Eigen - balancing matrix for eigenvalue
My experience (as well as some others: How do I get given eigenvectors from a generalized Schur factorization of a pair of matrices using LAPACK? ) Is that the eigenvalues obtained from Eigen (I don't care about eigenvectors) are not as reliable as those which are derived from numpy, matlab, etc. when the matrix is badly conditioned.
The internet ( https://www.mathworks.com/help/matlab/ref/balance.html ) suggests balancing is a solution, but I can't figure out how to do it in Eigen, Can anyone help?
I currently have an annoying two-tier solution that includes python and C ++ and I would like to push everything to C ++; the eigenvalue solver is the only part holding me back.
source to share
As it turns out, this fantastic little article on arxiv gives a nice clear description of balancing: https://arxiv.org/pdf/1401.5766.pdf . When I implement this balancing, the eigenvalues are almost consistent with numpy. It would be great if Eigen would balance the matrix before accepting the 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);
}
source to share