ENH: Vectorized computations in approx_derivative

[nmayorov]

Jacobian is estimated numerically for a function fun(x). If fun is implemented such that it accepts and returns 2-d arrays (it is very easy to achieve using vstack when fun is coded in terms of numpy functions), then fun will be called only once to compute all necessary values.

The idea is simple, but it can give substantial speed-ups for estimating Jacobians. Currently only least_squares uses approx_derivative. I will fill an issue that least_squares should be adjusted to allow vectorized Jacobian estimation (when this PR is in).

@ev-br @pv, please look when you have time.

[ev-br]

Do we really need to support non-vectorized computations?

[nmayorov]

@ev-br missed your comment. I'm glad that you looked into this.

Someone may have a complicated function (integration of ODE for example), which can't be vectorized efficiently (and vectorizing it is just a tedious and pointless job). So I guess it's necessary to provide both variants, although I agree it would be good to decrease the number of options. Another thing is that non-vectorized version is already used in least_squares.

So what do you think of this pull request?

[ev-br]

I think it's a worthy feature. If there's a function which can handle vectorized evaluations, why taking a perf hit of python looping.

However, I'm a bit wary of adding this much of almost-but-not-completely duplicated code.
Maybe it'd be possible to default to a vectorized evaluation, and recommend users wrap their function into np.vectorize otherwise. Or even bake it in, and use this new flag for deciding whether to np.vectorize or not?

[mike1808]

what is the status of this PR?

[nmayorov]

I think it is a useful addition and generally the right way to do things in scipy (i. e. allowing for "doubly" vectorized functions). I don't know, I think it even can be merged right now, then we can include it into least_squares. A general idea is to make this function public and replace calls of other implementations in optimizeby this function. I guess we will have num. diff. functions in a good shape after this year's GSoC (either in scipy.diff or in scipy.diff and scipy.optimize, the latter being used for numerical differentiation in optimization problems, i. .e Jacobian and Hessians including sparse versions).

By the way, why are you asking?

[mdhaber]

However, I'm a bit wary of adding this much of almost-but-not-completely duplicated code. Maybe it'd be possible to default to a vectorized evaluation, and recommend users wrap their function into np.vectorize otherwise. Or even bake it in...

Yes, I'd suggest the same. If the user function can't do vectorized calculations, then we np.vectorize it. I'd also consider going all the way with the vectorization to support Nd arrays unless it's particularly difficult for some reason.

Since this has been inactive for a while and there are merge conflicts, I'll go ahead and close, but all should feel free to reopen to complete this enhancement.