Solving "transposed" Lyapunov equation AX+X'B=Q

I need to solve AX+X'B=Q which is almost like Lyapunov equation, except X is transposed in one of the terms.

• Is there a name for this equation?
• Is there a way to massage lyap2 or some other function to work for this?

I like diagonalization-based approach of lyap2 because it works for nearly singular A,B (can modify lyap2 to ignore near zero-eigenvalues)

Edit: I did some digging to find that it's called "T-Sylvester" or "Sylvester-transpose" equation

[Matt J]

And Q is symmetric?

[Yaroslav Bulatov]

Q is not symmetric. This equation comes up when trying to invert operation Q=F(X) written in einsum notation as Q_{ij}=T_{ijkl}X_{kl} where T factors as T_{ijkl}=A_{ik}B_{jl}+C_{jk}D_{il}

Answer 1

[Matt J]

I don't know how large your system is, but if not super-large you could use func2mat,

https://www.mathworks.com/matlabcentral/fileexchange/44669-func2mat-convert-linear-function-to-matrix

to recast it as a regular linear system,

[nx,mx]=size(A); 

X=reshape( func2mat( @(X) A*X+X'*B, ones(mx,nx)    )\Q(:),  [mx,nx]);

[Yaroslav Bulatov]

My n-by-n matrices are 1024-by-1024 so I'm kind of limited to O(n^3) algorithms

[John D'Errico]

It is always important to tell specifics like the size of the problem. But really, your problem is not one of size 1024, since there really are 1048576 unknowns.

Answer 2

[Bruno Luong]

Sove to real equation.

It might take a couple of minutes for matrix of size 1024, but this code using iterative solver seems able to give back a solution.

Depending on A and B, this problem might be ill-posed, so you might run into all sort of issue of solving ill-posed linear system. But this is another topic.

% A = some n x n real matrix
% B = some n x n real matrix
% Q is A*X + X'*B real matrix

% Solve for X
maxiter = 10000; % increase if needed
x = lsqr(@(x,flag) Lyap(x,flag,A,B), Q(:), 1e-6, maxiter);
sz = size(A);
X = reshape(x,sz);

% Check residual
norm(A*X+X'*B-Q,'fro') / norm(Q,'fro')

%
function y = Lyap(x,flag,A,B)
n = sqrt(length(x));
sz = [n,n];
X = reshape(x,sz);
if strcmp(flag,'notransp') 
    Q = A*X+X'*B;
elseif strcmp(flag,'transp')
    Q = A'*X+B*X';
end
y = Q(:);
end

I believe the equation is no-longer linear for complex case.