数値的な例

Wolfram LibraryLink でアクセスできる関数を使うと,Mathematica の柔軟性と汎用性を維持しながら数値計算を最適化することが可能となる.既存の大規模な数値コードがある場合は,MathLinkLibraryLink の両方とも,Mathematica から駆動するコードを繋ぎ合わせるよい方法となる.一方数値計算を開発しているなら,Mathematica でプロトタイプを作成できるし,どこか一部が障害となっている場合はLibraryFunctionを使ってその部分の高速化を図ることができる.LibraryFunctionMathematica の数値計算関数とも直接繋がっているため,別のサンプル点で評価されるコードを書くことで,計算が高速化できることもある.このチュートリアルでは,LibraryLink の後者2つの使用方法に焦点を当てる.

このチュートリアルで示される例題のソースはドキュメントパクレットにある.

demo_numerical.cコンパイルすると共有ライブラリlibrary demo_numericalになる数値例題のソースコード

これは次の入力を評価すると見付けられる.

In[9]:=
Click for copyable input
Out[9]=
ファイルのソースを示す.

基本的な二次関数の例

以下は を計算するだけの非常に簡単な関数例である.

DLLEXPORT int parabola(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res)
{
    mreal x, a, f;
    if (Argc != 2) return LIBRARY_FUNCTION_ERROR;
    x = MArgument_getReal(Args[0]);
    a = MArgument_getReal(Args[1]);
    f = x*x - a;
    MArgument_setReal(Res, f);
    return LIBRARY_NO_ERROR;
}
関数をロードする.
In[15]:=
Click for copyable input
Out[15]=

関数がロードされたら,プロットしたり,サンプルしたり等,通常の関数と同じように使える.

関数の面プロットを作成する.
In[24]:=
Click for copyable input
Out[24]=
の関数として根をプロットする.
In[27]:=
Click for copyable input
Out[27]=

このような簡単な関数については,関数の評価のコストと比較して,Mathematica から結果を得たり Mathematica に結果を得たりするオーバーヘッドが大きいため,LibraryFunctionを使うことについて実際の利点はない.

マンデルブロ(Mandelbrot)集合

一方,関数が評価するのにより高価な場合は,大きな利点があることがある.その一例がマンデルブロ集合の画像作成に使える関数である.

以下に示す関数mandelbrotから始めて,列が有界でないTemplateBox[{{z, _, n}}, Abs]=4の外に至るまで, の繰返し数を計算するものである.mandelbrot_max_iterationsに等しくなるまで繰返しが続くと,点 は集合内であると想定され,0が返される

static mint mandelbrot_max_iterations = 1000;

DLLEXPORT int mandelbrot(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res)
{
    mint n = 0;
    mcomplex z = {0.,0.};
    mcomplex c = MArgument_getComplex(Args[0]);
    mreal rr = mcreal(z)*mcreal(z);
    mreal ii = mcimag(z)*mcimag(z);
    while ((n < mandelbrot_max_iterations) && (rr + ii < 4)) {
        mcimag(z) = 2.*mcreal(z)*mcimag(z) + mcimag(c);
        mcreal(z) = rr - ii + mcreal(c);
        rr = mcreal(z)*mcreal(z);
        ii = mcimag(z)*mcimag(z);
        iteration++;
    }

    if (n == mandelbrot_max_iterations)
        n = 0;

    MArgument_setInteger(Res, n);
    return LIBRARY_NO_ERROR;
}

Cコードには複素演算が組み込まれていないため,数学演算は実数演算を使って実装されている.

関数をロードする.
In[1]:=
Click for copyable input
Out[1]=

これで関数が複素平面上の点群を評価する Mathematica の通常の関数のように使えるようになった.

関数の面プロットを作成する.
In[2]:=
Click for copyable input
Out[2]=

関数は非常に不連続であるので,これは集合を可視化する最良の方法ではないだろう.より一般的で格段に速い方法は,通常の点の格子上で関数をサンプリングして画像を作成することである.

関数を通常の点の格子上でサンプリングする.
In[3]:=
Click for copyable input
Out[3]=
二値画像を作成する.
In[4]:=
Click for copyable input
Out[4]=

関数値から色チャンネルへのマップを作成しRasterを使うことでカラー画像が作れる.Imageを使うこともできるが,下に示すようにインタラクティブに使う場合はRasterの方が格段に速く表示できる.

カラー画像を作成する.
In[5]:=
Click for copyable input
Out[6]=

ランダムなカラーマップで,逃避時間の差をハイライトする.

In[7]:=
Click for copyable input
Out[8]=

このようなもののゴールは,画像からできるだけ素早くデータを取り出すことであることがよくある.上の例では,C関数が行う計算の主要部分は非常に素早く行われるが,残りの部分は Mathematica コンパイラを使うと高速化できる.Mathematica コンパイラとLibraryFunctionは緊密に統合されており,これにより非常に効率的になる.

CompiledFunctionがカラー画像に色を与えるようにする.
In[9]:=
Click for copyable input

ここでWithを使うのは,Compileは自身の引数を評価せず,は使われるとすぐに外部評価によってランタイムで解決する必要のあるシンボルであると想定されるからである.それらを明示することで,コンパイラは可能な限り最も効率的な評価を行う.

Listable属性とはリストで与えられたいくつかの複素点について評価できるということである.

虚軸上のいくつかの点についてCompiledFunctionを評価する.
In[10]:=
Click for copyable input
Out[10]=

これはつまり,色チャンネルデータを取得するのに必要なのは,この関数に必要な全サンプル点を与えるこどだけであるということである.これを素早く行う方法の一つは,コンパイル済み関数を書くということである.

複素平面の通常の格子状のサンプル点を与えるCompiledFunctionを定義する.画像の中で上昇することがの増加となるよう,の方向を逆にする.
In[11]:=
Click for copyable input
からまでの平面の長方形上にあるマンデルブロ集合の×の画像を与える関数を定義する.
In[12]:=
Click for copyable input
マンデルブロ集合の12×512の画像を作成するのに必要な全時間を与える.
In[13]:=
Click for copyable input
Out[13]=

これらのコンパイル済み関数を使って画像全体を取得するのに必要な時間は,単にサンプル点で関数を評価するのより速くなる.これは上記サンプル点の評価に示される計算時間の多くはの代りにListとしてTableを構築するの使われているからである.コンパイル済み関数では全体を通してパックアレーが使われる.

複数のコアを持つマシンを使用しているなら,並列化によってさらに計算速度を増加させることができる.

並列化に使うために,色チャンネルの計算のためのCompiledFunctionを再定義する.
In[14]:=
Click for copyable input
マンデルブロ集合の512×512の画像を並列評価で行ったときの所要時間を示す(ここで示されている時間は8コアマシンのものである).
In[15]:=
Click for copyable input
Out[15]=

高速化は顕著ではあるが,プロセッサ数に比例はしない.これはサンプル点の構築や画像の生成等,1つのプロセッサ上で行われる部分もあるからである.

このように速い更新では,集合のいくつかの側面が調べられるインタラクティブなインターフェースを設定しておくと便利である.EventHandlerおよびDynamicを使った設定はやや複雑であるが,効率的に動作させるようにする決め手は,上で設定したmimage関数のスピードである.

In[21]:=
Click for copyable input

マシンがかなり速ければ,各方向512ポイントでこれを実行することができる.そうでなくても200ポイントで適切な画像が得られる.子rを実行するためには

を使うとよい.ここで n は各方向の点の数を示す.ズームの倍率をかなり大きくすると,2つの理由で鮮明さが失われる.その一つは,境界を鮮明にするにはCコードでエンコードされた1000回の繰返し上限以上が必要になるという点である.2つ目は場合によっては近隣の点の軌道を実際に解像するのに,より高度な精度が要求されることがあるという点である.

Mathematica コンパイラは,コンパイルするとLibraryFunctionとしてロードできる共有ライブラリになるCコードを自動的に生成する.例えば,生成されたCコードはより遅いが,手で書いた関数と比べるとそれほど遅くない.

マンデルブロ集合の繰返し数を評価するCompiledFunctionを定義する.
In[17]:=
Click for copyable input
関数をマンデルブロ集合の画像を生成する関数に結合する.
In[18]:=
Click for copyable input
生成されたコードを使って並列評価でマンデルブロ集合の512×512の画像を作成するのに必要な合計時間を求める(ここで示されている時間は8コアマシンのものである).
In[20]:=
Click for copyable input
Out[20]=

手書きのCコードの約2倍の時間がかかっている.

ブリュセレータ(Brusselator)偏微分方程式

化学反応のモデル化のためのブリュセレータ偏微分方程式問題([HNW93]と[HW96]に記載されている)は,記号前処理のコンポーネントをいくつか含むNDSolve等の関数で使われるLibraryFunctionを定義するよい例である.

これから解く偏微分方程式問題のパラメータと方程式を設定する.
In[25]:=
Click for copyable input

問題はNDSolveで直接解くことができる.この例の論点の一部は,NDSolveと部分的に手作業でコードを作成した離散化との解像の速度を比較するということである.

NDSolveで問題を解く.
In[28]:=
Click for copyable input
Out[28]=
解の面プロットを示す.
In[29]:=
Click for copyable input
Out[29]=

参考文献[HW96]では二次の差分と 点の離散化を使って問題が設定されており,大きい については「硬い問題の全特徴を有している」と言っている.これと同じ離散化は,NDSolveを使って二次差分を強制し,線の方法のメソッドオプションを使って空間格子点の数を設定することで行える.を使い,結果のオブジェクトから関数を抽出してNDSolveで効率的に解き,常微分方程式の右辺を得ることもできる.

離散化方程式のオブジェクトを得る..
In[30]:=
Click for copyable input
Out[30]=
関数と初期条件をオブジェクトから得る.
In[31]:=
Click for copyable input

demo_numerical.cの関数は同じ離散化が返されるように定義されている.

をライブラリ関数としてロードする.
In[33]:=
Click for copyable input
Out[33]=

は離散化された右辺に明示的には現れないが,関数はそれを第1引数に取るように定義されている.これによりNDSolveは非常に非効率的なインターフェースを作ることができる.この関数が与える値はNDSolveが生成した関数が与える値とほぼ同じである.

関数の値を比較する.
In[34]:=
Click for copyable input
Out[34]=

わずかな差分は単に計算操作の順における小さな差のためである.

評価速度を比較する.
In[35]:=
Click for copyable input
Out[36]=

LibraryFunctionの方が約20倍速い.

NDSolveは偏微分方程式を解くのに,常微分方程式の一次系の右辺に関数を使用した.ベクトルの初期条件を使うと,これはNDSolveに直接与えることができる.

NDSolveの右辺関数を使って常微分方程式系を得.AccuracyGoalPrecisionGoalオプションは,偏微分方程式に対するNDSolveのものと矛盾しないように使われる.
In[37]:=
Click for copyable input
Out[37]=
LibraryFunction右辺関数を使って常微分方程式系を解く.
In[38]:=
Click for copyable input
Out[38]=

LibraryFunctionは評価するのに20倍近く速くても,解には3倍もかかる.これはヤコビアン(Jacobian)が必要とする硬さと,ヤコビアンが疎であるということに起因する.この問題を回避する方法は2つある.最も簡単なのは,ヤコビアンの疎のパターンを与え,比較的少ない評価で有限の差分が使えるようにし,疎な線形代数が使えるようにすることである.

疎のパターンを得るのが難しいことはよくあり,また非零であろう要素をすべて得ることが重要である.以下では,右辺が に依存していないという事実を利用し,有限差分のみ使用する.

ランダムなテストベクトルの有限差分を使ってヤコビアンの疎のパターンを取得する.
In[39]:=
Click for copyable input
Out[42]=
有限差分を使って疎になるように計算されたヤコビヤンにLibraryFunction右辺関数を使って常微分方程式系を解く.
In[43]:=
Click for copyable input
Out[43]=

ヤコビアンの計算を制御すると,格段に速くなる.

解析的ヤコビアンを定義するのは,パターンを正しく得なくてはならない上に,そのパターンを使って値のベクトルを与える関数から疎配列を構築しなければならないので,1ステップ複雑になる.

LibraryFunctionをロードして非零の位置と対応するヤコビアンの値を取得する.
In[44]:=
Click for copyable input
ヤコビアンをSparseArrayとして返す純関数を作る.
In[49]:=
Click for copyable input

ここではが一度だけ評価され,定義に位置が含まれるようにするためにWithが使われている.位置が時間やベクトルに依存する場合は,定義の中で明示的にを使えばよい.記号的評価を避けるためにPatternTest が使われているが,これにより構成が正しくなくなる.

解析的ヤコビアンの値と有限差分から計算された値を比較する.
In[50]:=
Click for copyable input
Out[50]=

有限差分は機械精度の半分までしか正確でないので,ここは差の大きさが想定される.

解析的ヤコビアンのLibraryFunction右辺関数を使って常微分方程式系を解く.
In[51]:=
Click for copyable input
Out[51]=

速度は有限差分を使った場合と同程度である.計算によってはヤコビアンの方が他のに比べてエラーになりやすいことがある.NDSolveについては,多くの場合ある程度近いヤコビアンで十分であることが多い.一方でFindRootFindMinimumのような関数については,導関数の確度は解の確度に大きく関係する.

NDSolveを使って特定の値の に対する右辺関数のためのCコードを生成することもできる.メソッドオプションを"ExpandEquationsSymbolically" -> Trueとすると,ベクトル関数を使う代りに線の方法が系を記号的に拡張する.Compiledオプションでは,その中で与えられたオプションをコンパイラに渡す.

解析的ヤコビアンのLibraryFunction右辺関数を使って常微分方程式の系を解く.
In[59]:=
Click for copyable input
Out[59]=

この処理を行うには長い時間がかかる.これは生成されたコードが非常に大きく,Cコンパイラがコードをコンパイルして最適化するのに時間がかかるからである.ヤコビアンが最初に評価されたときに,Cコードも生成されてコンパイルされるので,最初の評価には時間がかかる.

解析的ヤコビアンの値と生成されたコードから得られた値を比較する.
In[60]:=
Click for copyable input
Out[60]=
生成されたコードを使って常微分方程式の系を解く.
In[61]:=
Click for copyable input
Out[61]=

そのためコンパイルのための顕著な時間コストにおいて,生成されたコードの速度は手作業でコードした関数とほぼ同じである.生成されたコードはあまり柔軟ではないことに注意されたい. の別の値についてすべて再び生成しなければならない.

ダフィング(Duffing)方程式の折りたたみ

ダフィング方程式は定期的に駆動する固定梁であると考えることができる.これは非線形の復元力を持つ振動子を表す.

パラメータを持つダフィング方程式を定義する.
In[1]:=
Click for copyable input

これらのパラメータの値について,軌道は確定的にカオス的である.NDSolveが計算した解からこれを表示することができる.

から始まる近似解の位相面プロットを表示する.
In[3]:=
Click for copyable input
Out[3]=

カオス軌道は一般に初期条件に非常に敏感である.この例では,初期条件の線について軌跡の変化をどのように追跡できるかが示されている.

これを行う一つの方法に,曲線を特定の時間までに解決できるように線に十分近い空間から始めるということである.この方法は比較的簡単だが計算的に集中し,さらにいくつの点で開始すべきかを予測するのが難しい.

別の方法は,軌跡が進むにつれて特徴が自動的に解決されるように適応的に点を追加するというものである.これはすべての軌跡を時間を増やしながら進め,ある基準に応じて点を追加して行う.この例で使われた基準は単に点と点の間の空間である.曲率等,より高度なテストも使えるが,単純に距離を使うと物事がより簡単にすむ.

解を前進させる最も効率のよい方法は,ベクトルの初期条件で記述された一次系でNDSolveを使うものである.

解を一次系として時間 から まで前進させる関数を定義する.
In[4]:=
Click for copyable input
で始まる軌跡を,時間を数回だけ 増やして前進させ,軌跡の連続プロット上に点を示す.
In[25]:=
Click for copyable input
Out[8]=

NDSolveでは初期条件をベクトルにすることができるため,関数は多くの軌跡を一度に扱うことができる.初期値 xy は別々のベクトルなので,を使う代りにを使って xy の座標を分けたことに注意されたい.

時間を ずつ数回増やして異なる時間から始まる複数の軌跡を進め,それを示す.
In[9]:=
Click for copyable input
Out[12]=

点と点の間の距離が において大きく増加したことが分かる.この微調整方法は距離が大きくなりすぎる前に新しい点を加えるというものである.

指定の点集合について,点と点の間の最大距離が tol である,...を含む新しい点集合を与える関数を定義する.
In[13]:=
Click for copyable input
In[14]:=
Click for copyable input

関数にはいくつでも点が使えるため,指定の初期条件集合を前進させるためにしなければいけないことは,を使って点が十分に近いようにし,これらの点をすべて進め,それを繰り返すことである.

次の例は初期条件の線について,点と点の間の距離を0.025に保ちながら解を微調整し,前進させる.
In[15]:=
Click for copyable input
Out[19]=

の値をより大きくすると,計算に長い時間がかかることがある.これは隣接する2点は指数関数的に発散することが多く,加えられる点の数も時間に対して指数関数的に増えるからである.

計算時間を減らすには,まずどこに多くの時間がかけられているかを見極める必要がある.この計算では,考えなければいけない要素は,微調整と解の前進の2つだけである.時間の多くがどこに使われているかを見るには,軌跡が延びるにしたがって各ステップが必要とする時間をプロットすればよい.

初期条件の線を までの刻み幅で進めるときに,が各ステップで必要とする時間のプロットを作成する.
In[20]:=
Click for copyable input
Out[20]=

プロットを見ると,各ステップの微調整部分で半分以上の時間が使われているのが分かる.したがって,その部分を早くすると,計算全体の速度に最も大きく影響があるであろうと考えられる.

微調整関数はパックアレーとしてコンパイルや表現ができない異なる長さの配列を扱うため,Mathematica ではスピードアップが難しい.したがって,この部分を行うためのLibraryFunctionを設定するのが理にかなっている.

上で定義した関数と同じことを行うLibraryFunctionをロードする.
In[21]:=
Click for copyable input
初期条件の線を までの刻み幅で進めるときに,LibraryFunctionが各ステップで必要とする時間のプロットを作成する.
In[22]:=
Click for copyable input
Out[22]=

微調整にLibraryFunctionを使うと,計算が3倍早くなり,微分方程式の解の計算においては,計算時間が大半を占めることになる.

多くの場合,NDSolveに組み込まれたより高度な適応時間ステップ法を使って微分方程式を数値的に解くと,より効率的である.しかしこの例では,次の2つの事柄が組み合さって,簡単な固定ステップサイズのソルバを書くことが望ましい状態になっている.1つ目はメソッドによらず,局所誤差は最終的に指数関数的に大きくなるため,計算された解は実際の解を暗示しているにすぎず,長時間実行する場合は局所誤差を微調整しても,その誤差が大きすぎない限りは大きな違いは生じない.2つ目は,比較的高次の方法において,微調整間のステップは常微分方程式ソルバに対して適切であるということである.したがって,微調整はオーバーヘッドなしで非常に素早く結果を得るためのソルバの一ステップのみで置き換えることが可能である.ここに適したソルバは,古典的な4次のルンゲクッタ(Runge-Kutta)法であるが,より低い次数の方法で時間をさらに削減することも可能であるかもしれない.

軌跡を古典的な四次のルンゲクッタ法の1ステップで軌跡を前進させる関数をロードする.
In[23]:=
Click for copyable input
Out[23]=
LibraryFunctionを使って前進させ,の計算を行う.
In[24]:=
Click for copyable input
Out[24]=

単一ステップの前進LibraryFunctionは格段に速いことが分かる.時間の差は小さすぎてあまり意味を成さないので,プロットしなかった.さらに先の時間まで行き,時間の差をプロットするなら,時間の大半が微分方程式の解法に使われていることが分かるだろう.

繰返しが十分に最適化されたので,さらに大きな計算ができるようになった.比較的小さい列から始めて展開していくことで得られるプロットの列は,どのように近くの初期条件が伸ばされ折りたたまれていくかを示す.

初期条件の線を まで進め,時間区間ごとにプロットを作成する.
In[25]:=
Click for copyable input
Out[25]=
10番目ごとのプロットを示す.
In[26]:=
Click for copyable input
Out[26]=

非常に小さな線分がやがて引き伸ばされ,折れ曲がって元に戻ってくることが分かる.これを見るのに適した方法は,ListAnimate[plots]を入力することである.サイズの都合上,出力はこのノートブックには含めない.

New to Mathematica? Find your learning path »
Have a question? Ask support »