О рекуррентно-матричной форме метода оптимального исключения Текст научной статьи по специальности «Математика»

УДК 519.612

О РЕКУРРЕНТНО-МАТРИЧНОЙ ФОРМЕ МЕТОДА ОПТИМАЛЬНОГО ИСКЛЮЧЕНИЯ

© 2009 А. И. Жданов, Л. В. Яблокова Самарский государственный аэрокосмический университет

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

Прямой рекуррентный метод, прямой проекционный метод, метод оптимального исключения.

§ 1. Введение

В данной работе рассматривается задача решения систем линейных алгебраических уравнений (СЛАУ) с невырожденными матрицами, т.е. задача вычисления решения СЛАУ:

Ах = Ь, А е ЯпХп , Ь е Я” , гапк А = п . (1)

Прямые методы решения СЛАУ оцениваются по трем важнейшим качествам: числу выполняемых арифметических операций, требованию к объему оперативной памяти и максимальной точности, достижимой при их помощи.

Наилучшими по вышеперечисленным показателям методами решения невырожденных СЛАУ продолжают оставаться методы типа исключения Г аусса. Как известно [1] -[3], число операций этих методов оценивается величиной 2п3/з + о(п2), где в число операций включаются как умножения/деления, так и сложения/вычитания.

Гауссовское исключение существует во многих вариантах, которые алгебраически тождественны. Из всего множества вариантов исключения Гаусса метод оптимального исключения [1], [3] требует для своей реализации минимального объема оперативной памяти. Точно такие же характеристики имеет метод окаймления [1], [3], хотя он и не относится, вообще говоря, к методам исключе-

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

п величиной - п2 /4 машинных слов, а при нечетном п данный показатель будет близок к этому.

Как известно [4], в методе оптимального исключения попеременно выполняются операции прямого и обратного хода метода Г аусса. Это позволяет не вводить в оперативное запоминающее устройство очередную строку СЛАУ до тех пор, пока не преобразованы предыдущие строки, т.е. используется последовательный ввод строк матрицы и правых частей системы. Таким образом, метод оптимального исключения имеет рекуррентную структуру.

В настоящей работе рассматривается прямой рекуррентный метод решения СЛАУ, основанный на использовании формулы малоранговой модификации обратной матрицы (формулы Шермана-Моррисона [5]). Показана эквивалентность этого метода методу оп-

тимального исключения, хотя по своей алгоритмической структуре полученный рекуррентный метод не является методом исключения. Этот метод относится к классу прямых проекционных методов [6] - [8].

Указанная эквивалентность рассматриваемого прямого рекуррентного метода, прямого проекционного метода [8] и метода оптимального исключения дает возможность исследовать последний как конечный рекуррентный процесс, а также представлять его в компактной рекуррентно-матричной форме. Это позволяет получить оптимальные по требуемому объему оперативной памяти машинные алгоритмы решения СЛАУ со многими правыми частями, которые особенно эффективны в алгоритмах итерационного уточнения.

§ 2. Прямой рекуррентный алгоритм

Рассматриваемый в статье прямой рекуррентный метод решения СЛАУ [9] основан на формуле Шермана - Моррисона [5].

Л е м м а. Пусть дана матрица

В =

икк ик ,п-к

^°п -к ,к Оп -к ,п - к ^

где е Ккхк , к < п и гапкикк = к ,

ик п_к е Якх(п-к), Оп-кк и °п-кп-к - нулевые матрицы соответствующих размерностей.

Я

пхп

Введем обозначения.

Пусть і = 1, 2, • , п - номера строк матрицы А, аі = (аі1, ..., аіп) - строки матрицы А, т - знак транспонирования,

Аі =

аі

аі

V /

К'

іхп

Ап = А,

Ьі - элементы вектора Ь , Ь = (Ьі, ..., Ьп )6 є Кп,

А =

а11 • а1і

, аі1 • іа

ЯіХі, Аі = ап.

Ап = А.

Теорема. Пусть все главные миноры матрицы А отличны от нуля, т.е.

ёй Аі ф 0 "і = 1, 2, • , п,

(3)

и

пусть векторы Хі є Яп и матрицы

р. є ЯпХп , і = 1, 2, • , п удовлетворяют системе рекуррентных уравнений

Х+і = Х + (Ъ+1 - аг+1 Х )/ Щ+1, Хо = о, (4)

Тогда для всех |°| < о^іп ІЦкк) матрица р = р - а6 р Щ р = Е

' ' Аі+1 Аі 6і+1Ыі+11 і/ші+1 ’ ±0 ^п

(5)

рв = ііт(В + 8Е„ )8

о ®0

существует и равна

рв =

Окк

О

-и к!ик, п—к

•-к ,к

Е,

-к

(2)

где 8ш;пикк - минимальное сингулярное число матрицы икк, Еп - единичная матрица порядка п, к = 1, 2 , ..., п -1.

где р=ы^]е япхп, щ+1=а^

; = 1, 2, • , п -1.

Тогда условия (2) являются необходимыми и достаточными для того, чтобы

Щ+1 ^ 0 "; = 1, 2, • , п -1.

При этом хп = А_1Ь, т.е. хп является решением СЛАУ (1).

Пусть Аг = [Аг Аг,п-і ], где А є К

Лх(п-і)

1Х1 .

и

Аі п-і є V , тогда из условий (5) на ос-

новании леммы получаем, что "г = 1, 2, • , п матрицы Рг = Нш Р5 существуют и равны

5 ®0

О- А—1А.

О Е

,^-^П—І, І п—і у

(6)

Из (6) видно, что в матрице Р+1 первые г +1 векторов равны нулю, т.е. ё1+ = 0 "] = 1, 2, • , I -1. Следовательно, формулу (5)

можно записать в виде

І+1 і

ё , = ё ,

од

-ё

І+1 5

І+1

где

У, = а]+1ё), , = І + 2, • , п, І = 0, 1, -

(7)

а ёк = ек "к = 1, 2,

Из (7) непосредственно видно, что алгоритм (4), (5) является альтернативной формулировкой прямого проекционного метода [8].

Известно, что условия (2) не являются сильно ограничительными. Чтобы от них освободиться, достаточно выполнить перестановку, основанную на какой-либо стратегии выбора ведущих элементов. Как отмечается в [8], наиболее естественной и практичной для данного метода при решении СЛАУ с плотными матрицами продолжает оставаться стратегия с выбором максимального по модулю ведущего элемента по строке.

Формально эта стратегия на каждом г - м шаге состоит из следующих этапов:

1) определяется номер к из условия: к = шт N г , где

N = [к':

аІ+1 ёк'

■ тах

І+1<1 <п

аІ+1 ё1

2) выполняется перестановка векторов

ё\+1 и ёк .

В [8] отмечается, что вычислительные эксперименты свидетельствуют о том, что прямой проекционный метод с рассмотрен-

ной стратегией выбора ведущих элементов лишь незначительно уступает по точности методу исключения с использованием Ьи -разложения и частичным выбором (по строке) ведущих элементов.

Для СЛАУ с разреженными матрицами в [8] предложена новая стратегия выбора ведущих элементов. В этой стратегии номер к выбирается на г-м шаге из условия:

к = шт N., где

=[к': \а°+1ёк] > с ■ тгк \а-+1ё\ , 0 < с < 1}.

і і і +\<г<п ' і

В [8] показано, что эта стратегия в случае СЛАУ с разреженными матрицами при с < 1 позволяет значительно уменьшить заполнение матриц в процессе вычислений при незначительном снижении точности алгоритмов вида (6), (7).

§ 3. Анализ вычислительной схемы рекуррентного алгоритма

Нетрудно показать, что общее число арифметических операций, требуемых для реализации алгоритма (4), (5), оценивается

величиной 4п3 + О (п2). Однако, учитывая

специальную структуру (6) матриц Р , можно представить алгоритм (4), (5) в эквивалентном виде, для которого общее число арифметических операций оценивается величиной 2п3/3 + О (п2).

Для этого предварительно введем следующие обозначения:

а

Оп

Еп

Оі = (••• чП—і ) = — аі 1 Аі,п—і є &

{іх( п—і)

як'

& , к = 1,2,...,п-

я2 яз ••• Чп —і є &(і+1)х(п—і—1)

0 0 ... 0

для всех г = 1,2,...,п - 2.

Из определения матриц Qi и следу-

ет, что матрицы Qi+1 и (^. имеют одинаковые размерности, т.е.

(і+1)х( п —і—1)

а+ь (°і є я

для всех г = 1,2,...,п - 2.

Представим строки а° матрицы А в виде двух компонентов:

(аг )и е ^ и (аг )1 е Кп~1,

аІ =[(аХ°>+а У,]

а также введем вектор yi, состоящий из первых г компонент вектора хг,

Уі = (х{, •••, х-) є я1, где х = (х;, •••, х'п )0

"і = 1, 2, • , п .

С>1 =—( а12,....а1п )/

гУ і

(8)

(9)

і +1

і = 1,2,...,п — 1, п > 2,

У і+1 =

+

г(Ьі+1 — а 1+1 Уі)

од

І+1

(10)

(11)

і = 1,2,...,п — 1, п > 2, где

г =

Ч1

1

V /

я

і+1

аі+1 = (а/+1,1,-.>аі+1,і) є &

Уі+1 =(+1 )) ) +(+1)

++1 У.

Тогда, используя введенные обозначения, непосредственно получаем, что алгоритм (4), (5) можно записать в эквивалентном виде:

Таким образом, если выполняются условия теоремы, то у =х .

Г 5 * п п

Рекуррентная формула (5) эквивалентна рекуррентным формулам (8)-(9), а рекуррентная формула (4) - формулам (10)-(11).

Из формул (8)-(11) непосредственно видно, что вычислительная схема предлагаемого алгоритма практически совпадает с вычислительной схемой метода оптимального исключения [3]. Следовательно, число его арифметических операций оценивается величиной 2п3^3 + О (п2) , а для решения СЛАУ п - го порядка этим методом достаточно иметь

оперативную память машины порядка п2 /4 слов. Необходимая для машинной реализации алгоритма (8)-(11) оперативная память определяется максимальным размером мат-

рщы ^.

Таким образом, формулы (4), (5) (или, что эквивалентно, формулы (8) — (11)) представляют рекуррентно-матричную форму метода оптимального исключения.

§ 4. Решение систем со многими правыми частями

В [3] отмечалось, что единственным серьезным недостатком метода оптимального исключения по сравнению с Ьи - разложением является то, что он не позволяет эффективно решать системы со многими правыми частями, так как для его реализации требуется достаточно большой объем оперативной памяти, позволяющий запомнить все преобразования с матрицей системы. Использова-

ние рекуррентно - матричной формы алгоритма (формул (8), (9) и (10), (11)) дает возможность в значительной мере устранить этот недостаток классической вычислительной схемы метода оптимального исключения.

Подробно рассмотрим лишь общий случай СЛАУ с плотными матрицами.

Формулы (8), (9) и (10), (11) позволяют полностью разделить вычисления в методе оптимального исключения на два этапа аналогично использованию Ьи-разложения в методе исключения Гаусса. Первый этап — это обработка элементов матрицы А по формулам (8), (9). Второй этап — обработка правых частей и вычисление решений по формулам (10), (11). Точно так же, как и при Ьи-разложении, основная масса вычислений

(примерно 2п У 3 + О (п2)) приходится здесь

на вычисления по формулам (8), (9). Из формул (10), (11) видно, что для выполнения второго этапа нет необходимости запоминать на

первом этапе все матрицы Qi. Достаточно

запомнить лишь векторы д1, ...,^1п"1 и числа

щ, , где щ = а11. Они могут быть разме-

щены в верхней треугольной матрице

Я = (г ... Гп )е Я

пхп

где

Г = Ц,0,...,0 )6,

Г2 =

\о

Гп-1

-П—2

о

41 ) , )—1,0

Гп =

41 *) , щ

\о

е Яп.

Таким образом, для их хранения требуется оперативная память объемом п (п +1)12 .

В этом же массиве в процессе вычислений на первом этапе (по формулам (8), (9))

может храниться и текущая матрица Qi.

Для реализации второго этапа (вычислений по формулам (10), (11)) кроме верхней треугольной матрицы Я требуются элементы исходной матрицы А, расположенные ниже главной диагонали, т. е. элементы матрицы

Ьд =(. >,) ■

где

1у =

ау, апёе / > у,

0, апёе I < у.

Для хранения элементов матрицы Ьд

требуется п (п — 1)/2 машинных слов.

При этом возможны два варианта. В первом варианте элементы матрицы Ьд хранятся в оперативной памяти, т. е. при обработке элементов матрицы А на первом этапе происходит запоминание в оперативной памяти части элементов матрицы А в матрице Ьд.

Во втором варианте на первом этапе не происходит запоминания в оперативной памяти элементов матрицы Ьд и они сохраняются лишь на внешнем запоминающем устройстве.

В первом варианте для реализации ре-куррентно - матричной формы метода оптимального исключения требуется хранить в оперативной памяти машины элементы матрицы Р = ЬА + Я , т.е. требуется оперативная

память объемом п2 машинных слов. Следовательно, в этом варианте рекуррентно - матричная форма алгоритма не имеет преимуществ перед классической вычислительной схемой, основанной на Ьи-разложении матрицы А.

Во втором варианте требуется хранить в оперативной памяти лишь матрицу Я и, следовательно, требуемый объем оперативной

памяти — п (п +1)12 машинных слов. Таким

образом, в этом варианте рекуррентно - мат-

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

что снижает быстродействие этого метода по сравнению с аналогичным вариантом, основанным на Ьи-разложении.

Известно [2], что при итерационном уточнении возникает необходимость решения СЛАУ с несколькими правыми частями. При этом, если в итерационном уточнении используется Ьи-разложение, то предъявляются следующие требования к объему оперативной памяти машин: 2п2 машинных слов при хранении матрицы А в оперативной памяти и п2 при хранении матрицы А во внешней памяти. Это связано с тем, что для итерационного уточнения кроме Ьи-разложения требуется на каждой итерации еще и исходная матрица А.

Рассмотрим, как изменятся соответствующие характеристики, если в итерационном уточнении вместо Ьи-разложения использовать рекуррентно - матричную форму метода оптимального исключения, т.е. формулы (8)-(9). Если матрица А хранится в оперативной памяти, то отпадает необходимость

в матрице Ьд, и, следовательно, достаточно хранить только матрицы А и Я.

Таким образом, требуемый в этом случае объем оперативной памяти равен

п2 + п(п +1)12 = 3п2 /2 + п/ 2 машинных слов.

Если матрица А хранится во внешней памяти, то соответственно для итерационного уточнения требуется объем оперативной

памяти, равный п2 /2 + п/ 2 машинных слов.

Следовательно, использование в итерационном уточнении рекуррентно-матричной формы метода оптимального исключения предпочтительнее Ьи-разложения с точки зрения требуемой оперативной памяти при любом способе хранения исходной матрицы

А. При этом, как отмечалось в § 2, почти оди-

наково ведут себя оба метода с точки зрения устойчивости к ошибкам округления (при соответствующих стратегиях выбора ведущих элементов).

Очевидно, что в случае СЛАУ с разреженными матрицами данный алгоритм будет еще более эффективен по сравнению с LU-разложением с точки зрения требуемой оперативной памяти.

Библиографический список

1. Воеводин, В. В. Численные методы алгебры (теория и алгорифмы) [Текст] /

В. В. Воеводин. - М.: Наука, 1966. - 248с.

2. Воеводин, В. В. Вычислительные основы линейной алгебры [Текст] / В. В. Воеводин. - М.: Наука, 1977. - 303с.

3. Воеводин, В. В. Матрицы и вычисления [Текст] / В. В. Воеводин, Ю. А. Кузнецов. - М.: Наука, 1984. - 318 с.

4. Беклемишев, Д.В. Дополнительные главы линейной алгебры [Текст] / Д. В. Беклемишев. - М.: Наука, 1983. - 336 с.

5. Ортега, Дж. Итерационные методы решения нелинейных систем уравнений со многими неизвестными [Текст] / Дж. Ортега, В. Рейнболдт. - М.: Мир, 1975. - 558 с.

6. Фаддеев, Д. К. Вычислительные основы линейной алгебры [Текст] / Д. К. Фаддеев, В. Н.Фаддеева. - М.: Физматгиз, 1960. -656 с.

7. Abaffy J., Broyden C., Spedicato E. A class of direct methods for linear equations [Текст] / J.Abaffy, C.Broyden, E.Spedicato // Numer. Math. 1984. V.45. P.361-376.

8. Benzi М., Meyer C.D. A direct projection method for sparse linear systems [Текст] / М. Benzi, C. D. Meyer // SIAM J. Sci. Comput. 1995. V.16. N.5. Р1159-1176.

9. Жданов, А. И. Прямой последовательный метод решения систем линейных уравнений [Текст] / А. И. Жданов // Докл. РАН. - 1997. - Т.356, № 5. - С. 442 - 444.

10. Малышев, А. Н. Введение в вычислительную линейную алгебру [Текст] / А. Н. Малышев. - Новосибирск: Наука, 1991. - 228 с.

11. Гантмахер, Ф. Р. Теория матриц [Текст] / Ф. Р. Гантмахер. - М.: Наука, 1966. -576 с.

References

1. Voyevodin, V. V. Numerical methods in algebra (theory and algorithms) / V. V. Voyevodin - Moscow: Nauka, 1966 - 248 pp.

2. Voyevodin, V. V. Computing foundations of linear algebra / V. V. Voyevodin - Moscow: Nauka, 1977 - 303 pS.

3. Voyevodin, V. V., Kuznetsov Yu. A. Matrices and computations / V. V. Voyevodin, Yu. A. Kuznetsov - Moscow: Nauka, 1984 - 318 pp.

4. Beklemishev, D. V. Additional chapters of linear algebra / D. V. Beklemishev - Moscow: Nauka, 1983 - 336 pp.

5. Ortega, J. Iteration methods of solving non-linear systems of equations with many unknowns / J. Ortega, V. Reinboldt - Moscow: Mir, 1975 - 558 pp.

6. Faddeyev, D. K. Computing foundations of linear algebra / D. K. Faddeyev, V. N.

Faddeyeva - Moscow: Fizmatgiz, 1960 - 656 pp.

7. Abaffy J., Broyden C., Spedicato E. A class of direct methods for linear equations / J. Abaffy, C. Broyden, E. Spedicato // Numer. Math. 1984. V.45. P.361-376.

8. Benzi M., Meyer C. D. A direct projection method for sparse linear systems /

i. Benzi, C. D. Meyer // SIAM J. Sci. Comput. 1995. V.16. N.5. B.1159-1176.

9. Zhdanov, A. I. Direct sequential method of solving linear equation systems / A. I. Zhdanov // Report to the Russian Academy of Sciences. -1997 - vol. 356, No. 5 - pp.442 - 444.

10.Malyshev, A. N. Introduction to computing linear algebra / A. N. Malyshev -Novosibirsk: Nauka, 1991 - 228 pp.

11. Gantmakher F. R. Theory of matrices / F. R. Gantmakher - Moscow: Nauka, 1966 - 576 pp.

RECURRENT MATRIX FORM OF THE OPTIMAL ELIMINATION METHOD

© 2009 A. I. Zhdanov, L. V. Yablokova Samara State Aerospace University

The paper deals with the direct recurrent method [1] which is an alternative formulation of the direct projection method of solving systems of linear algebraic equations. The direct recurrent method is shown to be equivalent to the method of optimal elimination, which makes it possible to treat the latter as the direct recurrent method. This approach provides the potential for creating machine algorithms of solving systems of linear algebraic equations with many right sides, especially efficient in iteration refinement algorithms. The algorithms created have an optimal volume of on-line storage.

Direct recurrent method, direct projection method, method of optimal elimination.

Информация об авторах

Жданов Александр Иванович, заведующий кафедрой прикладной математики, доктор физико-математических наук, профессор, Самарский государственный аэрокосмический университет; e-mail: zhdanov@ssau.ru. Область научных интересов: вычислительная математика, матричные вычисления.

Яблокова Людмила Вениаминовна, ассистент кафедры прикладной математики, Самарский государственный аэрокосмический университет; e-mail: lyablokova@gml.com Область научных интересов: вычислительная математика.

Zhdanov Alexander Ivanovitch, Samara State Aerospace University, head of the department of applied mathematics, doctor of physical and mathematical science, professor, e-mail: zhdanov@ssau.ru. Area of research: computing mathematics, matrix computations.

Yablokova Lyudmila Veniaminovna, Samara State Aerospace University, assistant of the department of applied mathematics, e-mail: lyablokova@gml.com. Area of research: computing mathematics.