ここでは、量子計算をアインシュタインの縮約記法を用いて計算をします。そのための基礎知識などを少しずつ備忘録として書いていきます。
*どうやらnumpyは52までしか文字コードに対応していないので、opt_einsumを代わりに使う必要がありました。
参考： 量子コンピュータの計算をシミュレーションする際には、通常状態ベクトルや量子ゲートを使います。アインシュタインの縮約記法を使って量子回路を計算する際には、途中でテンソルなどが出てきます。解き方を確認します。参考ページから下記のような例題を持ってきました。
import numpy as np
#腕が4本と3本のテンソルを二種類作ります。
A = np.random.rand(1, 2, 2, 2)
B = np.random.rand(2, 2, 2)
#これらのテンソルを縮約してCを作ります。
C = np.einsum('ijkl,klm->ijm', A, B)
#プリント
print ('A:', A)
print ('B:', B)
print ('C:', C)
A: [[[[0.33851206 0.13251935]
[0.46585545 0.04285461]]
[[0.4669879 0.15598267]
[0.59812002 0.46526585]]]]
B: [[[0.23390734 0.85684659]
[0.87213317 0.95330902]]
[[0.00741036 0.94095885]
[0.77381743 0.83828137]]]
C: [[[0.23136878 0.89065983]
[0.60973267 1.50166669]]]
このように腕の縮約によって新しいテンソルができました。これらテンソルの縮約はFLOPカウントと呼ばれる計算量をキチンと測定できるようです。行列Mとベクトルvを使って試してみたいと思います。optimizeアルゴリズムはgreedyと書いてありますので、基本的にはヒューリスティクスで探索するようです。
import numpy as np
M = np.random.rand(2, 2)
v = np.random.rand(2)
path_info = np.einsum_path('ij,j->i', M, v, optimize='greedy')
print(path_info[0])
print(path_info[1])
['einsum_path', (0, 1)]
Complete contraction: ij,j->i
Naive scaling: 2
Optimized scaling: 2
Naive FLOP count: 8.000e+00
Optimized FLOP count: 9.000e+00
Theoretical speedup: 0.889
Largest intermediate: 2.000e+00 elements
--------------------------------------------------------------------------
scaling current remaining
--------------------------------------------------------------------------
2 j,ij->i i->i
path_info[1]を見ることで、scaling / FLOP count / Theoretical speedupなどが表示されています。正直なんだかわからないですが、数字を見ていけばよいようです。。今のところスケーリングは腕の本数に対応しているように見えます。
次は腕の本数が全部で5本です。
import numpy as np
A = np.random.rand(1, 2, 2, 2)
B = np.random.rand(2, 2, 2)
path_info = np.einsum_path('ijkl,klm->ijm', A, B, optimize='greedy')
print(path_info[0])
print(path_info[1])
['einsum_path', (0, 1)]
Complete contraction: ijkl,klm->ijm
Naive scaling: 5
Optimized scaling: 5
Naive FLOP count: 3.200e+01
Optimized FLOP count: 3.300e+01
Theoretical speedup: 0.970
Largest intermediate: 4.000e+00 elements
--------------------------------------------------------------------------
scaling current remaining
--------------------------------------------------------------------------
5 klm,ijkl->ijm ijm->ijm
import numpy as np
np.random.seed(123)
a = np.random.rand(2, 2, 2)
b = np.random.rand(2, 2, 2, 2)
c = np.random.rand(2, 2, 2)
d = np.random.rand(2, 2)
path_info = np.einsum_path('ijk,jlmn,klo,mo->in', a, b, c, d, optimize='greedy')
print(path_info[0])
print(path_info[1])
['einsum_path', (2, 3), (1, 2), (0, 1)]
Complete contraction: ijk,jlmn,klo,mo->in
Naive scaling: 7
Optimized scaling: 5
Naive FLOP count: 5.120e+02
Optimized FLOP count: 1.290e+02
Theoretical speedup: 3.969
Largest intermediate: 8.000e+00 elements
--------------------------------------------------------------------------
scaling current remaining
--------------------------------------------------------------------------
4 mo,klo->klm ijk,jlmn,klm->in
5 klm,jlmn->jkn ijk,jkn->in
4 jkn,ijk->in in->in
縮約の順番を工夫することで計算が速くなります。計算もoのエッジを最初に計算し、その次にmlを計算し、最後にjkを計算して答えが出ています。最初にbcを計算してしまうと、計算量のscalingが6となって増えてしまうので、順番に気を付ける必要があります。
現在NVIDIAのcuQunatumにはあるのかわかりませんが、縮約だけでなく分解を通じて計算量を減らせます。5ノードあり、1つがrank4のテンソルで腕が4本、ほかにベクトルが4つあります。
np.random.seed(123)
a = np.random.rand(2)
b = np.random.rand(2)
c = np.random.rand(2, 2, 2, 2)
d = np.random.rand(2)
e = np.random.rand(2)
path_info = np.einsum_path('i,j,ijkl,k,l', a, b, c, d, e, optimize='greedy')
print(path_info[0])
print(path_info[1])
['einsum_path', (0, 2), (0, 3), (0, 2), (0, 1)]
Complete contraction: i,j,ijkl,k,l->
Naive scaling: 4
Optimized scaling: 4
Naive FLOP count: 8.000e+01
Optimized FLOP count: 6.100e+01
Theoretical speedup: 1.311
Largest intermediate: 8.000e+00 elements
--------------------------------------------------------------------------
scaling current remaining
--------------------------------------------------------------------------
4 ijkl,i->jkl j,k,l,jkl->
3 jkl,j->kl k,l,kl->
2 kl,k->l l,l->
1 l,l-> ->
普通に考えれば最大で計算量はスケーリング4ですが、中心のテンソルを腕3本のテンソル二つに分解した後に計算をすると最適なパスで計算量が3におちます。
import numpy as np
np.random.seed(123)
a = np.random.rand(2)
b = np.random.rand(2)
c = np.random.rand(2, 2, 2)
d = np.random.rand(2, 2, 2)
e = np.random.rand(2)
f = np.random.rand(2)
path_info = np.einsum_path('i,j,imk,jml,k,l', a, b, c, d, e, f, optimize='greedy')
print(path_info[0])
print(path_info[1])
['einsum_path', (0, 2), (0, 1), (0, 2), (0, 1), (0, 1)]
Complete contraction: i,j,imk,jml,k,l->
Naive scaling: 5
Optimized scaling: 3
Naive FLOP count: 1.920e+02
Optimized FLOP count: 5.300e+01
Theoretical speedup: 3.623
Largest intermediate: 4.000e+00 elements
--------------------------------------------------------------------------
scaling current remaining
--------------------------------------------------------------------------
3 imk,i->km j,jml,k,l,km->
3 jml,j->lm k,l,km,lm->
2 km,k->m l,lm,m->
2 lm,l->m m,m->
1 m,m-> ->
このように、量子コンピュータの計算においてnumpy.einsumをつかって巨大な量子もつれを高速に計算できます。以上です。