Compute Only Diagonals of Matrix Multiplication in R

compute only diagonals of matrix multiplication in R

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

How to just calculate the diagonal of a matrix product in R

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.

compute only diagonal of product of matrices in Matlab

First of all, welcome! To be honest, this seems to be difficult. A little change is at least slightly increasing the speed:

N = 5000;
A = rand(N,N*2);
B = rand(N,N);

t = cputime;
diag(A'*B*A);
disp(['Elapsed cputime ' num2str(cputime-t)]);

t=cputime;
C = B*A;
sum(A'.*C',2);
disp(['Elapsed cputime ' num2str(cputime-t)]);

% slightly better...
t=cputime;
C = B*A;
sum(A.*C)';
disp(['Elapsed cputime ' num2str(cputime-t)]);

% slightly better than slightly better...
t=cputime;
sum(A.*(B*A))';
disp(['Elapsed cputime ' num2str(cputime-t)]);

Results:

Elapsed cputime 82.2593
Elapsed cputime 28.6106
Elapsed cputime 25.8338
Elapsed cputime 25.7714

Calculate the diagonal of a big.matrix in R

Ok, in fact, you don't need Rcpp here.
Just use the special matrix accessor of two-columns:

library(bigmemory)

X <- big.matrix(10, 10); X[] <- 1:100

d <- min(dim(X))
X[cbind(1:d, 1:d)]

X[cbind(1:d, 1:d)] will access X[1, 1], X[2, 2], ..., X[d, d].

How to compute only the diagonal of a matrix product in Octave?

The first element in the diagonal is the scalar product of the first row of A with the first column of B. The second element in the diagonal is the scalar product of the second row of A with the second column of B.

In other words:

vector = sum(A.*B',2);

Repeating a Matrix Multiplication with Different Blocks along its Diagonal

We can use

library(Matrix)
as.vector(as.matrix(bdiag(a[1:3, 1:3], a[4:6, 4:6], a[7:9, 7:9])) %*% b[1:9])
#[1]   38  200  362  866 1109 1352 2234 2558 2882

---

Or making it more dynamic

un1 <- unique(sub("\\d+", "", colnames(a)))
lst1 <- lapply(un1, function(x) a[grep(x, row.names(a)), 
          grep(x, colnames(a))])
as.matrix(bdiag(lst1)) %*% b[1:9]

---

Or if we want a matrix as output

out <- as.matrix(bdiag(Map(`%*%`, lst1, split(b[1:9], as.integer(gl(9, 3, 9))))))
row.names(out) <- row.names(a)
out
#     [,1] [,2] [,3]
#aaa1   38    0    0
#aaa2  200    0    0
#aaa3  362    0    0
#bbb1    0  866    0
#bb2    0 1109    0
#bbb3    0 1352    0
#ccc1    0    0 2234
#ccc2    0    0 2558
#ccc3    0    0 2882

NOTE: Using only base R installed packages

Loop diagonal multiplication - 7 * 7 matrix ... and so on

sapply(lapply(7:NCOL(df), function(i)
    df[, (i-6):i]), function(a)
        round(x = rev(cumprod(rev(diag(as.matrix(a))))), digits = 2))
#     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
#[1,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00
#[2,] 0.09 0.08 0.06 0.08 0.08 0.03 0.00
#[3,] 0.19 0.18 0.14 0.21 0.26 0.05 0.15
#[4,] 0.37 0.31 0.22 0.41 0.33 0.23 0.24
#[5,] 0.59 0.48 0.38 0.51 0.40 0.23 0.38
#[6,] 0.81 0.72 0.44 0.57 0.73 0.30 0.58
#[7,] 0.98 0.84 0.44 0.93 0.88 0.78 0.78

Let me know if the output is correct

DATA

df = structure(list(A = c(0.04, 0.44, 0.57, 0.61, 0.7, 0.76, 0.8), 
    B = c(0.03, 0.44, 0.57, 0.58, 0.65, 0.76, 0.84), C = c(0.03, 
    0.42, 0.52, 0.6, 0.66, 0.79, 0.89), D = c(0.04, 0.43, 0.59, 
    0.63, 0.71, 0.74, 0.84), E = c(0.04, 0.4, 0.62, 0.65, 0.73, 
    0.83, 0.82), F = c(0.07, 0.32, 0.51, 0.59, 0.66, 0.83, 0.83
    ), G = c(0.86, 0.64, 0.79, 0.81, 0.86, 0.86, 0.98), H = c(0.28, 
    0.02, 0.23, 0.83, 0.9, 1, 0.84), I = c(0.05, 0.33, 0.64, 
    1, 0.55, 0.61, 0.44), J = c(0.05, 0.36, 0.66, 0.63, 0.76, 
    0.83, 0.93), K = c(0.05, 0.3, 0.5, 0.57, 0.65, 0.38, 0.88
    ), L = c(0.04, 0.27, 0.55, 0.63, 0.66, 0.74, 0.78), M = c(0.04, 
    0.37, 0.6, 0.74, 0.74, 0.75, 0.78)), .Names = c("A", "B", 
"C", "D", "E", "F", "G", "H", "I", "J", "K", "L", "M"), class = "data.frame", row.names = c(NA, 
-7L))

vectoriced norm/matrix multiplication

You can do:

sigma <- matrix(c(1,0.5,0,0.5,1,0,0,0,1),3,3)
x <- t(matrix(rep(1:3, 10),3,10))

mynorm <- function(x, sig) t(x) %*% sig %*% x
apply(x, 1, mynorm, sig=sigma)

Here is a variant with tcrossprod():

mynorm <- function(x, sig) tcrossprod(x, sig) %*% x
apply(x, 1, mynorm, sig=sigma)

And here is the benchmark (including variants of the solution from compute only diagonals of matrix multiplication in R , thanks to @Benjamin for the link):

mynorm1 <- function(x, sig) t(x) %*% sig %*% x
mynorm2 <- function(x, sig) tcrossprod(x, sig) %*% x

microbenchmark(n1=apply(x, 1, mynorm1, sig=sigma), 
               n2=apply(x, 1, mynorm2, sig=sigma), 
               n3 = colSums(t(x) * (sigma %*% t(x))),
               n4 = rowSums(x * t(sigma %*% t(x))),
               n5 = rowSums(x * (x %*% t(sigma) )),
               n6 = rowSums(x * tcrossprod(x, sigma)),
               Eugen1 = diag(x %*% sigma %*% t(x)),
               Eugen2 = diag(x %*% tcrossprod(sigma, x)),
               unit="relative")

r: how to make matrix multiplication faster (special case)

I want to have a matrix d3 (m x n) whose each row is identical and
equal to a given vector (d0) with dimension n.

This is trivial to do using the matrix function and vector recycling.

m=4
n=5
set.seed(42)
d0=runif(n)
matrix(d0, nrow = m, ncol = n, byrow = TRUE)
#         [,1]      [,2]      [,3]      [,4]      [,5]
#[1,] 0.914806 0.9370754 0.2861395 0.8304476 0.6417455
#[2,] 0.914806 0.9370754 0.2861395 0.8304476 0.6417455
#[3,] 0.914806 0.9370754 0.2861395 0.8304476 0.6417455
#[4,] 0.914806 0.9370754 0.2861395 0.8304476 0.6417455

This should be the fastest solution.

---

Related Topics