子集非 NA
Subsetting non-NA
我有一个矩阵，其中每一行至少有一个 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.
作为菜鸟，我看不出有什么明显的方法可以做到这一点。谁能帮我出点主意？
我有一个解决方案，但它的扩展性不是很好：
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 密度，那么通过首先搜索最大的子矩阵，您可能会很快找到解决方案。
该算法的工作原理是计算所有 counts 行和 counts 列的笛卡尔积，这些乘积可以从原始索引中提取出来矩阵来产生子矩阵。然后，这组计数对按照每对计数产生的子矩阵的大小降序排列；换句话说，按两个计数的乘积排序。然后迭代这些对。对于每一对，它计算可用于该对计数的行索引和列索引的所有组合，并依次尝试每个组合，直到找到包含零个 NA 的子矩阵。找到这样的子矩阵后，它 returns 将那组行和列索引作为列表。
结果保证是正确的，因为它按递减顺序尝试子矩阵大小，所以它找到的第一个必须是满足条件的最大（或并列最大）可能的子矩阵。
## 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 的示例矩阵上运行得非常快，returns 一个正确的结果。
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 的随机输入矩阵，在 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 矩阵，可能是一个失败的原因。
等着看是否有人能想出一个绝妙的有效解决方案...
1) optim 在这种方法中，我们将问题放松为一个连续优化问题，我们用 optim 解决了这个问题。
objective 函数是 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 将 return 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 使用的 nr 和 Ainf）。下面我们试试不设起始值。
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
我有一个矩阵，其中每一行至少有一个 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.
作为菜鸟，我看不出有什么明显的方法可以做到这一点。谁能帮我出点主意？
我有一个解决方案，但它的扩展性不是很好：
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 密度，那么通过首先搜索最大的子矩阵，您可能会很快找到解决方案。
该算法的工作原理是计算所有 counts 行和 counts 列的笛卡尔积，这些乘积可以从原始索引中提取出来矩阵来产生子矩阵。然后，这组计数对按照每对计数产生的子矩阵的大小降序排列；换句话说，按两个计数的乘积排序。然后迭代这些对。对于每一对，它计算可用于该对计数的行索引和列索引的所有组合，并依次尝试每个组合，直到找到包含零个 NA 的子矩阵。找到这样的子矩阵后，它 returns 将那组行和列索引作为列表。
结果保证是正确的，因为它按递减顺序尝试子矩阵大小，所以它找到的第一个必须是满足条件的最大（或并列最大）可能的子矩阵。
## 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 的示例矩阵上运行得非常快，returns 一个正确的结果。
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 的随机输入矩阵，在 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 矩阵，可能是一个失败的原因。
等着看是否有人能想出一个绝妙的有效解决方案...
1) optim 在这种方法中，我们将问题放松为一个连续优化问题，我们用 optim 解决了这个问题。
objective 函数是 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 将 return 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 使用的 nr 和 Ainf）。下面我们试试不设起始值。
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