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

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

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

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

Если данные явно изогнуты — параболическая траектория, насыщающаяся кривая роста — прямая их не опишет. Тогда аппроксимируют полиномом степени $m$: $y = c_0 + c_1 x + c_2 x^2 + \ldots + c_m \cdot x^m$. Задача та же — минимизировать сумму квадратов отклонений, — но неизвестных теперь $m+1$. Приравняв к нулю частные производные по каждому $c_k$, получаем систему из $m+1$ линейных уравнений — нормальные уравнения. Решаем её любым методом из раздела про СЛАУ (например, Гауссом) — и готово.

Элементы матрицы нормальных уравнений — это суммы степеней x: элемент $(i,j)$ равен $\sum x_k^{i+j}$, а правая часть $i$ — $\sum y_k \cdot 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^2 + 1$: $c_0 = 1, c_1 = 0, c_2 = 1$. Тот же код с другой степенью подгонит кубику, прямую или константу — достаточно поменять параметр степень.

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

Полиномы — лишь частный случай. МНК работает с любым набором базисных функций $\varphi_0(x), \varphi_1(x), \ldots$: модель $y = \sum c_k \cdot \varphi_k(x)$. Хотите подогнать $y = c_0 + c_1 \cdot \sin(x) + c_2 \cdot \cos(x)$? Подставьте в роли $\varphi$ синусы и косинусы — структура нормальных уравнений та же. Это называется линейной по параметрам моделью: нелинейность по x допускается, лишь бы коэффициенты входили линейно. Именно так строят гармонический анализ, регрессию на признаках и базисные разложения.

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

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

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

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

Итоги

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