question about vectorization using indexes

Hello, I am trying to do the following operations in matlab but I have a problem with how to properly write my code using vectorization. This is just an example, m, n and the values of the vectors and matrices are just to illustrate my problem. In reality m and n can go up to 1000.

n=5; m=8;
a=4*ones(m,1); a(2)=2;a(n)=3;
b=2*ones(n,2);b(1,1)=5;b(3,1)=1;
ind=3*ones(n,2); 
ind(1,2)=0;ind(3,2)=0; b(1,2)=0;b(3,2)=0;

non=zeros(1,n);c=zeros(1,n);

for i=1:n
    non(i)=nnz(ind(i,:));
    c(i)=prod(a(ind(i,1:non(i)))'.^b(i,1:non(i)),2);
end

I tried the following but it does not give correct results.

i=1:n;c=prod(a(ind(i,1:non(i))).^b(i,1:non(i)),2);

Thank you in advance

[Adam Danz]

What's the goal with this line of code?

 c(i)=prod(a(ind(i,1:non(i)))'.^b(i,1:non(i)),2);

I'm having a hard time guessing. The rest can easily be vectorized and cleaned up.

n=5; 
m=8;
a=4*ones(m,1); 
a([2,n])=[2,3];
b=2*ones(n,2);
b([1,3],1)=[5,1];
ind=3*ones(n,2);
ind([1,3],2)=0;
b([1,3],2)=0;
c=zeros(1,n);
non = sum(ind ~= 0,2)'; 
for i=1:n
    c(i)=prod(a(ind(i,1:non(i)))'.^b(i,1:non(i)),2); % ????
end

[AM]

That is the line I am trying to vectorize, I do not want to use a loop for it.

If

a =

     4
     2
     4
     4
     3
     4
     4
     4

>> b

b =

     5     0
     2     2
     1     0
     2     2
     2     2

>> ind

ind =

     3     0
     3     3
     3     0
     3     3
     3     3
     

For i=2 for example

     c(i)=prod(a(ind(i,1:non(i)))'.^b(i,1:non(i)),2);
     

does a(3)^b(2,1) *a(3)^b(2,2)=4^2 * 4^2

the matrix ind gives me the indexes of a I need to raise the power of b. Sometimes there is a 0 in ind so i determine the number of non zero elements and only take those values of a in each row.

I hope it's a bit more clear.

[Guillaume]

Note that :

    non = nnz(ind(i,:));
    x = a(ind(i,1:non);

can be written more simply as:

    x = a(nonzeros(ind(i, :)));

which is a lot easier to understand (particularly given the poorly named variables).

--edit:--

I reiterate Adam's question, what is the intent of the line

    c(i)=prod(a(ind(i,1:non(i)))'.^b(i,1:non(i)),2)

We now that it's what you want to vectorise. It'd be a lot easier to do if we knew what you're trying to do with it.

In particular, I'll point out that with the example given, the above will always pick element a(3), regardless of i.

[AM]

Thank you both for taking time to answer.

The values taken are just examples (poorly chosen I see but the main focus was that line where c is calculated) , what I am trying to do is avoid using a loop seeing as in reality n is equal to 1300. I thought using vectorization would make my code run faster since I am solving odes (I used profile viewer realized that part of the code is run over 8000 times).

What I want is to write a code that allows me to detemine c (i've changed ind so it doesn't always pick element 3) :

n=5; m=8;
a=4*ones(m,1); a([2,3,4,5,m])=[2,1,3,1,0];
b=2*ones(n,2);b([1,3],1)=[5,1];
ind=3*ones(n,2); ind([1,2,3,5],1)=[1 2 4 5];
ind([1,3],2)=0; b([1,3],2)=0;
c=zeros(1,n);
non = sum(ind ~= 0,2)';
for i=1:n
    c(i)=prod(a(ind(i,1:non(i)))'.^b(i,1:non(i)),2);
end

I do not know if there is a simpler way to calculate c, I see it's not clear to either of you what that line does so maybe I chose a convoluted way to calculate it.

Basically, with this new example

>> a

a =

     4
     2
     1
     3
     1
     4
     4
     0

>> ind

ind =

     1     0
     2     3
     4     0
     3     3
     5     3

>> b

b =

     5     0
     2     2
     1     0
     2     2
     2     2
     

c(1)=a(1)^b(1,1)

c(2)=a(2)^b(2,1)*a(3)^b(2,2)

c(3)=a(4)^b(3,1)

and so on

     c(i)=prod(a(ind(i,1:non(i)))'.^b(i,1:non(i)),2);

I use a(ind(i,1:non(i))) as the element of a because otherwise if i coded it as a(ind(i,:)) it would try to read a(0) for i=1 for example and give the following error

Array indices must be positive integers or logical values.

Is there a simpler way to do it and one that does not require a loop?

Answer 1

[Stephen23]

Note that ind and b must be transposed for this to work:

>> a   = [4;2;1;3;1;4;4;0]; % must be column!
>> ind = [1,0;2,3;4,0;3,3;5,3].'; % transposed!
>> b   = [5,0;2,2;1,0;2,2;2,2].'; % transposed!
>> idx  = b~=0;
>> XC = ind(idx);
>> bC =   b(idx);
>> [~,idc] = find(idx);
>> out = accumarray(idc,a(XC).^bC,[],@prod)
out =
  1024
  4
  3
  1
  1

Answer 2

[Guillaume]

Another option is to append a 0 (or any finite value) to the start of a and increase ind by 1, so a(ind+1) is always valid. Assuming that b is 0 when ind is 0 as in your example (if not, it's trivially fixed), then anything.^0 is 1 and multiplying by 1 doesn't affect the result, so:

apadded = [0; a];
c = prod(apadded(ind + 1) .^ b, 2)

As a bonus, c is a column vector matching the rows of b.

If b can be non-zero when ind is 0:

c = prod(apadded(ind + 1) . ^ (b .* (ind ~= 0)), 2)

to compensate.

edit: actually, if b can be non-zero when ind is 0, the easiest is to pad a with a 1 instead of a zero. Since 1.^anything is 1, it doesn't affect anything:

apadded = [1; a];
c = prod(apadded(ind + 1) .^b, 2)   %b can be zero or non-zero where ind is 0. It'll result in 1.^something

[Guillaume]

Note: although both options are very fast even for very large inputs, this option is about twice as fast as the accepted answer. And much simpler and using less memory.

Answer 3

[Andrei Bobrov]

ind(ind == 0) = 1;
c = prod(a(ind).^b,2);

[Guillaume]

Yes, that will also works as long as b is 0 when ind is 0, since a(1).^0 is 1.

[Andrei Bobrov]

lo = ind == 0;
ind(lo) = 1;
c = prod(a(ind).^(b.*~lo),2);