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

Question

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.

Answer 1

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];
  solution = Thread[vars -> LinearProgramming[
      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}],
     Thread[r /@ rl == 1],
     Thread[c /@ cl == 1]
     };
  bm = CoefficientArrays[Equal @@@ constraints, vars];
  solution = Thread[vars -> Quiet[LinearProgramming[
       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

Answer 2

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}],
    Thread[0 <= vars <= 1]
    };

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.

Answer 3

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]

Answer 4

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

Part, Complement

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

---

Approach 2

Delete, Transpose

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

Answer 5

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.

Answer 6

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.

Answer 7

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.