question about vectorization using indexes

AM
2019 7 30
3 个回答
回答已采纳
更新时间：2019 7 31
1 次查看（30 天）

0 个投票

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更早的评论 隐藏 2更早的评论

Adam Danz
Adam Danz 2019-7-30
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
AM 2019-7-30
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
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

0 个投票

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

0 个投票

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 个评论
显示 -1更早的评论 隐藏 -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

0 个投票

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.
Andrei Bobrov
Andrei Bobrov 2019-7-31
lo = ind == 0;
ind(lo) = 1;
c = prod(a(ind).^(b.*~lo),2);

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Creating and Concatenating Matrices 的更多信息

产品

标签

Community Treasure Hunt

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

Start Hunting!

Translated by