[SciPy-user] Using SciPy/NumPy optimization

Alok Singhal as8ca at virginia.edu
Tue Feb 27 16:46:29 EST 2007


Hi,

Sorry about the long post.  I did not ask the original question, but I
tried the solution posted and could not get it to work.

Original question: given a set of data [[x1, y1], ..., [xn, yn]], how
to fit y = a*(1 + b*c*x) ** (-1.0/b) to it, subject to the
constraints:

 0 <  b <= 5
-1 <= c <= 1
a ~ >= max(x1, ... xn)

On 27/02/07: 14:25, Robert Kern wrote:
> There are a number of constrained optimizers in scipy.optimize .
> scipy.optimize.fmin_tnc seems most appropriate for simple bounds like you have.
> In order to get a function to minimize that depends on data, I usually like to
> use a class:
> 
>   import numpy as np
>   from scipy.optimize import fmin_tnc
> 
>   class LossFunction(object):
>     def __init__(self, x, y):
>       self.x = x
>       self.y = y
> 
>     def __call__(self, abc):
>       """ A function suitable for passing to the fmin() minimizers.
>       """
>       a, b, c = abc
>       y = a*(1.0 + b*c*self.x) ** (-1.0/b)
>       dy = self.y - y
>       return dy*dy
> 
>   x = np.array([...])
>   y = np.array([...])
>   lf = LossFunction(x, y)
>   abc0 = np.array([x.max(), 2.5, 0.0])  # or whatever
>   retcode, nfeval, abc_optimal = fmin_tnc(lf, abc0,
>     bounds=[(None, None), (0., 5.), (-1., 1.)])

I used optimize.leastsq to do the fitting, and it worked well me (see
http://www.scipy.org/Cookbook/FittingData).  But I was trying the
above method, and unfortunately it doesn't work as presented.  When I
define the variables x and y by:

from scipy.optimize import fmin_tnc
from scipy import rand
from numpy import mgrid, log

def f(abc, x):
  return abc[0]*(1.0+abc[1]*abc[2]*x)**(-1.0/abc[1])

a = 15.0
b = 2.5
c = 0.3
abc = (a, b, c)
num_points = 151
x = mgrid[1:10:num_points*1j]

y = f(abc, x) + rand(num_points) - 0.5
abc0 = [x.max(), 2.5, 0.0]
lf = LossFunction(x, y)

retcode, nfeval, abc_optimal = fmin_tnc(lf, abc0, bounds=[(None, None), (0., 5.), (-1., 1.)])

The call to fmin_tnc gives me an error:

File "/usr/lib/python2.4/site-packages/scipy/optimize/tnc.py", line 191, in func_and_grad
    f, g = func(x, *args)
ValueError: too many values to unpack

Looking at the documentation and the source file, the first parameter
to fmin_tnc should be a function that returns the function value and
the gradient vector (with respect to the parameters being estimated).
Also, it seems that I can tell fmin_tnc to estimate the gradient
itself, if I specify approx_grad as True:

retcode, nfeval, abc_optimal = fmin_tnc(lf, abc0, bounds=[(None, None), (0., 5.), (-1., 1.)], approx_grad=True)

But this gives me an error saying:

File "/usr/lib/python2.4/site-packages/scipy/optimize/optimize.py", line 576, in approx_fprime
    grad[k] = (apply(f,(xk+ei,)+args) - f0)/epsilon
ValueError: setting an array element with a sequence.

Then I tried redefining the __call__ in LossFunction to return the
gradient as well:

class LossFunction(object):
  def __init__(self, x, y):
    self.x = x
    self.y = y
  def __call__(self, abc):
    a, b, c = abc
    y = a*(1.0 + b*c*self.x) ** (-1.0/b)
    dy = self.y - y
    d = [0, 0, 0]
    d[0] = (1.0 + b*c*self.x) ** (-1.0/b)
    d[1] = y * (log(a*(1.0+b*c*self.x))/b**2 - c*self.x/(b*(1.0+b*c*self.x)))
    d[2] = -a/b*(1.0 + b*c*self.x) ** (-1.0/b) * b*self.x
    return dy*dy, d

When I call this function instead, I get another error:

File "/usr/lib/python2.4/site-packages/scipy/optimize/tnc.py", line 221, in fmin_tnc
    fmin, ftol, rescale)
ValueError: Bad return value from minimized function.

At this point, I gave up, but I would like to understand what I am
doing wrong.  To me, it seems like fmin_tnc will find the values of
parameters (a, b, c in this case) that minimize a given function that
depends *only* on a, b, c.  I don't understand how to make it work so
as to fit existing data.  I am sure I am missing something, but I
don't know what that is.

Thanks for all your help,
Alok

-- 
Alok Singhal                               *   *          
Graduate Student, dept. of Astronomy   *           *     *
University of Virginia                                    
http://www.astro.virginia.edu/~as8ca/              *    * 



More information about the SciPy-User mailing list