mutable の設計

設計目標

mutable はリポジトリの実行指向の半分です。行列を行優先の平坦な 1 本の配列に格納し、呼び出し側がインプレースで更新したり動的なビューを通じて操作したりできるようにし、浮動小数点の数値ルーチンを実装します。行列式、逆行列、階数（ランク）、行簡約、Cholesky 分解、対称行列の固有値、べき乗法です。公開 API は可能な限り値の操作のように読めるようにしてあります。変更は Unit を返すメソッドとビューに限られるので、呼び出し側はシグネチャから、呼び出しがレシーバーを変更するかどうかを判断できます。

数学的背景

以下では、 $u$ は丸め単位（Double では $2^{-53}$ 、Float では $2^{-24}$ ）、 $\mathrm{fl}(x \circ y) = (x \circ y)(1 + \delta)$ （ $|\delta| \le u$ ）、 $\gamma_n = n u / (1 - n u)$ とします。行列間の不等式は要素ごとに成り立つものとし、 $|A|$ は絶対値をとった行列です。

行列積

$(AB)_{ik} = \sum_j a_{ij} b_{jk}$ のコストは $rcn$ 回の積和です。カーネルによって総和の順序は異なります（展開したカーネルでは 1 ステップに 4 つの部分積、 $4 \times 16 \times 16$ 以上の積では $B$ の列をパックしたコピー）が、どの順序でも次を満たします。

\big|\mathrm{fl}(AB) - AB\big| \le \gamma_n |A|\,|B| ,

したがって、異なるターゲット上の結果はこの精度で一致しますが、ビット単位では一致しません。

部分ピボット選択付き LU 分解

正方行列 $A$ について、部分ピボット選択付き Gauss 消去は、置換行列 $P$ 、単位下三角行列 $L$ 、上三角行列 $U$ を計算し、次を満たします。

PA = LU .

ステップ $k$ では、 $|a^{(k)}_{pk}|$ が最大となる行 $p \ge k$ を選んで位置 $k$ に交換し、 $i > k$ について乗数 $l_{ik} = a^{(k)}_{ik} / a^{(k)}_{kk}$ を格納して $a^{(k+1)}_{ij} = a^{(k)}_{ij} - l_{ik} a^{(k)}_{kj}$ と更新します。ピボット選択により $|l_{ik}| \le 1$ が保証されます。コストは $\tfrac23 n^3$ flops です。

行列式。 $PA = LU$ の行列式をとり、 $\det L = 1$ と、交換が $s$ 回のときの $\det P = (-1)^{s}$ を使うと、

\det A = (-1)^{s} \prod_{k} u_{kk} .

求解。 $Ax = b$ は $L y = P b$ （前進代入）と $U x = y$ （後退代入）になり、それぞれ右辺 1 つあたり $n^2$ flops です。逆行列は $I$ の $n$ 本の列に対する解であり、 $\tfrac23 n^3 + n \cdot 2n^2 = \tfrac83 n^3$ flops です。

安定性。 計算された因子は後退誤差の限界（Wilkinson; Higham, Theorem 9.3）を満たします。¹1 N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002, 第 9 章（LU）、第 10 章（Cholesky）、第 8 章（三角行列系）。

\hat L \hat U = P(A + \Delta A), \qquad |\Delta A| \le \gamma_n |\hat L|\,|\hat U| .

$|l_{ik}| \le 1$ からこれは $\lVert \Delta A \rVert_\infty \le n \gamma_n \rho_n \lVert A \rVert_\infty$ を与えます。ここで $\rho_n = \max_{i,j,k} |a^{(k)}_{ij}| / \max_{i,j} |a_{ij}|$ は増大因子です。部分ピボット選択は $\rho_n \le 2^{n-1}$ と抑えます。この限界に達するのは人為的な行列だけで、実際には $\rho_n$ は小さいため、この方法は実用上後退安定です。求解の前進誤差は条件数によって支配されます: $\lVert \hat x - x \rVert / \lVert x \rVert \lesssim \kappa(A)\, n \gamma_n \rho_n$ 。

小さな行列式の閉じた公式

$n \le 4$ では、行列式は LU ではなく公式で評価します。規則 $ad - bc$ 、 $n = 3$ では第 1 行に沿った余因子展開、 $n = 4$ では最初の 2 行に沿った Laplace 展開による、相補的な $2 \times 2$ 小行列式の 6 つの積です。

\det A = \sum_{p<q} (-1)^{p+q+1} \det A_{\{0,1\},\{p,q\}} \det A_{\{2,3\},\overline{\{p,q\}}} .

これらは除算もピボット選択の判断も使わないので、小さな行列では許容誤差の判定を避けられます。LU の意味での後退安定性はありません。悪条件の入力では $ad - bc$ の桁落ちによって相対精度がすべて失われることがありますが、それはまさにそこで行列式自体が悪条件であるのと同じです。

階数（ランク）と被約行階段形

rank はコピーに対して部分ピボット選択付きの消去を行い、絶対値が許容誤差以上のピボットを数えます。これは絶対しきい値 $\tau$ による 数値的階数、つまり $\ge \tau$ のピボットの個数です。ピボットが $\tau$ から十分離れている行列では厳密ですが、そうでなければスケーリングに依存します。（信頼できる数値的階数 $\#\{\sigma_i > \tau\}$ は特異値分解で得られますが、実装されていません。）

reduce_row_elimination は Gauss–Jordan 消去です。各ピボット行をスケーリングしてピボットを 1 にし、ピボット列の上下を消去します。コストはおよそ $r c \min(r, c)$ flops で、インプレースで動作します。

Cholesky 分解

対称正定値（SPD）行列は、 $L$ が下三角で $l_{jj} > 0$ である一意の分解 $A = L L^{\mathsf T}$ を持ちます。 $i \ge j$ について $A = L L^{\mathsf T}$ の要素を比較すると、

a_{ij} = \sum_{k=0}^{j} l_{ik} l_{jk} \quad\Longrightarrow\quad l_{jj} = \sqrt{a_{jj} - \sum_{k<j} l_{jk}^2}, \qquad l_{ij} = \frac{1}{l_{jj}} \Big(a_{ij} - \sum_{k<j} l_{ik} l_{jk}\Big) .

コードはこれを行ごとに（Cholesky–Banachiewicz の順序で）評価し、 $\tfrac13 n^3$ flops かかります。ステップ $j$ の根号の中身は、ピボット選択なしの Gauss 消去のピボットに等しく、それは首座小行列式の比 $\det M_{j+1} / \det M_j$ です。Sylvester の判定法により、 $A$ が SPD であるのはすべての首座小行列式が正であるときに限るので、分解が成功するのはちょうど SPD 行列のときです。is_positive_definite が分解を試みることで実装されているのはこのためです。Cholesky 分解にはピボット選択が不要です。 $a_{jj} = \sum_k l_{jk}^2$ からすべての $|l_{jk}| \le \sqrt{a_{jj}}$ となり要素は増大しえず、計算された因子は $|\Delta A| \le \gamma_{n+1} |\hat L|\,|\hat L^{\mathsf T}|$ である $\hat L \hat L^{\mathsf T} = A + \Delta A$ を満たします（Higham, Theorem 10.3）。

対称固有値問題

実対称行列 $A$ について、スペクトル定理は $Q$ を直交行列、 $\Lambda$ を実対角行列として $A = Q \Lambda Q^{\mathsf T}$ を与えます。eigen はこれを 2 段階で計算します。

Householder 三重対角化。 Householder 鏡映 $H = I - 2 v v^{\mathsf T} / v^{\mathsf T} v$ は対称かつ直交で、ベクトルを座標ベクトルの定数倍に写すことができます。 $n - 2$ 個の鏡映を両側から適用すると、

Q_1^{\mathsf T} A Q_1 = T, \qquad Q_1 = H_1 H_2 \cdots H_{n-2},

となります。ここで $T$ は対角 $d$ 、副対角 $e$ を持つ対称三重対角行列です。コードは $Q_1$ を明示的に累積し、合わせておよそ $\tfrac83 n^3$ flops かかります。鏡映を作る前に各行をその絶対値の和でスケーリングすることで、 $\lVert x \rVert_2$ のオーバーフローとアンダーフローを防ぎます。

Wilkinson シフト付き陰的 QL 法。 三重対角行列 $T$ は平面回転によって対角化されます。 $l$ から始まる未簡約ブロックへの各掃き出しの前に、シフト $\sigma$ として、先頭の $2 \times 2$ ブロック $\begin{pmatrix} d_l & e_l \\ e_l & d_{l+1} \end{pmatrix}$ の固有値のうち $d_l$ に近いほうをとります。 $g = (d_{l+1} - d_l)/(2 e_l)$ とすると、このブロックの固有値は

\lambda_\pm = d_l + e_l \big(g \pm \sqrt{g^2 + 1}\big),

であり、 $d_l$ に近いほうは

\sigma = d_l + e_l\big(g - \operatorname{sgn}(g)\sqrt{g^2+1}\big) = d_l - \frac{e_l}{g + \operatorname{sgn}(g)\sqrt{g^2 + 1}} ,

です。コードが使う 2 番目の形は桁落ちを避けます。掃き出しでは、生じたバルジを追い出す Givens 回転を適用し、同じ回転で $Q$ を更新します。副対角要素は次の場合にゼロとみなします。

|e_m| \le \tau\,\big(|d_m| + |d_{m+1}| + 1\big),

この判定は、大きな対角要素に対しては相対的、ゼロ付近では絶対的です。Wilkinson シフトによる収束は、対称三重対角行列については実際上 3 次であり、有限の入力で失敗した例は観測されていません。それでもコードは、1 つの固有値について 60 回の掃き出し後に中断します。手続き全体は直交変換の積であり後退安定です。計算された固有値は $\lVert \Delta A \rVert_2 = O(u)\lVert A \rVert_2$ である $A + \Delta A$ に対して厳密なので、Weyl の不等式により

|\hat\lambda_i - \lambda_i| \le \lVert \Delta A \rVert_2 = O(u)\,\lVert A \rVert_2 .

固有ベクトルの精度は $u\lVert A \rVert / \text{gap}$ に比例します。ここでギャップとは、最も近い他の固有値までの距離です。

$2 \times 2$ の場合。 $A = \begin{pmatrix} a & b \\ c & d \end{pmatrix}$ について、特性多項式 $\lambda^2 - (a + d)\lambda + (ad - bc)$ から、 $m = (a + d)/2$ として次が得られます。

\lambda_{1,2} = m \pm \sqrt{m^2 - (ad - bc)} .

$b \ne 0$ のとき、ベクトル $(b, \lambda - a)^{\mathsf T}$ は固有ベクトルです。

\begin{pmatrix} a - \lambda & b \\ c & d - \lambda \end{pmatrix} \begin{pmatrix} b \\ \lambda - a \end{pmatrix} = \begin{pmatrix} 0 \\ bc - (\lambda - a)(\lambda - d) \end{pmatrix} = 0 ,

特性方程式により $(\lambda - a)(\lambda - d) = \lambda^2 - (a + d)\lambda + ad = bc$ だからです。コードはこれらのベクトルを正規化せずに返します。 $|\lambda_2| \ll |\lambda_1|$ のとき、減算 $m - \sqrt{\cdot}$ で桁落ちが起こり $\lambda_2$ の相対精度が失われます。安定な代替は $\lambda_2 = \det A / \lambda_1$ です。

べき乗法

$x_0 = (1, \dots, 1)$ （ $A x_0 = 0$ の場合は座標ベクトル）から出発し、この方法は次を反復します。

y = A x_k, \qquad x_{k+1} = y / \lVert y \rVert_\infty, \qquad \lambda_k = \frac{x_k^{\mathsf T} A x_k}{x_k^{\mathsf T} x_k} ,

そして $\lVert A x_k - \lambda_k x_k \rVert_\infty \le \tau$ となったら停止します。 $A$ が固有値 $|\lambda_1| > |\lambda_2| \ge \dots$ で対角化可能であり、 $x_0$ が固有ベクトル $v_1$ 方向に成分 $c_1 \ne 0$ を持つならば、

A^{k} x_0 = \lambda_1^{k}\Big(c_1 v_1 + \sum_{i \ge 2} c_i \big(\tfrac{\lambda_i}{\lambda_1}\big)^{k} v_i\Big),

統計量

variance は 2 パスで計算する母分散です。まず平均 $\bar a$ を求め、次に $\tfrac1N \sum (a_i - \bar a)^2$ を求めます。1 パスの公式 $\tfrac1N \sum a_i^2 - \bar a^2$ は、データの平均が大きく散らばりが小さいときにほぼ等しい 2 つの数を引き算することになり、負の値を返すことさえあります。2 パスの形は非負の項を足し合わせ、その丸め誤差は分散自体に対して相対的です（Chan, Golub and LeVeque, 1983）。個数 $N$ は T において 1 の和として累積されるので、Double では $2^{53}$ 個、Float では $2^{24}$ 個まで厳密です。

転置ビューと積

Transpose::mul は、ラップした行列に対して行列カーネルを再利用し、 $A^{\mathsf T} B^{\mathsf T}$ を $(BA)^{\mathsf T}$ として計算します。要素ごとに見ると、

\big(A^{\mathsf T} B^{\mathsf T}\big)_{ik} = \sum_j a_{ji} b_{kj}, \qquad \big((BA)^{\mathsf T}\big)_{ik} = \sum_j b_{kj} a_{ji} ,

であり、両者はスカラーが可換なときに一致します。Tolerance を持つスカラー型（Float、Double）はすべて可換ですが、Transpose::mul が要求するのは AddMonoid + Mul だけです。非可換なスカラー型では、結果は逆順の積になります（algebra の設計を参照）。

設計上の判断

行優先の平坦なストレージ

選択肢。 行配列の配列、永続的な構造、平坦な 1 本の配列。決定。 要素 $(i, j)$ を $i c + j$ に置く 1 本の Array[T] で、4 つのターゲットすべてで共通。理由。 座標ごとに 1 回の範囲検査で $O(1)$ のアクセスが得られ、消去法や積の内側のループで行が連続し、行と列のビューをコストなしで作れます。from_array は大きな入力をコピーしなくて済むように呼び出し側の配列を採用します。その代償はエイリアシングで、API でドキュメント化しています。

コピーではなくビュー

row_view、col_view、to_transpose は動的なビューを $O(1)$ で返します。ビューは（行列, インデックス）の組かラッパーなので、書き込みは共有ストレージに対して行われ、同期は不要です。実体化は常に明示的です（to_vector、materialize、transpose）。

意味論を共有するターゲット固有のカーネル

このパッケージは、行列・LU・ビュー・転置のコードについてターゲットごとに 1 つのソースファイルを持ちます。違いは各バックエンド向けに選んだループ構造（展開、パッキング、インデックス計算での除算の回避）だけで、範囲やエラーの振る舞いを含む公開された意味論は同一であり、テストは 4 つのターゲットすべてで実行されます。総和の順序が異なるため、浮動小数点の結果は最下位ビットで異なることがあります。

許容誤差に基づく判定

消去法では、計算されたピボットが「ゼロである」のはいつかを判定しなければなりません。このパッケージは Double と Float に対して 1 つの絶対しきい値 $\tau$ = Tolerance::tolerance() = $10^{-11}$ を使い、次のように適用します。

ルーチン	テスト
LU（ $n \ge 5$ の `determinant`、`inverse`、`is_invertible`）	ピボット $\lvert u_{kk}\rvert < \tau$ なら特異
`rank`	残りの最大の $\lvert a_{ik}\rvert < \tau$ ならピボットなし
`reduce_row_elimination`	$\le \tau$ の要素はゼロに設定
`cholesky_decomposition`	根号の中身が $\le \tau$ なら正定値でない
`is_symmetric`、高速経路の判別（単位、対角、置換、三角）	$\lvert a_{ij} - b_{ij}\rvert \le \tau$
`eigen` のデフレーション	$\lvert e_m\rvert \le \tau(\lvert d_m\rvert + \lvert d_{m+1}\rvert + 1)$
`power_method`	残差 $\lVert Ax - \lambda x\rVert_\infty \le \tau$

絶対しきい値は単純で予測しやすいものの、スケール不変ではありません。 $A$ を $10^{-12}$ 倍するとすべてのピボットが $\tau$ を下回るため、完全に良条件の行列でも特異と報告されます。 $10^{12}$ 倍すると、数値的に特異な行列が通ってしまいます。これらのルーチンを呼ぶ前に、データを 1 程度の大きさにスケーリングしてください。Float では $\tau = 10^{-11}$ は丸め単位 $\approx 6 \times 10^{-8}$ よりはるかに小さいので、判定は実質的に厳密なゼロの検査になり、ほぼ特異な Float 行列は検出されません。この trait は閉じている（pub(open) ではなく pub）ので、インスタンスはこの 2 つだけです。スケールを考慮した許容誤差の方針は今後の課題で、導入すれば破壊的変更になります。

高速経路

inverse は単位行列（コピーを返す）、対角行列（対角要素を反転する）、置換行列（転置を返す。列が互いに異なる座標ベクトルである行列では $P^{\mathsf T} P = I$ となるため）を判別します。determinant は $n \ge 5$ で三角行列を判別し、対角要素を掛け合わせます。各判定のコストは $O(n^2)$ で、 $O(n^3)$ の分解を節約します。判定には許容誤差を使うので、 $\tau$ の範囲で対角な行列は厳密に対角として扱われます。

カーネルを変えずに検査付き形式を提供する

すべての検査付きメソッドは検証（正方性、指数の符号、空でないこと、長さ）を行ってから、元の中断する振る舞いまたは Option の振る舞いを保つ検査なしの対を呼び出します。ここには検査付きの行列積はありません。* は検証して中断し、unchecked_matmul はまったく検証しません。@immut.Matrix::matmul とのこの非対称性は認識されています。検査付きの matmul は変更ではなく追加になるでしょう。

対称行列の固有値のみ

一般の実行列は複素固有値を持ちうるため、実数の T に対して Vector[T] を返す関数では表現できません。そのため eigen は、すべての固有値が実数で正規直交な固有基底が存在する対称行列だけを受け付け、それ以外では中断します。

正しさと不変条件

ストレージ。 常に data.length() == row * col です。すべての公開アクセサは行と列を別々に検査します。
値を返すメソッドは変更しない。 行列を変更するのは Unit を返すメソッド、ビューへの書き込み、reduce_row_elimination だけです。
検査付き・検査なしの法則。 前提条件が成り立つときは常に x.f() == Ok(x.unchecked_f()) です。inverse が Err(SingularMatrix) を返すのは、unchecked_inverse が None を返すときに限ります。
行列式の一貫性。 $n \ge 5$ で、 $\tau$ の範囲で三角でない行列について、determinant が 0 を返すのは LU 分解が $\tau$ を下回るピボットを報告するときに限り、それはまた is_invertible が false を返し inverse が失敗するときでもあります。 $n \le 4$ では determinant に閉じた公式を使う一方、is_invertible は引き続き LU を使うので、ごく小さな非ゼロの行列式を持つ行列は、determinant が非ゼロでありながら「逆行列を持たない」と判定されることがあります。
残差の保証。 cholesky_decomposition と eigen は再構成誤差が $O(u)\lVert A \rVert$ の因子を返します。power_method は残差が高々 $\tau$ の組だけを返します。
計算量。 *: $rcn$ 。determinant、inverse: $O(n^3)$ 。rank、reduce_row_elimination: $O(rc\min(r, c))$ 。cholesky_decomposition: $n^3/3$ 。eigen: $O(n^3)$ 。power_method: 反復ごとに $O(n^2)$ 。統計量: $O(rc)$ 。

却下した代替案

相対的な、またはノルムでスケーリングした許容誤差。 より頑健ですが、既存の呼び出し側の結果が変わってしまいます。許容誤差の方針を明示的に渡せるようになるまで見送ります。
複素固有値を返す。 このパッケージが複素数型に依存することになり、実対称な入力に対するシグネチャも変わります。
コピーオンライトのビュー。 書き込みのコストを隠し、ビューの存在意義であるインプレースの契約を壊します。
単一の移植可能なカーネル。 計測の結果、一部のターゲットで遅くなりました。ターゲットごとのファイルは、1 つの仕様を共有しながらコードサイズと引き換えに速度を得ています。

境界

mutable は、右辺に対する連立一次方程式の求解を公開メソッドとして提供せず（逆行列のみ）、QR、SVD、最小二乗ソルバー、非対称行列の固有値、疎なストレージ、条件数の推定、スケールを考慮した許容誤差も提供しません。algebra の trait も自身では実装しません。そのために backends/default がこれをラップしています。回帰や最適化のような分野固有のワークフローは下流のパッケージに属します。

N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002, 第 9 章（LU）、第 10 章（Cholesky）、第 8 章（三角行列系）。 ↩