寻找具有非零元素的最大大小子矩阵

Question

我有一个矩阵，其中行是基因名称，列是样本名称。该矩阵的元素为0或1（1表示该样本中的基因表达，而0表示未表达）。我想找到在最大子集样本中表达的基因的最大子集。换句话说，我想通过重新排列行和列来找到只有1的最大矩阵。

mat <- matrix(c(1,0,1,1,1,0,1,0,1),nrow = 3,byrow = T)
mat
      [,1] [,2] [,3]
[1,]    1    0    1
[2,]    1    1    0
[3,]    1    0    1
###first swap column2 and column3
mat1 <- mat
mat1[,2] <- mat[,3]
mat1[,3] <- mat[,2]
mat1
      [,1] [,2] [,3]
[1,]    1    1    0
[2,]    1    0    1
[3,]    1    1    0
###then swap row2 and row3
mat2 <- mat1
mat2[2,] <- mat1[3,]
mat2[3,] <- mat1[2,]
mat2
      [,1] [,2] [,3]
[1,]    1    1    0
[2,]    1    1    0
[3,]    1    0    1
###then the up-left is wanted matrix

- wutao

---

这个回答解决了你的问题吗？非方阵中的最大团问题 - undefined

Answer 1

问题

我理解我们被给定了一个矩阵M，其中有S行，对应于样本，G列，对应于基因，这样矩阵的每个元素m_sg等于1，如果基因g在样本s中，否则等于0，对于所有的s = 0,1,...,S-1和g = 0,1,...,G-1。

我们希望找出M的行和列的集合，形成一个全为1的矩阵，使得该矩阵的大小最大，大小通常定义为行数和列数的乘积。

解决方法

我建议解决S个问题，其中对于每个n = 0,1,...,S-1，选择的样本数固定为n。对于每个固定的n，最大化全为1的结果矩阵的大小(t)等价于最大化选择的列数(基因数)(c)，因为t = n*c。然后选择使得t最大的解作为最终解。

优化问题

一个线性规划模型是一个问题，其中变量被赋予非负实数值，以便以最大化（或最小化）给定线性函数的方式进行分配，同时受到不等式约束的限制。每个约束都被表达为变量的给定线性函数的值必须小于或等于给定值。附加约束可以将“小于或等于”替换为“等于”。

整数线性规划模型是一个线性规划模型，其附加要求是某些或所有变量受限于取整数值。

一种类型的整数线性规划模型是某些或所有变量受限于取零或一的值。这样的变量被称为二进制变量。二进制变量通常用于强制逻辑要求。下面我提议使用二进制变量解决一系列线性规划问题。

优化软件包（一些可以在线访问，一些免费）可用于解决所有这些类型的线性问题。然而，要求某些或所有变量受限于整数值使得问题变得更加困难和耗时。一个非常大的线性规划模型可能在不到一秒的时间内解决，而具有相同规模且所有变量都是二进制的线性规划问题可能需要数年时间来解决。

优化模型

如前所述，n以下等于一个常数，等于要选择的样本数。对于每个n = 0,1,...,S-1，将优化一个模型。

首先定义二进制变量：

• x_s = 1，如果选择样本s，否则为0，对于所有s = 0,1,...,S-1
• y_g = 1，如果选择基因g，否则为0，对于所有g = 0,1,...,G-1

现在陈述具有二进制变量的线性规划模型。

最大化y_g，对于所有g = 0,1,...,G-1

满足以下条件：

对于所有s = 0,1,...,S-1，x_s = 0或1

对于所有g = 0,1,...,G-1，y_g = 0或1

对于所有s = 0,1,...,S-1，x_s求和 ≤ n

x_s + y_g <= 1 对于所有满足 m_sg = 0 的 s 和 g

注释

最后一组约束条件防止选择与样本 s 对应的行（x_s = 1）和选择与基因 g 对应的列（y_g = 1），当 m_sg = 0。对于这样的一对 s 和 g，要么选择样本 s 而不选择基因 g，要么选择基因 g 而不选择样本 s，要么两者都不选择。

---

追求最优性应该强制要求将x_s在s = 0,1,...,S-1的范围内求和等于n，但如果在最优解中左侧小于n，那并不是一个问题。我将约束条件表达为不等式，希望能减少解决方案的时间。

---

我觉得有义务解决一系列问题，每个问题都有固定的物种数量，原因是结果矩阵的大小等于其行数和列数的乘积，这是模型变量的非线性函数。举个例子，假设我引入了（线性）约束条件：

nrows = y_s，其中s = 0,1,...,S-1

ncols = y_g，其中g = 0,1,...,G-1

目标函数将变为：

max nrows*ncols

这是非线性的，破坏了模型的线性结构。

---

解决优化问题可能不方便，但是对所有样本和基因的蛮力枚举可能不可行，而且除此之外的任何方法都会成为一种启发式（或者临时）算法，可能产生最优解或者接近最优解，也可能不会。

一些启发式算法可以产生界限，给出启发式解与最优解之间的最大差异。那些没有产生这种界限（或者没有产生“紧密”界限）的启发式算法的问题在于，我们永远不知道它们是否产生了一个“好”的解决方案还是一个“坏”的解决方案。

- Cary Swoveland

Answer 2

看，这是如何概括的。

mat[order(matrixStats::rowSums2(matrixStats::rowCumsums(mat))), 
    order(-matrixStats::colSums2(matrixStats::colCumsums(mat)))]
#      [,1] [,2] [,3]
# [1,]    1    1    0
# [2,]    1    1    0
# [3,]    1    0    1

然而，在下面的矩阵中，似乎更多是负的行和。

set.seed(42)
mat <- matrix(rbinom(25, 1, .5), 5, 5)
mat[order(-matrixStats::rowSums2(matrixStats::rowCumsums(mat))), 
    order(-matrixStats::colSums2(matrixStats::colCumsums(mat)))]
#      [,1] [,2] [,3] [,4] [,5]
# [1,]    1    1    1    0    1
# [2,]    1    1    1    1    0
# [3,]    1    1    1    0    0
# [4,]    1    1    0    1    0
# [5,]    0    0    0    1    1

但这可能只是一个开始。

如果还没有方法存在，你可能无法没有数学证明。

- jay.sf

---

这是一个有趣的尝试，通过对行或列进行累加来将“1”聚集成一个子矩阵，加一！ - undefined

Answer 3

从这个问答中：

找到由1组成的最大子矩阵，经过行/列的排列等同于找到最大的双全图。

library(igraph)

fMaxSubMat <- function(m) {
  cl <- max_cliques(complementer(graph_from_incidence_matrix(1 - m)))
  n <- vapply(cl, \(x) (s <- sum(x > nrow(m)))*(length(x) - s), 0)
  lapply(
    cl[n == max(n)], \(x) {
      r <- which(x <= nrow(m))
      m[x[r], x[-r] - nrow(m)]
    }
  )
}

测试OP的示例矩阵：

(mat <- matrix(c(1,0,1,1,1,0,1,0,1), 3, 3, 1, list(letters[1:3], LETTERS[1:3])))
#>   A B C
#> a 1 0 1
#> b 1 1 0
#> c 1 0 1

fMaxSubMat(mat)
#> [[1]]
#>   A C
#> a 1 1
#> c 1 1

测试一个更大的矩阵：

set.seed(1222560709)
(mat <- matrix(runif(12*20)%/%0.5, 12, 20, 1, list(letters[1:12], LETTERS[1:20])))
#>   A B C D E F G H I J K L M N O P Q R S T
#> a 1 1 0 0 0 1 0 0 0 0 1 0 1 1 0 0 1 1 1 1
#> b 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 0 1 1 0 1
#> c 1 0 0 0 0 0 0 1 1 1 0 1 1 0 1 0 1 1 0 1
#> d 1 1 1 0 1 1 0 0 1 1 0 0 0 0 0 1 1 0 0 0
#> e 1 0 1 0 0 0 0 1 0 0 1 0 0 1 1 1 1 1 1 0
#> f 0 1 0 0 0 1 0 1 0 1 1 1 1 0 0 0 0 1 0 1
#> g 0 1 1 0 1 0 0 1 0 1 1 1 1 1 0 0 1 0 0 0
#> h 0 1 1 1 0 0 0 0 1 1 0 1 1 1 0 0 0 1 1 1
#> i 1 0 1 1 0 0 1 0 0 0 1 1 1 0 0 0 1 0 1 0
#> j 0 0 0 1 1 0 1 1 1 0 1 0 1 0 1 1 0 1 1 1
#> k 1 1 0 1 1 0 0 1 0 0 1 1 0 0 0 1 1 1 1 0
#> l 1 1 0 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 1 1
fMaxSubMat(mat)
#> [[1]]
#>   R T M L J
#> b 1 1 1 1 1
#> h 1 1 1 1 1
#> c 1 1 1 1 1
#> f 1 1 1 1 1
#> 
#> [[2]]
#>   R T M B J L D C I N
#> b 1 1 1 1 1 1 1 1 1 1
#> h 1 1 1 1 1 1 1 1 1 1
#> 
#> [[3]]
#>   H R L J M T A I Q O
#> b 1 1 1 1 1 1 1 1 1 1
#> c 1 1 1 1 1 1 1 1 1 1

- jblood94

Answer 4

这是一种暴力方法，可以找到所有最大尺寸的子矩阵（正方形或非正方形）。它从矩阵的较小维度开始，例如行或列，即min(dim(mat))，并查看所有组合，以查看它们是否可以生成所需的全为1的子矩阵。 警告：如果min(dim(mat))很大，这种方法效率较低，因为组合的数量会迅速增加。

代码

f <- function(mat) {
    dimnames(mat) <- list(
        paste0("r", 1:nrow(mat)),
        paste0("c", 1:ncol(mat))
    )

    tf <- FALSE
    if (nrow(mat) < ncol(mat)) {
        mat <- t(mat)
        tf <- !tf
    }

    lst <- Filter(
        length,
        unlist(
            lapply(
                1:ncol(mat),
                \(m) {
                    combn(ncol(mat), m, \(k) {
                        v <- mat[, k, drop = FALSE]
                        u <- v[rowSums(v) == length(k), , drop = FALSE]
                        ifelse(tf, t, I)(u)
                    }, simplify = FALSE)
                }
            ),
            recursive = FALSE
        )
    )

    len <- lengths(lst)
    lst[len == max(len)]
}

输出

给定一个例子

> set.seed(0)

> (mat <- matrix(rbinom(80, 1, .5), 8, 10))
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]    1    1    0    0    1    0    0    0    0     1
[2,]    0    1    1    0    0    1    1    0    1     0
[3,]    0    0    1    0    0    1    1    1    0     0
[4,]    1    0    0    0    1    1    0    1    0     0
[5,]    1    0    1    0    1    1    1    0    1     1
[6,]    0    1    1    1    1    1    0    1    0     1
[7,]    1    0    0    0    0    1    0    0    1     0
[8,]    1    1    1    0    1    0    0    0    0     1

我们可以获得

> f(mat)
[[1]]
   c1 c5 c10
r1  1  1   1
r5  1  1   1
r8  1  1   1

[[2]]
   c2 c5 c10
r1  1  1   1
r6  1  1   1
r8  1  1   1

[[3]]
   c3 c6 c7
r2  1  1  1
r3  1  1  1
r5  1  1  1

[[4]]
   c3 c5 c10
r5  1  1   1
r6  1  1   1
r8  1  1   1

在行/列名称中，r*和c*分别表示行/列索引。

- ThomasIsCoding