R devel - May 2006 - optim "CG" bug w/patch proposal (PR#8786)

On Wed, 17 May 2006, Prof Brian Ripley wrote:
> 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.
Having spent some time with the reference and a debugger, as far as I can 
see what is happening here is that the second phase of the line search 
fails. That can only happen (in exactly this way) for a discontinuous 
function where the first phase has jumped over a discontinuity to an 
essentially flat region.  In those circumstances updating *Fmin seems 
sensible: probably ideally one should detect the lack of near-quadratic 
and force a restart.  But I don't think it is worth much effort to detect 
examples that are a very long way from the assumptions.  So we'll just 
update *Fmin.
>
>> 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