# Gradient Descent--Huber Loss
### Author: Jayanth Raman

**Date: 10 Dec 2019** v2

HomeWork2 of Machine Learning Basics, Tokyo DataScience

# Huber Loss
Point-wise Huber loss is defined piecewise as:
$$
L_\gamma (e_i) = 
\begin{cases}
 e_i^2                   & \text{for } |e_i| \le \gamma, \quad \gamma > 0\\
 \gamma (2 |e_i| - \gamma), & \text{otherwise.}
\end{cases}
$$
where, $\gamma > 0$ is a threshold and it defines the point where the Huber loss changes from a quadratic to a linear function or vice versa.  And $e_i$ are the individual data points.

The overall Huber loss of $M$ data points is the average:
$$
L_\gamma (e) = \frac{1}{M} \sum_{i=1}^{M} L_\gamma (e_i)
$$
for the $M$ points of $e_i$.

When $e$ is the error, $e = \hat{y} - y$, then the point-wise loss is
$$
L_\gamma (y_i, \hat{y_i}) = 
\begin{cases}
 (\hat{y_i} - y_i)^2                   & \text{for } |\hat{y_i} - y_i| \le \gamma, \quad \gamma > 0\\
 \gamma (2 |\hat{y_i} - y_i| - \gamma), & \text{otherwise.}
\end{cases}
$$
where, $\hat{y}$ is the prediction of $y$.

And the overall Huber loss is the average
$$
L_\gamma (y, \hat{y}) = \frac{1}{M} \sum_{i=1}^{M} L_\gamma (y_i, \hat{y_i})
$$

[From [Wikipedia](https://en.wikipedia.org/wiki/Huber_loss)]

For a derivation of the linear term so that it is continuous with the quadratic term, see [Appendix A](#Appendix-A----Huber-Loss-Derivation).

## Gradients
A polynomial model is defined as
$$
\hat{y}_i = \sum_{d=0}^{N} w_d x^d_i
$$
where $\hat{y}_i$ is the predicted value of $y_i$ at $x_i$ and $w_d$ are the weights.

For the Huber-loss model, the first gradient component (in the $w_0$ direction) is
$$
\frac{\partial L_\gamma}{\partial w_0} = 
\begin{cases}
\frac{2}{M} \sum_{i=1}^{M} (\hat{y_i} - y_i)  & |\hat{y_i} - y_i| < \gamma \\
\frac{2\gamma}{M} \sum_{i=1}^{M} \text{sgn}(\hat{y_i} - y_i) & \text{otherwise}
\end{cases}
$$
where $\text{sgn}$ is the signum or sign function.

Let $e_i = \hat{y_i} - y_i$.  Notice that, when $e_i < \gamma$, then each term in
the summation is $\text{sgn}(e_i)\,|e_i|$.  And, when $e_i \ge \gamma$, then each term is
$\text{sgn}(e_i)\,\gamma$.  Hence, each term in the summation is $\text{sgn}(e_i)\,\min(|e_i|, \gamma)$.

Hence the first component of the gradient can be simplified to:
$$
\frac{\partial L_\gamma}{\partial w_0} =
\frac{2}{M} \sum_{i=1}^{M} \text{sgn}(e_i) \, \min \left( |e_i|, \gamma \right)
$$
where $e_i = \hat{y_i} - y_i$.  And the above expression can be computed as an average.

The other gradient components are defined recursively as
$$
\frac{\partial L_\gamma}{\partial w_d} = x \frac{\partial L_\gamma}{\partial w_{d - 1}} \quad d=1, 2, ..., N
$$

For a derivation of the above, please see [Appendix B](#Appendix-B----Gradient-for-Huber-Loss).

# Summary of Findings
Synthetic data was generated with a few outliers.  The data, $(x_i, y_i)$, is linear:
$$
y_i = w_0 + w_1 x_i
$$

Then, two models were fit:
 * One with an MSE-loss function.
 * Another with a Huber-loss function.  The parameter, $\gamma$, was chosen to be
   twice the standard deviation of the error in the MSE-loss model.

The models were compared using MSE loss (which may be an advantage to the MSE-loss model when
all of the data is considered).  MSE loss was computed for all of the data and also
just for the unaffected data (i.e. data without the outliers).

With all of the data in consideration, both models performed about the same with the MSE-loss model slightly outperforming the Huber-loss model.

But, when only the unaffected data was considered, the Huber-loss model outperformed the MSE-loss model.  The MSE-loss for the two models is shown in the table below.  The MSE loss for the Huber-loss model is about 38% lower than that for the MSE-loss model, when only the unaffected data is considered.  For more details, see [here](#Comparison-of-MSE-Loss).

The models were also trained on the same synthetic data, but with the outliers suppressed.  Both models were similar and both had an MSE-loss of approximately 0.1170.

| Model            | Loss (unaffected data) | Loss (all data) | Loss (no outliers) |
|------------------|-----------------------:|----------------:|-------------------:|
| MSE-loss Model   |                 0.2468 |          1.9935 |             0.1170 |
| Huber-loss Model |                 0.1516 |          2.0246 |             0.1171 |

The code below is easily applicable to quadratic data or any higher degree polynomials.

# Calculations and Observations

What follows is the source code for performing the calculations and the comparisons.


```python
# Imports
import bokeh
import bokeh.plotting as bkp
import functools
import numpy as np
import pandas as pd   # for summary statistics
#
bkp.reset_output()
bkp.output_notebook()
```



<div class="bk-root">
    <a href="https://bokeh.org" target="_blank" class="bk-logo bk-logo-small bk-logo-notebook"></a>
    <span id="1001">Loading BokehJS ...</span>
</div>




The function to compute $y$ as a polynomial of $x$ is given below.


```python
# Function to generate `y` given `x` and a set of weights.
def polynomial_function(weights, x):
    '''
    weights = [a0, a1, a2, ...]
    output = a0 + a1 * x + a2 * x ** 2 + ...
    '''
    if isinstance(x, list):
        x = np.array(x)
    out = weights[0]
    xn = x
    for w in weights[1:]:
        out = out + w * xn
        xn = xn * x
    return out

def test_polynomial_function():
    x1 = 3.1415
    tmp1 = 4 + 3.2 * x1 + 2.1 * x1 * x1 + 8.53 * x1 * x1 * x1
    assert polynomial_function([6.24, 4.14], x1) == 6.24 + 4.14 * x1
    assert polynomial_function([4, 3.2, 2.1, 8.53], x1) == tmp1
    assert all(polynomial_function([4, 3.2, 2.1, 8.53], [x1, x1]) == np.array([tmp1, tmp1]))
test_polynomial_function()
```

A helper function to generate training data is given below.  Outliers are added using samples from a binomial distribution.  Gaussian noise is added to all data.  The function outputs $x$, $y$ with no outlier noise, and outlier noise.  In order to generate $y$ with outliers, the outlier noise vector is added.


```python
def gen_training_data(wt_true, npoints=50, gauss_noise_sd=0.4, binomp={'p': 0.03, 'mag': 6.0}, seed=414243):
    '''Generate (simulate) training data'''
    np.random.seed(seed)
    xtrain = np.random.random(npoints)
    # e.g. 6 * binomial(1, 0.03, size=n) - 6 * binomial(1, 0.03, size=n)
    binomial = np.random.binomial
    ntrials = 1
    outlier_noise = binomial(ntrials, binomp['p'], size=npoints) - binomial(ntrials, binomp['p'], size=npoints)
    outlier_noise = outlier_noise * binomp['mag']
    ytrain = polynomial_function(wt_true, xtrain) + gauss_noise_sd * np.random.randn(npoints)
    return xtrain, ytrain, outlier_noise
```

We define two update functions, one using an MSE loss function, and another using a Huber loss function.
We will refer to the former as the MSE-loss model and the latter as the Huber-loss model.


```python
def update_weights_mse(weights, x, y, learning_rate):
    y_predicted = polynomial_function(weights, x)
    derivative_of_loss = 2 * (y_predicted - y)
    weights[0] -= learning_rate * derivative_of_loss.mean()
    m1 = derivative_of_loss.mean()
    for ii in range(1, len(weights)):
        derivative_of_loss = x * derivative_of_loss
        weights[ii] -= learning_rate * derivative_of_loss.mean()
    return None

def update_weights_huber(weights, x, y, learning_rate, gamma):
    y_predicted = polynomial_function(weights, x)
    err = y_predicted - y
    derivative_of_loss = 2.0 * np.sign(err) * np.minimum(gamma, np.abs(err))
    weights[0] -= learning_rate * derivative_of_loss.mean()
    m1 = derivative_of_loss.mean()
    for ii in range(1, len(weights)):
        derivative_of_loss = x * derivative_of_loss
        weights[ii] -= learning_rate * derivative_of_loss.mean()
    return None
```

We also define helper functions to compute the MSE loss and the Huber loss.


```python
def mse_loss(y, yhat):
    return np.mean((y - yhat) ** 2)

def huber_loss(y, yhat, gamma):
    errabs = np.abs(y - yhat)
    loss = np.where(errabs <= gamma, (errabs ** 2), gamma * (2.0 * errabs - gamma))
    return np.mean(loss)

# test
def test():
    wt_true = [3.0, 6.0]
    xtrain, ytrain, outlier_noise = gen_training_data(wt_true, npoints=70)
    _ytrue = polynomial_function(wt_true, xtrain)
    ytrain = ytrain + outlier_noise
    m1, m2 = mse_loss(_ytrue, ytrain), mse_loss(ytrain, _ytrue)
    print('MSE loss:', m1, m2)
    assert m1 == m2
    h1, h2 = huber_loss(_ytrue, ytrain, np.inf), huber_loss(ytrain, _ytrue, np.inf)
    print('Huber loss, gamma=np.inf:', h1, h2)
    assert h1 == h2
    assert h1 == m1
    gamma = np.std(ytrain)
    h1, h2 = huber_loss(_ytrue, ytrain, gamma), huber_loss(ytrain, _ytrue, gamma)
    print(f'Huber loss, gamma={gamma:.4f}: {h1}, {h2}')
    assert h1 == h2
test()
```

    MSE loss: 2.1804075712825233 2.1804075712825233
    Huber loss, gamma=np.inf: 2.1804075712825233 2.1804075712825233
    Huber loss, gamma=1.9832: 1.252087574809613, 1.252087574809613


A helper function to do the training is defined below.  It takes as parameters the $x$ (`xtrain`) and $y$ (`ytrain`) data along with the learning rate and update functions and the number of iterations to run.  An optional loss function, if provided, results in the loss being computed and returned.  This is useful in checking if the training has progressed properly.  The iteration interval when the loss is sampled can also be provided as `loss_iter`.  Optional initial weights can also be provided.  This is useful in continuing a previous training run.  If not provided, random weights are generated.


```python
# Train
def train(xtrain, ytrain, nweights, update_func, niter, learning_rate=0.02,
          loss_func=None, init_weights=None, loss_iter=20):
    # initialize weights
    weights = np.random.randn(nweights) if init_weights is None else init_weights
    loss_values = []
    for ii in range(niter):
        update_func(weights, xtrain, ytrain, learning_rate)
        if loss_func and ii % loss_iter == (loss_iter - 1):
            loss_values.append(loss_func(polynomial_function(weights, xtrain), ytrain))
    return weights, loss_values
```

Training data is generated below using weights that define a linear relationship.  The number of outliers that were generated is shown below.


```python
wt_true = [3.0, 6.0]
xtrain, ytrain, outlier_noise = gen_training_data(wt_true, npoints=70)
print('Number of outliers:', sum(np.where(outlier_noise, 1, 0)))
```

    Number of outliers: 4


Let's look at the summary statistics of the generated data, with and without outliers.


```python
display(pd.DataFrame({'y (no outliers)': ytrain}).describe())
print()
display(pd.DataFrame({'y (with outliers)': ytrain + outlier_noise}).describe())
```


<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>y (no outliers)</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>count</th>
      <td>70.000000</td>
    </tr>
    <tr>
      <th>mean</th>
      <td>5.993224</td>
    </tr>
    <tr>
      <th>std</th>
      <td>1.806728</td>
    </tr>
    <tr>
      <th>min</th>
      <td>2.545869</td>
    </tr>
    <tr>
      <th>25%</th>
      <td>4.440494</td>
    </tr>
    <tr>
      <th>50%</th>
      <td>6.017799</td>
    </tr>
    <tr>
      <th>75%</th>
      <td>7.658137</td>
    </tr>
    <tr>
      <th>max</th>
      <td>8.887732</td>
    </tr>
  </tbody>
</table>
</div>


    



<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>y (with outliers)</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>count</th>
      <td>70.000000</td>
    </tr>
    <tr>
      <th>mean</th>
      <td>5.993224</td>
    </tr>
    <tr>
      <th>std</th>
      <td>1.997503</td>
    </tr>
    <tr>
      <th>min</th>
      <td>0.680286</td>
    </tr>
    <tr>
      <th>25%</th>
      <td>4.440494</td>
    </tr>
    <tr>
      <th>50%</th>
      <td>6.017799</td>
    </tr>
    <tr>
      <th>75%</th>
      <td>7.735788</td>
    </tr>
    <tr>
      <th>max</th>
      <td>9.859972</td>
    </tr>
  </tbody>
</table>
</div>


As expected, the standard deviation of the data with outliers is higher than the one without.  The minimum and maximum values are also wider.  Otherwise the mean and median are unchanged since the positive and negative outliers are balanced in the data that was produced.

We plot the data below.  The outliers are plotted in a different color.  The dashed line corresponds to the true weights that generated the data.


```python
p = bkp.figure(title='Training Data With Outliers', y_range=[-1, 11], width=600, height=400)
xin, xout = xtrain[outlier_noise == 0], xtrain[outlier_noise != 0]
yin, yout = ytrain[outlier_noise == 0], (ytrain + outlier_noise)[outlier_noise != 0]
p.circle(xin, yin, radius=0.01, legend_label='unaffected data')
p.circle(xout, yout, radius=0.01, color='darkorange', legend_label='outliers')
line = bokeh.models.Slope(gradient=wt_true[1], y_intercept=wt_true[0], 
                          line_color='gray', line_dash='dashed', line_width=3.5)
p.add_layout(line)
p.legend.location = 'bottom_left'
p.legend.border_line_alpha = 0.4
p.legend.border_line_color = 'black'
bkp.show(p)
```








<div class="bk-root" id="3f634038-ae32-48fa-bb89-466412f6296f" data-root-id="1002"></div>





## Predictions Without Outliers
First, we compare predictions using the MSE-loss and Huber-loss models using data without outliers i.e. we do not yet add the noise that will move some of the points and make them outliers.

### MSE-loss Model
First, we fit the MSE-loss model.


```python
# No outliers - Predict using MSE Loss
np.random.seed(424242)
weights_mse, loss_values_mse = train(xtrain, ytrain, nweights=2, update_func=update_weights_mse,
                                     loss_func=mse_loss, niter=250 * 20)
```


```python
print('loss values (last few):', loss_values_mse[-4:])
print('True weights:', wt_true)
print('Predicted weights:', weights_mse)
yhat = polynomial_function(weights_mse, xtrain)
print('MSE loss:', mse_loss(yhat, ytrain))
```

    loss values (last few): [0.11702782079848476, 0.11702782079843213, 0.11702782079838533, 0.11702782079834366]
    True weights: [3.0, 6.0]
    Predicted weights: [3.07659014 5.79207763]
    MSE loss: 0.11702782079834366


Notice that the loss values have stabilized.  Later, we also plot the loss against the iteration steps to drive home this point.

The predicted weights are close to the true weights.  And the MSE loss of the predictions is small compared to either the average value or the median value of $y$.

The error distribution is plotted below.  The two dashed lines are plotted at $\pm \gamma$.  The axis limits are set so that we can later directly compare it to the error generated with the outliers added.

In the Huber-loss model, $\gamma$ will be set to twice the standard deviation of the error.


```python
# Plot the errors
p = bkp.figure(title='Error', x_range=[-7, 7], y_range=[0.5, 0.8], width=600, height=200,
               x_axis_label='error=(yhat - y)')
err = polynomial_function(weights_mse, xtrain) - ytrain
p.circle(err, 0.6, fill_color='blue')
gamma = 2 * np.std(err)
p.line([-gamma, -gamma], [0, 1], line_dash='dashed', line_width=2)
p.line([gamma, gamma], [0, 1], line_dash='dashed', line_width=2)
p.yaxis.visible = False
bkp.show(p)
```








<div class="bk-root" id="252c45ef-9cd5-42f0-8f94-22792e3e93bc" data-root-id="1134"></div>





The error data is bunched up near the origin (again, the x-axis scale is set to account for the outliers which come later).  The lines at $\pm \gamma$ are at the data edges and there are really no outliers to be separated here.

### Huber-loss Model
Next we predict using the Huber-loss model.  We use twice the standard deviation of the error from the MSE-loss model as the $\gamma$ value.


```python
# No outliers - Predict using Huber Loss
np.random.seed(424242)
# gamma = np.std(ytrain + outlier_noise)
gamma = 2 * np.std(polynomial_function(weights_mse, xtrain) - ytrain)
print('gamma:', gamma)
update_func = functools.partial(update_weights_huber, gamma=gamma)
loss_func = functools.partial(huber_loss, gamma=gamma)
weights_huber, loss_values_huber = train(xtrain, ytrain, nweights=2, update_func=update_func,
                                         loss_func=loss_func, niter=250 * 20)
```

    gamma: 0.6841865850726339



```python
print('gamma:', gamma)
print('loss values (last few):', loss_values_huber[-4:])
print('True weights:', wt_true)
print('Predicted weights:', weights_huber)
yhat = polynomial_function(weights_huber, xtrain)
print('MSE loss:', mse_loss(yhat, ytrain))
```

    gamma: 0.6841865850726339
    loss values (last few): [0.11154368969877192, 0.11154368969833425, 0.11154368969794404, 0.11154368969759612]
    True weights: [3.0, 6.0]
    Predicted weights: [3.07095171 5.78240498]
    MSE loss: 0.11714691181634618


Again, the loss values have stabilized.

The predicted weights are close to the true weights and are *practically the same as that from the MSE-loss model.*

Next, we plot the predictions of both models against the data.


```python
p = bkp.figure(title='Predictions -- Data Without Outliers', y_range=[-1, 11], width=600, height=400)
p.circle(xtrain, ytrain, radius=0.01)
xtmp = np.linspace(0, 1)
ytrue = polynomial_function(wt_true, xtmp)
ymse = polynomial_function(weights_mse, xtmp)
yhuber = polynomial_function(weights_huber, xtmp)
p.line(xtmp, ytrue, color='gray', line_width=3, line_dash='dashed', legend_label='true weight')
p.line(xtmp, ymse, color='blue', line_width=3, legend_label='MSE')
p.line(xtmp, yhuber, color='darkorange', line_width=3, legend_label='Huber')
p.legend.location = 'bottom_right'
bkp.show(p)
```








<div class="bk-root" id="42e74734-4fcc-4423-8de3-3fad3f2c02a8" data-root-id="1277"></div>





From the above figure, the three predictions are practically indistinguishable.

Both the models yielded the same **MSE loss of 0.117**.

The true weights and the weights from the predictions are below:


```python
print(f'True weights      : {wt_true}')
print(f'MSE-loss weights  : {weights_mse}')
print(f'Huber-loss weights: {weights_huber}')
```

    True weights      : [3.0, 6.0]
    MSE-loss weights  : [3.07659014 5.79207763]
    Huber-loss weights: [3.07095171 5.78240498]


### Training Loss
For completeness, we plot the training loss against the training iterations for both MSE and Huber losses.  The training loss decreases with the number of steps as expected.  The Huber-loss model converges a bit more slowly than the MSE-loss model.


```python
p = bkp.figure(title='Training Loss vs Iterations (no outliers)', width=600, height=400)
p.line(np.arange(len(loss_values_mse)) * 20, loss_values_mse, line_width=2, legend_label='MSE')
p.line(np.arange(len(loss_values_huber)) * 20, loss_values_huber, line_width=2,
       color='darkorange', legend_label='Huber')
p.xaxis.axis_label = 'Training Steps'
p.yaxis.axis_label = 'Loss'
bkp.show(p)
```








<div class="bk-root" id="720acb7d-90a3-4151-aeb7-30da0ce4c6bb" data-root-id="1485"></div>





## Predictions With Outliers
We now add the outlier noise, thus creating the outliers, and then train the linear model using both the MSE-loss and Huber-loss models.

### MSE-loss Model
We train the MSE-loss model.


```python
# With outliers - Predict using MSE Loss
np.random.seed(424242)
weights_mse_ol, loss_values_mse_ol = train(xtrain, ytrain + outlier_noise, nweights=2,
                                           update_func=update_weights_mse,
                                           loss_func=mse_loss, niter=250 * 20)
```


```python
print('loss values (last few):', loss_values_mse_ol[-4:])
print('True weights:', wt_true)
print('Predicted weights:', weights_mse_ol)
yhat = polynomial_function(weights_mse_ol, xtrain)
print('MSE loss (all data):', mse_loss(yhat, ytrain + outlier_noise))
print('MSE loss (unaffected data):', mse_loss(yhat[outlier_noise == 0], ytrain[outlier_noise == 0]))
```

    loss values (last few): [1.9935106591383374, 1.993510659138311, 1.993510659138289, 1.9935106591382687]
    True weights: [3.0, 6.0]
    Predicted weights: [3.68645222 4.58096621]
    MSE loss (all data): 1.9935106591382687
    MSE loss (unaffected data): 0.24677995511981501


Notice that
 * the loss values have stabilized.
 * the predicted weights have diverged from the true weights, compared to the model build on the data without outliers.
 * the MSE loss for all of the data has increased to almost 2.0 from 0.117.
 * the MSE loss for the unaffected data has increased to 0.2468.

#### Error characteristics
Let's look at the error characteristics.  Ideally, one may want to base the value of the Huber-loss parameter, $\gamma$, based on the standard deviation of the error.


```python
yhat = polynomial_function(weights_mse_ol, xtrain)
err_mse_ol = yhat - ytrain - outlier_noise
pd.DataFrame({'error': err_mse_ol}).describe()
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>error</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>count</th>
      <td>7.000000e+01</td>
    </tr>
    <tr>
      <th>mean</th>
      <td>5.191414e-08</td>
    </tr>
    <tr>
      <th>std</th>
      <td>1.422112e+00</td>
    </tr>
    <tr>
      <th>min</th>
      <td>-5.657756e+00</td>
    </tr>
    <tr>
      <th>25%</th>
      <td>-2.994260e-01</td>
    </tr>
    <tr>
      <th>50%</th>
      <td>-1.110702e-02</td>
    </tr>
    <tr>
      <th>75%</th>
      <td>3.095180e-01</td>
    </tr>
    <tr>
      <th>max</th>
      <td>5.935515e+00</td>
    </tr>
  </tbody>
</table>
</div>



#### Choice of $\gamma$
For the Huber-loss function, we will choose $\gamma$ to be twice the standard deviation of the error.

The plot below shows the distribution of the error along with dashed lines at plus or minus one standard deviation of the error.  Observe that the lines at $\pm \gamma$ clearly separate the outliers from the unaffected data.  Compare this plot to the earlier plot of the error for when there were no outliers in the data.


```python
# Plot the errors
p = bkp.figure(title='Error', x_range=[-7, 7], y_range=[0.5, 1.0], width=600, height=200,
               x_axis_label='error=(yhat - y)')
p.circle(err_mse_ol[outlier_noise == 0], 0.6, fill_color='blue', legend_label='unaffected data')
p.circle(err_mse_ol[outlier_noise != 0], 0.6 + 0.01 * np.random.randn(4), color='darkorange',
         radius=0.1, legend_label='outlier')
estd = 2 * np.std(err_mse_ol)
p.line([-estd, -estd], [0, 1], line_dash='dashed', line_width=2)
p.line([estd, estd], [0, 1], line_dash='dashed', line_width=2)
p.legend.location = 'top_right'
p.yaxis.visible = False
bkp.show(p)
```








<div class="bk-root" id="3b7187fd-2523-47b8-9ac1-77a6b7221ed8" data-root-id="1652"></div>





### Huber-loss model
Next, we fit the Huber-loss model.


```python
# With outliers - Predict using Huber Loss
np.random.seed(424242)
gamma = 2 * np.std(err_mse_ol)
update_func = functools.partial(update_weights_huber, gamma=gamma)
loss_func = functools.partial(huber_loss, gamma=gamma)
weights_huber_ol, loss_values_huber_ol = train(xtrain, ytrain + outlier_noise, nweights=2, update_func=update_func,
                                               loss_func=loss_func, niter=250 * 20)
```


```python
print('gamma:', gamma)
print('loss values (last few):', loss_values_huber_ol[-4:])
print('True weights:', wt_true)
print('Predicted weights:', weights_huber_ol)
yhat = polynomial_function(weights_huber_ol, xtrain)
print('MSE loss (all data):', mse_loss(yhat, ytrain + outlier_noise))
yhat = polynomial_function(weights_huber_ol, xtrain[outlier_noise == 0])
print('MSE loss (unaffected data):', mse_loss(yhat, ytrain[outlier_noise == 0]))
```

    gamma: 2.8238347395966823
    loss values (last few): [1.5358121300535852, 1.5358121300532634, 1.5358121300529748, 1.5358121300527157]
    True weights: [3.0, 6.0]
    Predicted weights: [3.38589048 5.1598805 ]
    MSE loss (all data): 2.02456696435184
    MSE loss (unaffected data): 0.15162498018067325


Notice that
 * the loss values have stabilized.
 * the predicted weights have diverged from the true weights, but not as much as the MSE-loss model.
 * the MSE loss for all of the data has increased to over 2.0 from 0.117.
 * the MSE loss for the unaffected data has increased to 0.1516.  But this is lower than the loss for the MSE-loss model.


```python
p = bkp.figure(title='Predictions -- Huber-loss -- Data With Outliers', y_range=[-1, 11], width=600, height=400)
p.circle(xtrain, ytrain + outlier_noise, radius=0.01)
xtmp = np.linspace(0, 1)
ytrue = polynomial_function(wt_true, xtmp)
ymse = polynomial_function(weights_mse_ol, xtmp)
yhuber = polynomial_function(weights_huber_ol, xtmp)
p.line(xtmp, ytrue, color='gray', line_width=3, line_dash='dashed', legend_label='true weight')
p.line(xtmp, ymse, color='blue', line_width=3, legend_label='MSE')
p.line(xtmp, yhuber, color='darkorange', line_width=3, legend_label='Huber')
p.legend.location = 'bottom_right'
bkp.show(p)
```








<div class="bk-root" id="817b475b-f17b-40d5-83a6-774623edd370" data-root-id="1867"></div>





Observe that the Huber-loss model line is closer to the true line than the MSE-loss model line.

### Training Loss
For completeness, we plot the training loss against the training iterations for both MSE and Huber losses.  The training loss decreases with the number of steps as expected.  We can also observe that the training loss for the MSE-loss model is higher than that of the Huber-loss model.


```python
p = bkp.figure(title='Training Loss vs Iterations (data with outliers)', width=600, height=400)
p.line(np.arange(len(loss_values_mse_ol)) * 20, loss_values_mse_ol, line_width=2, legend_label='MSE')
p.line(np.arange(len(loss_values_huber_ol)) * 20, loss_values_huber_ol, line_width=2,
       color='darkorange', legend_label='Huber')
p.xaxis.axis_label = 'Training Steps'
p.yaxis.axis_label = 'Loss'
bkp.show(p)
```








<div class="bk-root" id="1414fbad-b461-4b40-bc29-40fbb79abe91" data-root-id="2113"></div>





## Comparison of Error Characteristics

We compare the characteristics of the error for the two models for all the data (including outliers) and for just the unaffected data.


```python
# Compute error for both models, with and without outliers, and display their statistics
y_mse_ol = polynomial_function(weights_mse_ol, xtrain)
err_mse_ol = (y_mse_ol - ytrain - outlier_noise)
#
y_huber_ol = polynomial_function(weights_huber_ol, xtrain)
err_huber_ol = (y_huber_ol - ytrain - outlier_noise)
#
tmp = pd.DataFrame({
    ('All Data', 'MSE-loss Model'): err_mse_ol,
    ('All Data', 'Huber-loss Model'): err_huber_ol,
    ('Unaffected Data', 'MSE-loss Model'): np.concatenate((err_mse_ol[outlier_noise == 0], np.nan * np.ones(4))),
    ('Unaffected Data', 'Huber-loss Model'): np.concatenate((err_huber_ol[outlier_noise == 0], np.nan * np.ones(4))),
}).describe()
display(tmp)
```


<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead tr th {
        text-align: left;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr>
      <th></th>
      <th colspan="2" halign="left">All Data</th>
      <th colspan="2" halign="left">Unaffected Data</th>
    </tr>
    <tr>
      <th></th>
      <th>MSE-loss Model</th>
      <th>Huber-loss Model</th>
      <th>MSE-loss Model</th>
      <th>Huber-loss Model</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>count</th>
      <td>7.000000e+01</td>
      <td>70.000000</td>
      <td>66.000000</td>
      <td>6.600000e+01</td>
    </tr>
    <tr>
      <th>mean</th>
      <td>5.191414e-08</td>
      <td>-0.009046</td>
      <td>0.006265</td>
      <td>1.993456e-07</td>
    </tr>
    <tr>
      <th>std</th>
      <td>1.422112e+00</td>
      <td>1.433117</td>
      <td>0.500536</td>
      <td>3.923744e-01</td>
    </tr>
    <tr>
      <th>min</th>
      <td>-5.657756e+00</td>
      <td>-5.893139</td>
      <td>-1.359935</td>
      <td>-1.332695e+00</td>
    </tr>
    <tr>
      <th>25%</th>
      <td>-2.994260e-01</td>
      <td>-0.216800</td>
      <td>-0.280476</td>
      <td>-2.026042e-01</td>
    </tr>
    <tr>
      <th>50%</th>
      <td>-1.110702e-02</td>
      <td>0.004658</td>
      <td>-0.011107</td>
      <td>4.657694e-03</td>
    </tr>
    <tr>
      <th>75%</th>
      <td>3.095180e-01</td>
      <td>0.277028</td>
      <td>0.306297</td>
      <td>2.218965e-01</td>
    </tr>
    <tr>
      <th>max</th>
      <td>5.935515e+00</td>
      <td>6.005147</td>
      <td>1.203468</td>
      <td>9.108528e-01</td>
    </tr>
  </tbody>
</table>
</div>


**Observations**
 * "All Data" refers to all the data including outliers, and "Unaffected Data" refers to
   data excluding the outliers.
 * *Mean*: Notice that the mean error is close to zero for both models, with and without outliers.
   This is probably because we have two outliers on either side of the unaffected data.
 * *Standard deviation--All Data*: The standard deviation for all the data is slightly higher for the Huber-loss
   model than for the MSE-loss model.  By this measure it would seem that both the models are similar,
   with the MSE-model being slightly better.
 * *Standard deviation--Unaffected Data*: But, when we look at the standard deviation for the
   unaffected data, then the Huber-loss model really shines.  It's error `std` is 0.037 compared to 0.5
   for the MSE-loss and is an order of magnitude better.
 * *Min-Max All Data*: Again, the minimum and maximum error values are slightly worse for the
   Huber-loss model--implying that it underperforms for the outliers themselves.  This makes sense
   given the position of the outliers in our data.  The outliers cause the Huber-loss line to be closer
   to the actual line which causes the errors at the outliers to be larger.
 * *Min-Max Unaffected Data*: For the unaffected data, the roles are reversed--the Huber-loss model
   slightly outperforms the MSE-loss model.  For example, the maximum error is 0.82 compared to 1.2.
   Again, for our data, this makes sense since the Huber-loss line is closer to the actual line.

Plots of the errors for the two models are shown below.  Also, the outliers are shown in
a different color than the unaffected data.


```python
# Plot the errors
p = bkp.figure(title='Predictions -- Data Without Outliers', x_range=[-7, 7], y_range=[0, 1], width=600, height=400)
p.circle(err_mse_ol[outlier_noise == 0], 0.6, fill_color='blue', legend_label='unaffected data')
p.circle(err_mse_ol[outlier_noise != 0], 0.6 + 0.01 * np.random.randn(4), color='darkorange',
         radius=0.1, legend_label='outlier')
p.add_layout(bokeh.models.Label(x=0, y=0.7, text='MSE-Loss Model', text_align='center', text_baseline='top'))
p.circle(err_huber_ol[outlier_noise == 0], 0.3, fill_color='blue', legend_label='unaffected data')
p.circle(err_huber_ol[outlier_noise != 0], 0.3 + 0.01 * np.random.randn(4), color='darkorange',
         radius=0.1, legend_label='outlier')
p.add_layout(bokeh.models.Label(x=0, y=0.2, text='Huber-Loss Model', text_align='center'))
bkp.show(p)
```








<div class="bk-root" id="74ac1640-5f8c-4f99-beb0-187a9981177d" data-root-id="2318"></div>





## Comparison of MSE Loss
The MSE loss using the weights of the following models is shown in the table below:
 1. Model using MSE loss, and
 2. Model using Huber loss

| Model            | Loss (unaffected data) | Loss (all data) | Loss (no outliers) |
|------------------|-----------------------:|----------------:|-------------------:|
| MSE-loss Model   |                 0.2468 |          1.9935 |             0.1170 |
| Huber-loss Model |                 0.1516 |          2.0246 |             0.1171 |

Legend:
 * **Loss (unaffected data)**: The fitted data contains outliers, but the outliers
   were left out when computing the MSE loss.
 * **Loss (all data)**: The fitted data contains outliers, and the MSE loss was
   computed on all the data including the outliers.
 * **Loss (no outliers)**: The fitted data does not contain outliers i.e. the
   outlier noise was not added to the data.  The MSE loss is for all of the data.

Observations:
 * When there are no outliers, there is practically no difference between the two models.
 * When the model is fit on data with outliers, and when the MSE loss is computed on the
   entire dataset, the MSE-loss model slightly outperforms the Huber-loss model.
 * When the model is fit on data with outliers, but the MSE loss is computed on the data
   excluding the outliers, then the Huber-loss model significantly outperforms the
   MSE-loss model.
 * Since the comparison is done using MSE-loss it is but natural that the MSE-loss model
   outperforms the Huber-loss model when all of the data is considered.

# Appendix A -- Huber Loss Derivation

Define the error, $e_i$, as
$$
e_i = \hat{y}_i - y_i
$$

The total L2 loss for all the $M$ datapoints is
$$
L_2 = \frac{1}{M} \sum_{i=1}^{M} e_i^2
$$

The quadratic L2 loss for each data point is
$$
L_2^{(i)} = e_i^2
$$

For the Huber loss function, we wish to convert from a quadratic to a linear function.
Hence, we define the L1 loss as
$$
L_1^{(i)} = c_1 + c_2 |e_i|
$$
for some constants $c_1$ and $c_2$.

For clarity we temporarily drop the $i$ subscript.  And we require that the Huber loss function
be continuous at the change point $\gamma > 0$ and also that the first derivative be continuous.
$$
\begin{align}
L_2 & = L_1 \\
\frac{\partial L_2}{\partial e} & = \frac{\partial L_1}{\partial e}
\end{align}
$$

Substituting for $L_1$ and $L_2$ at $e=\gamma$, we get
$$
\begin{align}
\gamma^2 & = c_1 + c_2 \gamma \\
2 \gamma & = c_2
\end{align}
$$

Which yields, $c_2 = 2 \gamma$ and $c_1 = -\gamma^2$.

Hence, the linear part of the Huber loss is
$$
\begin{align}
L_1^{(i)} & = c_1 + c_2 |e_i| \\
  & = -\gamma^2 + 2 \gamma |e_i| \\
  & = \gamma (2|e_i| - \gamma)
\end{align}
$$

And the point-wise Huber loss function is then
$$
L_\gamma (e_i) = 
\begin{cases}
 e_i^2                   & \text{for } |e_i| \le \gamma, \quad \gamma > 0\\
 \gamma (2 |e_i| - \gamma), & \text{otherwise.}
\end{cases}
$$


# Appendix B -- Gradient for Huber Loss

Define the point-wise error, $e_i$, as $e_i = \hat{y_i} - y_i$.

Where, for an $N$-degree polynomial model, the prediction, $\hat{y_i}$ for input $x_i$, is given by
$$
\hat{y_i} = w_0 + w_1 x + w_2 x^2 + ... + w_N x^N = \sum_{d=0}^{N} w_i x^i
$$

The point-wise Huber loss is defined as
$$
L_i = L_\gamma (e_i) = 
\begin{cases}
 e_i^2                   & \text{for } |e_i| \le \gamma, \quad \gamma > 0\\
 \gamma (2 |e_i| - \gamma), & \text{otherwise.}
\end{cases}
$$

To derive the gradient, we apply the following chain rule
$$
\frac{\partial L_i}{\partial w_d} = \frac{\partial L_i}{\partial e_i} \frac{\partial e_i}{\partial w_d}
$$

And,
$$
\frac{\partial e_i}{\partial w_d} = x^i
$$

Let's derive the gradients for the linear and quadratic terms separately.

For the quadratic term,
$$
\frac{\partial L_i}{e_i} = \frac{\partial e_i^2}{e_i} = 2 e_i \quad |e_i| \le \gamma
$$

For the linear term,
$$
\frac{\partial L_i}{e_i} = \frac{\partial \gamma (2 |e_i| - \gamma)}{e_i}
= 2 \gamma \text{sgn}(e_i) \quad |e_i| > \gamma
$$

The above two can be combined into one expression
$$
\frac{\partial L_i}{\partial e_i} = 2 \text{sgn}(e_i) \, \min(|e_i|, \gamma)
$$

By the chain rule and substituting for $\frac{\partial e_i}{w_i}$,
$$
\frac{\partial L_i}{\partial w_d} = 2 \text{sgn}(e_i) \, \min(|e_i|, \gamma) x^i
$$

Which can be defined recursively as follows.  This definition is suitable to
compute the gradients sequentially along each of the weights, $w_d$.
$$
\begin{align}
\frac{\partial L_i}{\partial w_0} & = 2 \, \text{sgn}(e_i) \, \min(|e_i|, \gamma) \\
\frac{\partial L_i}{\partial w_d} & = x \frac{\partial L_i}{\partial w_{d - 1}} \quad d = 1, 2, ..., N
\end{align}
$$



# Appendix C -- Variation of $\gamma$


```python
# different gamma
np.random.seed(424242)
gamma = 6.0
update_func = functools.partial(update_weights_huber, gamma=gamma)
loss_func = functools.partial(huber_loss, gamma=gamma)
weights, loss_values = train(xtrain, ytrain + outlier_noise, nweights=2, update_func=update_func,
                             loss_func=loss_func, niter=250 * 20)
print(loss_values[-10:])
print(wt_true)
print(weights)
```

    [1.9935106591385687, 1.993510659138518, 1.993510659138473, 1.9935106591384324, 1.993510659138397, 1.993510659138365, 1.9935106591383362, 1.993510659138311, 1.9935106591382883, 1.9935106591382683]
    [3.0, 6.0]
    [3.68645222 4.58096621]



```python

```

88be54cdad032f524185fb748a752b6cb3b9e316
ipynb
Jupyter Notebook
src/basics-ml/hw2_huber_loss.ipynb
jraman/tds-submissions
27fe4894c77dc4bc3fda7c53d4e530ecfb6fbca9
Apache-2.0
src/basics-ml/hw2_huber_loss.ipynb
jraman/tds-submissions
27fe4894c77dc4bc3fda7c53d4e530ecfb6fbca9
Apache-2.0
src/basics-ml/hw2_huber_loss.ipynb
jraman/tds-submissions
27fe4894c77dc4bc3fda7c53d4e530ecfb6fbca9
Apache-2.0
Qwen/Qwen-72B
1. YES
2. YES