Реализация алгоритма Бартельса – Стюарта в Eigen3? - PullRequest
Новая
       42

Реализация алгоритма Бартельса – Стюарта в Eigen3?

1 голос
/ 08 июля 2019

В прошлом, когда мне нужно было решить уравнение Сильвестра, AX + XB = C, я использовал функцию scipy, solve_sylvester [1], которая, очевидно, работает с использованием алгоритма Бартельс-Стюарта дляполучить вещи в верхней треугольной форме, а затем решить уравнение, используя 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

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

Ответы [ 2 ]

2 голосов
/ 15 июля 2019

Ваши матрицы A и 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.

Я не знаю, есть ли решатель Сильвестра, который также обрабатывает квазитреугольные матрицы, но самое простое решениеиспользование методов 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;
1 голос
/ 17 июля 2019

@ чтц правильный;это связано с тем, что матрица разложения Шурра является квазитреугольной, а не треугольной.Собственный решатель, который вы используете, не может иметь дело с такими матрицами.Однако, chtz не прав в том, что есть решатели Сильвестра, которые могут иметь дело с квазитреугольными решателями.Например, trsyl Лапака может работать с квазитреугольными матрицами.Это то, что называется scipy, что объясняет, почему реализация scipy OP работала нормально.

...