title: how to propery specify a gradient function for use in optim() or other optimizer
body: <p>I have an optimization problem that the <code>Nelder-Mead</code> method will solve, but that I would also like to solve using <code>BFGS</code> or Newton-Raphson, or something that takes a gradient function, for more speed, and hopefully more precise estimates. I wrote such a gradient function following (I thought) the example in the <code>optim</code> / <code>optimx</code> documentation, but when I use it with <code>BFGS</code> my starting values either don't move (<code>optim()</code>), or else the function outright doesn't run (<code>optimx()</code>, which returns <code>Error: Gradient function might be wrong - check it!</code>). I'm sorry there's a bit of code involved in reproducing this, but here goes:</p>
  
  <p>This is the function that I want to get parameter estimates for (this is for smoothing old-age mortality rates, where x is age, starting at age 80):</p>
  
  <pre><code>    KannistoMu &lt;- function(pars, x = .5:30.5){
        a &lt;- pars["a"]
        b &lt;- pars["b"]
        (a * exp(b * x)) / (1 + a * exp(b * x))
      }
  </code></pre>
  
  <p>And here's a log likelihood function for estimating it from observed rates (defined as deaths, <code>.Dx</code> over exposure, <code>.Exp</code>):</p>
  
  <pre><code>    KannistoLik1 &lt;- function(pars, .Dx, .Exp, .x. = .5:30.5){
        mu &lt;- KannistoMu(exp(pars), x = .x.)
        # take negative and minimize it (default optimizer behavior)
        -sum(.Dx * log(mu) - .Exp * mu, na.rm = TRUE) 
      }
  </code></pre>
  
  <p>you see <code>exp(pars)</code> in there because I give <code>log(pars)</code> to optimize over, in order to constrain the final <code>a</code> and <code>b</code> to be positive.</p>
  
  <p>Example data (1962 Japan females, if anyone is curious):</p>
  
  <pre><code>    .Dx &lt;- structure(c(10036.12, 9629.12, 8810.11, 8556.1, 7593.1, 6975.08, 
        6045.08, 4980.06, 4246.06, 3334.04, 2416.03, 1676.02, 1327.02, 
        980.02, 709, 432, 350, 217, 134, 56, 24, 21, 10, 8, 3, 1, 2, 
        1, 0, 0, 0), .Names = c("80", "81", "82", "83", "84", "85", "86", 
        "87", "88", "89", "90", "91", "92", "93", "94", "95", "96", "97", 
        "98", "99", "100", "101", "102", "103", "104", "105", "106", 
        "107", "108", "109", "110"))
      .Exp &lt;- structure(c(85476.0333333333, 74002.0866666667, 63027.5183333333, 
        53756.8983333333, 44270.9, 36749.85, 29024.9333333333, 21811.07, 
        16912.315, 11917.9583333333, 7899.33833333333, 5417.67, 3743.67833333333, 
        2722.435, 1758.95, 1043.985, 705.49, 443.818333333333, 223.828333333333, 
        93.8233333333333, 53.1566666666667, 27.3333333333333, 16.1666666666667, 
        10.5, 4.33333333333333, 3.16666666666667, 3, 2.16666666666667, 
        1.5, 0, 1), .Names = c("80", "81", "82", "83", "84", "85", "86", 
        "87", "88", "89", "90", "91", "92", "93", "94", "95", "96", "97", 
        "98", "99", "100", "101", "102", "103", "104", "105", "106", 
        "107", "108", "109", "110"))
  </code></pre>
  
  <p>The following works for the <code>Nelder-Mead</code> method:</p>
  
  <pre><code>    NMab &lt;- optim(log(c(a = .1, b = .1)), 
        fn = KannistoLik1, method = "Nelder-Mead",
        .Dx = .Dx, .Exp = .Exp)
      exp(NMab$par) 
      # these are reasonable estimates
         a         b 
      0.1243144 0.1163926
  </code></pre>
  
  <p>This is the gradient function I came up with:</p>
  
  <pre><code>    Kannisto.gr &lt;- function(pars, .Dx, .Exp, x = .5:30.5){
        a &lt;- exp(pars["a"])
        b &lt;- exp(pars["b"])
        d.a &lt;- (a * exp(b * x) * .Exp + (-a * exp(b * x) - 1) * .Dx) /
          (a ^ 3 * exp(2 * b * x) + 2 * a ^ 2 * exp(b * x) + a)
        d.b &lt;- (a * x * exp(b * x) * .Exp + (-a * x * exp(b * x) - x) * .Dx) /
          (a ^ 2 * exp(2 * b * x) + 2 * a * exp(b * x) + 1)
        -colSums(cbind(a = d.a, b = d.b), na.rm = TRUE)
      }
  </code></pre>
  
  <p>The output is a vector of length 2, the change with respect to the parameters <code>a</code> and <code>b</code>. I also have an uglier version arrived at by exploiting the output of <code>deriv()</code>, which returns the same answer, and which I don't post (just to confirm that the derivatives are right). </p>
  
  <p>If I supply it to <code>optim()</code> as follows, with <code>BFGS</code> as the method, the estimates do not move from the starting values:</p>
  
  <pre><code>    BFGSab &lt;- optim(log(c(a = .1, b = .1)), 
        fn = KannistoLik1, gr = Kannisto.gr, method = "BFGS",
        .Dx = .Dx, .Exp = .Exp)
      # estimates do not change from starting values:
      exp(BFGSab$par) 
        a   b 
      0.1 0.1
  </code></pre>
  
  <p>When I look at the <code>$counts</code> element of the output, it says that <code>KannistoLik1()</code> was called 31 times and <code>Kannisto.gr()</code> just 1 time. <code>$convergence</code> is <code>0</code>, so I guess it thinks it converged (if I give less reasonable starts they also stay put). I reduced the tolerance, etc, and nothing changes. When I try the same call in <code>optimx()</code> (not shown), I receive the waring I mentioned above, and no object is returned. I get the same results when specifying <code>gr = Kannisto.gr</code> with the <code>"CG"</code>. With the <code>"L-BFGS-B"</code> method I get the same starting values back as estimate, but it is also reported that both function and gradient were called 21 times, and there is an error message:
   <code>"ERROR: BNORMAL_TERMINATION_IN_LNSRCH"</code></p>
  
  <p>I'm hoping that there is some minor detail in the way the gradient function is written that will solve this, as this later warning and the <code>optimx</code> behavior are bluntly hinting that the function simply isn't right (I think). I also tried the <code>maxNR()</code> maximizer from the <code>maxLik</code> package and observed similar behavior (starting values don't move). Can anyone give me a pointer? Much obliged</p>
  
  <p>[Edit]
  @Vincent suggested I compare with the output from a numerical approximation:</p>
  
  <pre><code>    library(numDeriv)
      grad( function(u) KannistoLik1( c(a=u[1], b=u[2]), .Dx, .Exp ), log(c(.1,.1)) )
      [1] -14477.40  -7458.34
      Kannisto.gr(log(c(a=.1,b=.1)), .Dx, .Exp)
       a        b 
      144774.0  74583.4 
  </code></pre>
  
  <p>so different sign, and off by a factor of 10? I change the gradient function to follow suit:</p>
  
  <pre><code>    Kannisto.gr2 &lt;- function(pars, .Dx, .Exp, x = .5:30.5){
        a &lt;- exp(pars["a"])
        b &lt;- exp(pars["b"])
        d.a &lt;- (a * exp(b * x) * .Exp + (-a * exp(b * x) - 1) * .Dx) /
          (a ^ 3 * exp(2 * b * x) + 2 * a ^ 2 * exp(b * x) + a)
        d.b &lt;- (a * x * exp(b * x) * .Exp + (-a * x * exp(b * x) - x) * .Dx) /
          (a ^ 2 * exp(2 * b * x) + 2 * a * exp(b * x) + 1)
        colSums(cbind(a=d.a,b=d.b), na.rm = TRUE) / 10
      }
      Kannisto.gr2(log(c(a=.1,b=.1)), .Dx, .Exp)
      # same as numerical:
        a         b 
      -14477.40  -7458.34 
  </code></pre>
  
  <p>Try it in the optimizer:</p>
  
  <pre><code>    BFGSab &lt;- optim(log(c(a = .1, b = .1)), 
        fn = KannistoLik1, gr = Kannisto.gr2, method = "BFGS",
        .Dx = .Dx, .Exp = .Exp)
      # not reasonable results:
      exp(BFGSab$par) 
        a   b 
      Inf Inf 
      # and in fact, when not exp()'d, they look oddly familiar:
      BFGSab$par
        a         b 
      -14477.40  -7458.34 
  </code></pre>
  
  <p>Following Vincent's answer, I rescaled the gradient function, and used <code>abs()</code> instead of <code>exp()</code> to keep parameters positive. The most recent, and better performing objective and gradient functions:</p>
  
  <pre><code>    KannistoLik2 &lt;- function(pars, .Dx, .Exp, .x. = .5:30.5){
        mu &lt;- KannistoMu.c(abs(pars), x = .x.)
        # take negative and minimize it (default optimizer behavior)
        -sum(.Dx * log(mu) - .Exp * mu, na.rm = TRUE) 
      }
  
      # gradient, to be down-scaled in `optim()` call
      Kannisto.gr3 &lt;- function(pars, .Dx, .Exp, x = .5:30.5){
        a &lt;- abs(pars["a"])
        b &lt;- abs(pars["b"])
        d.a &lt;- (a * exp(b * x) * .Exp + (-a * exp(b * x) - 1) * .Dx) /
          (a ^ 3 * exp(2 * b * x) + 2 * a ^ 2 * exp(b * x) + a)
        d.b &lt;- (a * x * exp(b * x) * .Exp + (-a * x * exp(b * x) - x) * .Dx) /
          (a ^ 2 * exp(2 * b * x) + 2 * a * exp(b * x) + 1)
        colSums(cbind(a = d.a, b = d.b), na.rm = TRUE) 
      }
  
      # try it out:
      BFGSab2 &lt;- optim(
        c(a = .1, b = .1), 
        fn = KannistoLik2, 
        gr = function(...) Kannisto.gr3(...) * 1e-7, 
        method = "BFGS",
        .Dx = .Dx, .Exp = .Exp
      )
      # reasonable:
      BFGSab2$par
              a         b 
      0.1243249 0.1163924 
  
      # better:
      KannistoLik2(exp(NMab1$par),.Dx = .Dx, .Exp = .Exp) &gt; KannistoLik2(BFGSab2$par,.Dx = .Dx, .Exp = .Exp)
      [1] TRUE
  </code></pre>
  
  <p>This was solved much faster than I was expecting, and I learned more than a couple tricks. Thanks Vincent!</p>