针对多维矩阵数据优化二维区间划分

Question

有一个四维矩阵数据，以及每个元素对应的X坐标data_X和Y坐标data_Y。我想根据edge_X和edge_Y进行二维区间分割，或者说网格划分，得到一个新的数组data_class。当前的循环执行速度非常慢。如何提高计算速度？

clc;clear
rng("default")
data = randi([0, 100], 40,30,20, 10);
data_X =  randi([0, 100], 40,30,20, 10);
data_Y =  randi([0, 100], 40,30,20, 10);

% Define the edges for histograms in X and Y directions
edge_X = [0:10:100];
edge_Y = [0:20:100];

sz = size(data);
for ix = 1:sz(1)
    for iy = 1:sz(2)
        for ii = 1:sz(3)
            for ij = 1:sz(4)
               % Extract the current data, X, and Y values
                g_data = data(ix, iy, ii, ij);
                g_X = data_X(ix, iy, ii, ij);
                g_Y = data_Y(ix, iy, ii, ij);
    
                % Compute the 2D histogram count  
                g_id = histcounts2(g_X, g_Y, edge_X, edge_Y);
                g_new = g_id * g_data;
                data_class  (ix,iy,ii,ij,:,:) = g_new; 
            end
        end
    end
end

由于data的每个元素都是单独处理的，所以g_data在一个bin中只会被记为1，然后

g_new = g_id * g_data; 

g_new是重分类后g_data的位置

当前循环执行速度非常慢。如何提高计算速度？

Answer 1

您应该做的第一件事是查看编辑器中关于

data_class

随着每次循环迭代而增长的橙色下划线警告。在所有循环之前初始化它可以立即节省大部分运行时间

在循环之前预分配

data_class

：

szout = [sz, [numel(edge_X),numel(edge_Y)]-1];
data_class = NaN(szout);

但主要瓶颈是对每个值单独调用

histcounts2

，您可以通过将其置于

ii/ij

循环之外来节省大量时间，而不是使用第一个输入作为其所在容器的

1

计数在中，使用第四个和第五个输出来存储每个值应进入的容器。在您的情况下，这是行索引和列索引。

以下是一些时间安排：

• 基线代码：20.53秒
• 预分配

  data_class

  ：7.75秒（-62%总运行时间）
• 解除嵌套

  histcounts2

  ：0.53秒（-97%总运行时间）
• 仅分配非零值：0.10秒（-99.5％总运行时间）

解巢

histcounts2

sz = size(data);
szout = [sz, [numel(edge_X),numel(edge_Y)]-1];
data_class = NaN(szout);
g_new = zeros(szout(end-1:end));
for ix = 1:sz(1)
    for iy = 1:sz(2)

        g_X = data_X(ix, iy, :, :);
        g_Y = data_Y(ix, iy, :, :);
        [~,~,~,bx,by] = histcounts2(g_X, g_Y, edge_X, edge_Y);
        bx = squeeze(bx);
        by = squeeze(by);
        g_data = squeeze( data(ix, iy, :, :) );
            
        for ii = 1:sz(3)
            for ij = 1:sz(4)
                g_new(:) = 0;        
                g_new(bx(ii,ij),by(ii,ij)) = g_data(ii,ij);
                data_class(ix,iy,ii,ij,:,:) = g_new; 
            end
        end
    end
end

此外，您不需要仅为每个切片初始化零数组

g_new

来设置一个值。相反，将

data_class

初始化为零矩阵，然后覆盖您感兴趣的特定值：

仅分配非零值

sz = size(data);
szout = [sz, [numel(edge_X),numel(edge_Y)]-1];
data_class = zeros(szout);
for ix = 1:sz(1)
    for iy = 1:sz(2)

        g_X = data_X(ix, iy, :, :);
        g_Y = data_Y(ix, iy, :, :);
        [~,~,~,bx,by] = histcounts2(g_X, g_Y, edge_X, edge_Y);
        bx = squeeze(bx);
        by = squeeze(by);
        g_data = squeeze( data(ix, iy, :, :) );
            
        for ii = 1:sz(3)
            for ij = 1:sz(4)
                data_class(ix,iy,ii,ij,bx(ii,ij),by(ii,ij)) = g_data(ii,ij); 
            end
        end
    end
end

您可以使用

sub2ind

做一些事情来完全摆脱

ii/ij

循环，但在快速测试中，与仅使用循环相比，没有显示出任何性能。