Регрессия, выбросы и регуляризация

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

Одна строка поворачивает прямую

Возьмём реальный вечерний час велопроката: 60 наблюдений, температура в градусах и число поездок. Прямая наименьших квадратов даёт прирост $13{,}7$ поездки на градус. Теперь в одной строке — той, где было холодно, $+2{,}3$ °C, и город накатал скромные 159 поездок, — оператор ошибся регистром и вписал 2159. Больше ничего не менялось: те же 59 остальных строк, тот же метод. Наклон падает до $6{,}97$. Половина зависимости исчезла из-за одной ячейки.

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

Посмотрим на источник силы. Для одного наблюдения вклад в градиент равен

$$\frac{\partial}{\partial w}\,\tfrac12\bigl(y_i-w^\top x_i\bigr)^2 =-\bigl(y_i-w^\top x_i\bigr)x_i=-e_i\,x_i .$$

Тяга пропорциональна остатку $e_i$ и не имеет потолка. Удвоив ошибку в строке, мы удваиваем её голос; увеличив в сто раз — даём ей сто голосов. Наклон поворачивается ровно туда, куда тянет самый громкий крикун.

Рисунок шире экрана — проведите по немуОткрыть целиком ↗

Медиана против среднего

Разберём простейший случай: модель без признаков, одна константа $c$. Для суммы квадратов

$$\frac{d}{dc}\sum_i(y_i-c)^2=-2\sum_i(y_i-c)=0 \quad\Longrightarrow\quad c^\star=\bar y .$$

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

$$\frac{d}{dc}\sum_i|y_i-c|=\#\{y_i<c\}-\#\{y_i>c\},$$

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

Классическая иллюстрация — набор $0,0,0,0,M$. Среднее равно $M/5$ и уезжает в бесконечность вместе с $M$; медиана равна нулю при любом $M$. Долю выборки, которую оценка выдерживает без разрушения, называют точкой отказа: у среднего она нулевая, у медианы — половина.

Рисунок шире экрана — проведите по немуОткрыть целиком ↗

Хьюбер соединяет две геометрии

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

Проверим склейку в точке $e=\delta$. Значения совпадают: $\tfrac12\delta^2=\delta(\delta-\tfrac12\delta)$. Производные тоже:

$$\rho_\delta'(e)= \begin{cases} e,&|e|\leq\delta,\\ \delta\operatorname{sign}(e),&|e|>\delta \end{cases} \;=\;\operatorname{clip}(e,-\delta,\delta).$$

Производная — это и есть сила, с которой наблюдение тянет параметры. У квадрата она не ограничена, у Хьюбера упирается в потолок $\delta$. Далёкая точка не теряет голос совсем, но и не получает мегафон.

Рисунок шире экрана — проведите по немуОткрыть целиком ↗

Минимизировать $\sum_i\rho_\delta(y_i-w^\top x_i)$ удобно повторным взвешенным МНК. Приравняв градиент нулю, получаем

$$\sum_i \psi_\delta(e_i)\,x_i=0,\qquad \psi_\delta(e)=\operatorname{clip}(e,-\delta,\delta),$$

а это то же самое, что взвешенные нормальные уравнения с весами $u_i=\psi_\delta(e_i)/e_i=\min(1,\delta/|e_i|)$:

$$\bigl(X^\top U X\bigr)w=X^\top U y,\qquad U=\operatorname{diag}(u_i).$$

Веса зависят от остатков, остатки от весов, поэтому шаги повторяют до сходимости. Строка с промахом в пять порогов получает вес $0{,}2$: её не выбрасывают, её приглушают.

Порог измеряют в единицах разброса

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

Стандартное отклонение не годится: оно само построено на квадратах. Берут медиану абсолютных отклонений

$$\mathrm{MAD}=\operatorname{med}_i\bigl|e_i-\operatorname{med}_j e_j\bigr|, \qquad \hat\sigma=1{,}4826\cdot\mathrm{MAD}.$$

Множитель $1{,}4826$ подобран так, чтобы при точно нормальном шуме $\hat\sigma$ оценивала настоящее $\sigma$: для $N(0,\sigma^2)$ выполняется $\operatorname{med}|Z|=\Phi^{-1}(0{,}75)\,\sigma\approx0{,}6745\sigma$, а $1/0{,}6745\approx1{,}4826$.

На наших вечерних часах $\mathrm{MAD}=130$ поездок, значит $\hat\sigma\approx193$. Стандартный выбор $\delta=1{,}35\hat\sigma\approx261$ оставляет в линейной зоне около $13\%$ наблюдений — то есть подозрительной считается примерно каждая восьмая строка, остальные обрабатываются точно так же, как в обычном МНК.

Рисунок шире экрана — проведите по немуОткрыть целиком ↗

Выброс по ответу и рычаг по признаку

Робастная ошибка защищает от необычного $y$. Но бывает необычный $x$. Вспомним матрицу проектирования из урока 49: прогноз есть $\widehat y=Hy$, где

$$H=X\bigl(X^\top X\bigr)^{-1}X^\top,\qquad \frac{\partial\widehat y_i}{\partial y_i}=h_{ii},\qquad \sum_i h_{ii}=p .$$

Диагональ $h_{ii}$ и называют рычагом: она говорит, какую долю собственного ответа точка тянет на себя. При $p$ параметрах и $n$ строках средний рычаг равен $p/n$; строка с $h_{ii}$ в несколько раз выше среднего опасна независимо от того, какова её ошибка.

Добавим к нашим вечерним часам одну строку с температурой $45$ °C, какой в таблице нет вовсе, и правдоподобным на вид числом поездок. Её рычаг $h=0{,}14$ при медианном $0{,}025$ — почти в шесть раз больше типичного. Наклон МНК уезжает с $13{,}7$ до $16{,}4$. А Хьюбер? До $16{,}3$. Он не спас: прямая просто повернулась к новой точке, остаток стал умеренным, и приглушать оказалось нечего.

Рисунок шире экрана — проведите по немуОткрыть целиком ↗

Когда коэффициенты не определены данными

Вторая болезнь живёт в столбцах. Возьмём реальную таблицу диабета: 442 больных, десять признаков, среди них два показателя крови — общий холестерин s1 и липопротеины низкой плотности s2. Их корреляция после z-стандартизации равна $0{,}897$: почти одна и та же величина, записанная дважды.

Метод наименьших квадратов отвечает на такую таблицу вызывающе:

$$\widehat w_{s1}=-37{,}7,\qquad \widehat w_{s2}=+22{,}7 .$$

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

Чтобы понять причину, доведём ситуацию до предела. Пусть площадь дана в метрах $x_1$ и в футах $x_2=10{,}764\,x_1$. Тогда

$$\widehat y=w_1x_1+w_2x_2=(w_1+10{,}764\,w_2)\,x_1 ,$$

и любая пара с одной и той же комбинацией даёт одинаковый прогноз. Матрица $X^\top X$ вырождена, решений бесконечно много. При почти дублированных столбцах решение формально единственно, но ковариация оценок

$$\operatorname{Cov}(\widehat w)=\sigma^2\bigl(X^\top X\bigr)^{-1}$$

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

Рисунок шире экрана — проведите по немуОткрыть целиком ↗

Ridge поднимает диагональ

Добавим к ошибке штраф за длину вектора весов:

$$\widehat w_\lambda=\arg\min_w\Bigl\{\tfrac1n\|y-Xw\|_2^2+\lambda\|w\|_2^2\Bigr\}.$$

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

$$-\tfrac2n X^\top(y-Xw)+2\lambda w=0 \quad\Longrightarrow\quad \bigl(X^\top X+n\lambda I\bigr)\widehat w_\lambda=X^\top y .$$

Матрица $X^\top X$ симметрична и неотрицательно определена; прибавление $n\lambda I$ делает её строго положительно определённой, а значит обратимой при любом $\lambda>0$. Задача, у которой не было единственного решения, его получает. В библиотеках чаще пишут $\alpha=n\lambda$, и дальше мы держимся этой записи: $\alpha$ — то, что прибавляется к диагонали.

Что происходит с нашими двумя показателями крови:

$$\alpha=10:\quad w_{s1}=-11{,}3,\ w_{s2}=+1{,}8; \qquad \alpha=100:\quad w_{s1}=-2{,}1,\ w_{s2}=-3{,}7 .$$

Пара перестала воевать. При $\alpha=100$ оба коэффициента маленькие и отрицательные — теперь их можно читать как «повышенные липиды слегка связаны с худшим исходом», а не как гигантские силы разного знака.

Собственные значения и число обусловленности

Формула $(X^\top X+\alpha I)^{-1}$ становится прозрачной в собственном базисе. Пусть $X^\top X=V\Lambda V^\top$ с собственными значениями $d_j^2$ (это квадраты сингулярных чисел из урока 38). Тогда в координатах $V$

$$\widehat w_{\lambda,j}^{\,(V)} =\frac{d_j^2}{d_j^2+\alpha}\cdot\widehat w_{\mathrm{OLS},j}^{\,(V)} .$$

Каждое собственное направление сжимается своим множителем. Там, где данных много ($d_j^2\gg\alpha$), множитель близок к единице и ridge почти ничего не трогает. Там, где направление плохо освещено данными ($d_j^2\ll\alpha$), множитель близок к нулю, и штраф просто гасит шум. Это избирательное лечение, а не общее затупление модели.

Для нашей таблицы собственные значения лежат между $3{,}78$ и $1778{,}7$, то есть число обусловленности

$$\kappa\bigl(X^\top X\bigr)=\frac{d_{\max}^2}{d_{\min}^2} =\frac{1778{,}7}{3{,}78}\approx470 .$$

После добавки $\alpha$ оно становится $\kappa_\alpha=(d^2_{\max}+\alpha)/(d^2_{\min}+\alpha)$ и падает: до $129{,}8$ при $\alpha=10$, до $18{,}1$ при $\alpha=100$ и до $1{,}40$ при $\alpha=4420$ (это $\lambda=10$ в записи со средним по $n$). Геометрия из урока 36 здесь буквальна: $\kappa$ — вытянутость эллипса, в который преобразование переводит окружность, и ridge делает этот эллипс круглее.

Рисунок шире экрана — проведите по немуОткрыть целиком ↗

Ограничение вместо штрафа: форма Иванова

Тот же метод можно записать иначе — не как штраф, а как ограничение:

$$\min_w\ \tfrac1n\|y-Xw\|_2^2 \quad\text{при условии}\quad \|w\|_2\leq c .$$

Две записи связаны множителем Лагранжа из урока 23: каждому $c$ отвечает своё $\lambda$, и обратно. Но геометрически формулировка с ограничением честнее: мы прямо говорим, в каком множестве ищем ответ, и рисуем это множество — шар для $L_2$, ромб для $L_1$.

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

Рядом работала вся советская школа некорректных задач: Андрей Николаевич Тихонов, чей стабилизирующий функционал мы разбирали в уроке 23, и Михаил Михайлович Лаврентьев, изучавший условную устойчивость обратных задач геофизики. Наш $\lambda\|w\|^2$ — это тихоновский штраф, а $\|w\|\leq c$ — ивановское ограничение; одна и та же идея с двух сторон.

Lasso создаёт нули

Заменим шар ромбом, то есть $\|w\|_2^2$ на $\|w\|_1$:

$$\widehat w_\lambda=\arg\min_w \Bigl\{\tfrac1{2n}\|y-Xw\|_2^2+\lambda\sum_j|w_j|\Bigr\}.$$

Разберём один коэффициент при ортогональных столбцах. Нужно минимизировать

$$f(w)=\tfrac12(w-a)^2+\lambda|w| .$$

При $w>0$ производная равна $w-a+\lambda=0$, откуда $w=a-\lambda$ — но это годится, только если $a>\lambda$. При $w<0$ симметрично $w=a+\lambda$ при $a<-\lambda$. В нуле производной нет, есть субдифференциал $[-\lambda,\lambda]$, и ноль оптимален ровно когда $|a|\leq\lambda$. Собирая:

$$w^\star=\operatorname{sign}(a)\max\bigl(|a|-\lambda,\,0\bigr) \;=:\;S_\lambda(a).$$

Эта операция — мягкий порог. У ridge аналог был бы $a/(1+\lambda)$: множитель, который никогда не даёт ровного нуля. Разница не в силе сжатия, а в её форме: $L_1$ вычитает постоянную величину, $L_2$ делит.

Путь lasso: кто представляет группу

В уроке 23 мы уже строили путь lasso на этой же таблице и смотрели, какие признаки выживают дольше всех. Ответ был: bmi, s5 и bp — индекс массы тела, показатель липидного обмена и давление. Здесь нас занимает другой вопрос, который тогда остался за кадром: что происходит внутри группы похожих признаков и насколько устойчив выбор представителя.

Проследим пару s1/s2 вдоль пути. Первым из них в модель входит s1 при $\alpha\approx3{,}2$; s2 молчит почти до самого конца и появляется лишь при $\alpha\approx0{,}25$ — предпоследним из всех десяти признаков. При $\alpha=1$ картина такая:

$$w_{s1}=-4{,}84,\qquad w_{s2}=0,\qquad \text{ненулевых весов }7 .$$

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

Рисунок шире экрана — проведите по немуОткрыть целиком ↗

Насколько устойчив такой выбор? Пересоберите выборку — уберите половину строк или добавьте немного шума измерения, — и делегатом может стать s2, а s1 уйти в ноль. Прогноз при этом почти не изменится, а список «важных признаков» изменится полностью. Отсюда практическое правило: разреженность читают по частоте попадания признака в модель на многих подвыборках, а не по одному запуску.

Elastic net удерживает группу

Если жаль терять группу целиком, штрафы складывают:

$$\widehat w=\arg\min_w \Bigl\{\tfrac1{2n}\|y-Xw\|_2^2 +\alpha\Bigl(\gamma\|w\|_1+\tfrac{1-\gamma}{2}\|w\|_2^2\Bigr)\Bigr\}, \qquad \gamma\in[0,1].$$

Квадратичная часть делает задачу строго выпуклой, поэтому у неё единственное решение даже при точно совпадающих столбцах, а модульная часть сохраняет способность зануления. У этой смеси есть свойство группировки: чем сильнее коррелированы два признака, тем ближе друг к другу их коэффициенты. В пределе $x_1=x_2$ решение elastic net удовлетворяет $w_1=w_2$, тогда как lasso допускает любой делёж.

На нашей таблице при $\alpha=1$ и $\gamma=0{,}5$

$$w_{s1}=-0{,}24,\qquad w_{s2}=-2{,}37,\qquad \text{ненулевых весов }10 ,$$

то есть оба показателя крови остались в модели с согласованными знаками — тогда как чистый lasso при том же $\alpha$ отбросил s2 и оставил s1 с весом $-4{,}84$.

Рисунок шире экрана — проведите по немуОткрыть целиком ↗

Штраф как априорное распределение

У обоих штрафов есть вероятностное прочтение. Запишем логарифм апостериорной плотности при нормальном шуме и нормальном априорном распределении весов $w\sim N(0,\tau^2I)$:

$$\log p(w\mid y)= -\frac{1}{2\sigma^2}\|y-Xw\|_2^2 -\frac{1}{2\tau^2}\|w\|_2^2+\mathrm{const}.$$

Максимум по $w$ (оценка MAP) — это минимум суммы квадратов плюс штраф, причём

$$\lambda=\frac{\sigma^2}{\tau^2}.$$

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

Для lasso роль prior играет распределение Лапласа $p(w_j)\propto\exp(-|w_j|/b)$: у него острый пик в нуле, поэтому апостериорный максимум охотно садится ровно на ноль.

MAP — только вершина апостериорного распределения, одна точка вместо всей картины неопределённости. Что теряется при таком сжатии и как выглядит полный ответ, разбирает урок 52. Полезно помнить: выбирая $\lambda$ по кросс-валидации, мы фактически подбираем prior по данным — приём законный, но уже не вполне байесовский.

Как честно выбирать штраф

Величину $\alpha$ не выводят из формулы — её выбирают по данным, которых модель не видела. Пятикратная кросс-валидация на таблице диабета даёт

$$R^2_{\mathrm{OLS}}=0{,}489,\qquad R^2_{\mathrm{ridge},\,\alpha=100}=0{,}482,\qquad R^2_{\mathrm{lasso},\,\alpha=1}=0{,}490 .$$

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

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

$$\alpha=0{,}1:\ \mathrm{MSE}=4623,\ \mathrm{Var}=523; \qquad \alpha=63:\ \mathrm{MSE}=3829,\ \mathrm{Var}=124 .$$

Разброс упал вчетверо, смещение выросло чуть-чуть, и полная ошибка снизилась на $17\%$. Дальше, при $\alpha=10^4$, смещение берёт своё и ошибка вырастает до $5971$. Вот та самая U-образная кривая, ради которой всё и затевалось.

Рисунок шире экрана — проведите по немуОткрыть целиком ↗

Три правила честного выбора. Первое: стандартизацию признаков обучают внутри каждого train-fold, иначе среднее и разброс validation-части просачиваются в обучение. Второе: сетку $\alpha$ берут логарифмическую и достаточно широкую, чтобы минимум оказался внутри, а не на краю. Третье: если сеток и моделей перебрано много, лучший validation-результат сам становится оптимистичным, и окончательную оценку берут на отложенной test-части ровно один раз — этой ловушке посвящён урок 63.

Робастность как свойство всей процедуры

Ничто не мешает лечить обе болезни разом:

$$\min_w\ \frac1n\sum_i\rho_\delta\bigl(y_i-w^\top x_i\bigr)+\lambda\|w\|_2^2 .$$

Первое слагаемое ограничивает влияние больших остатков, второе стабилизирует плохо определённые направления. Задача остаётся выпуклой, решается тем же взвешенным МНК с добавкой к диагонали.

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

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

Лаборатория сжатия

В первой сцене — та самая таблица диабета со z-стандартизованными признаками. Поставьте штраф в минимум и убедитесь, что пара s1/s2 показывает $-37{,}7$ и $+22{,}7$; затем ведите ползунок вправо и следите не за общей картинкой, а за двумя фиолетовыми столбиками. В режиме ridge они сходятся плавно; переключите на lasso — и увидите, как счётчик ненулевых весов падает ступеньками, а s2 исчезает раньше s1. Найдите штраф, при котором в модели остаётся ровно три признака, и сравните его $R^2$ с полным.

Во второй сцене возьмитесь за жёлтую точку и тащите её вверх. Красная прямая квадратичной ошибки поедет за ней сразу; зелёная прямая Хьюбера будет держаться, пока вы не увеличите порог $\delta$. Доведите $\delta$ до четырёх $\hat\sigma$ и убедитесь, что робастность исчезла: при большом пороге Хьюбер — это тот же МНК. Затем утащите точку далеко вправо по температуре и посмотрите, как обе прямые поворачиваются вместе.

Что делать с плохо определённым

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

Есть и третья ситуация, которую легко спутать с первыми двумя: новая популяция. Если строки перестали приходить из того же распределения, ни $\rho_\delta$, ни $\lambda\|w\|^2$ не помогут, потому что менять надо не оценку, а модель мира. Умение различить «опечатка», «дублированный признак» и «другой мир» стоит дороже любой формулы сжатия — а формулы, к счастью, выводятся за полстраницы.

Задачи