# How to use equalities as constrains with constrOptim in R

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?

• 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.
– whuber
Aug 16 '18 at 12:26
• 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 Aug 16 '18 at 19:31

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)

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  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).

• 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 Aug 16 '18 at 7:42
• @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.
– juod
Aug 16 '18 at 8:02
• 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 Aug 16 '18 at 9:24