クロネッカー積

Wikipediaを見よ。
ざっくりいうと、2つの行列やベクトルの各要素の組み合わせを作っていく計算で、AxB行列とCxD行列のクロネッカー積はACxBD行列になって、繰り返しでクロネッカー積を取っていけばゴリゴリと行列が巨大になっていきます。

何に使うの? って感じですが、量子コンピュータをシミュレーションする上ではかなり使います。
複数の系をつなげて、ひとつの系として見るには、2つの量子状態のクロネッカー積をとることになるからです。

戦いに挑むプレイヤー

彼らで、パフォーマンス天下一武闘会をやっていきます。

プレイヤー	説明
普通の2次元配列 (np.array)	ぼくら大好きnumpyの2次元配列
対角要素だけの1次元配列(np.array)	対角行列同士のクロネッカー積は対角行列なので、対角要素だけ1次元配列で取り出せばよくないですか?
1次元配列(np.array) + `scipy.linalg.kron`	numpyに`np.kron`関数が用意されているんですが、scipyにも似た関数があるので、そっちも使ってみました
`scipy.sparse`のdia行列	本命登場。scipyの疎行列のひとつ。対角圧縮格納行列という形式
`scipy.sparse`のcoo行列	scipyの疎行列のひとつ。座標とデータを格納する方式
`scipy.sparse`のlil行列	scipyの疎行列のひとつ。リストのリストとしてデータを格納
`scipy.sparse`のdok行列	scipyの疎行列のひとつ。dictでデータを格納

疎行列が一杯出てきますね。
「0がいっぱいある行列って、持つの無駄じゃん」ってことで、それをなんとかしたやつが疎行列なんですが。こいつら、種類が多すぎて大変なんですよ。

結果

つまり、2x2行列とのクロネッカー積をとっていって、行列の次元を倍々にしてみましょう。

from time import perf_counter
from contextlib import contextmanager
import numpy as np
import scipy.linalg
from scipy.sparse import eye, kron, dia_matrix, coo_matrix, lil_matrix, dok_matrix

@contextmanager
def timer(ls):
    t0 = perf_counter()
    yield
    ls.append(perf_counter() - t0)

fmt = ['dia', 'coo', 'lil', 'dok']
mat1 = [np.eye(2), np.ones(2), np.ones(2).reshape((-1, 1))]
for f in fmt:
    mat1.append(eye(2, format=f))
mat2 = [x.copy() for x in mat1]

print('n\tdim\tnp.kron (mat)\tnp.kron (diag)\tscipy.linalg.kron\tdia\tcoo\tlil\tdok')
for i in range(15):
    time = []
    with timer(time):
        mat1[0] = np.kron(mat1[0], mat2[0])
    with timer(time):
        mat1[1] = np.kron(mat1[1], mat2[1])
    with timer(time):
        mat1[2] = scipy.linalg.kron(mat1[2], mat2[2])
    for j, f in enumerate(fmt, start=3):
        with timer(time):
            mat1[j] = kron(mat1[j], mat2[j], format=f)
    print(i, len(mat1[1]), *time, sep='\t')

こんなコードを動かしてみました。

結果、こうなりました

dim     np.kron (mat)   np.kron (diag)  scipy.linalg.kron       dia     coo     lil     dok
4       7.621300028404221e-05   4.2740000935737044e-05  1.618100213818252e-05   0.00034239199885632843  0.00010936500621028244  0.00040860599983716384       0.0002625720007927157
8       3.6978999560233206e-05  2.091900387313217e-05   1.4868004654999822e-05  0.00031008200312498957  0.00010704099986469373  0.00037807899934705347       0.0002334279997739941
16      3.546700463630259e-05   1.7302001651842147e-05  1.6571000742260367e-05  0.000312356001813896    0.00011324300430715084  0.0003890799998771399        0.000256961997365579
32      5.516300007002428e-05   2.5807996280491352e-05  2.214199776062742e-05   0.00030120500014163554  0.00010780199954751879  0.0004031159987789579        0.0005835240008309484
64      0.00010105999535880983  2.7360998501535505e-05  2.941599814221263e-05   0.00030793799669481814  0.0001086229967768304   0.00042860399844357744       0.0002643059997353703
128     0.00020814999879803509  4.005499795312062e-05   5.749799311161041e-05   0.00031550200219498947  0.00010875400039367378  0.0005037450027884915        0.00027479499840410426
256     0.0007497959959437139   6.817800021963194e-05   7.954899774631485e-05   0.0003335859946673736   0.00011099899711553007  0.0014838429997325875        0.0003385250020073727
512     0.00291782299609622     0.00013850899995304644  0.00015007200272521004  0.0004066420005983673   0.00011842299863928929  0.000978336000116542 0.00039470100455218926
1024    0.009874209994450212    0.0002527940014260821   0.0002898239981732331   0.0010409530004835688   0.0001297330018132925   0.0016027960009523667        0.0006418839984689839
2048    0.0948937740031397      0.0004355569981271401   0.000531908000994008    0.0004898590050288476   0.00015769599849591032  0.0030208060052245855        0.0009902669989969581
4096    0.4371686599988607      0.0008414180047111586   0.0010229589970549569   0.0006520530005218461   0.00018501699378248304  0.005185828005778603 0.0017490499958512373
8192    1.5510790590051329      0.0017292539996560663   0.0019197299989173189   0.00098355500085745     0.00034314399817958474  0.014461956001468934 0.002987473999382928
16384   5.583846742003516       0.0033475189993623644   0.004184869001619518    0.0011626810010056943   0.0003905220000888221   0.01908096900297096  0.006053636003343854

メモリがだいぶつらくなってきました。
要素数が万のオーダーになってくると、2次元配列で持つのがつらくなります。
脱落させて次を見てみましょう。

dim     np.kron (diag)  scipy.linalg.kron       dia     coo     lil     dok
## 略
16384   0.003250737994676456    0.0040016369966906495   0.0011192689999006689   0.0003866760016535409   0.01951057599944761     0.005936435001785867
32768   0.006436383002437651    0.008126152002660092    0.00257632300053956     0.0006475750051322393   0.04585514999780571     0.012816850001399871
65536   0.013466698997945059    0.016394651000155136    0.0036876980011584237   0.0012796200026059523   0.09913738199975342     0.027271613995253574
131072  0.027018057997338474    0.03451431899884483     0.007149820994527545    0.002364355997997336    0.20790910699724918     0.07132508500217227
262144  0.054609892998996656    0.07083234099991387     0.013791570003377274    0.004519036003330257    0.4256131450019893      0.15278254599979846
524288  0.11130649899860146     0.14217567000014242     0.029465308994986117    0.008979752994491719    0.8125153309956659      0.3088614529988263
1048576 0.22045852500014007     0.2833445610012859      0.060802528001659084    0.018181282001023646    1.6133216590023949      0.6439245029978338
2097152 0.4423036929947557      0.5692931669982499      0.13147616000060225     0.037687891999667045    3.204140270005155       1.342386914002418
4194304 0.8839434900000924      1.1379078309983015      0.265317548000894       0.07555588000104763     6.303365648003819       2.8236029719992075
8388608 1.7452782509935787      2.242528515998856       0.5512261999974726      0.14783157800411573     12.572135554000852      5.75408922600036
16777216        3.614834731000883       4.505664156000421       1.1048625590046868      0.3020367110002553      25.316714471002342      11.952817262994358

lilとdokが厳しくなってきました。
脱落させて、他やっていきましょう。

dim     np.kron (diag)  scipy.linalg.kron       dia     coo
## 略
16777216        3.6877648160007084      4.646497107998584       1.0930974110015086      0.29926496299594874
33554432        7.348309517001326       9.31613551099872        2.176151612002286       0.591943232997437
67108864        14.420524740002293      18.4003505590008        4.722769480998977       1.2042434280010639

多分、numpy.kronがいちばんシンプルなことやるだけで済むはずなんですが、遅いです。
あと、scipy.linalg.kronも、numpy.kronとの違いがわからないので使ってみましたが、あんまりパフォーマンス上のメリットはなさそうです。

もう結果は見えてきてますが、決勝戦やりましょう。
diaは対角圧縮という名前がついていて、対角行列に強いんじゃないかと期待したんですが、残念ながら勝てるか怪しいですね。

dim     dia     coo
## 略
67108864        4.889362873000209       1.2848818929996924
134217728       11.382003072998486      2.5668652079984895

時間はもう少し頑張れそうですが、メモリが厳しくなってきました。
cooがだいぶ早いです。

優勝: coo行列

coo行列がだいぶ速い、という結果になりました。
正直、理由はよくわからないですが、ライブラリのうちCで書かれている比率が高いとか、そういうようなことだと思っています。