희소 행렬의 선형대수

희소 행렬의 선형대수 (Linear Algebra on Sparse Matrices)

Octave에는 희소 행렬을 위한 다형(polymorphic) 솔버가 들어 있어요. 행렬을 인수분해할 때 정확히 어떤 솔버를 쓸지는 희소 행렬 자체의 특성에 따라 결정됩니다. 대체로 행렬 타입을 알아내는 비용은 행렬을 인수분해하는 비용에 비해 작지만, 어쨌든 타입은 한 번 계산되면 캐시돼서 선형 방정식에 쓸 때마다 다시 판정하지 않아요. Octave 매뉴얼 22.2절을 정리해 볼게요.

출처: Linear Algebra on Sparse Matrices

본문

선형 방정식을 푸는 선택 트리는 다음과 같아요.

행렬이 대각(diagonal)이면 직접 풀고 8로 간다.
행렬이 치환된 대각(permuted diagonal)이면 치환을 고려해 직접 풀고 8로 간다.
행렬이 정사각이면서 띠(banded)이고 띠 밀도가 spparms ("bandden")이 준 값보다 작으면 계속 진행하고, 아니면 4로 간다.
행렬이 삼대각(tridiagonal)이고 우변(right-hand side)이 희소가 아니면 계속 진행하고, 아니면 3b로 간다.
행렬이 Hermitian이고 양의 실수 대각을 가지면 LAPACK xPTSV로 Cholesky 인수분해를 시도한다.
위가 실패했거나 행렬이 양의 실수 대각을 가진 Hermitian이 아니면, 피보팅을 포함한 가우스 소거법(LAPACK xGTSV)를 쓰고 8로 간다.
행렬이 양의 실수 대각을 가진 Hermitian이면 LAPACK xPBTRF로 Cholesky 인수분해를 시도한다.
위가 실패했거나 행렬이 양의 실수 대각을 가진 Hermitian이 아니면, 피보팅을 포함한 가우스 소거법(LAPACK xGBTRF)을 쓰고 8로 간다.
행렬이 상삼각 또는 하삼각이면 희소 전방/후방 대입을 수행하고 8로 간다.
행렬이 열 치환이 있는 상삼각 또는 행 치환이 있는 하삼각이면 희소 전방/후방 대입을 수행하고 8로 간다.
행렬이 정사각이고, 실수 양의 대각을 가진 Hermitian이면 CHOLMOD로 희소 Cholesky 인수분해를 시도한다.
희소 Cholesky 인수분해가 실패했거나 행렬이 실수 양의 대각을 가진 Hermitian이 아니고 정사각이면, UMFPACK으로 인수분해·해를 구하고 정제(refinement)를 한 번 수행한다.
행렬이 정사각이 아니거나 이전 솔버 중 하나가 특이 또는 특이에 가까운 행렬을 표시하면 CXSPARSE (Footnote 10)로 최소 노름(minimum norm) 해를 찾는다.

띠 밀도는 띠 안의 0이 아닌 값의 개수를 전체 띠의 값 개수로 나눈 값이에요. 띠 행렬 솔버는 spparms로 bandden을 1로 설정하면 완전히 끌 수 있어요 (즉 spparms ("bandden", 1)).

QR 솔버는 문제를 Dulmage-Mendelsohn 분해로 인수분해해, 과결정(over-determined) 블록 여러 개와 마지막 과결정 블록으로 처리할 수 있는 블록들로 나눠요. 강하게 연결된 노드 블록이 있는 행렬에서는 많은 블록에 LU 분해를 쓸 수 있어 큰 이점이 됩니다. 또한 NaN 벡터를 그냥 돌려주는 대신 과결정 문제의 해를 찾을 가능성을 크게 높여줘요.

위 솔버들은 모두 조건수(condition number)의 추정값을 계산할 수 있어요. 이는 해에서 수치 안정성 문제를 감지하고 최소 노름 해를 강제하는 데 쓸 수 있죠. 다만 좁은 띠, 삼각, 대각 행렬의 경우 조건수를 계산하는 비용이 크고 실제로 행렬을 인수분해하는 비용을 초과할 수도 있어요. 그래서 이런 경우에는 조건수를 계산하지 않고, Octave는 특이 행렬을 감지하는 더 단순한 기법이나(띠 행렬의 경우) 내부 LAPACK 코드에 의존합니다.

사용자는 matrix_type 함수로 행렬 타입을 강제할 수 있어요. 이러면 행렬 타입을 알아내는 비용을 피할 수 있죠. 다만 행렬 타입을 잘못 지정하면 예측 불가능한 결과가 나올 수 있으니, matrix_type은 신중히 사용해야 합니다.

normest

nest = normest (A)
nest = normest (A, tol)
[nest, iter] = normest (...)

멱급수(power series) 분석으로 행렬 A의 2-노름을 추정해요. norm (A)를 계산하는 비용이 너무 커서 2-노름의 근사값으로 충분한 큰 행렬에 주로 쓰입니다. tol은 2-노름을 계산할 허용 오차로, 기본값은 1e-6이에요. 선택적 출력 iter는 normest가 수렴하는 데 필요한 반복 횟수를 돌려줍니다.

See also: normest1, norm, cond, condest.

normest1

nest = normest1 (A)
nest = normest1 (A, t)
nest = normest1 (A, t, x0)
nest = normest1 (Afcn, t, x0, p1, p2, …)
[nest, v] = normest1 (A, …)
[nest, v, w] = normest1 (A, …)
[nest, v, w, iter] = normest1 (A, …)

블록 알고리즘으로 행렬 A의 1-노름을 추정해요. 노름의 추정값만 필요한 크고 희소한 행렬에 최적입니다. 작거나 중간 크기 행렬은 norm (A, 1)을 고려하세요. 또한 행렬-벡터 곱 A * x와 A' * x를 저렴하게 계산할 수 있을 때 선형 연산자 A의 1-노름 추정에도 쓸 수 있어요. 이 경우 행렬 A 대신 함수 Afcn (flag, x)를 쓰는데, 이 함수는 다음을 돌려줘야 해요:

flag가 "dim"이면 A의 차원 n
flag가 "real"이면 A가 실수 연산자인지 여부
flag가 "notransp"이면 A * x의 결과
flag가 "transp"이면 A' * x의 결과

대표적인 예는 A가 b ^ m으로 정의되는 경우예요. 이때 결과 A * x를 b ^ m을 명시적으로 만들지 않고도 이렇게 계산할 수 있어요:

y = x;
for i = 1:m
  y = b * y;
endfor

매개변수 p1, p2, …는 Afcn (flag, x, p1, p2, …)의 인자예요. t의 기본값은 2이고, 알고리즘은 n×n과 n×t 크기의 행렬-행렬 곱을 필요로 해요.

초기 행렬 x0는 열의 1-노름이 1이어야 해요. 기본 초기 행렬 x0는 첫 번째 열이 ones (n, 1) / n이고, t > 1이면 나머지 열은 -1 / n, 1 / n의 무작위 원소를 n으로 나눈 값이에요.

출력에서 nest는 원하는 추정값이고, v와 w는 w = A * v를 만족하는 벡터예요 (이때 norm (w, 1) = c * norm (v, 1)). iter는 iter(1)에 반복 횟수(최대는 5로 하드코딩), iter(2)에 알고리즘이 수행한 A * x 또는 A' * x 곱의 총 횟수를 담아요.

Algorithm Note: normest1은 평가 중 무작위 수를 사용해요. 따라서 일관된 결과가 필요하면 normest1을 호출하기 전에 난수 생성기의 "state"를 고정해야 해요.

Reference: N. J. Higham and F. Tisseur, "A block algorithm for matrix 1-norm estimation, with and application to 1-norm pseudospectra", SIAM J. Matrix Anal. Appl., Vol. 21, No. 4, pp. 1185–1201, 2000.

See also: normest, norm, cond, condest.

condest

cest = condest (A)
cest = condest (A, t)
cest = condest (A, Ainvfcn)
cest = condest (A, Ainvfcn, t)
cest = condest (A, Ainvfcn, t, p1, p2, …)
cest = condest (Afcn, Ainvfcn)
cest = condest (Afcn, Ainvfcn, t)
cest = condest (Afcn, Ainvfcn, t, p1, p2, …)
[cest, v] = condest (...)

t개의 테스트 벡터와 무작위 1-노름 추정기를 사용해 정사각 행렬 A의 1-노름 조건수를 추정해요. 선택적 입력 t는 테스트 벡터 개수(기본 5)를 지정해요.

입력은 행렬 A일 수 있는데, 알고리즘은 특히 크고 희소한 행렬에 적합합니다. 또는 함수로 행렬의 동작을 암시적으로 정의할 수도 있어요. 암시적 정의를 쓸 때 condest는 다음 함수들을 필요로 해요:

Afcn (flag, x) — 다음을 돌려줘야 함
- flag가 "dim"이면 A의 차원 n
- flag가 "real"이면 A가 실수 연산자인지 여부
- flag가 "notransp"이면 A * x의 결과
- flag가 "transp"이면 A' * x의 결과
Ainvfcn (flag, x) — 다음을 돌려줘야 함
- flag가 "dim"이면 inv (A)의 차원 n
- flag가 "real"이면 inv (A)가 실수 연산자인지 여부
- flag가 "notransp"이면 inv (A) * x의 결과
- flag가 "transp"이면 inv (A)' * x의 결과

매개변수 p1, p2, …는 Afcn (flag, x, p1, p2, …)와 Ainvfcn (flag, x, p1, p2, …)의 추가 인자예요.

주된 출력은 1-노름 조건수 추정값 cest이에요. 선택적 두 번째 출력 v는 근사적으로 영 공간(null space)에 있는 벡터로, norm (A*v, 1) == norm (A, 1) * norm (v, 1) / cest를 만족해요.

Algorithm Note: condest는 1-노름을 근사할 때 무작위 알고리즘을 사용해요. 따라서 일관된 결과가 필요하면 "state"를 고정해야 합니다.

References: N.J. Higham and F. Tisseur, "A Block Algorithm for Matrix 1-Norm Estimation, with an Application to 1-Norm Pseudospectra", SIAM Journal on Matrix Analysis and Applications, Vol. 21, Iss. 4, pp. 1185–1201, 2000, https://dx.doi.org/10.1137/S0895479899356080.

See also: cond, rcond, norm, normest1, normest.

spparms

spparms ()
vals = spparms ()
[keys, vals] = spparms ()
val = spparms (key)
spparms (vals)
spparms ("default")
spparms ("tight")
spparms (key, val)

희소 솔버와 인수분해 함수가 쓰는 매개변수를 조회하거나 설정해요. 위에 있는 처음 네 호출은 현재 설정에 대한 정보를 얻고, 나머지는 현재 설정을 바꿉니다. 매개변수는 키-값 쌍으로 저장되는데, 값은 모두 실수(float)이고 키는 다음 문자열 중 하나예요:

'spumoni' — 솔버의 디버깅 정보 출력 수준 (기본 0)
'ths_rel' — 호환성을 위해 포함됨. 사용 안 함. (기본 1)
'ths_abs' — 호환성을 위해 포함됨. 사용 안 함. (기본 1)
'exact_d' — 호환성을 위해 포함됨. 사용 안 함. (기본 0)
'supernd' — 호환성을 위해 포함됨. 사용 안 함. (기본 3)
'rreduce' — 호환성을 위해 포함됨. 사용 안 함. (기본 3)
'wh_frac' — 호환성을 위해 포함됨. 사용 안 함. (기본 0.5)
'autommd' — LU/QR와 '\', '/' 연산자가 희소성을 보존하는 mmd 함수를 자동으로 사용할지 (기본 1)
'autoamd' — LU와 '\', '/' 연산자가 희소성 보존 amd 함수를 자동으로 사용할지 (기본 1)
'piv_tol' — UMFPACK 솔버의 피벗 허용오차 (기본 0.1)
'sym_tol' — UMFPACK 대칭 솔버의 피벗 허용오차 (기본 0.001)
'bandden' — 띠 행렬이 LAPACK 띠 솔버로 처리되기 전의 0이 아닌 원소 밀도 (기본 0.5)
'umfpack' — LU, '\', '/' 연산에 UMFPACK을 쓸지 mmd 솔버를 쓸지 (기본 1)

개별 키의 값은 spparms (key, val)로 설정할 수 있어요. 기본값은 특별 키워드 "default"로 복원할 수 있고, 특별 키워드 "tight"를 쓰면 mmd 솔버가 더 희소한 해를 시도하도록 할 수 있는데 잠재적으로 실행 시간이 길어질 수 있어요.

See also: chol, colamd, lu, qr, symamd.

sprank

p = sprank (S)

희소 행렬 S의 구조적 계수(structural rank)를 계산해요. 이 계산은 Dulmage-Mendelsohn 치환을 통한 블록 삼각 형태에 기반한, 행렬의 구조만 사용한다는 점에 주의하세요. 따라서 행렬 S의 수치적 계수는 sprank (S) >= rank (S)로 제한됩니다. 부동소수점 오류를 무시하면 sprank (S) == rank (S)예요.

symbfact

[count, h, parent, post, R] = symbfact (S)
[…] = symbfact (S, typ)
[…] = symbfact (S, typ, mode)

희소 행렬 S의 기호적(symbolic) 인수분해 분석을 수행해요. 입력 변수는:

S — 실수 또는 복소수 희소 행렬.
typ — 인수분해 타입. 다음 중 하나:
- "sym" (기본) — S를 인수분해. S가 대칭이라 가정하고 상삼각 부분을 사용.
- "col" — S' * S를 인수분해.
- "row" — S * S'를 인수분해.
- "lo" — S'를 인수분해. S가 대칭이라 가정하고 하삼각 부분을 사용.
mode — mode를 지정하지 않으면 R에 Cholesky 인수분해가 반환돼요. mode가 "lower" 또는 "L"이면 하삼각 인자인 켤레 전치 R'가 반환됩니다. 켤레 전치 버전은 더 빠르고 메모리를 덜 쓰지만, 다른 출력(count, h, parent, post)에는 동일한 값을 돌려줘요.

출력 변수는:

count — typ에 따라 결정된 Cholesky 인수분해의 행 개수. chol로 실제 인수분해를 수행했을 때의 계산 난이도는 sum (count .^ 2)예요.
h — 제거 트리(elimination tree)의 높이.
parent — 제거 트리 자체.
post — typ에 따라 결정된 Cholesky 인수분해의 구조를 가진 희소 불리언 행렬.

See also: chol, etree, treelayout.

정사각이 아닌 행렬의 경우 사용자는 spaugment 함수로 선형 방정식의 최소제곱 해를 찾을 수도 있어요.

spaugment

s = spaugment (A, c)

A의 확대 행렬(augmented matrix)을 만들어요. 다음으로 주어집니다:

[c * eye(m, m), A;
            A', zeros(n, n)]

이것은 A \ b의 최소제곱 해와 관련돼 있어요. 다음 식으로요:

s * [ r / c; x] = [ b, zeros(n, columns(b)) ]

여기서 r은 잔차 오류예요:

r = b - A * x

행렬 s는 대칭 부정부호(indefinite)이므로 lu로 인수분해할 수 있고, 따라서 qr 인수분해 없이도 최소 노름 해를 찾을 수 있어요. 잔차 오류는 과소결정(underdetermined) 문제에서 zeros (m, m)이 될 거예요. 예시를 보면:

m = 11; n = 10; mn = max (m, n);
A = spdiags ([ones(mn,1), 10*ones(mn,1), -ones(mn,1)],
             [-1, 0, 1], m, n);
x0 = A \ ones (m,1);
s = spaugment (A);
[L, U, P, Q] = lu (s);
x1 = Q * (U \ (L \ (P  * [ones(m,1); zeros(n,1)])));
x1 = x1(end - n + 1 : end);

과결정 문제의 해를 찾으려면 잔차 오류 r의 추정값이 필요하므로, spaugment 함수로 최소 노름 해를 만드는 게 더 복잡해요.

일반적으로 좌나눗셈 연산자가 spaugment 함수를 쓰는 것보다 더 안정적이고 빠릅니다.

eigs

d = eigs (A)
d = eigs (A, k)
d = eigs (A, k, sigma)
d = eigs (A, k, sigma, opts)
d = eigs (A, B)
d = eigs (A, B, k)
d = eigs (A, B, k, sigma)
d = eigs (A, B, k, sigma, opts)
d = eigs (Af, n)
d = eigs (Af, n, k)
d = eigs (Af, n, k, sigma)
d = eigs (Af, n, k, sigma, opts)
d = eigs (Af, n, B)
d = eigs (Af, n, B, k)
d = eigs (Af, n, B, k, sigma)
d = eigs (Af, n, B, k, sigma, opts)
[V, D] = eigs (...)
[V, D, flag] = eigs (...)

선택 기준에 따라 제한된 개수의 고유값과 고유벡터를 계산해요. 기본적으로 eigs는 표준 고유값 방정식을 풀고, 양의 정부호 행렬 B가 주어지면 일반 고유값 방정식을 풀어요.

입력 A는 n×n 차원의 정사각 행렬이에요. 대개 A는 크고 희소합니다. 일반 고유값 문제의 입력 B는 A와 같은 크기(n×n)의 정사각 행렬이에요. 대개 B도 크고 희소하죠.

계산할 고유값과 고유벡터의 개수는 k로 주어지며 기본값은 6이에요.

인자 sigma는 어떤 고유값을 반환할지 결정해요. sigma는 스칼라일 수도 있고 문자열일 수도 있어요. 스칼라일 때는 sigma에 가장 가까운 k개의 고유값이 반환됩니다. 문자열일 때는 다음 값 중 하나여야 해요:

"lm" — Largest Magnitude (기본값).
"sm" — Smallest Magnitude.
"la" — Largest Algebraic (실수 대칭 문제에서만 유효).
"sa" — Smallest Algebraic (실수 대칭 문제에서만 유효).
"be" — Both Ends, k가 홀수면 상위 쪽에서 하나 더 (실수 대칭 문제에서만 유효).
"lr" — Largest Real part (복소수 또는 비대칭 문제에서만 유효).
"sr" — Smallest Real part (복소수 또는 비대칭 문제에서만 유효).
"li" — Largest Imaginary part (복소수 또는 비대칭 문제에서만 유효).
"si" — Smallest Imaginary part (복소수 또는 비대칭 문제에서만 유효).

opts가 주어지면 eigs가 쓸 수 있는 선택적 옵션을 정의하는 구조체예요. opts 구조체의 필드는:

issym — Af가 주어지면 이 플래그(true/false)가 함수 Af가 대칭 문제를 정의하는지 결정. 행렬 A가 주어지면 무시. 기본 false.
isreal — Af가 주어지면 이 플래그가 함수 Af가 실수 문제를 정의하는지 결정. 행렬 A가 주어지면 무시. 기본 true.
tol — 필요한 수렴 허용오차를 정의하며 tol * norm (A)로 계산. 기본 eps.
maxit — 최대 반복 횟수. 기본 300.
p — 사용할 Lanczos 기저 벡터의 수. 벡터가 많을수록 수렴이 빠르지만 메모리를 더 씀. 최적 값은 문제에 따라 달라지며 k + 1에서 n 사이여야 함. 기본값은 2 * k.
v0 — 알고리즘의 시작 벡터. 최종 벡터에 가까운 초기 벡터는 수렴을 가속함. 기본은 ARPACK이 시작 벡터를 무작위로 생성. 지정되면 v0는 n = rows (A)인 n×1 벡터여야 함.
disp — 진단 출력 수준 (0|1|2). disp가 0이면 진단 비활성화. 기본 0.
cholB — 일반 고유값 문제를 계산한다면 이 플래그가 입력 B가 chol (B)를 나타내는지 그냥 행렬 B인지 지정. 기본 false.
permB — cholB가 true일 때 B의 Cholesky 인수분해의 치환 벡터. [R, ~, permB] = chol (B, "vector")로 얻음. 기본 1:n.

또한 A를 Af로 표시된 함수로 나타낼 수도 있어요. Af 뒤에는 Af가 받아들이는 벡터 인자의 길이를 정의하는 스칼라 인자 n이 와야 해요. Af는 함수 핸들, 인라인 함수, 또는 문자열일 수 있어요. 문자열일 때는 사용할 함수 이름을 담고 있습니다.

Af는 y = Af (x) 형태의 함수인데, Af의 필수 반환값은 sigma의 값에 따라 결정돼요. 네 가지 가능한 형태는:

A * x — sigma가 주어지지 않거나 "sm"이 아닌 문자열일 때.
A \ x — sigma가 0 또는 "sm"일 때.
(A - sigma * I) \ x — sigma가 0이 아닌 스칼라일 때. I는 A와 같은 크기의 단위행렬.
(A - sigma * B) \ x — 일반 고유값 문제에 대해.

반환 인자와 그 형태는 요청한 반환 인자의 수에 따라 달라져요. 반환 인자가 하나이면 찾은 k개의 고유값을 담은 길이 k의 열 벡터 d가 반환됩니다. 반환 인자가 두 개이면 V는 반환된 고유값에 해당하는 k개의 고유벡터를 열로 가진 n×k 행렬이에요. 고유값 자체는 D에 k×k 행렬로 반환되며, 대각 원소가 고유값입니다.

세 번째 반환 인자 flag는 수렴 상태를 돌려줘요. flag가 0이면 모든 고유값이 수렴한 거고, 다른 값은 수렴 실패를 나타냅니다.

Programming Notes: 작은 문제(n < 500)에는 eig (full (A))를 고려하세요.

ARPACK이 수렴에 실패하면 Lanczos 벡터 수(opt.p)를 늘리거나, 반복 횟수(opt.maxiter)를 늘리거나, 허용오차(opt.tol)를 줄여 보세요.

Reference: 이 함수는 R. Lehoucq, K. Maschhoff, D. Sorensen, C. Yang이 작성한 ARPACK 패키지에 기반합니다. 자세한 내용은 http://www.caam.rice.edu/software/ARPACK/ 참고.

svds

s = svds (A)
s = svds (A, k)
s = svds (A, k, sigma)
s = svds (A, k, sigma, opts)
[u, s, v] = svds (...)
[u, s, v, flag] = svds (...)

행렬 A의 몇 개의 특이값을 구해요. 특이값은 다음으로 계산합니다:

[m, n] = size (A);
s = eigs ([sparse(m, m), A;
                     A', sparse(n, n)])

eigs가 돌려준 고유값은 A의 특이값에 해당해요. 계산할 특이값의 개수는 k로 주어지며 기본값은 6이에요.

인자 sigma는 어떤 특이값을 찾을지 지정해요. sigma가 문자열 'L'(기본값)이면 A의 가장 큰 특이값을 찾아요. 그 외에는 sigma가 실수 스칼라여야 하고 sigma에 가장 가까운 특이값을 찾아요. 결과적으로 sigma = 0은 가장 작은 특이값을 찾아요. sigma가 상대적으로 작은 값이면 요청한 특이값 개수를 찾지 못할 가능성이 있으니, 그 경우 sigma를 늘려야 해요.

opts는 svds가 eigs에 넘겨줄 옵션을 정의하는 구조체예요. 이 구조체의 가능한 필드는 eigs에 문서화돼 있어요. 기본적으로 svds는 다음 세 필드를 설정해요:

tol — 특이값의 필요한 수렴 허용오차. 기본 1e-10. eigs에는 tol / sqrt (2)가 전달됨.
maxit — 최대 반복 횟수. 기본 300.
disp — 진단 출력 수준 (0|1|2). 0이면 진단 비활성화. 기본 0.

출력이 하나보다 많이 요청되면 svds는 A의 특이값 분해의 근사값을 돌려줍니다:

A_approx = u*s*v'

여기서 A_approx는 크기가 A와 같지만 계수(rank) k인 행렬이에요. flag는 알고리즘이 성공적으로 수렴하면 0, 그렇지 않으면 1을 돌려줍니다. 수렴 테스트는:

norm (A*v - u*s, 1) <= tol * norm (A, 1)

svds는 크고 희소한 행렬에서 몇 개의 특이값만 찾을 때 가장 좋아요. 그렇지 않으면 svd (full (A))가 더 효율적일 가능성이 높아요.

Footnotes

(10) CHOLMOD, UMFPACK, CXSPARSE 패키지는 Tim Davis가 작성했으며 http://faculty.cse.tamu.edu/davis/suitesparse.html에서 구할 수 있어요.

더 알아보기

octave-sparse-matrices — 희소 행렬 개념과 챕터 구성
octave-creating-sparse-matrices — 희소 행렬 만들기 함수들
octave-storage-of-sparse-matrices — 압축 열 저장 방식의 원리