在Eigen3中实现Bartels–Stewart算法? [英] Implementing the Bartels–Stewart algorithm in Eigen3?
问题描述
过去当我需要解 Sylvester 方程时,AX + XB = C
,我使用了 scipy
的函数,solve_sylvester
[1],显然是通过使用Bartels-Stewart算法将事物分解为上三角形式,然后使用 lapack
求解方程来实现的.
In the past when I've needed to solve the Sylvester equation, AX + XB = C
, I've used scipy
's function, solve_sylvester
[1], which apparently works by using the Bartels-Stewart algorithm to get things into upper triangular form, and then solving the equation using lapack
.
我现在需要使用 eigen
求解方程. eigen
提供了一个函数 matrix_function_solve_triangular_sylvester
[2],在文档中看起来与 scipy
的 lapack
函数相似.代码>调用.我试图在 eigen3
中准确地翻译 scipy
的实现,但是最后我对 X
的值不满足该等式.这是我的实现:
I now need to solve the equation using eigen
. eigen
provides an function, matrix_function_solve_triangular_sylvester
[2], which seems by the documentation to be similar to the lapack
function which scipy
calls. I'm attempting to translate exactly scipy
's implementation in eigen3
, but in the end my value for X
doesn't satisfy the equation. Here's my implementation:
#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;
}
有什么想法我在做什么错吗?
Any ideas what I'm doing wrong?
[1] https://github.com/scipy/scipy/blob/v0.15.1/scipy/linalg/_solvers.py#L23
[2] 推荐答案
@chtz是正确的;这是由于Schurr分解矩阵是准三角形而不是三角形.您使用的特征求解器无法处理此类矩阵.但是,chtz是错误的,因为有些西尔维斯特(Sylvester)求解器可以处理准三角求解器.例如,lapack的 trsyl 可以交易与准三角矩阵.这就是 @chtz is correct; this is due to the Schurr decomposition matrix being quasi-triangular rather than triangular. The eigen solver you are using cannot deal with such matrices. However, chtz is wrong in that there are Sylvester solvers that can deal with quasi-triangular solvers. For example lapack's trsyl can deal with quasi-triangular matrices. This is what is called by 这篇关于在Eigen3中实现Bartels–Stewart算法?的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋! scipy
所称的,它解释了OP的scipy实现工作正常的原因.scipy
, which explains why the OP's scipy implementation worked fine.