Ошибка в задаче нелинейной оптимизации: бесконечные или отсутствующие значения в "x"

Question

Я должен рассмотреть проблему оптимизации в симуляционном исследовании. Пример приведен ниже:

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')

Это дает ошибку:

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

Я не уверен, как указать и преодолеть эту проблему. Если это происходит из-за неправильной спецификации начальной точки, как можно указать ее в симуляционной работе, чтобы программа могла сама установить подходящую начальную точку и дать правильное решение? Иначе, почему возникает эта ошибка и как от нее избавиться? Может кто-нибудь, пожалуйста, помогите. Спасибо!

Answer 1

Как уже было сказано выше, смотрите задачу Максимизация нелинейных ограничений с использованием пакета r nloptr:
Вы должны предотвратить попадание решателя в область, где ваша целевая функция не определена, здесь это означает p_i >= 0 для каждого индекса i, И если это так, пусть целевая функция возвращает некоторое конечное значение. Упрощение вашей функции (для q = 0.5) это выглядит, например, как

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

Лучше также обеспечить ограничение неравенства для p_i > 0 как

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

Теперь решатель возвращает правдоподобный результат:

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

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

и условия равенства все выполнены:

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

Все это может быть сделано "программно", с небольшой осторожностью, естественно.

Answer 2

Этот ответ является ДОПОЛНЕНИЕМ к первому ответу, особенно нацеленный на ваш второй вопрос о значительном ускорении всего процесса.

Чтобы сделать оценку времени выполнения воспроизводимой, мы исправим начальное значение; все остальные определения ваши.

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) 

Во-первых, давайте сделаем это с расширенным лагранжианом и optim() как внутренний решатель.

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

Эта проблема является выпуклой с линейными ограничениями. Поэтому мы можем применить эффективный выпуклый решатель, такой как ECOS. Для моделирования мы будем использовать пакет CVXR.

# 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

Мы видим, что выпуклый решатель примерно в 500 раз быстрее (!), Чем первый подход со стандартным нелинейным решателем. ВАЖНО: Нам не нужно начальное значение, потому что у выпуклой задачи есть только один оптимум.