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


Thanks.


**EDIT:**


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



Non-linear matrix equation
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](https://en.wikipedia.org/wiki/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](https://en.wikipedia.org/wiki/Generalized_singular_value_decomposition) 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).


