GitHub Repository: restrepo/ComputationalMethods
Path: blob/master/material/Minimization.ipynb
⁹³⁴ views

Kernel: Python 3 (ipykernel)

Minimization

Find the minumum of a function

Bibliography

Example

In [1]:

%pylab inline

Out[1]:

Populating the interactive namespace from numpy and matplotlib

In [2]:

import numpy as np
import scipy.optimize as optimize
import scipy.interpolate as interpolate
import pandas as pd

Example

Consider the following dataset

In [3]:

df=pd.DataFrame({'X':[2.5,3.1,4.5,5,5.9,6.2],
              'Y':[3,-1.01,1.5,0.7,2.8,1.5]})
df

Out[3]:

Laggrange interpolation

In [10]:

import numpy as np

In [15]:

p=np.poly1d([1,2,3])

In [16]:

print(p)

Out[16]:

   2
1 x + 2 x + 3

In [17]:

p(1)

Out[17]:

6

In [18]:

p([1,2,3])

Out[18]:

array([ 6, 11, 18])

In [19]:

p.roots

Out[19]:

array([-1.+1.41421356j, -1.-1.41421356j])

In [21]:

print(p.deriv())

Out[21]:

 
2 x + 2

In [4]:

pol=interpolate.lagrange(df.X,df.Y)

In [5]:

print(pol)

Out[5]:

         5         4         3         2
-0.8941 x + 19.97 x - 174.7 x + 746.7 x - 1557 x + 1264

df.X is a Pandas Series

In [22]:

df.X

Out[22]:

  2.5
  3.1
  4.5
  5.0
  5.9
  6.2
Name: X, dtype: float64

To obtain some specific value by using slices, we must use .iloc

Note that df.X[3]=df.X.loc[3]=df.X.iloc[3]

In [23]:

df.X[3],df.X.loc[3],df.X.iloc[3]

Out[23]:

(5.0, 5.0, 5.0)

In [24]:

df.X.iloc[-1],df.X.max()

Out[24]:

(6.2, 6.2)

works!

In [25]:

df

Out[25]:

In [26]:

x=np.linspace(df.X.iloc[0],df.X.iloc[-1],100)
plt.plot(x,pol(x))
plt.plot(df.X,df.Y,'ro')
plt.xlabel('X')
plt.ylabel('Y')
plt.grid()

Out[26]:

Hermite interpolation

The recommend degree for the Hermite polynomial is $n-1$ where $n$ is the number of data of the dataset

In [27]:

H=np.polynomial.hermite.Hermite.fit(df.X,df.Y,5)

In [29]:

print(H)

Out[29]:

5768090026661477 - 21.34734250551783·H₁(x) + 3.149481575369614·H₂(x) -
275256998110663·H₃(x) + 0.3812886076702684·H₄(x) -
6054571273241235·H₅(x)

In [56]:

plt.plot(x,H(x),'k-')
plt.plot(x,pol(x),'c--')
plt.plot(df.X,df.Y,'ro')
plt.grid()
plt.xlabel('$x$')
plt.ylabel('$H(x)$')

Out[56]:

Text(0,0.5,'$H(x)$')

Finding the local minimum of a function

Finding the first minimum close to 3 (which corresponds to the global minimum), and the second close to 5 (a local minimum)

In [31]:

help(optimize.fmin_powell)

In [32]:

min1=optimize.fmin_powell(pol,3,full_output=True)
min2=optimize.fmin_powell(pol,5.2,full_output=True)

Out[32]:

Optimization terminated successfully.
         Current function value: -1.374113
         Iterations: 2
         Function evaluations: 28
Optimization terminated successfully.
         Current function value: 0.699898
         Iterations: 2
         Function evaluations: 50

In [33]:

min1

Out[33]:

(array([2.92784955]), array(-1.37411316), array([[1.]]), 2, 28, 0)

In [34]:

min2

Out[34]:

(array([5.00492929]), array(0.69989815), array([[0.00045871]]), 2, 50, 0)

In [35]:

print('The global minimum is f(x)={} for x={};\n the local minimum is f(x)={} for x={}'.format(
    min1[1],min1[0],min2[1],min2[0]))

Out[35]:

The global minimum is f(x)=-1.3741131581291484 for x=[2.92784955];
 the local minimum is f(x)=0.6998981514389016 for x=[5.00492929]

Activity Find the maximum values of the Hermite interpolation function of degree 5 to the set of points: https://github.com/restrepo/ComputationalMethods/blob/master/data/hermite.csv

In [13]:

df=pd.read_csv('https://raw.githubusercontent.com/restrepo/ComputationalMethods/master/data/hermite.csv')
H=np.polynomial.hermite.Hermite.fit(df.X,df.Y,5)
x=np.linspace(df.X.iloc[0],df.X.iloc[-1],100)
plt.plot(x,-H(x),'k-')
plt.grid()
plt.xlabel('$x$')
plt.ylabel('$H(x)$')

Out[13]:

Text(0, 0.5, '$H(x)$')

fmin_powell try to search the global minimum

In [14]:

optimize.fmin_powell(-H,4.5)

Out[14]:

Optimization terminated successfully.
         Current function value: -2.803764
         Iterations: 2
         Function evaluations: 27

array([5.918213])

Find a local minumum

close minimum

In [15]:

min1=optimize.minimize(-H,x0=4.5)
min1['x']

Out[15]:

array([4.01502454])

In [16]:

min1=optimize.minimize(-H,5.2)
min1['x']

Out[16]:

array([5.91821428])

minimum in a range

In [17]:

min1=optimize.minimize(-H,x0=4,bounds=((3.5, 4.5), ) )
min1['x']

Out[17]:

array([4.01502381])

In [18]:

xmin_local = optimize.fminbound(-H, 3.5, 4.5)
xmin_local

Out[18]:

4.01502533002799

In [19]:

optimize.brute(-H,ranges=((3.5,4.5,0.1),))

Out[19]:

array([4.01503906])

Find a global minumum (alternative)

In [20]:

min1=optimize.basinhopping(-H,5.5)
min1['x']

Out[20]:

array([5.91821428])

In [24]:

%pylab inline

Out[24]:

Populating the interactive namespace from numpy and matplotlib

In [21]:

import pandas as pd

The Higgs potential

To write greek letter inside a cell use the $\rm \LaTeX$ macro and the tab, e.g: \mu<TAB>, to produce μ

In [1]:

%pylab inline

Out[1]:

Populating the interactive namespace from numpy and matplotlib

In [2]:

import numpy as np

In [3]:

α=2

In [4]:

m_proton = 1 #GeV/c²

In [5]:

m_H=126 # GeV/c^2 (c→1)
G_F=1.1663787E-5 #GeV^-2
v=1/np.sqrt(np.sqrt(2.)*G_F) # GeV
μ=np.sqrt(m_H**2/2) #GeV
λ=m_H**2/(2.*v**2)
μ,λ,v

Out[5]:

(89.09545442950498, 0.13093799079487806, 246.21965079413738)

In [6]:

ϕ=np.linspace(-300,300)
Vp=lambda ϕ: 0.5*μ**2*ϕ**2+0.25*λ*ϕ**4
plt.plot(ϕ, Vp(ϕ) )
plt.xlabel(r'$\phi$ [GeV]',size=20 )
plt.ylabel(r'$V(\phi)$ [GeV]',size=20)
plt.xlabel(r'$\phi$ [GeV]',size=20 )
plt.ylabel(r'$V(\phi)$ [GeV$^4$]',size=20)
plt.grid()

Out[6]:

In [7]:

μ=μ*1j
V=lambda ϕ: Vp(ϕ).real

In [8]:

Vp(ϕ)[0].real

Out[8]:

-92060568.64037189

In [9]:

ϕ=np.linspace(-380,380)
plt.plot(ϕ, V(ϕ) )
plt.xlabel(r'$\phi$ [GeV]',size=20 )
plt.ylabel(r'$V(\phi)$ [GeV]',size=20)
plt.xlabel(r'$\phi$ [GeV]',size=20 )
plt.ylabel(r'$V(\phi)$ [GeV$^4$]',size=20)
plt.grid()

Out[9]:

In [11]:

from scipy import optimize

In [13]:

fp=optimize.fmin_powell(V,200,ftol=1E-16,full_output=True)

Out[13]:

Optimization terminated successfully.
         Current function value: -120308559.069597
         Iterations: 4
         Function evaluations: 74

In [14]:

ϕ_min=fp[0]

In [15]:

print(ϕ_min,v)

Out[15]:

[246.21964988] 246.21965079413738

Minimization in higher dimensions

For a complex scalar field with potential $\begin{equation} V(\phi)=\mu^2\phi^*\phi + \lambda (\phi^*\phi)^2 \end{equation}$ with $\begin{equation} \phi=\frac{\phi_1+i\phi_2 }{\sqrt{2} } \end{equation}$ and $\mu^2<0$ , and $\lambda>0$ , find some of the infinite number of minimum values of $\phi$ , as illustrated in the figure, with the plane, $\phi_1-\phi_2$ , moved to the minimum to easy the visualization. Expanding in terms of the real and imaginary part of $\phi$ $\begin{equation} V(\phi)=\frac{\mu^2}{2}\left(\phi_1^2+\phi_2^2 \right) + \frac{\lambda}{4}\left( \phi_1^2+\phi_2^2\right)^2 \end{equation}$

In [16]:

%pylab inline
from scipy import optimize
import pandas as pd

Out[16]:

Populating the interactive namespace from numpy and matplotlib

In [17]:

def f(ϕ,m_H=126,G_F=1.1663787E-5):
    v=1/np.sqrt(np.sqrt(2.)*G_F) # GeV
    μ=np.sqrt(m_H**2/2)*1j #imaginary mass
    λ=m_H**2/(2.*v**2)
    #print(μ,λ)
    return ( 0.5*μ**2*(ϕ[0]**2+ϕ[1]**2)+0.25*λ*(ϕ[0]**2+ϕ[1]**2)**2 ).real

Check a point of the function

In [18]:

f([0,10])

Out[18]:

-396572.6550230127

Check the minimim obtained when an inizialization point at $\phi_0=(0,0)$

In [19]:

fmin=optimize.fmin_powell(f,x0=[0,10],ftol=1E-16,full_output=True)
fmin[0]

Out[19]:

Optimization terminated successfully.
         Current function value: -120308559.069597
         Iterations: 3
         Function evaluations: 111

array([246.21914011,  -0.50137284])

Check the proyection of the minimum in the plane $\phi_1-\phi_2$

In [20]:

print('V(ϕ)={} GeV^4'.format(fmin[1].round(1)))

Out[20]:

V(ϕ)=-120308559.1 GeV^4

In [21]:

np.sqrt( fmin[0][0]**2+fmin[0][1]**2 )

Out[21]:

246.21965057683713

For random initialization points, we can get several minima

In [22]:

np.random.uniform(-300,300,2)

Out[22]:

array([255.08701519, 178.66169419])

In [23]:

df=pd.DataFrame()
for i in range(1000):
    ϕ0=np.random.uniform(-300,300,2)
    ϕmin=optimize.fmin_powell(f,x0=ϕ0,ftol=1E-16,disp=False)
    df=df.append({'ϕ1':ϕmin[0],'ϕ2':ϕmin[1]},ignore_index=True)

Projection of the minima in the plane, $\phi_1,\phi_2$

In [24]:

df[:3]

Out[24]:

In [25]:

plt.figure( figsize=(6,6) )
plt.plot(df['ϕ1'],df['ϕ2'],'r.',label=r'$\phi_{{\rm min}}={}$ GeV'.format(
    np.sqrt(df.loc[0,'ϕ1']**2+df.loc[0,'ϕ2']**2).round(1)))
plt.plot(df.loc[0,'ϕ1'],df.loc[0,'ϕ2'],'k*',label='Our Universe',markersize=10)
plt.xlabel('$\phi_1$',size=15)
plt.ylabel('$\phi_2$',size=15)
plt.legend(loc='best')
plt.grid()