Implementing the Bartels–Stewart algorithm in Eigen3?

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

Solution

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