在 Eigen3 中实现 Bartels–Stewart 算法？

过去当我需要解西尔维斯特方程时AX + XB = C,我用过scipy的函数，solve_sylvester[1]，它显然是通过使用 Bartels-Stewart 算法将事物转化为上三角形式，然后使用lapack.

我现在需要使用以下方法求解方程eigen. eigen提供一个函数，matrix_function_solve_triangular_sylvester[2]，从文档来看，它类似于lapack函数其中scipy来电。我正在尝试准确翻译scipy的实现在eigen3，但最终我的价值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 https://github.com/scipy/scipy/blob/v0.15.1/scipy/linalg/_solvers.py#L23

[2] https://bitbucket.org/eigen/eigen/src/dbb0b1f3b07a261d01f43f8fb94e85ceede9fac7/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h?at=default#lines-274 https://bitbucket.org/eigen/eigen/src/dbb0b1f3b07a261d01f43f8fb94e85ceede9fac7/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h?at=default#lines-274

Your A and B矩阵具有非实特征值，因此它们RealSchur分解将是非三角形的（仅是“准三角形”，即它在对角线上包含一个 2x2 块）。如果你编译时没有-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.

我不知道是否有一个西尔维斯特求解器也可以处理准三角矩阵，但使用特征方法的最简单的解决方案是使用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;

I think X应始终为实数，因此您可以将最后两行替换为

Eigen::MatrixXd X = ((U * Y) * V.adjoint()).real();

Eigen::MatrixXd Q_calc = A * X + X * B;