Солверы гладких выпуклых функций как ключ к пониманию Data Science

Солверы гладких выпуклых функций как ключ к пониманию Data Science: Реверс-инжиниринг scikit-learn

Зачем Data Scientist'у нужны солверы?

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

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

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

Почему солверы — это основа Data Science

Солверы — это не просто математические алгоритмы; это фундаментальные инструменты, которые позволяют нам извлекать знания из данных и строить мощные предсказательные модели. Понимание их принципов работы, преимуществ и недостатков дает вам не только технические навыки, но и глубокое интуитивное понимание того, как машины «учатся».

Начиная изучение Data Science с солверов, вы закладываете прочный фундамент. Вы учитесь не просто применять готовые библиотеки, но и понимать, что происходит под капотом. Это позволяет вам:

1. Выбирать правильный инструмент: зная особенности каждого солвера, вы можете осознанно выбирать наиболее подходящий для вашей конкретной задачи и набора данных.
2. Отлаживать модели: если модель не сходится или ведет себя странно, понимание работы солвера поможет вам диагностировать проблему.
3. Оптимизировать производительность: вы сможете настраивать параметры солвера (например, скорость обучения) для достижения лучшей производимости и более быстрой сходимости.
4. Разрабатывать новые алгоритмы: глубокое понимание основ открывает двери для создания собственных оптимизационных решений.

Как найти минимум функции?

Наша основная задача в Data Science часто сводится к минимизации функции потерь. Функция потерь (или целевая функция) измеряет, насколько плохо наша модель справляется с задачей. Чем меньше значение функции потерь, тем лучше модель. Для задач регрессии, как уже упоминалось, часто используется MSE. Математически это выглядит так:

MSE = (1/n) * Σ(y_i - ŷ_i)^2

Где:

• n — количество наблюдений.
• y_i — истинное значение для i-го наблюдения.
• ŷ_i — предсказанное значение нашей моделью для i-го наблюдения.

Предположим, наша модель — это простая линейная регрессия: ŷ = wX + b, где w — веса, X — признаки, b — смещение. Наша задача — найти такие w и b, которые минимизируют MSE. Это и есть та самая «низкая точка», которую мы ищем.

Роадмап по миру солверов

Теперь давайте погрузимся в мир солверов и рассмотрим, как каждый из них помогает нам найти оптимальные параметры модели.

1. Solver: Градиентный спуск (Gradient Descent) и Стохастический градиентный спуск (SGD)

Принцип работы: представьте, что вы стоите на склоне горы и хотите спуститься в долину. Самый простой способ — сделать шаг в направлении наибольшего уклона вниз. Градиентный спуск делает именно это: он итеративно корректирует параметры модели, двигаясь в направлении, противоположном градиенту функции потерь (наибольшему подъему). Скорость, с которой мы делаем шаги, называется скоростью обучения (learning rate).

• Градиентный спуск (GD): на каждой итерации он использует все данные для вычисления градиента. Это обеспечивает точное направление, но может быть очень медленным для больших наборов данных.
• Стохастический градиентный спуск (SGD): вместо всех данных, SGD использует только один случайный пример (или небольшой батч примеров) для вычисления градиента на каждой итерации. Это делает его намного быстрее, особенно для больших данных, но путь к минимуму может быть более «шумным» и извилистым.

Когда использовать: GD подходит для небольших наборов данных, где важна точность. SGD — это рабочий конь для большинства задач глубокого обучения и больших наборов данных, где скорость важнее идеальной точности на каждой итерации.

Кейс: обучение простой линейной регрессии на большом датасете.

Пример на Python (SGD):

import numpy as np

# Генерируем синтетические данные
np.random.seed(0)
X = 2 * np.random.rand(100, 1)
y = 4 + 3 * X + np.random.randn(100, 1)

# Добавляем столбец для смещения (intercept)
X_b = np.c_[np.ones((100, 1)), X]

# Параметры SGD
learning_rate = 0.01
n_iterations = 1000
m = len(X_b)

# Инициализируем случайные веса
theta = np.random.randn(2, 1)

for iteration in range(n_iterations):
    random_index = np.random.randint(m)
    xi = X_b[random_index:random_index+1]
    yi = y[random_index:random_index+1]
    gradients = 2 * xi.T.dot(xi.dot(theta) - yi)
    theta = theta - learning_rate * gradients

print("Оптимальные параметры (SGD):", theta.ravel())

2. Solver: L-BFGS-B (Limited-memory Broyden–Fletcher–Goldfarb–Shanno)

Принцип работы: L-BFGS-B — это квазиньютоновский метод. В отличие от градиентного спуска, который использует только информацию о первом порядке (градиент), квазиньютоновские методы пытаются аппроксимировать информацию о втором порядке (гессиан — матрица вторых производных). Это позволяет им делать более интеллектуальные шаги к минимуму. L-BFGS-B особенно хорош тем, что он использует ограниченный объем памяти, что делает его применимым для задач с большим количеством параметров.

Когда использовать: L-BFGS-B часто используется в задачах, где функция потерь сложна, но градиент может быть вычислен. Он быстрее сходится, чем градиентный спуск, и часто является хорошим выбором для задач с умеренным количеством параметров, где требуется высокая точность.

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

Пример на Python (с использованием scipy.optimize.minimize):

from scipy.optimize import minimize
import numpy as np

# Генерируем синтетические данные для линейной регрессии
np.random.seed(0)
X = 2 * np.random.rand(100, 1)
y = 4 + 3 * X + np.random.randn(100, 1)

# Добавляем столбец для смещения
X_b = np.c_[np.ones((100, 1)), X]

# Функция потерь (MSE) для L-BFGS-B
def mse_loss(theta, X, y):
    predictions = X.dot(theta)
    return np.mean((predictions - y)**2)

# Градиент MSE
def mse_gradient(theta, X, y):
    predictions = X.dot(theta)
    return 2 * X.T.dot(predictions - y) / len(y)

# Начальные параметры
theta_initial = np.random.randn(2)

# Оптимизация с L-BFGS-B
result = minimize(
    mse_loss,
    theta_initial,
    args=(X_b, y),
    method='L-BFGS-B',
    jac=mse_gradient,
)

print("Оптимальные параметры (L-BFGS-B):", result.x)

3. Solver: Нормальное уравнение (Normal Equation) и его разложения (Cholesky, LU, QR, SVD)

Принцип работы: для линейной регрессии с функцией потерь MSE существует аналитическое решение, которое позволяет найти оптимальные параметры w и b без итераций. Это решение называется нормальным уравнением:

theta = (X^T X)^-1 X^T y

Где theta — это вектор оптимальных параметров, X — матрица признаков, y — вектор целевых значений. Вычисление (X^T X)^-1 может быть вычислительно затратным и численно нестабильным для больших или плохо обусловленных матриц. Поэтому часто используются различные методы разложения матрицы для более эффективного и стабильного решения.

• Разложение Холецкого (Cholesky Decomposition): применяется, когда матрица X^T X является симметричной и положительно определенной (что часто бывает в линейной регрессии). Оно разлагает X^T X на произведение нижней треугольной матрицы L и ее транспонированной L^T (X^T X = L L^T). Это позволяет решить систему уравнений более эффективно.
• LU-разложение (LU Decomposition): разлагает матрицу X^T X на произведение нижней треугольной матрицы L и верхней треугольной матрицы U (X^T X = L U). Также используется для эффективного решения систем линейных уравнений.
• QR-разложение (QR Decomposition) и LSQR (Golub-Kahan): QR-разложение разлагает матрицу X на ортогональную матрицу Q и верхнюю треугольную матрицу R (X = QR). Это позволяет переформулировать нормальное уравнение и решить его более стабильно, особенно когда X^T X плохо обусловлена. LSQR — это итерационный метод, основанный на QR-разложении, который особенно эффективен для больших и разреженных систем.
• Разложение по сингулярным значениям (Singular Value Decomposition, SVD): это мощный метод разложения матрицы X на три матрицы: U, Σ (сигма, диагональная матрица сингулярных значений) и V^T (X = U Σ V^T). SVD является наиболее численно стабильным способом решения нормального уравнения, особенно когда матрица X^T X вырождена или близка к вырожденной (т.е. имеет мультиколлинеарность). Он также позволяет работать с псевдообратной матрицей, что полезно в случаях, когда X^T X не имеет обратной.

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

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

Пример на Python (Нормальное уравнение и SVD):

import numpy as np
from numpy.linalg import inv, svd

# Генерируем синтетические данные
np.random.seed(0)
X = 2 * np.random.rand(100, 1)
y = 4 + 3 * X + np.random.randn(100, 1)

# Добавляем столбец для смещения
X_b = np.c_[np.ones((100, 1)), X]

# 1. Решение с помощью Нормального уравнения
theta_normal = inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y)
print("Оптимальные параметры (Нормальное уравнение):", theta_normal.ravel())

# 2. Решение с помощью SVD
U, s, Vt = svd(X_b, full_matrices=False)
s_inv = np.diag(1 / s)
theta_svd = Vt.T.dot(s_inv).dot(U.T).dot(y)
print("Оптимальные параметры (SVD):", theta_svd.ravel())

4. Solver: Сопряженные градиенты (Conjugate Gradient, Sparse_cg)

Принцип работы: метод сопряженных градиентов — это итерационный алгоритм, предназначенный для решения больших систем линейных уравнений, особенно когда матрица системы симметрична и положительно определена. Он не требует хранения всей матрицы в памяти, что делает его очень эффективным для разреженных матриц (матриц с большим количеством нулевых элементов). Sparse_cg — это реализация сопряженных градиентов, оптимизированная для разреженных систем.

В отличие от градиентного спуска, который может двигаться по извилистой траектории, сопряженные градиенты выбирают направления, которые «ортогональны» предыдущим градиентам в определенном смысле, что позволяет им быстрее сходиться к минимуму.

Когда использовать: идеально подходит для очень больших систем линейных уравнений, возникающих, например, при обучении линейных моделей на разреженных данных (много нулей), таких как в обработке естественного языка или рекомендательных системах. Когда нормальное уравнение становится слишком дорогим для вычисления из-за размера матрицы X^T X.

Кейс: решение задачи линейной регрессии с огромным количеством признаков, где матрица признаков X разрежена.

Пример на Python (с использованием scipy.sparse.linalg.cg):

import numpy as np
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import cg

# Генерируем синтетические разреженные данные
np.random.seed(0)
n_samples, n_features = 1000, 10000  # Много признаков
X_dense = np.random.rand(n_samples, n_features)
X_dense[np.random.rand(*X_dense.shape) < 0.99] = 0  # Делаем разреженной
X_sparse = csr_matrix(X_dense)
y = X_dense.dot(np.random.rand(n_features)) + np.random.randn(n_samples) * 0.1

# Для сопряженных градиентов решаем систему A x = b, где A = X.T @ X, b = X.T @ y
A = X_sparse.T @ X_sparse
b = X_sparse.T @ y

# Решаем с помощью сопряженных градиентов
theta_cg, info = cg(A, b, atol=1e-6)

print("Оптимальные параметры (Conjugate Gradient, первые 5):", theta_cg[:5])
print("Статус сходимости (0 = успешно):", info)

5. Solver: LSQR (Golub-Kahan) и QR-разложение (Givens, Gram-Schmidt)

Принцип работы: эти методы также направлены на решение систем линейных уравнений, но с акцентом на численную стабильность и эффективность, особенно когда матрица X может быть плохо обусловлена или прямоугольной (больше строк, чем столбцов, что типично для регрессии).

• QR-разложение (Gram-Schmidt, Givens): как уже упоминалось, X = QR. Это разложение позволяет преобразовать исходную систему X theta = y в QR theta = y, а затем R theta = Q^T y. Поскольку R является верхней треугольной матрицей, эту систему легко решить обратной подстановкой. Методы Гивенса и Грама-Шмидта — это разные алгоритмы для выполнения QR-разложения. Гивенс более численно стабилен, чем Грама-Шмидта.
• LSQR (Golub-Kahan): это итерационный метод, который использует бидиагонализацию Ланцоша (связанную с QR-разложением) для решения систем наименьших квадратов. Он особенно эффективен для больших, разреженных систем и обладает хорошей численной стабильностью. LSQR не требует явного формирования X^T X, что позволяет избежать проблем с обусловленностью, которые могут возникнуть при прямом использовании нормального уравнения.

Когда использовать: QR-разложение и LSQR — отличный выбор, когда вам нужно решить задачу линейных наименьших квадратов (OLS/SSE) для больших или плохо обусловленных матриц, где нормальное уравнение может быть нестабильным или слишком медленным. LSQR предпочтителен для очень больших и разреженных систем.

Кейс: линейная регрессия на больших наборах данных с потенциальной мультиколлинеарностью или когда требуется высокая численная стабильность.

Пример на Python (LSQR):

import numpy as np
from scipy.sparse.linalg import lsqr

# Генерируем синтетические данные
np.random.seed(0)
X = 2 * np.random.rand(1000, 50)  # Много наблюдений, умеренно признаков
y = X.dot(np.random.rand(50)) + np.random.randn(1000) * 0.5

# Решаем с помощью LSQR
theta_lsqr, istop, itn, r1norm, r2norm, anorm, acond, arnorm, xnorm, var = lsqr(X, y)

print("Оптимальные параметры (LSQR, первые 5):", theta_lsqr[:5])
print("Статус сходимости (1 = успешно):", istop)

6. Solver: Stochastic Average Gradient (SAG) & SAGA

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

• SAG: хранит градиенты для каждого примера и на каждой итерации обновляет один случайный градиент, а затем использует среднее всех градиентов для обновления параметров. Это позволяет ему сходиться быстрее, чем SGD, и с лучшей асимптотической скоростью.
• SAGA: улучшенная версия SAG, которая также хранит градиенты, но использует немного другую схему обновления, что делает ее более гибкой и часто более быстрой на практике. SAGA также имеет теоретические гарантии сходимости, аналогичные SAG.

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

Кейс: обучение линейных моделей (регрессия, классификация) на очень больших наборах данных, где SGD слишком шумен, а GD слишком медленен.

Пример на Python (концептуальный, так как прямых реализаций в scipy нет, но есть в sklearn для некоторых моделей):

# Пример использования SAGA в scikit-learn для логистической регрессии
from sklearn.linear_model import LogisticRegression
from sklearn.datasets import make_classification

# Генерируем синтетические данные для классификации
X, y = make_classification(n_samples=10000, n_features=20, random_state=42)

# Обучаем логистическую регрессию с солвером SAGA
# solver= поддерживает L1 и L2 регуляризацию и хорошо масштабируется
model = LogisticRegression(solver='saga', max_iter=1000, random_state=42)
model.fit(X, y)

print("Коэффициенты модели (SAGA, первые 5):", model.coef_[0][:5])
print("Смещение (SAGA):", model.intercept_[0])

7. Solver: Scalar-form Normal Equation

Принцип работы: это не отдельный солвер в том же смысле, что и другие, а скорее способ решения нормального уравнения для очень простых случаев, например, для одномерной линейной регрессии (y = wx + b). В этом случае можно вывести формулы для w и b в скалярной форме, без использования матричной алгебры. Это по сути упрощенная версия нормального уравнения.

Когда использовать: только для очень простых случаев, когда у вас всего один признак и вы хотите быстро получить аналитическое решение без использования сложных библиотек. В реальных задачах Data Science это редко применяется.

Кейс: простейшая линейная регрессия с одним признаком для демонстрационных целей.

Пример на Python:

import numpy as np

# Генерируем синтетические данные для одномерной регрессии
np.random.seed(0)
X = 2 * np.random.rand(100, 1)
y = 4 + 3 * X + np.random.randn(100, 1)

# Вычисляем параметры по формулам для одномерной регрессии
# w = (sum(xi*yi) - n*mean(x)*mean(y)) / (sum(xi^2) - n*mean(x)^2)
# b = mean(y) - w*mean(x)
x_mean = np.mean(X)
y_mean = np.mean(y)

# Числитель для w
numerator = np.sum(X * y) - len(X) * x_mean * y_mean

# Знаменатель для w
denominator = np.sum(X**2) - len(X) * x_mean**2

w = numerator / denominator
b = y_mean - w * x_mean

print("Оптимальные параметры (Scalar-form Normal Equation):", w[0], b[0])

Выбор солвера: Какой выбрать для вашей задачи?

Выбор правильного солвера — это искусство, которое приходит с опытом, но есть несколько общих рекомендаций:

В этой статье мы рассмотрели широкий спектр солверов, от простых и интуитивно понятных до сложных и высокоэффективных. Мы увидели, как они применяются для минимизации функции потерь MSE/SSE в задачах линейной регрессии, и как каждый из них имеет свою нишу применения. Надеемся, что это путешествие по миру солверов вдохновило вас на дальнейшее изучение и применение этих мощных инструментов в вашей практике Data Science.