# How can I express the solution to a discrete Lyapunov equation as an eigenvalue problem?

Given the discrete Lyapunov equation $$AXA^T - X + Q = 0$$ how can I solve for $$X$$ as a function of the eigenvectors of some matrix $$H$$?

More precisely, in the case of the continuous Lyapunov equation $$AX + XA^T + Q = 0$$ This equation can be re-written as \begin{align} AX + XA^T + Q &= 0 \\ \begin{bmatrix} X & -I\end{bmatrix} \begin{bmatrix} A^T \\ -Q - AX\end{bmatrix} &= 0 \\ \begin{bmatrix} X & -I\end{bmatrix} \begin{bmatrix} A^T & 0 \\ -Q & -A\end{bmatrix} \begin{bmatrix} I \\ X\end{bmatrix} &= 0 \\ \begin{bmatrix} X & -I\end{bmatrix} H \begin{bmatrix} I \\ X\end{bmatrix} &= 0 \tag{1} \label{eq:cont_lyap} \end{align} where $$H = \begin{bmatrix} A^T & 0 \\ -Q & -A\end{bmatrix}$$ is a $$2n \times 2n$$ Hamiltonian matrix. From \eqref{eq:cont_lyap}, we see that $$H \begin{bmatrix} I \\ X\end{bmatrix}$$ is in the null-space of $$\begin{bmatrix} X & -I\end{bmatrix}$$, which can be used to show that there exists a matrix $$W$$, such that its columns form a basis for the image of $$\begin{bmatrix} I \\ X\end{bmatrix}$$, so that $$H \begin{bmatrix} I \\ X\end{bmatrix} = \begin{bmatrix} I \\ X\end{bmatrix}W$$ By diagonalizing $$W$$, so that $$W = X_1 \Lambda X_1^{-1}$$, we get the equation $$H \begin{bmatrix} X_1 \\ X_2\end{bmatrix} = \begin{bmatrix} X_1 \\ X_2\end{bmatrix}\Lambda$$ where $$X = X_2X_1^{-1}$$. This means that we can compute $$X$$ by computing the $$2n$$ eigenvectors of $$H$$, picking out $$n$$ of them and forming the matrix $$\begin{bmatrix} X_1 \\ X_2\end{bmatrix}$$, and then computing $$X = X_2X_1^{-1}$$.

What I am looking for is a similar derivation to the one above, but for the discrete (not continuous) Lyapunov equation.

You can do it with symplectic matrix pencils instead of Hamiltonian matrices, even in the more general case of discrete-time algebraic Riccati equations. $$\begin{bmatrix} A & 0\\ -Q & I \end{bmatrix} \begin{bmatrix} I\\X \end{bmatrix} = \begin{bmatrix} I & G\\ 0 & A^T \end{bmatrix} \begin{bmatrix} I\\X \end{bmatrix}W$$ is equivalent to $$X - Q = A^TX(I+GX)^{-1}A$$, after eliminating $$W$$ from the two resulting block equations. Set $$G=0$$ to recover the case of a discrete-time Lyapunov equation.

This is all "folklore" in numerical linear algebra / control theory circles, but if you wish to have a reference you can take my review paper https://dx.doi.org/10.1002/gamm.202000018 (it's Equation 25 there).

• Thank you very much for this. Your reference in particular is brilliant. For others who don't have paid access to Wiley, the arXiv version of the paper is here. Nov 15, 2022 at 14:48
• Just one question: based on what you wrote in section 4.1 of the paper, it seems to me that the generalized eigenvalue problem in eq. 25 can only be derived as long as $A$ is invertible. You did mention that the results are still valid if $A$ is singular, but I'm wondering if it would still be possible to derive this eigenvalue problem if $A$ was singular. That is, we don't need to left-multiply both sides of eq. 25 by $$\begin{bmatrix} I & G\\0 & A^T\end{bmatrix}^{-1}$$ Nov 15, 2022 at 14:50
• @mhdadk Yes, for a suitable generalization of "eigenvalue problem". One can define and solve generalized eigenvalue problems of the form $\det(M-\lambda N) = 0$, without the need to invert $N$. In Matlab, it's eig(M, N), for instance; in Python, scipy.linalg.eig(M, N). Nov 15, 2022 at 18:28
• @mh Alternatively, you can use Cayley transforms, i.e., a change of variable of the form $\lambda = \frac{\mu + \gamma}{\mu - \gamma}$, for a suitable $\gamma > 0$; this should transform your problem into a standard eigenvalue problem (assuming a certain matrix depending on $\gamma$ is invertible) but that's trickier to handle as it changes the location of the eigenvalues wrt the unit circle, and it is less numerically stable. Nov 15, 2022 at 18:32
• Hey Federico, if you have the time, could you possibly have a look at my question Is it possible to express the solution of a matrix Riccati differential equation as an eigenvalue problem?? I think it is related to your answer here, although I could be wrong. Dec 18, 2022 at 7:52

For completeness, I'm including a full derivation for the solution of the discrete Lyapunov equation that is adapted from @FedericoPoloni's paper (specifically section 4.1). For the derivation for the more general discrete algebraic Riccati equation, see the paper.

Given the discrete Lyapunov equation $$AXA^T - X + Q = 0$$ Suppose that $$A$$ is invertible, such that $$A^{-1}$$ exists. As @FedericoPoloni mentioned in section 4.1 in their paper, this is necessary for the following derivation, although the results in section 4.1 are valid even if $$A$$ is singular.

Note that the discrete Lyapunov equation can be re-written as \begin{align} AXA^T - X + Q &= 0 \\ A^{-1}AXA^T - A^{-1}X + A^{-1}Q &= 0 \\ XA^T - A^{-1}X + A^{-1}Q &= 0 \\ \begin{bmatrix} X & -I\end{bmatrix} \begin{bmatrix} A^T \\ A^{-1}X - A^{-1}Q\end{bmatrix} &= 0 \\ \begin{bmatrix} X & -I\end{bmatrix} \begin{bmatrix} A^T & 0 \\ -A^{-1}Q & A^{-1}\end{bmatrix} \begin{bmatrix} I \\ X\end{bmatrix} &= 0 \\ \begin{bmatrix} X & -I\end{bmatrix} H \begin{bmatrix} I \\ X\end{bmatrix} &= 0 \tag{1} \label{eq:sympl_eig1} \end{align} where we defined $$H = \begin{bmatrix} A^T & 0 \\ -A^{-1}Q & A^{-1}\end{bmatrix}$$ which is a symplectic matrix that has the property that (lemma 2 in section 4.1 of @FedericoPoloni's paper), if all the eigenvalues of $$A$$ lie strictly inside the unit circle in the complex plane, then the $$2n$$ eigenvalues of $$H$$ consist of $$n$$ eigenvalues inside the unit circle, and $$n$$ eigenvalues outside the unit circle (to quickly see why this is true, note that $$H$$ is block lower-triangular, so its eigenvalues are the union of the eigenvalues of the matrices on its diagonal). This property will be useful later.

Furthermore, note from \eqref{eq:sympl_eig1} that $$H \begin{bmatrix} I \\ X\end{bmatrix}$$ is in the null space of $$\begin{bmatrix} X & -I\end{bmatrix}$$. Equivalently, because the null space of $$\begin{bmatrix} X & -I\end{bmatrix}$$ is orthogonal to the image of $$\begin{bmatrix} X & -I\end{bmatrix}^T = \begin{bmatrix} X \\ -I\end{bmatrix}$$ (see orthogonal complementarity) then $$H \begin{bmatrix} I \\ X\end{bmatrix}$$ is in the image of \begin{align} \begin{bmatrix} 0 & -I \\ I & 0\end{bmatrix} \begin{bmatrix} X \\ -I\end{bmatrix} &= \begin{bmatrix} I \\ X\end{bmatrix} \end{align} where the matrix $$\begin{bmatrix} 0 & -I \\ I & 0\end{bmatrix}$$ represents a 90-degree counter-clockwise rotation in $$2n$$-dimensional space. Let $$W$$ be an $$n \times n$$ matrix such that each of its columns span the $$n$$-dimensional image of $$\begin{bmatrix} I \\ X\end{bmatrix}$$. Then, because $$H \begin{bmatrix} I \\ X\end{bmatrix}$$ is in the image of $$\begin{bmatrix} I \\ X\end{bmatrix}$$, then $$H \begin{bmatrix} I \\ X\end{bmatrix} = \begin{bmatrix} I \\ X\end{bmatrix}W$$ which almost looks like an eigenvalue problem. The problem is that $$W$$ is not necessarily diagonal (or in Jordan normal form). Before we solve this problem, we can first determine what $$W$$ corresponds to. Recall that \begin{align} H \begin{bmatrix} I \\ X\end{bmatrix} &= \begin{bmatrix} A^T & 0 \\ -A^{-1}Q & A^{-1}\end{bmatrix} \begin{bmatrix} I \\ X\end{bmatrix} \\ &= \begin{bmatrix} A^T \\ A^{-1}X - A^{-1}Q\end{bmatrix} \end{align} Then, \begin{align} H \begin{bmatrix} I \\ X\end{bmatrix} &= \begin{bmatrix} I \\ X\end{bmatrix}W \\ \begin{bmatrix} A^T \\ A^{-1}X - A^{-1}Q\end{bmatrix} &= \begin{bmatrix} W \\ XW\end{bmatrix} \end{align} which implies that $$W = A^T$$.

Going back to our problem where $$W$$ (or $$A^T$$) is not necessarily diagonal, suppose that $$W$$ is diagonalizable, such that $$W = X_1 \Lambda X_1^{-1}$$, where $$\Lambda$$ is diagonal and contains the eigenvalues of $$A^T$$ (which are the same as the eigenvalues of $$A$$). Then, \begin{align} H \begin{bmatrix} I \\ X\end{bmatrix} &= \begin{bmatrix} I \\ X\end{bmatrix}W \\ H \begin{bmatrix} I \\ X\end{bmatrix} &= \begin{bmatrix} I \\ X\end{bmatrix}X_1 \Lambda X_1^{-1} \\ H \begin{bmatrix} I \\ X\end{bmatrix} &= \begin{bmatrix} X_1 \\ XX_1\end{bmatrix} \Lambda X_1^{-1} \\ H \begin{bmatrix} X_1 \\ XX_1\end{bmatrix} &= \begin{bmatrix} X_1 \\ XX_1\end{bmatrix} \Lambda \end{align} Moreover, let $$X_2 = XX_1$$, such that $$X = X_2X_1^{-1}$$. Then, \begin{align} H \begin{bmatrix} X_1 \\ X_2\end{bmatrix} &= \begin{bmatrix} X_1 \\ X_2\end{bmatrix} \Lambda \end{align} which represents an eigenvalue problem. Note that the $$n$$ eigenvalues of $$W$$ (or $$A^T$$) are a subset of the $$2n$$ eigenvalues of $$H$$. Also, note that $$\begin{bmatrix} X_1 \\ X_2\end{bmatrix}$$ has $$2n$$ rows and $$n$$ columns, while $$H$$ has $$2n$$ eigenvalues and eigenvectors. This implies that we will have to pick out $$n$$ eigenvalues and eigenvectors out of the $$2n$$ possible ones of $$H$$. Because $$W = A^T$$ and $$\Lambda$$ consists of its eigenvalues, and because $$A$$ is already given and assumed to be stable, then we should pick the $$n$$ eigenvectors that correspond to the $$n$$ eigenvalues in $$\Lambda$$ to form the matrix $$\begin{bmatrix} X_1 \\ X_2\end{bmatrix}$$. Then, we can compute $$X$$ as $$X = X_2X_1^{-1}$$.

• Maybe you should switch the accepted answer to yours; there is much more detail here. Nov 16, 2022 at 7:58
• @FedericoPoloni this answer would not exist without your help :-) Nov 16, 2022 at 11:00
• Note that if you only care about getting $H\begin{bmatrix}I\\X\end{bmatrix} = \begin{bmatrix}I\\X\end{bmatrix} W$ (or the equivalent equation in my answer without inverting the 2x2 block matrix in the RHS), then this relation is easier to prove directly: just write down the equations corresponding to the two blocks, and check that they are verified if $X$ solves the discrete Lyapunov equation and $W=A^T$. Nov 16, 2022 at 18:22
• @FedericoPoloni not sure I understand. Do you mean this? \begin{align} H \begin{bmatrix} I \\ X\end{bmatrix} &= \begin{bmatrix} I \\ X\end{bmatrix}W \\ \begin{bmatrix} A^T \\ A^{-1}X - A^{-1}Q\end{bmatrix} &= \begin{bmatrix} W \\ XW\end{bmatrix} \end{align} Nov 16, 2022 at 19:33
• Yes; from the first block you get $A^T=W$, and if you multiply the second one by $A$ you get $X-Q = AXA^T$. Since $X$ solves the Lyapunov equation this equality is verified; this proves that $H\begin{bmatrix}I\\X\end{bmatrix} = \begin{bmatrix}I\\X\end{bmatrix} W$. Nov 16, 2022 at 21:26