# Finding the largest rectangular submatrix

I have a sparse non-symmetric binary matrix with a block structure. The dimensions of the matrix are thousands rows and less than one hundred columns.

How do I identify the largest non-contiguous rectangular sub-matrix consisting only of 1-entries?

### Edit

The pattern is either very sparse, or there are very few large submatrices. However, the blocks are likely to overlap.

In the simplified example below, it would be elements {2, 1} towards {2,6} and {5,1} towards {6,6}.

$$\qquad mat=\left( \begin{array}{ccc} 1&1&1&1&0&0\\ 1&1&1&1&1&1\\ 0&0&0&1&1&1\\ 1&1&0&1&1&1\\ 1&1&1&1&1&1\\ 1&1&1&1&1&1\\ \end{array} \right)$$

• How is this related to graphs and networks? Might there be an additional structure behind the matrix? For example, if mat is an adjacency matrix of a graph, you might be looking for a maximal complete subgraph or a clique. Szabolcs' package "IGraphM" has tools for that... – Henrik Schumacher Jan 24 '19 at 9:14
• You still want to solve the problem for general binary matrices? (It is not unlikely that his discrete optimization problem is a very hard...) – Henrik Schumacher Jan 24 '19 at 9:22
• If you want non-contiguous too, then this is indeed the clique problem, which is NP-complete, and there's not going to be a simpler solution. Use FindClique, then filter for blocks that also have 1s on the diagonal. – Szabolcs Jan 24 '19 at 9:39
• Before FindClique, remove every row/column that has a 0 on the diagonal. – Szabolcs Jan 24 '19 at 9:40
• I added the size (thousands long, less than a hundred wide and that I am looking for non-contiguous solutions. – Sander Jan 24 '19 at 9:54

## Introduction

The problem you describe is called the maximum biclique problem in graph theory.

Definitions: A clique is a complete subgraph. A biclique is a complete bipartite subgraph (of a bipartite graph).

We can interpret your matrix $$A$$ as a bipartite incidence matrix: $$a_{ij}=1$$ means that vertex $$i$$ of the first partition is connected to vertex $$j$$ of the second partition of a bipartite graph.

A question remains: what does "largest" submatrix, or equivalently "largest" biclique mean? There are multiple interpretations, e.g.

1. Submatrix with most elements = biclique with most edges. The maximum edge biclique problem is NP-complete, meaning that there is no fast exact solution for large cases.

2. Submatrix with most #rows + #columns in total = biclique with most vertices. This problem can be solved in polynomial time by reducing it to bipartite matching.

## Submatrix with most elements (maximum edge biclique)

I believe you want the 1st interpretation, but this does not seem to have fast solutions. Other answers already give slow solutions.

## Submatrix with max #rows + #columns (maximum vertex biclique)

We can attack the 2nd problem in multiple ways. The easiest is to reduce it to simple clique finding: connect all vertex pairs within each of the two partitions, and look for maximum cliques. Clique finding is still computationally hard, but Mathematica already has a function for it, so the implementation will be easy.

### Map to simple clique finding

The problem you are trying to solve is equivalent to biclique finding, i.e. finding complete bipartite graphs. We can think of the matrix as the bipartite incidence matrix of a bipartite graph.

Let us take this bipartite graph, and connect all vertex pair within each of the two partitions. Then we just need to look for maximal cliques to find maximal bicliques.

The following code looks for all maximal cliques:

Clear[getCliques]
getCliques[mat_] :=
Module[{am, g, m, n},
{m, n} = Dimensions[mat];
am = ArrayFlatten[{ (* transform bipartite incidence matrix to ajdacency matrix *)
{ConstantArray[1, {m, m}], mat},
{Transpose[mat], ConstantArray[1, {n, n}]}
}
];
g = AdjacencyGraph[am];
FindClique[g, Length /@ FindClique[g], All]
]


Here's some code to convert these to row/column indices and highlight them in the graph:

Clear[cliqueToRowCol]
cliqueToRowCol[mat_][clique_] :=
With[{m = First@Dimensions[mat]}, {Select[clique, # <= m &],
Select[clique, # > m &] - m}]

Clear[highlight]
highlight[mat_][rowCol_] :=
MatrixForm@MapAt[Style[#, Red] &, mat, Tuples[rowCol]]


This is your matrix:

mat = {{1, 1, 1, 1, 0, 0}, {1, 1, 1, 1, 1, 1}, {0, 0, 0, 1, 1, 1}, {1, 1, 0, 1, 1, 1}, {1, 1, 1, 1, 1, 1}, {1, 1, 1, 1, 1, 1}};


And this is the solutions we can find:

cliqueToRowCol[mat] /@ getCliques[mat]
(* { {{2, 5, 6},    {1, 2, 3, 4, 5, 6}},
{{2, 4, 5, 6}, {1, 2, 4, 5, 6}} } *)

highlight[mat] /@ cliqueToRowCol[mat] /@ getCliques[mat]


Once again, this method does not find submatrices which have the most elements. It finds those for which #row + #columns is maximal. In this case, one of the solutions happens to also have the most element ($$4 \times 5 = 20$$). The other only has $$3 \times 6 = 18$$. For both sub-matrices, $$4+5 = 3+6 = 9$$.

### Map to maximum matching

The other way is as follows:

Maximum vertex biclique finding is equivalent to finding a maximum independent vertex set in the bipartite complement graph. This, in turn, is equivalent to finding a minimum vertex cover: the vertices not in the covert will form the independent vertex set. According to Kőnig's theorem, in bipartite graphs, the minimum vertex cover can be formed from a maximum matching.

There is a function in Mathematica for finding a maximum matching: FindIndependentEdgeSet. To transform it to a min vertex cover, we can use a relatively simple algorithm, which so far I was too lazy to implement (might update this answer: http://tryalgo.org/en/matching/2016/08/05/konig/

• I am not an expert but would like to understand why, based on your statement: "Maximum vertex biclique finding is equivalent to finding a maximum independent vertex set in the bipartite complement graph." I could not use FindIndependentVertexSet on the complement graph? – Sander Jan 29 '19 at 4:25
• @Sander That would be basically the same as the method I already provided in the previous section (find cliques). The point of finding an independent edge set (matching) is that this can be solved in polynomial time, and is equivalent (through a series of transformations) to out original problem in a bipartite graph only. – Szabolcs Jan 29 '19 at 8:26
• @Sander But you are right that the implementation would have been a bit shorter with FindIndependentVertexSet – Szabolcs Jan 29 '19 at 8:27
• Got it ... Thanks for support and background; still digesting your answer, a lot of knowledge there... – Sander Jan 29 '19 at 9:07
• @Sander Actually, I learn as I go ... your question was an opportunity to learn some new things. – Szabolcs Jan 29 '19 at 9:17

Update: Putting together pieces from several sources (links below) to identify the largest contiguous rectangle in a binary matrix:

ClearAll[poP, stutteringAccumulate, largestRectangleInHistogram, maxRectangle]
SetAttributes[poP, HoldAllComplete];
poP[a_] := Module[{b}, If[EmptyQ[a], False, b = Last[a]; Set[a, Most[a]]; b]]
stutteringAccumulate = FoldList[#2  #1 + #2 &, #] &;
largestRectangleInHistogram = Module[{max = 0, a = Join[{-1}, #, { -1}], n = 2 + Length@#,
stack = {1},  h, area, i, index = 1, height = 0},
For[i = 1, i <= n, ++i,
While[a[[i]] < a[[Last@stack]],
h = a[[poP[stack]]];
area = h (i - Last[stack] - 1); max = Max[max, area];
If[max > area, index = index; height = height, index = i; height = h];
]; AppendTo[stack, i]];
{height, {# - 1 - max/height, # - 2} &@index, max}] &;
maxRectangle[mat_] := Module[{lr = largestRectangleInHistogram /@ stutteringAccumulate[mat],
l = List /@ Range[Length@mat]},
{#4 - {# - 1, 0}, #2, #3} & @@ MaximalBy[Join[lr, l, 2], Last][[1]]];


Examples:

mat = {{1, 1, 1, 1, 0, 0}, {1, 1, 1, 1, 1, 1}, {0, 0, 0, 1, 1, 1}, {1, 1, 0, 1, 1, 1},
{1, 1, 1, 1, 1, 1}, {1, 1, 1, 1, 1, 1}};


Construct a matrix of histograms:

histograms = stutteringAccumulate @ mat


{{1, 1, 1, 1, 0, 0}, {2, 2, 2, 2, 1, 1}, {0, 0, 0, 3, 2, 2}, {1, 1, 0, 4, 3, 3}, {2, 2, 1, 5, 4, 4}, {3, 3, 2, 6, 5, 5}}

Find the largest rectangle for each row of histograms:

largestrecs = largestRectangleInHistogram /@ histograms


{{1, {1, 4}, 4}, {2, {1, 4}, 8}, {2, {4, 6}, 6}, {3, {4, 6}, 9}, {4, {4, 6}, 12}, {5, {4, 6}, 15}}

Pick from largestrecs the one with largest area:

{rows, cols, area} = maxRectangle[mat]


{{2, 6}, {4, 6}, 15}

Row[Labeled[##, Top] & @@@ Transpose[{MatrixForm /@ {mat, histograms,
MapAt[Style[#, Red, Bold] &, mat, Span @@@ maxRectangle[mat][[;; 2]]]},
{"mat", "histograms", "max rectangle"}}]]


Row[Labeled[BarChart[#, ImageSize -> 100,
Background -> If[maxRectangle[mat][[-1]] ==
largestRectangleInHistogram[#][[-1]], LightBlue, White]],
Style[largestRectangleInHistogram@#, 12], Top] & /@ histograms, Spacer[5]]


With SeedRandom[1]; mat = RandomInteger[1, {20, 40}]; as input

maxRectangle[mat]


{{17, 20}, {32, 33}, 8}

MatrixForm @ MapAt[Style[#, Red, Bold] &, mat, Span @@@ maxRectangle[mat][[;; 2]]]


Grid[Partition[Labeled[BarChart[#, ImageSize -> 100,
Background -> If[maxRectangle[mat][[-1]] == largestRectangleInHistogram[#][[-1]],
LightBlue, White]], Style[largestRectangleInHistogram@#, 10], Top] & /@
(stutteringAccumulate@mat), 10]]


SeedRandom[123]
dta = RandomInteger[10, 50];
lr = largestRectangleInHistogram[dta];
BarChart[dta, BarSpacing -> 0,
ChartLabels -> Placed[Range@Length@dta, Axis], ImageSize -> Large,
PlotLabel -> Style[lr, 16],
Epilog -> {EdgeForm[Red], FaceForm[Opacity[.5, Red]],
Rectangle @@ (Transpose[{lr[[2]], {0, lr[[1]]}}] + {{-1/2, 0}, {1/2, 0}}) }]


Sources:

The idea of using an increasing stack to find the largest rectangle in a histogram and implementation is from this answer by Pei. The function largestRectangleInHistogram above is a Mathematica implementation of Pei's python function largestRectangleArea which is modified to return the column indices and the height in addition to the area of the largest rectangle.

The function poP is a slightly modified version of Pop from rosettacode - Stack.

The function stutteringAccumulate is from the posts by ciao and by Chip Hurst.

Okkes's links to Tushar Roy's YouTube videos has been extremely useful; especially, Maximum Rectangular Area in Histogram and Maximum Size Rectangle of All 1's Dynamic Programming.

Update 2: Dealing with non-necessarily-contiguous case for small matrices:

sa = SparseArray[mat];
al = DeleteCases[sa["AdjacencyLists"], {}];
nzprows = Union@sa["NonzeroPositions"][[All, 1]];
rowindices = MaximalBy[Subsets[nzprows, {2, Infinity}],
Length[#] Length[Intersection @@ #] &@al[[#]] &, 10];
rowscols = {#, Intersection @@ al[[#]]} & /@ rowindices;

Grid[Prepend[{## & @@ #, Times @@ Length /@ #} & /@ rowscols,
{"rows", "columns", "area"}], Dividers -> All] // TeXForm


$$\begin{array}{|c|c|c|} \hline \text{rows} & \text{columns} & \text{area} \\ \hline \{2,4,5,6\} & \{1,2,4,5,6\} & 20 \\ \hline \{2,5,6\} & \{1,2,3,4,5,6\} & 18 \\ \hline \{1,2,5,6\} & \{1,2,3,4\} & 16 \\ \hline \{2,4,5\} & \{1,2,4,5,6\} & 15 \\ \hline \{2,4,6\} & \{1,2,4,5,6\} & 15 \\ \hline \{4,5,6\} & \{1,2,4,5,6\} & 15 \\ \hline \{1,2,4,5,6\} & \{1,2,4\} & 15 \\ \hline \{2,3,4,5,6\} & \{4,5,6\} & 15 \\ \hline \{2,5\} & \{1,2,3,4,5,6\} & 12 \\ \hline \{2,6\} & \{1,2,3,4,5,6\} & 12 \\ \hline \end{array}$$

Original answer:

A brute force approach:

mat = {{1, 1, 1, 1, 0, 0}, {1, 1, 1, 1, 1, 1}, {0, 0, 0, 1, 1, 1}, {1, 1, 0, 1, 1, 1},
{1, 1, 1, 1, 1, 1}, {1, 1, 1, 1, 1, 1}};
pairs = Transpose /@ MaximalBy[DeleteDuplicates[CoordinateBounds /@
Subsets[SparseArray[mat]["NonzeroPositions"], {2}]],
Min[#] Total[#, 2] &@mat[[## & @@ Span @@@ #]] &]


{{{2, 4}, {6, 6}}}

• A brute brain force approach as well, thanks. Will try to understand this ... do you think the Subsets expansion will become a constraint? – Sander Jan 24 '19 at 12:29
• @Sanders, although we are considering only 2-subsets of the non-zero positions Subsets is the brute part likely to cause pain for large and relatively dense input matrices. – kglr Jan 24 '19 at 12:43
• Apologies, there is a larger selection possible in my answer that I omitted. At least you helped me a lot by giving me a fun challenge to see if I can work with your answer as a basis. – Sander Jan 24 '19 at 23:28
• it's very insight- and usefull nonetheless. Thanks! – Sander Jan 27 '19 at 3:39

This finds ONLY SQUARE sub matrix. Based on this explanation.

    SeedRandom@8
{m, n} = {15, 20};
(initMat = RandomChoice[{20, 1} -> {1, 0}, {m, n}]) // MatrixForm;

mat = initMat;
MatrixForm[
mat2 = Table[
If[mat[[i + 1, j + 1]] == 0, mat[[i + 1, j + 1]] = 0,
mat[[i + 1, j + 1]] =
Min[{mat[[i, j]], mat[[i, j + 1]], mat[[i + 1, j]]}] + 1], {i,
m - 1}, {j, n - 1}]];

max = Max@mat2;

pos = Position[mat2, max] + 1;

pairs = Table[{pos[[i]] - max + 1, pos[[i]]}, {i, Length@pos}];

highlight[list_, position_] :=
Grid[list, Background -> {None, None, # -> Yellow & /@ position}]

Table[highlight[initMat,
Join @@ CoordinateBoundsArray[Transpose@pairs[[i]]]], {i,
Length@pairs}]


• There is an explanation for rectangle matrix which I don't quite understand youtube.com/watch?v=g8bSdXCG-lA – OkkesDulgerci Jan 24 '19 at 20:07
• Thanks, unfortunately .... it’s mostly rectangular squares I am dealing with here. – Sander Jan 25 '19 at 12:01

Here's a brute force for the non-contiguous rectangular submatrix. Use your matrix.

mat = {{1, 1, 1, 1, 0, 0},
{1, 1, 1, 1, 1, 1},
{0, 0, 0, 1, 1, 1},
{1, 1, 0, 1, 1, 1},
{1, 1, 1, 1, 1, 1},
{1, 1, 1, 1, 1, 1}};

numRows = mat//Length;
numCols = mat//First//Length;


Get all permutations of rows and cols, skipping the first which is just { }

rows = Subsets[Range[numRows]] // Rest
cols = Subsets[Range[numCols]] // Rest


Do an Outer of the possible combinations, capturing the sum of their elements but noting the ones that have a 0 in them.

res = Outer[Total@Total@mat[[#1, #2]]*Min[mat[[#1, #2]]] &, rows, cols, 1];
bestSub = Max[res]


20

Find where it is at, noting that there may be multiple solutions. Just take the first.

location = Position[res,bestSub]//First;


Go back and get the particular rows and columns.

bestRows = rows[[location//First]]


{2, 4, 5, 6}

bestCols = cols[[location//Last]]


{1, 2, 4, 5, 6}

mat[[bestRows,bestCols]]


$$\left( \begin{array}{ccccc} 1 & 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 & 1 \\ \end{array} \right)$$

There's an obvious dynamic programming approach, but this at least gets the ball rolling.

EDIT

Its's fairly easy to bound the area of the rectangular submatrix, in order to restrict the search quite a bit. Create a random matrix 8 rows by 5 columns

numRows = 8;
numCols = 5;
matran = RandomChoice[{0, 1}, {numRows, numCols}]


$$\left( \begin{array}{ccccc} 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ 1 & 1 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 & 1 \\ 1 & 0 & 0 & 1 & 1 \\ 0 & 0 & 0 & 0 & 1 \\ 1 & 1 & 1 & 1 & 1 \\ 0 & 1 & 1 & 1 & 1 \\ \end{array} \right)$$

It is helpful to sort each row large to small, and then the matrix large to small, to see how we can bound the problem.

(ms = -1 Sort[Sort /@ (-matran)]) // MatrixForm


$$\left( \begin{array}{ccccc} 1 & 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 & 0 \\ 1 & 1 & 1 & 0 & 0 \\ 1 & 1 & 1 & 0 & 0 \\ 1 & 1 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 \\ \end{array} \right)$$

For the upper bound, count the number of 1's in each row, and then sort on that, large to small.

rowSums = Total /@ matran // Sort // Reverse
(* {5, 4, 3, 3, 2, 1, 1, 1} *)


Look at the list of max possible values, and the max of that.

maxPossibleList = MapIndexed[#1*#2 &, rowSums] // Flatten
(* {5, 8, 9, 12, 10, 6, 7, 8} *)

upperBound = Max[maxPossibleList]
(* 12 *)


For the lower bound, it's a little quirky to calculate. Idea is that values in rows are ordered in the worst possible ordering.

accRowSums = Accumulate[rowSums];
accNumCols = numCols*(Range[numRows] - 1);
mins = accRowSums - accNumCols ;
minList = MapIndexed[Max[#1*#2,0] &, mins] // Flatten
(* {5, 8, 6, 0, 0, 0, 0, 0} *)


The zero values mean in the worst case, it is possible to order the values in the rows so that there is no overlap of all terms. The lower bound is

lowerBound = Max[minPossibleList]
(* 8 *)


So you know that your submatrix rectangle at least has an area of 8, and no more than 12. In fact, for this matrix, the ideal works out to 8. Potential shapes can be seen...

(res = MapIndexed[(temp = #1*Times @@ #2;
If[temp < lowerBound, 0, temp]) &, ms, {2}]) // MatrixForm


$$\left( \begin{array}{ccccc} 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 8 & 0 \\ 0 & 0 & 9 & 0 & 0 \\ 0 & 8 & 12 & 0 & 0 \\ 0 & 10 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 8 & 0 & 0 & 0 & 0 \\ \end{array} \right)$$

So for example, you need to see if combinations of 4 (from 8) rows and 3 (from 5) columns result in a successful submatrix. And you can get the subsets to search over (corresponding to the coordinates of the non-zero elements in the matrix above).

res2 = MapIndexed[Boole@Positive[#1]*#2 &, res, {2}];
res3 = Partition[Flatten[res2], 2];
res4 = Select[res3, # != {0, 0} &]
(* {{2, 4}, {3, 3}, {4, 2}, {4, 3}, {5, 2}, {8, 1}} *)


Lots fewer potential subsets to search through.

• Thanks Mike, great solution, many thanks ... I suspect I will run into constraints when applying the Subset expansions .. but so far so good. – Sander Jan 25 '19 at 12:00

Update: This answer is not correct but for referencing, because the Subsets does not give all the possible slices of the matrix.

mat = {{1, 1, 1, 1, 0, 0}, {1, 1, 1, 1, 1, 1}, {0, 0, 0, 1, 1, 1}, {1,
1, 0, 1, 1, 1}, {1, 1, 1, 1, 1, 1}, {1, 1, 1, 1, 1, 1}};

mat//MatrixForm//TeXForm


$$\left( \begin{array}{cccccc} 1 & 1 & 1 & 1 & 0 & 0 \\ 1 & 1 & 1 & 1 & 1 & 1 \\ 0 & 0 & 0 & 1 & 1 & 1 \\ 1 & 1 & 0 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 & 1 & 1 \\ \end{array} \right)$$

lst = Subsets@mat;

result = DeleteDuplicates@(MatrixForm /@ Select[lst, DeleteDuplicates@Flatten@# == {1} &]) // Sort;
result // TeXForm


$$\left\{\left( \begin{array}{cccccc} 1 & 1 & 1 & 1 & 1 & 1 \\ \end{array} \right),\left( \begin{array}{cccccc} 1 & 1 & 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 & 1 & 1 \\ \end{array} \right),\left( \begin{array}{cccccc} 1 & 1 & 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 & 1 & 1 \\ \end{array} \right)\right\}$$

Is this ok? I'm not sure to apply Transpose to the last one of the result.

• Thanks Jerry, there seem to be two issues: 1 the largest sub-matrix in the example is 5x3, your answer results in (after transpose) a 6x3 sub-matrix; Also, I am concerned the Subset will explode beyond memory capacity once we work with large matrices? 2. I would like to recover the coordinates of where the sub-matrix is residing. – Sander Jan 24 '19 at 11:00
• @Sander Here is a similar question in other language, see geeksforgeeks.org/… – Jerry Jan 24 '19 at 11:08

This method will omit those submatrices smaller than 2*2 :

hilight[mat_] :=
Module[{pmat, i = 0, originMat, last = 0},
originMat = pmat = Image[mat];
While[Total[Flatten[ImageData[pmat, "Byte"]]] != 0,
pmat = Erosion[pmat, 1, Padding -> 0]; i++;
If[AnyTrue[
Values[ComponentMeasurements[
pmat, {"BoundingBoxArea", "Count"}]], Apply[Equal]],
last = i]];
HighlightImage[originMat,
Dilation[Erosion[Image[originMat], last, Padding -> 0], last]]]


As the mat from origin question

hilight[mat]


As the initMat from Okkes Dulgerci's answer:

SeedRandom@8
{m, n} = {15, 20};
(initMat = RandomChoice[{20, 1} -> {1, 0}, {m, n}]);
hilight[initMat]


• nice!! It seems that you need a more stringent condition; if your current condition is satisfied for a small isolated submatrix of 1s, you miss a larger block in other parts of the input matrix. Try, for example, mat2 = ArrayPad[mat, 1, 1]`. – kglr Jan 26 '19 at 6:51
• oh... Thanks for your reminder...@kglr – yode Jan 26 '19 at 7:11