1611688847-5b7354cc83380cb6c671f7c9dd5f83f8 (826648), страница 14
Текст из файла (страница 14)
. . , Hn−1 , состоящая из матрицотражения и, возможно, единичных матриц, таких что матрицаHn−1 Hn−2 · · · H2 H1 A = Rявляется правой треугольной матрицей.ПредложениеДля любой квадратной матрицы A существует конечнаяпоследовательность H1 , H2 , . . . , Hn−1 , состоящая из матрицотражения и, возможно, единичных матриц, таких что матрицаHn−1 Hn−2 · · · H2 H1 A = Rявляется правой треугольной матрицей.Раздельное упоминание матриц отражения и единичных матрицвызвано здесь тем, что единичная матрица не является матрицейотражения.Удобно применять систему обозначений матрично-векторныхобъектов, укоренившуюся в языках программирования Fortran,Matlab, Scilab, и им подобных.A(p : q, r : s) обозначает сечение массива A — массив с тем жеколичеством измерений и имеющий элементы, которые стоятна пересечении строк с номерами с p по qи столбцов с номерами с r по s.То есть, запись A(p : q, r : s) указывает в индексах матрицы Aне отдельные значения, а целые диапазоны изменения индексовэлементов, из которых образуется новая матрица, как подматрицаисходной.Доказательство.Используя результат Предложения о действии матрицы отражения,возьмём в качестве H1 матрицу отражения, которая переводит 1-йстолбец A в вектор, коллинеарный (1, 0, .
. . , 0)⊤ , если хотя бы одиниз элементов a21 , a31 , . . . , an1 не равен нулю.Иначе полагаем H1 = I.Затем переходим к следующему шагу.В результате выполнения первого шагаматрица СЛАУ приводится, как и в методе Гаусса, к видуà = ×××···×0××0...×...×...···......0××···×......×.где «×» обозначает элементы, которые не равны нулю.В результате выполнения первого шагаматрица СЛАУ приводится, как и в методе Гаусса, к видуà = ×××···×0××0...×...×...···......0××···×......×.где «×» обозначает элементы, которые не равны нулю.Проделаем теперь то же самое с матрицей Ã(2 : n, 2 : n), обнулив унеё подходящим отражением поддиагональные элементы первогостолбца, который является вторым во всей большой матрице Ã.И так далее до (n − 1)-го столбца.Для формального описания алгоритма определим теперь матрицуHj = Hj (u),j = 2, 3, .
. . , n − 1,как n × n-матрицу отражения, порождаемую вектором Хаусхолдераu ∈ Rn , которыйимеет нулевыми первые j − 1 компонент,подобран так, чтобы Hj (u) аннулировала в матрицеà = Hj−1 · · · H2 H1 Aподдиагональные элементы j-го столбца, если среди нихсуществуют ненулевые.Иначе, если в преобразуемой матрице à = (ãij ) все элементы ãj+1,j ,ãj+2,j , . . . , ãnj — нулевые, то полагаем Hj = I — единичнойn × n-матрице.Можно положить в блочной формеHj =I00H̃j!,гдев верхнем левом углу стоит единичная (j − 1) × (j − 1)-матрица,H̃j — матрица размера (n − j + 1) × (n − j + 1), котораяпереводит вектор Ã(j : n, j) в (n − j + 1)-вектор, коллинеарный(1, 0, .
. . , 0)⊤ , т. е. обнуляет поддиагональные элементы j-гостолбца в Ã.Можно положить в блочной формеHj =I00H̃j!,гдев верхнем левом углу стоит единичная (j − 1) × (j − 1)-матрица,H̃j — матрица размера (n − j + 1) × (n − j + 1), котораяпереводит вектор Ã(j : n, j) в (n − j + 1)-вектор, коллинеарный(1, 0, . . . , 0)⊤ , т. е. обнуляет поддиагональные элементы j-гостолбца в Ã.Если хотя бы один из элементов ãj+1,j , ãj+2,j , . . .
, ãnj не равеннулю, то H̃j — матрица отражения, способ построения которойописывается в Предложении о действии матриц отражения.Иначе, если (ãj+1,j , ãj+2,j , . . . , ãnj )⊤ = 0, то H̃j единичная(n − j + 1) × (n − j + 1)-матрица.Из представленияHn−1 Hn−2 · · · H2 H1 A = Rвытекает равенствоA = QRс ортогональной матрицейQ = Hn−1 Hn−2 · · · H2 H1−1.Таким образом, мы получаем QR-разложение матрицы A.Из представленияHn−1 Hn−2 · · · H2 H1 A = Rвытекает равенствоA = QRс ортогональной матрицейQ = Hn−1 Hn−2 · · · H2 H1−1.Таким образом, мы получаем QR-разложение матрицы A.Теорема доказана.QR-разложение матрицы с помощью отражений ХаусхолдераDO FOR j = 1 TO n − 1I˜ ← единичная матрица размера (n − j + 1);IF вектор A((j +1) : n, j) ненулевой THENвычислить вектор Хаусхолдера u ∈ Rn−j+1 ,порождающий отражение, которое переводитвектор A(j : n, j) в вектор, коллинеарныйвектору (1, 0, .
. . , 0)⊤ ;H̃ ← I˜ − 2uu⊤ ;ELSEH̃ ← I˜ ;END IFA(j : n, j : n) ← H̃ A(j : n, j : n) ;END DOМетод отражений ХаусхолдераЕсли известно QR-разложение матрицы A, то система уравненийAx = bравносильна системе уравнений(Qy = b,Rx = y,составные части которой решаются просто.Метод отражений ХаусхолдераЕсли известно QR-разложение матрицы A, то система уравненийAx = bравносильна системе уравнений(Qy = b,Rx = y,составные части которой решаются просто.При практической реализации удобнее дополнить алгоритминструкциями, которые задают преобразования правой части.Тогда результатом работы алгоритма сразу будет праваятреугольная система Rx = y, которую можно решать с помощьюобратной подстановки.Согласно Предложению о действии матриц отражения, вычислениевектора Хаусхолдера u требует нахождения вектора ũ, в которомимеется неоднозначность выбора знака второго слагаемого.При вычислениях на цифровых ЭВМ в стандартной арифметике сплавающей точкой имеет смысл брать A(j : n, j) + kA(j : n, j)k2 e, если ajj ≥ 0,ũ = A(j : n, j) − kA(j : n, j)k e, если a ≤ 0,2jjгде e = (1, 0, .
. . , 0)⊤ .Согласно Предложению о действии матриц отражения, вычислениевектора Хаусхолдера u требует нахождения вектора ũ, в которомимеется неоднозначность выбора знака второго слагаемого.При вычислениях на цифровых ЭВМ в стандартной арифметике сплавающей точкой имеет смысл брать A(j : n, j) + kA(j : n, j)k2 e, если ajj ≥ 0,ũ = A(j : n, j) − kA(j : n, j)k e, если a ≤ 0,2jjгде e = (1, 0, .
. . , 0)⊤ .Тогда вычисление первого элемента в столбце A(j : n, j), т. е. тогоединственного элемена, который останется ненулевым, не будетсопровождаться вычитанием чисел одного знака и, как следствие,возможной потерей точности.Ещё одно соображение по практической реализации описанноговыше алгоритма QR-разложения с помощью матриц Хаусхолдера.В действительности даже не нужноформировать в явном виде матрицу отражения H̃:умножение на неё можно выполнить по экономичной формулеI − 2uu⊤ A(j : n, j : n)= A(j : n, j : n) − 2u u⊤ A(j : n, j : n) .Метод отражений ХаусхолдераОпределённым недостатком метода Хаусхолдера в сравнениис методом Гаусса является привлечение операции извлеченияквадратного корня, которая приводит к иррациональностям.Метод отражений ХаусхолдераОпределённым недостатком метода Хаусхолдера в сравнениис методом Гаусса является привлечение операции извлеченияквадратного корня, которая приводит к иррациональностям.Это не позволяет точно (без округлений) реализоватьсоответствующие алгоритмы в поле рациональных чисел.Как следствие, нельзя воспользоваться программными системамитак называемых «безошибочных вычислений» или языкамипрограммирования типа Ruby, которые могут точно оперироватьрациональными дробями.Матрицы вращения и метод вращенийМатрицей вращения называется n × n-матрица G(k, l, θ) вида1...cos θ···− sin θk-ая строка 1......,...1l-ая строка sin θ···cos θ...1где 1 ≤ k < l ≤ n, 0 ≤ θ < 2π.Остальные не выписанные явно элементы вне главной диагоналиравны нулю.Таким образом, G(k, l, θ) — это матрица, которая отличается отединичной матрицы лишь элементами, находящимися в позициях(k, k),(k, l),(l, k),(l, l).Таким образом, G(k, l, θ) — это матрица, которая отличается отединичной матрицы лишь элементами, находящимися в позициях(k, k),(k, l),(l, k),(l, l).Нетрудно проверить, что G(k, l, θ) — ортогональная:⊤⊤= I.G(k, l, θ) G(k, l, θ) = G(k, l, θ) G(k, l, θ)Матрицаcos θ − sin θsin θcos θ!задаёт, как известно, вращение двумерной плоскости 0x1 x2на угол θ вокруг начала координат.Матрицаcos θ − sin θsin θcos θ!задаёт, как известно, вращение двумерной плоскости 0x1 x2на угол θ вокруг начала координат.Матрица G(k, l, θ) также задаёт вращение пространства Rn на уголθ вокруг оси, проходящей через начало координат и ортогональнойгиперплоскости 0xk xl .Матрицы вращения G(k, l, θ) называют также матрицами Гивенса.Мы будем иногда обозначать их посредством G(k, l),если конкретная величина угла θ несущественна.Подходящим вращением можно занулить любую из компонентдвумерного вектора:x20x1Если вектор a = (a1 , a2 )⊤ — ненулевой, то, взявcos θ =a1,kak2sin θ =−a2,kak2где kak2 =qa21 + a22 ,мы можем с помощью матрицы двумерного вращения занулитьвторую компоненту этого вектора:!!!kak2cos θ − sin θa1.=0a2sin θcos θЕсли вектор a = (a1 , a2 )⊤ — ненулевой, то, взявcos θ =a1,kak2sin θ =−a2,kak2где kak2 =qa21 + a22 ,мы можем с помощью матрицы двумерного вращения занулитьвторую компоненту этого вектора:!!!kak2cos θ − sin θa1.=0a2sin θcos θАналогично может быть занулена первая компонента вектора a,путём домножения на такую матрицу вращения чтоcos θ =a2,kak2sin θ =a1.kak2В общем случае для любой матрицы A = (aij ) в произведенииà = (ãij ) := G(k, l, θ) Aстроки k-ая и l-ая становятся линейными комбинациями строк сэтими же номерами из A:ãkj ← akj cos θ − alj sin θ,ãlj ← akj sin θ + alj cos θ,j = 1, 2, .
. . , n.Остальные элементы матрицы à совпадают с элементами в A.В общем случае для любой матрицы A = (aij ) в произведенииà = (ãij ) := G(k, l, θ) Aстроки k-ая и l-ая становятся линейными комбинациями строк сэтими же номерами из A:ãkj ← akj cos θ − alj sin θ,ãlj ← akj sin θ + alj cos θ,j = 1, 2, . . .
, n.Остальные элементы матрицы Ã совпадают с элементами в A.Из рассуждений предшествующего абзаца вытекает, что путёмспециального подбора угла θ можно занулить любой элемент k-ойили l-ой строк в произведении Ã = G(k, l, θ) A.Следствие.Любая квадратная матрица A может быть приведена к правомутреугольному виду с помощью последовательности умноженийслева на матрицы вращения.Следствие.Любая квадратная матрица A может быть приведена к правомутреугольному виду с помощью последовательности умноженийслева на матрицы вращения.Доказательство.Мы можем один за другим занулить поддиагональные элементыпервого столбца, потом второго, третьего и т. д., аналогично тому,как это делалось в прямом ходе метода Гаусса.При этом зануление поддиагональных элементов второго ипоследующих столбцов никак не испортит полученные ранеенулевые элементы предшествующих столбцов, так как линейноекомбинирование нулей даст снова нуль.Формально можно сказать, что существует набор матриц вращенияG(1, 2), G(1, 3), .