I want to solve the following non linear matrix equation for $X\in\mathbb{R}^{N\times N}$:

\begin{equation} XX^{\top}+ABX^{\top}-A=0 \qquad (1) \end{equation}

For a given matrices $A\in\mathbb{R}^{N\times N}, B\in\mathbb{R}^{N\times N}$ and $A=A^{T}, A\succeq0$.

Is there a known stable numerical solution for (1)?

If not, assuming that the solution is symmetric, $X=X^{\top}$, we get the following matrix equation:

\begin{equation} X^{2}+ABX-A=0 \qquad (2)\\ X=X^{\top} \end{equation}

Is there a known stable numerical solution for (2)?

I know that eq.(1) and eq.(2) must have a solution (due to basic optimization considerations).



(1) We are willing to relax the stability requirement.

(2) B is a low rank matrix.

(3) My attempt: I used the necessary condition (As Robert suggested): \begin{equation} ABX^{\top}=XB^{\top}A \qquad (3) \end{equation}

to get the following equation:

\begin{equation} XX^{\top}+\frac{1}{2}ABX^{\top}+\frac{1}{2}XB^{\top}A-A=0 \qquad (4) \end{equation}

Now we can use Riccati to find a symmetric solution. However, not every solution of (4) is a solution of (1).


I'll give an explicit expression for a family of solutions to your first problem (the more general one without the symmetry constraint).

Let us use Robert Israel's suggestion as in your last edit, and start from \begin{equation} XX^{\top}+\frac{1}{2}ABX^{\top}+\frac{1}{2}XB^{\top}A-A=0 \tag{4} \end{equation}

You can factor it as \begin{equation} (X+\frac12 AB)(X^\top +\frac12 B^\top A) = A + \frac14 AB B^\top A. \tag{*} \end{equation} The matrix $C = A + \frac14 AB B^\top A$ is positive semidefinite, hence it can be factored as $C=LL^\top$ (you can take $L=C^{1/2}$, for instance, or compute a Cholesky-like factorization where you stop the algorithm if you encounter a zero submatrix).

Then, a family of solutions is $X = LQ - \frac12 AB$, where $Q$ is any orthogonal matrix.

If $C$ is invertible, then this is the full set of solutions, since you can rewrite the equation as $L^{-1}(X+\frac12 AB)(X^\top +\frac12 B^\top A)L^{-\top} = I$, and hence the equality holds iff $Q=L^{-1}(X+\frac12 AB)$ is orthogonal. If $C$ is singular, I suspect there are more, and probably it is possible to find them all with some more work (change basis so that $C = \begin{bmatrix}C_{11} & 0 \\ 0 & 0\end{bmatrix}$).

But in any case, at this point the question is: do these solutions solve your problem, or do you have more constraints on $X$?

Remark: the condition $A \succeq 0$ can be relaxed: solutions exist if and only if $C\succeq 0$ (the LHS of $(*)$ is always positive semidefinite so $C$ must be too.)

EDIT: as noted by @loupblanc in his/her answer, this solution does not guarantee that $ABX^\top$ is symmetric for all choices of $Q$, hence it solves (4) but not in general (1). The missing condition is that $XB^\top A = (LQ-\frac12AB)B^\top A$ is symmetric, which is equivalent to requiring $LQB^\top A$ symmetric.

In the case in which $L$ is invertible, this is equivalent to requiring $QB^\top A L^{-\top}$ to be symmetric. A valid choice for $Q$ is given by computing the polar decomposition $L^{-1}AB=ZS$, with $Z$ orthogonal and $S$ symmetric, and taking $Q=Z$.

In the general case, of course there is an even more arcane matrix decomposition. We take a generalized SVD of the pair $(L^\top,B^\top A)$, i.e., $L^\top = U\Sigma_1Y, \, B^\top A = V\Sigma_2 Y$, where $Y$ is invertible, $\Sigma_1,\Sigma_2$ are diagonal with $\Sigma_1^2+\Sigma_2^2=I$, and $U,V$ are orthogonal (all square). Then, direct verification shows that choosing $Q=UV^\top$ gives a symmetric $XB^\top A$.

So in the end a solution to (1) exists iff $A + \frac14 ABB^\top A$ is positive semidefinite, and we can compute it in $O(n^3)$ using some lesser-known matrix decompositions (that are implemented stably, for instance, in Matlab or Scipy/Numpy).

  • $\begingroup$ write "@loup blanc", otherwise I am not advised. On the other hand, it's "he", otherwise my pseudo would be "louve blanche". $\endgroup$ – loup blanc Sep 11 '18 at 13:21
  • $\begingroup$ @loupblanc The syntax without the space is the correct one (as you can confirm if you received a notification of this comment). The reason why you were not notified is that @-notifications currently work only in comments and not in questions/answers. See meta.mathoverflow.net/q/577/1898 . $\endgroup$ – Federico Poloni Sep 11 '18 at 13:29

Comparing equation (1) with its transpose, we see that $AB X^T = X B^T A$. I would start by solving this linear equation.

  • $\begingroup$ For reviewers: this post appeared in the "Low quality posts" review queue and I was going to press "Recommend deletion" when, on a second reading, I realized that, despite its brevity, it deserves to stay. $\endgroup$ – Alex M. Sep 3 '18 at 18:37

In this post, I assume that $A$ is symmetric $>0$ (in particular invertible) and $B$ is invertible.

$\textbf{Proposition}$. i) Equation (1) generically admits $2^n$ real solutions $X$.

ii) There are equations (1) that do not admit any stable solutions (in the sense: $X$ is stable iff for every $\lambda\in spectrum(X), Re(\lambda)<0$).

$\textbf{Proof}$. i) I use the Robert's idea which is clearly a good one.

Let $S=XB^TA$ (necessarily a symmetric matrix). (1) is equivalent to

(2) $SUS+S-A=0$ where $U=A^{-1}B^{-T}B^{-1}A^{-1}$ is symmetric $>0$.

Thus (2) is a symmetric Riccati equation; GENERICALLY, it has only symmetric real solutions. Beware, (for example) if $U=A=I_n$, then there are non-symmetric solutions.

Note that (2) is equivalent to

(3) $Y^2+Y-AU=0$ (where $Y=SU$) or (3') $(Y+1/2I)^2=1/4I+AU$.

Note that $spectrum(AU)\subset \mathbb{R}^+$ and that $AU $ is diagonalizable. Generically (3') admits $2^n$ real solutions and, consequently, (2),(1) too.

ii) A (randomly chosen) example with no stable solution is as follows for $n=3$.

$A= \begin{pmatrix}2901/250& -687/100& -5479/500\\-687/100& 967/200& 3041/500\\-5479/500& 3041/500& 5507/500\end{pmatrix}$

$B= \begin{pmatrix}63& 76& 59\\89& -41& 9\\-15& 50& -24\end{pmatrix}$.

EDIT. I just read the good answer of @Federico Poloni ; yet, if I understand the OP's question (that is difficult because he changes his mind easily), we must add the so called Robert's condition: $(R)$ "$XB^TA$ is symmetric"; then the problem becomes more complicated.

Assume that $A>0$; then Federico shew that the general solution of $(4)$ (OP's notation) is $X=LQ-1/2AB$ where $L=\sqrt{A+1/4ABB^TA}$ is invertible and $Q$ is an arbitrary orthogonal matrix. Then, when we add the condition $(R)$, putting $U=B^TAL^{-1}$ (note that $rank(U)=rank(B)$), it remains to solve the non-obvious system in the unknown $Q$

$(*)$ $QU=U^TQ^T,QQ^T=I$.

The previous system generically has $2^n$ real solutions, as mentioned in the first part of my post. Yet, the result using the Federico'post is better because we assume only that $A$ is invertible and we do not suppose anything about $B$.

Of course, when $rank(B)<n$, the system may have an infinity of solutions. For example, some numerical tests give for $n=4$ and $U$ a random matrix

if $rank(U)=3$ then $16$ solutions. If $rank(U)=2$, then an infinity of solutions that depend on $1$ parameter.

EDIT 2. We can solve $(*)$ as follows

We consider the SVD of $U$: $VUW=diag(\sigma_1I_{i_1},\cdots,\sigma_{k-1}I_{i_{k-1}},0_{i_k})$ where the $\sigma_i> 0$ are distinct.

We seek $P\in O(n)$ s.t. $PVUW$ is symmetric. We obtain $P=diag(P_1,\cdots,P_{k-1},P_k)$ where $P_1,\cdots,P_{k-1}$ are orthogonal symmetries ($P_j^2=I_{i_j}$) and $P_k\in O(I_{i_k})$. The number of parameters for $P_j,j<k$ is $floor({i_j}^2/4)$ and for $P_k$ is $i_k(i_k-1)/2$.

Note that $WPVU$ is symmetric and that the set of required $Q$ is the set of $WPV$.

Remark. When the singular values of $U$ are distinct, $P=diag(\pm 1,\cdots,\pm 1)$ and we find the $2^n$ solutions of the generic case.

  • $\begingroup$ Thank you for your answer, unfortunately B is a low rank matrix. $\endgroup$ – Chen Zeno Sep 4 '18 at 15:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.