Main Content

bicgstabl

선형 연립방정식 풀기 — 쌍켤레 기울기 안정 (l)법(Stabilized BiConjugate Gradients (l) Method)

설명

예제

x = bicgstabl(A,b)쌍켤레 기울기 (l) 안정법(Biconjugate Gradients Stabilized (l) Method) 방법을 사용하여 x에 대한 선형 연립방정식 A*x = b를 풉니다. 시도가 성공한 경우 bicgstabl은 수렴을 확인하는 메시지를 표시합니다. 최대 반복 횟수 이후에도 bicgstabl이 수렴하지 않거나 어떠한 이유로든 중단될 경우 상대 잔차 norm(b-A*x)/norm(b)와, 이 계산이 중지된 반복 횟수가 포함된 진단 메시지가 표시됩니다.

예제

x = bicgstabl(A,b,tol)은 이 방법의 허용오차를 지정합니다. 디폴트 허용오차는 1e-6입니다.

예제

x = bicgstabl(A,b,tol,maxit)는 적용할 최대 반복 횟수를 지정합니다. maxit회의 반복 내에서 수렴하지 않을 경우 bicgstabl은 진단 메시지를 표시합니다.

예제

x = bicgstabl(A,b,tol,maxit,M)은 선조건자 행렬 M을 지정하고 y에 대해 시스템 AM1y=b를 풀어 x를 계산합니다. 여기서 y=Mx입니다. 선조건자 행렬을 사용하면 문제의 수치적 속성과 계산 효율성을 개선할 수 있습니다.

예제

x = bicgstabl(A,b,tol,maxit,M1,M2)M = M1*M2가 되도록 선조건자 행렬 M의 인수를 지정합니다.

예제

x = bicgstabl(A,b,tol,maxit,M1,M2,x0)은 해 벡터 x에 대한 초기 추측값을 지정합니다. 디폴트 값은 0으로 구성된 벡터입니다.

예제

[x,flag] = bicgstabl(___)은 알고리즘이 성공적으로 수렴하는지 여부를 지정하는 플래그를 반환합니다. flag = 0이면 수렴이 성공한 것입니다. 이 출력 구문은 위에 열거된 구문에 나와 있는 입력 인수를 원하는 대로 조합하여 사용할 수 있습니다. flag 출력값을 지정할 경우 bicgstabl은 진단 메시지를 표시하지 않습니다.

예제

[x,flag,relres] = bicgstabl(___)은 상대 잔차 norm(b-A*x)/norm(b)도 반환합니다. flag0이면 relres <= tol입니다.

예제

[x,flag,relres,iter] = bicgstabl(___)x가 계산된 반복 횟수 iter도 반환합니다.

예제

[x,flag,relres,iter,resvec] = bicgstabl(___)은 첫 번째 잔차 norm(b-A*x0)을 포함하여 각 4분할 반복에서 잔차 노름으로 구성된 벡터도 반환합니다.

예제

모두 축소

디폴트 설정으로 bicgstabl을 사용하여 정사각 선형 시스템을 풀고 풀이 과정에 사용된 허용오차 및 반복 횟수를 조정합니다.

50%의 밀도를 갖는 희소 형식의 확률 행렬 A를 만듭니다. Ax=b의 우변에 해당하는 확률 벡터 b도 만듭니다.

rng default
A = sprand(400,400,.5);
A = A'*A;
b = rand(400,1);

bicgstabl을 사용하여 Ax=b를 풉니다. 출력 표시에는 상대 잔차 오차 b-Axb의 값이 포함됩니다.

x = bicgstabl(A,b);
bicgstabl stopped at iteration 20 without converging to the desired tolerance 1e-06
because the maximum number of iterations was reached.
The iterate returned (number 20) has relative residual 0.09.

기본적으로 bicgstabl은 20회 반복과 1e-6의 허용오차를 사용하는데 이 알고리즘으로는 이 행렬에 대해 20회의 반복 내에 수렴할 수 없습니다. 잔차는 아직 큽니다. 이는 더 많은 반복(또는 선조건자 행렬)이 필요함을 알려줍니다. 알고리즘이 쉽게 수렴할 수 있도록 큰 허용오차를 사용할 수도 있습니다.

1e-4의 허용오차와 100회의 반복으로 시스템을 다시 풉니다.

x = bicgstabl(A,b,1e-4,100);
bicgstabl stopped at iteration 100 without converging to the desired tolerance 0.0001
because the maximum number of iterations was reached.
The iterate returned (number 100) has relative residual 0.034.

허용오차가 더 느슨해지고 반복 횟수가 늘어나도 잔차 오차가 크게 개선되지 않았습니다. 반복 알고리즘이 이와 같이 교착 상태에 빠진다는 것은 선조건자 행렬이 필요함을 알려줍니다.

A의 불완전 촐레스키 분해를 계산하고, bicgstabl에 대한 선조건자 입력값으로 L' 인수를 사용합니다.

L = ichol(A);
x = bicgstabl(A,b,1e-4,100,L');
bicgstabl converged at iteration 15.2 to a solution with relative residual 7.7e-05.

선조건자를 사용한 결과 bicgstabl이 수렴할 수 있을 만큼 문제의 수치적 속성이 개선됩니다.

bicgstabl에 선조건자 행렬을 사용하여 선형 시스템을 푸는 효과를 검토합니다.

479×479 실수 비대칭 희소 행렬인 west0479를 불러옵니다.

load west0479
A = west0479;

Ax=b의 참(True) 해가 1로만 구성된 벡터가 되도록 b를 정의합니다.

b = sum(A,2);

허용오차와 최대 반복 횟수를 설정합니다.

tol = 1e-12;
maxit = 20;

bicgstabl을 사용하여 요청된 허용오차와 반복 횟수로 해를 구합니다. 풀이 과정에 대한 정보를 반환하는 5개의 출력값을 지정합니다.

  • xA*x = b의 계산된 해입니다.

  • fl0은 알고리즘의 수렴 여부를 나타내는 플래그입니다.

  • rr0은 계산된 답 x의 상대 잔차입니다.

  • it0x가 계산된 반복 횟수입니다.

  • rv0b-Ax의 잔차 내역으로 구성된 벡터입니다.

[x,fl0,rr0,it0,rv0] = bicgstabl(A,b,tol,maxit);
fl0
fl0 = 1
rr0
rr0 = 1
it0
it0 = 0

bicgstabl이 요청된 반복 횟수인 20회 이내에 요청된 허용오차 1e-12으로 수렴하지 않으므로 fl01이 됩니다. 사실상 bicgstabl의 동작이 좋지 않기 때문에 초기 추측값 x0 = zeros(size(A,2),1)이 최적의 해가 되고 it0 = 0에서 볼 수 있듯이 이 해가 반환됩니다.

이러한 느린 수렴을 개선할 수 있도록 선조건자 행렬을 지정할 수 있습니다. A는 비대칭 행렬이므로 ilu를 사용하여 선조건자 M=L U를 생성합니다. 대각선상에 있지 않으면서 값이 1e-6보다 작은 요소는 무시하도록 기각 허용오차를 지정합니다. LUbicgstabl에 대한 입력값으로 지정하여 선조건이 적용된 시스템 A M-1M x=b를 풉니다.

setup = struct('type','ilutp','droptol',1e-6);
[L,U] = ilu(A,setup);
[x1,fl1,rr1,it1,rv1] = bicgstabl(A,b,tol,maxit,L,U);
fl1
fl1 = 0
rr1
rr1 = 1.0249e-15
it1
it1 = 2

ilu 선조건자를 사용한 결과 2번째 반복에서 미리 정해진 허용오차 1e-12보다 작은 상대 잔차가 생성됩니다. 출력값 rv1(1)norm(b)가 되고, 출력값 rv1(end)norm(b-A*x1)이 됩니다.

각 반복마다 상대 잔차를 플로팅하면 bicgstabl의 진행률을 추적할 수 있습니다. 각 해의 잔차 내역을 플로팅하고, 지정된 허용오차에 대한 선을 추가합니다.

semilogy(0:length(rv0)-1,rv0/norm(b),'-o')
hold on
semilogy(0:length(rv1)-1,rv1/norm(b),'-o')
yline(tol,'r--');
legend('No preconditioner','ILU preconditioner','Tolerance','Location','East')
xlabel('Iteration number')
ylabel('Relative residual')

Figure contains an axes object. The axes object with xlabel Iteration number, ylabel Relative residual contains 3 objects of type line, constantline. These objects represent No preconditioner, ILU preconditioner, Tolerance.

bicgstabl에 해의 초기 추측값을 제공하는 효과를 검토합니다.

삼중대각 희소 행렬을 만듭니다. x의 예상 해가 1로 구성된 벡터가 되도록 각 행의 합을 Ax=b의 우변에 대한 벡터로 사용합니다.

n = 900;
e = ones(n,1);
A = spdiags([e 2*e e],-1:1,n,n);
b = sum(A,2);

bicgstabl을 사용하여 Ax=b를 두 번 풉니다. 한 번은 디폴트 초기 추측값을 사용하여 풀고, 다른 한 번은 해에 대한 양호한 초기 추측값을 사용하여 풉니다. 두 해 모두에 대해 20회의 반복과 디폴트 허용오차를 사용합니다. 모든 요소가 0.99로 이루어진 벡터로 두 번째 해의 초기 추측값을 지정합니다.

maxit = 20;
x1 = bicgstabl(A,b,[],maxit);
bicgstabl converged at iteration 9.2 to a solution with relative residual 9.5e-07.
x0 = 0.99*e;
x2 = bicgstabl(A,b,[],maxit,[],[],x0);
bicgstabl converged at iteration 2 to a solution with relative residual 5.4e-07.

이 예제에서는 초기 추측값을 제공한 결과 bicgstabl이 더 빨리 수렴되었습니다.

중간 결과 반환하기

for 루프에서 bicgstabl을 호출하여 중간 결과를 얻는 데도 초기 추측값을 사용할 수 있습니다. 솔버에 대한 각 호출은 몇 회의 반복을 수행한 후 계산된 해를 저장합니다. 그런 다음 저장된 해를 다음 반복 배치의 초기 벡터로 사용합니다.

예를 들어, 다음 코드는 100회의 반복을 4번 수행하며, for 루프를 통과한 후 매번 해 벡터를 저장합니다.

x0 = zeros(size(A,2),1);
tol = 1e-8;
maxit = 100;
for k = 1:4
    [x,flag,relres] = bicgstabl(A,b,tol,maxit,[],[],x0);
    X(:,k) = x;
    R(k) = relres;
    x0 = x;
end

X(:,k)는 for 루프의 반복 k회에서 계산된 해 벡터이고, R(k)는 이 해의 상대 잔차입니다.

bicgstabl에 계수 행렬 A 대신 A*x를 계산하는 함수 핸들을 제공하여 선형 시스템을 풉니다.

gallery에 의해 생성되는 윌킨슨 테스트 행렬 중 하나는 21×21 삼중대각 행렬입니다. 행렬을 미리 봅니다.

A = gallery('wilk',21)
A = 21×21

    10     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     1     9     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     1     8     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     1     7     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     1     6     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     1     5     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     1     4     1     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     1     3     1     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     1     2     1     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     1     1     1     0     0     0     0     0     0     0     0     0     0
      ⋮

윌킨슨 행렬은 특수한 구조를 갖고 있으므로 연산 A*x를 함수 핸들로 나타낼 수 있습니다. A에 벡터를 곱할 때 결과 벡터의 대부분의 요소는 0입니다. 결과에 있는 0이 아닌 요소는 A의 0이 아닌 삼중대각선 요소에 해당합니다. 또한 주대각선에만 0과 1이 아닌 요소가 있습니다.

표현식 Ax는 다음과 같이 표현됩니다.

Ax=[1010001910001810017100161001510014100130001000110][x1x2x3x4x5x21]=[10x1+x2x1+9x2+x3x2+8x3+x4x19+9x20+x21x20+10x21].

결과 벡터는 다음과 같이 세 벡터의 합으로 작성할 수 있습니다.

Ax=[0+10x1+x2x1+9x2+x3x2+8x3+x4x19+9x20+x21x20+10x21+0]=[0x1x20]+[10x19x210x21]+[x2x210].

MATLAB®에서 이러한 벡터를 만들어서 더하여 A*x의 값을 제공하는 함수를 작성합니다.

function y = afun(x)
y = [0; x(1:20)] + ...
    [(10:-1:0)'; (1:10)'].*x + ...
    [x(2:21); 0];
end

(이 함수는 예제 끝에 로컬 함수로 저장되어 있습니다.)

이제 A*x를 계산하는 함수 핸들을 bicgstabl에 제공하여 선형 시스템 Ax=b를 풉니다. 허용오차 1e-12 및 50회 반복을 사용합니다.

b = ones(21,1);
tol = 1e-12;  
maxit = 50;
x1 = bicgstabl(@afun,b,tol,maxit)
bicgstabl converged at iteration 5.2 to a solution with relative residual 2e-15.
x1 = 21×1

    0.0910
    0.0899
    0.0999
    0.1109
    0.1241
    0.1443
    0.1544
    0.2383
    0.1309
    0.5000
      ⋮

afun(x1)이 1로 구성된 벡터를 생성하는지 확인합니다.

afun(x1)
ans = 21×1

    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
      ⋮

로컬 함수(Local Function)

function y = afun(x)
y = [0; x(1:20)] + ...
    [(10:-1:0)'; (1:10)'].*x + ...
    [x(2:21); 0];
end

입력 인수

모두 축소

계수 행렬로, 정사각 행렬 또는 함수 핸들로 지정됩니다. 이 행렬은 선형 시스템 A*x = b의 계수 행렬입니다. 일반적으로 A는 큰 희소 행렬이거나 큰 희소 행렬과 열 벡터의 곱을 반환하는 함수 핸들입니다.

A를 함수 핸들로 지정

선택적으로 계수 행렬을 행렬이 아닌 함수 핸들로 지정할 수도 있습니다. 함수 핸들은 전체 계수 행렬을 구성하는 대신 행렬-벡터 곱을 반환하여 계산을 보다 효율적으로 만듭니다.

함수 핸들을 사용하려면 함수 시그니처 function y = afun(x)를 사용하십시오. 함수를 파라미터화하기에서는 필요한 경우 함수 afun에 추가 파라미터를 제공하는 방법을 설명합니다. 함수 호출 afun(x)A*x의 값을 반환해야 합니다.

데이터형: double | function_handle
복소수 지원 여부:

선형 방정식의 우변으로, 열 벡터로 지정됩니다. b의 길이는 size(A,1)이어야 합니다.

데이터형: double
복소수 지원 여부:

방법에 대한 허용오차로, 양의 스칼라로 지정됩니다. 이 입력값을 사용하여 계산에서 정확도와 런타임 간의 상호 절충 관계를 조정할 수 있습니다. bicgstabl 함수가 성공하려면 허용된 반복 횟수 내에서 허용오차를 충족해야 합니다. tol 값이 작을수록 계산을 완료하려면 답이 더 정밀해야 합니다.

데이터형: double

최대 반복 횟수로, 양의 정수 스칼라로 지정됩니다. bicgstabl 함수가 더 많은 반복을 거쳐 허용오차 tol을 충족할 수 있도록 하려면 maxit의 값을 늘리십시오. 일반적으로 tol의 값이 작을수록 계산이 성공하려면 더 많은 반복이 필요합니다.

선조건자 행렬로, 행렬 또는 함수 핸들의 개별 인수로 지정됩니다. 선조건자 행렬 M 또는 그 행렬 계수 M = M1*M2를 지정하여 선형 시스템의 수치적 특성을 개선하고 bicgstabl 함수가 빠르게 수렴하도록 할 수 있습니다. 불완전 행렬 분해 함수 iluichol을 사용하여 선조건자 행렬을 생성할 수 있습니다. 분해 전에 equilibrate를 사용하여 계수 행렬의 조건수를 개선할 수도 있습니다. 선조건자에 대한 자세한 내용은 선형 시스템을 위한 반복법 항목을 참조하십시오.

bicgstabl 함수는 지정되지 않은 선조건자를 단위 행렬로 처리합니다.

M을 함수 핸들로 지정

선택적으로 M, M1 또는 M2를 행렬이 아닌 함수 핸들로 지정할 수도 있습니다. 함수 핸들은 전체 선조건자 행렬을 구성하는 대신 함수-벡터 연산을 수행하여 계산을 보다 효율적으로 만듭니다.

함수 핸들을 사용하려면 함수 시그니처 function y = mfun(x)를 사용하십시오. 함수를 파라미터화하기에서는 필요한 경우 함수 mfun에 추가 파라미터를 제공하는 방법을 설명합니다. 함수 호출 mfun(x)M\x 또는 M2\(M1\x)의 값을 반환해야 합니다.

데이터형: double | function_handle
복소수 지원 여부:

초기 추측값으로, 길이가 size(A,2)인 열 벡터로 지정됩니다. bicgstabl에 디폴트 값인 0으로 구성된 벡터보다 더 합리적인 초기 추측값 x0을 제공할 수 있는 경우 이 초기 추측값은 계산 시간을 절약하고 알고리즘이 더 빨리 수렴하도록 할 수 있습니다.

데이터형: double
복소수 지원 여부:

출력 인수

모두 축소

선형 시스템 해로, 열 벡터로 반환됩니다. 이 출력값은 선형 시스템 A*x = b에 대한 근사해를 제공합니다. 계산이 성공한 경우(flag = 0), relrestol보다 작거나 같습니다.

계산이 성공하지 않으면(flag ~= 0), bicgstabl 함수가 반환하는 해 x는 모든 반복에 대해 최소 잔차 노름이 계산된 해입니다.

수렴 플래그로, 다음 표에 있는 스칼라 값 중 하나로 반환됩니다. 수렴 플래그는 계산의 성공 여부를 나타내며 여러 가지 형태의 실패를 구분합니다.

플래그 값

수렴

0

성공 — bicgstabl 함수가 maxit 반복 횟수 내에, 기대한 허용오차 tol로 수렴되었습니다.

1

실패 — bicgstabl 함수가 maxit회만큼 반복되었지만 수렴되지 않았습니다.

2

실패 — 선조건자 행렬 M 또는 M = M1*M2의 조건이 나쁩니다.

3

실패 — 연속된 2회의 반복이 동일한 결과를 반환한 후 bicgstabl 함수가 정체되었습니다.

4

실패 — bicgstabl 알고리즘으로 계산된 스칼라 수량 중 하나가 너무 작아지거나 커져 계산을 계속할 수 없습니다.

상대 잔차 오차로, 스칼라로 반환됩니다. 상대 잔차 오차 relres = norm(b-A*x)/norm(b)는 답의 정확도를 나타냅니다. 계산이 maxit회의 반복 내에서 허용오차 tol로 수렴하면 relres <= tol입니다.

데이터형: double

반복 횟수로, 스칼라로 반환됩니다. 이 출력값은 x의 답이 계산된 반복 횟수를 나타냅니다. bicgstabl의 각 외부 반복은 4개의 내부 반복을 포함하므로 iter은 소수 형식의 반복 번호로 반환될 수 있습니다.

데이터형: double

잔차 오차로, 벡터로 반환됩니다. 잔차 오차 norm(b-A*x)는 알고리즘이 주어진 x 값에 얼마나 가깝게 수렴하는지를 나타냅니다. resvec 내 요소 개수는 반복을 4분할한 개수와 동일합니다. resvec의 내용을 검사하여 tol 또는 maxit의 값을 변경할지 여부를 결정할 수 있습니다.

데이터형: double

세부 정보

모두 축소

쌍켤레 기울기 (l) 안정법(Biconjugate Gradients Stabilized (l) Method)

쌍켤레 기울기 안정 l(BiCGSTABL) 알고리즘은 BiCG 방법을 개선하는 과정에서 개발된 BiCGSTAB 방법을 개선하기 위해 개발되었습니다.

BiCGSTAB와 마찬가지로 BiCGSTABL 알고리즘도 GMRES 단계를 사용하여, BiCG에 나타나는 불규칙 수렴 동작을 최소화합니다. 그런데 BiCGSTABL은 BiCGSTAB의 GMRES(1) 단계가 아닌 GMRES(2) 단계를 사용하므로 보다 덜 빈번하게 정체되는 더 나은 개선을 제공할 수 있습니다[1].

  • 대부분의 반복법의 수렴 여부는 계수 행렬의 조건수 cond(A)에 따라 결정됩니다. A가 정사각 행렬인 경우 equilibrate를 사용하여 이 행렬의 조건수를 개선할 수 있으며, 이로 인해 대부분의 반복 솔버의 수렴이 보다 쉬워집니다. equilibrate를 사용하면 이후 equilibrate가 적용된 행렬 B = R*P*A*C를 분해할 때 더 나은 품질의 선조건자 행렬을 만들 수도 있습니다.

  • 계수 행렬을 인수 분해할 때 dissectsymrcm과 같은 행렬 재정렬 함수를 사용하여 계수 행렬의 행과 열을 치환하고 0이 아닌 요소의 개수를 최소화하여 선조건자를 생성할 수 있습니다. 이렇게 하면 이후 선조건이 적용된 선형 시스템을 푸는 데 필요한 메모리와 시간을 절약할 수 있습니다.

참고 문헌

[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.

확장 기능

버전 내역

R2006a 이전에 개발됨