[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!
複製代碼