0
$\begingroup$

I want to solve a matrix system which have several solutions (infinite since it is overdetermined - 6 equations with 8 unknowns). However, the way I want to do it is to a criterion for the variables, which should be minimized, to find a specific solution. I am working in R

My first through was to use constrOptim in R However, this requires for me to set constraints that are equalities, and not "larger than or equal to", which is the the default in constrOptim.

The problem is essentially that I have a matrix, where I know none of the elements, and I want the the rowSums to have specific values, and the colSums to have specific values. This however, should be done while minimizing criterias that the i,j element divided by the sum of the i'th row, should be as close as possible to the sum of the j'th column divided by the sum of the whole matrix.

Does anyone know how I could solve this problem using constrOptim, or by using another function?

$\endgroup$
  • $\begingroup$ You have too few criteria to make this interesting, because--provided none of the denominators is zero--there always exists an exact solution. Letting the matrix elements be $x_{ij}$ and writing the row sums as $x_{i\cdot},$ the column sums as $x_{\cdot j},$ and the sum as $x_{\cdot\cdot},$ the solution is $$\hat x_{ij} = \frac{x_{i\cdot} x_{\cdot j} }{x_{\cdot\cdot}}.$$ Since this requires no optimization and is fast and easy to compute, it suggests you might not have formulated the question you meant to ask. $\endgroup$ – whuber Aug 16 '18 at 12:26
  • 1
    $\begingroup$ Thanks for the reponse. No, I formulated what I wanted to ask. I, however, might not have been precise. This can obviously be solved much quicker. However, I want the algorithm, to go the optimization route, because I will later on get specific restriction for certain of the $x_{i,j}$'s, which have to hold as well $\endgroup$ – pkpkPPkafa Aug 16 '18 at 19:31
1
$\begingroup$

So I found a solution that works. This might be relevant for others who find the question in the future.

Here cat1 and cat2 are the sums which are requried from the matrices.

The code that works:

require(Rsolnp)


cat1 <- c(10000,3000,2000)
cat2 <- c(5000, 2500,7500)

sum(cat11) == sum(cat2)
elements <- length(cat2)*length(cat1)


cat1.share <- c(cat1/sum(cat1))
cat1.counter <- rep(1:length(cat1), 3)

#specify your function
min.func <- function(x) {
sum.holder <- cat1[cat1.counter[1]]/sum(cat1)*cat2[ceiling(1/length(cat1)]-x[1]
    for (i in 2:elements)  {
    sum.holder <- c(sum.holder,cat1[cat1.counter[i]]/sum(cat1)*cat2[ceiling(i/length(cat1))]-x[i])
}
  sum(abs(sum.holder))
}

#specify the equality function. The number 15 (to which the function is equal)
#is specified as an additional argument
equal <- function(x) {

  y <- x[1:elements]

  matrix.hold <- matrix(y, nrow = length(cat2), ncol = length(cat1), byrow = TRUE) 

  hold <- c(rowSums(matrix.hold),
            colSums(matrix.hold)
  ) 
  hold <- hold[-length(hold)]
  hold

}

x.start<- rep(min(cat2)/9, elements)
x.low <- rep(0, elements)
x.high <- rep(max(c(cat1,cat2)), elements)


#the optimiser
opt <- solnp(x.start, #starting values 
             min.func, #function to optimise
             eqfun=equal, #equality function 
             eqB=c(cat2[1:length(cat2)], cat1[1:(length(cat1)-1)]),       #equality constraint
             LB=x.low, #lower bound for parameters i.e. greater than zero
             UB=x.high,
             control = c(delta = 1.0e-8, tol = 1e-11)) #upper bound for parameters 

opt.matrix <- matrix(opt$pars, nrow = length(cat2), ncol = length(cat1), byrow=TRUE)

    opt.matrix
$\endgroup$
0
$\begingroup$

Boring way: use alabama package which allows equality constraints.

Fun way: refactor your problem to eliminate equality constraints.
For example, if you were searching for a $2 \times 2$ matrix $$ \begin{bmatrix} a & b \\ c & d\end{bmatrix} $$ with the constraints $r_1, r_2, c_1, c_2$ such that $$ a + b = r_1 \\ c + d = r_2 \\ a + c = c_1 \\ b + d = c_2$$

you don't need to optimize over four variables. Optimize over $a$ and then define the other variables inside the cost function as: $$ b = r_1 - a \\ c = c_1 - a \\ d = c_2 - b$$

(see this answer for more formal explanation.)

If $a\in \mathbf{R}$, you don't need any constrains at all and can use optim straight out of the box. In terms of programming, it would looks something like:

costfn <- function(a){
    # obtain other elements using predefined rowSum and colSum constants
    b = r1 - a
    c = c1 - a
    d = c2 - b
    # calculate actual cost
    sum(a/r1 - c1/(c1+c2), b/r1 - c2/(c1+c2), c/r2 - c1/(c1+c2), d/r2 - c1/(c1+c2))
}
optim(costfn, ...)

Essentially, the cost function would take only $a$ as argument and define everything else inside the function. Similarly, one can trick optim to take bounded search spaces:

cost <- function(a){ # a in reals
    b = exp(a) + 1   # now b>1
    # use b as the parameter to calculate cost...
}
optim(cost, ...)

Just note that it can affect the search parameters if your problem is sensitive to that (e.g. a uniform grid is no longer uniform after exp).

$\endgroup$
  • $\begingroup$ Thanks for the rapid response. However, I am unsure how that is eliminating the constraints. As I see it there, it is only one of the equality constraints that is eliminated $\endgroup$ – pkpkPPkafa Aug 16 '18 at 7:42
  • $\begingroup$ @pkpkPPkafa Edited with some more details. Not sure if I represented your cost function correctly, but the trick is to calculate all other variables inside the cost function, using the the given constants. That way you enforce equality without requiring a separate constrain. $\endgroup$ – juod Aug 16 '18 at 8:02
  • $\begingroup$ thanks so much. I found a solution using the solnp, which is pretty flexible. It works for if I have a 3x3 matrix. When I extend it beyond that I get a error message. See separate answer I have added so others might use it $\endgroup$ – pkpkPPkafa Aug 16 '18 at 9:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.