optim "CG" bug w/patch proposal (PR#8786)

"Duncan" == Duncan Murdoch <murdoch at stats.uwo.ca>
    on Tue, 16 May 2006 08:34:06 -0400 writes:

   Duncan> On 5/16/2006 4:56 AM, westfeld at inf.tu-dresden.de
   Duncan> wrote:

   >> Probably I included too much at once in my bug report. I
   >> can live with an unfulfilled wishlist and thank you for
   >> thinking about it. The "badly-behaved" function is just
   >> an example to demonstrate the bug I reported. I think it
   >> is a bug if optim returns (without any warning) an
   >> unmatching pair of par and value: f(par) != value. And it
   >> is easily fixed.

   >>  Andreas

   Duncan> I agree with you that on return f(par) should be
   Duncan> value.  I agree with Brian that changes to the
   Duncan> underlying strategy need much more thought.

I agree (to both).
However, isn't Andreas' patch just fixing the problem
and not changing the underlying strategy at all?
[No, I did not study the code in very much detail ...]

Martin Maechler

   >> Prof Brian Ripley wrote:
   >>

   >>> [Sorry for the belated reply: this came in just as I was leaving for a
   >>> trip.]
   >>>
   >>> I've checked the original source, and the C code in optim does
   >>> accurately reflect the published algorithm.
   >>>
   >>> Since your example is a discontinuous function, I don't see why you
   >>> expect CG to work on it.  John Nash reports on his extensive
   >>> experience that method 3 is the worst, and I don't think we should let
   >>> a single 2D example of a badly-behaved function override that.
   >>>
   >>> Note that no other optim method copes with the discontiuity here: had
   >>> your reported that it would have been clear that the problem was with
   >>> the example.
   >>>
   >>> On Fri, 21 Apr 2006, westfeld at inf.tu-dresden.de wrote:
   >>>

   >>>> Dear R team,
   >>>>
   >>>> when using optim with method "CG" I got the wrong $value for the
   >>>> reported $par.
   >>>>
   >>>> Example:
   >>>> f<-function(p) {
   >>>> if (!all(p>-.7)) return(2)
   >>>> if (!all(p<.7)) return(2)
   >>>> sin((p[1])^2)*sin(p[2])
   >>>> }
   >>>> optim(c(0.1,-0.1),f,method="CG",control=list(trace=0,type=1))
   >>>> $par 19280.68 -10622.32
   >>>> $value -0.2346207 # should be 2!
   >>>>
   >>>> optim(c(0.1,-0.1),f,method="CG",control=list(trace=0,type=2))
   >>>> $par 3834.021 -2718.958
   >>>> $value -0.0009983175 # should be 2!
   >>>>
   >>>> Fix:
   >>>> --- optim.c     (Revision 37878)
   >>>> +++ optim.c     (Arbeitskopie)
   >>>> @@ -970,7 +970,8 @@
   >>>> if (!accpoint) {
   >>>> steplength *= stepredn;
   >>>> if (trace) Rprintf("*");
   >>>> -                           }
   >>>> +                           } else
   >>>> +                               *Fmin = f;
   >>>> }
   >>>> } while (!(count == n || accpoint));
   >>>> if (count < n) {
   >>>>
   >>>> After fix:
   >>>> optim(c(0.1,-0.1),f,method="CG",control=list(trace=0,type=1))
   >>>> $par 0.6993467 -0.4900145
   >>>> $value -0.2211150
   >>>> optim(c(0.1,-0.1),f,method="CG",control=list(trace=0,type=2))
   >>>> $par 3834.021 -2718.958
   >>>> $value 2
   >>>>
   >>>> Wishlist:

   >>>

   >> [wishlist deleted]
   >>
   >>

   Duncan> ______________________________________________
   Duncan> R-devel at r-project.org mailing list
   Duncan> https://stat.ethz.ch/mailman/listinfo/r-devel