Математика / Линейная алгебра

QR-разложение отражениями Хаусхолдера: устойчивое зануление столбцов

Отражение Хаусхолдера умножением слева зануляет элементы под выбранным диагональным элементом. Последовательность ортогональных отражений превращает матрицу A в верхнетреугольную R, а произведение транспонированных отражений даёт Q.

Опубликовано: Обновлено:

Формула

$$H=I-2\frac{vv^T}{v^Tv},\qquad A=QR$$
Схема Отражение столбца (4;3)
Отражение вектора с сохранением нормыx=(4;3)Hx=(−5;0)0x₁x₂Длина обоих векторов: 5

Синий вектор (4;3) переходит в коричневый (−5;0); обе длины равны 5. Это первый шаг численного примера.

Обозначения

$A$
исходная матрица размера m×n, единицы исходных данных
$v$
вектор, задающий отражающую гиперплоскость, единицы столбца
$H$
ортогональная матрица отражения, безразмерная
$Q$
произведение отражений, QᵀQ=I, безразмерная
$R$
верхнетреугольный множитель, единицы A

Условия применения

Для каждого шага выбирают v так, чтобы vᵀv не было нулём; sign(0) принимают равным 1, а знак ведущего элемента выбирают, чтобы избежать вычитания близких чисел. При m≥n первые n столбцов можно занулить под диагональю. Для задачи наименьших квадратов матрица должна иметь полный столбцовый ранг. Ортогональность выполняется в точной арифметике, а в вычислениях проверяется нормой ||QᵀQ−I||.

Ограничения

QR не устраняет плохую обусловленность исходной задачи: он лишь не добавляет к ней лишнего неортогонального усиления. Формирование Q целиком дороже хранения отражателей. Если столбец уже нулевой, шаг пропускают; при переопределённой системе после получения R решают треугольную задачу по первым n строкам и проверяют остаток по всем m строкам.

Подробное объяснение

Отражение Хаусхолдера — это симметричное ортогональное преобразование. Для ненулевого вектора x выбирают v=x+sign(x₁)||x||e₁. Тогда Hx имеет только первую ненулевую координату: геометрически вектор отражается относительно гиперплоскости, перпендикулярной v, и его длина не меняется. Знак выбирают согласованно с x₁, чтобы сумма не была близка к нулю.

В матрице A отражение применяют слева к текущему подстолбцу. После первого шага все элементы ниже a₁₁ становятся нулём; после второго — ниже a₂₂ во втором столбце, и так далее. Уже обработанные строки не изменяются, потому что вектор v дополняют нулевыми координатами выше текущей позиции. Произведение H_n…H₂H₁A=R, а так как H_iᵀ=H_i и H_i²=I, то Q=H₁H₂…H_n.

Ортогональность даёт ||Qy||₂=||y||₂ и ||Hx||₂=||x||₂. Поэтому решение наименьших квадратов получают из Qᵀb: минимизация ||Ax−b||₂ равносильна минимизации ||Rx−Qᵀb||₂, после чего решают верхнетреугольную систему. Остаток считают в исходной матрице A, а не в округлённой R.

Для реального кода полезно хранить v и коэффициент τ=2/(vᵀv), применяя обновление подматрицы как A←A−v(τvᵀA). Это уменьшает память и число операций. Независимые проверки — восстановление QR и оценка ||QᵀQ−I||₂ или нормы Фробениуса.

Как пользоваться формулой

  1. Возьмите текущий подстолбец x и вычислите его евклидову норму без преждевременного округления.
  2. Сформируйте v=x+sign(x₁)||x||e₁ и коэффициент τ=2/(vᵀv); при x=0 шаг пропустите.
  3. Обновите текущую и все правые колонки по правилу A←A−v(τvᵀA), оставляя верхнюю часть неизменной.
  4. Повторите шаг для следующего диагонального элемента и соберите R из полученной верхнетреугольной части.
  5. Примените отражения к b в порядке H₁, H₂ и далее, чтобы получить Qᵀb; Q собирайте произведением H₁H₂….
  6. Проверьте ||QᵀQ−I|| и ||A−QR||, затем решите треугольную систему и вычислите исходную невязку.

Историческая справка

Отражения, переводящие вектор в координатное направление, были введены Альстоном Скоттом Хаусхолдером в статье 1958 года. В ней преобразование использовалось как вычислительный инструмент линейной алгебры и спектральных задач. В последующей вычислительной практике отражения стали стандартной альтернативой ортогонализации Грама—Шмидта: они дают компактное представление и лучше контролируют накопление ошибок. Современные реализации LAPACK применяют блоковые варианты, где несколько отражений объединяются для эффективной работы с кэш-памятью. Историческое имя относится к преобразованию, а алгоритм QR-разложения развивался благодаря работам многих авторов численного анализа. Первичная работа называется Unitary Triangularization of a Nonsymmetric Matrix и опубликована в Journal of the ACM, том 5, страницы 339–342. Современный QR использует ту же геометрию сохранения нормы, но организация блоковых операций относится к более поздним библиотечным реализациям.

Историческая линия формулы

Первичная атрибуция отражения принадлежит A. S. Householder (1958), однако современное QR-разложение и его устойчивые блочные реализации — результат дальнейшего развития вычислительной линейной алгебры. Поэтому имя автора относится к преобразованию, а не означает, что вся современная библиотечная реализация создана одним исследователем.

Пример

Разложим A=[[4,0],[3,−5]]. Первый столбец x=(4,3)ᵀ имеет норму 5. Берём v=(9,3)ᵀ, vᵀv=90 и H₁=I−(1/45)vvᵀ=[[-4/5,−3/5],[-3/5,4/5]]. Тогда H₁A=[[-5,3],[0,−4]]. Нижний элемент первого столбца занулён, а норма столбца сохранилась: √(4²+3²)=5. Здесь R=H₁A, Q=H₁, потому что один шаг достаточен. Прямая проверка QR=H₁(H₁A)=A использует H₁²=I. Для правой части b=(7,−1)ᵀ сначала Qᵀb=H₁b=(−5,−5)ᵀ. Обратная подстановка в R: −4x₂=−5, x₂=5/4; затем −5x₁+3x₂=−5, откуда x₁=7/4. Подстановка в исходную матрицу даёт Ax=(7,−1)ᵀ, то есть невязка равна нулю. Этот квадратный пример показывает полный цикл QR и проверку результата в исходной системе.

Частая ошибка

Нельзя брать v=x−||x||e₁ автоматически: при x₁>0 это вызывает катастрофическое вычитание. Нельзя считать Q ортогональной только потому, что R верхнетреугольная; проверьте QᵀQ. Нельзя решать задачу по R, забыв преобразовать правую часть b теми же отражениями. Нельзя проверять остаток только по округлённому R: восстановите Ax исходной матрицей. При нулевом текущем столбце деление на vᵀv запрещено, шаг должен быть пропущен.

Практика

Задачи с решением

Одно отражение

Условие. Для x=(3,4)ᵀ постройте v=x+5e₁ и найдите Hx.

Решение. v=(8,4)ᵀ, vᵀv=80. H=I−vvᵀ/40=[[-3/5,−4/5],[−4/5,3/5]]. Поэтому Hx=(-3/5·3−4/5·4, −4/5·3+3/5·4)=(-5,0).

Ответ. Hx=(-5,0)ᵀ.

Проверка нормы

Условие. Покажите для x=(3,4)ᵀ, что отражение не меняет евклидову норму.

Решение. До отражения ||x||₂=√(9+16)=5. После отражения Hx=(-5,0), поэтому ||Hx||₂=√25=5. Равенство также следует из HᵀH=I.

Ответ. Обе нормы равны 5.

Дополнительные источники

  • https://web.stanford.edu/class/cme324/classics/householder.pdf
  • https://www.netlib.org/lapack/lug/node69.html

Связанные формулы

Математика

Матричная форма системы линейных уравнений

$Ax=b$

Матричная форма Ax = b записывает систему линейных уравнений как произведение матрицы коэффициентов на столбец неизвестных. Такая запись позволяет решать систему не как набор отдельных строк, а как единый линейный объект.

Математика

Прямой ход метода Гаусса

$R_i\leftarrow R_i-\frac{a_{ik}}{a_{kk}}R_k$

Прямой ход метода Гаусса зануляет коэффициенты под ведущими элементами. В результате система приводится к ступенчатому виду, из которого решение находят обратной подстановкой.

Математика

Обратная подстановка в методе Гаусса

$x_i=\frac{b'_i-\sum_{j=i+1}^{n}u_{ij}x_j}{u_{ii}}$

Обратная подстановка находит неизвестные после прямого хода метода Гаусса. Она идет снизу вверх по ступенчатой системе: сначала последняя ведущая переменная, затем предыдущие.