Форсайт Дж., Малькольм М., Моулер К. - Машинные методы математических вычислений (1040536), страница 42
Текст из файла (страница 42)
Преобразование порождается первой строкой А,: огг (О 50 26! 27 642) 1 !' 71 =. Р, = 1 — У, о,ог, -1 1' Аг = А,У1, Строка матрицы А, с номером 1 получается по формуле т11, г (1,-1,,г! г ( — гого) г Таким образом, из каждой строки А, вычитается надлежащее кратное огг. Зто дает матрицу — 7.416 32.732 0 0 -3.133 0.3 19 0 0.861 — 0.088 0 4.856 -0.495 0 8.850 ' — 0.901 0 10;605 -1.080 0 0 0.000 0 0 0.000 0 0 0.000 Аг — ЦЛАг- Поскольку первая компонента о,' нулевая, то нули первого столбца А, сохраняются в А,. Так как 1/, ортогональная, то длина каждой строки в А, равна длине соответствующей строки в Аг Теперь можно добиться новых нулей под диагональю, не испортив полученных ранее: -7.416 32.732 0 9.е Вычисление РАзлОжения 245 Поскольку ранг этой матрицы равен лишь 2, то теперь третий столбец имеет на диагонали и под диагональю элементы порядка опц~бки округления.
Эти элементы обозначены в матрице через 0.000, чтобы отличить их от элементов, в точности равных нулю. Если бы матрица имела полный ранг, то нужно было бы выполнить еще одно преобразование, чтобы получить нули в третьем столбце: -7А!6 32.732 0 О 10.605 1.080 0 0 0.000 0 0 0 0 0 0 В данном примере третий диагональный элемент был бы точным нулем, если бы не ошибки округлений. Элементы под диагональю во всех столбцах указаны как точные нули, потому что преобразования так и строились, чтобы получить там нули.
Последнее преобразование ста в этом примере могло бы быть тождественным, однако тогда оио ие будет хаусхолдеровым отражением. Фактически использованное ст', попутно изменяет знак элемеита— 1.080 в матрице А,. Получилась искомая двухдиагональная матрица, и первый этап закончен. Прямое использование ортогональных преобразований не позволяет получить какие-либо новые нули. Для общего порядка п нужно гг преобразований с/ и п — 2 преобразований )У, чтобы достигнуть этого места.
Число преобразований не зависит от строчной размерности пг, но от гл зависит работа, затрачиваемая на выполнение каждого преобразования. Второй этап алгоритма О'ег0 представляет собой итерационный процесс. Каждый шаг уменьшает величину изгоняемых наддиагональных элементов. В конечном счете этп элементы становятся сравнимыми с ошибками округленяй, и их без ущерба можно отбросить. Фактическое число шагов зависит от конкретной матрицы тл точности машины, но снорость сходимости очень высокая, так что обычно на выполнение второго этапа требуется меньше времени, чем на выполнение первого.
Это особенно верно, если гп много больше, чем п, поскольку на втором этапе :охраняются все нули, полученные ниже первых и строк. Второй этап есть вариант 4;)К-алгоритма, разработанного Фрэнсисом ') около 1060 г. Этот метод стал одним из наиболее надежных и широко используемых численных алгоритмов. Полное ') г2п-алгоритм (под оазаанием метода односторонних прагдений) бмл однопременно и пезаеисимо предложен В. В. Кублапоеской.— Прим, перга. В НАИМЕНЬШИЕ КВАДРАТЫ 246 обсуждение его деталей и, в частности, удовлетворительное объяснение того, почему он работает, выходят за рамки этой книги.
Мы лишь кратко обрисуем поведение алгоритма для рассматриваемого примера. Интересующемуся читателю советуем за деталями обратиться к указываемой ниже литературе. Каждый итерационный шаг начинается с двухдиагональиой матрицы, которая преобразуется в другую двухдиагональную матрицу с меньшими внедиагональными элементами. Мы подробно рассмотрим один такой шаг в нашем примере. Поскольку две последние строки нулевые и остаются нулевыми в дальнейшем, то их можно опустить. Будем указывать для матричных элементов пять десятичных знаков после запятой вместо трех, как было до сих пор, Начинаем с матрицы — 7.41620 32.73169 0 0 10.60517 1,08016 0 0 0.00000.
Исходя из трех элементов правого нижнего угла, вычисляется некоторое преобразование Ъ'; когда А замещается матрицей А У, изменяются два первых столбца: 33.56105 0.13897 0 10.35262 — 2.30062 1.08016 0 0 0.00000. Заметим, что под диагональю появился ненулевой элемент 10.35202. Ь!тобы избавиться от этого нежелательного элемента, вычисляется преобразование 17; при переходе от А к !7А изменяются две первые строки: 35.12152 — 0.54535 0.31839 0 — 2.23937 1.03216 0 0 0.00000.
Требуемый нуль восстановлен, однако выше наддиагонали воз-. ник новый ненулевой элемент 0.3!839. Следующее преобразова' ние Р' комбинирует два последних столбца, чтобы исключ этот ненужный элемент: 35.! 2! 52 0.63149 0 0 2.45431 0.23771 0 0.00000 0.00000. Если бы правый нижний элемент не имел значение 0.00000 под диагональю появился бы еще один ненулевой элемент, т где в нашей матрице стоит второй указатель 0.00000.
Следую з.с Вычнслгнип Рллложгния преобразование У, которое в данном примере мало отличается от единичной матрицы,'оперирует с двумя последними строками '7: ' 35.12152 0.63149 0 0 2.45431 0.23771 0 0 0.00000. Мы вернулись к двухдиагональной матрице. Благодаря первоначальному 1г-преобразованию оба падднагопальпых элемента уменьшились. (Для других матриц некоторые надднагопальные элементы могут на отдельных шагах даже возрастать, но в конечном счете онн уменьшаются.) Процесс повторяется. Следующей будет получена двухдиагональная матрица 35.12722 0.00309 0 О 2.46540 — 0.00014 О 0 0.00000.
Заметим, что диагональные элементы изменились мало и что оба наддиагональпых элемента уменыпились. После еще одного шага получим 35.12722 0.00002 0 0 2.46540 0.00000 0 0 0.00000 Изменения диагональных элементов происходят лишь в шестом или седьмом десятичном знаке после запятой, которых мы уже не показываем. Последний наддиагональный элемент имеет теперь порядок ошибки округления и может рассматриваться как нуль. Заключительный шаг уменьшает до пренебрежимо малого значения другой надднагональный элемент: 35.12722 0 0 0 2.46540 О 0 0 0.00000. Три диагональных элемента и есть сингулярные числа.
Мы не указывали явно использованные преобразования, а лишь результаты их воздействия на матрицу А; дело в том, что процесс выполняется посредством операций со строками и столбцами А, а не фактических матричных умножений. Если нужно находить только сингулярные числа, то А — единственный дву- '1 Смысл этого нреобраэоаання — ию<лгояение элемента (3,21,— Прим. лерга. к нлимепьшие квйдялты 'мерный массив, участвующий в вычислениях.
(В самом деле, в подпрограмме ЬНО двухдиагональная матрица представляется двумя одномерными массивами, один нз которых на выходе содержит сингулярные числа.) Если нужно вычислять не только сингулярные числа, по и матрицы (7 и У сингулярных векторов, то запасаются еще два массива, начальным состоянием которых служат единичные матрицы. В дальнейшем всякий раз, когда выполняется операция над строками или столбцами А, та же операция производится с одним из двух других массивов. Когда А в конечном счете приобретает диагональный вид, в остальных массивах содержатся матрицы (7 и К.
Зтот алгоритм типичен для используемых в современных мап|инных подпрограммах решения различных матричных задач на собственные значения. Если задача не находигся уже в специальной форме, то сначала матрица переводится в другую, содержащую много нулей. Затем к приведенной матрице применяется итерационный процесс, позволяющий получить еще больше нулей. Обычно повсюду используются ортогональные преобразования.
Первоначальное приведение не только обеспечивает эффективность итерационного процесса; оно также ускоряет его сходимость. Для алгоритма ЬНО, например, приведенная матрица двухдиагональная, итерационный процесс сконструирован так, чтобы аннулировать наддиагональные элементы, и теорема Уилкинсона утверждает, что итерации всегда сходятся н сходимость очень быстрая. Теорема Уилкинсона обсуждается в применении к ЬНЭ в книге Лоусон, Хансов (!974). 9.5. Подпрограмма 5ЧЬ АЛГОЛ-процедура для вычисления сингулярного разложения (авторы Голуб и Райнш) была опубликована в собрании; матричных алгоритмов под редакцией Уилкинсона и Райнша.
' В рамках проекта МАТЬ в Национальных лабораториях в Аргон- не были разработаны фортран-подпрограммы для ряда матричных задач на собственные значения, в том числе трансляция процедуры ЬН1). КАТЬ вЂ” это сокращение от Ка1!опа1 Ас11гд1 1о Тез( Ьо((маге (Национальная программа тестирования мате матического обеспечения); соответствующая библиотека по программ получила название Е1ЬРАСК.
Была опубликована ии струкция для пользователя, относящаяся к части второй очер ди Е1ЬРАСК (Смит и др. (1976)). Хотя ЬНР включена во втору очередь, она не описывается в инструкции издания 1976 г: а будет описана лишь в последующем издании. у а пОдпРОГРАммА Бго 243 Наш вариант 5ЧО отличается от варианта из Е)5РАСК тестом на малость элементов. Вариант из Е)5РАСК требует включения оператора РАТА, содержащего машинное эпсилон.
Мы исключили машинное эпсилон и заменили тест на малость элементов оиераторами типа !Р(АВ5(Р)+АР(ОКМ .ЕО, Ай)ОКМ)ОО ТО 565. В подпрограмме имеются три таких теста. Мы решились на это изменение с некоторой неохотой, потому что подпрограммы, вошедшие в Е!5РАСК, широко и тщательно тестировались, и распространяются и гарантируются проектом Р)АТ5. Однако внесенное нзменение очень незначительно, и мы полагаем, что достигнутая этим машинная независимость стоила такой переделки. Входом в ЯЧ) являются тки-матрица А, информация о размерностях и две логические переменные, указывающие, нужно ли вычислять матрицы (7 и )г.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
