Finding roots with scipy.optimize.root
Question: I am trying to find the root y of a function called f using Python.

Here is my code:

    
    
     def f(y):
        w,p1,p2,p3,p4,p5,p6,p7 = y[:8] 
        t1 = w - 0.500371726*(p1**0.92894164) - (-0.998515304)*((1-p1)**1.1376649)
        t2 = w - 8.095873128*(p2**0.92894164) - (-0.998515304)*((1-p2)**1.1376649)
        t3 = w - 220.2054377*(p3**0.92894164) - (-0.998515304)*((1-p3)**1.1376649)
        t4 = w - 12.52760758*(p4**0.92894164) - (-0.998515304)*((1-p4)**1.1376649)
        t5 = w - 8.710859537*(p5**0.92894164) - (-0.998515304)*((1-p5)**1.1376649)
        t6 = w - 36.66350261*(p6**0.92894164) - (-0.998515304)*((1-p6)**1.1376649)
        t7 = w - 3.922692207*(p7**0.92894164) - (-0.998515304)*((1-p7)**1.1376649)       
        t8 = p1 + p2 + p3 + p4 + p5 + p6 + p7 - 1
        return [t1,t2,t3,t4,t5,t6,t7,t8]
    
    
    x0 = np.array([-0.01,0.3,0.1,0.2,0.1,0.1,0.1,0.1])
    sol = scipy.optimize.root(f, x0, method='lm')
    print sol 
    print 'solution', sol.x
    print 'success', sol.success
    

Python does not find the root whatever the method I try in
scipy.optimize.root.

However there is one, I found it with the function fsolve in Matlab.

It is:

[-0.0622, 0.5855, 0.087, 0.0028, 0.0568, 0.0811, 0.0188, 0.1679].

When I specify x0 close to the root, the python algorithm converges. The
problem is that I have no idea a priori on the root to specify x0. In reality
I am solving many equations of this type.

I really want to use Python. Can anyone help me converge with python?


Answer: OK, after some fooling around, we focus on another aspect of good
optimization/root finding algorithms. In the comments above we went back and
forth around which method in scipy.optimize.root() to use. An equally
important question for near-bulletproof 'automatic' root finding is zeroing in
on good initial guesses. Often times, good initial guesses are, in fact, not
near the real answer at all. Instead, they need to be arranged so that they
will naturally lead the solver in the right direction.

In your particular case, your guesses were, in fact, sending the algorithm off
in strange directions.

My toy reconstruction of your problem is:

    
    
    import numpy as np
    import scipy as sp
    import scipy.optimize
    def f(y):
        w,p1,p2,p3,p4,p5,p6,p7 = y[:8]
        def t(p,w,a):       
            b = -0.998515304
            e1 = 0.92894164
            e2 = 1.1376649
            return(w-a*p**e1 - b*(1-p)**e2)
        t1 = t(p1,w,1.0)
        t2 = t(p2,w,4.0) 
        t3 = t(p3,w,16.0)
        t4 = t(p4,w,64.0)
        t5 = t(p5,w,256.0)
        t6 = t(p6,w,512.0)
        t7 = t(p7,w,1024.0)       
        t8 = p1 + p2 + p3 + p4 + p5 + p6 + p7 - 1.0
        return(t1,t2,t3,t4,t5,t6,t7,t8)
    guess = 0.0001
    x0 = np.array([-1000.0,guess,guess,guess,guess,guess,guess,guess])
    sol = sp.optimize.root(f, x0, method='lm')
    print('w=-1000: ', sol.x, sol.success,sol.nfev,np.sum(f(sol.x)))
    

Note that I did not use your specific prefactors (I wanted to broaden the
range I explored), although I kept your particular exponents on the p terms.

The real secret is in the initial guess, which I made the same for all p
terms. Having it a 0.1 or above bombed much of the time, since some terms want
to go one way and some the other. Reducing it to 0.01 worked well for this
problem. (I will note that the w term is very robust - varying from -1000. to
+1000. had no effect on the solution). Reducing the initial guess even further
has no effect on this particular problem, but it does no harm either. I would
keep it very small.

Yes, you know that at least some terms will be much larger. But, you are
putting the solver in a position where it can clearly and directly proceed
towards the real solution.

Good luck.

