[R] multiple dimensional diag()

[Peter Wolf]

Here is a function that does what you want (perhaps):

<<*>>=
a.block.diag <- function(a,b) {
  # a.block.daig combines arrays a and b and builds a new array which has
  # a and b as blocks on its diagonal. pw 10/2004
  if(length(dim.a<-dim(a))!= length(dim.b<-dim(b))){
    stop("a and b must have identical number of dimensions")}
  s<-array(0,dim.a+dim.b)
  s<-do.call("[<-",c(list(s),lapply(dim.a,seq),list(a)))
  ind<-lapply(seq(dim.b),function(i)seq(dim.b[[i]])+dim.a[[i]])
  do.call("[<-",c(list(s),ind,list(b)))
}

@
Try:
<<*>>=
a=matrix(1,3,4); b=matrix(2,2,2)
a.block.diag(a,b)

@
output-start
Fri Oct  1 17:20:26 2004
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    1    1    1    0    0
[2,]    1    1    1    1    0    0
[3,]    1    1    1    1    0    0
[4,]    0    0    0    0    2    2
[5,]    0    0    0    0    2    2
output-end

and an another example:

<<*>>=
xx<-array(1:8,c(2,rep(2,2))); yy<-array(-1*(1:9),c(2,rep(3,2)))
a.block.diag(xx,yy)

@
output-start
Fri Oct  1 17:21:38 2004
, , 1
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    3    0    0    0
[2,]    2    4    0    0    0
[3,]    0    0    0    0    0
[4,]    0    0    0    0    0
, , 2
     [,1] [,2] [,3] [,4] [,5]
[1,]    5    7    0    0    0
[2,]    6    8    0    0    0
[3,]    0    0    0    0    0
[4,]    0    0    0    0    0
, , 3
     [,1] [,2] [,3] [,4] [,5]
[1,]    0    0    0    0    0
[2,]    0    0    0    0    0
[3,]    0    0   -1   -3   -5
[4,]    0    0   -2   -4   -6
, , 4
     [,1] [,2] [,3] [,4] [,5]
[1,]    0    0    0    0    0
[2,]    0    0    0    0    0
[3,]    0    0   -7   -9   -2
[4,]    0    0   -8   -1   -3
, , 5
     [,1] [,2] [,3] [,4] [,5]
[1,]    0    0    0    0    0
[2,]    0    0    0    0    0
[3,]    0    0   -4   -6   -8
[4,]    0    0   -5   -7   -9
output-end

a.block.diag is not always the best solution. In case of
x<-array(1:8,rep(2,3)); y<-array(-1,rep(2,3))
the function adiag will be a little bit faster.

Peter Wolf



Robin Hankin wrote:

> Hi
>
> I have two arbitrarily dimensioned arrays, "a" and "b", with
> length(dim(a))==length(dim(b)).  I want to form a sort of
> "corner-to-corner" version of abind(), or a multidimensional version
> of blockdiag().
>
> In the case of matrices, the function is easy to write and if
> a=matrix(1,3,4) and b=matrix(2,2,2), then adiag(a,b) would return:
>
>      [,1] [,2] [,3] [,4] [,5] [,6]
> [1,]    1    1    1    1    0    0
> [2,]    1    1    1    1    0    0
> [3,]    1    1    1    1    0    0
> [4,]    0    0    0    0    2    2
> [5,]    0    0    0    0    2    2
>
>
> I am trying to generalize this to two higher dimensional arrays.
> If x <- adiag(a,b) then I want all(dim(x)==dim(a)+dim(b)); and if
> dim(a)=c(a_1, a_2,...a_d) then x[1:a_1,1:a_2,...,1:a_d]=a, and
> x[(a_1+1):(a_1+b_1),...,(a_d+1):(a_d+b_d)]=b.  Other elements of x are
> zero.
>
> The fact that I'm having difficulty expressing this succinctly makes
> me think I'm missing something basic.
>
> If a and b have identical dimensions [ie all(dim(a)==dim(b)) ], the
> following ghastly kludge (which is one of many) works:
>
> adiag <- function(a,b) {
>   if(any(dim(a) != dim(b))){stop("a and b must have identical 
> dimensions")}
>   jj <- array(0,rep(2,length(dim(a))))
>   jj[1] <- 1
>   jj[length(jj)] <- 1
>   jj <- kronecker(jj,b)
>   f <- function(i){1:i}
>   do.call("[<-",c(list(jj),sapply(dim(a),f,simplify=FALSE),list(a)))
> }
>
> Then "adiag(array(1:8,rep(2,3)),array(-1,rep(2,3)))" is OK.  What is
> the best way to bind arbitrarily dimensioned arrays together
> corner-to-corner?
>
>
>