机器学习算法之线性回归

来源： https://www.biaodianfu.com

线性回归是统计学总最常用的算法之一。从根本上来说，当你想表示两个变量间数学关系时，就可以使用线性回归。当你使用它时，你首先假设输出变量（有时称为响应变量、因变量或标签）和预测变量（有时称为自变量、解释变量或特征）之间存在线性关系。当然这种线性关系也可能存在于一个输出变量和数个预测变量之间。输出变量于预测变量之间存在线性关系是一个大胆的假设，同时也是一个最简单的假设。从数学表示形式来看，线性函数比非线性函数更加简单。线性模型作为最简单的参数化方法，始终值得关注。这是因为很多问题，甚至本质是非线性的问题，也可以采用线性模型解决。

在线性回归中，数据采用线性函数的方式进行数据建模，对于模型中的未知参数也采用数据进行估计。对于一个多变量的线性回归模型可以表示为如下公式：

=0+11+22+…++,=1,2,…,

其中，Y是X的线性函数，是误差项。线性则是的条件均值，在参数里是线性的。有些函数模型看起来并不一定是线性回归 但是通过代数转换可以转换为线性回归模型。对于一个简单的线性回归，可以表示为如下：

=+

对于单变量的线性回归，是目前最容易理解的线性模型，其回归式如下：

=++

线性回归是回归分析中最为广泛使用的模型， 在结果预测及函数关系的应用中较为频繁。

对于线性回归算法，我们希望从训练数据中学习到线性回归方程，即：

=+∑=1

其中，b称为偏置，为回归系数。对于上式，令0=1则上式可以表示为：

=∑=0

在线性回归模型中，其目的是 求出线性模型回归方程，即求出线性回归方程中的回归系数。线性回归的评价是指如何度量预测值（Prediction）于标签（Label）之间的接近程序，线性回归模型的损失函数可以是绝对损失（Absolute Loss）或者平方损失（Squared Loss）。其中，绝对损失函数为：

=|−̂|

其中，̂为预测值，且：̂=∑=0

平方损失函数为：

=(−̂)2

由于平方损失处处可导，通常使用平方误差作为线性回归模型的损失函数。假设有m个训练样本，每个样本中有n-1个特征，则平方误差可以表示为：

=12∑=1(−∑=0−1)2

对于如上的损失函数，线性回归的求解是希望求得平方误差最小值。

最小二乘法

最小二乘法也称作最小平方法， 最常用的是普通最小二乘法(Ordinary Least Square) ,它是一种数学中的优化方法，试图找到一个或一组估计值，使得实际值与估计值尽可能相似， 距离最小，目的是通过已有的数据来预测未知数据。一般通过一条多元一次的直线方程来计算，在二维坐标中即二元一次方程。 例如，在二维坐标中， 非常多的点分散在其中，试图绘制一条直线，使得这些分散的点到直线上的距离最小。这里的距离最小并非点到直线的垂直距离最短，而是点到直接的y轴距离最短，即通过该点并与y轴平行的直线，点到该y轴平行线与直线交点的距离最短。

最小二乘法的核心思想是通过最小化误差的平方和，试图找到最可能的函数方程。 例如，在二维坐标系中存在五个数据点(10,20), (11,23), (12,25), (13,27), (14,26), 希望找出一条该五个点距离最短的直线，根据二元一次方程：=+。将五个点分别带入该二元方程得到如下：

20=10+

23=11+

25=12+

27=13+

26=14+

由于最小二乘法是尽可能使得等号两边的方差值最小，因此：

(,)=[20−(10+)]2+[23−(11+)]2+[25−(12+)]62+[27−(13+)]2+[26−(14+)]2

因此求最小值即可通过对(,)求偏导数获得，并得一阶倒数得值为0，因此：

∂∂=1460+120−2936

∂∂=12+10−242

即得到关于求解未知变量a、b的二元一次方程：

{1460+120=293612+10=242

通过计算上述二元一次方法即得到a=0.0243, b=24.1708。 因此，在上述五个点中，通过最小二乘法得到直线方程：y=24.l708+0.0243x是使得五个点到该直线距离最小的直线。

最小二乘解法

对于线性回归模型，假设训练集中有m个训练样本，每个训练样本中有n-1个特征，可以使用矩阵的表示方法，预测函数可以表示为：=，其损失函数表示为：(−)(−)。

其中，标签Y为m×1的矩阵，训练特征X为m×n的矩阵，回归系数W 为n×1的矩阵。在最小二乘法中，对W 求导，即：

(−)(−)=(−)

令其为0，得到：

̂=()−1

Python实现：

def least_square(feature, label):

input: feature(mat):特征

label(mat):标签

w = (feature.T * feature).I * feature.T * label

最小二乘法虽然看似是一个直线方程的问题，但是在实际应用中却应用非常广泛，因为它得到的方程可以视为一个函数模型，该函数模型可以对后续的工作带来极大的便利。上述过程均是通过线性问题的求解方式进行阐述，但是在更多时候，需要解决的问题不是一个线性问题，它需要通过多项式拟合的方式进行处理，但是原理和求解方式均一致。当样本数据不断增加后，计算量会明显增加，在阶数更高时计算扯则更为复杂，为解决更多问题，基千最小二乘法衍生出了移动最小二乘法、 加权最小二乘法及偏最小二乘法等。除最小二乘法外， 梯度下降也是线性回归中的方法之一。

梯度下降法

什么是梯度下降法

梯度下降法的基本思想可以类比为一个下山的过程。假设这样一个场景：一个人被困在山上，需要从山上下来(i.e. 找到山的最低点，也就是山谷)。但此时山上的浓雾很大，导致可视度很低。因此，下山的路径就无法确定，他必须利用自己周围的信息去找到下山的路径。这个时候，他就可以利用梯度下降算法来帮助自己下山。具体来说就是，以他当前的所处的位置为基准，寻找这个位置最陡峭的地方，然后朝着山的高度下降的地方走，同理，如果我们的目标是上山，也就是爬到山顶，那么此时应该是朝着最陡峭的方向往上走。然后每走一段距离，都反复采用同一个方法，最后就能成功的抵达山谷。

在微积分里面，对多元函数的参数求∂偏导数，把求得的各个参数的偏导数以向量的形式写出来，就是梯度。比如函数f(x,y),分别对x,y求偏导数，求得的梯度向量就是(∂∂,∂∂)，简称grad f(x,y)或者▽(,)。对于在点(0,0)的具体梯度向量就是(∂∂0,∂∂0)，或者▽(0,0)。如果是3个参数的向量梯度，就是(∂∂,∂∂,∂∂),以此类推。

那么这个梯度向量求出来有什么意义呢？他的意义从几何意义上讲，就是函数变化增加最快的地方。具体来说，对于函数f(x,y),在点(0,0)，沿着梯度向量的方向就是(∂∂0,∂∂0)的方向是f(x,y)增加最快的地方。或者说，沿着梯度向量的方向，更加容易找到函数的最大值。反过来说，沿着梯度向量相反的方向，也就是−(∂∂0,∂∂0)的方向，梯度减少最快，也就是更加容易找到函数的最小值。

对于凸优化问题来说，导数为0（梯度为0向量）的点，就是优化问题的解。为了找到这个解，我们沿着梯度的反方向进行线性搜索，每次搜索的步长为某个特定的数值，直到梯度与0向量非常接近为止。上面描述的这个方法就是梯度下降法。迭代算法的步骤如下：

1. 当=0，自己设置初始点0=(01,02,⋯,0)，设置步长（也叫做学习率），设置迭代终止的误差忍耐度tol。
2. 计算目标函数在点上的梯度∇。
3. 计算+1，公式如下：+1=–∇
4. 计算梯度∇+1。如果‖∇+1‖2≤则迭代停止，最优解的取值为+1；如果梯度的二范数大于tol，那么i=i+1，并返回第（3）步。

梯度下降法的类型

基于如何使用数据计算代价函数的导数，梯度下降法可以被定义为不同的形式（various variants）。确切地说，根据使用数据量的大小（the amount of data），时间复杂度（time complexity）和算法的准确率（accuracy of the algorithm），梯度下降法可分为：

• 批量梯度下降法（Batch Gradient Descent, BGD）；
• 随机梯度下降法（Stochastic Gradient Descent, SGD）；
• 小批量梯度下降法（Mini-Batch Gradient Descent, MBGD）。

批量梯度下降法BGD

这是梯度下降法的基本类型，这种方法使用整个数据集（the complete dataset）去计算代价函数的梯度。每次使用全部数据计算梯度去更新参数，批量梯度下降法会很慢，并且很难处理不能载入内存（don’t fit in memory）的数据集。

批量梯度下降法的损失函数为：

∂()∂=∑=1(ℎ()−)

进一步得到批量梯度下降的迭代式为：

′=–∂()∂==∑=1(ℎ()−)

其中，m是训练样本（training examples）的数量。

每迭代一步，都要用到训练集所有的数据，如果样本数目很大，那么可想而知这种方法的迭代速度！

优点：全局最优解；易于并行实现；

缺点：当样本数目很多时，训练过程会很慢。

• 如果训练集有3亿条数据，你需要从硬盘读取全部数据到内存中；
• 每次一次计算完求和后，就进行参数更新；
• 然后重复上面每一步；
• 这意味着需要较长的时间才能收敛；
• 特别是因为磁盘输入/输出（disk I/O）是系统典型瓶颈，所以这种方法会不可避免地需要大量的读取。

从迭代的次数上来看，BGD迭代的次数相对较少。其迭代的收敛曲线示意图可以表示如下：

上图是每次迭代后的等高线图，每个不同颜色的线表示代价函数不同的值。运用梯度下降会快速收敛到圆心，即唯一的一个全局最小值。

下面的Python代码实现了批量梯度下降法：

def gradient_descent(alpha, x, y, ep=0.0001, max_iter=10000):

m = x.shape[0]  # number of samples 

t0 = np.random.random(x.shape[1])

t1 = np.random.random(x.shape[1])

# total error, J(theta) 

J = sum([(t0 + t1 * x[i] - y[i])  2 for i in range(m)])

# for each training sample, compute the gradient (d/d_theta j(theta)) 

grad0 = 1.0 / m * sum([(t0 + t1 * x[i] - y[i]) for i in range(m)])

grad1 = 1.0 / m * sum([(t0 + t1 * x[i] - y[i]) * x[i] for i in range(m)])

# update the theta_temp 

temp0 = t0 - alpha * grad0

temp1 = t1 - alpha * grad1

# mean squared error 

e = sum([(t0 + t1 * x[i] - y[i])  2 for i in range(m)])

if abs(J - e) <= ep:

print('Converged, iterations: ', iter, '!!!')

converged = True

J = e  # update error 

iter += 1  # update iter 

if iter == max_iter:

print('Max interactions exceeded!')

converged = True

小批量梯度下降法MBGD

小批量梯度下降法是最广泛使用的一种算法，该算法每次使用m个训练样本（称之为一批）进行训练，能够更快得出准确的答案。小批量梯度下降法不是使用完整数据集，在每次迭代中仅使用m个训练样本去计算代价函数的梯度。一般小批量梯度下降法所选取的样本数量在50到256个之间，视具体应用而定。

• 减少了参数更新时的变化，能够更加稳定地收敛。
• 能利用高度优化的矩阵，进行高效的梯度计算。

随机初始化参数后，按如下伪码计算代价函数的梯度：

这里b表示一批训练样本的个数，m是训练样本的总数。

随机梯度下降法SGD

批量梯度下降法被证明是一个较慢的算法，所以，我们可以选择随机梯度下降法达到更快的计算。随机梯度下降法的第一步是随机化整个数据集。在每次迭代仅选择一个训练样本去计算代价函数的梯度，然后更新参数。即使是大规模数据集，随机梯度下降法也会很快收敛。随机梯度下降法得到结果的准确性可能不会是最好的，但是计算结果的速度很快。在随机化初始参数之后，使用如下方法计算代价函数的梯度：

′=+1(−ℎ())

批量梯度下降会最小化所有训练样本的损失函数，使得最终求解的是全局的最优解，即求解的参数是使得风险函数最小。随机梯度下降会最小化每条样本的损失函数， 虽然不是每次迭代得到的损失函数都向着全局最优方向， 但是大的整体的方向是向全局最优解的，最终的结果往往是在全局最优解附近。

随机梯度下降法的优化过程为：

随机梯度下降是通过每个样本来迭代更新一次，如果样本量很大的情况（例如几十万），那么可能只用其中几万条或者几千条的样本，就已经将theta迭代到最优解了，对比上面的批量梯度下降，迭代一次需要用到十几万训练样本，一次迭代不可能最优，如果迭代10次的话就需要遍历训练样本10次。但是，SGD伴随的一个问题是噪音较BGD要多，使得SGD并不是每次迭代都向着整体最优化方向。

优点：训练速度快；

缺点：准确度下降，并不是全局最优；不易于并行实现。

从迭代的次数上来看，SGD迭代的次数较多，在解空间的搜索过程看起来很盲目。其迭代的收敛曲线示意图可以表示如下：

参考链接：http://sofasofa.io/tutorials/python_gradient_descent/index.php

牛顿法

除了前面说的梯度下降法，牛顿法也是机器学习中用的比较多的一种优化算法。牛顿法的基本思想是利用迭代点处的一阶导数（梯度）和二阶导数（Hessen矩阵）对目标函数进行二次函数近似，然后把二次模型的极小点作为新的迭代点，并不断重复这一过程，直至求得满足精度的近似极小值。牛顿法下降的速度比梯度下降的快，而且能高度逼近最优值。牛顿法分为基本的牛顿法和全局牛顿法。

基本牛顿法

基本牛顿法是一种基于导数的算法，它每一步的迭代方向都是沿着当前点函数值下降的方向。对于一维的情形，对于一个需要求解的优化函数f(x)，求函数的极值的问题可以转化为求导函数′()=0。对函数f(x)进行泰勒展开到二阶，得到：

()=()+′()(−)+12()(−)2

对上式求导并令其为0，则为：

′()+”()(−)=0

即得到：

=–′()”()

这就是牛顿法的更新公式。

基本牛顿法的流程

1. 给定终止误差值0≤≤1，初始点0∈ℝ ，令k=0；
2. 计算=∇()，若‖‖≤，则停止，输出∗≈ ；
3. 计算=∇2()，并求解线性方程组=−得解；
4. 令+1=+,=+1，并转2。

全局牛顿法

牛顿法最突出的优点是收敛速度快，具有局部二阶收敛性，但是，基本牛顿法初始点需要足够“靠近”极小点，否则，有可能导致算法不收敛，此时就引入了全局牛顿法。全局牛顿法的流程为：

1. 给定终止误差值0≤≤1，∈(0,1)，∈(0,0.5)，初始点0∈ℝ ，令k=0；
2. 计算=∇()，若‖‖≤，则停止，输出∗≈ ；
3. 计算=∇2()，并求解线性方程组=−得解；
4. 设是不满足下列不等式的最小非负整数m：(+)≤()+
5. 令=,+1=+,=+1，并转2。

全局牛顿法的具体实现：

def newton(feature, label, iterMax, sigma, delta):

input: feature(mat):特征

label(mat):标签

iterMax(int):最大迭代次数

sigma(float), delta(float):牛顿法中的参数

n = np.shape(feature)[1]

w = np.mat(np.zeros((n, 1)))

g = first_derivativ(feature, label, w)  # 一阶导数

G = second_derivative(feature)  # 二阶导数

m = get_min_m(feature, label, sigma, delta, d, w, g)  # 得到最小的m

w = w + pow(sigma, m) * d

print "t---- itration: ", it, " , error: ", get_error(feature, label , w)[0, 0]

在程序中，函数newton利用全局牛顿法对线性回归模型中的参数进行学习，函数newton的输入为训练特征feature、训练的目标值label、全局牛顿法的最大迭代次数iterMax以及全局牛顿法的两个参数sigma和delta。函数newton的输出是线性回归模型的参数 w。在函数 newton 中需要计算损失函数的一阶导数，计算损失函数的二阶导数，，同时需要计算最小的m值，最终根据上述的值更新权重。

def get_min_m(feature, label, sigma, delta, d, w, g):

input: feature(mat):特征

label(mat):标签

sigma(float),delta(float):全局牛顿法的参数

d(mat):负的一阶导数除以二阶导数值

w_new = w + pow(sigma, m) * d

left = get_error(feature, label , w_new)

right = get_error(feature, label , w) + delta * pow(sigma, m) * g.T * d

if left <= right:

程序实现了全局牛顿法中最小m值的确定，在函数get_min_m中，其输入为训练数据的特征 feature，训练数据的目标值 label，全局牛顿法的参数 sigma、delta、d以及损失函数的一阶导数值g。其输出是最小的m值m。在计算的过程中，计算损失函数值时使用到了get_error函数，其具体实现：

def get_error(feature, label, w):

input: feature(mat):特征

label(mat):标签

w(mat):线性回归模型的参数

return (label - feature * w).T * (label - feature * w) / 2

Armijo搜索

给定∈(0,1)，∈(0,0.5)，令步长因子= ，其中是满足下列不等式的最小非负整数：(+)≤()+

利用全局牛顿法求解线性回归模型

假设有m个训练样本，其中，每个样本有n-1个特征，则线性回归模型的损失函数为：

=12∑=1(−∑=0−1)2

若是利用全局牛顿法求解线性回归模型，需要计算线性回归模型损失函数的一阶导数和二阶导数，其一阶导数为：

∂∂=−∑=1[(−∑=0−1⋅)∗]

损失函数的二阶导数为：

∂∂∂=∑=1(⋅)

相关代码：

def first_derivativ(feature, label, w):

input: feature(mat):特征

label(mat):标签

m, n = np.shape(feature)

g = np.mat(np.zeros((n, 1)))

err = label[i, 0] - feature[i, ] * w

for j in xrange(n):

g[j, ] -= err * feature[i, j]

def second_derivative(feature):

input: feature(mat):特征

m, n = np.shape(feature)

G = np.mat(np.zeros((n, n)))

x_left = feature[i, ].T

x_right = feature[i, ]

G += x_left * x_right

first_derivativ 实现了损失函数一阶导数值的求解。在函数first_derivativ 中，其输入为训练数据的特征feature 和训练数据的目标值 label，其输出为损失函数的一阶导数g，其中g是一个n×1的向量。

second_derivative函数实现了损失函数二阶导数值的计算。在函数second_derivative中，其输入为训练数据的特征feature，输出为损失函数的二阶导数G，其中G是一个n×n的矩阵。

局部加权线性回归

在线性回归中会出现欠拟合的情况，有些方法可以用来解决这样的问题。局部加权线性回归（LWLR）就是这样的一种方法。局部加权线性回归采用的是给预测点附近的每个点赋予一定的权重，此时的回归系数可以表示为：

̂=()−1

M为给每个点的权重。

LWLR使用核函数来对附近的点赋予更高的权重，常用的有高斯核，对应的权重为：

(,)=exp(∣∣()−∣∣−22)

这样的权重矩阵只含对角元素。

局部加权线性回归的具体实现：

def lwlr(feature, label, k):

input: feature(mat):特征

label(mat):标签

k(int):核函数的系数

output: predict(mat):最终的结果

m = np.shape(feature)[0]

predict = np.zeros(m)

weights = np.mat(np.eye(m))

for j in xrange(m):

diff = feature[i, ] - feature[j, ]

weights[j,j] = np.exp(diff * diff.T / (-2.0 * k  2))

xTx = feature.T * (weights * feature)

ws = xTx.I * (feature.T * (weights * label))

predict[i] = feature[i, ] * ws

使用Scikit-Learn进行线性回归分析

这里使用SkLearn自带得数据集进行测试：

from sklearn.datasets import load_boston

from sklearn.model_selection import train_test_split

from sklearn.linear_model import LinearRegression

print(X.shape)  # 样本个数及特征数

print(boston.feature_names)  # 特征标签名字

把数据分为训练数据集和测试数据集(20%数据作为测试数据集）

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=3)

model = LinearRegression()

model.fit(X_train, y_train)

train_score = model.score(X_train, y_train)  # 模型对训练样本得准确性

test_score = model.score(X_test, y_test)  # 模型对测试集的准确性

模型优化

多项式与线性回归

当线性回归模型太简单导致欠拟合时，我们可以增加特征多项式来让线性回归模型更好的拟合数据。比如两个特征1,2，可以增加两特征的乘积1⋅2作为新特征3。我们还可以增加21作为另外一个新特征4。

在scikit-learn里，线性回归是由类sklearn.linear_model.LinearRegression实现，多项式由类sklearn.preprocessing.PolynomialFeatures实现。那么要怎么添加多项式特征呢？我们需要用一个管道把两个类串起来，即用sklearn.pipeline.Pipeline把两个模型串起来。一个Pipeline可以包含多个处理节点，在scikit-learn里，除了最后一个节点外，其他的节点都必须实现’fit()’和’transform()’方法，最后一个节点只需要实现fit()方法即可。当训练样本数据送进Pipeline里进行处理时，它会逐个调用节点的fit()和transform()方法，然后调用最后一个节点的fit()方法来拟合数据。

数据归一化

当线性回归模型有多个输入特征时，特别是使用多项式添加特征时，需要对数据进行归一化处理。比如，特征1的范围在[1，4]之间，特征2的范围在[1，2000]之间，这种情况下，可以让1除以4来作为新特征1，同时让2除以2000来作为新特征2，该过程称为特征缩放（feature scaling）。可以使用特征缩放来对训练样本进行归一化处理，处理后的特征值范围在[0，1]之间。

归一化处理的目的是让算法收敛更快，提升模型拟合过程中的计算效率。进行归一化处理后，当有个新的样本需要计算预测值时，也需要先进行归一化处理，再通过模型来计算预测值，计算出来的预测值要再乘以归一化处理的系数，这样得到的数据才是真实的预测数据。

在scikit-learn里，使用LinearRegression进行线性回归时，可以指定normalize=True来对数据进行归一化处理。具体可查阅scikit-learn文档。

优化后代码：

from sklearn.datasets import load_boston

from sklearn.model_selection import train_test_split

from sklearn.linear_model import LinearRegression

from sklearn.preprocessing import PolynomialFeatures

from sklearn.pipeline import Pipeline

print(X.shape)  # 样本个数及特征数

print(boston.feature_names)  # 特征标签名字

把数据分为训练数据集和测试数据集(20%数据作为测试数据集）

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=3)

def polynomial_model(degree=1):

polynomial_features = PolynomialFeatures(degree=degree, include_bias=False)

linear_regression = LinearRegression(normalize=True)

pipeline = Pipeline([("polynomial_features", polynomial_features), ("linear_regression", linear_regression)])

model = polynomial_model(degree=2)

model.fit(X_train, y_train)

train_score = model.score(X_train, y_train)  # 模型对训练样本得准确性

test_score = model.score(X_test, y_test)  # 模型对测试集的准确性

使用二阶多项式训练后样本分数和测试分数都提高了，看来模型得到了优化，改为三阶后，针对训练样本的分数达到了1，而针对测试样本的得分数却是负数。说明这个模型已经过拟合。