# day06: Gradient Descent for Linear Regression

# Objectives

* Learn how to fit weight parameters of Linear Regression to a simple dataset via gradient descent
* Understand impact of step size
* Understand impact of initialization


# Outline
* [Part 1: Loss and Gradient for 1-dim. Linear Regression](#part1)
* [Part 2: Gradient Descent Algorithm in a few lines of Python](#part2)
* [Part 3: Debugging with Trace Plots](#part3)
* [Part 4: Selecting the step size](#part4)
* [Part 5: Selecting the initialization](#part5)
* [Part 6: Using SciPy's built-in routines](#part6)

# Takeaways


* Gradient descent is a simple algorithm that can be implemented in a few lines of Python
* * Practical issues include selecting step size and initialization
* Step size matters a lot
* * Need to select carefully for each problem

* Initialization of the parameters can matter too!

* scipy offers some useful tools for gradient-based optimization
* * scipy's toolbox cannot do scalable "stochastic" methods (requires a modest size dataset, not too big)
* * "L-BFGS-B" method is highly recommended if you have your loss and gradient functions available



```python
import numpy as np
```


```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.

Naturally, the best slope for the finite dataset of N=100 examples 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=100, 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
x_N, y_N = create_dataset(N=50, noise_stddev=0.3)
```


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

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

## Define 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$)

\begin{align}
    \min_{w \in \mathbb{R}} ~~ &\ell(w)
    \\
    \text{calc_loss}(w) = \ell(w) &= \sum_{n=1}^N (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):
    ''' Compute loss for slope-only least-squares linear regression
    
    Args
    ----
    w : float
        Value of slope parameter

    Returns
    -------
    loss : float
        Sum of squared error loss at provided w value
    '''
    yhat_N = x_N * w
    sum_squared_error = np.sum((y_N - yhat_N)**2) # todo compute the sum of squared error between y and yhat
    return sum_squared_error
```

# Define the gradient function

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

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


```python
def calc_grad(w):
    ''' Compute gradient for slope-only least-squares linear regression
    
    Args
    ----
    w : float
        Value of slope parameter

    Returns
    -------
    g : float
        Value of derivative of loss function at provided w value
    '''
    g = 2.0 * w * np.sum(np.square(x_N)) - 2.0 * np.sum(x_N * y_N)
    return g
```

## Plot 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
w_grid = np.linspace(-3, 8, 300) # create array of 300 values between -3 and 8
```


```python
loss_grid = np.asarray([calc_loss(w) for w in w_grid])
plt.plot(w_grid, loss_grid, 'b.-');
plt.xlabel('w');
plt.ylabel('loss(w)');
```

### Discussion 1b: 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)

### Exercise 1c: Write NumPy code to identify which entry in the w_grid array corresponds to the lowest entry in the loss_grid array

Hint: use np.argmin


```python
# TODO write code here
index = np.argmin(loss_grid)
print(index)
print(w_grid[index])
```

    112
    1.120401337792643


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


```python
grad_grid = np.asarray([calc_grad(w) for w in w_grid])
plt.plot(w_grid, grad_grid, 'b.-');
plt.xlabel('w');
plt.ylabel('grad(w)');
```

### Discussion 1d: Visually, at what value of $w$ does the gradient function cross zero? Is it the same place as the location of the minimum in the loss above?

TODO interpret the graph above and write your answer here, then discuss with your group

### Exercise 1d: Numerically, at which value of w does grad_grid cross zero?

We might try to estimate numerically where the gradient crosses zero.

We could do this in a few steps:

1) Compute the distance from each gradient in `grad_grid` to 0.0 (we could use just absolute distance)

2) Find the index of `grad_grid` with smallest distance (using `np.argmin`)

3) Plug that index into `w_grid` to get the $w$ value corresponding to that zero-crossing


```python
dist_from_zero_G = np.abs(grad_grid - 0.0)

zero_cross_index = 0 # TODO fix me for step 2 above

print("Zero crossing occurs at w = %.4f" % w_grid[0]) # TODO fix me for step 3 above
```

    Zero crossing occurs at w = -3.0000


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


### Define minimize_via_grad_descent algorithm

Can you understand what each step of this algorithm does?


```python
def minimize_via_grad_descent(calc_loss, calc_grad, init_w=0.0, step_size=0.001, max_iters=100):
    ''' Perform minimization of provided loss function via gradient descent
    
    Args
    ----
    calc_loss : function
    calc_grad : function
    init_w : float
    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
    '''
    w = 1.0 * init_w 
    grad = calc_grad(w)

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

    for iter_id in range(max_iters):
        if iter_id > 0:
            w = w - step_size * grad
        
        loss = calc_loss(w)
        grad = calc_grad(w)    

        print("  iter %5d/%d | w  % 13.5f | loss % 13.4f | grad % 13.4f" % (
            iter_id, max_iters, 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 GD with step_size = 0.001

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_grad_descent(calc_loss, calc_grad, step_size=0.001);
```

      iter     0/100 | w        0.00000 | loss       93.3197 | grad     -156.5566
      iter     1/100 | w        0.15656 | loss       70.5104 | grad     -134.8304
      iter     2/100 | w        0.29139 | loss       53.5926 | grad     -116.1192
      iter     3/100 | w        0.40751 | loss       41.0446 | grad     -100.0047
      iter     4/100 | w        0.50751 | loss       31.7376 | grad      -86.1265
      iter     5/100 | w        0.59364 | loss       24.8345 | grad      -74.1743
      iter     6/100 | w        0.66781 | loss       19.7144 | grad      -63.8807
      iter     7/100 | w        0.73169 | loss       15.9168 | grad      -55.0156
      iter     8/100 | w        0.78671 | loss       13.1001 | grad      -47.3808
      iter     9/100 | w        0.83409 | loss       11.0110 | grad      -40.8055
      iter    10/100 | w        0.87489 | loss        9.4614 | grad      -35.1427
      iter    11/100 | w        0.91004 | loss        8.3121 | grad      -30.2657
      iter    12/100 | w        0.94030 | loss        7.4597 | grad      -26.0656
      iter    13/100 | w        0.96637 | loss        6.8274 | grad      -22.4483
      iter    14/100 | w        0.98882 | loss        6.3584 | grad      -19.3331
      iter    15/100 | w        1.00815 | loss        6.0106 | grad      -16.6501
      iter    16/100 | w        1.02480 | loss        5.7526 | grad      -14.3395
      iter    17/100 | w        1.03914 | loss        5.5612 | grad      -12.3495
      iter    18/100 | w        1.05149 | loss        5.4193 | grad      -10.6357
      iter    19/100 | w        1.06212 | loss        5.3140 | grad       -9.1597
      iter    20/100 | w        1.07128 | loss        5.2360 | grad       -7.8886
      iter    21/100 | w        1.07917 | loss        5.1781 | grad       -6.7938
      iter    22/100 | w        1.08597 | loss        5.1351 | grad       -5.8510
      iter    23/100 | w        1.09182 | loss        5.1032 | grad       -5.0390
      iter    24/100 | w        1.09686 | loss        5.0796 | grad       -4.3397
      iter    25/100 | w        1.10120 | loss        5.0621 | grad       -3.7375
      iter    26/100 | w        1.10493 | loss        5.0491 | grad       -3.2188
      iter    27/100 | w        1.10815 | loss        5.0394 | grad       -2.7721
      iter    28/100 | w        1.11092 | loss        5.0323 | grad       -2.3874
      iter    29/100 | w        1.11331 | loss        5.0270 | grad       -2.0561
      iter    30/100 | w        1.11537 | loss        5.0231 | grad       -1.7708
      iter    31/100 | w        1.11714 | loss        5.0201 | grad       -1.5250
      iter    32/100 | w        1.11866 | loss        5.0180 | grad       -1.3134
      iter    33/100 | w        1.11998 | loss        5.0164 | grad       -1.1311
      iter    34/100 | w        1.12111 | loss        5.0152 | grad       -0.9742
      iter    35/100 | w        1.12208 | loss        5.0143 | grad       -0.8390
      iter    36/100 | w        1.12292 | loss        5.0136 | grad       -0.7225
      iter    37/100 | w        1.12364 | loss        5.0132 | grad       -0.6223
      iter    38/100 | w        1.12427 | loss        5.0128 | grad       -0.5359
      iter    39/100 | w        1.12480 | loss        5.0125 | grad       -0.4615
      iter    40/100 | w        1.12526 | loss        5.0123 | grad       -0.3975
      iter    41/100 | w        1.12566 | loss        5.0122 | grad       -0.3423
      iter    42/100 | w        1.12600 | loss        5.0121 | grad       -0.2948
      iter    43/100 | w        1.12630 | loss        5.0120 | grad       -0.2539
      iter    44/100 | w        1.12655 | loss        5.0119 | grad       -0.2187
      iter    45/100 | w        1.12677 | loss        5.0119 | grad       -0.1883
      iter    46/100 | w        1.12696 | loss        5.0119 | grad       -0.1622
      iter    47/100 | w        1.12712 | loss        5.0118 | grad       -0.1397
      iter    48/100 | w        1.12726 | loss        5.0118 | grad       -0.1203
      iter    49/100 | w        1.12738 | loss        5.0118 | grad       -0.1036
      iter    50/100 | w        1.12749 | loss        5.0118 | grad       -0.0892
      iter    51/100 | w        1.12757 | loss        5.0118 | grad       -0.0768
      iter    52/100 | w        1.12765 | loss        5.0118 | grad       -0.0662
      iter    53/100 | w        1.12772 | loss        5.0118 | grad       -0.0570
      iter    54/100 | w        1.12777 | loss        5.0118 | grad       -0.0491
      iter    55/100 | w        1.12782 | loss        5.0118 | grad       -0.0423
      iter    56/100 | w        1.12787 | loss        5.0118 | grad       -0.0364
      iter    57/100 | w        1.12790 | loss        5.0118 | grad       -0.0314
      iter    58/100 | w        1.12793 | loss        5.0118 | grad       -0.0270
      iter    59/100 | w        1.12796 | loss        5.0118 | grad       -0.0233
      iter    60/100 | w        1.12798 | loss        5.0118 | grad       -0.0200
      iter    61/100 | w        1.12800 | loss        5.0118 | grad       -0.0172
      iter    62/100 | w        1.12802 | loss        5.0118 | grad       -0.0149
      iter    63/100 | w        1.12804 | loss        5.0118 | grad       -0.0128
      iter    64/100 | w        1.12805 | loss        5.0118 | grad       -0.0110
      iter    65/100 | w        1.12806 | loss        5.0118 | grad       -0.0095
      iter    66/100 | w        1.12807 | loss        5.0118 | grad       -0.0082
      iter    67/100 | w        1.12808 | loss        5.0118 | grad       -0.0070
      iter    68/100 | w        1.12808 | loss        5.0118 | grad       -0.0061
      iter    69/100 | w        1.12809 | loss        5.0118 | grad       -0.0052
      iter    70/100 | w        1.12810 | loss        5.0118 | grad       -0.0045
      iter    71/100 | w        1.12810 | loss        5.0118 | grad       -0.0039
      iter    72/100 | w        1.12810 | loss        5.0118 | grad       -0.0033
      iter    73/100 | w        1.12811 | loss        5.0118 | grad       -0.0029
      iter    74/100 | w        1.12811 | loss        5.0118 | grad       -0.0025
      iter    75/100 | w        1.12811 | loss        5.0118 | grad       -0.0021
      iter    76/100 | w        1.12812 | loss        5.0118 | grad       -0.0018
      iter    77/100 | w        1.12812 | loss        5.0118 | grad       -0.0016
      iter    78/100 | w        1.12812 | loss        5.0118 | grad       -0.0014
      iter    79/100 | w        1.12812 | loss        5.0118 | grad       -0.0012
      iter    80/100 | w        1.12812 | loss        5.0118 | grad       -0.0010
      iter    81/100 | w        1.12812 | loss        5.0118 | grad       -0.0009
      iter    82/100 | w        1.12812 | loss        5.0118 | grad       -0.0007
      iter    83/100 | w        1.12812 | loss        5.0118 | grad       -0.0006
      iter    84/100 | w        1.12812 | loss        5.0118 | grad       -0.0006
      iter    85/100 | w        1.12812 | loss        5.0118 | grad       -0.0005
      iter    86/100 | w        1.12813 | loss        5.0118 | grad       -0.0004
      iter    87/100 | w        1.12813 | loss        5.0118 | grad       -0.0004
      iter    88/100 | w        1.12813 | loss        5.0118 | grad       -0.0003
      iter    89/100 | w        1.12813 | loss        5.0118 | grad       -0.0003
      iter    90/100 | w        1.12813 | loss        5.0118 | grad       -0.0002
      iter    91/100 | w        1.12813 | loss        5.0118 | grad       -0.0002
      iter    92/100 | w        1.12813 | loss        5.0118 | grad       -0.0002
      iter    93/100 | w        1.12813 | loss        5.0118 | grad       -0.0001
      iter    94/100 | w        1.12813 | loss        5.0118 | grad       -0.0001
      iter    95/100 | w        1.12813 | loss        5.0118 | grad       -0.0001
      iter    96/100 | w        1.12813 | loss        5.0118 | grad       -0.0001
      iter    97/100 | w        1.12813 | loss        5.0118 | grad       -0.0001
      iter    98/100 | w        1.12813 | loss        5.0118 | grad       -0.0001
      iter    99/100 | w        1.12813 | loss        5.0118 | grad       -0.0001


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

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

### Exercise 2d: What exactly is the gradient of the returned "optimal" value of w?

Use your `calc_grad` function to check the result. What is the gradient of the returned `wopt`?

Does this look totally converged? Can you find a $w$ value that would be even better?


```python
# TODO call calc_grad on the return value from above
```

## Part 3: Diagnostic plots for gradient descent

Let's look at some trace functions.

Whenever you run gradient descent, an *excellent* debugging strategy is the ability to plot the loss, the gradient magnitude, and the parameter of interest at every step of the algorithm.


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

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');

plt.xlim([0, 100]);
```

### Discussion 3a: What value do we expect the *loss* to converge to? Should it always be zero?

### Discussion 3b: What value do we expect the *gradient* to converge to? Should it always be zero?

# Part 4: Larger step sizes

## Try with larger step_size = 0.014


```python
wopt, info_dict = minimize_via_grad_descent(calc_loss, calc_grad, step_size=0.014);
```


```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 4a: What happens here? How is this step size different than in Part 3 above?

TODO discuss with your group

## Try with even larger step size 0.1


```python
wopt, info_dict = minimize_via_grad_descent(calc_loss, calc_grad, step_size=0.1, max_iters=25);
```

### Discussion 3b: What happens here with this even larger step size? Is it converging?

### Exercise 3c: What is the largest step size you can get to converge reasonably?


```python
# TODO try some other step sizes here
wopt, info_dict = minimize_via_grad_descent(calc_loss, calc_grad, step_size=0) # TODO fix step_size
```

# Part 5: Sensitivity to initial conditions



### Exercise 5a: Try to call the defined procedure with a different initial condition for $w$. What happens?

You could try $w = 5.0$ or something else.


```python
# TODO try some other initial condition for init_w
wopt2, info_dict2 = minimize_via_grad_descent(calc_loss, calc_grad, init_w=0, step_size=0.001, max_iters=10) # TODO fix step_size
```

### Exercise 5b: Try again with another initial value. 


```python
# TODO try some other initial condition for init_w
wopt3, info_dict3 = minimize_via_grad_descent(calc_loss, calc_grad, init_w=0, step_size=0.001, max_iters=10) # TODO fix
```

### Exercise 5c: Make a trace plot

Make a trace plot showing convergence from multiple different starting values for $w$. What do you notice?


```python
# TODO
```

# Part 6: Using scipy's built-in gradient optimization tools



```python
import scipy.optimize
```

Take a look at SciPy's built in minimization toolbox

<https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html#scipy.optimize.minimize>

We'll use "L-BFGS", a second-order method that uses the function and its gradient.

This is a "quasi-newton" method, which you can get an intuition for here:

https://en.wikipedia.org/wiki/Newton%27s_method_in_optimization


```python
result = scipy.optimize.minimize(calc_loss, 0.0, jac=calc_grad, method='L-BFGS-B')

# Returns an object with several fields, let's print the result to get an idea
print(result)
```


```python
print(str(result.message))
```


```python
best_w = result.x
print(best_w)
```


```python

```

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