梯度下降法
梯度下降法用來求解目標函數的極值。這個極值是給定模型給定數據之后在參數空間中搜索找到的。迭代過程為：
可以看出，梯度下降法更新參數的方式為目標函數在當前參數取值下的梯度值，前面再加上一個步長控制參數alpha。梯度下降法通常用一個三維圖來展示，迭代過程就好像在不斷地下坡，最終到達坡底。為了更形象地理解，也為了和牛頓法比較，這里我用一個二維圖來表示：
懶得畫圖了直接用這個展示一下。在二維圖中，梯度就相當於凸函數切線的斜率，橫坐標就是每次迭代的參數，縱坐標是目標函數的取值。每次迭代的過程是這樣：
- 首先計算目標函數在當前參數值的斜率（梯度），然后乘以步長因子后帶入更新公式，如圖點所在位置（極值點右邊），此時斜率為正，那么更新參數后參數減小，更接近極小值對應的參數。
- 如果更新參數后，當前參數值仍然在極值點右邊，那么繼續上面更新，效果一樣。
- 如果更新參數后，當前參數值到了極值點的左邊，然后計算斜率會發現是負的，這樣經過再一次更新后就會又向着極值點的方向更新。
根據這個過程我們發現，每一步走的距離在極值點附近非常重要，如果走的步子過大，容易在極值點附近震盪而無法收斂。解決辦法：將alpha設定為隨着迭代次數而不斷減小的變量，但是也不能完全減為零。
梯度下降法實戰
我們來求解Ax=b，程序如下：
%matplotlib inline
import copy
import numpy as np
import matplotlib.pyplot as plt
A = np.array([[4,2],[1,3]])
b = np.array([[3],[2]])
loss = []
x0 = [copy.deepcopy([]), copy.deepcopy([])]
step = 0.01
x = [[0.],[4]]
while (A.dot(x)-b).T.dot(A.dot(x)-b) > 0.1:
dx = 2*A.T.dot(A).dot(x) - 2*A.T.dot(b)
x = x - step*dx
x0[0].append(x[0])
x0[1].append(x[1])
loss.append(np.squeeze((A.dot(x)-b).T.dot(A.dot(x)-b)))
line = np.linspace(0,len(x0[0])-1,len(x0[0]))
fig,(ax0,ax1)=plt.subplots(2,1, figsize=(9,6))
ax1.plot(line, x0[0])
ax1.plot(line, x0[1])
ax0.plot(line, loss)
ax1.plot(line, np.ones(len(x0[0]))*0.5)
plt.show
上圖表示loss函數的降低趨勢，下圖反映了解的收斂趨勢：
牛頓法
首先得明確，牛頓法是為了求解函數值為零的時候變量的取值問題的，具體地，
一階方法：
當要求解 f(θ)=0時，如果 f可導，那么可以通過迭代公式
來迭代求得最小值。通過一組圖來說明這個過程。
二階方法：
當應用於求解最大似然估計的值時，變成ℓ′(θ)=0的問題。這個與梯度下降不同，梯度下降的目的是直接求解目標函數極小值，而牛頓法則變相地通過求解目標函數一階導為零的參數值，進而求得目標函數最小值。那么迭代公式寫作：
當θ是向量時，牛頓法可以使用下面式子表示：
其中H叫做海森矩陣，其實就是目標函數對參數θ的二階導數。
通過比較牛頓法和梯度下降法的迭代公式，可以發現兩者及其相似。海森矩陣的逆就好比梯度下降法的學習率參數alpha。牛頓法收斂速度相比梯度下降法很快，而且由於海森矩陣的的逆在迭代中不斷減小，起到逐漸縮小步長的效果。
牛頓法的缺點就是計算海森矩陣的逆比較困難，消耗時間和計算資源。因此有了擬牛頓法。
實踐練習
%matplotlib inline
import copy
import numpy as np
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
A = np.array([[4,2],[1,3]])
b = np.array([[3],[2]])
X, Y = np.meshgrid(np.linspace(-10, 10, 200),np.linspace(-10, 10, 200))
def get_loss(X,Y):
Z = (A[0,0]**2+A[1,0]**2)*X**2 + (A[0,1]**2+A[1,1]**2)*Y**2 - 2*(A[0,0]*A[0,1]+A[1,0]*A[1,1])*X*Y - 2*(A[0,0]*b[0]+A[1,0]*b[1])*X \
- 2*(A[0,1]*b[0]+A[1,1]*b[1])*Y
return Z
Z = get_loss(X,Y)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X,Y,Z,rstride=5, cstride=5, alpha=0.3)
ax.contour(X,Y,Z, alpha=0.3)
cset = ax.contour3D(X, Y, Z, 10, zdir='z', offset=-100, cmap=cm.coolwarm)
# 梯度下降法
step = 0.01
x_g = np.array([[0],[-4]])
for i in range(100):
dx_g = 2*A.T.dot(A).dot(x_g[:,-1].reshape(2,1)) - 2*A.T.dot(b)
x_g = np.concatenate((x_g,x_g[:,-1].reshape(2,1) - step*dx_g),axis=1)
z_g = get_loss(x_g[0],x_g[1])
ax.plot(x_g[0],x_g[1],z_g,'ro')
# 牛頓法
# 黑塞矩陣
H = np.matrix([[A[0,0]**2+A[1,0]**2, A[0,0]*A[0,1]+A[1,0]*A[1,1]],
[A[0,0]*A[0,1]+A[1,0]*A[1,1], A[0,1]**2+A[1,1]**2]],dtype='float64') * 2
x_n = np.array([[0.],[-4.]])
for i in range(100):
dx_n = 2*A.T.dot(A).dot(x_n[:,-1].reshape(2,1)) - 2*A.T.dot(b)
x_n = np.concatenate((x_n, x_n[:,-1] - np.linalg.inv(H).dot(dx_n)),axis=1)
z_n = get_loss(np.squeeze(np.asarray(x_n[0])),np.squeeze(np.asarray(x_n[1])))
ax.plot(np.squeeze(np.asarray(x_n[0])),np.squeeze(np.asarray(x_n[1])),z_n,'go')
f, (ax1, ax2) = plt.subplots(1, 2)
ax1.plot(np.linspace(0,len(z_g)-1,len(z_g)),z_g,'b',label='grad')
ax1.plot(np.linspace(0,len(z_n)-1,len(z_n)),z_n,'r',label='Newton')
ax1.legend()
cs = ax2.contour(X,Y,Z)
ax2.clabel(cs, inline=1, fontsize=5)
ax2.plot(np.squeeze(np.asarray(x_n[0])),np.squeeze(np.asarray(x_n[1])),'b')
ax2.plot(x_g[0],x_g[1],'r')
程序學習：
嘗試了平面等高線圖和線標注的繪制，
cs = ax2.contour(X,Y,Z)
ax2.clabel(cs, inline=1, fontsize=5)
注意到數組拼接方法都是不破壞原數組，單純返回新數組的，且axis=0是行拼接（行數增加），axis=1是列拼接（列數增加），
x_n = np.concatenate((x_n, x_n[:,-1] - np.linalg.inv(H).dot(dx_n)),axis=1)
學習了numpy中的矩陣類型：np.matrix()，在牛頓法中我用的是matrix，在梯度下降法中我用的是array：
matrix是array的子類，特點是有且必須只是2維，matrix.I()可以求逆，和線代的求逆方法一致，所以繪圖時我不得不才用np.sequeeze(np.asarray())操作來降維，而由於x[:, -1]這種操作對array會自動降維（由兩行變為一行），所以要么使用matrix，要么切片后reshape(2,1)，總之不消停。
結果分析：
從兩張角度截一下圖，紅色是梯度下降藍色是牛頓法，可以看到，牛頓法收斂速度很快（外圍只有3個點），不過這是建立在黑塞矩陣的基礎上（需要求解目標函數的二階偏導數），這是牛頓法快速收斂的原因，也是牛頓法的瓶頸，而且這個瓶頸很直觀：我計算黑塞矩陣的numpy矩陣表達方法時的確費了挺大勁（其實原因更多是我渣... ...）
