Найдите ближайшую положительно определенную матрицу с помощью Eigen - PullRequest
Новая
       71

Найдите ближайшую положительно определенную матрицу с помощью Eigen

0 голосов
/ 06 мая 2020

Я хочу найти ближайшую положительно определенную матрицу к некоторой матрице X.

Я видел здесь этот ответ: https://scicomp.stackexchange.com/questions/30631/how-to-find-the-nearest-a-near-positive-definite-from-a-given-matrix

Но у меня проблема с переводом этого в ответ на Eigen. Я пробовал следующее: 

Matrix<double, 9, 9> X, D, DPLUS, Z, Y;
X.setRandom();
Y = 0.5 * (X + X.transpose().eval());
EigenSolver<Matrix<double, 9, 9>> es;
es.compute(Y);
D = es.eigenvalues().asDiagonal();
DPlus = D.cwiseMax(0);
Z = es.eigenvectors() * DPLUS * es.eigenvectors().transpose().eval();

Это дает мне ошибки, связанные со сложной проблемой, но соответствует ли это тому, что предлагает связанный ответ?

1 Ответ

1 голос
/ 07 мая 2020

Проблемы, в которых вы не используете правильный решатель собственных значений.

Поскольку матрица Y симметрична c по конструкции, вам необходимо использовать SelfAdjointEigenSolver. Эта симметрия гарантирует, что ваша реальная матрица диагонализируется и имеет действительные собственные значения и собственные векторы, поэтому вам не нужно иметь дело с комплексными числами.

Общий EigenSolver не делает этого предположения и, следовательно, не в этом случае будет работать так, как вы собираетесь. Вы можете получить с ним тот же результат (вы получите комплексные собственные значения, которые будут фактически действительными, но вам нужно будет создать векторы комплексов et c. Кроме того, это намного менее эффективно).

Это даст вам собственные значения и собственные векторы, которые вы можете использовать для построения своей ортогональной базисной матрицы Q . 

const MatrixXd Y = 0.5 * (X + X.transpose());
const SelfAdjointEigenSolver<MatrixXd> solver(Y);
const VectorXd D = solver.eigenvalues();
const MatrixXd Q = solver.eigenvectors();
const VectorXd Dplus = D.cwiseMax(0);
const MatrixXd Z = Q * Dplus.asDiagonal() * Q.transpose();

Конечно, если вам нравятся однострочники, вы всегда можете сократить его до: 

const SelfAdjointEigenSolver<MatrixXd> solver(0.5 * (X + X.transpose()));
const MatrixXd Z = solver.eigenvectors() * solver.eigenvalues().cwiseMax(0).asDiagonal() * solver.eigenvectors().transpose();
...