векторизация вложенного цикла, где одна переменная цикла зависит от другой

Недавно я узнал, как векторизовать «простой» вложенный цикл в предыдущем вопрос, который я задал. Однако теперь я пытаюсь также векторизовать следующий цикл

A=rand(80,80,10,6,8,8);

I=rand(size(A1,3),1);
C=rand(size(A1,4),1);
B=rand(size(A1,5),1);

for i=1:numel(I)
    for v=1:numel(C)
        for j=1:numel(B)
            for k=1:j
                A(:,:,i,v,j,k)= A(:,:,i,v,j,k)*I(i)*C(v)*B(j)*((k-1>0)+1);               
            end
        end
    end
end

Итак, теперь k зависит от j... Что я пробовал до сих пор: комбинация терминов j и k (т.е. B(j)*((k-1>0)+1) дает треугольную матрицу, которую мне удается векторизовать независимо:

  B2=tril([ones(8,1)*B']');
  B2(2:end,2:end)=2*B2(2:end,2:end);

Но это дает мне матрицу (j,k) правильно, а не способ использовать ее для векторизации оставшегося цикла. Может быть, я тоже на неправильном пути... Итак, как я могу векторизовать этот тип петли?


for-loop matlab multidimensional-array vectorization
person Max    schedule 02.10.2014    source источник

Ответы (2)


arrow_upward
13
arrow_downward

В один из ваших комментариев к принятое решение предыдущего вопроса, вы упомянули, что последовательные коды на основе bsxfun(@times,..,permute..) были быстрее. Если это так, вы можете использовать аналогичный подход и здесь. Вот код, который использует такой шаблон вместе с tril -

B1 = tril(bsxfun(@times,B,[1 ones(1,numel(B)-1).*2]));
v1 = bsxfun(@times,B1, permute(C,[3 2 1]));
v2 = bsxfun(@times,v1, permute(I,[4 3 2 1]));
A = bsxfun(@times,A, permute(v2,[5 6 4 3 1 2]));
person Divakar    schedule 03.10.2014
comment
замечательно! это намного элегантнее и работает на 25% быстрее, чем решение @natan. - person Max; 06.10.2014
comment
@Макс Круто! Это тоже хорошо знать! - person Divakar; 06.10.2014
comment
Это решение напоминает мне Ramanujan. Я совершенно не понимаю, как, черт возьми, ты пришел к такому ответу. - person rayryeng; 26.11.2015

arrow_upward
2
arrow_downward

Вы были близки. Предложенная вами векторизация действительно следует логике (j,k), но выполнение tril добавляет нули в местах, куда не входит цикл. Использование решения для вашего предыдущего вопроса (@david's) не является полным, поскольку оно умножает все элементы, включая эти элементы с нулевым значением, в которые не входит цикл. Мое решение состоит в том, чтобы найти эти нулевые элементы и заменить их на 1 (так просто):

Начиная с вашего кода:

B2=tril([ones(8,1)*B']');
B2(2:end,2:end)=2*B2(2:end,2:end);

и после векторизации, показанной в предыдущем вопросе:

s=size(A);
[b,c,d]=ndgrid(I,C,B2);
F=b.*c.*d;
F(F==0)=1; % this is the step that is important for your case.
A=reshape(A,s(1),s(2),[]);
A=bsxfun(@times,A,permute(F(:),[3 2 1]));
A=reshape(A,s);

Для размера A, используемого в вопросе, это сокращает примерно 50% времени выполнения, неплохо...

person bla    schedule 02.10.2014