Linear Regression, Ridge, and Lasso¶
Fit and compare linear regression models with quadratic and sparse penalties. Keep training and test data separate while inspecting coefficient paths and prediction errors, so that improved fit is not confused with improved generalization.
Run this tour¶
Run the cells in order with a Python 3 kernel. The first cell locates the companion data and toolbox and installs missing dependencies when needed. All worked examples include their implementation directly in this notebook. Random seeds make comparisons reproducible; you can change them to explore other samples.
# Locate the companion toolbox locally, or fetch it for a standalone/Colab copy.
from pathlib import Path
import importlib.util
import os
import subprocess
import sys
working = Path.cwd()
candidates = [working, working / "python", working.parent / "python"]
python_dir = next((p for p in candidates if (p / "nt_toolbox").is_dir()), None)
if python_dir is None:
checkout = working / "numerical-tours-support"
if not checkout.exists():
subprocess.run(
[
"git",
"clone",
"--depth",
"1",
"--branch",
"master",
"https://github.com/gpeyre/numerical-tours.git",
str(checkout),
],
check=True,
)
python_dir = checkout / "python"
os.chdir(python_dir)
if str(python_dir) not in sys.path:
sys.path.insert(0, str(python_dir))
requirements = python_dir / "requirements.txt"
if any(
importlib.util.find_spec(name) is None
for name in [
"numpy",
"scipy",
"matplotlib",
"skimage",
"sklearn",
"pywt",
"ipywidgets",
"cvxpy",
"skfmm",
"autograd",
"progressbar",
"celer",
]
):
subprocess.run(
[sys.executable, "-m", "pip", "install", "-r", str(requirements)], check=True
)
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(0)
plt.rcParams.update(
{
"figure.figsize": (8, 4),
"figure.dpi": 100,
"axes.spines.top": False,
"axes.spines.right": False,
"font.size": 11,
"image.cmap": "gray",
}
)
%matplotlib inline
$\newcommand{\dotp}[2]{\langle #1, #2 \rangle}$ $\newcommand{\enscond}[2]{\lbrace #1, #2 \rbrace}$ $\newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} }$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\umax}[1]{\underset{#1}{\max}\;}$ $\newcommand{\uargmin}[1]{\underset{#1}{argmin}\;}$ $\newcommand{\norm}[1]{\|#1\|}$ $\newcommand{\abs}[1]{\left|#1\right|}$ $\newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. }$ $\newcommand{\pa}[1]{\left(#1\right)}$ $\newcommand{\diag}[1]{{diag}\left( #1 \right)}$ $\newcommand{\qandq}{\quad\text{and}\quad}$ $\newcommand{\qwhereq}{\quad\text{where}\quad}$ $\newcommand{\qifq}{ \quad \text{if} \quad }$ $\newcommand{\qarrq}{ \quad \Longrightarrow \quad }$ $\newcommand{\ZZ}{\mathbb{Z}}$ $\newcommand{\CC}{\mathbb{C}}$ $\newcommand{\RR}{\mathbb{R}}$ $\newcommand{\EE}{\mathbb{E}}$ $\newcommand{\Zz}{\mathcal{Z}}$ $\newcommand{\Ww}{\mathcal{W}}$ $\newcommand{\Vv}{\mathcal{V}}$ $\newcommand{\Nn}{\mathcal{N}}$ $\newcommand{\NN}{\mathcal{N}}$ $\newcommand{\Hh}{\mathcal{H}}$ $\newcommand{\Bb}{\mathcal{B}}$ $\newcommand{\Ee}{\mathcal{E}}$ $\newcommand{\Cc}{\mathcal{C}}$ $\newcommand{\Gg}{\mathcal{G}}$ $\newcommand{\Ss}{\mathcal{S}}$ $\newcommand{\Pp}{\mathcal{P}}$ $\newcommand{\Ff}{\mathcal{F}}$ $\newcommand{\Xx}{\mathcal{X}}$ $\newcommand{\Mm}{\mathcal{M}}$ $\newcommand{\Ii}{\mathcal{I}}$ $\newcommand{\Dd}{\mathcal{D}}$ $\newcommand{\Ll}{\mathcal{L}}$ $\newcommand{\Tt}{\mathcal{T}}$ $\newcommand{\si}{\sigma}$ $\newcommand{\al}{\alpha}$ $\newcommand{\la}{\lambda}$ $\newcommand{\ga}{\gamma}$ $\newcommand{\Ga}{\Gamma}$ $\newcommand{\La}{\Lambda}$ $\newcommand{\Si}{\Sigma}$ $\newcommand{\be}{\beta}$ $\newcommand{\de}{\delta}$ $\newcommand{\De}{\Delta}$ $\newcommand{\phi}{\varphi}$ $\newcommand{\th}{\theta}$ $\newcommand{\om}{\omega}$ $\newcommand{\Om}{\Omega}$ $\newcommand{\eqdef}{\equiv}$
This tour studies linear regression method in conjunction with regularization. It contrasts ridge regression and the Lasso.
We recommend that after doing this Numerical Tours, you apply it to your own data, for instance using a dataset from LibSVM or Kaggle.
Disclaimer: these machine learning tours are intended to be overly-simplistic implementations and applications of baseline machine learning methods. For more advanced uses and implementations, we recommend to use a state-of-the-art library, the most well known being Scikit-Learn.
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
# Use this code to read from a CSV file.
# import pandas as pd
# U = pd.read_csv('myfile.csv')
Usefull functions to convert to a column/row vectors.
# convert to a column vector
def MakeCol(y):
return y.reshape(-1, 1)
# convert to a row vector
def MakeRow(y):
return y.reshape(1, -1)
# find non zero/true elements
def find(x):
return np.nonzero(x)[0]
Dataset Loading¶
We test the method on the prostate dataset in $n=97$ samples with features $x_i \in \RR^p$ in dimension $p=8$. The goal is to predict the price value $y_i \in \RR$.
Load the dataset.
from scipy import io
name = "prostate"
Data = io.loadmat("nt_toolbox/data/ml-" + name)
Ay = Data["A"]
class_names = Data["class_names"]
Randomly permute it.
Ay = Ay[np.random.permutation(Ay.shape[0]), :]
Separate the features $X$ from the data $y$ to predict information.
A_full = Ay[:, 0:-2]
y_full = MakeCol(Ay[:, -2])
c = MakeCol(Ay[:, -1])
Split into training and testing.
I0 = find(c == 1) # train
I1 = find(c == 0) # test
n = I0.size
n1 = I1.size
A = A_full[I0, :]
y = y_full[I0]
A1 = A_full[I1, :]
y1 = y_full[I1]
$n$ is the total number of samples, $p$ is the dimensionality of the features,
[n, p] = A.shape
print(n, p)
Normalize the features by the mean and std of the training set. This is optional.
mA = A.mean(axis=0)
sA = A.std(axis=0)
A = (A - mA) / sA
A1 = (A1 - mA) / sA
Remove the mean (computed from the test set) to avoid introducing a bias term and a constant regressor. This is optional.
m = y.mean()
y = y - m
y1 = y1 - m
Dimenionality Reduction and PCA¶
In order to display in 2-D or 3-D the data, dimensionality reduction is needed. The simplest method is the principal components analysis (PCA), which performs an orthogonal linear projection on the principal axes (eigenvectors) of the covariance matrix.
Display the covariance matrix of the training set.
C = A.transpose().dot(A)
plt.imshow(C);
u = A.transpose().dot(y)
plt.clf()
plt.bar(np.arange(1, p + 1), u.flatten())
plt.axis("tight");
Compute PCA ortho-basis and the feature in the PCA basis.
U, s, V = np.linalg.svd(A)
Ar = A.dot(V.transpose())
Plot sqrt of the eigenvalues.
plt.plot(s, ".-");
Display the features.
pmax = min(p, 8)
k = 0
plt.clf()
for i in np.arange(0, pmax):
for j in np.arange(0, pmax):
k = k + 1
plt.subplot(pmax, pmax, k)
if i == j:
plt.hist(A[:, i], 6)
plt.axis("tight")
else:
plt.plot(A[:, j], A[:, i], ".")
plt.axis("tight")
if i == 1:
plt.title(class_names[0][j][0])
plt.tick_params(axis="x", labelbottom=False)
plt.tick_params(axis="y", labelleft=False)
Display the points cloud of feature vectors in 2-D PCA space.
plt.plot(Ar[:, 0], Ar[:, 1], ".")
plt.axis("equal");
1D plot of the function to regress along the main eigenvector axes.
plt.clf()
for i in np.arange(0, 3):
plt.subplot(3, 1, i + 1)
plt.plot(Ar[:, i], y, ".")
plt.axis("tight")
Linear Regression¶
We look for a linear relationship $ y_i \approx \dotp{x}{a_i} $ written in matrix format $ y= A x $ where the rows of $A \in \RR^{n \times p}$ stores the features $a_i \in \RR^p$.
Since here $n > p$, this is an over-determined system, which can solved in the least square sense $$ \umin{ x } \norm{Ax-y}^2 $$ whose solution is given using the Moore-Penrose pseudo-inverse $$ x = (A^\top A)^{-1} A^\top y $$
Compute the least square solution.
x = np.linalg.solve(A.transpose().dot(A), A.transpose().dot(y))
Prediction (along 1st eigenvector).
plt.clf()
plt.plot(A1.dot(x), ".-")
plt.plot(y1, ".-")
plt.axis("tight")
plt.legend(("$y_1$", "$X_1 w$"));
Mean-square error on testing set.
E = np.linalg.norm(A1.dot(x) - y1) / np.linalg.norm(y1)
print(("Relative prediction error: " + str(E)));
Although this is not an effective method to solve this problem (a more efficient approach is to use for instance conjugate gradient), one can do a gradient descent to minimize the function $$ \min_x f(x) = \frac{1}{2}\norm{A x-y}^2. $$
def f(x):
return 1 / 2 * np.linalg.norm(A.dot(x) - y) ** 2
The gradient of $f$ is $$ \nabla f(x) = A^\top (Ax - y). $$
def Gradf(x):
return A.transpose().dot(A.dot(x) - y)
The maxium step size allowable by the gradient descent is $$ \tau \leq \tau_\max \eqdef \frac{2}{\norm{AA^\top}_{op}} $$ where $\norm{\cdot}_{op}$ is the maximum singular eigenvalue.
tau = 1 / np.linalg.norm(A, 2) ** 2
Initialize the algorithm to $x=0$.
x = np.zeros((p, 1))
One step of gradient descent reads: $$ x \leftarrow x - \tau \nabla f(x). $$
x = x - tau * Gradf(x)
tau_mult = [0.1, 0.5, 1, 1.8]
Worked example 0: Display the evolution of the training error $f(x)$ as a function of the number of iterations. Test for diffetent values of $\tau$.
niter = 100
flist = np.zeros((niter, 1))
tau_mult = [0.1, 0.5, 1, 1.8, 1.98]
xopt = np.linalg.solve(A.transpose().dot(A), A.transpose().dot(y))
plt.clf()
fig, (ax1, ax2) = plt.subplots(2, 1)
for itau in np.arange(0, 5):
tau = tau_mult[itau] / np.linalg.norm(A, 2) ** 2
x = np.zeros((p, 1))
for i in np.arange(0, niter):
flist[i] = f(x)
x = x - tau * Gradf(x)
# plt.subplot(2,1,1)
ax1.plot(flist)
ax1.axis("tight")
plt.title("f(x_k)")
# plt.subplot(2,1,2)
e = np.log10(flist - f(xopt) + 1e-20)
ax2.plot(e - e[0], label=str(tau_mult[itau]))
ax2.axis("tight")
leg = ax2.legend()
# ax2.legend( str( tau_mult[itau] ) )
plt.title("$log(f(x_k)-min J)$")
The optimal step size to minimize a quadratic function $\dotp{C w}{w}$ is $$ \tau_{\text{opt}} = \frac{2}{\sigma_\min(C) + \sigma_\max(C)} $$
C = A.transpose().dot(A)
tau_opt = 2 / (np.linalg.norm(C, 2) + np.linalg.norm(C, -2))
print(("Optimal tau = " + str(tau_opt * np.linalg.norm(A, 2) ** 2)) + " / |AA^T|");
Stochastic Gradient Method¶
We use SGD (which is not actually a descent algorithm) to minimize the quadratic risk $$ \umin{x} f(x) \eqdef \frac{1}{n} \sum_{i=1}^n f_i(x) = \frac{1}{n} \sum_{i=1}^n \frac{1}{2} ( \dotp{x}{a_i}-y_i )^2 $$ where we used $$ f_i(x) \eqdef \frac{1}{2} ( \dotp{x}{a_i}-y_i )^2$$ The algorithm reads $$ x_{k+1} = x_k - \tau_k \nabla f_{i_k}(x_k)$$ where at each iteration $i_k$ is drawn in $\{1,\ldots,n\}$ uniformly at random.
x = np.zeros((p, 1))
Draw $i_k$ are random.
ik = int(np.floor(np.random.rand() * n))
Compute $\nabla f_{i_k}(x)$.
gk = (A[ik, :].dot(x) - y[ik]) * A[ik, :].reshape(-1, 1)
Set the step size $\tau_k$ (for convergence is should converge to 0) and perform the update.
tauk = 1 / np.linalg.norm(A, 2) ** 2
x = x - tauk * gk
Worked example SGD1: Test different fixed step size $\tau_k=\tau$.
xopt = np.linalg.solve(
A.transpose().dot(A), A.transpose().dot(y)
) # least square solution
niter = 5000
flist = np.zeros((niter, 1))
tau_mult = [0.05, 0.3, 0.8]
for itau in np.arange(0, len(tau_mult)):
tauk = tau_mult[itau] / np.linalg.norm(A, 2) ** 2
x = np.zeros((p, 1))
for i in np.arange(0, niter):
ik = int(np.floor(np.random.rand() * n))
gk = (A[ik, :].dot(x) - y[ik]) * A[ik, :].reshape(-1, 1) # stochastic gradient
# gk = A.transpose().dot( A.dot(x)-y ) # batch gradient
x = x - tauk * gk
flist[i] = f(x)
plt.plot(flist / f(xopt) - 1)
plt.legend(("$\\tau=.05$", "$\\tau=.3$", "$\\tau=.8$"))
Worked example SGD2: Average on different runs to see the impact of the step size.
niter = 8000
tau_mult = [0.05, 0.8]
nruns = 10 # number of runs to compute the average performance
for itau in np.arange(0, len(tau_mult)):
tauk = tau_mult[itau] / np.linalg.norm(A, 2) ** 2
flist = np.zeros((niter, 1))
for iruns in np.arange(0, nruns):
x = np.zeros((p, 1))
for i in np.arange(0, niter):
ik = int(np.floor(np.random.rand() * n))
gk = (A[ik, :].dot(x) - y[ik]) * A[ik, :].reshape(-1, 1)
x = x - tauk * gk
flist[i] = flist[i] + f(x)
plt.plot((flist / nruns) / f(xopt) - 1)
plt.legend(("$\\tau=.05$", "$\\tau=.8$"))
Worked example SGD3: Use a decaying step size $\tau_k=\frac{\tau_0}{1+k/k_0}$.
niter = 20000
flist = np.zeros((niter, 1))
x = np.zeros((p, 1))
for i in np.arange(0, niter):
tauk = 1 / np.linalg.norm(A, 2) ** 2 * 1 / (1 + i / 10)
ik = int(np.floor(np.random.rand() * n))
gk = (A[ik, :].dot(x) - y[ik]) * A[ik, :].reshape(-1, 1)
x = x - tauk * gk
flist[i] = f(x)
plt.plot(flist / f(xopt) - 1)
print(flist[-1] / f(xopt) - 1)
Worked example SGD4: Use a decaying step size $\tau_k=\frac{\tau_0}{1+\sqrt{k/k_0}}$ and average the iteration $\frac{1}{K}\sum_{k<K}x_k$.
niter = 20000
flist = np.zeros((niter, 1))
x = np.zeros((p, 1))
x1 = np.zeros((p, 1))
for i in np.arange(0, niter):
tauk = 1 / np.linalg.norm(A, 2) ** 2 * 1 / (1 + np.sqrt(i / 10.0))
ik = int(np.floor(np.random.rand() * n))
gk = (A[ik, :].dot(x) - y[ik]) * A[ik, :].reshape(-1, 1)
x = x - tauk * gk
x1 = 1 / (i + 1) * x + i / (i + 1) * x1
flist[i] = f(x1)
plt.plot(flist / f(xopt) - 1);
Ridge Regularization¶
Regularization is obtained by introducing a penalty. It is often called ridge regression, and is defined as $$ \umin{ x } \norm{Ax-y}^2 + \lambda \norm{x}^2 $$ where $\lambda>0$ is the regularization parameter.
The solution is given using the following equivalent formula $$ x = (A^\top A + \lambda \text{Id}_p )^{-1} A^\top y, $$ $$ x = A^\top ( AA^\top + \lambda \text{Id}_n)^{-1} y, $$ When $p<n$ (which is the case here), the first formula should be prefered.
In contrast, when the dimensionality $p$ of the feature is very large and there is little data, the second is faster. Furthermore, this second expression is generalizable to Kernel Hilbert space setting, corresponding possibly to $p=+\infty$ for some kernels.
Lambda = 0.2 * np.linalg.norm(A) ** 2
x = np.linalg.solve(A.transpose().dot(A) + Lambda * np.eye(p), A.transpose().dot(y))
u = np.linalg.solve(A.dot(A.transpose()) + Lambda * np.eye(n), y)
x1 = A.transpose().dot(u)
print(("Error (should be 0): " + str(np.linalg.norm(x - x1) / np.linalg.norm(x))))
Worked example 1: Display the evolution of the test error $E$ as a function of $\lambda$.
q = 50
lmax = np.linalg.norm(A, 2) ** 2
lambda_list = lmax * np.linspace(0.3, 1e-3, q)
X = np.zeros((p, q))
E = np.zeros((q, 1))
for i in np.arange(0, q):
Lambda = lambda_list[i]
x = np.linalg.solve(A.transpose().dot(A) + Lambda * np.eye(p), A.transpose().dot(y))
X[:, i] = x.flatten() # bookkeeping
E[i] = np.linalg.norm(A1.dot(x) - y1) / np.linalg.norm(y1)
# find optimal lambda
i = E.argmin()
lambda0 = lambda_list[i]
xRidge = X[:, i]
print("Ridge: " + str(E.min() * 100) + "%")
# Display error evolution.
plt.clf()
plt.plot(lambda_list / lmax, E)
plt.plot([lambda0 / lmax, lambda0 / lmax], [E.min(), E.max()], "r--")
plt.axis("tight")
plt.xlabel(r"$\lambda/|X|^2$")
plt.ylabel("$E$")
Worked example 2: Display the regularization path, i.e. the evolution of $w$ as a function of $\lambda$.
plt.clf()
for i in np.arange(0, p):
plt.plot(lambda_list / lmax, X[i, :], label=class_names[0][i])
plt.plot(
[lambda0 / lmax, lambda0 / lmax], [X.flatten().min(), X.flatten().max()], "r--"
)
plt.axis("tight")
plt.xlabel(r"$\lambda/|X|^2$")
plt.ylabel("$x_i$")
plt.legend()
Sparse Regularization¶
In order to perform feature selection (i.e. select a subsect of the features which are the most predictive), one needs to replace the $\ell^2$ regularization penalty by a sparsity inducing regularizer. The most well known is the $\ell^1$ norm $$ \norm{x}_1 \eqdef \sum_i \abs{x_i} . $$
The energy to minimize is $$ \umin{x} f(x) \eqdef \frac{1}{2}\norm{Ax-y}^2 + \lambda \norm{x}_1. $$
def f(x, Lambda):
return 1 / 2 * np.linalg.norm(A.dot(x) - y) ** 2 + Lambda * np.linalg.norm(x, 1)
The simplest iterative algorithm to perform the minimization is the so-called iterative soft thresholding (ISTA), aka proximal gradient aka forward-backward.
It performs first a gradient step (forward) of the smooth part $\frac{1}{2}\norm{X w-y}^2$ of the functional and then a proximal step (backward) step which account for the $\ell^1$ penalty and induce sparsity. This proximal step is the soft-thresholding operator $$ \Ss_s(x) \eqdef \max( \abs{x}-\lambda,0 ) \text{sign}(x). $$
def Soft(x, s):
return np.maximum(abs(x) - s, np.zeros(x.shape)) * np.sign(x)
The ISTA algorithm reads $$ x_{k+1} \eqdef \Ss_{\la\tau}( x_k - \tau A^\top ( A x_k - y ) ), $$ where, to ensure convergence, the step size should verify $ 0 < \tau < 2/\norm{A}^2 $ where $\norm{A}$ is the operator norm.
Display the soft thresholding operator.
t = np.linspace(-5, 5, 201)
plt.clf()
plt.plot(t, Soft(t, 2))
plt.axis("tight");
Descent step size.
tau = 1.5 / np.linalg.norm(A, 2) ** 2
Choose a regularization parameter $\la$.
lmax = abs(A.transpose().dot(y)).max()
Lambda = lmax / 10
Initialization $w_0$.
x = np.zeros((p, 1))
A single ISTA step.
C = A.transpose().dot(A)
u = A.transpose().dot(y)
def ISTA(x, Lambda, tau):
return Soft(x - tau * (C.dot(x) - u), Lambda * tau)
x = ISTA(x, Lambda, tau)
Worked example 3: Implement the ISTA algorithm, display the convergence of the energy.
niter = 400
flist = np.zeros((niter, 1))
x = np.zeros((p, 1))
for i in np.arange(0, niter):
flist[i] = f(x, Lambda)
x = ISTA(x, Lambda, tau)
ndisp = int(niter / 4)
plt.clf()
plt.subplot(2, 1, 1)
plt.plot(flist[0:ndisp])
plt.axis("tight")
plt.title("f(x_k)")
plt.subplot(2, 1, 2)
e = np.log10(flist[0:ndisp] - flist.min() + 1e-20)
plt.plot(e - e[0])
plt.axis("tight")
plt.title("$log(f(x_k)-min f)$")
Worked example 4: Compute the test error along the full regularization path. You can start by large $\lambda$ and use a warm restart procedure to reduce the computation time. Compute the classification error. ind optimal lambda isplay error evolution.
q = 200
lambda_list = lmax * np.linspace(0.6, 1e-3, q)
X = np.zeros((p, q))
E = np.zeros((q, 1))
x = np.zeros((p, 1))
niter = 500
for iq in np.arange(0, q):
Lambda = lambda_list[iq]
# ISTA #
for i in np.arange(0, niter):
x = ISTA(x, Lambda, tau)
X[:, iq] = x.flatten() # bookkeeping
E[iq] = np.linalg.norm(A1.dot(x) - y1) / np.linalg.norm(y1)
# find optimal Lambda
i = E.argmin()
lambda0 = lambda_list[i]
xSparse = X[:, i]
print("Lasso: " + str(E.min() * 100) + "%")
# Display error evolution.
plt.clf()
plt.plot(lambda_list / lmax, E)
plt.plot([lambda0 / lmax, lambda0 / lmax], [E.min(), E.max()], "r--")
plt.axis("tight")
plt.xlabel(r"$\lambda/|A^* y|_\infty$")
plt.ylabel("$E$")
Worked example 5: Display the regularization path, i.e. the evolution of $w$ as a function of $\lambda$.
plt.clf()
for i in np.arange(0, p):
plt.plot(lambda_list / lmax, X[i, :], label=class_names[0][i])
plt.plot(
[lambda0 / lmax, lambda0 / lmax], [X.flatten().min(), X.flatten().max()], "r--"
)
plt.axis("tight")
plt.xlabel(r"$\lambda/|A^* y|_\infty$")
plt.ylabel("$x_i$")
plt.legend()
Worked example 6: Compare the optimal weights for ridge and lasso.
plt.clf()
plt.bar(np.arange(1, p + 1), abs(xSparse))
plt.bar(np.arange(1, p + 1), -abs(xRidge))
plt.legend(("Lasso", "Ridge"))
References and further reading¶
Trevor Hastie, Robert Tibshirani, and Jerome Friedman. The Elements of Statistical Learning. 2009, 2nd ed., Springer. Regression, classification, regularization, and model assessment.
Gilbert Strang. Linear Algebra and Learning from Data. 2019, Wellesley-Cambridge Press. Matrix factorizations, least squares, and low-rank representations.
Fabian Pedregosa et al.. Scikit-learn: Machine Learning in Python. 2011, Journal of Machine Learning Research 12, 2825–2830. Practical estimators and reproducible model evaluation.
Stephen Boyd and Lieven Vandenberghe. Convex Optimization. 2004, Cambridge University Press. Convexity, duality, optimality conditions, and interior-point methods.
Robert Tibshirani. Regression Shrinkage and Selection via the Lasso. 1996, Journal of the Royal Statistical Society B 58(1), 267–288. Sparse regression through an absolute-value penalty.