牛顿（Newton）方法

1、基本Newton方法

设f(x)f(x)具有连续的二阶偏导数，当前迭代点是xk.f(x)x_k. f(x)在xkx_k处的TaylorTaylor展开式为

f(xk+d)=fk+gTkd+12dTGkd+o(||d||2),

f(x_k+d)=f_k+g_k^Td+\frac{1}{2}d^TG_kd+o(||d||^2),
其中 d=x−xk.d = x-x_k.在点 xkx_k的邻域内，用二次函数

qk(d)=Δfk+gTkd+12dTGkd

q_k(d)\overset{\Delta}{=}f_k+g_k^Td+\frac{1}{2}d^TG_kd
近似 f(xk+d)f(x_k+d),求解问题

min qk(d)=fk+gTk+12dTGkd.

min\ q_k(d)=f_k+g_k^T+\frac{1}{2}d^TG_kd.
若 GkG_k可逆，则方程组

Gkd=−gk

G_kd=-g_k
的解 dk=−G−1kgkd_k=-G_k^{-1}g_k为上面问题 min qk(d)min\ q_k(d)的唯一解.
迭代步骤：
步1　给出x0∈Rn,ϵ>0,k:=0;x_0\in\mathbf R^n,\epsilon>0,k:=0;
步2　若终止准则满足，则输出有关信息，停止迭代；
步3　计算dkd_k;
步4　xk+1:=xk+dk,k:=k+1,x_{k+1}:=x_k+d_k,k:=k+1,转步2.
下面是实现的代码

#include <math.h>
#include <stdio.h>#define delta 1E-03
#define epsilon 1E-06double funcValue(double s)
{return pow(s,4)-6*pow(s,3)+1;
}double derivation(double s)
{double diff = funcValue(s+delta)-funcValue(s);return diff/delta;
}double dederivation(double s)
{double diff = derivation(s+delta)-derivation(s);return diff/delta;
}
int main()
{double x = 1;while(abs(derivation(x))>epsilon){printf("%f\n", derivation(x));double d = -derivation(x)/dederivation(x);x = x+d;}printf("%f, %f", x, funcValue(x));return 0;
}

经过各种试验之后，发现基本的Newton法依赖于初始点的选择，当初始点接近极小点时，迭代序列很快收敛于极小点，否则会出现迭代序列收敛到鞍点或极大点的情形，这会使迭代失败。
由收敛性证明也可得知，它的收敛是局部的。

2、阻尼Newton方法

阻尼牛顿法比基本牛顿法多了一步一维搜索步长α\alpha的方法：

xk+1=xk+αkdk,

x_{k+1}=x_k+\alpha_kd_k,
此法能保证对正定的 Gk,{fk}G_k,\{f_k\}单调下降；即使 xkx_k离 x∗x^*稍远，该方法产生的点列 {xk}\{x_k\}仍可能收敛至 x∗x^*.
其中一维搜索步长的方法有黄金比例法、多项式插值法、信赖域法等。
这里介绍一下黄金比例法（0.618比例法）：
首先我们需要确定一个初始区间，使其包含

ϕ(α)=f(xk+αdk),α>0

\phi(\alpha)=f(x_k+\alpha d_k), \alpha>0
的极小点，这是建立在满足单峰函数的基础上的. 即使不是单峰函数，可以先用进退法得到一个满足单峰函数的区间. 这是迭代步骤：
步1　给定α0∈[0,∞),γ0>0,t>1,i:=0.\alpha_0\in[0,\infty),\gamma_0>0,t>1,i:=0.
步2　计算αi+1=αi+γi.\alpha_{i+1}=\alpha_i+\gamma_i.若αi+1≤0\alpha_{i+1}\le0,则令αi+1:=0,\alpha_{i+1}:=0,转步4；若ϕ(αi+1)≥ϕ(αi),\phi(\alpha_{i+1})\ge\phi(\alpha_i),则转步4.
步3　令γi+1=tγi,α:=αi,αi:=αi+1,i:=i+1\gamma_{i+1}=t\gamma_i,\alpha:=\alpha_i,\alpha_i:=\alpha_{i+1},i:=i+1,转步2.
步4　若i=0,i=0,令γi:=−γi,α:=αi+1\gamma_i:=-\gamma_i,\alpha:=\alpha_{i+1},转步2；否则

a=min{α,αi+1},b=max{α,αi+1},

a=min \{\alpha,\alpha_{i+1}\},\\b=max\{\alpha,\alpha_{i+1}\},
输出a,ba,b,迭代停止.
0.618方法的迭代步骤：
步1　给定a0>0,b0>0,i:=0,ϵ>0,τ:=0.618.a_0>0,b_0>0,i:=0,\epsilon>0,\tau:=0.618.
步2　若bi−ai<ϵ,b_i-a_i则α∗:=bi+ai2\alpha^*:=\frac{b_i+a_i}{2},输出α∗\alpha^*,停止线搜索.
步3　计算

αli:=ai+(1−τ)(bi−ai),αri:=ai+τ(bi−ai).

\alpha_i^l:=a_i+(1-\tau)(b_i-a_i),\\\alpha_i^r:=a_i+\tau(b_i-a_i).
步4　若ϕ(αli)<ϕ(αri),\phi(\alpha_i^l)则ai+1:=ai,bi+1:=αri;a_{i+1}:=a_i,b_{i+1}:=\alpha_i^r;否则ai+1:=αli,bi+1:=bi,i:=i+1a_{i+1}:=\alpha_i^l,b_{i+1}:=b_i,i:=i+1,转步2.
后面的编程实现可以参考这位作者写的文章： http://dataunion.org/20714.html，里面有python写的详细实现方式。

参考资料:
[1]高立.数值最优化方法.北京大学出版社.高立.\color{red}{数值最优化方法}.北京大学出版社.