0

I have to consider optimization problem in simulation study. An instance is given below:

library(mvtnorm)
library(alabama)

n = 200
q = 0.5
X <- matrix(0, nrow = n, ncol = 2)
X[,1:2] <- rmvnorm(n = n, mean = c(0,0), sigma = matrix(c(1,1,1,4),
                                                          ncol = 2))
x0 = matrix(c(X[1,1:2]), nrow = 1)
y0 = x0 - 0.5 * log(n) * (colMeans(X) - x0)
X = rbind(X, y0)

x01 = y0[1]
x02 = y0[2]
x1 = X[,1]
x2 = X[,2]

pInit = matrix(rep(1/(n + 1), n + 1), nrow = n + 1) 

f1 <- function(p) mean(((n + 1) * p ) ^ q)

heq1 <- function(p)
  c(sum(x1 * p) - x01, sum(x2 * p) - x02, sum(p) - 1)

sol <- alabama::auglag(pInit, fn = function(p) -f1(p), heq = heq1)
cat("The maximum objective value is:", -sol$value, '\n')

This gives error:

Error in eigen(a$hessian, symmetric = TRUE, only.values = TRUE) : 
  infinite or missing values in 'x'

I am not sure how to point out and overcome this problem. If this occurs due to mis-specification of initial point, how one can specify it in simulation work so that the program can itself sets suitable initial point and gives the right solution? Otherwise, why this error occurs and how to get rid of it? Can someone please help. Thanks!

Hans W.
  • 1,799
  • 9
  • 16
Jafer
  • 33
  • 5
  • 2
    Not reproducible and not runnable. See [how to make a reproducible example](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example). – Bhas Nov 03 '18 at 11:51
  • You should display `sol`. From `sol$convergence` and `sol$message` you will see that there is no convergence. The error message you get implies a hessian containing `Inf`, `Nan` or even `NA`. You'll have to closely check your system to determine if there is some sort of error somewhere. – Bhas Nov 03 '18 at 16:38
  • Forgot to mention this. You must add `,control.outer=list(kkt2.check=FALSE)` to the `auglag` function call. Then you can see the full value of `sol`. – Bhas Nov 03 '18 at 19:43
  • @Bhas Thank you for your tips. But how to fix the problem specifically in the given example is another problem. – Jafer Nov 04 '18 at 09:35

2 Answers2

1

This answer is an ADDENDUM to the first answer, especially targeting your second question about significantly speeding up the whole process.

To make the run time estimation reproducible, we will fix the seed; all other definitions are yours.

set.seed(4789)
n = 200
q = 0.5
X <- mvtnorm::rmvnorm(n = n, mean = c(0,0),
                      sigma = matrix(c(1,1,1,4), ncol = 2))
x0 = matrix(c(X[1,1:2]), nrow = 1)
y0 = x0 - 0.5 * log(n) * (colMeans(X) - x0)
X = rbind(X, y0)
x01 = y0[1]; x02 = y0[2]
x1 = X[,1]; x2 = X[,2]
pInit = matrix(rep(1/(n + 1), n + 1), nrow = n + 1)

First, let's do it with augmented Lagrangian and optim() as inner solver.

f1 <- function(p) sum(sqrt(pmax(0, p)))
heq1 <- function(p) c(sum(x1 * p) - x01, sum(x2 * p) - x02, sum(p) - 1)
hin1 <- function(p) p - 1e-06
system.time( sol <- alabama::auglag(pInit, fn = function(p) -f1(p), 
                           heq = heq1, hin = hin1) )
##    user  system elapsed 
##  24.631   0.054  12.324 
-1 * sol$value; heq1(sol$par)
## [1] 7.741285
## [1] 1.386921e-09 3.431108e-10 4.793488e-10

This problem is convex with linear constraints. Therefore we can apply an efficient convex solver such as ECOS. For modeling we will make use of the CVXR package.

# install.packages(c("ECOSolveR", "CVXR"))
library(CVXR)

p <- Variable(201)
obj <- Maximize(sum(sqrt(p)))
cons <- list(p >= 0, sum(p) == 1,
             sum(x1*p)==x01, sum(x2*p)==x02)
prbl <- Problem(obj, cons)
system.time( sol <- solve(prbl, solver="ECOS") )
##    user  system elapsed 
##   0.044   0.000   0.044 

ps <- sol$getValue(p)
cat("The maximum value is:", sum(sqrt(pmax(0, ps))))
## The maximum value is: 7.74226
c(sum(ps), sum(x1*ps) - x01, sum(x2*ps) - x02)
## [1]  1.000000e+00 -1.018896e-11  9.167819e-12

We see that the convex solver is about 500 times faster (!) than the first approach with a standard nonlinear solver. IMPORTANT: We do not need a starting value because a convex problem has only one optimum.

Hans W.
  • 1,799
  • 9
  • 16
0

As said before, see Maximizing nonlinear constraints problem using r package nloptr:
You have to prevent the solver to run into an area where your objective function is not defined, here this means p_i >= 0 for each index i. And if it does, let the objective function return some finite value. Simplifying your function (for q = 0.5) it looks, for instance, like

f1 <- function(p) sum(sqrt(pmax(0, p)))

Better also to provide an inequality constraint for p_i > 0 as

heq1 <- function(p) c(sum(x1 * p) - x01, sum(x2 * p) - x02, sum(p) - 1)
hin1 <- function(p) p - 1e-06

Now the solver returns a plausible result:

sol <- alabama::auglag(pInit, fn = function(p) -f1(p), 
                       heq = heq1, hin = hin1)

-1 * sol$value
## [1] 11.47805

and the equality conditions are all satisfied:

heq1(sol$par)
## [1] -4.969690e-09  5.906888e-09  1.808652e-08

All this can be done 'programmatically', with a bit of care, naturally.

Hans W.
  • 1,799
  • 9
  • 16
  • Thank you for this explanation. It works now but there is a problem. The procedure is too slow to apply in simulation study. For example, in my study unlike the case when q > 1 that takes around 3 minutes for one replication, this case (q = 0.5) takes 2.5 hours. How to speed it up? – Jafer Nov 05 '18 at 10:24
  • Yeah, this is what I supposed anyway. R with augmented Lagrangian and `optim()` as solver is too slow for doing simulation studies. There are maybe two or three ways out, but I think it will require a new question, comments are too small for that. You could also ask here or on the [Computational Science](https://scicomp.stackexchange.com/) Steckexchange Q&A site. And think about accepting my answer(s). – Hans W. Nov 05 '18 at 11:56
  • See in my ADDENDUM answer here on this page how solving this as a convex optimization problem will speed it up considerable. – Hans W. Nov 06 '18 at 12:51
  • Thank you Hans W. It works surprisingly much faster than previous one. – Jafer Nov 07 '18 at 07:16