Rust 高性能线性代数库 faer：矩阵分解、性能调优与工程避坑

线性代数是那种"平时不显山露水、一旦用上就要命"的东西。你在做图形学、做物理仿真、做机器学习里的某个小算子、做信号处理，或者只是想让一段拟合代码跑得别那么肉，最后都会绕到矩阵分解、最小二乘、特征值这一类运算上。在 Rust 生态里做这件事，绕不开的一个名字就是faer——一个纯 Rust 实现的高性能线性代数库。它的定位很明确：不依赖外部的 Fortran 计算后端，不用 FFI 去桥接系统里装没装的数学库，整个从内存布局、SIMD 微内核到多线程调度都由 Rust 代码自己完成。这篇文章想聊的是我实际把 faer 塞进项目里之后，对它的设计取舍、API 手感、性能调优手段和踩过的坑的理解，面向的是能用 Rust 写出能编译的代码、但对数值库内部不太熟的开发者，也包括那些正在纠结"到底要不要引入一个第三方矩阵库"的人。读完你应该能独立判断：你的场景该不该上 faer，上了之后怎么配、怎么写、怎么调，以及哪些地方最容易翻车。

1. faer 到底解决了谁的痛点

1.1 纯 Rust 实现的分量有多重

先说清楚"纯 Rust"到底指的是什么，因为这个词经常被当成营销话术。传统的科学计算栈是分层的：上层是 LAPACK 这样的一批 Fortran 子程序，负责分解、求解这些高层算法；下层是 BLAS，负责矩阵乘法、向量点积这类基础算子；再往下是各家针对具体指令集手写的汇编或 intrinsics 内核。这套体系跑了几十年，成熟得很，但它的分发形态对 Rust 项目其实相当不友好。

你的程序想用它们，通常有几条路：链接系统里已经装好的动态库（用户机器上有没有、版本对不对、符号名字对不对，全是运气）、在构建时编译一份静态库（构建环境里得有 Fortran 编译器，交叉编译到某些目标平台会直接卡死）、或者用某个预先打包好的产物（体积、许可、平台覆盖面又成了新问题）。这些麻烦在服务端可能还能忍，一旦到了桌面应用、移动端、嵌入式，或者你想产出一个"下载即用、零外部依赖"的单文件二进制，立刻就变成实打实的负担。

faer 的做法是把整条链路用 Rust 重写一遍，从微内核开始。带来的直接后果是：cargo build是它唯一需要的构建步骤，交叉编译跟普通 Rust 项目没有任何区别，最终产物不引入任何额外的动态库依赖。你在 CI 里给某个少见的目标平台出包，不会因为缺少某个 Fortran 运行时链接失败。这个好处听起来朴素，但真被它救过一次的人会很有感触。

另一层价值在于类型系统的贯穿。整个库的标量类型是泛型参数，f32、f64、单精度复数、双精度复数都能用同一套 API 走通，不需要为每种类型单独选一个函数名或换一套调用约定。当你写的算法需要在精度和性能之间做切换（比如编辑器里实时预览用单精度，导出时用双精度），这一点能省掉大量重复代码。

1.2 和传统计算栈的取舍对照

把话说透一点，faer 不是"免费的午餐"，它是拿生态成熟度换取了集成体验。我做了一张对照表，列的是我在选型时真正在意的那几个维度：

这里我要强调一个常见误解：很多人默认"手写 Rust 一定慢"。这话在五年前成立，现在不太成立。faer 的核心竞争力恰恰是它把传统实现里最值钱的两块东西——分块策略和 SIMD 微内核——认真做了一遍，再加上 Rust 本身在内存安全上的优势，让你不太可能因为越界写坏一个缓冲区而得到一堆莫名其妙的数值。至于性能是不是"接近甚至超过"，这高度依赖矩阵形状、标量类型、线程数和你编译时开的指令集，我不建议你直接采信任何声称的倍数，自己测一遍更靠谱，方法我在第 4 节里给。

1.3 什么情况下我会劝你先别用

实话实说，有几类场景我会建议先缓缓。

第一类是你只需要四则运算级别的矩阵操作。如果你的全部需求就是两个 3×3 矩阵相乘、算个 4×4 逆、做个简单的变换叠加，那引入一个功能齐全的线性代数库属于杀鸡用牛刀，手写一个固定尺寸的小模块，配上[[f64; 3]; 3]这种数组类型，编译器能帮你把循环完全展开，性能未必差，依赖还更干净。这类需求在图形学里极其常见，我见过不少项目为了几个 4×4 乘法拖进一个巨型依赖，纯属浪费。

第二类是你需要的是符号运算或者超大规模稀疏求解。faer 主要面向稠密数值计算，稀疏部分在持续完善中，如果你的主战场是几十万维、百万维的稀疏方程组，或者需要一套完整的稀疏直接法生态，那还是要谨慎评估，先看它的稀疏模块能不能覆盖你的场景，别一上来就押注。

第三类是你对极端病态输入的鲁棒性有极高要求。传统实现里沉淀了大量针对奇异矩阵、极大条件数、非规格化浮点数的防御性处理，这些是几十年真实故障喂出来的经验。faer 在这一块做得不差，但历史积累的厚度上还是有差距。如果你的输入来自不可控的外部数据，做好自己的数值健康检查（残差验证、条件数估计）比依赖库自动兜底更靠谱。

2. 环境配置：从第一行 Cargo 依赖到指令集开关

2.1 依赖声明与 feature 的取舍

基本接入没什么花活，cargo add或者手写一行依赖就行。但 feature 的选择值得说一下，因为它直接影响你的编译时间和运行时行为。

[dependencies] faer = { version = "0.22", default-features = false, features = ["std", "rayon"] }

上面这行是我的常用配置。default-features = false再手动勾选，目的是把控制权拿在手里。几个关键 feature 的含义大致是这样：std提供标准库环境下的分配器和部分依赖标准库的实现，做no_std场景才需要关掉；rayon打开多线程能力，大矩阵运算会用到它；linalg是线性代数算法层的总开关，分解、求解这些都在它后面。如果你只想要矩阵容器和基础的按元素运算，把linalg关掉能省一些编译时间——不过说实话，绝大多数人最后都会用上，不如一开始就开着。

注意：不要把rayon关掉又指望大矩阵跑得快。单线程的 dgemm 在几百阶以上就会被多线程版本甩开，这个差距是数量级边缘的，不是几个百分点。

另外，这个库曾经用过 nightly 专属特性来追求更激进的优化。如果你的项目能接受 nightly 工具链，可以尝试打开相关 feature 看看收益；如果追求稳定、长期跟进的工程环境，就用 stable，功能上不会有缺口，只是上限略低一点点。

版本方面我要提醒一句：这个库还在快速迭代期，API 在不同小版本之间确实发生过调整，尤其是分解器的构造方式（用构造函数还是用矩阵上的便捷方法）。你从网上抄来的示例代码编译不过，第一反应应该是查当前版本的文档，而不是怀疑自己写错了。我后面给的代码也是以某个版本为基准的，你在实际项目里对照一下文档，改动通常就是个方法名的事。

2.2 编译参数：不做这一步性能会差一大截

这是我认为最容易被忽略、也最影响体感的一环。Rust 在 x86-64 平台上默认的基准指令集非常保守，只保证 SSE2 级别可用，因为它要确保产出的二进制在很老的机器上也能跑。但 faer 的微内核里有针对更宽 SIMD 指令（比如 256 位、512 位的向量指令）的实现，这些实现只有在编译器知道你允许使用它们的时候才会被启用。如果你什么都不配，很可能跑的是最窄的那条路径。

最直接的办法是在项目的.cargo/config.toml里加一段：

[build] rustflags = ["-C", "target-cpu=native"]

target-cpu=native让编译器针对当前这台编译机的 CPU 生成代码。这在开发机上很好用，一键拉满。但请务必注意它的坑：用 native 编译出来的二进制，拿到指令集更老的机器上会直接崩溃（非法指令）。如果你是给自己机器用，随便；如果你要分发给用户，绝对不能用 native。

分发场景下的稳妥做法是指定一个明确的微架构等级。Intel 和 AMD 平台在 x86-64-v3 这个等级上已经覆盖得比较广了，它大致对应 AVX2 加 FMA 那一代能力：

[build] rustflags = ["-C", "target-cpu=x86-64-v3"]

同时建议把发布配置也调一下，减少跨 crate 内联的障碍：

[profile.release] lto = "fat" codegen-units = 1

lto = "fat"让链接器跨 crate 做全局优化，能把这层薄薄的封装函数彻底内联掉，代价是编译时间明显变长。codegen-units = 1牺牲并行编译换取更好的优化机会。这两个组合起来，我实测在计算密集的场景里能拿到可观的提升，但在编译时间上是从几十秒涨到一两分钟的量级，所以只建议在 release profile 上开，别动 dev。

2.3 矩阵类型体系：搞清楚所有权的三种形态

faer 的类型设计和 Rust 的借用体系咬合得很紧，理解这一点能让你少写很多没必要的内存分配。

核心是三个东西：Mat<T>是拥有数据的矩阵，内部是一段连续分配的内存加上行列维度；MatRef<'a, T>是对一块矩阵数据的不可变借用视图；MatMut<'a, T>是可变借用视图。它们之间的关系很像Vec<T>、&[T]、&mut [T]。

这个设计的意义在于：几乎所有的运算接口都接受MatRef，而你手上的Mat只要调一次.as_ref()就能借出去，不复制任何数据。反过来，如果你的矩阵数据本来就在一个大的缓冲区里，你完全可以用视图直接切出一块参与运算，避免"先拷出来再算"的浪费。

内存布局是列主序，这点从 Fortran 时代一路继承下来，也是数值计算领域的事实标准。它意味着(i, j)元素存放在j * nrows + i这个位置。对你写代码的影响是：按列遍历比按行遍历对缓存更友好。如果你的算法天然是行优先的，要么调整循环顺序，要么接受一点性能损失，别指望库帮你自动重排。

元素访问用元组下标，看起来很自然：

use faer::Mat; let mut m = Mat::<f64>::zeros(3, 3); m[(0, 0)] = 2.0; m[(1, 1)] = 2.0; m[(2, 2)] = 2.0; println!("{}", m[(1, 1)]);

Index的实现会在越界时 panic，如果你需要不 panic 的版本，用get(i, j)拿Option<&T>。在热循环里逐元素用下标赋值是可以的，但它毕竟带着边界检查，如果一段代码要填几十万个数，用from_fn一次性构造会干净很多，后面会讲到。

3. 核心 API 实操：从建矩阵到解方程

3.1 构造矩阵的几种姿势与选择依据

构造方式的选择取决于你的数据从哪来。最通用的入口是Mat::from_fn，它接受行列数和闭包，把生成逻辑交给你：

use faer::Mat; // 构造一个三对角矩阵，模拟一维热传导的离散化 let n = 1000; let a = Mat::<f64>::from_fn(n, n, |i, j| { if i == j { 2.0 } else if (i as isize - j as isize).abs() == 1 { -1.0 } else { 0.0 } });

这种写法我很喜欢，因为它是按需生成，闭包被逐个调用，不会先分配再填一遍。对于结构化矩阵（对角、三对角、带状、Toeplitz），一行逻辑就能表达清楚，可读性比手动双重循环加一堆if好得多。

需要注意的是，from_fn内部仍然要遍历全部n × n个位置。上面那个例子n = 1000意味着一百万次闭包调用，其中绝大多数返回 0。如果你真的需要稀疏存储，应该考虑稀疏矩阵类型而不是稠密矩阵——稠密矩阵不因为它"看起来是稀疏的"就少占内存。

如果需要从已有的连续缓冲区接手，用切片视图更合适：

let buf: Vec<f64> = vec![0.0; 6]; // 把 buf 按列主序解释成 2×3 的矩阵视图 let view = faer::mat::from_column_major_slice::<f64, _, _>(&buf, 2, 3);

这种免拷贝的视图在和外部接口交换数据的时候特别有用。列主序这一点在这里必须记牢：如果你手上是行优先的数据，要么自己转置一次，要么用带步长的视图接口描述布局，别硬塞。

还有一种情况是逐元素计算。库为矩阵和标量、矩阵和矩阵之间的加减乘提供了运算符重载，写起来像数学公式：

let b = &a * 2.0; // 标量乘 let c = &a + &b; // 逐元素加（这里 a、b 形状需一致）

要注意这里的*在矩阵与矩阵之间是逐元素乘还是矩阵乘，这是个语义陷阱。在很多数值库里*是逐元素，矩阵乘要显式调用函数；不要凭印象写，看文档确认，写错了编译器不会报错，你的结果会悄悄错掉。

3.2 矩阵乘法：三种输入输出形态怎么选

矩阵乘法是性能核心，接口设计也最能体现一个库的取向。faer 里典型的调用形态是这样的：

use faer::linalg::matmul::matmul; use faer::linalg::matmul::Accum; use faer::{Mat, Par}; let m = 512; let k = 512; let n = 512; let lhs = Mat::<f64>::from_fn(m, k, |i, j| ((i + j) % 7) as f64); let rhs = Mat::<f64>::from_fn(k, n, |i, j| ((i * j) % 5) as f64); let mut dst = Mat::<f64>::zeros(m, n); matmul( dst.as_mut(), Accum::Replace, lhs.as_ref(), rhs.as_ref(), 1.0, Par::Rayon(0), );

这段里有几个点值得展开。

Accum::Replace表示输出矩阵先被清零再累加结果，对应数学上的C = A * B。如果你要做的是C = C + A * B，换成Accum::Add就省掉了手动清零。这个区分在迭代算法里很实用——比如幂法求主特征向量、或者分块累加的算法，很多情况下你希望把本次计算的结果叠加到已有结果上，用Add能避免一次额外的读写。

Par::Rayon(0)里的 0 表示让 rayon 用默认线程数，通常就是 CPU 逻辑核心数。传Par::Seq就是强制单线程。什么时候该强制单线程？小矩阵。线程池的调度和任务拆分开销在几百微秒量级，如果你的矩阵乘法本身只要 50 微秒，开多线程大概率是负优化，反而更慢。我的经验分界线是：矩阵规模在百阶以下，先测单线程；几百阶往上，多线程基本稳赢；中间地带必须实测。

还有一个常被忽视的点：输出的dst必须尺寸匹配。这类接口不做自动扩容，你传一个尺寸不对的目标矩阵进去，等着你的就是 panic 或者更隐蔽的错误结果。写的时候把m、n、k三个变量算清楚再调用，别指望库帮你兜。

形状不匹配是最常见的问题，我在第 5 节做了个速查表，这里先记住一条判据：(m × k) * (k × n) = (m × n)，中间那个维度必须对上，输出是两个外维度。

3.3 解线性方程组：LU、Cholesky、QR 怎么选

这是绝大多数工程场景真正要解决的问题：A x = b，已知A和b，求解x。库本身不替你做选择，选错了不仅慢，还可能算不准。

最通用的是带部分主元选择的 LU 分解。它对非奇异方阵普遍适用，通过行交换来抑制数值放大：

use faer::Mat; use faer::linalg::solvers::PartialPivLu; let n = 4; let a = Mat::<f64>::from_fn(n, n, |i, j| if i == j { 4.0 } else { 1.0 }); let b = Mat::<f64>::from_fn(n, 1, |i, _| (i + 1) as f64); let lu = PartialPivLu::new(a.as_ref()); let x = lu.solve(b.as_ref());

solve接受的是视图，返回的是一个拥有数据的矩阵，形状跟b一致。这个接口的好处是可以复用分解结果：如果你有多个不同的b要解同一个A，分解只需要做一次，后面每次求解都是廉价的前代回代。这一点在参数扫描、多工况仿真里能省下大量时间——分解是 O(n³)，求解是 O(n²)，差一个数量级。

如果你的矩阵是对称正定的（来自最小二乘的正规方程、有限元刚度矩阵、协方差矩阵、核矩阵都属于这一类），务必换用 Cholesky 分解。它在数学上等价但在计算量上只要 LU 的一半，而且数值稳定性对这个矩阵类别来说是有保证的。判断标准很简单：矩阵对称（A[i][j] == A[j][i]）且所有特征值为正。实践中你不一定知道特征值，但"它是不是从某个最小二乘问题里来的"这个来源判断通常靠得住。

如果你不确定矩阵是否正定，又想要比 LU 更好的稳定性，可以走 QR 分解。它对任意形状的矩阵都能用（包括矩阵不是方阵的过定系统），代价是计算量大约是 LU 的两倍。做最小二乘拟合时，标准做法就是用 QR 解A x ≈ b，因为对正规方程直接求逆会让条件数平方，精度损失明显。

我把选择逻辑整理成一张表，按你的实际情况对号入座：

提示：判断"病态"不要靠肉眼看矩阵元素。条件数估算是更靠谱的手段，很多库提供了接口。如果你手头没有，一个粗糙但有效的替代方案是：把解回代算残差||A x - b||，残差相对量级大就说明有问题。

3.4 特征值、奇异值与最小二乘

再往上走一层就是谱分解。对称矩阵的特征值问题有专门的算法，比一般矩阵便宜很多，而且结果保证是实数：

use faer::linalg::solvers::SelfAdjointEigen; use faer::Side; let n = 5; let a = Mat::<f64>::from_fn(n, n, |i, j| { // 构造一个对称矩阵 ((i + 1) * (j + 1)) as f64 / (i + j + 2) as f64 }); let eig = SelfAdjointEigen::new(a.as_ref(), Side::Lower); let values = eig.eigen_values(); let vectors = eig.eigen_vectors();

这里Side::Lower告诉算法"只读下三角就够了，上三角是它的镜像"。这个参数的用意是省掉一半的读取和计算，你必须保证矩阵真的对称，否则算法读的是下三角，上三角的数据被无视，算出来的结果会和你的预期对不上，而且不会报错。这是我见过最隐蔽的坑之一——尤其是当你的矩阵是从某个近似对称的数值过程来的，浮点误差让A[i][j]和A[j][i]差了 1e-16，虽然影响不大，但如果你搞反了哪边是"真实数据"，问题就来了。

奇异值分解（SVD）是数值线性代数里的万能工具：求伪逆、判断秩、做降维、做数据压缩。它的计算量比特征值分解大，但适用面也更广——不要求方阵、不要求对称。

let a = Mat::<f64>::from_fn(6, 4, |i, j| ((i * 3 + j) % 5) as f64); let svd = a.svd(); let s = svd.s(); // 奇异值，按降序排列 let u = svd.u(); let v = svd.v();

奇异值默认降序排列，这个约定很好用：你想看"能量集中在哪几个方向上"，直接看前几个奇异值的衰减速度就行。如果要降维，取前 k 个奇异值对应的列构造低秩近似，这是 PCA 的核心步骤。

最小二乘我用得最多的场景是曲线拟合。给定一组观测点和一组基函数，构造设计矩阵后解一个过定系统，QR 分解是首选：

let a = Mat::<f64>::from_fn(m, n, |i, j| basis(j, x[i])); let b = Mat::<f64>::from_fn(m, 1, |i, _| y[i]); let qr = a.qr(); let coeffs = qr.solve_lstsq(b.as_ref());

我要专门提醒一个容易被忽略的问题：基函数的选择对数值稳定性的影响比求解器大得多。如果你打算拟合一个多项式，直接用1, x, x², x³...当基，当 x 的范围稍微大一点（比如 [0, 10]），设计矩阵的条件数会爆炸，哪怕你用的是 QR 分解，精度该丢还是丢。解决办法是用正交多项式基，或者先把 x 归一化到 [-1, 1]。这个坑我在做传感器标定曲线拟合时踩过一次，当时以为是求解器不行，折腾了半天换库，最后发现是基函数的问题。

4. 性能调优：它凭什么快，你怎么接得住

4.1 微内核与分块：性能的真正来源

很多人以为线性代数库快是因为"用了 SIMD"，这个理解只对了一小半。真正的性能来源是三层结构配合。

最底层是微内核，处理的是非常小的一块（比如 8×8 或 16×4 这样的寄存器级分块），它被手写成能塞进 CPU 寄存器、把乘加流水线打满的形式，中间不留任何函数调用和分支。这一层的效率决定了整个库的天花板。

中间层是分块。原因在于缓存：一个 2000×2000 的f64矩阵占 32MB，远超 L2 甚至部分 L3 缓存。如果你按最朴素的三重循环算C = A * B，每个元素都要去主存里取一次数据，缓存疯狂失效，实测性能可能只有理论峰值的十分之一。分块的做法是把大矩阵切成能装进缓存的块，在块内密集计算，让数据在被换出缓存之前尽可能多地被复用。块大小的选择直接决定了缓存命中率，这是一门要靠实测调参的手艺活。

最上层是并行。大矩阵乘法天然可并行——把输出矩阵切块，不同线程算不同的块，互不干扰。这里的调度策略很重要：任务粒度太细，线程间同步开销吃掉收益；粒度太粗，负载不均衡，最后一个线程拖后腿。好的实现会根据矩阵规模动态选择粒度。

理解了这三层，你就能理解为什么"矩阵太小的时候反而要关并行"：小矩阵连微内核那层的效率都填不满，哪还有空间给分块和并行。

4.2 并行阈值的调校思路

库通常有一个默认的并行阈值——小于这个规模就用单线程。这个默认值是根据通用场景调的，你的硬件和工作负载未必匹配。

如果你发现某段代码在中等规模矩阵上表现异常，可以先做个简单实验：固定矩阵尺寸，分别用Par::Seq和Par::Rayon(0)跑，看交叉点在哪。我一般写个小循环扫一遍尺寸：

use faer::Par; // 伪代码示意：扫描不同尺寸，比较单线程与多线程 for &sz in &[32usize, 64, 128, 256, 512, 1024, 2048] { let lhs = Mat::<f64>::from_fn(sz, sz, |i, j| ((i ^ j) % 11) as f64); let rhs = Mat::<f64>::from_fn(sz, sz, |i, j| ((i + j) % 13) as f64); let mut dst = Mat::<f64>::zeros(sz, sz); // 单线程计时 let t0 = std::time::Instant::now(); for _ in 0..repeats { matmul(dst.as_mut(), Accum::Replace, lhs.as_ref(), rhs.as_ref(), 1.0, Par::Seq); } let seq = t0.elapsed(); // 多线程计时 let t0 = std::time::Instant::now(); for _ in 0..repeats { matmul(dst.as_mut(), Accum::Replace, lhs.as_ref(), rhs.as_ref(), 1.0, Par::Rayon(0)); } let par = t0.elapsed(); println!("size={} seq={:?} par={:?}", sz, seq, par); }

这种扫描跑一遍就能看出你的机器上交叉点大概在哪一档。我的经验是，交叉点通常比多数人猜的要低一些——现代多核机器上，128 到 256 阶之间就开始有明显收益了。但这也跟你的矩阵内容有关：如果内存带宽是瓶颈，加线程带来的收益会被摊薄。

注意：做这类对比实验时，务必在 release 模式下跑。debug 模式下的性能数据毫无参考价值，浮点运算在未优化时会慢一到两个数量级。

4.3 减少分配：比调参数更实在的优化

调并行阈值是微调，减少不必要的内存分配是数量级的改进。我总结了三条在项目里反复用的规则。

第一条：分解结果能复用就复用。前面说过，同一矩阵多个右端项的场景，分解只做一次。更进一步，如果你的算法在一个迭代循环里反复求解稍微变化的矩阵，可以考虑是否能用增量更新或者预条件技术，而不是每次都重新分解。

第二条：避免在热路径上构造临时矩阵。let temp = &a * 2.0;这种写法每次都会分配一块新内存。如果它在循环里，每次迭代都触发一次malloc加一次memcpy加一次free，开销能吃掉可观的比例。解决办法是提前分配好缓冲区，然后调用所有参数都是视图的接口，把结果写进已有缓冲区。matmul就是这种设计，dst是传可变视图进去的，天然支持复用。

第三条：注意隐藏的拷贝。有些操作看起来没拷贝，实际上触发了。典型的是对一个视图调用需要拥有权的方法。养成习惯：凡是能用as_ref()传视图的地方就别传值；必要时显式标注参数类型，让编译器帮你把"本来想传视图结果传了拥有型"的地方暴露出来。

我在一个图像处理的小工具里做过对比：把三处热循环里的临时矩阵分配改掉、改成复用缓冲区，端到端时间从 4 秒多降到 1 秒出头。算法一行没改，纯粹是内存行为的差别。这个例子说明一个道理：先把内存行为理顺，再去纠结 SIMD 和线程数。

4.4 基准测试该怎么做才可信

最后说方法论，因为网上流传的性能数字大多不可信。

用工具而不是Instant::now()做正式测量。手写计时的抖动太大，尤其在测量微秒级操作时，一次意外的上下文切换就能让结果翻倍。用基准测试框架（criterion是 Rust 生态里的常规选择），它会自动做预热、多次采样、统计置信区间，还能检测结果是否统计显著。一次基准大概长这样：

use criterion::{criterion_group, criterion_main, BenchmarkId, Criterion}; fn bench_solve(c: &mut Criterion) { let mut group = c.benchmark_group("solve"); for &n in &[64usize, 128, 256, 512] { let a = Mat::<f64>::from_fn(n, n, |i, j| if i == j { 4.0 } else { 1.0 }); let b = Mat::<f64>::from_fn(n, 1, |i, _| 1.0); group.bench_with_input(BenchmarkId::from_parameter(n), &n, |bench, _| { bench.iter(|| { let lu = PartialPivLu::new(a.as_ref()); lu.solve(b.as_ref()) }); }); } group.finish(); } criterion_group!(benches, bench_solve); criterion_main!(benches);

固定输入数据，别用随机数。随机数据每次不同，一方面引入噪声，另一方面某些算法对输入值的分布敏感（比如全零矩阵会走捷径分支），用确定性生成的伪随机数更稳妥。上面例子里的(i ^ j) % 11就是一种确定性"看起来乱但不真乱"的填充。

测你要用的尺寸区间。别拿 4096 阶的漂亮数字去支持你那个 30 阶的实际场景，两者优化策略完全不同。

锁定线程数。多线程基准跑出来的数字，换个核数不同的机器就不可比。做对比实验时用RAYON_NUM_THREADS之类的环境变量固定住，结果才有意义。

5. 常见坑与排查速查表

5.1 数值层面的坑

对称矩阵只传了一半三角。前面提过，用Side::Lower时算法只读下三角。如果你的矩阵实际上不对称，或者数据在另一侧，结果就是错的，而且不会有任何警告。防御手段：在提交给求解器之前，写个断言检查(a[(i,j)] - a[(j,i)]).abs() < eps，在调试构建里开着，release 里关掉。

正定性没验证就用 Cholesky。Cholesky 在遇到非正定矩阵时可能返回错误，也可能因为浮点误差"勉强完成"并给出一个毫无意义的结果。判断正定最廉价的方法是尝试分解并检查返回值，别自己用顺序主子式去判断——那个判据在数值上极不稳定。

条件数被忽略。解出来的x看起来正常，回代残差却很大，说明问题本身病态。我养成的一个习惯是：任何关键求解之后都算一次相对残差||A x - b|| / ||b||，超过 1e-8 就开始警觉。这不是库的问题，是数学问题，换任何库都一样。

浮点累加的误差累积。如果你的算法是长链累加（比如迭代求解、或者把几万个小矩阵乘起来），误差会累积。可以考虑用 Kahan 求和、或者中间结果提升到更高精度，甚至用双精度算单精度的输入。这类问题在图形学的变换链里尤其常见。

5.2 编译与运行时层面的坑

性能莫名其妙只有预期的一小部分。第一件事检查是不是没配target-cpu。第二件事检查是不是用 debug 模式在跑。这两个原因占了我遇到过的性能问题的绝大多数。

分发的二进制在别人机器上崩溃。典型的target-cpu=native后遗症。改成明确指定的微架构等级重新编译。

链接时间突然变长。开了lto = "fat"和codegen-units = 1的代价。建议只对最终产物开，或者干脆只在 CI 的发布流水线里开。

多线程环境下结果偶尔不一致。浮点加法的结合律不成立，多线程归约的顺序不固定，所以结果在最低几位上可能不同。这不是 bug，是浮点的本性。如果你的下游逻辑对最后几位敏感（比如做相等比较来判断收敛），把收敛判据改成带容差的比较。

5.3 排查速查表

6. 落地实践：把 faer 包进真实项目

6.1 封装一层，别让调用方直接见库类型

库在快速迭代期，API 会变。如果你的代码里Mat、MatRef、各种分解器类型散布在几十个文件里，某次升级你就得改几十处。我的做法是在项目里建一个薄薄的封装层，对外只暴露自己的类型和函数签名。

封装层的价值不只是隔离版本变化，还有几个实际好处。它可以统一你的数值类型选择——比如整个项目固定用f64，那封装层里就不需要泛型，接口更简单。它可以集中做输入校验——形状检查、正定性检查、NaN 检查，只写一遍。它还可以统一错误处理——把库的 panic 转成你自己的Result，避免一个数值错误直接掀翻整个服务。

具体形态大概是：一个Solver结构体持有分解结果，一个solve方法接受切片形式的输入输出，内部负责视图转换。这样外部代码面对的是普通的&[f64]，跟矩阵库完全解耦。

一个额外的好处是测试好写。你可以在封装层外面挂一组已知答案的测试用例（比如用对角矩阵验证、用单位阵验证、用解析解验证），升级库版本时跑一遍就知道有没有回归。我建议至少准备三类用例：平凡解（单位阵、对角阵）、有解析解的小规模问题、以及一个规模稍大的真实数据快照。

6.2 大矩阵的内存规划

内存是这类库最容易出问题的地方，因为矩阵的尺寸是平方增长的。一个 10000×10000 的f64矩阵是 800MB，而且很多算法内部还要额外分配工作区。你在 8GB 内存的开发机上可能跑得动，在容器里就 OOM 了。

几条实操建议。第一，估算峰值内存再动手。一个 LU 分解大致需要n²的存储加上n的置换向量，Cholesky 大约n²/2（只存一半），QR 可能要到n²甚至更多。SVD 的开销比 LU 大不少。心里有个数，别在临界点上试探。

第二，能分批就分批。如果你的最终目标是"处理 N 个样本，每个样本做一个中等规模求解"，那完全没有必要把 N 个矩阵全装上一次性算，一个循环里逐个处理、复用缓冲区，峰值内存就是单个样本的量级。

第三，注意工作区分配。有些分解接口允许传入预分配的工作区，能省掉每次调用的分配开销。在迭代场景里这个收益不小。

第四，小心栈溢出。某些实现会用栈上的定长数组做微内核缓冲区，如果你把矩阵尺寸调到很极端，在某些平台上可能触发栈问题。遇到莫名其妙的段错误（在 unsafe 或者 FFI 边界附近），先怀疑栈。

6.3 稀疏矩阵与结构化数据的取舍

最后聊聊稀疏。当你的矩阵里非零元素占比很低（一般低于 5% 到 10% 才值得考虑），稠密存储就是纯浪费——你花内存和算力去处理一堆零。

稀疏矩阵的存储一般是三元组或者压缩列格式，只存非零元素及其位置。好处是内存从 O(n²) 降到 O(nnz)，矩阵向量乘法从 O(n²) 降到 O(nnz)。代价是随机访问变慢，而且算法复杂度大幅上升——稀疏 LU 分解需要重排序来减少填充（fill-in），填充控制得好不好直接决定了性能和内存。

faer 的稀疏支持在持续演进，覆盖范围不如稠密部分那么全面。我的建议是：如果你的问题规模在稠密方案能接受的范围内（比如一万阶以下、内存够），优先用稠密，简单可靠；只有当内存或时间明显不可接受时，再考虑稀疏，并且做好前期调研——把稀疏求解当默认选项的人，往往低估了它的复杂度。

还有一个中间地带值得知道：结构化矩阵。如果你的矩阵有特殊结构（对称、带状、Toeplitz、低秩加对角），很多问题可以不经过通用分解，用专门算法直接解，速度快几个数量级。这需要你在建模阶段就有意识地识别结构，而不是等矩阵建好了才想优化。

我个人在几个项目里用下来的体会是，这个库最舒服的地方不是某个具体数字，而是那种"不用为构建环境操心"的踏实感——cargo build完就是能跑的产物，交叉编译到哪个平台都一个样，CI 里不用为了一个数学库装一堆额外工具链。性能方面它已经足够好到让我不再想去找替代方案，真正需要盯的反而是我自己的代码：有没有在热路径上乱分配、有没有把不该分摊的分解重复做了、有没有在基函数和条件数上埋雷。最后分享一个我觉得最值的小习惯——任何一次关键求解之后都算一遍残差并打日志，这个成本可以忽略，但它救回来的排查时间是以天计的。