Eigen - Balancing matrix for eigenvalue

Question

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.

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?

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.

Answer 1

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