[Rd] optim "CG" bug w/patch proposal (PR#8786)
ripley at stats.ox.ac.uk
ripley at stats.ox.ac.uk
Wed May 17 18:46:45 CEST 2006
On Wed, 17 May 2006, maechler at stat.math.ethz.ch 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 ...]
The (minor) issue is that x is updated but not f(x). I think the intended
stategy was to update neither, so Andreas' patch was a change of stategy.
In particular, a question is if this should be marked as a convergence
failure. But people really need to read the reference before commenting,
and I at least need to find the time to do so in more 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
>
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>
>
--
Brian D. Ripley, ripley at stats.ox.ac.uk
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UK Fax: +44 1865 272595
More information about the R-devel
mailing list