信赖域求解Rosenbrock函数（LM求解子问题）

二维求解

import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3Dfrom pandas import DataFramedef draw(x1, x2):    '''    采用matplolib绘制曲面图    :param x: x轴坐标数组    :param y: y轴坐标数组    :param z: z轴坐标数组    :return:    '''    X = x1    Y = x2    Z = (1 - x1) ** 2 + 100 * (x2 - x1 ** 2) ** 2    fig = plt.figure()    ax = fig.add_subplot(111, projection='3d')    ax.plot_trisurf(X, Y, Z)    plt.show()def draw2(x1,x2):    x = x1    y = x2    #绘制Xk的离散点图形    plt.xlim(xmax = 1,xmin = 0)    plt.ylim(ymax = 1,ymin = 0)    plt.plot(x,y,'ro')    plt.xlabel('x')    plt.ylabel('y')    plt.show()def cal_rosenbrock(x):    # 计算rosenbrock函数的值    return sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0)def cal_rosenbrock_prax(x1, x2):    """    对x1求偏导    """    return -2 + 2 * x1 - 400 * (x2 - x1 ** 2) * x1def cal_rosenbrock_pray(x1, x2):    """    对x2求偏导    """    return 200 * (x2 - x1 ** 2)def rosen_hess(x):    x = np.asarray(x)    H = np.diag(-400*x[:-1],1) - np.diag(400*x[:-1],-1)    diagonal = np.zeros_like(x)    diagonal[0] = 1200*x[0]**2-400*x[1]+2    diagonal[-1] = 200    diagonal[1:-1] = 202 + 1200*x[1:-1]**2 - 400*x[2:]    H = H + np.diag(diagonal)    return Hdef check(T):    c = 0    for i in T[0:]:        if i <= 0:            c = 1            break    return cdef for_rosenbrock_func(max_iter_count=100000,step_size = 0.001):    pre_x = np.zeros((2,), dtype=np.float32)    X_all = np.zeros((2,), dtype=np.float32)    loss = 10    cnt = 0    I = np.eye(2)    miu = 1    loss_all = []    while loss > 0.001 and cnt < max_iter_count:        der = np.zeros((2,), dtype=np.float32)        der[0] = cal_rosenbrock_prax(pre_x[0], pre_x[1])        der[1] = cal_rosenbrock_pray(pre_x[0], pre_x[1])        Hess = np.zeros((2,), dtype=np.float32)        G_k = rosen_hess(pre_x) + miu * I        T_vector = np.linalg.eigvals(G_k)  #计算特征值        # 如果不正定，阻尼系数扩大四倍        while check(T_vector) :            miu = 4 * miu            G_k = rosen_hess(pre_x) + miu * I            T_vector = np.linalg.eigvals(G_k)        S_k = np.linalg.solve(G_k,-1 * der)        error1 = cal_rosenbrock(pre_x) - cal_rosenbrock(pre_x + S_k)        error2 = -1.0 * np.dot(der,S_k.transpose()) + 0.5 * np.dot(np.dot(S_k.transpose(),rosen_hess(pre_x)),S_k)        r = error1 / error2        #比较r值        #如果r大于0        if r > 0:            pre_x = pre_x + S_k            # 如果r值小于0.25，步子太大，缩小范围，令 μk+1=4μk            if r < 0.25:                miu = 4 * miu            #  如果r大于0.75 ，到信赖域边缘，步子有点小，令 μk+1=μk/2             elif r > 0.75:                miu = miu / 2.0            # 如果r大于等于0.25且小于等于0.75 ，介于信赖域与非信赖域之间，令 μk+1=μk            elif r >= 0.25 and r <= 0.75:                miu = miu        # r小于0，即沿反方向上升，应停步，按照r小于0.25缩小范围        elif r <= 0:                pre_x = pre_x                miu = 4 * miu        loss = cal_rosenbrock(pre_x)  # 最小值为0        X_all = np.vstack((X_all,pre_x))        print("count: ", cnt, "the loss:", loss,"μ:",miu)        loss_all.append(loss)        cnt += 1    x = X_all[:,0]    y = X_all[:,1]    z = loss_all    draw(x,y)    draw2(x,y)    return pre_xif __name__ == '__main__':    w = for_rosenbrock_func()      print(w)

11维函数求解

import numpy as np import matplotlib as plt def cal_rosenbrock(x):    # 计算rosenbrock函数的值    return sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0)def rosen_der(x):    xm = x[1:-1]    xm_m1 = x[:-2]    xm_p1 = x[2:]    der = np.zeros_like(x)    der[1:-1] = 200*(xm-xm_m1**2) - 400*(xm_p1 - xm**2)*xm - 2*(1-xm)    der[0] = -400*x[0]*(x[1]-x[0]**2) - 2*(1-x[0])    der[-1] = 200*(x[-1]-x[-2]**2)    return derdef rosen_hess(x):    x = np.asarray(x)    H = np.diag(-400*x[:-1],1) - np.diag(400*x[:-1],-1)    diagonal = np.zeros_like(x)    diagonal[0] = 1200*x[0]**2-400*x[1]+2    diagonal[-1] = 200    diagonal[1:-1] = 202 + 1200*x[1:-1]**2 - 400*x[2:]    H = H + np.diag(diagonal)    return Hdef check(T):    c = 0    for i in T[0:]:        if i <= 0:            c = 1            break    return cdef for_rosenbrock_func(max_iter_count = 1000):    pre_x = np.zeros((11,),dtype=np.float32)    loss = 10    cnt = 0    I = np.eye(11)    miu = 1    while loss > 0.001 and cnt < max_iter_count:        der = np.zeros((11,), dtype=np.float32)        Hess = np.zeros((11,), dtype=np.float32)        # 计算梯度,hessian矩阵        der = rosen_der(pre_x)        G_k = rosen_hess(pre_x) + miu * I        T_vector = np.linalg.eigvals(G_k)  #计算特征值        # 如果不正定，阻尼系数扩大四倍        while check(T_vector) :            miu = 4 * miu            G_k = rosen_hess(pre_x) + miu * I            T_vector = np.linalg.eigvals(G_k)        S_k = np.linalg.solve(G_k,-1 * der)        error1 = cal_rosenbrock(pre_x) - cal_rosenbrock(pre_x + S_k)        error2 = -1.0 * np.dot(der,S_k.transpose()) + 0.5 * np.dot(np.dot(S_k.transpose(),rosen_hess(pre_x)),S_k)        r = error1 / error2        #比较r值        #如果r大于0        if r > 0:            pre_x = pre_x + S_k            # 如果r值小于0.25，步子太大，缩小范围，令 μk+1=4μk            if r < 0.25:                miu = 6 * miu            #  如果r大于0.75 ，到信赖域边缘，步子有点小，令 μk+1=μk/2             elif r > 0.75:                miu = miu / 2.0            # 如果r大于等于0.25且小于等于0.75 ，介于信赖域与非信赖域之间，令 μk+1=μk            elif r >= 0.25 and r <= 0.75:                miu = miu        # r小于0，即沿反方向上升，应停步，按照r小于0.25缩小范围        elif r <= 0:                pre_x = pre_x                miu = 6 * miu        loss = cal_rosenbrock(pre_x)  # 最小值为0        print("count: ", cnt, "the loss:", loss,"μ:",miu)        cnt += 1    return pre_xif __name__ == '__main__':    w = for_rosenbrock_func()      print(w)