Learn Again! 러너게인

이번 강의에서는 대칭 행렬(Symmetric Matrix)에 대해 이야기 하도록 하겠다. 지난 강의 에서 간략히 배우긴 했으나, 이번 강의에선 고유값과 고유벡터의 관점에서 대칭 행렬의 특성에 관한 내용을 다룰 것이다. 대칭 행렬은 굉장히 좋은 특성을 가지고 있기 때문에 특이값 분해(Singular Value Decomposition), 주성분 분석(Principal Component Analysis)등 여러 분야에 응용될 수 있다. 

1. 대칭 행렬의 개념

- Basic concept of symmetric matrix

이미 지난 강의 Lecture 5-(1)에서 공부했지만, 대칭행렬에 대해 간략히 정리해 보도록 하자. 대칭 행렬은 다음의 조건을 만족하는 행렬이다. 

대칭 행렬(symmetric matrix)는 어떤 행렬 A가 있다고 했을 때, 자신의 전치(transpose)행렬이 원래의 자기 자신과 같은 행렬이다. 전치라는 것은 행렬의 임의의 원소의 row와 column의 인덱스가 서로 바뀌어도, 즉 a_ij와 a_ji가 서로 같아야 한다. 이것은 식 (1)의 1번 조건을 만족하는 것이며, 이 조건을 만족하기 위해선 반드시 정방 행렬(square matrix)의 형태이어야만 한다. 

식 (1)의 3x3 행렬을 보면 row와 column의 인덱스가 같은 대각 원소(diagonal element)=[8, 3, 6]가 있다. 이 대각 원소를 기준으로 우측 상단의 상삼각행렬의 원소들과 좌측 하단의 하삼각행렬 원소들의 값이 같은 것을 볼 수 있다. 이처럼 전치 이후에도 그 결과가 완전히 같은 행렬을 대칭행렬이라 한다. 

2. 대칭행렬과 스펙트럼 정리(Spectral theorem)

- Eigenvalues and Eigenvectors of symmetric matrix

대칭행렬은 보통의 행렬과는 다른 특별한 형태의 행렬이다. 지난 강의 Lecture 24에서 우리는 마코브 행렬(Markov matrix)이 보통의 다른 행렬과는 달리 특별한 형태임을 배웠다. 모든 원소가 0보다 크고, 각 column의 합이 1이 되는 등의 특징을 보였는데, 이와 같이 특별한 특성을 보이는 행렬은 고유값과 고유벡터 역시 어떤 특징을 보이기 마련이다. 마코브 행렬은 적어도 하나의 고유값은 1을 가져야 한다는 특징을 보였다. 이와 마찬가지로 대칭행렬 역시 고유값과 고유벡터가 어떤 특징을 보인다. 그렇다면 대칭행렬은 어떤 특징을 보일까? 아래를 보도록 하자. 

대칭행렬은 (2)에 정리된 것과 같은 특성을 보인다. 첫 번째 특성은 대칭행렬의 고유값(eigenvalue)은 전부 실수(Real number)라는 것이다. Lecture 21-(2)에서 회전 행렬(rotation matrix)의 고유값을 구했을 때 허수(complex number)가 나왔었다. 하지만 대칭행렬에서는 절대 이와 같은 허수 고유값이 나오지 않는다. 그 이유는 잠시 후에 설명하도록 하겠다. 

두 번째 특성은 대칭행렬의 고유벡터(eigenvector)는 직각(perpendicular)을 이룬다는 것이다. 대칭행렬은 실수인 고유값들을 가지며, 각각의 고유값에 대응되는 고유벡터들은 1차원의 고유공간(eigenspace)을 형성한다. 그리고 각 고유 벡터, 즉 각각의 고유공간들은 서로 직교(orthogonal)하다는 것이 두 번째 특성이다. 

- Diagonalization of symmetric matrix

그렇다면 위의 특성으로부터 알 수 있는 것은 무엇일까? 일단 모든 고유벡터들이 서로 수직하다는 특성 2에 주목해보자. 이 특성은 어떤 nxn크기의 대칭행렬(symmetric matrix)은 n개의 서로 수직인 고유벡터들을 가지고 있다는 의미이다. 이는 행렬을 대각화하여 분해(factorization)할 때 특별한 형태를 띄게 된다. 그 특별한 형태를 알아보도록 하자. 

먼저 행렬의 대각화(diagonalization) (Lecture 22)에서 배웠듯이 어떤 행렬은 고유벡터행렬과 고유값 대각 행렬의 조합으로 분해될 수 있다. 

식 (3)은 임의의 행렬 A를 대각화하여 분해한 모습이며, S는 고유벡터행렬을, 대문자 람다는 고유값행렬을 각각 나타낸다. 그렇다면 대칭행렬은 분해할 때 무엇이 특별하다는 것일까? 특별한 형태가 나타나는 부분은 고유벡터행렬이다. 앞서 언급했듯이 대칭행렬의 고유벡터들은 서로 수직이며, 수직인 벡터들이 column 벡터 형태로써 고유벡터행렬 S를 형성한다. 이때 이 고유벡터들은 정규화(normalization)과정을 거쳐 그 크기를 1로 만들 수 있다. 고유벡터는 크기보단 방향이 의미가 있기 때문에 단위 벡터(unit vector)로 만들어도 문제가 없으며 실제로 MATLAB의 내장 함수 eig()를 이용하여 고유값을 구해도 정규화된 고유벡터가 나온다. 

어쨋든 대칭행렬의 고유벡터들은 서로 수직(perpendicular)하며 방향 성분만을 나타내는 단위 벡터(unit vector)로 만들 수 있다. 이들은 정규직교벡터(orthonormal vector)이며, 이 정규직교벡터인 고유벡터들을 모아 행렬을 만들었기 때문에 결과적으로 대칭행렬의 고유벡터행렬은 정규직교행렬(orthonormal matrix)이 된다. 이제 대칭행렬의 분해는 아래와 같은 식으로 표현할 수 있다. 

대칭행렬의 경우 식 (4)와 같이 원래의 고유벡터행렬 S가 정규직교행렬을 나타내는 Q로 대체되었다. 정규직교행렬은 Lecture 17-(1)에서 배운바 있다. 정규직교행렬은 역행렬이 전치(transpose)와 같음을 이미 배웠다. 따라서 (4.1)의 우변의 Q의 역행렬은 (4.2)와 같이 Q의 전치 행렬로 간단히 표현할 수 있다. 결과적으로 대칭행렬은 식 (4.2)와 같이 정규직교행렬 Q와 고유값 행렬 람다, 그리고 Q의 전치의 곱으로 분해하여 표현할 수 있으며, 이는 선형대수에서 굉장히 중요한 식이므로 잘 알아두도록 하자. 또한 반대로 해석하면 어떤 행렬을 (4.2)와 같이 분해하여 표현할 수 있다면, 그 행렬은 대칭행렬이라고 할 수 있다. 

식 (4.2)가 대칭행렬인지를 증명할 수 있는 방법이 한 가지 있다. 바로 (4.2)를 전치(transpose)시켜보는 것이다. 

식 (4.2)를 전치시키면 (4.3)을 거쳐 (4.4)와 같이 정리할 수 있다. 우선 각 행렬의 순서가 뒤집혀서 곱해지는데, 첫 번째 Q는 전치에 전치를 거쳐서 원래의 Q로, 고유값행렬은 대각행렬이기 때문에 원래 행렬 그대로, 마지막 Q는 전치된 행렬로 각각 정리할 수 있다. 결과 식인 (4.5)를 보면 전치를 하기 이전의 식인 (4.2)와 같은 것을 볼 수 있다. 결국 대칭행렬의 분해식을 기준으로 전치를 시켜도 분해식의 결과가 같기 때문에 대칭행렬의 분해를 식 (4.2)와 같이 정리할 수 있다. 

아래 그림은 식 (4.2)를 검증하기 위한 MATLAB 코드와 그 결과이다.  

Fig. 1 대칭행렬의 분해. [Left] 실행 결과. [Right] 소스코드

Fig. 1의 프로그램 결과를 살펴보면 먼저 대칭 행렬 A에 대한 고유값과 고유벡터를 구한 뒤, 이를 식 (4.2)와 같이 다시 곱하였다. 결과값이 정확히 원래의 대칭행렬과 일치하는 것을 볼 수 있다. 

- Special case

앞서 우리는 대칭행렬이 전부 실수(real number)인 고유값을 가지고, 고유벡터들은 모두 수직(perpendicular)하다는 것을 배웠으며 이를 식 (2)에 정리하였다. 그런데 식 (2)의 특성 2번의 "대칭행렬의 고유벡터들은 서로 수직이다"는 어떤 특수한 형태의 대칭행렬인 경우엔 약간 다르게 표현되어야 한다. 여기서 특수한 형태의 대칭행렬은 바로 단위 행렬(Identity matrix)을 의미하며, 이 경우 특성 2는 다음과 같이 표현될 수 있다. 

"고유벡터들은 서로 수직한 벡터들로 선택될 수 있다"

왜 이와 같이 표현하는지 당장은 이해가 가지 않을 것이다. 다음의 2x2 단위행렬의 고유값과 고유벡터를 구해서 그 이유를 알아보자. 

식 (5)는 대칭행렬이 단위행렬일 때의 고유값과 고유벡터를 나타낸 것이다. 행렬이 단위행렬이기 때문에 고유값은 식 (5.1)과 같이 A가 없는 것과 마찬가지로 볼 수 있다. 따라서 식을 만족시키는 고유값은 오직 1이기 때문에 두 개의 고유값 모두 1이 된다. 여기서 중요한 포인트는 고유값이 반복(repeated eigenvalue)된다는 것이다. 이렇게 고유값이 반복되는 경우엔 어떠한 벡터라도 고유벡터가 될 수 있다. 식 (5.2)를 보면 고유벡터를 구하는 과정에서 행렬의 모든 원소가 0이 되고, 이때의 null space(고유벡터공간을 의미)는 어떠한 벡터이든 가능하게 된다. 즉 2x2행렬 기준으로 전체 평면이 고유벡터가 되는 것이다. 이 수많은 벡터들중 어떠한 벡터들을 고유벡터로 선택하든 식은 성립한다. 그러나 어떠한 벡터든지 고유벡터가 될 수 있다면, 기왕이면 서로 수직한 벡터(위의 경우엔 x1=[1 0]T, x2=[0 1]T )를 선택하는 것이 여러 모로 더 유리하다고 볼 수 있다. 그렇기 때문에 식 (2)의 특성 2가 "고유벡터들은 서로 수직한 벡터들로 선택될 수 있다" 로 정의할 수 있는 것이다. 

- Spectral Theorem

지금까지 우리는 대칭행렬을 대각화하여 분해하는 과정을 살펴보았다. 이쯤에서 스펙트럼 정리(spectral theorem)에 대해 알아보자. 사실 스펙트럼 정리라는 것은 별다를 것이 없다. 우리가 지금까지 대각행렬을 대각화하여 분해하는 과정이 전부 스펙트럼 정리에 기반하여 이루어진 것이다. 스펙트럼 정리는 아래와 같다. 

스펙트럼 정리(Spectral Theorem) :  nxn 크기의 에르미트 행렬(Hermitian matrix 또는 self-adjoint matrix)은 아래 식과 같이 실수로 이루어진 고유값행렬와 유니터리 행렬(Unitary matrix) U로 대각화 과정을 통해 분해할 수 있다.

일단 식을 보아하니 (4.2)와 유사하게 생기긴 했는데, 어려운 용어가 나오고 지수부에 이상한 기호도 보여서 헷갈리는 분들도 있을 것이다. 그러나 어렵게 생각할 필요 없다. 차근차근 알아보자. 

에르미트 행렬(Hermitian matrix)은 행렬이 가질 수 있는 값들을 실수와 함께 복소수까지 고려하여 확장하여 설명한 개념이다. 즉 복소수를 원소로 가지고 있는 정방 행렬이고, 자기 자신과 그의 켤레전치(conjugate transpose)행렬이 같은 행렬을 의미하며 다음의 식을 만족시킨다. 

에르미트 행렬은 일단 (6.2)와 같이 대각 원소들은 반드시 실수(Real number)여야 한다. 또한 (6.2)의 왼쪽의 행렬 A에서 먼저 전치(transpose)를 한 뒤, 켤레 복소를 해주면 결국 원래의 행렬과 같아진다. 이러한 조건들을 만족시키는 행렬을 에르미트 행렬이라 한다. 지금까지 우리는 실수로만 이루어진 행렬만 다루어 왔지만, 사실 행렬의 원소로 복소수(complex number)가 올 수 있기 때문에 에르미트 행렬에 대해 알아두면 좋다. 

유니터리 행렬(unitary matrix) 역시 에르미트 행렬과 정의가 거의 유사하다. 복소수를 원소로 가질 수 있으며 원래의 행렬과 켤레전치(conjugate transpose)가 같은 행렬이다. 

식 (7)과 같이 자기 자신과 켤레전치행렬을 곱하면 단위 행렬(Identity matrix)이 되므로 결국 켤레전치행렬은 원래 행렬의 역행렬과 같다. 이 외에도 정규화된 행렬(normalized matrix)이고, 대각화가 가능(diagonalizable)하고, 고유 공간이 서로 직교(orthogonal)하는 등 여러 특성들이 존재한다. 자세한 사항은 위키를 참고 하자. (※ 켤레전치행렬(conjugate transpose matrix)의 표현은 선형대수에서는 * 혹은 H 둘 다 사용 가능하다. 즉 와 는 같은 의미이다. 양자역학(quantum mechanics)에서는 dagger를 사용)

중요한것은 유니터리 행렬 U의 각 column은 U의 고유벡터(eigenvector)라는 것이다. 그런데 스펙트럼 정리에 의해 A가 U와 실수 고유값행렬 $\Lambda$로 분해할 수 있으므로 U의 고유벡터들은 곧 A의 고유벡터가 된다. 

여기서 다시 한 번 생각해보자. 스펙트럼 정리가 의미하는 것은 어떤 복소수를 가질 수 있는 nxn크기의 정방행렬인 에르미트 행렬 A가 있고, 이 행렬은 실수로 이루어진 고유값행렬 $\Lambda$와 유니터리 행렬인 U와 U의 켤레전치행렬 $U^H$로 분해할 수 있다. 또한 유니터리 행렬의 성질에 의해 U의 역행렬은 U의 켤레복소행렬과 같고, U의 각 column vector는 서로 직교(orthogonal)한다. 이때 U의 column vector들은 A의 고유벡터들이다. 

바로 위에서 정리한 내용들을 보고 이미 눈치채신 분들이 있을 것이다. 그렇다. 스펙트럼 정리는 대칭행렬의 대각화와 매우 유사하다. 식 (4.2)와 비교해보면 더욱 잘 이해할 수 있을 것이다. 에르미트 행렬이 정의하는 것이 복소행렬까지 포함한 것을 제외하면 스펙트럼 정리에서 에르미트 행렬을 대칭행렬로 놓고 이해해도 무리가 없을 정도이다. 결과적으로 우리는 실수대칭행렬(real symmetric matrix)의 대각화를 스펙트럼 정리로 설명할 수 있다. 

스펙트럼 정리는 마치 빛을 서로 다른 파장을 가지는 스펙트럼으로 쪼개서 나타낼 수 있는 것과 같다. 원래는 그냥 하얗게 보이는 빛을 프리즘에 통과시켜보면 여러 색깔로 이루어져 있음을 알 수 있다. 마찬가지로 스펙트럼 정리도 어떤 행렬을 순수 고유값(pure eigenvalue)과 순수고유벡터(pure eigenvector)로 분해하여 표현할 수 있다는 의미에서 빛의 스펙트럼과 그 의미가 같다고 할 수 있다. 

어떤 행렬이 선형 변환(Linear transformation)을 시킴에 있어서 얼마만큼(고유값), 어느 방향으로(고유벡터) 변환시키는 지에 대한 순수 값들로 분해하는 과정. 마치 하얀 빛을 여러 파장의 순수한 빛들로 분해하여 나타내는 것과 같은 느낌으로 받아들이면 좋을 것 같다. 그러한 측면에서 대칭행렬은 정규직교(orthonormal)한 고유벡터와, 역행렬이 전치와 같음으로 인한 계산의 편의성 등 유용한 특성을 지니고 있는 좋은 특성을 지닌 행렬이라고 할 수 있다. 

- Why real eigenvalues?

이제 처음 부분에서 공부했던 대칭행렬의 특성문제로 돌아가보자. 식 (2)의 특성 1은 대칭행렬이 실수인 고유값을 가진다고 했다. 왜 그런지 알아보도록 하자. 

식 (8)은 A가 실수 대칭행렬(real symmetric matrix)일 때를 가정하여 정리한 것이다. 식 (8.1)은 우리가 이미 잘 알고있는 고유값/고유벡터에 대한 식이다. 여기서 양변에 복소켤레(complex conjugate)를 취하면 (8.2)와 같이 된다. 복소켤레는 실수부(real part)는 그대로두고 허수부(imaginary part)의 부호가 반대로 바뀌는 것을 의미하며, 표기는 각 알파벳 기호 위에 bar를 써서 표현한다. 양변에 똑같이 복소켤레를 취해도 여전히 식은 성립하기 때문에 (8.2)와 같이 정리할 수 있다. 

사실 (8.1)에서 (8.2)로 넘어가는 과정이 의미하는 것이 있다. 어떤 실수 행렬 A에 대한 고유값이 복소수 형태라면, 원래의 고유값 와 그의 켤레값(conjugate)인 가 쌍(pair)으로 존재한다. 이는 고유벡터에도 마찬가지로 적용이 되어 와 의 형태로 나타난다. 

그러나 여기서 증명하고자 하는 것은 행렬 A가 실수 대칭행렬일 때 고유값은 항상 실수라는 것이다. 따라서 이를 보이기위해 양변을 전치시켜주면 전치의 규칙에의해 순서가 바뀌어 (8.3)과 같이 정리할 수 있다. 여기서 행렬 A는 실수로 이루어진 대칭행렬이라 가정하였으므로 켤레(conjugate)와 전치(transpose)가 의미가 없어진다. 결과적으로 A의 bar와 전치 기호가 제거된 (8.4)와 같은 형태로 정리할 수 있다. 사실 (8.4)에서 고유값인 λ는 상수이기 때문에 x앞으로 빼서 정리해도 되지만, 전치의 결과를 보이기 위해 바뀐 순서를 그대로 하여 정리하였다. 

이제 식 (8)로부터 고유값이 실수(real number)임을 증명해보자. 어떻게 할 수 있을까? 식 (8.1)과 (8.4)를 이용하면 된다. 아래 식을 보자. 

먼저 식 (8.1)의 양변의 좌측에 $\bar{\textbf{x}}^T$를 곱해주면 (9.1)과 같이 정리할 수 있다. 람다는 상수이기 때문에 앞으로 빼서 정리할 수 있음을 기억하자. 다음으로 (8.4)의 양변의 우측에 x를 곱해주면 (9.2)와 같이 정리가 된다. 이때 (9.1)과 (9.2)의 좌변에 $\bar{\textbf{x}}^T A\textbf{x}$가 존재하는 것을 볼 수 있다. 따라서 식 (9.3)과 같이 (9.1)의 좌변과 (9.2)의 좌변이 같다고 정리할 수 있고, 다시 (9.3)의 양변을 $\bar{\textbf{x}}^T \textbf{x}$로 나눠주게 되면 결과적으로 람다와 람다 bar는 같다고 정의할 수 있다. 즉 복소수인 고유값과 그의 켤레(conjugate)가 같다는 것은 허수부(imaginary part)의 값이 0과 같다는 의미이다. 따라서 이때의 고유값은 항상 실수(real number)라고 할 수 있다. 

정리하자면 실수 대칭행렬인 A가 있을 때, A에 대한 고유값/고유벡터 식을 (8.1)과 같이 정의하고 그의 켤레 전치(conjugate transpose)에 대한 식을 (8.4)와 같이 도출한다. 이때 행렬이 실수 대칭행렬이기 때문에 (8.4)의 A에는 bar와 transpose기호가 없어진다. 이렇게 만들어진 고유값에 대한 식 (8.1)과 (8.4)를 이용하여 양변에 적절한 수를 곱해주어 (9.3)의 관계를 도출하였고, 이를 통해 결과적으로 실수 대칭행렬(real symmetric matrix)은 고유값이 무조건 실수(real number)라는 식 (2)의 조건 1을 증명한 것이다. 

지금까지의 과정을 통해 우리는 실수대칭행렬은 식 (2)에서와 같이 고유값은 언제나 실수(real number)이며 고유벡터들은 서로 수직(perpendicular)임을 증명하였다. 대각화(diagonalization)과정을 통해 실수대칭행렬을 분해하면 일반적인 행렬과는 달리 정규직교(orthonormal)한 고유벡터행렬 Q로 분해할 수 있으며, 이때 Q는 역행렬을 전치(transpose)를 통해 간단히 구할 수 있다. 이처럼 실수 대칭행렬은 좋은 특성들을 가지고 있기 때문에 굳이 이름을 붙이자면 좋은 행렬(good matrix)이라고 하겠다. 바꿔 말하면 좋은 행렬이 되기 위해선 대칭행렬이 가지는 특성들, 실수인 고유값과 직교인 고유벡터의 특성들을 가지고 있어야하며 실은 대칭행렬 그 자체가 되면 된다. 

우리는 지금까지 실수로 이루어진 행렬만을 가정하고 문제를 풀어왔다. 하지만 행렬에는 허수가 포함된 경우도 있을 수 있다. 만약 허수가 포함된 행렬이 있다면, 이또한 좋은 행렬이 될 수 있을까? 다시 말하면 허수가 포함된 대칭행렬도 실수인 고유값을 가지고 고유벡터들이 서로 직교(orthogonal)할 수 있을까? 결론부터 말하자면 어떤 허수가 포함된 대칭행렬(complex symmetric matrix) A가 이와 같은 특성을 가지기 위해선 A의 켤레전치행렬(conjugate transpose matrix)이 원래의 A와 같아야 한다는 조건을 만족시켜야 한다. 

좋은 행렬(good matrix, 실수 고유값, 직교 고유벡터)이 되려면... 

위 테이블의 좌측은 우리가 지금까지 다뤄왔던 A가 실수로만 이루어진 경우에 좋은 행렬이 되기 위한 조건이다. 반면에 오른쪽은 A가 허수를 포함한 경우에 좋은 행렬이 되기 위한 조건을 나타낸다. 원래의 행렬과 켤레전치행렬이 같은 경우에 좋은 행렬이 될 수 있다. 그런데 테이블의 오른쪽 허수를 포함한 행렬은 이번 강의에서 이미 다룬 것이다. 바로 식(6)의 에르미트 행렬(Hermitian matrix)이다. 행렬의 원소를 복소수가 포함된 경우까지 확장하여 생각해보면 에르미트 행렬인 경우에 좋은 행렬이 될 수 있다는 결론이 나온다. 에르미트 행렬이 되기 위해선 켤레전치행렬이 같아야 하고, 켤레전치행렬이 같기 위해선 서로 대응되는 원소, 즉 위의 경우엔 (1, 2)의 1+i와 (2, 1)의 1-i가 서로 켤레(conjugate)의 관계여야 한다. 물론 실제 응용문제에서는 거의 대부분이 실수대칭행렬문제이긴 하지만 그래도 실수대칭행렬과 에르미트 행렬의 차이를 잘 알아두도록 하자. 

3. 대칭행렬의 그 외의 특징들

- projection matrices in symmetric matrix

대칭행렬을 해석하는 또 다른 관점이 있다. 이미 알다시피 대칭행렬은 대각화를 통해 Q와 Λ의 곱으로 분해할 수 있다. 이 식을 다시 써보자. 

식 (10.1)은 분해된 행렬의 원소들이 곱해지는 과정을 자세히 풀어서 작성한 것이며, Lecture 3에서 배웠던 행렬 곱셈에 따라 column * row의 곱의 조합으로 나타낸 것이다. (10.2)와 (10.3)은 (10.1)의 우변에 해당하는 각각의 항들을 보다 이해하기 쉽게 시각화한 것이다. 이미 배운대로 column * row의 순으로 벡터를 곱하면 하나의 행렬이 만들어지는데, 여기에선 같은 벡터끼리 곱해줬기 때문에 (10.3)과 같이 nxn크기의 정방행렬이 나올 것이다. 모든 대칭행렬들이 이와 같은 조합(combination)의 형태로 표현될 수 있다. 

그렇다면 (10.3)에서 표현된 각 람다와 곱해진 정방행렬들은 무엇일까? 바로 상호간 서로 독립인 투영행렬(mutually independent projection matrix)이다. 투영행렬에 대한 내용은 Lecture 15를 참고하자. (10.3)의 각 투영행렬은 n개 만큼 나오고, rank는 1이다. 이는 대칭행렬 A가 nxn의 정방행렬이고 역행렬의 계산이 가능하다면 이때의 대칭행렬은 full rank일 것이고, 따라서 n개의 투영행렬의 조합이 가능하다면 각 투영행렬은 rank=1이 될 것이다. 이는 결국 스펙트럼 정리의 측면에서 보자면 nxn크기의 대칭행렬 A는 n개의 투영행렬과 고유값의 조합으로 표현이 가능하다는 것을 알 수 있다. 

대칭행렬이 실제로 어떻게 투영행렬의 조합으로 표현될 수 있는지 그래프를 통해 알아보도록 하자. 먼저 임의의 대칭행렬의 고유값과 고유벡터를 구해보자. 

(11.1)은 대칭행렬, (11.2)는 고유값, (11.3)은 아직 정규화(normalization)되지 않은 상태의 고유벡터를 각각 나타낸다. 자세한 풀이 방법은 Lecture 21-(2)을 참고하자. 다음으로 고유벡터를 정규화 한 뒤, (10.3)과 같이 투영행렬의 조합으로 표현해보자. 

(12.1)은 (10.1)의 형태로, (12.2)는 (10.3)의 형태로 정리한 것이다. 여기서 고유값 바로 뒤에 곱해진 정방행렬들이 바로 투영행렬들이다. 이제 (11.1)의 대칭행렬 A의 선형변환을 그래프로 표현해보자.  

Fig. 2 식 (11)의 대칭행렬의 선형변환 그래프

파란색 벡터는 변환 전, 빨간색 점선 벡터는 A에 의해 변환된 벡터를 의미하며, 벡터의 끝 부분에는 변환 전 원래 벡터의 좌표가 표시되어 있다. 변환 전과 변환 후의 벡터의 모습을 보고 대략적으로 이 행렬이 어떤식으로 변환을 시키는지를 유추할 수 있다. 변환의 방향은 고유벡터, 변하는 정도는 고유값이 영향을 미친다는 것을 알아두자. 

그렇다면 (12.2)에서 정리한 각 고유벡터의 투영행렬을 이용하여 위의 파란색 벡터를 투영시키면 어떤 모습이 될까? 바로 아래 그림과 같은 모습이 될 것이다. 

Fig. 3 (12.2)의 첫 번째 투영행렬을 이용하여 투영시킨 모습

Fig. 4 (12.2)의 두 번째 투영행렬을 이용하여 투영시킨 모습

Fig. 3과 4는 각각 (12.2)의 첫 번째와 두 번째 투영행렬을 이용하여 원래의 파란벡터들을 투영시킨 모습이다. 결과적으로 파란 벡터들이 투영행렬 A의 고유공간(eigenspace)으로 투영(projection)된 것을 볼 수 있다. 이것을 그래프상에서 해석해보면 Fig. 3의 투영된 벡터와 Fig. 4에서 투영된 벡터를 더하면 Fig. 2의 대칭행렬로 변환된 벡터와 일치한다는 것이다. 즉 예를 들면 [1, 1]의 벡터를 첫 번째 투영행렬로 투영시킨 벡터와, 두 번째 투영행렬로 투영시킨 벡터를(고유값의 곱도 포함함) 서로 더하면 [1, 1]을 A로 변환시킨 벡터와 같다는 것이다. 아래 그림은 [1, 1]의 투영에 대한 예시이다. 

Fig. 5 (11.1)의 대칭행렬을 이용한 [1, 1]벡터의 투영 예시

Fig. 5의 파란색 벡터는 [1, 1]의 변환 전 벡터를 나타내고 빨간색 벡터는 첫 번째 투영 벡터, 두꺼운 자홍색 벡터는 두 번째 투영벡터를, 녹색 벡터는 대칭행렬 A에 의해 변환된 벡터를 각각 의미한다. 얇은 자홍색 벡터는 첫 번째 투영벡터와 두 번째 투영벡터를 더한 결과를 나타낸다. 더한 결과가 일치하는 것을 볼 수 있다. (대칭행렬일 때만 항상 성립함을 주의하자) 이것을 스펙트럼 정리에 적용시켜 이해해보면 각각의 투영행렬은 대칭행렬의 스펙트럼(spectrum)이라고 이해해도 무리가 없다. 

이를 통해 우리는 다음의 결론을 내릴 수 있다. 

아래는 MATLAB코드이다. 

- pivots and eigenvalues

마지막으로 살펴볼 대칭행렬의 특징은 피벗(pivot)에 관한 것이다. 정확히는 피벗과 고유값 사이의 관계에 관한 내용이다. 얼핏 생각하면 피벗과 고유값은 별 관계가 없을 것 같지만 한 가지 공통점이 있다. 바로 부호(sign)와 관련된 것인데, 바로 대칭행렬에서 피벗과 고유값은 같은 부호를 갖는 원소의 개수가 같다는 것이다. 쉽게 예를 들어 설명하면 100x100크기의 정방행렬이 있을 때, 피벗과 고유값의 개수도 똑같이 100개가 될 것이다. 여기서 100개의 피벗 중 56개의 피벗이 양수이고, 나머지 44개의 피벗이 음수라고 가정하자. 이때 고유값도 마찬가지로 100개 중에 56개의 고유값이 양수, 나머지 44개의 고유값이 음수가 된다. 이것이 대칭행렬에서 피벗과 고유값의 부호에 관한 관계이다. 아래의 식을 보자. 

식 (13.1)의 2x2 행렬은 양수와 음수 피벗이 한 개씩 존재한다. 마찬가지로 고유값도 양수와 음수 한 개씩 존재하는 것을 볼 수 있다. (13.2)의 경우엔 피벗이 둘 다 양수이고 고유값도 마찬가지로 둘 다 양수임을 볼 수 있다. 

그렇다면 위와 같은 피벗과 고유값 사이의 부호 관계는 왜 알아야할까? 일단 Lecture 23-(1)에서 공부했던 미분방정식에서 고유값의 부호가 해당 시스템의 안정성(stability)에 영향을 미치는 등 중요한 역할을 한다는 것을 이미 배운바있다. 고유값의 부호가 중요하다는 것은 알겠는데, 피벗과는 왜 연관지어서 알아야할까? 

지금까지 우리가 다루어왔던 행렬들은 그 크기가 2x2, 3x3, 아무리 커 봐야 4x4 정도였다. 하지만 어떤 시스템에서는 50x50, 100x100, 혹은 그 이상의 크기의 행렬을 다루어야 하는 경우도 발생할 수 있다. 우리가 배웠던 고유값을 구하는 방식은 det(A-λI)의 꼴로 놓고 λ에 대한 다항식(polynomial)을 만들어 푸는 방식이다. 그러나 차수가 증가함에따라 고차다항식을 풀어서 고유값을 구해야 하는데, 고차로 갈 수록 해(solution)의 불안정성이 커지게 된다. 즉 100x100의 고유값을 구하려면 100차 다항식을 풀어야 한다는 소리다. 

이와 같은 커다란 행렬의 고유값을 구하려면 기존의 방법보다는 수치선형대수(numerical linear algebra)의 방법을 이용해서 값을 구하는 것이 훨씬 안정적이다. 즉 피벗 같은 경우엔 MATLAB같은 프로그램으로 비교적 안정적으로 해를 구할 수 있으므로 우선 피벗을 구해서 양의 개수가 몇 개인지, 음의 개수가 몇 개인지를 안 다음, 거기서부터 수치해석적으로 고유값을 구해나가는 방식을 사용하는 것이다. 

자세한 것을 다루진 않겠지만 고유값을 수치적으로 구하기 위해서 사용되는 중요한 사항이므로 이러한 것이 있다는 것 정도만 알아두도록 하자. 

4. 마치며

이번 강의에선 대칭행렬의 여러 가지 특성과 스펙트럼 정리에 대해 알아봤다. 기본적으로 대칭행렬은 실수인 고유값, 그리고 서로 수직인 고유벡터들을 가지며 이를 통해 행렬을 대각화하여 분해할 때 전치연산으로 간단히 역행렬을 구할 수 있는 등 좋은 특성들이 있음을 배웠다. 이러한 좋은 특성들은 특이값 분해(SVD), 주성분 분석(PCA)등의 연산에 기본이 되는 성질이다. 또한 대칭행렬은 스펙트럼 정리로 설명할 수 있으며, 행렬 원소를 복소수까지 확장한 에르미트 행렬에 대해서도 공부하였다. 이번 대칭행렬은 앞으로의 강의에서 다룰 positive definite matrix에 대한 준비과정이므로 잘 공부하도록 하자. 

이번 포스팅에서는 마코브 행렬(Markov matrix)에 대해 다루도록 하겠다. 마코브 행렬은 이전 강의에서 다루었던 행렬의 대각화 Lecture 22와 관련이 깊기 때문에 앞선 포스팅을 먼저 학습하길 바란다. 

1. 마코브 행렬과 마코브 체인

- What is the Markov matrix and Markov chain?

마코브 행렬(Markov matrix)은 1906년에 러시아의 수학자 Andrey Markov에 의해 처음으로 언급된 개념이다.  마코브 행렬은 마코브 체인(Markov chain)을 기술하기 위한 수학적 도구인데, 확률적 방법을 기반으로 하기 때문에 확률 행렬(stochastic matrix or probability matrix)로 불리기도 한다. 

여기서 마코브 체인은 확률을 이용하여 어떤 객체 상태를 시간에 따라 어떻게 변화할지를 모델링(modeling)하는 것이다. 말이 약간 어려워 보이지만 쉽게 말하자면 날씨 예측, 인구 이동 예측 등과 같이 어떤 객체(날씨, 인구)의 상태(맑음, 흐림, 10만명, 5만명, ...)가 시간이 지남에 따라 어떻게 변화할지를 확률을 이용하여 예측하는 것이다. 체인(chain)이라는 단어가 의미하듯이 객체의 시간에 따른 서로 다른 상태를 어떻게 연결할지를 기술하는 것이 마코브 체인이며, 이들을 연결시켜주는 매개체 역할을 마코브 행렬이 하는 것이다. 

마코브 체인(Markov chain)은 마코브 프로세스(Markov process)와 혼용되기도 하는데, 이 둘을 구분 짓는 정확한 정의는 없다. 다만 일반적으로 객체의 상태 공간이 이산 시간(discrete time)이면 마코브 체인, 연속시간(continuous time)이면 마코브 프로세스로 정의한다. 

마코브 행렬은 다른 말로 전이 행렬(transition matrix), 또는 치환 행렬(substitution matrix)이라고도 한다. 마코브 행렬이 왜 이렇게 불리는지 아래의 그림을 통해 알아보자. 

Fig. 1 삼성 갤럭시 핸드폰과 애플 아이폰의 시장 점유율 예측에 관한 마코브 체인 상태 전이 다이어그램(Markov chain state transition diagram)

Fig. 1은 삼성 갤럭시 핸드폰과 애플의 아이폰에 대한 시장 점유율을 마코브 체인으로 표현한 것이다. 물론 수치는 임의로 정한 것들이다. Fig. 1과 같은 도표를 마코브 체인 상태 전이 다이어그램(Markov chain state transition diagram)이라고 한다. 마코브 체인은 기본적으로 각 상태들이 다른 상태로의 전이(transition)가 얼마의 확률로 이루어질 것인가를 기술한다. 

애플의 아이폰을 사용하는 것을 state 1, 삼성의 갤럭시 핸드폰을 사용하는 것을 state 2로 각각 정의했고, 각 상태가 다른 상태로 전이되는 행위를 화살표로 표시하였다. 상태전이는 자기 자신으로 될 수도 있는데, 아이폰의 경우 상태(state)가 자기 자신으로 전이되는 것은 현재 아이폰을 사용하고 있는 유저가 계속 아이폰을 사용하는 것을 의미하며, 그 확률을 71.4%로 정의하였다. 반면 아이폰을 쓰던 유저가 갤럭시폰으로 갈아탈 확률은 28.6%라고 정의하였다. 이것을 state 1에서 state 2로 상태가 전이 될 확률로 해석할 수 있다. 

마찬가지로 삼성의 갤럭시폰을 쓰던 유저가 계속 갤럭시 폰을 사용할 확률은 63.7%, 아이폰으로 갈아탈 확률은 36.3%로 각각 정의하였다. 이것이 마코브 체인 상태 전이 다이어그램이다. 

우리가 이와 같이 마코브 체인을 구성한 이유는 앞선 문제의 경우엔 결국 핸드폰 사업에서 두 회사의 앞으로의 시장 점유율을 예측하기 위함이다. 마코브 방법을 이용하여 앞으로의 상태를 예측하기 위해선 먼저 위 그림과 같이 마코브 체인을 정의하고, 현재 시간의 상태와 그 다음 시간의 상태 사이를 연결시켜줄 수 있는 매개체인 마코브 행렬을 구해야한다. 마코브 행렬은 마코브 체인에서 각 상태로 전이될 확률만 나오면 쉽게 구할 수 있다. Fig. 1의 마코브 체인에 대한 마코브 행렬은 아래와 같다. 

식 (1)의 좌측의 행렬 M이 바로 마코브 행렬(Markov matrix)이다. 행렬의 각 원소는 어떤 상태에서 자기 자신, 혹은 다른 상태로 전이될 확률값을 가진다. 

행렬을 해석할 땐 column의 인덱스를 시작 상태로 보고, row의 인덱스를 전이 되는 상태로 보면 된다. 즉 M의 첫 번째 column인 Apple을 시작 상태로 봤을 때, 다시 Apple이 될 확률을 row 1로 보고 이때의 값은 0.714가 된다. Apple에서 Samsung으로 전이 될 확률은 row 2가 되며, 이때의 값은 0.286이다. 마찬가지로 두 번째 column인 Sam을 시작 상태로 봤을 때, Apple로 상태가 전이 될 확률은 row 1이 되고 그 값은 0.363이 되며, 다시 자기 자신인 Sam으로 될 확률은 row 2이고 값은 0.637이다. 결국 어떤 시작 상태에서 다시 자기 자신이 될 확률은 column과 row가 같은 인덱스를 가질 때이며, 각 column의 인덱스를 시작 상태로 하여 다른 상태로 전이될 확률을 보고 싶으면 각 column에 해당하는 row의 원소들을 보면 된다. 

식 (1)의 우측에 있는 u는 Fig. 1에서 Apple과 Samsung의 현재의 시장 점유율을 의미한다. 주어진 마코브 행렬의 우측에 u를 곱한 결과는 다음 시점에서의 시장 점유율의 예측값이 된다. 

참고로 필자가 마코브 행렬을 표현한 방식은 column을 기준으로 만든 행렬인데, 어떤 사람들은 row를 기준으로 행렬을 구성하기도 한다. 즉 row인덱스가 시작 상태가 되고 column이 전이될 상태를 의미하는 것이다. 이 경우 마코브 행렬은 column 기준 마코브 행렬이 전치(transpose)된 형태가 되며, 현재 상태인 u도 역시 전치된 row 벡터 형태로써 M의 좌측에 곱해진다. 어떤 방법이 더 옳다고 할 순 없지만, 필자는 column 기준 방식을 선호하기 때문에 이번 강의에서는 column 기준 마코브 행렬을 사용하도록 하겠다. 

이제 마코브 체인(Markov chain)과 마코브 행렬(Markov matrix)이 어떤 것을 의미하는지 대략적으로 파악했을 것이다. 자세한 풀이는 이후에 하도록 하고 먼저 마코브 행렬이 가지는 특성에 대해서 알아보도록 하자. 

2. 마코브 행렬의 특성

- Properties of Markov matrix

아래 행렬은 임의의 3x3크기의 마코브 행렬을 만든 것이다. 

식 (2)에 표현된 마코브 행렬에서 어떤 특징이 보이는가?  

마코브 행렬의 첫 번째 특징은 모든 원소의 값이 0보다 크거나 같다는 것이다. 즉 오직 양수만 허용된다. 마코브 행렬은 기본적으로 상태들의 전이 확률(transition probability)을 나타내기 위한 행렬이므로 확률값을 원소로 가진다. 확률 값에는 음수가 있을 수 없기 때문에 마코브 행렬의 원소값들은 반드시 0보다 크거나 같은 양수만 올 수 있다. 또한 마코브 행렬은 시간이 지남에 따라 변화하는 상태를 기술 및 예측하기 위한 행렬이므로 행렬의 거듭제곱이 발생한다. 이때 거듭제곱의 결과값도 모두 양수가 된다.

두 번째 특징은 각 column의 원소의 합은 1이 된다는 것이다. 식 (2)의 column 1의 원소들을 모두 더해보자. 그 결과값은 1이 되고, column 2, column 3도 마찬가지이다. 이 두 번째 특징은 마코브 행렬을 거듭제곱해도 마찬가지로 성립한다. 즉 M2=MxM 이라고 했을 때(M은 마코브 행렬을 의미), M2의 각 column 원소들의 합은 역시 1이 되는 것이다. column의 원소의 합이 1이 되어야 하는 이유는 하나의 column은 어떤 상태에 대한 모든 전이 확률을 나타내기 때문이다. 즉 column 1이 자기 자신이 될 확률(row1), 상태 2(row2)나 상태 3(row3)이 될 확률을 각각 기술하고 있으며, 따라서 하나의 상태에 대한 모든 전이 확률을 나타내는 것이기 때문에 column의 합은 1이 되어야한다. 

이 두 가지 특성을 정리하면 아래와 같다. 

모든 원소들이 0보다 크거나 같아야 한다. 모든 각 column원소들은 더했을 때 1이 되도록 해야한다. 이것이 마코브 행렬(Markov matrix)이 가지는 기본적인 특성이다. 

- Eigenvalues of Markov matrix and steady state

마코브 행렬은 고유값에 대한 한 가지 특성을 가진다. 이 특성을 이해하기 위해선 정상 상태(steady state)에 대해서 먼저 이야기해야 한다. 정상상태는 지난 미분방정식과 선형대수에 대한 강의 Lecture 23-(1)에서 언급된 바 있다. 해당 강의에서 우리는 미분방정식에서 해가 시간이 지남에 따라 크게 세 가지 형태가 될 수 있다고 하였으며, 그 형태는 각각 안정 상태(stability), 정상 상태(steady state), 발산(divergence)이다. 안정 상태가 되려면 행렬의 모든 고유값이 0보다 작아야하고, 정상 상태일때는 하나의 고유값이 0이고 나머지가 0보다 작아야 한다. 그리고 어느 하나의 고유값이라도 0보다 크면 그 해는 발산한다. 미분방정식의 해가 고유값에 따라 세 가지 형태를 띄는 이유는 미분방정식의 일반해(general solution)의 형태를 보면 쉽게 유추할 수 있다. 해당 강의의 일반해에 대한 식을 다시 써보자. 

식 (4)은 행렬로 정의된 미분방정식의 일반해에 대한 식을 나타낸다. 보다시피 일반해는 계수값 c1, c2, 지수 함수(exponential function)의 지수부에는 고유값이, 그리고 고유벡터 x1, x2의 선형 조합(linear combination)으로 구성되어 있다. 계수값 c와 고유벡터 x는 상수이기 때문에 고정되고, 변수 t와 곱해지는 고유값이 해(solution)에 영향을 미치게 되는데, 고유값이 0인 경우 지수 함수가 1이 되기 때문에 극한의 시간으로 가도 그 값은 결국 상수 c와 x의 곱으로 유지가 된다. 고유값이 0보다 작을 경우엔 시간이 지날 수록 값이 작아져서 결국엔 0이 된다. 이것이 하나의 고유값이 0이고, 나머지 고유값이 0보다 작을 경우 정상 상태가 되는 이유이다. 

이에 대한 자세한 사항은 해당 강의를 참조하자. 

여기서 우리가 관심을 가져야 하는 상태는 바로 정상 상태(steady state)이다. 그 이유는 마코브 행렬(Markov matrix)이 정상 상태(steady state)의 특성을 가지기 때문이다. 여기서 "마코브 행렬이 정상 상태의 특성을 가진다고? 그렇다면 마코브 행렬의 고유값 중 하나는 0의 값을 가지겠군!" 라고 생각하는 사람들이 있을 것이다. 그러나 이건 잘못된 생각이다. 미분방정식의 경우엔 일반해에서 고유값의 위치가 지수함수의 지수부이다. 따라서 지수함수를 시간에 상관없이 1로 만들기 위해선 고유값이 반드시 0이어야 한다. 

그러나 마코브 행렬의 경우엔 미래의 일을 예측하기 위해서는 행렬을 거듭제곱(power)해야 한다. 행렬을 거듭 제곱 함에 있어서 행렬의 고유값이 0이라면, 그리고 나머지 고유값들의 크기가 1보다 작다면 이 행렬의 거듭제곱의 결과값은 얼마 못가 0이 되고 만다. 행렬의 거듭제곱을 함에 있어 해당 행렬이 정상 상태의 특성을 가지기 위해선 아래의 두 가지 조건을 충족시켜야 한다. 

행렬의 거듭제곱에서 정상 상태의 특성을 보이기 위해선 적어도 하나의 고유값은 1이어야 한다. 그리고 다른 모든 고유값들의 크기는 1보다 작아야한다. 마코브 행렬은 이 두 가지 조건을 충족시키는 행렬이고, 결과적으로 행렬의 거듭제곱(power)에 있어서 정상 상태(steady state)의 특성을 갖는다. 조건 (5)에 대한 내용을 MATLAB을 통해 검증해보자. 

Fig. 2 고유값에 따른 행렬의 거듭제곱 결과. [Upper left] 고유값=1일 때. [Lower right] 고유값=0일 때

Fig. 2는 고유값에 따라 행렬의 거듭제곱 결과가 어떻게 나오는지를 MATLAB을 통해 검증한 결과이다. 왼쪽 상단(파란색 영역)은 하나의 고유값이 1이고 나머지 고유값의 크기가 1보다 작을 때 행렬을 100번을 거듭제곱 한 결과이다. rst=로 표현된 값이 행렬을 거듭제곱한 결과인데, 어떤 특정 값으로 수렴하는 모습을 보인다. 따라서 조건 (5)를 만족시키는 마코브 행렬은 정상 상태(steady state)의 성질을 가짐을 증명하였다. 

행렬의 거듭제곱의 결과값인 rst=로 표시된 값 아래에 LC1=로 표시된 결과가 나오는데, 앞서 계산한 거듭제곱의 결과(rst=)와 그 값이 같은 것을 볼 수 있다. LC1=로 표현된 값은 소스코드의 맨 아래 줄에 의해 계산된 값인데, 어디서 많이 본 기억이 날 것이다. 바로 지난 강의 Lecture 22의 행렬의 대각화에서 배운 내용이다. 행렬의 대각화 강의에서 우리는 계차방정식(difference equation)의 해를 행렬의 대각화를 통해 구했는데, 그때 행렬을 일일이 거듭제곱을 하지 않고서도 행렬의 고유값과 고유벡터, 그리고 계수의 선형조합의 형태로 해를 효율적으로 구하는 방법을 배웠다. 잘 기억나지 않는 분들을 위해 식을 다시 써보면 아래와 같다. 

식 (6)과 같이 행렬의 거듭제곱을 일일이 계산할 필요가 없다. 단지 행렬의 고유값과 고유벡터를 알면 훨씬 적은 계산량으로 구할 수 있으며, LC1이 그 결과값이다. 여기서 중요한 포인트는 계수 c1, c2와 고유벡터 x1, x2는 상수이기 때문에 고정된 값이고, 우리가 주목해야 할 값은 고유값이라는 것이다. 즉 고유값이 어떻게 나오느냐에 따라서 해당 행렬이 발산하는지, 정상 상태(특정 값으로 수렴)인지, 0으로 수렴하는지가 결정이 되는데 그 이유는 식 (6)과 같이 나머지는 상수이고 고유값만 k의 거듭제곱으로 표현이 되기 때문이다. Fig. 2에 나온 문제를 예로 들자면 첫 번째 고유값은 1이고, 100제곱을 하던 1000제곱을 하던 그 결과는 항상 1이다. 반면 두 번째 고유값은 그 크기가 1보다 작은 0.351이기 때문에 거듭제곱을 할 수록 그 값은 작아져서 0에 가까워지게 된다. 

식 (7)을 통해 말할 수 있는 것은 Fig. 2의 M1 마코브 행렬의 정상 상태가 바로 (7.2)인 c1x1이라는 것이다. 결국 마코브 행렬의 정상 상태(steady state)는 λ=1에 대응되는 고유벡터(eigenvector)라 할 수 있다. 

반면 Fig. 2의 오른쪽 하단(붉은색 영역)의 경우를 보면 마코브 행렬 M2가 0인 고유값을 가지는 것을 볼 수 있다. 또한 나머지 고유값의 크기도 1보다 작기 때문에 결과적으로 0행렬의 결과를 보인다. 이를 통해 행렬을 거듭제곱 함에 있어 정상 상태(steady state)의 특성을 보이기 위해선 조건 (5)를 만족시켜야 함을 알 수 있다. 

물론 마코브 행렬이면서 위의 조건을 충족시키지 않는 예외 적인 경우도 있는데, 바로 단위 행렬(identity matrix)이 그 경우다. 조건 (3)을 만족시키기 때문에 마코브 행렬이라 할 수 있지만, 정상 상태를 위한 조건 (5)는 만족시키지 않는다. 하지만 거듭 제곱(power)을 해도 발산하거나 0행렬이 되지 않기 때문에 결과적으로는 정상 상태라고 할 수 있다. 예외 적인 경우이니 잘 알아두도록 하자. 

- Eigenvalue 1

앞서 우리는 마코브 행렬이 정상 상태의 특성을 가지기 위해선 반드시 하나의 고유값이 1이어야 함을 배웠고, 이를 계차방정식을 이용하여 이를 설명하였다. 하지만 계차방정식 이외에도 고유값이 1이어야 하는 것을 설명할 수 있다. 바로 조건 (3)-2의 각 column의 합은 1이 된다는 조건으로부터 고유값이 1이어야 한다는 것을 설명할 수 있다. column의 합이 1이 된다는 조건과 고유값이 1인 것이 무슨 관계지? 라는 의문이 들 것이다. 하나씩 풀어나가보자. 앞서 만들었던 임의의 마코브 행렬 (2)를 다시 써보자. 

식 (2)는 각 원소의 크기가 0보다 크고 각 column의 합이 1이 되는 조건을 충족하는 마코브 행렬이다. 이제 이 행렬의 고유값을 구한다고 생각해보자. 고유값을 구하기 위해선 det(A-λI)를 계산해야 하는데, 우리는 A의 고유값이 λ=1임을 가정했기 때문에 (A-I)를 계산해보도록 하겠다. 

식 (8)은 A-I를 계산한 모습을 나타낸다. 원래 행렬 A의 특징 중 하나가 각 column의 합이 1이 되는 것이었는데, 단위 행렬을 빼고 나니 각 column의 합이 0이 되는 특징을 보인다. 이때 이 column의 합이 0이 되는 특성이 의미하는 것은 곧 A-I가 특이 행렬(singular matrix)임을 의미한다. 어째서 그럴까? 바로 0벡터 이외의 null space가 존재하기 때문이다. 이 null space가 존재하는지를 계산해보지도 않고 어떻게 바로 알 수 있을까? 그 힌트는 바로 식 (8)에 이미 나와있다. 

식 (8)의 행렬 A-I는 모든 각 column이 더했을 때 0이 된다. 이는 바꿔말하면 A-I 행렬의 row들이 dependent함을 의미한다. 즉 A-I는 transpose를 시켰을 때 null space가 존재한다는 의미다. Lecture 10을 공부했다면 A-I의 transpose의 null space가 Left Null Space라는 것을 알 것이다. 식 (8)의 Left null space를 정리하면 아래와 같다. 

식 (9)에 보이는 것 처럼 0벡터 이외의 $(A-I)^T \boldsymbol{x}=0$를 만족시키는 벡터가 존재하므로 (A-I)는 Left null space가 존재한다. Left null space의 값은 column을 다 더했을 때 0이 된다는 사실로부터 쉽게 유추할 수 있는데, 이 상태에서 전치(transpose)를 시켰을 때 모든 column을 그대로 더하면 0벡터가 나옴을 알 수 있고, 따라서 각 column에 계수 1을 곱해주면 된다는 결론이 나온다. 이렇게 되면 Left null space의 값은 식 (9)에서와 같이 $[1 \; 1 \; 1]^T$ 가 되어야만 한다. 

한편으로 Left null space가 존재한다는 의미는 곧 (A-I)는 full rank가 아니라는 의미가 되며, (A-I)의 column이 dependent하다는 결론을 내릴 수 있다. 따라서 아래와 같이 (A-I)의 null space도 존재하게 된다. 

식 (10)에서 (A-I)의 null space를 계산하였다. 그런데 생각해보면 식 (10)의 null space를 구하는 식은 사실 우리가 가정했던 "마코브 행렬 A는 고유값 λ=1을 가진다"는 사실로부터 출발하였고, (10)에서 구한 고유벡터는 A의 λ=1에 대응되는 고유벡터이다. 이것이 곧 A의 정상 상태(steady state)라고 할 수 있다. 결과적으로 마코브 행렬의 각 column은 더했을 때 그 크기가 1이 된다는 성질로 부터 마코브 행렬의 최소 하나의 고유값은 1이 된다는 사실을 연결지어 설명하였다. 설명이 조금 길어져서 잘 이해가 안될 수 있기 때문에 다시 한 번 정리해보자. 

또 다른 관점에서는 A의 행렬식의 성질로부터 이를 유추할 수 있다. Lecture 18에서 행렬식의 특성 10번으로 부터 A와 A의 transpose는 같은 고유값을 가진다는 것을 알 수 있다. A transpose가 고유값 1을 가질 때, 앞에서와 같이 null space [1, 1, 1]T를 가진다. 즉 A입장에서는 left null space를 가지는 것이다. 그런데 A는 A transpose와 같은 고유값을 가지므로 λ=1을 가질 것이고, 그에 따른 고유벡터도 존재하게 된다. 여기서 고유값 1에 대응되는 고유벡터가 마코브 행렬 A의 정상 상태이다. A와 A transpose는 같은 고유값을 가질 수는 있으나, 같은 고유값에 대응되는 각각의 고유벡터는 다를 수 있음에 유의하자. 

3. 마코브 행렬의 응용

- i-Phone vs Galaxy S

이제 앞서 제시했던 문제를 풀어보도록 하자. Fig. 1에서 우리는 애플의 아이폰과 삼성의 갤럭시 핸드폰의 시장 점유율 예측에 대한 가상의 문제를 마코브 체인으로 모델링 하였다. 또한 마코브 체인에 대한 마코브 행렬을 식 (1)에 정리하였다. 일단 이 마코브 행렬을 다시 써보자. 

앞서 언급한 것과 같이 식 (1)의 M은 아이폰과 갤럭시폰 유저들의 상태가 전이될 확률을 표현한 마코브 행렬, u는 현재의 상태(시장점유율)를 나타낸다. 만약 상태 전이에 대한 예측 단위가 한 달 단위로 이루어진다고 가정해보자. 올해 1월의 시장 점유율은 아이폰이 49.2%, 갤럭시가 50.8%를 차지하고 있다. 그렇다면 2월의 시장점유율을 예측하려면 어떻게 해야할까? 이미 모델링 해놓은 마코브 행렬을 u의 좌변에 곱해주면 된다. 실제 계산을 해보면 2월의 시장점유율 예측치는 아래와 같이 계산할 수 있다. 

Fig. 3 아이폰과 갤럭시폰의 마코브 체인

Fig. 3은 시장점유율의 계산을 위한 마코브 체인을 나타낸다. 1월의 시장점유율 벡터인 [0.492, 0.508]T를 마코브 행렬에 곱해주면 그 결과는 2월의 시장점유율이 되고 값은 [0.535, 0.464]T가 된다. 다시 2월의 시장점유율을 마코브행렬과 곱해주면 3월의 시장점유율에 대한 예측치가 나온다. 아래 식은 계산 과정을 나타낸다. 

식 (11)에 표현된 식들을 보고 약간 익숙하게 느끼는 분들이 있을 것이다. 바로 계차방정식 $u_{k+1}=Au_k$ 와 같은 형태이다. 이 말은 마코브 체인을 계차방정식으로 생각할 수 있고, 결국 고유값과 고유벡터의 선형조합의 형태로 표현하여 해를 구할 수 있다는 의미다. 우리는 이미 식 (6)에서 계차방정식의 해에 대한 식을 보았으므로 여기에 필요한 고유값과 고유벡터를 구해보자. 

기존에 고유값을 구할 때 행렬식(determinant)을 이용했는데, 식 (12.1)에서는 행렬식을 이용하지 않고 고유값의 합은 행렬의 trace와 같다는 성질을 이용하여 간단히 구했다. 행렬 M이 마코브 행렬임을 이미 알고있기 때문에 첫 번째 고유값은 1임을 미리 알 수 있으며, 나머지는 M의 trace에서 1을 빼주면 간단히 구할 수 있다. 물론 이 trace성질은 2x2행렬에 국한된 것임을 유의하자.  

식 (12.2)와 (12.3)에서는 각각 고유벡터를 구한 뒤, 정규화(normalization)를 해주었다. 이제 필요한 값들을 모두 구했으니 해를 구해보자. 

식 (13.2)에서는 계수값 c를 계산하였다. 이는 식 (13.1)에서 k=0으로 놓고 계산한 것이다. 식 (13.3)에서는 계수, 고유값, 고유벡터의 실제 값을 대입하여 정리하였다. 이제 이 식을 기반으로 향후 시장 전망을 예측할 수 있다. 물론 몇 step까지는 행렬을 일일이 곱하여 계산할 수도 있다. 하지만 이는 행렬의 크기가 커질 경우 속도가 매우 느리고 비효율적이다. 만약 어떤 마코브 체인 문제에서 마코브 행렬의 크기가 10000 x 10000인데 극한(infinite)에서의 결과값을 알고싶다고 하자. 실제 극한을 구할 수는 없으니 예를 들어 10,000,000 step이후의 상황을 예측한다고 가정해보자. 이 경우 직접 행렬을 곱해서 계산을 한다는 것은 매우 무식한 방법이다. 식 (13.3)과 같이 해를 계산한다면 단지 고유값의 거듭제곱만 계산하면 되기 때문에 훨씬 빠르고 효율적으로 답을 구할 수 있게 된다. 

이제 위의 식을 이용하여 아이폰과 갤럭시폰의 한 해 동안의 시장점유율을 예측해보자. 아래 그림은 마코브 체인으로 모델링한 애플과 삼성의 핸드폰 시장점유율 예측 그래프이다. 

Fig. 4 마코브 체인을 이용한 애플과 삼성의 휴대폰 시장점유율 예측 그래프

초기 점유율은 삼성이 더 높았지만, 마코브 체인을 통한 시뮬레이션 결과 급격하게 점유율이 역전하면서 약 4월달 정도에는 정상상태(steady state)로 수렴하는 모습을 보인다. 

한 가지 알아둬야 할 것은 마코브 체인은 그 한계점이 극명하게 존재한다는 것이다. 위의 문제에선 매달 시장점유율을 예측함에 있어서 그 달의 점유율에 마코브 행렬을 곱하여 다음달의 점유율에 대한 예측을 하였다. 중요한 점은 그 과정에서 마코브 행렬 자체는 한 번도 변하지 않았다는 것이다. 위 문제에서 마코브 행렬은 유저들의 이동변화율을 기술하였다. 그러나 이 변화율은 실제 세계에서는 시시각각 변하고 또한 변수도 너무나도 많다. 따라서 좀 더 정교한 예측을 하기 위해서는 마코브 행렬의 값을 스텝을 거칠 때마다 업데이트를 해야할 것이다. 하지만 마코브 행렬의 지속적인 업데이트가 필요없는 분야, 예를 들면 전자회로에서의 응용 문제 등에는 비교적 잘 들어맞을 때가 많다. 

이 외에도 마코브 체인의 단골 손님으로 나오는 문제인 날씨 예측, 인구이동 예측 문제 등이 있으니 직접 문제를 만들어 풀어보기 바란다. 문제 설정 및 풀이를 하는 과정은 애플-삼성 휴대폰 시장 점유율 예측 문제와 다르지않다. 아래는 소스코드이다. 

4. 마치며

이번 강의에서는 마코브 체인과 마코브 행렬에 대해 공부하였다. 마코브 행렬은 기본적으로 확률을 기반으로 한 상태 전이 행렬이며, 어떤 상태의 변환 이후에도 상태의 총량은 유지된다는 특성이 있다. 예를 들면 인구이동 문제의 경우 이동에 따른 상태 값의 변화(인구수)는 있겠으나 그 총량은 유지된다. 이는 확률을 기반으로 하며, 각 column의(혹은 row) 합은 1이 된다는 마코브 행렬의 특성에서도 유추해 볼 수 있다. 그 특성의 결과로 항상 하나의 고유값은 1을 가진다는 특성도 공부하였다. 또한 계차방정식과 연결지어 행렬의 거듭제곱 문제를 고유값의 거듭제곱 문제로 바꾸어 보다 빠르고 효율적인 계산방법을 배웠다. 

마지막으로 부가적인 설명을 하자면 애플과 삼성의 핸드폰 시장 점유율 예측에 대한 문제를 마코브 행렬로 모델링하고 이를 계산하는 과정에서 알아챈 분들이 있을지 모르겠으나, 이번 강의에서 설명한 마코브 체인은 1차 마코브 체인(first order Markov chain)이다. 이것이 의미하는 것은 상태 전이 과정에서 다른 시점에서의 상태와는 무관하고 오직 바로 이전상태에만 영향을 받는다는 것이다. 예측을 함에 있어 그 전전 상태까지 고려하는 2차 마코브 체인 등도 있지만 1차가 가장 보편적으로 사용된다. 또한 일반적으로 마코브 체인이라 하면 1차 마코브 체인을 의미한다. 

앞선 두 강의 Lecture 23-(1)과 Lecture 23-(2)에 이은 미분방정식과 선형대수의 마지막 강의이다. 지난 강의에서 우리는 u'=Au의 형태로 식이 주어졌을 때 그 풀이방법에 대해 주로 공부하였다. 이번에는 N차의 미분방정식이 주어졌을 때 이를 행렬의 형태로 만들어 푸는 방법에 대해 공부해보자. 앞의 두 강의를 충분히 이해했다면 이번에는 어렵지 않게 이해할 수 있을 것이다. 

6. 선형대수를 이용한 N차 미분방정식의 풀이

- 3rd order linear differential equation

미분방정식의 차수가 올라갈 수록 그 해를 구하기가 까다로워진다. 이때 선형대수를 이용하면 고차 미분방정식을 비교적 쉽게 풀 수 있다. 아래의 미분방정식을 보자. 

식 (38)은 3차 미분방정식을 나타낸다. y는 t의 함수이고, y의 최고차항은 3차항이며 Derivative form으로 표현한 식이다. 이제 이 미분방정식을 선형대수의 행렬을 이용하여 일반해(general solution)를 구해보도록 하자. 

해를 구하기 위해선 우선 식 (38)을 u'=Au의 형태로 만들어야 한다. 즉 u와 행렬 A를 식 (38)로 부터 만들어내야 한다. 어떻게 만들 수 있을까? 

가장 먼저 해야할 일은 벡터 u를 만드는 것이다. u에 정의되어야 할 것은 식 (38)에서 계수(coefficient)들을 제외한 나머지 미분항들을 넣으면 된다. u'=Au를 보면 u앞에 행렬 A를 곱해주면 u가 미분이 된다. 따라서 이 미분을 고려해보면 u에는 최고차항이 들어가면 안된다는 것을 알 수 있다. 그러므로 u에 들어갈 것은 최고차항 바로 전의 항까지 들어가면 되고, u'(u prime)에는 최고차항부터 최저차항 바로 전까지 들어가면 된다. 정리하면 아래와 같다. 

식 (40)과 같이 미분방정식을 u'=Au의 꼴로 정리하였고, u와 u'은 정의를 하였다. 이제 남은 과정은 행렬 A를 채워넣는 일이다. A를 어떻게 정리할 수 있을까? 힌트는 바로 미분방정식의 계수(coefficient)로부터 찾을 수 있다. A를 row별로 채워 넣는다고 생각하고 A의 row1과 u(t)가 곱해졌을 때를 가정해보면 아래와 같이 쓸 수 있다. 

식 (41.1)은 식 (40)의 u'(t) = row1 x u(t)를 가정한 것이다. 이때 ?로 표시된 계수들을 찾아야 하는데, 이 계수들은 사실 식 (38)에서 y'''을 제외한 나머지 항들을 우변으로 넘겨서 y'''=의 꼴로 정리를 했을 때와 같다. 이렇게 정리를 했을 때 식 (41.2)와 같이 쓸 수 있고, 이때의 계수들을 행렬 A에 삽입해주면 된다. 

row2부터는 굉장히 간단해진다. y''=?y''+?y'+?y에서 ?인 계수들을 찾으면 되는데, 사실 y''항에만 1을 곱해주고 나머지 항에는 0을 곱해주면 간단히 만들어지기 때문이다. 따라서 식을 다시 써보면 y''=1y''+0y'+0y와 같이 정리가 되고 row2는 [1 0 0]이 된다. 마찬가지 방법으로 row3까지 정리하면 row3=[0 1 0]이 될 것이다. 식 (40)을 행렬을 완성하여 다시 정리하면 아래와 같이 만들 수 있다. 

이렇게 하여 우리는 식 (38)의 3차 상미분방정식(3rd order ordinary differential equation)을 u(t)에 대한 1차 미분방정식으로 만들었다. u(t)가 벡터라고 해도 원래 미분방정식에 관련된 모든 정보를 담고있기 때문에 u(t)의 해를 구하면 y에 대한 해도 구할 수 있다. 이제 남은 것은 지금까지 배웠던 것처럼 식 (42)를 푸는 것이고, 그 핵심엔 A의 고유값(eigenvalue)과 고유벡터(eigenvector)가 있다. 이제 Lecture 21-(1)과 Lecture 21-(2)를 참고하여 행렬 A의 고유값과 고유벡터를 구해보자. 

고유값은 det(A-λI)의 계산을 통해서 람다에 대한 방정식을 만들어 풀 수 있다. (43)과 같이 행렬식(determinant)을 풀어서 고유값을 계산하면 각각 λ1=4, λ2=-1, λ3=1로 계산할 수 있다. 람다의 순서는 정하기 나름이므로 중요하진 않다. 이제 이 고유값을 가지고 고유벡터를 구해보자. 고유벡터는 알다시피 (A-λI)x=0에 대한 null space를 구하면 된다. 

식 (43)에서 고유값을 구하기 위해 세웠던 행렬식에 각각의 람다값을 대입하여 각 고유값에 해당하는 고유벡터들을 구하였다. 고유벡터를 구하는 과정에서 null space를 계산하야 하는데, 가우스 소거(Gauss elimination)를 이용하여 계산하는 방법도 있지만, 식이 간단할 경우 직관적으로 값을 대입해보면 쉽게 계산할 수 있다. 식 (44.1), (44.2), (44.3)은 각 고유값에 해당하는 고유벡터를 구하기 위한 null space 식을 나타내며, 소거법을 이용하지 않아도 Lecture 1-(1)에서 배웠던 column picture의 방식으로 생각하면 각 고유벡터의 원소값을 쉽게 구할 수 있다. 행렬 A의 고유값과 고유벡터를 구했으니 이를 이용하여 식 (42)의 일반해(general solution)를 아래와 같이 정의할 수 있다. 

식 (45.1)은 u(t)에 대한 일반해의 형태를 나타낸다.  C1, C2, C3는 임의의 계수이고 지수함수(exponential function)의 지수부에는 고유값이, 그리곡 각 항에는 고유벡터 x가 곱해진다. (45.2)는 (45.1)의 기본형태에 맞게 고유값과 고유벡터를 대입하여 일반해를 정리한 것이다. 만약 y에 대한 일반해를 정의하고자 한다면 고유벡터의 마지막 세 번째 원소만 고려하여 정리하면 된다. 이는 u(t)의 마지막 원소(component)가 y이기 때문이다. (45.2)의 각 고유벡터들의 마지막 원소값은 모두 1이기 때문에 y(t)에 대해서 해를 정리하면 (45.3)과 같이 된다. 

이렇게하여 미분방정식 (38)에 대한 일반해를 (45)와 같이 도출하였다. 아직 미지수로 남아있는 계수 C1, C2, C3는 초기 조건(initial condition)이 주어져 있다면 Lecture 23-(1)에서와 같이 구할 수 있다. 그러나 앞서 배웠듯이 이렇게 구한 일반해 말고도 또 하나의 방법이 있다. 바로 행렬지수함수(Matrix exponential) exp(At)를 이용하여 미분방정식의 해를 구하는 것이다. 지난 강의 Lecture 23-(2)에서 배운대로 행렬지수함수를 이용한 해를 구해보자. 

식 (46.1)과 (46.2)는 행렬지수함수로 정의하는 해 u(t)를 나타낸다. (46.3)의 S는 고유벡터행렬(eigenvector matrix)로써 각 고유벡터들이 column이 되는 행렬이다. (46.4)는 고유값 행렬(eigenvalue matrix)의 행렬지수함수이고 고유값 행렬이 대각 원소만 존재하기 때문에 행렬지수함수도 대각 원소만 존재한다. (46.5)는 S의 역행렬(inverse matrix)이고 S의 행렬식(determinant)과 cofactor matrix의 전치(transpose)를 이용하여 계산할 수 있다. 자세한 사항은 Lecture 20-(1)를 참고하기 바란다. 

식 (46.6), (46.7)은 행렬지수함수의 해의 정의에 따라 구한 해 u(t)를 나타낸다. 식이 약간 복잡해보이지만 계산 과정을 나타내기 위해서 자세히 풀어서 나타내었다. 실제 계산은 컴퓨터를 이용하여 하는 경우가 많기 때문에 이렇게 계산하는구나 정도만 알면 된다. 결과적으로 미분방정식의 해를 행렬지수함수를 이용하여 나타내었고 이는 고유벡터행렬 S와 고유값행렬의 행렬지수함수 e(λt), 그리고 고유벡터행렬 S의 역함수의 곱으로 정의됨을 다시 확인하였다. u(t)의 최종 결과값을 확인하려면 위의 식의 우측에 초기값 u(0)를 곱해주고 알고자하는 시점의 값 t를 대입하면 위의 미분방정식에 대한 최종 값을 구할 수 있다. 

- N-th order linear differential equation with linear algebra

앞서 우리는 3차 미분방정식을 행렬을 이용하여 1차의 미분식으로 만들어서 그 해를 구하는 방법을 배웠다. 그렇다면 4차, 5차, 나아가 N차 미분방정식일 때도 가능한 이야기일까? 그렇다. N차 미분방정식도 행렬을 이용하여 1차의 미분식으로 연결시켜 해를 구할 수 있다. 이에 대한 설명을 위해 앞서 다루었던 식을 일반화시켜 정리해보자. 

식 (47.1)은 3차 미분방정식이고 (47.2)는 이를 1차의 미분식으로 만든것이다. 여기서 나름의 규칙을 볼 수 있는데, 행렬 A의 row1은 미분방정식의 두 번째 고차항부터 그 아래 항의 계수들, 즉 (47.1)에선 a, b, c들이 부호가 바뀌어서 채워진다는 것을 알 수 있다. 그리고 뒤이어 row2, row3부터는 가장 왼쪽 끝 원소부터 대각선 방향으로 차례로 1이 채워지고 나머지 부분은 0이 된다. 이 규칙을 N차 방정식으로 정리해보면 아래와 같다. 

N차 미분방정식을 행렬의 형태로 만들면 식 (48)과 같이 nxn크기의 행렬이 만들어진다. 이때 첫 번째 row는 최고차항을 제외한 나머지 항들의 계수가 부호가 바뀌어 채워지고, 나머지 row들은 왼쪽부터 1이 대각선 방향으로 차례로 채워지게 된다. 결국 N차 미분방정식도 위와 같이 1차의 nxn크기의 행렬로 이루어진 식으로 바꾸어 풀 수 있다. 

7. 마치며

총 세 편에 걸쳐 미분방정식과 선형대수를 배웠다. 고차의 미분방정식이라도 선형대수를 이용하면 1차의 식으로 바꿀 수 있다는 것은 정말 매력적이지 않을 수 없다. 또한 해를 구하기 위한 핵심에는 고유값과 고유벡터가 빠지지 않고 등장한다. 이 고유값/고유벡터를 통해 행렬지수함수를 정의할 수 있었으며, 그 배경에는 테일러 급수라는 훌륭한 근사법이 존재했다. 이쯤 되면 고유값과 고유벡터의 깊이는 도대체 어디까지인가 하는 생각이 들기도 한다. 아무쪼록 세 번에 걸친 이번 강의를 통해 미분방정식과 선형대수의 연결고리를 이해하는데 도움이 되길 바랍니다. 

이번 포스팅은 지난 강의 Lecture 23-(1)에 이은 미분방정식과 선형대수의 두 번째 이야기이다. 앞선 강의와 이어지는 내용이니만큼 앞 강의 내용을 잘 이해하고 보는 것을 추천한다. 

4. 행렬 지수함수(Matrix exponential)

- Diagonalization for differential equation

지난 강의에서 계속 다루어왔던 아래의 식을 다시 생각해보자. 

식 (18)의 행렬 A는 벡터 u와 그의 미분값 u'을 연결시키는(couple)역할을 한다. 그런데 우리는 앞선 강의에서 A의 고유값과 고유벡터를 구하여 지수함수(exponential function)로 이루어진 일반해(general solution)를 도출한바 있다. 일반해는 각각의 상수와 고유값, 그리고 고유벡터들로 이루어진 식이며, 선형조합(Linear combination)의 형태로 분리가 되어 있다. 이는 기존의 행렬 A가 고유값/고유벡터에 의해서 분리(uncouple)가 되는 것이다. 이것이 무슨 말인지 이해하기 위해 식 (18)을 고유벡터의 관점에서 다시 정리해보자. 

식 (19)는 (18)의 u를 고유벡터와 가중치(weight)곱의 선형조합에 대한 식으로 정리한 것이다. 이를 행렬에 대한 식으로 정의한 것이 (19.1)이다. 여기서 v는 고유값에 대한 지수함수를, S는 고유벡터들로 이루어진 고유벡터행렬을 의미한다. u=Sv를 풀어서 써보면 식 (19.2)와 같이 일반해로 정리되기 때문에 식에 무리는 없다. 이제 식 (19)의 내용을 (18)에 대입하여 A를 분해하여 정리해보자. 

(20.1)은 (18)에서 u대신 Sv를 대입하여 정리한 것이다. 이때 S는 상수이기 때문에 좌변의 미분식에서 S는 앞으로 뺄 수 있다. 양변의 좌측에 S의 역함수를 곱해주면 (20.2)와 같이 정리할 수 있는데, 우변의 S^-1AS가 어디서 많이 본 형태이다. 바로 Lecture 22의 행렬의 대각화(diagonalization)에서 배웠던 형태인데, S^-1AS는 Λ(capital lambda)와 같음을 배웠다. 이를 통해 (20.2)를 다시 (20.3)과 같이 정리할 수 있다. 

식 (20.3)을 보면 dv/dt=Λv가 되는데, 식 자체의 의미를 놓고 보면 벡터 변수 v를 시간으로 미분한 것은(시간의 변화에 따른 변수 변화량의 비율) 변수에 어떤 상수값을 곱한 것과 같다이다. 이 의미가 성립하려면 v는 식 (19.1)과 같이 지수함수로 정의되어야 한다. 지수함수는 미분한 값이 바로 자기 자신이 되기 때문이다(지수함수의 지수부의 상수가 미분하면 지수함수 앞에 곱해지고 이것이 람다가 되는 꼴). v는 벡터이기 때문에 여러 단일 변수를 동시에 정의할 수 있으며 아래 식은 이에 대한 내용을 정리한 것이다. 

식 (21)은 dv/dt=Λv를 2x2크기로 가정하고 식을 풀어 정리한 것이다. 실제로 위의 식과 같이 v에 지수함수를 대입하여 t로 미분해보면 식이 성립함을 알 수 있다. 여기서 중요한 것은 위의 식을 활용하여 미분방정식을 푸는것이다. 우리가 결국 구하고자하는 것은 u이고, v는 식 (19)와 같이 u=Sv의 꼴로 나타냈을 때 S의 각 고유벡터들에 곱해지는 계수(coefficient)에 해당한다. 그러므로 먼저 v의 시간 t에 대한 함수를 정의해보자. 지난 강의 Lecture 23-(1) 에서 배웠듯이 미분방정식의 해(solution)는 지수함수(exponential function)를 도입하여 정의할 수 있다. 단일 변수 대신 다변수에 대해 정의한 행렬에 지수함수를 도입한다는 것이 다를 뿐 기본적인 개념은 같다. 

(22.1)은 식 (21)에서 도출한 식이고, 여기서 양변에 dt를 곱해주고 1/v를 곱해준 뒤 적분(integral)을 해준 것이 (22.2)이다. 적분을 하면 (22.3)이 되고, 좌변의 C1을 우변으로 이항하고 자연로그를 없애주기 위해 양변에 지수함수(exponential function)를 취하면 (22.4)가 된다. 절대값은 +, 혹은 -값이 되지만 임의의 값이 될 수 있으므로 +값을 선택하여 e^(c2-c1)자체를 상수 C로 치환하여 (22.5)를 정의할 수 있다. 

여기서 우리는 식을 미분방정식의 해를 구하기 위한 함수의 형태로 정리해야 한다. 즉 특정 시간 t에서 해당 미분방정식의의 값이 무엇인지 구할 수 있도록 v의 t에 대한 함수 v(t)의 꼴로 정의해야 한다. v(t)의 꼴로 정의하는 것은 반드시 초기값이 있어야 한다. v에서 t=0인 경우는 지수함수부분이 1이 되므로 상수 c1, c2만 남게 된다. 이 상수 벡터가 v의 초기값이 되고 (22.5)의 C를 대체하여 생각할 수 있다. C를 v의 초기값으로 대체하여 (22.6)과 같이 v(t)를 정의하였다. 

v(t)를 구했으니 이제 u(t)를 구해보자. 어떻게 구할 수 있을까? 바로 식 (19)의 u=Sv로부터 구할 수 있다. 아래의 식을 보자. 

u=Sv를 u(t)의 꼴로 정리하면 시간의 변수 t를 포함하고 있는 v도 역시 v(t)의 형태로 변해서 (23.2)와 같이 정리할 수 있다. 그 다음 (23.2)의 v(t)를 (23.1)의 v(t)의 우변으로 식을 치환하면 (23.3)과 같이 정리가 된다. u(t)도 마찬가지로 초기값 u(0)가 포함되어야 하는데, u(0)=Sv(0)이고 이로부터 (23.4)와 같이 v(0)를 정의할 수 있다. 여기서 정리한 v(0)를 (23.3)에 대입하여 정리하면 최종적으로 (23.5)와 같이 정리할 수 있다. 

정리하자면 우리가 최초에 풀고자 했던 문제는 식 (18)에서 정의한 것과 같이 다변수(multi-variable) 함수에 대한 미분방정식을 푸는 것이었다. 여기서 A는 u와 u의 미분값 사이를 정의하며, 이들을 연결해주는(couple) 역할을 한다. 특정 시간 t에서의 미분방정식의 해(solution)를 구하기 위해선 u를 시간 t에 대한 함수 u(t)의 꼴로 정의를 해줘야 하는데, 이때 사용한 방법이 이 행렬 A의 고유값과 고유벡터를 이용하여 정리하는 것이다. A의 고유값과 고유벡터를 이용하여 대각화를 하고, 여기에 지수함수를 적용하여 식 (23.5)와 같이 미분방정식의 해를 구할 수 있는 식을 만든 것이다. 결국 A의 고유값과 고유벡터를 이용하여 최초의 식 A를 분해하고 식 (19)와 같이 지수함수를 도입하여 다변수 미분방정식의 해(solution)를 구하기 위한 식을 (23.5)와 같이 도출하였다. 

그런데 식 (23.5)의 가운데 부분은 다시 아래와같이 정리할 수 있다. 

u(t)의 가운데 대각화된 식을 e^At라고 할 수 있는데, 이와 같이 지수함수를 밑으로 하는 형태의 행렬을 행렬지수함수(matrix exponential)라 한다. 뒤이어 알아볼 내용은 식 (24.1)의 행렬지수함수이다. 

행렬지수함수를 공부하기 전에 먼저 알아야 할 내용이 있다. 바로 테일러 급수(Taylor series)이다. 인터넷에 테일러 급수에 대한 수많은 자료가 있지만, 매우 중요한 개념이니만큼 간략히 정리하고 넘어가도록 하자. 

- Taylor series

테일러 급수(Taylor series)란 어떤 함수를 특정 위치 x=a에서 근사(approximation)하는 방법이다. 테일러 급수의 정의는 아래의 식과 같다. 

식 (25)를 풀어서 설명해보자면 n번 미분이 가능한 어떤 함수 f(x)가 있을 때, 이 함수의 특정 위치에서의 값, 즉 x=a에서의 값은 f(x)와 그의 미분값들로 계산한 항들의 무한대의 합으로 표현할 수 있다는 것이다. 이때 각 항은 미분차수값의 factorial로 나누어주어야 한다. 즉 1차 미분이면 1 factorial(1!), 2차 미분이면 2 factorial(2!)로 나누어줘야 한다. 테일러 급수의 식이 왜 이와 같이 표현되었는지를 이해하려면 멱급수(power series)의 전개부터 이해를 해야하는데, 이 부분은 본 포스팅의 범위를 넘어가기 때문에 다음의 동영상을 보는 것을 추천한다(Power series/Euler's great formula). 

위에서 설명한 테일러 급수에서 한 가지 중요한 것은  테일러 급수가 근사하는 것은 함수 전체가 아니라 함수의 특정 위치(x=a)에서의 근사만을 의미한다. 물론 어떤 함수의 경우엔 근사하는 위치 주변에서의 값도 잘 맞아떨어질 순 있겠지만, 기본적으로는 x=a에서의 근사를 하기 위함이다. 이와 같은 설명만으론 이해하기 어려울 수 있으므로 예를 한 번 들어보자. 

식 (26)의 함수 f(x)는 3차 함수이며 3번 미분이 가능하기 때문에 3차 미분까지만 고려하면 된다. 그 이상은 0이 되기 때문에 의미가 없다. 위의 함수를 x=0일때와 x=-2일때 각각 근사를 해보자. 참고로 x=0일때의 테일러 급수를 매클로린 급수(Maclaurin series)라 한다. 

식 (27)은 (26)의 함수 f(x)를 x=0에서 근사한 것이다. (27.1)은 f(x=0)일때의 값과 3차 미분까지의 값들을 나열하였다. (27.2)는 원래의 함수 f(x)를 테일러 급수로 근사한 결과이다. 각 미분항에 곱해지는 (x-a)^n에서 x=0에서의 근사이기 때문에 a에는 0이 대입된다. 정리한 결과 원래의 함수 f(x)가 그대로 복원되고 0을 대입하면 마찬가지로 2가 나오는 것을 볼 수 있다. 

이제 x=-2일 때의 근사(approximation)를 수행해보자. 

식 (28)은 x=-2에서 테일러 급수를 통해 f(x)를 근사한 모습이다. 마찬가지로 원래의 함수를 잘 복원하는 것을 볼 수 있다. 

그런데 여기서 어떤 사람들은 이렇게 생각할 수 있다. "이렇게 그대로 다시 복원할거면 도대체 테일러 급수를 왜 사용하는거지? 그냥 그대로 사용하면 될텐데..". 사실 위에서 예로든 함수 f(x)는 다항식의 유한번의 사칙연산으로 구할 수 있는 대수함수(algebraic function)이고, 따라서 테일러 급수를 통해 원래의 함수를 그대로 복원할 수 있다. 테일러 급수가 진정 빛을 발하는 경우는 대수적 연산만으로는 정확한 값을 구할 수 없고 극한의 과정을 거쳐야만 정확한 값을 구할 수 있는 함수, 이를 테면 지수 함수(exponential function), sin, cos, log등 초월 함수(transcendental function)를 구할 때이다. 초월 함수의 대표적인 테일러 급수의 예로 지수함수(exponential series)와 기하급수(geometric series)에 대한 식을 살펴보자. 

식 (29.1)과 (29.2)는 각각 지수함수와 기하 급수를 나타내며, (29.2)의 기하급수는 x의 크기가 1보다 작을 때에만 좌변과 우변의 식이 성립(수렴)한다. 이 외에도 sin, cos, log등의 다양한 초월 함수들이 존재하며 자세한 사항은 테일러 급수 위키(Taylor series)를 참고하기 바란다. 이번 강의를 위해선 일단 위의 두 함수에 대한 테일러 전개(Taylor expansion)를 알면 된다. 참고로 위의 지수함수와 기하급수 두 개의 테일러 급수는 각각 x=0일 때의 매클로린 급수이다. 

그렇다면 테일러 급수가 이러한 초월함수에서 힘을 발휘하는 이유는 무엇일까? 그것은 바로 근사(approximation)하는 능력에 있다. 제아무리 성능 좋은 수퍼컴퓨터라해도 무한번의 연산을 수행할 수는 없을 것이다. 따라서 우리는 해당 함수의 값을 원하는 정도의 정밀도가 나올 때까지만 근사하여 구하면 된다. 이때 이러한 근사의 과정을 테일러 급수를 통해서 수행하게 된다. 실제로 공학용 계산기에서 sin함수의 값은 테일러 급수를 통해 어느 정도로 근사된 값이 계산되어 나오는 것이다. 테일러 급수라는 개념이 없었다면 이러한 초월함수들의 계산 과정에 큰 어려움이 있었을 것이다. 근사(approximation)에 대한 이해를 돕기 위해 (29.2)의 기하 급수의 수렴에 대한 예를 살펴보자. 

식 (30.1)은 x=1/2일 때의 기하급수의 계산 결과를 나타낸다. 기하급수의 특성상 |x|<1일 때 더해가는 항, 즉 n을 늘려나갈수록 특정 값에 수렴하게 되는데, (30.2)는 두 번째 항까지, (30.3), (30.4)는 각각 세 번째, 네 번째 항까지 더한 결과를 나타낸다. 항이 늘어날 수록 수렴값 2에 가깝게 되는 것을 볼 수 있다. 사실 무한대로 더해야지만 2에 완전히 수렴하겠지만, 무한대로 더하는 것은 불가능하기 때문에 만족할만한 수준의 정밀도를 설정해놓고 그에 상응하는 항까지 더하여 근사를 하면 된다. 

이 외에도 테일러 급수는 어떤 복잡한 모델을 단순화 시켜서 풀거나, 미분방정식의 해를 구할 때, 정적분을 구할 때 등 과학과 공학에 있어서 다양하게 활용되는 굉장히 중요한 개념이므로 잘 이해하도록 하자. 

- Matrix exponential

앞서 테일러 급수를 간략히 공부했으니 이제 행렬지수함수(Matrix exponential)를 공부해보도록 하자. 미리 말하지만 행렬지수함수를 공부하는 이유는 식 (24)에서와 같이 미분방정식의 해(solution)에 대한 식을 도출하는 과정에서 행렬지수함수가 나오기 때문에 이를 정의하는 법을 알기 위함이다. 우리가 알아볼 행렬지수함수의 형태는 아래의 식과 같다. 

우리가 학창시절에 공부하면서 봤던 지수함수의 형태는 대부분 지수함수 e의 지수부에 상수, 혹은 변수가 있는 형태였다. 그런데 식 (31)에 정의된 행렬지수함수는 지수함수의 지수부에 행렬(matrix)이 있는 경우이다. 여기서 행렬 A는 고정된 값이고 t가 기존의 x와 같은 변수이다. 이러한 행렬지수함수는 어떻게 정의할 수 있을까? 그 힌트는 바로 앞서 공부했던 테일러 급수에 있다. 식 (29.1)의 지수함수의 테일러 전개가 (31)의 행렬지수함수에도 그대로 적용되는 것이다. 다만 단일 변수 대신 행렬의 형태로 전개가 되는 것이다. 아래의 식을 살펴보자.  

비교해보기 쉽도록 원래의 식 (29.1)과 함께 정리하였다. (29.1)의 1이 (32)에선 단위 행렬 I가 되고, x는 At가 되며 나머지 항들도 같은 패턴을 보이는 것을 알 수 있다. 이때 (29.1)의 1+x 부분이나 (32)의 I + At등의 부분은 사실은 n^0 !, n^1 !로 각각 나누어 준 것이다. 지수함수와 같은 패턴으로 행렬지수함수(Matrix exponential)이 계산된다는건 알겠는데, 정확히 어떤 형태로 나온다는건지 감이 오지 않을 수 있다. 아래의 행렬을 이용하여 행렬지수함수를 계산해 보도록 하자. 

식 (33.1)의 행렬 A와 t=1임을 가정하여 행렬지수함수를 계산하였다. 계산 항은 총 10차 까지 계산하였다. 계산 결과값은 (33.3)과 같은 결과를 보인다. 이 값이 과연 맞는 값인지 MATLAB을 통해 확인해보자. MATLAB에서 행렬지수함수를 계산하기 위해선 expm()함수를 이용하면 된다. 아래 그림은 계산 결과와 소스코드이다. 

Fig. 1 Matrix Exponential계산 결과와 소스코드. MATLAB 내장함수결과(ans =), 테일러급수 n=10차 계산결과(M =)

Fig. 1의 좌측은 MATLAB을 이용한 계산결과를, 우측은 소스코드를 나타낸다. 좌측 결과에서 ans = 는 MATLAB의 내장함수인 expm()함수를 이용한 계산결과이고 M = 은 테일러급수를 이용하여 10차까지 근사한 결과이다. 근사값이 내장함수값과 거의 일치함을 보이며, 차수를 더해갈 수록 오차는 줄어들게 되고 보다 정확한 근사를 할 수 있다. 이와 같이 행렬지수함수가 테일러 급수로 표현될 수 있음을 증명하였다. 

위에서 행렬지수함수를 테일러 급수를 통해 정의하는 법을 배웠는데, 우리는 앞서 지수함수와 더불어 식 (29.2)의 기하급수(geometric series)에 대해서도 배웠다. 기하급수 역시 식 (29.2)의 기존의 형태에 행렬을 대입하여 정리할 수 있다. 아래 식을 보자. 

비교를 위해 역시 기존의 식 (29.2)와 함께 정리했다. (29.2)의 1/1-x를 다시 쓰면 (1-x)^-1와 같이 쓸 수 있는데, 식 (34)와 같이 행렬에 대해 정리했을 때 이는 역함수로 해석할 수 있다. 즉 (34)의 (I-At)의 역함수로 생각할 수 있다. 1/(I-At)가 아니다. (29.2)에서 x의 크기가 1보다 작은 경우에만 식이 성립한다고 했다. 마찬가지로 행렬의 형태로 정리한 (34)에서도 이와 같은 조건이 있는데 바로 At의 모든 고유값이 1보다 작은 경우에만 식이 성립한다는 것이다. 이때 중요한것은 행렬 A에 대한 고유값이 아니라 상수 t가 곱해진 At의 모든 고유값이 1보다 작아야 한다는 것이다. 실제 값을 대입하여 계산 결과를 살펴보자. 

위의 계산 결과가 실제로 (I-At)의 역함수를 계산했을 때와 같은지 MATLAB을 통해 확인해보자. 

Fig. 2  기하급수(Geometric series)의 Matrix Exponential계산 결과와 소스코드. MATLAB 계산결과(ans =), 기하급수 n=100차 계산결과(M =)

Fig. 2의 왼쪽은 MATLAB 계산결과이고 오른쪽은 소스코드이다. V와 D는 각각 At의 고유벡터와 고유값을, mag1과 mag2는 At의 λ1과 λ2의 크기를 각각 의미한다. 

다음으로 ans = 는 (I-At)의 역함수를 계산한 결과이며, M은 t=2.2일 때 기하급수를 테일러 급수를 이용하여 계산한 결과이다. 이때 차수는 n=100차 까지 계산하였다. 결과적으로 역함수를 계산한 결과와 급수를 통해 근사한 결과인 M을 비교해보면 거의 같음을 알 수 있다. 물론 차수를 높일 수록 값은 더욱 같아진다. 이로써 기하급수를 행렬지수함수의 형태로 정리한 식 (34)가 조건이 맞을 경우 식이 성립함을 증명하였다. 위의 MATLAB코드에서 t와 A의 값을 바꾸어가며 테스트 해보면 mag1과 mag2가 1이 넘는 경우엔 식이 성립하지 않음을 볼 수 있으니 테스트해보기 바란다. 

급수(series)의 관점으로 봤을 땐 식 (32)가 (34)보다는 훨씬 뛰어난 특성을 보인다. (32)의 경우엔 A나 t가 얼마나 크던 간에 항을 더해갈수록 큰 수로 나누어주기 때문에 나중으로 갈 수록 0에 가까워지고 결국 수렴하게 된다. 반면에 (34)의 경우엔 At의 고유값이 1보다 큰 경우 발산할 수 있다. 기하급수의 행렬지수함수를 알아봤지만, 이번 강의에서 우리가 눈여겨봐야할 급수는 바로 식 (32)의 행렬지수함수이다. 이 행렬지수함수를 이용하여 앞에서 다루었던 식 (24)의 행렬지수함수를 증명해보도록 하자. 그 사이에 많은 내용을 다루었기 때문에 편의를 위해서 식 (24)를 다시 써보도록 하자. 

우리가 알고자 하는 것은 식 (24.2)가 어떻게 성립하는지, 즉 $Se^{\Lambda t}S^{-1}$를 어떻게 $e^{At}$로 정의할 수 있는지이다. 식 (32)의 행렬지수함수를 이용하여 이를 정의할 수 있다. 아래의 식을 보자. 

식 (36.2)는 (36.1)의 행렬지수함수에서 A를 $S \Lambda S^{-1}$로 치환한 것이다. $A=S \Lambda S^{-1}$는 Lecture 22의 행렬의 대각화 부분에서 배웠던 내용이다. (36.2)의 첫 번째 항은 t의 0제곱과 zero factorial이기 때문에 지수함수 e가 사라져서 고유벡터행렬 S만 남은 것이다. (36.2)에서 고유벡터행렬 $S$와 $S^{-1}$는 상수이기 때문에 양쪽으로 뺄 수 있다. 이렇게 S들을 빼서 정리한것이 (36.3)이고, 괄호안에는 고유값 행렬(Λ)에 대한 식이 남는데, 이는 $e^{\Lambda t}$로 볼 수 있다. 따라서 (36.4)와 같이 정의할 수 있고 식 (24.1)을 증명한 것이다. 

한 가지 기억해야 할 것은 식 (36.1)은 언제나 성립하는 반면 식 (36.4)는 항상 성립하는 것은 아니라는 것이다. (36.4)가 성립하기위해선 A가 대각화가 가능해야 한다는 것이다. 다시말하면 행렬 A는 n개의 독립(independent)인 고유벡터(eigenvector)를 가져야 한다. 이는 행렬의 대각화를 배울 때 강조했던 내용이며 (36.4)가 성립하기 위해선 반드시 충족시켜야 할 조건이다. 중요한 내용이므로 다시 한 번 강조하니 꼭 기억하도록 하자. 

이즈음에서 다시 한 번 생각해보자. 우리가 이렇게 행렬지수함수를 계산하는 법을 공부하는 이유는 무엇인가? 바로 선형연립미분방정식(simultaneous linear differential equation)의 해(solution)를 구하기 위함이며, 그 해는 식 (24)와 같이 u(t)=의 형태이다. u(t)의 식에서 결국 구해야 할 것은 주어진 행렬 A를 기반으로 식 (24.1)의 행렬지수함수를 구하는 것이다. 사실 (36.1)과 같이 테일러급수의 형태로 근사하여 해를 구할 수도 있다. 그러나 계산과정에 있어 비용이 많이 드는 것이 사실이다. 따라서 고유값과 고유벡터를 기반으로 한 (36.4)와 같은 형태로 해를 구하는 것이 훨씬 효율적이다(n개의 독립인 고유벡터가 존재한다는 가정하에서). 

행렬 A의 고유값과 고유벡터를 구하는 것은 어렵지 않지만, 문제는고유값 행렬(Λ)에 대한 행렬지수함수이다. 결국 람다에 대한 행렬지수함수를 얻기 위해 (36.1)의 과정을 거쳐야 한다면 (36.4)의 형태가 큰 장점은 없어보인다. 하지만 고유값 행렬(Λ)은 대각행렬이고, 이에 대한 행렬지수함수는 그 계산방법이 훨씬 간단하다는 이점이 있다. 아래의 식을 보자. 

식 (37)은 고유값 행렬(Λ)의 행렬지수함수를 나타낸 것이다. 보다시피 고유값만 알면 쉽게 계산할 수 있는 형태이다. 과연 위의 형태대로 계산하는 것이 맞는지 MATLAB을 통해 확인해보자. 

Fig. 3 고유값 행렬(Λ)의 행렬지수함수 계산결과(Left)와 MATLAB소스코드(Right)

Fig. 3은 고유값 행렬의 행렬지수함수의 계산 검증 결과를 나타낸다. 왼쪽의 계산결과에서 V와 D는 각각 행렬 A의 고유벡터행렬과 고유값 행렬이고, ans는 MATLAB 내장함수인 expm()으로 고유값 행렬의 행렬지수함수를 계산한결과를, 그리고 expD는 식 (37)의 형태로 직접 계산한 결과이다. 비교해보면 정확히 일치하는 것을 볼 수 있다. 이로써 계산이 훨씬 간편해졌다. 물론 코드가 더 추가되긴 했지만 행렬이 커짐에 따라 이 계산방식이 빛을 발할 것이다. (내장함수인 expm()의 계산 알고리즘이 이미 이와 같은 방식을 사용한다거나, 혹은 더 좋은 방법을 사용한다면 이야기는 달리진다)

정리하자면 최초에 식 (18)과 같이 행렬의 형태로 선형연립미분방정식이 있고 초기값(initial condition)이 주어졌을 때 우리는 식 (24)와 같은 형태로 해를 구할 수 있다. 해를 구하기 위해선 행렬 A의 고유값과 고유벡터를 기반으로 식 (24.2)와 같이 해를 정의할 수 있으며, 이때 고유값 행렬의 행렬지수함수는 (37)의 형태로 계산할 수 있다. 

마지막으로 중요한 포인트 한 가지를 다시 짚어보자. 바로 지난 강좌에서 일부 다루었던 내용이지만 중요한 부분이므로 다시 언급하도록 하겠다.

일반적으로 어떤 시스템의 해를 구한다는 것은 시간(t)이 지날수록 그 값이 수렴하는 형태일 것이다. 그렇다면 언제 수렴하는 것일까? 그 해답은 식 (37)의 고유값 행렬에 있다. 식 (24.2)에서 u(t)는 고유벡터행렬 S와 고유값 행렬로 정의되어 있는데, 여기서 고유벡터행렬은 상수이기 때문에 고정되어 있다. 움직이는 것은 t가 곱해져서 영향을 받는 고유값 행렬의 행렬지수함수인데, 식 (37)을 보면 대각원소가 각 고유값에 대한 지수함수로 정의되어 있다. 이것이 의미하는 것은 t가 양의 방향으로 나아간다고 가정했을 경우 해가 수렴하기 위해선 행렬 A의 모든 고유값들이 0보다 작아야 한다는 것이다. 보다 더 정확히 얘기하자면 행렬 A의 모든 실수부(real part)의 고유값이 0보다 작아야지만 시간 t가 지날 수록 해가 안정상태(stable)가 된다. 이는 어떤 시스템 행렬 A를 설계하는데에 있어서 중요한 포인트이므로 잘 기억하도록 하자. 

- Test

지금까지 배운 행렬지수함수(Matrix exponential)를 이용한 미분방정식의 해를 구하는 방법이 과연 맞는지 확인해보도록 하자. 지난 강의 Lecture 23-(1)  에서 직접 구했던 해를 이번 강좌에서 배운 행렬지수함수를 이용하여 구해서 결과가 같은지 확인해보자. 아래는 MATLAB을 이용해 두 방식을 테스트 한 결과이다. 

Fig. 4 직접 구한 해와 행렬지수함수를 이용한 해의 비교

Fig. 4의 왼쪽이 계산 결과, 오른쪽이 소스코드이다. 왼쪽의 결과에서 A는 행렬, t는 임의의 시간 4.75를 가정하였고 초기값인 v0는 동일하게 [1 0]T으로 설정하였다. Sol=로 표시된 값이 행렬지수함수를 이용하여 계산한 값이고 u=로 표시된 값이 기존에 직접 계산했던 값의 결과이다. 계산 결과가 동일한 것을 볼 수 있다. 

5. 마치며(continue)

이번 강좌에선 지난 시간에 이어 미분방정식과 선형대수에 대한 내용을 다루었다. 어떤 시스템의 미분방정식을 세우고 그 해를 u(t)의 꼴로 구하는 과정을 배웠다. 핵심은 결국 시스템 행렬 A의 고유값과 고유벡터에 있으며, 이들을 기반으로한 행렬지수함수를 통해 그 해를 정의할 수 있다. 또한 행렬지수함수를 테일러 급수를 이용하여 정의하였으며, 행렬의 대각화(Diagonalization)를 통해 행렬지수함수를 정의하는 과정을 배웠다. 다음 포스팅에서 미분방정식과 선형대수의 마지막 강좌를 하도록 하겠습니다. 

이번에 포스팅할 내용은 미분방정식(Differential equation)을 선형대수를 이용하여 푸는 방법이다. 이번 포스팅의 내용을 이해하기 위해선 고유값/고유벡터 강의와 대각화(Diagonalization)내용을 먼저 학습하고 오는 것을 추천한다. 특히 지난 강의 Lecture 22에서 다루었던 대각화와 행렬의 거듭제곱, 그리고 차분방정식(Difference equation)과 관련이 깊기 때문에 이를 확실하게 이해한 뒤 이번 강좌를 보는 것이 좋다. 본격적인 내용을 학습하기전에 미분방정식이란 무엇이며, 이번 강의에 필요한 관련 내용들에 대해서 간략히 살펴보고 가도록 하자. 

1. 미분방정식의 개요(Introduction to Differential equation)

- What is a differential equation?

미분방정식(Differential equation)이란 어떤 방정식(equation)에 도함수, 즉 미분(derivative)이 포함된 것이다. 이는 물리학에서의 운동방정식, 가령 스프링에 매달린 추의 운동 분석에 대한 방정식 등을 세우는데에 사용되기도 하고 그밖에 경제학, 생물학 등 다양한 분야에 걸쳐 사용될 수 있는 수학적 방법이다. 미분방정식은 일반적인 방정식보다 식을 세우고 이해하거나 해(solution)를 구하기가 굉장히 까다롭다. 미분방정식에 대해서 본격적으로 알아보기 전에 다른 형태의 방정식들은 어떤 것들이 있는지 살펴보자. 아래는 미분방정식이 아닌 일반적인 형태의 방정식들을 나열해논 것이다. 

식 (1.1)부터 차례로 선형(일차) 방정식(Linear equation), 이차 방정식(Quadratic equation), 삼차 방정식(Cubic equation)을 각각 나타내며, 최고차항에 따라 구분되는 형태이다. 이들 방정식은 미지수(unknown)와 그들의 앞에 곱해져있는 상수인 계수(coefficient)로 구성되어 있으며, 덧셈, 뺄셈 연산을 통한 조합으로 구성되어 있다. 이와 같이 어떤 미지수의 거듭제곱과 계수들의 조합으로 구성되어있는 식을 대수방정식(algebraic equation)이라고 한다. 미지수의 최고차항에 따라 식을 구분할 때 차수가 일차인 식을 특히 선형방정식이라고 하는데, 이는 식 자체가 선형성(Linearity)을 가지기 때문이다. 2차식부터는 비선형성(non-linearity)을 가진다. 

사실 엄밀히 따지면 식 (1.1)은 대수방정식은 아닌데, 미지수가 x와 y 두 개로 구성되어있기 때문이다. (1.1)처럼 미지수가 한 개 이상으로 구성되어 있는 경우 이들을 다항방정식(polynomial equation)이라고 하며 종종 다항식(polynomial)이라고 쓰인다. 다항방정식은 대수방정식을 포괄하기 때문에 보통 대수방정식이라는 말 보다는 다항식이라는 용어를 사용한다. 

이제 미분방정식이 어떤 형태인지 살펴보자. 미분방정식의 형태는 식 (1)의 다항식(polynomial)형태의 식에 미분(derivative)의 개념이 더해진 식이기 때문에 훨씬 복잡해진다. 아래의 식을 보자. 

식 (2)는 하나의 미분방정식을 여러 가지 형태로 표현한 것이다. 식 (2.1)은 derivative form이고, y의 미분식을 prime을 이용하여 표현하였다. prime하나당 한 번의 미분을 나타내며, 이차미분은 y의 double prime으로 표현한다. 

식 (2.2)는 differential form이고 미분을 dy/dx(delta x, delta y로 읽음)로 표현한다. 여기서 분모에 존재하는 미지수 y는 종속 변수(dependent variable)이고, 분자인 x는 독립 변수(independent variable)라고 한다. 이는 (2.2)에서 y는 x에 대한 함수이기 때문이다. 즉 y는 x가 변함에 따라 그 변화량이 얼마나 되는지를 봐야 하는 x에 종속되는 변수이고, x는 어디에도 종속되지 않는 독립적인 변수이다. 예를 들어 어떤 공장의 온도에 따른 가스의 압력에 대한 관계를 미분방정식으로 만들었다고 생각해보자. 여기서 온도를 x, 가스의 압력을 y라고 했을 때 온도의 변화에 따른 가스의 압력을 dy/dx와 같은 미분식을 넣어서 만들 수 있을 것이다. x로 표기된 온도는 독립적으로 우리가 조정하는 것이고, y로 표시된 가스의 압력은 x를 임의로 변화시킴에 따라 변화된 결과를 관찰하는 종속변수라고 이해하면 된다. 

마지막 (2.3)은 differential operator이고 미분을 D문자로 표기한다. 많이 사용되지는 않는 표현법이지만 이러한 표현법이 있다는 것 정도만 알아두자. 

이러한 미분방정식은 크게 두 가지 형태로 분류할 수 있다. 첫 번째는 ODE(Ordinary differential equation), 두 번째는 PDE(Partial differential equation)이다. 각각 한글로 상미분방정식, 편미분방정식으로 부른다. ODE는 1개 혹은 여러 개의 종속 변수(dependent variable)를 가질 수 있으나, 오직 단 한개의 독립 변수(independent variable)만을 가진다. 반면에 PDE는 ODE와 달리 여러 개의 독립 변수(independent variable)를 가질 수 있다. 단 한 개의 종속 변수를 가지느냐, 아니면 여러 개를 가질 수 있냐, 이것이 가장 큰 차이점이다. 아래 식을 보자. 

식 (3.1)은 (2.2)와 같은 식이고, 여러 개의 dependent variable을 가질 수 있으나, 오직 하나의 독립 변수(independent variable)인 dx만을 가지는 형태의 식을 ODE라고 한다. 일반적으로 미분방정식이라고 하면 상미분방정식(ODE)을 의미한다. 반면 식 (3.2)와 같이 여러 개의 독립변수를 가지는 형태를 편미분방정식이라 하며, 종속 변수를 서로 다른 독립변수로 각각 편미분하기 때문에 delta인 d대신에 round인 ∂를 이용하여 표기한다. 식 (3.2)에선 하나의 종속 변수 u를 각각 x와 y로 2차 편미분 하였다. 이번 강의에서는 편미분방정식은 따로 다루지 않고 상미분방정식만을 다룬다. 

- Relationship between Differential equation and Exponential function

미분방정식(Differential equation)은 지수함수(Exponential function)와 밀접한 관련이 있다. 이번 강좌를 학습하기 위해선 이 둘의 관계를 아는 것이 좋다. 아래의 식을 보자. 

식 (4.1)은 t의 변화에 대한 y의 변화율, 즉 y를 t로 미분한것을 나타낸다. 이 변화율은 y의 값에 어떤 상수 k를 곱한 것으로 나타낼 수 있다. 이 미분방정식에서 y의 값을 구하기위해 우리는 지수함수(exponential function)를 이용할 수 있다. 식 (4.2)는 지수함수를 이용한 y의 일반해(general solution)이고, (4.2)에서 C는 y의 초기값(initial value)이며 k는 비례 상수(proportionality constant)이다. 

그렇다면 식 (4.2)는 어떻게 (4.1)의 일반해가 될 수 있는 것일까? (4.1)의 식을 약간만 변형시킨다면 (4.2)와 같이 만들 수 있다. 아래 식을 보도록 하자. 

위의 식 (5.1)은 (4.1)과 같은 식이다. 양변을 y로 나누고 dt를 곱해준 뒤, 적분(integral)을 해주면 식 (5.2)와 같이 되고, 적분을 수행한 결과가 (5.3)이다. 좌변의 적분 상수 C1을 우변으로 넘겨주고 C1+C2=C를 만든 다음 양변에 지수함수(exponential function)를 취해주면 (5.4)와 같이 된다. 지수법칙에 의해 e^c를 분리한 뒤 C로 치환해주면 식 (5.5)와 같이 정리할 수 있다. 이때 (5.4)의 y의 절대값에 의해 C는 +C와 -C의 두 가지 경우가 발생하게 되는데, C자체가 임의의 값이 될 수 있으므로 +C를 선택하여 (5.5)와 같이 정리할 수 있다. 이렇게 하여 (5.5)와 같이 미분방정식의 일반해(general solution)를 도출하였다. 

이해를 돕기위해 아래의 연습문제를 풀어보도록 하자. 

Ex) 당신은 연이율 2.4% 복리의 거치식 적금에 가입하였고 100만원을 예금하였다. 이 적금이 10년후, 50년 후엔 각각 얼마가 되어있을까? 

식 (5.5)에 이 문제를 적용하면 각 변수는 다음과 같이 설정할 수 있다. 우선 예금금 100만원은 초기값으로 C=100만원이 되고, 연이율 2.4%는 비례 상수(proportionality constant)가 되어 k=0.024가 된다. 문제에서 나타난 수치들을 식 (5.5)의 수식으로 만들어보면 다음과 같다. 

식 (6.1)은 연습 문제의 내용을 식 (4.1)의 미분방정식의 형태로 나타낸 것이다. 해(t)가 지날 때마다 내가 저축한 원금(y)이 얼마만큼씩 변하는지가 (6.1)의 좌변이 나타내는 것이고, 이 변화율은 현재의 시간 t에서의 원금 y에 이율 2.4%를 곱한 것과 같다는 것이 우변이 의미하는 것이다. 문제의 내용을 미분방정식의 형태로 정리하긴 했는데, 이를 어떻게 풀지가 막막할 것이다. 이때 활용할 수 있는 것이 식 (5.5)의 일반해(general solution)이다. 식 (6.2)는 (6.1)의 미분방정식을 일반해로 나타내어 푼 것이다. 

계산 결과 10년 후에는 약 127만원, 50년 후에는 약 332만원이 통장에 남게된다. 이 값은 근사값(approximate value)으로 실제 값과는 약간 차이가 있으나 거의 근접한 값을 구할 수 있다. 그러나 t가 커질 수록 오차는 벌어지게 된다. 이와 같이 미분방정식(differential equation)은 지수함수(exponential function)가 포함된 일반해(general solution)의 형태로 표현할 수 있으며 이를 통해 미분방정식의 해(solution)를 구할 수 있다. 이 관계를 잘 숙지하도록 하자. 

2. 미분방정식과 선형대수

- Differential equations and Linear algebra

미분방정식은 식 자체가 굉장히 복잡한 만큼 여러 가지 풀이 방법이 존재한다. 이번 챕터에서 배울 내용은 이러한 미분방정식을 선형대수의 행렬을 이용해서 푸는 방법이다. 이곳에서 다룰 내용은 상미분방정식(ODE)만 해당한다. 우선 아래의 두 개의 선형미분방정식(Linear Differential Equation)을 살펴보자. 

식 (7)은 u1과 u2에 대한 선형미분방정식을 differential form으로 각각 나타낸 것이다. 이 표기법이 익숙치 않다면 u1, u2를 x, y로 각각 생각하면 된다. u1과 u2의 미분값은 각각 서로를 포함하고 있기 때문에 이 둘은 서로 관계가 있고 따라서 하나의 벡터로 만들어서 행렬을 이용해 함께 풀어낼 수 있다. 행렬을 이용하여 풀기 위해 이들을 벡터로 정의하면 아래와 같다.  

여기서 우리가 구하고자 하는 것은 식 (7)의 두 개의 미분방정식을 행렬을 이용하여 1차 식으로 만들어 푸는 것이다. 이를 위해 u1과 u2를 원소(component)로 하는 벡터 u(t)를 식 (8.1)과 같이 정의할 수 있다. (8.2)는 u를 미분한 du/dt를 나타낸 것이고, (8.3)은  t=0일때의 초기 조건(initial condition)에서의 값을 의미한다. 

u와 u'을 정의했으니 이제 남은 행렬 A를 정의하면 우리가 구하고자 하는 미분방정식의 행렬 형태를 정의할 수 있다. A는 어떻게 정할 수 있을까? 이를 위해 먼저 아래와 같이 식을 나열해보자. 

식 (9)의 u1'은 (7.1)의 좌변에 표현된 u1의 미분값에 해당한다. 이 u1'이 되기 위해선 우변의 u1과 u2에 각각 얼마씩이 곱해져야하는지를 생각해보면 A의 row1을 알 수 있다. -1과 2를 각각 u1과 u2에 곱하면 u1'이 된다. 따라서 A의 row1은 -1, 2가 되고, 마찬가지로 식 (7.2)를 참고하여 row2를 구해보면 1, -2가 되는 것을 알 수 있다. 

이렇게 하여 식 (9)를 다시 써보면 아래와 같다. 

위와 같이 식 (7)의 두 개의 미분방정식을 행렬에 대한 식으로 표현하였다. 이 식을 통해 두 미분방정식에 대하여 무엇을 알 수 있을까? 우리는 식 (7)에 표현된 u1과 u2에 대한 미분방정식이 시간이 지남에 따라 어떻게 변하는지를 식 (9)를 통해 알 수 있다. A의 고유값(eigenvalue)과 고유벡터(eigenvector)를 통해서 말이다. A의 고유값과 고유벡터를 구해보자. 

식 (10.1)을 통해 행렬 A의 고유값이 각각 0과 -3임을 알았다. A가 2x2이면서 특이행렬(singular matrix)이기 때문에 두 개의 고유값중 하나가 0이 나오는 것이다. 이 특성을 알고 2x2행렬에선 trace가 모든 고유값들의 합과 같다는 성질을 이용하여 계산과정 없이 고유값을 빠르게 구할 수도 있다. 다음으로 고유벡터를 구하기 위해 (A-λI)x=0 의 해인 null space를 구해야 한다. 이때의 행렬이 특이 행렬이므로 x2가 free variable이다. 이 free variable을 1로 설정한 뒤 나머지 해를 구하면 (10.2)와 (10.3)과 같이 고유벡터를 구할 수 있다. 

- General Solution

식 (9)의 행렬 A에 대한 고유값(eigenvalue)과 고유벡터(eigenvector)를 구했으니 이제 이들을 활용하여 u(t)에 대한 해를 구해보도록 하자. A는 2x2크기의 행렬이므로 두 개의 고유값과 고유벡터를 가진다. 따라서 u(t)에 대한 일반해(general solution)는 두 개의 special solution, 즉 두 개의 exponential solution의 조합으로 표현할 수 있다. 아래의 식을 보자. 

식 (11)은 (9)에 나타난 u(t)의 일반해(general solution)를 나타낸 것이다. 식 (5)에서 보인것과 같이 미분방정식의 일반해를 exponential solution의 형태로 보이긴 했는데, 과연 이들이 미분방정식의 해의 형태가 될 수 있는가에 대한 의문이 생길 수 있다. 우리가 최종적으로 구하고자 하는 du/dt = Au에 exponential solution을 대입해서 이를 확인해보도록 하자. 

(12.1)에서와 같이 초기값을 고려하지 않은 순수 지수함수(pure exponential)를 u로 여겨 원래의 식에 대입하여 정리하면 (12.2)와 같이 정리할 수 있다. 좌변은 u의 pure exponential을 t로 미분한 결과이며, 따라서 exponential의 지수부인 람다(λ)가 앞으로 나와 계수로 곱해진 것이다. (12.2)의 양변을 지수함수로 나누어주면 결과적으로 (12.3)과 같이 정리할 수 있으며, 이는 앞서 공부했던 고유값/고유벡터의 정의와 같은 것을 볼 수 있다. 마찬가지로 λ2에 대한 pure exponential을 대입 해도 해를 구할 수 있다. 이를 통해 우리가 구하고자 하는 (12.1)의 행렬로 정의된 미분방정식의 해를 지수함수(exponential function)를 통해 구할 수 있음을 증명한 셈이다. 

식 (11)에 정의된 지수함수(exponential function)로 정의된 미분방정식의 일반해는 지난 강의(Lecture 22)에서 배웠던 대각화(diagonalization)에서의 일반해와 같은 꼴로 생각할 수 있다. 단지 t로 정의된 연속적인(continuous) 개념인지, k로 정의된 이산적인(discrete) 개념인지가 다르다. 아래의 두 식은 2x2크기 행렬의 미분방정식(differential equation)과 차분방정식(difference equation)의 일반해를 각각 나타낸다.  

이제 지수함수와 벡터는 알겠는데, 각 항의 앞에 곱해진 상수 c1과 c2는 어떻게 알 수 있을까? 이들 상수는 초기 조건(initial condition)으로부터 계산할 수 있다. 이는 앞부분에 미분방정식의 개념을 설명할 때 배웠던 식 (4)의 일반해의 초기값(initial value) C에 해당하는 값이다. 식 (8.3)의 초기값으로부터 c1과 c2를 구해보도록 하자.  

우선 (14.1)은 앞서 구한 고유값과 고유벡터를 대입하여 일반해의 식을 정리한 것이다. 고유값이 각각 λ1=0, λ2=-3이므로 첫 번째 지수 함수는 1이 되고, 두 번째 지수함수의 지수부는 -3t가 되어 (14.1)과 같이 정리할 수 있다. 이때 t=0일 때의 초기값 u가 u(0)=[1 0]T이므로 t에 0을 대입하여 정리하면 (14.2)와 같이 정리할 수 있다. c1과 c2는 각각 column vector에 곱해진 계수가 되기 때문에 이 column vector(고유 벡터들)들을 행렬로 정리하면 (14.3)과 같이 나타낼 수 있으며, 이때의 행렬을 고유벡터 행렬(eigenvector matrix) S라고 한다. S를 가우스 소거(Gauss elimination)를 통해 c1과 c2를 계산하면 각각 c1=1/3과 c2=-1/3이 된다. (가우스 소거법은 Lecture 2참조)

고유벡터 행렬 S를 통해 우리는 해(solution)의 계수(coefficient) c1, c2, ...를 구할 수 있으며, 이 계수들은 해에서 각 순수 지수함수(pure exponential)들이 얼마만큼인지를 나타낸다. 이 계수들은 t=0일 때의 u인 u(0)로부터 구할 수 있으며, 이렇게 구한 계수들은 t가 무한대로 갔을 때 경우에 따라 극한값이 되기도 한다. 

c1과 c2를 구했으니 이들을 일반해 식에 대입하여 정리하면 아래와 같다. 

식 (15)가 (9)의 미분방정식의 일반해(general solution)이다. 일반해를 도출하기까지의 과정을 살펴보면 먼저 미분방정식을 행렬의 형태로 만든 다음, 초기 조건(initial condition)에 따른 초기 값을 구하고 행렬의 고유값과 고유벡터를 구한 뒤, 고유벡터행렬을 통해 각 term의 계수값들(c1, c2, ...)을 구한다. 행렬 A의 고유값, 고유벡터, 계수 등 모든 값들을 구하여 일반해 형태의 식에 대입하고 정리하여 최종적인 해를 구할 수 있다. 

식 (15)와 같이 일반해를 구했기 때문에 시간이 지남에 따라 미분방정식의 값이 어떻게 변하는지를 쉽게 알아낼 수 있다. 시간에 대한 변수인 t값만 대입해주면 해당 시간에서의 값을 알 수 있다. 그런데 우리가 세운 미분방정식이 만약 어떤 기계의 진동에 따른 기계 수명 예측과 같은 문제라고 한다면, 경영진 입장에서는 당연히 아주 오랜 시간이 지났을 경우 기계의 기대 수명이 궁금할 것이다. 즉 우리가 여기서 알고자 하는 것은 시간이 무한대로 지났을 때 미분방정식의 값이다. 이는 앞서 구했던 미분방정식 행렬의 고유값과 관련이 깊다. 다음 파트에서 고유값과 미분방정식의 무한대 시간에서의 해(infinite solution)의 관계에 대해 알아보자. 

- Stability

우리는 지난 대각화 강의에서 고유값을 통해 해당 차분방정식의 특성을 파악할 수 있음을 언급했다. 마찬가지로 이번 강의에서 다뤘던 미분방정식 행렬의 경우에도 고유값을 통해서 관련 특성을 파악할 수 있다. 즉 식 (9)의 A의 고유값을 통해 우리는 행렬 A에 내포된 미분방정식들이 시간이 지남에 따라 어떻게 변하는지를 예측할 수 있다. 그리고 궁극적으로 시간이 무한대로 지났을 때 그 극값이 어떻게 되는 지를 알 수 있다. 

식 (10)에서 계산한 첫 번째 고유값은 0이고, 이 고유값은 식 (11)에서 일반해의 첫 번째 term인 c1*exp(0*t)*x1의 지수 함수의 지수부가 된다. 따라서 지수함수는 t가 아무리 증가해도 0의 거듭제곱 꼴이 되기때문에 항상 1이 되므로 식 (15)에서 첫 번째 term에서 지수함수 부분이 사라진 것이다. 이는 결국 미분방정식이 시간이 지날수록 정상상태(steady state)로 수렴함을 의미한다. 

두 번째 고유값은 -3이고 식 (11)에서 c2*exp(-3*t)*x2의 꼴이 되기 때문에 t가 커질수록 지수함수의 지수부는 작아지고 점점 작은 값이 된다. 따라서 두 번째 term의 식 자체가 시간이 증가할수록 굉장히 작은 값이 된다. 결과적으로 우리가 구한 일반해 (15)의 극한값(t가 무한대로 갔을 때의 값)은 아래와 같다. 

식 (16.2)가 극한 값이며, 실제 값은 [0.666... , 0.333...]T이다. 식 (7.1)의 미분방정식의 t가 무한대일때의 극한값은 0.666..., 식 (7.2)의 극한값은 0.333...이다. 이들은 2개의 미분방정식 (7.1), (7.2)가 연립된 형태인 시스템 (9)에 대한 해 벡터(solution vector)의 극한값 u(t=∞)=[0.666..., 0.333...]T이다. 

미분방정식은 대수방정식(Algebraic equation)처럼 어떤 값이나 값들의 집합의 형태로 해가 나오는 것이 아니다. 미분방정식의 해는 (15)와 같이 식의 형태로 나오며, MATLAB과 같은 프로그램을 통해 시간이 지남에 따라, 혹은 종속 변수(dependent variable)가 변함에 따라 해가 어떤 식으로 변하는지를 알 수 있다. 아래 그래프는 (15)의 해가 시간이 지남에 따라 변하는 과정을 나타낸 것이다. 

Fig. 1 미분방정식의 해가 시간이 지남에 따라 변화하는 과정

Fig. 1에서 빨간색 선과 파란색 선이 각각 일반해의 u1과 u2를 의미한다. t=0에서부터 시간이 지남에 따라 특정 값에 수렴하는 것을 볼 수 있는데, 최종적으로 수렴하는 값이 u1=0.666..., u2=0.333... 임을 알 수 있다. 이는 (16)에서 구한 극값과 일치한다. uf(first)라고 표시된 벡터는 일반해에서 첫 번째 고유값인 0이 속해있는 식을 의미한다. 녹색과 자홍색 그래프가 uf의 u1과 u2에 각각 해당하는데, 시간이 지나도 어떤 특정값에서 변화가 없음을 확인할 수 있다. 이는 고유값이 0이기 때문에 지수함수 부분이 항상 1이 되기 때문이다. 

us(second)라고 표기된 벡터는 일반해에서 두 번째 고유값인 -3이 속해있는 식을 의미한다. 지수부의 고유값이 음수이기 때문에 시간이 지날수록 값이 작아지는 것을 볼 수 있다. 하늘색과 검은색 그래프가 초기 값에서 시간이 지날수록 0으로 수렴하는 것을 볼 수 있다. 

이와 같이 미분방정식을 행렬을 이용하여 풀 때 행렬의 고유값에 의해서 우리는 연립미분방정식의 해가 시간이 지남에 따라 정상상태(steady state)수렴할 수도, 값이 0이 되어 사라질 수도 혹은 값이 무한대로 커져서 발산할수도 있다. 그렇다면 고유값이 어떤 상태일 때 발산 혹은 수렴이 될까? 아래의 정리를 살펴보자. 

(1) Stability	(2) Steady State	(3) Divergence

미분방정식의 해는 크게 세 가지 상태가 될 수 있는데, 위의 표에 나타난 것과 같이 0으로 수렴하는 안정화 상태인 stability, 어떤 특정 값으로 수렴하는 정상 상태(steady state) 그리고 아예 값 자체가 무한정 커지는 divergence상태가 있다. 

만약 모든 고유값이 0보다 작다면 stability 상태가 된다. 이때는 해(solution) u(t)가 시간이 지날수록 0으로 수렴하게 되는데, 모든 고유값이 0보다 작은 음수인 경우 모든 지수함수의 값들이 점점 작아지기 때문이다. 여기서 고유값을 실수부(real part)라고 지정한 것은 해가 허수부(imaginary part)를 포함할 수도 있기 때문이다. 아래의 예를 보자. 

식 (17.1)은 허수부를 포함한 지수함수이다. 지수법칙에 의해 실수부와 허수부를 나눌 수 있는데, 이때 허수부의 크기는 오일러 공식에 의해 항상 1이 되기 때문에 t가 커진다고 해도 해에 실질적인 영향은 주지 않는다. 허수부는 오일러 공식 위키(https://en.wikipedia.org/wiki/Euler%27s_formula)에서 볼 수 있듯이 t가 증가해도 그저 복소 평면에서 단위 원(unit circle)을 반복해서 돌 뿐이다. 결국 Stability는 고유값의 실수부(real part)에 의해서만 결정된다. 

Stability를 만족시키는 조건은 2x2 크기의 행렬에서는 trace와 행렬식(determinant)로 쉽게 구할 수 있다. 2x2행렬에선 trace(A)=λ1+λ2, det(A)=λ1*λ2 임을 이용하여 trace < 0와 det > 0를 만족한다면 두 개의 고유값이 모두 음수이므로 Stability상태가 됨을 쉽게 알 수 있다. 

두 번째로 정상 상태(Steady state)의 경우는 해가 어떤 특정 값으로 수렴하는 경우를 의미하는데, (15), (16)에서 구한 해가 바로 이 경우 이다. 즉 어떤 하나의 고유값이 0이고 나머지 고유값이 0보다 작은 경우에 해당하는데, (16)의 경우에 첫 번째 고유값이 0이고, 두 번째 고유값이 0보다 작은 -3이므로 시간이 지날수록 두 번째 term은 0으로 수렴하여 사라지고 첫 번째 term의 상수인 c1=1/3과 첫 번째 고유값 x1=[2 1]T의 곱의 값으로 수렴하게 된다. Fig. 1의 빨간색과 파란색 그래프가 바로 이 특정 값으로 수렴하는 모습을 보여준다. 

마지막의 경우는 어떤 고유값이 하나라도 0보다 큰 경우 시간이 지날 수록 무한대로 커져서 발산하게 된다. 

종합해보면 복소수까지 고려했을 때 오직 실수부의 고유값만이 미분방정식의 해에 영향을 미치고, 고유값(eigenvalue)이 0보다 크고 작음, 혹은 같음에 따라 해의 상태가 결정된다. 

아래 그림은 Fig. 1을 plot하기 위한 MATLAB 소스코드이다. 

3. 마치며(continue)

이번 강의에선 미분방정식의 간단한 소개와 함께 선형대수를 이용하여 미분방정식을 푸는 방법을 배웠다. 기본적으로 미분방정식에 대한 식을 세운 뒤 고유값과 고유벡터를 구하여 일반해에 대한 식을 세운 다음 초기 조건으로부터 일반해의 계수를 구하는 순서로 해를 구하게 된다. 해는 지수함수의 고유값의 거듭제곱꼴로 나타나는데, 고유값이 0보다 크거나 같음 또는 작음에 따라 해의 상태가 결정됨을 배웠다. 다음 포스팅에서는 이번 강의에 이어서 미분방정식과 선형대수의 나머지 이야기를 하도록 하겠다.