线性回归之最小二乘法公式推导和原理介绍

1.最小二乘法的原理

最小二乘法的主要思想是通过确定未知参数θ\thetaθ（通常是一个参数矩阵），来使得真实值和预测值的误差（也称残差）平方和最小，其计算公式为E=∑i=0nei2=∑i=1n(yi−yi^)2E=\sum_{i=0}^ne_i^2=\sum_{i=1}^n(y_i-\hat{y_i})^2E=∑i=0nei2=∑i=1n(yi−yi^)2，其中yiy_iyi是真实值，yi^\hat{y_i}yi^是对应的预测值。如下图所示（来源于维基百科，Krishnavedala的作品），就是最小二乘法的一个示例，其中红色为数据点，蓝色为最小二乘法求得的最佳解，绿色即为误差。

图1

图中有四个数据点分别为：(1, 6), (2, 5), (3, 7), (4, 10)。在线性回归中，通常我们使用均方误差来作为损失函数，均方误差可以看作是最小二乘法中的E除以m（m为样本个数），所以最小二乘法求出来的最优解就是将均方误差作为损失函数求出来的最优解。对于图中这些一维特征的样本，我们的拟合函数为hθ(x)=θ0+θ1xh_\theta(x)=\theta_0+\theta_1xhθ(x)=θ0+θ1x，所以损失函数为 J(θ0,θ1)=∑i=0m(y(i)−hθ(x(i)))2=∑i=0m(y(i)−θ0−θ1x(i))2J(\theta_0,\theta_1)=\sum_{i=0}^m(y^{(i)}-h_\theta(x^{(i)}))^2=\sum_{i=0}^m(y^{(i)}-\theta_0-\theta_1x^{(i)})^2J(θ0,θ1)=i=0∑m(y(i)−hθ(x(i)))2=i=0∑m(y(i)−θ0−θ1x(i))2（这里损失函数使用最小二乘法，并非均方误差），其中上标(i)表示第i个样本。

2.最小二乘法求解

要使损失函数最小，可以将损失函数当作多元函数来处理，采用多元函数求偏导的方法来计算函数的极小值。例如对于一维特征的最小二乘法，J(θ0,θ1)J(\theta_0,\theta_1)J(θ0,θ1)分别对θ0\theta_0θ0，θ1\theta_1θ1求偏导，令偏导等于0得：

∂J(θ0,θ1)∂θ0=−2∑i=1m(y(i)−θ0−θ1x(i))=0(2.1)\frac{\partial J(\theta_0,\theta_1)}{\partial\theta_0}=-2\sum_{i=1}^{m}(y^{(i)}-\theta_0-\theta_1x^{(i)}) = 0\tag{2.1}∂θ0∂J(θ0,θ1)=−2i=1∑m(y(i)−θ0−θ1x(i))=0(2.1)

∂J(θ0,θ1)∂θ1=−2∑i=1m(y(i)−θ0−θ1x(i))x(i)=0(2.2)\frac{\partial J(\theta_0,\theta_1)}{\partial\theta_1}=-2\sum_{i=1}^{m}(y^{(i)}-\theta_0-\theta_1x^{(i)})x^{(i)} = 0\tag{2.2}∂θ1∂J(θ0,θ1)=−2i=1∑m(y(i)−θ0−θ1x(i))x(i)=0(2.2)

联立两式，求解可得：

θ0=∑i=1m(x(i))2∑i=1my(i)−∑i=1mx(i)∑i=1mx(i)y(i)m∑i=1m(x(i))2−(∑i=1mx(i))2(2.3)\theta_0 =\frac{\sum_{i=1}^m(x^{(i)})^2\sum_{i=1}^my^{(i)}-\sum_{i=1}^mx^{(i)}\sum_{i=1}^mx^{(i)}y^{(i)}}{m\sum_{i=1}^m(x^{(i)})^2-(\sum_{i=1}^mx^{(i)})^2} \tag{2.3}θ0=m∑i=1m(x(i))2−(∑i=1mx(i))2∑i=1m(x(i))2∑i=1my(i)−∑i=1mx(i)∑i=1mx(i)y(i)(2.3)

θ1=m∑i=1mx(i)y(i)−∑i=1mx(i)∑i=1my(i)m∑i=1m(x(i))2−(∑i=1mx(i))2(2.4)\theta_1 =\frac{m\sum_{i=1}^mx^{(i)}y^{(i)}-\sum_{i=1}^mx^{(i)}\sum_{i=1}^my^{(i)}}{m\sum_{i=1}^m(x^{(i)})^2-(\sum_{i=1}^mx^{(i)})^2} \tag{2.4}θ1=m∑i=1m(x(i))2−(∑i=1mx(i))2m∑i=1mx(i)y(i)−∑i=1mx(i)∑i=1my(i)(2.4)

对于图1中的例子，代入公式(2.3)(2.3)(2.3)和(2.4)(2.4)(2.4)进行结算得，θ0=3.5,θ1=1.4,J(θ)=4.2\theta_0 = 3.5, \theta_1=1.4,J(\theta) = 4.2θ0=3.5,θ1=1.4,J(θ)=4.2。

对于n维特征的样本，同样可以采用这种方式来求解。对于特征维度(x1,x2,⋯,xn)(x_1,x_2, \cdots,x_n)(x1,x2,⋯,xn)，我们增加一个第0维x0=1x_0=1x0=1，这样增广特征向量x=(x0,x1,⋯,xn)x = (x_0,x_1,\cdots,x_n)x=(x0,x1,⋯,xn)，增广权向量为θ=(θ0,θ1,…,θn)\theta = (\theta_0, \theta_1,\dots,\theta_n)θ=(θ0,θ1,…,θn).

此时我们的拟合函数变为：

hθ(x)=∑i=0nθixi=θ0+θ1x1+⋯+θnxnh_\theta(x) = \sum_{i=0}^n\theta_ix_i =\theta_0+ \theta_1x_1 + \cdots+\theta_nx_n hθ(x)=i=0∑nθixi=θ0+θ1x1+⋯+θnxn

损失函数变为：

J(θ)=∑j=1m(hθ(x(j))−y(j))2=∑j=1m(∑i=0nθixi(j)−y(j))2J(\theta)=\sum_{j=1}^m(h_\theta(x^{(j)})-y^{(j)})^2=\sum_{j=1}^m(\sum_{i=0}^n\theta_ix_i^{(j)}-y^{(j)})^2J(θ)=j=1∑m(hθ(x(j))−y(j))2=j=1∑m(i=0∑nθixi(j)−y(j))2

损失函数J(θ)J(\theta)J(θ)分别对θi(i=0,1,…,n)\theta_i(i=0,1,\dots,n)θi(i=0,1,…,n)求偏导，得：

∂J(θ)∂θi=2∑j=1m(hθ(x(j))−y(j))x(j)=2∑j=1m(∑i=0nθixi(j)−y(j))x(j)(i=0,1,…,n)\frac{\partial J(\theta)}{\partial\theta_i} = 2\sum_{j=1}^m(h_\theta(x^{(j)})-y^{(j)})x^{(j)}=2\sum_{j=1}^m(\sum_{i=0}^n\theta_ix_i^{(j)}-y^{(j)})x^{(j)}\quad (i=0,1,\dots,n) ∂θi∂J(θ)=2j=1∑m(hθ(x(j))−y(j))x(j)=2j=1∑m(i=0∑nθixi(j)−y(j))x(j)(i=0,1,…,n)

令偏导等于0，则有：

∑j=1m(∑i=0nθixi(j)−y(j))x(j)=0(i=0,1,…,n)\sum{j=1}^m(\sum{i=0}^n\theta_ix_i^{(j)}-y^{(j)})x^{(j)}=0\quad (i=0,1,\dots,n) ∑j=1m(∑i=0nθixi(j)−y(j))x(j)=0(i=0,1,…,n)

这样最终得到的结果就是一个线性方程组，未知数的个数为n+1，方程的个数也为n+1，这样就可以通过高斯消元法解出θi(i=0,1,…,n)\theta_i(i=0,1,\dots,n)θi(i=0,1,…,n)，具体可参见：详解最小二乘法原理和代码。

对于线性回归问题，我们可以依据拟合函数的形式进行特征空间变换，即广义线性回归。例如，hθ(x)=θ0+θ1x1+θ2x22h_\theta(x)=\theta_0+\theta_1x_1+\theta_2x_2^2hθ(x)=θ0+θ1x1+θ2x22，我们可以令x2:=x22x_2:=x_2^2x2:=x22，这里:=:=:=表示赋值，即将右边的值赋给左边。这样又变成了我们所熟悉的拟合函数形式。

对于非线性回归问题，最小二乘法的思想同样适用，只不过函数形式有所变化。例如，对于拟合函数hθ(x)=θ0+θ1x+θ2lnxh_\theta(x)=\theta_0+\theta_1x+\theta_2l nxhθ(x)=θ0+θ1x+θ2lnx，此时J(θ)=∑j=1m(hθ(x(j))−y(j))2J(\theta)=\sum_{j=1}^m(h_\theta(x^{(j)})-y^{(j)})^2J(θ)=∑j=1m(hθ(x(j))−y(j))2，求偏导的结果为：

∂J(θ)∂θi=2∑j=1m(hθ(x(j))−y(j))∂hθ(x)θi(i=0,1,2);其中∂hθ(x)θ0=1,∂hθ(x)θ2=x,∂hθ(x)θ2=lnx\frac{\partial J(\theta)}{\partial\theta_i}=2\sum_{j=1}^{m}(h_\theta(x^{(j)})-y^{(j)})\frac{\partial h_\theta(x)}{\theta_i}\quad (i=0,1,2);其中\frac{\partial h_\theta(x)}{\theta_0} = 1, \frac{\partial h_\theta(x)}{\theta_2} = x, \frac{\partial h_\theta(x)}{\theta_2} = lnx∂θi∂J(θ)=2j=1∑m(hθ(x(j))−y(j))θi∂hθ(x)(i=0,1,2);其中θ0∂hθ(x)=1,θ2∂hθ(x)=x,θ2∂hθ(x)=lnx

同样可以构造线性方程组，用高斯消元法求解。

3.矩阵求解最小二乘法

对于函数hθ(x)=θ0+θ1x1+⋯+θnxnh_\theta(x)=\theta_0+\theta_1x_1+\dots+\theta_nx_nhθ(x)=θ0+θ1x1+⋯+θnxn，我们将其用矩阵表示为：

Xθ=Y(3.1)X\theta = Y \qquad \tag{3.1}Xθ=Y(3.1)

其中，

X={(x(1))T(x(2))T⋮(x(m))T},Y={y(1)y(2)⋮y(m))},x(j)={x0(j))x1(j)⋮xn(j)},θ={θ0θ1⋮θn}X = \left\{\begin{matrix} (x^{(1)})^T \\ (x^{(2)})^T \\ \vdots \\(x^{(m)})^T \end{matrix} \right\} , Y = \left\{\begin{matrix} y^{(1)} \\ y^{(2)} \\ \vdots \\y^{(m)}) \end{matrix} \right\},x^{(j)}=\left\{\begin{matrix} x_0^{(j)}) \\ x_1^{(j)} \\ \vdots \\ x_n^{(j)} \end{matrix} \right\}, \theta = \left\{\begin{matrix} \theta_0 \\ \theta_1 \\ \vdots \\ \theta_n \end{matrix} \right\}X=⎩⎨⎧(x(1))T(x(2))T⋮(x(m))T⎭⎬⎫,Y=⎩⎨⎧y(1)y(2)⋮y(m))⎭⎬⎫,x(j)=⎩⎨⎧x0(j))x1(j)⋮xn(j)⎭⎬⎫,θ=⎩⎨⎧θ0θ1⋮θn⎭⎬⎫

m表示样本个数，n为特征维度，x0(i)=1(i=0,1,…,m)x_0^{(i)}=1\quad(i = 0,1,\dots,m)x0(i)=1(i=0,1,…,m)，即XXX的第一列全为1，xi(j)x_i^{(j)}xi(j)表示第j个样本的第i个特征，XXX为增广样本矩阵((1+n)*m维)，YYY为真实值组成的列向量。

损失函数表示为：

J(θ)=∑j=1m(hθ(x(j))−y(j))2=(Xθ−Y)T(Xθ−Y)(3.2)J(\theta)=\sum_{j=1}^m(h_\theta(x^{(j)})-y^{(j)})^2=(X\theta−Y)^T(X\theta−Y) \tag{3.2}J(θ)=j=1∑m(hθ(x(j))−y(j))2=(Xθ−Y)T(Xθ−Y)(3.2)

根据最小二乘法，利用矩阵求导得：（具体推导参见线性回归矩阵推导和线性回归相关向量求导）

∂J(θ)∂θ=2XT(Xθ−Y)\frac{\partial J(\theta)}{\partial\theta}=2X^T(X\theta-Y)∂θ∂J(θ)=2XT(Xθ−Y)

令求导结果等于0矩阵，可得：

XTXθ=XTY⇒θ=(XTX)−1XTY(3.3)X^TX\theta = X^TY\quad\Rightarrow \quad \theta = (X^TX)^{-1}X^TY \tag{3.3}XTXθ=XTY⇒θ=(XTX)−1XTY(3.3)

对于图1中的例子，利用公式(3.3)(3.3)(3.3)计算得：θ={3.51.4}\theta = \left\{\begin{matrix} 3.5 \\1.4\end{matrix} \right\}θ={3.51.4}

4.总结

最小二乘法可以直接求解参数矩阵，在计算时可以直接套入公式。但是仍有一定的局限性，主要体现在：

1.XTXX^TXXTX的逆矩阵可能不存在，这个在Matlab中，可以通过求伪逆来进行计算。

2.对于(3.1)(3.1)(3.1)式，可以将其看成一个线性方程组（假设各方程线性无关）。如果样本个数m小于特征维数n，那么此方程组有无穷多个解。如果m = n，有唯一解。如果m大于n，无解（即存在矛盾解）。最小二乘法一般是在m大于n的时候使用，此时求出来的解是最优近似解。

3.最小二乘法的时间复杂度为O(n3)O(n^3)O(n3)，当n特别大的时候（一般大于10000），求逆矩阵的过程非常复杂。此时采用最小二乘法，会非常耗时。

参考链接：

最小二乘法小结

半小时学习最小二乘法