[R] More block diagonal matrix construction code
Berton Gunter
gunter.berton at gene.com
Fri Sep 2 01:08:58 CEST 2005
Folks:
In answer to a query, Andy Liaw recently submitted some code to construct a
block diagonal matrix. For what seemed a fairly straightforward task, the
code seemed a little "overweight" to me (that's an American stock analyst's
term, btw), so I came up with a slightly cleaner version (with help from
Andy):
bdiag<-function(...){
mlist<-list(...)
## handle case in which list of matrices is given
if(length(mlist)==1)mlist<-unlist(mlist,rec=FALSE)
csdim<-rbind(c(0,0),apply(sapply(mlist,dim),1,cumsum ))
ret<-array(0,dim=csdim[length(mlist)+1,])
add1<-matrix(rep(1:0,2),nc=2)
for(i in seq(along=mlist)){
indx<-apply(csdim[i:(i+1),]+add1,2,function(x)x[1]:x[2])
## non-square matrix
if(is.null(dim(indx)))ret[indx[[1]],indx[[2]]]<-mlist[[i]]
## square matrix
else ret[indx[,1],indx[,2]]<-mlist[[i]]
}
ret
}
I doubt that there's any noticeable practical performance difference, of
course.
The strategy is entirely basic: just get the right indices for replacement
of the arguments into a matrix of 0's of the right dimensions. About the
only thing to notice is that the apply() construction returns either a list
or matrix depending on whether a matrix is square or not (a subtlety that
tripped me up in my first version of this code).
I would be pleased if this stimulated others to come up with cleverer/more
elegant approaches that they would share, as it's the sort of thing that
I'll learn from and find useful.
Cheers to all,
Bert Gunter
More information about the R-help
mailing list