R help - Sep 2005 - More block diagonal matrix construction code

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

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