Компьютерное решение систем линейных уравнений с предельной точностью тема автореферата и диссертации по математике, 01.01.07 ВАК РФ

9; п

АКАДЕМИЯ НАУК СССР СИБИРСКОЕ ОТДЕЛЕНИЕ ВЫЧИСЛИТЕЛЬНЫЙ ЦЕНТР

На правах рукописи УДК 519.612.2

Кирилюк Олег Петрович

Компьютерное решение систем линейных уравнений с предельной точностью

01.01.07 - вычислительная математика

АВТОРЕФЕРАТ

Диссертация на соискание ученой степени кандидата физико-математических наук

Новосибирск-1991

Работа выполнена в Институте математики СО АН СССР.

Научный руководитель - член-корреспондент АН СССР профессор С.К.Годунов.

Официальные оппоненты: доктор физико-математических наук профессор А.Н.Коновалов, кандидат физико-математических наук доцента А.Д.Митченко

Ведущая организация: Институт вычислительных технологий

/ С >-1

Защита состоится 15 января 1992г. в 1 ч. I- мин. на заседании специализированного совета К 002.10.01 по присуждению ученой степени кандидата наук в Вычислительном центре СО АН СССР по адресу: 630093, Новосибирск-90, пр. академика Лаврентьева, 6.

С диссертационной работой можно ознакомиться в читальном зале Вычислительного центра (Новосибирск-90, пр. академика Лаврентьева, 6).

Автореферат разослан " 76 " • 1991г.

Ученый секретарь специализированного совета кандидат физико-математических наук

1

Ю.И.Кузнецов

... < I. Общая характеристика работы

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

Цель работы. Построение эффективных численных алгоритмов решения систем линей.чых алгебраических уравнений с полностью заполненными матрицами коэффициентов дающих гарантированную сценку точности результата. Прл этом требуется, чтобы гарангаровалвзя оценка точности результата работы алгоритма была близка к экстремально возможно;». Научная новизна:

- Получены эффективные алгоритмы решения систем линейных алгебраических уравнений с квадратными и прямоугольными матрицами коэффициентов лоляого ранга. Гарантируемая точность вычисляемого решения получается близкой к точности представления вектора в памяти компьютера.

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

- Получен алгоритм вычисления г-решения системы в случае матрицы коэффициентов системы неполного ранга.

- Получены новые теоремы об упрощении структуры ленточных матриц специальными ортогональными преобразованиями.

Апробация работы. Основные результаты диссертации докладывались на семинарах ОВМ АН СССР (рук. чл.-хорр. АН СССР В.В.Воеводин, д.ф.-м.н. Ю.Л.Кузкецов), ВЦ СО АН СССР (рук. д.ф.-м.н. В.П.Ильин, д.ф.-м.н. А.Н.Коновалов), ИМ СО АН СССР (рук. чл.-корр. АН СССР С.К.Годунов).

Публикации. Основные результаты диссертации опубликованы в трех работах, приведенных в конце автореферата.

Структура и объем работы. Диссертационвая работа изложена на 112

страницах и состоит из введения и двух глав. Работа включает 4 таблицы и список литературы, содержащий 22 наименования.

2. Содержание работы

Во введении (§1) излагается постановка задачи решения с гарантированной точность систем линейных алгебраических уравнений на компьютере. Приводится кркткий обзор основных результатов полученных в данной области. Постановка задачи состоит в следующем. В нашем распоряжении имеется некоторый конкретный компьютер разрядная сетка которого и его арифметические операции характеризуются известным набором констант. Пусть теребуется решить систему линейных алгебраических уравнений на данном компьютере. В диссертации описываются и обосновываются алгоритмы решения систем, которые анализируют ту систему уравнений которую требуется решить и в зависимости от результата этого анализа получают один из следующих двух ответов: 1) приближенное решение системы с гарантированной оценкой его точности; 2) отказ от решения исходной системы как системы полного ранга с выдачей числа обусловленности системы или параметра ее несовместности как количественного подтверждения плохой поставленное™ решаемой задачи. При этом, в отличае от ранее существовавших алгоритмов, гарантируемая точность решения в случае варианта 1 близка к экстремальной точности решения, т.е. к точности представления вектора в памяти компьютера.

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

Известно, что произвольную (Ы х М>матрицу А можнопредставить в виде Л - У-К-и\ (1)

где Ц и V - ортогональные матрицы, а К имеет вид

К '

[ лг] при N > М I при N - М, X [X ; 0] при N < М

о

о

s aL, L— mln(N,M)-

Определение 1.2. Разложение вида (1) матрицы А в произведение двух ортогональных V и V и одной "диагональной" К будем называть сингулярным разложением А. (Слово диагональный взято в кавычки в силу прямоугольности матрицы Л).

Определение 1.3. Числа а , ... , oL будем называть сингулярными числами матрицы А.

Определение 1.4. Пару векторов и<0 и будем называть соответственно правым и левым сингулярными векторами А соответствующими ее сингулярному числу <т(А) если существует такое сингулярное разложение матрицы А, при котором и(" и \fn будут коллннеарны соответственно í-му столбцу матрицы U и г-му же столбцу V, а на fZw'+l^-м месте диагонали ÍT стоит «гМЛ Если i > L, то для N > М вектор (или и^ для N < М) будем также называть сингулярным вектором Л, если существует такое сингулярное разложение (1), что он коллинеарен соответствующему столбцу V (или V), несмотря на то, что ему не соответствует парного вектора и(" (и(") и сингулярного числа а(А).

Пусть задан некоторый положительный параметр г, 1 s г S min(N,M) и пусть сингулярное разложение матрицы А имеет вид (1). Определим матрицы ранга г.

О

£-г+/

О

О

к

ítv] 2

Г

[2?г!0]

при N > М при N « М при N < М

А= V-K -U .

Определение 1.5. Обобщенным нормальным r-решением системы

(1) (или просто г-решением) будем называть обобщенное нормальное решение системы

Ах'= /

г

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

В первой главе (§§2-7) приводятся и обосновываются алгоритмы решения систем линейных уравнений и вычисления базисов ядер и г-ядер матриц.'

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

Пусть требуется решить систему

Ах - / (2)

с квадратной ЫхМ матрицей А ранга N и некоторой правой частью /. И пусть в нашем распоряжении имеется некоторый метод приближенного решения этой системы с гарантированной точностью. Приближенное решение системы (2), полученное этим методом, обозначим через х ш Ав/. В случае, когда этот метод гарантирует оценку погрешности

Их -х\\ 5 г 1М1, IIЛе/ - А-'/[{ * С \\А-% причем константа е не зависит от правой части f системы (2), в диссертации описывается простой итерационный процесс уточнена! решения х.

Введем несколько обозначений. Символами о, е и в будем обозначать соответственно машинные операции умножения матрицы на вектор, вычитания и сложения двух векторов, допускающие следующие оценки точности |И©* - ¿Ы з И1-М, 1К* в у) - (X - 5 с, II* - >11 + сЛГМ + ЫЛ

||Г* ® у) - (X + УЛ1 5 «;||х - у|| + с2е0(\\х\\ + !|у||л где с! и с2, е/ и £д - константы, причем с/ и с2 зависят от размерности N матрицы А и вехторов х и у, а к ед - от разрядной сетки ЭВМ, на которой проводятся вычисления (число г1 характеризует относительную точность вычислений, а ед - абсолютную). Тогда предлагаемый алгоритм решения системы (2) выглядит следующим образом.

Алгоритм 2

1. Приводим матрицу А к двухдиагональному виду (или другому виду, позволяющему быстро вычислять Авц для любого вектора г), вычисляем и, ||Л-'й (под нормой матрицы здесь и ниже подразумевается ее операторная норма), определяем характеристику точности решения системы г, величины а п /2 по формулам

а = {е + /"г + с^П+е+с^; + с2еоМА)(Ш)}х(Не+с2ео) + 2 сд,

/? = « + /с^П+е + 2сд;дМ;хГНе+с2£о/)Г1+г> + сд.

а также <5 согласно

3 = с;«/1+£ + сЛ; + сд 1 В ъ

Если а > ^ или > ^'Е^ то процесс итерационного уточнения

заканчивается отказом от вычисления решения. Иначе полагаем

/ог/2сг 0/(1-а))

к =

и переходим к шагу 2.

2. ;о= лэ/, 1=о.

3. у = <4 о л: .

| _1

4. г = /ву(.

И^'Н',. II/!!*,.!

5. р----.

' 1-е,- 1К'II•«",№,« -г-лм;

6. Если р > 2«;, то перейти к шагу 7, иначе считаем расчет законченным с результатом

7. 2 = Л /• .

> (

8. <7.+/= •+ /?.

р- -:-—:-;-•

1 -е,- (\\z.W ■(1-е1)/(1-е)+ е 1\\А-,\\-\\гр/\\х!\\- З-^СА)

9. х = х ® г .

и-1 1 I

10. Если г+1 <к или 9 +;>2£(, то перейти к шагу 2, заменив / на ¿+1, иначе считаем расчет закопченным с результатом +

Третий параграф посвящен задаче решения систем линейных уравнений с

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

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

в случае М>М и к

Г"* Л А' о

О А

3

(3)

(4)

А {О)

при Ы<М. При этом параметр р входящий в матрицы расширенных систем выбирается равным <г{А)Н 2 ', где <г;М,М/||Л~'||. Далее для расширенных систем применяется алгоритм изложенный в §2. Даказано, что итоговая точность х-ых компонент полученного решения будет не хуже чем

в; - 4 * V ыз+ ы2 *

в случае Ы<М и

II* " *11 5 и \ 1 + 2М^ЛГ.М - 2е /1 + 2р'М/>"М

'4 о\(А)-||*И

в случае Н>Ы. Величина « вычисляется в процессе

решения системы.

В четвертом параграфе проводится детальный анализ погрешностей решения расширенных систем (3), (4) и формулируется итоговый алгоритм решения систем с прямоугольной матрицей коэффициентов полного ранга.

В пятом параграфе предлагается для вычисления обобщенного нормального г-решения хг использовать квадратную расширенную систему вида

(5)

в которой столбцами матриц 1/д Уд являются ортогональные базисы из М-г правых и Ы-г левых сингулярных векторов соответствующих наименьшим сингулярным числам матрицы А. Эти базисы в дальнейшем будут называться базисами г-адер А и А, а матрица системы (5) - через В. Обозначим через

ЪУ0 А- Р1- /

. 0 К м 0

сингулярные числа матрицы А

(а[А)-/ А ¿А А)). Положим (5) решить систему

+ а )

Ьт+1 1у-г

упорядоченные по возрастанию Тогда если вместо системы

о ри

где Уд и ио - приближенно вычисленные базисы К0 и 1/0, такие, что

а

IV У0И 5 «1*;1 и «V ".и 5 £11уо»

то ,--

|я -х\&< 2в|х||*41 + »УЛ/) ,

где через и (А^ обозначена величина

аналогичная параметру о введенному в §3 характеризует несовместность г-решения системы (2). Таким образом задача вычисления г-решения системы сводится к решению расширенной системы (6) и к вычислению с по ВОЗМОЖНОСТИ ВЫСОКОЙ ТОЧНОСТЬЮ матриц Ч^У столбцы которых являются приближенными ортогональными базисами из Ы-т правых и Ы-г левых сингулярных векторов соответствующих наименьшим сингулярным числам матрицы А.

Е седьмом параграфе приводятся н обосновывается алгоритм итерационного уточнения базиса г-ядра вообще говоря произвольной матрицы, а в шестом параграфе подробно разбирается важный частный случай т"т1п(Ы,М). Алгоритм уточнения г-ядра в общем случае имеет следующий ввд.

Алгоритм б

1) Пусть известен почтя ортогональный базис подпространства близкого к г-ядру матрицы А. Матрицу столбцами которой являются векторы этого базиса обозначим через У (17 размера Мх(М-г)). Пусть точность этого базиса задается двумя константами у^ и у^ смысл которых следующий: 1а) существует матрица О с ортогональными столбцами (17*17 ■ Г) такая, что

II" - щ * с

16) верна оценка

* г'2°>

где через П1 обозначен оператор ортогонального проектирования на г-ядро матрицы А.

Вычислим сингулярные числа с, , аг. а, матрицы А и определим по ним следующие величины:

р - <^/3. УГ^Фу-

а - {V [[^ + 1] 4 - 1) г,]}/1 +

' ' [№ " ') + 1 + < - 2у + а, <1 - 3 + у?,

,--£_

1 + 1 -«- <4°

1 - { - й-у

+ <5 ,

,<о> О'

С1 + с - о - /а + й-дд- с - </

4Л(а +

2(1

+ а-гв- с - о + /п 4 с - - 4а/а +

уз . _

ще 8„ и 6 „ характеризуют точность проведения (Зй-разлож'зния и зависят V 2к

от размеров матрицы А - М и # и от параметра е(, характеризующего разрядную сетку компьютера на котором проводятся вычисления. Проверим выполнение неравенств

г + й-у'0)< 1, п < 1, /"-'< у

>2 ' ' '< _>2 2 тах

с + ( + 2ЛГУ а + 5 (1 - О - Л д0< 1-

Если хотя бы одно из неравенств не выполнено - отказываемся от вычисления базиса г-ядра ввиду плохой обусловленности этой задачи. Если же все неравенства выполнены, то дополнительно проверяем условие

у<°>> у '2 'г ш(л

Если последнее условие не выполнено, го начальное приближение уже достаточно хорошее. Поэтому заканчиваем алгоритм и в качестве решений выдаем (/ и две константы у'^ и

Если же все приведенные выше неравенства выполнены, то выберем к равным наименьшему целому I среди /2:0 для которых выполнено неравенство

С

гЬ

) £

Положим, далее, р»1; 2) Пока выполняем следующее:

V и-

ж; т рр 1 н 1-1

■ л -Р1» У> о.

2б) Выполним (Зй-разложение матрицы X. Ортогональную матрицу этого разложения обозначим через а треугольную - через Я\ А"= /^й, -прямоугольная матрица тех же размеров что и X. причем I).

I II М-Г

2в) / увеличим на 1.

3) Результатом работы алгоритма считаем Р и константы у*= <5. и • . * ' о

у " 2{/Д - ц) характеризующие его точность.

Значения величин 6 и 6 зависят от реализации (Ж-разложения. В с С*

частности, нетрудно предложить реализацию при которой

С

Г

2 ты

5 = S = QR Q

23-"ml/^ 'A1 - + ной точности

использовании удвоен-

(4-п +26)-п У п £,

max irtin ml п 1

при использовании обычной точности

где п = min(N,M), п " max(N,M).

mm max

Вторая глава посвящена задаче упрощения структуры ленточных матриц ортогональными преобразованиями специального вида называемыми цепЬчками двумерных вращений.

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

Пусть квадратная (Nx ^-матрица А с вещественными элементами а^ имеет следующую структуру:

о

А =

Ы-к+1,1

N-k+I ,N

(7)

Л7 NN

т.е. элементы матрицы А удовлетворяют условию а= 0 при ¡'-/Н г 0. Определение 8.1. Матрицу А, имеющую структуру (7), будем называть Л-нижнегессенберговой.

Определение 8.2. Матрицу А, сопряженная к которой А является Л-нижнегессенберговой, будем называть ¿-верхнегессенберговой. Пусть матрицы С (2й г'зЛг) имеют вид

О

о

¡-я строка,

где с2+ sj= 1.

Определение 8.3. Матрицу

С = <У с>

будем называть цепочкой двумерных вращений, определяемой параметрами с,, 5, (25/'^Л0.

Определение 8.4. Будем говорить, что цепочка С поляризует

вектор ц~Си^ иг, ... , иы)т, если параметры с, определяющие цепочку

двумерных вращений С, подобраны так, что имеет место равенство

О "

Си

О

" ±1Нк"

О

±ии

Определение 8.5. Преобразование матрицы А по формуле САС , где С и С - цепочки двумерных вращений, поляризующие соответственно левый и правый сингулярные векторы А, будем называть сингулярным исчерпыванием матрицы А. ^

В девятом параграфе доказываются три основные теоремы об упрощении структуры Л-ншкнегессенберговых матриц.

Теорема 9.1. Пусть дана к-нижнегессенберговая матрица А, и векторы и и V связаны равенством

Аи-у, (8)

причем

+ < + -+

Подберем иепочку С, поляризующую вектор и, и С, поляризующую V. Тогда матрица А«<МС?, будет к-нижнегессенберговой и, кроме того, первые (УУ-1) элементов последнего столбца'А будут нулевыми.

Теорема 9.2. Пусть дана г-нижнегессенберговая матрица А и заданы

некоторые векторы и и V, причем

и2 + и2 + ... + и1 # 0. 1 2 м

Подберем цепочку С поляризующий, вектор и и С поляризующую V. Предположим, кроме того, что выполнено п (ЛЧ:+1 йпйМ первых соотношений (8). Тогда матрица А будет к-нижнегессенберговой и первые п-1 компоненты ее последнего столбца - нулевые.

Теорема 9.3. Пусть дана к-нижнегессенберговая матрица А (к> 1), и вектор и удовлетворяющий условию

и2 + и2 + ... + и2 * 0. Подберем цепочку С, поляризующую вектор и. Тогда матрица А-Ж?, будет (к-\)-нижнегессенберговой. Если же дополнительно п-я (Ы-к+1<п£Ы) компонента вектора Ли нулевая, то нулевой будет и п-я компонента последнего столбца Д.

В десятом параграфе теоремы об упрощении структуры к-нижнегессенберговых матриц применяются к задаче сингулярного исчерпывания двухдиагональных матриц. Теорема 10.1. Пусть матрица Р имеет вид

о -

Ч Ь2

<1 ъ

з з

О

о

а ь

Я-Г У

(9)

а ее элементы удовлетворяют условию

¡¡¡Ф О (1 Ь^ О (1(10)

Тогда, если цепочки С и С поляризуют векторы и и V удовлетворяющие равенству

Ли " V,

а элементы векторов и и у удовлетворяют условиям

и, + и* Ф О, V * О, то матрица ТН2ЕЮ примет вид

где

о - к+

2 и

К'

•О п, О п.

О

О

0 V, о о

. А

X

X X

X X X X

а п и £ задаются формулами 1 у

Ъ.

(г-2,...Л-1), к'

й с

Н У

о

X

X О

(11)

Теорема 10.2. Пусть цепочки С и С поляризуют векторы и и V удовлетворяющие равенству

б'у - и,

где элементы даухдиагональной матрицы О вида (9) удовлетворяют условиям (10), причем

v¡Ф О.

Тогда матрица имеет вид

где

КГ

к.

О

о

Я-1

о

, 7'

О X X О X

а элементы £ задаются формулами

X X X X X X

О X X

п

</ 5

I+1 1+1

(¡=1,2.....ЛЧ).

(12)

Теорема 10.3. Пусть цепочки С и С поляризуют правый и левый сингулярные векторы, соответствующие сингулярному числу а двухдиагоналыюй матрицы О элементы которой удовлетворяют условиям (10). Тогда матрица Ь=СЪс имееп вид

Ь = кг

0

О

с элементами к, >1 задаваемыми (11), (12). Причем

' кн '

Теорема 10.4. Пусть элементы двухдиагональной матрицы £> вида (9) удовлетворяют условиям

Эг, ¿=0; Ь.ФО (г2.....ЛО.

Тогда существует исчерпывание нулевого сингулярного числа матрицы 1) такое, что

Г*,

к2 П3

о = сое =

О

О

V/0

где

к =

Л * , ¡+1 ¡+1

Ь 'с

Ш Г

если с! & О, ¡+1

если (! = О, ¡+1

(¿=1,2,... ,N-1)

Ь я

¡+1 I

(¿=2,...,Лг-1)

Одиннадцатый параграф посвящен анализу погрешностей процеду исчерпывания двухдиагональных матриц. Эта процедура является базовой задаче построения базиса г-ядра произвольной матрицы. Теорема 11.1. Пусть цепочки С и С поляризуют векторы и и удм.и'ткорнннцне уравнениям

1)и = сп', Ь V = а и

и.+ и, £ о, V* 0.

12 N

Причем элементы двухдиагональаых матриц О и В вида

гг. У

D

■а, ъ2 hh

о

О

, ъ

О

О

а ъ

N-1 N

таковы, что

tf- (l+S^d i

fl+SlJJl

(lsi

N

йЛО,

ai-(i+fil)dl <2sj's«,

maxC | <5,1 , !<5(1,Ц1,Ц|}

Л элементы матрицы 3 задаются согласно 3« (\ +<?.)• -11LJ-L! (1 s isN-l),

d с

NN

l* 1

b s

где max{ | I , | I } S р. Тогда

. _ 2/2^ (/TT + Ua + VIT?

BCDC-DJ S -—-И|.

1 - V 2 (<x + q>)

ПоэтЬму итоговую оценку точности исчерпывания можно преобразовать к виду ||ÜDC*- ВЦ S (ISN/JF + 15/V - s/~N + SMjßH.

Список работ опубликованных по теме диссертации

1. Кирилюк О.П. Итерационное уточнение решения систем линейных уравнений. -Новосибирск, 1989. (Препринт /АН СССР. Саб. огд-кне, Ия-т математики; а 7).

2. Кирилюк О.П. Упрощение структуры ленточных матриц. / Численный анализ, тр. Ин-та математики, т. 15. - Новосибирск: Наука, 1989.

3. Киршпок О,П. Итерационное уточнение решения систем линейных уравнений. / Вычислительные методы линейной алгебра, тр. Ин-та математики, т. 17. -Новосибирск: Наука, 1990.

1