기초 개념

• 고유값 $\lambda$, 고유벡터 $x$: $Ax = \lambda x$를 만족하는 값.

• 스펙트럼: 역행렬이 있는 행렬 A는 n개의 고유값과 n개의 고유벡터를 가지는데, 이 때 모든 고유값의 집합 $\{ \lambda_1, \lambda_2, \cdots, \lambda_n \}$을 스펙트럼이라고 함.

• 최소고유값, 최대고유값: $|\lambda_1| \leq |\lambda_2| \leq \cdots \leq |\lambda_n|$으로 정렬했을 때 $\lambda_1$을 최소고유값, $\lambda_n$을 최대고유값이라고 함.

• 스펙트럼 반지름: 최대고유값을 스펙트럼 반지름이라고도 함.

• 대수적다중도(algebraic multiplicity): k개의 고유값이 서로 같은 값을 가질 때 k를 대수적다중도라고 함.

• 기하적다중도(geometric multiplicity): 고유값이 k개의 고유벡터를 가질 때 k를 기하적다중도라고 함.

• 특성행렬식: $\det(A - \lambda I) = 0$

• 특성다항식: $\det(A - \lambda I)$를 $\lambda$에 관한 n차다항식 $D(\lambda)$로 본 경우

$$\begin{align*} \det(A - \lambda I) &= (a_{11} - \lambda)(a_{22} - \lambda) \cdots (a_{nn} - \lambda) + \\ & (a_{ii} - \lambda) \text{꼴의 항이 최대 n-2개 있는 항들} \\ &= (\lambda_1 - \lambda)(\lambda_2 - \lambda) \cdots (\lambda_n - \lambda) \\ &= 0 \end{align*}$$

따라서 $\lambda_1 \lambda_2 \cdots \lambda_n = \det(A - 0 \cdot I) = \det A$,
$\lambda_1 + \lambda_2 + \cdots + \lambda_n = \lambda^{n-1} \text{의 계수} = a_{11} + a_{22} + \cdots + a_{nn} = \operatorname{trace} A$

대각화

$A^* = B^{-1}AB$를 A의 상사행렬(similar matrix)이라고 하고, 이러한 변환을 상사변환(similarity transformation)이라고 함.
이때 고유벡터를 열벡터로 갖는 행렬 X, 고유값을 주대각 원소로 갖는 대각행렬 D를 생각하자.
$AX = XD$이므로, $D = X^{-1}AX$가 되고, 이를 행렬의 대각화라고 한다.

고유값의 성질

• $A^k$의 고유값은 $\lambda_i^k$
• $\alpha A$의 고유값은 $\alpha \lambda_i$
• $A - \alpha I$의 고유값은 $\lambda_i - \alpha$
• 상사행렬 $A^* = B^{-1}AB$의 고유값은 $\lambda_i$ (즉 고유값은 똑같음)
• A가 삼각행렬이면 고유값은 주대각원소 $a_{ii}$
  이는 특성행렬식을 전개했을 때 $(a_{11} - \lambda)(a_{22} - \lambda) \cdots (a_{nn} - \lambda)$가 돼서 자명하다.

행렬식 탐색법

케일리-해밀튼 정리: 특성다항식 $D(\lambda) = 0$을 행렬다항식으로 봤을 때, 원래 행렬 A는 이 다항식의 근이 됨. 즉 $D(A) = 0$을 만족함.

즉 아래 두 식이 모두 성립한다.

$$\begin{align*} D(\lambda) &= -\lambda^n + c_1 \lambda^{n-1} + \cdots + c_{n-1} \lambda + 1 = 0 \\ D(A) &= -A^n + c_1 A^{n-1} + \cdots + c_{n-1} A + 1 = 0 \\ \end{align*}$$

Fadeev-Leverrier의 방법: 아래 점화식으로 특성방정식 계수를 모두 구할 수 있음

$$\begin{aligned} B_1 &= A & c_1 &= \operatorname{trace} B_1\\ B_2 &= A(B_1 - c_1 I) & c_2 &= \frac{1}{2} \operatorname{trace} B_2 \\ &\vdots & &\vdots \\ B_n &= A(B_{n-1} - c_{n-1} I) & c_n &= \frac{1}{n} \operatorname{trace} B_n \end{aligned}$$

거쉬고린(Gershgorin) 정리: eigenvalue의 범위를 특정할 수 있음

$$|z - a_{ii}| \leq \sum_{j=1, j \neq i}^{n} |a_{ij}|$$

각 주대각원소마다 거쉬고린 원(Gershgorin's circle)이 n개 나오는데, 모든 eigenvalue는 이 n개의 원의 합집합에 속해있다.

단점: 이렇게 열심히 해봤자 겨우 특성방정식만 구할 수 있음.. 특성방정식의 근은 따로 구해야 함

거듭제곱법

1. 고유벡터의 초기 추측값 $x^{(0)}$을 정함. (보통 1만 있는 벡터로 정함)
2. $y^{(k+1)} = Ax^{(k)}$를 구하고 $y^{(k+1)}$에서 제일 절대값이 큰 원소를 $\lambda^{(k+1)}$이라고 한다.
3. $x^{(k+1)} = \frac{1}{\lambda^{(k+1)}}y^{(k+1)}$으로 다시 절대값을 최대 1로 만든 다음 위 과정을 반복한다.
4. $\lambda, x$ 중 한 쪽이 허용오차 범위에 들어오면 종료한다.

위 방법으로 최대고유값 $\lambda_n$과 고유벡터 $x_n$을 쉽게 구할 수 있지만, $\lambda_{n-1}$과 비슷하거나 중근인 경우 (즉, 둘이 같은 경우) 수렴하지 않는다.

원리) 임의의 벡터를 고유벡터들의 선형결합으로 쓸 수 있으므로, $x^{(0)} = c_1​v_1​ + c_2​v_2​ + \cdots + c_n​v_n$으로 쓸 수 있다.

$$\begin{align*} \therefore A^kx^{(0)} &= A^k (c_1​v_1​ + c_2​v_2​ + \cdots + c_n​v_n) \\ &= c_1 \lambda_1^k v_1 + c_2 \lambda_2^k v_2 + \cdots + c_n \lambda_n^k v_n \\ &\approx c_1 \lambda_1^k v_1 \end{align*}$$

역거듭제곱법

$Ax = \lambda x \iff A^{-1}x = \frac{1}{\lambda}x$

따라서 $A^{-1}$에서 거듭제곱법을 쓰면 최소고유값을 (정확히는 최소고유값의 역수를) 구할 수 있다.
하지만 역행렬을 구하는 것은 비효율적이므로 대신 $Ay^{(k+1)} = x^{(k)}$를 LU분해나 가우스 소거법을 이용하여 구한다.

고유값 이동법

$Ax = \lambda x \iff (A - sI)x = (\lambda - s)x$

따라서 $A - sI$의 고유값은 $\lambda -s$이므로, $A - sI$로 고유값을 이동시켜 수렴속도를 가속시키거나 새로운 고유값을 구할 수 있다.
예를 들어 거듭제곱법으로 최대고유값 $\lambda_n$을 구했다면, $A - \lambda_n I$에서 거듭제곱법을 하면 다른 고유값을 구할 수 있다.

고유벡터의 결정

고유값 $\lambda$를 구해도 고유벡터를 구하려면 $(A - \lambda I) x = 0$을 풀어야 한다.
하지만 이 방정식은 무한히 많은 해를 가지므로, 미지상수가 필요하다.
e.g. 고유벡터가 $x = t(0, 1, 1) + s(1, 0, -1)$같은 식일 수도 있다.

$B = A - \lambda I$라고 하고, $Bx = 0$의 해를 구해보자.

1. 가우스 조단 소거 후, B의 k번째 열이 0벡터라면 $e_k$는 고유벡터이다. 만약 B의 n번째 열이 0벡터라면, 여기서 그만둔다.
2. 만약 $b_{nn} = 0$이라면 $x_n = 0, -1$ 두 가지 경우 모두 조사한다.

• $x_n = 0$인 경우, 마지막 행/열을 제거한 (n-1) x (n-1) 행렬에 대해서 다시 반복한다.
• $x_n = -1$인 경우, 마지막 열을 소스벡터로 두고 (n-1) x (n-1) 행렬에 대해서 연립방정식을 다시 푼다. (즉, $B_{n-1} x_{n-1} = b_{n}$)

1. 만약 $b_{nn} = 1$이라면 $x_n = 0$으로 두고 조사한다.
2. 최종 고유벡터가 영벡터라면 고유값이 틀렸다고 판단한다.

하우스홀더 변환

$$A^* = PAP, P = I - 2\frac{uu^T}{|u|^2}$$

여기서 P를 하우스홀더 행렬이라고 하고, $P = P^T = P^{-1}$을 만족한다.
하우스홀더 변환을 n-2번 하면 헤센버그 행렬(주대각 바로 아래의 대각원소를 제외하고 대각선 아래의 원소가 전부 0인 행렬)이 된다. 특히 A가 대칭이면 (A = A^T), 3대각행렬이 된다.
c.f. 상사행렬 꼴이므로 이 방법으로 기존 행렬과 고유값이 같은 헤센버그 행렬을 구할 수 있다.

하우스홀더 변환이 성립하는 u를 찾기 위해 $\alpha = \frac{2}{|u|^2}$, $\beta = u^T A u$라고 하면,

$$\begin{align*} A^* &= (I - \alpha uu^T)A(I - \alpha uu^T) \\ &= (I - \alpha uu^T)(A - \alpha Auu^T) \\ &= A - \alpha(uu^T A + Auu^T) + \alpha^2 uu^T Auu^T \\ &= A - \alpha(uu^T A + Auu^T) + \alpha^2 \beta uu^T \\ a_{ij}^* &= a_{ij} - \alpha \sum_{s=1}^n (u_iu_sa_{sj} + a_{is}u_su_j) + \alpha^2 \beta u_iu_j \end{align*}$$

A가 4x4 행렬인 경우를 가정하면, $u^T = \begin{bmatrix}0 & u_2 & u_3 & u_4\end{bmatrix}$로 잡고, $A^*$의 첫 열 아래쪽을 0으로 만들어야 한다. (즉, $a_{31}^* = a_{41}^* = 0$)
$uu^T$의 첫 행과 첫 열은 전부 0이므로, $A^*$의 첫 열과 관계가 없다.
따라서 $A - \alpha(uu^TA)$의 첫 열이 $A^*$의 첫 열이 된다.

$$\begin{align*} a^*_{31} &= 0 = a_{31} - \alpha u_3 (u_2 a_{21} + u_3 a_{31} + u_4 a_{41}) \\ a^*_{41} &= 0 = a_{41} - \alpha u_4 (u_2 a_{21} + u_3 a_{31} + u_4 a_{41}) \\ \therefore \frac{a_{31}}{u_3} &= \frac{a_{41}}{u_4} = \alpha (u_2 a_{21} + u_3 a_{31} + u_4 a_{41}) \end{align*}$$

따라서 $u_3 = a_{31}, u_4 = a_{41}$로 놓으면 정의에 의해 $\frac{1}{\alpha} = \frac{|u|^2}{2} = \frac{u_2^2 + a_{31}^2 + a_{41}^2}{2}$이므로 위 방정식은 $(u_2 - a_{21})^2 = a_{21}^2 + a_{31}^2 + a_{41}^2$가 된다.

$$\therefore u_2 = \gamma + a_{21}, \gamma = \operatorname{sign}(a_{21})\sqrt{a_{21}^2 + a_{31}^2 + a_{41}^2}$$

마무리오차를 줄이기 위해 절대값이 큰 쪽을 택한다.

일반적으로 n x n 행렬 A의 k번째 반복 과정에서 하우스홀더 변환을 위한 u벡터는

$$\gamma = \operatorname{sign}(a_{k+1, k})\sqrt{\sum_{j=k+1}^n a_{j,k}^2}, u_i = \begin{cases} 0 & 1 \leq i \leq k \\ a_{k+1, k} + \gamma & i = k+1 \\ a_{i, k} & i > k + 1 \end{cases}$$

이며, 하우스홀더 행렬은 아래와 같다.

$$\alpha = \frac{2}{|u|^2} = \frac{1}{\gamma^2 + \gamma a_{k+1, k}}, P = I - \alpha uu^T$$

QR 반복법

상블록 삼각행렬(Upper-block-triangular matrix)는 행렬 A를 여러 개의 블록으로 나누었을 때 (일반적으로 1x1이랑 2x2로만 나눔), 블록 단위로 삼각행렬인 행렬이다.

$$A = \begin{bmatrix} A_{11} & A_{12} & A_{13} & \cdots & A_{1k} \\ 0 & A_{22} & A_{23} & \cdots & A_{2k} \\ 0 & 0 & A_{33} & \cdots & A_{3k} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & A_{kk} \end{bmatrix}$$

상블록 삼각행렬의 경우 det를 간단하게 각 대각블록의 det의 곱으로 나타낼 수 있다.

$$\det(A) = \det(A_{11})\det(A_{22}) \cdots \det(A_{kk})$$

따라서 상블록 삼각행렬의 특성다항식도 1~2차 다항식들의 곱으로 나타낼 수 있으므로, 고유값을 곧바로 구할 수 있다.

이미 하우스홀더 변환으로 기존 행렬 A와 고유값이 같은 헤센버그 행렬 $A_0$을 구하는 방법을 배웠다.
따라서 헤센버그 행렬을 상블록 삼각행렬로 변환시키는 방법만 구하면 된다.

$A_k = Q_kR_k$로 QR분해 한 뒤, $A_{k+1} = Q_k^{-1}A_kQ_k$인 상사변환을 정의하면
$A_{k+1} = Q_k^{-1}Q_kR_kQ_k = R_kQ_k$고, 상사변환이니 고유값이 같고, $A_k$가 헤센버그 행렬이면 $A_{k+1}$도 헤센버그 행렬이 된다.
이때 이 과정, 즉 QR분해 후 행렬을 바꿔서 RQ로 곱하는 과정을 반복하면 $A_{k+1}$이 점점 상블록 삼각행렬로 수렴하게 된다.

고유값 이동

$A_k - s_kI = Q_kR_k, A_{k+1} = R_kQ_k + s_kI$ 처럼 고유값 이동을 했다가 취소하는 방법으로 수렴을 더 가속시킬 수 있다.
고유값이 $s_k$에 가까울수록 수렴이 가속된다.

$R_k = Q_k^T(A_k - s_kI)$를 대입하면 $A_{k+1} = Q_k^{-1}A_kQ_k$, 즉 고유값이 보존된다는 것을 보일 수 있다.

이중 QR 분해

일반적으로 고유값은 복소수이므로, 실수를 원소로 가지는 행렬에 고유값 이동을 그대로 적용하기엔 어려움이 따른다.
하지만 켤레복소수도 합과 곱이 실수라는 점을 사용하면 이러한 문제를 해결할 수 있다.

$$\begin{align*} A_k - s_kI &= Q_kR_k \\ A_{k+1} - s_{k+1}I &= Q_{k+1} R_{k+1} \\ A_{k+1} &= Q_k^{-1}A_kQ_k \\ \therefore A_k - s_{k+1}I &= Q_k(A_{k+1} - s_{k+1}I)Q_k^{-1} \\ &= Q_kQ_{k+1} R_{k+1}Q_k^{-1} \\ \therefore (A_k - s_{k+1}I)(A_k - s_kI) &= Q_kQ_{k+1} R_{k+1}R_k \end{align*}$$

이때 $A_{k+2} = Q_{k+1}^{-1}A_{k+1}Q_{k+1} = Q_{k+1}^{-1}Q_k^{-1}A_kQ_kQ_{k+1}$이다.
따라서 $(A_k - s_{k+1}I)(A_k - s_kI) = QR$로 QR분해 한 뒤, $A_{k+2} = Q^TA_kQ$로 QR반복법을 수행할 수 있다.

$s_k, s_{k+1}$의 값으로는 고유값의 합은 trace, 곱은 det인걸 이용하여, $A_k$의 오른쪽 아래 2x2 행렬을 $B_k$라고 하면
$(A_k - s_{k+1}I)(A_k - s_kI) = A_k^2 - \operatorname{trace} B_k \cdot A_k + \det B_k \cdot I$로 정할 수 있다.

프로베니우스 행렬

프로베니우스 행렬(Frobenius matrix)로 알려진 다음의 행렬

$$\mathbf{A} = \begin{bmatrix} -a_{n-1} & -a_{n-2} & \cdots & -a_1 & -a_0 \\ 1 & 0 & \cdots & 0 & 0 \\ 0 & 1 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & 1 & 0 \end{bmatrix}$$

의 특성방정식은 $\lambda^n + a_{n-1} \lambda^{n-1} + \cdots + a_1\lambda + a_0 = 0$임이 알려져있다.

따라서 이 행렬의 고유값을 구하는 문제와 임의의 n차 다항식의 근을 구하는 문제는 동치다.
프로베니우스 행렬은 이미 헤센버그 행렬이므로, QR 반복법을 이용하여 특성방정식의 근을 구할 수 있다.
이 방법이 n차 다항식의 근을 구하는 방법 중에서 제일 우수한 것으로 알려져 있다!

호텔링 수축법

(시험범위 아님!!)
점점 오래 걸리고 오차도 많이 쌓이지만 처음 몇 개 구할때는 괜찮음

대칭행렬 $A = A^T$의 경우 서로 다른 고유벡터는 수직(orthogonal)하다.
이를 이용해 대칭행렬의 고유값을 쉽게 구할 수 있다.

왠진 모르겠지만 이번엔 $|\lambda_1| \geq |\lambda_2| \geq \cdots \geq |\lambda_n|$으로 정렬해보자.
거듭제곱법으로 $\lambda_1, x_1$을 구할 수 있다. (고유벡터는 정규화되었다고 가정)

이때 $A_2 = A - \lambda_1x_1x_1^T$라고 하면,

$$\begin{align*} A_2x_1 &= Ax_1 - \lambda_1x_1x_1^Tx_1 = Ax_1 - \lambda_1x_1 = 0 \\ A_2x_j &= Ax_j - \lambda_1x_1x_1^Tx_j = Ax_j = \lambda_j x_j (j \neq 1) \end{align*}$$

따라서 $A_2$의 고유값은 $0, \lambda_2, \ldots, \lambda_n$이고, 고유벡터는 $x_1, x_2, \ldots, x_n$이다.
그런데 $\lambda_1$은 최대고유값이였으므로, $A_2$의 고유값은 $|\lambda_2| \geq \cdots \geq |\lambda_n| \geq 0$이 된다.
따라서 $A_2$에 거듭제곱법을 쓰면 $\lambda_2, x_2$를 구할 수 있고, 이를 반복해서 모든 고유값을 구할 수 있다.

하지만 이전 단계에서 구한 고유값과 고유벡터를 계속 사용하기 때문에 오차가 누적된다는 단점이 있어 대형 계산에는 적합하지 않다.

자코비 방법

대칭행렬을 회전변환을 사용한 상사변환만을 사용하여 대각화시키는 방법이다.

$$A = \begin{bmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{bmatrix}, Q = \begin{bmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{bmatrix}, A^* = Q^{-1}AQ = \begin{bmatrix} a_{11}^* & a_{12}^* \\ a_{21}^* & a_{22}^* \end{bmatrix}$$

라고 할 때, $a_{12}^* = a_{21}^* = a_{12}\cos 2\theta + (a_{22} - a_{11}) \frac{\sin 2\theta}{2} = 0$일 조건은

$$\begin{cases} \tan 2\theta = \frac{2a_{12}}{a_{11} - a_{22}} & a_{11} \neq a_{22} \\ \theta = \frac{\pi}{4} & a_{11} = a_{22} \end{cases}$$

이다. 이때 상사행렬 $A^*$는 원래의 행렬 $A$와 동일한 고유값을 가지므로, $\lambda_1 = a_{11}^*, \lambda_2 = a_{22}^*$이다.

n x n행렬로 확장하면 계속 2x2 행렬을 찾아서 대각화 시키는 것을 반복해 최종적으로 전체 행렬을 대각화시킨다.

1. $X_0 = I, A_0 = A$로 둔다.
2. r, s를 $|a_{rs}| = \max_{i \neq j} |a_{ij}|$로 잡는다. $|a_{rs}| < \epsilon$이면 종료한다.
3. $Q_k$를 I에서 r, s번째 행, 열만 $\begin{bmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{bmatrix}$로 바꾼 행렬로 정의한다.
4. $\theta$를 $\begin{cases} \tan 2\theta = \frac{2a_{rs}}{a_{rr} - a_{ss}} & a_{rr} \neq a_{ss} \\ \theta = \frac{\pi}{4} & a_{rr} = a_{ss} \end{cases}$로 정한다.
5. $A_{k+1} = Q_k^T A_k Q_k, X_{k+1} = X_k Q_k$로 놓고 반복한다.
   실제로 행렬곱을 계산할 필요는 없고 아래 값들만 바꾸면 된다.

$$\begin{align*} a_{rs}^* &= a_{sr}^* = 0 \\ a_{rr}^* &= a_{rr}\cos^2\theta + a_{ss}\sin^2\theta + 2a_{rs}\cos\theta\sin\theta \\ a_{ss}^* &= a_{rr}\sin^2\theta + a_{ss}\cos^2\theta - 2a_{rs}\cos\theta\sin\theta \\ a_{ir}^* &= a_{ri}^* = a_{ir}\cos\theta + a_{is}\sin\theta (i \neq r, s) \\ a_{is}^* &= a_{si}^* = a_{is}\cos\theta - a_{ir}\sin\theta (i \neq r, s) \\ x_{ir}^* &= x_{ir}\cos\theta + x_{is}\sin\theta \\ x_{is}^* &= x_{is}\cos\theta - x_{ir}\sin\theta \end{align*}$$

끝나면 고유값과 고유벡터가 한번에 구해지지만, 행렬의 크기가 커지면 연산이 과도하게 필요하다.

대칭행렬의 고유값, 고유벡터

대칭행렬의 고유값은 모두 실수이고, 하우스홀더 변환시 3대각행렬이 된다.

$$PAP = A^* = \begin{bmatrix} a_1 & b_1 & 0 & \cdots & 0 & 0 \\ b_1 & a_2 & b_2 & \cdots & 0 & 0 \\ 0 & b_2 & a_3 & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & a_{n-1} & b_{n-1} \\ 0 & 0 & 0 & \cdots & b_{n-1} & a_n \end{bmatrix}$$

위 특성방정식은 $p_k(\lambda) = (a_k - \lambda)p_{k-1}(\lambda) - b_{k-1}^2 p_{k-2}(\lambda)$라는 수열로 정의되고, 최종적으로 고유값은 $D(\lambda) = p_n(\lambda) = 0$으로 구한다. ($p_k(\lambda)$는 왼쪽 위 k x k 행렬의 특성방정식, $p_0(\lambda) = 1$)

위 수열의 근을

$$\begin{align*} p_{k-1}(\mu) &= 0, \mu_1 < \dots < \mu_i < \dots < \mu_{k-1} \\ p_k(\lambda) &= 0, \lambda_1 < \dots < \lambda_i < \dots < \lambda_k \end{align*}$$

이라고 하고, $\mu_0 = -\infty, \mu_k = \infty$라고 하면 항상 $\mu_{i-1} < \lambda_i < \mu_i$가 성립한다.
따라서 $p_k(\lambda) = 0$의 근을 이분법을 연속적으로 적용하여 구할 수 있다.
참고로 거쉬고린 정리를 사용해서 이분법의 범위를 줄일 수 있다.

여담으로 위 방식대로 구하면 행렬식도 $\det A = p_n(0)$으로 구할 수 있다.

고유벡터를 구할 때는 먼저 $A^*$의 고유벡터를 구한 뒤 $A$의 고유벡터를 구한다.
$A^* - \lambda I$를 전개하면 고유벡터는 아래 식으로 구할 수 있다.

$$x_1 = 1, x_i = \frac{(-1)^{i-1} p_{i-1}(\lambda)}{b_2b_3 \cdots b_i}$$

이 벡터 $x$는 $A^* = PAP$의 고유벡터이므로, $A$의 고유벡터는 $Px$로 구하면 된다.