Least Squares with {L}^{1} / {L}^{2} Regularization

Least Squares with \({L}^{1}\) / \({L}^{2}\) Regularization#

Notebook by:

Royi Avital RoyiAvital@fixelalgorithms.com

Revision History#

Version	Date	User	Content / Changes
1.0.000	24/03/2024	Royi Avital	First version

# Import Packages

# General Tools
import numpy as np
import scipy as sp
import pandas as pd

# Machine Learning
from sklearn.linear_model import Lasso, Ridge
from sklearn.linear_model import lars_path, lasso_path
from sklearn.metrics import r2_score
from sklearn.model_selection import train_test_split
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import PolynomialFeatures

# Miscellaneous
import math
import os
from platform import python_version
import random
import timeit

# Typing
from typing import Callable, Dict, List, Optional, Set, Tuple, Union

# Visualization
import matplotlib as mpl
import matplotlib.pyplot as plt
import seaborn as sns

# Jupyter
from IPython import get_ipython
from IPython.display import Image
from IPython.display import display
from ipywidgets import Dropdown, FloatSlider, interact, IntSlider, Layout, SelectionSlider
from ipywidgets import interact

Notations#

(?) Question to answer interactively.
(!) Simple task to add code for the notebook.
(@) Optional / Extra self practice.
(#) Note / Useful resource / Food for thought.

Code Notations:

someVar    = 2; #<! Notation for a variable
vVector    = np.random.rand(4) #<! Notation for 1D array
mMatrix    = np.random.rand(4, 3) #<! Notation for 2D array
tTensor    = np.random.rand(4, 3, 2, 3) #<! Notation for nD array (Tensor)
tuTuple    = (1, 2, 3) #<! Notation for a tuple
lList      = [1, 2, 3] #<! Notation for a list
dDict      = {1: 3, 2: 2, 3: 1} #<! Notation for a dictionary
oObj       = MyClass() #<! Notation for an object
dfData     = pd.DataFrame() #<! Notation for a data frame
dsData     = pd.Series() #<! Notation for a series
hObj       = plt.Axes() #<! Notation for an object / handler / function handler

Code Exercise#

Single line fill

vallToFill = ???

Multi Line to Fill (At least one)

# You need to start writing
????

Section to Fill

#===========================Fill This===========================#
# 1. Explanation about what to do.
# !! Remarks to follow / take under consideration.
mX = ???

???
#===============================================================#

# Configuration
# %matplotlib inline

seedNum = 512
np.random.seed(seedNum)
random.seed(seedNum)

# Matplotlib default color palette
lMatPltLibclr = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd', '#8c564b', '#e377c2', '#7f7f7f', '#bcbd22', '#17becf']
# sns.set_theme() #>! Apply SeaBorn theme

runInGoogleColab = 'google.colab' in str(get_ipython())

# Constants

FIG_SIZE_DEF    = (8, 8)
ELM_SIZE_DEF    = 50
CLASS_COLOR     = ('b', 'r')
EDGE_COLOR      = 'k'
MARKER_SIZE_DEF = 10
LINE_WIDTH_DEF  = 2

# Courses Packages
import sys
sys.path.append('../')
sys.path.append('../../')
sys.path.append('../../../')
from utils.DataVisualization import PlotRegressionData

# General Auxiliary Functions

def PlotPolyFit( vX: np.ndarray, vY: np.ndarray, vP: Optional[np.ndarray] = None, P: int = 1, numGridPts: int = 1001, 
                hA: Optional[plt.Axes] = None, figSize: Tuple[int, int] = FIG_SIZE_DEF, markerSize: int = MARKER_SIZE_DEF, 
                lineWidth: int = LINE_WIDTH_DEF, axisTitle: str = None ) -> None:

    if hA is None:
        hF, hA = plt.subplots(1, 2, figsize = figSize)
    else:
        hF = hA[0].get_figure()

    numSamples = len(vY)

    # Polyfit
    vW    = np.polyfit(vX, vY, P)
    
    # MSE
    vHatY = np.polyval(vW, vX)
    MSE   = (np.linalg.norm(vY - vHatY) ** 2) / numSamples
    
    # Plot
    xx  = np.linspace(np.floor(np.min(vX)), np.ceil(np.max(vX)), numGridPts)
    yy  = np.polyval(vW, xx)

    hA[0].plot(vX, vY, '.r', ms = 10, label = '$y_i$')
    hA[0].plot(xx, yy, 'b',  lw = 2,  label = '$\hat{f}(x)$')
    hA[0].set_title (f'$P = {P}$\nMSE = {MSE}')
    hA[0].set_xlabel('$x$')
    # hA[0].axis(lAxis)
    hA[0].grid()
    hA[0].legend()
    
    hA[1].stem(vW[::-1], label = 'Estimated')
    if vP is not None:
        hA[1].stem(vP[::-1], linefmt = 'C1:', markerfmt = 'D', label = 'Ground Truth')
    numTicks = len(vW) if vP is None else max(len(vW), len(vP))
    hA[1].set_xticks(range(numTicks))
    hA[1].set_title('Coefficients')
    hA[1].set_xlabel('$w$')
    hA[1].legend()

Linear Regression with LASSO (\({L}_{1}\)) / Ridge (\({L}_{2}\)) Regularization#

Regularization is a simple and effective way to balance the adaptation of the model to the data.
It allows us to use complex model yet to tune it to prevent overfit.

The models of optimization:

LASSO: \(\arg \min_{\boldsymbol{w}} \frac{1}{2} {\left\| X \boldsymbol{w} - \boldsymbol{y} \right\|}_{2}^{2} + \lambda {\left\| \boldsymbol{w} \right\|}_{1}\).
Promotes sparsity, modeled by a Laplace prior of the coefficients of \(\boldsymbol{w}\).
Ridge: \(\arg \min_{\boldsymbol{w}} \frac{1}{2} {\left\| X \boldsymbol{w} - \boldsymbol{y} \right\|}_{2}^{2} + \lambda {\left\| \boldsymbol{w} \right\|}_{2}^{2}\).
Promotes damped coefficients, modeled by a Gaussian prior of the coefficients of \(\boldsymbol{w}\).

(#) The LASSO model is implemented by the Lasso class in SciKit Learn.
(#) The Ridge model is implemented by the Ridge class in SciKit Learn.
(#) Both Lasso and Ridge (As well as LinearRegression) have the positive parameter to enforce non negative values for \(\boldsymbol{w}\).
(#) The SciKit Learn classes use \(\alpha\) instead of \(\lambda\) in their optimization function.

In this notebook we’ll show the effect of the regularization on the fitness level of the model.

# Parameters

# Data Generation
numSamples  = 30
noiseStd    = 17.5

vP = np.array([0.5, 0.25, 2, 5])

# Model
lPolyDeg = [3, 5, 7, 9]
lλ       = list(np.linspace(0, 200, 1001)) #<! Pay attention that for λ = 0 it is better to use LinearRegression model (Numerical reasons)
lλ       = lλ[1:] #<! Prevent issues with λ = 0
numSamplesTrain = 20
numSamplesTest  = numSamples - numSamplesTrain

maximumIter = 500_000

# Data Visualization
gridNoiseStd = 0.05
numGridPts   = 500

Generate / Load Data#

In the following we’ll generate data according to the following model:

\[ y_{i} = f \left( x_{i} \right) + \epsilon_{i} \]

Where

\[ f \left( x \right) = \frac{1}{2} {x}^{3} + \frac{1}{4} {x}^{2} + 2 x + 5 \]

# The Data Generating Function

def f( vX: np.ndarray, vP: np.ndarray ):
    
    return np.polyval(vP, vX)

hF = lambda vX: f(vX, vP)

# Generate Data

vX = np.linspace(-5, 5, numSamples, endpoint = True) + (gridNoiseStd * np.random.randn(numSamples))
vN = noiseStd * np.random.randn(numSamples)
vY = hF(vX) + vN

print(f'The features data shape: {vX.shape}')
print(f'The labels data shape: {vY.shape}')

The features data shape: (30,)
The labels data shape: (30,)

Plot Data#

# Plot the Data

PlotRegressionData(vX, vY)

plt.show()

../../../../_images/15fb3ffcf48255499a764ec9ae5abd596c58e5085d4f142fc92968bedd5f20b5.png

Train Test Split#

In order to show the importance of the regularization, we’ll apply a split on the data.

# Train Test Split

mX = np.reshape(vX, (-1, 1))
mXTrain, mXTest, vYTrain, vYTest = train_test_split(mX, vY, test_size = numSamplesTest, train_size = numSamplesTrain, shuffle = True)

(?) What does stratify mean in the context of splitting data for regression?

Train a Regularized Polyfit Regressor#

The regularized PolyFit optimization problem is given by:

\[ \arg \min_{\boldsymbol{w}} {\left\| \boldsymbol{X} \boldsymbol{w} - \boldsymbol{y} \right|}_{2}^{2} + \lambda R \left( \boldsymbol{w} \right) \]