Apply function on two matrices

Question

For matrices Beta and x1, I am trying to apply re.fn to find the maximum for the (colSums(Beta*b))^2. Here, b is the column of the matrix x1. I wonder how I can make this algorithm faster.

set.seed(1)
D=10000
M=1000; N=1000

Beta=matrix(rnorm(N*D),ncol=D)
x1=matrix(rnorm(N*M),N)

re.fn <- function(b) {
  sum1 <- colSums(Beta*b)
  T_nc1 <- sum1^2
  T_nc <- max(T_nc1)
  return(T_nc)
}
T_nc=apply(x1,2,re.fn)

Answer 1

Using crossprod should be substantially faster

T_nc2 <- apply(crossprod(Beta,x1)^2,2,max)

all.equal(T_nc,T_nc2)
# [1] TRUE