QR-разложение: ортогональность как суперсила

**Линейная регрессия на 100 000 примерах.** Учебник говорит: w = (XᵀX)⁻¹Xᵀy. Numpy делает иначе. `np.linalg.lstsq` внутри вычисляет QR-разложение и никогда не трогает XᵀX. Причина - числовая. Возведение в квадрат числа обусловленности матрицы X превращает решаемую задачу в нестабильную. QR этого не делает. Именно поэтому за 60 лет QR-разложение стало стандартом - не за красоту, а за устойчивость к реальным данным.

**Что этот урок открывает на самом деле**: QR - это не просто ещё одно разложение. Это способ мышления: любую задачу лучше решать через ортогональные преобразования, потому что ортогональные матрицы не усиливают погрешности. Три места, где это используется прямо сейчас: sklearn's linear regression, scipy's eigensolver, и инициализация нейросетей.

Что такое QR-разложение

Любую матрицу A размера m x n (m >= n) можно представить как произведение двух матриц:

**Ключевое свойство Q**: умножение на ортогональную матрицу не меняет длины векторов и углы между ними. Это и есть причина численной устойчивости - Q не усиливает погрешности округления, она их просто поворачивает.

Грам-Шмидт: строим Q вручную

Классический алгоритм построения Q - процесс Грама-Шмидта: берём столбцы A один за другим и последовательно делаем их ортонормированными.

Классический Грам-Шмидт теряет ортогональность при плохо обусловленных матрицах из-за накопления ошибок округления. На практике numpy использует отражения Хаусхолдера - они численно устойчивее.

Хаусхолдер: как numpy делает это на самом деле

Промышленная реализация QR использует отражения Хаусхолдера. Идея: серия отражений H1, H2, ... обнуляет элементы под диагональю столбец за столбцом, превращая A в верхнетреугольную R.

ML-применение 1: почему lstsq лучше нормального уравнения

Нормальное уравнение для линейной регрессии: w = (XᵀX)⁻¹Xᵀy. Элегантно, но опасно. Если матрица X имеет число обусловленности kappa, то XᵀX имеет kappa^2. Плохо обусловленные данные (мультиколлинеарность) превращаются в катастрофически плохие коэффициенты.

sklearn's LinearRegression тоже использует lstsq через scipy, который вызывает LAPACK's dgelsd (SVD-based) или dgelsy (QR-based). Нормальное уравнение через явное (XᵀX)⁻¹ - учебный артефакт, в производственном коде его не применяют.

ML-применение 2: QR-итерация для собственных векторов

scipy.linalg.eig находит все собственные значения квадратной матрицы. Под капотом - QR-итерационный алгоритм: последовательные QR-разложения превращают матрицу в квазитреугольную форму Шура, где собственные значения читаются с диагонали.

На практике используют QR с shifts (сдвигами) - это ускоряет сходимость с кубической до почти квадратичной для симметричных матриц. Метод реализован в LAPACK и применяется в scipy, numpy и MATLAB как стандарт.

ML-применение 3: ортогональная инициализация весов

При инициализации весов глубокой нейросети случайными числами из нормального распределения возникает проблема: матрицы перемножаются через слои, и если сингулярные значения > 1, градиенты взрываются; если < 1, затухают. Ортогональная инициализация решает это: Q * Q^T = I, все сингулярные значения равны 1.

Ортогональная инициализация особенно важна для рекуррентных сетей (LSTM, GRU), где матрица весов применяется сотни раз за последовательность. Стандартная инициализация Xavier/He не гарантирует ortogonality - только подходящую дисперсию. QR-инициализация даёт оба свойства.

LU против QR: когда что выбирать

Где QR-разложение работает прямо сейчас Components: numpy.linalg.lstsq, scipy.linalg.eig / eigh, torch.nn.init.orthogonal_, sklearn LinearRegression

Практика: Метод наименьших квадратов

- Почему numpy.linalg.lstsq не использует нормальное уравнение (X^TX)^{-1}X^Ty? - В чём смысл QR-итерации для поиска собственных значений? Почему перестановка R*Q вместо Q*R сходится к нужному? - Зачем инициализировать веса нейросети ортогональной матрицей, если Xavier/He уже работает?

A = Q * R A - исходная матрица m x n Q - ортогональная матрица m x m: Q^T * Q = I R - верхнетреугольная матрица m x n Пример для матрицы 3 x 2: A = | 1 1 | Q = | 1/sqrt(2) 1/sqrt(6) | R = | sqrt(2) 1/sqrt(2) | | 1 0 | | 1/sqrt(2) -1/sqrt(6) | | 0 sqrt(3/2) | | 0 1 | | 0 2/sqrt(6) | Q^T * Q = I (легко проверить: столбцы Q - ортонормированные)

Дано: столбцы a1, a2, ..., an матрицы A Шаг 1: e1 = a1 / ||a1|| Шаг 2: u2 = a2 - (a2 . e1) * e1 // вычитаем проекцию на e1 e2 = u2 / ||u2|| Шаг k: uk = ak - sum_{i<k} (ak . ei) * ei ek = uk / ||uk|| Матрица Q = [e1 | e2 | ... | en] Матрица R: R[i][j] = ai . ej (скалярные произведения) Почему это работает: каждый новый вектор - это старый столбец минус-проекции на уже построенные ортонормированные. Остаток всегда перпендикулярен всем предыдущим.

Матрица отражения: H = I - 2 * v*v^T / (v^T * v) Где v выбирается так, что H*a = -sign(a1)*||a||*e1 (отражает вектор a в направление первого базисного вектора) Процесс: H1 * A обнуляет всё под [1,1] H2 * H1 * A обнуляет всё под [2,2] (в подматрице 2:n, 2:n) ... Hk * ... * H1 * A = R - верхнетреугольная Тогда Q = H1 * H2 * ... * Hk (каждое H ортогонально) Главный бонус: каждый H ортогональный => накопление ошибок пропорционально eps*n, не eps*n^2 как в Граме-Шмидте.

Задача: min ||Xw - y||^2 Нормальное уравнение (ОПАСНО при мультиколлинеарности): X^T X w = X^T y => kappa(X^T X) = kappa(X)^2 Через QR (БЕЗОПАСНО): X = QR Q^T Q R w = Q^T y R w = Q^T y => kappa(R) = kappa(X) R верхнетреугольная - решается back substitution за O(n^2). Число обусловленности НЕ возводится в квадрат.

Дано: матрица A0 Итерация k: Ak = Qk * Rk // QR-разложение текущей матрицы A(k+1) = Rk * Qk // меняем порядок множителей При k -> infinity: Ak -> верхнетреугольная матрица T (форма Шура) Собственные значения = диагональные элементы T Собственные векторы = столбцы произведения Q1 * Q2 * ... * Qk Почему это работает: похоже на степенной метод, но сразу для всех собственных векторов. Каждая итерация сжимает ненужные направления быстрее, чем усиливает нужные.