immut の設計

設計目標

immut は値の意味論を持つ線形代数を提供します。行列やベクトルは一度構築されると決して変わらず、すべての演算は新しい値を返します。そのためプログラムは古いバージョンを保持でき、値をコンポーネント間で自由に共有でき、置換によってコードを推論できます。このパッケージはまた、スカラー型が許す限り 厳密に 計算することを目指しており、整数行列や多倍長整数行列では浮動小数点の近似ではなく、厳密な行列式とべき乗が得られます。

数学的背景

値の意味論と参照透過性

式が参照透過であるとは、プログラムを変えずにその式をその値で置き換えられることです。不変な行列では、

\texttt{let } B = A.\mathtt{set}(i, j, x) \;\Longrightarrow\; A \text{ is the same value before and after,}

となるので、 $A$ に言及する式は更新の前後どちらでも同じ意味を持ちます。そのため代数的な恒等式をそのままプログラム変換として使えます。厳密なスカラーでは $(A + B) + C$ と $A + (B + C)$ は同じ値を表し、中間結果は自由に再利用することも再計算することもできます。

永続的なストレージ

要素は行優先順で moonbitlang/core/immut/vector に格納されます。これは分岐因子 32 のトライとして実装された永続ベクトルです。1 つの要素を置き換えると、根から葉までのパスがコピーされ、それ以外のノードはすべて共有されるので、

\text{cost}(\mathtt{set}) = O(\log_{32} N), \qquad \text{extra memory} = O(\log_{32} N), \qquad N = rc ,

となり、古い行列は有効なままです。読み取りのコストも $O(\log_{32} N)$ で、高々 7 段で $2^{32}$ 個の要素を扱えます。行列全体に対する演算はトライを $O(N)$ で再構築します。

二乗の繰り返しによる行列のべき乗

半環上の正方行列について、 $k = \sum_t b_t 2^t$ と 2 進表記します。すると

A^{k} = \prod_{t : b_t = 1} A^{2^t}, \qquad A^{2^{t+1}} = \big(A^{2^t}\big)^2 ,

となり、 $k - 1$ 回の積の代わりに、高々 $\lfloor \log_2 k \rfloor$ 回の二乗と同じ回数の追加の積で済みます。積の並べ替えが正当なのは、行列の乗算が任意の半環上で結合的だからです（algebra の設計を参照）。因子はすべて同じ $A$ のべき乗なので、スカラーの可換性は必要ありません。pow は状態 $S$ 、指数 $e$ 、底 $B$ を保持し、次の不変条件を保ちます。

S \cdot B^{e} = A^{k},

初期値は $(I, k, A)$ です。各ステップでは、 $e$ が奇数なら $S$ に $B$ を掛け、続いて $e$ を半分にし、 $B$ を二乗します。不変条件は保たれ、 $e = 0$ になったとき状態は $A^k$ です。固定幅整数では、ラップアラウンドする算術がまさにその剰余環の環演算 そのもの なので、オーバーフローしても結果は $\mathbb{Z}/2^{32}\mathbb{Z}$ において厳密です。

分数を使わない行列式

体上の Gauss 消去は $\det A$ をピボットの積として計算しますが、各ステップで除算を行うため整数の範囲から外れます。Bareiss のアルゴリズムは中間値をすべて整数に保ちます。¹1 E. H. Bareiss, “Sylvester’s identity and multistep integer-preserving Gaussian elimination”, Mathematics of Computation 22 (1968), 565–578. $a^{(-1)}_{-1,-1} = 1$ , $a^{(0)}_{ij} = a_{ij}$ とし、 $k = 0, 1, \dots, n-2$ および $i, j > k$ について

a^{(k+1)}_{ij} = \frac{a^{(k)}_{kk}\, a^{(k)}_{ij} - a^{(k)}_{ik}\, a^{(k)}_{kj}}{a^{(k-1)}_{k-1,k-1}} .

Sylvester の行列式恒等式により、各 $a^{(k)}_{ij}$ は、行 $0, \dots, k-1, i$ と列 $0, \dots, k-1, j$ からなる $A$ の小行列式に等しくなります。

a^{(k)}_{ij} = \det \begin{pmatrix} a_{00} & \cdots & a_{0,k-1} & a_{0j} \\ \vdots & & \vdots & \vdots \\ a_{k-1,0} & \cdots & a_{k-1,k-1} & a_{k-1,j} \\ a_{i0} & \cdots & a_{i,k-1} & a_{ij} \end{pmatrix}.

ここから 2 つの帰結が得られます。第一に、漸化式の除算は 割り切れます。分子は直前のピボットの倍数なので、 $\mathbb{Z}$ のような整域上では商は再びその整域に属します。第二に、最後の値は行列式全体 $a^{(n-1)}_{n-1,n-1} = \det A$ です。恒等式の証明は下の添付資料にあります。

Bareiss 消去法: 厳密性と正しさ

ピボット選択を入れても厳密性は損なわれません。 $a^{(k)}$ の行 $i$ は $A$ の行 $i$ と行 $0, \dots, k-1$ にしか依存しないので、ステップ $k$ の前にインデックスが $\ge k$ の 2 つの行を交換することは、それらの行を交換した $A$ に対してアルゴリズムを実行することと同じであり、最終結果に $-1$ を掛けることになります。ピボットより下の列全体がゼロなら、それらの小行列式は消え、行は線形従属であり、 $\det A = 0$ です。

中間値はすべて $A$ の小行列式なので、Hadamard の不等式によってすべて抑えられます。

\big|a^{(k)}_{ij}\big| \le \prod_{r} \lVert \text{row}_r \rVert_2

ここで積は、関係する $k + 1$ 行にわたってとります。漸化式の分子はそうした小行列式の 2 つの積の差なので、その積の 2 乗の 2 倍で抑えられます。これにより Int のオーバーフロー判定基準が得られます。 $H = \prod_r \max(1, \lVert \text{row}_r \rVert_2)$ が $2H^2 < 2^{31}$ を満たせば、どの中間値もオーバーフローせず、結果は厳密です（Int64 では $2H^2 < 2^{63}$ ）。BigInt ではアルゴリズムは常に厳密で、中間の整数はすべて $2H^2$ で抑えられ、算術演算は $O(n^3)$ 回です。

$n \le 4$ では、このパッケージは代わりに閉じた公式を使います。 $n = 3$ では第 1 行に沿った余因子展開、 $n = 4$ では最初の 2 行に沿った Laplace 展開

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

で、互いに相補的な $2 \times 2$ 小行列式の 6 つの積からなります。これらの公式は環演算しか使わないので、どの環でも厳密であり、除算をまったく必要としません。

遅延行列

MatrixFn は $f : [r] \times [c] \to T$ を伴う組 $(\text{shape}, f)$ です。演算は関数を合成します。 $\mathtt{map}(g)$ は $g \circ f$ 、転置は $f \circ \mathrm{swap}$ 、積は

(f \cdot g)(i, k) = \sum_{j} f(i, j)\, g(j, k) ,

で、必要に応じて評価されます。何もキャッシュされないので、積の 1 要素のコストは、内側の次元に因子の要素のコストを掛けたものです。したがって $n \times n$ 行列に対する深さ $d$ の積の木は 1 要素あたり $O(n^{d})$ のコストがかかります。遅延したべき乗は構築するのは安価ですが、読み取るのは高価です。

設計上の判断

永続ベクトルによるストレージ

選択肢。 (a) 更新ごとにコピーする Array、(b) 永続トライ、(c) 関数のみ。決定。 Matrix と Vector には (b)、(c) は MatrixFn として別に提供。理由。 配列をコピーすると set のコストは $O(N)$ になりますが、トライなら読み取りを高速に保ち値を不変にしたまま $O(\log_{32} N)$ にできます。関数は構造的な行列や記号的な行列に有用ですが、読み取りのコストが演算の履歴に依存するため、コストモデルを隠さず明示した別の型としています。

スカラーが許すところでは厳密なアルゴリズム

determinant は体ではなく Compare + Num + Div だけを要求するので、Int、Int64、BigInt を受け付けます。Bareiss 消去法により、これらの型の上では結果は厳密です。Double 上ではスケーリングされたピボットによる消去法のように振る舞います。浮動小数点を対象とする mutable パッケージは、代わりに部分ピボット選択付き LU と許容誤差を使います。

短い名前は検査付き、検査なしは明示的な名前

matmul、trace、determinant、pow は Result を返し、対となる unchecked_* は中断します。演算子 +、-、* は Result を返せないので、形状が一致しない場合は中断します。検査付き形式は検証してから検査なし形式を呼ぶので、定義域上では構成上 checked(x) == Ok(unchecked(x)) という法則が成り立ちます。

`mutable` との整合

両パッケージが同じ操作を提供する場合、名前、引数の順序、検査付き・検査なしの慣例は @mutable と一致しており、consistency のテストがその結果を比較します。違いは意図的なもので、ここに列挙します。

操作	`immut`	`mutable`
単位行列	`Matrix::identity(n)`	トップレベルの `identity(n)`
更新	`set` は新しい行列を返す	`set` はインプレースで書き込む
`m[r][c] = x`	なし	あり
行列式	Bareiss、整域上で厳密	許容誤差付き LU、浮動小数点
分解、逆行列、統計量	なし	あり
`dot`	`Vector` にはない（`ImmutableDenseVector::dot` を参照）	`Vector::dot`

`Vector` には減算がない

Vector は Add、Mul、Neg を実装しますが Sub は実装しません。u - v は u + -v と書きます。これは以前のリリースから引き継いだ非対称性であり、数学的な主張ではありません。backends/default のラッパーは - を提供します。

正しさと不変条件

不変性。 既存の Matrix、Vector、MatrixFn を変更する公開操作はありません。set と swap_* は新しい値を返します。
形状の不変条件。 $r, c \ge 0$ であり、内部のベクトルはちょうど $rc$ 個の要素を持ちます。コンストラクタはそれ以外を中断によって拒否します。
範囲。 すべての公開読み取りは行と列を別々に検査するので、次の行の中に入ってしまう列インデックスは読み取られずに拒否されます。
厳密性。 BigInt 上では determinant と pow は厳密です。Int と Int64 上では、pow は $2^{w}$ を法として厳密で、determinant は小行列式がオーバーフローしない限り厳密です（Hadamard の限界を参照）。
空のケース。 $0 \times 0$ 行列の $\det$ と $\operatorname{tr}$ は $1$ と $0$ です。 $A^0 = I$ です。内側の次元が $0$ の積は零行列です。
計算量。 set: $O(\log_{32} N)$ 。map、+、transpose、swap_*: $O(N)$ 。matmul: $O(rcn)$ 。determinant: $O(n^3)$ 。pow: $O(n^3 \log k)$ 。

却下した代替案

immut の行列式に浮動小数点 LU を使う。 体と許容誤差が必要になり、このパッケージを使う主な理由である整数行列での厳密性が失われます。
Matrix の内部での実行時バックエンドセレクタ。 バックエンドは別々の型（backends/default）であり、Matrix の表現は 1 つです。
MatrixFn でのキャッシュ。 純粋な値が可変な状態を持つことになり、コストモデルが黙って変わります。代わりに Matrix::make で実体化してください。

境界

immut はインプレース更新、ビュー、逆行列、分解、固有値、統計量、許容誤差に基づく述語を提供しません。それらは mutable の担当です。整数オーバーフローを検出せず、検査付きの MatrixFn API も提供せず、algebra の trait も自身では実装しません（backends/default のラッパーが実装します）。

E. H. Bareiss, “Sylvester’s identity and multistep integer-preserving Gaussian elimination”, Mathematics of Computation 22 (1968), 565–578. ↩