Linear Algebra in R
Contents
Linear Algebra in R#
Here I briefly introduce the use of matrix algebra manipulations for R. Other programs are arguably better for pure linear algebra work (e.g. Matlab or Julia), but R is a very good environment for mixing modeling including running pre-packaged statistical commands and linear algebra. In my opinion stata’s linear algebra code is a lot more intuitive than R’s.
Creating vectors and matrices#
Suppose we want to define a row vector \(\mathbf{a}\) as
we can enter this in R as
a <- cbind(1,2,3)
print(a)
[,1] [,2] [,3]
[1,] 1 2 3
If instead, we wanted a column vector
we use this
a <- rbind(1,2,3)
print(a)
[,1]
[1,] 1
[2,] 2
[3,] 3
Creating Matrices and Vectors from DataFrames#
Sometimes we want to use information in R DataFrames for matrix operations. To convert part (or all) of a DataFrame to a matrix, use data.matrix(frame). Here is an example:
library(foreign)
tk.df = read.dta("http://rlhick.people.wm.edu/econ407/data/tobias_koop.dta")
tk.matrix = data.matrix(tk.df)
or if we only need a few columns and only some rows:
tk.matrix = data.matrix(tk.df[tk.df$time == 4,c("educ", "pexp", "ability")])
to construct a matrix with time period 4 data and the education, experience, and ability columns.
Referencing elements in arrays#
We can grab individual elements of this vector using slicing:
print(a[1:2])
print(a[3])
[1] 1 2
[1] 3
Note, since we are working with a column, we don’t need to refer to the row dimension, although we could:
print(a[1:2,1])
print(a[3,1])
[1] 1 2
[1] 3
Creating this matrix $\( B = \begin{bmatrix} 2&4&3\\1&5&7 \end{bmatrix} \)$ is done this way:
B <- matrix(c(2, 4, 3, 1, 5, 7), nrow = 2, byrow = TRUE)
print(B)
[,1] [,2] [,3]
[1,] 2 4 3
[2,] 1 5 7
slicing is also possible
print(B[,2:3])
[,1] [,2]
[1,] 4 3
[2,] 5 7
It is also to create an empty matrix (all zeros) and fill it in using slicing:
C <- matrix(0, 5, 5)
C[3:4, 5] = -999
print(C)
[,1] [,2] [,3] [,4] [,5]
[1,] 0 0 0 0 0
[2,] 0 0 0 0 0
[3,] 0 0 0 0 -999
[4,] 0 0 0 0 -999
[5,] 0 0 0 0 0
The identity matrix#
These are created using diag, having the arguments:
Value to place on the diagonal
Number of rows
Number of columns (if ommitted, columns=rows)
I = diag(1,5)
print(I)
[,1] [,2] [,3] [,4] [,5]
[1,] 1 0 0 0 0
[2,] 0 1 0 0 0
[3,] 0 0 1 0 0
[4,] 0 0 0 1 0
[5,] 0 0 0 0 1
or, if you only supply one argument (row, column dimension):
I = diag(5)
print(I)
[,1] [,2] [,3] [,4] [,5]
[1,] 1 0 0 0 0
[2,] 0 1 0 0 0
[3,] 0 0 1 0 0
[4,] 0 0 0 1 0
[5,] 0 0 0 0 1
Creating a column vector of ones#
print(cbind(rep(1, 5)))
[,1]
[1,] 1
[2,] 1
[3,] 1
[4,] 1
[5,] 1
Getting information about your matrices and vectors#
The most important information we need is the dimensions of our matrices. The function dim tells us rows and columns (for 2 dimensional objects):
dim(B)
- 2
- 3
We can extract row and column dimensions like this:
dim(B)[1]
dim(B)[2]
Linear Algebra Operations#
Scalar Addition#
a = 1
b = 1
print(a + b)
[1] 2
Matrix Addition#
a = rbind(2, 5, 8)
b = rbind(6, 4, 3)
print(b + a)
[,1]
[1,] 8
[2,] 9
[3,] 11
Note, conformability matters:
a = rbind(2, 5)
b = rbind(6, 4, 3)
print(b + a)
Error in b + a: non-conformable arrays
Traceback:
1. print(b + a)
Matrix Multiplication#
Matrix multiplication uses the %*% operator:
dim(B)
dim(b)
print(B%*%b)
- 2
- 3
- 3
- 1
[,1]
[1,] 37
[2,] 47
Again, comformability (and order) matters:
print(b%*%B)
Error in b %*% B: non-conformable arguments
Traceback:
1. print(b %*% B)
Matrix Transpose#
The transpose operator is t():
print(B)
print(t(B))
[,1] [,2] [,3]
[1,] 2 4 3
[2,] 1 5 7
[,1] [,2]
[1,] 2 1
[2,] 4 5
[3,] 3 7
Matrix Inversion#
A <- matrix(c(2, 4, 3, 4, 9, 1, 1, 5, 7), nrow = 3, byrow = TRUE)
A_inv = solve(A)
print(A_inv)
[,1] [,2] [,3]
[1,] 1.4146341 -0.3170732 -0.56097561
[2,] -0.6585366 0.2682927 0.24390244
[3,] 0.2682927 -0.1463415 0.04878049
and the matrix A_inv satisfies the properties of an inverse:
print(A_inv%*%A)
[,1] [,2] [,3]
[1,] 1.000000e+00 4.440892e-16 -8.881784e-16
[2,] -1.942890e-16 1.000000e+00 4.163336e-16
[3,] 6.245005e-17 1.457168e-16 1.000000e+00
Scalar Operations#
Scalar Addition and Subtraction
print(a)
print(a - 5)
[,1]
[1,] 2
[2,] 5
[,1]
[1,] -3
[2,] 0
Scalar Multiplication and Addition
print(a/5)
print(a*5)
[,1]
[1,] 0.4
[2,] 1.0
[,1]
[1,] 10
[2,] 25
Elementwise Operations#
Sometimes we want to combine arrays by performing arithmetic on the corresponding elements. For example, supposing that
and
we might want to calculate
by default R performs these operations using basic arithmetic operators so long as the arrays are of the same dimensions. Note in the first example below, \(\mathbf{b}\) and \(\mathbf{a}\) are not conformable.
print(dim(a))
print(dim(b))
print(a/b)
[1] 2 1
[1] 3 1
Error in a/b: non-conformable arrays
Traceback:
1. print(a/b)
However, if we redefine \(\mathbf{b}\) for having conformable dimensions, element-wise division exists
b=c(6,4)
print(a/b)
[,1]
[1,] 0.3333333
[2,] 1.2500000
Some other useful commands in R#
Sometimes we want to extract diagonal elements of a matrix or set off diagonal elements of a matrix to zero. In R, we use the diag command for these cases:
To grab the diagonal elements as a vector:
print(diag(B))
[1] 2 5
and to set off diagonal elements of \(\mathbf{B}\) to zero, apply diag twice:
print(diag(diag(B)))
[,1] [,2]
[1,] 2 0
[2,] 0 5