5 Введение в анализ данных и искусственный интеллект
В предыдущих главах мы учились находить корни уравнений и минимумы функций. На вид это два сюжета из вычислительной математики — однако именно они окажутся главными «инструментами в кармане», когда мы будем извлекать смысл из данных.
Эта глава вводит методологию анализа данных. Мы пройдём путь от простейшей задачи — определить долю голосов за кандидата по результатам небольшого опроса — до моделей, на которых работают современные большие языковые модели и системы распознавания изображений. Один и тот же приём — принцип максимума правдоподобия — окажется работающим и в первом, и в последнем сюжетах. Сменяются только модели: монета — нормальный шум — линейная регрессия — нейронная сеть.
Современный анализ данных уверенно стоит на трёх китах:
вероятностном — он превращает расплывчатое «извлечь информацию из данных» в строгую задачу: найти параметр модели, при котором наблюдаемые данные наиболее правдоподобны;
оптимизационном — он указывает, как именно найти этот параметр: чаще всего это сводится к минимизации некоторого функционала, и тут пригождаются методы из главы 4;
линейно-алгебраическом — он даёт язык, на котором записаны и сами данные (огромные таблицы — матрицы), и преобразования над ними (умножение матриц, малоранговые приближения, разложения).
В этой главе мы поочерёдно увидим все три кита в действии.
5.1 Простейшие примеры задач анализа данных. Принцип максимума правдоподобия
5.1.1 Опрос на втором туре выборов
На втором туре президентских выборов борьба идёт между двумя кандидатами, — назовём их условно A и B. До дня голосования социологическая служба проводит опрос: выбирает случайным образом n избирателей и спрашивает каждого, за кого тот собирается голосовать. Допустим, k человек из n ответили «за A».
Какую долю голосов в итоге наберёт кандидат A во всей стране?
Вероятностная модель. Предположим, что:
истинная (нам неизвестная!) доля сторонников A среди всех избирателей равна p\in(0,1);
опрашиваемые отбираются независимо друг от друга;
каждый отвечает честно — то есть с вероятностью p говорит «за A», с вероятностью 1-p говорит «за B».
Тогда ответ i-го опрошенного — это случайная величина X_i, принимающая значение 1 (за A) с вероятностью p и значение 0 (за B) с вероятностью 1-p. Такие величины называют бернуллиевскими, а описанную модель — схемой испытаний Бернулли.
Эксперимент, в котором n раз независимо повторяется испытание с двумя исходами («успех» с вероятностью p и «неудача» с вероятностью 1-p), называется схемой испытаний Бернулли. Число успехов S_n=X_1+\dots+X_n имеет биномиальное распределение:
\tag{5.1} \Pr(S_n=k)=\binom{n}{k}\,p^{k}(1-p)^{n-k},\qquad k=0,1,\dots,n.
Наблюдаемые данные — это вектор ответов X=(X_1,\dots,X_n) или, эквивалентно, число «единиц» в нём: S_n=X_1+\dots+X_n=k. Параметр p нам неизвестен — это и есть то, что мы хотим оценить по данным.
Идея метода максимума правдоподобия. Запишем вероятность того, что данные оказались именно такими, как мы их наблюдаем, — как функцию неизвестного параметра p:
\tag{5.2} L(p)\;=\;\Pr(X_1=x_1,\dots,X_n=x_n\mid p) \;=\;\prod_{i=1}^{n}p^{x_i}(1-p)^{1-x_i} \;=\;p^{k}(1-p)^{\,n-k},
где k=\sum_i x_i — наблюдаемое число «единиц». Функция L(p) называется функцией правдоподобия. Метод максимума правдоподобия (maximum likelihood estimation, сокращённо ОМП) состоит в выборе того \hat p, при котором L(p) достигает максимума:
\hat p\;=\;\operatorname*{arg\,max}_{p\in(0,1)} L(p).
Идея, на которой держится метод, очень житейская: из всех возможных значений параметра выбираем то, при котором наблюдаемые данные были бы наиболее правдоподобными. Если выпало 600 единиц из 1000, было бы странно объявить p=0{,}1 — ведь при таком p событие S_n=600 почти невозможно. Здравый смысл подсказывает: правдоподобнее всего гипотеза p\approx 0{,}6. Метод максимума правдоподобия и есть формализация этого здравого смысла.
Вычисление оценки. Удобнее искать максимум логарифма правдоподобия — это сводит произведение к сумме, не меняя точку максимума (логарифм монотонен):
\tag{5.3} \ln L(p)=k\,\ln p+(n-k)\,\ln(1-p).
Дифференцируем по p и приравниваем нулю (принцип Ферма):
\frac{d}{dp}\ln L(p)\;=\;\frac{k}{p}-\frac{n-k}{1-p}\;=\;0 \quad\Longleftrightarrow\quad k(1-p)=(n-k)p \quad\Longleftrightarrow\quad p=\frac{k}{n}.
Чтобы убедиться, что это именно максимум, посмотрим на вторую производную: \frac{d^{2}}{dp^{2}}\ln L(p)=-k/p^{2}-(n-k)/(1-p)^{2}<0 для всех p\in(0,1). Значит, \ln L — строго вогнутая функция, и найденная точка — её единственный максимум.
Оценкой максимума правдоподобия параметра p в схеме испытаний Бернулли по выборке X_1,\dots,X_n является выборочное среднее:
\tag{5.4} \boxed{\;\hat p \;=\; \frac{S_n}{n} \;=\; \overline{X} \;=\; \frac{X_1+\dots+X_n}{n}.\;}
Это очень житейский ответ: «доля сторонников A оценивается долей ответивших “за A”». Тем интереснее, что этот же ответ выводится из общего, абсолютно формального принципа — ОМП. На рис. 5.1 видно, что чем больше выборка, тем «острее» становится максимум функции правдоподобия — то есть тем точнее данные «локализуют» неизвестный параметр.
5.1.2 Почему ОМП хорош: теоретические результаты Фишера и Ле Кама
Естественный вопрос: является ли \hat p=\overline{X} лучшей возможной оценкой — или, может, существует ещё более точный способ извлечь p из данных?
Ответ на этот вопрос получили в первой половине XX века Р. Фишер и затем Л. Ле Кам. Изложим главные результаты на популярном уровне, без полных доказательств (они выходят за рамки школьной программы).
Рональд Эймлер Фишер (1890–1962) в работе 1922 г. «On the mathematical foundations of theoretical statistics» ввёл понятие функции правдоподобия и предложил метод её максимизации как универсальный способ оценивания. Он же ввёл термины «оценка», «параметр», «статистика» и доказал, что метод максимума правдоподобия обладает асимптотической эффективностью — его оценки имеют наименьшую возможную (по неравенству информации) дисперсию при больших размерах выборки.
Люсьен Ле Кам (1924–2000) в серии работ 1950-х—1970-х гг. построил общую теорию асимптотически нормальных статистических экспериментов и доказал, что оценка максимума правдоподобия является, в очень широком смысле, оптимальной: для гладких параметрических семейств никакая другая разумная оценка не может быть точнее.
Аналог теоремы об асимптотической эффективности в более общем бесконечномерном случае независимо изучал советский статистик Илья Александрович Ибрагимов (р. 1932) с соавторами (Ибрагимов, Хасьминский, ЛОМИ, 1970-е).
Главную мысль можно сформулировать так:
При n\to\infty оценка максимума правдоподобия \hat p_n для гладкого параметрического семейства распределений:
состоятельна: \hat p_n\to p при n\to\infty (в смысле сходимости по вероятности);
асимптотически нормальна: распределение нормированной разности \sqrt n\,(\hat p_n - p) при n\to\infty стремится к нормальному распределению с нулевым средним и определённой дисперсией;
асимптотически эффективна: эта предельная дисперсия — наименьшая среди всех «разумных» (несмещённых, регулярных) оценок. Соответственно, доверительный интервал, построенный по ОМП, оказывается кратчайшим среди всех таких интервалов, то есть наиболее точно локализует неизвестный параметр.
Иначе говоря, в очень широком классе задач ОМП — это асимптотически наилучшее из всего, на что можно надеяться при росте n. Этот результат лежит в основе огромной части современной математической статистики и машинного обучения.
5.1.3 Доверительные интервалы: Чебышёв и центральная предельная теорема
Получив точечную оценку \hat p=k/n, мы хотим понимать, насколько ей можно доверять. Социологи говорят: «46 % за кандидата A с погрешностью \pm 3\,\%». Откуда берётся это \pm 3\,\%?
Подход 1: неравенство Чебышёва
Поскольку S_n=X_1+\dots+X_n — сумма независимых одинаково распределённых случайных величин с математическим ожиданием p и дисперсией p(1-p), для \hat p_n=S_n/n:
\mathbb E[\hat p_n]=p,\qquad \mathbb D[\hat p_n]=\frac{p(1-p)}{n}.
Применим неравенство Чебышёва: для любого \varepsilon>0
\tag{5.5} \Pr\bigl(|\hat p_n - p| \geq \varepsilon\bigr)\;\leq\; \frac{\mathbb D[\hat p_n]}{\varepsilon^{2}}\;=\;\frac{p(1-p)}{n\,\varepsilon^{2}} \;\leq\;\frac{1}{4n\,\varepsilon^{2}},
где мы воспользовались тем, что p(1-p)\leq 1/4 для всех p\in[0,1].
Если потребовать, чтобы вероятность ошибки не превышала \alpha (например, \alpha=0{,}05), достаточно взять
\varepsilon\;=\;\frac{1}{2\sqrt{\alpha n}}.
Тогда с вероятностью 1-\alpha
\tag{5.6} \hat p_n - \frac{1}{2\sqrt{\alpha n}}\;\leq\; p \;\leq\; \hat p_n + \frac{1}{2\sqrt{\alpha n}}.
Подход 2: центральная предельная теорема
Неравенство Чебышёва справедливо при минимальных предположениях и потому очень грубое. Центральная предельная теорема (ЦПТ) — а её мы примем сейчас без доказательства — утверждает, что распределение нормированной суммы независимых одинаково распределённых случайных величин с конечной дисперсией стремится к нормальному.
Пусть X_1,X_2,\dots — независимые одинаково распределённые случайные величины с математическим ожиданием \mu и конечной дисперсией \sigma^2>0, S_n=X_1+\dots+X_n. Тогда при n\to\infty
\Pr\!\left(\frac{S_n-n\mu}{\sigma\sqrt n}\leq x\right)\;\longrightarrow\; \Phi(x)\;\stackrel{\mathrm{def}}{=}\;\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{x} e^{-t^{2}/2}\,dt
(функция стандартного нормального распределения).
Применим ЦПТ к нашей задаче: \mu=p, \sigma^{2}=p(1-p), и при больших n
\frac{\hat p_n - p}{\sqrt{p(1-p)/n}}\;\approx\;\mathcal N(0,1).
Если z_{1-\alpha/2} — квантиль уровня 1-\alpha/2 стандартного нормального распределения, то с вероятностью 1-\alpha
|\hat p_n-p|\;\leq\;z_{1-\alpha/2}\cdot\sqrt{p(1-p)/n}\;\leq\; \frac{z_{1-\alpha/2}}{2\sqrt n}.
Получаем доверительный интервал по ЦПТ:
\tag{5.7} \hat p_n - \frac{z_{1-\alpha/2}}{2\sqrt n}\;\leq\; p \;\leq\; \hat p_n + \frac{z_{1-\alpha/2}}{2\sqrt n}.
При \alpha=0{,}05 имеем z_{0{,}975}\approx 1{,}96.
Сравнение двух подходов
Сравним полуширины интервалов (5.6) и (5.7) при \alpha=0{,}05:
\varepsilon_\text{Чеб.}=\frac{1}{2\sqrt{0{,}05\,n}}\;=\;\frac{2{,}236}{2\sqrt n} \;=\;\frac{1{,}118}{\sqrt n}, \qquad \varepsilon_\text{ЦПТ}=\frac{1{,}96}{2\sqrt n}\;=\;\frac{0{,}98}{\sqrt n}.
Отношение \varepsilon_\text{Чеб.}/\varepsilon_\text{ЦПТ}\approx 2{,}28 — интервал Чебышёва шире более чем вдвое (рис. 5.2). Этим проиллюстрирована важная мысль: неравенство Чебышёва — это универсальная, но очень консервативная оценка; знание формы распределения (в нашем случае — асимптотической нормальности) позволяет её существенно улучшить.
Опрошено n=1000 избирателей; за кандидата A высказалось k=460. Точечная оценка — \hat p=0{,}46. Полуширина 95\,\% ДИ:
\varepsilon_\text{Чеб.}\;\approx\;\tfrac{1{,}118}{\sqrt{1000}}\;\approx\;0{,}0354, \qquad \varepsilon_\text{ЦПТ}\;\approx\;\tfrac{0{,}98}{\sqrt{1000}}\;\approx\;0{,}0310.
По ЦПТ: с уверенностью 95 % истинная доля сторонников A лежит в интервале [0{,}429,\,0{,}491], то есть точно меньше половины. Можно с большой уверенностью утверждать, что кандидат A проиграет — даже несмотря на то, что точечная оценка 0{,}46 не сильно ниже 0{,}5.
Из неравенства Чебышёва мы бы получили интервал [0{,}425,\,0{,}495] — его правый конец заметно ближе к 0{,}5, вывод об исходе выборов был бы менее уверенным.
5.1.4 Оценка площади множества: метод Монте-Карло
То, что мы сделали в задаче с выборами, имеет красивое геометрическое обобщение. Заметим: доля «успехов» k/n есть не что иное, как доля точек выборки, попавших в фиксированное множество. В нашем случае множество — это «множество всех избирателей, голосующих за A»; но точно так же мы можем оценивать, например, площадь криволинейной фигуры.
Пример: оценка числа \pi
Впишем в единичный квадрат Q=[0,1]^{2} круг D радиуса 1/2 с центром в точке (1/2,\,1/2). Площадь квадрата равна 1; площадь круга равна
\mathrm{S}(D)=\pi\cdot (1/2)^{2}=\pi/4.
Бросим в квадрат N случайных точек, равномерно распределённых по Q (для каждой точки независимо генерируем координаты x,y\sim U(0,1)). Пусть K — число точек, попавших в круг D. Тогда для каждой брошенной точки вероятность оказаться в D равна
p=\Pr\bigl((x,y)\in D\bigr)\;=\;\frac{\mathrm{S}(D)}{\mathrm{S}(Q)}\;=\;\pi/4.
Мы оказались в той же схеме испытаний Бернулли. По теореме 5.2 оценкой p является доля K/N, а оценкой числа \pi:
\tag{5.8} \boxed{\;\hat\pi_N\;=\;4\cdot\frac{K}{N}.\;}
Этот приём называется методом Монте-Карло.
Идея использовать случайные точки для приближённого вычисления определённых интегралов восходит к работам Энрико Ферми (конец 1930-х), а в современном виде была развита Станиславом Уламом, Джоном фон Нейманом и Николасом Метрополисом в Лос-Аламосе в 1946–1949 гг. при расчётах прохождения нейтронов сквозь вещество. Название «Монте-Карло» дал Метрополис — в честь известного казино, по аналогии: расчёты с розыгрышем случайных чисел напомнили коллегам игру в рулетку.
Стоит сразу отделить главное от второстепенного. На примере вычисления \pi метод Монте-Карло выглядит медленным: ошибка убывает лишь как 1/\sqrt N, и любая школьная формула численного интегрирования (трапеций, Симпсона) даёт куда быстрее куда более точный ответ. Так зачем же о нём говорить?
Сила метода — не в скорости, а в универсальности. Обычные квадратурные формулы плохо переносятся на функции многих переменных: чтобы посчитать интеграл по d-мерному кубу с точностью \varepsilon, узлов сетки требуется \sim (1/\varepsilon)^{d} — то самое «проклятие размерности», экспоненциальный рост с числом переменных. Монте-Карло же сохраняет свою скорость 1/\sqrt N независимо от размерности: для интеграла в 100 переменных он работает примерно так же хорошо (или плохо), как для одной.
Именно поэтому Монте-Карло сегодня применяется там, где размерность велика: в физике частиц (расчёты столкновений), квантовой химии (многоэлектронные интегралы), финансовой математике (опционы на много активов), байесовской статистике, а также при обучении нейронных сетей — например, при оценке градиента стохастическим спуском.
Численный эксперимент. На рис. 5.3 показан результат одного запуска при N=2000: красные точки попали вне круга, синие — внутри. По формуле (5.8), \hat\pi_{2000}\approx 3{,}154 — ошибка около 0{,}4\,\%.
Графика одной реализации, конечно, мало: при другом значении начального состояния генератора случайных чисел получится немного другой результат. Поэтому естественно изучать сходимость — как ведёт себя оценка \hat\pi_N при росте N. По формулам выше,
\mathbb E[\hat\pi_N]=\pi,\qquad \mathbb D[\hat\pi_N]=16\cdot\frac{p(1-p)}{N} =\frac{4\pi(4-\pi)}{N}\approx \frac{10{,}79}{N}.
Стандартное отклонение оценки убывает как 1/\sqrt N. На рис. 5.4 это видно отчётливо: точка пересекает горизонтальную прямую y=\pi туда-сюда, а 95\,\% «коридор» сужается как 1/\sqrt N.
Сколько точек нужно для трёх знаков после запятой? Чтобы получить \hat\pi_N с точностью \varepsilon=10^{-3} с уверенностью 95 %, нужно по (5.7), z=1{,}96: 1{,}96\cdot\sqrt{10{,}79/N}\leq 10^{-3}, откуда N\geq (1{,}96)^{2}\cdot 10{,}79\cdot 10^{6}\approx 4{,}14\cdot 10^{7}. Сорок миллионов точек ради трёх знаков — это и впрямь не самый эффективный способ вычислять \pi.
Из всех способов, доступных читателю, не выходя за рамки школьной программы, очень быстро сходится формула Мэчина (1706):
\frac{\pi}{4}\;=\;4\arctan\!\frac{1}{5}-\arctan\!\frac{1}{239}, \qquad \arctan x\;=\;x-\frac{x^{3}}{3}+\frac{x^{5}}{5}-\dots
Уже первые \sim\!10 членов рядов для \arctan(1/5) и \arctan(1/239) дают больше десяти верных знаков числа \pi. Сравним: метод Монте-Карло потребовал бы для тех же 10 знаков порядка 10^{20} точек — астрономическое число.
Формула Мэчина и её обобщения (формулы Эйлера, Гаусса, Шторма) были основным инструментом ручных вычислений \pi с XVIII по середину XX века. К началу 1970-х были найдены формулы, дающие десятки знаков на итерацию (Брент, Саламин); современные миллиарды знаков \pi считаются именно такими методами и алгоритмом БПФ — но это уже совсем другая история.
Итак, метод Монте-Карло — это плохой способ считать \pi, но отличный способ считать многомерные объёмы и интегралы. Принципиально важно: он опирается на тот же приём, что и опрос избирателей — выборочное среднее с границами, заданными ЦПТ.
5.1.5 Вторая задача: оценка скалярного параметра при гауссовом шуме
Перейдём ко второй классической задаче, в которой работает тот же принцип максимума правдоподобия. Теперь данные — не «нули и единицы», а вещественные числа с погрешностями.
В 1906 г. английский антрополог Фрэнсис Гальтон (1822–1911) посетил сельскохозяйственную выставку в Плимуте, где проходил традиционный конкурс: посетители угадывали вес выставленного быка после того, как его освежевали и взвесили. Гальтон собрал 787 заполненных бланков и обнаружил поразительный факт: среднее (а точнее, медиана) догадок участников совпало с истинным весом быка с точностью до фунта (1198 lb).
Этот эпизод — классический пример «мудрости толпы»: несмотря на то что каждый отдельный участник ошибался на десятки фунтов в обе стороны, их усреднение оказалось неожиданно точным. С точки зрения математической статистики, перед нами модель «истинное значение плюс шум», и среднее — это и есть оценка максимума правдоподобия при предположении о гауссовом шуме. Эту связь мы сейчас установим формально.
Модель. Пусть мы хотим оценить неизвестную величину \mu (вес быка; масса звезды; ускорение свободного падения — неважно). Мы располагаем n независимыми измерениями
\tag{5.9} X_i\;=\;\mu+\xi_i,\qquad i=1,\dots,n,
где \xi_i — ошибки измерения. О них естественно предположить, что:
они независимы (разные измерения не влияют друг на друга),
симметричны относительно нуля (нет систематического сдвига вверх или вниз) и
имеют одинаковое распределение.
Какое именно? Самым «канонически» правильным является нормальное (гауссовское) распределение (рис. 5.6).
Случайная величина \xi имеет нормальное распределение с параметрами \mu\in\mathbb{R} и \sigma^{2}>0, если её плотность вероятности задаётся формулой
\tag{5.10} \varphi_{\mu,\sigma^{2}}(x)\;=\;\frac{1}{\sigma\sqrt{2\pi}}\, \exp\!\left(-\frac{(x-\mu)^{2}}{2\sigma^{2}}\right),\qquad x\in\mathbb{R}.
Обозначение: \xi\sim\mathcal N(\mu,\sigma^{2}). Параметр \mu — математическое ожидание (среднее значение), \sigma^{2} — дисперсия; \sigma называют стандартным отклонением.
Принцип максимума правдоподобия для гауссова шума. Согласно модели (5.9), X_i\sim\mathcal N(\mu,\sigma^{2}) независимы. Совместная плотность выборки:
\begin{aligned} L(\mu) &=\prod_{i=1}^{n}\varphi_{\mu,\sigma^{2}}(x_i) =\frac{1}{(\sigma\sqrt{2\pi})^{n}} \exp\!\left(-\frac{1}{2\sigma^{2}}\sum_{i=1}^{n}(x_i-\mu)^{2}\right),\\[2pt] \ln L(\mu) &=-\,\frac{n}{2}\ln(2\pi\sigma^{2})\;-\;\frac{1}{2\sigma^{2}} \sum_{i=1}^{n}(x_i-\mu)^{2}. \tag{5.11} \end{aligned}
Параметр \sigma^{2} можно считать на этом этапе известным или несущественным (мы оцениваем только \mu). Тогда из всех слагаемых от \mu зависит лишь сумма квадратов:
\tag{5.12} \mu\mapsto S(\mu)\;\stackrel{\mathrm{def}}{=}\;\sum_{i=1}^{n}(x_i-\mu)^{2}.
Максимизация правдоподобия эквивалентна минимизации суммы квадратов отклонений. Это краеугольный камень всего, что последует дальше, поэтому отметим его в рамке:
\tag{5.13} \boxed{\;\;\text{ОМП при гауссовом шуме }\Longleftrightarrow\text{ метод наименьших квадратов (МНК)}\;\;}
Минимизация суммы квадратов. Найдём \hat\mu=\operatorname*{arg\,min}_{\mu} S(\mu). Дифференцируем (5.12) по \mu и приравниваем нулю:
\frac{d S}{d\mu}\;=\;-2\sum_{i=1}^{n}(x_i-\mu)\;=\;0 \quad\Longleftrightarrow\quad \sum_{i=1}^{n}x_i\;=\;n\mu.
Откуда
\tag{5.14} \boxed{\;\hat\mu\;=\;\frac{x_1+\dots+x_n}{n}\;=\;\overline{X}.\;}
Это снова выборочное среднее — теперь уже в задаче с непрерывными данными. Вторая производная S''(\mu)=2n>0, поэтому это действительно минимум, причём единственный (рис. 5.7).
5.1.6 Теорема Гаусса: симметричный шум должен быть нормальным
В рассуждении выше мы приняли нормальность шума как гипотезу. Поразительный факт — доказанный самим К. Ф. Гауссом, в чью честь и названо это распределение, — состоит в том, что эту гипотезу можно было бы и не принимать: она вытекает из требования, чтобы оценкой максимума правдоподобия было именно выборочное среднее.
Пусть случайные величины \xi_1,\dots,\xi_n независимы, одинаково распределены и имеют плотность f, симметричную относительно нуля (f(-x)=f(x)), дважды дифференцируемую, не обращающуюся в ноль. Если для любого n\geq 2 и любой реализации X_i=\mu+\xi_i оценкой максимума правдоподобия параметра \mu является выборочное среднее \overline X, то f есть плотность нормального распределения с нулевым средним.
Доказательство. Логарифмическая производная плотности \rho(x)=-f'(x)/f(x) играет роль «силы», возвращающей оценку к точке. Условие максимума ОМП
\sum_{i=1}^{n}\rho\bigl(X_i-\hat\mu\bigr)\;=\;0
должно совпадать с уравнением для среднего \sum_{i=1}^{n}(X_i-\hat\mu)=0 при любых наблюдениях. Отсюда вытекает, что \rho(x)=cx для некоторой константы c>0, то есть f'(x)/f(x)=-cx. Интегрирование даёт \ln f(x)=-cx^{2}/2+\mathrm{const}, то есть f(x) — плотность нормального распределения с нулевым средним и дисперсией 1/c. Полный вывод и обсуждение этого свойства Гаусса (характеризации нормальности через оптимальность среднего) можно найти в стандартных университетских курсах математической статистики.\square
Теорема Гаусса показывает: предположение о нормальности шума и использование среднего арифметического в качестве оценки — это одно и то же, две стороны одной медали. Если в практической задаче мы по каким-то соображениям усредняем измерения — мы неявно соглашаемся, что наши ошибки имеют гауссов характер.
В стилизованной реконструкции (рис. 5.5) среднее по n=787 ответам составило \overline X\approx 1196 фунтов при истинном весе \mu=1198 фунтов. Это не означает, что отдельный наблюдатель угадывает вес быка точно: разброс отдельных оценок составляет около \sigma\approx 75 фунтов. Однако стандартная погрешность среднего
\sigma_{\overline X}=\frac{\sigma}{\sqrt n}\;\approx\;\frac{75}{\sqrt{787}}\approx 2{,}7
фунта, и расхождение в 2 фунта вполне укладывается в одну стандартную погрешность. «Мудрость толпы» — следствие закона больших чисел, оформленного как ОМП для гауссовой модели (5.9).
5.1.7 Резюме параграфа
Что мы увидели. Две очень разные задачи — оценка доли голосов в опросе и оценка истинного значения по зашумлённым измерениям — решаются по одной и той же схеме:
Выписываем вероятностную модель порождения данных (Бернулли — в первой задаче; гауссов шум — во второй).
Записываем функцию правдоподобия L(\theta) — вероятность наблюдать имеющиеся данные при параметре \theta.
Находим точку максимума \hat\theta=\operatorname*{arg\,max} L(\theta) дифференцированием по параметру.
Строим вокруг \hat\theta доверительный интервал — с помощью неравенства Чебышёва (грубо, но универсально) или центральной предельной теоремы (точно, но при больших n).
В обеих задачах оценкой максимума правдоподобия оказалось одно и то же выражение — выборочное среднее. Это не случайно: за этим стоит теорема Гаусса (5.6), связывающая среднее арифметическое с гауссовым шумом, а в более общем плане — работы Фишера и Ле Кама об асимптотической оптимальности ОМП.
Что дальше. Возникающая задача минимизации суммы квадратов перестаёт быть тривиальной, как только величина \mu перестаёт быть скалярной и начинает зависеть от других переменных — скажем, мы хотим найти, как время падения связано с высотой (Галилей), как период обращения планеты связан с радиусом её орбиты (Кеплер) или, в общем случае, какова наилучшая линейная зависимость между y и вектором признаков x\in\mathbb{R}^{d}. Этот сюжет — линейная регрессия — будет в следующем параграфе. А затем мы шаг за шагом, не теряя нити, дойдём до глубоких нейронных сетей.
В социологическом опросе n=400, за кандидата A ответили k=220. Постройте 90\,\% доверительный интервал для доли p (а) по неравенству Чебышёва; (б) по ЦПТ. Можно ли утверждать, что A победит с вероятностью \geq 90\,\%?
Игральную кость подбросили n=600 раз. Выпало k=120 шестёрок. Найдите оценку максимума правдоподобия вероятности p выпадения шестёрки и постройте 95\,\% ДИ. Согласуется ли результат с гипотезой «кость честная»?
Запрограммируйте метод Монте-Карло для оценки числа \pi. Постройте график зависимости абсолютной ошибки |\hat\pi_N-\pi| от N в двойной логарифмической шкале и убедитесь, что наклон близок к -1/2.
Покажите, что для случайных величин X_1,\dots,X_n с плотностью Лапласа f(x;\mu)=\frac{1}{2}\exp(-|x-\mu|) оценкой максимума правдоподобия параметра \mu является выборочная медиана, а не среднее. (Указание: суммируйте -\ln f и минимизируйте.) Чем это объясняет известную «робастность» медианы по отношению к выбросам?
^{\star} Докажите теорему 5.6 полностью. Указание: подставьте X_2=\dots=X_n=0, тогда из условия \sum\rho(X_i-\hat\mu)=0 при \hat\mu=X_1/n выведите функциональное уравнение для \rho.
^{\star} Метод Монте-Карло для интеграла. Пусть требуется оценить I=\int_{0}^{1} g(x)\,dx, где g — известная функция. Покажите, что оценкой максимума правдоподобия в подходящей вероятностной модели является \hat I_N=\frac{1}{N}\sum_{i=1}^{N}g(U_i), U_i\sim U(0,1). Применяя ЦПТ, оцените сверху число точек, при котором ошибка приближённого вычисления \int_{0}^{1}e^{-x^{2}}dx не превзойдёт 10^{-3} с уверенностью 95\,\%.
5.2 Линейная регрессия и метод наименьших квадратов
В предыдущем параграфе мы оценивали скалярную величину \mu: вес быка, долю голосов, площадь множества. На практике куда чаще встречается другая постановка: мы хотим выяснить, как одна величина зависит от другой. Период обращения планеты от радиуса её орбиты. Время падения от высоты. Цена квартиры от её площади. Доход компании от числа сотрудников.
Простейшая модель такой зависимости — линейная. И именно её мы и будем сейчас изучать.
5.2.1 Постановка задачи
Дано: n пар чисел (x_1,y_1),\dots,(x_n,y_n) — результат измерений. Будем считать, что они связаны (приближённо!) линейной зависимостью
\tag{5.15} y_i\;=\;a\,x_i+b+\xi_i,\qquad i=1,\dots,n,
где
a,\,b — неизвестные коэффициенты, которые предстоит определить;
\xi_i — случайные ошибки измерений (или просто «отклонения от точной линейной зависимости», вызванные тем, что в природе идеально линейных зависимостей не бывает).
Задача линейной регрессии — найти такие \hat a,\,\hat b, чтобы прямая y=\hat a x+\hat b как можно лучше описывала наши данные.
5.2.2 Метод наименьших квадратов как ОМП при гауссовом шуме
Какой смысл вкладывать в слова «как можно лучше»? Эту неопределённость снимает принцип максимума правдоподобия из параграфа 5.1. Предположим, что ошибки \xi_i независимы и имеют нормальное распределение \mathcal N(0,\sigma^{2}). Тогда из модели (5.15) следует, что и каждое y_i нормально распределено — с математическим ожиданием a x_i+b и дисперсией \sigma^{2}. Совместная плотность выборки:
L(a,b)\;=\;\prod_{i=1}^{n}\frac{1}{\sigma\sqrt{2\pi}} \exp\!\left(-\frac{(y_i-a x_i-b)^{2}}{2\sigma^{2}}\right).
Логарифмируем и оставляем только слагаемые, зависящие от a,b:
\tag{5.16} \ln L(a,b)\;=\;C\;-\;\frac{1}{2\sigma^{2}}\sum_{i=1}^{n}(y_i-ax_i-b)^{2},
C — константа, нам не интересная. Максимизация \ln L по (a,b) эквивалентна минимизации функционала
\tag{5.17} \boxed{\;\;S(a,b)\;=\;\sum_{i=1}^{n}\bigl(y_i-a x_i-b\bigr)^{2}\longrightarrow\min_{a,b}.\;\;}
Это и есть знаменитый метод наименьших квадратов (МНК, англ. Ordinary Least Squares, OLS).
По теореме Гаусса (5.6), сумма квадратов — это единственный выбор функционала ошибки, при котором ОМП параметра \mu=a x+b соответствует гауссовскому шуму. Замените квадраты на модули |y_i-ax_i-b| — и вы будете неявно предполагать лапласовское распределение ошибок. Замените на четвёртую степень — и неявно предположите ещё более «острое» распределение. Гауссов шум привилегирован тем, что он естественно возникает в природе, когда ошибка — результат сложения большого числа независимых малых факторов.
5.2.3 Почему шум часто оказывается нормальным
Прежде чем считать формулы, разберёмся с предположением \xi_i\sim \mathcal N(0,\sigma^{2}). Откуда такая уверенность? Дело в двух свойствах нормального распределения, которые делают его «магнитом» в теории вероятностей.
Свойство 1: устойчивость относительно суммирования. Если \eta_1\sim\mathcal N(\mu_1,\sigma_1^{2}) и \eta_2\sim\mathcal N(\mu_2,\sigma_2^{2}) независимы, то их сумма снова нормальна: \eta_1+\eta_2\sim\mathcal N(\mu_1+\mu_2,\,\sigma_1^{2}+\sigma_2^{2}). Доказательство (через характеристические функции или прямой подсчёт свёртки плотностей) выходит за рамки школьной программы, но факт стоит запомнить: нормальное распределение — это «алгебраически самосогласованный» класс. Большинство других распределений после суммирования меняет форму, а нормальное — нет.
Свойство 2: ЦПТ. Гораздо более глубокий факт — центральная предельная теорема (5.3): сумма большого числа независимых случайных величин (любой природы, лишь бы с конечной дисперсией) распределена приблизительно нормально. В жизни ошибка измерения почти никогда не имеет одной-единственной природы: на неё влияют тепловые колебания прибора, вибрации, нечёткость зрения наблюдателя, погрешность шкалы, изменения освещённости и многое другое. Каждый из этих факторов даёт маленькую случайную поправку, а их сумма — по ЦПТ — получается почти нормальной (рис. 5.8).
Поэтому, когда мы пишем «шум измерения — нормальный», за этим стоят два обоснования: ЦПТ говорит, что природа благосклонна к гауссовскому шуму; а теорема Гаусса (5.6) говорит, что, выбрав МНК как метод оценки, мы автоматически предполагаем именно такой шум.
5.2.4 Решение задачи МНК. Принцип Ферма
Найдём минимум функции S(a,b) из (5.17). Это функция от двух переменных. По принципу Ферма (известному из главы 4), в точке минимума обе частные производные обращаются в ноль:
\tag{5.18} \frac{\partial S}{\partial a}\;=\;0,\qquad \frac{\partial S}{\partial b}\;=\;0.
Вычислим частные производные явно:
\begin{aligned} \frac{\partial S}{\partial a} &=\sum_{i=1}^{n}\frac{\partial}{\partial a}(y_i-a x_i-b)^{2} =\sum_{i=1}^{n}2(y_i-a x_i-b)\cdot(-x_i) =-2\sum_{i=1}^{n}x_i(y_i-a x_i-b),\\[2pt] \frac{\partial S}{\partial b} &=\sum_{i=1}^{n}2(y_i-a x_i-b)\cdot(-1) =-2\sum_{i=1}^{n}(y_i-a x_i-b). \end{aligned}
Приравнивание обеих производных нулю даёт систему нормальных уравнений:
\begin{aligned} a\sum_{i=1}^{n}x_i^{2}+b\sum_{i=1}^{n}x_i &=\sum_{i=1}^{n}x_i y_i,\\ a\sum_{i=1}^{n}x_i+b\,n &=\sum_{i=1}^{n}y_i. \end{aligned}
Обозначим выборочные средние \bar x=\frac{1}{n}\sum x_i, \bar y=\frac{1}{n}\sum y_i. Делим (5.20) на n: a\bar x+b=\bar y, откуда
\tag{5.21} \boxed{\;\hat b\;=\;\bar y-\hat a\,\bar x.\;}
Это означает, что прямая y=\hat a x+\hat b всегда проходит через «центр масс» данных (\bar x,\bar y).
Подставляем \hat b в (5.19) и после простых преобразований (умножаем (5.20) на \bar x и вычитаем из (5.19)) получаем
\tag{5.22} \boxed{\;\hat a\;=\;\dfrac{\sum_{i=1}^{n}(x_i-\bar x)(y_i-\bar y)}{\sum_{i=1}^{n}(x_i-\bar x)^{2}}.\;}
Числитель — так называемая ковариация выборок x и y (точнее, n раз ковариация), знаменатель — n раз выборочная дисперсия x. Формулы (5.21), (5.22) — основной рабочий инструмент простой линейной регрессии.
МНК минимизирует сумму квадратов вертикальных расстояний от точек данных до прямой — не перпендикулярных и не горизонтальных, а именно вертикальных (рис. 5.9). Это связано с тем, что в модели (5.15) мы считаем, что значение x_i — точное (точка на оси абсцисс), а вся неопределённость сосредоточена в y_i. Если, напротив, ошибки есть и в x, и в y, на смену МНК приходит так называемый общий МНК (Total Least Squares), минимизирующий перпендикулярные расстояния.
Возьмём шесть точек: (1,2{,}1), (2,3{,}9), (3,5{,}8), (4,8{,}3), (5,10{,}0), (6,11{,}9). Считаем по (5.21), (5.22):
\bar x=3{,}5,\quad \bar y=7{,}0,\quad \sum(x_i-\bar x)^{2}=17{,}5,\quad \sum(x_i-\bar x)(y_i-\bar y)=34{,}65.
Откуда \hat a=34{,}65/17{,}5\approx 1{,}980 и \hat b=7{,}0-1{,}98\cdot 3{,}5\approx 0{,}07. Регрессионная прямая: \hat y\approx 1{,}98\,x+0{,}07 — то есть данные «выложены вдоль» зависимости y\approx 2x.
5.2.5 Пример 1. III закон Кеплера по данным Тихо Браге
В качестве первого «настоящего» применения МНК разберём задачу, с которой началась современная физика. К XVI веку астрономы уже умели определять для каждой планеты:
период обращения T (за сколько лет планета возвращается на ту же точку звёздного неба), и
средний радиус орбиты a — большую полуось эллипса.
Тихо Браге (1546–1601), один из величайших астрономов-наблюдателей доньютоновской эпохи, посвятил всю жизнь точнейшим наблюдениям движения планет с помощью огромных угломерных инструментов на острове Вен (Дания). После его смерти эти таблицы достались его помощнику — молодому немцу Иоганну Кеплеру (1571–1630), который потратил два десятилетия на их интерпретацию.
Тихо Браге — величайший астроном-наблюдатель домикроскопной эпохи. Точность его измерений угловых положений планет составляла \sim\!1' (угловая минута) — предел для невооружённого глаза. К концу жизни он накопил каталог из почти тысячи звёзд и подробнейшие записи о движении пяти видимых планет за 21 год.
Кеплер унаследовал эти таблицы и в 1605 г. открыл первый закон движения планет: каждая планета движется по эллипсу, в одном из фокусов которого находится Солнце. В 1609 г. он добавил второй: радиус-вектор планеты заметает равные площади за равные времена. А в 1619 г. (в трактате Harmonices Mundi, «Гармонии мира») он опубликовал третий закон: квадраты периодов обращения относятся как кубы больших полуосей. Сам Кеплер выписал этот закон в форме T^{2}\sim a^{3}, не зная ни логарифмов в общеупотребительном виде, ни тем более понятия регрессии. Но именно на «эпохе Браге–Кеплера» данные впервые позволили установить количественный закон одного астрономического параметра как функции другого.
Постановка как задача регрессии. Подозреваем степенную зависимость T=C\cdot a^{\beta}. Логарифмируя обе части, получаем линейную зависимость в логарифмических осях:
\tag{5.23} \lg T\;=\;\beta\,\lg a\;+\;\lg C.
Обозначим x_i=\lg a_i, y_i=\lg T_i. Это уже стандартная задача линейной регрессии — ищем \hat\beta и \widehat{\lg C} по формулам (5.22), (5.21).
Данные. Шесть планет, известных Кеплеру (значения близки к тем, что использовал он сам, но приведены в современных единицах):
| Планета | a, а. е. | T, лет | \lg a | \lg T |
|---|---|---|---|---|
| Меркурий | 0,3871 | 0,2408 | -0{,}4122 | -0{,}6183 |
| Венера | 0,7233 | 0,6152 | -0{,}1407 | -0{,}2110 |
| Земля | 1,0000 | 1,0000 | 0 | 0 |
| Марс | 1,5237 | 1,8809 | 0,1829 | 0,2744 |
| Юпитер | 5,2034 | 11,862 | 0,7163 | 1,0742 |
| Сатурн | 9,5371 | 29,457 | 0,9794 | 1,4692 |
Расчёт МНК. Программа из десяти строк (или подсчёт «вручную» по формулам (5.22), (5.21)) даёт:
\tag{5.24} \hat\beta\;=\;1{,}4999628\dots,\qquad \widehat{\lg C}\;=\;-0{,}00003\dots,\qquad R^{2}\;=\;0{,}99999996.
Коэффициент \hat\beta совпадает с 3/2 с точностью до четвёртого знака после запятой. Это и есть III закон Кеплера:
\boxed{\;\;T^{2}\;=\;C\cdot a^{3},\qquad C\approx 1\;(\text{в единицах год/а.\,е.}^{3/2}).\;\;}
Это соотношение и графическое спрямление данных Кеплера в двойной логарифмической шкале показаны на рис. 5.10.
Кеплер не имел в распоряжении ни МНК (его открыли через 100 лет: Лежандр — в 1805 г., Гаусс — в 1809 г.), ни графиков логарифмов в современном виде. Он подбирал зависимость руками, сравнивая отношения. После многомесячных проб он заметил, что T^{2}_{\text{Юпитер}}/T^{2}_{\text{Земля}}\approx 140{,}7 и a^{3}_{\text{Юпитер}}/a^{3}_{\text{Земля}}\approx 140{,}9 — почти равенство. Сегодня та же зависимость устанавливается за десяток строк кода — но идею Кеплера это не уменьшает. Напротив, наглядно показывает, насколько мощным инструментом стал МНК для последующих поколений физиков.
5.2.6 Пример 2. Ускорение свободного падения по данным Галилея
Галилео Галилей (1564–1642) стоял у истоков физики как экспериментальной науки. Его трактат «Беседы и математические доказательства, касающиеся двух новых отраслей науки» (1638) изложил закон свободного падения: пройденный путь пропорционален квадрату времени падения,
\tag{5.25} h\;=\;\frac{g}{2}\,t^{2},
где g — ускорение свободного падения. Историки до сих пор спорят, на самом ли деле Галилей сбрасывал предметы с Пизанской башни (свидетельства его ученика Винченцо Вивиани, написанные через десятки лет, ставятся под сомнение). Но даже как «мысленный эксперимент» этот сюжет великолепно подходит, чтобы продемонстрировать МНК на ещё одной задаче (рис. 5.11).
Регрессионная модель. Записывать t через h напрямую нельзя (получится корень):
t\;=\;\sqrt{\frac{2 h}{g}}.
Поэтому регрессия выписывается для обратной зависимости — высота от квадрата времени:
\tag{5.26} h\;=\;\frac{g}{2}\,t^{2}+\xi.
Введя X=t^{2}, Y=h и A=g/2, получаем линейную модель Y=AX+\xi — частный случай (5.15) с b=0. (Свободный член можно оставить и формально — он покажет систематическую погрешность типа «реакции экспериментатора при пуске секундомера»; см. ниже.)
«Данные». Воспроизведём эксперимент в духе Discorsi. Высоты ярусов Пизанской башни (приблизительно): 7, 12, 19, 25{,}5, 32, 38, 44, 50 м над землёй. Для каждой высоты «измеряем» время падения — настоящий Галилей использовал водяные часы (взвешивал воду, вытекшую за время падения!) с типичной погрешностью около \sigma\sim 0{,}05 секунды на одно измерение. Получим набор пар (h_i,t_i):
| i | h_i, м | t_i, с (измерено) | t_i^{2}, с^{2} |
|---|---|---|---|
| 1 | 7,0 | 1,253 | 1,569 |
| 2 | 12,0 | 1,570 | 2,465 |
| 3 | 19,0 | 1,897 | 3,597 |
| 4 | 25,5 | 2,299 | 5,287 |
| 5 | 32,0 | 2,497 | 6,235 |
| 6 | 38,0 | 2,793 | 7,799 |
| 7 | 44,0 | 3,092 | 9,562 |
| 8 | 50,0 | 3,196 | 10,213 |
Расчёт МНК. Применяем формулы (5.22), (5.21) к парам (t_i^{2},h_i): \hat A\approx 4{,}757, \hat b\approx 0{,}653. Отсюда
\hat g\;=\;2\hat A\;\approx\;9{,}51\;\text{м/с}^{2},
при «истинном» значении g=9{,}81 м/с^{2}. Относительная погрешность \sim\!3\,\% — абсолютно реалистичный результат для экспериментатора XVII века с водяными часами. Коэффициент детерминации R^{2}\approx 0{,}9932 — то есть модель объясняет более 99\,\% дисперсии в данных (рис. 5.12).
И в задаче Кеплера, и в задаче Галилея исходная зависимость — не линейная (T=Ca^{3/2}, t=\sqrt{2h/g}). Прежде чем применить МНК, мы её линеаризуем: берём логарифм у Кеплера и квадрат t у Галилея. Это очень распространённый трюк: подбор показателя степени, экспоненциальный рост, степенные законы — все они сводятся к линейной регрессии заменой переменных.
Однако с этим связано важное наблюдение. После замены переменных меняется и распределение шума: гауссов шум на исходных y после взятия логарифма уже не гауссов на \lg y. Поэтому формально «линеаризованная» оценка слегка отличается от ОМП в исходной модели — но в большинстве практических задач отличие незначительно, и линеаризацией пользуются как удобным и универсальным приёмом.
5.2.7 Множественная линейная регрессия (краткий обзор)
В реальных задачах одной входной переменной обычно мало. Цена квартиры зависит одновременно от её площади, расстояния до метро, этажа, года постройки. Время на маршруте — от длины пути, числа светофоров, загруженности, времени суток.
Постановка. Пусть n — число наблюдений, d — число входных признаков. Для i-го наблюдения известны: вектор признаков \mathbf x_i=(x_{i,1},\dots,x_{i,d})\in\mathbb{R}^{d} и значение отклика y_i\in\mathbb{R}. Запишем общую множественную линейную модель:
\tag{5.27} y_i\;=\;w_1 x_{i,1}+w_2 x_{i,2}+\dots+w_d x_{i,d}+w_0+\xi_i, \qquad i=1,\dots,n.
Здесь w_1,\dots,w_d — коэффициенты при признаках, w_0 — свободный член, \xi_i\sim\mathcal N(0,\sigma^{2}) — независимый гауссов шум.
Матричная запись. Соберём матрицу плана X\in\mathbb{R}^{n\times(d+1)}, у которой i-я строка — это (x_{i,1},x_{i,2},\dots,x_{i,d},1) (последняя единица позволяет учесть свободный член w_0 как обычный коэффициент). Вектор параметров \mathbf w=(w_1,\dots,w_d,w_0)^{\top}\in\mathbb{R}^{d+1}, вектор откликов \mathbf y=(y_1,\dots,y_n)^{\top}. Тогда модель компактно записывается одной строкой:
\mathbf y\;=\;X\,\mathbf w+\boldsymbol\xi.
Метод наименьших квадратов. Сумма квадратов ошибок принимает вид
S(\mathbf w)\;=\;\sum_{i=1}^{n}\bigl(y_i-\mathbf x_i^{\top}\mathbf w\bigr)^{2} \;=\;\left\lVert \mathbf y-X\mathbf w \right\rVert^{2},
где для краткости \mathbf x_i^{\top}\mathbf w=x_{i,1}w_1+\dots+x_{i,d}w_d+w_0 (т. е. стандартное скалярное произведение). Приравнивание градиента S нулю даёт систему нормальных уравнений:
\tag{5.28} \boxed{\;\;(X^{\top}X)\,\hat{\mathbf w}\;=\;X^{\top}\mathbf y\;\;}\qquad(\text{нормальные уравнения}).
Если матрица X^{\top}X размера (d+1)\times(d+1) невырождена, отсюда \hat{\mathbf w}=(X^{\top}X)^{-1}X^{\top}\mathbf y. Эта формула — наиболее известная в линейной алгебре статистической оценки.
Условие применимости. Невырожденность X^{\top}X равносильна тому, что столбцы X линейно независимы: ни один признак не является линейной комбинацией других, и число наблюдений не меньше числа признаков (n\geq d+1). Если это не так — решение не единственно, и приходится либо отбирать признаки, либо использовать регуляризацию (см. 5.18).
Пусть для пяти квартир известны площадь x_{i,1} (м^{2}), расстояние до метро x_{i,2} (минуты пешком) и цена y_i (млн руб.):
\begin{array}{c|c|c|c} i & x_{i,1} & x_{i,2} & y_i \\ \hline 1 & 35 & 5 & 9{,}5 \\ 2 & 45 & 12 & 10{,}0 \\ 3 & 60 & 8 & 13{,}5 \\ 4 & 75 & 3 & 17{,}0 \\ 5 & 50 & 20 & 9{,}8 \end{array}
Подставим X\in\mathbb{R}^{5\times 3} (последний столбец из единиц) в формулу (5.28) и численно решим систему. Получим \hat{\mathbf w}\approx(0{,}183,\,-0{,}142,\,2{,}74), то есть
\hat y\;=\;0{,}183\,x_{1}\;-\;0{,}142\,x_{2}\;+\;2{,}74.
Содержательно: каждый дополнительный квадратный метр прибавляет к цене \approx\!183 тыс. руб., каждая дополнительная минута до метро уменьшает цену на \approx\!142 тыс. руб., а свободный член 2{,}74 млн — это «базовая стоимость», когда оба признака равны нулю.
Главная идея. С ростом размерности тип формул не меняется. МНК — универсальный инструмент, одинаково работающий и для одного-единственного признака (Кеплер, Галилей), и для тысячи признаков (модели ценообразования, прогнозирование спроса, физика частиц). Мы вернёмся к этим формулам в параграфе 5.4 при обсуждении малоранговых аппроксимаций матриц.
5.2.8 Резюме параграфа
Метод наименьших квадратов — первый реальный инструмент анализа данных, который мы получили. Его эквивалентность принципу максимума правдоподобия в гауссовой модели делает его универсальным.
На примере Кеплера МНК даёт III закон с точностью до семи знаков: наклон \hat\beta=1{,}49996\approx 3/2. На примере Галилея — ускорение свободного падения \hat g\approx 9{,}51 м/с^{2} при истинном 9{,}81.
Главные идеи, которые мы должны вынести из этого параграфа и которые будут проявляться дальше в более сложных моделях:
Принцип максимума правдоподобия, применённый к модели «параметр плюс гауссов шум», даёт минимизацию суммы квадратов как универсальный функционал ошибки.
Принцип Ферма (приравнивание частных производных к нулю) приводит к системе линейных уравнений, которая называется системой нормальных уравнений.
Зависимость, не являющаяся линейной, часто спрямляется заменой переменных (логарифмирование, возведение в квадрат) и после этого тоже становится задачей МНК.
В следующем параграфе мы шагнём от линейных моделей к существенно нелинейным — нейронным сетям. Принцип максимума правдоподобия и тут останется методологической отправной точкой, но функционал станет сложнее, а вместо аналитического решения нормальных уравнений мы будем использовать численные методы — градиентный спуск и его варианты.
Для точек (1,1{,}2), (2,2{,}3), (3,3{,}5), (4,4{,}3), (5,5{,}1) найдите \hat a,\hat b методом наименьших квадратов. Постройте график «данные + регрессионная прямая».
Выведите формулу (5.22) в явном виде, проделав все преобразования между (5.19) и (5.22).
Кеплер. По данным таблицы выше посчитайте вручную (или с помощью калькулятора) \bar x=\overline{\lg a}, \bar y=\overline{\lg T}, ковариацию и дисперсию, и убедитесь, что МНК даёт \hat\beta\approx 1{,}5.
Галилей. Зачем в (5.26) мы добавили свободный член \xi? Какое физическое явление он моделирует?
Постройте регрессию веса человека (y) от его роста (x) на выборке из 10 ваших знакомых (можно вымышленных). Сравните ваш результат с известной эмпирической формулой Брока y\approx x-100 (см. массой в кг, рост в см).
^{\star} Покажите, что МНК-оценка \hat a из (5.22) несмещённая: \mathbb E[\hat a]=a, если на самом деле y_i=a x_i+b+\xi_i с независимыми ошибками с нулевым математическим ожиданием. Найдите \mathbb D[\hat a].
^{\star} Метод наименьших квадратов с весами. Пусть нам известно, что точность i-го измерения характеризуется собственной дисперсией \sigma_i^{2} (разные приборы — разная погрешность). Покажите, что ОМП-оценка минимизирует взвешенную сумму квадратов \sum_{i}\sigma_i^{-2}(y_i-a x_i-b)^{2}. Найдите для этой задачи аналог формул (5.22), (5.21).
5.3 Задача классификации и нейронные сети
До сих пор мы решали задачи, в которых ответ модели — вещественное число. В задаче выборов это была доля голосов; у Гальтона — масса быка; у Кеплера — период обращения. Такие задачи называют задачами регрессии.
Но огромный класс практических задач устроен иначе. Алгоритму показывают фотографию — он должен сказать «кошка» или «собака». Письмо приходит в почтовый ящик — система решает: «спам» или «не спам». Изображение с томографа поступает в больницу — модель должна определить, есть ли на нём патология. Это задачи классификации: каждому входному объекту нужно сопоставить дискретную метку из конечного списка возможных классов.
В этом параграфе мы:
познакомимся с классической задачей классификации — распознаванием рукописных цифр;
покажем, как из принципа максимума правдоподобия выводится функционал кросс-энтропии — стандартная мера качества классификатора;
построим простейшие модели — полносвязную сеть и свёрточную сеть — и сравним их;
обнаружим важнейшее явление — переобучение (обнаружат его и сами вычислительные эксперименты);
обсудим теоретическое обоснование того, что нейронные сети могут приближать любую разумную зависимость — теоремы Колмогорова–Арнольда и Цыбенко;
изучим основной численный метод обучения — градиентный спуск — и его эффективную реализацию для нейронных сетей — обратное распространение ошибки.
5.3.1 Задача классификации MNIST и UCI Digits
Самый знаменитый «учебный» датасет для классификации — это MNIST (Modified National Institute of Standards and Technology database): 70 000 чёрно-белых изображений рукописных цифр размера 28\times 28 пикселей. Его собрал в 1998 г. Янн ЛеКун, который позже получил премию Тьюринга за работы по глубокому обучению. Все обсуждения и итоговые выводы этой главы относятся именно к MNIST как к каноническому ориентиру — именно его читатель встретит во всех профессиональных курсах и публикациях по машинному обучению.
Почему в численных экспериментах книги — не MNIST. Чтобы все приводимые ниже эксперименты можно было полностью воспроизвести на школьном ноутбуке за несколько секунд, мы используем UCI Optical Digits — меньший близкий по структуре датасет: те же 10 классов цифр (0–9), но изображения 8\times 8 (всего 1 797 примеров вместо 70 000). На полном MNIST обучение глубокой сети занимает минуты или десятки минут даже на современном железе и плохо ложится в формат «прочитал страницу учебника — запустил код — получил картинку». Качественные эффекты — переобучение, U-образная кривая ошибки, преимущество свёрточных сетей перед полносвязными — проявляются и на «миниатюре» точно так же; в подписях к рисункам сравнения с полноразмерным MNIST и ImageNet будут приведены отдельно (рис. 5.13).
Математическая постановка. Изображение размера 8\times 8 — это таблица из 64 чисел (значений яркости от 0 до 16). «Развернём» эту таблицу в один длинный вектор:
\mathbf x \;=\; (x_1,\,x_2,\,\dots,\,x_{64}) \;\in\; \mathbb{R}^{64}.
Так каждая цифра превращается в точку 64-мерного пространства. Для полноразмерного MNIST — 784-мерного (рис. 5.14).
Метка — одно из 10 значений: y\in\{0,1,\dots,9\}. Датасет — это набор пар
\{(\mathbf x_i,\,y_i)\}_{i=1}^{N},\qquad \mathbf x_i\in\mathbb{R}^{64},\quad y_i\in\{0,\dots,9\}.
Задача классификации — построить функцию-классификатор \hat y=\hat y(\mathbf x), которая по новому, ранее не виденному вектору \mathbf x предсказывает правильную метку.
5.3.2 Чем классификация отличается от регрессии
В задаче регрессии (5.2) функция, которую мы искали, отображала \mathbb{R}^{d}\to\mathbb{R} — то есть выдавала одно вещественное число. В задаче классификации мы хотим, чтобы алгоритм выдавал категориальный ответ. Принципиальные отличия — два.
Отличие 1: ответ не число, а распределение вероятностей. Хорошая модель должна выдавать не «эта цифра — семёрка», а более информативно: «семёрка с вероятностью 92 %, единица — 5 %, четвёрка — 2 %, остальное — меньше процента». Соответственно выход модели — это вектор из 10 чисел, неотрицательных и суммирующихся в 1:
\mathbf p(\mathbf x) \;=\; \bigl(p_0(\mathbf x),\,p_1(\mathbf x),\,\dots,\,p_9(\mathbf x)\bigr), \qquad p_c \geq 0,\quad \sum_{c=0}^{9} p_c=1.
Здесь p_c(\mathbf x) — предсказанная моделью вероятность того, что \mathbf x изображает цифру c. Окончательное предсказание метки — \hat y(\mathbf x)=\operatorname*{arg\,max}_{c} p_c(\mathbf x).
Отличие 2: функционал ошибки — не сумма квадратов, а кросс-энтропия. В регрессии мы минимизировали \sum_i (y_i-\hat y_i)^{2} — это была ОМП при гауссовом шуме. Сейчас «гауссова шума» нет: правильная метка y_i — это дискретная категория. Что вместо суммы квадратов? Принцип максимума правдоподобия — даст ответ.
Вывод кросс-энтропии из ОМП
Модель порождения данных. Считаем, что наблюдённые пары (\mathbf x_i,y_i) порождены так: для каждого i метка y_i — случайная величина со значениями в \{0,1,\dots,9\} и распределением, которое зависит от \mathbf x_i:
\Pr(y_i=c\,|\,\mathbf x_i,\boldsymbol\theta)\;=\;p_c(\mathbf x_i;\boldsymbol\theta).
Здесь \boldsymbol\theta — параметры модели (веса и смещения нейронной сети), которые предстоит определить. Совместное правдоподобие выборки размера N — это произведение вероятностей:
L(\boldsymbol\theta)\;=\;\prod_{i=1}^{N} p_{y_i}(\mathbf x_i;\boldsymbol\theta).
Берём отрицательный логарифм (минимизировать удобнее, чем максимизировать) и нормируем на размер выборки:
\tag{5.29} \boxed{\;\; \mathcal L(\boldsymbol\theta)\;\stackrel{\mathrm{def}}{=}\; -\frac{1}{N}\sum_{i=1}^{N}\ln p_{y_i}(\mathbf x_i;\boldsymbol\theta)\;\longrightarrow\;\min_{\boldsymbol\theta}. \;\;}
Этот функционал называется кросс-энтропией1 (англ. cross-entropy, categorical log-loss). Его минимизация по \boldsymbol\theta — это в точности максимизация правдоподобия. Мы снова применили один и тот же универсальный принцип — и получили совершенно другой функционал, естественный для своей задачи.
В теории информации Клод Шеннон (1948) ввёл понятие энтропии распределения \mathbf q:
H(\mathbf q)\;=\;-\sum_c q_c\ln q_c,
характеризующее «непредсказуемость» исхода. Кросс-энтропия двух распределений \mathbf q (истинное) и \mathbf p (предсказанное) определяется как
H(\mathbf q,\mathbf p)\;=\;-\sum_c q_c\ln p_c.
В нашей задаче \mathbf q — это «one-hot» вектор истины (q_c=1 только при c=y_i, иначе 0), а \mathbf p — предсказание модели. Тогда H(\mathbf q_i,\mathbf p_i)=-\ln p_{y_i}, и формула (5.29) оказывается средней по выборке кросс-энтропией.
5.3.3 Полносвязная нейронная сеть
Нам осталось определить семейство функций p_c(\mathbf x;\boldsymbol\theta), по которому будет проводиться минимизация. Самое известное — полносвязная нейронная сеть, она же многослойный персептрон (MLP, multilayer perceptron).
Идея «искусственного нейрона» — линейная функция от входов, пропущенная через нелинейность — восходит к работам Уоррена Маккаллока и Уолтера Питтса (1943) и обучаемому персептрону Фрэнка Розенблатта (1958). Розенблатт построил физическую машину Mark I Perceptron на лампах и фоторезисторах; пресса 1958 г. восторженно писала, что «машина скоро будет ходить, разговаривать и осознавать своё существование».
Однако в 1969 г. Марвин Минский и Сеймур Пейперт в книге Perceptrons строго доказали, что один-единственный слой не может реализовать даже простейшую функцию — исключающее ИЛИ (XOR). Это охладило интерес к нейросетям почти на 15 лет.
Возрождение пришло в 1986 г., когда Дэвид Румельхарт, Джеффри Хинтон и Рональд Уильямс представили универсальный алгоритм обучения многослойных сетей — метод обратного распространения ошибки (backpropagation). С этого момента и начинается эпоха современных нейронных сетей. Подчеркнём, что общая математическая формула того, что вычислительная сложность подсчёта градиента сравнима со сложностью самого вычисления функции, в общем случае была строго установлена советскими математиками: Юрий Геннадьевич Ким, Юрий Евгеньевич Нестеров, Александр Скоков, Виктор Черкасский (работы начала 1980-х). Об этом подробнее ниже, в разделе 5.22.
В 2018 г. премию Тьюринга — «нобелевскую» премию по информатике — получили Янн ЛеКун, Джеффри Хинтон и Йошуа Бенжио — «за концептуальные и инженерные прорывы, благодаря которым глубокие нейронные сети стали важнейшим компонентом современной информатики».
Строгое определение. Нам понадобится понятие матрично-векторного умножения.
Пусть A=(a_{ij}) — матрица размера m\times n (то есть m строк и n столбцов), и \mathbf z=(z_1,\dots,z_n)^{\top} — вектор из \mathbb{R}^{n}. Произведение A\mathbf z — это вектор из \mathbb{R}^{m}, i-я координата которого равна
\tag{5.30} (A\mathbf z)_i\;=\;\sum_{j=1}^{n} a_{ij}\,z_j,\qquad i=1,\dots,m.
Иначе говоря, мы скалярно умножаем i-ю строку матрицы на вектор и получаем i-ю координату результата.
A=\begin{pmatrix}1 & 2 & -1\\ 3 & 0 & 4\end{pmatrix},\qquad \mathbf z=\begin{pmatrix}1\\2\\3\end{pmatrix}.
Тогда
A\mathbf z=\begin{pmatrix}1\cdot 1+2\cdot 2+(-1)\cdot 3\\ 3\cdot 1+0\cdot 2+4\cdot 3\end{pmatrix}= \begin{pmatrix}2\\15\end{pmatrix}.
Теперь можем определить полносвязную сеть (рис. 5.15).
Полносвязная нейронная сеть (MLP) глубины L — это функция \mathbf x\mapsto \mathbf p(\mathbf x), заданная как композиция L слоёв:
$$ \begin{aligned} \mathbf z^{(\ell)} &= W^{(\ell)} \mathbf a^{(\ell-1)} + \mathbf b^{(\ell)},\qquad \mathbf a^{(\ell)} = \sigma\bigl(\mathbf z^{(\ell)}\bigr), \qquad \ell=1,\dots,L-1, \\ \mathbf z^{(L)} &= W^{(L)} \mathbf a^{(L-1)} + \mathbf b^{(L)},\qquad \mathbf p(\mathbf x) = \mathrm{softmax}\bigl(\mathbf z^{(L)}\bigr). \end{aligned}$$
Здесь \mathbf a^{(0)}=\mathbf x (вход), W^{(\ell)}\in\mathbb{R}^{d_\ell\times d_{\ell-1}} — матрица весов \ell-го слоя, \mathbf b^{(\ell)}\in\mathbb{R}^{d_\ell} — вектор смещений (bias). Функция \sigma применяется покомпонентно — это функция активации (см. ниже). Числа d_1,\dots,d_{L-1} — ширины скрытых слоёв; d_L — число классов (у нас — 10). Функция softmax превращает вектор произвольных вещественных чисел в распределение вероятностей:
\tag{5.33} \mathrm{softmax}(\mathbf z)_c\;=\;\frac{e^{z_c}}{\sum_{c'} e^{z_{c'}}}.
Совокупность всех матриц W^{(\ell)} и векторов \mathbf b^{(\ell)} объявляется обучаемыми параметрами \boldsymbol\theta модели.
Содержательная интерпретация. Без функций активации \sigma композиция линейных преобразований оставалась бы линейной (произведение матриц — снова матрица), и вся сеть свелась бы к простой линейной модели. Именно нелинейность делает сеть выразительной — позволяет ей описывать сложные, непрямолинейные зависимости.
Подсчёт числа параметров. Для сети 64 \to 32 \to 10 (то есть один скрытый слой шириной 32):
\#\boldsymbol\theta\;=\;\underbrace{64\cdot 32 + 32}_{\text{слой 1}}\;+\; \underbrace{32\cdot 10+10}_{\text{слой 2}} \;=\;2080+330\;=\;2410.
Сравните: датасет UCI Digits содержит около 1800 изображений, то есть параметров уже больше, чем обучающих примеров. Эта диспропорция — типичная для современных нейросетей — ставит важный вопрос о переобучении, к которому мы скоро вернёмся.
5.3.4 Обучение MLP и кривые потерь
Минимизация функционала (5.29) в нейронных сетях — это большая задача оптимизации в пространстве \boldsymbol\theta размерности от тысяч до миллиардов. Аналитического решения нет (как было в МНК), поэтому применяют численные методы — прежде всего градиентный спуск. Подробно к нему мы вернёмся в разделе 5.21; пока примем как данность и посмотрим на поведение функции потерь и точности в процессе обучения.
В этом эксперименте модель достигает \approx 92\,\% точности на тесте — вполне приличный результат для столь простой архитектуры.
5.3.5 Переобучение и U-образная кривая
Что произойдёт, если взять модель «помощнее»? Интуитивно: чем больше параметров, тем больше выразительная сила, тем лучше модель должна работать. Это верно до определённого предела — и катастрофически неверно после него. Это и называется переобучением.
Переобучение (англ. overfitting) — ситуация, когда модель отлично запоминает обучающую выборку (ошибка на train близка к нулю), но плохо работает на ранее не виденных данных (ошибка на test велика). Модель «выучила шум», а не закономерность.
Простейший пример — полиномиальная регрессия. Возьмём 12 точек, порождённых функцией y=\sin 2\pi x плюс небольшой гауссов шум, и будем приближать их полиномами растущей степени.
Три выборки: train, validation, test. Стандартный приём борьбы с переобучением — разбить имеющиеся данные на три непересекающиеся части.
Обучающая (train) — на ней непосредственно подбирают параметры \boldsymbol\theta (веса сети).
Валидационная (validation) — на ней подбирают гиперпараметры: глубину сети, ширину скрытого слоя, скорость обучения, степень полинома и т. п. Это не обучение в строгом смысле — мы лишь сравниваем разные варианты модели.
Контрольная (test) — её модель видит только один раз, в самом конце, и только для того, чтобы честно оценить качество финальной модели. Если test-выборка используется неоднократно для подбора, она перестаёт быть «честной» — происходит так называемая утечка информации.
Типичные пропорции: 60–80\,\% на обучение, по 10–20\,\% на остальные две.
5.3.6 Теоремы об универсальной аппроксимации
Возникает естественный вопрос: что в принципе может выучить нейронная сеть? Может ли она приблизить произвольную непрерывную функцию — или есть классы зависимостей, недоступные ей в принципе?
Ответ дают два глубоких математических результата, открытые с разрывом в почти 30 лет.
Всякая непрерывная функция f\colon[0,1]^{n}\to\mathbb{R} от n переменных представима в виде суперпозиции непрерывных функций одной переменной и операции сложения:
\tag{5.34} f(x_1,\dots,x_n)\;=\;\sum_{q=0}^{2n}\Phi_q\!\left(\sum_{p=1}^{n}\psi_{q,p}(x_p)\right),
где \Phi_q, \psi_{q,p} — некоторые непрерывные функции одной переменной.
Эта теорема была доказана А. Н. Колмогоровым в 1957 г. и усилена В. И. Арнольдом — в виде, эквивалентном решению 13-й проблемы Гильберта (доказательство того, что любая непрерывная функция трёх переменных представима суперпозицией функций двух переменных). С точки зрения нейронных сетей утверждение (5.34) говорит: двух слоёв и нескольких десятков аккуратно выбранных одномерных функций уже хватает для представления любой непрерывной функции от n переменных. Это, конечно, теоретическое утверждение: функции \Phi_q,\psi_{q,p} из доказательства Колмогорова чрезвычайно сложно устроены и непрактичны для прямого обучения. Но идейно теорема заложила фундамент.
Современная формулировка — с обычными функциями активации, такими как сигмоида или ReLU — была получена Дж. Цыбенко.
Пусть \sigma\colon\mathbb{R}\to\mathbb{R} — непрерывная, ограниченная, монотонно возрастающая функция (например, сигмоида). Тогда для любой непрерывной функции f на компакте K\subset\mathbb{R}^{n} и любого \varepsilon>0 существуют число N, числа a_j,b_j,w_{j,i}, j=1,\dots,N, i=1,\dots,n, такие что функция
\tag{5.35} g(\mathbf x)\;=\;\sum_{j=1}^{N} a_j\,\sigma\!\left(\sum_{i=1}^{n}w_{j,i}\,x_i+b_j\right)
приближает f равномерно на K с точностью \varepsilon: \sup_{\mathbf x\in K}|g(\mathbf x)-f(\mathbf x)|<\varepsilon.
Иначе говоря, одного скрытого слоя достаточной ширины достаточно, чтобы приблизить любую непрерывную функцию с заданной точностью. Это теорема об универсальной аппроксимации (рис. 5.18).
Теорема Цыбенко требует очень широких сетей: чтобы достичь точности \varepsilon, ширина N может расти экспоненциально с размерностью n или с 1/\varepsilon. На практике куда выгоднее глубокие сети — те, что используют несколько слоёв средней ширины. Это утверждение было математически уточнено в работах 2010-х годов: для некоторых классов функций глубокая сеть требует экспоненциально меньше параметров, чем эквивалентная по точности неглубокая.
5.3.7 Свёрточные нейронные сети
Полносвязная сеть «не знает» о том, что её вход — это изображение. Она одинаково обработала бы цифру и любую случайную перестановку её 64 пикселей. Это и теоретически, и практически — неоптимально: в изображении соседние пиксели связаны (вместе образуют чёрточки, дуги, кружочки — из которых и состоит цифра), а далёкие — почти независимы.
Идея использовать локальные «свёрточные» фильтры пришла из биологии. В 1962 г. нейрофизиологи Дэвид Хьюбел и Торстен Визель показали, что зрительная кора кошки устроена иерархически: на первом уровне — нейроны, реагирующие на простые признаки (края, отрезки прямых под определённым углом); на следующих — на их комбинации (углы, дуги, контуры); и так до целых форм. За это открытие они получили Нобелевскую премию по физиологии в 1981 г.
В 1980 г. японский исследователь Кунихико Фукусима реализовал эту идею в виде неокогнитрона — иерархической сети для распознавания символов. А в 1989 г. Янн ЛеКун — тогда ещё аспирант у Хинтона — добавил к неокогнитрону обучение методом обратного распространения и предложил универсальную архитектуру свёрточной нейронной сети (Convolutional Neural Network, CNN). Его сеть LeNet-5 (1998) распознавала рукописные цифры с точностью 99\,\% и применялась для чтения почтовых индексов и банковских чеков.
В 2012 г. работа Алекса Крижевского, Ильи Суцкевера и Джеффри Хинтона ImageNet Classification with Deep Convolutional Neural Networks (сеть AlexNet) продемонстрировала, что свёрточные сети, обученные на графических процессорах (GPU), могут с разгромным счётом превзойти все предыдущие подходы в задаче классификации изображений. С этого момента и началась «революция глубокого обучения».
Что такое свёртка. Рассмотрим изображение X размера H\times W и небольшой фильтр K размера k\times k (обычно k=3 или 5). Свёрткой изображения с фильтром называется новое изображение Y, в каждой точке которого стоит «скалярное произведение» фильтра и маленького окошка изображения:
\tag{5.36} Y_{r,c}\;=\;\sum_{u=1}^{k}\sum_{v=1}^{k}K_{u,v}\,X_{r+u-1,\,c+v-1}.
Фильтр «скользит» по всему изображению; в итоге получается новая карта признаков, чувствительная к тому, что «видит» фильтр (см. рис. 5.19 — четыре фильтра, обученные на UCI Digits, очевидно выделяют локальные «уголки» и «полосы»).
Свёрточная нейронная сеть — это нейронная сеть, в которой некоторые слои представляют собой не полносвязные матрично-векторные умножения, а свёртки вида (5.36) с обучаемыми фильтрами K. Обычно после свёртки применяется функция активации и операция пулинга (выборки максимума по окрестности — для уменьшения разрешения). После нескольких свёрточных слоёв полученная «карта признаков» разворачивается в вектор и передаётся в обычный полносвязный слой с softmax.
Сравнение MLP и CNN на UCI Digits. Чтобы оценить выигрыш от свёрточной структуры, обучим две сети с примерно одинаковым числом параметров: обычную MLP и CNN с одним свёрточным слоем. Результат показан на рис. 5.20.
5.3.8 Градиент и градиентный спуск
Перейдём к ключевому вопросу: как же подобрать параметры \boldsymbol\theta, минимизирующие функционал (5.29)? Аналитическое решение, как было в МНК (5.10), здесь невозможно: функция \mathcal L(\boldsymbol\theta) очень сложная, многомерная, нелинейная. Поэтому применяют итерационные численные методы, и важнейший из них — градиентный спуск.
Производная по направлению и градиент
Пусть F\colon\mathbb{R}^{d}\to\mathbb{R} — дифференцируемая функция (в нашем случае F=\mathcal L, а d — число параметров сети). В точке \boldsymbol\theta зададим какое-нибудь направление \mathbf v\in\mathbb{R}^{d} единичной длины. Производной функции F в точке \boldsymbol\theta по направлению \mathbf v называется
\partial_{\mathbf v} F(\boldsymbol\theta)\;=\; \lim_{t\to 0}\frac{F(\boldsymbol\theta+t\mathbf v)-F(\boldsymbol\theta)}{t}.
Чтобы вычислить её, рассмотрим функцию одной переменной \varphi(t)=F(\boldsymbol\theta+t\mathbf v) и применим правило дифференцирования сложной функции:
\tag{5.37} \partial_{\mathbf v} F\;=\;\varphi'(0)\;=\;\sum_{j=1}^{d}\frac{\partial F}{\partial\theta_j}\cdot v_j \;=\;\langle \nabla F,\,\mathbf v\rangle,
где \nabla F=\bigl(\partial F/\partial\theta_1,\dots,\partial F/\partial\theta_d\bigr) — градиент функции F в точке \boldsymbol\theta, а \langle\cdot,\cdot\rangle — скалярное произведение в \mathbb{R}^{d}.
Среди всех единичных векторов \mathbf v\in\mathbb{R}^{d} производная \partial_{\mathbf v} F(\boldsymbol\theta) принимает минимальное значение при
\tag{5.38} \mathbf v^{\star}\;=\;-\frac{\nabla F(\boldsymbol\theta)}{\left\lVert \nabla F(\boldsymbol\theta) \right\rVert}.
При этом \partial_{\mathbf v^{\star}} F=-\left\lVert \nabla F(\boldsymbol\theta) \right\rVert.
Доказательство. По неравенству Коши–Буняковского \langle \nabla F,\mathbf v\rangle\geq -\left\lVert \nabla F \right\rVert\cdot\left\lVert \mathbf v \right\rVert=-\left\lVert \nabla F \right\rVert, с равенством в точности тогда, когда \mathbf v сонаправлен с -\nabla F.\square
Этот простой факт — сердце всей теории численной оптимизации. Локально, для уменьшения функции F, надо двигаться в сторону антиградиента. Именно так устроен метод градиентного спуска:
\tag{5.39} \boxed{\; \boldsymbol\theta_{k+1}\;=\;\boldsymbol\theta_k\;-\;\alpha\,\nabla F(\boldsymbol\theta_k) \;}
Здесь \alpha>0 — шаг обучения (learning rate), гиперпараметр. Слишком маленький — сходимость медленная; слишком большой — итерации могут «перескакивать» через минимум и расходиться (рис. 5.21).
Стохастический градиентный спуск (SGD)
Когда обучающая выборка велика, вычислять полный градиент \nabla\mathcal L(\boldsymbol\theta)=\frac{1}{N}\sum_{i=1}^{N}\nabla\ell_i(\boldsymbol\theta) (где \ell_i=-\ln p_{y_i}(\mathbf x_i;\boldsymbol\theta) — ошибка на i-м примере) на каждом шаге — очень дорого. Поэтому используют стохастический градиентный спуск (SGD): на каждой итерации случайно выбирается батч (англ. batch, «пачка») — подмножество B_k\subset\{1,2,\dots,N\} из m примеров, — и шаг делается по приближённому градиенту, вычисленному только по этому батчу:
\tag{5.40} \boldsymbol\theta_{k+1}\;=\;\boldsymbol\theta_k\;-\;\alpha\,\widehat{\nabla\mathcal L}(\boldsymbol\theta_k), \qquad \widehat{\nabla\mathcal L}=\frac{1}{m}\sum_{i\in B_k}\nabla\ell_i.
Здесь B_k — случайное подмножество индексов («мини-батч» на шаге k), |B_k|=m, обычно m\in\{32,64,128,256\}. Такой приближённый градиент — несмещённая оценка истинного градиента (в духе ОМП!), но вычисляется в N/m раз быстрее.
Дополнительный бонус: «шум» от случайного выбора батча помогает итерациям выбираться из неглубоких локальных минимумов и седловых точек. Поэтому SGD — де-факто стандартный оптимизатор глубоких сетей.
5.3.9 Метод обратного распространения ошибки
Чтобы реализовать градиентный спуск, надо уметь вычислять градиент функции потерь \mathcal L по параметрам \boldsymbol\theta. В нейронной сети это нетривиально: параметры \boldsymbol\theta запрятаны в матрицы W^{(\ell)} нескольких слоёв, и каждая участвует в формуле потерь через цепочку преобразований
\mathbf x\;\to\;\mathbf z^{(1)}\;\to\;\mathbf a^{(1)}\;\to\;\mathbf z^{(2)}\;\to\;\dots\;\to\;\mathbf p\;\to\;\mathcal L.
Прямой подсчёт частных производных «в лоб» — через определение производной как предела — потребовал бы для каждого параметра отдельного перевычисления всей сети. Для сети с миллионом параметров это означало бы миллион полных прогонов вместо одного — полностью неприемлемо.
Обратное распространение ошибки (backpropagation) — алгоритм, который вычисляет градиент сразу по всем параметрам с вычислительной стоимостью, сравнимой со стоимостью одного прогона функции. Идея алгоритма — системное применение цепного правила дифференцирования.
Backpropagation — классический пример переоткрытия в науке. Идея многократно появлялась с 1960-х:
Сеппо Линнайнмаа (1970) — финский магистр, в диссертации описал общий алгоритм автоматического дифференцирования в обратном режиме.
Ким Юрий Геннадьевич, Нестеров Юрий Евгеньевич, Скоков Александр, Черкасский Виктор — работы начала 1980-х в СССР: в общей форме строго установлен фундаментальный результат о том, что вычислительная сложность подсчёта градиента дифференцируемой функции, заданной программой, ограничена сверху константой, умноженной на сложность вычисления самой функции. Этот результат — математическая основа того, что современные «миллиарды параметров» вообще возможно обучать.
Поль Вербос (1974, докторская диссертация Гарварда) и Дэвид Румельхарт, Джеффри Хинтон, Рональд Уильямс (1986) — именно последняя работа в журнале Nature популяризовала backpropagation в применении к обучению многослойных нейронных сетей.
Сегодня backpropagation реализован в каждой библиотеке глубокого обучения (PyTorch, TensorFlow, JAX) как automatic differentiation.
Идея алгоритма для MLP.
Рассмотрим сеть из определения 5.12. Ради прозрачности выкладок ограничимся одним обучающим примером (\mathbf x,y); функция потерь \ell=-\ln p_y.
Прямой проход. Сначала вычисляем все промежуточные значения \mathbf z^{(\ell)},\mathbf a^{(\ell)} ровно так, как они определены в формулах (5.31), (5.32). Это — один обычный прогон сети.
Обратный проход. Хотим вычислить градиенты потерь по всем параметрам. Для \ell-го слоя обозначим
\tag{5.41} \boldsymbol\delta^{(\ell)}\;\stackrel{\mathrm{def}}{=}\;\frac{\partial \ell}{\partial \mathbf z^{(\ell)}}\;\in\;\mathbb{R}^{d_\ell}.
Это вспомогательный вектор; его называют «ошибкой \ell-го слоя». Алгоритм состоит из трёх формул, получаемых строгим применением цепного правила:
(B1) Ошибка на последнем слое. Для пары «softmax + кросс-энтропия» прямой подсчёт даёт красивую формулу:
\tag{5.42} \boldsymbol\delta^{(L)}\;=\;\mathbf p\;-\;\mathbf e_y,
где \mathbf e_y — one-hot вектор истины (1 в позиции y, 0 в остальных). То есть ошибка — это просто разница между предсказанным распределением вероятностей и истинным.
(B2) Распространение ошибки назад по слоям. Связь между ошибками двух соседних слоёв даётся формулой
\tag{5.43} \boldsymbol\delta^{(\ell-1)}\;=\;\bigl(W^{(\ell)}\bigr)^{\top}\boldsymbol\delta^{(\ell)}\;\odot\;\sigma'\!\bigl(\mathbf z^{(\ell-1)}\bigr),
где \odot — покомпонентное произведение векторов, а \sigma' применяется покомпонентно. Откуда она берётся? Распишем цепное правило для одного перехода между слоями: \mathbf z^{(\ell)}=W^{(\ell)}\mathbf a^{(\ell-1)}+\mathbf b^{(\ell)} и \mathbf a^{(\ell-1)}=\sigma(\mathbf z^{(\ell-1)}). Тогда
\frac{\partial\ell}{\partial z^{(\ell-1)}_j} =\sum_{i}\frac{\partial\ell}{\partial z^{(\ell)}_i}\cdot \frac{\partial z^{(\ell)}_i}{\partial a^{(\ell-1)}_j}\cdot \frac{\partial a^{(\ell-1)}_j}{\partial z^{(\ell-1)}_j} =\sum_{i}\delta^{(\ell)}_i\cdot W^{(\ell)}_{i,j}\cdot \sigma'\!\bigl(z^{(\ell-1)}_j\bigr).
Собирая всё по j в один вектор, получаем именно (5.43). Это и есть «распространение ошибки назад»: ошибка более глубокого слоя \ell «протаскивается» через транспонированную матрицу весов в более ранний слой \ell-1.
(B3) Градиенты по параметрам. Зная \boldsymbol\delta^{(\ell)}, градиенты по параметрам \ell-го слоя выписываются сразу:
\begin{aligned} \frac{\partial \ell}{\partial W^{(\ell)}} &= \boldsymbol\delta^{(\ell)}\bigl(\mathbf a^{(\ell-1)}\bigr)^{\top}, \\ \frac{\partial \ell}{\partial \mathbf b^{(\ell)}} &= \boldsymbol\delta^{(\ell)}. \end{aligned}
(Первое — матрица того же размера, что W^{(\ell)}; второе — вектор.)
Стоимость алгоритма. Каждая из формул (5.43), (5.44) — это одно матрично-векторное умножение, по стоимости такое же, как и прямой проход для того же слоя. Итого: полный обратный проход стоит не дороже одного прямого — то есть мы получаем все градиенты по всем параметрам ценой одного дополнительного прохода. Это и есть тот фундаментальный результат, который в общей форме был установлен в работах Кима–Нестерова–Скокова–Черкасского начала 1980-х: для любой дифференцируемой функции, заданной программой,
\boxed{\;\text{стоимость подсчёта градиента}\;\leq\;C\cdot\text{(стоимость подсчёта функции)},\;}
с константой C\leq 5 в наихудшем случае. Без этого результата обучение современных сетей с миллиардами параметров было бы вычислительно невозможно.
В практике глубокого обучения один из любимых рецептов отладки backpropagation — численная проверка градиента:
\frac{\partial\ell}{\partial\theta_j}\;\approx\;\frac{\ell(\theta_j+h)-\ell(\theta_j-h)}{2h}, \qquad h\sim 10^{-5}.
Это очень дорого (по одному вычислению функции на каждый параметр), поэтому делается только на десяти-двадцати случайных параметрах маленьких сетей, и только при отладке. Совпадение с аналитическим градиентом до 6-го знака означает, что backprop реализован правильно.
5.3.10 Резюме параграфа
Мы прошли длинный путь от формулировки задачи классификации до полноценного «рецепта» обучения нейронной сети. Ключевые идеи:
Кросс-энтропия — естественный функционал для задач классификации, выводимый из принципа максимума правдоподобия (как МНК выводился из ОМП для гауссова шума).
Нейронная сеть — параметризованное семейство функций в виде композиции линейных и нелинейных слоёв. Теоремы Колмогорова–Арнольда и Цыбенко гарантируют, что это семейство достаточно богато, чтобы приблизить любую непрерывную функцию.
Переобучение — ключевая опасность: с ростом сложности модели train-ошибка убывает монотонно, а test-ошибка ведёт себя по U-образной кривой. Для борьбы с этим выделяют валидационную и контрольную выборки.
Свёрточные сети учитывают пространственную структуру изображения и работают лучше полносвязных там, где локальность существенна.
Градиентный спуск и его стохастический вариант — универсальный численный метод. Двигаться надо в сторону антиградиента; обоснование — неравенство Коши–Буняковского.
Backpropagation — вычислительно эффективная реализация градиента для нейронных сетей; работает за время, сопоставимое с одним прямым проходом сети.
В следующем параграфе мы перейдём к задачам обучения без учителя, в которых правильные ответы y_i для обучающих объектов не известны. Удивительным образом и там удастся применить часть тех же идей.
Покажите формулу (5.42): вычислите производную кросс-энтропии \ell=-\ln p_y=-\ln\frac{e^{z_y}}{\sum_c e^{z_c}} по z_c и убедитесь, что она равна p_c-[c=y].
Сколько обучаемых параметров в сети с архитектурой 784\to 256\to 128\to 10 (стандартная MLP для полного MNIST)? А если добавить ещё один скрытый слой шириной 64?
Покажите, что произведение двух матриц A и B можно определить так, чтобы для любого вектора \mathbf v выполнялось (AB)\mathbf v=A(B\mathbf v). Какое тогда условие на размеры A и B?
Реализуйте на Python градиентный спуск для функции F(x,y)=\tfrac12(x^{2}+10y^{2}). При каких значениях шага \alpha метод сходится? При каком из них — быстрее всего?
Запрограммируйте простейшую полносвязную сеть с одним скрытым слоем на UCI Digits и обучите её. Постройте кривые потерь и точности — как на рис. 5.16.
^{\star} Возьмите 12 точек (x_i,y_i) из любой гладкой функции (например, y=\sin x при x\in[0,2\pi]) с добавленным гауссовым шумом. Постройте полиномы степеней 1,2,\dots,11 и нарисуйте U-кривую, как на рис. 5.17.
^{\star} Докажите формулу (5.43) полностью, начиная с цепного правила \partial\ell/\partial z^{(\ell-1)}_j=\sum_k (\partial\ell/\partial z^{(\ell)}_k)\cdot(\partial z^{(\ell)}_k/\partial z^{(\ell-1)}_j).
5.4 Задачи обучения без учителя
Во всех задачах предыдущих параграфов каждому обучающему объекту \mathbf x_i сопутствовал правильный ответ y_i: голос за кандидата, период обращения планеты, рукописная цифра. Такая задача называется обучением с учителем (supervised learning): учителем выступает разметка, кем-то заранее проставленная.
На практике разметка часто отсутствует. У онлайн-кинотеатра — терабайты данных о просмотрах и оценках, но никто не объяснил алгоритму, что «Иван — любитель боевиков, а Маша — драм». У поисковой системы — миллиарды веб-страниц, но никто не пометил, какие из них «про физику», а какие «про поэзию». В таких ситуациях говорят об обучении без учителя (unsupervised learning). Цель — самостоятельно обнаружить структуру в данных, выделить «главные оси», на которых эти данные расположены, и научиться компактно представлять каждый объект.
Этот параграф устроен так. Мы разберём два центральных сюжета — восстановление матриц по неполным данным (matrix completion, она же ставшая знаменитой Netflix problem) и автокодировщики. Затем покажем, что оба сводятся к универсальной алгебраической конструкции — сингулярному разложению матриц (SVD). В завершении обсудим эмбеддинги и контрастивное обучение, которые лежат в основе современных систем распознавания лиц, поиска и больших языковых моделей.
5.4.1 Задача восстановления матрицы. Netflix problem
В октябре 2006 г. американская компания Netflix — тогда ещё сервис проката DVD по почте — объявила публичный конкурс с призом в 1 миллион долларов: первый, кто построит алгоритм, угадывающий оценки пользователей хотя бы на 10 % точнее, чем собственный алгоритм Netflix, получит всю сумму.
Компания выложила открытый датасет: 480\,189 пользователей, 17\,770 фильмов и \approx 100 млн оценок (по 5-балльной шкале). Каждая оценка сопровождалась датой; нужно было предсказывать неизвестные оценки.
Конкурс длился почти три года. Победила команда BellKor’s Pragmatic Chaos (объединение исследователей из AT&T, Yahoo, Big Chaos и Pragmatic Theory), и выиграла она 21 сентября 2009 г. с разницей в 20 минут над командой Ensemble. Сердцевиной победившего метода была именно малоранговая матричная факторизация, о которой мы сейчас расскажем. Этот сюжет принципиально изменил индустрию рекомендательных систем — сегодня описанные подходы лежат в основе «вам может понравиться» в Кинопоиске, Окко, YouTube, Spotify, Amazon и сотнях других сервисов.
Формальная постановка. Пусть имеется m пользователей и n фильмов. Полная информация об оценках — это матрица A\in\mathbb{R}^{m\times n}:
A_{ij}\;=\;\text{оценка, которую пользователь $i$ поставил фильму $j$.}
Однако в реальности большинство элементов этой матрицы неизвестны — средний пользователь Netflix посмотрел и оценил лишь несколько десятков фильмов из почти 18 тысяч. Обозначим
\Omega\;\subset\;\{1,\dots,m\}\times\{1,\dots,n\}
множество наблюдённых пар (i,j). Задача восстановления матрицы по неполным данным: по известным A_{ij}, (i,j)\in\Omega, восстановить значения A_{ij} для (i,j)\notin\Omega (рис. 5.22).
Идея: малоранговая факторизация. Без дополнительных предположений задача безнадёжна: восстановить неизвестное число можно как угодно. Но если предположить, что есть скрытые признаки — например, «насколько фильм комедийный», «насколько драматичный», «сколько в нём действия» — и описать каждого пользователя i парой чисел, отражающей его пристрастия, а каждый фильм j — такой же парой чисел, то оценка естественным образом получается как «скалярное произведение пристрастий и свойств».
Формализуем. Пусть r — небольшое число (например, r=2, 5 или 20). Введём
матрицу U\in\mathbb{R}^{m\times r}: i-я её строка \mathbf u_i — вектор пристрастий пользователя i;
матрицу V\in\mathbb{R}^{n\times r}: j-я её строка \mathbf v_j — вектор свойств фильма j.
Тогда модель оценок:
\tag{5.46} \boxed{\;A_{ij}\;\approx\;\langle \mathbf u_i,\,\mathbf v_j\rangle \;=\;\sum_{k=1}^{r} u_{ik}\,v_{jk},\;}\qquad\text{или в матричной форме}\qquad A\;\approx\;U\,V^{\top}.
Сколько параметров в модели? В исходной матрице A — m\cdot n чисел (например, 480\,189\cdot 17\,770\approx 8{,}5 млрд). В факторизованной модели — m\cdot r+n\cdot r=(m+n)r чисел (для r=20, m\approx 5\cdot 10^{5}, n\approx 2\cdot 10^{4} — около 10^{7}, то есть на три порядка меньше). Если истинная структура данных — действительно малоранговая, то этой существенно меньшей информации достаточно для восстановления.
Рангом матрицы A называется наименьшее r, при котором A можно представить в виде A=UV^{\top} с матрицами размера m\times r и n\times r соответственно. Эквивалентное определение: ранг r — это максимальное число линейно независимых столбцов (а также строк) матрицы. Чем ниже ранг, тем больше избыточности (повторяющихся закономерностей) в данных — и тем более «сжимаема» матрица.
Постановка как задачи оптимизации. Параметры \mathbf u_i,\,\mathbf v_j подбираются так, чтобы модель наилучшим образом соответствовала наблюдённым данным:
\tag{5.47} \sum_{(i,j)\in\Omega}\bigl(A_{ij}-\langle\mathbf u_i,\,\mathbf v_j\rangle\bigr)^{2}\;\longrightarrow\;\min_{U,V}.
Это снова метод наименьших квадратов — но теперь сумма квадратов только по известным элементам матрицы. Если интерпретировать данные как «истинные оценки плюс гауссов шум», функционал (5.47) получается из принципа максимума правдоподобия — ровно тем же выводом, что в 5.8.
Этот пример хорош тем, что наглядно показывает водораздел. В задачах математической статистики (как в 5.1, 5.2) распределение данных с точностью до неизвестных параметров известно, и функционал ошибки выводится из принципа максимума правдоподобия.
В задачах машинного обучения и ИИ распределение чаще всего неизвестно; функционал ошибки выбирают исходя из здравого смысла. В нашем примере — квадратичная невязка (A_{ij}-\langle\mathbf u_i,\mathbf v_j\rangle)^{2} — это естественный выбор, но он не выводится из строгой модели порождения данных. Можно, конечно, постулировать, что A_{ij}=\langle\mathbf u_i,\mathbf v_j\rangle+\xi_{ij} с \xi_{ij}\sim\mathcal N(0,\sigma^{2}), и тогда квадраты возникают из ОМП — но это уже не вывод из данных, а наше дополнительное предположение о структуре шума. Кроме того, в ИИ часто используются нелинейные модели — нейронные сети как универсальные аппроксиматоры.
Неотрицательность. В практической рекомендательной системе оценки лежат в шкале от 1 до 5, а пристрастия и свойства принципиально «качественные»: фильм скорее в большей или меньшей степени принадлежит жанру, чем «в антинаправленной». Поэтому естественно ограничиться неотрицательной малоранговой факторизацией:
\tag{5.48} A\;\approx\;UV^{\top},\qquad U\geq 0,\;V\geq 0,
где неравенство понимается покомпонентно. Такая факторизация называется неотрицательным матричным разложением (non-negative matrix factorization, NMF); её свойства глубоко изучали Д. Ли и Х. Сын в работах 1999 г.
Рассмотрим частично заполненную матрицу 4 пользователей и 5 фильмов (оценки от 1 до 5, ‘?’ — неизвестно):
A=\begin{pmatrix} 5 & 4 & ? & 1 & ? \\ ? & 5 & 5 & ? & 2 \\ 1 & ? & 2 & 5 & 5 \\ ? & 1 & 1 & 4 & 5 \\ \end{pmatrix}.
Видна закономерность: первые два пользователя любят первые три фильма; третий и четвёртый — последние два. Если попытаться найти факторизацию A\approx UV^{\top} с рангом r=2:
U=\begin{pmatrix}1{,}9 & 0{,}1\\1{,}9 & 0{,}2\\0{,}1 & 1{,}7\\0{,}0 & 1{,}7\end{pmatrix}, \qquad V=\begin{pmatrix}2{,}6 & 0{,}5\\2{,}3 & 0{,}3\\2{,}5 & 0{,}5\\0{,}6 & 2{,}8\\0{,}4 & 3{,}0\end{pmatrix},
то UV^{\top} почти точно воспроизводит наблюдённые элементы и «предсказывает» неизвестные:
UV^{\top}\approx\begin{pmatrix}5{,}0 & 4{,}4 & \boxed{4{,}8} & 1{,}4 & \boxed{1{,}1}\\ \boxed{5{,}1} & 4{,}4 & 4{,}9 & \boxed{1{,}7} & 1{,}4\\ \boxed{1{,}1} & \boxed{0{,}7} & 1{,}1 & 4{,}8 & 5{,}1\\ 0{,}0 & 0{,}5 & 0{,}9 & 4{,}8 & 5{,}1\end{pmatrix}.
Обведены восстановленные элементы. Видно: пользователи 1 и 2 в этой модели — любители «жанра 1» (первая координата пристрастий велика); пользователи 3 и 4 — любители «жанра 2». А столбцы V показывают, какому жанру какие фильмы соответствуют.
Алгоритм решения. Задача (5.47) не имеет аналитического решения и относится к нелинейной невыпуклой оптимизации. На практике её решают методом переменных наименьших квадратов (Alternating Least Squares, ALS): фиксируем V и оптимизируем U (это уже выпуклая задача и решается по формулам линейной регрессии), потом фиксируем U и оптимизируем V. Чередуем до сходимости. Для очень больших матриц используется стохастический градиентный спуск (5.21).
5.4.2 Три кита и роль малоранговой аппроксимации в ИИ
Разобранный пример хорошо демонстрирует обещание начала параграфа: современный анализ данных стоит на трёх китах.
Вероятностный кит. Хотя ни одного распределения мы явно не вводили, постановка (5.47) имеет естественную вероятностную интерпретацию — «истинные значения плюс гауссов шум». Эта вероятностная подкладка — то, что превращает расплывчатое «найти хорошее восстановление» в строгую задачу.
Оптимизационный кит. Полученная задача оптимизации (5.47) — невыпуклая нелинейная, она требует подходящих численных методов и (для матриц размером в миллиарды элементов) распределённых вычислений на сотнях GPU. Современная индустрия рекомендательных систем — это в первую очередь огромная оптимизационная машинерия.
Линейно-алгебраический кит. Сама постановка — «представить матрицу в виде произведения двух маленьких» — это вопрос линейной алгебры. И тут возникает важная общая идея, далеко выходящая за рамки рекомендательных систем.
Идея малоранговой аппроксимации играет ключевую роль в современных больших языковых моделях. Когда нужно дообучить (fine-tune) готовую модель с миллиардами параметров под узкую задачу (юриспруденция, медицина, диалект), обновлять все веса — слишком дорого. Метод LoRA (Low-Rank Adaptation, Хью и др., 2021) предлагает: оставить исходные веса W_0 замороженными, а поверх каждой матрицы добавить малоранговую поправку \Delta W=A B^{\top}, где A,B — матрицы с очень небольшим рангом r (обычно r=8 или 16).
Обновляются только параметры A,B — их в сотни раз меньше, чем параметров полной модели. При этом результат дообучения почти неотличим от полного. LoRA-адаптеры (один на каждую узкую задачу) весят считанные мегабайты при том, что базовая модель — сотни гигабайтов; их можно «навинчивать» и снимать на лету. Это и доказывает фундаментальную важность идеи малоранговой аппроксимации в современном машинном обучении.
5.4.3 Автокодировщики: общая идея
Второй классический сюжет обучения без учителя — автокодировщик (autoencoder). Его базовая идея проста и красива: научить нейронную сеть воспроизводить свой собственный вход, но сделать это, проходя через узкое «горло» из небольшого числа нейронов. Тогда сеть будет вынуждена сжать всю существенную информацию о входе в этом горле — получится содержательное представление, или эмбеддинг (рис. 5.23).
Автокодировщик — это пара функций f\colon\mathbb{R}^{n}\to\mathbb{R}^{r} (энкодер) и g\colon\mathbb{R}^{r}\to\mathbb{R}^{n} (декодер), r<n, обученные на выборке \{\mathbf x_1,\dots,\mathbf x_N\}\subset\mathbb{R}^{n} минимизацией функционала
\tag{5.49} \mathcal L(f,g)\;=\;\frac{1}{N}\sum_{i=1}^{N}\left\lVert \mathbf x_i-g(f(\mathbf x_i)) \right\rVert^{2}.
Вектор \mathbf z_i=f(\mathbf x_i)\in\mathbb{R}^{r} называется эмбеддингом (embedding, скрытым представлением) объекта \mathbf x_i. Размерность r — гиперпараметр модели.
5.4.4 Линейный автокодировщик и связь с PCA
Самый простой случай — когда f и g линейны: f(\mathbf x)=W_e\mathbf x, g(\mathbf z)=W_d\mathbf z, где W_e\in\mathbb{R}^{r\times n} и W_d\in\mathbb{R}^{n\times r} — обучаемые матрицы (для простоты считаем, что данные предварительно центрированы, \overline{\mathbf x}=\mathbf 0).
Постановка. Функционал (5.49) принимает вид
\tag{5.50} \mathcal L(W_d,W_e)\;=\;\frac{1}{N}\sum_{i=1}^{N}\left\lVert \mathbf x_i-W_d W_e \mathbf x_i \right\rVert^{2}\;\longrightarrow\;\min_{W_d,W_e}.
Решение и его геометрия. Обозначим X\in\mathbb{R}^{n\times N} — матрицу, столбцы которой — центрированные обучающие вектора \mathbf x_1,\dots,\mathbf x_N. Введём выборочную ковариационную матрицу C=\tfrac{1}{N}X X^{\top}\in\mathbb{R}^{n\times n}. Эта матрица симметрична и неотрицательно определена.
Матрица C\in\mathbb{R}^{n\times n} называется симметричной, если C=C^{\top}, и неотрицательно определённой, если \langle \mathbf v,C\mathbf v\rangle\geq 0 для любого \mathbf v\in\mathbb{R}^{n}. Спектральная теорема линейной алгебры утверждает: для любой симметричной матрицы C существует ортонормированный базис \mathbf e_1,\dots,\mathbf e_n из её собственных векторов, C\mathbf e_k=\lambda_k\mathbf e_k. Для неотрицательно определённой матрицы все \lambda_k\geq 0. Геометрически: оператор с матрицей C в этом базисе — просто покоординатное растяжение в \lambda_k раз вдоль k-й оси.
Упорядочим собственные значения C по убыванию: \lambda_1\geq\lambda_2\geq\dots\geq\lambda_n\geq 0, с соответствующими ортонормированными собственными векторами \mathbf e_1,\dots,\mathbf e_n.
Минимум функционала (5.50) (при условии, что f и g линейные, r<n) достигается, когда подпространство, порождённое строками W_e, совпадает с подпространством, натянутым на первые r собственных векторов \mathbf e_1,\dots,\mathbf e_r ковариационной матрицы C. При этом действие энкодера f есть ортогональная проекция на это подпространство.
(Доказательство опирается на спектральное разложение и теорему Эккарта–Янга–Мирского, к которой мы перейдём в 5.29.)
Что это означает по содержанию. Линейный автокодировщик ровно тот же объект, что метод главных компонент (Principal Component Analysis, PCA), открытый К. Пирсоном в 1901 г. задолго до всех нейронных сетей.
Геометрически: энкодер — ортогональная проекция в подпространство максимальной дисперсии; декодер — обратное вложение из этого подпространства назад. Эмбеддинг точки — её координаты в этом подпространстве.
5.4.5 Нелинейные автокодировщики и их разновидности
Принципиальный шаг вперёд — разрешить f и g быть произвольными нейронными сетями. Тогда эмбеддинг \mathbf z может «свернуть» в свои несколько координат нелинейные комбинации исходных признаков, что несоизмеримо мощнее линейной проекции.
Свёрточные автокодировщики. Если вход — изображение, естественно использовать свёрточные слои в энкодере (5.20) и обратные операции (deconvolution, или transposed convolution) в декодере. Так строятся свёрточные автокодировщики, лежащие в основе ранних систем сжатия изображений и шумоподавления.
Разреженный автокодировщик. Добавим к функционалу (5.49) штрафное слагаемое, требующее, чтобы эмбеддинг был разрежён (т. е. большинство его координат близки к нулю):
\mathcal L_{\text{sparse}}\;=\;\mathcal L+\lambda\sum_{i,k}|z_{ik}|, \qquad z_{ik}=f(\mathbf x_i)_k.
Это заставляет сеть «активировать» лишь несколько координат эмбеддинга для каждого входа — получается интерпретируемое представление, в котором за каждым активным нейроном «закрепляется» свой признак.
Шумоподавляющий автокодировщик. Берём \mathbf x_i, портим его шумом \tilde{\mathbf x}_i=\mathbf x_i+\boldsymbol\eta, а сеть учим восстанавливать чистое \mathbf x_i по шумному \tilde{\mathbf x}_i:
\mathcal L_{\text{denoise}}\;=\;\frac{1}{N}\sum_{i=1}^{N}\left\lVert \mathbf x_i-g(f(\tilde{\mathbf x}_i)) \right\rVert^{2}.
Это вариант, предложенный П. Венсаном и Й. Бенжио в 2008 г., — очень полезен на этапе предобучения глубоких сетей. Эта идея является концептуальным прародителем современных диффузионных моделей (Stable Diffusion, DALL-E — 2022 г. и позднее): они обучаются шаг за шагом удалять гауссов шум из изображения и таким образом порождают совершенно новые картинки из «случайной соли и перца».
Вариационный автокодировщик (VAE). Самая глубокая идейная модификация автокодировщика — вариационный автокодировщик (Кингма, Веллинг, 2013). Энкодер выдаёт не сам эмбеддинг, а параметры распределения над эмбеддингами:
f(\mathbf x)\;=\;\bigl(\boldsymbol\mu(\mathbf x),\,\boldsymbol\sigma(\mathbf x)\bigr),
после чего сам вектор \mathbf z выбирается случайно из нормального распределения \mathcal N\!\bigl(\boldsymbol\mu(\mathbf x),\boldsymbol\sigma(\mathbf x)^{2}\bigr) и подаётся в декодер. Обучаемая функция потерь складывается из двух слагаемых:
\mathcal L\;=\;\underbrace{\left\lVert \mathbf x-g(\mathbf z) \right\rVert^{2}}_{\text{невязка восстановления}} \;+\;\beta\,\underbrace{\mathrm{KL}\!\bigl(\mathcal N(\boldsymbol\mu,\boldsymbol\sigma^{2})\,\|\,\mathcal N(0,1)\bigr)}_{\text{регуляризатор}}.
Здесь \mathrm{KL} — расхождение Кульбака–Лейблера, или «расстояние» между двумя вероятностными распределениями (определяется в курсе теории вероятностей). В нашем случае оно штрафует отклонение распределения эмбеддингов \mathcal N(\boldsymbol\mu,\boldsymbol\sigma^{2}) от «эталонного» стандартного нормального \mathcal N(0,1). Содержательно это означает: похожие на вид изображения должны давать близкие \mathbf z, а всё пространство \mathbf z должно быть «равномерно заполненным».
Главное преимущество VAE — возможность генерации новых объектов: достаточно выбрать \mathbf z случайно из \mathcal N(0,1) и прогнать через декодер; благодаря регуляризатору почти любой такой \mathbf z декодируется в осмысленное изображение.
5.4.6 Эмбеддинги в современных моделях
Идея эмбеддинга вышла далеко за рамки автокодировщиков. Сегодня любая глубокая модель ИИ построена вокруг представления своих входов в виде векторов в каком-то многомерном пространстве — их называют общим словом эмбеддинги.
Эмбеддинги в компьютерном зрении. Рассмотрим обученную свёрточную сеть для классификации изображений (например, ResNet-50 на ImageNet с 1000 классов). Архитектурно её последние два слоя устроены так:
\underbrace{\;\mathbf x\;\longrightarrow\;\dots\;\longrightarrow\;\mathbf z\in\mathbb{R}^{2048}\;}_{\text{свёрточный «костяк»}} \;\xrightarrow{\;W\in\mathbb{R}^{1000\times 2048}\;}\; \mathbf p=\mathrm{softmax}(W\mathbf z)\in\mathbb{R}^{1000}.
Вектор \mathbf z на предпоследнем слое — это и есть эмбеддинг изображения: сжатое представление в \mathbb{R}^{2048}, из которого последнее линейное преобразование W восстанавливает вероятности классов. Эмпирически оказалось, что \mathbf z переносим: его можно использовать для совсем других задач — поиска похожих фотографий (сравнение по косинусной мере), классификации классов, которых не было в обучении (transfer learning), распознавания одного и того же лица на разных снимках. В этом и состоит главная практическая ценность эмбеддингов: одна обученная сеть даёт универсальное «представление изображения», пригодное для целого спектра приложений.
Эмбеддинги в больших языковых моделях. Каждое слово, а точнее — каждый токен (подпоследовательность символов), современная LLM кодирует вектором фиксированной длины (обычно 1024–4096 чисел). Эти эмбеддинги — обучаемые параметры модели; в процессе обучения семантически близкие слова приобретают близкие эмбеддинги (рис. 5.25).
Современная большая языковая модель (LLM; примеры — GPT-серии от OpenAI, Claude от Anthropic, Llama от Meta, GigaChat от Сбера, YandexGPT от Яндекса) построена на архитектуре трансформер (Васвани и др., 2017). На входе модель получает последовательность токенов; каждый токен превращается в эмбеддинг.
Дальше эти эмбеддинги многократно прогоняются через слои двух типов:
механизм внимания (self-attention): каждый токен «смотрит» на другие токены последовательности и обновляет свой эмбеддинг, учитывая контекст. Это и есть то, что делает модель «контекстно-зависимой»: то же слово в разных контекстах получает разный финальный эмбеддинг.
полносвязный блок (FFN, feedforward network) — обычный полносвязный слой большой ширины. Именно в этих блоках «хранятся знания» модели — факты, выученные из обучающих текстов.
В конце последний эмбеддинг конвертируется в распределение вероятностей над всеми токенами словаря — и модель предсказывает следующий токен. Так, генерируя по одному токену за раз, LLM пишет тексты, отвечает на вопросы, сочиняет код — всё это, по сути, последовательное предсказание следующего токена.
5.4.7 Сиамские сети и контрастивное обучение
В автокодировщике мы заставляли сеть восстанавливать вход целиком. Часто это избыточно: для практических задач нам нужно не «восстановить картинку до пикселя», а узнать, похожи ли две картинки между собой. Этим занимается контрастивное (или контрастное) обучение (contrastive learning).
Сиамская сеть. Сиамская сеть (Siamese network, Бромли с соавторами, 1993) состоит из двух копий одной и той же нейронной сети f с общими весами. На вход подают пару объектов; каждый прогоняется через свою копию f, и в результате получаются два эмбеддинга f(\mathbf x) и f(\mathbf y). Их близость или удалённость характеризует «похожесть» объектов (рис. 5.26).
Контрастивный функционал. Пусть имеется набор обучающих пар \{(\mathbf x_i,\mathbf y_i,t_i)\}, где t_i=1, если пара — «положительная» (похожие), t_i=0, если «отрицательная» (непохожие). Простейший контрастивный функционал (Хадселл, Шопра, ЛеКун, 2006):
\tag{5.51} \mathcal L\;=\;\sum_{i}\Bigl[\,t_i\cdot d_i^{2}\;+\;(1-t_i)\cdot\bigl(\max(0,\,m-d_i)\bigr)^{2}\,\Bigr],
где d_i=\left\lVert f(\mathbf x_i)-f(\mathbf y_i) \right\rVert, а m>0 — отступ (margin), гиперпараметр. Смысл: для похожих пар (t_i=1) минимизируется d_i^{2} — расстояние «прижимается» к нулю; для непохожих (t_i=0) штрафуется только то, что d_i меньше m — то есть от непохожих пар требуется минимальное расстояние m, а дальше уже всё равно.
Применения.
Распознавание лиц. Сиамская сеть FaceNet (Google, 2015) обучается на миллионах фотографий: эмбеддинг каждого лица — вектор в \mathbb{R}^{128}. Положительные пары — два снимка одного человека (часто в разных условиях); отрицательные — разных. После обучения по тестовому снимку можно определить личность простым поиском ближайшего эмбеддинга в базе. Точность — более 99{,}6\,\% на стандартном бенчмарке LFW. Такая система — основа всех современных систем распознавания лиц.
Поиск изображений по тексту (и обратно). Модель CLIP (OpenAI, 2021) обучает два энкодера: один — для изображений (свёрточная или трансформерная сеть), другой — для текстов (трансформер). Положительные пары — картинка и её настоящая подпись; отрицательные — картинка и случайная подпись из батча. После обучения эмбеддинги изображения и эмбеддинги его осмысленного текстового описания оказываются близкими в общем многомерном пространстве. Это даёт колоссальные возможности: искать картинки по произвольным описаниям («рыжий кот на пианино»), классифицировать изображения в произвольные классы, не имея размеченной выборки (zero-shot classification), и многое другое.
Дубликаты в поиске. Тот же контрастивный подход используют поисковые системы (Яндекс, Google) для дедупликации страниц, рекомендательные системы — для поиска похожих товаров, антивирусы — для группировки похожих образцов вредоносного кода.
5.4.8 Сингулярное разложение матриц^{\star}
Этот раздел — со звёздочкой: материал более продвинутый, и при первом чтении его можно пропустить. Однако именно здесь раскрывается глубокая связь между двумя сюжетами параграфа — матричной факторизацией и линейным автокодировщиком — через одну и ту же алгебраическую конструкцию.
Матрицы как линейные преобразования. Каждая матрица A\in\mathbb{R}^{m\times n} задаёт линейное преобразование \mathbb{R}^{n}\to\mathbb{R}^{m}, \mathbf x\mapsto A\mathbf x (см. определение 5.10). Столбцы матрицы A — это образы базисных ортов \mathbf e_1,\dots,\mathbf e_n: j-й столбец равен A\mathbf e_j. Произведение матриц C=AB соответствует суперпозиции преобразований: сначала действует B, потом A.
Простейшие линейные преобразования. Среди всех линейных преобразований \mathbb{R}^{n}\to\mathbb{R}^{n} выделяются два простейших класса.
Повороты — матрицы Q, сохраняющие длины: \left\lVert Q\mathbf x \right\rVert=\left\lVert \mathbf x \right\rVert. Они называются ортогональными и характеризуются равенством Q^{\top}Q=I.
Растяжения по осям — диагональные матрицы \Sigma: \Sigma\mathbf e_k=\sigma_k\mathbf e_k, \sigma_k\geq 0. Они растягивают k-ю ось в \sigma_k раз.
Сингулярное разложение. Знаменитый факт линейной алгебры (Бельтрами, 1873; Йордан, 1874; Сильвестр, 1889 для квадратных матриц; Экарт и Янг, 1936 для прямоугольных) — любое линейное преобразование можно представить как комбинацию поворота, растяжения и снова поворота:
Для любой матрицы A\in\mathbb{R}^{m\times n} существует представление
\tag{5.52} \boxed{\;A\;=\;U\,\Sigma\,V^{\top},\;}
где U\in\mathbb{R}^{m\times m} и V\in\mathbb{R}^{n\times n} — ортогональные матрицы (U^{\top}U=I_m, V^{\top}V=I_n), а \Sigma\in\mathbb{R}^{m\times n} — «диагональная» (т. е. \Sigma_{ij}=0 при i\neq j) с неотрицательными элементами на диагонали \sigma_1\geq\sigma_2\geq\dots\geq\sigma_{\min(m,n)}\geq 0. Эти числа называются сингулярными значениями матрицы A; столбцы U и V — сингулярными векторами.
Геометрически: преобразование с матрицей A можно представить последовательностью — сначала поворот в исходном пространстве (матрица V^{\top}), затем покоординатное растяжение (матрица \Sigma), затем поворот в пространстве-образе (матрица U) (рис. 5.27).
Полезный способ читать SVD: столбцы \mathbf v_1,\dots,\mathbf v_n матрицы V — это «входные направления», на которых матрица A действует чище всего (как простое растяжение, без поворота). Столбцы \mathbf u_1,\dots,\mathbf u_m матрицы U — соответствующие «выходные направления»: A\mathbf v_i=\sigma_i\,\mathbf u_i. Числа \sigma_i показывают, во сколько раз растягивается каждое такое направление; по соглашению они упорядочены по убыванию.
Связь с симметричной задачей. Доказательство SVD сводится к спектральной теореме для симметричной матрицы. Заметим:
A^{\top}A=(U\Sigma V^{\top})^{\top}(U\Sigma V^{\top})=V\Sigma^{\top}U^{\top}U\Sigma V^{\top}=V\Sigma^{\top}\Sigma V^{\top} — это спектральное разложение симметричной неотрицательной матрицы A^{\top}A. Её собственные значения — \sigma_k^{2}; собственные векторы — столбцы V.
Аналогично, AA^{\top}=U\Sigma\Sigma^{\top}U^{\top}. Собственные векторы AA^{\top} — столбцы U.
Так одно общее утверждение для произвольных матриц редуцируется к симметричному случаю.
Теорема о наилучшей малоранговой аппроксимации. SVD замечательно тем, что позволяет получить наилучшее приближение матрицы матрицей фиксированного ранга.
Пусть A\in\mathbb{R}^{m\times n} и A=U\Sigma V^{\top} — её сингулярное разложение. Для любого r\leq \min(m,n) определим обрезанное разложение
A_r\;=\;\sum_{k=1}^{r}\sigma_k\,\mathbf u_k\,\mathbf v_k^{\top}\;=\;U_r\,\Sigma_r\,V_r^{\top},
где U_r,V_r — первые r столбцов матриц U,V, а \Sigma_r=\mathrm{diag}(\sigma_1,\dots,\sigma_r). Тогда матрица A_r имеет ранг r и реализует наилучшую среди всех матриц ранга \leq r аппроксимацию исходной матрицы:
\tag{5.53} \min_{\mathrm{rank}\,B\leq r}\left\lVert A-B \right\rVert_{F}\;=\;\left\lVert A-A_r \right\rVert_{F}\;=\;\sqrt{\sigma_{r+1}^{2}+\dots+\sigma_{\min(m,n)}^{2}},
где \left\lVert \cdot \right\rVert_{F} — норма Фробениуса: \left\lVert M \right\rVert_{F}^{2}=\sum_{i,j}M_{ij}^{2}.
Иными словами, отбрасывая в разложении (5.52) все сингулярные значения, начиная с (r+1)-го, мы получаем оптимальную малоранговую аппроксимацию. Если первые r сингулярных значений велики, а последующие малы, то A_r — почти такая же хорошая, как A, но в \frac{mn}{(m+n)r} раз более компактна по числу параметров (рис. 5.28).
Связь с автокодировщиком и матричной факторизацией.
Линейный автокодировщик из 5.25 — это в точности SVD-факторизация центрированных данных. Если столбцы X — обучающая выборка, и X=U\Sigma V^{\top} её SVD, то оптимальный энкодер — f(\mathbf x)=U_r^{\top}\mathbf x, оптимальный декодер — g(\mathbf z)=U_r\mathbf z, где U_r — первые r столбцов матрицы U. Сингулярные значения \sigma_k показывают, насколько важна k-я главная компонента (большое \sigma_k — много дисперсии в этом направлении).
Матричная факторизация A\approx UV^{\top} из 5.23 — это (без ограничения на неотрицательность) ровно A\approx A_r из теоремы Эккарта–Янга–Мирского. Если бы все элементы A были известны, идеальной стратегией было бы взять обрезанное SVD-разложение ранга r. В реальной Netflix problem известны лишь \sim 1\,\% элементов матрицы, и приходится решать задачу (5.47), которая является обобщением SVD на случай пропусков.
Школа Е. Е. Тыртышникова и тензорные разложения. Идея малоранговой аппроксимации не ограничивается матрицами. Тензором называется многомерный массив чисел: матрица — двумерный тензор, а матрицы с тремя и более индексами — тензоры высших порядков. Тензоры возникают повсюду в современных приложениях: изображение видеоряда — 4-индексный тензор (время \times высота \times ширина \times цвет); веса трансформера — тысячи матриц, организованных в одну сложную структуру; квантово-механические задачи — тензоры с десятками индексов.
Аналог малоранговой аппроксимации для тензоров — значительно более тонкая задача. Современный мощный инструмент — разложение в виде тензорного поезда (tensor train decomposition), предложенное в 2009 г. российским математиком, профессором РАН Иваном Валерьевичем Оселедцем, учеником академика Евгения Евгеньевича Тыртышникова (МГУ, ИВМ РАН). Тензорный поезд позволяет хранить тензор с d индексами при размере N по каждому индексу с помощью \sim d\,N\,r^{2} параметров вместо N^{d} — то есть избавляет от знаменитого «проклятия размерности». Этот результат лежит в основе многих современных алгоритмов в квантовой химии, машинном обучении и решении многомерных дифференциальных уравнений.
За эти и связанные работы Е. Е. Тыртышников — один из мировых лидеров в области тензорных и матричных разложений — является лауреатом премий Сбербанка в области искусственного интеллекта.
5.4.9 Резюме параграфа
Обучение без учителя — это семейство задач, в которых модель должна самостоятельно обнаружить структуру в данных без подсказок-меток. В этом параграфе мы разобрали два центральных сюжета.
Восстановление матрицы (Netflix problem):
Гипотеза малоранговой факторизации A\approx UV^{\top} — из неё естественным образом получается рекомендательная система.
Функционал ошибки \sum_{(i,j)\in\Omega}(A_{ij}-\langle\mathbf u_i,\mathbf v_j\rangle)^{2} — здравый смысл, который при желании можно вывести из ОМП в предположении гауссова шума.
Идея малоранговой коррекции легла в основу LoRA — стандарта дообучения больших языковых моделей.
Автокодировщики:
Энкодер f сжимает вход в эмбеддинг \mathbf z малой размерности; декодер g восстанавливает вход из этого эмбеддинга. Минимизируется \left\lVert \mathbf x-g(f(\mathbf x)) \right\rVert^{2}.
В линейном случае оптимальное решение — ортогональная проекция на подпространство максимальной дисперсии (PCA).
Эмбеддинги — универсальное представление в современных моделях ИИ: компьютерное зрение, языковые модели, рекомендации.
Контрастивное обучение через сиамские сети позволяет обучать качественные эмбеддинги, не требуя восстанавливать вход целиком — достаточно знать, что одни пары похожи, а другие нет.
Оба сюжета фундаментально связаны через сингулярное разложение A=U\Sigma V^{\top} и теорему Эккарта–Янга–Мирского о наилучшей малоранговой аппроксимации. Это — мощный универсальный инструмент линейной алгебры, лежащий и за рекомендательными системами, и за PCA, и за современными тензорными разложениями (школа Тыртышникова–Оселедца).
Проверьте, что в примере 5.18 матрицы U и V действительно дают UV^{\top}, близкое к A во всех наблюдённых клетках.
Сколько параметров в модели A\approx UV^{\top}, если m=10^{6} пользователей, n=10^{4} фильмов, ранг r=20? Во сколько раз это меньше числа элементов исходной матрицы?
Покажите, что для центрированной выборки матрица C=\frac{1}{N}XX^{\top} действительно симметрична и неотрицательно определена.
Сделайте свой PCA: возьмите датасет \mathbf x_i\in\mathbb{R}^{2} вашего сочинения (хотя бы 30 точек), посчитайте C, найдите собственные значения и векторы (можно с помощью
numpy.linalg.eig). Постройте график как на рис. 5.24.Возьмите изображение в градациях серого. Получите его SVD (\texttt{numpy.linalg.svd}) и постройте малоранговые аппроксимации с r=1,5,20,50. Постройте график зависимости ошибки в норме Фробениуса от r. Найдите минимальное r, при котором ошибка не превышает 5\,\% от \left\lVert A \right\rVert_{F}.
^{\star} Докажите формулу (5.53), опираясь на SVD и ортогональность U,V в норме Фробениуса.
^{\star} Сиамская сеть для рукописных цифр. Возьмите MNIST/UCI Digits. Постройте сиамскую сеть, обучающуюся на парах «одна и та же цифра / разные цифры». После обучения покажите, что эмбеддинги одной цифры из разных написаний оказываются близкими.
5.5 Как искусственный интеллект меняет мир
В этом завершающем параграфе мы поднимемся над частными алгоритмами и посмотрим на современный ИИ как на явление, — разберёмся, как последние десять лет изменили картину мира, какие силы движут революцию больших моделей, что ИИ умеет, а что — пока нет, и как с большой осторожностью можно говорить о ближайшем будущем (рис. 5.29).
5.5.1 AlphaGo и второй прорыв в обучении с подкреплением
С точки зрения общественного восприятия, ИИ стал «настоящим ИИ» 9 марта 2016 г., в первый день матча AlphaGo против Ли Седоля. Если поражение Гарри Каспарова шахматному компьютеру Deep Blue в 1997 г. многие готовы были списать на «грубый перебор», то с игрой го такое объяснение не проходило: число позиций (\sim 10^{170}, что превосходит число атомов в наблюдаемой Вселенной) исключает любой полный перебор. Считалось, что компьютер достигнет человеческого уровня в го не раньше чем через 30 лет. AlphaGo, разработанная компанией DeepMind, обыграла Ли Седоля со счётом 4:1.
Архитектура AlphaGo. Под капотом находились две глубокие свёрточные нейронные сети:
сеть стратегии (policy network): по позиции выдавала распределение вероятностей следующего хода;
сеть оценки (value network): по позиции выдавала оценку шансов на победу.
Они комбинировались с поиском Монте-Карло по дереву (Monte Carlo Tree Search, MCTS) и обучались на партиях, сыгранных людьми, — своего рода имитационное обучение.
AlphaGo Zero: ИИ, который учился сам. Поразительным продолжением стала система AlphaGo Zero (октябрь 2017): она не использовала никаких человеческих партий, а училась только играя сама с собой — так называемое обучение с подкреплением через self-play. За 40 дней самообучения она достигла уровня, превосходящего AlphaGo 2016-го года, и в матче из 100 партий выиграла у предыдущей версии со счётом 100:0.
AlphaGo Zero показала: ИИ может учиться, взаимодействуя сам с собой, не нуждаясь ни в каких внешних данных — лишь имея формализованную задачу и среду, в которой можно играть. Это принципиально меняет масштаб возможного: данных не хватает — ИИ их сам производит. Этот принцип лежит в основе многих современных направлений: автоматического доказательства теорем, поиска новых материалов, проектирования лекарств. Развитие идеи привело к концепции мультиагентного ИИ — системы, где много обучаемых агентов взаимодействуют друг с другом, образуя экосистему, в которой коллективное знание растёт быстрее индивидуального.
5.5.2 Революция трансформеров и появление ChatGPT
Вторая, ещё более масштабная революция произошла в области обработки естественного языка — и она опирается на одну архитектурную идею.
В июне 2017 г. группа исследователей Google (Васвани, Шазир, Парма́р и др.) опубликовала статью Attention Is All You Need, в которой предложила архитектуру трансформер (Transformer). Идея была одновременно простой и революционной: отказаться от последовательной (рекуррентной) обработки текста и заменить её на механизм внимания, который позволяет каждому слову последовательности «смотреть» на любые другие слова одновременно. Это дало модели возможность учитывать сколь угодно далёкие зависимости в тексте и принципиально ускорило обучение за счёт параллельных вычислений на GPU.
Каждое слово (точнее — его эмбеддинг) превращается в три вектора: ключ (key), значение (value) и запрос (query). Затем для каждого слова считаются веса внимания — скалярные произведения его запроса с ключами всех остальных слов, пропущенные через softmax. Эти веса говорят, насколько каждое другое слово важно для текущего. Финальный эмбеддинг слова получается как взвешенное среднее значений (с этими весами).
Многократно повторяя эту операцию слоями, модель формирует всё более и более «смысловое» представление каждого слова в контексте.
От Transformer до ChatGPT.
BERT (Google, 2018) — 110 млн параметров; первая big-scale демонстрация мощи трансформеров на задачах NLP.
GPT-3 (OpenAI, 2020) — 175 млрд параметров; впервые показала, что увеличение масштаба ведёт к новым, ранее невидимым способностям (emergent abilities) — перевод, программирование, элементарная арифметика «из коробки», без специального обучения этим задачам.
ChatGPT (OpenAI, ноябрь 2022) — интерактивный чат-интерфейс к модели GPT-3.5, дообученной на диалогах с обратной связью (RLHF). За 5 дней — миллион пользователей; за 2 месяца — 100 миллионов (самый быстрорастущий сервис в истории интернета). С этого момента ИИ стал массовым потребительским явлением.
Отечественные модели. Россия — одна из немногих стран, имеющих собственную линейку больших языковых моделей промышленного уровня. ruGPT-3 (Сбер, 2020) была одной из первых крупных русскоязычных LLM в мире — 13 млрд параметров. GigaChat (Сбер, начиная с 2023) — семейство мультимодальных диалоговых моделей. YandexGPT (Яндекс, 2023+) — модель, оптимизированная под русский язык и интегрированная в продукты Яндекса (Алиса, Поиск, Браузер).
5.5.3 Законы скейлинга: формулы
К 2020-му году стало понятно, что качество LLM — удивительно предсказуемая функция трёх параметров: размера модели N (числа обучаемых параметров), размера датасета D (числа токенов в обучающей выборке) и вычислительного бюджета C (числа операций FLOPS, затраченных на обучение).
Закон Каплана (2020). Группа Каплана из OpenAI обнаружила, что кросс-энтропия (метрика качества языковой модели) убывает по степенному закону:
\tag{5.54} L(N)\;\approx\;L_\infty+\frac{A}{N^{\alpha}},
с \alpha\approx 0{,}076 (фиксируя D достаточно большим), где L_\infty — асимптотическая «непредсказуемая шумовая часть» естественного языка.
Закон Чинчилла (2022). Команда DeepMind (Хоффман и др., работа Training Compute-Optimal Large Language Models) уточнила картину: для оптимального обучения размер модели и размер датасета должны расти вместе, в определённой пропорции. Совместный закон:
\tag{5.55} \boxed{\;L(N,D)\;=\;L_\infty\;+\;\frac{A}{N^{0{,}34}}\;+\;\frac{B}{D^{0{,}28}},\;}
с константами L_\infty\approx 1{,}69, A\approx 406, B\approx 411 (в естественных нормировках).
Главный практический вывод — правило 20\colon 1. При фиксированном вычислительном бюджете C\sim 6\cdot N\cdot D (приближённая формула для трансформера) минимум функции L(N,D) достигается, когда
D \;\approx\; 20\cdot N.
Иными словами, на каждый параметр модели должно приходиться примерно 20 токенов обучающего текста. Это правило — неожиданное и критически важное: до его обнаружения почти все большие модели (включая GPT-3) были недообучены — слишком велики при слишком маленькой обучающей выборке (рис. 5.30).
Закон скейлинга объясняет, почему именно сейчас, в 2020-х, ИИ совершил качественный скачок. До 2017-го у нас не было архитектуры (трансформера), способной эффективно обучаться при таких масштабах. До 2010-х не было таких массивов данных (всемирная паутина накопила миллиарды документов). И не было дешёвых GPU, способных параллельно крутить миллиарды умножений матриц. Качество современных LLM определяется тремя факторами:
грамотными архитектурными идеями (внимание, residual connections, layer normalization, RLHF и пр.);
размером модели N;
объёмом и качеством обучающих данных D.
И всё это поверх классической оптимизации (5.21) и backpropagation (5.22).
5.5.4 AlphaFold и Нобелевская премия 2024 г.
Долгое время ИИ оставался занятием академическим. Прорывом, после которого никто уже не мог говорить о «модной игрушке», стала AlphaFold (DeepMind, 2018; AlphaFold 2 — 2020). Эта система решает задачу, считавшуюся одной из главных «50-летних» проблем биологии: по последовательности аминокислот предсказать трёхмерную структуру белка.
Почему это важно. Белки — работники клетки. От их трёхмерной формы зависит, как они работают: блокируют ли вирусную инфекцию, переносят ли кислород, разлагают ли пластик, вызывают ли болезнь Альцгеймера. Долгое время структуру белка определяли экспериментально — методами рентгеновской кристаллографии или криоэлектронной микроскопии. Расшифровка одного белка могла занимать годы, и стоила сотни тысяч долларов.
К концу 2010-х было известно \sim 170\,000 структур (база PDB), при том что в человеческом организме одних только белков несколько сотен тысяч.
Что сделала AlphaFold. В 2020 г. команда DeepMind во главе с Демисом Хассабисом и Джоном Джампером представила AlphaFold 2 на конкурсе CASP14 (Critical Assessment of protein Structure Prediction) — международном двухгодичном соревновании по предсказанию структур. Результат был поразительным: средняя точность \approx 87 по шкале GDT_TS, что сравнимо с экспериментальной точностью кристаллографии (см. рис. 5.31).
В 2021 г. DeepMind открыто опубликовала структуры \sim 200 миллионов белков — практически все известные белки в природе. К концу 2024 г. базой AlphaFold пользовались более 2 миллионов учёных из 190 стран.
Нобелевская премия 2024. 9 октября 2024 г. Нобелевская премия по химии была присуждена Демису Хассабису и Джону Джамперу (DeepMind) за создание AlphaFold, а также Дэвиду Бейкеру (Вашингтонский университет) за работы по компьютерному дизайну белков. День спустя — Нобелевская премия по физике была присуждена Джону Хопфилду и Джеффри Хинтону «за основополагающие открытия и изобретения, позволившие машинное обучение с помощью искусственных нейронных сетей». Впервые в истории Нобелевские премии по двум естественнонаучным номинациям в один и тот же год были связаны с ИИ.
Связь с трансформерами. В сердце AlphaFold лежит модификация архитектуры трансформера. Её основное отличие: вместо последовательности слов модель обрабатывает «последовательность аминокислот» (\sim 20 типов «букв»); аналогом внимания между словами выступает внимание между парами аминокислот в свёрнутом белке. Это и есть глубокая идея: задачу структурной биологии удалось редуцировать к языковой архитектуре. Аналогичный подход — редукция новой задачи к языку — сегодня применяется к десяткам новых областей: химия, материаловедение, геномика, программирование, математика.
5.5.5 Шкалирование во время вывода и цепочки рассуждений
К 2023–2024 гг. стало понятно: одного только увеличения размера модели для следующих качественных скачков может не хватить. Возникла новая идея — масштабирование на этапе вывода (inference-time scaling).
Цепочки рассуждений (Chain-of-Thought). Классическая большая языковая модель сразу выдаёт ответ. Оказалось, что если попросить модель сначала рассуждать пошагово вслух — выписать промежуточные шаги, потом сделать вывод — качество ответов резко возрастает, особенно в задачах математики, программирования, формальной логики. Этот феномен описан в работе Wei et al. (2022) и получил название цепочек рассуждений (chain-of-thought, CoT).
Reasoning-модели (o1, o3 и далее). В сентябре 2024 г. OpenAI представила модель o1 — reasoning model, специально обученную «думать дольше» перед ответом. На олимпиадных задачах по математике (AIME, IMO уровень), программированию (Codeforces) и физике она показала результаты, сравнимые с лучшими школьниками-олимпиадниками. В 2025-м — модель o3, более продвинутая. У других разработчиков появились свои аналоги — Claude Sonnet с extended thinking (Anthropic), Gemini Deep Think (Google DeepMind), DeepSeek-R1 (Китай) и многие другие.
Reasoning-модель тратит существенно больше вычислений в момент ответа на запрос, а не на этапе обучения. Технически это делается двумя способами: (а) обучить модель генерировать длинные цепочки промежуточных рассуждений (десятки или сотни тысяч токенов) перед финальным ответом; (б) явно искать в дереве возможных рассуждений (как MCTS в AlphaGo), оценивая разные ветки и выбирая лучшую.
Возникает новая ось скейлинга: время на запрос. Это принципиально — значит, дальнейший прогресс может идти не только за счёт увеличения моделей, но и за счёт более длительных рассуждений.
5.5.6 Воплощённый ИИ и физическая картина мира
При всех успехах языковых моделей у современного ИИ есть очевидные слабости. Главная — отсутствие у моделей того, что Янн ЛеКун (главный исследователь Meta AI и один из лауреатов премии Тьюринга 2018 г.) называет «здравым смыслом физического мира» (common sense). Даже самая большая LLM не знает простых вещей, известных любому двухлетнему ребёнку: если положить чашку на край стола и сдвинуть — она упадёт; если объект скрылся за дверью — он не перестал существовать; если предмет тяжёлый, его трудно поднять.
ЛеКун в работах 2022–2025 гг. предлагает альтернативную к LLM парадигму — world models (модели мира). Идея: ИИ должен учиться не на текстах (которые описывают мир), а на видео и сенсорных данных — так же, как учатся дети. Внутреннее представление мира получает структуру «причинно-следственного графа»: действие \to изменение состояния. Такой ИИ называется воплощённым (embodied AI). Он, возможно, не будет поражать литературными способностями, но сможет планировать, прогнозировать и взаимодействовать с физическим миром — то, чего нынешним моделям заметно не хватает.
Реализация этой идеи — одна из «святых граалей» современного ИИ. Если она удастся, мы получим ИИ, способный быть оператором робота, а не просто собеседником.
5.5.7 Исчерпание данных и пределы скейлинга
Закон Чинчилла говорит: чтобы получить более качественную модель, нужно одновременно увеличивать N и D. Но публично доступный интернет конечен. По разным оценкам, к 2026–2028 гг. передовые модели будут обучаться на всём качественном тексте, написанном человечеством в открытом доступе. Что дальше?
Возможные пути:
Синтетические данные — ИИ сам генерирует обучающую выборку. Уже сейчас часть обучающих данных моделей синтезируется другими моделями. Это путь к self-improving AI — системе, обучающейся самой на себе (по аналогии с AlphaGo Zero).
Мультимодальные данные — видео и аудио содержат на порядки больше информации, чем текст. Видео в сети (YouTube и др.) хватит на десятилетия обучения.
Алгоритмические улучшения — эффективные архитектуры могут учиться на меньших объёмах данных. Например, mixture of experts (MoE), retrieval-augmented generation (RAG), новые виды внимания, более эффективные оптимизаторы.
Вычислительные пределы. Закон Мура замедляется; энергетика обучения крупнейших моделей сравнима с потреблением небольших городов. Это уже инженерное ограничение.
5.5.8 Прогнозы и сингулярность по Курцвейлю
Рэй Курцвейль — американский инженер, футуролог, лауреат Национальной медали технологий США — знаменит своими долгосрочными прогнозами в области ИИ. Его взгляды изложены в двух книгах с почти одинаковыми названиями — они отличаются всего одной буквой:
The Singularity Is Near («сингулярность близка»), опубликована в 2005 году. В ней Курцвейль предсказал, что компьютер сравняется с человеком в задачах общего интеллекта к 2029 году, а сингулярность — момент кардинального ускорения технологического развития — наступит в 2045 году.
The Singularity Is Nearer («сингулярность ещё ближе»), опубликована почти через двадцать лет, в 2024 году. В этой новой книге Курцвейль подтвердил оба прогноза 2005 года, ссылаясь на ChatGPT и GPT-4 как на «доказательство по дороге к 2029-му».
Большинство экспертов считают эти прогнозы оптимистичными, но смены парадигмы — быстрого и кардинального изменения возможностей ИИ — сегодня не отрицает никто.
Artificial General Intelligence (AGI) — ИИ, превосходящий человека в любой интеллектуальной задаче. Сегодня все системы — узкие: AlphaGo гениально играет в го, но не знает, что такое кошка; GPT-5 пишет талантливые стихи, но не способен поработать поваром. AGI означал бы единую систему, которая может всё.
Прогнозы насчёт AGI разнятся от 2027 (Сэм Альтман) до конца 2040-х (многие академические эксперты). Само понятие AGI часто критикуют: проблема в том, что граница между «узким» и «общим» интеллектом расплывчата, а отдельные интеллектуальные навыки появляются у современных моделей постепенно, а не одновременно.
5.5.9 Меняющийся ландшафт профессий
ИИ принципиально меняет очень многие профессии. Главный тренд — не замена человека ИИ, а симбиоз: человек делает то, что у него хорошо получается (постановка задач, оценка результата, этические решения, общение, креативные прорывы), а рутинная интеллектуальная работа автоматизируется. Этот принцип хорошо называется «centaur» — кентавр, союз человека и машины.
Несколько примеров.
Программирование. Уже в 2025 г. значительная часть производственного кода в IT-компаниях пишется при участии или полностью с помощью моделей-кодеров (GitHub Copilot, Cursor, Claude Code, Codeium). Программист всё больше превращается в архитектора и ревьюера; рутинные функции, обвязки, тестовый код, документация — автоматизируются.
Наука. Уже сейчас существуют системы, способные самостоятельно формулировать научную гипотезу, ставить эксперимент (в симуляции), писать черновик статьи. Это AI Scientist — класс систем, появившийся в 2024–2025 гг. В материаловедении, химии, биологии счёт открытым новым материалам, структурам, реакциям с участием ИИ идёт на сотни тысяч в год.
Финансовый сектор. Уже более 70 % алгоритмической торговли ведётся ИИ-системами (включая системы на основе нейросетей). Анализ рисков, кредитный скоринг, обнаружение мошенничества, прогноз рынков — всё это давно автоматизировано.
Медицина. Системы помощи в диагностике: компьютерная томография анализируется ИИ; патология (распознавание раковых клеток на гистологических снимках); чтение ЭКГ. Уровень точности часто превосходит средневзвешенного врача.
Дизайн, киноиндустрия и музыка. Генеративные модели изображений (Midjourney, Stable Diffusion, DALL-E, FLUX) изменили рекламу, графический дизайн, начальные стадии разработки игр и кино. Музыка (Suno, Udio): системы, генерирующие профессионального уровня композиции по текстовому запросу. Эти технологии в основном работают вместе с человеком — как мощные кисти, а не как замена художника.
Юриспруденция. Анализ контрактов, поиск прецедентов, подготовка проектов документов — области, где ИИ за минуту делает работу, на которую тратился час или день.
Государственное управление и общественные сервисы. Цифровизация госуслуг (МФЦ, портал Госуслуг в России) — лишь начало. ИИ-ассистенты для граждан, помощь в обработке обращений, переводы официальных документов — всё это уже внедряется.
5.5.10 Учиться учиться: главный навык XXI века
Каков же универсальный совет тем, кто только начинает свой профессиональный путь? Парадоксально, ИИ-революция предъявляет к человеку гуманистические требования.
Что точно не заменится.
Постановка задач и целеполагание. ИИ может найти оптимальный путь к цели, но цель должен ставить человек.
Глубокое понимание области. ИИ многое знает поверхностно; специалист, который понимает, почему та или иная идея работает, остаётся незаменим как заказчик и ревьюер ИИ-инструментов.
Творчество и нетривиальные идеи. Революционные научные прорывы, художественные открытия, организация людей — области, где смелые, рисковые, неожиданные решения остаются за человеком.
Этика и ответственность. Решения, влияющие на жизни людей, должны принимать люди.
Что точно полезно.
Уметь быстро учиться. В мире, где появляются новые ИИ-инструменты каждый месяц, выигрывает тот, кто умеет быстро осваивать новое и применять его в своей работе.
Уметь работать с ИИ-инструментами. «Промпт-инженерия» (искусство задавать правильные вопросы моделям) и навык критической оценки результатов ИИ становятся базовыми.
Понимать математику. Парадоксально, но именно математика — наиболее «незаменимый» навык в эпоху ИИ. Кто понимает, как устроены модели, кто умеет читать их внутреннее устройство, кто может улучшать алгоритмы — тому открыты самые интересные позиции в этой отрасли.
ИИ не заменяет человека — он становится инструментом усиления человеческого интеллекта. Самая успешная команда XXI века — это человек, мастерски использующий ИИ, более продуктивный, чем команда того же размера без ИИ — но в десять или сто раз менее продуктивный, чем человек, который понимает, что и как делают эти системы. Книга, которую вы держите в руках, — небольшой шаг в эту сторону.
5.5.11 Что мы прошли в этой главе
Вернёмся к карте маршрута, нарисованного в начале главы. Мы вышли из точки «опрос на выборах» и пришли в точку «AGI и нобелевские лауреаты». На пути — один и тот же универсальный методологический принцип, проходящий красной нитью:
Принцип максимума правдоподобия: запиши вероятностную модель, выпиши функцию правдоподобия данных, максимизируй по параметрам — и получишь оптимальную оценку.
В простейших задачах ОМП даёт выборочное среднее (5.1). В линейной модели — метод наименьших квадратов (5.2). В классификации — кросс-энтропию (5.3). В рекомендательных системах и автокодировщиках — квадратичную невязку в матричных моделях (5.4). А в современных языковых моделях — то же самое: предсказание следующего токена — это ОМП для условного распределения \Pr(\text{токен}_t\mid\text{токены}_{1..t-1}).
Вокруг этого универсального принципа вращаются три кита анализа данных:
Вероятностный кит даёт постановку задачи: что мы хотим оценить и каковы шумы.
Оптимизационный кит даёт инструмент: градиентный спуск (5.21) с эффективной реализацией градиента через backpropagation (5.22).
Линейно-алгебраический кит даёт язык: матрично-векторные умножения внутри нейронных сетей (5.10); малоранговые аппроксимации (SVD, 5.21); тензорные разложения.
Этот аппарат, изложенный в одной главе на двух десятках страниц, — один и тот же — лежит и в основе академической математической статистики XX века, и в основе многомиллиардной индустрии ИИ XXI века. Между ними — лишь шкала: данных стало в миллиарды раз больше, моделей в миллиарды раз сложнее, вычислений в миллиарды раз быстрее. Но идеи — те же.
Заключение главы
Современный ИИ — не магия, а аккуратная математика, помноженная на инженерию, помноженную на огромные вычислительные ресурсы. Все ключевые идеи — от ОМП до backpropagation, от SVD до трансформеров — доступны старшекласснику. Чем глубже понимать эти идеи, тем меньше страха перед ИИ и тем больше возможностей применять его себе на пользу.
Мы живём в очень интересное время. Возможно, ваше поколение — первое поколение в истории, для которого партнёром в каждой интеллектуальной работе будет ИИ. Желаем вам быть в этом партнёрстве старшим партнёром — тем, кто понимает, кто ставит цели и кто принимает решения. А ИИ — помощник, инструмент усиления, как в своё время паровая машина, электричество, компьютер, интернет. Каждая из этих технологий перекраивала мир — но не отменяла главного, что делает нас людьми.
Посмотрите запись 2-й партии матча AlphaGo–Ли Седоль (Wikipedia: AlphaGo vs Lee Sedol). 37-й ход AlphaGo — знаменитый «нечеловеческий» ход. Прочитайте о нём и подумайте, почему он удивил даже мировых экспертов го.
Вычислите по закону Чинчилла (5.55) ожидаемую кросс-энтропию для модели на N=10^{11} параметров, обученной на D=2\cdot 10^{12} токенах. Сравните с моделью того же размера, но обученной на D=10^{11} токенах (типичный случай 2020 г.).
Найдите в открытых источниках три современных применения AlphaFold — например, в разработке новых лекарств. Опишите кратко.
Поразмышляйте над выбором профессии. Какие задачи вашей планируемой будущей деятельности можно делегировать ИИ уже сегодня? Какие — останутся за вами навсегда? Сформулируйте ответ в виде короткого эссе.
^{\star} Прочитайте оригинальную статью Attention Is All You Need (2017, Васвани и др.; arXiv:1706.03762). Попробуйте понять, как именно устроен механизм внимания \mathrm{softmax}(QK^{\top}/\sqrt{d})\,V. Запишите это формулой, объяснив роль матриц Q, K, V.
^{\star} Запустите локально маленькую LLM (например, через Ollama: модели Llama, Qwen, GigaChat-Mini, YandexGPT-Lite). Сравните качество её ответов на ваши вопросы. Какие задачи она выполняет хорошо? Где ошибается? Замечаете ли вы предсказуемые сценарии её слабости?
Также встречается название logloss или отрицательное лог-правдоподобие.↩︎