Mathematicaにおけるコンパイラの効用

松田裕幸
計算環境 Mathematica 14.0、 MacOS Sonoma 14.3.1, CPU Apple M1, メモリ 16GB
Mathematica関数 Compile, AbsoluteTiming, Clear

はじめに

Mathematica にコンパイル機能があることは昔からよく知られていますが、いま一つ使い方がわからず、あまり利用されていません。本稿では具体的な例 https://reference.wolfram.com/language/ref/Compile.html.ja 内[おもしろい例題]を参考に、コンパイルすると高速化(この例の場合は約30倍に高速化)されることを示します。

基本事項の確認

コンパイラを知らない人はいないと思いますが、コンパイルの基本がメモリ参照をCPU内メモリ(普通はレジスタ)に切り替えると理解している方は少ないかもしれません。メモリアクセス速度を低い順に並べると ファイル << メインメモリ < CPUレジスタとなります。Mathematica は通常、インタプリタとして機能し、変数が登場するたびにメインメモリへのアクセスが発生します。これに対し、コンパイラは変数をCPUレジスタに割り付けるために、参照速度が上がります。​なおコンパイルにおいて当該データが何バイト使用するか事前に知る必要があり、Mathematica のコンパイル関数 Compile では引数に必ず型を加えます。たとえば、
{x,_Real
} 「引数
x
の型は実数」はその例です。​したがってコンパイルの対象変数は実行時、整数や実数等の数値でなければなりません。ただし、それらの変数が Sin 等の組み込み関数の引数あるいはリスト構造の一部であっても、コンパイルの効果はあります。

コンパイル関数 Compile の使い方

式
Sin[x]+x^2-1/(1+x)
を本体に持つ関数 cf を定義します。
In[]:=
Clear[cf]
In[]:=
cf[x_]:=Sin[x]+x^2-1/(1+x);​​cf[N[Pi]]
Out[]=
9.62815
この関数の本体部分をコンパイルするには引数部分
{{x,_Real}}
(二重括弧は複数の引数を想定)を含め、本体を Compile で囲みます。結果の CompileFunciton が代入された変数 cf は通常の関数名として使用できます:
\bcf[Pi]
。なお、引数が数値を仮定しているため、シンボル
Pi
を
N
によって強制的に数値化する必要はありません。
In[]:=
cf=Compile[{{x,_Real}},Sin[x]+x^2-1/(1+x)]​​cf[Pi]
Out[]=
CompiledFunction
Argument count: 1
Argument types:
{_Real}

Out[]=
9.62815
ここでコンパイルの効果をみるために、かなり不自然な関数(同じ計算を100万回
繰り返す)を定義します。現在、Mathematica の Do 関数はかなり最適化されていますが、それでもコンパイルするのとしないとでは約10倍の速度差が出ることがわかります。
In[]:=
Clear[cf]
In[]:=
cf[x_Real]:=Do[Sin[x]+x^2-1/(1+x),{1000000}]​​AbsoluteTiming[cf[N[Pi]];]
Out[]=
{1.1818,Null}
In[]:=
cf=Compile[{{x,_Real}},Do[Sin[x]+x^2-1/(1+x),{1000000}]];​​AbsoluteTiming[cf[N[Pi]];]
Out[]=
{0.10731,Null}

実例

以下の各小節に「原コード」と「コンパイル化コード」を挙げますが、その前に関数引数部分に関し、コンパイル化前と後の比較を行います。なお With は一時的な値を設定するのに使用し、Block のように内部の変数を局所化しません。つまり With 内で関数は大域的に定義されます。
比較1.
fade[t_]:=t*t*t*(t*(t*6.0-15.0)+10.0);​​fade=Compile[{{t,_Real}},t*t*t*(t*(t*6.0-15.0)+10.0)];
比較2.
lerp[x_,y_,t_]:=(1.0-t)*x+t*y;​​lerp=Compile[{{x,_Real},{y,_Real},{t,_Real}},(1.0-t)*x+t*y];
比較3.
With[{grad={{1,1,0},{-1,1,0},{1,-1,0},{-1,-1,0},{1,0,1},{-1,0,1},{1,0,-1},{-1,0,-1},{0,1,1},{0,-1,1},{0,1,-1},{0,-1,-1}}},​​dot[gradIdx_,x_,y_,z_]:=​​grad[[gradIdx+1]][[1]]*x+grad[[gradIdx+1]][[2]]*y+grad[[gradIdx+1]][[3]]*z​​];​​dot=With[{grad={{1,1,0},{-1,1,0},{1,-1,0},{-1,-1,0},{1,0,1},{-1,0,1},{1,0,-1},{-1,0,-1},{0,1,1},{0,-1,1},{0,1,-1},{0,-1,-1}}},​​Compile[{{gradIdx,_Integer},{x,_Real},{y,_Real},{z,_Real}},​​grad[[gradIdx+1]][[1]]*x+grad[[gradIdx+1]][[2]]*y+grad[[gradIdx+1]][[3]]*z​​]​​];
比較4.
With[{permutations=Join[(permutations=RandomSample[Range[0,255]]),permutations]},​​signedNoise[x0_,y0_,z0_]:=​​(*以降省略*)]​​signedNoise=With[{permutations=Join[(permutations=RandomSample[Range[0,255]]),permutations]},​​Compile[{{x0,_Real},{y0,_Real},{z0,_Real}},(*以降省略*)]]
比較5.
With[{octaves=8},​​classicPerlin[xIndex_,yIndex_,amplitude_,frequency_,gain_,lacunarity_,scale_,increment_,width_,height_]:=(*以降省略*)]​​classicPerlin=​​With[{octaves=8},Compile[{{xIndex,_Integer},{yIndex,_Integer},{amplitude,_Real},{frequency,_Real},{gain,_Real},{lacunarity,_Real},{scale,_Real},{increment,_Real},{width,_Integer},{height,_Integer}},(*以降省略*)]]​​
比較4.比較5.ではコンパイラオプション
CompilationOptions{"InlineExternalDefinitions"True},"CompilationTarget""WVM"
が指定されています。
CompilationTarget はデフォルトでは Cコードを想定します。ただし signedNoise のように計算よりリスト処理が中心の場合、Mathematicaの仮想マシン用コード WVM を指定し、高速化します。
InlineExternalDefinitions を指定した場合、signedNoise 内で使用される外部関数(たとえば dot 等)がインライン展開され、コンパイルされます。これによりコンパイルコードは長くなりますが、関数呼び出しが減る分、処理時間は高速化されます。

原コード

In[]:=
Clear[dot,fade,lerp,signedNoise,classicPerlin]
In[]:=
​​With[{grad={{1,1,0},{-1,1,0},{1,-1,0},{-1,-1,0},{1,0,1},{-1,0,1},{1,0,-1},{-1,0,-1},{0,1,1},{0,-1,1},{0,1,-1},{0,-1,-1}}},​​dot[gradIdx_,x_,y_,z_]:=​​grad[[gradIdx+1]][[1]]*x+grad[[gradIdx+1]][[2]]*y+grad[[gradIdx+1]][[3]]*z​​];​​​​fade[t_]:=t*t*t*(t*(t*6.0-15.0)+10.0);​​lerp[x_,y_,t_]:=(1.0-t)*x+t*y;​​​​With[{permutations=Join[(permutations=RandomSample[Range[0,255]]),permutations]},​​signedNoise[x0_,y0_,z0_]:=​​Module[{x,y,z,ix,iy,iz,g000,g001,g010,g011,g100,g101,g110,g111,n000,n100,n010,n110,n001,n101,n011,n111,u,v,w,nx00,nx01,nx10,nx11,nxy0,nxy1,nxyz},​​ix=IntegerPart[x0];​​iy=IntegerPart[y0];​​iz=IntegerPart[z0];​​​​x=x0-ix;​​y=y0-iy;​​z=z0-iz;​​​​ix=Mod[ix,255]+1;​​iy=Mod[iy,255]+1;​​iz=Mod[iz,255]+1;​​​​g000=Mod[permutations[[ix+permutations[[iy+permutations[[iz]]]]]],12];​​g001=Mod[permutations[[ix+permutations[[iy+permutations[[iz+1]]]]]],12];​​g010=Mod[permutations[[ix+permutations[[iy+1+permutations[[iz]]]]]],12];​​g011=Mod[permutations[[ix+permutations[[iy+1+permutations[[iz+1]]]]]],12];​​g100=Mod[permutations[[ix+1+permutations[[iy+permutations[[iz]]]]]],12];​​g101=Mod[permutations[[ix+1+permutations[[iy+permutations[[iz+1]]]]]],12];​​g110=Mod[permutations[[ix+1+permutations[[iy+1+permutations[[iz]]]]]],12];​​g111=Mod[permutations[[ix+1+permutations[[iy+1+permutations[[iz+1]]]]]],12];​​​​n000=dot[g000,x,y,z];​​n100=dot[g100,x-1,y,z];​​n010=dot[g010,x,y-1,z];​​n110=dot[g110,x-1,y-1,z];​​n001=dot[g001,x,y,z-1];​​n101=dot[g101,x-1,y,z-1];​​n011=dot[g011,x,y-1,z-1];​​n111=dot[g111,x-1,y-1,z-1];​​​​u=fade[x];​​v=fade[y];​​w=fade[z];​​​​nx00=lerp[n000,n100,u];​​nx01=lerp[n001,n101,u];​​nx10=lerp[n010,n110,u];​​nx11=lerp[n011,n111,u];​​​​nxy0=lerp[nx00,nx10,v];​​nxy1=lerp[nx01,nx11,v];​​​​nxyz=lerp[nxy0,nxy1,w];​​​​nxyz]];​​With[{octaves=8},​​classicPerlin[xIndex_,yIndex_,amplitude_,frequency_,gain_,lacunarity_,scale_,increment_,width_,height_]:=​​Module[{noiseVal=0.0,x,y,z,freq=frequency,amp=amplitude},​​x=xIndex*frequency/scale;​​y=yIndex*frequency/scale;​​z=1.0*frequency/scale;​​Do[​​noiseVal+=signedNoise[x*freq,y*freq,z*freq]*amp;​​freq*=lacunarity;​​amp*=gain,​​{octaves}​​];​​Min[Max[noiseVal,0.0],1.0]​​]];​​createImage[img_,r_,g_,b_]:=Image[{r*img,g*img,b*img},InterleavingFalse]​​perlin[width_Integer,height_Integer]:=Table[classicPerlin[ii,jj,amplitude,frequency,gain,lacunarity,scale,increment,width,height],{ii,0,width},{jj,0,height}];

コンパイル化コード