Rambler's Top100 Service
Поиск   
 
Обратите внимание!   Посмотрите новые поступления ... Обратите внимание!
 
  Наука >> Математика >> Алгебра, математическая логика и теория чисел | Популярные статьи
 Написать комментарий  Добавить новое сообщение
 См. также

Популярные статьиОтветы на все вопросы

Линейная алгебра: от Гаусса до суперкомпьютеров будущего

В.П.Ильин, доктор физико-математических наук, профессор,
главный научный сотрудник Института вычислительной математики и математической геофизки СО РАН.
Опубликовано в журнале "Природа", N 6, 1999 г.
Содержание

Об авторе

Валерий Павлович Ильин, доктор физико-математических наук, профессор, главный научный сотрудник Института вычислительной математики и математической геофизки СО РАН. Область научных интересов - вычислительная математика, информатика, математическое моделирование. Автор более 250 научных статей и девяти монографий (часть из них написана в соавторстве).

...Поверил
Я алгеброй гармонию
А.С.Пушкин

Введение: история, предпосылки и проблемы

Задача численного решения систем линейных алгебраических уравнений (СЛАУ) имеет незапамятную историю. Классический метод исключения, активно развиваемый и изучаемый даже в наши дни, был открыт К.Гауссом в 1849 г. Однако еще за 2000 лет до этого в Древнем Китае изданы "Девять книг о математическом искусстве", где этот алгоритм уже изложен в характерной для своего времени "натуральной" форме, но фактически с использованием матричных преобразований.

Вплоть до начала XX в. алгебра оставалась "наукой о решении уравнений", после чего произошло ее разделение на высшую алгебру (операции с абстрактными объектами различной природы) и линейную, основа которой - матричное исчисление.

Становление современных вычислительных методов линейной алгебры можно считать состоявшимся после выхода в 1960 г. книги Д.К. и В.Н. Фаддеевых с одноименным названием, по которой училось не одно поколение российских и зарубежных математиков. Необходимо подчеркнуть, что актуальные проблемы вычислительной алгебры имеют фундаментальный характер не только потому, что текущая компьютеризация различных областей знаний в значительной степени сводится к векторно-матричным процедурам. Изучение матриц, являющихся операторами в простейших конечномерных пространствах, позволяет обнаружить наиболее глубокие и тонкие свойства математических объeктов, имеющих свое значение для функционального анализа, теории аппроксимации, дифференциальных уравнений и так далее

Последние десятилетия ознаменованы бурным развитием численных методов, нашедших отражение в книгах Г.И.Марчука, А.А.Самарского, С.К.Годунова, Дж.Голуба, О.Аксельсона и других. В рамках итерационных методов решения СЛАУ значительные продвижения достигнуты в таких современных направлениях, как алгоритмы сопряженных направлений, последовательной и симметричной верхней релаксации, неявные алгоритмы переменных направлений и методы неполной факторизации.

Несмотря на имеющиеся в этой области теоретические достижения, а также на бурный рост мировых вычислительных мощностей, проблема конструирования и исследования новых быстрых алгоритмов решения СЛАУ не теряет своей актуальности, и поток публикаций по этой тематике только увеличивается. Дело объясняется тем, что возникающие практические задачи наряду с усложнением вычислительных процедур требуют от методов решения все большей точности и надежности. Как говорится, аппетит приходит во время еды, и с этим связан тот уникальный феномен, что, невзирая на появление все новых поколений ЭВМ, за последние 50 лет решение "больших" (на текущий момент) задач всегда занимало примерно одинаковое астрономическое время.

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

Главный источник сверхбольших СЛАУ - сеточные методы (конечных разностей, конечных объемов и конечных элементов) решения многомерных краевых задач, возникающих при моделировании процессов или явлений. Для оценки потребностей вычислительных ресурсов рассмотрим сеточную трехмерную расчетную область в форме параллелепипеда, образованную координатными линиями
$x=x_i, i=1,...,I; y=y_j, j=1,...,J; z=z_k, k=1,...,K$ Если мы решаем одно дифференциальное уравнение второго порядка эллиптического типа (например, уравнение стационарной диффузии или теплопроводности), то в результате его простейшей аппроксимации для каждого узла (i,j,k) получается уравнение вида
$(Au)_{i,j,k} = p_0 u_{i,j,k} - p_1u_{i-1,j,k} - p_2u_{i,j-1,k} - p_3u_{i+1,j,k} - p_4u_{i,j+1,k} - p_5u_{i,j,k-1} - p_6u_{i,j,k+1} = f_{i,j,k'}$, (1)
где коэффициенты $p_k$ выражаются через шаги сетки и параметры исходной задачи. Такой системе уравнений отвечает квадратная матрица А порядка IJK с числом ненулевых элементов на каждой строке не более 7.

Пусть ради простоты $I=J=K \gg 1$. Тогда для реализации метода исключения Гаусса требуется выполнить примерно $Q=K^7$ арифметических операций и запомнить $P=K^5$ промежуточных чисел. На сегодня задаче средней сложности соответствует K=100, что приводит к значениям $Q=10^{14}$ и $P=10^{10}$. На обиходном языке это означает около трех часов работы компьютера быстродействием 10 гигафлоп, то есть $10^{10}$ "флоп" - арифметических операций над "нормальными" вещественными числами в секунду, при наличии оперативной памяти около 10 гигабайт. Это намного превосходит ресурсы персональных компьютеров и рабочих станций, но на мировом рынке такая ЭВМ - среднего класса. Однако если перейти к "большой" задаче с K=1000, то она уже оказывается крепким орешком даже для рекордного суперкомпьютера со скоростью в один терафлоп ($10^{12}$) и памятью в один терабайт.

Двумерные задачи гораздо легче. Однако и они требуют огромных ресурсов, если нужно находить оптимальные параметры искомых решений (чтобы удовлетворить поставленным ограничениям и выбранным критериям качества), рассчитывать длительные по времени или, наоборот, быстроменяющиеся процессы, учитывать сильнонелинейные эффекты. Это приводит к многократному решению промежуточных линейных систем. Важно отметить, что это - наиболее трудоемкая часть крупномасштабного вычислительного эксперимента, так как остальные расчетные этапы требуют объема арифметических действий, пропорционального общему количеству узлов сетки.

Кардинальное средство повышения производительности компьютеров (по профессиональным прогнозам - до петафлопного быстродействия, т.е. до $106{15}$ операций в секунду к 2015 г.) - создание многопроцессорных вычислительных систем (до $10^3$ или $10^6$ процессоров). Отсюда одной из первых по актуальности проблем в вычислительной алгебре становится распараллеливание алгоритмов и их отображение на архитектуру многопроцессорных ЭВМ. Насущность этого вопроса, в том числе при решении СЛАУ, проистекает из того общеизвестного факта, что эффективность загрузки всех устройств - это больное место различных суперкомпьютеров, и реализация алгоритмов без учета, например, коммуникационных потерь может снизить фактическую производительность многопроцессорной системы до нескольких процентов от ее паспортных данных. Один из самых эффективных методов решения СЛАУ высокого порядка базируется на итерационных алгоритмах неполной факторизации, получивших сейчас самое широкое распространение в мировой вычислительной практике1. Данная статья представляет последние результаты автора, касающиеся новых, перспективных подходов: обобщенного принципа компенсации в приближенной матричной факторизации, адаптивных итерационных методов с оптимальными упорядоченностями и, наконец, проблемы распараллеливания алгоритмов на основе декомпозиции областей.

Совсем не алгебраическая, но очень плодотворная идея

Идея метода неполной факторизации для решения алгебраической системы
$Au=f$ (2)
заключается в построении такой матрицы B, называемой предобусловливающей, которая легко бы обращалась и была бы близкой к исходной матрице A в том смысле, что собственные числа матрицы-произведения B-1A не сильно различаются между собой. Тогда последовательные приближения un строятся с помощью итерационного процесса, который в канонической записи имеет следующий вид:
$B(u^n-u^{n-1})= \omega_n (f-Au^{n-1})$ (3)
Здесь $\omega_n$ - вычисляемые параметры, а число итераций $n(\varepsilon)$ для достижения требуемой точности $\varepsilon \ll 1$ зависит от отношения
$k=max \mid \lambda \mid /min \mid \lambda \mid $,
называемого числом обусловленности матрицы B-1A. Эта величина k оказывается магическим числом в линейной алгебре, характеризуя главные вычислительные качества любой матрицы, в том числе чувствительность к погрешностям округлений.

Вводя обозначения
$v^n=u-u^n$, $r^n=f-Au^n$
для векторов коррекции и невязки, соотношение (3) можно разбить на два:
$Bv^n = r^{n-1}$, $u^n = u^{n-1}+\omega_n v^n$,
из которых видно, что наиболее трудоемкая процедура на каждой итерации - решение вспомогательной системы с матрицей B (конечно, если последняя не слишком тривиальна). Первооткрыватель методов неполной факторизации Н.И.Булеев предложил (середина 50-х годов, закрытый ядерный центр в г.Обнинске) искать B как приближенную (неполную) факторизацию исходной матрицы A:
$B = LU \cong A$, (4)
где L, U - легко вычисляемые (намного проще, чем в точной факторизации) и легко обратимые нижняя и верхняя треугольные матрицы.

Гениальное все просто: помимо конструктивного построения множителей L, U, Булеев предложил в качестве критерия близости матриц требование
$Be = Ae$, (5)
где e - вектор с постоянными (единичными) компонентами; он назвал его принципом компенсации. Эта совсем не алгебраическая идея исходила из тех сермяжных соображений, что в задачах матфизики решения в основном гладкие, а тогда при условии (5) на каждой итерации "лишние" члены (обусловленные неточностью факторизации (4)) будут взаимно компенсироваться.

Новые методы на многих практических задачах обнаруживали потрясающее быстродействие, но временами не все было гладко. Это говорило о необходимости более глубокого изучения нового класса алгоритмов, теории которых тогда не было. Но, как известно, нет пророка в своем отечестве - ведущие специалисты в СССР идею неполной факторизации не подхватили, продолжая увлеченно развивать другие (пришедшие из-за рубежа) модные направления.

На Западе же методы Булеева были несколько раз переоткрыты и стали активно усовершенствоваться. В частности, эвристическое условие (5), которое на алгебраическом языке не что иное, как принцип согласования строчных сумм матриц B и A, было объяснено и эффективно использовано в сочетании с ускорением итераций универсальными методами сопряженных градиентов. Поразительное совпадение, но практически одновременно (в 1983 г.) и независимо в работах автора с О.Аксельсоном (Голландия) и Дж.Голуба (США) с коллегами были предложены неявные (блочные) методы неполной факторизации, в которых матричные множители L, U - суть блочно-треугольные матрицы, причем диагональные блоки - ленточные матрицы, ненулевые элементы которых сосредоточены только в окрестности ленты заданной ширины d около главной диагонали. Такие алгоритмы позволили уточнить саму приближенную факторизацию и значительно сократить число итераций. Однако здесь возникла новая нетрадиционная задача вычислительной алгебры: найти ленточную часть матрицы, обратной к ленточной же матрице. Для этой задачи удалось найти изящное решение (естественно, трудоемкость процедуры увеличивается с ростом ширины ленты d).

В этом месте заложена принципиальная методологическая проблема, ожидающая до сих пор своего решения. Пусть, например, C - невырожденная матрица. Если мы вычисляем ленточную часть ширины $d \gt 1$ от обратной матрицы - обозначим ее через $(C^{-1})^{(d)}$, - то находим в определенном смысле приближение к обратной матрице $(C^{-1})$ (которая в общем случае плотная, т.е. не содержит нулей). Действительно, если d увеличивается вплоть до порядка матрицы N, то $(C^{-1})^{(d)} \rightarrow (C^{-1})$ При d=N матричная аппроксимация становится точной и метод решения СЛАУ из итерационного превращается в прямой, дающий точное решение за конечное число действий. Фактически в этом случае метод приближенной факторизации переходит в блочный метод Гаусса, который для сеточных систем уравнений слишком дорог.

Здесь и возникает вопрос: какова оптимальная ширина$1 \leq\le d \leq\le N$ ленточной аппроксимации обратной матрицы $(C^{-1})^{(d)}$, чтобы обеспечить итоговую эффективность итерационного метода? Или, другими словами, достичь заданной точности за минимальное число арифметических действий? Понятно, что лучшее - враг хорошего, и в полном объеме решение проблемы оптимизации алгоритмов - это недостижимая голубая мечта, к которой надо стремиться, но только дробными шагами, расширяя жизненное пространство конкретными результатами для отдельных типов систем.


1Il'in V.P. Iterative incomplete factorization methods. Singapore, 1992; Ильин В.П. Методы неполной факторизации для решения алгебраических систем. М., 1995.

Назад | Вперед


Написать комментарий
 Copyright © 2000-2015, РОО "Мир Науки и Культуры". ISSN 1684-9876 Rambler's Top100 Яндекс цитирования