torch.randn이 주는 것은 좌표끼리 아무 관계가 없는 둥근 구름입니다. 실제 데이터는 기운 타원입니다. 그 기울기를 재는 도구가 공분산 행렬이고, 촐레스키 분해로 그 행렬을 되돌려 상관 있는 표본을 직접 만들어 봅니다.
PALDYN Team//20 MIN READ
23MATH
중급
공분산 · 다변량 정규분포
확산 모델의 학습 코드를 열면 noise = torch.randn_like(x) 한 줄이 나옵니다. 이미지에 얹을 노이즈를 뽑는 자리입니다. 이 함수가 주는 것은 모든 좌표가 서로 아무 관계도 없는 노이즈입니다 — 픽셀 하나가 크게 나왔다고 옆 픽셀이 덩달아 커지지 않습니다.
그런데 진짜 이미지의 픽셀은 그렇지 않습니다. 옆 픽셀은 거의 같은 값입니다. 임베딩도 마찬가지여서, 512차원 벡터를 몇만 개 모아 놓으면 512개 방향으로 고르게 퍼지지 않고 몇몇 방향으로만 길쭉하게 퍼져 있습니다.
지난 글에서 정규분포의 손잡이가 μ 와 σ 둘뿐이라고 했습니다. 좌표가 d 개로 늘어나면 «퍼진 정도»는 수 하나로는 부족합니다. 방향마다 다르게 퍼져 있으니까요. 이 글은 그 «방향마다 다른 퍼짐»을 담는 행렬을 세우고, 마지막에는 그 행렬을 거꾸로 써서 randn의 둥근 구름을 원하는 모양의 타원으로 바꾸는 법까지 갑니다.
분산 하나로는 부족하다
좌표가 둘인 자료를 생각해 봅시다. x 의 분산과 y 의 분산을 각각 구하면 «각 축 방향으로 얼마나 퍼졌는가»는 알 수 있습니다. 하지만 그것만으로는 아래 세 그림을 구별하지 못합니다 — 셋 다 축 방향 퍼짐은 같을 수 있습니다.
빠진 것은 «두 좌표가 같이 움직이는가»입니다. 그것을 재는 값이 공분산입니다.
정의. 두 확률변수 X,Y 의 공분산은
Cov(X,Y)=E[(X−μX)(Y−μY)]
말로 풀면 이렇습니다. 각 표본에서 «x 가 자기 평균보다 얼마나 큰가»와 «y 가 자기 평균보다 얼마나 큰가»를 곱하고, 그 곱의 평균을 냅니다. 둘이 함께 평균 위에 있거나 함께 평균 아래에 있으면 곱이 양수, 엇갈리면 음수입니다. 그러니 양수면 같이 움직이고 음수면 반대로 움직입니다.
Y 자리에 X 를 그대로 넣으면 Cov(X,X)=E[(X−μX)2], 즉 분산입니다. 공분산은 분산을 두 변수로 늘린 것이고 분산은 자기 자신과의 공분산입니다.
다섯 점으로 손으로 세어 보기
자료가 다섯 개 있다고 합시다: (5,3), (5,11), (10,10), (15,9), (15,17).
평균은 xˉ=(5+5+10+15+15)/5=10, yˉ=(3+11+10+9+17)/5=10 으로 둘 다 10입니다.
편차를 적고 곱합니다.
점
x−xˉ
y−yˉ
곱
(5,3)
−5
−7
+35
(5,11)
−5
+1
−5
(10,10)
0
0
0
(15,9)
+5
−1
−5
(15,17)
+5
+7
+35
곱의 합은 35−5+0−5+35=60 입니다. 표본에서 잰 값이므로 지난 글에서 본 이유로 n 이 아니라 n−1=4 로 나눕니다.
sxy=460=15
같은 방식으로 각 축의 분산도 나옵니다. x 편차의 제곱합은 25+25+0+25+25=100 이므로 sx2=100/4=25 이고, y 편차의 제곱합도 49+1+0+1+49=100 이라 sy2=25 입니다.
공분산 행렬 — 값 네 개를 한 표에
값 세 개(sx2, sy2, sxy)가 나왔으니 표로 묶습니다.
정의. 벡터 확률변수 x∈Rd 의 공분산 행렬은
Σ=E[(x−μ)(x−μ)⊤],Σij=Cov(xi,xj)
(x−μ) 는 d×1, 그 전치는 1×d 이므로 곱하면 d×d 행렬이 나옵니다. 그 (i,j) 칸이 정확히 i 번 좌표와 j 번 좌표의 공분산입니다.
우리 자료라면 이렇습니다.
Σ=(25151525)
대각선에는 분산이, 나머지 자리에는 공분산이 앉습니다. 그리고 Cov(xi,xj)=Cov(xj,xi) 이므로 Σ 는 언제나 대칭 행렬입니다. 대칭 행렬을 다룬 글에서 「AI에서 만나는 행렬은 왜 대부분 대칭인가」를 물었는데, 그 답의 절반이 여기 있습니다.
대칭인 것으로 끝이 아니다 — 양반정치성
Σ 에는 성질이 하나 더 붙습니다. 아무 방향 a 나 잡아 이차형식을 계산해 봅시다.
a⊤Σa=a⊤E[(x−μ)(x−μ)⊤]a=E[(a⊤(x−μ))2]=Var(a⊤x)
가운데 단계는 a⊤v⋅v⊤a=(a⊤v)2 이라는 것뿐입니다. 그런데 오른쪽은 어떤 실수 확률변수의 분산이므로 절대 음수일 수 없습니다.
a⊤Σa=Var(a⊤x)≥0모든a에대해
모든 방향에서 이차형식이 0 이상인 대칭 행렬을 양반정치(positive semi-definite)라고 합니다. 즉 공분산 행렬은 대칭이면서 양반정치입니다. 스펙트럼 정리에 따라 대칭 행렬은 고유값이 전부 실수이고 서로 직교하는 고유벡터를 가지는데, 여기에 양반정치성이 붙으면 고유값이 전부 0 이상이라는 것까지 따라옵니다.
이 사실은 곱씹을 값어치가 있습니다. a 를 단위벡터로 잡으면 a⊤Σa 는 «그 방향으로 자료를 사영했을 때의 분산»입니다. 그러니 Σ 의 가장 큰 고유값에 딸린 고유벡터가 자료가 가장 넓게 퍼진 방향이고, 가장 작은 고유값 쪽이 가장 납작한 방향입니다. 임베딩 벡터를 모아 Σ 를 만들고 고유값을 정렬해 보면 앞쪽 몇 개가 유난히 큰 것 — 서두에 말한 «몇몇 방향으로만 길쭉하다»가 바로 이 숫자로 나타납니다.
상관계수 — 단위를 지운 공분산
공분산에는 불편한 점이 하나 있습니다. 단위에 딸려 다닙니다. 길이를 센티미터로 재다가 밀리미터로 바꾸면 공분산이 100배가 됩니다. 관계의 세기는 그대로인데 숫자만 커진 것입니다.
그래서 각자의 표준편차로 나눠 단위를 지웁니다.
정의.상관계수는 ρxy=σXσYCov(X,Y) 이고 항상 −1≤ρ≤1 이다.
우리 자료라면 ρ=15/(5×5)=0.6 입니다.
범위가 [−1,1] 인 것은 코시-슈바르츠 부등식에서 나옵니다. 앞에서 본 a⊤Σa≥0 에 a=(t,1) 을 넣으면 t 에 대한 이차식 sx2t2+2sxyt+sy2≥0 이 되고, 이차식이 항상 0 이상이려면 판별식이 0 이하여야 하므로 sxy2≤sx2sy2 입니다. 양변에 제곱근을 씌우면 그대로 ∣ρ∣≤1 입니다.
ρ=±1 은 판별식이 0인 경우, 즉 두 좌표가 완전히 직선 위에 놓인 경우입니다. ρ 가 재는 것은 «직선 관계»뿐이라는 뜻이기도 합니다 — 원 둘레를 따라 놓인 점들은 분명히 관계가 있지만 ρ 는 0에 가깝게 나옵니다.
다변량 정규분포
이제 분포로 갑니다. 좌표가 하나일 때 정규분포의 밀도는 exp(−(x−μ)2/2σ2) 에 상수를 곱한 것이었습니다. 지수 안은 «중심에서 떨어진 거리를 표준편차로 잰 값의 제곱»입니다.
d 차원으로 늘릴 때 바뀌는 것은 그 «거리를 재는 자» 하나뿐입니다.
정의. 평균 μ, 공분산 Σ 인 다변량 정규분포의 밀도는
f(x)=(2π)d/2∣Σ∣1/21exp(−21(x−μ)⊤Σ−1(x−μ))
∣Σ∣ 는 행렬식이고, 앞의 분수는 전체 적분이 1이 되도록 높이를 맞추는 정규화 상수입니다. 볼 곳은 지수 안입니다.
(x−μ)2/σ2 자리에 (x−μ)⊤Σ−1(x−μ) 가 들어갔습니다. 나누기가 역행렬 곱하기로 바뀐 것입니다. 이 값의 제곱근을 마할라노비스 거리라고 하고, «Σ 라는 자로 잰 중심까지의 거리»라는 뜻입니다.
등고선이 왜 타원인가
밀도가 같은 점들을 모으면 지수 안이 같은 점들을 모으는 것과 같습니다. 즉 등고선은
(x−μ)⊤Σ−1(x−μ)=c
를 만족하는 점들입니다. 그런데 이 왼쪽은 대칭 행렬 Σ−1 의 이차형식이고, 스펙트럼 정리는 대칭 행렬의 이차형식이 고유기저에서 보면 교차항 없이 제곱의 합이 된다고 말합니다. 고유벡터 방향의 좌표를 u1,u2 라 하면 Σ 의 고유값이 λ1,λ2 일 때 Σ−1 의 고유값은 1/λ1,1/λ2 이므로
λ1u12+λ2u22=c
가 됩니다. 이것이 타원의 방정식이고, 반지름은 λ1c 와 λ2c 입니다. 등고선이 타원인 이유는 그러니까 하나입니다 — Σ 가 대칭이라서입니다.
우리 Σ 로 확인해 봅시다. (25151525) 의 고유값은 대각합이 50이고 행렬식이 625−225=400 이므로 λ2−50λ+400=0 에서 λ=40,10 입니다. 고유벡터는 (1,1)/2 와 (1,−1)/2 입니다.
c=1 인 등고선의 반지름이 40≈6.32 와 10≈3.16 입니다. 두 좌표의 분산이 25로 같았기 때문에 타원의 축이 정확히 45도로 섰습니다 — 분산이 서로 다르면 축도 그만큼 기웁니다.
c 를 키우면 같은 모양이 그대로 커집니다. 그래서 다변량 정규분포의 등고선은 닮은 타원들이 겹겹이 놓인 모양입니다.
촐레스키 분해 — 행렬의 제곱근
여기까지가 «주어진 자료의 퍼짐을 읽는» 방향이었습니다. 이제 반대로 갑니다. Σ 를 정해 놓고 그 공분산을 가진 표본을 만들려면?
한 차원일 때는 쉬웠습니다. z∼N(0,1) 을 뽑아 σ 를 곱하면 분산이 σ2 입니다. 여러 차원에서도 같은 발상을 쓰려면 «Σ 의 제곱근»에 해당하는 행렬 L, 즉 LL⊤=Σ 인 L 이 필요합니다.
먼저 그런 L 이 있으면 실제로 통한다는 것부터 확인합니다. z 의 공분산이 단위행렬 I 일 때 x=Lz 의 공분산은
Cov(Lz)=E[Lzz⊤L⊤]=LE[zz⊤]L⊤=LIL⊤=LL⊤=Σ
입니다. L 은 상수 행렬이라 기댓값 밖으로 나오고, 남은 E[zz⊤] 이 정의상 z 의 공분산 행렬인 I 입니다. 게다가 정규분포는 선형변환에 닫혀 있으므로 Lz 는 여전히 정규분포입니다. 필요한 것은 L 을 구하는 방법 하나뿐입니다.
정의. 대칭 양정치 행렬 Σ 를 Σ=LL⊤ 으로 쓰는 것을 촐레스키 분해라 하고, 이때 L 은 대각 성분이 양수인 아래삼각행렬로 유일하게 정해진다.
아래삼각으로 잡는 것이 요령입니다. 미지수가 줄어 한 칸씩 차례로 풀리기 때문입니다. 2×2 로 직접 해 봅시다.
세 걸음에서 전부 «앞에서 구한 값을 대입하고 제곱근을 하나 뽑는» 일만 했습니다. d 차원에서도 똑같이 위에서 아래로, 왼쪽에서 오른쪽으로 한 칸씩 채워 나갑니다. 3단계에서 제곱근 안이 음수가 되지 않는 것이 바로 양정치성이 보장해 주는 것이고, 그래서 공분산 행렬에는 촐레스키 분해가 언제나 존재합니다.
이제 z=(z1,z2) 를 randn으로 뽑고 곱하면 됩니다.
x=Lz=(5z13z1+4z2)
둘째 좌표가 첫째 좌표와 z1 을 나눠 쓰고 있습니다 — 상관이 생기는 자리가 눈에 보입니다.
코드로 확인하기
import numpy as npSigma = np.array([[25.0, 15.0], [15.0, 25.0]])mu = np.array([10.0, 10.0])L = np.linalg.cholesky(Sigma) # 아래삼각 Lprint(L)# [[5. 0.]# [3. 4.]]rng = np.random.default_rng(0)z = rng.standard_normal((200_000, 2)) # 독립인 둥근 구름x = z @ L.T + mu # 행이 표본이라 L.T 를 오른쪽에 곱한다print(np.cov(x.T).round(2))# [[25.02 15.03]# [15.03 25.03]]print(np.corrcoef(x.T)[0, 1].round(3)) # 0.601print(np.linalg.eigvalsh(np.cov(x.T)).round(2)) # [ 9.99 40.06]
표본에서 되돌려 잰 공분산이 넣어 준 Σ 와 맞고, 고유값도 손으로 구한 10과 40입니다. L 을 그대로 찍어 보면 손계산과 같은 [[5,0],[3,4]] 입니다.
반대 방향도 한 줄입니다. z=L−1(x−μ) 로 되돌리면 공분산이 I 인 표본이 나옵니다 — 기운 타원을 다시 둥근 구름으로 펴는 이 조작을 백색화(whitening)라고 합니다.
w = np.linalg.solve(L, (x - mu).T).Tprint(np.cov(w.T).round(2))# [[1. 0.]# [0. 1.]]
정리
공분산은 두 좌표가 같이 움직이는 정도이고, 자기 자신과의 공분산이 분산이다.
공분산 행렬 Σ 는 언제나 대칭이고 양반정치다. a⊤Σa=Var(a⊤x)≥0 한 줄에서 나온다.
그래서 가장 큰 고유값의 방향이 자료가 가장 넓게 퍼진 방향이다. 임베딩의 이방성이 이 숫자로 보인다.
상관계수는 단위를 지운 공분산이고 범위 [−1,1] 은 코시-슈바르츠에서 나온다. 재는 것은 직선 관계뿐이다.
다변량 정규분포에서 바뀐 것은 거리를 재는 자 하나다. (x−μ)2/σ2 이 (x−μ)⊤Σ−1(x−μ) 이 되었을 뿐이다.
등고선이 타원인 것은 Σ 가 대칭이기 때문이고, 반지름은 고유값의 제곱근이다.
촐레스키 분해 Σ=LL⊤ 이 행렬의 제곱근 노릇을 한다. Lz 로 상관을 만들고 L−1(x−μ) 로 지운다.
서두의 torch.randn_like로 돌아가 봅시다. 그 한 줄이 주는 것은 Σ=I 인 표본입니다. 확산 모델이 그것으로 충분한 이유는 노이즈를 일부러 좌표마다 독립으로 넣기 때문이고, 반대로 데이터의 실제 분포를 흉내 내야 하는 자리에서는 L 을 곱하는 한 걸음이 반드시 끼어듭니다. VAE가 대각 성분만 학습하는 것도, 전체 Σ 대신 «촐레스키의 아래삼각 L 을 학습한다»는 설계가 종종 등장하는 것도 이 글의 식에서 곧장 읽힙니다.
다음은 분포를 «읽는» 쪽에서 «만드는» 쪽으로 한 번 더 옮겨 갈 차례입니다. 모델이 내놓는 것은 점수 벡터일 뿐인데 그것을 확률로 바꿔 주는 함수가 하나 있고, 그 함수의 모양이 왜 하필 그것인지를 유도합니다.
실수를 2^b개 격자에 사상할 때 오차의 분산이 왜 Δ²/12인지 유도하고, 그것이 비트당 6.02dB라는 SNR로 번역되는 과정을 실측과 대조했습니다. 이상치 하나가 나머지 값의 유효 비트를 어떻게 먹는지, 그리고 int4에서 성능이 무너지는 지점을 오차 예산으로 미리 계산하는 법까지.
최댓값 빼기, 로그 공간, log1p·expm1, 분산의 두 공식, 정규화의 ε, fp32 누산, 역행렬 대신 solve — 프레임워크가 몰래 해 주는 일곱 가지를 하나씩 꺼내 각각 어떤 고장을 막는지 직접 재 봤습니다. 수식을 그대로 옮긴 코드가 왜 라이브러리보다 나쁜지에 대한 목록입니다.
0.1 + 0.2가 0.3이 아닌 이유부터 시작해 머신 엡실론을 유도하고, 같은 16비트인데 fp16과 bf16이 서로 다른 지점에서 터지는 이유, 비슷한 수를 뺄 때 유효자리가 사라지는 파괴적 상쇄, 그리고 1,000만 개를 순서만 바꿔 더했을 때 오차가 백만 배 갈리는 실험까지 직접 재 봤습니다.