More block diagonal matrix construction code
hi Bert, List
well now seems a good time to adduce adiag() of package magic.
Function adiag binds together arrays of arbitrary dimension
corner-to-corner. Sensible interpretation is made for
arguments at the "edge" of acceptability (eg one array
being a scalar).
The "meat" of the code is as follows:
adiag <- function(a,b,pad=0){
s <- array(pad, dim.a + dim.b)
s <- do.call("[<-", c(list(s), lapply(dim.a, seq.new), list(a)))
ind <- lapply(seq(dim.b), function(i) seq.new(dim.b[[i]]) +
dim.a[[i]])
out <- do.call("[<-", c(list(s), ind, list(b)))
return(out)
}
where
seq.new <- function(i) { if (i == 0) { return(NULL)} else { return
(1:i) } }
[NB: untested].
so it creates an array "s" of the right size filled with "pad", and then
fills one corner with "a", then fills the other corner with "b".
Note the absence of any for() loops.
Hope this is useful
rksh
On 2 Sep 2005, at 00:08, Berton Gunter wrote:
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
______________________________________________ R-help at stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting- guide.html
-- Robin Hankin Uncertainty Analyst National Oceanography Centre, Southampton European Way, Southampton SO14 3ZH, UK tel 023-8059-7743