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

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

Перейти к: навигация, поиск


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

Содержание

Введение

Рассмотрим ортогональное разложение 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 в следующей итерации.

Литература


Данная статья была создана в рамках учебного задания.
Студент: Участник:DValov
Преподаватель: Участник:Константин Воронцов
Срок: 7 января 2010


В настоящее время задание завершено и проверено. Данная страница может свободно правиться другими участниками проекта MachineLearning.ru.

См. также методические указания по использованию Ресурса MachineLearning.ru в учебном процессе.

Личные инструменты