一文看懂最小二乘法

• 参考：
  1. 如何理解最小二乘法？
  2. 最小二乘法小结
  3. 最小二乘法的原理理解

---

1 回归与函数拟合

• 回归分析：指研究两组相互依赖的随机变量 $Y_1,Y_2,...Y_n$（常称为因变量） 和 $X_1,X_2,...X_m$（常称为自变量）之间关系的统计分析方法
  1. 按照自变量的多少，分为一元回归和多元回归
  2. 按照因变量的多少，可分为简单回归和多重回归
  3. 按照自变量和因变量之间的关系类型，可分为线性回归和非线性回归
• 回归建立了自变量和因变量间的关系，因此常用于拟合函数曲线。机器学习中，回归属于监督学习，是一种预测性建模技术，说白了其实就是拟合一条曲线，利用它可以做回归任务，如线性回归；也能做分类任务，如逻辑回归（此时看作拟合分类平面曲线）
• 补充几个和函数曲线拟合相关的概念
  1. 拟合函数/估值函数：给定来自真实分布的一批样本，拟合后得到的函数
  2. 损失函数 和 目标函数：做函数拟合时，常将其转化为一个优化问题，认为和数据匹配越好的函数越优，这时需要一个指示拟合程度的指标来指导优化方向。给定训练样本，该指标是从模型参数到拟合程度的映射，称之为损失函数或目标函数。损失函数越小/目标函数越大，拟合得就越好

     注：能列出损失函数/目标函数，就意味着模型形式（函数参数数量，线性关系还是二次关系…）已经确定

  3. 残差：针对某个样本而言， 拟合值与真实值的差。通常把关于的单个样本残差推广到整个训练集来构造损失函数或目标函数，有以下选择

     1. 残差和： 指所有训练样本拟合值与真实值的差值的和，有正有负会存在抵消的情况，不能反应真实误差
     2. 残差绝对值和：这可以解决残差和有正有负的问题，但是绝对值在后续的求导会异常麻烦
     3. 残差的平方和：这可以避免上述两个问题，运用广泛，它再取个均值就是机器学习中常用的 MSE 最小均方损失。其问题在于对异常值非常敏感，会被脏数据影响

  4. 最优化问题：指给定目标函数/损失函数时，通过决定某些可选择的变量应该取何值，使所给定目标函数最大化/损失函数最小化的问题。根据有无约束条件，又分为 无约束最优化问题 和 约束最优化问题
  5. 最优化方法：指解决最优化问题的方法

2 最小二乘法

• 最小二乘法/最小平方法：一种最优化方法，在模型给定的情况下，使用样本残差的平方和作为损失函数，通过令偏导数为 0 直接求解无约束最优化问题，得到最优模型参数的解析解。最小二乘法可用于曲线拟合，也能通过最小化能量或最大化熵来表达其他一些优化问题

2.1 适用模型

• 最小二乘法是针对 线性模型 的 多元简单回归 任务设计的方法，即要求标签为 1 维，且模型对样本特征向量 $x_i=[x_1,x_2,...,x_m{]}^{\top }$ 给出的预测为 $$f(x,w)=w_0+w_1{x}_1+w_2{x}_2+...+w_m{x}_m$$ 为了表示简便，通常增设一个 $x_0=1$，模型表示为
  $$f(x,w)=\sum _{i=0}^m{w}_i{x}_i=x^{\top }w$$
• 非线性模型无法使用最小二乘法估计参数，这是因为非线性的拟合函数里面可能有很复杂的结构，比如求和，积分，级数等，此时就算直接令导数为 0 得到方程，也可能很难求出方程的解析解
• 虽然针对线性模型设计，但并不意味着不能用最小二乘法优化非线性模型参数，可以先把非线性模型表示为非线性部分的线性组合，再使用最小二乘法，如 $$f(x,w)=w_0+w_1{f}_1(x)+w_2{f}_2(x)+...+w_m{f}_m(x)$$ 其中每个 $f_i(x)$ 都是非线性部分，比如多项式回归就属于这种情况
• 最小二乘法在求解时，只要有解，一定是全局最优解，参考：最小二乘法求解方程组为什么不会陷入局部最优解

2.2 求解方法

• 先从最常见的多元简单回归角度给出形式化描述：给定样本集
  $$\mathcal{D}=\{(x_1,y_1),(x_2,y_2),...,(x_n,y_n)\}$$设模型为 $f$，样本特征向量 $x_i=[1,x_{i1},x_{i1},...,x_{im}{]}^{\top }$，参数向量 $w=[w_0,w_1,w_2,...,w_m{]}^{\top }$，为表示简便，设 $p=m+1,x_i\in {\mathbb{R}}^p,y_i\in \mathbb{R}$，样本集 $X_{n\times p}=[x_1^{\top },x_2^{\top },...,x_n^{\top }{]}^{\top }$，标记集 $y=[y_1,y_2,...,y_n{]}^{\top }$，模型对于样本 $x_i$ 给出预测 ${\hat{y}}_i=f(x_i,w)=x_i^{\top }w$，对整个样本集有
  $$Xw=\left[\begin{array}{c}{x}_1^{\top }\\ x_2^{\top }\\ ⋮\\ x_n^{\top }\end{array}\right]w=\left[\begin{array}{c}{\hat{y}}_1\\ {\hat{y}}_2\\ ⋮\\ {\hat{y}}_n\end{array}\right]=\hat{y}$$ 使用残差平方和作为损失函数，即 $l(w)={\sum }_{i=1}^n(f(x_i,w)-y_i{)}^2$，求以下最优化问题
  $$\underset{\mathbf{w}}{\mathrm{argmin}}l(w)=\underset{\mathbf{w}}{\mathrm{argmin}}\sum _{i=1}^n(f(x_i,w)-y_i{)}^2=\underset{\mathbf{w}}{\mathrm{argmin}}(Xw-y{)}^{\top }(Xw-y)$$ 求解时，通过令偏导数为 0 直接计算，即令 $\frac{\partial l(\mathbf{w})}{\partial w_i}=0$ 直接解出参数 $w_i$

• 我们可以对每个参数令偏导数为 0 来求解，但是这样计算量太大了，所以直接对参数向量求导。首先展开目标函数
  $$\begin{array}{rl}(Xw-y{)}^{\top }(Xw-y)& =(w^{\top }{X}^{\top }-y^{\top })(Xw-y)\\ & =w^{\top }{X}^{\top }Xw-w^{\top }{X}^{\top }y-y^{\top }Xw+y^{\top }y\\ & =w^{\top }{X}^{\top }Xw-2w^{\top }{X}^{\top }y+y^{\top }y\end{array}\text{\hspace{1em}(1)}$$ 注意这里 $w^{\top }{X}^{\top }y$ 和 $y^{\top }Xw$ 两项都是 $1\times 1$ 的常数，可以合并。接下来求偏导数
  $$\begin{array}{rl}\frac{\partial (w^{\top }{X}^{\top }Xw-2w^{\top }{X}^{\top }y+y^{\top }y)}{\partial w}& =\frac{\partial (w^{\top }{X}^{\top }Xw-2w^{\top }{X}^{\top }y)}{\partial w}\\ & =\frac{\partial (w^{\top }{X}^{\top }Xw)}{\partial w}-2\frac{\partial (w^{\top }{X}^{\top }y)}{\partial w}\end{array}\text{\hspace{1em}(2)}$$ 这里用到向量内积对列向量求导公式 $\frac{\partial f({\mathbf{u}}^{\top }\mathbf{v})}{\partial \mathbf{x}}=\frac{\partial f({\mathbf{u}}^{\top })}{\partial \mathbf{x}}v+\frac{\partial f({\mathbf{v}}^{\top })}{\partial \mathbf{x}}u$，其中 $\mathbf{u},\mathbf{v}$ 是向量，根据这个变形有

  1. $\frac{\partial {\mathbf{x}}^{\top }\mathbf{x}}{\partial \mathbf{x}}=2\mathbf{x}$
  2. $\frac{\partial {\mathbf{x}}^{\top }\mathbf{Ax}}{\partial \mathbf{x}}=\frac{\partial {\mathbf{x}}^{\top }}{\partial \mathbf{x}}\mathbf{Ax}+\frac{\partial {\mathbf{x}}^{\top }{\mathbf{A}}^{\top }}{\partial \mathbf{x}}\mathbf{x}=({\mathbf{A}}^{\top }+\mathbf{A})\mathbf{x}$

  利用这两个结论展开公式2，得到
  $$\begin{array}{rl}\frac{\partial (w^{\top }{X}^{\top }Xw)}{\partial w}-2\frac{\partial (w^{\top }{X}^{\top }y)}{\partial w}& =(X^{\top }X+X^{\top }X)w-2X^{\top }y\\ & =2(X^{\top }Xw-X^{\top }y)\end{array}\text{\hspace{1em}(3)}$$ 令公式 3 为 0，有
  $$X^{\top }X\hat{w}=X^{\top }y\Rightarrow \hat{w}=(X^{\top }X{)}^{-1}{X}^{\top }y\text{\hspace{1em}(4)}$$ 这就是最小二乘法得到的解，其中 $(X^{\top }X{)}^{-1}{X}^{\top }$ 是 $X$ 的 Moore-Penrose 广义逆，也叫伪逆，记为 $X^{†}$

2.3 理解最小二乘法

2.3.1 函数图像角度

• 最小二乘法要最小化所有样本的预测值和真实值间误差的平方和，其主要适用于简单回归任务，线性模型因变量通常只有一维。一维特征空间中回归拟合一条曲线，样本误差可以看作样本特征点距离拟合函数曲线的 $y$ 轴方向上的差距；二维特征空间中回归拟合一个曲面，样本误差可以看作样本特征点距离拟合曲面的 $z$ 轴方向上的差距
  1. 常数函数模型：对于最简单的情况，试想我们要测量一根铅笔的长度，为了保证精准，使用五把不同的尺子多次测量。这种情况下模型为常数函数 $$f(x,w)=w_0$$ 拟合函数 $y=w_0$ 是一条平行于 $x$ 轴的直线，给定五个样本，各样本误差如下图绿色线所示所示
     
     按最小二乘法原理，求 ${\mathrm{argmin}}_{w_0}{\sum }_{i=1}^5(y_i-w_0{)}^2$，通过令偏导为 0 可以解出 $$\widehat{w_0}=\frac{y_1+y_2+...+y_5}{5}$$ 恰是算数平均值，因此多次测量取平均来减小误差的方式是符合最小二乘原理的
  2. 简单线性模型：当真实分布是一条直线时，对于线性模型 $$f(x,w)=w_0+w_{1}x$$ 拟合函数 $y=w_0+w_{1}x$ 是一条直线，给定五个样本，如下图所示
     
     紫色点代表真实标记 $b=[b_1,b_2,...,b_5{]}^{\top }$，绿色点代表拟合函数预测的样本标记 $p=[p_1,p_2,...,p_5{]}^{\top }$，二者间的误差 $e=[e_1,e_2,...,e_5{]}^{\top }=b-p$，我们要通过优化 $\hat{w}$ 来最小化 $||e|{|}_2$，即 $$\underset{\mathbf{w}}{\mathrm{argmin}}(p-b{)}^{\top }(p-b)=\underset{\mathbf{w}}{\mathrm{argmin}}{e}^{\top }e=\underset{\mathbf{w}}{\mathrm{argmin}}\sum _{i=1}^5{e}_i^2$$
  3. 多项式函数：对于更复杂数据分布情况，可以考虑做多项式回归，使用多项式模型 $$f(x,w)=w_0+w_{1}x+w_2{x}^2+...+w_m{x}^m$$ 在做最小二乘计算前，先做数据预处理把 $x$ 的各次方都算出来，这样就回到线性模型，另外通过选择不同的多项式函数（模型），能得到不同的拟合结果
     
     也可以扩展到多维特征空间，比如拟合得到以下二维多项式模型
     $$f(x,w)=0+x_1+2x_1^2+3x_2+4x_1{x}_2+5x_1^2{x}_2+6x_2^2+7x_1{x}_2^2+8x_1^2{x}_2^2$$
     
  4. 任意函数：从理论上讲，最小二乘法可以拟合任何函数，因为我们可以通过泰勒展开，使用多项式函数逼近将任何函数曲线。泰勒展开式如下
     $$f(x)=f(x_0)+(x-x_0)f^{\prime }(x_0)+(x-x_0{)}^2\frac{f^{\prime \prime }(x_0)}{2!}+...=\sum _{n=0}^{\infty }\frac{f^{(n)}(x_0)}{n!}(x-x_0{)}^n$$ 为了简化表示，在原点 $x_0=0$ 处展开，有
     $$f(x)=f(0)+x{f}^{\prime }(0)+x^2\frac{f^{\prime \prime }(0)}{2!}+...=\sum _{n=0}^{\infty }\frac{f^{(n)}(0)}{n!}{x}^n$$ 于是问题又变形成了多项式回归的形式

2.3.2 向量空间投影角度

• 从向量空间投影的角度看，模型 $X\hat{w}=\hat{y}$ 是对样本集 $X$ 的列向量按 $\hat{w}$ 做线性组合，得到 $X$ 列空间中的 $\hat{y}$。为了使得误差最小，这个 $\hat{y}$ 应当是真实标记向量 $y$ 在 $X$ 列空间中的投影
• 现实生活中常有这样的应用：根据一些测量值求解另一些值，比如根据卫星与几个基站的通讯延时测量其位置。这就是一个解方程的问题，假设一组观测数据可以构造一个线性方程（约束）。为了统一符号，假设每个测量样本为向量 $x\in {\mathbb{R}}^p$，参数向量 $w\in {\mathbb{R}}^p$，每个方程为 $x_i^{\top }w=b_i$，进而 $n$ 组观测数据有 $X_{n\times p}=[x_1^{\top },x_2^{\top },...,x_n^{\top }{]}^{\top }$ 和 $b=[b_1,b_2,...,b_n{]}^{\top }$，联立到得到方程组
  $$Xw=b$$为了提高准确性，一般情况下会收集很多数据，使得观测样本量（方程个数） $n$ 远大于未知数个数 $p$，除非其中含有多达 $n-p$ 个无效约束，否则这样的超定方程很难有解析解，因为我们没法对非方阵的 $X$ 求逆来得到 $w=X^{-1}b$

  • 上一节中 “常数函数模型” 这个示例中，未知数个数 $p=1$，由于每一个测量结果都有微小误差，除非 $n$ 个观测中有 $n-p=n-1$ 个都无效，否则把他们联立起来就会出现冲突导致无解，反之有解就意味着所有测量结果一模一样
  • 上一节中 “简单线性模型” 这个示例中，未知数有 $w_0,w_1$ 两个，$p=2$，只需要两个观测样本就能确定一条直线，而把这 $n=5$ 个样本联立在一起就无解（找不到一条可以通过全部五个样本点的直线），除非其中 $5-2=3$ 个样本无效（有四个样本共线）

• 这时我们可以不断地去除方程，直到剩下的方程组有解为止，但这会造成数据浪费。另一种更好的方法是找一个和已知数据 “误差最小” 的解，设这个近似解为 $\hat{w}$，我们有
  $$X\hat{w}=p$$ 这样得到的 $p$ 处于 $X$ 的列空间中，为了保证 “误差最小”，$p$ 应当是 $b$ 在 $X$ 列空间 $C(X)$ 中的投影，而误差 $e=b-p$ 在 $X$ 的左零空间 $N(X)$ 中，与 $C(X)$ 正交，如下图所示
  
  图中 $a$ 换成 $x$，灰色平面代表 $X$ 的列空间 $C(X)$，由于 $e\perp C(X)$，有
  $$\begin{array}{rl}& X^{\top }e=0\\ & \Rightarrow X^{\top }(b-p)=0\\ & \Rightarrow X^{\top }(b-X\hat{w})=0\\ & \Rightarrow X^{\top }b=X^{\top }X\hat{w}\\ & \Rightarrow \hat{w}=(X^{\top }X{)}^{-1}{X}^{\top }b\end{array}$$ 可见得到了和最小二乘法相同的解。这里我们其实推导了一遍子空间投影矩阵 $P=X(X^{\top }X{)}^{-1}{X}^{\top }$，任意向量左乘此投影矩阵都会被投影到 $C(X)$ 中，此情境中是 $Pb=p$。关于向量空间投影，可以参考 线性代数拾遗（6）—— 向量空间投影与投影矩阵
• 这里其实可以联合 2.3.1 节对比看一下
  1. 2.3.1 节的视角下，我们是把总误差分到各个样本上，对应于 $X$ 的行向量
  2. 本节视角下，我们是把总误差分到特征空间的各个维度上，对应于 $X$ 的列向量

2.3.3 概率角度

• 如果把误差值看做来自一个概率分布的采样，那么当分布是正态分布时，使用最大似然估计得到的估计参数和使用最小二乘法解出的估计参数相同
• 我们可以把每个样本的 $x_i$ 的真实标记值看做是预测标记值 $\widehat{y_i}$ 加上一个噪音 $e_i$ 得到，比如下图这个简单线性模型
  
  假设误差 $e_i\sim N(0,{\sigma }^2)$ 服从均值为 0 的高斯分布，则样本的真实标签 $$y_i=\widehat{y_i}+e_i=f(x_i,w)+e_i\sim N(f(x_i,w),{\sigma }^2)$$ 也服从高斯分布，有
  $$p(y_i|x_i,w)=\frac{1}{\sqrt{2\pi }\sigma }\text{exp}(-\frac{(y_i-f(x_i,w){)}^2}{2{\sigma }^2})$$ 由于样本独立同分布，整个样本集的似然函数为
  $$L(w)=p(y_1|x_1,w)p(y_2|x_2,w)....p(y_m|x_m,w)=\prod _{i=1}^{m}p(y_i|x_i,w)$$ 通过最大化对数似然函数的方式得到参数 $w$ 的估计值，即
  $$\begin{array}{rl}\hat{w}& =\underset{\hat{\mathbf{w}}}{\mathrm{argmax}}L(w)\\ & =\underset{\hat{\mathbf{w}}}{\mathrm{argmax}}\prod _{i=1}^{m}p(y_i|x_i,w)\\ & =\underset{\hat{\mathbf{w}}}{\mathrm{argmax}}\sum _{i=1}^m\log p(y_i|x_i,w)\\ & =\underset{\hat{\mathbf{w}}}{\mathrm{argmax}}\sum _{i=1}^m(\log \frac{1}{\sqrt{2\pi }\sigma }+\log \text{exp}(-\frac{(y_i-f(x_i,w){)}^2}{2{\sigma }^2}))\\ & =\underset{\hat{\mathbf{w}}}{\mathrm{argmax}}\sum _{i=1}^m-\frac{(y_i-f(x_i,w){)}^2}{2{\sigma }^2}\\ & =\underset{\hat{\mathbf{w}}}{\mathrm{argmin}}\sum _{i=1}^m(y_i-f(x_i,w){)}^2\end{array}$$ 得到了最小二乘法的目标，因此二者等价
• 我们认为，误差是由于随机的、无数的、独立的、多个因素造成的，因此根据中心极限定理，预测误差在大样本量的情况下确实服从正态分布，所以最小二乘法确实是合理的

2.4 $A^{\top }A$ 的可逆性

• 综上所述，面对方程组 $Ax=b$ 我们如下处理

  1. $A$ 可逆，直接 $x=A^{-1}b$ 为精确解析解
  2. $A$ 非方阵，两边同时左乘 $A^{\top }$ 变为 $A^{\top }A\hat{x}=A^{\top }b$，这时 $A^{\top }A$ 是一个方阵了，若其可逆则有 $\hat{x}=(A^{\top }A{)}^{-1}{A}^{\top }b$。注意这里 $x$ 和 $\hat{x}$ 不同，若 $x$ 有解则 $\hat{x}=x$；若 $x$ 无解则 $\hat{x}$ 是最小误差解

  可见计算时要求方阵 $A^{\top }A$ 必须可逆。本节说明以下结论

  若 $A$ 的列向量线性无关，则 $A^{\top }A$ 可逆

• 证明需要用到结论：对于任意矩阵 $A$，$A$ 列向量线性无关 $\iff$ $Ax=0$ 只有零解。这个等价可从两个角度理解

  1. $Ax$ 是对 $A$ 的列向量做线性组合，由于线性无关，要线性组合得到零向量只能是各项系数为 0
  2. 若 $A$ 是方阵，可逆矩阵的零空间只有零向量

• 证明流程

  1. 首先证明 $A^{\top }Ax=0$ 和 $Ax=0$ 同解
     $$\begin{array}{rl}& A^{\top }Ax=0\\ & \iff x^{\top }{A}^{\top }Ax=0\\ & \iff (Ax{)}^{\top }(Ax)=0\\ & \iff Ax=0\end{array}$$
  2. $A$ 列向量线性无关 $\Rightarrow Ax=0$ 只有零解 $\Rightarrow A^{\top }Ax=0$ 只有零解 $\Rightarrow A^{\top }A$ 列向量线性无关 $\Rightarrow A^{\top }A$ 可逆

• 注意这里 $A$ 放在回归任务的语境下就是上面的 $X$，其每一行是一个样本的特征向量，每一列对应一个特征，$A$ 的列向量不是线性无关说明特征设计有问题，出现了冗余特征

2.5 局限性

1. 要求 $X^{\top }X$ 必须可逆（梯度下降无此限制）

   若不可逆，可以通过对样本数据进行整理，去掉冗余特征，使其变为可逆矩阵

2. 样本特征维数很大时，$(X^{\top }X{)}^{-1}$ 计算量可能太大甚至不可行（梯度下降无此限制）

   通常超过 10000 个特征就不适合用最小二乘法了，此时可以改用梯度下降等迭代算法，或者通过主成分分析降低特征维度后再用最小二乘法

3. 模型必须是线性模型（梯度下降无此限制）

   如果拟合函数不是线性的，可以先将其转换为非线性成分的线性组合并对应地处理样本数据，再使用最小二乘法

4. 由于考虑误差的二次方，因此受离群点影响很大，使用前应做好数据预处理排除错误数据
5. 当样本量 $m$ 小于特征数 $n$ 的时候，这时拟合方程组是欠定的，常用的优化方法都无法去拟合数据。当样本量 $m$ 等于特征数 $n$ 的时候，用方程组求解就可以了。当 $m$ 大于 $n$ 时，拟合方程是超定的，也就是我们常用与最小二乘法的场景了