Полиномиальная аппроксимация и нормальные уравнения

Урок про обобщение МНК на полиномы и произвольные базисы через систему нормальных уравнений.

Нормальные уравнения — линейная система AᵀA·c = Aᵀy, дающая коэффициенты c МНК-модели; для полинома степени m это система размера (m+1)×(m+1).

Когда прямой мало

Если данные явно изогнуты — параболическая траектория, насыщающаяся кривая роста — прямая их не опишет. Тогда аппроксимируют полиномом степени m: y = c₀ + c₁x + c₂x² + ... + c_m·x^m. Задача та же — минимизировать сумму квадратов отклонений, — но неизвестных теперь m+1. Приравняв к нулю частные производные по каждому c_k, получаем систему из m+1 линейных уравнений — нормальные уравнения. Решаем её любым методом из раздела про СЛАУ (например, Гауссом) — и готово.

Элементы матрицы нормальных уравнений — это суммы степеней x: элемент (i,j) равен Σ x_k^{i+j}, а правая часть i — Σ y_k·x_k^i. Звучит громоздко, но в коде это два вложенных цикла.

def гаусс(A, b):
    n = len(A)
    M = [row[:] + [b[i]] for i, row in enumerate(A)]
    for k in range(n):
        p = max(range(k, n), key=lambda i: abs(M[i][k]))
        M[k], M[p] = M[p], M[k]
        for i in range(k+1, n):
            f = M[i][k] / M[k][k]
            for j in range(k, n+1):
                M[i][j] -= f * M[k][j]
    x = [0.0]*n
    for i in range(n-1, -1, -1):
        x[i] = (M[i][n] - sum(M[i][j]*x[j] for j in range(i+1, n))) / M[i][i]
    return x

def полином_мнк(xs, ys, степень):
    n = степень + 1
    # матрица нормальных уравнений: A[i][j] = sum x^(i+j)
    A = [[sum(x**(i+j) for x in xs) for j in range(n)] for i in range(n)]
    b = [sum(ys[k]*xs[k]**i for k in range(len(xs))) for i in range(n)]
    return гаусс(A, b)

# Данные точно на параболе y = x^2 + 1 -> ждём коэффициенты [1, 0, 1]
xs = [0, 1, 2, 3, 4]
ys = [1, 2, 5, 10, 17]
коэф = полином_мнк(xs, ys, 2)
print("коэффициенты [c0, c1, c2]:", [round(v, 4) for v in коэф])
print(f"модель: y = {коэф[0]:.2f} + {коэф[1]:.2f}·x + {коэф[2]:.2f}·x²")

Вывод:

коэффициенты [c0, c1, c2]: [1.0, -0.0, 1.0]
модель: y = 1.00 + -0.00·x + 1.00·x²

МНК точно восстановил y = x² + 1: c₀=1, c₁=0, c₂=1. Тот же код с другой степенью подгонит кубику, прямую или константу — достаточно поменять параметр степень.

Обобщение: любой набор базисных функций

Полиномы — лишь частный случай. МНК работает с любым набором базисных функций φ₀(x), φ₁(x), ...: модель y = Σ c_k·φ_k(x). Хотите подогнать y = c₀ + c₁·sin(x) + c₂·cos(x)? Подставьте в роли φ синусы и косинусы — структура нормальных уравнений та же. Это называется линейной по параметрам моделью: нелинейность по x допускается, лишь бы коэффициенты входили линейно. Именно так строят гармонический анализ, регрессию на признаках и базисные разложения.

Как работает под капотом

Прямое решение нормальных уравнений AᵀA·c = Aᵀy просто, но имеет коварство: матрица AᵀA для полиномов высокой степени чудовищно плохо обусловлена (близка к матрице Гильберта). Уже при степени 6–8 на равномерных данных число обусловленности взлетает, и коэффициенты «плывут». Поэтому в промышленных пакетах нормальные уравнения не решают напрямую: используют QR-разложение или SVD матрицы A, которые решают задачу МНК устойчиво, без возведения обусловленности в квадрат. А ещё применяют ортогональные полиномы (Чебышёва, Лежандра) вместо степеней 1, x, x² — у них матрица почти диагональна. В библиотеке: numpy.polyfit(x, y, deg) внутри использует именно устойчивое разложение, а не наивные нормальные уравнения.

Частые ошибки

Брать слишком высокую степень. Полином высокой степени переобучается (повторяет шум) и страдает от плохой обусловленности нормальных уравнений; обычно степень 2–4 достаточно.
Решать нормальные уравнения напрямую при высокой степени. AᵀA возводит число обусловленности в квадрат; для устойчивости берите QR/SVD или ортогональные базисы.
Экстраполировать полиномом. За пределами данных полином высокой степени улетает в бесконечность — прогнозам вне диапазона доверять нельзя.

Итоги

Полиномиальный МНК сводится к нормальным уравнениям AᵀA·c = Aᵀy размера (m+1)×(m+1).
Их решают методами из раздела про СЛАУ; элементы матрицы — суммы степеней x.
МНК работает с любым базисом (полиномы, тригонометрия), лишь бы модель была линейна по параметрам.
Прямые нормальные уравнения плохо обусловлены при высокой степени — на практике берут QR/SVD или ортогональные полиномы.