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 的 trie 实现的持久化向量。替换一个元素时，会复制从根到叶子的路径并共享其他所有节点，因此

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

而旧矩阵仍然有效。读取的开销同样是 $O(\log_{32} N)$；最多七层即可覆盖 $2^{32}$ 个元素。整矩阵运算以 $O(N)$ 重建 trie。

用反复平方计算矩阵幂

对半环上的方阵，把 $k = \sum_t b_t 2^t$ 写成二进制形式。于是

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

这最多需要 $\lfloor \log_2 k \rfloor$ 次平方和同样多次额外乘积，而不是 $k - 1$ 次乘积。乘积的重排之所以合法，是因为矩阵乘法在任意半环上都满足结合律（见 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}$ 中也是精确的，因为回绕算术正是该商环的环算术。

无分数行列式

域上的高斯消元把 $\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}.$$

由此得到两个推论。其一，递推中的除法是精确的：分子是前一个主元的倍数，因此在 $\mathbb{Z}$ 这样的整环上，商仍在该整环中。其二，最后一个值就是完整的行列式，$a^{(n-1)}_{n-1,n-1} = \det A$。该恒等式的证明见下方附件。

Bareiss 消元：精确性与正确性

选主元可以融入而不破坏精确性。$a^{(k)}$ 的第 $i$ 行只依赖于 $A$ 的第 $i$ 行以及第 $0, \dots, k-1$ 行，因此在第 $k$ 步之前交换两个索引 $\ge k$ 的行，等同于对交换了这两行的 $A$ 运行算法，这会使最终结果乘以 $-1$。如果主元下方整列都为零，这些子式全部为零，各行线性相关，于是 $\det A = 0$。

每个中间值都是 $A$ 的子式，因此 Hadamard 不等式为它们全部给出了界：

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

（对所涉及的 $k + 1$ 行取积）。递推的分子是两个此类子式乘积之差，因此其界为该乘积平方的两倍。这为 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$ 时是沿第一行的余子式展开；$n = 4$ 时是沿前两行的 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$ 子式的乘积。这些公式只使用环运算，因此对任意环都是精确的，完全不需要除法。

惰性矩阵

MatrixFn 是一个二元组 $(\text{shape}, f)$，其中 $f : [r] \times [c] \to T$。运算即函数复合：$\mathtt{map}(g)$ 是 $g \circ f$，转置是 $f \circ \mathrm{swap}$，乘积是

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

按需求值。没有任何缓存，因此乘积中一个元素的开销等于内维度乘以各因子元素的开销。于是，$n \times n$ 矩阵上深度为 $d$ 的乘积树每个元素的开销为 $O(n^{d})$：惰性幂构造起来很便宜，读取却很昂贵。

设计决策

持久化向量存储

选项。 (a) 每次更新复制一个 Array，(b) 持久化 trie，(c) 只用函数。决定。 Matrix 和 Vector 采用 (b)，(c) 作为 MatrixFn 单独提供。理由。 复制数组会使 set 的开销为 $O(N)$；trie 使其为 $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 只有一种表示。
• 在 MatrixFn 中缓存。 这会让纯值持有可变状态，并悄悄改变开销模型。请改用 Matrix::make 物化。

边界

immut 不提供原地更新、视图、逆矩阵、分解、特征值、统计或基于容差的谓词；这些属于 mutable。它不检测整数溢出，不提供受检的 MatrixFn API，自身也不实现 algebra trait（由 backends/default 中的包装类型实现）。

Footnotes

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