R 语言实现牛顿下降法

凸是一个很好的性质.如果已经证明了某个问题是凸的,那这个问题基本上算是解决了.
最近在解决一个多目标优化的问题.多目标的问题往往是非凸的.好在能够知道这个问题的近似解大概是多少.这样这个多目标优化的问题至少能够在局部运用凸优化的方法来解决了.解决凸优化的方法有很多,比如梯度下降法,内点法.在梯度下降法中,牛顿下降法是一种重要的方法,也容易实现.更好的是牛顿下降法的收敛速度是二次的,比通常的下降法的收敛速度要快很多.

牛顿算法
$ x(n+1) = x(n) - H(x(n))^{-1} grad f(x(n)) \( \)H(x)表示hessian矩阵$

从这里可以看出,如果hessian矩阵是奇异的,那么牛顿下降法将会失效.这是后就需要运用其他的算法了.比如拟牛顿法.

R语言实现(代码)

newton <- function(func = objfun, x0, tol = 1e-5, n.max = 100,...){
    x <- x0
    g <- grad(func, x, ...)
    h <- hessian(func, x, ...)

    n <- 0
    while( max(abs(g))>tol && n<n.max ){
        x <- x-solve(h,g)
        g <- grad(func, x, ...)
        h <- hessian(func, x, ...)
        n <- n+1
    }
    if(n == n.max){
        cat('newton failed to converge\n')
        return(x)
    }
    return(x)
}

依赖包

numDeriv
如果需要安装,在R控制台里键入:

install.package(numDeriv)

说明

func : 目标函数.
x0: 目标函数的极小化初始值.
tol:控制精度,越接近零越精确,代表梯度已经是0.
n.max:迭代次数.
... : 目标函数的其他参数.

例子

testfun <- function(x, a,b){
    y <- a*x[1]^2 + b*x[2]^2
    return (y)
}
library(numDeriv)
y <- newton(func = testfun, x0=c(1,1), a = 1,b = 1)

输出

y
[1] 3.596345e-12 3.596345e-12

原文链接