Почему numpy.polyfit и SVD не сходятся в задаче линейных наименьших квадратов?

Линейные наименьшие квадраты (МНК) – фундаментальный метод в статистике и машинном обучении, используемый для нахождения наилучшей линейной зависимости между переменными. numpy.polyfit и сингулярное разложение (SVD) – два подхода для решения задач МНК, но они могут давать разные результаты. Понимание причин этих расхождений критически важно для правильного выбора и применения метода.

Постановка задачи линейных наименьших квадратов

Задача МНК состоит в минимизации суммы квадратов разностей между наблюдаемыми значениями и значениями, предсказанными линейной моделью. Математически, дана матрица A (матрица признаков) и вектор b (вектор целевых значений), требуется найти вектор x, который минимизирует ||Ax - b||₂², где ||.||₂ обозначает евклидову норму.

Решение с использованием numpy.polyfit: полиномиальная регрессия

numpy.polyfit используется для полиномиальной регрессии. Он предполагает, что зависимость между переменными может быть аппроксимирована полиномом заданной степени. Функция принимает на вход массивы x и y (значения независимой и зависимой переменных) и степень полинома deg. Она возвращает коэффициенты полинома, минимизирующего сумму квадратов ошибок.

import numpy as np
from typing import Tuple


def polyfit_regression(x: np.ndarray, y: np.ndarray, degree: int) -> np.ndarray:
    """Выполняет полиномиальную регрессию с использованием numpy.polyfit."""
    coefficients: np.ndarray = np.polyfit(x, y, degree)
    return coefficients

# Пример использования
x_data: np.ndarray = np.array([1, 2, 3, 4, 5])
y_data: np.ndarray = np.array([2, 4, 5, 4, 5])
degree: int = 2
coefficients: np.ndarray = polyfit_regression(x_data, y_data, degree)
print(f"Коэффициенты полинома: {coefficients}")

Решение с использованием SVD: общее решение наименьших квадратов

SVD – мощный инструмент линейной алгебры, позволяющий разложить матрицу A на три матрицы: U, Σ, и Vᵀ, где U и V – унитарные матрицы, а Σ – диагональная матрица с сингулярными числами на диагонали. Решение задачи МНК с использованием SVD находится как x = V Σ⁺ Uᵀ b, где Σ⁺ – псевдообратная матрица к Σ.

import numpy as np
from numpy.linalg import svd, solve
from typing import Tuple


def svd_least_squares(A: np.ndarray, b: np.ndarray) -> np.ndarray:
    """Решает задачу наименьших квадратов с использованием SVD."""
    U, s, V = svd(A)
    # Создаем диагональную матрицу из сингулярных чисел
    Sigma = np.zeros((A.shape[0], A.shape[1]))
    Sigma[:A.shape[1], :A.shape[1]] = np.diag(s)
    # Вычисляем псевдообратную матрицу Sigma
    Sigma_inv = np.linalg.pinv(Sigma)
    # Решение x = V * Sigma_inv * U^T * b
    x = V @ Sigma_inv @ U.T @ b
    return x

# Пример использования
A: np.ndarray = np.array([[1, 1], [1, 2], [1, 3]])
b: np.ndarray = np.array([2, 4, 5])
x: np.ndarray = svd_least_squares(A, b)
print(f"Решение, полученное с помощью SVD: {x}")

Причины расхождения результатов numpy.polyfit и SVD

Различные цели и подходы: полиномиальная регрессия vs. общая аппроксимация

numpy.polyfit предназначен для полиномиальной регрессии, в то время как SVD предоставляет общее решение задачи наименьших квадратов. polyfit строит матрицу Вандермонда неявно, в то время как SVD работает с произвольной матрицей A. Если матрица A не соответствует структуре полиномиальной регрессии (например, содержит категориальные признаки), результаты будут различаться.

Масштабирование данных: влияние на SVD и numpy.polyfit

Масштабирование данных может существенно влиять на результаты SVD. Большие значения в матрице A могут привести к доминированию одних сингулярных чисел над другими, что повлияет на решение. numpy.polyfit менее чувствителен к масштабированию, поскольку он внутренне обрабатывает данные для полиномиальной регрессии.

Особенности реализации numpy.polyfit и возможные проблемы с обусловленностью матрицы

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

Реклама

Учет ранга матрицы при SVD и его влияние на решение

SVD позволяет явно учитывать ранг матрицы A. Если ранг матрицы меньше числа столбцов, это означает, что система уравнений является недоопределенной. SVD позволяет найти решение с минимальной нормой. numpy.polyfit не предоставляет такой явный контроль над рангом матрицы.

Практический пример расхождения и способы диагностики

Генерация синтетических данных для демонстрации расхождения

import numpy as np
import matplotlib.pyplot as plt


np.random.seed(42)
x = np.linspace(0, 10, 50)
y = 2*x + 1 + np.random.normal(0, 2, 50)

# Добавим нелинейность
y[25:] += (x[25:] - 5)**2

A = np.vstack([x, np.ones(len(x))]).T

# polyfit
coeffs_polyfit = np.polyfit(x, y, 1)
y_polyfit = np.poly1d(coeffs_polyfit)(x)

# SVD
U, s, V = np.linalg.svd(A)
s_inv = np.linalg.pinv(np.diag(s))
x_svd = V.T @ s_inv @ U.T @ y
y_svd = A @ x_svd

plt.scatter(x, y, label='Данные')
plt.plot(x, y_polyfit, color='red', label=f'Polyfit (наклон={coeffs_polyfit[0]:.2f}, сдвиг={coeffs_polyfit[1]:.2f})')
plt.plot(x, y_svd, color='green', label=f'SVD (наклон={x_svd[0]:.2f}, сдвиг={x_svd[1]:.2f})')
plt.legend()
plt.xlabel('x')
plt.ylabel('y')
plt.title('Сравнение Polyfit и SVD')
plt.show()

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

Визуализация результатов, полученных numpy.polyfit и SVD

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

Анализ сингулярных чисел для оценки обусловленности матрицы

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

Выбор оптимального метода решения в зависимости от задачи

Если задача заключается в аппроксимации данных полиномом, следует использовать numpy.polyfit. Если требуется общее решение задачи наименьших квадратов, особенно для плохо обусловленных матриц, следует использовать SVD.

Методы сглаживания расхождений

Предварительное масштабирование данных для улучшения результатов SVD и numpy.polyfit

Масштабирование данных (например, с использованием StandardScaler из scikit-learn) может улучшить обусловленность матрицы и повысить точность решений, полученных с помощью SVD и numpy.polyfit.

Регуляризация в SVD для повышения устойчивости решения (например, Tikhonov regularization)

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

Использование других функций numpy для полиномиальной аппроксимации (например, numpy.polynomial)

numpy.polynomial предоставляет более гибкие инструменты для полиномиальной аппроксимации, чем numpy.polyfit, включая возможность работы с различными базисами полиномов и управления обусловленностью.

Заключение

Ключевые выводы и рекомендации по выбору метода решения

numpy.polyfit и SVD – полезные инструменты для решения задач линейных наименьших квадратов, но они имеют разные цели и подходы. numpy.polyfit предназначен для полиномиальной регрессии, а SVD – для общего решения задачи наименьших квадратов. Выбор метода зависит от конкретной задачи и свойств данных.

Ограничения numpy.polyfit и SVD при решении задач линейных наименьших квадратов

numpy.polyfit может быть неустойчивым при высокой степени полинома. SVD может быть чувствителен к масштабированию данных и требовать регуляризации для плохо обусловленных матриц.

Перспективы развития и альтернативные подходы

Развитие алгоритмов решения задач наименьших квадратов продолжается, появляются новые методы, устойчивые к шуму и хорошо работающие с большими объемами данных. Примерами могут служить итеративные методы, методы случайного проецирования и методы на основе градиентного спуска.


Добавить комментарий