最佳化演算法推導&手寫實現——牛頓法和擬牛頓法1
在上一篇文章《機器學習演算法推導&實現——邏輯斯蒂迴歸》中，我們分別使用了
梯度下降法
和
牛頓法
來求解對數似然函式的最最佳化問題，從例項中我們也能發現
牛頓法的收斂速度遠快於梯度下降法
，但是在上文中，直接使用了牛頓法的結果，並沒有進行相應推導，故本文一方面是
補上牛頓法的推導
，另一方面是
展開討論下擬牛頓法。
牛頓法推導
牛頓法的推導得先從泰勒公式說起：
泰勒公式
假設現在函式f（x）迭代了k次的值為X（k），則在X（k）上進行二階泰勒展開可近似得到以下公式：
泰勒二階展開
我們要求得f（x）的極小值，則必要條件是f（x）在極值點處的一階導數為0，即：
因為我們把每輪迭代求得的滿足目標函式極小值的x作為下一輪迭代的值，因此我們可以假設第k+1輪的值就是最優解：
代入二階泰勒展開並求導可得：
令：
一階導數
二階導數
可得最終的最佳化公式為：
牛頓法迭代
雖然根據推導，我們已經知道了牛頓法的迭代方法，但是在實際應用過程中，我們會發現海塞矩陣的逆矩陣往往計算比較複雜，於是又有了擬牛頓法來簡化這一過程。
擬牛頓法的原理
在擬牛頓法中，我們考慮優化出一個n階矩陣D來代替海塞矩陣的逆矩陣，首先我們得先來看看替代矩陣D要滿足什麼條件：
首先根據二階泰勒展開，我們令：
代入f(x)，並求導，可得：
泰勒二階展開
再令：
最終可得矩陣D需滿足的條件：
替代矩陣Dk需滿足的條件
另外在每次迭代中可以更新矩陣D，故可以假設以下更新條件：
由此我們可以發現海塞矩陣逆矩陣的近似矩陣D（x）的選擇條件比較靈活，可以有多種具體的實現方法。
1、DFP演算法（Davidon-Fletcher-Powell）推導
在DFP演算法中，我們假設：
*在這裡我們知道凸函式的二階導數海塞矩陣（Hessian）是對稱矩陣，因此我們假設想要替代的D， P， Q矩陣也是對稱矩陣，即矩陣的轉置等於矩陣本身。
兩邊乘以
可得：
為了滿足上式要求，我們不妨令：
我們
先來求P矩陣
，最簡單的就是令：
接著兩邊乘以
可得：
求得
，並代入P矩陣，可求得：
同樣的方法我們
求解下Q矩陣
：
最終，
DFP演算法的迭代公式為：
DFP演算法替代海塞矩陣逆矩陣的迭代公式
python實現DFP演算法
我們繼承上一篇文章邏輯迴歸演算法的類，並且增加擬牛頓法的最佳化演算法。
我們先看下演算法的主函式部分，在繼承邏輯迴歸類的基礎上主要有2點改動：
一是把"method"引數提前，在例項化的時候就要求定位好最佳化演算法；二是重寫了train方法，衍生出其他擬牛頓演算法。
from
ML_LogisticRegression
import
LogisticRegressionSelf
as
logit
import
numpy
as
np
class
allLogitMethod
（
logit
）：
def
__init__
（
self
，
method
）：
“”“init class”“”
super
（）
。
__init__
（）
self
。
method
=
method
self
。
graList
=
［］
#記錄梯度模長的列表
#重寫下train方法
def
train
（
self
，
X
，
y
，
n_iters
=
1000
，
learning_rate
=
0。01
）：
“”“fit model”“”
if
X
。
ndim
<
2
：
raise
ValueError
（
“X must be 2D array-like！”
）
self
。
trainSet
=
X
self
。
label
=
y
if
self
。
method
。
lower
（）
==
“gradient”
：
self
。
_LogisticRegressionSelf__train_gradient
（
n_iters
，
learning_rate
）
elif
self
。
method
。
lower
（）
==
“newton”
：
self
。
_LogisticRegressionSelf__train_newton
（
n_iters
）
elif
self
。
method
。
lower
（）
==
“dfp”
：
self
。
__train_dfp
（
n_iters
，
learning_rate
）
else
：
raise
ValueError
（
“method value not found！”
）
return
在寫演算法主體前，我們先看下DFP的
簡易最佳化過程：
1。 初始化引數W0，初始化替代矩陣Dk0，計算初始梯度Gk0，迭代次數k；
2。 While （ k < n_iters ） and （ ||Gk0|| > e ）：
- 2。1 計算出W0需要最佳化的方向Pk = Dk*Gk0
-
2。2 更新引數W1如下：
-
for n in N：
- 嘗試用一維搜尋方法改變Pk的學習率，
- 求得最小似然值，
- 求出W1和deltaW（W1 - W0）
- 2。3 更新梯度Gk1，並求deltaG（Gk1 - Gk0）
- 2。4 更新DK1，利用上述推導的DFP的迭代公式
- 2。5 重新賦值W0，Dk0，Gk0
- 2。6 k += 1
3。 end。
我們可以大致看到DFP演算法的主體部分有外迴圈和內迴圈兩個部分組成，所以我們也是將內外迴圈分開來寫。
我們先看
內迴圈部分
。內迴圈是已經知道了引數W的最佳化方向的基礎上，想要找到最合適的學習率，使得更新後的W求得的似然值最小。實際有很多方法，在這裡為了簡單，筆者只是對學習率進行了i次方的操作。具體函式如下：
#一維搜尋法求出最優lambdak，更新W後，使得似然值最小
def
__updateW
（
self
，
X
，
Y
，
lambdak
，
W0
，
Pk
）：
“”“此處對lambdak的處理僅簡單用1~i次方來逐步變小，以選取到最小似然值的lambdak”“”
min_LLvalue
=
np
。
inf
W1
=
np
。
zeros
（
W0
。
shape
）
for
i
in
range
（
10
）：
Wi
=
W0
-
（
lambdak
**
i
）
*
Pk
Ypreprob
，
LLvalue
=
self
。
PVandLLV
（
X
，
Y
，
Wi
）
if
LLvalue
<
min_LLvalue
：
min_LLvalue
=
LLvalue
W1
=
np
。
copy
（
Wi
）
deltaW
=
-
（
lambdak
**
i
）
*
Pk
bestYpreprob
=
Ypreprob
return
W1
，
deltaW
，
min_LLvalue
，
bestYpreprob
再看
外迴圈部分
。
#新增擬牛頓法-DFP最佳化演算法
def
__train_dfp
（
self
，
n_iters
，
learning_rate
）：
“”“Quasi-Newton Method DFP（Davidon-Fletcher-Powell）”“”
n_samples
，
n_features
=
self
。
trainSet
。
shape
X
=
self
。
trainSet
y
=
self
。
label
#合併w和b，在X尾部新增一列全是1的特徵
X2
=
np
。
hstack
（（
X
，
np
。
ones
（（
n_samples
，
1
））））
#將y轉置變為（n_samples，1）的矩陣
Y
=
np
。
expand_dims
（
y
，
axis
=
1
）
#初始化特徵係數W，初始化替代對稱矩陣
W
=
np
。
zeros
（（
1
，
n_features
+
1
））
Dk0
=
np
。
eye
（
n_features
+
1
）
#計算初始的預測值、似然值，並記錄似然值
Ypreprob
，
LL0
=
self
。
PVandLLV
（
X2
，
Y
，
W
）
self
。
llList
。
append
（
LL0
）
#根據初始的預測值計算初始梯度，並記錄梯度的模長
Gk0
=
self
。
_LogisticRegressionSelf__calGradient
（
X2
，
Y
，
Ypreprob
）
graLength
=
np
。
linalg
。
norm
（
Gk0
）
self
。
graList
。
append
（
graLength
）
#初始化迭代次數
k
=
0
while
（
k
<
n_iters
）
and
（
graLength
>
self
。
tol
）：
#計算最佳化方向的值Pk=Gk0。Dk0
Pk
=
np
。
dot
（
Gk0
，
Dk0
）
#一維搜尋更新引數，並儲存求得的最小似然值
W
，
deltaW
，
min_LLvalue
，
Ypreprob
=
self
。
__updateW
（
X2
，
Y
，
learning_rate
，
W
，
Pk
）
self
。
llList
。
append
（
min_LLvalue
）
#更新梯度Gk和deltaG，同時求得梯度的模長和更新前後的模長差值
Gk1
=
self
。
_LogisticRegressionSelf__calGradient
（
X2
，
Y
，
Ypreprob
）
graLength
=
np
。
linalg
。
norm
（
Gk1
）
self
。
graList
。
append
（
graLength
）
deltaG
=
Gk1
-
Gk0
Gk0
=
Gk1
#更新替代矩陣Dk
Dk1
=
Dk0
+
np
。
dot
（
deltaW
。
T
，
deltaW
）
/
np
。
dot
（
deltaW
，
deltaG
。
T
）
-
\
np
。
dot
（
np
。
dot
（
np
。
dot
（
Dk0
，
deltaG
。
T
），
deltaG
），
Dk0
）
/
np
。
dot
（
np
。
dot
（
deltaG
，
Dk0
），
deltaG
。
T
）
Dk0
=
Dk1
k
+=
1
self
。
n_iters
=
k
self
。
w
=
W
。
flatten
（）［：
-
1
］
self
。
b
=
W
。
flatten
（）［
-
1
］
Ypre
=
np
。
argmax
（
np
。
column_stack
（（
1
-
Ypreprob
，
Ypreprob
）），
axis
=
1
）
self
。
accurancy
=
sum
（
Ypre
==
y
）
/
n_samples
（
“第{}次停止迭代，梯度模長為{}，似然值為{}，準確率為{}”
。
format
（
self
。
n_iters
，
self
。
graList
［
-
1
］，
self
。
llList
［
-
1
］，
self
。
accurancy
））
（
“w：{}；
\n
b：{}”
。
format
（
self
。
w
，
self
。
b
））
return
在這裡還有一點小小的提示：
筆者在實測過程中，發現DFP演算法經常會導致求解似然值和P（y=1|X）的時候報值溢位的錯誤。主要是python的numpy。exp（wx）中的wx過大導致的，比如看下圖的報錯：
值溢位報錯資訊
所以在這裡我們對求解似然值和P（y=1|X）的公式進行了一些調整。
當wx為負數時，使用原公式；當wx為正數時，使用轉換公式替換。
其實質是一致的。
原公式和替換公式
最後我們用擬牛頓法-DFP演算法和牛頓法實測比較下。
牛頓法：
迭代了7次，迭代時長0。18s，似然值176。81，準確率0。94，模型引數如下圖所示：
if
__name__
==
“__main__”
：
import
matplotlib。pyplot
as
plt
from
sklearn。datasets
import
make_classification
import
time
X
，
y
=
make_classification
（
n_samples
=
1000
，
n_features
=
4
）
#1、自編的牛頓法進行擬合
time_nt
=
time
。
time
（）
logit_nt
=
allLogitMethod
（
“newton”
）
logit_nt
。
train
（
X
，
y
，
n_iters
=
100
）
（
“迭代時長：”
，
time
。
time
（）
-
time_nt
）
plt
。
plot
（
range
（
logit_nt
。
n_iters
+
1
），
logit_nt
。
llList
）
plt
。
show
（）
牛頓法迭代結果
擬牛頓法-DFP演算法：
迭代了18次，迭代時長0。038s，似然值176。81，準確率0。94，模型引數如下圖所示：
#3、自編的擬牛頓法-DFP演算法進行擬合
time_dfp
=
time
。
time
（）
logit_dfp
=
allLogitMethod
（
“DFP”
）
logit_dfp
。
train
（
X
，
y
，
n_iters
=
20
，
learning_rate
=
0。5
）
（
“迭代時長：”
，
time
。
time
（）
-
time_dfp
）
fig
=
plt
。
figure
（
figsize
=
（
10
，
6
））
ax1
=
fig
。
add_subplot
（
1
，
2
，
1
）
ax1
。
plot
（
range
（
logit_dfp
。
n_iters
+
1
），
logit_dfp
。
llList
）
ax2
=
fig
。
add_subplot
（
1
，
2
，
2
）
ax2
。
plot
（
range
（
logit_dfp
。
n_iters
+
1
），
logit_dfp
。
graList
）
plt
。
show
（）
擬牛頓法迭代結果
經過兩種方法的對比，我們發現兩者最終訓練出來的
模型引數基本一致，似然值也下降到同一水平
；
DFP雖然迭代18次比牛頓法多，但是訓練總時間只有0.038s，遠小於牛頓法的訓練總時間。
可見，DFP演算法中的Dk矩陣完全可以逼近於牛頓法中的海塞矩陣的逆矩陣，而且DFP演算法中的Dk矩陣的訓練也比二階導數矩陣的訓練來得方便與快速。
以上就是本次的全部內容，謝謝閱讀。（今天先寫到這，寫的好累，擬牛頓法的其他演算法我們下次再接著推導與實現。）
全部程式碼可前往以下地址下載：
往期回顧：
機器學習演算法推導&實現——邏輯斯蒂迴歸
機器學習演算法推導&實現——線性迴歸
學無止境，歡迎關注筆者公眾號、知乎號，互相學習！
