Universal algorithms for solving discrete stationary Bellman equations

Full Text

Введение Вычислительные алгоритмы строятся на основе некоторых простых операций. Эти операции манипулируют данными, которые описывают ѕчислаї. Эти ѕчислаї являются элементами ѕчислового доменаї (numerical domain), то есть такого математического объекта, как поле вещественных чисел, кольца целых чисел, различные полукольца и т. д. В практических задачах числовые домены заменяются их компьютерными пред- ставлениями, т. е. элементами конкретных конечных моделей этих доменов. Примерами моделей, которые удобно использовать для компьютерного представления веществен- ных чисел, являются различные модификации арифметики с плавающей запятой, при- ближенная арифметика рациональных чисел [37], интервальная арифметика и другие. Отличие между математическими объектами (ѕидеальнымиї числами) и их конечны- ми моделями (машинными представлениями) выражается в вычислительных ошибках (например, при округлении). Алгоритм называется универсальным, если он не зависит от конкретного числового домена и/или его компьютерного представления [32, 33, 36, 40]. Типичным примером уни- версального алгоритма является вычисление скалярного произведения (x; y) векторов x = (x1; : : : ; xn) и y = (y1; : : : ; yn) по формуле (x; y) = x1y1 +: : :+xnyn: Этот алгоритм (формула) не зависит от конкретного домена и его машинной реализации; формула имеет смысл для любого полукольца. Очевидно, что один алгоритм может быть более универсальным, чем другой. Например, простейшая формула Ньютона-Котеса (форму- ла прямоугольников) является самым универсальным алгоритмом численного интегри- рования. В частности, эта формула верна также для идемпотентного интегрирования (т. е. над любым идемпотентным полукольцом, см., например, [6, 31]). Другие формулы интегрирования (например, формула трапеций или формула Симпсона) не зависят от машинной арифметики и могут использоваться (например, в итерационной форме) для вычислений с произвольной точностью. С другой стороны, алгоритмы, основанные на формулах Гаусса-Якоби, созданы для вычислений с заданной точностью: они содержат константы (коэффициенты этих формул), определенные с конечной точностью (конеч- но, алгоритмы такого типа можно сделать более универсальными, включив процедуры вычисления констант, однако, это влечет усложнение алгоритмов). Современные достижения в математике и разработке программного обеспечения позволяют рассматривать численные алгоритмы и их классификацию с новой точки зрения. Обычные численные алгоритмы ориентированы на программную (или аппарат- ную) реализацию на основе арифметики с плавающей или фиксированной запятой. Тем не менее, часто желательно выполнять вычисления с изменяющейся (и произвольной) УНИВЕРСАЛЬНЫЕ АЛГОРИТМЫ 395 точностью. Для этого от алгоритма требуется независимость от точности вычислений и от особенностей машинного представления чисел. На самом деле многие алгоритмы также не зависят не только от машинного представления чисел, но и от конкретных математических (алгебраических) операций над данными. В этом случае операции са- ми по себе могут рассматриваться как переменные. Такие алгоритмы реализуются в виде обобщенных программ, написанных на основе абстрактных типов данных, кото- рые определяются пользователем в дополнение ко встроенным в язык типам данных. Соответствующие программные средства появились сначала в языке Симула-67, но со- временные объектно-ориентированные языки (такие как C++, см., например, [41, 47]) более удобны для обобщенного программирования. Алгоритмы компьютерной алгеб- ры, используемые в системах Mathematica, Maple, REDUCE и других, также высоко универсальны. Другой формой алгоритмов являются итерационные алгоритмы для решения диф- ференциальных уравнений (например, методы Эйлера, Эйлера-Коши, Рунге-Кутта, Адамса, несколько вариантов метода разностных приближений и подобные), методы вычисления элементарных и некоторых специальных функций, основанных на разло- жении в ряды Тейлора и на непрерывные дроби (приближения Паде). Эти алгоритмы не зависят от машинного представления чисел. Понятие обобщенных программ было введено многими авторами: например, в [30] такие программы называются ѕпрограммными схемамиї. В настоящей статье мы об- судим универсальные алгоритмы, реализованные в виде обобщенных программ, и их особые возможности. Эта статья тесно связана с работами [8, 31-33, 35, 36, 40], в кото- рых определяется понятие универсального алгоритма, и обсуждается программная и аппаратная реализации таких алгоритмов в связи с задачами идемпотентной матема- тики, см., например, [6, 39, 42, 53, 54]. Так называемый идемпотентный принцип соответствия (см. [32, 33]), связывающий идемпотентную математику с обычной математикой над полями, будет рассмотрен ни- же. В двух словах, существует соответствие между интересными, полезными и важны- ми конструкциями и результатами над полем вещественных (или комплексных) чисел и похожими конструкциями над идемпотентыми полукольцами. Это соответствие мо- жет быть сформулировано в духе хорошо известного принципа соответствия Н. Бора в квантовой механике: на самом деле, оба принципа тесно связаны (см. [31-33]). В извест- ном смысле, традиционную математику над числовыми полями можно рассматривать как ѕквантовуюї теорию, в то время как идемпотентную математику можно рассмат- ривать как ѕклассическуюї тень (или двойник) традиционной. Важно, что идемпо- тентный принцип соответствия верен для алгоритмов, компьютерных программ и их аппаратных реализаций. В квантовой механике принцип суперпозиции означает, что уравнение Шрјдингера (основное уравнение теории) линейно. Подобным образом, в идемпотентной математике принцип суперпозиции, сформулированный В. П. Масловым, означает, что некоторые важные и основные задачи и уравнения, которые являются нелинейными в обычном смысле (например, уравнение Гамильтона-Якоби, одно из основных уравнений клас- сической механики, которое также появляется во многих задачах оптимизации, или уравнение Беллмана и его разновидности и обобщения), могут рассматриваться как 396 Г.Л. Литвинов, А. Я. Родионов, С. Н. Сергеев, А. Н. Соболевский линейные над подходящими идемпотентными полукольцами, см. [4, 5, 31]. Заметим, что численные алгоритмы для бесконечномерных линейных задач над идемпотентными полукольцами (например, идемпотентное интегрирование, интеграль- ные операторы и преобразования, уравнения Гамильтона-Якоби и обобщенные уравне- ния Беллмана) имеют дело с соответствующими конечномерными приближениями. Сле- довательно, идемпотентная линейная алгебра является основой идемпотентного числен- ного анализа и, в особенности, теории дискретной оптимизации. Б. А. Карре [18, 19] (см. также [24-26]) использовал идемпотентную линейную алгеб- ру, чтобы показать, что различные задачи оптимизации для конечных графов могут быть сформулированы единым образом и сведены к решению дискретных матричных уравнений Беллмана, то есть некоторых систем линейных алгебраических уравнений над идемпотентными полукольцами. Он также обобщил основные алгоритмы вычисли- тельной линейной алгебры на идемпотентный случай и показал, что некоторые из них совпадают с алгоритмами, разработанными ранее для решения задач оптимизации. На- пример, метод Беллмана решения задачи нахождения кратчайшего пути соответствует разновидности метода Якоби решения системы линейных уравнений, тогда как алго- ритм Форда соответствует методу Гаусса-Зейделя. Уравнения Беллмана являются ос- новной темой данной статьи, а идеи Карре используются для получения новых универ- сальных алгоритмов. Подчеркнем, что универсальные алгоритмы, описанные в данной статье, могут быть интерпретированы как демонстрация принципа идемпотентной су- перпозиции и идемпотентного принципа соответствия между алгоритмами линейной алгебры и соответствующими задачами оптимизации на графах. Заметим, что многие алгоритмы решения матричных уравнений Беллмана можно найти в [8, 9, 13, 18, 19, 21, 24, 26, 35, 36, 48]. Другие задачи тропической линейной алгебры рассмотрены, например, в [16]. Также мы вкратце обсудим интервальный анализ над идемпотентными и поло- жительными полукольцами. Идемпотентный интервальный анализ возник в работах [3, 10, 39], где он применялся к матричному уравнению Беллмана. Позднее различные задачи, возникающие в интервальной идемпотентной линейной алгебре, были рассмот- рены, например, в работах [20, 22, 28, 44, 45]. Важно заметить, что интервалы над идем- потентными полукольцами образуют новое идемпотентное полукольцо. Следовательно, универсальные алгоритмы могут применяться к элементам этого нового полукольца и порождают интервальное расширение исходных алгоритмов. Данная статья посвящена программным реализациям универсальных алгоритмов решения матричных уравнений Беллмана над полукольцами. Раздел 1. содержит вве- дение в математику полуколец и, особенно, в тропическую (идемпотентную) матема- тику. В разделе 2. мы описываем несколько хорошо известных и новых универсальных алгоритмов линейной алгебры над полукольцами, связанных с дискретным матричным уравнением Беллмана и алгебраической задачей о кратчайшем пути. Эти алгоритмы тесно связаны с их прототипами в традиционной линейной алгебре, описанными, на- пример, в знаменитой книге Дж. Голуба и Ч. Ван Лоуна [1], которая служит основным источником таких прототипов. Следуя стилю [1], мы записываем их на языке MATLAB. Опыт программной реализации универсальных алгоритмов и дальнейшие перспективы их развития также рассматриваются. УНИВЕРСАЛЬНЫЕ АЛГОРИТМЫ 397 1. Математика полуколец 1.1. Основные определения Широкий класс универсальных алгоритмов связан с понятием полукольца. Напом- ним соответствующее определение (см., например, [23]). Рассмотрим множество S; снабженное двумя ассоциативными операциями: сложением и умножением ; та- кими, что сложение коммутативно, умножение дистрибутивно относительно сложения с обеих сторон, 0 (соответственно 1 ) нейтральный элемент по сложению (соответствен- но, по умножению), 0 x = x 0 = 0 для всех x 2 S; и 0 6= 1: Пусть полукольцо S частично упорядочено таким отношением ; что 0 наименьший элемент, и из нера- венства x y следует, что xz yz; xz yz и zx zy для всех x; y; z 2 S ; в этом случае полукольцо S называется положительным (см., например, [23]). Полукольцо S называется полуполем, если каждый его ненулевой элемент обратим. Полукольцо S называется идемпотентным, если xx = x для всех x 2 S: В этом случае сложение задает канонический частичный порядок на полукольце S по правилу: x y тогда и только тогда, когда x y = y: Любое идемпотентное коль- цо положительно относительно этого порядка. Заметим также, что x y = supfx; yg относительно канонического порядка. Впоследствии мы полагаем все идемпотентные кольца упорядоченными канонически. Мы говорим, что положительное (например, идемпотентное) полукольцо S полно, если оно полно как упорядоченное множество. Это означает, что для каждого подмно- жества T S существуют элементы sup T 2 S и inf T 2 S: Самые известные и важные примеры положительных полуколец это ѕчисловыеї полукольца, состоящие из (подмножества) вещественных чисел, упорядоченных обыч- ным линейным порядком 6 на множестве R: полукольцо R+ с обычными опера- циями = +; = и нейтральными элементами 0 = 0; 1 = 1; полукольцо Rmax = R [ f 1g с операциями = max; = + и нейтральными элементами 0 = 1; 1 = 0; полукольцо bR max = Rmax [ f1g; где x 1; x 1 = 1 для всех x; x 1 = 1 x = 1; если x 6= 0; и 0 1 = 1 0 и полукольцо S[a;b] max;min = [a; b]; где 1 6 a < b 6 +1 с операциями = max; = min и нейтральными элементами 0 = a; 1 = b: Полукольца Rmax; bR max; и S[a;b] max;min = [a; b] идемпотентны. Полуколь- ца bR max; S[a;b] max;min; bR+ = R+ S f1g полны как упорядоченные множества. Вспомним, что каждое частично упорядоченное множество может быть вложено в свое пополнение (минимальное полное множество, содержащее исходное). Полукольцо Rmin = R S f1g с операциями = min и = + и нейтральными элементами 0 = 1; 1 = 0 изоморф- но полукольцу Rmax: Полукольцо Rmax называют алгеброй макс-плюс. Полуполя Rmax и Rmin также называют тропическими алгебрами. Термин ѕтропическийї впервые появился в [51] для целочисленной версии алгебры макс-плюс, см. также [27, 42, 54]. Обозначим Matmn(S) множество всех матриц A=(aij) с m строками и с n столбцами, элементы которых принадлежат полукольцу S: Сумма AB матриц A;B 2 Matmn(S) и произведение AB матриц A 2 Matlm(S) и B 2 Matmn(S) определяются в соответ- 398 Г.Л. Литвинов, А. Я. Родионов, С. Н. Сергеев, А. Н. Соболевский ствии с обычными правилами линейной алгебры: A B = (aij bij) 2 Matmn(S) и AB = Mm k=1 aij bkj ! 2 Matln(S); где A 2 Matlm(S) и B 2 Matmn(S): Если полукольцо S положительно, то множество Matmn(S) упорядочено посред- ством отношения A = (aij) B = (bij) тогда и только тогда, когда aij bij в S для всех 2 6 i 6 m; 1 6 j 6 n: Умножение матриц согласуется с порядком в следующем смысле: если A;A0 2 Matlm(S); B;B0 2 Matmn(S) и A A0; B B0; то AB A0B0 в Matln(S): Более того, множество Matnn(S) квадратных матриц (n n) над (положительным идемпотентным) полукольцом S образует [положительное идемпотентное] полукольцо с нулевым элементом O = (oij); где oij = 0; 1 6 i; j 6 n; и единичным элементом I = (ij); где ij = 1; если i = j; и ij = 0 в противном случае. Множество Matnn типичный пример некоммутативного полукольца (при n > 1 ). 1.2. Операция замыкания Пусть положительное полукольцо S снабжено частичной унарной операцией замы- кания такой, что из a b следует, что a b и a = 1 (a a) = 1 (a a) во всей ее области определения. В частности, 0 = 1 по определению. Из этих аксиом следует, что a = 1 a a2 (a an); если n > 1: Таким образом, x можно рассматривать как ѕрегуляризованную суммуї ряда a = 1 a a2 : : : : В положительном полукольце, при условии, что оно замкнуто относительно взятия ограниченных упорядоченных верхних граней и операций и ; дистрибутивных относительно таких верхних граней, мы можем определить a := sup k>0 1 a : : : ak; (1.1) если последовательность правых частей ограничена. В этом случае a наименьшее решение уравнений x = ax1 и x = xa1; и ab наименьшее решение уравнений x = ax b и x = xa b: Если S полно, то замыкание очевидным образом определено для любого элемента x 2 S: В случае идемпотентного сложения (1.1) становится особенно приятным: a = M i>0 ai = sup i>0 ai: (1.2) В числовых полукольцах операция замыкания обычно реализуется очень просто: x = (1 x) 1; если x < 1 в R+; или в bR+; и x = 1; если x > 1 в bR+ ; x = 1; если x 1 в Rmax и bR max; x = 1; если x 1 в bR max; x = 1 для всех x в Smax;min [a;b]: Во всех остальных случаях x не определено. Операция замыкания в матричных полукольцах над положительным полукольцом S может быть определена УНИВЕРСАЛЬНЫЕ АЛГОРИТМЫ 399 по индукции: во-первых, имеем A = (a11) = (a 11) в Mat11(S): Во-вторых, для любого целого числа n > 1 и любой матрицы A = A11 A12 A21 A22 ; где A112Matkk(S); A122Matk n k(S); A212Matn k k(S); A222Matn k n k(S); 16k6 n; определим A = 0 @ A 11 A 11A12DA21A 11 A 11A12D DA21A 11 D 1 A; (1.3) где D = A22 A21A 11A12: Можно доказать, что из определения A следует, что равен- ства A = AA I = AA I выполняются, и, таким образом, A является ѕрегуля- ризованной суммойї ряда I A A2 : : : : Более того, в случае когда A определено как наименьшее решение уравнения A = AA I и A = AA I; можно показать, что оно удовлетворяет уравнению (1.3). Заметим, что это рекуррентное соотношение совпадает с формулами эскалаторно- го метода обращения матриц в обычной линейной алгебре над полями вещественных и комплексных чисел с точностью до используемых алгебраических операций. Следо- вательно, этот алгоритм матричного замыкания требует полиномиального от n числа операций, подробнее см. ниже. Пусть S положительное полукольцо. Матричное (дискретное, стационарное) урав- нение Беллмана имеет вид X = AX B; (1.4) где A 2 Matnn(S); X;B 2 Matns(S); а матрица X неизвестна. Пусть A замыкание матрицы A: Из тождества A = AA I следует, что матрица AB удовлетворяет этому уравнению. Также как и в скалярном случае, можно показать (при некотором дополнительном предположении) , что для положительных полуколец, если матрица A определена как в уравнении (1.1), то AB является наименьшим во множестве решений уравнения (1.4) относительно частично порядка на множестве Matns(S): Вспомним, что в идемпотентном случае A = M i>0 Ai = sup i>0 Ai: (1.5) Рассмотрим также случай, когда матрица A = (aij) размерности n n является строго верхнетреугольной (такой, что aij = 0 при i > j ) или строго нижнетреуголь- ной (такой, что aij = 0 при i 6 j ). В этом случае An = 0; нулевой матрице, и можно показать с помощью итераций X = AX I; что это уравнение имеет единственное решение, а именно A = I A : : : An 1: (1.6) Интересно, что в случае алгебры макс-плюс формула (1.6) работает для более широкого класса матриц. Ряд (1.5) сходится там тогда и только тогда, когда он может быть усечен до (1.6). Это тесно связано с интерпретацией матрицы A как ѕматрицы оптимальных путейї (см. ниже). 400 Г.Л. Литвинов, А. Я. Родионов, С. Н. Сергеев, А. Н. Соболевский 1.3. Взвешенные ориентированные графы и матрицы над полукольцами Пусть S полукольцо с нулем 0 и единицей 1: Хорошо известно, что любая квадратная матрица A = (aij) 2 Matnn(S) задает взвешенный ориентированный граф. Эта геометрическая конструкция включает три вида объектов: множество X из n элементов x1; : : : ; xn; называемых вершинами, множество всех упорядоченных пар (xi; xj) таких, что aij 6= 0; называется дугами, и отображение A: ! S такое, что A(xi; xj) = aij : Элементы aij из полукольца S называются весами дуг. Наоборот, любой взвешенный ориентированный граф с n вершинами задает некоторую (един- ственную) матрицу A 2 Matnn(S): Это определение допускает, что некоторые пары вершин могут быть не соединены, если соответствующий элемент матрицы A равен 0 и некоторые дуги могут быть петлями, если у матрицы A есть ненулевые элементы на главной диагонали. Вспомним, что последовательность вершин вида: p = (y0; y1; : : : ; yk); где k > 0 и (yi; yi+1) 2 ; i = 0; : : : ; k 1; называется путем длины k; соединяю- щим вершину y0 с вершиной yk: Обозначим множество всех таких путей Pk(y0; yk): Вес A(p) пути p 2 Pk(y0; yk) определяется как произведение весов дуг, соединяющих последовательные вершины пути: A(p) = A(y0; y1) A(yk 1; yk): По определению, вес ѕпутиї p 2 P0(xi; xj) длины k = 0 равен 1; если i = j; и 0 в противном случае. Для каждой матрицы A 2 Matnn(S) определим A0 = I = (ij) (где ij = 1; если i = j; и ij = 0 в противном случае) и Ak = AAk 1; k > 1: Пусть a[k] ij (i; j) -й элемент матрицы Ak: Легко проверяется, что a[k] ij = M i0=i; ik=j 16i1;:::;ik 16n ai0i1 aik 1ik : Таким образом, a[k] ij точная верхняя грань множества весов, отвечающих всем путям длины k; соединяющим вершину xi0 = xi с вершиной xik = xj : УНИВЕРСАЛЬНЫЕ АЛГОРИТМЫ 401 Пусть матрица A определена в соответствии с формулой (1.5). Обозначим элемен- ты матрицы A как a ij ; i; j = 1; : : : ; n ; тогда a ij = M 06k<1 M p2Pk(xi;xj ) A(p): Матрица A является решением известного алгебраического обобщения задачи на нахождение оптимального пути: для каждой пары (xi; xj) требуется вычислить точ- ную верхнюю грань весов всех путей (произвольной длины), соединяющих вершину xi с вершиной xj : Операция замыкания в полукольце матриц хорошо изучена (см., на- пример, [2, 6, 13, 16, 18, 19, 21, 23, 25, 26, 39] и ссылки в них). П р и м е р 1.1 (Задача о кратчайшем пути). Пусть S = Rmin; т. е. веса являются вещественными числами. В этом случае A(p) = A(y0; y1) + A(y1; y2) + + A(yk 1; yk): Если элемент aij задает длину дуги (xi; xj) в некоторой метрике, тогда a ij длина кратчайшего пути, соединяющего xi с xj : П р и м е р 1.2 (Задача о пути максимальной ширины). Пусть S = R [ f0; 1g с = max; = min : Тогда a ij = max p2 S k>1 Pk(xi;xj ) A(p); A(p) = min(A(y0; y1); : : : ;A(yk 1; yk)): Если элемент aij задает ѕширинуї дуги (xi; xj); тогда ширина пути p определяется как наименьшая ширина дуги, принадлежащей этому пути, а элемент a ij является максимальной шириной всех путей, соединяющих xi с xj : П р и м е р 1.3 (Простая задача динамического программирования). Пусть S = Rmax; и положим, что aij задает выплаты при переходе от xi к xj : Определим вектор B = (bi) 2 Matn1(Rmax); элементы которого bi задают терминальные выпла- ты, соответствующие выходу из графа в узле xi: Конечно, отрицательные выплаты (то есть убытки) также допускаются. Пусть m итоговая выплата, соответствующая пути p 2 Pk(xi; xj); т. е. m = A(p) + bj : Тогда легко проверить, что наибольшая выплата, которая может быть достигнута на путях длины k; начинающихся в узле xi; равна (AkB)i; и наибольшая выплата, до- стижимая без ограничений на длину пути, равна (AB)i: Заметим, что если длина пути неограничена, то наибольшая выплата может быть равна бесконечности, как и соответствующее значение (AB)i: П р и м е р 1.4 (Задача обращения матрицы). Заметим, что в формулах этого раз- дела мы используем дистрибутивность умножения относительно сложения ; но не 402 Г.Л. Литвинов, А. Я. Родионов, С. Н. Сергеев, А. Н. Соболевский используем аксиому идемпотентности. Таким образом, алгебраическое обобщение зада- чи об оптимальном пути может быть сформулировано также и для неидемпотентного полукольца S (см., например, [48]). Например, если S = R; то A = I + A + A2 + = (I A) 1: Если kAk > 1; но матрица обратима, то это выражение определяет регуляризованную сумму расходящегося матричного степенного ряда P i>0 Ai: Подчеркнем, что эта связь между матричной операцией замыкания и решениями уравнения Беллмана порождает несколько различных алгоритмов для вычисления мат- ричного замыкания. Все эти алгоритмы являются адаптациями хорошо известных алго- ритмов традиционной вычислительной линейной алгебры, таких как метод исключения Гаусса-Жордана, различные итерационные и эскалаторные схемы и т. д. Это особый случай идемпотентного принципа суперпозиции (см. ниже). На самом деле, теория дискретного стационарного уравнения Беллмана может быть развита с использованием равенства A = AA I как дополнительной аксиомы без какой-либо существенной интерпретации (так называемые замкнутые полукольца, см., например, [23, 30, 48]). Такая теория может быть основана на следующих равенствах, которые в случае идемпотентных полуколец можно проверить, используя связь с за- дачей на нахождение оптимального пути. Эти равенства также хорошо известны и в случае вещественных чисел с обычной арифметикой: (A B) = (AB)A; (AB)A = A(BA): (1.1) 1.4. Интервальный анализ над положительными полукольцами Традиционный интервальный анализ это нетривиальная и популярная математи- ческая область, см., например, [12, 22, 29, 43, 46]. ѕИдемпотентнаяї версия интервально- го анализа (и кроме того, интервального анализа над положительными полукольцами) появилась в [3, 10, 39]. Довольно много работ по этой теме появились позднее, см., на- пример, [20, 22, 28, 44, 45]. Интервальный анализ над положительным полукольцом R+ обсуждается в [15]. Пусть множество S частично упорядочено отношением : Замкнутый интервал в S - это подмножество вида x = [x; x] = f x 2 S j x x x g; где элементы x x называются нижней и верхней границей интервала x: Порядок порождает частичное упорядочение на множестве всех замкнутых интервалов в S : x y тогда и только тогда, когда x y и x y: Слабым интервальным расширением I(S) положительного полукольца S назы- вается множество всех замкнутых интервалов в S; снабженное операциями и ; определенными следующим образом: x y = [x y; x y]; x y = [x y; x y]; и частичным порядком, порожденным порядком на S: Операция замыкания в I(S) опре- деляется как x = [x; x]: Существуют другие интервальные расширения (в том числе, так называемые сильные интервальные расширения [39]), но слабое расширение более удобно. УНИВЕРСАЛЬНЫЕ АЛГОРИТМЫ 403 Расширение I(S) положительное; I(S) идемпотентно, если S идемпотентное полукольцо. Универсальный алгоритм над S может применяться к I(S); и мы полу- чим интервальную версию исходного алгоритма. Обычно обе версии имеют одинако- вую сложность. Для дискретного стационарного уравнения Беллмана и соответству- ющей оптимизационной задачи на графах интервальный анализ подробно рассматри- вается в [3, 39]. Другие задачи идемпотентной линейной алгебры были рассмотрены в [20, 22, 28, 44, 45]. Идемпотентная математика оказывается проще ее традиционного аналога. Напри- мер, в традиционной интервальной арифметике перемножение интервалов не дистри- бутивно относительно сложения, тогда как в идемпотентной интервальной арифметике эта дистрибутивность сохраняется. Более того, в традиционном интервальном анали- зе множество всех квадратных интервальных матриц заданного порядка не образует даже полугруппы по матричному умножению: эта операция не ассоциативна, посколь- ку дистрибутивность теряется в традиционной интервальной арифметике. Наоборот, в идемпотентном (и положительном) случае ассоциативность сохраняется. Наконец, в обычном интервальном анализе некоторые задачи линейной алгебры, такие как реше- ние линейных систем интервальных уравнений, может быть очень сложным (точнее говоря, они NP -полные, см. [29] и ссылки там же). В [3, 39] показано, что в идем- потентном случае решение интервальных линейных систем требует полиномиального количества операций (так же как и для обычного метода исключения Гаусса). Два свой- ства, которые делают идемпотентную интервальную арифметику такой простой, это монотонность арифметических операций и положительность всех элементов идемпо- тентного полукольца. Отметим, что интервальные оценки в идемпотентной математике обычно являются точными. В традиционной теории такие оценки, как правило, слишком пессимистиче- ские. 1.5. Идемпотентный принцип соответствия Существует нетривиальная аналогия между математикой полуколец и квантовой механикой. Например, поле вещественных чисел можно рассматривать как ѕкванто- вый объектї относительно идемпотентных полуколец, которые рассматриваются как ѕклассическиеї и ѕполуклассическиеї объекты относительно поля вещественных чи- сел. Пусть R поле вещественных чисел, а R+ подмножество всех неотрицательных чисел. Рассмотрим следующую замену переменных: u 7! w = h ln u; где u 2 R+ n f0g; h > 0 ; таким образом, u = ew=h; w 2 R: Обозначим через 0 дополнительный элемент 1 и через S расширенную вещественную прямую R[f0g: Эта замена переменных имеет естественное расширение Dh на все S с Dh(0) = 0 ; также мы положим Dh(1) = 0 = 1: Обозначим Sh множество S; снабженное двумя операциями h (обобщенное сло- жение) и h (обобщенное умножение), такое что Dh гомоморфизм из fR+;+; g в 404 Г.Л. Литвинов, А. Я. Родионов, С. Н. Сергеев, А. Н. Соболевский fS;h;hg: Это означает, что Dh(u1 + u2) = Dh(u1) h Dh(u2); Dh(u1 u2) = Dh(u1) h Dh(u2); то есть w1 h w2 = w1 + w2 и w1 h w2 = h ln(ew1=h + ew2=h): Легко показать, что w1 h w2 ! maxfw1;w2g при h ! 0: Алгебраически R+ и Sh изоморфны, поэтому мы получаем Rmax как результат предельной деформации R+: Мы подчеркиваем очевидную аналогию с процедурой деквантования, где h аналог постоянной Планка. Таким образом, R+ (или R) иг- рает роль ѕквантового объектаї, а Rmax ѕклассическогої или ѕполуклассическогої объекта, который получается как результат деквантования квантового объекта R+: В случае Rmin соответствующая процедура деквантования порождается заменой пере- менных u 7! w = h ln u: В случае поля вещественных или комплексных чисел, идемпотентное деквантование в полукольцо Rmax (или Rmin ) определяется с помощью композиции отображения x 7! jxj и деформации, описанной выше. В то же время идемпотентное деквантование играет важную роль в методологии и в философии идемпотентной математики как обоснование возможности переноса с основного поля (как правило, вещественных или комплексных чисел) на идемпотентное полукольцо многих математических понятий, конструкций и результатов, см. [31, 33]. Идемпотентное деквантование приводит к следующей формулировке идемпотентного принципа соответствия [31, 33]: Существует эвристическое соответствие между интересными, полезны- ми и важными результатами над полем вещественных (или комплексных) чисел и похожими конструкциями и результатами над идемпотентными полукольцами в духе принципа соответствия Н. Бора в квантовой меха- нике. Квантовая механика Тртрвараварывпкуецкецуе Поля вещественных и комплексных чисел Принцип соответствия Н. Бора Традиционная математика Идемпотентная математика Классическая механика Идемпотентные полукольца и полуполя Идемпотентный принцип соответствия Идемпотентная математика может рассматриваться как ѕклассическая тень (или двойник)ї традиционной математики над полями. Систематическое приложение этого принципа соответствия ведет ко множеству теоретических и прикладных результатов, см., например, [31, 33, 36, 39, 42, 53, 54]. Связь с квантовой физикой подробно обсуждается, например, в [31]. УНИВЕРСАЛЬНЫЕ АЛГОРИТМЫ 405 Теперь мы обсудим идею о том, как идемпотентный принцип соответствия мог бы привести к унифицированному подходу к разработке аппаратного обеспечения. Подроб- нее об этой идее и ее осуществлении см. в [35, 40]. Наиболее важные и давно известные численные алгоритмы имеют множество аппа- ратных реализаций в виде технических устройств или специальных процессоров. Часто эти устройства могут использоваться как прототипы для новых аппаратных устройств, в результате простой замены обычных арифметических операций на их аналоги в полу- кольцах (и дополнительные средства для порождения нейтральных элементов 0 и 1 ). Разумеется, случай числовых полуколец, состоящих из вещественных чисел (возможно, кроме нейтральных элементов) и полуколец числовых интервалов, является самым про- стым и естественным. Заметим, что для полуколец (включая Rmax и Rmin ) операция деления также определена. Эффективные технические идеи и решения могут быть заимствованы из прототипов и реализованы в новых аппаратных устройвах. Таким образом, принцип соответствия порождает эвристический метод для разработки аппаратных средств. Заметим, что для получения патента необходимо представить так называемую ѕформулу изобретенияї, где требуется указать прототип предлагаемого устройства и разницу между устрой- ством и его прототипом. Рассмотрим (как типичный пример) самый важный и давно известный алгоритм вычисления скалярного произведения двух векторов: (x; y) = x1y1 + x2y2 + + xnyn: (1.1) Универсальная версия формулы (1.1) для любого полукольца A очевидна: (x; y) = (x1 y1) (x2 y2) (xn yn): (1.2) В случае A = Rmax эта формула принимает следующий вид: (x; y) = maxfx1 + y1; x2 + y2; ; xn + yng: (1.3) Это вычисление стандартно для многих алгоритмов оптимизации, поэтому полез- но создать аппаратное средство для вычисления (1.3). Существует множество различ- ных устройств и патентов для вычисления (1.1), и каждое такое устройство можно использовать в качестве прототипа для создания нового устройства для вычисления (1.3) и даже (1.2). Многие процессоры для перемножения матриц и других алгорит- мов линейной алгебры основаны на вычислении скалярного произведения и на соответ- ствующих ѕэлементарныхї устройствах. Используя современные технологии, можно создать дешевые специализированные многопроцессорные чипы и систолические мас- сивы элементарных процессоров, реализующих универсальные алгоритмы. См., напри- мер, [35, 40, 48], где обсуждаются систолические массивы и вопросы параллельного вы- числения в рамках задачи алгебраического пути. В частности, существует систоли- ческий массив из n(n + 1) элементарных процессоров, которые выполняют алгоритм исключения Гаусса-Жордана и могут решать задачу алгебраического пути за 5n 2 временных шагов. 406 Г.Л. Литвинов, А. Я. Родионов, С. Н. Сергеев, А. Н. Соболевский 2. Некоторые универсальные алгоритмы линейной алгебры В этом разделе мы рассмотрим универсальные алгоритмы вычисления A и AB: Мы начнем с простых методов: эскалаторного метода и метода исключения Гаусса- Жордана в подразделе 2.1., а потом рассмотрим модификацию этих методов для слу- чая систем Тјплица (подраздел 2.2.). Универсальное LDM-разложение для уравнения Беллмана описано в подразделе 2.3., а потом приводится адаптация этого метода для случаев симметричной и ленточной матрицы (подраздел 2.4.) . Далее обсуждаются две итерационные схемы, а также реализация универсальных алгоритмов (подраздел 2.7.). Сами алгоритмы будут записаны на языке MATLAB, следуя книге Голуба и ван Лоуна. Это сделано по следующим причинам: 1) желание упростить сравнение алгорит- мов с их аналогами, взятых преимущественно из [1]; 2) поскольку язык MATLAB раз- работан для матричных вычислений и удобен для их записи. Мы не будем формально описывать правила нашего MATLAB-образного языка, предпочитая только выделить следующие важные особенности: 1. Основные арифметические операции: a b; a b и a: 2. Эти же операции для векторов следуют правилам MATLAB. 3. Мы используем основные ключевые слова MATLAB: ѕforї, ѕifї и ѕendї, похожие на ключевые слова других языков программирования (C++, Java, ...). Дадим несколько примеров универсальных матричных вычислений на нашем языке. П р и м е р 2.1. v(1 : j) = a(1 : j; k) означает, что результат умножения (ска- лярного) первых j компонент k -й колонки матрицы A на замыкание присваивается первым j компонентам вектора v: П р и м е р 2.2. a(i; j) = a(i; j)a(i; 1 : n)a(1 : n; j) означает, что мы добавляем () к элементу aij матрицы A результат скалярного умножения (универсального) i -й строки и j -го столбца матрицы A (предполагается, что A - матрица n n ). П р и м е р 2.3. a(1 : n; k)b(l; 1 : m) означает внешнее произведение k -го столбца матрицы A и l -й строки матрицы B: Элементы матрицы-результата C =(cij) равны cij =aik blj для всех i = 1; : : : ; n и j = 1; : : : ; m: П р и м е р 2.4. x(1 : n) y(n : 1 : 1) - скалярное произведение вектора x и вектора y; чьи компоненты берутся в обратном порядке: правильное алгебраическое выражение Ln i=1 xi yn+1 i: П р и м е р 2.5. Следующий цикл вычисляет тот же результат, что и предыдущий пример: s = 0 for i = 1 : n s = s x(i) x(n + 1 i) end УНИВЕРСАЛЬНЫЕ АЛГОРИТМЫ 407 2.1. Эскалаторная схема и исключение Гаусса-Жордана Сначала мы проанализируем базовый эскалаторный метод, основанный на определе- нии замыкания матрицы (1.3). Пусть A квадратная матрица. Замыкания ее главных подматриц Ak могут быть найдены по индукции, начиная с A1 = (a11); замыкания первого диагонального элемента. В общем случае мы выражаем Ak+1 как Ak+1 = Ak gk hTk ak+1 ; предполагая, что мы уже нашли замыкание для Ak: В этой формуле gk и hk столбцы с k элементами и ak+1 скаляр. Мы также выражаем A k+1 как A k+1 = Uk vk wT k uk+1 : Используя (1.3), мы получаем, что uk+1 = (hTk A kgk ak+1); vk = A kgkuk+1; wT k = uk+1hTk A k; Uk = A kgkuk+1hTk A k A k: (2.1) Алгоритм, основанный на (2.1), записывается следующим образом. Алгоритм 1. Эскалаторный метод вычисления для A Вход: матрица A размером n n с элементами a(i; j); также используется для хранения результата вычисления и промежуточных результатов вычисления. a(1; 1) = (a(1; 1)) for i = 1 : n 1 Ag = a(1 : i; 1 : i) a(1 : i; i + 1) hA = a(i + 1; 1 : i) a(1 : i; 1 : i) a(i + 1; i + 1) = a(i + 1; i + 1) a(i + 1; 1 : i) Ag(1 : i; 1) a(i + 1; i + 1) = (a(i + 1; i + 1)) a(1 : i; i + 1) = a(i + 1; i + 1) Ag a(i + 1; 1 : i) = a(i + 1; i + 1) hA a(1 : i; 1 : i) = a(1 : i; 1 : i) Ag a(i + 1; i + 1) hA end В полной аналогии со своим прототипом из линейной алгебры алгоритм требует n3 + O(n2) операций сложения ; n3 + O(n2) операций умножения и n операций взятия алгебраического замыкания. Линейно-алгебраический прототип метода, запи- санного выше, также называется в литературе [11, 18] методом окаймления. В качестве альтернативы, мы можем получить решение X = AX B в результате процесса исключения, формальное объяснение которого дано ниже. Если A опреде- ляется как L i>0 Ai (включая скалярный случай), тогда AB наименьшее решение 408 Г.Л. Литвинов, А. Я. Родионов, С. Н. Сергеев, А. Н. Соболевский X = AX B для всех A и B подходящих размеров. В этом случае, решение, найден- ное с помощью процесса исключения, совпадает с AB: Для матрицы A = (aij) и вектор-столбцов x = (xi) и b = (bi) (ограничимся без потери общности случаем, когда x и b вектор-столбцы) уравнение Беллмана x = Ax b может быть записано как 0 BBB@ x1 x2 ... xn 1 CCCA = 0 BBB@ a11 a12 : : : a1n a21 a22 : : : a2n ... ... . . . ... an1 an2 : : : ann 1 CCCA 0 BBB@ x1 x2 ... xn 1 CCCA 0 BBB@ 1 0 : : : 0 0 1 : : : 0 ... ... . . . ... 0 0 : : : 1 1 CCCA 0 BBB@ b1 b2 ... bn 1 CCCA : (2.2) После выражения x1 через x2; : : : ; xn из первого уравнения и подстановки этого выражения для x1 во все другие уравнения от второго до n -го мы получим: 0 BBB@ x1 x2 ... xn 1 CCCA = 0 BBB@ 0 (a11)a12 : : : (a11)a1n 0 a22 (a21(a11)a12) : : : a2n (a21(a11)a1n) ... ... . . . ... 0 an2 (an1(a11)a12) : : : ann (an1(a11)a1n) 1 CCCA 0 BBB@ x1 x2 ... xn 1 CCCA 0 BBB@ (a11) 0 : : : 0 a21(a11) 1 : : : 0 ... ... . . . ... an1(a11) 0 : : : 1 1 CCCA 0 BBB@ b1 b2 ... bn 1 CCCA (2.3) Заметим, что нетривиальные элементы в обеих матрицах занимают дополнительные по отношению друг к другу позиции, поэтому во время вычислений обе матрицы могут храниться в одной квадратной матрице C(k): Обозначим ее элементы через c(k) ij ; где k число исключенных переменных. После l 1 исключения мы имеем xl = (c(l 1) ll )bl; c(l) il = c(l 1) il (c(l 1) ll ); i = 1; : : : ; l 1; l + 1; : : : ; n c(l) ij = c(l 1) ij c(l 1) il (c(l 1) ll )c(l 1) lj ; i; j = 1; : : : ; l 1; l + 1; : : : ; n c(l) li = (c(l 1) ll )c(l 1) li ; i = 1; : : : ; l 1; l + 1; : : : ; n (2.4) После n исключений имеем x = C(n)b: Выбирая в качестве b любой вектор с одной координатой, равной 1; а остальными, равными 0; получим C(n) = A: Запишем следующий алгоритм, основанный на рекурсии (2.4). Алгоритм 2. Метод исключения Гаусса-Жордана для вычисления A: Вход: матрица A размеров n n с элементами a(i; j); также используется для хранения окончательного результата и промежуточные результаты вычислений. for i = 1 : n a(i; i) = (a(i; i)) УНИВЕРСАЛЬНЫЕ АЛГОРИТМЫ 409 for k = 1 : n if k 6= i a(k; i) = a(k; i) a(i; i) end end for k = 1 : n for j = 1 : n if k 6= i & j 6= i a(k; j) = a(k; j) a(k; i) a(i; j) end end for j = 1 : n if j 6= i a(i; j) = a(i; i) a(i; j) end end end Замечание. Алгоритм 2 может рассматриваться как ѕуниверсальный алгоритм Флойда-Уоршаллаї, обобщающий хорошо известные алгоритмы Уоршалла и Флойда для вычисления транзитивного замыкания графа и всех оптимальных путей на графе, см., например, [7]. В свою очередь, эти методы могут рассматриваться как частные случаи алгоритма 2 для тропических и булевых полуколец. Алгоритм 2 также близок к методу ѕпополненияї Ершова для обращения матриц и решения систем вида Ax = b в классической линейной алгебре, подробности см. [11, глава 2]. 2.2. Системы Тјплица Запишем эскалаторный метод для нахождения решения x = Ab уравнения x = Ax b; где x и b вектор-столбцы, аналогичный описанному в предыдущем разделе эскалаторному методу для нахождения A: Начнем с x(1) = A1 b1: Пусть x(k) - вектор, найденный после (k 1) шага. Представим его в виде x(k+1) = z xk+1 : Используя (2.1), мы получаем xk+1 = uk+1(hTk x(k) bk+1); z = x(k) A kgkxk+1: (2.1) Требуется вычислить A kgk: Далее мы покажем, что это вычисление может быть про- ведено гораздо эффективнее, когда A - симметричная тјплицева матрица. Матрица A 2 Matnn(S) называется тјплицевой (или матрицей Тјплица), если существуют скаляры r n+1; : : : ; r0; : : : ; rn 1; такие, что Aij = rj i для всех i и j: Другими словами, тјплицевы матрицы обладают тем свойством, что их элементы по- стоянны вдоль любой линии, параллельной главной диагонали (и вдоль самой главной диагонали). Например, матрица A = 0 BB@ r0 r1 r2 r3 r 1 r0 r1 r2 r 2 r 1 r0 r1 r 3 r 2 r 1 r0 1 CCA 410 Г.Л. Литвинов, А. Я. Родионов, С. Н. Сергеев, А. Н. Соболевский тјплицева. Такие матрицы не обязательно симметричные, однако они всегда персим- метричные, то есть симметричные относительно побочной диагонали. Это свойство алгебраически выражается как A = EnATEn; где En = [en; : : : ; e1]: Через ei мы обо- значаем столбец, чей i -й элемент равен 1; а другие элементы равны 0: Из свойства E2n = In (где In единичная матрица n n ) следует, что произведение двух персим- метричных матриц является персимметричным. Следовательно, любая степень пер- симметричной матрицы персимметрична, как и замыкание персимметричной матрицы. Таким образом, если A персимметрична, то EnA = (A)TEn: (2.2) Далее мы рассматриваем только симметричные тјплицевы матрицы. Рассмотрим уравнение y = Tny r(n); где r(n) = (r1; : : : ; rn)T и Tn определяется скалярами r0; r1; : : : ; rn 1 так, что Tij = rjj ij для всех i и j: Это обобщение задачи Юла- Уокера [1]. Предположим, что мы получили наименьшее решение y(k) системы y = Tky r(k) для некоторого k такого, что 1 6 k 6 n 1; где Tk главная k k подматрица матрицы Tn: Тогда Tk+1 можно записать в следующем виде: Tk+1 = Tk Ekr(k) r(k)TEk r0 ; а y(k+1) и r(k+1) - как y(k+1) = z k ; r(k+1) = r(k) rk+1 : Используя (2.1), (2.2) и тождество T k r(k) = y(k); мы получаем, что k = (r0 r(k)T y(k))(r(k)TEky(k) rk+1); z = Eky(k)k y(k): Обозначим k = r0r(k)T y(k): Следующая выкладка показывает, что k может быть найден рекурсивно, если ( k 1) 1 существует: k = r0 [r(k 1)T rk] Ek 1y(k 1)k 1 y(k 1) k 1 = r0 r(k 1)T y(k 1) (r(k 1)TEk 1y(k 1) rk)k 1 = k 1 ( k 1) 1 2k 1: (2.3) Существование ( k 1) 1 не гарантировано, и поэтому мы записываем две версии нашего алгоритма, первая включает (2.3), а вторая не включает. Мы запишем обе эти версии в одной программе и выделим выражения, которые относятся только к первой или только к второй, с помощью комментариев языка MATLAB %1 и %2 соответ- ственно. Собирая выражения для k; k и z; мы получим следующее рекурсивное УНИВЕРСАЛЬНЫЕ АЛГОРИТМЫ 411 выражение: k = r0 r(k)T y(k); %2 k = k 1 ( k 1) 1 2k 1; %1 k = (k) (r(k)TEky(k) rk+1); y(k+1) = Eky(k)k y(k) k : (2.4) Рекурсия (2.4) это обобщенная версия метода Дурбина для задачи Юла-Уокера, про- тотип алгоритма см. в [1, алгоритм 4.7.1]. Алгоритм 3. Метод Дурбина для задачи Юла-Уокера для уравнений Беллмана с сим- метричной матрицей Тјплица. Вход: r0 : скаляр, r : n 1 1 вектор; y(1) = r 0 r(1) = r0 %1 = r 0 r(1) for k = 1 : n 1 = r0 r(1 : k) y(1 : k) %2 = () 1 2 %1 = (r(k : 1 : 1) y(1 : k) r(k + 1)) z(1 : k) = y(1 : k) y(k : 1 : 1) y(1 : k) = z(1 : k) y(k + 1) = end Выход: вектор y: В общем случае алгоритм требует 3=2n2 + O(n) операций и каждой и только n2 + O(n) операций и ; если допускаются инверсии алгебраических замыканий (как обычно, требуется только n таких замыканий в обоих случаях). Рассмотрим задачу нахождения x(n) = T nb(n); где Tn заданы как выше, а век- тор b(n) = (b1; : : : ; bn) произволен. Также мы введем вектор-столбец y(k); являющийся решением задачи Юла-Уокера: y(k) = T k r(k): Основная идея состоит в нахождении выражения для x(k+1) = T k+1b(k+1); включающих x(k) и y(k): Запишем x(k+1) и b(k+1) следующим образом: x(k+1) = v k ; b(k+1) = b(k) bk+1 : Используя персимметрию T k и тождества T k bk = x(k) и T k rk = y(k); преобразуем выражения (2.1) и получим, что k = (r0 r(k)T y(k)) (r(k)TEkx(k) bk+1); v = Eky(k)k x(k): 412 Г.Л. Литвинов, А. Я. Родионов, С. Н. Сергеев, А. Н. Соболевский Коэффициент r0 r(k)T y(k) = k может быть выражен рекурсивно: k = k 1 ( k 1) 1 (k 1)2; если замыкание (k 1) обратимо. Учитывая эту возможность, получаем следующие выражения: k = r0 r(k)T y(k); %2 k = k 1 ( k 1) 1 2k 1; %1 k = k (r(k)TEkx(k) bk+1); x(k+1) = Eky(k)k x(k) k : (2.5) Выражения (2.4) и (2.5) дают следующее обобщение алгоритма Левинсона решения линейных симметричных систем Тјплица, см. прототип в [1, алгоритм 4.7.2]. Алгоритм 4. Метод Левинсона для систем Беллмана с симметричной матрицей Тјплица Вход: r0 : скаляр, r : 1 n 1 вектор-строка; b : n 1 вектор-столбец. y(1) = r 0 r(1) ; x(1) = r 0 b(1) ; = r0 %1 = r 0 r(1) for k = 1 : n 1 = r0 r(1 : k) y(1 : k) %2 = () 1 2 %1 = (r(k : 1 : 1) x(1 : k) b(k + 1)) v(1 : k) = x(1 : k) y(k : 1 : 1) x(1 : k) = v(1 : k) x(k + 1) = if k < n 1 = (r(k : 1 : 1) y(1 : k) r(k + 1)) z(1 : k) = y(1 : k) y(k : 1 : 1) y(1 : k) = z(1 : k) y(k + 1) = end end Выход: вектор x: В общем случае алгоритм требует (5=2)n2 + O(n) операций и каждая, и только 2n2 + O(n) операций и ; если алгебраические замыкания обратимы (как обычно, требуется только n замыканий в обоих случаях). УНИВЕРСАЛЬНЫЕ АЛГОРИТМЫ 413 2.3. Разложение LDM Над классическими числовыми полями матричные уравнения вида AX = B можно решать с помощью разложения матрицы A в произведение A = LDM; где L и M нижняя и верхняя треугольные матрицы с единичной диагональю соответственно, а D диагональная матрица. Мы построим похожее разложение для решения уравнения Беллмана X = AX B над полукольцами. В случае системы AX = B над классическими числовыми полями, разложение A = LDM приводит к следующему разложению AX = B : LZ = B; DY = Z; MX = Y: (2.1) Следовательно, A 1 = M 1D 1L 1; (2.2) если матрица A обратима. По существу, достаточно найти матрицы L; D и M; так как после этого линейная система AX = B решается путем комбинации прямой под- становки Z (решение системы уравнений с нижне-треугольной матрицей), тривиаль- ной инверсии диагональной матрицы Y и обратной подстановки X (решение системы уравнений с верхне-треугольной матрицей). Используя разложение LDM для AX = B как прототип, мы можем записать аналогичное разложение для X = AX B : Z = LZ B; Y = DY Z; X = MX Y: (2.3) Тогда A = MDL: (2.4) Тройка (L;D;M); состоящая из нижней треугольной, диагональной и верхней тре- угольной матриц, называется LDM -разложением матрицы A; если выполнены соот- ношения (2.3) и (2.4). Заметим, что в этом случае главные диагонали матриц L и M нулевые. Другими словами, матрица L является строго нижне-треугольной, а M является строго верхне-треугольной. Наша универсальная модификация LDM -разложения схожа с универсальной мо- дификацией LU -разложения, предложенной Карре в [18, 19], см. также ниже. Универ- сальные алгоритмы LDM и LU -разложения можно применять и в случае обычных числовых полей, однако лишь в том случае, если не требуется операция ѕвыбора ве- дущего элементаї (pivoting). В этом случае они вычисляют классические LDM и LU разложения для обращения матрицы I A: Если A симметричная матрица над полукольцом с коммутативным умножением, то количество вычислений может быть уменьшено вдвое, так как M и L переходят друг в друга при транспозиции. Соответствующие алгоритмы будут подробно рассмот- рены ниже. Мы начнем со случая треугольной матрицы A = L (или A = M ). Тогда нахождение X сводится к прямой (или обратной) замене. В этих случаях уравнение X = AX B имеет единственное решение, которое может быть найдено с помощью простых алго- ритмов, приведенных ниже. В этих алгоритмах B вектор (обозначенный через b ), 414 Г.Л. Литвинов, А. Я. Родионов, С. Н. Сергеев, А. Н. Соболевский однако их (после соответствующей модификации) можно применять и в случае, когда B матрица любого подходящего размера. Нас интересует случай строго нижнетре- угольной матрицы, соответственно строго верхнетреугольной, когда aij = 0 для i 6 j; соответственно aij = 0 для i > j: Алгоритм 5. Прямая подстановка. Вход: Строго нижнетреугольная n n -матрица l ; n 1 вектор b: y = b for k = 2 : n y(k) = l(k; 1 : k 1) y(1 : k 1) end Выход: вектор y: Алгоритм 6. Обратная подстановка. Вход: Строго верхнетреугольная n n -матрица m; n 1 вектор b: y = b for k = n 1 : 1 : 1 y(k) = m(k; k + 1 : n) y(k + 1 : n) end Выход: вектор y: Оба алгоритма требуют n2=2 + O(n) операций и и никаких алгебраических за- мыканий. После LDM -разложения нам нужно вычислить замыкание диагональной матрицы. Это делается поэлементно. Сформулируем алгоритм LDM -разложения, т. е. вычисление матриц L; D и M; удовлетворяющих (2.3) и (2.4). Доказательство будет приведено ниже. Алгоритм 7. LDM-разложение (версия 1). Вход: n n -матрица A с элементами a(i; j); также используется для хранения результата вычислений и промежуточных результатов вычислений. for j = 1 : n 1 v(j) = (a(j; j)) a(j + 1 : n; j) = a(j + 1 : n; j) v(j) a(j + 1 : n; j + 1 : n) = a(j + 1 : n; j + 1 : n) a(j + 1 : n; j) a(j; j + 1 : n) a(j; j + 1 : n) = v(j) a(j; j + 1 : n) end Алгоритм требует n3=3 + O(n2) операций и и n 1 операций алгебраического замыкания. УНИВЕРСАЛЬНЫЕ АЛГОРИТМЫ 415 Строго треугольная матрица L записывается в нижнем треугольнике, строго верх- нетреугольная матрица M в верхнем треугольнике, и диагональная матрица D на диагонали матрицы, вычисленной алгоритмом 7. Теперь покажем, что A = MDL: Заметим, что приведенное ниже доказа- тельство является модификацией доказательства [14], данного для универсального LU -разложения. Идея заключается в том, чтобы определить универсальный аналог ѕматриц Гауссаї нижне-треугольных или верхне-треугольных, состоящих из одного столбца или одной строки, соответственно). A = a11 h(1) g(1) B(1) : (2.5) Можно проверить, что A = 1 h(1)a 11 On 11 In 1 a 11 O1n 1 On 11 (h(1)a 11g(1) B(1)) 1 O1n 1 a 11g(1) In 1; ; (2.6) где умножение в правой части ведет к выражению, полностью аналогичному (2.1), где (h(1)a 11g(1) B(1)) играет роль uk+1: Здесь и далее Okl означает k l -матрицы, состоящие только из нулей, а Il означает единичную матрицу размера l: Это также можно переписать как A = M 1D 1(A(2))L 1; (2.7) M1 = O h(1)a 11 O(n 1)1 O(n 1)(n 1) ; D1 = a11 O1(n 1) O(n 1)1 O(n 1)(n 1) ; A(2) = O11 O1(n 1) O(n 1)1 R(2) ; L1 = O11 O1(n 1) a 11g(1) O(n 1)(n 1) ; R(2) = h(1)a 11g(1) B(1): (2.8) Здесь мы используем, в частности, что L21 = 0 и M2 1 = 0; и, следовательно, L1 = IL1 и M 1 = I M1: Первый шаг алгоритма 7 ( k = 1 ) вычисляет a11 h(1)a 11 a 11g(1) R(2) = A(2) L1 M1 D1; (2.9) что содержит всю нужную информацию. 416 Г.Л. Литвинов, А. Я. Родионов, С. Н. Сергеев, А. Н. Соболевский Теперь мы продолжим с подматрицы R(2) матрицы A(2); раскладывая ее на мно- жители по аналогии с (2.7), и т. д. Опишем формально k -й шаг этой конструкции, отвечающий k -му шагу алгоритма 7. На этом общем шаге мы работаем с A(k) = O(k 1)(k 1) O(k 1)(n k+1) O(n k+1)(k 1) R(k) ; (2.10) где R(k) = h(k 1)(a(k 1) k 1;k 1)g(k 1) B(k 1) = a(k) kk h(k) g(k) B(k) ! : (2.11) Как и на первом шаге, мы представляем (A(k)) = M kD k(A(k+1))L k; (2.12) где Mk = 0 @ O(k 1)(k 1) O(k 1)1 O(k 1)(n k) O1(k 1) O11 h(k)(a(k) kk ) O(n k)(k 1) O(n k)1 O(n k)(n k) 1 A; Dk = 0 @ O(k 1)(k 1) O(k 1)1 O(k 1)(n k) O1(k 1) a(k) kk O1(n k) O(n k)(k 1) O(n k)1 O(n k)(n k) 1 A; Lk = 0 @ O(k 1)(k 1) O(k 1)1 O(k 1)(n k) O1(k 1) O11 O1(n k) O(n k)(k 1) (a(k) kk )g(k) O(n k)(n k) 1 A; A(k+1) = Okk Ok(n k) O(n k)k R(k+1) ; R(k+1) = h(k)(a(k) kk )g(k) B(k): (2.13) Заметим, что имеется следующее рекуррентное соотношение на элементы A(k) : a(k+1) ij = ( 0; если i 6 k или j 6 k; a(k) ij a(k) ik (a(k) kk )a(k) kj ; если i > k и j > k: (2.14) Эта рекурсия сразу видна в алгоритме 7. Более того, можно показать по индукции, что матрица, вычисленная на k -м шаге алгоритма, равна A(k+1) Mk i=1 Li Mk i=1 Mi Mk i=1 Di: (2.15) Иначе говоря, эта матрица составлена из h(1)a 11; ..., h(k)(a(k) kk ) (в верхнем треуголь- нике), a 11g(1); ..., (a(k) kk )g(k) (в нижнем треугольнике), a11; : : : ; a(k) kk (на диагонали) и R(k+1) (в юго-восточном углу). УНИВЕРСАЛЬНЫЕ АЛГОРИТМЫ 417 После раскрытия всех выражений (2.12) для A(k); где k = 1; : : : ; n; мы получаем A = M 1D 1 M nD nL n L 1 (2.16) (на самом деле, Mn = Ln = 0 и, следовательно, M n = L n = I ). Заметив, что D i и M j коммутируют при i < j; мы можем переписать A = M 1 M nD 1 D nL n L 1: (2.17) Рассмотрим равенства (D1 : : : Dn) = D 1 D n; (L1 : : : Ln) = L n L 1; (M1 : : : Mn) = M 1 M n: (2.18) Первое из этих равенств очевидно. По поводу остальных двух заметим, что M2 k = L2 k = 0 для всех k; следовательно, M k = I Mk и L k = I Lk: Далее, LiLj = 0 при i < j и MiMj = 0 при i > j: Используя эти равенства, можно показать, что (L1 : : : Ln) = Mn 1 i=0 (L1 : : : Ln)i = = (I Ln) (I L1) = L n L 1; (M1 : : : Mn) = Mn 1 i=0 (M1 : : : Mn)i = = (I M1) (I Mn) = M 1 M n; (2.19) что дает два последних равенства из (2.18). Заметим, что в (2.19) мы использовали нильпотентность L1 : : : Ln и M1 : : : Mn; что позволяет применить (1.6). Можно проверить, что матрицы M := M1 : : : Mn; L := L1 : : : Ln и D := D1 : : : Dn содержатся в верхнем треугольнике, в нижнем треугольнике и, соответственно, на диагонали матрицы, вычисленной алгоритмом 7. Эти матрицы удовлетворяют LDM -разложению A = MDL: Этим завершается объяснение ал- горитма 7. В терминах матричных вычислений алгоритм 7 есть версия LDM -разложения с внешним произведением. Этот алгоритм может быть переделан в форму, почти иден- тичную с [1, алгоритм 4.1.1]: Алгоритм 8. LDM-разложение (версия 2). Вход: n n -матрица A с элементами a(i; j); также используется для хранения результата вычислений и промежуточных результатов вычислений. for j = 1 : n v(1 : j) = a(1 : j; j) for k = 1 : j 1 418 Г.Л. Литвинов, А. Я. Родионов, С. Н. Сергеев, А. Н. Соболевский v(k + 1 : j) = v(k + 1 : j) a(k + 1 : j; k) v(k) end for i = 1 : j 1 a(i; j) = (a(i; i)) v(i) end a(j; j) = v(j) for k = 1 : j 1 a(j + 1 : n; j) = a(j + 1 : n; j) a(j + 1 : n; k) v(k) end d = (v(j)) a(j + 1 : n; j) = a(j + 1 : n; j) d end Этот алгоритм выполняет в точности те же операции, что и алгоритм 7, вычисляя по- следовательно один за другим столбцы результата. А именно, в первой части основного цикла он вычисляет элементы a(i) ij для i = 1; : : : ; j; сначала под видом элементов v и окончательно в присваивании ѕ a(i; j) = (a(i; i)) v(i) ї. Сложность этого алгоритма такая же, как у алгоритма 7. 2.4. LDM -разложение с симметрией и ленточной структурой Когда матрица A симметрична, т. е. aij = aji для всех i; j; естественно ожидать, что LDM -разложение тоже должно быть симметрично, т. е. M = LT : Действительно, повторяя рассуждения предыдущего раздела, можно показать, что все промежуточные матрицы A(k) симметричны, следовательно, Mk = LTk для всех k и M = LT : Те- перь мы представим две версии симметричного LDM -разложения, отвечающие двум версиям LDM -разложения, данным в предыдущем разделе. Заметим, что объем вы- числений в этих алгоритмах примерно вдвое меньше, чем в их полных версиях. В обоих случаях они требуют n3=6+O(n2) операций и (каждый) и n 1 операций взятия алгебраического замыкания. Алгоритм 9. Симметричное LDM -разложение (версия 1). Вход: n n -матрица A с элементами a(i; j); также используется для хранения результата вычислений и промежуточных результатов вычислений. for j = 1 : n 1 v(j) = (a(j; j)) for k = j + 1 : n for l = j + 1 : k a(k; l) = a(k; l) a(k; j) v(j) a(l; j) end end a(j + 1 : n; j) = a(j + 1 : n; j) v(j) end УНИВЕРСАЛЬНЫЕ АЛГОРИТМЫ 419 Строго треугольная матрица L содержит нижний треугольник результата, а матрица D на диагонали. Следующая версия обобщает [1, алгоритм 4.1.2]. Алгоритм 10. Симметричное LDM-разложение (версия 2). Вход: n n -матрица A с элементами a(i; j); также используется для хранения результата вычислений и промежуточных результатов вычислений. for j = 1 : n for i = 1 : j 1 v(i) = (a(i; i)) 1a(i; j) end v(j) = v(j) a(j; 1 : j 1) v(1 : j 1) a(j; j) = v(j) for k = 1 : j 1 a(j + 1 : n; j) = a(j + 1 : n; j) a(j + 1 : n; k) v(k) end d = (v(j)) a(j + 1 : n; j) = a(j + 1 : n; j) d end Заметим, что эта версия требует обратимости замыканий a(i; i); вычисленных алго- ритмом. Замечание. В случае идемпотентного полукольца мы имеем (D)2 = D; следо- вательно, A = (MD)(DL): Когда A симметрична, мы можем записать A = (G)TG; где G = DL: Очевидно, эта идемпотентное разложение Хо- лецкого приводит лишь к небольшой модификации алгоритма 9 (или 10). См. также классическое разложение Холецкого в [1, алгоритм 4.2.2]. A = (aij) называется ленточной матрицей с верхней шириной полосы q и нижней шириной полосы p; если aij = 0 для всех j > i+p и всех i > j+p: Ленточная матрица с p = q = 1 называется тридиагональной. Для обоснования универсальной версии LDM -разложения ленточных матриц нам нужно показать, что параметры лент матриц A = A(2); : : : ;A(k); вычисленных в процессе LDM -разложения, не больше параметров матрицы A(1) = A: Положим по индукции, что A = A(1); : : : ;A(k) имеют требуемые параметры ленты и рассмотрим элементы a(k+1) ij для i > j + p: Если i 6 k или j 6 k; то a(k+1) ij = 0; поэтому мы можем предположить, что i > k и j > k: В этом случае i > k + p; следовательно, a(k) ik = 0 и a(k+1) ij = a(k) ij a(k) ik (a(k) kk )a(k) kj = 0: Таким образом, мы показали, что нижняя полоса A(k) не больше p: Таким же образом можно показать, что ее верхняя полоса не превышает q: Мы используем это для построения ленточной версии LDM -разложения, еј прототип см. в [1, алгоритм 4.3.1]. 420 Г.Л. Литвинов, А. Я. Родионов, С. Н. Сергеев, А. Н. Соболевский Алгоритм 11. LDM-разложение ленточной матрицы. A ленточная n n -матрица с элементами a(i; j); нижняя ширина ленты p и верх- няя ширина ленты q также используется для хранения окончательного результата и промежуточных результатов вычислений. for j = 1 : n 1 v(j) = (a(j; j)) for i = j + 1 : min(j + p; n) a(i; j) = a(i; j) v(j) end for k = j + 1 : min(j + q; n) for i = j + 1 : min(j + p; n) a(k; j) = a(k; j) a(k; i) a(i; j) end end for k = j + 1 : min(j + q; n) a(j; k) = v(j) a(j; k) end end Когда p и q фиксированы и n >> p; q переменная, легко видеть, что алгоритм вы- полняет примерно npq каждой из операций и : Замечание. Есть некоторые еще более специальные виды ленточных матриц, на- пример, матрицы Хессенберга и трјхдиагональные матрицы. Матрицы Хессенберга определяются как ленточные матрицы с p = 1 и q = n; а в случае трјхдиагональ- ных матриц p = q = 1: Несложно записать дальнейшие адаптации алгоритма 11 для этих случаев. 2.5. Итерационные схемы По настоящему универсальную итерационную схему трудно найти, так как условия, при которых такие схемы работают, а также условия остановки для таких схем зависят от полукольца и от представления данных. Тем не менее, основные операции таких схем бывают универсальны и очень просты. Рассмотрим следующий итерационный процесс решения уравнения Беллмана X = AX B : X(k+1) = AX(k) B: (2.1) Итерируя выражения (2.1) для всех k от 0 до m 1; мы получаем X(m) = AmX(0) mM 1 i=0 AiB: (2.2) Таким образом, результат зависит от поведения AmX(0): Алгоритм может быть переписан следующим образом (для случая, когда B вектор-столбец). Алгоритм 12. Итерации Якоби Вход: n n -матрица A с элементами a(i; j) ; n 1 вектор-столбцы b и x УНИВЕРСАЛЬНЫЕ АЛГОРИТМЫ 421 situation=’вычисление’ while situation==’вычисление’ x = A x b situation=newsituation(...) if situation==’расходимость’ disp(’Итерации Якоби расходятся.’) exit end if situation==’сходимость’ disp(’Решение найдено:’) exit end end Выход: ситуация, x: Далее мы кратко обсудим поведение итераций Якоби над обычной арифметикой с неот- рицательными вещественными числами и над полукольцом Rmax: Для простоты в обо- их случаях мы ограничимся случаем неприводимой матрицы A; то есть когда ассоци- ированный граф сильно связан. Хорошо известно, что над обычной арифметикой (в неприводимом неотрицательном случае) итерации Якоби сходятся тогда и только тогда, когда наибольшее собственное значение r(A) матрицы A строго меньше 1: Это связано, в частности, с тем, как ведет себя последовательность fAmx(0)gm>0: Заметим, что итерации Якоби, даже в случае сходимости, как правило, дают лишь приближјнное решение уравнения x = Ax + b: В случае Rmax поведение последовательности fAmx(0)gm>0 отличается от случая обычной неотрицательной алгебры. Однако, как и в случае неотрицательной алгеб- ры, это поведение зависит от r(A); наибольшего собственного значения A в алгебре макс-плюс (то есть относительно собственной задачи A x = x в этой алгебре). А именно Amx(0) ! 0 и, следовательно, итерации сходятся если r(A) < 1: Более того, в этом случае A = (I A: : :An 1) и, следовательно, итерации доставляют точное решение уравнения Беллмана после конечного числа шагов. Наоборот, Amx(0) ! +1 и, следовательно, итерации расходятся при r(A) > 1: Подробности см., например, [18]. На границе r(A) = 1 степени Am достигают периодического режима после ко- нечного числа шагов. Следовательно, последовательность fAb Amx(0)gm>0 также становится периодической после конечного числа шагов. Если период равен единице, то есть если эта последовательность стабилизируется, то метод сходится к некоторому общему решению уравнения x = Axb; которое описывается как суперпозиция Ab и некоторого собственного вектора матрицы A [2, 17]. Вектор Ab может доминировать, в этом случае метод сходится к Ab как ѕожидалосьї. Однако, период последователь- ности (после выхода на периодический режим) может быть и больше единицы, в этом случае итерации Якоби не доставляют никакого решения уравнения x = Ax b непо- средственно, однако такое решение может быть найдено с помощью суммирования по периоду. Степени матриц и теория собственных векторов и собственных значений над алгеброй макс-плюс подробно описаны в монографии Бутковича [16], см. также [49, 50]. 422 Г.Л. Литвинов, А. Я. Родионов, С. Н. Сергеев, А. Н. Соболевский В более сложной схеме итераций Гаусса-Зейделя мы также можем использовать предварительно найденные координаты X(k): В этом случае матрица A записывается в виде L U; где L строго нижнетреугольная часть A; а U верхнетреугольная часть с диагональю. Итерации записываются как X(k) = LX(k) UX(k 1) B = LUX(k 1) LB: (2.3) Заметим, что преобразования правой части однозначно, так как L строго нижне- треугольная матрица, и L однозначно определяется как I L : : : Ln 1 (где n размерность A). Другими словами, мы только применяем формальную подстановку. Итерируя выражения (2.3) для всех k вплоть до m; мы получим X(m) = (LU)mX(0) mM 1 i=0 (LU)iLB: (2.4) Правая часть напоминает формулу (L U) = (LU)L; см. (1.1), поэтому естествен- но ожидать, что эти итерации сходятся к AB при хорошем выборе X(0): Результат зависит от поведения (LU)mX(0): Алгоритм может быть записан следующим образом (мы снова полагаем, что B вектор-столбец). Алгоритм 13. Итерации Гаусса-Зейделя Вход: n n -матрица A с элементами a(i; j) ; n 1 -вектор-столбцы b и x situation=’вычисление’ while situation==’вычисление’ for i = 1 : n y(i) = a(i; i : n) x(i : n) b(i) end for i = 2 : n x(i) = a(i; 1 : i 1) x(1 : i 1) ‘ end situation=newsituation(...) if situation==’расходимость’ disp(’Итерации Гаусса-Зейделя расходятся.’) exit end if situation==’сходимость’ disp(’Решение найдено:’) exit end end Выход: ситуация, x: Анализ работы схемы Гаусса-Зейделя в случаях тропической (макс-плюс) алгебры и неотрицательной линейной алгебры мы здесь не приводим, однако ясно, что его можно провести по аналогии с итерациями Якоби. УНИВЕРСАЛЬНЫЕ АЛГОРИТМЫ 423 2.6. LU -разложение Обсудим более кратко универсальную версию разложения A = LU; которая бы- ла описана в работах Б. А. Карре [18]. Этот раздел содержит лишь формулировки соответствующих алгоритмов, так как доказательства были даны в работах Карре и Бекхауса [14, 18, 19] и могут быть получены (как и сами алгоритмы) путем небольшой модификации соответствующих рассуждений в разделе про LDM -разложение (см. вы- ше). Напомним, что в случае классической системы AX = B над числовыми полями разложение A = LU; где L нижняя унитреугольная (с диагональными элементами равными 1 ), а U верхняя треугольная матрица, приводит к следующему разложе- нию исходного уравнения: LZ = B; UX = Z: Таким образом, мы получаем A 1 = U 1L 1 при условии, что матрица A обратима. Достаточно найти L и U; так как после этого линейная система (2.6.) может быть решена комбинацией прямой подстановки для Z и обратной подстановки для X: В случае системы X = AX B над полукольцами, следуя [18], мы можем записать Z = LZ B; X = UX Z: Таким образом, A = UL; где U верхне-треугольная, а L строго нижне-треугольная матрица. Далее мы приводим алгоритмы прямой и обратной подстановки, а также две версии универсального LU -разложения (то есть одновременного вычисления матриц L и U ). Алгоритм 14. Прямая подстановка. Вход: Строго нижнетреугольная n n -матрица l ; n 1 вектор b: y = b for k = 2 : n y(k) = l(k; 1 : k 1) y(1 : k 1) end Выход: вектор y: Алгоритм 15. Обратная подстановка. Вход: Верхнетреугольная n n -матрица u ; n 1 вектор b: y = b y(n) = (u(n; n)) y(n) for k = n 1 : 1 : 1 424 Г.Л. Литвинов, А. Я. Родионов, С. Н. Сергеев, А. Н. Соболевский y(k) = u(k; k + 1 : n) y(k + 1 : n) y(k) = (u(k; k)) y(k) end Выход: вектор y: Отметим, что обратная подстановка отличается от случая LDM -разложения, так как матрица U (в отличие от M ) является лишь треугольной, а не строго треугольной. Прямые подстановки для обеих схем, разумеется, идентичны. Алгоритм 16. LU -разложение (версия 1). Вход: n n -матрица A с элементами a(i; j); также используется для хранения ре- зультата вычислений и промежуточных результатов вычислений. for j = 1 : n 1 v(j) = (a(j; j)) a(j + 1 : n; j) = a(j + 1 : n; j) v(j) a(j + 1 : n; j + 1 : n) = a(j + 1 : n; j + 1 : n) a(j + 1 : n; j) a(j; j + 1 : n) end Алгоритм 17. LU -разложение (версия 2). Вход: n n -матрица A с элементами a(i; j); также используется для хранения ре- зультата вычислений и промежуточных результатов вычислений. for j = 1 : n v(1 : j) = a(1 : j; j) for k = 1 : j 1 v(k + 1 : j) = v(k + 1 : j) a(k + 1 : j; k) v(k) end a(1 : j; j) = v(1 : j) for k = 1 : j 1 a(j + 1 : n; j) = a(j + 1 : n; j) a(j + 1 : n; k) v(k) end d = (v(j)) a(j + 1 : n; j) = a(j + 1 : n; j) d end Эти алгоритмы основаны на следующих определениях, аналогичных (2.12) и (2.13). Полагая A(1) := A = a11 h(1) g(1) B(1) ; для любого k = 1; : : : ; n получаем (A(k)) = U k (A(k+1))L k; (2.1) УНИВЕРСАЛЬНЫЕ АЛГОРИТМЫ 425 где Uk = 0 @ O(k 1)(k 1) O(k 1)1 O(k 1)(n k) O1(k 1) a(k) kk h(k) O(n k)(k 1) O(n k)1 O(n k)(n k) 1 A; Lk = 0 @ O(k 1)(k 1) O(k 1)1 O(k 1)(n k) O1(k 1) O11 O1(n k) O(n k)(k 1) (a(k) kk )g(k) O(n k)(n k) 1 A; A(k+1) = Okk Ok(n k) O(n k)k R(k+1) ; R(k+1) = h(k)(a(k) kk )g(k) B(k) = a(k+1) k+1;k+1 h(k+1) g(k+1) B(k+1) ! : (2.2) 2.7. Программная реализация универсальных алгоритмов Программные реализации универсальных алгоритмов на полукольцах не могут быть такими же эффективными, как и стандартные программные реализации (относитель- но скорости вычислений), но они намного более гибкие. Программные модули могут иметь дело с абстрактными (и переменными) операциями и типами данных. Конкрет- ные значения этих операций и типов данных могут быть заданы соответствующими входными данными. В этом случае конкретные операции и типы данных порождаются дополнительными программными модулями. В программах, написанных в таком духе, удобно использовать специальные методы так называемого объектно-ориентированного (и функционального) дизайна, см., например, [41, 47, 52]. К счастью, недавно появились мощные инструменты, поддерживающие объектно-ориентированный дизайн программ, включающие компиляторы для реальных и удобных языков программирования (на- пример, C++ и Java) и современные системы компьютерной алгебры. С недавних пор, этот тип программных методов был назван обобщенным программированием (см., на- пример, [47]). Реализация на C++. Используя шаблоны и объектно-ориентированное программи- рование, А. В. Чуркин и С. Н. Сергеев [9] создали программу (на Visual C++), которая демонстрирует, как некоторые универсальные алгоритмы вычисляют замыкания мат- риц A и решают уравнения Беллмана x = Axb в различных полукольцах. Програм- ма также может решать обычную систему Ax = b в обычной арифметике, преобразо- вывая ее в форму ѕБеллманаї. Перед нажатием кнопки ѕSolveї пользователь должен выбрать полукольцо, задачу и алгоритм. После этого начальные данные записывают- ся в матрице (для удобства отображения размерность матрицы ограничена 10 10 ). Результат может быть и вектором, и матрицей в зависимости от решаемой задачи. Объектно-ориентированный подход позволяет реализовать различные полукольца как объекты с различными определениями основных операций, тогда как сама программа алгоритма записана компактно и в единственном ѕэкземпляреї (шаблоне). Примеры полуколец. Полукольцо соответствует объекту, используемому универсаль- ным алгоритмом, и определяет конкретную реализацию этого алгоритма. Пользователь может выбрать одно из следующих полуколец: 426 Г.Л. Литвинов, А. Я. Родионов, С. Н. Сергеев, А. Н. Соболевский 1) = + и = : обычная арифметика над действительными числами; 2) = max и = +: арифметика макс-плюс над R [ f 1g; 3) = min и = +: арифметика мин-плюс над R [ f+1g; 4) = max и = : арифметика макс-умножить над неотрицательными числами; 5) = max и = min : арифметика макс-мин над вещественным отрезком [a; b] (концы a и b может выбирать пользователь); 6) =OR и =AND: булева логика над двухэлементном множеством f0; 1g: Алгоритмы. Пользователь может выбрать один из следующих методов: 1) Метод исключения Гаусса, включающий универсальные реализации эскала- торного метода (алгоритм 1), Флойда-Уоршалла (алгоритм 2, алгоритм Ершова (основанный на прототипе из [11, гл. 2]) и универсальный алгоритм Ротэ [48]; 2) Методы для систем Тјплица, включающие универсальные реализации схем Дурбина и Левинсона (алгоритмы 3 и 4); 3) LDM -разложение (алгоритм 7) и его адаптации к симметрическим (ал- горитм 9), ленточным (алгоритм 11), трјхдиагональным матрицам и матрицам Хессенберга. 4) Итерационные схемы Якоби и Гаусса-Зейделя. Как было упомянуто выше, эти схемы не универсальны, т. к. критерии остановки и условия, при которых они сходятся и дают правильный ответ, различны для обычной арифметики и идем- потентных полуколец (в частности, макс-плюс). Типы матриц. Пользователь может выбрать работу с общими матрицами или с мат- рицами специального вида, например, симметричные, симметричные тјплицевы, лент- ночные, Хессенберга и тридиагональные. Отображение. В случае идемпотентных полуколец матрица может быть отображе- на как взвешенный направленный граф. После выполнения вычислений пользователь может захотеть найти оптимальный путь между данной парой вершин или отобразить дерево оптимальных путей. Эти задачи могут быть решены, используя родительские ссылки, как в случае с классическим методом Флойда-Уоршалла, вычисляющим все оптимальные пути. См., например, [7]. В нашем случае механизм родительских ссылок может быть реализован прямо в классе, описывающем идемпотентную арифметику. Другие арифметики и интервальные расширения. Возможна также реализация раз- личные типов арифметики как типов данных и их сочетание с выбором полуколец. Более того, все реализованные полукольца могут быть расширены до их интерваль- ных версий. Такие возможности не были реализованы в программе А. В. Чуркина и С. Н. Сергеева [9], они отложены до следующей версии. Список этих арифметик вклю- чает целые числа, дробную арифметику с цепными дробями и управляемой точностью. Реализации в MATLAB. Вся работа (кроме инструментов отображения) была по- вторена в MATLAB [9], который также обладает возможностями объектно- ориентированного программирования. Разумеется, универсальные алгоритмы, запи- санные в MATLAB, очень близки к описанным в настоящей статье. УНИВЕРСАЛЬНЫЕ АЛГОРИТМЫ 427 Перспективы на будущее. По-видимому, основная задача заключается в том, что- бы, используя идемпотентный принцип соответствия и вычисления над полукольцами, разработать программную систему широкого назначения, основанную на коллекции универсальных алгоритмов. Разработка такой системы на основе универсальных алго- ритмов может гарантировать сокращение рабочего времени для программистов и поль- зователей благодаря унификации программного обеспечения. Также может быть гаран- тирована произвольная точность и безопасность численных вычислений. При этом надо учесть, что высокоуровневые инструменты, такие как STL [47, 52], обладают и очевид- ными преимуществами, и некоторыми недостатками, поэтому должны использоваться с осторожностью. Система должна содержать несколько уровней (включая уровни пользователя и программиста) и большое количество модулей. Грубо говоря, она должна быть раз- делена на три части. Первая часть должна содержать модули, которые реализуют модули доменов (конечные представления основных математических объектов). Вто- рая часть должна реализовывать универсальные (инвариантные) методы вычисления. Третья часть должна содержать модули, реализующие алгоритмы, которые зависят от модели. Эти модули могут быть использованы в пользовательских программах на C++, Java, Maple, MATLAB и др. В частности, система должна содержать следующие модули: Доменные модули: - целые числа произвольной длины; - рациональные числа; - рациональные числа конечной точности (см. [8]); - комплексные рациональные числа конечной точности; - рациональные числа с фиксированной и плавающей запятой; - комплексные рациональные числа; - вещественные числа с плавающей запятой произвольной точности; - комплексные числа произвольной точности; - p -адические числа; - интервальные числа; - кольца многочленов над различными кольцами; - идемпотентные полукольца; - интервальные идемпотентные полукольца; - и другие. Алгоритмы: - линейная алгебра; - численное интегрирование; - корни многочленов; - интерполяция сплайнами и аппроксимации; - рациональные и полиномиальные интерполяции и аппроксимации; - вычисление специальных функций; - дифференциальные уравнения; - оптимизации и оптимальное управление; 428 Г. Л. Литвинов, А. Я. Родионов, С. Н. Сергеев, А. Н. Соболевский - идемпотентный функциональный анализ; - и другие. Эта программная система может быть особенно полезна разработчикам алгоритмов, программистам, студентам и математикам.

About the authors

Grigory L. Litvinov

Institute for Information Transmission Problems of the Russian Academy of Sciences

Email: glitvinov@gmail.com
Doctor of Physics and Mathematics, Professor 19 Bolshoy Karetny Pereulok, Moscow 127051, Russian Federation

Аnatoliy Ya. Rodionov

Moscow Center for Continuous Mathematical Education

Email: ayarodionov@yahoo.com
Teacher 11 Bolshoy Vasil’evsky Pereulok, Moscow 119002, Russian Federation

Serge˘ı Sergeev

University of Birmingham, School of Mathematics

Email: sergeevs@maths.bham.ac.uk
Associate Professor of the Mathematics School Watson Building, Birmingham B15 2TT, United Kingdom

Andrei N. Sobolevsky

Institute for Information Transmission Problems of the Russian Academy of Sciences

Email: sobolevski@iitp.ru
Doctor of Physics and Mathematics, Professor of RAS, Director of the Institute for Information Transmission Problems of the Russian Academy of Sciences (Kharkevich Institute) 19 Bolshoy Karetny Pereulok, Moscow 127051, Russian Federation

References

Supplementary files