机器学习笔记之无约束优化问题——(阶段性收尾)共轭方向法与Wolfe准则优化方法Python示例

机器学习笔记之无约束优化问题——基于共轭方向法与Wolfe准则优化方法的Python示例

• 引言
• • 小插曲：画图——非标准二次型的等值线
  • 算法在图像中的表示
  • 附：共轭梯度法完整代码

引言

本节使用 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=(11​12​);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​)(11​12​)(x1​x2​​)=(x1​+x2​x1​+2x2​)(x1​x2​​)=x12​+2x1​x2​+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​+2x1​x2​+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=±2C ​C>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)=21​xTQx+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​)]Tdk​C1​∈(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)