스텀 수열

1. 개요[편집]

스텀 수열 Sturm sequence
원형	실계수 다항식의 실근 개수 세기 (Sturm, 1829)
세는 법	$V(a) - V(b)$ — 양 끝에서의 부호변화 수 차이
행렬판	대칭 삼중대각의 선행 주소행렬식 점화식
실무 구현	$T - \lambda I = LDL^{\top}$ 의 음수 피벗 개수 (관성)
비용	계수 함수 1회 $O(n)$, 나눗셈 $n$ 번
LAPACK	dstebz (커널 dlaebz) + dstein

근을 구하지 말고 세기만 해라. 세는 건 나누기 $n$ 번이면 되고, 세는 것만 할 줄 알아도 원하는 근을 원하는 정밀도로 잡아낼 수 있다.

스텀 수열(Sturm sequence)은 어떤 점 $x$ 에서의 부호변화 개수가 그 점보다 작은 근의 개수를 세어 주도록 구성한 함수열이다. 1829년 샤를 스텀이 실계수 다항식의 구간 내 실근 개수를 정확히 세는 정리로 도입했고, 수치선형대수에서는 이것이 대칭 삼중대각행렬의 고유값을 원하는 것만 골라 뽑는 알고리즘의 근거가 된다.

핵심은 다음 한 줄이다. 스칼라 함수

$$\nu(\lambda) \;=\; \#\{\,i : \lambda_i(T) < \lambda \,\}$$

를 $O(n)$ 에 계산할 수 있으면, 이분법만으로 임의의 고유값을 임의의 정밀도로 격리할 수 있다. 고유값 문제를 근 찾기가 아니라 개수 세기로 환원하는 것이 이 문서의 전부이며, “5000차 행렬에서 300번째부터 320번째 고유값만 필요하다”는 요구에 답할 수 있는 유일한 고전적 경로다.

2. 다항식에서의 스텀 정리[편집]

제곱인수가 없는 실계수 다항식 $p$ 에 대해 유클리드 호제법을 나머지의 부호를 뒤집어 돌린다.

$$p_0 = p,\quad p_1 = p',\quad p_{k+1} = -\operatorname{rem}(p_{k-1},\,p_k)$$

상수가 나올 때까지 반복한 열 $p_0,p_1,\dots,p_s$ 가 스텀 연쇄(Sturm chain)다. $V(x)$ 를 이 열의 값에서 0을 빼고 센 부호변화 횟수라 하면, $p(a)\ne 0$, $p(b)\ne 0$ 인 $a<b$ 에 대해

$$\#\{\,\text{$(a,b]$ 안의 서로 다른 실근}\,\} \;=\; V(a) - V(b)$$

가 성립한다. 제곱인수가 있어도 호제법이 $\gcd(p,p')$ 에서 멈출 뿐 공식은 그대로 살아 있고, 이때도 세어지는 것은 중복도를 무시한 서로 다른 근의 개수다. 근사도 확률도 아닌 정확한 정수를 준다는 점이 데카르트 부호 규칙 같은 부등식 도구와 결정적으로 다른 지점이다.¹

작동 원리는 국소적으로 확인된다. $x$ 가 증가하며 중간 항 $p_k$ 의 근을 지날 때는 이웃 $p_{k-1}, p_{k+1}$ 의 부호가 서로 반대라 부호변화 총수가 변하지 않고, $p_0$ 의 근을 지날 때만 정확히 하나씩 줄어든다. 그래서 $V$ 는 계단처럼 단조감소하고, 그 낙차가 곧 근의 개수다.

3. 삼중대각행렬에서는 공짜로 나온다[편집]

이제 본론이다. 대칭 삼중대각행렬

$$T = \begin{bmatrix} \alpha_1 & \beta_1 & & \\ \beta_1 & \alpha_2 & \ddots & \\ & \ddots & \ddots & \beta_{n-1} \\ & & \beta_{n-1} & \alpha_n \end{bmatrix}$$

의 앞쪽 $k \times k$ 선행 주부분행렬을 $T_k$ 라 하고 $p_k(\lambda) = \det(T_k - \lambda I)$ 로 두자. 마지막 행으로 행렬식을 전개하면 3항 점화식이 그대로 튀어나온다.

$$p_k(\lambda) = (\alpha_k - \lambda)\,p_{k-1}(\lambda) - \beta_{k-1}^{2}\,p_{k-2}(\lambda),\qquad p_0 = 1,\ p_{-1} = 0$$

모든 $\beta_i \ne 0$ 인 비축약(unreduced) $T$ 에서는 $T_k$ 의 고유값이 $T_{k+1}$ 의 고유값을 강하게 교대(strict interlacing)하고, 그 결과 $p_0,p_1,\dots,p_n$ 이 스텀 수열의 조건을 만족한다. 결론은 딱 떨어진다.

$$\text{$\{p_0(\lambda),\dots,p_n(\lambda)\}$ 의 부호변화 수} \;=\; \nu(\lambda) = \#\{\lambda_i(T) < \lambda\}$$

호제법도, 다항식 나눗셈도 없다. 삼중대각 구조가 스텀 연쇄를 이미 들고 있다. 극단값을 넣어 보면 바로 납득된다. $\lambda \to -\infty$ 면 모든 $p_k > 0$ 이라 부호변화 0개이고 실제로 그보다 작은 고유값은 0개, $\lambda \to +\infty$ 면 $p_k \sim (-\lambda)^k$ 라 부호가 매 항 뒤집혀 $n$ 개다.

4. 행렬식 대신 관성 — 실무 구현[편집]

점화식을 그대로 쓰면 $n$ 이 조금만 커져도 $p_k$ 가 오버플로/언더플로로 죽는다. $\lambda^k$ 규모로 자라기 때문이다. 그래서 실제 코드는 비를 굴린다.

$$q_k = \frac{p_k}{p_{k-1}}\ :\qquad q_1 = \alpha_1 - \lambda,\qquad q_k = (\alpha_k - \lambda) - \frac{\beta_{k-1}^{2}}{q_{k-1}}$$

그러면 $p_{k-1}$ 과 $p_k$ 의 부호가 다른 것과 $q_k < 0$ 인 것이 같으므로

$$\nu(\lambda) = \#\{\,k : q_k < 0\,\}$$

이다. 그런데 이 $q_k$ 는 정확히 $T - \lambda I = LDL^{\top}$ 의 피벗 $d_k$ 다. 즉 음의 피벗 개수를 세는 것이고, 이건 실베스터 관성 법칙이 보장하는 불변량 — 합동변환은 관성(양·음·영 고유값의 개수)을 보존한다는 그 정리 — 을 세는 것과 같다. 부호변화라는 다항식 언어가 관성이라는 선형대수 언어로 번역되는 지점이고, 이 형태가 오버플로도 없고 $n$ 번의 나눗셈으로 끝난다.

남은 문제는 $q_{k-1} = 0$ 일 때의 0으로 나누기다. 처방이 둘 있다.

• 미세 섭동. $q_{k-1}$ 이 0이거나 지나치게 작으면 적절한 부호의 아주 작은 값(보통 $\beta_{k-1}^2 \cdot u$ 규모)으로 갈아치운다. LAPACK dlaebz 가 이 방식이다.
• 예외 산술 그대로 두기. IEEE 754 에서 $q_{k-1}=0$ 이면 $q_k = -\infty$ 가 되고, 다음 항은 $\beta^2/(-\infty) = -0$ 이라 $q_{k+1} = \alpha_{k+1}-\lambda$ 로 자동 복구된다. 부호까지 맞게 나온다. 데믈–리(1994)가 “예외 처리를 끄지 말고 이용하라”는 취지로 정리한 이 트릭은 분기 없는 내부 루프를 만들어 실측을 크게 앞당긴다.²

여기에 하나 더 필요한 것이 부동소수점에서의 단조성이다. 계산된 계수 함수 $\hat\nu$ 가 $\lambda$ 에 대해 단조가 아니면 이분법이 모순된 구간을 만들어 무한루프에 빠질 수 있다. 데믈–딜론–렌(1995)이 표준적인 부동소수점 연산 가정 아래 $\hat\nu$ 의 단조성을 증명해 두었고, 그래서 이분법 기반 고유값 코드는 “대체로 잘 돈다”가 아니라 증명된 정확성을 갖는 드문 사례가 됐다.

5. 이분법 — 원하는 고유값만[편집]

계수 함수가 준비되면 나머지는 초등학교 수준이다.

1. 게르슈고린 원판으로 전체 스펙트럼을 담는 초기 구간 $[a,b]$ 를 잡는다.
2. $k$ 번째 고유값을 원하면 $\nu(a) < k \le \nu(b)$ 인 구간을 유지하며 이분법으로 반씩 줄인다.
3. 구간 폭이 목표 정밀도 이하가 되면 멈춘다. 배정도 기준 대략 50~60회.

비용은 고유값 하나당 $O(n\log\frac{b-a}{\varepsilon})$ 이고, 고유값 $k$ 개를 원하면 $O(kn)$ 규모다. 전체를 다 구하면 $O(n^2)$ 로 QR 알고리즘과 같지만 상수가 커서 느리다. 이 방법의 존재 이유는 속도가 아니라 선택성이다.

• 구간 지정(RANGE='V'): ”$[0, \omega_{\max}^2]$ 안의 모든 모드” — 모드 해석에서 매일 쓰는 질의다.
• 번호 지정(RANGE='I'): ”$k$ 번째부터 $\ell$ 번째까지”.
• 완벽한 병렬성: 고유값마다 계산이 완전히 독립이고 통신이 0이다. 분산 메모리에서 이보다 편한 구조가 없다.

고유벡터는 별도다. 격리된 $\hat\lambda_i$ 를 시프트로 삼아 역반복법(dstein)을 돌리는데, 여기서 약점이 드러난다. 고유값이 뭉쳐 있으면 역반복이 같은 방향으로 수렴해 나온 벡터들이 직교하지 않고, 그람-슈미트 재직교화를 붙여야 한다. 무리 크기가 커지면 그 비용이 $O(n^3)$ 으로 튀고 직교성 보장도 없다. 이 결함을 정면으로 없앤 것이 MRRR(dstemr)이며, 재미있게도 MRRR 역시 시프트된 $LDL^{\top}$ 표현 위에서 돌아간다 — 관성이라는 도구는 그대로 두고 고유벡터 쪽만 갈아 끼운 셈이다.

6. 어디서 또 만나는가[편집]

• 일반화 문제와 스텀 검사. $K - \lambda M = LDL^{\top}$ 의 음의 피벗 개수는 $\lambda$ 보다 작은 고유진동수 제곱의 개수를 그대로 준다(촐레스키 분해 없이도 된다). 구조해석 상용 코드가 모드 추출 후 “스텀 시퀀스 체크”를 돌리는 이유가 이것이다. 모드를 빠뜨렸는지 여부를 독립적으로 검증할 수 있는 유일한 값싼 수단이라, 반복법으로 뽑은 란초스 알고리즘 결과의 누락(missed mode) 진단에 표준으로 쓰인다.
• 관성 기반 스펙트럼 분할. 구간을 관성으로 재귀 분할해 각 조각을 다른 프로세스에 던지는 것이 대규모 병렬 고유값 코드(예: 시프트-역변환 기반 솔버)의 부하 분산 전략이다.
• 다항식 실근 격리. 원래 용도인 이쪽에서는 정작 스텀이 밀렸다. 정확 산술로 스텀 연쇄를 만들면 계수가 기하급수적으로 불어나 큰 차수에서 실용성이 떨어져, 현대 컴퓨터 대수 시스템은 데카르트 부호 규칙 기반 분할(Descartes/VAS)을 기본으로 쓴다. 개수를 정확히 세는 정리가, 개수를 세는 실무에서는 지고 개수를 세지 않는 실무(고유값)에서는 이겼다는 얄궂은 결말.³

7. 관련 문서[편집]

• 삼중대각화 · 고유값 문제 · 분할 정복 고유값 알고리즘 · 역반복법
• 실베스터 관성 법칙 · 행렬식 · 최소다항식 · 이분법
• LU 분해 · 촐레스키 분해 · 그람-슈미트
• 부동소수점 연산 · 후진 오차 해석 · LAPACK · 병렬 컴퓨팅
• 모드 해석 · 란초스 알고리즘 · 뉴턴-랩슨법

8. Footnotes[편집]

1. 스텀이 이 정리를 내놓았을 때 푸리에가 극찬했다는 일화가 전해진다. 당시는 대수방정식의 근을 “구하는” 것이 수학의 중심 과제였는데, 근을 하나도 구하지 않고 개수만 정확히 세는 방법이 등장했으니 발상 자체가 신선했던 것. 참고로 스텀은 광속 측정과 리우빌과 함께 만든 스텀-리우빌 이론으로도 유명한데, 셋 다 완전히 다른 분야다. ↩

2. Demmel & Li (1994), “Faster numerical algorithms via exception handling”. 보통은 $\infty$ 와 NaN 을 보면 코드가 잘못됐다고 배우는데, 여기서는 반대로 분기문을 넣는 쪽이 느리고 위험하다. 물론 이 트릭은 컴파일러가 -ffast-math 로 IEEE 준수를 꺼 버리면 조용히 망가진다. 성능 플래그가 정확성을 먹는 고전적 사례. ↩

3. 계수 폭발이 어느 정도냐면, 정수 계수 다항식의 스텀 연쇄를 유리수 산술로 만들면 항의 계수 자릿수가 차수에 대해 이차적으로 늘어난다. 부분종결식(subresultant) PRS 로 억누르는 기법이 있긴 하지만, 그래도 20차만 넘어가면 “이럴 거면 그냥 근을 구하자”는 소리가 나온다. ↩