라쏘 회귀 보기¶

import numpy as np

def soft_threshold(rho, lam):
    """연성 문턱 연산자. 좌표하강의 갱신식에 쓰인다.

    |rho| 가 lam 보다 작으면 0 으로 보내고, 크면 그만큼 깎아 0 쪽으로 당긴다.
    라쏘가 계수를 정확히 0 으로 만들 수 있는 까닭이 이 평평한 구간에 있다.
    능형의 갱신식에는 이런 구간이 없어 0 이 될 수 없다.
    """
    if rho > lam:
        return rho - lam
    elif rho < -lam:
        return rho + lam
    return 0.0

def lasso_cd(X, y, lam, max_iter=1000, tol=1e-6):
    """좌표하강으로 라쏘를 푼다.

    L1 벌점은 0 에서 미분이 되지 않아 정규방정식 같은 닫힌 해가 없다.
    대신 계수를 하나씩 돌아가며 나머지를 고정한 채 최적화하면, 각 단계가
    연성 문턱 한 줄로 끝난다. 라쏘의 표준적인 푸는 법이다.

    매개변수
    --------
    X   : (n, p) 설계행렬. 미리 표준화해야 한다.
    y   : (n,)   반응벡터
    lam : L1 벌점 모수

    돌려주는 값
    ----------
    beta : (p,) 계수벡터
    """
    n, p = X.shape
    beta = np.zeros(p)
    for _ in range(max_iter):
        beta_old = beta.copy()
        for j in range(p):
            # j 번째 변수의 몫만 되살린 부분잔차. 나머지 변수의 설명은
            # 이미 빼 놓은 상태이므로, 여기에 j 만 단순회귀하는 셈이 된다.
            r_j = y - X @ beta + X[:, j] * beta[j]
            rho_j = X[:, j] @ r_j / n
            beta[j] = soft_threshold(rho_j, lam)
        if np.max(np.abs(beta - beta_old)) < tol:
            break
    return beta

# --- (1) 의 조건과 결과를 확인한다 ---
from sklearn.linear_model import Lasso

np.random.seed(42)
n, p = 150, 10
X_raw = np.random.randn(n, p)
beta_true = np.array([4.0, -3.0, 2.0] + [0.0] * 7)
y = X_raw @ beta_true + np.random.randn(n) * 2
X = (X_raw - X_raw.mean(axis=0)) / X_raw.std(axis=0)

print("||x_j||^2 / n =", np.round((X ** 2).sum(axis=0) / n, 6))
for lam in (1.0, 0.2212, 0.05):
    mine = lasso_cd(X, y, lam)
    sk = Lasso(alpha=lam, max_iter=100000, tol=1e-12).fit(X, y).coef_
    print(f"lam = {lam:<7} 최대 차이 = {np.abs(mine - sk).max():.3e}")

||x_j||^2 / n = [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
lam = 1.0     최대 차이 = 2.753e-08
lam = 0.2212  최대 차이 = 4.225e-08
lam = 0.05    최대 차이 = 6.010e-08

def lasso_path(X, y, lambdas):
    """lambda 격자 위에서 계수 경로를 구한다.

    lambda 를 크게 잡으면 모든 계수가 0 이고, 줄여 갈수록 하나씩 살아난다.
    먼저 살아나는 변수가 그만큼 중요하다는 뜻이라, 경로 그림 자체가
    변수 선택의 이야기를 담는다.
    """
    coefs = []
    for lam in lambdas:
        beta = lasso_cd(X, y, lam)
        coefs.append(beta.copy())
    return np.array(coefs)

# --- (1) 의 lambda_max 와 꺾임점을 확인한다 ---
lambdas = np.logspace(1, -2, 60)
path = lasso_path(X, y, lambdas)
n_nz = (np.abs(path) > 1e-8).sum(axis=1)

lam_max = np.abs(X.T @ (y - y.mean())).max() / n
print(f"이론 lambda_max = {lam_max:.4f}")
print(f"격자에서 계수가 모두 0 인 가장 작은 lambda = "
      f"{lambdas[n_nz == 0][-1]:.4f},  바로 다음 격자점 = {lambdas[(n_nz > 0).argmax()]:.4f}")

from sklearn.linear_model import lars_path
knots, active, _ = lars_path(X, y - y.mean(), method='lasso')
print("정확한 꺾임점 앞 넷 =", np.round(knots[:4], 4))
print("그 자리에서 들어온 변수 =", [f"x{j+1}" for j in active[:4]])
print("격자가 준 진입 lambda  =",
      [f"{lambdas[(np.abs(path[:, j]) > 1e-8).argmax()]:.4f}" for j in range(4)])

이론 lambda_max = 4.2017
격자에서 계수가 모두 0 인 가장 작은 lambda = 4.4062,  바로 다음 격자점 = 3.9194
정확한 꺾임점 앞 넷 = [4.2017 3.3635 2.1672 0.1892]
그 자리에서 들어온 변수 = ['x1', 'x2', 'x3', 'x4']
격자가 준 진입 lambda  = ['3.9194', '3.1012', '1.9415', '0.1867']

def cv_lasso(X, y, lambdas, folds=5):
    """lambda 마다 k겹 교차검증 MSE 를 구한다.

    벌점이 강할수록 훈련오차는 반드시 커진다. 그러니 lambda 는 훈련자료가
    아니라 떼어 놓은 자료에서 재야 고를 수 있다.
    """
    n = len(y)
    indices = np.arange(n)
    np.random.shuffle(indices)
    fold_size = n // folds
    cv_mse = np.zeros(len(lambdas))

    for k in range(folds):
        val_idx = indices[k * fold_size:(k + 1) * fold_size]
        train_idx = np.setdiff1d(indices, val_idx)
        X_tr, y_tr = X[train_idx], y[train_idx]
        X_va, y_va = X[val_idx], y[val_idx]
        for i, lam in enumerate(lambdas):
            beta = lasso_cd(X_tr, y_tr, lam)
            pred = X_va @ beta
            cv_mse[i] += np.mean((y_va - pred) ** 2)
    return cv_mse / folds

# --- (1) 의 단조성을 확인한다 ---
np.random.seed(42)
grid = np.logspace(1, -2, 30)
train_mse = np.array([np.mean((y - X @ lasso_cd(X, y, l)) ** 2) for l in grid])
print(f"훈련 MSE 가 lambda 에 대해 단조증가인가? "
      f"{bool(np.all(np.diff(train_mse[::-1]) >= -1e-12))}")
print("훈련 MSE (lambda = 10, 1, 0.2212, 0.01):",
      [f"{np.mean((y - X @ lasso_cd(X, y, l)) ** 2):.3f}" for l in (10, 1, 0.2212, 0.01)])

cv = cv_lasso(X, y, grid)
print(f"CV MSE  최솟값 {cv.min():.3f} at lambda = {grid[cv.argmin()]:.4f}")
print(f"CV MSE  양 끝: lambda = {grid[0]:.2f} -> {cv[0]:.3f},  "
      f"lambda = {grid[-1]:.2f} -> {cv[-1]:.3f}")

훈련 MSE 가 lambda 에 대해 단조증가인가? True
훈련 MSE (lambda = 10, 1, 0.2212, 0.01): ['35.154', '6.848', '4.398', '4.157']
CV MSE  최솟값 4.622 at lambda = 0.1743
CV MSE  양 끝: lambda = 10.00 -> 35.154,  lambda = 0.01 -> 4.964

import matplotlib.pyplot as plt

np.random.seed(42)

# 열 변수 중 참으로 쓰이는 것은 앞의 셋뿐이다. 교차검증으로 고른 라쏘가
# 그 셋을 찾아내는지 보는 것이 목적이다.
n, p = 150, 10
X_raw = np.random.randn(n, p)
beta_true = np.array([4.0, -3.0, 2.0, 0, 0, 0, 0, 0, 0, 0])
y = X_raw @ beta_true + np.random.randn(n) * 2

# L1 벌점은 계수의 크기를 그대로 재므로, 변수의 단위가 다르면 벌점이
# 불공평하게 걸린다. 표준화는 선택이 아니라 필수다.
X_mean = X_raw.mean(axis=0)
X_std = X_raw.std(axis=0)
X = (X_raw - X_mean) / X_std

# lambda 격자는 로그 눈금으로 잡는다. 벌점이 곱셈으로 작동하기 때문이다.
lambdas = np.logspace(1, -2, 60)
path = lasso_path(X, y, lambdas)

# 교차검증으로 lambda 고르기
lambdas_cv = np.logspace(1, -2, 30)
mse_cv = cv_lasso(X, y, lambdas_cv)
best_idx = int(np.argmin(mse_cv))
best_lam = lambdas_cv[best_idx]

beta_best = lasso_cd(X, y, best_lam)
n_nonzero = np.sum(np.abs(beta_best) > 1e-8)

print(f"Best lambda (5-fold CV):  {best_lam:.4f}")
print(f"Non-zero coefficients:    {n_nonzero}  (true: 3)")
print(f"Min CV MSE:               {mse_cv[best_idx]:.3f}")

# --- 직교 근사가 예측하는 축소량과 견준다 ---
from sklearn.linear_model import LinearRegression
ols = LinearRegression().fit(X, y).coef_
pred = np.sign(ols) * np.maximum(np.abs(ols) - best_lam, 0)
print(f"X'X/n 이 단위행렬에서 벗어난 최대 크기 = {np.abs(X.T @ X / n - np.eye(p)).max():.4f}")
print(f"  최소제곱          {np.round(ols[:3], 4)}")
print(f"  연성문턱 예측      {np.round(pred[:3], 4)}")
print(f"  라쏘 실제          {np.round(beta_best[:3], 4)}")
print(f"잡음 바닥 sigma^2 = 4,  최소 CV MSE = {mse_cv[best_idx]:.3f}  "
      f"(차이 {mse_cv[best_idx] - 4:.3f})")

Best lambda (5-fold CV):  0.2212
Non-zero coefficients:    3  (true: 3)
Min CV MSE:               4.774
X'X/n 이 단위행렬에서 벗어난 최대 크기 = 0.1857
  최소제곱          [ 3.7252 -3.129   1.8683]
  연성문턱 예측      [ 3.504  -2.9078  1.6471]
  라쏘 실제          [ 3.5463 -2.9119  1.6434]
잡음 바닥 sigma^2 = 4,  최소 CV MSE = 4.774  (차이 0.774)

def lasso_cd_warm(X, y, lam, beta_init=None, max_iter=1000, tol=1e-6):
    n, p = X.shape
    beta = beta_init.copy() if beta_init is not None else np.zeros(p)
    for _ in range(max_iter):
        beta_old = beta.copy()
        for j in range(p):
            r_j = y - X @ beta + X[:, j] * beta[j]
            rho_j = X[:, j] @ r_j / n
            beta[j] = soft_threshold(rho_j, lam)
        if np.max(np.abs(beta - beta_old)) < tol:
            break
    return beta

import numpy as np

np.random.seed(0)
n, p = 100, 200
X = np.random.randn(n, p)
X = (X - X.mean(0)) / X.std(0)
beta_true = np.zeros(p)
beta_true[:5] = [3, -2, 4, -1, 2]
y = X @ beta_true + np.random.randn(n) * 1.5

lambdas_cv = np.logspace(1, -2, 30)
mse_cv = cv_lasso(X, y, lambdas_cv)
best_lam = lambdas_cv[np.argmin(mse_cv)]
beta_hat = lasso_cd(X, y, best_lam)

selected = np.abs(beta_hat) > 1e-8
true_support = np.abs(beta_true) > 0

tp = np.sum(selected & true_support)
fp = np.sum(selected & ~true_support)
print(f"True positives:  {tp}/5")
print(f"False positives: {fp}")

True positives:  5/5
False positives: 11