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

Question

我在 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)]

这很好用，但现在我想处理非常大的矩阵。因此，我希望避免同时存储 res.upper 和 res（填充为 0）。有什么办法可以直接将 res.upper 转换为对称矩阵，而不必先初始化矩阵 res？

Answer 1

我认为这里有两个问题。

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

那么不要使用 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 进行寻址会导致性能下降。但我相信您可以将乘法移到内部循环之外，而只将加法留在内部。