I'd like to find the set of positive definite matrices $W$ which satisfy the equation
$$V = \Psi \left(W - \Psi^{-1/2}\Phi^{1/2}\right)\left(W - \Phi^{1/2}\Psi^{-1/2}\right)W^{-1}$$
where $V, \Psi,$ and $\Phi$ are known positive definite matrices. I'd appreciate any help or references on solving matrix equations of this sort.
Some context: I'm trying to simulate from a particular matrix variate probability distribution. It's not obvious how to simulate from the distribution for $W,$ but (for particular parameter values) the distribution of $V$ (a transformation of $W$) is easy to simulate from. If there were just one positive definite solution to this equation, one could simulate $V$ and then solve for $W$ and be done. If there is more than one solution, one can still sometimes use a similar approach, but assign a probability to each solution. The article https://www.jstor.org/stable/pdf/2683801.pdf has more details. My hope is that there are a finite number of positive definite solutions so that I can apply this sort of trick.
I assume that "positive definite" in your post means "self-adjoint strictly positive definite". Then the solution is, indeed, unique (in the class of positive definite matrices). To show that, just rewrite the equation as $(W-A)(W-A^*)=UW$ where $A=\Psi^{-1/2}\Phi^{1/2}, U=\Psi^{-1}V$. Since the LHS is self-adjoint, we must have $UW=WU^*$, so, replacing $UW$ by $\frac 12[UW+WU^*]$ and opening the parentheses, we see that $W$ satisfies $$ W^2=B^*W+WB+Q $$ with $B=\frac 12U^*+A^*$ and positive definite $Q=AA^*$, both known. That is the classical Riccati equation (just google the term to see a vast literature devoted to it) but since some rudiments of the theory relevant to your question seem to be omitted from the Wikipedia article and from most research papers you can see online as well-known, I'll include them here.
One possible approach is to notice that the span of the columns of $\begin{bmatrix}I\\W\end{bmatrix}$ is an eigenspace of the $2n\times 2n$ matrix $M=\begin{bmatrix}B&-I\\-Q&-B^*\end{bmatrix}$ (this is just the claim that the linear space of all vectors of the kind $\begin{bmatrix}x\\Wx\end{bmatrix}$ is $M$-invariant, which immediately follows from honest multiplication and observing that the equation can be rewritten as $-Q-B^*W=W(B-W)$. Now, if $\lambda\in C$ is an eigenvalue of $M$ with the corresponding eigenvector $\begin{bmatrix}x\\y\end{bmatrix}$, then $-\lambda$ is an eigenvalue of $M^*$ with the corresponding eigenvector $\begin{bmatrix}-y\\x\end{bmatrix}$ (direct verification), so $-\bar\lambda$ is also an eigenvalue of $M$.
The eigenvector equation for $M$ is equivalent to the pair of equations $Bx-y=\lambda x, -Qx-B^*y=\lambda y$. Note now that they can be rewritten as $y=(B-\lambda I)x, -Qx-2\Re\lambda y=(B^*-\bar\lambda I)y$, so $-Qx-2\Re\lambda y=(B^*-\bar\lambda I)(B-\lambda I)x$. Assume now that we are in our eigenspace, i.e., $y=Wx$. Taking the scalar product with $x$, we see that the right hand side is non-negative (the operator is positive semidefinite) while the left hand side is $-\langle Qx,x\rangle-2\Re\lambda\langle Wx,x\rangle$. Since both scalar products are strictly positive, we conclude that it may happen only for $\Re\lambda<0$, which tells us that if a positive definite solution $W$ exists at all (and, if I understand you right, the existence follows from some other considerations in your setting), then we can use only eigenvalues with the negative real part (which is what people usually do anyway though for a different reason), in which case, due to the mentioned above property that all eigenvalues except purely imaginary ones come in pairs $\lambda,-\bar\lambda$ with corresponding eigenspaces having the same dimension, we see that we must have exactly $n$-dimensional eigenspace corresponding to all eigenvalues with negative real part for your problem to be solvable and that eigenspace produces a unique solution.
To be completely honest, one needs also to discuss the possibility of non-trivial Jordan cells here. I'll just restrict myself to the observation that $M$ and $M^*$ have the same Jordan structure and the relation between the generalized eigenvectors of $M$ and $M^*$ is the same as between the eigenvectors (to check the latter, just observe that if $(M-\lambda I)\begin{bmatrix}x\\y\end{bmatrix}=\begin{bmatrix}u\\v\end{bmatrix}$, then $(M^*+\lambda I)\begin{bmatrix}-y\\x\end{bmatrix}=\begin{bmatrix}v\\-u\end{bmatrix}$), so the cells corresponding to the eigenvalues $\lambda$ and $-\bar\lambda$ in $M$ are of the same dimensions and that if some generalized eigenvector is in an $M$-invariant subspace, then some eigenvector with the same eigenvalue is in that subspace too, so the conclusion is still the same: we should use all eigenvalues with negative real part and to include the full Jordan blocks corresponding to them into our $n$-dimensional eigenspace.