НОВОЕ РЕШЕНИЕ:
Это более новое решение, построенное на более старом решении, которое решило ранее заданную формулу.Код в вопросе на самом деле является модификацией этой формулы, в которой перекрытие между двумя матрицами в третьем измерении многократно смещается (это похоже на свертку в третьем измерении данных).Предыдущее решение, которое я дал, вычисляло только результат для последней итерации кода в вопросе (т.е. k = kend).Итак, вот полное решение, которое должно быть намного более эффективным, чем код в вопросе для kend порядка 1000:
kend = size(A,3); %# Get the value for kend
C = zeros(3,3,kend); %# Preallocate the output
Anew = reshape(flipdim(A,3),3,[]); %# Reshape A into a 3-by-3*kend matrix
Bnew = reshape(permute(B,[1 3 2]),[],3); %# Reshape B into a 3*kend-by-3 matrix
for k = 1:kend
C(:,:,k) = Anew(:,3*(kend-k)+1:end)*Bnew(1:3*k,:); %# Index Anew and Bnew so
end %# they overlap in steps
%# of three
Даже при использовании только kend = 100для меня это решение оказалось примерно в 30 раз быстрее, чем рассматриваемое, и примерно в 4 раза быстрее, чем решение на основе чистых циклов (которое включало бы 5 циклов !).Обратите внимание, что приведенное ниже обсуждение точности с плавающей точкой по-прежнему применимо, поэтому вполне нормально и ожидается, что вы увидите небольшие различия между решениями порядка относительной точности с плавающей точкой .
СТАРЫЕ РЕШЕНИЯ:
Исходя из этой формулы, на которую вы ссылаетесь в комментарии:
кажется, что вы на самом делехочу сделать что-то отличное от кода, который вы указали в вопросе.Предполагая, что A и B являются матрицами 3 на 3 на k, результат C должен быть матрицей 3 на 3, а формула из вашей ссылки, записанная как набор вложенных циклов for, будетвыглядит следующим образом:
%# Solution #1: for loops
k = size(A,3);
C = zeros(3);
for i = 1:3
for j = 1:3
for r = 1:3
for l = 0:k-1
C(i,j) = C(i,j) + A(i,r,k-l)*B(r,j,l+1);
end
end
end
end
Теперь можно выполнить эту операцию без циклов for, изменив и реорганизовав A и B соответствующим образом:
%# Solution #2: matrix multiply
Anew = reshape(flipdim(A,3),3,[]); %# Create a 3-by-3*k matrix
Bnew = reshape(permute(B,[1 3 2]),[],3); %# Create a 3*k-by-3 matrix
C = Anew*Bnew; %# Perform a single matrix multiply
Вы могли бы даже переработать код, который у вас есть в вашем вопросе, чтобы создать решение с одним циклом, который выполняет матричное умножение ваших подматриц 3 на 3:
%# Solution #3: mixed (loop and matrix multiplication)
k = size(A,3);
C = zeros(3);
for l = 0:k-1
C = C + A(:,:,k-l)*B(:,:,l+1);
end
Итак, теперь вопросКакой из этих подходов быстрее / чище?
Ну, «чище» очень субъективно, и я, честно говоря, не могу сказать вам, какой из приведенных выше фрагментов кода облегчает понимание того, что за операцияделается.Все циклы и переменные в первом решении усложняют отслеживание происходящего, но оно явно отражает формулу.Второе решение разбивает все это на простую матричную операцию, но трудно понять, как это соотносится с исходной формулой.Третье решение похоже на золотую середину между этими двумя.
Итак, давайте сделаем скорость тай-брейка.Если я рассчитываю вышеупомянутые решения для ряда значений k, я получаю эти результаты (в секундах, необходимых для выполнения 10 000 итераций данного решения, MATLAB R2010b):
k | loop | matrix multiply | mixed
-----+--------+-----------------+--------
5 | 0.0915 | 0.3242 | 0.1657
10 | 0.1094 | 0.3093 | 0.2981
20 | 0.1674 | 0.3301 | 0.5838
50 | 0.3181 | 0.3737 | 1.3585
100 | 0.5800 | 0.4131 | 2.7311 * The matrix multiply is now fastest
200 | 1.2859 | 0.5538 | 5.9280
Ну, получаетсяиз этого следует, что при меньших значениях k (около 50 или менее) решение for-loop на самом деле выигрывает, еще раз демонстрируя, что циклы for не так «злы», как это считалось в более старых версиях MATLAB.При определенных обстоятельствах они могут быть более эффективными, чем умная векторизация.Однако, когда значение k больше, чем около 100, векторизованное решение умножения матриц начинает выигрывать, масштабируясь с увеличением k намного лучше, чем решение для цикла.Смешанное решение для цикла / матрицы-умножения масштабно жестоко по причинам, в которых я не совсем уверен.
Так что, если вы ожидаете, что k будет большим, я быпойти с векторизованным решением умножения матрицы.Следует иметь в виду, что результаты, которые вы получаете от каждого решения (матрица C), будут незначительно отличаться (на уровне точности с плавающей запятой ) с момента добавления иумножения, выполняемые для каждого решения, различны, что приводит к разнице в накоплении ошибок округления .Короче говоря, разница между результатами для этих решений должна быть незначительной, но вы должны знать об этом.