Sylvester-like solver
Submitted by Chen-Pang He
Assigned to Chen-Pang He
Link to original bugzilla bug (#840)
Description
This is an announcement for whom is interested in the same topic.
I'd like to write a module that provides solvers for Sylvester-like equations.
- Lyapunov equation (CT) AX + XA* = Q
- Lyapunov equation (DT) AXA* - X = Q
- Sylvester equation (CT) AX + XB = C (1)
- Sylvester equation (DT) AXB - X = C
More equations may be added after beta.
Algorithm
A classical solution is Bartels–Stewart algorithm. Take (1) for example. We compute the Schur forms of A and B, viz. UTU* and VSV*, transforming the equation into (2).
UTUX + XVSV = C (2)
Apply U* on the left and V on the right to transform (2) into the equation below.
T(UXV) + (UXV)S = U*CV
Parentheses are added to aid reading. A general SYCT (Sylvester CT) is reduced to triangular SYCT.
However, if A and B are to be used only once, a better way is not to compute the Schur form of the larger. We stop at cheaper Hessenberg form. Assume the larger is A without losing generality. Set the Hessenberg form of A is UHU*.
H(UXV) + (UXV)S = U*CV
Now we see the trade-off: we have to solve Hessenberg instead of triangular system.
Challenges
- Blocking (different depths on-the-left and on-the-right)
- Do we need HessenbergView or PartialPivLU is enough?
- RealSchur or ComplexSchur for real matrices?
- Do we need QuasiTriangularView, or HessenbergView or even PartialPivLU is enough?
- Blocking is more tricky