在Eigen3中实现Bartels–Stewart算法? [英] Implementing the Bartels–Stewart algorithm in Eigen3?

查看:76
本文介绍了在Eigen3中实现Bartels–Stewart算法?的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!

问题描述

过去当我需要解 Sylvester 方程时,AX + XB = C,我使用了 scipy 的函数,solve_sylvester [1],显然是通过使用Bartels-Stewart算法将事物分解为上三角形式,然后使用 lapack 求解方程来实现的.

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.

我现在需要使用 eigen 求解方程. eigen 提供了一个函数 matrix_function_solve_triangular_sylvester [2],在文档中看起来与 scipy lapack 函数相似.代码>调用.我试图在 eigen3 中准确地翻译 scipy 的实现,但是最后我对 X 的值不满足该等式.这是我的实现:

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] 推荐答案

@chtz是正确的;这是由于Schurr分解矩阵是准三角形而不是三角形.您使用的特征求解器无法处理此类矩阵.但是,chtz是错误的,因为有些西尔维斯特(Sylvester)求解器可以处理准三角求解器.例如,lapack的 trsyl 可以交易与准三角矩阵.这就是 scipy 所称的,它解释了OP的scipy实现工作正常的原因.

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

这篇关于在Eigen3中实现Bartels–Stewart算法?的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!

查看全文
登录 关闭
扫码关注1秒登录
发送“验证码”获取 | 15天全站免登陆