close

Вход

Забыли?

вход по аккаунту

?

Коларов В.Г. - Переодоление ошибок при вычислениях с плавающей запятой в линейном программировании (2006).pdf

код для вставкиСкачать
Переодоление ошибок при вычислениях с
плавающей запятой в линейном
программировании 1
(неотредактированный авторский перевод на русском языке,
версия 1 сентября 2006 г.)
Васил Георгиев Коларов2
1 сентября 2006 г.
1 обзор,
опубликованный на болгарском языке в журнале “Системи и управление”, изд. Высшего после-дипломного института при Университете национальной и
мировой экономики в городе Софии (бывший Высший экономический институт имени Карла Маркса), 1974, № 1 (стр. 41-51) и № 2 (стр. 33 - 41)
2 окончил Мехмат МГУ в 1966-ом году, кафедра Вычислительной математики
СОДЕРЖАНИЕ
1
Содержание
1 Решение задач ЛП
2
2 Представление чисел с плавающей запятой
3
3 Особености вычислений с плавающей запятой
3
4 Предварительное обсуждение проблемы ошибок
5
5 Задачи ЛП
5
6 Симплекс-метод
8
7 Эмпирическое исследование ошибок
12
8 Логическая часть симплекс-метода
14
9 Выбор метода
15
10 Симплекс-метод с LU−разложением
17
11 Толерансы
20
12 Реинверсии
22
13 Управление реинверсий. Проверки ошибок
27
14 Масштабирование
30
15 Заключение
33
16 Приложение
34
1 РЕШЕНИЕ ЗАДАЧ ЛП
1
2
Решение задач линейного программирования
Линейное программирование (ЛП) наиболее часто применяемый в практике
математический аппарат решения оптимизационных задач. Это обуславливается прежде всего следующими обстоятельствами. У линейного программирования своя, построенная в основном, теория (например [1], [6]). Разработаны эффективные методы решения задач. Современные компьютеры располагают программами, которые быстро решают задач линейного программирования с большими размерами. Например система математического программирования ОФФЕЛИЯ II, разработанная во Франции для компьютера
CDC 6600 решает примерно за 29 минут типичную задачу с 1 040 строками
ограничений, 8 056 переменными и 73 184 ненулевыми элементами в матрици
ограничений [12]. Есть сведения [16], что для компьютера ИЛЛИАК IV, при
котором экспериментируется идея о двухмерном массиве из 256 автономных
процессоров, разрабатывается (1969 год) система линейного программирования, которая благодаря паралельной работы процессоров решает за минут
две задачу с 4 000 строками, 10 000 переменными и 1% заполнености матрицы. Та же самая задача решается часов за семь мощным (замечание при
переводе: но однопроцессорным) компьютером IBM 360/75.
Наряду с этим важным фактором, помагающий широкому использованию
линейного программирования в человеческой практике решения оптимизационных задач — это возможность при помощи его моделировать множество
классов важных прикладных задач. Эта сторона линейного программирования обект многочисленных монографий, сборников, журналов, которые доступни в стране и может быть предметом самостоятельных обзоров и соответсвующих библиографий. Достаточно подчеркнуть, что ежедневно компьютеры в мире решают огромное число задач линейного программирования для
самых разных применениях. В экономике, военном деле и т.д. они обеспечивают возможность искать более эффективных вариантов, для раскрытия ряда значительных резервов. При этом существенная часть машинного времени
для подобных задач поглащается решением больших по размеру задач линейного программирования [44]. Обычайное явление задачи с 500 - 1 000 строками ограничений, но решаются и отдельные задачи с 1 000 - 10 000 строками
ограничений. В [18] Данциг сообщает, что при одном применении для “Нэшенел бискуит кампан” в США успешно решалась задача с около 32 000
ограничениями и почти с миллион переменными.
Объем вычислений при решении даже небольших задач линейного программирования однозначно определяет необходимость использования современных компьютеров.
2 ПРЕДСТАВЛЕНИЕ ЧИСЕЛ С ПЛАВАЮЩЕЙ ЗАПЯТОЙ
2
3
Представление чисел с плавающей запятой
При вычислениях такого рода обычно числа представляются в компьютере с
плавающей запятой, т.е. в виде
±m.β p ,
где m мантисса, β −1 6 m < 1, записанная в компьютере цифрами, число которых определено и конечно, а p - порядок числа, −n 6 p 6 N, N > 0, n > 0.
β - основание выбранной системы представления чисел (чаще всего β =
2, 16 или 10). Если
m = 0, m1 m2 m3 . . . mt
представление мантиссы в системе счисления с основой β , то цифры mi удовлетворяют условиям
0 6 mi < β
(i = 1, 2, . . .t).
Если еще m1 всегда отлично от нуля, говорят что представление чисел с плавающей запятой нормализовано.
Обычно в компьютере записываются последовательно знак мантисы, цифры mi (i = 1, 2, . . .t). и цифры "смещенного” порядка q = p + b , где например
b = n (т.е. знак порядка p не записывается отдельно, а помнится неявно).
Представление чисел с плавающей запятой привлекательно, так как позволяет работать с числами из широкого диапазона, при этом для их записывании всегда необходима порция памяти компьютера, которая имеет фиксированный объем, который не зависит от абсолютной величины числа в этом
диапазоне.
3
Особености вычислений с плавающей запятой
Вычисления с плавающей запятой обладают ряд особеностей, которых надо
учитывать при вычислениях с большим объемом. Прежде всего, дискретное
представление чисел означает, что в упомянутом диапазоне можно представит только конечное число чисел. Следовательно, если число, которое хотим
ввести в компьютере, не одно из конечного числа чисел, представимые точно
с плавающей запятой, несмотря на то, что оно из вопросного диапазона, то
в памяти компьютера записивается обычно одно из этих чисел, которое достаточно близко к оригинальному числу. Этим уже допускается определенная
неточность еще при вводе этого числа.
Из за фиксированной значности t мантисы числа с плавающей запятой,
произведение двух чисел, которое в общем случае нуждается в 2.t цифр, записывается или после отброса цифр, которые остаются после первых t значущих, или после соответствующего округления. В обоих случаях допускается
определенная погрешность по сравнению с точным результатом.
3 ОСОБЕНОСТИ ВЫЧИСЛЕНИЙ С ПЛАВАЮЩЕЙ ЗАПЯТОЙ
4
Вычитание двух близких чисел с плавающей запятой порождает другой
проблемы. Например, пусть нужно вычислить
z = x−y,
где x, y и z числа с плавающей запятой с восьмицифренной мантиси и с основании 10:
x = +0, 99999876.102
y = +0, 99999765.102 .
Тогда до нормализации
z = +0, 00000111.102 ,
а после нормализации
z = x − y = +0, 11100000.10−3 .
Беда в том, что у x и y по восемь значущих цифр, а у результата z всего три.
Остальные пять нулей генерированы компьютером при нормализации.
Если нужен результат с большей точности, необходимо [24, стр. 230231] или использовать плавающей арифметики с увеличенной точности (т.е.
с большим t), или переформулировать задачу в эквивалентной, но без вычитания близких чисел. К примеру, вместо y = 1 − cos x при маленких x можно
вычислять 2. sin2 ( 2x ) . При b2 ≫ a.c меньший по абсолютной стоимости корень
квадратного уравнения ax2 +bx+c = 0 лучше не вычислять по известной формуле
√
−b + b2 − 4ac
x1 =
,
2a
а по эквивалентной ей
2c
√
.
x1 =
b + b2 − 4ac
Из за этих особенностей арифметика с плавающей запятой не удовлетворяет некоторых основных законов - например ассоциативный [33, стр. 128 129] :
a = −0, 56324.105 ;
b = +0, 56323.105 ;
c = +0, 99216.100 ;
a + b = −0, 00001.105 = −0, 10000.101 ;
(a + b) + c = −0, 000784.101 = −0, 78400.10−2 ;
b + c = +0, 56323.105 ;
(b + c) + a = −0, 00001.105 = −0, 10000.101 ;
a + c = +0, 56323.105 ;
(a + c) + b = +0, 00000.105 = +0, 00000.100 .
(Числа записаны петцифренной мантиссой, а операции исполняются во вспомагательном регистре с емкостью десети цифр).
4 ПРЕДВАРИТЕЛЬНОЕ ОБСУЖДЕНИЕ ПРОБЛЕМЫ ОШИБОК 5
4
Предварительное обсуждение проблемы ошибок
Проблема ошибок вычислений с плавающей запятой при решении больших
по размеру задач линейного программирования не тривиальна [45, стр. 271].
При решении практических задач бывает, что соответствующие компьютерные программы не могут произвести разумные результаты от разумных входных данных [7, стр. 414]. При подобных ситуациях программы обычно не сообщают достаточно информации, наводящей на улучшение вещей.
Создается впечатление, что в специальной литературе эта проблема обсуждается мало и поверхностно. Исключение составляют быть может только публикации [45], [7], [34] и [16], которые обсуждаются подробно далее.
Вобщем можно сказать, что не создана целостная теория, описывающая в количественном и в качественном аспекте накопление ошибок в процессе этих
вычислений. Вместо этого молчаливо используются ряд недоказанных, более
или менее испитанных в практике средств для предотвращения, индикации и
уменьшения соответствующих ошибок [45], [34], [33] и [16]. Большая часть
из них описывается только в документации компьютерных программ и часто
представляют строжайший фирменный секрет.
В каждом конкретном пакете программ для оптимизации вопрос ошибок
вычислений решается определенным способом, учитивая научно-прикладных
исследований больших колективах специалистов (создание современного пакета программ линейной оптимизации оценивается на около 150 человеко-лет
высококвалифицированного труда). Вероятно одна из причин ограниченного
числа явных публикаций в этой актуальной области стремление сохранить от
конкурирующих фирм тайны любого исследования и любого решения вопросной проблемы.
Все это делает исследования в этой области актуальными и интересными
как с теоретической, так и с практической точки зрения.
5
Задачи линейного программирования
Общая задача линейного программирования, записанная в стандартной форме, может иметь следный вид: Требуется найти x , вектор-столбец с размерностью n , который удовлетворяет ограничениям
A.x = b,
(5.1)
где A матрица m × n, а вектор b размерностью m (m < n),
x > 0,
(5.2)
который минимизирует линейный функционал
L (x) = cT · x,
(5.3)
5 ЗАДАЧИ ЛП
6
где cT = (c1 , c2 , . . . , cn ) n−мерный вектор-строка, полученный транспонированием вектора-столбца c.
Здесь и далее в этом обзоре используются следующие обозначения. Если
A матрица, a ji — ее элемент в i−той строке и в j−том столбце; A j — j−той
столбец; A i — i−тая строка.
В основной литературе по линейному программированию описано как произвольную задачу линейного программирования можно привести в только что
указанной, скажем “стандартной” форме и рассмотрение ее не ограничивает
общность.
Каждый вектор, который удовлетворяет условий (5.1) и (5.2), называется допустимым вектором, а каждый допустимый вектор, минимизирующий
функционала (5.3) — оптимальным вектором или решением задачи (5.1),
(5.2) и (5.3).
Пусть A j1 , A j2 , . . . , A jm m линейно независимых столбцов матрицы A. Будем
говорить, что они формируют базис
B = (B1 , B2 , . . . , Bm ) = (A j1 , A j2 , . . . , A jm ),
а соответствующие x j1 , x j2 , . . . , x jm будем называть базисными переменными.
Не ограничивая общность будем предполагать, что столбцы A переупорядочены так, что j1 = 1, j2 = 2, . . . , jm = m, т.е. что
A=(B
..
. C ),
где C матрица с размерами m × (n − m) из остальных (замечание при переводе:
находящиеся вне базиса) столбцов A. Пусть yT = (x1 , x2 , . . . , xm ) m−мерный
вектор базисных переменных. Пусть еще
x T = ( yT
..
.
zT ),
где zT = (xm+1 , xm+2 , . . . , xn ) (n − m)-мерный вектор остальных, называемых
внебазисных (свободных) переменных, и
cT = ( d T
..
.
fT )
соответственно разложение c на m−мерный и (n − m)-мерный вектор, отвечающее на разложение x на y и z. Тогда выражения (5.1), (5.2) и (5.3) переписываются как следует:
Необходимо найти x
.
xT = ( yT .. zT ),
удовлетворяющий
B · y +C · z = b,
y>0
z > 0,
(5.4)
(5.5)
и для которого минимален функционал
L (x) = d T · y + f T · z .
(5.6)
5 ЗАДАЧИ ЛП
7
Так как по предположению B невырождена, то уравнение (5.4) может быть
решено относительно базисных переменных:
y = B−1 · b − B−1 ·C · z.
(5.7)
Функционал (5.6) так же может быть выражен посредством внебазисных переменных:
L (x) = d T · B−1 · b − (d T · B−1 ·C − f T ) · z.
(5.8)
Каждый вектор x, который удовлетворяет (5.4) при z = 0, называется базисным вектором, соответствующий базису B. Очевидно, что если x — базисный вектор, то
y = B−1 · b,
(5.9)
L (x) = d T · B−1 · b.
(5.10)
Из теории линейного программирования известно (например посмотрите
[31, стр. 25]), что если у задачи (5.1), (5.2) и (5.3) единственное решение x, то x
является базисным вектором при некотором базисе B. В случае, если решение
задачи — более чем один базисный вектор, то любая их выпуклая комбинация
тоже решение задачи (5.1), (5.2) и (5.3). Естествено, любое решение задачи
должно быть допустимым вектором.
Задачу нахождения допустимого базисного вектора можно свести (посмотрите например [1, стр. 100]) к решению вспомагательной задачи линейного программирования, у которой налицо очевидный допустимый вектор. В
результат ее решения или получается искомый допустимый вектор для (5.1) и
(5.2) , или устанавливается, что не существует такой вектор и следовательно
первоначальная задача не имеет решения. Ввиду этого достаточно рассматривать решение задачи (5.1), (5.2) и (5.3) при предположении, что существует
базис и соответствующий допустимый базисный вектор.
Интересно отметить, что обично программы современных компьютеров
для линейного программирования не используют этот способ нахождения
допустимого базисного вектора, а реализируют менее известного в нашей
стране так называемого композитного симплексного алгоритма (посмотрите [46], [33, стр. 47], [27, стр. 3], [23, стр. 6 - 8] и [28, стр. 2]), при котором
базисные векторы в отдельных итерациях не являются непременно допустимыми и минимизируется например общая сумма “неудовлетворения” отдельных ограничений.
Можно предположить также, что для допустимого базисного вектора сответствующий y > 0. В этом случае начальная задача (5.1), (5.2) и (5.3) называется невырожденной. Известен ряд способов переодоления возникшей
вырожденности при решении задач линейного программирования. (посмотрите [47], [1, стр. 231] и [33, стр. 131]), но следуя Уульфа [44] не будем рассматривать этот вопрос, так как он не имеет прямое отношение к проблеме об
ошибках.
6 СИМПЛЕКС-МЕТОД
6
8
Симплекс-метод
В разных своих модификациях симплекс-метод [1] — основной использованный метод решения задач линейного программирования. В [45] весьма метко
отмечено, что пока по меньшей мере “нет другой доступной альтернативы, которая была бы так эффективной”, как симплекс-метод.
После сделанных оговорок задача, подлежащая решению, формулируется
следующим образом:
Необходимо найти вектор x
.. T
. z ),
x T = ( yT
удовлетворяющий ограничений
B · y +C · z = b,
(6.1)
y>0
(6.2)
z > 0,
и который минимизирует функционал
L (x) = d T · y + f T · z .
(6.3)
при предположении, что B базис, а
T
T
x(0) = ( y(0)
..
.
T
z(0) ) ,
соответствующий базисный вектор, который предполагается допустимым, т.е.
(при невырожденной задаче)
следовательно
y(0) = B−1 · b > 0,
(6.4)
z(0) = 0,
(6.5)
L (x) = d T · B−1 · b.
(6.6)
Симплекс-метод состоит из последовательного выполнения ряда итераций. При каждой из них исходят из определенного базиса и соответствующего базисного вектора. Строится новый базис, различный от текущего в точно
одном столбце. Этому новому базису соответсвует новый базисный вектор с
меньшей стоимости функционала L (x). За конечное число итераций или находится оптимальный вектор, или устанавливается, что L (x) не ограничен снизу
на множестве допустимых векторов задачи.
Далее симплекс-метод рассмотрен более подробно, после чего кратко
формулируются основные правила алгоритма.
После решения уравнения (6.1) относительно базисных переменных и выражения формулы (6.3) через внебазисных переменных и после упрощающих
полаганий
b̄ = B−1 · b,
(6.7)
6 СИМПЛЕКС-МЕТОД
9
C̄ = B−1 ·C,
(6.8)
f¯T = d T · B−1 ·C − f T = d T · C̄ − f T ,
(6.9)
y = b̄ − C̄ · z,
(6.10)
L (x) = d T · b̄ − f¯T · z .
(6.11)
y(0) = b̄,
(6.12)
L (x(0) ) = d T · b̄ .
(6.13)
f¯T 6 0,
(6.14)
получается для каждого допустимого x
Более конкретно для x(0)
Если при текущем базисе B
то соответствующий допустимый базисный вектор x(0) — решение задачи, так
как у любого другого допустимого базисного вектора хотя бы одна координата вектора z была бы не ноль. Тогда ввиду условия (6.14) из формулы (6.11)
следует, что функционал имел бы большее значение.
В противном случае, т.е. если условие (6.14) не удовлетворено, существует
индекс s, для которого
f¯sT > 0.
(6.15)
Если соответствующая внебазисная переменная zs (до сих пор ноль в x(0) )
примет положительное значение, то для полученного вектора, назовем его
x(1) , будет верна формула
L (x(1) ) = d T · b̄ − f¯sT · zs
.
(6.16)
Но zs > 0 и f¯sT > 0, следовательно
L (x(1) ) = L (x(0) ) − f¯sT · zs < L (x(0) ) .
С точки зрения уменьшения L (x) вектор x(1) “лучше” по сравнению с x(0) .
При этом L (x(1) ) зависим линейно от zs и чем больше zs , тем меньше значение
L (x(1) ).
Конечно, по возможности увеличение zs должно произойти так, чтобы x(1)
был бы допустимым и базисным, как это было и при x(0) . Так как
z(1) = zs · es
(6.17)
(здесь es s−тый единичный вектор), то из выражения (6.10) получается для
x = x(1)
y(1) = b̄ − C̄ · z(1) =
= b̄ − C̄ · zs · es =
= b̄ − zs · C̄s
.
(6.18)
6 СИМПЛЕКС-МЕТОД
10
Если в выражении (6.18)
C̄s 6 0
(6.19)
,
то L (x) не ограничен снизу на множестве всех допустимых векторов, так как
построенный вектор x(1) представляет однопараметрическую фамилию допустимых векторов (параметр zs ), для которого L (x) уменьшается неограничено
с увеличением zs .
Допустимость x(1) следует из выражений (6.18) и zs > 0, а неограниченность L (x) — из (6.16), учитывая условий (6.15) и zs > 0.
В противоположном случае (6.19) существует хотя бы один индекс i, для
которого
c̄is > 0 .
(6.20)
Пусть
I{i : 1 6 i 6 m, i − целое, c̄si > 0} .
θ=
b̄r
b̄i
=
min
.
i∈I c̄is
c̄rs
(6.21)
(6.22)
zs не должно увеличиваться принимая значение больше чем θ , так как в
противном случае yr , т.е. r−тая координата (6.18), будет отрицательной ввиду
(6.22) и тогда x(1) не будет допустимым вектором. Полагая
zs = θ ,
(6.23)
получается, с одной стороны, наибольшее возможное уменьшение L (x), а с
другой — делает и
yr = 0.
(6.24)
Выражения (6.24) и (6.23) допускают естественое толкование — zs превращается в базисной переменной, а yr — во внебазисной переменной. Обновленный базис получается из столбцов старого базиса B заменой столбца
Br на Cs , т.е.
B(1) = (B1 , . . . , Br−1 ,Cs , Br+1 , . . . , Bm ).
(6.25)
Вместо того, чтобы заново решать (6.1) по отношение к новому базису,
ввиду “близости” B и B(1) достаточно в (5.7) сделать исключение по Гауссу
с ключевым элементом (pivot) c̄rs . В результат этого и после несущественной
размены местами переменных yr и zs , получается
−1
y(1) = (B(1) )
где
−1
· b − (B(1) )
·C(1) · z(1) ,
(6.26)
(y(1) )T = (y1 , . . . , yr−1 , zs , yr+1 , . . . , ym ),
(6.27)
(z(1) )T = (z1 , . . . , zs−1 , yr , zs+1 , . . . , yn−m ),
(6.28)
C(1) = (C1 , . . . , Cs−1 , Br ,Cs+1 , . . . , Cn−m ).
(6.29)
Одно доказательство утверждения, содержимое в последном предложении,
приведено в приложении в конце этого обзора.
6 СИМПЛЕКС-МЕТОД
11
Естествено, (6.26) можно записать более кратко снова как (6.10), но здесь
черточки сверху b и C уже означают (посмотрите так же (6.7), (6.8) и (6.9))
преобразование в отношение к новому базису B(1) .
Подобным образом, y(1) из (6.26) можно исключить из (6.3), получая таким путем
−1
L (x(1) ) = d T · (B(1) ) · b+
T
(1) −1
T
+ d · (B ) ·C − f · z(1)
(6.30)
при предположениях, что
d T = (d1 , . . . , dr−1 , fs , dr+1 , . . . , dm ),
(6.31)
f T = ( f1 , . . . , fs−1 , dr , fs+1 , . . . , fn−m ).
(6.32)
(6.30) можно получить и путем распространения операции исключения с
ключевым элементом c̄rs над уравнении (5.8).
Кратко, основные правила симплекс-алгоритма следующие:
I. Если при текущем базисе
f¯T = (d T · B−1 ·C − f T ) 6 0,
(6.33)
то соответствующий допустимый базисный вектор является оптимальным вектором.
II. Если существует s такое, что
f¯sT = (d T · B−1 ·C − f T )s > 0,
(6.34)
то s индекс столбца - кандидата для вхождения в базис.
а) если
C̄s = B−1 ·Cs 6 0 ,
(6.35)
то L (x) не ограничен снизу на множестве из допустимых векторов;
б) в противоположном случае, т.е. если не пусто множество
I{i : 1 6 i 6 m, i − целое, c̄si > 0} ,
то
(6.36)
i
b̄r
b̄i
(B−1 · b)
θ = r = min i = min
.
i∈I c̄s
c̄s
c̄is
(6.37)
Тогда r является индекс строки ключевого элемента c̄rs .
III. Выполняется исключение по Гауссу при помощи определенного ключового элемента. После размены местами между yr и zs , система ограничений
снова принимает вид (6.10) в отношение к новому базису. Подобным образом
трансформируется выражение для L (x). Новый допустимый базисный вектор получается из (6.10) и для него функционал принимает значение, которое
определяется из (6.10), полагая в обоих случаях z = 0. Переходится снова к
первому правилу для выполнения следующей итерации.
7 ЭМПИРИЧЕСКОЕ ИССЛЕДОВАНИЕ ОШИБОК
7
12
Эмпирическое исследование ошибок
При отсуствии соответствующих теоретических результатов относительно накоплении ошибок от вычислений с плавающей запятой в симплекс-методе,
представляют интерес эмпирические исследования в этой области. Студия
Мьюллер-Мербаха [32] едва ли не единственная публикация, в которой обстоятельно обсуждаются результаты представительного эмпиричекого изучения ошибок при решении задач линейного программирования. Исследования
проводились на 14 примерных задачах линейного программирования, размеры которых в интервалах 5 6 m 6 48, 9 6 n 6 75 и процент заполнености матриц ограничений p
q
p=
· 100 ,
m·n
где q — число ненулевых a ij , лежит в интервале 9, 9 6 p 6 84, 5. Большинство задач взяты из проекта “Стандартизованные численные эксперименты в
математическом программировании” [45].
Программы для компьютера Ай Би Эм 7090 (с фактической точностью
операций с плавающей запятой 26 битов, т.е. немного больше чем 8 десятичных цифр) реализируют трех модификаций симплекс-метода. Через каждые k
итераций (k = 1 для маленких и k = 5 для больших задач) вычисляются величины, характеризирующие накопленной ошибки. Более точно, пусть B матрица базиса, B−1 ее точная обратная матрица, а B−1
c — вычисленная обратная,
которая отличается от точной ошибками вычислений. Известно (посмотрите
например [42], [5]), что если
∆ = B · B−1
c −E .
(7.1)
и вычисление ∆ выполнить с повышенной точности (например используя арифметики с плавающей запятой с двойной точностью), то B−1 получается из приближенного равенства
−1
B−1 ≈ B−1
(7.2)
c − Bc · ∆ ,
при том умножение B−1
c · ∆ можно делать с обычной точностью. Следовательно матрица
(ε ji ) = B−1
c ·∆
с некоторым приближением представляет
мерой допущенных ошибок. Если
i
−1
i
еще обозначить (β j ) = Bc , то матрица e j , где
εi j
e ij = i
,
i
βj − εj (7.3)
можно интерпретировать, учитывая (7.2), как матрица относительных ошибок
при вычислении B−1 .
7 ЭМПИРИЧЕСКОЕ ИССЛЕДОВАНИЕ ОШИБОК
13
При этих обозначениях в [32] описаны результаты от вычислений следующих величин:
log e ij
,
(7.4)
M=∑
w
i, j
s
2
1
(7.5)
S=
· ∑ log e ij − M ,
w i, j
Emax = max log e ij .
(7.6)
i, j
В (7.4) и (7.5) w число тех элементов матрицы e ij , для которых
e ij > 0 .
(7.7)
Использование логаритма в приведенных формулах позволяет интерпретировать величин более непосредствено как порядок (позиция) первой ненулевой цифры в мантиссе соответствующего числа, представленное в десятичной
системе счисления.
Основные результаты можно попробовать формулировать кратко так:
• средний порядок относительной ошибки возрастает (не обязательно монотонно) начиная приблизительно с M = −8 до приблизительно M = −6
при 30-50 итерациях;
• стандартное отклонение обычно в порядке S = ±1;
• максимальный порядок ошибки варирует без видимой закономерности
обычно в интервале −6 6 Emax 6 −4, но в отдельных случаях достигает
и Emax = −3.
Наряду с этим в [32] обсуждаются и другие проблемы как например: испитание некоторых простых проверок индикации накопленных ошибок; проверка и елиминирование незначительных элементов; переодоление плохо обусловленных вершин многостена допустимых векторов; и другие.
В целом, можно оценить, что вопросное исследование, хотя и эмпирическое, обладает особой ценности, так как формирует количественное и качественное представление о процессе накопления ошибок. Результат Emax = −3
означает, что в конкретном случае при восемьцифренной мантиссе решение
было надеждно лишь до третьей цифры мантиссы. Это вывод еще раз подчеркивает важность проблемы анализа и управления ошибок, особенно учитивая
небольших размеров экспериментальных задач в [32].
8 ЛОГИЧЕСКАЯ ЧАСТЬ СИМПЛЕКС-МЕТОДА
8
14
Логическая часть симплекс-метода
Вероятно объязяны Филипу Уульфу за замечание [45, стр. 197] и [44, стр. 271],
которое очень полезно при анализе ошибок в симплекс-методе. Решение задачи линейного программирования при помощи симплекс-метода можно условно разделить на две части:
• логическая часть: в ней ищется и находится базис, которому соответствует оптимальный вектор;
• вычислительная часть: для найденного базиса вычисляются величины базисных переменных в соответствующем оптимальном векторе.
После того как в логической части правильно определен оптимальный базис, вычислительная часть должна решить обычной системы линейных уравнений. Эта задача хорошо изучена с точки зрения ошибок при вычислениях с плавающей запятой. Описаны соответствующие алгоритмы (например
разложение матрицы системы, A, при помощи исключений по Гауссу с полным или частичным выбором ключевого элемента, на произведение нижней
и верхной трехугольной матрицы A = L · U и последующее решение двух систем с трехугольными матрицами). Для этих алгоритмов выведены априорные
оценки относительной ошибки. Исследованы соответствующие итеративные
процедуры для возможного последующего уточнения найденного при помощи
точного метода решения. Существует ряд програм для современных компьютеров, реализирующие упомянотые алгоритмы. Посмотрите [42], [5], [3], [41],
[38], [37], [11], [43].
Важная отличительная характеристика больших по размеру задач линейного программирования ниская степень заполнености матрицы ограничений,
т.е. малый процент числа ненулевых элементов. Для больших задач эта величина как правило меньше 5%. Разработан ряд эффективных методов для решения систем линейных уравнений с “редкими” матрицами (и для обращения
“редких” матриц) — посмотрите прежде всего [40], [31, стр. 311-318], [30].
Основное применение этих методов в линейном программировании выполнение так называемых реинверсий — проблема, которая ввиду своего исключительного значения, обсуждается подробно далее. Конечно, в частности, вопросные алгоритмы позволяют успешно осуществлят упомянутой выше второй части симплекс-алгоритма, именно — при известном оптимальном базисе вычислить координаты соответствующего оптимального вектора.
Следовательно, можно сделать вывод, что для целей этого обсуждения,
подлинная часть симплекс-алгоритма его логическая часть, т.е. “открытие”
того базиса, которому соответствует оптимальный вектор. Накопление ошибок работы с плавающей запятой в часты вычислений может самое большое
изменить в определенных пределах некоторые значения чисел базисных переменных в оптимальном векторе. Ошибки при вычислениях в логической части
симплекс-метода могут по меньшей мере привести до неправильного объявления, что определенный вектор оптимален. Гораздо хуже случай, когда на-
9 ВЫБОР МЕТОДА
15
копленные ошибки в процессе решения могут привести к подмене первоначальной (исходной) задачи другой.
Вот почему далее рассматриваются основные способы при помощи которых предотвращается, выявляется, уменьшается или устраняется эффект
ошибок вычислений в логической части симплекс-алгоритма. Многие из способов, которые будут упомянуты, “изобретены” для практических целей, встроенны в современных компьютерных программ для линейного программирования, доказали в практике свою эффективность, хотя часто отсуствуют соответствующие обосновывающие теоретические разработки или скорее всего
публикации о них. Другие из способов имеют преимущественно теоретикоисследовательское значение и если будут упомянуты, это чтобы попробовать
создать более полную картину состояния вещей во вопросной области.
9
Выбор метода
В литературе по линейному программированию описано множество разных
модификаций симплекс-метода [6], [1]. Еще больше число модификаций, реализованное в многочисленных компьютерских программах для линейного
программирования.
Естествено и здесь выбор метода влияет чувствительно на накопление
ошибок при вычислениях с плавающей запятой, так же, как и при решении
других математических задач. Если судить по публикациям в специальной литературе, то прямые исследования теоретического или хотя бы эмпирического
характера этой интересной проблемы весьма скудны. Едва ли не единственное прямое исследование, которое можно найти, это работа Бартелса [7].
Опосредственно можно делать выводы на основе ряда сравнительных исследований разных модификаций симплекс-метода, в которых более конкретно оценивается число необходимых основных арифметических операций.
В этом отношение замечательными являются: монография Зойтендейка [2];
исследовение Уульфа и Катлера [45], в которой сообщаются и обсуждаются
результаты системного и очень широкого эмпиричного исследования; как и
публикация Смита и Орчад-Хейза [39].
Из упомянутых и других разработок (например [31]), ровно как и из анализа некоторых действующих компьютерских программ для линейного программирования, можно сделать вывод, что основной метод для решения типичных больших по размерам задач линейного программирования с “редкими” матрицами это так называемый мультипликативный модифицированный симплекс-метод (ММСМ) в разных своих вариациях. (Если матрица ограничений обладает специальной структурой, возможно преимущество будут иметь методы декомпозиции, при которых итеративно решаются
одна главная и множество подчиненных подзадач, все со существенно меньшими размерами по сравнению с первоначальной задачей, при чем осуществляется взаимный обмен информации между главной и подчиненними задачами).
9 ВЫБОР МЕТОДА
16
В ММСМ обратная матрица B−1 к матрице текущего базиса B представляется как произведение элементарных матриц, соответствующие последным
p операциям исключения
где
B−1 = E (p) · E (p−1) · . . . · E (1) ,
(9.1)
E k = (e1 , . . . , erk −1 , η k , erk +1 , . . . , em ) ,
(9.2)
ei i−тый m−мерный единичный вектор, а η (k) m−мерный вектор с координатами
(
c̄ i
i
i 6= r(k)
− c̄ sr
k
s
η
=
(9.3)
1
i = r(k) .
c̄ r
s
Благодаря специальной структуре матриц E k (только r(k) тый столбец различен от соответствующего столбца единичной матрицы), то достаточно в
памяти компьютеров хранить для данной E k только число r(k) , определяющее номер осбого столбца, и еще числа этого столбца, т.е. так называемый
η −вектор:
r(k) ; η1 , η2 , . . . , ηm .
Обычно даже η −вектор хранится плотно, т.е. только ненулевыми элементами, представленые парами: (ненулевой элемент и его номер)
r(k) ; (i1 , ηi1 ), (i2 , ηi2 ), . . . , (iq , ηiq ) .
Это целесообразно, так как η −векторы относительно разряжены (хотя достигается ценой специальных усилий).
При ММСМ правила симплекс-алгоритма из конца п.6 применяются учитывая (9.1). Более конкретно для правила I нужно вычислить координаты
(компоненты) вектора в левой части (6.33). Основная проблема здесь вычисление вектора так называемых симплекс-множителей π :
π T = d T · B−1
.
(9.4)
Учитывая (9.1) получается
π T = ( . . . ((d T · E (p) ) · E (p−1) ) . . . ) · E (1) .
(9.5)
Каждая из операций в скобках представляет более обобщенно взятое оперция типа vT · E (k) , для которой
где
vT · E (k) = (v1 , . . . , vr−1 , δ , vr+1 , . . . , vm ),
(9.6)
δ = v T · η = ∑ v i · ηi .
(9.7)
C̄s = E (p) · ( . . . (E (2) · (E (1) ·Cs ) ) . . . ),
(9.8)
Подобно, необходимый для правила II симплекс-алгоритма вектор из левой части (6.35) получается из
10 СИМПЛЕКС-МЕТОД С LU−РАЗЛОЖЕНИЕМ
17
где снова каждая из операций в скобках имеет специальный вид E (k) · v и выражается как
E (k) · v = (α1 , α2 , . . . , αm )T
(9.9)
при
i 6= r
v i + ηi v r
(9.10)
i=r.
ηr v r
Операция исключения в правило III ММСМ состоит просто из добавления текущего η −вектора, определенного на основе C̄s , к уже записанным в
памети компьютера η −векторов.
Хотя в литературе нет доступного прямого анализа проблемы ошибок в
ММСМ, косвенный анализ при помощью числа арифметических операций
и упомянутые уже эмпирические исследования [45] потверждают утверждения, что ММСМ более стабилен с точки зрения ошибок от работы с плавающей запятой по сравнению с другими версиями симплекс-метода. Учитывая
необходимость соответствующих теоретических разработок, предварительно можно высказать следующие соображения, подкрепляющие высказанной
выше тезы.
ММСМ (так же как и модифицированний симплекс метод с явной форме
обратной матрицы) [1] вычисляет на каждой итерации f¯T и C̄s используя первоначальных данных задачи, что вероятно уменьшает накопление ошибок от
вычислений с плавающей запятой. Далее, использование оригинальных данных, которые редки, и спецальные принимаемые меры, которые будут рассмотрени кратко дальше, приводят к тому, что и η −векторы обратной матрицы в степени и по мере возможности были бы тоже редкими. Тогда арифметическая операция (например в (9.7) и (9.10)) совершается только если оба
соответсвующие операнды ненулевые. Это тоже уменьшает число операций и
ошибки от работы с плавающей запятой.
αi =
10
Симплекс-метод с LU−разложением
Исследование [7], рядом с [8] и [9], представляет интересной попыткой анализировать проблемы ошибок при симплекс-методе и на этой основе предложить “стабилизированной” в отношении этих ошибок модификации алгоритма.
Объект внимания в начале [7] в основном операция исключения. Пусть
T
w = (w1 , w2 , . . . , wl ) l−мерный вектор, который преобразуется в vT = (v1 , v2 , . . . , vl )
при помощью операции исключения с ключевым элементом pβ — координата
с номером β вектора pT = (p1 , p2 , . . . , pl ). Используя так называемый обратный анализ ошибок, в [7] показано, что
v = f l(E β · w) = E β · (w + δ w),
(10.1)
где f l(α ) означение для результата вычислений выражения в арифметике с
плавающей запятой α , а E β элементарная матрица вида (9.2) и (9.3), соответствующая операции исключения по вектору p с ключевым элементом pβ .
10 СИМПЛЕКС-МЕТОД С LU−РАЗЛОЖЕНИЕМ
18
Другими словами, (10.1) означает, что вычисленный результат v точный результат исключения измененного вектора w + δ w. При этом
|δ wi | 6 |wi | · ε1 +
+
|pi |
· |wβ | · (2.ε1 + 3.ε12 + ε13 ).
|pβ |
(10.2)
Здесь ε1 оценка погрешности от совершения одной единственной элементарной арифметической операции с плавающей запятой.
Из (10.2) видно, что если |pβ | достаточно малое число по сравнению с
|pi | · |wβ |, то правая сторона (10.2) может быть произвольно большой. В [7]
приводятся простые примеры, которые показывают, что симплекс-метод, если в нем не наложены дополнительные ограничения снизу на абсолютную величину ключевого элемента, действительно приводит к операции исключения
со значительными ошибками при вычислениях.
В чем состоит предложение о стабилизировании симплекс-алгоритма?
При любой итерации в симплекс-методе необходимо решать системы линейных уравнений:
B·y = b,
(10.3)
π T · B = dT ,
(10.4)
B · C̄s = Cs .
(10.5)
B = L ·U.
(10.6)
На практике, так как в (10.3) b постоянно, система (10.3) решается на каждой итерации просто расспространением исключения по Гауссу и на b̄. Основная идея предложения в том, чтобы эти системы не решались с помощью
B−1 (сохраненная явно или как в ММСМ элементарными матрицами, которые умножаются), которая в конце каждой итерации обновляется исключением по Гауссу. Вместо этого три системы линейных уравнений решают разложением B на произведение нижней и верхней треухгольной матриц, так что
Следовательно, например B·v = q решается путем решения двух систем с треухгольными матрицами
L · t = q,
(10.7)
U · v = t.
(10.8)
Известно, что при определенных предположениях, которые в случае легко могут быть исполнены, решение системы линейных уравнений способом
LU−разложения и решение соответствующих треухгольных систем сопровождается априори ограниченными (при этом с приемливым пределом) ошибок от вычислений с плавающей запятой (посмотрите [42], [5], [3], [41]). Тогда
единственной проблемы остается нахождение простого и стабильного с точки зрения ошибок от арифметики с плавающей запятой споособа обновления LU−разложения в конце любой итерации, учитывая близкой связи между
двумя последовательными базисами в симплекс-методе.
10 СИМПЛЕКС-МЕТОД С LU−РАЗЛОЖЕНИЕМ
19
В [7] проблема решается, как следует. Пусть
B = (B1 , B2 , . . . , Bm ).
(10.9)
старый базис, для которого в силе (10.6), а
B(1) = (B1 , . . . , Br−1 , Br+1 , . . . , Bm ,Cs )
(10.10)
новый базис, полученный после выхода вектора Br и вхождения Cs на последнем месте в B(1) . Тогда
L−1 · B(1) =
(L−1 B1 , . . . , L−1 Br−1 , L−1 Br+1 , . . . , L−1 Bm , L−1Cs ) =
= (U1 , . . . ,Ur−1 ,Ur+1 , . . . ,Um , L−1Cs ) =
= H (1) .
(10.11)
Очевидно, что H (1) верхняя матрица Гессенберга (т.е. hij = 0 для i > j + 1),
в которой даже hij = 0 ниже главного диагонала в первых (r − 1) столбцов.
Она строится легко, так как вектор L−1 ·Cs уже вычислен при решении (10.7) с
q = Cs в правиле II симплекс-алгоритма из п. 6. Предлагается, при посредстве
ряда исключений по Гауссу, свести H (1) к верхней треухгольной матрице U (1) ,
т.е.
U (1) = E (m−1) · E (m−2) · . . . · E (r+1) · E (r) · H (1) ,
(10.12)
после чего новый базис снова разложен
B(1) = L ·U (1)
(10.13)
и можно перейти к следующей итерации симплекс-метода.
До сих пор с целью упрощения изложения не говорили о важной возможности: при начальном разложении в виде (10.6) и особенно при последовательности исключений (10.12) для сведения матрицы Гессенберга к верхней
треухгольной, можно менять местами строк с целью выбирать более хороший
ключевой элемент.
Далее в [7] сделан подробный анализ ошибок при предложенной модификации симплекс-метода. В результат этого определена одна, хотя и весьма
большая, но все таки априорная оценка того изменения первоначальной задачи, для которого “измененная” задача имеет как точное решение вычисленное
с плавающей запятой решение пероначальной (исходной) задачи.
В заключении [7] сообщаются результаты вычислительного эксперимента с реализацией модификации для компьютера Си Ди Си 6600, по сравнению со стандартным симплекс-алгоритмом. На указанном небольшом примере LU−модификация показывает большей точности, требуя около 1,6 раз
больше время для вычислений.
В [8] и [9] тот же автор вместье с Голубом предлагает процедуру на АЛГОЛе, которая реализирует идею, описанной в [7], с той разницей, что в конце
каждой итерации вместо более экономное приведение матрицы Гессенберга к
11 ТОЛЕРАНСЫ
20
верхней треухгольной, совершается совершенно новое LU−разложение нового базиса, конечно — только для его части вправо от (r − 1)-того столбца.
Строгий вывод (10.2) и проистекающие предупреждения относительно
абсолютной величины ключевого элемента при исключениях, в частности при
симплекс-методе, сделанные в 1971 году в [7], могут показаться и немного
опоздалымы, так как на практике в почти всех компьютерских программ для
линейного программирования, созданных после 1950 года, кандидат ключевой элемент ограничивается по абсолютной величине снизу нарочным прагом (толерансом), примерно равный 10−5 при ε1 = 10−15 (посмотрите (10.2)).
(Проблема толерансов обсуждается далее).
Что касается предлагаемого стабилизирования, которое впрочем тоже уязвимо выбором малого по абсолютной величине ключевого элемента, кажется
что главная проблема в том как использовать малой степени заполнености B
для поддержания редких L и U, что ускорило бы резко вычислений и между прочим уменьшило бы еще больше накапливаемые ошибки. Проблемой
представляет так же и неизменность и обновление L и U. В предлагаемой модификации L остается неизменной с начало и до конца. (В частности, если
начальный базис единичный базис, то L будет до конца единичной матрицой).
Интересно изучить есть ли другие LU−модификации, при которых обновляются и L, и U, но которые требуют меньше арифметических операций для одного обновления.
11
Толерансы
В логической части симплекс-метода ключевую роль играет ряд сравнений
вычисленных величин с нулем. Отнесение проверяемой величины к некоторому классу чисел, например положительные, неотрицательные и т.д., определяет в конечном счете какой быть замене в текущем базисе. В условиях
вычислений с плавающей запятой, при которых допускаются и накопляются
ошибки, вопросные сравнения с нулем могут привести к неправильному протеканию логической части алгоритма. Уульф отмечает в [45, стр. 281] “ефект,
который может не проявится пока не пройдет несколько итераций, будет сигнализироваться явлениями, против которых правила алгоритма должни были
бы предохранять как: ухудшение функционала; или появление отрицательных
чисел в текущем решении; или завершение с вектором, который не решение;
или с последующей невозможностью совершить реинверсии ввиду незамеченной вовремя линейной зависимости в данных задачи”.
Кажется единственный способ решить поставенной проблемы сравнивать
величины вместо с нулем — с малыми числами, так называемые толерансы
[44], [33, стр. 127], [32], [15], [31, стр. 318].
Например, проверка оптимальности (6.33) смягчается следующим образом:
f¯T = (d T · B−1 ·C − f T ) 6 ε1 > 0.
(11.1)
11 ТОЛЕРАНСЫ
21
Естествено, при этом существует опасность объявить какого-то вектора преждевременно оптимальным. Во всяком случае это намного лучше, чем продолжать вычислений с “надуманным” из-за ошибок f¯T .
Проверка (6.34) приобретает вид
f¯sT = d T · B−1 ·C − f T s > ε1 > 0.
(11.2)
Подобным образом (6.35) может перейти в
C̄s = B−1 ·Cs 6 ε2 > 0 .
(11.3)
Возможно L (x) будет объявлен неограниченным снизу, хотя некоторые компоненты C̄s положительны. Снова, это более приемливо, чем возможность
неправильно выбрать ключевой элемент, который малое положительное число только из за ошибок вычислений с плавающей запятой.
Сужение множества индексов (6.36) до множества
I{i : 1 6 i 6 m, i − целое, c̄si > ε2 }
(11.4)
может привести к появлению малых отрицательных значений некоторых базисных переменных, которые обычно замещаются автоматически нулями.
Так как выбор ключевого элемента для исключения по Гауссу имеет особое значение для точности вычислений, обычно накладываются более сильные требования к кандrдатам в ключевые элементы. Если можно ключевого
элемента искать в другие столбцы, то кандидат отбрасывается при
c̄sr < ε3 ,
(11.5)
где примерно ε3 = 10−5 . Наряду с этим при
c̄sr < ε4 · max |csi | ,
i
(11.6)
где, скажем ε4 = 10−2 , кандидат тоже отбрасывается, так как иначе породил
бы весьма большие по абсолютной величине элементы в соответном η −векторе.
Последняя проверка обосновывается например в [15] и [31].
Чтобы не порождались примерно в η −векторах ненулевые числа, которые
практически незначительны, обычно пользуются проверками
|η i | < ε5 ,
(11.7)
результат
< ε6 .
max |операнд|
(11.8)
где примерно ε5 = 10−12 и
(11.8) например сделает нулем результат вычитания двух весьма близких чисел, вместо того, чтобы допустить порождение относительно малого и “паразитического” (из за того, что не ноль) число.
12 РЕИНВЕРСИИ
22
Подобным образом, толерансы используются и при проверках как выполняются ограничения задачи.
Одна или другая система подобных толерансов предусмотрена в каждой
из современных компьютерских программных систем для решения задач линейного программирования. Конкретные величины толерансов, определение
которых считается искуством [44, стр. 283], зависит прежде всего от разрядности мантиссы чисел с плавающей запятой. Обычно тот, кто пользует соответствующих программ, может изменять стандартных толерансов, но
Орчард-Хейз в [33, стр. 128] отмечает, что “это все равно пробовать регулировать карбуратора современного автомобиля: лучше быть хорошим техником и хорошо знать двигателя”.
12
Реинверсии
Если для задачи линейного программирования, которой нужно решать, предварительно не известен начальный базис, вместо первоначальной задачи рассматривается другая, вспомагательная задача, в матрицы которой имеется
единичная под-матрица с размерами m × m. Естествено столбцы этой подматрицы формируют базис.
Пусть после p (p < m) итераций ММСМ
B(−1) = E (p) · E (p−1) · . . . · E (2) E (1)
(12.1)
и число базисных векторов, которые не единичные векторы, l. Очевидно l 6 p,
так как за одной итерации в базис может войти не более один существенный,
т.е. неединичный вектор. Обычно дааже число l чувствительно меньше p, что
объясняется возможностью некоторых существенных векторов входит в базис на некоторой итерации, выходить позднее из него, повторно входить и т.д..
Каждая итерация добавляет к списку η −векторов в памяти компьютера
новый η −вектор. Увеличение их числа замедляет все более основные операции симплекс-алгоритма, в которых участвует B(−1) . С нарастанием числа
η −векторов возрастает еще и их заполненность, т.е. число ненулевых компонент. Все это в конечном итоге увеличивает и ошибок от вычислений с плавающей запятой.
Реинверсии помагают решать эти проблемы. В случае, описанном выше,
возможно вычислить новое разложение B(−1) как произведение точно l элементарных матриц исключения:
B(−1) = Ê (l) · Ê (l−1) · . . . · Ê (2) · Ê (1) .
(12.2)
Другими словами, так как после p итераций в базисе только l (l < p) новых,
т.е. существенных векторов, то B(−1) может получиться более кратко после l
исключений (для единичных векторов исключения не нужны). Единственный
результат из предыдущих итераций, необходимый для реинверсии, последний
список номеров базисных векторов. В отличие от шагов симплекс-алгоритма,
12 РЕИНВЕРСИИ
23
при реинверсии налицо свобода выбора ключевого элемента. Следовательно,
реинверсия уменьшает ошибки при следующих вычислениях, с одной стороны, так как уменьшается число η −векторов, и с другой стороны, потому что
исключения при реинверсии совершаются с выбором ключевого элемента.
Оказывается, что число ненулевых компонент η −векторов зависит от порядка, в котором порождаются эти векторы. И здесь важен факт, что хотя при
определенном базисе его обратная матрица единствена, то ее разложение на
элементарные матрицы не единствено. В [33, стр. 74] приведен как иллюстрацию следующий пример. Пусть B имеет следующий


.
1

1

0

0
0
0
1
1
0
0
0
0
1
1
0
0
0
0
1
1
1

0

0

0
1
Если исключения совершаются сверху вниз по главному диагоналу, получаются следующие η −векторы:
η5
η
η
η
η
4
3
2
1
−1/2 0
0
0
1
1/2
0
0
1 −1
−1/2 1
1 −1 0 1/2
1 −1 0
0 1/2 −1 0
0
0
.
Место ключевого элемента показано жырным шрифтом.
Если порядок исключений изменится снизу вверх по главному диагоналу,
то η −векторы уже:
η5
η
η
η
η
4
3
2
1
−1/2 1 −1 1 −1
−1/2 1
0
0
0
1/2 −1 1
0
0 −1/2 1 −1 1
0 1/2 −1 1 −1 1 .
Второе представление обратной матрицы содержит уже 19 ненулевых элементов против 13 во первом.
12 РЕИНВЕРСИИ
24
Обратная матрица B−1 , записанная явно, заполнена целиком:


.
 1/2
1/2 −1/2 1/2 −1/2


−1/2 1/2
1/2 −1/2 1/2 



 1/2 −1/2 1/2
1/2
−1/2


−1/2 1/2 −1/2 1/2
1/2 
1/2 −1/2 1/2 −1/2 1/2
Более общим образом задача о реинверсии может быть поставлена так:
Пусть после определенного числа итераций текущий базис имеет вид
B = (B1 , B2 , . . . , Bm ).
(12.3)
Необходимо найти такое разложение B−1 на произведение элементарных матриц
B(−1) = E (q) · E (q−1) · . . . · E (2) E (1) ,
(12.4)
при котором общее число ненулевых элементов в матрицах E (k) (k = 1, . . . q)
минимално.
Тьюарсон пишет в [40], что не известно, чтобы поставленная задача решена в общем случае, хотя для многих и важных для практики специальных
случаев предложены точные решения и/или удовлетворительные эвристические решения.
Одна из первых простых и эффективных процедур для реинверсии в линейном программировании предложена Ларсеном [30]. Она устраняет приблизительно 40% из числа ненулевых элементов в η −векторах и из числа
необходимых арифметических операций для решения задач.
Современные компьютерские программы для линейного программирования используют более эффективные программы для реинверсии. Как видно
из [33, стр. 73] и [31, стр. 314], главные общие особености соответствующих
алгоритмов следующие.
12 РЕИНВЕРСИИ
25
Обычно матрица базиса B рассматривается как матрица специального
типа, как показано на следующей схеме.
Рис. 1: Структура матрицы базиса
В ней вне заштрихованной области элементы матрицы равны нулю, а
внутри заштрихованной области матрица редка, т.е. плотность ненулевых элементов мала. Использованы следующие обознечения:
• UT (Upper Triangle) Верхняя треухгольная часть;
• C (Core) Ядро;
• LT (Lower Triangle) Нижняя треухгольная часть.
Подматрица верхнего треухгольника BUT
UT нижняя треухгольная матрица,
следовательно последовательные исключения в столбцах BUT не влияют на
остальных столбцов матрицы. Эти исключения осуществляются в начале, при
чем необходимые η −векторы посылаются во внешней памяти компьютера.
Далее рассматриваются только остальные столбцы базиса, т.е. подматрица
C ∪ LT .
BC
∪ LT
Исключения в столбцах другой нижней трехугольной матрицы BLT
LT (нижC
ний треухгольник) не касают подматрицы B и так же могли бы быть совершены сразу, но это увеличить число ненулевых элементов в подматрице BCLT .
Это так, потому что строки из BLT с некоторой картиной расположения ненулевых элементов умножаются на определенные числа и добавляют к другим
12 РЕИНВЕРСИИ
26
строкам из BLT , которые в общем случае имеют другое расположение ненулевых элементов. Как увидится скоро, исключения в столбцах BC направлены
на преобразование подматрицы BC
C (ядро базисной матрицы) в треухгольную.
Вот почему, если временно отложить исключений в столбцов BLT , то после
превращении ядра в треухгольную подматрицу, то вся матрица, без уже обработанных векторов BUT , превращается в нижней треухгольной и соответствующие η −векторы могут быть произведены сразу, так как это было сделано в
начале для столбцов BUT . И так, исключения в столбцов BLT откладываются,
чтобы не породились бы лишные ненулевые элементы в BCLT и чтобы расшиC ∪ LT
рилась уже превращенная в трехгольную BC
C до нижней треухгольной BC ∪ LT ,
для которой непосредствено формируются соответствующие η −векторы.
Превращение ядра в треухгольной совершается так. Ключевой элемент b̄ij
выбирается таким образом, чтобы выполнились требования, подобные (11.5)
и (11.6), и еще, чтобы для него было минимально число
uij = (ri − 1)(c j − 1),
(12.5)
где ri число ненулевых элэментов в i−той строке, а c j — соответственно число ненулевых элементов в j−том столбце ядра. Этот выбор основывается на
том, что при исключении по j−том столбце с ключевым элементом b̄ij в самом плохом случае, т.е. при самом плохом полном несовпадении картин расположения ненулевых элементов, порождаются uij ненулевых элементов. На
самом деле, i−тая строка c j − 1 раз прибавляется, умноженная на константу к строке матрицы, а в ней ri − 1 ненулевых чисел, которые в общем случае
порождают ненулевые элементы. Генерированные ненулевые элементы могут
быть и меньше в своем числе, чем uij , так как на соответствующие места уже
могут быть ненулевые элементы.
На практике редко ищут абсолютный минимум uij в столбцах ядра — работают с минимумом в некотором подмножестве столбцов и строк из BC
C . После того как найден ключевой элемент, формируют необходимый η −вектор,
но только для триангуляции ядра (не для полного исключения). Все остальные столбцы подматрицы BC умножаются слева только что генерированной
элементарной матрицей. Строка и столбец выбранного ключевого элемента
отмечаятся как обработанными, снова ищут ключевой элемент в ядре, и так
далее.
После триангуляции ядра, последовательно производят η −векторы для
C ∪ LT .
исключений по Гауссу во всей треухгольной матрицы BC
∪ LT
Естествено, чтобы базисная матрица получила вид, показанный в схеме, не нужно делать фактическое физическое переупорядочивание ее строк и
столбцов. Соответствующие программы могут быть организованны так, чтобы обрабатывали матрицу, как будто она уже переупорядочена.
В заключении, реинверсии — основное профилактическое средство протев накопления ошибок от вычислений с плавающей запятой в логической
части симплекс-метода. Они уменьшают число η −векторов обратной матрицы и уменьшают число ненулевых элементов в них. В следующих итерациях
13 УПРАВЛЕНИЕ РЕИНВЕРСИЙ. ПРОВЕРКИ ОШИБОК
27
необходимо меньше арифметических операций и накопляется меньше ошибок от вычислений. Дополнительное уменьшение ошибки достигается за счет
возможности более свободно выбирать ключевой элемент.
13
Управление реинверсий. Проверки ошибок
Реинверсии отнимают определенное дополнительное время, которое при больших задачах может быть значительным. Вот почему они не могут осуществлятся весьма часто.
В некоторых современных программных системах для линейного программирования, например [23], предусмотрено, что тот, кто пользует программ,
определяет когда делать реинверсии. Классен в [15] считает целесообразным
делать реинверсии периодически, на каждые k итераций, где m/2 < k < m. В
[44] Уульф отмечает, что “частота реинверсий может быть задана предварительно, или ее осуществление может зависить от наблюдений хода решения
задач, которые совершаются программой или даже человеком”.
Смит и Орчард-Хейз в [39] описывают простой и изящный спооб как бы
“оптимального” управления периодических реинверсий. Цель — минимизировать общее среднее время для одной итерации, включительно время для
реинверсий. Если начать вести учет времени, прошедшее с начала последней реинверсии, то после каждой итерации можно вычислять общее среднее
время для совершения итераций до настоящего мамента. Если не делать новых реинвесий, то очевидно общее среднее время для итерации сначало будет сравнительно большым и будет уменьшаться в последующих итерациях
до достижения определенного минимума. Это обясняется тем, что на общее
среднее время в первых итерациях будет больше тяготеть время реинверсии.
Начиная с определенной итерации общее среднее время для одной итерации
начинает снова возрастать, так как обратная матрица базиса, представленная
η −векторами, становиться все больше и все более заполненой ненулевыми
элементами.
Тогда, предлагается следующий простой механизм управления реинверсий. Когда общее среднее время для итерации увеличиться более чем скажем
на 1% и при этом число итераций после последной реинверсии относительно
близко к числу итераций до вопросной реинверсии (с толерансом, примерно, 25%), то принимается, что наступило время для новой реинверсии. При
такой политики честоты реинверсий оказывается, что среднее время для реинверсий занимает около 15% общего времени для решения задач, включая
времени на реинверсий. Есть основания предполагать, что подобное управление реинверсий реализовано кроме в системе LP/90 для Ай Би Эм 7090,
как описано в [39], еще и в ряд других современных программных систем для
линейного программирования.
Описанное касается периодических реинверсий, которые являются профилактическим средством против накоплении ошибок вычислений с плавающей запятой. Наряду с этим, однако, реинверсии используются и как “те-
13 УПРАВЛЕНИЕ РЕИНВЕРСИЙ. ПРОВЕРКИ ОШИБОК
28
рапевтическое” средство, когда возникла ситуация (быть может много до
наступления время для профилактической реинверсии), дающая оснований
предполагать, что накопленные ошибки уже весьма большие. Сериозна проблема: каким быть проверкам, которые сигнализируют о попадении в подобной ситуации. С одной стороны, они должны иметь сильной коррелации с накопленной ошибкой, с другой стороны — должны быть простыми, чтобы не
отнимать весьма много времени.
При выполнении итераций симплекс-метода, если были другие возможности, то ключевой элемент отбрасывался при условий типа (11.5) или (11.6).
Если исчерпаны все возможности, т.е. вообще нет настолько хороший кандидат для ключевого элемента, то обычно требования смягчаются. Если и после
этого нет подходящего ключевого элемента, то это принимается как сигнал,
что реинверсия нужна.
В [32] изучены эмпирически пять разных проверок возможного сигнализирования накопленных ошибок. Их коррелация со средней ошибкой оказывается слабой, из за чего они не рассматриваются в этом обзоре.
Орчард-Хейз описывает в [33] следующие две естественные проверки, которые впрочем реализованны по меньшей мере как финальные проверки в системе [23]:
• вычисляется вектор:
T
(B · y − b̄) =
∑
k
bik · yk − b̄i
!
,
(13.1)
i=1,2, ... ,m
который в идеальном случае должен быть нулевым, т.е. проверяется насколько текущий базисный вектор допустим;
• вычисляется
π · B − dT =
∑ πk · bkj − d j
k
!
,
(13.2)
j=1,2, ... ,m
что тоже должно быть нулевым вектором.
Эти проверки на самом деле просты, так как требуют одно чтение из внешней памяти компьютера базисных строк (столбцов) основной матрицы и соответственно их умножение на векторы, которые постоянно находятся в оперативной памяти компьютера. Проверки (13.2) не нужно программировать
отдельно, так как соответствующие вычисления (так называемая оценка)
столбцов вне базиса, является частью симплекс-алгоритма.
Более сложны, но за то и более индикативны, другие две проверки, реализированые по меньшей мере в системах линейного программирования [27],
[28] и [29]. Основные операции в каждой итерации ММСМ — вычисление
симплекс-множителей
π T = d T · B−1
(13.3)
13 УПРАВЛЕНИЕ РЕИНВЕРСИЙ. ПРОВЕРКИ ОШИБОК
29
и преобразованного столбца, входящего в базис:
C̄s = B−1 ·Cs ,
(13.4)
При вычислении π T (посмотрите
(9.5)) нужны в обратном порядке η −векторы,
−1
представляющие B = β ji . Если d T состоял бы из одных единиц, то очевидно будет верно тождество
d T · B−1 = (1, 1, . . . , 1) · B−1 =
!
∑ β1i , ∑ β2i , . . . , ∑ βmi
=
i
i
(13.5)
.
i
После умножения справа на B получается:
∑ β1i , ∑ β2i , . . . , ∑ βmi
(1, 1, . . . , 1) =
i
i
i
!
·B.
(13.6)
Первая проверка состоит в умножении слева обратной матрицы базиса на
вектор из одных единиц, после чего полученным вектором умножают слева
матрицы базиса:
(ε1 , ε2 , . . . , εm ) = ((1, 1, . . . , 1) · B−1 ) · B .
(13.7)
Согласно (13.5) и (13.6) результат должен быть снова вектор из одных единиц. Максимальное по абсолютной величине отклонение
εB = max |ε j − 1|
(13.8)
j
является одной мерой ошибок при умножении на B−1 , когда η −векторы считываются в обратном порядке.
Вторая проверка опирается на тождестве
!T
B−1 ·
∑ b1j , ∑ b2j , . . . , ∑ bmj
j
= (1, 1, . . . , 1)T ,
(13.9)
j
j
которое проверяется сразу умножением слева на B.
Вычисляются суммы строк базисной матрицы, после чего обратная матрица умножается справа на полученный вектор. При этом (посмотрите (9.8))
обратная матрица, представленная η −векторами, используется в прямом
порядке. Если результат умножения будет
(ε1′ , ε2′ , . . . , εm′ )T =
= B−1 ·
∑ b1j , ∑ b2j , . . . , ∑ bmj
j
j
j
!
T
,
(13.10)
14 МАСШТАБИРОВАНИЕ
то
30
εF = max |ε ′j − 1|
(13.11)
j
мера ошибок при умножении на B−1 , когда η −векторы считываются в прямом
порядке.
Обычно
ε = max(εB , εF )
называется максимальной ошибки вычислений (maximum processing error),
которая сравнивается с определенным толерансом, например 10−2 при 31−
битовой мантиссе [28] и [29]. Если ошибка превзойдет этот толлеранс, то начинается черезвичайная реинверсия. Если и после ней ошибка больше толеранса, то решение задачи прекращается.
14
Масштабирование
Прежде чем обсуждать масштабирование в задачах линейного программирования, полезно обрисовать кратко состояние аналогичной проблемы при
решение сродных задач, например некоторые из задач линейной алгебры.
Известно (посмотрите например [3], [4], [5], [20]), что число
cond(A) = kAk · kA−1 k ,
(14.1)
называемое обусловленость матрицы A, где k.k некоторая из матричных
норм, например спектральная, играет важную роль для оценивания ошибок
при решении системы
A·x = b
,
(14.2)
как и при задач нахождения A−1 . Более конкретно, если A известна точно, b
содержит некоторые ошибки δ b по сравнению с b из (14.2) и x + δ x решение
A · (x + δ x) = b + δ b
то
(14.3)
,
kδ xk
kδ bk
6 cond(A)
kxk
kbk
,
(A + δ A) · (x + δ x) = b
,
(14.4)
где kxk эвклидова норма вектора x. Если A изменена по сравнению с (14.2) на
δ A, b точный, а x + δ x соответствующее решение
(14.5)
то по подобию с (14.4) верно
kδ Ak
kδ xk
6 cond(A)
kx + δ xk
kAk
.
(14.6)
14 МАСШТАБИРОВАНИЕ
31
Ввиду этого, естествено до того как решать (14.2) поискать такое масштабирование A, т.е. такие невырожденные диагональные матрицы D1 и D2 ,
чтобы обусловленость масштабированной A, т.е. D1 AD2
cond(D1 AD2 )
(14.7)
.
была бы минимальна.
Эта проблема изучается Бауером [10] , но ее решение зависить от матриц
|A| · |A−1 | и |A−1 | · |A|. Здесь |A| матрица абсолютных величин элементов матрицы A, т.е |A| = (|aij |). Как справедливо отмечает в [20] Форсайт, это уменьшает практической стоимости соответствующих результатов, так как первоначальная (исходная) задача (14.2), которую нужно решать, проще задачи нахождения A−1 , с помощью которой могли бы искать оптимального масштабирования.
Другой результат Бауера имеет фундаментальное значение для оценки
эффекта масштабирования при решении систем линейных уравнений при помощью исключений по Гауссу [20, стр. 505]. Оказывается, что если упорядоченное множество ключевых элементов фиксированно заранее, любое масштабирование A при помощью множителей, которые являются степени основы представления чисел с плавающей запятой, не меняет ни одной значущей цифры междинных и окончательных результатов исключений по Гауссу.
(Масштабирование степенями основы системы с плавающей запятой предпочитаются, так как при нем не вносятся ошибки от вычислений с плавающей
запятой).
Этот результат можно толковать так. Масштабирование может влиять
только на выбор и порядок ключевых элементов.
Примеры показывают, что ряд простых видов масштабирования, которые
уменьшают обусловленности одных матриц, могут увеличивать обусловленность других.
Интересны исследования ван дер Слуиса в [37] и [38]. Они рассматрывают
эффект балансированного масштабирования строк (столбцов) матрицы A в
задаче (14.2) при векторной норме
k x k∞ = max | xi |
i
.
(14.8)
Точнее, матрица называется балансированно масштабированной по строкам, если ее строки с равными ∞− нормами, т.е. если
k Ai k∞ = const
(i = 1, 2, . . . , m).
(14.9)
Аналогично, говорят, что матрица балансированно масштабированна по
столбцам, когда
k A j k∞ = const ( j = 1, 2, . . . , m).
(14.10)
Основной результат в [38] определение как “наилучшего” в определенном смысле то масштабирование, при котором сначала матрица масштабиру-
14 МАСШТАБИРОВАНИЕ
32
ется балансированно по столбцам, после чего полученная матрица масштабируетса балансированно по строкам, и в конце совершается преобразование обратное балансированому масштабированию по столбцам, т.е. матрица
“расмаштабируется” по столбацам.
Что касается задачи нахождения собственных чисел и векторов матрицы A , то Озборн обосновывает и предлагает [34] подходящее масштабирование, которое с определенными модификациями реализованно в процедурах на
АЛГОЛе [35]. Масштабированная матрица — с минимальной нормой и еще
балансированна в смысле:
(D1 AD2 )i = k(D1 AD2 )i k (i = 1, 2, . . . , m),
(14.11)
т.е. i−тая строка и i−тый столбец с одинаковыми нормами.
Масштабирование задач линейного программирования изучено значительно более слабо. В [44] Уульф пишет “ нам не известнd количественная оценка
важности масштабирования; единственное то, что использованные программы часто решают подходящим образом масштабированных задач, которые не
могли быть решены до этого”.
В [33] приводится следный практический критерий для масштабирования
матрицы ограничений
max
|a ij |
min
|a ij |
i, j
i, j
6R
(a ij 6= 0)
,
(14.12)
где R примерно 106 .
Простой алгоритм для нахождения множителей по столбцам и по строкам,
при которых левая сторона (14.11) минимизируется, описан Фулкерсоном и
Уульфом в [22].
В [33] Орчард-Хейз упоминает другой критерий при определении масштабирования матрицы ограничений, который кажется весьма подходящим.
По возможности все ненулевые элементы матрицы следует “центрировать”
к единице. Это можно например осуществит подбирая такие множители по
столбцам S j и строкам Ri , для которых функционал
∑i
(i, j :a j 6=0)
1 − S j · aij · Ri
2
,
(14.13)
наименьший.
Стоит отметить, что масштабирование при задачах линейного программирования кажется хотя бы настолько необходимым, как и при решении систем
линейных уравнений, потому что кроме всего прочего масштабирование делает чисел задачи более однородными и следовательно более равноправными
при множестве сопоставлений с фиксированными толерансами.
15 ЗАКЛЮЧЕНИЕ
15
33
Заключение
Состояние вещей в рассмотренной области делает актуальным как с теоретической, так и с практической точке зрения решение ряда проблем, например:
• анализирование ошибок от вычислений с плавающей запятой в ММСМ;
• исследование эффекта масштабирования матриц задач линейного программирования;
• исследование возможностей ускорения LU−модификации симплекс-метода.
Первую из указанных проблем, к примеру, можно рассмотреть с разных
сторон. Прежде всего надо решить какой тип анализа ошибок необходим.
Априорный анализ оценил бы скажем относительную ошибку независимо
от промежуточных и конечных результатов при решение конкретной задачи. Этот вид анализа ценен при выборе соответствующего метода. Уульф и
Катлер в [45] отмечают, однако, что “так как поведение (разных вариаций
симплекс-метода) зависит сильно от особеностей процесса, которые не могут быть известными заранее, априорные оценки их эффективности внушают
немного доверия”. Есть осонования предполагать, что подобное утверждение
верно и для оценок ошибок.
В определенном смысле полную противоположность априорному анализу предлагает динамически апостериорный анализ, идея которого оценивать
ошибки в процессе решения задачи путем оценки ошибок, возникающие при
выполнении каждой элементарной арифметической операциив компьютере.
В специальной литературе это направление представлено в общем плане или
по отношению других математических задач [25], [20], [14] и [21]. Может быть
понастоящем наиболее сериозный недостаток этого типа анализа, при отсуствии соответствующей специализированной аппаратуры (хардуер), несколькократное замедление решении задачи, которое при и без того значительном
объеме вычислений в линейном программировании, было бы немыслимым
люксом даже и для небольших задач.
Компромиссом между упомянутых двух полюсов быть может предлагает
подход, примененный к решении систем линейных уравнений в [13]. Выводятся апостериорные оценки ошибок, которые зависят от некоторых данных,
полученные в процессе решения задачи. Оценки могут вычислятся компьютером и выводиться вместе с решением. С одной стороны, они не такие общие,
как обычно бывают априорные оценки, и в этом смысле — полезны. С другой стороны, они не результат возможно самого подробного прослеживания
накопленных ошибок и благодаря этому вычисляются почти “бесплатно”.
Независимо от того какое решение будет принято относительно типа анализа, вероятнее всего метод анализа будет так называемый обратный анализ ошибок, при котором ошибки вычисления “перебрасываются назад” как
изменения в первоначальных данных, после чего ищут оценки этих изменений. Этот метод, введенным Гивенсом (посмотрите насчет этого [20, стр. 506])
16 ПРИЛОЖЕНИЕ
34
и разработанный в основном Уилкинсоном (посмотрите например [42], [43]),
применен успешно при анализе ошибок решая основных задач линейной алгебры.
Целесообразным кажется как объект анализирования ошибок принять
прежде всего самую развитую модификацию симплекс-метода, ММСМ, которая и чаще всего используется в практике. Желательно исследовать ММСМ в его самом современном виде, включающий к примеру эффективные алгоритмы реинверсий; учитывающий всюду, где это необходимо, малой плотности ненулевых данных, и так далее.
Если окажется возможным определить вычислимых апостериорных оценок ошибок в ММСМ, подобные оценкам при системах линеных уравнений
[13], они могли бы вычислятся периодически при решении задачи компьютером и при возможности служить более верному индицированию накопленных ошибок. Подобные оценки могли бы далее сравниваться с определенными пределами изменения начальных данных оптимизационной задачи, в рамках которых допустимый базисный вектор продолжает оставаться допустимым базисным вектором. Пределы этого типа определяются к примеру в [36].
Для проведения этого исследования и оформления обзора автор благодарен ряду коллег, которые создали благоприятных условий, прочли внимательно рукопись и помогли своими материалами, обсуждениями и критическими
замечаниями, особенно члену корреспонденту БАН проф. д-ру Благовесту
Сендову, доц. к.э.н. К. Габровскому, н.с. В. Спиридонову (рецензент), доц. Д.
Токареву (рецензент), н.с. Л. Филипову, н.с. Р. Карадимовой, н.с. К. Мариновой, н.с. Мл. Карадимову.
16
Приложение
Требуется доказать утверждение из п. 6, что исключение в выражении (5.7) с
ключевым элементом c̄rs приводит к (6.26).
Для этого используется факт, что операция исключения с ключевым элементом c̄rs эквивалентна умножению слева на специальную элементарную матрицу E (k) , определенной в (9.2) и (9.3). Легко проверяется равенство
E (k) = E −
1
(C̄s − er ) · eTr
c̄rs
,
(16.1)
где E единичная матрица.
Умножается (5.7) слева на E (1) и требуется доказать, что из него следует
(6.26).
E (1) · y = E (1) · B−1 · b − E (1) · B−1 ·C · z .
Представляется z как сумма z − zs · es и zs · es , т.е. “почти” вектора z, но с
нулевой s−той координатой, и “почти” нулевого вектора, но с s−той координатой zs :
E (1) · y = E (1) · B−1 · b−
16 ПРИЛОЖЕНИЕ
35
−E (1) · B−1 ·C · (z − zs · es ) − E (1) · B−1 ·C · zs · es
.
(16.2)
Последный член преобразуется
E (1) · B−1 ·C · zs · es = zs · E (1) · B−1 ·C · es =
= zs · E (1) · B−1 ·Cs = zs · E (1) · C̄s = zs · er
,
(16.3)
т.е. после исключения с ключевым элементом c̄rs из C̄s , этот столбец превращается в единичном. После подставления (16.3) в (16.2) получается
E (1) · y = E (1) · B−1 · b−
−E (1) · B−1 ·C · (z − zs · es ) − zs · er
.
(16.4)
Чтобы показать, что E (1) · B−1 обратная матрица новому базису B(1) , достаточно вычислить
E (1) · B−1 · B(1) =
= E (1) · B−1 · (B1 , . . . , Br−1 ,Cs , Br+1 , . . . , Bm ) =
= E (1) · (e1 , . . . , er−1 , B−1 ·Cs , er+1 , . . . , em ) =
= (e1 , . . . , er−1 , E (1) · B−1 ·Cs , er+1 , . . . , em ) =
= (e1 , . . . , er−1 , er , er+1 , . . . , em ) = E .
−1
= E (1) · B−1 так:
Тогда (16.4) переписывается с использованием B(1)
−1
· b−
E (1) · y · B(1)
−1
·C · (z − zs · es ) − zs · er .
− B(1)
Левая сторона (16.5) перерабатывается
T
1
(1)
E · y = E − r C̄s − er er · y =
c̄s
1
C̄s − er · yr =
r
c̄s
yr
= (y − yr · er ) + yr · er − r C̄s − er =
c̄s
C̄s − er
= (y − yr · er ) + yr · er −
=
c̄rs
= y−
= (y − yr · er ) + yr · η =
= (y − yr · er ) + yr · E (1) · er =
= (y − yr · er ) + yr · E (1) · B−1 · Br =
(16.5)
СПИСОК ЛИТЕРАТУРЫ
36
−1
(1)
· B r · yr .
= (y − yr · er ) + B
(16.6)
Подставляя (16.6) вместо левой стороны (16.5) получается
(y − yr · er ) + B
(1)
−1
· B r · yr =
−1
−1
·C · (z − zs · es ) − zs · er
· b − B(1)
= B(1)
или окончательно получается
−1
· b−
(y − yr · er ) + zs · er = B(1)
−1
−1
· B r · yr
·C · (z − zs · es ) − B(1)
− B(1)
.
(16.7)
Если теперь использовать (6.27), (6.28) и (6.29), то (16.7) превращается в
желанное (6.26).
Список литературы
[1] G. B. Dantzig, Linear Programming and Extentions, Princeton University
Press, Princeton, New Jersey, 1963
(та же книга: на русском языке, Издательство “Прогресс”, Москва,
1966)
[2] G. Zoutendijk, Methods of feasible directions, Elsevier Publishing
Company, New York, 1960
(та же книга: на русском языке,“Иностранная литература”, Москва,
1963)
[3] A. Ralston, A First Course in Numerical Analysis, McGrow-Hill Book
Company, New York, 1965
(эта же книга: на болгарском языке, “Наука и изкуство”, София, 1972)
[4] Д. К. Фаддеев, В. Н. Фаддеева, Вычислительные методы линейной алгебры, Госсиздат Физико-математической литературы, МоскваЛенинград, 1963
[5] G. E. Forsythe, C. B. Moler, Computer Solution of Linear Algebraic
Systems, Prentice-Hall, Inc., Englewood Cliffs, N. J., 1967
(та же книга: на русском языке,“Мир”, Москва, 1969)
[6] Д. Б. Юдин, Е. Г. Гольштейн, Линейное программирование -теория, методы, применения, “Наука”, Москва, 1969
СПИСОК ЛИТЕРАТУРЫ
37
[7] R. H. Bartels, A stabilisation of the simplex method, Numerishe
Mathematik, 16, 5, 1971, 414-434
[8] R. H. Bartels, G. H. Golub, The simplex method of linear programming
using LU decomposition, Communications of the ACM, 12, 5(May 1969),
266-268
[9] R. H. Bartels, G. H. Golub, Algorithm 350: Simplex Method procedure
employing LU decomposition, Communications of the ACM, 12, 5(May
1969), 275-278
[10] F. L. Bauer, Optimaly scaled matrices, Numerishe Mathematik, 5, 1, 1963,
73-87
[11] H. J. Bowdler, J. H. Wilkinson, . . . Solution of real and complex systems of
linear equations, Numerishe Mathematik, 8, 5, 1966, 217-234
[12] CDC 6600 OPHELIE II Mathematical Programming System
[13] B. A. Chartres, J. C. Geudor, Computable error bounds for direct solution
of linear equations, Journal of the ACM, 14, 1(Jan. 1967), 66-71
[14] B. A. Chartres, Automatic controled precision calculations, Journal of the
ACM, 13, 1966, 66-71, 386-403
[15] R. J. Clasen, Techniques for automatic tolerance control in linear
programming, Communications of the ACM, 9, 11(Nov.
1966), 802-803
[16] Review 18,939; Computing Reviews 11, 4(April 1970)
[17] G. B. Dantzig, Linear Programming and Extentions, Princeton University
Press, Princeton, New Jersey, 1963
[18] G. B. Dantzig, Linear programming and its progeny, Naval Research
Reviews, June 1966
[19] R. H. Dargel, F. R. Loscalzo, T. H. Witt, Automatic error bounds on real
zeros of rational functions, Communications of the ACM, 9, 11 (Nov. 1966),
806-809
[20] G. E. Forsyte, Today’s computational methods of linear algebra, SIAM
Review, 9, 3(July 1967), 489-515
[21] M. Fraser, N. Metropolis, Algorithms in unnormalized arithmetic III. Matrix
inversion, Numerishe Mathematik, 12, 5, 1968, 416-428
[22] D. R. Fulkerson, Philip Wolfe, An algorithm for scaling matrices, SIAM
Review, 4, 2 1962), 142-148
СПИСОК ЛИТЕРАТУРЫ
38
[23] Functional Mathematical Programming System CDC 6600, Bonner &
Moore Associates, Inc., Houston, Texas, 1971
[24] J. I. Golden FORTRAN IV programming and computing, Prentice-Hall,
Inc., Englewood Cliffs, New Jersey, 1965
[25] D. I. Good, R. L. London, Computer interval arithmetic: definition and proof
of correct implementation, Journal of the ACM, 17, 4(Oct. 1970), 603-612
[26] R. L. Graves, Philip Wolfe (Editors), Recent advances in mathematical
programming, McGraw-Hill Book Company, Inc., New York, 1963
[27] IBM Linear Programming System/360 (360A-CO-18X), Program
Description Manual (GH20-0607)
[28] IBM Program Product: Linear Programming System/1130 (LPS/1130)
Program Description (SH20-0633-1)
[29] IBM Application Program: 1130 Linear Programming - Mathematical
Optimization Subroutine System (1130 LP-MOSS)
[30] L. J. Larsen, A modified inversion procedure for product form of the inverse
- linear programming codes, Communications of the ACM, 5, 7(July 1962),
382-383
[31] L. S. Lasdon, Optimization theory for large systems, The Macmillian
Company, London, 1970
[32] H. Mueller-Merbach, On round-off errors in linear programming, SpringerVerlag, Berlin, 1970
[33] W. Orchard-Hays, Advanced linear-proggramming computing techniques,
McGraw-Hill Book Company, New York, 1968
[34] E. E. Osborne, On pre-conditioning of matrices, Journal of the ACM, 7,
4(1960), 338-345
[35] B. N. Parlett, C. Reinsch, Balancing a matrix for calculation of eigenvalues
and eigenvectors, Numerishe Mathematik, 13, 4, 1969, 293-304
[36] W. Sautter, Schwankungsbereiche linearer
Numerishe Mathematik, 21, 3, 1973, 264-269
Optimierungsaufgaben,
[37] A. van der Sluis, Stabilty of solutions of linear algebraic systems, Numerishe
Mathematik, 14, 3, 1970, 246-251
[38] A. van der Sluis, Condition, equilibrium and pivoting in linear algebraic
systems, Numerishe Mathematik, 15, 1, 1970, 74-86
СПИСОК ЛИТЕРАТУРЫ
39
[39] D. M. Smith, W. Orchard-Hays, Computational efficiency in product form
LP codes, in [25, 211-218]
[40] Computations with sparce matrices, SIAM Review, 12, 1970, 527-543
[41] J. Walish (Editor), Numerical analysis: An introduction, Academic Press,
London, 1966
[42] J. H. Wilkinson, Rounding errors in algebraic processes, H. M. S. O., 1963
[43] J. H. Wilkinson, Error analysis of direct methods of matrix inversion,
Journal of the ACM, 8, 3, 1961, 281-330
[44] Philip Wolfe, Error in the solution of linear programming problems, in L. B.
Rall (Editor), Errors in digital computation, Vol. 1, Willey, New York, 1964,
271-284
[45] Philip Wolfe, L. Cutler, Experiments in linear programming, in [25, 177200]
[46] Philip Wolfe, The composite simplex algorithm, SIAM Review, 7, 1, 1965,
42-54
[47] Philip Wolfe, A technique for resolving degeneracy in linear programming,
Journal of the ACM, 11, 1963, 205-211
Поступила в редакцию 20-ого мая 1974 года
Документ
Категория
Без категории
Просмотров
7
Размер файла
257 Кб
Теги
запятой, коларов, вычисления, переодоление, линейной, pdf, ошибок, программирование, 2006, плавающей
1/--страниц
Пожаловаться на содержимое документа