在 Eigen3 中实现 Bartels-Stewart 算法 -- 仅实矩阵?

基于这个问题和解决方案 -- 在 Eigen3 中实现 Bartels-Stewart 算法？ -- 我正在尝试使用特征库求解李雅普诺夫方程 (AX + XA^T = C(，但仅限于实矩阵。

下面的 R(使用 c++(代码有效，但涉及复数。它绝对可以简化(因为在这个框架中，没有 B 矩阵(，但主要困难是对复数的依赖。在这种情况下，真正的 schur 形式似乎是标准的替代方案，但特征函数matrix_function_solve_triangular_sylvester不起作用，因为输入矩阵不是上三角形，而是上块三角形。我很高兴看到a(消除对复数的需求的建议，然后如果可能的话，b(任何效率改进。

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