常用最优化算法-Python实现

computeGradient 梯度计算方法见文章
理查德外推法计算偏导数近似值-python实现

梯度法

梯度下降法的优化思想是用当前位置负梯度方向作为搜索方向，因为该方向为当前位置的最快下降方向，所以也被称为是”最速下降法“。最速下降法越接近目标值，步长越小，前进越慢

流程

实现

from matplotlib import pyplot
from computeGradient import computeGrad

def fun(x):
    return 100 * (x[1] - x[0] ** 2) ** 2 + (1 - x[0]) ** 2

def gradient(x):
    grad = computeGrad(fun, x, 1e-3)
    return grad

def armijo(fun,xk,dk,gk,beta,sigma,maxm=20):
    """
    非精确线搜索Armijo准则
    :param beta:(0,1)
    :param sigma:(0,0.5)
    :param maxm:最大m值(m 非负整数0,1,2,3,...)
    :return:mk
    """
    m = 0
    mk = 0
    while m < maxm:
        if fun(xk + beta ** m * dk) <= fun(xk) + sigma * beta ** m * np.dot(gk, dk):
            mk = m
            break
        m += 1

    return mk

def Gradient(fun, gfun, x0, epsilon):
    """
    Gradient 梯度法
    :param fun:目标函数
    :param gfun:梯度
    :param x0:初始点
    :param epsilon:容许误差
    :return:k 迭代次数，x 近似最优点，val 近似最优值
    """
    maxIteration = 5000
    beta = 0.5
    sigma = 0.4
    k = 0

    graphArr = []
    graphArr.append(x0)

    while k < maxIteration:
        gk = gfun(x0)
        dk = -1.0 * gk

        if np.linalg.norm(gk) < epsilon:
            break

        mk = armijo(fun,x0,dk,gk,beta,sigma,maxm=20)
        x0 = x0 + beta ** mk * dk
        graphArr.append(x0)
        k += 1

    x = x0
    val = fun(x)
    return k, x, val, graphArr

k, x, val, graphArr = Gradient(fun, gradient, np.array([-1.2, 1.0]), 1e-5)
print "minPoint:", x, "minValue:", int(val), "iteration:", k

arg_x = np.array(graphArr)[:, 0]
arg_y = np.array(graphArr)[:, 1]
pyplot.plot(arg_x, arg_y)
pyplot.show()

结果

minPoint: [ 0.99998906  0.99997808] minValue: 0 iteration: 1435

阻尼牛顿法

基础牛顿法初始点需要靠近极小点，否则算法可能不收敛。为了克服初始点选取困难的问题，阻尼牛顿法引入了线搜索方法以得到大范围收敛。

流程

实现

import numpy as np
from matplotlib import pyplot
from computeGradient import computeGrad
from computeHessian import computeHes

def fun(x):
    return 100 * (x[1] - x[0] ** 2) ** 2 + (1 - x[0]) ** 2

def gradient(x):
    grad = computeGrad(fun, x, 1e-3)
    return grad

def hessianMatrix(x):
    Hes = computeHes(fun, x, 1e-3)
    return Hes

def armijo(fun,xk,dk,gk,beta,sigma,maxm=20):
    """
    非精确线搜索Armijo准则
    :param beta:(0,1)
    :param sigma:(0,0.5)
    :param maxm: 最大m值(m 非负整数0,1,2,3,...)
    :return:mk
    """
    m = 0
    mk = 0
    while m < maxm:
        if fun(xk + beta ** m * dk) <= fun(xk) + sigma * beta ** m * np.dot(gk, dk):
            mk = m
            break
        m += 1

    return mk

def dampingNewton(fun, gfun, Hess, x0, epsilon):
    maxIteration = 5000
    beta = 0.5
    sigma = 0.4
    k = 0

    graphArr = []
    graphArr.append(x0)

    while k < maxIteration:
        gk = gfun(x0)
        Gk = Hess(x0)
        dk = -1.0 * np.linalg.solve(Gk, gk)
        if np.linalg.norm(gk) < epsilon:
            break
        mk = armijo(fun,x0,dk,gk,beta,sigma,maxm=20)
        x0 = x0 + beta ** mk * dk
        k = k + 1
        graphArr.append(x0)

    x = x0
    val = fun(x)

    return k, x, val, graphArr


k, x, val, graphArr = dampingNewton(fun, gradient, hessianMatrix, np.array([-1.2, 1.0]), 1e-5)
print "minPoint:", x, "minValue:", int(val), "iteration:", k

arg_x = np.array(graphArr)[:, 0]
arg_y = np.array(graphArr)[:, 1]

pyplot.plot(arg_x, arg_y)
pyplot.show()

结果

minPoint: [ 0.9994168   0.99883101] minValue: 0 iteration: 5000

拟牛顿法

解决无约束优化问题，在计算过程中不需要计算目标函数的Hesse阵，却在某种意义下具有使用Hesse阵时的功效。其基本思想是，用迭代点处的一阶导数和拟牛顿条件来构造目标函数的曲率近似，然后把极小值作为新的迭代点，并重复这一过程，直至求得满足精度的近似极小值点。

流程

实现: 对称秩1

import numpy as np
from matplotlib import pyplot

from computeGradient import computeGrad

def armijo(fun,xk,dk,gk,beta,sigma,maxm=20):
    """
    非精确线搜索Armijo准则
    :param beta:(0,1)
    :param sigma:(0,0.5)
    :param maxm:最大m值(m 非负整数0,1,2,3,...)
    :return:mk
    """
    m = 0
    mk = 0
    while m < maxm:
        if fun(xk + beta ** m * dk) <= fun(xk) + sigma * beta ** m * np.dot(gk, dk):
            mk = m
            break
        m += 1

    return mk

def sr1Newton(fun, gfun, x0, epsilon, maxIteration):
    """
    基于Armijo的对称秩1算法
    :param fun:目标函数
    :param gfun:梯度
    :param x0:初始点
    :param epsilon:容许误差
    :param maxIteration:最大迭代次数
    :return:k 迭代次数，x 近似最优点，val 近似最优值
    """
    beta = 0.55
    sigma = 0.4
    k = 0
    n = len(x0)
    Hk = np.eye(n)

    graphArr = []
    graphArr.append(x0)

    while k < maxIteration:
        gk = gfun(x0)
        if np.linalg.norm(gk) < epsilon:
            break
        dk = -1.0 * np.dot(Hk, gk)

        mk = armijo(fun, x0, dk, gk, beta, sigma, maxm=20)
        x = x0 + beta ** mk * dk
        sk = x - x0
        yk = gfun(x0) - gk

        sHy = sk - np.dot(Hk, yk)
        Hk = Hk + 1.0 * sHy * sHy.reshape((n, 1)) / 1.0 * sHy.reshape((n, 1)) * yk
        k = k + 1
        x0 = x
        graphArr.append(x0)

    x = x0
    val = fun(x0)

    return k, x, val, graphArr

k, x, val, graphArr = sr1Newton(fun, gradient, np.array([-1.0, 1.0]), 1e-5, 10000)
print "minPoint:", x, "minValue:", int(val), "iteration:", k

arg_x = np.array(graphArr)[:, 0]
arg_y = np.array(graphArr)[:, 1]
pyplot.plot(arg_x, arg_y)
pyplot.show()

结果

minPoint: [ 0.99998982  0.99997959] minValue: 0 iteration: 8121

实现: BFGS

# coding=utf-8
import numpy as np
from matplotlib import pyplot
from computeGradient import computeGrad


def fun(x):
    return 100 * (x[1] - x[0] ** 2) ** 2 + (1 - x[0]) ** 2

def gradient(x):
    # grad = np.array([-400 * x[0] * (x[1] - x[0] ** 2) - 2 * (1 - x[0]), 200 * (x[1] - x[0] ** 2)])
    grad = computeGrad(fun, x, 1e-5)
    return grad

def armijo(fun,xk,dk,gk,beta,sigma,maxm=20):
    """
    非精确线搜索Armijo准则
    :param beta:(0,1)
    :param sigma:(0,0.5)
    :param maxm:最大m值(m 非负整数0,1,2,3,...)
    :return:mk
    """
    gamma = 0.7
    m = 0
    mk = 0
    while m < maxm:
        # gk1 = gfun(x0 + beta ** m * dk)
        # if fun(xk + beta ** m * dk) <= fun(xk) + sigma * beta ** m * np.dot(gk, dk) and np.dot(gk1.T,dk) >= gamma * np.dot(gk.T, dk):
        if fun(xk + beta ** m * dk) <= fun(xk) + sigma * beta ** m * np.dot(gk, dk):
            mk = m
            break
        m += 1

    return mk

def bfgs(fun, gfun, x0, epsilon, maxIteration=1000):
    beta = 0.55
    sigma = 0.4
    gamma = 0.7
    k = 0
    n = np.shape(x0)[0]

    Bk = np.eye(n)

    graphArr = []
    graphArr.append(x0)

    while k < maxIteration:
        gk = gfun(x0) # 计算梯度
        if np.linalg.norm(gk) < epsilon: # 检验终止准则
            break
        # 通过解方程组 Bk * dk = -gk ，计算搜索方向 dk = -Bk/gk
        dk = -1.0 * np.linalg.solve(Bk, gk)
        mk = armijo(fun, x0, dk, gk, beta, sigma, maxm=20)

        # BFGS
        x = x0 + beta ** mk * dk
        # print "Iteration" + str(k) + ":" + str(x)
        sk = x - x0
        yk = gfun(x) - gk

        if np.dot(sk, yk) > 0:
            Bs = np.dot(Bk, sk)
            ys = np.dot(yk, sk)
            sBs = np.dot(np.dot(sk, Bk), sk)

            Bk = Bk - 1.0 * Bs.reshape((n, 1)) * Bs / sBs + 1.0 * yk.reshape((n, 1)) * yk / ys

        k += 1
        x0 = x
        graphArr.append(x0)

    return x0, fun(x0), k, graphArr


x0, fun0, k, graphArr = bfgs(fun, gradient, np.array([-1.2, 1.0]), 1e-5, 10000)
print "minPoint:", x0, "minValue:", int(fun0), "iteration:", k

arg_x = np.array(graphArr)[:, 0]
arg_y = np.array(graphArr)[:, 1]

pyplot.plot(arg_x, arg_y)
pyplot.show()

结果

minPoint: [ 0.99999999  0.99999999] minValue: 0 iteration: 32