Linear algebra

Matrix basics

  • dimensions of matrix
   seed=123
   X=matrix(1:12,ncol=3,nrow=4)
   X
##      [,1] [,2] [,3]
## [1,]    1    5    9
## [2,]    2    6   10
## [3,]    3    7   11
## [4,]    4    8   12
   dim(X)
## [1] 4 3
   dim(X)[1]
## [1] 4
   ncol(X)
## [1] 3
  • change dimensions of matrix
dim(X)=c(2,6)
   X
##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    1    3    5    7    9   11
## [2,]    2    4    6    8   10   12
  • change names of matrix
  a <- matrix(1:6,nrow=2,ncol=3,byrow=FALSE)
  b <- matrix(1:6,nrow=3,ncol=2,byrow=T)
  c <- matrix(1:6,nrow=3,ncol=2,byrow=T,dimnames=list(c("A","B","C"),c("boy","girl")))
  c
##   boy girl
## A   1    2
## B   3    4
## C   5    6
   rownames(c)
## [1] "A" "B" "C"
   colnames(c)
## [1] "boy"  "girl"
   rownames(c)=c("E","F","G")
   c
##   boy girl
## E   1    2
## F   3    4
## G   5    6
  • replace elements of matrix
   X=matrix(1:12,nrow=4,ncol=3)
   
   X[2,3]
## [1] 10
   X[2,3]=1000
      
   X
##      [,1] [,2] [,3]
## [1,]    1    5    9
## [2,]    2    6 1000
## [3,]    3    7   11
## [4,]    4    8   12
  • extract diagonal elements and replace them
   diag(X)
## [1]  1  6 11
   diag(X)=c(0,0,1)

   X
##      [,1] [,2] [,3]
## [1,]    0    5    9
## [2,]    2    0 1000
## [3,]    3    7    1
## [4,]    4    8   12
  • create diagonal identity matrix
  diag(c(1,2,3))
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    2    0
## [3,]    0    0    3
  diag(3)
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1
  • extract lower/upper triangle elements
  X
##      [,1] [,2] [,3]
## [1,]    0    5    9
## [2,]    2    0 1000
## [3,]    3    7    1
## [4,]    4    8   12
  X[lower.tri(X)]
## [1]  2  3  4  7  8 12
  X[upper.tri(X)]
## [1]    5    9 1000
  • create lower/upper triangle matrix
X[lower.tri(X)]=0
X
##      [,1] [,2] [,3]
## [1,]    0    5    9
## [2,]    0    0 1000
## [3,]    0    0    1
## [4,]    0    0    0

Operations

  • transform
 t(X)
##      [,1] [,2] [,3] [,4]
## [1,]    0    0    0    0
## [2,]    5    0    0    0
## [3,]    9 1000    1    0
  • summary by row or column
  A=matrix(1:12,3,4)
  A 
##      [,1] [,2] [,3] [,4]
## [1,]    1    4    7   10
## [2,]    2    5    8   11
## [3,]    3    6    9   12
  rowSums(A)
## [1] 22 26 30
  rowMeans(A)
## [1] 5.5 6.5 7.5
  • determinant of matrix
  X=matrix(rnorm(9),nrow=3,ncol=3)
  det(X)
## [1] -1.284583
  • Addition
  A=matrix(1:12,nrow=3,ncol=4)
  B=matrix(1:12,nrow=3,ncol=4)
  A+B #same dimensions (non-conformable arrays)
##      [,1] [,2] [,3] [,4]
## [1,]    2    8   14   20
## [2,]    4   10   16   22
## [3,]    6   12   18   24
  • addition by scale
A=matrix(1:12,nrow=3,ncol=4)
2+A
##      [,1] [,2] [,3] [,4]
## [1,]    3    6    9   12
## [2,]    4    7   10   13
## [3,]    5    8   11   14
  • multiple by scale
A=matrix(1:12,nrow=3,ncol=4)
2*A
##      [,1] [,2] [,3] [,4]
## [1,]    2    8   14   20
## [2,]    4   10   16   22
## [3,]    6   12   18   24
  • dot multiple
   A=matrix(1:12,nrow=2,ncol=4)
## Warning in matrix(1:12, nrow = 2, ncol = 4): data length differs from size of
## matrix: [12 != 2 x 4]
   B=matrix(1:12,nrow=4,ncol=3)
   
   A%*%B
##      [,1] [,2] [,3]
## [1,]   50  114  178
## [2,]   60  140  220
  • kronecker multiple
   A=matrix(1:4,2,2)
   B=matrix(rep(1,4),2,2)
   
   kronecker(A,B)
##      [,1] [,2] [,3] [,4]
## [1,]    1    1    3    3
## [2,]    1    1    3    3
## [3,]    2    2    4    4
## [4,]    2    2    4    4
  • inverse matrix

must be square

   A=matrix(rnorm(9),nrow=3,ncol=3) 
   A
##            [,1]      [,2]       [,3]
## [1,] 0.26588723 1.4689811 -0.8353634
## [2,] 0.10596427 0.4447248  0.7555293
## [3,] 0.02644485 0.3020064  1.4258913
   AI=solve(A)
   AI
##            [,1]       [,2]        [,3]
## [1,] -3.9966334 23.1052182 -14.5840782
## [2,]  1.2908175 -3.9499953   2.8491925
## [3,] -0.1992752  0.4081026   0.3683305
   # identity matrix    
   (A%*%AI)
##               [,1]          [,2]          [,3]
## [1,]  1.000000e+00  0.000000e+00 -1.110223e-16
## [2,] -1.110223e-16  1.000000e+00  0.000000e+00
## [3,]  0.000000e+00 -2.220446e-16  1.000000e+00
library(matlib)
inv(A)
##            [,1]       [,2]        [,3]
## [1,] -3.9966334 23.1052182 -14.5840782
## [2,]  1.2908175 -3.9499953   2.8491925
## [3,] -0.1992752  0.4081026   0.3683305
  • generalized inverse when it is not square
   library(MASS)
   A2=matrix(1:12,nrow=3,ncol=4)
   
   A%*%ginv(A)
##              [,1]          [,2]         [,3]
## [1,] 1.000000e+00 -3.996803e-15 2.831069e-15
## [2,] 6.383782e-16  1.000000e+00 1.276756e-15
## [3,] 1.110223e-16 -7.771561e-16 1.000000e+00
   A%*%ginv(A)%*%A
##            [,1]      [,2]       [,3]
## [1,] 0.26588723 1.4689811 -0.8353634
## [2,] 0.10596427 0.4447248  0.7555293
## [3,] 0.02644485 0.3020064  1.4258913
  • crossprod
B=matrix(1:12,nrow=4,ncol=3)
B
##      [,1] [,2] [,3]
## [1,]    1    5    9
## [2,]    2    6   10
## [3,]    3    7   11
## [4,]    4    8   12
crossprod(B)
##      [,1] [,2] [,3]
## [1,]   30   70  110
## [2,]   70  174  278
## [3,]  110  278  446
t(B)%*%B
##      [,1] [,2] [,3]
## [1,]   30   70  110
## [2,]   70  174  278
## [3,]  110  278  446
   B=matrix(1:12,nrow=4,ncol=3)
   ginv(B)
##        [,1]        [,2]        [,3]    [,4]
## [1,] -0.375 -0.14583333  0.08333333  0.3125
## [2,] -0.100 -0.03333333  0.03333333  0.1000
## [3,]  0.175  0.07916667 -0.01666667 -0.1125
   ginv(crossprod(B))
##             [,1]        [,2]        [,3]
## [1,]  0.26649306  0.07638889 -0.11371528
## [2,]  0.07638889  0.02222222 -0.03194444
## [3,] -0.11371528 -0.03194444  0.04982639
   # solve(crossprod(B)) #is singular
   # solve((B)) #is not square
  • eigen values (decomposition) mxn matrix A=UΛU
   A=matrix(1:9,nrow=3,ncol=3)
   
   Aeigen=eigen(A)
   
   Aeigen$values
## [1]  1.611684e+01 -1.116844e+00 -5.700691e-16
   val <- diag(Aeigen$values)
   val
##          [,1]      [,2]          [,3]
## [1,] 16.11684  0.000000  0.000000e+00
## [2,]  0.00000 -1.116844  0.000000e+00
## [3,]  0.00000  0.000000 -5.700691e-16
  • eigen vectors
Aeigen$vectors
##            [,1]       [,2]       [,3]
## [1,] -0.4645473 -0.8829060  0.4082483
## [2,] -0.5707955 -0.2395204 -0.8164966
## [3,] -0.6770438  0.4038651  0.4082483
 round(Aeigen$vectors%*%val%*%t(Aeigen$vectors))
##      [,1] [,2] [,3]
## [1,]    3    4    5
## [2,]    4    5    6
## [3,]    5    6    7
 A
##      [,1] [,2] [,3]
## [1,]    1    4    7
## [2,]    2    5    8
## [3,]    3    6    9

Advanced operations

  • Choleskey factor

positive definite matrix (symmetric square), A=P’P

covariance matrix is semi positive definite matrix

   A=diag(3)+1
   A  
##      [,1] [,2] [,3]
## [1,]    2    1    1
## [2,]    1    2    1
## [3,]    1    1    2
   chol(A)
##          [,1]      [,2]      [,3]
## [1,] 1.414214 0.7071068 0.7071068
## [2,] 0.000000 1.2247449 0.4082483
## [3,] 0.000000 0.0000000 1.1547005
   t(chol(A))%*%chol(A)
##      [,1] [,2] [,3]
## [1,]    2    1    1
## [2,]    1    2    1
## [3,]    1    1    2

inverse using Choleshey

chol2inv(chol(A))
##       [,1]  [,2]  [,3]
## [1,]  0.75 -0.25 -0.25
## [2,] -0.25  0.75 -0.25
## [3,] -0.25 -0.25  0.75
inv(A)
##       [,1]  [,2]  [,3]
## [1,]  0.75 -0.25 -0.25
## [2,] -0.25  0.75 -0.25
## [3,] -0.25 -0.25  0.75
  • singular value (svd) decomposition

m×n matrix

A=UDV

   A=matrix(1:18,3,6)
   
   svd(A)   $d
## [1] 4.589453e+01 1.640705e+00 1.366522e-15
   t(svd(A)   $u)%*%svd(A)   $u
##              [,1]         [,2]         [,3]
## [1,] 1.000000e+00 3.330669e-16 1.665335e-16
## [2,] 3.330669e-16 1.000000e+00 5.551115e-17
## [3,] 1.665335e-16 5.551115e-17 1.000000e+00
   t(svd(A)   $v)%*%svd(A)   $v
##              [,1]          [,2]          [,3]
## [1,] 1.000000e+00  2.775558e-17  2.775558e-17
## [2,] 2.775558e-17  1.000000e+00 -2.081668e-16
## [3,] 2.775558e-17 -2.081668e-16  1.000000e+00
   svd(A)   $u %*%diag(svd(A)   $d)%*% t(svd(A)   $v)
##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    1    4    7   10   13   16
## [2,]    2    5    8   11   14   17
## [3,]    3    6    9   12   15   18
   A
##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    1    4    7   10   13   16
## [2,]    2    5    8   11   14   17
## [3,]    3    6    9   12   15   18
  • QR decomposition

m×n matrix

A=QR，where Q’Q=I, Q is orthogonal matrix

   A=matrix(1:12,4,3)
   qr(A)
## $qr
##            [,1]        [,2]          [,3]
## [1,] -5.4772256 -12.7801930 -2.008316e+01
## [2,]  0.3651484  -3.2659863 -6.531973e+00
## [3,]  0.5477226  -0.3781696  1.601186e-15
## [4,]  0.7302967  -0.9124744 -5.547002e-01
## 
## $rank
## [1] 2
## 
## $qraux
## [1] 1.182574 1.156135 1.832050
## 
## $pivot
## [1] 1 2 3
## 
## attr(,"class")
## [1] "qr"
qr.Q(qr(A))
##            [,1]          [,2]       [,3]
## [1,] -0.1825742 -8.164966e-01 -0.4000874
## [2,] -0.3651484 -4.082483e-01  0.2546329
## [3,] -0.5477226 -1.665335e-16  0.6909965
## [4,] -0.7302967  4.082483e-01 -0.5455419
qr.R(qr(A))
##           [,1]       [,2]          [,3]
## [1,] -5.477226 -12.780193 -2.008316e+01
## [2,]  0.000000  -3.265986 -6.531973e+00
## [3,]  0.000000   0.000000  1.601186e-15

Solve linear equations

   X=matrix(c(2, 2, 2, 0, 1, 1  ,2, 2, 0),ncol=3,nrow=3)
   X
##      [,1] [,2] [,3]
## [1,]    2    0    2
## [2,]    2    1    2
## [3,]    2    1    0
   b=1:3
   b
## [1] 1 2 3
   solve(X,b) # whether it is singular
## [1]  1.0  1.0 -0.5

Summary

  • Eigen values (values and vectors, UΛU), singular values (s v “" d), and QR (orthogonal and upper triangle, QR) for any matrix.

  • Choleskey (P’ P, lower and upper triangle) factor for positive definite matrix (cov matrix is semi).