过去,当我需要求解西尔维斯特方程时,AX + XB = C
,我使用了 scipy
的函数,solve_sylvester
[1],这显然是通过使用 Bartels-Stewart 算法将事物转化为上三角形式,然后使用 lapack
求解方程。
我现在需要使用 eigen
求解方程。 eigen
提供了一个函数,matrix_function_solve_triangular_sylvester
[2],它在文档中看起来类似于 lapack
函数,scipy
调用。我试图在 eigen3
中准确翻译 scipy
的实现,但最终我对 X
的值不满足等式。这是我的实现:
#include <iostream>
#include <Eigen/Core>
#include <Eigen/Eigenvalues>
#include <unsupported/Eigen/MatrixFunctions>
int main()
{
Eigen::Matrix<double, 3, 3> A;
A << -17, -6, 0,
-15, 6, 14,
9, -12, 19;
Eigen::Matrix<double, 5, 5> B;
B << 5, -17, -12, 16, 11,
-4, 19, -1, 9, 13,
1, 3, 5, -5, 2,
8, -15, 5, 14, -12,
-2, -4, 13, -8, -17;
Eigen::Matrix<double, 3, 5> Q;
Q << 6, 5, -17, 12, 4,
-11, 15, 8, 1, 7,
15, -3, 9, -19, -10;
Eigen::RealSchur<Eigen::MatrixXd> SchurA(A);
Eigen::MatrixXd R = SchurA.matrixT();
Eigen::MatrixXd U = SchurA.matrixU();
Eigen::RealSchur<Eigen::MatrixXd> SchurB(B.transpose());
Eigen::MatrixXd S = SchurB.matrixT();
Eigen::MatrixXd V = SchurB.matrixU();
Eigen::MatrixXd F = (U.transpose() * Q) * V;
Eigen::MatrixXd Y =
Eigen::internal::matrix_function_solve_triangular_sylvester(R, S, F);
Eigen::MatrixXd X = (U * Y) * V.transpose();
Eigen::MatrixXd Q_calc = A * X + X * B;
std::cout << Q_calc - Q << std::endl;
// Should be all zeros, but instead getting:
// 421.868 193.032 -208.273 42.7449 -3.57527
//-1651.66 -390.314 2043.59 -1611.1 -1843.91
//-67.4093 207.414 1168.89 -1240.54 -1650.48
return EXIT_SUCCESS;
}
知道我做错了什么吗?
[1] https://github.com/scipy/scipy/blob/v0.15.1/scipy/linalg/_solvers.py#L23
最佳答案
您的 A
和 B
矩阵具有非实数特征值,因此它们的 RealSchur
分解将是非三角形的(仅“准三角形” ",即它在对角线上包含一个 2x2 block )。如果你编译时没有使用 -DNDEBUG
,你应该得到这样的断言:
../eigen/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h:277: MatrixType Eigen::internal::matrix_function_solve_triangular_sylvester(const MatrixType&, const MatrixType&, const MatrixType&) [with MatrixType = Eigen::Matrix<double, -1, -1>]: Assertion `A.isUpperTriangular()' failed.
我不知道,如果有一个 Sylvester 求解器也可以处理类三角矩阵,但使用 Eigen 方法的最简单的解决方案是使用 ComplexSchur
分解(也使用 adjoint()
而不是 transpose()
—— 并且不要转置 B
:
Eigen::ComplexSchur<Eigen::MatrixXd> SchurA(A);
Eigen::MatrixXcd R = SchurA.matrixT();
Eigen::MatrixXcd U = SchurA.matrixU();
Eigen::ComplexSchur<Eigen::MatrixXd> SchurB(B);
Eigen::MatrixXcd S = SchurB.matrixT();
Eigen::MatrixXcd V = SchurB.matrixU();
Eigen::MatrixXcd F = (U.adjoint() * Q) * V;
Eigen::MatrixXcd Y =
Eigen::internal::matrix_function_solve_triangular_sylvester(R, S, F);
Eigen::MatrixXcd X = (U * Y) * V.adjoint();
Eigen::MatrixXcd Q_calc = A * X + X * B;
我认为 X
应该总是真实的,所以你可以将最后两行替换为
Eigen::MatrixXd X = ((U * Y) * V.adjoint()).real();
Eigen::MatrixXd Q_calc = A * X + X * B;
关于linear-algebra - 在 Eigen3 中实现 Bartels–Stewart 算法?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/56929966/