파이썬 표준 라이브러리만으로 선형 회귀: 경사 하강법이 발산하는 순간 실측
numpy 없이 선형 회귀를 닫힌 형태 해, 배치 GD, 미니배치, SGD 네 가지로 구현했다. 학습률이 크면 발산하는 손실 값과 스케일링 없이 140,538 에포크가 걸린 실측, SGD 의 seed 의존성까지 담았다.
머신러닝 면접이나 스터디에서 "경사 하강법 업데이트 식을 유도해 보세요", "학습률이 너무 크면 어떻게 되나요"라는 질문은 단골입니다. 답은 외워도, 실제로 손실이 몇 에포크 만에 얼마나 커지는지 본 적은 드뭅니다. 이 글은 그 빈칸을 채우려는 개발자를 위해, numpy 없이 파이썬 표준 라이브러리만으로 선형 회귀를 네 가지 방법으로 구현하고 실제로 망가지는 모습을 숫자로 보여 줍니다. 결과부터 말하면 닫힌 형태 해는 0.47ms 만에 참값 기울기를 오차 0.01 이내로 찾았고, 경사 하강법은 스케일링을 하면 2 에포크, 안 하면 140,538 에포크가 걸렸습니다. 학습률을 한계보다 6% 크게 잡자 손실은 100 에포크 만에 5,014 에서 9.3×10¹² 로 치솟았습니다.
이 글은 시리즈 「머신러닝 밑바닥부터 (파이썬 표준 라이브러리만)」의 1강입니다. 뒤의 강의도 손실, 루프, 멈추는 조건이라는 같은 틀 위에 쌓입니다: k-평균 군집화, 나이브 베이즈 텍스트 분류기, 퍼셉트론에서 신경망까지, TF-IDF/BM25 검색 엔진.
실행 환경: Python 3.9.6 표준 라이브러리만 사용(math, random, time, sys), Apple Silicon Mac(M1 Pro). 사용 중단 예정 API 를 쓰지 않았으므로 Python 3.12 이상에서도 그대로 돌아갈 것입니다. 난수는 모두 seed 를 고정했습니다.
데이터: 정답을 아는 생성기
학습 알고리즘을 검증하려면 정답을 아는 데이터가 필요합니다. 부동산 앱에서 매물을 고를 때 흔히 보는 세 가지 항목을 흉내 내, 전용면적(20~150㎡), 연식(0~40년), 역까지 거리(0.1~3.0km)를 가진 매물 600건과 가격을 만듭니다.
가격 = 5.0 + 0.9 × 면적 − 0.4 × 연식 − 6.0 × 역 거리 + 가우스 노이즈(σ = 4)
행은 처음부터 무작위 순서이므로 앞의 480건을 학습용, 나머지 120건을 테스트용으로 씁니다. σ = 4 이므로 노이즈만으로 분산이 16 이고, 어떤 모델도 MSE 를 16 근처보다 크게 낮출 수 없습니다. 이것이 판단의 기준선입니다.
"""Linear regression from scratch: closed form vs gradient descent (Python stdlib only)."""
import math, random, sys, time
TRUE_B, TRUE_W = 5.0, [0.9, -0.4, -6.0] # the answer we hope to recover
def make_data(n=600, seed=42):
rng = random.Random(seed)
X, y = [], []
for _ in range(n):
x = [rng.uniform(20, 150), rng.uniform(0, 40), rng.uniform(0.1, 3.0)] # area, age, km
X.append(x)
y.append(TRUE_B + sum(w * v for w, v in zip(TRUE_W, x)) + rng.gauss(0, 4))
return X, y
def add_bias(X):
return [[1.0] + row for row in X]
def predict(row, w):
return sum(a * b for a, b in zip(row, w))
def mse(A, y, w):
return sum((predict(r, w) - t) ** 2 for r, t in zip(A, y)) / len(y)
add_bias 는 각 행 맨 앞에 상수 1.0 을 붙입니다. 절편을 "0번 가중치"로 다룰 수 있어서, 아래 식들에서 절편을 따로 처리할 필요가 없습니다.
닫힌 형태 해: 연립방정식을 한 번 풀기
손실을 L(w) = (1/n)·Σ(xᵢ·w − yᵢ)² 로 쓰고, 가중치 wⱼ 로 편미분합니다.
∂L/∂wⱼ = (2/n)·Σ(xᵢ·w − yᵢ)·xᵢⱼ
모든 j 를 묶으면 ∇L = (2/n)·Xᵀ(Xw − y) 입니다. 이 한 줄이 곧 경사 하강법의 전부이기도 하니 기억해 두세요. 최솟점에서는 기울기가 0 이므로 정규방정식 XᵀX w = Xᵀy 가 나옵니다. 가중치가 4개라 XᵀX 는 4×4 행렬이고, 가우스 소거법으로 풉니다. 역행렬은 만들지 않습니다. 푸는 쪽이 계산도 적고 정확도도 좋기 때문입니다.
# ---------- 1. closed form: solve (X^T X) w = X^T y ----------
def solve(M, v):
"""Gaussian elimination with partial pivoting."""
d = len(v)
M = [row[:] + [v[i]] for i, row in enumerate(M)]
for c in range(d):
p = max(range(c, d), key=lambda r: abs(M[r][c]))
M[c], M[p] = M[p], M[c]
for r in range(c + 1, d):
f = M[r][c] / M[c][c]
for k in range(c, d + 1):
M[r][k] -= f * M[c][k]
w = [0.0] * d
for r in range(d - 1, -1, -1):
w[r] = (M[r][d] - sum(M[r][k] * w[k] for k in range(r + 1, d))) / M[r][r]
return w
def closed_form(A, y):
d = len(A[0])
XtX = [[sum(r[i] * r[j] for r in A) for j in range(d)] for i in range(d)]
Xty = [sum(r[i] * t for r, t in zip(A, y)) for i in range(d)]
return solve(XtX, Xty)
출력(가중치 순서는 [절편, 면적, 연식, 역 거리]):
Python 3.9.6 train 480 test 120
true [5.000, 0.900, -0.400, -6.000]
closed [4.063, 0.909, -0.394, -6.005] train 16.6759 test 14.6532
true w train 16.8045 test 14.6847
target: train MSE <= 16.6926
세 기울기는 참값에서 0.01 이내입니다. 절편은 1 가까이 어긋났지만 버그가 아닙니다. 절편은 "면적 0㎡, 연식 0년, 역까지 0km"인 매물의 예측 가격인데, 데이터에는 20㎡ 미만 매물이 하나도 없습니다. 그만큼 바깥으로 외삽하면 기울기의 작은 오차가 크게 불어납니다. 회귀는 예측으로 평가해야지, 절편이 그럴듯한지로 평가하면 안 됩니다.
학습 MSE 16.68 은 노이즈 하한 16 바로 위입니다. 테스트 MSE 가 더 낮은(14.65) 이유는 true w 줄에 있습니다. 생성기의 참 가중치를 그대로 넣어도 학습용은 16.80, 테스트용은 14.68 이었습니다. 테스트 120건에 우연히 작은 노이즈가 붙었을 뿐입니다. 또 학습 데이터에서는 적합한 가중치가 참 가중치보다 더 좋은 값(16.68 < 16.80)을 냅니다. 최소제곱법이 노이즈를 조금 따라가기 때문이고, 테스트에서는 거의 동점(14.65 대 14.68)으로 돌아옵니다.
target 줄은 아래 모든 경사 하강법의 목표입니다. 학습 MSE 가 닫힌 형태 해의 최적값에서 0.1% 이내에 들어오면 "도달"로 봅니다.
배치 경사 하강법: 같은 미분에서 출발
경사 하강법은 방금 구한 기울기를 0 으로 놓고 푸는 대신, 비탈을 조금씩 내려갑니다.
w ← w − η·(2/m)·Xᵀ(Xw − y)
η 는 학습률, m 은 한 번 갱신에 쓰는 행 수입니다. 모든 행을 쓰면 배치 GD, 무작위로 고른 32행이면 미니배치, 한 행씩이면 확률적 경사 하강법(SGD)입니다. 그래서 batch 인자만 다른 함수 하나로 씁니다.
# ---------- 2/3. batch, mini-batch and stochastic gradient descent ----------
def gd(A, y, lr, max_epochs, batch=None, target=None, seed=0, trace=None):
rng = random.Random(seed)
n, d = len(A), len(A[0])
w, idx = [0.0] * d, list(range(n))
bs = batch or n # batch=None -> full batch
for ep in range(1, max_epochs + 1):
if batch:
rng.shuffle(idx)
for s in range(0, n, bs):
rows = idx[s:s + bs]
g = [0.0] * d
for i in rows:
e = predict(A[i], w) - y[i] # residual
for j in range(d):
g[j] += e * A[i][j]
k = 2.0 / len(rows) # d/dw of mean squared error
w = [wj - lr * k * gj for wj, gj in zip(w, g)]
try:
loss = mse(A, y, w)
except OverflowError:
loss = math.inf
if trace is not None and ep in trace:
print(f" epoch {ep:>5} loss {loss:.4g}")
if not math.isfinite(loss):
return w, ep, loss, "diverged"
if target is not None and loss <= target:
return w, ep, loss, "reached"
return w, max_epochs, loss, "not reached"
try/except OverflowError 는 실제로 프로그램이 죽어서 넣었습니다. 순수 파이썬에서는 발산한 가중치로 1e308 근처의 float 를 제곱하면 OverflowError: (34, 'Result too large') 가 나며 멈춥니다. numpy 처럼 조용히 inf 가 되지 않습니다. 에포크는 학습 데이터를 한 바퀴 도는 것이고, 에포크당 갱신 횟수는 배치 GD 1번, 미니배치 32 는 15번, SGD 는 480번입니다.
실패 1: 너무 큰 학습률
다음 절에서 설명할 표준화 특성에 풀 배치 GD 를 학습률 여섯 가지로 돌렸습니다.
[learning-rate sweep, scaled features, full batch]
lr 0.01 reached epoch 316 loss 16.69
lr 0.1 reached epoch 29 loss 16.69
lr 0.5 reached epoch 2 loss 16.69
lr 0.9 reached epoch 57 loss 16.69
lr 1.0 diverged epoch 2917 loss inf
lr 1.1 diverged epoch 1193 loss inf
눈에 띄는 점이 두 가지입니다. 크다고 빠르지 않습니다. 0.5 는 2 에포크, 0.9 는 57 에포크가 걸렸습니다. 그리고 0.9 와 1.0 사이에 절벽이 있습니다. trace 인자로 발산한 실행의 손실을 찍으면 이렇습니다.
[lr 1.1]
epoch 1 loss 7219
epoch 2 loss 1.037e+04
epoch 3 loss 1.507e+04
epoch 5 loss 3.29e+04
epoch 10 loss 2.858e+05
epoch 20 loss 5.169e+07
epoch 30 loss 1.616e+10
[lr 1.0]
epoch 1 loss 5014
epoch 10 loss 8174
epoch 100 loss 9.325e+12
epoch 500 loss 3.354e+54
epoch 1000 loss 2.954e+106
epoch 2000 loss 2.292e+210
epoch 2500 loss 2.019e+262
위 표가 숨긴 것이 있습니다. 표에서는 학습률 1.0 의 발산이 "2917 에포크"였지만, 손실은 100 에포크에 이미 9.3×10¹² 였습니다. inf 나 NaN 만 보는 검사는 수천 에포크 늦게 반응합니다. 실제 학습 코드에는 "첫 에포크 손실보다 커지면 멈춘다"도 넣어야 합니다. 이 검사였다면 두 실행 모두 10 에포크 안에 멈췄습니다.
절벽은 어디에 있을까요. MSE 손실은 그릇 모양이고 그 곡률은 헤세 행렬 H = (2/n)XᵀX 입니다. 최대 고윳값 λmax 방향에서는 스텝마다 오차가 (1 − η·λmax) 배가 됩니다. 이 값이 −1 보다 작아지면 오차는 부호를 바꾸며 커집니다. 즉 η > 2/λmax 이면 경사 하강법은 발산합니다. 거듭제곱법으로 λ 를 구했습니다.
# ---------- why a learning rate is "too high" ----------
def curvature(A):
"""largest / smallest eigenvalue of the MSE Hessian (2/n) X^T X (power iteration)"""
n, d = len(A), len(A[0])
H = [[2 / n * sum(r[i] * r[j] for r in A) for j in range(d)] for i in range(d)]
out = []
for step in (lambda v: [predict(row, v) for row in H], lambda v: solve(H, v)):
v = [1.0] * d
for _ in range(300):
u = step(v)
nrm = math.sqrt(sum(x * x for x in u))
v = [x / nrm for x in u]
out.append(predict(v, [predict(row, v) for row in H]))
return out # [lambda_max, lambda_min]
unscaled lambda_max 1.745e+04 lambda_min 0.1563 condition 111,625.8 max stable lr 2/lambda_max = 0.0001146
scaled lambda_max 2.127 lambda_min 1.878 condition 1.1 max stable lr 2/lambda_max = 0.9403
이론과 실험이 정확히 맞습니다. 0.9 < 0.9403 은 수렴했고 1.0 > 0.9403 은 발산했습니다. 0.9 가 느렸던 이유도 보입니다. 1 − 0.9 × 2.127 = −0.91 이므로, 가중치가 골짜기를 좌우로 오가며 오차를 에포크당 9% 밖에 줄이지 못합니다. 가장 좋은 고정 학습률은 2/(λmax + λmin) = 0.499 이고, 0.5 가 2 에포크 만에 끝난 이유가 이것입니다.
실패 2: 특성마다 스케일이 다를 때
위 출력의 unscaled 줄이 두 번째 실패입니다. 면적은 최대 150 이고 절편 열은 늘 1 이라서, 그릇이 한 방향으로 다른 방향보다 111,626배 가파릅니다. 학습률은 가파른 방향에 맞춰 0.000115 미만이어야 하고, 그 학습률에서는 완만한 방향이 거의 움직이지 않습니다. 한계를 조금 넘긴 1.2e-4 에서는 원래 특성도 발산합니다.
[unscaled, lr 1.2e-4]
epoch 1 loss 6047
epoch 2 loss 7201
epoch 3 loss 8584
epoch 5 loss 1.222e+04
epoch 10 loss 2.972e+04
epoch 50 loss 3.782e+07
epoch 100 loss 2.883e+11
epoch 200 loss 1.676e+19
epoch 400 loss 5.661e+34
해결책은 각 특성을 평균 0, 표준편차 1 로 맞추는 표준화입니다. 평균과 표준편차는 학습 데이터로만 계산하고 테스트 데이터에도 같은 값을 씁니다. 600건 전체로 계산하면 테스트 정보가 학습에 새어 들어갑니다(데이터 누수). 학습 뒤에는 가중치를 원래 단위로 되돌립니다.
# ---------- 4. feature scaling (standardize with TRAIN statistics only) ----------
def fit_scaler(X):
cols = list(zip(*X))
mean = [sum(c) / len(c) for c in cols]
std = [math.sqrt(sum((v - m) ** 2 for v in c) / len(c)) for c, m in zip(cols, mean)]
return mean, std
def scale(X, mean, std):
return [[(v - m) / s for v, m, s in zip(r, mean, std)] for r in X]
def unscale(w, mean, std):
"""weights learned on scaled features -> weights in original units"""
ws = [wj / s for wj, s in zip(w[1:], std)]
return [w[0] - sum(wj * m for wj, m in zip(ws, mean))] + ws
같은 목표까지의 에포크 수입니다. 스케일링 없는 쪽은 발산하지 않는 한계 바로 아래인 1.1e-4 로 돌렸습니다.
[epochs to target: scaled vs unscaled, full batch]
scaled lr 0.5 reached at epoch 2 [4.312, 0.906, -0.394, -5.980]
scaled lr 0.9 reached at epoch 57
epoch 1000 loss 34.73
epoch 10000 loss 18.55
epoch 20000 loss 17.73
epoch 100000 loss 16.74
unscaled lr 1.1e-4 reached at epoch 140538 loss 16.6926 100.7 s [3.610, 0.911, -0.389, -5.913]
전처리 하나로 2 에포크 대 140,538 에포크입니다. 스케일링 없는 실행에는 함정이 하나 더 있습니다. --long 없이 실행하면 20,000 에포크에서 끊기는데, 그때 손실 17.73 은 최적값보다 6% 높을 뿐이라 거의 수렴한 것처럼 보입니다. 하지만 가중치는 다른 이야기를 합니다.
unscaled lr 1.1e-4 not reached at epoch 20000 loss 17.7341 14.1 s [0.483, 0.925, -0.357, -5.196]
역 거리 가중치가 −6.0 이어야 하는데 −5.20, 절편은 약 4 여야 하는데 0.48 입니다. 0.1% 목표에 닿은 140,538 에포크 뒤에도 역 거리 가중치는 −5.913 으로, 닫힌 형태 해의 −6.005 와 여전히 차이가 있었습니다. 손실이 거의 줄지 않게 되어도 완만한 방향의 가중치에는 큰 오차가 남아 있을 수 있습니다. "손실이 더 안 움직이면 멈춘다"는 방식으로는 틀린 모델에서 멈출 수 있습니다.
미니배치와 SGD
같은 함수에 batch=32 와 batch=1, 표준화 특성, 셔플 seed 1 로 학습률을 몇 가지 바꿔 봤습니다.
[mini-batch and SGD, scaled, shuffle seed 1]
batch 32 lr 0.02 reached epoch 11 loss 16.6874 test 14.7088 [4.066, 0.908, -0.395, -5.974]
batch 32 lr 0.05 reached epoch 4 loss 16.6845 test 14.6414 [4.045, 0.908, -0.388, -6.005]
batch 32 lr 0.1 reached epoch 9 loss 16.6906 test 14.7560 [3.942, 0.910, -0.399, -5.992]
batch 32 lr 0.2 reached epoch 290 loss 16.6893 test 14.6619 [4.195, 0.907, -0.390, -6.058]
batch 1 lr 0.001 reached epoch 7 loss 16.6910 test 14.7085 [4.080, 0.908, -0.392, -6.007]
batch 1 lr 0.002 reached epoch 4 loss 16.6845 test 14.6064 [4.116, 0.910, -0.391, -6.045]
batch 1 lr 0.005 reached epoch 184 loss 16.6837 test 14.6878 [4.111, 0.909, -0.392, -6.053]
batch 1 lr 0.01 reached epoch 3 loss 16.6805 test 14.6572 [3.868, 0.910, -0.389, -5.972]
batch 1 lr 0.02 not reached epoch 2000 loss 17.7942 test 15.5188 [7.181, 0.882, -0.411, -6.256]
batch 1 lr 0.05 not reached epoch 2000 loss 17.8677 test 15.5023 [7.453, 0.885, -0.399, -6.600]
[SGD lr 0.005, shuffle seeds 1-5]
epochs: [184, 15, 17, 163, 68]
예상과 다른 결과였지만, 이 절에서 가장 쓸모 있는 부분이 여기입니다.
- 에포크 수가 학습률 순서대로 나오지 않습니다. SGD 에서 0.002 는 4 에포크, 0.005 는 184 에포크, 0.01 은 3 에포크였습니다. 한 행짜리 기울기는 전체 기울기의 노이즈 섞인 추정이라, 가중치가 최적점에 자리 잡지 못하고 그 주변에서 계속 흔들립니다. "목표 도달"은 흔들림 중 하나가 우연히 0.1% 폭 안에 들어왔다는 뜻입니다. 같은 0.005 도 셔플 seed 만 바꾸면 184, 15, 17, 163, 68 에포크로 들쭉날쭉했습니다. SGD 설정을 seed 하나로 비교하면 안 됩니다.
- 학습률이 너무 큰 SGD 는 폭발하지 않고 높은 곳에 머뭅니다. 0.02 와 0.05 에서는 손실이 유한한 채로 17.8 근처에 머물렀고, 절편은 7.2~7.5 까지 밀렸습니다. 흔들림 폭이 목표 폭보다 커서 끝내 들어가지 못합니다. 발산보다 조용해서 놓치기 쉬운 실패입니다.
- 미니배치는 노이즈를 평균으로 상쇄합니다. 한 번에 32행(에포크당 15번 갱신)이면 0.02~0.1 의 모든 학습률이 4~11 에포크에 도달했습니다. 0.2 에서는 흔들림이 커져 290 에포크가 걸렸습니다.
480행 정도에서는 SGD 의 이점이 없습니다. SGD 가 빛나는 것은 데이터 전체를 메모리에 올릴 수 없을 때, 새 행이 들어올 때마다 학습을 이어 가야 할 때입니다. 실무의 SGD 는 학습률을 시간에 따라 줄여 흔들림을 작게 만들지만, 여기서는 효과를 보여 주려고 고정했습니다.
실행 시간과 테스트 데이터
각각 5번 실행한 최솟값이며, 시작부터 목표 손실에 닿을 때까지의 시간입니다. SGD 줄은 학습률 0.005, seed 1, 즉 184 에포크가 걸린 운 나쁜 경우입니다.
[wall time, best of 5]
closed form 0.47 ms
batch GD lr 0.5 1.39 ms
batch GD lr 0.9 39.69 ms
mini-batch 32 3.28 ms
SGD (batch 1) 212.59 ms
test MSE closed 14.6532 batch GD 14.6333
GD 시간에는 멈출 때를 판단하려고 매 에포크 계산하는 학습 데이터 전체 손실도 들어 있습니다. 이 규모에서는 닫힌 형태 해가 압승입니다. 4×4 행렬을 한 번 만들고 작은 방정식을 한 번 풀면 끝나기 때문입니다. 테스트 MSE 차이는 0.02 로 노이즈 범위 안입니다. 학습 최적값의 0.1% 앞에서 멈춰도 손해가 없었고, 이번 테스트 데이터에서는 오히려 조금 나았습니다.
| 닫힌 형태 해 | GD (배치 / 미니배치) | SGD | |
|---|---|---|---|
| 비용 | XᵀX 만들기 O(n·d²) + 풀기 O(d³) | 에포크당 O(n·d) × 에포크 수 | 갱신당 O(d) |
| 쓸 때 | 특성 수 d 가 수천 이하, 데이터가 메모리에 들어감 | d 가 큼(d³ 가 부담), L1 처럼 닫힌 해가 없는 규제 | 데이터가 계속 들어옴, 메모리에 안 들어감, 온라인으로 계속 갱신 |
| 스케일링 | 필요 없음 | 필수(이번 2 대 140,538 에포크) | 필수 |
| 하이퍼파라미터 | 없음 | 학습률, 멈추는 조건 | 학습률(과 감쇠), 배치 크기, seed |
| 실패하는 경우 | XᵀX 가 특이(같은 특성 중복 등). 릿지 회귀나 최소제곱 솔버로 | η > 2/λmax 면 발산, 조건수가 나쁘면 느림 | 학습률이 크면 목표 위에서 계속 흔들림 |
솔직한 한계
- numpy 가 없어서 느립니다. 곱셈 하나마다 파이썬 인터프리터가 돕니다. 480행 × 가중치 4개는 괜찮지만 100만 행 × 특성 100개는 현실적이지 않습니다. 스케일링 없는 140,538 에포크는 100.7초가 걸렸는데, 같은 문제를 닫힌 형태 해는 1ms 도 안 걸려 풉니다. 벡터화된 라이브러리라면 에포크 하나가 훨씬 가벼워집니다.
- 라이브러리는 정규방정식으로 풀지 않습니다. XᵀX 를 만들면 문제의 조건수가 제곱이 됩니다. float64 에서 조건수 10⁵ 정도는 문제없지만, 거의 공선인 특성에서는 유효 숫자를 잃습니다. scikit-learn 이라면 대신
LinearRegression(XᵀX 를 만들지 않고 행렬 분해로 푸는 최소제곱 솔버), 스케일링에StandardScaler, SGD 판에SGDRegressor를 씁니다. 이 시리즈의 목적은 원리를 보는 것이라 이름만 적고 실행하지는 않았습니다. - 합성 데이터 한 종류뿐입니다. 특성은 독립이고 균등분포, 노이즈는 가우스 분포라 아주 얌전한 문제입니다. 상관된 특성은 표준화한 뒤에도 조건수를 키우고, 이상치는 제곱 오차를 크게 흔듭니다. 애초에 선형이 아닌 관계는 여기의 어떤 방법으로도 고칠 수 없습니다.
- 학습률 탐색은 작게 했습니다. 학습률 몇 개와 seed 1~5개만 시도했고, 학습률 감쇠도 쓰지 않았습니다. 실제 학습 코드에서는 스케줄을 두고, 손실이 오르면 일찍 멈춥니다.
직접 돌려 보기
위 코드 블록을 위에서부터 이어 붙이면 linreg.py 가 됩니다. 마지막에 아래 드라이버를 붙이세요. python3 linreg.py 는 이 환경에서 27초, python3 linreg.py --long 은 스케일링 없는 실행을 목표까지 돌리므로 100초쯤 더 걸립니다.
def fmt(w):
return "[" + ", ".join(f"{v:.3f}" for v in w) + "]"
def best_of(fn, k=5):
ts = []
for _ in range(k):
t = time.perf_counter(); fn(); ts.append(time.perf_counter() - t)
return min(ts)
if __name__ == "__main__":
LONG = "--long" in sys.argv # unscaled run until it converges (minutes)
X, y = make_data()
cut = int(len(X) * 0.8) # rows are already random: 480 / 120
Xtr, ytr, Xte, yte = X[:cut], y[:cut], X[cut:], y[cut:]
Atr, Ate = add_bias(Xtr), add_bias(Xte)
print(f"Python {sys.version.split()[0]} train {len(Xtr)} test {len(Xte)}")
print("true ", fmt([TRUE_B] + TRUE_W))
w_cf = closed_form(Atr, ytr)
best = mse(Atr, ytr, w_cf)
print(f"closed {fmt(w_cf)} train {best:.4f} test {mse(Ate, yte, w_cf):.4f}")
w_true = [TRUE_B] + TRUE_W # noise only: the floor for any model
print(f"true w train {mse(Atr, ytr, w_true):.4f} test {mse(Ate, yte, w_true):.4f}")
target = best * 1.001 # within 0.1% of the optimum
print(f"target: train MSE <= {target:.4f}")
mean, std = fit_scaler(Xtr)
Str, Ste = add_bias(scale(Xtr, mean, std)), add_bias(scale(Xte, mean, std))
for name, A in (("unscaled", Atr), ("scaled", Str)):
hi, lo = curvature(A)
print(f"{name:<9} lambda_max {hi:.4g} lambda_min {lo:.4g} "
f"condition {hi/lo:,.1f} max stable lr 2/lambda_max = {2/hi:.4g}")
print("\n[learning-rate sweep, scaled features, full batch]")
for lr in (0.01, 0.1, 0.5, 0.9, 1.0, 1.1):
w, ep, loss, st = gd(Str, ytr, lr, 5000, target=target)
print(f"lr {lr:<5} {st:<9} epoch {ep:>5} loss {loss:.4g}")
print("[lr 1.1]")
gd(Str, ytr, 1.1, 30, trace={1, 2, 3, 5, 10, 20, 30})
print("[lr 1.0]")
gd(Str, ytr, 1.0, 5000, trace={1, 10, 100, 500, 1000, 2000, 2500})
print("[unscaled, lr 1.2e-4]")
gd(Atr, ytr, 1.2e-4, 400, trace={1, 2, 3, 5, 10, 50, 100, 200, 400})
print("\n[epochs to target: scaled vs unscaled, full batch]")
w_s, ep, loss, st = gd(Str, ytr, 0.5, 5000, target=target)
print(f"scaled lr 0.5 {st} at epoch {ep} {fmt(unscale(w_s, mean, std))}")
w_s9, ep, loss, st = gd(Str, ytr, 0.9, 5000, target=target)
print(f"scaled lr 0.9 {st} at epoch {ep}")
t = time.perf_counter()
w_u, ep, loss, st = gd(Atr, ytr, 1.1e-4, 600000 if LONG else 20000, target=target,
trace={1000, 10000, 20000, 100000, 200000, 300000})
print(f"unscaled lr 1.1e-4 {st} at epoch {ep} loss {loss:.4f} "
f"{time.perf_counter()-t:.1f} s {fmt(w_u)}")
print("\n[mini-batch and SGD, scaled, shuffle seed 1]")
for bs, lrs in ((32, (0.02, 0.05, 0.1, 0.2)), (1, (0.001, 0.002, 0.005, 0.01, 0.02, 0.05))):
for lr in lrs:
w, ep, loss, st = gd(Str, ytr, lr, 2000, batch=bs, target=target, seed=1)
print(f"batch {bs:<3} lr {lr:<6} {st:<11} epoch {ep:>4} loss {loss:.4f} "
f"test {mse(Ste, yte, w):.4f} {fmt(unscale(w, mean, std))}")
print("[SGD lr 0.005, shuffle seeds 1-5]")
print(" epochs:", [gd(Str, ytr, 0.005, 2000, batch=1, target=target, seed=sd)[1]
for sd in range(1, 6)])
print("\n[wall time, best of 5]")
print(f"closed form {best_of(lambda: closed_form(Atr, ytr))*1000:8.2f} ms")
for lbl, lr, kw in (("batch GD lr 0.5", 0.5, {}), ("batch GD lr 0.9", 0.9, {}),
("mini-batch 32", 0.05, {"batch": 32, "seed": 1}),
("SGD (batch 1)", 0.005, {"batch": 1, "seed": 1})):
ms = best_of(lambda: gd(Str, ytr, lr, 5000, target=target, **kw)) * 1000
print(f"{lbl:<21} {ms:8.2f} ms")
print(f"\ntest MSE closed {mse(Ate, yte, w_cf):.4f} batch GD {mse(Ste, yte, w_s):.4f}")
손실, 에포크 수, 가중치는 seed 가 같으면 항상 같습니다. 실행 시간은 기기마다 다르겠지만 비율은 크게 달라지지 않을 것입니다.
다음 강의: k-평균 군집화 밑바닥부터. "기울기를 구해 한 걸음"이 아니라 "배정하고 갱신"하는 루프입니다. 같은 방식으로 실측 비교한 글로는 C 정렬 알고리즘 실측도 있습니다.
자주 묻는 질문
학습률이 크면 경사 하강법은 왜 발산하나요?
평균제곱오차 손실은 그릇 모양이고, 가장 가파른 방향의 곡률은 (2/n)XᵀX 의 최대 고윳값 λmax 입니다. 학습률이 2/λmax 를 넘으면 매 스텝 골짜기를 넘어가며 오차가 커집니다. 이번 표준화 데이터에서는 2/λmax = 0.94 였고, 0.9 는 수렴했지만 1.0 은 손실이 100 에포크 만에 5,014 에서 9.3×10¹² 까지 커졌습니다.
특성 스케일링을 하면 경사 하강법이 얼마나 빨라지나요?
같은 480행에서 표준화한 특성은 학습률 0.5 로 2 에포크 만에 목표 손실에 도달했습니다. 원래 특성은 발산하지 않는 최대 수준의 학습률로도 140,538 에포크가 필요했습니다. 방향별 곡률 차이가 111,626배였기 때문입니다.
선형 회귀에는 닫힌 형태 해와 경사 하강법 중 무엇을 써야 하나요?
특성이 적고 데이터가 메모리에 들어가면 닫힌 형태 해(최소제곱 솔버)입니다. 이번에는 0.47ms 로, 가장 잘 맞춘 경사 하강법의 1.39ms 보다 빨랐습니다. 특성이 매우 많거나, 데이터가 계속 들어오거나, 닫힌 해가 없는 규제를 붙일 때는 경사 하강법이나 SGD 를 씁니다.