我的经验(像其他一些：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)表明平衡是解决方案,但我无法弄清楚如何在Eigen中做到这一点.有人可以帮忙吗？

目前我有一个讨厌的双层解决方案,涉及python和C,我想把所有东西都推到C;特征值求解器是阻碍我的唯一部分.

解决方法:

事实证明,这篇关于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);
}