1

Converting Sympy derived equations to be compatible with scipy.optimize.newton for multivariate problem

I referred to Error using 'exp' in sympy -TypeError and Attribute Error is displayed. I found out that using lambdify will make the symbolic equations compatible with other packages such as numpy and scipy.

This is a continuation of my previous post: Error using 'exp' in sympy -TypeError and Attribute Error is displayed 

import numpy as np
import sympy as sym
import scipy.optimize
from sympy import symbols, diff, exp, log, power
from sympy.utilities.lambdify import lambdify

data = [3, 33, 146, 227, 342, 351, 353, 444, 556, 571, 709, 759, 836, 860, 968, 1056, 1726, 1846, 1872, 1986, 2311, 2366, 2608, 2676, 3098, 3278, 3288, 4434, 5034, 5049, 5085, 5089, 5089, 5097, 5324, 5389,5565, 5623, 6080, 6380, 6477, 6740, 7192, 7447, 7644, 7837, 7843, 7922, 8738, 10089, 10237, 10258, 10491, 10625, 10982, 11175, 11411, 11442, 11811, 12559, 12559, 12791, 13121, 13486, 14708, 15251, 15261, 15277, 15806, 16185, 16229, 16358, 17168, 17458, 17758, 18287, 18568, 18728, 19556, 20567, 21012, 21308, 23063, 24127, 25910, 26770, 27753, 28460, 28493, 29361, 30085, 32408, 35338, 36799, 37642, 37654, 37915, 39715, 40580, 42015, 42045, 42188, 42296, 42296, 45406, 46653, 47596, 48296, 49171, 49416, 50145, 52042, 52489, 52875, 53321, 53443, 54433, 55381, 56463, 56485, 56560, 57042, 62551, 62651, 62661, 63732, 64103, 64893, 71043, 74364, 75409, 76057, 81542, 82702, 84566, 88682]
n = len(data)
tn = data[n-1]


b, c = sym.symbols('b c', real=True)

f = -(-n +sum([sym.log(b*c*(num**(c-1))*sym.exp(-b*(num**c))) for num in data]))

bh = lambdify((b,c),diff(f,b),"numpy")
ch = lambdify((b,c),diff(f,c),"numpy")

sol = scipy.optimize.newton([bh,ch],(0.00404,1.0))
print(sol)

I can't seem to get the Newton's method working. Any information or resources is appreciated.

SuperKogito
  • 2,998
  • 3
  • 16
  • 37
Lucky
  • 347
  • 6
  • 15
  • `newton` works for scalar functions. For multidimensional functions see [here](https://docs.scipy.org/doc/scipy/reference/optimize.html#id2) – Tarifazo May 21 '19 at 13:03

2 Answers2

1

According to the scipy.optimize.newton your x0 should be a scalar and not an array or tuple (which is what you are passing to scipy.optimize.newton() in your code). Fortunately, scipy.optimize.fsolve can also find the zero position and accepts arrays as x0.

From what I see, you are making life harder by using sympy, when you can simply use numpy, which speeds your code (you can avoid the for loop). The following example shows you how it can be used (hope it helps): 

import warnings 
import numpy as np
import scipy.optimize as opt
# filter warnings
warnings.filterwarnings("ignore")

data = np.array([3, 33, 146, 227, 342, 351, 353, 444, 556, 571, 709, 759, 836, 860, 968,
        1056, 1726, 1846, 1872, 1986, 2311, 2366, 2608, 2676, 3098, 3278, 3288,
        4434, 5034, 5049, 5085, 5089, 5089, 5097, 5324, 5389,5565, 5623, 6080, 
        6380, 6477, 6740, 7192, 7447, 7644, 7837, 7843, 7922, 8738, 10089,
        10237, 10258, 10491, 10625, 10982, 11175, 11411, 11442, 11811, 12559, 
        12559, 12791, 13121, 13486, 14708, 15251, 15261, 15277, 15806, 16185,
        16229, 16358, 17168, 17458, 17758, 18287, 18568, 18728, 19556, 20567,
        21012, 21308, 23063, 24127, 25910, 26770, 27753, 28460, 28493, 29361,
        30085, 32408, 35338, 36799, 37642, 37654, 37915, 39715, 40580, 42015,
        42045, 42188, 42296, 42296, 45406, 46653, 47596, 48296, 49171, 49416,
        50145, 2042,  52489, 52875, 53321, 53443, 54433, 55381, 56463, 56485,
        56560, 57042, 62551, 62651, 62661, 63732, 64103, 64893, 71043, 74364,
        75409, 76057, 81542, 82702, 84566, 88682])

n    = len(data)
tn   = data[n-1]

f        = lambda b, c: n - np.sum(np.log(b * c * data**(c-1) * np.exp(-b * data**c)))
obj_func = lambda x   : np.array([f(x[0], x[1]) - x[0], f(x[0], x[1]) - x[1]])


x0  = np.array([0.00404, 1.0])
sol = opt.fsolve(obj_func, x0 = x0)
print('b, c = ', sol)

Output: 

b, c =  [0.01737803 0.4300348 ]
SuperKogito
  • 2,998
  • 3
  • 16
  • 37
  • Thank you. For the objective function, I have specified in my code, it is easier to use fsolve or root. But, I need to use Newton's method for much complicated objective function so I am trying to get this simple function working before moving on to a complex one, – Lucky May 21 '19 at 14:48
  • Well the problem is that `scipy.optimize.newton` does not support multidimensional functions so you will need to consider finding a different approach. That or solve for each variable at a time using a loop but that would be an ugly slow solution imo. Good luck. – SuperKogito May 21 '19 at 14:52
0

Instead of defining symbols, we can use the 'DeferredVector' in sympy to define the symbols. Instead of 

b, c = sym.symbols('b c', real=True)
f = -(-n +sum([sym.log(b*c*(num**(c-1))*sym.exp(-b*(num**c))) for num in data]))
bh = lambdify((b,c),diff(f,b),"numpy")
ch = lambdify((b,c),diff(f,c),"numpy"

do the following: x = DeferredVector('x') f = -(-n +sum([sym.log(x[0]x[1](num**(c-1)) for num in data]))

Now use diff w.r.t x[0] and x[1] and run the scipy.optimize.newton().

Lucky
  • 347
  • 6
  • 15