# Constrained Optimization in R - again

Note: I have asked a similar question earlier here, I got some good advice which helped me to create a solution using constrOptimfunction. It worked but it was very slow. So I asked seeked help to optimize it here , and someone suggested to use auglag function from Alabama package. The solution is much faster, however I am asking it as a new question as the equality constraint (# 3) was not there and I am unable implement this additional constraint using auglag

Question: I have an optimization problem which can be explained as:

$F =$ {$f_1, f_2, ... f_n$} ; where $f_i$ are scalars, is a column vector of size $n$. $n$ is generally between 35 and 55

$S =$ {$s_1, s_2, ... s_m$} ; where $s_j$ are scalars, is a row vector of size $m$. $m$ is generally between 100 and 200

$C =$ {$c_{ij}$} ; is a $n \times m$ matrix

Here $F$ and $S$ are fixed, and $C$ is variable.

The constraints on $C$ are:

1. $0 <= c_{ij} <= 1$
2. $\sum_{i=1}^n c_{ij} <= 1$
3. $c_{ij}=0$ for some particular $i,$ and $j$ based on research data

I am predicting $f_i$ as:

$\hat{f_i} = \sum_{j=1}^m c_{ij}s_j$

I want to minimize $E = \sum_{i=1}^n (f_i - \hat{f_i})^2$

In matrix notation $E = (F-CS^T)^T.(F-CS^T)$

The question is which package/function should I use to effectively solve this problem in R. By effectively I mean, run time < 30 minutes, and an approximate global optimal solution. Till now I have tried constrOptim, DEoptim, and auglag (package alabama), but I either couldn't solve the problem effectively or couldn't solve it at all. For auglag here is a reproducible version of my code:

library(alabama)

F = c(10,10,5)
S = c(8,8,9,8,4)

n = length(F)
m = length(S)

# Creating a dummy matrix to determine the coefficients under constraint #3
dat_mat = c(rep(0, n*m))
x = sample.int(n*m,round(n*m/2)) # In real data, around 80% to 90% of the values will be NA, but each column will have atleast one non-NA value
dat_mat[x] <-NA
dat_mat = matrix(dat_mat,nrow=n, ncol = m, byrow = T) #whenever any element of dat_mat is NA, the corresponding coefficient is 0.

#Initial solution
P_init = c(rep(0, n*m))
for (i in 1:n) {
for (j in 1:m) {
if (!is.na(dat_mat[i,j])) P_init[m*(i-1)+j] = 0.0001 #choice of .0001 is random, is there a better choice
}
}

loss_fun <- function(P){
T_mat <- t(matrix(P, nrow = m,byrow = FALSE)*S) #Created a matrix with m rows to avoid double transpose and save time
F2 = rowSums(T_mat) # Predicted values of F
E = F - F2  # Vector containing error for each of the n functions
return(sum(E*E))
}

# Creating inequality constraint function - works fine
hin <- function(P){
P_mat <- matrix(P, nrow = n,byrow = TRUE)
c(colSums(P_mat) * -1 +1, P)
}

# Creating equality constraint function - Including this Throws an error
heq <- function(P) {
P[which(P==0)]
}

# Creating gradient function - Including this gives suboptimal solution. Wanted to use gradient as document says it can decrease the run-time.
gr <- function(P){
P_mat = matrix(P,nrow=n,byrow=TRUE)
gr = ((F - P_mat %*% S)) %*% t(S)
gr = gr / max(gr)
}

system.time(aug_res <- auglag(P_init, loss_fun, hin = hin, control.outer = list(kkt2.check = FALSE))) # works as expected
system.time(aug_res <- auglag(P_init, loss_fun, hin = hin, gr = gr, control.outer = list(kkt2.check = FALSE))) #gives suboptimal solution
system.time(aug_res <- auglag(P_init, loss_fun, hin = hin, heq = heq, control.outer = list(kkt2.check = FALSE))) #throws an error

Error in colSums(lam[-inactive] * ij) :
'x' must be an array of at least two dimensions
In addition: There were 13 warnings (use warnings() to see them)
Timing stopped at: 0.01 0 0.01

> warnings()
Warning messages:
1: In func(dx, ...) - f :
longer object length is not a multiple of shorter object length
9: In lam * d0 : longer object length is not a multiple of shorter object length.


To validate/see the solution, I use this:

#Validation
P_final = matrix(aug_res$par,nrow=n,byrow=TRUE) T = t(S*t(P_final)) #proportion matrix * S, transpose to ensure multiplication is by row. (F2 = rowSums(T)) # Predicted values of F F # Actual values of F  In real problem,$F$,$S\$ and dat_mat comes from the real data.

• It is generally hard to solve optimization with equality constraints. Is there any way you can relax the "=0" constraint? – forecaster Nov 24 '16 at 13:36
• Ya, I can relax it like "<.0001". Your comment is food for thought. How to apply it without creating a huge matrix is I have to think now – Gaurav Singhal Nov 24 '16 at 16:13
• Maybe we can apply some maximum limit like .01 on colSums or rowSums of all these coefficients. Or on total sum of all these coefficients as well – Gaurav Singhal Nov 24 '16 at 16:31