R:将矩阵的上三角部分转换为对称矩阵 [英] R: Convert upper triangular part of a matrix to symmetric matrix

问题描述

我在R中有矩阵的上三角部分(无对角线)，并且想从上三角部分(在对角线上有1，但可以稍后进行调整)生成对称矩阵.我通常是这样的:

I have the upper triangular part of matrix in R (without diagonal) and want to generate a symmetric matrix from the upper triangular part (with 1 on the diagonal but that can be adjusted later). I usually do that like this:

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?

This works fine but now I want to work with very large matrices. Thus, I would want to avoid having to store res.upper and res (filled with 0) at the same time. Is there any way I can directly convert res.upper to a symmetric matrix without having to initialize the matrix res first?

推荐答案

我认为这里有两个问题.

I think there are two issues here.

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

now I want to work with very large matrices

然后不要使用R代码执行此工作. R将使用比您预期更多的内存.尝试以下代码:

Then do not use R code to do this job. R will use much more memory than you expect. Try the following code:

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中覆盖对象.

In R, you have to use 5 * (100 * 100) + 4950 double words to finish these operations. While in C, you only need at most 4950 + 100 * 100 double words (In fact, 100 * 100 is all that is needed! Will talk about it later). It is difficult to overwrite object directly in R without extra memory assignment.

有什么方法可以直接将res.upper转换为对称矩阵，而无需先初始化矩阵res?

Is there any way I can directly convert res.upper to a symmetric matrix without having to initialize the matrix res first?

您确实必须为res分配内存，因为最终您将获得此内存.但是无需为res.upper分配内存. 您可以初始化上方的三角形，同时填充下方的三角形.请考虑以下模板:

You do have to allocate memory for res because that is what you end up with; but there is no need to allocate memory for res.upper. You can initialize the upper triangular, while filling in the lower triangular at the same time. Consider the following template:

#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进行寻址将导致性能下降.但是我相信您可以将乘法移到内部循环之外，而只在内部保留加法.

The code has not been optimized, as using integer multiplication j * n + i for addressing in the innermost loop will result in performance penalty. But I believe you can move multiplication outside the inner loop, and only leave addition inside.

这篇关于R:将矩阵的上三角部分转换为对称矩阵的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！