将矩阵的上三角部分转换为对称矩阵

Question

将矩阵的上三角部分转换为对称矩阵

8

我有 R 中矩阵的上三角部分（不包括对角线），想要从上三角部分生成一个对称矩阵（对角线上为 1，但稍后可以进行调整）。我通常会这样做：

res.upper <- rnorm(4950)
res <- matrix(0, 100, 100)
res[upper.tri(res)] <- res.upper
rm(res.upper)
diag(res) <- 1
res[lower.tri(res)]  <- t(res)[lower.tri(res)]

这个代码能正常工作，但现在我想处理非常大的矩阵。因此，我希望能避免同时存储填满0的res.upper和res矩阵。有没有办法能直接将res.upper转换为对称矩阵而不必先初始化res矩阵？

- Lila

我可以编写已编译的代码，有时也会这样做以加快我的函数速度。然而，我并不真正理解这如何避免使用额外的内存。在C/C++代码中，我也会首先初始化像res这样的对象。那么这不也会使用额外的内存吗？或者说，在使用C/C++时这不是问题，因为这些语言中的内存分配更加“智能”吗？这可能是一个愚蠢的问题，但我是一名统计学家，而不是计算机科学家，所以我真的不知道内存分配的内部工作原理。 - Lila

你不必为我编写函数，这不是问题。我熟悉编译和内联包。我只是缺乏足够的背景来理解它如何解决我的内存问题。但如果你向我保证它可以解决问题，我会编写该函数（如果你将其作为答案而不是评论给出，我会接受你的答案）。 - Lila

你尝试过使用bigmemory包中的big.matrix吗？这可能是绕过内存限制的一种方法。 - konvas

最终的矩阵必须转换为kernlab库中的kernelMatrix类。我不知道这是否适用于big.matrix，但我一定会看看！谢谢你的提示！我不知道这个包，这是我第一次处理如此大的矩阵。 - Lila

你可能还会发现 "Matrix" 包很有用 -- 即 sparseMatrix(i = sequence(1:99), j = rep(2:100, 1:99), x = res.upper, symmetric = TRUE, dims = c(100, 100))，避免了多次复制和使用（可能）更小的对象。 - alexis_laz

@alexis_laz：谢谢你的提示！我想我会采用 Zheyuan Li 提供的 C 解决方案。但我还会看一下 sparseMatrix，了解它对于未来的问题也无妨。 - Lila

2个回答

2

要从上三角矩阵或下三角矩阵得到对称矩阵，可以将该矩阵与其转置相加并减去对角线元素。以下是相关方程式。

diag(U)是一个由U的对角线元素组成的对角矩阵。

Equation 1

ultosymmetric=function(m){
  m = m + t(m) - diag(diag(m))
return (m)}

如果您希望对角线元素为1，可以这样做。

ultosymmetric_diagonalone=function(m){
  m = m + t(m) - 2*diag(diag(m)) + diag(1,nrow=dim(m)[1])
return (m)}

Equation 2

- Yehya Abdelmohsen

为了简化起见，在您的函数中不需要使用 m = 和 return(m)。R默认返回最后一个被评估的表达式。 - Martin Gal

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Zheyuan Li · Accepted Answer

我认为这里有两个问题。

现在我想要处理非常大的矩阵

那么请不要使用R代码来完成此工作。R将占用比您预期更多的内存。请尝试以下代码：

res.upper <- rnorm(4950)
res <- matrix(0, 100, 100)
tracemem(res)  ## trace memory copies of `res`
res[upper.tri(res)] <- res.upper
rm(res.upper)
diag(res) <- 1
res[lower.tri(res)]  <- t(res)[lower.tri(res)]

这是您将获得的内容：

> res.upper <- rnorm(4950)  ## allocation of length 4950 vector
> res <- matrix(0, 100, 100)  ## allocation of 100 * 100 matrix
> tracemem(res)
[1] "<0xc9e6c10>"
> res[upper.tri(res)] <- res.upper
tracemem[0xc9e6c10 -> 0xdb7bcf8]: ## allocation of 100 * 100 matrix
> rm(res.upper)
> diag(res) <- 1
tracemem[0xdb7bcf8 -> 0xdace438]: diag<-  ## allocation of 100 * 100 matrix
> res[lower.tri(res)]  <- t(res)[lower.tri(res)]
tracemem[0xdace438 -> 0xdb261d0]: ## allocation of 100 * 100 matrix
tracemem[0xdb261d0 -> 0xccc34d0]: ## allocation of 100 * 100 matrix

在R中，您需要使用5 * (100 * 100) + 4950个双字来完成这些操作。而在C中，您只需要最多4950 + 100 * 100个双字（实际上，只需要100 * 100个双字！稍后会讨论）。在R中，直接覆盖对象而不进行额外的内存分配是困难的。

有没有办法可以直接将res.upper转换为对称矩阵，而无需先初始化矩阵res？

您确实需要为res分配内存，因为这就是您最终得到的；但是，没有必要为res.upper分配内存。您可以同时初始化上三角形，并填写下三角形。考虑以下模板：

#include <Rmath.h>  // use: double rnorm(double a, double b)
#include <R.h>  // use: getRNGstate() and putRNGstate() for randomness
#include <Rinternals.h>  // SEXP data type

## N is matrix dimension, a length-1 integer vector in R
## this function returns the matrix you want
SEXP foo(SEXP N) {
  int i, j, n = asInteger(N);
  SEXP R_res = PROTECT(allocVector(REALSXP, n * n));  // allocate memory for `R_res`
  double *res = REAL(R_res);
  double tmp;  // a local variable for register reuse
  getRNGstate();
  for (i = 0; i < n; i++) {
    res[i * n + i] = 1.0;  // diagonal is 1, as you want
    for (j = i + 1; j < n; j++) {
      tmp = rnorm(0, 1);  
      res[j * n + i] = tmp; // initialize upper triangular
      res[i * n + j] = tmp;  // fill lower triangular
      }
    }
  putRNGstate();
  UNPROTECT(1);
  return R_res;
  }

这段代码还没有被优化，因为在最内层循环中使用整数乘法 j * n + i 来寻址将导致性能损失。但我相信您可以将乘法移到内层循环之外，只留下加法。