Estimates Hessian of a scalar-valued function, and returns it
in a sparse Matrix format. The sparsity pattern must be known in advance. The
algorithm is especially efficient for hierarchical models with a large number of
heterogeneous units. See Braun, M. (2017)