Anticommutator linear system

[SangalliDavide]

Dear developers,
I need to solve the following linear system:
AX+XA^dag=M
where A and M are known "n x n" matricies and X is the unknown "n x n" matrix

There are two possible strategy to solve such problem:
a) Recast it as AA vec[X] = vec[M]
where
AA = (1 x A) + (A* x 1),
i.e. an "n^2 x n^2" matrix and vec[X] and vec[M] are two "n^2 x 1" vectors

b) Diagonalize A and rewrite the problem in the rotated basis set
a_i = eigenalues of A
Y=U^dag X U new unknown matrix
M'=U^dag M U
Y_{ij}(a_i+a_j) = M_{ij}

However they are both problematic for me, because:
a) works fine, but it is too slow, needing to deal with an "n^2 x n^2" matrix
b) is not numerically sensitive, since the A matrix is of the form 1+x, with |x|<<1, and doing the diagonalization there is a lot of noise

I was thus wondering if it exist, or if it would be possible to have, a lapack subroutine which solves the commutator linear system of equations in a faster way than a) without needing to diagonalize (or invert) the A (or A^dag) matrix.

Kind regards,
D.

[Julien Langou]

Hi Davide,

I think you want to solve a continuous Lyapunov equation. Please google this and you will find methods and software.

We have a `rough` Sylvester equation solver in LAPACK with DTRSYL, but the matrix needs to be first pre-processed with DHSEQR. This would work for you, but I think you want to go with a continuous Lyapunov equation solver.

Please have a look. (Not my expertise.)

Cheers,
Julien.

Note: I think you mean `not numerical sensible` as opposed to `not numerical sensitive` in your initial post.

[sven]

SLICOT (http://slicot.org/) has Lyapunov solvers and utilises LAPACK.

Best wishes,

Sven Hammarling.

[SangalliDavide]

Dear Julien and Sven thank you very much for your answers.

These lapack subroutine seems to be ok for what I need (I have complex matricies):
- http://www.netlib.org/lapack/explore-ht ... 6b54427981
- http://www.netlib.org/lapack/explore-ht ... 2073dc7ec2
I'll have a look to the needed pre-processing.

From the SLICOT library I found this:
http://slicot.org/objects/software/shar ... B03TD.html
However I would need a version which deals with complex matricies .. ?
Also my "M" (C in the language of the subroutine) is anti-hermitian, not symmetric/hermitian.
Does it exist a more general version.

Best,
D.

P.S.: ok, I guess I meant "numerical sensible" ... the point is my matricies have big (b1,b2) and small (s1,s2) numbers, and I need to avoid operations like (s2=b1-b2+s1) otherwise I get a lot of noise in the results. Recasting the problem as a Matrix * vector operation worked find in this sense.