2016-03-31 13 views
1

を帯状より効率的な方法は、私は次のコードを使用マトリックス

 [,1] [,2] [,3]  [,4] [,5] [,6] [,7] [,8] [,9] 
[1,] 4 0.8 0.4 0.2666667 0.3 0.0 0.0 0.0 0.00 
[2,] 0 4.0 0.8 0.4000000 0.4 0.3 0.0 0.0 0.00 
[3,] 0 0.0 4.0 0.8000000 0.6 0.4 0.3 0.0 0.00 
[4,] 0 0.0 0.0 4.0000000 1.2 0.6 0.4 0.3 0.00 
[5,] 0 0.0 0.0 0.0000000 9.0 1.8 0.9 0.6 0.45 
[6,] 0 0.0 0.0 0.0000000 0.0 9.0 1.8 0.9 0.60 
[7,] 0 0.0 0.0 0.0000000 0.0 0.0 9.0 1.8 0.90 
[8,] 0 0.0 0.0 0.0000000 0.0 0.0 0.0 9.0 1.80 
[9,] 0 0.0 0.0 0.0000000 0.0 0.0 0.0 0.0 9.00 

しかし、このコードはケースを計算するには余りにも遅いです大きなnの効率的なソリューションはありますか?

答えて

1

私はあなただけではなく引き算の任意の機能を可能にdiffのような第1のヘルパー機能、(5

の帯域幅を持つスパース帯行列をしたいようにあなたのnは、非常に大きいと仮定するつもりです-)。この場合

fdiff <- function(x,lag,f) { 
    i1 <- -seq_len(lag) 
    f(x[i1],x[-length(x):-(length(x)-lag+1L)]) 
} 

は、私たちが望む機能は、最初の対角がスペアバンド行列を取り込むために、我々は「bandSparse」を使用

x <- c(rep(4,4),rep(9,5)) 
0.2*fdiff(x,1,gm)/1 
# [1] 0.8 0.8 0.8 1.2 1.8 1.8 1.8 1.8 

によって与えられること

gm <- function(x,y) sqrt(x*y) 

そうですMatrixライブラリ

library(Matrix) 
x <- c(rep(4,4),rep(9,5)) 
bandSparse(n,k=0:4,diagonals= 
    c(list(x),lapply(1:4,function(lag) 0.2*fdiff(x,lag,gm)/lag))) 

出力:あなたが使用してどのような

9 x 9 sparse Matrix of class "dgCMatrix" 

[1,] 4 0.8 0.4 0.2666667 0.3 . . . . 
[2,] . 4.0 0.8 0.4000000 0.4 0.3 . . . 
[3,] . . 4.0 0.8000000 0.6 0.4 0.3 . . 
[4,] . . . 4.0000000 1.2 0.6 0.4 0.3 . 
[5,] . . . .   9.0 1.8 0.9 0.6 0.45 
[6,] . . . .   . 9.0 1.8 0.9 0.60 
[7,] . . . .   . . 9.0 1.8 0.90 
[8,] . . . .   . . . 9.0 1.80 
[9,] . . . .   . . . . 9.00 
0

、あなたは行列の上/下三角を入力する必要があります。この場合、パターンは非常に単純であり、すなわちmat1 <- rho/(col(mat)-row(mat));diag(mat1)=mat;#Cov matである。相関行列を計算するには、列にRを格納する要素に気づき、v=sqrt(diag(mat1));mat1=mat1/v;mat1=t(t(mat1)/v);を実行します。これを1行にまとめると他の人のコメントにmat1=mat1/vコピーを避けることができます。

しかし、各エントリのパターンがはるかに複雑な場合は、inline::cxxfunctionの使用を検討することができます。元のRコードをほとんど変更する必要はありません。メモリコピーを避けるために、行列のインデックス(実際には二重ベクトル)も使用できます。

library(inline) 
include=" 
#include <math.h> 
#include <vector> 
" 

body=" 
NumericMatrix x(X); 
int nrow = x.nrow(); 
int ncol = x.ncol(); 
std::vector<double> diag(nrow); 
for (int i=0;i<nrow;i++){ 
    diag[i] = sqrt(x(i,i)); 
} 
double rho = .2; 
for(int j=1;j<ncol;j++){ 
     for(int i=0; i<(j-1);i++){ 
       if (j < i + nrow-4){// Change to your version 
        x(i,j) = rho/double(j-i)*diag[i]*diag[j]; 
        x(j,i) = x(i,j); 
} 
    } 
} 
return(x); 
" 
f1 <- cxxfunction(signature(X='matrix'),body,plugin='Rcpp',include=include) 

使用法:

> dim(C) 
[1] 3000 3000 
> system.time(f1(C)) 
    user system elapsed 
    0.086 0.000 0.087 

は、したがって、コンパイル時間を無視し、それも私のクロームブックにかなり速いです。