固有値の固有ベクトル行列

Question

私の経験（他のいくつかと同様：How do I get specified Eigenvectors from the generalized Schur factorization of a matrix pair using LAPACK?）は、固有値から得られる固有値（私は固有ベクトルを気にしません）がnumpy、matlabなどから得られるものほど信頼性が低いということです。行列が悪条件調整されているとき。

インターネット（https://www.mathworks.com/help/matlab/ref/balance.html）は、バランシングが解決策であることを示唆していますが、Eigenでこれを行う方法を理解できません。誰も助けることができますか？

現時点では、私はpythonとC++を含む迷惑な2層ソリューションを持っており、すべてを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); 
}