algebra 设计

设计目标

algebra 为泛型线性代数代码提供一套描述整体向量和矩阵对象的词汇，并如实反映每个对象能做什么。只有满足其定律的类型才可以实现某个 trait，而算法应当能够只要求它真正用到的结构：形状、阿贝尔群、Hadamard 环、转置或矩阵乘积。因此，该包定义的是一组按包含关系排序的小 trait，而不是单个 Matrix 或 VectorSpace trait。API 页面列出了这些 trait；本页解释它们所编码的数学以及由此得出的设计选择。

数学背景

模与向量空间

设 $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)$。

矩阵构成范畴

矩阵乘法并不是单个集合上的运算。它的形状规则

$$(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}$$

两个和式逐项相等，当且仅当 $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$-模，但模有两个载体：向量和标量。

选项。 (a) 一个标量类型固定（例如 Double）的 Module trait。(b) 一个带关联标量类型的 trait。(c) 只描述加法结构，把标量留给具体方法。

决定。 (c)：AdditiveVector 和 AdditiveMatrix 只陈述阿贝尔群 $(V, +, -)$，不涉及标量。

理由。 MoonBit 的 trait 只有一个 Self 参数，没有关联类型，因此无法表达选项 (b)。选项 (a) 会使 Int、BigInt 或用户标量类型上的模都无法表示，并把某一种浮点类型硬编码进公开 trait。因此，标量乘法仍是具体类型的方法（scale、left_scale、right_scale），在那里标量类型是已知的。

向量 trait 中没有零元

运行时确定形状的向量类型没有唯一的零元：$(R^n, +)$ 的单位元取决于 $n$。trait 方法 zero() -> Self 不得不选定一个长度。因此这些 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$，超出了向量范畴；它们是具体后端的方法（DenseVector::dot），而不是结构 trait。

矩阵乘法是部分运算，以文档说明而非类型编码

问题。 在运行时确定形状的矩阵类型上，* 只对可复合的形状有定义。返回 Self 的 trait 方法无法报告失败。

选项。 (a) 把 Mul 放进最小的矩阵 trait。(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 要么为空，要么只有一个观察方法，因此正确性是每个实现自身的性质。实现所承诺的定律是：

级别	定律（只要两边都有定义）
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}$

对具体的标量类型有两点需要注意。

定宽整数。 Int 是环 $\mathbb{Z}/2^{32}\mathbb{Z}$，而不是 $\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。对一次内积中的 $n - 1$ 次加法和 $n$ 次乘法分别应用该模型即可得到此界。 两次应用该界，$\mathrm{fl}((AB)C)$ 与 $\mathrm{fl}(A(BC))$ 各自与精确乘积相差不超过 $2\gamma_k\,|A|\,|B|\,|C| + O(u^2)$，其中 $k = \max(n, p)$，因此两种求值顺序的结果最多可相差该量的两倍。这就是基于容差的比较必须允许的量级。

被否决的方案

• VectorSpace 或 Module trait。 在能够不固定标量类型地表达标量关联之前予以否决；见第一个决策。
• 依赖形状的 zero(rows, cols) 方法。 它会以不同含义重复上游的 Zero trait，而且仍无法服务于不知道形状的泛型代码。
• 把内积和范数作为结构 trait。 它们映射到标量，依赖于标量类型（取值于 Double 的范数并不是取值于 Int 的范数），而且有多种合理选择。它们仍作为后端方法。
• 包含所有运算的单个 Matrix trait。 它会把部分乘法和 Hadamard 积强加给不具备它们的类型。

边界

algebra 不定义标量 trait（来自 luna-generic 和 arithmetic），不涉及存储、元素访问、修改或构造（属于 container），也不涉及受检错误报告、分解、求解器、范数或内积。它不为具体的 @immut 和 @mutable 类型提供实现；这由 backends/default 通过自有的包装类型完成。

Footnotes

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 的原因。 ↩

2. N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002，§3.1 与 §3.5。对一次内积中的 $n - 1$ 次加法和 $n$ 次乘法分别应用该模型即可得到此界。 ↩