MFDS-NM-09: Pendahuluan Optimasi Numerik

Video Penjelasan: https://youtu.be/wXiqOUx8u4c

Open Code in Google Colaboratory

Metode Numerik

</center>

Optimasi Numerik
(C) Taufik Sutanto - 2022
taudata Analytics

https://taudata.blogspot.com/2020/04/mfdsnm-06.html¶

In [1]:

print("Detecting environment: ", end=' ')
try:
    import google.colab
    IN_COLAB = True
    print("Running the code in Google Colab. Installing and downloading dependencies.\nPlease wait...")
    !pip install plotly
except:
    IN_COLAB = False
    print("Running the code locally, please make sure all modules installed in your environment.")
# Please visit https://github.com/taudataid/PINN-DCAI for further detail such as requirements.txt file.

Detecting environment:  Running the code locally, please make sure all modules installed in your environment.

In [2]:

# importing modules
import numpy as np, matplotlib.pyplot as plt, pylab
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from numpy import exp,arange
from pylab import meshgrid,cm,imshow,contour,clabel,colorbar,axis,title,show
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter

np.random.seed(0)
"Modules loaded"

Out[2]:

'Modules loaded'

Outline:¶

• Minimalisasi Fungsi
• Rasio Emas
• Metode Nelder mead
• Minimalisasi dengan Turunan (gradient Descent)

Pendahuluan Aplikasi Turunan

• Statistik (machine learning/data science): Regresi, Klasifikasi, Clustering
• mechanical engineering
• Oceanography, dll

In [3]:

# Fungsi untuk plot
def contourFungsi(X,Y,Z, judul):
    im = imshow(Z,cmap=cm.RdBu, extent=[-3.0,3.0,-3.0,3.0], origin='lower') # drawing the function    
    cset = contour(X,Y,Z,arange(Z.min(),Z.max(),0.1),linewidths=2,cmap=cm.Set2)# adding the Contour lines with labels
    clabel(cset,inline=True,fmt='%1.1f',fontsize=10)
    colorbar(im) # adding the colobar on the right
    title(judul); show()

def plotFs3D(X, Y, Z):
    fig = plt.figure()
    ax = fig.gca(projection='3d')
    surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.RdBu,linewidth=0, antialiased=False)
    ax.zaxis.set_major_locator(LinearLocator(10))
    ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
    ax.set_xlabel('x-axis');ax.set_ylabel('y-axis');ax.set_zlabel('z-axis')
    ax.view_init(elev=25, azim=-120)
    fig.colorbar(surf, shrink=0.5, aspect=5);plt.show()

In [4]:

# Fungsi dan domain
def z_func(x,y):
    return (1-(x**2+y**3))*exp(-(x**2+y**2)/2) 

x = y = arange(-3.0,3.0,0.1)
X,Y = meshgrid(x, y) # grid of point
Z = z_func(X, Y) # evaluation of the function on the grid

In [5]:

%matplotlib qt
# interactive Plot
plotFs3D(X, Y, Z) # Plot Fungsi

C:\Users\taufi\AppData\Local\Temp\ipykernel_18840\1077025168.py:11: MatplotlibDeprecationWarning: Calling gca() with keyword arguments was deprecated in Matplotlib 3.4. Starting two minor releases later, gca() will take no keyword arguments. The gca() function should only be used to get the current axes, or if no axes exist, create new axes with default keyword arguments. To create a new axes with non-default arguments, use plt.axes() or plt.subplot().
  ax = fig.gca(projection='3d')

In [6]:

%matplotlib inline 
# Back to inline plot
contourFungsi(X,Y,Z, 'Contour fungsinya') # Contour Fungsi (Biru Maximum, Merah Minimum)

Meminimumkan Fungsi

In [7]:

%matplotlib inline 
def f(x):
    #return 3*x**2-2*x+5 
    return x**3+x**2-x+1 

X = arange(-2.0, 2.0, 0.1)
#X = arange(0,1.0,0.1)
Y = [f(x) for x in X]
plt.plot(X, Y, label = "Contoh 8.1")
plt.stem([-1, 1/3], [f(-1), f(1/3)], label = "Local Extremum points", linefmt='C2:', markerfmt = 'C2o')
plt.legend(); plt.show()

Diskusi¶

• Kita bisa mencari extremum points, tapi bagaimana kita tau jika titik itu extremum global/tidak?

Golden Ratio (section) Search

Ok ... tapi apanya yang golden?¶

• Source: https://www.youtube.com/watch?v=h0g6_dXO0SM
• Channel: drtedcoe - https://www.youtube.com/channel/UCuZQZbk2Nt7TsZRADcYAyuA

In [8]:

(np.sqrt(5)+1)/2

Out[8]:

1.618033988749895

Nature by Numbers: Golden Ratio, Fibonacci, dan Ciptaan Ilahi¶

• Source: https://www.youtube.com/watch?v=kkGeOWYOFoA
• Channel: Cristóbal Vila - https://www.youtube.com/channel/UCQEgFYXVoiS5RlfsMoAexSA

Algoritma Golden Search

In [9]:

def f(x):
    #return 3*x**2+2*x+5
    #return x+np.cos(x)
    return x**2 - np.sin(x)
X = arange(0.0,1.0,0.01)
Y = [f(x) for x in X]
plt.plot(X, Y, label = "Contoh 8.2")
xPts = [0, 0.450, 0.618, 1]
yPts = [f(x) for x in xPts]
plt.stem(xPts, yPts, label = "Interesting points", linefmt='C2:', markerfmt = 'C2o')
plt.legend(); plt.show()

In [10]:

# Simple version from Wikipedia
def gss(a, b, tol=1e-6):
    gr = (np.sqrt(5) + 1) / 2 # Gold Ratio
    c = b - (b - a) / gr
    d = a + (b - a) / gr 
    output = [(a,b,c,d)]
    while abs(c - d) > tol:        
        if f(c) < f(d):
            b = d
        else:
            a = c
        output.append((a,b,c,d))
        c = b - (b - a) / gr
        d = a + (b - a) / gr
    return output

In [11]:

# Versi Mathews 
def GS(a,b,n): # n = jumlah iterasi    
    r = (np.sqrt(5)-1)/2
    c = a + (1-r)*(b-a)
    d = a + r*(b-a)
    output = [(a,b,c,d)]
    for i in range(n): # or abs(c - d) > tol
        if f(c)<f(d):
            b = d
        else:
            a = c
        c = a + (1-r)*(b-a)
        d = a + r*(b-a)    
        output.append((a,b,c,d))
    return output

In [12]:

A, B = 0.0, 1.0
output = gss(A,B) # versi wikipedia
print('iter, \ta,      \tc,      \td,      \tb')
for i,(a,b,c,d) in enumerate(output):
    print(i, end = ', \t')
    print('{}, \t{}, \t{}, \t{}'.format(a,c,d,b))

iter, 	a,      	c,      	d,      	b
0, 	0.0, 	0.3819660112501052, 	0.6180339887498948, 	1.0
1, 	0.0, 	0.3819660112501052, 	0.6180339887498948, 	0.6180339887498948
2, 	0.2360679774997897, 	0.2360679774997897, 	0.3819660112501051, 	0.6180339887498948
3, 	0.3819660112501051, 	0.3819660112501051, 	0.4721359549995794, 	0.6180339887498948
4, 	0.3819660112501051, 	0.47213595499957933, 	0.5278640450004206, 	0.5278640450004206
5, 	0.3819660112501051, 	0.4376941012509463, 	0.4721359549995794, 	0.4721359549995794
6, 	0.41640786499873816, 	0.41640786499873816, 	0.4376941012509463, 	0.4721359549995794
7, 	0.4376941012509464, 	0.4376941012509464, 	0.4508497187473712, 	0.4721359549995794
8, 	0.4376941012509464, 	0.45084971874737123, 	0.45898033750315453, 	0.45898033750315453
9, 	0.44582472000672974, 	0.44582472000672974, 	0.4508497187473712, 	0.45898033750315453
10, 	0.44582472000672974, 	0.45084971874737123, 	0.45395533876251304, 	0.45395533876251304
11, 	0.4489303400218716, 	0.4489303400218716, 	0.4508497187473712, 	0.45395533876251304
12, 	0.4489303400218716, 	0.45084971874737123, 	0.4520359600370134, 	0.4520359600370134
13, 	0.4489303400218716, 	0.45011658131151383, 	0.4508497187473712, 	0.4508497187473712
14, 	0.449663477457729, 	0.449663477457729, 	0.4501165813115138, 	0.4508497187473712
15, 	0.449663477457729, 	0.4501165813115138, 	0.4503966148935864, 	0.4503966148935864
16, 	0.44994351103980157, 	0.44994351103980157, 	0.45011658131151383, 	0.4503966148935864
17, 	0.45011658131151383, 	0.45011658131151383, 	0.45022354462187414, 	0.4503966148935864
18, 	0.45011658131151383, 	0.4502235446218742, 	0.45028965158322604, 	0.45028965158322604
19, 	0.45011658131151383, 	0.4501826882728657, 	0.45022354462187414, 	0.45022354462187414
20, 	0.45015743766052224, 	0.45015743766052224, 	0.4501826882728657, 	0.45022354462187414
21, 	0.45015743766052224, 	0.45018268827286567, 	0.4501982940095307, 	0.4501982940095307
22, 	0.4501730433971872, 	0.4501730433971872, 	0.4501826882728657, 	0.4501982940095307
23, 	0.4501730433971872, 	0.45018268827286567, 	0.45018864913385226, 	0.45018864913385226
24, 	0.4501790042581738, 	0.4501790042581738, 	0.45018268827286567, 	0.45018864913385226
25, 	0.4501790042581738, 	0.4501826882728657, 	0.45018496511916034, 	0.45018496511916034
26, 	0.4501812811044685, 	0.4501812811044685, 	0.45018268827286567, 	0.45018496511916034

In [13]:

a = 0.0
b = 0.23606797749978967
(a+b)/2, f((a+b)/2)

Out[13]:

(0.11803398874989483, -0.1038280817430638)

In [14]:

output = GS(A,B,24) # versi Matthews
print('iter, \ta,      \tc,      \td,      \tb')
for i,(a,b,c,d) in enumerate(output):
    print(i, end = ', \t')
    print('{}, \t{}, \t{}, \t{}'.format(a,c,d,b))

iter, 	a,      	c,      	d,      	b
0, 	0.0, 	0.3819660112501051, 	0.6180339887498949, 	1.0
1, 	0.0, 	0.2360679774997897, 	0.3819660112501052, 	0.6180339887498949
2, 	0.2360679774997897, 	0.38196601125010515, 	0.47213595499957944, 	0.6180339887498949
3, 	0.38196601125010515, 	0.4721359549995794, 	0.5278640450004206, 	0.6180339887498949
4, 	0.38196601125010515, 	0.4376941012509464, 	0.4721359549995794, 	0.5278640450004206
5, 	0.38196601125010515, 	0.41640786499873816, 	0.4376941012509464, 	0.4721359549995794
6, 	0.41640786499873816, 	0.4376941012509464, 	0.45084971874737123, 	0.4721359549995794
7, 	0.4376941012509464, 	0.45084971874737123, 	0.45898033750315453, 	0.4721359549995794
8, 	0.4376941012509464, 	0.44582472000672974, 	0.4508497187473712, 	0.45898033750315453
9, 	0.44582472000672974, 	0.45084971874737123, 	0.45395533876251304, 	0.45898033750315453
10, 	0.44582472000672974, 	0.4489303400218716, 	0.4508497187473712, 	0.45395533876251304
11, 	0.4489303400218716, 	0.45084971874737123, 	0.4520359600370134, 	0.45395533876251304
12, 	0.4489303400218716, 	0.45011658131151383, 	0.4508497187473712, 	0.4520359600370134
13, 	0.4489303400218716, 	0.449663477457729, 	0.4501165813115138, 	0.4508497187473712
14, 	0.449663477457729, 	0.4501165813115138, 	0.4503966148935864, 	0.4508497187473712
15, 	0.449663477457729, 	0.44994351103980157, 	0.45011658131151383, 	0.4503966148935864
16, 	0.44994351103980157, 	0.45011658131151383, 	0.45022354462187414, 	0.4503966148935864
17, 	0.45011658131151383, 	0.4502235446218742, 	0.45028965158322604, 	0.4503966148935864
18, 	0.45011658131151383, 	0.4501826882728657, 	0.45022354462187414, 	0.45028965158322604
19, 	0.45011658131151383, 	0.45015743766052224, 	0.4501826882728657, 	0.45022354462187414
20, 	0.45015743766052224, 	0.45018268827286567, 	0.4501982940095307, 	0.45022354462187414
21, 	0.45015743766052224, 	0.4501730433971872, 	0.4501826882728657, 	0.4501982940095307
22, 	0.4501730433971872, 	0.45018268827286567, 	0.45018864913385226, 	0.4501982940095307
23, 	0.4501730433971872, 	0.4501790042581738, 	0.45018268827286567, 	0.45018864913385226
24, 	0.4501790042581738, 	0.4501826882728657, 	0.45018496511916034, 	0.45018864913385226

In [15]:

Fc = [f((a+b)/2) for a,b,_,_ in output]
plt.plot(np.linspace(a,b,len(Fc)), Fc, label = "Plot f(c)")
plt.legend(); plt.show()
Fc[-1]

Out[15]:

-0.23246557515815913

Latihan¶

Tentukan f(optimal) dan x optimal setelah 3 iterasi (c+d)/2 yang ke-3

$$f(x) = 7x^2-x+3;\quad [0,2]$$¶

In [16]:

def f(x):
    return 7*x**2-x+3#+np.cos(x)**2

A, B = 0, 2
output = GS(A,B,4)
print('iter, \ta,      \tc,      \td,      \tb')
for i,(a,b,c,d) in enumerate(output):
    print(i, end = ', \t')
    print('{}, \t{}, \t{}, \t{}'.format(a,c,d,b))
    
0.763/2

iter, 	a,      	c,      	d,      	b
0, 	0, 	0.7639320225002102, 	1.2360679774997898, 	2
1, 	0, 	0.4721359549995794, 	0.7639320225002104, 	1.2360679774997898
2, 	0, 	0.2917960675006309, 	0.4721359549995795, 	0.7639320225002104
3, 	0, 	0.1803398874989485, 	0.29179606750063103, 	0.4721359549995795
4, 	0, 	0.11145618000168246, 	0.18033988749894858, 	0.29179606750063103

Out[16]:

0.3815

Extreme Values of f(x,y) ~ Multivariate Function¶

Metode Nelder-Mead

1. Mulai dari sembarang segitiga $(x_k, y_k)$, k=1, 2, 3
2. Hitung $z_k = f(x_k, y_k)$, k=1, 2, 3
3. urutkan $z_k$ (dan urutkan ulang indeks k) sedemikian sehingga $z_i \leq z_2 \leq z_3$
4. Namakan $B = z_1, G = z_2, W = z_3$
5. Hitung (MidPoint of "goodSide") = $M=\frac{B+G}{2}$
6. Hitung titik : $ R = 2M - W $
7. Hitung titik : E = 2R-M (alternatif R)
8. Hitung $C_1 = \frac{W+M}{2}$ dan $C_2 = \frac{M+R}{2}$ (kandidat segitiga BGC)
9. Secara terurut kandidat terbaik pengganti W adalah (E, R, C2, C1)
10. Dengan syarat kandidat tersebut < W -- Lihat Gambar
11. Iterasikan dengan BGW yang baru (setelah diurutkan lagi), hingga |B-G|<toleransi

    

Animasinya¶

• Source: https://www.youtube.com/watch?v=HUqLxHfxWqU
• Channel: brainbrian123 - https://www.youtube.com/channel/UC79jrmp203LSCOxhiJ1wong

Latihan ... = Contoh 8.4¶

$f(x,y)=x^2-4x+y^2-y-xy$

• v1=(0,0), v2=(1.2, 0), v3=(0.0,0.8)

In [17]:

def F(v):
    x, y = v[0], v[1]
    return x**2-4*x+y**2-y-x*y

v1, v2, v3 = (0,0), (1.2,0), (0.0,0.8)
v1 = np.array(v1); v2 = np.array(v2); v3 = np.array(v3)

print(F(v1), F(v2), F(v3))
B, G, W = v3, v1, v2
M = (B+G)/2
R = 2*M-W
E = 2*R-M
print(F(E),F(B),F(G))

B, G, W = E, B, v2
M = (B+G)/2
R = 2*M-W
E = 2*R-M

print(F(E),F(B),F(G))
B

0 -3.36 -0.15999999999999992
18.48 -0.15999999999999992 0
83.99999999999999 18.48 -0.15999999999999992

Out[17]:

array([-2.4,  1.2])

Minimization using Derivative¶

• Konsepnya adalah Newton (Secant) pada permasalahan mencari akar f'(x) = 0
• Bracketing di Mathews (hal 410) ... HHhmmm... do you see why we're not discussing it further?
• .
• $\min -2x^4 + 3x-7$ jika $X_0=1$ hitung $x_3^*$ menggunakan metode Newton.

In [18]:

def F(x):
    return -2*x**4 + 3*x - 7
def g(x):
    return -8*x**3 + 3
def gg(x):
    return -24*x**2

# Newton (X3, xo=1)
xo =1
x1 = xo - g(xo)/gg(xo)
x2 = x1 - g(x1)/gg(x1)
x3 = x2 - g(x2)/gg(x2)
x3

Out[18]:

0.7211757921789095

(Stochastic-Mini batch) Gradient Descent with momentum¶

$$ w_{n+1} = w_n - \eta \nabla f(w_n) - \alpha \Delta w_n$$¶

• $\eta$ adalah learning rate, bayangkan sebagai "langkah" menuju nilai optimal fungsi $f(w_n)$
• $\alpha$ disebut sebagai momentum, ia menstabilkan learning rate dan "smoothing" osilasi bergeraknya $w_n$. Silahkan simak penjelasan Andrew Ng disini: https://www.youtube.com/watch?v=k8fTYJPd3_I
• Referensi lain tentang GD dan variasinya: https://www.analyticsvidhya.com/blog/2021/06/guide-to-gradient-descent-and-its-variants-with-python-implementation/
• Pada dasarnya Stochastic atau mini batch GD hanya bergantung penggunaan sampel data untuk update weight saja.
• Mari kita simulasikan pada fungsi sederhana berikut:

$$ f(x_1, x_2) = 1.5 x_1^2.0 + 2 x_1 x_2 + 3 x_2^2 - 2 x_1 + 8 x_2 $$¶

• Nilai optimal (eksak): $(x_1, x_2) = (2, -2)$, verify!

In [19]:

def f(x1, x2):
    return 1.5*x1**2.0+2.0*x1*x2+3.0*x2**2-2.0*x1+8*x2

x1 = np.linspace(-10, 10, 100)
x2 = np.linspace(-10, 10, 100)
X, Y = np.meshgrid(x1, x2) #x = np.column_stack((X, Y))
Z = f(X, Y)

fig = make_subplots(rows=1, cols=1, specs=[[{'type': 'surface'}]])
fig.add_trace(go.Surface(x=X, y=Y, z=Z, colorscale='Viridis', showscale=False), row=1, col=1)
fig.update_layout(title_text='Surface Plot Fungsi', height=500, width=500)
fig.update_layout(scene=dict(xaxis_title='X', yaxis_title='Y', zaxis_title='$Z=f(x_1, x_2)$'))
fig.show()

In [20]:

def df(x1, x2, h=10**-4): # Misal Forward difference
    return np.array([(f(x1+h, x2)-f(x1, x2))/h, (f(x1, x2+h)-f(x1, x2))/h])

x0 = np.array([1, -.2])
x1 = x0.copy() # X1 initialization for the momentum
lr, alfa_ = 0.2, 0.1 # Learning Rate & momentum
tol_ = 10**-4 # Tolerance
hist_f = [] # History of f(x) values from the iterations
iter_, maxIter_ = 0, 10
while iter_<=maxIter_ or abs(y0-y1)>tol_:
    y0 = f(x0[0], x0[1])
    hist_f.append(y0)
    x1 = x0 - lr*df(x0[0], x0[1]) - alfa_*(x1-x0)
    y1 = f(x1[0], x1[1])
    x0 = x1.copy()
    iter_ += 1
        
fig = plt.figure()
pylab.plot(range(len(hist_f)), hist_f,'blue', label='History function cost') 
plt.xlabel('Number of iteration')
plt.ylabel('Value of cost Function')
plt.legend(loc="upper left")
plt.show()

print("Optimal X = ", x1, " after ", iter_ , " iterations")
"min: ", f(x1[0], x1[1]), "surrounding value: ", f(x1[0]+0.05, x1[1]+0.05)

Optimal X =  [ 1.99444224 -1.99732922]  after  11  iterations

Out[20]:

('min: ', -9.999961954932456, 'surrounding value: ', -9.984033082540044)

End of Module¶

---

Referensi:

1. Johansson, R. (2018). Numerical Python: Scientific Computing and Data Science Applications with Numpy, SciPy and Matplotlib. Apress.
2. John H. Mathews, Numerical Methods for Mathematics, Prentice Hall, 1992 [Refernsi Utama]
3. Heath, M. T. (2018). Scientific computing: an introductory survey (Vol. 80). SIAM.
4. Conte, S. D., u0026amp; De Boor, C. (2017). Elementary numerical analysis: an algorithmic approach (Vol. 78). SIAM.
5. https://www.cs.princeton.edu/courses/archive/fall09/cos323/lectures/cos323_s06_lecture03_optimization.ppt
6. https://eclass.aueb.gr/modules/document/file.php/STAT263/W.%20John%20Braun%2C%20Duncan%20J.%20Murdoch%20-%20A%20First%20Course%20in%20Statistical%20Programming%20with%20R%20%282016%2C%20Cambridge%29.pdf
7. http://www.mee.tcd.ie/~sigmedia/pmwiki/uploads/Main.Tutorials/Introduction_to_Numerical_Optimization.pdf
8. Materi Bahasa Indonesia

Referensi Video

• https://www.youtube.com/watch?v=ms3aKKW_iRc
• https://www.youtube.com/watch?v=hTslJ9A3peA
• https://www.youtube.com/watch?v=tnOOiSjTI8I&list=PL6EA0722B99332589&index=4