特征值-特征值平衡矩阵 [英] Eigen - Balancing matrix for eigenvalue

问题描述

我的经验（和其他一些经验一样：如何使用LAPACK从矩阵对的广义Schur分解中获得特定的特征向量？）是从特征值获得的特征值（我不在乎特征向量）不如获得的可靠numpy，matlab等从矩阵中调出条件。

My experience (like some others: How do I get specified Eigenvectors from the generalized Schur factorization of a matrix pair using LAPACK?) is that the eigenvalues obtained from Eigen (I don't care about the eigenvectors) are not nearly as reliable as those obtained from numpy, matlab, etc. when the matrix is ill-conditioned.

互联网（ https://www.mathworks.com/help/matlab/ref/balance.html ）建议平衡是解决方案，但我可以请弄清楚如何在Eigen中执行此操作。有人可以帮忙吗？

The internet (https://www.mathworks.com/help/matlab/ref/balance.html) suggests that balancing is the solution, but I can't figure out how to do this in Eigen. Can anyone help?

目前，我有一个烦人的涉及Python和C ++的两层解决方案，我想将所有内容都推送到C ++中。

At the moment I have an annoying two-layer solution that involves python and C++ and I would like to push everything into C++; the eigenvalue solver is the only part that is holding me back.

推荐答案

事实证明，这份关于arxiv的奇妙小论文给出了一个很好的平衡说明： https://arxiv.org/pdf/1401.5766.pdf 。当我实现这种平衡时，特征值几乎与numpy一致。如果Eigen在获取特征值之前先平衡矩阵。

As it turns out, this fantastic little paper on arxiv gives a nice clear description of balancing: https://arxiv.org/pdf/1401.5766.pdf. When I implement this balancing, the eigenvalues agree almost perfectly with numpy. It would be great if Eigen would balance the matrix prior to taking eigenvalues.

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);
}

这篇关于特征值-特征值平衡矩阵的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！