Non-linear matrix equation 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).
 A: 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).
A: Comparing equation (1) with its transpose, we see that $AB X^T = X B^T A$.
I would start by solving this linear equation.
A: 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.
