[SciPy-User] fsolve with restriction on variable
Ashley DaSilva
I am using fsolve to solve a function, f(v), v=[x,y,z] is a list of three variables. However, I have a factor in f which contains (1-x**2)**(7./2). So, when I do the following,
fsolve(f,x0)
the code eventually tries x= -1.57, which clearly produces an error in my function due to the power 7./2.
I know that the solution for x should be between 0 and 1, is there a way to put this restriction on x while using fsolve?
Chris Colbert
http://docs.scipy.org/doc/scipy/reference/optimize.html
> _______________________________________________
> SciPy-User mailing list
> SciPy...@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
>
_______________________________________________
SciPy-User mailing list
SciPy...@scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user
Ashley DaSilva
I see that there are some constrained optimization techniques, but I don't want to find the minimum/maximum of my function, I want to find the root. I also see that there are some scalar root finding techniques that allow constraints (brentq,brenth,ridder,bisect). But my function is multivariate, so these will not work for me.
Harald Schilly
> t I don't
> want to find the minimum/maximum of my function, I want to find the root.
Very generally speaking, you can always find the root by minimization,
if you square your function (simply because no negative values are
possible)!
H
Ashley DaSilva
Ahh, of course. I should have thought of this :) I will try it, thanks for the suggestion.
Ashley
Warren Weckesser
quadratic convergence of Newton's method, which is at the heart (in some
form or another) of most good root-finding algorithms.
Warren
--
Warren Weckesser
Enthought, Inc.
515 Congress Avenue, Suite 2100
Austin, TX 78701
512-536-1057
Ashley DaSilva
Do you know a way to add constraints to fsolve, or some other root finding technique? If there's no other option, I'll have to go with Harald's suggestion, even if it is slow to converge.
Also, does anyone know what the input format is for these minimization techniques (fmin_l_bfgs_b, fmin_tnc, fmin_cobyla), I tried
def f(x):
Ashley DaSilva
Warren, thanks for your input.
Do you know a way to add constraints to fsolve, or some other root finding technique? If there's no other option, I'll have to go with Harald's suggestion, even if it is slow to converge.
Also, does anyone know what the input format is for these minimization techniques (fmin_l_bfgs_b, fmin_tnc, fmin_cobyla), I tried:
return [ n-P(x), n*u-Q(x), n-N(x)]
where P(x), Q(x), N(x) are defined functions and n,u are global constants. I get the error
TypeError: unsupported operand type(s) for -: 'list' and 'list'
which I think is referring to the numerical approximation to the gradient, when it tries to subtract f(x)-f0
Ashley
Robert Kern
> Sorry I pressed send by accident.
>
>> Warren, thanks for your input.
>>
>> Do you know a way to add constraints to fsolve, or some other root finding
>> technique? If there's no other option, I'll have to go with Harald's
>> suggestion, even if it is slow to converge.
>>
>> Also, does anyone know what the input format is for these minimization
>> techniques (fmin_l_bfgs_b, fmin_tnc, fmin_cobyla), I tried:
>
> def f(x):
> return [ n-P(x), n*u-Q(x), n-N(x)]
>
> where P(x), Q(x), N(x) are defined functions and n,u are global constants. I
> get the error
> TypeError: unsupported operand type(s) for -: 'list' and 'list'
> which I think is referring to the numerical approximation to the gradient,
> when it tries to subtract f(x)-f0
Always copy-and-paste the traceback, not just the final message. For
the fmin_cobyla constraints, you don't pass a function that returns a
list. You pass a list of functions.
--
Robert Kern
"I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth."
-- Umberto Eco
Ashley DaSilva
I am using fmin_tnc, where I have: import fmin_tnc as fmin
Traceback (most recent call last):
File "hf_improved.py", line 190, in <module>
solution=fmin(f2,x0,args=(eE,),bounds=[(0,1),(0,None),(0,None)],approx_grad=True)
File "/usr/lib/python2.5/site-packages/scipy/optimize/tnc.py", line 246, in fmin_tnc
fmin, ftol, xtol, pgtol, rescale)
File "/usr/lib/python2.5/site-packages/scipy/optimize/tnc.py", line 200, in func_and_grad
g = approx_fprime(x, func, epsilon, *args)
File "/usr/lib/python2.5/site-packages/scipy/optimize/optimize.py", line 617, in approx_fprime
grad[k] = (f(*((xk+ei,)+args)) - f0)/epsilon
TypeError: unsupported operand type(s) for -: 'list' and 'list'
def f2(v,E):
u=v[0]
return [(E-P(v))**2, (E*u-Q(v))**2, (1-N(v))**2]
where P, Q, N, are defined functions.
Ashley
Robert Kern
> Sorry, I will send the traceback.
> I am using fmin_tnc, where I have: import fmin_tnc as fmin
>
> Traceback (most recent call last):
> File "hf_improved.py", line 190, in <module>
>
> solution=fmin(f2,x0,args=(eE,),bounds=[(0,1),(0,None),(0,None)],approx_grad=True)
> File "/usr/lib/python2.5/site-packages/scipy/optimize/tnc.py", line 246,
> in fmin_tnc
> fmin, ftol, xtol, pgtol, rescale)
> File "/usr/lib/python2.5/site-packages/scipy/optimize/tnc.py", line 200,
> in func_and_grad
> g = approx_fprime(x, func, epsilon, *args)
> File "/usr/lib/python2.5/site-packages/scipy/optimize/optimize.py", line
> 617, in approx_fprime
> grad[k] = (f(*((xk+ei,)+args)) - f0)/epsilon
> TypeError: unsupported operand type(s) for -: 'list' and 'list'
>
>
> Other information that might be useful:
> def f2(v,E):
> u=v[0]
> return [(E-P(v))**2, (E*u-Q(v))**2, (1-N(v))**2]
>
> where P, Q, N, are defined functions.
Oh, you mean for the function itself. Return an array, not a list.
Sebastian Walter
systems by reformulating it as a NLP.
But if the function F(x) =0 is uniquely solvable and f is twice
continuously differentiable, then
f(x) := |F(x)|^2 is also C^2 and should have a strict local minimum.
Then Newton's method should locally converge quadratically.
Am I missing something here?
Ernest Adrogué
I ran into the same problem.
One thing you can try is add a couple of lines to your function
so that it returns a constant value (which is not a solution to your
problem) when (1-x**2) < 0, instead of the actual calculation.
In my case it did the trick.
Ernest
Warren Weckesser
> well, I agree that is not such a clever idea to solve nonlinear
> systems by reformulating it as a NLP.
> But if the function F(x) =0 is uniquely solvable and f is twice
> continuously differentiable, then
> f(x) := |F(x)|^2 is also C^2 and should have a strict local minimum.
>
Yes.
> Then Newton's method should locally converge quadratically.
>
No, because the derivative of the function being minimized is zero at
the root. In this case the convergence of Newton's method is only linear.
Warren
Sebastian Walter
definition of a KKT point.
I think you confuse Newton's method to solve nonlinear systems and
Newton's method for nonlinear optimization.
Yes, if you have a f(x) = f'(x) = 0 then you destroy the local
quadratic convergence of Newton's method.
But in nonlinear programming you do:
f(x) = | F(x)|^2
x_* = argmin_x f(x)
Then you define the first order necessary conditions for optimality:
0 = df(x)/dx
Therefore
0 =G(x) := 2 dF/dx F
I.e. you get a nonlinear system that you can solve with Newton's method.
If you differentiate once more you get
H = 2 (d^2/dx F) F + 2 (d/dx F)^2
and use the update rule
x_+ = x - H(x)^-1 G(x)
in our case H is symmetric positive definite and therefore Newton
should converge quadratically.
Sebastian Walter
x_* = argmin_x f(x)
subject to g(x) <= 0
then you can try to do:
x_* = argmin_x f(x) + \rho (max(g(x), 0))^p
where e.g. p =2, and make \rho large
problem: badly conditioned for \rho very large.
But that's still much better than adding a constant when g(x) > 0!
2009/8/4 Ernest Adrogué <eadr...@gmx.net>:
Ernest Adrogué
> what you can try to do is using a penalty method:
>
> x_* = argmin_x f(x)
> subject to g(x) <= 0
>
> then you can try to do:
>
> x_* = argmin_x f(x) + \rho (max(g(x), 0))^p
> where e.g. p =2, and make \rho large
>
> problem: badly conditioned for \rho very large.
>
> But that's still much better than adding a constant when g(x) > 0!
But... this won't prevent the function from being evaluated
outside of its domain, which is the real issue here, will it?
Sebastian Walter
Then, interior points algorithms as implemented in IPOPT should be
method of choice.
They have the nice property to stay feasible (granted you can find a
feasible starting value ...)
2009/8/4 Ernest Adrogué <eadr...@gmx.net>:
Warren Weckesser
> well, of course the derivative is zero at the root, that's the
> definition of a KKT point.
> I think you confuse Newton's method to solve nonlinear systems and
> Newton's method for nonlinear optimization.
>
Yes. :)
Ernest Adrogué
> uhmm, you are completely right... penalty methods are useless in that case.
> Then, interior points algorithms as implemented in IPOPT should be
> method of choice.
i would like to try ipopt, unfortunately it depends on third-party
libraries with draconian licenses. only mumps seems to actually be
"free", and yet in order to get it you have to fill a form and wait
until they send you the software by mail!!
Sebastian Walter
restrictive.
I use the ma27 (I think it is called) as academic user. It's not as
bad as you wrote above:
You register, you get a confirmation mail with a password, then you
can log in an download the code as a couple of fortran source code
files.
So, not such a big deal.
2009/8/4 Ernest Adrogué <eadr...@gmx.net>:
Ernest Adrogué
> Yeah, the licences of the linear solvers used in IPOPT are really too
> restrictive.
>
> I use the ma27 (I think it is called) as academic user. It's not as
> bad as you wrote above:
> You register, you get a confirmation mail with a password, then you
> can log in an download the code as a couple of fortran source code
> files.
> So, not such a big deal.
In all fairness, a little while later I received an e-mail with the link to
download MUMPS.
The README file states that the software is free of charge and in the public
domain. Why they don't put the download link directly on the web I don't know.
Sebastian Walter
of interfaces are defined.
Some collegues of mine use UMFPACK to solve sparse linear systems, as
far as i know.
2009/8/6 Ernest Adrogué <eadr...@gmx.net>:
Ernest Adrogué
> well, is MUMPS directly supported by IPOPT? I thought, only a couple
> of interfaces are defined.
> Some collegues of mine use UMFPACK to solve sparse linear systems, as
> far as i know.
>From the ipopt website:
"""
Currently, the following linear solvers can be used:
* MA27 from the Harwell Subroutine Library
(see http://www.cse.clrc.ac.uk/nag/hsl/).
* MA57 from the Harwell Subroutine Library
(see http://www.cse.clrc.ac.uk/nag/hsl/).
* MUMPS (MUltifrontal Massively Parallel sparse direct Solver)
(see http://graal.ens-lyon.fr/MUMPS/)
* The Parallel Sparse Direct Solver (PARDISO)
(see http://www.computational.unibas.ch/cs/scicomp/software/pardiso/).
Note: The Pardiso version in Intel's MKL library does not yet support the features necessary for IPOPT.
* The Watson Sparse Matrix Package (WSMP)
(see http://www-users.cs.umn.edu/~agupta/wsmp.html)
You should include at least one of the linear solvers above in order to run IPOPT, and if you want to be able to switch easily between different alternatives, you can compile IPOPT with all of them.
"""
So, I guess, yes, it is supported. Can't tell you for certain
because the Debian distribution that I'm using lacks a certain header
file from the metis/scotch library which is required to compile MUMPS,
so haven't had a chance to try it out yet.