2 Loops of same size, 1 goes very quick but one is very slow.

Hello,

I have 2 loops with 4 index each. Every loop is of this form

X1 =zeros (10,11);
s= zeros(1,10^2*11^2);
for i = 1:10
for j = 1:10
for k = 1:11
for l = 1:11
    
counter = counter + 1
s(counter) = 2;
X1(i,k) = X1(i,k) - 1;
X1(j,l) = X1(j,l) - 1;
X = reshape(X1',1,[]);
T1 = {s,X} ;                 
V1 = horzcat(T1{:}) ;        

A(counter,:) = V1(:);
X1 =zeros(10,11);
s= zeros(1,10^2*11^2);

end
end
end
end

The first loop does what it should do very quickly. And the second one is exactly of the same form but this runs very slowly.

I watch the variable "counter" to see how fast it runs and the first loop is finished after a few seconds, but the second one needs already 2 minutes for the first 500 cycles.

Why is it so slow and what can I do to make it faster ?

Best regards !

[Guillaume]

Why don't you explain what you're trying to do with that code?

As far as I can tell, it's just a massively overcomplicated way of creating vectors full of zeros with just one element set to 2 and two elements set to -1.

One example of overcomplication:

T1 = {s,X} ;                 
V1 = horzcat(T1{:}) ;

this is just:

V1 = [s, X];

Similarly,

X1(i,k) = X1(i,k) - 1;
X1(j,l) = X1(j,l) - 1;

is simply:

X1(i,k) = -1;
X1(j,l) = -1;

since X1 is always a matrix of 0s before that.

[Kernel7364]

I want to create a binary linear optimization problem and this will be the matrix for the constraints Ay . So y is a vector and I have constraints to s and x and thats why I want to create this matrix.

Thank you very much for the advice with V1.

But your second advice is not what I want. It has to be -2 if i=j=k=l and thats why I have done it this way. My other loop is of the same form, but just with other values for x1(i,j) and s. It looks like this:

X1 =zeros (10,11);
s= zeros(1,10^2*11^2);
for i = 1:10
for j = 1:10
for k = 1:11
for l = 1:11
    
counter = counter + 1
s(counter) = -1;
X1(i,k) = X1(i,k) + 1;
X1(j,l) = X1(j,l) + 1;
X = reshape(X1',1,[]);                
V1 = [s,X] ;        
A(counter,:) = V1(:);
X1 =zeros(10,11);
s= zeros(1,10^2*11^2);
end
end
end
end

I am just a beginner with Matlab, so please excuse me if I have very easy and big mistakes. I know that there are many zeros. But in the first loop which I postet first it takes Matlab just a few seconds. And I dont understand why the second one takes so long. I even did this :

clearvars i j k l s X1 T1 V1

to make it faster but that didnt work too.

Thank you very much for your quick answer !

[Guillaume]

Are you sure that counter is set to 0 before you start any of these loops? If not, you'll be growing s slowly in the loops to sizes far greater than 12100. I note that there are no

counter = 0;

in any of the code you've posted.

But your second advice is not what I want

The whole thing can be constructed without any loop whatsoever, so please explain exactly what your constraints are. In particular, at the moment your inner loops iterate over columns and the outer loops iterate over rows, which is less efficient than iterating over the rows in the inner loops. So is the order of iteration critical?

[Kernel7364]

Sorry for my mistake, i set the counter = 0 in front of the loop, but I forgot it to copy in the code to my question.

I want to minimize where

I have a transformation :

All in all I have to minimize my vector , since I have constraints to s and x.

Then i have the following constraints:

The first loop that I postet in my question is to fulfil the first to constraints to .

The second loop is to fulfil the last inequality.

Furthermore I have a constraint to my . But these one is at the end of my matrix A and it is very easy and fast to generate. This are just 10 more rows in my matrix A.

I had no other idea to make the matrix A since the constraints should hold to every and I have 10*10*11*11 of these variables and I have 3 constraints to it.

Do you have an idea to create the matrix A faster.

I am very glad that you try to help me. Thank you very much !!!

Answer 1

[Guillaume]

Here is a simple way to create your A (the one with +1 and +2)

[v1, v2] = ndgrid(num2cell(eye(110), 2));  %110 is 10*11
vsum = cellfun(@plus, v1, v2, 'UniformOutput', false);
vsum = cell2mat(vsum(:));
Apos = [2*eye(size(vsum, 1)), vsum];

Your A with -1 and -2 is simply:

Aneg = [2*eye(size(vsum, 1)), -vsum];

Note that the order of the rows is going to be different from what your loops generate but I assume it doesn't matter.

[Kernel7364]

This is absolutely brilliant. It is very good and fast, thank you very much.

There is just one problem, the algorithm can not found a solution, because the order how I have done it is important.

Because my vector is my variable which will be optimized.

And it is important that the first variable is only 1 if .

DO you know what I mean ?

Because and thats why it is important that s is 1 if and only if . I hope that you undersrtand what I tried to explain.

[Kernel7364]

For example if I want to minimize

I am in the 8. row if i=1, j=1, k=3, l=2.

So it has to be that

And if then s(8) hast to be 1.

The 8. row in my matrix Aneg should look like

to make sure that the substitution holds.

But at the moment with your construction it is

and that would basicly mean that

but that is not correct. I hope that you know what I mean and hope that you can help me. The first part with the generation of A is brilliant !!!

[Guillaume]

I don't believe the order of the rows in A is important. Your original code consider i and j as the row index of a 10x11 matrix and k and l as the column index.

My implementation swap them, with i and j as the column index, and k and l as the row index, because of the way linear indexing works in matlab.

If you want s to be 1 instead of 2 when there's a 2 (or -2) in vsum (so when ), then:

[v1, v2] = ndgrid(num2cell(eye(110), 2));  %110 is 10*11
vsum = cellfun(@plus, v1, v2, 'UniformOutput', false);
vsum = cell2mat(vsum(:));
s = eye(size(vsum, 1)) .* (2 - any(abs(vsum) == 2, 2));
Apos = [s, vsum];
Aneg = [s, -vsum];

[Kernel7364]

But when I try to solve it, it says that there is no inner point with my constraints. When I generated the matrix in my loops it gave me the correct solution. Is it possible to generate the rows in the matrix in the order I did it ?

The rest you did with your first code is absolutely briliiant and works very good. But I believe the order of the rows is important because otherwise the substitution is not valid anymore when I change the order.

I am very sorry if I did not explain it in a good way or does not understand it because I am just a beginner with matlab.

Thank you so much for your help !

If it is possible to generate the part in A where the x are is like this

X1(i,k) = X1(i,k) + 1;
X1(j,l) = X1(j,l) + 1;

in the right order of the index you would be the hero of my week ! That would be brilliant.

[Kernel7364]

I mean the part in A of x has completly the same structur in Apos and Aneg, just with different numbers (in Apos with 1 and 2 and in Aneg with -1 and -2). I just have to generate it once like you did it and than put an minus in front of this part in Aneg and thats it. It should not be that much effort for Matlab hopefully

[Guillaume]

I'm afraid I don't understand everything you ask.

To change the order of the rows of vsum to reflect your original order, you could do this:

vsum2 = reshape(vsum.', 11, 10, []);  %vsum as calculated before
vsum2 = permute(vsum2, [2, 1, 3]);
vsum2 - reshape(vsum2, 110, []).';

vsum2 will have exactly the same order as your loop generated.

[Kernel7364]

I did it with smaller dimensions not 11 and 10 but 2 and 3. So i goes to 2 and j goes to 2 and k goes to 3 and l goes to 3. So in every row of my Matrix A are 36 elements from s (one is 2 and the rest is zero) and the last elements describe the .

and if I do it your way I have

-2	0	0	0	0	0
-1	0	-1	0	0	0
-1	0	0	0	-1	0
-1	-1	0	0	0	0
-1	0	0	-1	0	0
-1	0	0	0	0	-1
-1	0	-1	0	0	0
0	0	-2	0	0	0

But when I do it with the loop I did it before it looks like this

-2	0	0	0	0	0
-1	-1	0	0	0	0
-1	0	-1	0	0	0
-1	-1	0	0	0	0
0	-2	0	0	0	0
0	-1	-1	0	0	0
-1	0	-1	0	0	0
0	-1	-1	0	0	0

I am really bad in explaining but I hope that you know what I mean.

Is it possible to get it in this structure ? If not I have to generate it with my loops, but if you know a tric it would be amazing. In my code here I just showed the first 8 rows as an example.

Thank you very much for your help.

[Kernel7364]

If my index i and j goes to 2 and k and l to 3 this is my matrix when I generate it with my loops.

This is just the part of your vsum. The part where the 2 is is correct and the rest should look like this:

-2	0	0	0	0	0
-1	-1	0	0	0	0
-1	0	-1	0	0	0
-1	-1	0	0	0	0
0	-2	0	0	0	0
0	-1	-1	0	0	0
-1	0	-1	0	0	0
0	-1	-1	0	0	0
0	0	-2	0	0	0
-1	0	0	-1	0	0
-1	0	0	0	-1	0
-1	0	0	0	0	-1
0	-1	0	-1	0	0
0	-1	0	0	-1	0
0	-1	0	0	0	-1
0	0	-1	-1	0	0
0	0	-1	0	-1	0
0	0	-1	0	0	-1
-1	0	0	-1	0	0
0	-1	0	-1	0	0
0	0	-1	-1	0	0
-1	0	0	0	-1	0
0	-1	0	0	-1	0
0	0	-1	0	-1	0
-1	0	0	0	0	-1
0	-1	0	0	0	-1
0	0	-1	0	0	-1
0	0	0	-2	0	0
0	0	0	-1	-1	0
0	0	0	-1	0	-1
0	0	0	-1	-1	0
0	0	0	0	-2	0
0	0	0	0	-1	-1
0	0	0	-1	0	-1
0	0	0	0	-1	-1
0	0	0	0	0	-2

[Guillaume]

Ok, I completely missed that the two inner loops were both iterating over the same coordinate. There's probably a way to generate the result in that order with my original algorithm with some reshape and permute but I can't be bothered to figure it out. So new tack:

numrows = 3;  %your maximum k and l
numcols = 2;  %your maximum i and j

destcoords = repelem({1:numrows, 1:numcols}, 2); %ranges of k, l, i, j
[destcoords{:}] = ndgrid(destcoords{:});         %cartesian product of ranges in the exact order you want
destcol1 = sub2ind([numrows, numcols], destcoords{1}, destcoords{3});  %convert coordinates into linear indices. This is going to be the destination column in A
destcol2 = sub2ind([numrows, numcols], destcoords{2}, destcoords{4});
destrow = repmat(1:numel(destcol1), 1, 2).';   %destination row in A
Aneg = accumarray([destrow, [destcol1(:); destcol2(:)]], -1);
Apos = accumarray([destrow, [destcol1(:); destcol2(:)]], 1);

S can be constructed the same as before, with eye.

[Kernel7364]

Thank you very much, it works perfectly !!!

That helped me so much.

Best wishes and have a nice day !

[Kernel7364]

The next problem I have is to solve this optimization problem... The Cplex solver needed a long time and if I do it for a bigger one it did not solve it after 2 hours.

Do you know how cplex.solve() works ?

I read that I can also use this one and this function should be good for optimization problems.

But it did not work because it says something that

The class Cplex has no Constant property or Static method named 'solve'.

I have no idea whats wrong. I can use other function from cplex, but not cplex.solve()

Maybe you know something about cplex.solve() and you can rescue me again.

Best regards !

[Guillaume]

The Cplex solver needed a long time and if I do it for a bigger one it did not solve it after 2 hours

I'm afraid this is not something I can help with. It's completely outside my domain of expertise, I don't even know what CPlex is. However, I wouldn't be surprised if it took a long time to find a solution to whatever it is you're trying to solve. You have 12100 unknown!

The class Cplex has no Constant property or Static method named 'solve'.

Sounds like solve is a normal method of the class, not a static method. If so, you need to call it on an instance of the class, not on the class itself. Something like:

myobj = CPlex(???);  %construct an instance of the CPlex class. Replace the ??? by the required inputs if any.
myobj.solve;   %call solve on the instance.

[Kernel7364]

Okay, thank you very much for all your help !

I will try it tomorrow and see if I can do something. I hop that Cplex can solve it a bit faster but I can imagine that it will take quiet a while since there are more than 12000 unknown.

Thank you so much for all your help !