Solve for matrix X in XAX'=B

Question

I have a problem of the form $XAX^\top=B$, where $A$ and $B$ are symmetric matrices and $X$ need not be symmetric. I'd like to solve for $X$, e.g. by expressing the problem as $C\text{vec}(X)=D$ (where vec is the vectorization operation) and solving by least squares. I'm looking for ways to express the problem in a solvable way. I'm familiar with the vec trick but that just seems to make things worse in this case.

Answer 1

There is an algorithmic way to solve $P^T AP = D_1$ which also handles the inverse $Q = P^{-1}$ and $Q^T D_1 Q = A.$ With a few extra steps, you may demand that the diagonal entries of $D_1$ are in decreasing order, $d_{11} \geq d_{22} ... \geq d_{nn}$

Pause: method discussed at http://math.stackexchange.com/questions/1388421/reference-for-linear-algebra-books-that-teach-reverse-hermite-method-for-symmetr

Do the same with $B,$ with $R^T B R = D_2,$ $S^T D_2 S = B.$ Again, force the diagonal entries in order.

IIIFFFF the counts of positive elements, zero elements, and negativ elements agree for $D_1$ and $D_2,$ you can make a diagonal matrix $M$ with some square root ( of the positive entry ratios) entries, once again nonsingular, and $M D_1M = D_2.$ Once again, $M^{-1}$ will be obvious, $M$ is diagonal with strictly positive elements... If such $M$ is not possible, by Sylvester's Law of Inertia, you are out of luck

Answer 2

As was clarified in the comments, we know that $A$ and $B$ are positive semidefinite. Let us further assume that they are of the same size.

A solution will exist to the equation if and only if $\operatorname{rank}(A) \geq \operatorname{rank}(B)$.

In this case, let $p$ be the rank of $A$ and let $q$ be the rank of $B$. Begin by finding $P,Q$ such that $A = PP^T$ and $B = QQ^T$, where $P$ has $p$ columns and $Q$ has $q$ columns. These can be found, for instance, with a Cholesky decomposition.

Rewrite the equation as $$ XAX^T = B \iff (XP)(XP)^T = QQ^T. $$ It now suffices to find any $X$ solving $XP = QU$, where $U$ is any choice of orthogonal matrix. This equation can be solved using your "vec trick", or you can separately solve the systems determined by each row of $X$.

I believe that every solution $X$ solves the above equation for some choice of $U$, but I don't have a proof off the top of my head.

---

I'll focus on the case where $A$ and $B$ are invertible and of the same size for now, but this method can be extended.

First of all, we note that your equation will have a solution if and only if $A$ and $B$ have equal signatures. Under this assumption, suppose that $p$ is the number of positive eigenvalues and $n$ is the number of negative eigenvalues associated with both $A$ and $B$.

There exist invertible matrices $P$ and $Q$ such that $A = PDP^T$ and $B = QDQ^T$ where $$ D = \pmatrix{I_p&0\\0&-I_n}. $$ Now, rewrite $$ XAX^T = B \iff\\ XP D P^TX^T = Q D Q^T \iff\\ (Q^{-1}XP) D (Q^{-1}XP)^T = D. $$ Now, let $Y = Q^{-1}XP$; it suffices to solve the system $YDY^T = D$. Just like $X$, $Y$ must be square with size $n + p$. Break $Y$ into blocks, i.e. take $$ Y = \pmatrix{Y_1 & Y_2} $$ where $Y_{1}$ has $p$ columns. By block-matrix multiplication, we rewrite the equation as $$ Y_1Y_1^T - Y_2Y_2^T = D. $$ It suffices to consider the solutions to this equation. Note that $Y_1Y_1^T$ must be positive semidefinite with rank at most $p$, and $Y_2Y_2^T$ is positive semidefinite with rank at most $n$. Note that corresponding to any solution we will have a decomposition $D$ into $D = D_+ - D_-$, where both $D_+$ and $D_-$ are positive semidefinite matrices where $Y_1Y_1^T = D_+$ and $Y_2Y_2^T = D_-$.

I believe that the only such decomposition where $D_+$ has rank at most $p$ and $D_-$ rank at most $n$ comes from $$ D_+ = \pmatrix{I_n&0\\0&0}, \quad D_- = \pmatrix{0&0\\0&I_p}. $$

Corresponding this case, we have the solutions have the form $$ Y_1 = \pmatrix{U_1\\0}, \quad Y_2 = \pmatrix{0\\U_2} $$ where $U_1$ and $U_2$ are orthogonal matrices of sizes $p$ and $n$ respectively.