고유값 문제

오류 제보
수치해석 계산물리 마지막 수정: 2026-10-02 04:22:11
○ 공개 검수 대기 AI 보조 작성 표시 기준 보기

1. 개요

고유값 문제
Eigenvalue Problem
표준형Ax = λx
일반화형Ax = λBx
조밀행렬QR 알고리즘(LAPACK)
희소행렬거듭제곱법, 란초스, 아놀디(ARPACK)
섭동 한계대칭 |Δλ| ≤ ‖ΔA‖ / 비대칭 κ(V)‖ΔA‖
등장 분야모드 해석, 좌굴, 밴드 구조

행렬을 통과시켜도 방향이 안 바뀌는 벡터가 있다. 걔가 그 행렬의 본질이다.

고유값 문제(eigenvalue problem)는 정방행렬 AA에 대해 Ax=λxA\mathbf{x} = \lambda\mathbf{x}를 만족하는 스칼라 λ\lambda(고유값, eigenvalue)와 영벡터가 아닌 벡터 x\mathbf{x}(고유벡터, eigenvector)를 찾는 문제다. 기하학적으로 고유벡터는 선형변환 AA를 가해도 방향이 바뀌지 않고 오직 λ\lambda배만큼 늘어나거나 줄어드는 특별한 방향이며, 고유값은 그 배율이다. 이 단순한 정의가 모드 해석, 좌굴, 밴드 구조 계산, 주성분 분석, 양자역학의 에너지 준위까지 온갖 수치해석 문제의 심장에 자리한다.

수학적으로는 특성방정식 det⁡(A−λI)=0\det(A - \lambda I) = 0의 근이 고유값이다. 하지만 이 다항식을 손으로 푸는 건 4×4만 넘어가도 지옥이고1, 수치적으로는 아예 하면 안 되는 짓이다(아래 § 검산에서 숫자로 확인한다). 실무에서는 전적으로 반복 알고리즘에 맡긴다. “고유값은 신에게 맡긴다”는 짤이 괜히 나온 게 아니다.

2. 표준 고유값 문제와 일반화 고유값 문제

가장 기본형은 표준 고유값 문제다.

Ax=λxA\mathbf{x} = \lambda \mathbf{x}

그런데 유한요소법 기반 구조해석에서는 두 개의 행렬이 등장하는 일반화 고유값 문제(generalized eigenvalue problem)를 만난다.

Ax=λBxA\mathbf{x} = \lambda B\mathbf{x}

진동 모드 해석에서는 AA가 강성행렬 KK, BB가 질량행렬 MM이 되어 Kϕ=ω2MϕK\boldsymbol{\phi} = \omega^2 M \boldsymbol{\phi} 형태로 나타난다. 여기서 고유값 λ=ω2\lambda = \omega^2은 고유진동수의 제곱, 고유벡터 ϕ\boldsymbol{\phi}는 진동 모드 형상이다. 만약 BB가 대칭 양정치라면 촐레스키 분해 B=LLTB = LL^\mathsf{T}로 L−1AL−Ty=λyL^{-1}AL^{-\mathsf{T}}\mathbf{y} = \lambda\mathbf{y} 꼴의 표준형으로 바꿀 수 있지만, 실무 솔버는 대개 일반화형을 직접 다루는 알고리즘(LAPACK의 dsygv, QZ 알고리즘)을 쓴다.

손으로 끝까지 따라갈 수 있는 최소 예제가 하나 있다. 질량 mm 두 개가 스프링 kk로 벽-질량-질량 순서로 묶인 2자유도 계는 K=k[2−1−11]K = k\begin{bmatrix}2 & -1\\ -1 & 1\end{bmatrix}, M=mIM = mI 이고, 고유값은

λ=3±52km={2.618034, 0.381966}km\lambda = \frac{3 \pm \sqrt{5}}{2}\frac{k}{m} = \{2.618034,\ 0.381966\}\frac{k}{m}

이다. 즉 ω={1.618034, 0.618034}k/m\omega = \{1.618034,\ 0.618034\}\sqrt{k/m} — 황금비 φ\varphi와 그 역수가 그대로 튀어나온다2. 이 예제는 아래 § 검산에서 알고리즘 수렴률을 재는 눈금자로 쓴다.

대칭 행렬의 고유값은 항상 실수이고 고유벡터는 서로 직교한다. 구조·양자 문제의 행렬이 대부분 대칭(에르미트)인 것은 축복인데, 덕분에 안정적이고 빠른 전용 알고리즘을 쓸 수 있다.

3. 알고리즘: QR 반복부터 란초스까지

고유값을 구하는 수치 방법은 크게 조밀행렬용과 희소행렬용으로 나뉜다.

**QR 알고리즘**은 조밀행렬의 모든 고유값을 구하는 표준 도구다. 행렬을 QR 분해하고 순서를 뒤집어 곱하는 Ak+1=RkQkA_{k+1} = R_k Q_k를 반복하면, 행렬이 점점 상삼각(대칭이면 대각) 형태로 수렴하며 대각 성분에 고유값이 드러난다. 실전에서는 먼저 하우스홀더 변환으로 헤센베르크(대칭이면 삼중대각) 형태로 줄인 뒤 시프트를 걸어 수렴을 가속한다 — 1961년 프랜시스(Francis)가 정리한 이 조합이 60여 년을 버틴 구조다. LAPACK이 이걸 극도로 최적화해 놨기 때문에, 수천 차원까지는 그냥 라이브러리를 부르면 끝난다3.

큰 희소행렬에서 고유값 몇 개만 필요할 때는 크리로프 부분공간법 계열이 답이다.

  • 거듭제곱법(power iteration): xk+1=Axk/∥Axk∥\mathbf{x}_{k+1} = A\mathbf{x}_k / \|A\mathbf{x}_k\|를 반복하면 절댓값이 가장 큰 고유값의 고유벡터로 수렴한다. 구글의 초기 페이지랭크가 바로 이 원리였다.
  • 란초스 알고리즘: 대칭(에르미트) 행렬 전용. 크리로프 부분공간을 만들며 원 행렬을 작은 삼중대각 행렬로 사영해, 극단 고유값 몇 개를 빠르게 뽑는다. 밴드 구조 계산과 대규모 진동 해석의 주력.
  • 아놀디 알고리즘: 비대칭 행렬용 일반화. 삼중대각이 아니라 헤센베르크 행렬로 사영되므로 3항 점화식이 깨지고 직교화 비용을 다 물어야 한다. ARPACK의 심장이 이것(정확히는 암시적 재시작 아놀디, IRAM)이며, 유체 유동의 선형 안정성 해석 등에 쓰인다.
  • 역반복법과 shift-invert: (A−σI)−1(A - \sigma I)^{-1}에 거듭제곱법을 걸면 σ\sigma에 가장 가까운 고유값이 가장 큰 고유값으로 바뀌어 먼저 수렴한다. 스펙트럼 내부의 특정 대역만 뽑는 표준 트릭이고, 대가는 매 반복 선형계 한 번(LULU 또는 전처리 반복)이다.

행렬 전체를 다루지 않고 “행렬-벡터 곱 AvA\mathbf{v}만 계산할 수 있으면 된다”는 점이 크리로프 방법의 강력함이다. 수백만 차원 희소행렬에서 저주파 모드 몇 개만 뽑는 문제가 이 방식으로 실용화된다.

4. 어디서 튀어나오는가: 모드·좌굴·밴드구조

고유값 문제가 특히 자주 출몰하는 세 현장이 있다.

진동 모드 해석: 구조물이 스스로 진동하려는 고유진동수와 그 모양(모드 형상)을 Kϕ=ω2MϕK\boldsymbol{\phi} = \omega^2 M\boldsymbol{\phi}의 고유값·고유벡터로 얻는다. 외력 주파수가 고유진동수에 근접하면 공진(resonance)이 일어나 구조물이 박살 나므로, 설계 단계에서 반드시 확인한다. 타코마 다리 붕괴가 교과서적 반면교사다4.

좌굴 해석: 압축 하중을 받는 기둥·판이 어느 순간 갑자기 옆으로 휘어 무너지는 좌굴 임계하중은 (K+λKG)ϕ=0(K + \lambda K_G)\boldsymbol{\phi} = 0 형태의 일반화 고유값 문제로 구한다. 여기서 KGK_G는 기하강성행렬, 최소 고유값 λ\lambda가 좌굴 안전계수, 고유벡터가 좌굴 모드 형상이다.

밴드 구조 계산: 결정 속 전자의 에너지 준위는 주기적 퍼텐셜에 대한 슈뢰딩거 방정식을 이산화한 일반화 고유값 문제 Hψ=ESψH\psi = E S\psi로 얻는다. 각 파수벡터 k\mathbf{k}마다 고유값을 풀어 에너지-운동량 관계를 그린 것이 밴드 구조이며, 반도체가 도체인지 절연체인지가 여기서 판가름 난다. 밀도범함수이론 계산의 알맹이가 바로 이 고유값 문제 반복이다.

5. 검산

말로 “불안정하다”고 하면 안 믿는 사람이 많으니 숫자를 박아 둔다. 아래는 모두 배정밀도로 직접 돌린 결과다.

특성방정식을 풀면 안 되는 이유. 윌킨슨(Wilkinson)의 고전적 반례5는 W(x)=∏i=120(x−i)W(x) = \prod_{i=1}^{20}(x-i) 다. 고유값(=근)은 정수 1부터 20까지로 완벽하게 떨어져 있다. 여기서 x19x^{19} 계수 −210-210에 2−23≈1.19×10−72^{-23} \approx 1.19\times10^{-7} 만 더한다 — 상대 섭동 5.7×10−105.7\times10^{-10}, 배정밀도 반올림보다 겨우 몇 자리 큰 양이다. 결과:

원래 근섭동 후
1 ~ 6변화 없음 (소수 4자리까지 동일)
8, 98.0072, 8.9173
10, 1110.0953 ± 0.6436i (복소쌍)
16, 1716.7307 ± 2.8126i
2020.8469

정수 근 열 개가 복소 평면으로 날아갔다. 반면 같은 고유값을 갖는 대각행렬 diag(1,…,20)\mathrm{diag}(1,\dots,20)에 같은 크기의 섭동을 주면 고유값은 10−710^{-7} 수준으로만 움직인다. 즉 나쁜 것은 고유값이 아니라 다항식 계수라는 표현이다.

대칭이 얌전한 이유. 대칭(에르미트)이면 바일(Weyl)의 부등식으로 ∣λi(A+ΔA)−λi(A)∣≤∥ΔA∥2|\lambda_i(A+\Delta A) - \lambda_i(A)| \le \|\Delta A\|_2 가 보장된다. 증폭 계수가 1이다. 일반(대각화 가능) 행렬에서는 바우어-파이크 정리가 min⁡i∣μ−λi∣≤κ(V) ∥ΔA∥2\min_i|\mu - \lambda_i| \le \kappa(V)\,\|\Delta A\|_2 를 줄 뿐이고, 여기서 κ(V)\kappa(V)는 고유벡터 행렬의 조건수다. [1101+δ]\begin{bmatrix}1 & 1\\ 0 & 1+\delta\end{bmatrix} 에서 δ=10−6\delta = 10^{-6} 이면 κ(V)=2×106\kappa(V) = 2\times10^{6} — 고유값이 가까워질수록 고유벡터가 평행해지고 그만큼 증폭된다.

결함행렬의 ε1/N\varepsilon^{1/N} 법칙. 극단은 고유벡터가 아예 부족한 경우다. N×NN\times N 조르당 표준형 블록의 모서리에 ε\varepsilon 을 넣으면 고유값은 ε1/N\varepsilon^{1/N} 만큼 흩어진다. N=10N=10, ε=10−16\varepsilon = 10^{-16}(배정밀도 반올림!)이면 ∣Δλ∣=2.5×10−2|\Delta\lambda| = 2.5\times10^{-2} — 입력 오차를 2.5×10142.5\times10^{14}배 증폭한 것이다. 비정규 행렬에서 고유값 하나하나를 신뢰하지 말라는 경고가 여기서 나온다.

수렴률은 고유값 간격이 정한다. 앞의 2자유도 예제 K/kK/k 에 거듭제곱법을 걸고 레일리 몫 ρ=xTAx/xTx\rho = \mathbf{x}^\mathsf{T}A\mathbf{x}/\mathbf{x}^\mathsf{T}\mathbf{x} 으로 고유값을 추정하면, 대칭이므로 오차가 비 (λ2/λ1)2=0.1458982=0.021286(\lambda_2/\lambda_1)^2 = 0.145898^2 = 0.021286 로 줄어야 한다. 실측: 1스텝 1.80×10−21.80\times10^{-2} → 2스텝 3.87×10−43.87\times10^{-4}(비 0.0215) → 5스텝 3.7×10−93.7\times10^{-9} → 10스텝에 배정밀도 한계. 이론과 실측이 세 자리까지 맞는다. 고유값이 뭉쳐 있으면(λ2/λ1→1\lambda_2/\lambda_1 \to 1) 이 비가 1에 붙어 수렴이 기어간다는 뜻이기도 하다 — shift-invert가 필요한 진짜 이유다.

6. 수치적 주의사항

  • 특성방정식을 직접 풀지 마라. 위 윌킨슨 예제가 답이다. 반드시 QR·란초스 같은 반복법을 쓴다. 역설적으로 QR 알고리즘은 컴패니언 행렬을 거쳐 다항식 근 찾기에 쓰이는 쪽이 더 안전하다.
  • 결함행렬(defective matrix). 고유값이 중복되면 고유벡터가 부족해 대각화가 불가능하고, 이 경우 조르당 표준형을 써야 한다. 다만 조르당 형태는 수치적으로 계산 불가능에 가까워, 실무에서는 슈어 분해로 우회한다.
  • 고유벡터는 고유값보다 민감하다. 고유값 자체는 후진 오차 해석 관점에서 안정적으로 나와도, 두 고유값이 가까우면 그 둘이 펼치는 부분공간만 신뢰할 수 있고 개별 고유벡터는 둘 사이에서 자유롭게 섞인다. 모드 형상 하나를 들고 와서 해석하기 전에 고유값 간격을 먼저 본다.
  • 수렴 판정. 크리로프 방법은 요구하는 고유값 개수와 전처리기, 재시작 전략에 따라 성능이 천차만별이다. 잔차 ∥Ax−λx∥\|A\mathbf{x} - \lambda\mathbf{x}\| 는 싸게 구할 수 있으니 반드시 확인하고, 내부 고유값이 필요하면 ARPACK의 shift-invert 모드를 쓴다.

7. 관련 문서

8. Footnotes

  1. 5차 이상 다항식은 근의 공식이 존재하지 않는다(아벨-루피니 정리). 즉 5×5 행렬부터는 고유값을 유한 번의 사칙연산과 거듭제곱근으로 표현하는 것 자체가 원리적으로 불가능하다. 반복법이 선택이 아니라 필연인 이유. ↩

  2. 우연이 아니다. 특성방정식이 λ2−3λ+1=0\lambda^2 - 3\lambda + 1 = 0 이고, φ2=φ+1\varphi^2 = \varphi + 1 을 쓰면 λ=φ2, φ−2\lambda = \varphi^2,\ \varphi^{-2} 임이 바로 나온다. 2자유도 숙제에서 1.6181.618 이 보이면 계산을 맞게 한 것이니 안심해도 된다. ↩

  3. numpy.linalg.eig 한 줄이 사실은 수십 년간 다져진 LAPACK의 dgeev/dsyev를 부르는 것이다. 이 한 줄 뒤에 하우스홀더 변환과 시프트 QR과 밸런싱이 숨어 있다. 편함에 감사하며 살자. ↩

  4. 흔히 공진으로만 설명되지만 실제 타코마 붕괴는 플러터(flutter)라는 유체-구조 연성 불안정에 가깝다. 그래도 “고유진동수 근처는 위험하다”는 교훈 자체는 유효하다. 학부 강의의 단골 낚시. ↩

  5. 윌킨슨 본인이 이 다항식을 “perfidious polynomial”(배신하는 다항식)이라 불렀다. 근이 깔끔한 정수인데도 계수 표현이 배신한다는 뜻이다. 교훈: 수학적으로 동등한 두 표현이 수치적으로 동등하다는 보장은 전혀 없다. ↩

검증 자료

이 문서의 주요 출처

  1. Francis (1961), The QR Transformation: A Unitary Analogue to the LR Transformation - Part 1, Comput. J. 4(3), 265-271
  2. Lanczos (1950), An Iteration Method for the Solution of the Eigenvalue Problem of Linear Differential and Integral Operators, J. Res. Natl. Bur. Stand. 45, 255-282
  3. Arnoldi (1951), The Principle of Minimized Iterations in the Solution of the Matrix Eigenvalue Problem, Quart. Appl. Math. 9(1), 17-29
  4. Bauer & Fike (1960), Norms and Exclusion Theorems, Numer. Math. 2, 137-141
  5. Anderson et al., LAPACK Users' Guide, 3rd ed. (SIAM, 1999) - 2.3 대칭/비대칭 고유문제 드라이버 루틴
  6. Lehoucq, Sorensen & Yang, ARPACK Users' Guide (SIAM, 1998) - 암시적 재시작 아놀디와 shift-invert 모드
이 문서가 문제를 푸는 데 도움이 됐나요? 한 번의 클릭이 다음 검수 순서를 정합니다.

의견을 반영할게요. 감사합니다.

다음으로 읽기 희소행렬 약 5분 더 탐색 →