Как я уже упоминал в комментарии , векторизация больше не всегда является огромным преимуществом. Поэтому существуют методы векторизации, которые замедляют код, а не ускоряют его. Вы должны всегда рассчитывать свои решения. Векторизация часто включает создание больших временных массивов или копирование больших объемов данных, которых избегают в циклическом коде. Это зависит от архитектуры, размера ввода и многих других факторов, если такое решение будет быстрее.
Тем не менее, в этом случае кажется, что подходы векторизации могут дать большое ускорение.
Первое, что следует заметить в исходном коде, это то, что X(:, i, 1) .* X(:, j, 2) пересчитывается во внутреннем цикле, хотя это постоянное значение. Переписав внутренний цикл, это сэкономит время:
Y = X(:, i, 1) .* X(:, j, 2);
for k = 1:p
T(i, j, k) = sum(Y .* X(:, k, 3));
end
Теперь мы заметили, что внутренний цикл является точечным произведением и может быть записан следующим образом:
Y = X(:, i, 1) .* X(:, j, 2);
T(i, j, :) = Y.' * X(:, :, 3);
Транспонирование .' в Y не копирует данные, так как Y является вектором. Далее мы замечаем, что X(:, :, 3) индексируется повторно. Давайте переместим это из внешнего цикла. Теперь я остался со следующим кодом:
T = zeros(p, p, p);
X1 = X(:, :, 1);
X2 = X(:, :, 2);
X3 = X(:, :, 3);
for i = 1:p
for j = 1:p
Y = X1(:, i) .* X2(:, j);
T(i, j, :) = Y.' * X3;
end
end
Вероятно, что удаление цикла над j столь же легко, что оставило бы один цикл над i. Но на этом я и останавливаюсь.
Это время, которое я вижу (R2017a, 3-летний iMac с 4 ядрами). Для n=10, p=20:
original: 0.0206
moving Y out the inner loop: 0.0100
removing inner loop: 0.0016
moving indexing out of loops: 7.6294e-04
Luis' answer: 1.9196e-04
Для большего массива с n=50, p=100:
original: 2.9107
moving Y out the inner loop: 1.3488
removing inner loop: 0.0910
moving indexing out of loops: 0.0361
Luis' answer: 0.1417
"Ответ Луиса" - , этот . Это намного быстрее для небольших массивов, но для больших массивов он показывает стоимость перестановки. Перемещение вычисления первого продукта из внутреннего цикла экономит чуть более половины стоимости вычислений. Но удаление внутреннего цикла значительно снижает стоимость (чего я не ожидал, я предполагаю, что один матричный продукт может использовать параллелизм лучше, чем многие маленькие элементные продукты). Затем мы получаем дополнительное сокращение времени за счет сокращения количества операций индексирования в цикле.
Это временной код:
function so()
n = 10; p = 20;
%n = 50; p = 100;
X = randn(n,p,3);
T1 = method1(X);
T2 = method2(X);
T3 = method3(X);
T4 = method4(X);
T5 = method5(X);
assert(max(abs(T1(:)-T2(:)))<1e-13)
assert(max(abs(T1(:)-T3(:)))<1e-13)
assert(max(abs(T1(:)-T4(:)))<1e-13)
assert(max(abs(T1(:)-T5(:)))<1e-13)
timeit(@()method1(X))
timeit(@()method2(X))
timeit(@()method3(X))
timeit(@()method4(X))
timeit(@()method5(X))
function T = method1(X)
p = size(X,2);
T = zeros(p, p, p);
for i = 1:p
for j = 1:p
for k = 1:p
T(i, j, k) = sum(X(:, i, 1) .* X(:, j, 2) .* X(:, k, 3));
end
end
end
function T = method2(X)
p = size(X,2);
T = zeros(p, p, p);
for i = 1:p
for j = 1:p
Y = X(:, i, 1) .* X(:, j, 2);
for k = 1:p
T(i, j, k) = sum(Y .* X(:, k, 3));
end
end
end
function T = method3(X)
p = size(X,2);
T = zeros(p, p, p);
for i = 1:p
for j = 1:p
Y = X(:, i, 1) .* X(:, j, 2);
T(i, j, :) = Y.' * X(:, :, 3);
end
end
function T = method4(X)
p = size(X,2);
T = zeros(p, p, p);
X1 = X(:, :, 1);
X2 = X(:, :, 2);
X3 = X(:, :, 3);
for i = 1:p
for j = 1:p
Y = X1(:, i) .* X2(:, j);
T(i, j, :) = Y.' * X3;
end
end
function T = method5(X)
T = sum(permute(X(:,:,1), [2 4 5 3 1]) .* permute(X(:,:,2), [4 2 5 3 1]) .* permute(X(:,:,3), [4 5 2 3 1]), 5);