algebra 設計

設計目標

algebra は、汎用的な線形代数コードに、ベクトルや行列のオブジェクト全体を扱うための語彙を与えます。この語彙は、各オブジェクトに何ができるかについて正直です。trait はその法則を満たす型にしか実装してはならず、アルゴリズムは自身が使う構造だけを正確に要求できるべきです。つまり、形状、アーベル群、Hadamard 環、転置、行列積のいずれかです。そのため、このパッケージは 1 つの Matrix trait や VectorSpace trait ではなく、包含関係で順序づけられた小さな trait 群を定義します。trait の一覧は API ページにあります。このページでは、それらが表す数学と、そこから導かれる選択を説明します。

数学的背景

加群とベクトル空間

$R$ を単位元をもつ環とします。（左） $R$ 加群とは、アーベル群 $(V, +, 0, -)$ にスカラー作用 $R \times V \to V$ , $(r, v) \mapsto r v$ を合わせたもので、すべての $r, s \in R$ と $u, v \in V$ について次を満たすものです。

\begin{aligned} r(u + v) &= r u + r v, & (r + s) v &= r v + s v, \\ (r s) v &= r (s v), & 1 v &= v . \end{aligned}

ベクトル空間とは体上の加群です。線形代数の大部分はベクトル空間について述べられますが、コード の大部分が必要とするものはずっと少なくて済みます。行列の加算にはアーベル群が、行列の乗算にはスカラーの半環が必要で、除算を必要とするのは消去法だけです。たとえば整数行列は $\mathbb{Z}$ 加群をなしますがベクトル空間ではありません。体を要求すると、除算を一切しないアルゴリズムからも整数行列が締め出されてしまいます。

座標

$R$ 上の長さ $n$ の密ベクトルは $R^n$ の元であり、密な $m \times n$ 行列は $R^{m \times n}$ の元です。要素ごとの加算によってどちらもアーベル群となり、スカラー作用 $r \cdot (a_{ij}) = (r a_{ij})$ によって $R$ 加群となります。 $n$ を固定すると、要素ごとの積

(u \odot v)_i = u_i v_i

によって $R^n$ は環（直積環）にもなり、その単位元は $(1, \dots, 1)$ です。

圏としての行列

行列の乗算は 1 つの集合上の演算ではありません。その形状規則

(m \times n) \cdot (n \times p) = m \times p

は、行列が自然数を対象とする圏 $\mathbf{Mat}_R$ の射であることを示しています。 $m \times n$ 行列は射 $n \to m$ であり、乗算は合成、単位行列 $I_n$ は $n$ 上の恒等射です。積 $AB$ が定義されるのは、射が合成可能なときに限ります。固定した対象 $n$ 上の射 $n \to n$ （正方行列）は環をなし、これがおなじみの $M_n(R)$ です。

合成の結合律は、 $R$ が半環であれば常に成り立ちます。 $A \in R^{m \times n}$ , $B \in R^{n \times p}$ , $C \in R^{p \times q}$ について:

\begin{aligned} \big((AB)C\big)_{il} &= \sum_{k=1}^{p} (AB)_{ik} C_{kl} = \sum_{k=1}^{p} \Big(\sum_{j=1}^{n} A_{ij} B_{jk}\Big) C_{kl} \\ &= \sum_{k=1}^{p} \sum_{j=1}^{n} A_{ij} B_{jk} C_{kl} && \text{right distributivity, associativity of } \cdot \\ &= \sum_{j=1}^{n} \sum_{k=1}^{p} A_{ij} B_{jk} C_{kl} && \text{associativity and commutativity of } + \\ &= \sum_{j=1}^{n} A_{ij} \Big(\sum_{k=1}^{p} B_{jk} C_{kl}\Big) = \big(A(BC)\big)_{il} && \text{left distributivity.} \end{aligned}

乗算の可換性は一度も使っていないので、この法則は任意の半環上の行列について成り立ちます。四元数のような非可換なものも含みます。 $C = I$ として同じ計算をすれば $AI = A$ が示され、行列の加算に対する積の分配律は $R$ における分配律から要素ごとに従います。

転置

転置は $R^{m \times n}$ から $R^{n \times m}$ への写像 $(A^{\mathsf T})_{ij} = A_{ji}$ です。これは全域的であり、次を満たします。

(A^{\mathsf T})^{\mathsf T} = A, \qquad (A + B)^{\mathsf T} = A^{\mathsf T} + B^{\mathsf T} .

積の規則にはさらに条件が必要です。合成可能な $A$ と $B$ について:

\begin{aligned} \big((AB)^{\mathsf T}\big)_{ik} &= (AB)_{ki} = \sum_j A_{kj} B_{ji}, \\ \big(B^{\mathsf T} A^{\mathsf T}\big)_{ik} &= \sum_j (B^{\mathsf T})_{ij} (A^{\mathsf T})_{jk} = \sum_j B_{ji} A_{kj} . \end{aligned}

2 つの和が項ごとに一致するのは、 $A_{kj} B_{ji} = B_{ji} A_{kj}$ 、すなわちスカラーが可換であるときに限ります。したがって $(AB)^{\mathsf T} = B^{\mathsf T} A^{\mathsf T}$ は可換環上の行列については定理ですが、一般には偽です。¹1 積の順序を反転させる対合 $x \mapsto \bar{x}$ 、つまり $\overline{xy} = \bar{y}\,\bar{x}$ を満たすものを備えた環では、共役転置 $A^{*} = \overline{A}^{\mathsf T}$ は可換性なしでも $(AB)^{*} = B^{*} A^{*}$ を満たします。各項が $\overline{A_{kj} B_{ji}} = \overline{B_{ji}}\,\overline{A_{kj}}$ となるからです。これが、具体的な行列型が transpose と並んで adjoint を提供している理由です。転置が反変関手 $\mathbf{Mat}_R \to \mathbf{Mat}_R$ となるのは可換な場合だけです。

設計上の判断

加群ではなくアーベル群

問題。 「ベクトル」を表す自然な trait は $R$ 加群ですが、加群には 2 つの台集合、すなわちベクトルとスカラーがあります。

選択肢。 (a) スカラー型を固定した（たとえば Double）Module trait。(b) 関連型としてスカラー型を持つ trait。(c) 加法構造だけを扱い、スカラーは具体的なメソッドに任せる。

決定。 (c)。AdditiveVector と AdditiveMatrix はアーベル群 $(V, +, -)$ を表し、スカラーについては何も述べません。

理由。 MoonBit の trait は Self パラメータを 1 つしか持たず、関連型もないため、選択肢 (b) は表現できません。選択肢 (a) では、Int、BigInt、あるいはユーザー定義のスカラー型上の加群がすべて表現できなくなり、さらに 1 つの浮動小数点型を公開 trait にハードコードすることになります。そのためスカラー倍は、スカラー型がわかっている具体的な型のメソッド（scale、left_scale、right_scale）にとどめています。

ベクトル trait にゼロはない

実行時に形状が決まるベクトル型には、単一のゼロがありません。 $(R^n, +)$ の単位元は $n$ に依存します。trait メソッド zero() -> Self を設けると、長さを 1 つ選ばなければならなくなります。そのため trait は Add、Neg、Sub だけを要求し、群の法則は名前付きのゼロを使わずに述べます。

(u + (-u)) + v = v \quad\text{for all } v \text{ of the same length as } u .

元 $u + (-u)$ が、ちょうどよい長さのゼロになります。単一のゼロが存在するスカラー型については、luna-generic の上流 Zero trait を引き続き使えます。

独立した層としての Hadamard 積

ベクトル型の Mul は Hadamard 積、内積、外積（クロス積）のいずれをも意味しえますが、閉じているのは最初のものだけです（すべての $n$ について $R^n \times R^n \to R^n$ ）。VecMulVector はその意味を Hadamard 積に固定し、AdditiveVector の上に置かれます。そのため、ベクトルを足すだけのアルゴリズムが積を持たない型を排除することはありません。内積のようなスカラー値の積は、ベクトルの圏の外への写像 $R^n \times R^n \to R$ であり、構造 trait ではなく具体的なバックエンドのメソッド（DenseVector::dot）です。

行列の乗算は部分的であり、型ではなくドキュメントで扱う

問題。 実行時に形状が決まる行列型では、* は合成可能な形状に対してのみ定義されます。Self を返す trait メソッドは失敗を報告できません。

選択肢。 (a) 最小の行列 trait に Mul を入れる。(b) 形状を型にエンコードする。(c) 形状・転置・加算を乗算から分離し、定義域の外での振る舞いは各実装にドキュメント化させる。

決定。 (c)。TransposeMatrix と AdditiveMatrix は乗算なしでも有用であり、MatMulMatrix は独立した層で、その実装者は失敗時の振る舞いをドキュメント化しなければなりません。固定の $1 \times 1$ や $3 \times 3$ 行列のように形状が静的にわかっている型は、全域的な * で MatMulMatrix を実装します。

理由。 MoonBit には型レベルの自然数がないため、動的な次元に対して (b) は使えません。選択肢 (a) では、乗算を一切しない行列型も含め、すべての行列型に部分的な演算を押し付けることになります。Result を返す検査付きの乗算は、エラー型が固定されている具体的な型（@immut.Matrix::matmul）で引き続き利用できます。

転置は `Self` を返す

TransposeMatrix::transpose は同じ型を返します。転置はあらゆる形状で全域的なので閉じた演算であり、trait メソッドに適しています。結果が実体化されたものかビューかは契約に含まれません。対照的に、container 層の転置は異なるターゲット型を生成することがあります。

正しさと不変条件

trait は空か、観測メソッドを 1 つ持つだけなので、正しさは各実装の性質です。実装が約束する法則は次のとおりです。

レベル	法則（両辺が定義されるとき）
`AdditiveVector`, `AdditiveMatrix`	アーベル群の法則。演算は形状を保つ
`VecMulVector`	$R^n$ は `+` と `*` について環: 結合律、分配律
`TransposeMatrix`	$\operatorname{shape}(A^{\mathsf T}) = \operatorname{shape}(A)$ の反転。 $(A^{\mathsf T})^{\mathsf T} = A$
`AdditiveMatrix`	$(A + B)^{\mathsf T} = A^{\mathsf T} + B^{\mathsf T}$
`MatMulMatrix`	$(AB)C = A(BC)$ 。分配律。スカラーが可換なら $(AB)^{\mathsf T} = B^{\mathsf T}A^{\mathsf T}$

具体的なスカラー型については 2 つの注意点があります。

固定幅整数。 Int は $\mathbb{Z}$ ではなく環 $\mathbb{Z}/2^{32}\mathbb{Z}$ です。ラップアラウンドする加算と乗算は $\mathbb{Z}/2^{32}\mathbb{Z}$ の環演算そのものなので、上の法則はその環で厳密に成り立ちます。オーバーフローする行列積も、 $2^{32}$ を法とした正しい積です。

浮動小数点。 Float と Double は環ではありません。丸めのもとでは加算が結合的でないからです。法則は丸め誤差の範囲でしか成り立たず、テストでは許容誤差付きで比較しなければなりません。任意の順序で計算した長さ $n$ の内積について、標準モデル $\mathrm{fl}(x \circ y) = (x \circ y)(1 + \delta)$ , $|\delta| \le u$ から、要素ごとに

\big|\mathrm{fl}(AB) - AB\big| \le \gamma_n\, |A|\,|B|, \qquad \gamma_n = \frac{n u}{1 - n u},

が得られます。ここで Double では $u = 2^{-53}$ です。²2 N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002, §3.1 および §3.5。この限界は、1 つの内積に含まれる $n - 1$ 回の加算と $n$ 回の乗算のそれぞれにモデルを適用することで得られます。この限界を 2 回適用すると、 $\mathrm{fl}((AB)C)$ と $\mathrm{fl}(A(BC))$ はどちらも正確な積から $2\gamma_k\,|A|\,|B|\,|C| + O(u^2)$ 以内（ $k = \max(n, p)$ ）にあるので、2 つの評価順序の差はその 2 倍まで生じえます。許容誤差に基づく比較が許容すべきなのは、この程度の大きさです。

却下した代替案

VectorSpace または Module trait。 スカラー型を固定せずにスカラーとの対応を表現できるようになるまでは却下です。最初の判断を参照してください。
形状に依存する zero(rows, cols) メソッド。 上流の Zero trait と意味の異なる重複になり、しかも形状を知らない汎用コードには役立ちません。
構造 trait としての内積とノルム。 これらはスカラーへの写像であり、スカラー型に依存し（Double への値をとるノルムは Int への値をとるノルムではありません）、妥当な選択肢が複数あります。バックエンドのメソッドのままにします。
すべての演算を持つ 1 つの Matrix trait。 部分的な乗算や Hadamard 積を、それらを持たない型にまで押し付けることになります。

境界

algebra が定義しないものは次のとおりです。スカラー trait（luna-generic と arithmetic から来ます）、ストレージ・要素アクセス・変更・構築（container の担当です）、検査付きのエラー報告、分解、ソルバー、ノルム、内積。また、具体的な @immut 型と @mutable 型に対する実装も提供しません。それは backends/default が、自身が所有するラッパー型を通じて行います。

積の順序を反転させる対合 $x \mapsto \bar{x}$ 、つまり $\overline{xy} = \bar{y}\,\bar{x}$ を満たすものを備えた環では、共役転置 $A^{*} = \overline{A}^{\mathsf T}$ は可換性なしでも $(AB)^{*} = B^{*} A^{*}$ を満たします。各項が $\overline{A_{kj} B_{ji}} = \overline{B_{ji}}\,\overline{A_{kj}}$ となるからです。これが、具体的な行列型が transpose と並んで adjoint を提供している理由です。 ↩
N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002, §3.1 および §3.5。この限界は、1 つの内積に含まれる $n - 1$ 回の加算と $n$ 回の乗算のそれぞれにモデルを適用することで得られます。 ↩