R: Isotonic regression Minimisation

Question

I want minimize the following equation:

F=SUM{u 1:20}sum{w 1:10}   Quw(ruw-yuw)

with the following constraints:

yuw >= yu,w+1
yuw >= yu-1,w
y20,0 >= 100
y0,10 >= 0

I have a 20*10 ruw and 20*10 quw matrix, I now need to generate a yuw matrix which adheres to the constraints. I am coding in R and am familiar with the lpsolve and optimx packages, but don't know how to use them for this particular question.

Answer 1

Because Quw and ruw are both data, all constraints as well as the objective are linear in the yuw decision variables. As a result, this is a linear programming problem that can be approached by the lpSolve package.

To abstract this out a bit, let's assume R=20 and C=10 describe the dimensions of the input matrices. Then there are R*C decision variables and we can assign them order y11, y21, ... yR1, y12, y22, ... yR2, ..., y1C, y2C, ..., yRC, reading down the columns of the matrix of variables.

The coefficient of each yuw variable in the objective is -Quw; note that the Quw*ruw terms in the summation are constants (aka not affected by the values we select for the decision variables) and therefore not inputted to the linear programming solver. Interestingly, this means that ruw actually has no effect on the optimization model solution.

The first R*(C-1) constraints correspond to the yuw >= yu,w+1 constraints, and the next (R-1)*C constraints correspond to the yuw >= yu-1,w constraints. The final two constraints correspond to the y20,1 >= 100 and y1,10 >= 0 constraints.

We can input this model into the lpsolve package with the following R command, taking as input a simple Q matrix where each entry is -1 (the resulting solution should have all decision variables set to 0 except the bottom left corner, which should be 100):

# Sample data
Quw <- matrix(-1, nrow=20, ncol=10)
R <- nrow(Quw)
C <- ncol(Quw)

# Build constraint matrix
part1 <- matrix(0, nrow=R*(C-1), ncol=R*C)
part1[cbind(1:(R*C-R), 1:(R*C-R))] <- 1
part1[cbind(1:(R*C-R), (R+1):(R*C))] <- -1
part2 <- matrix(0, nrow=(R-1)*C, ncol=R*C)
pos2 <- as.vector(sapply(2:R, function(r) r+seq(0, R*(C-1), by=R)))
part2[cbind(1:nrow(part2), pos2)] <- 1
part2[cbind(1:nrow(part2), pos2-1)] <- -1
part3 <- rep(0, R*C)
part3[R] <- 1
part4 <- rep(0, R*C)
part4[(C-1)*R + 1] <- 1
const.mat <- rbind(part1, part2, part3, part4)

library(lpSolve)
mod <- lp(direction = "min",
          objective.in = as.vector(-Quw),
          const.mat = const.mat,
          const.dir = rep(">=", nrow(const.mat)),
          const.rhs = c(rep(0, nrow(const.mat)-2), 100, 0))

We can now access the model solution:

mod
# Success: the objective function is 100
matrix(mod$solution, nrow=R)
#       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#  [1,]    0    0    0    0    0    0    0    0    0     0
#  [2,]    0    0    0    0    0    0    0    0    0     0
#  [3,]    0    0    0    0    0    0    0    0    0     0
#  [4,]    0    0    0    0    0    0    0    0    0     0
#  [5,]    0    0    0    0    0    0    0    0    0     0
#  [6,]    0    0    0    0    0    0    0    0    0     0
#  [7,]    0    0    0    0    0    0    0    0    0     0
#  [8,]    0    0    0    0    0    0    0    0    0     0
#  [9,]    0    0    0    0    0    0    0    0    0     0
# [10,]    0    0    0    0    0    0    0    0    0     0
# [11,]    0    0    0    0    0    0    0    0    0     0
# [12,]    0    0    0    0    0    0    0    0    0     0
# [13,]    0    0    0    0    0    0    0    0    0     0
# [14,]    0    0    0    0    0    0    0    0    0     0
# [15,]    0    0    0    0    0    0    0    0    0     0
# [16,]    0    0    0    0    0    0    0    0    0     0
# [17,]    0    0    0    0    0    0    0    0    0     0
# [18,]    0    0    0    0    0    0    0    0    0     0
# [19,]    0    0    0    0    0    0    0    0    0     0
# [20,]  100    0    0    0    0    0    0    0    0     0

Note that your model could easily become infeasible if Quw changed (for instance if we filled it with 1 instead of -1). In these cases the model will exit with status 3 (you can see this by running mod or mod$status).