Руководство по PyTorch для новичков: создаём модель множественной регрессии с нуля

Практический PyTorch: строим трёхслойную нейросеть для множественной регрессии

Незадолго до того, как большие языковые модели (LLM) стали объектом всеобщего хайпа, между фреймворками для «классического» машинного обучения и фреймворками для глубокого обучения проходила почти видимая граница.

Про машинное обучение обычно говорили в контексте scikit-learn, XGBoost и им подобных, тогда как PyTorch и TensorFlow доминировали на сцене, когда речь заходила о глубоком обучении.

После взрывного роста интереса к ИИ я всё чаще вижу, что PyTorch заметно опережает TensorFlow по популярности. Оба фреймворка очень мощные и позволяют дата-сайентистам решать самые разные задачи, включая обработку естественного языка, что вновь подогрело интерес к глубокому обучению.

Однако в этой статье я не собираюсь говорить про NLP. Вместо этого мы разберём задачу множественной линейной регрессии с двумя целями:

• Изучить, как создать модель с использованием PyTorch.

• Поделиться нюансами линейной регрессии, о которых не всегда рассказывают в других туториалах.

Приступим.

Подготовка данных

Давайте я избавлю вас от ещё одной вычурной формальной дефиниции линейной регрессии. Скорее всего, вы уже видели её десятки раз в бесконечных туториалах по всему Интернету. Достаточно сказать так: когда у вас есть переменная Y, которую вы хотите предсказывать, и другая переменная X, которая может объяснить изменение Y с помощью прямой линии, — в сущности, это и есть линейная регрессия.

Набор данных

В этом упражнении мы будем использовать набор данных Abalone [1].

Nash, W., Sellers, T., Talbot, S., Cawthorn, A., & Ford, W. (1994). Abalone [Dataset]. UCI Machine Learning Repository. https://doi.org/10.24432/C55C7W

Согласно документации к набору данных, возраст морского ушка (абалон) определяют, разрезая раковину по конусу, затем окрашивая срез и считая количество колец под микроскопом — занятие скучное и очень трудоёмкое. Поэтому для предсказания возраста используют другие, более простые для измерения характеристики.

Итак, давайте загрузим данные. Дополнительно мы применим one-hot-кодирование к переменной Sex, так как это единственная категориальная переменная в наборе.

# Загрузка данных
from ucimlrepo import fetch_ucirepo
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style('darkgrid')
from feature_engine.encoding import OneHotEncoder

# Загружаем набор данных
abalone = fetch_ucirepo(id=1)

# данные (в виде pandas DataFrame)
X = abalone.data.features
y = abalone.data.targets

# One-Hot кодирование признака Sex
ohe = OneHotEncoder(variables=['Sex'])
X = ohe.fit_transform(X)

# Объединяем признаки и целевую переменную в один DataFrame
df = pd.concat([X, y], axis=1)

Вот набор данных.

Чтобы построить более качественную модель, давайте исследуем данные.

Исследуем данные

Первые шаги, которые я обычно выполняю при разведочном анализе данных, такие:

1. Проверить распределение целевой переменной.

# Анализ распределения целевой переменной
plt.hist(y)
plt.title('Rings [Target Variable] Distribution')

График показывает, что целевая переменная не имеет нормального распределения. Это может повлиять на качество регрессии, но обычно проблему можно сгладить с помощью степенных преобразований, например логарифмирования или преобразования Бокса–Кокса.

2. Посмотреть статистическое описание.

Статистика позволяет заметить важные характеристики — среднее значение, стандартное отклонение, а также быстро увидеть аномалии в минимальных или максимальных значениях. Объясняющие переменные выглядят вполне корректно: находятся в небольшом диапазоне и на одной шкале. Целевая переменная (Rings) — на другой шкале.

# Статистическое описание данных
df.describe()

Теперь проверим корреляции.

# Анализ корреляций
(df
 .drop(['Sex_M', 'Sex_I', 'Sex_F'], axis=1)
 .corr()
 .style
 .background_gradient(cmap='coolwarm')
)

Объясняющие переменные имеют среднюю и сильную корреляцию с Rings. Также видно, что есть мультиколлинеарность между Whole_weight и Shucked_weight, Viscera_weight и Shell_weight. Length и Diameter тоже коллинеарны. Позже можно проверить, улучшится ли модель, если их убрать.

sns.pairplot(df);

Когда мы строим диаграммы рассеяния в парах и смотрим на связь переменных с Rings, можно быстро заметить несколько проблем.

• Нарушено предположение гомоскедастичности. Это означает, что дисперсия ошибки непостоянна.

• Посмотрите, как точки на графиках образуют форму конуса: по мере роста X увеличивается разброс значений Y. Поэтому при прогнозировании Rings для больших значений X погрешность будет выше.

• Переменная Height имеет минимум два явных выброса при Height > 0.3.

Если удалить выбросы и преобразовать целевую переменную через логарифм, получим следующую диаграмму. Она лучше, но проблема гомо��кедастичности всё равно остаётся.

Есть ещё один простой способ изучить данные — построить графики связи переменных, сгруппированные по Sex.

Переменная Diameter выглядит наиболее линейной при Sex = I, но на этом всё.

# Создаём FacetGrid с диаграммами рассеяния
sns.lmplot(x="Diameter", y="Rings", hue="Sex", col="Sex", order=2, data=df)

Переменная Shell_weight, напротив, при больших значениях показывает слишком сильный разброс, и это искажает линейную зависимость.

# Создаём FacetGrid с диаграммами рассеяния
sns.lmplot(x="Shell_weight", y="Rings", hue="Sex", col="Sex", data=df)

Все эти наблюдения показывают, что задача линейной регрессии для этого набора данных очень непростая и модель, вероятнее всего, покажет слабые результаты. Но мы всё равно попробуем.

Кстати, я не припомню статей, где подробно разбирается, что именно пошло не так. Так что, проходя через это, мы тоже чему-то научимся.

Моделирование: используем Scikit-Learn

Давайте запустим модель из sklearn и оценим её по RMSE (корню из среднеквадратичной ошибки).

from sklearn.linear_model import LinearRegression
from sklearn.metrics import root_mean_squared_error

df2 = df.query('Height < 0.3 and Rings > 2 ').copy()
X = df2.drop(['Rings'], axis=1)
y = np.log(df2['Rings'])

lr = LinearRegression()
lr.fit(X, y)

predictions = lr.predict(X)

df2['Predictions'] = np.exp(predictions)
print(root_mean_squared_error(df2['Rings'], df2['Predictions']))

2.2383762717104916

Если посмотреть на первые строки таблицы (head), можно убедиться, что модели сложно даются оценки для больших значений целевой переменной (например, строки 0, 6, 7 и 9).

Шаг назад: пробуем другие трансформации

Итак, что мы можем сделать дальше?

Скорее всего, стоит удалить ещё часть выбросов и попробовать снова. Попробуем применить несупервизируемый алгоритм, чтобы найти дополнительные выбросы. Мы воспользуемся методом Local Outlier Factor и отбросим 5% наблюдений, признанных выбросами.

Заодно уберём мультиколлинеарность, удалив признаки Whole_weight и Length.

from sklearn.neighbors import LocalOutlierFactor
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline

# загрузим набор данных
abalone = fetch_ucirepo(id=1)

# данные (как pandas DataFrame)
X = abalone.data.features
y = abalone.data.targets

# One Hot Encode для признака Sex
ohe = OneHotEncoder(variables=['Sex'])
X = ohe.fit_transform(X)

# удалим Whole_weight и Length (мультиколлинеарность)
X.drop(['Whole_weight', 'Length'], axis=1, inplace=True)

# посмотрим на датафрейм
df = pd.concat([X,y], axis=1)

# создадим Pipeline для масштабирования данных и поиска выбросов методом Local Outlier Factor (на основе k ближайши�� соседей)
steps = [
('scale', StandardScaler()),
('LOF', LocalOutlierFactor(contamination=0.05))
]
# обучим модель и получим предсказания (метки выбросов)
outliers = Pipeline(steps).fit_predict(X)

# добавим столбец с метками выбросов
df['outliers'] = outliers

# моделирование
df2 = df.query('Height < 0.3 and Rings > 2 and outliers != -1').copy()
X = df2.drop(['Rings', 'outliers'], axis=1)
y = np.log(df2['Rings'])

lr = LinearRegression()
lr.fit(X, y)

predictions = lr.predict(X)

df2['Predictions'] = np.exp(predictions)
print(root_mean_squared_error(df2['Rings'], df2['Predictions']))

2.238174395913869

Результат тот же. Хм…

Ладно, мы можем продолжать экспериментировать с переменными и заниматься генерацией признаков (feature engineering) — и увидим небольшие улучшения. Например, если добавить квадраты Height, Diameter и Shell_weight. В сочетании с обработкой выбросов это снизит RMSE до 2.196.

# признаки второй степени
X['Diameter_2'] = X['Diameter'] ** 2
X['Height_2'] = X['Height'] ** 2
X['Shell_2'] = X['Shell_weight'] ** 2

Важно отметить: каждый новый признак, добавляемый в модель линейной регрессии, влияет на R² и иногда искусственно «раздувает» его, создавая ложное ощущение, что модель стала лучше, хотя это не так. В данном случае модель действительно улучшается, потому что мы добавляем в неё нелинейные компоненты за счёт признаков второй степени. Это можно подтвердить, посчитав скорректированный R²: он вырос с 0.495 до 0.517.

# Скорректированный R²
from sklearn.metrics import r2_score


r2 = r2_score(df2['Rings'], df2['Predictions'])
n= df2.shape[0]
p = df2.shape[1] - 1
adj_r2 = 1 - (1 - r2) * (n - 1) / (n - p - 1)
print(f'R²: {r2}')
print(f'Adjusted R²: {adj_r2}')

С другой стороны, если вернуть признаки Whole_weight и Length, метрики можно немного улучшить, но я бы не рекомендовал этого делать. В этом случае мы снова добавляем мультиколлинеарность и завышаем значимость коэффициентов некоторых переменных, что может привести к ошибкам оценок в будущем.

Моделирование: использование PyTorch

Итак. Теперь, когда у нас есть базовая модель, наша цель — построить линейную модель с использованием глубокого обучения и попробовать превзойти RMSE 2.196.

Для начала сразу оговорюсь: модели глубокого обучения лучше работают с масштабированными данными. Однако наши признаки X уже находятся на одной шкале, поэтому об этом можно не беспокоиться и двигаться дальше.

import torch
import torch.nn as nn
import torch.optim as optim
from torch.utils.data import DataLoader, TensorDataset

Нужно подготовить данные для моделирования в PyTorch. Здесь нам потребуется преобразовать данные в формат, который PyTorch «понимает», так как он не работает напрямую с pandas DataFrame.

• Возьмём тот же датафрейм, что и для базовой модели.

• Разделим X и Y.

• Преобразуем целевую переменную Y с помощью логарифма.

• Преобразуем обе части в массивы numpy, так как PyTorch не принимает DataFrame.

df2 = df.query('Height < 0.3 and Rings > 2 and outliers != -1').copy()
X = df2.drop(['Rings', 'outliers'], axis=1)
y = np.log(df2[['Rings']])


# Преобразуем X и Y в Numpy
X = X.to_numpy()
y = y.to_numpy()

Далее, используя TensorDataset, мы превращаем X и Y в тензоры и выводим пример.

# Подготовка с помощью TensorDataset
# TensorDataset помогает преобразовать набор данных в объект Tensor
dataset = TensorDataset(torch.tensor(X).float(), torch.tensor(y).float())


input_sample, label_sample = dataset[0]
print(f'** Input sample: {input_sample}, \n** Label sample: {label_sample}')

** Input sample: tensor([0.3650, 0.0950, 0.2245, 0.1010, 0.1500, 1.0000,
 0.0000, 0.0000, 0.1332, 0.0090, 0.0225]),
 ** Label sample: tensor([2.7081])

Затем, с помощью функции DataLoader, мы можем разбить данные на батчи. Это означает, что нейросеть будет обрабатывать данные порциями размера batch_size.

# Теперь воспользуемся DataLoader
batch_size = 500
dataloader = DataLoader(dataset, batch_size=batch_size, shuffle=True)

Модели в PyTorch лучше всего определять в виде классов.

• Класс наследуется от nn.Module — базового класса PyTorch для нейросетей.

• В методе init мы задаём слои модели. 

  • Вызов super().__init__() гарантирует, что класс будет вести себя как torch-объект.

• Метод forward описывает, что происходит с входными данными при прохождении через модель.

Здесь мы пропускаем их через линейные слои, определённые в init, и используем функции активации ReLU, чтобы добавить нелинейность при прямом проходе.

# 2. Создаём класс
class AbaloneModel(nn.Module):
  def __init__(self):
    super().__init__()
    self.linear1 = nn.Linear(in_features=X.shape[1], out_features=128)
    self.linear2 = nn.Linear(128, 64)
    self.linear3 = nn.Linear(64, 32)
    self.linear4 = nn.Linear(32, 1)


  def forward(self, x):
    x = self.linear1(x)
    x = nn.functional.relu(x)
    x = self.linear2(x)
    x = nn.functional.relu(x)
    x = self.linear3(x)
    x = nn.functional.relu(x)
    x = self.linear4(x)
    return x


# Создаём экземпляр модели
model = AbaloneModel()

Далее попробуем модель в деле, используя скрипт, который имитирует Random Search по гиперпараметрам.

• Создадим функцию ошибки для оценки модели.

• Создадим список для хранения данных о лучших вариантах модели и зададим best_loss большим числом, чтобы его можно было постепенно заменять более низкими значениями в ходе итераций.

• Зададим диапазон для шага обучения: будем использовать степени от −2 до −4 (то есть от 0.01 до 0.0001).

• Зададим диапазон для параметра momentum от 0.9 до 0.99.

• Получим данные.

• Обнулим градиенты, чтобы очистить результаты вычислений с предыдущих шагов.

• Обучим модель.

• Посчитаем функцию потерь и запишем показатели для лучшей модели.

• Посчитаем градиенты весов и смещений на обратном проходе.

• Повторим эти шаги N раз и в конце выведем лучшую модель.

# Среднеквадратичная ошибка (MSE) — стандартная функция потерь для задач регрессии
criterion = nn.MSELoss()


# Random Search
values = []
best_loss = 999
for idx in range(1000):
  # Случайным образом выбираем показатель степени для шага обучения в диапазоне от 2 до 4
  factor = np.random.uniform(2,5)
  lr = 10 ** -factor


  # Случайным образом выбираем momentum в диапазоне от 0.90 до 0.99
  momentum = np.random.uniform(0.90, 0.99)


  # 1. Получаем данные
  feature, target = dataset[:]
  # 2. Обнуляем градиенты, чтобы очистить результаты прошлых обратных проходов
  optimizer = optim.SGD(model.parameters(), lr=lr, momentum=momentum)
  optimizer.zero_grad()
  # 3. Прямой проход: считаем предсказание
  y_pred = model(feature)
  # 4. Считаем функцию потерь
  loss = criterion(y_pred, target)
  # 4.1. Регистрируем лучшую (минимальную) ошибку
  if loss < best_loss:
    best_loss = loss
    best_lr = lr
    best_momentum = momentum
    best_idx = idx


  # 5. Обратный проход: считаем градиенты функции потерь по W и b
  loss.backward()
  # 6. Обновляем параметры: корректируем W и b с учётом вычисленных градиентов
  optimizer.step()
  values.append([idx, lr, momentum, loss])


print(f'n: {idx},lr: {lr}, momentum: {momentum}, loss: {loss}')

n: 999,lr: 0.004782946959508322, momentum: 0.9801209929050066, loss: 0.06135804206132889

Когда мы нашли лучший шаг обучения и momentum, можно двигаться дальше.

# --- 3. Функция потерь и оптимизатор ---


# Среднеквадратичная ошибка (MSE) — стандартная функция потерь для задач регрессии
criterion = nn.MSELoss()


# Стохастический градиентный спуск (SGD) с небольшим шагом обучения (lr)
optimizer = optim.SGD(model.parameters(), lr=0.004, momentum=0.98)

Дальше мы переобучим модель, используя те же шаги, что и раньше, но теперь фиксируя шаг обучения и momentum.

Обучение модели в PyTorch требует более развёрнутого кода, чем простой вызов метода fit() в Scikit-Learn. Но это не проблема: структура почти всегда будет такой:

1. перевести модель в режим обучения model.train();

2. создать цикл по числу итераций (эпох), которое вы хотите выполнить;

3. в каждой итерации обнулять градиенты через optimizer.zero_grad();

4. получать батчи данных из dataloader;

5. считать предсказания y_pred = model(X);

6. вычислить функцию потерь loss = criterion(y_pred, target);

7. выполнить обратный проход loss.backward(), чтобы посчитать градиенты весов и смещений;

8. обновить веса и смещения с помощью optimizer.step().

Мы будем обучать эту модель в течение 1000 эпох (итераций). Дополнительно добавим шаг по сохранению лучшей версии модели (с минимальной ошибкой), чтобы в конце использовать именно её.

# 4. Обучение
torch.manual_seed(42)
NUM_EPOCHS = 1001
loss_history = []
best_loss = 999


# Переводим модель в режим обучения
model.train()


for epoch in range(NUM_EPOCHS):
  for data in dataloader:


    # 1. Получаем данные
    feature, target = data


    # 2. Обнуляем градиенты перед обратным проходом
    optimizer.zero_grad()


    # 3. Прямой проход: считаем предсказание
    y_pred = model(feature)


    # 4. Считаем функцию потерь
    loss = criterion(y_pred, target)
    loss_history.append(loss)


    # Сохраняем лучшую модель
    if loss < best_loss:
      best_loss = loss
      best_model_state = model.state_dict()  # сохраняем состояние лучшей модели


    # 5. Обратный проход: считаем градиенты функции потерь по W и b
    loss.backward()


    # 6. Обновляем параметры: корректируем W и b с учётом вычисленных градиентов
    optimizer.step()


    # Загружаем лучшее состояние модели перед тем, как использовать её для предсказаний
    model.load_state_dict(best_model_state)


  # Печатаем прогресс каждые 200 эпох
  if epoch % 200 == 0:
    print(epoch, loss.item())
    print(f'Best Loss: {best_loss}')

0 0.061786893755197525
Best Loss: 0.06033024191856384
200 0.036817338317632675
Best Loss: 0.03243456035852432
400 0.03307393565773964
Best Loss: 0.03077109158039093
600 0.032522525638341904
Best Loss: 0.030613820999860764
800 0.03488151729106903
Best Loss: 0.029514113441109657
1000 0.0369877889752388
Best Loss: 0.029514113441109657

Отлично. Модель обучена. Теперь пора её оценить.

Оценка качества

Давайте проверим, стала ли эта модель лучше по сравнению с обычной линейной регрессией. Для этого я переведу модель в режим оценки с помощью model.eval(), чтобы PyTorch понял, что нужно изменить поведение: из режима обучения перейти в режим инференса. Например, в этом режиме отключаются слои Dropout и batch normalization переводится в режим оценки.

# Получаем признаки и целевую переменную
features, targets = dataset[:]


# Получаем предсказания
model.eval()
with torch.no_grad():
  predictions = model(features)


# Добавляем предсказания в датафрейм
df2['Predictions'] = np.exp(predictions.detach().numpy())


# RMSE
print(root_mean_squared_error(df2['Rings'], df2['Predictions']))
2.1108551025390625

2.1108551025390625

Улучшение получилось скромным — примерно 4%.

Давайте посмотрим на некоторые предсказания обеих моделей.

Обе модели дают очень похожие результаты. Им становится сложнее справляться с предсказаниями по мере роста значения Rings. Причина — в «конусной» форме распределения целевой переменной.

Если немного об этом подумать:

• По мере увеличения числа колец (Rings) растёт и разброс, вносимый объясняющими переменными.

• Морское ушко с 15 кольцами может иметь гораздо более широкий диапазон значений признаков, чем экземпляр с 4 кольцами.

• Это сбивает модель с толку, потому что ей приходится проводить одну-единственную линию посередине данных, которые по сути не особенно линейны.

Заключение

В этом проекте мы разобрали много всего:

• Как исследовать данные.

• Как понять, подходит ли линейная модель для конкретной задачи.

• Как создать модель на PyTorch для задачи многомерной линейной регрессии.

В итоге мы увидели, что неоднородная целевая переменная, даже после степенных преобразований, может приводить к низкой эффективности модели. Наша модель всё равно лучше, чем просто предсказывать одно и то же среднее значение для всех наблюдений, но ошибка остаётся высокой — порядка 20 % от среднего значения.

Мы попытались использовать глубокое обучение, чтобы улучшить результат, но всей этой мощности оказалось недостаточно, чтобы заметно снизить ошибку. Я бы, вероятно, выбрал модель на Scikit-Learn: она проще и лучше интерпретируется.

Другой п��ть улучшения качества — собрать кастомный ансамбль, например Random Forest + линейная регрессия. Но эту часть я оставляю вам, если захотите поэкспериментировать.

Научиться работать с важнейшими моделями машинного обучения, NLP, DL, рекомендательными системами на практике с реальными данными можно на курсе "Machine Learning. Professional". Пройдите вступительный тест, чтобы узнать, подойдет ли вам программа курса.

А если не хватает «живых» кейсов и хочется увидеть, как ML ведёт себя в продакшене, приходите на демо-уроки:

• 11 декабря 20:00 — Ускоряем работу с видеопотоками при помощи TorchCodec. Записаться

• 17 декабря 18:00 — Оптимизируем построение модели через Pipeline. Записаться

• 22 декабря 18:00 — ML как основа современного AI. Записаться