How to vectorize most constraint loops in CVX

1

Using permute and repmat, most constraint loops in CVX can be eliminated. For example, the following:

variables X(N) Y(M) Z(N);

for n=1:N
    for m=1:M
        if n~=m && n^3>m^2
           X(n)>=Y(m)+Z(n); % (1)
        end
    end
end

The first step is to expand. (1) actually means there are M*N constraints, so X is expanded to repmat(X,1,M), Y is expanded to repmat(Y,1,N), Z(n) is expanded to repmat(Z,1,M), according to their original dimensions. Now, repmat(X,1,M) is a N-by-M matrix, repmat(Y,1,N) is M-by-N, repmat(Z,1,M) is N-by-M. Note that some matrixs are M-by-N whlie others are N-by-M.

Before doing operations between them (in this case, add them), their dimensions need to be permuted to a uniform order. We choose the order N-by-M so that only repmat(Y,1,N) need to be permuted by using permute(repmat(Y,1,N),[2 1]). Then the whole loop can be easily replaced by the following

repmat(X,1,M) >= permute(repmat(Y,1,N),[2 1]) + repmat(Z,1,M); %(2)

if the condition “if n~=m && n^3>m^2” is removed.

To handle the condition “if n~=m && n^3>m^2”, we create a N-by-M filter matrix F, in which F(n,m)=1 if the condition is satisfied, otherwise 0. First of all, create a matrix for n using repmat((1:N)',1,M) , and for m using permute(repmat((1:M)',1,N),[2 1]), permute is used here to adjust order like above. Then obviously condition1 “n~=m” can be expressd as repmat((1:N)',1,M)~=permute(repmat((1:M)',1,N),[2 1]), condition2 “n^3>m^2” can be repmat((1:N)',1,M).^3>permute(repmat((1:M)',1,N),[2 1]).^2, the product of them is F

F=(repmat((1:N)',1,M)~=permute(repmat((1:M)',1,N),[2 1])).*(repmat((1:N)',1,M).^3>permute(repmat((1:M)',1,N),[2 1]).^2);

since it is a “and” relationship. To filter (2), we multiply it by F. So finally, the whole loop can be expressed as

F.*(repmat(X,1,M) - permute(repmat(Y,1,N),[2 1]) - repmat(Z,1,M))>=0;

With these three steps, I think most CVX constraint loops can be vectorized, so a lot of time can be saved especially when variables are too many. Obviously what’s described here can be easily generalized to other examples.

4 Likes
2

I declare @jackfsuia the Wizard of Vectorization.

3

It makes me recall my days as APL programmer.

4

It seems a good manner to save the running time of the CVX-based code. But I wonder whether the mentioned vectorization method is still feasible or not if the Y and Z become convex expressions? This case may be more meaningful for us.

5

The convex expressions does not refer to the linear terms, it could be the quadratic terms, logarithms, or other complex convex terms.

6

Yes, as long as the operator (the function) itself has been vectorized, and of course, they can be composite. For example, norm is vectorized as norms, inv_pos is vectorized, max is vectorized and some others. The idea above only applies to the loop of operators that have been vectorized by the developer , if it hasn’t, I don’t know what to do in general, so it depends on your problem. Besides, maybe you can

  1. Vectorize the operator yourself. The reshape command, which can reduce the dimensions, and some linear transform you do by hand might help. I don’t know.

  2. Use the idea to directly model the problem by yourself. I don’t know if this help, haven’t tried. If you do so, this book by MOSEK might help.

7

I would like to vectorize these constraints as much as possible. Ideally, full vectorization of everything, but the big bang for the buck would be to reduce the double for loop to no more than a single loop.

Input values for array sizes are provided just to make this code self-contained, but the actual values of n1 and n2 will be larger. Hopefully I haven’t made a mistake in the below code, but I certainly want to improve the formulation speed.

n1=3;n2=4;m=5;
cvx_begin
variable x_1(n1,m) binary
variable x_2(n2,m) binary
variable x(n1*n2,m) binary
for i=1:n1
  x((i-1)*n2+1:(i-1)*n2+n2,:) <= repmat(x_1(i,:),n2,1)
end
for j=1:n2
  x(j:n2:j+(n1-1)*n2,:) <= repmat(x_2(j,:),n1,1)
end
for i=1:n1
  for j=1:n2
    x((i-1)*n2+j,:) >= x_1(i,:)+x_2(j,:)-1
  end
end
% other portions of code not shown here
cvx_end

These are the linearization constraints for
x((i-1)*n2+j,:) = x1(i,:)*x2(j,:)
for from 1 to n1 and j from 1 to n2.
So there are n1*n2*m product terms to be linearized.

8

Please test the following 3 lines for the corresponding 3 loops you have, they might work

x <= reshape(permute(repmat(x_1,1,1,n2),[3 1 2]),n1*n2,m);

reshape(x,n2,n1,m) <= permute(repmat(x_2,1,1,n1),[1 3 2]);

permute(reshape(x,n2,n1,m),[2 1 3])>=permute(repmat(x_1,1,1,n2),[1 3 2])+permute(repmat(x_2,1,1,n1),[3 1 2])-1;

The definition of the x makes it not easy, why not just variable x(n1,n2,m) binary or something like that?

9

Thanks.

I think your proposed constraints might be correct, but have not put in the effort to verify. That is because your suggestion of using a 3D binary variable is a good one, and results in a better and easier to understand formulation. Fortunately .* works when multiplying 3D matrices element by element.

10

A higher-than-two-dimension example is here:


some of other examples are here:

12

Here are some explanations for permute and repmat, which I thought to be unneccessary, since they are a Matlab thing, people can just check its own website and easily learn (This thought seems to be just wrong), and my english is not good enough.

They are important and can be used to vectorize arbitrary layers of loops as I have showed. permute is adjusting the order of dimensions. For example, array X(I,J,K,N), if you want it to become X(N,K,J,I) in the second step above, use permute(X, [4 3 2 1]). Yes, it is actually simple, 4 is the original position of the dimension N in X(I,J,K,N), 3 is that of K, 2 is that of J, 1 is that of I. repmat is for expansion, if you have X(A,B,C), but you need X(A,B,C,D,E,F,G), use repmat(X, 1,1,1,D,E,F,G). Yes, the beginning three ones correspond to the dimensions you dont want to change. the rest DEFG is just DEFG in its order.OK, what if you need X(E,D,F,G,A,C,B), but given the X(A,B,C)?? Yes, repmat, then permute.

I haven’t used Matlab/CVX for a while since I graduated, hope these explanations make sense.

1 Like