8

I have two matrix A and B, so what's the fastest way to just calculate diag(A%*%B), i.e., the inner-product of the ith row of A and ith column of B, and the inner-product of other terms are not concerned.

supplement: A and B have large row and column numbers respectively.

Matthew Lundberg
  • 42,009
  • 6
  • 90
  • 112
Bs He
  • 717
  • 1
  • 10
  • 22
  • Possible duplicate of http://stackoverflow.com/questions/21708489/compute-only-diagonals-of-matrix-multiplication-in-r – akrun Mar 03 '17 at 03:04
  • It's similar but clearer as it involves only one matrix multiplication. – Andrey Shabalin Mar 03 '17 at 03:08
  • @akrun They are very similar, yes. But the linked question is a pattern used in mathematics, and this does not fit that pattern. That said, the answer to the linked question is in fact a near duplicate to the answer to this question. The other question may admit more answers. – Matthew Lundberg Mar 03 '17 at 03:14

1 Answers1

19

This can be done without full matrix multiplication, using just multiplication of matrix elements.

We need to multiply rows of A by the matching columns of B and sum the elements. Rows of A are columns of t(A), which we multiply element-wise by B and sum the columns.

In other words: colSums(t(A) * B)

Testing the code we first create sample data:

n = 5
m = 10000;

A = matrix(runif(n*m), n, m);
B = matrix(runif(n*m), m, n);

Your code:

diag(A %*% B)
# [1] 2492.198 2474.869 2459.881 2509.018 2477.591

Direct calculation without matrix multiplication:

colSums(t(A) * B)
# [1] 2492.198 2474.869 2459.881 2509.018 2477.591

The results are the same.

Andrey Shabalin
  • 4,389
  • 1
  • 19
  • 18
  • @AndreyShabalin Really a smart way. I would try this first. I think probably this is the fastest way with respect to software computing, but I think there should be better way with respect to algebra. – Bs He Mar 03 '17 at 03:13
  • It is numerically equivalent to `diag(A %*% B)`, the issue is there is no good name in linear algebra for `colSums` to phrase it my way. – Andrey Shabalin Mar 03 '17 at 03:16
  • It is not numerically equivalent to `diag(A %*% B)` because `A %*% B` performs many more multiplications. The linear algebra `colSums` would be `rep(1, m) %*% (t(A) * B)`. – IceCreamToucan Mar 20 '18 at 20:20
  • It performs many more multiplications, but always produces the same result. – Andrey Shabalin Mar 20 '18 at 21:41