[R] block diagonal matrix function

[Dimitris Rizopoulos]

Dear Robin,

you could use the following function which has been submitted to
s-news some time ago:

bdiag <- function(x){
     if(!is.list(x)) stop("x not a list")
     n <- length(x)
     if(n==0) return(NULL)
     x <- lapply(x, function(y) if(length(y)) as.matrix(y) else
stop("Zero-length component in x"))
     d <- array(unlist(lapply(x, dim)), c(2, n))
     rr <- d[1,]
     cc <- d[2,]
     rsum <- sum(rr)
     csum <- sum(cc)
     out <- array(0, c(rsum, csum))
     ind <- array(0, c(4, n))
     rcum <- cumsum(rr)
     ccum <- cumsum(cc)
     ind[1,-1] <- rcum[-n]
     ind[2,] <- rcum
     ind[3,-1] <- ccum[-n]
     ind[4,] <- ccum
     imat <- array(1:(rsum * csum), c(rsum, csum))
     iuse <- apply(ind, 2, function(y, imat) imat[(y[1]+1):y[2],
(y[3]+1):y[4]], imat=imat)
     iuse <- as.vector(unlist(iuse))
     out[iuse] <- unlist(x)
     return(out)
}

#########################

mats <- list(matrix(1:20, 5, 4), matrix(1:12, 4, 3), matrix(1:25, 5,
5))
bdiag(mats)

Regarding the dimnames of the matrices, I think you can easily create
a function thas calls bdiag and gives the appropriate names to the
final matrix.

I hope this helps.

Best,
Dimitris

----
Dimitris Rizopoulos
Doctoral Student
Biostatistical Centre
School of Public Health
Catholic University of Leuven

Address: Kapucijnenvoer 35, Leuven, Belgium
Tel: +32/16/396887
Fax: +32/16/337015
Web: http://www.med.kuleuven.ac.be/biostat/
     http://www.student.kuleuven.ac.be/~m0390867/dimitris.htm


----- Original Message ----- 
From: "Robin Hankin" <rksh at soc.soton.ac.uk>
To: <r-help at stat.math.ethz.ch>
Sent: Thursday, May 27, 2004 10:24 AM
Subject: [R] block diagonal matrix function


> Hello List
>
> I have just written a little function that takes two matrices as
> arguments and returns a large matrix that is composed of the two
input
> matrices in upper-left position and lower-right position with a
padding
> value everywhere else. (function definition and toy example below).
I
> need nonsquare matrices and rowname() and colname() inherited
appropriately.
>
> Two questions:
>
> (1) Is there a better way to do this? (kronecker() isn't applicable
here)
>
> (2) How do I generalize it to take an arbitrary number of matrices
as
>      inputs?
>
>
> TIV
>
> Robin
>
>
>
> "blockdiag" <-
> function (m1, m2, p.tr = 0, p.ll = 0)
> {
>      ## p.tr and p.ll are padding values
>      topleft <- m1
>      topright <- matrix(p.tr, nrow(m1), ncol(m2))
>      colnames(topright) <- colnames(m2)
>      lowleft <- matrix(p.ll, nrow(m2), ncol(m1))
>      lowright <- m2
>      rbind(cbind(topleft, topright), cbind(lowleft, lowright))
> }
>
> m1 <-
> structure(c(1, 1, 3, 1, 3, 4), .Dim = as.integer(c(2, 3)), .Dimnames
= list(
>      c("a", "b"), c("x", "y", "z")))
>
> m2 <-
> structure(c(2, 1, 1, 0, 3, 2, 2, 2, 0), .Dim = as.integer(c(3,
> 3)), .Dimnames = list(c("I", "II", "III"), c("A", "B", "C")))
>
> R> m1
>    x y z
> a 1 3 3
> b 1 1 4
>
> R> m2
>      A B C
> I   2 0 2
> II  1 3 2
> III 1 2 0
>
> R> blockdiag(m1,m2)
>      x y z A B C
> a   1 3 3 0 0 0
> b   1 1 4 0 0 0
> I   0 0 0 2 0 2
> II  0 0 0 1 3 2
> III 0 0 0 1 2 0
> R>
> -- 
> Robin Hankin
> Uncertainty Analyst
> Southampton Oceanography Centre
> SO14 3ZH
> tel +44(0)23-8059-7743
> initialDOTsurname at soc.soton.ac.uk (edit in obvious way; spam
precaution)
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html