docs: def eval_hessian(self, ordered_parameters=[], **parameters):
          """
          Hessian of :math:`S` in the
          :class:`~symfit.core.argument.Parameter`'s (:math:`\\nabla_\\vec{p}^2 S`).
  
          :param parameters: values of the
              :class:`~symfit.core.argument.Parameter`'s to evaluate :math:`\\nabla_\\vec{p} S` at.
          :return: ``np.array`` of length equal to the number of parameters..
          """
          evaluated_func = super(LeastSquares, self).__call__(
              ordered_parameters, **parameters
          )
          evaluated_jac = super(LeastSquares, self).eval_jacobian(
              ordered_parameters, **parameters
          )
          evaluated_hess = super(LeastSquares, self).eval_hessian(
              ordered_parameters, **parameters
          )
  
          result = 0
          for var, f, jac_comp, hess_comp in zip(self.model.dependent_vars,
                                                 evaluated_func, evaluated_jac,
                                                 evaluated_hess):
              y = self.dependent_data[var]
              sigma_var = self.model.sigmas[var]
              if y is not None:
                  sigma = self.sigma_data[sigma_var]
                  p1 = hess_comp * ((y - f) / sigma**2)[np.newaxis, np.newaxis, ...]
                  # Outer product
                  p2 = np.einsum('i...,j...->ij...', jac_comp, jac_comp)
                  p2 = p2 / sigma[np.newaxis, np.newaxis, ...]**2
                  # We sum away everything except the matrices in the axes 0 & 1.
                  axes = tuple(range(2, len(p2.shape)))
                  result += np.sum(p2 - p1, axis=axes, keepdims=False)
          return np.atleast_2d(np.squeeze(np.array(result)))