Test case:

data = RandomReal[1, {1000, 1000}];
pos = RandomInteger[{1, 1000}, {5*10^5, 2}];
r1 = ReplacePart[data, pos -> 0]; // AbsoluteTiming
(data[[##]] = 0) & @@@ pos; // AbsoluteTiming
data == r1

{3.433196, Null}
{1.892108, Null}

The similar matlab code only takes about 0.16 second:

pos= randperm(5e5);
(*Elapsed time is 0.166717 seconds.*)

I am looking for a more efficient way.

  • 6
    $\begingroup$ I doubt that you can make it any faster for as long as you insist on using ReplacePart (as you seem to be). It should be possible to speed it up with other methods. Here's one: SparseArray[pos -> 0., Dimensions[data], 1.] data; // AbsoluteTiming. It is a critical but not so obvious point that I used 0. and 1. (machine precision inexact numbers) instead of 0 and 1, to avoid unpacking data. $\endgroup$ – Szabolcs Sep 4 '13 at 16:18
  • 2
    $\begingroup$ I second @Szabolcs's comment. The big problem with ReplacePart is that it necessarily copies the entire expression. So, for large expressions the large overhead is pretty much inevitable with this method, particularly for comparatively small number of elements to be modified. $\endgroup$ – Leonid Shifrin Sep 4 '13 at 17:19

I think it's worth noting that the two random position codes you list aren't the same. If you want to do something like randperm in Mathematica you should use RandomSample. In addition to this, what you are trying to do is effectively to simultaneously address a large number of positions in the matrix. In Mathematica I do not believe you can do this for matrixes and other things with larger dimensions, however for a single list you can do this with speeds similar to MATLAB:

data = Flatten[RandomReal[1., {1000, 1000}], 1];
pos = RandomSample[Range[1, 5*10^5]];

Module[{temp = Flatten[data, 1]},
temp[[pos]] = 0.;
data = Partition[temp, 1000]
]; // AbsoluteTiming // First
(* 0.014001 *)

Where MATLAB returns Elapsed time is 0.058448 seconds. on my system. I believe the root of the difference here is that MATLAB treats your matrix like it was just a long list, and only uses the dimensions to translate between indexing using two coordinates and the singular index, while Mathematica has a more general structure, where you can't just assume that each row is the same length for the purpose of simultaneous indexing, which might be why you can't just do something like for instance data[[index[1,2],index[2,3]]]=0 for a matrix, even though you can do data[[{1,2}]]=0 for a list.

data = RandomReal[1, {1000, 1000}];
pos = RandomInteger[{1, 1000}, {5*10^5, 2}];

r1 = ReplacePart[data, pos -> 0]; // AbsoluteTiming

{3.145831, Null}

In your particular example you have many duplicate positions. If this can actually occur, then eliminate the duplicates before doing the replacements:

r2 = ReplacePart[data, Union[pos] -> 0]; // AbsoluteTiming

{1.447928, Null}

(data[[##]] = 0) & @@@ pos; // AbsoluteTiming

{1.670169, Null}

data == r1 == r2

  • 4
    $\begingroup$ Hi Bob, its great for you to join! This community will certainly benefit a lot from your presence here. $\endgroup$ – Leonid Shifrin Sep 4 '13 at 17:17
  • $\begingroup$ Wow.. Bob Hanlon here! Welcome, I'm big fan of your Mathematica way of think. $\endgroup$ – Murta Sep 4 '13 at 22:51

Your code using Part can be made more efficient by changing data one whole row at a time, like this:

With[{parts = {#[[1, 1]], #[[All, 2]]} & /@ GatherBy[pos, First]},
   (data[[##]] = 0) & @@@ parts];
  • $\begingroup$ Neat idea +1... $\endgroup$ – Rojo Sep 4 '13 at 17:57

...and you can always compile. Given that the MATLAB code you compare it to actually works with machine precision I assume that you don't mind the 0 you assign to be 0.

newData = Compile[{{data, _Real, 2}, {pos, _Integer, 2}},
    Module[{d = data}, Do[d[[i[[1]], i[[2]]]] = 0., {i, pos}];
     d]][data, pos];
  • 1
    $\begingroup$ I find the word "sped" weird. $\endgroup$ – Rojo Sep 4 '13 at 23:40

You can also flatten, change all at once, and then reshape

newData = 
  With[{flatPos = 
     Dimensions@data /. {width_, _} :> pos.{width, 1} - width},
   Module[{flatData = Flatten@data},
    flatData[[flatPos]] = 0.;

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.