Can ReplacePart be sped up?

Question

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}
True
*)

The similar matlab code only takes about 0.16 second:

tic;
A=rand(1000,1000);
pos= randperm(5e5);
A(pos)=0;
toc;
(*Elapsed time is 0.166717 seconds.*)

I am looking for a more efficient way.

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. — Szabolcs, Sep 04 '13 at 16:18
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. — Leonid Shifrin, Sep 04 '13 at 17:19

score 13 · Answer 1 · edited Sep 04 '13 at 17:17

13

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

True

edited Sep 04 '13 at 17:17

Leonid Shifrin

114,335
15
329
420

answered Sep 04 '13 at 17:11

Bob Hanlon

157,611
7
77
198

4

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

score 7 · Answer 2 · answered Sep 04 '13 at 17:52

7

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];

answered Sep 04 '13 at 17:52

Simon Woods

84,945
8
175
324

Neat idea +1... – Rojo Sep 04 '13 at 17:57

score 7 · Answer 3 · answered Sep 04 '13 at 17:58

7

...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];

answered Sep 04 '13 at 17:58

Rojo

42,601
7
96
188

1

I find the word "sped" weird. – Rojo Sep 04 '13 at 23:40

score 4 · Accepted Answer · answered Sep 04 '13 at 21:18

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:

SeedRandom[322112432]
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.

score 2 · Answer 5 · answered Sep 05 '13 at 02:06

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.;
    flatData~ArrayReshape~Dimensions@data
    ]
   ];

Can ReplacePart be sped up?

5 Answers5

Linked