In jcorrett/stats230.Rpackage: STATS230 R Functions

knitr::opts_chunk$set(echo = TRUE)
library(bench)
library(ggplot2)
library(tidyr)
is_available <- require(stats230.Rpackage)
if (!is_available) {
  library(devtools)
  install_github("jcorrett/stats230.Rpackage")
}
library(stats230.Rpackage) 

Github

The public Github repository for this R package, jcorrett/stats230.Rpackage, can be found online here.

Knitting this Rmd file should automatically install the package.

Test Function

For our multiplication function, we want to verify that we compute the products correctly. For this to be the case, we know that $$A\cdot (B\cdot\boldsymbol{x}) = (A\cdot B)\cdot\boldsymbol{x}.$$ So, to check that our function is working correctly, we compute $$||A\cdot (B\cdot\boldsymbol{x}) - (A\cdot B)\cdot\boldsymbol{x}||_\infty = \varepsilon.$$ If our function is working correctly, we should find that $\varepsilon \approx 1e-16$, the machine's working precision. We compute this value as follows on dense matrices $A,B \in \mathbb{R}^{10\times 10}$ and vector $\boldsymbol{x}\in\mathbb{R}^{10}$, where each entry of these matrices and vector is a random number drawn from a standard normal distribution. This does not guarantee that our matrices are full rank, but does ensure that they are dense and that little to no trivial computations are performed by the algorithm.

n <- 10
Atest <- matrix(rnorm(n^2),nrow=n,ncol=n)
Btest <- matrix(rnorm(n^2),nrow=n,ncol=n)
xtest <- matrix(rnorm(n^2),nrow=n,ncol=1)
prod1 <- matrixvectormultiply(Atest,Btest,xtest,TRUE)
prod2 <- matrixvectormultiply(Atest,Btest,xtest,FALSE)
error <- norm(prod1-prod2,type="I")
error

We see above that our error is in fact on the order of the machine precision, so we conclude that our algorithm is working correctly and we can continue with microbenchmarking.

Microbenchmarking

To microbenchmark this code, we will record the median computation time for a call with variable sized inputs, where $A,B \in \mathbb{R}^{n\times n}$, where $n \in [1,30]$, giving us 30 mean computation times for each way of computing $A\cdot B\cdot\boldsymbol{x}$. We can summarize the function performance with just the median computation time since each function call has only two matrix-matrix or matrix-vector multiplications, so the structure is simple enough to use one time measurement to summarize the whole function call. Then, using the bench package, we compute the median computation time as a function of $n$ for each algorithm as follows.

nmax <- 30
A <- matrix(rnorm(nmax^2),nrow=nmax,ncol=nmax)
B <- matrix(rnorm(nmax^2),nrow=nmax,ncol=nmax)
x <- matrix(rnorm(nmax^2),nrow=nmax,ncol=1)
method1_median <- c()
method2_median <- c()
nvec <- 1:nmax
for (n in nvec) {
  lb <- mark(matrixvectormultiply(A[1:n,1:n],B[1:n,1:n],matrix(x[1:n]),TRUE),matrixvectormultiply(A[1:n,1:n],B[1:n,1:n],matrix(x[1:n]),FALSE),min_time = Inf,max_iterations = 10000)
  method1_median[n] <- lb$median[1]
  method2_median[n] <- lb$median[2]
}

df <- data.frame(nvec,method1_median,method2_median)

Using the above data, we then summarize our results in the following plot:

cols = c("A*(B*x)"="red","(A*B)*x"="blue")
bench_plot <- ggplot(data = df, aes(x=nvec)) + geom_line(aes(y = method1_median,color="A*(B*x)"), size=1) + geom_line(aes(y = method2_median,color="(A*B)*x"),size=1) + labs(x = "n",y = "Median Computation Time (s)",color = "Legend",title = "Multiplication Algorithm Computation Time Summary")+ scale_color_manual(values = cols)+theme(plot.title = element_text(hjust = 0.5))
bench_plot 

From the above, we see that both algorithms compute our product on the order of microseconds, which is relatively quick. However, we see that the first algorithm, computing the product as $(A\cdot B)\cdot\boldsymbol{x}$, has a computation time that increases with $n$ much more rapidly than our other method. This can be attributed to multiplying two matrices scales $\sim O(n^3)$, but matrix-vector multiplication scales $\sim O(n^2)$. This difference leads to the significant difference in computation times, as expected.

To analyze further, we summarize the function call for $n=30$ with the following plot.

lb <- mark(matrixvectormultiply(A,B,x,TRUE),matrixvectormultiply(A,B,x,FALSE),min_time = Inf,max_iterations = 10000)
bench_plot <- plot(lb) + labs(color = "Legend",title = "Multiplication Algorithm Computation Time for n = 30")+theme(plot.title = element_text(hjust = 0.5))+theme(legend.position = "none")
bench_plot 

From this plot, we can see that the median is a good measure for this function, since the vast majority of the density of each distribution of computation times is near a constant value, with only a small handful of rare events exceeding the expected computation time. However, these rare events do not significantly violate the conclusion that the matrix-matrix multiplication algorithm is much slower than the matrix-vector algorithm due to its $n^3$ complexity.

jcorrett/stats230.Rpackage documentation built on March 21, 2022, 5:36 a.m.