[python][机器学习][回归]线性回归

时间 2019-11-08

原文原文链接

最后一次更新日期： 2019/5/12html

机器学习系列教程，旨在解释算法的基本原理并提供算法的入门级实现思路。python

须要引入如下模块：算法

import numpy as np数组

import pandas as pdbash

import matplotlib.pyplot as plt网络

from mpl_toolkits.mplot3d import Axes3Dapp

import mathdom

主要使用的计算与绘图库可查阅如下资料：机器学习

numpy使用指南ide

pandas使用指南

matplotlib使用指南

微积分和线代的入门，推荐3blue1brown的视频

微积分的本质

线性代数的本质

点击下方连接可前往各章节

一. 问题原型与假设函数

二. 代价函数与模型评估

三. 梯度降低法优化

四. 多元回归与矩阵运算

五. 正规方程法优化

六. 非线性映射

七. 过拟合问题

一. 问题原型与假设函数

返回目录

1. 问题原型

一般存在这样的一个场景，咱们但愿给程序输入一系列咱们持有的数值，而程序输出另外一个咱们所须要的数值，好比老生常谈的房价预测问题，输入一个房屋的市中心距离、房间数、空气质量等属性（这些属性通常被称为特征），输出房屋的价格预测（这种连续的输出数值也被称为回归值）。

这样多个输入到一个输出的映射关系，能够是简单的线性关系，即每一个输入与输出的函数图像都是直线（固定其余输入），也能够是复杂的非线性关系，即函数图像是曲线。

2. 假设函数

不妨就从最简单的场景开始，只有一个特征，且假设它与目标回归值为简单的线性关系，这种关系能够表示成这样一个函数： $h_\theta (x_1)=\theta_0+\theta_1x_1$ ， $\theta_0$ 和 $\theta_1$ 是影响输出的两个参数。能够很容易的得出，函数的完整形式应该是： $h_\theta (x)=\theta_0+\theta_1x_1+\theta_2x_2+\cdots+\theta_nx_n=\sum_{j=0}^m{\theta_j x_j}, x_0=1$ ，但此处先不考虑多个特征的状况。

绘制出二维图像后能够发现，咱们只要找出一条直线使其可以最优地贴合已知的一些数据点就好了，因此线性回归是一种监督学习方式，须要咱们提供一些已知结果的数据，不然咱们没法在训练过程当中判断模型输出的优劣，也就没法进行优化了。

为方便测试，可用numpy构造用于线性回归的简单数据集：

x1=np.linspace(0,10,20)
y=3*x1+2+np.random.randn(20)*5

plt.scatter(x1,y)
plt.plot(x1,3*x1+2,c='r')
plt.xlabel('x1')
plt.ylabel('y')
复制代码

假设函数的代码实现：

#假设函数
def linear(theta0,theta1,x1):
    return theta0+theta1*x1
复制代码

如今咱们有了一个将输入映射到输出的函数，可是尚未办法评估输出结果的好坏。

二. 代价函数与模型评估

返回目录

1. 代价函数：均方偏差

之因此叫作代价函数，能够理解为由于模型输出不正确的结果，咱们须要它为此付出代价，接受惩罚，错得越离谱，惩罚越重。因此代价函数至少应该有以下性质：当结果彻底正确时，代价值为0；偏离正确结果越多，代价值越大。

能够很快想到的就是将正确的输出与错误的输出相减（简称残差）再取绝对值就能获得衡量单条记录输出质量的数值，取绝对值是为了消除正负的影响，由于在当前场景下，不管是小于仍是大于正确的回归值，都不是咱们想要的，应该同等视之，而在衡量多条记录的总体输出误差时，每每会对代价值求平均，若是有正负会相互抵消，这也不是咱们想要的。

残差绝对值是一个可用的代价函数，可是却不便于优化，因为绝对值函数非连续可导，使用该函数做为代价函数会致使没法使用梯度降低法这类依赖于梯度信息的优化方法，因此此处不采用，改用残差平方和（也称均方偏差），一样能够正确衡量输出质量，公式为： $J(\theta)=\frac{1}{2n}\sum_{i=0}^n(h_\theta (x_{i})-y_i)^2$ ，其中是样本数量。公式中之因此额外除以了2只是出于计算上的习惯，在后面进行求导时分母上的2会被约去。

代价函数的代码：

#代价函数
#p_y是predict y，即模型的输出
def cost(p_y,y):
    return 0.5*np.power(y-p_y,2).mean()
复制代码

在只有一条记录时，代价函数的变化特色：

def cost1(p_y,y):
    return np.power(y-p_y,2)

x=np.full((100,),2)
y=np.full((100,),3)
p_y=y+np.linspace(-5,5,100)
theta1=p_y/x
theta0=p_y-x
cost_=cost1(p_y,y)

fig=plt.figure(figsize=(12,3))
ax=fig.add_subplot(131,title='y=3')
ax.plot(p_y,cost_,c='r')
ax.set_xlabel('p_y')
ax.set_ylabel('cost')
ax=fig.add_subplot(132,title='y=3, x=2, theta1=1')
ax.plot(theta0,cost_,c='r')
ax.set_xlabel('theta0')
ax.set_ylabel('cost')
ax=fig.add_subplot(133,title='y=3, x=2, theta0=0')
ax.plot(theta1,cost_,c='r')
ax.set_xlabel('theta1')
ax.set_ylabel('cost')
复制代码

关于代价函数和损失函数的区别：初学者每每会对这两个相似的概念感到困惑，若是追求严谨，能够用损失函数表示单条记录的输出误差，而代价函数则是数据集中全部记录损失的平均，或者只用两个名词中的一个进行统一的表述，本文中倾向于后者，统一称为代价。

2. 评估指标：R2

代价函数是定义给程序看的，程序可以经过代价函数计算输出的误差，并在偏离正确值越远时以更快地速度矫正。但对于用户来讲，代价函数就不那么直观了，由于不一样类数据集的代价函数值域不同，咱们很难从中看出离最优结果还差多少，因此须要一个更直观的统计指标，理想的状况是值域为[0,1]，越接近彻底准确的预测则输出越接近1，能够视做对模型优劣的评分。

R2（或叫R方）就是这样的一个指标，公式为 $R^2=\frac{SSR}{SST}=\frac{SST-SSE}{SST}$ ，其中 $SSE=\sum_{i=0}^n(y_i-\hat{y}_i)^2,SST=\sum_{i=0}^n(y_i-\overline{y})^2,\hat{y_i}=h_\theta(x_i),\overline{y}=\frac{1}{n}\sum_{i=0}^ny_i$ 。

SSE即sum of error squares（残差平方和），SST即sum of total squares（总平方和），SSR即sum of regression squares（回归平方和），等于SST-SSE。这其中的含义能够观察下面这段代码绘制的图理解：

x1=np.array([1,2,3,10])
y=3*x1+2
p_y=2*x1+6
y_mean=np.full(4,y.mean())

plt.plot(x1,y,c='b')
plt.text(9,31,'y')

plt.plot(x1,p_y,c='g')
plt.text(9,26,'p_y')

plt.plot(x1,y_mean,c='r')
plt.text(9,12,'y_mean')

x1=[8,8,8,8]
y=[3*8+2,2*8+6,y.mean(),0]
plt.scatter(x1,y,s=40,c='k')
plt.text(8.2,1,'x_i')
plt.text(8.2,24,'SSE')
plt.text(8.2,17,'SSR')

plt.plot([x1[0],x1[2]],[y[0],y[2]],c='k');
plt.plot([x1[2],x1[3]],[y[2],y[3]],c='k',linestyle='--');

plt.title('R2')
plt.xlabel('x1')
plt.ylabel('y')
plt.ylim([0,33])
复制代码

图中绘制了, $\hat{y}$ (p_y), $\overline{y}$ (y_mean)对应的直线，拿其中一个数据点举例，SSE对应的是蓝绿线之间的差值，当预测值与真实值越接近，该差值越小，SST对应的是蓝红线之间的差值，SST的公式相似于方差，能够描述的离散程度，SSR则是SST-SSE，至关于绿红线之间的差值（注意，由于实际运算时取平方消除了正负，这种对应并不严格，当 $\hat{y_i}$ 处在蓝线上方的区域时，应当参考蓝线下方的等距点）。

能够看出，R2就是在计算 $\hat{y}$ 所表达的数据分布与在多大程度上符合。R2的计算等价于 1 - 代价值/方差，评分须要基准值，而此处就是以方差做为基准值，当代价值为0， $\hat{y}$ 与彻底一致时，R2达到最大值1，而当模型输出很是糟糕，代价值甚至大于方差时，R2会是负数，不过在多数状况下R2仍是保持在[0,1]范围内。

R2是回归模型最经常使用的评估指标，sklearn中也是使用该指标进行评估。在R2之上，还有调整R2，公式 $R_{adj}^2=1-\frac{(n-1)(1-R^2)}{n-p-1}$ ，其中是样本数，是特征数。该指标将样本数归入了考虑范围，在相同R2的状况下，会倾向于样本数更少的模型，但易出现除零错误，通常可将调整R2与R2一同求出，在R2相同时再以调整R2做为进一步选择的参考。

代码实现示例：

def r2(y,p_y):
  SSE=(np.power(y-p_y,2)).sum()
  SST=(np.power(y-y.mean(),2)).sum()
  return 1-SSE/SST

def r2_adjusted(r2,p,n):
    return 1-(1-r2)*(n-1)/(n-p-1)
复制代码

三. 梯度降低法优化

返回目录

1. 批量梯度降低

观察一下上一章节代价函数的函数图像，能够发现，函数图像老是只有一个最低点，咱们但愿模型在已知数据上的输出尽量正确，即代价值要尽量低，只须要求出代价函数的最小值所在位置，该位置会肯定全部已知数据的指望输出，而后根据假设函数便能计算出正确的 $\theta$ 。

固然，上图被限制为只有一条记录，但能够肯定的是，即便应用在一个数据集上，代价函数依旧是只有一个最小值的凸函数，有兴趣的能够查找相关资料了解如何证实。

咱们实现代价函数的代码时输入是模型的预测输出，但实际进行优化时，咱们能改变的是模型的参数，因此须要将假设函数考虑进去，完整的表示出优化目标函数与优化项之间的关系。

下面一段代码将使用以前生成的数据集，并将 $\theta_0$ 和 $\theta_1$ 的变化置于同一个图中：

theta0=np.linspace(-8,12,100)
theta1=np.linspace(-7,13,100)
theta0,theta1=np.meshgrid(theta0,theta1)
cost_=[]
for t0,t1 in zip(theta0.flat,theta1.flat):
    cost_.append(cost(linear(t0,t1,x1),y))
cost_=np.array(cost_).reshape((100,100))

fig=plt.figure(figsize=(9,3))
ax=fig.add_subplot(121)
levels=[1,2,4,8,16,32,64,128]
c=ax.contour(theta0,theta1,cost_,levels,colors='k',linewidths=0.5)
ax.clabel(c,fontsize=10)
cf=ax.contourf(theta0,theta1,cost_,levels,cmap='YlOrRd')
ax.set_xlabel('theta0')
ax.set_ylabel('theta1')
ax=fig.add_subplot(122,projection='3d')
ax.plot_surface(theta0,theta1,cost_,cmap='coolwarm',alpha=0.8)
ax.contour(theta0,theta1,cost_,levels,colors='k',linewidths=0.5)
ax.set_xlabel('theta0')
ax.set_ylabel('theta1')
ax.set_zlabel('cost')
ax.view_init(45,135)
复制代码

在线性回归中，有这么一种直接求最优解的优化方法：正规方程法，后面会提到，但在此处先介绍另外一种更为通用的优化方法：梯度降低法。

梯度降低法具有更强的泛用性，在不少其余机器学习算法中也会用到，由于比起只有一个最小值，更多的状况是优化目标函数有多个极小值，并且特征每每也远不止一个，这种状况下直接求最优解会变得几乎不可能。所以，经过屡次迭代优化逐步逼近极值点是更经济可行的作法，这种作法虽然不能保证获得最优解，但不少时候极优解带来的准确率已经足够让人满意了。这类优化算法不止梯度降低法一种，还有坐标降低法、牛顿法、拟牛顿法等，此处不作介绍。

审视一下上面从新绘制的函数图像，梯度降低法的思想很简单：选择合适的方向和速度，一步步走下山坡，直到局部最低点。代价函数的图像在参数与代价构成的数据空间中表现为低一维的“坡面”，具备至少一个局部最低点（在线性回归中只有一个全局最低点），当选取一个起始点后，可经过数学方法估计可以较快逼近局部最低点的移动方向和速度，而后不断迭代优化。

在只有一个参数的简单场景下（仅 $\theta_1$ ），经过微积分知识计算代价关于参数的导数，咱们能够得出该参数应该增长仍是减少，以及变化多少的参考量。当处于最低点左侧时，导数是负的，参数应该增大；当处于最低点右侧时，导数是正的，参数值应该减少；同时越接近最低点，导数的绝对值越小，这对应着一开始远离最低点时，参数值须要更快的改变，而接近最低点后，速度应该放缓以期得到稳定的结果。导数的绝对值只能给出一个变化量的参考值，由于导数所表示是自变量在趋向于0的极小变化量下因变量的变化趋势，而咱们在进行优化时不可能一次只更新一个极小的量，那样太慢，因此通常会将负导数乘上一个学习率 $\alpha$ 做为每次的更新量，经过控制学习率的大小来调整准确度与变化速度之间的平衡，学习率过大会致使变化过分，难以准确地收敛至极值点甚至离极值点原来越远，而学习率太小会致使优化速度很慢。

对 $J(\theta)=\frac{1}{2n}\sum_{i=0}^n(h_\theta (x_i)-y_i)^2$ 求导，获得 $\frac{\partial J(\theta)}{\partial\theta}=\frac{1}{n}\sum_{i=0}^n(h_\theta (x_i)-y_i)x_i$ ，此处使用了复合函数求导的链式法则，其中是 $h_\theta(x_i)$ 关于 $\theta$ 的导数。而后每次迭代执行以下更新： $\theta^{new}=\theta-\alpha\frac{\partial J(\theta)}{\partial\theta}$ ，持续优化直至收敛到极值点（到什么程度算是收敛其实没有一个很是明确的标准，通常认为变化量小于必定的阈值，即继续更新对结果的影响微乎其微时，就能够认为是收敛了）。

学习率对优化效果的影响能够观察下面的示例：

x=np.linspace(0,10,20)
y=3*x1+np.random.randn(20)*5
theta_range=np.linspace(-7,13,21)
cost_range=[]
for t in theta_range:
    cost_range.append(cost(linear(0,t,x),y))
cost_range=np.array(cost_range)
theta_=10

#单次优化
def optimize(theta,x,y,alpha):
    p_y=linear(0,theta,x)
    theta-=alpha*((p_y-y)*x).sum()/len(x)
    return theta

#迭代优化
def train(theta,x,y,alpha,iters=10):
    optimize_h=[theta]
    for i in range(iters):
        theta=optimize(theta,x,y,alpha)
        optimize_h.append(theta)
    return optimize_h

#绘制变化过程
def plot_change(alpha,axes):
    axes.plot(theta_range,cost_range)
    optimize_h=train(theta_,x,y,alpha)
    for i in range(len(optimize_h)-1):
        theta_1=optimize_h[i]
        theta_2=optimize_h[i+1]
        cost_1=cost(linear(0,theta_1,x),y)
        cost_2=cost(linear(0,theta_2,x),y)
        axes.plot([theta_1,theta_2],[cost_1,cost_2],c='r')
        axes.scatter([theta_1,theta_2],[cost_1,cost_2],c='r',s=20)
    axes.set_xlim([-8,14])
    axes.set_ylim([-100,1900])
        
#不一样学习率的更新
fig=plt.figure(figsize=(14,3))
ax=fig.add_subplot(141,title='alpha=0.001')  
plot_change(0.001,ax)
ax=fig.add_subplot(142,title='alpha=0.01')  
plot_change(0.01,ax)
ax=fig.add_subplot(143,title='alpha=0.05')  
plot_change(0.05,ax)
ax=fig.add_subplot(144,title='alpha=0.06')  
plot_change(0.06,ax)
复制代码

由图可知，该测试用例中学习率取0.001更新太慢，0.01比较合适，0.05偏大出现了振荡但依旧能收敛，0.06过大没法收敛。

如今再来看两个参数的状况（ $\theta_0$ 和 $\theta_1$ ），多元函数中因变量关于每一个自变量的导数被称为偏导数，而全部偏导数组成的向量便是因变量关于自变量的梯度： $\nabla J(\theta)=\left[\begin{matrix} \frac{\partial J(\theta)}{\partial \theta_0}\\ \frac{\partial J(\theta)}{\partial \theta_1}\\ \end{matrix}\right]$ ，其中 $\frac{\partial J(\theta)}{\partial \theta_0}=\frac{1}{n}\sum_{i=0}^n(h_\theta (x_i)-y_i), \frac{\partial J(\theta)}{\partial \theta_1}=\frac{1}{n}\sum_{i=0}^n(h_\theta (x_i)-y_i)x_i$ 。多元函数在某点处的梯度向量指出了该点处函数值增加最快的方向，同时其模长可以体现增加量。在参数变多的状况下，咱们要作的事情并无复杂多少，只要根据当前位置的偏导数更新每个对应的参数便可，每次迭代执行以下更新： $\theta_j^{new}=\theta_j-\alpha\frac{\partial J(\theta)}{\partial\theta_j}$ 。

梯度降低的实现代码以下：

def grad_desc(theta0,theta1,x1,y,alpha):
  p_y=linear(theta0,theta1,x1)
  theta0-=alpha*(p_y-y).sum()/len(y)
  theta1-=alpha*((p_y-y)*x1).sum()/len(y)
  return theta0,theta1
复制代码

2. 随机梯度降低

上一小节中介绍的梯度降低法在每次计算梯度时都使用所有的训练样本，这被称为批量梯度降低，这种作法在数据量较大时计算梯度会变得很慢，迭代优化的速度也会所以变慢。

应对这种状况的一种方法就是使用随机梯度降低，每次优化只取用训练样本的一个随机子集，牺牲必定的准确度以大幅提升迭代速度（随机梯度降低一开始是指每次只使用一个训练样本的梯度降低，但后来都被用来指代使用随机子集的小批量梯度降低）。

为保证每一个训练样本都被均匀的覆盖到，不采用每次都随机取n个样本的作法，而是每次将训练集随机排序并分割为大小n的子集，使用完全部子集算做一轮迭代，而后重复。子集的大小，经常使用的默认设置为200左右。

牺牲必定的准确度还带来另外一个好处，在有多个极值点的数据空间中，每每还会存在不少鞍点，当批量梯度降低的起始点正好选取在鞍点的稳定方向上，就会致使最终收敛至鞍点而不是极值点，而随机梯度降低的优化方向会产生一些偏移，导致从不稳定方向上跌落鞍部。

至此，模型已经可以完成训练和预测，完整的实现代码和使用示例以下：

#线性回归
class LinearRegression:
    
    def __init__(self,learning_rate=0.001,iter_max=100,batch_size=256):
        self.learning_rate_=learning_rate
        self.iter_max_=iter_max
        self.batch_size_=batch_size
        
    #初始化参数
    def init_theta(self):
        self.theta0_=0
        self.theta1_=0
    
    #假设函数
    def linear(self,theta0,theta1,x1):
        return theta0+theta1*x1
    
    #代价函数
    def cost(self,p_y,y):
        return 0.5*np.power(y-p_y,2).mean()
    
    #预测
    def predict(self,x1):
        return self.linear(self.theta0_,self.theta1_,x1)
    
    #梯度降低
    def grad_desc(self,theta0,theta1,x1,y,alpha):
        p_y=self.linear(theta0,theta1,x1)
        theta0-=alpha*(p_y-y).sum()/len(y)
        theta1-=alpha*((p_y-y)*x1).sum()/len(y)
        return theta0,theta1
    
    #训练
    def train(self,x1,y):
        self.init_theta()
        #子集数量
        batches_n=math.ceil(len(y)/self.batch_size_)
        #优化历史
        self.theta_h_=[[self.theta0_,self.theta1_]]
        self.cost_h_=[self.cost(self.predict(x1),y)]
        #迭代
        for i in range(self.iter_max_):
            #随机排序
            idx=np.random.permutation(len(y))
            x1,y=x1[idx],y[idx]
            #遍历子集
            for j in range(batches_n):
                x1_=x1[j*self.batch_size_:(j+1)*self.batch_size_]
                y_=y[j*self.batch_size_:(j+1)*self.batch_size_]
                self.theta0_,self.theta1_=self.grad_desc(
                    self.theta0_,self.theta1_,x1_,y_,self.learning_rate_
                    )
            #记录当前参数和代价
            self.theta_h_.append([self.theta0_,self.theta1_])
            self.cost_h_.append(self.cost(self.predict(x1),y))
        self.theta_h_=np.array(self.theta_h_)
        self.cost_h_=np.array(self.cost_h_)
    
    #R方
    def r2(self,p_y,y):
        SSE=(np.power(y-p_y,2)).sum()
        SST=(np.power(y-y.mean(),2)).sum()
        return 1-SSE/SST
    
    #评分
    def score(self,x1,y):
        p_y=self.predict(x1)
        return self.r2(p_y,y)

#使用示例
x1=np.linspace(0,10,50)
y=3*x1+2+np.random.randn(50)*5

model=LinearRegression()
model.train(x1,y)

plt.title("score: %f"%model.score(x1,y))
plt.scatter(x1,y)
plt.plot(x1,model.predict(x1),c='r')
plt.xlabel('x1')
plt.ylabel('y')
复制代码

四. 多元回归与矩阵运算

返回目录

1. 改进假设函数

如今咱们的模型有一个很大的缺陷：只支持一个特征输入，若是须要支持更多的特征，按本来的写法须要增长更多独立的参数，但咱们不肯定具体使用时会有几个特征因此这种硬编码的方式是不可取的，更有效的作法是将参数向量化，简单地说就是以数组的形式存储。将参数向量化后，本来的许多运算就能够改写为矩阵运算的形式以进一步提高性能，numpy底层对矩阵运算作了充足的优化，在使用numpy时混杂大量对数组的循环遍历会发挥不出完整的性能。

能够说，凡是各元素相乘再总体求和的运算很大一部分均可以转换为矩阵运算。

假设函数的原写法：

#假设函数
    def linear(self,theta0,theta1,x1):
        return theta0+theta1*x1
复制代码

原完整公式： $h_\theta (x)=\sum_{j=0}^m \theta_i x_i, x_0=1$ ，如今能够改写为矩阵形式： $h_\theta (X)=X\theta$ ，其中 $\theta = \left [ \begin{matrix} \theta_0 \\ \vdots \\ \theta_m \\ \end{matrix} \right ]$ ， $X = \left [ \begin{matrix} 1 & x_{0 1} & \cdots & x_{0 m} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n 1} & \cdots & x_{n m} \\ \end{matrix} \right]$ ，对应样本数，对应特征数+1（额外增长了首列常数列），计算时按行取出每一个样本的特征向量，与 $\theta$ 做点积，整个运算返回一个长度的预测值向量。值得一提的是，此处的写法和常见的写法略有不一样，可能 $h_\theta (x)=\theta x$ 这样的形式更常见，区别在于此处给出的是做用于数据集的公式，且为了保持行对应记录、列对应属性这样的数据格式，没有对进行转置，而相应的，公式中做矩阵乘法的先后两项就颠倒了。

另外一种写法是 $h_\theta (X)= Xw+b$ ，其中 $w = \left [ \begin{matrix} w_0 \\ \vdots \\ w_m \\ \end{matrix} \right ]$ ， $X = \left [ \begin{matrix} x_{0 0} & \cdots & x_{0 m} \\ \vdots & \ddots & \vdots \\ x_{n 0} & \cdots & x_{n m} \\ \end{matrix} \right ]$ ，对应样本数，对应特征数，是标量，和分别是权重和偏置的意思，这种写法在神经网络中广泛使用，在优化时须要将和分开处理。

此处为了和后面正规方程法的处理方式统一，将采用第一种形式，须要对填充常量列。

改进后的代码以下：

#初始化参数
    def init_theta(self,m):
        self.theta_=np.zeros(m)
        
    #填充常量列
    def fill_x0(self,X):
        return np.insert(X,0,1,axis=1)
    
    #假设函数
    def linear(self,X,theta):
        return np.dot(X,theta)
复制代码

代价函数以预测值 $\hat{y}$ 和真实值做为输入，于是不须要改动，但能够稍做调整以适应批量计算。

2. 改进梯度降低

原先梯度计算公式的通常形式为： $\frac{\partial J(\theta)}{\partial \theta_j}=\frac{1}{n}\sum_{i=0}^n(h_\theta (x_{i})-y_{i})x_{i j}$ ，改写为矩阵形式为 $\nabla J(\theta) = \frac{1}{n} X^T(h_\theta (X)-y)$ 。

关于矩阵式如何写有个较为简单的思考方式，就是观察形状，形状的变化遵循以下规则 $(m,k) \cdot (k,n)=(m,n)$ （两个矩阵相乘时，左矩阵的第二个维度和右矩阵的第一个维度会做为抽取向量的维度，最终会有个组合，而每组向量在执行点积后会收缩为标量），好比此处 $\nabla J(\theta)$ 应当是一个长度的向量（与 $\theta$ 长度相同，是特征数+1），而是的矩阵， $h_\theta (X)$ 和均为长度的向量，如此一来，应当沿长度的维度抽取向量， $(1,n) \cdot (n,m)=(1,m)$ 或 $(m,n) \cdot (n,1)=(m,1)$ 都是可行的转换，上面便是采用了第二种。

每次迭代执行以下优化： $\theta^{new} = \theta -\alpha \nabla J(\theta)$ ，改进后的完整代码和使用示例以下：

#线性回归
class LinearRegression:
    
    def __init__(self,learning_rate=0.001,iter_max=100,batch_size=256):
        self.learning_rate_=learning_rate
        self.iter_max_=iter_max
        self.batch_size_=batch_size
        
    #初始化参数
    def init_theta(self,m):
        self.theta_=np.zeros(m)
        
    #填充常量列
    def fill_x0(self,X):
        return np.insert(X,0,1,axis=1)
    
    #假设函数
    def linear(self,X,theta):
        if X.shape[1]==theta.shape[0]-1:
            return np.dot(X,theta[1:])+theta[0]
        else:
            return np.dot(X,theta)
    
    #代价函数
    def cost(self,p_y,y):
        return 0.5*np.power(y-p_y,2).mean(axis=0)
    
    #预测
    def predict(self,X):
        return self.linear(X,self.theta_)
    
    #梯度降低
    def grad_desc(self,X,y,theta,alpha):
        p_y=self.linear(X,theta)
        theta-=alpha*np.dot(X.T,p_y-y)/len(y)
    
    #训练
    def train(self,X,y):
        X=self.fill_x0(X)
        self.init_theta(X.shape[1])
        #子集数量
        batches_n=math.ceil(len(y)/self.batch_size_)
        #优化历史
        self.theta_h_=[self.theta_]
        self.cost_h_=[self.cost(self.predict(X),y)]
        #迭代
        for i in range(self.iter_max_):
            #随机排序
            idx=np.random.permutation(len(y))
            X,y=X[idx],y[idx]
            #遍历子集
            for j in range(batches_n):
                X_=X[j*self.batch_size_:(j+1)*self.batch_size_]
                y_=y[j*self.batch_size_:(j+1)*self.batch_size_]
                self.grad_desc(
                    X_,y_,self.theta_,self.learning_rate_
                    )
            #记录当前参数和代价
            self.theta_h_.append(self.theta_)
            self.cost_h_.append(self.cost(self.predict(X),y))
        self.theta_h_=np.array(self.theta_h_)
        self.cost_h_=np.array(self.cost_h_)
    
    #R方
    def r2(self,p_y,y):
        SSE=(np.power(y-p_y,2)).sum()
        SST=(np.power(y-y.mean(),2)).sum()
        return 1-SSE/SST
    
    #评分
    def score(self,X,y):
        p_y=self.predict(X)
        return self.r2(p_y,y)

#使用示例
X=np.random.randn(1000,10)
y=np.full(1000,10)+np.random.randn(1000)*5
for i in range(X.shape[1]):
    y=y+(i+1)*X[:,i]
    
model=LinearRegression(learning_rate=0.1)
model.train(X,y)

plt.title("score: %f"%model.score(X,y))
plt.plot(range(len(model.cost_h_)),model.cost_h_,c='r')
plt.xlabel('iter')
plt.ylabel('cost')
复制代码

3. 特征归一化

如今咱们改用真实数据集看一下模型的效果如何，这里将会使用波士顿房价数据集，该数据集经过sklearn也可加载。

#使用示例：波士顿房价数据集
f = open('D:\\training_data\\used\\boston_house_price.txt')
buf = pd.read_table(f,header=None,delim_whitespace=True)
buf.columns=['CRIM','ZN','INDUS','CHAS','NOX','RM','AGE','DIS',
             'RAD','TAX','PTRATIO','B','LSTAT','MEDV']
X,y=buf.values[:,:-1],buf.values[:,-1]

model=LinearRegression(learning_rate=1e-6,iter_max=10000)
model.train(X,y)

plt.title("score: %f"%model.score(X,y))
plt.plot(range(len(model.cost_h_)),model.cost_h_,c='r')
plt.xlabel('iter')
plt.ylabel('cost')
复制代码

而后发现结果并不如预期，学习率很是很差控制，要么没法收敛，要么优化极其缓慢，不管怎么调整，评分都无法达到较高的水准，这是为何？

实际上是不一样特征取值范围的差别致使的。回想一下梯度降低的公式，在计算 $\frac{\partial J(\theta)}{\partial \theta_j}$ 时每一个样本的偏差都乘上了 $x_{i j}$ ，这样获得的梯度在 $\theta_j$ 方向上的份量大小会与特征的取值范围相关，取值范围大的特征对应的梯度份量也老是较大。

咱们能够选取数据集中两个取值范围差别较大的特征来观察一下cost的分布特色和优化过程：

#选取两个特征
X_=X[:,[0,-2]]

#cost分布（这里的计算使用了一些矩阵运算的技巧）
theta1=np.linspace(-20,40,100)
theta2=np.linspace(-20,20,100)
theta1,theta2=np.meshgrid(theta1,theta2)
theta=np.r_[theta1.reshape((1,-1)),theta2.reshape((1,-1))]
p_y=model.linear(X_,theta)
cost_=model.cost(p_y.T,y).reshape((100,100))

def mixed_plot(ax,learning_rate,theta_start=[30.,-15.]):
    ax.set_title('learning_rate: %f'%learning_rate)
    #cost等高线图
    c=ax.contour(theta1,theta2,cost_,colors='k',linewidths=0.5)
    ax.set_xlabel('theta1')
    ax.set_ylabel('theta2')
    #cost最小值
    theta1_min=theta1.ravel()[cost_.argmin()]
    theta2_min=theta2.ravel()[cost_.argmin()]
    ax.scatter(theta1_min,theta2_min,s=100,c='r',marker='*')
    #cost变化
    theta=np.array(theta_start)
    for i in range(100):
        theta_old=theta.copy()
        model.grad_desc(X_,y,theta,alpha=learning_rate)
        ax.plot([theta_old[0],theta[0]],
                [theta_old[1],theta[1]])
        
fig=plt.figure(figsize=(9,3))
ax=fig.add_subplot(121)
mixed_plot(ax,1e-6)
ax=fig.add_subplot(122)
mixed_plot(ax,1.4e-5)
复制代码

两个特征和，的取值范围更大，从等高线图中能够看出，对应的 $\theta_2$ 方向上坡面更陡峭，带来的结果是 $\theta_2$ 方向上的梯度份量老是偏大。若是咱们的起点在 $\theta_1$ 方向上距最低点较远时，问题就来了：一开始取用较小的学习率，咱们可能会发如今 $\theta_1$ 方向上优化得太慢，而 $\theta_2$ 方向很快地收敛了，而后开始调大学习率， $\theta_1$ 方向上的优化速度还没有达到满意值，学习率对 $\theta_2$ 而言却已通过大了，再增长下去 $\theta_2$ 方向上就没法正常收敛。

而解决这个问题最简单好用的办法就是将特征缩放到一样的取值范围，通常长度为1的区间，也称特征归一化，这样梯度在各方向上的份量就不容易出现差别过大的状况。归一化改变了特征的量纲，但却保留了数值的大小关系，不会所以致使模型不能正常训练。

归一化有多种方式，此处会采用最大最小值归一化的方式，公式为 $z_j=\frac{x_j-min(x_j)}{max(x_j)-min(x_j)}$ ，是全部样本的特征组成的向量，和是缩放的参考值，通常是由训练集数据计算获得，应当保存下来，预测时也用该参考值进行归一化，由于以后若是基于用于预测的数据从新计算参考值，可能会出现一样的数值缩放后和训练时不一致的状况。

代码实现以及上方示例的从新测试：

#缩放参照统计量
def scaling_ref(X):
    ref=pd.DataFrame()
    ref['min']=X.min(axis=0)
    ref['max']=X.max(axis=0)
    return ref
  
#最大最小值缩放
def minmax_scaling(X,ref):
    return (X-ref['min'].values)/(ref['max']-ref['min']).values

#执行特征缩放
ref=scaling_ref(X_)
X_=minmax_scaling(X_,ref)

theta1=np.linspace(-40,80,100)
theta2=np.linspace(-15,65,100)
theta1,theta2=np.meshgrid(theta1,theta2)
theta=np.r_[theta1.reshape((1,-1)),theta2.reshape((1,-1))]
p_y=model.linear(X_,theta)
cost_=model.cost(p_y.T,y).reshape((100,100))

fig=plt.figure(figsize=(9,3))
ax=fig.add_subplot(121)
mixed_plot(ax,1,[50.,0.])
ax=fig.add_subplot(122)
mixed_plot(ax,2.2,[50.,0.])
复制代码

能够看到，学习率相对容易设置了，优化的速度也快了不少，虽然未能达到最理想的状态（即等高线轮廓是圆形）。这主要是由于特征的取值范围不是惟一影响梯度总体大小的缘由，特征对预测值的重要性也会影响。由图能够看出，特征1的重要性是不如特征2的，若是将优化后的 $\theta$ 输出，会发现 $\theta_1$ 的绝对值要小于 $\theta_2$ ，二者相吻合。

从新测试房价数据：

X,y=buf.values[:,:-1],buf.values[:,-1]
ref=scaling_ref(X)
X=minmax_scaling(X,ref)

model=LinearRegression(learning_rate=0.3,iter_max=1000)
model.train(X,y)

plt.title("score: %f"%model.score(X,y))
plt.plot(range(len(model.cost_h_)),model.cost_h_,c='r')
plt.xlabel('iter')
plt.ylabel('cost')
复制代码

这一次咱们获得了理想的结果。

五. 正规方程法优化

返回目录

1. 矩阵求解线性方程组

下面介绍以前提到过的，用于线性回归一次性求解的优化方法。该方法其实是基于用矩阵解线性方程组的思想，一种线性代数的常见应用方式。

有以下一个方程组，能够将其转换为矩阵表达：

\left\{
\begin{array}{c}
4x+2y=1 \\
3y+z=2 \\
x+5y+2z=3 \\
\end{array}
\right.
\rightarrow
\left\{
\begin{array}{c}
4x+2y+0z=1 \\
0x+3y+1z=2 \\
1x+5y+2z=3 \\
\end{array}
\right.
\rightarrow
\left[
\begin{matrix}
4 & 2 & 0 \\
0 & 3 & 1 \\
1 & 5 & 2 \\
\end{matrix}
\right ]
\left [
\begin{matrix}
x \\
y \\
z \\
\end{matrix}
\right ]
=
\left [
\begin{matrix}
1 \\
2 \\
3 \\
\end{matrix}
\right ]

那么该如何求解呢？经过线性代数的知识咱们能够知道，若是一个向量通过一个矩阵变换后获得另外一个向量，那么只要对获得的向量进行逆矩阵变换就能够还原为原向量。

\left [
\begin{matrix}
x \\
y \\
z \\
\end{matrix}
\right ]
=
\left [
\begin{matrix}
4 & 2 & 0 \\
0 & 3 & 1 \\
1 & 5 & 2 \\
\end{matrix}
\right ]^{-1}
\left [
\begin{matrix}
1 \\
2 \\
3 \\
\end{matrix}
\right ]

2. 正规方程法求解参数

正规方程法之因此可用是基于两个事实：线性回归的代价函数仅有一个最小值；极值点处全部自变量的偏导数都为0。基于以上事实能够构建以下线性方程组：

\left\{
\begin{array}{c}
\frac{\partial J(\theta)}{\partial \theta_0}=0 \\
\vdots \\
\frac{\partial J(\theta)}{\partial \theta_j}=0 \\
\end{array}
\right.
\rightarrow
\frac{1}{n}X^T(h_\theta (X)-y)=\vec{0}

梯度的矩阵计算式前面有给出过，只要让其等于0向量就好了。下面是求解的方式：

\frac{1}{n}X^T(h_\theta (X)-y)=\vec{0}
\rightarrow
X^T(X \theta-y) = \vec{0}

\rightarrow
X^T X\theta = X^T y
\rightarrow
\theta=(X^T X)^{-1}X^Ty

正规方程求解的开销主要在矩阵的求逆，稠密矩阵求逆的复杂度为，其中在此处为特征数，也就是在特征数量不少的状况下，使用正规方程求解开销会很大，并且另外一个问题是矩阵的逆可能不存在，这对应着方程组是无解的，这种状况大可能是因为存在无用的特征。在特征数很少的状况下，用正规方程求解是一个很好的选择。

正规方程法的代码实现：

#添加进线性回归类
class LinearRegression:

    #正规方程
    def norm_equa(self,X,y):
        return np.dot(np.linalg.inv(np.dot(X.T,X)),np.dot(X.T,y))
    
    #训练
    def train_ne(self,X,y):
        X=self.fill_x0(X)
        self.theta_=self.norm_equa(X,y)

model=LinearRegression()
model.train_ne(X,y)
print("score: %f"%model.score(X,y))

#与sklearn对照
from sklearn.linear_model import LinearRegression
sk_model = LinearRegression()
sk_model.fit(X,y)
print("sklearn score: %f"%sk_model.score(X,y))

score: 0.740643
sklearn score: 0.740643
复制代码

与sklearn的线性回归对照，结果是同样的。使用正规方程求解时无需进行特征归一化。

六. 非线性映射

返回目录

1. 将非线性转换为线性

到如今为止，咱们一直都在假设特征与预测值之间的关系是线性的，但实际状况是，非线性关系每每远多于线性关系，那咱们的模型可以处理非线性问题吗？

线性回归模型自己是没有处理非线性问题的能力的，可是咱们可使用数学上常见的手法：转换，将非线性问题转换为线性的。譬如以下的关系，很明显，与的关系是非线性的，它的函数图像不会是一条直线的，可是，若是将和看做自变量和，那自变量与因变量的关系就成了，显然是线性的，那咱们只要基于这个式子求解线性回归参数，以后模型进行预测时也采用一样的手法转换就好了。

2. 多项式特征映射

那么该如何去选择映射的函数呢？非线性函数有无数种，不可能将全部可能都尝试一遍，而是要选择尽量简单通用的映射规则，咱们没法去猜想真实的关系中存在哪些非线性函数，但应该保证咱们所选的规则可以随着映射规模的拓展，实现对原关系愈来愈好的近似。

多项式映射是最经常使用的一种规则，其有效性从泰勒公式就能够看出，理论上是能够无限近似任意非线性关系的。该映射的规则就是计算出最高次的多项式的全部项，映射后的特征数量为 $m^1+m^2+\cdots+m^h=\frac{m(m^h-1)}{m-1}$ ，好比有特征，最高次为的映射应为。

多项式特征映射的代码以下：

#多项式映射
def polynomial_mapping(X,h,cross=True,size_limit=1e8):
    #样本数和特征数
    n,m=X.shape[0],X.shape[1]
    #映射后特征数
    if cross==True:
        m_new=m*(m**h-1)//(m-1)
    else:
        m_new=m*h
    #检查映射数量是否过于庞大
    if m_new*X.shape[0]>size_limit:
        raise Exception("The array size after mapping is over limit.")
    #根据是否含组合项用不一样的方式处理
    mapping=np.zeros((n,m_new))
    if cross==True:
        #指数序列
        exponent=[np.zeros(m)]
        j_new=0
        while len(exponent)>0:
            buf=exponent.pop(0)
            #达到次数上限
            if buf.sum()>=h:
                continue
            #每次在上一轮的指数序列上取一位+1
            for j in range(m):
                exponent_=buf.copy()
                exponent_[j]+=1
                exponent.append(exponent_)
                mapping[:,j_new]=np.power(X,exponent_).prod(axis=1)
                j_new+=1
    else:
        for j in range(m):
            for k in range(h):
                mapping[:,k*m+j]=np.power(X[:,j],k+1)
    return mapping
复制代码

多项式特征映射最大的问题在于随着最高次的提高，完整映射后的特征数量会爆炸式的增加，这个缺陷限制了线性回归模型在处理非线性问题上的能力。为提升可用性，能够提供一个不生成不一样特征之间组合项的选项，用于生成不完整的映射，由于可能刚好当前数据就不侧重组合项关系，这时若是还不得不生成完整的映射就会带来没必要要的困难。

为了更清楚的看到特征映射的效果，下面会使用随机生成的简单数据集测试：

x=(np.random.rand(50)-0.5)*6
x=np.sort(x)
y=np.sin(x)+np.cos(x)+np.random.randn(50)*0.2
X=x.reshape((-1,1))

def subplot(ax,h):
    X_=polynomial_mapping(X,h,False)
    model.train_ne(X_,y)
    ax.set_title('h=%d, score=%f'%(h,model.score(X_,y)))
    ax.scatter(x,y)
    ax.plot(x,model.predict(X_),c='r')
    ax.set_xlabel('x')
    ax.set_ylabel('y')

fig=plt.figure(figsize=(13.5,3))
ax=fig.add_subplot(131)
subplot(ax,h=1)
ax=fig.add_subplot(132)
subplot(ax,h=2)
ax=fig.add_subplot(133)
subplot(ax,h=3)
复制代码

原函数是由三角函数组成的，并不包含多项式，可是却能经过多项式去近似，随着映射最高次的不断提升，预测结果也愈来愈接近原函数。

注意，若是在多项式映射后用梯度降低求解，须要进行特征归一化，由于特征的取值范围差距在映射后会变得更大。

七. 过拟合问题

返回目录

1. 交叉验证

前面的内容中，咱们在训练和测试模型时都使用了同一个数据集，这会陷入一个误区，咱们会误觉得模型的效果已经足够好了，但咱们的模型在投入使用后处理的都是新的数据，若是在新的数据上表现得不好，就没有意义了。咱们须要将训练和测试用的数据拆分开来，用训练集去训练，而后用测试集去评估，这样才能正确地评判模型的好坏，这种作法称为交叉验证。

最经常使用的拆分方式是8:2拆分，原始数据抽取4/5做为训练集，1/5做为测试集。固然，为了不偶然性，有更严谨的方式：k折交叉验证，即将数据集平均分为k份，每次取其中一份做为测试集，其他做为训练集，而后训练模型，最后在全部组合中选择效果最好的。实际使用时用哪一种方式取决于项目的时间和对准确度的要求。

下面是八二分的代码实现（numpy和pandas两个版本）：

#numpy
def split_train_test(X,y,frac=0.8):
    idx=np.arange(y.shape[0])
    test_size=int(y.shape[0]*(1-frac))
    test_idx=np.random.choice(idx,size=test_size,replace=False)
    train_idx=idx[~np.isin(idx,test_idx)]
    train_X,train_y=X[train_idx],y[train_idx]
    test_X,test_y=X[test_idx],y[test_idx]
    return train_X,train_y,test_X,test_y

#pandas
def split_train_test(X,y,frac=0.8):
    train_X=X.sample(frac=frac)
    test_X=X[~X.index.isin(train_X.index)]
    train_y,test_y=y[train_X.index],y[test_X.index]
    return train_X,train_y,test_X,test_y
复制代码

从新在房价数据集上测试：

X,y=buf.values[:,:-1],buf.values[:,-1]
train_X,train_y,test_X,test_y=split_train_test(X,y)
ref=scaling_ref(train_X)
train_X=minmax_scaling(train_X,ref)
test_X=minmax_scaling(test_X,ref)

model.train_ne(train_X,train_y)
print("train score: %f"%model.score(train_X,train_y))
print("test score: %f"%model.score(test_X,test_y))

train score: 0.753238
test score: 0.656748
复制代码

可见，以前大于0.7的score只是假象，在新数据上并无那么高。

2. 欠拟合与过拟合

在前一章节提到了多项式映射，当准备对一个数据集执行映射时，咱们其实并不知道将最高次设为多少才是合适的，为了追求更好的非线性近似，可能会把最高次设置得尽量大。而结合上一小节的交叉验证，会发现另外一个问题的存在：

x=(np.random.rand(10)-0.5)*6
x=np.sort(x)
y=np.sin(x)+np.cos(x)+np.random.randn(10)*0.4
X=x.reshape((-1,1))

x_range=np.linspace(-3,3,100)
X_range=x_range.reshape((-1,1))

test_x=(np.random.rand(10)-0.5)*6
test_y=np.sin(test_x)+np.cos(test_x)+np.random.randn(10)*0.4
test_X=test_x.reshape((-1,1))

def subplot(ax,h):
    X_=polynomial_mapping(X,h,False)
    test_X_=polynomial_mapping(test_X,h,False)
    X_range_=polynomial_mapping(X_range,h,False)
    model.train_ne(X_,y)
    ax.set_title('h=%d\ntrain score=%f\ntest score=%f'%(
            h,model.score(X_,y),model.score(test_X_,test_y)))
    ax.scatter(x,y)
    ax.plot(x_range,model.predict(X_range_),c='r')
    ax.set_xlabel('x')
    ax.set_ylabel('y')

fig=plt.figure(figsize=(13.5,3))
ax=fig.add_subplot(131)
subplot(ax,h=1)
ax=fig.add_subplot(132)
subplot(ax,h=3)
ax=fig.add_subplot(133)
subplot(ax,h=5)
复制代码

三张图的状况分别对应了：欠拟合、正常拟合、过拟合。

能够看到，一开始增长最高次时，训练集和测试集的score都在提高，而到某个临界点后，再增长，只有训练集的score提高，测试集的score反而降低了（因为数据集是随机生成的，若是没有看到上图的状况能够多执行几回）。至于这其中的缘由，从图中能够看出端倪，在时，曲线尝试穿过每个训练集的点，彷佛在强行记住每个点的位置，这种作法显然有些极端了，不具有容错性，时的曲线反而更合适。当训练一个模型时，应当指望它去寻找数据中的广泛规律，使得在面对新的数据时具有足够强的泛化能力，而不是仅仅记住训练数据。

避免过拟合有多种方式：第一种是增长训练集大小，数据点越多，曲线就越难穿过每个点，迫使模型去寻找一个更普适的函数，但数据集并非那么好收集的，每每条件受限无法作到这点；第二种是根据模型的复杂度在代价函数中增长惩罚项，以在准确率和复杂度之间寻求平衡；第三种是设置优化中止阈值，只适用于迭代优化，当代价变化量小于一个指定值时或是测试集评分持续多轮没有提高就提早结束优化，迭代优化到过拟合是有一个过程的，只要在此以前中止就好了。

注意，特征映射越复杂，越可能出现过拟合，越不容易欠拟合，在不执行映射的线性模式下，过拟合的可能性很低。

3. L2正则化

正则化是一种在代价函数中增长惩罚项以免过拟合的技巧，经常使用的有L1和L2两种正则化方式。L1正则化（或称L1范数）是对参数的绝对值求和，公式为 $L_1=\sum_{j=0}^m |\theta_j|$ ，可是因为绝对值是非连续可导函数，使用L1正则化后会致使没法使用梯度降低法优化，因此这里不采用；L2正则化（或称L2范数）是对参数的平方求和，公式为 $L2=\sum_{j=0}^m \theta_j^2$ 。

不管是L1仍是L2正则化，都是在把参数的总体大小做为惩罚项，并经过一个额外的 $\lambda$ 变量来控制该惩罚项的重要性，但为何参数的总体大小可以表示模型的复杂度呢？对此本人也没有深刻的研究，但简单些说，当一个模型的非线性表达能力很强时，它每每须要更多、更大的参数来完成这种表达，多项式映射模式下，当模型优先采用高次项时，不光高次项要采用相对较大的参数，低次项也要使用更大的参数去微调函数的形态。上一小节中若是将不一样下的最终正则化值输出就会发现，越大，正则化值也越大。

加入L2正则化后因为代价函数发生了改变，优化算法也要从新推导。代价函数： $J(\theta)=\frac{1}{2n}\sum_{i=0}^n(h_\theta (x_{i})-y_i)^2+\frac{\lambda}{2} \sum_{j=0}^m \theta_j^2$

梯度降低： $\theta^{new}=(1-\alpha\lambda)\theta-\frac{\alpha}{n} X^T(X\theta-y)$

正规方程：

\frac{1}{n}X^T(h_\theta (X)-y)+\lambda \theta=\vec{0}
\rightarrow
X^T(X \theta-y)+n\lambda \theta = \vec{0}

\rightarrow
(X^T X + n\lambda I)\theta = X^T y
\rightarrow
\theta=(X^T X + n\lambda I)^{-1}X^Ty

改进代码以下：

#修改线性回归类
class LinearRegression:
    
    def __init__(self,learning_rate=0.001,iter_max=100,batch_size=256,l2_lambda=0.001):
        self.learning_rate_=learning_rate
        self.iter_max_=iter_max
        self.batch_size_=batch_size
        self.l2_lambda_=l2_lambda

    #代价函数
    def cost(self,p_y,y):
        return 0.5*np.power(p_y-y,2).mean(axis=-1)\
            +0.5*self.l2_lambda_*np.power(self.theta_,2).sum()

    #梯度降低
    def grad_desc(self,X,y,theta,alpha):
        p_y=self.linear(X,theta)
        theta-=alpha*self.l2_lambda_*theta
        theta-=alpha*np.dot(X.T,p_y-y)/len(y)
        
    #正规方程
    def norm_equa(self,X,y):
        I=np.eye(X.shape[1])
        return np.dot(np.linalg.inv(np.dot(X.T,X)+self.l2_lambda_*I*len(y)),np.dot(X.T,y))
复制代码

从新进行测试：

model=LinearRegression(l2_lambda=0.0)
fig=plt.figure(figsize=(10,3))
ax=fig.add_subplot(131)
subplot(ax,h=1)
ax=fig.add_subplot(132)
subplot(ax,h=3)
ax=fig.add_subplot(133)
subplot(ax,h=5)

model=LinearRegression(l2_lambda=0.01)
fig=plt.figure(figsize=(10,3))
ax=fig.add_subplot(131)
subplot(ax,h=1)
ax=fig.add_subplot(132)
subplot(ax,h=3)
ax=fig.add_subplot(133)
subplot(ax,h=5)
复制代码

显然，L2正则化对避免过拟合是有帮助的，时的测试集准确率提高了。可是要注意的是，这不意味着高次+正则化就必定获得更好的结果，就此处的例子而言，获得的结果最好，模型还没有过拟合时，使用正则化反而可能会下降准确率。正则化的做用在于你能够更放心地假设数据存在复杂的非线性关系而不容易陷入过拟合的陷阱。

4. 提早终止

另外一种避免过拟合的方法，只能用于梯度降低这类迭代优化算法中，就是设置提早终止优化的条件。这也是一种节省训练时间的方式，能够减小模型在低效优化上耗费的时间。通常可针对两种指标设置结束条件：一是代价变化量，当该变化量小于必定值时咱们能够认为它已经接近极值点了因此梯度逐渐平缓；二是测试集的评分连续几轮迭代都没有提高。这两种状况再继续优化对模型表现的提高颇有限还容易陷入过拟合，提早结束反而能节省大量时间。

由于第一种方法须要频繁的调整阈值，我的更喜欢第二种，如下是实现代码：

#在线性回归类中修改
class LinearRegression:

    def __init__(self,learning_rate=0.001,iter_max=100,batch_size=256, l2_lambda=0.001,early_stop=True):
        self.learning_rate_=learning_rate
        self.iter_max_=iter_max
        self.batch_size_=batch_size
        self.l2_lambda_=l2_lambda
        self.early_stop_=early_stop

    #训练
    def train(self,X,y,test_X=None,test_y=None):
        #提早终止机制须要提供test_X,test_y
        if self.early_stop_&(type(test_X)==type(None))|(type(test_y)==type(None)):
            raise Exception('test_X and test_y can not be None for early stopping')
        #初始化
        X=self.fill_x0(X)
        self.init_theta(X.shape[1])
        #子集数量
        batches_n=math.ceil(len(y)/self.batch_size_)
        #优化历史
        self.theta_h_=[self.theta_.copy()]
        if self.early_stop_:
            not_improved=0
            theta_best=self.theta_.copy()
            score_best=self.score(test_X,test_y)
            self.cost_h_=[[self.cost(self.predict(X),y),
                           self.cost(self.predict(test_X),test_y)]]
        else:
            self.cost_h_=[self.cost(self.predict(X),y)]
        #迭代
        for i in range(self.iter_max_):
            #随机排序
            idx=np.random.permutation(len(y))
            X,y=X[idx],y[idx]
            #遍历子集
            for j in range(batches_n):
                X_=X[j*self.batch_size_:(j+1)*self.batch_size_]
                y_=y[j*self.batch_size_:(j+1)*self.batch_size_]
                self.grad_desc(
                    X_,y_,self.theta_,self.learning_rate_
                    )
            #记录当前参数和代价
            self.theta_h_.append(self.theta_.copy())
            #提早中止
            if self.early_stop_:
                self.cost_h_.append([self.cost(self.predict(X),y),
                                     self.cost(self.predict(test_X),test_y)])
                score_new=self.score(test_X,test_y)
                if score_new<=score_best:
                    not_improved+=1
                else:
                    not_improved=0
                    theta_best=self.theta_.copy()
                    score_best=score_new
                if not_improved>=10:
                    self.theta_=theta_best
                    print("Early stopping at iter %d!"%(i+1))
                    break
            else:
                self.cost_h_.append(self.cost(self.predict(X),y))
        self.theta_h_=np.array(self.theta_h_)
        self.cost_h_=np.array(self.cost_h_)
复制代码

最后再在房价数据集上作一次测试：

X,y=buf.values[:,:-1],buf.values[:,-1]
ref=scaling_ref(X)
X=minmax_scaling(X,ref)
train_X,train_y,test_X,test_y=split_train_test(X,y)

model=LinearRegression(learning_rate=0.3,iter_max=1000,
                       l2_lambda=0.001,early_stop=True)
model.train(train_X,train_y,test_X,test_y)

plt.title("train score: %f\ntest score: %f"%(
        model.score(train_X,train_y),
        model.score(test_X,test_y)
        ))
plt.plot(range(len(model.cost_h_)),model.cost_h_[:,0],label='train')
plt.plot(range(len(model.cost_h_)),model.cost_h_[:,1],label='test')
plt.xlabel('iter')
plt.ylabel('cost')

Early stopping at iter 164!
复制代码

因为房价数据的非线性关系不明显，因此就没有作特征映射。即便抛开防止过拟合的做用，提早终止在节省计算时间上的效果也是很重要的，没有过拟合的风险不表明该机制就没有用。