Main Content

ldl

에르미트 부정부호 행렬(Hermitian Indefinite Matrix)의 블록 LDL 분해

구문

L = ldl(A)
[L,D] = ldl(A)
[L,D,P] = ldl(A)
[L,D,p] = ldl(A,'vector')
[U,D,P] = ldl(A,'upper')
[U,D,p] = ldl(A,'upper','vector')
[L,D,P,S] = ldl(A)
[L,D,P,S] = LDL(A,THRESH)
[U,D,p,S] = LDL(A,THRESH,'upper','vector')

설명

L = ldl(A)는 2출력 형식에서와 같이 치환된 하부 삼각 행렬 L만 반환합니다. 블록 대각선 인수 D처럼 치환 정보가 손실됩니다. 기본적으로 ldl은 A의 대각 행렬과 하부 삼각 행렬만 참조하며, 상부 삼각 행렬이 하부 삼각 행렬의 켤레 복소수 전치라고 가정합니다. 따라서 [L,D,P] = ldl(TRIL(A))와 [L,D,P] = ldl(A)는 모두 정확히 같은 인수를 반환합니다. 참고로, 이 구문은 희소 행렬 A에는 유효하지 않습니다.

[L,D] = ldl(A)는 블록 대각 행렬(Block Diagonal Matrix) D와 치환된 하부 삼각 행렬을 L에 저장하여 A = L*D*L'이 되도록 합니다. 블록 대각 행렬 D는 대각선에 1×1 블록과 2×2 블록을 포함합니다. 참고로, 이 구문은 희소 행렬 A에는 유효하지 않습니다.

[L,D,P] = ldl(A)는 단위 하부 삼각 행렬 L, 블록 대각 행렬 D와 치환 행렬 P를 반환하여 P'*A*P = L*D*L'이 되도록 합니다. 이는 [L,D,P] = ldl(A,'matrix')와 동일합니다.

[L,D,p] = ldl(A,'vector')는 치환 정보를 행렬이 아닌 벡터 p로 반환합니다. p 출력값은 A(p,p) = L*D*L'을 충족하는 행 벡터입니다.

[U,D,P] = ldl(A,'upper')는 A의 대각 행렬과 상부 삼각 행렬만 참조하며, 하부 삼각 행렬이 상부 삼각 행렬의 켤레 복소수 전치라고 가정합니다. 이 구문은 단위 상부 삼각 행렬 U를 반환하여 P'*A*P = U'*D*U가 되도록 하며, A를 상부 삼각 행렬이 아닌 에르미트 행렬으로 가정합니다. 마찬가지로, [L,D,P] = ldl(A,'lower')도 디폴트 동작을 제공합니다.

[L,D,p] = ldl(A,'lower','vector')에서처럼 [U,D,p] = ldl(A,'upper','vector')도 치환 정보를 벡터 p로 반환합니다. A는 비희소 행렬(Full Matrix)이어야 합니다.

[L,D,P,S] = ldl(A)는 단위 하부 삼각 행렬 L, 블록 대각 행렬 D, 치환 행렬 P 및 스케일링 행렬 S를 반환하여 P'*S*A*S*P = L*D*L'이 되도록 합니다. 이 구문은 실수 희소 행렬에만 사용할 수 있으며, A의 하부 삼각 행렬만 참조합니다.

[L,D,P,S] = LDL(A,THRESH)는 THRESH를 알고리즘의 피벗 허용오차로 사용합니다. THRESH는 구간 [0, 0.5]에 있는 double형 스칼라여야 합니다. THRESH의 디폴트 값은 0.01입니다. 더 작은 THRESH 값을 사용하면 행렬 분해 시간이 빨라지고 요소가 적어질 수 있지만, 행렬 분해의 안정성이 저하될 수도 있습니다. 이 구문은 실수 희소 행렬에만 사용할 수 있습니다.

[U,D,p,S] = LDL(A,THRESH,'upper','vector')는 위에서 설명된 대로 피벗 허용오차를 설정하고 상부 삼각 행렬 U와 치환 벡터 p를 반환합니다.

예제

다음 예제에서는 1출력, 2출력, 3출력 형식을 포함하여 ldl 함수의 다양한 형식을 사용하는 방법과 vector 및 upper 옵션의 사용 방법을 보여줍니다. 여기에서 다루는 항목은 다음과 같습니다.

이 예제들을 실행하기 전에 다음과 같은 양의 정부호 행렬과 부정부호 에르미트 행렬을 생성해야 합니다.

A = full(delsq(numgrid('L', 10)));
B = gallery('uniformdata',10,0);
M = [eye(10) B; B' zeros(10)];

여기에서 M의 구조는 최적화 및 유체 흐름 문제에서 매우 일반적이며, M은 사실상 부정부호 행렬입니다. 참고로, ldl은 희소 인수를 받지 않으므로 양의 정부호 행렬 A는 비희소 행렬이어야 합니다.

예제 1 — ldl의 2출력 형식

ldl의 2출력 형식은 L과 D를 반환하여, A-(L*D*L')이 작아지고 L이 치환된 단위 하부 삼각 행렬이 되며 D가 2×2 블록 대각 행렬이 되도록 합니다. 참고로, A는 양의 정부호 행렬이므로 D의 대각선은 모두 양수입니다.

[LA,DA] = ldl(A); 
fprintf(1, ...
'The factorization error ||A - LA*DA*LA''|| is %g\n', ...
norm(A - LA*DA*LA'));
neginds = find(diag(DA) < 0)

a b가 주어지면, LA, DA를 사용하여 Ax=b를 풉니다.

bA = sum(A,2);
x = LA'\(DA\(LA\bA));
fprintf(...
'The absolute error norm ||x - ones(size(bA))|| is %g\n', ...
norm(x - ones(size(bA))));

예제 2 — ldl의 3출력 형식

3출력 형식은 치환 행렬도 반환하므로, L은 사실상 단위 하부 삼각 행렬이 됩니다. 

[Lm, Dm, Pm] = ldl(M);
fprintf(1, ...
'The error norm ||Pm''*M*Pm - Lm*Dm*Lm''|| is %g\n', ...
norm(Pm'*M*Pm - Lm*Dm*Lm'));
fprintf(1, ...
'The difference between Lm and tril(Lm) is %g\n', ...
norm(Lm - tril(Lm)));

b가 주어지면, Lm, Dm 및 Pm을 사용하여 Mx=b를 풉니다.

bM = sum(M,2);
x = Pm*(Lm'\(Dm\(Lm\(Pm'*bM))));
fprintf(...
'The absolute error norm ||x - ones(size(b))|| is %g\n', ...
norm(x - ones(size(bM))));

예제 3 — D의 구조

D는 1×1 블록과 2×2 블록이 있는 블록 대각 행렬입니다. 따라서 삼중대각 행렬(Tridiagonal Matrix)의 특수한 경우가 됩니다. 입력 행렬이 양의 정부호 행렬이면 D는 거의 항상 대각 행렬입니다(행렬의 정부호 정도에 따라 다름). 그러나 행렬이 부정부호 행렬이면 D는 대각 행렬이거나 블록 구조를 표현할 수 있습니다. 예를 들어, A가 위와 같을 때 DA는 대각 행렬입니다. 그러나 A를 조금만 이동하면 부정부호 행렬이 되며, 블록 구조를 갖는 D를 계산할 수 있습니다.

figure; spy(DA); title('Structure of D from ldl(A)');
[Las, Das] = ldl(A - 4*eye(size(A)));
figure; spy(Das); 
title('Structure of D from ldl(A - 4*eye(size(A)))');

예제 4 — 'vector' 옵션 사용

lu 함수와 같이, ldl은 함수가 치환 벡터 또는 치환 행렬을 반환하는지를 확인하는 인수를 받습니다. ldl은 기본적으로 치환 행렬을 반환합니다. 'vector'를 선택하면 함수가 더 빨리 실행되며 메모리를 더 적게 사용합니다. 이런 이유로, 'vector' 옵션을 지정하는 것이 좋습니다. 참고로, 이러한 연산에서는 인덱싱이 곱셈보다 일반적으로 더 빠릅니다.

[Lm, Dm, pm] = ldl(M, 'vector');
fprintf(1, 'The error norm ||M(pm,pm) - Lm*Dm*Lm''|| is %g\n', ...
  norm(M(pm,pm) - Lm*Dm*Lm'));

% Solve a system with this kind of factorization.
clear x;
x(pm,:) = Lm'\(Dm\(Lm\(bM(pm,:))));
fprintf('The absolute error norm ||x - ones(size(b))|| is %g\n', ...
  norm(x - ones(size(bM))));

예제 5 — 'upper' 옵션 사용

chol 함수와 같이, ldl은 입력 행렬의 어느 삼각형이 참조되는지와, ldl이 하부(L) 또는 상부(L') 삼각 인수 중 어느 것을 반환하는지를 확인하는 인수를 받습니다. 조밀 행렬의 경우 하부 삼각 버전이 아닌 상부 삼각 버전을 사용한다고 해서 실제적으로 얻는 이득이 없습니다. 

Ml = tril(M);
[Lml, Dml, Pml] = ldl(Ml, 'lower'); % 'lower' is default behavior.
fprintf(1, ...
'The difference between Lml and Lm is %g\n', norm(Lml - Lm));
[Umu, Dmu, pmu] = ldl(triu(M), 'upper', 'vector');
fprintf(1, ...
'The difference between Umu and Lm'' is %g\n', norm(Umu - Lm'));

% Solve a system using this factorization.
clear x;
x(pm,:) = Umu\(Dmu\(Umu'\(bM(pmu,:))));
fprintf(...
'The absolute error norm ||x - ones(size(b))|| is %g\n', ...
norm(x - ones(size(bM))));

'upper' 옵션과 'vector' 옵션을 둘 다 지정할 경우 인수 목록에서 'upper'가 'vector'보다 앞에 나와야 합니다.

예제 6 — linsolve와 에르미트 부정형 솔버

linsolve 함수 사용 시, 시스템에 대칭 행렬이 있다는 것을 이용하여 향상된 성능을 얻을 수도 있습니다. 위 예제에서 사용된 행렬은 성능 향상의 정도를 느끼기에는 작기 때문에, 이 예제에서는 더 큰 행렬을 생성해 봅니다. 이 행렬은 양의 정부호 대칭 행렬이며, 아래를 보면 행렬에 대해 알고 있는 각 사항에 따라 그에 대응하는 속도 향상이 있다는 것을 알 수 있습니다. 즉, 대칭 솔버는 일반적인 솔버보다 더 빠른 반면, 양의 정부호 대칭 솔버는 대칭 솔버보다 더 빠릅니다.

Abig = full(delsq(numgrid('L', 30)));
bbig = sum(Abig, 2);
LSopts.POSDEF = false;
LSopts.SYM = false;
tic; linsolve(Abig, bbig, LSopts); toc;
LSopts.SYM = true;
tic; linsolve(Abig, bbig, LSopts); toc;
LSopts.POSDEF = true;
tic; linsolve(Abig, bbig, LSopts); toc;

참고 문헌

[1] Ashcraft, Cleve, Roger G. Grimes, and John G. Lewis. “Accurate Symmetric Indefinite Linear Equation Solvers.” SIAM Journal on Matrix Analysis and Applications 20, no. 2 (January 1998): 513–561. https://doi.org/10.1137/S0895479896296921.

[2] Anderson, E., ed. LAPACK Users’ Guide. 3rd ed. Software, Environments, Tools. Philadelphia: Society for Industrial and Applied Mathematics, 1999. https://doi.org/10.1137/1.9780898719604.

[3] Duff, Iain S. “MA57---a Code for the Solution of Sparse Symmetric Definite and Indefinite Systems.” ACM Transactions on Mathematical Software 30, no. 2 (June 2004): 118–144. https://doi.org/10.1145/992200.992202.

확장 기능

참고 항목

| |