This method does not fully parameterize the space of orthogonal matrices. It instead defines a chart over a Steifel manifold that will necessarily always have a slice of parameters that it cannot describe (i.e. the chart is not isomorphic to the entire space). All of the issues with reflections and reorientations are consequences of this.
It’s the same story that occasionally pops up about circular parameters – if you try to parameterize them as [0, 2pi] or [-pi, pi] then you run into consistent issues at the boundaries. If the posterior is sufficiently concentrated then it is possible to orient the coordinates so that the posterior doesn’t have appreciable mass across the boundary, but it’s tricky to verify when the posterior is tight enough and even then you have to know where it’s concentrating. This typically results in an awkward, iterative process.
That same process is technically possible for orthogonal (equivalently orthonormal) matrices as implied by the paper but it becomes extraordinarily harder as you now have to understand how the posterior concentrates in a multidimensional space.
Models over the space of orthogonal matrices are becoming more and more relevant, but the only way that they can principly be implemented in Stan (and all of these identification issues avoided) is to seriously change the algorithmic backend to enable the geodesics as described by Simon and Mark’s paper.
Everyone is welcome to experiment with this code, but please be very careful as failures many not be immediately obvious. Again, see how subtle one-dimensional circular parameters are to fit using similar techniques already.
</methodological pedant>