Part I
벡터 / 행렬 / 기본연산 / 선형결합
transpose / trace
단위행렬 / 대칭행렬
선형방정식
행렬식
선형독립 / rank
역행렬
Part II
eigenvalue / eigenvector
spectral decomposition
quadratic form / pd / nnd
projection / idempotency
partitioned matrix
Part I 의 주요 내용에 대한 간략한 소개
자습자료: [Garth Tarr]
참고서적: [Searle 도서관에 있음 ], [박흥선]
자습 필수: 강의를 하지 않은 내용이라도 강의를 한 것으로 간주
part II 에 대한 설명
개념의 의미를 설명함
수치적인 예를 제시함
특성, 변환, 분해법, 표준형 사이의 관계 :
행렬이 가진 특정한 특성의 동일성을 보기 위하여,
행렬을 변환하고, 그 과정을 행렬을 분해하는 형태로 표현하기도 하고,
그 표현 중 특정한 형태를 표준형으로 규정하게 되므로
행렬의 변환법, 분해법, 표준형 이라는 용어는 서로 밀접하게 관련되어 있다
행동치 row equivalence : [wiki]
가우스 소거법을 적용하여 행단위로 변환할 수 있는 행렬들을 행동치로 설정
행동치 관계인 행렬들은 해가 동일한 선형방정식을 만들게 됨
동치 equivalence : [wiki]
행동치와 열동치 개념을 동시 적용함
동치관계인 행렬들은 동일한 rank 를 가짐
닯음분해 similarity : [wiki]
정방행렬을 대상으로 하여, 행렬의 대각화를 위하여
행렬의 고유값과 고유벡터를 구하기 위한 과정
직교닮음 분해 orthogoanl similarity :
대칭행렬의 경우에 대한 닮음분해. 대칭행렬의 고유값과 고유벡터를 구하기 위한 과정
스펙트럴 분해법을 행렬에 적용하는 경우와 동일한 결과
닮음 similarity 와 합동성 congruence 이 동시에 만족되는 경우
특이값 분해 singular value decomposition : [wiki]
기타 많은 수의 행렬 분해법과 표준형이 사용되고 있음. [wiki]
선형방정식 관련 : LU / Cholesky / QR …
고유값 관련 : Jordan / Schur …
정방행렬을 대각행렬로 변환하는 방법
대각행렬 계산의 편리성 :
예를 들어, \(D = \begin{pmatrix} d_1 & 0 & 0 \\ 0 & d_2 & 0 \\ 0 & 0 & d_3 \end{pmatrix}\) 라 하자
만약 \(D^5\) 을 계산할 필요가 있다면, 대각행렬의 경우 행렬의 곱셈을 5번 하지 않고서도
\(D^5 = \begin{pmatrix} d_1^5 & 0 & 0 \\ 0 & d_2^5 & 0 \\ 0 & 0 & d_3^5 \end{pmatrix}\) 을 쉽게 얻을 수 있다.
이와 같은 대각행렬의 특성을 이용하여, 효과적으로 비대각행렬의 계산을 수행하는 방법을 생각해 보자
대각행렬의 특징 :
여기서, \(\mathbf{u}_1 =\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}\), \(\mathbf{u}_2 =\begin{pmatrix} 0 \\ 1 \\ 0 \end{pmatrix}\), \(\mathbf{u}_3 =\begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix}\) 라 하면,
\(k=1,2,3\) 에 대하여, 식 \(D \, \mathbf{u}_k = d_k \mathbf{u}_k\) 이 만족된다
(정방)행렬 \(A\)의 고유값 eigenvalue, 고유벡터 eigenvector,
식 \(A \mathbf{u} = d \, \mathbf{u}\) 을 만족하는 값 \(d\)와 벡터 \(\mathbf{u}\) 가 있다
고유값 : 위 식을 만족하는 값 \(d\) (실수 혹은 복소수)
고유벡터 : 영벡터가 아니면서,위 식을 만족하는 벡터 \(\mathbf{u}\)
예, \(A = \left( \begin{array}{rrr} 2 & 2 & 0 \\ 2 & 1 & 1 \\ -7 & 2 & -3 \end{array} \right)\) 의 경우
고유값 : \(d_1 = 1\), \(d_2 = 3\), \(d_3 = -4\)
고유벡터 : \(\mathbf{u}_1 =\left( \begin{array}{r} -2 \\ 1 \\ 4 \end{array} \right)\), \(\mathbf{u}_2 =\left( \begin{array}{r} -2 \\ -1 \\ 2 \end{array} \right)\), \(\mathbf{u}_3 = \left( \begin{array}{r} 1\\ -3 \\ 13 \end{array} \right)\)
\(A \, \mathbf{u}_k = d_k \mathbf{u}_k , \quad k=1,2,3\)
\(\mathbf{u}= -1 \, \mathbf{u}_1 + 2 \, \mathbf{u}_2\) 이면, \(A \, \mathbf{u}= (-1)\, (1) \, \mathbf{u}_1 + (2) \, (3) \, \mathbf{u}_2\)
\({A^2 \mathbf{u}_3 = A \left( A \mathbf{u}_3 \right) = A \left( (-4) \mathbf{u}_3 \right) = (-4) \, A \mathbf{u}_3 = (-4)^2 \,\mathbf{u}_3 }\) 이므로,
\(A^2\)의 고유값과 고유벡터 : \(d_k^2\), \(\mathbf{u}_k\), \(k=1,2,3\)
고유벡터 행렬 \(U\), 고유값 대각행렬 \(D\)
\(U\) : 고유벡터를 열벡터로 구성한 행렬,
\(D\) : 고유벡터의 순서에 따라 대응하는 고유값을 대각원소로 배치한 대각행렬
\(A \mathbf{u} = d \, \mathbf{u}\)인 관계를 행렬 형태로 표현하면, \(AU = UD\) 임
위의 예에서, \(U=\left( \begin{array}{rrr} -2& -2& 1\\ 1& -1& -3\\ 4& 2& 13\end{array} \right)\), \(D= \left( \begin{array}{rrr} 1 & 0& 0\\0 & 3& 0\\0 & 0& -4\end{array} \right)\)
행렬의 닮음 표준형 similar canonical form : \(\quad AU=UD\)
\(U\) 의 역행렬이 존재하면, \(A = UDU^{-1}\) 로 표현됨.
\(U\) 의 역행렬이 존재하지 않는 경우, \(A\) 를 deficient 행렬, defective 행렬이라 함.
다음의 경우에는 \(U\) 의 역행렬이 존재함이 증명되어 있음
고유값이 모두 다른 실수인 경우
\(A\) 가 대칭행렬인 경우
행렬 \(A\) 가 deficient 행렬이 아니고, \(A = UDU^{-1}\) 와 같이 표현되면,
\(A^2 = U D U^{-1} U D U^{-1}= U D^2 U^{-1}\)
이를 일반화 하면, \(A^c = U D^c U^{-1}\), 여기서 \(c\) 는 음수가 아닌 모든 정수.
즉, \(A^3 = U D^3 U^{-1}\) 이 됨.
이 때, \(D\)는 대각행렬이므로, \(D^3\) 은 고유값의 3제곱 값을 대각원소로 배치한 행렬.
\(A+3A^2\) 은 마찬가지로 \(U \left(D+3D^2\right) U^{-1}\) 으로 표현됨.
\(\exp(x)=1+x +x^2/2 +x^3 /3 + \cdots\) 이므로,
\(\exp(A)=U \left(I+D +D^2/2 + \cdots \right)U^{-1}\) 이고, 결국 \(\exp(A)=U \, \exp(D) \, U^{-1}\) 임.
여기서 \(\exp(D)\) 은 대각원소에 고유값 \(d_k\) 들에 대하여 \(\exp(d_k)\) 를 계산한 값을 배치한 대칭행렬.
만약 고유값이 0이 아니면, \(A^{-1} = U D^{-1} U^{-1}\) 임.
만약, 모든 고유값이 양수라면, \(c\) 는 모든 실수로 확장할 수 있음.
즉, \(A^{1/2} = U D^{1/2} U^{-1}\) 라고 할 수 있음 ( 말하자면, \(\sqrt{A}\) 를 구할 수 있음).
A= t(matrix(c(2,2,0, 2,1,1,-7, 2, -3),3,3))
u1 = matrix(c(-2, 1, 4), ncol=1)
u2 = matrix(c(-2, -1, 2), ncol=1)
u3 = matrix(c(1, -3, 13), ncol=1)
u = -1*u1 + 2*u2
A## [,1] [,2] [,3]
## [1,] 2 2 0
## [2,] 2 1 1
## [3,] -7 2 -3cbind(u1, u2, u3, u)## [,1] [,2] [,3] [,4]
## [1,] -2 -2 1 -2
## [2,] 1 -1 -3 -3
## [3,] 4 2 13 0cbind(A %*% u1, A %*% u2, A %*% u3, A %*% u, (-1)*(1)*u1+(2)*(3)*u2 )## [,1] [,2] [,3] [,4] [,5]
## [1,] -2 -6 -4 -10 -10
## [2,] 1 -3 12 -7 -7
## [3,] 4 6 -52 8 8( D=diag(c(1,3,-4)) )## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 3 0
## [3,] 0 0 -4( U=cbind(u1,u2,u3) )## [,1] [,2] [,3]
## [1,] -2 -2 1
## [2,] 1 -1 -3
## [3,] 4 2 13A %*% A %*% A # A^3## [,1] [,2] [,3]
## [1,] 14 26 0
## [2,] 26 1 13
## [3,] -91 26 -51U %*% diag( c(1,3,-4)^3 ) %*% solve(U) ## [,1] [,2] [,3]
## [1,] 14 26 4.440892e-16
## [2,] 26 1 1.300000e+01
## [3,] -91 26 -5.100000e+01solve(A) # A^{-1}## [,1] [,2] [,3]
## [1,] 0.41666667 -0.5 -0.1666667
## [2,] 0.08333333 0.5 0.1666667
## [3,] -0.91666667 1.5 0.1666667U %*% diag( c(1,3,-4)^{-1} ) %*% solve(U)## [,1] [,2] [,3]
## [1,] 0.41666667 -0.5 -0.1666667
## [2,] 0.08333333 0.5 0.1666667
## [3,] -0.91666667 1.5 0.1666667exp(A) # What ? Why ? Don't confuse of vector operation in R.## [,1] [,2] [,3]
## [1,] 7.389056099 7.389056 1.00000000
## [2,] 7.389056099 2.718282 2.71828183
## [3,] 0.000911882 7.389056 0.04978707U %*% diag( exp(c(1,3,-4)) ) %*% solve(U)## [,1] [,2] [,3]
## [1,] 14.892038 15.04050 2.326753
## [2,] 6.896868 9.69854 1.703370
## [3,] -15.413716 -12.88053 -1.768444eigen(A)## eigen() decomposition
## $values
## [1] -4 3 1
##
## $vectors
## [,1] [,2] [,3]
## [1,] 0.07474351 -0.6666667 -0.4364358
## [2,] -0.22423053 -0.3333333 0.2182179
## [3,] 0.97166562 0.6666667 0.8728716unit.vector <- function(u) { u/sqrt(sum(u*u))}
cbind(unit.vector(u3), unit.vector(u2), unit.vector(u1))## [,1] [,2] [,3]
## [1,] 0.07474351 -0.6666667 -0.4364358
## [2,] -0.22423053 -0.3333333 0.2182179
## [3,] 0.97166562 0.6666667 0.87287162차 형식 : \(\quad Q(A) = \mathbf{x} ' A \mathbf{x}\)
\(\mathbf{x}=(x_1, \ldots, x_p)'\) 라고 하고, \(A_{p \times p}\) 이면,
2차 형식 \(Q(A)\) 를 통하여, \(\mathbf{x}=(x_1, \ldots, x_p) '\) 에 대한 2차 함수 표현이 가능하다.
2차 형식에 사용되는 행렬 \(A\) 는 대칭행렬이라 본다.
양수 비음수 개념의 행렬로의 확장
비음수 : ‘실수 \(a\)가 비음수’ \(~\Leftrightarrow~\) ‘모든 \(x\) 에 대하여, \(ax^2 \ge 0\)’
양수 : ‘실수 \(a\)가 양수’ \(~\Leftrightarrow~\) ‘0 아닌 모든 \(x\) 에 대하여, \(ax^2 >0\)’
행렬 \(A\)는 비음정치 행렬 \(~\Leftrightarrow~\) 임의의 \(\mathbf{x}\) 에 대하여, \(Q(A) \ge 0\)
행렬 \(A\)는 양정치 행렬 \(~\Leftrightarrow~\) 영벡터가 아닌 임의의 \(\mathbf{x}\) 에 대하여, \(Q(A) > 0\)
양정치 행렬의 역행렬은 양정치 행렬이다. (why ?)
어떤 행렬 \(A\)에 대하여, \(A'A\) 형태로 주어지는 행렬은 대칭행렬이다.
어떤 행렬 \(A\)에 대하여, \(A'A\) 형태로 주어지는 행렬은 비음정치 행렬이다.
행렬 \(A\) 의 열벡터가 모두 선형독립이면, \(A'A\)은 양정치 행렬이다.(why ?)
( A <- matrix(c(2,-1,-1,2),2,2))## [,1] [,2]
## [1,] 2 -1
## [2,] -1 2egnA<- eigen(A)
egnA## eigen() decomposition
## $values
## [1] 3 1
##
## $vectors
## [,1] [,2]
## [1,] -0.7071068 -0.7071068
## [2,] 0.7071068 -0.7071068library(plotly, warn.conflicts = FALSE, quietly=TRUE)
Q <- function(xx){ t(xx) %*% A %*% xx }
x1sn <- seq(-2,2, l=41)
x2sn <- seq(-2,2, l=41)
ss <- expand.grid(x1=x1sn, x2=x2sn)
Q1 <- apply( ss, 1, Q)
# Q1[abs(ss$x1-ss$x2)>2] <- NA
Q <- t(matrix(Q1, ncol=41 ))
# Q[Q>4]<-NA
pp <- plot_ly(z=~Q,x=x1sn, y=x2sn) %>%
add_surface(showscale=FALSE,
contours = list(z = list(show=TRUE,
usecolormap=TRUE,
highlightcolor="#ff0000",
project=list(z=TRUE)
)
)
)
pp %>% layout(title="Quadratic form with a positive definite matrix",
scene=list(xaxis=list(title="x1",range=c(-3,3)), yaxis=list(title="x2",range=c(-3,3)))
)직교닮음 표준형 : \(\quad A=UD U'\) (여기서 \(A\)는 대칭행렬)
행렬의 닮음 표준형 \(AU=UD\) 표현에서, 실수 대칭행렬의 경우 \(U^{-1}\) 이 존재한다.
이 경우, \(A=UD U^{-1}\) 라고 표현이 가능하고,
\(A\) 가 대칭행렬이므로, \(UDU^{-1} = A = A'= (U^{-1})' D U'\) 인 관계가 성립한다.
\(U^{-1} = U'\) 인 경우, 이 조건을 만족한다. 즉, \(A=UD U'\) 이다.
\(U'\) 가 \(U\)의 역행렬이 되는 경우이므로, \(U'U = UU' =I\) (단위행렬)이라는 의미이고,
단위행렬의 비대각 원소는 0 이고 대각원소는 1이므로,
\(j\neq k\) 일 때, \(\mathbf{u}_j ' \mathbf{u}_k =0\) 이고, \(\mathbf{u}_k ' \mathbf{u}_k = 1\) 이라는 의미이다.
이는 \(U\) 의 모든 열벡터의 길이(norm)는 1이고, 서로 간에 직교한다는 의미이다.
즉, \(U\)가 직교행렬이라는 의미다.
이 경우, 즉, \(U^{-1} = U'\) 일 때, \(A=UD U'\) 를 직교닮음 표준형 표현이라고 한다.
스펙트럴 분해법 (\(A\) 는 대칭행렬) :
고유값과 고유벡터를 이용하여 \(A =\sum_k d_k \mathbf{u}_k \mathbf{u}_k '\) 와 같이 표현되고,
앞서와 마찬가지로, \(A^c =\sum_k d_k^c \mathbf{u}_k \mathbf{u}_k '\) 이다. 여기서 \(c\) 는 비음인 정수
만약 모든 고유값 \(d_k\) 가 0이 아니라면, \(A^{-1} =\sum_k d_k^{-1} \mathbf{u}_k \mathbf{u}_k '\) 가 성립함.
만약 모든 고유값 \(d_k\) 가 양수라면, \(A^{1/2} =\sum_k d_k^{1/2} \mathbf{u}_k \mathbf{u}_k '\) 가 성립함.
( A <- matrix(c(9,-5,2,-5,9,2,2,2,2),ncol=3) )## [,1] [,2] [,3]
## [1,] 9 -5 2
## [2,] -5 9 2
## [3,] 2 2 2( d<- eigen(A)$values ) # ill conditioned matrix## [1] 1.400000e+01 6.000000e+00 1.065814e-14( U <- eigen(A)$vectors )## [,1] [,2] [,3]
## [1,] 0.7071068 0.5773503 -0.4082483
## [2,] -0.7071068 0.5773503 -0.4082483
## [3,] 0.0000000 0.5773503 0.8164966round( t(U) %*% U , 12)## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] 0 0 1u1<- U[,1]; u2<- U[,2]; u3 <- U[,3]
u1 %*% t(u1)## [,1] [,2] [,3]
## [1,] 0.5 -0.5 0
## [2,] -0.5 0.5 0
## [3,] 0.0 0.0 0# check spectral decompsition
d[1]*u1 %*% t(u1) + d[2]*u2 %*% t(u2) + d[3]*u3 %*% t(u3)## [,1] [,2] [,3]
## [1,] 9 -5 2
## [2,] -5 9 2
## [3,] 2 2 2양정치 행렬 \(~\Leftrightarrow~\) 모든 고유값들이 양수
양정치 행렬의 역행열도 양정치 행렬이다 (역행렬의 고유값은 역수)
비음정치 행렬 \(~\Leftrightarrow~\) 모든 고유값들이 비음이다.
행렬식의 값은, 고유값들의 곱과 같다 (why ?)
행렬의 rank 는 0 아닌 고유값의 개수와 같다
행렬의 대각원소의 합 trace 는, 고유값의 합과 같다
행렬 \(X_{n \times k}\) 와 \(n\) 차원 벡터 \(\mathbf{y}\) 가 주어져 있을 때, 선형방정식 \(X \mathbf{b} = \mathbf{y}\) 의 해 \(\mathbf{b}\) :
\(\mathbf{y} \in Col(X)\) 이면, 일반해는 \(\mathbf{b}= X^{-} \mathbf{y} +(I-X^{-}X) \mathbf{z}, ~~ \mathbf{z} \in I\!\!R^k\) 이다
\(\mathbf{y} \not\in Col(X)\) 이면, 해가 없고 근사해를 구할 수 있다.
선형방정식 \(X \mathbf{b} = \mathbf{y}\) 에 대한 최소자승법에 의한 근사해 \(\mathbf{b}\) 는,
정규방정식 \(X' X \mathbf{b} = X' \mathbf{y}\) 으로부터 구한다
\(X' \mathbf{y} \in Col(X')\) 이므로, 정규방정식의 해는 항상 존재한다
정규방정식 \(X' X \mathbf{b} = X' \mathbf{y}\) 에서,
행렬 \(X\)의 열벡터들이 모두 선형독립이면 ( \(X\) 가 full column rank 이면),
\((X'X)^{-1}\) 이 (유일하게) 존재하고, \((X'X)^{-1}=(X'X)^{-}\) 이다.
정규방정식의 해는, \(\mathbf{b} = (X'X)^{-1} X' \mathbf{y}\) 로 유일하다
행렬 \(X\)의 열벡터들이 선형독립이 아니면,
일반해는 \(\mathbf{b}= (X'X)^{-} X'\mathbf{y} +(I-(X'X)^{-}X'X) \mathbf{z}, ~~ \mathbf{z} \in I\!\!R^k\) 이다
\((X'X)^{-}\)은 존재하고(유일하지 않음), 그 선택에 따라 일반해 \(\mathbf{b}\) 도 다양한 값이 될 수 있다
그러나 모든 일반해에 대하여 \(X \mathbf{b}\)는 동일한 값이다
인공신경망의 특성을 이해하는 데 있어서 이 점이 매우 중요함
참고 1: 일반화 역행렬 \((X'X)^{-}\) 의 성질
\(G=(X'X)^{-}~\) 일 때, 전치행렬 \(~G'~\) 도 \(~X'X~\)의 일반화 역행렬이다
\((X'X)(X'X)^{-}(X'X) = (X'X)\) 이다
\(X(X'X)^{-}(X'X) = X\) 이다. 즉, \((X'X)^{-}X' = X^{-}\) 이다. (참고 2 이용) (why ?)
\(X(X'X)^{-}X'\) 는 불변이다. 즉, \((X'X)^{-}\) 의 선택에 영향을 받지 않고 동일하다 (why?)
\(X(X'X)^{-}X'\) 는 대칭행렬이다 (\(~(X'X)^{-}~\)는 대칭이 아닐 수 있지만 )
참고 2: \(PX' X = Q X'X\) \(~\Longrightarrow~\) \(PX'= Q X'\) 에 대한 증명
\((PX'X -QX'X) (P' -Q') = (PX' -QX')(PX' -QX')'\)
\(A'A=\mathbf{0} ~~\Rightarrow~~ A=\mathbf{0}\)
사영행렬: 행렬 \(P\) 가 다음 두 조건을 만족하면 사영행렬이라 한다.
대칭성 symmetricity : \(P=P'\)
멱등성 idempotency : \(P=P^2\)
멱등성의 의미: 두 번 사영을 한 이미지는 한 번 사영한 것과 동일하다
사영의 직교성 :
\(0 = P-P^2 = P(I-P) = P' (I-P)\) 이고, \(\mathbf{y}' P' (I-P)\mathbf{y}=0\) 이므로
이미지 벡터 \(P\mathbf{y}\) 와 잔차 벡터 \((I-P)\mathbf{y}\) 는 직교한다
사영행렬의 고유값은, 0 혹은 1 이다.
\(P\mathbf{u} =d \mathbf{u}\) 이고, \(P^2\mathbf{u}= d \, P \mathbf{u} = d^2 \mathbf{u}\) 이고, \(P=P^2\) 이므로
(영벡터가 아닌) 고유벡터 \(\mathbf{u}\)에 대하여, \(d \mathbf{u}= d^2 \mathbf{u}\) 이므로, \(d^2=d\) 이고, \(d=0,1\) 임.
사영행렬 \(P\)에 대하여, \(trace(P)=rank(P)\) 이다. (why?)
\(P=X(X'X)^{-} X'\) 일 때,
행렬 \(P_{n \times n}\) 는 대칭성과 멱등성을 만족하여 사영행렬이 된다.
\(P\) 의 이미지 공간은 행렬 \(X\)의 열공간이다. 즉, \(P\mathbf{y} \in Col(X)\)
\(rank(P)=rank(X)\) 이다.
\(P_{n \times n}\) 일 때, \(rank(I-P)=n-rank(X)\) 이다.
행렬 \(X\) 가 full column rank 라면,
\(rank(X_{n \times k})=rank(X'X)=k\) 이고(why?),
\(P=X(X'X)^{-} X' = X(X'X)^{-1} X'\) 임.(why?)
모자행렬 hat matrix
보통 (열)벡터 \(~\mathbf{y}~\) 를 사영하여 얻은 결과를 \(~\hat{\mathbf{y}}~\) 으로 표현한다. 즉, \(~\hat{\mathbf{y}}= P \mathbf{y}~\) 이다
사영행렬 \(P_{n \times n}\) 은 \(\mathbf{y}\) 에 모자(hat) 를 씌워 \(\hat{\mathbf{y}}\) 을 만들게 되므로,
사영행렬을 다른 이름으로 hat matrix 라고 부르고, \(P\) 대신 \(H\) 라는 기호로 나타내기도 함
hat matrix (사영행렬)에는 \(n\) 개의 대각원소들이 있고 이를 hat values 라고 함.
행렬 \(X_{n \times k}\) 에는, \(n\) 개의 \(k\) 차원의 행벡터가 들어 있음.
\(n\) 개의 \(k\) 차원 행벡터 중, 다른 벡터들과 멀리 떨어진 값을 갖는 경우,
대응하는 hat value 가 다른 경우들 보다 큰 값을 갖게 됨.
사영행렬 계산 예 :
trees 데이터는 31 개의 체리나무에 대한 측정 자료이다
아래에서는, 간단한 예로, trees 자료의 앞 부분 8개의 관측값으로 구성된 행렬 \(X\) 로부터
사영행렬 \(P\) 를 얻고, 대칭성과 멱등성이 만족되는 지를 시험해 본다.
X <- as.matrix(trees[1:8,]); X## Girth Height Volume
## 1 8.3 70 10.3
## 2 8.6 65 10.3
## 3 8.8 63 10.2
## 4 10.5 72 16.4
## 5 10.7 81 18.8
## 6 10.8 83 19.7
## 7 11.0 66 15.6
## 8 11.0 75 18.2P = X %*% solve(t(X) %*% X ) %*% t(X) # projection matrix
round(P,5)## 1 2 3 4 5 6 7 8
## 1 0.52038 0.33821 0.26144 -0.01295 0.09670 0.09447 -0.20415 -0.08199
## 2 0.33821 0.30104 0.29559 0.05227 -0.00573 -0.03988 0.06398 -0.01436
## 3 0.26144 0.29559 0.32597 0.08237 -0.06468 -0.11726 0.19777 0.01252
## 4 -0.01295 0.05227 0.08237 0.17247 0.12735 0.12473 0.25112 0.19546
## 5 0.09670 -0.00573 -0.06468 0.12735 0.33424 0.39250 -0.06254 0.18647
## 6 0.09447 -0.03988 -0.11726 0.12473 0.39250 0.46767 -0.12384 0.19935
## 7 -0.20415 0.06398 0.19777 0.25112 -0.06254 -0.12384 0.63368 0.25484
## 8 -0.08199 -0.01436 0.01252 0.19546 0.18647 0.19935 0.25484 0.24455all.equal(P,t(P)) # checking symmetricity## [1] TRUEall.equal(P %*% P, P) # checking idempotency ## [1] TRUE수치 오차에 의한 대칭성의 상실
이론적으로는 당연히 대칭행렬이 되어야 하는 경우이지만,
막상 계산 과정에서, 수치 오차가 발생하여, 행렬의 대칭성이 깨질 수도 있다.
대칭성 상실의 효과
행렬의 대칭성이 깨지게 되면, 예를 들어, 이후 고유값을 구하는 과정에서
실수인 고유값이 나와야 함에도, 고유값이 복소수가 되어
후속 계산이 불가능해지는 경우와 같은 문제가 발생하기도 한다.
대칭성 복구
이런 경우 원래의 행렬 \(A\) 를 사용하는 대신
\((A+A')/2\) 를 사용하여 대칭성을 복구할 수 있다.
힐버트 행렬에서의 대칭성 상실
아래의 예에서는, 크기가 4인 힐버트 행렬, 크기가 11인 힐버트 행렬에 대하여
역행렬을 구하여 그 역행렬들의 대칭성을 시험해 본다.
두 경우 모두 대칭행렬이 되어야 하지만, 크기가 11인 힐버트 행렬의 경우
수치오차로 인하여 대칭성이 깨진 것으로 나타난다.
이 경우 대칭화 변환을 통하여, 대칭성을 복구할 수 있다
# 'symmetric matrix' may not be symmetric in computing
hilbert <- function(n) { i <- 1:n; 1 / outer(i - 1, i, "+") }
H4 <- hilbert(4)
iH4 <- solve(H4)
M <- iH4[1:4,1:4] ; M## [,1] [,2] [,3] [,4]
## [1,] 16 -120 240 -140
## [2,] -120 1200 -2700 1680
## [3,] 240 -2700 6480 -4200
## [4,] -140 1680 -4200 2800all.equal(M,t(M))## [1] TRUEH11 <- hilbert(11)
iH11 <- solve(H11)
M <- iH11[1:4,1:4] ; M## [,1] [,2] [,3] [,4]
## [1,] 120.9175 -7251.294 141342.6 -1318762
## [2,] -7251.2918 579881.054 -12717294.3 126576692
## [3,] 141342.4750 -12717291.722 297519276.8 -3084834612
## [4,] -1318760.5400 126576660.825 -3084834543.5 32900801617all.equal(M,t(M)) ## [1] "Mean relative difference: 3.213717e-08"S <- (M+t(M))/2 # giving symmetricity
all.equal(S,t(S))## [1] TRUE\(\qquad \qquad \qquad Q = \begin{pmatrix} A & B \\ C & D \end{pmatrix}\)
\(\qquad \qquad \qquad Q^{-1}=\begin{pmatrix} A^{-1}+A^{-1} B M^{-1} C A^{-1} & - A^{-1} B M^{-1} \\ - M^{-1} C A^{-1} & M^{-1} \end{pmatrix}\)
\(\qquad \qquad \qquad Q^{-1}=\begin{pmatrix} W^{-1} & - W^{-1} B D^{-1} \\ - D^{-1} C W^{-1} & D^{-1} + D^{-1} C W^{-1} B D^{-1} \end{pmatrix}\)
\(\qquad \qquad \qquad (A-B D^{-1} C)^{-1} = A^{-1}+A^{-1} B (D-C A^{-1} B )^{-1} C A^{-1}\)
\(\qquad \qquad \qquad \begin{pmatrix} A & B & |& I & \mathbf{0} \\ C & D & |& \mathbf{0} & I \end{pmatrix}\)
\(\qquad \qquad \qquad ~~\approx ~~ \begin{pmatrix} I & A^{-1} B & | & A^{-1} & \mathbf{0} \\ \mathbf{0} & D-C A^{-1} B & | & - C A^{-1} & I \end{pmatrix}\)
\(\qquad \qquad \qquad ~~\approx ~~ \begin{pmatrix} I & \mathbf{0}& | & A^{-1}+A^{-1} B M^{-1} C A^{-1} & - A^{-1} B M^{-1} \\ \mathbf{0} & I & | & - M^{-1} C A^{-1} & M^{-1} \end{pmatrix}\)