I have a two matrices $A$ and $B$ in $\mathbb{R}^{m \times n }$ ($m \gg $ n) such that there exists an orthonormal matrix $X \in \mathbb{R}^{n \times n }$, such that:

$$AX = B$$

Given that $X$ is orthonormal this is also true:

$$A = BX^T$$

**How to find $X$?**

I tried to use Moore-Penrose inversion $A^+$ and got non-orthonormal result $Y = A^+B$, that works only one way:

$$AY = B$$

but not in another. The problem is that both matrices $A$ and $B$ are not absolutely accurate (obtained in numerical calculations). So the non-orthonormal solution $Y$ arises. It is slightly more accurate than the exact solution $X$: $\|AY - B\| \lt \|AX - B\| < 0.001 $. But in the inverse case of course it doesn't work at all: $\|BY^T - A$$\| > 100.0$. Whilst the exact solution is good enough in both ways: $\|BX^T - A\| < 0.001$.

$A$ is full-rank, *i.e.* $A^+A = I$ **(Not really! See P.S.)**

So the question is how to find orthonormal solution of the overdetermined linear equations system?

P.S.

It turned out that my problems were caused by presence small of singular singular values of A and B. So the overall problem was ill-conditioned. I have accepted @Federico-Poloni's answer because it directly addresses the original question and has a reference. James's answer seems to be working either, but I have to choose only one.