SymGEigsShiftSolver与SymGEigsSolver选型及性能对比技术咨询

核心求解目标差异 #

• SymGEigsSolver：直接求解广义特征值问题 (Ax = \lambda Bx) 中两端极值的特征值（最大模、最小模、最大实部、最小实部等），无需指定偏移值。
• SymGEigsShiftSolver：通过移位-逆变换求解靠近特定偏移值σ的特征值，将原问题转化为求解 ((A - \sigma B)^{-1}Bx = \mu x) 的主导特征值，原特征值λ满足 (\mu = 1/(\lambda - \sigma))。

编程实操差异 #

1. 初始化参数

   • SymGEigsSolver需指定特征值排序规则（如SortRule::LargestMagnitude找最大模特征值），无需额外偏移参数。
   • SymGEigsShiftSolver必须传入偏移值σ，同时需结合排序规则明确目标特征值与σ的相对位置（如SortRule::LargestReal对应找比σ大的特征值，SortRule::SmallestReal对应找比σ小的）。

2. 矩阵算子要求

   • SymGEigsSolver仅需实现矩阵-向量乘法 (Av) 和 (Bv)，逻辑简单，无需处理逆运算。
   • SymGEigsShiftSolver需要实现移位逆算子：即计算 ((A - \sigma B)^{-1}Bv)。通常需先对 (A - \sigma B) 做分解（如Cholesky、LU），再通过解线性方程组完成逆运算，代码复杂度更高。

适用场景 #

• 优先用SymGEigsSolver的场景：

  • 需要获取最大/最小极值特征值（如结构力学固有频率、PCA主成分分析）。
  • 特征值分布稀疏，两端特征与其他特征差距明显，无需额外移位即可快速收敛。
  • 不想处理逆运算，追求代码简洁性。

• 优先用SymGEigsShiftSolver的场景：

  • 需要获取中间区域的特征值（如量子力学中特定能量附近的能级、优化问题中特定区间的特征值）。
  • 特征值分布密集，两端特征之外的目标特征用普通Solver收敛极慢，移位变换可将目标特征转化为迭代的主导项。

稳定性与收敛速度 #

• 收敛速度：

  • 当目标特征值靠近偏移值σ时，Shift Solver收敛速度远快于普通Solver——移位逆变换会放大目标特征值对应的迭代分量，减少迭代次数。
  • 若目标是两端极值特征，普通Solver更高效，无需额外逆运算的开销。

• 稳定性：

  • SymGEigsSolver的稳定性依赖于矩阵A、B的条件数，当B正定且A对称时，整体数值稳定性较好。
  • SymGEigsShiftSolver的稳定性受σ的选择影响极大：若σ离某个特征值过近，(A - \sigma B) 接近奇异，逆运算的数值误差会急剧增大，导致结果不稳定；若σ选择合理（与目标特征值距离适中），稳定性与普通Solver相当。

代码示例片段 #

SymGEigsSolver 实现 #

#include <Spectra/SymGEigsSolver.h>

// 定义A的矩阵-向量乘算子
class AOp {
public:
    using Scalar = double;
    int rows() const { return 100; }
    int cols() const { return 100; }
    void perform_op(const double* x, double* y) const {
        // 实现 y = A * x 的计算逻辑
    }
};

// 定义B的矩阵-向量乘算子
class BOp {
public:
    using Scalar = double;
    int rows() const { return 100; }
    int cols() const { return 100; }
    void perform_op(const double* x, double* y) const {
        // 实现 y = B * x 的计算逻辑
    }
};

int main() {
    AOp op;
    BOp bop;
    // 找最大的3个特征值，子空间维度设为6
    SymGEigsSolver<AOp, BOp, Spectra::SortRule::LargestMagnitude> solver(op, bop, 3, 6);
    solver.init();
    int nconv = solver.compute();
    // 获取结果
    auto vals = solver.eigenvalues();
    return 0;
}

SymGEigsShiftSolver 实现 #

#include <Spectra/SymGEigsShiftSolver.h>
#include <Eigen/Cholesky> // 假设用Eigen做矩阵分解

using Matrix = Eigen::MatrixXd;
const Matrix A = ...; // 初始化对称矩阵A
const Matrix B = ...; // 初始化正定矩阵B
const double sigma = 5.0; // 偏移值
const Matrix A_sigma_B = A - sigma * B;
const Eigen::LLT<Matrix> llt(A_sigma_B); // 对A-σB做Cholesky分解

// 定义移位逆算子：(A-σB)^{-1}Bv
class ShiftOp {
public:
    using Scalar = double;
    int rows() const { return A.rows(); }
    int cols() const { return A.cols(); }
    void perform_op(const double* x, double* y) const {
        Eigen::Map<const Eigen::VectorXd> vec_x(x, rows());
        Eigen::Map<Eigen::VectorXd> vec_y(y, rows());
        vec_y = llt.solve(B * vec_x); // 计算 (A-σB)^{-1} * Bx
    }
};

int main() {
    ShiftOp shift_op;
    // 找靠近sigma=5.0的3个特征值，子空间维度设为6
    SymGEigsShiftSolver<ShiftOp, Spectra::SortRule::LargestMagnitude> solver(shift_op, 3, 6, sigma);
    solver.init();
    int nconv = solver.compute();
    auto vals = solver.eigenvalues();
    // 原问题特征值为 sigma + 1/vals[i]
    return 0;
}

内容的提问来源于stack exchange，提问作者user8469759