선형 시스템을 위한 반복법
수치 선형 대수의 가장 중요하고 일반적인 응용 사례 중 하나는 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].
설명 | 참고 |
---|---|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
반복 솔버 선택하기
다음은 MATLAB의 반복 솔버를 표시한 순서도로, 각 솔버가 유용하게 사용되는 상황을 대략적으로 나타냅니다. 일반적으로 거의 모든 정사각 비대칭 문제에 gmres
를 사용할 수 있습니다. gmres
보다 쌍켤레 기울기법 알고리즘(bicg
, bicgstab
, cgs
등)이 더 효율적인 경우도 있지만, 예상할 수 없는 수렴 동작으로 인해 처음부터 gmres
를 선택하는 것이 더 적절할 때가 많습니다.
선조건자
반복법의 수렴 속도는 계수 행렬의 스펙트럼(고유값)에 따라 결정됩니다. 따라서 더 유용한 스펙트럼(군집화된 고유값 또는 1에 가까운 조건수)을 갖도록 선형 시스템을 변환하면 반복법 대부분의 수렴과 안정성을 높일 수 있습니다. 이 변환은 선조건자라는 두 번째 행렬을 시스템에 적용하여 수행됩니다. 이 과정에서 다음 선형 시스템이
다음과 같은 상응하는 시스템으로 변환됩니다.
이상적인 선조건자를 사용할 경우 모든 반복법이 1회 반복에서 수렴하기 때문에, 이러한 선조건자는 계수 행렬 A를 단위 행렬로 변환합니다. 실제로 좋은 선조건자를 찾으려면 상호 절충이 필요합니다. 변환은 세 가지 방법 즉, 왼쪽 선조건 지정, 오른쪽 선조건 지정 또는 분할 선조건 지정 중 하나로 수행됩니다.
첫 번째 경우는 선조건자 행렬 M이 A의 왼쪽에 위치하는 왼쪽 선조건 지정입니다.
왼쪽 선조건 지정을 사용하는 반복 솔버는 다음과 같습니다.
오른쪽 선조건 지정에서는 M이 A의 오른쪽에 위치합니다.
오른쪽 선조건 지정을 사용하는 반복 솔버는 다음과 같습니다.
마지막으로, 대칭 계수 행렬 A
에 대해 분할 선조건 지정을 사용하면 변환된 시스템이 계속 대칭을 유지할 수 있습니다. 선조건자 가 분할되고 인수들이 A의 서로 다른 쪽에 나타납니다.
분할 선조건 적용 시스템의 솔버 알고리즘은 위 방정식을 바탕으로 하지만, 실제로는 H를 계산할 필요가 없습니다. 솔버 알고리즘은 M
으로 직접 곱하여 풉니다.
분할 선조건 지정을 사용하는 반복 솔버는 다음과 같습니다.
모든 경우에서 선조건자 M은 반복법에 대한 수렴을 가속화하기 위해 선택되었습니다. 반복 해의 잔차 오차가 반복 간에 정체되거나 별다른 진척이 없는 경우 대개 선조건자 행렬을 생성해서 문제에 포함시켜야 합니다.
MATLAB에서 반복 솔버를 사용하면 단일 선조건자 행렬 M 또는 두 개의 선조건자 행렬 인수(예: M = M1M2)를 지정할 수 있습니다. 이렇게 하면 선조건자를 분해된 형태로 쉽게 지정할 수 있습니다(예: M = LU). 분할 선조건이 적용된 경우, M = HHT도 성립하는데 M1
, M2
입력값과 H 인수 사이에 관계가 없다는 사실에 유의하십시오.
경우에 따라 주어진 문제의 수학 모델에는 기본적으로 선조건자가 존재합니다. 기본적인 선조건자가 없는 경우 다음 표에 나와 있는 불완전 행렬 분해 중 하나를 사용하여 선조건자 행렬을 생성할 수 있습니다. 본질적으로 불완전 행렬 분해는 계산이 빠른 불완전한 직접법입니다.
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
과 같은 희소 행렬 재정렬 함수를 사용하여 재정렬합니다.
다음은 희소 계수 행렬에 대한 선조건자를 생성하기 위해 평형화 기법과 재정렬을 사용하는 예제입니다.
계수 행렬
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.