mutable 设计

设计目标

mutable 是仓库中面向执行的那一半。它把矩阵存放在一个扁平的行主序数组中，允许调用者原地更新并通过实时视图操作，并实现了浮点数值例程：行列式、逆矩阵、秩、行化简、Cholesky 分解、对称矩阵特征值以及幂法。公开 API 在可能的地方仍然读起来像值运算；修改仅限于返回 Unit 的方法和视图，因此调用者可以从签名判断一次调用是否会改变其接收者。

数学背景

本文中，$u$ 表示单位舍入误差（Double 为 $2^{-53}$，Float 为 $2^{-24}$），$\mathrm{fl}(x \circ y) = (x \circ y)(1 + \delta)$ 且 $|\delta| \le u$，$\gamma_n = n u / (1 - n u)$。矩阵之间的不等式按元素成立，$|A|$ 表示由绝对值构成的矩阵。

矩阵乘积

$(AB)_{ik} = \sum_j a_{ij} b_{jk}$ 需要 $rcn$ 次乘加。各内核的求和顺序不同（展开的内核每步计算四个部分积；对至少 $4 \times 16 \times 16$ 的乘积，会把 $B$ 的列打包复制），但每种顺序都满足

$$\big|\mathrm{fl}(AB) - AB\big| \le \gamma_n |A|\,|B| ,$$

因此不同目标上的结果在这一精度内一致，而不是逐位相同。

部分选主元的 LU 分解

对方阵 $A$，部分选主元的高斯消元计算出置换矩阵 $P$、单位下三角矩阵 $L$ 和上三角矩阵 $U$，使得

$$PA = LU .$$

在第 $k$ 步，它选取 $|a^{(k)}_{pk}|$ 最大的行 $p \ge k$，将其交换到位置 $k$，并对 $i > k$ 存储乘数 $l_{ik} = a^{(k)}_{ik} / a^{(k)}_{kk}$，更新 $a^{(k+1)}_{ij} = a^{(k)}_{ij} - l_{ik} a^{(k)}_{kj}$。选主元保证 $|l_{ik}| \le 1$。开销为 $\tfrac23 n^3$ flops。

行列式。 对 $PA = LU$ 两边取行列式，利用 $\det L = 1$ 以及 $s$ 次交换下的 $\det P = (-1)^{s}$，

$$\det A = (-1)^{s} \prod_{k} u_{kk} .$$

求解。 $Ax = b$ 化为 $L y = P b$（前代）和 $U x = y$（回代），每个右端项各需 $n^2$ flops。逆矩阵就是对 $I$ 的 $n$ 列求解：$\tfrac23 n^3 + n \cdot 2n^2 = \tfrac83 n^3$ flops。

稳定性。 计算得到的因子满足向后误差界（Wilkinson；Higham，定理 9.3）¹1 N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002，第 9 章（LU）、第 10 章（Cholesky）和第 8 章（三角方程组）。

$$\hat L \hat U = P(A + \Delta A), \qquad |\Delta A| \le \gamma_n |\hat L|\,|\hat U| .$$

由 $|l_{ik}| \le 1$ 可得 $\lVert \Delta A \rVert_\infty \le n \gamma_n \rho_n \lVert A \rVert_\infty$，其中 $\rho_n = \max_{i,j,k} |a^{(k)}_{ij}| / \max_{i,j} |a_{ij}|$ 是增长因子。部分选主元保证 $\rho_n \le 2^{n-1}$；只有刻意构造的矩阵才会达到该界，实践中 $\rho_n$ 很小，因此该方法在实践中是向后稳定的。求解的前向误差则由条件数决定：$\lVert \hat x - x \rVert / \lVert x \rVert \lesssim \kappa(A)\, n \gamma_n \rho_n$。

小规模行列式的闭式公式

当 $n \le 4$ 时，行列式用公式而非 LU 计算：$ad - bc$ 规则；$n = 3$ 时沿第一行做余子式展开；$n = 4$ 时沿前两行做 Laplace 展开，得到六个互补 $2 \times 2$ 子式的乘积，

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

它们不使用除法，也不做选主元判断，从而对极小的矩阵避免了容差检验。它们在 LU 意义下不是向后稳定的：对病态输入，$ad - bc$ 中的相消可能丢失全部相对精度，这与行列式本身在那里是病态的完全一致。

秩与简化行阶梯形

rank 在副本上做部分选主元的消元，并统计绝对值不小于容差的主元个数。这是带绝对阈值 $\tau$ 的数值秩：即 $\ge \tau$ 的主元个数。对于主元与 $\tau$ 明显分离的矩阵，它是精确的，否则依赖于缩放。（奇异值分解给出可靠的数值秩 $\#\{\sigma_i > \tau\}$；目前未实现。）

reduce_row_elimination 是 Gauss–Jordan 消元：每个主元行被缩放使主元为一，并将主元列的上方和下方都消为零。开销约为 $r c \min(r, c)$ flops，原地进行。

Cholesky 分解

对称正定（SPD）矩阵有唯一的分解 $A = L L^{\mathsf T}$，其中 $L$ 为下三角矩阵且 $l_{jj} > 0$。对 $i \ge j$ 比较 $A = L L^{\mathsf T}$ 的元素，

$$a_{ij} = \sum_{k=0}^{j} l_{ik} l_{jk} \quad\Longrightarrow\quad l_{jj} = \sqrt{a_{jj} - \sum_{k<j} l_{jk}^2}, \qquad l_{ij} = \frac{1}{l_{jj}} \Big(a_{ij} - \sum_{k<j} l_{ik} l_{jk}\Big) .$$

代码逐行计算这些值（Cholesky–Banachiewicz 顺序），需 $\tfrac13 n^3$ flops。第 $j$ 步的被开方数等于不选主元的高斯消元的主元，即顺序主子式之比 $\det M_{j+1} / \det M_j$。由 Sylvester 判据，$A$ 是 SPD 当且仅当所有顺序主子式都为正，因此分解恰好对 SPD 矩阵成功；这就是 is_positive_definite 通过尝试分解来实现的原因。Cholesky 不需要选主元：由 $a_{jj} = \sum_k l_{jk}^2$ 可知每个 $|l_{jk}| \le \sqrt{a_{jj}}$，因此元素不会增长，计算得到的因子满足 $\hat L \hat L^{\mathsf T} = A + \Delta A$，且 $|\Delta A| \le \gamma_{n+1} |\hat L|\,|\hat L^{\mathsf T}|$（Higham，定理 10.3）。

对称特征值问题

对实对称矩阵 $A$，谱定理给出 $A = Q \Lambda Q^{\mathsf T}$，其中 $Q$ 正交，$\Lambda$ 为实对角矩阵。eigen 分两个阶段计算它。

Householder 三对角化。 Householder 反射子 $H = I - 2 v v^{\mathsf T} / v^{\mathsf T} v$ 是对称正交的，可以把一个向量映射为某个坐标向量的倍数。从两侧应用 $n - 2$ 个反射子，

$$Q_1^{\mathsf T} A Q_1 = T, \qquad Q_1 = H_1 H_2 \cdots H_{n-2},$$

其中 $T$ 是对角线为 $d$、次对角线为 $e$ 的对称三对角矩阵。代码显式累积 $Q_1$；合计约 $\tfrac83 n^3$ flops。在构造反射子之前，每一行都按其绝对值之和缩放，以避免 $\lVert x \rVert_2$ 的上溢和下溢。

带 Wilkinson 位移的隐式 QL。 三对角矩阵 $T$ 通过平面旋转对角化。在从 $l$ 开始的未约化块上进行每轮扫描之前，位移 $\sigma$ 取前导 $2 \times 2$ 块 $\begin{pmatrix} d_l & e_l \\ e_l & d_{l+1} \end{pmatrix}$ 中更接近 $d_l$ 的特征值。令 $g = (d_{l+1} - d_l)/(2 e_l)$，该块的特征值为

$$\lambda_\pm = d_l + e_l \big(g \pm \sqrt{g^2 + 1}\big),$$

其中更接近 $d_l$ 的那个是

$$\sigma = d_l + e_l\big(g - \operatorname{sgn}(g)\sqrt{g^2+1}\big) = d_l - \frac{e_l}{g + \operatorname{sgn}(g)\sqrt{g^2 + 1}} ,$$

代码使用的第二种形式避免了相消。一轮扫描应用 Givens 旋转追赶由此产生的凸起，并用同样的旋转更新 $Q$。当满足以下条件时，次对角元被视为零：

$$|e_m| \le \tau\,\big(|d_m| + |d_{m+1}| + 1\big),$$

这一检验对较大的对角元是相对的，在零附近是绝对的。实践中，对称三对角矩阵在 Wilkinson 位移下呈三次收敛，在有限输入上从未观察到失败；尽管如此，代码在某个特征值经过 60 轮扫描后仍会中止。整个过程是正交变换的乘积，是向后稳定的：计算出的特征值是 $A + \Delta A$ 的精确特征值，其中 $\lVert \Delta A \rVert_2 = O(u)\lVert A \rVert_2$，因此由 Weyl 不等式

$$|\hat\lambda_i - \lambda_i| \le \lVert \Delta A \rVert_2 = O(u)\,\lVert A \rVert_2 .$$

特征向量的精度与 $u\lVert A \rVert / \text{gap}$ 成比例，其中 gap 是到最近的其他特征值的距离。

$2 \times 2$ 情形。 对 $A = \begin{pmatrix} a & b \\ c & d \end{pmatrix}$，令 $m = (a + d)/2$，特征多项式 $\lambda^2 - (a + d)\lambda + (ad - bc)$ 给出

$$\lambda_{1,2} = m \pm \sqrt{m^2 - (ad - bc)} .$$

当 $b \ne 0$ 时，向量 $(b, \lambda - a)^{\mathsf T}$ 是特征向量：

$$\begin{pmatrix} a - \lambda & b \\ c & d - \lambda \end{pmatrix} \begin{pmatrix} b \\ \lambda - a \end{pmatrix} = \begin{pmatrix} 0 \\ bc - (\lambda - a)(\lambda - d) \end{pmatrix} = 0 ,$$

因为由特征方程有 $(\lambda - a)(\lambda - d) = \lambda^2 - (a + d)\lambda + ad = bc$。代码返回这些向量时不做归一化。当 $|\lambda_2| \ll |\lambda_1|$ 时，减法 $m - \sqrt{\cdot}$ 会发生相消，$\lambda_2$ 丢失相对精度；$\lambda_2 = \det A / \lambda_1$ 才是稳定的替代方案。

幂法

从 $x_0 = (1, \dots, 1)$ 出发（若 $A x_0 = 0$ 则从某个坐标向量出发），该方法迭代

$$y = A x_k, \qquad x_{k+1} = y / \lVert y \rVert_\infty, \qquad \lambda_k = \frac{x_k^{\mathsf T} A x_k}{x_k^{\mathsf T} x_k} ,$$

并在 $\lVert A x_k - \lambda_k x_k \rVert_\infty \le \tau$ 时停止。若 $A$ 可对角化，特征值满足 $|\lambda_1| > |\lambda_2| \ge \dots$，且 $x_0$ 在特征向量 $v_1$ 方向上有分量 $c_1 \ne 0$，则

$$A^{k} x_0 = \lambda_1^{k}\Big(c_1 v_1 + \sum_{i \ge 2} c_i \big(\tfrac{\lambda_i}{\lambda_1}\big)^{k} v_i\Big),$$

因此方向以比率 $|\lambda_2 / \lambda_1|$ 线性收敛；对对称矩阵 $A$，Rayleigh 商以比率 $|\lambda_2 / \lambda_1|^2$ 收敛。当 $|\lambda_1| = |\lambda_2|$ 且 $\lambda_1 \ne \lambda_2$（例如 $\pm 1$）时，方向会振荡，残差检验永远无法通过；当 $A$ 幂零时，迭代向量会变为零；这两种情况都返回 None。

统计量

variance 是用两遍算法计算的总体方差：先求均值 $\bar a$，再求 $\tfrac1N \sum (a_i - \bar a)^2$。单遍公式 $\tfrac1N \sum a_i^2 - \bar a^2$ 在数据均值大、离散度小时会让两个几乎相等的数相减，甚至可能返回负值；两遍形式对非负项求和，其舍入误差相对于方差本身（Chan、Golub 和 LeVeque，1983）。计数 $N$ 以 T 中 1 的累加和得到，对 Double 在 $2^{53}$ 个元素以内、对 Float 在 $2^{24}$ 个元素以内是精确的。

转置视图与乘积

Transpose::mul 在被包装的矩阵上复用矩阵内核，把 $A^{\mathsf T} B^{\mathsf T}$ 计算为 $(BA)^{\mathsf T}$。逐元素来看，

$$\big(A^{\mathsf T} B^{\mathsf T}\big)_{ik} = \sum_j a_{ji} b_{kj}, \qquad \big((BA)^{\mathsf T}\big)_{ik} = \sum_j b_{kj} a_{ji} ,$$

二者在标量可交换时相等。所有实现了 Tolerance 的标量类型（Float、Double）都可交换，但 Transpose::mul 只要求 AddMonoid + Mul；对非交换标量类型，结果是顺序错误的乘积（见 algebra 设计）。

设计决策

扁平的行主序存储

选项。 行数组的数组、持久化结构，或一个扁平数组。决定。 一个 Array[T]，元素 $(i, j)$ 位于 $i c + j$，四个目标共用。理由。 它提供 $O(1)$ 访问，每个坐标只做一次边界检查；消元和乘积的内层循环可以访问连续的行；行视图和列视图零开销。from_array 直接接管调用者的数组，使大型输入无需复制；代价是别名，API 文档对此作了说明。

用视图代替副本

row_view、col_view 和 to_transpose 在 $O(1)$ 内返回实时视图。视图是一个（矩阵，索引）对或一个包装，因此写入会落到共享存储上，无需同步。物化总是显式进行的（to_vector、materialize、transpose）。

语义共享的目标专用内核

该包为矩阵、LU、视图和转置代码在每个目标上各保留一个源文件。它们只在为各后端选择的循环结构上有所不同（展开、打包、在索引计算中避免除法）；公开语义（包括边界和错误行为）完全相同，测试在全部四个目标上运行。由于求和顺序不同，浮点结果的最后几位可能有差异。

基于容差的判断

消元必须判断计算出的主元何时“为零”。该包对 Double 和 Float 使用同一个绝对阈值 $\tau$ = Tolerance::tolerance() = $10^{-11}$，用法如下：

例程	测试
LU（$n \ge 5$ 时的 determinant、inverse、is_invertible）	主元 $\lvert u_{kk}\rvert < \tau$ 即视为奇异
rank	剩余最大 $\lvert a_{ik}\rvert < \tau$ 即无主元
reduce_row_elimination	$\le \tau$ 的元素置零
cholesky_decomposition	被开方数 $\le \tau$ 即非正定
is_symmetric、快速路径检测（单位、对角、置换、三角）	$\lvert a_{ij} - b_{ij}\rvert \le \tau$
eigen 收缩	$\lvert e_m\rvert \le \tau(\lvert d_m\rvert + \lvert d_{m+1}\rvert + 1)$
power_method	残差 $\lVert Ax - \lambda x\rVert_\infty \le \tau$

绝对阈值简单且可预测，但不具有尺度不变性。把 $A$ 缩放 $10^{-12}$ 倍会使每个主元都低于 $\tau$，于是一个条件极好的矩阵会被报告为奇异；缩放 $10^{12}$ 倍则会让数值上奇异的矩阵通过检验。调用这些例程之前，请先把数据缩放到量级为一。对于 Float，$\tau = 10^{-11}$ 远小于其单位舍入误差 $\approx 6 \times 10^{-8}$，因此这些检验实际上只是在检查严格的零，接近奇异的 Float 矩阵不会被检测出来。该 trait 是封闭的（pub 而非 pub(open)），因此只有这两个实例；感知尺度的容差策略是未来的工作，并且将是一次破坏性变更。

快速路径

inverse 能识别单位矩阵（返回副本）、对角矩阵（对对角线求逆）和置换矩阵（返回转置，因为对列为互不相同坐标向量的矩阵有 $P^{\mathsf T} P = I$）。当 $n \ge 5$ 时，determinant 能识别三角矩阵并直接计算对角线乘积。每项检测的开销为 $O(n^2)$，可省去一次 $O(n^3)$ 的分解。这些检测使用容差，因此在 $\tau$ 范围内为对角的矩阵会被视为严格对角。

不改动内核的受检形式

每个受检方法都先校验（是否方阵、指数符号、是否非空、长度），再调用其非受检对应版本，后者保留原来的中止或 Option 行为。这里没有受检的矩阵乘积：* 会校验并中止，unchecked_matmul 则完全不校验。与 @immut.Matrix::matmul 的这种不对称是已知的；增加受检的 matmul 将是一次新增，而不是修改。

仅支持对称特征值

一般实矩阵可能有复特征值，而对实数 T，返回 Vector[T] 的函数无法表示它们。因此 eigen 只接受对称矩阵（其特征值全为实数，且存在标准正交特征基），对其他矩阵则中止。

正确性与不变量

• 存储。 任何时候都有 data.length() == row * col；每个公开访问器都分别检查行和列。
• 返回值的方法不做修改。 只有返回 Unit 的方法、视图写入和 reduce_row_elimination 会改变矩阵。
• 受检/非受检定律。 只要前置条件成立，就有 x.f() == Ok(x.unchecked_f())；inverse 返回 Err(SingularMatrix) 当且仅当 unchecked_inverse 返回 None。
• 行列式一致性。 当 $n \ge 5$ 且矩阵在 $\tau$ 范围内不是三角矩阵时，determinant 返回 0 当且仅当 LU 分解报告了低于 $\tau$ 的主元，这也正是 is_invertible 返回 false、inverse 失败的情形。当 $n \le 4$ 时，determinant 使用闭式公式，而 is_invertible 仍使用 LU，因此行列式非零但极小的矩阵可能“不可逆”，同时 determinant 非零。
• 残差保证。 cholesky_decomposition 和 eigen 返回的因子，其重构误差为 $O(u)\lVert A \rVert$；power_method 只返回残差不超过 $\tau$ 的特征对。
• 复杂度。 *：$rcn$；determinant、inverse：$O(n^3)$；rank、reduce_row_elimination：$O(rc\min(r, c))$；cholesky_decomposition：$n^3/3$；eigen：$O(n^3)$；power_method：每次迭代 $O(n^2)$；统计：$O(rc)$。

被否决的方案

• 相对容差或按范数缩放的容差。 更稳健，但会改变现有调用者的结果；推迟到可以显式传入容差策略时再考虑。
• 返回复特征值。 这会让本包依赖复数类型，并改变实对称输入的签名。
• 写时复制的视图。 这会隐藏写入的开销，并破坏视图存在的意义——原地契约。
• 单一的可移植内核。 实测在某些目标上更慢；按目标拆分的文件以代码体积换取速度，同时共享同一份规范。

边界

mutable 不提供以公开方法形式针对右端项的线性求解（只有逆矩阵），也不提供 QR、SVD 或最小二乘求解器、非对称矩阵的特征值、稀疏存储、条件数估计或感知尺度的容差。它自身不实现 algebra trait；backends/default 为此对它进行了包装。回归或优化这类领域工作流属于下游包。

Footnotes

1. N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002，第 9 章（LU）、第 10 章（Cholesky）和第 8 章（三角方程组）。 ↩