Python机器学习算法入门之梯度下降法实现线性回归

專 欄

❈

ZZR ，Python中文社区专栏作者，OpenStack工程师，曾经的NLP研究者。主要兴趣方向：OpenStack、Python爬虫、Python数据分析。

Blog：http://skydream.me/

CSDN：http://blog.csdn.net/titan0427/article/details/50365480

❈

—

—

1. 背景

    文章的背景取自

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

    线性回归的目标很简单，就是用一条线，来拟合这些点，并且使得点集与拟合函数间的误差最小。如果这个函数曲线是一条直线，那就被称为线性回归，如果曲线是一条二次曲线，就被称为二次回归。数据来自于

GradientDescentExample 中的data.csv文件，共100个数据点，如下图所示：

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

function）：

误

差

函

数

E r r o

r

( b , m )

=

1 N

∑

1

N

( ( b

m

x i

) −

y i

)

2

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

1. # y = b + mx

2. def compute_error_for_line_given_points(b, m, points):

3. totalError = sum((((b + m * point[0]) - point[1]) ** 2 for point in points))

4. return totalError / float(len(points))

    现在问题被转化为，寻找参数b和m，使得误差函数Error(b,m)有最小值。在这里，xi和yi都被视为已知值。从下图看，最小二乘法所做的是通过数学推导直接计算得到最低点；而梯度下降法所做的是从图中的任意一点开始，逐步找到图的最低点。
   

2. 多元线性回归模型

    从机器学习的角度来说，以上的数据只有一个feature，所以用一元线性回归模型即可。这里我们将一元线性模型的结论一般化，即推广到多元线性回归模型。这部分内部参考了

机器学习中的数学(1)-回归(regression)、梯度下降(gradient descent) 。假设有x1，x2，…, xn共n个feature，θ为x的系数，则

拟

合

函

数

h θ

θ 0

θ 1

x 1

. . .

θ n

x n

=

θ T

x

，

其

中

x 0

= 1

误

差

函

数

1 2

∑

1

m

(

h θ

(

x

( i )

) −

y

( i )

) 2

，

m

代

表

有

m

组

样

本

    更一般地，我们可以得到广义线性回归。ϕ(x)可以换成不同的函数，从而得到的拟合函数就不一定是一条直线了。

广

义

线

性

函

数

h θ

θ T

θ 0

∑

1

n

θ i

ϕ i

(

x i

)

2.1 误差函数的进一步思考

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

y

( i )

=

θ T

x

( i )

ε

( i )

    假定误差ε(i)(1⩽i⩽m)是独立同分布的，由中心极限定理可得，ε(i)服从均值为0，方差为σ2的正态分布（均值若不为0，可以归约到θ0上）。进一步的推导来自从

@邹博_机器学习 的机器学习课件。

    所以求maxL(θ)的过程，就变成了求minJ(θ)的过程，从理论上解释了误差函数J(θ)的由来。

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

    最小二乘法（normal equation）相信大家都很熟悉，这里简单进行解释并提供python实现。首先，我们进一步把J(θ)写成矩阵的形式。X为m行n列的矩阵（代表m个样本，每个样本有n个feature），θ和Y为m行1列的矩阵。所以   

1 2

∑

1

m

(

h θ

(

x

( i )

) −

y

( i )

) 2

=

1 2

( X θ − Y

) T

( X θ − Y )

    所以θ的最优解为：θ=(XTX)−1XTY。




    当然这里可能遇到一些问题，比如X必须可逆，比如求逆运算时间开销较大。具体解决方案待补充。

3.1 python实现最小二乘法

    这里的代码仅仅针对背景里的这个问题。部分参考了

回归方法及其python实现 。

1. # 通过最小二乘法直接得到最优系数，返回计算出来的系数b, m

2. def least_square_regress(points):

3. x_mat = np.mat(np.array([np.ones([len(points)]), points[:, 0]]).T) # 转为100行2列的矩阵，2列其实只有一个feature，其中x0恒为1

4. y_mat = points[:, 1].reshape(len(points), 1) # 转为100行1列的矩阵

5. xT_x = x_mat.T * x_mat

6. if np.linalg.det(xT_x) == 0.0:

7. print('this matrix is singular,cannot inverse') # 奇异矩阵，不存在逆矩阵

8. return

9. coefficient_mat = xT_x.I * (x_mat.T * y_mat)

10. return coefficient_mat[0, 0], coefficient_mat[1, 0] # 即系数b和m

    程序执行结果如下：   
    

b = 7.99102098227, m = 1.32243102276, error = 110.257383466, 相关系数 = 0.773728499888

    拟合结果如下图：

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

    有了最小二乘法以后，我们已经可以对数据点进行拟合。但由于最小二乘法需要计算X的逆矩阵，计算量很大，因此特征个数多时计算会很慢，只适用于特征个数小于100000时使用；当特征数量大于100000时适合使用梯度下降法。最小二乘法与梯度下降法的区别见

最小二乘法和梯度下降法有哪些区别？ 。

4.1. 梯度

    首先，我们简单回顾一下微积分中梯度的概念。这里参考了

方向导数与梯度 ，具体的证明请务必看一下这份材料，很短很简单的。

    讨论函数z=f(x,y)在某一点P沿某一方向的变化率问题。设函数z=f(x,y)在点P(x,y)的某一邻域U(P)内有定义，自点P引射线l到点P′(x+Δx,y+Δy)且P′∈U(P)，如下图所示。










    定义函数z=f(x,y)在点P沿方向l的方向导数为：

∂ f

∂ l

=

lim

ρ → 0

f ( x

Δ x , y

Δ y ) − f ( x , y )

ρ

，

其

中

( Δ x

) 2

( Δ y

) 2

    方向导数可以理解为，函数z=f(x,y)沿某个方向变化的速率。可以类比一下函数y=kx+b的斜率k=dydx。斜率越大，函数yy增长得越快。那么现在问题来了，函数z=f(x,y)在点P沿哪个方向增加的速度最快？而这个方向就是梯度的方向

∂ f

∂ x

i →

∂ f

∂ y

j →

    从几何角度来理解，函数z=f(x,y)表示一个曲面，曲面被平面z=c截得的曲线在xoy平面上投影如下图，这个投影也就是我们所谓的等高线。










    函数z=f(x,y)在点P(x,y)处的梯度方向与点P的等高线f(x,y)=c在这点的法向量的方向相同，且从数值较低的等高线指向数值较高的等高线。

4.2 梯度方向计算

    理解了梯度的概念之后，我们重新回到

1. 背景 中提到的例子。
2. 背景 提到，梯度下降法所做的是从图中的任意一点开始，逐步找到图的最低点。那么现在问题来了，从任意一点开始，b和m可以往任意方向"走"，如何可以保证我们走的方向一定是使误差函数Error(b,m)减小且减小最快的方向呢？回忆 4.1.梯度 中提到的结论，梯度的方向是函数上升最快的方向，那么函数下降最快的方向，也就是梯度的反方向。有了这个思路，我们首先计算梯度方向，

∂ E r r o

r

( b , m )

∂ m

=

∑

1

N

x i

( ( b

m

x i

) −

y i

)

∂ E r r o

r

( b , m )

∂ b

=

∑

1

N

( ( b

m

x i

) −

y i

)

，

x 0

恒

为

1

    有了这两个结果，我们就可以开始使用梯度下降法来寻找误差函数Error(b,m)的最低点。我们从任意的点(b,m)开始，逐步地沿梯度的负方向改变b和m的值。每一次改变，Error(b,m)都会得到更小的值，反复进行该操作，逐步逼近Error(b,m)的最低点。




    回到更一般的情况，对于每一个向量θ的每一维分量θi，我们都可以求出梯度的方向，也就是错误函数J(θ)下降最快的方向：   

∂

∂

θ j

∂

∂

θ j

1 2

∑

1

m

(

h θ

(

x

( i )

) −

y

( i )

) 2

=

∑

1

m

(

h θ

(

x

( i )

) −

y

( i )

)

x j

( i )

4.3 批量梯度下降法

    从上面的公式中，我们进一步得到特征的参数θj的迭代式。因为这个迭代式需要把m个样本全部带入计算，所以我们称之为批量梯度下降

θ j ′

=

θ j

− α

∂ J ( θ )

∂

θ j

=

θ j

− α

∑

1

m

(

h θ

(

x

( i )

) −

y

( i )

)

x j

( i )

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

1. def step_gradient(b_current, m_current, points, learningRate):

2. b_gradient = 0

3. m_gradient = 0

4. N = float(len(points))

5. for i in range(0, len(points)):

6. x = points[i, 0]

7. y = points[i, 1]

8. m_gradient += (2 / N) * x * ((b_current + m_current * x) - y)

9. b_gradient += (2 / N) * ((b_current + m_current * x) - y)

10. new_b = b_current - (learningRate * b_gradient) # 沿梯度负方向

11. new_m = m_current - (learningRate * m_gradient) # 沿梯度负方向

12. 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文件，在

1. # this is the instance used by the matplotlib classes

2. rcParams = rc_params()

    后面加上：
   

3. # fix a bug by ZZR

4. 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 随机梯度下降法

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

θ j ′

=

θ j

− α (

h θ

(

x

( i )

) −

y

( i )

)

x j

( i )

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




    这里还有一些重要的细节没有提及，比如如何确实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动画指南

ARTICLES

近期热门文章

⊙ 生成器 ：

关于生成器的那些事儿

⊙ 爬虫代理 ：

如何构建爬虫代理服务

⊙ 地理编码 ：

怎样用Python实现地理编码

⊙ nginx日志 ：

使用Python分析nginx日志

⊙ 淘宝女郎 ：

一个批量抓取淘女郎写真图片的爬虫

⊙ IP代理池 ：

突破反爬虫的利器——开源IP代理池

⊙ 布隆去重 ：

基于Redis的Bloomfilter去重（附代码）

⊙ QQ空间爬虫 ：

QQ空间爬虫最新分享，一天 400 万条数据

⊙ 匿名代理池 ：

进击的爬虫：用Python搭建匿名代理池

在公众号底部回复上述关键词可直接打开相应文章

我 们 终 将 改 变 潮 水 的 方 向

§

§

Python中文社区

www.python-cn.com

致力于成为

国内最好的Python社区

QQ群：152745094

专栏作者申请邮箱

sinoandywong@gmail.com

— Life is short,we use Python —

点击阅读原文可直接进入作者博客