# Given a large binary matrix, find the largest submatrix containing non-zero elements

I'm looking for a way to reduce a binary matrix containing zeros at some positions into a matrix that contains no zeros by deleting rows and columns of the original matrix until only non-zero values are left.

There is an important constraint on this procedure in that I also want to minimize the total number of elements from the original matrix that are removed to arrive at the solution.

For example, consider this matrix:

SeedRandom[0];
(initmat = RandomChoice[{8, 1} -> {1, 0}, {10, 12}]) // MatrixForm


The optimal solution in this case is to remove rows 1,2 and 4 and columns 10 and 12 to end up with a 7x10 matrix of ones. As you can see, the resulting submatrix need not have been a contiguous block of the original matrix.

My first stab at a solution is along the lines of:

NestWhile[Transpose[Rest[Sort[#]]] &, initmat, Cases[#, 0, 2] != {} &]


However, this is neither efficient nor optimal in most cases, and could probably be improved, though the general version of the problem I'm describing has been shown to be NP-hard.

• Can you maybe give example inputs and expected outputs, to help us further specify what you need? Commented Feb 25, 2016 at 8:53
• One example is initmat, above. Another, trivial, example is a matrix which has one row and one column consisting of nothing but zeros, and no zeros anywhere else. The expected output in that case would be a matrix with the zero-containing row and column deleted. The nestwhile approach works for the trivial case, but not for initmat: it deletes more of the original matrix than it needs to (it produces a matrix that is 7x9, whereas the optimal solution is 6x11)
– tavr
Commented Feb 25, 2016 at 9:12
• That's what I meant by my request. Can you supply a particular initmat (RandomChoice[] will of course give different results on different machines), and show what your "optimal solution" looks like for that particular input. Commented Feb 25, 2016 at 10:28
• A particular initmat is provided, as the code is preceded by SeedRandom[0]; now RandomChoice will always give the same output first time is called. So far so good. What we do need now is the output, preferable step by step to understand unambiguously what is the procedure you mean. Commented Feb 25, 2016 at 11:36
• I wouldn't say it's a duplicate of that. To your earlier question, I don't really have an algorithm in place that delivers the optimal output. From eyeballing the initmat, I can see that the solution is a 7x10 matrix of 1's (rows 1,2,4 deleted; columns 10 and 12 deleted). What I need is a formal approach that gets me to that solution.
– tavr
Commented Feb 25, 2016 at 16:43

Because you want to minimize the number of elements removed, a natural way to do this is with constrained optimization functions like Minimize and LinearProgramming.

In another answer I used the the easy-to-implement, and slow, method with Minimize; here the faster method of LinearProgramming.

LinearProgramming require that the objective function and constraints are linear in the decision variables. So we need to rewrite the objective function in some linear form. This is possible adding some more (many more!) decision variables and constraints.

Given:

• the matrix $A$, with dimensions $(m,n)$ and elements $a_{i,j}$

you want to determine the value of variables ($i=1,\ldots,m$, $j=1,\ldots,n$):

• $e_{i,j} \in \{0,1\}$ (it's $1$ if we want/need to remove the element $(i,j)$, and $0$ otherwise)
• $r_i \in \{0,1\}$ (it's $1$ if we remove row $i$, and $0$ otherwise)
• $c_j \in \{0,1\}$ (it's $1$ if we remove column $j$, and $o$ otherwise)

subject to:

• $e_{i,j}=1$ for every $(i,j)$ such that $a_{i,j}=0$
• $r_i+c_j \ge e_{i,j}$ i.e. $e_{i,j}=1 \Rightarrow (r_i=1 \vee c_j=1)$
• $e_{i,j} \ge r_i$ i.e. $r_i = 1 \Rightarrow e_{i,j}=1$
• $e_{i,j} \ge c_j$ i.e. $c_j = 1 \Rightarrow e_{i,j}=1$

minimizing the objective function

• $\sum_{i=1}^m \sum_{j=1}^n e_{i,j}$.

Note you can also enforce $e_{i,j} \ge 1$ instead of $e_{i,j} = 1$ because of the objective function. Indeed it's this we do in the following code. But we need to build the arguments of LinearProgramming:

• the cost vector
• the constraints matrix
• the right-hand side of the constraints
• the range of the decision variables
• the domain of the decision variables

I packaged the process in the following function:

zeroFreeSubmatrix[A_?MatrixQ] :=
Module[{m, n, e, r, c, vars, constraints, bm, solution},
{m, n} = Dimensions[A];
vars = Flatten@{Array[e, {m, n}], Array[r, m], Array[c, n]};
constraints = Flatten@{
Thread[e @@@ Position[A, 0, {2}] == 1],
Table[{r[i] + c[j] >= e[i, j], e[i, j] >= r[i],
e[i, j] >= c[j]}, {i, m}, {j, n}]
};
bm = CoefficientArrays[Equal @@@ constraints, vars];
vars /. {_e -> 1, (_r | _c) -> 0},
bm[[2]], -bm[[1]],
Table[{0, 1}, Length@vars], Integers
]];

Sort@Cases[solution, (#[i_] -> 0) :> i] & /@ {r, c}
]


The usage:

SeedRandom[0];
A = RandomChoice[{8, 1} -> {1, 0}, {15, 20}];
{rows, cols} =
zeroFreeSubmatrix[A] // Timing // EchoFunction["Timing:", First] //
Last
A[[rows, cols]] // MatrixForm


During evaluation of LinearProgramming::lpip: Warning: integer linear programming will use a machine-precision approximation of the inputs. >>

Timing: 0.234375

{{1, 3, 4, 6, 8, 9, 10, 11, 12, 13, 14}, {2, 3, 4, 7, 9, 12, 14, 15, 17, 18, 19, 20}}

Faster than the simpler approach with Minimize.

There is also a way to directly build the arguments of LinearProgramming with SparseArray but I don't think in this case it deserve the effort.

# Update

If you are sure that some rows or columns will be removed in any optimal solution you can use the following generalization of the previous routine.

zeroFreeSubmatrix[A_?MatrixQ, rl : _?VectorQ : {},
cl : _?VectorQ : {}] :=
Module[{m, n, e, r, c, vars, constraints, bm, solution},
{m, n} = Dimensions[A];
vars = Flatten@{Array[e, {m, n}], Array[r, m], Array[c, n]};
constraints = Flatten@{
Thread[e @@@ Position[A, 0, {2}] == 1],
Table[{r[i] + c[j] >= e[i, j], e[i, j] >= r[i],
e[i, j] >= c[j]}, {i, m}, {j, n}],
};
bm = CoefficientArrays[Equal @@@ constraints, vars];
vars /. {_e -> 1, (_r | _c) -> 0},
bm[[2]], -bm[[1]],
Table[{0, 1}, Length@vars], Integers
], LinearProgramming::lpip]];

Sort@Cases[solution, (#[i_] -> 0) :> i] & /@ {r, c}
]


For example if you are sure that the rows $2,5$ and colum $1,5$ are to be removed you can use:

{rows, cols} = zeroFreeSubmatrix[A, {2, 5}, {1, 5}]
A[[rows, cols]] // Dimensions // Apply[Times]


{{1, 3, 4, 6, 8, 9, 10, 11, 12, 13, 14}, {2, 3, 4, 7, 9, 12, 14, 15, 17, 18, 19, 20}}

132

and you get the (same) optimal solution.

But if you are wrong you can get a different, maybe sub-optimal, solution (the optimal solution where you enforce the removal of that rows/columns). For example:

{rows, cols} = zeroFreeSubmatrix[A, {1, 3}, {2, 3}]
A[[rows, cols]] // Dimensions // Apply[Times]


{{4, 5, 6, 9, 10, 11, 13, 15}, {1, 4, 5, 6, 7, 9, 12, 13, 14, 15, 16, 17, 18, 19}}

112

• Suppose I had some heuristic that told me with near-certainty that some rows and columns will be removed in any optimal solution. How could I inject that information into this approach?
– tavr
Commented Feb 27, 2016 at 21:34
• @tavr If you add a constraint like r[i]==1 you get the optimal solution between those where the row $i$ is removed... Commented Feb 27, 2016 at 21:46
• @tavr I added an updated version to force the removal of some rows/columns Commented Feb 27, 2016 at 23:00
• Thanks. I can see that this works, and with enough hints about removals, it's quite fast. Probably the best answer so far.
– tavr
Commented Feb 27, 2016 at 23:44
• Genious, +1. I have to read up on LinearProgramming Commented Mar 3, 2016 at 9:42

Because you want to minimize the number of elements removed, a natural way to do this is with constrained optimization functions like Minimize and LinearProgramming.

In this answer I use the basic idea, the easy-to-implement, and the slow method with Minimize; in another answer the faster method of LinearProgramming which is a bit less straightforward.

Given:

• a matrix $A$, with dimensions $(m,n)$ and elements $a_{i,j}$

you want to determine the value of variables:

• $r_i\in\{0,1\}$, $i=1,\ldots,m$, (it's $1$ if the row $i$ of $A$ has to be removed, and $0$ otherwise)
• $c_j\in\{0,1\}$, $j=1,\ldots,n$, (it's $1$ if the column $j$ of $A$ has to be removed, and $0$ otherwise)

subject to:

• $r_i+c_j>0$ for every $(i,j)$ s.t. $a_{i,j}=0$, i.e. $a_{i,j}=0 \Rightarrow (r_i > 0 \vee c_j > 0)$

minimizing the objective function:

• $n \sum_{i=1}^m r_i + (m-\sum_{i=1}^m r_i) \sum_{j=1}^n c_j$.

This is easy to translate with with Minimize:

SeedRandom[0];
(A = RandomChoice[{8, 1} -> {1, 0}, {10, 12}]) // MatrixForm

{m, n} = Dimensions[initmat];

vars = Flatten@{
Table[r[i], {i, m}],
Table[c[j], {j, n}]
};

constraints = Flatten@{
r@#1 + c@#2 > 0 & @@@ Position[A, 0, {2}],
};

objective =
n Sum[r[i], {i, m}] + (m - Sum[r[i], {i, m}])*Sum[c[j], {j, n}] //
Expand;

({opt, solution} =
Minimize[{objective, constraints}, vars, Integers]) // Timing

A[[Cases[solution, (r[i_] -> 0) :> i],
Cases[solution, (c[i_] -> 0) :> i]]] // MatrixForm


{28.3125, {50, {r2 -> 1, r2 -> 1, r[3] -> 0, r[4] -> 1, r[5] -> 0, r[6] -> 0, r[7] -> 0, r[8] -> 0, r[9] -> 0, r[10] -> 0, c2 -> 0, c2 -> 0, c[3] -> 0, c[4] -> 0, c[5] -> 0, c[6] -> 0, c[7] -> 0, c[8] -> 0, c[9] -> 0, c[10] -> 1, c[11] -> 0, c[12] -> 1}}}

As expected, it's not fast, mostly because of the objective function, which is not linear in the decision variables.

• Can you describe in more detail how to do this faster? As it is, the approach is not very scalable (# elements in objective function \[Proportional] to mxn)
– tavr
Commented Feb 27, 2016 at 20:51
• @tavr Please check my other answer with LinearProgramming but don't expect some magic... Commented Feb 27, 2016 at 21:03

There is a brute force method, but it actually avoids matrix operations. The problem is, that once there are too many rows with zeros, the list returned by Subset gets far too large. Anyway, Taking initmat as the matrix we're trying to reduce:

pos = GatherBy[Position[initmat, 0], First] /. {a : {{x_, _} ..} :> {x, Last /@ a}}


Gets us the list of positions of zeros in the form

{1, {5, 6}}, {2, {1, 9, 10}} ... }


meaning "row 1 has zeros in columns 5 and 6, row 2 has zeros in columns 1, 9, and 10..."

l = Length[pos]; (* 6 *)


For the example matrix given in the OP

SeedRandom[0]; (initmat = RandomChoice[{8, 1} -> {1, 0}, {10, 12}])


if we delete row 1, for example, we do not need to delete column 6 anymore. Let's find out, which columns are left to delete after we delete certain rows containing zeros:

rowscols =
{pos[[#, 1]], Union @@ (Last /@ Delete[pos, Transpose@{#}])} & /@ Subsets[Range[l]]


This will give a list in the form

{{rows1,cols1}, {rows2, cols2}...}


Where rowsi and colsi are list of row and column numbers that should be deleted. Each combination is "optimal" in the sense that after a row with a zero is deleted, the column of that zero which does not contain any more zeros will not be deleted.

The size of the matrix (number of retained elements) after each deleting a given combination is given by

Times@@(Dimensions[initmat] - {Length[rowsi],Length[colsi]})


So we can find the possible sizes like so:

sizes = (Times@@(Dimensions[initmat] - {Length[#1],Length[#2]}))&@@@rowscols


Then

max = Position[sizes,Max[sizes]];
optimal = rowscols[[Flatten[max]]];


and the reduced matrix (or matrices, if there's more than one possibility) will be

out = Transpose[
Delete[Transpose[Delete[initmat, List /@ First@#]],
List /@ Last@#]] & /@ optimal


And here's the entire code in a single Module. It's quite procedural and can likely be better written.

reduceMatrix[initmat_] :=
Module[{pos, l, rowscols, sizes, max, optimal, out},
pos = GatherBy[Position[initmat, 0],
First] /. {a : {{x_, _} ..} :> {x, Last /@ a}};
l = Length[pos];
rowscols = {pos[[#, 1]],
Union @@ (Last /@ Delete[pos, Transpose@{#}])} & /@
Subsets[Range[l]];
sizes = (Times @@ (Dimensions[initmat] - {Length[#1],
Length[#2]})) & @@@ rowscols;
max = Position[sizes, Max[sizes]];
optimal = rowscols[[Flatten[max]]];
out = Transpose[
Delete[Transpose[Delete[initmat, List /@ First@#]],
List /@ Last@#]] & /@ optimal]

• I think this is a general solution, but it is way too expensive, computationally. I'm dealing with 500x500 arrays typically, with something like 15% of the elements zeros, so need a solution that will be faster. I was hoping that the Subsets[Range[l]] could be reduced to a small set, either through heuristics, or otherwise.
– tavr
Commented Feb 27, 2016 at 2:46
• @tavr I believe, a start could be to partition the matrix into groups of a small number of rows and reduce each separately, then somehow find the optimal union of reductions. Can't say much more from mobile right now though Commented Feb 27, 2016 at 16:39
• @LLlAMnYP See The Consecutive Ones Submatrix Problem for Sparse Matrices: goo.gl/bR6wls for a formal description of the problem. So the question could say "Given a 500x500 matrix to efficiently find the biggest submatrix containing non-zero elements". Commented Feb 27, 2016 at 17:35

This does not answer the OP's question.

I take all the positions where there are zeros using Position

Then get a sorted list of unique rows and columns where zeros were found.

Map[Sort@*DeleteDuplicates, Transpose[Position[m, s]]]


Then give only the columns and rows that are not in that list

## Approach 1

f[m_, s_] := Part[m,
Sequence @@
Apply[
Complement,
Transpose@{
Map[Range, Dimensions[m]]
, Map[Sort@*DeleteDuplicates, Transpose[Position[m, s]]]
}, 1]
]


## Approach 2

rDelete[m_, row_, col_] := Delete[
Transpose[
Delete[
Transpose[m]
, List /@ col
]
], List /@ row
]

g[m_, s_] := rDelete[m,
Sequence @@ Map[Sort@*DeleteDuplicates, Transpose[Position[m, s]]]
]


## Results

• @tavr, is this what you needed? Commented Feb 25, 2016 at 14:39
• OP seems to want deleting minimal number of rows and columns possible (so in comment OP claims for 6X11 output, not 4X5). I.e. to choose order of deleting so as to eliminate maximum 0 during each iteration. IMHO Commented Feb 25, 2016 at 15:51
• @garej is right. I need to delete the minimum number of rows and columns necessary to achieve a matrix with no zeros. So, for example, if you've already deleted a column containing a zero, there is no longer a need to delete the row that contained that zero (unless that row contained another zero).
– tavr
Commented Feb 25, 2016 at 16:26
• @tavr Please add all necessary clarification in the question, as it is, is about to be closed for been unclear. Commented Feb 25, 2016 at 16:29
• @tavr Ok, so now that implies there is more than one way of doing it. How do you define the correct way? Please edit and explain in your question Commented Feb 25, 2016 at 16:30

I have not rigorously tested this procedure, but I believe, it might just work.

In fact, this approach does not work!
I leave it here, because some snippets of code may be useful in the correct algorithm, but this appears to be a sort of traveling salesman kind of problem and I'm not sure, if there's always a perfect algorithm.

Let your matrix be of size n*m and stored as mat. Then initialize a padded matrix like so:

m = SparseArray[{{i_, j_} -> Max[i, j] - 1}, {n+1, m+1}];
m[[2;;,2;;]] = mat;
m = Normal@m;


That will basically get you a matrix with a header row and column, numbering it from 1 to n (m). Then repeated execution of the block of code below should get you to where you want

(m = Transpose@m);
m[[2 ;;]] = SortBy[-Counts[Rest@#][0] &]@m[[2 ;;]];
(m = Transpose@m);
m[[2 ;;]] = SortBy[-Counts[Rest@#][0] &]@m[[2 ;;]];
m // MatrixForm
If[(Length[m[[2]]] - 1 - Counts[m[[2]]][0])/Counts[m[[2]]][0] <= (
Length[m[[All, 2]]] - 1 - Counts[m[[All, 2]]][0])/
Counts[m[[All, 2]]][0], m = Delete[m, 2],
m = Transpose[Delete[Transpose[m], 2]], Print["Done!"];


Once the zeros are gone, you can use the header rows and columns to return to the original ordering.

There is no need to use pure mathematical methods to do that. Since this is a zero-one matrix, the morphological method for image processing will be unbelievably faster than others, because you can use GPU power to deal with super large matrices. Here, I simplified this problem to find the largest all-one square matrix to show this idea. Firstly, I added a zero boundary to this matrix, and transfer it into an image, which show one as white pixel and zero as black one. Then I iteratively use HissMissTransform to shrink the all-one square matrices in this image, until no white pixels found in this figure. Then use the last second figure Here is the code:

    Module[{image = ImagePad[Image@initmat, 1, 0], d = 0, imagelast},
While[Max@ImageData@image == 1, d++; imagelast = image;
image = HitMissTransform[image, ConstantArray[1, {2, 2}]];
Print@Image@image]; {d,
imagelast, (# - {1, 1}) & /@ Position[ImageData@imagelast, 1]}]


Then you will get the results, the size of the largest all-one square matrix (d=6), images that will show this algorithm and the positions of those matrices.

All the intermediate outputs are as below:

So to find the largest all-one matrix, I continued to extend the mask along the row or column direction. Here is the code:

    Module[{image0 = ImagePad[Image@initmat, 1, 0], image, d = 0,
imagelast, dx = 0, dy = 0, imagex, imagey}, image = image0;
Print["find largest square matrix"];
While[Max@ImageData@image == 1, d++; imagelast = image;
Print@Image@image;
image = HitMissTransform[image, ConstantArray[1, {2, 2}]];];
imagey =
imagex = HitMissTransform[image0, ConstantArray[1, {d, d}]];
Print["find largest matrix: row < column"];
While[Max@ImageData@imagex == 1, dx++;
imagex = HitMissTransform[image0, ConstantArray[1, {d, d + dx}]];
Print@Image@imagex;];
dx--;
Print["find largest matrix: row > column"];
While[Max@ImageData@imagey == 1, dy++;
imagey = HitMissTransform[image0, ConstantArray[1, {d + dy, d}]];
Print@Image@imagey;];
dy--; If[dx > dy, {d + dx, d}, {d, d + dy}]]


The last vector showed the size of the largest matrix. To find its position, you just need to re-do the HitMissTransform.

• I admit this is interesting and I don't fully understand. But doesn't this find the largest contiguous (square) submatrix? As of OP: "As you can see, the resulting submatrix need not have been a contiguous block of the original matrix."... Commented Mar 2, 2016 at 7:29
• Second what unlikely said in his comment. This might not seem like a big deal, and "close enough", but for large binary matrices with relatively uniform density of zeros, it will spell the difference between a solution that leaves only 1% of the original elements (in the contiguous case), and one that leaves 25% (in the non-contiguous). I do like the approach though. I wonder if it can be adapted to the problem at hand.
– tavr
Commented Mar 2, 2016 at 7:38

This question and its answers have been really useful for me, especially @unlikely's Linear Programming strategy (I implemented it in python for my needs). I also thought about this : what if the right permutations on columns and rows (based on the number of 0s by row/column) reduce the original problem to the search of the largest (contiguous) rectangle in the processed matrix:

I tested it and it seems to give the same results as the linear programming solution.

• Can you add the Mathematica code corresponding to your approach? Commented May 9, 2023 at 15:20
• I haven't coded in Mathematica for almost ten years, but I can show my python code if you wish (I know it is a bit off tropic). Commented May 10, 2023 at 10:54
• Since my last comment, I have been able to refine and assess my approach (which is actually a heuristic) and created a package (unfortunately, in Python). I present my algorithm in this notebook (github.com/CyrilJl/OptiMask/blob/main/notebooks/Optimask.ipynb). Commented Nov 26, 2023 at 17:36