Implementing the Bartels–Stewart algorithm in Eigen3?

Question

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.

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] https://bitbucket.org/eigen/eigen/src/dbb0b1f3b07a261d01f43f8fb94e85ceede9fac7/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h?at=default#lines-274

Answer 1

Your A and B matrices have non-real eigenvalues, therefore their RealSchur decomposition will be non-triangular (only "quasi-triangular", i.e., it contains a 2x2 block on the diagonal). If you compile without -DNDEBUG, you should get an assertion like this:

../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.

I don't know, if there is a Sylvester-solver which also handles quasi-triangular matrices, but the easiest solution using Eigen methods would be to use the ComplexSchur decomposition (also use adjoint() instead of transpose() -- and 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;

I think X should always be real, so you can replace the last two lines by

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

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

Answer 2

@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 scipy, which explains why the OP's scipy implementation worked fine.