最佳化演算法推導&手寫實現——牛頓法和擬牛頓法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矩陣的訓練也比二階導數矩陣的訓練來得方便與快速。
以上就是本次的全部內容,謝謝閱讀。(今天先寫到這,寫的好累,擬牛頓法的其他演算法我們下次再接著推導與實現。)
全部程式碼可前往以下地址下載:
往期回顧:
機器學習演算法推導&實現——邏輯斯蒂迴歸
機器學習演算法推導&實現——線性迴歸
學無止境,歡迎關注筆者公眾號、知乎號,互相學習!