# Least squares solution to underdetermined Lyapunov equation

I need to solve an underdetermined Lyapunov equation for unknown $$n\times n$$ matrix $$X$$.

$$AX + XA = B$$

The naive method is to vectorize $$x=\operatorname{vec}(X)$$ and use a least squares solver on the following equation to find the least-squares solution.

$$(I\otimes A + A\otimes I)x = \operatorname{vec}B$$

Experimentally I found that when $$A$$ and $$B$$ are positive semidefinite, I get the same solution by expressing $$X$$ in terms of $$U,s$$, the eigenvectors and eigenvalues of $$A$$ as below (proof) and modifying pointwise (Hadamard) division to skip division by zero, like how pseudo-inverse implementations do it.

$$X=U \left( \frac{U' BU}{s + s'} \right) U'$$

Does this appear in the literature, or does someone see a way to prove that this recovers least-squares solution?

Example

$$\text{A=}\left( \begin{array}{ccc} 8 & -8 & -8 \\ -8 & 9 & 8 \\ -8 & 8 & 8 \\ \end{array} \right)$$

$$\text{B=}\left( \begin{array}{ccc} 5 & 5 & -5 \\ 5 & 9 & -3 \\ -5 & -3 & 6 \\ \end{array} \right)$$

This equation is underdetermined because $$A$$ and $$B$$ are singular, hence standard Lyapunov solver fails. However, both least squares and truncated spectral decomposition methods succeed with the same answer:

$$X=\frac{1}{640}\left( \begin{array}{ccc} 1789 & 2928 & -1329 \\ 2928 & 4672 & -1968 \\ -1329 & -1968 & 869 \\ \end{array} \right)$$

Notebook

• Since $U\ {\rm and}\ s$ both depend on $X$, the new equation is implicit. How do you calculate the RHS, i.e. eigenvectors/eigenvalues of an unknown matrix?
– greg
Commented Nov 16, 2023 at 20:25
• They depend on A not X, heres constructive proof this works math.stackexchange.com/a/3355045/998 Commented Nov 16, 2023 at 21:03

To be frank, this is obvious. Let $$USU^\ast$$ be a unitary diagonalisation of $$A$$ and write $$S=\operatorname{diag}(s_1,s_2,\ldots,s_n)$$. Let $$C=U^\ast BU$$ and $$Y=U^\ast XU$$. Then $$\|AX+XA-B\|_F^2 =\|SY+YS-C\|_F^2 =\sum_{i,j}|(s_i+s_j)y_{ij}-c_{ij}\big|^2.\tag{1}$$ As each $$|(s_i+s_j)y_{ij}-c_{ij}\big|$$ contains a unique element of $$Y$$, minimising the above is equivalent to minimising each residual individually. When $$s_i+s_j>0$$, the minimum value of $$|(s_i+s_j)y_{ij}-c_{ij}\big|$$ is zero and it is attained by $$y_{ij}=(s_i+s_j)^{-1}c_{ij}$$. When $$s_i+s_j=0$$, the residual is equal to the constant $$|c_{ij}|$$. Hence $$y_{ij}$$ can assume any value and the least-norm least-squares solution is $$y_{ij}=0$$. So, in summary, the least-norm solution to the least-squares problem can be written as $$y_{ij} =\begin{cases} \frac{c_{ij}}{s_i+s_j}&\text{if } s_i+s_j\ne0,\\ 0&\text{otherwise}, \end{cases}$$ i.e., $$y_{ij}=(s_i+s_j)^+c_{ij}$$ where $$(s_i+s_j)^+$$ denotes the Moore-Penrose pseudo-inverse of the number $$s_i+s_j$$. Since $$\|Y\|_F=\|UYU^\ast\|_F$$, the matrix $$X=UYU^\ast$$ constructed from $$Y$$ is also the least-norm solution to the least-squares problem of minimising $$(1)$$.

• thanks for the clear explanation Commented Nov 20, 2023 at 17:58

$$\def\o{{\tt1}} \def\LR#1{\left(#1\right)} \def\op#1{\operatorname{#1}} \def\diag#1{\op{diag}\LR{#1}} \def\Diag#1{\op{Diag}\LR{#1}} \def\trace#1{\op{Tr}\LR{#1}} \def\frob#1{\left\| #1 \right\|_F} \def\qiq{\quad\implies\quad} \def\p{\partial} \def\grad#1#2{\frac{\p #1}{\p #2}}$$Write the Lyapunov/Sylvester equation as \eqalign{ F(X) = \LR{AX + XA - B} \;\approx\; 0 \\ } The associated Least-Squares problem is
\eqalign{ \min_X\:\tfrac12\,\frob{F}^2 \\ } Calculate the differential and gradient of the LS function \eqalign{ \phi &= \tfrac12\,\LR{F:F}\\ d\phi &= F:dF \\ &= F:\LR{A\:dX+dX\:A} \\ &= \LR{AF+FA}:dX \\ \grad{\phi}{X} &= AF+FA \\ &= A^2X + 2AXA + XA^2 - AB - BA \\ } To solve the zero-gradient condition, use the technique from your previous posts and substitute the eigenvalue decomposition \eqalign{ A &= USU^T, \qquad s=\diag S,\quad S=\Diag s \\ } and use this to define some new matrix variables \eqalign{ C &= U^TBU, \qquad Y &= U^TXU, \qquad P &= S^2 = \Diag p \\ } Pre/post multiplying the zero-gradient condition by $$U^T\ {\rm and}\ U\:$$ yields \eqalign{ \LR{PY + 2SYS + YP} &= \LR{SC+CS} \\ } Now employ standard tricks involving Hadamard products, diagonal matrices, and the all-ones vector $$\o$$ \eqalign{ SMP = sp^T\odot M \\ SM = s\o^T\odot M \\ MS = \o s^T\odot M \\ } to rewrite the condition as \eqalign{ &Y\odot\LR{p\o^T+ 2ss^T +\o p^T} = C\odot\LR{s\o^T+\o s^T} \\ &Y = C\odot\LR{s\o^T+\o s^T}\oslash\LR{p\o^T+ 2ss^T +\o p^T} \\ &Y = C\oslash\LR{s\o^T+\o s^T} \\ } So the least-squares solution to the problem is indeed \eqalign{ X &= UYU^T \\\\ }

In the above, $$\{\oslash,\,\odot\}$$ denote Hadamard division and multiplication, while a colon is used to denote the matrix inner product (aka the Frobenius product) $$A:B = \trace{A^TB}$$