ilu
불완전 LU 분해
설명
[___] = ilu(는 구조체 A,options)options로 지정된 옵션으로 A에 대한 불완전 LU 분해를 수행합니다.
예를 들어, options의 type 필드를 "ilutp"로 설정하여 피벗 연산으로 불완전 LU 분해를 수행할 수 있습니다. 그런 다음 milu 필드를 "row" 또는 "col"로 설정하여 행의 합계 또는 열의 합계를 유지하는 수정된 불완전 LU 분해를 지정할 수 있습니다. 이 옵션 조합은 치환 행렬 P를 반환합니다. 여기서 "row" 옵션의 경우 L과 U는 A*P의 불완전 인수이거나(이때 U는 열 치환됨), "col" 옵션의 경우 L과 U는 P*A의 불완전 인수입니다(이때 L은 행 치환됨).
예제
ilu 함수는 3가지 유형의 불완전 LU 분해, 즉 0 채우기 분해(ILU(0)), 크라우트 버전(ILUC) 및 임계값 기각과 피벗 연산을 사용한 분해(ILUTP)를 제공합니다.
기본적으로 ilu는 희소 행렬 입력값에 대해 0 채우기 불완전 LU 분해를 수행합니다. 예를 들어, 7840개의 0이 아닌 요소를 갖는 희소 행렬에 대한 완전 및 불완전 행렬 분해를 구합니다. 완전 LU 인수는 126,478개의 0이 아닌 요소를 가지며, 0 채우기가 있는 불완전 LU 인수는 A와 같은 개수인 7840개의 0이 아닌 요소를 가집니다.
A = gallery("neumann",1600) + speye(1600);
n = nnz(A)n = 7840
n = nnz(lu(A))
n = 126478
n = nnz(ilu(A))
n = 7840
0 채우기 분해는 LU 인수에 입력 행렬의 희소성 패턴을 유지하므로, 분해의 상대 오차는 A의 0이 아닌 요소의 패턴에 대해 본질적으로 0입니다.
[L,U] = ilu(A); e = norm(A-(L*U).*spones(A),"fro")/norm(A,"fro")
e = 4.8874e-17
그러나 이러한 0 채우기 인수의 곱은 원래 행렬의 적절한 근사가 아닙니다.
e = norm(A-L*U,"fro")/norm(A,"fro")
e = 0.0601
정확도를 높이기 위해 필인(fill-in)이 있는 다른 유형의 불완전 LU 분해를 사용할 수 있습니다. 예를 들어, 1e-6 기각 허용오차와 함께 크라우트 버전을 사용합니다.
options.droptol = 1e-6;
options.type = "crout";
[L,U] = ilu(A,options);크라우트 버전의 불완전 행렬 분해는 LU 인수에 (완전 행렬 분해보다 적은) 51,482개의 0이 아닌 요소를 갖습니다. 필인을 사용하면 불완전 LU 인수의 곱이 원래 행렬의 더 적절한 근사가 됩니다.
n = nnz(ilu(A,options))
n = 51482
e = norm(A-L*U,"fro")./norm(A,"fro")
e = 9.3040e-07
비교를 위해, 동일한 입력 행렬 A에 대해 임계값 기각과 피벗 연산을 사용한 불완전 행렬 분해를 수행할 경우 크라우트 버전과 유사한 결과가 산출됩니다.
options.type = "ilutp";
[L,U,P] = ilu(A,options);
n = nnz(ilu(A,options))n = 51541
norm(P*A-L*U,"fro")./norm(A,"fro")
ans = 9.4960e-07
희소 행렬을 인수 분해하도록 불완전 LU 분해의 기각 허용오차를 변경합니다.
실수 값의 479×479 비대칭 희소 행렬인 west0479 행렬을 불러옵니다. condest를 사용하여 행렬의 조건수를 추정합니다.
load west0479
A = west0479;
c1 = condest(A)c1 = 1.4244e+12
equilibrate를 사용하여 행렬의 조건수를 개선합니다. 원래 행렬 A를 치환하고 다시 스케일링하여 새 행렬 B = R*P*A*C를 만듭니다. 이 행렬은 더 나은 조건수를 가지며 대각선 요소가 1과 –1로만 구성되어 있습니다.
[P,R,C] = equilibrate(A); B = R*P*A*C; c2 = condest(B)
c2 = 5.1036e+04
B에 대해 행의 합계를 유지하는 임계값 기각과 피벗 연산을 사용한 불완전 LU 분해를 수행하는 옵션을 지정합니다. 비교를 위해, 먼저 필인의 기각 허용오차를 0으로 지정합니다. 이 경우 완전 LU 분해가 생성됩니다.
options.type = "ilutp"; options.milu = "row"; options.droptol = 0; [L,U,P] = ilu(B,options);
이 분해는 입력 행렬 B를 매우 정확하게 근사하지만, 인수는 B보다 훨씬 더 조밀합니다.
e = norm(B*P-L*U,"fro")/norm(B,"fro")
e = 1.0345e-16
nLU = nnz(L)+nnz(U)-size(B,1)
nLU = 19118
nB = nnz(B)
nB = 1887
기각 허용오차를 변경하여 더 희소하지만 B를 근사하는 데 덜 정확한 불완전 LU 인수를 얻을 수 있습니다. 예를 들어, 다음 플롯은 기각 허용오차에 대해 플로팅된 입력 행렬의 밀도에 대한 불완전 인수의 밀도의 비율과 불완전 행렬 분해의 상대 오차를 보여줍니다.
ntols = 20; tau = logspace(-6,-2,ntols); e = zeros(1,ntols); nLU = zeros(1,ntols); for k = 1:ntols options.droptol = tau(k); [L,U,P] = ilu(B,options); nLU(k) = nnz(L)+nnz(U)-size(B,1); e(k) = norm(B*P-L*U,"fro")/norm(B,"fro"); end figure semilogx(tau,nLU./nB,LineWidth=2) title("Ratio of nonzeros of LU factors with respect to B") xlabel("drop tolerance") ylabel("nnz(L)+nnz(U)-size(B,1)/nnz(B)",FontName="FixedWidth")
figure loglog(tau,e,LineWidth=2) title("Relative error of the incomplete LU factorization") xlabel("drop tolerance") ylabel("$||BP-LU||_F\,/\,||B||_F$",Interpreter="latex")
이 예제에서 임계값 기각을 사용한 불완전 LU 분해의 상대 오차는 기각 허용오차와 동일한 차수(항상 동일하지 않음)입니다.
불완전 LU 분해를 bicgstab에 대한 선조건자로 사용하여 선형 시스템을 풉니다.
479×479 실수 비대칭 희소 행렬인 west0479를 불러옵니다.
load west0479
A = west0479;의 참(True) 해가 1로만 구성된 벡터가 되도록 b를 정의합니다.
b = full(sum(A,2));
허용오차와 최대 반복 횟수를 설정합니다.
tol = 1e-12; maxit = 20;
bicgstab를 사용하여 요청된 허용오차와 반복 횟수로 해를 구합니다. 풀이 과정에 대한 정보를 반환하는 5개의 출력값을 지정합니다.
x는A*x = b의 계산된 해입니다.fl0은 알고리즘의 수렴 여부를 나타내는 플래그입니다.rr0은 계산된 답x의 상대 잔차입니다.it0은x가 계산된 반복 횟수입니다.rv0은 에 대한 각 절반의 반복에서 잔차 내역으로 구성된 벡터입니다.
[x,fl0,rr0,it0,rv0] = bicgstab(A,b,tol,maxit); fl0
fl0 = 1
rr0
rr0 = 1
it0
it0 = 0
bicgstab가 요청된 반복 횟수인 20회 이내에 요청된 허용오차 1e-12로 수렴하지 않으므로 fl0은 1이 됩니다. 사실상 bicgstab의 동작이 좋지 않기 때문에 초기 추측값 x0 = zeros(size(A,2),1)이 최적의 해가 되고 it0 = 0에서 볼 수 있듯이 이 해가 반환됩니다.
이러한 느린 수렴을 개선할 수 있도록 선조건자 행렬을 지정할 수 있습니다. A는 비대칭 행렬이므로 ilu를 사용하여 선조건자 를 생성합니다. 대각선상에 있지 않으면서 값이 1e-6보다 작은 요소는 무시하도록 기각 허용오차를 지정합니다. L과 U를 bicgstab에 대한 입력값으로 지정하여 선조건이 적용된 시스템 를 풉니다.
options = struct("type","ilutp","droptol",1e-6); [L,U] = ilu(A,options); [x_precond,fl1,rr1,it1,rv1] = bicgstab(A,b,tol,maxit,L,U); fl1
fl1 = 0
rr1
rr1 = 5.9892e-14
it1
it1 = 3
ilu 선조건자를 사용한 결과 3번째 반복에서 요청된 허용오차 1e-12보다 작은 상대 잔차 rr1이 생성됩니다. 출력값 rv1(1)은 norm(b)가 되고, 출력값 rv1(end)는 norm(b-A*x1)이 됩니다.
각 절반의 반복에서 상대 잔차를 플로팅하면 bicgstab의 진행률을 추적할 수 있습니다. 각 해의 잔차 내역을 플로팅하고, 지정된 허용오차에 대한 선을 추가합니다.
semilogy((0:numel(rv0)-1)/2,rv0/norm(b),"-o") hold on semilogy((0:numel(rv1)-1)/2,rv1/norm(b),"-o") yline(tol,"r--"); legend("No preconditioner","ILU preconditioner","Tolerance",Location="East") xlabel("Iteration number") ylabel("Relative residual")
입력 인수
출력 인수
참고 문헌
[1] Saad, Y. "Preconditioning Techniques", chap. 10 in Iterative Methods for Sparse Linear Systems. The PWS Series in Computer Science. Boston: PWS Pub. Co, 1996.
