켤레기울기법

"오늘의AI위키"는 AI 기술로 일관성 있고 체계적인 최신 지식을 제공하는 혁신 플랫폼입니다.
"오늘의AI위키"의 AI를 통해 더욱 풍부하고 폭넓은 지식 경험을 누리세요.

1. 개요
2. 방법의 설명
3. 직접법으로서의 켤레 기울기법
4. 반복법으로서의 켤레 기울기법
- 4.1. 알고리즘
  - 4.1.1. 알파와 베타의 계산
5. 전처리 켤레 기울기법
6. 유연한 전처리 켤레 기울기법
7. 정규 방정식에 대한 켤레 기울기법 (CGNR)
8. 수렴성
9. 복소수 에르미트 행렬에 대한 켤레 기울기법
10. 장단점
참조

1. 개요

켤레 기울기법은 대칭 양의 정부호 행렬을 계수로 갖는 선형 연립방정식을 풀기 위한 반복법으로, 직접법으로 유도되었으나 반복법으로도 활용된다. 이 방법은 해를 구하기 위해 A-켤레성을 만족하는 탐색 방향 벡터를 구성하며, 특히 n이 매우 커서 직접법으로 풀기 어려운 문제에 유용하다. 켤레 기울기법은 수렴 속도가 행렬의 조건수에 영향을 받으므로 전처리를 통해 속도를 향상시킬 수 있으며, 유연한 전처리 켤레 기울기법과 CGNR(Conjugate Gradient on Normal Equations) 등의 변형된 방법이 존재한다. 이론적으로는 유한한 반복 후에 정확한 해를 생성하지만, 실제로는 반올림 오차로 인해 근사해를 구하며, 대규모 희소 행렬 시스템에 효율적인 해법이지만, 행렬의 조건수가 크거나 전처리기가 부적절하면 수렴이 느려질 수 있다. 켤레 기울기법은 복소수 에르미트 행렬에도 적용 가능하다.

더 읽어볼만한 페이지

최적화 알고리즘 - 기댓값 최대화 알고리즘
최적화 알고리즘 - 유전 알고리즘
유전 알고리즘은 생물 진화 과정을 모방하여 최적해를 탐색하는 계산 최적화 기법으로, 해를 유전자로 표현하고 적합도 함수를 통해 평가하며 선택, 교차, 돌연변이 연산을 반복하여 최적해에 근접하는 해를 생성하며, 성능은 매개변수에 영향을 받고 초기 수렴 문제 해결을 위한 다양한 변형 기법과 관련 기법이 존재한다.
수치선형대수학 - 가우스 소거법
가우스 소거법은 연립 일차 방정식의 해를 구하기 위해 행렬을 사다리꼴로 변환하는 알고리즘이며, 기본 행 연산을 통해 전진 소거와 후퇴 대입 단계를 거쳐 해를 계산하고, 행렬식 계산 등 다양한 분야에 응용된다.
수치선형대수학 - LINPACK
LINPACK은 부동소수점 연산 성능을 평가하는 벤치마크 프로그램이자 FORTRAN 라이브러리로, 슈퍼컴퓨터 성능 측정 기준으로 사용되는 HPLinpack 벤치마크의 기반이 되었으며, TOP500 목록에서 고성능 컴퓨터 순위를 결정하는 데 기여한다.

켤레기울기법
개요
종류	반복적인 방법
문제 유형	선형 방정식 Ax = b 풀기, 2차 함수 최소화
발견자	매그너스 헤스테네스, 에두아르트 슈티펠
개발 시기	1952년
설명
conjugate gradient method (CG method)	공액 기울기법 (共役傾斜法)
biconjugate gradient method (BiCG method)	쌍공액 기울기법 (双共役勾配法)
풀고자 하는 문제	실수 대칭인 정부호 행렬 A에 대한 선형 방정식 Ax = b의 해를 구하는 문제 2차 함수의 최소값을 찾는 문제
특징
장점	A가 희소 행렬인 경우 계산 효율성이 높음 추가적인 저장 공간 필요량이 적음
단점	A가 정부호 행렬이어야 함 전처리가 중요함
관련 알고리즘
선형 공액 기울기법	Ax = b 꼴의 선형 시스템에 적용
비선형 공액 기울기법	비선형 최적화 문제에 적용
참고 문헌
논문	Hestenes, Magnus R.; Stiefel, Eduard (December 1952). Methods of Conjugate Gradients for Solving Linear Systems. Journal of Research of the National Bureau of Standards. 49 (6): 409. doi:10.6028/jres.049.044
학위 논문	Straeter, T. A. (1971). On the Extension of the Davidon–Broyden Class of Rank One, Quasi-Newton Minimization Methods to an Infinite Dimensional Hilbert Space with Applications to Optimal Control Problems. NASA Technical Reports Server. North Carolina State University.
서적	Polyak, Boris (1987). Introduction to Optimization Greenbaum, Anne (1997). Iterative Methods for Solving Linear Systems Ambros Speiser, Konrad Zuse und die ERMETH: Ein weltweiter Architektur-Vergleich. In Hans Dieter Hellige (ed.). Geschichten der Informatik. Visionen, Paradigmen, Leitmotive. Berlin: Springer. 2004. ISBN 3-540-00217-0 야마모토 테츠로(山本哲朗) (2003-06). 수치 해석 입문(数値解析入門) (증정판(増訂版) ed.). 사이언스사(ja). 사이언스 라이브러리 현대 수학 입문 14. ISBN 4-7819-1038-6.

2. 방법의 설명

다음과 같은 선형 연립 방정식을 푸는 것을 목표로 한다.

:'''Ax''' = '''b'''

여기서 '''A'''는 알려진 ''n'' × ''n'' 대칭 행렬(즉, '''A'''^T = '''A''')이고, 양의 정부호 행렬(즉, 0이 아닌 모든 벡터 '''x''' ∈ '''R'''^''n''에 대해 '''x'''^T'''Ax''' > 0)이며, 실수 행렬이다. '''b'''는 알려진 ''n'' × 1 크기의 실수 열벡터이고, '''x'''는 구하고자 하는 ''n'' × 1 크기의 미지 실수 열벡터이다.

이 연립 방정식의 유일한 해를 $\mathbf{x}_*$ 라고 표기한다.

두 개의 0이 아닌 벡터 '''u'''와 '''v'''가 다음 조건을 만족할 때, 행렬 '''A'''에 대해 켤레(conjugate)라고 정의한다.

: $\mathbf{u}^\mathrm{T} \mathbf{A} \mathbf{v} = 0$

행렬 '''A'''는 대칭이고 양의 정부호이므로, 다음과 같은 내적을 정의할 수 있다.

: $\langle \mathbf{u},\mathbf{v} \rangle_\mathbf{A} := \langle \mathbf{A} \mathbf{u}, \mathbf{v}\rangle = \langle \mathbf{u}, \mathbf{A}^\mathrm{T} \mathbf{v}\rangle = \langle \mathbf{u}, \mathbf{A}\mathbf{v} \rangle = \mathbf{u}^\mathrm{T} \mathbf{A} \mathbf{v}$

따라서 두 벡터가 이 내적 $\langle \cdot, \cdot \rangle_\mathbf{A}$ 에 대해 직교하면, 두 벡터는 '''A'''에 대해 켤레 관계에 있다.

켤레 관계는 대칭적이다. 즉, '''u'''가 '''v'''에 대해 켤레이면, '''v'''도 '''u'''에 대해 켤레이다.

이제 $P = \{\mathbf{p}_1, \mathbf{p}_2, \dots, \mathbf{p}_n\}$ 가 ''n''개의 서로 켤레인 방향 벡터들의 집합이라고 가정하자. 즉, 모든 $i \neq k$ 에 대해 $\langle\mathbf{p}_i,\mathbf{p}_k\rangle_{A} = \mathbf{p}_i^\mathrm{T} \mathbf{A} \mathbf{p}_k = 0$ 이다. 이러한 벡터 집합 $P$ 는 $\mathbb{R}^n$ 공간의 기저를 이룬다.

따라서, 연립 방정식의 해 $\mathbf{x}_*$ 는 기저 벡터 $\mathbf{p}_i$ 들의 선형 결합으로 다음과 같이 유일하게 표현될 수 있다.

: $\mathbf{x}_* = \sum^{n}_{i=1} \alpha_i \mathbf{p}_i$

여기서 계수 $\alpha_i$ 를 구하기 위해, 원래 방정식 '''Ax''' = '''b'''에 위 식을 대입한다.

: $\mathbf{A}\mathbf{x}_* = \mathbf{A} \left( \sum^{n}_{i=1} \alpha_i \mathbf{p}_i \right) = \sum^{n}_{i=1} \alpha_i (\mathbf{A} \mathbf{p}_i) = \mathbf{b}$

이제 이 식의 양변에 $\mathbf{p}_k^\mathrm{T}$ (단, $1 \le k \le n$ )를 왼쪽에서 곱한다.

: $\mathbf{p}_k^\mathrm{T} \left( \sum^{n}_{i=1} \alpha_i \mathbf{A} \mathbf{p}_i \right) = \sum^{n}_{i=1} \alpha_i \left( \mathbf{p}_k^\mathrm{T} \mathbf{A} \mathbf{p}_i \right) = \mathbf{p}_k^\mathrm{T} \mathbf{b}$

벡터 $\mathbf{p}_i$ 들이 서로 켤레이므로, $i \neq k$ 일 때는 $\mathbf{p}_k^\mathrm{T} \mathbf{A} \mathbf{p}_i = \langle \mathbf{p}_k, \mathbf{p}_i \rangle_\mathbf{A} = 0$ 이다. 따라서 합계 기호 안에서는 $i = k$ 인 항만 남게 된다.

: $\alpha_k (\mathbf{p}_k^\mathrm{T} \mathbf{A} \mathbf{p}_k) = \mathbf{p}_k^\mathrm{T} \mathbf{b}$

$\mathbf{p}_k$ 는 0벡터가 아니고 '''A'''는 양의 정부호 행렬이므로 $\mathbf{p}_k^\mathrm{T} \mathbf{A} \mathbf{p}_k = \langle \mathbf{p}_k, \mathbf{p}_k \rangle_\mathbf{A} = \|\mathbf{p}_k\|_\mathbf{A}^2 > 0$ 이다. 따라서 계수 $\alpha_k$ 는 다음과 같이 계산된다.

: $\alpha_k = \frac{\mathbf{p}_k^\mathrm{T} \mathbf{b}}{\mathbf{p}_k^\mathrm{T} \mathbf{A} \mathbf{p}_k} = \frac{\langle \mathbf{p}_k, \mathbf{b}\rangle}{\langle \mathbf{p}_k, \mathbf{p}_k\rangle_\mathbf{A}} = \frac{\langle \mathbf{p}_k, \mathbf{b}\rangle}{\|\mathbf{p}_k\|_\mathbf{A}^2}$

이 결과는 위에서 정의한 내적 $\langle \cdot, \cdot \rangle_\mathbf{A}$ 을 고려할 때 직관적으로 이해될 수 있다. 즉, 해 $\mathbf{x}_*$ 의 $\mathbf{p}_k$ 방향 성분은 $\mathbf{b}$ 를 $\mathbf{p}_k$ 방향으로 투영한 것과 관련된다.

결론적으로, '''Ax''' = '''b''' 형태의 선형 연립 방정식을 푸는 한 가지 방법은 다음과 같다.

# ''n''개의 서로 켤레인 방향 벡터 $\mathbf{p}_1, \dots, \mathbf{p}_n$ 을 찾는다.

# 각 방향 벡터 $\mathbf{p}_k$ 에 대해 계수 $\alpha_k = \frac{\mathbf{p}_k^\mathrm{T} \mathbf{b}}{\mathbf{p}_k^\mathrm{T} \mathbf{A} \mathbf{p}_k}$ 를 계산한다.

# 해 $\mathbf{x}_* = \sum_{k=1}^{n} \alpha_k \mathbf{p}_k$ 를 구한다.

켤레 기울기법은 이러한 켤레 방향 벡터들을 효율적으로 생성하고 동시에 해를 점진적으로 개선하는 방법을 제공한다.

3. 직접법으로서의 켤레 기울기법

켤레 기울기법은 원래 선형 연립 방정식 '''Ax''' = '''b'''를 풀기 위한 직접법으로 고안되었다. 이 방법의 핵심 아이디어는 행렬 '''A'''에 대해 서로 켤레(conjugate) 또는 '''A'''-직교(orthogonal)인 ''n''개의 탐색 방향 벡터 '''p'''_''k''를 구성하여 해 '''x'''_*를 찾는 것이다.^[4]^[21]^[22]^[23]

두 개의 영벡터가 아닌 벡터 '''u'''와 '''v'''가 다음 조건을 만족할 때, 이 두 벡터는 행렬 '''A'''에 대해 켤레라고 한다.

: $\mathbf{u}^\mathrm{T} \mathbf{A} \mathbf{v} = 0$

행렬 '''A'''가 대칭이고 양의 정부호이기 때문에, 다음과 같이 내적을 정의할 수 있다.

: $\langle \mathbf{u},\mathbf{v} \rangle_\mathbf{A} := \langle \mathbf{A} \mathbf{u}, \mathbf{v}\rangle = \langle \mathbf{u}, \mathbf{A}^\mathrm{T} \mathbf{v}\rangle = \langle \mathbf{u}, \mathbf{A}\mathbf{v} \rangle = \mathbf{u}^\mathrm{T} \mathbf{A} \mathbf{v}$

따라서 두 벡터 '''u'''와 '''v'''가 '''A'''-켤레라는 것은 위에서 정의한 내적 `< · , · >_'''A'''`에 대해 서로 직교한다는 의미와 같다. 켤레 관계는 대칭적이어서, '''u'''가 '''v'''에 대해 켤레이면 '''v'''도 '''u'''에 대해 켤레이다.

이제 ''n''개의 벡터로 이루어진 집합 $P = \{\mathbf{p}_1, \dots, \mathbf{p}_n\}$ 가 있고, 이 벡터들이 서로 '''A'''-켤레라고 가정하자. 즉, 모든 $i \neq k$ 에 대해 $\langle\mathbf{p}_i, \mathbf{p}_k\rangle_{A} = \mathbf{p}_i^\mathrm{T} \mathbf{A} \mathbf{p}_k = 0$ 이 성립한다. 이 ''n''개의 상호 켤레인 벡터들의 집합 $P$ 는 $\mathbb{R}^n$ 공간의 기저를 형성한다. 따라서 방정식 '''Ax''' = '''b'''의 유일한 해 '''x'''_*는 이 기저 벡터들의 선형 결합으로 다음과 같이 표현될 수 있다.

: $\mathbf{x}_* = \sum^{n}_{i=1} \alpha_i \mathbf{p}_i$

이 식의 양변에 행렬 '''A'''를 곱하면 다음과 같다.

: $\mathbf{A}\mathbf{x}_* = \sum^{n}_{i=1} \alpha_i \mathbf{A} \mathbf{p}_i = \mathbf{b}$

이제 위 식의 양변에 특정 기저 벡터 '''p'''_k^T를 왼쪽에 곱하면,

: $\mathbf{p}_k^\mathrm{T} \mathbf{A}\mathbf{x}_* = \sum^{n}_{i=1} \alpha_i\mathbf{p}_k^\mathrm{T} \mathbf{A} \mathbf{p}_i = \mathbf{p}_k^\mathrm{T} \mathbf{b}$

벡터 '''p'''_i와 '''p'''_k는 $i \neq k$ 일 때 서로 '''A'''-켤레이므로, $\mathbf{p}_k^\mathrm{T} \mathbf{A} \mathbf{p}_i = \langle \mathbf{p}_k, \mathbf{p}_i \rangle_\mathbf{A} = 0$ 이다. 따라서 위 합계에서 $i = k$ 인 항만 남게 된다.

: $\alpha_k \mathbf{p}_k^\mathrm{T} \mathbf{A} \mathbf{p}_k = \mathbf{p}_k^\mathrm{T} \mathbf{b}$

이를 통해 계수 α_k를 다음과 같이 구할 수 있다.

: $\alpha_k = \frac{\mathbf{p}_k^\mathrm{T} \mathbf{b}}{\mathbf{p}_k^\mathrm{T} \mathbf{A} \mathbf{p}_k} = \frac{\langle \mathbf{p}_k, \mathbf{b}\rangle}{\langle \mathbf{p}_k, \mathbf{p}_k\rangle_\mathbf{A}} = \frac{\langle \mathbf{p}_k, \mathbf{b}\rangle}{\|\mathbf{p}_k\|_\mathbf{A}^2}$

여기서 $\langle \cdot, \cdot \rangle$ 는 표준 유클리드 내적을 나타낸다.

결론적으로, 선형 시스템 '''Ax''' = '''b'''를 풀기 위해 먼저 ''n''개의 서로 '''A'''-켤레인 방향 벡터 '''p'''_k들을 찾고, 그 다음 위의 공식을 사용하여 계수 α_k들을 계산하면 해 '''x'''_*를 구할 수 있다. 이것이 켤레 기울기법을 직접법으로 이해하는 방식이다.^[4]

4. 반복법으로서의 켤레 기울기법

선형 연립 방정식 $\mathbf{Ax} = \mathbf{b}$ 에서 행렬 $\mathbf{A}$ 의 크기 ''n''이 매우 크면 직접법으로 해를 구하는 데 많은 시간이 소요될 수 있다. 이런 경우 켤레 기울기법은 해 $\mathbf{x}_*$ 를 효율적으로 근사하는 반복법으로 사용될 수 있다.^[21]^[22]^[23] 켤레 벡터 $\mathbf{p}_k$ 들을 적절히 선택하면, 모든 ''n''개의 벡터를 사용하지 않고도 해에 충분히 가까운 근사값을 얻을 수 있기 때문이다.

반복법으로서 켤레 기울기법은 초기 추측값 $\mathbf{x}_0$ (특별한 정보가 없다면 보통 $\mathbf{x}_0 = \mathbf{0}$ 으로 설정)에서 시작하여, 반복적인 계산을 통해 실제 해 $\mathbf{x}_*$ 에 점진적으로 근접해 나간다. 각 반복 단계에서 해에 얼마나 가까워졌는지를 판단하는 기준이 필요한데, 이는 해 $\mathbf{x}_*$ 가 다음 이차 함수 $f(\mathbf{x})$ 를 유일하게 최소화하는 벡터라는 사실을 이용한다.

: $f(\mathbf{x}) = \frac12 \mathbf{x}^\mathrm{T} \mathbf{A}\mathbf{x} - \mathbf{x}^\mathrm{T} \mathbf{b} , \quad \mathbf{x}\in\mathbf{R}^n.$

행렬 $\mathbf{A}$ 가 대칭 양의 정부호 행렬이므로 이 함수는 유일한 최소점을 가지며, 그 최소점에서 함수의 기울기는 0이 된다. 즉, $\nabla f(\mathbf{x}_*) = \mathbf{A} \mathbf{x}_* - \mathbf{b} = \mathbf{0}$ 이므로, $f(\mathbf{x})$ 를 최소화하는 것은 원래의 선형 연립 방정식을 푸는 것과 같다. 따라서 반복 과정에서 $f(\mathbf{x}_k)$ 의 값이 계속 감소하면 해 $\mathbf{x}_*$ 에 가까워지고 있다고 볼 수 있다.

''k''번째 반복 단계에서의 잔차는 다음과 같이 정의된다.

: $\mathbf{r}_k = \mathbf{b} - \mathbf{Ax}_k.$

이 잔차 $\mathbf{r}_k$ 는 현재 위치 $\mathbf{x}_k$ 에서 함수 $f(\mathbf{x})$ 의 음의 기울기( $-\nabla f(\mathbf{x}_k)$ )와 같다. 경사 하강법에서는 이 기울기의 반대 방향, 즉 잔차 $\mathbf{r}_k$ 방향으로 해를 업데이트하지만, 켤레 기울기법은 다른 전략을 사용한다.

첫 번째 탐색 방향 $\mathbf{p}_0$ 는 초기 잔차 $\mathbf{r}_0 = \mathbf{b} - \mathbf{Ax}_0$ 로 설정한다. 이후의 탐색 방향 $\mathbf{p}_k$ ( $k \ge 1$ )는 현재 잔차 $\mathbf{r}_k$ 와 이전 탐색 방향들을 이용하여, 이전의 모든 탐색 방향 $\mathbf{p}_0, \dots, \mathbf{p}_{k-1}$ 들과 '''A'''-켤레(A-orthogonal) 관계, 즉 $\mathbf{p}_k^\mathrm{T} \mathbf{A} \mathbf{p}_i = 0$ ( $i < k$ )를 만족하도록 결정된다. 이 때문에 "켤레 기울기법"이라는 이름이 붙었다. 이 과정은 그람-슈미트 과정과 유사하게, 잔차 $\mathbf{r}_k$ 에서 이전 탐색 방향 $\mathbf{p}_i$ ( $i < k$ )들의 성분을 '''A'''-켤레 조건에 맞게 제거하여 새로운 탐색 방향 $\mathbf{p}_k$ 를 만드는 것으로 볼 수 있다. 이론적으로는 다음과 같이 표현된다.

: $\mathbf{p}_{k} = \mathbf{r}_{k} - \sum_{i < k}\frac{\mathbf{p}_i^\mathrm{T} \mathbf{A} \mathbf{r}_{k}}{\mathbf{p}_i^\mathrm{T}\mathbf{A} \mathbf{p}_i} \mathbf{p}_i$

새로운 탐색 방향 $\mathbf{p}_k$ 가 결정되면, 이 방향으로 얼마나 이동할지를 결정하는 스텝 크기 $\alpha_k$ 를 계산하여 해를 업데이트한다.

: $\mathbf{x}_{k+1} = \mathbf{x}_k + \alpha_k \mathbf{p}_k$

스텝 크기 $\alpha_k$ 는 $\mathbf{x}_{k+1}$ 에서 함수 $f(\mathbf{x})$ 의 값을 최소화하도록, 즉 $\mathbf{p}_k$ 방향으로의 직선 탐색을 통해 결정된다. 이는 다음 식으로 계산된다.

: $\alpha_{k} = \frac{\mathbf{p}_k^\mathrm{T} \mathbf{r}_k}{\mathbf{p}_k^\mathrm{T} \mathbf{A} \mathbf{p}_k}$

위의 $\mathbf{p}_k$ 계산식만 보면 모든 이전 탐색 방향 벡터들을 저장하고 계산해야 하므로 비효율적으로 보일 수 있다. 그러나 켤레 기울기법의 중요한 성질인 잔차 벡터들의 직교성( $\mathbf{r}_i^\mathrm{T} \mathbf{r}_j = 0, i \neq j$ )과 탐색 방향 벡터들의 '''A'''-켤레성( $\mathbf{p}_i^\mathrm{T} \mathbf{A} \mathbf{p}_j = 0, i \neq j$ ) 덕분에, 실제 알고리즘에서는 단지 이전 단계의 $\mathbf{x}_k, \mathbf{r}_k, \mathbf{p}_k$ 값만을 이용하여 $\mathbf{x}_{k+1}, \mathbf{r}_{k+1}, \mathbf{p}_{k+1}$ 을 효율적으로 계산할 수 있다. 이 과정은 해 $\mathbf{x}_k$ 를 크릴로프 부분 공간 $\mathcal{K}_k(\mathbf{A}, \mathbf{r}_0) = \operatorname{span}\{ \mathbf{r}_0, \mathbf{Ar}_0, \dots, \mathbf{A}^{k-1}\mathbf{r}_0 \}$ 내에서 찾는 과정과 밀접하게 관련되어 있다.^[6] 구체적인 계산 알고리즘과 $\alpha_k$ , $\beta_k$ 의 유도는 하위 섹션 '알고리즘'에서 자세히 설명한다.

4. 1. 알고리즘

선형 연립 방정식

\mathbf{Ax} = \mathbf{b}

의 해

\mathbf{x}_*

를 구하기 위한 반복법으로서의 켤레 기울기법은 다음과 같이 정리될 수 있다. 이 방법은 특히 행렬

\mathbf{A}

가 크고 희소 행렬일 때 효과적이다. 알고리즘은 실수 값을 가지며, 대칭이고 양의 정부호 행렬인

\mathbf{A}

에 대해 정의된다.

초기 추측값

\mathbf{x}_0

가 주어졌을 때 (보통

\mathbf{x}_0 = \mathbf{0}

으로 설정), 알고리즘은 다음 단계를 반복적으로 수행한다.

1. 초기 잔차(residual)와 검색 방향(search direction)을 계산한다:

$\mathbf{r}_0 := \mathbf{b} - \mathbf{A x}_0$
$\mathbf{p}_0 := \mathbf{r}_0$

2. 반복 시작 (

k = 0, 1, 2, \dots

):

스텝 크기(step size) $\alpha_k$ 를 계산한다. 이는 현재 검색 방향 $\mathbf{p}_k$ 로 얼마나 이동하여 잔차를 최소화할지 결정한다.

\alpha_k := \frac{\mathbf{r}_k^\mathrm{T} \mathbf{r}_k}{\mathbf{p}_k^\mathrm{T} \mathbf{A p}_k}

해의 근사값 $\mathbf{x}_k$ 를 업데이트한다.

\mathbf{x}_{k+1} := \mathbf{x}_k + \alpha_k \mathbf{p}_k

잔차 $\mathbf{r}_k$ 를 업데이트한다.

\mathbf{r}_{k+1} := \mathbf{r}_k - \alpha_k \mathbf{A p}_k

종료 조건을 확인한다. 잔차 $\mathbf{r}_{k+1}$ 의 노름 $||\mathbf{r}_{k+1}||$ 이 충분히 작아지면 (미리 정한 허용 오차보다 작아지면) 반복을 종료하고 $\mathbf{x}_{k+1}$ 를 해로 반환한다.

\hbox{if } ||\mathbf{r}_{k+1}|| \hbox{ is sufficiently small then exit loop}

다음 검색 방향을 만들기 위한 계수 $\beta_k$ 를 계산한다. 이는 새로운 검색 방향이 이전 검색 방향과 '''A'''-켤레 관계를 가지도록 보장한다.

\beta_k := \frac{\mathbf{r}_{k+1}^\mathrm{T} \mathbf{r}_{k+1}}{\mathbf{r}_k^\mathrm{T} \mathbf{r}_k}

다음 검색 방향 $\mathbf{p}_{k+1}$ 을 계산한다. 이는 새로운 잔차와 이전 검색 방향의 조합으로 만들어지며, 이전의 모든 검색 방향과 '''A'''-켤레 관계를 유지한다.

\mathbf{p}_{k+1} := \mathbf{r}_{k+1} + \beta_k \mathbf{p}_k

반복 횟수 $k$ 를 1 증가시킨다. ( $k := k + 1$ )

이 알고리즘은 각 반복 단계에서 이전 단계의 값들(

\mathbf{x}_k, \mathbf{r}_k, \mathbf{p}_k

)만을 필요로 하므로 저장 공간 측면에서 효율적이다. 또한, 각 반복에서 주로 요구되는 계산은 행렬-벡터 곱셈(

\mathbf{Ap}_k

) 한 번과 몇 번의 벡터 내적 및 벡터 덧셈/스칼라 곱셈이므로 계산 비용도 상대적으로 낮다.^[4]^[23]

\alpha_k

와

\beta_k

의 구체적인 유도 과정은 하위 섹션에서 자세히 다룬다. 이 알고리즘의

\beta_k

계산 방식은 비선형 최적화 문제에 적용되는 플레처-리브스(Fletcher–Reeves) 켤레 기울기법에서도 사용된다.

4. 1. 1. 알파와 베타의 계산

알고리즘에서 스칼라 값

\alpha_k

와

\beta_k

는 특정 조건을 만족하도록 계산된다.

'''

\alpha_k

의 계산'''

\alpha_k

는 새로운 잔차 벡터

\mathbf{r}_{k+1}

가 이전 잔차 벡터

\mathbf{r}_k

에 직교하도록 선택된다. 즉, 두 벡터의 내적이 0 (

\mathbf{r}_{k+1}^\mathrm{T} \mathbf{r}_k = 0

)이 되도록 한다.

알고리즘의 잔차 업데이트 식

\mathbf{r}_{k+1} = \mathbf{r}_k - \alpha_k \mathbf{A} \mathbf{p}_k

의 양변에

\mathbf{r}_k^\mathrm{T}

를 곱하면 다음과 같다.

\mathbf{r}_k^\mathrm{T} \mathbf{r}_{k+1} = \mathbf{r}_k^\mathrm{T} \mathbf{r}_k - \alpha_k \mathbf{r}_k^\mathrm{T} \mathbf{A p}_k

좌변이 0이 되어야 하므로, 위 식을

\alpha_k

에 대해 정리하면 다음과 같다.

\alpha_k = \frac{\mathbf{r}_{k}^\mathrm{T} \mathbf{r}_{k}}{\mathbf{r}_{k}^\mathrm{T} \mathbf{A p}_k}

여기서 검색 방향

\mathbf{p}_k

는

\mathbf{p}_k = \mathbf{r}_k + \beta_{k-1} \mathbf{p}_{k-1}

로 정의되고, 켤레기울기법의 성질에 의해

\mathbf{p}_{k-1}

과

\mathbf{r}_k

는 직교(

\mathbf{p}_{k-1}^\mathrm{T} \mathbf{r}_k = 0

)하며,

\mathbf{p}_{k-1}

과

\mathbf{p}_k

는 행렬 '''A'''에 대해 켤레(A-직교,

\mathbf{p}_{k-1}^\mathrm{T} \mathbf{A p}_k = 0

) 관계이다. 이 성질들을 이용하면 분모를 다음과 같이 변형할 수 있다.

\mathbf{p}_k^\mathrm{T} \mathbf{A p}_k = (\mathbf{r}_k + \beta_{k-1} \mathbf{p}_{k-1})^\mathrm{T} \mathbf{A p}_k = \mathbf{r}_k^\mathrm{T} \mathbf{A p}_k + \beta_{k-1} \underbrace{\mathbf{p}_{k-1}^\mathrm{T} \mathbf{A p}_k}_{0} = \mathbf{r}_k^\mathrm{T} \mathbf{A p}_k

따라서 최종적으로 알고리즘에서 사용되는

\alpha_k

식은 다음과 같다.

\alpha_k = \frac{\mathbf{r}_k^\mathrm{T} \mathbf{r}_k}{\mathbf{p}_k^\mathrm{T} \mathbf{A p}_k}

'''

\beta_k

의 계산'''

\beta_k

는 새로운 검색 방향

\mathbf{p}_{k+1}

가 이전 검색 방향

\mathbf{p}_k

에 대해 행렬 '''A'''에 대해 켤레(A-직교)가 되도록 선택된다. 즉,

\mathbf{p}_{k+1}^\mathrm{T} \mathbf{A} \mathbf{p}_k = 0

이 되도록 한다.

알고리즘의 검색 방향 업데이트 식

\mathbf{p}_{k+1} = \mathbf{r}_{k+1} + \beta_k \mathbf{p}_k

의 양변에

\mathbf{p}_k^\mathrm{T} \mathbf{A}

를 곱하면 다음과 같다.

\mathbf{p}_k^\mathrm{T} \mathbf{A} \mathbf{p}_{k+1} = \mathbf{p}_k^\mathrm{T} \mathbf{A} \mathbf{r}_{k+1} + \beta_k \mathbf{p}_k^\mathrm{T} \mathbf{A} \mathbf{p}_k

좌변이 0이 되어야 하므로, 위 식을

\beta_k

에 대해 정리하면 다음과 같다. (대칭행렬 '''A'''의 성질(

\mathbf{A}^T = \mathbf{A}

)에 의해

\mathbf{p}_k^\mathrm{T} \mathbf{A} \mathbf{r}_{k+1} = (\mathbf{A} \mathbf{p}_k)^\mathrm{T} \mathbf{r}_{k+1} = \mathbf{r}_{k+1}^\mathrm{T} \mathbf{A} \mathbf{p}_k

이다.)

\beta_k = - \frac{\mathbf{r}_{k+1}^\mathrm{T} A \mathbf{p}_k}{\mathbf{p}_k^\mathrm{T} A \mathbf{p}_k}

분자를 변형하기 위해 잔차 업데이트 식

\mathbf{r}_{k+1} = \mathbf{r}_k - \alpha_k \mathbf{A p}_k

에서 유도되는

\mathbf{A p}_k = \frac{1}{\alpha_k} (\mathbf{r}_k - \mathbf{r}_{k+1})

관계를 이용한다.

\mathbf{r}_{k+1}^\mathrm{T} A \mathbf{p}_k = \mathbf{r}_{k+1}^\mathrm{T} \frac{1}{\alpha_k} (\mathbf{r}_k - \mathbf{r}_{k+1}) = \frac{1}{\alpha_k} (\underbrace{\mathbf{r}_{k+1}^\mathrm{T} \mathbf{r}_k}_{0} - \mathbf{r}_{k+1}^\mathrm{T} \mathbf{r}_{k+1})

앞서

\alpha_k

를 정의할 때

\mathbf{r}_{k+1}

과

\mathbf{r}_k

는 직교하도록 했으므로 (

\mathbf{r}_{k+1}^\mathrm{T} \mathbf{r}_k = 0

), 분자는 다음과 같이 간단해진다.

\mathbf{r}_{k+1}^\mathrm{T} A \mathbf{p}_k = - \frac{1}{\alpha_k} \mathbf{r}_{k+1}^\mathrm{T} \mathbf{r}_{k+1}

분모는 위에서

\alpha_k

를 계산할 때 유도된 관계

\mathbf{p}_k^\mathrm{T} \mathbf{A p}_k = \mathbf{r}_k^\mathrm{T} \mathbf{A p}_k

와

\mathbf{A p}_k = \frac{1}{\alpha_k} (\mathbf{r}_k - \mathbf{r}_{k+1})

를 이용한다.

\mathbf{p}_k^\mathrm{T} A \mathbf{p}_k = \mathbf{r}_k^\mathrm{T} \mathbf{A p}_k = \mathbf{r}_k^\mathrm{T} \frac{1}{\alpha_k} (\mathbf{r}_k - \mathbf{r}_{k+1}) = \frac{1}{\alpha_k} (\mathbf{r}_k^\mathrm{T} \mathbf{r}_k - \underbrace{\mathbf{r}_k^\mathrm{T} \mathbf{r}_{k+1}}_{0})

마찬가지로

\mathbf{r}_k

와

\mathbf{r}_{k+1}

는 직교하므로 (

\mathbf{r}_k^\mathrm{T} \mathbf{r}_{k+1} = 0

), 분모는 다음과 같이 간단해진다.

\mathbf{p}_k^\mathrm{T} A \mathbf{p}_k = \frac{1}{\alpha_k} \mathbf{r}_k^\mathrm{T} \mathbf{r}_k

이제

\beta_k

식에 변형된 분자와 분모를 대입하면

\alpha_k

가 소거되고 최종적인

\beta_k

식을 얻는다.

\beta_k = - \frac{- \frac{1}{\alpha_k} \mathbf{r}_{k+1}^\mathrm{T} \mathbf{r}_{k+1}}{\frac{1}{\alpha_k} \mathbf{r}_k^\mathrm{T} \mathbf{r}_k} = \frac{\mathbf{r}_{k+1}^\mathrm{T} \mathbf{r}_{k+1}}{\mathbf{r}_k^\mathrm{T} \mathbf{r}_k}

이 계산 방식은 비선형 최적화 방법 중 하나인 플레처-리브스(Fletcher–Reeves) 방법에서도 사용된다.

5. 전처리 켤레 기울기법

켤레기울기법의 수렴 속도는 계수 행렬 '''A'''의 조건수에 크게 의존한다. 조건수가 클수록, 즉 행렬이 특이행렬에 가까울수록 수렴 속도는 느려진다. 전처리는 기존의 선형 시스템 '''Ax''' = '''b'''를 동일한 해를 가지면서도 조건수가 더 작은, 즉 더 빨리 수렴하는 새로운 시스템으로 변환하는 기법이다.^[10] 대부분의 경우, 켤레기울기법의 빠른 수렴을 위해서는 적절한 전처리가 필수적이다.

전처리 켤레 기울기법은 전처리 행렬 '''M'''을 도입하여 적용한다. 이때 '''M'''은 대칭 양의 정부호 행렬이어야 하며, '''M'''^-1'''A'''의 조건수가 원래 행렬 '''A'''의 조건수보다 작아야 한다. 전처리 켤레 기울기법의 알고리즘은 다음과 같다:^[10]

: $\mathbf{r}_0 := \mathbf{b} - \mathbf{A x}_0$

: $\text{Solve } \mathbf{M}\mathbf{z}_0 := \mathbf{r}_0$ (즉, $\mathbf{z}_0 = \mathbf{M}^{-1}\mathbf{r}_0$ 계산)

: $\mathbf{p}_0 := \mathbf{z}_0$

: $k := 0 \,$

:'''반복'''

:: $\alpha_k := \frac{\mathbf{r}_k^\mathsf{T} \mathbf{z}_k}{\mathbf{p}_k^\mathsf{T} \mathbf{A p}_k}$ (스텝 길이 계산)

:: $\mathbf{x}_{k+1} := \mathbf{x}_k + \alpha_k \mathbf{p}_k$ (근사 해 업데이트)

:: $\mathbf{r}_{k+1} := \mathbf{r}_k - \alpha_k \mathbf{A p}_k$ (잔차 업데이트)

::'''만약''' '''r'''_''k''+1 이 충분히 작으면 (수렴 판정) '''반복 종료''' '''만약'''

:: $\text{Solve } \mathbf{M}\mathbf{z}_{k+1} := \mathbf{r}_{k+1}$ (전처리된 잔차 계산)

:: $\beta_k := \frac{\mathbf{r}_{k+1}^\mathsf{T} \mathbf{z}_{k+1}}{\mathbf{r}_k^\mathsf{T} \mathbf{z}_k}$ (개선 방향 계수 계산)

:: $\mathbf{p}_{k+1} := \mathbf{z}_{k+1} + \beta_k \mathbf{p}_k$ (새로운 검색 방향 계산)

:: $k := k + 1 \,$

:'''반복 종료'''

:결과는 '''x'''_''k''+1이다.

이 알고리즘은 다음과 같이 변형된 시스템에 일반적인 켤레 기울기법을 적용한 것과 수학적으로 동일하다.^[11]

: $\mathbf{E}^{-1}\mathbf{A}(\mathbf{E}^{-1})^\mathsf{T}\mathbf{\hat{x}}=\mathbf{E}^{-1}\mathbf{b}$

여기서 $\mathbf{M} = \mathbf{EE}^\mathsf{T}$ ('''E'''는 '''M'''의 촐레스키 분해 인자 등)이고 $\mathbf{\hat{x}}=\mathbf{E}^\mathsf{T}\mathbf{x}$ 이다. 이 변환은 시스템의 대칭성과 양의 정부호성을 유지한다. 중요한 점은 촐레스키 분해 $\mathbf{M} = \mathbf{EE}^\mathsf{T$ 를 실제로 계산할 필요 없이, 각 반복 단계에서 $\mathbf{Mz = r}$ 형태의 선형 시스템을 푸는 것으로 충분하다는 것이다. 변환된 행렬 $\mathbf{E}^{-1}\mathbf{A}(\mathbf{E}^{-1})^\mathsf{T}$ 는 $\mathbf{M}^{-1}\mathbf{A}$ 와 동일한 스펙트럼(고유값 분포)을 가지므로, 조건수가 개선되는 효과를 얻는다.

전처리 행렬 '''M'''은 반드시 대칭 양의 정부호 행렬이어야 하며, 반복 과정 동안 변하지 않아야 한다. 만약 이 조건들이 만족되지 않으면 알고리즘의 수렴성이 보장되지 않거나 예측 불가능한 동작을 보일 수 있다.

다양한 종류의 전처리자가 존재한다.

불완전 촐레스키 분해: 가장 널리 사용되는 전처리 방법 중 하나이다.^[12] 행렬 '''A'''의 촐레스키 분해를 근사적으로 계산하여 '''M'''을 구성한다.
야코비 전처리 또는 대각 스케일링: 가장 간단한 형태로, '''M'''을 '''A'''의 대각 성분만으로 이루어진 대각 행렬로 선택한다. 구현이 매우 쉽고 메모리 사용량이 적지만, 성능 향상은 제한적일 수 있다.
가우스-자이델법 기반 전처리: 가우스-자이델 반복법의 아이디어를 활용하여 전처리 행렬을 구성한다.
대칭 SOR법 기반 전처리: SOR 방법을 대칭적으로 적용하여 전처리 행렬을 만든다.

효과적인 전처리 행렬 '''P'''(여기서는 '''M'''과 동일한 역할)는 원래 시스템 '''Ax'''='''b'''를 푸는 것보다 변환된 시스템 '''P'''^-1'''A''' ('''P'''^T)^-1 '''x'''′ ='''P'''^-1'''b'''′ (여기서 '''P''' = '''L''' '''L'''^T 형태이고 '''x'''′ = '''P'''^T'''x''')를 푸는 것이 더 쉽도록, 즉 변환된 시스템의 조건수가 원래 시스템보다 작도록 선택된다.^[23] 전처리 행렬을 선택할 때는 조건수 개선 효과와 전처리 시스템(

\mathbf{Mz = r}

)을 푸는 데 드는 계산 비용 사이의 균형을 고려해야 한다.

알고리즘 구현 시, 각 반복 단계에서

\mathbf{Mz = r}

를 풀어야 한다. 이때 전처리 행렬 '''M'''의 역행렬

\mathbf{M}^{-1}

을 직접 계산하는 것은 매우 비효율적이며, 종종 켤레 기울기법 자체보다 더 많은 계산 시간을 요구한다. 대신, '''M'''의 구조를 이용하여 이 시스템을 효율적으로 풀어야 한다.

예를 들어, 불완전 촐레스키 분해를 통해

\mathbf{M} = \mathbf{LL}^\mathsf{T}

('''L'''은 하삼각 행렬) 형태의 전처리 행렬을 얻었다면,

\mathbf{Mz = r}

즉

\mathbf{LL}^\mathsf{T}\mathbf{z} = \mathbf{r}

는 다음 두 단계로 나누어 풀 수 있다.

1.

\mathbf{La = r}

: 전진 대입법을 사용하여 중간 벡터 '''a'''를 구한다. '''L'''이 하삼각 행렬이므로 이 계산은 빠르다.

2.

\mathbf{L}^\mathsf{T}\mathbf{z = a}

: 후진 대입법을 사용하여 최종 벡터 '''z'''를 구한다. '''L'''^T는 상삼각 행렬이므로 이 계산 역시 빠르다.

이 방식을 사용하면 '''M'''이나 '''L'''의 역행렬을 명시적으로 계산할 필요 없이, 두 번의 삼각 시스템 풀이를 통해 필요한 벡터 '''z'''를 효율적으로 얻을 수 있다.

6. 유연한 전처리 켤레 기울기법

수치적으로 어려운 문제를 풀 때 정교한 전처리기가 필요할 수 있는데, 이는 반복 과정에서 전처리 행렬이 바뀌는 가변 전처리 상황을 만들 수 있다. 전처리기가 매 반복마다 대칭 양의 정부호(SPD)이더라도, 전처리기가 바뀔 수 있다는 점 때문에 기존 켤레기울기법의 수렴 속도가 실제 계산에서는 상당히 느려질 수 있다.

이러한 가변 전처리를 허용하는 방법으로 유연한(flexible) 전처리 켤레 기울기법이 있다.^[14] 이 방법은 기존의 Fletcher–Reeves 공식

: $\beta_k := \frac{\mathbf{r}_{k+1}^\mathsf{T} \mathbf{z}_{k+1}}{\mathbf{r}_k^\mathsf{T} \mathbf{z}_k}$

대신 Polak–Ribière 공식

: $\beta_k := \frac{\mathbf{r}_{k+1}^\mathsf{T} \left(\mathbf{z}_{k+1}-\mathbf{z}_{k}\right)}{\mathbf{r}_k^\mathsf{T} \mathbf{z}_k}$

을 사용하여 β_k를 계산한다. 이 공식을 사용하면 가변 전처리 상황에서 수렴 속도가 크게 향상될 수 있다.^[13]

유연한 버전은 전처리기가 SPD가 아닌 경우에도 안정적으로 작동하는 것으로 알려져 있다.^[15] 다만, 유연한 버전을 구현하려면 추가적인 벡터를 저장해야 한다.

만약 고정된 SPD 전처리기를 사용한다면, 정확한 산술 연산(즉, 반올림 오차가 없는 경우)에서는 $\mathbf{r}_{k+1}^\mathsf{T} \mathbf{z}_{k}=0$ 이므로 Fletcher–Reeves 공식과 Polak–Ribière 공식은 동일한 값을 계산한다.

Polak–Ribière 공식을 사용했을 때 수렴이 더 잘 되는 이유는 이 방법이 해당 상황에서 국소적으로 최적(locally optimal)이기 때문이며, 특히 국소적으로 최적인 최급강하법보다 느리게 수렴하지 않는다는 수학적 설명이 있다.^[16]

7. 정규 방정식에 대한 켤레 기울기법 (CGNR)

켤레 기울기법은 임의의 ''n'' × ''m'' 행렬 '''A'''에도 적용할 수 있다. 이는 정규 방정식 '''A'''^T'''Ax''' = '''A'''^T'''b'''에 켤레 기울기법을 적용하는 방식으로 이루어지는데, 모든 '''A'''에 대해 '''A'''^T'''A'''가 대칭 양의 준정부호 행렬이기 때문이다.^[36] 이 방법을 정규 방정식에 대한 켤레 기울기법(Conjugate Gradient on Normal equations, CGNR)이라고 부른다.

:'''A'''^T'''Ax''' = '''A'''^T'''b'''

CGNR은 반복적인 방법이므로, 행렬 '''A'''^T'''A'''를 메모리에 명시적으로 저장할 필요 없이 행렬-벡터 곱셈과 전치 행렬-벡터 곱셈만 수행하면 된다. 따라서 '''A'''가 희소 행렬일 때 특히 유용하며, 이러한 연산은 일반적으로 매우 효율적이다.

그러나 정규 방정식을 사용하는 것의 단점은 조건수 κ('''A'''^T'''A''')가 κ('''A''')²와 같다는 점이다. 이 때문에 CGNR의 수렴 속도가 느려질 수 있으며, 근사 해의 품질이 반올림 오류에 민감해질 수 있다. 따라서 좋은 전처리자를 찾는 것이 CGNR 방법을 효과적으로 사용하는 데 중요한 요소가 된다.

이러한 단점을 보완하기 위해 여러 알고리즘이 제안되었는데, 대표적으로 CGLS(Conjugate Gradient Least Squares)^[37]와 LSQR이 있다. 특히 [https://web.stanford.edu/group/SOL/software/lsqr/ LSQR] 알고리즘은 행렬 '''A'''의 조건이 나쁠 때, 즉 조건수가 클 때 가장 수치적으로 안정적인 해법으로 알려져 있다.^[38]^[39]

8. 수렴성

켤레 기울기법은 이론적으로 직접적인 방법으로 볼 수 있는데, 이는 반올림 오차가 없을 경우 행렬의 크기( $n$ )보다 크지 않은 유한한 수의 반복 후에 정확한 해 $\mathbf{x}_*$ 를 생성하기 때문이다. 그러나 실제 계산에서는 반올림 오차의 영향과 작은 섭동에도 불안정한 특성 때문에 유한 번 반복으로 정확한 해를 얻는 경우는 드물다. 예를 들어, 계산 과정에서 생성되는 방향 벡터들이 이론과 달리 실제로는 크릴로프 부분 공간을 제대로 생성하지 못하고 서로 켤레(conjugate) 관계를 잃기 쉽다.

반복법으로서 켤레 기울기법은 반복 횟수 $k$ 가 증가함에 따라 근사해 $\mathbf{x}_{k}$ 가 실제 해 $\mathbf{x}_*$ 에 단조적으로 수렴하는 특성을 가진다. 즉, 오차 $\mathbf{e}_k = \mathbf{x}_k - \mathbf{x}_*$ 의 크기는 에너지 노름 $\|\mathbf{e}_k\|_\mathbf{A} = \sqrt{\mathbf{e}_k^\mathrm{T} \mathbf{A} \mathbf{e}_k}$ 기준으로 계속 감소한다. 수렴 속도는 시스템 행렬 $\mathbf{A}$ 의 조건수 $\kappa(\mathbf{A}) = \|\mathbf{A}\| \|\mathbf{A}^{-1}\|$ (보통 $\|\cdot\|_2$ 노름 사용 시 $\kappa(\mathbf{A}) = \lambda_{\max}/\lambda_{\min}$ )에 크게 의존하며, 조건수가 클수록 수렴은 느려진다.^[8]

만약 $\kappa(\mathbf{A})$ 가 크다면, 수렴 속도를 높이기 위해 전처리(preconditioning) 기법을 사용한다. 이는 원래 선형 시스템 $\mathbf{A x} = \mathbf{b}$ 를 푸는 대신, 이와 동등하면서도 조건수가 더 작은 새로운 시스템 $\mathbf{M}^{-1}\mathbf{A x} = \mathbf{M}^{-1}\mathbf{b}$ 를 푸는 방식이다. 여기서 $\mathbf{M}$ 은 전처리기(preconditioner)라 불리는 행렬로, $\kappa(\mathbf{M}^{-1}\mathbf{A}) < \kappa(\mathbf{A})$ 가 되도록 선택한다.

수렴 속도를 더 자세히 분석하기 위해, 최대 차수가 $k$ 이고 $p(0)=1$ 인 다항식의 집합을 $\Pi_k^* := \left\lbrace p \in \Pi_k \ : \ p(0)=1 \right\rbrace$ 로 정의하자. 켤레 기울기법의 $k$ 번째 반복에서의 오차 $\mathbf{e}_k$ 는 에너지 노름 $\|\cdot\|_\mathbf{A}$ 에 대해 다음 부등식을 만족한다.^[4]^[9]

: $\begin{align}\left\| \mathbf{e}_k \right\|_\mathbf{A}&= \min_{p \in \Pi_k^*} \left\| p(\mathbf{A}) \mathbf{e}_0 \right\|_\mathbf{A}\\&\leq \min_{p \in \Pi_k^*} \left( \max_{\lambda \in \sigma(\mathbf{A})} |p(\lambda)| \right) \left\| \mathbf{e}_0 \right\|_\mathbf{A}\\&\leq 2 \left( \frac{ \sqrt{\kappa(\mathbf{A})}-1 }{ \sqrt{\kappa(\mathbf{A})}+1 } \right)^k \left\| \mathbf{e}_0 \right\|_\mathbf{A}\\&\leq 2 \exp\left(\frac{-2k}{\sqrt{\kappa(\mathbf{A})}}\right) \left\| \mathbf{e}_0 \right\|_\mathbf{A}\end{align}$

여기서 $\mathbf{e}_0 = \mathbf{x}_0 - \mathbf{x}_*$ 는 초기 오차이고, $\sigma(\mathbf{A})$ 는 행렬 $\mathbf{A}$ 의 고유값 스펙트럼을 나타낸다.

이 부등식은 오차 $\|\mathbf{e}_k\|_\mathbf{A}$ 가 반복 횟수 $k$ 에 따라 지수적으로 감소하며, 그 감소율이 조건수 $\kappa(\mathbf{A})$ 의 제곱근에 반비례함을 보여준다. 특히, 초기 오차를 $2\varepsilon$ ( $\varepsilon > 0$ ) 수준으로 줄이기 위해서는 대략 $k \approx \frac{1}{2}\sqrt{\kappa(\mathbf{A})} \log\left(\frac{\|\mathbf{e}_0\|_\mathbf{A}}{\varepsilon}\right)$ 번의 반복이 필요함을 알 수 있다.

조건수 $\kappa(\mathbf{A})$ 가 매우 클 때 ( $\kappa(\mathbf{A}) \gg 1$ ), 수렴 인자 $\frac{ \sqrt{\kappa(\mathbf{A})}-1 }{ \sqrt{\kappa(\mathbf{A})}+1 }$ 는 $1 - \frac{2}{\sqrt{\kappa(\mathbf{A})}}$ 에 가깝다. 이는 자코비나 가우스-자이델과 같은 고전적인 반복법의 수렴 인자( $\approx 1 - \frac{C}{\kappa(\mathbf{A})}$ 형태)보다 1에 더 가깝지 않으므로, 켤레 기울기법이 훨씬 빠른 수렴 속도를 가짐을 의미한다.

이론적인 수렴 분석은 반올림 오차가 없다고 가정하지만, 실제 계산에서는 반올림 오차가 수렴에 영향을 미친다. 그럼에도 불구하고, 앤 그린바움(Anne Greenbaum) 등이 보였듯이 위에서 제시된 수렴 경계는 실제 계산에서도 대체로 유효한 경향을 보인다.^[5]

실제 계산에서 켤레 기울기법의 수렴은 종종 다음과 같은 단계적 양상을 보인다.^[5]

초기 단계: 처음 몇 번의 반복 동안에는 수렴이 매우 빠르다. 이는 오차 성분 중 행렬 $\mathbf{A}$ 의 큰 고유값에 해당하는 성분들이 먼저 제거되면서, 마치 작은 유효 조건수를 가진 문제처럼 동작하기 때문이다.
중간 단계: 수렴 속도가 이론적인 예측치인 $\approx 1 - \frac{2}{\sqrt{\kappa(\mathbf{A})}}$ 에 가까워진다. 때로는 행렬 $\mathbf{A}$ 의 고유값 분포나 초기 오차의 특성에 따라 이론적인 속도보다 더 빠른 초선형적(superlinear) 수렴을 보이기도 한다.
최종 단계: 반복이 계속되면 반올림 오차가 누적되어 더 이상 오차가 줄어들지 않고 정체되거나, 심한 경우 오차가 다시 증가하며 발산할 수도 있다.

따라서 배정밀도 부동 소수점 형식을 사용하는 일반적인 과학 계산에서는 보통 미리 정해진 허용 오차(tolerance)에 도달하면 초기 또는 중간 단계에서 반복을 멈추는 중지 기준을 사용한다.

9. 복소수 에르미트 행렬에 대한 켤레 기울기법

켤레 기울기법은 약간의 수정을 통해 복소수 값을 갖는 행렬 '''A'''와 벡터 '''b'''가 주어졌을 때, 선형 방정식 시스템 $\mathbf {A} \mathbf {x} =\mathbf {b}$ 를 푸는 데 사용될 수 있다. 이 방법이 적용되려면 행렬 '''A'''가 에르미트 행렬(즉, A' = A, 여기서 '는 켤레 전치를 의미)이고 양의 정부호 행렬이어야 한다. 이때 해 벡터 '''x'''도 복소수 값을 가진다. 필요한 수정은 단순히 실수 전치 연산을 켤레 전치 연산으로 대체하는 것이다.

10. 장단점

켤레 기울기법은 이론적으로 반올림 오차가 없을 경우, 행렬의 크기보다 크지 않은 유한한 수의 반복 후에 정확한 해를 생성하기 때문에 직접 해법으로 볼 수 있다. 그러나 실제로는 작은 섭동에도 불안정하여 정확한 해를 얻는 경우는 드물다. 이는 계산 과정에서 생성되는 방향 벡터들이 크릴로프 부분 공간을 생성하는 퇴행적인 성질 때문에 실제로는 켤레가 아니게 되는 경우가 많기 때문이다.
장점:

반복법으로서, 켤레 기울기법은 정확한 해에 대한 근사값 $\mathbf{x}_{k}$ 을 (에너지 노름에서) 단조적으로 개선하며, 문제 크기에 비해 비교적 적은 수의 반복 후에 필요한 공차에 도달할 수 있다.
전처리를 통해 수렴 속도를 크게 향상시킬 수 있다. 원래 시스템 $\mathbf{A x}-\mathbf{b} = 0$ 을 $\mathbf{M}^{-1}(\mathbf{A x}-\mathbf{b}) = 0$ 으로 대체하는 방식인데, 적절한 전처리기 $\mathbf{M}$ 을 사용하면 변환된 시스템의 조건수 $\kappa(\mathbf{M}^{-1}\mathbf{A})$ 가 원래 시스템의 조건수 $\kappa(\mathbf{A})$ 보다 작아져 수렴이 빨라진다.

단점:

실제 계산에서는 반올림 오차나 작은 섭동에 불안정하여, 이론적인 장점인 유한 반복 내 정확한 해 도달은 보장되지 않는다.
수렴 속도는 시스템 행렬 $A$ 의 조건수 $\kappa(A)$ 에 의해 결정되며, 일반적으로 선형적이다. 조건수가 클수록 수렴 속도는 느려진다.^[8]
조건수가 매우 크거나 부적절한 전처리 행렬을 사용하면 수렴이 매우 느려지거나 불안정해질 수 있다.

정규 방정식에 대한 켤레 기울기 (CGNR)켤레 기울기법은 정규 방정식

\mathbf{A}^\mathrm{T}\mathbf{A x} = \mathbf{A}^\mathrm{T}\mathbf{b}

에 적용하여 임의의 ''n''×''m'' 행렬

\mathbf{A}

에 대한 최소제곱 문제를 풀 수도 있다. 이 방법을 CGNR(Conjugate Gradient on Normal Residual) 또는 CGN이라고 부른다.^[36] 행렬

\mathbf{A}^\mathrm{T}\mathbf{A}

는 항상 대칭 양의 정부호 행렬(또는 준정부호 행렬)이므로 켤레 기울기법을 적용할 수 있다.

장점 (CGNR):
반복법이므로 행렬 $\mathbf{A}^\mathrm{T}\mathbf{A}$ 를 메모리에 명시적으로 구성할 필요 없이, 행렬-벡터 곱셈( $\mathbf{A}\mathbf{v}$ )과 전치 행렬-벡터 곱셈( $\mathbf{A}^\mathrm{T}\mathbf{u}$ ) 연산만 수행하면 된다.
따라서 $\mathbf{A}$ 가 희소 행렬일 때 이러한 연산은 일반적으로 매우 효율적이며, CGNR은 특히 유용하다.
단점 (CGNR):
정규 방정식을 구성하는 것의 단점은 조건수가 $\kappa(\mathbf{A}^\mathrm{T}\mathbf{A}) = \kappa(\mathbf{A})^2$ 로 제곱이 되어, 원래 문제보다 조건수가 훨씬 커진다는 점이다.
이로 인해 CGNR의 수렴 속도가 느려지는 경향이 있다.
근사 해의 품질이 반올림 오차에 민감할 수 있다.
따라서 좋은 전처리기를 찾는 것이 CGNR 방법을 사용하는 데 종종 중요한 부분이다.

CGNR의 수렴 속도 및 안정성 문제를 개선하기 위해 CGLS (Conjugate Gradient Least Squares)^[37]나 [https://web.stanford.edu/group/SOL/software/lsqr/ LSQR]과 같은 여러 알고리즘이 제안되었다. 특히 LSQR은 행렬

\mathbf{A}

가 잘 조절되지 않은 경우(ill-conditioned), 즉 조건수가 큰 경우에 가장 수치적으로 안정적인 방법으로 알려져 있다.^[38]^[39]

켤레 기울기법의 장단점에 대한 자세한 내용은 네미로프스키와 벤탈의 강의 노트에 요약되어 있다.^[18]

참조

_[1] 논문 Methods of Conjugate Gradients for Solving Linear Systems http://nvlpubs.nist.[...] 1952-12
_[2] 학위논문 On the Extension of the Davidon–Broyden Class of Rank One, Quasi-Newton Minimization Methods to an Infinite Dimensional Hilbert Space with Applications to Optimal Control Problems North Carolina State University 1971
_[3] 서적 Geschichten der Informatik. Visionen, Paradigmen, Leitmotive Springer
_[4] 서적 Introduction to Optimization https://www.research[...]
_[5] 서적 Iterative Methods for Solving Linear Systems
_[6] 논문 Universality for the Conjugate Gradient and MINRES Algorithms on Sample Covariance Matrices https://onlinelibrar[...] 2023-03
_[7] 서적 An Introduction to the Conjugate Gradient Method Without the Agonizing Pain http://www.cs.cmu.ed[...]
_[8] 서적 Iterative methods for sparse linear systems https://archive.org/[...] Society for Industrial and Applied Mathematics
_[9] 서적 Iterative solution of large sparse systems of equations Springer 2016-06-21
_[10] 서적 Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods http://www.netlib.or[...] SIAM 2020-03-31
_[11] 서적 Matrix Computations Johns Hopkins University Press 2013
_[12] 논문 Block Preconditioning for the Conjugate Gradient Method https://escholarship[...]
_[13] 논문 Inexact Preconditioned Conjugate Gradient Method with Inner-Outer Iteration
_[14] 논문 Flexible Conjugate Gradients
_[15] 논문 Nonsymmetric Preconditioning for Conjugate Gradient and Steepest Descent Methods 1
_[16] 논문 Steepest Descent and Conjugate Gradient Methods with Variable Preconditioning
_[17] arXiv An Optimal Control Theory for Accelerated Optimization 2019
_[18] 웹사이트 Optimization III: Convex Optimization http://www2.isye.gat[...] 2023
_[19] 웹사이트 Random Matrix Theory and Machine Learning Tutorial https://random-matri[...] 2023-12-05
_[20] 서적 数値解析入門サイエンス社 2003-06
_[21] 서적 数値解析第2版共立出版
_[22] 서적 数値線形代数の数理とHPC 共立出版 2018-08
_[23] 서적 UNIX & Informatioin Science-5 C 言語による数値計算入門 2005
_[24] 서적 偏微分方程式の数値解析岩波書店 2010
_[25] 서적 偏微分方程式の数値シミュレーション東京大学出版会 2003
_[26] 논문 Numerical linear algebra and solvability of partial differential equations 2002
_[27] 서적 Numerical linear algebra and optimization Addison-Wesley 1991
_[28] 논문 Global convergence properties of conjugate gradient methods for optimization 1992
_[29] 논문 The conjugate gradient method and trust regions in large scale optimization 1983
_[30] 웹사이트 Biconjugate Gradient Method https://mathworld.wo[...]
_[31] 간행물 Nonlinear conjugate gradient methods Wiley Encyclopedia of Operations Research and Management Science 2010
_[32] 논문 A survey of nonlinear conjugate gradient methods 2006
_[33] 논문 Convergence properties of nonlinear conjugate gradient methods 2000
_[34] 논문 Efficient implementation of a class of preconditioned conjugate gradient methods 1981
_[35] 논문 Preconditioned conjugate gradients for solving singular systems 1988
_[36] 웹사이트 Conjugate Gradient Method on the Normal Equations. https://mathworld.wo[...]
_[37] 서적 Numerical methods for least squares problems SIAM 1996
_[38] 논문 LSQR: An Algorithm for Sparse Linear Equations and Sparse Least Squares. 1982
_[39] 논문 Algorithm 583: LSQR: Sparse linear equations and least squares problems. 1982

본 사이트는 AI가 위키백과와 뉴스 기사,정부 간행물,학술 논문등을 바탕으로 정보를 가공하여 제공하는 백과사전형 서비스입니다.
모든 문서는 AI에 의해 자동 생성되며, CC BY-SA 4.0 라이선스에 따라 이용할 수 있습니다.
하지만, 위키백과나 뉴스 기사 자체에 오류, 부정확한 정보, 또는 가짜 뉴스가 포함될 수 있으며, AI는 이러한 내용을 완벽하게 걸러내지 못할 수 있습니다.
따라서 제공되는 정보에 일부 오류나 편향이 있을 수 있으므로, 중요한 정보는 반드시 다른 출처를 통해 교차 검증하시기 바랍니다.

문의하기 : help@durumis.com