# day 12: Stochastic Gradient Descent

# Objectives

* See how Stochastic Gradient Descent can be used on a simple model (1-dim. linear regression)
* See how SGD can be used on a complex model (MLPClassifier)
* Understand impact of batch size and learning rate (aka step size)


# Outline
* [Review: Loss and Gradient for 1-dim. Linear Regression](#part1)
* [Part 1: Stochastic Estimates of the Loss](#part1)
* [Part 2: Stochastic Estimates of the Gradient](#part2)
* [Part 3: Stochastic Gradient Descent Algorithm in a few lines of Python](#part3)
* [Part 4: Using sklearn code to train MLPClassifier with SGD](#part4)

# Takeaways

* Stochastic estimates of loss functions are possible when the function is *additive* over training examples

* Stochastic gradient descent is a simple algorithm that can be implemented in a few lines of Python
* * Practical issues include selecting step size and batch size $B$

* Selecting batch size trades off two things:
* * Runtime cost of computing each gradient estimate (scales with $O(B)$, so smaller is better)
* * Quality of the estimate (larger $B$ leads to less variance)


```python
import numpy as np
```


```python
import sklearn.neural_network
```


```python
# import plotting libraries
import matplotlib
import matplotlib.pyplot as plt

%matplotlib inline
plt.style.use('seaborn') # pretty matplotlib plots

import seaborn as sns
sns.set('notebook', font_scale=1.25, style='whitegrid')
```

# Create simple dataset:   y = 1.234 * x + noise

We will *intentionally* create a toy dataset where we know that a good solution has slope near 1.234.

We'll generate N = 1000 examples.

Naturally, the best slope for the finite dataset we create won't be exactly 1.234 (because of the noise added plus the fact that our dataset size is limited).


```python
def create_dataset(N=1000, slope=1.234, noise_stddev=0.1, random_state=0):
    random_state = np.random.RandomState(int(random_state))

    # input features
    x_N = np.linspace(-2, 2, N)
    
    # output features
    y_N = slope * x_N + random_state.randn(N) * noise_stddev
    
    return x_N, y_N
```


```python
xtrain_N, ytrain_N = create_dataset(N=1000, noise_stddev=0.3)
```


```python
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(5,5))
plt.plot(xtrain_N, ytrain_N, 'k.');
plt.xlabel('x');
plt.ylabel('y');
```

# Review: Gradient Descent for 1-dim. Linear Regression

## Define prediction model

Consider the *simplest* linear regression model. A single weight parameter $w \in \mathbb{R}$ representing the slope of the prediction line. No bias/intercept.

To make predictions, we just compute the weight multiplied by the input feature
$$
\hat{y}(x) = w \cdot x
$$

## Define loss function

We want to minimize the total *squared error* across all N observed data examples (input features $x_n$, output responses $y_n$)

Given a full dataset of $N$ examples, we can compute the loss as:

\begin{align}
    \min_{w \in \mathbb{R}} ~~ &\ell(w)
    \\
    \text{calc_loss}(w) = \ell(w) &= \frac{1}{N} \sum_{n=1}^N \frac{1}{2} (y_n - w x_n)^2
\end{align}

Given a random minibatch of $B$ examples, we can *estimate* the loss as:

\begin{align}
    \text{estimate_loss}(w) = &= \frac{1}{B} \sum_{n=1}^B \frac{1}{2} (y_n - w x_n)^2
\end{align}

### Exercise 1A: Complete the code below

You should make it match the math expression above.


```python
def calc_loss(w, xbatch_B, ybatch_B):
    ''' Compute loss for slope-only least-squares linear regression
    
    Will compute the estimate of the loss given a minibatch of x,y values
    
    Args
    ----
    w : float
        Value of slope parameter

    Returns
    -------
    loss : float
        Mean of squared error loss at provided w value
    '''
    B = xbatch_B.shape[0]
    yhat_B = xbatch_B * w
    half_mean_squared_error = 0.5 * np.mean(np.square(ybatch_B - yhat_B))
    return half_mean_squared_error
```

### Define the gradient function

Given a full dataset of $N$ examples, we can compute the gradient as:

\begin{align}
\text{calc_grad}(w) = \ell'(w) &= \frac{1}{N} \frac{\partial}{\partial w} [ \sum_{n=1}^N \frac{1}{2} (y_n - w x_n)^2] 
\\
&= \frac{1}{N}  \sum_{n=1}^N (y_n - w x_n) (-x_n)
\\
&= \frac{1}{N} \sum_{n=1}^N (w x_n - y_n) (x_n)
\\
&= w \left( \frac{1}{N} \sum_{n=1}^N x_n^2 \right) - \frac{1}{N} \sum_{n=1}^N y_n x_n
\end{align}

Given a random minibatch of $B$ examples, we can *estimate* the gradient using:

\begin{align}
\text{estimate_grad}(w) &= \frac{1}{B} \frac{\partial}{\partial w} [ \sum_{n=1}^B \frac{1}{2} (y_n - w x_n)^2] 
\\
&= w \left( \frac{1}{B} \sum_{n=1}^B x_n^2 \right) - \frac{1}{B} \sum_{n=1}^B y_n x_n
\end{align}

Below, we've implemented the gradient calculation in code for you


```python
def calc_grad(w, xbatch_B, ybatch_B):
    ''' Compute gradient for slope-only least-squares linear regression
    
    Will compute a deterministic estimate of the gradient given a minibatch of x,y values

    Args
    ----
    w : float
        Value of slope parameter

    Returns
    -------
    g : float
        Value of derivative of loss function at provided w value
    '''
    g = w * np.mean(np.square(xbatch_B)) - np.mean(xbatch_B * ybatch_B)
    return g
```

# Part 1: Stochastic estimates of the loss

### Plot whole-dataset loss evaluated at each w from -3 to 8

We should see a "bowl" shape with one *global* minima, because our optimization problem is "convex"


```python
G = 101
w_grid = np.linspace(-3, 8, G) # create array of 300 values between -3 and 8
print(w_grid)
```

    [-3.   -2.89 -2.78 -2.67 -2.56 -2.45 -2.34 -2.23 -2.12 -2.01 -1.9  -1.79
     -1.68 -1.57 -1.46 -1.35 -1.24 -1.13 -1.02 -0.91 -0.8  -0.69 -0.58 -0.47
     -0.36 -0.25 -0.14 -0.03  0.08  0.19  0.3   0.41  0.52  0.63  0.74  0.85
      0.96  1.07  1.18  1.29  1.4   1.51  1.62  1.73  1.84  1.95  2.06  2.17
      2.28  2.39  2.5   2.61  2.72  2.83  2.94  3.05  3.16  3.27  3.38  3.49
      3.6   3.71  3.82  3.93  4.04  4.15  4.26  4.37  4.48  4.59  4.7   4.81
      4.92  5.03  5.14  5.25  5.36  5.47  5.58  5.69  5.8   5.91  6.02  6.13
      6.24  6.35  6.46  6.57  6.68  6.79  6.9   7.01  7.12  7.23  7.34  7.45
      7.56  7.67  7.78  7.89  8.  ]



```python
loss_grid = np.zeros(G)
for gg in range(G):
    loss_grid[gg] = calc_loss(w_grid[gg], xtrain_N, ytrain_N)

plt.plot(w_grid, loss_grid, 'b.-');
plt.xlabel('w');
plt.ylabel('loss(w)');
plt.ylim([0, 35]);
```

### Discussion 1a: Visually, at what value of $w$ does the loss function have a minima? Is it near where you would expect (hint: look above for the "true" slope value used to generate the data)


```python

```

### Sampling a minibatch of size B

We have provided some starter code that samples a minibatch of B examples from a training set of N examples

Uses a provided `random_state` pseudo-random number generator, which defaults to numpy's if not specified.


```python
def draw_minibatch(xtrain_N, ytrain_N, batch_size=100, random_state=np.random):
    ''' Sample a minibatch of desired size from provided training set
    
    Returns
    -------
    xbatch_B : 1D array, size (B,)
        x values of minibatch
    ybatch_B : 1D array, size (B,)
        y values of minibatch
    '''
    N = ytrain_N.size
    selected_row_ids = random_state.choice(np.arange(N), size=batch_size, replace=False)
    xbatch_B = xtrain_N[selected_row_ids].copy()
    ybatch_B = ytrain_N[selected_row_ids].copy()
    return xbatch_B, ybatch_B
```

### Show several minibatches of size 1


```python
draw_minibatch(xtrain_N, ytrain_N, 1)
```




    (array([-1.76376376]), array([-2.28530683]))




```python
draw_minibatch(xtrain_N, ytrain_N, 1)
```




    (array([1.42342342]), array([1.62548]))




```python
draw_minibatch(xtrain_N, ytrain_N, 1)
```




    (array([-1.35935936]), array([-1.82685918]))



### Show example minibatch of size 500



```python
xbatch_B, ybatch_B = draw_minibatch(xtrain_N, ytrain_N, 500)

print("Showing first 5 entries of x array of shape %s" % str(xbatch_B.shape))
print(xbatch_B[:5])
print("Showing first 5 entries of y array of shape %s" % str(ybatch_B.shape))
print(ybatch_B[:5])
```

    Showing first 5 entries of x array of shape (500,)
    [ 1.67967968  0.41041041 -0.42642643 -0.04204204 -0.96296296]
    Showing first 5 entries of y array of shape (500,)
    [ 2.0268924   0.2231359  -1.01115247  0.01427242 -1.00518248]


### Exercise 1b: Would you expect the loss on the minibatch to be better at $w=1.0$ or at $w = -1.0$?

TODO write answer here to guess, then use code below to figure it out.

Show the loss at this random minibatch when $w = 1.0$


```python
calc_loss(1.0, xbatch_B, ybatch_B)
```




    0.08034805534465543



Show the loss at this random minibatch when w = -2.0


```python
calc_loss(-1.0, xbatch_B, ybatch_B)
```




    3.290188505364923



### Exercise 1c: Can you draw a random minibatch of size 25 and display *all* of it?


```python
xbatch_25, ybatch_25 = draw_minibatch(xtrain_N, ytrain_N, 25)

print("Showing x array")
print(xbatch_25)
print("Showing y array")
print(ybatch_25)

```

    Showing x array
    [ 1.83983984 -0.03803804  0.90690691 -0.63863864  1.54754755 -1.25125125
      1.11911912 -0.8028028   1.35935936  1.83583584  1.75575576  1.23523524
      0.02602603  1.62362362 -0.47047047 -1.04304304 -0.37037037  1.003003
     -0.50650651  1.27127127  0.8028028   1.57957958  0.68268268  1.38338338
      0.25025025]
    Showing y array
    [ 1.88661403 -0.35591952  0.83928641 -0.46752726  2.19059324 -1.52839452
      1.36085275 -0.87964189  1.77088357  2.96941553  1.9289166   1.01781127
     -0.45767424  2.55036863 -1.25722983 -1.0840427  -0.43935936  1.46951416
     -0.6998666   2.07332527  0.5571766   1.52262524  0.67871684  1.64914208
      0.36732972]


### Exercise 1d: What is the loss at your minibatch at $w=1.0$ and $w = 1.5$?


```python
calc_loss(1.0, xbatch_25, ybatch_25)
```




    0.09975696099948157




```python
calc_loss(-1.0, xbatch_25, ybatch_25)
```




    3.3594047098831754



## Plot: Compute the *stochastic estimate* of loss as a function of $w$, using batch size 50

We'll make 3 lines, for 3 separate *trials* of this procedure, so we can see how much each trial's curve might vary


```python
for trial in range(3):
    random_state = np.random.RandomState(trial) # set seed based on the trial id

    stoch_loss_grid = np.zeros(G)
    for gg in range(G):
        xbatch_B, ybatch_B = draw_minibatch(xtrain_N, ytrain_N, batch_size=50, random_state=random_state)
        stoch_loss_grid[gg] = calc_loss(w_grid[gg], xbatch_B, ybatch_B)

    plt.plot(w_grid, stoch_loss_grid, '.-', label='Trial %d' % trial);
plt.xlabel('w');
plt.ylabel('loss(w)');
plt.ylim([0, 35]);
```

## Exercise 1e: Compute stochastic estimate of loss as function of w, for batch_size = 5


```python
for trial in range(3):
    random_state = np.random.RandomState(trial) # set seed based on the trial id

    stoch_loss_grid = np.zeros(G)
    for gg in range(G):
        xbatch_B, ybatch_B = draw_minibatch(xtrain_N, ytrain_N, batch_size=5, random_state=random_state) # TODO use draw_minibatch
        stoch_loss_grid[gg] = calc_loss(w_grid[gg], xbatch_B, ybatch_B)             # TODO use calc_loss at the current batch

    plt.plot(w_grid, stoch_loss_grid, '.-', label='Trial %d' % trial);
    
    
plt.xlabel('w');
plt.ylabel('loss(w)');
```

## Exercise 1f: Compute stochastic estimate of loss as function of w, for batch_size = 1

Repeat the above for batch_size = 1


```python
for trial in range(3):
    random_state = np.random.RandomState(trial) # set seed based on the trial id

    stoch_loss_grid = np.zeros(G)
    for gg in range(G):
        xbatch_B, ybatch_B = draw_minibatch(xtrain_N, ytrain_N, batch_size=1, random_state=random_state) # TODO use draw_minibatch
        stoch_loss_grid[gg] = calc_loss(w_grid[gg], xbatch_B, ybatch_B)             # TODO use calc_loss at the current batch

    plt.plot(w_grid, stoch_loss_grid, '.-', label='Trial %d' % trial);
    
    
plt.xlabel('w');
plt.ylabel('loss(w)');
```

## Discussion 1f: What can you say about the *variance* of the stochastic estimates

How does variance change as a function of batch size ?

## Part 2: Stochastic estimates of the gradient

## Sanity check: plot whole dataset gradient evaluated at each w from -3 to 8


```python
grad_grid = np.zeros(G)
for gg in range(G):
    grad_grid[gg] = calc_grad(w_grid[gg], xtrain_N, ytrain_N)

plt.plot(w_grid, grad_grid, 'b.-');
plt.xlabel('w');
plt.ylabel('grad(w)');
plt.ylim([-7, 9]);
```

## Plot: Compute the *stochastic estimate* of grad as a function of $w$, using batch size 50


```python
for trial in range(3):
    random_state = np.random.RandomState(trial) # set seed based on the trial id

    stoch_grad_grid = np.zeros(G)
    for gg in range(G):
        xbatch_B, ybatch_B = draw_minibatch(xtrain_N, ytrain_N, batch_size=50, random_state=random_state)
        stoch_grad_grid[gg] = calc_grad(w_grid[gg], xbatch_B, ybatch_B)
        
    plt.plot(w_grid, stoch_grad_grid, '.-', label='Trial %d' % trial);
    

plt.xlabel('w');
plt.ylabel('grad(w)');
plt.ylim([-7, 9]);
```

## Exercise 2a: Repeat the above plot at batch_size = 5



```python
for trial in range(3):
    random_state = np.random.RandomState(trial) # set seed based on the trial id

    stoch_grad_grid = np.zeros(G)
    for gg in range(G):
        xbatch_B, ybatch_B = draw_minibatch(xtrain_N, ytrain_N, batch_size=5, random_state=random_state)
        stoch_grad_grid[gg] = calc_grad(w_grid[gg], xbatch_B, ybatch_B)
        
    plt.plot(w_grid, stoch_grad_grid, '.-', label='Trial %d' % trial);
    

plt.xlabel('w');
plt.ylabel('grad(w)');
plt.ylim([-7, 9]);
```

## Exercise 2c: Repeat the above plot at batch_size = 1


TODO copy code from above.


```python
for trial in range(3):
    random_state = np.random.RandomState(trial) # set seed based on the trial id

    stoch_grad_grid = np.zeros(G)
    for gg in range(G):
        xbatch_B, ybatch_B = draw_minibatch(xtrain_N, ytrain_N, batch_size=1, random_state=random_state)
        stoch_grad_grid[gg] = calc_grad(w_grid[gg], xbatch_B, ybatch_B)
        
    plt.plot(w_grid, stoch_grad_grid, '.-', label='Trial %d' % trial);
    

plt.xlabel('w');
plt.ylabel('grad(w)');
plt.ylim([-7, 9]);
```

## Discussion 2d: What happens to the variance of the grad estimate as batch_size increases?


## Discussion 2d: What happens to the cost of computing the estimate as batch_size increases?


## Part 3: Stochastic Gradient Descent (GD) as an algorithm in Python

### Define minimize_via_sgd algorithm

Can you understand what each step of this algorithm does?


```python
def minimize_via_sgd(xtrain_N, ytrain_N, init_w=0.0, batch_size=10, step_size=0.001, max_iters=100, random_state=0):
    ''' Perform minimization of provided loss function via gradient descent
    
    Each "iteration" performs one or more gradient updates, until total training set has been "seen"
    
    Args
    ----
    xtrain_N : numpy array, shape (N,)
    ytrain_N : numpy array, shape (N,)
    init_w : float
    batch_size : int
    step_size : float
    max_iters : positive int
    
    Return
    ----
    wopt: float
        array of optimized weights that approximately gives the least error
    info_dict : dict
        Contains information about the optimization procedure useful for debugging
        Entries include:
        * trace_loss_list : list of loss values
        * trace_grad_list : list of gradient values
    '''
    N = int(ytrain_N.size)
    B = int(batch_size)
    
    if isinstance(random_state, int):
        random_state = np.random.RandomState(random_state)

    w = 1.0 * init_w 
    
    # Create some lists to track progress over time (for debugging)
    trace_loss_list = []
    trace_w_list = []
    trace_grad_list = []

    step_id = 0
    for iter_id in range(max_iters):
        
        n_examples_seen_this_iter = 0
        while n_examples_seen_this_iter < N:
            xbatch_B, ybatch_B = draw_minibatch(xtrain_N, ytrain_N, batch_size, random_state=random_state)
            n_examples_seen_this_iter += batch_size
            
            loss = calc_loss(w, xbatch_B, ybatch_B)    
            grad = calc_grad(w, xbatch_B, ybatch_B)    
            w = w - step_size * grad
            step_id += 1

        print("  iter %5d/%d done | step %5d | w  % 13.5f | loss % 13.4f | grad % 13.4f" % (
            iter_id, max_iters, step_id, w, loss, grad))
    
        trace_loss_list.append(loss)
        trace_w_list.append(w)
        trace_grad_list.append(grad)
    
    wopt = w
    info_dict = dict(
        trace_loss_list=trace_loss_list,
        trace_w_list=trace_w_list, 
        trace_grad_list=trace_grad_list)
    
    return wopt, info_dict
```

### Discussion 2a: Which line of the above function does the *parameter update* happen?

Remember, in math, the parameter update of gradient descent is this:
$$
w \gets w - \alpha \nabla_w \ell(w)
$$

where $\alpha > 0$ is the step size.

In words, this math says *move* the parameter $w$ from its current value a *small step* in the "downhill" direction (indicated by gradient).

TODO write down here which line above *you* think it is, then discuss with your group


```python

```

### Try it! Run SGD with step_size = 0.01 and batch_size = 200

Running the cell below will have the following effects:

1) one line will be printed for every iteration, indicating the current w value and its associated loss

2) the "optimal" value of w will be stored in the variable named `wopt` returned by this function

3) a dictionary of information useful for debugging will be stored in the `info_dict` returned by this function


```python
wopt, info_dict = minimize_via_sgd(xtrain_N, ytrain_N, step_size=0.01, batch_size=200);
```

      iter     0/100 done | step     5 | w        0.08015 | loss        1.0488 | grad       -1.7202
      iter     1/100 done | step    10 | w        0.15328 | loss        0.8576 | grad       -1.4773
      iter     2/100 done | step    15 | w        0.22259 | loss        0.7444 | grad       -1.3535
      iter     3/100 done | step    20 | w        0.28735 | loss        0.6190 | grad       -1.2232
      iter     4/100 done | step    25 | w        0.35058 | loss        0.5804 | grad       -1.2300
      iter     5/100 done | step    30 | w        0.40936 | loss        0.5670 | grad       -1.2455
      iter     6/100 done | step    35 | w        0.46134 | loss        0.4426 | grad       -1.0366
      iter     7/100 done | step    40 | w        0.51053 | loss        0.3971 | grad       -0.9791
      iter     8/100 done | step    45 | w        0.55889 | loss        0.3367 | grad       -0.8819
      iter     9/100 done | step    50 | w        0.60231 | loss        0.3295 | grad       -0.8519
      iter    10/100 done | step    55 | w        0.64268 | loss        0.2927 | grad       -0.8156
      iter    11/100 done | step    60 | w        0.68140 | loss        0.2557 | grad       -0.7506
      iter    12/100 done | step    65 | w        0.71875 | loss        0.2472 | grad       -0.7419
      iter    13/100 done | step    70 | w        0.75133 | loss        0.1835 | grad       -0.5930
      iter    14/100 done | step    75 | w        0.78277 | loss        0.1722 | grad       -0.5639
      iter    15/100 done | step    80 | w        0.81064 | loss        0.1731 | grad       -0.5796
      iter    16/100 done | step    85 | w        0.83925 | loss        0.1716 | grad       -0.6222
      iter    17/100 done | step    90 | w        0.86445 | loss        0.1464 | grad       -0.5074
      iter    18/100 done | step    95 | w        0.88863 | loss        0.1269 | grad       -0.4665
      iter    19/100 done | step   100 | w        0.91101 | loss        0.1231 | grad       -0.4666
      iter    20/100 done | step   105 | w        0.93204 | loss        0.1004 | grad       -0.4243
      iter    21/100 done | step   110 | w        0.95136 | loss        0.0874 | grad       -0.3560
      iter    22/100 done | step   115 | w        0.97035 | loss        0.0987 | grad       -0.3765
      iter    23/100 done | step   120 | w        0.98704 | loss        0.0808 | grad       -0.2952
      iter    24/100 done | step   125 | w        1.00346 | loss        0.0852 | grad       -0.3169
      iter    25/100 done | step   130 | w        1.01760 | loss        0.0784 | grad       -0.3205
      iter    26/100 done | step   135 | w        1.03255 | loss        0.0701 | grad       -0.2819
      iter    27/100 done | step   140 | w        1.04555 | loss        0.0659 | grad       -0.2295
      iter    28/100 done | step   145 | w        1.05768 | loss        0.0687 | grad       -0.2576
      iter    29/100 done | step   150 | w        1.06847 | loss        0.0614 | grad       -0.2003
      iter    30/100 done | step   155 | w        1.07939 | loss        0.0661 | grad       -0.2469
      iter    31/100 done | step   160 | w        1.08989 | loss        0.0582 | grad       -0.2022
      iter    32/100 done | step   165 | w        1.10010 | loss        0.0630 | grad       -0.2164
      iter    33/100 done | step   170 | w        1.10851 | loss        0.0565 | grad       -0.1732
      iter    34/100 done | step   175 | w        1.11634 | loss        0.0552 | grad       -0.1632
      iter    35/100 done | step   180 | w        1.12382 | loss        0.0565 | grad       -0.1650
      iter    36/100 done | step   185 | w        1.12964 | loss        0.0494 | grad       -0.1079
      iter    37/100 done | step   190 | w        1.13630 | loss        0.0470 | grad       -0.1289
      iter    38/100 done | step   195 | w        1.14222 | loss        0.0478 | grad       -0.1077
      iter    39/100 done | step   200 | w        1.14832 | loss        0.0590 | grad       -0.1447
      iter    40/100 done | step   205 | w        1.15326 | loss        0.0419 | grad       -0.1074
      iter    41/100 done | step   210 | w        1.15853 | loss        0.0460 | grad       -0.0721
      iter    42/100 done | step   215 | w        1.16395 | loss        0.0509 | grad       -0.1238
      iter    43/100 done | step   220 | w        1.16902 | loss        0.0426 | grad       -0.1063
      iter    44/100 done | step   225 | w        1.17258 | loss        0.0459 | grad       -0.0702
      iter    45/100 done | step   230 | w        1.17556 | loss        0.0466 | grad       -0.0392
      iter    46/100 done | step   235 | w        1.17890 | loss        0.0428 | grad       -0.0586
      iter    47/100 done | step   240 | w        1.18149 | loss        0.0433 | grad       -0.0558
      iter    48/100 done | step   245 | w        1.18422 | loss        0.0477 | grad       -0.0653
      iter    49/100 done | step   250 | w        1.18750 | loss        0.0428 | grad       -0.0637
      iter    50/100 done | step   255 | w        1.19068 | loss        0.0525 | grad       -0.0635
      iter    51/100 done | step   260 | w        1.19410 | loss        0.0478 | grad       -0.0637
      iter    52/100 done | step   265 | w        1.19610 | loss        0.0371 | grad       -0.0345
      iter    53/100 done | step   270 | w        1.19875 | loss        0.0444 | grad       -0.0442
      iter    54/100 done | step   275 | w        1.20074 | loss        0.0458 | grad       -0.0389
      iter    55/100 done | step   280 | w        1.20295 | loss        0.0467 | grad       -0.0212
      iter    56/100 done | step   285 | w        1.20470 | loss        0.0410 | grad       -0.0131
      iter    57/100 done | step   290 | w        1.20604 | loss        0.0386 | grad       -0.0209
      iter    58/100 done | step   295 | w        1.20676 | loss        0.0387 | grad       -0.0180
      iter    59/100 done | step   300 | w        1.20872 | loss        0.0407 | grad       -0.0393
      iter    60/100 done | step   305 | w        1.20974 | loss        0.0391 | grad       -0.0075
      iter    61/100 done | step   310 | w        1.21074 | loss        0.0448 | grad       -0.0445
      iter    62/100 done | step   315 | w        1.21162 | loss        0.0469 | grad        0.0055
      iter    63/100 done | step   320 | w        1.21239 | loss        0.0385 | grad        0.0217
      iter    64/100 done | step   325 | w        1.21317 | loss        0.0404 | grad       -0.0250
      iter    65/100 done | step   330 | w        1.21514 | loss        0.0422 | grad       -0.0427
      iter    66/100 done | step   335 | w        1.21667 | loss        0.0505 | grad       -0.0682
      iter    67/100 done | step   340 | w        1.21734 | loss        0.0441 | grad       -0.0089
      iter    68/100 done | step   345 | w        1.21843 | loss        0.0487 | grad       -0.0520
      iter    69/100 done | step   350 | w        1.21968 | loss        0.0426 | grad       -0.0217
      iter    70/100 done | step   355 | w        1.22061 | loss        0.0451 | grad       -0.0162
      iter    71/100 done | step   360 | w        1.22148 | loss        0.0494 | grad       -0.0131
      iter    72/100 done | step   365 | w        1.22243 | loss        0.0432 | grad       -0.0476
      iter    73/100 done | step   370 | w        1.22310 | loss        0.0359 | grad       -0.0035
      iter    74/100 done | step   375 | w        1.22355 | loss        0.0445 | grad       -0.0733
      iter    75/100 done | step   380 | w        1.22387 | loss        0.0475 | grad       -0.0229
      iter    76/100 done | step   385 | w        1.22515 | loss        0.0418 | grad       -0.0054
      iter    77/100 done | step   390 | w        1.22519 | loss        0.0457 | grad        0.0189
      iter    78/100 done | step   395 | w        1.22512 | loss        0.0411 | grad        0.0211
      iter    79/100 done | step   400 | w        1.22509 | loss        0.0423 | grad       -0.0424
      iter    80/100 done | step   405 | w        1.22611 | loss        0.0435 | grad       -0.0184
      iter    81/100 done | step   410 | w        1.22761 | loss        0.0454 | grad       -0.0352
      iter    82/100 done | step   415 | w        1.22761 | loss        0.0444 | grad        0.0206
      iter    83/100 done | step   420 | w        1.22788 | loss        0.0434 | grad        0.0002
      iter    84/100 done | step   425 | w        1.22852 | loss        0.0393 | grad       -0.0144
      iter    85/100 done | step   430 | w        1.22787 | loss        0.0438 | grad       -0.0081
      iter    86/100 done | step   435 | w        1.22866 | loss        0.0451 | grad       -0.0259
      iter    87/100 done | step   440 | w        1.22834 | loss        0.0410 | grad        0.0073
      iter    88/100 done | step   445 | w        1.22975 | loss        0.0406 | grad       -0.0768
      iter    89/100 done | step   450 | w        1.23034 | loss        0.0427 | grad       -0.0344
      iter    90/100 done | step   455 | w        1.23006 | loss        0.0496 | grad       -0.0083
      iter    91/100 done | step   460 | w        1.23052 | loss        0.0438 | grad       -0.0056
      iter    92/100 done | step   465 | w        1.23073 | loss        0.0441 | grad        0.0013
      iter    93/100 done | step   470 | w        1.23094 | loss        0.0471 | grad        0.0267
      iter    94/100 done | step   475 | w        1.23106 | loss        0.0392 | grad        0.0114
      iter    95/100 done | step   480 | w        1.23164 | loss        0.0459 | grad        0.0033
      iter    96/100 done | step   485 | w        1.23252 | loss        0.0481 | grad       -0.0301
      iter    97/100 done | step   490 | w        1.23304 | loss        0.0368 | grad       -0.0023
      iter    98/100 done | step   495 | w        1.23335 | loss        0.0429 | grad       -0.0010
      iter    99/100 done | step   500 | w        1.23343 | loss        0.0513 | grad       -0.0131



```python
fig, axes = plt.subplots(nrows=1, ncols=3, sharex=True, sharey=False, figsize=(12,3))

axes[0].plot(info_dict['trace_loss_list'], '.-');
axes[0].set_title('loss');
axes[1].plot(info_dict['trace_grad_list'], '.-');
axes[1].set_title('grad');
axes[2].plot(info_dict['trace_w_list'], '.-');
axes[2].set_title('w');
```

### Discussion 2b: Does it appear from the *parameter* values in trace above that the SGD procedure converged?

### Discussion 2c: Does it appear from the *loss* values in trace above that the SGD procedure converged?

### Discussion 2d: Does it appear from the *grad* values in trace above that the SGD procedure converged?


```python

```


```python

```

## Try with smaller batch_size = 10


```python
wopt, info_dict = minimize_via_sgd(xtrain_N, ytrain_N, step_size=0.01, batch_size=10);
```

      iter     0/100 done | step   100 | w        0.91625 | loss        0.1191 | grad       -0.3511
      iter     1/100 done | step   200 | w        1.14951 | loss        0.0377 | grad        0.0505
      iter     2/100 done | step   300 | w        1.20424 | loss        0.0991 | grad        0.1091
      iter     3/100 done | step   400 | w        1.23031 | loss        0.0381 | grad        0.0042
      iter     4/100 done | step   500 | w        1.23886 | loss        0.0511 | grad       -0.2389
      iter     5/100 done | step   600 | w        1.24595 | loss        0.0859 | grad       -0.1187
      iter     6/100 done | step   700 | w        1.25205 | loss        0.0334 | grad        0.1019
      iter     7/100 done | step   800 | w        1.23007 | loss        0.0334 | grad       -0.0147
      iter     8/100 done | step   900 | w        1.22628 | loss        0.0573 | grad        0.0265
      iter     9/100 done | step  1000 | w        1.23810 | loss        0.0412 | grad        0.0452
      iter    10/100 done | step  1100 | w        1.22357 | loss        0.0227 | grad        0.0010
      iter    11/100 done | step  1200 | w        1.23342 | loss        0.0515 | grad       -0.0504
      iter    12/100 done | step  1300 | w        1.24008 | loss        0.0484 | grad       -0.2045
      iter    13/100 done | step  1400 | w        1.24485 | loss        0.0738 | grad       -0.0568
      iter    14/100 done | step  1500 | w        1.23728 | loss        0.0317 | grad        0.0817
      iter    15/100 done | step  1600 | w        1.22575 | loss        0.0404 | grad        0.0695
      iter    16/100 done | step  1700 | w        1.22235 | loss        0.0293 | grad        0.1196
      iter    17/100 done | step  1800 | w        1.22595 | loss        0.0330 | grad       -0.1241
      iter    18/100 done | step  1900 | w        1.23493 | loss        0.0149 | grad        0.0567
      iter    19/100 done | step  2000 | w        1.23858 | loss        0.0361 | grad       -0.0242
      iter    20/100 done | step  2100 | w        1.23371 | loss        0.0423 | grad        0.0176
      iter    21/100 done | step  2200 | w        1.23280 | loss        0.0428 | grad       -0.0151
      iter    22/100 done | step  2300 | w        1.24703 | loss        0.0544 | grad       -0.0138
      iter    23/100 done | step  2400 | w        1.24159 | loss        0.0544 | grad        0.0364
      iter    24/100 done | step  2500 | w        1.24189 | loss        0.0278 | grad       -0.0217
      iter    25/100 done | step  2600 | w        1.23713 | loss        0.0509 | grad        0.0451
      iter    26/100 done | step  2700 | w        1.23292 | loss        0.0502 | grad        0.1502
      iter    27/100 done | step  2800 | w        1.23131 | loss        0.0266 | grad       -0.0233
      iter    28/100 done | step  2900 | w        1.23236 | loss        0.0267 | grad        0.0038
      iter    29/100 done | step  3000 | w        1.22569 | loss        0.0350 | grad       -0.0734
      iter    30/100 done | step  3100 | w        1.22976 | loss        0.0349 | grad        0.0029
      iter    31/100 done | step  3200 | w        1.23033 | loss        0.0173 | grad        0.0927
      iter    32/100 done | step  3300 | w        1.23814 | loss        0.0387 | grad        0.0252
      iter    33/100 done | step  3400 | w        1.23415 | loss        0.0182 | grad       -0.0323
      iter    34/100 done | step  3500 | w        1.23714 | loss        0.0549 | grad        0.1117
      iter    35/100 done | step  3600 | w        1.23755 | loss        0.0294 | grad       -0.0801
      iter    36/100 done | step  3700 | w        1.23231 | loss        0.0303 | grad        0.1700
      iter    37/100 done | step  3800 | w        1.23270 | loss        0.0449 | grad        0.0136
      iter    38/100 done | step  3900 | w        1.22366 | loss        0.0325 | grad        0.0768
      iter    39/100 done | step  4000 | w        1.22075 | loss        0.0216 | grad       -0.0362
      iter    40/100 done | step  4100 | w        1.23983 | loss        0.0784 | grad       -0.1227
      iter    41/100 done | step  4200 | w        1.22596 | loss        0.0632 | grad        0.1253
      iter    42/100 done | step  4300 | w        1.23038 | loss        0.0280 | grad       -0.1657
      iter    43/100 done | step  4400 | w        1.23377 | loss        0.0390 | grad       -0.1455
      iter    44/100 done | step  4500 | w        1.22033 | loss        0.0512 | grad       -0.1874
      iter    45/100 done | step  4600 | w        1.22807 | loss        0.0153 | grad        0.0763
      iter    46/100 done | step  4700 | w        1.23235 | loss        0.0493 | grad       -0.2164
      iter    47/100 done | step  4800 | w        1.22289 | loss        0.0420 | grad       -0.1002
      iter    48/100 done | step  4900 | w        1.22318 | loss        0.0355 | grad        0.0695
      iter    49/100 done | step  5000 | w        1.22875 | loss        0.0320 | grad        0.0094
      iter    50/100 done | step  5100 | w        1.21265 | loss        0.0573 | grad        0.0005
      iter    51/100 done | step  5200 | w        1.23306 | loss        0.0531 | grad       -0.1133
      iter    52/100 done | step  5300 | w        1.23758 | loss        0.0500 | grad        0.0407
      iter    53/100 done | step  5400 | w        1.23312 | loss        0.0659 | grad        0.0761
      iter    54/100 done | step  5500 | w        1.23100 | loss        0.0128 | grad       -0.0745
      iter    55/100 done | step  5600 | w        1.23669 | loss        0.0392 | grad        0.0078
      iter    56/100 done | step  5700 | w        1.24015 | loss        0.0410 | grad       -0.0476
      iter    57/100 done | step  5800 | w        1.22657 | loss        0.0248 | grad        0.0298
      iter    58/100 done | step  5900 | w        1.22985 | loss        0.0573 | grad        0.0474
      iter    59/100 done | step  6000 | w        1.22832 | loss        0.0259 | grad       -0.0250
      iter    60/100 done | step  6100 | w        1.23766 | loss        0.0662 | grad       -0.0484
      iter    61/100 done | step  6200 | w        1.23839 | loss        0.0862 | grad       -0.2737
      iter    62/100 done | step  6300 | w        1.23734 | loss        0.0323 | grad       -0.0365
      iter    63/100 done | step  6400 | w        1.23474 | loss        0.0286 | grad       -0.0304
      iter    64/100 done | step  6500 | w        1.24034 | loss        0.0578 | grad        0.1138
      iter    65/100 done | step  6600 | w        1.23524 | loss        0.0212 | grad        0.0872
      iter    66/100 done | step  6700 | w        1.22606 | loss        0.0438 | grad        0.1876
      iter    67/100 done | step  6800 | w        1.23196 | loss        0.0306 | grad       -0.0178
      iter    68/100 done | step  6900 | w        1.23354 | loss        0.0525 | grad       -0.0704
      iter    69/100 done | step  7000 | w        1.22663 | loss        0.0442 | grad       -0.0648
      iter    70/100 done | step  7100 | w        1.22977 | loss        0.0419 | grad        0.0042
      iter    71/100 done | step  7200 | w        1.23995 | loss        0.0629 | grad        0.0045
      iter    72/100 done | step  7300 | w        1.22573 | loss        0.0197 | grad        0.0462
      iter    73/100 done | step  7400 | w        1.22610 | loss        0.0434 | grad        0.1277
      iter    74/100 done | step  7500 | w        1.22869 | loss        0.0187 | grad       -0.0309
      iter    75/100 done | step  7600 | w        1.22590 | loss        0.0532 | grad        0.0942
      iter    76/100 done | step  7700 | w        1.22963 | loss        0.0493 | grad        0.1311
      iter    77/100 done | step  7800 | w        1.23985 | loss        0.0852 | grad        0.0897
      iter    78/100 done | step  7900 | w        1.23384 | loss        0.0514 | grad        0.0544
      iter    79/100 done | step  8000 | w        1.23430 | loss        0.0295 | grad       -0.1374
      iter    80/100 done | step  8100 | w        1.22765 | loss        0.0141 | grad       -0.0433
      iter    81/100 done | step  8200 | w        1.23420 | loss        0.0216 | grad       -0.0395
      iter    82/100 done | step  8300 | w        1.22563 | loss        0.0665 | grad       -0.1015
      iter    83/100 done | step  8400 | w        1.24080 | loss        0.0364 | grad        0.0033
      iter    84/100 done | step  8500 | w        1.23774 | loss        0.0632 | grad       -0.0353
      iter    85/100 done | step  8600 | w        1.23484 | loss        0.0781 | grad       -0.0488
      iter    86/100 done | step  8700 | w        1.23275 | loss        0.0117 | grad       -0.0342
      iter    87/100 done | step  8800 | w        1.23411 | loss        0.0321 | grad        0.0455
      iter    88/100 done | step  8900 | w        1.22806 | loss        0.0688 | grad       -0.0674
      iter    89/100 done | step  9000 | w        1.23307 | loss        0.0222 | grad        0.1727
      iter    90/100 done | step  9100 | w        1.23491 | loss        0.0251 | grad       -0.0312
      iter    91/100 done | step  9200 | w        1.23106 | loss        0.0519 | grad        0.0098
      iter    92/100 done | step  9300 | w        1.23303 | loss        0.0437 | grad       -0.1625
      iter    93/100 done | step  9400 | w        1.23277 | loss        0.0708 | grad       -0.2593
      iter    94/100 done | step  9500 | w        1.22882 | loss        0.0636 | grad        0.0411
      iter    95/100 done | step  9600 | w        1.24039 | loss        0.0208 | grad       -0.0708
      iter    96/100 done | step  9700 | w        1.23499 | loss        0.0839 | grad        0.0251
      iter    97/100 done | step  9800 | w        1.23123 | loss        0.0467 | grad        0.0922
      iter    98/100 done | step  9900 | w        1.23677 | loss        0.0182 | grad        0.0165
      iter    99/100 done | step 10000 | w        1.22925 | loss        0.0492 | grad        0.0620



```python
fig, axes = plt.subplots(nrows=1, ncols=3, sharex=True, sharey=False, figsize=(12,3))

axes[0].plot(info_dict['trace_loss_list'], '.-');
axes[0].set_title('loss');
axes[1].plot(info_dict['trace_grad_list'], '.-');
axes[1].set_title('grad');
axes[2].plot(info_dict['trace_w_list'], '.-');
axes[2].set_title('w');
```

### Try with even smaller batch_size of 1


```python
# TODO write code here
```

### Discussion 3b: What happens here with this smaller batch size? Is it converging?

# Part 5: SGD for MLPClassifier



Let's revisit the XOR dataset from our previous lab, and try SGD for it.


```python
def make_xor_dataset(n_per_blob=50, stddev=0.4, random_state=0):
    random_state = np.random.RandomState(random_state)
    cov_22 = np.square(stddev) * np.eye(2)
    x_00 = random_state.multivariate_normal([-1, -1], cov_22, size=n_per_blob)
    x_01 = random_state.multivariate_normal([-1, +1], cov_22, size=n_per_blob)
    x_10 = random_state.multivariate_normal([+1, -1], cov_22, size=n_per_blob)
    x_11 = random_state.multivariate_normal([+1, +1], cov_22, size=n_per_blob)

    N = n_per_blob * 4
    x_N2 = np.vstack([x_00, x_11, x_01, x_10])
    assert x_N2.shape == (N, 2)

    y_N = np.hstack([np.ones(N//2), np.zeros(N//2)]).astype(np.int32)
    assert y_N.shape == (N,)

    # Shuffle the order
    perm_ids = random_state.permutation(N)
    x_N2 = x_N2[perm_ids].copy()
    y_N = y_N[perm_ids].copy()

    return x_N2, y_N
```


```python
x_tr_N2, y_tr_N = make_xor_dataset(n_per_blob=50)
```


```python
plt.figure(figsize=(4,4))
plt.plot(x_tr_N2[y_tr_N==0,0], x_tr_N2[y_tr_N==0,1], 'rx', label='y=0', mew=2);
plt.plot(x_tr_N2[y_tr_N==1,0], x_tr_N2[y_tr_N==1,1], 'b+', label='y=1', mew=2);

plt.xlabel('x_1');
plt.ylabel('x_2');
plt.legend(bbox_to_anchor=(1.0, 0.5));
```

# Setup: Create Utility function for visualizing classifier predictions

You do NOT need to understand the details of this function. We'll just use it as is.


```python
def plot_pretty_probabilities_for_clf(
        clf,
        ax=None,
        x1_grid=(-2.1, 2.1, 50), x2_grid=(-2.1, 2.1, 50),
        x1_ticks=[-1, 0, 1], x2_ticks=[-1, 0, 1],
        do_show_colorbar=False,
        c_ticks=np.asarray([0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]),
        c_num_levels=100,
        ):
    ''' Display predicted probabilities from classifier as color contour plot

    Args
    ----
    clf : sklearn object that implements classifier API
    ax : matplotlib axis handle, or None
        If provided, will use axis handle as primary axis to plot on.
        If None, will use the current handle, or make new one as necessary.
    x1_grid : tuple-like or array-like
        If tuple of length 3, interpreted as args to np.linspace
        Otherwise, cast to array and assumed to be a 1d grid of x1 values
    x2_grid : tuple-like or array-like
        If tuple of length 3, interpreted as args to np.linspace
        Otherwise, cast to array and assumed to be a 1d grid of x2 values

    '''
    # Activate the current axis, if necessary
    if ax is None:
        cur_ax = plt.gca()
    else:
        cur_ax = ax
        plt.sca(cur_ax)

    # Create dense grid of x1 and x2 values
    # useful for visualizing probabilities
    if isinstance(x1_grid, tuple) and len(x1_grid) == 3:
        x1_grid = np.linspace(x1_grid[0], x1_grid[1], x1_grid[2])
    if isinstance(x2_grid, tuple) and len(x2_grid) == 3:
        x2_grid = np.linspace(x2_grid[0], x2_grid[1], x2_grid[2])
    x1_grid = np.asarray(x1_grid).flatten()
    x2_grid = np.asarray(x2_grid).flatten()

    c_levels = np.linspace(0.0, 1.0, c_num_levels)

    # Get regular grid of G x H points, where each point is an (x1, x2) location
    G = x1_grid.size
    H = x2_grid.size    
    x1_GH, x2_GH = np.meshgrid(x1_grid, x2_grid)
    
    # Combine the x1 and x2 values into one array
    # Flattened into M = G x H rows
    # Each row of x_M2 is a 2D vector [x_m1, x_m2]
    x_M2 = np.hstack([x1_GH.flatten()[:,np.newaxis], x2_GH.flatten()[:,np.newaxis]])
    
    # Predict proba for each point in the flattened grid
    yproba1_M = clf.predict_proba(x_M2)[:,1]
    # Reshape the M probas into the GxH 2D field
    yproba1_GH = np.reshape(yproba1_M, x1_GH.shape)
    # Contour plot
    cmap = plt.cm.RdYlBu
    my_contourf_h = plt.contourf(
            x1_GH, x2_GH, yproba1_GH,
            levels=c_levels, vmin=0, vmax=1.0,
            cmap=cmap, alpha=0.5)

    # Edit the ticks observed
    if x1_ticks is not None:
        plt.xticks(x1_ticks, x1_ticks);
    if x2_ticks is not None:
        plt.yticks(x2_ticks, x2_ticks);
    if do_show_colorbar:
        left, bottom, width, height = plt.gca().get_position().bounds
        cax = plt.gcf().add_axes([left + 1.1*width, bottom, 0.03, height])
        plt.colorbar(my_contourf_h, orientation='vertical', cax=cax, ticks=c_ticks);
        plt.sca(cur_ax);
```

## Try an MLP with 2 hidden units, batch size of 25

This time, we'll use SGD as built in solver


```python
mlp_2hidden_run5 = sklearn.neural_network.MLPClassifier(
    hidden_layer_sizes=[2],
    activation='relu',
    solver='sgd',
    learning_rate_init=0.1,
    random_state=5,
    batch_size=25,
    )
```


```python
# Fit the model to training data

mlp_2hidden_run5.fit(x_tr_N2, y_tr_N)
```




    MLPClassifier(batch_size=25, hidden_layer_sizes=[2], learning_rate_init=0.1,
                  random_state=5, solver='sgd')



### Visualize the results using our utility function



```python
plt.figure(figsize=(4,4))
plot_pretty_probabilities_for_clf(mlp_2hidden_run5, do_show_colorbar=True, ax=plt.gca());

plt.plot(x_tr_N2[y_tr_N==0,0], x_tr_N2[y_tr_N==0,1], 'rx', label='y=0', mew=2);
plt.plot(x_tr_N2[y_tr_N==1,0], x_tr_N2[y_tr_N==1,1], 'b+', label='y=1', mew=2);
```


```python
### Visualize the trace of the loss

plt.plot(mlp_2hidden_run5.loss_curve_, 'k.-')
```

# Try several trials with batch_size of 25


```python
for trial in [1, 2, 3, 4, 5]:
    mlp_2hidden_batchsize25 = sklearn.neural_network.MLPClassifier(
        hidden_layer_sizes=[2],
        activation='relu',
        solver='sgd',
        learning_rate_init=0.05,
        random_state=trial,
        batch_size=25,
        )
    
    # Fit the model to training data
    mlp_2hidden_batchsize25.fit(x_tr_N2, y_tr_N)
    
    # Visualize the trace of the loss
    plt.plot(mlp_2hidden_batchsize25.loss_curve_, '.-', label='trial %d' % trial)
    
plt.legend();
```

### Now try several trials with a batch_size of 1


```python
for trial in [1, 2, 3, 4, 5]:
    mlp_2hidden_batchsize1 = sklearn.neural_network.MLPClassifier(
        hidden_layer_sizes=[2],
        activation='relu',
        solver='sgd',
        learning_rate_init=0.05,
        random_state=trial,
        batch_size=1,
        )
    
    # Fit the model to training data
    mlp_2hidden_batchsize1.fit(x_tr_N2, y_tr_N)
    
    # Visualize the trace of the loss
    plt.plot(mlp_2hidden_batchsize1.loss_curve_, '.-', label='trial %d' % trial)
    
plt.legend();
```

### Discussion 4a: What is happening here?


```python

```


```python

```

12db9be85d85700d46ad654c87de4d0fe171a198
ipynb
Jupyter Notebook
labs/day12_StochasticGradientDescent.ipynb
brawnerquan/comp135-20f-assignments
9570c17b872b7334b0e5b86160d868e3b854c71a
MIT
labs/day12_StochasticGradientDescent.ipynb
brawnerquan/comp135-20f-assignments
9570c17b872b7334b0e5b86160d868e3b854c71a
MIT
labs/day12_StochasticGradientDescent.ipynb
brawnerquan/comp135-20f-assignments
9570c17b872b7334b0e5b86160d868e3b854c71a
MIT
Qwen/Qwen-72B
1. YES
2. YES