【问题标题】:Subsetting non-NA子集非 NA
【发布时间】:2016-07-27 14:46:27
【问题描述】:

我有一个矩阵,其中每一行至少有一个 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。

作为菜鸟,我认为没有明显的方法可以做到这一点。谁能帮我出点主意?

【问题讨论】:

  • A[c(2,4,5),3:5] 不是最好的解决方案,有 9 个单元格吗?
  • 对于 77x132 的矩阵,您正在考虑大约 2^(77+132) ~ 8.2E62 个可能的子矩阵。我很想知道如何解决这个问题......
  • @bgoldst 或就此而言A[2:4, c(2, 4, 5)]
  • @Frank 我怀疑我们可以通过首先识别所有NAs 来显着降低维度......但除此之外也感到困惑
  • 如果您只允许连续矩阵,则更易于管理。允许行或列跳过的问题要困难得多

标签: r matrix subset na


【解决方案1】:

1) 优化 在这种方法中,我们将问题放松为一个连续优化问题,我们使用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 参数提供起始值,这在某些情况下可能会改进解决方案)因此代码要短得多,但由于它使用模拟退火,因此预计需要更长的时间跑步。使用(1)中的f(以及f 使用的nrAinf)。下面我们试试不带起始值。

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

【讨论】:

  • 您能详细说明seriate 的作用吗?帮助文件对我来说太行话了
  • 它根据method 参数排列输出排列的行和列。我们使用默认值。您可以使用该参数尝试不同的方法。它运行得非常快,但不一定会给你想要的东西,所以你必须玩一下它。这更像是一个起点,而不是一个完成的解决方案,尽管它似乎确实适用于问题中的小问题。
  • 谢谢!这对于示例矩阵和一些测试矩阵非常有效且非常快,但并非全部。使用我的实际矩阵和许多随机测试矩阵,它返回一个单元格。我将在下一条评论中举例说明。
  • A &lt;- structure( c( NA, NA, 4, 1, 4, NA, 2, NA, 3, 9, 9, NA, NA, 9, NA, 8, 1, 8, 9, 6, NA, 3, 3, 3, 1, 9, NA, NA, 8, 8, 1, 5, 8, 6, 9, 0, 3, 5, NA, 3, NA, 3, 6, 2, 5, NA, NA, 8, 6, NA, 6, 3, NA, 5, NA, NA, NA, 2, 4, 5, 4, NA, 4, 8, NA, 2, 2, NA, 9, 7, 1, 2, 9, 2, NA, 6, 1, 5, 2, 4, 4, 4, 7, NA, 0, NA, 5, 0, 7, NA, 9, NA, 1, 6, 9, 8, 8, 7, 9, 2 ), .Dim = c(10L, 10L), .Dimnames = list(paste0("R", 1:10), paste0("C", 1:10)) ) 结果是:[[1]] [1] "R10" [[2]] [1] "C5"
  • 已完全修改答案。
【解决方案2】:

我有一个解决方案,但不能很好地扩展:

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

我写了上面的函数是为了让测试更容易。我们可以调用它来测试NC 大小为NR 的随机输入矩阵,在probNA 的任何给定单元格中选择NA 的概率。


这里有一些简单的测试:

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 矩阵,可能会失败。

正在等待看看是否有人能想出一个出色的高效解决方案......

【讨论】:

    猜你喜欢
    • 1970-01-01
    • 1970-01-01
    • 2018-09-08
    • 1970-01-01
    • 2020-03-26
    • 1970-01-01
    • 2021-12-15
    • 2022-06-22
    • 2016-09-06
    相关资源
    最近更新 更多