58、无约束优化框架

import math

def AB(A, B):
    return [[sum([A[i][k]*B[k][j] for k in range(len(A[0]))]) for j in range(len(B[0]))] for i in range(len(A))]

def Ax(A, x):
    return [sum([A[i][j]*x[j] for j in range(len(A[0]))]) for i in range(len(A))]

def xy(x, y):
    return sum([x[i]*y[i] for i in range(len(x))])

def mod(v):
    return math.sqrt(sum(v[i]*v[i] for i in range(len(v))))

def x_next(x, di, ai):
    return [x[i]+di[i]*ai for i in range(len(x))]

# 函数为f(x) = 2(x-2)^2+(y-10)^2
def f(x):
    return 2*(x[0]-2)*(x[0]-2) + (x[1]-10)*(x[1]-10)

# 确定下降方向
# 梯度方向
def d_gradient(x, **kwargs):
    return [-2*x[0]+4, -x[1]+10]

def d_newton(x, **kwargs):
    return [2-x[0], 10-x[1]]

def d_conjugate(x, **kwargs):
    H = [[4, 0], [0, 2]]
    g = [4 * x[0] - 8, 2 * x[1] - 20]

    if kwargs["loop"] == 0:
        return [-g[i] for i in range(len(g))]

    di = kwargs["di"]
    b = xy(Ax(H, di), g)/xy(Ax(H, di), di)
    return [-g[i]+b*di[i] for i in range(len(g))]



# 确定步长
# 三分法步长
def a_three(x, f, di, **kwargs):
    left = 0
    right = 1

    # 放缩法确定区间
    while f(x_next(x, di, right)) < f(x):
        right = right * (2 - kwargs['r'])

    # 三分法确定步长
    mid1 = left * 2 / 3 + right / 3
    mid2 = left / 3 + right * 2 / 3
    while abs(left - right) * mod(di) > kwargs['step_min']:
        if f(x_next(x, di, mid1)) < f(x_next(x, di, mid2)):
            right = mid2
        else:
            left = mid1
        mid1 = left * 2 / 3 + right / 3
        mid2 = left / 3 + right * 2 / 3

    return left

# 牛顿法步长
def a_newton(x, f, di, **kwargs):
    f1 = 4*(x[0]-2)*di[0] + 2*(x[1]-10)*di[1]
    f2 = 4*di[0]*di[0] + 2*di[1]*di[1]
    return -f1/f2 if f2 > 0 else 0

def a_conjugate(x, f, di, **kwargs):
    H = [[4,0],[0,2]]
    g = [4*x[0]-8, 2*x[1]-20]
    return -xy(g, di)/xy(Ax(H, di), di) if xy(Ax(H, di), di) > 0 else 0


def opt(x0, f, d, a, **kwargs):

    x = x0

    print("\nf(x):")
    print("loop  value")
    loop = 0
    di = [0 for i in range(len(x))]

    while True:
        di = d(x, loop = loop, di = di, **kwargs)
        ai = a(x, f, di, **kwargs)
        if abs(ai) * mod(di) < kwargs['step_min']:
            break
        x = x_next(x, di, ai)
        loop += 1
        print(loop, ":  ", f(x))

    print("\nx:\n", x)
    print("\ndi:\n", di)
    print("\nai:\n", ai)
    print("\nmod(di):\n", mod(di))

    return x

#x为优化向量，f为价值函数，这里默认是最小化，d为优化方向，a为优化步长
if __name__ == '__main__':

    x0 = [0, 0]

    print("\nf(x0):\n", f(x0))
    # x = opt(x0, f, d_gradient, a_three, r=0.8, step_min=0.00001)
    # x = opt(x0, f, d_gradient, a_newton, r=0.8, step_min=0.00001)
    # x = opt(x0, f, d_gradient, a_conjugate, r=0.8, step_min=0.00001)
    #
    # x = opt(x0, f, d_newton, a_three, r=0.8, step_min=0.00001)
    # x = opt(x0, f, d_newton, a_newton, r=0.8, step_min=0.00001)
    # x = opt(x0, f, d_newton, a_conjugate, r=0.8, step_min=0.00001)
    #
    x = opt(x0, f, d_conjugate, a_three, r=0.8, step_min=0.00001)
    # x = opt(x0, f, d_conjugate, a_newton, r=0.8, step_min=0.00001)
    # x = opt(x0, f, d_conjugate, a_conjugate, r=0.8, step_min=0.00001)

    print("\nf(x):\n", f(x))