机器学习笔记之无约束优化问题——基于共轭方向法与Wolfe准则优化方法的Python示例
- 引言
- 小插曲：画图——非标准二次型的等值线
- 算法在图像中的表示
- 基于精确搜索的共轭梯度法
- 基于Wolfe准则的共轭梯度法
- 附：共轭梯度法完整代码
引言
本节使用 Python \text{Python} Python对共轭梯度法的精确搜索与非精确搜索进行示例。
本人数学水平与代码水平有限，欢迎小伙伴一起讨论~关联文章：
小插曲：画图——非标准二次型的等值线
非标准二次型——这意味着：对应函数
f
(
x
)
=
x
T
Q
x
\begin{aligned}f(x) = x^T \mathcal Q x\end{aligned}
f(x)=xTQx中的正定矩阵
Q
\mathcal Q
Q不是对角阵。本节以凸二次函数：
f
(
x
)
=
x
T
Q
x
Q
=
(
1
1
1
2
)
;
x
∈
R
2
f(x) = x^T \mathcal Q x \quad \mathcal Q = \begin{pmatrix}1 & 1 \\ 1 & 2\end{pmatrix};x \in \mathbb R^2
f(x)=xTQxQ=(1112);x∈R2
的等值线为例，使用
Python
\text{Python}
Python做出基于
Wolfe Condition
\text{Wolfe Condition}
Wolfe Condition与精确搜索的共轭梯度法效果。代码见文章最下方，后面不再赘述。如果使用二元函数进行表示，可以表示为如下形式：
f
(
x
1
,
x
2
)
=
(
x
1
x
2
)
(
1
1
1
2
)
(
x
1
x
2
)
=
(
x
1
+
x
2
x
1
+
2
x
2
)
(
x
1
x
2
)
=
x
1
2
+
2
x
1
x
2
+
2
x
2
2
\begin{aligned} f(x_1,x_2) & = (x_1 \text{ } x_2)\begin{pmatrix} 1 & 1 \\ 1 & 2 \end{pmatrix}\begin{pmatrix} x_1 \\ x_2 \end{pmatrix} \\ & = (x_1 + x_2 \quad x_1 + 2x_2) \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} \\ & = x_1^2 + 2 x_1 x_2 + 2x_2^2 \end{aligned}
f(x1,x2)=(x1 x2)(1112)(x1x2)=(x1+x2x1+2x2)(x1x2)=x12+2x1x2+2x22
很明显：该函数中不仅包含平方项，并且包含交叉项。因而无法将
x
1
,
x
2
x_1,x_2
x1,x2进行独立表示。对于等值线函数表示如下：
基于
C
>
0
\mathcal C>0
C>0的不同取值，可以得到相应大小的等值线结果。
x
1
2
+
2
x
1
x
2
+
2
x
2
2
=
C
x_1^2 + 2 x_1 x_2 + 2x_2^2 = \mathcal C
x12+2x1x2+2x22=C
针对上述问题，这里的思路是：在给定等值线
C
\mathcal C
C的条件下，确定
x
x
x的范围。判断
x
1
x1
x1是否为边缘的条件是：将
x
1
,
C
x_1,\mathcal C
x1,C均视作常数，此时上述函数就是关于
x
2
x_2
x2的一元二次方程，只需要求解根判别公式：
Δ
=
b
2
−
4
a
c
≜
0
\Delta = b^2 - 4ac \triangleq 0
Δ=b2−4ac≜0即可找到该范围：
因为Δ ≜ 0 \Delta \triangleq 0 Δ≜0说明一元二次方程有唯一解,意味着随机变量 x 1 x_1 x1的正交基在该位置与函数图像相切,见下示例图:
其中a = 2 , b = 2 x 1 , c = x 1 2 − C a = 2,b=2 x_1,c = x_1^2 - \mathcal C a=2,b=2x1,c=x12−C,带入即可。
0 ≜ Δ = b 2 − 4 a c = 4 x 1 2 − 8 ( x 1 2 − C ) ⇒ x = ± 2 C C > 0 \begin{aligned} 0 \triangleq \Delta & = b^2 - 4 ac \\ & = 4 x_1^2 - 8(x_1^2 - \mathcal C) \\ & \Rightarrow x = \pm \sqrt{2 \mathcal C} \quad \mathcal C > 0 \end{aligned} 0≜Δ=b2−4ac=4x12−8(x12−C)⇒x=±2CC>0
基于该范围，范围边缘的
±
2
C
\pm \sqrt{2 \mathcal C}
±2C只有唯一解，范围内的其他点均对应两个不相同的解。使用求根公式：
−
b
±
b
2
−
4
a
c
2
a
\begin{aligned}\frac{-b \pm \sqrt{b^2 - 4ac}}{2a}\end{aligned}
2a−b±b2−4ac将大值与小值分开作图：
一次画半个椭圆~
对应代码表示如下：
import math
import matplotlib.pyplot as plt
def CalxLimits(C):
# 根判别式
return math.sqrt(2 * C),-1 * math.sqrt(2 * C)
# def f(x,y):
# 没有用到;只是描述一下这个非标准二次型函数。
# return 2 * (y ** 2) + 2 * x * y + (x ** 2)
def GetAnalytical(x,C):
# 求根公式
return 0.25 * (-2 * x + math.sqrt(8 * C - (4 * (x ** 2)))),0.25 * (-2 * x - math.sqrt(8 * C - (4 * (x ** 2))))
def DrawContourPicture(CList):
for C in CList:
UpperPointList,LowerPointList = list(),list()
Upper,Lower = CalxLimits(C)
xList = list(np.linspace(Lower,Upper,10000))
for x in xList:
Upperx,Lowerx = GetAnalytical(x,C)
UpperPointList.append(Upperx)
LowerPointList.append(Lowerx)
plt.plot(xList, UpperPointList,'--',c="tab:blue")
plt.plot(xList, LowerPointList,'--',c="tab:blue")
plt.show()
if __name__ == '__main__':
CList = [0.01, 0.1, 0.5, 2, 4.5]
DrawContourPicture(CList)
对应图像效果表示如下：
算法在图像中的表示
基于精确搜索的共轭梯度法
由于目标函数
f
(
⋅
)
f(\cdot)
f(⋅)是凸二次函数，那么该函数一定存在全局最优解。因而可以使用基于精确搜索的共轭梯度法获取最优解。回顾无约束优化问题——共轭梯度法，其算法步骤表示如下：
初始化操作：
- 给定初始点 x 0 = ( 3 2 ) T x_0 = (3 \quad 2)^T x0=(32)T，记 d 0 = − ∇ f ( x 0 ) d_0 = - \nabla f(x_0) d0=−∇f(x0)；设置阈值 ϵ = 0.05 ; k = 0 \epsilon = 0.05;k=0 ϵ=0.05;k=0
算法过程：
- 首先判断范数
∥
∇
f
(
x
k
)
∥
≤
ϵ
\|\nabla f(x_k)\| \leq \epsilon
∥∇f(xk)∥≤ϵ是否成立
?
?
?若成立，则算法终止；
个人理解：这里利用范数来侧面描述梯度向量∇ f ( x k ) \nabla f(x_k) ∇f(xk)的大小。当∥ ∇ f ( x k ) ∥ ⇒ 0 \|\nabla f(x_k)\| \Rightarrow 0 ∥∇f(xk)∥⇒0意味着向量∇ f ( x k ) \nabla f(x_k) ∇f(xk)趋于零向量。 - 计算当前迭代步骤的最优步长
α
k
\alpha_k
αk：
需要注意一下：上面链接文章中对目标函数的定义为f ( x ) = 1 2 x T Q x + C T x \begin{aligned}f(x) =\frac{1}{2}x^T \mathcal Q x + \mathcal C^T x\end{aligned} f(x)=21xTQx+CTx,而系数1 2 \begin{aligned}\frac{1}{2}\end{aligned} 21只是为了方便求导。在本节中的f ( x ) = x T Q x f(x) = x^T \mathcal Q x f(x)=xTQx没有系数，因而需要在相应α k \alpha_k αk化简结果中填上一个系数1 2 \begin{aligned}\frac{1}{2}\end{aligned} 21。
α k = − [ ∇ f ( x k ) ] T d k 2 ( d k ) T Q d k \alpha_k = - \frac{[\nabla f(x_k)]^T d_k}{2(d_k)^T \mathcal Q d_k} αk=−2(dk)TQdk[∇f(xk)]Tdk - 计算新位置点：
x
k
+
1
=
x
k
+
α
k
⋅
d
k
x_{k+1} = x_k + \alpha_k \cdot d_k
xk+1=xk+αk⋅dk，并计算共轭方向
d
k
+
1
d_{k+1}
dk+1：
d k + 1 = − ∇ f ( x k + 1 ) + β k ⋅ d k , β k = [ ∇ f ( x k + 1 ) ] T Q d k ( d k ) T Q d k d_{k+1} = - \nabla f(x_{k+1}) + \beta_k \cdot d_k,\beta_k = \frac{[\nabla f(x_{k+1})]^T \mathcal Q d_k}{(d_k)^T \mathcal Q d_k} dk+1=−∇f(xk+1)+βk⋅dk,βk=(dk)TQdk[∇f(xk+1)]TQdk - 令 k = k + 1 k = k+1 k=k+1，转步骤 1 1 1重新进行判断。
相应代码表示如下：
def ConjugateGradient():
def f(PointInput, Q):
# 二次型
return np.dot(np.dot(PointInput, Q), PointInput)
def Nablaf(PointInput, Q):
# 二次型求导
return 2 * np.dot(Q, PointInput)
def GetNorm(ArrayInput):
return np.linalg.norm(ArrayInput)
Epsilon = 0.05
InitPoint = np.array([3., 2.])
Q = np.array([[1., 1.], [1., 2.]])
PointList = list()
PointList.append(InitPoint)
ConjugateStart = -1 * Nablaf(InitPoint, Q)
while True:
if GetNorm(Nablaf(InitPoint, Q)) <= Epsilon:
break
else:
alpha = -0.5 * (np.dot(Nablaf(InitPoint, Q), ConjugateStart) / np.dot(np.dot(ConjugateStart, Q),ConjugateStart))
NewPoint = InitPoint + alpha * ConjugateStart
Beta = np.dot(np.dot(Nablaf(NewPoint, Q), Q), ConjugateStart) / np.dot(np.dot(ConjugateStart, Q),ConjugateStart)
NewConjugate = -1 * Nablaf(NewPoint, Q) + Beta * ConjugateStart
PointList.append(NewPoint)
print("[info] Iterations: ", len(PointList))
InitPoint = NewPoint
ConjugateStart = NewConjugate
return PointList
对应图像表示如下：
很明显，可以发现：使用精确搜索，它仅需要
1
1
1次线搜索：
n
=
2
n=2
n=2次迭代必然能够找到最优解。由于下降方向是共轭方向
d
k
d_k
dk，并且
Q
\mathcal Q
Q不是对角阵，因此上面迭代产生的下降方向之间并不满足垂直关系。
只是看起来优点迷惑人~
当然，如果在计算
α
k
\alpha_k
αk过程中没有乘以
1
2
\begin{aligned}\frac{1}{2}\end{aligned}
21，会得到下面的结果：
循环无法停止，它只会在这几个点上无限循环下去。不要问我是怎么知道的~
基于Wolfe准则的共轭梯度法
在非精确搜索——
Wolfe Condition
\text{Wolfe Condition}
Wolfe Condition一节中介绍了这种方法，它主要通过参数
C
1
\mathcal C_1
C1来描述所选
α
\alpha
α的上界；以及参数
C
2
\mathcal C_2
C2来描述上界以下范围内
α
\alpha
α满足的梯度范围：
其中
ϕ
(
α
)
=
f
(
x
k
+
1
)
=
f
(
x
k
+
α
⋅
d
k
)
\phi(\alpha) = f(x_{k+1}) = f(x_k + \alpha \cdot d_k)
ϕ(α)=f(xk+1)=f(xk+α⋅dk)
{
ϕ
(
α
)
≤
f
(
x
k
)
+
C
1
[
∇
f
(
x
k
)
]
T
d
k
⋅
α
ϕ
′
(
α
)
≥
C
2
⋅
[
∇
f
(
x
k
)
]
T
d
k
C
1
∈
(
0
,
1
)
C
2
∈
(
C
1
,
1
)
\begin{cases} \phi(\alpha) \leq f(x_k) + \mathcal C_1 [\nabla f(x_k)]^T d_k \cdot \alpha \\ \phi'(\alpha) \geq \mathcal C_2 \cdot [\nabla f(x_k)]^T d_k \\ \mathcal C_1 \in (0,1) \\ \mathcal C_2 \in (\mathcal C_1,1) \end{cases}
⎩
⎨
⎧ϕ(α)≤f(xk)+C1[∇f(xk)]Tdk⋅αϕ′(α)≥C2⋅[∇f(xk)]TdkC1∈(0,1)C2∈(C1,1)
线搜索过程中，每次迭代只选择一个满足上述条件的优质结果参与迭代，从而得到近似最优解。关于
Wolfe Condition
\text{Wolfe Condition}
Wolfe Condition的收敛性证明详见传送门。对应代码描述如下：
def WolfeConditionOperation(C1,C2,PointInput,Conjugate,Q,UpperLimits):
def phi(alpha,PointInput,Conjugate,Q):
return np.dot(np.dot(PointInput + alpha * Conjugate,Q),PointInput + alpha * Conjugate)
def Derphi(alpha,PointInput,Conjugate,Q):
# phi()函数关于alpha的导函数
return np.dot(np.dot(Conjugate,Q),PointInput) + np.dot(np.dot(PointInput,Q),Conjugate) \
+ 2 * alpha * np.dot(np.dot(Conjugate,Q),Conjugate)
assert 0 < C1 < 1 and C1 < C2 < 1
while True:
alpha = random.uniform(0.0,UpperLimits)
if phi(alpha,PointInput,Conjugate,Q) <= f(PointInput,Q) + C1 * alpha * np.dot(Nablaf(PointInput,Q),Conjugate):
if Derphi(alpha,PointInput,Conjugate,Q) >= C2 * np.dot(Nablaf(PointInput,Q),Conjugate):
if not alpha:
continue
else:
UpperLimits /= 1.2
break
return alpha
基于
Wolfe Condition
\text{Wolfe Condition}
Wolfe Condition的共轭梯度法收敛效果描述如下：
需要说明的是，每次迭代产生的
α
\alpha
α不是固定的，对应图像也不是固定的。甚至有些时候选择出的
α
\alpha
α结果不满足迭代条件甚至卡死。多试几次~
附：共轭梯度法完整代码
import math
import random
import numpy as np
import matplotlib.pyplot as plt
def CalxLimits(C):
return math.sqrt(2 * C),-1 * math.sqrt(2 * C)
# def f(x,y):
# 没有用到;只是描述一下这个非标准型函数。
# return 2 * (y ** 2) + 2 * x * y + (x ** 2)
def GetAnalytical(x,C):
return 0.25 * (-2 * x + math.sqrt(8 * C - (4 * (x ** 2)))),0.25 * (-2 * x - math.sqrt(8 * C - (4 * (x ** 2))))
def DrawContourPicture(CList):
for C in CList:
UpperPointList,LowerPointList = list(),list()
Upper,Lower = CalxLimits(C)
xList = list(np.linspace(Lower,Upper,10000))
for x in xList:
Upperx,Lowerx = GetAnalytical(x,C)
UpperPointList.append(Upperx)
LowerPointList.append(Lowerx)
plt.plot(xList, UpperPointList,'--',c="tab:blue")
plt.plot(xList, LowerPointList,'--',c="tab:blue")
# plt.show()
def ConjugateGradient(mode="WolfeCondition"):
# 使用精确搜索(步长是最优解)的迭代效果。
def f(PointInput,Q):
# 二次型
return np.dot(np.dot(PointInput,Q),PointInput)
def Nablaf(PointInput,Q):
# 二次型求导
return 2 * np.dot(Q,PointInput)
def GetNorm(ArrayInput):
return np.linalg.norm(ArrayInput)
def WolfeConditionOperation(C1,C2,PointInput,Conjugate,Q,UpperLimits):
def phi(alpha,PointInput,Conjugate,Q):
return np.dot(np.dot(PointInput + alpha * Conjugate,Q),PointInput + alpha * Conjugate)
def Derphi(alpha,PointInput,Conjugate,Q):
# phi()函数关于alpha的导函数
return np.dot(np.dot(Conjugate,Q),PointInput) + np.dot(np.dot(PointInput,Q),Conjugate) \
+ 2 * alpha * np.dot(np.dot(Conjugate,Q),Conjugate)
assert 0 < C1 < 1 and C1 < C2 < 1
while True:
alpha = random.uniform(0.0,UpperLimits)
if phi(alpha,PointInput,Conjugate,Q) <= f(PointInput,Q) + C1 * alpha * np.dot(Nablaf(PointInput,Q),Conjugate):
if Derphi(alpha,PointInput,Conjugate,Q) >= C2 * np.dot(Nablaf(PointInput,Q),Conjugate):
if not alpha:
continue
else:
UpperLimits /= 1.2
break
return alpha
assert mode in ["Exact","WolfeCondition"]
Epsilon = 0.05
InitPoint = np.array([3.,2.])
Q = np.array([[1.,1.],[1.,2.]])
PointList = list()
PointList.append(InitPoint)
UpperLimits = 2.0
C1 = 0.5
C2 = 0.8
ConjugateStart = -1 * Nablaf(InitPoint,Q)
while True:
if GetNorm(Nablaf(InitPoint,Q)) <= Epsilon:
break
else:
if mode == "Exact":
alpha = -0.5 * (np.dot(Nablaf(InitPoint,Q),ConjugateStart) / np.dot(np.dot(ConjugateStart,Q),ConjugateStart))
else:
alpha = WolfeConditionOperation(C1,C2,InitPoint,ConjugateStart,Q,UpperLimits)
NewPoint = InitPoint + alpha * ConjugateStart
Beta = np.dot(np.dot(Nablaf(NewPoint,Q),Q),ConjugateStart) / np.dot(np.dot(ConjugateStart,Q),ConjugateStart)
NewConjugate = -1 * Nablaf(NewPoint,Q) + Beta * ConjugateStart
PointList.append(NewPoint)
InitPoint = NewPoint
ConjugateStart = NewConjugate
return PointList
def DrawPicture(PointList):
CList = [0.01,0.1,0.5,2,4.5]
DrawContourPicture(CList)
plotList = list()
for (x,y) in PointList:
plotList.append((x,y))
plt.scatter(x, y, s=40, facecolor="none", edgecolors="tab:red", marker='o')
if len(plotList) < 2:
continue
else:
plt.plot([plotList[0][0], plotList[1][0]], [plotList[0][1], plotList[1][1]], c="tab:red")
plotList.pop(0)
plt.show()
if __name__ == '__main__':
PointList = ConjugateGradient()
# PointList = ConjugateGradient(mode="Exact")
DrawPicture(PointList)
