矩阵分解 — LU, Cholesky, QR, SVD¶

算法流程¶

LU 分解（Doolittle 算法）¶

目标：将方阵 \(A\) 分解为 \(A = P \cdot L \cdot U\)，其中 \(L\) 是单位下三角，\(U\) 是上三角，\(P\) 是置换矩阵。

算法流程（带部分主元）：

1. 复制 \(A\) 到工作矩阵 \(U\)，初始化 \(L = I\), \(P = I\)

2. 对每一列 \(k = 0, \dots, n-2\)：

3. 选主元：在 \(U[k..n-1, k]\) 中找绝对值最大的元素，记录行号 \(p\)

4. 换行：交换 \(U[k, :]\) 与 \(U[p, :]\)，\(P[k, :]\) 与 \(P[p, :]\)，\(L[k, 0..k-1]\) 与 \(L[p, 0..k-1]\)

5. 消元：对每一行 \(i = k+1, \dots, n-1\)：

   • \(L_{ik} = U_{ik} / U_{kk}\)
   • \(U_{i, k..n-1} \;\;{-}{=}\;\; L_{ik} \cdot U_{k, k..n-1}\)
   • 返回 \(L, U, P\)

复杂度：\(O(n^3)\) 次浮点运算。

选型建议：

• 对一般方阵，始终启用主元（use_pivoting = true）
• 若已知矩阵对角占优，可省略主元以节省约 20% 时间

Cholesky 分解¶

目标：将对称正定矩阵 \(A\) 分解为 \(A = L \cdot L^T\)，其中 \(L\) 是下三角且对角元素 \(> 0\)。

算法流程（原地计算版）：

对每一列 \(j = 0, \dots, n-1\)：

1. 计算对角元素：\(L_{jj} = \sqrt{A_{jj} - \sum_{k=0}^{j-1} L_{jk}^2}\)

2. 对每一行 \(i = j+1, \dots, n-1\)：

3. \(L_{ij} = \left( A_{ij} - \sum_{k=0}^{j-1} L_{ik} L_{jk} \right) / L_{jj}\)

关键特性：

• 速度：约 LU 的一半运算量（利用对称性）
• 稳定性：SPD 矩阵上数值稳定，无需选主元
• 存储：只需存储 \(L\)（约 LU 的一半内存）
• 前置条件：必须先确认 \(A\) 对称正定（调用 is_symmetric() + is_positive_definite()）

复杂度：\(O(n^3 / 3)\) 次浮点运算。

QR 分解（Householder 反射）¶

目标：将任意 \(m \times n\) 矩阵 \(A\) 分解为 \(A = Q \cdot R\)，其中 \(Q\) 正交，\(R\) 上三角。

算法流程：

对每一列 \(k = 0, \dots, n-1\)（假设 \(m \geq n\)）：

1. 取当前子列 \(\mathbf{x} = A[k..m-1, k]\)

2. 构造 Householder 向量 \(\mathbf{v}\)：

3. \(\alpha = -\text{sign}(x_0) \cdot \|\mathbf{x}\|_2\)

4. \(\mathbf{v} = \mathbf{x} - \alpha \cdot \mathbf{e}_0\)，\(v_0 = x_0 - \alpha\)
5. \(\beta = 2 / \|\mathbf{v}\|_2^2\)
6. 应用反射：\(A[k..m-1, k..n-1] \;\;{=}\;\; A - \beta \cdot \mathbf{v} \cdot (\mathbf{v}^T \cdot A)\)
7. 累积 \(Q\)：\(Q = Q \cdot (I - \beta \cdot \mathbf{v} \cdot \mathbf{v}^T)\)

对比 MGS（用于 C 库 tiny_decomp）：

复杂度：\(O(2mn^2)\) 次浮点运算。

SVD 分解（双边 Jacobi + QR 精化）¶

目标：将任意 \(m \times n\) 矩阵 \(A\) 分解为 \(A = U \cdot \Sigma \cdot V^T\)，其中 \(U, V\) 正交，\(\Sigma\) 对角（奇异值 \(\sigma_1 \geq \dots \geq \sigma_r \geq 0\)）。

算法流程：

Phase 1 — 构造对称矩阵：

1. 计算 Gram 矩阵 \(G = A^T \cdot A\)（大小 \(n \times n\)）
2. 问题转化为求 \(G\) 的特征值与特征向量

Phase 2 — Jacobi 旋转：

1. 初始化 \(V = I_n\)，取工作矩阵 \(B = G\)

2. 扫描所有非对角对 \((p, q)\)，其中 \(p < q\)：

3. 若 \(|B_{pq}|\) 大于容差，计算 Givens 旋转角 \(\theta\)：

   • \(\tau = (B_{qq} - B_{pp}) / (2 \cdot B_{pq})\)
   • \(t = \text{sign}(\tau) / (|\tau| + \sqrt{1 + \tau^2})\)
   • \(c = 1 / \sqrt{1 + t^2}\), \(s = t \cdot c\)
   • 应用旋转：\(B = J^T \cdot B \cdot J\)，\(V = V \cdot J\)
   • 重复扫描直到 \(\max |B_{pq}| < \text{tolerance}\)

Phase 3 — 排序与截断：

1. 将特征值（= 奇异值平方）降序排列，同步排列 \(V\) 的列
2. 取平方根得到奇异值：\(\sigma_i = \sqrt{|\lambda_i|}\)
3. 恢复 \(U\)：\(U = A \cdot V \cdot \Sigma^{-1}\)
4. 根据容差确定数值秩：\(\text{rank} = \#\{i \mid \sigma_i > \text{tolerance}\}\)

复杂度：约 \(O(n^3)\)，其中 \(n\) 是列数。

坑点：

• 通过 \(A^T A\) 间接计算会损失精度（条件数平方倍放大）；备选方案是双向 Jacobi 直接作用于 \(A\)
• 容差选择直接影响秩和伪逆结果

矩阵属性与分解¶

矩阵属性检查¶

检查对称性¶

bool Mat::is_symmetric(float tolerance = 1e-6f) const;

描述:

检查矩阵在给定容差内是否对称。矩阵A是对称的，如果A = A^T，即对于所有i, j，A(i,j) = A(j,i)。

数学原理:

对于对称矩阵，所有特征值都是实数，并且可以选择特征向量使其正交。对称矩阵在许多应用中都是基础的，特别是在结构动力学和优化中。

参数:

• float tolerance : 允许的最大差值 |A(i,j) - A(j,i)|（默认：1e-6）。

返回值:

bool - 如果近似对称返回true，否则返回false。

使用建议:

• 特征分解: 对称矩阵可以使用更高效和稳定的特征分解方法（例如Jacobi方法）。

• Cholesky分解: 只有对称正定矩阵可以使用Cholesky分解进行分解。

• 结构动力学: 结构分析中的刚度和质量矩阵通常是对称的。

检查正定性¶

bool Mat::is_positive_definite(float tolerance = 1e-6f, int max_minors_to_check = -1) const;

描述:

使用Sylvester准则检查矩阵是否正定。对称矩阵A是正定的，如果对于所有非零向量x，x^T A x > 0，或者等价地，所有特征值都是正的。

数学原理:

Sylvester准则指出，对称矩阵是正定的当且仅当所有前导主子式都是正的。为了效率，函数检查前几个前导主子式和对角元素。

参数:

• float tolerance : 数值检查的容差（默认：1e-6）。
• int max_minors_to_check : 要检查的前导主子式个数。-1表示检查全部；>0表示仅检查前若干个。

返回值:

bool - 如果矩阵是正定的返回true，否则返回false。

使用建议:

• Cholesky分解: 正定矩阵可以使用Cholesky分解进行分解，这比LU分解更快、更稳定。

• 优化: 正定Hessian矩阵在优化问题中表示局部最小值。

• 稳定性分析: 在控制系统中，某些矩阵的正定性确保系统稳定性。

矩阵分解结构¶

LU分解结构¶

struct Mat::LUDecomposition
{
    Mat L;                 // 下三角矩阵（单位对角线）
    Mat U;                 // 上三角矩阵
    Mat P;                 // 置换矩阵（如果使用主元）
    bool pivoted;          // 是否使用主元
    tiny_error_t status;   // 计算状态

    LUDecomposition();
};

描述:

LU分解结果的容器。分解A = P * L * U（带主元）或A = L * U（不带主元），其中L是单位对角线的下三角矩阵，U是上三角矩阵，P是置换矩阵。

数学原理:

LU分解将矩阵分解为下三角和上三角矩阵，实现线性系统的高效求解。使用主元时，可以更好地处理近奇异矩阵。

Cholesky分解结构¶

struct Mat::CholeskyDecomposition
{
    Mat L;                 // 下三角矩阵
    tiny_error_t status;   // 计算状态

    CholeskyDecomposition();
};

描述:

Cholesky分解结果的容器。对于对称正定矩阵，A = L * L^T，其中L是下三角矩阵。

数学原理:

Cholesky分解是用于对称正定矩阵的专用LU分解。它只需要LU分解的一半存储和计算。

QR分解结构¶

struct Mat::QRDecomposition
{
    Mat Q;                 // 正交矩阵 (Q^T * Q = I)
    Mat R;                 // 上三角矩阵
    tiny_error_t status;   // 计算状态

    QRDecomposition();
};

描述:

QR分解结果的容器。A = Q * R，其中Q是正交的（Q^T * Q = I），R是上三角矩阵。

数学原理:

QR分解将矩阵表示为正交矩阵和上三角矩阵的乘积。它在数值上稳定，是最小二乘问题的基础。

SVD分解结构¶

struct Mat::SVDDecomposition
{
    Mat U;                 // 左奇异向量（正交矩阵）
    Mat S;                 // 奇异值（对角矩阵或向量）
    Mat V;                 // 右奇异向量（正交矩阵，V^T）
    int rank;              // 矩阵的数值秩
    int iterations;        // 执行的迭代次数
    tiny_error_t status;   // 计算状态

    SVDDecomposition();
};

描述:

SVD分解结果的容器。A = U * S * V^T，其中U和V是正交矩阵，S在对角线上包含奇异值。

数学原理:

SVD是最通用的矩阵分解。奇异值揭示矩阵的秩、条件数，并能够计算秩亏矩阵的伪逆。

矩阵分解方法¶

LU分解¶

Mat::LUDecomposition Mat::lu_decompose(bool use_pivoting = true) const;

描述:

计算LU分解：A = P * L * U（带主元）或A = L * U（不带主元）。对于求解具有相同系数矩阵的多个系统很高效。

数学原理:

• 不带主元: A = L * U，其中L具有单位对角线

• 带主元: P * A = L * U，其中P是置换矩阵

分解通过求解Ly = Pb（前向替换）然后Ux = y（后向替换）来求解Ax = b。

参数:

• bool use_pivoting : 是否使用部分主元以提高数值稳定性（默认：true）。

返回值:

LUDecomposition，包含L、U、P矩阵和状态。

使用建议:

• 多个右端项: 一旦分解，使用solve_lu()高效求解具有不同右端项的多个系统。

• 行列式: det(A) = det(P) * det(L) * det(U) = det(P) * det(U)（因为det(L) = 1）。

• 逆矩阵: 可以通过为每个单位向量eᵢ求解LUx = eᵢ来计算A^(-1)。

• 性能: 分解为O(n³)，分解后每次求解为O(n²)。

坑点:

• 实际工程里不要省略主元: 不带主元的LU会让某些本来可解的矩阵变得非常不稳定。

• L 是单位对角线: 不要把 L 的对角元素当作缩放信息，缩放信息在 U 里。

• 行列式符号会变: 置换矩阵 P 会改变行列式符号。

Cholesky分解¶

Mat::CholeskyDecomposition Mat::cholesky_decompose() const;

描述:

计算对称正定矩阵的Cholesky分解：A = L * L^T。对于SPD矩阵比LU更快，常用于结构动力学。

数学原理:

对于对称正定矩阵A，存在唯一的具有正对角元素的下三角矩阵L，使得A = L * L^T。这本质上是利用对称性的专用LU分解。

参数:

无（矩阵必须是对称正定的）。

返回值:

CholeskyDecomposition，包含L矩阵和状态。

使用建议:

• 效率: 需要大约LU分解的一半计算和存储。

• 稳定性: 对于对称正定矩阵比LU更稳定。

• 应用:

• 结构动力学：质量和刚度矩阵通常是SPD

• 优化：Newton方法中的Hessian矩阵

• 统计：协方差矩阵

• 错误处理: 如果矩阵不对称或不是正定的，返回错误。

坑点:

• 必须同时满足对称和正定: 只有对称还不够，必须是正定。

• 浮点对称性: 数值计算里微小不对称可能导致失败，最好先做容差判断。

• 先验已知时再用: 如果你已经知道矩阵是SPD，这就是正确选择；否则先用 is_symmetric() 和 is_positive_definite() 检查。

QR分解¶

Mat::QRDecomposition Mat::qr_decompose() const;

描述:

计算QR分解：A = Q * R，其中Q是正交的，R是上三角的。数值稳定，用于最小二乘和正交化。

数学原理:

QR分解将任何矩阵表示为正交矩阵Q（Q^T * Q = I）和上三角矩阵R的乘积。使用改进的Gram-Schmidt过程和重新正交化计算分解。

参数:

无。

返回值:

QRDecomposition，包含Q和R矩阵和状态。

使用建议:

• 最小二乘: 对于超定系统Ax ≈ b，最小化||Ax - b||₂的解是x = R^(-1) * Q^T * b。

• 数值稳定性: QR分解对于最小二乘问题比正规方程更稳定。

• 特征分解: QR算法迭代使用QR分解来查找特征值。

• 秩揭示: A的秩等于R的非零对角元素数。

坑点:

• 正交性可能漂移: 基于 Gram-Schmidt 的 QR 方便，但对病态矩阵不如 Householder QR 稳。

• 不是通用逆矩阵工具: QR 适合最小二乘，但不是方阵精确求解的首选，除非你需要它的稳定性特征。

• 秩亏要看容差: R 对角线上很小的值往往意味着近线性相关，必须结合容差判断。

SVD分解¶

Mat::SVDDecomposition Mat::svd_decompose(int max_iter = 100, float tolerance = 1e-6f) const;

描述:

计算奇异值分解：A = U * S * V^T。最通用的分解，用于秩估计、伪逆、降维。使用基于特征分解的迭代方法。

数学原理:

SVD将任何m × n矩阵A分解为：

• U: m × m正交矩阵（左奇异向量）

• S: m × n对角矩阵（奇异值σ₁ ≥ σ₂ ≥ ... ≥ σᵣ ≥ 0）

• V: n × n正交矩阵（右奇异向量）

奇异值揭示矩阵的基本属性：秩、条件数和数值行为。

参数:

• int max_iter : 最大迭代次数（默认：100）。

• float tolerance : 收敛容差（默认：1e-6）。

返回值:

SVDDecomposition，包含U、S、V矩阵、秩和状态。

使用建议:

• 秩估计: 数值秩是高于容差阈值的奇异值数量。

• 伪逆: A⁺ = V * S⁺ * U^T，其中S⁺对于非零σᵢ有1/σᵢ。

• 降维: 截断SVD（仅保留最大奇异值）提供低秩近似。

• 条件数: κ(A) = σ₁ / σᵣ，其中σᵣ是最小的非零奇异值。

• 应用:

• 秩亏系统的最小二乘

• 主成分分析（PCA）

• 图像压缩

• 降噪

坑点:

• 代价高: SVD信息最丰富，但也是最昂贵的分解之一。

• 容差决定秩: 数值秩不是绝对值，它会随着阈值变化。

• 内存占用: 在嵌入式平台上，完整 U 和 V 的存储开销很大。

• 解释要谨慎: 奇异值可以反映条件数和能量集中程度，但不能自动替你完成建模选择。

使用分解求解线性系统¶

使用LU分解求解¶

static Mat Mat::solve_lu(const LUDecomposition &lu, const Mat &b);

描述:

使用预计算的LU分解求解线性系统Ax = b。当求解具有相同系数矩阵的多个系统时，比solve()更高效。

数学原理:

给定A = P * L * U，通过以下方式求解Ax = b：

1. 求解Ly = Pb（前向替换）
2. 求解Ux = y（后向替换）

参数:

• const LUDecomposition &lu : 预计算的LU分解。

• const Mat &b : 右端项向量 (N×1)。

返回值:

Mat - 解向量 (N×1)。

使用建议:

• 多个右端项: 计算一次LU分解后，高效求解多个系统。

• 性能: 每次求解O(n²) vs 完整求解O(n³)，对于多个右端项节省显著。

• 内存: 重用分解，避免重复计算。

坑点:

• 适合多右端项，不适合频繁换A: 如果 A 每次都变，LU 的收益会很快缩小。

• 置换矩阵不能忽略: 求解时必须把 P 考虑进去，否则答案会错。

使用Cholesky分解求解¶

static Mat Mat::solve_cholesky(const CholeskyDecomposition &chol, const Mat &b);

描述:

使用预计算的Cholesky分解求解线性系统Ax = b。对于对称正定矩阵比LU更高效。

数学原理:

给定A = L * L^T，通过以下方式求解Ax = b：

1. 求解Ly = b（前向替换）
2. 求解L^T x = y（后向替换）

参数:

• const CholeskyDecomposition &chol : 预计算的Cholesky分解。

• const Mat &b : 右端项向量 (N×1)。

返回值:

Mat - 解向量 (N×1)。

使用建议:

• 效率: 对于SPD矩阵，在分解和求解方面都比LU更快。

• 稳定性: 对于SPD矩阵在数值上更稳定。

• 应用: 结构动力学、优化、统计。

坑点:

• 仅限SPD: 如果矩阵不是对称正定，这个方法就不合法。

• 先检查再分解: 先做对称性/正定性检查，比事后排查分解失败更省时间。

使用QR分解求解（最小二乘）¶

static Mat Mat::solve_qr(const QRDecomposition &qr, const Mat &b);

描述:

使用QR分解求解线性系统。为超定系统（方程数多于未知数）提供最小二乘解。

数学原理:

对于Ax ≈ b（超定），最小二乘解最小化||Ax - b||₂。使用A = Q * R：

• x = R^(-1) * Q^T * b

这避免了数值不稳定的正规方程A^T * A * x = A^T * b。

参数:

• const QRDecomposition &qr : 预计算的QR分解。

• const Mat &b : 右端项向量 (M×1，其中M ≥ N)。

返回值:

Mat - 最小二乘解向量 (N×1)。

使用建议:

• 超定系统: 处理方程数多于未知数的情况。

• 数值稳定性: 比直接求解正规方程更稳定。

• 应用:

• 曲线拟合

• 数据回归

• 信号处理

坑点:

• 这是最小二乘，不是精确求解: 如果系统不相容，QR给出的是最佳拟合，而不是“精确解”。

• 仍要看秩: R 的对角线上若出现很小的值，说明问题可能病态或秩亏。

伪逆¶

static Mat Mat::pseudo_inverse(const SVDDecomposition &svd, float tolerance = 1e-6f);

描述:

使用SVD分解计算Moore-Penrose伪逆A⁺。适用于秩亏或非方矩阵，其中常规逆不存在。

数学原理:

对于A = U * S * V^T，伪逆为A⁺ = V * S⁺ * U^T，其中S⁺对于奇异值σᵢ > tolerance有1/σᵢ，否则为0。

伪逆的性质:

• A * A⁺ * A = A

• A⁺ * A * A⁺ = A⁺

• (A * A⁺)^T = A * A⁺

• (A⁺ * A)^T = A⁺ * A

参数:

• const SVDDecomposition &svd : 预计算的SVD分解。

• float tolerance : 奇异值阈值（默认：1e-6）。低于此值的奇异值被视为零。

返回值:

Mat - 伪逆矩阵。

使用建议:

• 秩亏系统: 为A不是满秩的系统提供解。

• 最小范数解: 对于欠定系统，给出具有最小||x||₂的解。

• 最小二乘: 对于超定系统，给出最小二乘解。

• 应用:

• 控制系统

• 信号处理

• 机器学习（正则化）

坑点:

• 容差是建模选择: 伪逆结果依赖奇异值截断阈值，阈值一变，实际解就可能变。

• 不一定省内存: 密集伪逆可能破坏结构并增加内存占用。

• 能直接 solve 就别伪逆: 如果矩阵方正且条件良好，直接求解通常更便宜。