chapter1 Matrix Algebra_2. etc...

건국대학교 응용통계학과 3-2 김형문 교수님의 다변수함수론의 R코드를 사용하였습니다.

 

Rank of a matrix

 

아래는 해당 내용 증명이다.

install.packages("matrixcalc")
library(matrixcalc)
A <- matrix(c(1,2,2,4),2,2); B <- diag(1:2)
A%*%B

     [,1] [,2]
[1,]    1    4
[2,]    2    8

matrix.rank(A%*%B)<=min(matrix.rank(A),matrix.rank(B))

[1] TRUE

matrix.rank(A+B)<=matrix.rank(A)+matrix.rank(B)

[1] TRUE

all.equal(matrix.rank(A),matrix.rank(t(A)),matrix.rank(A%*%t(A)),matrix.rank(A%*%t(A))) # 이 모든게 갖은가?

[1] TRUE

 

Trace of a matrix

library(matrixcalc)
A <- matrix(c(1,2,2,4),2,2); B <- diag(1:2)

matrix.trace(A%*%B) == matrix.trace(B%*%A)

[1] TRUE

해당 내용 증명은 아래를 참고하면 된다. (tr(AB) = tr(BA))

linear algebra - How to prove $\operatorname{Tr}(AB) = \operatorname{Tr}(BA)$? - Mathematics Stack Exchange

How to prove $\operatorname{Tr}(AB) = \operatorname{Tr}(BA)$?

 

Eigenvalues and Eigenvectors

library(matrixcalc)
A <- matrix(c(1,9,4,1),2,2) 
eigen(A)

eigen() decomposition
$values
[1]  7 -5

$vectors
          [,1]       [,2]
[1,] 0.5547002 -0.5547002
[2,] 0.8320503  0.8320503

# eigen value가 7 ,-5가 나온다는 말이고, 7에 해당되는 eigen vector는 1열에, -5에 해당되는 eigen vector는 2열에 존재한다.

sum((eigen(A)$vectors[,1])^2)

[1] 1

sum((eigen(A)$vectors[,2])^2)

[1] 1

# 이때 vecotor는 정규백터 (norm이 1)로 만들어졌다는 이야기다.

matrix.trace(A) == sum(eigen(A)$values) # eigen value의 첫 번째 성질이다. tr(A) = sum of eigenvalues (수학적 귀납법으로 간단히 증명가능하다.)

[1] TRUE

det(A)

[1] -35

prod(eigen(A)$values) # eigen value의 두 번째 성질이다. det(A) = product of eigenvalues (이 또한, 수학적 귀납법으로 증명 가능하다.)

[1] -35

det(A) == prod(eigen(A)$values) # 왜 False가 나올까?? 위에서는 같았는데 말이다.

[1] FALSE

det(A) -prod(eigen(A)$values) # 거의 0가 나온다. 컴퓨터의 구동방식에 의해 요렇게 나온다..

[1] -7.105427e-15

 

Gram-schmidt process of orthogonalization

install.packages("pracma")
library(pracma)
x1 <- c(2,0,0)
x2 <- c(2,2,0)
x3 <- c(2,0,2)
X <- cbind(x1,x2,x3)
gramSchmidt(X) # QR decomposition 즉, X = QR이며 Q는 직교행렬, Q는 상삼각행렬이다.

$Q
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1

$R
     [,1] [,2] [,3]
[1,]    2    2    2
[2,]    0    2    0
[3,]    0    0    2

더 깊은 QR decompostition은 아래 링크를 참조하길 바란다.

QR decomposition - Wikipedia

QR decomposition - Wikipedia

Quadratic forms and definite matrices

A <- matrix(c(3,5,1,5,13,0,1,0,1),3,3)
all(eigen(A)$values >0)

[1] TRUE

eigen(A)$values

[1] 15.08151247  1.88327962  0.03520791

뭐,, all of eigenvalues 가 양수라면, definite matrice라는 것인데.. 이렇게 보일 필요있나.. 근데 교수님 사랑해요.

 

Special matrices - Orthogonal matrix

P <- matrix(c(-1,0,0,0,-1,0,0,0,1),3,3)

t(P) == solve(P) # 교재에는 sum(t(P) - solve(P))로 되어있는데, 내 코드가 조금더 증명이 명확하다.

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

det(P)

[1] 1

all(abs(diag(P))<=1)

[1] TRUE

Special matrices - Idempotent matrix

library(matrixcalc)
A <- matrix(c(2,-1,1,-2,3,-2,-4,4,-3),3,3)
A%*%A == A # 교재에서는 sum(A%*%A -A)로 나와있지만, 이렇게 보이는게 더 명확하다고 생각한다.!

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

matrix.rank(A) == matrix.trace(A)

[1] TRUE

Decomposition of a matrix 

행렬의 대각화는 유명하니깐~~증명은 생략하겠슴다!

A <- matrix(c(3,5,1,5,13,0,1,0,1),3,3)
D <- diag(eigen(A)$values)
D # 대각행렬

         [,1]    [,2]       [,3]
[1,] 15.08151 0.00000 0.00000000
[2,]  0.00000 1.88328 0.00000000
[3,]  0.00000 0.00000 0.03520791

p1 <- eigen(A)$vectors[,1]
p2 <- eigen(A)$vectors[,2]
p3 <- eigen(A)$vectors[,3]
P <- cbind(p1,p2,p3)
P # A = PDP_t 가능하게 해주는 P

              p1         p2         p3
[1,] -0.38418583  0.6344832  0.6706954
[2,] -0.92285258 -0.2853734 -0.2586603
[3,] -0.02728299  0.7183266 -0.6951709

P1 <- P%*%D%*%t(P)
P11 <- D[1,1]*p1%*%t(p1) +D[2,2]*p2%*%t(p2) +D[3,3]*p3%*%t(p3)

abs(sum(A-P1))

[1] 1.573915e-14

abs(sum(A-P11))

[1] 1.307461e-14

abs(sum(P1-P11))

[1] 2.664535e-15 # 동일하다는 것을 알 수 있다.

P2 <- P%*%solve(D)%*%t(P) # A^(-1)도 마찬가지로 할 수 있다는 것을 보여준다. 단, 고유값의 역수가 활용된다!
P22 <- (1/D[1,1])*p1%*%t(p1) +(1/D[2,2])*p2%*%t(p2) +(1/D[3,3])*p3%*%t(p3)

abs(sum(solve(A)-P2))

[1] 1.44329e-14

abs(sum(solve(A)-P22))

[1] 1.798561e-14

abs(sum(P2-P22))

[1] 3.552714e-15 # 동일하다는 것을 알 수 있다.

Kronecker(direct) product of matrices

아직 필자가 제대로 활용하는 것을 못봐서,,ㅠㅠ 그냥 코드만 적어놓겠다.

A <- matrix(c(1,2,3,4),2,2)
B <- matrix(c(2,3,4,5),2,2)
C <- matrix(c(3,4,5,6),2,2)
D <- matrix(c(4,5,6,7),2,2)
kronecker(A,B)

     [,1] [,2] [,3] [,4]
[1,]    2    4    6   12
[2,]    3    5    9   15
[3,]    4    8    8   16
[4,]    6   10   12   20
kronecker(A,B)%*%kronecker(C,D)

     [,1] [,2] [,3] [,4]
[1,]  420  600  644  920
[2,]  555  795  851 1219
[3,]  616  880  952 1360
[4,]  814 1166 1258 1802

kronecker(A%*%C,B%*%D)

     [,1] [,2] [,3] [,4]
[1,]  420  600  644  920
[2,]  555  795  851 1219
[3,]  616  880  952 1360
[4,]  814 1166 1258 1802 

# 즉, kronecker(A,B)%*%kronecker(C,D)과 kronecker(A%*%C,B%*%D)이 같다.

a <- 1:3 ; b <- 4:6
all.equal(a%*%t(b),kronecker(a,t(b)),kronecker(t(b),a))

[1] TRUE