Solve for matrix X in XAX'=B 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.
 A: 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 
A: 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.
