Dual Ridge Regression When P Exceeds N

Why Solve Ridge Regression in Sample Space?

In high-dimensional genomic data, the number of measured features can be far larger than the number of samples. For example, an expression matrix might have thousands of gene columns but only a few dozen sample rows. Ridge regression remains well-defined in this setting, but its usual feature-space formula involves a \(p\)-by-\(p\) matrix. The dual formulation instead solves a system involving an \(n\)-by-\(n\) matrix, which can be much smaller when \(p\) greatly exceeds \(n\).

The earlier material on the penalized least-squares objective and ridge shrinkage establishes the optimization setup and its relationship to the design matrix. Here the focus is the algebra that moves the computation from feature space to sample space, and what that change does—and does not—save.

Setup and Assumptions

Let \(X\) be an \(n\)-by-\(p\) design matrix, with one row per sample and one column per predictor, and let \(y\) be the \(n\)-vector of responses. For the centered formulation below, each column of \(X\) and the response vector \(y\) have mean zero. If a model includes an intercept, center the predictors and response, estimate the slopes from the centered data, and recover the intercept as \(\bar y-\bar x^\mathsf{T}\hat\beta\). The intercept is not penalized.

Assume that the ridge penalty \(\lambda\) is strictly positive. No full-rank condition on \(X\) is needed. The derivation is an optimization identity for a fixed design matrix; it does not require normally distributed errors. The ridge objective used here is the unaveraged sum of squares plus penalty. If an objective instead uses the mean squared error, its penalty parameter must be scaled accordingly to represent the same fit.

The penalty makes the minimizer unique: \(X^\mathsf{T}X\) is positive semidefinite, so for every nonzero \(v\), \(v^\mathsf{T}(X^\mathsf{T}X+\lambda I_p)v=\|Xv\|_2^2+\lambda\|v\|_2^2>0\). Thus \(X^\mathsf{T}X+\lambda I_p\) is invertible even when \(p>n\) or predictors are collinear. Likewise, \(K+\lambda I_n\) is invertible.

Deriving the Dual Solution

Differentiating the objective with respect to \(\beta\) and setting the gradient to zero gives the ridge normal equations. Their solution is the feature-space formula:

To obtain the sample-space form, first solve the \(n\)-dimensional system \((XX^\mathsf{T}+\lambda I_n)\alpha=y\), then set \(\hat\beta=X^\mathsf{T}\alpha\). This is valid because substituting that coefficient vector into the feature-space normal equations yields:

The feature-space normal equations have a unique solution, so \(X^\mathsf{T}\alpha\) must be that solution. This proves the identity

The fitted response vector has an especially direct sample-space expression. With \(K=XX^\mathsf{T}\), it is \(\hat y=X\hat\beta=K(K+\lambda I_n)^{-1}y\). Since \(K\) commutes with \(K+\lambda I_n\), the residual vector is \(y-\hat y=\lambda(K+\lambda I_n)^{-1}y\). In computation, solve the linear system for \(\alpha\); do not explicitly calculate a matrix inverse.

Dual Coordinates and Predictions

The vector \(\alpha\) has one entry per training sample, rather than one entry per predictor. It is a set of weights that expresses the fitted coefficient vector as a linear combination of the training predictor columns: \(\hat\beta=\sum_{i=1}^n\alpha_i x_i\), where \(x_i^\mathsf{T}\) is row \(i\) of \(X\). This representation explains why prediction can also be carried out without first storing the \(p\) coefficients.

For a new centered predictor vector \(x_{\mathrm{new}}\) of length \(p\), its predicted response is \(x_{\mathrm{new}}^\mathsf{T}\hat\beta\). Substitution of the dual expression gives \(x_{\mathrm{new}}^\mathsf{T}X^\mathsf{T}\alpha=(Xx_{\mathrm{new}})^\mathsf{T}\alpha\). The vector \(Xx_{\mathrm{new}}\) contains the inner products between the new sample and each training sample. If an intercept was fitted, add the recovered intercept to this centered prediction.

When the Sample-Space Solve Is Cheaper

The feature-space formula requires solving a system with \(p\) rows and columns; the dual formula requires solving one with \(n\) rows and columns. For dense matrices, a direct factorization has cubic cost in the system dimension. Forming the Gram matrix also has a cost: \(XX^\mathsf{T}\) takes on the order of \(n^2p\) arithmetic operations, while forming \(X^\mathsf{T}X\) takes on the order of \(p^2n\). The dual approach is typically attractive when \(p\) greatly exceeds \(n\), but forming the matrix and producing predictions still require work. Both approaches must access the design matrix, and obtaining all \(p\) coefficients requires an additional multiplication by \(X^\mathsf{T}\).

Checks, Alternatives, and Limitations

The feature-space and sample-space formulas are algebraically identical for the same \(X\), \(y\), and \(\lambda\). Choose the system with the smaller dimension as a first computational guide, then consider whether the Gram matrix can be formed and stored efficiently. The feature-space route may be preferable when \(p\) is not large relative to \(n\), or when coefficients are the primary object and the feature-space solve is manageable. The dual route is often preferable when \(p\gg n\), especially when predictions can be computed from sample inner products.

Neither representation removes the modeling choices behind ridge regression. The value of \(\lambda\) still governs shrinkage, predictor scaling still affects how a common penalty acts on coefficients, and a training fit alone does not establish predictive performance. Keep the penalty convention consistent: changing from a sum of squares to an average changes the numerical scale of \(\lambda\). If features or the response were centered for fitting, apply the same training-derived centering to new observations.

Useful diagnostics include checking that the Gram matrix is symmetric up to numerical precision, verifying that the computed solution satisfies \((K+\lambda I_n)\alpha\approx y\), and comparing \(X^\mathsf{T}\alpha\) with a feature-space solution on a small test problem. Large discrepancies usually indicate inconsistent centering, mismatched penalty scaling, a transposed design matrix, or a faulty linear solve. A very small positive \(\lambda\) can also make the solve numerically sensitive when the data contain nearly redundant directions; use a stable factorization and inspect the residual rather than explicitly forming an inverse.