Consider the "positive part" of the hyperboloid model $$ \mathbb H_+^d := \{ (x_1, \ldots, x_d, x_{d + 1}) \in \mathbb R^{d} \times \mathbb R_{> 0}: \sum_{k = 1}^{d} x_k^2 = x_{d + 1}^2 - 1 \}. $$ How does the orthogonal projection $\text{proj}_{\mathbb H_+^d}(p)$ of a point $p \in \mathbb R^{d + 1}$ onto $\mathbb H_+^d$ look like?
My work. We have $\text{proj}_{\mathbb H_+^d}(p) = \text{argmin}_{x \in \mathbb H_+^d} \| p - x \|_2^2$, where $\| \cdot \|_2$ denotes the Euclidean norm. Setting up the Lagrange problem $$ \min_{x \in \mathbb R^{d + 1}} \| p - x \|_2^2 + \lambda \left( \sum_{k = 1}^{d} x_k^2 - x_{d + 1}^2 + 1\right) + \mu x_{d + 1} $$ for some Lagrange multipliers $\mu,\lambda \in \mathbb R$ yields the optimality condition (setting the gradient with respect to $x$ to zero) $$ x_k = \frac{p_k}{1 - \lambda} \; \forall k \in \{ 1, \ldots, d \} \qquad \text{and} \qquad x_{d + 1} = \frac{p_{d + 1}}{1 + \lambda} + \mu. $$ Plugging this into the objective and trying to solve for $\lambda$ yields an equation which is quartic in $\lambda$ and has many and very complicated solutions.
Is there another way to find this projection?