Подмножество не-NA

1) optim В этом подходе мы сводим проблему к задаче непрерывной оптимизации, которую мы решаем с помощью optim.

Целевой функцией является f, а входными данными для нее является вектор 0-1, первые nrow(A) элемента которого соответствуют строкам, а остальные элементы соответствуют столбцам. f использует матрицу Ainf, полученную из A путем замены NA большим отрицательным числом, а не-NA на 1. С точки зрения Ainf отрицательное число элементов в прямоугольнике строк и столбцов, соответствующих x, равно -x[seq(6)] %*% Ainf %*$ x[-seq(6)] который мы минимизируем как функцию x при условии, что каждый компонент x лежит между 0 и 1.

Хотя это смягчение исходной задачи до непрерывной оптимизации, кажется, что мы в любом случае получаем целочисленное решение, как и хотелось.

На самом деле большая часть кода ниже предназначена только для получения начального значения. Для этого мы сначала применяем сериализацию. Это переставляет строки и столбцы, создавая более блочную структуру, а затем в переставленной матрице мы находим наибольшую квадратную подматрицу.

В случае конкретного A в вопросе самая большая прямоугольная подматрица оказывается квадратной, и начальные значения уже достаточно хороши, чтобы дать оптимум, но мы все равно выполним оптимизацию, так что в целом это работает. Вы можете поиграть с различными начальными значениями, если хотите. Например, измените k с 1 на некоторое большее число в largestSquare, и в этом случае largestSquare вернет k столбцов, дающих k начальные значения, которые можно использовать в k прогонах optim, выбирая лучшие.

Если начальные значения достаточно хороши, это должно привести к оптимуму.

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")

давая:

> 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 Другой вариант — повторить (1), но вместо optim использовать GenSA из пакета GenSA. Он не требует начальных значений (хотя вы можете указать начальное значение с помощью аргумента par, и это может улучшить решение в некоторых случаях), поэтому код значительно короче, но, поскольку он использует имитацию отжига, можно ожидать, что его выполнение займет значительно больше времени. . Используя f (а также nr и Ainf, которые использует f) из (1). Ниже мы попробуем это без начального значения.

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

давая:

> 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

У меня есть решение, но оно не очень хорошо масштабируется:

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()

Он основан на идее, что если матрица имеет достаточно малую плотность NA, то, выполнив поиск сначала самых больших подматриц, вы, скорее всего, найдете решение довольно быстро.

Алгоритм работает путем вычисления декартова произведения всех количеств строк и количеств столбцов, которые могут быть проиндексированы из исходной матрицы для получения подматрицы. Затем набор пар отсчетов упорядочивается по убыванию размера подматрицы, которая будет создана каждой парой отсчетов; другими словами, упорядоченный произведением двух счетов. Затем он перебирает эти пары. Для каждой пары он вычисляет все комбинации индексов строк и индексов столбцов, которые могут быть взяты для этой пары счетчиков, и пробует каждую комбинацию по очереди, пока не найдет подматрицу, содержащую нулевые NA. Найдя такую ​​подматрицу, он возвращает этот набор индексов строк и столбцов в виде списка.

Результат гарантированно будет правильным, потому что он пробует размеры подматриц в порядке убывания, поэтому первая найденная должна быть самой большой (или привязанной к самой большой) возможной подматрице, которая удовлетворяет условию.

## 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

Мы видим, что функция работает очень быстро на примере матрицы OP и возвращает правильный результат.

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()

Я написал эту функцию, чтобы упростить тестирование. Мы можем вызвать его для проверки случайной входной матрицы размера NR на NC с вероятностью выбора NA в любой заданной ячейке probNA.

Вот несколько тривиальных тестов:

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

Я не знаю простого способа проверить правильность приведенного выше результата, но мне он кажется хорошим. Но для получения этого результата потребовалось почти 9 секунд. Запуск функции на матрицах среднего размера, особенно на матрице 77x132, вероятно, является безнадежным делом.

Ожидание, может ли кто-нибудь придумать блестящее эффективное решение...