1 最小二乘法概述
自從開始做畢設以來，發現自己無時無刻不在接觸最小二乘法。從求解線性透視圖中的消失點，m元n次函數的擬合，包括后來學到的神經網絡，其思想歸根結底全都是最小二乘法。
1-1 “多線→一點”視角與“多點→一線”視角
最小二乘法非常簡單，我把它分成兩種視角描述：
（1）已知多條近似交匯於同一個點的直線，想求解出一個近似交點：尋找到一個距離所有直線距離平方和最小的點，該點即最小二乘解；
（2）已知多個近似分布於同一直線上的點，想擬合出一個直線方程：設該直線方程為y=kx+b，調整參數k和b，使得所有點到該直線的距離平方之和最小，設此時滿足要求的k=k0，b=b0，則直線方程為y=k0x+b0。
1-2 思維拓展
這只是舉了兩個簡單的例子，其實在現實生活中我們可以利用最小二乘法解決更為復雜的問題。比方說有一個未知系數的二元二次函數f(x,y)=w0x^2+w1y^2+w2xy+w3x+w4y+w5，這里w0~w5為未知的參數，為了確定下來這些參數，將會給定一些樣本點(xi,yi,f(xi,yi))，然后通過調整這些參數，找到這樣一組w0~w5，使得這些所有的樣本點距離函數f(x,y)的距離平方之和最小。至於具體用何種方法來調整這些參數呢？有一種非常普遍的方法叫“梯度下降法”，它可以保證每一步調整參數，都使得f(x,y)朝比當前值更小的方向走，只要步長α選取合適，我們就可以達成這種目的。
而這里不得不提的就是神經網絡了。神經網絡其實就是不斷調整權值w和偏置b，來使得cost函數最小，從這個意義上來講它還是屬於最小二乘法。更為可愛的一點是，神經網絡的調參用到的仍是梯度下降法，其中最常用的當屬隨機梯度下降法。而后面偉大的bp算法，其實就是為了給梯度下降法做個鋪墊而已，bp算法的結果是cost函數對全部權值和全部偏置的偏導，而得知了這些偏導，對於各個權值w和偏置b該走向何方就指明了方向。
因此，最小二乘法在某種程度上無異於機器學習中基礎中的基礎，且具有相當重要的地位。至於上面所說的“梯度下降法”以及“利用最小二乘法求解二元二次函數的w0~w5”，我將會在后面的博客中進行更加詳細的探討。
2 scipy庫中的leastsq函數
當然，最小二乘法本身實現起來也是不難的，就如我們上面所說的不斷調整參數，然后令誤差函數Err不斷減小就行了。我們將在下一次博客中詳細說明如何利用梯度下降法來完成這個目標。
而在本篇博客中，我們介紹一個scipy庫中的函數，叫leastsq，它可以省去中間那些具體的求解步驟，只需要輸入一系列樣本點，給出待求函數的基本形狀（如我剛才所說，二元二次函數就是一種形狀——f(x,y)=w0x^2+w1y^2+w2xy+w3x+w4y+w5，在形狀給定后，我們只需要求解相應的系數w0~w6），即可得到相應的參數。至於中間到底是怎么求的，這一部分內容就像一個黑箱一樣。
2-1 函數形為y=kx+b
這一次我們給出函數形y=kx+b。這種情況下，待確定的參數只有兩個：k和b。
此時給出7個樣本點如下：
1 Xi=np.array([8.19,2.72,6.39,8.71,4.7,2.66,3.78]) 2 Yi=np.array([7.01,2.78,6.47,6.71,4.1,4.23,4.05])
則使用leastsq函數求解其擬合直線的代碼如下：
1 ###最小二乘法試驗### 2 import numpy as np 3 from scipy.optimize import leastsq 4 5 ###采樣點(Xi,Yi)### 6 Xi=np.array([8.19,2.72,6.39,8.71,4.7,2.66,3.78]) 7 Yi=np.array([7.01,2.78,6.47,6.71,4.1,4.23,4.05]) 8 9 ###需要擬合的函數func及誤差error### 10 def func(p,x): 11 k,b=p 12 return k*x+b 13 14 def error(p,x,y,s): 15 print s 16 return func(p,x)-y #x、y都是列表，故返回值也是個列表 17 18 #TEST 19 p0=[100,2] 20 #print( error(p0,Xi,Yi) ) 21 22 ###主函數從此開始### 23 s="Test the number of iteration" #試驗最小二乘法函數leastsq得調用幾次error函數才能找到使得均方誤差之和最小的k、b 24 Para=leastsq(error,p0,args=(Xi,Yi,s)) #把error函數中除了p以外的參數打包到args中 25 k,b=Para[0] 26 print"k=",k,'\n',"b=",b 27 28 ###繪圖，看擬合效果### 29 import matplotlib.pyplot as plt 30 31 plt.figure(figsize=(8,6)) 32 plt.scatter(Xi,Yi,color="red",label="Sample Point",linewidth=3) #畫樣本點 33 x=np.linspace(0,10,1000) 34 y=k*x+b 35 plt.plot(x,y,color="orange",label="Fitting Line",linewidth=2) #畫擬合直線 36 plt.legend() 37 plt.show()
我把里面需要注意的點提點如下：
1、p0里放的是k、b的初始值，這個值可以隨意指定。往后隨着迭代次數增加，k、b將會不斷變化，使得error函數的值越來越小。
2、func函數里指出了待擬合函數的函數形狀。
3、error函數為誤差函數，我們的目標就是不斷調整k和b使得error不斷減小。這里的error函數和神經網絡中常說的cost函數實際上是一回事，只不過這里更簡單些而已。
4、必須注意一點，傳入leastsq函數的參數可以有多個，但必須把參數的初始值p0和其它參數分開放。其它參數應打包到args中。
5、leastsq的返回值是一個tuple，它里面有兩個元素，第一個元素是k、b的求解結果，第二個元素我暫時也不知道是什么意思，先留下來。
其擬合效果圖如下：
2-2 函數形為y=ax^2+bx+c
這一次我們給出函數形y=ax^2+bx+c。這種情況下，待確定的參數有3個：a，b和c。
此時給出7個樣本點如下：
1 Xi=np.array([0,1,2,3,-1,-2,-3]) 2 Yi=np.array([-1.21,1.9,3.2,10.3,2.2,3.71,8.7])
這一次的代碼與2-1差不多，除了把待求參數再增加一個，換了一下訓練樣本，換了一下func中給出的函數形，幾乎沒有任何變化。
1 ###最小二乘法試驗### 2 import numpy as np 3 from scipy.optimize import leastsq 4 5 ###采樣點(Xi,Yi)### 6 Xi=np.array([0,1,2,3,-1,-2,-3]) 7 Yi=np.array([-1.21,1.9,3.2,10.3,2.2,3.71,8.7]) 8 9 ###需要擬合的函數func及誤差error### 10 def func(p,x): 11 a,b,c=p 12 return a*x**2+b*x+c 13 14 def error(p,x,y,s): 15 print s 16 return func(p,x)-y #x、y都是列表，故返回值也是個列表 17 18 #TEST 19 p0=[5,2,10] 20 #print( error(p0,Xi,Yi) ) 21 22 ###主函數從此開始### 23 s="Test the number of iteration" #試驗最小二乘法函數leastsq得調用幾次error函數才能找到使得均方誤差之和最小的a~c 24 Para=leastsq(error,p0,args=(Xi,Yi,s)) #把error函數中除了p以外的參數打包到args中 25 a,b,c=Para[0] 26 print"a=",a,'\n',"b=",b,"c=",c 27 28 ###繪圖，看擬合效果### 29 import matplotlib.pyplot as plt 30 31 plt.figure(figsize=(8,6)) 32 plt.scatter(Xi,Yi,color="red",label="Sample Point",linewidth=3) #畫樣本點 33 x=np.linspace(-5,5,1000) 34 y=a*x**2+b*x+c 35 plt.plot(x,y,color="orange",label="Fitting Curve",linewidth=2) #畫擬合曲線 36 plt.legend() 37 plt.show()
不過我們發現，它依舊能夠非常順利地解出待求的三個參數。其擬合情況如圖所示：
2-3 leastsq擬合y=kx+b可視化
本部分內容是建立在2-1代碼的基礎上，用Mayavi繪3D圖，以簡單地說明最小二乘法到底是怎么一回事。該部分知識用到了mgrid函數，具體是如何實施的請移步《Python閑談（一）mgrid慢放》。
step 1：創建一個k矩陣和b矩陣。在mgrid擴展后，有：
（1）k=[k1,k2,k3,...,kn]
mgrid(k)（朝右擴展）=
[k1,k1,k1,...,k1]
[k2,k2,k2,...,k2]
[k3,k3,k3,...,k3]
...
[kn,kn,kn,...,kn]
（2）b=[b1,b2,b3,...,bn]
mgrid(b)（朝下擴展）=
[b1,b2,b3,...,bn]
[b1,b2,b3,...,bn]
[b1,b2,b3,...,bn]
...
[b1,b2,b3,...,bn]
其中k矩陣和b矩陣等大（皆為n維向量,或者說1*n的矩陣），且這兩個矩陣里面的元素都非常密集。舉個例子以說明什么叫矩陣中的元素很密集：a是個矩陣，假設aij 為a矩陣中第i行第j列元素，則aij 和 a{i+1}j 的差值很小，aij 和 ai{j+1} 的差值也很小。也就是同一行或者同一列中相鄰的兩個元素的值非常接近。為什么要讓矩陣元素如此密集呢？因為我們的根本目的是用“密集的離散”來逼近“連續”，這里的思想就像微積分一樣。
而放在這里，就是ku和k{u+1}很接近，bv和b{v+1}也很接近。
step 2：令k矩陣和b矩陣中的元素按照其位置一一對應。對應后的結果為：
Combine_kb=
[(k1,b1),(k1,b2),(k1,b3)...,(k1,bn)]
[(k2,b1),(k2,b2),(k2,b3)...,(k2,bn)]
[(k3,b1),(k3,b2),(k3,b3)...,(k3,bn)]
...
[(kn,b1),(kn,b2),(kn,b3)...,(kn,bn)]
step 3：對矩陣中每一個(ku,bv)，我們分別求出該種情況下每一個訓練樣本點的誤差平方之和，即有：
Err{(ku,bv)}=∑{i=1~m}（(yi-（ku*xi+bv）)**2）
其中m為給定的訓練樣本點的個數。例如在這里：
1 Xi=np.array([8.19,2.72,6.39,8.71,4.7,2.66,3.78]) 2 Yi=np.array([7.01,2.78,6.47,6.71,4.1,4.23,4.05])
則有m=7。
什么意思呢？舉個例子，當i=1的時候，這個時候把(x1,y1)（=（8.19,7.01））代入（(yi-（ku*xi+bv）)**2）式子里面，由於此時已經鎖定(ku,bv)，因此式中所有的數都是常數，我們可以解出一個常數（(y1-（ku*x1+bv）)**2）。然后依次令i=2,3,4,...,7,可以分別求解出一個（(yi-（ku*xi+bv）)**2）值來，這7個（(yi-（ku*xi+bv）)**2）值加起來即Err{(ku,bv)}。
注意了，最終我們算出的那個Err{(ku,bv)}將會存放到ku、bv對應的那個位置，比方說u=3,v=2：
mgrid(k)=
[k1,k1,k1,...,k1]
[k2,k2,k2,...,k2]
[k3,k3,k3,...,k3]
...
[kn,kn,kn,...,kn]
mgrid(b)=
[b1,b2,b3,...,bn]
[b1,b2,b3,...,bn]
[b1,b2,b3,...,bn]
...
[b1,b2,b3,...,bn]
則剛才算出來的Err{(k3,b2)}應該放在這個位置：
Err=
[Err11,Err12,Err13,...,Err1n]
[Err21,Err22,Err23,...,Err2n]
[Err31,Err32,Err33,...,Err3n]
...
[Errn1,Errn2,Errn3,...,Errnn]
如此這般對於每一對(ku,bv)都這樣算，則上方的Err矩陣中每一個元素的值都可以算出來；將計算出的結果正確地放在Err矩陣中對應位置，即得到Err矩陣。
step 4：繪制曲面。
截至目前我們已經得到了兩個重要矩陣Combine_kb和Err，其中Combine_kb提供點的x、y軸坐標，Err矩陣提供點的z軸坐標。
Combine_kb=
[(k1,b1),(k1,b2),(k1,b3)...,(k1,bn)]
[(k2,b1),(k2,b2),(k2,b3)...,(k2,bn)]
[(k3,b1),(k3,b2),(k3,b3)...,(k3,bn)]
...
[(kn,b1),(kn,b2),(kn,b3)...,(kn,bn)]
Err=
[Err11,Err12,Err13,...,Err1n]
[Err21,Err22,Err23,...,Err2n]
[Err31,Err32,Err33,...,Err3n]
...
[Errn1,Errn2,Errn3,...,Errnn]
我們再將這兩個矩陣合並一下得到Combine_kbErr矩陣：
Combine_kbErr=
[(k1,b1,Err11),(k1,b2,Err12),(k1,b3,Err13)...,(k1,bn,Err1n)]
[(k2,b1,Err21),(k2,b2,Err22),(k2,b3,Err23)...,(k2,bn,Err2n)]
[(k3,b1,Err31),(k3,b2,Err32),(k3,b3,Err33)...,(k3,bn,Err3n)]
...
[(kn,b1,Errn1),(kn,b2,Errn2),(kn,b3,Errn3)...,(kn,bn,Errnn)]
在三維空間直角坐標系下繪制出Combine_kbErr中的每一個點，然后將這些點與其各自相鄰的點連起來，則得到我們想要的Err(k,b)函數曲面。
step 5：本部分代碼如下：
1 """part 2""" 2 ###定義一個函數，用於計算在k、b已知時∑((yi-(k*xi+b))**2)### 3 def S(k,b): 4 ErrorArray=np.zeros(k.shape) #k的shape事實上同時也是b的shape 5 for x,y in zip(Xi,Yi): #zip(Xi,Yi)=[(8.19,7.01),(2.72,2.78),...,(3.78,4.05)] 6 ErrorArray+=(y-(k*x+b))**2 7 return ErrorArray 8 9 ###繪制ErrorArray+最低點### 10 from enthought.mayavi import mlab 11 12 #畫整個Error曲面 13 k,b=np.mgrid[k0-1:k0+1:10j,b0-1:b0+1:10j] 14 Err=S(k,b) 15 face=mlab.surf(k,b,Err/500.0,warp_scale=1) 16 mlab.axes(xlabel='k',ylabel='b',zlabel='Error') 17 mlab.outline(face) 18 19 #畫最低點（即k，b所在處） 20 MinErr=S(k0,b0) 21 mlab.points3d(k0,b0,MinErr/500.0,scale_factor=0.1,color=(0.5,0.5,0.5)) #scale_factor用來指定點的大小 22 mlab.show()
對要點說明如下：
1、為了讓最小二乘法求解的結果出現在繪制曲面的范圍內，我們以最終leastsq求得的k0、b0為中心創建k向量和b向量。
2、傳入S函數的是k向量和b向量mgrid后的結果。
3、S函數中的ErrorArray+=(y-(k*x+b))**2 操作里，k、b皆為矩陣（是k、b向量mgrid后的結果），而x、y皆為常數，故這里的操作實際上是對矩陣的操作。這個ErrorArray就是上面我說的Err矩陣。
4、在繪圖時之所以對Err除以500，是因為Err和k、b的差距不是一般的大，直接繪圖會導致什么都看不出來。舉一個最簡單的例子就是比如我們要畫個二維直角坐標系下的圖，x的取值范圍是0~1，y的取值范圍是0~1000，而兩個坐標軸卻都按一個單位△x=△y=0.1來畫，想想看結果會成什么樣子？
這里也是同樣的道理，於是得給Err除以一個大數才能讓圖像正常顯示。
其實matplotlib畫三維坐標系下的圖會幫你調整到合適，只有Mayavi才會出現這種情況，反正注意一下比例問題就好了。
5、該程序除過繪制Err曲面外，還把(k0,b0)也畫出來了，見灰色小球。
step 6：整個程序的全部代碼如下，其中part1與2-1的代碼是完全一樣的。
1 ###【最小二乘法試驗】### 2 import numpy as np 3 from scipy.optimize import leastsq 4 5 ###采樣點(Xi,Yi)### 6 Xi=np.array([8.19,2.72,6.39,8.71,4.7,2.66,3.78]) 7 Yi=np.array([7.01,2.78,6.47,6.71,4.1,4.23,4.05]) 8 9 """part 1""" 10 ###需要擬合的函數func及誤差error### 11 def func(p,x): 12 k,b=p 13 return k*x+b 14 15 def error(p,x,y): 16 return func(p,x)-y #x、y都是列表，故返回值也是個列表 17 18 p0=[1,2] 19 20 ###最小二乘法求k0、b0### 21 Para=leastsq(error,p0,args=(Xi,Yi)) #把error函數中除了p以外的參數打包到args中 22 k0,b0=Para[0] 23 print"k0=",k0,'\n',"b0=",b0 24 25 """part 2""" 26 ###定義一個函數，用於計算在k、b已知時，∑((yi-(k*xi+b))**2)### 27 def S(k,b): 28 ErrorArray=np.zeros(k.shape) #k的shape事實上同時也是b的shape 29 for x,y in zip(Xi,Yi): #zip(Xi,Yi)=[(8.19,7.01),(2.72,2.78),...,(3.78,4.05)] 30 ErrorArray+=(y-(k*x+b))**2 31 return ErrorArray 32 33 ###繪制ErrorArray+最低點### 34 from enthought.mayavi import mlab 35 36 #畫整個Error曲面 37 k,b=np.mgrid[k0-1:k0+1:10j,b0-1:b0+1:10j] 38 Err=S(k,b) 39 face=mlab.surf(k,b,Err/500.0,warp_scale=1) 40 mlab.axes(xlabel='k',ylabel='b',zlabel='Error') 41 mlab.outline(face) 42 43 #畫最低點（即k，b所在處） 44 MinErr=S(k0,b0) 45 mlab.points3d(k0,b0,MinErr/500.0,scale_factor=0.1,color=(0.5,0.5,0.5)) #scale_factor用來指定點的大小 46 mlab.show()
3 結語
本次博客給出了最小二乘法的Python實現方法，它用到了scipy庫中的leastsq函數。在上面我們給出了兩個實例，分別實現了對一元一次函數的擬合和一元二次函數的擬合，而事實上，對於函數並不一定得是一元函數，對於更多元的函數也同樣能夠利用最小二乘法完成擬合工作，不過隨着元和次的增加，待求參數也就越來越多了，比方說二元二次函數就有6個待求參數w0~w6。
然為了更好地理解神經網絡的訓練算法，並不建議直接使用leastsq函數完成對未知參數的求解，因此在以后的博客中我會詳細說明如何利用梯度下降法來求解誤差函數的最小值。
4 下面要寫的博客
1、梯度下降法（什么是梯度下降法，如何使用梯度下降法求一元二次函數最小值，如何使用梯度下降法求二元二次函數的最小值）
2、最小二乘法擬合二元二次函數（即求解w0~w5，相當於對梯度下降法的一個應用）
以上兩個內容我會放在同一個博客中。
2016.5.14
by 悠望南山