backends/default 設計

設計目標

algebra の trait は、実際の密な型がそれを実装していなければ役に立ちません。backends/default は、具体的な @mutable 型と @immut 型をラップすることでその参照バックエンドを提供しつつ、それらの具体的なパッケージが実験的な algebra 層に一切依存しないようにしています。これはエコシステムの中心ではなく、1 つのバックエンドにすぎません。汎用アルゴリズムは trait に依存し、このパッケージはそれを満たす方法の 1 つです。

数学的背景

インスタンスは法則の根拠である

ある型に @algebra.MatMulMatrix を実装することは、algebra の設計 に挙げた法則、つまり定義される範囲での * の結合律と + に対する分配律を主張することです。ラッパーでは、これらの法則はラップされた型から受け継がれます。どの演算子も「アンラップし、内側の演算子を適用し、ラップする」として定義されているからです。

$$\mathrm{wrap}(A) \cdot \mathrm{wrap}(B) = \mathrm{wrap}(A \cdot B) .$$

したがって wrap はすべての演算について準同型であり、内側の型で成り立つ等式はラッパーでも成り立ちます。

密行列積の部分性

実行時に形状が決まる密行列では、積は $\{(A, B) : \operatorname{cols}(A) = \operatorname{rows}(B)\}$ 上で定義されます。trait メソッドは Self を返すので、この集合の外ではラッパーは何らかの処理をしなければなりません。内側の型と同様に中断します。これは、MatMulMatrix の契約が実装に明示を求めている、ドキュメント化された実行時の前提条件です。

内積とその丸め誤差

dot は $s_n = \sum_{i=1}^{n} u_i v_i$ を漸化式 $s_0 = 0$, $s_i = s_{i-1} + u_i v_i$ で計算します。浮動小数点では、各ステップがそれまでのすべての項の誤差にさらに因子 $(1 + \delta)$ を掛けるため、標準的な議論から次が得られます。

$$\big|\mathrm{fl}(s_n) - s_n\big| \le \gamma_n \sum_{i=1}^{n} |u_i v_i|, \qquad \gamma_n = \frac{n u}{1 - n u} .$$

したがって相対誤差は、項の符号がそろっているときには小さく、$\sum |u_i v_i| \gg |s_n|$ のとき（桁落ち）には大きくなりえます。matvec はこのような内積を $m$ 個計算するもので、行ごとに同じ限界を受け継ぎます。

設計上の判断

所有するラッパー型

問題。 MoonBit では impl Trait for Type を書けるのは、その trait か型を所有するパッケージだけです。@mutable.Matrix に @algebra.MatMulMatrix を実装するには、algebra（具体的な型を知ってはならない）か mutable（そうすると実験的な algebra 層に依存してしまう）のどちらかで行う必要があります。

決定。 このパッケージで新しい型 DenseMatrix[T]、DenseVector[T] とその不変版を定義し、それぞれを公開フィールド inner を 1 つ持つ構造体とし、ここで trait を実装します。

理由。 ラッパーは trait を実装するパッケージが所有するので規則を満たし、依存の向きは backends/default → algebra, immut, mutable のままです。代償は型における 1 段の間接参照ですが、inner() と from_backend によってコピーなしで取り除けます。

スカラー値の写像はバックエンドのメソッドにする

dot、scale、axpy、matvec は trait ではなくラッパーのメソッドです。これらはスカラー型を明示的に扱いますが、algebra の trait はスカラー型を名指しできません（algebra の設計を参照）。メソッドにしておくことで、各バックエンドが独自の制約を選べます。たとえば可変ベクトルの dot は AddMonoid + Mul、不変ベクトルでは Zero + Add + Mul です。

axpy は新しいベクトルを返す

BLAS のルーチン axpy は $y \leftarrow a x + y$ をインプレースで更新します。ここでは x.axpy(a, y) はどちらのラッパーでも $x a + y$ を新しいベクトルとして返すため、可変と不変のバックエンドが値を返す同じ契約を共有します。インプレースの更新は内側の @mutable.Vector で引き続き利用できます。

スカラーは右から掛ける

scale は right_scale を使い、$v_i a$ を計算します。可換なスカラーではこれは $a v_i$ と同じですが、非可換なスカラーではこの選択が結果に現れるため、隠さずにドキュメント化しています。

正しさと不変条件

• 準同型。 すべての演算子について inner(a op b) == inner(a) op inner(b) が成り立つので、ラッパーはラップされた型の法則をそのまま満たします。
• ラップ時にコピーしない。 from_backend(x).inner() は物理的に x そのものです。可変な内側の値への書き込みはラッパーを通じても見えます。
• 転置。 transpose は実体化します。ビューを返すことはないので、TransposeMatrix の法則 $(A^{\mathsf T})^{\mathsf T} = A$ は値として成り立ちます。
• 計算量。 +、-、scale: $O(n)$。dot: $n$ 回の積和。matvec: $mn$。*: $rcn$ 回の積和。transpose: $O(rc)$ 回のコピー。

却下した代替案

• 1 つの行列型の内部での実行時バックエンドセレクタ。 すべての演算がバックエンドで分岐することになり、どのカーネルが実行されるかが見えなくなります。バックエンドは型で選びます。
• immut と mutable で直接 trait を実装する。 安定した具体的パッケージを実験的な algebra 層に縛り付けることになります。
• インプレースの axpy。 同じ名前に対して 2 つのラッパーで異なる意味論を持たせることになります。

境界

backends/default は新しい trait も新しい数値アルゴリズムも定義しません。mutable の分解、逆行列、統計量には inner() を通じてアクセスします。疎、遅延、静的サイズ、GPU のバックエンドは提供しません。そうしたバックエンドは、これらのラッパーに変換するのではなく、自身の型に algebra の trait を実装すべきです。かつてこのパッケージの隣にあったネイティブ OpenBLAS バックエンドは本リリースで取り下げられ、contrib/openblas_backend に保存されています。