Solving multivariate nonlinear system of equations

[Fmwyso null]

To finish some of my research, I'd like to use sympy to compute a multivariate nonlinear system of equations. This system is under-determined but I can guarantee it is solvable. 

The simplest case I would like to solve has 21 unknowns and 12 constraint equations; see following for all constraint equations. Ideally, I would like to have a solution set for the 21 unknowns (I'll call them "a" - "u") that has a small norm2. 

1 = a*j*p + d*l*q + g*n*r
0 = b*j*p + e*l*q + h*n*r
0 = c*j*p + f*l*q + i*n*r
0 = a*k*p + d*m*q + g*o*r
1 = b*k*p + e*m*q + h*o*r
0 = c*k*p + f*m*q + i*o*r
0 = a*j*s + d*l*t + g*n*u
1 = b*j*s + e*l*t + h*n*u
0 = c*j*s + f*l*t + i*n*u
0 = a*k*s + d*m*t + g*o*u
0 = b*k*s + e*m*t + h*o*u
1 = c*k*s + f*m*t + i*o*u

Currently, I am trying to use sympy's "solve" for all of these equations and it hasn't even completed computation in 24 hours.

Is my problem simply too complicated or am I using the wrong function for this? I appreciate any suggestions you may have.

Thanks for your time!

[Denis Akhiyarov]

Are you trying to find set of symbolic solutions by elimination or a unique numerical solution by optimization, since it is indeterminate?

[Fmwyso null]

Sorry for the vagueness of "solution set", I should have specified that all I require is a unique numerical solution. 

If optimization can provide this optimal solution (with minimized norm2), that would be great.

Otherwise, if Sympy could provide a set of symbolic solutions, I could construct my own minimized norm2 solution.

[Denis Akhiyarov]

something like this:

~ $ ptipython
Python 3.5.1 |Anaconda 4.0.0 (64-bit)| (default, Dec  7 2015, 11:16:01)
Type "copyright", "credits" or "license" for more information.


IPython 4.1.2 -- An enhanced Interactive Python.
?         -> Introduction and overview of IPython's features.
%quickref -> Quick reference.
help      -> Python's own help system.
object?   -> Details about 'object', use 'object??' for extra details.


In [1]: def equations(p):
   ...:     (a, b, c, d, e, f,
   ...:     g, h, i, j, k, l,
   ...:     m, n, o, p, q, r, s, t, u)=p
   ...:     return (1 - a*j*p + d*l*q + g*n*r,
   ...:              0 - b*j*p + e*l*q + h*n*r,
   ...:              0 - c*j*p + f*l*q + i*n*r,
   ...:              0 - a*k*p + d*m*q + g*o*r,
   ...:              1 - b*k*p + e*m*q + h*o*r,
   ...:              0 - c*k*p + f*m*q + i*o*r,
   ...:              0 - a*j*s + d*l*t + g*n*u,
   ...:              1 - b*j*s + e*l*t + h*n*u,
   ...:              0 - c*j*s + f*l*t + i*n*u,
   ...:              0 - a*k*s + d*m*t + g*o*u,
   ...:              0 - b*k*s + e*m*t + h*o*u,
   ...:              1 - c*k*s + f*m*t + i*o*u)


In [2]: from scipy.optimize import fmin


In [3]: import numpy as np


In [4]: def f(p):return np.abs(np.sum(np.array(equations(p))**2)-0)


In [5]: fmin(f,(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1))
Warning: Maximum number of function evaluations has been exceeded.
Out[5]:
array([-0.89996132,  3.24921453,  2.05127943, -2.03527783,  1.69544993,
        1.40438226,  1.92348097,  1.06625688,  0.24227183,  0.08901109,
        0.09085688,  2.91819752,  0.6803978 ,  4.09234832,  1.90642669,
        2.95941855,  0.17154083, -0.02346012,  1.44161738, -0.02367669,
       -0.03451558])


In [6]: from scipy.optimize import fmin_l_bfgs_b


In [7]: fmin_l_bfgs_b(f,(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1))
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-7-b040792986e3> in <module>()
----> 1 fmin_l_bfgs_b(f,(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1))


/home/dta/anaconda3/lib/python3.5/site-packages/scipy/optimize/lbfgsb.py in fmin_l_bfgs_b(func, x0, fprime, args, approx_grad, bounds, m, factr, pgtol, epsilon, iprint, maxfun, maxiter, disp, callback, maxls)
    191
    192     res = _minimize_lbfgsb(fun, x0, args=args, jac=jac, bounds=bounds,
--> 193                            **opts)
    194     d = {'grad': res['jac'],
    195          'task': res['message'],


/home/dta/anaconda3/lib/python3.5/site-packages/scipy/optimize/lbfgsb.py in _minimize_lbfgsb(fun, x0, args, jac, bounds, disp, maxcor, ftol, gtol, eps, maxfun, maxiter, iprint, callback, maxls, **unknown_options)
    328                 # minimization routine wants f and g at the current x
    329                 # Overwrite f and g:
--> 330                 f, g = func_and_grad(x)
    331         elif task_str.startswith(b'NEW_X'):
    332             # new iteration


/home/dta/anaconda3/lib/python3.5/site-packages/scipy/optimize/lbfgsb.py in func_and_grad(x)
    276     else:
    277         def func_and_grad(x):
--> 278             f = fun(x, *args)
    279             g = jac(x, *args)
    280             return f, g


/home/dta/anaconda3/lib/python3.5/site-packages/scipy/optimize/optimize.py in function_wrapper(*wrapper_args)
    287     def function_wrapper(*wrapper_args):
    288         ncalls[0] += 1
--> 289         return function(*(wrapper_args + args))
    290
    291     return ncalls, function_wrapper


/home/dta/anaconda3/lib/python3.5/site-packages/scipy/optimize/optimize.py in __call__(self, x, *args)
     62         self.x = numpy.asarray(x).copy()
     63         fg = self.fun(x, *args)
---> 64         self.jac = fg[1]
     65         return fg[0]
     66


IndexError: invalid index to scalar variable.


In [8]: from scipy.optimize import fmin_bfgs


In [9]: fmin_bfgs(f,(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1))
Optimization terminated successfully.
         Current function value: 2.500000
         Iterations: 12
         Function evaluations: 483
         Gradient evaluations: 21
Out[9]:
array([ 0.99750085,  1.13125434,  0.99750097,  0.84397729,  0.75589231,
        0.84398644,  0.8439863 ,  0.75589223,  0.84397754,  1.07736067,
        1.07736081,  0.73335887,  0.73335281,  0.73335274,  0.73335887,
        1.07736068,  0.73335893,  0.73335289,  1.07736063,  0.73335281,
        0.73335887])

[Denis Akhiyarov]

with differential global optimization one of the roots is found:

In [10]: differential_evolution(f,((-2,2), (-2,2), (-2,2), (-2,2), (-2,2), (-2,2
   ....: ), (-2,2), (-2,2), (-2,2), (-2,2), (-2,2), (-2,2), (-2,2), (-2,2), (-2,
   ....: 2), (-2,2),(-2,2), (-2,2),(-2,2), (-2,2), (-2,2)))
Out[10]:
     fun: 1.8444985066496224e-09
     jac: array([  8.23037829e-06,  -3.25546909e-05,  -1.20137119e-05,
        -3.40010940e-06,  -3.18740293e-05,  -3.19651128e-05,
         1.18561859e-05,  -5.92874802e-05,  -4.87322296e-05,
         1.17806110e-06,   2.53454759e-05,   6.31250916e-06,
        -1.75807827e-06,   1.92369061e-06,   6.78304593e-05,
         1.76279553e-05,   1.43592260e-05,   4.05940515e-05,
        -4.35746205e-05,  -2.54559468e-05,   3.37541106e-05])
 message: 'Maximum number of iterations has been exceeded.'
    nfev: 315667
     nit: 1000
 success: False
       x: array([ 0.62224365,  1.46881062, -1.47343122, -0.85611247,  0.09195341,
        0.423449  , -0.20369088, -0.76820833, -0.41410515,  0.66288976,
       -1.01724249, -1.97129718,  0.63188554,  0.69403993,  0.91437698,
       -0.2956973 , -0.60242029,  0.74482681,  0.45377653,  0.19310465,
        0.98123914])


In [11]: res=_


In [12]: res.x
Out[12]:
array([ 0.62224365,  1.46881062, -1.47343122, -0.85611247,  0.09195341,
        0.423449  , -0.20369088, -0.76820833, -0.41410515,  0.66288976,
       -1.01724249, -1.97129718,  0.63188554,  0.69403993,  0.91437698,
       -0.2956973 , -0.60242029,  0.74482681,  0.45377653,  0.19310465,
        0.98123914])

In [13]: equations(res.x)
Out[13]:
(-2.8165037928157277e-06,
 -9.5617945984893815e-06,
 -1.4916310316442916e-05,
 -4.1531045657239307e-06,
 -5.3874567714773391e-06,
 -1.5715540758742819e-05,
 2.8150894887946087e-06,
 7.9860279372789833e-06,
 7.2868002708448287e-06,
 9.2345080527356238e-06,
 -2.9513712621720423e-05,
 -1.2176049990486604e-05)