在R中高效地比较矩阵

Question

我有一个包含一些矩阵的数组a。现在我需要高效地检查我有多少个不同的矩阵以及它们在数组中的索引（按升序排列）。我的方法是：将矩阵的列粘贴为字符向量，然后查看频率表，像这样：

n <- 10 #observations
a <- array(round(rnorm(2*2*n),1),
           c(2,2,n))

paste_a <- apply(a, c(3), paste, collapse=" ") #paste by column
names(paste_a) <- 1:n
freq <- as.numeric( table(paste_a) ) # frequencies of different matrices (in ascending order)
indizes <- as.numeric(names(sort(paste_a[!duplicated(paste_a)]))) 
nr <- length(freq) #number of different matrices

然而，当您将 n 增加到很大的数字时，这会变得非常低效（主要是 paste() 变得越来越慢）。有没有更好的解决方案？这里有一个包含 100 个观测值的“真实”数据集，其中一些矩阵是实际重复的（与上面的示例不同）：https://pastebin.com/aLKaSQyF。非常感谢。

- yrx1702

---

可以把这些矩阵存储在一个列表里吗？将它们存储在列表里更容易处理和应用函数。 - AntoniosK

最好提供一些实际的重复/复制样本。另外，由于您没有提供种子，所有结果都是不同的。 - talat

Answer 1

既然您的实际数据由整数0、1、2、3组成，为什么不利用4进制呢？整数比整个矩阵对象更容易比较。 （下面所有出现的a都是从链接中找到的真实数据集中的数据。）

Base4Approach <- function() {
    toBase4 <- sapply(1:dim(a)[3], function(x) {
        v <- as.vector(a[,,x])
        pows <- which(v > 0)
        coefs <- v[pows]
        sum(coefs*(4^pows))
    })
    myDupes <- which(duplicated(toBase4))
    a[,,-(myDupes)]
}

既然问题涉及效率，让我们进行基准测试：

MartinApproach <- function() {
    ### commented this out for comparison reasons
    # dimnames(a) <- list(1:dim(a)[1], 1:dim(a)[2], 1:dim(a)[3])
    a <- a[,,!duplicated(a, MARGIN = 3)]
    nr <- dim(a)[3]
    a
}

identical(MartinApproach(), Base4Approach())
[1] TRUE

microbenchmark(Base4Approach(), MartinApproach())
Unit: microseconds
            expr     min       lq      mean    median       uq      max neval
 Base4Approach() 291.658  303.525  339.2712  325.4475  352.981  636.361   100
MartinApproach() 983.855 1000.958 1160.4955 1071.9545 1187.321 3545.495   100

@d.b.的方法实际上与前两种方法并不完全相同（它仅仅是识别重复项而不移除它们）。

DBApproach <- function() {
    a[, , 9] = a[, , 1]

    #Convert to list
    mylist = lapply(1:dim(a)[3], function(i) a[1:dim(a)[1], 1:dim(a)[2], i])
    temp = sapply(mylist, function(x) sapply(mylist, function(y) identical(x, y)))
    temp2 = unique(apply(temp, 1, function(x) sort(which(x))))

    #The indices in 'a' where the matrices are same
    temp2[lengths(temp2) > 1]
}

然而，Base4Approach 仍然占主导地位:

    microbenchmark(Base4Approach(), MartinApproach(), DBApproach())
Unit: microseconds
            expr      min         lq       mean    median         uq       max neval
 Base4Approach()  298.764   324.0555   348.8534   338.899   356.0985   476.475   100
MartinApproach() 1012.601  1087.9450  1204.1150  1110.662  1162.9985  3224.299   100
    DBApproach() 9312.902 10339.4075 11616.1644 11438.967 12413.8915 17065.494   100

由@alexis_laz提供的更新

正如@alexis_laz在评论中提到的，我们可以做得更好。

AlexisBase4Approach <- function() {
    toBase4 <- colSums(a * (4 ^ (0:(prod(dim(a)[1:2]) - 1))), dims = 2)
    myDupes <- which(duplicated(toBase4))
    a[,,-(myDupes)]
}

microbenchmark(Base4Approach(), MartinApproach(), DBApproach(), AlexisBase4Approach(), unit = "relative")
Unit: relative
                 expr       min        lq       mean     median         uq        max neval
      Base4Approach()  11.67992  10.55563   8.177654   8.537209   7.128652   5.288112   100
     MartinApproach()  39.60408  34.60546  27.930725  27.870019  23.836163  22.488989   100
         DBApproach() 378.91510 342.85570 262.396843 279.190793 231.647905 108.841199   100
AlexisBase4Approach()   1.00000   1.00000   1.000000   1.000000   1.000000   1.000000   100

## Still gives accurate results
identical(MartinApproach(), AlexisBase4Approach())
[1] TRUE

- Joseph Wood

---

我还看到一些2s离子在那里。 - Martin Schmelzer

@MartinSchmelzer，谢谢你发现了这个问题。回答已更新。 - Joseph Wood

没关系。你的方法很好。我也修改了我的答案。 - Martin Schmelzer

谢谢你们的努力，非常感激！ - yrx1702

"toBase4"也可以计算为colSums(a * (4 ^ (0:(prod(dim(a)[1:2]) - 1))), dims = 2)，除非我遗漏了什么。 - alexis_laz

我认为 toBase4 <- packBits(intToBits(a)[c(T,T,rep(F,30))], type = 'integer') 应该很快。它比你的原始函数要快，但 @alexis_laz 的版本仍然胜出。 - dww

Answer 2

我的第一次尝试实际上是非常慢的。所以这里稍微改变了你的版本：

  dimnames(a) <- list(1:dim(a)[1], 1:dim(a)[2], 1:dim(a)[3])
  a   <- a[,,!duplicated(a, MARGIN = 3)]
  nr  <- dim(a)[3] #number of different matrices
  idx <- dimnames(a)[[3]] # indices of left over matrices

- Martin Schmelzer

Answer 3

我不确定这是否完全符合你的要求，但以下是一种提取矩阵相同位置索引的方法。为得到你想要的结果可能需要进行更多处理。

#DATA
n <- 10
a <- array(round(rnorm(2*2*n),1), c(2,2,n))
a[, , 9] = a[, , 1]

temp = unique(apply(X = sapply(1:dim(a)[3], function(i)
    sapply(1:dim(a)[3], function(j) identical(a[, , i], a[, , j]))),
    MARGIN = 1,
    FUN = function(x) sort(which(x))))
temp[lengths(temp) > 1]
#[[1]]
#[1] 1 9

- d.b