c++ - Eigen - 特征值的平衡矩阵

Question

我的经验(和其他人一样:How do I get specified Eigenvectors from the generalized Schur factorization of a matrix pair using LAPACK?)是从 Eigen 获得的特征值(我不关心特征向量)不如从 numpy、matlab 等获得的特征值可靠。当矩阵是 body 状况不佳。

互联网 ( https://www.mathworks.com/help/matlab/ref/balance.html ) 表明平衡是解决方案，但我不知道如何在 Eigen 中做到这一点。谁能帮忙？

目前我有一个烦人的两层解决方案，涉及 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);
}