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

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:

1. Please make type=3 the default in optim (it is more robust).

2. The $par reported for type=2 is still not satisfactory. I found out
that this can be improved by limiting G3 to a maximum of about 2000
(maybe even smaller). However, I'm not a "CG" expert and can live
with a
suboptimal result.

--- optim.c     (Revision 37878)
+++ optim.c     (Arbeitskopie)
@@ -946,6 +946,8 @@
                        G3 = G1 / G2;
                    else
                        G3 = 1.0;
+                   if (G3 > 2e3)
+                       G3 = 2e3;
                    gradproj = 0.0;
                    for (i = 0; i < n; i++) {
                        t[i] = t[i] * G3 - g[i];

Andreas

-- 
Andreas Westfeld, 0432 01CC F511 9E2B 0B57 5993 0B22 98F8 4AD8 EEEA
<westfeld at inf.tu-dresden.de> http://www.inf.tu-dresden.de/~aw4
TU Dresden Fakult?t Informatik, Institut f?r Systemarchitektur
Datenschutz und Datensicherheit, Tel. +49-351-463-37918

[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:
>
> 1. Please make type=3 the default in optim (it is more robust).
>
> 2. The $par reported for type=2 is still not satisfactory. I found out
> that this can be improved by limiting G3 to a maximum of about 2000
> (maybe even smaller). However, I'm not a "CG" expert and can
live with a
> suboptimal result.
>
> --- optim.c     (Revision 37878)
> +++ optim.c     (Arbeitskopie)
> @@ -946,6 +946,8 @@
>                        G3 = G1 / G2;
>                    else
>                        G3 = 1.0;
> +                   if (G3 > 2e3)
> +                       G3 = 2e3;
>                    gradproj = 0.0;
>                    for (i = 0; i < n; i++) {
>                        t[i] = t[i] * G3 - g[i];
>
> Andreas
>
>
-- 
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

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

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]


-- 
Andreas Westfeld, 0432 01CC F511 9E2B 0B57 5993 0B22 98F8 4AD8 EEEA
<westfeld at inf.tu-dresden.de> http://www.inf.tu-dresden.de/~aw4
TU Dresden Fakult?t Informatik, Institut f?r Systemarchitektur
Datenschutz und Datensicherheit, Tel. +49-351-463-37918

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

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

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