Linear Algebra / Matrix Operations in R

Informal notes about common linear algebra / matrix operations in R
Author

Max Rohde

Published

April 29, 2022

The following material should be understandable to anyone with a basic knowledge of R and linear algebra. I plan to update this page with more material in the future.

Creating vectors

Vectors can be created with the c() function.

c(1,3,7)
[1] 1 3 7

For vectors that are a continuous sequence of numbers, you can use the : operator.

1:8
[1] 1 2 3 4 5 6 7 8

For vectors of repeated numbers, such as a zero vector, use the rep() function.

# Create a zero vector of length 10
rep(0, times=10)
 [1] 0 0 0 0 0 0 0 0 0 0
# Or to be more concise
# rep(0, 10)

More complicated sequences can be created with the seq() function.

# Create a sequence from 1 to 3 with a step size of 0.5
seq(1,3, by = 0.5)
[1] 1.0 1.5 2.0 2.5 3.0
# Create 10 equal spaced numbers between 1 and 3
seq(1,3, length.out = 10)
 [1] 1.000000 1.222222 1.444444 1.666667 1.888889 2.111111 2.333333 2.555556
 [9] 2.777778 3.000000

Creating matrices

There are two common ways to create matrices in R.

The first method is using the matrix() function. You pass the elements of the matrix into matrix() and specify the number of rows and columns. Note that R fills in the numbers going down each column. This can be unintuitive.

# Put the numbers 1 to 12 in 3 rows and 4 columns
matrix(1:12,
       nrow=3,
       ncol=4)
     [,1] [,2] [,3] [,4]
[1,]    1    4    7   10
[2,]    2    5    8   11
[3,]    3    6    9   12
# Put the numbers 1 to 12 in 2 rows and 6 columns
matrix(1:12,
       nrow=2,
       ncol=6)
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    3    5    7    9   11
[2,]    2    4    6    8   10   12

The second method is using rbind() and cbind(). Use rbind() to combine vectors into a matrix by row, and cbind() to combine vectors into a matrix by column.

x1 <- c(1,2,3,4)
x2 <- c(5,6,7,8)
x3 <- c(9,10,11,12) 
rbind(x1, x2, x3)
   [,1] [,2] [,3] [,4]
x1    1    2    3    4
x2    5    6    7    8
x3    9   10   11   12
cbind(x1, x2, x3)
     x1 x2 x3
[1,]  1  5  9
[2,]  2  6 10
[3,]  3  7 11
[4,]  4  8 12

To create a \(n \times n\) identity matrix, use the diag(n) function.

diag(5)
     [,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

Once a matrix is created, the dim() can be used to obtain the dimensions.

m <- cbind(x1, x2, x3)

# See that m has 4 rows and 3 columns
dim(m)
[1] 4 3

If you are only interested in the number of rows and columns, use nrow() and ncol() respectively.

nrow(m)
[1] 4
ncol(m)
[1] 3

Matrix indexing

Once you have a matrix in R, how do you subset parts of the matrix? Let’s use this matrix, m as an example.

m <- matrix(1:12,
            nrow=3,
            ncol=4)

print(m)
     [,1] [,2] [,3] [,4]
[1,]    1    4    7   10
[2,]    2    5    8   11
[3,]    3    6    9   12

The syntax m[i] selects the ith element. Recall that R counts going down columns.

m[2]
[1] 2

The syntax m[i, ] selects the ith row

m[2,]
[1]  2  5  8 11

The syntax m[,i] selects the ith column.

m[,2]
[1] 4 5 6

The syntax m[i,j] selects the element in the ith row and jth column.

m[2,2]
[1] 5

Matrix operations

Matrix multiplication uses the %*% operator.

# Define matrices
A <- matrix(1:9, nrow=3, ncol=3)
B <- matrix(1:6, nrow=3, ncol=2)

# Matrix multiplication
A %*% B
     [,1] [,2]
[1,]   30   66
[2,]   36   81
[3,]   42   96

Remember that the order of matrix multiplication is important!

B %*% A
Error in B %*% A : non-conformable arguments

To find the inverse of a matrix, we use the (somewhat unintuitively named) solve() function.

A <- matrix(runif(9),
            nrow=3,
            ncol=3)
A
          [,1]      [,2]      [,3]
[1,] 0.3784592 0.3445372 0.6930376
[2,] 0.2448913 0.4762224 0.2208319
[3,] 0.2830958 0.6282155 0.5572572
          [,1]       [,2]      [,3]
[1,]  3.553595  6.8289748 -7.125670
[2,] -2.074969  0.4125499  2.417066
[3,]  0.533901 -3.9343129  2.689617

Multiplying A by its inverse, we see that we obtain the identity matrix. The numbers very close to zero in the off-diagonal elements are due to floating-point error.

solve(A) %*% A
             [,1]          [,2]         [,3]
[1,] 1.000000e+00  0.000000e+00 0.000000e+00
[2,] 2.220446e-16  1.000000e+00 2.220446e-16
[3,] 1.110223e-16 -2.220446e-16 1.000000e+00

To take the transpose of a matrix, use the t() function.

B
     [,1] [,2]
[1,]    1    4
[2,]    2    5
[3,]    3    6
t(B)
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    4    5    6