3
votes

I'm trying to build a constraint matrix that I'm going to use with linprog, and I'm struggling to build it efficiently without using for loops. Here is an example of what I'm trying to achieve:

A = zeros(3, 3); % Constraint matrix
i = [1 3]; % row indexes
j = [1 2 ; 1 3]; % column indexes
x = [1 2 ; 3 4]; % values to assign

Expected result after assignment:

A = [1 2 0 ; 0 0 0 ; 3 0 4]

I want to perform the following operations:

A(i(1), j(1,:)) = x(1,:)
A(i(2), j(2,:)) = x(2,:)

For now, I'm using a for loop:

for k=1:length(i)
    A(i(k), j(k,:)) = x(k,:);
end

Is there any better way to do this for any i and j? I could use a for loop everywhere, but the number of constraint depends on the number of variables, so my code gets filled with for-loops. What's the standard way of defining the constraint matrix to use with linprog?

3

3 Answers

4
votes

What you are looking for is the sub2ind function. Here's how you would improve your matrix creations:

>> indx = sub2ind(size(A),[1 3 1 3]',j(:))
indx =

   1
   3
   4
   9

>> A(indx)=x(:)
A =

   1   2   0
   0   0   0
   3   0   4

Note that you would have to tweak your i definition a little bit so i and j have the same number of elements.

3
votes

One vectorized approach with bsxfun -

A(bsxfun(@plus,ii(:),(jj-1)*size(A,1))) = x

You need the expansion with bsxfun as the number of row indices don't match up with the number of column indices.

Also, please note that I have replaced the variable names i with ii and j with jj as i and j are also used with complex numbers, as that might cause some conflict otherwise.

Sample run -

>> ii(:)    %// Row indices
ans =
     1
     3
>> jj       %// Column indices
jj =
     1     2     4
     1     3     5
>> x        %// Values to assign
x =
     1     2     6
     3     4     8
>> A        %// Output
A =
     1     2     0     6     0
     0     0     0     0     0
     3     0     4     0     8

It has two benefits:

  1. Avoiding the call to sub2ind with a raw-version of it, for being more efficient with runtime.

  2. The expansion being done internally without using any loop or repmat saves looping or another function call respectively.

1
votes

Here's another way, making use of sparse:

A = full(sparse(repmat(ii,size(jj,1),1).', jj ,x));

I'm using ii and jj as variable names instead of i and j, as in Divakar's answer.