牛顿迭代法
newton_iteration_method
阿贝尔—鲁菲尼定理指出，五次及更高次的代数方程没有一般的代数解法，即这样的方程不能由方程的系数经有限次四则运算和开方运算求根。
即对于低次n≤4多项式方程，根的存在性是可以通过构造的方式得到的，可以得到求根公式。但n≥5时一般没有求根公式。
高斯的代数基本定理：任何复系数一元n次多项式 方程在复数域上至少有一根(n≥1)，由此推出，n次复系数多项式方程在复数域内有且只有n个根（重根按重数计算）。代数基本定理在代数乃至整个数学中起着基础作用。 据说，关于代数学基本定理的证明，现有200多种证法。
一元函数
牛顿迭代法（Newton's method）又称为牛顿-拉夫逊（拉弗森）方法（Newton-Raphson method）
多数方程不存在求根公式，牛顿迭代法使用函数的泰勒级数的前几项来寻找方程的根
设\(r\)是\(f(x)=0\)的根，选取\(x_0\)作为\(r\)的初始近似值；
过点\((x_0，f(x_0))\)做曲线\(y=f(x)\)的切线\(L_1\)，切线\(L_1\)的方程为\(y=f(x_0)+f'(x_0)(x-x_0)\)，切线L_1与x轴有交点\(x_1\)，则\(x_1=x_0-f(x_0)/f'(x_0)\)；
过点\((x1，f(x1))\)做曲线\(y=f(x)\)的切线\(L_2\)，切线\(L_2\)的方程为\(y=f(x_1)+f'(x_1)(x-x_1)\)，切线\(L_2\)与\(x\)轴有交点\(x_2\)，则\(x_2=x_1-f(x_1)/f'(x_1)\)；
……
\(n+1\)次迭代公式：\(x_{n+1}=x_n-f(x_n)/f'(x_n)\)
\(n→∞，x_{n+1}→x_n→r\)，故\(f(x_n)/f'(x_n)=0\)，进而\(f(x_n)=0\)
1.以\(f(x)=x^2-x-1\)为例，求其零点。
\(x_{n+1}=x_n-f(x_n)/f'(x_n)=（x_n^2+1)/（2*x_n-1)\)
\(x_0\)选取值2和-1的情况，迭代5次，误差\(10^{-8}\)
| \(i\) | \(x_i\) | \(x_{i+1}\) | \(x_i\) | \(x_{i+1}\) | |
| 0 | 2.0000 | 1.6667 | -1.0000 | -0.6667 | |
| 1 | 1.6190 | 1.6190 | -0.6667 | -0.6190 | |
| 2 | 1.6180 | 1.6180 | -0.6190 | -0.6180 | |
| 3 | 1.6180 | 1.6180 | -0.6180 | -0.6180 | |
| 4 | 1.6180 | 1.6180339887 | -0.6180 | -0.6180339887 | |
| 根 | (1+√5)/2=1.61803398875 | (1-√5)/2=-0.61803398875 |
选取初始值会影响最后根的结果
2.以\(f(x)=e^x-x^2\)为例，求其零点。\(（e=2.7182818284590452353602874713527)\)
\(f'(x)=e^x-2x\)
\(x_{n+1}=x_n-f(x_n)/f'(x_n)=x_n-(e^x-x^2)/(e^x-2x)\)
\(x_0\)选取值-1和0.5的情况，都只能迭代到-0.70346附近，即只能找到这一个根
| \(i\) | \(x_i\) | \(x_{i+1}\) | \(x_i\) | \(x_{i+1}\) | |
| 0 | -1.0000 | -0.7330 | 0.5 | -1.6561 | |
| 1 | -0.7330 | -0.7038 | -1.6561 | -0.9277 | |
| 2 | -0.7038 | -0.7035 | -0.9277 | -0.721 | |
| 3 | -0.7035 | -0.7034674683 | -0.721 | -0.7036 | |
| 4 | -0.7036 | -0.7034674283 | |||
| \(e^x-x^2\) | -0.000000087166031 | -0.000000010998992 |
显然\(f(x)=0\)只有1个根，零点位于（-1，0)
3.以\(f(x)=e^x-x^6\)为例，求其零点。
\(f'(x)=e^x-6x^5\)
\(x_{n+1}=x_n-f(x_n)/f'(x_n)=x_n-(e^x-x^6)/(e^x-6x^5)\)
显然\(f(x)=0\)只有\(3\)个根，分别位于\(（-1，0)（0，2)（2，100)\)；初始值依次选取0，2，20；通过定时器记录1000次运行的总时间，得到平均每一次运行的时间
| \(X\) | \(n\) | \(erro=e^X-X^6\) | Time |
| -0.8656497053 | 6 | -0.000 000 003237437873 | 0.00012441060000000000 |
| 1.2268886960 | 8 | -0.000 000 000006451728 | 0.00016655110000000001 |
| 16.9988873523 | 9 | 0.000 000 070780515671 | 0.00020211259999999999 |
3.1优化算法，对\(e^x\)泰勒展开到8阶
\(e^x=1+x+x^2/2+x^3/6+x^4/24+x^5/120+x^6/720+x^7/7/720+x^8/56/720\)
带入迭代公式中，初始值依次选取0，2，20
| \(X\) | \(n\) | \(erro=e^X-X^6\) | Time |
| -0.8656499126 | 6 | -0.000 000 695075708823 | 0.00012678300000000000 |
| 1.226887208 | 8 | 0.000 019750994960432 | 0.00018222190000000000 |
| 1.2268872075 | 21 | 0.000 019751001404167 | 0.00047578460000000003 |
二者都是调用函数形式
显然，也能取到三个近似根，精度，迭代次数和运行时间上不占优势，这是没有采用指数形式，精确度的不够，但是或许在某些场景会更快？可能吗，目前还不知道matlab的底层算法是怎样的（挖坑）
牛顿迭代法的特点和优缺点
特点：
- 牛顿法收敛速度为二阶，对于正定二次函数一步迭代即达最优解
- 牛顿法是局部收敛的，当初始点选择不当时，往往导致不收敛
- 牛顿法不是下降算法，当二阶海塞矩阵非正定时，不能保证产生方向是下降方向
- 二阶海塞矩阵必须可逆，否则算法求解过程进行困难
- 对函数要求苛刻（二阶连续可微，海塞矩阵可逆），而且运算量大
优点：
- 对于二次连续可导的凸函数，二阶收敛很强，收敛速度快
- 求解多元非线性方程：牛顿迭代法可以同时求解多元非线性方程
- 可以求解高精度解：牛顿迭代法可以求得高精度的解，因此适用于精确计算
缺点：
- 函数要求显性：牛顿迭代法只适用于二次函数及其幂级数展开，对其他形式的函数不适用
- 可能不收敛：牛顿迭代法有时可能不收敛，尤其是在函数的二阶导数不连续或不存在的情况下
- 函数零点数量不确定
- 初始点选择影响：牛顿迭代法的结果与初始点选择有关，如果初始点选择不当，可能导致不收敛；若初始点离零点不够近，牛顿法可能会发散，或者变得十分低效；Carmark平方根倒数算法对牛顿迭代法初值的选择可微经典
快速收敛式
拉马努金圆周率公式
\(\frac{1}{\pi}=\frac{2 \sqrt{2}}{9801} \sum_{k=0}^{\infty} \frac{(4 k) !(1103+26390 k)}{(k!)^{4} 396^{4 \mathrm{k}}}\)
\(e^ {x} =1+x+\frac {x^ {2}}{2!} + \frac {x^ {3}}{3!} + \cdots\)
\(\mathrm{e}=\lim _{x \rightarrow \infty}\left(1+\frac{1}{x}\right)^{x}\)
\(\lim _ {n\rightarrow \infty } (1+\frac {x}{n})^ {n} = e^ {x}\)
\(e^x=(e^{x/1000})^{1000}\)
代码实现
三个.m文件，如下
function y=primitive_function(x)
e = 2.7182818284590452353602874713527;
x = x-0; %x0处展开
% e = 2.718281828459;
% y=x*x-x-1; %函数f(x)的表达式
% y=(1+x+x^2/2+x^3/6+x^4/24+x^5/120+x^6/720+x^7/7/720+x^8/56/720)-x^6; %函数e^x-x^2的x=x0处泰勒展开表达式
y=e^x-x^6;
end
function z = derived_function(x)
e = 2.7182818284590452353602874713527;
x = x-0; %x0处展开
% e = 2.718281828459;
% z = x*2-1; %函数f(x)的导数表达式z
% z = (1+x+x^2/2+x^3/6+x^4/24+x^5/120+x^6/720+x^7/7/720)-6*x^5; %2种方法，导数后泰勒展开，或者泰勒展开后求导;
% 太慢了，牛顿法改进
z = e^x-6*x^5;
end
%主程序，计算e^X-X*X+1的零点，只能一个
tic; %启动一个定时器
e = 2.7182818284590452353602874713527;
l=0; %记录循环程序次数
maxl=1000; %循环程序总次数
while l<maxl
X=0; %迭代初值
% x0=4; %泰勒展开处的值
i=1; %迭代次数计算
I=1000; %最大迭代次数
while i
Xn=X-primitive_function(X)/derived_function(X); %牛顿迭代格式
% Xn=X-(X*X-X-1)/(X*2-1); %或者直接函数表达
% Xn=X-(e^X-X^6)/(e^X-6*X^5);
disp([X Xn]); %显示每次迭代值
if abs(Xn-X)>0.00000001 && i<I && abs(e^X-X*X)>0.0000000001 %收敛判断
X=Xn;
else
break
end
i=i+1;
end
l=l+1;
end
t=toc;%显示自定时器tic开启到当前所经历的时间。若定时器没有运行，则toc命令返回0
fprintf('\n %s %.20f','time=',t/maxl)
if i==I
disp('增加迭代次数或者减小误差或者增加展开级数'); %输出结果
else
fprintf('\n %s %.10f \t %s %d \t %s %.18f %d','X=',X,'i=',i,'erro=',e^X-X^6) %输出结果
end
多元函数的牛顿迭代
函数\(y=g(x)\)的零点\(g(x)=0\)，极小值点\(g'(x)=0\)（黑塞矩阵正定），极大值点\(g'(x)=0\)（黑塞矩阵负定）；只要令\(f(x)=g(x)\)或者\(g'(x)\)即可
迭代公式：\(x_{n+1}=x_n-f(x_n)/f'(x_n)\)
\(n→∞，x_{n+1}→x_n→r\)，故\(f(r)/f'(r)=0\)，分母不为0进而\(f(r)=0\)，得到零点
\(n\)元函数\(g(x_1,x_2,……,x_n)\)的极值点，对进行偏导，\(g\)对\(x_i\)的偏导\(∂g/∂x_i\)记作\(g_i\)，即获得\(n\)个\(g_i(x_1,x_2,……,x_n)\)，联立偏导函数\(g_i(x_1,x_2,……,x_n)=0\)时变元\(x_i\)的值就是极小值点。\(（i=1，2，3，……,n）（g：R^n→R）\)。
\(g_1（x_1,x_2,……,x_n)=0\)
\(g_2（x_1,x_2,……,x_n)=0\)
……
\(g_n（x_1,x_2,……,x_n)=0\)
矩阵形式，令\(X=[x_1,x_2,……,x_n]\)，\(f(X)=[g_1(X)，g_2(X)，……，g_n(X)]\)
所以\(f(X)=0（f：R^n→R^n)\)
牛顿迭代法：
\(X_k=X_{k-1}-f(X_k)/f`(X_{k-1})=X_{k-1}-f(X_k)·[f`(X_{k-1})]^{-1}=X_{k-1}-g`(X_k)·[g``(X_{k-1})]^{-1}\)
其中\([g''(X_{k-1})]^{-1}\)为\(n*n\)矩阵，并且极小值要求是正定的，即每个特征值恒为正（也极大值时要求负定矩阵，反正特征值不能为正负分布)
以求\(g(x，y)=x^4+y^4\)的极小值点为例
Sol
\(Let X=（x，y)，then g'(X)=（g_x，g_y)=（4x^3，4y^3),\)
\(g''(X_{k-1})=▽g'=dg' (X)/dX=d（g_x，g_y)/d（x，y)=(\{g_{xx}，g_{xy}\}，\{g_{yx}，g_{yy}\})=({12x^2，0}，{0，12y^2})\)
\(\begin{bmatrix} g_{xx} & g_{xy} \\ g_{yx} & g_{yy}\end{bmatrix}=\begin{bmatrix} 12x^2 & 0 \\ 0 & 12y^2\end{bmatrix}\)
\(Let X_0=（x_0，y_0)\),其中\(x_0=a，y_0=b\)，则
\(X_1=X_0-g'(X_0)·[g''(X_0)]^{-1}=（a，b)-（4a^3，4b^3)({1/（12a^2)，0}，{0，1/（12b^2)})=2/3（a，b)=2/3·X_0\)
\(k→∞\)时，\(X_k=（2/3)^k·X_0=（0，0)\)
\(g（x，y)\)极小值点为\(g（0，0)=0\)，黑塞矩阵是正定的
偏导矩阵
雅可比矩阵是函数的一阶偏导数以一定方式排列成的矩阵，从一个 n 维的欧式空间转换到 m 维的欧氏空间，可以不是方阵
使用雅克比矩阵和{任意点与p点的距离} 来估算x所对应的f值，pytorch和tensorflow利用雅可比矩阵实现梯度下降法
黑塞矩阵（Hessian Matrix）为函数的二阶偏导数形成的矩阵，为方阵
若\(f\)在\(D\)域内连续可导，则其黑塞矩阵为对称阵
诸如牛顿法等梯度方法中，使用海森矩阵的正定性可以非常便捷的判断函数是否有凸性，也就是是否可收敛到局部/全局的最优解。
仅仅使用梯度信息的优化算法称之为 一阶优化算法，如梯度下降算法等。
使用 Hessian 矩阵的优化算法称之为 二阶优化算法，如牛顿法等。
雅可比矩阵，Jacobi matrix 或者 Jacobian，是向量值函数（\(f : {R}^n → {R}^m\))的一阶偏导数按行排列所得的矩阵。
黑塞矩阵，又叫海森矩阵，Hesse matrix，是多元函数（\(f : {R}^n → R\))的二阶偏导数组成的方阵。
A为黑塞矩阵
- 当\(H\)正定矩阵时, \(f(x_1，x_2，……，x_n)\)在\(M_0(a_1，a_2，……，a_n)\)处是极小值
- 当\(H\)负定矩阵时，\(f(x_1，x_2，……，x_n)\)在\(M_0(a_1，a_2，……，a_n)\)处是极大值
- 当\(H\)不定矩阵时，\(M_0(a_1，a_2，……，a_n)\)不是极值点
- 当\(H\)为半正定矩阵或半负定矩阵时，M_0(a_1，a_2，……，a_n)是“可疑"极值点,尚需要利用其他方法来判定。
\(f1（x1，x2，……，xn)=0\)
\(f2（x1，x2，……，xn)=0\)
\(……\)
\(fm（x1，x2，……，xn)=0\)
最速下降法、牛顿法、共轭梯度法原理及对比
拟牛顿法
信赖域法
浙公网安备 33010602011771号