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

原文鏈接