선형 시스템을 위한 반복법

수치 선형 대수의 가장 중요하고 일반적인 응용 사례 중 하나는 A*x = b 형식으로 표현할 수 있는 선형 시스템의 해입니다. A가 큰 희소 행렬인 경우 계산 실행 시간과 해의 정밀도 사이를 상호 절충하는 반복법을 사용하여 선형 시스템을 풀 수 있습니다. 여기에서는 방정식 A*x = b를 풀 때 MATLAB^®에서 사용할 수 있는 반복법에 대해 설명합니다.

직접법과 반복법 비교

선형 방정식 A*x = b를 푸는 데 사용할 수 있는 두 가지 유형의 방법이 있습니다.

직접법은 가우스 소거법(Gaussian Elimination)의 변형입니다. 이러한 방법은 LU 분해, QR 분해 또는 촐레스키 분해와 같은 행렬 연산을 통해 개별 행렬 요소를 직접 사용합니다. 직접법을 사용하면 높은 정밀도로 선형 방정식을 풀 수 있습니다. 하지만 이러한 방법은 큰 희소 행렬에 대한 연산을 수행하는 경우 속도가 느릴 수 있습니다. 직접법을 사용하여 선형 시스템을 푸는 속도는 계수 행렬의 밀도 및 채우기 패턴에 따라 크게 달라집니다.
예를 들어 다음 코드는 작은 선형 시스템을 풉니다.
```
A = magic(5);
b = sum(A,2);
x = A\b;
norm(A*x-b)
```
```
ans =

   1.4211e-14
```
MATLAB은 decomposition, lsqminnorm, linsolve 같은 함수뿐만 아니라 행렬 나눗셈 연산자 / 및 \를 통해서도 직접법을 구현합니다.
반복법은 유한한 횟수의 단계를 거친 후 선형 시스템에 대한 근사해(approximate solution)를 산출합니다. 이러한 방법은 실행 시간 단축을 위해 정밀도를 낮추는 것이 가능한 큰 연립방정식에 유용합니다. 반복법은 행렬-벡터 곱이나 추상 선형 연산자를 통해 간접적으로만 계수 행렬을 사용합니다. 반복법은 모든 행렬에 사용할 수 있지만 일반적으로 직접법 사용 시 속도가 느린 큰 희소 행렬에 적용됩니다. 간접법을 사용하여 선형 시스템을 푸는 속도는 직접법만큼 계수 행렬의 채우기 패턴에 크게 영향을 받지 않습니다. 그러나 반복법을 사용하려면 대개 각 특정 문제에 대해 파라미터를 조정해야 합니다.
예를 들어, 다음 코드는 양의 정부호 대칭 계수 행렬을 갖는 큰 희소 선형 시스템을 풉니다.
```
A = delsq(numgrid('L',400));
b = ones(size(A,1),1);
x = pcg(A,b,[],1000);
norm(b-A*x)
```
```
pcg converged at iteration 796 to a solution with relative residual 9.9e-07.

ans =

   3.4285e-04
```
MATLAB은 계수 행렬 A의 속성에 따라 다른 장단점을 갖는 다양한 반복법을 구현합니다.

직접법은 대개 이를 수행하는 데 사용할 수 있는 저장 공간이 충분한 경우 간접법보다 더 빠르고 더 일반적으로 적용 가능합니다. 일반적으로는 x = A\b를 먼저 사용해 보아야 합니다. 직접법이 너무 느리면 반복법을 사용해 볼 수 있습니다.

일반적인 반복 알고리즘

선형 방정식을 푸는 대부분의 반복 알고리즘은 대개 다음과 같은 절차를 따릅니다.

해 벡터 x0에 대한 초기 추측값으로 시작합니다. 사용자가 더 나은 추측값을 지정한 경우가 아니라면 초기 추측값은 대개 0으로 구성된 벡터입니다.
잔차 노름 res = norm(b-A*x0)을 계산합니다.
잔차와 지정된 허용오차를 비교합니다. res <= tol인 경우 계산을 끝내고 x0에 대해 계산된 답을 반환합니다.
잔차 값과 계산된 다른 수량에 따라 A*x0을 적용하고 벡터 x0의 크기와 방향을 업데이트합니다. 이 단계에서 대부분의 계산이 수행됩니다.
x0의 값이 허용오차를 충분히 충족할 때까지 2단계~4단계를 반복합니다.

반복법마다 4단계에서 x0의 크기와 방향을 업데이트하는 방식에 차이가 있으며 일부 반복법의 경우 2단계와 3단계의 수렴 조건이 약간 다르지만 이는 모든 반복 솔버가 따르는 기본 절차를 나타냅니다.

반복법 요약

MATLAB에는 선형 연립방정식을 풀기 위해 반복법을 구현하는 여러 개의 함수가 있습니다. 이러한 방법은 Ax = b의 해를 구하거나 노름(Norm) ||b – Ax||를 최소화하기 위해 설계되었습니다. 이러한 방법 중 몇몇은 유사성이 있고 동일한 기본 알고리즘에 기반하지만, 알고리즘마다 특정 상황에서 갖는 이점이 있습니다 [1], [2].

설명	참고
`pcg`(선조건 적용 켤레 기울기법)(Preconditioned Conjugate Gradient)	계수 행렬은 양의 정부호 대칭 행렬이어야 합니다. 제한된 개수의 벡터용 저장 공간만 필요하기 때문에, 양의 정부호 대칭 시스템에 가장 효과적인 솔버입니다.
`lsqr`(최소제곱법)(Least Squares)	사각 시스템에 사용 가능한 유일한 솔버입니다. 정규 방정식 `(A'A)x = A'*b`에 적용되는 켤레 기울기(PCG)법과 해석적으로 동일합니다.
`minres`(최소 잔차법)(Minimum Residual)	계수 행렬은 대칭이어야 하지만 양의 정부호 행렬일 필요는 없습니다. 반복할 때마다 2-노름에서 잔차 오차가 최소화되기 때문에, 알고리즘이 단계를 거치면서 더 나은 결과를 보입니다. 알고리즘이 해를 구하는 작업을 더 이상 진행할 수 없어 중단되는 경우에도 실패로 인한 영향을 받지 않습니다.
`symmlq`(대칭 LQ 계산법)(Symmetric LQ)	계수 행렬은 대칭이어야 하지만 양의 정부호 행렬일 필요는 없습니다. 투영 시스템에 대한 해를 구하고 잔차를 이전 모든 잔차에 직교하도록 유지합니다. 알고리즘이 해를 구하는 작업을 더 이상 진행할 수 없어 중단되는 경우에도 실패로 인한 영향을 받지 않습니다.
`bicg`(쌍켤레 기울기법)(Biconjugate Gradient)	계수 행렬은 정사각 행렬이어야 합니다. `bicg`는 리소스를 적게 사용하여 계산하지만, 수렴이 불규칙적이고 불안정합니다. `bicg`는 이 방법에 대한 개선 방법으로 다른 많은 반복 알고리즘이 개발되었기 때문에 역사적으로 중요합니다.
`bicgstab`(쌍켤레 기울기 안정법)(Biconjugate Gradients Stabilized)	계수 행렬은 정사각 행렬이어야 합니다. 안정성을 높이기 위해 BiCG 단계를 GMRES(1) 단계와 번갈아가며 사용합니다.
`bicgstabl`(쌍켤레 기울기 안정법(l))	계수 행렬은 정사각 행렬이어야 합니다. 안정성을 높이기 위해 BiCG 단계를 GMRES(2) 단계와 번갈아가며 사용합니다.
`cgs`(켤레 기울기 제곱법)(Conjugate Gradients Squared)	계수 행렬은 정사각 행렬이어야 합니다. 반복마다 `bicg`와 동일한 횟수의 연산을 필요로 하지만, 잔차 제곱을 사용한 전치 사용을 피합니다.
`gmres`(일반화된 최소 잔차법)(Generalized Minimum Residual)	계수 행렬은 정사각 행렬이어야 합니다. 각 반복에서 잔차 노름이 최소화되기 때문에, 가장 신뢰할 수 있는 알고리즘 중 하나입니다. 작업과 필요한 저장 공간이 반복 횟수에 따라 선형적으로 늘어납니다. 불필요한 작업과 저장 공간을 방지하기 위해 적절한 `restart` 값을 선택하는 것이 중요합니다.
`qmr`(준최소 잔차법)(Quasi-Minimal Residual)	계수 행렬은 정사각 행렬이어야 합니다. 반복당 오버헤드가 `bicg`보다 약간 더 높지만, 더 안정적입니다.
`tfqmr`(비전치 준최소 잔차법)(Transpose-free Quasi-Minimal Residual)	계수 행렬은 정사각 행렬이어야 합니다. 메모리가 제한적일 때 부정부호 대칭 시스템에 사용해 볼 수 있는 최적의 솔버입니다.

반복 솔버 선택하기

다음은 MATLAB의 반복 솔버를 표시한 순서도로, 각 솔버가 유용하게 사용되는 상황을 대략적으로 나타냅니다. 일반적으로 거의 모든 정사각 비대칭 문제에 gmres를 사용할 수 있습니다. gmres보다 쌍켤레 기울기법 알고리즘(bicg, bicgstab, cgs 등)이 더 효율적인 경우도 있지만, 예상할 수 없는 수렴 동작으로 인해 처음부터 gmres를 선택하는 것이 더 적절할 때가 많습니다.

Workflow to choose an appropriate iterative solver to use for a given problem.

선조건자

반복법의 수렴 속도는 계수 행렬의 스펙트럼(고유값)에 따라 결정됩니다. 따라서 더 유용한 스펙트럼(군집화된 고유값 또는 1에 가까운 조건수)을 갖도록 선형 시스템을 변환하면 반복법 대부분의 수렴과 안정성을 높일 수 있습니다. 이 변환은 선조건자라는 두 번째 행렬을 시스템에 적용하여 수행됩니다. 이 과정에서 다음 선형 시스템이

$A x = b$

다음과 같은 상응하는 시스템으로 변환됩니다.

$\tilde{A} \tilde{x} = \tilde{b} .$

이상적인 선조건자를 사용할 경우 모든 반복법이 1회 반복에서 수렴하기 때문에, 이러한 선조건자는 계수 행렬 A를 단위 행렬로 변환합니다. 실제로 좋은 선조건자를 찾으려면 상호 절충이 필요합니다. 변환은 세 가지 방법 즉, 왼쪽 선조건 지정, 오른쪽 선조건 지정 또는 분할 선조건 지정 중 하나로 수행됩니다.

첫 번째 경우는 선조건자 행렬 M이 A의 왼쪽에 위치하는 왼쪽 선조건 지정입니다.

$(M^{- 1} A) x = (M^{- 1} b) .$

왼쪽 선조건 지정을 사용하는 반복 솔버는 다음과 같습니다.

오른쪽 선조건 지정에서는 M이 A의 오른쪽에 위치합니다.

$(A M^{- 1}) (M x) = b .$

오른쪽 선조건 지정을 사용하는 반복 솔버는 다음과 같습니다.

마지막으로, 대칭 계수 행렬 A에 대해 분할 선조건 지정을 사용하면 변환된 시스템이 계속 대칭을 유지할 수 있습니다. 선조건자 $M = H H^{T}$ 가 분할되고 인수들이 A의 서로 다른 쪽에 나타납니다.

$(H^{- 1} A H^{- T}) H^{T} x = (H^{- 1} b)$

분할 선조건 적용 시스템의 솔버 알고리즘은 위 방정식을 바탕으로 하지만, 실제로는 H를 계산할 필요가 없습니다. 솔버 알고리즘은 M으로 직접 곱하여 풉니다.

분할 선조건 지정을 사용하는 반복 솔버는 다음과 같습니다.

모든 경우에서 선조건자 M은 반복법에 대한 수렴을 가속화하기 위해 선택되었습니다. 반복 해의 잔차 오차가 반복 간에 정체되거나 별다른 진척이 없는 경우 대개 선조건자 행렬을 생성해서 문제에 포함시켜야 합니다.

MATLAB에서 반복 솔버를 사용하면 단일 선조건자 행렬 M 또는 두 개의 선조건자 행렬 인수(예: M = M₁M₂)를 지정할 수 있습니다. 이렇게 하면 선조건자를 분해된 형태로 쉽게 지정할 수 있습니다(예: M = LU). 분할 선조건이 적용된 경우, M = HH^T도 성립하는데 M1, M2 입력값과 H 인수 사이에 관계가 없다는 사실에 유의하십시오.

경우에 따라 주어진 문제의 수학 모델에는 기본적으로 선조건자가 존재합니다. 기본적인 선조건자가 없는 경우 다음 표에 나와 있는 불완전 행렬 분해 중 하나를 사용하여 선조건자 행렬을 생성할 수 있습니다. 본질적으로 불완전 행렬 분해는 계산이 빠른 불완전한 직접법입니다.

함수 행렬 분해 설명

함수	행렬 분해	설명
`ilu`	$A \approx LU$	정사각 행렬 또는 직사각 행렬에 대한 불완전 LU 분해.
`ichol`	$A \approx L L^{*}$	양의 정부호 대칭 행렬에 대한 불완전 촐레스키 분해.

ilu

$A \approx LU$

정사각 행렬 또는 직사각 행렬에 대한 불완전 LU 분해.

ichol

$A \approx L L^{*}$

양의 정부호 대칭 행렬에 대한 불완전 촐레스키 분해.

ilu 및 ichol에 대한 자세한 내용은 불완전 행렬 분해 항목을 참조하십시오.

선조건자 예제

정사각 2차원 영역에서 라플라스 방정식(Laplace's equation)에 대한 5개 점을 가지는 유한 차분 근삿값이 있다고 가정해 보겠습니다. 다음 명령은 선조건자를 M = L*L'으로 하는 선조건 적용 켤레 기울기(PCG)법을 사용합니다. 여기서 L은 A에 대한 0 채우기(Zero-fill) 불완전 촐레스키 인수입니다. 이 시스템의 경우 선조건자 행렬을 지정하지 않으면 pcg는 해를 구할 수 없습니다.

A = delsq(numgrid('S',250));
b = ones(size(A,1),1);
tol = 1e-3;
maxit = 100;
L = ichol(A);
x = pcg(A,b,tol,maxit,L,L');

pcg converged at iteration 92 to a solution with relative residual 0.00076.

pcg는 지정된 허용오차 달성을 위해 92회의 반복을 필요로 합니다. 그러나 다른 선조건자를 사용하여 더 나은 결과를 얻을 수 있습니다. 예를 들어, ichol을 사용하여 수정된 불완전 촐레스키를 생성하는 경우 pcg는 단지 39회의 반복 후에 지정된 허용오차를 충족할 수 있습니다.

L = ichol(A,struct('type','nofill','michol','on'));
x = pcg(A,b,tol,maxit,L,L');

pcg converged at iteration 39 to a solution with relative residual 0.00098.

평형화 기법과 재정렬

힘든 계산을 요하는 문제의 경우 ilu 또는 ichol로 직접 생성된 선조건자보다 더 나은 선조건자가 필요할 수 있습니다. 예를 들어, 더 나은 품질의 선조건자를 생성하거나 수행되는 계산량을 최소화하고자 할 수 있습니다. 이러한 경우 평형화 기법을 사용하여 계수 행렬을 좀 더 대각선 우위 행렬로 만들고(이로 인해 선조건자 품질이 향상될 수 있음), 재정렬을 사용하여 행렬 인수에서 0이 아닌 요소의 개수를 최소화할 수 있습니다(이로 인해 메모리 요구 사항이 줄어들고 이후 계산의 효율성이 향상될 수 있음).

평형화 기법과 재정렬을 모두 사용하여 선조건자를 생성하는 경우 절차는 다음과 같습니다.

계수 행렬에 equilibrate를 사용합니다.
평형화 기법이 적용된 행렬을 dissect 또는 symrcm과 같은 희소 행렬 재정렬 함수를 사용하여 재정렬합니다.
ilu 또는 ichol을 사용하여 최종 선조건자를 생성합니다.

다음은 희소 계수 행렬에 대한 선조건자를 생성하기 위해 평형화 기법과 재정렬을 사용하는 예제입니다.

계수 행렬 A와 선형 방정식의 우변에 해당하는 1로 구성된 벡터 b를 만듭니다. A의 조건수 추정값을 계산합니다.
```
load west0479;
A = west0479;
b = ones(size(A,1),1);
condest(A)
```
```
ans =

   1.4244e+12
```
equilibrate를 사용하여 계수 행렬의 조건수를 개선합니다.
```
[P,R,C] = equilibrate(A);
Anew = R*P*A*C;
bnew = R*P*b;
condest(Anew)
```
```
ans =

   5.1042e+04
```
평형화 기법이 적용된 행렬을 dissect를 사용하여 재정렬합니다.
```
q = dissect(Anew);
Anew = Anew(q,q);
bnew = bnew(q);
```
불완전 LU 분해를 사용하여 선조건자를 생성합니다.
```
[L,U] = ilu(Anew);
```
선조건자 행렬, 허용오차 1e-10, 최대 외부 반복 50회, 내부 반복 30회를 사용하여 gmres로 선형 시스템을 풉니다.
```
tol = 1e-10;
maxit = 50;
restart = 30;
[xnew, flag, relres] = gmres(Anew,bnew,restart,tol,maxit,L,U);
x(q) = xnew;
x = C*x(:);
```
이제 gmres로 반환된 상대 잔차 relres(선조건자가 포함됨)를 선조건자가 없는 상대 잔차 resnew와 평형화 기법이 없는 상대 잔차 res와 비교합니다. 결과를 보면, 선형 시스템이 모두 동일하더라도 방법마다 각 요소에 서로 다른 가중치를 적용하기 때문에 잔차 값에 큰 영향을 미칠 수 있음을 알 수 있습니다.
```
relres
resnew = norm(Anew*xnew - bnew) / norm(bnew)
res = norm(A*x - b) / norm(b)
```
```
relres =
   8.7537e-11
resnew =
   3.6805e-08
res =
   5.1415e-04
```

행렬 대신 선형 연산자 사용하기

MATLAB에서 반복 솔버를 사용할 경우 A에 대한 숫자형 행렬을 제공할 필요가 없습니다. 솔버에서 수행하는 계산에는 행렬-벡터 곱셈 A*x 또는 A'*x의 결과가 사용되기 때문에, 그러한 선형 연산의 결과를 계산하는 함수를 대신 제공하면 됩니다. 이러한 수량을 계산하는 함수를 흔히 선형 연산자라고 합니다.

계수 행렬 A 대신에 선형 연산자를 사용할 뿐 아니라, 선조건자 M에 대한 행렬 대신에도 선형 연산자를 사용할 수 있습니다. 이 경우 함수는 솔버에 대한 함수 도움말 페이지에 표시된 대로 M\x 또는 M'\x를 계산해야 합니다.

선형 연산자를 사용하면 A 또는 M의 패턴을 활용하기 때문에 솔버가 명시적으로 행렬을 사용하여 비희소 행렬-벡터 곱셈을 수행할 때보다 선형 연산의 값을 보다 효율적으로 계산할 수 있습니다. 또한 대개 선형 연산자가 행렬을 전혀 구성하지 않은 상태로 행렬-벡터 곱셈의 결과를 계산하기 때문에, 계수 또는 선조건자 행렬을 저장할 메모리가 필요 없습니다.

예를 들어, 다음과 같은 계수 행렬이 있다고 가정해 보겠습니다.

A = [2 -1  0  0  0  0;
    -1  2 -1  0  0  0;
     0 -1  2 -1  0  0;
     0  0 -1  2 -1  0;
     0  0  0 -1  2 -1;
     0  0  0  0 -1  2];

A에 벡터를 곱할 때 결과 벡터의 대부분의 요소는 0입니다. 결과에 있는 0이 아닌 요소는 A의 0이 아닌 삼중대각 요소에 해당합니다. 따라서 지정된 벡터 x에 대해 A*x의 값을 계산하려면 선형 연산자 함수는 3개의 벡터를 함께 더하기만 하면 됩니다.

function y = linearOperatorA(x)
y = -1*[0; x(1:end-1)] ...
  + 2*x ...
  + -1*[x(2:end); 0];
end

대부분의 반복 솔버에서는 A에 대한 선형 연산자 함수가 A*x 값을 반환해야 합니다. 마찬가지로, 선조건자 행렬 M에 대해 일반적으로 M\x를 계산해야 합니다. 솔버 lsqr , qmr 및 bicg의 경우 선형 연산자 함수는 A'*x 또는 M'\x의 값도 반환해야 합니다(요청될 경우). 선형 연산자 함수의 예제와 설명은 반복 솔버 함수 도움말 페이지를 참조하십시오.

참고 문헌

[1] Barrett, R., M. Berry, T. F. Chan, et al., Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, SIAM, Philadelphia, 1994.

[2] Saad, Yousef, Iterative Methods for Sparse Linear Equations. PWS Publishing Company, 1996.