-
一开始没太搞明白为什么是教材上的迭代方式,其中最为困惑的就是为什么就直接假设以 f ( x 0 ) f(x_0) f(x0)展开的 f ( x ) f(x) f(x)处的导数,也就是 f ′ ( x ) f'(x) f′(x)为0,并且以这个条件引出了迭代公式呢?
-
归根结底,这里直接使用以 x 0 x_0 x0表示出来的 f ( x ) f(x) f(x)还是要从泰勒展开进行理解。 x 0 x_0 x0未必就是实际的极值点处,只是用之作为一个初值表示出 x x x处,最后利用极值点的必要条件,也即 f ′ ( x ) = 0 f'(x)=0 f′(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...