15

I need only the diagonal elements from a matrix multiplication:

enter image description here,

in R. As Z is huge I want to avoid the full out multiplication....

Z <- matrix(c(1,1,1,2,3,4), ncol = 2)
Z
#     [,1] [,2]
#[1,]    1    2
#[2,]    1    3
#[3,]    1    4

X <- matrix(c(10,-5,-5,20), ncol = 2)
X
#     [,1] [,2]
#[1,]   10   -5
#[2,]   -5   20

Z %*% D %*% t(Z)
#     [,1] [,2] [,3]
#[1,]   70  105  140
#[2,]  105  160  215
#[3,]  140  215  290

diag(Z %*% D %*% t(Z))
#[1]  70 160 290

X is always a small square matrix (2x2 , 3x3 or 4x4), where Z will have the number of columns equal to the dimension of X. Is there a function available to do this?

guyabel
  • 8,014
  • 6
  • 57
  • 86

1 Answers1

19

I don't think you can avoid the first matrix multiplication (i.e. ZX), but you can the second one, which is the expensive one:

rowSums((Z %*% X) * Z)
# [1]  70 160 290

The second multiplication is NOT a matrix multiply. This is much faster:

library(microbenchmark)
set.seed(1)
X <- matrix(c(10,-5,-5,20), ncol = 2)
Z <- matrix(sample(1:1000), ncol=2)    # made Z a little bigger    

microbenchmark(
  res.new <- rowSums((Z %*% X) * Z),   # this solution
  res.old <- diag(Z %*% X %*% t(Z))    # original method
)
# Unit: microseconds
#                               expr     min      lq       mean   median        uq       max neval
#  res.new <- rowSums((Z %*% X) * Z)  20.956  23.233   34.77693  29.6150   44.0025    67.852   100
#  res.old <- diag(Z %*% X %*% t(Z)) 571.214 699.247 1761.08885 760.4295 1188.4485 47488.543   100     

all.equal(res.new, res.old)
# [1] TRUE
BrodieG
  • 51,669
  • 9
  • 93
  • 146
  • @Marcinthebox, thanks. Until you upvoted I was wondering if anyone would notice! – BrodieG Feb 12 '14 at 21:59
  • I had to try my own trial and error explorations for a while before I realized what you long discovered. Props – Marc in the box Feb 12 '14 at 22:48
  • This is brilliant. I am using it for SVD, so as my middle matrix is diagonal I can even call `rowSums(t(t(u) * d) * v)`. – Daniel Sparing Aug 01 '14 at 08:33
  • Why is the second matrix multiplication "the expensive one?" – Matthew Lundberg Mar 03 '17 at 03:15
  • 1
    @MatthewLundberg because the fist multiplication is between the big matrix and the small matrix, but the second multiplication is between the result of the first multiplication, which will be a big matrix, and another big matrix. – BrodieG Mar 03 '17 at 13:42