Mouk
Mouk

Reputation: 93

How could I build a function that extracts the diagonal block matrices of a larger one in R

How could I build a function that extracts the diagonal blocks matrices of a larger one? The problem is as follows. The function takes a centred matrix as argument, computes the full error covariance matrix and extracts the blocks on the leading diagonal? I tried the following, but not working.

    err_cov <- function(x){
       m <- nrow(x)
       n <- ncol(x)
       #compute the full error covariance matrix as the inner product
       #of vec(x) and its transpose. Note that, omega is a mnxmn matrix
       vec <- as.vector(x)
       omega <- vec%*%t(vec)
       sigmas <- list()
       for(i in 0:n-1){
         #here the blocks have to be m nxn matrices along the
         #leading diagonal
         for (j in 1:m)
           sigmas[[j]] <- omega[(n*i+1):n*(i+1), (n*i+1):n*(i+1)]
       }
       return(sigmas)
    }

So, for instance for

A [,1] [,2] [,3] [,4] [1,] 1 4 7 10 [2,] 2 5 8 11 [3,] 3 6 9 12

    > B<-as.vector(A)
    > B
     [1]  1  2  3  4  5  6  7  8  9 10 11 12

    > C<-B%*%t(B)
    > C
          [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
     [1,]    1    2    3    4    5    6    7    8    9    10    11    12
     [2,]    2    4    6    8   10   12   14   16   18    20    22    24
     [3,]    3    6    9   12   15   18   21   24   27    30    33    36
     [4,]    4    8   12   16   20   24   28   32   36    40    44    48
     [5,]    5   10   15   20   25   30   35   40   45    50    55    60
     [6,]    6   12   18   24   30   36   42   48   54    60    66    72
     [7,]    7   14   21   28   35   42   49   56   63    70    77    84
     [8,]    8   16   24   32   40   48   56   64   72    80    88    96
     [9,]    9   18   27   36   45   54   63   72   81    90    99   108
    [10,]   10   20   30   40   50   60   70   80   90   100   110   120
    [11,]   11   22   33   44   55   66   77   88   99   110   121   132
    [12,]   12   24   36   48   60   72   84   96  108   120   132   144

The function should return:

    > C1
         [,1] [,2] [,3]
    [1,]    1    2    3
    [2,]    2    4    6
    [3,]    3    6    9
    > C2
         [,1] [,2] [,3]
    [1,]   16   20   24
    [2,]   20   25   30
    [3,]   24   30   36
    > C3
         [,1] [,2] [,3]
    [1,]   49   56   63
    [2,]   56   64   72
    [3,]   63   72   81
    > C4
         [,1] [,2] [,3]
    [1,]  100  110  120
    [2,]  110  121  132
    [3,]  120  132  144

Thanks for answering.

Upvotes: 1

Views: 152

Answers (1)

Neal Fultz
Neal Fultz

Reputation: 9687

I think a clearer solution is to reset the dimensions and then let R do the index calculations for you:

err_cov <- function(x){                                                          
  m <- nrow(x)                                                                   
  n <- ncol(x)                                                                   
  #compute the full error covariance matrix as the inner product                 
  #of vec(x) and its transpose                                                   
  vec <- as.vector(x)                                                            
  omega <- tcrossprod(vec)                                                       
  dim(omega) <- c(n,m,n,m)                                                       
  sigmas <- list()                                                               
  for (j in 1:m)                                                                 
    sigmas[[j]] <- omega[,j,,j]                                                  

  return(sigmas)                                                                 
}    

Here is an example:

> x
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6
> tcrossprod(vec)
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    2    3    4    5    6
[2,]    2    4    6    8   10   12
[3,]    3    6    9   12   15   18
[4,]    4    8   12   16   20   24
[5,]    5   10   15   20   25   30
[6,]    6   12   18   24   30   36
> err_cov(x)
[[1]]
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    2    4    6
[3,]    3    6    9

[[2]]
     [,1] [,2] [,3]
[1,]   16   20   24
[2,]   20   25   30
[3,]   24   30   36

Upvotes: 1

Related Questions