Least squares solution to underdetermined Lyapunov equation

Question

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

Answer 1

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)$.

Answer 2

$ \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} $$