A^K for Matrix Multiplication in R

A^k for matrix multiplication in R?

Although Reduce is more elegant, a for-loop solution is faster and seems to be as fast as expm::%^%

m1 <- matrix(1:9, 3)
m2 <- matrix(1:9, 3)
m3 <- matrix(1:9, 3)
system.time(replicate(1000, Reduce("%*%" , list(m1,m1,m1) ) ) )
#   user  system elapsed 
#  0.026   0.000   0.037 
mlist <- list(m1,m2,m3)
m0 <- diag(1, nrow=3,ncol=3)
system.time(replicate(1000, for (i in 1:3 ) {m0 <- m0 %*% m1 } ) )
#   user  system elapsed 
#  0.013   0.000   0.014 

library(expm)  # and I think this may be imported with pkg:Matrix
system.time(replicate(1000, m0%^%3))
# user  system elapsed 
#0.011   0.000   0.017 

On the other hand the matrix.power solution is much, much slower:

system.time(replicate(1000, matrix.power(m1, 4)) )
   user  system elapsed 
  0.677   0.013   1.037 

@BenBolker is correct (yet again). The for-loop appears linear in time as the exponent rises whereas the expm::%^% function appears to be even better than log(exponent).

> m0 <- diag(1, nrow=3,ncol=3)
> system.time(replicate(1000, for (i in 1:400 ) {m0 <- m0 %*% m1 } ) )
   user  system elapsed 
  0.678   0.037   0.708 
> system.time(replicate(1000, m0%^%400))
   user  system elapsed 
  0.006   0.000   0.006 

Matrix multiplication by loop in R

You can do:

cbind(A[, 1]*B, A[,2]*B) # or 
matrix(apply(A, 2, function(x) x*B), 2)

data

A <- matrix(1:4, 2)
B <- matrix(11:14, 2)

Matrix Multiplication in r

You could just use apply() to multiply each column of mat1 by mat2. (The "*" will carry out R's typically vectorized element-wise multiplication of the two equal-length vectors).

apply(mat1, 2, "*", mat2)
           [,1]       [,2]
[1,]  0.1785476  0.4175557
[2,]  0.2644247 -0.3745997
[3,] -0.5328542  0.8945527
[4,] -2.7351502 -0.7715341
[5,] -0.9719129 -0.1346929

Or better yet, convert mat1 to a vector to take advantage of R's recycling rules:

mat2 <- matrix(1:10, ncol=2)
mat1 <- matrix(1:5, ncol=1)

as.vector(mat1)*mat2
     [,1] [,2]
[1,]    1    6
[2,]    4   14
[3,]    9   24
[4,]   16   36
[5,]   25   50

Manually multiplying two matrices in R

Okay. I found the solution.

matrixMul <- function(mat1)
{
  rows <- nrow(mat1)
  cols <- ncol(mat1)

  matOut <- matrix(nrow = rows, ncol = cols) # empty matrix

  for (i in 1:rows) 
  {
    for(j in 1:cols)
    {
      vec1 <- mat1[i,]
      vec2 <- mat1[,j]

      mult1 <- vec1 * vec2

      matOut[i,j] <- sum(mult1)
    }
  }

  return(matOut) 
}

mat1 <- matrix(c(1,2,3,4), nrow = 2, ncol=2)
mult1 <- matrixMul(mat1)
mult1

Dimensions of matrix multiplication

You have 2 problems. One is that R eagerly converts single-column matrices to vectors, so K[1:2,1] is interpreted as a vector, not a matrix. We can fix this with the drop = FALSE argument to [.

There's also an order-of-operations problem, %*% has higher precedence than *, so we need parentheses to make the * happen first:

(K[1:2, 1, drop = F] * FF[1]) %*% K1[1, 1:2]
#       [,1]  [,2]
# [1,] 4e+07 2e+07
# [2,] 2e+07 1e+07

Unrelated to your problem, but it seems strange to make FF an array, seems like a length-100 vector would be simpler.

Matrix Multiplication in R without using %*% or crossprod

The only way I could think of solving it, is using embedded for loops (my first ever since 2013...)

The function

MatMult <- function(x, y){
  res <- matrix(NA, dim(x)[1], dim(y)[2])
    for(i in seq_along(y[1, ])){
     for(j in seq_along(x[, 1])){
      res[j, i] <- sum(x[j, ] * y[, i])
     }
    }
  res
 }

Your matrices

x <- matrix(1:4, ncol = 2)
y <- matrix(5:8, ncol = 2)

Testing

MatMult(x, y)
##      [,1] [,2]
## [1,]   23   31
## [2,]   34   46

x%*%y
##      [,1] [,2]
## [1,]   23   31
## [2,]   34   46

fast large matrix multiplication in R

There are many ways to approach this depending upon your code, effort, and hardware.

1. Use the 'best' function for the job

The simplest is to use crossprod which is the same as t(a)%*% b (Note - this will only be a small increase in speed)

crossprod(a,b) 

2. Use Rcpp (and likely RcppEigen/RcppArmadillo).

C++ will likely greater increase the speed of your code. Using linear algebra libraries will likely also help this further (hence Eigen and Armadillo). This however assumes you are willing to write some C++.

3. Use a better backend

After this point you are looking at BLAS backends such as OpenBLAS, Atlas, etc. Hooking these up to R varies depending upon your OS. It is quite easy if you are using a Debian system like Ubuntu. You can find a demo here. These can sometimes be leveraged further by libraries such as Armadillo and Eigen.

4. GPU Computing

If you have a GPU (e.g. AMD, NVIDIA, etc.) you can leverage the many cores within to greatly speed up your computations. There are a few that could be useful including gpuR, gputools, and gmatrix

EDIT - to address @jenesaiquoi comment on benefit of Rcpp

test.cpp

// [[Rcpp::depends(RcppArmadillo, RcppEigen)]]

#include <RcppArmadillo.h>
#include <RcppEigen.h>

// [[Rcpp::export]]
SEXP armaMatMult(arma::mat A, arma::mat B){
    arma::mat C = A * B;

    return Rcpp::wrap(C);
}

// [[Rcpp::export]]
SEXP eigenMatMult(Eigen::MatrixXd A, Eigen::MatrixXd B){
    Eigen::MatrixXd C = A * B;

    return Rcpp::wrap(C);
}

// [[Rcpp::export]]
SEXP eigenMapMatMult(const Eigen::Map<Eigen::MatrixXd> A, Eigen::Map<Eigen::MatrixXd> B){
    Eigen::MatrixXd C = A * B;

    return Rcpp::wrap(C);
}

test.R

library(Rcpp)

A <- matrix(rnorm(10000), 100, 100)
B <- matrix(rnorm(10000), 100, 100)

library(microbenchmark)
sourceCpp("test.cpp")
microbenchmark(A%*%B, armaMatMult(A, B), eigenMatMult(A, B), eigenMapMatMult(A, B))

Unit: microseconds
                  expr     min       lq     mean   median       uq      max neval
               A %*% B 885.846 892.1035 933.7457 901.1010 938.9255 1411.647   100
     armaMatMult(A, B) 846.688 857.6320 915.0717 866.2265 893.7790 1421.557   100
    eigenMatMult(A, B) 205.978 208.1295 233.1882 217.0310 229.4730  369.369   100
 eigenMapMatMult(A, B) 192.366 194.9835 207.1035 197.5405 205.2550  366.945   100

multiplying matrices with reference matrix n times

If you only want the last matrix we can use a recursive function, as suggested by Sotos:

A <- matrix(runif(9), 3) 
X <- matrix(runif(9), 3) 

repmult <- function(A,x,reps)
{
  if(reps==0){
    return(A)
  }
  else
  {
    repmult(A%*%x,x,reps-1)
  }
}

repmult(A,X,15)

---

If you want all intermediate results in a list as well, we can modify the function from the answer on this SO question (although you may want to change its name):

Mpow <- function(A,x, n) {
  L <- list(A%*%x)
  if (n==1) return(L)
  P <- A
  for (i in 1:n) L[[i]] <- (P <- P %*% x)
  return(L)
}

Mpow(A,X,3)

Matrix multiplication with integers or doubles in data table

In the first case, the column values of 'a' act as column index to subset the 'b', while in the second case it wont because the values are 0.1 to 0.5. So, instead of looping through the columns, loop through the sequence of the columns, and then subset the 'b' and the column values, which would work in both cases

a[, lapply(seq_along(.SD),function(x)(.SD[[x]]*b[,x]))]

---

Or just replicate the rows of 'a' and then multiply with 'b'

a[rep(seq_len(.N), nrow(b))] * b

---

Related Topics