👨‍💻Всем сап! В прошлом посте разбирали градиентный спуск, теперь самое интересное ➡️ КОД (пост оказался сложнее, чем я думал, приходится сидеть и думать 🤔🤔)

🔵Функцию берем вот эту 🟰 f(x) = x² - 5x + 5 🟰

🔵Теперь производная 🟰 f'(x) = 2x - 5 🟰


0⃣ Библиотеки:

import time
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np

1⃣ Функции, которые я прописал выше (просто переписываем):

def f(x): # Функция f(x)
    return x * x - 5 * x + 5

def df(x): # Проивзодная функции f(x)
    return 2 * x - 5

2⃣ Теперь вводим некоторые переменные:


🔵(Число итераций (N); — проще говоря, сколько шагов градиентного спуска мы делаем

🔵Начальное значение (xx); — стартовая точка, c нее мы начнем ползти по антиградиенту (это вниз если что, градиент указывает направление наискорейшего роста функции, антиградиент — направление наискорейшего убывания)

🔵Шаг сходимости (lmd) — проще говоря, это learning rate (длина шага).

🟰 xₖ₊₁ = xₖ - λf'(xₖ) 🟰 — вот по такой формуле делает шаг градиентный спуск. Если мы подставим 🟰 f'(x) = 2x - 5 🟰то получим 2.5

‼️ Кстати, эту формулу также записывают как 🟰 xₖ₊₁ = xₖ - λ∇ f'xₖ) 🟰Разница лишь в том, что в первом случае f'(x) — это обычная производная по одной переменной, а в другом — ∇ f(x) — вектор частных производных (как итог, формула одна и та же)

N = 20 
xx = 0 
lmd = 0.1

👨‍💻 Итак, мы хотим найти x, где функция потерь f(x) минимальная, формулу шага градиентного спуска мы уже разобрали, хотя, можно еще один момент разобрать:

∇f(x) — стрелка вверх (т.е куда быстрее всего растет f) — это направление наискорейшего роста

-∇f(x) — направление наискорейшего убывания (в комментах материал)

Можем кстати подставить в нашу формулу пример 🟰 xₖ₊₁ = xₖ - λf'(xₖ) 🟰 ➡️ получаем 🟰 xₖ₊₁ = xₖ - λ * (2x - 5)) 🟰 ➡️ 🟰(1 - 2λ) * xₖ + 5λ🟰 — это правило обновления итераций (наша функция)

‼️САМА СУТЬ — НАЙТИ ТОЧКУ, ГДЕ ДВИЖЕНИЕ ПРЕКРАЩАЕТСЯ (итерации прекращаются) — еще проще, найти глобальный минимум.

🔵А так, я нашел, что глобальный минимум можно найти кучей способов, самый простой — вершина параболы: напомню, у нас функция 🟰 f(x) = x² - 5x + 5 🟰

🔵 x (вершина параболы) = -b / 2a = 2.5 (просто подставляем наши значения, вот вам и глобальный минимум)

Совсем скоро продолжим 💤

👨‍💻 Итак, продолжаем:

⚙️Наш код на данный момент выглядит так:

import numpy as np
import matplotlib.pyplot as plt

def f(x):   # f(x) = x^2 - 5x + 5 # парабола 
    return x*x - 5*x + 5

def df(x):  # f'(x) = 2x - 5 # наклон графика в точке x 
    return 2*x - 5

N   = 20    # число итераций
x   = 0.0   # стартовая точка x0
lr  = 0.1   # шаг обучения (λ)

🔵N — это количество обновлений

🔵lr — это длина шага

❓ Возможно, вы заметили, что я убрал библиотеку time, ее я заменил привычной нам matplotlib, суть не поменялась (все это не влияет на наш алгоритм, исключительно для плавной анимации спуска)

3⃣ Создаем массив точек из 400 точек, по этим точкам и будем рисовать кривую f(x)

xs = np.linspace(-1, 5, 400)

4⃣ Включаем интерактивный режим Matplotlib (для анимации в цикле). ‼️ Важно, анимация у меня в колабе не отобразилась, мб это особенность блокнотов, в VS все нормально.

🔵Что касаемо остальных строк, fig, ax = plt.subplots() — будет создавать окошко Figure и систему координат Axes (сетка, подписи и две оси — x, y)

🔵ax.grid(True) — включает сетку на осях

🔵ax.plot(xs, f(xs), label="f(x)") — рисует линию функции f(x) по сетке xs. ‼️ В общем, мы создавали массив из 400 точек, так вот когда мы вызываем f(xs), Numpy поэлементно будет считать y = f(x) для каждого из этих 400 точек

🔵pt = ax.scatter([x], [f(x)], c='red', s=50) — рисует красную точку в текущем положении (x, f(x)) для визуализации текущего x на каждом шаге спуска.

🔵ax.set_xlim(xs.min(), xs.max()) — фиксирует границы оси X от минимума до максимума нашей сетки (таким образом, будем отображаться ровно тот интервал, где мы нарисовали кривую + масштаб во время анимации не будет прыгать)

🔵ax.set_ylim(f(xs).min()-1, f(xs).max()+1) — фиксирует границы оси Y по значениям функции. Почему берём min/max именно от xs и f(xs) ❓ Потому что это левые/правые границы, где мы считали и рисовали функцию.

🔵plt.show(block=False) — показывает окно не блокируя выполнение кода.

plt.ion()
fig, ax = plt.subplots()
ax.grid(True)
ax.plot(xs, f(xs), label="f(x)")
pt = ax.scatter([x], [f(x)], c='red', s=50)
ax.set_xlim(xs.min(), xs.max())
ax.set_ylim(f(xs).min()-1, f(xs).max()+1)
plt.show(block=False)

5⃣ Самое важное: градиентный спуск и анимация

🔵Это и есть обновление по формуле 🟰 xₖ₊₁ = xₖ - λ * (2x - 5)) 🟰Если lr = 0.5 ➡️ сразу попадаем в глобальный минимум

🔵pt — это объект, который вернул 🟰 pt = ax.scatter([x], [f(x)], c='red', s=50) 🟰

for _ in range(N):
    x = x - lr * df(x)   # шаг: x_{k+1} = x_k - λ f'(x_k)
    pt.set_offsets([[x, f(x)]])   # кладет новые (x, y) в объект точки (грубо говоря, это обновление уже нарисованной точки, будет переносить маркер в новые корды)
    fig.canvas.draw_idle() # Рисует содержимое окна ; перерисовывает
    plt.pause(0.03) # фигура перерисовывается + пазуа

6⃣ Фиксируем результаты:

plt.ioff()   # выключить интерактив
ax.scatter([x], [f(x)], c='blue', s=50, label="финал")  # финальная точка (это куда мы пришли, глобальный минимум)
ax.axvline(2.5, ls='--', label="x* = 2.5")   # вертикальная линия минимума
ax.legend()
print(x)
plt.show()

‼️ Теперь пройдемся по логике кода (САМОЕ САМОЕ ВАЖНОЕ): мы взяли 400 точек по оси x ➡️подставляем каждую точку в функцию f(x) ➡️ ax.plot соединяет точки линией (xs — это x корды, f(xs) — это y корды) ➡️ далее, у нас есть текущая точка x, считаем производную f'(x) и делаем шаг вниз (это антиградиент) ➡️ f'(x) задает направление и масштаб, λ — длина шага

👨‍💻 Всем ку! Заждались? Сегодня разбираем интересное задание, в последующем разберем несколько решений.

⬇️ ВЕСЬ МАТЕРИАЛ В КОММЕНТАРИЯХ ⬇️

⚙️ Итак, надеюсь, задание вы прочитали, теперь приступаем:

0⃣ Что нам вообще нужно ❓

‼️ Начнем с формулировки "восстановления/аппроксимации" функции. Это значит, что нам нужно заменить (ИЛИ ПРИБЛИЗИТЬ) нашу функцию🟰 f(x) = 0.1x² - sin(x) + 5 🟰 кубическим полиномом (это который 🟰a(x) = w₀ + w₁... 🟰таким образом, чтобы подобранные веса w делали среднеквадратическую ошибку MSE между a(x) — (тот самый полином) и f(x) МИНИМАЛЬНОЙ НА СЕТКЕ [-5, 5]

‼️ Если не поняли, могу сказать по-другому: нам нужно подобрать веса w в кубическом полиноме a(x) таким образом, чтобы MSE между a(x) и f(x) на отрезке [-5, 5] была МИНИМАЛЬНЫМ.

1⃣ В общем, подбираем коэффициенты w, но все по порядку:

🔵Нам нужны данные для обучения и сама функция: у нас по условию уже есть точки x на отрезке [-5. 5] и тру-значения y = f(x), функция также дана

import numpy as np

# исходная функция, которую нужно аппроксимировать моделью a(x)
def func(x):
    return 0.1 * x**2 - np.sin(x) + 5


# здесь объявляйте необходимые функции


coord_x = np.arange(-5.0, 5.0, 0.1) # значения по оси абсцисс [-5; 5] с шагом 0.1
coord_y = func(coord_x) # значения функции по оси ординат

sz = len(coord_x)  # количество значений функций (точек)
eta = np.array([0.1, 0.01, 0.001, 0.0001]) # шаг обучения для каждого параметра w0, w1, w2, w3
w = np.array([0., 0., 0., 0.]) # начальные значения параметров модели
N = 200 # число итераций градиентного алгоритма

👨‍💻 Все что вы видите сверху — это все, что у нас есть, теперь будем думать:

2⃣ Наша следующая функция — features, она будет строить матрицу признаков S с колонками [1, x, x², x³].

def features(x):
    x = np.asarray(x)
    return np.stack([np.ones_like(x), x, x**2, x**3], axis=1)

🔵Функция np.asarray() в библиотеке NumPy — используется для преобразования входных данных в объект ndarray.

🔵ndarray (N-dimensional array) — класс многомерных массивов в библиотеке NumPy

🔵Входные данные могут быть в виде списков, кортежей или других объектов, похожих на массивы. Если входные данные уже являются массивом NumPy, функция не создаёт копию

🔵В нашем коде делает из входа (список/массив) numpy-вектор формы (n,).
А нужно это, чтобы дальше работали векторные операции x**2, x**3

🔵np.ones_like(x) — делает массив единиц такой же формы и типа, как x. Это столбец для свободного члена w₀ — форма также (n, ) — т.е вектор из n единиц. В нашей функции признаков этот вектор из единиц идёт как первый столбец

Что такое формы (n,), (n,1), (n,4) ❓

Я уже как-то разбирал данный момент, но вернуться стоит:

🟢 (n,) — одномерный вектор длины n — это просто линейный массив...

💤 СКОРО ПРОДОЛЖИМ

🔤🔤🌹🌷 5️⃣

👨‍💻 В дополнении к формам массивов:

🟢(n,) — одномерный вектор длины n. Нет строк/столбцов, просто линейный массив, форма (3, ) буде выглядеть как то так:

x = np.array([10., 20., 30.])
x.shape # (3, )
------------
[10. 20. 30.]

🟢(n, 1) — двумерная матрица: n строк, 1 столбец, массив форму (3, 1) буде выглядеть вот так:

x_col = x.reshape(-1, 1)
x_col.shape # (3, 1)
------------
[[10.]
 [20.]
 [30.]]

🟢(n, 4) — матрица: n строк, 4 столбца

S = np.column_stack([np.ones_like(x), x, x**2, x**3])
S.shape
-------------
[[  1.  10. 100. 1000.]
 [  1.  20. 400. 8000.]
 [  1.  30. 900. 27000.]]

⬇️ ПОЛЕЗНЫЕ ССЫЛОЧКИ ВНИЗУ ⬇️

👨‍💻 Сегодня продолжим разбирать ГРАДИЕНТНЫЙ СПУСК 📉
➿➿➿➿➿➿
⬇️ ПРОШЛЫЕ ПОСТЫ ДЛЯ ОЗНАКОМЛЕНИЯ ⬇️
➿➿➿➿➿➿
0⃣ ОБЩИЕ ПОНЯТИЯ И ФУНЦИЯ ПОТЕРЬ

1⃣ ВЫПУКЛЫЕ И НЕВЫПУКЛЫЕ ФУНКЦИИ

2⃣ ПРАКТИКА НА КОДЕ

3⃣ РАЗБОР ФОРМУЛ И СТОХАСТИЧЕСКИЙ ГРАДИЕНТНЫЙ СПУСК

4⃣ РАБОТА АЛГОРИТМА И ПРИМЕРЫ

5⃣ ПОДКЛЮЧАЕМ КОД

6⃣ ПРОДОЛЖАЕМ ПИСАТЬ КОД

7⃣ РАЗБИРАЕМ ПОЛНОЦЕННОЕ ЗАДАНИЕ

8⃣ ФОРМЫ МАССИВОВ

👨‍💻 Следующая функция в нашем коде — predict()

def predict(w, x):
    S = features(x)
    return S @ w    

❓ Итак, че тут происходит? Функция вычисляет 🟰 a(x) = Sw 🟰. Откуда взялась эта запись? Это просто компактная запись кубического полинома 🟰 a(x) = w₀ + w₁.. 🟰. S — это по сути вектор признаков [1, x, x², x³] — он дан в условии.

    S = features(x)

🔵Выше мы вызываем функцию features(), которую разбирали в нескольких постах — ЗДЕСЬ МИНИ-НАВИГАЦИЯ/ В общем, features превращала каждый скаляр xᵢ в вектор признаков в вектор признаков [1, x, x², x³] — воот, и все эти векторы складываются в матрицу признаков S.

    return S @ w

⚙️ Строка выше — тут матричное умножение. Берем вектор весов w и умножаем его на матрицу S, на выходе predict() реализует формулу:

🟰 a(x) = w₀ + w₁x + w₂x² + w₃x³ = Sw 🟰

👨‍💻 Теперь код выглядит так:

import numpy as np

def func(x): # Функция
    return 0.1 * x ** 2 - np.sin(x) + 5

def features(x): # Матрица признков
    x = np.asarray(x)
    return np.stack([np.ones_like(x), x, x ** 2, x ** 3], axis = 1)

def predict(w, x): # a(x) = Sw
    S = features(x)
    return S @ w

🔤🔤🌹🌷 6️⃣

👨‍💻 Итак, следующая функция — risk() — по сути это просто MSE (среднеквадратичная) между предсказаниями и тру-значениями

def risk(w, x, y):
    y_hat = predict(w, x)  # считаем предсказания: S @ w
    n = len(x) # количество точек (должно равняться len(y))
    return (1.0 / n) * np.sum((y_hat - y) ** 2)  # MSE

⚙️ Теперь разбираем построчно:

y_hat = predict(w, x)

🔵Так ну строка выше это вектор предсказаний, мы вызываем функцию predict, которую разбирали — ТУТ.

❓ Что делает predict(w, x) внутри risk? Функция predict() строит матрицу признаков S = features(x) и умножает ее на w.

‼️ Возьмем например 'свои' w и x, пусть x = [0, 1, 2], w = [1, 2, 0, 0], истинные значения возьмем из функции func(x) = 0.1*x**2 - sin(x) + 5

🔵Признаки:

x           = [0, 1, 2]
ones_like   = [1, 1, 1]
x**2        = [0, 1, 4]
x**3        = [0, 1, 8]

🔵Матрица признаков S:

S = stack(ones, x, x**2, x**3, axis=1) =
[[1, 0, 0, 0],
 [1, 1, 1, 1],
 [1, 2, 4, 8]] # shape (3, 4) - форма массива

🔵Наши веса:

w = [1, 2, 0, 0]        # shape (4,)

🟢Теперь y_hat = S @ w (скалярно по строкам):

🔵для x=0: [1,0,0,0]·[1,2,0,0] = 1

🔵для x=1: [1,1,1,1]·[1,2,0,0] = 1 + 2 = 3

🔵для x=2: [1,2,4,8]·[1,2,0,0] = 1 + 4 = 5

Предсказания равны:

y_hat = [1, 3, 5]       # shape (3,)

🔵 Теперь считаем функцию поэлементно:

func(0) = 0.1*0 - sin(0) + 5 = 5.00000

func(1) = 0.1*1 - sin(1) + 5 ≈ 0.1 - 0.84147 + 5 = 4.25853

func(2) = 0.1*4 - sin(2) + 5 ≈ 0.4 - 0.90930 + 5 = 4.49070

🔵Ниже наши истинные y = func(x)

y ≈ [5.00000, 4.25853, 4.49070]

🔵 Что делает risk(w, x, y)?

Ошибка e = y_hat - y ≈ [1-5, 3-4.25853, 5-4.49070] ≈ [-4.00000, -1.25853, 0.50930]

Квадраты: e**2 ≈ [16.00000, 1.5839, 0.2594]

Сумма: ≈ 16 + 1.5839 + 0.2594 = 17.8433

Среднее (делим на n=3): Q ≈ 17.8433 / 3 ≈ 5.9478

‼️ То есть risk вернул MSE ≈ 5.95 для выбранных w, но позже я еще разъясню, всем сладких 💤

Друзья, занимался делами, теперь допишем пост для кода:

🧑‍💻 На данный момент наш код выглядит как-то так:

import numpy as np

def func(x):
    return 0.1 * x ** 2 - np.sin(x) + 5

def features(x):
    x = np.asarray(x)
    return np.stack([np.ones_like(x), x, x ** 2, x ** 3], axis = 1)

def predict(w, x):
    S = features(x)
    return S @ w

def risk(w, x, y):
    y_hat = predict(w, x)
    n = len(x)
    return (1.0 / n) * np.sum((y_hat - y) ** 2)

🟢 Давайте еще раз быстро пробежимся и закончим:

🔵 Функция? — Есть, за нее отвечает функция func()

🔵 Признаки мы тоже предусмотрели, за нее отвечает функция features(), входные признаки мы преобразовываем в numpy массив, и возвращаем матрицу признаков S (в матрице у нас столбец единиц — для bias, x, x², x³ ), короче говоря, features превращала каждый скаляр xᵢ в вектор признаков в вектор признаков.

🔵 Теперь нужна функция для прогноза — predict(). Функция вычисляет 🟰 a(x) = Sw 🟰— по сути, умножение нашей матрицы S (переменной мы присваиваем значение матрицу) на вектор признаков w.

🔵Теперь функция risk() — возвращает по-сути MSE — среднеквадратическая между предсказаниями и истинными значениями. Здесь подаем в функцию аж 3 аргумента — w (параметры модели, и это единственное, что мы меняем, risk() будет зависеть от w ⏺ x (входные данные по которым мы считаем признаки) и делаем предсказания y_hat — это вектор предсказаний (результат матричного умножения S @ w) ⏺ y — это истинные значения. Логика примерно такая: минимизируем функцию потерь по w, а x, y — это данные.