从MATLAB矩阵中提取数据而不进行for循环

Question

在MATLAB中,假设我有一个10 x 100矩阵,称为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)).

我在for循环中有这个,但我想看看是否有一种方法可以在没有for循环的情况下提取这些数据.谢谢.

Answer 1

提出的解决方案

我终于破解了它,用于真正的矢量化解决方案,用于logical indexing从输入矩阵中选择要求的元素.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

由此获得的基准结果 -

结论

1. 对于行数小于列数的数据,迭代次数是一个小数,循环代码似乎占上风.

2. 随着我们增加行数,真正的矢量化方法的好处变得清晰.您还会注意到bsxfun on CPU基于方法对第3组做得很好,直到12000 x 300对非非向量化方法进行标记,其背后的原因是,bsxfun创建了这个巨大的逻辑掩码,此时内存带宽需求变得太高而无法应对计算能力的bsxfun.这是有道理的,因为定义的矢量化操作意味着一次性对许多元素执行操作,因此内存带宽是必不可少的.因此,如果你有一台更好的RAM机器,那么这个12000 x 300标记应该进一步延伸.

3. 如果可以进一步扩展行数,只要保持内存带宽,矢量化解决方案的好处就会变得更加清晰.

基准代码

最后这里是基准测试代码,如果有人想在他们的系统上测试它 -

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