선형계를 최적화 문제로 바꾸기
일반적으로, 매우 복잡한 다변수 함수의 해를 찾는 것은 거의 불가능에 가깝습니다.
다변수 뉴턴법(Newton's method) 등을 써서 근삿값을 구한다 하더라도 여전히
힘들죠. 하지만, 다변수 함수의 극솟값을 구하는 최적화 문제를 푸는 것은
상대적으로 쉽습니다. 기울기를 구해서 함수값이 감소하는 방향으로 계속 가다보면
충분한 횟수 이후에는 극솟값으로 수렴하겠죠.
그렇다면 선형계(linear system)를 최적화 문제로 바꿔서 풀면 어떨까요? 예를
들어, 선형계
\begin{equation*} \begin{bmatrix} 4 & -1 \\ -1 & 2 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} = \begin{bmatrix} 1 \\ 5\end{bmatrix} \end{equation*}
의 해 $(x_1, x_2)=(1, 3)$은 다음 함수의 극소점에 해당합니다.
\begin{equation*} f(x_1, x_2) = 2x_1^2-x_1x_2+x_2^2-x_1-5x_2 \end{equation*}
일반적으로 다음이 성립합니다.
대칭인 양의 정부호(positive definite) 행렬 $A$에 대해 선형계 $A\x=\b$의 해는
아래 함수 $f(\x)$의 최소점이고, 이때 $\nabla f=A\x-\b$가 성립한다.
\[ f(\x) = \frac{1}{2}\x^TA\x-\b^T\x \]
먼저 최소점임을 증명합시다. $A\x=\b$의 해를 $\x^*$라
하고 $f$에 $\x=\x^*+\h$를 대입합니다. 여기서 $\h$는 영이 아닌 임의의
벡터입니다.
\begin{align*}
f(\x^*+\h) &= \frac{1}{2}(\x^{*T}+\h^T)A(\x^*+\h)-\b^T(\x^*+\h) \\
&=\frac{1}{2}\x^{*T}A\x^*+\frac{1}{2}\x^{*T}A\h+\frac{1}{2}\h^TA\x^*+\frac{1}{2}\h^TA\h-\b^T\x^*-\b^T\h \\
&=\left(\frac{1}{2}\x^{*T}A\x^*-\b^T\x^*\right)+\h^T\cancel{(A\x^*-\b)}+\frac{1}{2}\h^TA\h \\
&=f(\x^*)+\frac{1}{2}\h^TA\h \\
&>f(\x^*)
\end{align*}
따라서 $\x^*$는 $f$의 최소점입니다. 증명 과정에서 다음 성질을 이용했습니다.
\begin{gather*}
\u^T\v=\v^T\u \\
\u^TA\v=(A\v)^T\u=\v^TA\u
\end{gather*}
그 다음으로 $\nabla f$를 계산해야 합니다. $f$를 $A_{ij}$, $x_i$, $b_i$로
풀어 써봅시다.
\begin{equation*} f(\x) = \frac{1}{2}\sum_i\sum_jx_iA_{ij}x_j-\sum_i b_ix_i \end{equation*}
$x_k$에 대해 미분하면
\begin{equation*} (\nabla f)_k = \frac{\partial f}{\partial x_k} = \frac{1}{2}\sum_j A_{kj}x_j+\frac{1}{2}\sum_ix_iA_{ik}-b_k = \sum_iA_{ki}x_i-b_k = (A\x-\b)_k \end{equation*}
이로부터 $\nabla f = A\x-\b$도 증명되었습니다.
반복법으로 최소점 찾기
보통 함수의 최소점을 찾는다고 하면 실제 최소점 $\x^*$에 가까울 것 같은 초기값
$\x_{(0)}$를 적절히 고르고(이건 순전히 본인의 감입니다), 이를 $\x_{(1)}$,
$\x_{(2)}$, $\x_{(3)}$...으로 업데이트하면서 최종적으로 $\x^*$에 수렴시키는
반복법을 씁니다. 제일 많이 쓰는 업데이트 방법은 다음과 같이 적절한 탐색 방향
$\d$를 골라서 $\d$의 배수만큼 업데이트해주는 겁니다.
\[ \x_{(i+1)} = \x_{(i)} + \alpha_{(i)}\d_{(i)} \]
계수 $\alpha_{(i)}$는 어떻게 구할까요? 이왕이면 $f(\x_{(i+1)})$이 최소가
되도록 정하면 좋겠죠. 정리 1을 증명할 때처럼 하면
\[ f(\x_{(i+1)}) = f(\x_{(i)})+\alpha_{(i)}\d_{(i)}^T(A\x_{(i)}-\b)+\frac{1}{2}\alpha_{(i)}^2 \d_{(i)}^TA\d_{(i)} \]
이게 최소가 되려면 $\alpha_{(i)}$로 미분한 값이 0이 되어야 합니다.
\begin{gather}
\d_{(i)}^TA\d_{(i)}\alpha_{(i)} + \d_{(i)}^T(A\x_{(i)}-\b) = 0 \nonumber \\
\therefore \alpha_{(i)} = \frac{\d_{(i)}^T(\b-A\x_{(i)})}{\d_{(i)}^TA\d_{(i)}}
\end{gather}
잔차(residual)를 $\r = \b-A\x$로 정의하여 다시 쓰면
\[ \alpha_{(i)} = \frac{\d_{(i)}^T\r_{(i)}}{\d_{(i)}^TA\d_{(i)}} \]
최속 강하법
그럼 탐색 방향 $\d$는 어떻게 정해야 할까요? 가장 쉽게 생각할 수 있는 방법은
$f$가 가장 빨리 감소하는 방향, 즉 $-\nabla f$로 정하는 겁니다. 그런데
$\nabla f=A\x-\b$니까 $-\nabla f=\b-A\x$, 즉 잔차 $\r$이 됩니다. 따라서 이때
$\d_{(i)}=\r_{(i)}$입니다.
\[ \x_{(i+1)}=\x_{(i)}+\alpha_{(i)}\r_{(i)}, \qquad \alpha_{(i)} = \frac{\r_{(i)}^T\r_{(i)}}{\r_{(i)}^TA\r_{(i)}}, \qquad \r_{(i)}=\b-A\x_{(i)} \]
이 알고리즘은 최속 강하법(method of steepest descent)라고 부르는 것으로,
계산량이 아주 적으며 조금만 변형하면 비선형계에서도 쓸 수 있는 알고리즘입니다.
하지만 느린 수렴 속도 때문에 선형계에서는 굳이 쓰지 않습니다.
켤레 방향법
두 벡터 $\u$와 $\v$가 $\u^TA\v=\v^TA\u=0$을 만족할 때 $\u$와 $\v$를 서로
켤레(conjugate)라고 정의합시다. 일종의 내적을 정의하는 거죠. 그러면
그람-슈미트 과정 등의 방법을 써서 서로 켤레인 $n$차원 영이 아닌 벡터 $n$개의
집합을 만들어낼 수 있고 이 집합은 이 내적에 대해 $n$차원 직교 기저를 이룹니다.
또한 당연히 이 집합은 일차 독립 집합입니다. 이를 써서 아래의 재미있는 정리를
증명할 수 있습니다.
서로 켤레인 $n$개의 영이 아닌 벡터 $\d_{(0)}$, $\d_{(1)}$, …,
$\d_{(n-1)}$을 탐색 방향으로 쓰면 $\x_{(n)}=\x^*$이다. 즉,
$A\x_{(n)}=\b$이다.
$A\x_{(n)}$을 풀어 써봅시다.
\[ \begin{aligned}
A\x_{(n)} &= A\x_{(n-1)} + \alpha_{(n-1)}A\d_{(n-1)} \\
&= A\x_{(n-2)} + \alpha_{(n-2)}A\d_{(n-2)} + \alpha_{(n-1)}A\d_{(n-1)} \\
&=\cdots \\
&= A\x_{(0)} + \alpha_{(0)}A\d_{(0)} + \alpha_{(1)}A\d_{(1)} + \cdots + \alpha_{(n-1)}A\d_{(n-1)}
\end{aligned} \label{eq:Axn} \]
양변에서 $\b$를 빼고 잔차로 바꾸면
\begin{equation*} -\r_{(n)} = -\r_{(0)} + \alpha_{(0)}A\d_{(0)} + \alpha_{(1)}A\d_{(1)} + \cdots + \alpha_{(n-1)}A\d_{(n-1)} \end{equation*}
양변의 왼쪽에 $\d_{(i)}^T$를 곱합니다.
\begin{equation*} -\d_{(i)}^T\r_{(n)} = -\d_{(i)}^T\r_{(0)} + \alpha_{(0)}\d_{(i)}^TA\d_{(0)} + \alpha_{(1)}\d_{(i)}^TA\d_{(1)} + \cdots + \alpha_{(n-1)}\d_{(i)}^TA\d_{(n-1)} \end{equation*}
각 $\d$는 서로 켤레라고 했으므로 켤레의 정의에 의해서
$\d_{(i)}^TA\d_{(j)}$는 $i\neq j$일 때 0이 됩니다. 따라서,
\begin{equation*} -\d_{(i)}^T\r_{(n)} = -\d_{(i)}^T\r_{(0)} + \alpha_{(i)}\d_{(i)}^TA\d_{(i)} \end{equation*}
$\alpha_{(i)}$의 정의를 대입합시다.
\[ -\d_{(i)}^T\r_{(n)} = -\d_{(i)}^T\r_{(0)} + \d_{(i)}^T\r_{(i)} \label{eq:diTrn} \]
그런데 식 \eqref{eq:Axn}처럼 $\r_{(i)}$를 풀어쓰면
\begin{align*}
\r_{(i)} &= \b-A\x_{(i)} \\
&= \b-A\x_{(0)}-\alpha_{(0)}A\d_{(0)}-\alpha_{(1)}A\d_{(1)}-\cdots-\alpha_{(i-1)}A\d_{(i-1)} \\
&= \r_{(0)}-\alpha_{(0)}A\d_{(0)}-\alpha_{(1)}A\d_{(1)}-\cdots-\alpha_{(i-1)}A\d_{(i-1)}
\end{align*}
여기에도 양변에 $\d_{(i)}^T$를 곱합니다.
\begin{equation*} \d_{(i)}^T\r_{(i)} = \d_{(i)}^T\r_{(0)}-\alpha_{(0)}\d_{(i)}^TA\d_{(0)}-\alpha_{(1)}\d_{(i)}^TA\d_{(1)}-\cdots-\alpha_{(i-1)}\d_{(i)}^TA\d_{(i-1)}=\d_{(i)}^T\r_{(0)} \end{equation*}
그러므로 식 \eqref{eq:diTrn}은
\begin{equation*} \d_{(i)}^T\r_{(n)} = 0 \end{equation*}
$\d$의 집합은 일차 독립 집합이라고 했죠? 정의에 의해 $\r_{(n)}$은 다음과
같이 유일하게 표현할 수 있습니다.
\begin{equation*} \r_{(n)} = \sum_ic_i\d_{(i)} \end{equation*}
그렇다면
\begin{equation*} \r_{(n)}^T\r_{(n)} = \sum_i c_i\d_{(i)}^T\r_{(n)} = 0 \end{equation*}
이고 결과적으로 $\r_{(n)}=\mathbf{0}$, 또는 $A\x_{(n)}=\b$입니다.
이렇게 서로 켤레인 $n$개의 영이 아닌 벡터들을 탐색 방향으로 쓰는 알고리즘을
켤레 방향법(conjugate direction method)라고 합니다.
아래 행렬 $A$는 대칭이고 양의 정부호 행렬입니다.
\begin{equation*}
A = \begin{bmatrix}
3 & -1 & 2 \\
-1 & 7 & 0 \\
2 & 0 & 5
\end{bmatrix}
\end{equation*}
이 행렬 $A$에 대하여 다음 세 벡터는 서로 켤레입니다. 켤레의 정의대로 직접 곱해서 확인하면 됩니다.
\begin{equation*}
\d_{(0)}=\begin{bmatrix}1\\0\\0\end{bmatrix},\quad
\d_{(1)}=\begin{bmatrix}1\\1\\-1\end{bmatrix},\quad
\d_{(2)}=\begin{bmatrix}-1\\1\\2\end{bmatrix}
\end{equation*}
이제 위 세 벡터를 탐색 방향으로 써서 선형계
\begin{equation*}
\begin{bmatrix}
3 & -1 & 2 \\
-1 & 7 & 0 \\
2 & 0 & 5
\end{bmatrix} \x = \b = \begin{bmatrix}7\\3\\-2\end{bmatrix}
\end{equation*}
의 해를 구해봅시다. 초기값은 $\x_{(0)}=\mathbf{0}$으로 둡니다. 먼저
첫 번째 단계에서는
\begin{equation*}
\r_{(0)} = \b-A\x_{(0)} = \begin{bmatrix}7\\3\\-2\end{bmatrix}, \quad
\alpha_{(0)} = \frac{\d_{(0)}^T\r_{(0)}}{\d_{(0)}^TA\d_{(0)}} = \frac{7}{3}, \quad
\x_{(1)} = \x_{(0)}+\alpha_{(0)}\d_{(0)} = \begin{bmatrix}\frac{7}{3}\\0\\0\end{bmatrix}
\end{equation*}
두 번째는
\begin{equation*}
\r_{(1)} = \b-A\x_{(1)} = \begin{bmatrix}0\\\frac{16}{3}\\-\frac{20}{3}\end{bmatrix}, \quad
\alpha_{(1)} = \frac{\d_{(1)}^T\r_{(1)}}{\d_{(1)}^TA\d_{(1)}} = \frac{4}{3}, \quad
\x_{(2)} = \x_{(1)}+\alpha_{(1)}\d_{(1)} = \begin{bmatrix}\frac{11}{3}\\\frac{4}{3}\\-\frac{4}{3}\end{bmatrix}
\end{equation*}
마지막으로
\begin{equation*}
\r_{(2)} = \b-A\x_{(2)} = \begin{bmatrix}0\\-\frac{8}{3}\\-\frac{8}{3}\end{bmatrix}, \quad
\alpha_{(2)} = \frac{\d_{(2)}^T\r_{(2)}}{\d_{(2)}^TA\d_{(2)}} = -\frac{1}{3}, \quad
\x_{(3)} = \x_{(2)}+\alpha_{(2)}\d_{(2)} = \begin{bmatrix}4\\1\\-2\end{bmatrix}
\end{equation*}
확인해보면 $\x_{(3)}$는 해가 맞습니다.
그람-슈미트 과정으로 탐색 방향 계산하기
서로 켤레인 벡터 집합은 어떻게 구하냐고요? 위에서 말한 대로 내적
$\langle\u,\v\rangle=\u^TA\v$에 대한 직교 기저를 그람-슈미트 과정을 써서
구하면 됩니다. 일차 독립 집합 $\{\u_0, \u_1, \cdots, \u_{n-1}\}$이 주어지면
이 내적에 대한 그람-슈미트 과정은 다음과 같습니다.
\[
\begin{aligned}
\d_{(0)} &= \u_0 \\
\d_{(1)} &= \u_1-\frac{\d_{(0)}^TA\u_1}{\d_{(0)}^TA\d_{(0)}}\d_{(0)} \\
\d_{(2)} &= \u_2-\frac{\d_{(0)}^TA\u_2}{\d_{(0)}^TA\d_{(0)}}\d_{(0)}-\frac{\d_{(1)}^TA\u_2}{\d_{(1)}^TA\d_{(1)}}\d_{(1)} \\
&\ \vdots \\
\d_{(n-1)} &= \u_{n-1}-\sum_{k=0}^{n-2} \frac{\d_{(k)}^TA\u_{n-1}}{\d_{(k)}^TA\d_{(k)}}\d_{(k)}
\end{aligned}
\label{eq:gram-schmidt}
\]
예제 1의 행렬 $A$에 대해 그람-슈미트 과정으로 직교 기저를 계산해봅시다.
일차 독립 집합 $\{\u_0, \u_1, \u_2\}$는 $\mathbb{R}^3$의 표준 기저로
두겠습니다.
\begin{align*}
\d_{(0)} &= \u_0 = \begin{bmatrix}1\\0\\0\end{bmatrix} \\
\d_{(1)} &= \u_1-\frac{\d_{(0)}^TA\u_1}{\d_{(0)}^TA\d_{(0)}}\d_{(0)} = \begin{bmatrix}0\\1\\0\end{bmatrix}+\frac{1}{3}\begin{bmatrix}1\\0\\0\end{bmatrix}=\begin{bmatrix}\frac{1}{3}\\1\\0\end{bmatrix} \\
\d_{(2)} &= \u_2-\frac{\d_{(0)}^TA\u_2}{\d_{(0)}^TA\d_{(0)}}\d_{(0)}-\frac{\d_{(1)}^TA\u_2}{\d_{(1)}^TA\d_{(1)}}\d_{(1)} = \begin{bmatrix}0\\0\\1\end{bmatrix}-\frac{2}{3}\begin{bmatrix}1\\0\\0\end{bmatrix}-\frac{1}{10}\begin{bmatrix}\frac{1}{3}\\1\\0\end{bmatrix} = \begin{bmatrix}-\frac{7}{10}\\-\frac{1}{10}\\1\end{bmatrix}
\end{align*}
이 세 벡터가 서로 켤레인 것은 쉽게 검증할 수 있습니다.
하지만 이렇게 그람-슈미트 과정으로 탐색 방향을 계산하는 건 효율적이지 않습니다. 뒤쪽 벡터로 갈수록 계산량이 크게 늘어나기 때문이죠.
켤레 기울기법
최속 강하법은 계산량이 적지만 수렴 속도가 느리고, 켤레 방향법은 $n$번 이내에
수렴함이 보장되지만 미리 켤레인 탐색 방향을 구해놓아야 시작할 수 있습니다.
켤레 기울기법(conjugate gradient method, CG)은 이 둘의 절충으로, 최속
강하법처럼 탐색 방향을 잔차로 넣어주는 대신 켤레 방향법의 그람-슈미트 과정에서
필요한 일차 독립 집합을 각 단계에서의 잔차로 넣어줍니다. 즉
$\u_i=\r_{(i)}$이고,
\[ \d_{(i)}=\r_{(i)}-\sum_{k=0}^{i-1}\frac{\d_{(k)}A\r_{(i)}}{\d_{(k)}^TA\d_{(k)}}\d_{(k)} \label{eq:conjgrad-d} \]
이게 왜 좋은 선택일까요? 먼저 정리 하나를 더 보고 가겠습니다.
켤레 방향법에서 $i<j$이면 $\d_{(i)}^T\r_{(j)}=0$이다.
정리 2에서 $\d_{(i)}^T\r_{(n)}=0$을 증명한 것과 똑같이 하면 됩니다.
\begin{align*}
\r_{(j)} &= \r_{(0)}-\alpha_{(0)}A\d_{(0)}-\cdots-\alpha_{(j-1)}A\d_{(j-1)} \\
\therefore\d_{(i)}^T\r_{(j)} &= \d_{(i)}^T\r_{(0)}-\alpha_{(0)}\d_{(i)}A\d_{(0)}-\cdots-\alpha_{(j-1)}\d_{(i)}^TA\d_{(j-1)} \\
&= \d_{(i)}^T\r_{(0)}-\alpha_{(i)}\d_{(i)}^TA\d_{(i)} \\
&= \d_{(i)}^T\r_{(0)}-\d_{(i)}^T\r_{(i)} \\
&= 0
\end{align*}
예제 1에서 나온 벡터들을 대입해보면 잘 맞습니다.
정리 3에 따라 켤레 기울기법에서 잔차는 서로 직교함을 증명할 수 있습니다.
켤레 기울기법에서 $i<j$이면 $\r_{(i)}^T\r_{(j)}=0$이다.
식 \eqref{eq:conjgrad-d}을 전치하고 양변에 $\r_{(j)}$를 곱합시다. 여기서
$i<j$입니다.
\begin{equation*} \d_{(i)}^T\r_{(j)} = \r_{(i)}^T\r_{(j)}-\sum_{k=0}^{i-1}\frac{\d_{(k)}^TA\r_{(i)}}{\d_{(k)}^TA\d_{(k)}}\d_{(k)}^T\r_{(j)} \end{equation*}
$k<i<j$이므로 정리 3에 의해 $\d_{(i)}^T\r_{(j)}$와
$\d_{(k)}^T\r_{(j)}$는 모두 0이 되어 위 식에서
$\r_{(i)}^T\r_{(j)}=0$입니다. 덧붙여, 만약 $i=j$이면 우변 두 번째 항만
0이 되므로 이 때에는 $\d_{(i)}^T\r_{(i)}=\r_{(i)}^T\r_{(i)}$입니다.
잔차의 집합 $\{\r_{(0)},\cdots,\r_{(n-1)}\}$에 그람-슈미트 과정을 적용해서
켤레인 방향들의 집합을 얻을 수 있는 것도 정리 4 덕분입니다. 이 집합이 직교
집합이라서 선형 독립이기 때문입니다. 물론 이건 잔차 중에 영벡터가 없다는
보장이 있어야 하긴 한데, 애당초 잔차가 $\mathbf{0}$이면 해를 구한 거니 거기서
끝내면 됩니다.
이제 정리 4에 의해,
\begin{align*}
\r_{(i+1)} &= \r_{(i)}-\alpha_{(i)}A\d_{(i)} \\
\Rightarrow \ \r_{(i+1)}^T\r_{(j)} &= \r_{(i)}^T\r_{(j)}-\alpha_{(i)}\d_{(i)}^TA\r_{(j)} \\
\therefore \ \d_{(i)}^TA\r_{(j)} &= \frac{1}{\alpha_{(i)}}\left[\r_{(i)}^T\r_{(j)}-\r_{(i+1)}^T\r_{(j)}\right] \\
&= \begin{cases}
\displaystyle -\frac{\r_{(j)}^T\r_{(j)}}{\alpha_{(j-1)}}, & \text{for } i=j-1 \\
0, & \text{for } i < j-1
\end{cases}
\end{align*}
식 \eqref{eq:gram-schmidt}에 대입하면 합에서 $k<i-1$일 때가 모두
소거되므로 다음과 같이 아주 간단한 식이 됩니다.
\begin{align*}
\d_{(i)}&=\r_{(i)}+\frac{\r_{(i)}^T\r_{(i)}}{\alpha_{(i-1)}\d_{(i-1)}^TA\d_{(i-1)}}\d_{(i-1)}\\
&=\r_{(i)}+\frac{\r_{(i)}^T\r_{(i)}}{\d_{(i-1)}^T\r_{(i-1)}}\d_{(i-1)}\\
&=\r_{(i)}+\frac{\r_{(i)}^T\r_{(i)}}{\r_{(i-1)}^T\r_{(i-1)}}\d_{(i-1)}
\end{align*}
최종적으로 켤레 기울기법은
\[
\begin{gathered}
\d_{(0)}=\r_{(0)}=\b-A\x_{(0)}, \\
\x_{(i+1)}=\x_{(i)}+\alpha_{(i)}\d_{(i)}, \qquad
\alpha_{(i)}=\frac{\r_{(i)}^T\r_{(i)}}{\d_{(i)}^TA\d_{(i)}}, \\
\d_{(i+1)}=\r_{(i+1)}+\frac{\r_{(i+1)}^T\r_{(i+1)}}{\r_{(i)}^T\r_{(i)}}\d_{(i)}, \qquad
\r_{(i+1)}=\r_{(i)}-\alpha_{(i)}A\d_{(i)}
\end{gathered}
\]
잔차를 정의대로 계산하지 않고 이전 단계 잔차를 이용해서 업데이트하는데, 이는
$A\x_{(i)}$를 굳이 계산하지 않고 $A\d_{(i)}$를 또 사용해서 계산량을 줄이기
위해서입니다.