r - 子集非NA

Question

我有一个矩阵，其中每一行至少有一个 NA 单元格，每一列也至少有一个 NA 单元格。我需要的是找到这个矩阵中不包含 NA 的最大子集。

例如，对于这个矩阵 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

有两种解决方案(8 个单元格): A[c(2, 4), 2:5]和 A[2:5, 4:5] ，尽管找到一个有效的解决方案就足以满足我的目的。我的实际矩阵的尺寸是 77x132。

作为一个菜鸟，我认为没有明显的方法可以做到这一点。有人可以帮我一些想法吗？

Answer 1

1) 优化 在这种方法中，我们将问题放松为一个连续优化问题，我们用 optim 解决了这个问题。 .

目标函数是 f并且它的输入是一个 0-1 向量，其第一个 nrow(A)条目对应于行，其余条目对应于列。 f使用矩阵 Ainf源自 A通过用大的负数替换 NAs，用 1 替换非 NAs。就 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