Handling sparse matrix?

[jshong]

Hi.

I wonder whether I can express the following constraints more efficiently. 

N = 5, M = 1000

z = binvar(N,M);

for i = 1:M

   

   for k = 1:size(idx{i},2)

                            

           const_z = [const_z, z(:,i) >= z(:,idx{i}(k))];

                

  end

  

end

  

Basically, what I want to express is the inequality constraints of two binary variables in some predetermined pair. 

For example, z_1_1 >= z_1_10, z_1_1 >= z_1_15 , z_1_2 > z_1_17 and so on. Here, 10, 15 and 17 are predetermined. 

When I run with above loop expression, it never starts to solve the problem. 

I think this is due to sparsity of constraints.

Is there any way I can make above expression efficient?

Interestingly, when M = 100, it starts to solve in less than 2-3 minutes. 

Thanks.

-JS

[Johan Löfberg]

Is the code really correct. It doesn't match your description.

When i = k = 1, you generate 5 constraints

[z(1,1);z(2,1)... ;z(5,1)] >=[z(1,r);z(2,r)... ;z(5,r)]

where r is idx{i}(k)). Is that what you meant? In other words, your matrix is a bunch of column vectors [z1 z2...z1000] and you are generating a large amount of constraint of the vectors z1>=z?, z1>=z?, z2>=z?, etc etc.

[jshong]

You are right. I meant 

z(1,i);z(2,i)... ;z(5,i)] >=[z(1,idx(i));z(2,idx(i))... ;z(5,idx(i))]

for some predetermined index of i,i.e. idx(i). 

And I want to impose these constraints all over i =1:1000.

Thanks. 

[Johan Löfberg]

I don't understand. That is exactly what I just wrote, no difference.

To summarize: Your model wants column i to be larger than column idx(i)?

[jshong]

Yes. 

[Johan Löfberg]

So you want to do this more efficiently (here saying row i larger than row i+1. How long is your list idx?

N = 5, M = 1000

z = binvar(N,M);

idx = [2:M];

const_z = [];

for i = 1:M   

   const_z = [const_z, z(:,i) >= z(:,idx(i))];

end

[jshong]

Actually idx is a cell matrix with different length.

I want const_z = [const_z, z(:,i) >= z(:,idx{i}(k))]; for k = 1:length(idx{i})

idx{1} = [1 3 5];

idx{2} = [3 5 10 .... 990];

...

idx{100} = [2 5 ... 500];

The longest idx{i} is 990 and the shortest one is 1.

Thanks for your help. 

[Johan Löfberg]

So the first column should be larger than column 1 3 and 5? Why is 1 included, column1 >= column1?

[jshong]

That's a good point. 

1 is not needed. 

Yes. First column should be larger than 3 and 5. 

And a similar constraint with different length of idx(i) goes until the last column which is 1000th. 

[Johan Löfberg]

So this would be the slow version

N = 5, M = 1000

% Random indicies

for i = 1:M

    idx{i} = randi(M,5,1);

end

z = binvar(N,M,'full');

const_z = [];

% Slow extraction of columns

tic

for i = 1:M

    for k = 1:length(idx{i})

        const_z = [const_z, z(:,i) >= z(:,idx{i}(k))];

    end

end

toc

The only way is to vectorize the indexations

tic

indiciesLeft = [];

indiciesRight = [];

for i =1:length(idx)

    for k = 1:length(idx{i})

        indiciesLeft = [indiciesLeft sub2ind([N M],1:N,repmat(i,1,N))];

        indiciesRight = [indiciesRight sub2ind([N M],1:N,repmat(idx{i}(k),1,N))];

    end

end

const_z = z(indiciesLeft) >= z(indiciesRight)

toc

[Johan Löfberg]

Eaiser

tic

ii = [];

kk = [];

for i = 1:M

    for k = 1:length(idx{i})

        ii = [ii i];

        kk = [kk idx{i}(k)];       

    end

end

const_z = [const_z, z(:,ii) >= z(:,kk)];

toc

[jshong]

Thanks. It helps me a lot.