mxmult.r

# mxmult.r         common pitfall in matrix multiplication
#                  September 2008, revised September 2009
#
# clean up workspace
rm( list=ls() )
# slr design matrix X
X <- cbind( rep(1,6), (1:6) )
X

##      [,1] [,2]
## [1,]    1    1
## [2,]    1    2
## [3,]    1    3
## [4,]    1    4
## [5,]    1    5
## [6,]    1    6

# weight vector w
w <- sqrt( (1:6) )          # just so they're different
w

## [1] 1.000 1.414 1.732 2.000 2.236 2.449

# weight matrix
W <- diag(w)
# now (X'inv(W) X) four ways:
#   first          correct, but slow, W takes lots of space
W <- diag(w)
xwx1 <- t(X) %*% solve(W) %*% X      
xwx1

##       [,1]  [,2]
## [1,]  3.64 10.83
## [2,] 10.83 42.90

#   second         correct, less slow, takes lots of space 
xwx2 <- t(X) %*% diag(1/w) %*% X    
xwx2

##       [,1]  [,2]
## [1,]  3.64 10.83
## [2,] 10.83 42.90

#   third          wrong, but faster -- recycles w wrong
xwx3 <- ( t(X) / w ) %*% X           # *****wrong*****
xwx3

##        [,1]  [,2]
## [1,]  4.049 13.07
## [2,] 10.710 44.89

#   fourth         correct, but looks wrong
xwx4 <- t(X) %*% (X/w)               # uses recycling correctly, fast
xwx4

##       [,1]  [,2]
## [1,]  3.64 10.83
## [2,] 10.83 42.90

#                  looks different but same execution
crossprod(X,X/w)                     # correct, fast

##       [,1]  [,2]
## [1,]  3.64 10.83
## [2,] 10.83 42.90

# done
rm( list=ls() )
# q()