question about vectorization using indexes

1 view (last 30 days)
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
  4 Comments
Guillaume
Guillaume on 30 Jul 2019
Edited: Guillaume on 30 Jul 2019
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
AM on 30 Jul 2019
Edited: AM on 30 Jul 2019
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?

Sign in to comment.

Accepted Answer

Stephen23
Stephen23 on 30 Jul 2019
Edited: Stephen23 on 30 Jul 2019
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

More Answers (2)

Guillaume
Guillaume on 30 Jul 2019
Edited: Guillaume on 30 Jul 2019
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
  1 Comment
Guillaume
Guillaume on 31 Jul 2019
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.

Sign in to comment.


Andrei Bobrov
Andrei Bobrov on 31 Jul 2019
ind(ind == 0) = 1;
c = prod(a(ind).^b,2);
  2 Comments
Guillaume
Guillaume on 31 Jul 2019
Yes, that will also works as long as b is 0 when ind is 0, since a(1).^0 is 1.
Andrei Bobrov
Andrei Bobrov on 31 Jul 2019
lo = ind == 0;
ind(lo) = 1;
c = prod(a(ind).^(b.*~lo),2);

Sign in to comment.

Categories

Find more on Programming in Help Center and File Exchange

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!