소프트맥스 회귀 구현¶

import numpy as np

def softmax(z):
    """수치적으로 안정한 소프트맥스.

    가장 큰 값을 빼고 나서 exp 를 씌운다. 지수가 커지면 exp 가 넘쳐
    inf 가 되는데, 모든 항에서 같은 값을 빼면 분자와 분모에서 약분되어
    결과는 그대로이면서 넘침만 막을 수 있다.
    """
    z_shifted = z - np.max(z, axis=1, keepdims=True)
    exp_z = np.exp(z_shifted)
    return exp_z / np.sum(exp_z, axis=1, keepdims=True)

# 하필 3x3 인 로짓. keepdims 를 빠뜨려도 오류가 나지 않는다.
Z = np.array([[2.0, 1.0, -1.0],
              [0.0, 3.0,  1.0],
              [5.0, 5.0,  5.0]])

good = softmax(Z)
shift_wrong = np.exp(Z - np.max(Z, axis=1))
bad = shift_wrong / np.sum(shift_wrong, axis=1, keepdims=True)

print("keepdims=True :\n", np.round(good, 6))
print("keepdims 뺀 것:\n", np.round(bad, 6))
print("틀린 쪽의 행 합", bad.sum(axis=1), "  두 결과가 같은가", np.allclose(good, bad))

# 정사각이 아니면 그 자리에서 오류가 난다.
try:
    Zr = np.array([[2.0, 1.0, -1.0], [0.0, 3.0, 1.0]])
    np.exp(Zr - np.max(Zr, axis=1))
except ValueError as e:
    print("정사각이 아니면:", type(e).__name__, e)

keepdims=True :
 [[0.705385 0.259496 0.035119]
 [0.04201  0.843795 0.114195]
 [0.333333 0.333333 0.333333]]
keepdims 뺀 것:
 [[0.878878 0.118943 0.002179]
 [0.11731  0.866813 0.015876]
 [0.705385 0.259496 0.035119]]
틀린 쪽의 행 합 [1. 1. 1.]  두 결과가 같은가 False
정사각이 아니면: ValueError operands could not be broadcast together with shapes (2,3) (2,)

# 로짓의 차이가 확률의 비를 정한다. 2 와 1 의 차이가 1 이므로 첫 확률이
# 둘째의 e 배쯤 된다.
z = np.array([[2.0, 1.0, -1.0]])
print(softmax(z))
# [[0.7054  0.2595  0.0351]]

p = softmax(z)[0]
print(f"분모 e^2 + e + e^-1 = {np.exp(2) + np.exp(1) + np.exp(-1):.6f}")
print(f"p1/p2 = {p[0] / p[1]:.9f}   e   = {np.e:.9f}")
print(f"p2/p3 = {p[1] / p[2]:.9f}   e^2 = {np.e**2:.9f}")

[[0.70538451 0.25949646 0.03511903]]
분모 e^2 + e + e^-1 = 10.475217
p1/p2 = 2.718281828   e   = 2.718281828
p2/p3 = 7.389056099   e^2 = 7.389056099

def cross_entropy_loss(Y, Y_hat, eps=1e-12):
    """평균 교차엔트로피 손실.

    원-핫 이름표를 곱하므로 실제로는 "참 범주에 준 확률의 로그"만 더해진다.
    참 범주에 0 에 가까운 확률을 주면 손실이 무한대로 치솟는다.
    eps 는 로그가 발산하는 것을 막는 안전장치다.

    매개변수
    --------
    Y : (n, C) 원-핫 이름표 행렬
    Y_hat : (n, C) 예측확률 행렬
    """
    n = Y.shape[0]
    return -np.sum(Y * np.log(Y_hat + eps)) / n

Y = np.eye(3)            # 세 관측치의 참 범주가 각각 0, 1, 2
print(f"완벽      {cross_entropy_loss(Y, Y):.6e}")
print(f"균등      {cross_entropy_loss(Y, np.full((3, 3), 1 / 3)):.6f}"
      f"   log 3 = {np.log(3):.6f}")
print(f"가장 나쁨 {cross_entropy_loss(Y, 1 - Y):.6f}"
      f"   -log(1e-12) = {-np.log(1e-12):.6f}")

완벽      -1.000089e-12
균등      1.098612   log 3 = 1.098612
가장 나쁨 27.631021   -log(1e-12) = 27.631021

def compute_gradients(X, Y, Y_hat):
    """교차엔트로피의 W, b 에 대한 기울기.

    소프트맥스와 교차엔트로피를 함께 쓰면 미분이 Y_hat - Y 라는 아주
    간단한 꼴로 떨어진다. 로지스틱 회귀의 기울기와 같은 모양이며,
    일반화선형모형 전체에서 되풀이되는 구조다.
    """
    n = X.shape[0]
    error = Y_hat - Y                   # (n, C)
    dW = X.T @ error / n                # (d, C)
    db = np.mean(error, axis=0)         # (C,)
    return dW, db

def one_hot_tmp(y, C):                  # 보기 5 에서 다시 다룬다
    Y = np.zeros((y.shape[0], C))
    Y[np.arange(y.shape[0]), y] = 1.0
    return Y

rng = np.random.default_rng(1)
n, d, C = 7, 4, 3
X = rng.normal(size=(n, d))
y = rng.integers(0, C, n)
Y = one_hot_tmp(y, C)
W = rng.normal(size=(d, C)) * 0.5
b = rng.normal(size=C) * 0.5

dW, db = compute_gradients(X, Y, softmax(X @ W + b))
print("dW 의 행 합", dW.sum(axis=1))
print("db 의 합   ", db.sum())

# 중심차분으로 모든 성분을 검사한다.
eps, worst = 1e-6, 0.0
for P, G in ((W, dW), (b, db)):
    for idx in np.ndindex(P.shape):
        o = P[idx]
        P[idx] = o + eps
        lp = cross_entropy_loss(Y, softmax(X @ W + b))
        P[idx] = o - eps
        lm = cross_entropy_loss(Y, softmax(X @ W + b))
        P[idx] = o
        num = (lp - lm) / (2 * eps)
        worst = max(worst, abs(num - G[idx]) / max(abs(num), abs(G[idx]), 1e-12))
print(f"{W.size + b.size}개 성분 모두 검사,  최대 상대오차 {worst:.2e}")

dW 의 행 합 [0.00000000e+00 0.00000000e+00 0.00000000e+00 2.77555756e-17]
db 의 합    -1.3877787807814457e-17
15개 성분 모두 검사,  최대 상대오차 9.27e-08

def one_hot(y, C):
    """정수 이름표를 원-핫 행렬로 바꾼다.

    범주에 매긴 번호를 그대로 쓰면 "2가 1보다 크다" 같은 뜻이 없는 순서가
    생긴다. 원-핫은 그 순서를 지운다.
    """
    n = y.shape[0]
    Y = np.zeros((n, C))
    Y[np.arange(n), y] = 1.0
    return Y

y = np.array([0, 2, 1, 2, 2])
Y = one_hot(y, 3)
print(Y)
print("행 합", Y.sum(axis=1), " 열 합", Y.sum(axis=0),
      " 범주별 개수", np.bincount(y, minlength=3))
print("Y^T Y =\n", Y.T @ Y)

print("이름표 -1 을 넣으면:", one_hot(np.array([-1]), 3))
try:
    one_hot(np.array([3]), 3)
except IndexError as e:
    print("이름표 3 을 넣으면:", type(e).__name__)

[[1. 0. 0.]
 [0. 0. 1.]
 [0. 1. 0.]
 [0. 0. 1.]
 [0. 0. 1.]]
행 합 [1. 1. 1. 1. 1.]  열 합 [1. 1. 3.]  범주별 개수 [1 1 3]
Y^T Y =
 [[1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 3.]]
이름표 -1 을 넣으면: [[0. 0. 1.]]
이름표 3 을 넣으면: IndexError

from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

# --- 자료 준비 ---
iris = load_iris()
X, y = iris.data, iris.target
C = len(np.unique(y))

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.3, random_state=42)

# 표준화는 훈련자료로만 적합하고 시험자료에는 변환만 적용한다. 시험자료의
# 평균과 표준편차까지 보고 맞추면 정보가 새어 들어간다.
scaler = StandardScaler()
X_train = scaler.fit_transform(X_train)
X_test = scaler.transform(X_test)

Y_train = one_hot(y_train, C)
Y_test = one_hot(y_test, C)

# --- 모수 초기화 ---
d = X_train.shape[1]
np.random.seed(0)
W = np.random.randn(d, C) * 0.01
b = np.zeros(C)
row_sums_0 = W.sum(axis=1).copy()

# --- 학습 ---
# 전체 자료로 한 번에 기울기를 구하는 순수 경사하강법이다. 자료가 105건뿐이라
# 묶음으로 나눌 까닭이 없다.
lr = 0.5
epochs = 200
loss_history = []

for epoch in range(epochs):
    Z = X_train @ W + b                 # 로짓 (n, C)
    Y_hat = softmax(Z)                  # 확률 (n, C)
    loss = cross_entropy_loss(Y_train, Y_hat)
    loss_history.append(loss)
    dW, db = compute_gradients(X_train, Y_train, Y_hat)
    W -= lr * dW
    b -= lr * db

print(f"Final training loss: {loss_history[-1]:.4f}")
print(f"첫 세대 손실 {loss_history[0]:.6f},  log 3 = {np.log(3):.6f},"
      f"  차이 {loss_history[0] - np.log(3):.6f}")
print(f"손실이 매 세대 줄었는가: {bool(np.all(np.diff(loss_history) < 0))}")
print(f"W 의 행 합  처음 {np.round(row_sums_0, 8)}")
print(f"           나중 {np.round(W.sum(axis=1), 8)}")
print(f"최대 차이 {np.abs(row_sums_0 - W.sum(axis=1)).max():.2e},  b 의 합 {b.sum():.1e}")

Final training loss: 0.1326
첫 세대 손실 1.100430,  log 3 = 1.098612,  차이 0.001817
손실이 매 세대 줄었는가: True
W 의 행 합  처음 [0.03142948 0.03131173 0.00695512 0.02008916]
           나중 [0.03142948 0.03131173 0.00695512 0.02008916]
최대 차이 1.87e-15,  b 의 합 0.0e+00

from scipy import stats

# 시험자료에서의 정확도. 가장 큰 확률을 가진 범주를 고른다.
Z_test = X_test @ W + b
Y_hat_test = softmax(Z_test)
y_pred = np.argmax(Y_hat_test, axis=1)
accuracy = np.mean(y_pred == y_test)
print(f"Test accuracy: {accuracy:.4f}")

n_test = len(y_test)
z = stats.norm.ppf(0.975)
print(f"n = {n_test},  z = {z:.6f},  z^2 = {z**2:.6f}")
print(f"윌슨 아래끝 1/(1 + z^2/n) = {1 / (1 + z**2 / n_test):.6f}")
print(f"클로퍼-피어슨 아래끝 0.025^(1/n) = {0.025 ** (1 / n_test):.6f}")
for m in (10, 45, 100, 1000):
    print(f"  n = {m:4d} 을 모두 맞혔을 때의 아래끝 {1 / (1 + z**2 / m):.4f}")

Test accuracy: 1.0000
n = 45,  z = 1.959964,  z^2 = 3.841459
윌슨 아래끝 1/(1 + z^2/n) = 0.921348
클로퍼-피어슨 아래끝 0.025^(1/n) = 0.921295
  n =   10 을 모두 맞혔을 때의 아래끝 0.7225
  n =   45 을 모두 맞혔을 때의 아래끝 0.9213
  n =  100 을 모두 맞혔을 때의 아래끝 0.9630
  n = 1000 을 모두 맞혔을 때의 아래끝 0.9962

from sklearn.linear_model import LogisticRegression

# sklearn 과 맞춰 본다. 직접 구현한 것과 비슷하게 나오면 제대로 짠 것이다.
clf = LogisticRegression(solver='lbfgs', max_iter=1000)
clf.fit(X_train, y_train)
print(f"scikit-learn accuracy: {clf.score(X_test, y_test):.4f}")

print(f"우리 W 의 행 합      {np.round(W.sum(axis=1), 6)}")
print(f"sklearn 의 범주 합   {np.round(clf.coef_.sum(axis=0), 12)}")
print(f"프로베니우스 노름   우리 {np.linalg.norm(W):.4f}"
      f"   sklearn {np.linalg.norm(clf.coef_):.4f}")

for j, k in ((0, 1), (0, 2), (1, 2)):
    a = W[:, j] - W[:, k]
    c = clf.coef_[j] - clf.coef_[k]
    cos = a @ c / np.linalg.norm(a) / np.linalg.norm(c)
    print(f"  w_{j} - w_{k}:  코사인유사도 {cos:.6f}"
          f"   길이비 {np.linalg.norm(a) / np.linalg.norm(c):.4f}")

ours = np.argmax(softmax(X_test @ W + b), axis=1)
print(f"45개 검정점에서 예측이 모두 같은가: {bool(np.all(ours == clf.predict(X_test)))}")
print(f"예측확률의 최대 차이 {np.abs(softmax(X_test @ W + b) - clf.predict_proba(X_test)).max():.4f}")

scikit-learn accuracy: 1.0000
우리 W 의 행 합      [0.031429 0.031312 0.006955 0.020089]
sklearn 의 범주 합   [-0.  0. -0. -0.]
프로베니우스 노름   우리 5.1935   sklearn 4.2964
  w_0 - w_1:  코사인유사도 0.999029   길이비 1.2699
  w_0 - w_2:  코사인유사도 0.999719   길이비 1.2030
  w_1 - w_2:  코사인유사도 0.997450   길이비 1.1921
45개 검정점에서 예측이 모두 같은가: True
예측확률의 최대 차이 0.0656

lam = 0.1
np.random.seed(0)
W_reg = np.random.randn(d, C) * 0.01
b_reg = np.zeros(C)

for epoch in range(200):
    Z = X_train @ W_reg + b_reg
    Y_hat = softmax(Z)
    loss = cross_entropy_loss(Y_train, Y_hat) + 0.5 * lam * np.sum(W_reg ** 2)
    dW, db = compute_gradients(X_train, Y_train, Y_hat)
    dW += lam * W_reg       # regularization gradient
    W_reg -= lr * dW
    b_reg -= lr * db

y_pred_reg = np.argmax(softmax(X_test @ W_reg + b_reg), axis=1)
print(f"Regularized accuracy: {np.mean(y_pred_reg == y_test):.4f}")

Regularized accuracy: 0.8889

z = np.array([[1000.0, 1001.0, 999.0]])

# 순진한 소프트맥스. 최댓값을 빼지 않았다
exp_z_naive = np.exp(z)
softmax_naive = exp_z_naive / np.sum(exp_z_naive, axis=1, keepdims=True)
print("Naive:", softmax_naive)
# 출력: [[nan nan nan]] — np.exp(1001) 이 inf 가 되기 때문이다

# 안정한 소프트맥스. 최댓값을 빼고 계산한다
print("Stable:", softmax(z))
# 출력: [[0.2447  0.6652  0.0900]]

Naive: [[nan nan nan]]
Stable: [[0.24472847 0.66524096 0.09003057]]

import matplotlib.pyplot as plt

np.random.seed(0)
W_sched = np.random.randn(d, C) * 0.01
b_sched = np.zeros(C)
W_const = W_sched.copy()
b_const = b_sched.copy()

eta0 = 1.0
gamma = 0.95
total_epochs = 500
loss_sched, loss_const = [], []

for epoch in range(total_epochs):
    # 학습률을 세대에 따라 줄여 간다
    lr_t = eta0 * (gamma ** (epoch // 50))

    # --- Scheduled ---
    Z = X_train @ W_sched + b_sched
    Yh = softmax(Z)
    loss_sched.append(cross_entropy_loss(Y_train, Yh))
    dW, db = compute_gradients(X_train, Y_train, Yh)
    W_sched -= lr_t * dW
    b_sched -= lr_t * db

    # --- Constant ---
    Z = X_train @ W_const + b_const
    Yh = softmax(Z)
    loss_const.append(cross_entropy_loss(Y_train, Yh))
    dW, db = compute_gradients(X_train, Y_train, Yh)
    W_const -= eta0 * dW
    b_const -= eta0 * db

fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(loss_sched, label="Scheduled LR")
ax.plot(loss_const, label="Constant LR", linestyle="--")
ax.set_xlabel("Epoch")
ax.set_ylabel("Cross-Entropy Loss")
ax.set_title("Training Loss: Scheduled vs Constant Learning Rate")
ax.legend()
plt.tight_layout()
plt.show()