켤레기울기법

이 문서의 내용은 출처가 분명하지 않습니다. 이 문서를 편집하여, 신뢰할 수 있는 출처를 표기해 주세요. 검증되지 않은 내용은 삭제될 수도 있습니다. 내용에 대한 의견은 토론 문서에서 나누어 주세요. (2013년 5월)

켤레기울기법 또는 공역기울기법(영어: Conjugate Gradient Method, 일본어: 共役勾配法)이란 수학에서 대칭인 양의 준정부호행렬(陽-準定符號行列, 영어: positive-semidefinite matrix)을 갖는 선형계의 해를 구하는 수치 알고리즘이다. 이는 보통 반복알고리즘으로 풀리고, 따라서 숄레스키 분해와 같은 방법이나 직접 풀기에 너무 큰 계가 갖는 희소행렬 등에 사용하기 적합하다. 그러한 계는 편미분 방정식들이나 최적화 문제들을 수치적으로 풀 때 자주 등장한다.

켤레기울기법은 또한 에너지 최소화같은 제약조건이 없는 최적화문제를 풀 수 있다. 이것은 Magnus Hestenes와 Eduard Stiefel에 의해 개발되었다.

우리는 다음과 같은 선형 계의 방정식을 풀기 원한다고 가정하자.

Ax = b

A는 n x n 대칭행렬이고 (i.e. A^T = A), 양의 정부호행렬 positive definite (i.e. Rⁿ에서 모두 0이 아닌 벡터들 x에 관하여 x^TAx > 0 )이고, 실수이고, x, b는 n x 1 실수 인 열벡터이다.

여기서 우리는 이 계의 유일한 해를 ${\displaystyle \scriptstyle \mathbf {x} _{*}}$라고 하자.

다음과 같이

${\displaystyle \mathbf {u} ^{\mathrm {T} }\mathbf {A} \mathbf {v} =0.}$ 일 때,

두개의 영이 아닌 벡터들을 u 과 v를 켤레라고 하자. (A에 대하여) A는 대칭이고 양의 정부호를 갖는 행렬이기 때문에 왼쪽 변이 내적으로 정의 된다.

${\displaystyle \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에 켤레이면 v가 u에 켤레이다.

${\displaystyle P=\{\mathbf {p} _{k}:\forall i\neq k,i,k\in [1,n],\langle \mathbf {p} _{i},\mathbf {p} _{k}\rangle _{A}=0\}}$

는 n개 상호적인 켤레 방향들의 집합이라고 가정하자. 그러면 ${\displaystyle P}$는 ${\displaystyle \mathbb {R} ^{n}}$의 기저이고, 그래서 ${\displaystyle P}$내에서 우리는 그 해 ${\displaystyle \mathbf {Ax} =\mathbf {b} }$의 ${\displaystyle \mathbf {x} _{*}}$를 확장시킬 수 있다:

${\displaystyle \mathbf {x} _{*}=\sum _{i=1}^{n}\alpha _{i}\mathbf {p} _{i}}$

그리고 우리가 보면

${\displaystyle \mathbf {b} =\mathbf {A} \mathbf {x} _{*}=\sum _{i=1}^{n}\alpha _{i}\mathbf {A} \mathbf {p} _{i}.}$

어떤 ${\displaystyle \mathbf {p} _{k}\in P}$에 관하여,

${\displaystyle \mathbf {p} _{k}^{\mathrm {T} }\mathbf {b} =\mathbf {p} _{k}^{\mathrm {T} }\mathbf {A} \mathbf {x} _{*}=\sum _{i=1}^{n}\alpha _{i}\mathbf {p} _{k}^{\mathrm {T} }\mathbf {A} \mathbf {p} _{i}=\alpha _{k}\mathbf {p} _{k}^{\mathrm {T} }\mathbf {A} \mathbf {p} _{k}.}$

(왜냐하면 ${\displaystyle \forall i\neq k,p_{i},p_{k}}$는 상호적인 켤레이기 때문에)

${\displaystyle \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}}}.}$

이 결과는 위에 정의된 내적을 고려하면서 이 결과는 아마도 매우 명백하다.

이것은 Ax = b 방정식을 풀 때 다음의 방법을 쓴다. : n개의 켤레 방향들의 수열을 찾고, 그런뒤 그 계수들 ${\displaystyle \scriptstyle \alpha _{k}}$을 계산한다.

만일 우리가 켤레 벡터들을 p_k 신중하게 선택하고, 우리는 그 해 ${\displaystyle \scriptstyle \mathbf {x} _{*}}$ 의 좋은 근사 값을 얻기 위해서 그들의 모두가 필요하지 않을 수 있다. 그래서 우리는 반복적인 방법으로서 공역 구배법을 보기를 원한다. n이 매우 커서 직접적인 방법으로 풀기에 시간이 매우 많이 걸리는 계들의 경우 이 방법을 거의 정확하게 풀게 해준다.

${\displaystyle \scriptstyle \mathbf {x} _{*}}$의 초기 추측값을 x₀라고 하자. 우리는 일반성을 잃지 않고 x₀ = 0 (마찬가지로, 대신에 그 계를 Az = b − Ax₀ ) 라고 가정할 수 있다. x₀와 함께 시작하여 우리는 그 해를 찾고, 각각의 반복에서 우리는 그 해 ${\displaystyle \scriptstyle \mathbf {x} _{*}}$(우리에게 잘알려지지 않은 해) 에 근접한게 무엇인지 알려주는 측도가 필요하다. 이 측도는 그 해 ${\displaystyle \scriptstyle \mathbf {x} _{*}}$가 또한 다음의 이차함수의 유일하게 최소치라는 사실로부터 나왔다. 그래서 만일 f(x)가 반복에 있어서 작아진다면 그것은 우리가 ${\displaystyle \scriptstyle \mathbf {x} _{*}}$에 가까이 간다는 것을 의미한다.

${\displaystyle f(\mathbf {x} )={\frac {1}{2}}\mathbf {x} ^{\mathrm {T} }\mathbf {A} \mathbf {x} -\mathbf {x} ^{\mathrm {T} }\mathbf {b} ,\quad \mathbf {x} \in \mathbf {R} ^{n}.}$

x = x₀에서 f의 구배(영어:gradient)의 음의 값이 되도록 처음 기저 p₀를 잡는다.

그 f의 구배(영어:gradient)가 Ax−b와 같고, 추측 해x₀를 가지고 시작한다. (우리는 만일 우리가 다른 어떤 것을 추측할 이유가 없다면 ${\displaystyle \scriptstyle \mathbf {x} _{*}}$ 가 0 이고 그 집합을 x₀ 0으로 항상 추측할 수 있다. ) 이것을 p₀ = b−Ax₀라고 두자. 그 기저안에 다른 벡터들은 그 구배(영어:gradient)에 켤레가 될 것이다. 그래서 이 방법의 이름이 공역 구배법(영어:conjugate gradient method)이다.

우선 k번째 단계에서 r_k를 유수(영어:residue)라고 하자: ${\displaystyle \mathbf {r} _{k}=\mathbf {b} -\mathbf {Ax} _{k}.\,}$

r_k는 x = x_k에서 음의 f의 구배(gradient)라는 점을 염두에 두고, 그래서 경사 하강법은 r_k 방향에서 움직일 것이다. 여기, 우리는 그 방향들 p_k이 각각 다른것과 켤레라고 주장한다. 이 시행은 충분히 합리적인데, 우리는 또한 다음 찾는 방향은 최근의 유수와 이전의 모든 검색방향들에 의해서 정하게 한다.

그 공역 제약조건은 직교-유형의 제약조건이고 그래서, 그 알고리즘은 그램 -슈미트와 유사하다.

이것은 다음의 표현과 같다:

${\displaystyle \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}}$

(수렴에 켤레성 제약조건의 영향에 관한 문서의 상단에 있는 그림을 참조) 이 방향에 따르면, 그 다음 최적의 위치는 다음에 의해 주어진다.

${\displaystyle \mathbf {x} _{k+1}=\mathbf {x} _{k}+\alpha _{k}\mathbf {p} _{k}}$

${\displaystyle \alpha _{k}={\frac {\mathbf {p} _{k}^{\mathrm {T} }\mathbf {b} }{\mathbf {p} _{k}^{\mathrm {T} }\mathbf {A} \mathbf {p} _{k}}}={\frac {\mathbf {p} _{k}^{\mathrm {T} }(\mathbf {r} _{k-1}+\mathbf {Ax} _{k-1})}{\mathbf {p} _{k}^{\mathrm {T} }\mathbf {A} \mathbf {p} _{k}}}={\frac {\mathbf {p} _{k}^{\mathrm {T} }\mathbf {r} _{k-1}}{\mathbf {p} _{k}^{\mathrm {T} }\mathbf {A} \mathbf {p} _{k}}},}$

를 포함함

(p_k 와x_k-1 가 켤레이기 때문에 마지막 등식이 성립한다.)

그 위에 알고리즘은 공역 구배법의 가장 간단한 설명을 제공한다. 일견에, 그 알고리즘에 명시된대로 이전의 검색 방향들과 유수 벡터들의 저장이 필요할 뿐만 아니라 많은 행렬 벡터의 곱셉들이 필요로 하고, 따라서 계산하는데 비용이 많이 들 수 있다. 그러나 그 알고리즘의 근접한 분석은 다음 과 같이 볼 수 있다.

모든 i < k에 관하여 r_k+1는 p_i과 켤레이다. (예를 들어 귀납법에 의해 증명될 수 있다) 그러므로 오직 r_k, p_k, 그리고 x_k 이 r_k+1, p_k+1, 그리고 x_k+1를 세우는데 필요로 한다. 게다가 오직 하나의 행렬-벡터 곱셈이 각각의 반복적인 계산마다 필요하다.

실수이고, 대칭행렬이고 양의 정부호 행렬인 A에 관하여 Ax = b를 푸는 알고리즘 방법은 밑에 상세히 설명되어 있다.

그 입력 벡터 x₀는 거의 정확한 초기 해나 0이 될 수 있다. 이것은 위에 명시된 정확한 절차의 다른 식이다.

${\displaystyle {\begin{aligned}&\mathbf {r} _{0}:=\mathbf {b} -\mathbf {Ax} _{0}\\&\mathbf {p} _{0}:=\mathbf {r} _{0}\\&k:=0\\&{\hbox{repeat}}\\&\qquad \alpha _{k}:={\frac {\mathbf {r} _{k}^{\mathrm {T} }\mathbf {r} _{k}}{\mathbf {p} _{k}^{\mathrm {T} }\mathbf {Ap} _{k}}}\\&\qquad \mathbf {x} _{k+1}:=\mathbf {x} _{k}+\alpha _{k}\mathbf {p} _{k}\\&\qquad \mathbf {r} _{k+1}:=\mathbf {r} _{k}-\alpha _{k}\mathbf {Ap} _{k}\\&\qquad {\hbox{if }}r_{k+1}{\hbox{ is sufficiently small then exit loop}}\\&\qquad \beta _{k}:={\frac {\mathbf {r} _{k+1}^{\mathrm {T} }\mathbf {r} _{k+1}}{\mathbf {r} _{k}^{\mathrm {T} }\mathbf {r} _{k}}}\\&\qquad \mathbf {p} _{k+1}:=\mathbf {r} _{k+1}+\beta _{k}\mathbf {p} _{k}\\&\qquad k:=k+1\\&{\hbox{end repeat}}\\&{\hbox{The result is }}\mathbf {x} _{k+1}\end{aligned}}}$

이것은 일반적으로 자주 사용되는 알고리즘이다. 이와 같은 식은${\displaystyle \beta _{k}}$은 비선형 공역 구배법인 플레처 리브스에서 또한 사용된다.

알고리즘에서, ${\displaystyle \alpha _{k}}$ 는 ${\displaystyle \mathbf {r} _{k+1}}$가 ${\displaystyle \mathbf {r} _{k}}$에 직교가 되도록 선택된다. 그 분모는 다음과 같이 단순화된다.

${\displaystyle \alpha _{k}={\frac {\mathbf {r} _{k}^{\mathrm {T} }\mathbf {r} _{k}}{\mathbf {r} _{k}^{\mathrm {T} }\mathbf {Ap} _{k}}}={\frac {\mathbf {r} _{k}^{\mathrm {T} }\mathbf {r} _{k}}{\mathbf {p} _{k}^{\mathrm {T} }\mathbf {Ap} _{k}}}}$

${\displaystyle \mathbf {r} _{k}=\mathbf {p} _{k}-\mathbf {\beta } _{k-1}\mathbf {p} _{k-1}}$때문에 그 ${\displaystyle \beta _{k}}$는 ${\displaystyle \mathbf {p} _{k+1}}$ 가 ${\displaystyle \mathbf {p} _{k}}$에 켤레가 되도록 선택된다. 처음에, ${\displaystyle \beta _{k}}$ 는

${\displaystyle \beta _{k}=-{\frac {\mathbf {r} _{k+1}^{\mathrm {T} }A\mathbf {p} _{k}}{\mathbf {p} _{k}^{\mathrm {T} }A\mathbf {p} _{k}}}}$

${\displaystyle \mathbf {r} _{k}=\mathbf {r} _{k-1}-\alpha _{k-1}A\mathbf {p} _{k-1}}$를 이용하고, 마찬가지로 ${\displaystyle A\mathbf {p} _{k-1}={\frac {1}{\alpha _{k-1}}}(\mathbf {r} _{k-1}-\mathbf {r} _{k})}$, ${\displaystyle \beta _{k}}$의 분자는 다음과 같이 다시 쓰면,

${\displaystyle \mathbf {r} _{k+1}^{\mathrm {T} }A\mathbf {p} _{k}={\frac {1}{\alpha _{k}}}\mathbf {r} _{k+1}^{\mathrm {T} }(\mathbf {r} _{k}-\mathbf {r} _{k+1})=-{\frac {1}{\alpha _{k}}}\mathbf {r} _{k+1}^{\mathrm {T} }\mathbf {r} _{k+1}}$

${\displaystyle \mathbf {r} _{k+1}}$ 과 ${\displaystyle \mathbf {r} _{k}}$ 처음 구상에 의해 직교이기 때문에 그분모는 다음과 같이 다시 쓰면,

${\displaystyle \mathbf {p} _{k}^{\mathrm {T} }A\mathbf {p} _{k}=(\mathbf {r} _{k}+\beta _{k-1}\mathbf {p} _{k-1})^{\mathrm {T} }A\mathbf {p} _{k}={\frac {1}{\alpha _{k}}}\mathbf {r} _{k}^{\mathrm {T} }(\mathbf {r} _{k}-\mathbf {r} _{k+1})={\frac {1}{\alpha _{k}}}\mathbf {r} _{k}^{\mathrm {T} }\mathbf {r} _{k}}$

검색 방향들 ${\displaystyle \mathbf {p} _{k}}$은 켤레이고 다시 그 유수들과 직교라는 사실을 이용하여 이 알고리즘에서 ${\displaystyle \alpha _{k}}$를 제거한 후${\displaystyle \beta }$를 얻는다.

function [x] = conjgrad(A,b,x)
    r=b-A*x;
    p=r;
    rsold=r'*r;

    for i=1:1e6
        Ap=A*p;
        alpha=rsold/(p'*Ap);
        x=x+alpha*p;
        r=r-alpha*Ap;
        rsnew=r'*r;
        if sqrt(rsnew)<1e-10
              break;
        end
        p=r+rsnew/rsold*p;
        rsold=rsnew;
    end
end

공역 구배법을 설명하기 위해서, 우리는 간단한 예를 들겠다. 다음과 같은 선형 계인 Ax = b을 고려하면

${\displaystyle \mathbf {A} \mathbf {x} ={\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}x_{1}\\x_{2}\end{bmatrix}}={\begin{bmatrix}1\\2\end{bmatrix}},}$

그 계의 거의 정확한 해를 찾기 위해서 우리는 초기 추측값을 설정하여 공역 구배법의 두 단계만에 수행할 것이다.

${\displaystyle \mathbf {x} _{0}={\begin{bmatrix}2\\1\end{bmatrix}}}$

참고에 따르면 정확한 해는 :: ${\displaystyle \mathbf {x} ={\begin{bmatrix}{\frac {1}{11}}\\\\{\frac {7}{11}}\end{bmatrix}}}$

우리의 첫번째 단계는 x₀과 연관된 유수 벡터 r₀를 계산 하는 것이다.

이 유수는 식 r₀ = b - Ax₀으로부터 계산되고, 그리고 우리에 경우에는 다음과 같다.

${\displaystyle \mathbf {r} _{0}={\begin{bmatrix}1\\2\end{bmatrix}}-{\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}2\\1\end{bmatrix}}={\begin{bmatrix}-8\\-3\end{bmatrix}}.}$

이것은 첫번째 반복 계산이기 때문에, 우리는 우리의 초기 검색방향p₀ 으로써 유수 벡터 r₀를 사용할 것이다.

앞으로의 반복 계산에서 p_k를 선택하는 방법은 변할 것이다.

우리는 이제 위의 관계를 이용하여 이제 스칼라 α₀를 계산한다.

${\displaystyle \alpha _{0}={\frac {\mathbf {r} _{0}^{\mathrm {T} }\mathbf {r} _{0}}{\mathbf {p} _{0}^{\mathrm {T} }\mathbf {Ap} _{0}}}={\frac {{\begin{bmatrix}-8&-3\end{bmatrix}}{\begin{bmatrix}-8\\-3\end{bmatrix}}}{{\begin{bmatrix}-8&-3\end{bmatrix}}{\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}-8\\-3\end{bmatrix}}}}={\frac {73}{331}}.}$

우리는 그 식을 이용하여 x₁을 계산 할 수 있다.

${\displaystyle \mathbf {x} _{1}=\mathbf {x} _{0}+\alpha _{0}\mathbf {p} _{0}={\begin{bmatrix}2\\1\end{bmatrix}}+{\frac {73}{331}}{\begin{bmatrix}-8\\-3\end{bmatrix}}={\begin{bmatrix}0.2356\\0.3384\end{bmatrix}}.}$

이 결과는 첫번째 반복 계산을 완료되고, 그 계에 "향상된" 거의 정확한 해x₁가 된다. 우리는 이제 이동하여 이 식을 이용하여 새로운 유수 벡터 r₁를 계산할 것이다.

${\displaystyle \mathbf {r} _{1}=\mathbf {r} _{0}-\alpha _{0}\mathbf {A} \mathbf {p} _{0}={\begin{bmatrix}-8\\-3\end{bmatrix}}-{\frac {73}{331}}{\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}-8\\-3\end{bmatrix}}={\begin{bmatrix}-0.2810\\0.7492\end{bmatrix}}.}$

그 과정속에서 우리의 다음 단계는 다음 검색 방향 p₁를 결정하는데 결과적으로 사용될 스칼라 β₀를 계산하는 것이다.

${\displaystyle \beta _{0}={\frac {\mathbf {r} _{1}^{\mathrm {T} }\mathbf {r} _{1}}{\mathbf {r} _{0}^{\mathrm {T} }\mathbf {r} _{0}}}={\frac {{\begin{bmatrix}-0.2810&0.7492\end{bmatrix}}{\begin{bmatrix}-0.2810\\0.7492\end{bmatrix}}}{{\begin{bmatrix}-8&-3\end{bmatrix}}{\begin{bmatrix}-8\\-3\end{bmatrix}}}}=0.0088.}$

이제, 이 스칼라 β₀를 이용하여, 우리는 그 관계를 이용하여 다음 검색 방향 p₁을 계산 할 수 있다.

${\displaystyle \mathbf {p} _{1}=\mathbf {r} _{1}+\beta _{0}\mathbf {p} _{0}={\begin{bmatrix}-0.2810\\0.7492\end{bmatrix}}+0.0088{\begin{bmatrix}-8\\-3\end{bmatrix}}={\begin{bmatrix}-0.3511\\0.7229\end{bmatrix}}.}$

우리는 α₀에서 사용했던 같은 방식으로 새로 얻은 p₁를 이용하여 이제 스칼라 α₁를 계산 할 수 있다.

${\displaystyle \alpha _{1}={\frac {\mathbf {r} _{1}^{\mathrm {T} }\mathbf {r} _{1}}{\mathbf {p} _{1}^{\mathrm {T} }\mathbf {Ap} _{1}}}={\frac {{\begin{bmatrix}-0.2810&0.7492\end{bmatrix}}{\begin{bmatrix}-0.2810\\0.7492\end{bmatrix}}}{{\begin{bmatrix}-0.3511&0.7229\end{bmatrix}}{\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}-0.3511\\0.7229\end{bmatrix}}}}=0.4122.}$

마지막으로, 우리는 x₁을 찾았던 같은 방법으로 x₂를 찾는다.

${\displaystyle \mathbf {x} _{2}=\mathbf {x} _{1}+\alpha _{1}\mathbf {p} _{1}={\begin{bmatrix}0.2356\\0.3384\end{bmatrix}}+0.4122{\begin{bmatrix}-0.3511\\0.7229\end{bmatrix}}={\begin{bmatrix}0.0909\\0.6364\end{bmatrix}}.}$

그 결과, x₂는 그 계의 해에 x₁ 과 x₀보다 "더나은" 근사치가 된다.

만일 제한된 정밀도 대신에 이 예에서 정확한 계산이 되었다면, 그 정확한 해는 이론적으로 n=2 반복수행 후 도달 할 것이다. (n은 그 계의 크기임).

• 신뢰전파
• 반복법