矩阵分解与数值线性代数：LU、Cholesky、QR、特征分解与 SVD¶

一篇会成长的笔记。卡尔曼/组合导航里，"解方程"和"保持数值稳定"无处不在：求增益要解线性方程、协方差要平方根化、可观测性要数奇异值、姿态确定要用 SVD。这篇把最常用的矩阵分解一次讲透：LU / Cholesky / QR / 特征分解 / SVD，按"定义 → 推导 → 几何直觉 → 课程里用在哪 → 常见坑"。每节末尾附外链 Wiki。

核心要点¶

• 矩阵分解 = 把 \(A\) 拆成"好算的因子连乘"，一次性换来三样东西：解方程快、数值稳、结构可读。
• LU：高斯消元的矩阵写法，解 \(Ax=b\) 的最朴素工具（记得配部分主元）；
• Cholesky：对称正定矩阵的特权 \(A=LL^T\)——卡尔曼 \(P\) 阵的平方根形式（平方根滤波）、高斯采样的发动机；
• QR：最小二乘的数值稳定解法（\(A=QR\)），避免正规方程 \(A^TA\) 的条件数平方爆炸；
• 特征分解（对称阵谱分解 \(A=Q\Lambda Q^T\)）：误差椭球主轴、\(e^{At}\) 的优雅表达；
• SVD：对任意矩阵都成立、最稳的分解 \(A=U\Sigma V^T\)——可观测性分析、伪逆、Wahba 姿态确定全用它；
• 选型口诀：对称正定 → Cholesky；最小二乘 → QR；要特征/奇异值信息 → 特征分解/SVD；只解方程 → LU。

一、为什么需要矩阵分解¶

一句话：解线性方程 \(Ax=b\) 有三种境界——直接高斯消元（一步步磨）、求逆 \(A^{-1}b\)（贵且不稳）、把 \(A\) 拆成好算的因子（一次拆、多次用、数值稳）。

分解带来的三个直接好处：

1. 快：\(PA=LU\) 拆一次，\(Ax=b\) 退化成两次三角回代，复杂度 \(O(n^2)\) 而不是 \(O(n^3)\)；卡尔曼增益每步都要解这种方程。
2. 稳：直接求逆会放大浮点误差；分解后乘的顺序可控、误差不倍增（数值线性代数的全部意义）。
3. 读结构：特征值/奇异值直接告诉你怎么缩、怎么转、哪些方向"看不见"（可观测性）——这是组合导航设计要的东西。

二、LU 分解（高斯消元的矩阵版本）¶

2.1 定义¶

其中 \(L\) 是单位下三角（对角线全 1），\(U\) 是上三角。

例（\(2\times2\)）：

2.2 怎么用（解两次三角方程替代一次高斯消元）¶

求解 \(Ax=b\) 分两步：

两次都是三角回代，各 \(O(n^2)\)，比每步重新消元快得多。

2.3 数值注意：部分主元¶

当主元很小时，直接消元会放大误差。实际代码用

\(P\) 是置换矩阵（把当前列最大元素换到主元位置）——部分主元让算法对病态系统稳得多（对应深坑 §八.1）。

直觉：LU 就是"把高斯消元的每一行操作记下来，变成矩阵乘"，这样同一个 \(A\) 可以反复配不同的 \(b\)——卡尔曼里 \(F,H\) 每步变但结构固定，就适合这种"拆一次、用多次"。

三、Cholesky 分解（对称正定的特权）¶

3.1 定义与存在性¶

若 \(A\) 是对称正定矩阵（卡尔曼的 \(P,R,Q\) 都是），则存在唯一的下三角 \(L\)（对角线>0），使得

存在且唯一的充要条件就是对称正定——这是其他分解没有的"专属特权"。

3.2 计算（逐列递推）¶

对 \(i=1,\dots,n\)：

本质上是"把 \(A=\sum_k L_{*k}L_{*k}^T\) 逐项剥出来"，复杂度 \(O(n^3/3)\)（比 LU 少一半）。

3.3 几何直觉：它是"标准差"的矩阵版¶

• 方差的正平方根 \(\sigma=\sqrt{\sigma^2}\) 是"不确定性的长度"；
• Cholesky \(L\) 是"协方差的平方根"：\(P=LL^T\) 把不确定性"开方"成一个三角因子。

3.4 课程里用在哪（两大用途）¶

用途 1：平方根滤波（数值稳定地保持 \(P\) 正定）

直接传播 \(P^+=FPF^T+Q\) 长时间会因浮点误差失去对称正定性；改为存 \(S\)（\(P=SS^T\)）并更新 \(S\)，就保证 \(P\) 永远半正定——这就是 Joseph 形式 / 平方根滤波（UD/Cholesky 分解）的动机（概率统计基础 §九坑 5 的解法）。

用途 2：从高斯分布采样（Monte-Carlo 仿真 / 双轨验证）

要生成 \(x\sim\mathcal N(\mu,\,P)\)：先采样 \(z\sim\mathcal N(0,I)\)，再

因为 \(\mathrm{Cov}(Lz)=L\,\mathrm{Cov}(z)L^T=LL^T=P\)。你的知识库里 MATLAB/Python 双轨对拍（PSINS 演示）生成随机轨迹就是这条。

四、QR 分解（最小二乘的稳定解法）¶

4.1 定义¶

其中 \(Q\) 是正交矩阵（列正交归一，\(Q^TQ=I\)），\(R\) 是上三角。对任意（含矩形）矩阵都存在。

4.2 为什么最小二乘要用 QR¶

超定最小二乘 \(Ax\approx b\) 的教科书解法是正规方程：

但它有个致命伤：\(A^TA\) 的条件数是 \(A\) 的条件数平方，病态时数值上直接废掉。

改用 QR：

条件数保持不变，稳得多。

4.3 怎么算¶

• Gram-Schmidt 正交化（教科书，数值不稳）；
• Householder 反射（实际首选：一次把一个子列变成标准基）；
• Givens 旋转（逐步旋转清零，适合稀疏/并行）。

直觉：QR 把 \(A\) 的列"掰成"一组正交基 \(Q\)，剩下的系数全塞进 \(R\)——相当于"换到标准正交系里再看方程"，代价最低又最稳。

五、特征分解与谱分解（椭球主轴与 \(e^{At}\)）¶

5.1 一般矩阵：\(A=V\Lambda V^{-1}\)¶

5.2 对称矩阵：谱分解（卡尔曼场景基本都是对称的）¶

\(A\) 对称 ⟹ 可正交对角化：

特征值 \(\lambda_i\) 全实数。

5.3 几何直觉¶

对称矩阵的变换 = "先转到一个特殊坐标系（\(Q^T\)），各轴独立拉伸 \(\lambda_i\)，再转回来（\(Q\)）"。所以：

• 特征向量 = 椭球主轴方向，特征值 = 各轴方差——协方差矩阵的误差椭圆就是这么看出来的（数学基础 §二）；
• \(e^{At}\) 有了优雅表达：\(e^{At}=Qe^{\Lambda t}Q^T\)，对角指数逐元素算——接上 微分方程与状态空间离散化 §三的矩阵指数。

直觉：谱分解把"旋转变换"翻译成"沿主轴的中点拉伸"；\(e^{At}\) 的每个模式就是一根特征值指数。

六、SVD 奇异值分解（万能且最稳）¶

6.1 定义¶

对任意 \(m\times n\) 矩阵 \(A\)：

• \(U\)（\(m\times m\) 正交）、\(V\)（\(n\times n\) 正交）——列分别是 \(AA^T\)、\(A^TA\) 的特征向量；
• \(\Sigma\)：\(\sigma_1\ge\sigma_2\ge\cdots\ge\sigma_r\ge0\)（奇异值，\(r=\mathrm{rank}(A)\)），其余为 0。

6.2 与特征分解的区别（别混）¶

6.3 几何直觉¶

奇异值 \(\sigma_i\) = 各方向的压缩倍数：大的方向"信息强"、小的方向"几乎压没了"。

6.4 课程里用在哪（三大战场）¶

① 伪逆与最小二乘：\(A^\dagger = V\Sigma^\dagger U^T\)（奇异值取倒数），解 \(Ax\approx b\) 的最小范数解，比正规方程稳。

② 可观测性分析（组合导航的灵魂问题）：可观测性矩阵 \(Ob=[H;\ HF;\ HF^2;\dots]\) 做 SVD：

• 零奇异值 = 该方向状态从量测中看不见（比如静基座对方位失准角不可观——经典结论就是这么来的）；
• 奇异值大小 = 可观测方向上的"可见强度"；
• 这是 09 纯惯导误差传播 / 12 观测模型 里"哪些误差能估、哪些不能"的数值判据。

③ Wahba 姿态确定（初始对准/双矢量定姿）：求使 \(\sum |b_i-Rr_i|^2\) 最小的旋转 \(R\)，解就是 \(B=\sum b_ir_i^T\) 的 SVD 结构：\(R=U\,\mathrm{diag}(1,1,\det(UV^T))V^T\)——08 初始对准 里 TRIAD 的"更好版本"。

④ 条件数诊断：\(\mathrm{cond}(A)=\sigma_{\max}/\sigma_{\min}\)——矩阵病态程度，数值滤波/可观测性分析里的体检指标。

七、选型总表（什么时候用哪个）¶

八、常见坑¶

1. 忘配部分主元直接 LU：主元为 0 或极小就崩；代码里永远用 \(PA=LU\)。
2. 对非对称正定的矩阵用 Cholesky：出现负的平方根就是上当了——先用 \(\lambda_{\min}\) 检查条件。
3. 正规方程当默认解最小二乘：condition 平方爆炸，病态数据用 QR/SVD。
4. 混淆特征分解与 SVD：特征分解只对"可对角化的方阵"，SVD 对任意矩阵且更稳；不要把非对称阵硬塞特征分解。
5. 直接传播 \(P\) 导致失正定：长跑用 \(P=SS^T\)（Cholesky/UD）或 Joseph 形式，别等 \(P\) 变负了才后悔（对应 概率统计基础 §九坑 5）。
6. 忽略小奇异值：可观测性分析里"小但不为零"的奇异值不是噪声就是弱可观，要先归一化再定阈值，别一眼看成"不可观"。
7. 条件数凭空估：让计算机算 \(\sigma_{\max}/\sigma_{\min}\)，别用行列式判断病态（行列式会骗人）。

九、一句话总结¶

矩阵分解是数值线性代数的"乐高"：LU 拆方程、Cholesky 开平方、QR 保稳定、特征分解看主轴、SVD 通吃一切——卡尔曼/组合导航的稳定性、可观测性、姿态确定，全站在这些分解上。

十、数值案例：手算验证（每个分解一个可验算的例子）¶

下面每个例子都小到能笔算，目的是让你"看见"分解到底在算什么、结果对不对。建议拿纸笔跟一遍。

10.1 LU（连带解一次方程）¶

取

第一步消元：第二行减第一行的 2 倍 → 消元矩阵 \(L_1=\begin{bmatrix}1&0\\-2&1\end{bmatrix}\)，得 \(L_1A=\begin{bmatrix}2&1\\0&1\end{bmatrix}=U\)。于是

其中 \(L=\begin{bmatrix}1&0\\2&1\end{bmatrix}\)（单位下三角）、\(U=\begin{bmatrix}2&1\\0&1\end{bmatrix}\)（上三角）。验算 \(LU=A\) 成立。

用 LU 解 \(Ax=b\)，取 \(b=\begin{bmatrix}5\\13\end{bmatrix}\)：

1. 前向 \(Ly=b\)：\(y_1=5,\ 2y_1+y_2=13\Rightarrow y_2=3\)，得 \(y=\begin{bmatrix}5\\3\end{bmatrix}\)；
2. 后向 \(Ux=y\)：\(x_2=3,\ 2x_1+x_2=5\Rightarrow x_1=1\)，得 \(x=\begin{bmatrix}1\\3\end{bmatrix}\)。

验算 \(A\begin{bmatrix}1\\3\end{bmatrix}=\begin{bmatrix}5\\13\end{bmatrix}=b\)（成立）。注意解里没有出现过 \(A^{-1}\)——这就是分解的意义。

10.2 Cholesky（协方差开平方）¶

取对称正定

逐列递推：

得

物理含义：\(L\) 就是"协方差 \(P\) 的平方根"——要生成 \(x\sim\mathcal N(0,P)\)，算 \(x=Lz,\ z\sim\mathcal N(0,I)\) 即可（§三.4 用途 2）。

10.3 QR（把列向量掰正交）¶

取

列 \(a_1=\begin{bmatrix}1\\1\\0\end{bmatrix},\ a_2=\begin{bmatrix}1\\0\\1\end{bmatrix}\)。Gram–Schmidt：

于是 \(Q=[q_1\ q_2]\)，\(R=\begin{bmatrix}\sqrt2&1/\sqrt2\\0&\sqrt{1.5}\end{bmatrix}\)，验算 \(QR=A\)（成立）。可见 \(R\) 把"原始列在正交基下的坐标"记下来了。

10.4 SVD 与条件数（为什么不能直接求逆）¶

取

它本身已是对角 SVD：\(U=V=I,\ \Sigma=\begin{bmatrix}1&0\\0&0.1\end{bmatrix}\)，条件数 \(\kappa=\sigma_{\max}/\sigma_{\min}=10\)。

若换成近奇异的

则 \(\kappa=10^{10}\)，而 \(A^{-1}=\begin{bmatrix}1&0\\0&10^{10}\end{bmatrix}\)——元素直接爆到 \(10^{10}\)。这就是"直接求逆 → 放大误差 → 滤波器发散"的数值根源：条件数一大，求逆把微小浮点误差放大成天文数字。SVD 的 \(\sigma_{\min}\) 一眼就告诉你矩阵有多病态。

十一、物理意义速查：每个分解在 ESKF 里对应什么¶

矩阵分解不是数学炫技，是保命工具：直接求逆在浮点误差下极易让 \(P\) 失去正定、滤波器瞬间发散（"飞点"/NaN）。下表把五个分解一一对应到 ESKF 15 态里你真正在乎的物理量。

一句话：对称正定（\(P,R,Q\)）→ Cholesky 开方；最小二乘 / 求增益 → QR；要"哪些能估哪些不能" → SVD；只解普通方程 → LU。

十二、代码对比：MATLAB ↔ NumPy / SciPy¶

12.1 函数对照表¶

注意 NumPy svd 的坑：MATLAB 的 svd 返回 \(V\)，NumPy 返回 \(V^H\)（共轭转置）。所以 MATLAB 里 \(A=U S V^T\)，NumPy 里 \(A=U\,\mathrm{diag}(S)\,V_h\)。做 Wahba / 伪逆时别把 \(V_h\) 当 \(V\) 直接用，要先 .T。

12.2 求卡尔曼增益：三种写法（MATLAB + NumPy）¶

设 \(S=H P H^T+R\)，增益 \(K=P H^T S^{-1}\)。前两种都不直接求 \(S^{-1}\)。

写法 A：直接求逆（不稳，仅作对照）

K = P * H' * inv(S);              % MATLAB

K = P @ H.T @ np.linalg.inv(S)    # NumPy

写法 B：Cholesky 解（对称正定首选）

% MATLAB：解 S*K^T = (P*H')^T
L = chol(S, 'lower');
K = (L' \ (L \ (P*H')'))';

# NumPy
L = np.linalg.cholesky(S)
K = np.linalg.solve(L.T, np.linalg.solve(L, (P @ H.T).T)).T

写法 C：QR 解（S 非对称 / 奇异时也稳）

% MATLAB：解 R*K^T = Q'*(P*H')^T
[Qq, Rq] = qr(S);
K = (Rq \ (Qq' * (P*H')'))';

# NumPy
Qq, Rq = np.linalg.qr(S)
K = np.linalg.solve(Rq, Qq.T @ (P @ H.T).T).T

三种写法在 \(S\) 良态时结果一致；\(S\) 条件数一大，写法 A 的数值误差会明显大于 B / C——这正是组合导航里坚持用 B / C 的原因（对应 §八 坑 3、§十一 物理意义）。

关联¶

• 系列内：惯性导航系列首页
• 姊妹篇：数学地基（矩阵/协方差/特征值入口）· 概率统计基础（\(P\) 正定性与平方根化）· 微分方程与状态空间离散化（\(e^{At}\) 与状态转移）
• 应用篇：08 初始对准（Wahba/TRIAD）· 09 纯惯导误差传播（可观测性）· 12 观测模型 · 10 卡尔曼滤波基础（\(K\) 涉及求逆/条件数）
• 想深入数值：MATLAB chol/qr/lu/svd 文档 + 数值线性代数 — 维基

外链 Wiki¶

• 矩阵分解 — 维基百科
• Cholesky 分解 — 维基百科
• QR 分解 — 维基百科
• 特征分解 — 维基百科
• 奇异值分解 — 维基百科
• LU 分解 — 维基百科
• 条件数 — 维基百科
• Wahba 问题 — Wikipedia