Building on the blueprint from Goemans and Williamson (1995) for the Max-Cut problem, we construct a polynomial-time approximation algorithm for orthogonally constrained quadratic optimization problems. First, we derive a semidefinite relaxation and propose a randomized rounding algorithm to generate feasible solutions from the relaxation. Second, we derive constant-factor approximation guarantees for our algorithm. When optimizing for m orthonormal vectors in dimension n, we leverage strong duality and semidefinite complementary slackness to show that our algorithm achieves a 1/3-approximation ratio. For any m of the form 2^q for some integer q, we also construct an instance where the performance of our algorithm is exactly (m+2)/(3m), which can be made arbitrarily close to 1/3 by taking m large enough, hence showing that our analysis is tight.