如何用Eigen(C++)复现C库的矩阵Hessenberg分解与QR特征值计算？

差异根源 #

• 分解定义不一致：原C库的Hessenberg_Form_Elementary可能实现的是下Hessenberg矩阵，而Eigen的HessenbergDecomposition默认生成上Hessenberg矩阵；另外变换矩阵的构造逻辑也可能不同——比如原库是左乘变换、Eigen是右乘，或者原库返回的是变换的逆矩阵（转置/伴随）。
• 反射向量符号差异：Eigen在Householder反射向量的构造中，会强制第一个元素为正以保证数值稳定性，但旧C库可能没有这个处理，导致变换矩阵和Hessenberg矩阵的对应符号完全相反。
• 变换矩阵累积顺序：原库的SSS可能是按从右到左的顺序累积Householder变换，而Eigen的matrixQ()是从左到右累积，最终矩阵的乘积顺序相反。
• QR迭代细节不同：QR_Hessenberg_Matrix的后续运算中，Eigen的优化QR实现（含Wilkinson移位、收敛判断）和旧库的迭代步骤、阈值参数不一致，导致结果偏差。

Eigen的适配实现步骤 #

• 切换Hessenberg类型：若原库输出下Hessenberg矩阵，可通过矩阵转置适配：

  Eigen::MatrixXd A = ...; // 输入原矩阵
  Eigen::HessenbergDecomposition<Eigen::MatrixXd> hess_decomp(A.transpose());
  Eigen::MatrixXd H = hess_decomp.matrixH().transpose(); // 得到下Hessenberg矩阵
  Eigen::MatrixXd Q = hess_decomp.matrixQ().transpose(); // 对应变换矩阵
  

• 修正变换矩阵形式：如果原库的SSS是正交变换的逆（即Q^T），直接调用hess_decomp.matrixQ().adjoint()（实矩阵等价于转置）即可。
• 对齐反射向量符号：遍历Hessenberg分解的反射向量，强制符号和原库一致：

  Eigen::HessenbergDecomposition<Eigen::MatrixXd> hess_decomp(A);
  auto& householder_q = hess_decomp.householderQ();
  for (int i = 0; i < A.rows() - 2; ++i) {
      Eigen::VectorXd vec = householder_q.householderVector(i);
      // 根据原库的符号调整，比如原库第一个元素为负则取反
      if (vec(0) * original_sign < 0) vec = -vec;
      // 重新构造变换矩阵（可参考Eigen HouseholderSequence的逻辑）
  }
  

• 手动对齐QR迭代：若后续QR运算差异大，可基于Eigen的Hessenberg矩阵手动实现QR迭代，照搬原库的移位策略（比如固定移位、Wilkinson移位）和收敛阈值，确保每一步运算逻辑对齐。

调用原C库的取舍 #

• 如果对数值一致性要求极高（比如遗留系统兼容），且Eigen的适配无法完全匹配，可通过extern "C"封装原库函数，对接Eigen矩阵：

  extern "C" {
      void Hessenberg_Form_Elementary(double* A, int n, double* E, double* SSS);
  }
  
  void wrap_hess_c(Eigen::MatrixXd& A, Eigen::MatrixXd& E, Eigen::MatrixXd& SSS) {
      int n = A.rows();
      Hessenberg_Form_Elementary(A.data(), n, E.data(), SSS.data());
  }
  

• 优先推荐用Eigen实现，它的数值稳定性更强，长期维护成本更低，只要对齐分解的核心定义和细节，完全可以匹配原库的结果。

内容的提问来源于stack exchange，提问作者Oláh Balázs