LASSO的解法

Posted 2021-06-17 分析101

LASSO非常实用，但由于它的惩罚项不可以常规地进行求导，使得很多人以为它无法显式地求出解析解。但其实并不是这样的。

1 单变量情形：软阈值法

1.1 软阈值的分类讨论

将\\(N\\)个样本的真实值记为\\(N\\)维向量\\(y\\)，将\\(N\\)个样本的自变量记为\\(z\\)，假设我们已经将自变量做过标准化，即\\(z\' \\ell_n=0\\)，\\(z\'z/N=1\\)，这也意味着在LASSO模型中截距项为\\(0\\)。系数\\(\\beta\\)是要优化的参数，惩罚项参数为\\(\\lambda\\gt 0\\)。

LASSO就是要求解

\\[\\argmin_\\beta \\dfrac{1}{2N}(y-z\\beta)\'(y-z\\beta)+\\lambda |\\beta| \\tag{1} \\]

忽略常数项后，上式等价于

\\[\\argmin_\\beta \\dfrac{1}{2}\\beta^2 -\\dfrac{y\'z}{N}\\beta+ \\lambda |\\beta| \\]

将损失函数写成分段函数形式：

\\[f(\\beta)=\\begin{cases} f_1(\\beta) = \\dfrac{1}{2}\\beta^2 -\\left(\\dfrac{y\'z}{N} + \\lambda\\right)\\beta , \\beta \\lt 0\\\\ f_2(\\beta) = \\dfrac{1}{2}\\beta^2 -\\left(\\dfrac{y\'z}{N}- \\lambda\\right)\\beta, \\beta \\geq 0 \\end{cases} \\]

分类讨论：

• 若\\(\\dfrac{y\'z}{N}\\gt \\lambda\\)，则\\(f_1(\\beta) \\gt 0\\)，\\(f_2(\\beta)\\)在\\(\\hat\\beta=\\dfrac{y\'z}{N}- \\lambda\\)处取到最小值\\(f_2(\\hat\\beta)\\lt 0\\)，因此解为\\(\\hat\\beta=\\dfrac{y\'z}{N}- \\lambda\\)；
• 若\\(\\left|\\dfrac{y\'z}{N}\\right|\\leq \\lambda\\)，则\\(f_1(\\beta) \\geq 0\\)，\\(f_2(\\beta) \\geq 0\\)，且在\\(\\hat\\beta=0\\)处有\\(f_1(\\hat\\beta)=f_2(\\hat\\beta)=0\\)，因此解为\\(\\hat\\beta=0\\)；
• 若\\(\\dfrac{y\'z}{N}\\lt -\\lambda\\)，则\\(f_2(\\beta) \\gt 0\\)，\\(f_1(\\beta)\\)在\\(\\hat\\beta=\\dfrac{y\'z}{N}+\\lambda\\)处取到最小值\\(f_1(\\hat\\beta)\\lt 0\\)，因此解为\\(\\hat\\beta=\\dfrac{y\'z}{N}+\\lambda\\)。

利用软阈值算子（soft-thresholding operator）\\(S_\\lambda(x)=\\text{sign}(x)(|x|-\\lambda)_+\\)，可将以上三种解统一为

\\[\\hat\\beta=S_\\lambda(y\'z/N) \\]

其实在我们的设定下，OLS估计量为\\(\\tilde\\beta=y\'z/N\\)，因此，将OLS估计量通过一个软阈值算子的操作，就变成了LASSO估计量。

1.2 次梯度

如果引入次梯度（subgradient）的概念，可以更直接地求解\\((1)\\)式。设\\(|\\beta|\\)的次梯度为\\(s\\in \\text{sign}(\\beta)\\)，它的形式是，当\\(\\beta \\neq 0\\)时有\\(s= \\text{sign}(\\beta)\\)，当\\(\\beta = 0\\)时有\\(s\\in [-1,1]\\)。根据凸优化（convex optimization）理论，求解\\((1)\\)相当于求解

\\[-\\dfrac{1}{N}z\'(y-z\\beta) +\\lambda s =0 \\]

的解\\((\\hat\\beta,\\hat\\lambda)\\)。化简后得到\\(y\'z/N = \\beta+\\lambda s\\in\\beta+\\lambda \\text{sign}(\\beta)\\)，最终同样可以解出\\(\\hat\\beta=S_\\lambda(y\'z/N)\\)。比如\\(\\beta=0\\)时，就意味着\\(y\'z/N \\in[-\\lambda,\\lambda]\\)。

2 多变量情形：循环坐标下降法

我们来看多变量的完整版LASSO问题。将自变量排成\\(N\\times p\\)的矩阵\\(X\\)，我们要求解的是

\\[\\argmin_{\\beta\\in \\mathbb{R}^p} \\dfrac{1}{2N}\\left\\Vert y-X\\beta\\Vert_2^2+\\lambda \\right\\Vert\\beta\\Vert_1 \\]

在这里，我们使用循环坐标下降法（cyclic coordinate descent），它的思想是，按一定顺序依次对\\(p\\)个参数进行优化，比如按\\(j=1,\\ldots,p\\)的顺序，在第\\(j\\)次优化时，保持其他所有系数不变，变动\\(\\beta_j\\)使损失函数最小化。

根据以上思想，我们将第\\(j\\)次的最优化目标写为

\\[\\argmin_{\\beta_j\\in \\mathbb{R}^p} \\dfrac{1}{2N}\\left\\Vert y-\\sum_{k\\neq j}x_{\\cdot k}\\beta_k-x_{\\cdot j}\\beta_j\\right\\Vert^2_2+\\lambda \\sum_{k\\neq j}|\\beta_k| +\\lambda |\\beta_j| \\]

记\\(r^{(j)}=y-\\sum_{k\\neq j}x_{\\cdot k}\\hat{\\beta}_k\\)，这称为partial residual，那么根据第1.1节中的结果，我们可以得出

\\[\\hat\\beta_j = S_\\lambda (x_{\\cdot j}\' r^{(j)}/N) \\]

记\\(r = r^{(j)}-x_{\\cdot j}\\hat\\beta_j\\)，上式相当于更新规则

\\[\\hat\\beta_j \\leftarrow S_\\lambda (\\hat\\beta_j +x_{\\cdot j}\' r/N) \\]

由于目标函数是凸的，没有局部的极小值，因此，每次更新一个坐标，最终可以收敛到全局最优解。

Pathwise coordinate descent（逐路径坐标下降）：可以先取一个使\\(\\hat\\beta=0\\)的最小的\\(\\lambda\\)，然后，略微减小\\(\\lambda\\)的值，并以上一次得到的\\(\\hat\\beta\\)作为“warm start”，用坐标下降法迭代直到收敛，不断重复这个过程，最终就可以得到在\\(\\lambda\\)的一系列变化范围内的解。

那么，怎样才能使\\(\\hat\\beta=0\\)？利用次梯度，我们可以知道，对于\\(\\hat\\beta_j=0\\)，必有\\(x_{\\cdot j}\'y /N \\in [-\\lambda,\\lambda]\\)，即要求\\(\\lambda \\geq |x_{\\cdot j}\'y| /N\\)，若要使整个\\(\\hat\\beta=0\\)，则可取\\(\\lambda =\\max_j |x_{\\cdot j}\'y| /N\\)，这就是使\\(\\hat\\beta=0\\)的最小的\\(\\lambda\\)。

3 其他解法

求解LASSO还有其他的解法，如homotopy method，它可以从\\(0\\)开始，得到序列型的解的路径，路径是分段线性的。

还有LARS（least angle regression）算法，这是homotopy method之一，可以有效得到分段线性路径。

这里不作展开。

4 正交基

在上面的过程中，如果将自变量正交化，可以大大简化计算。如在第2节中，如果自变量之间是正交的，则\\(x_{\\cdot j}\' r^{(j)}/N = x_{\\cdot j}\' y/N\\)，此时\\(\\hat\\beta_j\\)就是将\\(y\\)对\\(x_{\\cdot j}\\)做回归的OLS解，通过软阈值算子后的值。