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

[murdoch at stats.uwo.ca]

On 5/17/2006 11:07 AM, Martin Maechler wrote:
>>>>>> "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 ...]

Brian and I only quoted part of his message.  The patch we quoted isn't 
bad, but I'm not sure it's the best:  in particular, with the patch 
optim() returns a function value that is larger than f at the starting 
value (see below).  I think this means it would be better to change 
optim$par rather than changing optim$value to achieve consistency, but a 
quick look at optim.c made me think it would take more time than I had 
to do this without messing up something else.

I don't think I'll have a chance to look at this before 2.3.1, so if 
nobody else takes it on, I'd prefer to leave this as an unresolved bug 
report for now.

> 
> 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!

I think this is f(0.1, -0.1), so really $par should be 0.1, -0.1 in this 
case.  In the one above, it appears to have made a little progress 
before it went off track, but -0.234 is better than 2, so it should be 
returned if it's really an f(p) value.

Duncan Murdoch


>     >>>> 
>     >>>> 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