【问题标题】:Implement Levenberg-Marquardt algorithm in R在 R 中实现 Levenberg-Marquardt 算法
【发布时间】:2021-03-25 20:44:37
【问题描述】:

我被告知在 R studio 中实现Levenberg-Marquardt algorithm,考虑到 lambda 的初始值等于 10。当梯度的范数低于公差时,算法必须停止。我还需要打印 x1、x2、λ、∇f(x)、d1 和 d2 每次迭代的值。关于如何做的任何想法?非常感谢提前

这就是我所拥有的:

library(pracma)
library(matlib)
MetodeLM<-function(f,xi,t)
{
  l=10
  stop=FALSE
  x<-xi
  k=0
  while (stop==FALSE){
    dk<- inv(hessian(f,x)+l*diag(diag(hessian(f,x))))
    x1<-x+dk
    if (Norm(grad(f,x1))<t){
      stop<-TRUE
    }
    else{
      if (f(x1) < f(x)){
        l<-l/10
        k<-k+1
        stop<-FALSE
      }
      else{
        l<-l*10
        stop<-FALSE
      }
    }
  }
}

【问题讨论】:

  • 欢迎来到 SO!不是每个人都知道 Levenberg Marquadt 算法是什么,所以如果不定义它,就会限制可以帮助你的人数。此外,由于没有展示您自己的任何努力(您尝试了什么?为什么没有成功?您是否收到任何错误消息?)您不会激励人们帮助您。
  • 我刚刚添加了我所拥有的,对不起

标签: r algorithm optimization mathematical-optimization nonlinear-optimization


【解决方案1】:

纠正你代码中的一些错误,下面的 Levenberg Marquadt 算法的实现应该可以工作(注意该算法的更新规则如下图所示):

library(pracma)

# tolerance = t, λ = l 
LM <- function(f, x0, t, l=10, r=10) { 
    
    x <- x0
    k <- 0
    while (TRUE) {
      H <- hessian(f, x)
      G <- grad(f, x)
      dk <- inv(H + l * diag(nrow(H))) %*% G   # dk <- solve(H + l * diag(nrow(H)), G)
      x1 <- x - dk   # update rule
      print(k)  # iteration
      # print(l) # λ
      print(x1) # x1, x2
      print(G)  # ∇f(x)
      print(dk) # d1, d2
      if (Norm(G) < t) break
      l <- ifelse(f(x1) < f(x), l / r, l * r)
      k <- k + 1
      x <- x1 # update the old point 
    }
}

例如,使用以下函数,非线性优化算法将快速找到一个局部最小值点(在第 10 次迭代中),如下所示

f <- function(x) {
   return ((x[1]^2+x[2]-25)^2 + (x[1]+x[2]^2-25)^2)
}

x0 <- rep(0,2)
LM(f, x0, t=1e-3, l=400, r=2)
# [1] 0
#      [,1]
# [1,] 0.165563
# [2,] 0.165563
# [1] -50 -50
#      [,1]
# [1,] -0.165563
# [2,] -0.165563
# [1] 1
#      [,1]
# [1,] 0.7986661
# [2,] 0.7986661
# [1] -66.04255 -66.04255
#      [,1]
# [1,] -0.6331031
# [2,] -0.6331031
# ...
# [1] 10
#     [,1]
# [1,] 4.524938
# [2,] 4.524938
# [1] 0.0001194898 0.0001194898
#         [,1]
# [1,] 5.869924e-07
# [2,] 5.869924e-07

下面的动画展示了函数收敛到局部最小点:

以下是带Log功能的

【讨论】:

  • 您很可能希望将inv(H + l * diag(nrow(H))) %*% G 替换为solve(H + l * diag(nrow(H)), G)。出于数字原因,您几乎从不想求逆和乘法而不是求解。
  • 是的,只是不想更改OP的代码,因为他正在使用pracma的函数
  • 我明白了。很好的答案。
  • 我不是 OP。我给了它一个 +1。
猜你喜欢
  • 1970-01-01
  • 1970-01-01
  • 2011-06-04
  • 1970-01-01
  • 2012-12-12
  • 2018-11-26
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
相关资源
最近更新 更多