# Vectorised root finding in Python

I have an array of size (254, 80) which I am trying to use Scipy's fsolve on. I have found that the speed of using fsolve on a vector is quicker than it is in a for loop but only for vectors upto about 100 values long. After this, the speed quickly drops off and becomes very slow, sometimes completely stopping.

I'm currently looping through one dimension of the array and using a vectorised fsolve on the smaller dimension but it's still taking longer than I would expect/like.

Does anyone have a good work around for this or a know of a similar function which will be happy handling a vector of a larger size? I have tried Scipy's Newton Krylov function which, upon testing, is really fast for simple functions but it's slower than fsolve on my function. - Will this be faster if I use the analytical jacobian? If so, how do I implement it?

I wondered if it is a problem with my initial solutions but the difference is in the range of ~0.001 - 0.8 which doesn't seem too bad. Although perhaps it is this?

Here's the current code:

for i in range(array.shape[0]):
f = lambda y: a[i] - m[i]*y - md[i]*(( y**4 + 2*(y**2)*np.cos(Thetas[i,:]) )**0.25)
ystar[i,:] = fsolve(f, y0[i])


(Variables a, m and md are all the same size i.e. (254,80))

Digging in to this further, I have found that a function such as

f = lambda y: y*np.tanh(y) - a0/(m**2)


is faster to solve than

f = lambda y: (m**2)y*np.tanh(y) - a0


where m and a0 are large 2D np arrays.

Can anyone explain why this is?

• Finding the zeros over one of the dimensions is not the same as finding the zeros over the whole array. In one case the you might have variables that are not zero. Regarding your question at the end, it might be faster to solve one case over the other because of the scaling, that leads to better behaved functions Dec 2 '17 at 15:04
• I want a zero for each entry so it doesn't matter what order it is done in - each value is independent Dec 3 '17 at 14:33