  We establish small energy H\"{o}lder bounds for minimizers $u_\varepsilon$ of
\[E_\varepsilon (u):=\int_\Omega W(\nabla u)+ \frac{1}{\varepsilon^2}
\int_\Omega f(u),\] where $W$ is a positive definite quadratic form and the
potential $f$ constrains $u$ to be close to a given manifold $\mathcal N$. This
implies that, up to subsequence, $u_\varepsilon$ converges locally uniformly to
an $\mathcal N$-valued $W$-harmonic map, away from its singular set. We treat
general energies, covering in particular the 3D Landau-de Gennes model for
liquid crystals, with three distinct elastic constants. Similar results are
known in the isotropic case $W(\nabla u)=\vert \nabla u\vert^2$ and rely on
three ingredients: a monotonicity formula for the scale-invariant energy on
small balls, a uniform pointwise bound, and a Bochner equation for the energy
density. In the level of generality we consider, all of these ingredients are
absent. In particular, the lack of monotonicity formula is an important reason
why optimal estimates on the singular set of $W$-harmonic maps constitute an
open problem. Our novel argument relies on showing appropriate decay for the
energy on small balls, separately at scales smaller and larger than
$\varepsilon$: the former is obtained from the regularity of solutions to
elliptic systems while the latter is inherited from the regularity of
$W$-harmonic maps. This also allows us to handle physically relevant boundary
conditions for which, even in the isotropic case, uniform convergence up to the
boundary was open.
