ENH: optimize._hessian: use _jacobian to compute accurate Hessian

[mdhaber]

Reference issue

Follow-up to gh-20630

What does this implement/fix?

Adds a private function to compute accurate Hessians using a finite difference approximation. Vectorized to compute Hessian at an arbitrary number of points in a single call.

Additional information

To do eventually:

• Test (and fix, where applicable) attributes of return object (for this and _jacobian). It doesn't look like nfev is correct. Fixed: only sum the inner nfev counts for the first nit outer iterations.
• Improve error estimate by considering error of _jacobian. Just tighten tolerance of inner Jacobian. In the future we can use a legit 2-D stencil.
• Consider adding features to help user vectorize with respect to multiple parameters. (Separate PR; not essential.)
• Consider using first Jacobian evaluation to estimate appropriate step size. This could reduce number of iterations required in subsequent evaluations. (Separate PR; not essential.)

[mdhaber]

I think the failure occurs when the errors are equal; there would need to be an additional condition to ensure that corresponding off-diagonals of i are mutually exclusive.

An unrelated issue that is probably rare but could happen: what if the lower error corresponds with a nonzero status. Probably need to make the condition contingent on success.

I'm conflicted about what to do here, so I'll wait for feedback.