From maechler at stat.math.ethz.ch  Wed Mar  8 09:02:19 2017
From: maechler at stat.math.ethz.ch (Martin Maechler)
Date: Wed, 8 Mar 2017 09:02:19 +0100
Subject: [Rd] Bug in nlm()
In-Reply-To: <22713.42179.765118.892734@stat.math.ethz.ch>
References: <0051D1E3C4DC574698F0EDFC7FFA6ADBEBF89E8E@HERMES01.demogr.mpg.de>
	<22713.42179.765118.892734@stat.math.ethz.ch>
Message-ID: <22719.47755.541134.331752@stat.math.ethz.ch>

       {This was sent to me, MM, only, but for completeness should
        have gone back to R-devel.

	Further: I now *have* added Marie B to the members'of "R bugzilla"
	-- M.Maechler}


I had already read the R bug reporting guide and I'm sure it is a bug.
The bug occurs when the user provides not only the analytic gradient but also the analytic Hessian of the objective function. In that case, the algorithm does not converge due to an erroneous implementation of the modified Cholesky decomposition of the Hessian matrix. It is actually a bug in the C-code called by nlm(), therefore it is hard to show that the non-convergence of the algorithm is really due to this bug with only a MRE.
However, a short example (optimizing the Rosenbrock banana valley function with and without analytic Hessian) is:

fg <- function(x){  
  gr <- function(x1, x2) c(-400*x1*(x2 - x1*x1)-2*(1-x1), 200*(x2 - x1*x1)) 
  x1 <- x[1]; x2 <- x[2]
  res<- 100*(x2 - x1*x1)^2 + (1-x1)^2 
  attr(res, "gradient") <- gr(x1, x2)
  return(res)
} 
nlm.fg <- nlm(fg, c(-1.2, 1))

fgh <- function(x){ 
  gr <- function(x1, x2) c(-400*x1*(x2 - x1*x1) - 2*(1-x1), 200*(x2 - x1*x1)) 
  h <- function(x1, x2){
    a11 <- 2 - 400*x2 + 1200*x1*x1
    a21 <- -400*x1 
    matrix(c(a11, a21, a21, 200), 2, 2)
  } 
  x1 <- x[1];  x2 <- x[2];  res<- 100*(x2 - x1*x1)^2 + (1-x1)^2 
  attr(res, "gradient") <- gr(x1, x2)
  attr(res, "hessian") <- h(x1, x2) 
  return(res)
}
nlm.fgh <- nlm(fgh, c(-1.2,1))

I have almost finished a more detailed bug report, which I would like to submit.

Best,
Marie Boehnstedt

>>>>> Martin Maechler <maechler at stat.math.ethz.ch>
>>>>>     on Fri, 3 Mar 2017 18:15:47 +0100 writes:

>>>>> Boehnstedt, Marie <Boehnstedt at demogr.mpg.de>
>>>>>     on Fri, 3 Mar 2017 10:23:12 +0000 writes:

    >> Dear all, I have found a bug in nlm() and would like to
    >> submit a report on this.  Since nlm() is in the
    >> stats-package, which is maintained by the R Core team,
    >> bug reports should be submitted to R's Bugzilla. However,
    >> I'm not a member of Bugzilla. Could anyone be so kind to
    >> add me to R's Bugzilla members or let me know to whom I
    >> should send the bug report?

    > Dear Marie,

    > I can do this ... but are you really sure?  There is
    > https://www.r-project.org/bugs.html which you should spend
    > some time reading if you haven't already.

    > I think you would post a MRE (Minimal Reproducible
    > Example) here {or on stackoverflow or ...} if you'd follow
    > what the 'R bugs' web page (above) recommends and only
    > report a bug after some feedback from "the public".

    > Of course, I could be wrong.. and happy if you explain /
    > tell me why.

    > Best, Martin Maechler

    >> Thank you in advance.

    >> Kind regards, Marie B?hnstedt


    >> Marie B?hnstedt, MSc Research Scientist Max Planck
    >> Institute for Demographic Research Konrad-Zuse-Str. 1,
    >> 18057 Rostock, Germany
    >> www.demogr.mpg.de<http://www.demogr.mpg.de/>



[Rd] Bug in nlm()
Martin Maechler