Finding geodesics: Lagrangian vs Hamiltonian I have a question referring to how to compute geodesics of a given spacetime (say, Kerr).
I know that the direct way is via the geodesic equation
$$\frac{d^{2}x^{\mu}}{d\lambda^{2}}+\Gamma^{\mu}_{\nu\kappa}\dot{x}^{\mu}\dot{x}^{\kappa}~=~0.$$
But also read that one can write down the geodesic equations using the Lagrangian formalism. From what I have seen so far, there are two approaches: either to write down the Euler-Lagrange equations with Lagrangian 
$$L~=~\frac{1}{2}g_{\mu\nu}\dot{x}^{\mu}\dot{x}^{\nu}.$$
and then Euler-Lagrange equations are exactly the geodesic equations.
OR: write down the Hamiltonian equations using the Lagrangian above, and then use those equations as the geodesic equations.
My question: are Euler-Lagrange and Hamiltonian approaches fully equivalent when it comes to writing down the geodesic equations? Does one have advantage over the other?
 A: They are all equivalent. The answer to your other question is: the Hamiltonian approach usually works best. 
Geodesics can be defined a few ways, since the connection of spacetime is taken to be Levi-Civita. 
Let $(\mathcal{M},g)$ denote spacetime $\mathcal{M}$ with metric $g$ and $\gamma:\mathbb{R}\to M,t\mapsto \gamma(t)$ be a curve on $\mathcal{M}$. If $\nabla$ is the Levi-Civita connection on spacetime, then the geodesic equation is
$$\nabla_{\dot\gamma}\dot\gamma=0,$$
where $\dot\gamma$ is the tangent vector of $\gamma$. Some manipulations can bring this into the form given in the OP. On the other hand, we can define the length functional
$$\ell [\gamma]:=\int_a^b\sqrt{-g(\dot\gamma,\dot\gamma)}\,\mathrm{d}t.$$
(The negative appears because $\gamma$ is usually taken to be timelike for GR purposes.) Then it may be shown that
$$\frac{\delta\ell}{\delta\gamma}=0\Longleftrightarrow\nabla_{\dot\gamma}\dot\gamma=0.$$
So  the geodesic problem of GR is actually an Euler-Lagrange type problem with Lagrangian $L=\sqrt{-g(\dot\gamma,\dot\gamma)}$. 
We can also get a Hamiltonian approach. We first note that the energy functional
$$E[\gamma]=\frac{1}{2}\int_a^bg(\dot\gamma,\dot\gamma)\,\mathrm{d}t$$
has identical Euler-Lagrange equations as the length functional. Proof: Let $D$ be the operator $$Df=\frac{\mathrm{d}}{\mathrm{d}t}\frac{\partial f}{\partial \dot x}-\frac{\partial f}{\partial x}$$
where $f=f(x,\dot x,t)$, so that $Df=0$ is the Euler-Lagrange equation for $f$. Then a short computation shows that $D\ell=\ell\cdot DE$, so $D\ell=0\implies DE=0$. We note that $H:=\frac{1}{2}g(\dot\gamma,\dot\gamma)$ is reminiscent of the $\frac{1}{2}m\mathbf{v}^2$ term of the classical free particle Hamiltonian. (Hence the name energy functional.) 
We can do this more fancily on the cotangent bundle $T^*\mathcal{M}$. Locally trivialize the cotangent bundle in a chart, so we have coordinates $(x^\mu,p_\mu)$. Then put
$$H(x,p):=\frac{1}{2}g^{\mu\nu}(x)p_\mu p_\nu.$$
The Hamiltonian equations
$$\dot x^\mu=\frac{\partial H}{\partial p_\mu},\quad \dot p^\mu=-\frac{\partial H}{\partial x^\mu}$$
are equivalent to the geodesic equation. 
The flow obtained from these equations is a Hamiltonian flow. 
It turns out that the first Hamilton approach is very useful. It is quite easy to vary the integral $E[\gamma]$. This, combined with the Killing method for obtaining first integrals of the geodesic equation, make for a better strategy than computing the Christoffel symbols and solving the equations outright. 
It should be noted that there is another method for solving the geodesic equations: apply the methods of Hamilton-Jacobi theory to the Hamiltonian system described above. This method is used in N. Straumann, General Relativity (2013) to find the geodesics of Kerr spacetime. See also V.I. Arnold, Mathematical Methods of Classial Mechanics (1989) for more on Hamiltonian systems in general. 
It should be noted that the Lagrangian and Hamiltonian described above are not actually Legendre transforms of each other. 
A: *

*The Hessian matrix
$$\tag{1} \frac{\partial^2 L_0}{\partial\dot{x}^{\mu}\partial\dot{x}^{\nu}}~=~2g_{\mu\nu}(x), \qquad \det g_{\mu\nu}(x)~\neq~0,$$
of OP's Lagrangian 
$$\tag{2} L_0~=~g_{\mu\nu}(x)~\dot{x}^{\mu}\dot{x}^{\nu}$$
is non-singular. This means that the Legendre transformation is regular, and therefore that the Lagrangian and Hamiltonian formulations are equivalent, cf. e.g. this Phys.SE post.

*However, since OP mentions geodesics, which extremize lengths of curves, it should be mentioned that the underlying Lagrangian is actually the square root Lagrangian
$$\tag{3} L~=~-\sqrt{-g_{\mu\nu}(x)~\dot{x}^{\mu}\dot{x}^{\nu}}$$
rather than OP's Lagrangian (2).

*The equivalence between the two Lagrangians (2) and (3) is discussed in e.g. this Phys.SE post in the Riemannian case.

*Interestingly, the Hessian matrix of the square root Lagrangian (3) is in fact singular, i.e. the Legendre transformation is singular, and there appears a constraint. (This is related to the worldline reparametrization invariance of the corresponding action for (3).) However, one may still show the equivalence between the Lagrangian and Hamiltonian formulations via e.g. Dirac-Bergmann analysis. This is e.g. done in this Phys.SE post.
