Вычислительные методы алгебры и оценивания. И.В. Семушкин (2011) (1185350), страница 21
Текст из файла (страница 21)
. . , 2, 1Для r = 1 до nnPλ = αkuk (i − k + 1)yr (i)i=k←-aT1aT2 aT3 s4Для i = k до nyr (i) := yr (i) + λuk (i − k + 1)↓Y := A−1Рис. 7.5. Вверху: Сохранение преобразования T и вычисление вектора y = T z, ∀y ∈∈ Rm . Внизу: Вычисление матрицы A−1 после сохранения преобразования TСледовательно, легко найти (рис. 7.6), что координаты y10 , y20 повернутого вектора yr = (y10 y20 )T определяются в виде y10 = y1 cos θ + y2 sin θ,y20 = −y1 sin θ + y2 cos θ.Записывая это в матричной форме и требуя, чтобы поворот P1,2 в плоскости (e1 , e2) происходил до совмещения с первой координатной осью, получим qcsrc , cos θ = y1 /ryr =y = P1,2y =,, r , y12 + y22 ,−s c0s , sin θ = y2 /rгде, очевидно, матрица P1,2 плоского вращения ортогональна при любом θ.1217 Ортогональные преобразованияe2e02ye01e1θ0Рис.
7.6. Геометрия вращенийТриангуляризация матрицы преобразованиями Гивенса. Выбор θтакой, что вторая координата вектора yr становится равной нулю, используем для триангуляризации матрицы A(m, n). На первом шаге нужны преобразования, делающие равными нулю все элементы ниже первого диагонального элемента. Для этого, очевидно, нужно выполнить последовательно элементарные вращения P1,2, P1,3, . .
. , P1,m . Так определенные преобразованиявоздействуют на все столбцы матрицы, но только первый столбец, которыйуместно назвать ведущим столбцом в преобразовании Гивенса, приобретаетжелаемый вид.Лемма 7.2. Пусть дана матрица A(m, n) и y — ее ведущий (левый)столбец. Тогда существует ортогональное преобразование Гивенса, задаваемое матрицей P1 = P1 (m, m), такое, что1n−1z}|{ z }| {1{r à ,P1 A =m−10 (7.34)P1 = P1,m · · · P1,3P1,2 ,и матрицы P1,j определяются по алгоритму на рис. 7.7.Теорема 7.2 (Триангуляризация матрицы по методу Гивенса).Пусть A1 := A(m, n) и для каждого j = 1, 2, . . .
, k, k ≤ min (m − 1, n)серия элементарных преобразований Гивенса, задаваемая матрицей Pj размера (m + 1 − j) × (m + 1 − j), выбрана, как сказано ниже.1227.7 Преобразование ГивенсаyP1,i (i = 2, 3, . . . , m)11ic1,is1,i....0.c1,i−s1,iy1r1,1 := y1.1im←10...m1Для i = 2, 3, . . . , mq2r1,i := r1,i−1+ yi2rc1,i := 1,i−1r1,iys1,i := r i1,imPr , r1,m =yi2 = kyk2i=1←yiymРис. 7.7. Вычисление матрицы P1,jДля ведущего (левого) столбца yj матрицы Aj эта Pj выбрана так, что1n−jz}|{ z }| {T1{rajjPj Aj = , j = 1, .
. . , k ;m−j0 Aj+1k ≤ min (m−1, n). (7.35)Тогда после k повторных применений леммы 7.2 имеем следующий промежуточный результат триангуляризации матрицы A:aT1r1aT2r2...P (k) A =0rk...aTk(7.36)Ak+1с отвечающей этому моменту процесса итоговой матрицей преобразований I0I0k−1P (k) =··· 1P1 , Pj = Pj,m−j+1 · · · Pj,3 Pj,2 , (7.37)0Pk0P21237 Ортогональные преобразованияyj⇓Pj,t (t = 2, 3, . . . , l), l , m − j + 1, j = 1, 2, . . . , k, i = j + t − 1.11tcj,isj,i....0.cj,i−sj,i←10y1,jrj,1 := y1,j.1tl...1lДля t = 2, . .
. , lq22rj,t := rj,t−1+ yt,jrcj,j+t−1 := j,t−1rj,tysj,j+t−1 := rt,jj,tlP2rj , rj,l =yt,j= kyj k2←t=1Формула (7.37) имеет рекуррентный вид произведенияI0j−1P (j) =P (j−1) , P (1) = P1 0Pj, N = min (m − 1, n).P =P· · · P P , j = 2, . . . , N jj,m−j+1yt,jyl,j(7.38)j,3 j,2Все участвующие здесь матрицы являются ортогональными, поэтому финальная матрица P , P (N ) также ортогональна.
Общее число используемыхпри этом элементарных матриц вращения равно (m − 1) + (m − 2) + . . . ++ (m − N ) = (2m − N − 1)N/2. В результате (в случае m > n) получимR, Rne ,(7.39)PA = ··· , R =00где индекс ne подчеркивает, что в треугольной матрице R заполненной частью может быть только «северо-восточная» (north-by-east) часть. ПолагаяQ = P T , при m = n имеем QR-разложение матрицы A = A(n, n), т.
е.A = QR. Матрицы в (7.37) и (7.38) непосредственно не вычисляются.Для алгоритма Гивенса — так же, как и для других матричных алгоритмов, — имеются две схемы вычислений: (1) строчно ориентированнаясхема Гивенса и (2) столбцово ориентированная схема Гивенса (рис. 7.8).Как и в алгоритме преобразований Хаусхолдера (см. рис. 7.5), здесь обычнотребуется сохранять информацию о произведенных по ходу алгоритма (7.38)элементарных преобразованиях.1247.7 Преобразование ГивенсаyjP (j) As⇓j⇓j ⇒ rji ⇒ ζj,iajsyjjaisyijms⇓j⇓j ⇒ rjyjjajsyijaisi ⇒ ζj,im(е)4352−|r|6178|r|(а)(в)Для j = 1 до min (m − 1, n)rj := ajjДля i = j + 1 до ma := rj ; b := aijполучить (r, c, s, ζ)из (a, b)(ajj , rj ) := r(aij , ζj,i) := ζcj,i := c; sj,i := sДля s = j + 1 до nДля i = j + 1 до mα = cj,iajs + sj,iaisais = −sj,iajs + cj,iaisajs := αДля j = 1 до min (m − 1, n)Для i = j + 1 до ma := ajj ; b := aijполучить (r, c, s, ζ)из (a, b)(ajj , rj ) := r(aij , ζj,i) := ζДля s = j + 1 до nα = cajs + saisais = −sajs + caisajs := α(б).
. . . . . . .Процедура . . . . . . . .получить (r, c, s, ζ) из (a, b)sgn [a],|a| > |b|σ=sgn [b],|a| ≤ |b|22 1/2r = σ(a + b )a/r,r 6= 0c=1,r=0b/r,r 6= 0s=0,r=0s,|a| > |b||a| ≤ |b| &ζ=1/c,& c 6= 01,c=0(г). . . . . . . .Процедура . . . . . . . .получить (c, s) из (ζ)c := 0ζ =1 ⇒s := 1d := |1 − ζ 2 |1/2c := d|ζ| < 1 ⇒s := ζc := 1/ζ|ζ| > 1 ⇒s := d/|ζ|(д)Для j = 1 до min (m − 1, n)Для i = j + 1 до mполучить (c, s) из ζj,iα = cyjj + syijyij = −syjj + cyijyjj := αРис.
7.8. Преобразование Гивенса: (a) столбцово ориентированная схема вычисления матрицы P A, где P = P (j) при j = min (m − 1, n) (нижняя матрица слева);(б) вычисление координаты r вектора (a, b)T , повернутого до совмещения с первойосью, а также косинуса и синуса угла поворота и рабочего признака ζ; (в) строчноориентированная схема вычисления матрицы P A (верхняя матрица (г) восстановление косинуса и синуса угла поворота из признака ζ; (д) получение вектора y темипреобразованиями Pj,i произвольного вектора z ∈ Rm , которые сохранены в рабочих признаках ζj,i и восстанавливаются из них; (е) вследствие п. (б) векторы 1, 2,3 и 4 поворачиваются к положительному направлению первой координатной оси, авекторы 5, 6, 7 и 8 — к отрицательному направлению этой оси.1257 Ортогональные преобразованияСохранение этой информации позволит впоследствии решать системыуравнений Ax = z (совместные или несовместные, в последнем случае —по методу наименьших квадратов, см.
подразд. 7.3) или же находить обратную матрицу A−1 (когда m = n).Необходимая информация — это значения косинуса и синуса, однако ихсохранение было бы неэффективным решением. Gentleman (1973) предложил способ [11], включенный в рис. 7.8(б) и (г) с геометрической иллюстрацией его действия на рис. 7.8(е). Введенный им рабочий признак ζ — этоодно число, которое можно хранить в позиции (i, j) как ζj,i вместо нулевогоэлемента, появляющегося в позиции (i, j) матрицы (7.39) в момент преобразования Pj,t (t = i + 1 − j) в (7.37).
Как и с преобразованиями Хаусхолдера,нахождение A−1 после преобразований Гивенса требует такой же последовательности процедур: сначала находят R−1 (см. рис. 7.5 (внизу)), затем к R−1применяют с правой стороны финальное преобразование P , P (N ) (7.38),так как A−1 = R−1P . Для этого для рис. 7.5 (внизу) надо взять алгоритм изрис.
7.8(д), который также отыскивает P z при решении уравнения Ax = z.7.8Варианты заполнения матрицы RТрадиционно ортогональные преобразования (выше рассмотрены T —преобразование Хаусхолдера и P — преобразование Гивенса) приводят матрицу к виду, показанному на рис. 7.4 или в выражении (7.39). Однако выбортого угла матрицы, который должен остаться треугольно заполненым, естественно, произволен. Предпочтения диктуются целями использования, т.
е.предназначением преобразования. Преследуя здесь учебно-тренировочныецели, включим в проект (см. подразд. 7.16) все четыре возможных варианта заполнения матрицы R. Вариант 1 показан в (7.39), другие три имеютследующий вид. Вариант 2:RQA = · · · ,0R=0, Rnw ,(7.40)где Q обозначает либо T (преобразование Хаусхолдера), либо P (преобразование Гивенса), индекс nw подчеркивает, что в треугольной матрице Rзаполненной частью может быть только «северо-западная» (north-by-west)1267.9 Правосторонние ортогональные преобразованиячасть.
Вариант 3:0QA = · · · ,RR=0, Rse,(7.41)где индекс se подчеркивает, что в треугольной матрице R заполненнойчастью может быть только «юго-восточная» (south-by-east) часть. Вариант 4:00PA = ··· , R =, Rsw ,(7.42)Rгде индекс sw подчеркивает, что в треугольной матрице R заполненнойчастью может быть только «юго-западная» (south-by-west) часть.
Вполнеочевидно, что эти варианты получаются простым изменение порядка действий в алгоритмах преобразований.7.9Правосторонние ортогональные преобразованияС правосторонними ортогональными преобразованиями мы уже сталкивались (см.
подразд. 7.6); тогда для квадратной матрицы A после T A = Rвычисляли A−1 = R−1 T . Однако, можно начинать с правостороннего преобразования матрицы A; тогда отыскание A−1 потребует, соответственно,левостороннего преобразования.Пусть A = A(n, n) — квадратная невырожденная матрица. Будем рассматривать ее строки как векторы в Rn . Преобразования вектора как матрицы-строки в n-мерном линейном пространстве задаются умножением ее напреобразующую матрицу справа.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
