Introduction to linear algebra in R

Søren Højsgaard,
Department of Mathematical Sciences,
Aalborg University, Denmark
https://people.math.aau.dk/~sorenh/

Contents

1  Introduction

2  Vectors

2.1  Vectors

• In what follows all vectors are column vectors unless otherwise stated.
• However, writing column vectors takes up more space than row vectors. Therefore we shall frequently write vectors as row vectors, but with the understanding that it really is a column vector.

Figure 1: Two 2-vectors

> a <- c(1, 3, 2)
> a
 

[1] 1 3 2
 

The vector a is in R printed "in row format" but can really be regarded as a column vector, cfr. the convention above.

2.2  Transpose of vectors

Transposing a vector means turning a column (row) vector into a row (column) vector. The transpose is denoted by "^T ".



⎡ ⎢ ⎢ ⎢ ⎣ 1 3 2 ⎤ ⎥ ⎥ ⎥ ⎦ T   = [1,3,2]    and   [1,3,2]T = ⎡ ⎢ ⎢ ⎢ ⎣ 1 3 2 ⎤ ⎥ ⎥ ⎥ ⎦

Hence transposing twice takes us back to where we started:


a = (aT )T

> t(a)
 

     [,1] [,2] [,3]
[1,]    1    3    2
 

2.3  Multiplying a vector by a number

If a is a vector and α is a number then αa is the vector


αa = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ αa1 αa2 : αan ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

See Figure 2.



7 ⎡ ⎢ ⎢ ⎢ ⎣ 1 3 2 ⎤ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎣ 7 21 14 ⎤ ⎥ ⎥ ⎥ ⎦


Figure 2: Multiplication of a vector by a number

> 7 * a
 

[1]  7 21 14
 

2.4  Sum of vectors

Let a and b be n-vectors. The sum a+b is the n-vector


a + b = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ a1 a2 : an ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ + ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ b1 b2 : bn ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ a1 + b1 a2 + b2 : an + bn ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = b+a

See Figure 3 and 4. Only vectors of the same dimension can be added.



⎡ ⎢ ⎢ ⎢ ⎣ 1 3 2 ⎤ ⎥ ⎥ ⎥ ⎦ + ⎡ ⎢ ⎢ ⎢ ⎣ 2 8 9 ⎤ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎣ 1+2 3+8 2+9 ⎤ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎣ 3 11 11 ⎤ ⎥ ⎥ ⎥ ⎦


Figure 3: Addition of vectors

Figure 4: Addition of vectors and multiplication by a number

> a <- c(1, 3, 2)
> b <- c(2, 8, 9)
> a+b
 

[1]  3 11 11
 

2.5  Inner product of vectors

Let a=(a₁,..., a_n) and b=(b₁,...,b_n). The inner product of a and b is


a·b = a1b1 + ... + anbn

Note, that the inner product is a number - not a vector.

> sum(a * b)
 

[1] 44
 

2.6  The length (norm) of a vector

The length (or norm) of a vector a is


|| a || = √   a ·a   =   ⎛ √ n ∑ i=1  ai2

> sqrt(sum(a * a))
 

[1] 3.741657
 

2.7  The 0-vector and 1-vector

The 0-vector (1-vector) is a vector with 0 (1) on all entries. The 0-vector (1-vector) is frequently written simply as 0 (1) or as 0_n (1_n) to emphasize that its length n.

> rep(0, 5)
 

[1] 0 0 0 0 0
 

> rep(1, 5)
 

[1] 1 1 1 1 1
 

2.8  Orthogonal (perpendicular) vectors

Two vectors v₁ and v₂ are orthogonal if and only if their inner product is zero, written


v1 ⊥v2 ⇔ v1·v2 = 0

Note that any vector is orthogonal to the 0-vector.

> v1 <- c(1, 1)
> v2 <- c(-1, 1)
> sum(v1 * v2)
 

[1] 0
 

3  Matrices

3.1  Matrices

> A <- matrix(c(1, 3, 2, 2, 8, 9), ncol=3)
> A
 

     [,1] [,2] [,3]
[1,]    1    2    8
[2,]    3    2    9
 

Note that the numbers 1,3,2,2,8,9 are read into the matrix column-by-column. To get the numbers read in row-by-row do

> A2 <- matrix(c(1, 3, 2, 2, 8, 9), ncol=3, byrow=T)
> A2
 

     [,1] [,2] [,3]
[1,]    1    3    2
[2,]    2    8    9
 

3.2  Multiplying a matrix with a number

For a number α and a matrix A, the product αA is the matrix obtained by multiplying each element in A by α.




7 ⎡ ⎢ ⎢ ⎢ ⎣ 1 2 3 8 2 9 ⎤ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎣ 7 14 21 56 14 63 ⎤ ⎥ ⎥ ⎥ ⎦

> 7 * A
 

     [,1] [,2] [,3]
[1,]    7   14   56
[2,]   21   14   63
 

3.3  Transpose of matrices

A matrix is transposed by interchanging rows and columns and is denoted by "^T ".




⎡ ⎢ ⎢ ⎢ ⎣ 1 2 3 8 2 9 ⎤ ⎥ ⎥ ⎥ ⎦ T   = ⎡ ⎢ ⎢ ⎣ 1 3 2 2 8 9 ⎤ ⎥ ⎥ ⎦

Note that if A is an r×c matrix then A^T is a c× r matrix.

> t(A)
 

     [,1] [,2]
[1,]    1    3
[2,]    2    2
[3,]    8    9
 

3.4  Sum of matrices

Let A and B be r×c matrices. The sum A+B is the r×c matrix obtained by adding A and B elementwise.
Only matrices with the same dimensions can be added.



⎡ ⎢ ⎢ ⎢ ⎣ 1 2 3 8 2 9 ⎤ ⎥ ⎥ ⎥ ⎦ + ⎡ ⎢ ⎢ ⎢ ⎣ 5 4 8 2 3 7 ⎤ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎣ 6 6 11 10 5 16 ⎤ ⎥ ⎥ ⎥ ⎦

> B <- matrix(c(5, 8, 3, 4, 2, 7), ncol=3, byrow=T)
> A + B
 

     [,1] [,2] [,3]
[1,]    6   10   11
[2,]    7    4   16
 

3.5  Multiplication of a matrix and a vector

Let A be an r×c matrix and let b be a c-dimensional column vector. The product Ab is the r×1 matrix


Ab = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ a11 a12 ... a1c a21 a22 ... a2c : : ··· : ar1 ar2 ... arc ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ b1 b2 : bc ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ a11b1 + a12b2 + ... + a1cbc a21b1 + a22b2 + ... + a2cbc : ar1b1 + ar2b2 + ... + arcbc ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦




⎡ ⎢ ⎢ ⎢ ⎣ 1 2 3 8 2 9 ⎤ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎣ 5 8 ⎤ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎣ 1 ·5 + 2 ·8 3 ·5 + 8 ·8 2 ·5 + 9 ·8 ⎤ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎣ 21 79 82 ⎤ ⎥ ⎥ ⎥ ⎦

> A %*% a
 

     [,1]
[1,]   23
[2,]   27
 

Note the difference to

> A * a
 

     [,1] [,2] [,3]
[1,]    1    4   24
[2,]    9    2   18
 

Please figure out yourself what goes on!

3.6  Multiplication of matrices

Let A be an r×c matrix and B a c×t matrix, i.e.  B=[b₁:b₂:...:b_t]. The product AB is the r×t matrix given by:


AB = A [b1:b2:...:bt] = [Ab1:Ab2:...:Abt]




⎡ ⎢ ⎢ ⎢ ⎣ 1 2 3 8 2 9 ⎤ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎣ 5 4 8 2 ⎤ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎣ ⎡ ⎢ ⎢ ⎢ ⎣ 1 2 3 8 2 9 ⎤ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎣ 5 8 ⎤ ⎥ ⎥ ⎦  :   ⎡ ⎢ ⎢ ⎢ ⎣ 1 2 3 8 2 9 ⎤ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎣ 4 2 ⎤ ⎥ ⎥ ⎦ ⎤ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎣ 1 ·5 + 2 ·8  1 ·4 + 2 ·2 3 ·5 + 8 ·8 3 ·4 + 8 ·2 2 ·5 + 9 ·8 2 ·4 + 9 ·2 ⎤ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎣ 21 8 79 28 82 26 ⎤ ⎥ ⎥ ⎥ ⎦

Note that the product AB can only be formed if the number of rows in B and the number of columns in A are the same. In that case, A and B are said to be conforme.
In general AB and BA are not identical.
A is :


⎡ ⎢ ⎢ ⎢ ⎣ 1 2 3 8 2 9 ⎤ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎣ 5 4 8 2 ⎤ ⎥ ⎥ ⎦ = 5 4 8 2 1 2 1 ·5 + 2 ·8 1 ·4 + 2 ·2 3 8 3 ·5 + 8 ·8 3 ·4 + 8 ·2 2 9 2 ·5 + 9 ·8 2 ·4 + 9 ·2 = ⎡ ⎢ ⎢ ⎢ ⎣ 21 8 79 28 82 26 ⎤ ⎥ ⎥ ⎥ ⎦

> A <- matrix(c(1, 3, 2, 2, 8, 9), ncol=2)
> B <- matrix(c(5, 8, 4, 2), ncol=2)
> A %*% B
 

     [,1] [,2]
[1,]   21    8
[2,]   79   28
[3,]   82   26
 

3.7  Vectors as matrices

One can regard a column vector of length r as an r×1 matrix and a row vector of length c as a 1× c matrix.

3.8  Some special matrices

• An n×n matrix is a
• A matrix A is if A=A^T .
• A matrix with 0 on all entries is the and is often written simply as 0.
• A matrix consisting of 1s in all entries is often written J.
• A square matrix with 0 on all off-diagonal entries and elements d₁, d₂, ..., d_n on the diagonal a and is often written diag{d₁, d₂, ..., d_n}
• A diagonal matrix with 1s on the diagonal is called the and is denoted I. The identity matrix satisfies that IA=AI=A. Likewise, if x is a vector then Ix=x.
• 0-matrix and 1-matrix

  > matrix(0, nrow=2, ncol=3)
   
  

       [,1] [,2] [,3]
  [1,]    0    0    0
  [2,]    0    0    0
   
  

  > matrix(1, nrow=2, ncol=3)
   
  

       [,1] [,2] [,3]
  [1,]    1    1    1
  [2,]    1    1    1
   
  
• Diagonal matrix and identity matrix

  > diag(c(1, 2, 3))
   
  

       [,1] [,2] [,3]
  [1,]    1    0    0
  [2,]    0    2    0
  [3,]    0    0    3
   
  

  > diag(1, 3)
   
  

       [,1] [,2] [,3]
  [1,]    1    0    0
  [2,]    0    1    0
  [3,]    0    0    1
   
  

  Note what happens when diag is applied to a matrix:

  > diag(diag(c(1, 2, 3)))
   
  

  [1] 1 2 3
   
  

  > diag(A)
   
  

  [1] 1 8
   
  

  3.9  Inverse of matrices

  In general, the inverse of an n×n matrix A is the matrix B (which is also n×n) which when multiplied with A gives the identity matrix I. That is,
  

  AB = BA = I.

  One says that B is A's inverse and writes B=A⁻¹. Likewise, A is Bs inverse.
  Let
  

  A= ⎡ ⎢ ⎢ ⎣ 1 3 2 4 ⎤ ⎥ ⎥ ⎦    B= ⎡ ⎢ ⎢ ⎣ −2 1.5 1 −0.5 ⎤ ⎥ ⎥ ⎦

  Now AB = BA = I so B=A⁻¹.
  If A is a 1 ×1 matrix, i.e. a number, for example A=4, then A⁻¹=1/4.
  Some facts about inverse matrices are:
  • Only square matrices can have an inverse, but not all square matrices have an inverse.
  • When the inverse exists, it is unique.
  • Finding the inverse of a large matrix A is numerically complicated (but computers do it for us).
  Finding the inverse of a matrix in R is done using the solve() function:

  > A <- matrix(c(1, 3, 2, 4), ncol=2,byrow=T)
  > A
   
  

       [,1] [,2]
  [1,]    1    3
  [2,]    2    4
   
  

  > #M2 <- matrix(c(-2,1.5,1,-0.5),ncol=2,byrow=T)
  > B <- solve(A)
  > B
   
  

       [,1] [,2]
  [1,]   -2  1.5
  [2,]    1 -0.5
   
  

  > A %*% B
   
  

       [,1] [,2]
  [1,]    1    0
  [2,]    0    1
   
  

  3.10  Solving systems of linear equations

  Matrices are closely related to systems of linear equations. Consider the two equations
  

  x1 + 3x2 = 7 2x1 + 4 x2 = 10

  The system can be written in matrix form
  

  ⎡ ⎢ ⎢ ⎣ 1 3 2 4 ⎤ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎣ x1 x2 ⎤ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎣ 7 10 ⎤ ⎥ ⎥ ⎦   i.e. Ax=b

  Since A⁻¹A=I and since Ix=x we have
  

  x = A−1b = ⎡ ⎢ ⎢ ⎣ −2 1.5 1 −0.5 ⎤ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎣ 7 10 ⎤ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎣ 1 2 ⎤ ⎥ ⎥ ⎦

  A geometrical approach to solving these equations is as follows: Isolate x₂ in the equations:
  

  x2 = 7 3 − 1 3 x1    x2 = 1 0 4 − 2 4 x1

  These two lines are shown in Figure 5 from which it can be seen that the solution is x₁=1, x₂=2.
  

  Figure 5: Solving two equations with two unknowns.
  From the Figure it follows that there are 3 possible cases of solutions to the system
  1. Exactly one solution - when the lines intersect in one point
  2. No solutions - when the lines are parallel but not identical
  3. Infinitely many solutions - when the lines coincide.

  > A <- matrix(c(1, 2, 3, 4), ncol=2)
  > b <- c(7, 10)
  > x <- solve(A) %*% b
  > x
   
  

       [,1]
  [1,]    1
  [2,]    2
   
  

3.11  Some additional rules for matrix operations

3.12  Details on inverse matrices*

3.12.1  Inverse of a 2×2 matrix*

3.12.2  Inverse of diagonal matrices*

3.12.3  Generalized inverse*

> A <- matrix(c(1, 2, 3, 2, 3, 4, 3, 5, 7), ncol=3)
> A # 3rd column is sum of the two first; the inverse does not exist
 

     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    2    3    5
[3,]    3    4    7
 

> det(A) 
 

[1] 0
 

> G <- MASS::ginv(A)
> A %*% G        # Not identity
 

           [,1]      [,2]       [,3]
[1,]  0.8333333 0.3333333 -0.1666667
[2,]  0.3333333 0.3333333  0.3333333
[3,] -0.1666667 0.3333333  0.8333333
 

> A %*% G %*% A  # This is A
 

     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    2    3    5
[3,]    3    4    7
 

3.12.4  Inverting an n×n matrix*

> A <- matrix(c(2, 2, 3, 3, 5, 9, 5, 6, 7), ncol=3)
> A
 

     [,1] [,2] [,3]
[1,]    2    3    5
[2,]    2    5    6
[3,]    3    9    7
 

> AB <- cbind(A, diag(c(1, 1, 1)))
> AB
 

     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    2    3    5    1    0    0
[2,]    2    5    6    0    1    0
[3,]    3    9    7    0    0    1
 

1. Multiply a row by a (non-zero) constant.
2. Multiply a row by a (non-zero) constant and add the result to another row.
3. Interchange two rows.

• First, we make sure that AB[1,1]=1. Then we subtract a constant times the first row from the second to obtain that AB[2,1]=0, and similarly for the third row:

  > AB[1,] <- AB[1,] / AB[1,1]
  > AB[2,] <- AB[2,] - 2 * AB[1,]
  > AB[3,] <- AB[3,] - 3 * AB[1,]
  > AB
   
  

       [,1] [,2] [,3] [,4] [,5] [,6]
  [1,]    1  1.5  2.5  0.5    0    0
  [2,]    0  2.0  1.0 -1.0    1    0
  [3,]    0  4.5 -0.5 -1.5    0    1
   
  
• Next we ensure that AB[2,2]=1. Afterwards we subtract a constant times the second row from the third to obtain that AB[3,2]=0:

  > AB[2,] <- AB[2,] / AB[2,2]
  > AB[3,] <- AB[3,] - 4.5 * AB[2,]
   
  
• Now we rescale the third row such that AB[3,3]=1:

  > AB[3,] <- AB[3,] / AB[3,3]
  > AB
   
  

       [,1] [,2] [,3]       [,4]      [,5]       [,6]
  [1,]    1  1.5  2.5  0.5000000 0.0000000  0.0000000
  [2,]    0  1.0  0.5 -0.5000000 0.5000000  0.0000000
  [3,]    0  0.0  1.0 -0.2727273 0.8181818 -0.3636364
   
  

  Then AB has zeros below the main diagonal.
• We then work our way up to obtain that AB has zeros above the main diagonal:

  > AB[2,] <- AB[2,] - 0.5 * AB[3,]
  > AB[1,] <- AB[1,] - 2.5 * AB[3,]
  > AB
   
  

       [,1] [,2] [,3]       [,4]        [,5]       [,6]
  [1,]    1  1.5    0  1.1818182 -2.04545455  0.9090909
  [2,]    0  1.0    0 -0.3636364  0.09090909  0.1818182
  [3,]    0  0.0    1 -0.2727273  0.81818182 -0.3636364
   
  

  > AB[1,] <- AB[1,] - 1.5 * AB[2,]
  > AB
   
  

       [,1] [,2] [,3]       [,4]        [,5]       [,6]
  [1,]    1    0    0  1.7272727 -2.18181818  0.6363636
  [2,]    0    1    0 -0.3636364  0.09090909  0.1818182
  [3,]    0    0    1 -0.2727273  0.81818182 -0.3636364
   
  

> B <- AB[,4:6]
> A%*% B
 

              [,1]         [,2]         [,3]
[1,]  1.000000e+00 3.330669e-16 1.110223e-16
[2,] -4.440892e-16 1.000000e+00 2.220446e-16
[3,] -2.220446e-16 9.992007e-16 1.000000e+00
 

4  Least squares

Figure 6: Regression

> y
 

[1] 3.7 4.2 4.9 5.7 6.0
 

> X
 

       x
[1,] 1 1
[2,] 1 2
[3,] 1 3
[4,] 1 4
[5,] 1 5
 

> beta.hat <- solve(t(X) %*% X) %*% t(X) %*% y
> beta.hat
 

  [,1]
  3.07
x 0.61
 

> as.numeric(X %*% beta.hat)
 

[1] 3.68 4.29 4.90 5.51 6.12
 

4.1  Least squares with generalized inverse

> X2 <- cbind(X, rowSums(X))
> X2
 

       x  
[1,] 1 1 2
[2,] 1 2 3
[3,] 1 3 4
[4,] 1 4 5
[5,] 1 5 6
 

> G <- MASS::ginv(t(X2) %*% X2)
> G
 

           [,1]       [,2]        [,3]
[1,]  0.6333333 -0.4333333  0.20000000
[2,] -0.4333333  0.3000000 -0.13333333
[3,]  0.2000000 -0.1333333  0.06666667
 

> beta2.hat <- G %*% t(X2) %*% y
> beta2.hat
 

           [,1]
[1,]  1.8433333
[2,] -0.6166667
[3,]  1.2266667
 

> beta22.hat <- c(beta.hat, 0)
 

> as.numeric(X2 %*% beta2.hat)
 

[1] 3.68 4.29 4.90 5.51 6.12
 

> as.numeric(X2 %*% beta22.hat)
 

[1] 3.68 4.29 4.90 5.51 6.12
 

5  A neat little exercise - from a bird's perspective

• Write up two equations with two unknowns.
• Solve these equations using the methods you have learned from linear algebra.
• Simply finding the solution by trial-and-error is considered cheating.