How to solve Lyapunov equation via Matlab/Octave

Question

Let's say that I have my Lyapunov equation for finding the controllbility gramian $P$.

The equation is: $$APA^T - P + BB^T = 0$$

The matrecies $A, B $ are known.

Question:

Can I use the command fzero or fsolve to find $P $ ?

Answer 1

Since $P$ only appears linearly in the equation you can solve it as a system of linear equations. But Matlab and Octave are only good at solving this when the unknown is in vector form. By reshaping $P$, and eventually the rest of the equation as well, by stacking each column of $P$ on top of each other it is possible to transform the initial equation into a form that is easily solvable by Octave and Matlab.

So if the elements of the matrix $P$ can be defined as

$$ P = \begin{bmatrix} p_{11} & p_{12} & \cdots & p_{1n} \\ p_{21} & p_{22} & \cdots & p_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ p_{n1} & p_{n2} & \cdots & p_{nn} \end{bmatrix} \tag{1} $$

then the vector form of $P$ will defined as

$$ \vec{p} = \begin{bmatrix} p_{11} & p_{21} & \cdots & p_{n1} & p_{12} & p_{22} & \cdots & p_{n2} & \cdots \cdots & p_{1n} & p_{2n} & \cdots & p_{nn} \end{bmatrix}^\top. \tag{2} $$

When starting with a matrix product $X\,P$ or $P\,Y$, with all matrices in $\mathbb{R}^{n \times n}$, then their results can also be converted into a vector. The vector form of the product can also be expressed as $\tilde{X}\,\vec{p}$ and $\tilde{Y}\,\vec{p}$ respectively. It can be shown that $\tilde{X}$ and $\tilde{Y}$ are both in $\mathbb{R}^{n^2 \times n^2}$ and can be expressed as follows

$$ \tilde{X} = \begin{bmatrix} X & 0 & \cdots & 0 \\ 0 & X & \ddots & \vdots \\ \vdots & \ddots & \ddots & 0 \\ 0 & \cdots & 0 & X \end{bmatrix} \tag{3} $$

$$ \tilde{Y} = \begin{bmatrix} y_{11}\,I & y_{21}\,I & \cdots & y_{n1}\,I \\ y_{12}\,I & y_{22}\,I & \cdots & y_{n2}\,I \\ \vdots & \vdots & \ddots & \vdots \\ y_{1n}\,I & y_{2n}\,I & \cdots & y_{nn}\,I \end{bmatrix} \tag{4} $$

where $0$ and $I$ are the zero- and identity-matrix in $\mathbb{R}^{n \times n}$. In the case of $X\,P\,Y$ it can be shown that it is equivalent to $\tilde{X}\,\tilde{Y}\,\vec{p}$ or $\tilde{Y}\,\tilde{X}\,\vec{p}$, since $\tilde{X}$ and $\tilde{Y}$ commute.

In the case of the discrete Lyapunov equation, which is the slightly more general case of your infinite time controllability Gramian

$$ A\,P\,A^\top - P + Q = 0 \tag{5} $$

the matrices $X$ and $Y$ are the each others transpose. Using this property then $\tilde{X}\,\tilde{Y}$ can also be written as follows

$$ \tilde{A}\,\tilde{A^\top} = \begin{bmatrix} a_{11}\,A & a_{12}\,A & \cdots & a_{1n}\,A \\ a_{21}\,A & a_{22}\,A & \cdots & a_{2n}\,A \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1}\,A & a_{n2}\,A & \cdots & a_{nn}\,A \end{bmatrix} \tag{6} $$

Therefore the discrete Lyapunov equation from $(5)$ can also be written as

$$ \left(\tilde{A}\,\tilde{A^\top} - I_{n^2 \times n^2}\right)\,\vec{p} = - \vec{q} \tag{7} $$

where $\vec{q}$ is derived from $Q$ in the same way as $\vec{p}$ in $(2)$ from $(1)$. The solution to this equation is simply

$$ \vec{p} = -\left(\tilde{A}\,\tilde{A^\top} - I_{n^2 \times n^2}\right)^{-1}\,\vec{q} \tag{8} $$

In order to obtain $P$ you just have to reshape $\vec{p}$ into a square matrix again. One example of an implementation of this might look like:

function P = discrete_lyapunov_equation(A, Q)
% This function solves A * P * A' - P + Q = 0

n = length(A);
M = zeros(n^2);
for k = 1 : n
    for l = 1 : n
        M((1:n)+n*(k-1),(1:n)+n*(l-1)) = A * A(k,l);
    end
end
P = reshape(-(M - eye(n^2)) \ Q(:), n, n);

end

When comparing this against the Matlab function dlyap() it usually achieves the same accuracy or sometimes a little better for small $n$ (the accuracy is measured my taking the 2-norm of the left hand side of $(5)$ using the $P$ value obtained from each function). However my function is slower to calculate when $n$ is big. Namely my function is around 5 times slower when $n = 10$, but is around 17 times slower when $n = 20$. For more details look at the first figure below. Above $n = 10$ Matlab sometimes also complains when running my function that the matrix M is close the singular. However when comparing the accuracy of my function to that of dlyap() then mine is much better above $n = 13$ as can be seen in the second figure below. I think I would prefer accuracy over computation time, especially if you do not have any control over how accurate the result is going to be.

When looking into the description of dlyap() I looked up one of the references, which mentions the use of Bartel-Stewart's algorithm, which is also mentioned here which states: "...which consists of transforming A and B into Schur form by a QR algorithm, and then solving the resulting triangular system via back-substitution. This algorithm, whose computational cost is $\mathcal{O}(n^3)$ arithmetical operations". So that might also be a place to look into if you would like to know how Matlab potentially does it. However this would go way to deep into numerical methods for the scope of this question, so I will leave it as this. Also this method might give results which have really poor accuracy at large $n$.