Модифицированная ортогонализация Грама-Шмидта

Материал из MachineLearning.

Модифицированная ортогонализация Грама-Шмидта - известный алгоритм ортогонализации, который строит ортогональные векторы $g_1, . . . , g_n$ , линейная оболочка которых совпадает с линейной оболочкой $f_1, . . . , f_n$ .

Содержание

1 Введение
2 Метод ортогонализации Грама-Шмидта
3 Вспомогательные утверждения
4 Алгоритм 1.1. Ортогонализация Грама-Шмидта
5 Модифицированная ортогонализация Грама-Шмидта
6 Алгоритм 1.2. Решение линейной задачи наименьших квадратов путём ортогонализации Грама-Шмидта с последовательным добавлением признаков
7 Литература

Введение

Рассмотрим ортогональное разложение $F = GR$ , где $R$ - верхняя треугольная матрица размера $n$ × $n$ , $G$ - ортогональная $l$ × $n$ - матрица, $G^TG = I_n$ . Для любой матрицы $F$ существует бесконечно много разложений указанного вида. Имея одно из них, легко выразить псевдообратную матрицу $F^+$ через $G$ и $R$ :

$F^+ = (R^TG^TGR)^{-1}R^TG^T = R^{-1}R^{-1T}R^TG^T = R^{-1}G^T$ .

Для вычисления псевдообратной $F^+$ достаточно построить какое-нибудь ортогональное разложение матрицы $F$ , обратить верхнюю треугольную матрицу $R$ и умножить её на $G^T$ . Этот метод во многом удобнее явного обращения матрицы.

Метод ортогонализации Грама-Шмидта

Для построения разложения $F = GR$ воспользуемся процессом ортогонализации Грама-Шмидта. Запишем матрицы $F$ и $G$ по столбцам:

$F = (f_1, . . . , f_n)$ ;

$G = (\tilde{g}_1, . . . , \tilde{g}_n)$ .

Волной здесь и далее обозначается операция нормирования вектора:

$\tilde{v} = v/||v||$ .

Процесс ортогонализации Грама-Шмидта строит ортогональные векторы $g_1, . . . , g_n$ , линейная оболочка которых совпадает с линейной оболочкой $f_1, . . . , f_n$ :

$g_1 = f_1$ ;

$g_2 = f_2 - \tilde{g}_1\tilde{g}^T_1f_2$ ;

· · ·

$g_j = f_j - \tilde{g}_1\tilde{g}^T_1f_j - ... - \tilde{g}_{j-1}\tilde{g}^T_{j-1}f_j$ .

Легко проверить, что для всех $k, j$ из $\{1, . . . , n\}$ , $k \ne j$ , векторы $\tilde{g}_k$ и $\tilde{g}_j$ ортогональны.

Вспомогательные утверждения

Лемма 1.1. На $j$ -м шаге процесса, $j = 1, . . . , n$ , матрица $F_j = (f_1, . . . , f_j)$ представима в виде ортогонального разложения $F_j = G_jR_j$ , где

$G_j = (\tilde{g}_1, . . . , \tilde{g}_j)$ - ортонормированная матрица;

$R_j = \begin{pmatrix} r_{11} & \cdots & r_{1j} \\ & \ddots & \vdots \\ 0 & & r_{jj} \end{pmatrix}$ - верхняя треугольная матрица, $r_{ij} = \left\{ \begin{eqnarray} \tilde{g}^T_if_j, & i < j;\\ ||g_j||, & i=j. \end{eqnarray}$

По окончании процесса ортогонализации получаем ортогональное разложение $F = G_nR_n$ .

С вычислительной точки зрения процесс Грама-Шмидта удобен тем, что на каждом шаге матрицы $G_j$ и $R_j$ получаются путём дописывания справа по одному столбцу к матрицам $G_{j-1}$ и $R_{j-1}$ соответственно. При этом предыдущие столбцы не изменяются (если не считать изменением дописывание нулей снизу к матрице $R_j$ - при разумной программной реализации эти нули всё равно не хранятся).

В следующей лемме утверждается, что обратная матрица $T_j = R^{-1}_j$ также является верхней треугольной и формируется путём дописывания столбцов справа.

Лемма 1.2. Пусть матрицы $R_j$ невырождены и в блочной записи имеют вид

$R_1 = (r_{11})$ ;

$R_j = \begin{pmatrix} R_{j-1} & & r_{j} \\ 0 & & r_{jj} \end{pmatrix}$ , $j = 2, . . . , n$ ,

где $r_{jj}$ - скаляр, $r_j$ - вектор-столбец размера $(j-1)$ . Тогда матрицы $T_j = R^{-1}_j$ могут быть вычислены по рекуррентной формуле

$T_1 = (t_{11})$ ;

$T_j = \begin{pmatrix} T_{j-1} & & t_{j} \\ 0 & & t_{jj} \end{pmatrix}$ , $j = 2, . . . , n$ ,

где $t_{jj} = 1/r_{jj}$ - скаляр, $t_j = -t_{jj}T_{j-1}r_j$ - вектор-столбец размера $(j - 1)$ .

Замечание 1.1. Обеспечить невырожденность матрицы $R_j$ в процессе ортогонализации очень просто. Допустим, матрица $R_{j-1}$ невырождена. Поскольку $R_j$ - верхняя треугольная, вырожденность может возникнуть только в том случае, если $r_{jj} = 0$ . Такое возможно только при $g_j = 0$ , а это означает, что вектор $f_j$ линейно зависит от векторов $\{f_1, . . . , f_{j-1}\}$ . Если в ходе процесса $r_{jj}$ оказывается равным нулю, то коэффициент $\alpha_j$ обнуляется и $j$ -й признак далее не учитывается, как будто его вообще не существовало. Если $r_{jj}$ не равен, но близок к нулю, может возникнуть проблема неустойчивости решения, поскольку на $r_{jj}$ приходится делить. На практике признак $f_j$ исключают, например, по такому условию: $r_{jj} < \delta\max_{i<j}r_{ii}$ , где $\delta$ имеет порядок $10^{-2}..10^{-5}$ .

Назовём вектор $\alpha_j = F^+_jy$ текущим вектором коэффициентов на $j$ -м шаге. Этот вектор имеет размерность $j$ . По окончании процесса $\alpha_n = \alpha^*$ .

Лемма 1.3. Пусть выполняются условия предыдущей леммы. Тогда на $j$ -м шаге процесса вектор $\alpha_j$ может быть вычислен по рекуррентной формуле

$\alpha_1 = t_{11}(y^T\tilde{g}_1)$ ,

$\alpha_j = \begin{pmatrix} \alpha_{j-1} + t_{j}(y^T\tilde{g}_j) \\ t_{jj}(y^T\tilde{g}_j) \end{pmatrix}$ , $j = 2, . . . , n$ .

Назовём величину $Q_j = \min_\alpha ||y - F_j\alpha||^2 = ||y - F_j\alpha_j||^2$ текущим значением функционала $Q$ на $j$ -м шаге. Оно равно кратчайшему расстоянию от $y$ до линейной оболочки столбцов $F_j$ . По окончании процесса $Q_n = Q(\alpha^*)$ . Следующая лемма показывает, что текущее значение $Q_j$ от шага к шагу только уменьшается.

Лемма 1.4. Значения $Q_j$ могут быть вычислены в ходе ортогонализации по рекуррентной формуле

$Q_0 = ||y||^2$ ;

$Q_j = Q_{j-1} - (y^T\tilde{g}_j)^2, \ j = 1, . . . , n$ .

Лемма 1.5. Текущий вектор невязок $\varepsilon_j = y - F_j\alpha_j$ на $j$ -м шаге процесса ортогонализации вычисляется по рекуррентной формуле

$\varepsilon_0 = y$ ;

$\varepsilon_j = \varepsilon_{j-1} - \tilde{g}_j(y^T\tilde{g}_j), \ j = 1, . . . , n$ .

Алгоритм 1.1. Ортогонализация Грама-Шмидта

1: инициализация: $g_j$ := $f_j, \ j$ := $1, . . . , n$ ;

2: для $j$ := $1, . . . , n$

3: $\ \$ для $i$ := $1, . . . , j - 1$

4: $\ \ \ \ \ r_{ij}$ := $\tilde{g}^T_i g_j$ ; (вычисление $i$ -й компоненты вектор-столбца $r_j \in R^{j-1}$ );

5: $\ \ \ \ \ g_j$ := $g_j - \tilde{g}_ir_{ij}$ ; (ортогонализация $g_j$ относительно $g_i$ );

6: $\ \ r_{jj}$ := $||g_j||$ ;

Модифицированная ортогонализация Грама-Шмидта

Если вместо $r_{ij}$ := $\tilde{g}^T_i f_j$ вычислять $r_{ij}$ := $\tilde{g}^T_i g_j$ , то формально результат не изменится, поскольку $\tilde{g}^T_i g_j = \tilde{g}^T_i(f_j - \sum_{s=0}^{i-1}\tilde{g}_s r_{sj}) = \tilde{g}^T_i f_j - \sum_{s=0}^{i-1}\tilde{g}^T_i \tilde{g}_s r_{sj} = \tilde{g}^T_i f_j$ .

Данная модификация повышает численную устойчивость алгоритма. Это объясняется тем, что вектор $g_j$ обладает минимальной нормой среди всех векторов вида $f_j - \sum_{s=0}^{i-1}\beta_s\tilde{g}_s$ , где $\beta_s$ - произвольные коэффициенты. Поэтому при скалярном умножении на $g_j$ ошибки накапливаются существенно медленнее.

Прежде чем переходить к следующей модификации, запишем основную часть алгоритма ортогонализации, вычисляющую $G_j$ и $R_j$ .

Изменим порядок ортогонализации столбцов. До сих пор мы ортогонализовали столбец $g_j$ относительно предыдущих столбцов $g_1, . . . , g_{j-1}$ . Но можно сделать и подругому - ортогонализовать все последующие столбцы $g_{j+1}, . . . , g_n$ относительно $g_j$ :

$g_i$ := $g_i - \tilde{g}_j(\tilde{g}^T_jg_i), \ i = j + 1, . . . , n$ .

Тогда в начале $j$ -го шага все столбцы $g_j , . . . , g_n$ по построению будут ортогональны всем столбцам $g_1, . . . , g_{j-1}$ . При этом подматрицы $G_j ,\ R_j ,\ T_j$ и вектор $\alpha_j$ останутся такими же, как и до модификации.

Описанная модификация обладает важным преимуществом. Теперь на каждом шаге можно выбрать столбец $g_m \in \{g_j , . . . , g_n\}$ , добавление которого наиболее выгодно. Чтобы не менять обозначений, будем полагать, что перед добавлением столбцы $g_j$ и $g_m$ переставляются местами (при реализации придётся сохранять соответствие между старой и новой нумерацией признаков, но мы не будем останавливаться на столь мелких технических деталях).

Возможны альтернативные критерии выбора добавляемого столбца:

1) столбец с максимальной нормой $||g_m|| \to \max_m$ , что соответствует выбору столбца $f_m$ , максимально некоррелированного с $g_1, . . . , g_{j-1}$ ; применение этого критерия решает проблему вырожденности $R_j$ (см. Замечание 1.1); здесь существенно, чтобы матрица $F$ была заранее стандартизована;

2) столбец, наиболее коррелированный с вектором ответов: $\frac{y^Tg_m}{||g_m||} \to \max_m$ ; его добавление ведёт к скорейшему убыванию функционала $Q$ ;

3) столбец, добавление которого ведёт к наименьшему увеличению нормы вектора коэффициентов $||\alpha_j||$ , что соответствует применению регуляризации;

4) столбец, после добавления которого значение функционала качества $Q$ на независимой контрольной выборке $X^k = \{x'_1, . . . , x'_k\}$ окажется минимальным, что соответствует применению скользящего контроля (hold-out CV).

Алгоритм 1.2. Решение линейной задачи наименьших квадратов путём ортогонализации Грама-Шмидта с последовательным добавлением признаков

Вход:

$F = (f_1, . . . , f_n)$ - матрица информации;

$y$ - вектор ответов;

$\delta Q$ - параметр критерия останова.

Выход:

$\alpha_j$ - вектор коэффициентов линейной комбинации;

$Q_j$ - минимальное значение функционала.

1: инициализация:

$\ \ Q_0$ := $||y||^2; \ g_j$ := $f_j; \ Z_j$ := $||g_j||^2; \ D_j$ := $y^Tg_j; \ \ j$ := $1, . . . , n;$

2: для $j$ := $1, . . . , n$

3: $\ \$ выбор $m \in {j, . . . , n}$ по критериям $Z_m \to \max_m$ и $(D^2_m /Z_m) \to \max_m;$

4: $\ \$ перестановка местами столбцов:

$\ \ \ \ \ \ g_j \rightleftharpoons g_m, \ f_j \rightleftharpoons f_m, \ r_j \rightleftharpoons r_m;$

5: $\ \ r_{jj}$ := $\sqrt{Z_m}$ ; нормировка: $\tilde{g}_j$ := $g_j/r_{jj};$

6: $\ \$ вычисление текущего значения функционала:

$\ \ \ \ \ \ d_j$ := $D_j/r_{jj}$ ; (эффективное вычисление $d_j$ := $y^T\tilde{g}_j$ );

$\ \ \ \ \ \ Q_j$ := $Q_{j-1} - d^2_j;$

7: $\ \$ обращение верхней треугольной матрицы $T_j = R^{-1}_j$ :

$\ \ \ \ \ \ t_{jj} = 1/r_{jj}; \ \ t_j = -t_{jj}T_{j-1}r_j$ (вектор-столбец $t_j$ длины $j - 1$ );

$\ \ \ \ \ \ T_j = \begin{pmatrix} T_{j-1} & & t_{j} \\ 0 & & t_{jj} \end{pmatrix};$

8: $\ \$ вычисление текущего вектора коэффициентов:

$\ \ \ \ \ \ \alpha_j = \begin{pmatrix} \alpha_{j-1} + t_{j}d_j \\ t_{jj}d_j \end{pmatrix};$

9: $\ \ \$ если $Q_j <\delta Q$ то

10: $\ \ \ \ \$ прекратить добавление столбцов; выход;

11: $\ \$ для $i$ := $j + 1, . . . , n$

12: $\ \ \ \ r_{ji}$ := $\tilde{g}^T_jg_i$ ; (компоненты вектор-столбца $r_i$ );

13: $\ \ \ \ g_i$ := $g_i - \tilde{g}_jr_{ji}$ ; (ортогонализация $g_i$ относительно $g_j$ );

14: $\ \ \ \ Z_i$ := $Z_i - r^2_{ji}$ ; (теперь $Z_i = ||g_i||^2$ );

15: $\ \ \ \ D_i$ := $D_i - d_jr_{ji}$ ; (теперь $D_i = y^Tg_i$ );

16: конец цикла по $j$ .

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

В Алгоритме 1.2 применяются первые два критерия.

Алгоритм состоит из основного цикла по $j$ , в котором поочерёдно добавляются столбцы. На шаге 3 принимается решение, какой из оставшихся столбцов $f_m$ добавить, затем он меняется местами со столбцом $f_j$ . Шаги 5–8 обновляют текущие значения функционала $Q_j$ , обратной матрицы $T_j$ и коэффициентов $\alpha_j$ . На шаге 9 проверяется условие останова: если достаточная точность аппроксимации уже достигнута, то добавлять оставшиеся столбцы не имеет смысла. Таким образом, Алгоритм 1.2 осуществляет отбор информативных признаков (features selection). Шаги 11–15 реализуют вложенный цикл, в котором все последующие столбцы ортогонализуются относительно только что добавленного столбца. Заодно обновляются значения квадратов норм столбцов $Z_j = ||g_j||^2$ и скалярных произведений $D_j = y^Tg_j$ , необходимые для эффективного выбора признака $f_m$ на шаге 3 в следующей итерации.