Руководство по PyTorch: Множественная регрессия с нуля

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

В этой статье акцент делается не на обработке естественного языка, а на задаче множественной линейной регрессии с двумя целями: продемонстрировать процесс построения модели в 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, поскольку это единственная категориальная переменная.

# Data Load 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 # fetch dataset abalone = fetch_ucirepo(id=1) # data (as pandas dataframes) X = abalone.data.features y = abalone.data.targets # One Hot Encode Sex ohe = OneHotEncoder(variables=['Sex']) X = ohe.fit_transform(X) # View df = pd.concat([X,y], axis=1)

Вот структура набора данных.

Чтобы построить эффективную модель, проведем анализ данных.

Анализ данных

Первоначальные шаги при изучении набора включают:

1. Оценку распределения целевой переменной.

# Looking at our Target variable plt.hist(y) plt.title('Rings [Target Variable] Distribution');

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

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

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

# Statistical description df.describe()

Далее оценим корреляции.

# Looking at the correlations (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, выявляют несколько проблем:

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

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

Дополнительный анализ — построение графиков для проверки связей переменных, сгруппированных по Sex.

Переменная Diameter проявляет наиболее линейную связь при Sex=I, но это единственный заметный случай.

# Create a FacetGrid with scatterplots sns.lmplot(x="Diameter", y="Rings", hue="Sex", col="Sex", order=2, data=df);

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

# Create a FacetGrid with scatterplots sns.lmplot(x="Shell_weight", y="Rings", hue="Sex", col="Sex", data=df);

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

Построение модели: Использование Scikit-Learn

Запустим модель из Scikit-Learn и оценим ее с помощью корня из среднеквадратичной ошибки.

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

Анализ заголовка подтверждает, что модель испытывает трудности с предсказаниями для высоких значений (например, строки 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 # fetch dataset abalone = fetch_ucirepo(id=1) # data (as pandas dataframes) X = abalone.data.features y = abalone.data.targets # One Hot Encode Sex ohe = OneHotEncoder(variables=['Sex']) X = ohe.fit_transform(X) # Drop Whole Weight and Length (multicolinearity) X.drop(['Whole_weight', 'Length'], axis=1, inplace=True) # View df = pd.concat([X,y], axis=1) # Let's create a Pipeline to scale the data and find outliers using KNN Classifier steps = [ ('scale', StandardScaler()), ('LOF', LocalOutlierFactor(contamination=0.05)) ] # Fit and predict outliers = Pipeline(steps).fit_predict(X) # Add column df['outliers'] = outliers # Modeling 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.

# Second Order Variables 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.

# Adjusted 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. Фреймворк не принимает обычные DataFrame, поэтому внесем корректировки.

• Используем тот же DataFrame из базовой модели.
• Разделим на X и Y.
• Преобразуем Y в логарифм.
• Переведем в numpy-массивы.

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 and Y to Numpy X = X.to_numpy() y = y.to_numpy()

С помощью TensorDataset преобразуем X и Y в тензоры и выведем пример.

# Prepare with TensorData # TensorData helps us transforming the dataset to Tensor object dataset = TensorDataset(torch.tensor(X).float(), torch.tensor(y).float()) input_sample, label_sample = dataset[0] print(f'** Input sample: {input_sample}, 
** 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.

# Next, let's use 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. Creating a class 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 # Instantiate model model = AbaloneModel()

Протестируем модель с помощью скрипта, имитирующего случайный поиск.

• Определим критерий ошибки для оценки.
• Создадим список для хранения данных лучшей модели и инициализируем best_loss высоким значением.
• Зададим диапазон learning rate от -2 до -4 (от 0.01 до 0.0001).
• Диапазон momentum от 0.9 до 0.99.
• Получим данные.
• Обнулим градиенты для очистки от предыдущих итераций.
• Обучим модель.
• Вычислим ошибку и зафиксируем параметры лучшей модели.
• Выполним обратное распространение для весов и смещений.
• Проведем N итераций и выведем лучшую модель.

# Mean Squared Error (MSE) is standard for regression criterion = nn.MSELoss() # Random Search values = [] best_loss = 999 for idx in range(1000): # Randomly sample a learning rate factor between 2 and 4 factor = np.random.uniform(2,5) lr = 10 ** -factor # Randomly select a momentum between 0.85 and 0.99 momentum = np.random.uniform(0.90, 0.99) # 1. Get Data feature, target = dataset[:] # 2. Zero Gradients: Clear old gradients before the backward pass optimizer = optim.SGD(model.parameters(), lr=lr, momentum=momentum) optimizer.zero_grad() # 3. Forward Pass: Compute prediction y_pred = model(feature) # 4. Compute Loss loss = criterion(y_pred, target) # 4.1 Register best Loss if loss < best_loss: best_loss = loss best_lr = lr best_momentum = momentum best_idx = idx # 5. Backward Pass: Compute gradient of the loss w.r.t W and b' loss.backward() # 6. Update Parameters: Adjust W and b using the calculated gradients 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

Получив оптимальные learning rate и momentum, продолжим.

# --- 3. Loss Function and Optimizer --- # Mean Squared Error (MSE) is standard for regression criterion = nn.MSELoss() # Stochastic Gradient Descent (SGD) with a small learning rate (lr) optimizer = optim.SGD(model.parameters(), lr=0.004, momentum=0.98)

Переобучим модель, следуя тем же шагам, но с фиксированными learning rate и momentum.

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

1. Активировать режим model.train().
2. Создать цикл для заданного числа итераций (эпох).
3. Обнулить градиенты предыдущих проходов с optimizer.zero_grad().
4. Получить батчи из dataloader.
5. Вычислить предсказания с model(X).
6. Рассчитать ошибку с criterion(y_pred, target).
7. Выполнить обратное распространение для весов и смещений: loss.backward().
8. Обновить веса и смещения: optimizer.step().

Обучим модель на 1000 эпохах. Добавим шаг для сохранения лучшей модели по ошибке.

# 4. Training torch.manual_seed(42) NUM_EPOCHS = 1001 loss_history = [] best_loss = 999 # Put model in training mode model.train() for epoch in range(NUM_EPOCHS): for data in dataloader: # 1. Get Data feature, target = data # 2. Zero Gradients: Clear old gradients before the backward pass optimizer.zero_grad() # 3. Forward Pass: Compute prediction y_pred = model(feature) # 4. Compute Loss loss = criterion(y_pred, target) loss_history.append(loss) # Get Best Model if loss < best_loss: best_loss = loss best_model_state = model.state_dict() # save best model # 5. Backward Pass: Compute gradient of the loss w.r.t W and b' loss.backward() # 6. Update Parameters: Adjust W and b using the calculated gradients optimizer.step() # Load the best model before returning predictions model.load_state_dict(best_model_state) # Print status every 50 epochs 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.

# Get features features, targets = dataset[:] # Get Predictions model.eval() with torch.no_grad(): predictions = model(features) # Add to dataframe df2['Predictions'] = np.exp(predictions.detach().numpy()) # RMSE print(root_mean_squared_error(df2['Rings'], df2['Predictions']))

2.1108551025390625

Улучшение скромное, около 4%.

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

Обе модели дают схожие результаты, испытывая большие трудности при росте числа колец. Это обусловлено конической формой распределения целевой переменной.

Размышляя об этом:

• С увеличением колец растет дисперсия от объясняющих переменных.
• Абалон с 15 кольцами попадает в более широкий диапазон значений, чем с 4 кольцами.
• Модель вынуждена проводить единую линию посреди нелинейных данных.

Итоги

В проекте рассмотрены ключевые аспекты:

• Методы анализа данных.
• Проверка пригодности линейной модели.
• Создание модели множественной линейной регрессии в PyTorch.

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

Глубокое обучение не дало значительного снижения ошибки. Предпочтительнее модель Scikit-Learn за простоту и интерпретируемость.

Для улучшения можно создать ансамбль из Random Forest и линейной регрессии, но это оставляем на усмотрение читателя.