【问题标题】:Finding sums of array elements that exceed a given value?查找超过给定值的数组元素的总和?
【发布时间】:2017-03-21 16:35:05
【问题描述】:

为此,我正在尝试利用 MATLAB 中的矢量化,但我可能不得不求助于 for 循环。我真的不想那样做!是时候学习算法了。

给定这个(11×3)数组:

x = [...
4.9000   -0.1000   -5.1000
4.6000   -0.4000   -5.4000
3.0000   -2.0000   -7.0000
2.9000   -2.1000   -7.1000
2.9000   -2.1000   -7.1000
2.9000   -2.1000   -7.1000
2.8000   -2.2000   -7.2000
2.7000   -2.3000   -7.3000
2.7000   -2.3000   -7.3000
2.2000   -2.8000   -7.8000
1.8000   -3.2000   -8.2000
];

我想找到数组中 11 个元素的所有 3^11 = 177147 个可能的总和,其中 11 个元素中的每一个都来自不同的行。然后,我想将超过阈值 16.0 的总和以及构成这些总和的 11 个元素存储在一个(12 乘?)数组中。

有什么让我开始的想法吗?感谢您的帮助。

【问题讨论】:

    标签: matlab matrix subset vectorization combinatorics


    【解决方案1】:

    以下是如何以矢量化方式进行操作:

    TR = 16;
    
    sets = num2cell(single(x),2);
    
    c = cell(1, numel(sets));
    [c{:}] = ndgrid( sets{:} );
    cartProd = cell2mat( cellfun(@(v)v(:), c, 'UniformOutput',false) );
    
    validRows = cartProd(sum(cartProd,2) > TR,:); % output is [353x11]
    

    请注意我如何使用single 来节省空间并使计算速度稍快。

    上述解决方案是this答案的改编。


    经过进一步的思考,我想我想出了一种既更快又更节省内存的方法。我们通过索引 x 来做到这一点,然后在索引上执行前面的过程。你可能会问,为什么这样更好?这是因为我们可以将索引存储为uint8,这比double 甚至single 消耗的内存要少得多。我们还可以保持x 的完整double 精度。因此:

    function validRows = q42933114(x,thresh)
    %% Input handling
    if nargin < 2
      thresh = 16;
    end
    if nargin < 1
      x = [...
        4.9000   -0.1000   -5.1000
        4.6000   -0.4000   -5.4000
        3.0000   -2.0000   -7.0000
        2.9000   -2.1000   -7.1000
        2.9000   -2.1000   -7.1000
        2.9000   -2.1000   -7.1000
        2.8000   -2.2000   -7.2000
        2.7000   -2.3000   -7.3000
        2.7000   -2.3000   -7.3000
        2.2000   -2.8000   -7.8000
        1.8000   -3.2000   -8.2000
      ];
    end
    
    I = reshape(uint8(1:numel(x)),size(x));
    
    sets = num2cell(I,2);
    
    c = cell(1, numel(sets));
    [c{:}] = ndgrid( sets{:} );
    cartProd = cell2mat( cellfun(@(v)v(:), c, 'UniformOutput',false) );
    validRows = x(cartProd(sum(x(cartProd),2) > thresh,:));
    

    内存消耗对比:

    方法1(旧):

    >> whos
      Name                Size              Bytes  Class     Attributes
    
      c                   1x11            7795700  cell                
      cartProd       177147x11            7794468  single              
      sets               11x1                1364  cell                
      validRows         353x11              15532  single              
    

    方法 2(新):

    >> whos
      Name                Size              Bytes  Class     Attributes
    
      c                   1x11            1949849  cell                
      cartProd       177147x11            1948617  uint8               
      sets               11x1                1265  cell                
      validRows         353x11              31064  double              
    

    我们看到内存消耗确实更小(大约 4 倍),正如预期的那样。

    运行时对比:

    Method 1 -- 0.0110
    Method 2 -- 0.0186
    

    这里我们看到 2nd 方法实际上有点慢。在对此进行分析时,我们看到原因是x(...),它相对昂贵。

    【讨论】:

    • 啊,谢谢!我也在研究一个解决方案,我刚刚得到它,尽管效率要低得多。我最终使用了ndgridspalloc
    • @abscissa 您可能希望将您的解决方案发布为另一个答案。
    【解决方案2】:

    我是这样做的。变量名显然还有改进的余地。

    注意有353 匹配的行,这与@Dev-iL 的答案一致。

    p = 11;
    [a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11] = ...
        ndgrid(x(1,:),x(2,:),x(3,:),x(4,:),x(5,:),x(6,:),x(7,:),x(8,:),x(9,:),x(10,:),x(11,:));
    a = a1+a2+a3+a4+a5+a6+a7+a8+a9+a10+a11;
    y = spalloc(p+1,3^p,(p+1)*3^p);
    for i = 1:3^p
        if a(i) >= 16.1
            y(:,i) = [a1(i),a2(i),a3(i),a4(i),a5(i),a6(i),a7(i),a8(i),a9(i),a10(i),a11(i),a(i)];
        end
    end
    nnz(y(p+1,:)); % 353 rows matching the criteria
    

    【讨论】:

      【解决方案3】:

      我认为没有比使用 for 循环更好的运气了。可能有一个 Matlab 函数用于生成所有 3^11 组合,并将其用作一种索引,但那样会消耗大量内存。

      代码也很难阅读。

      然而,最新版本的 Matlab 并没有表现得那么糟糕,因为它们 JIT 代码。在它刚刚被解释之前,或者 JIT-ing 被用于特定目的之前。因此,您不会想在 Matlab 中重新实现矩阵例程,但是对于像这样的简单代码,它应该可以很好地执行。

      【讨论】:

      • 这应该是一条评论。
      猜你喜欢
      • 2022-01-09
      • 2020-02-22
      • 1970-01-01
      • 2018-11-03
      • 1970-01-01
      • 2014-05-26
      • 2021-05-25
      • 1970-01-01
      • 1970-01-01
      相关资源
      最近更新 更多