optimize.approx_fprime & jacobians

[MartinSavko]

Hello,
While trying to run scipy.optimize.check_grad function in the following code

!/usr/bin/python

'''Fitting the function y = c + r*sin(phi + alpha)'''
from math import sin, cos, pi
import scipy.optimize as op
import numpy as np

def f(als, c=1., r=1., phi=15*pi/180.):
return np.array([y(a, c, r, phi) for a in als])

def y(a, c=1., r=1., phi=15*pi/180.):
return c + r * sin(phi + a)

def y_grad(als, c=1., r=1., phi=15*pi/180):
return np.array([y_prime(a, c, r, phi) for a in als])

def y_prime(a, c=1., r=1., phi=15*pi/180):
return -r * cos(phi + a)

def diffsquares(y_exp, y_fit):
return (y_exp - y_fit)**2

step = 2_pi / 360.
a = np.arange(0., 2_pi + step, step)

result = op.check_grad(f, y_grad, a)

print 'result', result

on Ubuntu 12.04 (scipy version '0.9.0') I get the following error:

Traceback (most recent call last):
File "./fit.py", line 28, in
result = op.check_grad(f, y_grad, a)
File "/usr/lib/python2.7/dist-packages/scipy/optimize/optimize.py", line 382, in check_grad
return sqrt(sum((grad(x0,_args)-approx_fprime(x0,func,_epsilon,_args))*2))
File "/usr/lib/python2.7/dist-packages/scipy/optimize/optimize.py", line 377, in approx_fprime
grad[k] = (f(((xk+ei,)+args)) - f0)/epsilon
ValueError: setting an array element with a sequence.

looking in the source code I think there is an obvious typo there, the numerator is not indexed, I have checked the code on gitHub and the bug seems to be still present (line 610 https://github.com/scipy/scipy/blob/master/scipy/optimize/optimize.py) .

line 610
grad[k] = (f(*((xk + d,) + args)) - f0) / d[k]

should probably read
grad[k] = (f(*((xk + d,) + args))[k] - f0[k]) / d[k]

Cheers,
Martin

[josef-pkt]

f is supposed to return a scalar, see docstring

From what I remember when I looked at the function last time, everything worked as expected.

[pv]

Yes, it works correctly. It evaluates only gradients, however, not Jacobians.

[ev-br]

Can be seen as a documentation issue. From reading the docstring of check_gradient for the first time, I didn't figure out whether f is supposed to be a a scalar function of a vector argument (and what is being checked is indeed a gradient), or is the requirement that x0 is the array just a way of vectorizing a computation of a scalar derivative of a scalar function at a series of points.

The docstring of approx_fprime is way clearer. A quick fix could be to add an example to the docstring. Something as simple as

import scipy.optimize as op
import numpy as np

def fv(x):
    return x*x

def fv_prime(x):
    return 2.*x

def fs(x):
    return np.sum(x*x)

print op.check_grad(fs, fv_prime, np.array([3]))

On a related note, I wonder is there's a pressing reason for check_grad refusing to deal with scalar arguments, not arrays:

op.check_gradient(fs, fv_prime, 3.)

[pv]

Ok, yep, the documentation is a bit confusing here.

@MartinSavko: in principle it would be OK to extend approx_fprime to evaluate also Jacobians of vector-valued functions, and not only gradients of scalar functions. If you decide to do that, please file a second pull request, thanks!

[andyfaff]

At the present time It would be better to expose optimize._differentiable_functions.VectorFunction, or optimize._numdiff.approx_derivative as a public interface.
The latter function can do either gradients (scalar functions) or Jacobians (vector valued functions).
I'm in the process of making approx_fprime rely on approx_derivative for calculation.

[mdhaber]

Based on the comments, looks like the doc fix is a good first issue.

[Jacob-Stevens-Haas]

So it looks like approx_fprime does rely on approx_derivative now. Moreover, it only performs two main steps: enforcing that x is an array and enforcing that f(x) is scalar. Then it calls approx_derivative, which can handle vector-valued functions.

So the question is, I could (a) make check_grad call approx_derivative directly, or (b) remove the step in approx_fprime that restricts f(x) to be a scalar. I would imagine that there's a reason it was made to be a scalar to support some other higher-level functions. On the other hand, not sure if you want check_grad to directly call the non-public optimize._numdiff.approx_derivative.

[andyfaff]

It might be better to keep approx_derivative private at this point.

approx_fprime could probably do with an overhaul. For example, the function itself is no longer used in minimisation, so the last sentence in the notes section is out of date.

It should be fairly straightforward to allow approx_fprime to calculate Jacobians (R^m -> R^n) instead of solely being used to calculate gradients (R^m -> 1), would you like to add that functionality, it'll be a worthwhile PR.

check_grad could be changed to call approx_derivative directly as well, but it's not clear to me that check_grad is actually all that useful.

[ev-br]

One thing to keep in mind if someone's willing to work on finite differencing: it might be better to separate scalar derivatives, gradients of vector functions and jacobians early on because of broadcasting.