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

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

5 голосов
/ 26 мая 2020

На основе этого вопроса и решения - Реализация алгоритма Бартельса – Стюарта в Eigen3? - Я пытаюсь решить уравнения Ляпунова (AX + XA ^ T = C) с помощью библиотеки Eigen , но ограничусь реальными матрицами.

Приведенный ниже код R (с C ++) работает, но включает комплексные числа. Его определенно можно упростить (поскольку в этом кадре нет B-матрицы), но главная трудность заключается в использовании комплексных чисел. Реальная форма Шура кажется стандартной альтернативой в этом случае, но функция Eigen matrix_function_solve_triangular_sylvester не работает, потому что входная матрица - это не верхний треугольник angular, а верхний блок tri angular. Я был бы рад увидеть предложения по а) устранению необходимости в комплексных числах, а затем, если это возможно, б) любых улучшениях эффективности. 

library(expm)
library(Rcpp)
library(RcppEigen)
library(inline)
# R -----------------------------------------------------------------------
d<-6 #dimensions 
A<-matrix(rnorm(d^2),d,d) #continuous time transition
G <- matrix(rnorm(d^2),d,d)
C<-G %*% t(G) #continuous time pos def error
AHATCH<-A %x% diag(d) + diag(d) %x% A
Xtrue<-matrix(-solve(AHATCH,c(C)), d) #asymptotic error from continuous time


# c++ in R ---------------------------------------------------------------------

sylcpp <- '

    using Eigen::Map;
    using Eigen::MatrixXd;

  // Map the double matrix A from Ar
   const Map<MatrixXd> A(as<Map<MatrixXd> >(Ar));

   // Map the double matrix Q from Qr
   const Map<MatrixXd> Q(as<Map<MatrixXd> >(Qr));


  Eigen::MatrixXd B = A.transpose();

  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::MatrixXd X = ((U * Y) * V.adjoint()).real();

  return wrap(X);
'



syl <- cxxfunction(signature(Ar = "matrix",Qr='matrix'), sylcpp, plugin = "RcppEigen")

X=syl(A,-C)
X-Xtrue #approx zero
...