Как найти функцию регрессии

Основы линейной регрессии

Здравствуй, Хабр!

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

! Осторожно, трафик! В статье присутствует заметное число изображений для иллюстраций, часть в формате gif.

Содержание

• Введение
• Метод наименьших квадратов
  • Математический анализ
  • Статистика
  • Теория вероятностей
• Мультилинейная регрессия
  • Линейная алгебра
• Произвольный базис
• Заключительные замечания
  • Проблема выбора размерности
  • Численные методы
• Реклама и заключение

Введение

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

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

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

Регрессия — способ выбрать из семейства функций ту, которая минимизирует функцию потерь. Последняя характеризует насколько сильно пробная функция отклоняется от значений в заданных точках. Если точки получены в эксперименте, они неизбежно содержат ошибку измерений, шум, поэтому разумнее требовать, чтобы функция передавала общую тенденцию, а не точно проходила через все точки. В каком-то смысле регрессия — это «интерполирующая аппроксимация»: мы хотим провести кривую как можно ближе к точкам и при этом сохранить ее максимально простой чтобы уловить общую тенденцию. За баланс между этими противоречивыми желаниями как-раз отвечает функция потерь (в английской литературе «loss function» или «cost function»).

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

Цель регрессии — найти коэффициенты этой линейной комбинации, и тем самым определить регрессионную функцию

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

Регрессия с нами уже давно: впервые метод опубликовал Лежандр в 1805 году, хотя Гаусс пришел к нему раньше и успешно использовал для предсказания орбиты «кометы» (на самом деле карликовой планеты) Цереры. Существует множество вариантов и обобщений линейной регрессии: LAD, метод наименьших квадратов, Ridge регрессия, Lasso регрессия, ElasticNet и многие другие.

Метод наименьших квадратов

Начнём с простейшего двумерного случая. Пусть нам даны точки на плоскости

и мы ищем такую аффинную функцию

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

.

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

поэтому измерять нужно вдоль оси

.

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

. Простейший вариант — сумма модулей отклонений

приводит к Least Absolute Distance (LAD) регрессии.

Впрочем, более популярная функция потерь — сумма квадратов отклонений регрессанта от модели. В англоязычной литературе она носит название Sum of Squared Errors (SSE)

Метод наименьших квадратов (по англ. OLS) — линейная регрессия c

в качестве функции потерь.

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

.

Математический анализ

Простейший способ найти

— вычислить частные производные по

и

, приравнять их нулю и решить систему линейных уравнений

Значения параметров, минимизирующие функцию потерь, удовлетворяют уравнениям

которые легко решить

Мы получили громоздкие и неструктурированные выражения. Сейчас мы их облагородим и вдохнем в них смысл.

Статистика

Полученные формулы можно компактно записать с помощью статистических эстиматоров: среднего

, вариации

(стандартного отклонения), ковариации

и корреляции

Перепишем

как

где

это нескорректированное (смещенное) стандартное выборочное отклонение, а

— ковариация. Теперь вспомним, что коэффициент корреляции (коэффициент корреляции Пирсона)

и запишем

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

Во-первых, это уравнение сразу указывает на два свойства регрессионной прямой:

Во-вторых, теперь становится понятно, почему метод регрессии называется именно так. В единицах стандартного отклонения

отклоняется от своего среднего значения меньше чем

, потому что

. Это называется регрессией(от лат. regressus — «возвращение») по отношению к среднему. Это явление было описано сэром Фрэнсисом Гальтоном в конце XIX века в его статье «Регрессия к посредственности при наследовании роста». В статье показано, что черты (такие как рост), сильно отклоняющиеся от средних, редко передаются по наследству. Характеристики потомства как бы стремятся к среднему — на детях гениев природа отдыхает.

Возведя коэффициент корреляции в квадрат, получим коэффициент детерминации

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

, равный

, означает что функция идеально ложится на все точки — данные идеально скоррелированны. Можно доказать, что

показывает какая доля вариативности в данных объясняется лучшей из линейных моделей. Чтобы понять, что это значит, введем определения

— вариация исходных данных (вариация точек

).

— вариация остатков, то есть вариация отклонений от регрессионной модели — от

нужно отнять предсказание модели и найти вариацию.

— вариация регрессии, то есть вариация предсказаний регрессионной модели в точках

(обратите внимание, что среднее предсказаний модели совпадает с

).

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

или

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

Мы стремимся избавиться от вариативности, связанной с шумом и оставить лишь вариативность, которая объясняется моделью, — хотим отделить зерна от плевел. О том, насколько это удалось лучшей из линейных моделей, свидетельствует

, равный единице минус доля вариации ошибок в суммарной вариации

или доле объясненной вариации (доля вариации регрессии в полной вариации)

равен косинусу угла в прямоугольном треугольнике

. Кстати, иногда вводят долю необъясненной вариации

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

Теория вероятностей

Ранее мы пришли к функции потерь

из соображений удобства, но к ней же можно прийти с помощью теории вероятностей и метода максимального правдоподобия (ММП). Напомню вкратце его суть. Предположим, у нас есть

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

и

). Для этого нужно вычислить вероятность получить

датапоинтов в предположении постоянных, но пока неизвестных параметров. Благодаря независимости измерений, мы получим произведение вероятностей реализации каждого измерения. Если мыслить полученную величину как функцию параметров (функция правдоподобия) и найти её максимум, мы получим оценку параметров. Зачастую вместо функции правдоподобия используют ее логарифм — дифференцировать его проще, а результат — тот же.

Вернемся к задаче простой регрессии. Допустим, что значения

нам известны точно, а в измерении

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

где

— нормально распределенная случайная величина

Исходя из предположений выше, запишем функцию правдоподобия

и ее логарифм

Таким образом, максимум правдоподобия достигается при минимуме

что дает основание принять ее в качестве функции потерь. Кстати, если

мы получим функцию потерь LAD регрессии

которую мы упоминали ранее.

Подход, который мы использовали в этом разделе — один из возможных. Можно прийти к такому же результату, используя более общие свойства. В частности, свойство постоянства распределения можно ослабить, заменив на свойства независимости, постоянства вариации (гомоскедастичность) и отсутствия мультиколлинеарности. Также вместо ММП эстимации можно воспользоваться другими методами, например линейной MMSE эстимацией.

Мультилинейная регрессия

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

, однако обычно регрессор — это

-мерный вектор

. Другими словами, для каждого измерения мы регистрируем

фич, объединяя их в вектор. В этом случае логично принять модель с

независимыми базисными функциями векторного аргумента —

степеней свободы соответствуют

фичам и еще одна — регрессанту

. Простейший выбор — линейные базисные функции

. При

получим уже знакомый нам базис

.

Итак, мы хотим найти такой вектор (набор коэффициентов)

, что

Знак «

» означает, что мы ищем решение, которое минимизирует сумму квадратов ошибок

Последнее уравнение можно переписать более удобным образом. Для этого расположим

в строках матрицы (матрицы информации)

Тогда столбцы матрицы

отвечают измерениям

-ой фичи. Здесь важно не запутаться:

— количество измерений,

— количество признаков (фич), которые мы регистрируем. Систему можно записать как

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

которую мы намерены минимизировать

Продифференцируем финальное выражение по

(если забыли как это делается — загляните в Matrix cookbook)

приравняем производную к

и получим т.н. нормальные уравнения

Если столбцы матрицы информации

линейно независимы (нет идеально скоррелированных фич), то матрица

имеет обратную (доказательство можно посмотреть, например, в видео академии Хана). Тогда можно записать

где

псевдообратная к

. Понятие псевдообратной матрицы введено в 1903 году Фредгольмом, она сыграла важную роль в работах Мура и Пенроуза.

Напомню, что обратить

и найти

можно только если столбцы

линейно независимы. Впрочем, если столбцы

близки к линейной зависимости, вычисление

уже становится численно нестабильным. Степень линейной зависимости признаков в

или, как говорят, мультиколлинеарности матрицы

, можно измерить числом обусловленности — отношением максимального собственного значения к минимальному. Чем оно больше, тем ближе

к вырожденной и неустойчивее вычисление псевдообратной.

Линейная алгебра

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

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

, образ которого

ближе остальных к

. Напомню, что множество образов или колоночное пространство

— это линейная комбинация вектор-столбцов матрицы

—

-мерное линейное подпространство (мы считаем фичи линейно независимыми), линейная оболочка вектор-столбцов

. Итак, если

принадлежит

, то мы можем найти решение, если нет — будем искать, так сказать, лучшее из нерешений.

Если в дополнение к векторам

мы рассмотрим все вектора им перпендикулярные, то получим еще одно подпространство и сможем любой вектор из

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

, тогда

равен нулю в том и только в том случае, если

перпендикулярен всем

, а значит и целому

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

. Иногда это обозначают c помощью символа ортогональной прямой суммы

где

. В каждое из подпространств можно попасть с помощью соответствующего оператора проекции, но об этом ниже.

Теперь представим

в виде разложения

Если мы ищем решение

, то естественно потребовать, чтобы

была минимальна, ведь это длина вектора-остатка. Учитывая перпендикулярность подпространств и теорему Пифагора

но поскольку, выбрав подходящий

, я могу получить любой вектор колоночного пространства, то задача сводится к

а

останется в качестве неустранимой ошибки. Любой другой выбор

сделает ошибку только больше.

Если теперь вспомнить, что

, то легко видеть

что очень удобно, так как

у нас нет, а вот

— есть. Вспомним из предыдущего параграфа, что

имеет обратную при условии линейной независимости признаков и запишем решение

где

уже знакомая нам псевдообратная матрица. Если нам интересна проекция

, то можно записать

где

— оператор проекции на колоночное пространство.

Выясним геометрический смысл коэффициента детерминации.

Заметьте, что фиолетовый вектор

пропорционален первому столбцу матрицы информации

, который состоит из одних единиц согласно нашему выбору базисных функций. В RGB треугольнике

Так как этот треугольник прямоугольный, то по теореме Пифагора

Это геометрическая интерпретация уже известного нам факта, что

Мы знаем, что

а значит

Красиво, не правда ли?

Произвольный базис

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

и её результатом есть модель

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

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

. Как можно было заметить, на самом деле ни вид

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

ложатся на параболу, а не на прямую, то стоит выбрать базис

. Количество базисных функций может быть как меньшим, так и большим, чем количество изначальных фич.

Если мы определились с базисом, то дальше действуем следующим образом. Мы формируем матрицу информации

записываем функцию потерь

и находим её минимум, например с помощью псевдообратной матрицы

или другим методом.

Заключительные замечания

Проблема выбора размерности

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

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

монотонно растет с ростом размерности базиса. Поэтому вводят скорректированный коэффициент

где

— размер выборки,

— количество независимых переменных. Следя за

, мы можем вовремя остановиться и перестать добавлять дополнительные степени свободы.

Вторая группа подходов — регуляризации, самые известные из которых Ridge(

/гребневая/Тихоновская регуляризация), Lasso(

регуляризация) и Elastic Net(Ridge+Lasso). Главная идея этих методов: модифицировать функцию потерь дополнительными слагаемыми, которые не позволят вектору коэффициентов

неограниченно расти и тем самым воспрепятствуют переобучению

где

и

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

Численные методы

Скажу пару слов, как минимизировать функцию потерь на практике. SSE — это обычная квадратичная функция, которая параметризируется входными данными, так что принципиально ее можно минимизировать методом скорейшего спуска или другими методами оптимизации. Разумеется, лучшие результаты показывают алгоритмы, которые учитывают вид функции SSE, например метод стохастического градиентного спуска. Реализация Lasso регрессии в scikit-learn использует метод координатного спуска.

Также можно решить нормальные уравнения с помощью численных методов линейной алгебры. Эффективный метод, который используется в scikit-learn для МНК — нахождение псевдообратной матрицы с помощью сингулярного разложения. Поля этой статьи слишком узки, чтобы касаться этой темы, за подробностями советую обратиться к курсу лекций К.В.Воронцова.

Реклама и заключение

Эта статья — сокращенный пересказ одной из глав курса по классическому машинному обучению в Киевском академическом университете (преемник Киевского отделения Московского физико-технического института, КО МФТИ). Автор статьи помогал в создании этого курса. Технически курс выполнен на платформе Google Colab, что позволяет совмещать формулы, форматированные LaTeX, исполняемый код Python и интерактивные демонстрации на Python+JavaScript, так что студенты могут работать с материалами курса и запускать код с любого компьютера, на котором есть браузер. На главной странице собраны ссылки на конспекты, «рабочие тетради» для практик и дополнительные ресурсы. В основу курса положены следующие принципы:

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

Если хотите посмотреть на результат — загляните на страничку курса на GitHub.

Надеюсь вам было интересно, спасибо за внимание.

Содержание:

Регрессионный анализ:

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

Пусть между случайными величинами X и Y существует линейная корреляционная зависимость. Это означает, что математическое ожидание Y линейно зависит от значений случайной величины X. График этой зависимости (линия регрессии Y на X) имеет уравнение

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

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

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

Такие значения , согласно необходимым условиям экстремума, находятся из системы уравнений:

Решения этой системы уравнений дают оценки называемые оценками по методу наименьших квадратов.

и

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

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

Пример:

По данным наблюдений двух случайных величин найти коэффициент корреляции и уравнение линии регрессии Y на X 

Решение. Вычислим величины, необходимые для использования формул (3.7.1)–(3.7.3):

 

По формулам (3.7.1) и (3.7.2) получим

Итак, оценка линии регрессии имеет вид Так как то по формуле (3.1.3)

Аналогично, Поэтому в качестве оценки коэффициента корреляции имеем по формуле (3.7.3) величину 

Ответ.  

Пример:

Получена выборка значений величин X и Y

Для представления зависимости между величинами предполагается использовать модель Найти оценки параметров

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

Необходимые условия экстремума приводят к системе из двух уравнений:

Откуда

Решения системы уравнений (3.7.4) и (3.7.5) и будут оценками по методу наименьших квадратов для параметров

На основе опытных данных вычисляем:

В итоге получаем систему уравнений (?????) и (?????) в виде

Эта система имеет решения

Ответ. 

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

В этой таблице равно числу наблюдений, для которых X находится в интервале а Y – в интервале Через обозначено число наблюдений, при которых а Y произвольно. Число наблюдений, при которых а X произвольно, обозначено через

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

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

Пример:

Проделано 80 наблюдений случайных величин X и Y. Результаты наблюдений представлены в виде таблицы. Найти линию регрессии Y на X. Оценить коэффициент корреляции.

Решение. Представителем каждого интервала будем считать его середину. Перенесем начало координат и изменим масштаб по каждой оси так, чтобы значения X и Y были удобны для вычислений. Для этого перейдем к новым переменным Значения этих новых переменных указаны соответственно в самой верхней строке и самом левом столбце таблицы.

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

Нанесем эти значения на координатную плоскость, соединив для наглядности их отрезками прямой (рис. 3.7.1).

По виду полученной ломанной линии можно предположить, что линия регрессии Y на X является прямой. Оценим ее параметры. Для этого сначала вычислим с учетом группировки данных в таблице все величины, необходимые для использования формул (3.31–3.33): 

Тогда

В новом масштабе оценка линии регрессии имеет вид График этой прямой линии изображен на рис. 3.7.1.

Для оценки по корреляционной таблице можно воспользоваться формулой (3.1.3):

Подобным же образом можно оценить величиной Тогда оценкой коэффициента корреляции может служить величина

Вернемся к старому масштабу:

 

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

Ответ. 

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

Регрессионный анализ

Основные положения регрессионного анализа:

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

Предпосылки регрессионного анализа:

1. Y — независимые случайные величины, имеющие постоянную дисперсию;
2. X— величины наблюдаемого признака (величины не случайные);
3. условное математическое ожидание можно представить в виде 

Выражение (2.1), как уже упоминалось в п. 1.2, называется функцией регрессии (или модельным уравнением регрессии) Y на X. Оценке в этом выражении подлежат параметры называемые коэффициентами регрессии, а также — остаточная дисперсия.

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

В данном вопросе рассмотрен линейный регрессионный анализ. Линейным он называется потому, что изучаем лишь те виды зависимостей которые линейны по оцениваемым параметрам, хотя могут быть нелинейны по переменным X. Например, зависимости

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

Важное место в линейном регрессионном анализе занимает так называемая «нормальная регрессия». Она имеет место, если сделать предположения относительно закона распределения случайной величины Y. Предпосылки «нормальной регрессии»:

1. Y — независимые случайные величины, имеющие постоянную дисперсию и распределенные по нормальному закону;
2. X— величины наблюдаемого признака (величины не случайные);
3. условное математическое ожидание можно представить в виде (2.1).

В этом случае оценки коэффициентов регрессии — несмещённые с минимальной дисперсией и нормальным законом распределения. Из этого положения следует что при «нормальной регрессии» имеется возможность оценить значимость оценок коэффициентов регрессии, а также построить доверительный интервал для коэффициентов регрессии и условного математического ожидания M(YX=x).

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

Рассмотрим простейший случай регрессионного анализа — модель вида (2.1), когда зависимость линейна и по оцениваемым параметрам, и

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

Составим систему нормальных уравнений: первое уравнение

откуда  

второе уравнение

откуда

Итак,

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

Остаётся получить оценку параметра . Имеем

где т — количество наблюдений.

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

где — частоты повторений соответствующих значений переменных. В дальнейшем часто используется этот наглядный приём вычислений.
 

Нелинейная регрессия

Рассмотрим случай, когда зависимость нелинейна по переменным х, например модель вида
   

На рис. 2.1 изображено поле корреляции. Очевидно, что зависимость между Y и X нелинейная и её графическим изображением является не прямая, а кривая. Оценкой выражения (2.6) является уравнение регрессии

где —оценки коэффициентов регрессии

Принцип нахождения коэффициентов тот же — метод наименьших квадратов, т.е.

или

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

В общем случае нелинейной зависимости между переменными Y и X связь может выражаться многочленом k-й степени от x:

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

Вычислив коэффициенты системы, её можно решить любым известным способом.
 

Оценка значимости коэффициентов регрессии. Интервальная оценка коэффициентов регрессии

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

имеет распределение Стьюдента с к= n-2 степенями свободы (b — оценка коэффициента регрессии, — оценка среднеквадратического отклонения

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

Оценки среднеквадратического отклонения коэффициентов регрессии вычисляют по следующим формулам:

где   — оценка остаточной дисперсии, вычисляемая по
формуле (2.5).

Доверительный интервал для значимых параметров строят по обычной схеме. Из условия

где а — уровень значимости, находим

 

Интервальная оценка для условного математического ожидания

Линия регрессии характеризует изменение условного математического ожидания результативного признака от вариации остальных признаков.

Точечной оценкой условного математического ожидания является условное среднее    Кроме точечной оценки для можно
построить доверительный интервал в точке

Известно, что имеет распределение
Стьюдента с k=n—2 степенями свободы. Найдя оценку среднеквадратического отклонения для условного среднего, можно построить доверительный интервал для условного математического ожидания

Оценку дисперсии условного среднего вычисляют по формуле

или для интервального ряда

Доверительный интервал находят из условия

где а — уровень значимости. Отсюда

Доверительный интервал для условного математического ожидания можно изобразить графически (рис, 2.2).

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

Проверка значимости уравнения регрессии

Оценить значимость уравнения регрессии — значит установить, соответствует ли математическая, модель, выражающая зависимость между Y и X, экспериментальным данным. Для оценки значимости в предпосылках «нормальной регрессии» проверяют гипотезу Если она отвергается, то считают, что между Y и X нет связи (или связь нелинейная). Для проверки нулевой гипотезы используют основное положение дисперсионного анализа о разбиении суммы квадратов на слагаемые. Воспользуемся разложением — Общая сумма квадратов отклонений результативного признака

разлагается на (сумму, характеризующую влияние признака

X) и (остаточную сумму квадратов, характеризующую влияние неучтённых факторов). Очевидно, чем меньше влияние неучтённых факторов, тем лучше математическая модель соответствует экспериментальным данным, так как вариация У в основном объясняется влиянием признака X.

Для проверки нулевой гипотезы вычисляют статистику  которая имеет распределение Фишера-Снедекора с А степенями свободы (в п — число наблюдений). По уровню значимости а и числу степеней свободы находят по таблицам F-распределение для уровня значимости а=0,05 (см. табл. 3 приложений) критическое значение удовлетворяющее условию . Если нулевую гипотезу отвергают, уравнение считают значимым. Если то нет оснований отвергать нулевую гипотезу.

Многомерный регрессионный анализ

В случае, если изменения результативного признака определяются действием совокупности других признаков, имеет место многомерный регрессионный анализ. Пусть результативный признак У, а независимые признаки Для многомерного случая предпосылки регрессионного анализа можно сформулировать следующим образом: У -независимые случайные величины со средним и постоянной дисперсией — линейно независимые векторы . Все положения, изложенные в п.2.1, справедливы для многомерного случая. Рассмотрим модель вида 

Оценке подлежат параметры и остаточная дисперсия.

Заменив параметры их оценками, запишем уравнение регрессии

Коэффициенты в этом выражении находят методом наименьших квадратов.

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

Как и в двумерном случае, составляют систему нормальных уравнений

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

Подставим в остальные уравнения системы вместо полученное выражение:

Пусть С — матрица коэффициентов при неизвестных параметрах — матрица, обратная матрице С; — элемент, стоящий на пересечении i-Й строки и i-го столбца матрицы    — выражение
. Тогда, используя формулы линейной алгебры,

запишем окончательные выражения для параметров:

Оценкой остаточной дисперсии является

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

Если выборка получена из нормально распределенной генеральной совокупности, то, аналогично изложенному в п. 2.4, можно проверить значимость оценок коэффициентов регрессии, только в данном случае статистику вычисляют для каждого j-го коэффициента регрессии

где —элемент обратной матрицы, стоящий на пересечении i-й строки и j-
го столбца; —диагональный элемент обратной матрицы.

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

Для значимых коэффициентов регрессии целесообразно построить доверительные интервалы по формуле (2.10). Для оценки значимости уравнения регрессии следует проверить нулевую гипотезу о том, что все коэффициенты регрессии (кроме свободного члена) равны нулю:  — вектор коэффициентов регрессии). Нулевую гипотезу проверяют, так же как и в п. 2.6, с помощью статистики , где — сумма квадратов, характеризующая влияние признаков X; — остаточная сумма квадратов, характеризующая влияние неучтённых факторов; Для уровня значимости а и числа степеней свободы по табл. 3 приложений находят критическое значение Если то нулевую гипотезу об одновременном равенстве нулю коэффициентов регрессии отвергают. Уравнение регрессии считают значимым. При нет оснований отвергать нулевую гипотезу, уравнение регрессии считают незначимым.

Факторный анализ

Основные положения. В последнее время всё более широкое распространение находит один из новых разделов многомерного статистического анализа — факторный анализ. Первоначально этот метод

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

Один из основоположников факторного анализа Л. Терстоун приводит такой пример: несколько сотен мальчиков выполняют 20 разнообразных гимнастических упражнений. Каждое упражнение оценивают баллами. Можно рассчитать матрицу корреляций между 20 упражнениями. Это большая матрица размером 20><20. Изучая такую матрицу, трудно уловить закономерность связей между упражнениями. Нельзя ли объяснить скрытую в таблице закономерность действием каких-либо обобщённых факторов, которые в результате эксперимента непосредственно, не оценивались? Оказалось, что обо всех коэффициентах корреляции можно судить по трём обобщённым факторам, которые и определяют успех выполнения всех 20 гимнастических упражнений: чувство равновесия, усилие правого плеча, быстрота движения тела.

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

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

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

Представление, информации в факторном анализе

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

Основная модель факторного анализа. Основная модель факторного анализа имеет вид

где -j-й признак (величина случайная); — общие факторы (величины случайные, имеющие нормальный закон распределения); — характерный фактор; — факторные нагрузки, характеризующие существенность влияния каждого фактора (параметры модели, подлежащие определению); — нагрузка характерного фактора.

Модель предполагает, что каждый из j признаков, входящих в исследуемый набор и заданных в стандартной форме, может быть представлен в виде линейной комбинации небольшого числа общих факторов и характерного фактора

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

Термин «характерный фактор» показывает, что он относится только к данному j-му признаку. Это специфика признака, которая не может быть, выражена через факторы

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

Для j-го признака и i-го объекта модель (2.19) можно записать в. виде

где значение k-го фактора для i-го объекта.

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

Если общие и характерные факторы не коррелируют между собой, то дисперсию j-го признака можно представить в виде

где —доля дисперсии признака приходящаяся на k-й фактор.

Полный вклад k-го фактора в суммарную дисперсию признаков

Вклад общих факторов в суммарную дисперсию
 

Факторное отображение

Используя модель (2.19), запишем выражения для каждого из параметров:

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

Факторный анализ позволяет получить не только матрицу отображений, но и коэффициенты корреляции между параметрами и

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

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

Рассмотрим связь между элементами структуры и коэффициентами отображения. Для этого, учитывая выражение (2.19) и определение выборочного коэффициента корреляции, умножим уравнения системы (2.21) на соответствующие факторы, произведём суммирование по всем n наблюдениям и, разделив на n, получим следующую систему уравнений:

где — выборочный коэффициент корреляции между j-м параметром и к-
м фактором; — коэффициент корреляции между к-м и р-м факторами.

Если предположить, что общие факторы между собой, не коррелированы, то уравнения    (2.22) можно записать в виде

, т.е. коэффициенты отображения равны
элементам структуры.

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

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

Матрица остаточных коэффициентов корреляции называется остаточной матрицей или матрицей остатков

где — матрица остатков; R — матрица выборочных парных коэффициентов корреляции, или полная матрица; R’— матрица вычисленных по отображению коэффициентов корреляции.

Результаты факторного анализа удобно представить в виде табл. 2.10.

Здесь суммы квадратов нагрузок по строкам — общности параметров, а суммы квадратов нагрузок по столбцам — вклады факторов в суммарную дисперсию параметров. Имеет место соотношение

Определение факторных нагрузок

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

Введём понятие редуцированной корреляционной матрицы или просто редуцированной матрицы. Редуцированной называется матрица выборочных коэффициентов корреляции у которой на главной диагонали стоят значения общностей :

Редуцированная и полная матрицы связаны соотношением

где D — матрица характерностей.

Общности, как правило, неизвестны, и нахождение их в факторном анализе представляет серьезную проблему. Вначале определяют (хотя бы приближённо) число общих факторов, совокупность, которых может с достаточной точностью аппроксимировать все взаимосвязи выборочной корреляционной матрицы. Доказано, что число общих факторов (общностей) равно рангу редуцированной матрицы, а при известном ранге можно по выборочной корреляционной матрице найти оценки общностей. Числа общих факторов можно определить априори, исходя из физической природы эксперимента. Затем рассчитывают матрицу факторных нагрузок. Такая матрица, рассчитанная методом главных факторов, обладает одним интересным свойством: сумма произведений каждой пары её столбцов равна нулю, т.е. факторы попарно ортогональны.

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

Максимум должен быть найден при условии

где —общностьпараметра

Затем рассчитывают матрицу коэффициентов корреляции с учётом только первого фактора Имея эту матрицу, получают первую матрицу остатков:

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

Максимум находят из условия

где — коэффициент корреляции из первой матрицы остатков; — факторные нагрузки с учётом второго фактора. Затем рассчитыва коэффициентов корреляций с учётом второго фактора и вычисляют вторую матрицу остатков:

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

Адекватность факторной модели оценивается по матрице остатков (если величины её коэффициентов малы, то модель считают адекватной).

Такова последовательность шагов для нахождения факторных нагрузок. Для нахождения максимума функции (2.24) при условии (2.25) используют метод множителей Лагранжа, который приводит к системе т уравнений относительно m неизвестных

Метод главных компонент

Разновидностью метода главных факторов является метод главных компонент или компонентный анализ, который реализует модель вида

где m — количество параметров (признаков).

Каждый из наблюдаемых, параметров линейно зависит от m не коррелированных между собой новых компонент (факторов) По сравнению с моделью факторного анализа (2.19) в модели (2.28) отсутствует характерный фактор, т.е. считается, что вся вариация параметра может быть объяснена только действием общих или главных факторов. В случае компонентного анализа исходной является матрица коэффициентов корреляции, где на главной диагонали стоят единицы. Результатом компонентного анализа, так же как и факторного, является матрица факторных нагрузок. Поиск факторного решения — это ортогональное преобразование матрицы исходных переменных, в результате которого каждый параметр может быть представлен линейной комбинацией найденных m факторов, которые называют главными компонентами. Главные компоненты легко выражаются через наблюдённые параметры.

Если для дальнейшего анализа оставить все найденные т компонент, то тем самым будет использована вся информация, заложенная в корреляционной матрице. Однако это неудобно и нецелесообразно. На практике обычно оставляют небольшое число компонент, причём количество их определяется долей суммарной дисперсии, учитываемой этими компонентами. Существуют различные критерии для оценки числа оставляемых компонент; чаще всего используют следующий простой критерий: оставляют столько компонент, чтобы суммарная дисперсия, учитываемая ими, составляла заранее установленное число процентов. Первая из компонент должна учитывать максимум суммарной дисперсии параметров; вторая — не коррелировать с первой и учитывать максимум оставшейся дисперсии и так до тех пор, пока вся дисперсия не будет учтена. Сумма учтённых всеми компонентами дисперсий равна сумме дисперсий исходных параметров. Математический аппарат компонентного анализа полностью совпадает с аппаратом метода главных факторов. Отличие только в исходной матрице корреляций.

Компонента (или фактор) через исходные переменные выражается следующим образом:

где — элементы факторного решения:— исходные переменные; .— k-е собственное значение; р — количество оставленных главных
компонент.

Для иллюстрации возможностей факторного анализа покажем, как, используя метод главных компонент, можно сократить размерность пространства независимых переменных, перейдя от взаимно коррелированных параметров к независимым факторам, число которых р

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

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

Приложения

Значение t — распределения Стьюдента

Понятие о регрессионном анализе. Линейная выборочная регрессия. Метод наименьших квадратов (МНК)

Основные задачи регрессионного анализа:

•  Вычисление выборочных коэффициентов регрессии
•  Проверка значимости коэффициентов регрессии
•  Проверка адекватности модели
•  Выбор лучшей регрессии
•  Вычисление стандартных ошибок, анализ остатков

Построение простой регрессии по экспериментальным данным.

Предположим, что случайные величины  связаны линейной корреляционной зависимостью  для отыскания которой проведено  независимых измерений 

Диаграмма рассеяния (разброса, рассеивания)

 — координаты экспериментальных точек.

Выборочное уравнение прямой линии регрессии  имеет вид

Задача: подобрать  таким образом, чтобы экспериментальные точки как можно ближе лежали к прямой 

Для того, что бы провести прямую  воспользуемся МНК. Потребуем,

чтобы 

Постулаты регрессионного анализа, которые должны выполняться при использовании МНК.

1. подчинены нормальному закону распределения.
2. Дисперсия  постоянна и не зависит от номера измерения.
3. Результаты наблюдений  в разных точках независимы.
4. Входные переменные  независимы, неслучайны и измеряются без ошибок.

Введем функцию ошибок  и найдём её минимальное значение

Решив систему, получим искомые значения 

 является несмещенными оценками истинных значений коэффициентов 

 где 

 несмещенная оценка корреляционного момента (ковариации),
 несмещенная оценка дисперсии 

 выборочная ковариация,

   выборочная дисперсия 

 — выборочный коэффициент корреляции

Коэффициент детерминации

 — наблюдаемое экспериментальное значение  при 

 — предсказанное значение  удовлетворяющее уравнению регрессии

 — средневыборочное значение 

 — коэффициент детерминации, доля изменчивости  объясняемая  рассматриваемой регрессионной моделью. Для парной линейной регрессии

Коэффициент детерминации принимает значения от 0 до 1. Чем ближе значение коэффициента к 1, тем сильнее зависимость. При оценке регрессионных моделей это используется для доказательства адекватности модели (качества регрессии). Для приемлемых моделей предполагается, что коэффициент детерминации должен быть хотя бы не меньше 0,5 (в этом случае коэффициент множественной корреляции превышает по модулю 0,7). Модели с коэффициентом детерминации выше 0,8 можно признать достаточно хорошими (коэффициент корреляции превышает 0,9). Подтверждение адекватности модели проводится на основе дисперсионного анализа путем проверки гипотезы о значимости коэффициента детерминации.

 регрессия незначима

 регрессия значима

 — уровень значимости 

 — статистический критерий

Критическая область — правосторонняя; 

Если  то нулевая гипотеза отвергается на заданном уровне значимости, следовательно, коэффициент детерминации значим, следовательно, регрессия адекватна.

Мощность статистического критерия. Функция мощности

Определение. Мощностью критерия  называют вероятность попадания критерия в критическую область при условии, что справедлива конкурирующая гипотеза.

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

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

Пример:

По паспортным данным автомобиля расход топлива на 100 километров составляет 10 литров. В результате измерения конструкции двигателя ожидается, что расход топлива уменьшится. Для проверки были проведены испытания 25 автомобилей с модернизированным двигателем; выборочная средняя расхода топлива по результатам испытаний составила 9,3 литра. Предполагая, что выборка получена из нормально распределенной генеральной совокупности с математическим ожиданием  и дисперсией  проверить гипотезу, утверждающую, что изменение конструкции двигателя не повлияло на расход топлива.

3) Уровень значимости 

4) Статистический критерий

5) Критическая область — левосторонняя

  следовательно  отвергается на уровне значимости 

Пример:

В условиях примера 1 предположим, что наряду с  рассматривается конкурирующая гипотеза  а критическая область задана неравенством  Найти вероятность ошибок I рода и II рода.

 автомобилей имеют меньший расход топлива)

  автомобилей, имеющих расход топлива 9л на 100 км, классифицируются как автомобили, имеющие расход 10 литров).

Определение. Пусть проверяется  — критическая область критерия с заданным уровнем значимости  Функцией мощности критерия  называется вероятность отклонения  как функция параметра  т.е.

 — ошибка 1-ого рода

 — мощность критерия

Пример:

Построить график функции мощности из примера 2 для 

 попадает в критическую область.

Пример:

Какой минимальный объем выборки следует взять в условии примера 2 для того, чтобы обеспечить 

Лемма Неймана-Пирсона.

При проверке простой гипотезы  против простой альтернативной гипотезы  наилучшая критическая область (НКО) критерия заданного уровня значимости  состоит из точек выборочного пространства (выборок объема  для которых справедливо неравенство:

 — константа, зависящая от

 — элементы выборки;

 — функция правдоподобия при условии, что соответствующая гипотеза верна.

Пример:

Случайная величина  имеет нормальное распределение с параметрами  известно. Найти НКО для проверки  против причем 

Решение:

Ошибка первого рода: 

НКО: 

Пример:

Для зависимости заданной корреляционной табл. 13, найти оценки параметров  уравнения линейной регрессии  остаточную дисперсию; выяснить значимость уравнения регрессии при 

Решение. Воспользуемся предыдущими результатами

Согласно формуле (24), уравнение регрессии будет иметь вид  тогда 

Для выяснения значимости уравнения регрессии вычислим суммы Составим расчетную таблицу:

Из (27) и (28) по данным таблицы получим 

 по табл. П7 находим  

Вычислим статистику

Так как  то уравнение регрессии значимо. Остаточная дисперсия равна

Мы начнем с самых простых и понятных моделей машинного обучения: линейных. В этой главе мы разберёмся, что это такое, почему они работают и в каких случаях их стоит использовать. Так как это первый класс моделей, с которым вы столкнётесь, мы постараемся подробно проговорить все важные моменты. Заодно объясним, как работает машинное обучение, на сравнительно простых примерах.

Почему модели линейные?

Представьте, что у вас есть множество объектов $mathbb{X}$, а вы хотели бы каждому объекту сопоставить какое-то значение. К примеру, у вас есть набор операций по банковской карте, а вы бы хотели, понять, какие из этих операций сделали мошенники. Если вы разделите все операции на два класса и нулём обозначите законные действия, а единицей мошеннические, то у вас получится простейшая задача классификации. Представьте другую ситуацию: у вас есть данные геологоразведки, по которым вы хотели бы оценить перспективы разных месторождений. В данном случае по набору геологических данных ваша модель будет, к примеру, оценивать потенциальную годовую доходность шахты. Это пример задачи регрессии. Числа, которым мы хотим сопоставить объекты из нашего множества иногда называют таргетами (от английского target).

Таким образом, задачи классификации и регрессии можно сформулировать как поиск отображения из множества объектов $mathbb{X}$ в множество возможных таргетов.

Математически задачи можно описать так:

• классификация: $mathbb{X} to {0,1,ldots,K}$, где $0, ldots, K$ – номера классов,
• регрессия: $mathbb{X} to mathbb{R}$.

Очевидно, что просто сопоставить какие-то объекты каким-то числам — дело довольно бессмысленное. Мы же хотим быстро обнаруживать мошенников или принимать решение, где строить шахту. Значит нам нужен какой-то критерий качества. Мы бы хотели найти такое отображение, которое лучше всего приближает истинное соответствие между объектами и таргетами. Что значит «лучше всего» – вопрос сложный. Мы к нему будем много раз возвращаться. Однако, есть более простой вопрос: среди каких отображений мы будем искать самое лучшее? Возможных отображений может быть много, но мы можем упростить себе задачу и договориться, что хотим искать решение только в каком-то заранее заданном параметризированном семействе функций. Вся эта глава будет посвящена самому простому такому семейству — линейным функциям вида

$$
y = w_1 x_1 + ldots + w_D x_D + w_0,
$$

где $y$ – целевая переменная (таргет), $(x_1, ldots, x_D)$ – вектор, соответствующий объекту выборки (вектор признаков), а $w_1, ldots, w_D, w_0$ – параметры модели. Признаки ещё называют фичами (от английского features). Вектор $w = (w_1,ldots,w_D)$ часто называют вектором весов, так как на предсказание модели можно смотреть как на взвешенную сумму признаков объекта, а число $w_0$ – свободным коэффициентом, или сдвигом (bias). Более компактно линейную модель можно записать в виде

$$y = langle x, wrangle + w_0$$

Теперь, когда мы выбрали семейство функций, в котором будем искать решение, задача стала существенно проще. Мы теперь ищем не какое-то абстрактное отображение, а конкретный вектор $(w_0,w_1,ldots,w_D)inmathbb{R}^{D+1}$.

Замечание. Чтобы применять линейную модель, нужно, чтобы каждый объект уже был представлен вектором численных признаков $x_1,ldots,x_D$. Конечно, просто текст или граф в линейную модель не положить, придётся сначала придумать для него численные фичи. Модель называют линейной, если она является линейной по этим численным признакам.

Разберёмся, как будет работать такая модель в случае, если $D = 1$. То есть у наших объектов есть ровно один численный признак, по которому они отличаются. Теперь наша линейная модель будет выглядеть совсем просто: $y = w_1 x_1 + w_0$. Для задачи регрессии мы теперь пытаемся приблизить значение игрек какой-то линейной функцией от переменной икс. А что будет значить линейность для задачи классификации? Давайте вспомним про пример с поиском мошеннических транзакций по картам. Допустим, нам известна ровно одна численная переменная — объём транзакции. Для бинарной классификации транзакций на законные и потенциально мошеннические мы будем искать так называемое разделяющее правило: там, где значение функции положительно, мы будем предсказывать один класс, где отрицательно – другой. В нашем примере простейшим правилом будет какое-то пороговое значение объёма транзакций, после которого есть смысл пометить транзакцию как подозрительную.

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

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

Ответ (не открывайте сразу; сначала подумайте сами!)Линейные зависимости не так просты, как кажется. Пусть мы решаем задачу регрессии. Если мы подозреваем, что целевая переменная $y$ не выражается через $x_1, x_2$ как линейная функция, а зависит ещё от логарифма $x_1$ и ещё как-нибудь от того, разные ли знаки у признаков, то мы можем ввести дополнительные слагаемые в нашу линейную зависимость, просто объявим эти слагаемые новыми переменными и добавив перед ними соответствующие регрессионные коэффициенты

$$y approx w_1 x_1 + w_2 x_2 + w_3log{x_1} + w_4text{sgn}(x_1x_2) + w_0,$$

и в итоге из двумерной нелинейной задачи мы получили четырёхмерную линейную регрессию.

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

Ответ (не открывайте сразу; сначала подумайте сами!)В линейную модель можно подать только численные признаки, так что категориальную фичу придётся как-то закодировать. Рассмотрим для примера вот такой датасет

Здесь два категориальных признака – pet_type и color. Первый принимает четыре различных значения, второй – пять.

Самый простой способ – использовать one-hot кодирование (one-hot encoding). Пусть исходный признак мог принимать $M$ значений $c_1,ldots, c_M$. Давайте заменим категориальный признак на $M$ признаков, которые принимают значения $0$ и $1$: $i$-й будет отвечать на вопрос «принимает ли признак значение $c_i$?». Иными словами, вместо ячейки со значением $c_i$ у объекта появляется строка нулей и единиц, в которой единица стоит только на $i$-м месте.

В нашем примере получится вот такая табличка:

Можно было бы на этом остановиться, но добавленные признаки обладают одним неприятным свойством: в каждом из них ровно одна единица, так что сумма соответствующих столбцов равна столбцу из единиц. А это уже плохо. Представьте, что у нас есть линейная модель

$$y sim w_1x_1 + ldots + w_{D-1}x_{d-1} + w_{c_1}x_{c_1} + ldots + w_{c_M}x_{c_M} + w_0$$

Преобразуем немного правую часть:

$$ysim w_1x_1 + ldots + w_{D-1}x_{d-1} + underbrace{(w_{c_1} — w_{c_M})}_{=:w’_{c_1}}x_{c_1} + ldots + underbrace{(w_{c_{M-1}} — w_{c_M})}_{=:w’_{C_{M-1}}}x_{c_{M-1}} + w_{c_M}underbrace{(x_{c_1} + ldots + x_{c_M})}_{=1} + w_0 = $$

$$ = w_1x_1 + ldots + w_{D-1}x_{d-1} + w’_{c_1}x_{c_1} + ldots + w’_{c_{M-1}}x_{c_{M-1}} + underbrace{(w_{c_M} + w_0)}_{=w’_{0}}$$

Как видим, от одного из новых признаков можно избавиться, не меняя модель. Больше того, это стоит сделать, потому что наличие «лишних» признаков ведёт к переобучению или вовсе ломает модель – подробнее об этом мы поговорим в разделе про регуляризацию. Поэтому при использовании one-hot-encoding обычно выкидывают признак, соответствующий одному из значений. Например, в нашем примере итоговая матрица объекты-признаки будет иметь вид:

Конечно, one-hot кодирование – это самый наивный способ работы с категориальными признаками, и для более сложных фичей или фичей с большим количеством значений оно плохо подходит. С рядом более продвинутых техник вы познакомитесь в разделе про обучение представлений.

Помимо простоты, у линейных моделей есть несколько других достоинств. К примеру, мы можем достаточно легко судить, как влияют на результат те или иные признаки. Скажем, если вес $w_i$ положителен, то с ростом $i$-го признака таргет в случае регрессии будет увеличиваться, а в случае классификации наш выбор будет сдвигаться в пользу одного из классов. Значение весов тоже имеет прозрачную интерпретацию: чем вес $w_i$ больше, тем «важнее» $i$-й признак для итогового предсказания. То есть, если вы построили линейную модель, вы неплохо можете объяснить заказчику те или иные её результаты. Это качество моделей называют интерпретируемостью. Оно особенно ценится в индустриальных задачах, цена ошибки в которых высока. Если от работы вашей модели может зависеть жизнь человека, то очень важно понимать, как модель принимает те или иные решения и какими принципами руководствуется. При этом не все методы машинного обучения хорошо интерпретируемы, к примеру, поведение искусственных нейронных сетей или градиентного бустинга интерпретировать довольно сложно.

В то же время слепо доверять весам линейных моделей тоже не стоит по целому ряду причин:

• Линейные модели всё-таки довольно узкий класс функций, они неплохо работают для небольших датасетов и простых задач. Однако, если вы решаете линейной моделью более сложную задачу, то вам, скорее всего, придётся выдумывать дополнительные признаки, являющиеся сложными функциями от исходных. Поиск таких дополнительных признаков называется feature engineering, технически он устроен примерно так, как мы описали в вопросе про «полиномиальные модели». Вот только поиском таких искусственных фичей можно сильно увлечься, так что осмысленность интерпретации будет сильно зависеть от здравого смысла эксперта, строившего модель.
• Если между признаками есть приближённая линейная зависимость, коэффициенты в линейной модели могут совершенно потерять физический смысл (об этой проблеме и о том, как с ней бороться, мы поговорим дальше, когда будем обсуждать регуляризацию).
• Особенно осторожно стоит верить в утверждения вида «этот коэффициент маленький, значит, этот признак не важен». Во-первых, всё зависит от масштаба признака: вдруг коэффициент мал, чтобы скомпенсировать его. Во-вторых, зависимость действительно может быть слабой, но кто знает, в какой ситуации она окажется важна. Такие решения принимаются на основе данных, например, путём проверки статистического критерия (об этом мы коротко упомянем в разделе про вероятностные модели).
• Конкретные значения весов могут меняться в зависимости от обучающей выборки, хотя с ростом её размера они будут потихоньку сходиться к весам «наилучшей» линейной модели, которую можно было бы построить по всем-всем-всем данным на свете.

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

Линейная регрессия и метод наименьших квадратов (МНК)

Мы начнём с использования линейных моделей для решения задачи регрессии. Простейшим примером постановки задачи линейной регрессии является метод наименьших квадратов (Ordinary least squares).

Пусть у нас задан датасет $(X, y)$, где $y=(y_i)_{i=1}^N in mathbb{R}^N$ – вектор значений целевой переменной, а $X=(x_i)_{i = 1}^N in mathbb{R}^{N times D}, x_i in mathbb{R}^D$ – матрица объекты-признаки, в которой $i$-я строка – это вектор признаков $i$-го объекта выборки. Мы хотим моделировать зависимость $y_i$ от $x_i$ как линейную функцию со свободным членом. Общий вид такой функции из $mathbb{R}^D$ в $mathbb{R}$ выглядит следующим образом:

$$color{#348FEA}{f_w(x_i) = langle w, x_i rangle + w_0}$$

Свободный член $w_0$ часто опускают, потому что такого же результата можно добиться, добавив ко всем $x_i$ признак, тождественно равный единице; тогда роль свободного члена будет играть соответствующий ему вес:

$$begin{pmatrix}x_{i1} & ldots & x_{iD} end{pmatrix}cdotbegin{pmatrix}w_1\ vdots \ w_Dend{pmatrix} + w_0 =
begin{pmatrix}1 & x_{i1} & ldots & x_{iD} end{pmatrix}cdotbegin{pmatrix}w_0 \ w_1\ vdots \ w_D end{pmatrix}$$

Поскольку это сильно упрощает запись, в дальнейшем мы будем считать, что это уже сделано и зависимость имеет вид просто $f_w(x_i) = langle w, x_i rangle$.

Сведение к задаче оптимизации

Мы хотим, чтобы на нашем датасете (то есть на парах $(x_i, y_i)$ из обучающей выборки) функция $f_w$ как можно лучше приближала нашу зависимость.

Для того, чтобы чётко сформулировать задачу, нам осталось только одно: на математическом языке выразить желание «приблизить $f_w(x)$ к $y$». Говоря простым языком, мы должны научиться измерять качество модели и минимизировать её ошибку, как-то меняя обучаемые параметры. В нашем примере обучаемые параметры — это веса $w$. Функция, оценивающая то, как часто модель ошибается, традиционно называется функцией потерь, функционалом качества или просто лоссом (loss function). Важно, чтобы её было легко оптимизировать: скажем, гладкая функция потерь – это хорошо, а кусочно постоянная – просто ужасно.

Функции потерь бывают разными. От их выбора зависит то, насколько задачу в дальнейшем легко решать, и то, в каком смысле у нас получится приблизить предсказание модели к целевым значениям. Интуитивно понятно, что для нашей текущей задачи нам нужно взять вектор $y$ и вектор предсказаний модели и как-то сравнить, насколько они похожи. Так как эти вектора «живут» в одном векторном пространстве, расстояние между ними вполне может быть функцией потерь. Более того, положительная непрерывная функция от этого расстояния тоже подойдёт в качестве функции потерь. При этом способов задать расстояние между векторами тоже довольно много. От всего этого разнообразия глаза разбегаются, но мы обязательно поговорим про это позже. Сейчас давайте в качестве лосса возьмём квадрат $L^2$-нормы вектора разницы предсказаний модели и $y$. Во-первых, как мы увидим дальше, так задачу будет нетрудно решить, а во-вторых, у этого лосса есть ещё несколько дополнительных свойств:

• $L^2$-норма разницы – это евклидово расстояние $|y — f_w(x)|_2$ между вектором таргетов и вектором ответов модели, то есть мы их приближаем в смысле самого простого и понятного «расстояния».

• Как мы увидим в разделе про вероятностные модели, с точки зрения статистики это соответствует гипотезе о том, что наши данные состоят из линейного «сигнала» и нормально распределенного «шума».

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

$$L(f, X, y) = |y — f(X)|_2^2 = $$

$$= |y — Xw|_2^2 = sum_{i=1}^N(y_i — langle x_i, w rangle)^2$$

Такой функционал ошибки не очень хорош для сравнения поведения моделей на выборках разного размера. Представьте, что вы хотите понять, насколько качество модели на тестовой выборке из $2500$ объектов хуже, чем на обучающей из $5000$ объектов. Вы измерили $L^2$-норму ошибки и получили в одном случае $300$, а в другом $500$. Эти числа не очень интерпретируемы. Гораздо лучше посмотреть на среднеквадратичное отклонение

$$L(f, X, y) = frac1Nsum_{i=1}^N(y_i — langle x_i, w rangle)^2$$

По этой метрике на тестовой выборке получаем $0,12$, а на обучающей $0,1$.

Функция потерь $frac1Nsum_{i=1}^N(y_i — langle x_i, w rangle)^2$ называется Mean Squared Error, MSE или среднеквадратическим отклонением. Разница с $L^2$-нормой чисто косметическая, на алгоритм решения задачи она не влияет:

$$color{#348FEA}{text{MSE}(f, X, y) = frac{1}{N}|y — X w|_2^2}$$

В самом широком смысле, функции работают с объектами множеств: берут какой-то входящий объект из одного множества и выдают на выходе соответствующий ему объект из другого. Если мы имеем дело с отображением, которое на вход принимает функции, а на выходе выдаёт число, то такое отображение называют функционалом. Если вы посмотрите на нашу функцию потерь, то увидите, что это именно функционал. Для каждой конкретной линейной функции, которую задают веса $w_i$, мы получаем число, которое оценивает, насколько точно эта функция приближает наши значения $y$. Чем меньше это число, тем точнее наше решение, значит для того, чтобы найти лучшую модель, этот функционал нам надо минимизировать по $w$:

$$color{#348FEA}{|y — Xw|_2^2 longrightarrow min_w}$$

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

МНК: точный аналитический метод

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

Пусть $x^{(1)},ldots,x^{(D)}$ – столбцы матрицы $X$, то есть столбцы признаков. Тогда

$$Xw = w_1x^{(1)}+ldots+w_Dx^{(D)},$$

и задачу регрессии можно сформулировать следующим образом: найти линейную комбинацию столбцов $x^{(1)},ldots,x^{(D)}$, которая наилучшим способом приближает столбец $y$ по евклидовой норме – то есть найти проекцию вектора $y$ на подпространство, образованное векторами $x^{(1)},ldots,x^{(D)}$.

Разложим $y = y_{parallel} + y_{perp}$, где $y_{parallel} = Xw$ – та самая проекция, а $y_{perp}$ – ортогональная составляющая, то есть $y_{perp} = y — Xwperp x^{(1)},ldots,x^{(D)}$. Как это можно выразить в матричном виде? Оказывается, очень просто:

$$X^T(y — Xw) = 0$$

В самом деле, каждый элемент столбца $X^T(y — Xw)$ – это скалярное произведение строки $X^T$ (=столбца $X$ = одного из $x^{(i)}$) на $y — Xw$. Из уравнения $X^T(y — Xw) = 0$ уже очень легко выразить $w$:

$$w = (X^TX)^{-1}X^Ty$$

Вопрос на подумать Для вычисления $w_{ast}$ нам приходится обращать (квадратную) матрицу $X^TX$, что возможно, только если она невырожденна. Что это значит с точки зрения анализа данных? Почему мы верим, что это выполняется во всех разумных ситуациях?

Ответ (не открывайте сразу; сначала подумайте сами!)Как известно из линейной алгебры, для вещественной матрицы $X$ ранги матриц $X$ и $X^TX$ совпадают. Матрица $X^TX$ невырожденна тогда и только тогда, когда её ранг равен числу её столбцов, что равно числу столбцов матрицы $X$. Иными словами, формула регрессии поломается, только если столбцы матрицы $X$ линейно зависимы. Столбцы матрицы $X$ – это признаки. А если наши признаки линейно зависимы, то, наверное, что-то идёт не так и мы должны выкинуть часть из них, чтобы остались только линейно независимые.

Другое дело, что зачастую признаки могут быть приближённо линейно зависимы, особенно если их много. Тогда матрица $X^TX$ будет близка к вырожденной, и это, как мы дальше увидим, будет вести к разным, в том числе вычислительным проблемам.

Вычислительная сложность аналитического решения — $O(N^2D + D^3)$, где $N$ — длина выборки, $D$ — число признаков у одного объекта. Слагаемое $N^2D$ отвечает за сложность перемножения матриц $X^T$ и $X$, а слагаемое $D^3$ — за сложность обращения их произведения. Перемножать матрицы $(X^TX)^{-1}$ и $X^T$ не стоит. Гораздо лучше сначала умножить $y$ на $X^T$, а затем полученный вектор на $(X^TX)^{-1}$: так будет быстрее и, кроме того, не нужно будет хранить матрицу $(X^TX)^{-1}X^T$.

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

Проблемы «точного» решения

Заметим, что для получения ответа нам нужно обратить матрицу $X^TX$. Это создает множество проблем:

1. Основная проблема в обращении матрицы — это то, что вычислительно обращать большие матрицы дело сложное, а мы бы хотели работать с датасетами, в которых у нас могут быть миллионы точек,
2. Матрица $X^TX$, хотя почти всегда обратима в разумных задачах машинного обучения, зачастую плохо обусловлена. Особенно если признаков много, между ними может появляться приближённая линейная зависимость, которую мы можем упустить на этапе формулировки задачи. В подобных случаях погрешность нахождения $w$ будет зависеть от квадрата числа обусловленности матрицы $X$, что очень плохо. Это делает полученное таким образом решение численно неустойчивым: малые возмущения $y$ могут приводить к катастрофическим изменениям $w$.

Пара слов про число обусловленности.Пожертвовав математической строгостью, мы можем считать, что число обусловленности матрицы $X$ – это корень из отношения наибольшего и наименьшего из собственных чисел матрицы $X^TX$. Грубо говоря, оно показывает, насколько разного масштаба бывают собственные значения $X^TX$. Если рассмотреть $L^2$-норму ошибки предсказания, как функцию от $w$, то её линии уровня будут эллипсоидами, форма которых определяется квадратичной формой с матрицей $X^TX$ (проверьте это!). Таким образом, число обусловленности говорит о том, насколько вытянутыми являются эти эллипсоиды.ПодробнееДанные проблемы не являются поводом выбросить решение на помойку. Существует как минимум два способа улучшить его численные свойства, однако если вы не знаете про сингулярное разложение, то лучше вернитесь сюда, когда узнаете.

1. Построим $QR$-разложение матрицы $X$. Напомним, что это разложение, в котором матрица $Q$ ортогональна по столбцам (то есть её столбцы ортогональны и имеют длину 1; в частности, $Q^TQ=E$), а $R$ квадратная и верхнетреугольная. Подставив его в формулу, получим

   $$w = ((QR)^TQR)^{-1}(QR)^T y = (R^Tunderbrace{Q^TQ}_{=E}R)^{-1}R^TQ^Ty = R^{-1}R^{-T}R^TQ^Ty = R^{-1}Q^Ty$$

   Отметим, что написать $(R^TR)^{-1} = R^{-1}R^{-T}$ мы имеем право благодаря тому, что $R$ квадратная. Полученная формула намного проще, обращение верхнетреугольной матрицы (=решение системы с верхнетреугольной левой частью) производится быстро и хорошо, погрешность вычисления $w$ будет зависеть просто от числа обусловленности матрицы $X$, а поскольку нахождение $QR$-разложения является достаточно стабильной операцией, мы получаем решение с более хорошими, чем у исходной формулы, численными свойствами.

2. Также можно использовать псевдообратную матрицу, построенную с помощью сингулярного разложения, о котором подробно написано в разделе про матричные разложения. А именно, пусть

   $$A = Uunderbrace{mathrm{diag}(sigma_1,ldots,sigma_r)}_{=Sigma}V^T$$

   – это усечённое сингулярное разложение, где $r$ – это ранг $A$. В таком случае диагональная матрица посередине является квадратной, $U$ и $V$ ортогональны по столбцам: $U^TU = E$, $V^TV = E$. Тогда

   $$w = (VSigma underbrace{U^TU}_{=E}Sigma V^T)^{-1}VSigma U^Ty$$

   Заметим, что $VSigma^{-2}V^Tcdot VSigma^2V^T = E = VSigma^2V^Tcdot VSigma^{-2}V^T$, так что $(VSigma^2 V^T)^{-1} = VSigma^{-2}V^T$, откуда

   $$w = VSigma^{-2}underbrace{V^TV}_{=E}Sigma U^Ty = VSigma^{-1}Uy$$

   Хорошие численные свойства сингулярного разложения позволяют утверждать, что и это решение ведёт себя довольно неплохо.

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

Полностью вылечить проблемы мы не сможем, но никто и не обязывает нас останавливаться на «точном» решении (которое всё равно никогда не будет вполне точным). Поэтому ниже мы познакомим вас с совершенно другим методом.

МНК: приближенный численный метод

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

Как известно, градиент функции в точке направлен в сторону её наискорейшего роста, а антиградиент (противоположный градиенту вектор) в сторону наискорейшего убывания. То есть имея какое-то приближение оптимального значения параметра $w$, мы можем его улучшить, посчитав градиент функции потерь в точке и немного сдвинув вектор весов в направлении антиградиента:

$$w_j mapsto w_j — alpha frac{d}{d{w_j}} L(f_w, X, y) $$

где $alpha$ – это параметр алгоритма («темп обучения»), который контролирует величину шага в направлении антиградиента. Описанный алгоритм называется градиентным спуском.

Посмотрим, как будет выглядеть градиентный спуск для функции потерь $L(f_w, X, y) = frac1Nvertvert Xw — yvertvert^2$. Градиент квадрата евклидовой нормы мы уже считали; соответственно,

$$
nabla_wL = frac2{N} X^T (Xw — y)
$$

Следовательно, стартовав из какого-то начального приближения, мы можем итеративно уменьшать значение функции, пока не сойдёмся (по крайней мере в теории) к минимуму (вообще говоря, локальному, но в данном случае глобальному).

Алгоритм градиентного спуска

w = random_normal()             # можно пробовать и другие виды инициализации
repeat S times:                 # другой вариант: while abs(err) > tolerance
   f = X.dot(w)                 # посчитать предсказание
   err = f - y                  # посчитать ошибку
   grad = 2 * X.T.dot(err) / N  # посчитать градиент
   w -= alpha * grad            # обновить веса

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

• Поскольку задача выпуклая, выбор начальной точки влияет на скорость сходимости, но не настолько сильно, чтобы на практике нельзя было стартовать всегда из нуля или из любой другой приятной вам точки;
• Число обусловленности матрицы $X$ существенно влияет на скорость сходимости градиентного спуска: чем более вытянуты эллипсоиды уровня функции потерь, тем хуже;
• Темп обучения $alpha$ тоже сильно влияет на поведение градиентного спуска; вообще говоря, он является гиперпараметром алгоритма, и его, возможно, придётся подбирать отдельно. Другими гиперпараметрами являются максимальное число итераций $S$ и/или порог tolerance.

Иллюстрация.Рассмотрим три задачи регрессии, для которых матрица $X$ имеет соответственно маленькое, среднее и большое числа обусловленности. Будем строить для них модели вида $y=w_1x_1 + w_2x_2$. Раскрасим плоскость $(w_1, w_2)$ в соответствии со значениями $|X_{text{train}}w — y_{text{train}}|^2$. Тёмная область содержит минимум этой функции – оптимальное значение $w_{ast}$. Также запустим из из двух точек градиентный спуск с разными значениями темпа обучения $alpha$ и посмотрим, что получится:

Заголовки графиков («Round», «Elliptic», «Stripe-like») относятся к форме линий уровня потерь (чем более они вытянуты, тем хуже обусловлена задача и тем хуже может вести себя градиентный спуск).

Итог: при неудачном выборе $alpha$ алгоритм не сходится или идёт вразнос, а для плохо обусловленной задачи он сходится абы куда.

Вычислительная сложность градиентного спуска – $O(NDS)$, где, как и выше, $N$ – длина выборки, $D$ – число признаков у одного объекта. Сравните с оценкой $O(N^2D + D^3)$ для «наивного» вычисления аналитического решения.

Сложность по памяти – $O(ND)$ на хранение выборки. В памяти мы держим и выборку, и градиент, но в большинстве реалистичных сценариев доминирует выборка.

Стохастический градиентный спуск

На каждом шаге градиентного спуска нам требуется выполнить потенциально дорогую операцию вычисления градиента по всей выборке (сложность $O(ND)$). Возникает идея заменить градиент его оценкой на подвыборке (в английской литературе такую подвыборку обычно именуют batch или mini-batch; в русской разговорной терминологии тоже часто встречается слово батч или мини-батч).

А именно, если функция потерь имеет вид суммы по отдельным парам объект-таргет

$$L(w, X, y) = frac1Nsum_{i=1}^NL(w, x_i, y_i),$$

а градиент, соответственно, записывается в виде

$$nabla_wL(w, X, y) = frac1Nsum_{i=1}^Nnabla_wL(w, x_i, y_i),$$

то предлагается брать оценку

$$nabla_wL(w, X, y) approx frac1Bsum_{t=1}^Bnabla_wL(w, x_{i_t}, y_{i_t})$$

для некоторого подмножества этих пар $(x_{i_t}, y_{i_t})_{t=1}^B$. Обратите внимание на множители $frac1N$ и $frac1B$ перед суммами. Почему они нужны? Полный градиент $nabla_wL(w, X, y)$ можно воспринимать как среднее градиентов по всем объектам, то есть как оценку матожидания $mathbb{E}nabla_wL(w, x, y)$; тогда, конечно, оценка матожидания по меньшей подвыборке тоже будет иметь вид среднего градиентов по объектам этой подвыборки.

Как делить выборку на батчи? Ясно, что можно было бы случайным образом сэмплировать их из полного датасета, но даже если использовать быстрый алгоритм вроде резервуарного сэмплирования, сложность этой операции не самая оптимальная. Поэтому используют линейный проход по выборке (которую перед этим лучше всё-таки случайным образом перемешать). Давайте введём ещё один параметр нашего алгоритма: размер батча, который мы обозначим $B$. Теперь на $B$ очередных примерах вычислим градиент и обновим веса модели. При этом вместо количества шагов алгоритма обычно задают количество эпох $E$. Это ещё один гиперпараметр. Одна эпоха – это один полный проход нашего сэмплера по выборке. Заметим, что если выборка очень большая, а модель компактная, то даже первый проход бывает можно не заканчивать.

Алгоритм:

 w = normal(0, 1)
 repeat E times:
   for i = B, i <= n, i += B
      X_batch = X[i-B : i]       
      y_batch = y[i-B : i]
      f = X_batch.dot(w)                 # посчитать предсказание
      err = f - y_batch                  # посчитать ошибку
      grad = 2 * X_batch.T.dot(err) / B  # посчитать градиент
      w -= alpha * grad

Сложность по времени – $O(NDE)$. На первый взгляд, она такая же, как и у обычного градиентного спуска, но заметим, что мы сделали в $N / B$ раз больше шагов, то есть веса модели претерпели намного больше обновлений.

Сложность по памяти можно довести до $O(BD)$: ведь теперь всю выборку не надо держать в памяти, а достаточно загружать лишь текущий батч (а остальная выборка может лежать на диске, что удобно, так как в реальности задачи, в которых выборка целиком не влезает в оперативную память, встречаются сплошь и рядом). Заметим, впрочем, что при этом лучше бы $B$ взять побольше: ведь чтение с диска – намного более затратная по времени операция, чем чтение из оперативной памяти.

В целом, разницу между алгоритмами можно представлять как-то так:

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

Существует определённая терминологическая путаница, иногда стохастическим градиентным спуском называют версию алгоритма, в которой размер батча равен единице (то есть максимально шумная и быстрая версия алгоритма), а версии с бОльшим размером батча называют batch gradient descent. В книгах, которые, возможно, старше вас, такая процедура иногда ещё называется incremental gradient descent. Это не очень принципиально, но вы будьте готовы, если что.

Вопрос на подумать. Вообще говоря, если объём данных не слишком велик и позволяет это сделать, объекты лучше случайным образом перемешивать перед тем, как подавать их в алгоритм стохастического градиентного спуска. Как вам кажется, почему?

Также можно использовать различные стратегии отбора объектов. Например, чаще брать объекты, на которых ошибка больше. Какие ещё стратегии вы могли бы придумать?

Ответ (не открывайте сразу; сначала подумайте сами!)Легко представить себе ситуацию, в которой объекты как-нибудь неудачно упорядочены, скажем, по возрастанию таргета. Тогда модель будет попеременно то запоминать, что все таргеты маленькие, то – что все таргеты большие. Это может и не повлиять на качество итоговой модели, но может привести и к довольно печальным последствиям. И вообще, чем более разнообразные батчи модель увидит в процессе обучения, тем лучше.

Стратегий можно придумать много. Например, не брать объекты, на которых ошибка слишком большая (возможно, это выбросы – зачем на них учиться), или вообще не брать те, на которых ошибка достаточно мала (они «ничему не учат»). Рекомендуем, впрочем, прибегать к этим эвристикам, только если вы понимаете, зачем они вам нужны и почему есть надежда, что они помогут.

Неградиентные методы

После прочтения этой главы у вас может сложиться ощущение, что приближённые способы решения ML задач и градиентные методы – это одно и тоже, но вы будете правы в этом только на 98%. В принципе, существуют и другие способы численно решать эти задачи, но в общем случае они работают гораздо хуже, чем градиентный спуск, и не обладают таким хорошим теоретическим обоснованием. Мы не будем рассказывать про них подробно, но можете на досуге почитать, скажем, про Stepwise regression, Orthogonal matching pursuit или LARS. У LARS, кстати, есть довольно интересное свойство: он может эффективно работать на выборках, в которых число признаков больше числа примеров. С алгоритмом LARS вы можете познакомиться в главе про оптимизацию.

Регуляризация

Всегда ли решение задачи регрессии единственно? Вообще говоря, нет. Так, если в выборке два признака будут линейно зависимы (и следовательно, ранг матрицы будет меньше $D$), то гарантировано найдётся такой вектор весов $nu$ что $langlenu, x_irangle = 0 forall x_i$. В этом случае, если какой-то $w$ является решением оптимизационной задачи, то и $w + alpha nu $ тоже является решением для любого $alpha$. То есть решение не только не обязано быть уникальным, так ещё может быть сколь угодно большим по модулю. Это создаёт вычислительные трудности. Малые погрешности признаков сильно возрастают при предсказании ответа, а в градиентном спуске накапливается погрешность из-за операций со слишком большими числами.

Конечно, в жизни редко бывает так, что признаки строго линейно зависимы, а вот быть приближённо линейно зависимыми они вполне могут быть. Такая ситуация называется мультиколлинеарностью. В этом случае у нас, всё равно, возникают проблемы, близкие к описанным выше. Дело в том, что $Xnusim 0$ для вектора $nu$, состоящего из коэффициентов приближённой линейной зависимости, и, соответственно, $X^TXnuapprox 0$, то есть матрица $X^TX$ снова будет близка к вырожденной. Как и любая симметричная матрица, она диагонализуется в некотором ортонормированном базисе, и некоторые из собственных значений $lambda_i$ близки к нулю. Если вектор $X^Ty$ в выражении $(X^TX)^{-1}X^Ty$ будет близким к соответствующему собственному вектору, то он будет умножаться на $1 /{lambda_i}$, что опять же приведёт к появлению у $w$ очень больших по модулю компонент (при этом $w$ ещё и будет вычислен с большой погрешностью из-за деления на маленькое число). И, конечно же, все ошибки и весь шум, которые имелись в матрице $X$, при вычислении $ysim Xw$ будут умножаться на эти большие и неточные числа и возрастать во много-много раз, что приведёт к проблемам, от которых нас не спасёт никакое сингулярное разложение.

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

Для того, чтобы справиться с этой проблемой, задачу обычно регуляризуют, то есть добавляют к ней дополнительное ограничение на вектор весов. Это ограничение можно, как и исходный лосс, задавать по-разному, но, как правило, ничего сложнее, чем $L^1$- и $L^2$-нормы, не требуется.

Вместо исходной задачи теперь предлагается решить такую:

$$color{#348FEA}{min_w L(f, X, y) = min_w(|X w — y|_2^2 + lambda |w|^k_k )}$$

$lambda$ – это очередной параметр, а $|w|^k_k $ – это один из двух вариантов:

$$color{#348FEA}{|w|^2_2 = w^2_1 + ldots + w^2_D}$$

или

$$color{#348FEA}{|w|_1^1 = vert w_1 vert + ldots + vert w_D vert}$$

Добавка $lambda|w|^k_k$ называется регуляризационным членом или регуляризатором, а число $lambda$ – коэффициентом регуляризации.

Коэффициент $lambda$ является гиперпараметром модели и достаточно сильно влияет на качество итогового решения. Его подбирают по логарифмической шкале (скажем, от 1e-2 до 1e+2), используя для сравнения моделей с разными значениями $lambda$ дополнительную валидационную выборку. При этом качество модели с подобранным коэффициентом регуляризации уже проверяют на тестовой выборке, чтобы исключить переобучение. Более подробно о том, как нужно подбирать гиперпараметры, вы можете почитать в соответствующей главе.

Отдельно надо договориться о том, что вес $w_0$, соответствующий отступу от начала координат (то есть признаку из всех единичек), мы регуляризовать не будем, потому что это не имеет смысла: если даже все значения $y$ равномерно велики, это не должно портить качество обучения. Обычно это не отображают в формулах, но если придираться к деталям, то стоило бы написать сумму по всем весам, кроме $w_0$:

$$|w|^2_2 = sum_{color{red}{j=1}}^{D}w_j^2,$$

$$|w|_1 = sum_{color{red}{j=1}}^{D} vert w_j vert$$

В случае $L^2$-регуляризации решение задачи изменяется не очень сильно. Например, продифференцировав новый лосс по $w$, легко получить, что «точное» решение имеет вид:

$$w = (X^TX + lambda I)^{-1}X^Ty$$

Отметим, что за этой формулой стоит и понятная численная интуиция: раз матрица $X^TX$ близка к вырожденной, то обращать её сродни самоубийству. Мы лучше слегка исказим её добавкой $lambda I$, которая увеличит все собственные значения на $lambda$, отодвинув их от нуля. Да, аналитическое решение перестаёт быть «точным», но за счёт снижения численных проблем мы получим более качественное решение, чем при использовании «точной» формулы.

В свою очередь, градиент функции потерь

$$L(f_w, X, y) = |Xw — y|^2 + lambda|w|^2$$

по весам теперь выглядит так:

$$
nabla_wL(f_w, X, y) = 2X^T(Xw — y) + 2lambda w
$$

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

Вопрос на подумать. Рассмотрим стохастический градиентный спуск для $L^2$-регуляризованной линейной регрессии с батчами размера $1$. Выберите правильный вариант шага SGD:

(а) $w_imapsto w_i — 2alpha(langle w, x_jrangle — y_j)x_{ji} — frac{2alphalambda}N w_i,quad i=1,ldots,D$;

(б) $w_imapsto w_i — 2alpha(langle w, x_jrangle — y_j)x_{ji} — 2alphalambda w_i,quad i=1,ldots,D$;

(в) $w_imapsto w_i — 2alpha(langle w, x_jrangle — y_j)x_{ji} — 2lambda N w_i,quad i=1,ldots D$.

Ответ (не открывайте сразу; сначала подумайте сами!)Не регуляризованная функция потерь имеет вид $mathcal{L}(X, y, w) = frac1Nsum_{i=1}^Nmathcal{L}(x_i, y_i, w)$, и её можно воспринимать, как оценку по выборке $(x_i, y_i)_{i=1}^N$ идеальной функции потерь

$$mathcal{L}(w) = mathbb{E}_{x, y}mathcal{L}(x, y, w)$$

Регуляризационный член не зависит от выборки и добавляется отдельно:

$$mathcal{L}_{text{reg}}(w) = mathbb{E}_{x, y}mathcal{L}(x, y, w) + lambda|w|^2$$

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

$$nabla_wmathcal{L}_{text{reg}}(w) = mathbb{E}_{x, y}nabla_wmathcal{L}(x, y, w) + 2lambda w,$$

Градиент по батчу – это тоже оценка градиента идеальной функции потерь, только не на выборке $(X, y)$, а на батче $(x_{t_i}, y_{t_i})_{i=1}^B$ размера $B$. Он будет выглядеть так:

$$nabla_wmathcal{L}_{text{reg}}(w) = frac1Bsum_{i=1}^Bnabla_wmathcal{L}(x_{t_i}, y_{t_i}, w) + 2lambda w.$$

Как видите, коэффициентов, связанных с числом объектов в батче или в исходной выборке, во втором слагаемом нет. Так что верным является второй вариант. Кстати, обратите внимание, что в третьем ещё и нет коэффициента $alpha$ перед производной регуляризационного слагаемого, это тоже ошибка.

Вопрос на подумать. Распишите процедуру стохастического градиентного спуска для $L^1$-регуляризованной линейной регрессии. Как вам кажется, почему никого не волнует, что функция потерь, строго говоря, не дифференцируема?

Ответ (не открывайте сразу; сначала подумайте сами!)Распишем для случая батча размера 1:

$$w_imapsto w_i — alpha(langle w, x_jrangle — y_j)x_{ji} — frac{lambda}alpha text{sign}(w_i),quad i=1,ldots,D$$

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

Отметим, что $L^1$- и $L^2$-регуляризацию можно определять для любой функции потерь $L(w, X, y)$ (и не только в задаче регрессии, а и, например, в задаче классификации тоже). Новая функция потерь будет соответственно равна

$$widetilde{L}(w, X, y) = L(w, X, y) + lambda|w|_1$$

или

$$widetilde{L}(w, X, y) = L(w, X, y) + lambda|w|_2^2$$

Разреживание весов в $L^1$-регуляризации

$L^2$-регуляризация работает прекрасно и используется в большинстве случаев, но есть одна полезная особенность $L^1$-регуляризации: её применение приводит к тому, что у признаков, которые не оказывают большого влияния на ответ, вес в результате оптимизации получается равным $0$. Это позволяет удобным образом удалять признаки, слабо влияющие на таргет. Кроме того, это даёт возможность автоматически избавляться от признаков, которые участвуют в соотношениях приближённой линейной зависимости, соответственно, спасает от проблем, связанных с мультиколлинеарностью, о которых мы писали выше.

Не очень строгим, но довольно интуитивным образом это можно объяснить так:

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

1. Линии уровня $L^1$-нормы – это $N$-мерные октаэдры. Точки их касания с линиями уровня лосса, скорее всего, лежат на грани размерности, меньшей $N-1$, то есть как раз в области, где часть координат равна нулю:

Заметим, что данное построение говорит о том, как выглядит оптимальное решение задачи, но ничего не говорит о способе, которым это решение можно найти. На самом деле, найти такой оптимум непросто: у $L^1$ меры довольно плохая производная. Однако, способы есть. Можете на досуге прочитать, например, вот эту статью о том, как работало предсказание CTR в google в 2012 году. Там этой теме посвящается довольно много места. Кроме того, рекомендуем посмотреть про проксимальные методы в разделе этой книги про оптимизацию в ML.

Заметим также, что вообще-то оптимизация любой нормы $L_x, 0 < x leq 1$, приведёт к появлению разреженных векторов весов, просто если c $L^1$ ещё хоть как-то можно работать, то с остальными всё будет ещё сложнее.

Другие лоссы

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

MAE

Mean absolute error, абсолютная ошибка, появляется при замене $L^2$ нормы в MSE на $L^1$:

$$color{#348FEA}{MAE(y, widehat{y}) = frac1Nsum_{i=1}^N vert y_i — widehat{y}_ivert}$$

Можно заметить, что в MAE по сравнению с MSE существенно меньший вклад в ошибку будут вносить примеры, сильно удалённые от ответов модели. Дело тут в том, что в MAE мы считаем модуль расстояния, а не квадрат, соответственно, вклад больших ошибок в MSE получается существенно больше. Такая функция потерь уместна в случаях, когда вы пытаетесь обучить регрессию на данных с большим количеством выбросов в таргете.

Иначе на эту разницу можно посмотреть так: MSE приближает матожидание условного распределения $y mid x$, а MAE – медиану.

MAPE

Mean absolute percentage error, относительная ошибка.

$$MAPE(y, widehat{y}) = frac1Nsum_{i=1}^N left|frac{widehat{y}_i-y_i}{y_i}right|$$

Часто используется в задачах прогнозирования (например, погоды, загруженности дорог, кассовых сборов фильмов, цен), когда ответы могут быть различными по порядку величины, и при этом мы бы хотели верно угадать порядок, то есть мы не хотим штрафовать модель за предсказание 2000 вместо 1000 в разы сильней, чем за предсказание 2 вместо 1.

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

$$mathcal{L}(f, X, y) = sum_{i=1}^Nh_{delta}(y_i — langle w_i, xrangle),mbox{ где }h_{delta}(z) = begin{cases}
frac12z^2, |z|leqslantdelta,\
delta(|z| — frac12delta), |z| > delta
end{cases}$$

Число $delta$ является гиперпараметром. Сложная формула при $vert zvert > delta$ нужна, чтобы функция $h_{delta}(z)$ была непрерывной. Попробуйте объяснить, зачем может быть нужна такая функция потерь.

Ответ (не открывайте сразу; сначала подумайте сами!)Часто требования формулируют в духе «функция потерь должна слабее штрафовать то-то и сильней штрафовать вот это». Например, $L^2$-регуляризованный лосс штрафует за большие по модулю веса. В данном случае можно заметить, что при небольших значениях ошибки берётся просто MSE, а при больших мы начинаем штрафовать нашу модель менее сурово. Например, это может быть полезно для того, чтобы выбросы не так сильно влияли на результат обучения.

Линейная классификация

Теперь давайте поговорим про задачу классификации. Для начала будем говорить про бинарную классификацию на два класса. Обобщить эту задачу до задачи классификации на $K$ классов не составит большого труда. Пусть теперь наши таргеты $y$ кодируют принадлежность к положительному или отрицательному классу, то есть принадлежность множеству ${-1,1}$ (в этой главе договоримся именно так обозначать классы, хотя в жизни вам будут нередко встречаться и метки ${0,1}$), а $x$ – по-прежнему векторы из $mathbb{R}^D$. Мы хотим обучить линейную модель так, чтобы плоскость, которую она задаёт, как можно лучше отделяла объекты одного класса от другого.

В идеальной ситуации найдётся плоскость, которая разделит классы: положительный окажется с одной стороны от неё, а отрицательный с другой. Выборка, для которой это возможно, называется линейно разделимой. Увы, в реальной жизни такое встречается крайне редко.

Как обучить линейную модель классификации, нам ещё предстоит понять, но уже ясно, что итоговое предсказание можно будет вычислить по формуле

$$y = text{sign} langle w, x_irangle$$

Почему бы не решать, как задачу регрессии?Мы можем попробовать предсказывать числа $-1$ и $1$, минимизируя для этого, например, MSE с последующим взятием знака, но ничего хорошего не получится. Во-первых, регрессия почти не штрафует за ошибки на объектах, которые лежат близко к *разделяющей плоскости*, но не с той стороны. Во вторых, ошибкой будет считаться предсказание, например, $5$ вместо $1$, хотя нам-то на самом деле не важно, какой у числа модуль, лишь бы знак был правильным. Если визуализировать такое решение, то проблемы тоже вполне заметны:

Нам нужна прямая, которая разделяет эти точки, а не проходит через них!

Сконструируем теперь функционал ошибки так, чтобы он вышеперечисленными проблемами не обладал. Мы хотим минимизировать число ошибок классификатора, то есть

$$sum_i mathbb{I}[y_i neq sign langle w, x_irangle]longrightarrow min_w$$

Домножим обе части на $$y_i$$ и немного упростим

$$sum_i mathbb{I}[y_i langle w, x_irangle < 0]longrightarrow min_w$$

Величина $M = y_i langle w, x_irangle$ называется отступом (margin) классификатора. Такая фунция потерь называется misclassification loss. Легко видеть, что

• отступ положителен, когда $sign(y_i) = sign(langle w, x_irangle)$, то есть класс угадан верно; при этом чем больше отступ, тем больше расстояние от $x_i$ до разделяющей гиперплоскости, то есть «уверенность классификатора»;

• отступ отрицателен, когда $sign(y_i) ne sign(langle w, x_irangle)$, то есть класс угадан неверно; при этом чем больше по модулю отступ, тем более сокрушительно ошибается классификатор.

От каждого из отступов мы вычисляем функцию

$$F(M) = mathbb{I}[M < 0] = begin{cases}1, M < 0,\ 0, Mgeqslant 0end{cases}$$

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

Вопрос на подумать. Допустим, мы как-то обучили классификатор, и подавляющее большинство отступов оказались отрицательными. Правда ли нас постигла катастрофа?

Ответ (не открывайте сразу; сначала подумайте сами!)Наверное, мы что-то сделали не так, но ситуацию можно локально выправить, если предсказывать классы, противоположные тем, которые выдаёт наша модель.

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

Ответ (не открывайте сразу; сначала подумайте сами!)На первый взгляд кажется, что первая модель действительно лучше: ведь она предсказывает «увереннее», но на самом деле всё не так однозначно: во многих случаях модель, которая умеет «честно признать, что не очень уверена в ответе», может быть предпочтительней модели, которая врёт с той же непотопляемой уверенностью, что и говорит правду. В некоторых случаях лучше может оказаться модель, которая, по сути, просто отказывается от классификации на каких-то объектах.

Ошибка перцептрона

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

$$F(M) = max(0, -M)$$

Давайте запишем такой лосс с $L^2$-регуляризацией:

$$L(w, x, y) = lambdavertvert wvertvert^2_2 + sum_i max(0, -y_i langle w, x_irangle)$$

Найдём градиент:

$$
nabla_w L(w, x, y) = 2 lambda w + sum_i
begin{cases}
0, & y_i langle w, x_i rangle > 0 \
— y_i x_i, & y_i langle w, x_i rangle leq 0
end{cases}
$$

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

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

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

Hinge loss, SVM

Для таких случаев, как на картинке выше, возникает логичное желание не только найти разделяющую прямую, но и постараться провести её на одинаковом удалении от обоих классов, то есть максимизировать минимальный отступ:

Это можно сделать, слегка поменяв функцию ошибки, а именно положив её равной:

$$F(M) = max(0, 1-M)$$

$$L(w, x, y) = lambda||w||^2_2 + sum_i max(0, 1-y_i langle w, x_irangle)$$

$$
nabla_w L(w, x, y) = 2 lambda w + sum_i
begin{cases}
0, & 1 — y_i langle w, x_i rangle leq 0 \
— y_i x_i, & 1 — y_i langle w, x_i rangle > 0
end{cases}
$$

Почему же добавленная единичка приводит к желаемому результату?

Интуитивно это можно объяснить так: объекты, которые проклассифицированы правильно, но не очень «уверенно» (то есть $0 leq y_i langle w, x_irangle < 1$), продолжают вносить свой вклад в градиент и пытаются «отодвинуть» от себя разделяющую плоскость как можно дальше.

К данному выводу можно прийти и чуть более строго; для этого надо совершенно по-другому взглянуть на выражение, которое мы минимизируем. Поможет вот эта картинка:

Если мы максимизируем минимальный отступ, то надо максимизировать $frac{2}{|w|_2}$, то есть ширину полосы при условии того, что большинство объектов лежат с правильной стороны, что эквивалентно решению нашей исходной задачи:

$$lambda|w|^2_2 + sum_i max(0, 1-y_i langle w, x_irangle) longrightarrowminlimits_{w}$$

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

Итоговое положение плоскости задаётся всего несколькими обучающими примерами. Это ближайшие к плоскости правильно классифицированные объекты, которые называют опорными векторами или support vectors. Весь метод, соответственно, зовётся методом опорных векторов, или support vector machine, или сокращённо SVM. Начиная с шестидесятых годов это был сильнейший из известных методов машинного обучения. В девяностые его сменили методы, основанные на деревьях решений, которые, в свою очередь, недавно передали «пальму первенства» нейросетям.

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

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

Строгий вывод постановки задачи SVM можно прочитать тут или в лекции К.В. Воронцова.

Логистическая регрессия

В этом параграфе мы будем обозначать классы нулём и единицей.

Ещё один интересный метод появляется из желания посмотреть на классификацию как на задачу предсказания вероятностей. Хороший пример – предсказание кликов в интернете (например, в рекламе и поиске). Наличие клика в обучающем логе не означает, что, если повторить полностью условия эксперимента, пользователь обязательно кликнет по объекту опять. Скорее у объектов есть какая-то «кликабельность», то есть истинная вероятность клика по данному объекту. Клик на каждом обучающем примере является реализацией этой случайной величины, и мы считаем, что в пределе в каждой точке отношение положительных и отрицательных примеров должно сходиться к этой вероятности.

Проблема состоит в том, что вероятность, по определению, величина от 0 до 1, а простого способа обучить линейную модель так, чтобы это ограничение соблюдалось, нет. Из этой ситуации можно выйти так: научить линейную модель правильно предсказывать какой-то объект, связанный с вероятностью, но с диапазоном значений $(-infty,infty)$, и преобразовать ответы модели в вероятность. Таким объектом является logit или log odds – логарифм отношения вероятности положительного события к отрицательному $logleft(frac{p}{1-p}right)$.

Если ответом нашей модели является $logleft(frac{p}{1-p}right)$, то искомую вероятность посчитать не трудно:

$$langle w, x_irangle = logleft(frac{p}{1-p}right)$$

$$e^{langle w, x_irangle} = frac{p}{1-p}$$

$$p=frac{1}{1 + e^{-langle w, x_irangle}}$$

Функция в правой части называется сигмоидой и обозначается

$$color{#348FEA}{sigma(z) = frac1{1 + e^{-z}}}$$

Таким образом, $p = sigma(langle w, x_irangle)$

Как теперь научиться оптимизировать $w$ так, чтобы модель как можно лучше предсказывала логиты? Нужно применить метод максимума правдоподобия для распределения Бернулли. Это самое простое распределение, которое возникает, к примеру, при бросках монетки, которая орлом выпадает с вероятностью $p$. У нас только событием будет не орёл, а то, что пользователь кликнул на объект с такой вероятностью. Если хотите больше подробностей, почитайте про распределение Бернулли в теоретическом минимуме.

Правдоподобие позволяет понять, насколько вероятно получить данные значения таргета $y$ при данных $X$ и весах $w$. Оно имеет вид

$$ p(ymid X, w) =prod_i p(y_imid x_i, w) $$

и для распределения Бернулли его можно выписать следующим образом:

$$ p(ymid X, w) =prod_i p_i^{y_i} (1-p_i)^{1-y_i} $$

где $p_i$ – это вероятность, посчитанная из ответов модели. Оптимизировать произведение неудобно, хочется иметь дело с суммой, так что мы перейдём к логарифмическому правдоподобию и подставим формулу для вероятности, которую мы получили выше:

$$ ell(w, X, y) = sum_i big( y_i log(p_i) + (1-y_i)log(1-p_i) big) =$$

$$ =sum_i big( y_i log(sigma(langle w, x_i rangle)) + (1-y_i)log(1 — sigma(langle w, x_i rangle)) big) $$

Если заметить, что

$$
sigma(-z) = frac{1}{1 + e^z} = frac{e^{-z}}{e^{-z} + 1} = 1 — sigma(z),
$$

то выражение можно переписать проще:

$$
ell(w, X, y)=sum_i big( y_i log(sigma(langle w, x_i rangle)) + (1 — y_i) log(sigma(-langle w, x_i rangle)) big)
$$

Нас интересует $w$, для которого правдоподобие максимально. Чтобы получить функцию потерь, которую мы будем минимизировать, умножим его на минус один:

$$color{#348FEA}{L(w, X, y) = -sum_i big( y_i log(sigma(langle w, x_i rangle)) + (1 — y_i) log(sigma(-langle w, x_i rangle)) big)}$$

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

$$
nabla_w L(y, X, w) = -sum_i x_i big( y_i — sigma(langle w, x_i rangle)) big)
$$

Вывод формулы градиентаНам окажется полезным ещё одно свойство сигмоиды:

$$
frac{d log sigma(z)}{d z} = left( log left( frac{1}{1 + e^{-z}} right) right)’ = frac{e^{-z}}{1 + e^{-z}} = sigma(-z)
$$ $$
frac{d log sigma(-z)}{d z} = -sigma(z)
$$

Отсюда:

$$
nabla_w log sigma(langle w, x_i rangle) = sigma(-langle w, x_i rangle) x_i
$$ $$
nabla_w log sigma(-langle w, x_i rangle) = -sigma(langle w, x_i rangle) x_i
$$

и градиент оказывается равным

$$
nabla_w L(y, X, w) = -sum_i big( y_i x_i sigma(-langle w, x_i rangle) — (1 — y_i) x_i sigma(langle w, x_i rangle)) big) =
$$ $$
= -sum_i big( y_i x_i (1 — sigma(langle w, x_i rangle)) — (1 — y_i) x_i sigma(langle w, x_i rangle)) big) =
$$ $$
= -sum_i big( y_i x_i — y_i x_i sigma(langle w, x_i rangle) — x_i sigma(langle w, x_i rangle) + y_i x_i sigma(langle w, x_i rangle)) big) =
$$ $$
= -sum_i big( y_i x_i — x_i sigma(langle w, x_i rangle)) big)
$$

Предсказание модели будет вычисляться, как мы договаривались, следующим образом:

$$p=sigma(langle w, x_irangle)$$

Это вероятность положительного класса, а как от неё перейти к предсказанию самого класса? В других методах нам достаточно было посчитать знак предсказания, но теперь все наши предсказания положительные и находятся в диапазоне от 0 до 1. Что же делать? Интуитивным и не совсем (и даже совсем не) правильным является ответ «взять порог 0.5». Более корректным будет подобрать этот порог отдельно, для уже построенной регрессии минимизируя нужную вам метрику на отложенной тестовой выборке. Например, сделать так, чтобы доля положительных и отрицательных классов примерно совпадала с реальной.

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

Вопрос на подумать. Проверьте, что, если метки классов – это $pm1$, а не $0$ и $1$, то функцию потерь для логистической регрессии можно записать в более компактном виде:

$$mathcal{L}(w, X, y) = sum_{i=1}^Nlog(1 + e^{-y_ilangle w, x_irangle})$$

Вопрос на подумать. Правда ли разделяющая поверхность модели логистической регрессии является гиперплоскостью?

Ответ (не открывайте сразу; сначала подумайте сами!)Разделяющая поверхность отделяет множество точек, которым мы присваиваем класс $0$ (или $-1$), и множество точек, которым мы присваиваем класс $1$. Представляется логичным провести отсечку по какому-либо значению предсказанной вероятности. Однако, выбор этого значения — дело не очевидное. Как мы увидим в главе про калибровку классификаторов, это может быть не настоящая вероятность. Допустим, мы решили провести границу по значению $frac12$. Тогда разделяющая поверхность как раз задаётся равенством $p = frac12$, что равносильно $langle w, xrangle = 0$. А это гиперплоскость.

Вопрос на подумать. Допустим, что матрица объекты-признаки $X$ имеет полный ранг по столбцам (то есть все её столбцы линейно независимы). Верно ли, что решение задачи восстановления логистической регрессии единственно?

Ответ (не открывайте сразу; сначала подумайте сами!)В этот раз хорошего геометрического доказательства, как было для линейной регрессии, пожалуй, нет; нам придётся честно посчитать вторую производную и доказать, что она является положительно определённой. Сделаем это для случая, когда метки классов – это $pm1$. Формулы так получатся немного попроще. Напомним, что в этом случае

$$L(w, X, y) = -sum_{i=1}^Nlog(1 + e^{-y_ilangle w, x_irangle})$$

Следовательно,

$$frac{partial}{partial w_{j}}L(w, X, y) = sum_{i=1}^Nfrac{y_ix_{ij}e^{-y_ilangle w, x_irangle}}{1 + e^{-y_ilangle w, x_irangle}} = sum_{i=1}^Ny_ix_{ij}left(1 — frac1{1 + e^{-y_ilangle w, x_irangle}}right)$$

$$frac{partial^2L}{partial w_jpartial w_k}(w, X, y) = sum_{i=1}^Ny^2_ix_{ij}x_{ik}frac{e^{-y_ilangle w, x_irangle}}{(1 + e^{-y_ilangle w, x_irangle})^2} =$$

$$ = sum_{i=1}^Ny^2_ix_{ij}x_{ik}sigma(y_ilangle w, x_irangle)(1 — sigma(y_ilangle w, x_irangle))$$

Теперь заметим, что $y_i^2 = 1$ и что, если обозначить через $D$ диагональную матрицу с элементами $sigma(y_ilangle w, x_irangle)(1 — sigma(y_ilangle w, x_irangle))$ на диагонали, матрицу вторых производных можно представить в виде:

$$nabla^2L = left(frac{partial^2mathcal{L}}{partial w_jpartial w_k}right) = X^TDX$$

Так как $0 < sigma(y_ilangle w, x_irangle) < 1$, у матрицы $D$ на диагонали стоят положительные числа, из которых можно извлечь квадратные корни, представив $D$ в виде $D = D^{1/2}D^{1/2}$. В свою очередь, матрица $X$ имеет полный ранг по столбцам. Стало быть, для любого вектора приращения $une 0$ имеем

$$u^TX^TDXu = u^TX^T(D^{1/2})^TD^{1/2}Xu = vert D^{1/2}Xu vert^2 > 0$$

Таким образом, функция $L$ выпукла вниз как функция от $w$, и, соответственно, точка её экстремума непременно будет точкой минимума.

А теперь – почему это не совсем правда. Дело в том, что, говоря «точка её экстремума непременно будет точкой минимума», мы уже подразумеваем существование этой самой точки экстремума. Только вот существует этот экстремум не всегда. Можно показать, что для линейно разделимой выборки функция потерь логистической регрессии не ограничена снизу, и, соответственно, никакого экстремума нет. Доказательство мы оставляем читателю.

Вопрос на подумать. На картинке ниже представлены результаты работы на одном и том же датасете трёх моделей логистической регрессии с разными коэффициентами $L^2$-регуляризации:

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

Как вам кажется, какие картинки соответствуют самому большому коэффициенту регуляризации, а какие – самому маленькому? Почему?

Ответ (не открывайте сразу; сначала подумайте сами!)Коэффициент регуляризации максимален у левой модели. На это нас могут натолкнуть два соображения. Во-первых, разделяющая прямая проведена достаточно странно, то есть можно заподозрить, что регуляризационный член в лосс-функции перевесил функцию потерь исходной задачи. Во-вторых, модель предсказывает довольно близкие к $frac12$ вероятности – это значит, что значения $langle w, xrangle$ близки к нулю, то есть сам вектор $w$ близок к нулевому. Это также свидетельствует о том, что регуляризационный член играет слишком важную роль при оптимизации.

Наименьший коэффициент регуляризации у правой модели. Её предсказания достаточно «уверенные» (цвета на верхнем графике сочные, то есть вероятности быстро приближаются к $0$ или $1$). Это может свидетельствовать о том, что числа $langle w, xrangle$ достаточно велики по модулю, то есть $vertvert w vertvert$ достаточно велик.

Многоклассовая классификация

В этом разделе мы будем следовать изложению из лекций Евгения Соколова.

Пусть каждый объект нашей выборки относится к одному из $K$ классов: $mathbb{Y} = {1, ldots, K}$. Чтобы предсказывать эти классы с помощью линейных моделей, нам придётся свести задачу многоклассовой классификации к набору бинарных, которые мы уже хорошо умеем решать. Мы разберём два самых популярных способа это сделать – one-vs-all и all-vs-all, а проиллюстрировать их нам поможет вот такой игрушечный датасет

Один против всех (one-versus-all)

Обучим $K$ линейных классификаторов $b_1(x), ldots, b_K(x)$, выдающих оценки принадлежности классам $1, ldots, K$ соответственно. В случае с линейными моделями эти классификаторы будут иметь вид

$$b_k(x) = text{sgn}left(langle w_k, x rangle + w_{0k}right)$$

Классификатор с номером $k$ будем обучать по выборке $left(x_i, 2mathbb{I}[y_i = k] — 1right)_{i = 1}^{N}$; иными словами, мы учим классификатор отличать $k$-й класс от всех остальных.

Логично, чтобы итоговый классификатор выдавал класс, соответствующий самому уверенному из бинарных алгоритмов. Уверенность можно в каком-то смысле измерить с помощью значений линейных функций:

$$a(x) = text{argmax}_k left(langle w_k, x rangle + w_{0k}right) $$

Давайте посмотрим, что даст этот подход применительно к нашему датасету. Обучим три линейных модели, отличающих один класс от остальных:

Теперь сравним значения линейных функций

и для каждой точки выберем тот класс, которому соответствует большее значение, то есть самый «уверенный» классификатор:

Хочется сказать, что самый маленький класс «обидели».

Проблема данного подхода заключается в том, что каждый из классификаторов $b_1(x), dots, b_K(x)$ обучается на своей выборке, и значения линейных функций $langle w_k, x rangle + w_{0k}$ или, проще говоря, «выходы» классификаторов могут иметь разные масштабы. Из-за этого сравнивать их будет неправильно. Нормировать вектора весов, чтобы они выдавали ответы в одной и той же шкале, не всегда может быть разумным решением: так, в случае с SVM веса перестанут являться решением задачи, поскольку нормировка изменит норму весов.

Все против всех (all-versus-all)

Обучим $C_K^2$ классификаторов $a_{ij}(x)$, $i, j = 1, dots, K$, $i neq j$. Например, в случае с линейными моделями эти модели будут иметь вид

$$b_{ij}(x) = text{sgn}left( langle w_{ij}, x rangle + w_{0,ij} right)$$

Классификатор $a_{ij}(x)$ будем настраивать по подвыборке $X_{ij} subset X$, содержащей только объекты классов $i$ и $j$. Соответственно, классификатор $a_{ij}(x)$ будет выдавать для любого объекта либо класс $i$, либо класс $j$. Проиллюстрируем это для нашей выборки:

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

$$a(x) = text{argmax}_ksum_{i = 1}^{K} sum_{j neq i}mathbb{I}[a_{ij}(x) = k]$$

Для нашего датасета получается следующая картинка:

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

Многоклассовая логистическая регрессия

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

В логистической регрессии для двух классов мы строили линейную модель

$$b(x) = langle w, x rangle + w_0,$$

а затем переводили её прогноз в вероятность с помощью сигмоидной функции $sigma(z) = frac{1}{1 + exp(-z)}$. Допустим, что мы теперь решаем многоклассовую задачу и построили $K$ линейных моделей

$$b_k(x) = langle w_k, x rangle + w_{0k},$$

каждая из которых даёт оценку принадлежности объекта одному из классов. Как преобразовать вектор оценок $(b_1(x), ldots, b_K(x))$ в вероятности? Для этого можно воспользоваться оператором $text{softmax}(z_1, ldots, z_K)$, который производит «нормировку» вектора:

$$text{softmax}(z_1, ldots, z_K) = left(frac{exp(z_1)}{sum_{k = 1}^{K} exp(z_k)},
dots, frac{exp(z_K)}{sum_{k = 1}^{K} exp(z_k)}right).$$

В этом случае вероятность $k$-го класса будет выражаться как

$$P(y = k vert x, w) = frac{
exp{(langle w_k, x rangle + w_{0k})}}{ sum_{j = 1}^{K} exp{(langle w_j, x rangle + w_{0j})}}.$$

Обучать эти веса предлагается с помощью метода максимального правдоподобия: так же, как и в случае с двухклассовой логистической регрессией:

$$sum_{i = 1}^{N} log P(y = y_i vert x_i, w) to max_{w_1, dots, w_K}$$

Масштабируемость линейных моделей

Мы уже обсуждали, что SGD позволяет обучению хорошо масштабироваться по числу объектов, так как мы можем не загружать их целиком в оперативную память. А что делать, если признаков очень много, или мы не знаем заранее, сколько их будет? Такое может быть актуально, например, в следующих ситуациях:

• Классификация текстов: мы можем представить текст в формате «мешка слов», то есть неупорядоченного набора слов, встретившихся в данном тексте, и обучить на нём, например, определение тональности отзыва в интернете. Наличие каждого слова из языка в тексте у нас будет кодироваться отдельной фичой. Тогда размерность каждого элемента обучающей выборки будет порядка нескольких сотен тысяч.
• В задаче предсказания кликов по рекламе можно получить выборку любой размерности, например, так: в качестве фичи закодируем индикатор того, что пользователь X побывал на веб-странице Y. Суммарная размерность тогда будет порядка $10^9 cdot 10^7 = 10^{16}$. Кроме того, всё время появляются новые пользователи и веб-страницы, так что на этапе применения нас ждут сюрпризы.

Есть несколько хаков, которые позволяют бороться с такими проблемами:

• Несмотря на то, что полная размерность объекта в выборке огромна, количество ненулевых элементов в нём невелико. Значит, можно использовать разреженное кодирование, то есть вместо плотного вектора хранить словарь, в котором будут перечислены индексы и значения ненулевых элементов вектора.
• Даже хранить все веса не обязательно! Можно хранить их в хэш-таблице и вычислять индекс по формуле hash(feature) % tablesize. Хэш может вычисляться прямо от слова или id пользователя. Таким образом, несколько фичей будут иметь общий вес, который тем не менее обучится оптимальным образом. Такой подход называется hashing trick. Ясно, что сжатие вектора весов приводит к потерям в качестве, но, как правило, ценой совсем небольших потерь можно сжать этот вектор на много порядков.

Примером открытой библиотеки, в которой реализованы эти возможности, является vowpal wabbit.

Parameter server

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

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

Подытожим

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

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