Python: Реализовать линейную регрессию с нуля

Решение

1. Аналитическое решение (Нормальное уравнение)

import numpy as np
import matplotlib.pyplot as plt

class LinearRegressionNormalEquation:
    """Линейная регрессия через нормальное уравнение
    
    w = (X^T X)^(-1) X^T y
    """
    
    def __init__(self):
        self.w = None  # веса
        self.b = None  # смещение
    
    def fit(self, X, y):
        """Обучение модели
        
        Args:
            X: (n_samples, n_features) - признаки
            y: (n_samples,) - целевая переменная
        """
        n_samples, n_features = X.shape
        
        # Добавляем столбец единиц для смещения (bias)
        X_with_bias = np.column_stack([np.ones(n_samples), X])
        
        # Нормальное уравнение: w = (X^T X)^(-1) X^T y
        # w[0] = bias, w[1:] = коэффициенты
        w_full = np.linalg.inv(X_with_bias.T @ X_with_bias) @ X_with_bias.T @ y
        
        self.b = w_full[0]
        self.w = w_full[1:]
        
        return self
    
    def predict(self, X):
        """Предсказание"""
        return X @ self.w + self.b
    
    def mse(self, y_true, y_pred):
        """Mean Squared Error"""
        return np.mean((y_true - y_pred) ** 2)

# Пример
print("=== НОРМАЛЬНОЕ УРАВНЕНИЕ ===")

# Синтетические данные
np.random.seed(42)
n_samples = 100
n_features = 3

X = np.random.randn(n_samples, n_features)
y = 2 * X[:, 0] + 3 * X[:, 1] - X[:, 2] + 5 + np.random.randn(n_samples) * 0.1

print(f"Размер данных: X {X.shape}, y {y.shape}")
print(f"Истинные коэффициенты: [2, 3, -1], смещение: 5")

# Обучение
model_ne = LinearRegressionNormalEquation()
model_ne.fit(X, y)

print(f"\nОбученные коэффициенты: {model_ne.w}")
print(f"Обученное смещение: {model_ne.b}")

# Предсказание
y_pred = model_ne.predict(X)
mse = model_ne.mse(y, y_pred)
print(f"MSE на обучающем наборе: {mse:.6f}")

print(f"""
=== ОБЪЯСНЕНИЕ ===

Нормальное уравнение: w = (X^T X)^(-1) X^T y

Математика:
1. Ищем минимум функции потерь: L = ||Xw - y||^2
2. Берем производную по w и приравниваем к нулю
3. dL/dw = 2X^T(Xw - y) = 0
4. X^T Xw = X^T y
5. w = (X^T X)^(-1) X^T y

Преимущества:
- Аналитическое решение (точно за O(n^3))
- Не нужна настройка learning rate
- Гарантированная сходимость

Недостатки:
- O(n^3) временная сложность (медленно для больших n)
- Нужна инверсия матрицы (может быть численно неустойчива)
- Требует много памяти
""")

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

class LinearRegressionGradientDescent:
    """Линейная регрессия через градиентный спуск
    
    w := w - alpha * grad(L)
    где grad(L) = 2/n * X^T(Xw - y)
    """
    
    def __init__(self, learning_rate=0.01, n_iterations=1000, verbose=False):
        self.learning_rate = learning_rate
        self.n_iterations = n_iterations
        self.verbose = verbose
        self.w = None
        self.b = None
        self.loss_history = []
    
    def fit(self, X, y):
        """Обучение модели"""
        n_samples, n_features = X.shape
        
        # Инициализируем веса нулями (или случайными)
        self.w = np.zeros(n_features)
        self.b = 0
        
        # Градиентный спуск
        for iteration in range(self.n_iterations):
            # Предсказание
            y_pred = X @ self.w + self.b
            
            # Ошибка
            error = y_pred - y
            
            # Градиенты
            dw = (2 / n_samples) * (X.T @ error)
            db = (2 / n_samples) * np.sum(error)
            
            # Обновление весов
            self.w -= self.learning_rate * dw
            self.b -= self.learning_rate * db
            
            # Сохраняем функцию потерь
            mse = np.mean(error ** 2)
            self.loss_history.append(mse)
            
            if self.verbose and (iteration + 1) % 100 == 0:
                print(f"Iteration {iteration + 1}/{self.n_iterations}, Loss: {mse:.6f}")
        
        return self
    
    def predict(self, X):
        """Предсказание"""
        return X @ self.w + self.b
    
    def mse(self, y_true, y_pred):
        """Mean Squared Error"""
        return np.mean((y_true - y_pred) ** 2)

# Пример
print("\n=== ГРАДИЕНТНЫЙ СПУСК ===")

model_gd = LinearRegressionGradientDescent(
    learning_rate=0.01,
    n_iterations=1000,
    verbose=True
)

model_gd.fit(X, y)

print(f"\nОбученные коэффициенты: {model_gd.w}")
print(f"Обученное смещение: {model_gd.b}")

y_pred_gd = model_gd.predict(X)
mse_gd = model_gd.mse(y, y_pred_gd)
print(f"MSE на обучающем наборе: {mse_gd:.6f}")

# Визуализация сходимости
plt.figure(figsize=(10, 5))
plt.plot(model_gd.loss_history)
plt.xlabel('Iteration')
plt.ylabel('Loss (MSE)')
plt.title('Gradient Descent Convergence')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('gradient_descent_convergence.png')
print(f"\nГрафик сходимости сохранен")

print(f"""
=== ОБЪЯСНЕНИЕ ГРАДИЕНТНОГО СПУСКА ===

1. Инициализируем веса (обычно нулями)
2. На каждой итерации:
   a) Предсказываем: y_pred = X @ w + b
   b) Вычисляем ошибку: error = y_pred - y
   c) Вычисляем градиенты:
      dw = 2/n * X^T @ error
      db = 2/n * sum(error)
   d) Обновляем веса:
      w := w - learning_rate * dw
      b := b - learning_rate * db
3. Повторяем шаг 2 до сходимости

Преимущества:
- O(n) на итерацию (быстро для больших данных)
- Легко параллелизируется
- Работает для больших датасетов

Недостатки:
- Нужно выбрать learning rate и количество итераций
- Может зависнуть или дивергировать
- Сходимость медленнее чем нормальное уравнение
""")

3. Варианты градиентного спуска

class LinearRegressionSGD:
    """Stochastic Gradient Descent - обновление весов по одному примеру"""
    
    def __init__(self, learning_rate=0.01, n_epochs=10, batch_size=1):
        self.learning_rate = learning_rate
        self.n_epochs = n_epochs
        self.batch_size = batch_size
        self.w = None
        self.b = None
    
    def fit(self, X, y):
        n_samples, n_features = X.shape
        self.w = np.zeros(n_features)
        self.b = 0
        
        for epoch in range(self.n_epochs):
            # Перемешиваем данные
            indices = np.random.permutation(n_samples)
            X_shuffled = X[indices]
            y_shuffled = y[indices]
            
            # Минибатчи
            for i in range(0, n_samples, self.batch_size):
                X_batch = X_shuffled[i:i+self.batch_size]
                y_batch = y_shuffled[i:i+self.batch_size]
                
                y_pred = X_batch @ self.w + self.b
                error = y_pred - y_batch
                
                dw = (2 / len(X_batch)) * (X_batch.T @ error)
                db = (2 / len(X_batch)) * np.sum(error)
                
                self.w -= self.learning_rate * dw
                self.b -= self.learning_rate * db
        
        return self
    
    def predict(self, X):
        return X @ self.w + self.b

print("\n=== STOCHASTIC GRADIENT DESCENT ===")

model_sgd = LinearRegressionSGD(
    learning_rate=0.01,
    n_epochs=100,
    batch_size=10
)

model_sgd.fit(X, y)
y_pred_sgd = model_sgd.predict(X)
mse_sgd = np.mean((y - y_pred_sgd) ** 2)
print(f"MSE (SGD): {mse_sgd:.6f}")

4. Сравнение подходов

print("\n=== СРАВНЕНИЕ МЕТОДОВ ===")

# Сравнение предсказаний
comparison = pd.DataFrame({
    'Method': ['Normal Equation', 'Gradient Descent', 'SGD'],
    'Coefficients': [model_ne.w, model_gd.w, model_sgd.w],
    'Bias': [model_ne.b, model_gd.b, model_sgd.b],
    'MSE': [mse, mse_gd, mse_sgd]
})

print(comparison)

print(f"""
=== СРАВНЕНИЕ ХАРАКТЕРИСТИК ===

                 | Normal Eq | Grad Desc | SGD
             Time| O(n^3)    | O(n*iter)| O(n*iter)
          Memory | O(n^2)    | O(n)     | O(batch)
       Stability | Численная | Хорошая  | Хорошая
    Big Data     | Плохо     | Хорошо   | Отлично
    Learning Rate| -         | Нужна    | Нужна
  Convergence    | Точная    | Прибли.. | Прибли..

=== КОГДА ИСПОЛЬЗОВАТЬ ===

✓ Normal Equation:
  - Маленькие датасеты (< 10k примеров)
  - Когда нужна точность
  - Когда нет требований к скорости

✓ Batch Gradient Descent:
  - Средние датасеты
  - Когда нужен контроль сходимости
  - GPU training

✓ SGD / Mini-batch SGD:
  - Большие датасеты (> 100k примеров)
  - Online learning
  - Когда памяти не хватает на весь датасет
  - Production системы
""")

# Визуализация сравнения
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Сравнение функции потерь
axes[0].plot(model_gd.loss_history, label='GD', linewidth=2)
axes[0].set_xlabel('Iteration')
axes[0].set_ylabel('Loss')
axes[0].set_title('Convergence Comparison')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
axes[0].set_yscale('log')

# Сравнение коэффициентов
coeff_comparison = np.array([
    model_ne.w,
    model_gd.w,
    model_sgd.w
])

x_pos = np.arange(n_features)
width = 0.25

axes[1].bar(x_pos - width, coeff_comparison[0], width, label='Normal Eq')
axes[1].bar(x_pos, coeff_comparison[1], width, label='GD')
axes[1].bar(x_pos + width, coeff_comparison[2], width, label='SGD')
axes[1].set_ylabel('Coefficient Value')
axes[1].set_title('Learned Coefficients Comparison')
axes[1].set_xticks(x_pos)
axes[1].set_xticklabels([f'w{i}' for i in range(n_features)])
axes[1].legend()
axes[1].grid(True, alpha=0.3, axis='y')

plt.tight_layout()
plt.savefig('methods_comparison.png')
print(f"\nСравнение визуализировано в methods_comparison.png")