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

【问题标题】：Most Efficient way to create a symmetric matrix创建对称矩阵的最有效方法
【发布时间】：2015-12-22 13:37:08
【问题描述】：

我有以下矩阵/数据框：

> 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] 等）。是否有一种最有效的方法来填充缺失值（在我的情况下，矩阵大小 N 很大）？有没有比嵌套循环更好的方法？

【问题讨论】：

标签： r loops

【解决方案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)) )  
}

然后测试：

> 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)] 慢得多，但至少您有一个自动对称矩阵的函数（默认情况下取上半部分并创建下半部分），以防万一有一个很大的矩阵，内存是个问题，这可能会派上用场。

附：无论您在矩阵中的何处看到.，这都表示零。通过使用不同的系统，稀疏矩阵是一种“压缩”对象，使其内存效率更高。

【讨论】：

您需要转置矩阵t(e)[lower.tri(e)] 的上三角部分，就像我在回答中所做的那样，否则您将不会得到与e + t(e) 相同的结果
@mpalanco 是的，你是对的。我没有注意到这一点。 Avi 如果我是你，我会改变接受的答案。如果我要切换到 mpalanco 的，这里没有这个答案
你很善良。一开始我也是这样做的，但我意识到它不是对称的。发布基准的好主意。谢谢。
我将答案更改为内存高效解决方案，并创建了一个将矩阵转换为对称矩阵的函数。我仍然认为@mpalanco 的想法最适合中小型矩阵，但如果您有内存限制或大型矩阵，这应该会很好。

【解决方案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

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

【讨论】：

在Matrix 包中，x == symmpart(x) + skewpart(x) 用于所有方阵
这是一个很好的答案，实际上是最快的（我赞成），但这仅适用于下部或上部全部为零的情况，只有那时。正如描述中所说，矩阵计算为(x + t(x))/2。

【解决方案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"))

【讨论】：

【解决方案4】：

e + t(e)

添加矩阵和该矩阵的转置，这是你想要的吗？

【讨论】：

是的。但我想以最有效的方式（最短时间）来做。
我不确定这种方式是否是最有效的方式，但它应该比我猜的嵌套循环更好？
如果我没记错的话，这是对 BLAS 的直接吸引力......所以如果你想要更快的速度，请转到 C/Julia/Fortran
R 使用 BLAS。 basic search 建议您通过更改默认 BLAS may get some improvement ATLAS
供参考，here 是在 C 中实现t.default 的源代码（如do_transpose）。