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

[ripley at stats.ox.ac.uk]

[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