[R] Can R solve this optimization problem?

[Paul Smith]

On Jan 8, 2008 2:58 AM, Gabor Grothendieck <ggrothendieck at gmail.com> wrote:
> > Thanks again, Gabor. Apparently, something is missing in your
> > approach, as I tried to apply it to the original problem (without sin)
> > and the result seems not being correct::
> >
> > f <- function(z) {
> >  a <- sum(-(999:1)*z)
> >  return(a)
> > }
> >
> > result <- optim(rep(0,999), f,
> > NULL,lower=rep(-1/1000,999),upper=rep(1/1000,999),method="L-BFGS-B")
> > b <- result$par
> >
> > x <- 0
> >
> > for (i in 1:999) {
> >  x[i+1] <- b[i] + x[i]
> > }
> >
> > The vector x contains the solution, but it does not correspond to the
> > analytical solution. The analytical solution is
> >
> > x(t) = t, if t <= 1/2 and
> >
> > x(t) = 1 - t, if t >= 1/2.
> >
> > Am I missing something?
> >
>
> No, but I was -- the x[n] = 0 constraint.  Add it through
> a penalty term and try this:
>
> theta <- 1000
> n <- 100
> x <-  rep(1/n, n)
> f <- function(x) sum(n:1 * x) - theta * sum(x)^2
> optim(x, f, lower = rep(-1, n), upper = rep(1, n), method = "L-BFGS-B",
>   control = list(maxit = 1000, fnscale = -1, lmm = 25))
> plot(cumsum(res$par))

Thanks, Gabor. Inspired by your code, I have written the following:

theta <- 50000000
n <- 1000
x <-  rep(1/n, n)
f <- function(x) sum(n:1 * x) - theta * sum(x)^2
res <- optim(x, f, lower = rep(-1/n, n), upper = rep(1/n, n), method =
"L-BFGS-B",control = list(maxit = 2000, fnscale = -1))

b <- res$par

x <- 0

for (i in 1:n) {
  x[i+1] <- b[i] + x[i]
}

plot(seq(0,1,by=1/n),x,type="l")

to solve the same problem. It works fine unless one replaces

theta <- 50000000

by

theta <- 1e100

According to a book that I consulted meanwhile, the larger the value
of theta the better the solution should be. Why does not it happen
with my example? Any clues?

Paul