[R] Relatively Simple Maximization Using Optim Doesnt Optimize

Duncan Murdoch murdoch@dunc@n @end|ng |rom gm@||@com
Fri Mar 13 11:47:05 CET 2020


On 12/03/2020 8:52 p.m., Abby Spurdle wrote:
>> L= 1006.536
>> L= 1006.537
>> L= 1006.535
>> It appears to have chosen step size 0.001, not 0.00001.  It should be
>> getting adequate accuracy in both 1st and 2nd derivatives.
>> Those little ripples you see in the plot are not relevant.
> 
> I'm impressed.
> But you're still wrong.
> 
> Try this:
> ---------
> #not good R code!
> v = numeric ()
> 
> production3 <- function(L){
>      #store in vector
>      v <<- c (v, L)
> 
>     budget=100000
>     Lcost=12
>     Kcost=15
>     K=(budget-L*Lcost)/Kcost
>     machines=0.05*L^(2/3)*K^(1/3)
>     return(machines)
> }
> 
> optim.sol <- optim (1001, production3 ,method="CG", control = list(fnscale=-1) )
> 
> n = length (v)
> print (n)
> 
> plot (1:n ,v, type="l")
> ---------
> 
> After 401 iterations (on my computer), the algorithm hasn't converged.
> And I note it's converging extremely slowly, so I don't see any
> argument for increasing the number of iterations.
> 
> And try this:
> (The first 30 steps).
> ---------
> plot (1:30 ,v [1:30], type="l")
> ---------
> 
> Little ripples aren't going anywhere...
> 

That's the bug.

It is correctly signalling that it hasn't converged (look at 
optim.sol$convergence, which "indicates that the iteration limit maxit 
had been reached".)  But CG should be taking bigger steps.  On a 1D 
quadratic objective function with no errors in the derivatives, it 
should take one step to convergence.  Here we're not quadratic (though 
it's pretty close), and we don't have exact derivatives (your ripples), 
so the fact that it is sticking to one step size is a sign that it is 
not working.   If those ripples are big enough to matter (and I'm not 
convinced of that), it should take highly variable steps.

The fact that it doesn't give a warning() when it knows it has failed to 
converge is also a pretty serious design flaw.

Duncan Murdoch



More information about the R-help mailing list