Why Penalization Is Needed When P Exceeds N

Why the Number of Predictors Matters

In genomic prediction, a study may measure expression for thousands of genes in a cohort of only dozens or hundreds of samples. If the aim is to predict a quantitative outcome from those measurements, ordinary least squares (OLS) may seem like a natural starting point. But when the predictor count \(p\) exceeds the sample count \(n\), the data cannot determine a unique coefficient for every predictor. Even when a unique solution exists, strong dependence among gene-expression profiles can make the coefficients highly sensitive to small changes in the data.

These are distinct problems. Nonuniqueness means that different coefficient vectors give exactly the same fitted values on the observed samples. Instability means that small changes in the data can produce large changes in the estimated coefficients. Either can make gene-by-gene interpretations unreliable. Penalized regression methods, including ridge, lasso, and elastic net, are designed to constrain or stabilize coefficient estimates, although they do not remove every modeling difficulty.

OLS and the Design Matrix

Let \(X\) be the \(n \times p\) design matrix of predictor measurements, \(y\) the vector of \(n\) observed outcomes, and \(\beta\) a vector of \(p\) regression coefficients. For this discussion, the intercept is omitted or handled separately; the dimension argument applies to the predictor columns in either case. OLS chooses coefficients to minimize the residual sum of squares, the sum of squared differences between observed and fitted outcomes.

Definition: The rank of \(X\) is the number of linearly independent columns (equivalently, the dimension of its column space). The null space of \(X\) is the set of vectors \(v\) for which \(Xv=0\). If the rank is less than \(p\), the null space contains a nonzero vector, and the coefficients are not uniquely determined by the fitted values.

The rank of an \(n \times p\) matrix cannot exceed \(n\). Thus, when \(p>n\), the rank is strictly less than \(p\), so there is at least one nonzero vector \(v\) satisfying \(Xv=0\). This is a consequence of linear algebra, not a claim about the biological importance of any particular gene. It guarantees a coefficient ambiguity: moving along such a vector does not change the predictions for the observed samples.

\operatorname{rank}(X) \leq n < p \quad\Longrightarrow\quad \dim\{\text{null space of }X\}=p-\operatorname{rank}(X)>0

The OLS fitted values themselves are the projection of \(y\) onto the column space of \(X\). That projection is unique, even when the coefficients producing it are not. Consequently, saying “OLS has no solution” is generally imprecise: a least-squares fit exists, but its coefficient vector may not be unique. The distinction matters because predictions on the training samples can be well defined while the estimated contribution assigned to each gene is not.

Worked Example: Different Coefficients, Identical Training Fits

Consider a hypothetical expression dataset with two samples and three predictors. The intercept is omitted, and the design matrix and outcome vector are

X= \begin{pmatrix} 1&0&1\\ 0&1&1 \end{pmatrix}, \qquad y= \begin{pmatrix} 2\\3 \end{pmatrix}.

For coefficients \(\beta=(\beta_1,\beta_2,\beta_3)^\mathsf{T}\), an exact fit requires

\beta_1+\beta_3=2, \qquad \beta_2+\beta_3=3.

Setting \(\beta_3=t\) gives the full family of exact fits:

\beta= \begin{pmatrix} 2-t\\3-t\\t \end{pmatrix}, \qquad X\beta= \begin{pmatrix} (2-t)+t\\(3-t)+t \end{pmatrix} = \begin{pmatrix} 2\\3 \end{pmatrix}.

For \(t=0\), the coefficients are \((2,3,0)\); for \(t=1\), they are \((1,2,1)\). Both have residual sum of squares \(0\), so the observed data cannot distinguish them. Yet for a new sample with predictor vector \((0,0,1)^\mathsf{T}\), the two coefficient vectors predict \(0\) and \(1\), respectively. An exactly matching training fit does not by itself establish that the coefficients or predictions for new samples are stable.

Near-Collinearity and Instability

Having \(p>n\) guarantees nonuniqueness, but instability can occur even when the number of predictors does not exceed the number of samples. If two gene-expression profiles are almost linear combinations of one another, their columns contribute nearly redundant information. When \(X\) has full column rank, the OLS estimate is \(\hat{\beta}=(X^\mathsf{T}X)^{-1}X^\mathsf{T}y\). If \(X^\mathsf{T}X\) has a very small eigenvalue, its inverse has a very large eigenvalue in the corresponding direction, amplifying changes in \(y\).

For example, under the standard linear model with independent errors of mean zero and common variance \(\sigma^2\), and with full column rank, the coefficient covariance is \(\sigma^2(X^\mathsf{T}X)^{-1}\). This formula shows why small eigenvalues matter: they can create large sampling variance for coefficient combinations in poorly informed directions. In high-dimensional data, this difficulty can accompany exact nonuniqueness rather than merely occur near it.

Worked Example: A Small Outcome Change, a Larger Coefficient Change

Suppose two hypothetical predictors are measured on two samples:

X= \begin{pmatrix} 1&1\\ 1&1.01 \end{pmatrix}.

The columns are almost identical. With outcome \(y=(2,2)^\mathsf{T}\), the equations for an exact OLS fit are \(\beta_1+\beta_2=2\) and \(\beta_1+1.01\beta_2=2\). Subtracting the first equation from the second gives \(0.01\beta_2=0\), so \(\hat{\beta}=(2,0)^\mathsf{T}\).

Now change only the second outcome by \(0.001\), giving \(y'=(2,2.001)^\mathsf{T}\). The equations become \(\beta_1+\beta_2=2\) and \(\beta_1+1.01\beta_2=2.001\). Subtraction gives \(0.01\beta_2=0.001\), so \(\hat{\beta}'=(1.9,0.1)^\mathsf{T}\). A \(0.001\)-unit change in one outcome has changed the coefficient vector by \((-0.1,0.1)^\mathsf{T}\), whose length is about \(0.1414\). The fitted training values change only by the \(0.001\) added to the second outcome, but the allocation of the fitted signal between predictors changes much more.

The example is deliberately small and exactly solvable. Its lesson is not that every small outcome change causes this degree of coefficient movement, but that nearly redundant columns can make the allocation unstable. Coefficients may then be poor summaries of separate gene effects even when fitted values are relatively stable.

What Penalization Is Intended to Change

Penalized regression modifies the OLS fitting problem by discouraging some coefficient vectors more than others. The broad purpose is to choose a more constrained fit when the data alone cannot reliably identify all \(p\) coefficients. The specific objectives and their mathematical properties require separate treatment; here, the key is to understand the problem each family is intended to address.

Ridge regression penalizes the overall size of the coefficient vector. It can give a unique, stabilized solution in settings where the unpenalized coefficients are not unique, and it tends to retain many predictors with shrunken coefficients.
Lasso regression can set some coefficients exactly to zero, producing a sparse fitted model. When predictors are highly correlated, however, which predictor is retained can be sensitive to the data; sparsity does not guarantee that the chosen genes are uniquely or biologically identified.
Elastic net combines ridge-like and lasso-like penalties. It provides a way to balance coefficient shrinkage with sparsity, and can be useful when correlated predictors occur in groups. It does not guarantee that every correlated group will be selected or interpreted correctly.

Worked Example: What Redundant Expression Profiles Leave Undetermined

Suppose two hypothetical genes have identical expression values across two samples, \(x_1=(1,2)^\mathsf{T}\) and \(x_2=(1,2)^\mathsf{T}\), and the outcome is \(y=(2,4)^\mathsf{T}\). An exact fit requires

\beta_1 x_1+\beta_2 x_2 =(\beta_1+\beta_2) \begin{pmatrix}1\\2\end{pmatrix} = \begin{pmatrix}2\\4\end{pmatrix}, \qquad \beta_1+\beta_2=2.

Thus, \((2,0)\), \((1,1)\), and \((0,2)\) all give identical fitted values. Among these equally fitting choices, the squared coefficient sizes are \(4\), \(2\), and \(4\), respectively, so a ridge-type preference for smaller squared coefficients favors the balanced allocation \((1,1)\). The absolute coefficient sums are \(2\), \(2\), and \(2\), so a lasso-type preference for smaller absolute coefficients does not distinguish these nonnegative allocations by that criterion alone. A lasso fit can therefore be nonunique in this exact duplicate-predictor case. Elastic net includes a squared-size component as well as an absolute-size component, so its squared-size component favors balance among equally fitting allocations.

This example motivates the different penalty families without implying that a penalty discovers which gene is causally relevant. The two measured profiles contain no information that could separate their individual contributions.

Assumptions, Diagnostics, and Failure Modes

The linear algebra conclusions require no assumption that gene expression causes the outcome. They concern the geometry of \(X\) and the fitted linear model. For statistical interpretations based on a linear model, one often assumes that the conditional mean of the outcome is linear in the predictors and that errors have mean zero; classical variance calculations additionally require specified error-variance conditions. These assumptions should be considered separately from the rank problem. A penalty cannot repair a misspecified mean relationship, unmeasured confounding, batch effects, or biased sampling.

Useful checks include comparing \(p\) with \(n\), examining the rank or singular values of the design matrix, and checking how coefficients or selected predictors change under resampling or modest perturbations. In genomic data, large groups of highly correlated genes, batch-related structure, and small cohorts can all complicate interpretation. Cross-validation can assess predictive performance on held-out folds, but it does not prove that individual coefficients are stable or biologically meaningful.

If \(p>n\), the model can sometimes interpolate the training outcomes. Specifically, if \(X\) has rank \(n\), its columns span all of \(\mathbb{R}^n\), so every outcome vector \(y\) can be fitted exactly. In that case the training residual sum of squares is zero, leaving no residual degrees of freedom for the usual estimate based on training residuals. Exact training fit is not evidence of strong prediction on new samples.

A practical mistake is to treat a unique numerical answer returned by software as proof that the data uniquely identify the OLS coefficients. Software may use a generalized inverse or another convention to return one solution among many. That output can be computationally well defined while the underlying OLS minimizer is not unique. Conversely, when columns are nearly dependent but full rank, the coefficients can be unique yet still unstable. Always distinguish uniqueness from numerical stability.

The alternatives also have limitations. Ridge generally does not produce a sparse model; lasso may choose one of several correlated predictors in a way that changes across samples; and elastic net still depends on its penalty settings and the data available. No one method is automatically best merely because \(p>n\). The appropriate aim might be prediction, a compact set of candidate predictors, or stable description, and these aims need not select the same method. Penalization is a response to underdetermination and instability, not a substitute for careful study design, preprocessing, validation, or cautious interpretation.

Key takeaway: When \(p>n\), OLS coefficients cannot be uniquely determined from the observed design matrix. Near-collinearity can make coefficients unstable even without \(p>n\). Ridge, lasso, and elastic net constrain the fitting problem in different ways, but none can create information that the measurements do not contain.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions below.

With \(n=12\) samples and \(p=5000\) predictors, what is the minimum possible null-space dimension if the intercept is omitted?
If an \(8 \times 10\) design matrix has rank \(7\), how many independent coefficient directions leave all fitted values unchanged?
Does a full-rank design with nearly identical predictor columns have nonunique OLS coefficients? What problem can remain?
Why does a zero training residual sum of squares not establish that an OLS model will predict well for new genomic samples?
In one sentence each, state a basic motivation for ridge, lasso, and elastic net.

Solutions. (1) The rank is at most \(12\), so the null-space dimension is at least \(5000-12=4988\). (2) It is \(p-\operatorname{rank}(X)=10-7=3\). (3) No: full column rank gives unique OLS coefficients, but near-collinearity can make them unstable and sensitive to data changes. (4) The model may fit the observed outcomes without capturing patterns that generalize; if \(X\) spans \(\mathbb{R}^n\), any observed outcome vector can be fitted exactly. (5) Ridge shrinks overall coefficient size; lasso can produce sparse coefficients; elastic net combines shrinkage with a sparsity-oriented component and can be useful for correlated predictors.