在周志华老师的《机器学习》中的线性回归章节内，通过简单的微积分给出了最小二乘
\[\min_{\pmb w}\quad\Vert\pmb{y}-\pmb{Xw}\Vert_2^2\]的解析解，用类似的方法，我们也可以求出带$L_2$正则化项的最小二乘
\[\min_{\pmb w}\quad\Vert\pmb{y}-\pmb{Xw}\Vert_2^2+\lambda\Vert\pmb{w}\Vert_2^2\]的解，上面两个模型对应sklearn.linear_model中的LinearRegression和RidgeRegression，但线性回归模型中，还存在以$L_1$正则化项的最小二乘，也就是LassoRegression：
由于后面一项不可微，导致我们无法用简单的求导去求解。笔者当时虽然有实现它的想法，但囿于能力水平的限制，未能付诸实施。本文将介绍LassoRegression的实现方法。
近端梯度下降法用于求解
\[\min_{x}\quad f(x)=g(x)+h(x)\]其中$g$是可微凸函数，$h$是不可微凸函数，可以发现*式符合该问题的形式。此外，$\nabla f$需满足$L$-Lipschitz条件，也就是存在正常数$L$，对任意的$\pmb x_1,\pmb x_2$，有
\[\Vert\nabla f(\pmb x_1)-\nabla f(\pmb x_2)\Vert_2^2\leq L\Vert\pmb x_1-\pmb x_2\Vert_2^2\]我们来证明线性回归下的均方误差函数满足该条件：定义
\[F(\pmb x_1,\pmb x_2)=\Vert\nabla f(\pmb x_1)-\nabla f(\pmb x_2)\Vert_2^2-L\Vert\pmb x_1-\pmb x_2\Vert_2^2\]那么我们有
\[\begin{aligned} F(\pmb w_1,\pmb w_2)&=\bigg\Vert \frac1n\pmb X^T(\pmb{Xw}_1-\pmb{y})-\frac1n\pmb X^T(\pmb{Xw}_2-\pmb{y}) \bigg\Vert_2^2-L\Vert\pmb{w}_1-\pmb{w}_2\Vert_2^2\\ &=\dfrac{1}{n^2}\bigg\Vert\pmb{X}^T\pmb X(\pmb w_1-\pmb w_2)\bigg\Vert_2^2-L\Vert\pmb{w}_1-\pmb{w}_2\Vert_2^2\\ &=L\Vert\pmb{w}_1-\pmb{w}_2\Vert_2^2\bigg(\dfrac{\Vert\pmb{X}^T\pmb X(\pmb w_1-\pmb w_2)\Vert_2^2}{n^2L\Vert\pmb{w}_1-\pmb{w}_2\Vert_2^2}-1\bigg) \end{aligned}\]只要证明
\[\dfrac{\Vert\pmb{X}^T\pmb X(\pmb w_1-\pmb w_2)\Vert_2^2}{\Vert\pmb{w}_1-\pmb{w}_2\Vert_2^2}\]是有界的，那么只要$L$足够大，$F(\pmb w_1,\pmb w_2)$就会非负。现在证有界，实际上是证明，对任意的$\pmb x$和实对称矩阵$A$，我们总有
\[\dfrac{\Vert A\pmb x\Vert_2^2}{\Vert\pmb x\Vert_2^2}\]有界。由实对称矩阵的特征值性质，我们可以将$\pmb x$分解为$A$的所有特征向量的线性组合，因为这些特征向量都是正交的：
\[\begin{aligned} \dfrac{\Vert A\pmb x\Vert_2^2}{\Vert\pmb x\Vert_2^2} &=\dfrac{\Vert A\sum_{i=1}^n k_i\pmb v_i\Vert_2^2}{\Vert\sum_{i=1}^n k_i\pmb v_i\Vert_2^2}\\ &=\dfrac{\Vert \sum_{i=1}^n k_i A\pmb v_i\Vert_2^2}{\Vert\sum_{i=1}^n k_i\pmb v_i\Vert_2^2}\\ &=\dfrac{\Vert \sum_{i=1}^n k_i\lambda_i\pmb v_i\Vert_2^2}{\Vert\sum_{i=1}^n k_i\pmb v_i\Vert_2^2}\\ &=\dfrac{\sum_{i=1}^n k_i^2\lambda_i^2\Vert v_i\Vert_2^2}{\sum_{i=1}^n k_i^2\Vert v_i\Vert_2^2}\\ &=\dfrac{\sum_{i=1}^n k_i^2\lambda_i^2}{\sum_{i=1}^n k_i^2}\\ &\leq \max_i(\lambda_i^2) \end{aligned}\]因此有界，此外我们可以求符合条件的$L$:
\[\begin{aligned} F(\pmb w_1,\pmb w_2) &=\dfrac{\Vert\pmb{X}^T\pmb X(\pmb w_1-\pmb w_2)\Vert_2^2}{n^2L\Vert\pmb{w}_1-\pmb{w}_2\Vert_2^2}-1\\ &\leq\dfrac{\max_i(\lambda_i^2)}{n^2L}-1\\ &\leq0\\ L&\geq\dfrac{\max_i(\lambda_i^2)}{n^2} \end{aligned}\]其中$\lambda_i$是矩阵$\pmb X^T\pmb X$的特征值。
在函数具备$L$-Lipschitz条件后，我们对函数进行泰勒展开：
\[\begin{aligned} \hat{f}(\pmb x)&\approx f(\pmb x_k)+\nabla f(\pmb x_k)^T(\pmb x-\pmb x_k)+\dfrac{L}{2}\Vert\pmb x-\pmb x_k\Vert_2^2\\ &=\dfrac{L}{2}\bigg\Vert\pmb{x}-\big(\pmb x_k-\dfrac1L\nabla f(\pmb x_k) \big) \bigg\Vert_2^2+\text{const} \end{aligned}\]是一个二次函数，令
\[\pmb{x}=\pmb x_k-\dfrac1L\nabla f(\pmb x_k)\]可以实现最大程度的最小化，这也是梯度下降的一种理解方式，在这里我们也可以看到，$\frac1L$是一个合理的步长。将上式推广到Lasso问题，就可以得到迭代式：
\[\pmb x_{k+1}=\mathop{\arg\min}\limits_{\pmb x}\dfrac{L}{2}\bigg\Vert\pmb{x}-\big(\pmb x_k-\dfrac1L\nabla f(\pmb x_k) \big) \bigg\Vert_2^2+\lambda\Vert\pmb x\Vert_1\]将上式重组一下：
\[\begin{aligned} \dfrac{L}{2}\bigg\Vert\pmb{x}-\big(\pmb x_k-\dfrac1L\nabla f(\pmb x_k) \big) \bigg\Vert_2^2+\lambda\Vert\pmb x\Vert_1&=\frac L2\sum_{i=1}^n\bigg(x^{(i)}-(x_k^{(i)}-\frac1L\nabla_i f(\pmb x_k))\bigg)^2+\lambda\sum_{i=1}^n\vert x_i\vert \end{aligned}\]因此，对参数$\pmb x$的每一个维度分别去优化，结果就是最优的，也就是求
\[\mathop{\arg\min}\limits_x\dfrac L2(x-z)^2+\lambda\vert x\vert\]上式的求解可以用分类讨论，但我们这里通过次梯度进行求解。首先看次梯度，对于一可微凸函数$f$，我们总有（一阶特性）
\[f(y)\ge f(x)+g^T(y-x)\]但当函数不可微时，我们也可以根据该性质来定义其次梯度。一个凸函数$f$在$x$的次梯度$g$定义为：
\[f(y)\ge f(x)+g^T(y-x)\]如果$x$是$f$上的完全可微点，其次梯度只有一个，就是其梯度。以$f(x)=\vert x\vert$为例，其在$x=0$处不可微，那么根据次梯度定义：
\[f(y)=\vert y\vert\ge 0+gy=gy\]解得$g\in[-1,1]$，也就是函数在该点的次梯度为$[-1,1]$，而在其他点，函数在其上都是可微的，所以次梯度就是其梯度，为-1或1.
再引入次微分的概念，微分$\partial f$是一个确定值，而凸函数$f$的次微分是其所有次梯度的集合。次微分和次梯度将函数的最优条件进行扩展：对于任意$f$，都有
\[f(x^*)=\min_x f(x)\leftrightarrows0\in\partial f(x^*)\]还是以绝对值函数为例，在$x<0,x=0,x>0$三种情况下，其次梯度分别是
\[\{-1\},[-1,1],\{1\}\]因此，只有$x=0$才是最优点。我们回到前面的问题，要找到
\[h(z)=\mathop{\arg\min}\limits_x\dfrac L2(x-z)^2+\lambda\vert x\vert\]其中$\lambda>0$，那么就有
\[\begin{aligned} 0&\in\partial \dfrac L2(x-z)^2+\lambda\vert x\vert\\ 0&\in L(x-z)+\lambda\partial\vert x\vert\\ &\to\begin{cases} L(x-z)+\lambda\cdot\text{sgn}(x)=0&x\neq 0\\ \vert L(x-z)\vert\leq\lambda&x=0 \end{cases} \end{aligned}\]分类讨论，当$x<0$，此时最优的$x$：
\[\begin{aligned} L(x-z)-\lambda&=0\\ x&=z+\dfrac{\lambda}{L} \end{aligned}\]显然要满足
\[z+\dfrac{\lambda}{L}<0\\ z<-\dfrac{\lambda}{L}\]类似的，当$x>0$，最优的$x$:
\[x=z-\dfrac{\lambda}{L}\]且满足
\[z>\dfrac{\lambda}{L}\]当$x=0$，则需要满足
\[\begin{aligned} -\lambda\leq&-Lz\leq\lambda\\ \vert z\vert&\leq\frac{\lambda}L \end{aligned}\]因此，在计算出$z$后，我们就可以利用上面的规则进行更新。
我们可以将PGD算法封装成一个优化器，从而实现端到端（end to end）的优化过程。
class PGD:
def __init__(self, *params, lr=0.1, alpha=1):
self.lr = lr
self.params = list(params)
self.n_params = len(self.params)
self.bound = self.lr * alpha
def step(self, *grad):
assert len(grad) == len(self.params)
for i in range(self.n_params):
z = self.params[i] - self.lr * grad[i]
self.params[i][np.abs(z) <= self.bound] = 0
self.params[i][z > self.bound] = z[z > self.bound] - self.bound
self.params[i][z < -self.bound] = z[z < -self.bound] + self.bound
def get_params(self):
return self.params
由此，我们就可以用PGD算法去求解Lasso回归，以及一系列使用$L_1$正则化的优化问题。我们这里使用sklearn中的数据集进行测试：
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.datasets import load_boston
from sklearn.metrics import mean_squared_error
# 省略PGD optimizer的定义
X, y = load_boston(return_X_y=True)
X = (X - X.mean(0)) / X.std(0) # 标准化，否则梯度会overflow
X = np.hstack((X, np.ones((len(X), 1)))) # 右边加上全1列，因为要考虑偏置b
train_X, test_X, train_y, test_y = train_test_split(X, y, train_size=0.7)
lr, n_iter, alpha = 0.1, 1000, 1
w = np.zeros(X.shape[1]) # 初始化参数
optim = PGD(
w,
lr=lr,
alpha=alpha,
) # 定义优化器
error = mean_squared_error(train_y, train_X @ w) # 初始误差
for i in range(n_iter):
grad = train_X.T @ (train_X @ w - train_y) / len(train_X) # 计算梯度
optim.step(grad) # 优化
# 如果误差下降过小，选择退出迭代
new_error = mean_squared_error(train_y, train_X @ w)
if abs(error - new_error) < 1e-8:
break
error = new_error
test_error = mean_squared_error(test_y, test_X @ w)
print("迭代次数 : {}\n训练误差 : {}\n测试误差 : {}\n参数 :\n{}".format(
i,
new_error,
test_error,
w,
))
输出：
迭代次数 : 414
训练误差 : 28.503625887591152
测试误差 : 29.24557875464323
参数 :
[ 0. 0. 0. 0.06968633 0. 2.87566021
0. 0. 0. 0. -1.29345419 0.04994334
-3.74374931 21.71602387]
可以发现$L_1$范数正则化下，参数的很多项都是0，意味着这些项对应的特征对标签值的没有贡献或贡献极小，有效帮助我们进行特征选择。
