Eigen - 特徵值的平衡矩陣

Question

我的經驗（像其他一些：How do I get specified Eigenvectors from the generalized Schur factorization of a matrix pair using LAPACK?）是從Eigen（我不關心特徵向量）獲得的特徵值幾乎不像​​從numpy，matlab等獲得的特徵值一樣可靠。當矩陣病態時。

互聯網（https://www.mathworks.com/help/matlab/ref/balance.html）認爲平衡是解決方案，但我想不出如何徵做到這一點。誰能幫忙？

目前，我有一個涉及Python和C++，我想一切都推到C++一個惱人的雙層解決方案;特徵值解算器是阻止我回歸的唯一部分。

Answer 1

事實證明，上的arXiv這個夢幻般的小給出了平衡的一個很好的清晰描述：https://arxiv.org/pdf/1401.5766.pdf。當我實現這種平衡時，特徵值與numpy幾乎完全一致。如果Eigen在採用特徵值之前平衡矩陣，那將是非常好的。

void balance_matrix(const Eigen::MatrixXd &A, Eigen::MatrixXd &Aprime, Eigen::MatrixXd &D) { 
    // https://arxiv.org/pdf/1401.5766.pdf (Algorithm #3) 
    const int p = 2; 
    double beta = 2; // Radix base (2?) 
    Aprime = A; 
    D = Eigen::MatrixXd::Identity(A.rows(), A.cols()); 
    bool converged = false; 
    do { 
     converged = true; 
     for (Eigen::Index i = 0; i < A.rows(); ++i) { 
      double c = Aprime.col(i).lpNorm<p>(); 
      double r = Aprime.row(i).lpNorm<p>(); 
      double s = pow(c, p) + pow(r, p); 
      double f = 1; 
      while (c < r/beta) { 
       c *= beta; 
       r /= beta; 
       f *= beta; 
      } 
      while (c >= r*beta) { 
       c /= beta; 
       r *= beta; 
       f /= beta; 
      } 
      if (pow(c, p) + pow(r, p) < 0.95*s) { 
       converged = false; 
       D(i, i) *= f; 
       Aprime.col(i) *= f; 
       Aprime.row(i) /= f; 
      } 
     } 
    } while (!converged); 
}