What's the solution to $X = AXA$? I have given a general $n \times d$ matrix $B$, and need to compute two similarity matrices $R \in \mathbb{R}^{n \times n}$ and $S \in \mathbb{R}^{d \times d}$ as follows:
$R = BSB^\top$ and $S = B^\top RB$.
Inserting one equation into the other yields $S = B^\top B S B^\top B$. I assume a solution doesn't always exist. In this case, I'm interested in an approximate solution that minimizes $||S - B^\top B S B^\top B||_F$.
Is there a name for this kind of problem? And is there an efficient way to solve it (possibly using some iterative algorithm)?
Note: The zero matrix would satisfy above equation, but I need similarity matrices, i.e. symmetric matrices for which the diagonal elements must be equal to 1.
 A: Given $A$ is $d \times d$ and real-valued:
The formula $X = AXA$ can be rewritten with $X$ reordered as $\operatorname{vec}(X)$, a column vector of $d^2$ elements:
$$
\operatorname{vec}(X) = \left( A^T \otimes A \right) \operatorname{vec}(X)
$$
... where $\otimes$ is the Kronecker product,  a $d^2 \times d^2$ matrix.
So, there are non-zero solutions iff $\left(A^T \otimes A\right)$ has any eigenvalues equal to 1, and the solution(s) are the corresponding eigenvectors.
The $d^2$ eigenvalues of  $\left( A^T \otimes A \right)$ can be obtained by multiplying all pairings of the $d$ eigenvalues of $A$. So if any of the eigenvalues of $A$ have an absolute value of 1 (including complex values), or if any pair of them satisfies $\lambda_j \lambda_k =1$, you have a solution to $X=AXA$.
[In your case $A=B^TB$ is symmetric so all the eigenvalues of $A$ will be real.]
Of course, it may not be possible to scale it to that $X$ is a similarity matrix. Clearly, all solutions $X$ will be singular, unless the determinant of $A$ is 1 or -1. In fact, any $X$ obtained directly from an eigenvector of the Kronecker product will be of rank 1, being the outer product of eigenvectors of $A$ (this is discussed more below), so a non-singular $X$ solution could be obtained only via a sum of multiple such solutions.

For the case where $A$ is symmetric:
let $P=\left( A^T \otimes A \right)$ 


*

*if $X$ is a solution then so is $X^T$

*For every eigenvalue of $P$ which is 1, the corresponding eigenvector gives a solution for $X$.

*If there is more than one solution, you can arbitrarily combine them linearly to get others.

*If there is a unit eigenvalue of $P$ where the corresponding eigenvector gives $X$ which is not symmetric, there will be 
another one giving $X^T$ so you can add these and get a symmetric solution. These are the unit eigenvalues of $P$ which are products of two eigenvalues of $A$ and thus occur in pairs. Those which arise from squares of {-1,1} eigenvalues of $A$ give symmetric solutions directly.


So, the more unit eigenvalues of $P$ you have, the more symmetric solutions you have; if you have enough of these you can solve for how to combine them to get a result with 1's along the diagonal.
For instance, suppose $d=5$ and the eigenvalues of $A$ are $[ 1,1, {1 \over 2} , 2,2]$ you'll have 8 unit eigenvalues of $P$ - two from squaring these, 3+3 from multiplying them in pairs; and you'll thus get 5 symmetric results for $X$. I'm not sure that all 5 will be independent on their diagonals, though.
And each $X$ candidate -- the eigenvector of $P$ corresponding to an eigenvalue of $P$ which is the product of two given eigenvalues of $A$ -- can be found as the outer product of the two corresponding eigenvectors of $A$, reordered to a vector. Or left in matrix form since you want $X$. So you don't need to find $P$ or all of its eigenvectors.
In a practical application, you may need to do a fitting to the best symmetrical solutions available, including ones corresponding to eigenvalues of $P$ which are 'close enough' to 1.

Or, if you need to adjust the eigenvalues of $A$ a bit to get an exact solution for $X$, it's relatively simple to reflect those changes back to $B$, i.e. find a new $B$ close to the original but yielding an exact $X$ solution. The new $B$ will have the same singular value decomposition as the original but with slightly different SV's:
(1) find eigen decomposition of $A$ :
$$ B^T B = A = Q L Q^T
$$
... where $L$ is a diagonal matrix with the eigenvalues of $A$; and $Q$ is an orthonormal matrix with the eigenvectors in columns
(2) Find $X$ solution(s) as above, adjusting eigenvalues as needed; This leads to a new set of eigenvalues $L'$ and new value for $A$, call it $A'=Q L' Q^T$
(3) find a diagonal matrix $R$, whose diagonal elements are the square roots of the factors by which you multiplied the eigenvalues. So:
$$\begin{align}  
L' &= RLR \\ 
A' &= Q R L R^T Q^T 
\end{align}$$
(4) find new value $B' = BQRQ^T$.
Now
$$\begin{align}
 {B'}^T B' &= QR Q^T B^T B Q R Q^T\\
        &= QR Q^T A Q R^T Q^T \\
        &= QR Q^T Q L Q^T Q R^T Q^T \\
        &= QR L R^T Q^T \\
        &= A'
\end{align}
$$
So, assuming $R$ is close to an identity matrix, $B'$ is close to the original $B$ and gives the exact $X$ solution(s) selected by adjusting eigenvalues at (2)
