Biostatistical Computing, PHC 6068

Matrix operation

Zhiguang Huo (Caleb)

Wednesday Oct 07, 2020

Outline

Algebra notation

\[ \textbf{x} = (x_1, x_2, \ldots, x_p)^\top \]

\[ \textbf{X}= \begin{pmatrix} x_{11} & x_{12} & \ldots & x_{1p} \\ x_{21} & x_{22} & \ldots & x_{2p} \\ \ldots \\ x_{n1} & x_{n2} & \ldots & x_{np} \end{pmatrix} \]

Cross product of two vectors

\[ \textbf{a}^\top \textbf{b} = \textbf{b}^\top \textbf{a} = a_1b_1 + a_2b_2 + \ldots + a_p b_p = \sum_{i=1}^p a_ib_i \]

Matrix notation

\[ \textbf{a}^\top\textbf{X} = \textbf{a}^\top (\textbf{X}_1, \textbf{X}_2, \ldots, \textbf{X}_p) = (\textbf{a}^\top \textbf{X}_1, \textbf{a}^\top\textbf{X}_2, \ldots, \textbf{a}^\top\textbf{X}_p) \in \mathbb{R}^{1 \times p} \]

Matrix notation

\[\textbf{X} \textbf{Y} = \begin{pmatrix} \sum_{i=1}^p x_{1i}y_{i1} & \sum_{i=1}^p x_{1i}y_{i2} & \ldots & \sum_{i=1}^p x_{1i}y_{im} \\ \sum_{i=1}^p x_{2i}y_{i1} & \sum_{i=1}^p x_{2i}y_{i2} & \ldots & \sum_{i=1}^p x_{2i}y_{im} \\ \ldots \\ \sum_{i=1}^p x_{ni}y_{i1} & \sum_{i=1}^p x_{ni}y_{i2} & \ldots & \sum_{i=1}^p x_{ni}y_{im} \end{pmatrix} \]

Matrix notation

\[\textbf{X}= \begin{pmatrix} \textbf{x}_{1}^\top\\ \ldots \\ \textbf{x}_{n}^\top \end{pmatrix} \]

\[\textbf{Y}= (\textbf{y}_{1}, \ldots, \textbf{y}_{m}) \]

\[\textbf{X} \textbf{Y} = \begin{pmatrix} \textbf{x}_{1}^\top \textbf{y}_{1} & \textbf{x}_{1}^\top \textbf{y}_{2} & \ldots & \textbf{x}_{1}^\top \textbf{y}_{m} \\ \textbf{x}_{2}^\top \textbf{y}_{1} & \textbf{x}_{2}^\top \textbf{y}_{2} & \ldots & \textbf{x}_{2}^\top \textbf{y}_{m} \\ \ldots \\ \textbf{x}_{n}^\top \textbf{y}_{1} & \textbf{x}_{n}^\top \textbf{y}_{2} & \ldots & \textbf{x}_{n}^\top \textbf{y}_{m} \\ \end{pmatrix} \]

Inverse of a matrix

\[\textbf{Y}^{-1} \textbf{Y} = I\] \[\textbf{Y} \textbf{Y}^{-1} = I\]

\[\textbf{Y}^{-1} = \frac{I}{\textbf{Y}}\]

Create a matrix

mdat <- matrix(c(1,2,3, 11,12,13), nrow = 2, ncol = 3, byrow = FALSE,
               dimnames = list(c("row1", "row2"),
                               c("C.1", "C.2", "C.3")))
mdat
##      C.1 C.2 C.3
## row1   1   3  12
## row2   2  11  13
mdat <- c(1,2,3, 11,12,13)
dim(mdat) <- c(2,3)
rownames(mdat) <- c("row1", "row2")
colnames(mdat) <- c("C.1", "C.2", "C.3")
mdat
##      C.1 C.2 C.3
## row1   1   3  12
## row2   2  11  13

Appending a column or a row to a matrix

dim(mdat)
## [1] 2 3
mat2 <- cbind(mdat, c(3,4))
dim(mat2)
## [1] 2 4
mat3 <- rbind(mdat, c(2,5,8))
dim(mat3)
## [1] 3 3

Subset of a matrix

mdat[1,1]
## [1] 1
mdat[1,1:2]
## C.1 C.2 
##   1   3
mdat[1,]
## C.1 C.2 C.3 
##   1   3  12
mdat[1,-1]
## C.2 C.3 
##   3  12
mdat[,]
##      C.1 C.2 C.3
## row1   1   3  12
## row2   2  11  13
dimnames(mdat); rownames(mdat); colnames(mdat)
## [[1]]
## [1] "row1" "row2"
## 
## [[2]]
## [1] "C.1" "C.2" "C.3"
## [1] "row1" "row2"
## [1] "C.1" "C.2" "C.3"
mdat["row1", "C.1"]
## [1] 1
mdat["row2", c("C.2", "C.3")]
## C.2 C.3 
##  11  13
mdat["row1",-2] ## mixed usage
## C.1 C.3 
##   1  12

Transponse of a Matrix

mat1 <- matrix(1:6,nrow=2,ncol=3)
mat1
##      [,1] [,2] [,3]
## [1,]    1    3    5
## [2,]    2    4    6
t(mat1)
##      [,1] [,2]
## [1,]    1    2
## [2,]    3    4
## [3,]    5    6

Element-wise Multiplication by a Scalar/Vector/Matrix

mdat * 3
##      C.1 C.2 C.3
## row1   3   9  36
## row2   6  33  39
mdat * c(1,2)
##      C.1 C.2 C.3
## row1   1   3  12
## row2   4  22  26
mdat * matrix(1:6,nrow=2,ncol=3)
##      C.1 C.2 C.3
## row1   1   9  60
## row2   4  44  78

Addition/Subtraction by a Scalar/Vector/Matrix

mdat + 3
##      C.1 C.2 C.3
## row1   4   6  15
## row2   5  14  16
mdat - c(1,2)
##      C.1 C.2 C.3
## row1   0   2  11
## row2   0   9  11
mdat + matrix(1:6,nrow=2,ncol=3)
##      C.1 C.2 C.3
## row1   2   6  17
## row2   4  15  19

Matrix Multiplication

mat1 <- matrix(1:6,nrow=2,ncol=3)
mat1
##      [,1] [,2] [,3]
## [1,]    1    3    5
## [2,]    2    4    6
mat2 <- matrix(1:6,nrow=3,ncol=2)
mat2
##      [,1] [,2]
## [1,]    1    4
## [2,]    2    5
## [3,]    3    6
mat1 %*% mat2
##      [,1] [,2]
## [1,]   22   49
## [2,]   28   64
mat2 %*% mat1
##      [,1] [,2] [,3]
## [1,]    9   19   29
## [2,]   12   26   40
## [3,]   15   33   51

crossProd

A <- matrix(c(1,2), ncol =1)
B <- matrix(c(2,2), ncol=1)

A 
##      [,1]
## [1,]    1
## [2,]    2
B
##      [,1]
## [1,]    2
## [2,]    2
crossprod(A,B)
##      [,1]
## [1,]    6
crossprod(t(A),t(B))
##      [,1] [,2]
## [1,]    2    2
## [2,]    4    4
crossprod(A)
##      [,1]
## [1,]    5
crossprod(t(A))
##      [,1] [,2]
## [1,]    1    2
## [2,]    2    4

Diagonal Matrix

S <- matrix(1:9,3,3)
S
##      [,1] [,2] [,3]
## [1,]    1    4    7
## [2,]    2    5    8
## [3,]    3    6    9
diag(S)
## [1] 1 5 9
diag(c(1,4,7))
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    4    0
## [3,]    0    0    7
diag(3)
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1

Inverse of a Matrix

S <- matrix(c(4,4,-2,2,6,2,2,8,4),3,3)
SI <- solve(S)
SI
##      [,1] [,2] [,3]
## [1,]  1.0 -0.5  0.5
## [2,] -4.0  2.5 -3.0
## [3,]  2.5 -1.5  2.0
S %*% SI
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1
SI %*% S
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1

More about solve

set.seed(32608)
(A <- matrix(rnorm(9),3,3))
##             [,1]       [,2]       [,3]
## [1,]  0.07699352  0.3381365 -1.9363229
## [2,] -0.56610906 -0.3786229 -0.5720703
## [3,]  0.80627284  0.7446219  0.7198154
(b <- runif(3))
## [1] 0.5756702 0.8367205 0.2964995
(x <- solve(A,b))
## [1] -7.2976775  7.5871849  0.7374597
A %*% x
##           [,1]
## [1,] 0.5756702
## [2,] 0.8367205
## [3,] 0.2964995

Moore-Penrose Generalized Inverse

library(MASS)
set.seed(32608)
A <- matrix(rnorm(9),3,3)
ginv(A)
##            [,1]      [,2]       [,3]
## [1,]  0.7013869 -7.703428 -4.2355174
## [2,] -0.2457031  7.389853  5.2121071
## [3,] -0.5314604  0.984167  0.7417649
(A2 <- A[,1:2])
##             [,1]       [,2]
## [1,]  0.07699352  0.3381365
## [2,] -0.56610906 -0.3786229
## [3,]  0.80627284  0.7446219
ginv(A2)
##           [,1]      [,2]      [,3]
## [1,] -2.572491 -1.640808 0.3338684
## [2,]  3.079399  1.232375 0.5712269
A2 %*% ginv(A2) %*% A2
##             [,1]       [,2]
## [1,]  0.07699352  0.3381365
## [2,] -0.56610906 -0.3786229
## [3,]  0.80627284  0.7446219

Eigen value decomposition

A0 <- c(3,1,1,
       1,3,2,
       1,2,3)
A <- matrix(A0,3,3)
eigen(A)
## eigen() decomposition
## $values
## [1] 5.732051 2.267949 1.000000
## 
## $vectors
##            [,1]       [,2]       [,3]
## [1,] -0.4597008  0.8880738  0.0000000
## [2,] -0.6279630 -0.3250576 -0.7071068
## [3,] -0.6279630 -0.3250576  0.7071068
Q <- eigen(A)$vectors
lambda <- eigen(A)$values

Q %*% t(Q)
##               [,1]          [,2]          [,3]
## [1,]  1.000000e+00 -5.551115e-17  0.000000e+00
## [2,] -5.551115e-17  1.000000e+00 -2.220446e-16
## [3,]  0.000000e+00 -2.220446e-16  1.000000e+00
Q %*% diag(lambda) %*% ginv(Q)
##      [,1] [,2] [,3]
## [1,]    3    1    1
## [2,]    1    3    2
## [3,]    1    2    3
A %*% Q[,1]
##           [,1]
## [1,] -2.635029
## [2,] -3.599516
## [3,] -3.599516
lambda[1] * Q[,1]
## [1] -2.635029 -3.599516 -3.599516
A %*% Q[,2]
##            [,1]
## [1,]  2.0141063
## [2,] -0.7372141
## [3,] -0.7372141
lambda[2] * Q[,2]
## [1]  2.0141063 -0.7372141 -0.7372141

Singular value decomposition (SVD)

M <- matrix(1:6,2,3) 
M
##      [,1] [,2] [,3]
## [1,]    1    3    5
## [2,]    2    4    6
svd(M)
## $d
## [1] 9.5255181 0.5143006
## 
## $u
##            [,1]       [,2]
## [1,] -0.6196295 -0.7848945
## [2,] -0.7848945  0.6196295
## 
## $v
##            [,1]       [,2]
## [1,] -0.2298477  0.8834610
## [2,] -0.5247448  0.2407825
## [3,] -0.8196419 -0.4018960
svdRes <- svd(M)
svdRes$u %*% diag(svdRes$d) %*% t(svdRes$v)
##      [,1] [,2] [,3]
## [1,]    1    3    5
## [2,]    2    4    6
t(svdRes$u) %*% svdRes$u
##              [,1]         [,2]
## [1,] 1.000000e+00 2.220446e-16
## [2,] 2.220446e-16 1.000000e+00
t(svdRes$v) %*% svdRes$v
##              [,1]         [,2]
## [1,] 1.000000e+00 1.110223e-16
## [2,] 1.110223e-16 1.000000e+00

Determinant and trace of a Matrix

S <- matrix(c(4,4,-2,2,6,2,2,8,4),3,3)
det(S)
## [1] 8
S <- matrix(c(4,4,-2,2,6,2,2,8,4),3,3)
sum(diag(S))
## [1] 14

Rank of a Matrix

S <- matrix(c(4,4,-2,2,6,2,2,8,4),3,3)
QR <- qr(S)
QR$rank
## [1] 3