Principled approach to arrive at geodesic Hamiltonian $H = g^{\mu \nu} p_\mu p_\nu$? The background:
If we have a spacetime path $x^\mu(t)$ parameterized by arbitrary parameter $t$, the proper time along the path between $t_1$ and $t_2$ is
$$
\int_{t_1}^{t_2} (g_{\mu \nu} \dot x^\mu \dot x^\nu)^{1/2}dt.
$$
This action has a reparameterization symmetry. If we take the Lagrangian to be
$$L = (g_{\mu \nu} \dot x^\mu \dot x^\nu)^{1/2}$$
then the Euler Lagrange equation ends up being (after multiplying by the inverse metric)
$$
\ddot x^\mu = - \Gamma^\mu_{\; \beta \gamma} \dot x^\beta \dot x^\gamma + \dot x^\mu \frac{d}{dt} \ln(L).
$$
If we use an affine parameter $t$ such that $L$ is constant along the path, then this is just the regular geodesic equation.
We can see that the reparameterization (gauge) symmetry is a huge pain. For example, it seems to render the Hamiltonian to be $0$.
\begin{align*}
H &= \Big( \frac{\partial L}{\partial\dot x^\mu} \Big)\dot x^\mu - L \\
&= 2\frac{1}{2}\frac{g_{\mu \nu} \dot x^\nu}{(g_{\alpha \beta} \dot x^\alpha \dot x^\beta )^{1/2}} \dot x^\mu - (g_{\mu \nu} \dot x^\mu \dot x^\nu)^{1/2} \\
& = 0.
\end{align*}
This seems to be related to the reparameterization symmetry because $H$ should generate time translations. However, if time evolution is not deterministic, then there's nothing for $H$ to reasonably be.
We would therefore like to find a less "pathological" Lagrangian and Hamiltonian with the same equations of motion, but in an automatically affinely parameterized form. The answer is to take
$$
L = \frac{1}{2} g_{\mu \nu} \dot x^\mu \dot x^\nu \implies H = \frac{1}{2} g^{\mu \nu} p_{\mu} p_{\nu}.
$$
Note that $L_{\rm new} = \tfrac{1}{2}L_{\rm old}^2$. The Euler Lagrange equations are indeed the affinely parameterized geodesic equations. The same goes for Hamilton's equations
$$
\dot x^\mu = \frac{\partial H}{\partial p_{\mu}} , \hspace{1 cm} \dot p_{\mu} = -  \frac{\partial H}{\partial x^{\mu}} \implies \ddot x^\mu = - \Gamma^\mu_{\; \beta \gamma} \dot x^\beta \dot x^\gamma.
$$
The question:
It seems like a complete coincidence that squaring the Lagrangian gives us sort of "gauge fixed" version of our original Lagrangian. Is there a principled, philosophical approach we can take to move from $L_{\rm old}$ to $L_{\rm new}$, rather than just seeing that the equations of motion give us what we want? Is this a special case of a more general procedure?
 A: *

*Physicists conventionally normalize the square root Lagrangian a bit differently, namely as$^1$
$$L_0~:=~ -m\sqrt{-\dot{x}^2}, \qquad \dot{x}^2~:=~g_{\mu\nu}(x)~ \dot{x}^{\mu}\dot{x}^{\nu}~<~0,\tag{1}$$
where $m>0$ is the mass. (OP assumes that $m=1$.)

*It is important to note that the Legendre transformation of the square root Lagrangian (1) is singular: the momentum has to satisfy the mass-shell constraint
$$  p^2+m^2~\approx~0 , \qquad p^2~:=~g^{\mu\nu}(x)~ p_{\mu}p_{\nu}.\tag{2} $$
Therefore, even though OP is correct that the original Hamiltonian is zero because of worldline (WL) reparametrization invariance, the full Hamiltonian 
$$ H~=~ \frac{e}{2}(p^2+m^2)\tag{3} $$
becomes a Lagrange multiplier $e$ times the mass-shell constraint (2).

*The inverse Legendre transformation of the Hamiltonian (3) leads to
the Lagrangian
$$L~=~\frac{\dot{x}^2}{2e}-\frac{e m^2}{2}.\tag{4}$$
If we integrate out the $e$ field, the Lagrangian (4) becomes the square root Lagrangian (1), cf. e.g. this related Phys.SE post.

*The main point is now that OP's quadratic Hamiltonian and Lagrangian are the $e=1$ gauge of the Hamiltonian (3) and the Lagrangian (4), respectively, up to irrelevant constant terms. In that sense they follow systematically from the Dirac-Bergmann algorithm for constrained systems.

*For further details, see e.g. this related Phys.SE post.
--
$^1$ We use the sign convention $(-,+,+,+)$, and put the speed-of-light $c=1$.
