Mathematicaにおけるプログラムの高速化手法
Mathematicaにおいてプログラムの実行速度を最適化する際の項目を思いつく限り挙げてみた。
関数型パラダイムで書く
必然的に組み込み関数を多く使い、リストをまとめて操作することになるので手続き型で書くより速くなることが多い。 コード量も少なくなって読みやすくなるので、よほどのことでない限りMathematicaでは関数型で書く。
具体的には、手続き型ループ構文(Do, For, Whileなど)をやめて、MapやThreadを使うようにする。
出来る限り組み込みの関数を使い、呼び出し回数を減らす
組み込み関数でできることは出来る限りやらせる。 ドキュメントを探すと、Mathematicaは意外と多くのものが組み込みでできるようになっている。
例えば、整数商を求める場合はFloor[a / b]とするより組み込み関数一発でQuotient[a, b]として関数呼び出し回数を減らすと高速化する。
v = RandomInteger[{1, 10}, {10000, 2}];
{Timing[Floor[#[[1]]/#[[2]]] & /@ v;],
Timing[Quotient[#[[1]], #[[2]]] & /@ v;]}
{{0.037664, Null}, {0.00116, Null}}
リストは全体をまとめて扱う
大概のループの中身はリストへの関数の適用や、テンソルの演算などに帰着する。
四則演算はリストにそのまま適用可能。
内積はDot、外積はCrossがあるし、それらを一般化したInner、Outerといった関数もある。
Mathematicaでは個々の要素に対する操作は記述しなくてもよい場合が多い。
myDifferences1[v_] := Table[v[[i + 1]] - v[[i]], {i, Length[v] - 1}]; (* 個々の要素を操作する *)
myDifferences2[v_] := Most[RotateLeft[v] - v]; (* リスト全体をまとめて操作する *)
{Timing[myDifferences1[#];],
Timing[myDifferences2[#];],
Timing[Differences[#];]} &[RandomReal[{0, 5}, 5000000]]
{{1.919, Null}, {0.296, Null}, {0.187, Null}}
リストの要素数変更を避ける
手続き型で要素数が分からないリストを生成するには、ReapとSowを用いる。AppendToやPrependToはかなり遅く、厳禁。
myPrimeChoice1[n_] := Module[{r = {}, i}, Do[If[PrimeQ[i], AppendTo[r, i]], {i, n}]; r];
myPrimeChoice2[n_] := Module[{i}, Reap[Do[If[PrimeQ[i], Sow[i]], {i, n}]][[2, 1]]];
{Timing[myPrimeChoice1[200000];], Timing[myPrimeChoice2[200000];]}
{{3.01, Null}, {0.297, Null}}
無駄な式の生成を避ける
返り値を使わないのであればMapではなくScanを使う(効果は小さい)。
無駄な型変換を避ける
最終結果が機械精度実数や機械精度複素数でほしいのならば、最初からその精度で計算する。桁落ちや情報落ち、丸め誤差の蓄積には注意。
{Timing[N[Total[Table[1/i^2, {i, 1, 100000}]]]],
Timing[Total[Table[1./i^2, {i, 1, 100000}]]],
Timing[Total[Table[1./i^2, {i, 1., 100000.}]]]}
{{2.418, 1.64492}, {0.344, 1.64492}, {0.046, 1.64492}}
無駄な遅延評価を避ける
SetDelayed(:=)やRuleDelayed(:>)ではなくSet(=)やRule(->)でも良い場合はかならずSetやRuleを使う。
SetDelayedの右辺やTableの中身などを予め記号的に評価できるときはEvaluateを使って先に評価してしまう。
func1[x_] := Integrate[Exp[-x0^2] + x0^2 + x0 + 1, {x0, 0, x}];
func2[x_] := Evaluate[Integrate[Exp[-x0^2] + x0^2 + x0 + 1, {x0, 0, x}]];
{Timing[Do[func1[i], {i, 0, 100}]], Timing[Do[func2[i], {i, 0, 100}]]}
{{7.161, Null}, {0., Null}}
無駄な記号的評価を避ける
PlotやPlot3Dなどグラフを生成する関数や、
FindRoot、NIntegrate、NDSolveなど数値解析を行う関数では引数を記号的に処理することによってより適切なアルゴリズムを選択しようとするが、
式が複雑で記号的評価に時間がかかってしまう場合や解析的でないことが分かっている場合は関数の定義を引数が数値の場合に限ることで記号的評価を抑制できる。
Takagi[x_] := Evaluate[Sum[Abs[x 2^n - Round[x 2^n]]/2^n, {n, 0, 16}]]; (* 高木関数 *)
Takagi2[x_?NumberQ] := NSum[Abs[x 2^n - Round[x 2^n]]/2^n, {n, 0, 16}];
{Timing[Plot[Takagi[x], {x, 0, 1}];], Timing[Plot[Takagi2[x], {x, 0, 1}];]}
{{1.256, Null}, {0.604, Null}}
Compileであらかじめコンパイルする
引数やモジュール変数、戻り値がすべて機械範囲整数、機械精度実数、機械精度複素数、真偽値およびこれらの完全配列で、単純な演算を大量に行う関数の場合は予めCompileしておくことで高速化できる。
sin1[x_] := Evaluate[N@Normal@Series[Sin[x], {x, 0, 50}]];
sin2 = Compile[{x}, Evaluate[N@Normal@Series[Sin[x], {x, 0, 50}]]];
{Timing[Do[sin1[i], {i, 0, 1, 0.00001}]],
Timing[Do[sin2[i], {i, 0, 1, 0.00001}]]}
{{3.697, Null}, {0.624, Null}}
CompileでCコードを生成してコンパイルする
Jan. 8 2015追記:
Mathematica 8からの新機能で、CompileでCompilationTarget -> "C"を指定することで、Cのコードを生成して機械語にコンパイルすることができる。
さらにParallelization -> Trueを指定すると複数のスレッドを使うようになるので、さらに高速化する場合がある。
使用するCコンパイラによっても性能が変わるので、CCompilerDriver`$CCompileを変更して使用するCコンパイラを変更できる。
exp1[x_] := Module[{c = 1., f = 1., x2 = 1.}, Do[x2 = x2*x; c += x2/f; f *= i + 1, {i, 1, 100}]; c];
exp2 = Compile[x, Module[{c = 1., f = 1., x2 = 1.}, Do[x2 = x2*x; c += x2/f; f *= i + 1, {i, 1, 100}]; c]];
exp3 = Compile[x, Module[{c = 1., f = 1., x2 = 1.}, Do[x2 = x2*x; c += x2/f; f *= i + 1, {i, 1, 100}]; c], CompilationTarget -> "C"];
exp4 = Compile[x, Module[{c = 1., f = 1., x2 = 1.}, Do[x2 = x2*x; c += x2/f; f *= i + 1, {i, 1, 100}]; c], CompilationTarget -> "C", Parallelization -> True];
{Timing[Do[exp1[i], {i, 0, 1, 0.00001}]],
Timing[Do[exp2[i], {i, 0, 1, 0.00001}]],
Timing[Do[exp3[i], {i, 0, 1, 0.00001}]],
Timing[Do[exp4[i], {i, 0, 1, 0.00001}]]}
{{23.6656, Null}, {0.757096, Null}, {0.145064, Null}, {0.096116, Null}}
Dispatchで呼び出しテーブルを作っておく
大量の変換規則をReplaceAll(/.)で適用するときには予めDispatchで呼び出し表を作成しておく。
rule1 = MapIndexed[x[#2[[1]]] -> #1 &, RandomInteger[10000, 10000]];
rule2 = Dispatch[rule1];
{Timing[Do[x[i] /. rule1, {i, 10000}]],
Timing[Do[x[i] /. rule2, {i, 10000}]]}
{{10.124, Null}, {0.016, Null}}
SparseArrayを使う
ほとんどの成分が0である疎行列であればSparseArrayを用いると、疎行列用のコードで計算されるので高速化する。
matrixa = SparseArray[{Band[{1, 1}] -> 2, Band[{2, 1}] -> 1, Band[{1, 2}] -> 1}, {500, 500}];
matrixb = Normal[matrixa];
{Timing[matrixa.matrixa;], Timing[matrixb.matrixb;]}
{{0., Null}, {2.901, Null}}
パックアレーにする
すべての要素が機械範囲整数か機械精度実数、機械精度複素数のいずれか1種類の完全配列であればMathematicaはパックアレーを使う(こともある)。 パックアレーとそうでない通常のリストの間の変換には時間がかかるので、できるかぎりパックアレーの状態を保つようにする。
unpacked = Developer`FromPackedArray[Table[i + 2 j, {i, 1000}, {j, 1000}]];
packed = Developer`ToPackedArray[unpacked];
{Timing[Developer`PackedArrayQ[Nest[Transpose, unpacked, 50]]],
Timing[Developer`PackedArrayQ[Nest[Transpose, packed, 50]]]}
{{2.558, False}, {0.375, True}}
パターンマッチではPatternTest(?)、Condition(/;)を避ける
パターンマッチではCondition(/;)よりPatternTest(?)のほうが若干高速であるが、これらを使わずにすむのであればそちらのほうが高速である。
{Timing[Count[RandomInteger[{1, 100}, 1000000]/RandomInteger[{1, 100}, 1000000], x_ /; IntegerQ[x]]],
Timing[Count[RandomInteger[{1, 100}, 1000000]/RandomInteger[{1, 100}, 1000000], _?IntegerQ]],
Timing[Count[RandomInteger[{1, 100}, 1000000]/RandomInteger[{1, 100}, 1000000], _Integer]]}
{{1.903, 47767}, {1.809, 48343}, {1.42, 48168}}
Listable属性を活用する
Listable属性をうまく使えば、MapやMapThreadを使わずにリストへ関数を適用でき、高速化できる。
h1[x_, y_] := Boole[x^2 + y^2 <= 1.];
h2[x_, y_] := Boole[myLessEqual[x^2 + y^2, 1]];
myLessEqual[a_, b_] := a <= b;
SetAttributes[myLessEqual, Listable];
{Timing[Total[MapThread[h1, RandomReal[{}, {2, 400000}]]]],
Timing[Total[h2 @@ RandomReal[{}, {2, 400000}]]]}
{{2.418, 314145}, {1.155, 314415}}
求める精度や確度を下げる
精度が必要ない場面ではNIntegrateやFindRootのPrecisionGoalやAccuracyGoalを小さく設定することで関数の評価回数や再帰の深さを小さくできる。
Quiet[{Timing[Do[NIntegrate[Sin[1/x], {x, 10^-i, 1}], {i, 3, 200}]],
Timing[Do[NIntegrate[Sin[1/x], {x, 10^-i, 1}, AccuracyGoal -> 3, PrecisionGoal -> 3], {i, 3, 200}]]}]
{{2.059, Null}, {0.952, Null}}
あらかじめInterpolatingFunctionにしておく
- 複雑な計算をして数値を返す
- 何回も呼び出す
- 引数の範囲が限定されている
- 精度が必要ない
という関数はInterpolatingFunctionにしてしまう手もある。
ある関数の中でInterpolatingFunctionを呼び出して使うような場合は精度がそのInterpolatingFunctionによって決まってしまうので、
FunctionInterpolationをかけて単一のInterpolatingFunctionにしてしまってもよい。
目安としてInterpolatingFunctionのメッシュ点数よりも関数の呼び出し回数の方が多ければ全体としてみて高速化できるだろう。
並列計算させる
Mathematica 7から標準でついてくるようになった自動並列化を用いる。 個々の計算の粒度が粗いほうが並列化の効果があがる。
副作用のない関数型でプログラムを書いていればDistributeDefinitionsしたあと、
適当にParallelizeを入れたり、MapをParallelMapに置換するだけで並列化できる。
とにかく少しでも高速化するための小技
労力の割に報われずプログラムが読みにくくなる、やらないほうがいいことも多い小技。
[[1]]や[[-1]]ではなくFirstやLastを使う- 0との大小比較には
Positive、Negative、NonPositive、NonNegativeを使う - 割り算を掛け算にする
- 割り算(
i/j)をDivide[i, j]にする(遅延評価されないと無意味) - 整数倍を足し算にする(遅延評価されないと無意味)
- 整数乗を掛け算にする(遅延評価されないと無意味)
- 機械精度実数の分数は
1./6.ではなく6^^0.1(6進数の0.1)などと入力する(0.16666666666666666`よりは読みやすい) Andの中はFalseになる確率が高い順番に、OrではTrueになる確率が高い順番に条件式を記述する
番外編1:Cで書いたプログラムをMathLink経由で呼び出す
これが出来る人はこんなページ見なくても大丈夫だろうなぁ･･･。 ドキュメントセンターの該当ページを参照。
Jan. 8 2015 追記:
Mathematica 8からCompileでCコードを生成して自動でMathLink経由で呼び出せるようになった(上記CompileでCコードを生成してコンパイルするを参照)。
番外編2:MathCode C++かMathCode F90を買う
お金があるならどうぞ。