1611688890-f641c9ec8276824e4686da772eb56520 (826652), страница 53
Текст из файла (страница 53)
. . , 0)⊤ ;H̃ ← I − 2uu⊤ ;ELSEH̃ ← I ;END IFA(j : n, j : n) ← H̃ A(j : n, j : n) ;END DOвытекает равенство A = QR с ортогональной матрицейQ = Hn−1 Hn−2 · · · H2 H1−1.Таким образом, мы получаем QR-разложение матрицы A, т. е. Предложения 3.7.2 и 3.7.3 дают в совокупности ещё одно, конструктивное,доказательство Теоремы 3.7.1.
Соответствующий псевдокод алгоритмадля вычисления QR-разложения матрицы приведён в Табл. 3.2.Как следствие, исходная система уравнений Ax = b становится равносильной системе уравнений(Qy = b,Rx = y,с несложно решаемыми составными частями. При практической реализации удобнее дополнить алгоритм Табл. 3.2 инструкциями, которые3.7. Методы на основе ортогональных преобразований323задают преобразования вектора правой части СЛАУ, и тогда результатом работы нового алгоритма будет правая треугольная СЛАУ Rx = y.Её можно решать с помощью обратной подстановки (3.58).Согласно Предложению 3.7.2 вычисление вектора Хаусхолдера u вкачестве первого шага требует нахождения из (3.85) вектора ũ, в котором имеется неоднозначность выбора знака второго слагаемого.
Привычислениях на цифровых ЭВМ в стандартной арифметике с плавающей точкой имеет смысл брать(A(j : n, j) + kA(j : n, j)k2 e, если ajj ≥ 0,ũ =A(j : n, j) − kA(j : n, j)k2 e, если ajj ≤ 0,где e = (1, 0, . . . , 0)⊤ . Тогда вычисление первого элемента в столбцеA(j : n, j), т. е. того единственного элемена, который останется ненулевым, не будет сопровождаться вычитанием чисел одного знака и, какследствие, возможной потерей точности.Ещё одно соображение по практической реализации описанного вПредложении 3.7.3 алгоритма состоит в том, что в действительностидаже не нужно формировать в явном виде матрицу отражения H̃:умножение на неё можно выполнить по экономичной формулеI − 2uu⊤ A(j : n, j : n)= A(j : n, j : n) − 2u u⊤ A(j : n, j : n) .Определённым недостатком метода Хаусхолдера и описываемого вследующем пункте метода вращений в сравнении с методом Гаусса является привлечение неарифметической операции извлечения квадратного корня, которая приводит к иррациональностям. Это не позволяетточно (без округлений) реализовать соответствующие алгоритмы в поле рациональных чисел, к примеру, в программных системах так называемых «безошибочных вычислений» или языках программированиятипа Ruby [83], которые могут оперировать рациональными дробями счислителем и знаменателем в виде целых чисел.3.7дМатрицы вращения и метод вращенийПусть даны натуральные числа k, l, не превосходящие n, т.
е. размерности пространства Rn , и задано значение угла θ, 0 ≤ θ < 2π.3243. Численные методы линейной алгебрыМатрицей вращения называется n × n-матрица G(k, l, θ) вида1...cosθ···−sinθk-ая строка 1......,...1l-ая строкаsin θ···cos θ...(3.86)1где все не выписанные явно элементы вне главной диагонали равнынулю. Таким образом, G(k, l, θ) — это матрица, которая отличаетсяот единичной матрицы лишь элементами, находящимися в позициях(k, k), (k, l), (l, k) и (l, l) для данных индексов k, l. Нетрудно проверить,что она ортогональна.Матрица!cos θ − sin θ(3.87)sin θcos θзадаёт, как известно, вращение двумерной плоскости 0x1 x2 на угол θвокруг начала координат.18 Матрица G(k, l, θ) также задаёт вращениепространства Rn на угол θ вокруг оси, проходящей через начало координат и ортогональной гиперплоскости 0xk xl . Матрицы вращенияG(k, l, θ) называют также матрицами Гивенса, и мы будем иногда обозначать их посредством G(k, l), если конкретная величина угла θ несущественна.Если вектор a = (a1 , a2 )⊤ — ненулевой, то, взявqa1−a2cos θ =,sin θ =,где kak2 = a21 + a22 ,kak2kak2мы можем с помощью матрицы двумерного вращения (3.87) занулитьвторую компоненту этого вектора:!!!kak2a1cos θ − sin θ.=0a2sin θcos θ18 Напомним, что положительным направлением вращения плоскости считаетсявращение «против часовой стрелки».3.7.
Методы на основе ортогональных преобразований325x20x1Рис. 3.18. Подходящим вращением можно занулитьлюбую из компонент двумерного вектора.Аналогично может быть занулена первая компонента вектора a, путёмдомножения на такую матрицу вращения (3.87), чтоcos θ =a2,kak2sin θ =a1.kak2В общем случае умножение любой матрицы A = (aij ) слева наматрицу вращения G(k, l, θ) приводит к тому, что в их произведенииà = (ãij ) := G(k, l, θ) A строки k-ая и l-ая становятся линейными комбинациями строк с этими же номерами из A:ãkj ← akj cos θ − alj sin θ,ãlj ← akj sin θ + alj cos θ,(3.88)j = 1, 2, .
. . , n. Остальные элементы матрицы Ã совпадают с элементами матрицы A. Из рассуждений предшествующего абзаца вытекает,что путём специального подбора угла θ можно занулить любой элементk-ой или l-ой строк матрицы Ã = G(k, l, θ) A.Как следствие, любая квадратная матрица A может быть приведена к правому треугольному виду с помощью последовательности умножений слева на матрицы вращения. Более точно, мы можем один задругим занулить поддиагональные элементы первого столбца, потомвторого, третьего и т.
д., аналогично тому, как это делалось в прямомходе метода Гаусса. При этом зануление поддиагональных элементов3263. Численные методы линейной алгебрывторого и последующих столбцов никак не испортит полученные ранеенулевые элементы предшествующих столбцов, так как линейное комбинирование нулей согласно формулам (3.88) даст снова нуль. На формальном языке можно сказать, что существует набор матриц вращенияG(1, 2), G(1, 3), .
. . , G(1, n), G(2, 3), . . . , G(n − 1, n), таких чтоG(n − 1, n) · · · G(2, 3) G(1, n) · · · G(1, 3) G(1, 2) A = R— правая треугольная матрица. ОтсюдаA = G(1, 2)⊤ G(1, 3)⊤ · · · G(1, n)⊤ G(2, 3)⊤ · · · G(n − 1, n)⊤ R,и мы получили QR-разложение матрицы A, поскольку произведениетранспонированных матриц вращения также является ортогональнойматрицей.Использование преобразований вращения — ещё один конструктивный способ получения QR-разложения, технически даже более простой, чем метод отражений Хаусхолдера. При его реализации организовывать полноценные матрицы вращения G(k, l, θ) и матричные умножения с ними, конечно, нецелесообразно, так как большинством элементов в G(k, l, θ) являются нули. Результат умножения слева на матрицу вращения разумно находить путём перевычисления лишь ненулевых элементов всего двух строк по формулам (3.88).Для плотно заполненных матриц использование вращений в полтора раза более трудоёмко, чем получение QR-разложения с помощьюматриц отражения, но зато вращения более предпочтительны для разреженных матриц в силу своей большей гибкости при занулении отдельных элементов.3.7еПроцессы ортогонализацииОртогонализацией называют процесс построения по заданному базису линейного пространства некоторого ортогонального базиса, который имеет ту же самую линейную оболочку.
Ввиду удобства ортогональных базисов для представления решений разнообразных задач и,как следствие, их важности во многих приложениях (см., к примеру,§2.10г) огромное значение имеют и процессы ортогонализации.Исторически первым процессом ортогонализации был алгоритм, который по традиции связывают с именами Й. Грама и Э. Шмидта.19 По19 Иногдаэтот процесс называют «ортогонализацией Сонина-Шмидта».3.7. Методы на основе ортогональных преобразований327конечной линейно независимой системе векторов {v1 , v2 , . . . , vn } процесс Грама-Шмидта строит ортогональный базис {q1 , q2 , . . .
, qn } линейной облочки векторов {v1 , v2 , . . . , vn }.Возмём в качестве первого вектора q1 конструируемого ортогонального базиса вектор v1 , первый из исходного базиса. Далее для построения q2 можно использовать v2 «как основу», но откорректировав его сучётом требования ортогональности к q1 и принадлежности линейнойоболочке векторов q1 = v1 и v2 . Естественно положить q2 = v2 − α12 q1 ,где коэффициент α12 подлежит определению из условия ортогональностиhq1 , v2 − α12 q1 i = 0.Отсюдаhq1 , v2 i.hq1 , q1 iДалее аналогичным образом находится q3 = v3 − α13 q1 − α23 q2 , и т. д.В целом ортогонализация Грама-Шмидта выполняется в соответствии со следующими расчётными формулами:α12 =qj ← vj −j−1Xhqk , vj iqk ,hqk , qk ij = 1, 2, .
. . , n(3.89)k=1(суммирование не выполняется, если его верхний предел меньше нижнего). В случае, когда очередной вычисленный вектор qj — нулевой,ясно, что v1 , v2 , . . . , vj линейно зависимы. Тогда построение для линейной оболочки lin {v1 , v2 , . . . , vn } ортогонального базиса, состоящегоиз того же количества векторов, невозможно.Если получающиеся векторы ортогонального базиса дополнительнонормируются сразу после их нахождения, то hqk , qk i = 1 в (3.89), чтоупрощает выражения для поправочных коэффициентов αkj . Псевдокод соответствующего варианта ортогонализации Грама-Шмидта данв Табл. 3.3.Каково матричное представление процесса ортогонализации ГрамаШмидта? Пусть векторы v1 , v2 , .
. . , vn заданы своими координатнымипредставлениями в некотором базисе, и из вектор-столбцов этих координатных представлений мы организуем матрицу W . В результате ортогонализации мы должны получить ортогональную матрицу, в которой первый столбец — это нормированный первый вектор, второй столбец — это нормированная линейная комбинация первых двух векторстолбцов, и т.