ENH: optimize._jacobian: use _differentiate to compute accurate Jacobian

[mdhaber]

Reference issue

gh-20063
Towards gh-17059

What does this implement/fix?

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

Additional information

Example: evaluate the gradient of the Rosenbrock function in three dimensions at 10 points.

import numpy as np
from scipy.optimize import rosen, rosen_der
from scipy.optimize._differentiate import _jacobian
rng = np.random.default_rng(259824658972546)
x = rng.random(size=(3, 10))
res = _jacobian(rosen, x)
ref = rosen_der(x)
res.df.shape  # (3, 10)
np.testing.assert_allclose(res.df, ref)

(Tests demonstrate evaluation of true Jacobians; i.e., for functions $\mathbf{R}^m \rightarrow \mathbf{R}^n$.)

Yes, this requires that the callable be vectorized. If it is not, it is trivial for the user to vectorize it. For example,

def non_vectorized_func(x1, x2, x3):
    assert np.isscalar(x1)
    assert np.isscalar(x2)
    assert np.isscalar(x3)
    return rosen([x1, x2, x3])

func = lambda x: np.vectorize(non_vectorized_func)(*x)
res2 = _jacobian(rosen, x)
np.testing.assert_allclose(res2.df, ref)

We can include an example in the documentation if need be.

Adding a batch parameter to reduce memory usage (to control the number of partial derivatives that are taken in a single call) can be a separate PR.

Implementation of a _hessian function in terms of _jacobian is simple. The thing that might take a bit of thought there is the status messages, error estimates, function eval counts, etc.. since there are nested calls.

import numpy as np
from scipy.optimize import rosen, rosen_hess
from scipy.optimize._differentiate import _jacobian
rng = np.random.default_rng(259824658972546)

def hessian(f, x):
    def df(x):
        return _jacobian(f, x).df
    return _jacobian(df, x)

x = rng.random(3)
ref = rosen_hess(x)
np.testing.assert_allclose(hessian(rosen, x).df, ref, atol=1e-9)  # atol needed only for zeros

[steppi]

I think this is great. I'd be happy to merge today if you think that's fine @mdhaber. Agreed that documenting how to vectorize a callable and adding a batch parameter can each be done in follow-ups.

Edit: Oh, I guess we need to sort out the mypy stuff. I'm not well versed in type hinting, is this something that could be solved by adding a .pyi file?

[mdhaber]

OK this commit adds "finite real" and should address mypy's complaints.