创建对称矩阵的最有效方法

Question

I have the following matrix/ dataframe:

> e
  V1 V2 V3 V4 V5
1  0  2  3  4  5
2  0  0  6  8 10
3  0  0  0 12 15
4  0  0  0  0 20
5  0  0  0  0  0

在这种情况下，N=5（行数等于列数）。我想填补这个对称矩阵中的缺失值（例如e[1,2]=e[2,1]）。是否有一种更有效的方法来填充缺失值（在我的情况下矩阵大小很大）？是否有比嵌套循环更好的方法？

- Avi

Answer 1

为了完整起见，我想展示一下这个技巧。如果矩阵的下半部分（对角线以下）已经填入了值，只是使用转换将不起作用，因为它会将它们添加到矩阵的上半部分。

使用Matrix包，我们可以创建一个稀疏矩阵，如果创建大矩阵的对称矩阵，则需要更少的内存甚至加快速度。

为了从矩阵e创建对称稀疏矩阵，我们应该：

library(Matrix)
rowscols <- which(upper.tri(e), arr.ind=TRUE)
sparseMatrix(i=rowscols[,1],    #rows to fill in
             j=rowscols[,2],    #cols to fill in
             x=e[upper.tri(e)], #values to use (i.e. the upper values of e)
             symmetric=TRUE,    #make it symmetric
             dims=c(nrow(e),nrow(e))) #dimensions

输出：

5 x 5 sparse Matrix of class "dsCMatrix"

[1,] .  2  3  4  5
[2,] 2  .  6  8 10
[3,] 3  6  . 12 15
[4,] 4  8 12  . 20
[5,] 5 10 15 20  .

微基准测试:

让我们编写一个函数，将矩阵转换为对称矩阵（默认情况下将矩阵的上半部分复制到下半部分）：

symmetrise <- function(mat){
  rowscols <- which(upper.tri(mat), arr.ind=TRUE)
  sparseMatrix(i=rowscols[,1], 
               j=rowscols[,2], 
               x=mat[upper.tri(mat)], 
               symmetric=TRUE, 
               dims=c(nrow(mat),ncol(mat)) )  
}

And test:

> microbenchmark(
e + t(e),
symmetrise(e),
e[lower.tri(e)] <- t(e)[lower.tri(e)],
times=1000
)
Unit: microseconds
                                  expr      min       lq      mean   median        uq       max neval cld
                              e + t(e)   75.946   99.038  117.1984  110.841  134.9590   246.825  1000 a  
                         symmetrise(e) 5530.212 6246.569 6950.7681 6921.873 7034.2525 48662.989  1000   c
 e[lower.tri(e)] <- t(e)[lower.tri(e)]  261.193  322.771  430.4479  349.968  395.3815 36873.894  1000  b 

正如您所看到的，symmetrise 比 e + t(e) 或 df[lower.tri(df)] <- t(df)[lower.tri(df)] 慢得多，但至少您有一个自动对称矩阵的函数（默认情况下，它取上半部分并创建下半部分），如果您有一个内存问题的大矩阵，这可能会派上用场。

附注：在矩阵中任何地方看到 . 表示零。通过使用不同的系统，稀疏矩阵是一种“压缩”的对象，使其更节省内存。

- LyzandeR

---

你需要转置矩阵t(e)[lower.tri(e)]的上三角部分，就像我在我的答案中所做的那样，否则你将无法得到与e + t(e)相同的结果。 - mpalanco

@mpalanco 是的，你说得对。我没有注意到这一点。如果我是 Avi，我会更改已接受的答案。如果我要切换到 mpalanco 的答案，那么这里就没有这个答案的意义了。 - LyzandeR

你非常友善。我一开始也是这样做的，后来才意识到它不对称。发布基准测试结果的想法很好。谢谢。 - mpalanco

我改变了我的答案，提供了一个内存高效的解决方案，并创建了一个将矩阵转换为对称矩阵的函数。我仍然认为@mpalanco的想法对于小到中等大小的矩阵是最好的，但如果您有内存限制或大型矩阵，则这个解决方案应该非常有效。 - LyzandeR

Answer 2

另外需要考虑速度：

2*symmpart(as.matrix(e))

这是一个基准测试：

Unit: microseconds
                                      expr      min       lq        mean    median        uq       max neval
                                  e + t(e)  572.505  597.194  655.132028  611.5420  628.4860  8424.902  1000
                             symmetrise(e) 1128.220 1154.562 1215.740071 1167.0020 1185.6585 10656.059  1000
 e[lower.tri(e)] <- e[upper.tri(e, FALSE)]  285.013  311.191  350.846885  327.1335  339.5910  8106.006  1000
                2 * symmpart(as.matrix(e))   78.392   93.953  101.330522  102.1860  107.9215   153.628  1000

它之所以能达到这样的速度，是因为它直接创建对称矩阵。

- Neal Fultz

---

在Matrix包中，对于所有方阵，都有x == symmpart(x) + skewpart(x)。 - Neal Fultz

这是一个很好的答案，实际上也是最快的（我已经点赞了），但是仅在下部或上部全部为零且仅在这种情况下才有效。正如描述中所说，矩阵被计算为 (x + t(x))/2。 - LyzandeR

Answer 3

df[lower.tri(df)] <- t(df)[lower.tri(df)]

输出：

  V1 V2 V3 V4 V5
1  0  2  3  4  5
2  2  0  6  8 10
3  3  6  0 12 15
4  4  8 12  0 20
5  5 10 15 20  0

数据：

df <- structure(list(V1 = c(0L, 0L, 0L, 0L, 0L), V2 = c(2L, 0L, 0L, 
0L, 0L), V3 = c(3L, 6L, 0L, 0L, 0L), V4 = c(4L, 8L, 12L, 0L, 
0L), V5 = c(5L, 10L, 15L, 20L, 0L)), .Names = c("V1", "V2", "V3", 
"V4", "V5"), class = "data.frame", row.names = c("1", "2", "3", 
"4", "5"))

- mpalanco

Answer 4

e + t(e)

你想要的是将矩阵和它的转置相加吗？

- CarlAH

---

是的。但我想以最有效率的方式完成它（最短时间）。 - Avi

我不确定那种方法是否是最有效的方式，但我猜应该比嵌套循环要好。 - CarlAH

如果我没记错的话，这是对BLAS的直接调用...所以如果你想要更快的速度，就转向C/Julia/Fortran。 - MichaelChirico

有没有办法在R中使用BLAS来处理这种情况？ - Avi

R使用BLAS。基本搜索表明，通过更改默认的BLAS（ATLAS），您可能会获得一些改进。 - MichaelChirico

参考此处，其中实现了t.default的C源代码（作为do_transpose）。 - MichaelChirico