```python
import numpy as np

import numpy as np
import scipy as sp
from scipy import linalg
from scipy import optimize
from scipy import interpolate
import sympy as sm

%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('seaborn-whitegrid')
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
```

# Gradient descent

Let $\boldsymbol{x} = \left[\begin{array}{c}
x_1 \\
x_2\\
\end{array}\right]$ be a two-dimensional vector. Consider the following algorithm:

**Algorithm:** `gradient_descent()`

**Goal:** Minimize the function $f(\boldsymbol{x})$.

1. Choose a tolerance $\epsilon>0$, a scale factor $ \Theta > 0$, and a small number $\Delta > 0$
2. Guess on $\boldsymbol{x}_0$ and set $n=1$
3. Compute a numerical approximation of the jacobian for $f$ by

    $$
    \nabla f(\boldsymbol{x}_{n-1}) \approx \frac{1}{\Delta}\left[\begin{array}{c}
    f\left(\boldsymbol{x}_{n-1}+\left[\begin{array}{c}
    \Delta\\
    0
    \end{array}\right]\right)-f(\boldsymbol{x}_{n-1})\\
    f\left(\boldsymbol{x}_{n-1}+\left[\begin{array}{c}
    0\\
    \Delta
    \end{array}\right]\right)-f(\boldsymbol{x}_{n-1})
    \end{array}\right]
    $$

4. Stop if the maximum element in $|\nabla f(\boldsymbol{x}_{n-1})|$ is less than $\epsilon$
5. Set $\theta = \Theta$ 
6. Compute $f^{\theta}_{n} = f(\boldsymbol{x}_{n-1} - \theta \nabla f(\boldsymbol{x}_{n-1}))$
7. If $f^{\theta}_{n} < f(\boldsymbol{x}_{n-1})$ continue to step 9
8. Set $\theta = \frac{\theta}{2}$ and return to step 6     
9. Set $x_{n} = x_{n-1} - \theta \nabla f(\boldsymbol{x}_{n-1})$
10. Set $n = n + 1$ and return to step 3

**Question:** Implement the algorithm above such that the code below can run.

Jacobian function


```python
def _rosen(x1,x2):
    f = (1.0-x1)**2 + 2*(x2-x1**2)**2

x1 = sm.symbols('x_1')
x2 = sm.symbols('x_2')
f = (1.0-x1)**2 + 2*(x2-x1**2)**2

Df = sm.Matrix([sm.diff(f,i) for i in [x1,x2]])
Df
```




    Matrix([
    [-8*x_1*(-x_1**2 + x_2) + 2*x_1 - 2.0],
    [                   -4*x_1**2 + 4*x_2]])




```python

def rosen(x):
    return _rosen(x[0],x[1])
def rosen_jac(x):
    return np.array([-(2.0-2*x[0])-8*x[0]*(x[1]-x[0]**2),4*(x[1]-x[0]**2)])

def gradient_descent(f,x0,epsilon=1e-6,Theta=0.1,Delta=1e-8,max_iter=10_000):

    """ minimize function with gradient descent
        
    Args:

        f (callable): function
        x0 (np.ndarray): initial values
        jac (callable): jacobian - MADE MYSELF IN FUNCTION
        alpha (list): potential step sizes - CHANGED TO Theta
        max_iter (int): maximum number of iterations
        tol (float): tolerance
        
    Returns:
    
        x (np.ndarray): minimum
        n (int): number of iterations used
        
    """
    
    # step 1: initialize
    x = x0
    fx = f(x0)
    n = 1
    
    # step 2-6: iteration
    while n < max_iter:
            
        x_prev = x
        fx_prev = fx
        
        # step 2: evaluate gradient
        jacx = rosen_jac(x)
        
        # step 3: find good step size (line search)
        fx_ast = np.inf
        theta_ast = Theta
        theta = Theta / 2
        x = x_prev - theta*jacx
        fx = f(x)
        if fx < fx_ast:
            fx_ast = fx
            theta_ast = theta
        
        # step 4: update guess
        x = x_prev - theta_ast*jacx
                            
        # step 5: check convergence
        fx = f(x)
        if abs(fx-fx_prev) < epsilon:
            break
            
        # d. update i
        n += 1
        
    return x,n
    pass
```

**Test case:**


```python
def rosen(x):
    return (1.0-x[0])**2+2*(x[1]-x[0]**2)**2

x0 = np.array([1.1,1.1])
try:
    x,it = gradient_descent(rosen,x0)
    print(f'minimum found at ({x[0]:.4f},{x[1]:.4f}) after {it} iterations')
    assert np.allclose(x,[1,1])
except:
    print('not implemented yet')
```

    minimum found at (1.0049,1.0109) after 131 iterations
    not implemented yet



```python
gradient_descent(rosen,x0)
```




    (array([1.0049464 , 1.01092185]), 131)




```python
def minimize_gradient_descent(f,x0,jac,alphas=[0.01,0.05,0.1,0.25,0.5,1],max_iter=500,tol=1e-8):
    """ minimize function with gradient descent
        
    Args:

        f (callable): function
        x0 (np.ndarray): initial values
        jac (callable): jacobian
        alpha (list): potential step sizes
        max_iter (int): maximum number of iterations
        tol (float): tolerance
        
    Returns:
    
        x (np.ndarray): minimum
        n (int): number of iterations used
        
    """
    
    # step 1: initialize
    x = x0
    fx = f(x0)
    n = 1
    
    # step 2-6: iteration
    while n < max_iter:
            
        x_prev = x
        fx_prev = fx
        
        # step 2: evaluate gradient
        jacx = jac(x)
        
        # step 3: find good step size (line search)
        fx_ast = np.inf
        alpha_ast = np.nan
        for alpha in alphas:
            x = x_prev - alpha*jacx
            fx = f(x)
            if fx < fx_ast:
                fx_ast = fx
                alpha_ast = alpha
        
        # step 4: update guess
        x = x_prev - alpha_ast*jacx
                            
        # step 5: check convergence
        fx = f(x)
        if abs(fx-fx_prev) < tol:
            break
            
        # d. update i
        n += 1
        
    return x,n
```


```python
minimize_gradient_descent(rosen,x0,rosen_jac)
```




    (array([1.00020707, 1.00046163]), 51)



Can use Nelder-Mead without analytical hessian

Newton if we have the analytical hessian

BFGS is the best without analytical

Maybe potential step sizes are the theta function with theta/2

##

459856d3a2a2da20784797fa3a1408076480a448
ipynb
Jupyter Notebook
Exam project/Marcus3.ipynb
notnasobe666/BlackHatGang
e6ec842cfaedaa3752003d7dc1d50ef66091fe6d
MIT
Exam project/Marcus3.ipynb
notnasobe666/BlackHatGang
e6ec842cfaedaa3752003d7dc1d50ef66091fe6d
MIT
Exam project/Marcus3.ipynb
notnasobe666/BlackHatGang
e6ec842cfaedaa3752003d7dc1d50ef66091fe6d
MIT
Qwen/Qwen-72B
1. YES
2. YES