NonLinear Regression & Lasso & Homotopy Method

Part One : Coordinate Descent for Lasso

The Lasso optimization problem can be formulated as
$$\hat{w}\in argmi{n}_{w\in R^d}\sum _{i=1}^m(h_w(x_i)-y_i{)}^2+\lambda \Vert w{\Vert }_1,$$
where $h_w(x)=w^{T}x$, and $\Vert w{\Vert }_1={\sum }_{i=1}^d|w_i|$. Note that to align with Murpy’s formulation below, and for historical reasons, we are using the total square loss, rather than the average square loss, in the objective function.

Since the ${\ell }_1$-regularization term in the objective function is non-differentiable, it’s not immediately clear how gradient descent or SGD could be used to solve this optimization problem directly.

Another approach to solving optimization problems is coordinate descent, in which at each step we optimize over one component of the unknown parameter vector, fixing all other components. The descent path so obtained is a sequence of steps, each of which is parallel to a coordinate axis in $R^d$, hence the name. It turns out that for the Lasso optimization problem, we can find a closed form solution for optimization over a single component fixing all other components. This gives us the following algorithm, known as the shooting algorithm.

The "soft thresholding’’ function is defined as
$$\text{soft}\left(a,\delta \right)=\text{sign}(a){\left(\left|a\right|-\delta \right)}_{+},$$
for any $a,\delta \in R$.

NOTE: The algorithm does not account for the case that $a_j=c_j=0$, which occurs when the $j$th column of $X$ is identically $0$. One can either eliminate the column (as it cannot possibly help the solution), or you can set $w_j=0$ in that case since it is, as you can easily verify, the coordinate minimizer. Note also that Murphy is suggesting to initialize the optimization with the ridge regession solution. Although theoretically this is not necessary (with exact computations and enough time, coordinate descent will converge for lasso from any starting point), in practice it’s helpful to start as close to the solution as we’re able.

There are a few tricks that can make selecting the hyperparameter $\lambda$ easier and faster. First, as we’ll see in a later problem, you can show that for any $\lambda \ge 2\Vert X^T(y-\bar{y}){\Vert }_{\infty }$, the estimated weight vector $\hat{w}$ is entirely zero, where $\bar{y}$ is the mean of values in the vector $y$, and $\Vert \cdot {\Vert }_{\infty }$ is the infinity norm (or supremum norm), which is the maximum over the absolute values of the components of a vector. Thus we need to search for an optimal $\lambda$ in $\left[0,{\lambda }_{\text{max}}\right]$, where ${\lambda }_{\text{max}}=2\Vert X^T(y-\bar{y}){\Vert }_{\infty }$. (Note: This expression for ${\lambda }_{\text{max}}$ assumes we have an unregularized bias term in our model. That is, our decision functions are of the form $h_{w,b}(x)=w^{T}x+b$. In our the experiments, we do not have an unregularized bias term, so we should use ${\lambda }_{\text{max}}=2\Vert X^{T}y{\Vert }_{\infty }$.)

The second trick is to use the fact that when $\lambda$ and ${\lambda }^{\prime }$ are close, the corresponding solutions $\hat{w}(\lambda )$ and $\hat{w}({\lambda }^{\prime })$ are also close. Start with $\lambda ={\lambda }_{\text{max}}$, for which
we know $\hat{w}({\lambda }_{\text{max}})=0$. You can run the optimization anyway, and initialize the optimization at $w=0$. Next, $\lambda$ is reduced (e.g. by a constant factor close to $1$), and the optimization problem is solved using the previous optimal point as the starting point. This is called warm starting the optimization. The
technique of computing a set of solutions for a chain of nearby $\lambda$'s is called a continuation or homotopy method. The resulting set of parameter values $\hat{w}(\lambda )$ as $\lambda$ ranges over $\left[0,{\lambda }_{\text{max }}\right]$ is known as a regularization path.

If we do not regularize, the loss function should be :
$$\begin{gathered} L(w)=(Y-Xw{)}^T(Y-Xw) \\ \frac{\partial L(w)}{\partial W_j}=-2\sum _{i=1}^N{X}_{ij}(y_i-w^T{X}_i) \\ =-2\sum _{i=1}^N{X}_{ij}(y_i-w_j{X}_{ij}-(w^T{X}_i-w_j{X}_{ij})) \\ =w_{j}2\sum _{i=1}^N{X}_{ij}^2-2\sum _{i=1}^N{X}_{ij}(y_i-w^T{X}_i+w_j{X}_{ij}) \end{gathered}$$

Let $a_j = 2\sum_{i=1}^N X_{ij}^{2}, c_j = 2\sum_{i=1}^N X_{ij}(y_i - w^TX_i + w_jX_{ij}) $ we can get :
$$\frac{\partial L(w)}{\partial W_j}=a_j{w}_j-c_j$$

Now, we add the regularization term :
$$L(w,\lambda )=L(w)+\lambda \Vert w{\Vert }_1$$

So now we get :
KaTeX parse error: Expected & or \\ or \cr or \end at position 204: … \ & w_j > 0 \ \̲ ̲\end{cases}

$$w_j=\{\begin{array}{ll}(c_j+\lambda )/a_j,& c_j<-\lambda ,\\ 0,& c_j\in [-\lambda ,\lambda ]\\ (c_j-\lambda )/a_j,& c_j>-\lambda \end{array}$$

We write the expression of $a_j , c_j $ in matrix form :
$$\begin{gathered} a_j=2X_{.j}^T{X}_{.j} \\ {c}_j=2X_{.j}^{T}y-2X_{.j}^{T}Xw+2w_j{X}_{.j}^T{X}_{.j} \end{gathered}$$

The let’s write a function that computes the Lasso solution for a given $\lambda$ using the shooting algorithm. For convergence criteria, continue coordinate descent until a pass through the coordinates reduces the objective function by less than $10^{-8}$, or have taken 1000 passes through the coordinates.

Import Modules

import numpy as np
import pandas as pd
from random import shuffle

Cyclic coordinate Descent for Lasso

def soft_f(a, delta) :
    if a >0 :
        result = max(abs(a)-delta, 0)
    elif a == 0 :
        result = 0
    else :
        result = -max(abs(a)-delta, 0)
    
    return result

def Lasso_Loss_fun(X, Y, w, lambda_reg) :
    temp = Y - np.dot(X, w)
    return np.dot(temp.T, temp) + lambda_reg*np.sum(np.absolute(w))

def Lasso_Co_Descent(X, Y, lambda_reg, init) :
    num_Co = X.shape[1]
    if init == 1 :
        w = np.dot(np.dot(np.linalg.inv(np.dot(X.T, X) +  lambda_reg*np.identity(num_Co)), X.T), Y)
    else :
        w = np.ones(num_Co)
    
    count = 0
    obj_fun_result = []
    obj_fun_result.append(Lasso_Loss_fun(X,Y,w,lambda_reg))
    
    while count <1000 :
        for i in range(num_Co) :
            if w[i] != 0 :
                a = 2*np.dot(X.T[i], X.T[i])
                c = 2*np.dot(X.T[i], Y) - 2*np.dot(np.dot(X.T[i], X), w) + 2*w[i]*np.dot(X.T[i], X.T[i])
                if(a == 0 or c==0) :
                    w[i] = 0
                else :
                    w[i] = soft_f(c/a, lambda_reg/a)
            else :
                w[i] = 0
                 
        count = count + 1
        obj_fun_result.append(Lasso_Loss_fun(X, Y, w, lambda_reg))
        if abs(obj_fun_result[count]-obj_fun_result[count-1]) < 1e-8  :
            break
            
    
    return w, np.array(obj_fun_result)

Randomized Coordinate Descent for Lasso

def Lasso_Co_Randomized_Descent(X, Y, lambda_reg, init) :
    num_Co = X.shape[1]
    if init == 1 :
        w = np.dot(np.dot(np.linalg.inv(np.dot(X.T, X) +  lambda_reg*np.identity(num_Co)), X.T), Y)
    else :
        w = np.ones(num_Co)
    
    count = 0
    obj_fun_result = []
    obj_fun_result.append(Lasso_Loss_fun(X,Y,w,lambda_reg))
    sequence = list(range(0, num_Co))
    
    while count<1000 :
        shuffle(sequence)
        for i in sequence :
            if w[i] != 0 :
                a = 2*np.dot(X.T[i], X.T[i])
                c = 2*np.dot(X.T[i], Y) - 2*np.dot(np.dot(X.T[i], X), w) + 2*w[i]*np.dot(X.T[i], X.T[i])
                if(a == 0 or c == 0) :
                    w[i] = 0
                else :
                    w[i] = soft_f(c/a, lambda_reg/a)
                    
            else :
                w[i] = 0
                    
        count = count + 1
        obj_fun_result.append(Lasso_Loss_fun(X, Y, w, lambda_reg))
        if abs(obj_fun_result[count]-obj_fun_result[count-1]) < 1e-8 :
            break
            
    
    return w, np.array(obj_fun_result)

Compare these two method :

from setup_problem import load_problem
from ridge_regression import RidgeRegression,plot_prediction_functions
#Load the training and validation data
lasso_data_fname = "lasso_data.pickle"
x_train, y_train, x_val, y_val, target_fn, coefs_true, featurize = load_problem(lasso_data_fname)

# Generate features
X_train = featurize(x_train)
X_val = featurize(x_val)

If we use Murphy Initial

w1, Or_Loss = Lasso_Co_Descent(X_train, y_train, 0.01, 1)
w2, Ra_Loss = Lasso_Co_Randomized_Descent(X_train, y_train, 0.01, 1)

w3, Or_Loss_new = Lasso_Co_Descent(X_train, y_train, 0.1, 1)
w4, Ra_Loss_new = Lasso_Co_Randomized_Descent(X_train, y_train, 0.1, 1)

import matplotlib.pyplot as plt

fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(Or_Loss, label = 'Cyclic $\lambda = 0.01$')
ax.plot(Ra_Loss, label = 'Random $\lambda = 0.01$')
plt.grid()
plt.xlabel('Number of Steps')
plt.ylabel('Training Lasso Loss Value')
plt.title('Training Lasso Loss: Cyclic Descsent and Murphy Initial')
plt.legend()

<matplotlib.legend.Legend at 0x1a1b408710>

fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(Or_Loss_new, label = 'Cyclic $\lambda = 0.1$')
ax.plot(Ra_Loss_new, label = 'Random $\lambda = 0.1$')
plt.grid()
plt.xlabel('Number of Steps')
plt.ylabel('Training Lasso Loss Value')
plt.title('Training Lasso Loss: Cyclic Descsent and Murphy Initial')
plt.legend()

<matplotlib.legend.Legend at 0x1a219a7f98>

If we do not use Murphy Initial

w5, Or_Loss_non = Lasso_Co_Descent(X_train, y_train, 0.01, 0)
w6, Ra_Loss_non = Lasso_Co_Randomized_Descent(X_train, y_train, 0.01, 0)

w7, Or_Loss_new_non = Lasso_Co_Descent(X_train, y_train, 0.1, 0)
w8, Ra_Loss_new_non = Lasso_Co_Randomized_Descent(X_train, y_train, 0.1, 0)

fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(Or_Loss_non[50:], label = 'Cyclic $\lambda = 0.01$')
ax.plot(Ra_Loss_non[50:], label = 'Random $\lambda = 0.01$')
plt.grid()
plt.xlabel('Number of Steps')
plt.ylabel('Training Lasso Loss Value')
plt.title('Training Lasso Loss: Cyclic Descsent and Non Initial')
plt.legend()

<matplotlib.legend.Legend at 0x1a22682ba8>

fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(Or_Loss_new_non[50:], label = 'Cyclic $\lambda = 0.1$')
ax.plot(Ra_Loss_new_non[50:], label = 'Random $\lambda = 0.1$')
plt.grid()
plt.xlabel('Number of Steps')
plt.ylabel('Training Lasso Loss Value')
plt.title('Training Lasso Loss: Cyclic Descsent and non Initial')
plt.legend()

<matplotlib.legend.Legend at 0x1a22965c50>

Grid Search for the best parameter

def compute_square_loss(X_train, Y_train, x_val, y_val, reg) :
    w, loss_list = Lasso_Co_Descent(X_train, y_train, reg, 1)
    residual = np.dot(x_val, w) - y_val
    square_loss = np.dot(residual.T, residual)
    return square_loss

square_loss_list = []
para_array = np.arange(0.1, 3, 0.05)

for i in range(len(para_array)) :
    square_loss_list.append(compute_square_loss(X_train, y_train, X_val, y_val, para_array[i]))

square_loss_array = np.array(square_loss_list)
square_loss_df = pd.DataFrame({'square loss' : square_loss_array}, index = para_array)
square_loss_df.plot()

<matplotlib.axes._subplots.AxesSubplot at 0x1a22cc7240>

Plot the Best Lasso Model together with Ridge and Target Function

w, loss_list = Lasso_Co_Descent(X_train, y_train, 1, 1)

from sklearn.linear_model import Ridge

plt.clf()
fig, ax = plt.subplots(figsize = (10,6))
plt.scatter(x_train, y_train, s=3, color='k', label='Training data')

x = np.arange(0,1,.001)
ax.plot(x, target_fn(x), 'r', label='Target function (i.e. Bayes prediction function)')

Ridge_estimator = Ridge(alpha = 2.35)
Ridge_estimator.fit(X_train, y_train)
x = np.sort(np.concatenate([np.arange(0,1,.001), x_train]))
X = featurize(x)
ax.plot(x,  Ridge_estimator.predict(X))
ax.plot(x,  np.dot(X, w))

plt.legend(['Target Function', 'Prediction with Ridge $\lambda$ = 2.35', 'Prediction with Lasso $\lambda = 1$', 'Training Data'], loc =1)
plt.title('Prediction Functions')
plt.xlabel('')
plt.show()

<Figure size 432x288 with 0 Axes>

Part Two : Homotopy Method

def avg_square_loss(X, Y, w) :
    residual = np.dot(X, w) - Y
    return np.dot(residual.T, residual)/X.shape[0] 

def Lasso_Co_Descent_Homotopy(X, Y, lambda_reg, W) :

    count = 0
    w = W
    obj_fun_result = []
    obj_fun_result.append(Lasso_Loss_fun(X,Y,w,lambda_reg))
    
    while count <1000 :
        for i in range(num_Co) :
            if w[i] != 0 :
                a = 2*np.dot(X.T[i], X.T[i])
                c = 2*np.dot(X.T[i], Y) - 2*np.dot(np.dot(X.T[i], X), w) + 2*w[i]*np.dot(X.T[i], X.T[i])
                if(a == 0 or c==0) :
                    w[i] = 0
                else :
                    w[i] = soft_f(c/a, lambda_reg/a)
            else :
                w[i] = 0
                 
        count = count + 1
        obj_fun_result.append(Lasso_Loss_fun(X, Y, w, lambda_reg))
        if abs(obj_fun_result[count]-obj_fun_result[count-1]) < 1e-8  :
            break
            
    
    return w

def homotopy_method(X_train, Y_train, x_val, y_val, lambda_max) :
    pre_w = np.ones(X_train.shape[1])
    w_list = []
    lambda_list = []
    avg_sqr_loss = []
    
    for i in range(30) :
        lambda_reg = lambda_max*((0.8)**(i))
        lambda_list.append(lambda_reg)
        curr_w = Lasso_Co_Descent_Homotopy(X_train, Y_train, lambda_reg, pre_w)
        w_list.append(curr_w)
        loss = avg_square_loss(x_val, y_val, curr_w)
        avg_sqr_loss.append(loss)
        pre_w = curr_w
        
    return lambda_list, avg_sqr_loss

lambda_list , avg_sqr_loss= homotopy_method(X_train, y_train, X_val, y_val, 2)

square_loss_array_H = np.array(avg_sqr_loss)
square_loss_H_df = pd.DataFrame({'avg square loss' : square_loss_array_H}, index = np.array(lambda_list))
square_loss_H_df.plot()

<matplotlib.axes._subplots.AxesSubplot at 0x1a23556b70>