Rsolnp (version 1.14)

gosolnp: Random Initialization and Multiple Restarts of the solnp solver.

Description

When the objective function is non-smooth or has many local minima, it is hard to judge the optimality of the solution, and this usually depends critically on the starting parameters. This function enables the generation of a set of randomly chosen parameters from which to initialize multiple restarts of the solver (see note for details).

Usage

gosolnp(pars = NULL, fixed = NULL, fun, eqfun = NULL, eqB = NULL, ineqfun = NULL, 
ineqLB = NULL, ineqUB = NULL, LB = NULL, UB = NULL, control = list(), 
distr = rep(1, length(LB)), distr.opt = list(), n.restarts = 1, n.sim = 20000, 
cluster = NULL, rseed = NULL, ...)

Arguments

pars

The starting parameter vector. This is not required unless the fixed option is also used.

fixed

The numeric index which indicates those parameters which should stay fixed instead of being randomly generated.

fun

The main function which takes as first argument the parameter vector and returns a single value.

eqfun

(Optional) The equality constraint function returning the vector of evaluated equality constraints.

eqB

(Optional) The equality constraints.

ineqfun

(Optional) The inequality constraint function returning the vector of evaluated inequality constraints.

ineqLB

(Optional) The lower bound of the inequality constraints.

ineqUB

(Optional) The upper bound of the inequality constraints.

LB

The lower bound on the parameters. This is not optional in this function.

UB

The upper bound on the parameters. This is not optional in this function.

control

(Optional) The control list of optimization parameters. The eval.type option in this control list denotes whether to evaluate the function as is and exclude inequality violations in the final ranking (default, value = 1), else whether to

distr

A numeric vector of length equal to the number of parameters, indicating the choice of distribution to use for the random parameter generation. Choices are uniform (1), truncated normal (2), and normal (3).

distr.opt

If any choice in distr was anything other than uniform (1), this is a list equal to the length of the parameters with sub-components for the mean and sd, which are required in the truncated normal and normal distributions.

n.restarts

The number of solver restarts required.

n.sim

The number of random parameters to generate for every restart of the solver. Note that there will always be significant rejections if inequality bounds are present. Also, this choice should also be motivated by the width of the upper and lower bound

cluster

If you want to make use of parallel functionality, initialize and pass a cluster object from the parallel package (see details), and remember to terminate it!

rseed

(Optional) A seed to initiate the random number generator, else system time will be used.

...

(Optional) Additional parameters passed to the main, equality or inequality functions

Value

• A list containing the following values:
• pars
  Optimal Parameters.
• convergence
  Indicates whether the solver has converged (0) or not (1).
• values
  Vector of function values during optimization with last one the value at the optimal.
• lagrange
  The vector of Lagrange multipliers.
• hessian
  The Hessian at the optimal solution.
• ineqx0
  The estimated optimal inequality vector of slack variables used for transforming the inequality into an equality constraint.
• nfuneval
  The number of function evaluations.
• elapsed
  Time taken to compute solution.
• start.pars
  The parameter vector used to start the solver

Details

Given a set of lower and upper bounds, the function generates, for those parameters not set as fixed, random values from one of the 3 chosen distributions. Depending on the eval.type option of the control argument, the function is either directly evaluated for those points not violating any inequality constraints, or indirectly via a penalty barrier function jointly comprising the objective and constraints. The resulting values are then sorted, and the best N (N = random.restart) parameter vectors (corresponding to the best N objective function values) chosen in order to initialize the solver. Since version 1.14, it is up to the user to prepare and pass a cluster object from the parallel package for use with gosolnp, after which the parLapply function is used. If your function makes use of additional packages, or functions, then make sure to export them via the clusterExport function of the parallel package. Additional arguments passed to the solver via the ...option are evaluated and exported by gosolnp to the cluster.

References

Y.Ye, Interior algorithms for linear, quadratic, and linearly constrained non linear programming, PhD Thesis, Department of EES Stanford University, Stanford CA. Hu, X. and Shonkwiler, R. and Spruill, M.C. Random Restarts in Global Optimization, 1994, Georgia Institute of technology, Atlanta.

Examples

# [Example 1]
# Distributions of Electrons on a Sphere Problem:
# Given n electrons, find the equilibrium state distribution (of minimal Coulomb 
# potential) of the electrons positioned on a conducting sphere. This model is 
# from the COPS benchmarking suite. See http://www-unix.mcs.anl.gov/~more/cops/.	
gofn = function(dat, n)
{
	
	x = dat[1:n]
	y = dat[(n+1):(2*n)]
	z = dat[(2*n+1):(3*n)]
	ii = matrix(1:n, ncol = n, nrow = n, byrow = TRUE)
	jj = matrix(1:n, ncol = n, nrow = n)
	ij = which(ii<jj, arr.ind = TRUE)
	i = ij[,1]
	j = ij[,2]
	#  Coulomb potential
	potential = sum(1.0/sqrt((x[i]-x[j])^2 + (y[i]-y[j])^2 + (z[i]-z[j])^2))
	potential
}

goeqfn = function(dat, n)
{
	x = dat[1:n]
	y = dat[(n+1):(2*n)]
	z = dat[(2*n+1):(3*n)]
	apply(cbind(x^2, y^2, z^2), 1, "sum")
}

n = 25
LB = rep(-1, 3*n)
UB = rep(1,  3*n)
eqB = rep(1, n)
ans = gosolnp(pars  = NULL, fixed = NULL, fun = gofn, eqfun = goeqfn, eqB = eqB, 
LB = LB, UB = UB, control = list(outer.iter = 100, trace = 1), 
distr = rep(1, length(LB)), distr.opt = list(), n.restarts = 2, n.sim = 20000, 
rseed = 443, n = 25)
# should get a function value around 243.813

# [Example 2]
# Parallel functionality for solving the Upper to Lower CVaR problem (not properly
# formulated...for illustration purposes only).

mu =c(1.607464e-04, 1.686867e-04, 3.057877e-04, 1.149289e-04, 7.956294e-05)
sigma = c(0.02307198,0.02307127,0.01953382,0.02414608,0.02736053)
R = matrix(c(1, 0.408, 0.356, 0.347, 0.378,  0.408, 1, 0.385, 0.565, 0.578, 0.356, 
0.385, 1, 0.315, 0.332, 0.347, 0.565, 0.315, 1, 0.662, 0.378, 0.578, 
0.332, 0.662, 1), 5,5, byrow=TRUE)
# Generate Random deviates from the multivariate Student distribution
set.seed(1101)
v = sqrt(rchisq(10000, 5)/5)
S = chol(R)
S = matrix(rnorm(10000 * 5), 10000) %*% S
ret = S/v
RT = as.matrix(t(apply(ret, 1, FUN = function(x) x*sigma+mu)))
# setup the functions
.VaR = function(x, alpha = 0.05)
{
	VaR = quantile(x, probs = alpha, type = 1)
	VaR
}

.CVaR = function(x, alpha = 0.05)
{
	VaR = .VaR(x, alpha)
	X = as.vector(x[, 1])
	CVaR = VaR - 0.5 * mean(((VaR-X) + abs(VaR-X))) / alpha
	CVaR
}
.fn1 = function(x,ret)
{
	port=ret%*%x
	obj=-.CVaR(-port)/.CVaR(port)
	return(obj)
}

# abs(sum) of weights ==1
.eqn1  = function(x,ret)
{
	sum(abs(x))
}

LB=rep(0,5)
UB=rep(1,5)
pars=rep(1/5,5)
ctrl = list(delta = 1e-10, tol = 1e-8, trace = 0)
cl = makePSOCKcluster(2)
# export the auxilliary functions which are used and cannot be seen by gosolnp
clusterExport(cl, c(".CVaR", ".VaR"))
ans = gosolnp(pars, fun = .fn1, eqfun = .eqn1, eqB = 1, LB = LB, UB = UB, 
control = list(n.restarts = 2, n.sim=500), cluster = cl, ret = RT)
ans
# don't forget to stop the cluster!
stopCluster(cl)