特征值 - 特征值的平衡矩阵

Eigen - Balancing matrix for eigenvalue

本文关键字：特征值 平衡      更新时间：2023-10-16

我的经验(像其他人一样：如何使用LAPACK从矩阵对的广义舒尔分解中获得指定的特征向量？(是，当矩阵条件不佳时，从特征向量获得的特征值(我不关心特征向量(几乎不如从numpy，matlab等获得的特征值可靠。

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

目前，我有一个烦人的两层解决方案，涉及python和C++，我想将所有内容都推到C++中;特征值求解器是唯一阻碍我的部分。

事实证明，这篇关于arxiv的精彩小论文对平衡进行了很好的清晰描述：https://arxiv.org/pdf/1401.5766.pdf。 当我实现这种平衡时，特征值几乎与 numpy 完全一致。 如果特征在获取特征值之前平衡矩阵，那就太好了。

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);
}
相关文章：