# Constrained optimization in R

My data is categorized by two different parameters (say F having n groups and S having m groups) and I want to get a relationship between the two. For example $F =${$f_1 , f_2 , f_3$} = {$10,10,5$} and $S =${$s_1, s_2 , s_3, s_4$} = {$8,8,8,4$}. (Read $f_1$ has 10 elements in it, $s_1$ has 8 elements in it.)

The problem is to get a relationship like $f_1 = p_{11}.s_1 + p_{12}.s_2 + p_{13}.s_3$ for all groups in F. The exact relationship is not possible (due to constraint 1 below), so we have to find the most approximate solution. The constraints are:

1. Sum of all the numbers in F (10+10+5) < sum of all the numbers in S (8+8+8+4). Note that, this is a default property of data and we don't have to apply this.

2. Total contribution of any $s_i$ can not be more than 100%. That is if $s_1$ contribute 40% in all 3 equations its total contribution will be 120% which should not happen.

3. Contribution can not be negative.

I framed this as an optimization problem as

Given:

$F = (f_1, f_2, ... f_n)$

$S = (s_1, s_2, ... s_m)$

Each $f_i$ can be represented as: $f_i = p_{i1} s_1 + p_{i2}s_2 + ... + p_{im}s_m + \epsilon_i$

That is: $P =${$p_{ij}$} is a $n * m$ matrix of row vectors $P_1, P_2, ..., P_n$ then $f_i = P_i.S+ \epsilon_i$

where . is the dot product (element by element multiplication).

I have to minimize $\sum_{i} \epsilon_i^2$, with following constraints:

1. $p_{11} + p_{12}+ ... + p_{1m}<=1$ That is some of each row of matrix P is less than 1.

2. $p_{ij} >=0$ for all i,j

I am trying to use constrOptim in R to solve this problem, however I am getting stuck with following:

1. Framing the constraints

2. How to create a code in which I just pass F and S and it gives me the matrix P and error vector as output.

Below is my code:

F = c(10,10,5)
S = c(8,8,8,4)
n = length(F)
m = length(S)
P_init = matrix(rep(0,m*n),nrow=n, ncol=m)

loss_fun <- function(P){

T = S*P #proportion matrix * S
F2 = rowSums(T) # Predicted values of F
E = F - F2  # Error
return(sum(E*E))
}
x = loss_fun(P_init)
z = constrOptim(P_init,loss_fun,NULL,ui=c(rep(-1,n)), ci=rep(-1,m))


Since the constraints are defined by: ui %*% theta - ci >= 0, I believe in my case ui = {-1,-1,-1}, theta = P_init ( a 3 x 4 matrix) and ci = {-1,-1,-1,-1}. However I get the following if I run the code. Error in ui %\*% theta : non-conformable arguments. Is theta = P_init is not true? Or there is some other error. As this is just one approach, I can also explore other approaches like optim or any other function.

I'm not sure if I'm right but help constrOptim {stats} say theta: numeric (vector) starting value (of length p): must be in the feasible region, and your theta is a matrix.

I think here is near to your idea (Input is a vector and the function makes it a matrix).
[EDITED] NOT NEAR because of difference constraints. Please refer to OP's ANSWER.

F = c(10,10,5)
S = c(8,8,8,4)
n = length(F)
m = length(S)
P_init2 = rep(0, 12)

loss_fun2 <- function(a){
P <- matrix(a, nrow=n, ncol=m)
T = S*P
F2 = rowSums(T)
E = F - F2
return(sum(E*E))
}

x = loss_fun2(rep(0, 12))
x  # 225
z = constrOptim(P_init2, loss_fun2, NULL, ui = rep(-1, n * m), ci = -1)

• Thanks for the answer, it clears the mistake I was doing, but it does not constraints the variables as I wanted. First constraint: $p_{i1} + p_{i2} + p_{i3} + p_{i4} < 0$ and 2nd: $p_{ij}>0$. z$par  0.3908330 0.7895228 0.4255702 -0.6770411 0.8194554 0.7260158  1.8119389 -3.9910594 -0.3568335 -0.6142514 1.6365515 -0.3395050. Certainly it violates the 2nd constraint. In this example, it does not violate first constraint as rowSums(matrix(z$par, nrow = 3, byrow=T)) gives  0.929 -0.634 0.326, so all rowSums are < 1, but I can't figure out if it will be true always – Gaurav Singhal Jun 15 '16 at 16:02
• @GauravSinghal How on earth is the sum of the $pi$s negative but each is positive? – AdamO Jun 15 '16 at 16:07
• @AdamO, sorry my bad, this was a typo in the comment, it should be $p_{i1}+p_{i2}+p_{i3}+p_{i4}<1$ as in the question. Thanks for pointing it out. – Gaurav Singhal Jun 15 '16 at 17:21
• @crayfish44, I used your advice to define P inside the function using P_init2, and to put constraints I used help from here stackoverflow.com/questions/16345271/…. Finally I was able to put the constraints as I wanted. – Gaurav Singhal Jun 15 '16 at 19:40
• @GauravSinghal I'm glad I could be of help. – cuttlefish44 Jun 15 '16 at 23:51

Using answer given by @crayfish and a detailed answer on how to put construct constraints here, I was able to come up with a solution.

#define F and S here
F = c(10,10,5)
S = c(8,8,9,8,4)

#loss_fun: to be minimized
loss_fun <- function(A){
P = matrix(A, nrow = n,ncol = m, byrow=T)
T = S*P #proportion matrix * S
F2 = rowSums(T) # Predicted values of F
E = F - F2  # Error
return(sum(E*E))
}

n = length(F)
m = length(S)
#Initial solution (theta)
P_init = c(rep(0.1,n*m))

# Creating Constraints

vi = matrix(rep(0,n*m*n*(m+1)),ncol = n*m, byrow = T)
for (i in 1:n){
for (j in 1:(n*m)) {
if (j <= m *  i & j > m * (i-1)) vi[i,j] = -1
}
}
for (i in (n+1):(n*m+n)) {
for (j in 1:(n*m)) {

if ((i-n) == j)  vi[i,j] = 1
}
}

myci = c(rep(-1,n),rep(0,n*m))

# check if initial value is in feasible region
vi %*% P_init - myci

#run the optimization module
z = constrOptim(P_init,loss_fun,NULL,ui=vi, ci=myci)

#result
P_final = matrix(z\$par,nrow=n,byrow=T)


Since creation of constraints may not be very easy to understand from the code, here is a visualization of constraints when n = 3 and m = 4

-1. p11 + -1. p12 + -1. p13 + -1. p14                                                                               >= -1
-1. p21 + -1. p22 + -1. p23 + -1. p24                                         >= -1
-1. p31 + -1. p32 + -1. p33 + -1. p34     >= -1
p11                                                                                                             >= 0
p12                                                                                                   >= 0
p13                                                                                         >= 0
p14                                                                               >= 0
p21                                                                       >= 0
p22                                                              >= 0
p23                                                   >= 0
p24                                          >= 0
p31                                   >= 0
p32                         >= 0
p33               >= 0
p34     >= 0