R语言实现牛顿迭代算法

• 2019 年 11 月 4 日
• 笔记

我们今天给大家介绍一个用来迭代的算法牛顿迭代法（Newton's method）。单变量下又称为切线法。它是一种在实数域和复数域上近似求解方程的方法。首先我们看下牛顿迭代算法的公式：

其中，Xn+1和Xn主要是指的在n+1和n这个位置的X值。

以上公式的推导可以用泰勒展开公式进行推导。当然上面只是一元一阶的情况的，如果更复杂的那就是多元多阶的情况需要用到更复杂的Hessian矩阵获得以上的通用公式模型。

关于推导的部分，我们就不展开了。接下来我们直接用一个R语言的实例来看下，牛顿迭代是如何工作的。我们看下下面这个例题：

以上就是简单的一元函数求解，当然我们基于我们数学的基础也可以人工展开计算，但是当次幂升到很高，那我们就无从下手了，这时候就可以直接通过牛顿迭代进行获取根。我们就基于上面的例子进行程序设计：

funs=function(x){    f=x    J=(2*x^3-4*x^2+3*x-6)/(6*x^2-8*x+3)    list(f=f,J=J);                              }    #Newton迭代法  Newtons=function(fun,x,ep=1e-5,it_max=1000){   index=0;k=1   while(k<=it_max){   norm=abs(fun(x)$J)   x=x-fun(x)$J         if(norm<ep){ index=1;break}         k=k+1        }     list(root=x,it=k,index=index)  }  Newtons(funs,1)

通过上面的程序计算就得到我们的根：

上面root就是我们得到的根，it指的迭代的次数，index指的最后的结果1代表找到根；0代表没找到根。

上面是一个简单的例子，那么如果多元的函数我们如何去计算它的根呢，那么就需要我们的Hessian矩阵方法帮我们进行一步转化，然后再生成我们想要的函数。接下来看一下实例：

上面这个例子呢，其实是已经帮大家转化好了对应的矩阵方程，一般我们遇到的那可就不这么简单了，可能是上面的一个方程式组，更可能是只有一个函数方程，需要我们进行下一步的转化，一直到最后转化为我们想要的矩阵方程式。然后构建我们的程序（以下程序引自网友PPT）：

funs=function(x){    f=c(x[1]^2+x[2]^2-5,(x[1]+1)*x[2]-(3*x[1]+1));    J=matrix(c(2*x[1],2*x[2],x[2]-3,x[1]+1),nrow=2,byrow=T);    list(f=f,J=J);                              }    #Newton迭代法  Newtons=function(fun,x,ep=1e-5,it_max=100){   index=0;k=1   while(k<=it_max){   x1=x;obj=fun(x);x=x-solve(obj$J,obj$f);norm=sqrt((x-x1)%*%(x-x1));         if(norm<ep){ index=1;break}         k=k+1        }     obj=fun(x)     list(root=x,it=k,index=index,Funval=obj$f)  }    Newtons(funs,c(0,1))

上图结果多了个Funval，也就是对应的根的f(x)的值。由结果可以看出，的确可以迭代到非常接近根的位置。

当然还有其他的迭代算法梯度下降法、拟牛顿法，三者并称是机器学习中最常见的三大类迭代法。

具体在真实世界的应用，大家可以去探索发现。