Regularized regression

Regularized regression is an important topic in modern statistics. For illustration purposes we consider the fixed design regression problem

y=Xβ+ε, (1.1)

with X ∈ R^n×d, β ∈ R^d and ε is an n-dimensional random vector with independent and identically distributed components. We assume throughout this chapter thatn≥d, i.e., we have

more observations than variables. Assuming that the columns ofX have mean zero and thaty is centred we denote the sample covariance matrix withA=n⁻¹X^TXand the cross covariance withb=n⁻¹X^Ty. The ordinary least squares estimatorβb_OLS =A⁻¹bis the minimizer (inβ) of the squared Euclidean distance betweenyandXβ. This estimator is unbiased and has several other important properties, see, e.g., Rao and Toutenburg (1999), chapter 3.

On the other hand it is obvious that the variance ofβb_OLS is high whenAis ill-conditioned. This problem is closely related to high collinearity in the columns ofX and thusAwill have small eigenvalues. As was mentioned in the previous section this occurs in the modelling of protein dynamics.

This complication can lead to unstable estimates of the coefficients ofβand, although the data used for model building can be estimated exactly, can lead to a poor generalization error also known as model overfitting (Hawkins, 2004).

When the quality of an estimatorβbis measured via the mean squared error the well known bias-variance decomposition can be used to analyze its behaviour. A biased estimator can improve upon βb_OLS in this sense if the variance is significantly lowered and at the same time the bias increases only slightly.

We consider estimators of the formβb_fθ = f_θ(A)b for a functionf_θ : [0,∞)→ Rthat depends on a parameter θ ∈ Θ ⊂ R. Usually fθ is chosen such that fθ(A) is better conditioned than A⁻¹. Here f_θ(A)is to be understood as the functional calculus of A, i.e., applying f_θ to the eigenvalues ofA. Of course forβb_OLS we have f_θ(x) = x⁻¹, x > 0. Typicallyf_θ has the role of a function that regularizesx⁻¹ and the degree of regularization depends on the regularization parameterθ.

We will first consider linear methods, that is, f_θ does not depend on y. Among this class of

methods are two of the most well known regression techniques, ridge regression and principal component regression. Partial least squares is a nonlinear regression technique and will be the focus of the next section.

Ridge regression (Hoerl and Kennard, 1970) is a biased method that is frequently used by statis-ticians when the regressor matrix is ill conditioned. It is also known in the literature of ill-posed problems as Tikhonov regularization (Tikhonov and Arsenin, 1977).

The regularization function isf_θ(x) = (x+θ)⁻¹,x≥0, for a parameterθ >0. For anyθ >0 the matrix A+θI_d, I_dbeing the d×didentity matrix, is invertible. Furthermore for small θ the perturbation of the original problem might be small enough that βb_θ^RR = fθ(A)b is a good estimator forβwith low variance.

It can be shown that the ridge estimator is the solution to the optimization problem min_v∈R^d∥Xv −y∥² +θ∥v∥² and thus large choices ofθ shrink the coefficients towards zero.

This hinders the regression estimates from blowing up like they can in the ordinary least squares estimator.

The simple description makes the theoretical analysis of ridge regression attractive. The opti-mality under a rotational invariant prior distribution on the coeffientsβ in Bayesian statistics (Frank and Friedman, 1993), make ridge regression a strong regularized regression technique when there is no prior belief on the size of the coefficientsβ. A major disadvantage is the need for the inversion of ad×dmatrix that can be quite cumbersome ifdis large. The choice ofθ is crucial, see Khalaf and Shukur (2005) for an overview of approaches.

Principal component regression is a technique that is based on principal component analysis (Pearson, 1901). Let us denote the eigenvalues of A withλ₁ ≥ λ₂ ≥ · · · ≥ λ_d ≥ 0 and the corresponding eigenvectorsv , . . . , v ∈R^d. Denote withIthe indicator function. For principal

component regression the functionf_a(x) = x⁻¹I(x ≥ λ_a),x > 0, is used with regularization parameter θ = a ∈ {1, . . . , d}, i.e., all eigenvalues that are smaller than λ_a are ignored for the inversion of A. This leads to the estimatorβb_a^{P CR} = f_a(A)b, a = 1, . . . , d, that avoids the collinearity problem ifais not chosen too large.

The principal component regression estimators can also be written as βb_a^{P CR} = W_a(W_a^TAW_a)⁻¹W_a^Tb. The matrix W_a = (w₁, . . . , w_a) ∈ R^d×a is calculated as follows. In the first step the aim is to find a vectorw₁ ∈R^dthat maximizes the empirical variance ofXv₁ and has unit norm, yieldingw1 = v1. Subsequent principal component vectors are calculated in the same way under the additional constraint that they are orthogonal tow₁, . . . , w_i−1. This givesw_i =v_i. See Jolliffe (2002) for details on the method.

Thus principal component regression also solves the problem of dimensionality reduction as we restrict our estimator to the space spanned by the first aeigenvectors. These eigenvectors are the ones that contribute most to the variance inX. For proteins this corresponds to the largest collective motions. To compute the principal component estimator it is necessary to calculate the first a eigenvectors of the matrix A, which, similar to the inversion in ridge regression, can be time intensive for large matrices. The number of used eigenvalues is crucial for the regularization properties of principal component regression and there are several ways to choose them, e.g., cross validation or the explained variance in the model.

We will mention some other methods, that are not necessarily linear iny, only shortly: the least absolute shrinkage and selection operator (Tibshirani, 1996), factor analysis (Gorsuch, 1983), least angle regression (Efron et al., 2004) and variable subset selection (Guyon and Elisseeff, 2003), to name only a few. We refer to Hastie et al. (2009) for an overview of the mentioned as well as other regularized regression methods.