R help - Mar 2020 - Relatively Simple Maximization Using Optim Doesnt Optimize

I am trying to familiarize myself with optim() with a relatively simple
maximization.

Description:
L and K are two terms which are constrained to add up to a total 100000
(with respective weights to each). To map this constraint I plugged K into
the function (to make this as simple as possible.)

Together these two feed into one nonlinear function which is the product of
two monotonic (on the positive interval) functions. Then that numbers is
returned in a function fed to optim, which should maximize the output by
adjusting L. The whole code is:

production1 <- function(L){
  budget=100000
  Lcost=12
  Kcost=15
  K=(budget-L*Lcost)/Kcost
  machines=0.05*L^(2/3)*K^(1/3)
  return(machines)
}

# production1(6000) #example of number with much higher output vs optim
result
S1=optim(1001,production1,method="CG",control=list(fnscale=-1))
S1

Output:
$par
[1] 1006.536

$value
[1] 90.54671

$counts
function gradient
     201      101

$convergence
[1] 1

$message
NULL


For some reason this never explores the problem space and just spits out
some answer close to the initial condition. What am I doing wrong?

Thanks,
Skyler S.

	[[alternative HTML version deleted]]

The help file points out that CG is "fragile" ... and I would expect
that failing to define a gradient function will exacerbate that.

I think you should use a different algorithm or specify a gradient function. You
might also consider working with the more recent optimr package contributed by
Dr Nash, author of the original optim function in R.

On March 12, 2020 2:30:26 AM PDT, Skyler Saleebyan <skylerbsaleebyan at
gmail.com> wrote:
>I am trying to familiarize myself with optim() with a relatively simple
>maximization.
>
>Description:
>L and K are two terms which are constrained to add up to a total 100000
>(with respective weights to each). To map this constraint I plugged K
>into
>the function (to make this as simple as possible.)
>
>Together these two feed into one nonlinear function which is the
>product of
>two monotonic (on the positive interval) functions. Then that numbers
>is
>returned in a function fed to optim, which should maximize the output
>by
>adjusting L. The whole code is:
>
>production1 <- function(L){
>  budget=100000
>  Lcost=12
>  Kcost=15
>  K=(budget-L*Lcost)/Kcost
>  machines=0.05*L^(2/3)*K^(1/3)
>  return(machines)
>}
>
># production1(6000) #example of number with much higher output vs optim
>result
>S1=optim(1001,production1,method="CG",control=list(fnscale=-1))
>S1
>
>Output:
>$par
>[1] 1006.536
>
>$value
>[1] 90.54671
>
>$counts
>function gradient
>     201      101
>
>$convergence
>[1] 1
>
>$message
>NULL
>
>
>For some reason this never explores the problem space and just spits
>out
>some answer close to the initial condition. What am I doing wrong?
>
>Thanks,
>Skyler S.
>
>	[[alternative HTML version deleted]]
>
>______________________________________________
>R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
>https://stat.ethz.ch/mailman/listinfo/r-help
>PLEASE do read the posting guide
>http://www.R-project.org/posting-guide.html
>and provide commented, minimal, self-contained, reproducible code.
-- 
Sent from my phone. Please excuse my brevity.

It is possible to work out this problem explicitly. Playing with a few
different calls to optim shows that the method="L-BFGS-B" gives the
correct
answer.
I don't have particular insight into why method="CG" is
problematic.


On Thu, Mar 12, 2020 at 4:12 PM Jeff Newmiller <jdnewmil at
dcn.davis.ca.us>
wrote:
> The help file points out that CG is "fragile" ... and I would
expect that
> failing to define a gradient function will exacerbate that.
>
> I think you should use a different algorithm or specify a gradient
> function. You might also consider working with the more recent optimr
> package contributed by Dr Nash, author of the original optim function in R.
>
> On March 12, 2020 2:30:26 AM PDT, Skyler Saleebyan <
> skylerbsaleebyan at gmail.com> wrote:
> >I am trying to familiarize myself with optim() with a relatively simple
> >maximization.
> >
> >Description:
> >L and K are two terms which are constrained to add up to a total 100000
> >(with respective weights to each). To map this constraint I plugged K
> >into
> >the function (to make this as simple as possible.)
> >
> >Together these two feed into one nonlinear function which is the
> >product of
> >two monotonic (on the positive interval) functions. Then that numbers
> >is
> >returned in a function fed to optim, which should maximize the output
> >by
> >adjusting L. The whole code is:
> >
> >production1 <- function(L){
> >  budget=100000
> >  Lcost=12
> >  Kcost=15
> >  K=(budget-L*Lcost)/Kcost
> >  machines=0.05*L^(2/3)*K^(1/3)
> >  return(machines)
> >}
> >
> ># production1(6000) #example of number with much higher output vs optim
> >result
>
>S1=optim(1001,production1,method="CG",control=list(fnscale=-1))
> >S1
> >
> >Output:
> >$par
> >[1] 1006.536
> >
> >$value
> >[1] 90.54671
> >
> >$counts
> >function gradient
> >     201      101
> >
> >$convergence
> >[1] 1
> >
> >$message
> >NULL
> >
> >
> >For some reason this never explores the problem space and just spits
> >out
> >some answer close to the initial condition. What am I doing wrong?
> >
> >Thanks,
> >Skyler S.
> >
> >       [[alternative HTML version deleted]]
> >
> >______________________________________________
> >R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> >https://stat.ethz.ch/mailman/listinfo/r-help
> >PLEASE do read the posting guide
> >http://www.R-project.org/posting-guide.html
> >and provide commented, minimal, self-contained, reproducible code.
>
> --
> Sent from my phone. Please excuse my brevity.
>
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
	[[alternative HTML version deleted]]

It looks like a bug in the CG method.  The other methods in optim() all 
work fine.  CG is documented to be a good choice in high dimensions; why 
did you choose it for a 1 dim problem?

Duncan Murdoch

On 12/03/2020 2:30 a.m., Skyler Saleebyan wrote:
> I am trying to familiarize myself with optim() with a relatively simple
> maximization.
> 
> Description:
> L and K are two terms which are constrained to add up to a total 100000
> (with respective weights to each). To map this constraint I plugged K into
> the function (to make this as simple as possible.)
> 
> Together these two feed into one nonlinear function which is the product of
> two monotonic (on the positive interval) functions. Then that numbers is
> returned in a function fed to optim, which should maximize the output by
> adjusting L. The whole code is:
> 
> production1 <- function(L){
>    budget=100000
>    Lcost=12
>    Kcost=15
>    K=(budget-L*Lcost)/Kcost
>    machines=0.05*L^(2/3)*K^(1/3)
>    return(machines)
> }
> 
> # production1(6000) #example of number with much higher output vs optim
> result
> S1=optim(1001,production1,method="CG",control=list(fnscale=-1))
> S1
> 
> Output:
> $par
> [1] 1006.536
> 
> $value
> [1] 90.54671
> 
> $counts
> function gradient
>       201      101
> 
> $convergence
> [1] 1
> 
> $message
> NULL
> 
> 
> For some reason this never explores the problem space and just spits out
> some answer close to the initial condition. What am I doing wrong?
> 
> Thanks,
> Skyler S.
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

As author of CG (at least the code that was used to build it), I can say I was
never happy with that code. Rcgmin is the replacement I wrote, and I believe
that
could still be improved.

BUT:
  - you have a 1D optimization. Use Brent method and supply bounds.
  - I never intended CG (or BFGS or Nelder-Mead or ...) to work for 1D
    problems
  - as Jeff points out, you need the gradient. I stop Rcgmin and Rvmmin if
    user hasn't supplied one, as numerical approximations need to be very
    good for these gradient methods

JN

On 2020-03-12 10:03 a.m., Jeff Newmiller wrote:
> The help file points out that CG is "fragile" ... and I would
expect that failing to define a gradient function will exacerbate that.
> 
> I think you should use a different algorithm or specify a gradient
function. You might also consider working with the more recent optimr package
contributed by Dr Nash, author of the original optim function in R.
> 
> On March 12, 2020 2:30:26 AM PDT, Skyler Saleebyan <skylerbsaleebyan at
gmail.com> wrote:
>> I am trying to familiarize myself with optim() with a relatively simple
>> maximization.
>>
>> Description:
>> L and K are two terms which are constrained to add up to a total 100000
>> (with respective weights to each). To map this constraint I plugged K
>> into
>> the function (to make this as simple as possible.)
>>
>> Together these two feed into one nonlinear function which is the
>> product of
>> two monotonic (on the positive interval) functions. Then that numbers
>> is
>> returned in a function fed to optim, which should maximize the output
>> by
>> adjusting L. The whole code is:
>>
>> production1 <- function(L){
>>  budget=100000
>>  Lcost=12
>>  Kcost=15
>>  K=(budget-L*Lcost)/Kcost
>>  machines=0.05*L^(2/3)*K^(1/3)
>>  return(machines)
>> }
>>
>> # production1(6000) #example of number with much higher output vs optim
>> result
>>
S1=optim(1001,production1,method="CG",control=list(fnscale=-1))
>> S1
>>
>> Output:
>> $par
>> [1] 1006.536
>>
>> $value
>> [1] 90.54671
>>
>> $counts
>> function gradient
>>     201      101
>>
>> $convergence
>> [1] 1
>>
>> $message
>> NULL
>>
>>
>> For some reason this never explores the problem space and just spits
>> out
>> some answer close to the initial condition. What am I doing wrong?
>>
>> Thanks,
>> Skyler S.
>>
>> 	[[alternative HTML version deleted]]
>>
>> ______________________________________________
>> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>

I'm sorry, Duncan.
But I disagree.

This is not a "bug" in optim function, as such.
(Or at least, there's nothing in this discussion to suggest that there's
a bug).
But rather a floating point arithmetic related problem.

The OP's function looks simple enough, at first glance.
But it's not.

Plotting a numerical approximation of the derivative, makes the
problem more apparent:
----------
plot_derivative <- function (f, a = sol - offset, b = sol + offset,
sol, offset=0.001, N=200)
{   FIRST <- 1:(N - 2)
    LAST <- 3:N
    MP <- 2:(N - 1)

    x <- seq (a, b, length.out=N)
    y <- f (x)
    dy <- (y [LAST] - y [FIRST]) / (x [LAST] - x [FIRST])

    plot (x [MP], dy, type="l", xlab="x", ylab="dy/dx
(approx)")
}

optim.sol <- optim (1001, production1 ,method="CG", control = list
(fnscale=-1) )$par
plot_derivative (production1, sol=optim.sol)
abline (v=optim.sol, lty=2, col="grey")
----------

So, I would say the optim function (including the CG method) is doing
what it's supposed to do.

And collating/expanding on Nash's, Jeff's and Eric's comments:
(1) An exact solution can be derived quickly, so using a numerical
method is unnecessary, and inefficient.
(2) Possible problems with the CG method are noted in the documentation.
(3) Numerical approximations of the function's derivative need to be
well-behaved for gradient-based numerical methods to work properly.


On Fri, Mar 13, 2020 at 3:42 AM Duncan Murdoch <murdoch.duncan at
gmail.com> wrote:
>
> It looks like a bug in the CG method.  The other methods in optim() all
> work fine.  CG is documented to be a good choice in high dimensions; why
> did you choose it for a 1 dim problem?
>
> Duncan Murdoch

It seems CG is having problems with the cube root.  This converges while
still using CG:

S1 <- optim(1001,function(x) (production1(x)^3), method = "CG",
  control = list(fnscale=-1))

On Thu, Mar 12, 2020 at 9:34 AM Skyler Saleebyan
<skylerbsaleebyan at gmail.com> wrote:
>
> I am trying to familiarize myself with optim() with a relatively simple
> maximization.
>
> Description:
> L and K are two terms which are constrained to add up to a total 100000
> (with respective weights to each). To map this constraint I plugged K into
> the function (to make this as simple as possible.)
>
> Together these two feed into one nonlinear function which is the product of
> two monotonic (on the positive interval) functions. Then that numbers is
> returned in a function fed to optim, which should maximize the output by
> adjusting L. The whole code is:
>
> production1 <- function(L){
>   budget=100000
>   Lcost=12
>   Kcost=15
>   K=(budget-L*Lcost)/Kcost
>   machines=0.05*L^(2/3)*K^(1/3)
>   return(machines)
> }
>
> # production1(6000) #example of number with much higher output vs optim
> result
> S1=optim(1001,production1,method="CG",control=list(fnscale=-1))
> S1
>
> Output:
> $par
> [1] 1006.536
>
> $value
> [1] 90.54671
>
> $counts
> function gradient
>      201      101
>
> $convergence
> [1] 1
>
> $message
> NULL
>
>
> For some reason this never explores the problem space and just spits out
> some answer close to the initial condition. What am I doing wrong?
>
> Thanks,
> Skyler S.
>
>         [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.


-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

######################################################
I ran <while (TRUE) 0> before posting, and waited a while...
(Re: The posting guide, which I'm going to start putting a lot more weight
on).

Noting, I was wondering if the posting guide has a mistake, because
<4*runif(1)> doesn't do anything special...
(Hopefully, Martin is reading this...).

And I'm planning to make this my last post in this thread...
######################################################
Here's a much simpler example of the problem:

optim (4, dnorm, method="CG", control = list (fnscale=-1) )$par

This problem isn't limited to objective functions of one variable.
I tried similar problems with functions of two, three and four variables.
And the same thing happened.

I'm not sure if this is a bug in the R code, or not.

If it's a bug...

And if it's not a bug, I'm struggling to see why anyone would want to
use this method in an applied setting...


On Fri, Mar 13, 2020 at 3:42 AM Duncan Murdoch <murdoch.duncan at
gmail.com> wrote:
>
> It looks like a bug in the CG method.  The other methods in optim() all
> work fine.  CG is documented to be a good choice in high dimensions; why
> did you choose it for a 1 dim problem?
>
> Duncan Murdoch