R : 행렬의 위 삼각형 부분을 대칭 행렬로 변환

Question

행렬의 위 삼각형 부분이 R (대각선 없음)이고 위 삼각형 부분에서 대칭 행렬을 생성하려고합니다. 나중에 조정할 수 있음). 나는 보통 다음과 같이 이것을한다 :

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를 초기화하지 않고도 res.upper를 대칭 매트릭스로 직접 변환 할 수있는 방법이 있습니까?

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

이 당신이 얻을 것이다 무엇 : R에서

> 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 

를 사용하면 이러한 작업을 완료 5 * (100 * 100) + 4950 두 단어를 사용합니다. C에서, 당신은 오직 4950 + 100 * 100 더블 단어 (사실, 100 * 100 전부가 필요합니다! 나중에 얘기 할 것입니다). 여분의 메모리 할당없이 R에서 객체를 직접 덮어 쓰는 것은 어렵습니다.

먼저 res 행렬을 초기화하지 않고 res.upper을 대칭 행렬로 직접 변환 할 수있는 방법이 있습니까?

결국 res에 메모리를 할당해야합니다. 그러나 res.upper에 메모리를 할당 할 필요는 없습니다. 하단 삼각형을 동시에 채우면서 상단 삼각형을 초기화 할 수 있습니다. 다음 템플릿 고려해 성능이 저하 될 것이다 최 루프 어드레싱 정수 승산 j * n + i을 사용하는 것으로,

#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; 
    } 

코드가 최적화되지 않았습니다. 그러나 내부 루프 외부로 곱셈을 이동하고 내부에 추가 만 남길 수 있다고 생각합니다.