Optimization of values in a matrix

Question

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

score 1 · Accepted Answer · answered Feb 23 '17 at 18:01

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

I thought that might be the only way but I wanted to make sure. It seems like it could slow the optimization down a lot. The matrices I am dealing with are pretty large and gosolnp has to call these functions many times just to do a single optimization. — user3390169, Feb 24 '17 at 01:05

Optimization of values in a matrix

1 Answers1

Linked