Практика 37 · Жорсткість і сплайн¶
📖 Теорія: відкрий
lecture.htmlу цій же теці. 📝 Домашнє завдання:homework.md· 🧪 Тест:quiz.html⏮ Перша частина теми:../36-double-descent/
У першій частині ми побачили подвійний спуск. Тут ми його порахуємо — на моделі, де все можна перевірити руками: шість точок і кусково-лінійна функція з випадковими зламами.
Що зробимо:
- Зберемо кусково-лінійну модель на базисі ½·|x − t| — рівно те, що вміє ReLU-мережа
- Реалізуємо розвʼязок мінімальної норми, що не штрафує пряму α + βx
- Покажемо, як один невдалий вузол розганяє розмах розвʼязку в сотні разів
- Порахуємо розподіл помилки по сотнях розкладів вузлів і побачимо важкий хвіст біля порогу
- Переконаємось, що норма спадає як 1/√K
- Порівняємо розвʼязок мінімальної норми з
scipy.interpolate.CubicSpline(bc_type="natural")
1. Дані: шість точок параболи, без жодного шуму¶
Дані навмисно крихітні. Відсутність шуму принципова: ми хочемо показати, що пік на порозі виникає навіть на ідеально чистих даних — тобто це властивість самої моделі, а не забруднення міток.
In [1]:
import numpy as np
import matplotlib.pyplot as plt
X_ДАНІ = np.array([0.0, 1.0, 2.0, 3.0, 4.0, 5.0])
КІЛЬКІСТЬ_ТОЧОК = len(X_ДАНІ) # це і є N
ЛІВА_МЕЖА, ПРАВА_МЕЖА = 0.0, 5.0
def істинна_функція(x):
"""Парабола з вершиною в 2.5. Саме її модель має відновити між точками."""
return 4.0 * (np.asarray(x, dtype=float) - 2.5) ** 2
Y_ДАНІ = істинна_функція(X_ДАНІ)
print(f"N = {КІЛЬКІСТЬ_ТОЧОК} точок")
print(f"x = {X_ДАНІ}")
print(f"y = {Y_ДАНІ}")
print("\nШуму немає взагалі: y — це точні значення параболи.")
N = 6 точок x = [0. 1. 2. 3. 4. 5.] y = [25. 9. 1. 1. 9. 25.] Шуму немає взагалі: y — це точні значення параболи.
2. Модель: ламана зі зламами у випадкових вузлах¶
Модель — кусково-лінійна функція, у якої злами дозволені лише в заздалегідь кинутих точках t₁ … t_K (їх називають вузлами, knots):
f(x) = α + β·x + γ₁·φ₁(x) + … + γ_K·φ_K(x), де φ_k(x) = ½·|x − t_k|
Чому саме модуль? Бо |x − t| лінійна всюди, крім самої точки t, де вона зламується. Сума таких доданків плюс пряма α + βx — це ламана з дозволеними зламами рівно у вузлах.
Вузли не навчаються. Вони кидаються випадково й далі стоять нерухомо. Навчаються тільки коефіцієнти α, β, γ. Тому задача лишається звичайною лінійною регресією — просто в нестандартному базисі.
Параметрів у моделі p = K + 2: два на пряму й по одному на кожен вузол. Отже поріг інтерполяції p = N досягається при K = N − 2 = 4 вузли.
In [2]:
def базисна_функція(x, вузол):
"""½·|x − t| — та сама «галочка», яка ламає нахил рівно на одиницю."""
return 0.5 * np.abs(np.asarray(x, dtype=float) - вузол)
def матриця_плану(x, вузли):
"""Стовпці: [1, x, φ₁(x), …, φ_K(x)]. Розмір: len(x) × (K + 2)."""
x = np.asarray(x, dtype=float)
стовпці = [np.ones_like(x), x]
for вузол in вузли:
стовпці.append(базисна_функція(x, вузол))
return np.column_stack(стовпці)
# як виглядає одна базисна функція і одна ламана з неї
сітка = np.linspace(ЛІВА_МЕЖА, ПРАВА_МЕЖА, 600)
fig, (ліва, права) = plt.subplots(1, 2, figsize=(12, 3.8))
for вузол in [1.5, 3.0]:
ліва.plot(сітка, базисна_функція(сітка, вузол), lw=2, label=f"φ(x) при t = {вузол}")
ліва.set_title("Базисна функція ½·|x − t|"); ліва.legend(); ліва.grid(alpha=.25)
приклад_вузлів = [1.2, 2.4, 3.9]
приклад_коефіцієнтів = np.array([2.0, -1.0, 3.0, -5.0, 4.0]) # α, β, γ₁, γ₂, γ₃
права.plot(сітка, матриця_плану(сітка, приклад_вузлів) @ приклад_коефіцієнтів, lw=2.2, color="teal")
for вузол in приклад_вузлів:
права.axvline(вузол, color="crimson", ls=":", lw=1.4)
права.set_title("Ламана: злами рівно у вузлах (пунктир)"); права.grid(alpha=.25)
plt.tight_layout(); plt.show()
print("Кусково-лінійні функції — це рівно те, що вміє будувати мережа з ReLU.")
print("Тому висновки цієї практики переносяться далі, ніж здається.")
Кусково-лінійні функції — це рівно те, що вміє будувати мережа з ReLU. Тому висновки цієї практики переносяться далі, ніж здається.
3. Правило підгонки¶
Так само, як у першій частині, але з однією важливою деталлю:
- p < N — звичайний метод найменших квадратів.
- p ≥ N — інтерполянтів багато; беремо той, у якого мінімальна ‖γ‖². Коефіцієнти α і β у норму не входять: пряма — «безкоштовна» частина моделі, штрафувати її немає сенсу.
Ця деталь заважає просто написати pinv: він мінімізував би довжину всього вектора,
включно з α і β. Тому спершу приберемо пряму з рівнянь.
Як прибрати. Помножимо систему Aθ = y на матрицю, рядки якої ортогональні
до стовпців «1» і «x». Тоді доданки α·1 і β·x зникають, і залишається система
тільки на γ — розміром (N − 2) × K. Її вже можна розвʼязувати через pinv,
і мінімальна норма буде саме тією, що потрібно.
In [3]:
# ортонормований базис доповнення до площини span{1, x} у просторі R^N
пряма_на_точках = np.column_stack([np.ones(КІЛЬКІСТЬ_ТОЧОК), X_ДАНІ])
ліві_вектори, _, _ = np.linalg.svd(пряма_на_точках, full_matrices=True)
ДОПОВНЕННЯ = ліві_вектори[:, 2:] # N × (N − 2) = 6 × 4
print(f"матриця доповнення: {ДОПОВНЕННЯ.shape}")
print(f"перевірка ортогональності до «1» і до «x»: "
f"{np.max(np.abs(ДОПОВНЕННЯ.T @ пряма_на_точках)):.2e}")
def навчити(вузли):
"""Коефіцієнти (α, β, γ₁ … γ_K) за правилом «МНК до порогу, мінімальна ‖γ‖ після»."""
A = матриця_плану(X_ДАНІ, вузли)
кількість_параметрів = len(вузли) + 2
if кількість_параметрів < КІЛЬКІСТЬ_ТОЧОК:
коефіцієнти, *_ = np.linalg.lstsq(A, Y_ДАНІ, rcond=None)
return коефіцієнти
# 1) прибираємо α і β з рівнянь — лишається система тільки на γ
система_на_гамма = ДОПОВНЕННЯ.T @ A[:, 2:]
права_частина = ДОПОВНЕННЯ.T @ Y_ДАНІ
гамма = np.linalg.pinv(система_на_гамма) @ права_частина # мінімальна ‖γ‖
# 2) знаючи γ, добираємо α і β так, щоб крива пройшла крізь усі точки
залишок = Y_ДАНІ - A[:, 2:] @ гамма
альфа_бета, *_ = np.linalg.lstsq(пряма_на_точках, залишок, rcond=None)
return np.concatenate([альфа_бета, гамма])
def прогноз(коефіцієнти, вузли, x):
return матриця_плану(x, вузли) @ коефіцієнти
def норма_гамма(коефіцієнти):
"""Довжина вектора коефіцієнтів при вузлах — без α і β."""
return np.linalg.norm(коефіцієнти[2:])
def сума_залишків(коефіцієнти, вузли):
"""Наскільки крива промахнулась повз шість точок даних."""
return np.sum(np.abs(Y_ДАНІ - прогноз(коефіцієнти, вузли, X_ДАНІ)))
def похибка_проти_істини(коефіцієнти, вузли, кроків=401):
"""Середнє відхилення від параболи на всьому відрізку, а не лише в шести точках."""
сітка_оцінки = np.linspace(ЛІВА_МЕЖА, ПРАВА_МЕЖА, кроків)
відхилення = прогноз(коефіцієнти, вузли, сітка_оцінки) - істинна_функція(сітка_оцінки)
return np.mean(np.abs(відхилення))
матриця доповнення: (6, 4) перевірка ортогональності до «1» і до «x»: 9.99e-16
Перевірка: чи справді від K = 4 модель інтерполює¶
Порогу відповідає K = 4. Починаючи з нього сума залишків має стати нулем — але лише якщо вузли розкидані по відрізку. Чому так, зрозуміти легко: між двома сусідніми вузлами функція строго пряма, а пряма не пройде крізь три точки параболи, які не лежать на одній лінії. Тому в кожній «дірці» шириною у дві сусідні точки має бути хоч один злам.
In [4]:
# розкидані вузли: по одному в кожному проміжку між точками даних
ПРИКЛАДИ_ВУЗЛІВ = {
0: np.array([]),
1: np.array([2.5]),
2: np.array([1.5, 3.5]),
3: np.array([1.5, 2.5, 3.5]),
4: np.array([0.5, 1.5, 2.5, 3.5]),
6: np.array([0.4, 1.2, 1.8, 2.6, 3.4, 4.2]),
12: np.linspace(0.2, 4.8, 12),
40: np.linspace(0.05, 4.95, 40),
}
print(f"{'K':>4} {'p':>4} {'сума залишків':>16} {'‖γ‖':>10} {'похибка проти істини':>22}")
for K, вузли in ПРИКЛАДИ_ВУЗЛІВ.items():
коеф = навчити(вузли)
мітка = " ← ПОРІГ" if K == КІЛЬКІСТЬ_ТОЧОК - 2 else ""
print(f"{K:>4} {K + 2:>4} {сума_залишків(коеф, вузли):16.2e} {норма_гамма(коеф):10.2f} "
f"{похибка_проти_істини(коеф, вузли):22.3f}{мітка}")
assert сума_залишків(навчити(ПРИКЛАДИ_ВУЗЛІВ[4]), ПРИКЛАДИ_ВУЗЛІВ[4]) < 1e-8
assert сума_залишків(навчити(ПРИКЛАДИ_ВУЗЛІВ[12]), ПРИКЛАДИ_ВУЗЛІВ[12]) < 1e-8
print("\n✅ від K = 4 модель проходить крізь усі шість точок точно")
K p сума залишків ‖γ‖ похибка проти істини 0 2 5.33e+01 0.00 7.308 1 3 1.07e+01 24.00 1.988 2 4 5.33e-14 22.63 0.931 3 5 2.53e-14 22.63 0.931 4 6 3.55e-14 22.63 0.931 ← ПОРІГ 6 8 3.73e-14 14.76 0.251 12 14 2.53e-14 10.80 0.153 40 42 5.62e-14 5.97 0.133 ✅ від K = 4 модель проходить крізь усі шість точок точно
А якщо вузли збилися в купу¶
Тепер зворотний бік. Кинемо чотири вузли так, щоб усі вони опинились праворуч від x = 2. Тоді на відрізку [0, 2] функція — одна пряма, а точки (0, 25), (1, 9), (2, 1) на одній прямій не лежать. Інтерполянта просто не існує, і модель змушена промахнутись.
In [5]:
ЗБИТІ_ВУЗЛИ = np.array([2.7, 3.0, 3.6, 4.7])
коеф_збиті = навчити(ЗБИТІ_ВУЗЛИ)
print(f"вузли зліплені праворуч: {ЗБИТІ_ВУЗЛИ}")
print(f"сума залишків = {сума_залишків(коеф_збиті, ЗБИТІ_ВУЗЛИ):.3f} — не нуль!")
print(f"розкидані вузли: {ПРИКЛАДИ_ВУЗЛІВ[4]}")
print(f"сума залишків = {сума_залишків(навчити(ПРИКЛАДИ_ВУЗЛІВ[4]), ПРИКЛАДИ_ВУЗЛІВ[4]):.2e}")
# наскільки часто випадковий розклад узагалі не дає інтерполянта
rng = np.random.default_rng(0)
print(f"\n{'K':>4} {'розкладів без інтерполянта':>28}")
for K in [4, 5, 6, 8, 12, 20]:
невдач = 0
for _ in range(400):
вузли = np.sort(rng.uniform(ЛІВА_МЕЖА, ПРАВА_МЕЖА, K))
if сума_залишків(навчити(вузли), вузли) > 1e-6:
невдач += 1
print(f"{K:>4} {100 * невдач / 400:27.1f}%")
assert сума_залишків(коеф_збиті, ЗБИТІ_ВУЗЛИ) > 1.0
print("\n✅ на самому порозі більш ніж половина випадкових розкладів не дає інтерполянта")
print("Це та сама жорсткість, тільки в найгострішій формі: свободи настільки мало,")
print("що іноді її не вистачає навіть на те, щоб пройти крізь дані.")
вузли зліплені праворуч: [2.7 3. 3.6 4.7] сума залишків = 5.333 — не нуль! розкидані вузли: [0.5 1.5 2.5 3.5] сума залишків = 3.55e-14 K розкладів без інтерполянта 4 58.0%
5 35.0% 6 20.8% 8 7.5%
12 0.8% 20 0.0% ✅ на самому порозі більш ніж половина випадкових розкладів не дає інтерполянта Це та сама жорсткість, тільки в найгострішій формі: свободи настільки мало, що іноді її не вистачає навіть на те, щоб пройти крізь дані.
4. Класичний режим: від прямої до порогу¶
При K = 0 модель — просто пряма, і на симетричній параболі вона вироджується в горизонталь. Додаємо вузли по одному й дивимось, як ситуація виправляється.
In [6]:
fig, axes = plt.subplots(1, 5, figsize=(17, 3.6))
rng = np.random.default_rng(11)
for вісь, K in zip(axes, [0, 1, 2, 3, 4]):
вузли = np.sort(rng.uniform(ЛІВА_МЕЖА, ПРАВА_МЕЖА, K))
коеф = навчити(вузли)
вісь.plot(сітка, істинна_функція(сітка), color="gray", ls="--", lw=1.8)
вісь.plot(сітка, прогноз(коеф, вузли, сітка), color="teal", lw=2.2)
вісь.scatter(X_ДАНІ, Y_ДАНІ, s=28, color="crimson", zorder=3)
for вузол in вузли:
вісь.axvline(вузол, color="darkorange", ls=":", lw=1.2)
заголовок = f"K = {K}, p = {K + 2}"
if K == 4:
заголовок += " ← поріг"
вісь.set_title(f"{заголовок}\nпохибка = {похибка_проти_істини(коеф, вузли):.2f}", fontsize=10)
вісь.set_ylim(-6, 30); вісь.grid(alpha=.2)
plt.tight_layout(); plt.show()
print("Поки вузлів мало, похибка чесно спадає: додавання нелінійності допомагає,")
print("як і обіцяє класична теорія. І з цим конкретним розкладом вузлів модель")
print("на порозі виглядає цілком пристойно. Але вузли — випадкові, і ми")
print("подивились лише на один щасливий жереб.")
Поки вузлів мало, похибка чесно спадає: додавання нелінійності допомагає, як і обіцяє класична теорія. І з цим конкретним розкладом вузлів модель на порозі виглядає цілком пристойно. Але вузли — випадкові, і ми подивились лише на один щасливий жереб.
5. Пастка порогу: один вузол вирішує все¶
Ось у чому суть. При K = 4 кількість рівнянь дорівнює кількості невідомих, і вибору немає взагалі. Прослідкуймо ланцюжок.
Якщо між x = 0 і x = 1 вузла не випало, то на цьому проміжку функція обовʼязково пряма — а пряма крізь (0, 25) і (1, 9) єдина, її нахил уже зафіксовано: −16. Ця пряма тягнеться до першого вузла, там ламається рівно настільки, щоб потрапити в наступну точку даних. І так далі — свободи немає на жодному кроці.
Тепер посунемо найлівіший вузол ближче до x = 2 і подивимось на розмах кривої.
In [7]:
БАЗОВІ_ВУЗЛИ = np.array([0.5, 2.6, 3.4, 4.4])
print(f"{'t₁':>6} {'розмах |f|':>14} {'‖γ‖':>12} {'сума залишків':>16}")
розмахи = []
позиції = [0.5, 1.0, 1.4, 1.7, 1.9, 1.95, 1.99]
for t1 in позиції:
вузли = БАЗОВІ_ВУЗЛИ.copy()
вузли[0] = t1
коеф = навчити(вузли)
крива = прогноз(коеф, вузли, сітка)
розмах = крива.max() - крива.min()
розмахи.append(розмах)
print(f"{t1:>6.2f} {розмах:14.1f} {норма_гамма(коеф):12.1f} {сума_залишків(коеф, вузли):16.2e}")
assert розмахи[-1] > 20 * розмахи[0], "розмах мав злетіти щонайменше в 20 разів!"
print("\n✅ розмах розвʼязку злетів у "
f"{розмахи[-1] / розмахи[0]:.0f} разів — і це при нульових залишках весь час")
print("Модель ідеально інтерполює дані й одночасно поводиться жахливо.")
t₁ розмах |f| ‖γ‖ сума залишків 0.50 28.8 31.9 8.88e-14 1.00 28.8 28.7 5.86e-14 1.40 25.6 18.6 3.82e-14 1.70 30.4 56.6 6.00e-14 1.90 77.0 269.2 4.48e-13 1.95 171.6 591.7 7.53e-13 1.99 934.5 3173.6 3.52e-12 ✅ розмах розвʼязку злетів у 32 разів — і це при нульових залишках весь час Модель ідеально інтерполює дані й одночасно поводиться жахливо.
In [8]:
fig, (ліва, права) = plt.subplots(1, 2, figsize=(13, 4.2))
for t1, колір in zip([0.5, 1.7, 1.95], ["teal", "darkorange", "crimson"]):
вузли = БАЗОВІ_ВУЗЛИ.copy()
вузли[0] = t1
коеф = навчити(вузли)
ліва.plot(сітка, прогноз(коеф, вузли, сітка), lw=2, color=колір, label=f"t₁ = {t1}")
ліва.plot(сітка, істинна_функція(сітка), color="gray", ls="--", lw=1.8, label="істина")
ліва.scatter(X_ДАНІ, Y_ДАНІ, s=28, color="black", zorder=3)
ліва.set_ylim(-80, 80)
ліва.set_title("Три однаково «ідеальні» інтерполянти"); ліва.legend(); ліва.grid(alpha=.25)
права.plot(позиції, розмахи, marker="o", lw=2.2, color="crimson")
права.set_yscale("log")
права.set_xlabel("положення лівого вузла t₁"); права.set_ylabel("розмах |f| (лог. шкала)")
права.set_title("Що ближче вузол до x = 2, то дикіша крива"); права.grid(alpha=.25, which="both")
plt.tight_layout(); plt.show()
print("Перший сегмент — це пряма крізь (0, 25) і (1, 9), тобто y = 25 − 16x.")
print("У вузлі t₁ вона вже впала до 25 − 16·t₁ — при t₁ біля двійки це майже −7.")
print("Щоб звідти дістатись до точки (2, 1), наступному сегменту потрібен нахил")
print("близько 8 / (2 − t₁) — і коли t₁ → 2, він летить у нескінченність.")
Перший сегмент — це пряма крізь (0, 25) і (1, 9), тобто y = 25 − 16x. У вузлі t₁ вона вже впала до 25 − 16·t₁ — при t₁ біля двійки це майже −7. Щоб звідти дістатись до точки (2, 1), наступному сегменту потрібен нахил близько 8 / (2 − t₁) — і коли t₁ → 2, він летить у нескінченність.
6. Наскільки погано буває насправді¶
Одна невдала конфігурація ще нічого не доводить. Порахуймо чесно: кинемо вузли 300 разів для кожного K і подивимось на розподіл помилки.
Середнє арифметичне тут рахувати немає сенсу — теорія каже, що на порозі воно нескінченне, а на практиці воно щоразу визначається одним найгіршим розкладом і не стабілізується, скільки б спроб ми не робили. Тому дивимось на стійкі величини: медіану (типовий випадок), 90-й процентиль (як буває в одному випадку з десяти) і частку катастроф.
In [9]:
КІЛЬКІСТЬ_РОЗКЛАДІВ = 300
КАТАСТРОФА = 10.0 # похибка більша за цю вважається провалом
СІТКА_K = [0, 1, 2, 3, 4, 5, 6, 8, 12, 20, 40, 80]
rng = np.random.default_rng(0)
статистика = {}
for K in СІТКА_K:
похибки, норми = [], []
for _ in range(КІЛЬКІСТЬ_РОЗКЛАДІВ):
вузли = np.sort(rng.uniform(ЛІВА_МЕЖА, ПРАВА_МЕЖА, K))
коеф = навчити(вузли)
похибки.append(похибка_проти_істини(коеф, вузли))
норми.append(норма_гамма(коеф))
похибки = np.array(похибки)
статистика[K] = {
"медіана": np.median(похибки),
"p90": np.percentile(похибки, 90),
"найгірша": похибки.max(),
"катастроф": np.mean(похибки > КАТАСТРОФА),
"середнє": похибки.mean(),
"норма": np.median(норми),
}
print(f"{'K':>4} {'p':>4} {'медіана':>10} {'90-й проц.':>12} {'найгірша':>12} "
f"{'середнє':>12} {'катастроф':>11}")
for K in СІТКА_K:
s = статистика[K]
мітка = " ← ПОРІГ" if K == 4 else ""
print(f"{K:>4} {K + 2:>4} {s['медіана']:10.3f} {s['p90']:12.3f} {s['найгірша']:12.3g} "
f"{s['середнє']:12.3g} {100 * s['катастроф']:10.1f}%{мітка}")
K p медіана 90-й проц. найгірша середнє катастроф 0 2 7.308 7.308 7.31 7.31 0.0% 1 3 4.136 6.786 7.32 4.27 0.0% 2 4 2.118 5.215 18 2.76 1.3% 3 5 1.703 5.168 2.94e+03 13.3 5.0% 4 6 1.384 5.001 1.69e+13 5.65e+10 6.3% ← ПОРІГ 5 7 1.174 5.203 1.32e+04 52.4 6.3% 6 8 0.952 3.148 110 2.47 3.3% 8 10 0.699 1.783 101 1.4 1.3% 12 14 0.413 0.749 18.5 0.542 0.3% 20 22 0.302 0.441 1.05 0.31 0.0% 40 42 0.210 0.299 0.469 0.219 0.0% 80 82 0.167 0.231 0.409 0.177 0.0%
In [10]:
медіани = np.array([статистика[K]["медіана"] for K in СІТКА_K])
процентилі = np.array([статистика[K]["p90"] for K in СІТКА_K])
катастрофи = np.array([статистика[K]["катастроф"] for K in СІТКА_K])
fig, (ліва, права) = plt.subplots(1, 2, figsize=(13, 4.4))
ліва.plot(СІТКА_K, медіани, marker="o", lw=2.2, color="teal", label="медіана")
ліва.plot(СІТКА_K, процентилі, marker="o", lw=2.2, color="crimson", label="90-й процентиль")
ліва.axvline(4, color="gray", ls="--", lw=1.5)
ліва.set_xscale("symlog"); ліва.set_yscale("log")
ліва.set_xlabel("кількість вузлів K"); ліва.set_ylabel("похибка (лог. шкала)")
ліва.set_title("Медіана спадає монотонно, хвіст — ні"); ліва.legend(); ліва.grid(alpha=.25, which="both")
права.plot(СІТКА_K, 100 * катастрофи, marker="o", lw=2.2, color="darkorange")
права.axvline(4, color="gray", ls="--", lw=1.5)
права.set_xscale("symlog")
права.set_xlabel("кількість вузлів K"); права.set_ylabel("% розкладів із похибкою > 10")
права.set_title("Частка катастрофічних розкладів"); права.grid(alpha=.25)
plt.tight_layout(); plt.show()
медіана_на_порозі = статистика[4]["медіана"]
найгірша_на_порозі = статистика[4]["найгірша"]
print(f"на порозі K = 4: типовий розклад дає похибку {медіана_на_порозі:.2f},")
print(f"а найгірший із {КІЛЬКІСТЬ_РОЗКЛАДІВ} розкладів — {найгірша_на_порозі:.3g}.")
print(f"Середнє арифметичне при цьому дорівнює {статистика[4]['середнє']:.3g} —")
print("це число не має жодного сенсу: воно цілком визначене одним розкладом")
print("і при наступному запуску буде іншим. Саме це й означає «нескінченне сподівання».")
assert найгірша_на_порозі > 20 * медіана_на_порозі, "хвіст мав бути важким!"
assert статистика[4]["катастроф"] > 10 * статистика[20]["катастроф"] or \
статистика[20]["катастроф"] == 0.0
assert медіани[-1] < медіани[СІТКА_K.index(4)], "медіана мала спадати за порогом!"
print("\n✅ біля порогу медіана пристойна, а хвіст важкий: катастрофи трапляються")
print(" у кількох відсотках розкладів. Від K = 12 вони зникають узагалі.")
на порозі K = 4: типовий розклад дає похибку 1.38, а найгірший із 300 розкладів — 1.69e+13. Середнє арифметичне при цьому дорівнює 5.65e+10 — це число не має жодного сенсу: воно цілком визначене одним розкладом і при наступному запуску буде іншим. Саме це й означає «нескінченне сподівання». ✅ біля порогу медіана пристойна, а хвіст важкий: катастрофи трапляються у кількох відсотках розкладів. Від K = 12 вони зникають узагалі.
7. За порогом: усі проходять крізь точки, але по-різному¶
При K = 12 параметрів p = 14, а рівнянь усього шість. Множина розвʼязків, що проходять крізь усі точки, — це вже не одна точка, а восьмивимірний плоский шматок простору.
Подивимось на кілька законних інтерполянтів із цієї множини і порівняємо їх із розвʼязком мінімальної норми.
In [11]:
rng = np.random.default_rng(4242)
вузли_12 = np.sort(rng.uniform(ЛІВА_МЕЖА, ПРАВА_МЕЖА, 12))
коеф_мін = навчити(вузли_12)
# напрямки, уздовж яких можна рухатись, не змінюючи жодного значення в точках даних
система_на_гамма = ДОПОВНЕННЯ.T @ матриця_плану(X_ДАНІ, вузли_12)[:, 2:]
_, сингулярні, права_матриця = np.linalg.svd(система_на_гамма)
ядро = права_матриця[len(сингулярні):]
fig, ax = plt.subplots(figsize=(10, 4.6))
ax.plot(сітка, істинна_функція(сітка), color="gray", ls="--", lw=1.8, label="істина")
print(f"розмірність множини інтерполянтів: {ядро.shape[0]}")
print(f"\n{'варіант':>22} {'‖γ‖':>10} {'сума залишків':>16}")
print(f"{'мінімальна норма':>22} {норма_гамма(коеф_мін):10.3f} "
f"{сума_залишків(коеф_мін, вузли_12):16.2e}")
for спроба in range(4):
зсув_гамма = rng.normal(size=ядро.shape[0]) @ ядро * 8
інший = коеф_мін.copy()
інший[2:] = інший[2:] + зсув_гамма
# α і β добираємо заново, щоб крива знову пройшла крізь точки
залишок = Y_ДАНІ - матриця_плану(X_ДАНІ, вузли_12)[:, 2:] @ інший[2:]
інший[:2], *_ = np.linalg.lstsq(пряма_на_точках, залишок, rcond=None)
print(f"{'інший інтерполянт':>22} {норма_гамма(інший):10.3f} "
f"{сума_залишків(інший, вузли_12):16.2e}")
ax.plot(сітка, прогноз(інший, вузли_12, сітка), color="darkorange", lw=1.4, alpha=.8)
ax.plot(сітка, прогноз(коеф_мін, вузли_12, сітка), color="teal", lw=3, label="мінімальна норма")
ax.scatter(X_ДАНІ, Y_ДАНІ, s=32, color="crimson", zorder=4, label="дані")
ax.set_ylim(-40, 70)
ax.set_title("Усі криві проходять крізь ті самі шість точок")
ax.legend(); ax.grid(alpha=.25)
plt.tight_layout(); plt.show()
print("\nЗалишки нульові в усіх — дані не дають жодної підстави когось із них обрати.")
print("Вибір робить критерій мінімальної норми, і саме він витягує найспокійнішу криву.")
розмірність множини інтерполянтів: 8
варіант ‖γ‖ сума залишків
мінімальна норма 10.296 4.26e-14
інший інтерполянт 24.133 4.84e-14
інший інтерполянт 25.167 2.49e-14
інший інтерполянт 17.323 3.93e-14
інший інтерполянт 21.276 3.91e-14
Залишки нульові в усіх — дані не дають жодної підстави когось із них обрати. Вибір робить критерій мінімальної норми, і саме він витягує найспокійнішу криву.
8. Норма спадає як 1/√K¶
У лекції ми вивели: ‖γ‖² ≈ (L/K)·∫g″(x)²dx. Множник L/K від функції не залежить, тому норма має спадати як 1/√K. Перевіримо це прямою прямою на логарифмічних осях: нахил має вийти близьким до −0.5.
In [12]:
# щоб побачити чистий закон 1/√K, треба відійти від порогу: там норма ще роздута
K_ДЛЯ_НОРМИ = [10, 20, 40, 80, 160, 320]
rng = np.random.default_rng(3)
норми_медіанні = []
for K in K_ДЛЯ_НОРМИ:
по_розкладах = [норма_гамма(навчити(np.sort(rng.uniform(ЛІВА_МЕЖА, ПРАВА_МЕЖА, K))))
for _ in range(40)]
норми_медіанні.append(np.median(по_розкладах))
норми_медіанні = np.array(норми_медіанні)
# нахил прямої в координатах (log K, log ‖γ‖) — це і є показник степеня
нахил, зсув = np.polyfit(np.log(K_ДЛЯ_НОРМИ), np.log(норми_медіанні), 1)
fig, ax = plt.subplots(figsize=(9, 4.4))
біля_порогу = [K for K in СІТКА_K if K > 0]
ax.plot(біля_порогу, [статистика[K]["норма"] for K in біля_порогу], marker="o", lw=2.2,
color="crimson", label="медіанна ‖γ‖ (зона порогу)")
ax.plot(K_ДЛЯ_НОРМИ, норми_медіанні, marker="s", lw=2.2, color="darkorange",
label="медіанна ‖γ‖ (далеко за порогом)")
ax.plot(K_ДЛЯ_НОРМИ, np.exp(зсув) * np.array(K_ДЛЯ_НОРМИ, dtype=float) ** нахил, ls="--",
lw=2, color="teal", label=f"пряма з нахилом {нахил:.2f}")
ax.axvline(4, color="gray", ls="--", lw=1.5)
ax.set_xscale("log"); ax.set_yscale("log")
ax.set_xlabel("кількість вузлів K (лог. шкала)"); ax.set_ylabel("‖γ‖ (лог. шкала)")
ax.set_title("Максимум норми — у зоні порогу, далі спад як 1/√K")
ax.legend(fontsize=9); ax.grid(alpha=.25, which="both")
plt.tight_layout(); plt.show()
print(f"{'K':>6} {'медіанна ‖γ‖':>15}")
for K, норма in zip(K_ДЛЯ_НОРМИ, норми_медіанні):
print(f"{K:>6} {норма:15.3f}")
print(f"\nвиміряний нахил = {нахил:.3f}, теорія обіцяє −0.5")
print(f"норма на порозі K = 4: {статистика[4]['норма']:.1f} — максимум усієї кривої")
assert abs(нахил + 0.5) < 0.1, "спад норми розійшовся з теорією 1/√K!"
print("\n✅ норма справді спадає як 1/√K — рівно те, що дає викладка з лекції")
K медіанна ‖γ‖
10 14.133
20 9.522
40 6.275
80 4.370
160 3.028
320 2.116
виміряний нахил = -0.548, теорія обіцяє −0.5
норма на порозі K = 4: 35.6 — максимум усієї кривої
✅ норма справді спадає як 1/√K — рівно те, що дає викладка з лекції
9. Границя: натуральний кубічний сплайн¶
Головний результат другої частини. Мінімізація ‖γ‖² — це насправді мінімізація ∫g″(x)²dx, тобто пошук найменш вигнутої кривої крізь точки. А функція, яка проходить крізь задані точки й мінімізує цей інтеграл, давно відома з чисельних методів: це натуральний кубічний сплайн.
Отже теорема звучить так: наші моделі мінімальної норми при K → ∞ збігаються
до натурального кубічного сплайна. scipy вміє будувати його однією командою —
звіримось.
In [13]:
from scipy.interpolate import CubicSpline
сплайн = CubicSpline(X_ДАНІ, Y_ДАНІ, bc_type="natural")
сітка_порівняння = np.linspace(ЛІВА_МЕЖА, ПРАВА_МЕЖА, 401)
значення_сплайна = сплайн(сітка_порівняння)
rng = np.random.default_rng(77)
K_ДЛЯ_ЗБІЖНОСТІ = [4, 8, 16, 32, 64, 128, 256, 512]
розбіжності = []
print(f"{'K':>5} {'p':>5} {'макс. розбіжність зі сплайном':>32} {'похибка проти параболи':>25}")
for K in K_ДЛЯ_ЗБІЖНОСТІ:
по_розкладах, похибки = [], []
for _ in range(9): # медіана по кількох розкладах вузлів
вузли = np.sort(rng.uniform(ЛІВА_МЕЖА, ПРАВА_МЕЖА, K))
коеф = навчити(вузли)
по_розкладах.append(np.max(np.abs(прогноз(коеф, вузли, сітка_порівняння) - значення_сплайна)))
похибки.append(похибка_проти_істини(коеф, вузли))
розбіжності.append(np.median(по_розкладах))
print(f"{K:>5} {K + 2:>5} {np.median(по_розкладах):32.4f} {np.median(похибки):25.4f}")
похибка_сплайна = np.mean(np.abs(значення_сплайна - істинна_функція(сітка_порівняння)))
print(f"\nпохибка самого сплайна проти параболи: {похибка_сплайна:.4f}")
assert розбіжності[-1] < 0.5, "модель не зійшлася до сплайна!"
assert розбіжності[-1] < розбіжності[0] / 5, "розбіжність мала помітно спадати!"
print("\n✅ розбіжність зі сплайном упала з "
f"{розбіжності[0]:.2f} до {розбіжності[-1]:.3f} — моделі злилися")
K p макс. розбіжність зі сплайном похибка проти параболи
4 6 4.9011 1.4139
8 10 2.7348 0.8919
16 18 0.8405 0.2960
32 34 0.6386 0.2224
64 66 0.2404 0.1688
128 130 0.3301 0.1566
256 258 0.1718 0.1480
512 514 0.1741 0.1303
похибка самого сплайна проти параболи: 0.1330
✅ розбіжність зі сплайном упала з 4.90 до 0.174 — моделі злилися
In [14]:
fig, (ліва, права) = plt.subplots(1, 2, figsize=(13, 4.4))
rng = np.random.default_rng(77)
for K, колір, товщина in [(4, "crimson", 1.6), (16, "darkorange", 1.8), (400, "teal", 2.6)]:
вузли = np.sort(rng.uniform(ЛІВА_МЕЖА, ПРАВА_МЕЖА, K))
коеф = навчити(вузли)
ліва.plot(сітка_порівняння, прогноз(коеф, вузли, сітка_порівняння), lw=товщина,
color=колір, label=f"мінімальна норма, K = {K}")
ліва.plot(сітка_порівняння, значення_сплайна, color="black", ls="--", lw=2.4,
label="натуральний кубічний сплайн")
ліва.scatter(X_ДАНІ, Y_ДАНІ, s=32, color="black", zorder=4)
ліва.set_ylim(-8, 32)
ліва.set_title("Що більше вузлів, то ближче до сплайна"); ліва.legend(fontsize=9)
ліва.grid(alpha=.25)
права.plot(K_ДЛЯ_ЗБІЖНОСТІ, розбіжності, marker="o", lw=2.2, color="teal")
права.set_xscale("log"); права.set_yscale("log")
права.set_xlabel("кількість вузлів K (лог. шкала)")
права.set_ylabel("макс. розбіжність (лог. шкала)")
права.set_title("Розбіжність зі сплайном тане"); права.grid(alpha=.25, which="both")
plt.tight_layout(); plt.show()
print("Модель із сотнями параметрів на шести точках даних робить рівно те саме,")
print("що робить сплайн — давно вивчений, гарантовано розумний інтерполятор.")
print("Ось тут подвійний спуск і перестає бути парадоксом: уся загадка зводиться")
print("до того, що кубічний сплайн наближає параболу краще, ніж ламана з чотирьох сегментів.")
Модель із сотнями параметрів на шести точках даних робить рівно те саме, що робить сплайн — давно вивчений, гарантовано розумний інтерполятор. Ось тут подвійний спуск і перестає бути парадоксом: уся загадка зводиться до того, що кубічний сплайн наближає параболу краще, ніж ламана з чотирьох сегментів.
💻 Завдання¶
🟢 Рівень 1 — разом¶
- Заміни параболу на іншу гладку функцію (наприклад,
np.sin(x)на відрізку [0, 5]). Чи лишається пік на порозі? Чи так само модель збігається до сплайна? - Додай до
Y_ДАНІдрібний шум і перебудуй розділ 6. Що змінилось у медіані, а що — у хвості?
🟡 Рівень 2 — самостійно¶
- Зроби вузли рівномірними замість випадкових (
np.linspace). Куди подівся важкий хвіст на порозі й чому? Це і є відповідь на питання, що саме робить модель жорсткою. - Побудуй залежність похибки від кількості точок даних N (при фіксованому K = 20). Де тепер стоїть поріг?
🔴 Рівень 3 — виклик¶
- Реалізуй мінімізацію ∫g″(x)²dx напряму — через дискретну другу похідну на щільній
сітці й
scipy.optimize— і покажи, що результат збігається з нашим розвʼязком мінімальної норми при великому K. - Замість базису ½·|x − t| візьми базис із гладких хвиль (як у практиці 36). Що тепер означає «мінімальна норма»? До якої кривої збігається модель? Це і є ілюстрація індуктивного зміщення архітектури.
🧪 Самоперевірка¶
1. Чому α і β не входять у норму, яку ми мінімізуємо?
відповідь
Пряма α + βx — «безкоштовна» частина моделі: вона не має зламів і не робить криву вигнутою. Мінімізуємо ми, по суті, вигнутість (∫g″²), а в неї пряма не вносить нічого. Штрафувати α і β означало б без потреби притискати модель до нуля.2. На порозі отримали похибку 0.9, а сума залишків — 10⁻¹³. Модель хороша?
відповідь
Ні. Нульові залишки означають лише, що крива пройшла крізь шість точок — вони нічого не кажуть про те, що діється між точками. У розділі 5 ми бачили криві з нульовими залишками й розмахом у сотні одиниць. Дивитись треба на похибку проти істини, а не на залишки.3. Чому середню помилку на порозі рахувати безглуздо?