Как аналитически решается задача линейной регрессии?

Как аналитически решается задача линейной регрессии?

Линейная регрессия решается аналитически через нормальное уравнение (Normal Equation), которое находит оптимальные коэффициенты в закрытом виде, без итеративных алгоритмов оптимизации.

Формулировка задачи

Модель:

ŷ = β₀ + β₁x₁ + β₂x₂ + ... + βₖxₖ = Xβ

Матричная форма:

ŷ = Xβ

где:

• X — матрица признаков (n × (k+1), с колонкой единиц для β₀)
• β — вектор коэффициентов (k+1 × 1)
• n — количество примеров
• k — количество признаков

Целевая функция (Loss)

Минимизируем квадратичную ошибку (OLS - Ordinary Least Squares):

L(β) = ||y - Xβ||² = Σ(yᵢ - ŷᵢ)²

Аналитическое решение

Шаг 1: Берём производную по β

∂L/∂β = ∂/∂β (y - Xβ)ᵀ(y - Xβ)
      = ∂/∂β (yᵀy - 2βᵀXᵀy + βᵀXᵀXβ)
      = -2Xᵀy + 2XᵀXβ

Шаг 2: Приравниваем к нулю

-2Xᵀy + 2XᵀXβ = 0
XᵀXβ = Xᵀy

Шаг 3: Нормальное уравнение

β = (XᵀX)⁻¹Xᵀy

Это замкнутое решение — находим оптимальные коэффициенты за один шаг!

Практическая реализация

Вариант 1: С нуля (NumPy)

import numpy as np
from sklearn.datasets import make_regression

# Генерируем данные
X, y = make_regression(n_samples=100, n_features=5, noise=10, random_state=42)

# Добавляем колонку единиц для свободного члена (intercept)
X = np.column_stack([np.ones(X.shape[0]), X])

# Нормальное уравнение: β = (XᵀX)⁻¹Xᵀy
beta = np.linalg.inv(X.T @ X) @ X.T @ y

print(f"Коэффициенты (β): {beta}")
print(f"Свободный член (β₀): {beta[0]:.4f}")
print(f"Коэффициенты признаков: {beta[1:]}")

# Предсказания
y_pred = X @ beta

# Оценка качества
from sklearn.metrics import r2_score, mean_squared_error
r2 = r2_score(y, y_pred)
mse = mean_squared_error(y, y_pred)
print(f"\nR² = {r2:.4f}")
print(f"MSE = {mse:.4f}")

Вариант 2: Через NumPy lstsq (более численно устойчиво)

# np.linalg.lstsq решает XB = Y напрямую
# Более numerically stable, чем обращение матрицы

beta, residuals, rank, s = np.linalg.lstsq(X, y, rcond=None)
print(f"Коэффициенты: {beta}")
print(f"Остатки: {residuals}")
print(f"Ранг матрицы: {rank}")

Вариант 3: Через scikit-learn

from sklearn.linear_model import LinearRegression
from sklearn.datasets import make_regression

X, y = make_regression(n_samples=100, n_features=5, noise=10, random_state=42)

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

print(f"Коэффициенты: {model.coef_}")
print(f"Свободный член: {model.intercept_}")
print(f"R² на train: {model.score(X, y):.4f}")

Матричный расчёт (пошагово)

# Пример с 2 признаками
X = np.array([
    [1, 1, 2],      # Пример 1: intercept=1, x1=1, x2=2
    [1, 2, 3],      # Пример 2: intercept=1, x1=2, x2=3
    [1, 3, 4],      # Пример 3: intercept=1, x1=3, x2=4
])
y = np.array([5, 8, 11])  # Ответы

print("X матрица:")
print(X)
print("\ny вектор:")
print(y)

# Шаг 1: XᵀX
XTX = X.T @ X
print("\nXᵀX:")
print(XTX)

# Шаг 2: Xᵀy
XTy = X.T @ y
print("\nXᵀy:")
print(XTy)

# Шаг 3: (XᵀX)⁻¹
XTX_inv = np.linalg.inv(XTX)
print("\n(XᵀX)⁻¹:")
print(XTX_inv)

# Шаг 4: β = (XᵀX)⁻¹Xᵀy
beta = XTX_inv @ XTy
print("\nβ (коэффициенты):")
print(beta)
print(f"\nРешение: y = {beta[0]:.4f} + {beta[1]:.4f}*x1 + {beta[2]:.4f}*x2")

Сравнение методов

Численная нестабильность

# Проблема: если XᵀX близка к сингулярной матрице
# (мультиколлинеарность), то (XᵀX)⁻¹ может быть нестабильна

# Пример плохой обусловленности
X_bad = np.array([
    [1, 100, 100.01],
    [1, 200, 200.02],
    [1, 300, 300.03],
])
y = np.array([1, 2, 3])

# XᵀX будет почти сингулярной
XTX = X_bad.T @ X_bad
cond_number = np.linalg.cond(XTX)
print(f"Число обусловленности: {cond_number:.2e}")
print(f"Это плохо! (>1000 → численная нестабильность)")

# Решение: Использовать lstsq вместо inv
beta_lstsq = np.linalg.lstsq(X_bad, y, rcond=None)[0]
beta_inv = np.linalg.inv(XTX) @ X_bad.T @ y

print(f"\nТочное (lstsq): {beta_lstsq}")
print(f"Через inv: {beta_inv}")
print(f"Разница: {np.abs(beta_lstsq - beta_inv)}")

Геометрическая интерпретация

Векторное пространство:
- y — целевой вектор
- Xβ — проекция y на подпространство, натянутое столбцами X
- ŷ = Xβ — ближайшая точка в подпространстве
- (y - ŷ) — ортогональна Xβ

Для оптимума: Xᵀ(y - Xβ) = 0
              Xᵀy = XᵀXβ
              β = (XᵀX)⁻¹Xᵀy

Регуляризация (Ridge, Lasso)

# Ridge: добавляем штраф за большие коэффициенты
# β_ridge = (XᵀX + λI)⁻¹Xᵀy

lambda_reg = 0.1
I = np.eye(X.shape[1])
beta_ridge = np.linalg.inv(X.T @ X + lambda_reg * I) @ X.T @ y

print(f"OLS: {beta}")
print(f"Ridge (λ=0.1): {beta_ridge}")
print(f"Ridge коэффициенты меньше (регуляризация)")

Когда аналитическое решение работает лучше

• Маленькие k (< 1000) — аналитическое решение быстрее
• Полный батч — когда используются все данные сразу
• Нужна точность — аналитическое решение точнее градиента
• Нет необходимости в online learning — не нужно обновлять по частям

Аналитическое решение — золотой стандарт для линейной регрессии!