Given two origin-centered ellipsoids $E_0$ and $E_1$ in $\mathbb{R}^n$, I'd like to find an SPD (symmetric positive definite) transformation matrix $M$ that transforms $E_0$ into $E_1$.
Let's say $E_0$ and $E_1$ are specified by SPD matrices that take the unit sphere to the respective ellipsoid: $$ M_0 = R_0 D_0 {R_0}^{-1} \mathrm{\ takes\ unit\ sphere\ to\ }E_0\\ M_1 = R_1 D_1 {R_1}^{-1} \mathrm{\ takes\ unit\ sphere\ to\ }E_1 $$ where each $R_i$ is a rotation matrix and $D_i$ is a positive diagonal matrix whose diagonal entries are the respective ellipsoid's principal radii.
Then of course the matrix $M_1 {M_0}^{-1}$ will take $E_0$ to $E_1$, but it is in general not a symmetric matrix (even though each $M_i$ is symmetric), so that's not a solution.
I strongly suspect there is a unique SPD matrix that takes $E_0$ to $E_1$. How can it be computed?
More generally, if $R$ is any rotation matrix, then $M_1 R {M_0}^{-1}$ will take $E_0$ to $E_1$. In fact I suspect the matrices that take $E_0$ to $E_1$ are precisely the matrices of this form. So perhaps the question boils down to "Given SPD $M_0$, $M_1$, find a rotation $R$ such that $M_1 R {M_0}^{-1}$ is symmetric".
This comes from a statistics question: given a point cloud in $\mathbb{R}^n$ and a target covariance matrix, I want to find an SPD transformation such that the transformed point cloud has the target covariance. Intuitively, it seems like what's needed is an SPD transformation that transforms the point cloud's error ellipsoid (found by taking the square root of its covariance matrix) to the error ellipsoid of the target covariance matrix; so that's why I'm asking this question. However, even if I found such a transform, I'm not certain the transformed point cloud will have exactly the desired target covariance; if it turns out that it doesn't, I'll ask a different question for the underlying statistics problem.
Denote by $B$ the unit sphere. As $QB=B$ for every real orthogonal matrix $Q$, if we can solve the equation $$PM_0Q=M_1\tag{1}$$ for some positive definite $P$ and some real orthogonal $Q$, then $P$ will take $E_0$ to $E_1$, because $PE_0=PM_0B=PM_0QB=M_1B=E_1$. Now, $(1)$ is equivalent to $(M_0^\top PM_0)Q=M_0^\top M_1$. Hence we may simply perform a polar decomposition $\tilde{P}Q$ on $M_0^\top M_1$ and set $P=(M_0^\top)^{-1}\tilde{P}M_0^{-1}$.
Note that we do not require that $M_0$ and $M_1$ are positive definite in the above. All we need is that $M_0$ is invertible.