Vector 算術アルゴリズム

概要

sangi の Vector<T> および StaticVector<T,N> は、 要素型 $T$ の数値ベクトルを式テンプレートと SIMD パケット評価で扱うクラスである。 本記事ではユーザ視点ではなく、内部表現・誤差管理・SIMD 化・$O(\cdot)$ 上の漸近挙動など、 実装側の判断を扱う。

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

内部表現

動的版 Vector<T>

Vector<T> は std::vector<T> をストレージとする 1 次元連続配列である。 要素はメモリ上で隣接して並ぶため、ハードウェアプリフェッチャと SIMD ロード/ストアの恩恵を直接受ける。

  • 要素間ストライド: 常に $\mathrm{sizeof}(T)$ バイト (詰めて配置)
  • アラインメント: 標準アロケータ依存。SIMD パケットロードは未アラインド命令 (_mm_loadu_* 系) で読み出し、アラインメントを要求しない
  • サイズ管理: size() と capacity() を分離し、reserve で前置確保が可能

静的版 StaticVector<T,N>

StaticVector<T,N> は std::array<T,N> をストレージとする。 ヒープ確保が一切発生しないため、3 次元座標やクォータニオン等の小ベクトルで Vector<T> より定数倍速い。 サイズ $N$ がコンパイル時に確定するため、内部ループも完全アンロールされやすい。

両者は CRTP 基底 VecExpr<Derived> を継承し、式テンプレート上では同等に振る舞う。

要素ごとの算術

加算・減算・スカラー乗除算はすべて要素ごとの $O(n)$ 演算である。 式テンプレートが導入されていない場合、c = a + b は一時オブジェクト a + b を作って c に代入する 2 パスとなるが、 sangi では式テンプレートにより単一ループ単一パスで評価される。

$$c_i = \alpha a_i + \beta b_i - \gamma d_i \quad (i = 0, \ldots, n-1)$$

このような複合式 c = alpha*a + beta*b - gamma*d は中間バッファを作らず、 SIMD パケットを 1 個ずつ a, b, d から読み込んで合成し、 結果パケットを c に書き出す。 メモリ帯域律速の問題で 4-5 倍の高速化が得られる。

BLAS Level-1 互換の融合関数 axpy ($\mathbf{y} \mathrel{+}= \alpha \mathbf{x}$) と axpby ($\mathbf{y} = \alpha \mathbf{x} + \beta \mathbf{y}$) も、 単一ループで実装され、コンパイラの自動ベクタライズに好適な形になっている。

内積 (dot)

内積 $\langle \mathbf{a}, \mathbf{b} \rangle = \sum_{i=0}^{n-1} a_i b_i$ は最も頻出する縮約演算である。 単純実装ではループキャリー依存により浮動小数点パイプラインが詰まるため、 sangi は複数アキュムレータ展開を採用する。

$$s = \sum_{k=0}^{3} \left( \sum_{i \equiv k \pmod 4} a_i b_i \right)$$

4 つの独立した SIMD アキュムレータに分けて累積し、最後にそれらを横方向に縮約する。 Skylake 以降の x86 では FMA レイテンシ 4-5 サイクル・スループット 0.5 サイクルなので、 最低 4 並列のアキュムレータが必要であり、sangi の実装はこの理論限界に届く。

Kahan 補正の採否

sangi の dot は素朴な累積 (Kahan 補正なし) を採用している。理由は以下のとおり。

  • BLAS の ddot も Kahan 補正なしであり、業界標準と整合する
  • $\|\mathbf{a}\| \cdot \|\mathbf{b}\|$ でキャンセレーションが起きる場合は、利用側でアルゴリズム的に分割すべき問題
  • Kahan 補正は素朴版の 3-4 倍遅く、汎用 dot として常時負担するコストが大きい

誤差解析: 素朴な累積では相対誤差が $O(n \epsilon)$、4 アキュムレータでは並列累積効果により $O(\sqrt{n} \cdot \epsilon)$ 程度に収まる ($\epsilon$ は機械イプシロン)。 高精度が必要な利用者は、$\mathbf{a}$, $\mathbf{b}$ を Float 多倍長型に持ち上げて dot を呼ぶことで任意精度内積を得られる。

関連記事: 内積の定義と性質 / 桁落ち

外積 (cross)

外積は 3 次元ベクトル空間に固有の演算であり、 sangi では StaticVector<T,3> 専用のフリー関数として提供される。 動的サイズの Vector<T> には cross は存在しない。 サイズ不一致を実行時例外で投げる API は、線形代数の文脈では誤用の温床になりやすいためである。

$$\mathbf{a} \times \mathbf{b} = \begin{pmatrix} a_1 b_2 - a_2 b_1 \\ a_2 b_0 - a_0 b_2 \\ a_0 b_1 - a_1 b_0 \end{pmatrix}$$

SIMD 化のメリットは小さい (3 成分しかないため) が、シャッフル命令 1 回 + FMA 2 回で完了する。 実装は完全アンロールされた 3 行のスカラ式である。

7 次元の外積 (Cayley 数からの誘導) や $n$ 次元の Hodge dual は提供しない。 必要な場合は Matrix 経由の楔積を利用する。

関連記事: 外積 (3 次元)

ノルム L1 / L2 / L∞

ノルムは 3 種類とも内部的に縮約演算で実装されており、SIMD パケット評価の対象となる。

関数定義計算量
norm_l1$\|\mathbf{v}\|_1 = \sum_i |v_i|$$O(n)$
norm (= L2)$\|\mathbf{v}\|_2 = \sqrt{\sum_i v_i^2}$$O(n)$
norm_linf$\|\mathbf{v}\|_\infty = \max_i |v_i|$$O(n)$
norm_lp$\|\mathbf{v}\|_p = (\sum_i |v_i|^p)^{1/p}$$O(n)$

L2 ノルムの overflow / underflow 回避

$\sum v_i^2$ を素朴に計算すると、$|v_i|$ が大きい場合に $v_i^2$ が浮動小数点表現の上限を超え、 小さい場合に $v_i^2$ がアンダーフローしてゼロに丸められる。 $\|\mathbf{v}\|$ 自体は表現可能でも、途中計算で破綻する。

この問題は古典的に「hypot 風スケーリング」で回避される:

$$\|\mathbf{v}\| = m \cdot \sqrt{\sum_i (v_i / m)^2}, \quad m = \max_i |v_i|$$

2 パス必要 (まず max、次にスケーリング後の和) で、素朴版の約 2 倍のコストがある。 sangi の norm はデフォルトで素朴版を採用し、 オーバーフロー耐性が必要な場合は別名 safe_norm (将来追加予定) または hypot を利用するよう設計している。 ベクトル長 $n$ が大きく要素が中庸 ($|v_i| \approx 1$) な通常用途では素朴版で十分であり、 二重スケーリングのコストを常時支払うのは回避している。

L∞ ノルムの SIMD 縮約 (T = double / float の場合)

$T$ が double または float のとき、$\max_i |v_i|$ は _mm256_max_pd 系の lane-wise max 命令で SIMD 化される。 4 アキュムレータに分けた後、水平縮約で最終 max を得る。 絶対値は and マスク (符号ビット消去) で計算するため分岐なし。 IEEE 754 のビットレイアウトに依存した最適化であり、 多倍長型 (Int / Float / Rational) や Complex<T> ではこの SIMD パスを通らず、要素ごとに abs() を呼ぶスカラー実装になる (多倍長型は内部表現が可変長であり 256-bit レーンに収まらない)。

関連記事: ノルムの定義と性質 / ベクトルノルム (数値解析) / オーバーフロー・アンダーフロー

正規化と射影

normalize

normalized(v) は $\mathbf{v} / \|\mathbf{v}\|$ を返す。 $\|\mathbf{v}\| = 0$ の場合、例外を投げる代わりにゼロベクトルを返す (数値計算の中間結果としてゼロベクトルが頻出するため、例外伝播が制御フローを複雑にすることを避ける選択)。

In-place 版 v.normalize() は内部で norm を一度だけ計算してから、 ベクトル全体をその逆数倍する。$1/\|\mathbf{v}\|$ を先に計算してから乗算するため、除算は 1 回のみ。

射影 (projection)

ベクトル $\mathbf{a}$ を $\mathbf{b}$ 方向に射影するには:

$$\mathrm{proj}_{\mathbf{b}}(\mathbf{a}) = \frac{\langle \mathbf{a}, \mathbf{b} \rangle}{\langle \mathbf{b}, \mathbf{b} \rangle} \mathbf{b}$$

sangi では専用関数は提供しておらず、(dot(a,b) / dot(b,b)) * b として 式テンプレート経由で記述する。中間 Vector は生成されず、単一ループでスケーリングされる。 Gram-Schmidt 直交化や Householder 反射子の構築では、この式を内部で多用する。

$\langle \mathbf{b}, \mathbf{b} \rangle$ が小さい場合、商が数値的に不安定になる。 正規化済みの単位ベクトル $\hat{\mathbf{b}}$ を保持しておき、$\langle \mathbf{a}, \hat{\mathbf{b}} \rangle \hat{\mathbf{b}}$ で射影する方が安定。 QR 分解の Householder 経路ではこちらの定式化を採用している。

関連記事: 線形変換と射影 / 射影行列

VectorMap と ストライドビュー

VectorMap<T> と ConstVectorMap<T> は、 外部メモリ領域に対するゼロコピービューを提供する。 既存の C 配列、std::vector、Python numpy.ndarray 等の 生バッファを所有権を移さずに sangi の式テンプレート上で扱える。

ストライドの扱い

VectorMap はストライド付きアクセスをサポートする。たとえば、 $M \times N$ の行列の 1 列を抜き出してベクトルとして読みたい場合、 要素間のストライドは $M$ ワードになる。 ストライド付きビューは:

  • SIMD パケットロードは不可 — 要素が連続していないため gather 命令か逐次ロードに退化する
  • 式テンプレートには参加できる — ノードの addressing 関数経由でアクセスし、コンパイラに最適化を委ねる
  • 意味論的にはコピーと等価 — 結果が同じであることが保証されるだけで、性能は連続アクセスより劣る

ストライド 1 のビュー (連続アクセス) は連続バッファと同じ最適化が適用される。 head, tail, segment はストライド 1 の VectorView を返す。

SIMD と式テンプレート

PacketTraits<T>

sangi は型ごとに SIMD パケット幅を PacketTraits<T> で抽象化する。 現在の specialization は AVX2 と SSE2 の 2 系統で、 double は 4 (AVX2) / 2 (SSE2)、float は 8 (AVX2) / 4 (SSE2)、 Float (多倍長) や Int は 1 (パケットなし、スカラ評価) である。 パケット幅 $W$ で、$n$ 要素のベクトルは $\lfloor n/W \rfloor$ 個のパケット処理 + 末尾 $n \bmod W$ 個のスカラ tail で評価される。 AVX-512 への拡張は PacketTraits 層では未実装 (将来課題) だが、 BLAS 風の自由関数 (dot / axpy / scal 等の simd_backend) は別経路で コンパイル時に AVX-512 専用コードを選択する。

式テンプレートと auto の罠

式テンプレートにより auto v = a + b の v の型は Vector<T> ではなく、 VecBinOp<Plus, Vector<T>, Vector<T>> のような中間ノード型になる。 この中間型は元の a, b への参照を保持しているため、 $\mathbf{a}$, $\mathbf{b}$ のスコープを抜けるとダングリング参照となる。

API リファレンスでも警告しているとおり、計算結果を他の関数に渡したり長く保持する場合は、 Vector<T> v = a + b; のように具体型で受けるか、即時評価関数 (eval()) を経由する。

パケット評価の境界条件

SIMD ロード/ストアは未アラインドアクセス命令を用いるため、 利用者がメモリアラインメントを意識する必要はない。 境界ケースは以下の 2 系統で扱う:

  • $n \geq W$: パケットループ + スカラ tail
  • $n < W$: スカラ評価 (パケットセットアップのコストを避ける)

$n$ が極小 (例: $n = 1$) で大量に呼ばれる場合、パケットセットアップのオーバーヘッドが パケット並列性のメリットを上回るため、純粋スカラに切り替える分岐が含まれている。

設計判断

動的版と静的版を統合しない理由

Eigen のように単一テンプレート Matrix<T, Rows, Cols> で動的・静的の両方を扱う設計は 柔軟だが、ヘッダ展開が大きくコンパイル時間が長い。 sangi は読者がオーバーロード解決の結果を追いやすいよう、 動的版 (Vector<T>) と静的版 (StaticVector<T,N>) を別クラスにしている。 両者は CRTP 基底 VecExpr を共有するため、フリー関数 (dot, norm 等) は共通で書ける。

cross は 3 次元限定

$\mathbb{R}^3$ の外積は Cayley 数 (四元数の虚部) に由来する特異な代数構造であり、 7 次元を除いて他の次元には拡張されない。 本来の外積を $n$ 次元化したものは Hodge dual と楔積であるが、 これは Matrix 経由で別途扱う。 Vector<T> に cross(a,b) を提供してサイズ不一致時に例外を投げる API は誤用の温床になるため避けている。

Kahan 補正をデフォルトにしない理由

素朴な dot は BLAS 標準と一致しており、4 アキュムレータ展開時の経験誤差 $O(\sqrt{n}\epsilon)$ は 多くの応用で十分な精度を持つ。 Kahan 補正版を常時走らせると 3-4 倍遅くなる。 究極的に厳密な内積が必要な場合は、要素型を Float (sangi 多倍長浮動小数点) に置き換える方が、 Kahan 補正を SIMD で実装するより筋がよい。

参考文献

  • Higham, N. J. Accuracy and Stability of Numerical Algorithms. 2nd ed. SIAM, 2002. (浮動小数点内積の誤差解析)
  • Goldberg, D. (1991). "What Every Computer Scientist Should Know About Floating-Point Arithmetic". ACM Computing Surveys, 23 (1), 5–48.
  • Veldhuizen, T. (1995). "Expression Templates". C++ Report, 7 (5), 26–31.