Main Content

이 페이지의 최신 내용은 아직 번역되지 않았습니다. 최신 내용은 영문으로 볼 수 있습니다.

cgs

선형 연립방정식 풀기 — 켤레 기울기 제곱법(Conjugate Gradients Squared Method)

설명

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

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

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

예제

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

예제

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

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

예제

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

예제

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

예제

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

예제

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

예제

모두 축소

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

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

rng default
A = sprand(600,600,.5);
A = A'*A + speye(size(A));
b = rand(600,1);

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

x = cgs(A,b);
cgs 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.068.

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

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

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

허용오차가 더 느슨해지고 반복 횟수를 늘렸는데도 cgs가 수렴하지 않았습니다. 반복 알고리즘이 이런 식으로 정체 상태에 빠진다는 것은 선조건자 행렬이 필요하다는 의미입니다.

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

L = ichol(A);
x = cgs(A,b,1e-4,40,L);
cgs converged at iteration 13 to a solution with relative residual 7.1e-05.

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

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

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

load west0479
A = west0479;

참(True) 해가 1로만 이루어진 벡터가 되도록 b를 정의합니다.

b = sum(A,2);

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

tol = 1e-12;
maxit = 20;

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

  • x0A*x0 = b의 계산된 해입니다.

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

  • rr0은 계산된 답 x0의 잔차입니다.

  • it0x0이 계산된 반복 횟수입니다.

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

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

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

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

setup = struct('type','ilutp','droptol',1e-6);
[L,U] = ilu(A,setup);
[x1,fl1,rr1,it1,rv1] = cgs(A,b,tol,maxit,L,U);
fl1
fl1 = 0
rr1
rr1 = 4.3851e-14
it1
it1 = 3

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

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

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')

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

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

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

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

maxit = 200;
x1 = cgs(A,b,[],maxit);
cgs converged at iteration 17 to a solution with relative residual 8.8e-07.
x0 = 0.99*ones(size(A,2),1);
x2 = cgs(A,b,[],maxit,[],[],x0);
cgs converged at iteration 4 to a solution with relative residual 8e-07.

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

중간 결과 반환하기

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

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

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

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

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

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

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의 각 행에 x의 요소를 곱하므로 단지 몇 개의 결과(삼중대각의 0이 아닌 요소에 해당함)만 0이 아니게 됩니다. 또한 주대각선에만 0과 1이 아닌 요소가 있습니다. 표현식 Ax는 다음과 같이 표현됩니다.

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

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

[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

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

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

b = ones(21,1);
tol = 1e-12;  
maxit = 50;
x1 = cgs(@afun,b,tol,maxit)
cgs converged at iteration 11 to a solution with relative residual 1.3e-14.
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
복소수 지원 여부:

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

데이터형: double

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

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

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

M을 함수 핸들로 지정

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

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

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

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

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

출력 인수

모두 축소

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

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

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

플래그 값

수렴

0

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

1

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

2

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

3

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

4

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

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

데이터형: double

반복 횟수로, 스칼라로 반환됩니다. 이 출력값은 x의 답이 계산된 반복 횟수를 나타냅니다.

데이터형: double

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

데이터형: double

세부 정보

모두 축소

켤레 기울기 제곱법(Conjugate Gradients Squared Method)

켤레 기울기 제곱(CGS) 알고리즘은 쌍켤레 기울기(BiCG) 알고리즘을 개선하기 위해 개발되었습니다. CGS 알고리즘은 잔차와 그 켤레를 사용하는 대신 잔차 제곱을 사용함으로써 계수 행렬의 전치를 사용하는 것을 피합니다[1].

CGS는 BiCG와 거의 동일한 계산 비용으로 더 빠른 수렴을 달성하지만 불규칙 수렴 동작을 보일 수 있으며 특히 초기 추측값이 해에 가까울 때 그러합니다[1].

  • 대부분의 반복법의 수렴 여부는 계수 행렬의 조건수 cond(A)에 따라 결정됩니다. equilibrate를 사용하여 A의 조건수를 개선할 수 있으며, 이로 인해 대부분의 반복 솔버의 수렴이 보다 쉬워집니다. 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.

[2] Sonneveld, Peter, “CGS: A fast Lanczos-type solver for nonsymmetric linear systems,” SIAM J. Sci. Stat. Comput., January 1989, Vol. 10, No. 1, pp. 36–52.

확장 기능

R2006a 이전에 개발됨