§3. Вычисление главных компонент

Следующий простой алгоритм (степенной метод) позволяет вычислить для симметричной (не обязательно неотрицательно определённой) матрицы $\boldsymbol A$ несколько наибольших по модулю собственных значений $\lambda_l$ и соответствующих собственных векторов $\boldsymbol c_l$. Сначала разберём частный случай этого алгоритма — вычисление $\lambda_1$ и $\boldsymbol c_1$.

Возьмём произвольный ненулевой вектор $\boldsymbol x_0$ из $\mathbb R^m$. Пусть выполнены два условия: $|\lambda_1| > |\lambda_2| \ge \dots \ge |\lambda_m|$ и $\boldsymbol c_1^T\boldsymbol x_0 \ne 0$. Положим $\boldsymbol x_{i+1} = \boldsymbol A\boldsymbol x_i$ (то есть $\boldsymbol x_i = \boldsymbol A^i\boldsymbol x_0$, что объясняет название метода). Тогда последовательность значений $t_{i+1} = \boldsymbol x_i^T\boldsymbol x_{i+1}/|\boldsymbol x_i|^2$ сходится к $\lambda_1$ со скоростью геометрической прогрессии со знаменателем $\gamma = |\lambda_2/\lambda_1| < 1$.

Доказательство. Обозначим через $b_1, \dots, b_m$ координаты вектора $\boldsymbol x_0$ в базисе $\boldsymbol c_1, \dots, \boldsymbol c_m$, то есть $\boldsymbol x_0 = b_1\boldsymbol c_1 + \dots + b_m\boldsymbol c_m$. Так как $\boldsymbol A\boldsymbol c_l = \lambda_l\boldsymbol c_l$, то $\boldsymbol A^i\boldsymbol c_l = \lambda_l^i\boldsymbol c_l$. Отсюда

где $\boldsymbol\delta_i = b_2(\lambda_2/\lambda_1)^i\boldsymbol c_2 + \dots + b_m(\lambda_m/\lambda_1)^i\boldsymbol c_m$, причём в силу первого условия $|\boldsymbol\delta_i| = O(\gamma^i) \to 0$ при $i \to \infty$.

В частности, когда все $\lambda_l > 1$, это означает, что при каждой итерации длина вектора $\boldsymbol x_i$ увеличивается, но быстрее остальных растёт его координата по оси $\boldsymbol c_1$. При этом направление $\boldsymbol x_i$ приближается к направлению первой главной оси.

Формально это можно доказать с помощью известного из линейной алгебры неравенства Коши — Буняковского — Шварца $|\boldsymbol x^T\boldsymbol y| \le |\boldsymbol x||\boldsymbol y|$. Согласно второму условию алгоритма $b_1 = \boldsymbol c_1^T\boldsymbol x_0 \ne 0$. Отсюда $|\boldsymbol x_i| = (\boldsymbol x_i^T\boldsymbol x_i)^{1/2} = |\lambda_1|^i[|b_1| + O(\gamma^i)]$; $\boldsymbol e_i = \boldsymbol x_i/|\boldsymbol x_i| = (\operatorname{sign}\lambda_1)^i(\operatorname{sign}b_1)\boldsymbol c_1 + \boldsymbol r_i$, где $|\boldsymbol r_i| = O(\gamma^i)$; $t_{i+1} = \lambda_1 + O(\gamma^i)$. $\blacksquare$

По кирпичикам: каждое умножение на $\boldsymbol A$ растягивает координату вдоль $\boldsymbol c_l$ в $\lambda_l$ раз. После $i$ умножений координата по $\boldsymbol c_1$ обгоняет остальные в $(\lambda_1/\lambda_l)^i$ раз — экспоненциально. Остальные «тонут», и вектор разворачивается вдоль $\boldsymbol c_1$.

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

Замечание 4. Если $|\lambda_1| > 1$, то $|\boldsymbol x_i| \to \infty$, а если $|\lambda_1| < 1$, то $|\boldsymbol x_i| \to 0$ при $i \to \infty$. При вычислении на компьютере в первом случае возможно переполнение, во втором — исчезновение порядка. Предотвратить эти ситуации позволяет нормировка векторов $\boldsymbol x_i$ на каждом шаге.

Учитывая замечание 4, приходим к следующему алгоритму для вычисления $\lambda_1$ и $\boldsymbol c_1$.

1. Положим $i = 0$, $t_0 = 0$ и зададим произвольный ненулевой вектор $\boldsymbol x_0$ из $\mathbb R^m$. Нормируем $\boldsymbol x_0$: пусть $\boldsymbol e_0 = \boldsymbol x_0/|\boldsymbol x_0|$ (обычно сразу берут $\boldsymbol e_0 = (1, 0, \dots, 0)$).
2. Вычислим $\boldsymbol x_{i+1} = \boldsymbol A\boldsymbol e_i$ и $t_{i+1} = \boldsymbol e_i^T\boldsymbol x_{i+1}$.
3. Нормируем $\boldsymbol x_{i+1}$: $\boldsymbol e_{i+1} = \boldsymbol x_{i+1}/|\boldsymbol x_{i+1}|$.
4. Если $|t_{i+1} - t_i| \le \varepsilon$, где $\varepsilon$ — достаточно малое число, то положим $\lambda_1 = t_{i+1}$ и $\boldsymbol c_1 = \boldsymbol e_{i+1}$ и закончим итерационный процесс. Иначе увеличим $i$ на 1 и перейдём к шагу 2.

Замечание 5. Если при выборе $\boldsymbol x_0$ случится, что $\boldsymbol c_1^T\boldsymbol x_0 = 0$, то через несколько итераций за счёт погрешностей округления, как правило, появится ненулевая координата у $\boldsymbol x_i$ в направлении $\boldsymbol c_1$, и процесс, хотя и с некоторым запаздыванием, выйдет на первое собственное значение. Можно также стартовать с разных $\boldsymbol x_0$ и сравнивать результаты.

Следующие компоненты: ортогонализация

Пусть $|\lambda_1| > |\lambda_2| > \dots > |\lambda_k| > |\lambda_{k+1}| \ge \dots \ge |\lambda_m|$. Для вычисления $k > 1$ первых главных компонент кроме умножения на $\boldsymbol A$ и нормировки понадобится ещё ортогонализация. Именно, предположим, что $\lambda_l$ и $\boldsymbol c_l$ уже известны при всех $l \le k - 1$. Чтобы найти $\lambda_k$ и $\boldsymbol c_k$, выполним следующие шаги.

1. Положим $i = 0$, $t_0 = 0$ и зададим произвольный вектор $\boldsymbol x_0$ из $\mathbb R^m$ (надо только, чтобы $\boldsymbol c_k^T\boldsymbol x_0 \ne 0$).
2. Проведём ортогонализацию $\boldsymbol x_i$ к векторам $\boldsymbol c_1, \dots, \boldsymbol c_{k-1}$: $\boldsymbol y_i = \boldsymbol x_i - (\boldsymbol c_1^T\boldsymbol x_i)\boldsymbol c_1 - \dots - (\boldsymbol c_{k-1}^T\boldsymbol x_i)\boldsymbol c_{k-1}$.
3. Нормируем $\boldsymbol y_i$: $\boldsymbol e_i = \boldsymbol y_i/|\boldsymbol y_i|$.
4. Вычислим $\boldsymbol x_{i+1} = \boldsymbol A\boldsymbol e_i$ и $t_{i+1} = \boldsymbol e_i^T\boldsymbol x_{i+1}$.
5. Нормируем $\boldsymbol x_{i+1}$: $\boldsymbol z_{i+1} = \boldsymbol x_{i+1}/|\boldsymbol x_{i+1}|$.
6. Если $|t_{i+1} - t_i| \le \varepsilon$, то положим $\lambda_k = t_{i+1}$, $\boldsymbol c_k = \boldsymbol z_{i+1}$ и закончим процесс. В противном случае приравняем $\boldsymbol x_{i+1} = \boldsymbol z_{i+1}$, увеличим $i$ на 1 и перейдём к шагу 2.

Смысл по-пацански: самого громкого оратора уже нашли — выведем его из зала (вычтем проекцию на $\boldsymbol c_1$). Тогда толпа соберётся у второго по громкости. И так далее.

Для матрицы корреляций собак и волков $\gamma = |\lambda_2/\lambda_1| = 0{,}883/4{,}100 \approx 0{,}22$, и $\lambda_1$ находится за 8 шагов. Для $\lambda_2$ скорость определяется отношением $|\lambda_3/\lambda_2| = 0{,}639/0{,}883 \approx 0{,}72$ — шагов нужно уже около 40. Чем ближе соседние собственные значения, тем медленнее степенной метод.

(Матрица рис. 3 округлена до сотых, поэтому собственные значения по ней чуть отличаются от точных 4,100; 0,883; 0,639.)

Насколько доверять собственным значениям

В заключение обсудим, насколько следует доверять вычисленным значениям $\lambda_l$ и $\boldsymbol c_l$. Следующий пример показывает, что для числа объектов $n$ порядка нескольких десятков разброс этих значений может оказаться достаточно большим.

Пример 2. Статистика главных компонент для нормальной модели. Предположим, что объекты представляют собой выборку размера $n$ из закона $\mathcal N(\boldsymbol 0, \boldsymbol\Sigma)$ с невырожденной матрицей $\boldsymbol\Sigma$, у которой все собственные значения $\lambda_l$ ($l = 1, \dots, m$) различны. Тогда:

1. собственные значения $\hat\lambda_1, \dots, \hat\lambda_m$ и отвечающие им собственные векторы $\hat{\boldsymbol c}_1, \dots, \hat{\boldsymbol c}_m$ выборочной ковариационной матрицы $\hat{\boldsymbol\Sigma}$ являются оценками максимального правдоподобия для соответствующих $\lambda_1, \dots, \lambda_m$, $\boldsymbol c_1, \dots, \boldsymbol c_m$ (см. пример 2 главы 19) и, следовательно, обладают асимптотической эффективностью (см. §4 главы 9);
2. $\sqrt n(\hat\lambda_l - \lambda_l) \xrightarrow{d} \xi \sim \mathcal N(0, 2\lambda_l^2)$ при $n \to \infty$;
3. $\sqrt n(\hat{\boldsymbol c}_l - \boldsymbol c_l) \xrightarrow{d} \boldsymbol\zeta \sim \mathcal N(\boldsymbol 0, \boldsymbol V)$ при $n \to \infty$, где элементами матрицы $\boldsymbol V$ служат $v_{kl} = \lambda_l\sum_{k\ne l}|c_l|^2\lambda_k/(\lambda_k - \lambda_l)^2$.

Построим на основе свойства 2 асимптотический доверительный интервал (см. §1 главы 11) для $\lambda_1$ в примере 1. Так как свойство 2 равносильно сходимости $\sqrt{n/2}(\hat\lambda_l/\lambda_l - 1) \xrightarrow{d} \xi/(\lambda_l\sqrt2) \sim \mathcal N(0, 1)$, то при достаточно больших $n$ приближённо с вероятностью $1 - \alpha$ выполняется неравенство

где $x_{1-\alpha/2}$ — квантиль уровня $(1 - \alpha/2)$ закона $\mathcal N(0, 1)$ (см. §3 главы 7). Для $\alpha = 5\,\%$ по таблице T2 находим, что $x_{1-\alpha/2} = 1{,}96$. При $n = 43$ получаем 95 %-й доверительный интервал для $\lambda_1$ от 2,88 до 7,1.

По кирпичикам: $\sqrt{2/43} \approx 0{,}216$, $1{,}96\cdot0{,}216 \approx 0{,}423$. Тогда $4{,}100/1{,}423 \approx 2{,}88$ и $4{,}100/0{,}577 \approx 7{,}10$.

Смысл по-пацански: «первая компонента объясняет 68 %» звучит солидно, но на 43 черепах сама $\lambda_1$ гуляет от 2,9 до 7,1 — в два с половиной раза. Красивые проценты на малых выборках — с большим люфтом.

Проверь себя