MATRIX ALGEBRA

In this lesson we’re going to describe the language of linear model which is called “linear algebra”. Later, we’ll see how this notation is useful to describe models and perform computation on these models, to code them up, to describe them in papers.

VECTORS In R, write down a series of numbers as vectors:

#install.packages("UsingR")
#library(UsingR)
data("father.son",package="UsingR")
summary(father.son)

    fheight         sheight     
 Min.   :59.01   Min.   :58.51  
 1st Qu.:65.79   1st Qu.:66.93  
 Median :67.77   Median :68.62  
 Mean   :67.69   Mean   :68.68  
 3rd Qu.:69.60   3rd Qu.:70.47  
 Max.   :75.43   Max.   :78.36  

x<-father.son$fheight #from the dataset father.son we create the vector x that is the height of the sons.
length(x) #we can count the entrees of the vector with the function length().

[1] 1078

In mathematical notation a vector can be expressed as follows:

\[Y=\begin{pmatrix}y_{1}\\ y_{2} \\ \vdots\\y_{n}\end{pmatrix}\]

capital \(Y\) represents the vector. And it holds in the \(n\) observations \(y_{1}\) through \(y_{n}\). We can rewrite these vectors for the predictor of covariates as well.

So in the example of dropping object we showed in the intro section this would be the time on \(X_{1}\) and then time squared in \(X_{2}\).

\[\begin{align*}X_{1}&=\begin{pmatrix}x_{1,1}\\ x_{2,1} \\ \vdots\\x_{N,1}\end{pmatrix} &\text{and}& & X_{2}&=\begin{pmatrix}x_{1,2}\\ x_{2,2} \\ \vdots\\x_{N,2}\end{pmatrix}\end{align*}\]

It is useful to put these vectors join them together into what we call matrices. Matrices, we can thing of them as a series of vectors one next to each other.So here we can see a matrix that has \(n\) rows and 2 columns:

\[X=[X_{1}X_{2}]=\begin{pmatrix}x_{1,1}& x_{1,2}\\ x_{2,1} & x_{2,2} \\ \vdots & \vdots\\x_{N,1} &x_{N,2}\end{pmatrix}\]

The number of rows and columns define its dimensions. This is an \(N\times 2\) matrix.

n<- 25
t<- seq(0,3.4,len=n) #time in secs, t is a base function. function seq()  generates a sequence from 0 to 3.4
X<- cbind(X1=t,X2=t^2) #this functions combine the time variable and the time squared variable
X

             X1          X2
 [1,] 0.0000000  0.00000000
 [2,] 0.1416667  0.02006944
 [3,] 0.2833333  0.08027778
 [4,] 0.4250000  0.18062500
 [5,] 0.5666667  0.32111111
 [6,] 0.7083333  0.50173611
 [7,] 0.8500000  0.72250000
 [8,] 0.9916667  0.98340278
 [9,] 1.1333333  1.28444444
[10,] 1.2750000  1.62562500
[11,] 1.4166667  2.00694444
[12,] 1.5583333  2.42840278
[13,] 1.7000000  2.89000000
[14,] 1.8416667  3.39173611
[15,] 1.9833333  3.93361111
[16,] 2.1250000  4.51562500
[17,] 2.2666667  5.13777778
[18,] 2.4083333  5.80006944
[19,] 2.5500000  6.50250000
[20,] 2.6916667  7.24506944
[21,] 2.8333333  8.02777778
[22,] 2.9750000  8.85062500
[23,] 3.1166667  9.71361111
[24,] 3.2583333 10.61673611
[25,] 3.4000000 11.56000000

dim(X)

[1] 25  2

Matrices can have any number of rows and columns, so in general we write \(X\) to represent any matrix. For the general linear model we use:

\[X=\begin{pmatrix}x_{1,1}&...&x_{1,p}\\ x_{2,1}&...&x_{2,p}\\ \vdots &\vdots & \vdots\\ x_{N,1}&...&x_{N,p} \end{pmatrix}\]

\(X\) is an \(N\times p\) matrix.

#In R we have vectors and matrices. You can create your own vectors with the function c().
vector<-c(1,2,3)
#Vectors are also the output of many functions such as:
rnorm(10)
#You can turn vectors into matrices using functions such as rbind(), cbind() or matrix().

Matrix Notation Exercises #1:

Create the matrix from the vector 1:1000 like this:

#What is the entry in row 25, column 3?
X[25,3]

[1] 225

Matrix Notation Exercises #2

Using the function cbind(), create a 10 x 5 matrix with first column x=1:10. Then use 2x, 3x, 4x and 5x respectively in columns 2 through 5.

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

[1] 105

sum(M[7,])

[1] 105

Before we go into the use of matrix algebra to solve and operate on linear models, we have to learn a few operations. Linear algebra was developed to solve a system of equations, like the one I’m showing here. Solve this system of equations: \[\begin{align}a+b+c&=6\\ 3a-2b+c&=2\\ 2a+b-c&=1\end{align}\] Matrix algebra notation:

\[\begin{pmatrix}1 & 1 & 1\\ 3 & -2 & 1\\2 & 1 & -1\end{pmatrix}\begin{pmatrix}a\\b\\c\end{pmatrix}=\begin{pmatrix}6\\2\\1\end{pmatrix}\Longrightarrow\begin{pmatrix}a\\b\\c\end{pmatrix}=\begin{pmatrix}1 & 1 & 1\\ 3 & -2 & 1\\2 & 1 & -1\end{pmatrix}^{-1}\begin{pmatrix}6\\2\\1\end{pmatrix}\] it gives us the general solution to any system of equations.

Before going and solve systems of equations we’ll remind some basic operations of matrices:

• Multiplying by scalar.

If \(a\in\mathbb{R}\) is a scalar and \(X\) is a \(N\times N\) matrix, then \[aX=\begin{pmatrix}ax_{1,1}&...&ax_{1,N}\\ \vdots & \ddots & \vdots\\ ax_{N,1}&...&ax_{N,N}\end{pmatrix}\]

a<-2
X<-matrix(1:20,4,5)
mult<-a*X
print(X)

     [,1] [,2] [,3] [,4] [,5]
[1,]    1    5    9   13   17
[2,]    2    6   10   14   18
[3,]    3    7   11   15   19
[4,]    4    8   12   16   20

print(mult)

     [,1] [,2] [,3] [,4] [,5]
[1,]    2   10   18   26   34
[2,]    4   12   20   28   36
[3,]    6   14   22   30   38
[4,]    8   16   24   32   40

• The Transpose.

The transpose is an operation that simply change the columns to rows: \[\begin{align}X=\begin{pmatrix}x_{1,1}&...&x_{1,p}\\ \vdots & \ddots & \vdots\\ x_{N,1}&...&x_{N,p}\end{pmatrix} & & &\Longleftrightarrow& X^{T}=\begin{pmatrix}x_{1,1} & ... & x_{N,1}\\ \vdots & \ddots &\vdots\\x_{1,p}&...&x_{N,p} \end{pmatrix} \end{align}\]

To do so, we only take the columns of the matrix X and make them rows.

X<-matrix(1:24,6,4)
print(X)

     [,1] [,2] [,3] [,4]
[1,]    1    7   13   19
[2,]    2    8   14   20
[3,]    3    9   15   21
[4,]    4   10   16   22
[5,]    5   11   17   23
[6,]    6   12   18   24

t(X)

     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    2    3    4    5    6
[2,]    7    8    9   10   11   12
[3,]   13   14   15   16   17   18
[4,]   19   20   21   22   23   24

Let’s back to the system of equations we described before: \[\begin{align}a+b+c&=6\\ 3a-2b+c&=2\\ 2a+b-c&=1\end{align}\] - Matrix multiplication The multiplication of the matrix works as follows: Each entry in the matrix rows are multiplied by each entry in the second matrix column: \[\begin{pmatrix}1 & 1 & 1\\ 3 & -2 & 1\\2 & 1 & -1\end{pmatrix}\begin{pmatrix}a\\b\\c\end{pmatrix}=\begin{pmatrix}a+b+c\\3a-2b+c\\2a+b-c\end{pmatrix}=\begin{pmatrix}6\\ 2\\ 1\end{pmatrix}\]

multiplication

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

For now let’s define general matrix multiplication:

\[\begin{align}AX&=\begin{pmatrix}a_{11} & ... & a_{1,N}\\ \vdots & \ddots & \vdots\\ a_{M,1} & ... & a_{M,N}\end{pmatrix} \begin{pmatrix}x_{1,1} &...& x_{1,p}\\ \vdots &\ddots & \vdots\\ x_{N,1}&...&x_{N,p}\end{pmatrix}\\ \\ &=\begin{pmatrix}\sum_{i=1}^{N}a_{1,i}x_{i,1} &...& \sum_{i=1}^{N}a_{1,i}x_{i,p}\\ &\vdots&\\ \sum_{i=1} ^{N}a_{M,i}x_{i,1} &...& \sum_{i=1}^{N}a_{M,i}x_{i,p} \end{pmatrix}\end{align}\] Note that the number of columns of first must match with the number of rows of the second.

• Identity Matrix The identity matrix (\(I\))has the property that \(AI=A\), it has this form: \[I=\begin{pmatrix}1&0&0&0&...&0&0\\0&1&0&0&...&0&0\\ 0&0&1&0&...&0&0\\ \vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\0&0&0&0&...&1&0\\0&0&0&0&...&0&1\end{pmatrix}\]

diag(5) #function diag(dimension) generates an identity matrix of any dimension. Note that the identity matrix must be square

-The Inverse

The inverse is very useful. If you remember the system of equations that we showed before, the solution can be written as the inverse of a matrix times a vector of numbers. That was going to give us a solution to what a, b, and c are. Some square matrices have inverses, not all of them. The inverse of matrix \(X\) is denoted by \(X^{-1}\). And \(XX^{-1}=I\).

Let’s back to our system of equations:

\[\begin{align}a+b+c&=6\\ 3a-2b+c&=2\\ 2a+b-c&=1\end{align}\] The solution requires to invert the matrix: \[\begin{pmatrix}a\\b\\c\end{pmatrix}=\begin{pmatrix}1 & 1 & 1\\ 3 & -2 & 1\\2 & 1 & -1\end{pmatrix}^{-1}\begin{pmatrix}6\\2\\1\end{pmatrix}\]

X

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

Y<-matrix(c(6,2,1),3,1)
Y

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

solve(X)%*%Y

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

Matrix Operation Exercises #2

Solve the following system of equations using R:

\[\begin{align} 3a+4b-5c+d&=10\\2a+2b+2c-d&=5\\a-b+5c-5d&=7\\5a+d&=4\end{align}\] \[\begin{align}\begin{pmatrix} 3&4&-5&1\\ 2&2&2&-1\\ 1&-1&5&-5\\5&0&0&1\end{pmatrix}\begin{pmatrix}a\\b\\c\\d\end{pmatrix}&=\begin{pmatrix}10\\5\\7\\4\end{pmatrix}\\ \\ \Leftrightarrow \underbrace{\begin{pmatrix} 3&4&-5&1\\ 2&2&2&-1\\ 1&-1&5&-5\\5&0&0&1\end{pmatrix}^{-1}}_{A^{-1}}\underbrace{\begin{pmatrix}10\\5\\7\\4\end{pmatrix}}_{B}&=\begin{pmatrix}a\\b\\c\\d\end{pmatrix}\end{align}\]

A<-matrix(c(3,2,1,5,4,2,-1,0,-5,2,5,0,1,-1,-5,1),4,4)
A

     [,1] [,2] [,3] [,4]
[1,]    3    4   -5    1
[2,]    2    2    2   -1
[3,]    1   -1    5   -5
[4,]    5    0    0    1

B<-matrix(c(10,5,7,4),4,1)
X<-solve(A)%*%B
name<-c("a","b","c","d")
X<-data.frame(name,X)
X

a <- matrix(1:12, nrow=4)
b <- matrix(1:15, nrow=3)

Matrix Operation Exercises #3

What is the value in the 3rd row and the 2nd column of the matrix product of \(a\) and \(b\)?

N<-a%*%b
N

     [,1] [,2] [,3] [,4] [,5]
[1,]   38   83  128  173  218
[2,]   44   98  152  206  260
[3,]   50  113  176  239  302
[4,]   56  128  200  272  344

N[3,2]

[1] 113

Matrix Operation Exercises #4

Multiply the 3rd row of \(a\) with the 2nd column of \(b\), using the element-wise vector multiplication with *.

What is the sum of the elements in the resulting vector?

sum(a[3,]*b[,2])

[1] 113