CC BY
65
Вычислительные технологии
Том 8, № 1, 2003
ОПТИМАЛЬНОЕ ВНЕШНЕЕ ОЦЕНИВАНИЕ МНОЖЕСТВ РЕШЕНИЙ ИНТЕРВАЛЬНЫХ СИСТЕМ УРАВНЕНИЙ
Часть 2
С. П. ШАРЫй Институт вычислительных технологий СО РАН Новосибирск, Россия e-mail: shary@ict.nsc.ru
The work is devoted to developing parameter partitioning methods (PPS-methods) for optimal (exact) outer component-wise estimation of the solution sets to interval linear equations systems. The results of computational experiments and comparisons with the other known approaches to the problem are presented.
Введение
Настоящая работа является продолжением статьи [1] и посвящена задаче внешнего оценивания объединенного множества решений интервальных систем линейных алгебраических уравнений (ИСЛАУ) вида
anxi + ai2X2 + ... + ai ,nx,n = bi,
a2ixi + a22x2 + • • • + a2nxn = b2,
. .. . (1)
anixi + an2X2 + • • • + annxn = bn
с интервалами a^- и b или в краткой форме
Ax = b (2)
с квадратной интервальной n х n-матрицей A = (a^-) и интервальным n-вектором правой части b = (bi).
Напомним, что объединенное множество решений интервальной линейной системы Ax = b — это множество
EUni(A, b) := { x G Rn | (3A G A)(3b G b)( Ax = b) }, (( С.П. Шарый, 2003.
образованное всеми решениями точечных систем Ax = b с A £ A и b £ b. Нас интересует вычисление по возможности более точных внешних покоординатных оценок этого множества решений:
Для данной интервальной линейной системы уравнений Ax = b
найти оценки для величин min{ xv | x £ Suni(A, b) } снизу и (3)
для величин max{ xv | x £ Suni(A, b) } сверху, v = 1, 2,... ,n.
В постановке (3) можно ограничиться требованием вычисления только минимума min{ xv | x £ Suni(A, b) }, поскольку
max{ xv | x £ Euni(A, b) } = - min{ xv | x £ Euni(A, -b) }
для всех v = 1, 2,... , n. Ниже мы для краткости иногда будем называть задачу (3) внешней задачей для интервальной линейной системы Ax = b.
Одним из основных итогов нашей предшествующей статьи [1] является идея алгоритмов, использующих для уточнения оценок множеств решений интервальных систем уравнений адаптивное дробление интервалов их параметров. От красивой идеи путь к эффективным и практичным алгоритмам идет, как правило, через кропотливую работу по модификации исходной схемы и ее "скрещиванию" с более мелкими усовершенствованиями, обеспечивающими успех на реальных вычислительных задачах. Цель этой публикации — представить работу по практическому усовершенствованию методов дробления параметров в применении к интервальным линейным системам уравнений, приведшую к созданию лучших в своем классе алгоритмов для оптимального внешнего оценивания множеств решений ИСЛАУ.
1. Модификации методов дробления параметров
Приступая к построению более совершенных методов дробления параметров для задачи внешнего оценивания объединенного множества решений ИСЛАУ, мы возьмем алгоритм из табл. 2 работы [1] как основу, которая будет развиваться и дополняться рядом усовершенствований, уже стандартных для интервальных методов глобальной оптимизации, базирующихся на стратегии "ветвей и границ". Как правило, их перечень (не претендующий на полноту) включает в себя следующее (см., например, [2-7]):
1. Строят более качественную внешнюю оценивающую функцию (интервальное расширение) для целевой функции покоординатной оценки.
2. Выявлением монотонности целевой функции по тем или иным переменным на интервалах из рабочего списка L добиваются уменьшения размерности этих интервалов.
3. На основе специфических локальных свойств целевой функции в соответствующих интервалах применяют более эффективные, чем бисекция, процедуры минимизации (например, методы градиентного спуска там, где целевая функция гладкая и выпуклая).
4. Наряду с оцениванием целевой функции по целым интервалам вычисляют ее значения в каких-то точках этих интервалов, — эти значения доставляют верхнюю границу
искомого глобального минимума, и ее знание позволяет чистить рабочий список £ от тех записей, которые заведомо не могут быть ведущими.
Необходимость решительной модификации простейшего метода дробления параметров диктуется, в частности, его большой информационной сложностью: фактически каждый шаг алгоритма сопровождается ростом рабочего списка £, так что начиная с некоторого момента быстрая оперативная память ЭВМ может исчерпаться, а следующий за этим обмен с внешними носителями может стать фактором, который замедляет на порядки производительность процессора (см. § 2).
1.1. Тест на монотонность
Пусть дана интервальная система линейных алгебраических уравнений (х = г и нам известны
1 ( # | > . , . ,
дхи (Q, r) дхи (Q, r)
и
dqij
интервальные расширения соответствующих производных
дх^ (Q,r)
dq
и
dri вод
дх^ (Q,r)
j
дт-i
от V-й компоненты вектора решения системы уравнений Qx = г по г]-му элементу матрицы Q и ¿-му элементу вектора г. Если интервальные п х п-матрица Ш = (с^) и п-вектор Г = (г образованы из элементов
qij
[ q., q.] при
—ij —ij
[ qij, qij] пРи
qij
при int
дху (Q, r) дqij
дхи (Q, r) дц^
дхи(Q, r)
дц
> 0,
< 0,
0,
(4)
ij
[ ri, rj при [ ri, ri] при
дх„(Q, r) дт-i
дх„ (Q, r) дг,.
0,
0,
. дху (Q, r)
при int --- Э 0,
дт-i
(5)
где "int" обозначает внутренность интервала, то, очевидно,
min{ х„ | х £ Su„i(Q, r) } = min{ хи | х е Suni(Q, r) }.
А поскольку количество существенно интервальных (с ненулевой шириной) элементов в Q и r может быть, вообще говоря, значительно меньшим, чем в Q и r, то, переходя от
r
исходной ИСЛАУ Qx = г к системе Qx = г, мы тем самым добьемся упрощения задачи вычисления шт{ xV | х £ Surai(Q, г) }.
Как найти фигурирующие в (4), (5) интервальные расширения производных? Это можно сделать, например, следующим образом. Из курса математического анализа известно, что У = (у^) — обратная матрица для 3 = (^), поэтому производные решения вещественной линейной системы Qx = г по ее коэффициентам даются формулами
дхV
dqij dri
(см., к примеру, [8, 9]). Следовательно, в случае, когда Y = (yij) — "обратная интервальная матрица" для Q, т. е. внешняя интервальная оценка для множества обратных матриц
Y D{ Q-1 | Q е Q },
а Xj — j-я компонента некоторого интервального вектора x D Surai (Q, r), то мы можем принять следующие интервальные оценки производных:
дх^(Q, r) дх^(Q, r)
= -yviXj, = yvi- (6)
Вычисление "обратной интервальной матрицы" можно выполнить как решение задачи внешнего оценивания объединенного множества решений интервального матричного уравнения
AY = I,
где I — единичная матрица, применив n раз тот же самый метод внешнего оценивания End, который выбран базовым для всего алгоритма. Иногда для вычисления Y и x уместно использовать какие-нибудь дешевые приближенные алгоритмы (вроде метода Хансена [10] для локализации "обратной интервальной матрицы").
1.2. Стратегия дробления
Традиционно в интервальных методах глобальной оптимизации, которые основаны на схеме "ветвей и границ" и являются ближайшими родственниками методов дробления параметров, ведущие интервалы на каждом шаге дробятся по своей самой длинной грани. Этот способ дробления гарантирует, как известно (см. работы [7, 11]), сходимость использующих его алгоритмов, и мы тоже применили его в простейшем методе дробления параметров, хотя для него подобные проблемы не возникают из-за его конечности.
После того как собственно сходимость достигнута, на первый план выходит желание оптимизации метода, т.е. обеспечения наибольшей быстроты сходимости. В свою очередь, эту задачу далее обычно редуцируют к упрощенной постановке — обеспечить наиболее значительное улучшение целевой функции на каждом шаге алгоритма. Строгая и точная оптимизация методов дробления параметров даже в таком виде трудна и вряд ли целесообразна в полном объеме. Мы будем решать эту задачу, руководствуясь разумными эвристическими соображениями, на основе знания оценок производных целевой функции.
Если матрицы <3 = () и 3 = () отличаются друг от друга лишь в единственном (к,1)-м элементе, < , то в силу теоремы Лагранжа о среднем имеем
(3-1г)и - (3-1г)и = [ дк1,дк1 ]
для некоторой матрицы < £ [ <, <]. Аналогично, если векторы г = (г^) и г = (Г) отличаются только в единственной к-й компоненте и Гк < г к, то
(д-1^ - (<-1г)^ = дх;(д,г~^ [ гк , г ]
дгк
для некоторого вектора 5 € [г, г].
Предположим теперь, что интервальные матрицы (( и (( получены из интервальной матрицы ( заменой элемента чы на его левую ч и правую ^кг границы:
ч кг = кг = чкг.
Допустим также, что
шт{ Xv | х £ 2«Ы(С(, г) } и шт{ Xv | х £ ^«„¿((С, г) }
достигаются на одном и том же семействе концов интервальных элементов матрицы и вектора правой части, что почти всегда имеет место для "достаточно узких" интервальных систем. Тогда
шт{ х^ | х £ ~«т(С, г) } - шт{ XV | х £ , г) } = дхгу(<, г) ^
дды
для некоторой матрицы < £ ( и вектора г £ г. Сходным образом предположим, что Г и Г — это интервальные векторы, полученные из интервального вектора г заменой его к-й компоненты на левую гк и правую гк границы:
Гк = Тк, гк = гк •
При прежнем условии, что
шт{ Xv | х £ г) }, шт{ Xv | х £ г) }
достигаются на одном и том же множестве концов интервальных элементов матрицы и вектора правой части ИСЛАУ, снова получаем
дxv (<, г)
шт{ Xv | х £ г) } - шт{ Xv | х £ г) } = —^—-— wid гк
дгк
V
для некоторых матрицы < £ ( и вектора г £ г. Следовательно, величина произведения ширины (или радиуса) интервального элемента на модуль интервального расширения соответствующей производной решения может служить в некотором смысле локальной мерой того, как бисекция элементов из ( либо г влияет на шт{ xv | х £ ^„„¿(С, г) } и размеры объединенного множества решений ИСЛАУ.
Далее, оценки объединенного множества решений ИСЛАУ, получаемые с помощью большинства существующих методов, являются тем более точными, чем меньше размеры этого множества решений. Д. Гей, например, в [12] доказал для своих методов оценку погрешности, пропорциональную ширине интервальной оболочки множества решений. Аналогичная оценка доказана А. Ноймайером для метода Кравчика [13]. С подобными базовыми методами уменьшение размеров множества решений г) должно приво-
дить к аналогичному и сравнимому по величине изменению в оценке Т((, г). При этом,
следовательно, требование наиболее быстрого улучшения целевой функции за один шаг метода дробления параметров становится по существу эквивалентным условию наиболее быстрого уменьшения размеров множества решений при дроблении ведущей ИСЛАУ.
Учитывая сделанные выше (частично эвристические) заключения, мы рекомендуем поэтому рассекать ведущие ИСЛАУ по элементам, на которых достигается максимальная из величин
öxv(Q, r)
dq
wid q
öxv (Q, r)
ÖTi
wid г,
(7)
г,] £ { 1, 2,...,п }, т.е. по элементам, на которых достигается максимум произведения оценки производной решения на ширину интервала, а не просто по самым широким элементам ИСЛАУ.1
1.3. Модификация Рона
Результат следствия 2 из теоремы 2.3 предшествующей работы [1], который мы использовали с целью вывода метода дробления параметров, для объединенных множеств решений ИСЛАУ был значительно усилен в 80-е годы И. Роном, уточнившим множество вершин матрицы A и вектора правой части b, на которых достигаются минимальное и максимальное значения компонент точек из Suni(A, b) [17].
Чтобы дать математически строгую формулировку результата И. Рона, нам потребуются некоторые дополнительные обозначения. Пусть
E := { x G Rn | | x,| = 1 для всех i = 1, 2,... , n }
— множество n-векторов с компонентами ± 1. Для данной интервальной матрицы A и фиксированных о, т G E обозначим через ACTT = (aj) матрицу размера n х n, образованную элементами
!aij, если о7-т,- = — 1, , если a,Tj = 1.
Кроме того, через Ьа = (Щ) мы будем обозначать n-вектор, составленный из элементов
!b,, если о, = 1, b,, если о, = -1.
Матрица ACTT и вектор Ьа образованы, таким образом, наборами концов элементов матрицы A и вектора b соответственно, и всего имются 2n • 2n = 4n различных матрично-векторных пар вида (ACTT ,Ьа) для о и т, независимо пробегающих множество E.
Оказывается [17], что для невырожденной матрицы A как минимальное, так и максимальное значения компонент точек множества решений Suni(A, b) могут достигаться лишь на множестве 4n матриц ACTT и ассоциированных с ними векторов Ьст, т.е.
min{ xv | x G Su„i(A, b) } = min f( ACTT)-V) , (8)
ct,T ££ V / v
тах{ XV | X £ ^„¿(А, Ь) } = шаЛ (АСТТ) V) (9)
<Г,Т е£ V /V
хЗа рубежом имеется тенденция связывать такой способ выбора рассекаемой компоненты с именем Д. Ратца, рассматривавшего его в своих работах [14, 15]. Мы не следуем этому, поскольку независимо от Д. Ратца и даже раньше такая стратегия дробления, требующая максимизации величин (7), была предложена автором в статье [16].
для каждого индекса V = 1, 2,... ,п. Рассмотрим, как можно использовать этот факт в наших методах дробления параметров.
Важно осознавать, что приведенный выше результат не накладывает никаких ограничений на концы отдельно взятого интервального элемента матрицы А либо вектора правой части Ь, если в рассмотрение не привлечена информация о концах других элементов. Ограничение на ту или иную комбинацию концов интервалов, следующее из (8), (9), является существенно коллективным, и, чтобы принять его во внимание, мы должны отслеживать концы задействованных интервальных элементов по всей матрице А и всему вектору Ь. Для практической реализации этих идей с каждой интервальной линейной системой Qж = г, Q = (qij), г = (гг), порожденной в процессе дробления исходной ИСЛАУ (1), (2), мы связываем:
— вспомогательную п х п-матрицу Ш = (), элементы которой равны ±1 или 0, такую что
( -1, если = ,
0, если = ,
[ 1, если = '
— вспомогательные п-векторы в = (вг) и £ = &) с компонентами
для всех г, = 1, 2,... , п, и
( -1, если гг = Ьг,
вг := < 0, если гг = Ьг,
(10)
1, если гг = Ьг.
Кроме того, рабочий список £ будет состоять теперь из расширенных записей
Q, г, Х^, г), Ш, 5, £ ), (11)
имеющих шесть полей за счет добавления трех полей для хранения Ш, в и которые получаются с предыдущих шагов алгоритма.
Векторам в и £ предстоит быть "приближениями" соответственно к а и т, входящим в равенства (8), (9), в то время как Ш = в£т — это "приближение" к матрице (атт). В начале работы алгоритма мы инициализируем все элементы Ш, в и £ нулями, а затем перевычисляем их с целью заменить нулевые элементы (соответствующие неопределенности) на ненулевые (соответствующие тому или иному концу интервала). Матрица Ш и векторы в,£, взаимно влияя друг на друга и перевычисляясь в процессе работы алгоритма, предназначены, следовательно, для наблюдения над процессом дробления исходной интервальной линейной системы и его корректировки таким образом, чтобы порождать лишь варианты, допускаемые равенствами (8), (9). Это реализуется следующим образом.
На каждом шаге алгоритма при разбиении интервального элемента Ь ведущей ИСЛАУ Qж = г мы смотрим на соответствующее значение:
— матрицы Ш = (), если элемент Ь есть qы из матрицы Q;
— вектора в = (вг), если элемент Ь есть из вектора правой части г.
В случае Wki = 0 (sk = 0 соответственно) мы порождаем, согласно обычной процедуре бисекции из простейшего метода дробления параметров табл. 2 работы [1], две интервальных системы-потомка — Q'x = г' и Q''x = г'', которые отвечают обоим концам рассекаемого интервала h. Но если Wki = 0 (sk = 0 соответственно), то после бисекции ведущей ИСЛАУ мы можем оставить в рабочем списке только одного потомка от Qx = г. Какого именно, зависит от знака элемента Wki (sk соответственно). Говоря более точно, мы выполняем процедуру, представленную в табл. 1.
Почему это в принципе возможно? Другими словами, не может ли отбрасывание второй интервальной системы-потомка в вышеописанной процедуре нарушить то свойство ведущих оценок Y(Q, г), что они приближают искомый min{ xv | x £ Suni(A, b) } снизу?
Таблица 1
Порождение систем-потомков
IF ( рассекается qki ) THEN IF ( Wki = 0 ) THEN
порождаем системы Q'x = г' и Q''x = г'' так, что
qj := qij := qij для (i,j) = (М), qki := q^ qk'i := q^ г' :=г'' :=г;
ELSE
порождаем систему Q'x = г' так, что
г' :=г, qij := qij для (i,j) = (к,
q' .= \ qki для wki = 1, qki = \ — 1 { qki для wki = -1;
END IF END IF
IF ( рассекается Tk ) THEN IF ( sk = 0 ) THEN
порождаем системы Q'x = г' и Q''x = г'' так, что
Q' := Q'' := Q,
^ := гi для i = := Iki, г/k' := Tk,
ELSE
порождаем систему Q'x = г' так, что Q' := Q, гi := г для i = k,
' l !k для Sk = -1, г! := '
!k для Sk = 1;
END IF END IF
Для ответа на этот вопрос заметим, что в новой процедуре дробления из табл. 1 мы отбрасываем лишь те системы, которые:
— не принадлежат множеству точечных систем { Латх = ba | а,т £ E };
— не содержат систем из этого множества.
Следовательно, в силу свойства (C1) базисного оценивающего метода и, принимая во внимание равенство (8), мы получаем
min{ х„ | х £ Evni(A, b) } = min ((Лстт )-V) > minYM^ ,ba) >
<г,тe£ V 4 ' ) v <г,тe£
> min{ T(Q, r) | Q Э Лат и r Э ba для некоторых а,т £ E} >
> min{ T(Q, r) | система Qx = r присутствует в списке L} = = ведущая оценка Y(Q, r).
Таким образом, с новой процедурой дробления ведущая оценка действительно приближает min{ xv | х £ Svni(A, b) } снизу.
Обращение к описанию формальной вычислительной схемы "модификации Рона" мы начнем с установления правил перевычисления W, s и t в процессе работы алгоритма. Существует взаимно-однозначное соответствие между вектором s и правой частью интервальной системы Qx = r, в то время как дробление интервальной матрицы Q рассматриваемой ИСЛАУ влияет на векторы s и t лишь неявно — посредством матрицы W и условий (10). Последнее тем не менее позволяет организовать перевычисление W, s и t на каждом таком шаге алгоритма, который завершается дроблением интервального элемента ведущей системы на два конца. В противном случае, если ведущая интервальная система порождает только одного потомка Q^ = r' в соответствии с табл. 1, векторы s, t и матрица W остаются неизменными, так что
s'
t' := t, W' := W.
Итак, пусть ведущая интервальная система (х = г рассечена на два ИСЛАУ-потомка — Ш'х = г' и С'х = г", определенные в табл. 1. Каким должен быть закон, в соответствии с которым формируются матрицы Ш', Ш'' и векторы в', в'', Ь', Ь'', соответствующие системам-потомкам? Сначала мы можем положить
Ш'' := Ш' := в'' := в' := в, Ь'' := Ь' := Ь,
а затем выполнить следующую двухэтапную процедуру перевычисления: Во-первых, модифицируем Ш', Ш'' и в', в'' на основе информации о только что выполненной бисекции ведущей ИСЛАУ. Именно:
(I) если рассеченным элементом был Скг из матрицы Ш, то в матрицах Ш' = () и Ш'' = () мы полагаем
т'ы := 1 и т'ы := -1;
(II) если рассеченным элементом был г к из вектора правой части г,
то в векторах в' = (в[) и в'' = (в") мы полагаем
' ''
вк := —1 и вк := 1.
Во-вторых, перевычисляем каждое из двух семейств взаимосвязанных объектов, — Ш', в', Ь' и Ш'', в'', Ь'' соответственно, — используя соотношения (10), а именно:
s
(I) если s' или t' изменился, пытаемся перевычислить матрицу W'; (II) если W' или t' изменился, пытаемся перевычислить вектор s'; (III) если W' или s' изменился, пытаемся перевычислить вектор t'.
Инструкции (I) - (III) повторяются последовательно одна за другой в цикле до тех пор, пока изменения в W', s' и t' не прекратятся. Та же самая процедура применяется затем к W'', s'', t''. Полная алгоритмическая схема приведенной выше процедуры оказывается довольно сложной, так что имеет смысл снабдить читателя более строгим и детальным ее описанием. Соответствующий псевдокод представляет табл. 2. Ниже мы даем объяснения к нему. Булевские переменные
W' Change, s'Change, t'Change, W'C, s'C, t'C
и
W ''Change, s''Change, t''Change, W''C, s''C, t''C
— это "флаги", введенные с целью отображения текущего статуса изменений в W', s', t' и W'', s'', t'' соответственно. Значение true показывает, что объект, с которым связан флаг, был изменен на текущей итерации процесса перевычисления, тогда как значение false — "без изменений".
Для того чтобы завершить формальное описание модификации Рона, нам нужно лишь определить, что имеется в виду под "пытаемся перевычислить матрицу W' ", "пытаемся перевычислить вектор s' " и тому подобными инструкциями из табл. 2.
Обозначим прописными греческими буквами K', Л' и П' подмножества индексов элементов вектора s', вектора t' и матрицы W' соответственно, которые изменили свои значения (с 0 на ±1) на текущем шаге процедуры перевычисления табл. 2. K' и Л' — это, следовательно, подмножества множества первых n натуральных чисел
{ 1, 2,... , n },
тогда как П есть подмножество множества всех пар, составленных из первых n натуральных чисел, т. е. подмножество множества
{ (1,1),(1,2),..., (1,n),(2,1),...,(n,n) }.
Каждое из множеств K', Л', П' может быть пустым, но может содержать и более чем один элемент. Тогда наша "попытка перевычислить вектор s' " может быть организована следующим образом (табл. 3).
"Попытка перевычислить вектор t'" может быть выполнена аналогично тому, как это произведено выше с тем единственным отличием, что цикл "DO FOR l £ Л'" во втором IF-операторе должен быть заменен на "DO FOR k £ K'". Далее приводится процедура перевычисления W' (табл. 4).
То же следует произвести для перевычисления s'', t'' и W''. При этом нам необходимо организовать индексные подмножества K'', Л'' и П'' для представления индексов компонент векторов s'', t'' и матрицы W'' соответственно, которые изменились на текущем шаге перевычислительного процесса.
С. П. Шарый Таблица 2
Перевычисление W', W'', s', s", t', t''.
W'Change : = false; s'Change : = false; if Change : = false; W"Change := false; s''Change : = false; t"Change := false; IF ( ( рассекается qkl из Q ) AND ( qkl рассекается на два конца ) ) THEN
m'kl := 1; m' := -1; W'Change := true; W''Change := true; END IF
IF ( ( рассекается rk из r ) AND ( rk рассекается на два конца ) ) THEN
s'k := -1; s'k := 1; s'Change := true; s''Change := true; END IF
DO WHILE ( W 'Change OR s'Change OR if Change ) IF ( s'Change OR t'Change ) THEN
пытаемся перевычислить матрицу W'; IF ( W' изменилась ) THEN W'C := true ELSE W'C := false END IF END IF
IF ( W'Change OR t'Change ) THEN
пытаемся перевычислить вектор s'; IF ( s' изменился ) THEN s'C := true ELSE s'C := false END IF END IF
IF ( W'Change OR s'Change ) THEN
пытаемся перевычислить вектор t'; IF ( t' изменился ) THEN t'C := true ELSE t'C := false END IF END IF
W'Change := W'C; s'Change := s'C; t'Change := t'C; END DO
DO WHILE W''Change OR s''Change OR t''Change IF ( s''Change OR t''Change ) THEN
пытаемся перевычислить матрицу W''; IF ( W'' изменилась ) THEN W''C := true ELSE W''C := false END IF END IF
IF ( W''Change OR t''Change ) THEN
пытаемся перевычислить вектор s''; IF ( s'' изменился ) THEN s''C := true ELSE s''C := false END IF END IF
IF ( W''Change OR s''Change ) THEN
пытаемся перевычислить вектор t''; IF ( t'' изменился ) THEN t''C := true ELSE t''C := false END IF END IF
W''Change := W''C; s''Change := s''C; t''Change := t'C; END DO
Таблица 3
Перевычисление в'
1Р ( Ш'ОНаиде ) ТИЕЫ БОРОН (к, £ П'
1р ( ¿1 =0 ) вк := 'ш'ыМ ЕЫБ БО ЕЫБ 1Р
1Р ( ¿'ОНаиде ) ТИЕЫ БО РОИ I £ Л'
БО РОИ к = 1 ТО и
1Р ( вк = 0 АШ w'kl = 0 ) вк := /¿1 ЕЫБ БО ЕЫБ БО ЕЫБ 1Р
Таблица 4
Перевычисление Ш'
1Р ( в'ОНаиде ) ТИЕЫ БО РОИ к £ К'
БО РОИ I = 1 ТО и
1Р ( ¿1 = 0 ) wkl := вк¿1
ЕЫБ БО ЕЫБ БО ЕЫБ 1Р
1Р ( ¿'ОНаиде ) ТИЕЫ БО РОИ I £ Л'
БО РОИ к = 1 ТО и
1Р ( вк = 0 ) wkl := вк¿1 ЕЫБ БО ЕЫБ БО ЕЫБ 1Р
Наконец, мы должны упомянуть следующее замечательное свойство матрицы Ш: в каждой ее 2 х 2-подматрице любой элемент равен произведению остальных трех элементов. Чтобы убедиться в этом, обозначим через г1,г2 номера строк и через ]1, ]2 — номера столбцов, образующих рассматриваемую подматрицу. Тогда в соответствии с определением (10)
.1 ТЛ , ^¿1 Т.?2 , ^2Л а»2 ТЛ , ^2.2 а»2 Т.?2 •
Перемножая любые три из вышеприведенных равенств, к примеру, первое, второе и четвертое, получим
^¿1.1 ^¿1.2 ^¿2.2
аг1 ап Т.2 ^¿2 Т.2 •
Квадрат любой компоненты векторов а и т есть 1, так что
^¿1.1 ^¿1.2 ^¿2.2 = ТЛ а»2 = ^¿2.1 •
(12)
Таблица 5 Уточнение W посредством поиска 2 х 2-подматриц
DO FOR i = 1 TO n
DO FOR j = 1 TO n
IF ( i = k AND j = l ) THEN
IF ( wij = 0 AND wii = 0 AND wkj = 0 ) THEN Wki := WijWiiWkj; EXIT END IF END IF END DO END DO
То же самое с остальными элементами подматрицы:
Wiljl Wilj2 Wi2 jl Wilj2 Wi2 jl Wi2j2 = Wil jl Wi2 jl Wi2 j2 =
Иногда соотношения (12) - (15) могут оказаться полезными для дальнейшего уточнения матрицы W. Например, пусть мы намереваемся рассечь ведущую интервальную систему Qx = г по элементу qki, в то время как соответствующий элемент Wki матрицы W равен нулю. При этом согласно обычному правилу дробления мы должны были бы породить две системы-потомка
вместо ИСЛАУ Qx = г. Но разумно сделать попытку доопределить Wki, поискав 2 х 2-подматрицы в W, имеющие ненулевыми все элементы, за исключением Wki. Если такая подматрица в W найдется, то мы присваиваем элементу Wki значение произведения остальных трех элементов. Сказанное может быть реализовано в виде следующего алгоритма из табл. 5, где оператор EXIT означает выход из всех блоков и циклов выписанного псевдокода.
1.4. Отсев бесперспективных записей
Рассмотрим, наконец, модификацию простейшего метода дробления параметров, связанную с вычислением оценок Y(mid Q, mid г) для "середин" ведущих систем. Она позволит контролировать точность вычисления значения min{ xv | x £ Suni(A, b),}, а также удалять из списка L бесперспективные записи, т. е. такие, которые никогда не сделаются ведущими. Благодаря последнему обстоятельству несколько ограничивается рост длины рабочего списка L.
Действительно, пусть наряду с оценкой Y(Q, г) для интервальных систем Qx = г, порождаемых алгоритмом, вычисляются еще и величины Y(HQ, Иг), где символом "□" обозначена операция взятия какой-то фиксированной точки из интервала. Очевидно, что
Y(HQ, Иг) > Y(Q, г),
и значения Y(HQ, Иг) приближают min{ x | x £ Suni(A, b) } сверху. Если для каждого шага алгоритма определять величину
и = minYpQ, Иг) (16)
W
i2 j2
W
il jl
W
il j2
(13)
(14)
(15)
по всем интервальным линейным системам Qx = r, когда-либо побывавшим в списке L до этого шага, то
min{ ж^ | x E Su„i(A, b) } < w. С другой стороны, если Qx = r — ведущая ИСЛАУ, то
Y(Q, r) < min{ Xv | x E Su„i(A, b) },
и потому еще одним критерием естественной остановки алгоритма может служить достижение требуемой малости разности (w — Y(Q, r)) для ведущей оценки Y(Q, r).
Далее, интервальная система Qx = r, удовлетворяющая на некотором шаге условию
Y(Q, r) >w, (17)
при дальнейшем выполнении алгоритма уже никогда не станет ведущей, и удаление соответствующей записи из рабочего списка L не окажет никакого влияния на результат работы алгоритма. Вообще говоря, тестироваться неравенством (12) должны на каждом шаге алгоритма все вновь заносимые в список L записи, тогда как "чистки" списка — его просмотр и удаление всех удовлетворяющих (12) записей — имеет смысл проводить лишь после изменений (т. е. уменьшений) величины w.
Конечно, было бы идеальным выбирать (HQ, Иг) E Arg min{ Y(Q,r) | Q E Q,r E r }. Фактически мы так и будем поступать, если оценка Y(Q, r) — точная, полагая
Y(HQ, Hr) = Y(Q, r).
Но в общем случае подобный "удачный" выбор не менее труден, чем решение исходной задачи, и потому в целях минимизации возможных отклонений (HQ, Иг) от множества Argmin{ Y(Q,r) | Q E Q,r E r } целесообразней всего брать тогда в качестве HQ и Иг середины матрицы Q и вектора r, т.е. mid Q и mid r соответственно.
1.5. Влияние базового алгоритма
Обсудим теперь очень важный вопрос о влиянии базового метода Enal и, следовательно, способа получения оценки Y(Q, r) на вычислительную схему конкретных реализаций методов дробления параметров.
Для работы многих методов внешнего оценивания объединенного множества решений ИСЛАУ требуются некоторые начальные приближения, содержащие оцениваемое множество решений. Таковы, например, интервальный метод Гаусса — Зейделя, методы Кравчика, Гея и др. Нетрудно понять, что интервальный вектор внешней оценки объединенного множества решений ИСЛАУ Qx = r может служить начальным приближением для процедуры внешнего оценивания объединенных множеств решений систем-потомков Q'x = r' с Q' С Q и г' С г. То же самое относится и к вычислению "обратной интервальной матрицы", которая нужна в тесте на монотонность из § 1.1 и при выборе рассекаемого интервального элемента из § 1.2.
Имеет смысл поэтому при выборе соответствующих базовых методов хранить интервальные векторы внешней оценки множеств решений и "обратную интервальную матрицу" с предыдущих шагов алгоритма. Для этого потребуется расширить образующие рабочий список L записи еще на два поля, так что теперь в методе с модификацией Рона мы будем оперировать не тройками, как в простейшем методе дробления параметров, или шестерками (11), а восьмерками
(Q, r, Y(Q, r),W,s,t, Y, х), (18)
т. е. записями, состоящими из восьми членов, где
Q — интервальная n х n-матрица, Q С A;
r — интервальный n-вектор, r С b;
T(Q, r) — нижний конец v-й компоненты (для заданного фиксированного номера v) внешней интервальной оценки множества решений Suni(Q, r);
W, s, t — вспомогательные величины, необходимые для реализации модификации Рона и определенные в § 1.3;
Y — интервальная n х n-матрица, такая что Y D { Q-1 | Q G Q };
x =5 "«гаг (Q, r).
Другое соображение. Как правило, любая из методик решения внешней задачи для ИСЛАУ, удовлетворяющая условию (C2) из [1], дает точную оценку Y(Q, r) не только тогда, когда Q и r — точечные (неинтервальные), но и на некотором более широком множестве начальных данных, когда часть элементов в Q и r могут быть интервалами. Например, интервальный метод Гаусса — Зейделя [6, 13], методы Гея [12] и другие итерационные методы, основанные на теореме о сжимающем отображении [8], обеспечивают точную оценку Y(Q, r) в случае вещественной матрицы Q, а вектор правых частей r может быть при этом любым. Получаемая по интервальному методу Гаусса оценка Y(Q, r), очевидно, точна для верхних треугольных матриц Q и т. д. Возможны и более хитроумные условия на взаимное расположение элементов интервальной матрицы Q, их знак, ширину и т. п. Как правило, выявление всех практически наиболее значимых подобных ситуаций нетрудно алгоритмизовать.
Таким образом, для естественной остановки метода дробления параметров совсем не обязательно дожидаться полного "овеществления" ведущей ИСЛАУ Qx = r (условие завершения цикла "DO WHILE" в табл. 2 работы [1]). Достаточно того, чтобы ведущая оценка Y(Q,r) являлась точной.
Можно пойти еще дальше и предусмотреть возможность динамической смены в алгоритме базового метода End. Первоначально это может быть какой-нибудь метод с широкой сферой применимости, но и с большой погрешностью внешнего оценивания. Затем, по мере достижения в процессе дробления заданной узости интервальных систем-потомков, он заменяется на более точный специализированный алгоритм.
Вывод, вытекающий из вышеизложенного, прежний: для достижения наибольшей эффективности методов дробления параметров все составляющие реального алгоритма, применяемого к той или иной конкретной задаче, как-то: структура данных (т. е. записей в рабочем списке L); способ их обработки; способ дробления ведущих ИСЛАУ; другие характеристики алгоритма — должны быть тщательно увязаны с особенностями решаемой задачи. Общая схема методов дробления параметров предоставляет большую свободу разработчику, но нужно умело ей пользоваться.
1.6. Итоговая схема алгоритма
Приводимые ниже псевдокоды в табл. 6 и 7 суммируют развитые в предшествующих пунктах работы модификации методов дробления параметров (PPS-методов) для внешнего оценивания объединенного множества решений ИСЛАУ.
Теоретически модификация Рона позволяет уменьшить экспоненциальный множитель в верхней оценке трудоемкости PPS-методов с 2n +n до 4n, но необходимо отметить следующее:
— это сделано ценой существенного утяжеления алгоритма, так что его программная сложность становится весьма высокой;
— при решении больших практических задач размерности, превосходящей несколько десятков, как 2n +n, так и 4n делаются недостижимыми величинами, а методы дробления параметров следует рассматривать скорее как итерационную уточняющую процедуру, которая никогда не выполняется до своего естественного завершения.
На наш взгляд, в каждом конкретном случае пользователь, основываясь на информации о размере задачи, ее структуре, наличных вычислительных ресурсах, самостоятельно должен решать вопрос о целесообразности включения модификации Рона в ту или иную реализацию метода дробления параметров. По этой причине мы приводим ниже две итоговые вычислительные схемы методов дробления параметров: с модификацией Рона и без нее.
В табл. 6 и 7 предполагается, что в качестве базового метода внешнего оценивания End взят интервальный метод Гаусса — Зейделя, но это сделано лишь для определенности. Подчеркнем, что методы дробления параметров — это общая схема, а табл. 6 и 7 представляют лишь два ее возможных варианта. Развитые нами выше конструкции содержат ряд "свободных переменных", которые должны быть настраиваемы под конкретную задачу, так что мы можем, следовательно, говорить о целом классе методов, основанных на общей идее дробления интервальных параметров системы.
Алгоритм, приведенный в табл. 6, описывает метод дробления параметров без модификации Рона. Его рабочий список L образуется пятичленными записями вида
( Q, r, Y(Q, r), Y, x ).
Для начала работы этого алгоритма, который мы называем LinPPS1, нам нужно:
— найти грубые внешние оценки для объединенного множества решений исходной интервальной системы и "обратной интервальной матрицы", т. е. вычислить x D ^«ni(A, b) и Y D A-1;
— назначить точность e > 0;
— положить Y(A, b) := x и и :=
— инициализировать рабочий список L записью (A, b, x, Y, x).
Алгоритм, приведенный в табл. 7, названный нами LinPPS2, представляет метод дробления параметров с модификацией Рона, который оперирует с восьмичленными записями вида (18). Для начала его работы необходимо выполнить первые три пункта, как для алгоритма LinPPS1, далее положить
W := 0, s := 0, t := 0
для матрицы W и векторов s и t, введенных в § 2.3, и инициализировать рабочий список L записью (A, b, x, W, s,t, Y, x).
В заключение полезно проследить, в каких отношениях построенные нами алгоритмы находятся к другим методам решения ИСЛАУ. Прежде всего, отметим, что имеется глубокая связь (фактически являющаяся двойственностью) между развиваемыми в настоящей работе методами дробления параметров и предложенным автором ранее классом методов дробления решений [11, 18]. Но к идее построения методов дробления параметров для интервальных линейных систем можно прийти и несколько иным путем, доводя до логического завершения некоторые из давно и широко известных в интервальном анализе подходов к решению ИСЛАУ.
С. П. Шарый Таблица 6
Алгоритм LinPPSl
Ш1ЬЕ ( ( ведущая оценка T(Q, г) неточна ) ИЛИ (и — T(Q, г) > е) ) по формулам (6) вычисляем интервальные расширения для производных дж^(Q, г) дж^г)
дг '
соответствующие элементам qij и гг с ненулевой шириной;
"сжимаем" в соответствии с (4), (5) интервальные элементы в Q и г, на которых была выявлена монотонная зависимость от ^ и гг, и полученные в результате этой процедуры интервальную матрицу и интервальный вектор также обозначаем через Q и г;
ищем среди элементов ведущей ИСЛАУ Qж = г интервальный элемент б, которому соответствует наибольшее из произведений
dxv(Q, г)
ij
• wid qij,
dxy (Q, г) dri
wid гi, i, j £ { 1, 2,... , n };
порождаем интервальные системы-потомки Q'x = г' и Q''x = г": если s = qki для некоторых k, / £ { 1, 2,...,n }, то полагаем
qij := qij := qij для j) = (k0, qki := qki, qki :=
г' := г'' := г;
если s = ^ для некоторого k £{ 1, 2,...,n }, то полагаем г := гi' := гi для г = k, 4 := !k, ^ := rk,
Q' := Q'' := Q;
вычисляем интервальные векторы x' = End (Q', г') и x'' = End (Q'', г''); вычисляем оценки Y(Q', г') и Y(Q'', г'');
вычисляем внешние оценки для "обратных интервальных матриц" Y' Э (Q')-1 и Y' Э (Q')-1;
вычисляем оценки Y(mid Q', mid г') и Y(mid Q'', mid г'') и полагаем ^ := min{ Y(mid Q', mid г'), Y(mid Q'', mid г'') };
удаляем из L бывшую ведущей запись (Q, г, Y(Q, г),р);
если Y(Q', г') < и, то заносим запись (Q', г', Y(Q', г'), Y', x') в список L в порядке возрастания значений третьего поля;
если Y(Q'', г'') < и, то заносим запись (Q'', г'', Y(Q'', г''), Y'', x'') в список L в порядке возрастания значений третьего поля;
если и > то полагаем и := ^ и чистим список L:
удаляем из него все такие записи (Q, г, Y(Q, г), Y, x), что Y(Q, г) > и;
Таблица 7
Алгоритм ЫпРРБ2
DO WHILE ^ ( ведущая оценка Y(Q, r) неточна ) OR (и — Y(Q, r) > e) ^
используя формулы (6), вычисляем интервальные расширения производных дж^(Q, r) дж^(Q, r)
dqij dri '
которые соответствуют интервальным элементам qj и r с ненулевой шириной;
"сжимаем" в соответствии с (4), (5) в Q и r интервальные элементы , для которых выделена монотонность зависимости от qj и ri? и обозначаем полученные при этом интервальные матрицу и вектор снова через Q и r;
находим среди элементов системы Qx = r интервал h, который соответствует наибольшему из произведений
дж^(Q, r)
dq
wid qij,
джу(Q, r) dri
■ wid r^ g{ 1, 2,...,n };
пытаемся уточнить матрицу W в соответствии с процедурой табл. 5;
"дробим" элемент h и порождаем одну или две интервальные системы-потомки Q'x = r' и Q''x = r" в соответствии с процедурой табл. 1;
если были порождены две системы-потомка, вычисляем новые матрицы W', W'' и векторы s', s'', t', t'' в соответствии с процедурами табл. 2-4; иначе полагаем W' := W, s' := s, t' := t;
вычисляем интервальные вектор x' = End (Q', r') и, возможно, x'' = End (Q'', r''), беря x в качестве начального приближения;
присваиваем оценку и' := Y(Q', r') и, возможно, и'' := Y(Q'', r'');
уточняем внешнюю оценку для "обратной интервальной матрицы" Y' D (Q')-1 и, возможно, для Y'' D (Q'')-1, беря Y в качестве начального приближения;
вычисляем оценку Y(mid Q', mid r') и, возможно, Y(mid Q'', mid r'') и полагаем ^ := min{ Y(mid Q', mid r'), Y(mid Q'', mid r'') };
удаляем бывшую ведущую запись (Q, r, и, Y, W, s, t, x) из списка L;
если и' < и, то помещаем запись (Q', r', и', W', s', t', Y', x') в список L в порядке возрастания третьего поля;
если из ведущей ИСЛАУ были порождены две системы-потомка и и'' < и, то помещаем запись (Q'', r'',u'', W'', s'', t'', Y'', x'') в список L в порядке возрастания третьего поля;
если и > то полагаем и := ^ и чистим список L, т. е.
удаляем из него все такие записи (Q, r, и, W, s, t, Y, x), что и > и;
Одной из прародительниц алгоритма ЫпРРБ1 можно, по-видимому, считать процедуру Купермана — Хансена (см. [8]), в которой внешние интервальные оценки для объединенного множества решений ИСЛАУ ищутся на основе пассивного однократного использования информации о производных компонент решения по элементам матрицы и правой части системы, подобно тому как это делается в § 1.1. Но в реальных интервальных системах не всегда удается выявить определенный знак этих производных, а потому какие-то интервальные элементы все равно останутся интервальными и после применения теста на монотонность. И Куперман, и Хансен в этом месте останавливались и завершали процесс решения "внешней задачи". Дальнейшее уточнение искомой оценки на этом пути возможно лишь при добавлении активной процедуры измельчения элементов матрицы ИСЛАУ, хотя рассматривать будет нужно все получающиеся при таком измельчении-дроблении интервальные системы.
Фактически это и реализуется в алгоритме ЫпРРЭ1. Список £ хранит все возникшие в процессе работы варианты (за исключением тех заведомо бесперспективных, которые выявляются посредством теста из § 1.4), а обрабатывается на каждом шаге алгоритма "самый перспективный" вариант, если мерой этой перспективности считать значение оценки
т г).
2. Сравнения и численные эксперименты
Как мы уже упоминали выше, задача внешнего оценивания объединенного множества решений ИСЛАУ является на сегодня одной из наиболее популярных и наиболее разработанных интервальных постановок задач. За сорок лет, прошедших с момента ее формулировки, были предложены четыре принципиально различных подхода к вычислению оптимальных (точных) оценок объединенного множества решений для общих интервальных линейных систем. Один из них восходит к работе У. Оеттли [19], который первым обнаружил, что пересечение объединенного множества решений с ортантами пространства Кп является выпуклым полиэдром. Таким образом, точное значение шт{ хи | х £ ^ипг } может быть найдено путем решения в каждом из ортантов некоторой задачи линейного программирования и последующим взятием минимального из результатов. Этот подход, развитие которого рассмотрено, например, Х. Янссоном в [20], и который мы распространили в [1] на задачи оценивания обобщенных множеств решений, основывается на пассивной переборной стратегии, а трудоемкость его экспоненциально растет в зависимости от размерности п. Поэтому практическая значимость его очень ограничена.
Следующие два вычислительных подхода к оптимальному решению "внешней задачи" для интервальных линейных систем — это методы дробления решений и методы дробления параметров (называемые также РББ-методами и РРБ-методами), предложенные автором (см., например, работы [11, 16, 18, 21]). Хотя для п х п-системы сложность выполнения методов дробления решений в худшем случае пропорциональна 2п, а верхняя оценка сложности выполнения методов дробления параметров равна даже 2п +п, эти методы, имея в основе стратегию "ветвей и границ", являются адаптивными. При исполнении каждого последующего шага как в методах дробления решений, так и в методах дробления параметров мы полностью используем информацию о результатах предыдущих шагов.
Наконец, четвертый и в настоящий момент наиболее популярный в зарубежной литературе подход к оптимальному решению "внешней задачи" предложен И. Роном [17] (см. также [13, 22]). Отталкиваясь от характеризации Оеттли — Прагера для объединенно-
го множества решений, И. Рон показывает, что для случая квадратной невырожденной матрицы A искомые значения min{ xv | x £ Suni(A, b) } и max{ xv | x £ Suni(A, b) }, v = 1, 2,... , n, достигаются на множестве не более чем 2n решений так называемого уравнения Оеттли — Прагера
| mid A ■ x — mid b | = rad A ■ |x| + rad b. (19)
Вычисляя все возможные решения этого уравнения и сравнивая их между собой, мы получим оптимальные оценки множества решений за конечное число шагов. Поскольку процесс определения каждого последующего решения для (19) никак не зависит от решений, уже найденных раньше, алгоритм Рона в целом не является адаптивным (т. е. подобен переборным методикам), в то время как его трудоемкость пропорциональна 4n в худшем случае и 2n в лучшем.
Таким образом, все подходы, разработанные к настоящему моменту для вычисления оптимальных оценок объединенного множества решений интервальных линейных систем, имеют экспоненциальную трудоемкость в наихудшем случае. Как уже отмечалось, этот факт не является следствием "плохости" предложенных алгоритмов, а отражает глубокие свойства самого объединенного и других множеств решений интервальных линейных систем: задача вычисления его оптимальных покоординатных оценок оказалась NP-трудной [23-28]. Следовательно, экспоненциальная сложность всех перечисленных выше алгоритмов существенна и не может быть устранена [29].
Каковы же преимущества и недостатки методов дробления решений и методов дробления параметров в сравнении с другими подходами для нахождения оптимальных решений "внешней задачи"? Наиболее важная их особенность состоит в том, что они порождают последовательности приближенных оценок искомых величин "с нужной стороны", т. е. для min{ xv | x £ } снизу, а для max{ xv | x £ } сверху. Именно такие оценки и требуются в соответствии со смыслом "внешней задачи". Процесс выполнения метода дробления параметров, например, разбивается на ряд эффективно вычислимых этапов-шагов, в результате каждого из которых мы получаем некоторое решение "внешней задачи". После того как в таком алгоритме проработал хотя бы один из этих этапов, его прерывание в любой момент приведет к тому, что мы все равно будем иметь в своем распоряжении некоторое решение "внешней задачи". Иными словами, если у нас имеются достаточные вычислительные мощности, то, реализуя методы дробления параметров (как и методы дробления решений), мы можем быть вполне уверены, что некоторый ответ к задаче будет наверняка получен, хотя, возможно, и не оптимальный. Следовательно, как методы дробления решений, так и методы дробления параметров являются последовательно гарантирующими. С учетом труднорешаемости "внешней задачи" для интервальных линейных систем это свойство радикальным образом выделяет методы дробления решений и методы дробления параметров из всех других подходов для вычисления оптимального решения.
Напротив, подходы Оеттли и Рона к нахождению оптимальных решений "внешней задачи", имея экспоненциальную в худшем случае трудоемкость, дают желаемые "внешние" оценки объединенного множества решений лишь в финале, при естественном завершении своей работы, поскольку раньше нельзя гарантировать то, что вычисленная оценка в действительности не превосходит min{ xv | x £ } (соответственно, не меньше max{ xv | x £ }). Следовательно, эти алгоритмы являются финально гарантирующими. Если размерность интервальной линейной системы достаточно велика (больше
нескольких десятков), то в силу труднорешаемости "внешней задачи" количество арифметических и логических операций, необходимое для того, чтобы задача была наверняка решена, начинает превосходить количество операций, выполнимое на сколь угодно мощном компьютере за любое разумное время (час, день, год или даже столетие). В этих условиях нельзя быть уверенным, что финально гарантирующий алгоритм, будучи примененным к задаче, вообще завершит свою работу и, таким образом, будет получено решение поставленной задачи. Иначе говоря, применяя финально гарантирующий алгоритм, мы рискуем попусту растратить время и ресурсы без получения хоть какого-то ответа к нашей задаче.
Этот пессимистичный прогноз особенно справедлив для пассивных переборных алгоритмов. Если адаптивные (последовательные) алгоритмы достигают экспоненциальной трудоемкости лишь на наихудших вариантах, а в среднем для "не очень плохих задач" работают с приемлемыми трудозатратами, то, напротив, пассивные переборные алгоритмы, какими являются подходы Оеттли — Прагера и Рона, требуют экспоненциальных трудозатрат всегда, вне зависимости от характера решаемой задачи.
Отметим, что алгоритмы из [17, 19] для нахождения оптимальных решений "внешней задачи" — все-таки последовательно гарантирующие, но относительно так называемого "слабого внутреннего" оценивания. Фактически эти алгоритмы решают не "внешнюю", а некоторую другую задачу для интервальной линейной системы, требующую оценивания min{ xv | x G Suni(A, b) } сверху и max{ xv | x G Suni(A, b) } снизу.
Сказанное мы проиллюстрируем результатами вычислительных экспериментов на рабочей станции Sun ULTRA-10 (тактовая частота процессора 300 МГц, шины — 100 МГц, емкость ОЗУ 128 Мбайт), которые были проведены с несколькими версиями методов дробления параметров, реализованными на языке программирования Fortran фирмы Sun Microsystems.2
Для того чтобы продемонстрировать влияние различных факторов на поведение методов дробления параметров для внешнего оценивания объединенного множества решений ИСЛАУ, мы приводим в табл. 8 результаты тестовых расчетов с четырьмя методами:
A — простейший метод дробления параметров (т. е. алгоритм из табл. 2 работы [1]), но дополненный процедурой удаления из рабочего списка L бесперспективных записей;
B — тот же алгоритм, что и A, но дополненный процедурой сжатия элементов матрицы ИСЛАУ на основе информации о монотонности оценки (см. § 1.1);
C — тот же алгоритм, что и B, но дополненный процедурой сжатия компонент правой части ИСЛАУ на основе информации о монотонности оценки (см. § 1.1);
D — тот же алгоритм, что и C, но в котором рассекаемая компонента выбирается из условия максимизации произведения оценки производной на ширину интервала (как это описано в § 1.2).
Модификация Рона из § 1.3 нами в этих алгоритмах не использовалась. Фактически на основании приводимых ниже результатов можно сравнивать даже пять алгоритмов,
2Это фактически язык Fortran-90 со встроенными интервальными типами данных, операциями и отношениями между ними, см. http://www.sun.com.
включая самый простейший метод дробления параметров, поскольку он отличается от алгоритма А лишь длиной списка.
Таблица 8 4 х 4-система Ноймайера с диагональю Ь = 5.5
Алгоритм А В С Б
Количество итераций Максимальная длина списка Время счета 51 51 39 15 45 45 21 9 < 1 с < 1 с < 1 с < 1 с
5 х 5 система Ноймайера с диагональю Ь = 7
Алгоритм А В С Б
Количество итераций Максимальная длина списка Время счета 506 506 421 59 503 503 285 48 < 1 с 2 с 2 с < 1 с
6 х 6-система Ноймайера с диагональю Ь = 8.5
Алгоритм А В С Б
Количество итераций Максимальная длина списка Время счета 15760 15760 10655 441 15744 15744 4644 302 110 с 190 с 103 с 6 с
В качестве базового метода мы использовали интервальный метод Гаусса — Зейделя со стандартным предобусловливанием обратной средней матрицей, этот же метод использовался для внешнего оценивания "обратной интервальной матрицы" ИСЛАУ. Модельной задачей нам служила интервальная линейная система
( * [0,2] [0,2] *
V [0, 2] [0, 2]
[0, 2] \
[0, 2]
* )
х =
( [-1,1] \ [-1,1]
V [-1,1])
(20)
которая была предложена А. Ноймайером в книге [13]. Оказывается, что матрицы системы (20) четного порядка п невырождены при * > п, а матрицы нечетного порядка п невырождены при * > \/п2 — 1.
Некоторую информацию о строении объединенного множества решений 2 интервальной линейной системы (20) можно получить из соображений симметрии. Из-за уравновешенности вектора правой части рассматриваемая ИСЛАУ инвариантна относительно
изменения знаков всех компонент решения на противоположный, так что ее множество решений £ центрально симметричное относительно начала координат, в частности,
min{ x | x G £ } = — max{ x | x G £ }, i = 1, 2,...,n. (21)
Так как для любых i, j G {1, 2,... , n} после замены Xj на Xj и наоборот интервальная система (20) остается неизменной, множество £ симметрично относительно биссектрисы положительного и отрицательного ортантов пространства так что
min{ Xj | x G £ } = min{ Xj | x G £ }, max{ xj | x G £ } = max{ xj | x G £ }
для любых i,j G {1, 2, ...,n}. Сопоставляя эти соотношения с (21), можно заключить, что интервальная оболочка множества решений £ системы (20) является гиперкубом с центром в начале координат.
В случае приближения значения t к границам невырожденности (n для четных размерностей и Vn2 — 1 для нечетных) размеры объединенного множества решений ИСЛАУ (20) неограниченно возрастают. Варьируя t и n, легко получить из (20) набор тестов для проверки развитых нами алгоритмов оптимального решения "внешней задачи" для ИСЛАУ.
Результаты тестовых расчетов с интервальными линейными системами Ноймайера размера 4 х 4, 5х5 и 6х6 представлены в табл. 8. С одной стороны, подобные задачи являются все-таки не слишком тяжелыми и просчитываются рассматриваемыми нами алгоритмами "до конца", т. е. до того момента, когда ведущая оценка станет точной. А с другой стороны, комбинаторная сложность этих задач при решении "лобовым перебором" (количество крайних точечных систем уравнений, равное соответственно 220 , 230 и 242) достаточно велика и позволяет предметно судить о сравнительной эффективности вычислительных схем.
Наконец, для интервальной 7 х 7-системы Ноймайера с диагональю t =10 алгоритмы A, B и C вообще не просчитывают оптимальную оценку объединенного множества решений до конца даже за время порядка часов. Дело в том, что примерно после 50 000 итераций этих алгоритмов "быстрая" оперативная память использованной нами ЭВМ оказывалась исчерпанной списком L и начинавшийся обмен с "медленной" памятью на жестком магнитном диске практически сводил на нет производительность процессора, так что вычисления делались чрезвычайно медленными. Но алгоритм D все же позволяет успешно решить рассматриваемую задачу до конца за 112 с процессорного времени, сделав 5246 итераций-бисекций при максимальной длине списка 4050.
Из анализа табл. 8 можно сделать выводы об эффективности тех или иных модификаций методов дробления параметров. Довольно неожиданными являются данные о процессорном времени решения задач, согласно которым алгоритм B оказывается наименее эффективным. То есть затраты на исследование монотонности только по элементам матрицы оборачиваются (по крайней мере, для систем Ноймайера) лишь утяжелением алгоритма.
Интересно и поучительно сравнить по эффективности наши методы дробления параметров с методами Рона, предложенными в 80-е годы [17] и, как уже упоминалось, в настоящее время наиболее популярными среди зарубежных исследователей и пользователей интервального анализа [13, 22]. Методы Рона были реализованы и тщательно оттестированы К. Мадсеном и О. Тофтом [30, 31], которые использовали в качестве модельных
задач разнообразные интервальные линейные системы размерности от 2 до 30. Для нас особенную ценность имеют их результаты на ИСЛАУ предельной размерности 30 с тремя матрицами, которые являются интервализациями известных тестовых матриц вычислительной линейной алгебры из справочника [32]: матрицы
/10 0 0 1 0 0 0 1
V
000
1 2 3
0 0 0
1 2
3
1 п- 1
п — 1 п
(22)
/
матрицы конечно-разностной аппроксимации второй производном на симметричном шаблоне3
( 2 -1 -1 2 -1 -1 2 ••.
\
0
0
и матрицы
V
1 2
3
(23)
2 -1
-1 2 -1 -1 2
п - 1 п
п - 1 п
п — 1 п
п — 1 п — 1 п — 1 ••
п — 1 п
(24)
у п п п • • • п п у
Интервальный вектор правой части для всех трех случаев брался одинаковым и равным
([0.999,1.001], [0.999,1.001],..., [0.999,1.001] )Т.
Интересно, что еще в начале 90-х годов на персональном компьютере с процессором Intel 80486DX и математическим сопроцессором К. Мадсен и О. Тофт не смогли успешно просчитать до конца методом Рона интервальные линейные 30 х 30-системы:
1) с матрицей (22), все элементы которой равномерно уширены на интервал [-0.002, 0.002];
2) с матрицей (23), все элементы которой уширены на [-0.0003, 0.0003];
3) с матрицей (24), все элементы которой уширены на [-0.00001, 0.00001].
3Она является M-матрицей см., например, [13].
Вскоре задача 2 из этого списка все-таки была успешно решена, но для этого К. Мад-сену и О. Тофту потребовались привлечение транспьютера Meiko Transputer System и распараллеливание вычислений на 32 процессорах. Общее последовательное время работы всех процессоров составило при этом 2471.27 с (более 40 мин). О решении задач 1 и 3 даже на транспьютере в работах [30, 31] ничего не сообщается.
Что касается методов дробления параметров, то алгоритм, приведенный в табл. 8 под литерой D, успешно справился со всеми тремя упомянутыми задачами на однопроцессорной рабочей станции SUN Ultra-10. При этом:
Задача 1 была решена всего за две итерации-бисекции и время порядка одной секунды по каждой отдельной оцениваемой компоненте множества решений.
Задача 2 была решена всего за две итерации-бисекции и время порядка одной секунды по каждой из отдельных оцениваемых компонент множества решений.
Задача 3 даже с более широкой матрицей радиуса 0.0001 успешно решалась для отдельной оцениваемой компоненты множества решений за время, не превосходящее 14 мин и не более чем за 1646 итераций-бисекций. Этот худший результат был достигнут для 10-й компоненты множества решений, а для других компонент трудоемкость исполнения варьировалась в широких пределах — от 50-60 итераций-бисекций до нескольких сотен и тысячи с лишним.
Наконец, с помощью методов дробления параметров автором был получен рекордный результат. Алгоритм D позволил успешно вычислять оптимальные покоординатные оценки множества решений интервальной ИСЛАУ размера 300 х 300 с матрицей (22), уширенной на [-0.002,0.002] всего за две итерации-бисекции и время порядка 20 мин (по любой из компонент). Конечно, задача 1 с интервализацией матрицы (22) умеренно сложная, но тем не менее столь большие размеры решаемых систем в задаче оптимального внешнего оценивания множеств решений общих ИСЛАУ не рассматривал еще никто в мире!
В заключение можно резюмировать, что приведенные выше данные тестовых прогонов, а также свойство последовательной гарантии методов дробления параметров (см., например, [1, 33]) убедительно свидетельствуют в пользу их решительного преимущества перед методами Рона.
Список литературы
[1] Шарый С.П. Оптимальное внешнее оценивание множеств решений интервальных систем уравнений. Ч. 1 // Вычисл. технологии. 2000. Т. 7, №6. С. 90-113.
[2] Панков П. С. Алгоритм доказательного поиска экстремума с использованием миноранты по области // Изв. АН Киргизской ССР. 1979. №6. C. 12, 13.
[3] Панков П. С. Алгоритмы доказательства устойчивых утверждений и глобальной оптимизации в ограниченной области. Фрунзе, 1984. 13 с. Деп. в ВИНИТИ, №5250-84.
[4] Asaithambi N. S., Shen ZüHE, Moore R. E. On computing the range of values // Computing. 1982. Vol. 28, No 3. P. 225-237.
[5] Hansen E. Global Optimization Using Interval Analysis. N. Y.: Marcel Dekker, 1992.
[6] Kearfott R. B. Rigorous Global Search: Continuous Problems. Dordrecht: Kluwer, 1996.
[7] Ratschek H., ROKNE J. New Computer Methods for Global Optimization. Chichester, N. Y.: Ellis Horwood, Halsted Press, 1988.
[8] Алефельд Г., Херцбергер Ю. Введение в интервальные вычисления. М.: Мир, 1987.
[9] ЗОРИЧ В. А. Математический анализ. Т. 1. М.: Наука, 1981; T. 2. М.: Наука, 1984.
[10] Калмыков С. А., Шокин Ю. И., Юлдашев З. Х. Методы интервального анализа. Новосибирск: Наука, 1986.
[11] Shary S. P. On optimal solution of interval linear equations // SIAM J. Numer. Analysis. 1995. Vol. 32, No 2. P. 610-630.
[12] Gay D. M. Solving interval linear equations // SIAM J. Numer. Analysis. 1982. Vol. 19, No 4. P. 858-870.
[13] NEUMAIER A. Interval Methods for Systems of Equations. Cambridge: Cambridge Univ. Press, 1990.
[14] Ratz D. Automatische Ergebnisverifikation bei globalen Optimierungsproblemen. Ph.D. Dissertation. Karlsruhe: Univ. Karlsruhe, 1992.
[15] Ratz D., Csendes T. On the selection of subdivision directions in interval branch-and-bound methods for global optimization //J. Global Optimization. 1995. Vol. 7. P. 183207.
[16] Shary S. P. A new class of algorithms for optimal solution of interval linear systems // Interval Computations. 1992. No 2(4). P. 18-29.
[17] ROHN J. Systems of linear interval equations // Linear Algebra Appl. 1989. Vol. 126. P. 39-78.
[18] Shary S. P. Optimal solution of interval linear algebraic systems. I // Interval Computations. 1991. Vol. 1, No 2. P. 7-30.
[19] Oettli W. On the solution set of a linear system with inaccurate coefficients // SIAM J. Numer. Analysis. 1965. Vol. 2, No 1. P. 115-118.
[20] JANSSON C. Calculation of exact bounds for the solution sets of linear interval systems // Linear Algebra Appl. 1997. Vol. 251. P. 321-340.
[21] ШАРЫЙ С. П. Новый класс алгоритмов для оптимального решения интервальных линейных систем // Актуальные проблемы прикладной математики: Мат. конф., Саратов, 20-22 мая 1991. Саратов, 1991. С. 113-119.
[22] KOLEV L. V. Interval Methods for Circuit Analysis. Singapore: World Scientific, 1993.
[23] ЛАКЕЕВ А. В., Носков С. И. Описание множества решений линейного уравнения с интервально заданными оператором и правой частью // Докл. РАН. 1993. Т. 330, №4. С. 430-433.
[24] ЛАКЕЕВ А. В., Носков С. И. О множестве решений линейного уравнения с интер-вально заданными оператором и правой частью // Сиб. мат. журнал. 1994. Т. 35, №5. С. 1074-1084.
[25] KREINOVICH V., Lakeyev A. v., Noskov S. I. Optimal solution of interval linear systems is intractable (NP-hard) // Interval Computations. 1993. No 1. P. 6-14.
[26] KREINOVICH V., Lakeyev A. V, Noskov S. I. Approximate linear algebra is intractable // Linear Algebra Appl. 1996. Vol. 232. P. 45-54.
[27] KREINOVICH V., Lakeyev A., Rohn J., Kahl P. Computational Complexity and Feasibility of Data Processing and Interval Computations. Dordrecht: Kluwer Acad. Publ., 1997.
[28] ROHN J., KREINOVICH V. Computing exact componentwise bounds on solutions of linear system is NP-hard // SIAM J. Matrix Analysis Appl. 1995. Vol. 16. P. 415-420.
[29] ГЭРИ М., Джонсон Д. Вычислительные машины и труднорешаемые задачи. М.: Мир, 1982.
[30] MADSEN K., Toft O. A parallel method for linear interval equations // Interval Computations. 1994. No 3. P. 81-105.
[31] Toft O. Sequential and parallel solution of linear interval equations // Eksamensproject: NI-E-92-04, Numerisk Institute, Danmarks Tekniske Hojskole. Lyngby, 1992. 98 p.
[32] Gregory R.T., Karney D.L. A Collection of Matrices for Testing Computational Algorithms. N. Y.: Wiley Interscience, John Wiley, 1969.
[33] ШАРЫЙ С.П. Интервальные алгебраические задачи и их численное решение: Дис. ... д-ра физ.-мат. наук. Новосибирск, 2000.
[34] ПАПАДИМИТРИУ Х., СтАЙГЛИЦ К. Комбинаторная оптимизация. Алгоритмы и сложность. М.: Мир, 1985.
[35] ШАРЫЙ С. П. Алгебраический подход во "внешней задаче" для интервальных линейных систем // Вычисл. технологии. 1998. Т. 3, №2. С. 67-114.
[36] Heindl G., KREINOVICH V., Lakeyev A. Solving linear interval systems is NP-hard even if we exclude overflow and underflow // Reliable Computing. 1998. Vol. 4. P. 383388.
[37] MOORE R. E. Methods and Applications of Interval Analysis. Philadelphia: SIAM, 1979.
[38] Shary S. P. Algebraic approach in the "outer problem" for interval linear equations // Reliable Computing. 1997. Vol. 3, No 2. P. 103-135.
Поступила в редакцию 10 августа 2001 г.
