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

原文連接