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.