Given an $n \times n$ orthostochastic matrix $\mathbf{A}$, i.e., there exists an orthogonal matrix $\mathbf{O}$ with $A_{ij} = O_{ij}^2$ for all $1\leq i,j \leq n$. What is the fastest way to find $\mathbf{O}$ from $\mathbf{A}$?

As the missing information is only the signs, a complete search over all sign patterns would solve the problem. However, this involves testing $2^{N^2}$ sign patterns. The number of sign patterns can be reduced slightly by some equivalence operations like sign inverting a row. Also, there are $\textit{sign patterns allowing orthogonality}$, i.e., collection of all sign patterns which occur in orthogonal matrices. Besides that the characterization of sign patterns allowing orthogonality is still open for $n > 6$, it will probably not reduce the number of possible candidates to a polynomial magnitude.

There is an analogous case of reconstructing unitary matrices from unistochastic matrices. There is some literature on the analytic side like "The Importance of Being Unistochastic" by Bengtsson, I. (2004). However, I am more concerned with an efficient algorithmic solution to this problem. Any ideas?

A local improvement can be achieved by the following: given $\mathbf{A}$ and a random $\pm$ sign matrix $\mathbf{S}$. Then we can find the closest orthogonal matrix $\mathbf{P}$ by solving the Procrustes problem $\min \lVert \mathbf{P} - \mathbf{A} \circ \mathbf{S} \rVert_F$ with the singular value decomposition, where $\circ$ is the Hadamard product. If now, $\textrm{sign}(\mathbf{P}) \neq \mathbf{S}$ then we can improve by solving with replacing $\mathbf{S}$ with $\textrm{sign}(\mathbf{P})$, because $\mathbf{A}$ is non-negative. This works until the sign matrices are equal.