Ridge Regression and LASSO Regression.

岭回归与LASSO回归是将正则化思想引入线性回归模型后得到的方法。

岭回归对模型的参数加入L2正则化，能够有效的防止模型过拟合，解决非满秩下求逆困难的问题。
LASSO回归全称为Least Absolute Shrinkage and Selection Operator(最小绝对收缩选择算子, 套索算法)，对模型的参数加入L1正则化，能够稀疏矩阵，进行庞大特征数量下的特征选择。

1. Ridge Regression

岭回归（Ridge Regression）是引入了L2正则化的线性回归，求解问题如下：

\[L(w) = \sum_{i=1}^{N} {(w^Tx^{(i)}-y^{(i)})^2} + λ||w||_2 = (Xw-y)^T(Xw-y)+ λw^Tw\]

岭回归通过增加对权重的限制，使得模型把有限的权重放到更重要的特征维度上；并且每个权重都不会太大，否则自变量的微小变化将会引起输入的巨大变化。

令梯度为零可以得到：

\[\nabla_wL(w) = 2X^TXw - 2Xy + 2λw = 0\]

该问题的解析解为：

\[w = (X^TX+λI)^{-1}Xy\]

注意到正则化系数$λ$通常大于零，故矩阵$X^TX+λI$一定可逆。

正则化系数$λ$越大，不重要的维度的权重会减小；因此通过调整正则化系数可以实现初步的特征选择。下图给出了某岭回归模型中八个维度的权重随正则化系数的变化曲线。从图中可以看出，特征维度4和5的重要性较大，并且特征4起到正作用，特征5起到反作用。

求解岭回归问题的程序如下：

def RidgeRegression(X, Y):
    return np.dot(np.linalg.inv(np.dot(X.T,X)+lambd*np.eye(X.shape[1])), np.dot(X,Y))

⚪ 讨论：岭回归等价于噪声和先验服从高斯分布的最大后验估计

引入高斯噪声$ε$~$N(0,σ^2)$，对线性回归建模：

\[y = w^Tx + ε\]

贝叶斯角度认为参数$w$不再是常数，而是随机变量，假设其先验概率为$N(0,σ_0^2)$,

由贝叶斯定理，参数$w$的后验概率：

\[P(w | y) = \frac{P(y | w)P(w)}{P(y)} \propto P(y | w)P(w)\]

由最大后验估计：

\[\begin{aligned} \hat{w} &= \mathop{\arg \max}_{w}\log P(w | y) = \mathop{\arg \max}_{w}\log P(y | w)P(w) \\ &= \mathop{\arg \max}_{w} \log (\prod_{i=1}^{N} {\frac{1}{\sqrt{2\pi}σ}\exp(-\frac{(y^{(i)}-w^Tx^{(i)})^2}{2σ^2})\frac{1}{\sqrt{2\pi}σ_0}\exp(-\frac{w^Tw}{2σ_0^2})}) \\ &= \mathop{\arg \max}_{w} \sum_{i=1}^{N} {\log (\frac{1}{\sqrt{2\pi}σ}\exp(-\frac{(y^{(i)}-w^Tx^{(i)})^2}{2σ^2})\frac{1}{\sqrt{2\pi}σ_0}\exp(-\frac{w^Tw}{2σ_0^2}))} \\ &\propto \mathop{\arg \max}_{w} \sum_{i=1}^{N} {-\frac{(y^{(i)}-w^Tx^{(i)})^2}{2σ^2}-\frac{w^Tw}{2σ_0^2}} \\ &= \mathop{\arg \min}_{w} \sum_{i=1}^{N} {(y^{(i)}-w^Tx^{(i)})^2+\frac{σ^2}{σ_0^2}w^Tw} \end{aligned}\]

该问题等价于引入L2正则化的最小二乘法（岭回归）。

2. Kernel Ridge Regression

将核方法引入岭回归，即可得到核岭回归(Kernelized Ridge Regression)。

岭回归的目标函数：

\[L(w) = \frac{1}{N}\sum_{i=1}^{N} {(w^Tx^{(i)}-y^{(i)})^2} + \frac{λ}{N}w^Tw\]

根据表示定理，上述优化问题的最优解$w$可以表示为所有样本的线性组合：

\[w^* = \sum_{n=1}^{N} {β_nx^{(n)}}\]

代入目标函数：

\[L(β) = \frac{1}{N} \sum_{i=1}^{N} {(\sum_{n=1}^{N} {β_n{x^{(n)}}^Tx^{(i)}}-y^{(i)})^2} + \frac{λ}{N}\sum_{i=1}^{N} {\sum_{j=1}^{N} {β_iβ_j{x^{(i)}}^Tx^{(j)}}}\]

引入核函数来代替样本的特征转换和内积运算，记$K(x,x')={φ(x)}^Tφ(x')$，则核岭回归的损失函数为：

\[L(β) = \frac{1}{N} \sum_{i=1}^{N} {(\sum_{n=1}^{N} {β_nK(x^{(n)},x^{(i)})}-y^{(i)})^2} + \frac{λ}{N}\sum_{i=1}^{N} {\sum_{j=1}^{N} {β_iβ_jK(x^{(i)},x^{(j)})}}\]

上式写作矩阵形式：

\[L(β) = \frac{1}{N} (y-Kβ)^T(y-Kβ) + \frac{λ}{N}β^TKβ\]

对$β$求梯度，令其为零，(注意到$K$是对称矩阵)可得：

\[\nabla_βL(β) = \frac{1}{N}(2K^TKβ-2K^Ty) + \frac{λ}{N}2Kβ = \frac{2K^T}{N}(Kβ-y+λβ) = 0\]

可解得引入核岭回归的解析解：

\[β = (K+λI)^{-1}y\]

对矩阵$K+λI$的讨论：

由于$λ>0$，$K$是半正定矩阵，其逆矩阵一定存在。
该矩阵是一个稠密(dense)矩阵，大部分元素不是$0$，求逆过程计算量较大。

求得$β$后，可以得到回归函数：

\[y = \sum_{n=1}^{N} {β_nφ(x^{(n)})^Tφ(x)} = \sum_{n=1}^{N} {β_nK(x^{(n)},x)}\]

3. LASSO Regression

lasso (英/læˈsuː/ 美/ˈlæsoʊ/): n. (捕马、套牛等用的)套索

LASSO回归是引入了L1正则化的线性回归，求解问题如下：

\[L(w) = \sum_{i=1}^{N} {(w^Tx^{(i)}-y^{(i)})^2} + λ\sum_{d=1}^{D} {| w_d |}= (Xw-y)^T(Xw-y)+ λ||w||_1\]

选择L1范数使得学习得到的权重更加稀疏，不重要的权重接近$0$。

由于引入了L1正则化项，因此LASSO回归的目标函数不满足处处可到，没有闭式解；在实践中常采用梯度下降搜索法(坐标轴下降法)求解，即依次对权重$w_d$求局部最优，使其梯度趋近于0，不断迭代直至所有权重都不产生显著变化。

LASSO回归的目标函数可写作：

\[\begin{aligned} L(w) &= \sum_{i=1}^{N} {(w^Tx^{(i)}-y^{(i)})^2} + λ\sum_{d=1}^{D} {| w_d |} \\ &= \sum_{i=1}^{N} {(\sum_{d=1}^Dw_dx_d^{(i)}-y^{(i)})^2} + λ\sum_{d=1}^{D} {| w_d |} \\ &= \sum_{i=1}^{N} {((\sum_{d=1}^Dw_dx_d^{(i)})^2+(y^{(i)})^2-2(\sum_{d=1}^Dw_dx_d^{(i)})y^{(i)})} + λ\sum_{d=1}^{D} {| w_d |}\end{aligned}\]

求上式对权重$w_d$的梯度：

\[\begin{aligned} \frac{\partial L(w)}{w_d} &= \sum_{i=1}^{N} (2(\sum_{d=1}^Dw_dx_d^{(i)})x_d^{(i)}-2x_d^{(i)}y^{(i)}) + \lambda \text{sign}(w_d) \\ &= 2 \sum_{i=1}^{N} ((\sum_{d=1}^Dw_dx_d^{(i)}-y^{(i)})x_d^{(i)}) + \lambda \text{sign}(w_d) \\ &= 2 x_d^T(w^Tx^{(i)}-y^{(i)}) + \lambda \text{sign}(w_d) \end{aligned}\]

则可以通过梯度下降更新权重$w_d$：

\[w_d \leftarrow w_d - \alpha \cdot (2 x_d^T(w^Tx^{(i)}-y^{(i)}) + \lambda \text{sign}(w_d))\]

# 坐标轴下降法
def CoordinateDescent(X, Y, epochs, lr, lam):
    N, D= X.shape
    w = np.ones([D, 1])
    # 进行 epoches 轮迭代
    for k in range(epochs):
        # 保存上一轮的w
        pre_w = copy.copy(w)
        # 逐维度进行参数寻优
        for d in range(D):
            # 在每个维度上找到最优的w_d
            for j in range(epochs):
                Y_hat = X*w
                g_d = 2*X[:,d].T*(Y_hat-Y) + lam*np.sign(w[d])
                # 进行梯度下降
                w[d] = w[d] - g_d*lr
                if np.abs(g_d) < 1e-3:
                    break
        # 计算上一轮的w和当前轮w的差值，如果每个维度的w都没有什么变化则退出
        diff_w  = np.array(list(map(lambda x:abs(x)<1e-3,pre_w-w)))
        if diff_w.all():
            break
    return w

⚪ 讨论：LASSO回归等价于噪声服从高斯分布、参数服从拉普拉斯分布的最大后验估计。

引入高斯噪声$ε$~$N(0,σ^2)$，对线性回归建模：

\[y = w^Tx + ε\]

贝叶斯角度认为参数$w$不再是常数，而是随机变量，假设其先验概率为拉普拉斯分布$L(0,σ_0^2)$:

\[w \text{~} L(0,σ_0^2) = \frac{1}{2σ_0^2} \exp(-\frac{|w|}{σ_0^2})\]

由贝叶斯定理，参数$w$的后验概率：

\[P(w | y) = \frac{P(y | w)P(w)}{P(y)} \propto P(y | w)P(w)\]

由最大后验估计：

\[\begin{aligned} \hat{w} &= \mathop{\arg \max}_{w}\log P(w | y) = \mathop{\arg \max}_{w}\log P(y | w)P(w) \\ &= \mathop{\arg \max}_{w} \log (\prod_{i=1}^{N} {\frac{1}{\sqrt{2\pi}σ}\exp(-\frac{(y^{(i)}-w^Tx^{(i)})^2}{2σ^2}) \frac{1}{2σ_0^2} \exp(-\frac{|w|}{σ_0^2}) }) \\ &= \mathop{\arg \max}_{w} \sum_{i=1}^{N} {\log (\frac{1}{\sqrt{2\pi}σ}\exp(-\frac{(y^{(i)}-w^Tx^{(i)})^2}{2σ^2})\frac{1}{2σ_0^2} \exp(-\frac{|w|}{σ_0^2}))} \\ &\propto \mathop{\arg \max}_{w} \sum_{i=1}^{N} {-\frac{(y^{(i)}-w^Tx^{(i)})^2}{2σ^2} -\frac{|w|}{σ_0^2}} \\ &= \mathop{\arg \min}_{w} \sum_{i=1}^{N} {(y^{(i)}-w^Tx^{(i)})^2+\frac{2σ^2}{σ_0^2}|w|} \end{aligned}\]

该问题等价于引入L1正则化的最小二乘法（LASSO回归）。