Mathematicaにおけるプログラムの高速化手法

Mathematicaにおいてプログラムの実行速度を最適化する際の項目を思いつく限り挙げてみた。 高速化に関する公式のドキュメントはこちら。

関数型パラダイムで書く

必然的に組み込み関数を多く使い、リストをまとめて操作することになるので手続き型で書くより速くなることが多い。 コード量も少なくなって読みやすくなるので、よほどのことでない限りMathematicaでは関数型で書く。

具体的には、手続き型ループ構文(Do, For, Whileなど)をやめて、MapやThreadを使うようにする。

出来る限り組み込みの関数を使う

組み込み関数でできることは出来る限りやらせる。 ドキュメントを探すと、Mathematicaは意外と多くのものが組み込みでできるようになっている。

リストは全体をまとめて扱う

大概のループの中身はリストへの関数の適用や、テンソルの演算などに帰着する。

四則演算はリストにそのまま適用可能。 内積は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[[email protected]@Series[Sin[x], {x, 0, 50}]];
sin2 = Compile[{x}, Evaluate[[email protected]@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}}

これは http://moriyan.jp/karakuri/nokisita/BBS_MSG_031221003623.html を参照すること。

求める精度や確度を下げる

精度が必要ない場面では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では逆の順番に条件式を記述する

番外編1:Cで書いたプログラムをMathLink経由で呼び出す

これが出来る人はこんなページ見なくても大丈夫だろうなぁ・・・。 ドキュメントセンターの該当ページを参照。

Jan. 8 2015 追記: Mathematica 8からCompileでCコードを生成して自動でMathLink経由で呼び出せるようになった(上記CompileでCコードを生成してコンパイルするを参照)。

番外編2:MathCode C++かMathCode F90を買う

お金があるならどうぞ。