B Linear algebra
Reading materials: STA108_LinearAlgebra_S20.pdf.
B.1 Vector
x = c(1,2,4,2.5,-1.7,1) # creates a vector (you can use "<-" instead of "=")
length(x) # check the length of x
class(x) # class of the object; you can also check typeof()
### Try z=list(); class(z);
x[ 3 ] # extract one element in x..
x[c(1,2,3)]
x[2:4]
z<-numeric(10)
### Transpose:
x.t=t(x);
length(x.t)
x.t
x
y = x^2
### Q: what is the length of y?
length(y)
### The "default" operation in R is elementwise:
x/y # elementwise division
y-x*x # all elementwise..
### There are other operators..
x %o% y # outer product of x and y
x %*% y # Matrix product. Q: what is R doing here?
sum(x*y)B.2 Matrix
M = matrix(0,2,3) # matrix with 2 rows and 3 columns with all entries 0
# ?matrix
Xmat1 = matrix(x,2,3) # 2 x 3 matrix formed by elements of vector x,
### ordered in *column-wise* manner
Xmat2 = matrix(x,2,3,byrow=T) # ordered row-wise
M<- Xmat1;
dimnames(M)
rownames(M) = c("first","second")
colnames(M) = c("a","b","c")
dimnames(M)
typeof(dimnames(M))
### access elements through names of rows/columns
M[,"a"]
M[,1]
M[,c("a","c")]
M[,c(1,3)]
M %*% c(4,-3,1) # multiply matrix M to vector c(4,-3,1)
M + c(1,-2) # add the column vector c(1,-2) to each column of M
M2 = cbind(M,c(4,5)) # create matrix M2 by appending column c(4,5) to M
### How to append row to M2?
new_row = c(1,3,5,7)
#?cbind
#?rbind
### What if we call length() on a matrix?
length(M)
dim(M)
### Transpose of the matrix
M.t=t(M);
identical(t(M.t),M)
### Diagonal and identity matrix
## Identity matrix
diag(5)
## Diagonal matrix
diag(c(1,2,3,4,5))
## Symmatric matrix
A= matrix(rnorm(9),3,3);
A=A+t(A); # To make sure it is symmetric
identical(A,t(A))
### Trace of a matrix:
sum(diag(A))
## Define a tr() function to match the common convention
tr<-function(x){
out=sum(diag(x));
return(out)
}
tr(A)
## Verify some properties of the trace
B= matrix(rnorm(9),3,3);
tr(A)
tr(t(A))
tr(A+B)
tr(A)+tr(B)
tr(A%*%B)
tr(B%*%A) # if the dimensions check out..
tr(A)
sum(eigen(A)$values)
### Matrix-vector product:
Xmat1 = matrix(rnorm(6),2,3)
xvec1=rnorm(3)
Xmat1%*%xvec1
### elementwise product of Xmat1 and Xmat2
Xmat1*Xmat2
### Matrix multiplication
# Xmat1%*%Xmat2
dim(Xmat1)
dim(Xmat2)
Xmat1%*%t(Xmat2)
### Verify another property of transpose:
t( Xmat1%*%t(Xmat2) )
Xmat2%*%t(Xmat1)
# solve(Xmat1) # what does the error message say?### Inverse of a marix
# set.seed(1)
# A= matrix(rnorm(9),3,3);
# set.seed(): specify the random seed for generating random numbers
# Because a computer (if not a quantum one) can only generate pseudorandomness
# set.seed(1)
# B=matrix(rnorm(9),3,3);
# identical(A,B)
set.seed(1)
A= matrix(rnorm(9),3,3);
A.inv=solve(A);
A.inv%*%A
A%*%A.inv
B=matrix(rnorm(9),3,3);
identical(A,B)
B.inv=solve(B);
AB.inv=solve(A%*%B)
B.inv%*%A.inv
AB.inv
max(abs(B.inv%*%A.inv-AB.inv))
# Why?
# This is because
# (A%*%B)%*%AB.int= I if AB.inv is the inverse of AB=A%*%B
# and that
# A%*%B %*%B.inv%*%A.inv=A%*%(B%*%B.inv)%*%A.inv =A%*%A.inv=I
### Projection matrix
set.seed(12)
X = matrix(rnorm(6),3,2)
P=X%*%solve(t(X)%*%X)%*%t(X)
P%*%X
X
(diag(3)-P)%*%X
### Rank of a matrix:
qr(X)$rank # Check out the qr() function yourself!
# The rank of X is 2 when the dimension of X is 3 by 2
# This means that the matrix X is full rank
## How about x.transpose times X?
XtX=t(X)%*%X
dim(XtX)
qr(XtX)$rank
# This is a full rank matrix
## How about X times X.transpose?
XXt=X%*%t(X)
dim(XXt)
qr(XXt)$rank
# This is NOT a full rank matrix
# Rank and inverse
solve(XtX)
# solve(XXt)
### Determinants
det(A)
1/det(solve(A))
det(A%*%B)
det(A)%*%det(B)
det(3*A)
A.symm = (A+t(A))/2;
det(A.symm)
A.symm.eigenvalues=eigen(A.symm)$values
prod(A.symm.eigenvalues)
# Trace: summation of eigenvalues
# Determinant: product of eigenvalues (for symmetric matrix only!)
### Eigenvalue and eigenvector
A.symm.eigenvectors=eigen(A.symm)$vectors
A.symm%*%A.symm.eigenvectors[,1]
A.symm.eigenvalues[1]*A.symm.eigenvectors[,1]
A.symm.eigenvalues
qr(A.symm)$rank
qr(XXt)$rank
eigen(XXt) # the last eigenvalue equals zero, with some numerical error
### Matrix definiteness
eigen(XtX)
eigen(P)
### Eigenvalue decomposition
A.symm.eigenvectors%*%diag(A.symm.eigenvalues)%*%t(A.symm.eigenvectors)
A.symm
### Singular value decomposition
svd(X) # singular value decomposition of XB.3 Other operations on vectors and matrices
### How does R handle mathematically invalid operations?
### Example 1:
x = c(1,2,4,2.5,-1.7,1)
w =1:3 # short-handed for c(1,2,3)
### Q: What is the length of w?
x + w # no error displayed, but not quite a mathematically valid operation
### Q: what did R just do???
c(x[1:3]+w, x[4:6]+w)
### What if...
# Create a vector of length 4 (from 1 to 4)
n<-1:4;
x+n
### Example 2:
z = sin(x) + x^{1.5}* exp(-x^3)
### What is (-1.7)^{1.5}???
### How did R handle this?
x[5]
z[5]B.3.1 Array
### Arrays, a.k.a. tensors
### A matrix is a two-way array
# ?array # check for R documentation on array
A = array(rnorm(2*3*4),dim=c(2,3,4)) # creates an array with entries that are i.i.d. N(0,1)
#?rnorm
dim(A) # checks the dimension of A
A[,1,] # returns the 2 x 4 matrix formed by restricting the second dimension index to 1
### What is the dimension of A[,1,]?
A[2,,] # returns the 3 x 4 matrix formed by restricting the first dimension index to 2
A[2,1,2]
A[A<0]B.3.2 Datasets as matrices
library(AER) # install.packages('AER')
data("Fatalities")
#?Fatalities
data("CreditCard")
dat.CC=CreditCard
#?CreditCard
dim(dat.CC)
head(dat.CC)
Y.CC=dat.CC[,1]=='yes';
X.CC=as.matrix(dat.CC[,c('reports','age','income','share','expenditure','dependents','months','majorcards','active')])
# Check the rank of the covariates:
qr(X.CC)$rank
dim(X.CC) # full rank, i.e., no multicolinearity
# Linear regression:
fit.CC=lm(Y.CC~X.CC-1);
summary(fit.CC)
# Fitted coefficients:
beta.hat=solve(t(X.CC)%*%X.CC)%*%t(X.CC)%*%Y.CC;
fit.CC$coefficients-beta.hat
# Hat matrix:
P.CC=X.CC%*%solve(t(X.CC)%*%X.CC)%*%t(X.CC)
# Fitted Y
Yhat.CC=P.CC%*%Y.CC
# Compare with the fits from lm()
max(abs(fit.CC$fitted.values-Yhat.CC))