Three variables exponential equations in Python - solver options

Question

Asking here is my last resort, I do not know where else to look or what my next option would be, apart from changing languages.

I have a set of 3 heavy exponential equations, with 3 variables. All have the unknown variable in the exponent as well as outside.

I have tried sympy's solve, nsolve, and nonlinsolve as well as scipy's fsolve. Scipy's fsolve spits up one solution, but one containing negative values, which are not allowed, since what I am trying to calculate is the transitional intensity of a person getting sick, and it can not be negative. My guess is based on empirical data, to apply any other random guess doesn't have a sound ground behind the decision. Sympy's solve returns a [] or just keeps running, depending on if called from within a main or standalone (as copied later), nsolve spits out the error below, and nonlinsolve runs endlessly without returning anything.

The papers I have followed on the topic list MATLAB as the chosen tool, however, I would much rather stay in Python, as I am very rusty when it comes to MATLAB, and I have written (however badly, since I am not a real programmer, just a person using Python for a paper) a bunch of code to get to this unhappy point, I am now stuck at.

If it helps, this is what the equations are, and the example of a call to solve I am making (extracted from the other parts of the code):

import sympy as sm
from sympy import nsolve
from sympy import solve

sigma_R,sigma_SU,sigma_MU=sm.symbols("sigma_R,sigma_SU,sigma_MU", real=True)

f1 = -0.262617004746481*sigma_SU*sm.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/((-sigma_MU -   sigma_R - sigma_SU + 0.000931334226877499)*(0.260884896827677*sigma_MU*sm.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/(-sigma_MU - sigma_R - sigma_SU + 0.00170434902663555) + 0.256617389979587*sigma_R*sm.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/(-sigma_MU - sigma_R - sigma_SU + 0.00363094118807797) + 0.262617004746481*sigma_SU*sm.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/(-sigma_MU - sigma_R - sigma_SU + 0.000931334226877499) + 0.719582172729324*sm.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU))) + 1.54

f3 = -0.260884896827677*sigma_MU*sm.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/((-sigma_MU - sigma_R - sigma_SU + 0.00170434902663555)*(0.260884896827677*sigma_MU*sm.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/(-sigma_MU - sigma_R - sigma_SU + 0.00170434902663555) + 0.256617389979587*sigma_R*sm.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/(-sigma_MU - sigma_R - sigma_SU + 0.00363094118807797) + 0.262617004746481*sigma_SU*sm.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/(-sigma_MU - sigma_R - sigma_SU + 0.000931334226877499) + 0.719582172729324*sm.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU))) + 0.3

f2 = -0.256617389979587*sigma_R*sm.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/((-sigma_MU - sigma_R - sigma_SU + 0.00363094118807797)*(0.260884896827677*sigma_MU*sm.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/(-sigma_MU - sigma_R - sigma_SU + 0.00170434902663555) + 0.256617389979587*sigma_R*sm.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/(-sigma_MU - sigma_R - sigma_SU + 0.00363094118807797) + 0.262617004746481*sigma_SU*sm.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/(-sigma_MU - sigma_R - sigma_SU + 0.000931334226877499) + 0.719582172729324*sm.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU))) + 0.37

sol=solve((f1,f2,f3),(sigma_R,sigma_SU,sigma_MU),simplify=False)
print(sol)

My ask here is:

1. Is there a way I could debugged, step-in etc. within any of the solve functions, the fact that it is a black-box doesn't help, since if I am somehow sending in faulty data, I get no exceptions.
2. If somebody knows of any approach I have not tried, please share; I am out of ideas. Even a math suggestion would work, but I do not see any means to simplify the equations themselves.

As said above, I have tried:

sol = solve((f1,f2,f3),(sigma_R,sigma_SU,sigma_MU),simplify=False)
sol = nsolve((f1, f2,f3), (sigma_R,sigma_SU,sigma_MU), (0.0000813,0.0000202,0.000384))
sol = sm.nonlinsolve([y0,y1,y2],[sigma_R,sigma_SU,sigma_MU])

And the part with fsolve:

from scipy.optimize import fsolve
import numpy as np
import math
import datetime as dt

def equations(z):
    sigma_SU=z[0]
    sigma_MU=z[1]
    sigma_R=z[2]
    f=np.empty(3)
    f[0] = -0.262617004746481*sigma_SU*np.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/((-sigma_MU - sigma_R - sigma_SU + 0.000931334226877499)*(0.260884896827677*sigma_MU*np.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/(-sigma_MU - sigma_R - sigma_SU + 0.00170434902663555) + 0.256617389979587*sigma_R*np.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/(-sigma_MU - sigma_R - sigma_SU + 0.00363094118807797) + 0.262617004746481*sigma_SU*np.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/(-sigma_MU - sigma_R - sigma_SU + 0.000931334226877499) + 0.719582172729324*np.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU))) + 1.54
    f[1] = -0.260884896827677*sigma_MU*np.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/((-sigma_MU - sigma_R - sigma_SU + 0.00170434902663555)*(0.260884896827677*sigma_MU*np.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/(-sigma_MU - sigma_R - sigma_SU + 0.00170434902663555) + 0.256617389979587*sigma_R*np.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/(-sigma_MU - sigma_R - sigma_SU + 0.00363094118807797) + 0.262617004746481*sigma_SU*np.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/(-sigma_MU - sigma_R - sigma_SU + 0.000931334226877499) + 0.719582172729324*np.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU))) + 0.3
    f[2] =  -0.256617389979587*sigma_R*np.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/((-sigma_MU - sigma_R - sigma_SU + 0.00363094118807797)*(0.260884896827677*sigma_MU*np.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/(-sigma_MU - sigma_R - sigma_SU + 0.00170434902663555) + 0.256617389979587*sigma_R*np.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/(-sigma_MU - sigma_R - sigma_SU + 0.00363094118807797) + 0.262617004746481*sigma_SU*np.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/(-sigma_MU - sigma_R - sigma_SU + 0.000931334226877499) + 0.719582172729324*np.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU))) + 0.37
    return f

myGuess=np.array([0.0000813,0.0000202,0.000384])
z=np.array([0,0,0])
counter=0

ct1 = dt.datetime.now()
print(ct1)

z=fsolve(equations, myGuess)

ct2=dt.datetime.now()
print(ct2-ct1)
print(z)
print(counter)

Answer 1

You need to rewrite your equations for sanity and performance to use Numpy properly. Then, switch from fsolve (MINPACK hybr*) to minimize. Initial vector (0, 0, 0) seems to converge more easily than the guess you have provided. Of the following four methods that successfully converge with no additional "help",

• Nelder-Mead,
• Powell,
• L-BFGS-B, and
• trust-constr,

trust-constr is the method that converges the most quickly on a solution that is "exact" to machine precision:

from scipy.optimize import minimize
import numpy as np

c1 = np.array((0.262617004746481000, 0.260884896827677000, 0.256617389979587000))
c2 = np.array((0.000931334226877499, 0.001704349026635550, 0.003630941188077970))
c3 = np.array((1.54, 0.30, 0.37))


def equations(z: np.ndarray) -> np.ndarray:
    # sigma_SU, sigma_MU, sigma_R = z

    inner1 = -z.sum()
    inner2 = np.exp(45*inner1)
    inner3 = c1 * z * inner2/(inner1 + c2)
    f = c3 - inner3/(inner3.sum() + 0.719582172729324 * inner2)
    return f


def lstsq_cost(z: np.ndarray) -> float:
    f = equations(z)
    return f.dot(f)


def regression_tests() -> None:
    expected = (1.48037098, 0.29461646, 0.33099203)
    actual = equations(np.array((8.13e-5, 2.02e-5, 3.84e-4)))
    assert np.allclose(expected, actual)

    expected = (0.93253452, 0.21240985, 0.06805596)
    actual = equations(np.array((-0.49717498, -0.00991807,  0.50892057)))
    assert np.allclose(expected, actual)


def optimize() -> None:
    result = minimize(
        fun=lstsq_cost, method='powell',
        x0=(0, 0, 0), tol=1e-9,
    )
    assert result.success, result.message
    np.set_printoptions(precision=2)
    print('Function at roots:', equations(result.x))
    np.set_printoptions(precision=20)
    print('Roots:', result.x)


if __name__ == '__main__':
    regression_tests()
    optimize()

Function at roots: [-0. -0.  0.]
Roots: [ 0.0032663289158950093   0.00011197597802975169 -0.0015103347529418544 ]

The brute-force outer search loop:

    for method in (
        'Nelder-Mead',
        'Powell',
        'CG',
        'BFGS',
        'Newton-CG',
        'L-BFGS-B',
        'TNC',
        'COBYLA',
        'SLSQP',
        'trust-constr',
        'dogleg',
        'trust-ncg',
        'trust-exact',
        'trust-krylov',
    ):
        try:
            result = minimize(
                fun=lstsq_cost, method=method,
                x0=(0, 0, 0), tol=1e-9,
            )
            assert result.success, result.message
            y = equations(result.x)
            if y.dot(y) > 0.1:
                continue
            np.set_printoptions(precision=2)
            print(method, 'converged in nfev =', result.nfev)
            print('Function at roots:', y)
            np.set_printoptions(precision=20)
            print('Roots:', result.x)
            print()
        except Exception:
            pass

Nelder-Mead converged in nfev = 285
Function at roots: [-9.10e-08 -1.46e-07  3.42e-08]
Roots: [ 0.0032675556912995915   0.00011212905180447633 -0.001511370899412161  ]

Powell converged in nfev = 1135
Function at roots: [ 6.51e-14 -7.53e-14  5.76e-14]
Roots: [ 0.0032675561103559787   0.00011212902858510153 -0.0015113712768471598 ]

L-BFGS-B converged in nfev = 188
Function at roots: [2.39e-05 2.28e-05 7.52e-06]
Roots: [ 0.0032676040558970083   0.00011211616767734874 -0.0015114200558563654 ]

trust-constr converged in nfev = 292
Function at roots: [-0. -0.  0.]
Roots: [ 0.0032663289158950093   0.00011197597802975169 -0.0015103347529418544 ]

Answer 2

Since I mentioned it in the comments, I will also add an alternative method to @Reinderien's answer for completeness (though their answer is better and more complete, so it should be accepted over mine).

If you don't want to manually rewrite the equations and still get equivalent results, you can use sympy's lambdify function to produce vectorized functions for each of the three functions you're working with (these functions won't be the most efficient in terms of runtime, but it is quick and dirty way to get vectorized function). You can use the "lambdified" functions to create a vectorized equations function. Solving for the roots can be done using scipy.optimize.root (fsolve is legacy syntax; you should be using this function instead) with the initial guess of np.zeros(3) (for this case, root doesn't converge for the initial guess from the question).

import sympy as sp
import numpy as np
from scipy.optimize import root

sigma_R, sigma_SU, sigma_MU = sp.symbols("sigma_R,sigma_SU,sigma_MU", real=True)
variables = (sigma_SU, sigma_MU, sigma_R)

# define f1, f2, and f3 using sympy as in the question

f1_f = sp.lambdify((variables), f1)
f2_f = sp.lambdify((variables), f2)
f3_f = sp.lambdify((variables), f3)

def equations(z):
    return np.array([f1_f(*z), f2_f(*z), f3_f(*z)])

sol = root(equations, np.zeros(3))
print(sol)
print(sol.x)

Output:

message: The solution converged.
success: True
 status: 1
    fun: [-6.737e-12 -2.788e-11  9.857e-12]
      x: [ 3.268e-03  1.121e-04 -1.511e-03]
   nfev: 51
   fjac: [[-8.232e-01  2.522e-01 -5.087e-01]
          [ 2.705e-01 -6.136e-01 -7.418e-01]
          [ 4.992e-01  7.483e-01 -4.369e-01]]
      r: [ 2.049e+03 -2.050e+03  2.379e+03  1.686e+03  3.289e+01
          -4.480e+02]
    qtf: [-2.165e-13  2.192e-10 -1.184e-09]

array([ 0.0032675561102804376 ,  0.00011212902858868879,
       -0.0015113712767817957 ])

---

P.S. If you're timing the speed of a function by taking the time before and after, use time.perf_counter() from the time library instead of the datetime method you used.

from time import perf_counter

start = perf_counter()
# code here
end = perf_counter()
print(end - start)

If you want to get really serious about timing, you can look into the timeit library.