Skip to content

Data Manipulation - make diagonal matrix of each element of a matrix

4 messages · Clemontina Alexander, Rui Barradas

#
Dear R list,
I have the following data:

set.seed(1)
n  <- 5     # number of subjects
tt <- 3     # number of repeated observation per subject
numco <- 2  # number of covariates
x <- matrix(round(rnorm(n*numco),2), ncol=numco)   # the actual covariates
x
[,1]  [,2]
[1,] -0.63 -0.82
[2,]  0.18  0.49
[3,] -0.84  0.74
[4,]  1.60  0.58
[5,]  0.33 -0.31

I need to form a matrix X such that
X =      x11      0      0     x21      0      0
              0   x11      0        0   x21      0
              0      0   x11        0      0   x21
           x12      0      0     x22      0      0
              0   x12      0        0   x22      0
              0      0   x12        0      0   x22
                       ...
           x15      0      0     x25      0      0
              0   x15      0        0   x25      0
              0      0   x15        0      0   x25
where both tt and numco can change. (So if tt=5 and numco=4, then X
needs to have 20 columns and n*tt rows. Each diagonal matrix should be
5x5 and there will be 4 of them for the 4 covariates.) I wrote this
funky for loop:

idd <- length(diag(1,tt))    # length of intercept matrix
X <- matrix(numeric(n*numco*idd),ncol=tt*numco)
for(i in 1:numco){
      X[,((i-1)*tt+1):(i*tt)] <- matrix(
        c(matrix(rep(diag(1,tt),n),ncol=tt, byrow=TRUE))   *
rep(rep(x[,i],each=tt),tt)
       , ncol=tt)
}
X

It works fine, but is there an easier way when n, tt, and numco get
larger and larger?
Thanks,
Tina


--
Clemontina Alexander
Ph.D Student
Department of Statistics
NC State University
Email: ckalexa2 at ncsu.com
#
I'm sorry, the indices of my X matrix are wrong.
It should be:

 X =     x11      0      0     x12      0      0
              0   x11      0        0   x12      0
              0      0   x11        0      0   x12
           x21      0      0     x22      0      0
              0   x21      0        0   x22      0
              0      0   x21        0      0   x22
                       ...
           xn1      0      0     x52      0      0
              0   xn1      0        0   x52      0
              0      0   xn1        0      0   x52

or

 X =     -0.63        0          0    -0.82         0         0
              0   -0.630          0         0    -0.82         0
              0          0   -0.630         0          0   -0.82
           0.18        0          0     0.49          0         0
              0     0.18          0         0      0.49         0
              0         0      0.18         0          0     0.49
                       ...
            0.33      0          0    -0.31          0          0
              0    0.33          0          0    -0.31          0
              0        0      0.33          0         0     -0.31

Sorry for the confusion.
Tina








On Thu, Dec 15, 2011 at 10:02 AM, Clemontina Alexander
<ckalexa2 at ncsu.edu> wrote:

  
    
#
Hello,

I believe I can help, or at least, my code is simpler.
First, look at your first line:

idd <- length(diag(1,tt))   # length of intercept matrix
#
not needed: diag(tt) would do the job but it's not needed,
    why call 2 functions, and one of them, 'diag', uses memory(*), if the
    result is tt squared? It's much simpler!
    (*)like you say, "larger and larger" amounts of it

My solution to your problem is as follows (as a function, and yours).

fun2 <- function(n, tt, numco){
    M.Unit <- matrix(rep(diag(1,tt),n), ncol=tt, byrow=TRUE)
    M <- NULL
    for(i in 1:numco) M <- cbind(M, M.Unit*rep(x[,i], each=tt))
    M
}

fun1 <- function(n, tt, numco){
    idd <- length(diag(1,tt))    # length of intercept matrix
    X <- matrix(numeric(n*numco*idd),ncol=tt*numco)
    for(i in 1:numco){
          X[,((i-1)*tt+1):(i*tt)] <- matrix(
            c(matrix(rep(diag(1,tt),n),ncol=tt, byrow=TRUE))*
    		rep(rep(x[,i],each=tt),tt)
           , ncol=tt)
    }
    X
}

I' ve tested the two with larger values of 'n', 'tt' and 'numco'
using the following timing instructions


n  <- 1000
tt <- 50
numco <- 15
set.seed(1)
x <- matrix(round(rnorm(n*numco),2), ncol=numco)   # the actual covariates

Runs <- 10^1

t1 <- system.time(for(i in 1:Runs) a1 <- fun1(n, tt, numco))[c(1,3)]
t2 <- system.time(for(i in 1:Runs) a2 <- fun2(n, tt, numco))[c(1,3)]

rbind(t1, t2, t1/t2)

      user.self     elapsed
t1 23.210000   31.060000
t2 14.970000   22.540000
     1.550434    1.377995

As you can see, it's not a great speed improvement.
I hope it's at least usefull.

Rui Barradas


--
View this message in context: http://r.789695.n4.nabble.com/Data-Manipulation-make-diagonal-matrix-of-each-element-of-a-matrix-tp4200321p4201305.html
Sent from the R help mailing list archive at Nabble.com.
#
Thank you, that is much simpler!
On Thu, Dec 15, 2011 at 2:04 PM, Rui Barradas <ruipbarradas at sapo.pt> wrote: