ReplacePartの中身, Mathematica

 Mathematicaで以前からこれが遅くなる原因ではないか…と思っていた部分を検証したのでメモ．

概要

Mathematicaの代入

expr /. {replace1 -> value1, replace2 -> value2, ...}

は，代入されるもの・代入するものが多いと遅くなる (vectorizationできるならすべき)．

これの類推で，expr (ここでは行列) の一部を書き換える関数ReplacePart

    ReplacePart[expr, {{i1,j1} -> value1, {i2,j2} -> value2, ...}]

をTableを用いて

    ReplacePart[expr, Table[position[[i]] -> value[[i]], {i, Length@position}]]

のように書くのは遅くなる原因なのだろうか (7年ぐらい前に書いたコードはこれ)．

MapThreadを使うとベクトルをベクトル表記のまま同じことができる:

    ReplacePart[expr, MapThread[Rule, {position, value}]]

ので，Tableを使う場合とMapThreadを使う場合で差が出るか見てみた．

結果

MapThreadのほうが1割ぐらい速かった．実行環境はやや古いiMacのMathematica14.0.0.

コード

以下コード．

pattern01[nrow_, ncol_, nreplace_, nDo_] := 

 Module[{timebegin, mat, replacePosition, newValue},

  (* record computational time *)

  timebegin = SessionTime[] ;

  (*"set seed"*)

  ParallelEvaluate[SeedRandom[1]] (* enable to be replicated *);

  Do[

   (* a matrix that I want replace with newValue *)

   mat = RandomReal[1, {nrow, ncol}];

   (* positions of mat to be replaced. 

   the number of replace is at most nreplace. *)

   replacePosition = 

    DeleteDuplicates@

     Transpose@{RandomInteger[{1, nrow}, nreplace], 

       RandomInteger[{1, ncol}, nreplace]};

   (* new values (taking Exp is optional) *)

   newValue = Exp@RandomReal[1, Length@replacePosition];(* 

   replacepart with Table *)

   mat = 

    ReplacePart[mat, 

     Table[replacePosition[[i]] -> newValue[[i]], {i, 

       Length@replacePosition}] ];,

   nDo] (* end of Do *);

   Print[ToString@(SessionTime[] - timebegin) <> " sec"];

  Return@DateString[]

  ]

pattern02[nrow_, ncol_, nreplace_, nDo_] := 

 Module[{timebegin, mat, replacePosition, newValue},

  (* record computational time *)

  timebegin = SessionTime[] ;

  (*"set seed"*)

  ParallelEvaluate[SeedRandom[1]] (* enable to be replicated *);

  Do[

   (* a matrix that I want replace with newValue *)

   mat = RandomReal[1, {nrow, ncol}];

   (* positions of mat to be replaced. 

   the number of replace is at most nreplace. *)

   replacePosition = 

    DeleteDuplicates@

     Transpose@{RandomInteger[{1, nrow}, nreplace], 

       RandomInteger[{1, ncol}, nreplace]};

   (* new values (taking Exp is optional) *)

   newValue = Exp@RandomReal[1, Length@replacePosition];(* 

   replacepart with MapThread *)

   mat = 

    ReplacePart[mat, MapThread[Rule, {replacePosition, newValue}] ];,

   nDo] (* end of Do *);

   Print[ToString@(SessionTime[] - timebegin) <> " sec"];

  Return@DateString[]

  ]

 

"Parallel Karnels" (* 使ってないけど *)

LaunchKernels[];

"run (execution requires about 20min)"

 (* 16384x128行列のうち最大1024箇所を書き換える，を65536回実行 *) 

pattern01[2^14, 2^7, 2^10, 2^16] 

(* 1420 sec *)

pattern02[2^14, 2^7, 2^10, 2^16] 

(* 1318 sec *)

その他

Parallelize@ReplacePartのようにして並列化はできなさそう．