<jupyter_start><jupyter_text>Optimization Algorithms<NAME><jupyter_code>import numpy as np
from numpy import linalg as LA
import matplotlib.pyplot as plt
import torch as th
from torch.autograd.functional import hessian<jupyter_output><empty_output><jupyter_text>F function below is multi modal function, non-convex because there are multiple optimal solutions. (Convex functions have just one optimal solution)<jupyter_code>def f(x, y):
    return 2 * np.sum(np.log(LA.norm(x-y, axis=1)), axis=0)

y = np.array([[1., 0.], [-1., 1.], [-1.5, -1.]])
f([0,0], y)
#y = y.reshape(1,3,2)
x_ = np.linspace(-2, 2, 300)
y_ = np.linspace(-2, 2, 300)
xx, yy = np.meshgrid(x_, y_)
xy = np.array([xx.flatten(),yy.flatten()])
zz = f(xy, y.reshape(3,2,1))
_ = plt.contour(xx, yy, zz.reshape(xx.shape), levels=30)
def grad_f(x, y):
    return np.sum(2 * (x-y) / (LA.norm(x - y, axis=1).reshape(3,1))**2, axis=0)

x = np.array([0.5, 1.5])
y = np.array([[1., 0.], [-1., 1.], [-1.5, -1.]])

grad_f(x, y)<jupyter_output><empty_output><jupyter_text>Gradient Descent  Guess a solution $\vec x_0$ determine         $$ \vec \triangle_0 = - \nabla f(\vec x_0)  $$update estimate        $$ \vec x_1 = \vec x_0 + \eta \vec \triangle_0$$repeat        $$ \vec \triangle_t = - \nabla f(\vec x_t)  $$        $$ \vec x_{t+1} = \vec x_t + \eta \vec \triangle_t$$<jupyter_code>def grad_descent(x, y, lr=0.01):
    x_hist = np.array(x)

    for i in range(35):
        grad_x = grad_f(x, y)
        x = x - lr * grad_x
        x_hist = np.vstack((x_hist, x))
        
    return x_hist
        
lr = 0.05
x = np.array([1, 1.5])
y = np.array([[1., 0.], [-1., 1.], [-1.5, -1.]])
x_hist = grad_descent(x, y, lr=lr)
#print(x_hist)
plt.contour(xx, yy, zz.reshape(xx.shape), levels=30)
_ = plt.plot(x_hist[:,0],x_hist[:,1], "->")<jupyter_output><empty_output><jupyter_text>Gradient Descent with Line Searches Guess a solution $\vec x_0$ repeat        $$ \vec \triangle_t = - \nabla f(\vec x_t)  $$        $$ \eta_t = argmin_{\eta} \space f(\vec x_t + \eta \vec \triangle_t)    $$        $$ \vec x_{t+1} = \vec x_t + \eta_t \vec \triangle_t$$<jupyter_code>def grad_desc_line_search(x, y):
    
    x_hist = np.array(x)
    
    for i in range(100):
          
        change = -1 * grad_f(x, y)
        
        lr_list = np.linspace(0.,1,1000)
        min_value = f(x + lr_list[0] * change, y)
        min_lr = lr_list[0]

        for lr in lr_list:
            value = f(x + lr * change, y)
            if value < min_value:
                min_value = value
                min_lr = lr
        
        if min_lr!=0:
            x = x + min_lr * change
            x_hist = np.vstack((x_hist, x))
    
    return x_hist

x = np.array([1, 1.5])
y = np.array([[1., 0.], [-1., 1.], [-1.5, -1.]])

x_hist2 = grad_desc_line_search(x, y)
print(x_hist2[-1])
print(x_hist2.shape)
plt.contour(xx, yy, zz.reshape(xx.shape), levels=30)
_ = plt.plot(x_hist2[:,0], x_hist2[:,1], "-*r")<jupyter_output><empty_output><jupyter_text>Gradient Descent Using Conjugate Gradients Guess a solution $\vec x_0$         $$ \vec \triangle_0 = - \nabla f(\vec x_0)  $$        $$ \vec c_0 = \vec \triangle_0  $$        $$ \eta_0 = argmin_{\eta} \space f(\vec x_0 + \eta \vec c_0)  $$        $$ \vec x_1 = \vec x_0 + \eta_0 \vec c_0$$repeat        $$ \vec \triangle_t = - \nabla f(\vec x_t)  $$        $$ \vec c_t = \vec \triangle_t + \beta_t \vec c_{t-1}  $$        $$ \eta_t = argmin_{\eta} \space f(\vec x_t + \eta \vec c_t)  $$        $$ \vec x_{t+1} = \vec x_t + \eta_t \vec c_t$$        For beta:the Fletcher-Reeves formula:        $$ \beta_t = max \bigg \{ 0, \frac{\vec \triangle_t^T \vec \triangle_t}{\vec \triangle_{t-1}^T \vec \triangle_{t-1}} \bigg \}$$        the Polak-Ribière formula:              $$ \beta_t = max \bigg \{ 0, \frac{\vec \triangle_t^T (\vec \triangle_t - \vec \triangle_{t-1})}{\vec \triangle_{t-1}^T \vec \triangle_{t-1}} \bigg \}$$<jupyter_code>def f(x, y):
    return 2 * np.sum(np.log(LA.norm(x-y, axis=1)), axis=0)

def grad_f(x, y):
    return np.sum(2 * (x-y) / (LA.norm(x - y, axis=1).reshape(3,1))**2, axis=0)

def min_lr_search(x, y, change):
    
    lr_list = np.linspace(0.,1,1000)
    min_value = f(x + lr_list[0] * change, y)
    min_lr = lr_list[0]

    for lr in lr_list:
        value = f(x + lr * change, y)
        if value < min_value:
            min_value = value
            min_lr = lr
            
    return min_lr

"""
For beta:

Fletcher-Reeves formula
beta = np.max((0, (change[i].T @ change[i])/(change[i-1].T @ change[i-1])))

Polak-Ribière formula     
beta = np.max((0, (change[i].T @ (change[i] - change[i-1])) / (change[i-1].T @ change[i-1]) ))
"""
def grad_desc_conj_grad(x, y):
    
    c = []
    change = []
    x_hist = np.array(x)
    
    change.append(-1 * grad_f(x, y)) 
    c.append(change[0])
    min_lr = min_lr_search(x, y, c[0])

    x = x + min_lr * c[0]
    x_hist = np.vstack((x_hist, x))
    
    
    for i in range(1,30):
        
        change.append(-1 * grad_f(x, y))
        
        beta = np.max((0, (change[i].T @ (change[i] - change[i-1])) / (change[i-1].T @ change[i-1]) ))  #PR
        #beta = np.max((0, (change[i].T @ change[i])/(change[i-1].T @ change[i-1])))  #FR
          
        c.append(change[i] + beta * c[i-1])

        min_lr = min_lr_search(x, y, c[i])
        
        if min_lr != 0: 
            x = x + min_lr * c[i]
            x_hist = np.vstack((x_hist, x))
    
    return x_hist

x = np.array([1, 1.5])
y = np.array([[1., 0.], [-1., 1.], [-1.5, -1.]])

x_hist2 = grad_desc_conj_grad(x, y)
print("x_hist2.shape= ", x_hist2.shape)
print(x_hist2[-1])
plt.contour(xx, yy, zz.reshape(xx.shape), levels=30)
_ = plt.plot(x_hist2[:,0], x_hist2[:,1], "-*r")<jupyter_output><empty_output><jupyter_text>Gradient Descent Using Newton's Method Guess a solution $\vec x_0$ repeat        $$ \vec \triangle_t = - \nabla f(\vec x_t)  $$        $$ \vec x_{t+1} = \vec x_t + \eta H^{-1} \vec \triangle_t $$<jupyter_code>def f_th(x): 
    y = th.tensor([[1., 0.], [-1., 1.], [-1.5, -1.]])
    return 2 * th.sum(th.log(th.norm(x-y, dim=1)), dim=0)

def grad_f(x):
    y = np.array([[1., 0.], [-1., 1.], [-1.5, -1.]])
    return np.sum(2 * (x-y) / (LA.norm(x - y, axis=1).reshape(3,1))**2, axis=0)

def newtons_method(x, lr):

    x_hist = np.array(x)

    for i in range(5):
        change = -1 * grad_f(x)
        hessian = th.autograd.functional.hessian(f_th, th.tensor(x))

        x = x - lr * LA.inv(hessian.numpy()) @ change
        x_hist = np.vstack((x_hist, x))

    return x_hist    

lr = 0.3
x = np.array([1.,1.])

x_hist = newtons_method(x, lr)
x_hist[-1]
plt.contour(xx, yy, zz.reshape(xx.shape), levels=30)
_ = plt.plot(x_hist[:,0], x_hist[:,1], "-*r")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.20166443750928817, 'fraction_numerical': 0.034329023629068214, 'mean_word_length': 2.467284904688305, 'pattern_counts': {'":': 0, '<': 28, '<?xml version=': 0, '>': 27, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '278935', 'license': 'MIT', 'n_tokens_mistral': 2989, 'n_tokens_neox': 2893, 'n_words': 619, 'path': 'Numerical_Optimization.ipynb', 'repo_name': 'elifcansuyildiz/MachineLearningNotebooks', 'size': '6740'}