# How can I use Python's SciPy and NumPy functions in Mathematica to find the minimum of a function?

I have written my code in Mathematica, final aim of this code is to find the minimum of a function.

I found that FindMinimum and NMinimize are very slow (see this "FindMinimum doesn't increase step size when necessary") when compared to scipy.optimize.minimize function in Python (some times, for my case it is slow).

I would like to run Python's scipy.optimize.minimize in the Mathematica notebook to find the minimum of the function written in the same Mathematica notebook.

I have gone through the links "Run Python with package NumPy in Mathematica", "Importing Mathematica defined function in python to do further computations. [Mathematica-Python interfacing]" and "Is it possible to call back to Mathematica when using Python through ExternalEvaluate?". These links are not much helpful. How can I write code to achieve this task?

For example, I would like to find minimum of this function f[x, y] ( for my case very big code in Mathematica notebook) and use this function in scipy.optimize.minimize and get the output in the same notebook.

f[x_, y_] := x^2 + y*Sin[x + y] + Sin[5*x];


I agree that Mathematica can find the solution for the above function within a fraction of a second.

How can I proceed?

spython = StartExternalSession["Python"]
sol = ExternalEvaluate[spython,"from scipy.optimize import minimize;b=(-10,10);bnds=(b,b);xin=[-10,-10];sol=minimize("<>ToString[f]<>",xini,bounds=bnds,method=SLSQP);sol"]


How can I pass arguments and how will Python will pass an argument to a Mathematica function to find the minimum?

• Your python code involves typo, this isn't the key issue, of course. Jun 14, 2021 at 15:08

I believe you're not in the correct direction, optimizing your Mathematica code should be more practical and easier. (You've already learned the numeric capability of Mathematica under your previous question, don't you? )

Anyway, here's my attempt to minimize your toy function f using Python in Mathematica. I know little about Python, feel free to point out if I've done something improper:

spython = StartExternalSession["Python"];

f[x_, y_] = x^2 + y Sin[x + y] + Sin[5 x];
(* helper function is defined,
because fun of scipy.optimize.minimize seems to only accept 1D array: *)
func[x_] := Evaluate@f[x[[1]], x[[2]]]; // Quiet

ExternalEvaluate[spython,
{"from wolframclient.evaluation import WolframLanguageSession",
"from wolframclient.language import wlexpr",
"session = WolframLanguageSession()",
"from wolframclient.language import Global",
"from scipy.optimize import minimize"}];

ExternalEvaluate[spython,
{StringTemplate["session.evaluate(wlexpr('1'))"][
func // Definition // InputForm // ToString],
"f= session.function(Global.func)"}];

rst = ExternalEvaluate[spython,
"b=(-10,10);bnds=(b,b);xin=[-10,-10];minimize(f,xin,bounds=bnds,method='SLSQP')"]; \
// AbsoluteTiming
(* {0.0694416, Null} *)

rst /@ {"fun", "x"} // Normal
(* {-0.365518, {0.887446, -0.457956}} *)


Notice Python hasn't done a good job in this case:

NMinimize[{f[x, y], -10 < x < 10, -10 < y < 10}, {x, y},
Method -> "DifferentialEvolution"] // AbsoluteTiming
(* {0.200323, {-10.1225, {x -> 0.905711, y -> 10.}}} *)


And I won't be too surprised if it's hard or even impossible to adapt the method above to more complicated cases.

• thanks for the solution. I should have given my actual code . I got my result in 5 minutes in python, where as in Mathematica it took be days together, This might be because the way I might have written the code in Mathematica. I will keep my actual code here as an edit. I felt linking Mathematica and python might be of some use in some cases. I'm very thankful to your solution. Jun 14, 2021 at 15:35
• @GummalaNavneeth I'd suggest posting your actual problem as a new question, if it's not too localized. Jun 14, 2021 at 15:50
• I will make it a general problem and post my localized problem also , so that it would be of use to others as well. Thanks a lot again Jun 14, 2021 at 15:54

The first solution will use socket to connect python and Mathematica. We'll set up a server that returns the result of the function that has been defined in Mathematica and on the python side, we'll connect to the server and pass the function's inputs and receive the output.

We'll use SocketListen which was introduced in Mathematica 11.2.

To keep the first solution simple, we'll ignore some standards and use pure python without any third-party module.

### In Mathematica

We use SocketListen to set up a server in port 36000. On each request, it will pass the input data to the function f and return the result. (Note that the standard HTTP response was not used).

ClearAll[f, s];

f[x_, y_] := N[x^2 + y*Sin[x + y] + Sin[5*x]]

s = SocketListen[36000,
Function[{assoc},
WriteString[assoc["SourceSocket"],
ToString[
NumberForm[
Apply[f,
ToExpression[StringReplace[assoc["Data"], "e" -> "E"]]],
DefaultPrintPrecision -> 20]]]]]["Socket"];


### In Python

We'll connect to the port 36000 and define a function that receives an array (uses the first two elements), send them to the socket, and return the result.

import socket
from scipy import optimize

s = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
s.connect(('localhost', 36000))

def send_recieve(x):
s.send(bytes(f"{{{x[0]:.20f},{x[1]:.20f}}}","utf-8"))
temp = s.recv(4096).decode().strip()
if temp != '':
return float(temp)
return None

result = optimize.minimize(send_recieve,[0,1])

s.close()

print(result)


### Result:

Tooks 1 second.

      fun: -0.930090558146007
hess_inv: array([[ 0.03815335, -0.02060073],
[-0.02060073,  0.51689528]])
jac: array([-2.68220901e-07,  6.70552254e-08])
message: 'Optimization terminated successfully.'
nfev: 36
nit: 10
njev: 12
status: 0
success: True
x: array([-0.29631374,  0.14869772])


Notes:

• ToExpression was used instead of ImportString[..., "PythonExpression"]
• DefaultPrintPrecision is set to 20 digits
• If you want to re-evaluate the SocketListen, first it should be closed: Close[s]

## Solution 2 - using wolframclient

First, you should save your function definition to a file in which we could import it in python (default path for windows is Documents).

### In Mathematica

ClearAll[f];

f[x_, y_] := N[x^2 + y*Sin[x + y] + Sin[5*x]]

Save["function.m", f]


### In Python

from wolframclient.evaluation import WolframLanguageSession
from wolframclient.language import wl
from scipy.optimize import minimize

wolfSession = WolframLanguageSession()

def target(x):
return wolfSession.evaluate(wl.Global.f(x[0],x[1]))

res = minimize(target, [0,1])

wolfSession.terminate()

print(res)


### Result:

Tooks 5.6 seconds.

      fun: -0.9300905581460049
hess_inv: array([[ 0.0381547 , -0.02059617],
[-0.02059617,  0.51690923]])
jac: array([-2.68220901e-07,  1.11758709e-07])
message: 'Optimization terminated successfully.'
nfev: 36
nit: 10
njev: 12
status: 0
success: True
x: array([-0.29631374,  0.14869774])

• @ Ben Izd, Thanks a lot for providing different ways to connect Mathematica and python. " socket " Method is very interesting. I request you to Post your own Question and account for this "... ignore some standards and use pure python without any third-party module."(I don't even know what this is) , not that i'm facing a problem with present solution, Felt your method in complete form might be very useful to many us. Jun 14, 2021 at 15:44

This does not answer the question "How can I use Python's SciPy and NumPy functions in Mathematica to find the minimum of a function?". Rather, it should be shown how easily the minimum or maximum of a function can be found with Mathematica 12.3.

f = x^2 + y Sin[x + y] + Sin[5 x];
sol = Solve[{Grad[f, {x, y}] == 0, -10 < x < 10, -10 < y < 11}, {x,y}] // N;
min = MinimalBy[Thread@{f /. sol, sol}, First]
{{-10.3714, {x -> 0.830464, y -> 10.2622}}}


Proof that min is a (global) minimum in the given box:

Eigenvalues@D[f, {{x, y}, 2}] /. min[[1, 2]]
{6.46061, 37.3416}

Show[
Plot3D[f, {x, -5, 5}, {y, -10, 11}, Mesh -> 15,
MeshFunctions -> {#3 &}, PlotTheme -> "FilledSurface"],
Graphics3D[{Red, [email protected], Point[{x, y, f} /. min[[1, 2]]]}]
]


Correction

Solve delivers exact (root) objects. While NMinimize with a given region also calculates the global minimum as above.

reg = Rectangle[{-10, -10}, {10, 11}];
NMinimize[{f, {x, y} \[Element] reg}, {x, y}] // N // AbsoluteTiming
{0.0826161, {-10.3714, {x -> 0.830464, y -> 10.2622}}}


The calculation takes place within 0.08 second. Does this computing time take too long? So, there is no longer any need to reach for other stars.

• Looks like a local minimum to me. Jun 15, 2021 at 15:28
• @rmw, hi thanks for the response. many times we won't be able use Solve and "DifferentialEvolution" in FindMinimum or in NMinimize when function is total numerical values consisting of ?NumericQ. I'm not doubting the numerical capabilities of Mathematica. For my case I got the results linking my Mathematica code with python way faster. Jun 15, 2021 at 16:00
• @ Markus Roellig It is a global minimum in the range 0 < y < 11. If you expand to 0 < y < 11.5, you get a new minimum {-12.2574, {x -> -0.327794, y -> 11.4108}}. This new minimum can already be seen in the beginning in the plot.
– rmw
Jun 15, 2021 at 20:10
• @Gummala Navneeth A lot has changed from version 12.2 to 12.3 with regard to solving transcendent equations, also with numerical numbers (? NumericQ).
– rmw
Jun 15, 2021 at 20:10