Re: 3 Dimensional Cross Variance

[<Nick.Ellis@csiro.au>]

I take it that X is (n x p) and Q is (n x n x p). Then you can do the matrix 
operations by reshaping Q to (n x n*p), premultiplying by t(X), reshping the 
result to (p^2 x n), post-multiplying by X and reshaping back to a (p x p x p) 
array.

n <- 5;p <- 3
X <- matrix(1:(n*p),nrow=n,ncol=p)
Q <- array(1:(n*n*p),dim=c(n,n,p))

# use loop
tXQX.loop <- array(NA,dim=c(p,p,p))
for( i in 1 : p )
        for( j in 1 : p )
                for( k in 1 : p )                
                        tXQX.loop[ i, j, k ] = 0.5 * t( X[ , i ] ) %*% Q[ , , 
k] %*% X[ , j ] 

# use matrix multiplication
dim(Q) <- c(n,n*p)
tXQ <- t(X) %*% Q
dim(tXQ) <- c(p,n,p)
tXQ <- aperm(tXQ,perm=c(1,3,2))
dim(tXQ) <- c(p*p,n)
tXQX <- 0.5 * tXQ %*% X
dim(tXQX) <- c(p,p,p)
tXQX <- aperm(tXQX,perm=c(1,3,2))

range(tXQX.loop-tXQX)
[1] 0 0

Nick Ellis
CSIRO Marine Research   mailto:Nick.Ellis@csiro.au
PO Box 120                      ph    +61 (07) 3826 7260
Cleveland QLD 4163      fax   +61 (07) 3826 7222
Australia                       http://www.marine.csiro.au
  
-----Original Message-----
From: s-news-owner@lists.biostat.wustl.edu 
[mailto:s-news-owner@lists.biostat.wustl.edu]On Behalf Of Paul Lasky
Sent: Friday, 23 April 2004 9:26 AM
To: s-news
Subject: [S] 3 Dimensional Cross Variance


This is a pretty basic question but  . . . is there a better way, without 
explicit looping, to compute the  3 dimensional cross variance and a diagonal ? 
 The crossvariance defined as:
      crossvar = array( rep( NA,  p^3 ), dim=c ( p, p, p ) )
      for( i in 1 : p )
        {
          for( j in 1 : p )
           {
             for( k in 1 : p )
              {
                crossvar[ i, j, k ] = 0.5 * t( x[ , i ] ) % * % Q[ , , k] % * % 
x[ , j ] 
                   }
               }
         }
where x are the  bi-variates, and Q is the covariance array.

  The diag(crossvar) is defined as:

         crossvar.diag = matrix( rep( NA, p^2), nrow = p)
         for ( i in 1: p )
           {
             for( j in 1: p )
              {
                crossvar.diag[ i, j ] = 0.5 * t ( x[ , i] ) % * % Q[ , , j ] % 
* % x[ , i ]  
              }
          }

Paul Lasky
P & B Consultants
Rancho Santa Fe, CA .