Линейная регрессия++

Эта заметка написана несколько в другом стиле, чем многие предыдущие… Поскольку автор постоянно совершенствует курс по машинному обучению, здесь берётся самая простая и популярная тема классических курсов по ML, и показывается, о чём в ней можно / стоит ещё рассказать (хотя об этом часто забывают), какие здесь есть сложные и интересные вопросы (если Вы хотите проверить свои или чужие знания по линейной регрессии).

Line.pngСначала напомним стандартные вещи. Если в задаче регрессии у нас есть гипотеза о линейной зависимости целевой переменной от исходных признаков, то мы ищем решение в виде линейной комбинации признаков:

f1.png

(здесь мы применяем стандартный приём – вводим фиктивный признак X0=1, чтобы записать результат модели в виде скалярного произведения вектора весов на вектор признаков), тогда задача обучения по конечной выборке формализуется в виде системы линейных равенств:

f2.png

которую сразу можно записать в матричном виде Xw=y, вместо точного решения (которое может и не существовать) обычно решается задача

f4.png

Это означает, что мы настраиваемся на функцию ошибки MSE, мы уже писали при каких предположениях это корректно (хотя разбирать это в теме «Линейная регрессия» я не рекомендую, поскольку обоснование MSE и изучение конкретной модели алгоритмов логически не связано). Решение оптимизационной задачи выписывается в явном виде с помощью векторного дифференцирования

f5.png

Матрица, стоящая перед вектором y называется псевдообратной матрицей Мура-Пенроуза, это обобщение понятия «обратной матрицы» для неквадратных матриц (мнемоническое правило для запоминания: исходное уравнение Xw=y умножается слева на XT, после чего появляется квадратная матрица XTX для обращения). Основные проблемы линейной регрессии связаны с возможной вырожденностью (или просто плохой обусловленностью) матрицы XTX. Заметим, что эта матрица вырождена, например, когда число объектов в обучающей выборке меньше числа признаков. Есть следующие способы борьбы с вырожденностью:

Дальше всегда рассказывают про регуляризацию (про уменьшение размерности логично делать отдельную лекцию). Есть её байесовское обоснование, которое явно рановато в условиях изучения линейной регрессии (обычно это одна из первых тем в курсе ML). Есть несколько приёмов, которые позволяют показать, почему плохо, когда коэффициенты в модели большие без углубления в теорию. Например, объяснить, что значение εjwj может быть большим, когда два объекта отличаются лишь в j-м признаке на εj, а коэффициент wj огромный, что вызывает существенно различные ответы регрессии на похожих объектах. Формально есть много способов «заставить» коэффициенты быть небольшими, в каждом случае наша оптимизационная задача меняется другой:

f6.png

Слева здесь записана регуляризация по Иванову, справа — по Тихонову. Две описаные задачи эквивалентны, в том смысле, что решение одной является и решением другой при определённом выборе регуляризационного коэффициента. Но только в таком смысле! Формально их решения (так, как они записаны) не совпадают. А есть ещё регуляризация по Морозову, что это такое — «хороший» вопрос для желающих завалить кого-нибудь по теме «регуляризация»!

Задача оптимизации при регуляризации по Тихонову также имеет аналитическое решение:

f7.png

Здесь полезно вспомнить сингулярное разложение (SVD), применив его для матрицы X=UΛVT, получаем записанное выше (красным выделено то, что отностится к регуляризации). Эта формула позволяет продемонстрировать, что решение записывается в виде линейной комбинации столбцов матрицы V, коэффициенты зависят от скалярного произведения столбцов матрицы U и целевого вектора, кроме того, от сингулярных чисел матрицы X: когда эти числа близки к нулю соответствующие коэффициенты становятся аномально большими. Добавка λ «отодвигает» знаменатель от нуля. Эти выкладки позволяют ответить, например, на такой вопрос: верно ли, что Ridge-регрессии вектор оптимальных параметров всегда является линейной комбинацией строк  матрицы объект-признак, т.е. X (гребневой или Ridge-регрессией как раз называется линейная регрессией с регуляризацией по Тихонову с помощью L2-нормы).

Многие думают, что регуляризация борется со всеми негативными явлениями в машинном обучении, например, c выбросами. Рис. 1 показывает пример того, что регуляризация не борется с выбросами: чем больше коэффициент регуляризации (здесь он назван α — именно так он называется в библиотеке scikit-learn), тем ещё больше решение отклоняется в сторону выбросов.

h1.jpg
Рис. 1.  Линейная регрессия при разных коэфициентах в L2-регуляризации.

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

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

h2.png
Рис. 2. Коэффициенты линейной регрессии при разной силе L2-регуляризации.

Хороший вопрос: чем отличаются (с точки зрения проведения регуляризации) графики на рис 2 (слева и справа)? Это графики изменения коэффициентов в линейной регуляризации в зависимости от значения множителя при квадрате L2-нормы в одной и той же задаче! В большинстве учебников и слайдах курсов подобные рисунки похожи на левый, а должны быть (если Вы умеете правильно применять регуляризацию) похожи на правый (оговоримся, что в рассматриваемой модельной задаче все коэффициенты в «правильном ответе» положительные).

Регуляризацию можно выполнять выбирая в качестве штрафующего слагаемого не L2-норму весов, а L1 и их комбинации, популярные методы формализуются следующими  оптимизациоными задачами:

f8.png

Тут можно сделать замечание (это типичная ошибка), что «Лассо» – это не ковбойский аркан, а аббревиатура. Обычно поясняют, что L1-регуляризация отлично годится для автоматической селекции признаков. Это можно проиллюстрировать графиками коэффициентов:

f9
Рис. 3. Коэффициенты линейной регрессии при разной силе L1-регуляризации.

Здесь видно, что только несколько коэффициентов, в отличие от L2-регуляризации, отличны от нуля и графики стремительно «ныряют вниз». Кстати, на первом рисунке, когда два веса обнуляются – «разнуляются» другие. Это тоже интересный эффект. Хороший вопрос: когда происходит «разнуление весов» при L1-регуляризации и по каким причинам?

Также зануление часто иллюстрируют линиями уровня минимизируемой функции и соотвестствующими ограничениями (в форме регуляризации Иванова), см. рис. 4. При L2-регуляризации ограничения описываются окружностью (на плоскости), при L1 – ромбом, при L1+L2 – «полненький ромбом». У ромба большая вероятность коснуться линии уровня вершиной, чем у круга верхней / нижней / левой или правой точкой, это соответсвует тому, что в оптимальном решении есть нулевой коэффициент. Хороший вопрос: какая вероятность зануления коэффициента в задаче с двумя признаками и L1-регуляризацией?

h2.jpg
Рис. 4. Линии уровня и ограничения на коэффициенты в модельной задаче.

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

h3.jpg
Рис. 5. Применение метода стохастического градиента: изменнеие параметров(верхние рис.), измениние графика регрессии (левый нижний рис.), изменение ошибки (правый нижний).

Интересный факт, который почти никто не знает, в реализации sklearn.linear_model.Ridge есть разные параметры метода, но они зависимые! Например, если в модели «отключён» свободный член (fit_intercept=False), то игнорируется предварительная нормировка данных (параметр normalize), т.е. комбинация fit_intercept=False, normalize=True (пронормировать, но не использовать свободный член) не реализуема! Этому, кстати, есть логическое объяснение, хороший вопрос — какое?😉

Ещё интересное наблюдение и хороший вопрос по реализации методов: почему регуляризация в Ridge-регрессии контролитруется множителем при норме весов, а вот в SVM и логистической регресии — обратной величиной? Но тут, чтобы ответить, надо иметь опыт копания в книгах по ML…

Важно пояснить, что регрессии y(x) и x(y) (т.е. в зависимости от того, что считать целевым признаком в задаче с двумя признаками), естественно, дают разный результат, см. рис 6. Отличный вопрос: а когда их результат совпадает?

pic13.png
Рис. 6. Линейная регрессия y(x) (синяя) и x(y) (розовая).

Обычно, говоря про линейную регрессию, заточенную под функцию ошибки MSE, упоминают неустойчивость к выборосам. Здесь очень много тонкостей: выбросы, вообще говоря, даже можно оценить по степени влияния на регрессию. Пусть природа данных линейна: y=Xw*+ε, т.е. есть линейная зависимость с точностью до шума (даже не важно, какой это шум). Мы решаем задачу ya=Xw с точностью до оценки(!) шума

h.png

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

Во-вторых, мы ввели матрицу H, которая называется проекционной. Легко вывести, что a=Hy и по смыслу ответ нашего регрессора – это проекция на пространство линейных комбинаций L(X) столбцов матрицы X целевого вектора. Хороший вопрос – ортогональная ли эта проекция?

f2.jpg
Рис. 7. Геометрический смысл проекционной матрицы.

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

h2.png

На вопрос «что такое расстояние Кука?» уже мало кто ответит, а вот более интересный вопрос — почему приведённое выражение называется расстоянием? Расстояние Кука записывается через диагональные элементы проекционной матрицы:

h3.png

Кстати, на рис. 8 показано значения выражений, входящих в формулу в случае настройки обычной линейной регрессии на отрезке [0,1] (точки распределены равномерно, шум нормальный), скажем выброс при X=0.0 в четыре раза влиятельнее выброса при X=0.5.

f1.jpg
Рис. 8. Значения диагональных элементов в модельной задаче в зависимости от x.

Если говорить про линейную регрессию, устойчивую к выбросам, то следует упомянуть метод RANdom SAmple Consensus (RANSAC), простой, довольно эффективный и имеющий реализацию в scikit-learn. Метод заключается в выполнении следующих действий:

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

На рис. 8 показана одна итерация метом RANSAC. Строго говоря, описанный метод годится для любого регрессора, но чаще используется с линейным. Здесь много тонких вопросов по каждому из пунктов: насколько большое базовое множество выбирать, как обучать, как оценивать качество и т.п. Например, в scikit-learn под качеством понимается число базовых точек модели (после пополнения), если у двух моделей число базовых точек совпадает, то сравнивается средняя ошибка на всей выборке.

f3.jpg

На все описанные моменты я обращаю внимание на лекциях в рамках курса «Машинное обучение» программы Ozonmasters, поэтому можно оценить, насколько курс похож на другие и чем отличается на примере простейшей темы (ещё у меня эта тема есть в рамках потокового курса на ВМК МГУ, но там я уже не успеваю обращать на всё внимание — курс короткий).

Линейная регрессия++: 13 комментариев

  1. А как все эти знания помогут улучшать качество модели, например для решения Kaggle?

    У той же sklearn.linear_model.ElasticNet мало параметров, можно всё запихнуть в HyperOpt (параметры и признаки) и найти наилучшее решение за ночь.

    • Например, если уметь отвечать на вопрос «чем отличаются (с точки зрения проведения регуляризации) графики на рис 2 (слева и справа)? «, то значит уметь проводить правильно регуляризацию. Если хотите пример из Кэгла, то я со студентами решал задачу на кэгле — был в топ 10, а они во второй половине LB, хотя реализовывали один и тот же метод, но они не делали правильно регуляризацию.

      Ну, вот, а я ElasticNet настрою минут за 10… перебигать гипероптом параметры для ElasticNet это изумительно;)

      • Ну так hyperopt просто успеет перебрать прорву комбинаций столбцов, включенных в датасет + параметров самой ElasticNet с точки зрения качества на отложенной выборке, почему так удивительно?

        Если хочу только коэффициенты настроить — я просто по сетке смогу пробежать и отрисовать их. За 10 минут тоже управлюсь.

      • Есть, кстати, такая рекомендация (она не доказана формально, но раньше гуляла по ML-форумам): никогда не использовать для селекции признаков и для решения задачи один и тот же алгоритм.

  2. Хороший вопрос: чем отличаются (с точки зрения проведения регуляризации) графики на рис 2 (слева и справа)?

    — Слева график со всеми признаками. Справа выкинуты признаки, которые разнуляются при увеличении регуляризации.

    Хороший вопрос: когда происходит «разнуление весов» при L1-регуляризации и по каким причинам?

    — Одна из причин — скореллированные признаки, когда модели всё равно какой из признаков выбрать, она при занулении веса первого признака — разнуляет второй признак.

    Хороший вопрос: какая вероятность зануления коэффициента в задаче с двумя признаками и L1-регуляризацией?

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

    • > Слева график со всеми признаками.
      Нет, совсем не так.

      > Одна из причин — скореллированные признаки
      Да, так и есть. Более тонкий вопрос, какие признаки при этом разнуляются…

      > вероятность равна вероятности данным двум признакам быть скоррелированными.
      Тут надо, конечно, более подробно формализовать задачу. Но, как видно на рис., при возрастаннии коэф. L1-регуляризации признаки начинают зануляться, хотя они и не коррелировали (а при L2 — просто стремятся к нулю).

  3. // а должны быть (если Вы умеете правильно применять регуляризацию) похожи на правый

    А почему регуляризация должна уменьшать все коэффициэнты равномерно?

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

      • Хорошо, выражусь точнее — уменьшение, а не занулени.

        На правом графике коэффициенты уменьшаются одновременно, на правом одни быстрее/другие медленнее/третьи немного увеличиваются перед падением.

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

      • Давайте так, чтобы объяснить, что я имел в виду надо

        1) ответить на вопрос «Чем отличаются эти графики с точки зрения проведения регуляризации?»

        2) попробуйте рассмотреть типичную задачу с двумя признаками (см. рис 4) и ответить на вопрос, когда веса ведут себя как на правом рис, а когда как на левом (при увеличении коэф. регуляризации)

        3) Давайте возьмём самую примитивную модельную задачу
        Скажем Y = 1*X_1 + 2*X_2 + 3*X_3 (можно + слабый шум)
        Все X_i независимые, но одинаково распределённые с.в.
        Делаем линейную регрессию с L2-регуляризацией, повышаем коэффициент… по Вашей логике, что должно происходить? По моей — как на левом рисунке (можно провести эксперимент и убедиться)

  4. #### Давайте возьмём самую примитивную модельную задачу ####
    #### Скажем Y = 1*X_1 + 2*X_2 + 3*X_3 (можно + слабый шум) ####

    import pandas as pd
    import numpy as np

    import matplotlib.pyplot as plt
    import seaborn as sns
    import warnings; warnings.simplefilter(«ignore»)

    from sklearn.linear_model import Ridge#Lasso
    from sklearn.preprocessing import StandardScaler

    plt.style.use(‘seaborn-whitegrid’)
    sns.set_style(«white»)
    %matplotlib inline

    def make_s(n_rows):
    tmp = pd.DataFrame({‘x1’: (100*np.random.rand(n_rows)).astype(int),
    ‘x2’: (100*np.random.rand(n_rows)).astype(int),
    ‘x3′: (100*np.random.rand(n_rows)).astype(int),
    #’x4′: (100*np.random.rand(n_rows)).astype(int),
    #’x5’: (100*np.random.rand(n_rows)).astype(int),
    ‘noise’: (100*np.random.rand(n_rows)).astype(int)
    })
    #tmp[‘x5’] = tmp[‘x4’] * 2 + tmp[‘x3’] * 3 — 1* tmp[‘noise’]
    tmp[‘y’] = 1*tmp[‘x1’] + 2*tmp[‘x2’] + 3*tmp[‘x3’]
    #+8*tmp[‘x4’] + 9*tmp[‘x5’]
    + tmp[‘noise’]
    tmp.drop([‘noise’], axis = 1, inplace = True)

    return tmp

    data = make_s(500000)

    scaler = StandardScaler()
    X_train = scaler.fit_transform(data.drop([‘y’], axis = 1))
    y_train = data[‘y’]

    def plot_weights(alpha = 1):
    model = Ridge(alpha=alpha)
    model.fit(X_train, y_train)
    coef = model.coef_
    return list(np.append(np.abs(coef), np.round(alpha,2)))

    res = []

    N_ROWS_NEGATIVE = 100
    N_ROWS_POSITIVE = 100

    for i in reversed(range(N_ROWS_NEGATIVE)):
    res.append(plot_weights(alpha = 10 ** ((-i)/10) ))
    for i in range(N_ROWS_POSITIVE):
    res.append(plot_weights(alpha = 10 ** ((i+1)/10) ))

    res_df = pd.DataFrame(res)
    res_df.columns = data.columns
    res_df.rename(index=str, columns={‘y’: ‘alpha’}, inplace=True)

    plt.figure(figsize=(12,10), dpi= 80)
    plt.plot(res_df.drop([‘alpha’], axis = 1))
    plt.xticks(np.arange(res_df[‘alpha’].shape[0]), res_df[‘alpha’], rotation=90)
    plt.legend(res_df.columns);

    Не могу к сожалению вставить картинку 😦
    Получилось как на правом рисунке L1-регуляризаци

    • Вот картинка с Вашего кода: http://alexanderdyakonov.narod.ru/tmp_regul.png

      Я, кстати, описался… «По моей — как на левом рисунке» -> «По моей — как на правом рисунке»

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

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

      А я говорил, что если всё правильно делать — должно быть как на правом графике… и Вы всё правильно сделали!

      Естественно, что при решении более сложных задач должно быть точно также, как в простом идеальном случае;)

      • #### Тут получается, что надо просто после построения графика — выкидывать колонки с самыми большими коэффициентам, для оптимизации? Просто вообще не понятно что с этим дальше делать. Сделал функции, которые строят графики RMSE (MAE) в зависимости от рандомно задаваемых в различных распределениях данных. Не понимаю что с этим теперь делать дальше, потому что мой вариант в Hyperopt находит глобальный минимум, а этот с графиками — нет 😦
        ####

        import pandas as pd
        import numpy as np
        import matplotlib.pyplot as plt
        import warnings; warnings.simplefilter(«ignore»)

        from sklearn.linear_model import Lasso, Ridge
        from sklearn.preprocessing import StandardScaler
        from sklearn.metrics import mean_squared_error, mean_absolute_error
        from sklearn.model_selection import train_test_split

        %matplotlib inline

        def make_s(n_rows=5000, func=’rand’, n_rows_negative = 3, n_rows_positive = 40):
        if func == ‘rand’:
        tmp = pd.DataFrame({‘x1’: (100*np.random.rand(n_rows)).astype(int),
        ‘x2’: (100*np.random.rand(n_rows)).astype(int),
        ‘x3’: (100*np.random.rand(n_rows)).astype(int),
        ‘x4’: (100*np.random.rand(n_rows)).astype(int),
        ‘x5’: (100*np.random.rand(n_rows)).astype(int),
        ‘noise’: 2*(np.random.rand(n_rows)).astype(int)
        })
        if func == ‘norm’:
        tmp = pd.DataFrame({‘x1’: np.random.normal(0, 100, n_rows),
        ‘x2’: np.random.normal(0, 100, n_rows),
        ‘x3’: np.random.normal(0, 100, n_rows),
        ‘x4’: np.random.normal(0, 100, n_rows),
        ‘x5’: np.random.normal(0, 100, n_rows),
        ‘noise’: np.random.normal(0, 100, n_rows)
        })
        #tmp[‘x4’] = tmp[‘x3’] * 20 — 10* tmp[‘noise’]
        tmp[‘x5’] = 5*tmp[‘x1’] — 2*tmp[‘x2’] + 3*tmp[‘x3’] — 4*tmp[‘x4’] + 5*tmp[‘x5’] — 10* tmp[‘noise’]

        tmp[‘y’] = 1*tmp[‘x1’] + 2*tmp[‘x2’] + 3*tmp[‘x3’] +4*tmp[‘x4’] + 5*tmp[‘x5’] + tmp[‘noise’]

        tmp.drop([‘noise’], axis = 1, inplace = True)

        X_train, X_test, y_train, y_test = train_test_split(tmp.drop([‘y’], axis = 1), tmp[‘y’])

        scaler = StandardScaler()
        X_train = scaler.fit_transform(X_train)
        X_test = scaler.fit_transform(X_test)

        res = []
        def create_dataset(alpha=1):
        model = Ridge(alpha=alpha)#Lasso
        model.fit(X_train, y_train)
        coef = model.coef_
        y_pred = model.predict(X_test)
        rmse = mean_squared_error(y_test, y_pred)#mean_absolute_error
        res = list(np.append(np.abs(coef), [np.round(alpha,2), np.round(rmse)]))
        return res

        for i in reversed(range(n_rows_negative)):
        res.append(create_dataset(alpha = 10 ** ((-i)/10) ))

        for i in range(n_rows_positive):
        res.append(create_dataset(alpha = 10 ** ((i+1)/10) ))

        df = pd.DataFrame(res)
        df.columns = [‘x1’, ‘x2’, ‘x3’, ‘x4’, ‘x5’, ‘alpha’, ‘rmse’]

        return df

        def plot_df(res_df):
        fig, ax1 = plt.subplots(figsize=(12,10), dpi= 80)
        plt.xticks(np.arange(res_df[‘alpha’].shape[0]), res_df[‘alpha’], rotation=90)
        ax1.plot(res_df.drop([‘alpha’, ‘rmse’], axis = 1))
        ax1.legend(res_df.columns)

        ax2 = ax1.twinx()

        color = ‘tab:red’
        ax2.set_ylabel(‘rmse’, color=color)
        ax2.plot(res_df[‘rmse’], color=color, marker=’o’, linestyle=’dashed’,
        linewidth=3)
        ax2.tick_params(axis=’y’)
        plt.yscale(‘linear’);

        df = make_s(n_rows=5000, func=’rand’, n_rows_positive = 50)
        plot_df(df)
        df_norm = make_s(n_rows=5000, func=’norm’, n_rows_positive = 50)
        plot_df(df_norm)

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

Заполните поля или щелкните по значку, чтобы оставить свой комментарий:

Логотип WordPress.com

Для комментария используется ваша учётная запись WordPress.com. Выход /  Изменить )

Google photo

Для комментария используется ваша учётная запись Google. Выход /  Изменить )

Фотография Twitter

Для комментария используется ваша учётная запись Twitter. Выход /  Изменить )

Фотография Facebook

Для комментария используется ваша учётная запись Facebook. Выход /  Изменить )

Connecting to %s