A Deep Reading of L1

A Bird Building Its Nest

Sometimes a beginning is strong: it catches your attention and leaves you searching for an answer throughout the story. At other times it seems ordinary, giving no hint of what will follow or what will stay with you.

You may wonder what this introduction has to do with a symbol or a story, but a metaphor can make the explanation flow. So imagine a bird trying to build its nest; it cannot carry every branch it finds, and must choose enough to build with.

We will return to it when numbers become choices, but let us stay with three ordinary numbers for now.

The 39 Python examples share state and are intended to run in order. Download the complete script and its requirements.

What Is $L_{1}$

Suppose we have three numbers: 3, −4, and 2. Put them in a vector and call it w.

Geometrically, that vector might describe a direction or displacement; its entries might also be weights learned by a model. Our question is simpler for now: how large are these numbers together?

We need a norm to measure the vector, and $L_{1}$ does this by adding the absolute value of each entry: remove its sign, then add the magnitudes. The index j identifies an entry and d counts the entries, so our vector has $L_{1}$ magnitude 9.

Run the code blocks in order; the first imports NumPy; other tools appear when we need them.

1import numpy as np2w = np.array([3., -4., 2.])3print(np.abs(w).sum())  # 9.0

The Club of Norms

The reason for taking absolute values appears if we use $(100, -100)$. Ordinary addition gives zero although both entries are large; $L_{1}$ gives 200.

Mathematicians require conditions before admitting a measurement into the club of norms: it must be nonnegative, vanish only at the zero vector, scale by |c| when the vector is multiplied by c, and obey the triangle inequality. That last condition says that the magnitude of a sum cannot exceed the sum of the magnitudes.

In distance terms, a direct route cannot be longer than a route through an intermediate point when both use the same norm. To make the scaling example precise, w=$(3, -2)$ gives 2w=$(6, -4)$, whose magnitude is 10, twice 5.

1w = np.array([3., -2.])2u, v = np.array([2., -1.]), np.array([-1., 3.])3print(2*w, np.abs(2*w).sum(), 2*np.abs(w).sum())4print(np.abs(u+v).sum() <= np.abs(u).sum()+np.abs(v).sum())

From Streets to Shapes

Geometrically, $L_{1}$ measures Manhattan distance. Reaching $(3, 2)$ from $(0, 0)$ through horizontal and vertical streets costs 5 units, whereas the straight line has length $\sqrt{13}$.

Collect all points whose norm is at most one, and $L_{1}$ gives a diamond, $L_{2}$ a disk, and $L_{\infty}$ a square. These shapes are more than decoration: their boundaries will later help determine where a solution can settle.

Still, calculating a norm alone changes no model; regularization begins when we give that magnitude a cost inside an optimization problem.

1point = np.array([3., 2.])2print("Manhattan:", np.abs(point).sum())3print("Euclidean:", np.linalg.norm(point))

The Family and Its Distant Cousin

Our protagonist belongs to a wider family. For $p\ge 1$, add the pth powers of the magnitudes and take the corresponding root.

At $p=1$ we obtain $L_{1}$, at $p=2$ Euclidean length, and as p grows without bound we obtain the largest individual magnitude. Its distant cousin $L_{0}$ counts nonzero entries, so it fails the homogeneity required of a true norm: doubling the entries does not double their count.

For $(3, 4, 0)$, the values are $L_{0}=2$, $L_{1}=7$, $L_{2}=5$, and $L_{\infty}=4$. Each answers a different question; two vectors can therefore have the same $L_{1}$ norm while using very different numbers of entries.

1w = np.array([3., 4., 0.])2print(np.count_nonzero(w), np.abs(w).sum(),3      np.linalg.norm(w), np.max(np.abs(w)))4a, b = np.array([1., 0., 0., 0.]), np.full(4, .25)5print(np.abs(a).sum(), np.abs(b).sum(),6      np.count_nonzero(a), np.count_nonzero(b))

Where Did the Story Begin

We may now ask who created $L_{1}$. Absolute value has no single beginning that can be assigned to an inventor: absolute-deviation estimation predates machine learning, and the development of $L_{p}$ spaces is a separate mathematical story.

Specific methods are easier to attribute. Hoerl and Kennard developed Ridge regression in 1970 to address unstable coefficients; Williams studied Laplace-prior regularization and pruning in neural networks in 1995; Tibshirani introduced Lasso in 1996, combining regression shrinkage with variable selection.

Chen, Donoho, and Saunders developed Basis Pursuit in 1998, followed by Elastic Net in 2005 and major compressed-sensing results in 2006. Each landmark added a question rather than erasing earlier work.

How Many Branches Do We Need

Return to the bird: choosing concerns the number of branches, not just their combined length. A vector is sparse when most entries are zero, and its support S(w) is the set of positions occupied by the remaining entries.

The support of $(0, 4, 0, -2, 0, 0)$ is {2, 4}, so its $L_{0}$ count is 2. To charge a fixed cost for every retained entry, we can add $\lambda L_{0}$ to the fitting loss F.

Searching over subsets can be difficult, since d positions allow $2^{d}$ subsets. In the one-coordinate example, compare the zero solution’s cost $a^{2}/2$ with the cost $\lambda$ of keeping a.

Choose zero when its cost is lower and keep a unchanged when $\lambda$ is lower; at equality both are optimal. Many small entries are not exact sparsity, although a sparse approximation may represent them well.

1w = np.array([0., 4., 0., -2., 0., 0.])2print("Support, counted from 1:", np.flatnonzero(w)+1)3a, lam = .6, .34hard_solution = 0. if a*a/2 <= lam else a5print("Hard threshold:", hard_solution)

What Does Convex Mean

We now need the word convex, because ease of optimization is not merely a description. Join two points on a convex function’s graph with a straight line, and the graph stays below or on that line.

The following inequality states this when a and b are mixed in proportion t. Every local minimum of a convex function on a convex domain is global, but this alone guarantees neither uniqueness nor existence.

Absolute value is convex despite its corner at zero, and squared regression loss is convex in the coefficients, so their sum in Lasso is convex. $L_{0}$ instead jumps at zero: mixing 0 and 1 equally gives cost 1 at the midpoint, larger than the average cost 1/2.

Adding a convex penalty to a nonconvex neural-network loss does not make the entire problem convex; nonsmoothness and nonconvexity are different properties.

1a, b, t = -2., 3., .42print("Absolute value:", abs(t*a+(1-t)*b), t*abs(a)+(1-t)*abs(b))3count = lambda x: float(x != 0)4print("L0 counterexample:", count(.5), .5*count(0)+.5*count(1))

A Replacement That Makes the Search Easier

This is why we sometimes replace difficult counting with a convex measurement, calling the replacement a convex relaxation. It creates a new problem; it does not prove that the new solution equals the old one.

On the bounded interval [−M, M], |w|/M is the largest convex function below the nonzero indicator: its convex envelope. Between zero and M, a convex lower bound cannot rise above the line joining $(0, 0)$ and (M, 1); scaled absolute value attains that line, and the negative side works similarly.

On a box bounding every entry, the envelope is $L_{1}$/M. On the whole unbounded space it is zero, rather than $L_{1}$.

Even a small example gives different solutions: a=0.6 and $\lambda=0.3$ make $L_{0}$ choose zero while $L_{1}$ chooses 0.3. Their relationship is not equality.

1M = 2.2grid = np.linspace(-M, M, 501)3indicator = (grid != 0).astype(float)4envelope = np.abs(grid)/M5print("Lower bound:", np.all(envelope <= indicator+1e-12))6print("L0 / L1 solutions:", hard_solution, max(.6-.3, 0))

When the Solution Touches a Corner

We can see the replacement’s effect before differentiating it. Give two coefficients an $L_{1}$ budget t: the allowed region is a diamond with vertices on the axes.

Equal-loss contours are ellipses for a positive-definite quadratic loss. We seek the smallest loss contour touching the allowed region, and contact at a vertex makes one coefficient zero.

In the following example, the unconstrained solution is $(2, 0.3)$. A budget of 1.2 gives the $L_{1}$ solution (1.2, 0), whereas $L_{2}$ gives approximately (1.1792, 0.2225).

Yet a diamond has edges too, and an optimum on an edge can have two nonzero coefficients. Under suitable Lagrangian conditions, a budget can be related to a penalty $\lambda$, but there is no universal fixed conversion between t and $\lambda$.

1from scipy.optimize import brentq2t = 1.23lagrange_l1 = brentq(lambda z: max(2-z, 0)+max(.3-z/2, 0)-t, 0, 3)4l1_point = np.array([max(2-lagrange_l1, 0), max(.3-lagrange_l1/2, 0)])5lagrange_l2 = brentq(lambda z: np.hypot(2/(1+z), .6/(2+z))-t, 0, 10)6l2_point = np.array([2/(1+lagrange_l2), .6/(2+lagrange_l2)])7print(l1_point, l2_point)

Why Does a Model Need Regularization

Why impose a budget or penalty at all? A model learns from examples it has seen, but we want it to perform on unseen examples.

It may memorize training details and noise, leading to overfitting. The average observed loss is empirical risk; population risk is expected loss under the distribution of future examples.

Regularization adds a preference $\Omega$ to fitting quality, with $\lambda$ setting its price. If model A has loss 0.20 and magnitude 10, while B has loss 0.23 and magnitude 2, then $\lambda=0.01$ gives objectives 0.30 and 0.25, favoring B.

Sacrificing some training fit may improve generalization, but does not guarantee it. Loss, regularizer, weighted penalty, and the full objective need to remain distinguishable.

1loss = np.array([.20, .23])2magnitude = np.array([10., 2.])3print(loss + .01*magnitude)  # [0.30, 0.25]

The Price of Magnitude: $\lambda$

When $\Omega$ is $L_{1}$, we charge the combined magnitudes of the coefficients we choose to regularize, which need not be every model parameter. We may exclude an intercept or some biases, so a strong penalty does not necessarily make the whole model zero. $\lambda$ is distinct from learning rate $\eta$, pruning threshold $\tau$, and stopping tolerance: they change the preference, the search step, a later removal decision, and numerical accuracy respectively.

Its effective price also depends on units and loss normalization. Multiplying both mean loss and penalty by n preserves the solution; changing the mean to a sum without scaling the penalty changes its strength.

For centered squared-error Lasso, λmax is the largest absolute column–response correlation divided by n, at which all slope coefficients can be zero. The two printed values agree because the penalty is also multiplied by the observation count; leaving $\lambda$ unchanged when switching from a mean to a sum changes the problem.

1losses = np.array([.4, .2, .3]); penalty, lam = 2., .052print(len(losses)*(losses.mean()+lam*penalty),3      losses.sum()+len(losses)*lam*penalty)

Data We Can Recreate

Let us examine data we can reproduce, using the research report’s synthetic example. x2 is strongly correlated with x1, x3 is noise, and the response depends on x1 plus additional noise. Each row of X is an observation, each column a predictor, and y is what we want to predict.

Fix the seeds to repeat the experiment. Learn column means and scales from training alone, then apply them to an independent reference sample; subtract the training response mean to fit slopes without penalizing an intercept.

These numbers are for teaching. In an actual study, the reference sample should not become data used repeatedly to tune the model.

1def make_data(n, seed):2    rng = np.random.default_rng(seed)3    x1 = rng.normal(size=n)4    x2 = .985*x1 + np.sqrt(1-.985**2)*rng.normal(size=n)5    x3 = rng.normal(size=n)6    return np.column_stack([x1, x2, x3]), 2*x1+.6*rng.normal(size=n)7Xraw, yraw = make_data(180, 12)8Xtest_raw, ytest = make_data(800, 31)9mu, scale, ybar = Xraw.mean(0), Xraw.std(0), yraw.mean()10X, y = (Xraw-mu)/scale, yraw-ybar11Xtest = (Xtest_raw-mu)/scale

From Least Squares to Ridge and Lasso

Ordinary least squares, OLS, seeks $\beta$ making Xβ close to y. Differentiation gives the normal equations; independent columns allow an inverse formula, but in practice a stable numerical least-squares solver is preferable, and dependent columns can give nonunique coefficients.

Ridge adds half $\lambda$ times squared $L_{2}$, shrinking weakly determined directions; Lasso adds $\lambda L_{1}$ and can set some coefficients to zero. Our data give approximately (2.0563, −0.1672, −0.0386) for OLS, (1.0626, 0.7356, −0.0325) for Ridge at 0.10, and (1.8211, 0, 0) for Lasso at 0.07.

Reference errors are close: 0.3862, 0.4051, and 0.3869. Different coefficient allocations need not mean substantially different predictions.

The code leaves solver details to the library so the comparison stays clear; we will unpack thresholding shortly. To see when every coefficient disappears, try multiples of λmax; at and above it, only the mean response remains in this centered model.

1from sklearn.linear_model import Lasso2beta_ols = np.linalg.lstsq(X, y, rcond=None)[0]3beta_ridge = np.linalg.solve(X.T@X+len(y)*.10*np.eye(3), X.T@y)4beta_lasso = Lasso(alpha=.07, fit_intercept=False,5                   tol=1e-10, max_iter=50000).fit(X, y).coef_6for name, beta in [("OLS", beta_ols), ("Ridge", beta_ridge), ("Lasso", beta_lasso)]:7    error = np.mean((ytest-ybar-Xtest@beta)**2)8    print(name, np.round(beta, 4), "test MSE:", round(error, 4))

1lambda_max = np.max(np.abs(X.T@y))/len(y)2for ratio in [.1, .5, 1., 1.2]:3    b = Lasso(alpha=ratio*lambda_max, fit_intercept=False,4              tol=1e-10, max_iter=50000).fit(X, y).coef_5    print("lambda/lambda_max:", ratio, "active:", np.count_nonzero(np.abs(b)>1e-8))

How Does Zero Become a Solution

The absolute-value corner is what makes zero possible. Solve the one-number problem so the idea does not hide inside matrices: minimize (w−a)²/2+λ|w|, where a is the value preferred by the data alone.

If w is positive, the derivative is w−a+$\lambda$, giving a−$\lambda$ when a>$\lambda$. If w is negative, the derivative is w−a−$\lambda$, giving a+$\lambda$ when a<−$\lambda$.

At zero there is no single derivative, but there is a set of supporting slopes [−1, 1], the convex subdifferential. Zero is optimal exactly when |a|≤λ.

Soft-thresholding therefore removes values inside this interval and subtracts $\lambda$ from the magnitudes that survive. Ridge instead gives a/(1+$\lambda$), ordinarily shrinking a nonzero value without removing it at finite $\lambda$.

1def soft(v, threshold):2    return np.sign(v) * np.maximum(np.abs(v)-threshold, 0.0)3a = np.array([-2., -.5, 0., .5, 2.]); lam = .84print("OLS:", a)5print("Ridge:", a/(1+lam))6print("Lasso:", soft(a, lam))

The Balance at the Solution

Returning to a vector, we examine each coefficient’s relationship with the residual r=y−Xβ after joint fitting. The KKT conditions say that a column’s residual correlation equals $\lambda$ times the sign of its nonzero coefficient; for a zero coefficient, the correlation may lie anywhere in [−$\lambda$, $\lambda$].

This balances fitting pressure against penalty pressure and certifies optimality for the convex problem. Coordinate descent isolates the current column’s effect, forms its partial residual, and computes correlation c_j and average squared column magnitude q_j.

The update is Sλ(c_j)/q_j. The solver does not guess a zero: the thresholding subproblem creates it.

A printed value near zero means the balance conditions hold numerically; it does not by itself tell us how the model will perform on new data.

1r = y-X@beta_lasso2corr = X.T@r/len(y)3active = np.abs(beta_lasso) > 1e-84violations = np.where(active, np.abs(corr-.07*np.sign(beta_lasso)),5                      np.maximum(np.abs(corr)-.07, 0))6print("Largest KKT violation:", violations.max())

A Step That Reaches Zero

Following a slope directly may jump across the corner instead of landing on it. For |w|, starting at 0.03 with step 0.1 can alternate between 0.03 and −0.07.

Proximal gradient separates movement that reduces the smooth loss from handling the penalty: take a gradient step, then apply threshold ηλ. In the example below, (0.04, −0.20, 0) first becomes (0.03, −0.15, 0.06), then $(0, -0.11, 0.02)$.

One entry disappears, another shrinks, and a previously zero entry returns. Zero is not permanent when the fitting gradient changes.

For a smooth convex loss with L-Lipschitz gradient, the usual condition 0<$\eta\le 1$/L supports the appropriate convergence results.

1weights = np.array([.04, -.20, 0.])2gradient = np.array([.10, -.50, -.60])3intermediate = weights-.1*gradient4print(intermediate, "->", soft(intermediate, .1*.4))

How Accurate Is the Solution

After finding a solution, we need to know its numerical accuracy, a different question from prediction quality. In squared-error Lasso, a feasible dual value provides a lower bound for the original optimum; the primal–dual gap bounds remaining objective error.

The code builds a feasible dual vector from the residual and calculates that gap. Approaching zero means approaching the solution of this optimization problem, without certifying prediction quality.

Here both objectives approach approximately 0.3200; rounding can leave a tiny difference. The gap can also support safe screening of columns, but here we stay with our first question: has the numerical search finished?

1lam, n = .07, len(y)2r = y-X@beta_lasso3u = r/max(1., np.max(np.abs(X.T@r))/(n*lam))4primal = r@r/(2*n)+lam*np.abs(beta_lasso).sum()5dual = y@u/n-u@u/(2*n)6print("Primal / dual:", primal, dual)7print("Gap:", primal-dual)  # Near zero; tiny rounding errors are possible.

Choosing $\lambda$

Even an accurate solution can use an unsuitable $\lambda$, so choose regularization strength with observations kept out of fitting. Cross-validation divides training into folds, fits on some and evaluates on the rest.

Standardization must be learned within each training fold to prevent validation information leaking into preprocessing. Pipeline joins the steps and GridSearchCV searches $\lambda$, which Lasso calls alpha under the same normalization used here.

This example selects about 0.0297, but that number is not a general recommendation. If we estimate performance after tuning without an independent test sample, an outer evaluation such as nested validation is needed.

Temporal or grouped data need splits respecting their time or group structure.

1from sklearn.linear_model import Lasso2from sklearn.model_selection import GridSearchCV, KFold3from sklearn.pipeline import Pipeline4from sklearn.preprocessing import StandardScaler5procedure = Pipeline([("scale", StandardScaler()),6                      ("fit", Lasso(max_iter=50000, tol=1e-10))])7search = GridSearchCV(procedure,8    {"fit__alpha": np.geomspace(.003, .5, 30)},9    scoring="neg_mean_squared_error", cv=KFold(5, shuffle=True, random_state=20))10search.fit(Xraw, yraw)11print(search.best_params_, np.mean((ytest-search.predict(Xtest_raw))**2))