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.

enter image description here

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:

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.
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)

Do you know that a root exists? If you do, have you double- or triple-checked that your equations are typed in correctly? — jared, Jul 16 '23 at 18:45
Can you give some hint as to how these equations are derived? I'm in the middle of vectorising them, and it would be helpful to have more meaningful symbols. — Reinderien, Jul 17 '23 at 00:01
@Reinderien I was able to easily vectorize them by taking the sympy equations, running them through `lambdify` (with variables in the order of `(sigma_SU, sigma_MU, sigma_R)`, but that order doesn't matter as long as you're consistent), and then having the `equations(z)` function just return `np.array([f1_f(*z), f2_f(*z), f3_f(*z)])` (where `f#_f` is the lambdify result of `f#`). — jared, Jul 17 '23 at 00:16
Sure, though I don't trust the output from such a method. It wasn't too bad to do it by hand, it just doesn't have meaningful variable names. — Reinderien, Jul 17 '23 at 00:25
@Reinderien Why not? I compared the results of the automatically-vectorized equations to the sympy `.subs` results and they matched. I also tested your solution with the automatically-vectorized equations and got the same result (to machine precision). — jared, Jul 17 '23 at 00:44
I would consider it useful as a second method of regression test, but: does it provide a multi-stage expression for variable reuse? Partially I'm optimizing for legibility. — Reinderien, Jul 17 '23 at 00:47
@Reinderien I cannot argue that it returns the most efficient functions, but it is efficient if you consider the time it takes to rewrite what was given. — jared, Jul 17 '23 at 01:03

Reinderien · Answer 1 · 2023-07-17T01:31:28.417

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 ]

`scipy.optimize.root(equations, np.zeros(3))` also seems to succeed, though it does not when using the guess from the question. — jared, Jul 17 '23 at 01:03
Interesting. I think it is machine-dependent since I don't see any significant difference in the results when I use your code and use the zero vector as the initial guess (i.e. on my machine the function evaluation stays at the same order of 1e-13). — jared, Jul 17 '23 at 01:10
Hi! Thanx for the help, I have a new approach now, but all the root values still need to be positive, otherwise it makes no sense in the real world. I'll try to adjust the initial estimate. — nparadina, Jul 23 '23 at 15:49

score 1 · Answer 2 · answered Jul 17 '23 at 01:18

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.

Three variables exponential equations in Python - solver options

2 Answers2