5
$\begingroup$

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.

Improve this question
$\endgroup$

2 Answers 2

Reset to default
3
$\begingroup$

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
Improve this answer
$\endgroup$
1
3
$\begingroup$

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)
Improve this answer
$\endgroup$
5
  • $\begingroup$ 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 [1] 0.3908330 0.7895228 0.4255702 -0.6770411 0.8194554 0.7260158 [7] 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 [1] 0.929 -0.634 0.326, so all rowSums are < 1, but I can't figure out if it will be true always $\endgroup$
    – Gaurav Singhal
    Commented Jun 15, 2016 at 16:02
  • 2
    $\begingroup$ @GauravSinghal How on earth is the sum of the $pi$s negative but each is positive? $\endgroup$
    – AdamO
    Commented Jun 15, 2016 at 16:07
  • $\begingroup$ @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. $\endgroup$
    – Gaurav Singhal
    Commented Jun 15, 2016 at 17:21
  • $\begingroup$ @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. $\endgroup$
    – Gaurav Singhal
    Commented Jun 15, 2016 at 19:40
  • $\begingroup$ @GauravSinghal I'm glad I could be of help. $\endgroup$
    – cuttlefish44
    Commented Jun 15, 2016 at 23:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.