02. 一元线性回归的求解

来源：原始 Jupyter Notebook。内容已转换为网页阅读格式；下载原文件。

import numpy as np
import matplotlib.pyplot as plt

plt.xlim(30, 60)
plt.ylim(100, 600)

X = np.array([39.93, 42.05, 43.18, 44.68, 49.87, 53.57])
Y = np.array([199,   290,   298,   310,   399,   420])
plt.scatter(X, Y)

x = np.linspace(30, 60, 100)
y = 15*x-360
plt.plot(x, y)

plt.show()

运行结果

Notebook 运行结果

求解一元线性回归，实际上就是找一条直线，使得它和散点图中的所有点距离最近。直线和散点的距离一般通过损失函数来度量，损失函数（Loss Function）有时也称为代价函数（Cost Function），常用的损失函数有两种，一种是 绝对损失（Absolute Loss），另一种是 平方损失（Squared Loss）：

绝对损失函数：

loss = |y - \hat{y}|

平方损失函数：

loss = (y - \hat{y})^2

其中， $\hat{y}$ 表示预测值，由于它是一元线性的，我们可以记为： $\hat{y} = ax + b$ ，这条直线也叫 回归线。这条回归线就是我们要学习的对象，可以把它叫做模型（model），它反映了数据的某种潜在规律，因此也称为假设（hypothesis）。

在一元线性回归的实际求解中，因为平方损失处处可导，可以方便后面的计算，常常使用平方损失函数：

\begin{equation} \begin{aligned} loss &= (y - \hat{y})^2 \\ &= (y - (ax + b))^2 \\ &= \sum_{i=1}^{n}(y_i - (ax_i + b))^2\\ &= (y_1 - (ax_1 + b))^2 + (y_2 - (ax_2 + b))^2 + ... + (y_n - (ax_n + b))^2 \end{aligned} \end{equation}

将其完全展开可得：

loss = (y_1^2 + y_2^2 + ... + y_n^2) + a^2(x_1^2 + x_2^2 + ... + x_n^2) + nb^2 + 2ab(x_1 + x_2 + ... + x_n) - 2a(x_1y_1 + x_2y_2 + ... + x_ny_n) -2b(y_1 + y_2 + ... + y_n)

令 $\bar{x} = \frac{x_1 + x_2 + ... + x_n}{n}$ ，表示 x 的平均数，则：

loss = n\bar{y^2} + a^2n\bar{x^2} + nb^2 + 2abn\bar{x} - 2an\bar{xy} - 2bn\bar{y}

很显然，这是一个关于 a、b 变量的二元二次函数，我们记为：

F(a, b) = n\bar{y^2} + a^2n\bar{x^2} + nb^2 + 2abn\bar{x} - 2an\bar{xy} - 2bn\bar{y}

那么问题就转变成了求函数 F(a, b) 的最小值。为了求解这个问题，我们先看看一元函数的极值问题。

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(-2, 2, 100)
y = x*x*x*x - 2*x*x + 8
plt.plot(x, y)

plt.show()

运行结果

Notebook 运行结果

上面画出了一元函数的图像：

y = x^4 - 2x^2 + 8

从图中可以很直观的看出，当 x = 0 时，函数有极大值 8，当 x = 1 或 x = -1 时，函数有极小值 7。这里引申出两个数学概念：极值和最值。在数学分析中，某个函数在指定范围内的最大值和最小值被统称为极值，一般通过导数来求极值和最值：

y' = 4x^3 -4x

令 $y' = 0$ ，可以求得：

x^3 - x = 0 \to x = 0, x = \pm 1

根据极值点将函数的定义域划分为几个区间，并计算各个区间的单调性：

x	( $-\infty$ , -1)	-1	(-1, 0)	0	(0,1)	1	(1, $+\infty$ )
y'	-	0	+	0	-	0	+
y	$\searrow$	7	$\nearrow$	8	$\searrow$	7	$\nearrow$

检查极值点两侧导数的符号，如果是左正右负，则为极大值点，如果左负右正，则为极小值点。显然，由于函数在 ( $-\infty$ , -1) 区间内单调递减，在 (1, $+\infty$ ) 区间内单调递增，所以比较所有的极小值，最小的那个就是函数在整个定义域内的最小值。

这里介绍的是一元函数的极值求解方法，再回到我们刚刚的损失函数：

F(a, b) = n\bar{y^2} + a^2n\bar{x^2} + nb^2 + 2abn\bar{x} - 2an\bar{xy} - 2bn\bar{y}

这是一个关于 a、b 变量的二元函数，类似的，二元函数的极值也可以通过导数来求解，不过这里要用到偏导的概念，也就是只对其中一个变量求导，其余变量当作常数看待，那么二元函数就有两个偏导，并令其等于 0：

\left\{ \begin{aligned} \frac{\partial}{\partial a} F(a, b) &=& 2n\bar{x^2}a + 2bn\bar{x} - 2n\bar{xy} &=& 0 \\ \frac{\partial}{\partial b} F(a, b) &=& 2nb + 2an\bar{x} - 2n\bar{y} &=& 0 \\ \end{aligned} \right.

结合两个方程可以求得：

\left\{ \begin{aligned} a &=& \frac{\bar{xy} - \bar{x}\bar{y}}{\bar{x^2} - (\bar{x})^2} \\ b &=& \bar{y} - a\bar{x} \end{aligned} \right.

通过这个公式，就可以求出 a，b 的解，这种严格按照公式推导计算出来的解叫做 闭式解（closed-form solution）或 解析解（analytical solution）。

import numpy as np

X = np.array([39.93, 42.05, 43.18, 44.68, 49.87, 53.57])
Y = np.array([199,   290,   298,   310,   399,   420])

bar_x = np.sum(X) / X.size
bar_y = np.sum(Y) / Y.size
print("bar_x = {0}".format(bar_x))
print("bar_y = {0}".format(bar_y))

bar_xx = np.sum(np.multiply(X, X)) / X.size
bar_yy = np.sum(np.multiply(Y, Y)) / X.size
bar_xy = np.sum(np.multiply(X, Y)) / X.size
print("bar_xx = {0}".format(bar_xx))
print("bar_yy = {0}".format(bar_yy))
print("bar_xy = {0}".format(bar_xy))

a = (bar_xy - bar_x * bar_y) / (bar_xx - bar_x * bar_x)
b = bar_y - a * bar_x
print("a = {0}".format(a))
print("b = {0}".format(b))

运行结果

bar_x = 45.546666666666674
bar_y = 319.3333333333333
bar_xx = 2096.697333333333
bar_yy = 107367.66666666667
bar_xy = 14876.089999999997
a = 14.93439993913547
b = -360.878802561157

所以，求得一元线性回归的解为：

y = ax + b = 14.93x - 360.88

这种通过求解平方损失函数的最小值来求解线性回归的方法，我们称之为 最小二乘法（Least Squares，简称 LS），这是机器学习中最基础的一个算法。

通常我们会把损失函数写成下面的形式：

J(a,b) = \frac{1}{2m} \sum_{i=1}^{m}(h(x^{(i)}) - y^{(i)})^2

这个式子也叫做 均方误差（mean-square error, MSE），在统计学中， $h(x^{(i)}) - y^{(i)}$ 被称为残差，所以这个式子又叫 残差平方和（Sum of Squares due to Error，SSE）。其实它就是我们常用的欧几里得距离（欧式距离，Euclidean Distance），也就是残差的 $L_2$ 范数。让损失函数的值最小化，求参数 a b 的值可以写成：

(a^*, b^*) = arg min_{a,b} J(a,b)

求解 a 和 b 使损失函数最小化的过程，称为线性回归模型的最小二乘参数估计（parameter estimation）。

实际上，上面的损失函数就是三维空间里的一个曲面（图片来自 Andrew Ng 的在线课程）：

损失函数的最小值就是该曲面的最低点。