机器学习基础 线性回归与 Softmax 回归

机器学习基础 线性回归与 Softmax 回归

1. 最小二乘法

1.1 数据集定义

​ 给定数据集 $D=\{(x_1,y_1),\cdots ,(x_N,y_N)\}$，其中 $x_i\in {\mathbb{R}}^p$，$y_i\in \mathbb{R}$，$i=1,2,\cdots ,N$。对于输入样本定义为：
$$\begin{array}{r}X=(x_1{x}_2\cdots x_N{)}^T=\left[\begin{array}{c}{x}_1^T\\ x_2^T\\ ⋮\\ x_N^T\end{array}\right]={\left[\begin{array}{cccc}{x}_{11}& x_{12}& \cdots & x_{1p}\\ x_{21}& x_{22}& \cdots & x_{2p}\\ ⋮& ⋮& \ddots & ⋮\\ x_{N1}& x_{N2}& \cdots & x_{Np}\end{array}\right]}_{N\times p}.\end{array}\text{\hspace{1em}(1)}$$
​ 定义输出为：
$$\begin{array}{r}Y={\left[\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_N\end{array}\right]}_{N\times 1}.\end{array}\text{\hspace{1em}(2)}$$

1.2 最小二乘的矩阵推导

​ 最小二乘法的损失函数定义为：
$$\begin{array}{rl}L(w)& =\sum _{i=1}^N\Vert w^T{x}_i-y_i{\Vert }^2\\ & =\sum _{i=1}^N(w^T{x}_i-y_i{)}^2\\ & =(w^T{x}_1-y_1{w}^T{x}_2-y_2\cdots w^T{x}_N-y_N)\left[\begin{array}{c}{w}^T{x}_1-y_1\\ w^T{x}_2-y_2\\ ⋮\\ w^T{x}_N-y_N\end{array}\right],\end{array}\text{\hspace{1em}(3)}$$
第一项可继续化简：
$$\begin{array}{rl}& =(w^T{x}_1{w}^T{x}_2\cdots w^T{x}_N)-(y_1{y}_2\cdots y_N)\\ & =w^T(x_1{x}_2\cdots x_N)-(y_1{y}_2\cdots y_N)\\ & =w^T{X}^T-Y^T.\end{array}\text{\hspace{1em}(4)}$$
​ 因为 $x_i^T=[x_{i1}{x}_{i2}\cdots x_{ip}]$ 是 $1\times p$ 的行向量，那么 $x_i$ 是一个 $p\times 1$ 的列向量。$w$ 作为权重系数，是一个 $p\times 1$ 的列向量，转置后就变成了 $1\times p$ 的行向量。故第二项可以写成如下的矩阵形式：
$$\begin{array}{r}Xw-Y\end{array}.\text{\hspace{1em}(5)}$$
​ 因此，可以将损失函数表示为：
$$\begin{array}{rl}L(w)& =(w^T{X}^T-Y^T)(Xw-Y)\\ & =w^T{X}^{T}Xw-Y^{T}Xw-w^T{X}^{T}Y+Y^{T}Y\\ & =w^T{X}^{T}Xw-2Y^{T}Xw+Y^{T}Y.\end{array}\text{\hspace{1em}(6)}$$
​ 最小二乘的目标是求解损失函数 $L(w)$ 的最小值，即：
$$\begin{array}{r}\hat{w}=\underset{w}{\mathrm{argmin}}L(w).\end{array}\text{\hspace{1em}(7)}$$
​ 对矩阵求导有：
$$\begin{array}{rl}\frac{\partial L(w)}{\partial w}& =2X^{T}Xw-2X^{T}Y\stackrel{△}{=}0\\ & \Rightarrow X^{T}Xw=X^{T}Y\\ & \Rightarrow w=(X^{T}X{)}^{-1}{X}^{T}Y,\end{array}\text{\hspace{1em}(8)}$$
其中，$X^{T}X=X^{+}$ 表示伪逆。

1.3 最小二乘的几何解释

​ 公式 (3) 可以理解为线性回归总的误差，假设拟合的是二维平面上的一条直线，那么第一种几何解释就是每个样本点的拟合值到真实值的距离，可视化如图 1 所示。

图1 基于距离的几何解释.

​ 我们假设拟合的直线是 $f(w)=Xw$，如果横着看 $X$，每一行可以理解为是一个样本。如果竖着看，那么每一列都是一个列向量，这些列向量都存在于 $N$ 维的空间中。

​ 当列向量线性无关时，$p$ 个列向量会张成一个 $p$ 维的子空间。比如，两个列向量在 3 维空间中张成了一个平面。在一般情况下，观测向量 $Y$ 不会完美地被直线拟合，故 $Y$ 不会存在于 $X$ 的列空间中。

​ 最小二乘法的另一种几何解释即为，在 $p$ 维子空间中寻找一个向量 $f(w)$，使得 $f(w)$ 与 $Y$ 的距离最小。显然是当 $f(w)$ 为 $Y$ 的投影时距离最小，则满足 $Y-f(w)$ 与 $p$ 维子空间的基向量正交，故求解：
$$\begin{array}{rl}& X^T(Y-f(w))=0\\ & \Rightarrow X^{T}Y=X^{T}Xw\\ & \Rightarrow w=(X^{T}X{)}^{-1}{X}^{T}Y.\end{array}\text{\hspace{1em}(9)}$$

图2 基于子空间的几何解释.

1.4 概率视角下的最小二乘估计

​ 假设 $y=w^{T}x+ϵ$，这里的 $x$ 是 $p\times 1$ 维的列向量，$ϵ$ 是数据中的噪声，且 $ϵ\sim \mathcal{N}(0,{\sigma }^2)$。因此，对于 $y$ 有：
$$\begin{array}{r}y|x;w\sim \mathcal{N}(w^{T}x,{\sigma }^2).\end{array}\text{\hspace{1em}(10)}$$
​ 因为 $y$ 也服从正态分布，接下来可以用极大似然估计求解 $w$。假设 $N$ 个样本之间独立同分布，有以下似然函数：
$$\begin{array}{r}\mathcal{L}(w)=\text{log}P(Y|X,w)=\text{log}\prod _{i=1}^{N}P(y_i|x_i,w)=\sum _{i=1}^N\text{log}P(y_i|x_i,w),\end{array}\text{\hspace{1em}(11)}$$
其中，正态分布的概率密度函数：
P ( y ∣ x ; w ) = 1 2 π σ exp { − ( y − w T x ) 2 2 σ 2 } , (12) \begin{align} P(y|x;w)=\frac{1}{\sqrt{2\pi}\sigma}\text{exp}\{-\frac{(y-w^Tx)^2}{2\sigma^2}\} \end{align},\tag{12} P(y∣x;w)=2π ​σ1​exp{−2σ2(y−wTx)2​}​​,(12)
代回似然函数中可得：
L ( w ) = ∑ i = 1 N log P ( y i ∣ x i , w ) = ∑ i = 1 N log 1 2 π σ exp { − ( y i − w T x i ) 2 2 σ 2 } = ∑ i = 1 N ( log 1 2 π σ − ( y i − w T x i ) 2 2 σ 2 ) . \begin{align} \mathcal{L}(w)&=\sum_{i=1}^N \text{log}\; P(y_i|x_i,w)\\[4pt] &=\sum_{i=1}^N \text{log}\; \frac{1}{\sqrt{2\pi}\sigma}\text{exp}\{-\frac{(y_i-w^Tx_i)^2}{2\sigma^2}\}\\[4pt] &=\sum_{i=1}^N (\text{log}\; \frac{1}{\sqrt{2\pi}\sigma} - \frac{(y_i-w^Tx_i)^2}{2\sigma^2}).\tag{13} \end{align} L(w)​=i=1∑N​logP(yi​∣xi​,w)=i=1∑N​log2π ​σ1​exp{−2σ2(yi​−wTxi​)2​}=i=1∑N​(log2π ​σ1​−2σ2(yi​−wTxi​)2​).​(13)​
​ 极大似然估计有：
$$\begin{array}{rl}\hat{w}& =\underset{w}{\mathrm{argmax}}\mathcal{L}(w)\\ & =\underset{w}{\mathrm{argmax}}-\frac{1}{2{\sigma }^2}(y_i-w^T{x}_i{)}^2\\ & =\underset{w}{\mathrm{argmin}}(y_i-w^T{x}_i{)}^2.\end{array}\text{\hspace{1em}(14)}$$
​ 由此可见，最小二乘估计等价于噪声服从正态分布的极大似然估计。

2. 正则化

2.1 L1 范数与 L2 范数

​ L1 正则化和 L2 正则化可以看做是损失函数的惩罚项。所谓惩罚是指对损失函数中的某些参数做一些限制。对于线性回归模型，使用 L1 正则化的模型建叫做 Lasso 回归，使用 L2 正则化的模型叫做 Ridge 回归。其中，Lasso 回归的形式如下：
$$\underset{w}{\mathrm{argmin}}\frac{1}{2n}\Vert Xw-y{\Vert }_2^2+\alpha \Vert w{\Vert }_1.\text{\hspace{1em}(15)}$$
​ Ridge 回归的形式如下：
$$\underset{w}{\mathrm{argmin}}\Vert Xw-y{\Vert }_2^2+\alpha \Vert w{\Vert }_2^2.\text{\hspace{1em}(16)}$$

2.2 正则化的作用

​ 如图 3 所示，其目标在于在满足正则化约束的前提下，让损失函数最小。也就是说，让误差函数的等值线（椭圆）尽可能小，但又不能越过正则化约束线（菱形或圆），交点即为最优解。

​ 由于菱形有“尖角”正好在坐标轴上，故 L1 范数得到的最优解，恰好可能落在某条坐标轴上。进而某些参数会被压缩为 0，相当于自动做了“特征选择”。

​ L2 范数的作用是改善过拟合。过拟合是：模型训练时候的误差很小，但是测试误差很大，也就是说模型复杂到可以拟合到所有训练数据，但在预测新的数据的时候，结果很差。相对于 L1 范数，L2 范数可以使得 $w$ 的每个元素都很小，都接近于 0。

图3 L1 正则化比 L2 正则化更容易得到稀疏解.

2.3 Lasso 回归的求解

2.3.1 L-Lipschitz 条件

​ 较为常见的是使用近端梯度下降法对 Lasso 进行求解，近端梯度下降要求如下形式的 $\nabla f$ 满足 L-Lipschitz 条件
$$\underset{x}{\min }f(x)+h(x)\text{\hspace{1em}(17)}$$
​ 因此，首先证明 Lasso 表达式（公式 15）的第一部分的梯度满足 L-Lipschitz 条件，即证明线性回归模型的均方误差满足 L-Lipschitz 条件。这里引用知乎用户 cielo 的回答，说得非常形象。

L-Lipschitz 连续，要求函数图像的曲线上任意两点连线的斜率一致有界，就是任意的斜率都小于同一个常数，这个常数就是 L-Lipschitz 常数。

从局部看：我们可以取两个充分接近的点，如果这个时候斜率的极限存在的话，这个斜率的极限就是这个点的导数。也就是说函数可导，又是 L-Lipschitz 连续，那么导数有界。反过来，如果可导函数，导数有界，可以推出函数 L-Lipschitz 连续。

从整体看：L-Lipschitz 连续要求函数在无限的区间上不能有超过线性的增长，所以 $x^2$, $e^x$ 这些函数在无限区间上不是 L-Lipschitz 连续的。

​ 这里可以通过维基百科的示意图进行辅助解释：

图4 维基百科上对 L-Lipschitz 条件的示意图

2.3.2 L-Lipschitz 条件的证明

​ 定义线性回归的均方误差函数为：
$$f(w)=\frac{1}{2n}\Vert y-Xw{\Vert }_2^2.\text{\hspace{1em}(18)}$$
​ Lasso 回归的 $f(x)$ 部分与线性回归的均方误差函数相同，下面证明该部分满足 L-Lipschitz 条件，定义：
$$F(x_1,x_2)=\Vert \nabla f(x_1)-\nabla f(x_2){\Vert }_2^2-L\Vert x_1-x_2{\Vert }_2^2.\text{\hspace{1em}(19)}$$
​ 首先求解 $\nabla f(w)$，对式 (18) 进行展开得：
$$f(w)=\frac{1}{2n}(y-Xw{)}^T(y-Xw).\text{\hspace{1em}(20)}$$
​ 根据矩阵求导法则：
$$\frac{\partial }{\partial w}(a^{T}a)=2a^T\frac{\partial a}{\partial w}.\text{\hspace{1em}(21)}$$
​ 令 $a=(y-Xw)$，对式 (18) 求导后转置，故 $\nabla f(w)$ 为：
$$\nabla f(w)=-\frac{1}{n}{X}^T(y-Xw),\text{\hspace{1em}(22)}$$
代入式 (19) 得：
$$\begin{array}{rl}F(w_1,w_2)& =\Vert \frac{1}{n}{X}^T(X{w}_1-y)+\frac{1}{n}{X}^T(X{w}_2+y){\Vert }_2^2-L\Vert w_1-w_2{\Vert }_2^2\\ & =\frac{1}{n^2}\Vert X^{T}X(w_1-w_2){\Vert }_2^2-L\Vert w_1-w_2{\Vert }_2^2\end{array}.\text{\hspace{1em}(23)}$$
​ 现在要证明 $F(w_1,x_2)\le 0$，则要证明下式：
$$\frac{\Vert X^{T}X(w_1-w_2){\Vert }_2^2}{n^2\Vert w_1-w_2{\Vert }_2^2}\le L,\text{\hspace{1em}(24)}$$
记 $A=X^{T}X$，实际上是证明对于任意的 $x$，下式有界：
$$\frac{\Vert Ax{\Vert }_2^2}{\Vert x{\Vert }_2^2}.\text{\hspace{1em}(25)}$$
​ 由于 $A$ 是实对称矩阵，因此可以被对角化，即存在一个正交矩阵 $Q$ 和一个对角矩阵 $\Lambda$，使得：
$$A=Q\Lambda Q^T,\text{\hspace{1em}(26)}$$
其中，$Q$ 是 $A$ 的特征向量，$\Lambda$ 的对角线元素是 $A$ 的特征值。

​ 对于特征向量的正交性，如果 $v_i$ 和 $v_j $ 是两个不同的特征向量，分别对应不同的特征值 ${\lambda }_i\ne {\lambda }_j$，满足：
$$A{v}_i={\lambda }_i{v}_i,A{v}_j={\lambda }_j{v}_j,\text{\hspace{1em}(27)}$$
则它们之间满足正交关系：
$$v_i^T{v}_j=0.\text{\hspace{1em}(28)}$$
​ 除此之外，可以进一步选择特征向量使其单位化，即 $\Vert v_i\Vert =1$，使得它们形成一个标准正交基。

​ 基于此，$x$ 可以表示为 $A$ 的特征向量的线性组合：
$$x=\sum _{i=1}^n{k}_i{v}_i,\text{\hspace{1em}(29)}$$
其中，$v_i$ 是 $A$ 的特征向量，$k_i$ 是系数。

​ 因此，$Ax$ 可以化简为：
$$Ax=A(\sum _{i=1}^n{k}_i{v}_i)=\sum _{i=1}^n{k}_{i}A{v}_i=\sum _{i=1}^n{k}_i{\lambda }_i{v}_i,\text{\hspace{1em}(30)}$$
其中，$\Vert Ax{\Vert }_2^2$ 可以化简为：
$$\Vert Ax{\Vert }_2^2=\Vert \sum _{i=1}^n{k}_i{\lambda }_i{v}_i{\Vert }_2^2=\sum _{i=1}^n(k_i{\lambda }_i{)}^2\Vert v_i{\Vert }_2^2=\sum _{i=1}^n(k_i{\lambda }_i{)}^2.\text{\hspace{1em}(31)}$$
​ 同理，$\Vert x{\Vert }_2^2$可以化简为：
$$\Vert x{\Vert }_2^2=\sum _{i=1}^n{k_i}^2.\text{\hspace{1em}(32)}$$
​ 综上所述，可以将式 (25) 化简为：
$$\frac{\Vert Ax{\Vert }_2^2}{\Vert x{\Vert }_2^2}=\frac{{\sum }_{i=1}^n(k_i{\lambda }_i{)}^2}{{\sum }_{i=1}^n{k_i}^2}.\text{\hspace{1em}(33)}$$
​ 式 (33) 是一个加权平均式，因此有如下结论：
$$\underset{i}{\min }({\lambda }_i^2)\le \frac{{\sum }_{i=1}^n(k_i{\lambda }_i{)}^2}{{\sum }_{i=1}^n{k_i}^2}\le \underset{i}{\max }({\lambda }_i^2).\text{\hspace{1em}(34)}$$
​ 因此式 (25) 有界，即：
$$\frac{\Vert X^{T}X(w_1-w_2){\Vert }_2^2}{\Vert w_1-w_2{\Vert }_2^2}\le \underset{i}{\max }({\lambda }_i^2).\text{\hspace{1em}(35)}$$
​ 接下来求解符合条件的 $L$：
$$\begin{array}{rl}F(w_1,w_2)& =\frac{1}{n^2}\Vert X^{T}X(w_1-w_2){\Vert }_2^2-L\Vert w_1-w_2{\Vert }_2^2\\ & =L\Vert w_1-w_2{\Vert }_2^2(\frac{\Vert X^{T}X(w_1-w_2){\Vert }_2^2}{n^{2}L\Vert w_1-w_2{\Vert }_2^2}-1)\\ & \le L\Vert w_1-w_2{\Vert }_2^2(\frac{1}{n^{2}L}\underset{i}{\max }({\lambda }_i^2)-1)\le 0.\end{array}\text{\hspace{1em}(36)}$$
​ 故可得 L-Lipschitz 常数的上界：
$$L\ge \frac{{\max }_i({\lambda }_i^2)}{n^2}.\text{\hspace{1em}(37)}$$

2.4 近端梯度下降

2.4.1 近端梯度下降的定义

​ 在上一节已经证明了 Lasso 的 $\nabla f$ 满足 L-Lipschitz 条件，因此可以使用近端梯度下降算法进行求解。实际上，近端梯度下降的本质可以视为先对光滑的函数 $f(x)$ 进行梯度下降，再对正则项进行隐性的梯度下降。

​ 梯度的 L-Lipschitz 连续性要求梯度的变化不能过于剧烈，保证了存在合适的步长，使得 $f(x)$ 沿着梯度方向前进这个步长可以使函数值下降，即 descent lemma，这可用于证明算法收敛性。

​ 近端梯度算法（Proximal Gradient Descent, PGD）常用来解决如下形式的优化问题：
$$\underset{x}{\min }f(x)+h(x).\text{\hspace{1em}(38)}$$
​ 近端梯度下降的 $h(x)$ 有以下几种形式：

1. $h(x)=0$：则为一般的梯度下降算法。

2. $h(x)=I_C(x)$：则称为投影梯度下降（Projected gradient descent），其中：
   $$I_C(x)=\{\begin{array}{ll}0,& x\in C,\\ \infty ,& x\notin C.\end{array}\text{\hspace{1em}(39)}$$

3. $h(x)=\lambda \Vert x{\Vert }_1$：则被称为迭代收缩阈值算法（Iterative Shrinkage Thresholding Algorithm, ISTA），Lasso 回归为该种类型。

2.4.2 迭代方法

​ 对于函数 $f(x)$，在 $x_k$ 附近进行二阶泰勒展开得：
$$\hat{f}(x)\simeq f(x_k)+\nabla f(x_k{)}^T(x-x_k)+\frac{1}{2}(x-x_k{)}^T{\nabla }^{2}f(x_k)(x-x_k),\text{\hspace{1em}(40)}$$
上式中，${\nabla }^{2}f(x_k)$ 为 Hesse 矩阵，并且正定，L-Lipschitz 常数是 Hesse 矩阵的上界，在优化算法中常这样进行简化，故上式可以化简为如下形式：
$$\hat{f}(x)\le f(x_k)+⟨\nabla f(x_k),(x-x_k)⟩+\frac{L}{2}\Vert (x-x_k){\Vert }_2^2,\text{\hspace{1em}(41)}$$
上述化简中，向量内积写法为：$\nabla f(x_k{)}^T(x-x_k)=⟨\nabla f(x_k),x-x_k⟩$，范数平方写法为：$\Vert x-x_k{\Vert }_2^2=(x-x_k{)}^T(x-x_k)$。此外，由于 $\nabla f$ 满足 L-Lipschitz 条件，故 $\Vert {\nabla }^{2}f(x){\Vert }_2\le L,\forall x$，综合这些条件即可化简上式。

​ 进一步化简有：
$$\begin{array}{rl}\hat{f}(x)& \le f(x_k)+⟨\nabla f(x_k),(x-x_k)⟩+⟨\frac{L}{2}(x-x_k),x-x_k⟩\\ & =f(x_k)+⟨\nabla f(x_k)+\frac{L}{2}(x-x_k),x-x_k⟩\\ & =f(x_k)+\frac{L}{2}⟨\frac{2}{L}\nabla f(x_k)+(x-x_k),x-x_k⟩\\ & =f(x_k)+\frac{L}{2}⟨x-x_k+\frac{1}{L}\nabla f(x_k)+\frac{1}{L}\nabla f(x_k),x-x_k+\frac{1}{L}\nabla f(x_k)-\frac{1}{L}\nabla f(x_k)⟩\\ & =f(x_k)+\frac{L}{2}\Vert x-x_k+\frac{1}{L}\nabla f(x_k){\Vert }_2^2-\frac{1}{2L}\Vert \nabla f(x_k){\Vert }_2^2\\ & =\frac{L}{2}\Vert x-(x_k-\frac{1}{L}\nabla f(x_k)){\Vert }_2^2+\text{const},\end{array}\text{\hspace{1em}(42)}$$
上式中，第一步到第二步的化简利用了内积合并的线性性质：$⟨a,c⟩+⟨b,c⟩=⟨a+b,c⟩$，第四步到第五步的化简可以令 $a=x-x_k+\frac{1}{L}\nabla f(x_k)$，则 $⟨a+\delta ,a-\delta ⟩=\Vert a{\Vert }^2-\Vert \delta {\Vert }^2$，代入后乘上因子 $\frac{L}{2}$ 即可得到第五步。最后化简结果中的 $\text{const}$ 是常数。

​ 显然，式 (42) 的最小值在如下 $x_{k+1}$ 获得：
$$x_{k+1}=x_k-\frac{1}{L}\nabla f(x_k).\text{\hspace{1em}(43)}$$
​ 下面证明通过迭代上式可以使 $f(x)$ 逐步下降。

2.4.3 Descent Lemma 的证明

​ 对于函数 $f(x)$，由定积分的定义可得：
$$f(x_2)-f(x_1)={\int }_{x_1}^{x_2}{f}^{\prime }(x)dx.\text{\hspace{1em}(44)}$$
​ 设 $f(x)$ 二次连续可微，$x_k\in {\mathbb{R}}^n$，令 $g(t)=f(x_k+t(x-x_k))$，则：
$$\begin{array}{rl}g(0)& =f(x_k),\\ g(1)& =f(x).\end{array}\text{\hspace{1em}(45)}$$
​ 因此，由定积分的定义可得：
$$f(x)-f(x_k)=g(1)-g(0)={\int }_0^1{g}^{\prime }(t)dt.\text{\hspace{1em}(46)}$$
​ 对于 $g^{\prime }(t)$ 有：
$$\begin{array}{rl}{g}^{\prime }(t)& =\frac{\partial f(x_k+t(x-x_k))}{\partial t},\\ & =\nabla f(x_k+t(x-x_k){)}^T(x-x_k),\end{array}\text{\hspace{1em}(47)}$$
代入式 (46) 可得：
$$\begin{array}{rl}f(x)-f(x_k)& ={\int }_0^1\nabla f(x_k+t(x-x_k){)}^T(x-x_k)dt,\\ & ={\int }_0^1[\nabla f(x_k+t(x-x_k))-\nabla f(x_k)+\nabla f(x_k){]}^T(x-x_k)dt,\\ & ={\int }_0^1\nabla f(x_k{)}^T(x-x_k)dt+{\int }_0^1[\nabla f(x_k+t(x-x_k))-\nabla f(x_k){]}^T(x-x_k)dt,\end{array}\text{\hspace{1em}(48)}$$
对上式使用如下内积性质进行放缩：
$$x^{T}y=\Vert x\Vert \cdot \Vert y\Vert \cos \theta \le \Vert x\Vert \cdot \Vert y\Vert .\text{\hspace{1em}(49)}$$
​ 故式 (48) 可以放缩为如下形式：
$$f(x)-f(x_k)\le {\int }_0^1\nabla f(x_k{)}^T(x-x_k)dt+{\int }_0^1\Vert \nabla f(x_k+t(x-x_k))-\nabla f(x_k){\Vert }_2\cdot \Vert x-x_k{\Vert }_{2}dt.\text{\hspace{1em}(50)}$$
​ 接下来对上式使用 L-Lipschitz 常数进行放缩，可得到：
$$f(x)-f(x_k)\le {\int }_0^1\nabla f(x_k{)}^T(x-x_k)dt+\Vert x-x_k{\Vert }_2{\int }_0^{1}L\Vert t(x-x_k){\Vert }_{2}dt,\text{\hspace{1em}(51)}$$
对上式进行定积分，可化简为：
$$f(x)-f(x_k)=\nabla f(x_k{)}^T(x-x_k)+\frac{L}{2}\Vert (x-x_k){\Vert }_2^2,\text{\hspace{1em}(52)}$$
令
$$\hat{f}(x)=f(x_k)+\nabla f(x_k{)}^T(x-x_k)+\frac{L}{2}\Vert (x-x_k){\Vert }_2^2,\text{\hspace{1em}(53)}$$
则有 $f(x)\le \hat{f}(x)$，当 $x=x_k$ 时等号成立，该不等式被称为 descent lemma。

​ 由于 $f(x)$在$x=x_k-\frac{1}{L}\nabla f(x_k)$ 处取得最小值，因此以下不等式成立：
$$\hat{f}(x_k-\frac{1}{L}\nabla f(x_k))\le \hat{f}(x_k).\text{\hspace{1em}(54)}$$
​ 根据 descent lemma 的推导有 $f(x)\le \hat{f}(x)$ 恒成立，因此以下不等式也成立：
$$f(x_k-\frac{1}{L}\nabla f(x_k))\le \hat{f}(x_k-\frac{1}{L}\nabla f(x_k)),\text{\hspace{1em}(55)}$$
且当 $x=x_k$ 时，由式 (53) 可知 $f(x_k)=\hat{f}(x_k)$，因此：
$$f(x_k-\frac{1}{L}\nabla f(x_k))\le \hat{f}(x_k-\frac{1}{L}\nabla f(x_k))\le \hat{f}(x_k)=f(x_k).\text{\hspace{1em}(56)}$$
​ 综上所述，可以通过迭代下式使得 $f(x)$ 逐步下降：
$$x_{k+1}=x_k-\frac{1}{L}\nabla f(x_k).\text{\hspace{1em}(57)}$$

2.4.4 求解 Lasso

​ 上面的证明和推导都是针对 Lasso 的第一部分，Lasso 相对于线性回归只是加上了 L1 范数，故将式 (57) 应用在 Lasso 回归上，可得到如下一阶近似的迭代式：
$$x_{k+1}=\underset{x}{\mathrm{argmin}}\frac{L}{2}\Vert x-(x_k-\frac{1}{L}\nabla f(x_k)){\Vert }_2^2+\lambda \Vert x{\Vert }_1,\text{\hspace{1em}(58)}$$
令 $z=x_k-\frac{1}{L}\nabla f(x_k)$，$L=\frac{1}{t}$ 其实际上为 $f(x)$ 的梯度下降形式，可以将上式转换成如下形式：
$$x_{k+1}=\underset{x}{\mathrm{argmin}}\frac{1}{2}\Vert x-z{\Vert }_2^2+\lambda t\Vert x{\Vert }_1,\text{\hspace{1em}(59)}$$
上式也被称为近端算子（proximal operator），记为：
$$pro{x}_{t,h(\cdot )}(z)=\underset{x}{\mathrm{argmin}}\frac{1}{2}\Vert x-z{\Vert }_2^2+t\cdot h(x).\text{\hspace{1em}(60)}$$
​ 总体而言，首先对可导部分使用一阶泰勒展开和二次近似，然后再对整体目标函数应用近端操作。即先走一小步梯度下降，然后用近端算子来处理非光滑部分。特别地，当 $h(x)=\lambda \Vert x{\Vert }_1$时，近端算子有封闭解，称为 Soft Thresholding：
$$[pro{x}_{t\lambda \Vert \cdot {\Vert }_1}(z){]}_i=\{\begin{array}{ll}{z}_i-t\lambda ,& z_i>t\lambda ,\\ 0,& |z_i|\le t\lambda ,\\ z_i+t\lambda ,& z_i<-t\lambda ,\end{array}\text{\hspace{1em}(61)}$$
更简洁地，可以写为：
$$pro{x}_{t\lambda \Vert \cdot {\Vert }_1}(z)=\text{sign}(z)\cdot \max (|z|-t\lambda ,0).\text{\hspace{1em}(62)}$$
​ 对于式 (17) 的问题，近端梯度下降算法的一般迭代过程是先对 $f(x)$ 进行梯度下降求得 $z^{(k)}$，然后代入近端算子 $pro{x}_{t,h(\cdot )}(z^{(k)})$，然后迭代求解：
$$\begin{array}{rl}{z}^{(k)}& =x^{(k)}-t\nabla f(x^{(k)}),\\ x^{(k+1)}& =pro{x}_{t,h(\cdot )}(z^{(k)}).\end{array}\text{\hspace{1em}(63)}$$

2.4.5 具体求解例子

​ 假设数据设置有：
$$X=\left[\begin{array}{cc}1& 3\\ 3& 4\end{array}\right],y=\left[\begin{array}{c}1\\ 2\end{array}\right]w^{(0)}=\left[\begin{array}{c}0\\ 0\end{array}\right],\lambda =0.1,t=\mathrm{0.01.}\text{\hspace{1em}(64)}$$
​ 第一步需要计算梯度 $\nabla f(w)$：
$$\begin{array}{rl}f(w)& =\frac{1}{2n}\Vert Xw-y{\Vert }^2\\ & \Rightarrow \nabla f(w)=-\frac{1}{n}{X}^T(y-Xw)\\ & \Rightarrow \nabla f(w^{(0)})=\left[\begin{array}{c}-7\\ -10\end{array}\right].\end{array}\text{\hspace{1em}(65)}$$
​ 第二步进行梯度下降和近端操作：
$$z^{(0)}=w^{(0)}-t\nabla f(w^{(0)})=\left[\begin{array}{c}0.07\\ 0.10\end{array}\right],\text{\hspace{1em}(66)}$$
对每一维使用 Soft-thresholding：
$$\begin{array}{r}0.07\to \max (0.07-0.001,0)=0.069,\\ 0.10\to \max (0.10-0.001,0)=0.099,\end{array}\text{\hspace{1em}(67)}$$
因此第一轮更新的权重为：
$$w^{(1)}=prox(z^{(0)})=\left[\begin{array}{c}0.069\\ 0.099\end{array}\right].\text{\hspace{1em}(68)}$$

3. 最大后验估计

​ 最大后验估计是从贝叶斯的角度上的一种求解方式。假设 $w$ 有一个先验，满足 $w\sim \mathcal{N}(0,{\sigma }_0^2)$，噪声 $ϵ\sim \mathcal{N}(0,{\sigma }^2)$，目标变量 $y\sim \mathcal{N}(w^{T}x,{\sigma }^2)$。根据高斯分布的概率密度函数可写出似然函数和 $w$ 的先验分布：

p ( y ∣ w ) = ∏ i = 1 n 1 2 π σ exp ⁡ { − ( y i − w T x i ) 2 2 σ 2 } , (69) \begin{align} p(y|w)&=\prod_{i=1}^n\frac{1}{\sqrt{2\pi}\sigma}\exp\{-\frac{(y_i-w^Tx_i)^2}{2\sigma^2}\}, \end{align}\tag{69} p(y∣w)​=i=1∏n​2π ​σ1​exp{−2σ2(yi​−wTxi​)2​},​(69)
取对数并忽略常数：
$$\log p(y|w)\propto -\frac{1}{2{\sigma }^2}\sum _{i=1}^n(y_i-w^T{x}_i{)}^2=-\frac{1}{2{\sigma }^2}\Vert Xw-y{\Vert }^2.\text{\hspace{1em}(70)}$$
​ 先验项有：
p ( w ) = 1 ( 2 π σ 0 2 ) d / 2 exp ⁡ { − ∥ w ∥ 2 2 σ 0 2 } , (71) p(w)=\frac{1}{(2\pi \sigma_0^2)^{d/2}}\exp\{-\frac{\lVert w \rVert^2}{2\sigma_0^2}\},\tag{71} p(w)=(2πσ02​)d/21​exp{−2σ02​∥w∥2​},(71)
取对数并忽略常数：
$$\log p(w)\propto -\frac{1}{2{\sigma }_0^2}\Vert w{\Vert }^2.\text{\hspace{1em}(72)}$$
​ 最大后验等价于最小负对数后验：
$$\begin{array}{rl}\hat{w}& =\underset{w}{\mathrm{argmax}}\log p(y|w)+\log p(w)\\ & =\underset{w}{\mathrm{argmin}}\frac{1}{2{\sigma }^2}\Vert Xw-y{\Vert }^2+\frac{1}{2{\sigma }_0^2}\Vert w{\Vert }^2,\end{array}\text{\hspace{1em}(73)}$$
乘上 $2{\sigma }^2$ 不影响最小化目标：
$$\hat{w}=\underset{w}{\mathrm{argmin}}\Vert Xw-y{\Vert }^2+\lambda \Vert w{\Vert }^2,\lambda =\frac{{\sigma }^2}{{\sigma }_0^2}.\text{\hspace{1em}(74)}$$
​ 回顾最小二乘估计的求解结果，如果样本量 $N$ 没有远大于维度 $p$ 的时候，$X^{T}X$ 往往不可逆，也就是直观上理解的“过拟合”。在这种情况下，有三种方法：增加样本量、特征选择/降维、正则化。其中，正则化的框架为：

$$J(w)=L(w)+\lambda P(w).\text{\hspace{1em}(75)}$$
​ 在岭回顾中，其克服过拟合的思想为，一个半正定矩阵加上对角矩阵为正定矩阵，必定可逆。在上式框架中令 $P(w)=w^{T}w$，经最小二乘估计可求出 $w$ 为：
$$\hat{w}=\underset{x}{\mathrm{argmin}}J(w)=(X^{T}X+\lambda I{)}^{-1}{X}^{T}Y,\text{\hspace{1em}(76)}$$
上式等价于式子 (16)，观察式子 (16) 和式子 (74)，因此可以得出结论：噪声和先验为高斯分布的最大后验估计等价于正则化的最小二乘估计。

4. 自动求导

4.1 向量求导

图5 矩阵求导规则

​ 如果 $x$ 和 $y$ 都是标量，那么 $\partial y/\partial x$ 也是标量；如果 $y$ 是标量，$\mathbf{x}$ 是向量，那么 $\partial y/\partial \mathbf{x}$ 是行向量；如果 $\mathbf{y}$ 是向量，$x$ 是标量，那么 $\partial \mathbf{y}/\partial x$ 还是向量；如果 $\mathbf{x}$ 和 $\mathbf{y}$ 都是向量，那么 $\partial \mathbf{y}/\partial \mathbf{x}$ 是矩阵。具体而言，可以看作是向量 $\mathbf{y}$ 的每个元素对 $\mathbf{x}$ 求导，那么就会得到一个向量，向量中的每个元素为 $\partial y_i/\partial \mathbf{x}$，相当于标量对向量求导，变成一个个行向量，比如 $[\partial y_1/\partial x_1,\partial y_1/\partial x_2,\cdots ,\partial y_1/\partial x_n]$。

4.2 链式法则

图6 拓展到向量的链式法则

图7 链式法则例子

4.3 反向传播

图8 反向传播例子

​ 反向传播可以理解为神经网络加速梯度计算的方式。在图 8 中，为了求 $L$ 对 $w_1$ 的梯度，通过链式法则写成多个梯度相乘的形式。其中，计算的顺序是从后往前，黄色部分是已经计算的部分，不需要重新计算。

图9 正向传播计算图

​ 在计算机中，使用计算图来表示求导过程。如图 9 所示，这是一个正向传播的结构图，每个方块里面的计算都简单且可复制。

图10 反向传播计算图

​ 如图 10 所示，这是一个反向传播的计算图，黄色的式子就是通过之前的步骤得到的，绿色的式子就是在向前传播过程中得到的结果。

5. 线性回归代码实现

5.1 线性回归的从零实现

​ 代码来自为李沐《动手学深度学习》：https://zh.d2l.ai/chapter_linear-networks/linear-regression-scratch.html

import random
import torch
import matplotlib.pyplot as plt# 随机生成数据
def synthetic_data(w, b, num_examples):X = torch.normal(0, 1, (num_examples, len(w))) # 正态分布生成随机数y = torch.matmul(X, w) + b # 矩阵乘法y += torch.normal(0, 0.01, y.shape)return X, y.reshape((-1, 1)) # reshape 把一维张量转换成二维列向量true_w = torch.tensor([2, -3.4])
true_b = 4.2
features, labels = synthetic_data(true_w, true_b, 1000)# 获取批量数据
def data_iter(batch_size, features, labels):num_examples = len(features)indices = list(range(num_examples))random.shuffle(indices)for i in range(0, num_examples, batch_size):batch_indices = torch.tensor(indices[i:min(i + batch_size, num_examples)])yield features[batch_indices], labels[batch_indices] # 反复调用反复返回数据batch_size = 10
for X, y in data_iter(batch_size, features, labels):print(X,'\n', y)break# 初始化参数
w = torch.normal(0, 0.01, size=(2, 1), requires_grad = True) # 启用梯度跟踪
b = torch.zeros(1, requires_grad = True)# 定义线性回归模型
def linreg(X, w, b):return torch.matmul(X, w) + b# 定义损失函数
def square_loss(y_hat, y):return (y_hat - y.reshape(y_hat.shape)) ** 2 / 2# 定义优化算法
def sgd(params, lr, batch_size):# 小批量随机梯度下降with torch.no_grad(): # 禁用梯度计算for param in params:param -= lr * param.grad / batch_size # .grad 会储存张量在反向传播中计算的梯度param.grad.zero_() # 清零梯度，pytorch不会自动清零# 训练过程
lr = 0.03
num_epochs = 3
net = linreg
loss = square_lossfor epoch in range(num_epochs):for X, y in data_iter(batch_size, features, labels):l = loss(net(X, w, b), y) # 此时并不是一个标量l.sum().backward()'''对这个标量损失进行反向传播，计算所有 requires_grad=True 的参数（w 和 b）的梯度，并存储到它们的 .grad 属性中梯度计算基于链式法则，从损失函数反向传播到每个参数'''sgd([w, b], lr, batch_size)with torch.no_grad():train_l = loss(net(features, w, b), labels)print(f'epoch {epoch + 1}, loss {float(train_l.mean()):f}')print(f' w 的估计误差: {true_w - w.reshape(true_w.shape)}')
print(f' b 的估计误差: {true_b - b}')

5.2 线性回归的简洁实现

​ 代码来自为李沐《动手学深度学习》：https://zh.d2l.ai/chapter_linear-networks/linear-regression-concise.html

# 线性回归的简洁实现
import numpy as np
import torch
from torch.utils import data# 随机生成数据
def synthetic_data(w, b, num_examples):X = torch.normal(0, 1, (num_examples, len(w))) # 正态分布生成随机数y = torch.matmul(X, w) + b # 矩阵乘法y += torch.normal(0, 0.01, y.shape)return X, y.reshape((-1, 1)) # reshape 把一维张量转换成二维列向量true_w = torch.tensor([2, -3.4])
true_b = 4.2
features, labels = synthetic_data(true_w, true_b, 1000)# 调用框架中现有的 API 来读取数据
def load_array(data_arrays, batch_size, is_train = True):# 构造一个pytorch数据迭代器dataset = data.TensorDataset(*data_arrays) # 将传入的数据集转换成 Tensor 形式return data.DataLoader(dataset, batch_size, shuffle=is_train) # 随机抽取一些样本batch_size = 10
data_iter = load_array((features, labels), batch_size)
next(iter(data_iter))# 使用框架的预定义好的层
from torch import nn
net = nn.Sequential(nn.Linear(2, 1)) 
# 把包含一个线性层的网络放在 Sequential 容器里面，容器可以放置多个层
# 网络的输入特征维度为 2，输出维度为 1# 初始化模型参数
net[0].weight.data.normal_(0, 0.01) 
net[0].bias.data.fill_(0)# 计算均方误差
loss = nn.MSELoss()# 实例化 SGD 实例
trainer = torch.optim.SGD(net.parameters(), lr = 0.3)
# 把这个神经网络里所有要学习的参数都给net.parameters()，用优化器去更新它们# 训练部分
num_epochs = 3
for epoch in range(num_epochs):for X, y in data_iter:l = loss(net(X), y)trainer.zero_grad() # 梯度清零l.backward()trainer.step() # 用于根据当前参数的梯度更新模型参数l = loss(net(features), labels)print(f'epoch {epoch + 1}, loss {l:f}')

6. Softmax 回归

6.1 Softmax 与 交叉熵

​ 输出的匹配概率非负且和为 1：
$$\hat{\mathbf{y}}=\text{softmax}(\mathbf{o}),\text{\hspace{1em}(77)}$$
也可以以元素的方式进行表示：
$${\hat{y}}_i=\frac{\exp o_i}{{\sum }_k\exp (o_k)}.\text{\hspace{1em}(78)}$$
​ 交叉熵常用来衡量两个概率的区别：
$$H(\mathbf{p},\mathbf{q})=\sum _i-p_i\log (q_i),\text{\hspace{1em}(79)}$$
上式的 $i$ 不是样本索引，而是类别索引。将上式作为 Softmax 的损失：
$$\begin{array}{rl}l(\mathbf{y},\hat{\mathbf{y}})& =-\sum _i{y}_i\log {\hat{y}}_i\\ & =-y_y\log {\hat{y}}_y-\sum _{i\ne y}{y}_i\log {\hat{y}}_i\\ & =-\log {\hat{y}}_y.\end{array}\text{\hspace{1em}(80)}$$
​ 对于上式为什么可以这样化简，这是因为真实标签 $\mathbf{y}$ 是一个 ont-hot 向量：
$$y_i=\{\begin{array}{ll}1,& i=y,\\ 0,& \text{otherwise}.\end{array}\text{\hspace{1em}(81)}$$
​ 因此，公式 (80) 只有在 $i=y$ 时不为 0，$y_y$ 表示真实类别。

​ 进一步地，梯度是真实概率和预测概率的区别：
$$\begin{array}{rl}\frac{\partial }{\partial o_i}l(\mathbf{y},\hat{\mathbf{y}})& =-\sum _j{y}_j\cdot \frac{\partial }{\partial o_i}\log \widehat{y_j}\\ & =\sum _j{y}_j\cdot \frac{\partial }{\partial o_i}\log (\frac{\exp (o_j)}{{\sum }_k\exp (o_k)})\\ & =\sum _j{y}_j\cdot \frac{\partial }{\partial o_i}(o_j-\log (\sum _k\exp (o_k)))，\end{array}\text{\hspace{1em}(82)}$$
对 $o_j$ 求导的部分，如果 $i=j$，则结果为 1；否则结果为 0。对对数项求导的部分，使用链式法则有：
$$\begin{array}{rl}\frac{\partial }{\partial o_i}\log (\sum _k\exp (o_k))& =\frac{1}{{\sum }_k\exp (o_k)}\cdot \frac{\partial }{\partial o_i}(\sum _k\exp (o_k))\\ & =\frac{1}{{\sum }_k\exp (o_k)}\cdot \exp (o_i)\\ & ={\hat{y}}_i.\end{array}\text{\hspace{1em}(83)}$$
​ 综合上述两项的结果，求偏导的结果为：
$$\frac{\partial }{\partial o_i}\log \widehat{y_j}=\{\begin{array}{ll}1-\widehat{y_i},& \text{if}i=j,\\ -\widehat{y_i},& \text{if}i\ne j.\end{array}\text{\hspace{1em}(84)}$$
​ 将上式代入式 (82) 可得：
$$\begin{array}{rl}\frac{\partial }{\partial o_i}l(\mathbf{y},\hat{\mathbf{y}})& =-\sum _j{y}_j\cdot ({\delta }_{ij}-\widehat{y_i})\\ & =\sum _j{y}_j{\delta }_{ij}+\sum _j{y}_j\widehat{y_i}\\ & =-y_i+\widehat{y_i}.\end{array}\text{\hspace{1em}(85)}$$

6.2 损失函数

​ 下面给出几种损失函数的可视化，都是真实值为 0 的情况。

图11 均方损失：蓝线为损失函数，绿线为似然函数，橙色为梯度

​ 图 11 想表达的是，当真实值和预测值相差较大的时候，梯度比较大，下降比较快；反之，真实值和预测值相差较小的时候，梯度较小。

图12 绝对值损失函数

​ 绝对值损失函数想表达的是，无论误差有多大，梯度始终是常数，权重更新不会太大，带来稳定性，但也具有不平滑和不可导的缺陷。比如，当预测值和真实值较为接近时，在迭代后期，会变得不稳定。

图13 结合均方误差和绝对值误差的损失函数

参考

[1] Machine Learning基础：正则化项（Regularization）

[2] 【西瓜书】 第11章 特征选择与稀疏学习

[3] 近端梯度(Proximal Gradient)下降算法的过程以及理解|ISTA算法|LASSO问题 - 红泥小火炉的文章 - 知乎

[4] Lasso回归与近端梯度下降推导与实现Posted by Welt Xing on September 30, 2021

[5] 【西瓜书】 第11章 特征选择与稀疏学习笔记

[6] 利普希茨连续的几何意义是什么？怎么较好的理解它呢？ - cielo的回答 - 知乎

[7] L1范数中为什么需要证明目标函数的梯度满足Lipschitz（利普希茨）连续？ - Spectre的回答 - 知乎

[8] 【机器学习】【白板推导系列】【合集 1～33】

[9] 08 线性回归 + 基础优化算法【动手学深度学习v2】

[10] [5分钟深度学习] #02 反向传播算法