在 Eigen3 中实现 Bartels–Stewart 算法？

Question

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

求解方程

我现在需要使用 eigen 求解方程。 eigen 提供了一个函数，matrix_function_solve_triangular_sylvester [2]，从文档看来它类似于 scipy 调用的 lapack 函数。我试图准确翻译 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

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

Answer 1

您的 A 和 B 矩阵具有 non-real 个特征值，因此它们的 RealSchur 分解将是 non-triangular（仅 "quasi-triangular"，即，它在对角线上包含一个 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.

我不知道，如果有一个 Sylvester-solver 也处理 quasi-triangular 矩阵，但是使用 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;

Answer 2

@chtz 是正确的；这是因为 Schurr 分解矩阵是 quasi-triangular 而不是三角形。您使用的本征求解器无法处理此类矩阵。但是，chtz 是错误的，因为有 Sylvester 求解器可以处理 quasi-triangular 求解器。例如 lapack 的 trsyl 可以处理 quasi-triangular 矩阵。这就是 scipy 所称的，这解释了为什么 OP 的 scipy 实施工作正常。

在 Eigen3 中实现 Bartels–Stewart 算法？

Implementing the Bartels–Stewart algorithm in Eigen3?

linear-algebra

eigen3