浏览本文须要的背景知识点：线性回归算法、一丢丢编程常识
一、引言
上一节咱们学习了解决多重共线性的一种办法是对代价函数正则化，其中一种正则化的算法叫岭回归算法（Ridge Regression Algorithm）。上面咱们来学习另一种正则化的算法 - Lasso回归算法)1（Lasso Regression Algorithm），LASSO的残缺名称叫最小绝对值收敛和抉择算子算法（least absolute shrinkage and selection operator）。
二、模型介绍
先来回顾一下岭回归的代价函数，在原来规范线性回归代价函数上加上了一个带惩办系数 的 w 向量的L2-范数的平方：
$$\operatorname{Cost}(w) = \sum_{i = 1}^N \left( y_i - w^Tx_i \right)^2 + \lambda\|w\|_{2}^{2}$$
Lasso回归算法也同岭回归一样加上了正则项，只是改成加上了一个带惩办系数 的 w 向量的L1-范数作为惩办项（L1-范数的含意为向量 w 每个元素绝对值的和），所以这种正则化形式也被称为L1正则化。
$$\operatorname{Cost}(w) = \sum_{i = 1}^N \left( y_i - w^Tx_i \right)^2 + \lambda\|w\|_1$$
同样是求使得代价函数最小时 w 的大小：
$$w = \underset{w}{\operatorname{argmin}}\left(\sum_{i = 1}^N \left( y_i - w^Tx_i \right)^2 + \lambda\|w\|_1\right)$$
因为退出的是向量的L1-范数，其中存在绝对值，导致其代价函数不是处处可导的，所以就没方法通过间接求导的形式来间接失去 w 的解析解。上面介绍两种求解权重系数 w 的办法：坐标降落法2（coordinate descent）、最小角回归法3（Least Angle Regression，LARS）
三、算法步骤
坐标降落法：
坐标降落法的外围与它的名称一样，就是沿着某一个坐标轴方向，通过一次一次的迭代更新权重系数的值，来慢慢迫近最优解。
具体步骤：
（1）初始化权重系数 w，例如初始化为零向量。
（2）遍历所有权重系数，顺次将其中一个权重系数当作变量，其余权重系数固定为上一次计算的后果当作常量，求出以后条件下只有一个权重系数变量的状况下的最优解。
在第 k 次迭代时，更新权重系数的办法如下：
$$\begin{matrix}w_m^k 示意第k次迭代，第m个权重系数 \\w_1^k = \underset{w_1}{\operatorname{argmin}} \left( \operatorname{Cost}(w_1, w_2^{k-1}, \dots, w_{m-1}^{k-1}, w_m^{k-1}) \right) \\ w_2^k = \underset{w_2}{\operatorname{argmin}} \left( \operatorname{Cost}(w_1^{k}, w_2, \dots, w_{m-1}^{k-1}, w_m^{k-1}) \right) \\\vdots \\w_m^k = \underset{w_m}{\operatorname{argmin}} \left( \operatorname{Cost}(w_1^{k}, w_2^{k}, \dots, w_{m-1}^{k}, w_m) \right) \\\end{matrix}$$
（3）步骤（2）为一次残缺迭代，当所有权重系数的变动不大或者达到最大迭代次数时，完结迭代。
By Nicoguaro - Own work
如上图所示，每次迭代固定其余的权重系数，只朝着其中一个坐标轴的方向更新，最初达到最优解。
最小角回归法：
最小角回归法是一种特征选择的办法，计算出每个特色的相关性，通过数学公式逐步计算出最优解。
具体步骤：
（1）初始化权重系数 w，例如初始化为零向量。
（2）初始化残差向量 residual 为指标标签向量 y - Xw，因为此时 w 为零向量，所以残差向量等于指标标签向量 y
（3）抉择一个与残差向量相关性最大的特征向量 x_i，沿着特征向量 xi 的方向找到一组权重系数 w，呈现另一个与残差向量相关性最大的特征向量 x_j 使得新的残差向量 residual 与这两个特征向量的相关性相等（即残差向量等于这两个特征向量的角平分向量上），从新计算残差向量。
（4）反复步骤（3），持续找到一组权重系数 w，使得第三个与残差向量相关性最大的特征向量 x_k 使得新的残差向量 residual 与这三个特征向量的相关性相等（即残差向量等于这三个特征向量的等角向量上），以此类推。
（5）当残差向量 residual 足够小或者所有特征向量都已被抉择，完结迭代。
Least Angle Regression - Figure 2
上图演示了只有两个特征向量时的状况，初始残差向量为 y_2，其中 x_1 向量与残差向量的相关性最大（向量夹角最小），找到一个 _11 使得新的残差向量 y_2 - x_1 _11 是 x_1 和 x_2 的角平分线（图中为 u_2），此时 w1 = _11, w2 = 0。最初找到一组 _21，_22 使得新的残差向量 y_2 - x_1 _11 - (x_1 _21 + x_2 _22) 尽可能小， 此时 w1 = _11 + _21，w2 = _22。所有特征向量都已被抉择，所以完结迭代。
坐标降落法绝对最小角回归法绝对好了解一些，每次只须要关怀一个权重系数即可。最小角回归法令是通过一些奇妙的数学公式缩小了迭代次数，使其的最坏计算复杂度和最小二乘法相似。
四、原理证实
Lasso回归代价函数为凸函数
同样须要证实：
$$f\left(\frac{x_{1}+x_{2}}{2}\right) \leq \frac{f\left(x_{1}\right)+f\left(x_{2}\right)}{2}$$
不等式右边：
$$\operatorname{Cost}\left(\frac{w_{1}+w_{2}}{2}\right)=\sum_{i=1}^{N}\left[\left(\frac{w_{1}+w_{2}}{2}\right)^{T} x_{i}-y_{i}\right]^{2}+\lambda\left\|\frac{w_{1}+w_{2}}{2}\right\|_{1}$$
不等式左边：
$$\frac{\operatorname{Cost}\left(w_{1}\right)+\operatorname{Cost}\left(w_{2}\right)}{2}=\frac{\sum_{i=1}^{N}\left(w_{1}^{T} x_{i}-y_{i}\right)^{2}+\sum_{i=1}^{N}\left(w_{2}^{T} x_{i}-y_{i}\right)^{2}+\lambda\left\|w_{1}\right\|_{1}+\lambda\left\|w_{2}\right\|_{1}}{2}$$
（1）不等式两边的前半部分与规范线性回归统一，只须要证实剩下的差值大于等于零即可
（2）依据向量范数的正齐次性，常数项系数能够乘进去
（3）将 提到括号外
（4）依据向量范数的定义，满足三角不等式，必然是大于等于零
$$\begin{aligned}\Delta &=\lambda\left\|w_{1}\right\|_{1}+\lambda\left\|w_{2}\right\|_{1}-2 \lambda\left\|\frac{w_{1}+w_{2}}{2}\right\|_{1} & (1) \\&=\lambda\left\|w_{1}\right\|_{1}+\lambda\left\|w_{2}\right\|_{1}-\lambda\left\|w_{1}+w_{2}\right\|_{1} & (2) \\&=\lambda\left(\left\|w_{1}\right\|_{1}+\left\|w_{2}\right\|_{1}-\left\|w_{1}+w_{2}\right\|_{1}\right) & (3) \\& \geq 0 & (4)\end{aligned}$$
人为的管制 的大小，最初的后果在实数范畴内必然大于等于零，证毕。
最小回归法数学公式
求单位角平分向量：
（1）设单位角平分向量为 u_A，能够将其看成选中的特征向量的线性组合，下标 A 示意选中的特色的序号汇合
（2）因为每个选中的特征向量与角平分向量的夹角都雷同，所以特征向量与角平分向量的点积后的向量每一个元素必然雷同，1_A 为元素都为 1 的向量，z 为常数
（3）依据（2）式能够示意出线性组合的系数向量 _A
$$\begin{matrix} u_\mathcal{A} = X_\mathcal{A} \omega_\mathcal{A} & (1)\\ X_\mathcal{A}^Tu_\mathcal{A} = X_\mathcal{A}^TX_\mathcal{A} \omega_\mathcal{A} = z 1_\mathcal{A} & (2)\\\omega_\mathcal{A} = z(X_\mathcal{A}^TX_\mathcal{A})^{-1} 1_\mathcal{A} & (3)\end{matrix}$$
（1）角平分向量 u_A 为单位向量，所以长度为 1
（2）改写成向量的模式
（3）带入上一步的线性组合的系数向量 _A
（4）依据转置的性质第一项的括号能够化简，提出常数 z
（5）矩阵乘以其逆矩阵为单位矩阵，所以能够化简
（6）最初解得常数 z 的表达式
$$\begin{array}{cc}\left\|u_{\mathcal{A}}\right\|^{2}=1 & (1)\\\omega_{\mathcal{A}}^{T} X_{\mathcal{A}}^{T} X_{\mathcal{A}} \omega_{\mathcal{A}}=1 & (2) \\\left(z\left(X_{\mathcal{A}}^{T} X_{\mathcal{A}}\right)^{-1} 1_{\mathcal{A}}\right)^{T} X_{\mathcal{A}}^{T} X_{\mathcal{A}} z\left(X_{\mathcal{A}}^{T} X_{\mathcal{A}}\right)^{-1} 1_{\mathcal{A}}=1 & (3) \\z^{2} 1_{\mathcal{A}}^{T}\left(X_{\mathcal{A}}^{T} X_{\mathcal{A}}\right)^{-1}\left(X_{\mathcal{A}}^{T} X_{\mathcal{A}}\right)\left(X_{\mathcal{A}}^{T} X_{\mathcal{A}}\right)^{-1} 1_{\mathcal{A}}=1 & (4) \\z^21_\mathcal{A}^T\left(X_\mathcal{A}^TX_\mathcal{A}\right)^{-1}1_\mathcal{A} = 1 & (5) \\z= \frac{1}{\sqrt[]{1_\mathcal{A}^T\left(X_\mathcal{A}^TX_\mathcal{A}\right)^{-1}1_\mathcal{A}}} = \left(1_\mathcal{A}^T\left(X_\mathcal{A}^TX_\mathcal{A}\right)^{-1}1_\mathcal{A}\right)^{-\frac{1}{2} } & (6) \\\end{array}$$
（1）将特征向量的转置乘以特征向量用 G_A 示意
（2）带入 G_A，失去 z 的表达式
（3）带入 G_A ，失去 _A 的表达式
$$\begin{array}{c}G_{\mathcal{A}}=X_{\mathcal{A}}^{T} X_{\mathcal{A}} & (1)\\z=\left(1_{\mathcal{A}}^{T}\left(G_{\mathcal{A}}\right)^{-1} 1_{\mathcal{A}}\right)^{-\frac{1}{2}} & (2)\\\omega_{\mathcal{A}}=z\left(G_{\mathcal{A}}\right)^{-1} 1_{\mathcal{A}} & (3)\end{array}$$
将X_A 乘以 _A ，就失去了角平分向量 u_A 的表达式，这样就能够求出单位角平分向量。更加具体的证实请参考 Bradley Efron 的论文4。
求角平分向量的长度：
（1）_A 示意以后预测值，沿着角平分向量的方向更新一个 长度
（2）C 示意特征向量与残差向量的相关性，为两个向量的点积，带入（1）式，失去相关性的更新表达式
（3）当特征向量为被选中的特征向量时，因为每个特征向量与角平分向量的乘积都雷同，都等于 z
（4）计算每个特征向量与角平分向量的乘积
（5）当特征向量不是被选中的特征向量时，应用 a 来带入（2）式
（6）要想退出下一个特征向量，则（3）式与（5）式的绝对值必然相等，能力保障下一次更新后的相关性也是雷同的
（7）解得 的表达式
$$\begin{array}{c}\mu_{A^{+}}=\mu_{A}+\gamma u_{A} & (1)\\C_{+}=X^{T}\left(y-\mu_{A^{+}}\right)=X^{T}\left(y-\mu_{A}-\gamma u_{A}\right)=C-\gamma X^{T} u_{A} & (2)\\C_{+}=C-\gamma z \quad(j \in A) & (3)\\a=X^{T} u_{A} & (4)\\C_{+}=C_{j}-\gamma a_{j} \quad(j \notin A) & (5)\\|C-\gamma z|=\left|C-\gamma a_{j}\right| & (6)\\\gamma=\min _{j \notin A}^{+}\left\{\frac{C-C_{j}}{z-a_{j}}, \frac{C+C_{j}}{z+a_{j}}\right\} & (7)\end{array}$$
这样就求出了 的表达式，也就是角平分向量的长度。更加具体的证实请参考 Bradley Efron 的论文4。
五、代码实现
应用 Python 实现Lasso回归算法（坐标降落法）：
def lassoUseCd(X, y, lambdas=0.1, max_iter=1000, tol=1e-4): """ Lasso回归，应用坐标降落法（coordinate descent） args: X - 训练数据集 y - 指标标签值 lambdas - 惩办项系数 max_iter - 最大迭代次数 tol - 变动量容忍值 return: w - 权重系数 """ # 初始化 w 为零向量 w = np.zeros(X.shape[1]) for it in range(max_iter): done = True # 遍历所有自变量 for i in range(0, len(w)): # 记录上一轮系数 weight = w[i] # 求出以后条件下的最佳系数 w[i] = down(X, y, w, i, lambdas) # 当其中一个系数变动量未达到其容忍值，持续循环 if (np.abs(weight - w[i]) > tol): done = False # 所有系数都变动不大时，完结循环 if (done): break return wdef down(X, y, w, index, lambdas=0.1): """ cost(w) = (x1 * w1 + x2 * w2 + ... - y)^2 + ... + (|w1| + |w2| + ...) 假如 w1 是变量，这时其余的值均为常数，带入上式后，其代价函数是对于 w1 的一元二次函数，能够写成下式： cost(w1) = (a * w1 + b)^2 + ... + |w1| + c (a,b,c, 均为常数) => 开展后 cost(w1) = aa * w1^2 + 2ab * w1 + |w1| + c (aa,ab,c, 均为常数) """ # 开展后的二次项的系数之和 aa = 0 # 开展后的一次项的系数之和 ab = 0 for i in range(X.shape[0]): # 括号内一次项的系数 a = X[i][index] # 括号内常数项的系数 b = X[i][:].dot(w) - a * w[index] - y[i] # 能够很容易的失去开展后的二次项的系数为括号内一次项的系数平方的和 aa = aa + a * a # 能够很容易的失去开展后的一次项的系数为括号内一次项的系数乘以括号内常数项的和 ab = ab + a * b # 因为是一元二次函数，当导数为零时，函数值最小值，只须要关注二次项系数、一次项系数和 return det(aa, ab, lambdas)def det(aa, ab, lambdas=0.1): """ 通过代价函数的导数求 w，当 w = 0 时，不可导 det(w) = 2aa * w + 2ab + = 0 (w > 0) => w = - (2 * ab + ) / (2 * aa) det(w) = 2aa * w + 2ab - = 0 (w < 0) => w = - (2 * ab - ) / (2 * aa) det(w) = NaN (w = 0) => w = 0 """ w = - (2 * ab + lambdas) / (2 * aa) if w < 0: w = - (2 * ab - lambdas) / (2 * aa) if w > 0: w = 0 return w应用 Python 实现Lasso回归算法（最小角回归法）：
def lassoUseLars(X, y, lambdas=0.1, max_iter=1000): """ Lasso回归，应用最小角回归法（Least Angle Regression） 论文：https://web.stanford.edu/~hastie/Papers/LARS/LeastAngle_2002.pdf args: X - 训练数据集 y - 指标标签值 lambdas - 惩办项系数 max_iter - 最大迭代次数 return: w - 权重系数 """ n, m = X.shape # 已被抉择的特色下标 active_set = set() # 以后预测向量 cur_pred = np.zeros((n,), dtype=np.float32) # 残差向量 residual = y - cur_pred # 特征向量与残差向量的点积，即相关性 cur_corr = X.T.dot(residual) # 选取相关性最大的下标 j = np.argmax(np.abs(cur_corr), 0) # 将下标增加至已被抉择的特色下标汇合 active_set.add(j) # 初始化权重系数 w = np.zeros((m,), dtype=np.float32) # 记录上一次的权重系数 prev_w = np.zeros((m,), dtype=np.float32) # 记录特色更新方向 sign = np.zeros((m,), dtype=np.int32) sign[j] = 1 # 均匀相关性 lambda_hat = None # 记录上一次均匀相关性 prev_lambda_hat = None for it in range(max_iter): # 计算残差向量 residual = y - cur_pred # 特征向量与残差向量的点积 cur_corr = X.T.dot(residual) # 以后相关性最大值 largest_abs_correlation = np.abs(cur_corr).max() # 计算以后均匀相关性 lambda_hat = largest_abs_correlation / n # 当均匀相关性小于时，提前结束迭代 # https://github.com/scikit-learn/scikit-learn/blob/2beed55847ee70d363bdbfe14ee4401438fba057/sklearn/linear_model/_least_angle.py#L542 if lambda_hat <= lambdas: if (it > 0 and lambda_hat != lambdas): ss = ((prev_lambda_hat - lambdas) / (prev_lambda_hat - lambda_hat)) # 从新计算权重系数 w[:] = prev_w + ss * (w - prev_w) break # 更新上一次均匀相关性 prev_lambda_hat = lambda_hat # 当全副特色都被抉择，完结迭代 if len(active_set) > m: break # 选中的特征向量 X_a = X[:, list(active_set)] # 论文中 X_a 的计算公式 - (2.4) X_a *= sign[list(active_set)] # 论文中 G_a 的计算公式 - (2.5) G_a = X_a.T.dot(X_a) G_a_inv = np.linalg.inv(G_a) G_a_inv_red_cols = np.sum(G_a_inv, 1) # 论文中 A_a 的计算公式 - (2.5) A_a = 1 / np.sqrt(np.sum(G_a_inv_red_cols)) # 论文中 的计算公式 - (2.6) omega = A_a * G_a_inv_red_cols # 论文中角平分向量的计算公式 - (2.6) equiangular = X_a.dot(omega) # 论文中 a 的计算公式 - (2.11) cos_angle = X.T.dot(equiangular) # 论文中的 gamma = None # 下一个抉择的特色下标 next_j = None # 下一个特色的方向 next_sign = 0 for j in range(m): if j in active_set: continue # 论文中 的计算方法 - (2.13) v0 = (largest_abs_correlation - cur_corr[j]) / (A_a - cos_angle[j]).item() v1 = (largest_abs_correlation + cur_corr[j]) / (A_a + cos_angle[j]).item() if v0 > 0 and (gamma is None or v0 < gamma): gamma = v0 next_j = j next_sign = 1 if v1 > 0 and (gamma is None or v1 < gamma): gamma = v1 next_j = j next_sign = -1 if gamma is None: # 论文中 的计算方法 - (2.21) gamma = largest_abs_correlation / A_a # 选中的特征向量 sa = X_a # 角平分向量 sb = equiangular * gamma # 解线性方程（sa * sx = sb） sx = np.linalg.lstsq(sa, sb) # 记录上一次的权重系数 prev_w = w.copy() d_hat = np.zeros((m,), dtype=np.int32) for i, j in enumerate(active_set): # 更新以后的权重系数 w[j] += sx[0][i] * sign[j] # 论文中 d_hat 的计算方法 - (3.3) d_hat[j] = omega[i] * sign[j] # 论文中 _j 的计算方法 - (3.4) gamma_hat = -w / d_hat # 论文中 _hat 的计算方法 - (3.5) gamma_hat_min = float("+inf") # 论文中 _hat 的下标 gamma_hat_min_idx = None for i, j in enumerate(gamma_hat): if j <= 0: continue if gamma_hat_min > j: gamma_hat_min = j gamma_hat_min_idx = i if gamma_hat_min < gamma: # 更新以后预测向量 - (3.6) cur_pred = cur_pred + gamma_hat_min * equiangular # 将下标移除至已被抉择的特色下标汇合 active_set.remove(gamma_hat_min_idx) # 更新特色更新方向汇合 sign[gamma_hat_min_idx] = 0 else: # 更新以后预测向量 cur_pred = X.dot(w) # 将下标增加至已被抉择的特色下标汇合 active_set.add(next_j) # 更新特色更新方向汇合 sign[next_j] = next_sign return w六、第三方库实现
scikit-learn6 实现（坐标降落法）：
from sklearn.linear_model import Lasso# 初始化Lasso回归器，默认应用坐标降落法reg = Lasso(alpha=0.1, fit_intercept=False)# 拟合线性模型reg.fit(X, y)# 权重系数w = reg.coef_scikit-learn7 实现（最小角回归法）：
from sklearn.linear_model import LassoLars# 初始化Lasso回归器，应用最小角回归法reg = LassoLars(alpha=0.1, fit_intercept=False)# 拟合线性模型reg.fit(X, y)# 权重系数w = reg.coef_七、动画演示
上面动图展现了后面一节的工作年限与均匀月工资的例子，每次只朝这坐标轴的一个方向扭转权重系数，逐步迫近最优解的过程。
<center>坐标降落法</center>
上面动图展现了 对各个自变量权重系数的影响，横轴为惩办系数 ，纵轴为权重系数，每一个色彩示意一个自变量的权重系数（训练数据来源于sklearn diabetes datasets）：
<center> 对权重系数的影响</center>
能够看到当 逐步增大时（ 向左挪动），某些特色的权重系数会疾速变成零，通过这个性质阐明Lasso回归能够用来做特征选择，管制 的大小来抉择出要害特色。
八、思维导图
九、参考文献
- https://en.wikipedia.org/wiki...(statistics)
- https://en.wikipedia.org/wiki...
- https://en.wikipedia.org/wiki...
- https://web.stanford.edu/~has...
- https://zhuanlan.zhihu.com/p/...
- https://scikit-learn.org/stab...
- https://scikit-learn.org/stab...
残缺演示请点击这里
注：本文力求精确并通俗易懂，但因为笔者也是初学者，程度无限，如文中存在谬误或脱漏之处，恳请读者通过留言的形式批评指正
本文首发于——AI导图，欢送关注