[R] More block diagonal matrix construction code

[Robin Hankin]

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