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})$$

が得られる。加減算は、1回の分割あたり4回から18回へ増える（数えているのはスカラーではなくブロック行列の加減算で、各再帰レベルの分割ごとに18回必要になる）。ただし加減算は $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 は加減算を多用するため、桁落ちが累積しやすい。浮動小数点演算での誤差解析（Higham, Accuracy and Stability of Numerical Algorithms, 第23章。$\infty$ ノルムと、閾値で素朴乗算に切り替える標準的な再帰を仮定）によれば、素朴な行列乗算の前方誤差は

$$|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})$ であるが、定数項が改善されるため実用的にはやや高速になる。

Winograd 変種は、7回の乗算を保ったまま加減算を減らす代表的な改良であり、高速行列乗算を実装する際に検討されることが多い。ただし 5.5 で述べたとおり、主要な高性能 BLAS 実装はそもそも Strassen 系を採用していない。

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$ 台まで押し下げているが、隠れた定数係数と低次項が天文学的に大きいため、現実的な行列サイズでは Strassen どころか素朴な乗算にも勝てない。

このような「漸近的には速いが実用域では遅い (あるいは実装不能)」アルゴリズムは 銀河系アルゴリズム (galactic algorithm) と呼ばれる。Strassen はその意味で「銀河系外」の数少ない実用的な高速行列乗算アルゴリズムであり、実際の実装が試みられてきた数少ない例でもある。

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

関連記事

よくある質問

Q1. Strassen アルゴリズムとは何か

$n \times n$ 行列の乗算を通常の $O(n^3)$ から $O(n^{2.807})$ に高速化する分割統治アルゴリズムである。$2 \times 2$ ブロック行列の積を、8回の乗算の代わりに7回の乗算と18回のブロック加減算で計算する。1969 年に Volker Strassen が発表した。

Q2. 実用上の閾値はどの程度か

閾値は実装・ハードウェア・データ型、そして比較対象のベースケース乗算がどれだけ最適化されているかに強く依存するため、実測で決める必要がある。数十〜数百程度を用いる実装が多い。なお BLIS・Intel MKL・OpenBLAS といった主要な高性能 BLAS は Strassen を採用しておらず、キャッシュブロッキングと SIMD で素朴乗算を高速化している（5.5 参照）。

Q3. Strassen アルゴリズムの欠点は何か

主な欠点は、(1) 浮動小数点演算での誤差の増大、(2) 加減算が増えるため小さい行列ではオーバーヘッドが大きいこと、(3) 素朴な再帰実装ではメモリトラフィックと一時領域が増えること、(4) 奇数サイズの行列にはパディングが必要なこと、である。