question about vectorization using indexes

1 次查看（过去 30 天）
AM
AM 2019-7-30
评论： Andrei Bobrov 2019-7-31
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 个评论
显示 2更早的评论
Guillaume
Guillaume 2019-7-30
编辑：Guillaume 2019-7-30
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 2019-7-30
编辑：AM 2019-7-30
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?

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Stephen23
Stephen23 2019-7-30
编辑：Stephen23 2019-7-30
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

更多回答（2 个）

Guillaume
Guillaume 2019-7-30
编辑：Guillaume 2019-7-30
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 个评论
Guillaume
Guillaume 2019-7-31
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.

请先登录，再进行评论。

Andrei Bobrov
Andrei Bobrov 2019-7-31
ind(ind == 0) = 1;
c = prod(a(ind).^b,2);
  2 个评论
Guillaume
Guillaume 2019-7-31
Yes, that will also works as long as b is 0 when ind is 0, since a(1).^0 is 1.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Function Creation 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by