优化方法 - 牛顿法 以及Python实现

• 一开始没太搞明白为什么是教材上的迭代方式，其中最为困惑的就是为什么就直接假设以$f(x_0)$展开的$f(x)$处的导数，也就是$f^{\prime }(x)$为0，并且以这个条件引出了迭代公式呢？

• 归根结底，这里直接使用以$x_0$表示出来的$f(x)$还是要从泰勒展开进行理解。$x_0$未必就是实际的极值点处，只是用之作为一个初值表示出$x$处，最后利用极值点的必要条件，也即$f^{\prime }(x)=0$得到迭代公式。

• 几行代码实现牛顿法求解方程最小值，首先举例一维最简单的情况

# find the minimum of y
# y = x^2 + x + 1
# y' = 2x + 1
# y'' = 2

x0 = 91.91  # initial point

g = 2*x0 + 1  # f'(x0)
gg = 2  # f''(x0)

while abs(g) > 0.001:
    # update x0
    x0 = -g/gg + x0

    # update one order partial
    g = 2*x0 + 1

    # update two order partial, well, in this case it is constant
    gg = 2

    print(g)

y = x0**2 + x0 + 1
print('\nx_final is {0}, and min_y is {1}'.format(x0, y))

从实验现象来看，这个简单的例子仅仅迭代了一次就得到了和解析解相同的结果！为什么会这样？

答：因为我们找的迭代步长，直接可以达到一阶导为0的点，对于二次凸函数，这已经是极小值点！

• 一维情况下，没有看到所谓的海塞矩阵，接下来举例二维情况。程序中如何构造H矩阵？
• 答：事先根据偏导公式找到各变量的系数，计算时直接带入。

#  2-Dim Newton
# z = x^2 + y^2 + x + y + xy + 1
# partial x = 2*x + y + 1
# partial y = 2*y + x + 1

x0 = [1, 1]  # nit point

p_x = 2*x0[0] + x0[1] + 1
p_y = 2*x0[1] + x0[0] + 1

p_xx = 2
p_yy = 2
p_xy = 1
p_yx = 1

g0 = [p_x, p_y]
g0_norm = np.linalg.norm(g0)

H0 = np.array([[p_xx, p_xy],
               [p_yx, p_yy]
               ])

print('H0:\n', H0)

while g0_norm > 0.01:
    # print(g0_norm)

    x0 = x0 - 0.1 * np.dot(np.linalg.inv(H0), g0)  # update x0
    # print('x0:', x0)

    # update g0
    p_x = 2 * x0[0] + x0[1] + 1
    p_y = 2 * x0[1] + x0[0] + 1
    g0 = [p_x, p_y]
    g0_norm = np.linalg.norm(g0)  # don't forget to update end_condition!
    print('g0: ', g0)

    # this case, H is constant!

    # check y
    y = x0[0]**2 + x0[1]**2 + x0[0] + x0[1] + x0[0]*x0[1] + 1
    print('y: ', y)

    time.sleep(0.1)  # fastidium delay...

• 参考
  1.牛顿法与拟牛顿法学习笔记（一）牛顿法

  2.理解牛顿法