The function below may solve part of my probelm. It is named "mprod" vs. prod, similar to times vs. mtimes. With some reshaping, it uses multiprod recursively. In general, a recursive function call is slower than a loop. Multiprod claims to be >100x faster, so it should more than compensate.
function sqMat = mprod(M)
% Multiply *many* square matrices together, stored
% as 3D array M. Speed gain through recursive use
% of function 'multiprod' (Leva, 2010).
% check if M consists of multiple matrices
if size(M,3) > 1
% check for odd number of matrices
if mod(size(M,3),2)
siz = size(M,1);
M = cat(3,M,eye(siz));
end
% create two smaller 3D arrays
X = M(:,:,1:2:end); % odd pages
Y = M(:,:,2:2:end); % even pages
% recursive call
sqMat = mprod(multiprod(X,Y));
else
% create final 2D matrix and break recursion
sqMat = M(:,:,1);
end
end
I have not tested this function for speed or accuracy. I believe this is much faster than a loop. It does not 'vectorize' the operation since it cannot be used with higher dimensions; any repeated use of this function must be done within a loop.
EDIT Below is new code that seems to work fast enough. Recursive calls to functions are slow and eat up stack memory. Still contains a loop, but reduces the number of loops by log(n)/log(2). Also, added support for more dimensions.
function sqMats = mprod(M)
% Multiply *many* square matrices together, stored along 3rd axis.
% Extra dimensions are conserved; use 'permute' to change axes of "M".
% Speed gained by recursive use of 'multiprod' (Leva, 2010).
% save extra dimensions, then reshape
dims = size(M);
M = reshape(M,dims(1),dims(2),dims(3),[]);
extraDim = size(M,4);
% Check if M consists of multiple matrices...
% split into two sets and multiply using multiprod, recursively
siz = size(M,3);
while siz > 1
% check for odd number of matrices
if mod(siz,2)
addOn = repmat(eye(size(M,1)),[1,1,1,extraDim]);
M = cat(3,M,addOn);
end
% create two smaller 3D arrays
X = M(:,:,1:2:end,:); % odd pages
Y = M(:,:,2:2:end,:); % even pages
% recursive call and actual matrix multiplication
M = multiprod(X,Y);
siz = size(M,3);
end
% reshape to original dimensions, minus the third axis.
dims(3) = [];
sqMats = reshape(M,dims);
end
bsxfun
. – paddyA = [a1 a2 a3 ... am]
, then usingB = A(:,:,i)*B
will produceam*am-1*...*a2*a1
, while usingB = B*A(:,:,i)
would producea1*a2*a3*...*am
. Matrix products are non-commutative, so these outcomes are generally different. Which one do you want? – Rody Oldenhuis