Matrix 算術アルゴリズム

概要

Matrix<T> は sangi の動的サイズ行列であり、 StaticMatrix<T,M,N> はコンパイル時固定サイズの変種である。 ヘッダオンリーで、ヘッダは math/core/matrix.hpp。 本記事では算術演算 (加減乗、転置、要素積、Kronecker 積、トレース) の内部表現と漸近挙動を扱う。 分解・逆行列・行列式は 行列式と逆行列 で別途扱う。

API の使い方については Matrix API リファレンス を参照。

内部表現と行優先・列優先

row-major レイアウト

sangi の Matrix<T> は行優先 (row-major) で要素を格納する。 $M \times N$ 行列 $A$ の要素 $a_{ij}$ は線形バッファ上で位置 $i \cdot N + j$ に置かれる。

$$\mathrm{offset}(i, j) = i \cdot N + j \quad (0 \leq i < M, 0 \leq j < N)$$

これは C/C++ の 2 次元配列、Python NumPy のデフォルト、PyTorch / TensorFlow のテンソルと整合する。 一方で Fortran / LAPACK / BLAS は列優先 (column-major) である。 この差は転置 + leading dimension の指定で吸収できる:

$$A_{\mathrm{rowmajor}} = (A^\top)_{\mathrm{colmajor}}$$

すなわち row-major で格納された $A$ は、col-major で見ると $A^\top$ と等価である。 BLAS の dgemm を呼び出す際は CblasRowMajor オプションを渡すか、 $C = AB$ を $C^\top = B^\top A^\top$ に書き換えて転置 leading dimension を渡せばよい。

静的版 StaticMatrix<T,M,N>

StaticMatrix<T,M,N> は std::array<T, M*N> をストレージとする。 ヒープ確保がなく、$M, N$ がコンパイル時定数なので全ループが完全アンロール対象となる。 $4 \times 4$ 行列 (SE(3) 変換) や $3 \times 3$ (回転) で Matrix<T> より大幅に速い。

leading dimension とビュー

MatrixMap<T> 等のビュークラスは leading dimension (LD) を持つ。 $M \times N$ 行列の部分行列 (たとえば左上 $r \times c$ ブロック) を抜き出すとき、 部分行列の行と行のあいだに余白 (パディング) が生じる:

$$\mathrm{offset}(i, j) = i \cdot LD + j, \quad LD = N \geq c$$

このため部分行列のメモリ走査は連続ではなく、SIMD パケットロードは行内のみで有効となる。

加減算と転置

行列の加減算 $C = A \pm B$ は $M \cdot N$ 個の要素ごとの演算であり、 Vector の加減算と同じ機構 (式テンプレート + SIMD パケット評価) で実装される。 計算量は $O(MN)$、メモリ帯域律速で 4 アキュムレータ展開の恩恵を直接受ける。

転置

転置 $B = A^\top$ は数学的には要素の入れ替えだが、メモリ上では行と列が交差するため、 最も素朴な実装ではキャッシュミスを大量に発生させる。 sangi はキャッシュタイル転置を採用し、 L1 キャッシュに収まる $b \times b$ ブロック単位で in-place 転置を行う:

$$B[i][j] = A[j][i], \quad (i, j) \in [bi, bi+b) \times [bj, bj+b)$$

ブロックサイズ $b$ は L1 キャッシュサイズに基づいて選ぶ (L1 = 32 KiB、$T = $ double (8 B)、2 行列分で $b = 32$ あたりが目安)。 $b \times b$ 転置内では SIMD shuffle 命令 (_MM_TRANSPOSE4_PS 系) が使える。

転置は明示呼び出しでバッファに書き出すほか、 Lazy 評価で「転置ビュー」として扱う設計も可能だが、sangi は素直なエージャ評価を採用している。 後段で行列乗算等を行う際、トランスポーズ済みデータを用意する方が一貫した高速化が得られるためである。

関連記事: 線形空間の基礎

行列乗算 (一般 N×N)

$C = AB$、$A: M \times K$、$B: K \times N$、$C: M \times N$ の場合:

$$c_{ij} = \sum_{k=0}^{K-1} a_{ik} b_{kj}$$

素朴実装は 3 重ループで $O(MNK)$。行優先 $A$、行優先 $B$ の組み合わせでは $b_{kj}$ のアクセスがストライド $N$ になりキャッシュ局所性が悪い。 sangi の対応は次の三層構造である。

サイズ閾値選択アルゴリズム計算量
小 ($N \lesssim 64$)素朴 ijk ループ + SIMD$O(N^3)$
中 ($64 \lesssim N \lesssim 512$)ブロック乗算 + 内側 micro-kernel$O(N^3)$
大 ($N \gtrsim 1024$)Strassen 法 + 再帰下端でブロック$O(N^{\log_2 7}) \approx O(N^{2.807})$

BLAS バックエンド (OpenBLAS / Intel MKL) が利用可能な場合、中・大サイズで dgemm を呼ぶ経路もある。 BLAS なしでも素直なブロック実装でメモリ階層の恩恵をほぼ引き出せる構成になっている。

ブロック化と CPU キャッシュ

行列乗算のキャッシュ最適化はメモリ階層 (L1, L2, L3) に応じて 3 段のブロックサイズ $m_c, k_c, n_c$ を選ぶ Goto-Van de Geijn 流の枠組みが標準である。 sangi は最内 micro-kernel をハンドチューンせず、コンパイラの SIMD 自動ベクタライズに委ねる方針を取り、 その代わり外側のブロックサイズを実機ベンチで決定する。

$$C_{ij} \mathrel{+}= A_{ip} \cdot B_{pj}$$

($A_{ip}$ は $m_c \times k_c$、$B_{pj}$ は $k_c \times n_c$、$C_{ij}$ は $m_c \times n_c$)

  • $k_c$: $A$ のパネルが L1 キャッシュに収まる範囲。double なら $k_c \approx 256$
  • $m_c$: $A$ のパネル + $B$ のスライスが L2 に収まる。$m_c \approx 96$ 程度
  • $n_c$: $B$ のパネルが L3 に収まる。$n_c \approx 4096$

$B$ をブロック単位でパッキングすることで内側ループのアクセスを連続化できる。 パッキングコストは 1 回限りなので、$k_c, n_c$ が十分大きければ償却される。

BLAS バックエンドを使う場合はこれらの定数調整はライブラリ側に委ねる。 sangi 内蔵実装は移植性重視で $m_c = 96, k_c = 256, n_c = 4096$ をデフォルトとし、 ターゲットマシンで自動チューニングする仕組みを将来追加予定。

Strassen 法

Strassen 法 (1969) は行列乗算を再帰的に分割し、 $2 \times 2$ ブロック分割における 8 個の部分積を 7 個に削減する。 通常の $O(N^3)$ より漸近的に高速で、計算量は $O(N^{\log_2 7}) \approx O(N^{2.807})$。

$A, B, C$ を $2 \times 2$ のブロックに分割する:

$$A = \begin{pmatrix} A_{11} & A_{12} \\ A_{21} & A_{22} \end{pmatrix}, \quad B = \begin{pmatrix} B_{11} & B_{12} \\ B_{21} & B_{22} \end{pmatrix}$$

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

結果ブロックを線形結合で復元する:

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

閾値と数値安定性

Strassen は加減算の回数が増えるため、小サイズではブロック古典法より遅い。 sangi の閾値は$N \approx 1024$ ($T = $ double の場合) であり、 再帰下端では古典ブロック法に切り替わる。

数値安定性については、Strassen の誤差境界は古典法より緩いが、 実用上は IEEE 754 倍精度で $N \leq 10^4$ 程度なら検知できる差異は出ない。 sangi のデフォルトでは Strassen は有効化フラグ付きであり、明示的に要求された場合のみ走る。 多くの用途では古典ブロック法 + BLAS バックエンドの方が定数倍で勝つためである。

関連記事: Strassen 法 (数値解析・上級)

Hadamard 積

Hadamard 積 (要素ごとの積) は記号 $\circ$ で表され:

$$(A \circ B)_{ij} = a_{ij} \cdot b_{ij}$$

計算量は $O(MN)$、Matrix 加減算と同じ機構 (式テンプレート + SIMD) で実装される。 Neural Network のゲーティング (LSTM/GRU の forget gate 等) や 確率モデルの要素ごと積でよく使われる。

$A \circ B$ は行列乗算 $AB$ とは別演算であり混同に注意。 Matrix 型では現状フリー関数 hadamard(A, B) は Tensor 型に対してのみ提供しており、 Matrix の要素ごと積はループまたは式テンプレート経由で記述する。Matrix オーバーロードの追加は今後の課題である。

関連記事: Hadamard 積

Kronecker 積

Kronecker 積 $A \otimes B$ は $A$ が $M \times N$、$B$ が $P \times Q$ のとき $MP \times NQ$ の行列を生成する。

$$(A \otimes B)_{(i p + k)(j q + l)} = a_{ij} \cdot b_{kl}$$

計算量は $O(MNPQ)$、出力サイズは入力の積に膨れる。 テンソル積空間の表現、Sylvester 方程式 $AX + XB = C$ の vec 化、 群表現論などで頻出する。

sangi の実装では、出力行列を 4 重ループで埋める。 外側 2 重ループ ($i, j$ over $A$) ごとに B のブロックを スケール乗算しながらコピーする形になり、 最内ループ ($k, l$ over $B$) は連続メモリ書き込みで SIMD パケット化される。

応用上 $A \otimes B$ を陽に作らず、$\mathrm{vec}(BXA^\top) = (A \otimes B) \mathrm{vec}(X)$ の関係で 作用だけ評価する方が効率的なケースが多い。 sangi の Sylvester ソルバはこの陰的な形を使う。

関連記事: Kronecker 積

トレース

$$\mathrm{tr}(A) = \sum_{i=0}^{N-1} a_{ii}$$

正方行列の対角和。計算量は $O(N)$。 対角要素のアドレッシングは $i \cdot N + i$ で stride $N + 1$ のアクセスになり、SIMD パケットロードは効かない。 $N$ が大きい場合でも対角要素は連続でないため、純粋スカラループで十分。

$\mathrm{tr}(AB) = \sum_{i,j} a_{ij} b_{ji}$ の関係は frequenctly 利用される。 sangi の trace(A * B) は式テンプレートを通じて中間 $AB$ を作らずに済むが、 完全な融合は実装側で trace_product(A, B) 専用関数を提供する設計である。

不均衡サイズの乗算分割

$A: M \times K$ と $B: K \times N$ で $M, K, N$ が大きく異なる場合、 メモリアクセスパターンが極端に偏り、単一のブロックサイズでは最適化が効きにくい。 典型例:

  • 外積系 ($K \ll M, N$): $K$ が小さく $C$ が低 rank。各列ベクトルの和として直接書ける
  • 内積系 ($M, N \ll K$): $C$ が小さく $K$ 軸の縮約が支配的。各要素は長い dot 積
  • パネル系 ($M \approx K, N \ll K$): 縦長 panel との積。BLAS Level 2.5 相当

sangi は前段で次元比 $\rho = K / \min(M, N)$ を見て分岐する:

  • $\rho < 0.25$ → 外積系として処理 ($A$ の各列 × $B$ の各行 の rank-1 update 集合)
  • $\rho > 4$ → 内積系として処理 ($k_c$ ブロックの長い縮約を最適化)
  • $0.25 \leq \rho \leq 4$ → 通常のブロック乗算

結果として、極端な比率でも常識的なスループットが得られる。 BLAS では dgemm の transA/transB と leading dimension で同等のことを表現できるが、 sangi の内蔵実装ではこのディスパッチを陽に書いている。

設計判断

なぜ row-major なのか

C/C++ の慣習、Python NumPy/PyTorch のデフォルトと一致するため、 外部データとのやり取り (画像、テンソル、Pandas DataFrame 等) でコピーや転置が不要になる。 LAPACK は col-major だが、BLAS Level 3 の各ルーチンは CblasRowMajor オプションでフラットに row-major を扱える。 col-major が必須の場合は転置ビューを噛ませれば良いだけで、コストは漸近的にゼロ。

Strassen をデフォルト無効にする理由

BLAS バックエンドが利用可能な場合、$N \leq 10^4$ では BLAS の古典 GEMM (高度にチューンされた micro-kernel + 並列化) の方が Strassen より速いことが多い。 Strassen は再帰段ごとに $O(N^2)$ の加減算 + アロケーションが乗るため、 定数倍が大きいのである。 sangi は移植性のため BLAS なしでも動くが、Strassen は明示要求時のみ有効にして、デフォルトでは古典ブロック法 + 利用可能なら BLAS を使う。

並列化

行列乗算の外側ループ ($i, j$ blocks) は OpenMP で並列化可能だが、 sangi のデフォルトでは並列化を有効にしていない。 プロセス全体での並列性 (Python multiprocessing、TBB、外部スレッドプール等) と競合しないことを優先する判断である。 ユーザは num_threads パラメータで明示的に並列化を要求できる。

参考文献

  • Strassen, V. (1969). "Gaussian Elimination is not Optimal". Numerische Mathematik, 13, 354–356.
  • Goto, K. and van de Geijn, R. (2008). "Anatomy of High-Performance Matrix Multiplication". ACM Transactions on Mathematical Software, 34 (3), Article 12, 12:1–12:25.
  • Higham, N. J. Accuracy and Stability of Numerical Algorithms. 2nd ed. SIAM, 2002. (Strassen の数値誤差解析)