Как векторизовать следующий блок Matlab

Question

Предположим, что A имеет границы (1: 2,1: 2,1: numfoo), как векторизовать следующие строки:

W = zeros( 2, 2, numfoo );
for i = 1:numfoo
    temp(1:2,1:2) = inv( A(1:2,1:2,i) );
    W(1:2,1:2,i) = ( temp *  (temp') );
end

TYIA!

Answer 1

Поскольку размер вашей матрицы составляет всего 2, вы можете использовать явные выражения для векторизации вашего кода. https://en.wikipedia.org/wiki/Inverse_of_a_matrix#Inversion_of_2.C3.972_matrices

dets=A(1,1,:).*A(2,2,:)-A(1,2,:).*A(2,1,:);
temp=[A(2,2,:)./dets -A(1,2,:)./dets ; -A(2,1,:)./dets A(1,1,:)./dets];
W=[temp(1,1,:).^2+temp(1,2,:).^2,...
   temp(1,1,:).*temp(2,1,:)+temp(1,2,:).*temp(2,2,:);...
   temp(2,1,:).*temp(1,1,:)+temp(2,2,:).*temp(1,2,:),...
   temp(2,1,:).^2+temp(2,2,:).^2];

Я проверил, он дает тот же результат , а ускорение x100

Elapsed time is 1.070547 seconds.
Elapsed time is 0.012767 seconds.

Answer 2

Я нашел много примеров того, что векторизованный код не всегда быстрее или понятнее.

Если вы хотите векторизованный код, используйте массив ячеек без цикла for, но он медленнее и не так понятен, как ваш.

Acell = mat2cell(A,2,2,ones(1,1,numfoo));
temp = cellfun(@inv,Acell,'UniformOutput',0);
W = cellfun(@(x,y)x*x', temp,'UniformOutput',0);
W = cell2mat(W);

Результаты такие же, как у вашего кода.

Answer 3

Небольшое увеличение скорости возможно при использовании mldivide вместо inv.

clear
clc
tic
numfoo=10000;
W = zeros( 2, 2, numfoo );
A=rand(2,2,numfoo);
for i = 1:numfoo
    temp(1:2,1:2) = inv( A(:,:,i) );
    W(1:2,1:2,i) = temp *temp';
end
toc

tic
for i = 1:numfoo
    temp(1:2,1:2) = A(:,:,i)\[1 0;0 1];
    W(1:2,1:2,i) = temp * temp' ;
end
toc

На моей машине вы получаете, Истекшее время составляет 0,182324 секунды.Прошедшее время составляет 0,162933 секунды.

Answer 4

Прежде всего, этот код будет значительно упрощен:

W = zeros( 2, 2, numfoo );
for i = 1:numfoo
    temp = inv( A(:,:,i) );
    W(:,:,i) = temp *  temp';
end

Во-вторых, что вы подразумеваете под "векторизацией следующих строк"? Какова ваша цель? Вы хотите избавиться от цикла FOR? Как вы думаете, этот код будет работать быстрее, если он будет «векторизован»? Что вы ищете?