LASSO原作者的論文，來讀讀看

Regression Shrinkage and Selection via the lasso

眾所周知，Robert Tibshirani是統計領域的大佬，這篇文章在1996年提出了LASSO，之后風靡整個高維領域，並延伸出許多種模型。這篇文章截止2019.5.16已經獲得了27991的引用量（跪下）。

雖然LASSO是非常直觀且大家都很熟悉的模型，但重溫經典也無不可。了解一個模型就去讀原作者的文章，獲得的信息是最沒有損失的。

Background introduction

在回歸模型的場景中，我們使用最多的是最小二乘法估計參數（OLS）。

但OLS有很強的局限性，具體表現在：

• 預測精度（prediction accuracy）

  OLS估計量出於最小化均方誤差的目的，通常給出的估計結果有很低的偏差（擬合訓練集的誤差）但是有很高的方差（模型的泛化能力）

• 解釋能力（interpretation）

  保留大量的解釋變量會降低模型的可解釋性，因此需要找到一個具有很強影響力的解釋變量子集。

為了解決上述問題，目前有兩種方法：

• subset selection

  Cons：通過篩選變量增加了模型的解釋能力。

  Pron：但是模型並不穩健，因為選擇變量是一個離散的過程，變量只會面臨進入或退出兩種可能。數據的微小改變會導致模型非常大的改變，因此降低預測的准確性。

• ridge regression

  Cons：是一個連續的變量系數縮減過程，因此較為穩定。

  Pron：沒有設置任何變量的系數為0，不能形成解釋力強的模型。

基於此，作者提出LASSO模型，全稱為'least absolute shrinkage and selection operator'（最小絕對收斂和選擇算子），能夠同時保留兩種方法的優點。

一個類似的模型是Breiman於1993年提出的non-negative garotte.

\[(\hat{\alpha},\hat{\beta})=\arg min{\sum_{i=1}^N(y_i-\alpha-\sum_jc_j\hat{\beta_j}^ox_{ij})^2} \\subject \quad to\quad c_j\ge0\quad\sum_jc_j\le t \]

個人感覺這個模型是通過允許縮減系數的存在但限制縮減系數的大小來實現對傳統回歸模型的優化。該模型相較subset selection能夠減少預測誤差，效果和ridge regression差不多。但該模型直接使用了OLS估計，如果OLS估計量本身表現差，那么這個模型也會表現差。

LASSO避免了對OLS估計量的直接使用。

Basic ideas

接下來，作者以標准化的解釋變量為例。set \(x^i=(x_{i1},\dots,x_{ip})^T\) ,where \(\sum_ix_{ij}/N=0,\sum_ix_{ij}^2/N=1\), i.e, \(X^T\cdot X=I\) .

參數向量\(\hat{\beta}=(\hat{\beta_1},\cdots,\hat{\beta_p})^T\),lasso估計量通過求解下列優化問題得到：

\[(\hat{\alpha},\hat{\beta})=\arg min{\sum_{i=1}^N(y_i-\alpha-\sum_j\beta_jx_{ij})^2} \\subject \quad to\quad \sum_j|{\beta_j}|\le t \tag{1} \]

同時假設\(\bar{y}=0\),因此可以省略\(\alpha\)。估計最優系數不需要設計矩陣是滿秩的，這一點對於高維問題有很好的適應性。

參數估計及解釋（一元➡️二元）

方程（1）的解為

\[\hat{\beta_j}=sign(\hat{\beta^0_j})(|\hat{\beta^0_j|}-\gamma)^+ \]

Ridge Regression	Garotte Method
\(\frac{1}{1+\gamma}\hat{\beta_j^o}\)	\((1-\frac{\gamma}{\hat{\beta_j^{o2}}})^+\hat{\beta_j^o}\)

接下來從圖形上展示lasso方法，使用二維

空間p=2。由於lasso的限制條件是有棱角的，因此更容易切在某一維度的參數取零的點。

將變量個數拓展為兩個后，對系數的估計結果為（不失一般性，假設OLS估計的系數為正）：

\[\begin{array}{c}{\hat{\beta}_{1}=\left(\frac{t}{2}+\frac{\hat{\beta}_{1}^{\circ}-\hat{\beta}_{2}^{\circ}}{2}\right)^{+}} \\ {\hat{\beta}_{2}=\left(\frac{t}{2}-\frac{\hat{\beta}_{1}^{\circ}-\hat{\beta}_{2}^{\circ}}{2}\right)^{+}}\end{array} \]

系數的估計對變量之間的相關系數並不敏感。

然而，在二元情況下，嶺回歸的估計結果取決於變量之間的相關程度。當相關系數為0時，嶺回歸有縮減系數的效果，但隨着變量的相關性增強，系數縮減效果越來越低，甚至在增加。（后面的simulation中可以看到）

As pointed out by Jerome Friedman, this is due to the tendency of ridge regression to try to make the coefficients equal to minimize their square norm.

利用bootstrap估計參數標准差

首先將懲罰項寫成\(\Sigma \beta_{j}^{2} /\left|\beta_{j}\right|\)的形式，然后借助嶺回歸來估計系數的協方差矩陣：

\[ \left(\mathbf{X}^{\mathrm{T}} \mathbf{X}+\lambda \mathbf{W}^{-}\right)^{-1} \mathbf{X}^{\mathrm{T}} \mathbf{X}\left(\mathbf{X}^{\mathrm{T}} \mathbf{X}+\lambda \mathbf{X}^{-}\right)^{-1} \hat{\boldsymbol{\sigma}}^{2} \]

此外，還有一個通過迭代嶺回歸求解的方式，但在算法效率上較低。不過這一方式可以用於選擇lasso parameter t。

實證——prostate cancer data

x (after standardized)

log(cancer volume)

log(prostate weight)

age

log(benign prostatic hyperplasia amount)

seminal vesicle invasion

log(capsular penetration)

Gleason score

percentage Gleason scores 4 or 5

需要預測的y為log(prostate specific antigen)【前列腺特定抗原】

Lasso shrinkage of coefficients in the prostate cancer example: each curve represents a
coefficient (labelled on the right) as a function of the (scaled) lasso parameter \(s=t / \Sigma\left|\beta_{j}^{\circ}\right|\) (the intercept is not plotted); the broken line represents the model for \(\hat{s}\) = 0.44, selected by generalized cross-validation.

上圖是設定了s的閾值后，各參數的縮減程度曲線。對於嶺回歸和subset regression，變量的選擇不是單調的，即最后選出的5個最優變量和選出的4個最優變量不是包含關系。 \(\hat{s}\) = 0.44時，相當於我們只選擇了一半的變量。

Notice that the coefficients and Z-scores for the selected predictors from subset selection tend to be larger than the full model values: this is common with positively correlated predictors. However, the lasso shows the opposite effect, as it shrinks the coefficients and Z-scores from their full model values.

對於subset selection，參數估計值通常要高於OLS估計（這種情況在變量有正相關時較為普遍）。但LASSO仍然能夠保持縮減系數的效果。前面提到過，對於0系數標准差的估計，通過嶺回歸近似公式不能計算，但可以用bootstrap代替。

預測誤差和最優超參數t的選擇

作者使用了三種方法估計lasso parameter：cross-validation, generalized cross-validation and an analytical unbiased estimate of risk。前兩種適用於x為隨機變量的情形，后一種適用於x固定的情形。在實際應用中，兩種場景差別不大。

直接給出作者推導的結果。

1⃣️通過交叉驗證計算

給出PE（預測誤差）：

\[\mathrm{PE}=E\{Y-\hat{\eta}(\mathbf{X})\}^{2}=\mathbf{M E}+\sigma^{2} \]

作者采用了五折交叉驗證，選出一個最優的\(\hat{s}\)，使得PE最低。

2⃣️通過lasso估計量的線性估計來推導t。(需要假設正交的design matrix)

\[p(t)=\operatorname{tr}\left\{\mathbf{X}\left(\mathbf{X}^{\mathbf{T}} \mathbf{X}+\lambda \mathbf{W}^{-}\right)^{-1} \mathbf{X}^{\mathbf{T}}\right\}\\ \mathrm{GCV}(t)=\frac{1}{N} \frac{\mathrm{rss}(t)}{\{1-p(t) / N\}^{2}}\\ E_{\mu}\|\hat{\mu}-\mu\|^{2}=p+E_{\mu}\left(\|\mathbf{g}(\mathbf{z})\|^{2}+2 \sum_{1}^{p} \mathrm{d} g_{i} / \mathrm{d} z_{i}\right)\\ R\{\hat{\boldsymbol{\beta}}(\gamma)\} \approx \hat{\tau}^{2}\left\{p-2 \#\left(j ;\left|\hat{\beta}_{j}^{\circ} / \hat{\tau}\right|<\gamma\right)+\sum_{j=1}^{p} \max \left(\left|\hat{\beta}_{j}^{\circ} / \hat{\tau}\right|, \gamma\right)^{2}\right\}\\ \hat{\gamma}=\arg \min _{\gamma \geqslant 0}[R\{\hat{\boldsymbol{\beta}}(\gamma)\}]\\ \hat{t}=\sum\left(\left|\hat{\beta}_{j}^{\circ}\right|-\hat{\gamma}\right)^{+} \]

第二種方法在計算效率上遠勝第一種方法。generalized cross-validation介於上述兩種方法之間。

當參數\(\beta_j\)具有獨立雙指數分布的先驗假設時，也可以推導出LASSO估計量。由於指數分布的厚尾性質，可以在0點附近給更大的概率值。

P維情況下求解lasso估計量

Lawson and Hansen (1974) 提供了解決有不等式約束的線性最小二乘問題的方法。但對LASSO模型來說，\(m=2^p\)個約束條件太多了，實施起來很困難。但我們可以選擇將不等式約束逐一加入，尋找滿足KKT條件的解。

令\(g(\boldsymbol{\beta})=\boldsymbol{\Sigma}_{i=1}^{N}\left(y_{i}-\boldsymbol{\Sigma}_{j} \boldsymbol{\beta}_{j} \boldsymbol{x}_{i j}\right)^{2}\) ，將\(\Sigma\left|\beta_{j}\right| \leqslant t\) 拆成\(2^p\)個去掉絕對值的不等式約束\(\boldsymbol{\delta}_{i}^{\mathrm{T}} \boldsymbol{\beta} \leqslant t\)。

對於給定的\(\beta\)，設置等式集E和不等式集S:\(E=\left\{i : \boldsymbol{\delta}_{i}^{\mathrm{T}} \boldsymbol{\beta}=t\right\} \text { and } S=\left\{i : \delta_{i}^{\mathrm{T}} \boldsymbol{\beta}<t\right\}\) .

\(\text { Denote by } G_{E} \text { the matrix whose rows are } \delta_{i} \text { for } i \in E\).

Algorithm outline：

(a) Start with \(E=\left\{i_{0}\right\}\)where \(\boldsymbol{\delta {i{0}}=\operatorname{sign}\left(\hat{\boldsymbol{\beta}}^{\circ}\right), \hat{\boldsymbol{\beta}}^{\circ}}\) being the overall least squares estimate.

(b) Find $ \hat{\beta}$ to minimize \(g(\boldsymbol{\beta})\)subject to \(G_{E} \boldsymbol{\beta} \leqslant t \mathbf{1}\).（1是1向量）#添加第一個不等式，求解

(c) While \(\left\{\Sigma\left|\hat{\beta}_{j}\right|>t\right\}\).#檢查最終的條件

(d) add i to the set E where \(\boldsymbol{\delta}_{i}=\operatorname{sign}(\hat{\boldsymbol{\beta}})\). Find \(\hat{\beta}\) to minimize \(g(\boldsymbol{\beta})\)subject to

\(G_{E} \boldsymbol{\beta} \leqslant t \mathbf{1}\).#添加第二個不等式，並求解

該算法的收斂特性：

This procedure must always converge in a finite number of steps since one element
is added to the set E a t each step, and there is a total of 2p elements.

The average number of iterations required is in the range (0.5p, 0.75p).

另一種完全不同的算法由David Gay提出。將原始問題轉化為雖然具有2p個變量但限制條件減少為2p+1的問題。這種算法得到的解和原始問題的解一樣。

\(\begin{array}{l}{\text { We write each } \beta_{j} \text { as } \beta_{j}^{+}-\beta_{j}^{-}, \text { where } \beta_{j}^{+} \text { and } \beta_{j}^{+} \text { are non-negative. Then we solve the }} \\ {\text { least squares problem with the constraints } \beta_{j}^{+} \geqslant 0, \beta_{j}^{-} \geqslant 0 \text { and } \Sigma \beta_{j}^{+}+\Sigma_{j} \beta_{j}^{-} \leqslant t . }\end{array}\)

標准的二次規划也可用於解決該問題，2p+1步內能保證收斂。

Standard quadratic programming techniques can be applied, with the convergence
assured in 2p +1 steps.

Simulation

目的：比較普通最小二乘和lasso、non-negative garotte，best subset selection and ridge regression. 在每個case中都用五折驗證法選擇參數。

\[\operatorname{PE}(J)=\frac{1}{5} \sum_{v=1}^{5} \operatorname{PE}^{v}(J) \]

Zhang(1993) and Shao(1992)指出，如果只選擇誤差最小的模型，會導致估計的不一致性。

Example 1

估計了50個data set，每個data set都包含20個觀測值。

DGP為\(y=\boldsymbol{\beta}^{\mathrm{T}} \mathbf{x}+\sigma \epsilon\),where \(\boldsymbol{\beta}=(3,1.5,0,0,2,0,0,0)^{\mathrm{T}}\) and \(\epsilon\) 服從標准正態分布。變量之間的相關系數為\(\rho^{|i-j|}\)，這里\(\rho\)設置為0.5，\(\sigma=3\).

這里可以看出，系數中有5個無意義。signal-to-noise ratio大約為5.7。

在200個模擬中，lasso效果最好，其次是garotte和嶺回歸。subset selection正確篩選掉了0系數，但是它的估計結果很不穩健。

同時，對比lasso和subset regression可以發現，雖然lasso完全估計出正確的模型（1，2，5）的頻率很少，為2.5%，但它給出的模型包含正確參數的頻率高達95%。而subset regression盡管給出完全正確模型的頻率為24%，但所有模型中包含正確參數的頻率僅為53.5%。

Example 2

基本假設與example1一樣，但是參數改為\(\beta_j=0.85\)。signal-to-noise ratio大約為1.8.最終結果顯示LASSO和嶺回歸的在MSE意義上優於OLS，而其他幾種方法效果都較差。LASSO的效果在這個example中不如嶺回歸，原因是DGP的參數不再包含0參數，而lasso傾向於篩選掉一些變量，因此相比嶺回歸存在更大的信息損失。

Example 3

作者挑選了一個非常適合做subset selection的數據集，數據生成和example1一樣，但設置\(\boldsymbol{\beta}=(5,0,0,0,0,0,0,0)\) ，\(\sigma=2\),因此signal-to-noise ratio大約為7.

結果顯示garotte和subset selection效果最好，其次是lasso。嶺回歸的效果不佳，甚至不如OLS。

Example 4

通過模擬一個更大的模型來檢查lasso的表現。作者模擬了50個數據集，每個都有100個觀測值和40個變量。DGP為：

\(x_{i j}=z_{i j}+z_{i} \text { where } z_{i j} \text { and } z_{i} \text { are independent standard normal variates. }\)

\(\boldsymbol{\beta}=(0,0, \ldots, 0,2,2, \ldots, 2,0,0, \ldots, 0,2,2, \ldots, 2)\)，每個block重復了10次，\(\sigma=15\).

signal-to-noise ratio大約為9. 結果顯示，嶺回歸效果最好，lasso其次。

Application to Generalized Regression models

IRLS: models a quadratic approximation to \(\ l (\beta)\) and leads to an iteratively reweighed least squares procedure for computation of \(\beta\) .和牛頓迭代原理一樣。

Logistic regression

使用駝背數據kyphosis data（ analysed in Hastie and Tibshirani (1990)）。

description：

0=absent，1=present。

x1=age，x2=number of vertebrae levels，x3=starting vertebrae level.

number of oberservations=83.

由於解釋變量的影響不是線性的（可能影響越來越大or小，醫學上的經驗？）。作者還在模型中加入了二次項。變量都經過了標准化。

logistic regression：

\[-2.64+0.83 x_{1}+0.77 x_{2}-2.28 x_{3}-1.55 x_{1}^{2}+0.03 x_{2}^{2}-1.17 x_{3}^{2} \]

backward stepwise deletion：

\[-2.64+0.84 x_{1}+0.80 x_{2}-2.28 x_{3}-1.54 x_{1}^{2}-1.16 x_{3}^{2} \]

Lasso:

\[-1.42+0.03 x_{1}+0.31 x_{2}-0.48 x_{3}-0.28 x_{1}^{2} \]

大約5次后收斂。

Results on soft thresholding

lasso：\(\hat{\beta}_{j}=\operatorname{sign}\left(\hat{\beta}_{j}^{\circ}\right)\left(\left|\hat{\beta}_{j}^{\circ}\right|-\gamma\right)^{+}\)軟閾值。

subset selection:\(\tilde{\beta}_{j}=\hat{\beta}_{j}^{\circ} I\left(\left|\hat{\beta}_{j}^{\circ}\right|>\gamma\right)\)硬閾值。

可以求出Risk的下界。

Conclusion

LASSO的優點：

• 連續選擇變量
• 稀疏模型+變量篩選
• 效果優於subset selection和ridge regression
• 具有可解釋性

另外提到lasso是bridge（使用了Lp正則）的特例。