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

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

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

Возьмём реальный вечерний час велопроката — все $730$ наблюдений в 17:00, температура в градусах и число поездок. Прямая наименьших квадратов даёт прирост $14{,}74$ поездки на градус. Теперь в одной строке — той, где было холодно, $+2{,}34$ °C, и город накатал скромные 159 поездок, — оператор ошибся регистром и вписал 2159. На семистах тридцати строках наклон почти не шелохнулся: $14{,}24$, потеря $3{,}4\%$. Одна ячейка тонет в толпе.

Но таблицы редко бывают такими полными. Возьмём случайные $60$ вечерних часов (seed 51) — ровно тот размер, на котором работает медицинский пример из последнего раздела урока. Там та же прямая даёт $13{,}70$, а после той же опечатки — $6{,}97$: $49\%$ зависимости исчезло из-за одной ячейки. По шести seed падение наклона лежит между $30\%$ и $72\%$ — величина эффекта зависит от того, какие 60 часов достались, но не его знак и не его порядок. Дальше в этом разделе мы работаем с этими 60 строками, потому что на них видно глазом то, что на 730 строках прячется в третьем знаке.

В уроке 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\},$$

Эта производная существует всюду, кроме самих точек $y_i$; она кусочно постоянна, скачет на каждой точке выборки и меняет знак ровно в медиане. При чётном $n$ она обращается в ноль на всём отрезке между двумя средними порядковыми значениями — там и лежит множество минимумов. При нечётном $n$ в самой медиане производной нет: слева она отрицательна, справа положительна, — и минимум достигается именно там. Здесь важна не формула, а её геометрия: модуль знает только направление промаха, а не его величину.

Классическая иллюстрация — набор $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$.

На наших 60 вечерних часах $\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\ \ \text{(при полном ранге }X).$$

Равенство $\sum_i h_{ii}=p$ — это след проектора: $\operatorname{tr}H$ равен размерности того подпространства, на которое $H$ проектирует, и она равна $p$ ровно тогда, когда столбцы $X$ линейно независимы. Диагональ $h_{ii}$ и называют рычагом: она говорит, какую долю собственного ответа точка тянет на себя. При $p$ параметрах и $n$ строках средний рычаг равен $p/n$; строка с $h_{ii}$ в несколько раз выше среднего опасна независимо от того, какова её ошибка.

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

Посмотрим честно. Остаток этой строки после подгонки равен $483$ поездки при $2\hat\sigma=387$: он второй по модулю среди всех 61 и стоит выше порога. Хьюбер её и правда приглушает — вес $\delta/|e|=261/491=0{,}53$, ровно вдвое. Но вклад строки в наклон пропорционален не только весу, а произведению $u_i\,(x_i-\bar x)^2$, и второй множитель у неё в $13{,}9$ раза больше медианного. Приглушение вдвое против рычага в четырнадцать раз — этого мало. Робастность по остатку ограничивает голос строки, но не её позицию.

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

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

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

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

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

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

Чтобы понять причину, доведём ситуацию до предела. Пусть одна и та же площадь записана дважды: $x_1$ — в квадратных метрах, $x_2$ — в квадратных футах. Фут равен $0{,}3048$ метра, поэтому $1\ \text{м}^2=1/0{,}3048^2=10{,}764$ фут², и $x_2=10{,}764\,x_1$. Тогда

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

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

$$\operatorname{Cov}(\widehat w)=\sigma^2\bigl(X^\top X\bigr)^{-1} \qquad\text{при }\operatorname{Cov}(\varepsilon)=\sigma^2I \text{ и фиксированной }X$$

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

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

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{,}70}{3{,}78}\approx470 .$$

Здесь стоит быть точным в том, что именно усиливается. Число обусловленности самой матрицы признаков равно $\kappa(X)=d_{\max}/d_{\min}=\sqrt{470}=21{,}7$, и именно оно управляет тем, как шум в ответах $y$ переходит в шум в коэффициентах. Квадрат $\kappa(X^\top X)=\kappa(X)^2=470$ появляется, когда возмущают саму матрицу $X$ (или когда решают нормальные уравнения численно, в лоб). Спутав эти два числа, читатель завышает чувствительность вдвое по логарифмической шкале.

После добавки $\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\}.$$

Разберём один коэффициент. Пусть столбцы ортогональны и нормированы так, что $X^\top X=nI$ (после z-стандартизации $\|x_j\|_2^2=n$ автоматически, остаётся потребовать ортогональность). Тогда

$$\tfrac1{2n}\|y-Xw\|_2^2 =\mathrm{const}+\sum_j\Bigl[\tfrac12 w_j^2-w_j\,\tfrac{(X^\top y)_j}{n}\Bigr] =\mathrm{const}+\sum_j\tfrac12\bigl(w_j-\widehat w_{\mathrm{OLS},j}\bigr)^2,$$

потому что при $X^\top X=nI$ решение МНК равно $(X^\top y)_j/n$. Задача распадается на $p$ независимых скалярных, и каждая из них — это

$$f(w)=\tfrac12(w-a)^2+\lambda|w|,\qquad a=\widehat w_{\mathrm{OLS},j} .$$

(Если $X^\top X=\operatorname{diag}(c_j)$ без общего множителя, задача тоже распадается, но порог у каждой координаты свой и равен $n\lambda/c_j$.)

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

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

Эта операция — мягкий порог. Ridge в той же нормировке решает $\tfrac12(w-a)^2+\lambda w^2$; приравнивая производную нулю, получаем $(w-a)+2\lambda w=0$, то есть

$$w^\star=\frac{a}{1+2\lambda}$$

— множитель, который никогда не даёт ровного нуля. (Двойка здесь не украшение: она приходит из производной $\lambda w^2$. В записи с $1/n$ и без множителя $\tfrac12$ у квадрата ошибки тот же вывод даёт $a/(1+\lambda)$ — сравнивать формулы можно только внутри одной нормировки.) Разница не в силе сжатия, а в её форме: $L_1$ вычитает постоянную величину, $L_2$ делит.

Ноль делает метод, а не сам штраф

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

У модуля $|w|$ в нуле нет производной: слева наклон $-1$, справа $+1$. Первый способ — не обращать на это внимания: в нуле взять любое число между $-1$ и $+1$ и делать обычный шаг $w_{k+1}=w_k-\alpha_k g_k$. Такой шаг попадает ровно в ноль только чудом — это всё равно что бросать камешек и надеяться, что он остановится точно на черте.

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

$$w_{k+1}=S_{\alpha\lambda}\bigl(w_k-\alpha\nabla f(w_k)\bigr), \qquad S_\kappa(v)=\operatorname{sign}(v)\cdot\max\bigl(|v|-\kappa,\ 0\bigr).$$

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

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

Путь 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].$$

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

$$|\widehat w_i-\widehat w_j|\leq \frac{\|y\|_2\sqrt{2(1-\rho_{ij})}}{\alpha(1-\gamma)},$$

и у неё есть условия: столбцы должны быть стандартизованы, а корреляция $\rho_{ij}$ — положительна. Это верхняя граница, а не обещание близости: она говорит, что при $\rho\to1$ коэффициенты обязаны совпасть, но ничего не обещает при $\rho=0{,}9$. При $\rho\to-1$ сближаются не $w_i$ и $w_j$, а $w_i$ и $-w_j$.

Чистый предел виден на точно совпадающих столбцах. Продублируем s1 в таблице: elastic net при $\alpha=1,\ \gamma=0{,}5$ делит вес поровну ($-0{,}17$ и $-0{,}17$), а lasso раскладывает как придётся — $-0{,}12$ и $-4{,}72$.

Теперь сравним на настоящей таблице, причём при равной $L_1$-силе: elastic net с $\alpha=1,\ \gamma=0{,}5$ несёт перед $\|w\|_1$ множитель $0{,}5$, значит его честный соперник — lasso с $\alpha=0{,}5$, а не с $\alpha=1$.

$$\text{EN}:\ w_{s1}=-0{,}24,\ w_{s2}=-2{,}37,\ \text{ненулевых }10; \qquad \text{lasso}_{0{,}5}:\ w_{s1}=-7{,}76,\ w_{s2}=0,\ \text{ненулевых }8 .$$

Оба показателя крови остались в модели у смеси и с согласованными знаками — и это не артефакт сравнения: lasso той же $L_1$-силы всё равно выбрасывает s2 (вместе с age). Но заметьте, чего здесь нет: коэффициенты пары не сблизились — $-0{,}24$ и $-2{,}37$ отличаются в десять раз. Группировка при $\rho=0{,}897$ означает «оба живы», а не «оба равны». Для справки: lasso с $\alpha=1$ отбросил бы s2, age и s4, оставив 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) — это минимум ненормированной суммы квадратов $\|y-Xw\|_2^2$ плюс штраф с коэффициентом $\sigma^2/\tau^2$. Чтобы не путать буквы, вспомним обозначения урока: $\alpha$ — то, что прибавляется к диагонали $X^\top X$, а $\lambda=\alpha/n$ — множитель при $\|w\|^2$ в функционале со средним $\tfrac1n\|y-Xw\|^2$. В этих обозначениях

$$\alpha=\frac{\sigma^2}{\tau^2},\qquad \lambda=\frac{\sigma^2}{n\tau^2}.$$

Разница не косметическая: при $n=442$ одно и то же априорное распределение даёт $\lambda$ в 442 раза меньше, чем $\alpha$.

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

Для 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{,}483,\qquad R^2_{\mathrm{lasso},\,\alpha=1}=0{,}490 .$$

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

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

$$\alpha=0{,}1:\ \mathrm{MSE}=4623,\ \mathrm{Var}=523,\ \text{смещение}^2{+}\text{шум}=4100;$$ $$\alpha=63:\ \mathrm{MSE}=3829,\ \mathrm{Var}=124,\ \text{смещение}^2{+}\text{шум}=3704 .$$

Здесь надо остановиться и прочитать числа, а не учебник. Разброс упал в $4{,}2$ раза — это ожидаемо. Но смещение вместе с шумом тоже упало, с $4100$ до $3704$, примерно на десятую часть. Классическая картинка «платим смещением за устойчивость» на этом участке просто не работает: при $n=60$ и $p=10$ нерегуляризованная подгонка настолько неустойчива, что даже её средний по бутстрепу прогноз промахивается сильнее, чем средний прогноз слегка стянутой модели. Штраф здесь бесплатен.

Плата приходит позже. При $\alpha=10^4$ смещение с шумом вырастает до $5867$, а полная ошибка — до $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-стандартизованными признаками. Оставьте режим ridge, уведите ползунок в самый левый край ($\log_{10}\alpha=-4$) и убедитесь, что пара s1/s2 показывает $-37{,}7$ и $+22{,}7$ — те самые числа, с которых начался разговор о коллинеарности. Затем ведите ползунок вправо и следите не за общей картинкой, а за двумя фиолетовыми столбиками. В режиме ridge они сходятся плавно; переключите на lasso — и увидите, как счётчик ненулевых весов падает ступеньками, а s2 исчезает раньше s1. Найдите штраф, при котором в модели остаётся ровно три признака ($\log_{10}\alpha\approx1{,}25$, то есть $\alpha\approx18$), и сравните его $R^2$ на обучении с полным: $0{,}392$ против $0{,}518$.

Во второй сцене возьмитесь за жёлтую точку и тащите её вверх. Красная прямая квадратичной ошибки поедет за ней сразу; зелёная прямая Хьюбера будет держаться. Теперь увеличивайте порог $\delta$ и следите за наклоном Хьюбера в показаниях. Для исходного положения точки (строка поднята на 2000) при $\delta=1{,}35\hat\sigma$ он равен $12{,}9$, при $\delta=4\hat\sigma$ — ещё $10{,}3$, и только начиная с $\delta=8\hat\sigma$ он в точности совпадает с МНК-овскими $7{,}0$: порог должен перекрыть остаток выброса, а вообще говоря совпадение с МНК — это предел $\delta\to\infty$. Утащите точку выше, и граница уедет вместе с ней. Затем утащите её далеко вправо по температуре и посмотрите, как обе прямые поворачиваются вместе.

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

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

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

Задачи