【问题标题】:Optimization issue, Nonlinear: automatic analytical Jacobian/Hessian from objective and constraints in R?优化问题,非线性:从 R 中的目标和约束自动分析 Jacobian/Hessian?
【发布时间】:2014-02-21 05:30:24
【问题描述】:

在 R 中,当您只为优化问题提供目标函数和约束时,是否可以通过解析方式找到 Jacobian/Hessian/稀疏模式?

AMPL 可以做到这一点,据我所知,甚至 MATLAB 也可以做到这一点,但我不知道您是否需要 Knitro。

但是,R 的所有优化工具(例如 nloptr)似乎都需要 me 自己输入梯度和 Hessian,这非常困难,因为我正在处理复杂的模型。 p>

【问题讨论】:

  • optim 中有一个hessian 参数。这不是你想要的吗?
  • 粗麻布是用数值计算的,这对我正在做的工作类型来说是个问题。

标签: r automatic-differentiation


【解决方案1】:

您要查找的内容称为automatic differentiation。可悲的是,在我看来它在 R 中不可用。

attempts 5 years ago 来实施它,但我的简短调查表明这些尝试已经消失。

有一个相当新的 R 包 (Automatic Differentiation Model Builder),但我不清楚如何使用它,或者如何将其应用于您的情况。 (我自己不使用 R,所以我不知道。)

【讨论】:

  • 我刚刚用 R 写了一个广告包:github.com/shabbychef/madness。如果您有任何问题,请提交问题。
  • 请注意阅读本文的人,显然 Julia(语言)可以做到这一点。我将为此开始部分切换
  • @robertevansanders 其他语言中存在几个自动微分库,但该问题专门要求在 R 中提供解决方案。显然,上面链接的 shabbychef / madness R 包也可以做到这一点,但它在 2014 年不存在当我写下我的答案时。
【解决方案2】:

1)optim 中默认的 Nelder Mead 方法不需要导数,也不在内部计算它们。

2) Dderiv 和相关的 R 函数(参见 ?deriv)可以计算简单的符号导数。

3) Ryacas package 可以计算符号导数。

【讨论】:

    【解决方案3】:

    我为 R 中的自动微分编写了一个基本包,名为 madness。虽然它主要是为多元 delta 方法设计的,但它也应该可用于自动计算梯度。示例用法:

    require(madness)
    
    # the 'fit' is the Frobenius norm of Y - L*R
    # with a penalty for strongly negative R.
    compute_fit <- function(R,L,Y) {
        Rmad <- madness(R)
        Err <- Y - L %*% Rmad
        penalty <- sum(exp(-0.1 * Rmad))
        fit <- norm(Err,'f') + penalty
    }
    
    set.seed(1234)
    R <- array(runif(5*20),dim=c(5,20))
    L <- array(runif(1000*5),dim=c(1000,5))
    Y <- array(runif(1000*20),dim=c(1000,20))
    ftv <- compute_fit(R,L,Y)
    show(ftv)
    show(val(ftv))
    show(dvdx(ftv))
    

    我得到以下信息:

    class: madness 
            d (norm((numeric - (numeric  %*% R)), 'f') + (sum(exp((numeric * R)), na.rm=FALSE) + numeric))
     calc: ------------------------------------------------------------------------------------------------ 
                                                          d R
      val: 207.6 ...
     dvdx: 4.214 4.293 4.493 4.422 4.672 2.899 2.222 2.565 2.854 2.718 4.499 4.055 4.161 4.473 4.069 2.467 1.918 2.008 1.874 1.942 0.2713 0.2199 0.135 0.03017 0.1877 5.442 4.81
    1 5.472 5.251 4.674 1.933 1.62 1.79 1.902 1.665 5.232 4.435 4.789 5.183 5.084 3.602 3.477 3.419 3.592 3.376 4.109 3.937 4.017 3.816 4.2 1.776 1.784 2.17 1.975 1.699 4.365 4
    .09 4.475 3.964 4.506 1.745 1.042 1.349 1.084 1.237 3.1 2.575 2.887 2.524 2.902 2.055 2.441 1.959 2.467 1.985 2.494 2.223 2.124 2.275 2.546 3.497 2.961 2.897 3.111 2.9 4.44
    2 3.752 3.939 3.694 4.326 0.9582 1.4 0.8971 0.8756 0.9019 2.399 2.084 2.005 2.154 2.491 ...
     varx:  ...
    
          [,1]
    [1,] 207.6
    
          [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20]  [,21]  [,22] [,23]   [,24]  [,25] [,26] [,27]
    [1,] 4.214 4.293 4.493 4.422 4.672 2.899 2.222 2.565 2.854 2.718 4.499 4.055 4.161 4.473 4.069 2.467 1.918 2.008 1.874 1.942 0.2713 0.2199 0.135 0.03017 0.1877 5.442 4.811
         [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48] [,49] [,50] [,51] [,52] [,53] [,54]
    [1,] 5.472 5.251 4.674 1.933  1.62  1.79 1.902 1.665 5.232 4.435 4.789 5.183 5.084 3.602 3.477 3.419 3.592 3.376 4.109 3.937 4.017 3.816   4.2 1.776 1.784  2.17 1.975
         [,55] [,56] [,57] [,58] [,59] [,60] [,61] [,62] [,63] [,64] [,65] [,66] [,67] [,68] [,69] [,70] [,71] [,72] [,73] [,74] [,75] [,76] [,77] [,78] [,79] [,80] [,81]
    [1,] 1.699 4.365  4.09 4.475 3.964 4.506 1.745 1.042 1.349 1.084 1.237   3.1 2.575 2.887 2.524 2.902 2.055 2.441 1.959 2.467 1.985 2.494 2.223 2.124 2.275 2.546 3.497
         [,82] [,83] [,84] [,85] [,86] [,87] [,88] [,89] [,90]  [,91] [,92]  [,93]  [,94]  [,95] [,96] [,97] [,98] [,99] [,100]
    [1,] 2.961 2.897 3.111   2.9 4.442 3.752 3.939 3.694 4.326 0.9582   1.4 0.8971 0.8756 0.9019 2.399 2.084 2.005 2.154  2.491
    

    请注意,导数是标量相对于 5 x 20 矩阵的导数,但在这里被展平为向量。 (不幸的是一个行向量。)

    【讨论】:

      【解决方案4】:

      看看solnp,包Rsolnp。它是一个非线性规划求解器,不需要解析 Jacobian 或 Hessian:

      min f(x)
      s.t. g(x) = 0
      l[h] <= h(x) <= u[h]
      l[x] <= x <= u[x]
      

      【讨论】:

      • 如果它以数字方式计算雅可比矩阵和黑森矩阵,你知道吗?还是它能够像数学编程语言那样理解优化问题并找到它可以使用的解析导数?
      • 当然,它使用有限差分类型的数值近似。顺便说一下,后者(符号微分)在 R 中也是 available
      • 有限差异是我正在做的工作中的一个大问题。计算负担非常高,如果要在合理的时间内完成,我需要为求解器提供解析导数。您是说我必须始终自己在 R 中编写派生程序吗? (不像 MATLAB 和 AMPL?)
      • 我不能说,总是;在大多数情况下,是的。这是comprehensive list,您可能会在那里找到有用的东西。如果您的搜索成功,请告诉我。
      猜你喜欢
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 2014-03-13
      • 2013-02-15
      • 1970-01-01
      • 1970-01-01
      • 2018-03-15
      相关资源
      最近更新 更多