Strassen アルゴリズム

Strassen's Algorithm for Fast Matrix Multiplication

上級

1. はじめに

$n \times n$ 行列の積 $C = AB$ を素朴に計算すると、各成分 $c_{ij} = \displaystyle\sum_{k=1}^n a_{ik} b_{kj}$ を求めるのに $n$ 回の乗算と $n-1$ 回の加算が必要であり、全体で $O(n^3)$ の計算量がかかる。

1969年、Volker Strassen は $2 \times 2$ 行列の積を 8回ではなく7回の乗算で計算する方法を示し、これを再帰的に適用することで計算量を

$$O(n^{\log_2 7}) = O(n^{2.807\ldots})$$

に削減した。これは「行列乗算に $O(n^3)$ は必要ない」ことを初めて示した画期的な結果であり、その後の高速行列乗算研究の出発点となった。

2. 素朴な行列乗算の計算量

$2 \times 2$ 行列の積を考える。

$$\begin{pmatrix} c_{11} & c_{12} \\ c_{21} & c_{22} \end{pmatrix} = \begin{pmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{pmatrix} \begin{pmatrix} b_{11} & b_{12} \\ b_{21} & b_{22} \end{pmatrix}$$

通常の定義では

$$c_{ij} = \displaystyle\sum_{k=1}^{2} a_{ik} b_{kj}$$

であり、8回の乗算と4回の加算が必要である。$n \times n$ 行列を $n/2 \times n/2$ のブロックに分割して再帰的に適用すると、漸化式

$$T(n) = 8\, T(n/2) + O(n^2)$$

が得られ、マスター定理 (Master theorem; 分割統治型の再帰式 $T(n) = a\,T(n/b) + f(n)$ の漸近解を $f(n)$ と $n^{\log_b a}$ の比較で与える定理) より $T(n) = O(n^{\log_2 8}) = O(n^3)$ となる。

3. Strassen のアイデア — 7回の乗算

Strassen は以下の 7つの補助積 $M_1, \ldots, M_7$ を導入した。

Strassen の7つの乗算

$$\begin{aligned} M_1 &= (a_{11} + a_{22})(b_{11} + b_{22}) \\ M_2 &= (a_{21} + a_{22})\, b_{11} \\ M_3 &= a_{11}\, (b_{12} - b_{22}) \\ M_4 &= a_{22}\, (b_{21} - b_{11}) \\ M_5 &= (a_{11} + a_{12})\, b_{22} \\ M_6 &= (a_{21} - a_{11})(b_{11} + b_{12}) \\ M_7 &= (a_{12} - a_{22})(b_{21} + b_{22}) \end{aligned}$$

これら7つの $M_i$ から、$C$ の各成分は加減算だけで復元できる。

$C$ の成分の復元

$$\begin{aligned} c_{11} &= M_1 + M_4 - M_5 + M_7 \\ c_{12} &= M_3 + M_5 \\ c_{21} &= M_2 + M_4 \\ c_{22} &= M_1 - M_2 + M_3 + M_6 \end{aligned}$$

3.1 正しさの検証

例えば $c_{11}$ を展開すると、

$$\begin{aligned} c_{11} &= M_1 + M_4 - M_5 + M_7 \\ &= (a_{11}b_{11} + a_{11}b_{22} + a_{22}b_{11} + a_{22}b_{22}) \\ &\quad + (a_{22}b_{21} - a_{22}b_{11}) \\ &\quad - (a_{11}b_{22} + a_{12}b_{22}) \\ &\quad + (a_{12}b_{21} + a_{12}b_{22} - a_{22}b_{21} - a_{22}b_{22}) \\ &= a_{11}b_{11} + a_{12}b_{21} \end{aligned}$$

となり、正しい結果が得られる。他の成分も同様に確認できる。

3.2 計算量

$n \times n$ 行列を $n/2 \times n/2$ のブロックに分割し、$M_i$ の計算で Strassen を再帰的に適用すると、乗算の回数が8回から7回に減るため漸化式は

$$T(n) = 7\, T(n/2) + O(n^2)$$

マスター定理より

$$T(n) = O(n^{\log_2 7}) = O(n^{2.807\ldots})$$

が得られる。加減算は18回に増えるが（素朴な方法では4回）、加減算は $O(n^2)$ であり、$n$ が大きくなると乗算回数の削減が支配的になる。

3.3 具体例: $2 \times 2$ 行列での Strassen

$A = \begin{pmatrix} 1 & 2 \\ 3 & 4 \end{pmatrix}$, $B = \begin{pmatrix} 5 & 6 \\ 7 & 8 \end{pmatrix}$ で Strassen を実行する。

Step 1: 7 つの補助積を計算

$$\begin{aligned} M_1 &= (a_{11}+a_{22})(b_{11}+b_{22}) = (1+4)(5+8) = 5 \cdot 13 = 65 \\ M_2 &= (a_{21}+a_{22})\, b_{11} = (3+4) \cdot 5 = 35 \\ M_3 &= a_{11}(b_{12}-b_{22}) = 1 \cdot (6-8) = -2 \\ M_4 &= a_{22}(b_{21}-b_{11}) = 4 \cdot (7-5) = 8 \\ M_5 &= (a_{11}+a_{12})\, b_{22} = (1+2) \cdot 8 = 24 \\ M_6 &= (a_{21}-a_{11})(b_{11}+b_{12}) = (3-1)(5+6) = 22 \\ M_7 &= (a_{12}-a_{22})(b_{21}+b_{22}) = (2-4)(7+8) = -30 \end{aligned}$$

Step 2: 結果を復元

$$\begin{aligned} c_{11} &= M_1 + M_4 - M_5 + M_7 = 65 + 8 - 24 - 30 = 19 \\ c_{12} &= M_3 + M_5 = -2 + 24 = 22 \\ c_{21} &= M_2 + M_4 = 35 + 8 = 43 \\ c_{22} &= M_1 - M_2 + M_3 + M_6 = 65 - 35 - 2 + 22 = 50 \end{aligned}$$

検算: $AB = \begin{pmatrix} 1\cdot5+2\cdot7 & 1\cdot6+2\cdot8 \\ 3\cdot5+4\cdot7 & 3\cdot6+4\cdot8 \end{pmatrix} = \begin{pmatrix} 19 & 22 \\ 43 & 50 \end{pmatrix}$ ✓

4. 再帰構造の図解

図1: Strassen アルゴリズムの再帰構造 n×n行列をn/2×n/2の4ブロックに分割し、7回の再帰的乗算を行う様子を示す図 A₁₁ A₁₂ A₂₁ A₂₂ A (n×n) × B₁₁ B₁₂ B₂₁ B₂₂ B (n×n) 7 回の乗算（再帰） M₁ = (A₁₁+A₂₂)(B₁₁+B₂₂) M₂ = (A₂₁+A₂₂) B₁₁ M₃ = A₁₁ (B₁₂−B₂₂) M₄ = A₂₂ (B₂₁−B₁₁) M₅ = (A₁₁+A₁₂) B₂₂ M₆ = (A₂₁−A₁₁)(B₁₁+B₁₂) M₇ = (A₁₂−A₂₂)(B₂₁+B₂₂) 18 回の加減算で C を復元 C₁₁ = M₁+M₄−M₅+M₇ C₁₂ = M₃+M₅ C₂₁ = M₂+M₄ 素朴: 8乗算 → O(n³)　 Strassen: 7乗算 → O(n^2.807)
図1：Strassen アルゴリズムの再帰構造。$n \times n$ 行列を4ブロックに分割し、7回の $n/2 \times n/2$ 行列乗算で計算する。

5. 実用上の考慮事項

5.1 交差点（閾値）

Strassen は加減算の回数が18回に増えるため、小さい行列ではオーバーヘッドが支配的になる。実際の実装では、再帰の底で行列サイズが閾値 $n_0$ 以下になったら素朴な乗算に切り替える。

一般に $n_0$ は数十〜数百程度であり、ハードウェア（キャッシュサイズ、SIMD幅）やメモリアクセスパターンに依存する。

5.2 数値安定性

Strassen は加減算を多用するため、桁落ちが累積しやすい。浮動小数点演算での誤差解析によれば、素朴な行列乗算の前方誤差は

$$|c_{ij} - \hat{c}_{ij}| \leq n\, u\, \|A\|_\infty \|B\|_\infty + O(u^2)$$

であるのに対し、Strassen は

$$|c_{ij} - \hat{c}_{ij}| \leq O(n^{\log_2 12})\, u\, \|A\|_\infty \|B\|_\infty$$

（$u$ は単位丸め）程度に劣化する。このため、高い精度が求められる場面（逆行列計算、固有値分解の内部など）では注意が必要である。なお、多倍長整数の乗算（丸め誤差なし）では安定性の問題はない。

5.3 奇数サイズの処理

$n$ が奇数の場合、行列を等しく2分割できない。この場合はパディング（1行1列を追加して偶数にする）で対処する。追加された行・列は最後に除去する。

例: $n=3$ の場合、$3 \times 3$ 行列を $4 \times 4$ にゼロ埋めする。

$$A = \begin{pmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{pmatrix} \;\longrightarrow\; A' = \begin{pmatrix} a_{11} & a_{12} & a_{13} & 0 \\ a_{21} & a_{22} & a_{23} & 0 \\ a_{31} & a_{32} & a_{33} & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}$$

Strassen で $C' = A' B'$ を計算後、左上の $3 \times 3$ ブロックを取り出せば $AB$ が得られる。再帰的に呼ぶ場合は、奇数になるたびに各段で同様のパディングを適用する。

5.4 メモリ使用量

再帰の各段で補助積 $M_i$ のための一時領域が必要であり、素朴な乗算に比べてメモリ使用量が増加する。再帰の深さは $O(\log n)$ であるが、各段で $O(n^2)$ の一時領域が必要なため、合計 $O(n^2 \log n)$ の追加メモリを消費する。

5.5 キャッシュ効率

Strassen の再帰的な分割は、素朴な3重ループ (ijk / ikj 順) で達成されるキャッシュ局所性を壊しやすい。特に、補助積 $M_i$ ごとの一時行列の生成・破棄がキャッシュミスを増加させ、メモリ帯域がボトルネックとなる場面が多い。

このため、BLIS や Intel MKL、OpenBLAS などの高性能 BLAS は Strassen を採用せず、キャッシュブロッキング (block × block の階層的ブロック化) と SIMD を組み合わせた素朴乗算で $O(n^3)$ の理論計算量のまま実測性能を引き出している。Strassen は乗算回数では勝つが、キャッシュフレンドリーでないため $n$ が数千程度でも実時間で必ず勝てるわけではない。

5.6 実装上の工夫

  • 作業領域の使い回し: 再帰の各段で新規行列を確保するのではなく、深さ全体で必要なメモリを事前に確保し、補助積間で同じバッファを使い回す。これによりアロケータ呼び出しと GC 負荷を削減する。
  • コピー回避 (ストライド付きビュー): 部分行列をデータコピーで作るのではなく、leading dimension を持つ「ビュー」として表現する。BLAS の lda 引数に相当する設計で、サブブロック演算がインプレースで進む。
  • 底打ちでの素朴乗算 + SIMD: 再帰閾値 $n_0$ (典型的に 64〜256) 以下では Strassen をやめ、SIMD 化された素朴ブロック乗算 (例: $4\times4$ のレジスタブロック) に切り替えると最も速い。
  • 並列化: Strassen の 7 つの補助積は独立に計算できるため、スレッド・タスク並列で 7 並列まで容易にスケールする。

6. Winograd 変種

S. Winograd は Strassen と同じ7回の乗算で、加減算を18回から 15回に削減する変種を示した。漸近的な計算量は同じ $O(n^{2.807})$ であるが、定数項が改善されるため実用的にはやや高速になる。

現在の多くの高性能線形代数ライブラリ（BLAS Level 3 の内部実装など）で使われる Strassen 系アルゴリズムは、実質的に Winograd 変種を採用していることが多い。

7. 行列乗算の指数 $\omega$ の歴史

$n \times n$ 行列の乗算に必要な演算回数を $O(n^\omega)$ と書くとき、$\omega$ を行列乗算の指数と呼ぶ。理論的な下界は $\omega \geq 2$（出力が $n^2$ 個あるため）である。

年 研究者 指数 $\omega$
—素朴な方法$3.000$
1969Strassen$2.807$
1978Pan$2.796$
1981Schönhage$2.548$
1986Strassen$2.479$
1990Coppersmith–Winograd$2.376$
2012Williams$2.3727$
2024Duan–Wu–Zhou$2.371339$

$\omega = 2$ が達成可能かどうかは、計算複雑性理論の主要な未解決問題の一つである。

7.1 Strassen 以降の進展と「銀河系アルゴリズム」

Strassen 以降の改良はテンソル分解と laser method に基づくより精緻な構成を用いる。Coppersmith–Winograd 以降のアルゴリズムは指数を $\omega \approx 2.37$ 台まで押し下げているが、隠れた定数係数と低次項が天文学的に大きいため、現実的な行列サイズ ($n \le 10^{12}$ 程度) では Strassen どころか素朴な乗算にも勝てない。

このような「漸近的には速いが実用域では遅い (あるいは実装不能)」アルゴリズムは 銀河系アルゴリズム (galactic algorithm) と呼ばれる。Strassen はその意味で「銀河系外」の数少ない実用的な高速行列乗算アルゴリズムであり、これが今なお BLAS の高速化や多倍長乗算で生きている理由でもある。

8. まとめ

  • Strassen は $2 \times 2$ 行列の積を 7回の乗算で計算する方法を示した。
  • 再帰的に適用することで、$n \times n$ 行列の乗算を $O(n^{2.807})$ で実行できる。
  • 実用上は、行列サイズが閾値以下なら素朴な乗算に切り替えるハイブリッド戦略が用いられる。
  • 数値安定性はやや劣るが、多倍長整数演算では問題にならない。
  • Winograd 変種は加減算を18回から15回に削減する。
  • 行列乗算の指数 $\omega$ の最小化は、計算複雑性理論の活発な研究分野である。

参考資料

  • V. Strassen, "Gaussian elimination is not optimal," Numerische Mathematik, vol. 13, pp. 354–356, 1969. — DOI: 10.1007/BF02165411
  • T. H. Cormen, C. E. Leiserson, R. L. Rivest, C. Stein, Introduction to Algorithms, 4th ed., MIT Press, 2022. — Chapter 4.2: Strassen's algorithm for matrix multiplication.
  • N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002. — Chapter 23: Fast matrix multiplication.
  • Strassen algorithm — Wikipedia
  • Matrix multiplication algorithm — Wikipedia

関連記事