1. 背景

文章的背景取自An Introduction to Gradient Descent and Linear Regression，本文想在该文章的基础上，完整地描述线性回归算法。部分数据和图片取自该文章。没有太多时间抠细节，所以难免有什么缺漏错误之处，望指正。

线性回归的目标很简单，就是用一条线，来拟合这些点，并且使得点集与拟合函数间的误差最小。如果这个函数曲线是一条直线，那就被称为线性回归，如果曲线是一条二次曲线，就被称为二次回归。数据来自于GradientDescentExample中的data.csv文件，共100个数据点，如下图所示：

我们的目标是用一条直线来拟合这些点。既然是二维，那么 y=b+mx <script type="math/tex" id="MathJax-Element-1">y = b + mx</script>这个公式相信对于中国学生都很熟悉。其中 b <script type="math/tex" id="MathJax-Element-2">b</script>是直线在y轴的截距（y-intercept），m<script type="math/tex" id="MathJax-Element-3">m</script>是直线的斜率（slope）。寻找最佳拟合直线的过程，其实就是寻找最佳的 b <script type="math/tex" id="MathJax-Element-4">b</script>和m<script type="math/tex" id="MathJax-Element-5">m</script>的过程。为了寻找最佳的拟合直线，这里首先要定义，什么样的直线才是最佳的直线。我们定义误差（cost function）：

<script type="math/tex; mode=display" id="MathJax-Element-6">误差函数\ Error_{(b, m)}=\frac{1}{N}\sum_{1}^{N}((b+mx_i)-y_i)^{2}</script>

计算损失函数的python代码如下：

# y = b + mx
def compute_error_for_line_given_points(b, m, points):
    totalError = sum((((b + m * point[0]) - point[1]) ** 2 for point in points))
    return totalError / float(len(points))

现在问题被转化为，寻找参数 b <script type="math/tex" id="MathJax-Element-7">b</script>和m<script type="math/tex" id="MathJax-Element-8">m</script>，使得误差函数 Error(b,m) <script type="math/tex" id="MathJax-Element-9">Error_{(b, m)}</script>有最小值。在这里， xi <script type="math/tex" id="MathJax-Element-10">x_i</script>和 yi <script type="math/tex" id="MathJax-Element-11">y_i</script>都被视为已知值。从下图看，最小二乘法所做的是通过数学推导直接计算得到最低点；而梯度下降法所做的是从图中的任意一点开始，逐步找到图的最低点。

2. 多元线性回归模型

从机器学习的角度来说，以上的数据只有一个feature，所以用一元线性回归模型即可。这里我们将一元线性模型的结论一般化，即推广到多元线性回归模型。这部分内部参考了机器学习中的数学(1)-回归(regression)、梯度下降(gradient descent)。假设有 x1 <script type="math/tex" id="MathJax-Element-12">x_1</script>， x2 <script type="math/tex" id="MathJax-Element-13">x_2</script>， ... <script type="math/tex" id="MathJax-Element-14">...</script>, xn <script type="math/tex" id="MathJax-Element-15">x_n</script>共 n <script type="math/tex" id="MathJax-Element-16">n</script>个feature，θ<script type="math/tex" id="MathJax-Element-17">\theta</script>为 x <script type="math/tex" id="MathJax-Element-18">x</script>的系数，则

<script type="math/tex; mode=display" id="MathJax-Element-19">拟合函数\ h_\theta(x)=\theta_0+\theta_1x_1+...+\theta_nx_n=\theta^Tx，其中x_0=1</script>

<script type="math/tex; mode=display" id="MathJax-Element-20">误差函数\ J(\theta)=\frac{1}{2}\sum_{i=1}^{m}(h_\theta(x^{(i)})-y^{(i)})^2，m代表有m组样本</script>

更一般地，我们可以得到广义线性回归。 ϕ(x) <script type="math/tex" id="MathJax-Element-21">\phi (x)</script>可以换成不同的函数，从而得到的拟合函数就不一定是一条直线了。

<script type="math/tex; mode=display" id="MathJax-Element-22">广义线性函数\ h_\theta(x)=\theta^Tx=\theta_0+\sum_{i=1}^{n}\theta_i\phi _i(x_i)</script>

2.1 误差函数的进一步思考

这里有一个有意思的东西，就是误差函数为什么要写成这样的形式。首先是误差函数最前面的系数 12 <script type="math/tex" id="MathJax-Element-23">\frac{1}{2}</script>，这个参数其实对结果并没有什么影响，这里之所以取 12 <script type="math/tex" id="MathJax-Element-24">\frac{1}{2}</script>，是为了抵消求偏导过程中得到的 2 <script type="math/tex" id="MathJax-Element-25">2</script>。可以实验，把Error(b,m)<script type="math/tex" id="MathJax-Element-26">Error_{(b,m)}</script>最前面的 1N <script type="math/tex" id="MathJax-Element-27">\frac{1}{N}</script>修改或者删除并不会改变最终的拟合结果。那么为什么要使用平方误差呢？考虑以下公式：

<script type="math/tex; mode=display" id="MathJax-Element-28">y^{(i)}=\theta^Tx^{(i)}+\varepsilon ^{(i)}</script>

假定误差 ε(i)(1⩽i⩽m) <script type="math/tex" id="MathJax-Element-29">\varepsilon ^{(i)}(1\leqslant i\leqslant m)</script>是独立同分布的，由中心极限定理可得， ε(i) <script type="math/tex" id="MathJax-Element-30">\varepsilon ^{(i)}</script>服从均值为 0 <script type="math/tex" id="MathJax-Element-31">0</script>，方差为σ2<script type="math/tex" id="MathJax-Element-32">\sigma ^2</script>的正态分布（均值若不为0，可以归约到 θ0上 <script type="math/tex" id="MathJax-Element-33">\theta_0上</script>）。进一步的推导来自从@邹博_机器学习的机器学习课件。

所以求 maxL(θ) <script type="math/tex" id="MathJax-Element-34">max L(\theta)</script>的过程，就变成了求 minJ(θ) <script type="math/tex" id="MathJax-Element-35">minJ(\theta)</script>的过程，从理论上解释了误差函数 J(θ) <script type="math/tex" id="MathJax-Element-36">J(\theta)</script>的由来。

3 最小二乘法求误差函数最优解

最小二乘法（normal equation）相信大家都很熟悉，这里简单进行解释并提供python实现。首先，我们进一步把 J(θ) <script type="math/tex" id="MathJax-Element-37">J(\theta)</script>写成矩阵的形式。 X <script type="math/tex" id="MathJax-Element-38">X</script>为m<script type="math/tex" id="MathJax-Element-39">m</script>行 n <script type="math/tex" id="MathJax-Element-40">n</script>列的矩阵（代表m<script type="math/tex" id="MathJax-Element-41">m</script>个样本，每个样本有 n <script type="math/tex" id="MathJax-Element-42">n</script>个feature），θ<script type="math/tex" id="MathJax-Element-43">\theta</script>和 Y <script type="math/tex" id="MathJax-Element-44">Y</script>为m<script type="math/tex" id="MathJax-Element-45">m</script>行 1 <script type="math/tex" id="MathJax-Element-46">1</script>列的矩阵。所以

<script type="math/tex; mode=display" id="MathJax-Element-47">\ J(\theta)=\frac{1}{2}\sum_{i=1}^{m}(h_\theta(x^{(i)})-y^{(i)})^2=\frac{1}{2}(X\theta-Y)^T(X\theta-Y)</script>

所以 θ <script type="math/tex" id="MathJax-Element-48">\theta</script>的最优解为： θ=(XTX)−1XTY <script type="math/tex" id="MathJax-Element-49">\theta=(X^TX)^{-1}X^TY</script>。

当然这里可能遇到一些问题，比如 X <script type="math/tex" id="MathJax-Element-50">X</script>必须可逆，比如求逆运算时间开销较大。具体解决方案待补充。

3.1 python实现最小二乘法

这里的代码仅仅针对背景里的这个问题。部分参考了回归方法及其python实现。

# 通过最小二乘法直接得到最优系数，返回计算出来的系数b, m
def least_square_regress(points):
    x_mat = np.mat(np.array([np.ones([len(points)]), points[:, 0]]).T)  # 转为100行2列的矩阵，2列其实只有一个feature，其中x0恒为1
    y_mat = points[:, 1].reshape(len(points), 1)  # 转为100行1列的矩阵
    xT_x = x_mat.T * x_mat
    if np.linalg.det(xT_x) == 0.0:
        print('this matrix is singular,cannot inverse')  # 奇异矩阵，不存在逆矩阵
        return
    coefficient_mat = xT_x.I * (x_mat.T * y_mat)
    return coefficient_mat[0, 0], coefficient_mat[1, 0] # 即系数b和m

程序执行结果如下：
b = 7.99102098227, m = 1.32243102276, error = 110.257383466, 相关系数 = 0.773728499888

拟合结果如下图：

4. 梯度下降法求误差函数最优解

有了最小二乘法以后，我们已经可以对数据点进行拟合。但由于最小二乘法需要计算X<script type="math/tex" id="MathJax-Element-51">X</script>的逆矩阵，计算量很大，因此特征个数多时计算会很慢，只适用于特征个数小于100000时使用；当特征数量大于100000时适合使用梯度下降法。最小二乘法与梯度下降法的区别见最小二乘法和梯度下降法有哪些区别？。

4.1. 梯度

首先，我们简单回顾一下微积分中梯度的概念。这里参考了方向导数与梯度，具体的证明请务必看一下这份材料，很短很简单的。

讨论函数 z=f(x,y) <script type="math/tex" id="MathJax-Element-52">z = f(x, y)</script>在某一点 P <script type="math/tex" id="MathJax-Element-53">P</script>沿某一方向的变化率问题。设函数z=f(x,y)<script type="math/tex" id="MathJax-Element-54">z = f(x, y)</script>在点 P(x,y) <script type="math/tex" id="MathJax-Element-55">P(x, y)</script>的某一邻域 U(P) <script type="math/tex" id="MathJax-Element-56">U(P)</script>内有定义，自点 P <script type="math/tex" id="MathJax-Element-57">P</script>引射线l<script type="math/tex" id="MathJax-Element-58">l</script>到点 P′(x+Δx,y+Δy) <script type="math/tex" id="MathJax-Element-59">P'(x+\Delta x, y+\Delta y)</script>且 P′∈U(P) <script type="math/tex" id="MathJax-Element-60">P'\in U(P)</script>，如下图所示。

定义函数 z=f(x,y) <script type="math/tex" id="MathJax-Element-61">z = f(x, y)</script>在点 P <script type="math/tex" id="MathJax-Element-62">P</script>沿方向l<script type="math/tex" id="MathJax-Element-63">l</script>的方向导数为：

<script type="math/tex; mode=display" id="MathJax-Element-64">\frac{\partial f}{\partial l}=\lim_{\rho \to 0}\frac{f(x+\Delta x, y+\Delta y)-f(x,y)}{\rho}，其中\rho = \sqrt{(\Delta x)^2+(\Delta y)^2}</script>

方向导数可以理解为，函数 z=f(x,y) <script type="math/tex" id="MathJax-Element-65">z = f(x, y)</script>沿某个方向变化的速率。可以类比一下函数 y=kx+b <script type="math/tex" id="MathJax-Element-66">y=kx+b</script>的斜率 k=dydx <script type="math/tex" id="MathJax-Element-67">k=\frac{{\rm d}y}{{\rm d}x}</script>。斜率越大，函数 y <script type="math/tex" id="MathJax-Element-68">y</script>增长得越快。那么现在问题来了，函数z=f(x,y)<script type="math/tex" id="MathJax-Element-69">z = f(x, y)</script>在点 P <script type="math/tex" id="MathJax-Element-70">P</script>沿哪个方向增加的速度最快？而这个方向就是梯度的方向

<script type="math/tex; mode=display" id="MathJax-Element-71">gradf(x, y)=\frac{\partial f}{\partial x}\overrightarrow{i}+\frac{\partial f}{\partial y}\overrightarrow{j}</script>

从几何角度来理解，函数 z=f(x,y) <script type="math/tex" id="MathJax-Element-72">z = f(x, y)</script>表示一个曲面，曲面被平面 z=c <script type="math/tex" id="MathJax-Element-73">z=c</script>截得的曲线在 xoy <script type="math/tex" id="MathJax-Element-74">xoy</script>平面上投影如下图，这个投影也就是我们所谓的等高线。

函数 z=f(x,y) <script type="math/tex" id="MathJax-Element-75">z = f(x, y)</script>在点 P(x,y) <script type="math/tex" id="MathJax-Element-76">P(x, y)</script>处的梯度方向与点 P <script type="math/tex" id="MathJax-Element-77">P</script>的等高线f(x,y)=c<script type="math/tex" id="MathJax-Element-78">f(x, y)=c</script>在这点的法向量的方向相同，且从数值较低的等高线指向数值较高的等高线。

4.2 梯度方向计算

理解了梯度的概念之后，我们重新回到1. 背景中提到的例子。1. 背景提到，梯度下降法所做的是从图中的任意一点开始，逐步找到图的最低点。那么现在问题来了，从任意一点开始， b <script type="math/tex" id="MathJax-Element-79">b</script>和m<script type="math/tex" id="MathJax-Element-80">m</script>可以往任意方向”走”，如何可以保证我们走的方向一定是使误差函数 Error(b,m) <script type="math/tex" id="MathJax-Element-81">Error_{(b, m)}</script>减小且减小最快的方向呢？回忆4.1. 梯度中提到的结论，梯度的方向是函数上升最快的方向，那么函数下降最快的方向，也就是梯度的反方向。有了这个思路，我们首先计算梯度方向，

<script type="math/tex; mode=display" id="MathJax-Element-82">\frac{\partial Error_{(b, m)}}{\partial m}=\sum_{i=1}^{N}x_i((b+mx_i)-y_i)</script> <script type="math/tex; mode=display" id="MathJax-Element-83">\frac{\partial Error_{(b, m)}}{\partial b}=\sum_{i=1}^{N}((b+mx_i)-y_i)，x_0恒为1</script>

有了这两个结果，我们就可以开始使用梯度下降法来寻找误差函数 Error(b,m) <script type="math/tex" id="MathJax-Element-84">Error_{(b, m)}</script>的最低点。我们从任意的点 (b,m) <script type="math/tex" id="MathJax-Element-85">(b, m)</script>开始，逐步地沿梯度的负方向改变 b <script type="math/tex" id="MathJax-Element-86">b</script>和m<script type="math/tex" id="MathJax-Element-87">m</script>的值。每一次改变， Error(b,m) <script type="math/tex" id="MathJax-Element-88">Error_{(b,m)}</script>都会得到更小的值，反复进行该操作，逐步逼近 Error(b,m) <script type="math/tex" id="MathJax-Element-89">Error_{(b,m)}</script>的最低点。

回到更一般的情况，对于每一个向量 θ <script type="math/tex" id="MathJax-Element-90">\theta</script>的每一维分量 θi <script type="math/tex" id="MathJax-Element-91">\theta_i</script>，我们都可以求出梯度的方向，也就是错误函数 J(θ) <script type="math/tex" id="MathJax-Element-92">J(\theta)</script>下降最快的方向：

<script type="math/tex; mode=display" id="MathJax-Element-93">\frac{\partial}{\partial \theta_j}J(\theta)=\frac{\partial}{\partial \theta_j}\frac{1}{2}\sum_{i=1}^{m}(h_\theta(x^{(i)})-y^{(i)})^2=\sum_{i=1}^{m}(h_\theta(x^{(i)})-y^{(i)})x_j^{(i)}</script>

4.3 批量梯度下降法

从上面的公式中，我们进一步得到特征的参数 θj <script type="math/tex" id="MathJax-Element-94">\theta_j</script>的迭代式。因为这个迭代式需要把m个样本全部带入计算，所以我们称之为批量梯度下降

<script type="math/tex; mode=display" id="MathJax-Element-95">\theta'_j=\theta_j-\alpha \frac{\partial J(\theta)}{\partial \theta_j}=\theta_j-\alpha \sum_{i=1}^{m}(h_\theta(x^{(i)})-y^{(i)})x_j^{(i)}</script>

针对此例，梯度下降法一次迭代过程的python代码如下：

def step_gradient(b_current, m_current, points, learningRate):
    b_gradient = 0
    m_gradient = 0
    N = float(len(points))
    for i in range(0, len(points)):
        x = points[i, 0]
        y = points[i, 1]
        m_gradient += (2 / N) * x * ((b_current + m_current * x) - y)
        b_gradient += (2 / N) * ((b_current + m_current * x) - y)
    new_b = b_current - (learningRate * b_gradient)  # 沿梯度负方向
    new_m = m_current - (learningRate * m_gradient)  # 沿梯度负方向
    return [new_b, new_m]

其中learningRate是学习速率，它决定了逼近最低点的速率。可以想到的是，如果learningRate太大，则可能导致我们不断地最低点附近来回震荡； 而learningRate太小，则会导致逼近的速度太慢。An Introduction to Gradient Descent and Linear Regression提供了完整的实现代码GradientDescentExample。

这里多插入一句，如何在python中生成GIF动图。配置的过程参考了使用Matplotlib和Imagemagick实现算法可视化与GIF导出。需要安装ImageMagick，使用到的python库是Wand: a ctypes-based simple ImageMagick binding for Python。然后修改C:\Python27\Lib\site-packages\matplotlib__init__.py文件，在

# this is the instance used by the matplotlib classes
rcParams = rc_params()

后面加上：

# fix a bug by ZZR
rcParams['animation.convert_path'] = 'C:\Program Files\ImageMagick-6.9.2-Q16\convert.exe'

即可在python中调用ImageMagick。如何画动图参见Matplotlib动画指南，不再赘述。

learningRate=0.0001，迭代100轮的结果如下图：

After {100} iterations b = 0.0350749705923, m = 1.47880271753, error = 112.647056643, 相关系数 = 0.773728499888
After {1000} iterations b = 0.0889365199374, m = 1.47774408519, error = 112.614810116, 相关系数 = 0.773728499888
After {1w} iterations b = 0.607898599705, m = 1.46754404363, error = 112.315334271, 相关系数 = 0.773728499888
After {10w} iterations b = 4.24798444022, m = 1.39599926553, error = 110.786319297, 相关系数 = 0.773728499888

4.4 随机梯度下降法

批量梯度下降法每次迭代都要用到训练集的所有数据，计算量很大，针对这种不足，引入了随机梯度下降法。随机梯度下降法每次迭代只使用单个样本，迭代公式如下：

<script type="math/tex; mode=display" id="MathJax-Element-96">\theta'_j=\theta_j-\alpha (h_\theta(x^{(i)})-y^{(i)})x_j^{(i)}</script>

可以看出，随机梯度下降法是减小单个样本的错误函数，每次迭代不一定都是向着全局最优方向，但大方向是朝着全局最优的。

这里还有一些重要的细节没有提及，比如如何确实learningRate，如果判断何时递归可以结束等等。

---

参考文献

1. An Introduction to Gradient Descent and Linear Regression
2. 方向导数与梯度
3. 最小二乘法和梯度下降法有哪些区别？
4. GradientDescentExample
5. 机器学习中的数学(1)-回归(regression)、梯度下降(gradient descent)
6. @邹博_机器学习
7. 回归方法及其python实现
8. 使用Matplotlib和Imagemagick实现算法可视化与GIF导出
9. Wand: a ctypes-based simple ImageMagick binding for Python
10. Matplotlib动画指南