查找矩阵子矩阵的最有效方法 [matlab]
Most efficent way of finding submatrices of a matrix [matlab]
假设我们有一个由 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
我们想找到所有具有这些属性的子矩阵（我们只需要角的行索引和列索引）：
- 至少包含 L 个 1 和 L 个零
- 包含最大 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
很好，因为同时满足这两个属性。
然后objective就是求所有最小面积为2*L，最大面积为H且至少包含L个1和L个0的子矩阵
如果我们将矩阵视为矩形，通过查看从 H 到 2*L 的所有数字的约数，很容易找到最大面积 H 和最小面积 2*L 的所有可能的子矩形。
例如，对于 H=5 和 L=1，所有可能的 subrectangles/submatrices 都由
的除数给出
- 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)
首先，我建议您结合查找允许的子矩阵大小。
for smaller = 1:sqrt(H)
for larger = 2*L:H/smaller
# add smaller X larger and larger x smaller to your shapes list
接下来，从形状中的最小 个矩形开始。注意，任何一个小矩形的解都可以向任意方向扩展，达到H的面积限制，添加的元素不会使你找到的解失效。这将确定许多解决方案，而无需费心检查其中的人口。
跟踪您找到的解决方案。当您朝着更大的矩形前进时，您可以避免检查 solutions 集合中已有的任何内容。如果将其保存在散列 table 中，则检查成员资格是 O(1)。之后您需要检查的是 mostly-1 与 mostly-0 相邻的较大块。这应该会稍微加快处理速度。
这是否足以提供帮助？
这是一个以2D matrix convolution.为中心的解决方案粗略的想法是将每个子矩阵形状与第二个矩阵进行卷积P，使得结果矩阵的每个元素表示有多少个在其左上角位于所述元素的子矩阵。像这样，您可以一次获得单个形状的所有解决方案，而无需循环 rows/columns，大大加快了速度（在我 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
假设我们有一个由 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
我们想找到所有具有这些属性的子矩阵（我们只需要角的行索引和列索引）：
- 至少包含 L 个 1 和 L 个零
- 包含最大 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
很好，因为同时满足这两个属性。
然后objective就是求所有最小面积为2*L，最大面积为H且至少包含L个1和L个0的子矩阵
如果我们将矩阵视为矩形，通过查看从 H 到 2*L 的所有数字的约数，很容易找到最大面积 H 和最小面积 2*L 的所有可能的子矩形。
例如，对于 H=5 和 L=1，所有可能的 subrectangles/submatrices 都由
的除数给出- 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)
首先，我建议您结合查找允许的子矩阵大小。
for smaller = 1:sqrt(H)
for larger = 2*L:H/smaller
# add smaller X larger and larger x smaller to your shapes list
接下来，从形状中的最小 个矩形开始。注意，任何一个小矩形的解都可以向任意方向扩展，达到H的面积限制，添加的元素不会使你找到的解失效。这将确定许多解决方案，而无需费心检查其中的人口。
跟踪您找到的解决方案。当您朝着更大的矩形前进时，您可以避免检查 solutions 集合中已有的任何内容。如果将其保存在散列 table 中，则检查成员资格是 O(1)。之后您需要检查的是 mostly-1 与 mostly-0 相邻的较大块。这应该会稍微加快处理速度。
这是否足以提供帮助？
这是一个以2D matrix convolution.为中心的解决方案粗略的想法是将每个子矩阵形状与第二个矩阵进行卷积P，使得结果矩阵的每个元素表示有多少个在其左上角位于所述元素的子矩阵。像这样，您可以一次获得单个形状的所有解决方案，而无需循环 rows/columns，大大加快了速度（在我 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