Восстановление коэффициентов¶

Задача на CodeRun. Даны 1000 точек $(x, y)$:

$$y = (a\sin x + b\ln x)^2 + c\,x^2 + \text{шум}.$$

  • $x$, $y$ — данные из файла, они известны;
  • $a$, $b$, $c > 0$ — неизвестные коэффициенты, их нужно найти с точностью $10^{-2}$;
  • шум — отклонения коэффициентов не больше $0{,}001$. Это в 10 раз меньше нужной точности, поэтому шум просто не учитываем.

Как решаем. Берём модель $\hat y = (a\sin x + b\ln x)^2 + c\,x^2$ и подбираем $a$, $b$, $c$ так, чтобы предсказания $\hat y_i$ были как можно ближе к $y_i$. Для этого считаем ошибку, находим её производные по $a$, $b$, $c$ и запускаем градиентный спуск.

Автор: Алексей Чернуха. CC BY 4.0.

In [ ]:
import urllib.request
from pathlib import Path

import numpy as np

URL = "https://assets.contest.yandex.net/testsys/statement-file?hash=eyJhbGciOiJkaXIiLCJlbmMiOiJBMjU2R0NNIn0..S1G3epXKnK0cjv51.l6h0gXNaY9YG2YcFEE8q8DwFzfhZY2OQmoOu-Px4yuNzMK2FyfcFJhosWW9cZ5XUMGOUA0wVdiMZ0U3D5NGyk-OOl7LB37A.7Q6iR5TS6Y7SPME_yxW64w"
if not Path("data.csv").exists():
    urllib.request.urlretrieve(URL, "data.csv")

x, y = np.loadtxt("data.csv", delimiter=",", unpack=True)
print(len(x), "точек")

Если скачать не получилось, возьмите data.csv со страницы задачи и положите рядом с ноутбуком.

1. Производные¶

Ошибка — MSE, средний квадрат отклонения модели от данных:

$$L = \frac1n\sum_{i=1}^{n} \Bigl((a\sin x_i + b\ln x_i)^2 + c\,x_i^2 - y_i\Bigr)^2.$$

  • $i$ — номер точки, $n = 1000$ — число точек;
  • в большой скобке — на сколько модель ошиблась в точке $i$: предсказание минус правильный ответ;
  • $L$ зависит только от $a$, $b$, $c$: данные $x_i$, $y_i$ — просто числа.

Найдите на бумаге $\dfrac{\partial L}{\partial a}$, $\dfrac{\partial L}{\partial b}$, $\dfrac{\partial L}{\partial c}$.

💡 Подсказка

Производная суммы — сумма производных, а $\frac1n$ — просто множитель. Дальше цепное правило, причём дважды: снаружи большой квадрат, а внутри него ещё один квадрат $(a\sin x_i + b\ln x_i)^2$.

✅ Решение по шагам

Считаем производную по $a$. Числа $b$, $c$, $x_i$, $y_i$ при этом — константы.

Шаг 1. Заносим производную под знак суммы. Производная суммы равна сумме производных, множитель $\frac1n$ выносится:

$$\frac{\partial L}{\partial a} = \frac1n\sum_{i=1}^{n} \frac{\partial}{\partial a}\Bigl((a\sin x_i + b\ln x_i)^2 + c\,x_i^2 - y_i\Bigr)^2.$$

Шаг 2. Внешний квадрат. По цепному правилу $(u^2)' = 2u\cdot u'$, где $u$ — вся большая скобка:

$$\frac{\partial L}{\partial a} = \frac1n\sum_{i=1}^{n} 2\Bigl((a\sin x_i + b\ln x_i)^2 + c\,x_i^2 - y_i\Bigr)\cdot\frac{\partial}{\partial a}\Bigl((a\sin x_i + b\ln x_i)^2 + c\,x_i^2 - y_i\Bigr).$$

Шаг 3. Производная того, что внутри большой скобки. Берём по слагаемым:

  • $\dfrac{\partial}{\partial a}\,c\,x_i^2 = 0$ — в нём нет $a$;
  • $\dfrac{\partial}{\partial a}\,y_i = 0$ — тоже нет $a$;
  • $(a\sin x_i + b\ln x_i)^2$ — снова квадрат, снова цепное правило:

$$\frac{\partial}{\partial a}(a\sin x_i + b\ln x_i)^2 = 2(a\sin x_i + b\ln x_i)\cdot\frac{\partial}{\partial a}(a\sin x_i + b\ln x_i) = 2(a\sin x_i + b\ln x_i)\sin x_i.$$

Последний переход: у $a\sin x_i$ производная по $a$ равна $\sin x_i$ (как у $3a$ производная $3$), а $b\ln x_i$ от $a$ не зависит.

Шаг 4. Подставляем в шаг 2. Двойки перемножаются: $2\cdot 2 = 4$.

$$\frac{\partial L}{\partial a} = \frac4n\sum_{i=1}^{n}\Bigl((a\sin x_i + b\ln x_i)^2 + c\,x_i^2 - y_i\Bigr)(a\sin x_i + b\ln x_i)\sin x_i.$$

По $b$ — те же шаги. Отличие только в шаге 3: $\dfrac{\partial}{\partial b}(a\sin x_i + b\ln x_i) = \ln x_i$.

$$\frac{\partial L}{\partial b} = \frac4n\sum_{i=1}^{n}\Bigl((a\sin x_i + b\ln x_i)^2 + c\,x_i^2 - y_i\Bigr)(a\sin x_i + b\ln x_i)\ln x_i.$$

По $c$ шаги 1 и 2 те же, а шаг 3 проще: $(a\sin x_i + b\ln x_i)^2$ и $y_i$ не содержат $c$, а $\dfrac{\partial}{\partial c}\,c\,x_i^2 = x_i^2$. Второго квадрата нет, поэтому множитель остаётся $2$:

$$\frac{\partial L}{\partial c} = \frac2n\sum_{i=1}^{n}\Bigl((a\sin x_i + b\ln x_i)^2 + c\,x_i^2 - y_i\Bigr)\,x_i^2.$$

В коде $\frac1n\sum_i$ — это np.mean.

2. Градиентный спуск¶

Каждый шаг сдвигает коэффициенты против градиента: $a \leftarrow a - \eta\,\dfrac{\partial L}{\partial a}$, и так же для $b$ и $c$. Здесь $\eta$ — шаг (lr).

Одни и те же куски повторяются во всех трёх формулах, поэтому в коде они посчитаны заранее:

  • s — это $a\sin x_i + b\ln x_i$;
  • err — это большая скобка $(a\sin x_i + b\ln x_i)^2 + c\,x_i^2 - y_i$.

Оба — массивы сразу по всем точкам. Допишите градиенты по формулам из пункта 1 и запустите.

In [ ]:
a, b, c = 1.0, 1.0, 1.0
lr = 1e-4

for step in range(5000):
    s = a * np.sin(x) + b * np.log(x)
    err = s**2 + c * x**2 - y

    grad_a = ...  # ваш код
    grad_b = ...
    grad_c = ...

    a -= lr * grad_a
    b -= lr * grad_b
    c -= lr * grad_c

    if step % 1000 == 0:
        print(step, np.mean(err**2))

print(f"{a:.2f} {b:.2f} {c:.2f}")
✅ Решение 
grad_a = 4 * np.mean(err * s * np.sin(x))
grad_b = 4 * np.mean(err * s * np.log(x))
grad_c = 2 * np.mean(err * x**2)
🤔 Ошибка растёт или получается nan?

Скорее всего, шаг слишком большой. $x$ доходит до 10, значит $x^2$ — до 100, и градиент по $c$ огромный: уже при lr = 1e-3 спуск разлетается. Верните lr = 1e-4.

Если шаг 1e-4, а ошибка всё равно растёт, проверьте формулы градиентов: знаки и множители 4 и 2.

Последняя строка — ответ для CodeRun.