这是您发布的解决方案的稍微完善的版本，并进行了一些小的改进。
我们检查行数是否多于列数或相反，然后通过选择行与矩阵相乘或矩阵与列相乘来进行相应的乘法运算（从而执行最少的循环迭代）。
注意：这可能并不总是最好的策略（按行而不是按列），即使行数少于列数； MATLAB 数组存储在内存中的 column-major order 中这一事实使得按列切片更有效，因为元素是连续存储的。而访问行涉及通过strides 遍历元素（这对缓存不友好——想想spatial locality）。
除此之外，代码应处理双/单、实数/复数、完整/稀疏（以及不可能组合的错误）。它还尊重空矩阵和零维。
function C = my_mtimes(A, B, outFcn, inFcn)
% default arguments
if nargin < 4, inFcn = @times; end
if nargin < 3, outFcn = @sum; end
% check valid input
assert(ismatrix(A) && ismatrix(B), 'Inputs must be 2D matrices.');
assert(isequal(size(A,2),size(B,1)),'Inner matrix dimensions must agree.');
assert(isa(inFcn,'function_handle') && isa(outFcn,'function_handle'), ...
'Expecting function handles.')
% preallocate output matrix
M = size(A,1);
N = size(B,2);
if issparse(A)
args = {'like',A};
elseif issparse(B)
args = {'like',B};
else
args = {superiorfloat(A,B)};
end
C = zeros(M,N, args{:});
% compute matrix multiplication
% http://en.wikipedia.org/wiki/Matrix_multiplication#Inner_product
if M < N
% concatenation of products of row vectors with matrices
% A*B = [a_1*B ; a_2*B ; ... ; a_m*B]
for m=1:M
%C(m,:) = A(m,:) * B;
%C(m,:) = sum(bsxfun(@times, A(m,:)', B), 1);
C(m,:) = outFcn(bsxfun(inFcn, A(m,:)', B), 1);
end
else
% concatenation of products of matrices with column vectors
% A*B = [A*b_1 , A*b_2 , ... , A*b_n]
for n=1:N
%C(:,n) = A * B(:,n);
%C(:,n) = sum(bsxfun(@times, A, B(:,n)'), 2);
C(:,n) = outFcn(bsxfun(inFcn, A, B(:,n)'), 2);
end
end
end
比较
这个函数无疑是整个过程都比较慢，但是对于更大的尺寸，它比内置的矩阵乘法差几个数量级：
(tic/toc times in seconds)
(tested in R2014a on Windows 8)
size mtimes my_mtimes
____ __________ _________
400 0.0026398 0.20282
600 0.012039 0.68471
800 0.014571 1.6922
1000 0.026645 3.5107
2000 0.20204 28.76
4000 1.5578 221.51
这里是测试代码：
sz = [10:10:100 200:200:1000 2000 4000];
t = zeros(numel(sz),2);
for i=1:numel(sz)
n = sz(i); disp(n)
A = rand(n,n);
B = rand(n,n);
tic
C = A*B;
t(i,1) = toc;
tic
D = my_mtimes(A,B);
t(i,2) = toc;
assert(norm(C-D) < 1e-6)
clear A B C D
end
semilogy(sz, t*1000, '.-')
legend({'mtimes','my_mtimes'}, 'Interpreter','none', 'Location','NorthWest')
xlabel('Size N'), ylabel('Time [msec]'), title('Matrix Multiplication')
axis tight
额外
为了完整起见，下面是实现广义矩阵乘法的两种更简单的方法（如果要比较性能，请将my_mtimes 函数的最后一部分替换为其中任何一种）。我什至不会费心发布他们经过的时间:)
C = zeros(M,N, args{:});
for m=1:M
for n=1:N
%C(m,n) = A(m,:) * B(:,n);
%C(m,n) = sum(bsxfun(@times, A(m,:)', B(:,n)));
C(m,n) = outFcn(bsxfun(inFcn, A(m,:)', B(:,n)));
end
end
另一种方式（使用三重循环）：
C = zeros(M,N, args{:});
P = size(A,2); % = size(B,1);
for m=1:M
for n=1:N
for p=1:P
%C(m,n) = C(m,n) + A(m,p)*B(p,n);
%C(m,n) = plus(C(m,n), times(A(m,p),B(p,n)));
C(m,n) = outFcn([C(m,n) inFcn(A(m,p),B(p,n))]);
end
end
end
接下来要尝试什么？
如果您想获得更多性能，您将不得不迁移到 C/C++ MEX 文件以减少解释 MATLAB 代码的开销。您仍然可以通过从 MEX 文件中调用优化的 BLAS/LAPACK 例程来利用它们（例如，请参阅 the second part of this post）。 MATLAB 附带 Intel MKL 库，坦率地说，在 Intel 处理器上进行线性代数计算时，您无法击败它。
其他人已经在 File Exchange 上提到了一些将通用矩阵例程实现为 MEX 文件的提交（请参阅@natan 的答案）。如果您将它们与优化的 BLAS 库链接起来，它们会特别有效。