Практика 15 · Регуляризація

📖 Теорія: відкрий lecture.html у цій же теці. 📝 Домашнє завдання: homework.md. 🧠 Тест: quiz.html.

Лекція показала третій спосіб приборкати перенавчання: не спрощувати модель і не добувати даних, а додати до функції втрат штраф за величину ваг. Тут ми зробимо це руками й побачимо в числах усе, про що йшлося.

Що зробимо:

1. Відтворимо перенавчений многочлен і надрукуємо його ваги — ті самі мільйони
2. Напишемо гребеневу регресію самі й переконаємось, що вона збігається з Ridge
3. Подивимось у таблиці, як штраф стискає ваги — і що жодна з них не стає нулем
4. Замінимо квадрат на модуль і порахуємо, скільки ваг Lasso обнулив точно
5. Підберемо силу штрафу крос-валідацією через RidgeCV і LassoCV
6. Перевіримо, що без зрівнювання масштабу відповідь залежить від одиниць виміру

Числа тут не збігатимуться з лекцією до гривні: там дані породжував генератор браузера, тут — NumPy. Явище й усі висновки ті самі.

0. Ті самі оголошення про телефони

Дошка оголошень про вживані телефони рідкісної марки. Ціна залежить від року випуску: телефон дешевшає з віком, але не по прямій, плюс невеликий горб на 2019-му — тодішня серія вийшла вдалою. До закономірності додається шум: подряпина на корпусі, продавець поспішає, продавець поставив із запасом.

У житті істину ніхто не знає. Тут знаємо ми, бо самі її задали, — і тому зможемо показати пальцем, де модель вивчила закономірність, а де шум.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from numpy.polynomial import legendre

ПЕРШИЙ_РІК = 2010          # найстаріший телефон на дошці
ОСТАННІЙ_РІК = 2025        # найновіший
РОЗКИД_ЦІН = 750           # шум: стандартне відхилення в гривнях


def грн(число):
    """Число з пробілами між тисячами: 1 615 407 читається краще за 1615407."""
    return f"{число:,.0f}".replace(",", " ")


def справжня_ціна(рік):
    """Закономірність, якої модель не знає: здешевлення з віком плюс горб на 2019."""
    здешевлення = 1900 + 18500 * 0.80 ** (ОСТАННІЙ_РІК - рік)
    вдала_серія = 1500 * np.exp(-((рік - 2019) / 1.3) ** 2)
    return здешевлення + вдала_серія


def згенерувати_оголошення(rng, скільки):
    """Оголошення = справжня ціна свого року плюс випадкове відхилення."""
    роки = np.sort(rng.uniform(ПЕРШИЙ_РІК, ОСТАННІЙ_РІК, скільки))
    ціни = справжня_ціна(роки) + rng.normal(0, РОЗКИД_ЦІН, скільки)
    return роки, ціни


rng = np.random.default_rng(42)
роки_навчання, ціни_навчання = згенерувати_оголошення(rng, 16)
роки_тесту, ціни_тесту = згенерувати_оголошення(rng, 300)

print(f"навчальних оголошень: {len(роки_навчання)}")
print(f"тестових оголошень:   {len(роки_тесту)}")
print("\nперші чотири навчальні оголошення:")
for рік, ціна in zip(роки_навчання[:4], ціни_навчання[:4]):
    print(f"  {рік:.2f} → {грн(ціна):>7} ₴")

навчальних оголошень: 16
тестових оголошень:   300

перші чотири навчальні оголошення:
  2011.41 →   3 069 ₴
  2011.92 →   2 180 ₴
  2013.41 →   3 952 ₴
  2015.56 →   4 116 ₴

1. Одна ознака перетворюється на пʼятнадцять

Модель у нас та сама лінійна регресія, просто ознак у неї більше однієї: замість самого лише року ми даємо їй пʼятнадцять стандартних форм кривої. Форма 1 — похила пряма, форма 2 — одна дуга, форма 3 — дві дуги, і далі чим більший номер, тим більше вигинів. Модель підбирає, скільки гривень додає кожна форма, — це і є її ваги.

Тут є технічна пастка. Якщо будувати форми як сирі степені року — рік, рік², …, рік¹⁵ — то при роках близько 2020 значення рік¹⁵ має порядок 10⁴⁹, стовпці матриці стають майже однаковими, і замість перенавчання ми отримаємо чисельне сміття. Рятує це дві дії разом: стиснути роки у відрізок [-1, 1] і взяти не сирі степені, а многочлени Лежандра.

Нульову форму — горизонтальну пряму — ми з матриці ознак прибираємо. Її роль грає вільний член моделі, а вільний член штрафувати не можна: інакше модель каратимуть просто за те, що телефони коштують десять тисяч, а не нуль.

СТЕПІНЬ = 15


def у_відрізок(роки):
    """Стискаємо роки в [-1, 1]: без цього високі степені розвалюються чисельно."""
    return 2 * (роки - ПЕРШИЙ_РІК) / (ОСТАННІЙ_РІК - ПЕРШИЙ_РІК) - 1


def ознаки(роки):
    """Пʼятнадцять форм кривої для кожного року. Нульову (сталу) відкидаємо."""
    return legendre.legvander(у_відрізок(роки), СТЕПІНЬ)[:, 1:]


X_навчання = ознаки(роки_навчання)
X_тесту = ознаки(роки_тесту)

print("матриця ознак навчання:", X_навчання.shape, "— 16 оголошень × 15 форм")
print("\nрозмах кожної форми (стандартне відхилення):")
print(np.round(X_навчання.std(axis=0), 3))
print("\nформи майже одного розмаху — це знадобиться нам у розділі 7")

матриця ознак навчання: (16, 15) — 16 оголошень × 15 форм

розмах кожної форми (стандартне відхилення):
[0.549 0.403 0.348 0.318 0.274 0.253 0.224 0.188 0.176 0.211 0.201 0.195
 0.215 0.158 0.149]

форми майже одного розмаху — це знадобиться нам у розділі 7

2. Перенавчена модель зсередини

Навчимо модель без жодного штрафу й подивимось не на криву, а на числа, з яких вона складена. Пʼятнадцять ваг плюс вільний член — рівно стільки ж вільних чисел, скільки в нас оголошень, тож крива пройде точно через кожну точку.

Дивись на порядок величин.

from sklearn.linear_model import LinearRegression


def rmse(справжні_ціни, прогнози):
    """Середня квадратична помилка в гривнях."""
    return float(np.sqrt(np.mean((справжні_ціни - прогнози) ** 2)))


без_штрафу = LinearRegression().fit(X_навчання, ціни_навчання)

print("ваги моделі без штрафу:")
for номер, вага in enumerate(без_штрафу.coef_, start=1):
    print(f"  форма {номер:>2}: {грн(вага):>18} ₴")

print()
print(f"найбільша вага:      {грн(np.abs(без_штрафу.coef_).max())} ₴")
print(f"сума модулів ваг:    {грн(np.abs(без_штрафу.coef_).sum())} ₴")
print(f"помилка на навчанні: {грн(rmse(ціни_навчання, без_штрафу.predict(X_навчання)))} ₴")
print(f"помилка на тесті:    {грн(rmse(ціни_тесту, без_штрафу.predict(X_тесту)))} ₴")

ваги моделі без штрафу:
  форма  1:       -815 392 351 ₴
  форма  2:      1 259 912 020 ₴
  форма  3:     -1 575 066 777 ₴
  форма  4:      1 735 873 393 ₴
  форма  5:     -1 747 504 061 ₴
  форма  6:      1 625 020 889 ₴
  форма  7:     -1 407 702 978 ₴
  форма  8:      1 128 366 008 ₴
  форма  9:       -835 024 028 ₴
  форма 10:        557 842 723 ₴
  форма 11:       -333 230 013 ₴
  форма 12:        169 679 488 ₴
  форма 13:        -72 351 840 ₴
  форма 14:         22 435 160 ₴
  форма 15:         -4 770 957 ₴

найбільша вага:      1 747 504 061 ₴
сума модулів ваг:    13 290 172 687 ₴
помилка на навчанні: 0 ₴
помилка на тесті:    1 065 790 375 ₴

Ваги — мільярди гривень, а прогнозує ця модель телефони по десять тисяч. Сусідні доданки протилежних знаків майже повністю гасять одне одного, і відповідь моделі — це маленький залишок від різниці величезних чисел. Такий залишок нічим не контрольований: саме тому між точками крива пірнає в мінус мільйони.

Помилка на навчанні при цьому нульова. Якби ти дивився лише на неї, ти б доповів, що задачу розвʼязано ідеально.

3. Штраф: пишемо гребеневу регресію самі

Гребенева регресія (Ridge) мінімізує не просто суму квадратів промахів, а суму промахів плюс α помножене на суму квадратів ваг. Розвʼязок знаходиться так само, як звичайний МНК, тільки до діагоналі матриці додається α:

(XᵀX + αI) · w = Xᵀy

Щоб вільний член лишився без штрафу, ознаки й відповідь спершу центрують — віднімають середні. Напишемо це в чотири рядки й порівняємо з бібліотечним Ridge.

from sklearn.linear_model import Ridge


def наш_гребінь(X, y, альфа):
    """Ridge вручну. Центрування потрібне, щоб вільний член лишився без штрафу."""
    середні_ознаки = X.mean(axis=0)
    середня_ціна = y.mean()
    Xц = X - середні_ознаки
    yц = y - середня_ціна
    ваги = np.linalg.solve(Xц.T @ Xц + альфа * np.eye(X.shape[1]), Xц.T @ yц)
    вільний_член = середня_ціна - середні_ознаки @ ваги
    return ваги, вільний_член


наші_ваги, наш_вільний = наш_гребінь(X_навчання, ціни_навчання, альфа=1.0)
бібліотечна = Ridge(alpha=1.0).fit(X_навчання, ціни_навчання)

print(f"найбільше розходження ваг:  {np.abs(наші_ваги - бібліотечна.coef_).max():.2e} ₴")
print(f"розходження вільного члена: {abs(наш_вільний - бібліотечна.intercept_):.2e} ₴")

assert np.allclose(наші_ваги, бібліотечна.coef_), "розрахунок розійшовся!"
assert np.allclose(наш_вільний, бібліотечна.intercept_), "вільний член розійшовся!"
print("\n✅ збігається — усередині Ridge той самий доданок αI на діагоналі")

найбільше розходження ваг:  9.09e-13 ₴
розходження вільного члена: 9.09e-13 ₴

✅ збігається — усередині Ridge той самий доданок αI на діагоналі

4. Що штраф робить із вагами

Тепер прогонимо α від майже нуля до дуже великого й подивимось одразу на чотири речі: наскільки стиснулась найбільша вага, скільки ваг стало точно нулем, і як рухаються обидві помилки.

рядки = []
for альфа in [0.001, 0.01, 0.1, 1, 10, 100, 1000]:
    модель = Ridge(alpha=альфа).fit(X_навчання, ціни_навчання)
    рядки.append({
        "альфа": альфа,
        "найбільша вага": round(np.abs(модель.coef_).max()),
        "сума модулів": round(np.abs(модель.coef_).sum()),
        # саме точний нуль, а не «менше за 1e-6»: у цьому вся різниця L1 і L2
        "точних нулів": int((модель.coef_ == 0).sum()),
        "помилка навчання": round(rmse(ціни_навчання, модель.predict(X_навчання))),
        "помилка тесту": round(rmse(ціни_тесту, модель.predict(X_тесту))),
    })

таблиця_ridge = pd.DataFrame(рядки).set_index("альфа")
print(таблиця_ridge.to_string())

          найбільша вага  сума модулів  точних нулів  помилка навчання  помилка тесту
альфа                                                                                
0.001               7841         46257             0               198           2894
0.010               7530         26444             0               242           1416
0.100               7346         20253             0               291           1081
1.000               6264         15235             0               897           1241
10.000              2623          6250             0              3264           3205
100.000              388           916             0              4740           4554
1000.000              41            96             0              4970           4766

Три речі, які варто прочитати в цій таблиці:

• найбільша вага падає монотонно — з тисяч до десятків. Мільярдів немає й близько;
• колонка «точних нулів» усю дорогу показує 0. Похідна від w² дорівнює 2w і зникає біля нуля: чим менша вага, тим слабше на неї тиснуть, тож останній крок до самого нуля ніхто не робить;
• помилка на тесті йде літерою U: спершу падає в тисячі разів, потім розвертається вгору. Останній рядок — це вже недонавчання: обидві помилки високі й розриву між ними немає.

5. Ласо: модуль замість квадрата

Замінимо в штрафі w² на |w|. Похідна від модуля дорівнює ±1 завжди, скільки б вага не важила, — тиск не слабшає біля нуля й доводить вагу рівно до нього.

Одна технічна деталь: scikit-learn мінімізує для Lasso вираз (1/2n)·Σ(y−ŷ)² + α·Σ|w|, а для Ridge — Σ(y−ŷ)² + α·Σw². Через цей множник однакове alpha означає для них зовсім різну силу штрафу, тому й сітки різні.

from sklearn.linear_model import Lasso

рядки = []
for альфа in [0.1, 1, 10, 30, 100, 300, 1000]:
    # ознаки-многочлени сильно корельовані, тож покоординатному спуску
    # всередині Lasso потрібно більше проходів, ніж за замовчуванням
    модель = Lasso(alpha=альфа, max_iter=2_000_000, tol=1e-6).fit(X_навчання, ціни_навчання)
    рядки.append({
        "альфа": альфа,
        "найбільша вага": round(np.abs(модель.coef_).max()),
        "точних нулів": int((модель.coef_ == 0).sum()),
        "помилка навчання": round(rmse(ціни_навчання, модель.predict(X_навчання))),
        "помилка тесту": round(rmse(ціни_тесту, модель.predict(X_тесту))),
    })

таблиця_lasso = pd.DataFrame(рядки).set_index("альфа")
print(таблиця_lasso.to_string())

        найбільша вага  точних нулів  помилка навчання  помилка тесту
альфа                                                                
0.1               9269             3               190           3676
1.0               8396             5               246           1572
10.0              7757             6               317            884
30.0              7821            10               467            844
100.0             7546            12               679            849
300.0             6986            13              1097            991
1000.0            5185            14              2548           2146

Колонка «точних нулів» тепер росте. Це не «майже нуль» і не «дуже мало» — це нуль, який дорівнює нулю, тобто ознака з моделі зникла зовсім. Подивимось поіменно, які форми ласо лишило собі.

ласо = Lasso(alpha=30, max_iter=2_000_000, tol=1e-6).fit(X_навчання, ціни_навчання)

вижили = [номер for номер, вага in enumerate(ласо.coef_, start=1) if вага != 0]
викинуті = [номер for номер, вага in enumerate(ласо.coef_, start=1) if вага == 0]

print("форми, які ласо лишило: ", вижили)
print("форми, які ласо викинуло:", викинуті)
print(f"\nмодель тепер спирається на {len(вижили)} ознак замість {СТЕПІНЬ},")
print("і вирішила це вона сама — ми лише виставили ціну за складність")

форми, які ласо лишило:  [1, 2, 3, 6, 7]
форми, які ласо викинуло: [4, 5, 8, 9, 10, 11, 12, 13, 14, 15]

модель тепер спирається на 5 ознак замість 15,
і вирішила це вона сама — ми лише виставили ціну за складність

6. Скільки штрафу насипати: крос-валідація

Досі ми підглядали в тестову помилку, щоб знайти найкраще α. У справжній задачі так робити не можна: дані, за якими ти щось обираєш, більше не годяться, щоб оцінювати обране.

Правильний інструмент — крос-валідація на навчальних оголошеннях. Оголошень у нас шістнадцять, тож беремо вісім частин по два: у кожному розбитті модель бачить чотирнадцять прикладів, а міряється на двох, яких не бачила. Тестового конверта не торкаємось узагалі.

from sklearn.linear_model import RidgeCV, LassoCV

сітка_ridge = np.logspace(-3, 3, 7)     # 0.001, 0.01, … 1000 — крок у порядок
сітка_lasso = np.logspace(-1, 3, 9)     # своя сітка: штраф міряється інакше

гребінь_cv = RidgeCV(alphas=сітка_ridge, cv=8).fit(X_навчання, ціни_навчання)
ласо_cv = LassoCV(alphas=сітка_lasso, cv=8,
                  max_iter=2_000_000, tol=1e-6).fit(X_навчання, ціни_навчання)

print(f"RidgeCV обрала α = {гребінь_cv.alpha_:g}")
print(f"   помилка на тесті: {грн(rmse(ціни_тесту, гребінь_cv.predict(X_тесту)))} ₴")
print(f"LassoCV обрала α = {ласо_cv.alpha_:g}, обнуливши "
      f"{int((ласо_cv.coef_ == 0).sum())} ваг із {СТЕПІНЬ}")
print(f"   помилка на тесті: {грн(rmse(ціни_тесту, ласо_cv.predict(X_тесту)))} ₴")
print()
print(f"а без штрафу було:   {грн(rmse(ціни_тесту, без_штрафу.predict(X_тесту)))} ₴")

RidgeCV обрала α = 0.1
   помилка на тесті: 1 081 ₴
LassoCV обрала α = 31.6228, обнуливши 10 ваг із 15
   помилка на тесті: 842 ₴

а без штрафу було:   1 065 790 375 ₴

Модель не змінилась ані на коефіцієнт: той самий степінь 15, ті самі пʼятнадцять форм. Змінилось лише те, що тепер вона платить за їхню величину.

Намалюймо всі три криві на одній картинці.

сітка_років = np.linspace(ПЕРШИЙ_РІК, ОСТАННІЙ_РІК, 600)
X_сітки = ознаки(сітка_років)

fig, ax = plt.subplots(figsize=(9, 4.6))
ax.plot(сітка_років, справжня_ціна(сітка_років), "--", color="grey",
        label="справжня залежність")
ax.plot(сітка_років, без_штрафу.predict(X_сітки), color="crimson", lw=1.6,
        label="без штрафу")
ax.plot(сітка_років, гребінь_cv.predict(X_сітки), color="teal", lw=2, label="Ridge")
ax.plot(сітка_років, ласо_cv.predict(X_сітки), color="darkorange", lw=2, label="Lasso")
ax.scatter(роки_навчання, ціни_навчання, color="black", zorder=5, s=28,
           label="16 оголошень")
ax.set_ylim(0, 24000)
ax.set_xlabel("рік випуску")
ax.set_ylabel("ціна, ₴")
ax.set_title("Той самий многочлен 15 степеня — без штрафу і зі штрафом")
ax.legend(fontsize=9)
ax.grid(alpha=0.3)

print("крива без штрафу майже не видна на графіку — вона вилітає за його межі:")
print(f"її мінімум на цьому проміжку — {грн(без_штрафу.predict(X_сітки).min())} ₴")
plt.tight_layout()
plt.show()

крива без штрафу майже не видна на графіку — вона вилітає за його межі:
її мінімум на цьому проміжку — -16 574 858 ₴

7. Без однакового масштабу штраф стає лотереєю

Штраф α·Σw² однаковий для всіх ваг: він не питає, у чому виміряна ознака. Але ваги живуть у різних масштабах саме через одиниці. Запиши ознаку в одиницях, у тисячу разів дрібніших, — вага стане в тисячу разів меншою, а її внесок у штраф впаде в мільйон разів, бо в штрафі стоїть квадрат.

Перевіримо це. Візьмемо пʼять різних наборів множників — уявних «одиниць виміру» для наших пʼятнадцяти форм — і кожного разу підберемо α крос-валідацією. Двічі: без масштабування й з StandardScaler усередині конвеєра.

from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline

еталон = rmse(ціни_тесту, гребінь_cv.predict(X_тесту))   # ознаки як є, розділ 6
rng_одиниці = np.random.default_rng(0)

рядки = []
for варіант in range(1, 6):
    # множники від тисячних до сотень — так виглядають ознаки в різних одиницях
    множники = 10.0 ** rng_одиниці.uniform(-2.5, 2.5, СТЕПІНЬ)
    X_н, X_т = X_навчання * множники, X_тесту * множники

    без_шкали = RidgeCV(alphas=сітка_ridge, cv=8).fit(X_н, ціни_навчання)
    зі_шкалою = make_pipeline(StandardScaler(),
                              RidgeCV(alphas=сітка_ridge, cv=8)).fit(X_н, ціни_навчання)

    рядки.append({
        "варіант одиниць": варіант,
        "без StandardScaler": round(rmse(ціни_тесту, без_шкали.predict(X_т))),
        "зі StandardScaler": round(rmse(ціни_тесту, зі_шкалою.predict(X_т))),
    })

таблиця_масштабу = pd.DataFrame(рядки).set_index("варіант одиниць")
print(таблиця_масштабу.to_string())
print()
print(f"для порівняння, ознаки як є (розділ 6): {round(еталон)} ₴")
print()
print("без StandardScaler:", таблиця_масштабу["без StandardScaler"].nunique(),
      "різних відповідей із 5")
print("зі StandardScaler: ", таблиця_масштабу["зі StandardScaler"].nunique(),
      "різних відповідей із 5")

assert таблиця_масштабу["зі StandardScaler"].nunique() == 1, "масштаб має зникати!"
print("\n✅ StandardScaler робить відповідь незалежною від одиниць виміру")

                 без StandardScaler  зі StandardScaler
варіант одиниць                                       
1                              2576               2151
2                              5048               2151
3                              1557               2151
4                              1749               2151
5                              2382               2151

для порівняння, ознаки як є (розділ 6): 1081 ₴

без StandardScaler: 5 різних відповідей із 5
зі StandardScaler:  1 різних відповідей із 5

✅ StandardScaler робить відповідь незалежною від одиниць виміру

Читай цю таблицю уважно, бо висновок тут тонший, ніж «масштабування покращує якість».

Ліва колонка — пʼять різних чисел. Дані ті самі, модель та сама, α щоразу підібране крос-валідацією чесно. Змінились лише уявні одиниці виміру — те, що не має жодного стосунку до задачі, — і результат поїхав у рази. Який саме ти отримаєш, залежить від того, у чому колега записав ознаку, і дізнатись, пощастило тобі чи ні, можна лише за тестовою вибіркою, якої в реальній задачі немає.

Права колонка — одне число пʼять разів. StandardScaler віднімає середнє й ділить на розкид, тому будь-які множники зникають ще до навчання. Відповідь перестала залежати від одиниць — вона тепер визначається даними, а не випадковим рішенням.

І чесна деталь наостанок: StandardScaler не зобовʼязаний давати найкраще число. Наші пʼятнадцять форм і без нього були приблизно одного розмаху (ми дивились на це в розділі 1), і варіант «як є» тут виявився трохи вдалішим. Зрівнювання — це теж вибір, і різні способи зрівняти дають різні результати. Але кожен із них відтворюваний, а варіант «не зрівнювати нічого» — ні.

---

🎯 Завдання

🟢 Рівень 1

Побудуй графік «сила штрафу → найбільша вага моделі» для Ridge з alpha від 1e-3 до 1e3 (не менше двадцяти точок, вісь alpha логарифмічна). Поруч, другою лінією, — «сила штрафу → помилка на тесті».

Зроблено, якщо: на графіку видно, що найбільша вага спадає монотонно, а помилка на тесті має мінімум, і ти назвав alpha, при якому цей мінімум досягається.

🟡 Рівень 2

Додай до пʼятнадцяти форм ще десять ознак-сміття — стовпців випадкових чисел, ніяк не повʼязаних із ціною. Підбери alpha крос-валідацією для RidgeCV і для LassoCV і порахуй, скільки ознак-сміття обнулила кожна модель.

Зроблено, якщо: Lasso обнулив принаймні половину сміття, Ridge — жодної ознаки, і ти пояснив двома реченнями, звідки береться ця різниця.

🔴 Рівень 3

Реалізуй ласо з нуля покоординатним спуском: по черзі перебирай ваги й для кожної став w = знак(ρ)·max(|ρ| − поріг, 0) / нормаΔ, де ρ — скалярний добуток цієї ознаки із залишком, порахованим без неї. Порівняй результат із sklearn.linear_model.Lasso через assert np.allclose(...).

Зроблено, якщо: твої ваги збігаються з бібліотечними для трьох різних alpha, і ти показав, що саме операція max(|ρ| − поріг, 0) дає точні нулі, а її аналог для гребеневої регресії (ділення на 1 + α) — ні.