Finding multiple solutions to $a = Ra$ with Symmetric Orthogonal $R \in SO(3), R = R^T$ and $a \in Re(3)$ I have a problem where I wish to find multiple solutions to the equation:
$ \mathbf{a} = \mathbf{R} \mathbf{a}$
where I wish to calculate the solutions for $\mathbf{R}$ that is both symmetric ($\mathbf{R} = \mathbf{R}^T$) and special orthogonal ($\mathbf{R} \in SO(3)$) from a given $\mathbf{a}$ that is a real-valued 3-dimensional vector of non-zero magnitude.
For example, if I have the following vector:
$\mathbf{a} = \begin{bmatrix} 1.70349 && -1.67761 && 9.51419\end{bmatrix}^T$
Then two solutions are $\mathbf{R} = \mathbf{I}$ and
$\mathbf{R} = \begin{bmatrix} 
-0.939693 && -0.0593912  && 0.336824 \\
-0.0593912 && -0.941511 && -0.331707 \\
  0.336824 && -0.331707 &&  0.881204 
\end{bmatrix}$
This second solution can be generated from $\mathbf{R} = \mathbf{Q} \mathbf{R}_z \mathbf{Q}^T$ where $\mathbf{R}_z$ is a rotation about the z-axis by $\pi$ and
$\mathbf{Q} = \begin{bmatrix} 
0.984808 && 0.0301537  && -0.17101 \\
        0  && 0.984808 &&  0.173648 \\
 0.173648  && -0.17101 && 0.969846
\end{bmatrix}$
How would one compute the non-trivial solution? Does this approach generalize to an extended problem where the vectors are different, i.e. solving for the solution(s) for $\mathbf{R}$ in the equation $\mathbf{a} = \mathbf{R} \mathbf{b}$, given $\mathbf{a}$ and $\mathbf{b}$  ?
 A: (The following will be contingent on $a$ being a unit vector, but you can easily get to that case by normalizing.)
In an orthonormal basis $B$ containing $a$, this is easy. Let's say $a$ is the third basis vector, and that its orientation is the same as of the standard basis (so the basis change matrix is in $\operatorname{SO}(3)$, though that's only cosmetic), then the matrix representation $\tilde R$ wrt this basis must be of the form
$$\begin{pmatrix}\cos\varphi&-\sin\varphi&0\\
\sin\varphi&\cos\varphi&0\\
0&0&1\end{pmatrix}$$
for some $\varphi\in[0,2\pi)$. The condition ${\tilde R}^T=\tilde R$ enforces $\varphi\in\{0,\pi\}$, one of which results in the identity matrix, while the other result in the solution you described after a change of basis: Let $Q$ be the change of basis matrix from the standard basis to our orthonormal basis $B$. Then
$$R=Q^T\tilde RQ.$$
The case $a=Rb$ is harder. I was thinking of considering the matrix representation wrt two different orthonormal bases, one containing $a$, the other $b$, and then doing the same as above. But this way, $R^T=R$ might no longer be true after a change of bases. So I'm out of ideas for now.
A: A nice approach is to use the fact that every $n\times n$ real orthogonal matrix is generated by the product of at most $n$ reflection (Householder) matrices.
(all vectors have unit length )
i.e. $Q^{(k)}=I-2\mathbf v_k\mathbf v_k^T$
$\det\big(Q^{(k)}\big) = -1$
In your problem $n=3$ and $\det\big(R\big) = 1$ so $R$ must be the product of exactly 2 Householder matrices. (Note the involution $\big(Q^{(k)}\big)^2=I$ so even the identity matrix may be written as a product of 2 Householder matrices.)
Select the mutually orthonormal set  $\{\mathbf a, \mathbf v_1,\mathbf v_2\}$ so
$\Big(R\Big)\mathbf a=\Big(\big(I-2\mathbf v_1\mathbf v_1^T\big)\big(I-2\mathbf v_2\mathbf v_2^T\big)\Big)\mathbf a=\big(I-2\mathbf v_1\mathbf v_1^T\big)\mathbf a=\mathbf a$
You have an incremental constraint of $R=R^T$ (or equivalently the involution $R^2=I$).  This is satisfied in the above since we ensured $\mathbf v_1^T\mathbf v_2=0$
$R $
$=\big(I-2\mathbf v_1\mathbf v_1^T\big)\big(I-2\mathbf v_2\mathbf v_2^T\big) $
$= I -2\mathbf v_1\mathbf v_1^T - 2\mathbf v_2\mathbf v_2^T$
$= \big(I-2\mathbf v_2\mathbf v_2^T\big)\big(I-2\mathbf v_1\mathbf v_1^T\big)$
$=R^T$
as for the second question
$ R \mathbf{b}= \mathbf{a} $
First select the Householder matrix such that
$\big(I-2\mathbf v_2\mathbf v_2^T\big)\mathbf b=\mathbf a$
(check: $\mathbf v_2 \propto \mathbf a -\mathbf b$)
then select $\mathbf v_1$ such that $\mathbf v_1 \perp \mathbf a$ and $\mathbf v_1 \perp \mathbf v_2$,
Thus
$\big(I-2\mathbf v_1\mathbf v_1^T\big)\big(I-2\mathbf v_2\mathbf v_2^T\big)\mathbf b=\big(I-2\mathbf v_1\mathbf v_1^T\big)\Big(\big(I-2\mathbf v_2\mathbf v_2^T\big)\mathbf b\Big)=\big(I-2\mathbf v_1\mathbf v_1^T\big)\mathbf a = \mathbf a$
so
$R:=\big(I-2\mathbf v_1\mathbf v_1^T\big)\big(I-2\mathbf v_2\mathbf v_2^T\big)$ and
$R=R^T$
as desired
A: The answer can be found using a variation of Wahba's Problem.
Since R is symmetric, it implies that:
$$
\begin{align}
a &= Rb \\
b &= Ra
\end{align}
$$
Following Ref [1], we may then construct a symmetric matrix:
$$ B = \frac{1}{2} \left(ab^T + ba^T \right) $$
And compute the SVD:
$$ B = U S V^T $$
Since $B$ is symmetric, $U = V$.
If $B$ is of rank 2 (i.e. $a \neq b$), then we may compute the solution as:
$$
R = U_+ V_+^T
$$
where $U_+$ and $V_+$ are forced to be special orthogonal by:
$$
\begin{align}
U_+ &= U \operatorname{diag} \left( \begin{bmatrix} 1 && 1 && \det(U)\end{bmatrix} \right) \\
V_+ &= V \operatorname{diag} \left( \begin{bmatrix} 1 && 1 && \det(V)\end{bmatrix} \right) \\
\end{align}
$$
If $B$ is of rank 1 (i.e. $a=b$) then we may introduce a single axis rotation $W$ about the x-axis:
$$
R = U_+ W V_+^T
$$
To keep R symmmetric, we have two choices of W:
$$
\begin{align}
W &= I \\
W &= \operatorname{diag} \left( \begin{bmatrix} 1 && -1 && -1 \end{bmatrix} \right)
\end{align}
$$
Which corresponds as the two variants that we seek.
References:
[1] Markley, F. L. Attitude Determination using Vector Observations and the Singular Value Decomposition, Journal of the Astronautical Sciences, 1988, 38:245-258
