Раздел: Документация

0 ... 16 17 18 19 20 21 22 ... 78

нужных собственных векторов априорно известна; это большое достоинство. Многие химики работают на компьютерных системах средней производительности, и их трудности связаны, скорее, с особенностями операционной системы, чем с методом вычислений.

Программное обеспечение, разработанное большинством специалистов по собственным значениям, встроено в их специализированные пакеты программ, и трудно понять, как они могут воспользоваться более совершенной программой вычисления собственных значений больших матриц, разработанной для универсальных ЭВМ. Им придется, конечно, в этом подробно разобраться, а затем ввести эти программы в свои системы. Однако эту работу следует предпринимать только в том случае, если будут видны перспективы в улучшении на порядок надежности, быстродействия или требований к объему памяти.

6.8. МЕТОД ИТЕРИРОВАНИЯ ПОДПГОСТРАНСТВА ДЛЯ РЕШЕНИЯ УРАВНЕНИЯ (К - ХМ)z=0

Этот метод представляет собой любопытную версию блочного метода обратной итерации с перемежающимися сдвигами (более детально см. [12, гл. 14, 15 J). Пользователь должен выбрать множество изр ортонормирован-ных начальных векторов. Пусть это будут столбцы некоторой матрицы X размера пХр. Выбрать р не так легко, но пусть р<5(). Подпространство, о котором идет речь в названии раздела, является подпространством столбцов матрицы X. Внутренний цикл алгоритма SI, по существу, состоит в следующем:

1)решить матричное уравнение (К оМ) Y = MX дли (пХр) -матрицы Y;

2)сформировать (рХр)-матрицы К = YrKY и М = YMY;

3)определить все собственные значения малых матриц (К - аМ) G = MGD, где матрица D является диагональной и G - матрица размера рХр, составленная из собственных векторов;

4)положить X = YG;

5)проверить на сходимость элементы матрицы D.

Существует ряд существенных видоизменений этой схемы, Но они не важны для данного обзора (более подробно см. [1 и 12, гл. 14]).

Если элементы матрицы D сходятся к предельным значениям, соответствующие векторы матрицы X являются удовлетворительной аппроксимацией собственных векторов.

Типичными значениями являются п = 500, р = 20. Иногда значение р достаточно велико и приходится использовать оперативное ЗУ на первом шаге. Если р выбрано слишком малым, то число шагов для достижения сходимости резко возрастает.

В алгоритме SI в неявном виде предусмотрена обработка матрицы (К - оМ)~М, которая является несимметричной, но самосопряженной со скалярным произведением с матрицей М, что является достаточным для того, чтобы собственные значения были действительными. (Скалярное произведение с матрицей М определяется как (х, у) =yfMx.)

69 АЛГОРИТМ ЛАНЦОША [12, гл. .13]

Пусть А = Аг г - произвольный начальный вектор, q0 -0, Pi -11 II- Приведем внутренний цикл базисного алгоритма Ланцоша: для к = 1,2,.. - до сходимости, выполнить

1.q* «- «#*

2.r<-/lq*.

3.r-r-qpV [ 4. Занести qbl в оперативное ЗУ.

5.ак+-цкг.

6.г <- г - qk ак. 7- ft + ,*-r. 8. Проверку на сходимость.

Заметим основную особенность этого цикла. Он изобилует векторными операцГми, линейными комбинациями и скалярными произведениями ESTI и Р рассматриваются как элементы трехдиагональнои матрицы

Р2

Рк

которая с каждым шагом увеличивается на одну строку и один столбец. Пусть Q#~ (qi, q2, • - - , q#). Наиболее важным является соотношение

aq* - q А = о = рк+1чк+А + F*>

где Fk связано с погрешностями округления и всегда мало. Здесь вектор е£= (0, 0, - . . , 0,1) содержит к элементов; все другие векторы имеют длину п. При абсолютно точных вычислениях

<4Q* = i*

но в действительности такою не бывает. Тем не менее это несоответствие только замедляет выполнение алгоритма, но не препятствует сходимости.

На самом деле происходит то, что некоторые собственные значения матрицы Т,- при увеличении к стремятся к собственным значениям матрицы А. Если & и s представляют собой собственную пару матрицы Тк, т. е.

Tfcs=s#,s = 1,

119

---

то аппроксимируемый собственный вектор, соответствующий есть у = = Qks.

Существует граница погрешности вычислений д. Пусть s (/с) - последний элемент s. Существует собственное значение матрицы А, удовлетворяющее условию

Ну IMlVll-

В большинстве прикладных задач величина ,б/с никогда ие бывает малой, так что сходимость обусловлена уменъшешем is (к) по меря того, как к увеличивается. Для вычисления s (к) необходимо только найти собственный вектор трехдиш оналыюй (кЛк) -матрицы. На самом деле существуют методы даже более быстрого нахождения s (к). Здесь мы опять опустили многочисленные детали, чтобы выделить основные особенности.

Матрица А часто совпадает с матрицей, описанной в разд. 6.4, поэтому второй шаг (г *-Аа.) в основном определяет стоимость вычислений по алгоритму Ланцоша. Достоинство этого алгоритма состоит в том, что в отличие от алгоритма итерирования подпространства, при котором на каждом шаге вычислений одна матрица X заменяется на другую, здесь информация не отбрасывается. Теория сходимости алгоритма Ланцоша хорошо разработана, и крайние собственные значения матрицы А часто сходятся после выполнения примерно 2\[п шагов.

При реализации алгоритма Ланцоша в сверхоперативной памяти ЭВМ требуется хранить только три или четыре «-мерных вектора совместно с данными, необходимыми для выполнения процедуры г -<-Aq. Однако при необходимости нахождения собственных векторов векторы qi, 42 > • • • - Ча должны быть извлечены из оперативного ЗУ для образования вектора

У = QfcS = q,s(l) + q2s(2) +• •• + qks(k).

В настоящее время продолжаются работы, направленные на превращение алгоритма Ланцоша в надежную программу. В то же время известна хорошая документированная программа, носящая имя LAS02, разработанная Д. Скоттом в Национальном энергетическом центре программного обеспечения Национальной лаборатории Аргонны.

6.10. СИНГУЛЯРНЫЕ ЗНАЧЕНИЯ (ДЛЯ МАЛЫХ МАТРИЦ)

В большинстве вычислительных центров имеются очень хорошие программы для вычисления малых матриц. Наиболее известная из них находится в библиотеке программ UNPACK (но не EISPACK). Сингулярное разложение действительной матрицы В размера тХп есть В = UXV, где матрица U размера тХг и матрица V размера пХг являются ортонормированными, 120

г. е. UU = VV = I,. и I = diag(a,,ог), где оу >0. Здесь г - ранг матрицы В. Ясно, что

I BfB = V£2Vf,

откуда следует, что £с2]- - ненулевые собственные значения матрицы ВГВ, столбцы матрицы V являются соответствующими собственными векторами.

Нет необходимости записывать матрицу ВВ в явном виде. Программа LINPACK эквивалентна QR-методу Хаусхолдера применительно к матрице BfB, но фактически явно преобразуются только элементы матрицы В [4, гл. 11 ]. При этом программа выполняется в два этапа.

I 1. Последовательность п - 1 простых двусторонних ортогональных преобразований

B-»P*BQ

используется для того, чтобы преобразовать матрицу В к диагональному виду

2. Другая последовательность матриц Р и Q используется для уменьшения ненулевых внедиагональных элементов без искажения каких-либо нулевых элементов. Сингулярные векторы (так названы столбцы матриц U и V) определяются путем накопления произведений матриц Р и Q по мере их использования. Это устойчивый компактный алгоритм, реализуемый с гарантией. Сингулярное разложение дает наиболее информативный и надежный метод решения линейных задач по методу наименьших квадратов, но в то же время это наиболее дорогостоящий метод, мощность которого во много раз превышает потребности ряда прикладных задач. Если т>п, то число операций умножения, которое требуется для вычисления матриц О и V, а также расчета £, приблизительно равно

2т* (и - т(3) + 5т2п;

при п-т число этих операций составляет примерно 6л3 (что эквивалентно шести операциям умножения матрицы на матрицу).

Заметим, что матрицы U и V образованы с помощью операций вида М<~ UP, V <- VQ, где Р и Q практически ортогональны. Следовательно, матрицы U и V должны быть практически ортонормированы, а поэтому левые и правые сингулярные векторы должны быть ортогональны, как бы близко сингулярные значения ни располагались.

---

6.11. СИНГУЛЯРНЫЕ ЗНАЧЕНИЯ (ДЛЯ БОЛЬШИХ МАТРИЦ) Эта проблема по-прежнему пока что изучается. Не существует доступного стандартного программного обеспечения. Задачи, решаемые методом наименьших квадратов, характеризуются значениями п ~300ООО и т= 2000000, но при этом всего несколькими ненулевыми элементами в матрице В [10]. Для решения задач методом наименьших квадратов полное сингулярное разложение производить не надо. Поэтому используются другие, более дешевые методы. Один из таких методов состоит в вычислении QR-разложения (не следует путать с QR-алгоритмом) с перестановкой столбцов матрицы. В том случае, когда матрица имеет полный ранг, наиболее просто использовать множитель Холецкого С матрицы ВВ (т.е. В*В = СГС, где С - верхняя треугольная матрица размера пХп) для решения нормальных уравнений. При абсолютно точном вычислении матрица ВС-1 является ортогональной и ее столбцы дают ортонормированный базис пространства столбцов матрицы В.

Существует несколько хороших параллельных алгоритмов для вычисления множителя Холецкого положительно определенных матриц [21]. Однако будут ли они пригодны в случае, когда и 300000, не ясно! С возможностями применения параллельных алгоритмов для малых значений п можно познакомиться в работе [24].

6.12. О ПРОБЛЕМЕ СОБСТВЕННЫХ ЗНАЧЕНИЙ ДЛЯ БОЛЬШИХ НЕСИММЕТРИЧНЫХ МАТРИЦ

В настоящее время нет эффективного, надежного метода для вьгделения некоторых или вычисления всех собственных значений больших разреженных несимметричных матриц и очень мало хороших доступных средств программирования. Это является предметом активных исследований [13 - 15], но экономические побудительные мотивы здесь кажутся не столь важными, как в случае симметричных задач. Обзор этих исследований можно найти в работе [18]. Особо важный специальный случай, который оправдывает внимание к несимметричным задачам, связан с наличием возмущений в симметричных задачах.

6.13. ВОПРОСЫ НА БУДУЩЕЕ

Не удивительно, что в методах, которые оказались наиболее удачными для применения в последовательных компьютерах, в полной мере учитываются основные особенности этих компьютеров: легкий доступ процессора к ЗУПВ большой емкости и легкость реализации сложных циклов в программах. Было бы удивительным, если бы эти методы были пригодны и для параллельных процессоров.

Кажется, что все исследования в области параллельных алгоритмов сфокусированы на том, чтобы превзойти возможности последовательных ЭВМ при решении стандартных малых задач (при и<100), таких как умножение матрицы на матрицу или решение линейных уравнений. Однако именно эти 122

малые задачи не интересуют специалистов по вычислительным методам, которые считают, что все их требования в области малых задач удовлетворены. Их исследования направлены на решение больших задач: либо на совершенствование общих методов, либо на разработку специальных приемов для наиболее важных частных случаев, таких как решение квадратичных X матриц (л2А + лВ + С)х = 0, или на использование суперкомпьютеров типа Сгау-1 вместо последовательных компьютеров.

Все, что я читал о методах вычисления собственных значений с помощью параллельных алгоритмов, посвящено обработке заполненных матриц (см., например, [8]). Однако большинство больших матриц являются разреженными и мог бы оказаться пригодным любой алгоритм, учитывающий разреженность, если стараться сохранить емкость оперативной памяти.

Задачи обработки изображений и распознавания образов - это задачи обработки заполненных матриц. Рассмотрим кратко задачу вычисления 50 наибольших сингулярных значений (и векторов) матрицы В размера 1000Х Х500 элементов. Данные содержат 106/2 чисел. (Если эти числа известны только с точностью до одной десятичной цифры, то кажется не совсем удачным тратить 106/2 машинных слов для их запоминания.) Применение общих приемов, таких как использование нормальных уравнений или параллельных алгоритмов [7, 23], является "стрельбой из пушки по воробьям". Согласно им выполняется вся работа, необходимая для нахождения всех 500 сингулярных значений. Возможно ли здесь использовать априорные знания? Важной помощью было бы указание 10 или 20 столбцов матрицы В с наибольшей нормой.

Это позволило бы в последовательных ЭВМ использовать блочный алгоритм Панцоша с порядком блока 10 или 20 для оператора ВВ. Нет необходимости образовывать произведение матриц BfB для вычисления матричного произведения W (ВХ). Здесь достаточно 20 или 10 итераций для сходимости 50 наибольших собственных значений матричной формы ВВ, а число арифметических операций примерно равно тому, которое необходимо для образования ВВ. Каким образом в параллельных алгоритмах можно эффективно учесть изменение числа необходимых сингулярных значений? Очень малые изменения в приближениях, соответствующих большой задаче, могут сильно повлиять на эффективность других методов.

Здесь есть одна трудность. Как могут систолические алгоритмы с очень однородным потоком информации эффективно адаптироваться к изменениям в конкретной задаче или имеющейся информации? Если они не смогут адаптироваться, то могут ли они быть эффективными рдя больших задач?

СПИСОК ЛИТЕРАТУРЫ

[1] K.-J. Bathe, and S. Ramaswarmy, "An Accelerated Subspace Iteration Method," Comput. Methods Appl. Mech. Eng., 25:313-331 (1980).

[2] K.-J. Bathe, and E. Wilson, Numerical Methods in Finite Element Analysis, Prentice-Hall, Englewood Cliffs, NJ, 1976.

0 ... 16 17 18 19 20 21 22 ... 78