Most efficent way of finding combinations of submatrices of a matrix

Say we have a matrix P of zeros and ones and we want to find the best(*) combination of non-overlapping submatrices of P, so that each submatrix:

• contain at least L zeros and L ones (min area=2*L);
• contain at most H elements (max area=H).

(*) the best combination is the one which maximize the total number of elements in all the submatrices.

Notice that the best combination could be not unique, and that is not always possibile to cover all the elements of P with the submatrices of the best combination, i.e. with

L = 1
H = 6
P =
     0     0     0     0
     1     1     0     0
     0     0     0     0

one of the best combination is given by the two submatrices

     1     2     2     4   % top left corner 1 2 - bottom right corner 2 4
     1     3     1     1   % top left corner 1 1 - bottom right corner 3 1

which covers only 9 of the 12 elements of P.

The code I wrote to solve the problem is divided into 2 parts:

1. first it finds all the possibile submatrices of P which satisfy the two previous properties (at least L zeros and L ones, at maximum H elements) using the built-in function conv2 (more info in the code below); this part takes only few seconds and it is the easy one;
2. then it analyzes the set of all the submatrices using a recursive technique; this is the core of the problem and can take hours to find the solution (unless the computer freezes first).

The recursion is done in this way:

• call the recursive function using, as inputs, the set A of all the submatrices and a copy B of it
• fix the first submatrix of A and add it to the current combination
• find the set C of all the non-overlapping submatrices
• call the recursive function using, as inputs, A and C
• fix the first submatrix of C and add it to the current combination
• find the set D of all the non-overlapping submatrices
• call the recursive function using, as inputs, A and D
• and so on, until the set Z of non-overlapping submatrices is empty; the final combination is a vector containing all the submatrices, the counter (i.e. the number of submatrices in the combination) and the number of elements remained in the matrix after removing the elements contained in the submatrices of the combination
• then the second submatrix of Y is fixed, and the last submatrix of the previous combination is replaced with it
• if the set of non-overlapping submatrices is empty, the new final combination is compared with the previous one, if the set of non-overlapping submatrices is non empty then its first submatrix is added to the combination and so on
• the script is terminated when the perfect combination is found or when the first for cycle over all the submatrices of A is ended.

This method is fast only for small matrices (less than 10x10), and can be a nightmare for bigger matrices; I tested a 200x200 matrix and my computer freezed after few minutes. The problem is that if the set of all the submatrices contains thousands of elements, then the recursion will generate hundreds of nested for loops, consuming lot of RAM and CPU.

I wonder what is the most efficient method to achieve the goal, since my approach is very bad.

Here is my code:

%%PROBLEM
%
%  let P be a matrix whose elements are zeros and ones
%  find the best(*) combination of non-overlapping submatrices of P
%  so that each submatrix respect these properties:
%   - contains at least L zeros and L ones (min area=2*L)
%   - contains at most H elements (max area=H)
%
%  (*) the best is the one which maximize the total number of elements in all the submatrices
%
%  notices: the best combination could be not unique
%           is not always possibile to cover all the elements of P with the submatrices of the best combination
%
%%INPUT
P=round(rand(8,8)); L=1; H=5;
%P=dlmread('small.txt'); L=1; H=5;  % small can be found here https://pastebin.com/RTc5L8We
%P=dlmread('medium.txt'); L=2; H=8; % medium can be found here https://pastebin.com/qXJEiZTX
%P=dlmread('big.txt'); L=4; H=12;   % big can be found here https://pastebin.com/kBFFYg3K
%P=[0 0 0 0 0 1;0 0 0 0 0 1;0 1 0 1 0 1;0 0 0 0 0 0;0 0 0 0 0 0]; L=1; H=6;
P=[0 0 0 0 0;0 1 1 1 0;0 0 0 0 0]; L=1; H=6;
%P=[1,0,0,0,0;1,1,1,1,1;1,0,0,0,0]; L=1; H=5;
%
%%FIND ALL THE SUBMATRICES OF AREA >= 2*L & <= H
%
%  conv2(input_matrix,shape_matrix,'valid')
%  creates a matrix, where each element is the sum of all the elements contained in
%  the submatrix (contained in input_matrix and with the shape given by shape_matrix)
%  having its top left corner at said element
% 
%  ex.  conv2([0,1,2;3,4,5;6,7,8],ones(2,2),'valid')
%       ans =
%             8    12
%            20    24
%       where 8=0+1+3+4 12=1+2+4+5  20=3+4+6+7  24=4+5+7+8
%
s=[]; % will contain the indexes and the area of each submatrix
      % i.e.  1 3 2 5 9  is the submatrix with area 9 and corners in 1 2 and in 3 5 
for sH = H:-1:2*L
    div_sH = divisors(sH);
    fprintf('_________AREA %d_________\n',sH)
    for k = 1:length(div_sH)
        a = div_sH(k);
        b = div_sH(end-k+1);
        convP = conv2(P,ones(a,b),'valid');
        [i,j] = find((convP >= L) & (convP <= sH-L));
        if ~isempty([i,j])
            if size([i,j],1) ~= 1
%                        rows     columns           area
                s = [s;[i,i-1+a , j,j-1+b , a*b*ones(numel(i),1)]];
            else
                s = [s;[i',i'-1+a,j',j'-1+b,a*b*ones(numel(i),1)]];
            end
            fprintf('[%dx%d] submatrices: %d\n',a,b,size(s,1))
        end
    end
end
fprintf('\n')
s(:,6)=1:size(s,1);
%
%%FIND THE BEST COMBINATION
tic
[R,C]=size(P); % rows and columns of P
no_ones=sum(P(:)); % counts how many ones are in P
% a combination of submatrices cannot contain more than max_no_subm submatrices
if no_ones <= R*C-no_ones
    max_no_subm=floor(no_ones/L);
else
    max_no_subm=floor(R*C-no_ones/L);
end
comb(2,1)=R*C+1; % will contain the best combination
s_copy=s; % save a copy of s
[comb,perfect]=recursion(s,s_copy,comb,0,0,R,C,0,false,H,[],size(s,1),false,[0,0,0],0,0,0,0,0,0,max_no_subm);
fprintf('\ntime: %2.2fs\n\n',toc)
if perfect
    disp('***********************************')
    disp('***  PERFECT COMBINATION FOUND  ***')
    disp('***********************************')
end
%
%%PRINT RESULTS
if (R < 12 && C < 12)
    for i = 1:length(find(comb(2,3:end)))
        optimal_comb_slices(i,:)=s(comb(2,i+2),:);
    end
    optimal_comb_slices(:,1:5)
    P
end

and the function

function [comb,perfect,counter,area,v,hold_on,ijk,printed,first_for_i,second_for_i,third_for_i] = recursion(s,s_copy,comb,counter,area,R,C,k,hold_on,H,v,size_s,perfect,ijk,size_s_ovrlppd,size_s_ovrlppd2,printed,third_for_i,second_for_i,first_for_i,max_no_subm)
%
% OUTPUT (that is actually going to be used in the main script)
% comb [matrix] a matrix of two rows, the first one contains the current combination
%        the second row contains the best combination found
% perfect [boolean] says if the combination found is perfect (a combination is perfect if
%           the submatrices cover all the elements in P and if it is made up with
%           the minimum number of submatrices possible)
%
% OUTPUT (only needed in the function itself)
% counter [integer] int that keeps track of how many submatrices are in the current combination
% area [integer] area covered with all the submatrices of the current combination
% v [vector] keeps track of the for loops that are about to end
% hold_on [boolean] helps v to remove submatrices from the current combination
%
% OUTPUT (only needed to print results)
% ijk [vector] contains the indexes of the first three nested for loops (useful to see where the function is working)
% printed [boolean] used to print text on different lines
% first_for_i second_for_i third_for_i [integers] used by ijk
%
%
% INPUT
% s [matrix] the set of all the submatrices of P
% s_copy [matrix] the set of all the submatrices that don't overlap the ones in the current combination
%                 (is equal to s when the function is called for the first time)
% R,C [integers] rows and columns of P
% k [integer] area of the current submatrix
% H [integer] maximum number of cells that a submatrix can contains
% size_s [integer] number of rows of s i.e. number of submatrices in s
% size_s_ovrlppd [integer] used by ijk
% size_s_ovrlppd2 [integer] used by ijk
% max_no_subm [integer] maximum number of submatrices contained in a combination
%
%
%  the function starts considering the first submatrix (call it sub1) in the set 's' of all the submatrices
%  and adds it to the combination
%  then it finds 's_ovrlppd' i.e. the set of all the submatrices that don't overlap sub1
%  and the function calls itself considering the first submatrix (call it sub2) in the set 's_ovrlppd'
%  and adds it to the combination
%  then it finds the set of all the submatrices that don't overlap sub2 and
%  so on until there are no more non-overlapping submatrices
%  then it replaces the last submatrix in the combination with the second one of the last set of non-overlapping
%  submatrices and saves the combination which covers more elements in P
%  and so on for all the submatrices of the set 's'
%
%  DOWNSIDE OF THIS METHOD
%    if 's' contains thousands of submatrices, the function will create hundreds of nested for loops
%    so both time and space complexities can be really high and the computer might freeze
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%  SAVE AND RESET COMBINATIONS   %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   s_copy is empty when no more submatrices can be added to the current
%   combination, in this case we have to check if this combination is
%   better than the best combination previosly found, if so then we overwrite it
%
%   then we have to remove one or more submatrices from the combination (depending on
%   how many nested for loops are about to be closed)
%   and compute another combination
%   to 'remove one or more submatrices from the combination' it is necessary to do these things:
%    - reduce the area
%    - reduce the combination
%    - reduce the counter
%
    if isempty(s_copy)
        comb(1,2)=counter;  % final no of submatrices in the combination
        comb(1,1)=R*C-area; % no. of cells remained in P after removing the cells contained in the submatrices of the combination
%       if the combination just found is better than the previous overwrite it
        if comb(1,1)<comb(2,1) || (comb(1,1)==comb(2,1) && comb(1,2)<comb(2,2))
            comb(2,:)=comb(1,:);
            disp(['[area_left] ', num2str(comb(2,1)), ' [slices] ', num2str(comb(2,2))])
            printed=true;
        end
%        
        area=area-k; % tot area of the combination excluding the last sumatrix
        if ~isempty(v) && ~hold_on % more than one submatrix must be removed
            i=size(v,2);
            if i>length(find(v)) % if v contains at least one 0
                while v(i)==1 % find the index of the last 0
                    i=i-1;
                end
                last_i_counter=size(v(i+1:end),2); % no. of consecutive for loop that are about to end
                v=v(1:i-1);
            else
                last_i_counter=i;
                v=[];
            end
            for i=1:last_i_counter
                area=area-s(comb(1,2+counter-i),5); % reduce the area
            end
            comb(1,2+counter-last_i_counter:2+counter)=0; % remove submatrices from the combination
            counter=counter-(last_i_counter+1); % reduce the counter
            hold_on=true;
        else % exactly one submatrix must be removed
            comb(1,2+counter)=0;
            counter=counter-1;
        end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

FIND COMBINATIONS %%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    else
        for i = 1:size(s_copy,1) % fix the i-th submatrix
%
%           if the combination cover all elements of P & the no. of submatrices is the minimum possibile
            if comb(2,1)==0 && comb(2,2)==ceil(R*C/H)
                perfect=true;
                return
            end
%            
            pos=s_copy(i,6); % position of current submatrix in s
            comb(1,3+counter)=pos; % add the position to the combination
            k=s(pos,5); % save the area of the current submatrix
            area=area+k; % area covered with all the submatrices of the combination up now
            counter=counter+1;
%
%           if the area in P covered by the current combination could not
%           be bigger than the best combination found up now, then discard
%           the current combination and consider the next one
            if R*C-area-(max_no_subm-counter)*H > comb(2,1) && i < size(s_copy,1)
%              R*C-area-(max_no_subm-counter)*H can be negative if ceil(R*C/H) < max_no_subm
                counter=counter-1;
                comb(1,3+counter)=0;
                area=area-k;
            else
                s_ovrlppd=s_copy; % initializing the set of non-overlapping submatrices
                s_ovrlppd(s_copy(:,1)<=s_copy(i,2) & s_copy(:,2)>=s_copy(i,1) & s_copy(:,3)<=s_copy(i,4) & s_copy(:,4)>=s_copy(i,3),:)=[]; % delete submatrices that overlap the i-th one
                s_ovrlppd(s_ovrlppd(:,6)<pos,:)=[]; % remove submatrices that will generate combinations already studied
%               KEEP TRACK OF THE NESTED 'FOR' LOOPS ENDS REACHED
                if i==size(s_copy,1) % if i is the last cycle of the current for loop
                    v(size(v,2)+1)=1; % a 1 means that the code entered the last i of a 'for' loop
                    if size(s_ovrlppd,1)~=0 % hold on until an empty s_ovrlppd is found
                        hold_on=true;
                    else
                        hold_on=false;
                    end
                elseif ~isempty(v) && size(s_ovrlppd,1)~=0
                    v(size(v,2)+1)=0; % a 0 means that s_ovrlppd in the last i of a 'for' loop is not empty => a new 'for' loop is created
                end
%%%%%%%%%%%%%%%%%%%%%%%%

PRINT STATUS %%%

%%%%%%%%%%%%%%%%%%%%%%%%
                if size(s_copy,1)==size_s
                    ijk(1)=i;
                    ijk(2:3)=0;
                    fprintf('[%d,%d,%d]\n',ijk)
                    size_s_ovrlppd=size(s_ovrlppd,1);
                    first_for_i=i;
                    second_for_i=0;
                elseif size(s_copy,1)==size_s_ovrlppd
                    ijk(2)=i;
                    ijk(3)=0;
                    if ~printed
                        fprintf(repmat('\b',1,numel(num2str(first_for_i))+numel(num2str(second_for_i))+numel(num2str(third_for_i))+2+2+1)) % [] ,, return
                    else
                        printed=false;
                    end
                    fprintf('[%d,%d,%d]\n',ijk)
                    size_s_ovrlppd2=size(s_ovrlppd,1);
                    second_for_i=i;
                    third_for_i=0;
                elseif size(s_copy,1)==size_s_ovrlppd2
                    ijk(3)=i;
                    if ~printed
                        fprintf(repmat('\b',1,numel(num2str(first_for_i))+numel(num2str(second_for_i))+numel(num2str(third_for_i))+2+2+1))
                    else
                        printed=false;
                    end
                    fprintf('[%d,%d,%d]\n',ijk)
                    third_for_i=i;
                end
                [comb,perfect,counter,area,v,hold_on,ijk,printed,first_for_i,second_for_i,third_for_i]=recursion(s,s_ovrlppd,comb,counter,area,R,C,k,hold_on,H,v,size_s,perfect,ijk,size_s_ovrlppd,size_s_ovrlppd2,printed,third_for_i,second_for_i,first_for_i,max_no_subm);
            end
        end
    end
end