§6. Нелинейные методы понижения размерности

Круговое движение первее прямолинейного: оно проще и более совершенно.

— Аристотель

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

Пример 7. Проецирование вершин многомерного симплекса. Возьмём 16 групп точек по 4 точки в группе. Каждая группа генерировалась с помощью датчика случайных чисел вблизи одной из 16 вершин правильного симплекса в $\mathbb R^{15}$ (см. формулу (15) в §4). У Лагутина на первую (из пятнадцати) компоненту приходится 8,3 % следа ковариационной матрицы, на вторую — 7,6 %, на обе — лишь 15,9 %. Проецирование на плоскость двух первых главных компонент не позволяет сохранить разделение на группы.

По кирпичикам: у правильного симплекса нет «длинных» и «коротких» направлений — все 15 осей равноправны, на каждую приходится около $1/15 \approx 6{,}7\,\%$ разброса. Две любые оси видят лишь $2/15$ картины, а остальные 13 «проекций» схлопываются друг на друга.

Смысл по-пацански: 16 компаний, каждая на своём углу, и все углы одинаково далеко друг от друга. Сфотографировать такое в 2D честной проекцией невозможно — на фото они встанут в одну кучу. Нужна не проекция, а «нарисовать по памяти».

Метод Сэммона

Для анализа подобных сложных конфигураций применяются нелинейные методы понижения размерности. Они заключаются в минимизации некоторой функции $F$, выражающей суммарное расхождение между заданными различиями объектов $d_{ij}$ и расстояниями $\delta_{ij}$ между образами объектов в подпространстве небольшой размерности.

Наиболее простой из таких функций является

Широкое распространение получил метод Сэммона (предложен Дж. Сэммоном в 1969 г.):

который обладает свойством более точно передавать небольшие различия и менее точно — большие, так как при отображении бо́льших расстояний допустимы бо́льшие ошибки. Деление на константу $C_1$ нужно для инвариантности $F_1$ к изменению масштаба величин $d_{ij}$.

По кирпичикам: ошибка на паре делится на $d_{ij}$. Промахнулся на 1 при истинном расстоянии 2 — штраф $1/2$. Промахнулся на 1 при расстоянии 100 — штраф $1/100$. Сэммон бережёт ближних соседей: локальная структура важнее, чем точные дальние расстояния.

Метод Сэммона, стартовав с проекции на главные компоненты, «расталкивает» группы, и разделение на группы сохраняется на плоскости. В нашей модели (разброс вокруг вершин $\sigma = 0{,}05$) значение $F_1$ уменьшается примерно с 0,42 до 0,13. У Лагутина — с 0,44 до 0,13; разброс точек в книге не указан, поэтому мы его подобрали.

Для поиска минимума функции $F$ прежде всего строится начальная конфигурация точек в пространстве небольшой размерности $m$. Это могут быть просто проекции в подпространство, образованное какими-либо из $m$ координатных осей, или в подпространство $m$ первых главных компонент, или, наконец, случайная конфигурация, в которой $nm$ координат получены с помощью датчика случайных чисел.

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

Метод сопряжённых градиентов

1. Определяется направление $\boldsymbol p_t$ в пространстве $\mathbb R^{nm}$ по формулам

где $\boldsymbol g_t = \Bigl(\frac{\partial F}{\partial x_{11}}, \dots, \frac{\partial F}{\partial x_{nm}}\Bigr)$ обозначает градиент минимизируемой функции $F$, $\boldsymbol p_{t-1}$ — направление на предыдущем шаге, коэффициент $\beta_t = |\boldsymbol g_t|^2/|\boldsymbol g_{t-1}|^2$.

В задаче 2 предлагается проверить, что координаты градиента функции $F_1$, заданной формулой (18), удовлетворяют соотношениям

где $i = 1, \dots, n$; $l = 1, \dots, m$. Считается, что все объекты различны ($d_{ij} \ne 0$). В случае совпадения образов объектов в $\mathbb R^m$ ($\delta_{ij} = 0$) соответствующее слагаемое полагается равным некоторой константе $c \ne 0$, чтобы образы могли «разойтись» в процессе минимизации. 2. Производится перемещение до точки минимума по выбранному направлению (например, с помощью рассматриваемого ниже алгоритма золотого сечения): $\boldsymbol x_{t+1} = \boldsymbol x_t + \alpha_t\boldsymbol p_t$, где $F(\boldsymbol x_t + \alpha_t\boldsymbol p_t) = \min_{\alpha\ge0}F(\boldsymbol x_t + \alpha\boldsymbol p_t)$. 3. Если $|\boldsymbol g_{t+1}| \le \varepsilon$, где $\varepsilon$ — достаточно малое число (скажем, $\varepsilon = 0{,}001$), то поиск заканчивается, иначе $t$ увеличивается на 1, и происходит возвращение к пункту 1.

Впервые метод сопряжённых градиентов был опубликован в качестве численного алгоритма для решения системы линейных уравнений $\boldsymbol A\boldsymbol x = \boldsymbol b$ с положительно определённой матрицей $\boldsymbol A$. Эта задача равносильна минимизации по $\boldsymbol x$ квадратичной функции $F(\boldsymbol x) = \frac12\boldsymbol x^T\boldsymbol A\boldsymbol x - \boldsymbol b^T\boldsymbol x$. Оказывается, что (при отсутствии погрешностей округления) метод сопряжённых градиентов позволяет найти точку минимума $F(\boldsymbol x)$ не более чем за $k$ итераций, где $k$ — размерность $\boldsymbol x$. При этом направления одномерной минимизации $\boldsymbol p_1, \dots, \boldsymbol p_k$ являются взаимно сопряжёнными относительно матрицы $\boldsymbol A$: $\boldsymbol p_i^T\boldsymbol A\boldsymbol p_j = 0$ при всех $i \ne j$, что объясняет название метода.

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

1. на каждом шаге осуществляется минимизация по выбранному направлению, что ценно ввиду большого числа $nm$ оптимизируемых переменных;
2. в отличие от метода Ньютона не требуется память для хранения и не расходуется время на вычисление матрицы вторых производных размера $nm \times nm$;
3. в отличие от наискорейшего спуска (см. пример 8) метод сопряжённых градиентов обычно удовлетворительно работает при минимизации «овражных» функций.

Смысл по-пацански: градиентный спуск — как мяч, пущенный в узкий овраг: каждый шаг «вниз по склону», и его кидает от стенки к стенке. Сопряжённые градиенты помнят предыдущий шаг и потому идут вдоль дна.

Пример 8. Критика метода наискорейшего спуска

Иногда в качестве направления минимизации на очередном шаге рекомендуют выбирать антиградиент $-\boldsymbol g_t$, поскольку в этом направлении функция убывает быстрее всего. (В частности, для поверхности, «горы», задаваемой гладкой функцией двух переменных, антиградиент указывает направление локально наискорейшего спуска.) Однако эта «жадность» метода приводит к огромному числу смен направления при движении вдоль «дна оврага». Наглядный пример даёт поиск минимума функции Розенброка

которая определяет в трёхмерном пространстве поверхность типа серповидного ущелья с точкой минимума $(1, 1)$.

У Лагутина процесс минимизации $F(x, y)$ методом наискорейшего спуска из точки $(-6/5, 1)$ надолго застрял вблизи начала координат. Не помогла бы и счастливая случайность — на одной из итераций процедура одномерного поиска неожиданно нашла второй минимум по направлению. Затем было выполнено ещё несколько сотен шагов, но ощутимого изменения целевой функции это не дало, и счёт был прерван после 1000 итераций вдали от искомого решения. Метод сопряжённых градиентов, в свою очередь, позволил дойти до точки минимума за два десятка шагов. Хорошо виден его циклический характер в случае двух переменных: движение вдоль дна ущелья подправляется перемещением по антиградиенту.

В нашей реализации (одномерный поиск — удвоение шага и золотое сечение до машинной точности) наискорейший спуск за 1000 шагов всё же добирается до точки около $(0{,}91;\ 0{,}83)$, но так и не доходит до $(1, 1)$. Сопряжённым градиентам (с перезапуском по антиградиенту каждые 2 шага — отсюда и «циклический характер», о котором пишет книга) хватает 33 шагов. Насколько далеко застрянет спуск, зависит от деталей одномерного поиска, а вывод тот же: на «оврагах» наискорейший спуск безнадёжно медленный.

Метод золотого сечения

Для полноты изложения приведём один из алгоритмов минимизации по направлению. Алгоритм включает следующие пункты.

1. Необходимо проверить, что $\boldsymbol p_t$ — это направление убывания $F$. Для этого зададим малый шаг $\delta > 0$ (скажем, $\delta = 0{,}01$) и вычислим при $\alpha = \delta$ значение $\varphi(\alpha) = F(\boldsymbol x_t + \alpha\boldsymbol p_t)$. Если $\varphi(\delta) > \varphi(0)$, то заменим направление $\boldsymbol p_t$ на антиградиент $-\boldsymbol g_t$ и будем уменьшать величину шага делением пополам до тех пор, пока не добьёмся убывания функции.

2. Определим отрезок $[a, b]$, содержащий точку минимума функции $\varphi(\alpha)$. Положим $a = 0$, а для нахождения $b$ будем увеличивать шаг вдвое до тех пор, пока функция не прекратит убывать. (Для заданных выше $F_0$ или $F_1$ это обязательно произойдёт, поскольку они стремятся к $+\infty$ при $|\boldsymbol x| \to \infty$.)

   При поиске точки минимума $\alpha^*$ (непрерывной) функции $\varphi(\alpha)$ методом золотого сечения предполагается, что $\varphi$ унимодальна, то есть строго убывает слева от $\alpha^*$ и строго возрастает справа от $\alpha^*$. Для такой функции положение точки минимума можно уточнить, вычислив $\varphi$ в двух внутренних точках отрезка $[a, b]$. Например, если $\varphi(r) < \varphi(s)$ при $r < s$, то $\alpha^*$ обязана находиться на отрезке $[a, s]$. (Если нет унимодальности, то в результате применения метода золотого сечения будет найден один из локальных минимумов.)

3. На отрезке $[a, b]$ возьмём две пробные точки $r = a + (1 - \varkappa)(b - a)$ и $s = a + \varkappa(b - a)$, где $\varkappa = (\sqrt5 - 1)/2 \approx 0{,}618$ — золотое сечение (см. решение задачи 5 главы 1), и вычислим значения $\varphi(r)$ и $\varphi(s)$.

4. Если $\varphi(r) < \varphi(s)$, то полагаем $\alpha_t = r$, $a' = a$, $r' = a + (1 - \varkappa)(s - a)$, $s' = r$, $b' = s$ и вычисляем $\varphi(r')$. В противном случае полагаем $\alpha_t = s$, $a' = r$, $r' = s$, $s' = a + \varkappa(b - r)$, $b' = b$ и вычисляем $\varphi(s')$.

5. Если $b' - a' \le \varepsilon$ (скажем, $\varepsilon = 0{,}01$), то заканчиваем поиск, иначе убираем у всех переменных штрихи и возвращаемся к пункту 4.

Замечательно, что точки $r$ и $s$ осуществляют золотые сечения не только отрезка $[a, b]$, но и отрезков $[a, s]$ и $[r, b]$ соответственно. Это позволяет вычислять всего одно значение $\varphi$ при следующей итерации. При каждом выполнении пункта 4 длина отрезка локализации точки минимума уменьшается в $1 + \varkappa \approx 1{,}618$ раз. Всего за 10 повторов первоначальный отрезок сожмётся в $(1 + \varkappa)^{10} \approx 123$ раза.

По кирпичикам: почему одна точка «переиспользуется». $\varkappa$ — корень уравнения $\varkappa^2 = 1 - \varkappa$. Новый отрезок $[a, s]$ имеет длину $\varkappa(b - a)$, а старая точка $r$ отстоит от $a$ на $(1 - \varkappa)(b - a) = \varkappa^2(b - a)$, то есть ровно на долю $\varkappa$ новой длины — как раз там, где должна стоять новая точка $s'$.

Замечание 8. В отличие от линейного шкалирования Торгерсона (§4) метод Сэммона позволяет обрабатывать матрицы различий с пропусками. Для этого суммирование в формуле (18) достаточно проводить только для тех пар объектов, у которых различия измерены. Экспериментально установлено, что качество восстановления конфигурации будет почти таким же, как для полной матрицы, даже если доля пропусков составляет около 30 %.

Пропуски в таблицах данных

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

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

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

Эти недостатки влекут смещённость и несостоятельность оценок параметров моделей при обработке заполненной таблицы данных.

В качестве альтернативы рассмотрим

Метод локального заполнения

Пусть у $i$-й строки таблицы $\boldsymbol X = \|x_{il}\|_{n\times m}$ есть данные в столбцах с номерами $l_1, \dots, l_k$, $k < m$. Обозначим через $S_i$ множество строк $\boldsymbol X$, не имеющих пропусков в столбцах $l_1, \dots, l_k$. Для заполнения пропуска в $j$-й ячейке $i$-й строки выделим в $S_i$ подмножество $R_{ij}$ тех строк, у которых нет пропуска в $j$-й ячейке. Если $R_{ij}$ пусто, то заполнение невозможно ($i$-ю строку придётся исключить из обработки). В противном случае найдём в $R_{ij}$ строку, ближайшую к $i$-й в смысле, скажем, евклидова расстояния в $k$-мерном подпространстве, порождённом столбцами с номерами $l_1, \dots, l_k$ (если таких строк окажется несколько, то выберем одну из них случайно). Заполним пропуск значением из $j$-й ячейки найденной строки.

Смысл по-пацански: не знаешь рост пацана — не пиши «средний рост по району». Найди того, кто на него больше всех похож по остальным данным, и возьми рост у него. Так не создаётся толпы искусственных «середнячков».

Какими статистическими свойствами обладает этот метод? Будем рассматривать строки таблицы как реализации случайных векторов $\boldsymbol X_1, \dots, \boldsymbol X_n$. Пусть случайная величина $M_{il}$ — индикатор наличия пропуска в $l$-м столбце $i$-й строки. Положим $\boldsymbol M_i = (M_{i1}, \dots, M_{im})$, $i = 1, \dots, n$. Допустим, что:

1. случайные векторы $(\boldsymbol X_1, \boldsymbol M_1), \dots, (\boldsymbol X_n, \boldsymbol M_n)$ независимы и одинаково распределены;
2. $\boldsymbol M_1$ не зависит от $\boldsymbol X_1$ (то есть пропуски не зависят от значений);
3. распределение $\boldsymbol X_1$ имеет $m$-мерную плотность (см. П8);
4. $\mathbf P(M_{1l} = 1) = p_l < 1$ для всех $l = 1, \dots, m$.

Обозначим через $F(\boldsymbol x) = F(x_1, \dots, x_m)$ функцию распределения вектора $\boldsymbol X_1 = (X_{11}, \dots, X_{1m})$: $F(\boldsymbol x) = \mathbf P(X_{1l} \le x_l,\ l = 1, \dots, m)$, а через $\tilde F_n(\boldsymbol x)$ — эмпирическую функцию распределения, построенную по полученной методом локального заполнения таблице $\tilde{\boldsymbol X}$:

Тогда при выполнении указанных выше условий 1–4

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

При $\rho = 0{,}7$ и 30 % пропусков заполнение средним занижает коэффициент корреляции примерно с 0,70 до 0,60, а стандартное отклонение $Y$ — примерно с 0,97 до 0,82. Локальное заполнение сохраняет и то и другое.

Проверь себя