Вычислительные методы алгебры и оценивания. И.В. Семушкин (2011) (1185350), страница 23
Текст из файла (страница 23)
Обратим внимание на различия между двумятипами ортогональных преобразований, а именно: между преобразованиями Хаусхолдера и Гивенса (если говорить о левосторонних их версиях, хотявозможны и правосторонние) — это один тип преобразований, и ортогонализацией Грама–Шмидта — другой тип. Первый тип обеспечивает равенства(7.39), (7.40), (7.41) или (7.42), где преобразованная матрица QA имеет тотже размер, что и исходная матрица A, и в ее составе присутствует блок,появляющийся как треугольная матрица R в одном из четырех углов этойматрицы QA.
При этом матрица Q ортогонального преобразования — квадратная. В случае ортогонализации Грама–Шмидта имеем A = QR, где Q —матрица того же размера, что и A, но с ортонормированными столбцами, аR — строго квадратная матрица с треугольным заполнением.7.12Алгоритмы ортогонализации Грама–ШмидтаРассмотрим задачу QR-разложения матрицы A(m, n), m ≥ n, полногоранга на основе ортогонализации Грама–Шмидта.В данной задаче, рассматривая пример n = 3, имеемr11 r12 r13r22 r23 ,R=r331337 Ортогональные преобразованияA = [a1, a2 , a3 ] = [q1, q1, q3]R = [r11q1 , r12q1 + r22q2, r13q1 + r23q2 , r33q3].В результате получаем линейную системуr11q1 = a1 ,r12q1 + r22q2 = a2 ,r13q1 + r23q2 + r33q3 = a3 .Ясно, что каждое из этих выражений представляет собой разложение вектора a1 , a2 или a3 по системе ортов {q1, q2, q3}, при этом коэффициент rijесть алгебраическая проекция вектора aj на орт qi.В силу треугольности матрицы R эта система легко решается.
Из первогоуравнения находим орт q1 вдоль вектора a1 и координату (проекцию какчисло) первого вектора a1 вдоль орта q1 :q1 = a1 /ka1k, r11 = ka1 k.Второе уравнение есть разложение вектора a2 на сумму проекций вдольортов q1 и q2. Так как орт q1 уже найден, то координата r12 легко определяется в видеr12 = aT2 q1.После этого из второго уравнения имеемr22q2 = a2 − r12q1 ,следовательно,q2 = (a2 − r12q1 )/ka2 − r12q1k,r22 = ka2 − r12q1 k.Замечание 7.8. По предположению, rank A = n, т. е. ни один извекторов ai не является нулевым, и все ai образуют линейно независимуюсистему.
Поэтому r11 6= 0, r22 6= 0 и r33 6= 0, следовательно, существует R−1 .Продолжая решение системы, для третьего уравнения находимr13 = aT3 q1 , r23 = aT3 q2и затем определяемОтсюдаr33q3 = a3 − r13q1 − r23q2 .q3 = (a3 − r13q1 − r23q2 )/ka3 − r13q1 − r23q2 k,r33 = ka3 − r13q1 − r23q2 k.1347.12 Алгоритмы ортогонализации Грама–ШмидтаТаким образом, получили классический метод ГШО, отличающийся тем,что матрица R определяется по столбцам с номерами k = 1, 2, .
. . , n.В настоящее время существуют две более эффективные версии — ортогонализации Грама–Шмидта:•Алгоритм МГШО (модифицированая схема).•Алгоритм МГШО с выбором ведущего вектора.Первые две версии — классическая и модифицированная — показаныздесь — на стр. 135, а модифицированная с выбором ведущего вектора —на стр.
136.Алгоритм ГШО (классическая схема)Для k = 1 до nrik = aTk qi , i = 1, 2, . . . , k − 1,v = ak −k−1Prik qi,i=1rkk = (v T v)1/2,qk = v/rkk .Алгоритм МГШО (модифицированая схема)Для k = 1 до nrkk = kak k = (aTk ak )1/2,ak = ak /rkk ,Для j = k + 1 до nrkj = aTj ak ,aj = aj − rkj ak .1357 Ортогональные преобразованияАлгоритм МГШО с выбором ведущего вектораДля k = 1 до nq(k) = k,Для k = 1 до nДля s = k до nНайти № l, для которого kaq(l) k = max kaq(s) k,sПереставить номера: q(k) q(l),rq(k),q(k) = kaq(k) k = (aTq(k) aq(k) )1/2,aq(k) = aq(k) /rq(k),q(k) ,Для j = k + 1 до nrq(k),q(j) = aTq(j) aq(k) ,aq(j) = aq(j) − rq(k),q(j) aq(k) .Первая из более современных версий, называемая МГШО (Rice, 1966[11]), отличается порядком вычислений матрицы R.
В этой версии матрица Rопределяется по строкам с номерами k = 1, 2, . . . , n. Этот алгоритм требуетменьше оперативной памяти, так как в нем не используется промежуточныйвектор v. Кроме того, матрица A заменяется матрицей Q, потому что послеоперации деления имеем ak = qk . Одним из его преимуществ является то,что в него легко внедрить процедуру выбора ведущего столбца.Чтобы получить вторую из более современных версий — так называемую МГШО с выбором ведущего вектора, — нужно изменить алгоритмМГШО таким образом, чтобы очередным ортогонализируемым векторомоказался не k-й, а тот, чья норма наибольшая среди всех оставшихся s-хвекторов от s = k до s = n.
Как и в подразд. 2.2, реально переставляются нестолбцы матрицы A, а обмениваются значениями только элементы дополнительного вектора q, в котором фиксируются номера столбцов матрицыA. Доступ к элементам матрицы A осуществляется с использованием этоговектора. Перед началом работы основного алгоритма ортогонализации этотвектор перестановок q заполняется числами от 1 до n в естественном порядкенумерации столбцов матрицы A.1367.13 Решение систем после ортогонализации7.13Решение систем после ортогонализации1.
Пусть дана система линейных алгебраических уравнений с квадратнойневырожденной матрицей Ax = f . Тогда после ортогонального приведения матрицы с помощью одной из версий ортогонализации Грама–Шмидтаимеем представление этой матрицы в виде A = QR и, следовательно,QRx = f и Rx = QT f .2. Пусть дана система линейных алгебраических уравнений с прямоугольной матрицей A(m, n), m > n, полного ранга. Такая система называетсяпереопределенной системой. Нормальное псевдорешение x̄, найденное пометоду наименьших квадратов (МНК), удовлетворяет нормальным уравнениямAT Ax̄ = AT f.Поскольку A = QR и QT Q = I, эти уравнения эквивалентны уравнениюRx̄ = QT f,которое совпадает по виду с уравнением из п. 1.Чтобы вычислить x (для п.
1) или x̄ (для п. 2), находят вектор f 0 = QT f ,а затем решают систему с треугольной матрицей R (методом подстановки).7.14Обращение матриц после ортогонализацииДля матрицы A = A(n, n) имеем A = QR, где Q = Q(n, n). ОтсюдаA−1 = R−1Q−1 = R−1QT .Следовательно, A−1 есть решение матричного уравнения RX = QT .Чтобы найти i-й столбец матрицы A−1, надо в качестве правой части взятьi-й столбец матрицы QT и решить систему с треугольной матрицей R (какв подразд. 7.13 или подробнее в подразд. 7.6).7.15Задание на лабораторный проект № 6Написать и отладить программу, реализующую ваш вариант ортогонального преобразования для численного решения систем линейныхалгебраических уравнений Ax = f с квадратной матрицей A, вычисления1377 Ортогональные преобразования±(det(A)) и A−1.
Предусмотреть предупреждение о невозможности решения указанных задач из-за присутствия (почти) линейно зависимых векторов среди столбцов матрицы A (в пределах ошибок округления ЭВМ илидругого, заранее определенного критерия). Отделить основные части программы:а) подпрограмму факторизации матрицы A, отвечающую вашему варианту метода ортогонального приведения;б) подпрограмму решения систем линейных алгебраических уравнений;в) подпрограмму вычисления определителя матриц;г) подпрограмму обращения матриц;д) сервисные подпрограммы.Уделить особое внимание эффективности программы (в смысле экономииоперативной памяти и скорости решения указанных выше задач). Предусмотреть пошаговое выполнение алгоритма ортогонального приведения свыводом результата на экран.
Выполнить следующие пункты задания:1. Провести подсчет фактического количества операций, выполняемыхпри решении системы линейных алгебраических уравнений (отдельно числоопераций сложения, число операций умножения, число операций деленияи число операций извлечения квадратного корня) и сравнить эти числа стеоретическими (оценочными) числами.2.
Оценить скорость решения задач, т.е. определить время, затраченное нарешение системы линейных алгебраических уравнений, и время, затраченноена обращение матриц. Для этого спроектировать и провести эксперимент,который охватывает матрицы порядка от 10 до 100 (через 10 порядков).Представить результаты в виде таблицы и графика зависимости временивыполнения (в минутах и секундах) от порядка матриц. Таблицу и графиквывести на экран.3.
Оценить точность решения систем линейных алгебраических уравнений, имеющих тот же самый порядок, что и задачи из п. 2. Для этого сгенерировать случайные матрицы A, выбрать точное решение x∗ и образоватьправые части f = Ax∗. Провести анализ точности вычисленного решения xот порядка матрицы.
Результаты представить в виде таблицы и графика.Для заполнения матрицы A использовать случайные числа из диапазонаот −100 до 100. В качестве точного решения взять вектор x∗ = (1, 2, . . . , n),где n — порядок матрицы. Для оценки точности использовать норму вектораkxk∞ = max(|xi|).i1387.15 Задание на лабораторный проект № 6Таблица 7.2. Варианты задания на лабораторный проект № 6Вариантзаполненияматрицы R00a——c—d—e—bОтраженияХаусхолдераВращенияГивенсаОртогонализацияГрама–Шмидтаababcde, Rne12345670 , Rnw891011121314, Rse151617181920210 , Rsw22232425262728столбцово-ориентированный алгоритм;строчно-ориентированный алгоритм;классическая схема;модифицированая схема;модифицированая схема с выбором ведущего вектора.4.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
