Вычислительные методы алгебры и оценивания. И.В. Семушкин (2011) (1185350), страница 20

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 20)

7.3) для определения вектора u. Второй и далее столбцы,обозначенные на рис. 7.3 произвольно как вектор a, отражают от найденнойтаким образом гиперплоскости U⊥ , решая для этого задачу 1 (см. рис. 7.2)с y := a и тем самым получая блок Ã.Теорема 7.1 (Триангуляризация матрицы по методу Хаусхолдера).Пусть A1 := A(m, n) и для каждого j выбрано элементарное преобразованиеХаусхолдера Tj так, что1n−jz}|{ z }| {T1{saj jT j Aj = , j = 1, .

. . , k;m−j0 Aj+1114k ≤ min (m − 1, n). (7.19)7.4 Преобразование ХаусхолдераТогда в процессе после k повторных применений свойства 6 и леммы 7.1имеем следующий промежуточный результат триангуляризации матрицы A:aT1s1aT2s2...aTk...T (k) A =sk0(7.20)Ak+1с отвечающей этому моменту процесса итоговой матрицей преобразованийI0I0k−1T (k) =··· 1T1 .(7.21)0Tk0T2Замечание 7.4.

Важно подчеркнуть, что алгоритм триангуляризации (7.19) не требует вычисления или запоминания ортогональной матрицыT (k) , так как правая часть равенства (7.4) вычисляется непосредственно всоответствии с замечанием 7.3. Стоит также заметить, как неявно определяется Aj+1 рекурсией по j в алгоритме (7.19). Кроме Aj+1, на шаге j этойрекурсии определяются скаляр sj и (n−j) компонент вектор-строки aTj . Этинеявные соотношения для sj , aTj и Aj+1 и весь процесс вычислений (рис. 7.4)представлены в явном виде в подразд. 7.5.mn}|z{0|{zm}(a)nn}|z0(b){mzn}|0(c){nmn}|z{0| {z }k(d)Рис.

7.4. Представление возможных случаев применения теоремы 7.1 к матрицеA(m, n); (a) недоопределенная система: k = m − 1 ≤ n; (b) определенная система:k = n − 1, m = n; (c) переопределенная система: k = n < m; (d) k < n < m1157 Ортогональные преобразования7.5Шаг триангуляризации матрицы методом ХаусхолдераПусть матрица A = A(m, n) задана.

Тогда, согласно лемме 7.1, базоваяоперация процесса триангуляризации заключается в вычислении скаляра sи матрицы à = Ã(m, n − 1) таких, что1n−1z}|{ z }| {1{s à  .Tu A =m−10 (7.22)Алгоритм. Для вычисления s следует применять свойство 6 (см.стр. 113), т. е. выполнять (7.23) для всех k = 1 до min (m − 1, n). Затемдля вычисления à следует применять свойство 5 (см. стр. 112), т. е. последовательно в цикле по j = 2, . .

. , n для всех i = k, . . . , m выполнять (7.24).Здесь λ , −γ (см. (7.15)), α , −β (см. (7.14) и использовано свойство 6).Для k = 1 до min (m − 1, n)!1/2mX2[A(i, k)] ,sk = − sgn [A(k, k)]i=k(7.23)uk (1) = A(k, k) − sk ,uk (i) = A(k + i − 1, k), i = 2, . . . , m − k + 1, α = 1/(s u (1))(α < 0).kk kkДля j = k + 1, . . . , nmXuk (i − k + 1)A(i, j),λ := αk ·i=kДля i = k, k + 1, . .

. , mA(i, j) := A(i, j) + λuk (i − k + 1).(7.24)Приведенный алгоритм (7.23), (7.24) называют столбцово ориентированным алгоритмом Хаусхолдера, так как операции (7.24) вычисляют целикомкаждый j-й столбец матрицы, находящийся справа от ведущего, т. е. k-гостолбца. Альтернативная схема вычислений называется строчно ориентированным алгоритмом Хаусхолдера. Ее можно получить из выражения Tu == I − βuuT для матрицы Хаусхолдера следующим образом.1167.6 Решение треугольной системы и обращение матриц√Введем вспомогательные обозначения: µ , β, w , µu, чтобы записать Tu = I − wwT .

Тогда (TuA) = A − wz T , где z T , wT A = µv T ,PvT , mi=1 u(i)A(i, ·) и A(i, ·) есть i-я строка матрицы A = A(m, n). Введем обозначение λT = αv T , используя ранее введенное (см. (7.23)) α , −β.Отсюда получаем формулу для любой i-й строки (TuA)(i, ·) преобразованнойматрицы (TuA) в виде(TuA)(i, ·) = A(i, ·) − w(i)z T = A(i, ·) − µ2 u(i)v T = A(i, ·) + λT u(i).Алгоритм (строчно ориентированный), эквивалентный (7.23) и (7.24).Для k = 1 до min (m − 1, n)sk = − sgn [A(k, k)]mXi=k!1/2[A(i, k)]2 ,uk (1) = A(k, k) − sk ,uk (i) = A(k + i − 1, k), i = 2, . . . , m − k + 1,αk = 1/(sk uk (1))(αk < 0).Для j = k + 1, .

. . , nmXuk (i − k + 1)A(i, j),λk (j − k) := αk ·i=kДля i = k, k + 1, . . . , mДля j = k + 1, . . . ,nA(i, j) := A(i, j) + λk (j − k)uk (i − k + 1).7.6(7.25)(7.26)Решение треугольной системы и обращение матрицКак отмечено в подразд. 7.3, мы часто заинтересованы в решении уравненияRx = z,(7.27)где R = R(n, n) — верхняя треугольная невырожденная матрица.

Еслинужно иметь только решение x, то R−1 (для x = R−1z) вычислять не надо.Следующий алгоритм обратной подстановки позволяет вычислить решениеx непосредственно.1177 Ортогональные преобразованияАлгоритм. Для j = n, n − 1, . . . , 1 вычислятьnXR(j, k)x(k) /R(j, j).x(j) = z(j) −(7.28)k=j+1По сложности этот алгоритм почти такой же, как матричное умножение.Он допускает записывать x(j) поверх z(j), что очень удобно в приложениях.Если все же требуется иметь матрицу U , R−1, то ее можно вычислятьпо алгоритму окаймления, основанному на следующем легко проверяемомтождестве для треугольных матриц:−1 −1−1Rj−Rj−1yσj+1Rjy−1= Rj+1.(7.29)=−10σj+10σj+1Это соотношение позволяет вычислять матрицу U , R−1 рекуррентно, аименно: если Rj−1 = Uj , где Rj обозначает верхнюю левую часть матрицыR, тоUj −Uj [R(1, j + 1), .

. . , R(j, j + 1)]T σj+1Uj+1 =,(7.30)0σj+1где σj+1 = 1/R(j + 1, j + 1). Полагая U = R−1, этот результат представим валгоритмической форме.Алгоритм. Обращение верхней треугольной матрицы. Задать начальное значениеU (1, 1) = 1/R(1, 1) .(7.31)Для j = 2, . . . , n вычислять по формулам (7.32) и (7.33):U (k, j) = −j−1Xi=kU (j, j) = 1/R(j, j) ,!U (k, i)R(i, j) U (j, j) ,(7.32)k = 1, . . .

, j − 1 .(7.33)Замечание 7.5. R−1 вычисляется по столбцам, при этом элементыматрицы R−1 могут записываться поверх элементов исходной матрицы R.В справочных целях приведем примеры реализации данного алгоритмана языке FORTRAN. Эти примеры могут помочь студентам написать своисобственные программы на других языках высокого уровня при выполнениилабораторного проекта № 6, описание которого дано ниже в подразд.

7.16.1187.6 Решение треугольной системы и обращение матрицОбращение верхней треугольной матрицы: U := R−1. Реализуются формулы (7.31), (7.32) и (7.33). Если нужно, U может замещать R [15].R(N, N ),U (N, N ),R и U — верхние треугольные матрицыU (1, 1) = 1./R(1, 1)DO 20 J = 2, NU (J, J) = 1./R(J, J)JM1 = J − 1DO 20 K = 1, JM1SU M = 0.DO 10 I = K, JM110 SU M = SU M − U (K, I) ∗ R(I, J)20 U (K, J) = SU M ∗ U (J, J)В случаях, когда важно или нужно экономить память компьютера, матрицы в программе объявляют как одномерные массивы (см. подразд. 6.3).Хотя в компьютере даже многомерно объявленные массивы всегда хранятсякак одномерные, компилятор генерирует индексные выражения с операциями умножения и деления.

Операции сложения и вычитания в компьютерахвыполняются гораздо быстрее, поэтому индексы для доступа к элементамматриц следует программировать в рекуррентной инкрементной форме, экономя таким образом время процессора (табл. 7.1). В этой программе преобразование в треугольную форму выполняется отождествлением J(J − 1)/1 + Iс (I, J). Рекуррентное инкрементное вычисление KK, JJ и KK экономитвычисления.Как отмечалось на с. 118, иногда требуется вычислять R−1. Такая ситуация появляется, если требуется найти A−1, для которой уже выполненопреобразование T A = R, где T = T (n−1) по формуле (7.21), так как в теореме 7.1 для этого случая m = n и A−1 = R−1T . Последнее означает, что тоже самое ортогональное преобразование T теперь надо применить к строкам матрицы R−1 , но уже в обратном порядке следования элементарныхпреобразований, составляющих полное преобразование T = T (n−1) по формуле (7.21).

Таким образом, возникает проблема запоминания элементарныхпреобразований, составляющих полное преобразование T = T (n−1) , чтобыпозже можно было его применить в задаче отыскания A−1 или же для решения уравнения Ax = z с невырожденной матрицей A после преобразованияT A = R.1197 Ортогональные преобразованияТаблица 7.1.

Эффективное обращение верхней треугольной матрицы U := R−1 .R(N, N),U(1) = 1./R(1)JJ = 1DO 20 J = 2, NJJ = JJU(N, N),R и U — верхнетреугольные матрицы.c R(1) ≡ R(1, 1)c JJ = J(J − 1)/2 = (J −− 1, J − 1)c JJ = J(J + 1)/2 = (J, J)JJ = JJ + JU(JJ) = 1./R(JJ)JM1 = J − 1KK = 0DO 20 K = 1, JM1c KK = K(K + 1)/2KK = KK + KKI = KKSUM = 0.DO 10 I = K, JM1c KI = (K, I), JJ + 1 =SUM = SUM − U(KI) ∗ R(JJ + I) = (I, J)c KI = (K, I + 1)10 KI = KI + Ic JJ + K = (K, J)20 U(JJ + K) = SUM ∗ U(JJ)Реализуются формулы (7.31), (7.32) и (7.33). Верхние треугольные матрицы R и Uхранятся как векторы размерности N(N + 1)/2. Если нужно, U может замещать R.Как видно из (7.24), для отражения вектора y = T (k) z от гиперплоскости U⊥ , заданной вектором u, требуется иметь сохраненными две величины: вектор u и скаляр α.

Поскольку нули ниже диагонали, получающиесяв результате триангуляризации, хранить не нужно, это место можно отвестидля сохранения вектора u (вместе с диагональю, поэтому диагональные элементы sk приходится хранить отдельно). Данное предложение схематическипоказано на рис. 7.5 (вверху). Каким образом можно выполнить вычислениематрицы A−1, показано на рис. 7.5 (внизу) на конкретном примере размерностей, а именно: m = 4, n = 4.7.7Преобразование ГивенсаПреобразование Гивенса осуществляет вращение векторав плоскостиTдвух координат. Очевидно, поворот вектора y = (y1 y2) на угол θ почасовой стрелке эквивалентен повороту системы координат против часовойстрелки на тот же угол.1207.7 Преобразование ГивенсаT (k) A =k⇓T (k) zTzRykaTkz1ymz2α1 · · · αk · · · αnTA =s1 · · · sk · · · snY := R−1Au1aT2u2u3α1α2α3s1s2s3aT1aT3s4−→&—–—y1T —–— —–—y2T —–—  —–—y3T —–— —–—y4T —–—R=s1↑s2s30i=kДля i = k до myi := yi + λuk (i − k + 1)...ukДля k = 1 до min (m − 1, n)mPλ = αkuk (i − k + 1)yi,→R0}n, Rx = z1}m− nДля k = n − 1, .

Характеристики

Список файлов книги