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

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

(Различия между версиями)
Перейти к: навигация, поиск
(Новая: {{Задание|DValov|Константин Воронцов|6 января 2010}})
Строка 1: Строка 1:
 +
[[Категория:Непроверенные учебные задания]]
 +
 +
'''Модифицированная ортогонализация Грама-Шмидта''' - известный алгоритм ортогонализации, который строит ортогональные векторы
 +
<tex>g_1, . . . , g_n</tex>, линейная оболочка которых совпадает с линейной оболочкой <tex>f_1, . . . , f_n</tex>.
 +
 +
== Введение ==
 +
 +
Рассмотрим ортогональное разложение <tex>F = GR</tex>, где <tex>R</tex> - верхняя треугольная матрица размера <tex>n</tex> × <tex>n</tex>, <tex>G</tex> - ортогональная <tex>l</tex> × <tex>n</tex> - матрица, <tex>G^TG = I_n</tex>. Для любой матрицы <tex>F</tex> существует бесконечно много разложений указанного вида. Имея одно из них, легко выразить псевдообратную матрицу <tex>F^+</tex> через <tex>G</tex> и <tex>R</tex>:
 +
 +
<tex>F^+ = (R^TG^TGR)^{-1}R^TG^T = R^{-1}R^{-1T}R^TG^T = R^{-1}G^T</tex>.
 +
 +
Для вычисления псевдообратной <tex>F^+</tex> достаточно построить какое-нибудь ортогональное разложение матрицы <tex>F</tex>, обратить верхнюю треугольную матрицу <tex>R</tex> и умножить её на <tex>G^T</tex>. Этот метод во многом удобнее явного обращения матрицы.
 +
 +
== Метод ортогонализации Грама-Шмидта ==
 +
 +
Для построения разложения <tex>F = GR</tex> воспользуемся процессом ортогонализации Грама-Шмидта. Запишем матрицы <tex>F</tex> и <tex>G</tex> по столбцам:
 +
 +
<tex>F = (f_1, . . . , f_n)</tex>;
 +
 +
<tex>G = (\tilde{g}_1, . . . , \tilde{g}_n)</tex>.
 +
 +
Волной здесь и далее обозначается операция нормирования вектора:
 +
 +
<tex>\tilde{v} = v/||v||</tex>.
 +
 +
Процесс ортогонализации Грама-Шмидта строит ортогональные векторы <tex>g_1, . . . , g_n</tex>, линейная оболочка которых совпадает с линейной оболочкой <tex>f_1, . . . , f_n</tex>:
 +
 +
<tex>g_1 = f_1</tex>;
 +
 +
<tex>g_2 = f_2 - \tilde{g}_1\tilde{g}^T_1f_2</tex>;
 +
 +
· · ·
 +
 +
<tex>g_j = f_j - \tilde{g}_1\tilde{g}^T_1f_j - ... - \tilde{g}_{j-1}\tilde{g}^T_{j-1}f_j</tex>.
 +
 +
Легко проверить, что для всех <tex>k, j</tex> из <tex>\{1, . . . , n\}</tex>, <tex>k \ne j</tex>, векторы <tex>\tilde{g}_k</tex> и <tex>\tilde{g}_j</tex> ортогональны.
 +
 +
== Вспомогательные утверждения ==
 +
 +
'''Лемма 1.1.''' На <tex>j</tex>-м шаге процесса, <tex>j = 1, . . . , n</tex>, матрица <tex>F_j = (f_1, . . . , f_j)</tex> представима в виде ортогонального разложения <tex>F_j = G_jR_j</tex> , где
 +
 +
<tex>G_j = (\tilde{g}_1, . . . , \tilde{g}_j)</tex> - ортонормированная матрица;
 +
 +
<tex>R_j = \begin{pmatrix}
 +
r_{11} & \cdots & r_{1j} \\
 +
& \ddots & \vdots \\
 +
0 & & r_{jj}
 +
\end{pmatrix}</tex> - верхняя треугольная матрица, <tex>r_{ij} = \left\{
 +
\begin{eqnarray}
 +
\tilde{g}^T_if_j, & i < j;\\
 +
||g_j||, & i=j.
 +
\end{eqnarray}</tex>
 +
 +
 +
По окончании процесса ортогонализации получаем ортогональное разложение <tex>F = G_nR_n</tex>.
 +
 +
С вычислительной точки зрения процесс Грама-Шмидта удобен тем, что на каждом шаге матрицы <tex>G_j</tex> и <tex>R_j</tex> получаются путём дописывания справа по одному столбцу к матрицам <tex>G_{j-1}</tex> и <tex>R_{j-1}</tex> соответственно. При этом предыдущие столбцы не изменяются (если не считать изменением дописывание нулей снизу к матрице <tex>R_j</tex> - при разумной программной реализации эти нули всё равно не хранятся).
 +
 +
В следующей лемме утверждается, что обратная матрица <tex>T_j = R^{-1}_j</tex> также является верхней треугольной и формируется путём дописывания столбцов справа.
 +
 +
'''Лемма 1.2.''' Пусть матрицы <tex>R_j</tex> невырождены и в блочной записи имеют вид
 +
 +
<tex>R_1 = (r_{11})</tex>;
 +
 +
<tex>R_j = \begin{pmatrix}
 +
R_{j-1} & & r_{j} \\
 +
0 & & r_{jj}
 +
\end{pmatrix}</tex>, <tex>j = 2, . . . , n</tex>,
 +
 +
где <tex>r_{jj}</tex> - скаляр, <tex>r_j</tex> - вектор-столбец размера <tex>(j-1)</tex>. Тогда матрицы <tex>T_j = R^{-1}_j</tex> могут быть вычислены по рекуррентной формуле
 +
 +
<tex>T_1 = (t_{11})</tex>;
 +
 +
<tex>T_j = \begin{pmatrix}
 +
T_{j-1} & & t_{j} \\
 +
0 & & t_{jj}
 +
\end{pmatrix}</tex>, <tex>j = 2, . . . , n</tex>,
 +
 +
где <tex>t_{jj} = 1/r_{jj}</tex> - скаляр, <tex>t_j = -t_{jj}T_{j-1}r_j</tex> - вектор-столбец размера <tex>(j - 1)</tex>.
 +
 +
'''Замечание 1.1.''' Обеспечить невырожденность матрицы <tex>R_j</tex> в процессе ортогонализации очень просто. Допустим, матрица <tex>R_{j-1}</tex> невырождена. Поскольку <tex>R_j</tex> - верхняя треугольная, вырожденность может возникнуть только в том случае, если <tex>r_{jj} = 0</tex>. Такое возможно только при <tex>g_j = 0</tex>, а это означает, что вектор <tex>f_j</tex> линейно зависит от векторов <tex>\{f_1, . . . , f_{j-1}\}</tex>. Если в ходе процесса <tex>r_{jj}</tex> оказывается равным нулю, то коэффициент <tex>\alpha_j</tex> обнуляется и <tex>j</tex>-й признак далее не учитывается, как будто его вообще не существовало. Если <tex>r_{jj}</tex> не равен, но близок к нулю, может возникнуть проблема неустойчивости решения, поскольку на <tex>r_{jj}</tex> приходится делить. На практике признак <tex>f_j</tex> исключают, например, по такому условию: <tex>r_{jj} < \delta\max_{i<j}r_{ii}</tex>, где <tex>\delta</tex> имеет порядок <tex>10^{-2}..10^{-5}</tex>.
 +
 +
Назовём вектор <tex>\alpha_j = F^+_jy</tex> текущим вектором коэффициентов на <tex>j</tex>-м шаге.
 +
Этот вектор имеет размерность <tex>j</tex>. По окончании процесса <tex>\alpha_n = \alpha^*</tex>.
 +
 +
'''Лемма 1.3.''' Пусть выполняются условия предыдущей леммы. Тогда на <tex>j</tex>-м шаге
 +
процесса вектор <tex>\alpha_j</tex> может быть вычислен по рекуррентной формуле
 +
 +
<tex>\alpha_1 = t_{11}(y^T\tilde{g}_1)</tex>,
 +
 +
<tex>\alpha_j = \begin{pmatrix}
 +
\alpha_{j-1} + t_{j}(y^T\tilde{g}_j) \\
 +
t_{jj}(y^T\tilde{g}_j)
 +
\end{pmatrix}</tex>, <tex> j = 2, . . . , n</tex>.
 +
 +
 +
Назовём величину <tex>Q_j = \min_\alpha ||y - F_j\alpha||^2 = ||y - F_j\alpha_j||^2</tex> текущим значением функционала <tex>Q</tex> на <tex>j</tex>-м шаге. Оно равно кратчайшему расстоянию от <tex>y</tex> до линейной оболочки столбцов <tex>F_j</tex>. По окончании процесса <tex>Q_n = Q(\alpha^*)</tex>. Следующая лемма показывает, что текущее значение <tex>Q_j</tex> от шага к шагу только уменьшается.
 +
 +
'''Лемма 1.4.''' Значения <tex>Q_j</tex> могут быть вычислены в ходе ортогонализации по рекуррентной формуле
 +
 +
<tex>Q_0 = ||y||^2</tex>;
 +
 +
<tex>Q_j = Q_{j-1} - (y^T\tilde{g}_j)^2, \ j = 1, . . . , n</tex>.
 +
 +
'''Лемма 1.5.''' Текущий вектор невязок <tex>\varepsilon_j = y - F_j\alpha_j</tex> на <tex>j</tex>-м шаге процесса ортогонализации вычисляется по рекуррентной формуле
 +
 +
<tex>\varepsilon_0 = y</tex>;
 +
 +
<tex>\varepsilon_j = \varepsilon_{j-1} - \tilde{g}_j(y^T\tilde{g}_j), \ j = 1, . . . , n</tex>.
 +
 +
== Алгоритм 1.1. Ортогонализация Грама-Шмидта ==
 +
 +
1: инициализация: <tex>g_j</tex> := <tex>f_j, \ j</tex> := <tex>1, . . . , n</tex>;
 +
 +
2: для <tex>j</tex> := <tex>1, . . . , n</tex>
 +
 +
3: <tex> \ \ </tex> для <tex>i</tex> := <tex>1, . . . , j - 1</tex>
 +
 +
4: <tex> \ \ \ \ \ r_{ij}</tex> := <tex>\tilde{g}^T_i g_j</tex> ; (вычисление <tex>i</tex>-й компоненты вектор-столбца <tex>r_j \in R^{j-1}</tex>);
 +
 +
5: <tex> \ \ \ \ \ g_j</tex> := <tex>g_j - \tilde{g}_ir_{ij}</tex> ; (ортогонализация <tex>g_j</tex> относительно <tex>g_i</tex>);
 +
 +
6: <tex> \ \ r_{jj}</tex> := <tex>||g_j||</tex>;
 +
 +
 +
== Модифицированная ортогонализация Грама-Шмидта ==
 +
 +
Если вместо <tex>r_{ij}</tex> := <tex>\tilde{g}^T_i f_j</tex>
 +
вычислять <tex>r_{ij}</tex> := <tex>\tilde{g}^T_i g_j</tex> , то формально результат не изменится, поскольку
 +
<tex>\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</tex> .
 +
 +
Данная модификация повышает численную устойчивость алгоритма. Это объясняется тем, что вектор <tex>g_j</tex> обладает минимальной нормой среди всех векторов вида <tex>f_j - \sum_{s=0}^{i-1}\beta_s\tilde{g}_s</tex>, где <tex>\beta_s</tex> - произвольные коэффициенты. Поэтому при скалярном умножении на <tex>g_j</tex> ошибки накапливаются существенно медленнее.
 +
 +
Прежде чем переходить к следующей модификации, запишем основную часть
 +
алгоритма ортогонализации, вычисляющую <tex>G_j</tex> и <tex>R_j</tex>.
 +
 +
Изменим порядок ортогонализации столбцов. До сих пор мы ортогонализовали столбец <tex>g_j</tex> относительно предыдущих столбцов <tex>g_1, . . . , g_{j−1}</tex>. Но можно сделать и подругому - ортогонализовать все последующие столбцы <tex>g_{j+1}, . . . , g_n</tex> относительно <tex>g_j</tex> :
 +
 +
<tex>g_i</tex> := <tex>g_i - \tilde{g}_j(\tilde{g}^T_jg_i), \ i = j + 1, . . . , n</tex>.
 +
 +
Тогда в начале <tex>j</tex>-го шага все столбцы <tex>g_j , . . . , g_n</tex> по построению будут ортогональны всем столбцам <tex>g_1, . . . , g_{j-1}</tex>. При этом подматрицы <tex>G_j ,\ R_j ,\ T_j</tex> и вектор <tex>\alpha_j</tex> останутся такими же, как и до модификации.
 +
 +
Описанная модификация обладает важным преимуществом. Теперь на каждом шаге можно выбрать столбец <tex>g_m \in {g_j , . . . , g_n}</tex>, добавление которого наиболее выгодно. Чтобы не менять обозначений, будем полагать, что перед добавлением столбцы <tex>g_j</tex> и <tex>g_m</tex> переставляются местами (при реализации придётся сохранять соответствие между старой и новой нумерацией признаков, но мы не будем останавливаться на столь мелких технических деталях).
 +
 +
Возможны альтернативные критерии выбора добавляемого столбца:
 +
 +
1) столбец с максимальной нормой <tex>||g_m|| \to \max_m</tex>, что соответствует выбору столбца <tex>f_m</tex>, максимально некоррелированного с <tex>g_1, . . . , g_{j-1}</tex>; применение этого критерия решает проблему вырожденности <tex>R_j</tex> (см. Замечание 1.1); здесь существенно, чтобы матрица <tex>F</tex> была заранее стандартизована;
 +
 +
2) столбец, наиболее коррелированный с вектором ответов: <tex>\frac{y^Tg_m}{||g_m||} \to \max_m</tex>; его добавление ведёт к скорейшему убыванию функционала <tex>Q</tex>;
 +
 +
3) столбец, добавление которого ведёт к наименьшему увеличению нормы век-
 +
тора коэффициентов <tex>||\alpha_j||</tex>, что соответствует применению регуляризации;
 +
 +
4) столбец, после добавления которого значение функционала качества <tex>Q</tex> на независимой контрольной выборке <tex>X^k = \{x'_1, . . . , x'_k\}</tex> окажется минимальным, что соответствует применению скользящего контроля (hold-out CV).
 +
 +
== Алгоритм 1.2. Решение линейной задачи наименьших квадратов путём ортогонализации Грама-Шмидта с последовательным добавлением признаков ==
 +
 +
Вход:
 +
 +
<tex>F = (f_1, . . . , f_n)</tex> - матрица информации;
 +
 +
<tex>y</tex> - вектор ответов;
 +
 +
<tex>\delta Q</tex> - параметр критерия останова.
 +
 +
Выход:
 +
 +
<tex>\alpha_j</tex> - вектор коэффициентов линейной комбинации;
 +
 +
<tex>Q_j</tex> - минимальное значение функционала.
 +
 +
1: инициализация:
 +
 +
<tex>\ \ Q_0</tex> := <tex>||y||^2; \ g_j</tex> := <tex>f_j; \ Z_j</tex> := <tex>||g_j||^2; \ D_j</tex> := <tex>y^Tg_j; \ \ j</tex> := <tex>1, . . . , n;</tex>
 +
 +
2: для <tex>j</tex>:= <tex>1, . . . , n</tex>
 +
 +
3: <tex>\ \ </tex> выбор <tex>m \in {j, . . . , n}</tex> по критериям <tex>Z_m \to \max_m </tex> и <tex>(D^2_m
 +
/Z_m) \to \max_m;</tex>
 +
 +
4: <tex>\ \ </tex> перестановка местами столбцов:
 +
 +
<tex>\ \ \ \ \ \ g_j \rightleftharpoons g_m, \ f_j \rightleftharpoons f_m, \ r_j \rightleftharpoons r_m;</tex>
 +
 +
5: <tex>\ \ r_{jj}</tex>:= <tex>\sqrt{Z_m}</tex>; нормировка: <tex>\tilde{g}_j</tex>:= <tex>g_j/r_{jj};</tex>
 +
 +
6: <tex>\ \ </tex> вычисление текущего значения функционала:
 +
 +
<tex>\ \ \ \ \ \ d_j</tex> := <tex>D_j/r_{jj}</tex> ; (эффективное вычисление <tex>d_j</tex> := <tex>y^T\tilde{g}_j</tex>);
 +
 +
<tex>\ \ \ \ \ \ Q_j</tex> := <tex>Q_{j-1} - d^2_j;</tex>
 +
 +
7: <tex>\ \ </tex>обращение верхней треугольной матрицы <tex>T_j = R^{-1}_j </tex>:
 +
 +
<tex>\ \ \ \ \ \ t_{jj} = 1/r_{jj}; \ \ t_j = -t_{jj}T_{j-1}r_j</tex> (вектор-столбец <tex>t_j</tex> длины <tex>j - 1</tex>);
 +
 +
<tex>\ \ \ \ \ \ T_j = \begin{pmatrix}
 +
T_{j-1} & & t_{j} \\
 +
0 & & t_{jj}
 +
\end{pmatrix};</tex>
 +
 +
8: <tex>\ \ </tex>вычисление текущего вектора коэффициентов:
 +
 +
<tex>\ \ \ \ \ \ \alpha_j = \begin{pmatrix}
 +
\alpha_{j-1} + t_{j}d_j \\
 +
t_{jj}d_j
 +
\end{pmatrix};</tex>
 +
 +
9: <tex>\ \ \ </tex>если <tex>Q_j <\delta Q</tex> то
 +
 +
10: <tex>\ \ \ \ \ </tex>прекратить добавление столбцов; выход;
 +
 +
11: для <tex>i</tex> := <tex>j + 1, . . . , n</tex>
 +
 +
12: <tex>\ \ \ \ r_{ji}</tex> := <tex>\tilde{g}^T_jg_i</tex>; (компоненты вектор-столбца <tex>r_i</tex>);
 +
 +
13: <tex>\ \ \ \ g_i</tex> := <tex>g_i - \tilde{g}_jr_{ji}</tex>; (ортогонализация <tex>g_i</tex> относительно <tex>g_j</tex>);
 +
 +
14: <tex>\ \ \ \ Z_i</tex> := <tex>Z_i - r_{ji}</tex>; (теперь <tex>Z_i = ||g_i||^2</tex>);
 +
 +
15: <tex>\ \ \ \ D_i</tex> := <tex>D_i - d_jr_{ji}</tex>; (теперь <tex>D_i = y^Tg_i</tex>);
 +
 +
16: конец цикла по <tex>j</tex>.
 +
 +
 +
 +
 +
 +
 +
 +
 +
 +
{{Задание|DValov|Константин Воронцов|6 января 2010}}
{{Задание|DValov|Константин Воронцов|6 января 2010}}

Версия 17:40, 5 января 2010


Модифицированная ортогонализация Грама-Шмидта - известный алгоритм ортогонализации, который строит ортогональные векторы 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_{ji}; (теперь Z_i = ||g_i||^2);

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

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






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

До указанного срока статья не должна редактироваться другими участниками проекта MachineLearning.ru. По его окончании любой участник вправе исправить данную статью по своему усмотрению и удалить данное предупреждение, выводимое с помощью шаблона {{Задание}}.

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