线性回归与岭回归

线性方程扰动问题

在讨论一般的凸优化问题之前，我们先从一个基本的线性方程扰动问题开始研究。

  • 考虑一个系统：Ax⃗=y⃗A\vec{x}=\vec{y}，其中矩阵AA可逆（可看作系统内部属性），y⃗\vec{y}为观测值，而x⃗\vec{x}则是从观测值反推出的未知量（可看作系统的解）。
    • 我们需要考察当观测值有扰动时，对应的解扰动的度量。（也称为敏感性分析）具体而言，对于一个加入扰动的系统： A(x⃗+δ⃗x)=y⃗+δ⃗yA(\vec{x}+\vec{\delta}_x)=\vec{y}+\vec{\delta}_y 我们希望在给定∥y⃗∥2\|\vec{y}\|_2和∥δ⃗y∥2\|\vec{\delta}_y\|_2的前提下对解的变化率∥δ⃗x∥2∥x⃗∥2\dfrac{\|\vec{\delta}_x\|_2}{\|\vec{x}\|_2}进行评估。特别地，我们希望强化系统的稳健性，即让这个变化率尽量的小。
  • 下面我们对这个变化率给出一个上界：
    • 首先求∥δ⃗x∥2\|\vec{\delta}_x\|_2的上界： Aδ⃗x=δ⃗y⟹δ⃗x=A−1δ⃗y⟹∥δ⃗x∥2=∥A−1δ⃗y∥2≤max⁡z⃗∈Rn∥z⃗∥2=∥δ⃗y∥2∥A−1z⃗∥2=max⁡z⃗∈Rn∥z⃗∥2=1∥A−1z⃗∥2∥δ⃗y∥2=∥A−1∥2∥δ⃗y∥2\begin{aligned} A\vec{\delta}_x&=\vec{\delta}_y\Longrightarrow \vec{\delta}_x=A^{-1}\vec{\delta}_y\\ \Longrightarrow \|\vec{\delta}_x\|_2&=\|A^{-1}\vec{\delta}_y\|_2\leq\max_{\substack{\vec{z}\in\R^n\\\|\vec{z}\|_2=\|\vec{\delta}_y\|_2}}\|A^{-1}\vec{z}\|_2=\max_{\substack{\vec{z}\in\R^n\\\|\vec{z}\|_2=1}}\|A^{-1}\vec{z}\|_2\|\vec{\delta}_y\|_2=\|A^{-1}\|_2\|\vec{\delta}_y\|_2 \end{aligned} 其中∥A−1∥2\|A^{-1}\|_2表示矩阵的谱范数。
    • 然后求∥x⃗∥2\|\vec{x}\|_2的下界： ∥y⃗∥2=∥Ax⃗∥2≤∥A∥2∥x⃗∥2⟹∥x⃗∥2≥∥y⃗∥2∥A∥2\|\vec{y}\|_2=\|A\vec{x}\|_2\leq\|A\|_2\|\vec{x}\|_2\Longrightarrow \|\vec{x}\|_2\geq\frac{\|\vec{y}\|_2}{\|A\|_2}
    • 二者结合得到总体的上界： ∥δ⃗x∥2∥x⃗∥2≤∥A∥2∥A−1∥2∥δ⃗y∥2∥y⃗∥2.\frac{\|\vec{\delta}_x\|_2}{\|\vec{x}\|_2}\leq \|A\|_2\|A^{-1}\|_2\frac{\|\vec{\delta}_y\|_2}{\|\vec{y}\|_2}.
  • 由此可知解的扰动受到观测值扰动的控制。而由谱范数与SVD的结论，∥A∥2\|A\|_2可用σmax⁡{A}\sigma_{\max}\{A\}表示，而∥A−1∥2\|A^{-1}\|_2可用σmax⁡{A−1}=1σmin⁡{A}\sigma_{\max}\{A^{-1}\}=\dfrac{1}{\sigma_{\min}\{A\}}表示（AA与A−1A^{-1}奇异值互为倒数）。于是可以定义 κ(A)≔∥A∥2∥A−1∥2=σmax⁡{A}σmin⁡{A}\kappa(A)\coloneqq \|A\|_2\|A^{-1}\|_2=\frac{\sigma_{\max}\{A\}}{\sigma_{\min}\{A\}} 称κ(A)\kappa(A)为AA的条件数（condition number）。（当然，条件数的定义也可以推广到不可逆矩阵，此时值为无穷大）
    • 由定义可知，条件数刻画了观测值扰动对解扰动影响的相对比率。条件数越大，观测值的微小扰动就越可能导致解的大幅扰动。
  • 另外，由于现实情况中矩阵AA一般不是方阵（行数一般大于列数），此时我们使用正规方程（normal equation）替代： A⊤Ax⃗=A⊤y⃗A^\top A\vec{x}=A^\top\vec{y} 此时使用的条件数变为κ(A⊤A)\kappa(A^\top A)。又因为A⊤AA^\top A是半正定对称矩阵，其特征值与奇异值相等，所以条件数可以表示为 κ(A⊤A)=λmax⁡{A⊤A}λmin⁡{A⊤A}.\kappa(A^\top A)=\frac{\lambda_{\max}\{A^\top A\}}{\lambda_{\min}\{A^\top A\}}.

岭回归

  • 在最小二乘问题中，如果A⊤AA^\top A的条件数很大（接近奇异），那么这会导致计算出的参数估计值与真实值相去甚远。
  • 对此，我们考虑一种处理方法：对A⊤AA^\top A的所有特征值统一加一个常量，这样特征数就可以缩小（当λmin⁡{A⊤A}\lambda_{\min}\{A^\top A\}接近00时缩小的幅度非常可观）。具体而言，我们将原有的系统改造为 (A⊤A+λI)x⃗=A⊤y⃗(A^\top A+\lambda I)\vec{x}=A^\top\vec{y} 其中λ>0\lambda > 0。实际上，这恰好对应下述优化问题： min⁡x⃗∈Rn{∥Ax⃗−y⃗∥22+λ∥x⃗∥22}⟹x⃗∗=(A⊤A+λI)−1A⊤y⃗\min_{\vec{x}\in\R^n}\left\{\|A\vec{x}-\vec{y}\|_2^2+\lambda\|\vec{x}\|_2^2\right\}\Longrightarrow \vec{x}^*=(A^\top A+\lambda I)^{-1}A^\top\vec{y} 上述优化也被称为岭回归（Ridge Regression）。下面给出这个优化问题解的推导： ∇x⃗{∥Ax⃗−y⃗∥22+λ∥x⃗∥22}=∇x⃗{x⃗⊤A⊤Ax⃗−2y⃗⊤Ax⃗+y⃗⊤y⃗+λx⃗⊤x⃗}=2A⊤Ax⃗−2A⊤y⃗+2λx⃗=2(A⊤A+λI)x⃗−2A⊤y⃗.\begin{align*} \nabla_{\vec{x}} \left\{ \|A\vec{x} - \vec{y}\|_2^2 + \lambda \|\vec{x}\|_2^2 \right\} &= \nabla_{\vec{x}} \{\vec{x}^\top A^\top A \vec{x} - 2 \vec{y}^\top A \vec{x} + \vec{y}^\top \vec{y} + \lambda \vec{x}^\top \vec{x}\} \\ &= 2A^\top A \vec{x} - 2A^\top \vec{y} + 2\lambda \vec{x} \\ &= 2(A^\top A + \lambda I) \vec{x} - 2A^\top \vec{y}. \end{align*} 而A⊤A+λIA^\top A + \lambda I是正定矩阵，因此可逆。所以令梯度取零点就得到上述优化问题的解。【之后我们会证明凸函数的最值一定满足梯度为零】
    补充

    当然，还可以将岭回归问题转化为最小二乘问题：

    min⁡x⃗∈Rn{∥Ax⃗−y⃗∥22+λ∥x⃗∥22}⟹min⁡x⃗∈Rn{∥[AλI]x⃗−[y⃗0⃗]∥22}\min_{\vec{x}\in\R^n}\left\{\|A\vec{x}-\vec{y}\|_2^2+\lambda\|\vec{x}\|_2^2\right\}\Longrightarrow \min_{\vec{x}\in\R^n}\left\{\left\| \begin{bmatrix} A \\ \sqrt{\lambda}I \end{bmatrix} \vec{x} - \begin{bmatrix} \vec{y} \\ \vec{0} \end{bmatrix} \right\|_2^2\right\}

    然后套用最小二乘问题的解法得到：

    x⃗=([AλI]⊤[AλI])−1[AλI]⊤[y⃗0⃗]=([A⊤λI][AλI])−1[A⊤λI][y⃗0⃗]=(A⊤A+λI)−1A⊤y⃗.\begin{aligned} \vec{x} &= \left( \begin{bmatrix} A \\ \sqrt{\lambda}I \end{bmatrix}^\top \begin{bmatrix} A \\ \sqrt{\lambda}I \end{bmatrix} \right)^{-1} \begin{bmatrix} A \\ \sqrt{\lambda}I \end{bmatrix}^\top \begin{bmatrix} \vec{y} \\ \vec{0} \end{bmatrix} \\ &= \left( \begin{bmatrix} A^\top & \sqrt{\lambda}I \end{bmatrix} \begin{bmatrix} A \\ \sqrt{\lambda}I \end{bmatrix} \right)^{-1} \begin{bmatrix} A^\top & \sqrt{\lambda}I \end{bmatrix} \begin{bmatrix} \vec{y} \\ \vec{0} \end{bmatrix} \\ &= ( A^\top A + \lambda I )^{-1} A^\top \vec{y}. \end{aligned}

    这相当于在原有Ax⃗=y⃗A\vec{x}=\vec{y}的基础上增加了额外的信息：λIx⃗≈0⃗\sqrt{\lambda}I\vec{x}\approx\vec{0}（x⃗\vec{x}的每个分量大小都被控制）

    岭回归中的λ∥x⃗∥22\lambda\|\vec{x}\|_2^2也被称作正则化项，这相当于给优化问题增加了一个约束或惩罚。（在解的误差与解本身的范数大小之间进行权衡）
  • 下面我们讨论岭回归与SVD（PCA）之间的关系：设A=UΣV⊤A=U\Sigma V^\top，那么有 x⃗∗=(A⊤A+λI)−1A⊤y⃗=(VΣ⊤ΣV⊤+λI)−1VΣ⊤U⊤y⃗=(V(Σ⊤Σ+λI)V⊤)−1VΣ⊤U⊤y⃗=V(Σ⊤Σ+λI)−1V⊤VΣ⊤U⊤y⃗=V(Σ⊤Σ+λI)−1Σ⊤U⊤y⃗=V[σ1{A}σ1{A}2+λ⋱σn{A}σn{A}2+λ]U⊤y⃗=∑i=1rσi{A}σi{A}2+λ(u⃗i⊤y⃗)⋅v⃗i.\begin{aligned} \vec{x}^*&=(A^\top A+\lambda I)^{-1}A^\top\vec{y}\\ &= (V\Sigma^\top \Sigma V^\top + \lambda I)^{-1} V\Sigma^\top U^\top \vec{y} \\ &= (V (\Sigma^\top \Sigma + \lambda I) V^\top)^{-1} V\Sigma^\top U^\top \vec{y} \\ &= V (\Sigma^\top \Sigma + \lambda I)^{-1} V^\top V\Sigma^\top U^\top \vec{y} \\ &= V (\Sigma^\top \Sigma + \lambda I)^{-1} \Sigma^\top U^\top \vec{y} \\ &= V \begin{bmatrix}\dfrac{\sigma_1\{A\}}{\sigma_1\{A\}^2+\lambda}&&\\&\ddots&\\&&\dfrac{\sigma_n\{A\}}{\sigma_n\{A\}^2+\lambda}\end{bmatrix} U^\top \vec{y}\\ &=\sum_{i=1}^{r} \frac{\sigma_i\{A\}}{\sigma_i\{A\}^2 + \lambda} (\vec{u}_i^\top \vec{y}) \cdot \vec{v}_i. \end{aligned} 可以看到，λ\lambda的引入让原有最小二乘的解向零点拉近，不同奇异值对应的分量拉近的幅度不同。
    • 实际上，当λ\lambda增大时，大奇异值对应的分量（如v⃗1\vec{v}_1）影响会比小奇异值对应分量更小，这与主成分分析得到的效果类似。因此从某种程度上，岭回归可看作一种软化的PCA。

Tikhonov 回归

  • 上面的岭回归中，我们对∥x⃗∥22\|\vec{x}\|_2^2进行了惩罚。实际上，我们也可以考虑让最小二乘解尽量接近一个特定向量x⃗0\vec{x}_0，即对∥x⃗−x⃗0∥22\|\vec{x}-\vec{x}_0\|_2^2进行惩罚。此时目标优化函数就变为 min⁡x⃗∈Rn{∥Ax⃗−y⃗∥22+λ∥x⃗−x⃗0∥22}\min_{\vec{x}\in\R^n}\left\{\|A\vec{x}-\vec{y}\|_2^2+\lambda\|\vec{x}-\vec{x}_0\|_2^2\right\} 更一般地，我们可以对向量的不同分量实施不同程度的惩罚。引入对角矩阵W1∈Rm×m,W2∈Rn×nW_1\in\R^{m\times m},W_2\in\R^{n\times n}，优化问题变为 min⁡x⃗∈Rn{∥W1(Ax⃗−y⃗)∥22+∥W2(x⃗−x⃗0)∥22}\min_{\vec{x}\in\R^n}\left\{\|W_1(A\vec{x}-\vec{y})\|_2^2+\|W_2(\vec{x}-\vec{x}_0)\|_2^2\right\} 此时岭回归即为W1=I,W2=λI,x⃗0=0⃗W_1=I,W_2=\lambda I,\vec{x}_0=\vec{0}的特殊情况。上述问题也称为Tikhonov回归，其解为 x⃗∗=(A⊤W12A+W22)−1(A⊤W12y⃗+W22x⃗0).\vec{x}^* = (A^\top W_1^2 A + W_2^2)^{-1} (A^\top W_1^2 \vec{y} + W_2^2 \vec{x}_0).
  • 下面我们将其与极大似然估计/极大后验估计进行联系：设 A=[a⃗1⊤⋮a⃗m⊤],y⃗=[y1⋮ym]A=\begin{bmatrix}\vec{a}_1^\top\\\vdots\\\vec{a}_m^\top\end{bmatrix},\vec{y}=\begin{bmatrix}y_1\\\vdots\\y_m\end{bmatrix} 考虑概率模型y⃗=Ax⃗+w⃗\vec{y}=A\vec{x}+\vec{w}，其中 w⃗∼N(0⃗,Σw⃗),Σw⃗=diag(σ12,⋯ ,σm2).\vec{w}\sim\mathcal N(\vec{0},\Sigma_{\vec{w}}),\Sigma_{\vec{w}}=\text{diag}(\sigma_1^2,\cdots,\sigma_m^2). 另设px⃗(y⃗)p_{\vec{x}}(\vec{y})表示y⃗\vec{y}的在给定x⃗\vec{x}条件下的密度函数（即似然函数）。
    • 首先给出MLE的一个结论： arg max⁡x⃗∈Rnpx⃗(y⃗)=arg min⁡x⃗∈Rn∥Σw⃗−1/2(Ax⃗−y⃗)∥22.\argmax_{\vec{x} \in \mathbb{R}^n} p_{\vec{x}}(\vec{y}) = \argmin_{\vec{x} \in \mathbb{R}^n} \left\| \Sigma_{\vec{w}}^{-1/2} (A\vec{x} - \vec{y}) \right\|_2^2. 即极大似然估计与W1=Σw⃗−1/2,W2=0W_1=\Sigma_{\vec{w}}^{-1/2},W_2=\mathbf{0}的Tikhonov回归结果等价。其证明如下： arg max⁡x⃗∈Rnpx⃗(y⃗)=arg max⁡x⃗∈Rnlog⁡(px⃗(y⃗))=arg max⁡x⃗∈Rn∑i=1mlog⁡(px⃗(yi))=arg max⁡x⃗∈Rn∑i=1mlog⁡(12πσi2exp⁡(−(yi−a⃗i⊤x⃗)22σi2))=arg max⁡x⃗∈Rn{∑i=1mlog⁡(exp⁡(−(yi−a⃗i⊤x⃗)22σi2))}=arg min⁡x⃗∈Rn∑i=1m(yi−a⃗i⊤x⃗)2σi2=arg min⁡x⃗∈Rn∥Σw⃗−1/2(Ax⃗−y⃗)∥22.\begin{aligned} \argmax_{\vec{x} \in \mathbb{R}^n} p_{\vec{x}}(\vec{y}) &= \argmax_{\vec{x} \in \mathbb{R}^n} \log(p_{\vec{x}}(\vec{y})) \\ &= \argmax_{\vec{x} \in \mathbb{R}^n} \sum_{i=1}^m \log(p_{\vec{x}}(y_i)) \\ &= \argmax_{\vec{x} \in \mathbb{R}^n} \sum_{i=1}^m \log\left( \frac{1}{\sqrt{2\pi\sigma_i^2}} \exp\left( -\frac{(y_i - \vec{a}_i^\top \vec{x})^2}{2\sigma_i^2} \right) \right)\\ &= \argmax_{\vec{x} \in \mathbb{R}^n} \left\{ \sum_{i=1}^m \log \left( \exp \left( -\frac{(y_i - \vec{a}_i^\top \vec{x})^2}{2\sigma_i^2} \right) \right) \right\} \\ &= \argmin_{\vec{x} \in \mathbb{R}^n} \sum_{i=1}^m \frac{(y_i - \vec{a}_i^\top \vec{x})^2}{\sigma_i^2} = \argmin_{\vec{x} \in \mathbb{R}^n} \left\| \Sigma_{\vec{w}}^{-1/2} (A\vec{x} - \vec{y}) \right\|_2^2. \end{aligned} 特别地，如果w⃗\vec{w}各分量独立同分布，那么极大似然估计等价于最小二乘法。
    • 然后再给出极大后验估计的结论：在上述模型的基础上进一步假设x⃗=x⃗0+v⃗\vec{x}=\vec{x}_0+\vec{v}，其中 x⃗=[x1⋮xn],x⃗0=[(x⃗0)1⋮(x⃗0)n],v⃗∼N(0⃗,Σv⃗),Σv⃗=diag(τ12,⋯ ,τm2)\vec{x}=\begin{bmatrix}x_1\\\vdots\\x_n\end{bmatrix},\vec{x}_0=\begin{bmatrix}(\vec{x}_0)_1\\\vdots\\(\vec{x}_0)_n\end{bmatrix},\vec{v}\sim\mathcal{N}(\vec{0},\Sigma_{\vec{v}}),\Sigma_{\vec{v}}=\text{diag}(\tau_1^2,\cdots,\tau_m^2) 那么有 arg max⁡x⃗∈Rnp(x⃗∣y⃗)=arg min⁡x⃗∈Rn{∥Σw⃗−1/2(Ax⃗−y⃗)∥22+∥Σv⃗−1/2(x⃗−x⃗0)∥22}.\argmax_{\vec{x} \in \mathbb{R}^n} p(\vec{x}|\vec{y}) = \argmin_{\vec{x} \in \mathbb{R}^n} \left\{ \left\| \Sigma_{\vec{w}}^{-1/2} \left( A \vec{x} - \vec{y} \right) \right\|_2^2 + \left\| \Sigma_{\vec{v}}^{-1/2} \left( \vec{x} - \vec{x}_0 \right) \right\|_2^2 \right\}. 其中p(x⃗∣y⃗)p(\vec{x}|\vec{y})表示后验概率。其证明如下： arg max⁡x⃗∈Rnp(x⃗∣y⃗)=arg max⁡x⃗∈Rnlog⁡(p(x⃗∣y⃗))=arg max⁡x⃗∈Rnlog⁡(p(y⃗∣x⃗)p(x⃗)p(y⃗))=arg max⁡x⃗∈Rn{log⁡(p(y⃗∣x⃗))+log⁡(p(x⃗))}=arg max⁡x⃗∈Rn{∑i=1mlog⁡(p(yi∣x⃗))+∑j=1nlog⁡(p(xj))}=arg max⁡x⃗∈Rn{∑i=1mlog⁡(12πσi2exp⁡(−(yi−a⃗i⊤x⃗)22σi2))+∑j=1nlog⁡(12πτj2exp⁡(−(xj−(x⃗0)j)22τj2))}=arg max⁡x⃗∈Rn{∑i=1m(−(yi−a⃗i⊤x⃗)22σi2)+∑j=1n(−(xj−(x⃗0)j)22τj2)}=arg min⁡x⃗∈Rn{∑i=1m((yi−a⃗i⊤x⃗)2σi2)+∑j=1n((xj−(x⃗0)j)2τj2)}=arg min⁡x⃗∈Rn{∥Σw⃗−1/2(Ax⃗−y⃗)∥22+∥Σv⃗−1/2(x⃗−x⃗0)∥22}.\begin{aligned} \argmax_{\vec{x} \in \mathbb{R}^n} p(\vec{x}|\vec{y}) &= \argmax_{\vec{x} \in \mathbb{R}^n} \log(p(\vec{x}|\vec{y})) \\ &= \argmax_{\vec{x} \in \mathbb{R}^n} \log\left(\frac{p(\vec{y}|\vec{x})p(\vec{x})}{p(\vec{y})}\right) \\ &= \argmax_{\vec{x} \in \mathbb{R}^n} \left\{ \log(p(\vec{y}|\vec{x})) + \log(p(\vec{x})) \right\} \\ &= \argmax_{\vec{x} \in \mathbb{R}^n} \left\{ \sum_{i=1}^m \log(p(y_i|\vec{x})) + \sum_{j=1}^n \log(p(x_j)) \right\} \\ &= \argmax_{\vec{x} \in \mathbb{R}^n} \left\{ \sum_{i=1}^m \log\left(\frac{1}{\sqrt{2\pi\sigma_i^2}} \exp\left(-\frac{(y_i - \vec{a}_i^\top \vec{x})^2}{2\sigma_i^2}\right)\right) + \sum_{j=1}^n \log\left(\frac{1}{\sqrt{2\pi\tau_j^2}} \exp\left(-\frac{(x_j - (\vec{x}_0)_j)^2}{2\tau_j^2}\right)\right) \right\} \\ &= \argmax_{\vec{x} \in \mathbb{R}^n} \left\{ \sum_{i=1}^m \left(-\frac{(y_i - \vec{a}_i^\top \vec{x})^2}{2\sigma_i^2}\right) + \sum_{j=1}^n \left(-\frac{(x_j - (\vec{x}_0)_j)^2}{2\tau_j^2}\right) \right\} \\ &= \argmin_{\vec{x} \in \mathbb{R}^n} \left\{ \sum_{i=1}^m \left(\frac{(y_i - \vec{a}_i^\top \vec{x})^2}{\sigma_i^2}\right) + \sum_{j=1}^n \left(\frac{(x_j - (\vec{x}_0)_j)^2}{\tau_j^2}\right) \right\} \\ &= \argmin_{\vec{x} \in \mathbb{R}^n} \left\{ \left\| \Sigma_{\vec{w}}^{-1/2} \left( A \vec{x} - \vec{y} \right) \right\|_2^2 + \left\| \Sigma_{\vec{v}}^{-1/2} \left( \vec{x} - \vec{x}_0 \right) \right\|_2^2 \right\}. \end{aligned}