# 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? 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 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 Nov 24 '16 at 16:31