How do I preallocate memory for constraints

[Xu Wang84]

Hi there, 

When I add more and more constraints to the constraint vector, using format:

Constraint = [Constraint, a >= b];

it becomes time consuming, probably due to matrix resizing. So I want to pre-allocate memory for a large constraint vector (e.g., 10,000). How can I do that? 

What I have tried so far is: 

1. try to populate the last entry of the vector:

Constraint(1,10000) = [a >= b];

2. try to create a double vector and convert element to constraint:

Constraint = zeros(1,10000); Constraint(1,1) = [a >= b];

3. try to find some type identifier other than double:

Constraint = constraint(1,10000); 

Apparently, none works. Can somebody help me please? Thanks a lot.

Xu

[Johan Löfberg]

You cannot.

What does your code look like now?

[Xu Wang84]

The code constitutes several 'for' loop blocks. Within each block, there are many iterations, therefore generates many occurrence of resizing. I suppose this resizing is the main reason for long wait before the problem starts being solved. 

[Mark L. Stone]

This may not suit you needs, but if you have a large number of linear constraints, try to put them into a matrix formulation.  Fro example if x = sdpvar(n,1), for some n, and you have m constraints of the same type, formulate an m by n matrix A, and an n by 1 vector b, such that A*x == b, or A*x <= b or A*x >= b, as appropriate, and then specify that to YALMIP as your constraint. You can create A with spalloc if desired.  You can use vectorized assignments or even multiply nested for loops to quickly populate A and b.  Even if it is a little painful to figure out how to populate A and b the first time, once you do, you can can quickly generate A and b for other instances of the same problem. 

If m is very big, this will be much faster than having the constraint declarations processed individually through YALMIP.

[Mark L. Stone]

Upon seeing your attachment, it looks like it shouldn't be difficult to implement my recommendation above.

[Xu Wang84]

Good. I have thought about your method before, but ended up thinking that many constraints are limited to be element-wise. I will give it another try. Thanks. 

[Mark L. Stone]

Of course b should be m by 1, not n by 1 as I erroneously wrote.

Good Luck!!

[Johan Löfberg]

You should vectorize your code (as almost always in MATLAB)

Just picking a random example

%Spectral efficiency (inverse) limit
    for s = 1 : Nnodes - 1
        for d = s + 1 : Nnodes
            if Dsd(s,d) ~= 0,
                Constraints = [Constraints, SE(s,d) <= 1];
                Constraints = [Constraints, SE(s,d) >= 0.25];
            end
        end
    end

You could do (I am doing this in my head, most likely off slightly on the triu so you might need the second argument and tweak the code, but you get the idea)

index = ones(Nnodes);
index = find(triu(index) & Dsd~=0)
Constraints = [Constraints, SE(index) <= 1];


Profile your code and target the loop which takes the most time

BTW, why are you turning off the strict warnings, and use strict inequalities. The strict inequalities are treated as non-strict, there is no such thing as a strict inequality in YALMIP, and the possibility to turn off that warning is simply a patch for people with large old code bases.

[Xu Wang84]

Great, I will get back to you on my progress on reformatting the code. 

With regard to the turning off the strict warnings, I am not familiar with YALMIP, but what will happen if I have a constraint that requires strict inequality? (I will change the strict inequality in my constraints for now.) Also, is there side-effect if I turn the warning off?

[Johan Löfberg]

If x is integer, you write x<=1 instead of x<2 (of course,) and if x is continuous, you use x <= 2-smallnumber where smallnumber is a number you pick to make the constraint sufficiently (in your view) strictly satisfied, but still large enough for the solver to actually return a solution which is strict (the solver has a finite precision and thresholds so even if you say x<=0 you might get x=1e-12 etc)

The reason strict inequalities aren't used is, besides the fact that solver work with finite precision and thresholds etc., because they lead to ill-posed problems. As a trivial example, there is no solution to the problem min x s.t x>0

[Xu Wang84]

I see, thank you.

[Xu Wang84]

Hi guys, 

Thanks for the previous tips. I was able to vectorize my code: I get rid of all the for-loop in the constraints, and now, the creation of constraint matrix takes secs instead of minutes. But I also have a few new problem: 

1. It seems to me that if I declare a variable as a = sdpvar(3,3). The lower triangle of a has to be the same as the upper triangle (e.g. I can't get different values for a(1,3) and a(3,1)). Is this true? I tried a get-around by having two matrices one representing upper triangle and one representing lower triangle and use them in pair within the constraints. 

2. I was doing some matrix replicate using command 'repmat'. I want to expand a 2-D variable matrix to a 3-D variable matrix so that I can do some calculation with it. But it won't allow me to go over 2-D. Is this also true? What I did is to spread out the 3-D I was hoping to get to a 2-D version, which I expect to work but makes the code unclear. 

I attach my new file. If anyone can help me with the questions, I would be very appreciate that. 

Xu

[Johan Löfberg]

1. Yes. Square is symmetric by default

http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Tutorials.Basics

2. Don't understand the question. Give a small example of what you want to do, and what doesn't work

[Xu Wang84]

For second question, say, I have a variable matrix a = sdpvar(3,3), then I want to do following operation:

repmat(a, [1,1,3]); 

It won't let me. Error message says "N-dimensional indexing allowed for Full matrices only."

The reason I want to do that is because I have a constant 3*3*3 matrix (number 3 is just an example), and I want that 3-D variable matrix to create easy clean constraint expression. 

My get around is:

repmat(a, [1,3]);

which is doable, but requires a few more 'repmat' and 'reshape'. 

[Johan Löfberg]

Aha. Edit sdpvar/repmat and change the

catch
 error
end

with

catch


  X = varargin{1};

  Y = X;

  Y.basis = [];

  n = Y.dim(1);

  m = Y.dim(2);

  for i = 1:length(Y.lmi_variables)+1

    temp = repmatfixed(reshape(full(X.basis(:,i)),n,m),varargin{2:end});

    Y.basis(:,i) = temp(:);

  end

  Y.dim = size(temp); 

  Y = flush(Y);

  % Reset info about conic terms

  Y.conicinfo = [0 0];

  if length(temp)>2

      Y = ndsdpvar(Y);

  end

end

[Xu Wang84]

Works like a charm. Thanks.

[Mark L. Stone]

Leaving aside the general value of Yohan's fix to spdvar/repmat, if this a = sdpvar(3,3) is the same matrix variable which you want to be non-symmetric, then shouldn't it be declared as a full matrix  a=sdpvar(3,3,'full') , for which Yohan's fix is not needed?

[Xu Wang84]

Thank Yohan here first. Same question as Mark's. I thought declaring it in the form of a=sdpvar(3,3,'full') will save the trouble of using the fix. but it turned out not to be true.

[Johan Löfberg]

The repmat-bug is completely unrelated to symmetry of its argument. The problem was the step from 2D to 3D

[Pablo Luiz]

Hi Guys, 

I would like to now how can I improve this code?

          X = sdpvar(double(N));

           for i=1:N-2

               for j=i+1:N-1

                   for k=j+1:N

constraints= [constraints, X(i,j) + X(i,k) + X(j,k) >=-1, X(i,j) - X(i,k) - X(j,k) >= -1,-X(i,j) + X(i,k) - X(j,k) >= -1, -X(i,j) - X(i,k) + X(j,k) >= -1 ];
                     

                   end

               end 

           end

[Johan Löfberg]

Why don't you post this under a new topic? There is no such thing as pre-allocating memory for constraints.

What you want to do, as always in MATLAB; is to vectorize the code to get rid of the tripe for-loop

However, before any kind of optimization, I have to ask what N is? When I run this code with N=1000, which means you have a symmetric matrix with 500500 variables and thus can expect computational time in the order of minutes or hours, the setup of the constraints run in 1 second. Why even bother optimizing that?

BTW, what's with the double(N)? Why not simply sdpvar(N) (and note that the matrix is symmetric)

[Johan Löfberg]

Sorry, my timing was way off, so indeed you have to vectorize this

Trivial start is

for i=1:N-2

    for j=i+1:N-1

        k = j+1:N;        

        constraints= [constraints, X(i,j) + X(i,k) + X(j,k) >=-1, X(i,j) - X(i,k) - X(j,k) >= -1,-X(i,j) + X(i,k) - X(j,k) >= -1, -X(i,j) - X(i,k) + X(j,k) >= -1 ];                               

    end

end

However, you should do this more clever. I'm not saying this is exactly what you are doing, but it looks you are doing things like X + repmat(X(:,k),1,N) +  repmat(X(k,:),N,1) >= -1 etc.

[Pablo Luiz]

Im my computer (CPU intel core i5), If N is equal to 70 the setup of the constraints run in 48 second.

The value of N is assigned after reading file and the declaration of the variable X only worked with double (N)

[Johan Löfberg]

So what class does N have before casting it?

You will not be able to beat this, 0.1s for N=70. Completely unreadable though

N = 70;

X = sdpvar(N);

index = nchoosek(1:N,3);

use = index(:,1) <= N-2 & index(:,2) >= index(:,1)+1 & index(:,2) <= N-1 & index(:,3) >= index(:,2)+1;

index = index(find(use),:);

ij = sub2ind([N N],index(:,1),index(:,2));

ik = sub2ind([N N],index(:,1),index(:,3));

jk = sub2ind([N N],index(:,2),index(:,3));

tic

constraints=[X(ij) + X(ik) + X(jk) >=-1, X(ij) - X(ik) - X(jk) >= -1,-X(ij) + X(ik) - X(jk) >= -1, -X(ij) - X(ik) + X(jk) >= -1 ];

toc        

[Pablo Luiz]

Johan Löfberg 

I had managed to reduce the time to setup of the constraints, but not as much as your code. 

That's exactly what I would like to have accomplished.

Thank you very much.