3

在 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 循环的情况下提取这些数据。谢谢你。

4

3 回答 3

5

建议的解决方案

我终于能够破解它以获得真正的矢量化解决方案，该解决方案用于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 110：列数保持在问题中使用的输入矩阵中行数的倍数。
  • 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。这是有道理的，因为根据定义，向量化操作意味着一次对许多元素执行操作，因此内存带宽是必不可少的。所以，如果你有一台更好的机器和更多的内存，这个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
于 2014-11-15T18:30:26.363 回答
3

你可以试试这个单线

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.

对于 30×1000 矩阵

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.
于 2013-06-12T16:21:33.277 回答
1

由于 sparse/matrixmult 方法有一个有趣的性能影响，我将发布一些额外的结果：

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.

随着矩阵维数的增长，普通循环的效率更高。

于 2013-06-12T17:10:38.027 回答