CODEWITHSUNDEEP
roptimization
user3390169
user3390169

Reputation: 1045

Optimization of values in a matrix

I typically use Rsolnp for optimization but I am having trouble figuring out how to ask R to find values to fill a matrix (instead of a vector). Is that possible with Rsolnp or any other optimizer?

Here is a simplified example that is not working:

library(Rsolnp)

a<-matrix(rnorm(9), ncol=3)
b<-matrix(rnorm(9), ncol=3)

f1<-function(H) {
    return(sum(H*a))
}

f2<-function(H) {
    return(sum(H*b))
}

lH<-matrix(rep(0, 9), ncol=3)
uH<-matrix(rep(1, 9), ncol=3)
pars<-uH
target<-1.2

sol <- gosolnp(pars, fixed=NULL, fun=f1, eqfun=f2, eqB=target, LB=lH, UB=uH, distr=uH, n.restarts=10, n.sim=20000, cluster= NULL)

As you can see from the output, Rsolnp seems to be confused by the request:

> sol
$values
[1] 1e+10

$convergence
[1] 0

$pars
[1] NA NA NA NA NA NA NA NA NA

$start.pars
[1] 0.90042133 0.33262541 0.94586530 0.02083822 0.99953060 0.10720068 0.14302770 0.67162637 0.25463806

$rseed
[1] 1487866229

Upvotes: 0

Views: 1022

Answers (1)

Stibu
Stibu

Reputation: 15897

It seems that gosolnp() does not work with matrices. I went through the function in debugging mode and there is a call of solnp() that fails with the message:

Error in pb/cbind(vscale[(neq + 2):(neq + mm + 1)], vscale[(neq + 2):(neq + : non-conformable arrays

But since a matrix is just a vector with the dimension attribute set, you can always reformulate your problem in terms of vectors. In your case, this is very easy, because you never do something that actually requires a matrix (like, for instance, a matrix product). Just omitting matrix() everywhere works fine.

But I assume that this is just a property of your simplified problem and your actual problem indeed needs to be expressed in terms of matrices. You could get around the problem by converting your vectors into matrices only inside the functions f1() and f2() as follows:

f1 <- function(H) {
    return(sum(matrix(H, ncol = 3) * a))
}

f2 <- function(H) {
    return(sum(matrix(H, ncol = 3) * b))
}

You can then define a and b as matrices as before, but lH and uH must be vectors:

a <- matrix(rnorm(9), ncol=3)
b <- matrix(rnorm(9), ncol=3)

lH <- rep(0, 9)
uH <- rep(1, 9)
pars <- uH
target <- 1.2

And now you can call gosolnp():

sol <- gosolnp(pars, fixed = NULL, fun = f1, eqfun = f2,
               eqB = target, LB = lH, UB = uH, distr = uH,
               n.restarts = 10, n.sim = 20000, cluster = NULL)
sol$pars
## [1] 3.917819e-08 9.999997e-01 4.748336e-07 1.000000e+00 5.255060e-09 5.114680e-10
## [7] 4.899963e-01 1.000000e+00 9.260947e-08

Upvotes: 1

Related Questions