  We develop regularity theory for critical points of variational integrals
defined on Hessian spaces of functions on open, bounded subdomains of
$\mathbb{R}^n$, under compactly supported variations. The critical point solves
a fourth order nonlinear equation in double divergence form. We show that for
smooth convex functionals, a $W^{2,\infty}$ critical point with bounded Hessian
is smooth provided that its Hessian has a small bounded mean oscillation (BMO).
We deduce that the interior singular set of a critical point has Hausdorff
dimension at most $n-p_0$, for some $p_0 \in (2,3)$. We state some applications
of our results to variational problems in Lagrangian geometry. Finally, we use
the Hamiltonian stationary equation to demonstrate the importance of our
assumption on the a priori regularity of the critical point.
