• R语言
  • 三维数组操作---矩阵运算

我有个MCMC的程序,中间有一步涉及到加权最小二乘。需要计算如下的量:

β^=(i=1nXiR1Xi)1(i=1nXiR1Zi)\widehat{ \beta} = ( \sum_{i=1}^n X_i^\top R^{-1} X_i )^{-1} ( \sum_{i=1}^n X_i^\top R^{-1} Z_i)

其中,XiX_i 是一个T×p T\times p 的矩阵,RR 是一个 T×T T\times T 的矩阵, Zi Z_i 是一个 T×1T\times 1 的向量。
问题在于,我现在将 Xi X_i 堆起来放在一个 N×T×pN \times T\times p 的数组 X 里面,ZiZ_i 堆起来放在 N×TN\times T 的矩阵 ZZ里面。
那么每次计算,都需要做一个循环。样本量稍微大一点,计算的代价就太大。尤其是在mcmc中,每次迭代都要涉及这样一个显式的循环。
能否直接化成矩阵或者数组的操作而避免显式的循环呢?
多谢!

Compared to the use of an explicit loop, my guess is that some other problems deserve more attention in terms of speed. Examing profiling results might help.

If you want to keep the loop, then one thing to consider is to reorganize your NxTxP matrix to a TxPxN matrix, if your loop iterates over the index i. This will allow consecutive memory access and might make better use of cache.

If you want to remove the loop, then an NT x P matrix is possible, esp. useful when R is diagonal.

Is R diagonal? Will R change over MCMC? Will Z change over MCMC?

回复 第1楼 的 dingpeng:
对于第一个部分可以这么写(至少看是看不出循环的)

mat<-list(X1=matrix(1:4,2),X2=matrix(5:8,2))
R<-matrix(rnorm(4),2)
f<-function(i) {mat[[i]]%*%R%*%mat[[i]]}
LLL=1:2###2 matrix
part1<-solve(Reduce('+',lapply(LLL,f)));
#or use snow package:
cl <- makeCluster(2, type = "MPI")
mat<-list(X1=matrix(1:4,2),X2=matrix(5:8,2))
R<-matrix(rnorm(4),2)
f<-function(i,mat=mat,R=R) {mat[[i]]%*%R%*%mat[[i]]}
LLL=1:2###2 matrix
Reduce('+',parLapply(cl,LLL,f,mat=mat,R=R))
stopCluster(cl)

如果 R\mathbf{R}有一定规模,时间估计主要会消耗在矩阵求逆和相乘上吧,换掉默认的 BLAS 提升一下这部分的效率,也许会比改写显式循环收效更明显。

回复 第4楼 的 nan.xiao:R可以很大,但是一般来说,我考虑小型的R,比如维数小于5*5。复杂度,我觉得还是主要在循环上。

回复 第2楼 的 micro@: thanks for you comments. R will change during the mcmc iterations and R is not diagonal, because it is a general correlation matrix.

i will think about your suggestions, which are very helpful.

回复 第3楼 的 enthumelon:谢谢!我去试试看。

7 年 后
Cloud2016 更改标题为「三维数组操作---矩阵运算