c – 特征值的特征平衡矩阵
作者:互联网
我的经验(像其他一些:How do I get specified Eigenvectors from the generalized Schur factorization of a matrix pair using LAPACK?)是从Eigen(我不关心特征向量)获得的特征值几乎不像从numpy,matlab等获得的特征值那样,当矩阵病态时.
互联网(https://www.mathworks.com/help/matlab/ref/balance.html)表明平衡是解决方案,但我无法弄清楚如何在Eigen中做到这一点.有人可以帮忙吗?
目前我有一个讨厌的双层解决方案,涉及python和C,我想把所有东西都推到C;特征值求解器是阻碍我的唯一部分.
解决方法:
事实证明,这篇关于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);
}
标签:c,eigen,eigen3 来源: https://codeday.me/bug/20190828/1752934.html