```python
from plotly.offline import init_notebook_mode, iplot
from IPython.display import display, HTML
import numpy as np
```


```python
init_notebook_mode(connected=True)
```

# Gradient Descent

[Wiki: gradient descent](https://en.wikipedia.org/wiki/Gradient_descent)

In general, our optimization problem can be viewd as following: we are minimizing the **Objective** function with respect to $\theta$:

$$\text{argmin}_{\theta}\ L(\theta | x) \quad \textbf{Objective}\ \text{function},$$

for simplicity, we only consider the class of convex function. (Only one minimal point)

**Goal**

Generate a sequence of $\theta_i$ such that $L(\theta_1|x) \ge L(\theta_2|x) \ge \dots \ge L(\theta_i|x) \ge L(\theta_{i+1}|x) \ge \dots$

**Steps**

1. Calculate the gradient of objective function$\nabla L(\theta_i)$
2. Move to the negative direction of derivative $-\nabla L(\theta_i)$
3. Update $\theta_{i + 1} = \theta_i -\gamma\nabla L(\theta_i | x), \ \gamma > 0$
4. Check whether converge. Does $\nabla L(\theta_i)$ close to zero? Since if it is a minimal point, $\nabla L(\theta) = 0$



# Examples

The gradient descent algorithm is applied to find a local minimum of the function $f(x)=x^{4}-3x^{3}+2$, with derivative $f'{(x)}=4x^{3}-9x^{2}$. Here is an implementation is:


```python
# From calculation, it is expected that the local minimum occurs at x=9/4

cur_x = 3 # The algorithm starts at x=6
gamma = 0.01 # step size multiplier
precision = 0.00001
previous_step_size = 1 
max_iters = 10000 # maximum number of iterations
iters = 0 #iteration counter

ff = lambda x: x**4 - 3*x**3 + 2
df = lambda x: 4 * x**3 - 9 * x**2

xpath = []
ypath = []

while previous_step_size > precision and iters < max_iters:
    xpath.append(cur_x)
    ypath.append(ff(cur_x))
    prev_x = cur_x
    cur_x -= gamma * df(prev_x)
    previous_step_size = abs(cur_x - prev_x)
    iters+=1

print("The local minimum occurs at", cur_x)
```


```python
def plot_path(x, y):
    
    xm, xM = min(x) - 0.5, max(x) + 1
    ym, yM = min(y) - 0.5, max(y) + 1

    data = [dict(x = x, y = y, mode = 'lines', line = dict(width = 2, color = 'blue')),
           dict(x = x, y = y, mode = 'lines', line = dict(width = 2, color = 'blue'))]
    layout = dict(xaxis = dict(range = [xm, xM], autorange = False, zeroline = False, title = "X"), 
                  yaxis = dict(range = [ym, yM], autorange = False, zeroline = False, title = "Loss"),
                  title = "Plot of function curve",
                  hovermode = 'closest',
                  updatemenus= [{'type': 'buttons',
                           'buttons': [{'label': 'Play',
                                        'method': 'animate',
                                        'args': [None]}]}])
    N = len(x)
    mod = 1
    if N > 20:
        mod = int(N / 20)
    
    frames = [dict(data = [dict(x = [x[k]], y = [y[k]], 
                   mode = 'markers',
                   marker = dict(color='red', size=10))]) for k in range(N) if k % mod == 0]
    figure1 = dict(data = data, layout = layout, frames = frames)
    return iplot(figure1)
    
```


```python
plot_path(xpath, ypath)
```

## Examples

### Linear regression OLS

Let's first start with a simple example, the coefficient estimation in Linear regression (OLS)

\begin{align}
    Y &= X\beta + \epsilon, \ \epsilon \sim N(0, \sigma^2 I) \quad \text{// Statistical Model} \\
    \text{argmin}_{\theta}\ L(\beta | X) &= \frac{1}{2N}(Y - X\beta)^T(Y - X\beta) \quad \text{// Negative log-likelihood} \\
    &= \frac{1}{2N} (Y^TY - 2Y^TX\beta + \beta^T X^T X \beta)\\
\end{align}

$N$ is the sample size

- OLS result: $\hat{\beta}\ = (X^TX)^{-1}X^TY$
- Gradient descent:
\begin{align}
    &\nabla L(\beta_i | X) = \frac{-X^T(Y - X \beta_i)}{N} = \frac{-X^T Y + X^TX\beta_i}{N}\quad \text{// Gradient of } \beta_i \\
    &\theta_{i + 1} = \beta_i -\gamma\nabla L(\beta_i | x) \quad \text{// Update step}\\
    &\text{Check }\nabla L(\beta_i | X)
\end{align}

### Data

We use a simple example to perform the task. Here we set `N = 100000` and `p = 5`. 

$$\beta = \begin{bmatrix} 1\\ 2\\ 3\\ -2\\ -3\end{bmatrix} \quad \text{True}$$
$$X \sim N(0, I)$$


```python
N = 100000
p = 5
beta_true = np.array([1., 2.0, 3.0, -2.0, -3.0])
x = np.random.randn(N, p)
y = x @ beta_true + np.random.randn(N)
```


```python
def lm(X, y):
    return np.linalg.inv(X.T @ X) @ X.T @ y
```


```python
%time
lm(x, y)
```


```python
def loss(beta, x, y):
    N = len(y)
    xtx = x.T @ x
    xty = x.T @ y
    yty = y.T @ y
    return (yty - 2 * beta.T @ xty + beta.T @ xtx @ beta) / (2 * N)

def derivative(beta, x, y):
    N = len(y)
    xtx = x.T @ x
    xty = x.T @ y
    return (-xty + xtx @ beta) / N
```


```python
def gdescent_lm(init, x, y, derivative, loss,
                precision = 0.001, gamma = 0.001,
                previous_step_size = 1, max_iters = 10000):
    
    cur_x = init
    loss_path = [loss(cur_x, x, y)]
    
    iters = 1 #iteration counter
    while previous_step_size > precision and iters < max_iters:
        
        delta = # write yourself
        cur_x = # update step
        
        previous_step_size = np.linalg.norm(delta, ord=2)
        
        cur_y = loss(cur_x, x, y)
        
        if cur_y > loss_path[-1]:
            print("Current stepsize {} maybe too large, change to {}".format(gamma, gamma * 0.9))
            gamma *= 0.9
        
        iters += 1
        if iters % 50 == 0:
            print("Current iteration: {}".format(iters))
        
        loss_path.append(cur_y)

    return cur_x, loss_path, iters
```


```python
%time beta_gd, loss_path_gd, steps_gd = gdescent_lm(np.array([0., 0., 0., 0., 0.]), x, y, derivative, loss, gamma = 0.01)
```


```python
plot_path([i for i in range(steps)], loss_path)
```

# Stochastic Gradient Descent

`@` is a very expensive step, also `.T` in the gradient evaluation. Because it involves summation of all the data points. Stochastic Gradient Descent only evaluate one data point at once or a batch ($\ll N$) of data points at once. 

## Steps:

1. Set batch i = some numbers. The batch means a chunk of (x, y) where the batch size usually $\ll N$
2. Calculate the gradient $\nabla L(\theta_i | X_{\text{batch i}})$
3. Move the negative step of the $-\gamma\nabla L(\theta_i | X_{\text{batch i}})$
4. Update $\theta_{i + 1} = \theta_i -\gamma\nabla L(\theta_i | X_{\text{batch i}})$
5. Check convergence
6. get another batch, go back to step 1


```python
def sgd_lm(init, x, y, derivative, loss, 
           batchSize = 500, epoch = 10,
           precision = 0.001, gamma = 0.001,
           previous_step_size = 1):
    N = len(y)
    if N < batchSize:
        print("Current sample size: {}, batch size: {}\nUpdate batch size to 10% samples: {}".format(N, batchSize, N // 10))
        batchSize = N // 10
    
    cur_x = init
    e = 0 ## Start from first epoch
    indexN = np.arange(N)

    while e < epoch and previous_step_size > precision:
        
        np.random.shuffle(indexN) ## give a random sequence of index
        parts = np.array_split(indexN, N // batchSize)
        cur_y = 0
        
        preloss = loss_path[-1] 
        
        for index in parts:
            
            delta = # write yourself
            cur_x = # update step
            curloss = loss(cur_x, x[index,], y[index]) 
            cur_y += curloss * len(index)
            
            
            if curloss > preloss:
                # print("Current stepsize {} maybe too large, change to {}".format(gamma, gamma * 0.99))
                gamma *= 0.99
            loss_path.append(curloss)
        #loss_path.append(cur_y / N)
        
        previous_step_size = np.linalg.norm(delta, ord=2)
        print("Finished epoch: {}".format(e))
                
        e += 1
    return cur_x, loss_path, e
```


```python
%time beta_sgd, loss_path_sgd, steps_sgd = sgd_lm(np.array([0., 0., 0., 0., 0.]), x, y, derivative, loss, batchSize=100, gamma = 0.01)
```


```python
plot_path([i for i in range(len(loss_path))], loss_path)
```


```python
beta
```

# Practice

Consider use gd and sgd for ridge regression? or neural network

You basically only need to create both `derivative` and `loss` function for ridge. Also consider the tunning parameter.


```python
def ridge_derivative(beta, x, y, tunning = 1):
    """write you ridge derivative"""
def ridge_loss(beta, x, y, tunning = 1):
    """write you ridge loss"""
```


```python
beta_sgd, loss_path_sgd, steps_sgd = sgd_lm(np.array([0., 0., 0., 0., 0.]), x, y, ridge_derivative, ridge_loss, batchSize=100, gamma = 0.01)
```

8b7a238d7b34194d9cac035e601754e8ae9a0cba
ipynb
Jupyter Notebook
_site/assets/python/SGD-discussion.ipynb
BaconZhou/stat628
fd154bbdcd6070a70010c2940e0dadab36ccf125
MIT
2020-01-15T16:20:48.000Z
2020-01-15T16:20:48.000Z
assets/python/SGD-discussion.ipynb
BaconZhou/stat628
fd154bbdcd6070a70010c2940e0dadab36ccf125
MIT
assets/python/SGD-discussion.ipynb
BaconZhou/stat628
fd154bbdcd6070a70010c2940e0dadab36ccf125
MIT
Qwen/Qwen-72B
1. YES
2. YES