I rencently learned about the lagrangian/variational method for computing Christoffel symbols.
Let $\mathcal{M}$ be a $m-$dimensional manifold with $g_{ij}$ being the metric tensor components and let $\gamma(\lambda) = (x^1(\lambda)\dots x^m(\lambda))$ a geodesic in $\mathcal{M}$. By assigning the following Lagrangian to the metric tensor $g_{ij}$ : $$\mathcal{L} = g_{ij} \dot{x}^i \dot{x}^j,$$ with $\dot{x}^i = \dfrac{dx^i}{d\lambda}$ and using the least action principle one gets the Euler-Lagrange equations : $$\dfrac{d}{d\lambda}\dfrac{\partial \mathcal{L}}{\partial \dot{x}^i} - \dfrac{\partial\mathcal{L}}{\partial x^i} = 0 \tag{1}$$ Remebering the geodesic equations written by means of Christoffel symbols : $$\ddot{x}^k + \Gamma^k_{\ ij} \dot{x}^i \dot{x}^j = 0,\tag{2}$$ Christoffel symbols are obtained comparing $(1)$ and $(2)$.
How does the fact the basis one uses is holonomic or not is taken into account ?
Simply, your eq. (1) is dynamic equation of motion without constrains. The dynamic nonholonomic equation of motion reads:
$\frac{d}{dt}\frac{\partial L}{\partial \dot{x}^i} -\frac{\partial L}{\partial x^i}=\lambda_j a^j_i $,
where $\lambda_j$ - lagrange multipliers defining nonholonomic constrains