numpyとの比較

Juliaを利用する一つの同期として、Pythonでfor文書くと遅い！というのがあると思います。Python/Numpyのnp.einsumを利用した丁寧な解説記事があったので、それをベースにPythonとJuliaのEinsum（と、その他のパフォーマンス）について比較します。

• 【einsum】アインシュタインの縮約記法のように使えるnumpyの関数。性能と使い方を解説。 - プロクラシスト
• 比較に関したレポジトリ（ファイル置き場とも）: https://github.com/cocomoff/SampleEinsum

行列積の比較

dim1を30に固定して、dim2を変化させた場合の計算時間を計測します。Python側はtimeitを、Julia側はtic()してtoq()した結果を計算時間と見なしています（これで良いのでしょうか？）

using Einsum

function benchmark(;dim1=30, dim2=20)
    A = randn(dim1, dim2)
    B = randn(dim1, dim2)
    C = randn(dim1, dim2)
    D = randn(dim1, dim2)
    E = randn(dim1, dim2)
    X = zeros(dim2, dim2, dim2, dim2, dim2)
    @einsum X[a,b,c,d,e] = A[r,a] * B[r,b] * C[r,c] * D[r,d] * E[r,d]
end

結果です。噂通り、高次元に対してJuliaの優位性が出てきます(もちろん、だからといってJuliaが素晴らしい！というつもりはないですが)。

実装上の話

numpyではOKだったけど、JuliaではOKじゃなかった例を示します（僕の実装のせいかもしれませんが）。まずはpythonではOKだったコードです。このようにPythonのnp.einsumでは、A[0]・A[1]・A[2]のような例を暗黙的に(A[0, i]・A[1, i]・A[2, i]のこと)入力することが出来ます。

A = np.arange(9).reshape(3, 3)
det_numpy = np.linalg.det(A)
det_einsum = np.einsum("ijk,i,j,k", eijk, A[0], A[1], A[2])
print(det_numpy, det_einsum)

一方で、これに類することを実行するに、PythonでA[0]と実装するところをA1 = A[:, 1]などと切り出してから、それを指定する必要があります。これは実装上、numpyではeinsumにおいて"ij->k"など、計算を文字列で指定しているのに対して、Juliaでは与えられた式をExprで格納して、それをパースしているため、インデクス(Pythonだと0、Juliaだと1)の型がSymbol型じゃないから起きます。先に述べた通り、ベクトルで切り出すと一応回避できます(他にも方法がありそうです)。

A = randn(3, 3)
det_julia = det(A)
A1 = A[:, 1]
A2 = A[:, 2]
A3 = A[:, 3]
@einsum det_einsum := eijk[i,j,k] * A1[i] * A2[j] * A3[k]