Solving Lyapunov-like equation

Question

The following matrix equation might be a Lyapunov-like equation, but it seems hard for me to develop a simpler way to solve it. From the computation effort, I need some help for solving the special case of the following Lyapunov equation:

Let $X$ be an $n\times n$ symmetric matrix, and $I$ an identity matrix, and $A$ is a matrix whose entries are all between 0 and 1, and $A$ is invertible. I need to solve $X$ in the following equation:

$$AX+XA^T=I$$

Previously, I found some article discussing on using Krylov subspace to solve the following Lyapunov equation:

$$AX+XA^T=b \cdot b^T$$

where $b$ is a vector. Due to $b \cdot b^T$ being a rank-one matrix, Krylov subspace appraoch is highly efficient. Now in my case it is the identity matrix $I$, but $X$ in my case is symmetric. I found that in my equation $AX$ and $XA^T$ are symmetric. So by letting $Y=AX$, my equation can be reduced to :

$$Y+Y^T=I \quad \textrm{with } \ \ Y=AX$$.

I don't know how to continue this.

Another common way is to use tensor product to rewrite my equation as:

$$(I \otimes A + A \otimes I) vec(X) = vec(I)$$

but the LHS of the above equation is $n^2 \times n^2$ size, which is too large to solve.

Is there any other efficient way to solve this? Any advices are warmly welcome!

Answer 1

I will summarize everything, for future reference.

1. All Lyapunov equations $AX+XA^T=B$ have a unique, symmetric solution $X=X^T$, unless there is a $\lambda\in\mathbb{C}$ such that $\lambda$ and $-\lambda$ are both eigenvalues of $A$. This holds, for instance, when $A$ is Hurwitz stable, which is a common case.
2. As far as I know, reducing to $B=I$ (or to $Y+Y^T=I$ with a nontrivial constraint on $Y$) does not help in its numerical solution.
3. There are integral representations and reformulations as a system of $n^2$ linear equations in $n^2$ unknowns using Kronecker products, but they are also worthless from the point of view of numerical solution.
4. Up to $n\approx 1000$, the standard algorithm (also implemented in MATLAB's lyap) is the $O(n^3)$ Bartels-Stewart algorithm. The Hessenberg-Schur algorithm is a variation that helps for the more general Sylvester equation, but not in this case. For a general Sylvester equation $AX+XC=B$, the HS method requires one Schur decomposition of $A$ instead of one of $A$ and one of $C$, but if $C=A^T$ then the second decomposition comes for free.

5. For values of $n$ larger than $1000$, and in general for the case of $A$ large and sparse, the standard algorithms are either rational Krylov subspace methods (MATLAB code and papers on V. Simoncini's page, http://www.dm.unibo.it/~simoncin/software.html) or the ADI method (Matlab and C code on the page of Peter Benner's group, http://www.mpi-magdeburg.mpg.de/csc). Both can handle problems where $B$ has low rank.

6. The case in which $A$ is large and sparse and $B$ has full rank (for instance $B=I$) is more difficult; solving it efficiently is still an active research problem. This preprint seems to go in the right direction http://www.math.cts.nthu.edu.tw/download.php?filename=683_a93bf34d.pdf&dir=publish&title=prep2012-05-005, but the methods is much newer and I have never tried it. EDIT: I should rather recommend this approach by a collaborator of mine https://arxiv.org/abs/1711.05493 ; it should solve exactly this problem for the case of large matrices.

Answer 2

As Federico Poloni pointed out, the Hessenberg-Schur algorithm, used by MATLAB's lyap.m function is a much better choice. It is a refined version of the older Bartels-Stewart algorithm (which also works pretty well). Here's the original paper for Hessenberg-Schur, by Golub-Nash-van Loan:

https://www.cs.cornell.edu/cv/ResearchPDF/Hessenberg.Schur.Method.pdf

In your case, since $A=B^T$, things are even a bit simpler, in that only one matrix needs to be decomposed.

Answer 3

Let $A$ be an invertible $n\times n$ matrix which is antisymmetric: $A^T=-A$, e.g. the symplectic matrix $$ \begin{pmatrix} 0&1 \\\\ -1&0 \end{pmatrix}. $$ The equation $AX+XA^T=I$ cannot have a matrix solution $X$ since that would imply $$ AX-XA=I, $$ which is impossible since $trace (AX-XA)=0$. Of course that does not contradict the previous answer, but shows that some further conditions should be imposed on $A$.

Answer 4

If the $n\times n$ matrix $A$ is negatively stable (i.e. $\mbox{Re} \; \lambda_i <0$ for all $i=1,...,n$ where $\lambda_i$ are the eigenvalues of $A$), then for any $n\times n$ matrix $C$ there exists a unique $X$ such that

$$ AX+XA^T = C.$$

See Theorem 6.4.2 of Ortega ("Matrix Theory: A Second Course" 1987).

Answer 5

I hope the answer below is somewhat helpful.

Let me first summarize some basic facts.

It is known that the equation

\begin{equation*} AX + XA^T = B, \end{equation*}

has a unique solution if the matrix $A$ is positively stable (i.e., has spectrum in the right half plane). If $A$ is diagonal with entries $a_1,\ldots,a_n$, then the solution to the equation can be given in closed form

\begin{equation*} X = D \circ B, \end{equation*} where [EDIT:] $D$ is a matrix with entries $1/(\bar{a}_i+a_j)$.

In the more general case, for positively stable $A$, the solution to the above equation can be represented as

\begin{equation*} X = \int_0^\infty e^{-tA}B(e^{-tA})^Tdt \end{equation*}

But that does not seem to be computationally that nice.

If $n$ is largish, one can still solve the linear system written using tensor products by using an iterative algorithm for solving the linear system, as long as the iterative algorithm (e.g., conjugate gradient, or other related methods) depends on just "matrix-vector" products. Because you would need to only compute $(A \otimes I + I \otimes A)x$ several times, and that can be done using matrix multiply without actually forming the tensor products.