Solving "transposed" Lyapunov equation AX+X'B=Q
5 views (last 30 days)
Show older comments
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
2 Comments
Answers (2)
Matt J
on 31 Oct 2020
Edited: Matt J
on 31 Oct 2020
I don't know how large your system is, but if not super-large you could use func2mat,
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]);
2 Comments
John D'Errico
on 31 Oct 2020
Edited: John D'Errico
on 31 Oct 2020
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.
Bruno Luong
on 31 Oct 2020
Edited: Bruno Luong
on 31 Oct 2020
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.
0 Comments
See Also
Categories
Find more on Matrix Computations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!