Математика / Линейная алгебра
QR-разложение отражениями Хаусхолдера: устойчивое зануление столбцов
Отражение Хаусхолдера умножением слева зануляет элементы под выбранным диагональным элементом. Последовательность ортогональных отражений превращает матрицу A в верхнетреугольную R, а произведение транспонированных отражений даёт Q.
Формула
Синий вектор (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||₂ или нормы Фробениуса.
Как пользоваться формулой
- Возьмите текущий подстолбец x и вычислите его евклидову норму без преждевременного округления.
- Сформируйте v=x+sign(x₁)||x||e₁ и коэффициент τ=2/(vᵀv); при x=0 шаг пропустите.
- Обновите текущую и все правые колонки по правилу A←A−v(τvᵀA), оставляя верхнюю часть неизменной.
- Повторите шаг для следующего диагонального элемента и соберите R из полученной верхнетреугольной части.
- Примените отражения к b в порядке H₁, H₂ и далее, чтобы получить Qᵀb; Q собирайте произведением H₁H₂….
- Проверьте ||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 записывает систему линейных уравнений как произведение матрицы коэффициентов на столбец неизвестных. Такая запись позволяет решать систему не как набор отдельных строк, а как единый линейный объект.
Математика
Прямой ход метода Гаусса
Прямой ход метода Гаусса зануляет коэффициенты под ведущими элементами. В результате система приводится к ступенчатому виду, из которого решение находят обратной подстановкой.
Математика
Обратная подстановка в методе Гаусса
Обратная подстановка находит неизвестные после прямого хода метода Гаусса. Она идет снизу вверх по ступенчатой системе: сначала последняя ведущая переменная, затем предыдущие.