在eigen3中实施Bartels -Stewart算法

Question

过去我需要求解Sylvester方程时，AX + XB = C时，我已经使用了scipy的函数solve_sylvester [1]，显然可以使用Bartels-Stewart来工作算法将物体变成上层三角形形式，然后使用lapack。

求解方程

我现在需要使用eigen求解方程。eigen提供了一个函数matrix_function_solve_triangular_sylvester [2]，该文档似乎与scipy调用的lapack函数相似。我试图在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

[2] https://bitbucket.org/eigen/eigen/src/dbb0b1f3b07a261d01f43f43fb94e855cescee9fac7fac7/

您的 A和 B矩阵具有非真实特征值，因此其RealSchur分解将是非三角形的。如果您没有-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-solver，它也可以处理准三角矩阵，但是使用特征方法的最简单解决方案是使用ComplexSchur分解(也使用adjoint()代替transpose()-和DON''t 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;

@chtz是正确的；这是由于Schurr分解矩阵是准三角形而不是三角形的。您使用的特征求解器无法处理此类矩阵。但是，CHTZ是错误的，因为有一些Sylvester求解器可以处理准三角形求解器。例如，拉帕克(Lapack(的trsyl可以处理准三角矩阵。这就是scipy所谓的，它解释了为什么OP的Scipy实现效果很好。