Rank deficiency

It might happen that the columns of \(\mathbf{X}\) are not linearly independent, so that \(\mathbf{X}\) is not of full rank. This would occur, for example, if two of the inputs were perfectly correlated, e.g., \(\mathbf{x}_2=3\mathbf{x}_1\).

Then \(\mathbf{X}^T\mathbf{X}\) is singular and the least squares coefficients \(\hat\beta\) are not uniquely defined. However, the fitted values \(\hat{\mathbf{y}}=\mathbf{X}\hat\beta\) are still the projection of \(\mathbf{y}\) onto the \(\text{col}(\mathbf{X})\); there are just more than one way to express that projection in terms of the column vectors of \(\mathbf{X}\).

The non-full-rank case occurs most often when one or more qualitative inputs are coded in a redundant fashion.

There is usually a natural way to resolve the non-unique representation, by recording and/or dropping redundant columns in \(\mathbf{X}\). Most regression software packages detect these redundancies and automatically implement some strategy for removing them.

Rank deficientcies can also occur in signal and image analysis, where the number of inputs \(p\) can exceed the number of training cases \(N\). In this case, the features are typically reduced by filtering or else the fitting is controlled by regularization (\(\S\) 5.2.3 and Chapter 18).

The statement of the theorem

We focus on estimation of any linear combination of the parameters \(\theta=a^T\beta\). The least squares estimate of \(a^T\beta\) is

\[\begin{equation} \hat\theta = a^T\hat\beta = a^T\left(\mathbf{X}^T\mathbf{X}\right)^{-1}\mathbf{X}^T\mathbf{y}. \end{equation}\]

Considering \(\mathbf{X}\) to be fixed and the linear model is correct, \(a^T\beta\) is unbiased since

\[\begin{align} \text{E}(a^T\hat\beta) &= \text{E}\left(a^T(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}\right) \\ &= a^T(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{X}\beta \\ &= a^T\beta \end{align}\]

The Gauss-Markov Theorem states that if we have any other linear estimator \(\tilde\theta = \mathbf{c}^T\mathbf{y}\) that is unbiased for \(a^T\beta\), that is, \(\text{E}(\mathbf{c}^T\mathbf{y})=a^T\beta\), then

\[\begin{equation} \text{Var}(a^T\hat\beta) \le \text{Var}(\mathbf{c}^T\mathbf{y}). \end{equation}\]

The proof (Exercise 3.3) uses the triangle inequality.

For simplicity we have stated the result in terms of estimation of a single parameter \(a^T \beta\), but with a few more definitions one can state it in terms of the entire parameter vector \(\beta\) (Exercise 3.3).

Implications of the Gauss-Markov theorem

Consider the mean squared error of an estimator \(\tilde\theta\) of \(\theta\):

\[\begin{align} \text{MSE}(\tilde\theta) &= \text{E}\left(\tilde\theta-\theta\right)^2 \\ &= \text{E}\left(\tilde\theta-\text{E}(\tilde\theta)+\text{E}(\tilde\theta)-\theta\right)^2\\ &= \text{E}(\tilde\theta-\text{E}(\tilde\theta))^2+\text{E}(\text{E}(\tilde\theta)-\theta)^2+2\text{E}(\tilde\theta-\text{E}(\tilde\theta))(\text{E}(\tilde\theta)-\theta)\\ &= \text{Var}\left(\tilde\theta\right) + \left[\text{E}\left(\tilde\theta\right)-\theta\right]^2 \\ &= \text{Var} + \text{Bias}^2 \end{align}\]

The Gauss-Markov theorem implies that the least squares estimator has the smallest MSE of all linear estimators with no bias. However there may well exist a biased estimator with smaller MSE. Such an estimator would trade a little bias for a larger reduction in variance.

Biased estimates are commonly used. Any method that shrinks or sets to zero some of the least squares coefficients may result in a biased estimate. We discuss many examples, including variable subset selection and ridge regression, later in this chapter.

From a more pragmatic point of view, most models are distortions of the truth, and hence are biased; picking the right model amounts to creating the right balance between bias and variance. We go into these issues in more detail in Chapter 7.

Relation between prediction accuracy and MSE

MSE is intimately related to prediction accuracy, as discussed in Chapter 2.

Consider the prediction of the new response at input \(x_0\),

\[\begin{equation} Y_0 = f(x_0) + \epsilon. \end{equation}\]

Then the expected prediction error of an estimate \(\tilde{f}(x_0)=x_0^T\tilde\beta\) is

\[\begin{align} \text{E}(Y_0 - \tilde{f}(x_0))^2 &= \text{E}\left(Y_0 -f(x_0)+f(x_0) - \tilde{f}(x_0)\right)^2\\ &= \sigma^2 + \text{E}\left(x_o^T\tilde\beta - f(x_0)\right)^2 \\ &= \sigma^2 + \text{MSE}\left(\tilde{f}(x_0)\right). \end{align}\]

Therefore, expected prediction error and MSE differ only by the constant \(\sigma^2\), representing the variance of the new observation \(y_0\).

The first error component \(\sigma^2\) is unrelated to what model is used to describe our data. It cannot be reduced for it exists in the true data generation process. The second source of error corresponding to the term \(\text{MSE}(\tilde{f}(x_0))\) represents the error in the model and is under control of the statistician. Thus, based on the above expression, if we minimize the \(\text{MSE}\) of our estimator \(\tilde{f}(x_0)\) we are effectively minimizing the expected (quadratic) prediction error which is our ultimate goal anyway.

We will explore methods that minimize the mean square error. The mean square error can be broken down into two terms: a model variance term and a model bias squared term. We will explore methods that seek to keep the total contribution of these two terms as small as possible by explicitly considering the trade-offs that come from methods that might increase one of the terms while decreasing the other.

Building blocks for multiple linear regression

Suppose next that the columns of the data matrix \(\mathbf{X} = \left[\mathbf{x}_1,\cdots,\mathbf{x}_p\right]\) are orthogonal, i.e.,

\[\begin{equation} \langle \mathbf{x}_j,\mathbf{x}_k\rangle = 0\text{ for all }j\neq k. \end{equation}\]

Then it is easy to check that the multiple least squares estimates are equal to the univariate estimates:

\[\begin{equation} \hat\beta_j = \frac{\langle\mathbf{x}_j,\mathbf{y}\rangle}{\langle\mathbf{x}_j,\mathbf{x}_j\rangle}, \forall j \end{equation}\]

In other words, when the inputs are orthogonal, they have no effect on each other’s parameter estimates in the model.

Orthogonal inputs occur most often with balanced, designed experiments (where orthogonality is enforced), but almost never with observational data. Hence we will have to orthogonalize them in order to carry this idea further.

Gram-Schmidt procedure and QR decomposition

Algorithm 3.1 is known as the Gram-Schmidt procedure for multiple regression, and is also a useful numerical strategy for computing the estimates. We can obtain from it not just \(\hat\beta_p\), but also the entire multiple least squares fit (Exercise 3.4).

We can represent step 2 of Algorithm 3.1 in matrix form:

\[\begin{equation} \mathbf{X} = \mathbf{Z\Gamma}, \end{equation}\]

where

• \(\mathbf{Z}\) has as columns the \(\mathbf{z}_j\) (in order)
• \(\mathbf{\Gamma}\) is the upper triangular matrix with entries \(\hat\gamma_{kj}\).

Introducing the diagonal matrix \(\mathbf{D}\) with \(D_{jj}=\|\mathbf{z}_j\|\), we get

\[\begin{align} \mathbf{X} &= \mathbf{ZD}^{-1}\mathbf{D\Gamma} \\ &= \mathbf{QR}, \end{align}\]

the so-called QR decomposition of \(\mathbf{X}\). Here

• \(\mathbf{Q}\) is an \(N\times(p+1)\) orthogonal matrix s.t. \(\mathbf{Q}^T\mathbf{Q}=\mathbf{I}\),
• \(\mathbf{R}\) is a \((p+1)\times(p+1)\) upper triangular matrix.

The QR decomposition represents a convenient orthogonal basis for the \(\text{col}(\mathbf{X})\). It is easy to see, for example, that the ordinary least squares (OLS) solution is given by

\[\begin{align} \hat\beta &= \left(\mathbf{X}^T\mathbf{X}\right)^{-1}\mathbf{X}^T\mathbf{y}\\ &=\left(\mathbf{R}^T\mathbf{Q}^T\mathbf{Q}\mathbf{R}\right)^{-1}\mathbf{R}^T\mathbf{Q}^T\mathbf{y}\\ &=\left(\mathbf{R}^T\mathbf{R}\right)^{-1}\mathbf{R}^T\mathbf{Q}^T\mathbf{y}\\ &=\mathbf{R}^{-1}\mathbf{R}^{-T}\mathbf{R}^T\mathbf{Q}^T\mathbf{y}\\ &=\mathbf{R}^{-1}\mathbf{Q}^T\mathbf{y}, \\ \hat{\mathbf{y}} &= \mathbf{X}\hat\beta =\mathbf{Q}\mathbf{R}\mathbf{R}^{-1}\mathbf{Q}^T\mathbf{y}=\mathbf{QQ}^T\mathbf{y}. \end{align}\]

This last equation expresses the fact in ordinary least squares we obtain our fitted vector \(\mathbf y\) by first computing the coefficients of \(\mathbf y\) in terms of the basis spanned by the columns of \(\mathbf Q\) (these coefficients are given by the vector \(\mathbf Q^T\mathbf y\)). We next construct \(\hat{\mathbf y}\) using these numbers as the coefficients of the column vectors in \(\mathbf Q\) (this is the product \(\mathbf{QQ}^T\mathbf y\)).

Note that the triangular matrix \(\mathbf{R}\) makes it easy to solve (Exercise 3.4).

Correlated errors

If the errors \(\epsilon = (\epsilon_1,\cdots,\epsilon_K)\) are correlated with \(\text{Cov}(\epsilon)=\mathbf{\Sigma}\), then the multivariate weighted criterion

\[\begin{equation} \text{RSS}(\mathbf{B};\mathbf{\Sigma}) = \sum_{i=1}^N (y_i-f(x_i))^T \mathbf{\Sigma}^{-1} (y_i-f(x_i)) \end{equation}\]

arises naturally from multivariate Gaussian theory. Here

• \(f(x) = \left(f_1(x),\cdots,f_K(x)\right)^T\) is the vector function,
• \(y_i\) the vector of \(K\) responses for observation \(i\). However, the solution is again the same with ignoring the correlations as

\[\begin{equation} \hat{\mathbf{B}} = \left(\mathbf{X}^T\mathbf{X}\right)^{-1}\mathbf{X}^T\mathbf{Y}. \end{equation}\]

In Section 3.7 we pursue the multiple output problem and consider situations where it does pay to combine the regressions.

Forward-stepwise selection

Forward-stepwise selection starts with the intercept, and the sequentially adds into the model the predictor that most improves the fit. Clever updating algorithms can exploit the QR decomposition for the current fit to rapidly establish the next candidate. Like best-subset regression, the subset size \(k\) must be determined.

Forward-stepwise selection is a greedy algorithm, producing a nested sequence of models. In this sense it might seem sub-optimal compared to best-subset selection. However, there are several reasons why it might be preferred:

• Computational; we can always compute the forward stepwise sequence (even when \(p \gg N\)).
• Statistical;forward stepwise is a more constrained search than the best subset selection, and will have lower variance, but perhaps more bias.

Backward-stepwise selction

Backward-stepwise selection starts with the full model, and sequentially deletes the predictor that has the least impact on the fit. The candidate for dropping is the variable with the smallest Z-score. Backward selection can only be used when \(N>p\), while forward selection can always be used.

On the prostate cancer example, best-subset, forward and backward selection all gave exactly the same sequence of terms.

# there are 19 predictors
X.columns

Index(['AtBat', 'Hits', 'HmRun', 'Runs', 'RBI', 'Walks', 'Years', 'CAtBat',
       'CHits', 'CHmRun', 'CRuns', 'CRBI', 'CWalks', 'PutOuts', 'Assists',
       'Errors', 'League_N', 'Division_W', 'NewLeague_N'],
      dtype='object')

[p for p in X.columns if p not in []]

['AtBat',
 'Hits',
 'HmRun',
 'Runs',
 'RBI',
 'Walks',
 'Years',
 'CAtBat',
 'CHits',
 'CHmRun',
 'CRuns',
 'CRBI',
 'CWalks',
 'PutOuts',
 'Assists',
 'Errors',
 'League_N',
 'Division_W',
 'NewLeague_N']

def forward(predictors):

    # Pull out predictors we still need to process
    remaining_predictors = [p for p in X.columns if p not in predictors]
    
    tic = time.time()
    
    results = []
    
    for p in remaining_predictors:
        results.append(processSubset(predictors+[p]))
    
    # Wrap everything up in a nice dataframe
    models = pd.DataFrame(results)
    
    # Choose the model with the highest RSS
    best_model = models.loc[models['RSS'].argmin()]
    
    toc = time.time()
    print("Processed ", models.shape[0], "models on", len(predictors)+1, "predictors in", (toc-tic), "seconds.")
    
    # Return the best model, along with some other useful information about the model
    return best_model

models_fwd = pd.DataFrame(columns=["RSS", "model"])

tic = time.time()
predictors = []

for i in range(1,len(X.columns)+1):    
    models_fwd.loc[i] = forward(predictors)
    predictors = models_fwd.loc[i]["model"].model.exog_names

toc = time.time()
print("Total elapsed time:", (toc-tic), "seconds.")

Processed  19 models on 1 predictors in 0.03515005111694336 seconds.
Processed  18 models on 2 predictors in 0.033811330795288086 seconds.
Processed  17 models on 3 predictors in 0.02893376350402832 seconds.
Processed  16 models on 4 predictors in 0.029803991317749023 seconds.
Processed  15 models on 5 predictors in 0.028636455535888672 seconds.
Processed  14 models on 6 predictors in 0.027070999145507812 seconds.
Processed  13 models on 7 predictors in 0.025038480758666992 seconds.
Processed  12 models on 8 predictors in 0.025738239288330078 seconds.
Processed  11 models on 9 predictors in 0.022580862045288086 seconds.
Processed  10 models on 10 predictors in 0.02110886573791504 seconds.
Processed  9 models on 11 predictors in 0.01943826675415039 seconds.
Processed  8 models on 12 predictors in 0.016745567321777344 seconds.
Processed  7 models on 13 predictors in 0.014856338500976562 seconds.
Processed  6 models on 14 predictors in 0.013204336166381836 seconds.
Processed  5 models on 15 predictors in 0.011202096939086914 seconds.
Processed  4 models on 16 predictors in 0.009953975677490234 seconds.
Processed  3 models on 17 predictors in 0.007590293884277344 seconds.
Processed  2 models on 18 predictors in 0.005232572555541992 seconds.
Processed  1 models on 19 predictors in 0.0031538009643554688 seconds.
Total elapsed time: 0.427997350692749 seconds.

print(models_fwd.loc[1, "model"].summary())
print(models_fwd.loc[2, "model"].summary())

                                 OLS Regression Results                                
=======================================================================================
Dep. Variable:                 Salary   R-squared (uncentered):                   0.665
Model:                            OLS   Adj. R-squared (uncentered):              0.663
Method:                 Least Squares   F-statistic:                              519.2
Date:                Tue, 16 Feb 2021   Prob (F-statistic):                    4.20e-64
Time:                        16:41:14   Log-Likelihood:                         -1952.4
No. Observations:                 263   AIC:                                      3907.
Df Residuals:                     262   BIC:                                      3910.
Df Model:                           1                                                  
Covariance Type:            nonrobust                                                  
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Hits           4.8833      0.214     22.787      0.000       4.461       5.305
==============================================================================
Omnibus:                       90.075   Durbin-Watson:                   1.949
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              293.080
Skew:                           1.469   Prob(JB):                     2.28e-64
Kurtosis:                       7.256   Cond. No.                         1.00
==============================================================================

Notes:
[1] R² is computed without centering (uncentered) since the model does not contain a constant.
[2] Standard Errors assume that the covariance matrix of the errors is correctly specified.
                                 OLS Regression Results                                
=======================================================================================
Dep. Variable:                 Salary   R-squared (uncentered):                   0.761
Model:                            OLS   Adj. R-squared (uncentered):              0.760
Method:                 Least Squares   F-statistic:                              416.7
Date:                Tue, 16 Feb 2021   Prob (F-statistic):                    5.80e-82
Time:                        16:41:14   Log-Likelihood:                         -1907.6
No. Observations:                 263   AIC:                                      3819.
Df Residuals:                     261   BIC:                                      3826.
Df Model:                           2                                                  
Covariance Type:            nonrobust                                                  
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Hits           2.9538      0.261     11.335      0.000       2.441       3.467
CRBI           0.6788      0.066     10.295      0.000       0.549       0.809
==============================================================================
Omnibus:                      117.551   Durbin-Watson:                   1.933
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              654.612
Skew:                           1.729   Prob(JB):                    7.12e-143
Kurtosis:                       9.912   Cond. No.                         5.88
==============================================================================

Notes:
[1] R² is computed without centering (uncentered) since the model does not contain a constant.
[2] Standard Errors assume that the covariance matrix of the errors is correctly specified.

#Let's see how the models stack up against best subset selection:
print(models_best.loc[6, "model"].summary())
print(models_fwd.loc[6, "model"].summary())

                                 OLS Regression Results                                
=======================================================================================
Dep. Variable:                 Salary   R-squared (uncentered):                   0.795
Model:                            OLS   Adj. R-squared (uncentered):              0.790
Method:                 Least Squares   F-statistic:                              166.3
Date:                Tue, 16 Feb 2021   Prob (F-statistic):                    1.79e-85
Time:                        16:42:03   Log-Likelihood:                         -1887.6
No. Observations:                 263   AIC:                                      3787.
Df Residuals:                     257   BIC:                                      3809.
Df Model:                           6                                                  
Covariance Type:            nonrobust                                                  
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
AtBat         -1.5488      0.477     -3.248      0.001      -2.488      -0.610
Hits           7.0190      1.613      4.352      0.000       3.843      10.195
Walks          3.7513      1.212      3.095      0.002       1.364       6.138
CRBI           0.6544      0.064     10.218      0.000       0.528       0.781
PutOuts        0.2703      0.075      3.614      0.000       0.123       0.418
Division_W  -104.4513     37.661     -2.773      0.006    -178.615     -30.287
==============================================================================
Omnibus:                      106.414   Durbin-Watson:                   1.986
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              768.429
Skew:                           1.433   Prob(JB):                    1.37e-167
Kurtosis:                      10.869   Cond. No.                     1.29e+03
==============================================================================

Notes:
[1] R² is computed without centering (uncentered) since the model does not contain a constant.
[2] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[3] The condition number is large, 1.29e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
                                 OLS Regression Results                                
=======================================================================================
Dep. Variable:                 Salary   R-squared (uncentered):                   0.795
Model:                            OLS   Adj. R-squared (uncentered):              0.790
Method:                 Least Squares   F-statistic:                              166.3
Date:                Tue, 16 Feb 2021   Prob (F-statistic):                    1.79e-85
Time:                        16:42:03   Log-Likelihood:                         -1887.6
No. Observations:                 263   AIC:                                      3787.
Df Residuals:                     257   BIC:                                      3809.
Df Model:                           6                                                  
Covariance Type:            nonrobust                                                  
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Hits           7.0190      1.613      4.352      0.000       3.843      10.195
CRBI           0.6544      0.064     10.218      0.000       0.528       0.781
Division_W  -104.4513     37.661     -2.773      0.006    -178.615     -30.287
PutOuts        0.2703      0.075      3.614      0.000       0.123       0.418
AtBat         -1.5488      0.477     -3.248      0.001      -2.488      -0.610
Walks          3.7513      1.212      3.095      0.002       1.364       6.138
==============================================================================
Omnibus:                      106.414   Durbin-Watson:                   1.986
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              768.429
Skew:                           1.433   Prob(JB):                    1.37e-167
Kurtosis:                      10.869   Cond. No.                     1.29e+03
==============================================================================

Notes:
[1] R² is computed without centering (uncentered) since the model does not contain a constant.
[2] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[3] The condition number is large, 1.29e+03. This might indicate that there are
strong multicollinearity or other numerical problems.

For this data, the best one-variable through six-variable models are each identical for best subset and forward selection.

# Backward Selection
def backward(predictors):
    
    tic = time.time()
    
    results = []
    
    for combo in itertools.combinations(predictors, len(predictors)-1):
        results.append(processSubset(combo))
    
    # Wrap everything up in a nice dataframe
    models = pd.DataFrame(results)
    
    # Choose the model with the highest RSS
    best_model = models.loc[models['RSS'].argmin()]
    
    toc = time.time()
    print("Processed ", models.shape[0], "models on", len(predictors)-1, "predictors in", (toc-tic), "seconds.")
    
    # Return the best model, along with some other useful information about the model
    return best_model

models_bwd = pd.DataFrame(columns=["RSS", "model"], index = range(1,len(X.columns)))

tic = time.time()
predictors = X.columns

while(len(predictors) > 1):  
    models_bwd.loc[len(predictors)-1] = backward(predictors)
    predictors = models_bwd.loc[len(predictors)-1]["model"].model.exog_names

toc = time.time()
print("Total elapsed time:", (toc-tic), "seconds.")

Processed  19 models on 18 predictors in 0.05166292190551758 seconds.
Processed  18 models on 17 predictors in 0.03813934326171875 seconds.
Processed  17 models on 16 predictors in 0.034967660903930664 seconds.
Processed  16 models on 15 predictors in 0.03255939483642578 seconds.
Processed  15 models on 14 predictors in 0.03638648986816406 seconds.
Processed  14 models on 13 predictors in 0.0325469970703125 seconds.
Processed  13 models on 12 predictors in 0.0324857234954834 seconds.
Processed  12 models on 11 predictors in 0.02887129783630371 seconds.
Processed  11 models on 10 predictors in 0.026468515396118164 seconds.
Processed  10 models on 9 predictors in 0.02369976043701172 seconds.
Processed  9 models on 8 predictors in 0.021591901779174805 seconds.
Processed  8 models on 7 predictors in 0.020296335220336914 seconds.
Processed  7 models on 6 predictors in 0.017458677291870117 seconds.
Processed  6 models on 5 predictors in 0.017169475555419922 seconds.
Processed  5 models on 4 predictors in 0.014002323150634766 seconds.
Processed  4 models on 3 predictors in 0.01012420654296875 seconds.
Processed  3 models on 2 predictors in 0.005807399749755859 seconds.
Processed  2 models on 1 predictors in 0.0037310123443603516 seconds.
Total elapsed time: 0.4661753177642822 seconds.

For this data, the best one-variable through six-variable models are each identical for best subset and forward selection. However, the best seven-variable models identified by forward stepwise selection, backward stepwise selection, and best subset selection are different:

print("------------")
print("Best Subset:")
print("------------")
print(models_best.loc[7, "model"].params)

------------
Best Subset:
------------
Hits            1.680029
Walks           3.399961
CAtBat         -0.328835
CHits           1.347017
CHmRun          1.349373
PutOuts         0.248166
Division_W   -111.943760
dtype: float64

print("-----------------")
print("Foward Selection:")
print("-----------------")
print(models_fwd.loc[7, "model"].params)

-----------------
Foward Selection:
-----------------
Hits            7.277149
CRBI            0.652415
Division_W   -110.656338
PutOuts         0.259787
AtBat          -1.644651
Walks           3.684324
League_N       49.978410
dtype: float64

print("-------------------")
print("Backward Selection:")
print("-------------------")
print(models_bwd.loc[7, "model"].params)

-------------------
Backward Selection:
-------------------
AtBat         -1.601655
Hits           6.148449
Walks          5.866033
CRuns          1.097453
CWalks        -0.650614
PutOuts        0.310125
Division_W   -95.027171
dtype: float64

Cross-validation, briefly

1. Cross-validation works by dividing the training data randomly into ten equal parts.
2. The learning method is fit – for a range of values of the complexity parameter – to nine-tenths of the data.
3. The prediction error is computed on the remaining one-tenth.
4. Repeat step 2 - step 3 for each one-tenth of the data, and the ten prediction error estimates are averaged.
5. Then we obtain an estimated prediction error curve as a function of the complexity parameter so that a proper complexity parameter can be chosen.

We have used the “one-standard-error” rule – we pick the most parsimonious model within one standard error of the minimun. Such a rule acknowledges the fact that the tradeoff curve is estimated with error, and hence takes a conservative approach.

Note that we have already divided these data into a training set of size 67 and a test set of size 30. Cross-validation is applied to the training set, since selecting the shrinkage parameter is part of the training process. The test set is there to judge the performance of the selected model.

np.arange(67)

array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,
       17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33,
       34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50,
       51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66])

divmod(67, 10)

(6, 7)

div, mod = divmod(67, 10)
unit_sizes = [div for _ in range(10)]
for i in range(mod):
    unit_sizes[i] += 1
print(unit_sizes)
sum(unit_sizes)

[7, 7, 7, 7, 7, 7, 7, 6, 6, 6]





67

for k, unit_size in enumerate(unit_sizes):
    print(k, unit_size)

0 7
1 7
2 7
3 7
4 7
5 7
6 7
7 6
8 6
9 6

np.random.choice(np.arange(67), unit_size, replace=False)

array([35,  7, 64, 58,  4, 63])

"""FIGURE 3.7. Estimated prediction error curves and their standard errors
for the best-subset selection via tenfold cross-validation."""

def index_tenfold(n:int) -> np.ndarray:
    """Produce index array for tenfold CV with dataframe length n."""
    original_indices = np.arange(n)
    tenfold_indices = np.zeros(n)

    div, mod = divmod(n, 10)
    unit_sizes = [div for _ in range(10)]
    for i in range(mod):
        unit_sizes[i] += 1

    for k, unit_size in enumerate(unit_sizes):
        tenfold = np.random.choice(original_indices, unit_size,
                                   replace=False)
        tenfold_indices[tenfold] = k
        original_indices = np.delete(
            original_indices,
            [np.argwhere(original_indices == val) for val in tenfold],
        )
        # print(tenfold, original_indices)
    return tenfold_indices

print(index_tenfold(67))
len(index_tenfold(67))

[9. 7. 9. 8. 5. 1. 7. 5. 8. 3. 2. 4. 7. 4. 9. 9. 0. 0. 3. 4. 3. 4. 3. 0.
 1. 1. 8. 7. 1. 2. 3. 5. 6. 8. 2. 9. 2. 8. 2. 0. 2. 6. 5. 4. 5. 7. 4. 1.
 0. 1. 6. 6. 6. 5. 9. 6. 5. 0. 0. 2. 6. 4. 3. 7. 3. 8. 1.]





67

data_x_train = data_x_normalized[data['train'] == 'T']
data_y_train = data_y[data['train'] == 'T']
data_x_test = data_x_normalized[data['train'] == 'F']
data_y_test = data_y[data['train'] == 'F']
vec_y = data_y_train.values
vec_y_test = data_y_test.values

size_train = sum(data['train'] == 'T')
size_test = sum(data['train'] == 'F')
size_predictor = len(data_x_train.columns)

indices_tenfold = index_tenfold(size_train)
cv_results = collections.defaultdict(list)

for cv in range(10):
    mask_cv = indices_tenfold != cv
    size_cv_train = sum(mask_cv == True)

    df_x = data_x_train[mask_cv]
    v_y = vec_y[mask_cv]
    #All Subsets OLS
    for k in range(size_predictor+1):
        ols = ols_with_subset_size(df_x, v_y, k)
        ols_best = min(ols, key=op.itemgetter('rss'))
        cv_results[k].append(ols_best['rss']/size_cv_train)

rss_average = [sum(rss)/10 for _, rss in cv_results.items()]
rss_stderr = [math.sqrt((sum(np.array(rss)**2)-10*rss_average[k]**2)/9)
              for k, rss in cv_results.items()]

rss_average

[1.435705040941674,
 0.6627739101582181,
 0.55114165895448,
 0.5119002341038311,
 0.48542751353551605,
 0.4682757004380743,
 0.44811692523738395,
 0.43369589944674913,
 0.4329883352963549]

list(cv_results.keys())

[0, 1, 2, 3, 4, 5, 6, 7, 8]

fig37 = plt.figure(figsize=(10, 10))
ax = fig37.add_subplot(1, 1, 1)
ax.plot(list(cv_results.keys()), rss_average, 'o-', color='C1')
for k, (ave, stderr) in enumerate(zip(rss_average, rss_stderr)):
    ax.plot([k, k], [ave-stderr, ave+stderr], color='C0', linewidth=1)
    ax.plot([k-.1, k+.1], [ave-stderr, ave-stderr], color='C0', linewidth=1)
    ax.plot([k-.1, k+.1], [ave+stderr, ave+stderr], color='C0', linewidth=1)
ax.set_xlabel('Subset Size k')
ax.set_ylabel('RSS_average of cv')
plt.show()

png

FIGURE 3.7. Estimated prediction error curves and their standard errors for the various selection and shrinkage methods. Each curve is plotted as a function of the corresponding complexity parameter for that method. The horizontal axis has been chosen so that the model complexity increases as we move from left to right. The estimates of prediction error and their standard errors were obtained by tenfold cross-validation; full details are given in Section 7.10.

ols_k2 = ols_with_subset_size(data_x_train, vec_y, k=2)
ols_k2_best = min(ols_k2, key=op.itemgetter('rss'))
[list(ols_k2_best['column_names'])]

[['lcavol', 'lweight']]

ols_k2_best['beta']

array([2.47735734, 0.73589083, 0.31469341])

"""Table 3.3. Estimated coefficients ad test error result for the
best-subset selection applied to the prostate data"""

ols_k2 = ols_with_subset_size(data_x_train, vec_y, k=2)
ols_k2_best = min(ols_k2, key=op.itemgetter('rss'))

df_x_test = data_x_test[list(ols_k2_best['column_names'])]
mat_x_test = np.hstack((np.ones((size_test, 1)), df_x_test.values))
vec_y_test_fitted = mat_x_test @ ols_k2_best['beta']

print('{0:>15} {1:>15}'.format('Term', 'Best Subset'))
print('-'*31)
print('{0:>15} {1:>15.3f}'.format('Intercept', ols_k2_best['beta'][0]))
for idx, col_name in enumerate(ols_k2_best['column_names']):
    print('{0:>15} {1:>15.3f}'.format(col_name, ols_k2_best['beta'][idx+1]))
print('-'*31)
print('{0:>15} {1:>15.3f}'.format(
    'Test Error',
    sum((vec_y_test-vec_y_test_fitted)**2)/size_test))
print('{0:>15} {1:>15.3f}'.format(
    'Std Error',
    math.sqrt((sum(np.array(vec_y_test-vec_y_test_fitted)**2)-size_test*(np.mean(vec_y_test-vec_y_test_fitted))**2)/(size_test-1)),
))

           Term     Best Subset
-------------------------------
      Intercept           2.477
         lcavol           0.736
        lweight           0.315
-------------------------------
     Test Error           0.492
      Std Error           0.714

Scaling and centering

The ridge solutions are not equivariant under scaling of the inputs, and so one normally standardizes the inputs before solving the ridge problem.

Also notice that the intercept \(\beta_0\) has been left out of the penalty term. Penalization of the intercept would make the procedure depend on the origin chosen for \(Y\); i.e., adding a constant \(c\) to each of the targets \(y_i\) (i.e. simply shifting) would not simply result in a shift of the predictions by the same constant \(c\).

It can be shown that the ridge solution can be separated into two parts, after reparametrization using centered inputs: Each \(x_{ij}\) gets replaced by \(x_{ij}-\bar{x}_j\). 1. We estimate \(\beta_0\) by the mean response \(\bar{y} = \frac{1}{N}\sum_1^N y_i\). 2. The remaining coefficients get estimated by a ridge regression without intercept, using centered \(x_{ij}\).

Henceforth we assume that this centering has been done, so that the input matrix \(\mathbf{X}\) has \(p\) columns rather than \(p+1\).

From the Bayesian point of view

Ridge regression can also be derived as the mean or mode of a posterior distribution, with a suitably chosen prior distribution.

Suppose

\[\begin{align} y_i &\sim N(\beta_0+x_i^T\beta, \sigma^2) \\ \beta_j &\sim \text{ i.i.d. }N(0, \tau^2) \end{align}\]

Then the log-posterior density of \(\beta\), with \(\tau^2\) and \(\sigma^2\) assumed known, is equal to the expression

\[\begin{equation} \sum_{i=1}^N\left( y_i - \beta_0 - \sum_{j=1}^p x_{ij}\beta_j \right)^2 +\frac{\sigma^2}{\tau^2}\sum_{j=1}^p \beta_j^2, \end{equation}\]

which is the penalized residual sum of squares with \(\lambda = \sigma^2/\tau^2\). Thus the ridge estimate is the mode of the posterior distribution; since the distribution is Gaussian, it is also the posterior mean.

Matrix form

We write the criterion:

\[\begin{equation} \hat\beta^{\text{ridge}} = {\arg\min}_{\beta}\left\lbrace \sum_{i=1}^N\left( y_i - \beta_0 - \sum_{j=1}^p x_{ij}\beta_j \right)^2 +\lambda\sum_{j=1}^p \beta_j^2 \right\rbrace, \end{equation}\]

in matrix form: \[\begin{equation} \text{RSS}(\lambda) = (\mathbf{y}-\mathbf{X}\beta)^T(\mathbf{y}-\mathbf{X}\beta) + \lambda\beta^T\beta \end{equation}\]

This is a quadratic function in the \(p+1\) parameters. Differentiating w.r.t. \(\beta\) we obtain

\[\begin{align} \frac{\partial\text{RSS}(\lambda)}{\partial\beta} &= -2\mathbf{X}^T\left(\mathbf{y}-\mathbf{X}\beta\right) + 2\lambda\beta \\ \frac{\partial^2\text{RSS}(\lambda)}{\partial\beta\partial\beta^T} &= 2\mathbf{X}^T\mathbf{X} + 2\lambda \end{align}\]

we set the first derivative to zero

\[\begin{equation} \mathbf{X}^T(\mathbf{y}-\mathbf{X}\beta) - \lambda\beta = \mathbf{X}^T\mathbf{y} - (\mathbf{X}^T\mathbf{X}+\lambda\mathbf{I})\beta = 0 \end{equation}\]

to obtain the unique solution

\[\begin{equation} \hat\beta^{\text{ridge}} = \left( \mathbf{X}^T\mathbf{X} + \lambda\mathbf{I} \right)^{-1}\mathbf{X}^T\mathbf{y}, \end{equation}\]

where \(\mathbf{I}\) is the \(p\times p\) identity matrix. Notice that the ridge solution is again a linear function of \(\mathbf{y}\) by the choice of quadratic penalty \(\beta^T\beta\), resulting in addition of a positive constant to the diagonal of \(\mathbf{X}^T\mathbf{X}\) before inversion. This makes the problem nonsingular, even if \(\mathbf{X}^T\mathbf{X}\) is not of full rank. It was actually the main motivation for ridge regression when it was first introduced in statistics (Hoerl and Kennard, 1970).

The singular value decomposition (SVD)

The SVD of the centered input matrix \(\mathbf{X}\) gives us some additional insight into the nature of ridge regression. The SVD of the \(N\times p\) matrix \(\mathbf{X}\) has the form

\[\begin{equation} \mathbf{X} = \mathbf{UDV}^T, \end{equation}\]

where

• \(\mathbf{U}\) is \(N\times p\) orthogonal matrix, with the columns of \(\mathbf{U}\) spanning the \(\text{col}(\mathbf{X})\)
• \(\mathbf{V}\) is \(p\times p\) orthogonal matrix, with the columns of \(\mathbf{V}\) spanning the \(\text{row}(\mathbf{X})\)
• \(\mathbf{D}\) is a \(p\times p\) diagonal matrix, with diagonal entries
  \(d_1 \ge d_2 \ge \cdots \ge d_p \ge 0\) called the singular values of \(\mathbf{X}\).
• If one or more values \(d_j = 0\), \(\mathbf{X}\) is singular.

If we compute the singular value decomposition (SVD) of the \(N \times p\) centered data matrix \(\mathbf X\) as \[\mathbf X=\mathbf U\mathbf D\mathbf V^T\] where \(\mathbf U\) is a \(N × p\) matrix with orthonormal columns that span the column space of \(\mathbf X\), \(\mathbf V\) is a \(p\times p\) orthogonal matrix, and \(\mathbf D\) is a \(p\times p\) diagonal matrix with elements \(d_j\) ordered such that \(d_1 \ge d_2 \ge \cdots d_p \ge 0\). From this representation of \(\mathbf X\) we can derive a simple expression for \(\mathbf X^T\mathbf X\): \[\mathbf X^T\mathbf X=\mathbf V\mathbf D\mathbf U^T\mathbf U\mathbf D\mathbf V^T=\mathbf V\mathbf D^2\mathbf V^T\] Using this expression we can compute the least squares fitted values \(\hat{\mathbf y}^{ls} = \mathbf X \hat{\beta}^{ls}\) as \[\begin{align} \hat{y}^{ls} &= \mathbf X \hat{\beta}^{ls}\\ &=\mathbf X(\mathbf X^T\mathbf X)^{-1}\mathbf X^T\mathbf y\\ &=\mathbf U\mathbf D\mathbf V^T(\mathbf V\mathbf D^2\mathbf V^T)^{-1}\mathbf V\mathbf D\mathbf U^T\mathbf y\\ &=\mathbf U\mathbf D\mathbf V^T(\mathbf V^{-T}\mathbf D^{-2}\mathbf V^{-1})\mathbf V\mathbf D\mathbf U^T\mathbf y\\ &=\mathbf U\mathbf U^T\mathbf y\\ &=\sum_{j=1}^{p}u_ju_j^T\mathbf y \end{align}\]

Note

• \(\mathbf{U}^T\mathbf{y}\) are the coordinates of \(\mathbf{y}\) w.r.t. the orthonormal basis \(\mathbf{U}\).
• the similarity with QR decomposition; \(\mathbf{Q}\) and \(\mathbf{U}\) are generally different orthogonal bases for \(\text{col}(\mathbf{X})\).

To compare how the fitted values \(\hat{\mathbf y}^{ls}\) obtained in ridge regression compare with ordinary least squares, we next consider the SVD expression for \(\hat{\beta}^{ridge}\). In the same way as for least squares we find: \[ \begin{align} \hat\beta^{\text{ridge}} &= \left( \mathbf{X}^T\mathbf{X} + \lambda\mathbf{I} \right)^{-1}\mathbf{X}^T\mathbf{y}\\ &=(\mathbf V\mathbf D^2\mathbf V^T+\lambda\mathbf{V}\mathbf{V}^T)^{-1}\mathbf V\mathbf D\mathbf U^T\mathbf{y}\\ &=(\mathbf V(\mathbf D^2+\lambda\mathbf{I})\mathbf{V}^T)^{-1}\mathbf V\mathbf D\mathbf U^T\mathbf{y}\\ &=\mathbf{V}^{-T}(\mathbf D^2+\lambda\mathbf{I})^{-1}\mathbf D\mathbf U^T\mathbf{y}\\ &=\mathbf{V}(\mathbf D^2+\lambda\mathbf{I})^{-1}\mathbf D\mathbf U^T\mathbf{y}\\ \end{align} \] And \[ \begin{align} \hat y^{\text{ridge}} &=\mathbf X\hat\beta^{\text{ridge}}\\ &=\mathbf U\mathbf D\mathbf V^T\mathbf{V}(\mathbf D^2+\lambda\mathbf{I})^{-1}\mathbf D\mathbf U^T\mathbf{y}\\ &=\mathbf U\mathbf D(\mathbf D^2+\lambda\mathbf{I})^{-1}\mathbf D\mathbf U^T\mathbf{y}\\ \end{align} \] Now note that in this last expression \(\mathbf D(\mathbf D^2+\lambda\mathbf{I})^{-1}\mathbf D\) is a diagonal matrix with elements given by \[\frac{d_j^2}{d_j^2+\lambda}\] and the vector \(\mathbf U^T\mathbf y\) is the coordinates of the vector \(\mathbf y\) in the basis spanned by the \(p\)-columns of \(\mathbf U\). Thus writing the expression by summing columns we obtain: \[\begin{align} \hat y^{\text{ridge}} &= \mathbf X\hat\beta^{\text{ridge}}\\ &=\mathbf U\mathbf D(\mathbf D^2+\lambda\mathbf{I})^{-1}\mathbf D\mathbf U^T\mathbf{y}\\ &=\sum_{j=1}^{p}u_j\left(\frac{d_j^2}{d_j^2+\lambda}\right)u_j^T\mathbf y\\ \end{align}\]

Note that this result is similar to that found in Equation \[\sum_{j=1}^{p}u_ju_j^T\mathbf y\] derived for ordinary least squares regression but in ridge-regression the inner products \(u_j^T\mathbf y\) are now scaled by the factors \(\left(\frac{d_j^2}{d_j^2+\lambda}\right)\). This means that a greater amount of shrinkage is applied to the coordinates of basis vectors with smaller \(d_j^2\). Then what does a small value of \(d_j^2\) mean?

The SVD and the principal components

The SVD of the centered matrix \(\mathbf{X}\) is another way of expressing the principal components of the variables in \(\mathbf{X}\). The sample covariance matrix is given by

\[\begin{equation} \mathbf{S} = \frac{1}{N}\mathbf{X}^T\mathbf{X}, \end{equation}\]

and via the SVD,

\[\begin{equation} \mathbf{X}^T\mathbf{X} = \mathbf{VD}^2\mathbf{V}^T, \end{equation}\]

which is the eigen decomposition of \(\mathbf{X}^T\mathbf{X}\) (and of \(\mathbf{S}\), up to a factor \(N\)). The eigenvectors \(v_j\) (columns of \(\mathbf{V}\)) are also called the principal components (or Karhunen-Loeve) directions of \(\mathbf{X}\) (in the \(p\)-dimension space). The first principal component direction \(\mathbf v_1\) has the property that \(\mathbf{z}_1 = \mathbf{X}\mathbf v_1\) has the larger sample projected variance than any other direction. This sample variance is easily seen to be

\[\begin{equation} \mathbf{X}\mathbf{V}=\mathbf{U}\mathbf{D} \end{equation}\]

\[\begin{equation} \text{Var}(\mathbf{z}_1) = \text{Var}(\mathbf{X}\mathbf v_1) = \text{Var}(\mathbf{u_1}d_1)=d_1^2\text{Var}(\mathbf{u_1})=\frac{d_1^2}{N}(\mathbf{u_1}^T\mathbf{u_1})= \frac{d_1^2}{N}, \end{equation}\]

and in fact \(\mathbf{z}_1 = \mathbf{X}\mathbf v_1 = \mathbf{u}_1 d_1\). The derived variable \(z_1\) is called the first principal component of \(\mathbf{X}\), and hence \(\mathbf{u}_1\) is the normalized first principal component (in the \(N\)-dimension space). And subsequent principal components \(z_j\) have maximum variance \(d_j^2/N\), subject to being orthogonal to the earlier ones. Conversely the last principal component has minimum variance. Hence the small singular values \(d_j\) correspond to directions in the \(\text{col}(\mathbf{X})\) having small variance, and ridge regression shrinks these directions the most. In simple words, it does not care what seems not worth. See FIGURE 3.9 for the graphical representation of the principal components.

The 1st principal component compresses \(N\) points in \(p\)-dimention into \(N\) points in \(1\)-dimention which maintains the largest variance. The 2nd principal component compresses \(N\) points in \(p\)-dimention into \(N\) points in \(1\)-dimention which is perpendicular to the 1st principal component and maintains the largest remained variance. And so on with the \(N^{st}\) principal component.

Ridge regression protects against the potentially high variance of gradients estimated in the short directions. The implicit assumption is that the response will tend to vary most in the directions of high variance of the inputs. This is often a reasonable assumption, since predictors are often chosen for study because they vary with the response variable, but need not hold in general.

The effective degrees of freedom

\[\hat y^{\text{ridge}} =\mathbf X\hat\beta^{\text{ridge}}=\mathbf X\left( \mathbf{X}^T\mathbf{X} + \lambda\mathbf{I} \right)^{-1}\mathbf{X}^T\mathbf{y}\] In FIGURE 3.7 we have plotted the estimated prediction error versus the quantity. The definition of the effective degrees of freedom \(df(\lambda)\) in ridge regression is given by \[\begin{align} \text{df}(\lambda) &= \text{tr}\left( \mathbf{X}(\mathbf{X}^T\mathbf{X} + \lambda\mathbf{I})^{-1}\mathbf{X}^T \right)\\ &= \text{tr}(\mathbf{H}_\lambda) \\ &= \text{tr}\left(\mathbf U\mathbf D(\mathbf D^2+\lambda\mathbf{I})^{-1}\mathbf D\mathbf U^T\right)\\ &=\text{tr}\left(\sum_{j=1}^{p}u_j\left(\frac{d_j^2}{d_j^2+\lambda}\right)u_j^T\right)\\ &= \sum_{j=1}^p \frac{d_j^2}{d_j^2+\lambda}. \end{align}\] Since the trace of a matrix can be shown to equal the sum of its eigenvalues.

This monotone decreasing function of \(\lambda\) is the effective degrees of freedom of the ridge regression fit. Usually in a linear-regression fit, the degrees-of-freedom of the fit is \(p\), the number of free parameters. The idea is that although all \(p\) coefficients in a ridge fit will be non-zero, they are fit in a restricted fashion controlled by \(\lambda\). Note that

\[\begin{align} \text{df}(\lambda) &= p \text{ when }\lambda = 0, \\ \text{df}(\lambda) &\rightarrow 0 \text{ as }\lambda \rightarrow \infty. \end{align}\]

Of course there is always an additional one degree of freedom for the intercept, which was removed apriori. This definition is motivated in more detail in Section 3.4.4 and Sections 7.4-7.6.

One important consequence of this expression is that we can use it to determine the values of \(\lambda\) for which to use when applying cross validation. For example, the book discusses how to obtain the estimate of y when using ridge regression and it is given by Equation \[\hat y^{\text{ridge}} =\mathbf X\hat\beta^{\text{ridge}}=\mathbf X\left( \mathbf{X}^T\mathbf{X} + \lambda\mathbf{I} \right)^{-1}\mathbf{X}^T\mathbf{y}\] but no mention of the numerical values of \(\lambda\) we should use in this expression to guarantee that we have accurate coverage of all possible regularized linear models. The approach taken in generating the ridge regression results in Figure 3.8 is to consider \(df\) in Equation \[\text{df}(\lambda)=\sum_{j=1}^p \frac{d_j^2}{d_j^2+\lambda}\] a function of \(\lambda\). As such we set \(df(\lambda) = k\) for \(k = 1, 2, \cdots , p\) representing all of the possible values for the degree of freedom. We then use Newton’s root finding method to solve for \(\lambda\) in the expression \[\text{df}(\lambda)=\sum_{j=1}^p \frac{d_j^2}{d_j^2+\lambda}=k, (k = 1, 2, \cdots , p)\]

Newton’s method (Newton-Raphson method) uses tangent lines of the graph of \(y = ƒ(x)\) near the points where \(ƒ(x)=0\) to estimate the solution. The goal of Newton’s method for estimating a solution of an equation \(ƒ(x) = 0\) is to produce a sequence of approximations that approach the solution. The initial estimate, \(x_0\), may be found by just plain guessing. The method then uses the tangent to the curve \(y = ƒ(x)\) at \((x_0, ƒ(x_0))\) to approximate the curve, calling the point \(x_1\) where the tangent meets the \(x\)-axis. The number \(x_1\) is usually a better approximation to the solution than is \(x_0\). The point \(x_2\) where the tangent to the curve at \((x_1, ƒ(x_1))\) crosses the \(x\)-axis is the next approximation in the sequence. We continue on, using each approximation to generate the next, until we are close enough to the root to stop.
Given the approximation \(x_n\), the point-slope equation for the tangent to the curve at \((x_n,ƒ(x_n))\) is \[y = ƒ(x_n) + ƒ'(x_n)(x - x_n)\] We can find where it crosses the \(x\)-axis by setting \(y=0\), \[0 = ƒ(x_n) + ƒ'(x_n)(x - x_n)\] \[x - x_n=-\frac{ƒ(x_n)}{ƒ'(x_n)}\] \[x=x_n-\frac{ƒ(x_n)}{ƒ'(x_n)}\] This value of \(x\) is the next approximation \(x_{n+1}\). Then a solution of the equation \(f(x)=0\) is a fixed point of the function \[x \mapsto x-\frac{f(x)}{f'(x)}\] Newton’s method is the use of the fixed-point method to approximate a solution of the equation by finding a fixed point of the function \[x-\frac{f(x)}{f'(x)}\]

Thus we define a function \(d(\lambda)\) given by \[d(\lambda)=\sum_{j=1}^p \frac{d_j^2}{d_j^2+\lambda}-k\] and we want \(\lambda\) such that\(d(\lambda) = 0\). We use Newton’s algorithm for this where we iterate given a starting value of \(\lambda_0\) \[\lambda_{n+1}=\lambda_n-\frac{d(\lambda_n)}{d'(\lambda_n)}\] Thus we need the derivative of \(d(\lambda)\) which is given by \[d'(\lambda)=-\sum_{j=1}^{p}\frac{d_j^2}{(d_j^2+\lambda)^2}\] and an initial guess for \(\lambda_0\). Since we are really looking for \(p\) values of \(\lambda\) (one for each value of \(k\)) we will start by solving the problems for \(k = p, p−1, p−2, \cdots , 1\). When \(k = p\) the value of \(\lambda\) that solves \(df(\lambda) = p\) is seen to be \(\lambda = 0\). For each subsequent value of \(k\) we use the estimate of \(\lambda\) found in the previous Newton solve as the initial guess for the current Newton solve.

"""FIGURE 3.8. Profiles of ridge coefficients for the prostate cancer
example as tuning parameter lambda is varied.

Coefficients are plotted versus df(lambda), the effective degrees of
freedom. In case of orthonormal inputs, the ridge estimates are just a
scaled version of the least squares estimates;
beta_ridge = beta_ols/(1+lambda)
"""
import scipy
import scipy.stats
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

data = pd.read_csv('../../data/prostate/prostate.data', delimiter='\t',
                   index_col=0)
data_y = data.pop('lpsa')
mask_train = data.pop('train')

data_x_train = data[mask_train == 'T']
data_y_train = data_y[mask_train == 'T']
beta_intercept = data_y_train.mean()
# Centering for the training data y by subtracting the mean
data_y_train_centered = data_y_train.subtract(beta_intercept)
# Centering and scaling for the training data x using zscore
data_x_train_normalized = data_x_train.apply(scipy.stats.zscore)
vec_y = data_y_train_centered.values

data_x_test = data[mask_train == 'F']
# Centering and scaling for the test data x using zscore
data_x_test_normalized = data_x_test.apply(scipy.stats.zscore)
data_y_test = data_y[mask_train == 'F']
vec_y_test = data_y_test.values

size_train = sum(mask_train == 'T')
size_test = sum(mask_train == 'F')
size_predictor = data_x_train.columns.size

def lambdas_from_edf(singular_squared:np.ndarray, interval:int) ->np.ndarray:
    """Given squared singular values of data matrix, calculate the lambdas
    with `interval` parameter to split unit intervals s.t. the resulting
    effective degrees of freedom are equidistant with 1/interval, via the
    Newton-Raphson method. e.g., if interval = 10, it produces lambdas for
    0, 0.5, 0.6, 0.7, ..."""
    p = singular_squared.size
    edfs = np.linspace(.5, p-.5, (p-1)*interval+1)
    threshold = 1e-3
    lambdas = []
    for edf in edfs:
        # Newton-Raphson
        lambda0 = (p-edf)/edf
        lambda1 = 1e6
        diff = lambda1 - lambda0
        while diff > threshold:
            num = (singular_squared/(singular_squared+lambda0)).sum()-edf
            denom = (singular_squared/((singular_squared+lambda0)**2)).sum()
            lambda1 = lambda0 + num/denom
            diff = lambda1 - lambda0
            lambda0 = lambda1
        lambdas.append(lambda1)
    lambdas.append(0)
    edfs = np.concatenate(([0], edfs, [p]))
    return edfs, np.array(lambdas)

# singular value decomposition of data_x_train_normalized
u, s, vh = scipy.linalg.svd(data_x_train_normalized, full_matrices=False)
s2 = s**2
edfs, lambdas = lambdas_from_edf(s2, 10)
# beta_ols = vh.T @ scipy.diag(scipy.reciprocal(s)) @ u.T @ vec_y
# print(beta_ols)
beta_ridge_array = [np.zeros(size_predictor)]
for lamb in lambdas:
    mat_diag = np.diag(s/(s2+lamb))
    beta_ridge = vh.T @ mat_diag @ u.T @ vec_y
    beta_ridge_array.append(beta_ridge)
beta_ridge_array = np.array(beta_ridge_array)

fig38 = plt.figure(figsize=(10, 10))
ax = fig38.add_subplot(1, 1, 1)
ax.plot(edfs, beta_ridge_array, 'o-', markersize=2)
ax.legend(data_x_train.columns)
ax.set_xlabel(r'df ($\lambda$)')
ax.set_ylabel('Coefficients')
ax.set_xlim(-.2, 9)
ax.plot([-.2, 9], [0, 0], '--', color='0', linewidth=0.5)
ax.plot([5, 5], [-0.35, 0.75], '--', color='r', linewidth=0.5)
for x, y, s in zip(np.ones(8)*(8+0.1), beta_ridge_array[-1], data_x_train.columns):
    ax.text(x, y, s, color='0', fontsize=12)
ax.margins(0,0)
plt.show()

png

len(lambdas)

72

FIGURE 3.8. Profiles of ridge coefficients for the prostate cancer example, as the tuning parameter \(\lambda\) is varied. Coefficients are plotted versus df(\(\lambda\)), the effective degrees of freedom. A vertical line is drawn at df = 5.0, the value chosen by cross-validation.

import numpy as np
import pandas as pd
# 10-folds cross-validation

def index_tenfold(n:int) ->np.ndarray:
    """Produce index array for tenfold CV with dataframe length n."""
    original_indices = np.arange(n)
    tenfold_indices = np.zeros(n)

    div, mod = divmod(n, 10)
    unit_sizes = [div for _ in range(10)]
    for i in range(mod):
        unit_sizes[i] += 1

    for k, unit_size in enumerate(unit_sizes):
        tenfold = np.random.choice(original_indices, unit_size,
                                   replace=False)
        tenfold_indices[tenfold] = k
        original_indices = np.delete(
            original_indices,
            [np.argwhere(original_indices == val) for val in tenfold],
        )
        # print(tenfold, original_indices)
    return tenfold_indices

"""FIGURE 3.7. CV result of the ridge regression"""
import collections
import math

cv10_indices = index_tenfold(size_train)
cv_beta = collections.defaultdict(list)
cv_rss = collections.defaultdict(list)

for cv_idx in range(10):
    mask_cv = cv10_indices != cv_idx
    size_cv_train = (mask_cv == True).size
    # Use the zscore Centered and scaled training data x
    df_x = data_x_train_normalized[mask_cv]
    cv_x = df_x.values
    # Use the subtracting-mean Centered training data y
    v_y = vec_y[mask_cv]

    # Exceptional case for lambda = infinity, or edf = 0 (constant model)
    v_y_mean= v_y.mean()
    cv_rss[0].append(((v_y-v_y_mean)**2).sum()/size_cv_train)
    
    u, s, vh = np.linalg.svd(cv_x, full_matrices=False)
    s2 = s**2
    edfs, lambdas = lambdas_from_edf(s2, 2)

    for edf, lamb in zip(edfs[1:], lambdas):
        mat_diag = np.diag(s/(s2+lamb))
        beta_ridge = vh.T @ mat_diag @ u.T @ v_y
        cv_beta[edf].append(beta_ridge)
        cv_y_fitted = cv_x @ beta_ridge
        cv_rss[edf].append(((v_y-cv_y_fitted)**2).sum()/size_cv_train)

cv_rss_average = [sum(rss)/len(rss) for _, rss in cv_rss.items()]
cv_rss_stderr = [math.sqrt(((np.array(rss)**2).sum()-sum(rss)**2/10)/9)
                 for _, rss in cv_rss.items()]

# For the 17 edfs, each edfs has 10 cv_rss and 1 cv_rss_average
len(cv_rss_average)

17

fig27 = plt.figure(figsize=(10, 10))
ax = fig27.add_subplot(1, 1, 1)
ax.plot(edfs, cv_rss_average, 'o-', color='C1')
for k, (ave, stderr) in enumerate(zip(cv_rss_average, cv_rss_stderr)):
    ax.plot([k/2, k/2], [ave-stderr, ave+stderr], color='C0', linewidth=1)
    ax.plot([k/2-.1, k/2+.1], [ave-stderr, ave-stderr],
            color='C0', linewidth=1)
    ax.plot([k/2-.1, k/2+.1], [ave+stderr, ave+stderr],
            color='C0', linewidth=1)
ax.set_xlabel('Degrees of Freedom')
ax.set_ylabel('CV Error')
ax.set_title('Ridge Regression')
ax.set_xlim(-.2, 8.2)
ax.plot([-.2, 8.2], [cv_rss_average[-1]+cv_rss_stderr[-1], cv_rss_average[-1]+cv_rss_stderr[-1]], '--', color='purple', linewidth=0.5)
ax.plot([6, 6], [0.3, 1.4], '--', color='purple', linewidth=0.5)
ax.margins(0,0)
plt.show()

png

FIGURE 3.7. Estimated prediction error curves and their standard errors for the ridge shrinkage methods. The curve is plotted as a function of the corresponding complexity parameter for the method. The horizontal axis has been chosen so that the model complexity increases as we move from left to right. The estimates of prediction error and their standard errors were obtained by tenfold cross-validation; full details are given in Section 7.10. The least complex model within one standard error of the best (the horizontal dash line) is chosen, indicated by the purple vertical broken lines.

# The data_y is the 'lpsa' column of the original data 
data_y

1    -0.430783
2    -0.162519
3    -0.162519
4    -0.162519
5     0.371564
        ...   
93    4.385147
94    4.684443
95    5.143124
96    5.477509
97    5.582932
Name: lpsa, Length: 97, dtype: float64

# vec_y_test is the 'lpsa' column of the original data with train=='F'
data_y_test = data_y[mask_train == 'F']
vec_y_test = data_y_test.values
vec_y_test

array([0.7654678, 1.047319 , 1.047319 , 1.3987169, 1.6582281, 1.7316555,
       1.7664417, 1.8164521, 2.008214 , 2.0215476, 2.0856721, 2.3075726,
       2.3749058, 2.5687881, 2.5915164, 2.5915164, 2.6844403, 2.6912431,
       2.7047113, 2.7880929, 2.8535925, 2.8820035, 2.8820035, 2.8875901,
       3.0563569, 3.0750055, 3.5130369, 3.5709402, 5.1431245, 5.5829322])

"""Table 3.3. for the ridge regression. Use lambda = 5 as the textbook does.
"""
lambda_from_cv = lambdas[np.where(edfs == 5)[0]-1]  # lambda with edf = 5

u, s, vh = scipy.linalg.svd(data_x_train_normalized, full_matrices=False)
s2 = s**2

mat_diag = np.diag(s/(s2+lambda_from_cv))
beta_ridge = vh.T @ mat_diag @ u.T @ vec_y

# Use the zscore Centered and scaled test data x
mat_x_test = data_x_test_normalized.values
vec_y_test_fitted = mat_x_test @ beta_ridge

print('{0:>15} {1:>15}'.format('Term', 'Best Subset'))
print('-'*31)
print('{0:>15} {1:>15.3f}'.format('Intercept', beta_intercept))
for idx, col_name in enumerate(data.columns):
    print('{0:>15} {1:>15.3f}'.format(col_name, beta_ridge[idx]))
print('-'*31)
print('{0:>15} {1:>15.3f}'.format(
    'Test Error',
    #uncenter vec_y_test_fitted by adding beta_intercept
    ((vec_y_test-(vec_y_test_fitted+beta_intercept))**2).sum()/size_test),
)

           Term     Best Subset
-------------------------------
      Intercept           2.452
         lcavol           0.443
        lweight           0.255
            age          -0.051
           lbph           0.171
            svi           0.238
            lcp          -0.005
        gleason           0.041
          pgg45           0.137
-------------------------------
     Test Error           0.503

Singular Value Decomposition

from scipy.stats import multivariate_normal
# multivariate_normal(mean=array, cov=array)
rv = multivariate_normal([0, 0], [[10, 6], [6, 10]])
# .rvs() Draw random samples from a multivariate normal distribution.
X = rv.rvs(100)
sum(X)
# X is 100X2 array

array([-45.25513411, -54.32794596])

# np.linalg.svd() Singular Value Decomposition. mXn = (mXm)@(mXn)@(nXn)
# When a is a 2D array, it is factorized as u @ np.diag(s) @ vh = (u * s) @ vh, 
# where u and vh are 2D unitary arrays and s is a 1D array of a’s singular values. 
# When a is higher-dimensional, SVD is applied in stacked mode.
# full_matricesbool, optional If True (default), u and vh have the shapes (..., M, M) and (..., N, N), respectively. 
# Otherwise, the shapes are (..., M, K) and (..., K, N), respectively, where K = min(M, N).

U, d, Vt = np.linalg.svd(X, full_matrices=True)
# U is 100X100 array, d is 100X2 array, V is 2X2 array

V = Vt.T
d = np.sqrt(d)
print(U, '\n\n', d, '\n\n', V)

[[ 0.01  0.09 -0.07 ...  0.08  0.08  0.  ]
 [-0.03 -0.04  0.23 ...  0.02  0.11 -0.01]
 [ 0.06  0.24  0.94 ... -0.   -0.02  0.  ]
 ...
 [ 0.08 -0.02  0.   ...  0.99 -0.01 -0.  ]
 [ 0.12  0.06 -0.02 ... -0.01  0.98  0.  ]
 [-0.   -0.01  0.   ... -0.    0.    1.  ]] 

 [6.45 4.62] 

 [[ 0.71 -0.7 ]
 [ 0.7   0.71]]

V[0, 0]*d[0]

4.611849897053329

fig, ax1 = plt.subplots(figsize=(5, 5))
ax1.plot([0, V[0, 0]*d[0]], [0, V[1, 0]*d[0]])
ax1.plot([0, V[0, 1]*d[1]], [0, V[1, 1]*d[1]])
ax1.scatter(X[:,0], X[:,1],s=2)
ax1.set_xlim([-10, 10])
ax1.set_ylim([-10, 10])

(-10.0, 10.0)

png

FIGURE 3.9. Principal components of some input data points. The largest prin- cipal component is the direction that maximizes the variance of the projected data, and the smallest principal component minimizes that variance. Ridge regression projects y onto these components, and then shrinks the coefficients of the low– variance components more than the high-variance components.

Explicit solutions in the ideal situation

TABLE 3.4. In the case of an orthonormal input matrix \(\mathbf{X}\), the three procecures have explicit solutions as below. \(M\) and \(\lambda\) are constants chosen by the corresponding techniques; sign denotes the sign of its argument (±1), and \(x_+\) denotes “positive part” of \(x\). Below the table, estimators are shown by broken red lines. The \(45^\circ\) line in gray shows the unrestricted estimate for reference.

\[\begin{array}{l} \text{Estimator} & \text{Formula}\\ \hline \text{Best subset (size M)} &\hat{\beta}_j\cdot\mathbf I(\lvert\hat{\beta}_j\rvert\ge\lvert\hat{\beta}_{(M)}\rvert)\\ \text{Ridge} & \hat{\beta}_j\big/(1+\lambda)\\ \text{Lasso} & \text{sign}(\hat{\beta}_j)(\lvert\hat{\beta}_j\rvert-\lambda)_+\\ \hline \end{array}\]

• Best subset (size \(M\)) drops all variables with coefficients smaller than the \(M\)th largest; this is a form of “hard-thresholding.”
  \[\begin{equation} \hat\beta_j\cdot I\left( \lvert\hat\beta_j\rvert \ge \lvert\hat\beta_{(M)}\rvert \right) \end{equation}\]
• Ridge does a proportional shrinkage.
  \[\begin{equation} \frac{\hat\beta_j}{1+\lambda} \end{equation}\]
• Lasso translates each coefficient by a constant factor \(\lambda\), truncating at zero. This is called “soft-thresholding.”
  \[\begin{equation} \text{sign}\left( \hat\beta_j \right)\left( \lvert\hat\beta_j\rvert - \lambda \right)_+ \end{equation}\]

import scipy
import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(-1, 1, 100)

# Set \hat\beta_{M}=0.33
y_best_subset = x.copy()
y_best_subset[np.absolute(y_best_subset) <= .33] = 0

# Set \lambda=0.5
y_ridge = x/1.5

# Set \lambda=0.33
y_lasso = x.copy()
y_lasso[np.absolute(y_lasso) <= .33] = 0
y_lasso[y_lasso > .33] -= .33
y_lasso[y_lasso < -.33] += .33
        

fig34 = plt.figure(figsize=(9, 3))
ax1 = fig34.add_subplot(1, 3, 1)
ax1.plot(x, x, color='gray')
ax1.plot(x, y_best_subset, 'r--')
ax1.text(.4, .1, r'$|\hat\beta_{(M)}|$')
ax1.set_title('Best Subset', color = 'b')
ax1.grid()

ax2 = fig34.add_subplot(1, 3, 2)
ax2.plot(x, x, color='gray')
ax2.plot(x, y_ridge, 'r--')
ax2.set_title('Ridge', color = 'b')
ax2.grid()

ax3 = fig34.add_subplot(1, 3, 3)
ax3.plot(x, x, color='gray')
ax3.plot(x, y_lasso, 'r--')
ax3.plot([1, 1], [y_lasso[-1], 1.0], '--', color='b', linewidth=1)
ax3.text(1.0, .8, r'$\lambda$')
ax3.set_title('Lasso', color = 'b')
ax3.grid()
plt.show()

png

Back to the reality; the nonorthonormal case

Back to the nonorthogonal case; some pictures help understand their relationship. Figure 3.11 depicts the lasso (left) and ridge regression (right) when there are only two parameters. The residual sum of squares has elliptical contours, centered at the full least squares estimate. The constraint region for ridge regression is the disk \(\beta_1^2 + \beta_2^2 \leq t\), while that for lasso is the diamond \(|\beta_1| + |\beta_2| \leq t\). Both methods find the first point where the elliptical contours hit the constraint region. Unlike the disk, the diamond has corners; if the solution occurs at a corner, then it has one parameter \(\beta_j\) equal to zero. When \(p > 2\), the diamond becomes a rhomboid, and has many corners, flat edges and faces; there are many more opportunities for the estimated parameters to be zero.

import scipy
import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(-1, 1, 100)

# lasso
y_0 = x.copy()
y_1 = 1-np.absolute(y_0)
y_2 = -(1-np.absolute(y_0))


# ridge
y_3 = np.sqrt(1-np.absolute(y_0)**2)
y_4 = -(np.sqrt(1-np.absolute(y_0)**2))


fig35 = plt.figure(figsize=(8, 4))
ax1 = fig35.add_subplot(1, 2, 1)
ax1.fill(x, y_1, color='cyan')
ax1.fill(x, y_2, color='cyan')
ax1.axis('off')
ax1.set_xlim(-2, 5)
ax1.set_ylim(-2, 5)
ax1.arrow(-2, 0, dx=6, dy=0, color='0', head_width=0.1,width=.01, linewidth=0.5)
ax1.arrow(0, -2, dx=0, dy=6, color='0', head_width=0.1,width=.01, linewidth=0.5)
ax1.text(4.0, -.6, r'$\beta_1$')
ax1.text(-0.6, 4.0, r'$\beta_2$')

ax2 = fig35.add_subplot(1, 2, 2)
ax2.fill(x, y_3, color='cyan')
ax2.fill(x, y_4, color='cyan')
ax2.axis('off')
ax2.set_xlim(-2, 5)
ax2.set_ylim(-2, 5)
ax2.arrow(-2, 0, dx=6, dy=0, color='0', head_width=0.1,width=.01, linewidth=0.5)
ax2.arrow(0, -2, dx=0, dy=6, color='0', head_width=0.1,width=.01, linewidth=0.5)
ax2.text(4.0, -.6, r'$\beta_1$')
ax2.text(-0.6, 4.0, r'$\beta_2$')

plt.show()

png

FIGURE 3.11. Estimation picture for the lasso (left) and ridge regression (right). Shown are contours of the error and constraint functions. The solid blue areas are the constraint regions \(|\beta_1| + |\beta_2| \leq t\) and \(\beta_1^2 + \beta_2^2 \leq t^2\), respectively, while the red ellipses are the contours of the least squares error function.

Generalization to Bayes estimates

Consider the criterion, for \(q \ge 0\),

\[\begin{equation} \tilde\beta = \arg\min_\beta \left\lbrace \sum_{i=1}^N \left( y_i - \beta_0 - \sum_{j=1}^p x_{ij}\beta_j \right)^2 + \lambda\sum_{j=1}^p\lvert\beta_j\rvert^q\right\rbrace \end{equation}\]

See FIGURE 3.12 for the contours of \(\sum_j\lvert\beta_j\rvert^q\).

Thinking of \(\lvert\beta_j\rvert^q\) as the log-prior density for \(\beta_j\), these are also the equicontours of the prior distrobution of the parameters.

• \(q=0\) \(\Rightarrow\) variable subset selction, as the penalty simply counts the number of nonzero parameters,
• \(q=1\) \(\Rightarrow\) corresponds to the lasso, Laplace distribution with density

\[\begin{equation} \frac{1}{2\tau}\exp\left(-\lvert\beta\rvert/\tau\right), \end{equation}\] where \(\tau = 1/\lambda\). In this view, the lasso, ridge, and best subset selection are Bayes estimates with different priors. Note, however, that they are derived as posterior modes, that is, maximizers of the posterior. It is more common to use the mean of the posterior as the Bayes estimate. Ridge solution is also the posterior mean, but the lasso and best subset selection are not.

Note that the lasso case is the smallest \(q\) such that the constraint region is convex; non-convext constraint regions make the optimization problem more difficult.

In this view, the lasso, ridge regression and best subset selection are Bayes estimates with different priors. Note, however, that they are derived as posterior modes, that is, maximizers of the posterior. It is more common to use the mean of the posterior as the Bayes estimate. Ridge regression is also the posterior mean, but the lasso and best subset selection are not.

We might try using other values of \(q\) besides \(0\), \(1\), or \(2\). Although one might consider estimating \(q\) from data, our experience is that it is not worth the effort for the extra variance incurred.

Although this is the case, with \(q \gt 1\), \(\lvert\beta_j\rvert^q\) is differentiable at \(0\), and so does not share the ability of lasso (\(q=1\)) for setting coefficients exactly to zero. Partly for this reason as well as for computational tractability, Zou and Hastie (2005) introduced the elastic-net penalty, introduced in \(\S\) 18.4.

x = np.linspace(-1, 1, 100)

# q=4
y_0 = x.copy()
y_1 = np.sqrt(np.sqrt((1-np.absolute(y_0)**4)))
y_2 = -np.sqrt(np.sqrt((1-np.absolute(y_0)**4)))

# ridge q=2
y_3 = np.sqrt(1-np.absolute(y_0)**2)
y_4 = -(np.sqrt(1-np.absolute(y_0)**2))

# lasso q=1
y_0 = x.copy()
y_5 = 1-np.absolute(y_0)
y_6 = -(1-np.absolute(y_0))


# q=0.5
y_7 = (1-np.sqrt(np.absolute(y_0)))**2
y_8 = -(1-np.sqrt(np.absolute(y_0)))**2

# q=0.1
y_9 = (1-(np.absolute(y_0))**0.1)**10
y_10 = -(1-(np.absolute(y_0))**0.1)**10


fig36 = plt.figure(figsize=(10, 2))
ax1 = fig36.add_subplot(1, 5, 1)
ax1.plot(x, y_1, color='cyan')
ax1.plot(x, y_2, color='cyan')
ax1.axis('off')
ax1.plot([-1.2,1.2],[0,0], color='0', linewidth=1)
ax1.plot([0,0],[-1.2,1.2], color='0', linewidth=1)
ax1.set_xlim(-2, 2)
ax1.set_ylim(-2, 2)
ax1.set_title('q=4')

ax2 = fig36.add_subplot(1, 5, 2)
ax2.plot(x, y_3, color='cyan')
ax2.plot(x, y_4, color='cyan')
ax2.axis('off')
ax2.plot([-1.2,1.2],[0,0], color='0', linewidth=1)
ax2.plot([0,0],[-1.2,1.2], color='0', linewidth=1)
ax2.set_xlim(-2, 2)
ax2.set_ylim(-2, 2)
ax2.set_title('q=2')

ax3 = fig36.add_subplot(1, 5, 3)
ax3.plot(x, y_5, color='cyan')
ax3.plot(x, y_6, color='cyan')
ax3.axis('off')
ax3.plot([-1.2,1.2],[0,0], color='0', linewidth=1)
ax3.plot([0,0],[-1.2,1.2], color='0', linewidth=1)
ax3.set_xlim(-2, 2)
ax3.set_ylim(-2, 2)
ax3.set_title('q=1')

ax4 = fig36.add_subplot(1, 5, 4)
ax4.plot(x, y_7, color='cyan')
ax4.plot(x, y_8, color='cyan')
ax4.axis('off')
ax4.plot([-1.2,1.2],[0,0], color='0', linewidth=1)
ax4.plot([0,0],[-1.2,1.2], color='0', linewidth=1)
ax4.set_xlim(-2, 2)
ax4.set_ylim(-2, 2)
ax4.set_title('q=0.5')

ax5 = fig36.add_subplot(1, 5, 5)
ax5.plot(x, y_9, color='cyan')
ax5.plot(x, y_10, color='cyan')
ax5.axis('off')
ax5.plot([-1.2,1.2],[0,0], color='0', linewidth=1)
ax5.plot([0,0],[-1.2,1.2], color='0', linewidth=1)
ax5.set_xlim(-2, 2)
ax5.set_ylim(-2, 2)
ax5.set_title('q=0.1')

plt.show()

png

FIGURE 3.12. Contours of constant value of \(\sum_{j}|\beta_j|^q\) for given values of \(q\).

Algorithm 3.2. Least Angle Regression

1. Standardize the predictors to have mean zero and unit norm.
   Start with the residual \(r = y - \bar{y}\),and coefficients \(\beta_0 = \beta_1 = \cdots = \beta_p = 0\). Since with all \(\beta_j = 0\) and standardized predictors the constant coefficient \(\beta_0=\bar{y}\).

2. Set \(k = 1\) and begin start the \(k\)-th step. Since all values of \(\beta_j\) are zero the first residual is \(r_1 = y - \bar{y}\). Find the predictor \(x_j\) that is most correlated with this residual \(r_1\). Then as we begin this \(k = 1\) step we have the active step given by \(A_1 = \{x_j\}\) and the active coefficients given by \(\beta_{A_1} = [0]\).

3. Move \(\beta_j\) from \(0\) towards its least-squares coefficient \[\delta_1=\left(X_{A_1}^TX_{A_1}\right)^{-1}X_{A_1}^Tr_1=\frac{x_j^Tr_1}{x_j^Tx_j}=x_j^T r_1\] until some other competitor \(x_k\) has as much correlation with the current residual as does \(x_j\). The path taken by the elements in \(\beta_{A_1}\) can be parameterized by \[\beta_{A_1}(\alpha)\equiv\beta_{A_1}+\alpha\delta_1=0+\alpha x_j^T r_1=(x_j^T r_1)\alpha,\quad(0\leq\alpha\leq1)\] This path of the coefficients \(\beta_{A_1}(\alpha)\) will produce a path of fitted values given by \[\hat{f}_1(\alpha)=X_{A_1}\beta_{A_1}(\alpha)=(x_j^T r_1)\alpha x_j\] and a residual of \[r(\alpha)=y-\bar{y}-\hat{f}_1(\alpha)=y-\bar{y}-(x_j^T r_1)\alpha x_j=r_1-(x_j^T r_1)\alpha x_j\] Now at this point \(x_j\) itself has a correlation with this residual as \(\alpha\) varies given by \[x_j^T(r_1-(x_j^T r_1)\alpha x_j)=x_j^Tr_1-(x_j^T r_1)\alpha=(1-\alpha)x_j^T r_1\]

4. Move \(\beta_j\) and \(\beta_k\) in the direction defined by their joint least squares coefficient of the current residual on \((x_j,x_k)\), until some other competitor \(x_l\) has as much correlation with the current residual.

5. Continue in this way until all \(p\) predictors have been entered. After \(\min(N-1,p)\) steps, we arrive at the full least-squares solution.

The termination condition in step 5 requires some explanation. If \(p > N-1\), the LAR algorithm reaches a zero residual solution after \(N-1\) steps (the \(-1\) is because we have centered the data).

Suppose

• \(\mathcal{A}_k\) is the active set of variables at the beginning of the \(k\)th step,
• \(\beta_{\mathcal{A}_k}\) is the coefficient vector for these variables at this step. There will be \(k-1\) nonzero values, and the one just entered will be zero. If \(\mathbf{r}_k = \mathbf{y} - \mathbf{X}_{\mathcal{A}_k}\beta_{\mathcal{A}_k}\) is the current residual, then the direction for this step is

\[\begin{equation} \delta_k = \left( \mathbf{X}_{\mathcal{A}_k}^T\mathbf{X}_{\mathcal{A}_k} \right)^{-1}\mathbf{X}_{\mathcal{A}_k}^T \mathbf{r}_k\\ =\left(\mathbf{V}\mathbf D^2\mathbf{V}^T\right)^{-1}\mathbf{V}\mathbf D\mathbf U^T\mathbf{r}_k\\ =\mathbf{V}^{-T}\mathbf D^{-2}\mathbf D\mathbf U^T\mathbf{r}_k\\ \end{equation}\]

The coefficient profile then evolves as

\[\begin{equation} \beta_{\mathcal{A}_k}(\alpha) = \beta_{\mathcal{A}_k} + \alpha\delta_k, \end{equation}\]

and this direction keeps the correlations tied and decreasing (See Exercise 3.23).

If the fit vector at the beginning of this step is \(\hat{\mathbf{f}}_k\), then it evolves as

\[\begin{equation} \hat{\mathbf{f}}_k(\alpha) = \hat{\mathbf{f}}_k + \alpha\mathbf{u}_k, \end{equation}\]

where \(\mathbf{u}_k = \mathbf{X}_{\mathcal{A}_k}\delta_k\) is the new fit direction.

The name “least angle” arises from a geometrical interpretation of this process; \(\mathbf{u}_k\) makes the smallest (and equal) angle with each of the predictors in \(\mathcal{A}_k\) (Exercise 3.24).

Note that, by construction, the coefficients in LAR change in a piecewise linear fashion, and we do not need to take small steps and recheck the correlations in step 3; using knowledge of the covariance of the predictors and the piecewise linearity of the algorithm, we can work out the exact step length at the beginning of each step (Exercise 3.25).

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

fig = plt.figure(figsize=(5, 5))
ax = fig.add_subplot(1, 1, 1)

x = np.arange(0, 2, 0.01)
ax.plot([0,1], [0,1], color='0', linewidth=1)
ax.plot([1, 0.5],[0, 0.5], color='0', linewidth=1)
ax.set_xlim(0,1)
ax.text(0.5, -0.04, 'r', size=15)
ax.arrow(0.01, 0.05, dx=0.4, dy=0.4, color='0', head_width=0.05,width=.01, linewidth=0.5)
ax.text(0.3, 0.45, 'u', size=15)
ax.set_title('Animation of the LAR moves the coefficient of \n a variable continuously toward its least-sqaures value')

for n in range(50):
    plt.plot([1, x[n]],[0, x[n]], color='cyan')
    
plt.show()

png

"""FIGURE 3.10. LARS(lasso) coefficients profile on prostate cancer data"""
import scipy
import scipy.stats
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

df = pd.read_csv('../../data/prostate/prostate.data', delimiter='\t',
                 index_col=0)
dfy = df.pop('lpsa')
dfmask = df.pop('train')
size_predictor = df.columns.size

dfmask_train = dfmask == 'T'
dfx_train = df[dfmask_train]
dfy_train = dfy[dfmask_train]
size_train = dfmask_train.sum()

dfmask_test = dfmask == 'F'
dfx_test = df[dfmask_train]
dfy_test = dfy[dfmask_train]
size_test = dfmask_train.sum()

# Centering for the training data
meany_train = dfy_train.mean()
dfy_train_centered = dfy_train.subtract(meany_train)
dfx_train_centered = dfx_train.subtract(dfx_train.mean())

normx = np.sqrt((dfx_train_centered**2).sum())
dfx_train_scaled = dfx_train_centered.divide(normx)

matx = dfx_train_scaled.values
vecy = dfy_train_centered.values

"""Main loop for LAR

Reference:
Least Angle Regression, Efron et al., 2004, The Annals of Statistics
"""
# Initial data
signs = np.zeros(size_predictor)
betas = np.zeros(size_predictor)
indices_predictor = np.arange(size_predictor)
vecy_fitted = np.zeros_like(vecy)
beta_lars = [[0]*size_predictor]

for k in range(size_predictor):
    #vecc=X^T\bar{y}
    vecc = matx.T @ (vecy-vecy_fitted)
    vecc_abs = np.absolute(vecc)
    
    #maxc is x_j\bar{y} of the most correlated x_j
    maxc = vecc_abs.max()
    mask_maxc = np.isclose(vecc_abs, maxc)
    active = indices_predictor[mask_maxc]
    signs = np.where(vecc[active] > 0, 1, -1)
    
    #get the active column (predictor) of matx
    matx_active = signs*matx[:, active]

    u, s, vh = scipy.linalg.svd(matx_active, full_matrices=False)
    
    # X^TX = VD^2V^T
    matg = vh.T @ np.diag(s**2) @ vh
    # (X^TX)^-1 = VD^-2V^T
    matg_inv = vh.T @ np.diag(np.reciprocal(s**2)) @ vh
    vec1 = np.ones(len(active))
    scalara = (matg_inv.sum())**(-.5)
    
    #sum in the rows of matg_inv, equal to (X^TX)^-1 @ 1
    vecw = scalara * matg_inv.sum(axis=1)
    vecu = matx_active @ vecw

    veca = matx.T @ vecu

    if k < size_predictor-1:
        inactive = indices_predictor[np.invert(mask_maxc)]
        arr_gamma = np.concatenate([(maxc-vecc[inactive])/(scalara-veca[inactive]),
                                       (maxc+vecc[inactive])/(scalara+veca[inactive])])
        scalargamma = arr_gamma[arr_gamma > 0].min()
    else:
        scalargamma = maxc/scalara

    vecy_fitted += scalargamma*vecu
    betas[active] += scalargamma*signs
    beta_lars.append(list(betas))

fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(1, 1, 1)
l1length_ols_beta = sum(abs(coeff) for coeff in beta_lars[-1])
l1length_lars_beta = [sum(abs(coeff) for coeff in beta)/l1length_ols_beta for beta in beta_lars]
ax.plot(l1length_lars_beta, beta_lars, 'o-')
ax.grid()
ax.set_xlabel('Shrinkage Factor s')
ax.set_ylabel('Coefficients')
ax.legend(dfx_train.columns)
ax.set_xlim(-.2, 1.2)
ax.set_ylim(-2.2, 11.0)
for x, y, s in zip(np.ones(8)*(1+0.02), beta_lars[-1], df.columns):
    ax.text(x, y, s, color='0', fontsize=12)
ax.margins(0,0)
plt.show()

png

For the prostate data, the LAR coefficient profile turns out to be identical to the lasso profile in FIGURE 3.10, which never crosses zero. These observations lead to a simple modification of the LAR algorithm that gives the entire lasso path, which is also piecewise-linear.

Algorithm 3.2a. Least Angle Regression: Lasso Modification

4a. If a non-zero coefficient hits zero, drop its variable from the active set and recompute the current joint least squares direction.

Computational efficiency

The LAR (lasso) algorithm is extremely efficient, requiring the same order of computation as that of a single least squares fit using the \(p\) predictors. LAR always take \(p\) steps to get to the full least squares estimates. The lasso path can have more than \(p\) steps, although the two are often quite similar. Also the LAR algorithm with the lasso modification is efficient, especially when \(p \gg N\). Osborne et al. (2000) also discovered a piecewise-linear path for computing the lasso, which they called a homotopy algorithm.

Heuristic argument for why LAR and lasso are so similar

Although the LAR algorithm is stated in terms of correlations, if the input features are standardized, it is equivalent and easier to work with inner-products.

Let \(\mathcal{A}\) be the active set at some stage, tied in their absolute inner-product with the current residuals \(\mathbf{y}-\mathbf{X}\beta\), expressed as

\[\begin{equation} \mathbf{x}_j^T\left( \mathbf{y}-\mathbf{X}\beta \right) = \gamma\cdot s_j, \forall j\in\mathcal{A}, \end{equation}\]

where \(s_j\in\lbrace -1, 1\rbrace\) indicates the sign of the inner-product, and \(\gamma\) is the common value. Also note that

\[\begin{equation} \lvert \mathbf{x}_k^T \left( \mathbf{y}-\mathbf{X}\beta \right) \rvert \le \gamma, \forall k \notin \mathcal{A} \end{equation}\]

Now consider the lasso criterion

\[\begin{equation} \hat\beta^{\text{lasso}} = \arg\min_\beta \left\lbrace \frac{1}{2}\sum_{i=1}^N \left( y_i-\beta_0-\sum_{j=1}^p x_{ij}\beta_j \right)^2 + \lambda\sum_{j=1}^p |\beta_j| \right\rbrace, \end{equation}\]

which we write in the vector form

\[\begin{equation} R(\beta) = \frac{1}{2}\|\mathbf{y}-\mathbf{X}\beta\|^2_2 + \lambda\|\beta\|_1. \end{equation}\]

Let \(\mathcal{B}\) be the active set in the solution for a given value of \(\lambda\). For these variables \(R(\beta)\) is differentiable, and the stationary conditions give

\[\begin{align*} \frac{\partial R}{\partial\beta} &= -\mathbf{X}^T(\mathbf{y}-\mathbf{X}\beta) + \lambda \frac{\partial \|\beta\|_1}{\partial\beta} = 0 \\ \mathbf{x}^T_j\left( \mathbf{y}-\mathbf{X}\beta \right) &= \lambda\cdot\text{sign}(\beta_j), \forall j\in\mathcal{B} \end{align*}\]

This is identical to the above inner-product expression only if the sign of \(\beta_j\) matches the sign of the inner-product. This is why the LAR and lasso start to differ when an active coefficient passes through zero; the above condition is violated for that variable, and it is kicked out of the active set \(\mathcal{B}\). Exercise 3.23 shows that these equations imply a piecewise linear coefficient profile as \(\lambda\) decreases.

The stationary conditions for the non-active variables require that

\[\begin{equation} \lvert \mathbf{x}_k^T(\mathbf{y}-\mathbf{X}\beta) \rvert \le \lambda, \forall k\notin\mathcal{B}, \end{equation}\]

which again agrees with the LAR algorithm.

Degrees-of-freedom formula for LAR and lasso

Suppose that we fit a linear model via the LAR procedure, stopping at some number of steps \(k<p\), or equivalently using a lasso bound \(t\) that produces a constrained version of the full least squares fit. How many parameters, or “degrees of freedom” have we used?

In classical statistics, the number of linearly independent parameters is what is meant by “degrees of freedom.” So a least squares model with \(k\) features has degress of freedom to be \(k\). Alternatively, suppose that we carry out a best subset selection to determine the “optimal” set of \(k\) predictors. Then the resulting model has \(k\) parameters, but in some sense we have used up more than \(k\) degrees of freedom.

We need a more general defintion for the effective degrees of freedom of an adaptively fitted model. Define the degrees of freedom of the fitted vector \(\hat{\mathbf{y}}\) as

\[\begin{equation} \text{df}(\hat{\mathbf{y}}) = \frac{1}{\sigma^2}\sum_{i=1}^N \text{Cov}(\hat{y}_i, y_i), \end{equation}\]

where \(\text{Cov}(\hat{y}_i, y_i)\) refers to the sampling covariance between \(\hat{y}_i\) and \(y_i\). This makes intuitive sense: The harder that we fit to the data, the larger this covariance and hence \(\text{df}(\hat{\mathbf{y}})\). This is a useful notion of degrees of freedom, one that can be applied to any model prediction \(\hat{\mathbf{y}}\), including models adaptively fitted to the training data. This definition is motivated and discussed further in \(\S\) 7.4-7.6.

We relate this scalar expression into a vector inner product expression as \[\begin{align} \text{df}(\hat{\mathbf{y}}) &= \frac{1}{\sigma^2}\sum_{i=1}^N \text{Cov}(\hat{y}_i, y_i)\\ &=\frac{1}{\sigma^2}\sum_{i=1}^N \text{Cov}(\mathbf{e}_i^T\hat{\mathbf{y}}, \mathbf{e}_i^T\mathbf{y})\\ &=\frac{1}{\sigma^2}\sum_{i=1}^N \mathbf{e}_i^T\text{Cov}(\hat{\mathbf{y}}, \mathbf{y})\mathbf{e}_i\\ \end{align}\]

Now for ordinary least squares regression we have \[\hat{\mathbf{y}}=\mathbf{X}\hat{\beta}^{ls}=\mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}\] so that the above expression for \(\text{Cov}(\hat{\mathbf{y}}, \mathbf{y})\) becomes \[\text{Cov}(\hat{\mathbf{y}}, \mathbf{y})=\text{Cov}(\mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}, \mathbf{y})=\mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\text{Cov}(\mathbf{y}, \mathbf{y})=\mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\sigma^2\] Thus \[\text{Cov}(\hat{y}_i, y_i)=\mathbf{e}_i^T\text{Cov}(\hat{\mathbf{y}}, \mathbf{y})\mathbf{e}_i^T=\mathbf{e}_i^T\mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\sigma^2\mathbf{e}_i\]

Note that \(\mathbf{X}^T\mathbf{e}_i=\mathbf{x}_i\) is the \(i\)th samples feature vector for \(1\leq i\leq N\) and we have \[\text{Cov}(\hat{\mathbf{y}}, \mathbf{y})=\mathbf{x}_i^T(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{x}_i\sigma^2\] which when we sum for \(i = 1\) to \(N\) and divide by \(\sigma^2\) gives

\[\begin{align} \text{df}(\hat{\mathbf{y}}) &= \frac{1}{\sigma^2}\sum_{i=1}^N \text{Cov}(\hat{y}_i, y_i)\\ &=\frac{1}{\sigma^2}\sum_{i=1}^N \text{Cov}(\mathbf{e}_i^T\hat{\mathbf{y}}, \mathbf{e}_i^T\mathbf{y})\\ &=\frac{1}{\sigma^2}\sum_{i=1}^N \mathbf{e}_i^T\text{Cov}(\hat{\mathbf{y}}, \mathbf{y})\mathbf{e}_i\\ &=\sum_{i=1}^N \mathbf{x}_i^T(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{x}_i\\ &=\sum_{i=1}^N \text{trace}\left(\mathbf{x}_i^T(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{x}_i\right)\\ &=\sum_{i=1}^N \text{trace}\left(\mathbf{x}_i\mathbf{x}_i^T(\mathbf{X}^T\mathbf{X})^{-1}\right)\\ &=\text{trace}\left(\sum_{i=1}^N \left(\mathbf{x}_i\mathbf{x}_i^T\right)(\mathbf{X}^T\mathbf{X})^{-1}\right)\\ \end{align}\] Since \[ \sum_{i=1}^N \left(\mathbf{x}_i\mathbf{x}_i^T\right)=\begin{bmatrix} \mathbf{x}_1 & \mathbf{x}_2 & \cdots & \mathbf{x}_N \end{bmatrix}\begin{bmatrix} \mathbf{x}_1^T \\ \mathbf{x}_2^T \\ \vdots \\ \mathbf{x}_N^T\\ \end{bmatrix}=\mathbf{X}^T\mathbf{X} \]

Thus when there are \(k\) predictors we get \[\text{df}(\hat{\mathbf{y}})=\text{trace}\left(\mathbf{X}^T\mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1}\right)=\text{trace}\left(\mathbf{I}\right)=k\]

To do the same thing for ridge regression: \[\begin{equation} \hat{\mathbf{y}}=\mathbf{X}\hat\beta^{\text{ridge}} = \mathbf{X}\left( \mathbf{X}^T\mathbf{X} + \lambda\mathbf{I} \right)^{-1}\mathbf{X}^T\mathbf{y}, \end{equation}\]

so that \[\text{Cov}(\hat{y}_i, y_i)=\mathbf{e}_i^T\text{Cov}(\hat{\mathbf{y}}, \mathbf{y})\mathbf{e}_i^T=\mathbf{e}_i^T\mathbf{X}\left( \mathbf{X}^T\mathbf{X} + \lambda\mathbf{I} \right)^{-1}\mathbf{X}^T\sigma^2\mathbf{e}_i\\ =\sigma^2\left(\mathbf{X}^T\mathbf{e}_i\right)^T\left( \mathbf{X}^T\mathbf{X} + \lambda\mathbf{I} \right)^{-1}\mathbf{X}^T\mathbf{e}_i\\ =\sigma^2\mathbf{x}_i^T\left( \mathbf{X}^T\mathbf{X} + \lambda\mathbf{I} \right)^{-1}\mathbf{x}_i\] Then summing for \(i = 1, 2, \cdots ,N\) and dividing by \(\sigma^2\) to get \[\begin{align} \text{df}(\hat{\mathbf{y}})&=\text{trace}\left(\mathbf{X}^T\mathbf{X}\left(\mathbf{X}^T\mathbf{X} + \lambda\mathbf{I} \right)^{-1}\right)\\ &=\text{trace}\left(\mathbf{X}\left(\mathbf{X}^T\mathbf{X} + \lambda\mathbf{I} \right)^{-1}\mathbf{X}^T\right) \end{align}\]

Now

• for a linear regression with \(k\) fixed predictors, it is easy to show that \(\text{df}(\hat{\mathbf{y}}) = k\).
• Likewise for ridge regression, \(\text{df}(\hat{\mathbf{y}}) = \text{trace}(\mathbf{S}_\lambda)\).

In both these cases, \(\text{df}\) is simple to evaluate because the fit

\[\begin{equation} \hat{\mathbf{y}} = \mathbf{H}_\lambda \mathbf{y} \end{equation}\]

is linear in \(\mathbf{y}\). If we think about the above definition of degrees of freedom in the context of a best subset selection of size \(k\), it seems clear that \(\text{df}(\hat{\mathbf{y}}) > k\), and this can be verified by estimating \(\text{Cov}(\hat{y}_i,y_i)/\sigma^2\) directly by simulation. However there is no closed form method for estimating \(\text{df}\) for best subset selection.

For LAR and lasso, estimation of \(\text{df}\) is more tractable since these techniques are adaptive in a smoother way than best subset selection.

• After the \(k\)th step of the LAR procedure, \(\text{df} = k\).
• For the lasso, the modified LAR procedure often takes more than \(p\) steps, since predictors can drop out. Hence the definition is a little different; for the lasso, at any stage \(\text{df}\) approximately equals the number of predictors in the model.
  While this approximation works reasonably well anywhere in the lasso path, for each \(k\) it works best at the last model in the sequence that contains \(k\) predictors.

A detailed study of the effective degrees of freedom for the lasso may be found in Zou et al. (2007).

Comparison with ridge regression

If \(M=p\), since the columns of \(\mathbf{Z} = \mathbf{UD}\) span the \(\text{col}(\mathbf{X})\),

\[\begin{equation} \hat\beta^{\text{pcr}}(p) = \hat\beta^{\text{ls}}. \end{equation}\]

For \(M<p\) we get a reduced regression and we see that PCR is very similar to ridge regression: both operate via the principal components of the input matrix.

• Ridge regression shrinks the coefficients of the principal components (FIGURE 3.17), shrinking more depending on the size of the corresponding eigenvalue;
• PCR discards the \(p-M\) smallest eigenvalue components. FIGURE 3.17 illustrates this.

import scipy
import scipy.stats
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

data = pd.read_csv('../../data/prostate/prostate.data', delimiter='\t',
                   index_col=0)
data_y = data.pop('lpsa')
mask_train = data.pop('train')

data_x_train = data[mask_train == 'T']
data_y_train = data_y[mask_train == 'T']
beta_intercept = data_y_train.mean()
# Centering for the training data y by subtracting the mean
data_y_train_centered = data_y_train.subtract(beta_intercept)
# Centering and scaling for the training data x using zscore
data_x_train_normalized = data_x_train.apply(scipy.stats.zscore)
vec_y = data_y_train_centered.values

data_x_test = data[mask_train == 'F']
# Centering and scaling for the test data x using zscore
data_x_test_normalized = data_x_test.apply(scipy.stats.zscore)
data_y_test = data_y[mask_train == 'F']
vec_y_test = data_y_test.values

size_train = sum(mask_train == 'T')
size_test = sum(mask_train == 'F')
size_predictor = data_x_train.columns.size

def lambdas_from_edf(singular_squared:np.ndarray, interval:int) ->np.ndarray:
    """Given squared singular values of data matrix, calculate the lambdas
    with `interval` parameter to split unit intervals s.t. the resulting
    effective degrees of freedom are equidistant with 1/interval, via the
    Newton-Raphson method. e.g., if interval = 10, it produces lambdas for
    0, 0.5, 0.6, 0.7, ..."""
    p = singular_squared.size
    edfs = np.linspace(.5, p-.5, (p-1)*interval+1)
    threshold = 1e-3
    lambdas = []
    for edf in edfs:
        # Newton-Raphson
        lambda0 = (p-edf)/edf
        lambda1 = 1e6
        diff = lambda1 - lambda0
        while diff > threshold:
            num = (singular_squared/(singular_squared+lambda0)).sum()-edf
            denom = (singular_squared/((singular_squared+lambda0)**2)).sum()
            lambda1 = lambda0 + num/denom
            diff = lambda1 - lambda0
            lambda0 = lambda1
        lambdas.append(lambda1)
    lambdas.append(0)
    edfs = np.concatenate(([0], edfs, [p]))
    return edfs, np.array(lambdas)

# singular value decomposition of data_x_train_normalized
u, s, vh = scipy.linalg.svd(data_x_train_normalized, full_matrices=False)
s2 = s**2
edfs, lambdas = lambdas_from_edf(s2, 10)
# beta_ols = vh.T @ scipy.diag(scipy.reciprocal(s)) @ u.T @ vec_y
# print(beta_ols)
beta_ridge_array = [np.zeros(size_predictor)]
shrinkage_factors = []
for lamb in lambdas:
    shrinkage = s2/(s2+lamb)
    shrinkage_factors.append(shrinkage)
    shrinkage_factors_array = np.array(shrinkage_factors)

ridge = shrinkage_factors_array.mean(axis=0)
pcr = np.array((1,1,1,1,1,1,1,0,0))


fig = plt.figure(figsize=(10, 6))
ax = fig.add_subplot(1, 1, 1)
ax.plot(range(1, 9), ridge, 'o-', markersize=2)
ax.plot(np.array((1,2,3,4,5,6,7,7,8)), pcr, 'o-', markersize=2)
ax.legend(('ridge','pcr'))
ax.set_xlabel('Index')
ax.set_ylabel('Shrinkage Factor')
ax.set_ylim(-0.1,1.1)
plt.show()

png

FIGURE 3.17. Ridge regression shrinks the regression coefficients of the prin- cipal components, using shrinkage factors \(d_j^2 /(d_j^2 + \lambda)\) as in (3.47). Principal component regression truncates them. Shown are the shrinkage and truncation patterns corresponding to Figure 3.7, as a function of the principal component index.

Algorithm 3.3. Partial least squares.

1. Standardized each \(\mathbf{x}_j\) to have mean \(0\) and variance \(1\).
   Set \[\begin{align} \hat{\mathbf{y}}^{(0)} &= \bar{y}\mathbf{1} \\ \mathbf{x}_j^{(0)} &= \mathbf{x}_j, \text{ for } j=1,\cdots,p. \end{align}\]

2. For \(m = 1,2,\cdots,p\)

• \(\mathbf{z}_m = \sum_{j=1}^p \hat\varphi_{mj}\mathbf{x}_j^{(m-1)}\), where \(\hat\varphi_{mj} = \langle \mathbf{x}_j^{(m-1)},\mathbf{y}\rangle\).
• \(\hat\theta_m = \langle\mathbf{z}_m,\mathbf{y}\rangle \big/ \langle\mathbf{z}_m,\mathbf{z}_m\rangle\).
• \(\hat{\mathbf{y}}^{(m)} = \hat{\mathbf{y}}^{(m-1)} + \hat\theta_m \mathbf{z}_m\).
• Orthogonalize each \(\mathbf{x}_j^{(m-1)}\) w.r.t. \(\mathbf{z}_m\):
  \(\mathbf{x}_j^{(m)} = \mathbf{x}_j^{(m-1)} - \frac{\langle\mathbf{z}_m,\mathbf{x}_j^{(m-1)}\rangle}{\langle\mathbf{z}_m,\mathbf{y}\rangle}\mathbf{z}_m, \text{ for } j=1,2,\cdots,p\).

3. Output the sequence of fitted vectors \(\left\lbrace \hat{\mathbf{y}}^{(m)}\right\rbrace_1^p\).
   Since the \(\left\lbrace \mathbf{z}_l \right\rbrace_1^m\) are linear in the original \(\mathbf{x}_j\), so is
   \[\begin{equation} \hat{\mathbf{y}}^{(m)} = \mathbf{X}\hat\beta^{\text{pls}}(m). \end{equation}\] These linear coefficients can be recovered from the sequence of PLS transformations.

Gist of the PLS algorithm

PLS begins by computing the weights \[ \begin{equation} \hat\varphi_{1j} = \langle \mathbf{x}_j,\mathbf{y} \rangle, \text{ for each } j, \end{equation} \] which are in fact the univariate regression coefficients, since \(\mathbf{x}_j\) are standardized (only for the first step \(m=1\)).

From this we construct derived input \[ \begin{equation} \mathbf{z}_1 = \sum_j \hat\varphi_{1j}\mathbf{x}_j, \end{equation} \] which is the first PLS direction. Hence in the construction of each \(\mathbf{z}_m\), the inputs are weighted by the strength of their univariate effect on \(\mathbf{y}\).

The outcome \(\mathbf{y}\) is regressed on \(\mathbf{z}_1\) giving coefficient \(\hat\theta_1\), and then we orthogonalize \(\mathbf{x}_1,\cdots,\mathbf{x}_p\) w.r.t. \(\mathbf{z}_1\).

We continue this process, until \(M\le p\) directions have been obtained. In this manner, PLS produces a sequence of derived, orthogonal inputs or directions \(\mathbf{z}_1,\cdots,\mathbf{z}_M\).

• As with PCR, if \(M=p\), then \(\hat\beta^{\text{pls}} = \hat\beta^{\text{ls}}\).
• Using \(M<p\) directions produces a reduced regression.

Relation to the optimization problem

PLS seeks direction that have high variance and have high correlation with the response, in contrast to PCR with keys only on high variance (Stone and Brooks, 1990; Frank and Friedman, 1993).

Since it uses the response \(\mathbf{y}\) to construct its directions, its solution path is a nonlinear function of \(\mathbf{y}\).

In particular, the \(m\)th principal component direction \(v_m\) solves: \[ \begin{equation} \max_\alpha \text{Var}(\mathbf{X}\alpha)\\ \text{subject to } \|\alpha\| = 1, \alpha^T\mathbf{S} v_l = 0, \text{ for } l = 1,\cdots, m-1, \end{equation} \] where \(\mathbf{S}\) is the sample covariance matrix of the \(\mathbf{x}_j\). The condition \(\alpha^T\mathbf{S} v_l= 0\) ensures that \(\mathbf{z}_m = \mathbf{X}\alpha\) is uncorrelated with all the previous linear combinations \(\mathbf{z}_l = \mathbf{X} v_l\).

The \(m\)th PLS direction \(\hat\varphi_m\) solves: \[ \begin{equation} \max_\alpha \text{Corr}^2(\mathbf{y},\mathbf{S}\alpha)\text{Var}(\mathbf{X}\alpha)\\ \text{subject to } \|\alpha\| = 1, \alpha^T\mathbf{S}\hat\varphi_l = 0 \text{ for } l=1,\cdots, m-1. \end{equation} \] Further analysis reveals that the variance aspect tends to dominate, and so PLS behaves much like ridge regression and PCR. We discuss further in the next section.

If the input matrix \(\mathbf{X}\) is orthogonal, then PLS finds the least squares estimates after the first \(m=1\) step, and subsequent steps have no effect since the \(\hat\varphi_{mj} = 0\) for \(m>1\) (Exercise 3.14).

It can be also shown that the sequence of PLS coefficients for \(m=1,2,\cdots,p\) represents the conjugate gradient sequence for computing the least squares solutions (Exercise 3.18).

Algorithm 3.4 Incremental Forward Stagewise Regression—\(FS_{\varepsilon}\).

1. Start with the residual \(r\) equal to \(y\) and \(\beta_1, \beta_2, \cdots, \beta_p = 0\). All the predictors are standardized to have mean zero and unit norm.

2. Find the predictor \(x_j\) most correlated with \(r\)

3. Update \(\beta_j \gets \beta_j + \delta_j\) , where \[\delta_j=\varepsilon\cdot\text{sign}[\langle x_j , r\rangle]\] and \(\varepsilon > 0\) is a small step size, and set \[r \gets r − \delta_jx_j\].

4. Repeat steps 2 and 3 many times, until the residuals are uncorrelated with all the predictors.

If \[\delta_j=\langle x_j , r\rangle\] (the least-squares coefficient of the residual on \(j\)th predictor), then this is exactly the usual forward stagewise procedure (FS) outlined in Section 3.3.3.

Efron originally thought that the LAR Algorithm 3.2 was an implementation of \(FS_0\), allowing each tied predictor a chance to update their coefficients in a balanced way, while remaining tied in correlation. However, he then realized that the LAR least-squares fit amongst the tied predictors can result in coefficients moving in the opposite direction to their correlation, which cannot happen in Algorithm 3.4. The following modification of the LAR algorithm implements \(FS_0\):

Algorithm 3.2b Least Angle Regression: \(FS_0\) Modification.

4. Find the new direction by solving the constrained least squares problem \[\underset{b}{\text{min}}\lVert r-X_Ab\rVert_2^2\text{ subject to } b_js_j\ge 0,j\in A\] where \(s_j\) is the sign of \(\langle x_j, r\rangle\). The modification amounts to a non-negative least squares fit, keeping the signs of the coefficients the same as those of the correlations.