손으로 푸는 선형대수 15장 — 양정치와 이차형식: 공분산이 왜 대칭 양반정치인가

15 손으로 푸는 선형대수 · 제4부 고유값과 특이값

양정치와 이차형식 — 공분산이 왜 대칭 양반정치인가

행렬에 벡터를 양쪽에서 곱하면 숫자 하나가 나온다. 그 숫자가 항상 양수면 함수의 모양이 그릇이 되고, 그릇이면 최솟값이 하나뿐이다. 최적화가 잘 되는 조건, 공분산이 얌전한 이유, 그리고 17장 PCA가 성립하는 근거가 전부 이 한 성질에서 나온다.

1. 직관 — 그릇인가 안장인가 도랑인가
2. 수식 읽는 법 — 세 가지 동치 판정과 A^TA
3. 손으로 풀기 — G에 세 판정을 전부 적용
4. 코드 — 판정 세 개를 직접 구현하고 깨뜨려 본다
5. 시각화 — 이차형식의 등고면은 타원체다

약어 및 기호 정의

이차형식quadratic form q(x) = x^TSx

벡터를 넣으면 스칼라가 나오는 함수. 모든 항이 x_ix_j꼴, 즉 차수가 정확히 2인 다항식이다.

양정치positive definite S ≻ 0

대칭행렬 S가 모든 x ≠ 0에 대해 x^TSx > 0일 때. 그릇 모양.

양반정치positive semidefinite S ⪰ 0

모든 x에 대해 x^TSx ≥ 0. 어떤 x ≠ 0에서 0이 될 수 있다. 바닥이 평평한 도랑.

부정치indefinite

x^TSx가 양수도 되고 음수도 되는 경우. 안장 모양이고 최솟값이 없다.

선행주소행렬식 leading principal minor d_k

왼쪽 위 k×k 부분행렬의 행렬식. 7장에서 G의 값이 2, 27, 144임을 이미 구했다.

실베스터 판정법Sylvester's criterion

대칭행렬이 양정치일 필요충분조건은 모든 선행주소행렬식이 양수라는 정리.

콜레스키 분해Cholesky decomposition S = LL^T

L은 대각이 양수인 하삼각행렬. 양정치일 때만 존재하며 유일하다. "행렬의 제곱근"에 해당한다.

LDL^T 분해

L이 대각 1인 하삼각, D가 대각행렬. 5장 LU 소거의 대칭 버전이며 D가 곧 피벗이다.

중심화 행렬centering matrix H = I − 1n11^T

왼쪽에 곱하면 각 열에서 그 열의 평균을 빼 준다. H^T = H, H² = H — 10장의 투영행렬이다.

표본공분산 sample covariance Σ̂

Σ̂ = 1n−1 X_c^TX_c. 중심화한 자료행렬의 그람 행렬을 n−1로 나눈 것.

헤시안Hessian ∇²f

2계 편미분을 모은 대칭행렬. 이것이 양정치인 임계점이 곧 국소 최솟값이다. 19장에서 다시 만난다.

조건수 condition number κ₂

8장에서 본 σ_max/σ_min. 우리 A는 4, G = A^TA는 16.

STEP 01

직관 — 그릇인가 안장인가 도랑인가

지금까지 행렬은 늘 한쪽에서만 벡터를 만났다. Ax는 벡터를 벡터로 보내고(1장), Px는 벡터를 평면 위로 내린다(10장). 이번 장의 연산은 다르다. 같은 벡터를 양쪽에서 곱한다.

       행벡터        행렬        열벡터        스칼라
       (1×n)        (n×n)        (n×1)         (1×1)
        xᵀ     ·      S     ·      x      =    숫자 하나

모양을 따라가면 결과가 1×1, 즉 숫자 하나다.

즉 q(x) = x^TSx는 벡터를 넣으면 숫자가 나오는 함수다. 이런 함수를 이차형식quadratic form이라 부른다. 왜 "이차"인가. 성분으로 풀어 쓰면 모든 항이 x_ix_j 꼴이라 차수가 정확히 2이기 때문이다. 일차항도 상수항도 없다. n=1이면 q(x) = sx², 그냥 포물선이다.

포물선 sx²에서 s의 부호가 모든 것을 결정했던 것을 떠올리자. s > 0이면 위로 열린 포물선이라 최솟값이 x=0 하나뿐이고, s < 0이면 뒤집혀서 최솟값이 없다. n차원에서도 똑같은 분류가 있다. 다만 방향이 여러 개라서 "어떤 방향에서는 위로, 어떤 방향에서는 아래로"라는 제3의 경우가 생긴다.

도해 1. 이차형식 q(x) = x^TSx의 표면 세 가지. 고유값의 부호만으로 결정된다. 최적화하는 입장에서 관심은 왼쪽 하나다 — 그릇이면 아무 데서 굴려도 같은 바닥에 도착한다. 오른쪽 도랑은 바닥이 평평한 골짜기라, 최솟값은 존재하지만 어디인지가 유일하지 않다. 3장의 랭크결손, 12장의 다중공선성이 전부 이 그림이다.

이 분류가 왜 중요한가. 우리가 ML에서 푸는 문제는 거의 전부 "무엇인가를 최소화"하는 문제고, 손실함수를 최소점 근처에서 2차까지 근사하면 반드시 이차형식이 나온다. 테일러 전개를 보자.

f(x₀ + h) ≈ f(x₀) + ∇f(x₀)^Th + 12 h^T(∇²f) h 임계점에서는 ∇f = 0이므로, 모양을 결정하는 것은 마지막 이차형식 항뿐이다. (15.1)

임계점에서 헤시안 ∇²f가 양정치면 그 점은 최솟값이고, 부정치면 안장점이다. 신경망 학습이 안장점에서 헤매는 이야기, 뉴턴법이 헤시안의 역행렬을 쓰는 이야기가 전부 여기서 갈라진다. 그리고 최소제곱에서 헤시안은 정확히 2A^TA다(19장에서 유도한다). 그러니 우리의 G = A^TA가 양정치인지가 곧 "10장의 최소제곱 해 β̂ = (1,1,1)이 유일한 최솟값인가"라는 질문이다.

STEP 02

수식 읽는 법 — 세 가지 동치 판정

(가) 이차형식을 성분으로 펼치기

정의는 단순하다. S가 n×n 대칭행렬일 때

q(x) = x^TSx = nΣi=1 nΣj=1 s_ij x_i x_j 대각항은 s_iix_i², 비대각항은 대칭성 때문에 2s_ijx_ix_j로 두 번씩 들어온다. (15.2)

여기서 왜 대칭행렬만 다루는가에 답하고 넘어가자. S가 대칭이 아니어도 x^TSx는 계산된다. 그런데 결과는 스칼라이므로 자기 전치와 같다. x^TSx = (x^TSx)^T = x^TS^Tx. 둘을 더해 반으로 나누면

x^TSx = x^T S + S^T2 x 비대칭 부분은 이차형식에 아무 기여도 하지 않는다. 그래서 처음부터 대칭만 본다.

(나) 세 가지 판정 — 전부 같은 말이다

대칭행렬 S에 대해 다음 세 명제는 서로 동치다.

판정	내용	어디서 왔나	비용
① 고유값	모든 고유값 λi > 0	14장 스펙트럼 분해	O(n³), 상수 큼
② 실베스터	모든 선행주소행렬식 dk > 0	7장 행렬식	손계산에 좋음
③ 콜레스키	S = LLT가 존재 (L 대각 > 0)	5장 LU 소거	n³/3, 실무 표준

①이 왜 참인가. 14장에서 대칭행렬은 S = QΛQ^T로 쓸 수 있고 Q는 직교행렬이었다. y = Q^Tx로 두면

x^TSx = x^TQΛQ^Tx = y^TΛy = λ₁y₁² + λ₂y₂² + ⋯ + λ_ny_n² 직교 좌표계를 고유벡터 방향으로 돌려 놓으면 교차항이 전부 사라진다. 남는 것은 포물선 n개의 합이다. (15.3)

Q가 직교라 x ≠ 0 ⟺ y ≠ 0이다. 그러니 q > 0이 항상 성립할 필요충분조건은 모든 계수 λ_i가 양수인 것이다. 도해 1의 세 그림은 결국 (λ₁, λ₂)가 (+,+)인가 (+,−)인가 (+,0)인가의 차이다.

③이 왜 실무 표준인가. 고유값을 다 구하는 것은 비싸다. 콜레스키는 소거와 같은 비용으로 돌면서, 실패 자체가 판정이 된다 — 진행 도중 제곱근 안이 양수가 아니면 그 순간 양정치가 아니다. numpy와 scipy가 LinAlgError를 던지는 지점이 정확히 여기다.

(다) A^TA는 항상 양반정치다 — 한 줄 증명

이 장의 핵심 한 줄이다. 임의의 m×n 행렬 A와 임의의 x에 대해

x^T(A^TA)x = (Ax)^T(Ax) = ‖Ax‖² ≥ 0 괄호를 옮겨 묶었을 뿐이다. 그런데 오른쪽은 길이의 제곱이라 음수가 될 수 없다. (15.4)

증명 끝이다. 가정도 조건도 없다. A가 어떻게 생겼든 A^TA는 양반정치다. 그리고 등호가 언제 성립하는지를 보면 양정치 여부가 갈린다.

x ᵀ(AᵀA)x = 0   ⟺   ‖Ax‖ = 0   ⟺   Ax = 0   ⟺   x ∈ N(A)

6장의 영공간이 여기서 돌아온다.

따라서 A의 열이 일차독립이면 N(A) = {0}이므로 등호는 x = 0에서만 성립하고, A^TA는 양정치다. 열이 종속이면 0이 되는 x ≠ 0이 존재하므로 양반정치이지 양정치가 아니다.

이 한 줄이 이 시리즈에서 G가 늘 얌전했던 이유를 전부 설명한다. 7장에서 det G = 144 > 0이었던 것, 선행주소행렬식이 2, 27, 144로 전부 양수였던 것, 14장에서 고유값이 16, 9, 1로 전부 양수였던 것, 10장에서 (A^TA)⁻¹이 존재해 투영행렬을 만들 수 있었던 것 — 전부 식 (15.4)의 따름정리다. 우연이 하나도 없다.

도해 2. x^TGx를 두 가지 경로로 읽는다. 아래 경로는 "행렬 하나를 양쪽에서 곱한다"이고, 위 경로는 "A로 보낸 뒤 길이를 잰다"이다. 같은 숫자이므로 부호가 구조적으로 결정된다. 등호가 성립하려면 Ax 자체가 0이어야 하고, 그것은 x가 영공간에 있을 때뿐이다.

STEP 03

손으로 풀기 — G에 세 판정을 전부 적용한다

늘 쓰던 A와 G = A^TA다.

      ⎡ 1  1  0 ⎤        ⎡  2   1   3 ⎤
 A =  ⎢ 0  2  0 ⎥   G =  ⎢  1  14   3 ⎥      λ(G) = 16, 9, 1      det G = 144
      ⎢ 1  0  3 ⎥        ⎣  3   3  10 ⎦
      ⎣ 0  3  1 ⎦

(15.5)

(가) x^TGx를 세 방향에서 — 그리고 ‖Ax‖²와 맞춰 본다

먼저 식 (15.2)를 G에 대해 펼쳐 두자. 대각에서 2, 14, 10, 비대각에서 2·1, 2·3, 2·3이 나온다.

q(x) = 2x₁² + 14x₂² + 10x₃² + 2x₁x₂ + 6x₁x₃ + 6x₂x₃ (15.6)

방향 1 — x = (1, 0, 0). 식 (15.6)에 넣으면 첫 항만 살아남아 q = 2. 행렬로 보면 G의 (1,1) 성분을 그냥 집어 온 것이다. 일반적으로 e_i^TSe_i = s_ii이므로 양정치 행렬은 대각이 전부 양수여야 한다(필요조건이지 충분조건은 아니다).

검산: A(1,0,0) = c₁ = (1, 0, 1, 0), ‖c₁‖² = 1 + 0 + 1 + 0 = 2. 일치한다. 사실 G의 대각은 원래부터 열의 길이 제곱이었다(1장).

방향 2 — x = (1, 1, 1). 이 경우 q는 G의 아홉 성분을 전부 더한 값이다.

대각   : 2 + 14 + 10 = 26
비대각 : 2·(1 + 3 + 3) = 2·7 = 14
합     : 26 + 14 = 40

검산: A(1,1,1) = c₁ + c₂ + c₃ = (1+1+0, 0+2+0, 1+0+3, 0+3+1) = (2, 2, 4, 4). 길이 제곱은 4 + 4 + 16 + 16 = 40. 일치한다. 10장에서 Pb = (2,2,4,4) = A(1,1,1)이었으니, 이 40은 그때 본 ‖Pb‖² = 40과 같은 숫자다(76 = 40 + 36의 그 40이다).

방향 3 — 고유방향 x = v₁ = (1, 5, 3). 고유벡터에서는 계산이 두 배로 쉬워진다. Gv = λv이므로

v^TGv = v^T(λv) = λ‖v‖² = 16 · 35 = 560 ‖v₁‖² = 1 + 25 + 9 = 35. 고유값 16은 14장에서 특성다항식으로 구했다. (15.7)

정의대로도 해 보자. Gv₁을 먼저 계산한다.

Gv₁ = ( 2·1 + 1·5 + 3·3 ,  1·1 + 14·5 + 3·3 ,  3·1 + 3·5 + 10·3 )
    = ( 2 + 5 + 9 , 1 + 70 + 9 , 3 + 15 + 30 )
    = ( 16 , 80 , 48 ) = 16 · (1, 5, 3)   ✓ 고유벡터 확인
v₁ᵀ(Gv₁) = 1·16 + 5·80 + 3·48 = 16 + 400 + 144 = 560   ✓

검산: Av₁ = (1+5, 10, 1+9, 15+3) = (6, 10, 10, 18), ‖Av₁‖² = 36 + 100 + 100 + 324 = 560. 세 번째도 일치한다. 그리고 560 = σ₁²‖v₁‖² = 16 · 35이므로 특이값 σ₁ = 4가 여기서도 보인다 — 16장에서 이것이 SVD의 출발점이 된다.

x	xTGx	Ax	‖Ax‖²	일치
(1, 0, 0)	2	(1, 0, 1, 0)	1+0+1+0 = 2	✓
(1, 1, 1)	40	(2, 2, 4, 4)	4+4+16+16 = 40	✓
(1, 5, 3)	16·35 = 560	(6, 10, 10, 18)	36+100+100+324 = 560	✓

세 방향 모두 양수다. 물론 세 개를 확인했다고 모든 방향에서 양수라는 보장은 없다. 그것을 보장해 주는 것이 다음의 세 판정이다.

(나) 세 판정을 전부 적용한다

판정 ① 고유값. 14장에서 특성다항식 det(G − λI) = −(λ−16)(λ−9)(λ−1)을 손으로 전개해 λ = 16, 9, 1을 얻었다. 셋 다 양수다. G ≻ 0.

검산은 두 줄이면 된다. 16 + 9 + 1 = 26 = 2 + 14 + 10 = tr G, 16 · 9 · 1 = 144 = det G. 둘 다 맞는다.

판정 ② 실베스터. 7장에서 이미 구했다.

d₁ = det[ 2 ] = 2
d₂ = det ⎡ 2   1 ⎤ = 2·14 − 1·1 = 28 − 1 = 27
         ⎣ 1  14 ⎦
d₃ = det G = 144

2 > 0, 27 > 0, 144 > 0. 세 개가 전부 양수이므로 실베스터 판정법에 의해 G ≻ 0이다. 7장에서 "15장에서 이름이 붙는다"고 예고한 바로 그 조건이며, 거기서 본 기하적 이유 — d_k는 A의 앞 k개 열이 만드는 k차원 부피의 제곱 — 가 실베스터 판정법의 내용을 그대로 말해 준다. 부피가 0이 아니면 양정치다.

주의할 점 하나. "선행주소행렬식이 전부 ≥ 0"은 양반정치의 판정이 아니다. diag(0, −1)은 d₁ = d₂ = 0이지만 양반정치가 아니다. 실베스터는 양정치 전용이다.

판정 ③ 콜레스키. G = LL^T를 손으로 구한다. L을 하삼각으로 두고 양변의 성분을 왼쪽 위부터 하나씩 맞춰 간다.

     ⎡ ℓ₁₁   0    0  ⎤          ⎡ ℓ₁₁  ℓ₂₁  ℓ₃₁ ⎤     ⎡  2   1   3 ⎤
 L = ⎢ ℓ₂₁  ℓ₂₂   0  ⎥    Lᵀ =  ⎢  0   ℓ₂₂  ℓ₃₂ ⎥  =  ⎢  1  14   3 ⎥
     ⎣ ℓ₃₁  ℓ₃₂  ℓ₃₃ ⎦          ⎣  0    0   ℓ₃₃ ⎦     ⎣  3   3  10 ⎦

1

(1,1): ℓ₁₁² = 2 → ℓ₁₁ = √2. 여기서 제곱근 안이 양수여야 한다.

2

(2,1): ℓ₂₁ℓ₁₁ = 1 → ℓ₂₁ = 1√2 = √22

3

(3,1): ℓ₃₁ℓ₁₁ = 3 → ℓ₃₁ = 3√2 = 3√22

4

(2,2): ℓ₂₁² + ℓ₂₂² = 14 → ℓ₂₂² = 14 − 12 = 272 → ℓ₂₂ = 3√62 (√(27/2) = 3√3/√2 = 3√6/2)

5

(3,2): ℓ₃₁ℓ₂₁ + ℓ₃₂ℓ₂₂ = 3. ℓ₃₁ℓ₂₁ = 3√2·1√2 = 32 이므로 ℓ₃₂ℓ₂₂ = 32 → ℓ₃₂ = 3/23√6/2 = 1√6 = √66

6

(3,3): ℓ₃₁² + ℓ₃₂² + ℓ₃₃² = 10 → ℓ₃₃² = 10 − 92 − 16 = 60 − 27 − 16 = 326 = 163 → ℓ₃₃ = 4√3 = 4√33

세 번 모두 제곱근 안이 양수였다. 콜레스키가 완주했다 — 따라서 G는 양정치다.

     ⎡  √2        0        0    ⎤        대각:  √2 ,  √(27/2) ,  √(16/3)
 L = ⎢  √2/2    3√6/2      0    ⎥              ≈ 1.414214, 3.674235, 2.309401
     ⎣ 3√2/2    √6/6    4√3/3   ⎦        det L = √2·√(27/2)·√(16/3) = √144 = 12

(15.8) det G = (det L)² = 144. 7장에서 본 144 = 12²가 여기서 제곱근이 벗겨진다.

검산은 LL^T를 두 개만 확인하면 충분하다. (2,2): (√2/2)² + (3√6/2)² = 12 + 544 = 12 + 272 = 14 ✓. (3,3): 92 + 16 + 163 = 27 + 1 + 326 = 606 = 10 ✓.

LU와 콜레스키는 같은 계산이다 — L_chol = L·√D

L의 대각을 다시 보자. √2, √(27/2), √(16/3). 괄호 안의 2, 272, 163는 5장과 7장에서 소거로 구한 LU의 피벗 그 자체다. 우연이 아니다.

대칭행렬에 소거를 하면 G = LU에서 U의 각 행에서 피벗을 뽑아낼 수 있고, 남는 것이 L^T가 된다. 즉 대칭행렬의 LU는 자동으로 LDL^T가 된다.

      ⎡  1    0    0 ⎤       ⎡ 2    0     0   ⎤        ⎡ 1  1/2  3/2 ⎤
 G =  ⎢ 1/2   1    0 ⎥   ·   ⎢ 0  27/2    0   ⎥   ·    ⎢ 0   1   1/9 ⎥
      ⎣ 3/2  1/9   1 ⎦       ⎣ 0    0   16/3  ⎦        ⎣ 0   0    1  ⎦
          L (대각 1)              D = 피벗                   Lᵀ

여기서 L의 성분 1/2, 3/2, 1/9는 5장 소거에서 쓴 승수 그대로다. D가 전부 양수이므로 √D = diag(√2, √(27/2), √(16/3))를 정의할 수 있고, D = √D·√D이니 괄호를 다시 묶으면

G = L D Lᵀ = (L√D)(√D Lᵀ) = (L√D)(L√D)ᵀ = L_chol · L_cholᵀ

즉 L_chol = L·√D다. 열마다 그 열의 피벗의 제곱근을 곱한 것이다. 확인해 보자. 2열은 (0, 1, 1/9)에 √(27/2) = 3√6/2를 곱하면 (0, 3√6/2, √6/6) — 식 (15.8)의 2열과 정확히 같다.

여기서 양정치가 왜 필요한지가 선명해진다. √D를 만들려면 피벗이 전부 양수여야 한다. 그리고 7장의 관계 피벗_k = d_k/d_k−1에 의해, 피벗이 전부 양수인 것과 선행주소행렬식이 전부 양수인 것은 같은 말이다. 판정 ②와 판정 ③은 사실 하나의 계산을 두 언어로 부른 것이다. 그리고 판정 ①은 식 (15.3)이 말하듯 같은 사실을 직교 좌표계에서 본 것이다. 세 판정이 동치인 이유가 여기까지다.

실무에서 콜레스키를 쓰는 이유도 분명하다. LU가 2n³/3번 연산을 쓰는 데 비해 대칭성을 이용하는 콜레스키는 n³/3번이면 된다. 절반이다. 게다가 양정치 행렬에서는 행 교환(피벗팅)이 필요 없다는 것이 증명되어 있어, 5장에서 걱정했던 수치 안정성 문제도 없다. 가우스 과정, 능형회귀, 칼만 필터, 그래프 라플라시안 — 대칭 양정치가 나오는 곳에서는 거의 항상 콜레스키가 돈다.

(다) 등고면은 타원체다 — 조건수가 곧 납작한 정도

이차형식의 모양을 보는 가장 좋은 방법은 같은 값을 주는 점들을 모으는 것이다. 등고선을 그리는 것과 같다. x^TGx = 1인 점들의 집합을 보자. 식 (15.3)에서 좌표를 고유벡터 방향으로 돌리면

16y₁² + 9y₂² + y₃² = 1 y = Q^Tx, Q의 열은 정규화한 고유벡터 v₁/√35, v₂/√14, v₃/√10. (15.9)

이것은 표준형 타원체 방정식 y₁²/a² + y₂²/b² + y₃²/c² = 1에서 a² = 1/16, b² = 1/9, c² = 1인 경우다. 즉 반축 길이는

주축 방향            반축 길이 = 1/√λ        고유값
v₁ = (1, 5, 3)        1/√16 = 1/4             16    ← 가장 짧은 축
v₂ = (1, −2, 3)       1/√9  = 1/3              9
v₃ = (3, 0, −1)       1/√1  = 1                1    ← 가장 긴 축

주의해서 읽을 곳이 여기다. 고유값이 큰 방향일수록 축이 짧다. 직관적으로도 맞는다 — 그 방향으로는 q가 빨리 커지므로, q = 1이라는 문턱에 금방 도달한다. "가파른 방향 = 좁은 골짜기"다.

가장 긴 축과 가장 짧은 축의 비를 보자.

1/√λ_min1/√λ_max = √λ_max√λ_min = √16√1 = 4 = √κ₂(G) = κ₂(A) (15.10)

8장의 조건수가 타원체의 납작한 정도로 나타난 것이다. κ₂(G) = 16이지만 타원체의 축비는 그 제곱근인 4, 즉 κ₂(A) = 4다. 우리 A는 조건수가 4라 그림이 적당히 찌그러진 정도지만, 조건수가 1000인 문제라면 축비가 1000인 바늘 같은 타원체가 된다. 경사하강법이 그런 골짜기에서 지그재그로 튀는 이유가 이 그림 한 장에 들어 있다.

짝이 되는 그림도 기억해 두자. 8장·16장의 상 타원체 {Ax : ‖x‖=1}은 반축이 σ = 4, 3, 1이다. 이번 장의 x^TGx = 1은 ‖Ax‖ = 1과 같은 말이므로 반축이 1/σ인 역수 타원체다.

도해 3. 같은 A가 만드는 두 타원체를 가장 납작한 단면(v₁–v₃ 평면)에서 잘라 본 것. 왼쪽은 이번 장의 등고면 x^TGx = 1, 오른쪽은 8장·16장의 상 타원체다. 서로 90° 돌아간 모양인 것은 반축이 1/σ와 σ로 역수 관계이기 때문이다. 축비는 양쪽 다 4이며 이것이 조건수다.

STEP 04

코드 — 판정 세 개를 직접 구현하고, 깨뜨려 본다

세 판정을 각각 바닥부터 짜서 G에 적용하고, 콜레스키는 scipy와 소수점까지 대조한다.

positive_definite.py — 세 판정 직접 구현 + 콜레스키 대조

import numpy as np
from scipy.linalg import cholesky
np.set_printoptions(precision=6, suppress=True)
A = np.array([[1,1,0],[0,2,0],[1,0,3],[0,3,1]], float)
G = A.T @ A
def chol(S):                                   # 콜레스키 직접 구현
    n = len(S); L = np.zeros((n, n))
    for i in range(n):
        for j in range(i+1):
            s = S[i,j] - L[i,:j] @ L[j,:j]
            if i == j:
                if s <= 1e-10:                 # 대각이 양수가 아니면 양정치가 아니다
                    raise np.linalg.LinAlgError("%d번째 대각이 %+.2e" % (i+1, s))
                L[i,i] = np.sqrt(s)
            else:
                L[i,j] = s / L[j,j]
    return L
def verdict(name, S):                          # 세 판정을 각각 적용
    ev = np.linalg.eigvalsh(S)
    lpm = [np.linalg.det(S[:k,:k]) for k in range(1, len(S)+1)]
    print("[%s]" % name)
    print("  1 고유값         %-28s  모두>0 : %s" % (np.round(ev,6).tolist(), np.all(ev > 1e-10)))
    print("  2 선행주소행렬식 %-28s  모두>0 : %s" % (np.round(lpm,6).tolist(), np.all(np.array(lpm) > 1e-10)))
    try:
        L = chol(S); print("  3 콜레스키       성공, diag(L) =", np.round(np.diag(L),6).tolist())
    except np.linalg.LinAlgError as e:
        print("  3 콜레스키       LinAlgError:", e)
verdict("G = AᵀA", G)
print("  → 세 판정이 모두 '양정치'로 일치\n")
for x in ([1,0,0], [1,1,1], [1,5,3]):          # 이차형식 = ‖Ax‖²
    x = np.array(x, float); Ax = A @ x
    print("  x=%-10s xᵀGx =%7.1f   ‖Ax‖² =%7.1f   Ax = %s"
          % (tuple(int(t) for t in x), x @ G @ x, Ax @ Ax, [int(t) for t in Ax]))
L, Ls = chol(G), cholesky(G, lower=True)       # scipy 와 대조
print("\n직접 구현 L =\n", np.round(L, 6))
print("최대 차이(vs scipy) :", abs(L - Ls).max(), "   ‖LLᵀ − G‖max :", abs(L@L.T - G).max())
print("손계산 대각 √2, √(27/2), √(16/3) =", [round(float(np.sqrt(v)), 6) for v in (2, 27/2, 16/3)])

[G = AᵀA]
  1 고유값         [1.0, 9.0, 16.0]              모두>0 : True
  2 선행주소행렬식 [2.0, 27.0, 144.0]            모두>0 : True
  3 콜레스키       성공, diag(L) = [1.414214, 3.674235, 2.309401]
  → 세 판정이 모두 '양정치'로 일치
  x=(1, 0, 0)  xᵀGx =    2.0   ‖Ax‖² =    2.0   Ax = [1, 0, 1, 0]
  x=(1, 1, 1)  xᵀGx =   40.0   ‖Ax‖² =   40.0   Ax = [2, 2, 4, 4]
  x=(1, 5, 3)  xᵀGx =  560.0   ‖Ax‖² =  560.0   Ax = [6, 10, 10, 18]
직접 구현 L =
[[1.414214 0.       0.      ]
[0.707107 3.674235 0.      ]
[2.12132  0.408248 2.309401]]
최대 차이(vs scipy) : 0.0    ‖LLᵀ − G‖max : 4.440892098500626e-16
손계산 대각 √2, √(27/2), √(16/3) = [1.414214, 3.674235, 2.309401]

고유값 16, 9, 1, 선행주소행렬식 2, 27, 144, 콜레스키 대각 1.414214, 3.674235, 2.309401. 손계산의 √2, √(27/2), √(16/3)과 마지막 자리까지 같고, 직접 구현과 scipy의 차이는 정확히 0이다. 세 방향의 이차형식도 2, 40, 560으로 종이와 일치한다.

이제 깨뜨려 본다. 3장에서 만든 랭크결손 변형 A'(세 번째 열을 c₃' = c₁ + c₂로 바꾼 가상의 행렬)의 그람 행렬은 양반정치이지 양정치가 아니어야 한다. 세 판정이 전부 같은 방식으로 실패하는지 본다. 이어서 표본공분산을 중심화 행렬로 실제로 만들어 성질을 확인한다.

semidefinite.py — 양정치가 깨지는 경우와 표본공분산

import numpy as np
from scipy.linalg import cholesky
from positive_definite import A, chol, verdict   # 위 코드의 함수를 재사용
# ① 양정치가 깨지는 경우 — 3장의 랭크결손 변형 A' (c₃' = c₁ + c₂)
Ad = np.column_stack([A[:,0], A[:,1], A[:,0] + A[:,1]])
Gd = Ad.T @ Ad
print("A' 의 열 :", [[int(t) for t in c] for c in Ad.T], "   G' =", Gd.astype(int).tolist())
verdict("G' = A'ᵀA'  (열이 종속)", Gd)
try:
    cholesky(Gd, lower=True)
except np.linalg.LinAlgError as e:
    print("  scipy.linalg.cholesky →", type(e).__name__, ":", e)
z = np.array([1., 1., -1.])
print("  영방향 z = (1,1,−1) :  A'z =", [int(t) for t in Ad @ z],
      "  zᵀG'z =", z @ Gd @ z, " → 0 이 되는 z≠0 이 있다 = 양반정치\n")
# ② 표본공분산은 왜 항상 양반정치인가 — 중심화 행렬 H 로 실측
def cov(X):
    n = len(X); H = np.eye(n) - np.ones((n, n)) / n     # 중심화 행렬
    Xc = H @ X
    return Xc.T @ Xc / (n - 1), Xc
S, Xc = cov(A)
print("X_c 의 열합 :", np.round(Xc.sum(0), 12).tolist(), " (중심화 확인)")
print("3S = X_cᵀX_c =", (3*S).astype(int).tolist())
verdict("표본공분산 S   n=4 > d=3", S)
S3, _ = cov(A[:3])                                      # 표본 3개, 특징 3개
print("\nn=3 < d=3 → rank(X_c) ≤ n−1 = 2 이므로 0 고유값이 반드시 생긴다")
verdict("표본공분산 S₃  n=3 ≤ d=3", S3)
print("  rank(S₃) =", np.linalg.matrix_rank(S3), " 대칭성 ‖S₃−S₃ᵀ‖max =", abs(S3 - S3.T).max())

A' 의 열 : [[1, 0, 1, 0], [1, 2, 0, 3], [2, 2, 1, 3]]    G' = [[2, 1, 3], [1, 14, 15], [3, 15, 18]]
[G' = A'ᵀA'  (열이 종속)]
  1 고유값         [0.0, 2.577795, 31.422205]    모두>0 : False
  2 선행주소행렬식 [2.0, 27.0, 0.0]              모두>0 : False
  3 콜레스키       LinAlgError: 3번째 대각이 +3.55e-15
  scipy.linalg.cholesky → LinAlgError : 3-th leading minor of the array is not positive definite
  영방향 z = (1,1,−1) :  A'z = [0, 0, 0, 0]   zᵀG'z = 0.0  → 0 이 되는 z≠0 이 있다 = 양반정치
X_c 의 열합 : [0.0, 0.0, 0.0]  (중심화 확인)
3S = X_cᵀX_c = [[1, -2, 1], [-2, 5, -3], [1, -3, 6]]
[표본공분산 S   n=4 > d=3]
  1 고유값         [0.053091, 0.922708, 3.024201]  모두>0 : True
  2 선행주소행렬식 [0.333333, 0.111111, 0.148148]  모두>0 : True
  3 콜레스키       성공, diag(L) = [0.57735, 0.57735, 1.154701]
n=3 < d=3 → rank(X_c) ≤ n−1 = 2 이므로 0 고유값이 반드시 생긴다
[표본공분산 S₃  n=3 ≤ d=3]
  1 고유값         [0.0, 0.402832, 3.930501]     모두>0 : False
  2 선행주소행렬식 [0.333333, 0.083333, -0.0]    모두>0 : False
  3 콜레스키       LinAlgError: 3번째 대각이 -4.44e-16
  rank(S₃) = 2  대칭성 ‖S₃−S₃ᵀ‖max = 0.0

읽을 것이 많다. 하나씩 보자.

관측	무슨 뜻인가
G'의 고유값에 0이 하나	세 열이 2차원만 만든다(3장의 랭크 2). 고유값 0의 고유벡터가 바로 영방향 z = (1,1,−1)이다.
d₁, d₂는 2, 27 그대로, d₃만 0	앞 두 열은 건드리지 않았으므로 G와 왼쪽 위 2×2가 같다. 세 번째에서 부피가 무너진다.
콜레스키가 세 번째 대각에서 멈춤	내 구현은 +3.55e−15(수치적 0)에서, scipy는 "3-th leading minor"에서 같은 지점을 지목한다.
A'z = 0, zTG'z = 0	식 (15.4)의 등호 조건이 실제로 성립한다. 양반정치의 정의 그 자체다.
3Σ̂ = XcTXc가 정수 행렬	[[1,−2,1],[−2,5,−3],[1,−3,6]]. 중심화해도 우리 자료는 분수로 흩어지지 않는다.
Σ̂₃에서 고유값 하나가 0	n = 3 표본으로는 d = 3차원을 채울 수 없다. rank = 2 = n − 1.

마지막 줄이 실무에서 가장 자주 부딪히는 지점이다. 유전체 데이터처럼 n = 100 표본에 d = 20000 특징인 경우, 표본공분산은 반드시 특이행렬이다. 역행렬이 없으니 마할라노비스 거리도 선형판별분석도 그대로는 돌지 않는다. 그래서 대각에 λI를 더해 고유값을 밀어 올리는 축소 추정(shrinkage)을 쓰는데, 그것이 정확히 20장에서 볼 능형회귀의 A^TA + λI와 같은 조작이다 — 16, 9, 0을 16+λ, 9+λ, λ로 만들어 양반정치를 양정치로 바꾸는 것이다.

왜 공분산행렬은 항상 대칭 양반정치인가

자료행렬 X가 (n, d)일 때, 중심화 행렬 H = I − 1n11^T로 열평균을 빼면 X_c = HX이고, 표본공분산은

Σ̂ = X_cᵀ X_c / (n − 1)

분모 n−1은 양수이므로 부호에 아무 영향이 없다. 남는 것은 X_c^TX_c — A^TA 꼴이다. 그러므로 식 (15.4)가 그대로 적용된다.

1. 항상 대칭이다. (X_c^TX_c)^T = X_c^TX_c. 계산 순서를 뒤집어도 같은 행렬이다. 그래서 cov(x_i, x_j) = cov(x_j, x_i)가 정의가 아니라 구조에서 나온다.
2. 고유값이 전부 ≥ 0이다. w^TΣ̂w = ‖X_cw‖²/(n−1) ≥ 0. 통계 언어로 읽으면 왼쪽은 "가중치 w로 만든 합성변수 X_cw의 표본분산"이다. 분산이 음수일 수 없다는 사실이 곧 양반정치다. 「통계학 기본」 3장에서 분산의 정의로 배운 것이 여기서 행렬의 성질이 된다.
3. 고유벡터가 서로 직교한다. 대칭이므로 14장의 스펙트럼 정리가 적용되어 Σ̂ = QΛQ^T, Q^TQ = I로 쓸 수 있다.

17장 PCA가 성립하는 근거가 이 셋 전부다. PCA는 "분산이 가장 큰 방향"을 찾는 문제인데, 그 방향이 최대 고유값의 고유벡터라는 것은 (2)의 이차형식 최대화에서 나오고, 주성분들이 서로 무상관(직교)이라는 것은 (3)에서 나오며, 설명 분산 비율 λ_i/Σλ를 비율로 읽을 수 있는 것은 (2)에 의해 λ_i ≥ 0이라 분모가 절댓값 걱정 없이 합으로 쓰이기 때문이다. 고유값이 음수가 될 수 있다면 "설명 분산 120%" 같은 말이 나온다.

표본 수가 모자라면 양정치가 될 수 없다

rank(Σ̂) = rank(X_c) ≤ min(n − 1, d)다. n−1인 이유는 중심화가 각 열에서 평균을 빼면서 "모든 행의 합이 0"이라는 제약 하나를 만들기 때문이다 — H는 1 방향을 죽이는 투영행렬이고 랭크가 n−1이다(10장).

따라서 n ≤ d이면 rank(Σ̂) ≤ n − 1 < d이므로 반드시 0 고유값이 생긴다. 위 출력의 Σ̂₃가 그 예다: n = 3, d = 3, 랭크 2, 고유값 0, 0.402832, 3.930501. 양반정치이되 양정치가 아니다.

우리 A는 n = 4 > d = 3이고 중심화 후에도 열이 독립이라 Σ̂가 양정치였다(고유값 0.053091, 0.922708, 3.024201). 다만 가장 작은 고유값이 0.053으로 작아 Σ̂의 조건수는 약 57이다 — 원래 G의 16보다 나쁘다. 중심화가 열의 공통 성분을 빼내면서 남은 방향들이 더 가까워진 것이다. 「통계학 기본」 16장의 상관행렬, 「손으로 푸는 회귀분석」 7장의 VIF가 재는 것이 바로 이 "얼마나 가까운가"다.

STEP 05

시각화 — 콜레스키가 하는 일, 그리고 PCA로 가는 길

먼저 콜레스키를 그림으로 보자. S = LL^T는 "대칭 양정치 행렬을 삼각행렬 하나로 쪼개는 것"인데, 그 L이 무엇을 하는 사상인지가 핵심이다.

도해 4. 위: LU 소거의 피벗에 제곱근을 씌우면 콜레스키 인자가 된다. 아래: L은 상관 없는 단위 잡음을 원하는 공분산으로 바꾸는 변환이다. 몬테카를로에서 상관 있는 정규난수를 만들 때, 가우스 과정에서 사전분포를 표집할 때 쓰는 것이 정확히 이 L이다. 반대로 L⁻¹x는 표본을 되돌려 백색화한다 — 마할라노비스 거리 x^TG⁻¹x = ‖L⁻¹x‖²가 그 뜻이다.

마지막으로 이 장이 앞뒤로 무엇과 연결되는지를 한 장에 모은다. 출발점은 식 (15.4) 하나다.

도해 5. 식 (15.4) 한 줄에서 갈라져 나오는 결론들. 왼쪽 줄기는 최적화(해가 유일한가), 가운데 줄기는 분해(어떻게 쪼개는가), 오른쪽 줄기는 통계(공분산이 왜 얌전한가)다. 세 줄기가 16·17·19장에서 각각 한 장씩 차지한다. 이 장은 그 셋의 공통 뿌리다.

이것만 기억하자

1. A^TA는 언제나 양반정치다. x^T(A^TA)x = ‖Ax‖² ≥ 0 — 증명이 한 줄이다. 등호는 Ax = 0일 때뿐이므로, 열이 독립이면 양정치가 된다. 이 시리즈에서 G가 늘 얌전했던 이유(det = 144, 고유값 16,9,1, 피벗 전부 양수)가 전부 여기서 나온다.
2. 세 판정은 하나의 사실을 세 언어로 말한 것이다. 고유값 16, 9, 1 > 0(14장) = 선행주소행렬식 2, 27, 144 > 0(7장) = 콜레스키 L의 대각 √2, √(27/2), √(16/3)가 실수(5장의 피벗에 제곱근을 씌운 것). 실무에서는 셋 중 가장 싼 콜레스키가 판정기 노릇을 한다.
3. 공분산은 A^TA 꼴이므로 대칭 · 양반정치 · 직교 고유벡터를 공짜로 얻는다. Σ̂ = X_c^TX_c/(n−1)이고, w^TΣ̂w는 합성변수의 표본분산이라 음수가 될 수 없다. 17장 PCA의 근거가 전부 이 세 줄이며, n ≤ d이면 랭크가 n−1 이하라 반드시 0 고유값이 생긴다.

흔한 오해

1. "대각이 전부 양수면 양정치다." — 아니다. 필요조건일 뿐이다. [[1, 2], [2, 1]]은 대각이 둘 다 1이지만 d₂ = 1 − 4 = −3 < 0이라 부정치다. x = (1, −1)을 넣으면 1 − 4 + 1 = −2가 나온다. 대각만 보는 습관이 공분산 추정에서 사고를 만든다.
2. "행렬식이 양수면 양정치다." — 아니다. det는 고유값의 곱이라 음수가 짝수 개면 양수가 된다. diag(−1, −1)은 det = 1 > 0이지만 모든 방향에서 음수다. 실베스터가 모든 선행주소행렬식을 요구하는 이유가 이것이다.
3. "양반정치면 고유값이 0인 방향은 데이터에 없는 방향이다." — 방향이 없는 것이 아니라 그 방향의 변동이 0이다. 우리 A'의 z = (1,1,−1)은 멀쩡한 방향이지만 A'z = 0이라 아무 신호도 만들지 않는다. 회귀에서 이것이 "계수를 결정할 수 없다"는 뜻이 되고, 12장의 다중공선성과 같은 사건이다.
4. "콜레스키가 실패하면 데이터가 잘못된 것이다." — 대개는 정상이다. 표본 수가 특징 수보다 적거나(n ≤ d), 중복 특징이 있거나(원-핫 인코딩에서 범주를 다 넣는 흔한 실수), 반올림 오차로 마지막 고유값이 −10⁻¹⁶이 된 경우다. 위 출력의 −4.44e−16이 그 예다. 대각에 아주 작은 λ를 더하는 처방(지터·축소·능형)이 바로 양반정치를 양정치로 미는 조작이며, 20장에서 16, 9, 1 → 16+λ, 9+λ, 1+λ로 다시 본다.

이전 — 14장 대칭행렬의 스펙트럼 분해 다음 — 16장 SVD — 모든 행렬의 정규분해