Модифицированная ортогонализация Грама-Шмидта
Материал из MachineLearning.
 (Новая: {{Задание|DValov|Константин Воронцов|6 января 2010}})  | 
				 (→Литература)  | 
			||
| (4 промежуточные версии не показаны) | |||
| Строка 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>\ \ </tex>для <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^2_{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>.  | ||
| + | |||
| + | ----  | ||
| + | Наконец, можно использовать совокупность критериев, выбирая тот столбец, добавление которого выгодно с нескольких точек зрения.  | ||
| + | |||
| + | В Алгоритме 1.2 применяются первые два критерия.  | ||
| + | |||
| + | Алгоритм состоит из основного цикла по <tex>j</tex>, в котором поочерёдно добавляются столбцы. На шаге 3 принимается решение, какой из оставшихся столбцов <tex>f_m</tex> добавить, затем он меняется местами со столбцом <tex>f_j</tex> . Шаги 5–8 обновляют текущие  | ||
| + | значения функционала <tex>Q_j</tex> , обратной матрицы <tex>T_j</tex> и коэффициентов <tex>\alpha_j</tex>. На шаге 9 проверяется условие останова: если достаточная точность аппроксимации уже достигнута, то добавлять оставшиеся столбцы не имеет смысла. Таким образом, Алгоритм 1.2 осуществляет отбор информативных признаков (features selection). Шаги 11–15 реализуют вложенный цикл, в котором все последующие столбцы ортогонализуются относительно только что добавленного столбца. Заодно обновляются значения квадратов норм столбцов <tex>Z_j = ||g_j||^2</tex> и скалярных произведений <tex>D_j = y^Tg_j</tex> , необходимые для эффективного выбора признака <tex>f_m</tex> на шаге 3 в следующей итерации.  | ||
| + | |||
| + | == Литература ==  | ||
| + | *''К.В. Воронцов'', [[Машинное обучение (курс лекций, К.В.Воронцов)|Машинное обучение (курс лекций)]]  | ||
| + | |||
| + | {{ЗаданиеВыполнено|DValov|Константин Воронцов|7 января 2010}}  | ||
| + | |||
| + | [[Категория:Линейная регрессия]]  | ||
| + | [[Категория:Методы отбора признаков]]  | ||
Текущая версия
Модифицированная ортогонализация Грама-Шмидта - известный алгоритм ортогонализации, который строит ортогональные векторы 
, линейная оболочка которых совпадает с линейной оболочкой 
.
Введение
Рассмотрим ортогональное разложение , где 
 - верхняя треугольная матрица размера 
 × 
, 
 - ортогональная 
 × 
 - матрица, 
. Для любой матрицы 
 существует бесконечно много разложений указанного вида. Имея одно из них, легко выразить псевдообратную матрицу 
 через 
 и 
:
.
Для вычисления псевдообратной  достаточно построить какое-нибудь ортогональное разложение матрицы 
, обратить верхнюю треугольную матрицу 
 и умножить её на 
. Этот метод во многом удобнее явного обращения матрицы.
Метод ортогонализации Грама-Шмидта
Для построения разложения  воспользуемся процессом ортогонализации Грама-Шмидта. Запишем матрицы 
 и 
 по столбцам:
;
.
Волной здесь и далее обозначается операция нормирования вектора:
.
Процесс ортогонализации Грама-Шмидта строит ортогональные векторы , линейная оболочка которых совпадает с линейной оболочкой 
:
;
;
· · ·
.                              
Легко проверить, что для всех  из 
, 
, векторы 
 и 
 ортогональны.
Вспомогательные утверждения
Лемма 1.1. На -м шаге процесса, 
, матрица 
 представима в виде ортогонального разложения 
 , где
 - ортонормированная матрица;
 - верхняя треугольная матрица, 
По окончании процесса ортогонализации получаем ортогональное разложение .
С вычислительной точки зрения процесс Грама-Шмидта удобен тем, что на каждом шаге матрицы  и 
 получаются путём дописывания справа по одному столбцу к матрицам 
 и 
 соответственно. При этом предыдущие столбцы не изменяются (если не считать изменением дописывание нулей снизу к матрице 
 - при разумной программной реализации эти нули всё равно не хранятся).
В следующей лемме утверждается, что обратная матрица  также является верхней треугольной и формируется путём дописывания столбцов справа.
Лемма 1.2. Пусть матрицы  невырождены и в блочной записи имеют вид
;
,  
,
где  - скаляр, 
 - вектор-столбец размера 
. Тогда матрицы 
 могут быть вычислены по рекуррентной формуле
;
,  
,
где  - скаляр, 
 - вектор-столбец размера 
.
Замечание 1.1. Обеспечить невырожденность матрицы  в процессе ортогонализации очень просто. Допустим, матрица 
 невырождена. Поскольку 
 - верхняя треугольная, вырожденность может возникнуть только в том случае, если 
. Такое возможно только при 
, а это означает, что вектор 
 линейно зависит от векторов 
. Если в ходе процесса 
 оказывается равным нулю, то коэффициент 
 обнуляется и 
-й признак далее не учитывается, как будто его вообще не существовало. Если 
 не равен, но близок к нулю, может возникнуть проблема неустойчивости решения, поскольку на 
 приходится делить. На практике признак 
 исключают, например, по такому условию: 
, где 
 имеет порядок 
.
Назовём вектор  текущим вектором коэффициентов на 
-м шаге.
Этот вектор имеет размерность 
. По окончании процесса 
.
Лемма 1.3. Пусть выполняются условия предыдущей леммы. Тогда на -м шаге
процесса вектор 
 может быть вычислен по рекуррентной формуле
,
, 
.
Назовём величину  текущим значением функционала 
 на 
-м шаге. Оно равно кратчайшему расстоянию от 
 до линейной оболочки столбцов 
. По окончании процесса 
. Следующая лемма показывает, что текущее значение 
 от шага к шагу только уменьшается.
Лемма 1.4. Значения  могут быть вычислены в ходе ортогонализации по рекуррентной формуле
;
.
Лемма 1.5. Текущий вектор невязок  на 
-м шаге процесса ортогонализации вычисляется по рекуррентной формуле
;
.
Алгоритм 1.1. Ортогонализация Грама-Шмидта
1: инициализация:  := 
 := 
;
2: для  := 
3:   для 
 := 
4:  := 
 ; (вычисление 
-й компоненты вектор-столбца 
);
5:  := 
 ; (ортогонализация 
 относительно 
);
6:  := 
;
Модифицированная ортогонализация Грама-Шмидта
Если вместо  := 
вычислять 
 := 
 , то формально результат не изменится, поскольку
 .
Данная модификация повышает численную устойчивость алгоритма. Это объясняется тем, что вектор  обладает минимальной нормой среди всех векторов вида 
, где 
 - произвольные коэффициенты. Поэтому при скалярном умножении на 
 ошибки накапливаются существенно медленнее.
Прежде чем переходить к следующей модификации, запишем основную часть
алгоритма ортогонализации, вычисляющую  и 
.
Изменим порядок ортогонализации столбцов. До сих пор мы ортогонализовали столбец  относительно предыдущих столбцов 
. Но можно сделать и подругому - ортогонализовать все последующие столбцы 
 относительно 
 :
 := 
.
Тогда в начале -го шага все столбцы 
 по построению будут ортогональны всем столбцам 
. При этом подматрицы 
 и вектор 
 останутся такими же, как и до модификации.
Описанная модификация обладает важным преимуществом. Теперь на каждом шаге можно выбрать столбец , добавление которого наиболее выгодно. Чтобы не менять обозначений, будем полагать, что перед добавлением столбцы 
 и 
 переставляются местами (при реализации придётся сохранять соответствие между старой и новой нумерацией признаков, но мы не будем останавливаться на столь мелких технических деталях).
Возможны альтернативные критерии выбора добавляемого столбца:
1) столбец с максимальной нормой , что соответствует выбору столбца 
, максимально некоррелированного с 
; применение этого критерия решает проблему вырожденности 
 (см. Замечание 1.1); здесь существенно, чтобы матрица 
 была заранее стандартизована;
2) столбец, наиболее коррелированный с вектором ответов: ; его добавление ведёт к скорейшему убыванию функционала 
;
3) столбец, добавление которого ведёт к наименьшему увеличению нормы вектора коэффициентов , что соответствует применению регуляризации;
4) столбец, после добавления которого значение функционала качества  на независимой контрольной выборке 
 окажется минимальным, что соответствует применению скользящего контроля (hold-out CV).
Алгоритм 1.2. Решение линейной задачи наименьших квадратов путём ортогонализации Грама-Шмидта с последовательным добавлением признаков
Вход:
 - матрица информации;
 - вектор ответов;
 - параметр критерия останова.
Выход:
 - вектор коэффициентов линейной комбинации;
 - минимальное значение функционала.
1: инициализация:
 := 
 := 
 := 
 := 
 := 
2: для := 
3:  выбор 
 по критериям 
 и 
4:  перестановка местами столбцов:
5: := 
; нормировка: 
:= 
6:  вычисление текущего значения функционала:
 := 
 ; (эффективное вычисление 
 := 
);
 := 
7: обращение верхней треугольной матрицы 
:
 (вектор-столбец 
 длины 
);
8: вычисление текущего вектора коэффициентов:
9: если 
 то
10: прекратить добавление столбцов; выход;
11: для 
 := 
12:   := 
; (компоненты вектор-столбца 
);
13:  := 
; (ортогонализация 
 относительно 
);
14:  := 
; (теперь 
);
15:  := 
; (теперь 
);
16: конец цикла по .
Наконец, можно использовать совокупность критериев, выбирая тот столбец, добавление которого выгодно с нескольких точек зрения.
В Алгоритме 1.2 применяются первые два критерия.
Алгоритм состоит из основного цикла по , в котором поочерёдно добавляются столбцы. На шаге 3 принимается решение, какой из оставшихся столбцов 
 добавить, затем он меняется местами со столбцом 
 . Шаги 5–8 обновляют текущие
значения функционала 
 , обратной матрицы 
 и коэффициентов 
. На шаге 9 проверяется условие останова: если достаточная точность аппроксимации уже достигнута, то добавлять оставшиеся столбцы не имеет смысла. Таким образом, Алгоритм 1.2 осуществляет отбор информативных признаков (features selection). Шаги 11–15 реализуют вложенный цикл, в котором все последующие столбцы ортогонализуются относительно только что добавленного столбца. Заодно обновляются значения квадратов норм столбцов 
 и скалярных произведений 
 , необходимые для эффективного выбора признака 
 на шаге 3 в следующей итерации.
Литература
- К.В. Воронцов, Машинное обучение (курс лекций)
 
|   |  Данная статья была создана в рамках учебного задания.
 
 См. также методические указания по использованию Ресурса MachineLearning.ru в учебном процессе.  | 

