プロクラシスト

今日の寄り道 明日の近道

【注目記事】 Pythonでゼロから機械学習/データ分析を学ぶためのサイトマップ

【einsum】アインシュタインの縮約記法のように使えるnumpyの関数。性能と使い方を解説。

技術メモ 技術メモ-python

スポンサーリンク

f:id:imslotter:20170508210354p:plain

こんにちは、ほけきよです。

大学で物理*1を嗜んだ方ならわかるであろうEinsteinの縮約記号。 計算の上で色々省略できるしとにかく慣れれば色々便利な記法です！

物理学者以外には馴染みがなく微妙かもしれませんが、「便利そうだな〜〜」と思って使ってみたり試してみたりしたので、メモとしてまとめておきます。誰かの参考になれば幸いです。

einsum

Einstein(アインシュタイン)の縮約記法風に足し合わせができるnumpyの関数

縮約記法とは

たとえば、行列 A=a_{ij} と B=b_{jk} っていうのがあって、行列の演算をしたいときバカ正直に書くと AB=\sum_{j=i}^N a_{ij}b_{jk}という風に計算される。 こういう演算で、\Sigmaが大量に出てきすぎるので、「じゃあ略してしまえ」となった。

ルールは

  • 同じ記号のところは足し合わせる。

これだけ。さっきの足(a_{i}のiの部分)が

便利なのは、足の数が3つ以上ついている*2テンソル計算も簡単にできるということ。 物理学の世界(相対論など)では、この高次テンソルでの計算がたくさん出てくるため

「そんなに\Sigma書いてられるか！」

ということで、略される。

例

  • 行列演算(行列A=a_{ij}, B=b_{jk}) -> c_{ik} = a_{ij}b_{jk}
  • ベクトルa_{i}とb_{j}の内積 -> c = a_{i}b_{i}
  • 行列A_{ij}のトレース -> T = a_{ii}

文法

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) #60
np.einsum("ii->i",a) #[ 0,  6, 12, 18, 24]
np.einsum("ij->i",a) #[ 10,  35,  60,  85, 110]
np.einsum("ij->j",a) #[50, 55, 60, 65, 70]
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) #AB
np.einsum("ij,jk->ki",a,b)  #BA 足(ik)が変わったのに注意
np.einsum("ij,ij->ij",a,b) #いわゆるアダマール積

ベクトル演算その②

ベクトルの演算も、行列とほぼ同様にできる

v1 = np.arange(3)       #[0,1,2]
v2 = np.arange(3)+1   #[1,2,3]

np.einsum("i,i",v1,v2) #内積
np.einsum("i,j->ij",v1,v2) #直積

レビチビタ記号とその応用(行列式、逆行列、外積)

テンソル表記の時に気になるのが、行列式や逆行列、外積の計算ってどうするのよ？？ということ。 こういう時に便利なのがレビチビタ(Levi Civita)記号。定義は以下の通り

hogehoge

3次元以上空間の時には、この記号が計算をスッキリまとめられて便利。 これを使うと、einsumの枠組みでより多様な計算が実現できる

ただし、レビチビタ記号は調べた範囲では実装されてなさそう？*3なので自分でつくる

# using "Levi civita symbols" (As long as I search, there is no built-in function)
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)

c^{k} = \epsilon_{ijk}a_{i}b_{j}

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)

f:id:imslotter:20170508203507p:plain

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]) #0

※ ちょっとだけややこしい？でもスッキリはしているし、規則性があるから拡張も容易。

逆行列

あぁ、ややこしきCramer's fomula

こいつを使うと、Levi-Civita記号を使って行列演算で逆行列を求められる (4次元以上は高次のLevi-Civita記号への拡張が必要(後述))

# inverse(逆行列)
A = np.arange(9).reshape(3,3)
A[0,0] = 1
# numpyの実験コード http://d.hatena.ne.jp/sleepy_yoshi/20120513/p1
inv_numpy1 = np.linalg.inv(A) #numpy 
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倍くらい早い。 f:id:imslotter:20170508204445p:plain

コード*5

%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の圧勝 f:id:imslotter:20170508204718p:plain

行列演算2 : 要素が1のN×Nの行列の(AB)C(3つの積)

左が1-100次元。右が1-5次元。圧勝具合がさらにひどい f:id:imslotter:20170508204741p:plain

ちなみに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ですると、ソッチのほうが早いかも。

高次の行列演算がバンバン出てくる機械学習などには少し不向きかなと思う。 ただ、 ​ - 低次元では遜色ない計算量 - 直感的な表記での計算が可能 ​ なことを考えると、物理系の演算で用いるのはいいかもしれないです。。 高階のテンソルになると、どの足を消すんだっけ？とか面倒なことになるので。 ​ 慣れている表記法でそのまま計算できるのは、ミスも少なくなりますしね。 ​ 何かの物理計算をちょこっと確かめる時とかに、使ってみてはいかがでしょうか。ではではっ！ ​

*1:相対論とかテンソルとかかな

*2:この足の数が多いのを"高階"という。高次元とは別

*3:あったら教えて欲しい...

*4:しかし、4次元の外積ってどうなるんだ？？

*5:この後ほとんど同じコードなので、省略

imslotter

関連記事
スポンサーリンク


PROCRASIST をフォローする
'); $(".shwbtn").click(function(){ var $this = $(this); if($Contents.css('display') == 'none'){ $Contents.slideDown(400), $this.text("目次を隠す"); }else{ $Contents.slideUp(400), $this.text("目次をみる") }; }); });
コメントを書く

PROCRASIST