深入淺出最優(yōu)化(1) 最優(yōu)化問題概念與基本知識

1 最優(yōu)化問題

1.1 什么是最優(yōu)化問題

最優(yōu)化問題大體上分為連續(xù)最優(yōu)化問題和離散最優(yōu)化問題兩種祠丝。后者涉及到離散數(shù)學(xué)、組合數(shù)學(xué)等學(xué)科装哆，屬于計(jì)算機(jī)專業(yè)的專業(yè)課程伐厌，而前者的雛形在微積分課程中，甚至在高中褐望、初中勒庄、小學(xué)的數(shù)學(xué)課堂上就有所涉及。

我們還記得求一個(gè)定義域內(nèi)連續(xù)可微的一元實(shí)函數(shù)的最小值點(diǎn)的求解方法：求導(dǎo)瘫里，尋找導(dǎo)數(shù)為0的點(diǎn)实蔽，再求二階導(dǎo)，這些點(diǎn)當(dāng)中二階導(dǎo)大于0的點(diǎn)是我們需要的極值點(diǎn)谨读，再將這些極點(diǎn)對應(yīng)的函數(shù)值進(jìn)行對比局装，函數(shù)值最小的點(diǎn)就是我們需要的最小值點(diǎn)。

上述問題就是一個(gè)最基本的最優(yōu)化問題劳殖，但它已經(jīng)充分描述了最優(yōu)化問題的基本目標(biāo)：求解函數(shù)的最小值點(diǎn)铐尚。當(dāng)然，求解最大值點(diǎn)的方法也是相同的哆姻，因此我們再討論最優(yōu)化問題時(shí)只討論最小化問題宣增。

若函數(shù)定義域?yàn)?img class="math-inline" src="https://math.jianshu.com/math?formula=R" alt="R" mathimg="1">，則該最優(yōu)化問題為無約束最優(yōu)化問題矛缨。否則爹脾，就是約束最優(yōu)化問題。

下面箕昭，針對多元函數(shù)灵妨，我們提出最優(yōu)化問題規(guī)范化的描述。

1.2 名詞與符號

最優(yōu)點(diǎn) $x^*$
最優(yōu)值 $f(x^*)$
最優(yōu)解 $(x^*,f(x^*))^T$
可行域 $D$ 落竹，對應(yīng)于基本問題中的定義域
梯度 $g=\nabla f(x)$ 泌霍，對應(yīng)于基本問題中的一階導(dǎo)數(shù)
黑森矩陣 $G=\nabla^2f(x)$ ，對于n元函數(shù)筋量， $(i,j)\quad i,j\in n$ 位置的值為 $\frac{\partial^2 f(x)}{\partial x_i\partial x_j}$ 烹吵，對應(yīng)于基本問題中的二階導(dǎo)數(shù)
局部最優(yōu)解：在 $x^*$ 附近一個(gè)鄰域內(nèi) $f(x^*)$ 最小碉熄，對應(yīng)基本問題中的極值點(diǎn)
全局最優(yōu)解：在 $D$ 內(nèi) $f(x^*)$ 最小，對應(yīng)基本問題中的最小值點(diǎn)

1.3 最優(yōu)解條件

一階必要條件：f(x)一階連續(xù)可微且 $g^*=\nabla f(x^*)=0$ 是 $x^*$ 為f(x)的局部最優(yōu)解的必要條件肋拔。滿足 $g^*=\nabla f(x^*)=0$ 的點(diǎn)稱為駐點(diǎn)锈津。駐點(diǎn)分為：極小點(diǎn)，極大點(diǎn)凉蜂，鞍點(diǎn)

在這里插入圖片描述

僅從三維空間來看琼梆，梯度為0的幾何意義是該處曲面水平。當(dāng)然這個(gè)曲面還可以拓展到更高維的空間中窿吩，在此不再贅述

二階充分條件：f(x)二階連續(xù)可微且 $\nabla f(x^*)=0, \nabla^2f(x^*)$ 正定茎杂，是 $x^*$ 為f(x)的嚴(yán)格局部最優(yōu)解的充分條件。如果說一階條件刻畫了f(x)在x處切平面的法向纫雁，二階條件就刻畫了曲面的彎曲方向煌往，而 $\nabla^2f(x^*)$ 正定代表了此處曲面向上彎曲

以上兩個(gè)條件雖然拓展到了高維，但實(shí)際上和一元函數(shù)最優(yōu)化問題的條件是通過1.2中給出的對應(yīng)關(guān)系對應(yīng)的轧邪，結(jié)合起來就不難理解

2 用計(jì)算機(jī)求解問題

2.1 迭代搜索

我們最終的目標(biāo)是要實(shí)現(xiàn)讓計(jì)算機(jī)代替我們?nèi)ふ易顑?yōu)化問題的解刽脖。我們不能指望計(jì)算機(jī)具有和我們一樣的智能、對于一個(gè)方程可以使用人類總結(jié)出來的諸如分離變量之類的各種方法去求解忌愚。我們只能使用數(shù)值方法曲管，通過循環(huán)迭代解出一個(gè)逼近于最優(yōu)解的答案。這個(gè)過程稱為迭代搜索硕糊。

我們考慮如圖所示的一個(gè)三維曲面院水， $x_k$ 是當(dāng)前我們所在點(diǎn)。假如這個(gè)點(diǎn)不是我們需要的最優(yōu)解简十，我們就需要移動到下一個(gè)點(diǎn)檬某。移動需要確定方向和距離，因此我們有迭代搜索的基本元素如下：

步長：迭代搜索第k步的長度螟蝙，用 $\alpha_k$ 表示
下降方向：迭代搜索第k步的方向橙喘，用 $d_k$ 表示，沿該方向函數(shù)值需要下降胶逢，這樣才能保證算法的收斂性。其中可行方向表示沿著該方向搜索饰潜，存在使得下一個(gè)點(diǎn)不超出可行域的步長

在這里插入圖片描述

終止準(zhǔn)則：考慮一階必要條件初坠，即 $\nabla f(x)=0$ ，我們將這個(gè)條件作為迭代終止的準(zhǔn)則彭雾。如果我們當(dāng)前所在點(diǎn)滿足準(zhǔn)則碟刺，就不需要再尋找下一個(gè)點(diǎn)了。當(dāng)然薯酝，對于搜索算法來說半沽，我們不能保證我們一定可以通過有限步的搜索找到精確地滿足 $\nabla f(x)=0$ 的點(diǎn)爽柒，而我們在實(shí)際工程運(yùn)用中對于最小值點(diǎn)的要求也并非是絕對精確的。因此者填，引入一個(gè)精確度 $\epsilon$ 浩村，只要梯度滿足 $||g(x_k)||\leq \epsilon$ 即可認(rèn)為滿足了終止準(zhǔn)則。

根據(jù)基本元素占哟，我們可以得出n元函數(shù)迭代搜索的步驟：

給出 $x_0\in R^n,0\leq \epsilon<<1,k=0$
如果在 $x_k$ 點(diǎn)終止準(zhǔn)則被滿足心墅，停止迭代
確定 $x_k$ 點(diǎn)下降方向 $d_k$
確定步長因子 $\alpha_k$
令 $x_{k+1}=x_k+\alpha_kd_k$ ，轉(zhuǎn)步驟2

2.2 質(zhì)量評估

收斂性：對于迭代序列 $x^{(k)}$ 榨乎，嚴(yán)格的收斂定義要求k趨近于正無窮時(shí)等于最優(yōu)點(diǎn)怎燥，弱化的收斂定義要求k趨近于無窮時(shí)一階偏導(dǎo)向量的范數(shù)為0
收斂速度：為p階收斂，則
1. $p=1,0<\beta<1$ 蜜暑，則算法線性收斂
2. $p=1,\beta=0$ 铐姚，則算法超線性收斂
3. $p=2,\beta>0$ ，則算法二階收斂肛捍。二階收斂必定超線性收斂隐绵，且二階收斂的算法收斂速度高于線性收斂的算法
二次終止性：若某個(gè)算法對于任意的正定二次函數(shù)，從任意初始點(diǎn)出發(fā)篇梭，都能經(jīng)過有限步迭代達(dá)到極小值氢橙，則稱算法具有二次終止性

一個(gè)可行的收斂算法必須滿足收斂性，收斂速度越快的算法性質(zhì)越好恬偷，而具有二次終止性的算法可以在求解二次函數(shù)最優(yōu)化問題時(shí)具有廣泛的應(yīng)用價(jià)值悍手。但證明這些性質(zhì)并不作為本系列教程的重點(diǎn)。本教程會給出每一種算法性質(zhì)的結(jié)論與相關(guān)證明教程鏈接袍患，有興趣的讀者可以自行查閱

3 最小二乘問題——無約束最優(yōu)化問題實(shí)例

點(diǎn)列的曲線擬合是我們高中開始就接觸過的問題坦康。為了尋找一個(gè)待定系數(shù)的函數(shù)，可以以最小的誤差去描述點(diǎn)列诡延，我們需要用到最小二乘法滞欠。有關(guān)最小二乘法可以參閱：https://www.zhihu.com/question/37031188

最小二乘法是我們研究無約束最優(yōu)化問題的一個(gè)出色的實(shí)例。它具有廣泛的應(yīng)用價(jià)值肆良，而且目標(biāo)函數(shù)的局部最優(yōu)解就是全局最優(yōu)解筛璧。我們在接下來的無約束最優(yōu)化問題部分將針對一個(gè)最小二乘法例題進(jìn)行研究。

例題：

$min\displaystyle\sum_{i=1}^{11}(y_i-\hat{f}(x,t_i))^2,\hat{f}(x,t)=\frac{x_1(t^2+x_2 t)}{t^2+x_3t+x_4}$

在這里插入圖片描述

給出了11組 $(t_i,y_i)$ 的數(shù)據(jù)惹恃，要求我們通過最小二乘法解出 $x_i,i\in1,2,3,4$ 夭谤，使得函數(shù) $\hat{f}(x,t)$ 可以最好地?cái)M合 $y=f(t)$

代碼實(shí)現(xiàn)

本博客所有代碼在https://github.com/HarmoniaLeo/optimization-in-a-nutshell開源，如果幫助到你巫糙，請點(diǎn)個(gè)star朗儒，謝謝這對我真的很重要！

本次代碼實(shí)現(xiàn)內(nèi)容是梯度、梯度二范數(shù)醉锄、黑森矩陣等最優(yōu)化問題基本件的python數(shù)值解法乏悄，以及簡化的python線性代數(shù)庫

梯度、梯度二范數(shù)恳不、黑森矩陣求法參考自博客檩小，有修改：https://www.cnblogs.com/kisetsu/p/9145393.html

制作一個(gè)文件Fuction.py，在里面寫入：

import numpy as np

class Function:
    def __init__(self, _f):
        self.fun = _f

    def value(self, val):
        return self.fun(val)

    def part(self, var_index, val):
        a = self.fun(val)
        b = a + 1
        i = 0
        e = 2 ** 10 - 1
        e1 = 2 ** 10
        while 10 ** (-6) < e < e1 or i > -6:
            e1 = e
            a = b
            val_ = list(val)
            val_[var_index] += 10 ** i
            m = self.fun(val_)
            n = self.fun(val)
            b = (m - n) / 10 ** i
            i -= 2
            e = abs(b - a)
        return a

    def part_2(self, x_index, y_index, val):
        return self.__diff_(x_index).__diff_(y_index).value(val)

    def diff(self, val):
        a = self.fun(val)
        b = a + 1
        i = 0
        e = 2 ** 10 - 1
        e1 = 2 ** 10
        while 10 ** (-6) < e < e1 or i > -6:
            e1 = e
            a = b
            val_ = val + 10 ** i
            m = self.fun(val_)
            n = self.fun(val)
            b = (m - n) / 10 ** i
            i -= 2
            e = abs(b - a)
        return a

    def grad(self, val):
        g = np.array(val).astype('float')
        for i in range(0, g.size):
            g[i] = self.part(i, val)
        return np.array(g)

    def __diff_(self, index):
        def diff_f(vals):
            vals_ = list(vals)
            vals_[index] = vals_[index] + 10 ** (-6)
            m = self.fun(vals_)
            n = self.fun(vals)
            return (m - n) / 10 ** (-6)
        return Function(diff_f)

    def hesse(self, val):
        v = np.mat(val)
        G = np.mat(np.dot(v.T, v)).astype('float')
        for i in range(0, v.size):
            for j in range(0, v.size):
                p = self.part_2(i, j, val)
                G[i, j] = p
        G=np.array(G)
        return G
    
    def norm(self, val):
        s = 0
        for x in self.grad(val):
            s += x ** 2
        return np.sqrt(s)

在之后的每個(gè)方法實(shí)現(xiàn)中妆够，我們都會用到Fuction.py

from Function import Function

tar=Function(func)  #func是一個(gè)以列表為參數(shù)的多元函數(shù)
val=tar.val(para)   #傳入?yún)?shù)獲取值
diff=tar.diff(para) #一元函數(shù)求導(dǎo)
part=tar.part(index,para)   #對第index個(gè)參數(shù)求導(dǎo)
grad=tar.grad(para) #求梯度數(shù)組
hesse=tar.hesse(para)   #求hesse陣
norm=tar.norm(para) #求梯度范數(shù)

另外制作一個(gè)文件lagb.py识啦，在里面寫入：

import numpy as np

def turn(a):    #轉(zhuǎn)置
    if len(a.shape)==1:
        a=a.reshape((1,a.shape[0]))
        return a
    else:
        return np.array(np.mat(a).T)

def dot(a,b):   #點(diǎn)乘
    if len(a.shape)==1:
        a=a.reshape((a.shape[0],1))
    if len(b.shape)==1:
        b=b.reshape((b.shape[0],1))
    res=np.array(np.mat(a)*np.mat(b))
    if res.shape==(1,1):
        return res[0][0]
    else:
        if res.shape[1]==1:
            return res.squeeze()
        else:
            return res

def muldot(*args):  #連乘
    res=args[0]
    for i in range(1,len(args)):
        res=dot(res,args[i])
    return res

def rev(a): #求逆
    return np.array(np.mat(a).I)

在之后的每個(gè)方法實(shí)現(xiàn)中，我們也都會用到lagb.py

from lagb import *

#輸入?yún)?shù)可以為向量或矩陣神妹，即numpy創(chuàng)建的數(shù)組
#用numpy創(chuàng)建的一維數(shù)組會默認(rèn)為列向量
a=turn(a)   #求轉(zhuǎn)置
c=dot(a,b)  #點(diǎn)乘颓哮，視情況會返回向量、矩陣鸵荠、數(shù)值
d=dot([a,b,c])  #連續(xù)的點(diǎn)乘
b=rev(a)    #求逆矩陣