algorithm - 查找矩阵子矩阵的最有效方法 [matlab]

Question

假设我们有一个由 0 和 1 组成的矩阵

 0     1     1     1     0     0     0
 1     1     1     1     0     1     1
 0     0     1     0     0     1     0
 0     1     1     0     1     1     1
 0     0     0     0     0     0     1
 0     0     0     0     0     0     1

我们想找到所有具有这些属性的子矩阵(我们只需要角的行索引和列索引):

1. 至少包含L个1和L个0
2. 包含最多H个元素

即取上一个 L=1 和 H=5 的矩阵，子矩阵 1 2 1 4(行索引 1 2 和列索引 1 4)

 0     1     1     1
 1     1     1     1

满足属性1但是有8个元素(大于5)所以不好；

矩阵 4 5 1 2

 0     1
 0     0

很好，因为同时满足这两个属性。

目标是找到最小面积为 2*L、最大面积为 H 且至少包含 L 个 1 和 L 个零的所有子矩阵。

如果我们将矩阵视为矩形，通过查看从 H 到 2*L 的所有数字的约数，很容易找到最大面积 H 和最小面积 2*L 的所有可能的子矩形。

例如，当 H=5 且 L=1 时，所有可能的子矩形/子矩阵均由以下的除数给出

• H=5 -> 除数 [1 5] -> 区域 5 的可能矩形是 1x5 和 5x1
• 4 -> 除数 [1 2 4] -> 区域 4 的可能矩形是 1x4 4x1 和 2x2
• 3 -> 除数 [1 3] -> 区域 3 的可能矩形是 3x1 和 1x3
• 2*L=2 -> 除数 [1 2] -> 区域 2 的可能矩形是 2x1 和 1x2

我编写了这段代码，它为每个数字找到它的除数并循环它们以找到子矩阵。要找到子矩阵，它会这样做:以 1x5 子矩阵为例，代码所做的是固定矩阵的第一行并逐步(沿着矩阵的所有列)从子矩阵的左边缘移动子矩阵矩阵移动到矩阵的右边缘，然后代码固定矩阵的第二行并沿着所有列从左到右移动子矩阵，依此类推，直到到达最后一行。

它对所有 1x5 子矩阵执行此操作，然后考虑 5x1 子矩阵，然后是 1x4、4x1、2x2 等。

代码在 2 秒内完成工作(找到所有子矩阵)但是对于大矩阵，即 200x200，需要很多分钟才能找到所有子矩阵。所以我想知道是否有更有效的方法来完成这项工作，最终哪种方法最有效。

这是我的代码:

clc;clear all;close all

%% INPUT
P=  [0     1     1     1     0     0     0 ;
     1     1     1     1     0     1     1 ;
     0     0     1     0     0     1     0 ;
     0     1     1     0     1     1     1 ;
     0     0     0     0     0     0     1 ;
     0     0     0     0     0     0     1];
L=1; % a submatrix has to containg at least L ones and L zeros
H=5; % max area of a submatrix

[R,C]=size(P); % rows and columns of P
sub=zeros(1,6); % initializing the matrix containing the indexes of each submatrix (columns 1-4), their area (5) and the counter (6)
counter=1; % no. of submatrices found

%% FIND ALL RECTANGLES OF AREA >= 2*L & <= H
%
% idea: all rectangles of a certain area can be found using the area's divisors
%       e.g. divisors(6)=[1 2 3 6] -> rectangles: 1x6 6x1 2x3 and 3x2
tic
for sH = H:-1:2*L % find rectangles of area H, H-1, ..., 2*L
    div_sH=divisors(sH); % find all divisors of sH
    disp(['_______AREA ', num2str(sH), '_______'])

    for i = 1:round(length(div_sH)/2) % cycle over all couples of divisors        
        div_small=div_sH(i);
        div_big=div_sH(end-i+1);        

        if div_small <= R && div_big <= C % rectangle with long side <= C and short side <= R

            for j = 1:R-div_small+1 % cycle over all possible rows

                for k = 1:C-div_big+1 % cycle over all possible columns                    
                    no_of_ones=length(find(P(j:j-1+div_small,k:k-1+div_big))); % no. of ones in the current submatrix

                    if  no_of_ones >= L  &&  no_of_ones <= sH-L % if the submatrix contains at least L ones AND L zeros                        
                        %               row indexes    columns indexes         area         position
                        sub(counter,:)=[j,j-1+div_small , k,k-1+div_big , div_small*div_big , counter]; % save the submatrix
                        counter=counter+1;
                    end
                end
            end
            disp([' [', num2str(div_small), 'x', num2str(div_big), '] submatrices: ', num2str(size(sub,1))])
        end
        if div_small~=div_big % if the submatrix is a square, skip this part (otherwise there will be duplicates in sub)

            if div_small <= C && div_big <= R % rectangle with long side <= R and short side <= C

                for j = 1:C-div_small+1 % cycle over all possible columns

                    for k = 1:R-div_big+1 % cycle over all possible rows                        
                        no_of_ones=length(find(P(k:k-1+div_big,j:j-1+div_small)));

                        if no_of_ones >= L  &&  no_of_ones <= sH-L
                            sub(counter,:)=[k,k-1+div_big,j,j-1+div_small , div_big*div_small, counter];
                            counter=counter+1;
                        end
                    end
                end
                disp([' [', num2str(div_big), 'x', num2str(div_small), '] submatrices: ', num2str(size(sub,1))])
            end
        end
    end
end
fprintf('\ntime: %2.2fs\n\n',toc)

Answer 1

这是一个以 2D matrix convolution. 为中心的解决方案粗略的想法是将每个子矩阵形状的 P 与第二个矩阵进行卷积，使得结果矩阵的每个元素指示子矩阵中有多少元素，其左上角位于所述元素。像这样，您可以一次获得单个形状的所有解决方案，而无需遍历行/列，大大加快了速度(在我 8 岁的笔记本电脑上，200x200 矩阵需要不到一秒的时间)

P=  [0     1     1     1     0     0     0
    1     1     1     1     0     1     1
    0     0     1     0     0     1     0
    0     1     1     0     1     1     1
    0     0     0     0     0     0     1
    0     0     0     0     0     0     1];

L=1; % a submatrix has to containg at least L ones and L zeros
H=5; % max area of a submatrix
submats = [];
for sH = H:-1:2*L
    div_sH=divisors(sH); % find all divisors of sH
    for i = 1:length(div_sH) % cycle over all couples of divisors
        %number of rows of the current submatrix
        nrows=div_sH(i); 
        % number of columns of the current submatrix
        ncols=div_sH(end-i+1); 
        % perpare matrix to convolve P with
        m = zeros(nrows*2-1,ncols*2-1);
        m(1:nrows,1:ncols) = 1;
        % get the number of ones in the top left corner each submatrix
        submatsums = conv2(P,m,'same');
        % set values where the submatrices go outside P invalid
        validsums = zeros(size(P))-1;
        validsums(1:(end-nrows+1),1:(end-ncols+1)) = submatsums(1:(end-nrows+1),1:(end-ncols+1));
        % get the indexes where the number of ones and zeros is >= L
        topLeftIdx = find(validsums >= L & validsums<=sH-L);
        % save submatrixes in following format: [index, nrows, ncols]
        % You can ofc use something different, but it seemed the simplest way to me
        submats = [submats ; [topLeftIdx bsxfun(@times,[nrows ncols],ones(length(topLeftIdx),1))]];
    end
end