こんにちは、ほけきよです。
大学で物理*1を嗜んだ方ならわかるであろうEinsteinの縮約記号。
計算の上で色々省略できるしとにかく慣れれば色々便利な記法です！
物理学者以外には馴染みがなく微妙かもしれませんが、「便利そうだな〜〜」と思って使ってみたり試してみたりしたので、メモとしてまとめておきます。誰かの参考になれば幸いです。
einsum
Einstein(アインシュタイン)の縮約記法風に足し合わせができるnumpyの関数
縮約記法とは
たとえば、行列 
と 
っていうのがあって、行列の演算をしたいときバカ正直に書くと
という風に計算される。
こういう演算で、
が大量に出てきすぎるので、「じゃあ略してしまえ」となった。
ルールは
これだけ。さっきの足(
の
の部分)が
便利なのは、足の数が3つ以上ついている*2テンソル計算も簡単にできるということ。
物理学の世界(相対論など)では、この高次テンソルでの計算がたくさん出てくるため
「そんなに書いてられるか！」
ということで、略される。
例
- 行列演算(行列
) -> 
- ベクトルと
の内積 -> 
- 行列
のトレース -> 
文法
np.einsum('今ある足'->'残したい足')
これだけじゃいまいちわからないと思うので、実際に使ってみる。
行列(単体)
1つの行列で色々計算してみる
import numpy as np
a = np.arange(25).reshape(5,5)
b = np.arange(25).reshape(5,5)*2
print("A = ",a)
print("B =", b)
>>> A = [[ 0 1 2 3 4]
>>> [ 5 6 7 8 9]
>>> [10 11 12 13 14]
>>> [15 16 17 18 19]
>>> [20 21 22 23 24]]
>>>B = [[ 0 2 4 6 8]
>>> [10 12 14 16 18]
>>> [20 22 24 26 28]
>>> [30 32 34 36 38]
>>> [40 42 44 46 48]]
np.einsum("ii",a)
np.einsum("ii->i",a)
np.einsum("ij->i",a)
np.einsum("ij->j",a)
np.einsum("ij->ji",a)
>>> array([[ 0, 5, 10, 15, 20],
>>> [ 1, 6, 11, 16, 21],
>>> [ 2, 7, 12, 17, 22],
>>> [ 3, 8, 13, 18, 23],
>>> [ 4, 9, 14, 19, 24]])
行列演算その①
2つ以上の行列の演算を行う(使うのはさっきのA,B)
np.einsum("ij,jk->ik",a,b)
np.einsum("ij,jk->ki",a,b)
np.einsum("ij,ij->ij",a,b)
ベクトル演算その②
ベクトルの演算も、行列とほぼ同様にできる
v1 = np.arange(3)
v2 = np.arange(3)+1
np.einsum("i,i",v1,v2) #内積
np.einsum("i,j->ij",v1,v2) #直積
テンソル表記の時に気になるのが、行列式や逆行列、外積の計算ってどうするのよ？？ということ。
こういう時に便利なのがレビチビタ(Levi Civita)記号。定義は以下の通り
hogehoge
3次元以上空間の時には、この記号が計算をスッキリまとめられて便利。
これを使うと、einsumの枠組みでより多様な計算が実現できる
ただし、レビチビタ記号は調べた範囲では実装されてなさそう？*3なので自分でつくる
eijk = np.zeros((3, 3, 3))
eijk[0, 1, 2] = eijk[1, 2, 0] = eijk[2, 0, 1] = 1
eijk[0, 2, 1] = eijk[2, 1, 0] = eijk[1, 0, 2] = -1
外積(Cross Product)
v1 = np.arange(3)
v2 = np.arange(3)+1
cross_numpy = np.cross(v1,v2))
cross_einsum = np.einsum('ijk,i,j->k', eijk, v1, v2))
行列式(Determinant)
- Levi Civitaを使った行列式の公式(3×3)
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])
※ ちょっとだけややこしい？でもスッキリはしているし、規則性があるから拡張も容易。
あぁ、ややこしきCramer's fomula
こいつを使うと、Levi-Civita記号を使って行列演算で逆行列を求められる
(4次元以上は高次のLevi-Civita記号への拡張が必要(後述))
A = np.arange(9).reshape(3,3)
A[0,0] = 1
inv_numpy1 = np.linalg.inv(A)
inv_numpy2 = np.linalg.solve(A, np.identity(3))
det_einsum = np.einsum('ijk,i,j,k',eijk,A[0],A[1],A[2])
inv_einsum = np.einsum('ijk,pqr,qj,rk->ip',eijk,eijk,A,A)/(2.0*det_einsum)
高階/高次への拡張
高階(足の数が多い)かつ高次元のテンソルに対しても、イプシロンを定義することにより、計算が可能(Generalized Edinton's Epsilon)
Wikipediaに書いている(エディントンのイプシロン)
これで任意の次元で上記式が使える*4
性能比較
速度、実際のところbuilt-inのnumpyと比べてどうなのか、ちょっと測ってみる
今回測るのは
- 内積計算 : (1,1,....,1)同士の内積を10000回計算した時の時間(1-100次元)
- 行列演算1 : 要素が1のN×N行列の積を10000回計算
- 行列演算2 : 要素が1のN×Nの行列(AB)Cを100回計算
- 外積 : (1,1,1)×(1,1,1) を10000回計算
- 行列式 : np.arange(9).reshape(3,3)の行列式を10000回計算
- 逆行列 : np.arange(9).reshape(3,3)の[0,0]成分だけ1に変換した行列の逆行列を10000回計算
また、測ったPCのスペックは
- OS : Mac OSX
- プロセッサ : 1.3GHz Intel Core M
- メモリ 8GB 1600MHz DDR3
弱いスペックのマシンなので、参考程度にどうぞ
numpyの方が2倍くらい早い。
%matplotlib inline
import matplotlib.pyplot as plt
import time
time_num = []
time_ein = []
d = 101
iteration = 10000
for N in range(d):
a = np.ones(N)
b = np.ones(N)
np_start = time.time()
for _ in range(iteration):
np.dot(a,b)
time_num.append(time.time()-np_start)
ein_start = time.time()
for _ in range(iteration):
np.einsum("i,i",a,b)
time_ein.append(time.time()-ein_start)
plt.plot(range(d),time_num,label="builtin")
plt.plot(range(d),time_ein,label="einsum")
plt.xlabel("dimension")
plt.ylabel("time")
plt.xlim(0,100)
plt.savefig("inner.png")
plt.legend()
行列演算1 : 要素が1のN×N行列の積
左が1-100次元。右が1-20次元。次元が高いところではnumpyの圧勝
行列演算2 : 要素が1のN×Nの行列の(AB)C(3つの積)
左が1-100次元。右が1-5次元。圧勝具合がさらにひどい
ちなみに3つの行列の積をeinsumでかくとこう
np.einsum("ij,jk,kl->il",a,b,c)
外積 : (1,1,1)×(1,1,1)
numpy : 0.21898913383483887
einsum: 0.04619288444519043
einsumに軍配。
行列式 : np.arange(9).reshape(3,3)の行列式
numpy : 0.14093804359436035
einsum: 0.06131696701049805
einsumに軍配。
逆行列 : np.arange(9).reshape(3,3)の[0,0]成分だけ1に変換した行列の逆行列
np.linalg.inv : 0.15855097770690918
np.linalg.solve : 0.27668190002441406
einsum : 0.45436620712280273
numpyに軍配
まとめと感想
- 高次元になると圧倒的に計算時間がかかる
- 低次元(d<10)程度なら2アバイ程度の計算時間
- 外積と行列式だけは組み込みより早い
- ただし、レビチビタをつかって行列演算をnumpyですると、ソッチのほうが早いかも。
高次の行列演算がバンバン出てくる機械学習などには少し不向きかなと思う。
ただ、
- 低次元では遜色ない計算量
- 直感的な表記での計算が可能
なことを考えると、物理系の演算で用いるのはいいかもしれないです。。
高階のテンソルになると、どの足を消すんだっけ？とか面倒なことになるので。
慣れている表記法でそのまま計算できるのは、ミスも少なくなりますしね。
何かの物理計算をちょこっと確かめる時とかに、使ってみてはいかがでしょうか。ではではっ！
を丁寧に教えて、「わかった！！」って言ってたのに、
になると解けなくなる。
