Как ускорить процесс нелинейной оптимизации в R

Рассмотрим следующий пример задачи нелинейной оптимизации. Процедура слишком медленная, чтобы ее можно было применять в симуляционных исследованиях. Например, в случае с моими исследованиями всего на одну репликацию уходит 2,5 часа. Как ускорить процесс, чтобы можно было оптимизировать время обработки?

library(mvtnorm)
library(alabama)

n = 200
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(0.1, n + 1), nrow = n + 1)

outopt = list(kkt2.check=FALSE, "trace" = FALSE)
f1 <- function(p) sum(sqrt(pmax(0, p)))/sqrt(n+1)
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,
                       control.outer = outopt)
-1 * sol$value