matrix with row and column constraints

Question

I need to solve for a n x n (n usually <12) matrix subject to a few constraints:

1.Predetermined row and column sums are satisfied.

2.Each element in the matrix having a row number greater than column number must be zero (so basically the only nonzero elements must be in the top right portion).

3.For a given row, every element more than three columns to the right first nonzero element must also be zero.

So, a 4x4 matrix might look something like this (the row and column constraints will be much larger in practice, usually around 1-3 million):

|3 2 1 0| = 6
|0 2 1 1| = 4
|0 0 2 1| = 3
|0 0 0 4| = 4
 3 4 4 6

I have been trying to use some solver approaches to do this in excel and also have tried some R based optimization packages but have been so unsuccessful so far.

Any suggestions on how else I might approach this would be much appreciated.

Thanks!

Ben Bolker Ben Bolker · Accepted Answer · 2013-07-25T15:37:37

Test data:

x <- c(2,2,2,1,1,1,1)
rowVals <- c(6,4,3,4)
colVals <- c(3,4,4,6)

Function to construct the appropriate test matrix from (3N-5) parameters:

makeMat <- function(x,n) {
    ## first and last element of diag are constrained by row/col sums
    diagVals <- c(colVals[1],x[1:(n-2)],rowVals[n])
    ## set up off-diagonals 2,3
    sup2Vals <- x[(n-1):(2*n-3)]
    sup3Vals <- x[(2*n-2):(3*n-5)]
    ## set up matrix
    m <- diag(diagVals)
    m[row(m)==col(m)-1] <- sup2Vals
    m[row(m)==col(m)-2] <- sup3Vals
    m
}

Objective function (sum of squares of row & column deviations):

objFun <- function(x,n) {
    m <- makeMat(x,n)
    ## compute SSQ deviation from row/col constraints
    sum((rowVals-rowSums(m))^2+(colVals-colSums(m))^2)
}

Optimizing:

opt1 <- optim(fn=objFun,par=x,n=4)
## recovers original values, although it takes a lot of steps

opt2 <- optim(fn=objFun,par=rep(0,length(x)),n=4)
makeMat(opt2$par,n=4)
##      [,1]     [,2]      [,3]      [,4]
## [1,]    3 2.658991 0.3410682 0.0000000
## [2,]    0 1.341934 1.1546649 1.5038747
## [3,]    0 0.000000 2.5042858 0.4963472
## [4,]    0 0.000000 0.0000000 4.0000000
## 

## conjugate gradients might be better
opt3 <- optim(fn=objFun,par=rep(0,length(x)),n=4,
              method="CG")

It seems that there are multiple solutions to this problem, which isn't surprising (since there are 2N constraints on (N-2)+(N-1)+(N-2)= 3N-5 parameters).

You didn't say whether you need integer solutions or not -- if so you will need more specialized tools ...

matrix with row and column constraints

1 Answers