M。我想要做的是以向量化的方式立即提取该矩阵的特定索引，并基于行索引对它们进行操作。sum(M(1, 1:1:100))。然后对于第二行，我想要sum(M(2, 1:2:100))。对于第三行，我想要sum(M(3, 1:3:100))等等。对于第十行，我当然有sum(M(10, 1:10:100))。我最终成功破解了它，实现了一种真正的向量化解决方案，使用逻辑索引从输入矩阵中选择要求和的元素。这种神奇的方法是通过bsxfun及其可选的函数句柄@mod实现的。代码如下：
[m,n] = size(M);
mask = bsxfun(@mod,1:n,(1:m)')==1; %//'# all of magic happens here as it creates
%// a logical mask of 1's at places in input matrix
%// whose elements are to be summed and 0's elsewhere.
mask(1,:) = 1; %// set the first row as all ones as we need to sum all of those
sumvals = sum(mask.*M,2); %// finally get the sum values
在这个基准测试部分中，我包含了四种方法 - 早期在本篇文章中列出的方法及其GPU端口版本，以及在其他解决方案中列出的基于arrayfun和sparse的方法。
基准测试使用了三组输入数据 -
Set 1：将输入矩阵的行数保持为列数的10倍数，与问题中所使用的输入矩阵相同。Set 2：扩展行数，使行数现在是列数的10x倍。这将真正测试循环代码，因为在这种情况下迭代次数会更多。Set 3：这个集合是set2的扩展，进一步增加了行数，因此将是真正向量化方法与其他方法之间的另一个重要测试。用于基准测试的函数代码如下所示 -
function sumvals = sumrows_stepped_bsxfun(M)
//... same as the code posted earlier
return
function sumvals = sumrows_stepped_bsxfun_gpu(M)
gM = gpuArray(M);
[m,n] = size(gM);
mask = bsxfun(@mod,gpuArray.colon(1,n),gpuArray.colon(1,m)')==1; %//'
sumvals = gather(sum(mask.*gM,2));
sumvals(1) = sum(M(1,:));
return
function S = sumrows_stepped_arrayfun(M)
[m,n] = size(M);
S = arrayfun(@(x) sum(M(x,1:x:n)), 1:m);
return
function B = sumrows_stepped_sparse(M)
sz = size(M);
A=sparse(sz(1),sz(2));
for n=1:sz(1),
A(n, 1:n:end)=1;
end
B=full(sum(M.*A,2));
return
timeit 用于计时 CPU based 代码，而 gputimeit 用于 GPU based 代码。MATLAB Version: 8.3.0.532 (R2014a)
Operating System: Windows 7
RAM: 3GB
CPU Model: Intel® Pentium® Processor E5400 (2M Cache, 2.70 GHz)
GPU Model: GTX 750Ti 2GB
结论：clear all; clc; close all
outputfile = 'results.xlsx';
delete(outputfile); %// remove file, so that new results could be written into
base_datasize_array = 40:60:400;
methods = {'BSXFUN on GPU','BSXFUN on CPU','ARRAYFUN','SPARSE'};
num_approaches = numel(methods);
num_sets = 3;
timeall_all = zeros(num_approaches,numel(base_datasize_array),num_sets);
datasize_lbs = cell(numel(base_datasize_array),num_sets);
for set_id = 1:num_sets
switch set_id
case 1
N1_arr = base_datasize_array*2;
N2_arr = N1_arr*10;
case 2
N2_arr = base_datasize_array*2;
N1_arr = N2_arr*10;
case 3
N2_arr = base_datasize_array;
N1_arr = N2_arr*40;
end
timeall = zeros(num_approaches,numel(N1_arr));
for iter = 1:numel(N1_arr)
M = rand(N1_arr(iter),N2_arr(iter));
f = @() sumrows_stepped_bsxfun_gpu(M);
timeall(1,iter) = gputimeit(f); clear f
f = @() sumrows_stepped_bsxfun(M);
timeall(2,iter) = timeit(f); clear f
f = @() sumrows_stepped_arrayfun(M);
timeall(3,iter) = timeit(f); clear f
f = @() sumrows_stepped_sparse(M);
timeall(4,iter) = timeit(f); clear f
end
timeall_all(:,:,set_id) = timeall;
wp = repmat({' '},numel(N1_arr),1);
datasize_lbs(:,set_id) = strcat(cellstr(num2str(N1_arr.')),' x ',...
wp,cellstr(num2str(N2_arr.')));
end
for set_id=1:num_sets
out_cellarr = cell(numel(methods)+1,numel(N1_arr)+1);
out_cellarr(1,1) = {'Methods'};
out_cellarr(2:end,1) = methods;
out_cellarr(1,2:end) = datasize_lbs(:,set_id);
out_cellarr(2:end,2:end) = cellfun(@(x) num2str(x),...
num2cell(timeall_all(:,:,set_id)),'Uni',0);
xlswrite(outputfile, out_cellarr,set_id);
end
S=arrayfun(@(n) sum(M(n,1:n:100)), 1:10)
另外，您可以提前创建一个稀疏矩阵。
A=sparse(100,10);
for n=1:10,
A(1:n:100, n)=1;
end
通过以下方式找到总和
S=diag(M*A);
这可以通过定义 A=sparse(10,100) 进一步优化，特别是对于更大的矩阵。
S=sum(M.*A,2);
我的快速基准测试
M=rand(10,100);
sz = size(M);
tic;
for k=1:10000,
for n=1:sz(1),
B(n)=sum(M(n,1:n:end));
end
end
toc
tic;
for k=1:10000,
B=arrayfun(@(n) sum(M(n,1:n:end)), 1:sz(1));
end
toc
tic;
for k=1:10000,
A=sparse(sz(2), sz(1));
for n=1:sz(1),
A(1:n:end, n)=1;
end
B=diag(M*A);
end
toc
tic;
A=sparse(sz(2),sz(1));
for n=1:sz(1),
A(1:n:end, n)=1;
end
for k=1:10000,
B=diag(M*A);
end
toc
tic;
A=sparse(sz(1),sz(2));
for n=1:sz(1),
A(n, 1:n:end)=1;
end
for k=1:10000,
B=sum(M.*A,2);
end
toc
返回值
Elapsed time is 0.552470 seconds.
Elapsed time is 2.409102 seconds.
Elapsed time is 0.638072 seconds.
Elapsed time is 0.052246 seconds.
Elapsed time is 0.061893 seconds.
针对30x1000的矩阵
Elapsed time is 1.785664 seconds.
Elapsed time is 3.954034 seconds.
Elapsed time is 4.760436 seconds.
Elapsed time is 0.926118 seconds.
Elapsed time is 0.865330 seconds.
对于一个1000乘以100的矩阵
Elapsed time is 51.389322 seconds.
Elapsed time is 63.443414 seconds.
Elapsed time is 68.327187 seconds.
Elapsed time is 29.056304 seconds.
Elapsed time is 1.147215 seconds.
由于稀疏/矩阵乘法方法具有有趣的性能效果，我将发布一些额外的结果：
M = rand(1000,100);
sz = size(M);
% PLAIN LOOP
tic
out1 = zeros(sz(1),1);
for k = 1:10000
for n = 1:sz(1)
out1(n) = sum(M(n,1:n:100));
end
end
toc
% SPARSE MATRIXMULT
tic
A = sparse(sz);
for n = 1:sz(1)
A(1:n:sz(2),n) = 1;
end
for k = 1:10000
out2 = diag(M*A);
end
toc
isequal(out1,out2) % ok
Plain loop: 11.441380 seconds.
Sparse/matrixmult: 27.503829 seconds.
M 30x1000，我得到了普通循环的运行时间为0.439726秒，而稀疏/矩阵乘法的运行时间为0.858779秒，其中创建 A 只需要0.004381秒（我有一个加速它的想法，但是这是无意义的）。 - Oleg