数理解析研究所講究録 1295 巻 2002 年 213-219

213

Risa/Asir の Matrix 演算の新しい実装について
村尾裕一 \dag er
電気通信大学兵頭礼子 *

AlphaOmega Inc.

齋藤友克\ddagger

AlphaOmega lnc.

Abstract
近年数式処理の分野においても, 高速算法 (例えば FFT など) の利用は常識となっている. しかし ,では, 行列の演算には古典的計算方法が取られており, 特に実に高速化ははかられていない. 今回, 数値計算分野のアルゴリズムを$64\cross
64$
倍強の速度で乗算行列において 3
装したところ, 5 次多項式を各要素に持つを行うことができた. 他のサイズの行列でも, 古典的アルゴリズムよりも高速に乗算を行うことができる. しかし , 理論的な速度の向上と実際は異なっており今後の課題$\mathrm{R}\mathrm{i}\mathrm{s}\mathrm{a}/\mathrm{A}\mathrm{s}\mathrm{i}\mathrm{r}$$\mathrm{R}\mathrm{i}\mathrm{s}\mathrm{a}/\mathrm{A}\mathrm{s}\mathrm{i}\mathrm{r}$

である.

1

$\mathrm{R}\mathrm{i}\mathrm{s}\mathrm{a}/\mathrm{A}\mathrm{s}\mathrm{i}\mathrm{r}$の行列乗算の実情

は, 行列の積は定義にしたがった行列の内積演算によるアルゴリズムがの行列 , と $m\cross n$ の行列 $B$ の乗算を行う場合, 積の回数は実装されている.
の正方行列の積の, 和は $l\cdot n(m-1)$ 回発生する. よって, 計算量は楊合, $O(n^{3})$ である.
$\mathrm{R}\mathrm{i}\mathrm{s}\mathrm{a}/\mathrm{A}\mathrm{s}\mathrm{i}\mathrm{r}$$A$$l\mathrm{x}m$$n\mathrm{x}n$$l\cdot n\cdot m\fbox\Pi$

2

今回

2.1

Strassen-Winograd アルゴリズム$l\mathrm{x}m$$\mathrm{R}\mathrm{i}\mathrm{s}\mathrm{a}/\mathrm{A}\mathrm{s}\mathrm{i}\mathrm{r}$

行列 $A,$$m\mathrm{x}n$に実装した方法

行列 $B$ において $A\cdot B$ の演算を行う場合, 行列 $A,$ $B$ を$A=(\begin{array}{ll}A_{11} A_{12}A_{21} A_{22}\end{array})$,$B=(\begin{array}{ll}B_{11} B_{12}B_{21} B_{22}\end{array})$とすると$A_{11}B_{11}+A_{12}B_{21}$$A\cdot B=\{$$w+u+A_{22}(B_{21}-t_{2})$$\overline{l\mathrm{h}\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{k}\mathrm{o}@\mathrm{a}2\mathrm{z}.\mathrm{c}\mathrm{o}.\mathrm{j}\mathrm{p}}$\dagger murao@cs .uec.ac.jp

kaitO@a2z .co.jp

XXXII-I$w+v+(A_{12}-s_{2})B_{22}w+u+v)$214$s_{1}$$=$$A_{21}+A_{22}$$s_{2}$$=$$s_{1}-A_{11}$$t_{1}$$=$$B_{12}-B_{11}$$t_{2}$$=$$B_{22}-t_{1}$

$u$

$=$

$(A_{11}-A_{21})(B_{22}-B_{12})$

$v$

$=$

$w$

$=$

$=$

$-A_{11}+A_{21}+A_{22}$$=$$B_{11}-B_{12}+B_{22}$$s_{1}t_{1}$$=$$(A_{21}+A_{22})(B_{12}-B_{11})$$A_{11}B_{11}+s_{2}t_{2}$$=$$A_{11}B_{11}+$$(-A_{11}+A_{21}+A_{22})(B_{11}-B_{12}+B_{22})$となる.

Strasse-Winograd の方法での計算量は, 次数が半分の行列の乗算が 7 回, 加算が
回発生する
.
15
$m2^{k}$の行列の計算量を $T(m, k)$ とすると,次数$T(m, k)$(次数 $m2^{k-1}$ の行列の和 : $m2^{2(k-1)}t_{+}$ )$=$$7\cross T(m, k-1)+15\mathrm{x}$$=$$7^{k}\mathrm{x}T(m,0)+15m^{2}(2^{2k}-7^{k})/(2^{2}-7)t+$$=$$7^{k}(m^{3}t_{*}+m(m-1)t_{+})-5m^{2}(2^{2k}-7^{k})t_{+}$

よって,$O(7^{k})\approx O(7^{\log_{2}n})=O(n^{\log_{2}7})\approx O(n^{2.8})$となる.

このアルゴリズムを, 分割した行列のサイズが 4 行 4 列以下になるまで適用し, 4
行 4 列以下になったところで古典的アルゴリズムを適用して, 演算を行う. 行列が奇数行奇数列であった場合は, 偶数行偶数列の行列になるように要素が全て 0 の行または列を padding して演算を実行するよう実装した.

3

計算速度の比較

3.1

実験データ

計算量の比較のため次の行列を考える. 行列 $A,$ $B$ は大きさとする.
あり各成分はおのおの$a_{ij},$Case$n\mathrm{x}n$の正方行列で$b_{:,j}$$\mathrm{I}$

$\{$

$aij$

$=$

$(i+1)(j+1)(x^{5}+x^{4}+x^{3}+x^{2}+x)$ ,$b_{ij}$$=$$(i+1)(j+1)(x^{5}+x^{4}+x^{3}+x^{2}+x)$$1\leq i,j\leq n$

Case II$a_{ij}$$\{$$b_{ij}$$=$$(i+1)(j+1)(x^{5}+x^{4}+x^{3}+x^{2}+x)$ ,$=$$(i+1)(j+1)(y^{5}+y^{4}+y^{3}+y^{2}+y)$$1\leq i,j\leq n$とする. また実験は, AMD athron$1400\mathrm{M}\mathrm{h}\mathrm{z}(\mathrm{O}\mathrm{S}:\mathrm{F}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{B}\mathrm{S}\mathrm{D}4.4)$

算時間の計測は秒単位である.

XXXII-2

によって行った. 演

215
3.1.1

Case I の演算時間

古典的アルゴリズム (標準) と Strassen-Winograd アルゴリズムで Case I の場合の $A\cross B$ を実行し, 演算時間 (単位は秒) を計測 (表 1) する.

表 1: 演算時間の比較 (Case I)

lb

Strassen-Winograd
0.01526
0.01598
0.10260.181
0.5701
1.137
2.761
10.24
12.83
128 94.78
60.18
256 864.1

Size
8
16
32
64

$\mathrm{F}_{\overline{\mathrm{t}}}^{\sqrt[\backslash ]{}}\tau\not\leqq\backslash ffi\backslash \backslash$

1.047
1.764
1.994
3.709
7.387
14.36

3.12Case 垣の演算時間
Case I と同様に Case II の場合の古典的アルゴリズムと Strassen-Winograd ア
ルゴリズムで $A\cross B$ の演算を実行し , 計算速度を計測 (表 2) する. 表 1, 表 2 から,
表 2: 演算時間の比較 (Case

Size
8
16
32

$\ovalbox{\tt\small REJECT}^{\sqrt[\backslash ]{}}\backslash \#\backslash \mathrm{a}\backslash$$\mathrm{I}\mathrm{I}$)Strassen-Winograd

lb

0.011780.08877
0.5214
2.684
13.03
66.13

0.97961.245

0.011540.1105
1.025
9.235

64
128 75.95
256 953.1

1.966

3.441
5.829
14.41

Strassen-Winograd アルゴリズムを適用した場合, 古典的アルゴリズムで演算を行っ, 計算量の比とは明た場合よりも演算時間が短くなっている. しかし, 演算時間の$1\mathrm{t}fi\backslash ^{\backslash }$らかにかけ離れた値になっている. よって, 計算量以外の要因が, 演算時間の短縮につながっていると考えられる.

4

演算回数, 計算量の比較

4.1

Case I の場合の演算回数演算の高速化が計算量以上にはかられたため, 計算量以外の要因を考える. Case Iの場合の行列の乗算を行う際に, 行列要素どうしの演算の加算, 乗算の発生回数をカウントした. 以下の表 3 に示す. 行列要素どうしの演算回数では, 演算回数は古典的アルゴリズムに比べ減少しているものの, 計算量の理論通りの減少率となり, 古典的アルゴリズムと Strassen-Winograd アルゴリズムの演算速度の差につな力 るような関$1^{\backslash }$係は見付からなかった.
XXXII-3

216(Case I)$\mathrm{g}3:4\overline{\mathrm{T}}\partial’|\rfloor\ovalbox{\tt\small REJECT}\ovalbox{\tt\small REJECT}_{\backslash }\epsilon^{\backslash ^{\backslash }}\vee\grave{J}\mathrm{b}a)^{\backslash }\{,\ovalbox{\tt\small REJECT}\backslash \ovalbox{\tt\small REJECT}\Gamma\fbox \mathrm{J}\ovalbox{\tt\small REJECT}$

Size$\mathrm{F}\overline{\mathrm{T}}\epsilon\backslash \sqrt[\backslash ]{}ff\backslash ]7J\triangleright=\grave{\mathrm{f}}$$\not\supset\Pi\ovalbox{\tt\small REJECT}$

8
16

32
64
128
256

$|)$

$\ovalbox{\tt\small REJECT}\ovalbox{\tt\small REJECT}$

448

Strassen-Winograd$\lambda^{\theta}\mathrm{A}$512$\mathrm{A}^{-}\mathrm{p}_{\mathrm{Q}}^{\equiv+}$

960$\not\supset \mathrm{O}\ovalbox{\tt\small REJECT}$

624
5520

4096
7936
3840
32768
64512
31744
43248
2621442580485201923211682080768 2097152 4177920 2321904
16711680 16777216 33488896 16548240$\ovalbox{\tt\small REJECT} \mathrm{g}$$\mathrm{r}_{\equiv+}^{-}\mathrm{R}\mathrm{D}$

448
1072
3136
8656
21952
65200
153664474832
1975648 4297552
7529536 24077776

そこで, さらに行列要素の各項レベルでの演算回数を調べる. つまり, 加算が発生
した場合は, 加算した結果の項の数, 乗算の場合は掛け合わせる各要素の項の数を乗算した数をそれぞれカウントし, 実際の項レベルの演算回数 (表 4) を調べる. 各サイ

表 4: Strassen-Winograd アルゴリズムの項レベル演算回数 (Case I)Size$7J\triangleright\grave{-}$$|)$$X\mathrm{A}$$\mathrm{F}^{\sqrt[\backslash ]{}}\overline{*}-\not\in\backslash \backslash \#\backslash$

8

Strassen-Winograd$\Phi--\sqrt\subsetneqq\backslash \backslash \not\in\backslash$

16

Strassen-Winograd$\ovalbox{\tt\small REJECT}_{-}^{\sqrt}\backslash \not\in\backslash \ae\backslash$

32

Strassen-Winograd$ff_{\overline{\mathrm{s}}}^{\sqrt[\backslash ]{}}\neg.\not\leqq\backslash \ae\backslash$

64

Strassen-Winograd$\mathrm{J}\coprod\Leftrightarrow$

4608
5680
36864
39712
29491223684823592961290160

18874368
Strassen-Winograd
6641776
150994944
Strassen-Winograd 32973392$\otimes-\backslash \#\sqrt[\backslash ]{}\backslash \ae\backslash$128$\mathrm{F}^{\sqrt[\backslash ]{}}\mathrm{U}\backslash -\not\in\backslash \backslash \#\backslash$$\circ \mathrm{K}a$$A\cdot\cup$$\ovalbox{\tt\small REJECT}\ovalbox{\tt\small REJECT}$

A
B

$-$
$\equiv+$$1\mathrm{b}$$\mathrm{D}$

12800
17408
11200
16880 1.031
10240013926462400
102112 1.364
1114112819200302400539248 2.066
6553600
89128961339200
2629360 3.340
52428800 71303168
5572800 12214576 5.838
419430400 570425344
22161600 55134992 10.35ズの行列の乗算で, 旧アルゴリズムと Strassen-Winograd アルゴリズムの演算回数の比は, 8 $\mathrm{x}8,16\cross 16,32\cross 32,64\cross 64,128\cross 128,256\cross 256$ 行列の順に, 1031, 1364,

2066, 3340, 5838, 1035 となる.

4.2

Case II の場合の演算回数

Case 垣の場合の, 古典的アルゴリズムと Strassen-Winograd アルゴリズムでの演算を実行し , Case I の場合と同様に演算回数を計測 (表 5) し比較する. 各サイズの行列の乗算で, 標準的アルゴリズムと Strassen-Winograd アルゴリズムの計算量の$A\mathrm{x}B$

比は, $8\cross 8,16\cross 16,32\cross 32,64\cross 64,128\cross 128,256\cross 256$ 行列の順に, 09907, 1266,
1.859, 2.965, 4.990, 8.693 となる.
XXXII-4

217$\ovalbox{\tt\small REJECT} 5$

Size:Strassen-Winograd$7J\triangleright=\grave{l}$$|)$$\mathrm{X}^{\grave{\backslash }}\mathrm{A}$$\mathrm{F}\overline{\tau}^{\sqrt[\backslash ]{}^{\backslash }}\mathrm{E}\backslash \ae\backslash$

8

Strassen-Winograd$\Phi_{-}^{-\sqrt[\backslash ]{}^{\backslash }}\not\in\backslash ffi\backslash$

16

Strassen-Winograd$ff_{\mathrm{T}-}^{\sqrt[\backslash ]{}}\backslash \not\in\backslash \ae\backslash$

32

Strassen-Winograd$\Phi-\not\leqq\sqrt[\backslash ]{}^{\backslash }\backslash \ae\backslash$

64

Strassen-Winograd$\Phi-\backslash -\sqrt{}^{\backslash }\not\in^{\backslash }\backslash \#$128)$\vee 7J\triangleright\supset^{\backslash }|$$l\square \Leftrightarrow$

12800
14640

10240099360
8192005791206553600
308088052428800
15442800$i\mathrm{X}^{\grave{\backslash }}\mathrm{A}C^{)}\grave{\grave{;}},\mathrm{F}\Leftrightarrow\fbox \mathrm{c}\mathrm{J}\ovalbox{\tt\small REJECT}$$\ovalbox{\tt\small REJECT}\Leftrightarrow$

12800
11200
10240062400
8192003024006553600
1339200
52428800
5572800(Case$\mathrm{I}\mathrm{I}$$\infty\square \mathrm{n}\equiv-+$

)

kb

25600
25840 0.991

204800161760 1.266
1638400
881520 1.859
13107200
4420080 2.965
104857600

21015600 4.990
419430400 419430400 838860800
Strassen-Winograd 74336080 22161600 96497680 8.693

Strassen-Winograd$\mathrm{F}_{\tau\backslash -}^{\sqrt[\backslash ]{}}.\not\in^{\backslash }\backslash ffi$$\circ\propto a$$4\cup \mathrm{U}$以上 Case , Case 垣の結果から, 行列の各要素の項レベルでの演算回数が演算時$\mathrm{I}$間の比に比較的近いものとなり, 項レベルの演算回数が演算の高速化に影響を与えるものと考えられる.

4.3

分割終 7 条件

Strassen-Winograd アルゴリズムは, 行列を 4 分割し再帰的に Strassen-Winograd
アルゴリズムを呼び出し , 実行している. そこで, 各要素に Case I の場合の $128\cross 128$行列で, Strassen-Winograd アルゴリズムから古典的アルゴリズムに切り替わるサイズを 64,32,16,8,4 と変化させた場合の演算速度と演算回数を比較し , 分割を終了する
サイズを決定する. Strassen-Winograd アルゴリズムと旧アルゴリズムとの演算時間
表 6: 分割サイズによる計算量と演算速度の比較 (Case I)

Size
64
32
16
8
4

$\not\supset\Pi\Leftrightarrow$$\ovalbox{\tt\small REJECT}\ovalbox{\tt\small REJECT}$

18874368
16936960
12605440
90150407089088664177652428800
45875200
31948800
19353600
10713600
5572800$\bigwedge_{\mathrm{R}^{\frac{-}{\mathrm{Q}}+}}^{=}$$\mathrm{P}_{J}^{\rfloor\backslash }$71303168 102. 1
62812160 81.57
44554240 53.29
28368640 32.6
17802688 19.53
12214576 12.73の比 (表 6) は $64\cross 64,32\cross 32,16\cross 16,8\cross 8,4\cross 4$ に分割サイズを変化させるにつれて, 088, 071, 064, 063, 0686 となる. また, 演算回数の比 (図 3) は同様に, 080,
065, 061, 060, 065 となる.

XXXII-5218
行列に分割するまで滅少し続けており, また, 演算時間も減少し演算回数は,
分割は行列になるまで実行することとした.
続けている. よって,
$4\cross 4$$4\cross 4$5まとめに実装されていた, 行列の定義に基づく行列の乗算とStrassen-Winograd アルゴリズムの演算時間の比較を行った. これまでの計算量の理
今回, これまで$\mathrm{R}\mathrm{i}\mathrm{s}\mathrm{a}/\mathrm{A}\mathrm{s}\mathrm{i}\mathrm{r}$論値を越えた高速化が計られており, 直接の理由については今だ不明である. しかし,行列の要素に含まれる項レベルの演算量の比較から, 多項式の項レベルの計算量力吠
きな要因となっていることが推測できる. 数式処理の場合では, 従来の計算量ではな
くメモリアクセス等の空間計算量のような新たな尺度を用いる必要があるのではな
いかと思われる.

参考文献[Bailey, 1988] Bailey, D. H. (1988). Extra high speed matrix multiplication on the
Cray-2. SIAM J. Sci. Stat. Comp., $9(3):603-607$ .
[Brent, 1970] Brent, R. P. (1970). Algorithms for matrix multiplication. Technical, Stanford University.
Report[Bunch and Hopcroft, 1974] Bunch, J. R. and Hopcroft, J. E. (1974). Triangularfactorization and inversion by fast matrix multiplication. Mathematics of Computation, 28(125).
[Cohen and Roth, 1976] Cohen, J. and Roth, M. (1976). On the implementation ofStrassen’s fast multiplication algorithm. Acta Infomatica, 6:341-355.
[Coppersmith and Winograd, 1990] Coppersmith, D. and Winograd, S. (1990). Matrix multiplication via airthmetic progressions. Journal of Symbolic Computation,$\mathrm{C}\mathrm{S}157$

9:251-280.

[Cray SciLib, 1990] Cray SciLib (1990). (UNICOS) Math and Scientific Library, Release 6.0. Cray Research Inc.
Reference Manual[Douglas et al., 1994] Douglas, C. C., Heroux, M., Slishman, G., and Smith, R. M.
(1994). GEMMW: Aportable level 3BLAS Winograd variant of Strassen’s matrix
matrix multiply algorithm. Journal of Computational Physics, 10:1-10.
[Higham, 1990] Higham, N. J. (1990). Exploiting fast matrix multiplication within
the leve1-3 BLAS. $ACM$ Transactions on Mathematical Software, 16:352-368.
] Huss-Lederman, S., Jacobson, E. M., Johnson, J. R.,[Huss-Lederman et al.,Tsao, A., and Turnbull, T. (1996a). Implementation of Strassen’s algorithm for
matrix multiplication. In Supercomputing ’96 Conference Proceedings.
] Huss-Lederman, S., Jacobson, E. M., Johnson, J. R.,[Huss-Lederman et al.,Tsao, A., and Turnbull, T. (1996b). Strassen’s algorithm for matrix multiplication: Modeling, analysis, and implementation. Technical Report CCS-TR-96-147,
Center for Computing Sciences.
[IBM ESSL, 2000]
ESSL (-2000). Engineering and Scientific Subroutine Library: Guide and Reference. IBM.
.
.rs6000.
http:$(SR- \mathit{2}\mathit{0}\mathit{8}\mathit{1})$$1996\mathrm{a}$$1996\mathrm{b}$$\mathrm{I}\mathrm{B}\mathrm{M}$$//\mathrm{m}\mathrm{w}$$\mathrm{i}\mathrm{b}\mathrm{m}.\mathrm{c}\mathrm{o}\mathrm{m}/\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{o}\mathrm{u}\mathrm{r}\mathrm{c}\mathrm{e}/\mathrm{a}\mathrm{i}\mathrm{x}_{-}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{o}\mathrm{u}\mathrm{r}\mathrm{c}\mathrm{e}/\mathrm{s}\mathrm{p}-\mathrm{b}\mathrm{o}\mathrm{o}\mathrm{k}\mathrm{s}/\mathrm{e}\mathrm{s}\mathrm{s}1/$

XXXII-6

219[Pan, 1978] Pan, V. Y. (1978). Strassen algorithm is not optimal. trilinear techniqueof aggregating, uniting and canceling for constructing fast algorithms for matrix
multiplication. In Proc. 19th FOCS, pages 166-176.
[Strassen, 1969] Strassen, V. (1969). Gaussian elimination is not optimal. Numer.
Math., 13:354-356.
[Strassen, 1986] Strassen, V. (1986). The asymptotic spectrum of tensors and theexponent of matrix multiplication. In Proc. 27th FOCS, pages 49-54.
matrices. Linear[Winograd, 1971] Winograd, S. (1971). On multiplication of
Algebra and its Applicahons, 4:381-388.
$2\cross 2$

XXXII-7