Subsetting non-NA

Question

I have a matrix in which every row has at least one NA cell, and every column has at least one NA cell as well. What I need is to find the largest subset of this matrix that contains no NAs.

For example, for this matrix A

A <- 
  structure(c(NA, NA, NA, NA, 2L, NA,
              1L, 1L, 1L, 0L, NA, NA,
              1L, 8L, NA, 1L, 1L, NA, 
              NA, 1L, 1L, 6L, 1L, 3L, 
              NA, 1L, 5L, 1L, 1L, NA),
            .Dim = c(6L, 5L),
            .Dimnames = 
              list(paste0("R", 1:6),
                   paste0("C", 1:5)))

A
    C1  C2  C3  C4  C5
R1  NA  1   1   NA  NA
R2  NA  1   8   1   1
R3  NA  1   NA  1   5
R4  NA  0   1   6   1
R5  2   NA  1   1   1
R6  NA  NA  NA  3   NA

There are two solutions (8 cells): A[c(2, 4), 2:5] and A[2:5, 4:5], though finding just one valid solution is enough for my purposes. The dimensions of my actual matrix are 77x132.

Being a noob, I see no obvious way to do this. Could anyone help me with some ideas?

Answer 1

1) optim In this approach we relax the problem to a continuous optimization problem which we solve with optim.

The objective function is f and the input to it is a 0-1 vector whose first nrow(A) entries correspond to rows and whose remaining entries correspond to columns. f uses a matrix Ainf which is derived from A by replacing the NAs with a large negative number and the non-NAs with 1. In terms of Ainf the negative of the number of elements in the rectangle of rows and columns corresponding to x is -x[seq(6)] %*% Ainf %*$ x[-seq(6)] which we minimize as a function of x subject to each component of x lying between 0 and 1.

Although this is a relaxation of the original problem to continuous optimization it seems that we get an integer solution, as desired, anyways.

Actually most of the code below is just to get the starting value. To do that we first apply seriation. This permutes the rows and columns giving a more blocky structure and then in the permuted matrix we find the largest square submatrix.

In the case of the specific A in the question the largest rectangular submatrix happens to be square and the starting values are already sufficiently good that they produce the optimum but we will perform the optimization anyways so it works in general. You can play around with different starting values if you like. For example, change k from 1 to some higher number in largestSquare in which case largestSquare will return k columns giving k starting values which can be used in k runs of optim taking the best.

If the starting values are sufficiently good then this should produce the optimum.

library(seriation) # only used for starting values

A.na <- is.na(A) + 0

Ainf <- ifelse(A.na, -prod(dim(A)), 1) # used by f
nr <- nrow(A) # used by f
f <- function(x) - c(x[seq(nr)] %*% Ainf %*% x[-seq(nr)])

# starting values

# Input is a square matrix of zeros and ones.
# Output is a matrix with k columns such that first column defines the
# largest square submatrix of ones, second defines next largest and so on.
# Based on algorithm given here:
# http://www.geeksforgeeks.org/maximum-size-sub-matrix-with-all-1s-in-a-binary-matrix/
largestSquare <- function(M, k = 1) {
  nr <- nrow(M); nc <- ncol(M)
  S <- 0*M; S[1, ] <- M[1, ]; S[, 1] <- M[, 1]
  for(i in 2:nr) 
    for(j in 2:nc)
      if (M[i, j] == 1) S[i, j] = min(S[i, j-1], S[i-1, j], S[i-1, j-1]) + 1
  o <- head(order(-S), k)
  d <- data.frame(row = row(M)[o], col = col(M)[o], mx = S[o])
  apply(d, 1, function(x) { 
    dn <- dimnames(M[x[1] - 1:x[3] + 1, x[2] - 1:x[3] + 1])
    out <- c(rownames(M) %in% dn[[1]], colnames(M) %in% dn[[2]]) + 0
    setNames(out, unlist(dimnames(M)))
  })
}
s <- seriate(A.na)
p <- permute(A.na, s)
# calcualte largest square submatrix in p of zeros rearranging to be in A's  order
st <- largestSquare(1-p)[unlist(dimnames(A)), 1]

res <- optim(st, f, lower = 0*st, upper = st^0, method = "L-BFGS-B")

giving:

> res
$par
R1 R2 R3 R4 R5 R6 C1 C2 C3 C4 C5 
 0  1  1  1  0  0  0  1  0  1  1 

$value
[1] -9

$counts
function gradient 
       1        1 

$convergence
[1] 0

$message
[1] "CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL"

2) GenSA Another possibility is to repeat (1) but instead of using optim use GenSA from the GenSA package. It does not require starting values (although you can provide a starting value using the par argument and this might improve the solution in some cases) so the code is considerably shorter but since it uses simulated annealing it can be expected to take substantially longer to run. Using f (and nr and Ainf which f uses) from (1). Below we try it without a starting value.

library(GenSA)
resSA <- GenSA(lower = rep(0, sum(dim(A))), upper = rep(1, sum(dim(A))), fn = f)

giving:

> setNames(resSA$par, unlist(dimnames(A)))
R1 R2 R3 R4 R5 R6 C1 C2 C3 C4 C5 
 0  1  1  1  0  0  0  1  0  1  1 

> resSA$value
[1] -9

Answer 2

I have a solution, but it doesn't scale very well:

findBiggestSubmatrixNonContiguous <- function(A) {
    A <- !is.na(A); ## don't care about non-NAs
    howmany <- expand.grid(nr=seq_len(nrow(A)),nc=seq_len(ncol(A)));
    howmany <- howmany[order(apply(howmany,1L,prod),decreasing=T),];
    for (ri in seq_len(nrow(howmany))) {
        nr <- howmany$nr[ri];
        nc <- howmany$nc[ri];
        rcom <- combn(nrow(A),nr);
        ccom <- combn(ncol(A),nc);
        comcom <- expand.grid(ri=seq_len(ncol(rcom)),ci=seq_len(ncol(ccom)));
        for (comi in seq_len(nrow(comcom)))
            if (all(A[rcom[,comcom$ri[comi]],ccom[,comcom$ci[comi]]]))
                return(list(ri=rcom[,comcom$ri[comi]],ci=ccom[,comcom$ci[comi]]));
    }; ## end for
    NULL;
}; ## end findBiggestSubmatrixNonContiguous()

---

It's based on the idea that if the matrix has a small enough density of NAs, then by searching for the largest submatrices first, you'll be likely to find a solution fairly quickly.

The algorithm works by computing a cartesian product of all counts of rows and counts of columns that could be indexed out of the original matrix to produce the submatrix. The set of pairs of counts is then decreasingly ordered by the size of the submatrix that would be produced by each pair of counts; in other words, ordered by the product of the two counts. It then iterates over these pairs. For each pair, it computes all combinations of row indexes and column indexes that could be taken for that pair of counts, and tries each combination in turn until it finds a submatrix that contains zero NAs. Upon finding such a submatrix, it returns that set of row and column indexes as a list.

The result is guaranteed to be correct because it tries submatrix sizes in decreasing order, so the first one it finds must be the biggest (or tied for the biggest) possible submatrix that satisfies the condition.

---

## OP's example matrix
A <- data.frame(C1=c(NA,NA,NA,NA,2L,NA),C2=c(1L,1L,1L,0L,NA,NA),C3=c(1L,8L,NA,1L,1L,NA),C4=c(NA,1L,1L,6L,1L,3L),C5=c(NA,1L,5L,1L,1L,NA),row.names=c('R1','R2','R3','R4','R5','R6'));
A;
##    C1 C2 C3 C4 C5
## R1 NA  1  1 NA NA
## R2 NA  1  8  1  1
## R3 NA  1 NA  1  5
## R4 NA  0  1  6  1
## R5  2 NA  1  1  1
## R6 NA NA NA  3 NA
system.time({ res <- findBiggestSubmatrixNonContiguous(A); });
##    user  system elapsed
##   0.094   0.000   0.100
res;
## $ri
## [1] 2 3 4
##
## $ci
## [1] 2 4 5
##
A[res$ri,res$ci];
##    C2 C4 C5
## R2  1  1  1
## R3  1  1  5
## R4  0  6  1

We see that the function works very quickly on the OP's example matrix, and returns a correct result.

---

randTest <- function(NR,NC,probNA,seed=1L) {
    set.seed(seed);
    A <- replicate(NC,sample(c(NA,0:9),NR,prob=c(probNA,rep((1-probNA)/10,10L)),replace=T));
    print(A);
    print(system.time({ res <- findBiggestSubmatrixNonContiguous(A); }));
    print(res);
    print(A[res$ri,res$ci,drop=F]);
    invisible(res);
}; ## end randTest()

I wrote the above function to make testing easier. We can call it to test a random input matrix of size NR by NC, with a probability of choosing NA in any given cell of probNA.

---

Here are a few trivial tests:

randTest(8L,1L,1/3);
##      [,1]
## [1,]   NA
## [2,]    1
## [3,]    4
## [4,]    9
## [5,]   NA
## [6,]    9
## [7,]    0
## [8,]    5
##    user  system elapsed
##   0.016   0.000   0.003
## $ri
## [1] 2 3 4 6 7 8
##
## $ci
## [1] 1
##
##      [,1]
## [1,]    1
## [2,]    4
## [3,]    9
## [4,]    9
## [5,]    0
## [6,]    5

---

randTest(11L,3L,4/5);
##       [,1] [,2] [,3]
##  [1,]   NA   NA   NA
##  [2,]   NA   NA   NA
##  [3,]   NA   NA   NA
##  [4,]    2   NA   NA
##  [5,]   NA   NA   NA
##  [6,]    5   NA   NA
##  [7,]    8    0    4
##  [8,]   NA   NA   NA
##  [9,]   NA   NA   NA
## [10,]   NA    7   NA
## [11,]   NA   NA   NA
##    user  system elapsed
##   0.297   0.000   0.300
## $ri
## [1] 4 6 7
##
## $ci
## [1] 1
##
##      [,1]
## [1,]    2
## [2,]    5
## [3,]    8

---

randTest(10L,10L,1/3);
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
##  [1,]   NA   NA    0    3    8    3    9    1    6    NA
##  [2,]    1   NA   NA    4    5    8   NA    8    2    NA
##  [3,]    4    2    5    3    7    6    6    1    1     5
##  [4,]    9    1   NA   NA    4   NA   NA    1   NA     9
##  [5,]   NA    7   NA    8    3   NA    5    3    7     7
##  [6,]    9    3    1    2    7   NA   NA    9   NA     7
##  [7,]    0    2   NA    7   NA   NA    3    8    2     6
##  [8,]    5    0    1   NA    3    3    7    1   NA     6
##  [9,]    5    1    9    2    2    5   NA    7   NA     8
## [10,]   NA    7    1    6    2    6    9    0   NA     5
##    user  system elapsed
##   8.985   0.000   8.979
## $ri
## [1]  3  4  5  6  8  9 10
##
## $ci
## [1]  2  5  8 10
##
##      [,1] [,2] [,3] [,4]
## [1,]    2    7    1    5
## [2,]    1    4    1    9
## [3,]    7    3    3    7
## [4,]    3    7    9    7
## [5,]    0    3    1    6
## [6,]    1    2    7    8
## [7,]    7    2    0    5

I don't know an easy way of verifying if the above result is correct, but it looks good to me. But it took almost 9 seconds to generate this result. Running the function on moderately larger matrices, especially a 77x132 matrix, is probably a lost cause.

Waiting to see if someone can come up with a brilliant efficient solution...