跳转至

特征值分解 — 幂迭代、Jacobi、Francis QR

算法流程

幂迭代（Power Iteration）

目标：求矩阵 \(A\) 的模最大特征值 \(\lambda_1\) 及其对应特征向量 \(\mathbf{v}_1\)。

算法：

  1. 初始化随机向量 \(\mathbf{v}^{(0)}\)（本实现使用列绝对值之和）

  2. 迭代 \(k = 0, 1, \dots\)：

  3. \(\mathbf{w}^{(k+1)} = A \cdot \mathbf{v}^{(k)}\)

  4. \(\mathbf{v}^{(k+1)} = \mathbf{w}^{(k+1)} / \|\mathbf{w}^{(k+1)}\|_2\)
  5. \(\lambda^{(k+1)} = \left(\mathbf{v}^{(k+1)}\right)^T \cdot A \cdot \mathbf{v}^{(k+1)}\)（Rayleigh 商）
  6. 收敛判据：\(|\lambda^{(k+1)} - \lambda^{(k)}| < \text{tolerance} \cdot |\lambda^{(k+1)}|\)

收敛速率：线性，比率为 \(|\lambda_2 / \lambda_1|\)。若 \(\lambda_1 \approx \lambda_2\)，收敛极慢。

局限：

  • 只能得到一个特征对
  • 需要 \(|\lambda_1| > |\lambda_2|\)（主特征值唯一）
  • 初始向量不能与 \(\mathbf{v}_1\) 正交

反幂迭代（Inverse Power Iteration）

目标：求矩阵 \(A\) 的模最小特征值 \(\lambda_n\) 及其特征向量（结构动力学中的基频）。

原理：对 \(A^{-1}\) 应用幂迭代。\(A^{-1}\) 的最大特征值是 \(1/\lambda_n\)，因此收敛到 \(\lambda_n\)。

算法：

  1. 初始化 \(\mathbf{v}^{(0)}\)（使用交替符号避免与主特征向量对齐）

  2. 迭代 \(k = 0, 1, \dots\)：

  3. 求解线性系统 \(A \cdot \mathbf{y} = \mathbf{v}^{(k)}\)（即 \(\mathbf{y} = A^{-1} \mathbf{v}^{(k)}\)）

  4. \(\mathbf{v}^{(k+1)} = \mathbf{y} / \|\mathbf{y}\|_2\)
  5. \(\lambda^{(k+1)} = \left(\mathbf{v}^{(k+1)}\right)^T \cdot A \cdot \mathbf{v}^{(k+1)}\)
  6. 收敛判据同幂迭代

关键差异：

  • 每次迭代需要解线性系统（\(O(n^3)\)）vs 幂迭代仅需矩阵乘（\(O(n^2)\)）
  • 要求 \(A\) 可逆（非奇异）
  • 收敛通常更快（靠近零的特征值间隔相对更大）

Jacobi 特征分解（对称矩阵专用）

目标：对对称矩阵 \(A\) 求全部实特征值 \(\lambda_1, \dots, \lambda_n\) 和正交特征向量矩阵 \(V\)。

原理：通过一系列 Givens 旋转（正交相似变换）将 \(A\) 对角化：\(V^T \cdot A \cdot V = \Lambda\)。

算法：

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

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

  3. 若 \(|B_{pq}| > \text{tolerance}\)：

    • 计算旋转角度： [ \tau = \frac{B{qq} - B,\quad t = \frac{\text{sign}(\tau)}{|\tau| + \sqrt{1+\tau^2}},\quad c = \frac{1}{\sqrt{1+t^2}},\quad s = t \cdot c ]}}{2 B_{pq}
    • 应用左右旋转：\(B \leftarrow J^T \cdot B \cdot J\)，\(V \leftarrow V \cdot J\)
    • 重复完整扫描直至 \(\max |B_{pq}| < \text{tolerance}\)
    • \(\lambda_i = B_{ii}\)，特征向量为 \(V\) 的列

特性：

  • 精度高：对称矩阵上正交性保持极好
  • 复杂度：约 \(O(n^3)\) 每次扫描，通常需 \(3{-}5\) 次完整扫描
  • 适合：\(n \leq 10\) 的小型对称矩阵（如结构刚度/质量矩阵）

Francis QR 特征分解（一般矩阵）

目标：对任意方阵 \(A\) 求全部特征值（可能为复数，当前仅返回实部）。

算法流程（简化版：单移位 QR）：

Phase 1 — Hessenberg 约化：

  • 对每一列 \(k = 0, \dots, n-3\)：

  • 构造 Householder 反射将 \(A[k+1..n-1, k]\) 化为 \([\alpha, 0, \dots, 0]^T\)

  • 左右同时反射：\(A \leftarrow H^T \cdot A \cdot H\)
  • 结果：上 Hessenberg 矩阵（次对角线以下全零）

Phase 2 — QR 迭代：

对每个 \(m = n, n-1, \dots, 2\)（逐步 deflate）：

  1. 若 \(|H_{m,m-1}| < \text{tolerance} \cdot (|H_{mm}| + |H_{m-1,m-1}|)\)，视为零 -> deflate
  2. 计算 Wilkinson 移位 \(\mu\)（\(H[m-1..m, m-1..m]\) 更接近 \(H_{mm}\) 的特征值）

  3. 单步 QR：

  4. \([Q, R] = \text{qr}(H[0..m-1, 0..m-1] - \mu I)\)

  5. \(H[0..m-1, 0..m-1] = R \cdot Q + \mu I\)
  6. 累积特征向量：\(V \leftarrow V \cdot Q\)

Phase 3 — 提取特征值：

  • 1×1 块：\(\lambda = H_{ii}\)（实特征值）
  • 2×2 块：计算子块特征值（可能为复数共轭对）；当前实现返回实部

特性：

  • 通用性：可处理非对称矩阵
  • 收敛：带 Wilkinson 移位时通常二次收敛
  • 复杂度：每次迭代 \(O(n^2)\)（利用 Hessenberg 结构），共需 \(O(n^3)\)
  • 局限：复数特征值当前只保留实部；完整双移位版本需要时再扩展

方法选择总结：

场景 推荐方法
仅需最大特征对 power_iteration()
仅需最小特征对（基频） inverse_power_iteration()
对称矩阵，全部特征值 eigendecompose_jacobi()
非对称矩阵，全部特征值 eigendecompose_qr()
不确定矩阵类型 eigendecompose()（自动选择）

线性代数 - 特征值与特征向量

结构体：Mat::EigenPair

Mat::EigenPair::EigenPair();
// fields:
// float eigenvalue;      // 特征值（power_iteration 为最大模，inverse_power_iteration 为最小模）
// Mat eigenvector;       // 对应的特征向量（n x 1）
// int iterations;        // 迭代次数（若为迭代法返回）
// tiny_error_t status;   // 计算状态（TINY_OK / 错误码）

描述:

用于保存单一特征值/特征向量对及其计算信息。常由 power_iteration 或 inverse_power_iteration 返回。

说明:

这个结构刻意保持很小：一个特征值、一个特征向量、再加收敛信息。它适合迭代法只关心一个模态的场景。

结构体：Mat::EigenDecomposition

Mat::EigenDecomposition::EigenDecomposition();
// fields:
// Mat eigenvalues;    // n x 1 矩阵，存放特征值
// Mat eigenvectors;   // n x n 矩阵，通常每一列为对应的特征向量
// int iterations;     // 迭代次数（若为迭代法返回）
// tiny_error_t status; // 计算状态（TINY_OK / 错误码）

描述:

用于保存完整的特征值分解结果，包括全部特征值和对应的特征向量矩阵。

说明:

这是完整的谱信息快照。只有在矩阵规模可控、且你确实需要多个模态或基变换时，它才划算。

幂迭代（求主特征值/向量）

Mat::EigenPair Mat::power_iteration(int max_iter, float tolerance) const;

描述:

使用幂迭代法计算矩阵的主特征值（绝对值最大）及对应特征向量。快速方法，适用于实时SHM应用，可快速识别主频率。

数学原理:

幂迭代通过迭代地将矩阵应用于向量来找到绝对值最大的特征值：

  1. 从随机向量v₀开始

  2. 迭代：vₖ₊₁ = A * vₖ / ||A * vₖ||

  3. 特征值估计：λₖ = (vₖ^T * A * vₖ) / (vₖ^T * vₖ) (Rayleigh商)

收敛性:

方法收敛到主特征值，如果：

  • 主特征值是唯一的 (|λ₁| > |λ₂| ≥ ... ≥ |λₙ|)

  • 初始向量在主特征向量方向上有非零分量

参数:

  • int max_iter : 最大迭代次数（典型默认值：1000）。

  • float tolerance : 收敛容差（例如 1e-6）。收敛性通过 |λₖ - λₖ₋₁| < tolerance * |λₖ| 检查。

返回值:

EigenPair，包含 eigenvalue、eigenvector、iterations 和 status。

使用建议:

  • 实时应用: 对于分离良好的特征值快速收敛，适用于实时结构健康监测。

  • 初始化: 实现使用智能初始化策略（列绝对值之和）以避免收敛到较小的特征值。

  • 收敛速度: 收敛是线性的，速度为 |λ₂|/|λ₁|。当特征值接近时较慢。

  • 局限性:

  • 只找到一个特征值-特征向量对

  • 需要 |λ₁| > |λ₂|（主特征值必须唯一）

  • 如果特征值接近，可能收敛缓慢

  • 应用:

  • 主成分分析（第一主成分）

  • PageRank算法

  • 结构动力学（基频）

坑点:

  • 只能得到一个模态: 它不会告诉你其他特征值。

  • 需要明显主导特征值: 如果最大模特征值不够“突出”，收敛会慢，甚至含糊。

  • 缩放会影响收敛: 缩放很差的矩阵会让收敛速度和向量质量变差。

反幂迭代（求最小特征值/向量）

Mat::EigenPair Mat::inverse_power_iteration(int max_iter, float tolerance) const;

描述:

使用反幂迭代法计算矩阵的最小（最小模）特征值及其对应特征向量。对于系统识别至关重要——在结构动力学中找到基频/最低模态。该方法对于SHM应用至关重要，其中最小特征值对应于系统的基本频率。

数学原理:

反幂迭代将幂迭代应用于A^(-1)， 其特征值为1/λᵢ。由于1/λₙ是A^(-1) 的最大特征值，该方法收敛到A的最小特征值：

  1. 从向量v₀开始

  2. 迭代：求解A * yₖ = vₖ，然后vₖ₊₁ = yₖ / ||yₖ||

  3. 特征值估计：λₖ = (vₖ^T * A * vₖ) / (vₖ^T * vₖ) (Rayleigh商)

收敛性:

收敛到最小特征值，如果：

  • 最小特征值是唯一的 (|λₙ| < |λₙ₋₁| ≤ ... ≤ |λ₁|)

  • 矩阵A是可逆的（非奇异）

  • 初始向量在最小特征向量方向上有分量

参数:

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

  • float tolerance : 收敛容差（默认：1e-6）。使用相对容差：|λₖ - λₖ₋₁| < tolerance * max(|λₖ|, 1.0)。

返回值:

EigenPair，包含最小特征值、特征向量、迭代次数和状态。

算法步骤:

  1. 初始化归一化特征向量v（使用交替符号以避免与主特征向量对齐）

  2. 迭代：使用solve()求解A * y = v（等价于y = A^(-1) * v）

  3. 归一化y得到新的v

  4. 使用Rayleigh商计算特征值估计：λ = (v^T * A * v) / (v^T * v)

  5. 使用相对容差检查收敛性

使用建议:

  • 系统识别: 对于在结构动力学中查找基频至关重要，其中最小特征值对应于最低固有频率。

  • 数值稳定性: 实现包括对奇异矩阵的检查，并优雅地处理近奇异情况。

  • 初始化策略: 使用交替符号模式以避免收敛到较大的特征值，确保收敛到最小特征值。

  • 性能: 每次迭代需要求解线性系统（对于稠密矩阵为O(n³)），但通常比幂迭代收敛更快。

  • 与幂迭代互补:

  • 幂迭代：找到λ_max（最高频率）

  • 反幂迭代：找到λ_min（基频）

  • 两者一起提供系统的频率范围

  • 应用:

  • 结构健康监测（基频检测）

  • 模态分析（最低模态形状）

  • 系统识别

  • 稳定性分析（最小特征值表示稳定性裕度）

注意:

  • 要求矩阵为方阵且数据指针非空；否则返回错误状态。

  • 矩阵必须是可逆的（非奇异）才能使此方法工作。如果矩阵是奇异的或接近奇异的，该方法将优雅地失败。

  • 反幂迭代只返回最小特征值/向量对。若需全部特征值/向量请使用下面的分解函数。

坑点:

  • 要求可逆: 奇异或近奇异矩阵会让每一步求解变得不稳定甚至不可行。

  • 每次迭代都要解线性方程: 所以它只适合你只想要一个极值特征对的场景。

  • 最小模不等于最小数值: 对有正有负的谱来说，“绝对值最小”和“最负”不是一回事。

Jacobi 特征分解（对称矩阵）

Mat::EigenDecomposition Mat::eigendecompose_jacobi(float tolerance, int max_iter) const;

描述:

使用Jacobi方法计算完整的特征值分解。推荐用于对称矩阵（良好的精度和稳定性，适用于结构动力学应用）。稳健且准确，是结构动力学矩阵的理想选择。

数学原理:

Jacobi方法通过一系列正交相似变换（Givens旋转）对角化对称矩阵：

  1. 找到最大的非对角元素aₚq
  2. 计算旋转角θ以将该元素置零
  3. 应用旋转：A' = J^T * A * J，其中J是旋转矩阵
  4. 重复直到所有非对角元素低于容差

收敛性:

当最大非对角元素低于容差时方法收敛。每次旋转将一个非对角元素置零，过程持续直到矩阵对角化。

参数:

  • float tolerance : 收敛阈值（例如 1e-6）。允许的非对角元素最大幅度。

  • int max_iter : 最大迭代次数（例如 100）。对于n×n矩阵，通常需要O(n²)次迭代收敛。

返回值:

EigenDecomposition，包含 eigenvalues、eigenvectors、iterations 和 status。

使用建议:

  • 对称矩阵: 专为对称矩阵设计。对于非对称矩阵，使用QR方法。

  • 数值稳定性: 对于对称矩阵非常稳定，具有良好的正交性保持。

  • 精度: 高精度，适用于需要精确特征值/特征向量对的应用。

  • 性能: 每次迭代O(n³)，但对于对称矩阵通常比QR需要更少的迭代。

  • 应用:

  • 结构动力学：刚度和质量矩阵是对称的

  • 主成分分析（PCA）

  • 谱聚类

  • 二次型优化

注意:

如果矩阵不是近似对称，函数会发出警告，但仍可能运行。对于非对称矩阵，推荐使用QR方法。

坑点:

  • 对称性是前提: 该方法依赖正交相似变换。非对称输入会破坏算法假设。

  • 迭代次数有限: 很小的容差可能需要很多轮旋转，注意 max_iter 是否足够。

  • 对角化后仍是近似值: 浮点误差会留下极小的非对角项，要把结果当作近似对角化。

QR 特征分解（一般矩阵）

Mat::EigenDecomposition Mat::eigendecompose_qr(int max_iter, float tolerance) const;

描述:

使用QR算法计算特征值分解。适用于一般（可能非对称）矩阵。支持非对称矩阵，但可能产生复数特征值（仅返回实部）。

数学原理:

QR算法迭代应用QR分解：

  1. 从A₀ = A开始

  2. 对于k = 0, 1, 2, ...： 计算QR分解：Aₖ = Qₖ * Rₖ，然后更新：Aₖ₊₁ = Rₖ * Qₖ

  3. Aₖ收敛到上三角形式（Schur形式），特征值在对角线上

收敛性:

当Aₖ近似上三角（次对角元素 < 容差）时算法收敛。特征值出现在对角线上，特征向量从Q矩阵累积。

参数:

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

  • float tolerance : 收敛容差（例如 1e-6）。使用相对容差，比较次对角元素与对角元素。

返回值:

EigenDecomposition，包含特征值、特征向量、迭代次数和状态。

使用建议:

  • 一般矩阵: 可以处理非对称矩阵，与Jacobi方法不同。

  • 复数特征值: 非对称矩阵可能具有复数特征值；当前实现仅返回实部。

  • 数值稳定性: 使用改进的Gram-Schmidt和重新正交化以提高稳定性。

  • 性能: 每次迭代O(n³)。可能需要多次迭代才能收敛，特别是对于病态矩阵。

  • 收敛加速: 实现可以从移位（Wilkinson移位）中受益以加快收敛，但当前版本使用基本QR迭代。

  • 应用:

  • 一般矩阵特征值问题

  • 动力系统分析

  • 控制理论（系统极点）

注意:

QR在此实现中使用Gram–Schmidt构造Q/R；对于病态矩阵可能不太稳定。对于对称矩阵，由于更好的稳定性和精度，推荐使用Jacobi。

坑点:

  • 复谱限制: 当前接口只保留特征值实部，不是完整的复数特征值求解器。

  • 收敛可能较慢: 如果没有位移策略，困难矩阵上 QR 可能需要很多轮迭代。

  • 正交化质量很重要: Q 用 Gram-Schmidt 构造，病态输入会影响迭代质量。

自动特征分解（根据矩阵特性选择方法）

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

描述:

简便接口，会自动根据矩阵特性选择最优算法。先调用 is_symmetric(tolerance * 10.0f) 判断矩阵是否近似对称：

  • 若为对称，使用 eigendecompose_jacobi(tolerance, max_iter)（更稳定和精确）；
  • 否则使用 eigendecompose_qr(max_iter, tolerance)（处理一般矩阵）。

算法选择流程:

  1. 测试矩阵是否对称：is_symmetric(tolerance * 10.0f)
  2. 若对称 → 使用 eigendecompose_jacobi(tolerance, max_iter)
  3. 若非对称 → 使用 eigendecompose_qr(max_iter, tolerance)

参数:

  • float tolerance：用于对称性检测与分解收敛判断（默认 1e-6f）。
  • int max_iter：最大迭代次数（必须 > 0，默认值 = 100）。

返回值:

EigenDecomposition 结构，包含所有特征值和特征向量。

使用建议:

  • 自动优化：无需手动选择算法，同时仍能获得最佳性能。

  • 边缘计算：非常适合嵌入式系统，无需手动调优即可获得良好性能。

  • 鲁棒性：对称性测试使用放宽的容差（10×），以处理数值误差，确保正确识别对称矩阵。

实用技巧:

  • 已知对称性：若明确知道矩阵为对称矩阵（如刚度矩阵、质量矩阵），直接使用 eigendecompose_jacobi 可获得最佳稳定性和略好的性能。

  • 未知特性：对于一般矩阵或对称性未知的情况，使用 eigendecompose 进行自动选择。

  • 性能考虑：

  • 在嵌入式平台上，对大型矩阵进行特征分解的计算成本很高

  • 对于 n > 20，当只需要少数特征值时，考虑使用降阶方法或迭代方法（幂迭代）

  • 对于实时应用，使用 power_iteration() 或 inverse_power_iteration() 计算单个特征值

  • 内存使用：完全特征分解需要存储所有特征向量（n×n 矩阵），对于大型矩阵可能占用大量内存。

坑点:

  • 自动不等于免费: 对称性测试和分解本身都要耗时。若你已经知道矩阵类别，直接调用专用方法更好。

  • 容差有双重作用: 它同时影响对称性判定和收敛条件，所以调它会影响算法选择和数值停止条件。

  • 大矩阵不友好: 在嵌入式平台上，完整特征分解很容易成为时间和内存瓶颈。