library(tidyverse)
library(knitr)
library(broom)
options(digits = 3)
set.seed(1234)
theme_set(theme_minimal())Matrix algebra provides an elegant way of representing both the data the kind of operations on tables or arrays that frequently come up in data analysis, and when implemented numerically, matrix algebra also provides efficient an means of carrying those operations out. The following is a brief introduction to matrix algebra as implmented in R. For more on theory and mathematics behind matricies, see materials from the Computational Math/Statistics Camp.
Create the \(3\) row by \(2\) column matrix \(\mathbf{A}\). Note the use of the concatenation c() function to collect the individual matrix elements (the \(a_{ij}\)’s) together, and the default fill order (byrow = FALSE), which implies filling the matrix by columns:
A <- matrix(c(6, 9, 12, 13, 21, 5), nrow = 3, ncol = 2)
A## [,1] [,2]
## [1,] 6 13
## [2,] 9 21
## [3,] 12 5class(A)## [1] "matrix"The class() function indicates that A is indeed a matrix (as opposed to a data frame).
Create another matrix, \(\mathbf{B}\), with the same elements, only filled by row this time:
B <- matrix(c(6, 9, 12, 13, 21, 5), nrow = 3, ncol = 2, byrow = TRUE)
B## [,1] [,2]
## [1,] 6 9
## [2,] 12 13
## [3,] 21 5Individual matrix elements can be reference using the square-bracket selection rules.
A[1, 2] # row 1, col 2## [1] 13A[2, 1] # row 2, col 1## [1] 9A[2, ] # all elements in row 2## [1] 9 21A[, 2] # all elements in column 2## [1] 13 21 5For comparison, create a vector \(\mathbf{c}\):
c <- c(1,2,3,4,5,6,7,8,9)
c## [1] 1 2 3 4 5 6 7 8 9class(c)## [1] "numeric"At this point, \(\mathbf{c}\) is just a list of numbers. The as.matrix() function creates a 9 row by 1 column vector, which can be verified by the dim() function:
c <- as.matrix(c)
class(c)## [1] "matrix"dim(c)## [1] 9 1The vector \(\mathbf{c}\) has 9 rows and 1 column.
A vector can be reshaped into a matrix:
C <- matrix(c, nrow = 3, ncol = 3)
C## [,1] [,2] [,3]
## [1,] 1 4 7
## [2,] 2 5 8
## [3,] 3 6 9A vector can also be created from a single row or column of a matrix:
a1 <- as.matrix(A[, 1]) # vector from column 1
a1## [,1]
## [1,] 6
## [2,] 9
## [3,] 12dim(a1)## [1] 3 1\(\mathbf{a1}\) is a 3 row by 1 column column vector.
Transposition t() flips the rows and columns of a matrix:
A## [,1] [,2]
## [1,] 6 13
## [2,] 9 21
## [3,] 12 5t(A)## [,1] [,2] [,3]
## [1,] 6 9 12
## [2,] 13 21 5C## [,1] [,2] [,3]
## [1,] 1 4 7
## [2,] 2 5 8
## [3,] 3 6 9t(C)## [,1] [,2] [,3]
## [1,] 1 2 3
## [2,] 4 5 6
## [3,] 7 8 9Vectors and also be transposed, which simply turns a column vector, e.g. \(\mathbf{a1}\) into a row vector
a1t <- t(a1)
a1t## [,1] [,2] [,3]
## [1,] 6 9 12dim(a1t)## [1] 1 3Matrix algebra is basically analogous to scalar algebra (with the exception of division), and obeys most of the same rules that scalar algebra does.
Add two matrices, \(\mathbf{A}\) and \(\mathbf{B}\):
F <- A + B
F## [,1] [,2]
## [1,] 12 22
## [2,] 21 34
## [3,] 33 10Note that the individual elements of \(\mathbf{A}\) and \(\mathbf{B}\) are simply added together to produce the corresponding elements of \(\mathbf{F}\) (i.e. \(f_{ij} = a_{ij} + b_{ij}\)).
In order to be added together, the matrices have to be of the same shape (i.e. have the same number of rows and colums). The shape of a matrix can be verified using the dim() function:
dim(C)## [1] 3 3dim(A)## [1] 3 2Here, \(\mathbf{A}\) and \(\mathbf{C}\) are not the same shape, and the following code, if executed, would product an error message:
G <- A + C## Error in A + C: non-conformable arraysScalar multiplication involves muliplying each element of a matrix by a scalar value:
H <- 0.5 * A
H## [,1] [,2]
## [1,] 3.0 6.5
## [2,] 4.5 10.5
## [3,] 6.0 2.5Here, \(h_{ij} = 0.5 \times a_{ij}\). Element-by-element multiplication is also possible for identically shaped matrices, e.g., \(p_{ij} = a_{ij} \times b_{ij}\):
P <- A * B
P## [,1] [,2]
## [1,] 36 117
## [2,] 108 273
## [3,] 252 25Matrix multiplication results in a a set of sums and crossproducts, as opposed to element-by-element products. Matrix multiplication is symbolized by the %*% operator:
Q <- C %*% A
Q## [,1] [,2]
## [1,] 126 132
## [2,] 153 171
## [3,] 180 210dim(C)## [1] 3 3dim(A)## [1] 3 2Note that the matrices have to be conformable, as they are here (the number of columns of the first matrix must equal the number of rows of the second, and the product matrix \(\mathbf{Q}\) here has the number of rows of the first matrix and the number of columns of the second).
The matrices \(\mathbf{A}\) and \(\mathbf{B}\) are not conformable for multiplication; although they have the same shape, they are not square, and the following code would produce an error:
T <- A %*% B## Error in A %*% B: non-conformable argumentsdim(A)## [1] 3 2dim(B)## [1] 3 2There are a number of special matrices that come up in data analysis. Here \(\mathbf{D}\) is a diagonal matrix, with non-zero values along the principal diagonal, and zeros elsewhere:
D <- diag(c(6,2,1,3), nrow = 4, ncol = 4)
D## [,1] [,2] [,3] [,4]
## [1,] 6 0 0 0
## [2,] 0 2 0 0
## [3,] 0 0 1 0
## [4,] 0 0 0 3A special form of a diagonal matrix is the identity matrix \(\mathbf{I}\), which has ones along the principal diagonal, and zeros elsewhere:
I <- diag(1, nrow = 4, ncol = 4)
I## [,1] [,2] [,3] [,4]
## [1,] 1 0 0 0
## [2,] 0 1 0 0
## [3,] 0 0 1 0
## [4,] 0 0 0 1A special scalar that appears often in practice is the norm (or Euclidean norm) of a vector, which is simply the square root of the sum of squares of the elements of the vector:
anorm <- sqrt(t(a1) %*% a1)
anorm## [,1]
## [1,] 16.2It can be verified that sqrt(sum(a1^2)) = \(16.155\).
The matrix algebra equivalent of division is the multiplication of one matrix by the inverse of another. Invertible matricies must be square and have a non-zero determinant. Consider the followng matrix \(\mathbf{R}\):
R <- matrix(c(0, 9, 3, 2, 2, 1, 9, 4, 4), nrow = 3, ncol = 3)
R## [,1] [,2] [,3]
## [1,] 0 2 9
## [2,] 9 2 4
## [3,] 3 1 4The inverse of the matrix \(\mathbf{R}\), \(\mathbf{R^{-1}}\), is obtained using the solve() function:
Rinv <- solve(R)
Rinv## [,1] [,2] [,3]
## [1,] -0.190 -0.0476 0.476
## [2,] 1.143 1.2857 -3.857
## [3,] -0.143 -0.2857 0.857As in scalar division, where \(a \times \frac{1}{a} = 1\), postmuliplying a matrix by its inverse yields the Identity matrix, \(\mathbf{I}\):
D <- R %*% Rinv
D## [,1] [,2] [,3]
## [1,] 1.00e+00 0 0.00e+00
## [2,] 1.11e-16 1 -4.44e-16
## [3,] 0.00e+00 0 1.00e+00After rounding using the zapsmall() function \(\mathbf{D}\) indeed equals \(\mathbf{I}\):
zapsmall(D)## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] 0 0 1An important concept that comes up in multivariate analysis is the decomposition of a matrix, into another square matrix, \(\mathbf{E}\), and a diagonal matrix, \(\mathbf{V}\), that each have some desireable properties, and which make the following statement true: \(\mathbf{RE} = \mathbf{EV}\), where \(\mathbf{V}\) is a diagonal matrix with the elements \(v_i\) along the diagonal. The matrix \(\mathbf{E}\) contains the eigenvectors of \(\mathbf{R}\), while the \(v_i\)’s are the eigenvalues of \(\mathbf{R}\).
E <- eigen(R)
E## eigen() decomposition
## $values
## [1] 9.920 -4.401 0.481
##
## $vectors
## [,1] [,2] [,3]
## [1,] -0.506 -0.611 0.0666
## [2,] -0.771 0.781 -0.9733
## [3,] -0.387 0.125 0.2198devtools::session_info()## Session info -------------------------------------------------------------## setting value
## version R version 3.5.1 (2018-07-02)
## system x86_64, darwin15.6.0
## ui RStudio (1.1.456)
## language (EN)
## collate en_US.UTF-8
## tz America/Chicago
## date 2018-10-25## Packages -----------------------------------------------------------------## package * version date source
## assertthat 0.2.0 2017-04-11 CRAN (R 3.5.0)
## backports 1.1.2 2017-12-13 CRAN (R 3.5.0)
## base * 3.5.1 2018-07-05 local
## bindr 0.1.1 2018-03-13 CRAN (R 3.5.0)
## bindrcpp * 0.2.2 2018-03-29 CRAN (R 3.5.0)
## broom * 0.5.0 2018-07-17 CRAN (R 3.5.0)
## cellranger 1.1.0 2016-07-27 CRAN (R 3.5.0)
## cli 1.0.0 2017-11-05 CRAN (R 3.5.0)
## colorspace 1.3-2 2016-12-14 CRAN (R 3.5.0)
## compiler 3.5.1 2018-07-05 local
## crayon 1.3.4 2017-09-16 CRAN (R 3.5.0)
## datasets * 3.5.1 2018-07-05 local
## devtools 1.13.6 2018-06-27 CRAN (R 3.5.0)
## digest 0.6.15 2018-01-28 CRAN (R 3.5.0)
## dplyr * 0.7.6 2018-06-29 cran (@0.7.6)
## emo 0.0.0.9000 2017-10-03 Github (hadley/emo@9f2e0f2)
## evaluate 0.11 2018-07-17 CRAN (R 3.5.0)
## fansi 0.3.0 2018-08-13 CRAN (R 3.5.0)
## forcats * 0.3.0 2018-02-19 CRAN (R 3.5.0)
## ggplot2 * 3.0.0 2018-07-03 CRAN (R 3.5.0)
## ggthemes * 4.0.0 2018-07-19 CRAN (R 3.5.0)
## glue 1.3.0 2018-07-17 CRAN (R 3.5.0)
## graphics * 3.5.1 2018-07-05 local
## grDevices * 3.5.1 2018-07-05 local
## grid 3.5.1 2018-07-05 local
## gtable 0.2.0 2016-02-26 CRAN (R 3.5.0)
## haven 1.1.2 2018-06-27 CRAN (R 3.5.0)
## highr 0.7 2018-06-09 CRAN (R 3.5.0)
## hms 0.4.2 2018-03-10 CRAN (R 3.5.0)
## htmltools 0.3.6 2017-04-28 CRAN (R 3.5.0)
## httpuv 1.4.5 2018-07-19 CRAN (R 3.5.0)
## httr 1.3.1 2017-08-20 CRAN (R 3.5.0)
## jsonlite 1.5 2017-06-01 CRAN (R 3.5.0)
## knitr * 1.20 2018-02-20 CRAN (R 3.5.0)
## labeling 0.3 2014-08-23 CRAN (R 3.5.0)
## later 0.7.3 2018-06-08 CRAN (R 3.5.0)
## lattice 0.20-35 2017-03-25 CRAN (R 3.5.1)
## lazyeval 0.2.1 2017-10-29 CRAN (R 3.5.0)
## lubridate 1.7.4 2018-04-11 CRAN (R 3.5.0)
## magrittr 1.5 2014-11-22 CRAN (R 3.5.0)
## memoise 1.1.0 2017-04-21 CRAN (R 3.5.0)
## methods * 3.5.1 2018-07-05 local
## mime 0.5 2016-07-07 CRAN (R 3.5.0)
## miniUI 0.1.1.1 2018-05-18 CRAN (R 3.5.0)
## modelr 0.1.2 2018-05-11 CRAN (R 3.5.0)
## munsell 0.5.0 2018-06-12 CRAN (R 3.5.0)
## nlme 3.1-137 2018-04-07 CRAN (R 3.5.1)
## patchwork * 0.0.1 2018-09-06 Github (thomasp85/patchwork@7fb35b1)
## pillar 1.3.0 2018-07-14 CRAN (R 3.5.0)
## pkgconfig 2.0.2 2018-08-16 CRAN (R 3.5.1)
## plyr 1.8.4 2016-06-08 CRAN (R 3.5.0)
## promises 1.0.1 2018-04-13 CRAN (R 3.5.0)
## purrr * 0.2.5 2018-05-29 CRAN (R 3.5.0)
## R6 2.2.2 2017-06-17 CRAN (R 3.5.0)
## rcfss * 0.1.5 2018-05-30 local
## Rcpp 0.12.18 2018-07-23 CRAN (R 3.5.0)
## readr * 1.1.1 2017-05-16 CRAN (R 3.5.0)
## readxl 1.1.0 2018-04-20 CRAN (R 3.5.0)
## rlang 0.2.1 2018-05-30 CRAN (R 3.5.0)
## rmarkdown 1.10 2018-06-11 CRAN (R 3.5.0)
## rprojroot 1.3-2 2018-01-03 CRAN (R 3.5.0)
## rsconnect 0.8.8 2018-03-09 CRAN (R 3.5.0)
## rstudioapi 0.7 2017-09-07 CRAN (R 3.5.0)
## rvest 0.3.2 2016-06-17 CRAN (R 3.5.0)
## scales 1.0.0 2018-08-09 CRAN (R 3.5.0)
## shiny 1.1.0 2018-05-17 CRAN (R 3.5.0)
## stats * 3.5.1 2018-07-05 local
## stringi 1.2.4 2018-07-20 CRAN (R 3.5.0)
## stringr * 1.3.1 2018-05-10 CRAN (R 3.5.0)
## tibble * 1.4.2 2018-01-22 CRAN (R 3.5.0)
## tidyr * 0.8.1 2018-05-18 CRAN (R 3.5.0)
## tidyselect 0.2.4 2018-02-26 CRAN (R 3.5.0)
## tidyverse * 1.2.1 2017-11-14 CRAN (R 3.5.0)
## tools 3.5.1 2018-07-05 local
## utf8 1.1.4 2018-05-24 CRAN (R 3.5.0)
## utils * 3.5.1 2018-07-05 local
## withr 2.1.2 2018-03-15 CRAN (R 3.5.0)
## xml2 1.2.0 2018-01-24 CRAN (R 3.5.0)
## xtable 1.8-2 2016-02-05 CRAN (R 3.5.0)
## yaml 2.2.0 2018-07-25 CRAN (R 3.5.0)This work is licensed under the CC BY-NC 4.0 Creative Commons License.