Ctrl+K

Критерии переобучения и регуляризация

Критерии переобучения и регуляризация#

(Здесь и далее перезапускайте ячейку перед взаимодействием с графиками. Используются одинаковые глобальные переменные)

Вспомним основную проблему - переобучение. Характерный пример:

from bokeh.io import output_notebook
from bokeh.resources import INLINE
output_notebook(INLINE)
Loading BokehJS ...
Hide code cell source 
import numpy as np
from sklearn.preprocessing import PolynomialFeatures
from bokeh.layouts import column, row
from bokeh.models import Slider, ColumnDataSource, Div, CustomJS
from bokeh.plotting import figure, show, output_file

# Генерация данных
np.random.seed(42)
X = np.random.rand(50, 1) * 2 - 1  # Данные в диапазоне [-1, 1]
y = -4 + 3 * X - 2 * X**2 + 5 * X**3 + np.random.randn(50, 1) * 0.5  # Квадратичная зависимость с шумом

# Подготовка данных для Bokeh
X_flat = X.flatten()
y_flat = y.flatten()

# Источник данных
source_data = ColumnDataSource(data=dict(x=X_flat, y=y_flat))
source_pred = ColumnDataSource(data=dict(x=X_flat, y=np.zeros_like(y_flat)))
source_hist = ColumnDataSource(data=dict(top=[], left=[], right=[]))  # Данные для гистограммы

# График данных и предсказаний
p = figure(y_range=(-14, 2), x_range=(-1, 1), width=400, height=400, title="Polynomial Overfitting Visualization")
p.scatter('x', 'y', source=source_data, color="red", legend_label="Data")
pred_line = p.line('x', 'y', source=source_pred, color="blue", legend_label="Prediction")
p.legend.location = "bottom_right"

# График гистограммы весов
hist_fig = figure(width=400, height=400, title="Weights Histogram")
hist_fig.quad(top='top', bottom=0, left='left', right='right', source=source_hist, fill_color="green", line_color="black")

# Ползунок для степени полинома
slider = Slider(start=1, end=49, value=1, step=1, title="Degree of Polynomial")

# Текст для отображения MSE
mse_text = Div(text=f"MSE: 0")

# CustomJS для обновления предсказаний и пересчета MSE на клиентской стороне
callback = CustomJS(args=dict(source_data=source_data, source_pred=source_pred, source_hist=source_hist, slider=slider, mse_text=mse_text), code="""
    const degree = slider.value;  // Степень полинома
    let x = source_data.data['x'].slice();  // Копируем массив x
    let y_true = source_data.data['y'].slice();  // Копируем массив y
    const n = x.length;
    
    // Сортируем x и соответствующие y_true по возрастанию x
    const indices = Array.from(x.keys()).sort((a, b) => x[a] - x[b]);
    x = indices.map(i => x[i]);
    y_true = indices.map(i => y_true[i]);

    // Генерация более плотной сетки для предсказаний
    const grid_points = 200;  // Количество точек для более плотной сетки
    const x_min = Math.min(...x);
    const x_max = Math.max(...x);
    const x_dense = [];
    const step = (x_max - x_min) / (grid_points - 1);
    
    for (let i = 0; i < grid_points; i++) {
        x_dense.push(x_min + i * step);  // Заполняем массив точками на сетке
    }

    // Матрица X для полиномиальной регрессии на исходных данных
    const X_poly = [];
    for (let i = 0; i < n; i++) {
        const row = [];
        for (let d = 0; d <= degree; d++) {
            row.push(Math.pow(x[i], d));
        }
        X_poly.push(row);
    }
    
    // Простое решение линейной регрессии
    const XT = X_poly[0].map((_, colIndex) => X_poly.map(row => row[colIndex]));
    const XTX = XT.map(row => row.map((_, colIndex) => row.reduce((acc, curr, i) => acc + curr * X_poly[i][colIndex], 0)));
    const XTy = XT.map(row => row.reduce((acc, curr, i) => acc + curr * y_true[i], 0));
    
    function solve(A, b) {
        const n = A.length;
        const x = Array(n).fill(0);
        for (let i = 0; i < n; i++) {
            let maxEl = Math.abs(A[i][i]), maxRow = i;
            for (let k = i + 1; k < n; k++) {
                if (Math.abs(A[k][i]) > maxEl) {
                    maxEl = Math.abs(A[k][i]);
                    maxRow = k;
                }
            }
            [A[maxRow], A[i]] = [A[i], A[maxRow]];
            [b[maxRow], b[i]] = [b[i], b[maxRow]];
            for (let k = i + 1; k < n; k++) {
                const c = -A[k][i] / A[i][i];
                for (let j = i; j < n; j++) {
                    A[k][j] += c * A[i][j];
                }
                b[k] += c * b[i];
            }
        }
        for (let i = n - 1; i >= 0; i--) {
            x[i] = b[i] / A[i][i];
            for (let k = i - 1; k >= 0; k--) {
                b[k] -= A[k][i] * x[i];
            }
        }
        return x;
    }
    
    const coefs = solve(XTX, XTy);
    
    // Предсказания на плотной сетке
    const y_dense_pred = [];
    for (let i = 0; i < grid_points; i++) {
        let pred = 0;
        for (let d = 0; d <= degree; d++) {
            pred += coefs[d] * Math.pow(x_dense[i], d);
        }
        y_dense_pred.push(pred);
    }

    // Обновление предсказаний на графике (плотная сетка)
    source_pred.data['x'] = x_dense;  // Заменяем на плотную сетку
    source_pred.data['y'] = y_dense_pred;  // Обновляем предсказанные y
    source_pred.change.emit();
    
    // Вычисление MSE на исходных данных
    let mse = 0;
    for (let i = 0; i < n; i++) {
        let pred = 0;
        for (let d = 0; d <= degree; d++) {
            pred += coefs[d] * Math.pow(x[i], d);
        }
        mse += Math.pow(y_true[i] - pred, 2);
    }
    mse /= n;
    
    // Обновление текста с MSE
    mse_text.text = `MSE: ${mse.toFixed(5)}`;
    
    // Обновление гистограммы весов
    const top = coefs;  // Значения коэффициентов
    const left = Array.from({length: coefs.length}, (_, i) => i - 0.4);  // Индексы коэффициентов для left
    const right = Array.from({length: coefs.length}, (_, i) => i + 0.4);  // Индексы коэффициентов для right

    source_hist.data['top'] = top;
    source_hist.data['left'] = left;
    source_hist.data['right'] = right;
    source_hist.change.emit();
""")

# Привязка ползунка к callback
slider.js_on_change('value', callback)

# Компоновка
layout = column(row(p, hist_fig), slider, mse_text)

# Визуализация
show(layout)

Видим характерные признаки переобучения:

  1. Потеря предсказательной способности на тестовой выборке - а зачем нам такая модель?

  2. «Взрыв» весов - при большой сложности модели \(m\) почти все веса становятся очень большими

Возникает идея, что корень проблемы лежит именно в странном поведении весов модели при увеличении её сложности. Действительно, кажется маловероятным, что в какой-либо задаче регрессии мы будем получать коэффициенты порядка \(10^{8}\) и больше. Значит, логично давать модели штраф за «сложность» (величину коэффициентов).

На этом и основан метод борьбы с переобучением - регуляризация.

Регуляризация наивно#

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

\[ \mathcal{L} \rightarrow \mathcal{L}+R(\mathbf w) \]

Такой штраф должен ограничивать количество признаков (чтобы их не стало больше количества точек в train-е) или рост весов модели \(\mathbf w\) (чтобы функция не могла быстро осциллировать).

L2-регуляризация#

Таким образом, самой простой функцией регуляризации является просто квадрат нормы вектора весов модели \(w\), также называемая L2-регуляризация (ridge).

\[ R(\mathbf w)=\tau\|\mathbf w\|^2=\tau \cdot \mathbf w^{T} \mathbf w, \quad \tau > 0 \]

Здесь параметр \(\tau\) отвечает «силе» регуляризации. Чем он больше, теми сильнее мы «штрафуем» модель за большие признаки.

Таким образом, модифицированный MSE, например, будет выглядеть как

\[ \mathcal{L}=\frac{1}{N}\|X w-Y\|^2+\tau\|w\|^2 \]

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

\[ \frac{\partial \mathcal{L}}{\partial w}=2 X^{T}(X w-Y)+2 \tau w=0 \]
\[ 2\left(X^{T} X+\tau \cdot \mathbb{1}\right) w=2 X^{T} Y \]
\[ w_{best}=\left(X^{T} X+\tau \cdot \mathbb{1}\right)^{-1} X^{T} Y \]

Видим, что при росте \(\tau\) в обращаемой матрице растёт диагональ, что приводит к изменению её числа обусловленности. Как мы помним, число обусловленности напрямую влияет на эффективность и погрешность любых численных методов и чем оно меньше, тем лучше.

\[ \mu(X^TX)=\sqrt{\frac{\lambda_{\text {max }}}{\lambda_{\text {min }}}} \rightarrow \sqrt{\frac{\lambda_{\text {max }}+\tau}{\lambda_{\text {min }}+\tau}} \]

Таким образом, чем меньше это выражениe, тем проще и точнее происходит обращение матрицы. Тем не менее, наша основная цель не эффективность обращения матрицы, а хорошая модель! Это значит, что параметр \(\tau\) надо выбирать не из условия минимизации числа обусловленности, а из условия минимизации Loss-a.

Говорят, что \(\tau\) - гиперпараметр (внешний параметр), т.к. он не изменяется при обучении модели (процесса нахождения наилучшего \(\mathbf w\)). Они напрямую влияют на качество обучения и итоговую модель. Как подбирать наилучшие значения гиперпараметров обсудим в конце семинара.

L1-регуляризация#

Иным выбором функции регуляризации может служить L1-регуляризация (lasso). Вместо квадрата модуля, мы возьмём модуль вектора непосредственно.

\[ R(\mathbf w)=\mu\|\mathbf w\|_1=\mu \sum_{\alpha=1}^m\left|\text w_\alpha\right|, \quad \mu > 0 \]

Таким образом, новая функция ошибок будет выглядеть как

\[ \mathcal{L}=\frac{1}{N}\left\|X\mathbf w-Y\right\|^2+\mu \sum_\alpha\left|\text w_\alpha\right| \]

Отметим, что при данной регуляризации, функция ошибок недиффиренциируема в нуле и аналитическое решение получить можно, но крайне сложно и никто так не делает. Необходимо численно её минимизировать - как это делать обсудим позже в этом семинаре.

Тем не менее, данная регуляризация имеет существенное преимущество - вес признака \(w_k\) может в определённый момент обнулиться, что отвечает неважности этого признака для модели. Т.е. данная регуляризация приводит к отбору признаков. Это позволяет избавляться от линейно-зависимых или бессмысленных признаков.

Чем больше гиперпараметр \(mu\), тем больше признаков будет зануляться. При этом по мере роста \(mu\) самые неинформативные признаки будут зануляться первыми. Этой процедурой можно отбирать самые важные признаки (коэффициенты модели).

ElasticNet#

Иногда применяют комбинированный способ - сумма ridge и lasso. Такой регуляризатор называется ElasticNet. В среднем работает лучше.

\[ \mathcal{L} \rightarrow \mathcal{L}+\tau\|\mathbf w\|^2+\mu \sum_\alpha\left|\text w_\alpha\right| \]

Есть и другие регуляризации (см. документацию sklearn и лекцию).

Важность нормировки#

Отдельно отметим важность нормировки данных. Во всех приведённых регуляризаторах веса признаков одинаковы (\(\tau\) и \(\mu\)) - таким образом нам важно нормировать входные данные для модели (и иксы, и игреки) - т.е. вычитать среднее и масштабировать к единичной дисперсии, чтобы регуляризатор «одинаково относился» ко всем признакам. Естественно при таком раскладе, мы как бы «шифруем» данные для модели и при получении предсказаний их надо обратно «дешифровать».

\[ x_i^{\prime}=\left(x_i-\bar{x}\right) / \sigma_x \]

Реализация и визуализация#

Реализация на python с помощью sklearn:

import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import Ridge # Lasso
from sklearn.metrics import mean_squared_error, r2_score

# Генерация данных
np.random.seed(42)
X = np.random.rand(50, 1) * 2 - 1  # Данные в диапазоне [-1, 1]
y = -4 + 3 * X - 2 * X**2 + 5 * X**3 + np.random.randn(50, 1) * 0.5  # Кубическая зависимость с шумом

# Параметры
degrees = [1, 3, 5, 7]  # Степени полинома для проверки
alpha = 0.1  # Параметр регуляризации

plt.figure(figsize=(14, 10))

for i, degree in enumerate(degrees):
    # Создание полиномиальных признаков
    poly_features = PolynomialFeatures(degree=degree, include_bias=False)
    X_poly = poly_features.fit_transform(X)
    
    # Создание и обучение Ridge модели
    ridge_model = Ridge(alpha=alpha)
    ridge_model.fit(X_poly, y.ravel())
    
    # Предсказание на более плотной сетке для построения графика
    X_test = np.linspace(-1, 1, 300).reshape(-1, 1)
    X_test_poly = poly_features.transform(X_test)
    y_pred = ridge_model.predict(X_test_poly)
    
    # Вычисление метрик
    y_train_pred = ridge_model.predict(X_poly)
    mse = mean_squared_error(y, y_train_pred)
    r2 = r2_score(y, y_train_pred)
    
    # Построение графика
    plt.subplot(2, 2, i+1)
    plt.scatter(X, y, color='blue', label='Данные')
    plt.plot(X_test, y_pred, color='red', label='Модель')
    plt.title(f'Степень полинома = {degree}\nMSE = {mse:.4f}, R2 = {r2:.4f}')
    plt.xlabel('X')
    plt.ylabel('y')
    plt.legend()

plt.tight_layout()
plt.show()

# Вывод коэффициентов для лучшей модели (используем степень 3, так как знаем истинную зависимость)
best_poly_features = PolynomialFeatures(degree=3, include_bias=False)
X_best_poly = best_poly_features.fit_transform(X)
best_ridge_model = Ridge(alpha=alpha)
best_ridge_model.fit(X_best_poly, y.ravel())

coefficients = best_ridge_model.coef_
intercept = best_ridge_model.intercept_

print("Коэффициенты полинома степени 3:")
print(f"y = {intercept:.4f} + {coefficients[0]:.4f}x + {coefficients[1]:.4f}x^2 + {coefficients[2]:.4f}x^3")
../../_images/46d5a988ec2eae87f0119265e847118f69fbbad36961ea1f53a974609d41127a.png 
Коэффициенты полинома степени 3:
y = -4.0113 + 3.2281x + -1.9761x^2 + 4.4991x^3

Визуализация Ridge регуляризации#

Hide code cell source 
import numpy as np
from sklearn.preprocessing import PolynomialFeatures
from bokeh.layouts import column, row
from bokeh.models import Slider, ColumnDataSource, Div, CustomJS
from bokeh.plotting import figure, show, output_file

# Генерация данных
np.random.seed(42)
X = np.random.rand(50, 1) * 2 - 1  # Данные в диапазоне [-1, 1]
y = -4 + 3 * X - 2 * X**2 + 5 * X**3 + np.random.randn(50, 1) * 0.5  # Квадратичная зависимость с шумом

# Подготовка данных для Bokeh
X_flat = X.flatten()
y_flat = y.flatten()

# Источник данных
source_data = ColumnDataSource(data=dict(x=X_flat, y=y_flat))
source_pred = ColumnDataSource(data=dict(x=X_flat, y=np.zeros_like(y_flat)))
source_hist = ColumnDataSource(data=dict(top=[], left=[], right=[]))  # Данные для гистограммы

# График данных и предсказаний
p = figure(y_range=(-14, 2), x_range=(-1, 1), width=400, height=400, title="Polynomial Fitting with Ridge Regularization")
p.scatter('x', 'y', source=source_data, color="red", legend_label="Data")
pred_line = p.line('x', 'y', source=source_pred, color="blue", legend_label="Prediction")
p.legend.location = "bottom_right"

# График гистограммы весов
hist_fig = figure(width=400, height=400, title="Weights Histogram")
hist_fig.quad(top='top', bottom=0, left='left', right='right', source=source_hist, fill_color="green", line_color="black")

# Ползунки для степени полинома и параметра регуляризации
degree_slider = Slider(start=1, end=49, value=1, step=1, title="Degree of Polynomial")
tau_slider = Slider(start=0, end=50, value=0, step=0.1, title="Regularization Parameter (tau)")

# Текст для отображения MSE
mse_text = Div(text=f"MSE: 0")

# CustomJS для обновления предсказаний и пересчета MSE на клиентской стороне
callback = CustomJS(args=dict(source_data=source_data, source_pred=source_pred, source_hist=source_hist, 
                              degree_slider=degree_slider, tau_slider=tau_slider, mse_text=mse_text), code="""
    const degree = degree_slider.value;  // Степень полинома
    const tau = tau_slider.value;  // Параметр регуляризации
    let x = source_data.data['x'].slice();  // Копируем массив x
    let y_true = source_data.data['y'].slice();  // Копируем массив y
    const n = x.length;
    
    // Сортируем x и соответствующие y_true по возрастанию x
    const indices = Array.from(x.keys()).sort((a, b) => x[a] - x[b]);
    x = indices.map(i => x[i]);
    y_true = indices.map(i => y_true[i]);

    // Генерация более плотной сетки для предсказаний
    const grid_points = 200;  // Количество точек для более плотной сетки
    const x_min = Math.min(...x);
    const x_max = Math.max(...x);
    const x_dense = [];
    const step = (x_max - x_min) / (grid_points - 1);
    
    for (let i = 0; i < grid_points; i++) {
        x_dense.push(x_min + i * step);  // Заполняем массив точками на сетке
    }

    // Матрица X для полиномиальной регрессии на исходных данных
    const X_poly = [];
    for (let i = 0; i < n; i++) {
        const row = [];
        for (let d = 0; d <= degree; d++) {
            row.push(Math.pow(x[i], d));
        }
        X_poly.push(row);
    }
    
    // Решение линейной регрессии с Ridge регуляризацией
    const XT = X_poly[0].map((_, colIndex) => X_poly.map(row => row[colIndex]));
    const XTX = XT.map(row => row.map((_, colIndex) => row.reduce((acc, curr, i) => acc + curr * X_poly[i][colIndex], 0)));
    const XTy = XT.map(row => row.reduce((acc, curr, i) => acc + curr * y_true[i], 0));
    
    // Добавление регуляризации к XTX
    for (let i = 0; i < XTX.length; i++) {
        XTX[i][i] += tau;
    }
    
    function solve(A, b) {
        const n = A.length;
        const x = Array(n).fill(0);
        for (let i = 0; i < n; i++) {
            let maxEl = Math.abs(A[i][i]), maxRow = i;
            for (let k = i + 1; k < n; k++) {
                if (Math.abs(A[k][i]) > maxEl) {
                    maxEl = Math.abs(A[k][i]);
                    maxRow = k;
                }
            }
            [A[maxRow], A[i]] = [A[i], A[maxRow]];
            [b[maxRow], b[i]] = [b[i], b[maxRow]];
            for (let k = i + 1; k < n; k++) {
                const c = -A[k][i] / A[i][i];
                for (let j = i; j < n; j++) {
                    A[k][j] += c * A[i][j];
                }
                b[k] += c * b[i];
            }
        }
        for (let i = n - 1; i >= 0; i--) {
            x[i] = b[i] / A[i][i];
            for (let k = i - 1; k >= 0; k--) {
                b[k] -= A[k][i] * x[i];
            }
        }
        return x;
    }
    
    const coefs = solve(XTX, XTy);
    
    // Предсказания на плотной сетке
    const y_dense_pred = [];
    for (let i = 0; i < grid_points; i++) {
        let pred = 0;
        for (let d = 0; d <= degree; d++) {
            pred += coefs[d] * Math.pow(x_dense[i], d);
        }
        y_dense_pred.push(pred);
    }

    // Обновление предсказаний на графике (плотная сетка)
    source_pred.data['x'] = x_dense;  // Заменяем на плотную сетку
    source_pred.data['y'] = y_dense_pred;  // Обновляем предсказанные y
    source_pred.change.emit();
    
    // Вычисление MSE на исходных данных
    let mse = 0;
    for (let i = 0; i < n; i++) {
        let pred = 0;
        for (let d = 0; d <= degree; d++) {
            pred += coefs[d] * Math.pow(x[i], d);
        }
        mse += Math.pow(y_true[i] - pred, 2);
    }
    mse /= n;
    
    // Обновление текста с MSE
    mse_text.text = `MSE: ${mse.toFixed(5)}`;
    
    // Обновление гистограммы весов
    const top = coefs;  // Значения коэффициентов
    const left = Array.from({length: coefs.length}, (_, i) => i - 0.4);  // Индексы коэффициентов для left
    const right = Array.from({length: coefs.length}, (_, i) => i + 0.4);  // Индексы коэффициентов для right

    source_hist.data['top'] = top;
    source_hist.data['left'] = left;
    source_hist.data['right'] = right;
    source_hist.change.emit();
""")

# Привязка ползунков к callback
degree_slider.js_on_change('value', callback)
tau_slider.js_on_change('value', callback)

# Компоновка
layout = column(row(p, hist_fig), degree_slider, tau_slider, mse_text)

# Визуализация
show(layout)

Видим, что теперь веса модели не взрываются + смогли получить хорошую аппроксимацию даже при большом \(m\)! Главное - не переборщить с величиной параметра \(\tau\)

Визуализация Lasso регуляризации#

Hide code cell source 
import numpy as np
from sklearn.preprocessing import PolynomialFeatures
from bokeh.layouts import column, row
from bokeh.models import Slider, ColumnDataSource, Div, CustomJS
from bokeh.plotting import figure, show

# Генерация данных
np.random.seed(42)
X = np.random.rand(50, 1) * 2 - 1  # Данные в диапазоне [-1, 1]
y = -4 + 3 * X - 2 * X**2 + 5 * X**3 + np.random.randn(50, 1) * 0.5  # Кубическая зависимость с шумом

# Подготовка данных для Bokeh
X_flat = X.flatten()
y_flat = y.flatten()

# Источник данных
source_data = ColumnDataSource(data=dict(x=X_flat, y=y_flat))
source_pred = ColumnDataSource(data=dict(x=X_flat, y=np.zeros_like(y_flat)))
source_hist = ColumnDataSource(data=dict(top=[], left=[], right=[]))  # Данные для гистограммы

# График данных и предсказаний
p = figure(y_range=(-14, 2), x_range=(-1, 1), width=400, height=400, title="Polynomial Fitting with Lasso Regularization")
p.scatter('x', 'y', source=source_data, color="red", legend_label="Data")
pred_line = p.line('x', 'y', source=source_pred, color="blue", legend_label="Prediction")
p.legend.location = "bottom_right"

# График гистограммы весов
hist_fig = figure(width=400, height=400, title="Weights Histogram")
hist_fig.quad(top='top', bottom=0, left='left', right='right', source=source_hist, fill_color="green", line_color="black")

# Ползунки для степени полинома и параметра регуляризации
degree_slider = Slider(start=1, end=49, value=1, step=1, title="Degree of Polynomial")
tau_slider = Slider(start=0, end=0.5, value=0, step=0.001, title="Regularization Parameter (mu)")

# Текст для отображения MSE
mse_text = Div(text=f"MSE: 0")

# CustomJS для обновления предсказаний и пересчета MSE на клиентской стороне
callback = CustomJS(args=dict(source_data=source_data, source_pred=source_pred, source_hist=source_hist, 
                              degree_slider=degree_slider, tau_slider=tau_slider, mse_text=mse_text), code="""
    const degree = degree_slider.value;  // Степень полинома
    const tau = tau_slider.value;  // Параметр регуляризации
    let x = source_data.data['x'].slice();  // Копируем массив x
    let y_true = source_data.data['y'].slice();  // Копируем массив y
    const n = x.length;

    // Сортируем x и соответствующие y_true по возрастанию x
    const indices = Array.from(x.keys()).sort((a, b) => x[a] - x[b]);
    x = indices.map(i => x[i]);
    y_true = indices.map(i => y_true[i]);

    // Генерация более плотной сетки для предсказаний
    const grid_points = 200;  // Количество точек для более плотной сетки
    const x_min = Math.min(...x);
    const x_max = Math.max(...x);
    const x_dense = [];
    const step = (x_max - x_min) / (grid_points - 1);

    for (let i = 0; i < grid_points; i++) {
        x_dense.push(x_min + i * step);  // Заполняем массив точками на сетке
    }

    // Матрица X для полиномиальной регрессии на исходных данных
    const X_poly = [];
    for (let i = 0; i < n; i++) {
        const row = [];
        for (let d = 0; d <= degree; d++) {
            row.push(Math.pow(x[i], d));
        }
        X_poly.push(row);
    }
    let current_weights = null;  // Глобальная переменная для хранения весов между обновлениями
    
    function lassoCoordinateDescent(X, y, tau, max_iter=500, tol=1e-3) {
        const m = X.length;
        const n = X[0].length;
        
        // Инициализируем веса: если есть текущие веса, используем их (warm start)
        let w = current_weights ? [...current_weights] : new Float64Array(n);
        let w_old = new Float64Array(n);  // Старые веса
    
        // Предварительные вычисления
        const X_T = X[0].map((_, colIndex) => X.map(row => row[colIndex]));  // Транспонируем X
        const X_T_squared = X_T.map(col => col.reduce((acc, x_ij) => acc + x_ij * x_ij, 0));  // X^T * X для каждого столбца
    
        for (let iter = 0; iter < max_iter; iter++) {
            let max_change = 0;
            for (let j = 0; j < n; j++) {
                const residuals = y.map((y_i, i) => y_i - X[i].reduce((acc, x_ij, k) => acc + (k !== j ? x_ij * w[k] : 0), 0));
                const rho_j = X_T[j].reduce((acc, x_ij, i) => acc + x_ij * residuals[i], 0);
    
                const w_j_old = w[j];
                if (rho_j < -tau / 2) {
                    w[j] = (rho_j + tau / 2) / X_T_squared[j];
                } else if (rho_j > tau / 2) {
                    w[j] = (rho_j - tau / 2) / X_T_squared[j];
                } else {
                    w[j] = 0;
                }
    
                max_change = Math.max(max_change, Math.abs(w[j] - w_j_old));
            }
    
            // Условие остановки
            if (max_change < tol) break;
        }
    
        // Сохраняем текущие веса для использования в следующей итерации (warm start)
        current_weights = [...w];
    
        return w;
    }


    // Решение линейной регрессии без регуляризации (метод наименьших квадратов)
    function ordinaryLeastSquares(X, y) {
        const XT = X[0].map((_, colIndex) => X.map(row => row[colIndex]));  // Транспонирование X
        const XTX = XT.map(row => row.map((_, i) => row.reduce((acc, xk, j) => acc + xk * X[j][i], 0)));  // X^T * X
        const XTy = XT.map(row => row.reduce((acc, xk, i) => acc + xk * y[i], 0));  // X^T * y
    
        // Решение системы (X^T * X) w = X^T * y через метод Гаусса или через обратную матрицу
        function gaussJordan(A, b) {
            const n = A.length;
            for (let i = 0; i < n; i++) {
                let maxRow = i;
                for (let k = i + 1; k < n; k++) {
                    if (Math.abs(A[k][i]) > Math.abs(A[maxRow][i])) {
                        maxRow = k;
                    }
                }
    
                // Перестановка строк
                const temp = A[i];
                A[i] = A[maxRow];
                A[maxRow] = temp;
                const t = b[i];
                b[i] = b[maxRow];
                b[maxRow] = t;
    
                // Приводим ведущий элемент к 1
                const div = A[i][i];
                for (let j = i; j < n; j++) A[i][j] /= div;
                b[i] /= div;
    
                // Убираем остальные элементы в столбце
                for (let k = 0; k < n; k++) {
                    if (k !== i) {
                        const factor = A[k][i];
                        for (let j = i; j < n; j++) A[k][j] -= A[i][j] * factor;
                        b[k] -= b[i] * factor;
                    }
                }
            }
            return b;
        }
    
        return gaussJordan(XTX, XTy);  // Возвращаем веса
    }
    
    // Проверяем, какой метод использовать: lasso или OLS
    let coefs;
    if (tau === 0) {
        coefs = ordinaryLeastSquares(X_poly, y_true);
    } else {
        coefs = lassoCoordinateDescent(X_poly, y_true, tau);
    }

    // Предсказания на плотной сетке
    const y_dense_pred = [];
    for (let i = 0; i < grid_points; i++) {
        let pred = 0;
        for (let d = 0; d <= degree; d++) {
            pred += coefs[d] * Math.pow(x_dense[i], d);
        }
        y_dense_pred.push(pred);
    }

    // Обновление предсказаний на графике (плотная сетка)
    source_pred.data['x'] = x_dense;  // Заменяем на плотную сетку
    source_pred.data['y'] = y_dense_pred;  // Обновляем предсказанные y
    source_pred.change.emit();
    
    // Вычисление MSE на исходных данных
    let mse = 0;
    for (let i = 0; i < n; i++) {
        let pred = 0;
        for (let d = 0; d <= degree; d++) {
            pred += coefs[d] * Math.pow(x[i], d);
        }
        mse += Math.pow(y_true[i] - pred, 2);
    }
    mse /= n;
    
    // Обновление текста с MSE
    mse_text.text = `MSE: ${mse.toFixed(5)}`;

    // Обновление гистограммы весов
    const top = coefs;  // Значения коэффициентов
    const left = Array.from({length: coefs.length}, (_, i) => i - 0.4);  // Индексы коэффициентов для left
    const right = Array.from({length: coefs.length}, (_, i) => i + 0.4);  // Индексы коэффициентов для right

    source_hist.data['top'] = top;
    source_hist.data['left'] = left;
    source_hist.data['right'] = right;
    source_hist.change.emit();
""")

# Привязка ползунков к callback
degree_slider.js_on_change('value', callback)
tau_slider.js_on_change('value', callback)

# Компоновка
layout = column(row(p, hist_fig), degree_slider, tau_slider, mse_text)

# Визуализация
show(layout)

Видим, что также смогли получить относительно хорошую аппроксимацию при большой сложности модели.

Заметим, что теперь почти все признаки занулились! Более того, чем больше параметр \(\mu\), тем больше параметров зануляется. Так можно выявлять наиболее важные параметры в задаче регрессии.

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

Содержание