Solving non linear system containing a sum
<p>I have several functions which comes from sympy.lambdify:</p>

<pre><code>f_1 = sym.lambdify((z, m_1, m_2, s_1, s_2), expression_1, modules=['numpy', 'sympy'])
f_2 = sym.lambdify((z, m_1, m_2, s_1, s_2), expression_2, modules=['numpy', 'sympy'])
f_3 = sym.lambdify((z, m_1, m_2, s_1, s_2), expression_3, modules=['numpy', 'sympy'])
f_4 = sym.lambdify((z, m_1, m_2, s_1, s_2), expression_4, modules=['numpy', 'sympy'])
</code></pre>

<p>where <code>m_1</code>, <code>m_2</code>, <code>s_1</code>, <code>s_2</code> are scalars and <code>z</code> is a known 1D array (not necessary the same for each <code>f_i</code>). The output of each <code>f_i</code> is a scalar.</p>

<p>I would like to (numerically) find <code>m_1</code>, <code>m_2</code>, <code>s_1</code>, <code>s_2</code> so that,</p>

<pre><code>sum(f_1(z_i, m_1, m_2, s_1, s_2)) = 0
sum(f_2(z_i, m_1, m_2, s_1, s_2)) = 0
sum(f_3(z_i, m_1, m_2, s_1, s_2)) = 0
sum(f_4(z_i, m_1, m_2, s_1, s_2)) = 0
</code></pre>

<p>The sum being on i.</p>

<p>Using scipy.optimize, I do not know how to implement it (using fsolve or root).</p>



<p>To drive 4 functions F1, F2, F3, F4 all to 0,
minimize the sum of squares F1^2 + F2^2 + F3^2 + F4^2:
if the sum is small, each F must be small too.
(For example, sum &lt; 0.000001 &rArr; each |F| &lt; 0.001 .)
Use <a href="http://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html#least-squares-minimization-least-squares" rel="nofollow">scipy.optimize.least_squares</a>
to do this:</p>

<pre><code>import numpy as np
from scipy.optimize import least_squares
    # http://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html#least-squares-minimization-least-squares
    # http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html

# minimize F1**2 + F2**2 + F3**2 + F4**2 --
starting_guess = ...
ret = least_squares( [F1, F2, F3, F4], x0=starting_guess, max_nfev=20, verbose=2 )
xmin = ret.x
</code></pre>

<p>(Do NOT square the F s -- <code>least_squares</code> does that for you.)</p>

<p>In your case, you want to move the sum outside <code>lambdify</code>:</p>

<pre><code>f1 = sym.lambdify((z, m1, m2, s1, s2) ... ))
Z1 = numpy array whose rows are all the z_i for f1

def F1( m1, m2, s1, s2 ):
    sum = np.sum([ f( z, m1, m2, s1, s2 ) for z in Z1 ])
    print "F1: %-10.3g at %-10.3g %-10.3g %-10.3g %-10.3g" % (sum, m1, m2, s1, s2)
    return sum

# check a few values in ipython --
F1( 0, 0, 0, 0 )
for x in np.eye( 4 ):
    F1( *x )
</code></pre>

<p>When that works, F2 F3 F4 are similar. Check them too. Then optimize:</p>

<pre><code>starting_guess = ...
ret = least_squares( [F1, F2, F3, F4], x0=starting_guess, max_nfev=10, verbose=2 )
xmin = ret.x
</code></pre>

<p>or more iterations, or use <code>xtol</code> <code>ftol</code>.<br>
<code>least_squares</code> has lots of options, see
the <a href="http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html" rel="nofollow">full doc</a> .</p>

<p>There are fancier ways of doing the above for any number of functions, not just 4,
of any number of variables: more general, more obscure.</p>

<p>(Avuncular advice for optimization:</p>

<ul>
<li>start small</li>
<li>check every step, with print statements, and interactively in ipython.</li>
</ul>

python optimization scipy