# 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? – J. M. will be back soon Feb 25 '16 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 Feb 25 '16 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. – J. M. will be back soon Feb 25 '16 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. – rhermans Feb 25 '16 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 Feb 25 '16 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 Feb 27 '16 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... – unlikely Feb 27 '16 at 21:46
• @tavr I added an updated version to force the removal of some rows/columns – unlikely Feb 27 '16 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 Feb 27 '16 at 23:44
• Genious, +1. I have to read up on LinearProgramming – LLlAMnYP Mar 3 '16 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 Feb 27 '16 at 20:51
• @tavr Please check my other answer with LinearProgramming but don't expect some magic... – unlikely Feb 27 '16 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 Feb 27 '16 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 – LLlAMnYP Feb 27 '16 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". – Schopenhauer Feb 27 '16 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? – rhermans Feb 25 '16 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 – garej Feb 25 '16 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 Feb 25 '16 at 16:26
• @tavr Please add all necessary clarification in the question, as it is, is about to be closed for been unclear. – rhermans Feb 25 '16 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 – rhermans Feb 25 '16 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."... – unlikely Mar 2 '16 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 Mar 2 '16 at 7:38