Approximation approach to the solution of gravity and magnetic problems. 2. New methods of getting stable approximate solutions of linear algebraic equations with their approximated right sides Текст научной статьи по специальности «Математика»

Аппроксимационный подход к решению задач гравиметрии и магнитометрии.

II. Новые методы нахождения устойчивых приближенных решений систем линейных алгебраических уравнений с приближенно заданной правой частью

В. Н. Страхов, А. В. Страхов

Объединенный институт физики Земли им. О. Ю. Шмидта РАН

Введение

В первой части настоящей работы [Страхов, Страхов, 1999] были рассмотрены две темы:

1) сущность классического подхода [см. Воеводин, 1969; Иванов и др., 1978; Лаврентьев, 1962; Лисковец, 1981; Морозов, 1973, 1974, 1987; Танана, 1981; Тихонов, 1963a, 1963b, 1965; Тихонов, Арсенин, 1979; Тихонов и др., 1985] к регуляризации систем линейных алгебраических уравнений с приближенными данными;

2) сущность новой, разработанной В. Н. Страховым (частично в соавторстве, см. [Страхов, 1990a, 1990b, 1991a, 1991b, 1991c, 1992a, 1992b, 1995a, 1995b, 1996a, 1996b, 1996c, 1997a, 1997b, 1997c, 1997d, 1997e, 1997f, 1997g, 1997h, 1998a, 1998b, 1998c, 1998d, 1998e; Страхов, Страхов, 1998, 1900, 1900, 1900, 1900; Страхов, Тетерин, 1991a, 1991b, 1991c, 1993; Strakhov et al., 1995]) теории нахождения устойчивых приближенных решений систем линейных алгебраических уравнений с приближенными данными.

Однако вторая тема была рассмотрена в [Страхов, Страхов, 1999] только с самых общих позиций. Более глубокому раскрытию этой темы и посвящена вторая часть работы. Таким образом, рассмотрение собственно основной темы, вынесенной в заголовок работы, а именно - сущности аппроксимационного подхода

©1999 Российский журнал наук о Земле.

Статья N КЛВ99021.

Онлайновая версия этой статьи опубликована 15 ноября 1999. ИЯЬ: http://eos.wdcb.rssi.ru/rjes/RJE99021/RJE99021.htm

Это может показаться парадоксальным, но вся наука подчинена идее аппроксимации.

Бертран Рассел

Умный начинает с конца, дурак кончает в начале.

Дьердь Пойа

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

План изложения в настоящей, второй, части работы таков. Сначала в первом разделе рассматриваются вопросы ортогональных преобразований матриц и векторов, в частности - целый ряд новых конструкций ортогональных преобразований, разработанных В. Н. Страховым, частично в соавторстве с А. В. Страховым. Рассмотрение этих ортогональных преобразований имеет принципиальное значение, по крайней мере по двум причинам:

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

б) классические алгоритмы ортогональных преобразований, основанные на использовании гивенсовых вращений [см. Фаддеев, Фаддеева, 1963; Воеводин, Кузнецов, 1984; Уилкинсон, Райнш, 1976], как правило “отказывают” в случае линейных систем большой размерности - при использовании стандартного 64-битового машинного слова.

Во втором разделе работы обсуждаются наиболее перспективные (по мнению авторов) методы (алгоритмы) нахождения так называемых “пробных” приближенных решений линейных систем - для ситуации, когда приближенно известна лишь правая часть

354

СТРАХОВ, СТРАХОВ: АППРОКСИМАЦИОННЫй ПОДХОД К РЕШЕНИЮ ЗАДАЧ II

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

Далее в четвертом разделе статьи рассматривается проблема, которая для гравиметрии и магнитометрии имеет в некотором смысле центральное значение - каким именно способом находить разумные приближенные решения линейных систем (в особенности - систем большой размерности), в ситуации, когда отсутствует “привычная” информация о векторе 6/ аддитивной помехи в задании правой части системы, т.е. неизвестны вообще постоянные в неравенствах

5/

< 52

п — тах

либо величина

б2.

тт

62

тах

очень мала, условно

в2 < -.

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

§ 1. Начнем со следующего общего замечания -в данном разделе статьи рассматриваются ортогональные преобразования систем линейных алгебраических уравнений

Ах = / = / + 5/

(4)

(1)

(2)

(3)

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

Наконец, в “Заключении” дается общее резюме по первым двум частям работы и намечается общая программа действий по реализации новой теории нахождения устойчивых приближенных решений систем линейных алгебраических уравнений с приближенными данными.

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

где А суть (МхМ)-матрица с вещественными элементами а^, 1 < г < N, 1 < ] < М, (М - число неизвестных, М - число неизвестных), / - заданный М-вектор (/ - вектор полезного сигнала, 6/ - вектор помехи), х - искомый М-вектор; ортогональные преобразования всегда имеют цель упростить процедуру нахождения устойчивого приближенного решения, согласованного с имеющейся априорной информацией об х (или /) и 6/. При этом понятие “упрощения процедуры нахождения устойчивого приближенного решения” достаточно сложное. Именно, возможны и фактически реализуются два подхода:

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

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

Достаточно очевидными являются следующие моменты:

1) прежде всего надо дать общую классификацию ортогональных преобразований - по структурному принципу;

2) далее необходимо выделить смысловые (содержательные) классы ортогональных преобразований;

3) наконец, необходимо достаточно детально описать содержательные классы ортогональных преобразований.

В данном разделе статьи (имеющем в некотором смысле ключевое значение) последовательность подачи материала такова: в § 2 приводится структурная классификация ортогональных преобразований; в § 3 выделяются содержательные классы ортогональных преобразований, а в §§ 4-8 описываются преобразования из основных классов; наконец, в § 9 обсуждается принципиальная проблема устойчивости ортогональных преобразований.

2

2

в

§ 2. В рамках новой теории регуляризации систем линейных алгебраических уравнений рассматриваются три основных типа ортогональных преобразований систем линейных уравнений (4):

I тип: приведения системы (4) к виду

Bx

PS

где

B = VTA, ps = VTfS

где

Az = fs

A = AU, z = UJ

где и есть ортогональная по столбцам (МхМ)-мат-рица.

III тип: приведения системы (4) к виду

Cw = фs

A = ATj(ф), 1 < i<j < M ,

z = TT(ф)х, 1 < i < j < M .

(14)

(5)

(6)

Наконец, элементарное ортогональное преобразование системы (4) 3-го типа:

при этом

и V есть ортогональная по столбцам (МхМ)-матрица. II тип: приведения системы (4) к виду

(7)

(8)

A w = фs

° т

A = TT (p)ATj (ф) ,

1 < i < j < M, 1 < j' < j' < M

w = Tj (ф)х ,

^s = tT (p)fs .

(15)

(16)

(9)

В приведенных соотношениях ((12), (14), (16)) Т^-(р), Т^-(-0) и т.д., суть элементарные матрицы плоского вращения, а именно - матрицы со следующими элементами

1) Т^(а): элементы ,

где

C = VTAU, фs = VTfs,

= UT a

(10)

а V и и суть ортогональные по столбцам матрицы, размерностей (МхМ)- и (МхМ)-, соответственно.

Ясно, что это лишь классификация самого общего вида; в каждом из трех выделенных типов имеются подтипы. Здесь важно сразу же подчеркнуть, что ортогональные матрицы V и И, используемые в преобразованиях (5)-(6), (7)-(8) и (9)—(10) линейных систем, всегда конструируются в форме произведений элементарных матриц плоского вращения. Иначе говоря, общие ортогональные преобразования (5)-(6), (7)-(8) и (9)-(10) конструируются в виде последовательностей элементарных ортогональных преобразований трех типов.

Во всех случаях за исходную принимается система

(4).

Элементарное ортогональное преобразование системы (4) 1-го типа:

tp,q *

cos a, p = q, p = i, j ;

- sin a, p = i, q = j ; sin a, p = j, q = i ;

1 p = q, p = ьj ;

0, для всех остальных (p, q).

(17)

2) TT(a): элементы tp q ,

tp,q ^q,p

(18)

где суть элементы матрицы Т^-(а).

Пусть аи,т, о^т, ак т, а¿,,т суть элементы матриц

Л V ◦

А, А, А, А в соотношениях (4), (12), (14) и (16). Тогда: 1) в матрице А по сравнению с матрицей А изменены только г-ая и у-ая строки, согласно соотношениям

a(i) = v(i) cos p + v(i) sin p ,

a(j) = —v(i) sin p + v(j) cos p (19)

A x = fs

A = TT (p)A, 1 < i<j < N

fs = TT (p)fs,

1 < i < j < N .

(11)

(12)

Элементарное ортогональное преобразование системы (4) 2-го типа:

(где а(г) и суть г-ая и у-ая вектор-строки А, а(г)

и а^ - соответственно г-ая и у-ая вектор-строки А);

V

2) в матрице А по сравнению с матрицей А изменены только г-ый и у-ый столбцы, согласно соотношениям:

v i-i

a(í) = a(í) cos ф + vj) sin ф ,

v • /о. / (20)

a(j) = —v(i) sin ф + v(j) cos ф

v

Av z = fs

V V ^

(13) (где и а^-) суть і-ьій и ^-ый вектор-столбцы А, «(¿)

x

w

и - соответственно г-ый и у-ый вектор-столбцы

А);

о

3) в матрице А по сравнению с матрицей А изменены только г-ая и у-ая строки и ¿'-ый и у'-ый столбцы, согласно соотношениям:

° (*) * (j) * •

a = a(t)’ cos p + au)’ sin p,

° (j) (*) * • (j) *

a = —a(i)’ sin p + au)’ cos p

a**') = a*i') cos ф + aj sin ф ,

°

a*j') = —a**') sin ф + a*j') cos ф

(21)

(22)

* . * V V

в которых а(г)’ , а ^ , а(,), а*.,), соответственно

а(г)’ , а^’ , а*4,), а*,), суть векторы-строки и векторы-столбцы с выброшенными элементами, находящимися на пересечениях строк и столбцов (в позициях (г, г'), (г,у'), (г', г) и (г',у), соответственно); для элементов в указанных позициях выполняются следующие соотношения:

° ai,i' = (ai,i' cos p + aj,*' sin p) cos ф

+ (ai,j' cos p + a¿j' sin p) sin ф ,

° ai,j' = —(ai ,*' cos p + aj,i' sin p) sin ф

+ (ai,j' cos p + a¿j' sin p) cos ф ,

° aj,i' = (—ai ,*' sin p + aj,*' cos p) cos ф

+ (—ai ,j' sin p + ajj/ cos p) sin ф ,

° aj,j' = —(— a*,*' sin p + aj ,*' cos p) sin ф

+ (—ai ,j/ sin p + a^y cos p) cos ф •

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

Замечание 1. В общем случае ортогональные по столбцам матрицы V и и (размеров (МхМ)- и (МхМ)-, соответственно), фигурирующие в (12), (14), (16), представимы в форме произведения элементарных матриц плоского вращения,

Замечание 2. Существуют и другие виды ортогональных преобразований, например, осуществляемые с помощью так называемых матриц отражения [см. Воеводин, Кузнецов, 1984; Уилкинсон, Райнш, 1976; Фаддеев, Фаддеева, 1963]; однако в рамках развиваемой новой общей теории регуляризации они не рассматриваются, ибо подобные матрицы всегда могут быть представлены в форме конечного или бесконечного произведения элементарных матриц плоского вращения.

Замечание 3. Если в структурной классификации ортогональных преобразований основной упор делается на вопрос о структуре произведений элементарных матриц плоского вращения (24)-(24’) - является оно конечным или бесконечным, и каким именно образом задается последовательность пар индексов (i, j) у матриц элементарных плоских вращений, то при выделении содержательных классов преобразований основным является вопрос задания параметров преобразований (cosр, sinр, cosф, sinф и т.д.), из каких именно условий они находятся; последние же (условия на параметры) определяются общей целевой установкой на конечный итог преобразований. Замечание 4. Конечно, очень заманчива мысль объединить структурный и содержательный подходы в некоей единой классификации ортогональных преобразований. Но по мнению авторов эта цель почти не реализуема - настолько сложной окажется искомая единая классификация.

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

В рамках типа I выделяются три основных подтипа:

11) - в нем матрица V представима непосредственно в форме (24) с конечным числом сомножителей, без какого-либо дополнительного структурирования; и каждая из матриц Vp есть некоторое конечное произведение элементарных матриц плоского вращения;

12) - в нем матрица V имеет следующую структуру:

V = П V •

p=1

(25)

(24)

(24’

и каждая из матриц V есть некоторое конечное произведение элементарных матриц плоского вращения; ясно, что числа сомножителей в представлениях матриц могут быть различными;

13) - в нем матрица V имеет следующую структуру:

которое может быть конечным или бесконечным; в последнем случае важно понять структуру бесконечного произведения.

V = П V •

p=1

(26)

Иначе говоря, если ввести последовательность

преобразований систем (4)

А = Ао ,

Ар = УТАр-1 , (27)

р = 1, 2,... ,

то в случае ортогональных преобразований подтипа

12) индекс р имеет лишь значения, определяемые предельным значением Р, а в случае ортогональных преобразований подтипа 1з) значение р теоретически не-ограничено.

Ясно, что ортогональные преобразования подтипа

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

В случае ортогональных преобразований II типа аналогичным образом выделяются три подтипа:

Щ) - в нем матрица и представима в форме (25) с конечным числом сомножителей, без какого-либо дополнительного структурирования;

112) - в нем матрица И имеет структуру

р

И = П Ир , (28)

р=1

и каждая из матриц Ир есть конечное произведение элементарных матриц плоского вращения; ясно, что числа сомножителей в представлениях матриц Ир могут быть различными;

113) - в нем матрица И представима по (28) с Р = то, т.е. ортогональные преобразования фактически реализуют некоторый бесконечный процесс.

В случае ортогональных преобразований III типа аналогичным образом выделяются следующие подтипы (всего их 4):

III1) - в нем вначале осуществляется ортогональное преобразование с матрицей V (А = УТА), а затем - ортогональное преобразование с матрицей И

(А = АИ); при этом подтипы как первого, так и второго преобразований могут быть различными;

Ш2) - в нем вначале осуществляется ортогональ-

V

ное преобразование с матрицей И (А = АИ), а затем осуществляется преобразование с матрицей V

◦ Т V

(А = Vх А); при этом подтипы как первого, так и второго преобразований могут быть различными; Ш3) - в нем принимается

рр V= П V?, И= П Ир , (29)

р=1 р=1

т.е. используется последовательность преобразований

А = А о ,

О гр О

Ар = Vp Ар-1Ир , (30)

р = 1, 2,..., Р;

при этом каждая из матриц V? и Ир представима в виде произведений элементарных матриц плоского вращения, содержащих конечное (но, возможно, различное) число сомножителей в виде элементарных матриц плоского вращения; число сомножителей может быть различным как при разных значениях р, так и при одном и том же р - у матриц V? и Ир;

Ш4) - в этом подтипе остаются справедливыми формулы (29) и (30), но только при Р = то, т.е. фактически речь идет об итерационных процессах, порождаемых последовательностями ортогональных преобразований исходной системы (4); ясно, что и в этом случае матрицы V? и Ир представимы произведениями из конечного числа сомножителей в форме элементарных матриц плоского вращения и при этом число сомножителей при разных р (или у V? и Ир при одном и том же р) может быть различным.

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

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

А = Ао,

А? = VX Ар_1Ир , (31)

р = 1, 2,...

и соответствующие ортогональные преобразования независимых переменных:

х = х(0),

х(^ = Ир х(Р-1) , р = 1, 2,... , системы же линейных уравнений брать в форме:

А^) = VT /^ = $) ,

р = 1, 2,... , (33)

где

/(1) = /,

/^ = /^-1) - А„_1Х <*-1) , (34)

р = 2, 3,... ,

где ж^-1) суть приближенные решения систем (33) -с индексом (р — 1). При этом приближенное решение исходной системы (4) берется в форме

р

жр = Дх(?) , (35)

?=1

где

Дх(?) = ж^ , х^ = И1 (... (ир_1 (ирх(?)...) , (36)

р =1,2,...

и суть приближенное решение системы (33).

Таким образом, здесь налицо расширение области использования ортогональных преобразований, а именно - как основы для итерационных процессов нового типа.

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

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

Класс 1. Сюда относятся ортогональные преобразования, у которых параметры матриц элементарных плоских вращений (произведения которых опре-

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

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

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

Сразу же сделаем следующее принципиальное замечание: речь идет о классификации матриц V? и И? (в соотношениях (25), (30), (31) и т.д.); возможны (и фактически весьма важны) такие ортогональные преобразования систем линейных уравнений, в которых указанные ортогональные матрицы при различных р принадлежат к разным классам; тем более это справедливо для более общих ортогональных преобразований, порождающих итерационные процессы, описанные в конце предыдущего § 2 данного раздела. Ясно, что:

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

б) Наиболее широким (что может показаться неискушенному читателю странным) является класс 1; сюда относятся все ортогональные преобразования, как классические (традиционные для вычислительной линейной алгебры, [см. Воеводин, 1977; Воеводин, Кузнецов, 1984; Уилкинсон, Райнш, 1976; Фаддеев, Фаддеева, 1963]), так и новые, обеспечивающие приведения матриц к специальным формам:

- диагональной;

- ленточной (в частности, трехдиагональной и двухдиагональной);

- хессенберговой (верхней или нижней);

- треугольной (верхней или нижней);

- стреловидной (как правило нижней).

При этом классические ортогональные преобразования указанного вида образуют подкласс 1 1 ), новые -подклассы 12) и 13); кроме того, имеется еще три подкласса: 14), включающий ортогональные преобра-

зования последовательного выметания столбцов или строк матриц; 15), включающий ортогональные преобразования столбцов и строк матриц; 16), включающий ортогональные преобразования столбцов, обеспечивающих в итоге минимальное значение отноше-

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

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

31) - ортогональных преобразований, обеспечивающих некоррелированность заданного подмножества столбцов итоговой матрицы (либо заданного подмножества парциальных столбцов итоговой матрицы) с вектором правой части системы (с заданным парциальным столбцом вектора правой части системы);

32) - ортогональных преобразований, обеспечивающих максимально возможную коррелированность некоторого столбца итоговой матрицы (либо некоторой группы столбцов итоговой матрицы) с вектором правой части системы;

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

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

Ниже приводятся более подробные (хотя по необходимости и краткие) характеристики ортогональных преобразований из выделенных подклассов. Читателю настоятельно рекомендуется перед изучением последующих параграфов раздела освежить свои знания по вычислительной линейной алгебре - в части, касающейся классических ортогональных преобразований, например, по книгам [Воеводин, 1977; Воеводин, Кузнецов, 1984; Уилкинсон, Райнш, 1976; Фаддеев, Фаддеева, 1963].

§ 4. Напомним сначала сущность классических ортогональных преобразований, основанных на использовании последовательностей элементарных ортогональных преобразований с матрицами элементарных плоских вращений, с помощью которых матрица приводится к специальному виду (диагональному, треугольному, трехдиагональному и т.д.). Эти преобразования могут принадлежать ко всем трем типам (I, II, III), указанным в предыдущем параграфе:

1) ортогональные преобразования, реализующие

приведения к треугольным или хессенберговым формам, могут принадлежать к первому или второму типам;

2) ортогональные преобразования, реализующие приведения к диагональной, ленточной (например, трехдиагональной или двухдиагональной) и стреловидной формам, всегда принадлежат к типу III.

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

Пусть, например, матрицей элементарного плоского вращения преобразуются две вектор-строки матрицы А, см. соотношения (11)-(12) и (19), и пусть преобразование осуществляется под условием:

а

= 0, 1 < к < М

(37)

в этом случае параметры преобразования определяются соотношениями:

008 р =

81П р =

аг,к

аа + аЬ

(38)

а7,к

а’,й + У

Аналогичным образом, если матрицей элементарного плоского вращения преобразуются два вектор-столбца матрицы А, см. соотношения (13)—(14) и (20), и если преобразование осуществляется под условием

а ¡у =0, 1 < к < N ,

(39)

то в этом случае параметры преобразования определяются соотношениями:

008 р

81П р =

акг

¡у

(40)

аку

ак + аку

Элементарные ортогональные преобразования строк

- по соотношениям (37)-(38) - и столбцов - по соотношениям (39)-(40) - в литературе по проблемам вычислительной линейной алгебры именуются гивенсо-выми вращениями.

Ясно, каким именно образом в классических ортогональных преобразованиях формируются последовательности элементарных преобразований - из ги-венсовых вращений:

2

1) если матрица A с N>M приводится к верхней треугольной форме, то последовательно аннулируются элементы 1-го, 2-го, 3-го и т.д. подстолбцов (т.е. в позициях (k, i), k > i, i = 1, 2,... , M-1);

2) если матрица A с N>M приводится к верхней хессенберговой форме, то также последовательно аннулируются элементы в 1-ом, 2-ом, 3-ем, и т.д. под-столбцах - в позициях (k, i), k > i + 1, i = 1, 2,... , M-2;

3) если матрица A с N<M приводится к нижней треугольной форме, то последовательно аннулируются элементы 1-ой, 2-ой, 3-ей и т.д. строк, находящиеся в позициях (i, k), k > i, i = 1, 2,... , M-1;

4) если матрица A с N<M приводится к нижней хессенберговой форме, то последовательно аннулируются элементы 1-ой, 2-ой, 3-ей и т.д. строк, находящиеся в позициях (i, k), k > i +1, i = 1, 2,... , M-2.

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

Прежде всего укажем, что в случае матриц A с N>M приведение к диагональной форме (с неотрицательными элементами - сингулярными числами <7j = <Tj(A) - на диагонали), т.е. нахождение сингулярного разложения, осуществляется с помощью следующей последовательности ортогональных преобразований:

Ri = VTA, Li = Riüi ,

R-2 = VTLi, L2 = R2U2 , ......................... (41)

Rk Vk Lk — Ъ Lk Rkük ,

k = 1, 2,... ,

при этом в пределе

lim Lk = lim Rk = Sa , (42)

k—k—

где Sa - диагональная матрица с элементами а* = ffj(A).

В (41)-(42) Rk суть верхние треугольные матрицы, Lk - нижние треугольные матрицы.

В случае N<M процесс приведения к диагональной форме аналогичен, иной является только последовательность преобразований:

Li = Aüi, Ri =VT Li ,

L2 = Riü2, R2 = V2T L2 ,

......................... (43)

Lk Rk — i ük, Rk Vk Lk ,

k = 1, 2,...

и при этом в пределе снова справедливы соотношения

(42).

Нет необходимости подробно пояснять тот факт, что в (41) и (43) Ук и Ик суть ортогональные по столбцам матрицы, представимые в форме произведений элементарных матриц плоского вращения - реализующих гивенсовы вращения.

Приведение к ленточной форме осуществляется по близким алгоритмам; здесь следует выделять три типа ленточных матриц:

1) верхние (для которых все элементы равны нулю при г — у > 0 и г — у < —т);

2) нижние (для которых все элементы равны нулю при г — у < 0 и г — у >т);

3) двусторонние (для которых равны нулю все элементы при г — у < —т1 и г — у > т2).

Приведение к верхней ленточной форме при М>М осуществляется по схеме (ширина ленты равна т + 1):

А1 = УТ А, В1 = А1и1 ,

А2 = УТ В1, В2 = А2И2 ,

......................... (44)

А к = УТ В* _1 , В к = Ак и* , к = 1, 2,..., М — (т + 2) ,

и далее

АМ — (т + 1) = УМ — (т +1)ВМ-(т+2) ,

АЯ = уТ АЯ— 1, (45)

д = М — (т + 2), М — (т + 1),..., М .

При этом в (44) преобразование с ортогональной матрицей У1 аннулирует элементы 1-го подстолбца, преобразование с ортогональной матрицей И1 аннулирует элементы 1-ой строки, начиная с (т + 2)-го, преобразование с ортогональной матрицей У2 аннулирует элементы 2-го подстолбца, преобразование с ортогональной матрицей и2 аннулирует элементы

2-ой строки, начиная с (т + 3)-го, и т.д. После того, как матрица Ик при к = М — (т + 2) аннулирует всего один элемент (М — (т + 2))-ой строки, начинает осуществляться последовательное аннулирование только элементов в подстолбцах. Ясно, что аннулирование всех элементов осуществляется с помощью соответствующих гивенсовых вращений (с помощью элементарных преобразований, осуществляемых матрицами элементарных плоских вращений).

Ясно, что приведение матриц с N<М к нижней ленточной форме (с шириной ленты т + 1) осуществляется с помощью аналогичной последовательности ортогональных преобразований:

аа

А1 =АИ1, В1 =УТ А1 ,

V

А 2

БіИ

V

Б 2

ут А 2

матрица вида

и далее

V V гр

Ак = Бк-і ик, Бк = Ук Ак

к = 1, 2,...,Ж - (т + 2) ,

V V

Аы — (т+1) Аы — (т+2)иы — (т+1)

(46)

Я2 = N

.М

0

— к

Т

N - М

I

Т

М-к

I

Т

к

I

(49)

М к к

V V Ад = Ад — іи^

д = N — (т + 2),..., N ;

(47)

суть последовательности преобразований должна быть читателю совершенно ясна: с помощью матриц ик, к = 1, 2,..., N, последовательно аннулируются элементы строк (1-ой, 2-ой, и т.д.), а с помощью матриц Ук, к = 1, 2,..., N — (т + 2), аннулируются элементы в подстолбцах - таким образом, чтобы получалась нижняя ленточная матрица с шириной ленты (т +1).

Наконец, алгоритм приведения матрицы А к ленточной форме с шириной ленты (т-1 + т2 + 1), как в случае М>М, так и в случае М<М, должен быть читателю после всего сказанного совершенно ясен и поэтому мы предоставляем ему возможность выписать все формулы самостоятельно.

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

Первая форма (которую обозначим Б1) используется в случаях, когда у матрицы А М<М; это блочная матрица

<— N — к —> *— к ^

Б1 = N

0 —к

N — к

I

(48)

к

I

М

у которой есть диагональная матрица размера

(М-к)х(М-к), блок в левом верхнем углу размера (М-М)х(М-к) - нулевой, остальные четыре блока -с ненулевыми элементами.

Алгоритмы преобразования матриц А к стреловидным формам достаточно очевидны. Например, в случае М<М сначала осуществляется ортогональное преобразование

Б = Аи0

которое преобразует матрицу А к виду

(50)

Б

Т

0 И N —

I

Т

к

I

^ (М — N) ^ ^ N — к ^ ^ к ^

у которой есть диагональная матрица размера

(М-к)х(М-к), блок в левом верхнем углу размера ^-к)х(М-Ы) - нулевой, остальные четыре блока -с ненулевыми элементами.

Вторая форма (которую обозначим Я2) используется в случаях, когда у матрицы А М>М; это блочная

(51)

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

В случае же М>М сначала осуществляется ортогональное преобразование

С = УкА

(52)

362

СТРАХОВ, СТРАХОВ: АППРОКСИМАЦИОННЫй ПОДХОД К РЕШЕНИЮ ЗАДАЧ II

которое преобразует матрицу A к виду вляется в форме:

C =

0

L

Í

N — M

I

Í

M - к I Í к I

(53)

М-к

к

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

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

а) диагонализации симметричных положительно полуопределенных матриц;

б) приведения симметричных положительно полу-определенных матриц к ленточной форме (с шириной ленты 2т + 1);

в) приведения симметричных положительно полу-определенных матриц к стреловидной форме.

Заметим, прежде всего, что если требуется (пхп)-матрицу привести к стреловидной форме Б,

k 1 S сэ

Í

n — к I

(54)

(n — к)

к

где Вп-к есть диагональная матрица размера (п-к)х(п-к), то просто сначала матрица А предста-

A=

B

Í

n — к I Í к I

(n — к) ^ ^ к ^

и далее применяется описываемый ниже алгоритм диагонализации к матрице В; при этом два блока, смежных к В, также изменяются.

Далее очевидно, что приведение матрицы А = Ат > 0 к ленточной форме (ширина ленты 2т + 1) осуществляется с помощью последовательности ортогональных преобразований

A = Ао ,

Ak = VT Ak-iVk ,

к = 1, 2,.. ., n — (m + 1)

(56)

при этом каждое из ортогональных преобразований (с матрицами Ук) обеспечивает аннулирование элементов в к-ом столбце и к-ой строке, индексы которых не принадлежат тем, которые определяют ленту. Ясно, что У к представляются в форме произведений элементарных матриц плоского вращения. Ясно, что для приложений особенно важен случай т = 1, 2т + 1=3, т.е. А = Ат > 0 ортогональным преобразованием III типа приводится к трехдиагональной форме.

Рассмотрим теперь в некотором смысле самый важный классический алгоритм для матриц - алгоритм приведения к диагональной форме.

Этот алгоритм таков:

A=A

о

TT

Ak к = 1, 2

ik Зк (Pk) Ak-1Tikjk (Pk ) , . . , 1 < ik < jk < n .

(57)

Здесь, как всегда, Т^(р>к), 1 < 1к < ]к < п, суть элементарные матрицы плоского вращения, параметры которых выбираются из условия [см. Воеводин, 1977; Воеводин, Кузнецов, 1984; Уилкинсон, Райнш, 1976; Фаддеев, Фаддеева, 1963]:

ik jk

jk ik

(58)

здесь apq суть элементы матриц Ak. Ортогональные преобразования, фигурирующие в (57), параметры

- (cos ifk, sinpk) - которых выбираются из условий

(58), именуются якобиевыми вращениями. Важным

<----

і----

S

компонентом алгоритма является генерация последовательности (¿к, ук). Здесь существует целый ряд возможностей, [см. Воеводин, 1977; Воеводин, Кузнецов, 1984; Уилкинсон, Райнш, 1976; Фаддеев, Фад-деева, 1963]. Важнейшая из них - последовательное аннулирование максимальных по модулю элементов в каждой строке.

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

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

Замечание 1. Пусть произвольная матрица А х) преобразуется в матрицу А с помощью (12); тогда для элементов преобразуемых ¿-ой и у-ой строк справедливы соотношения (вытекающие из (19)):

(59)

Замечание 2. Пусть произвольная матрица А 2)

V

преобразуется в матрицу А с помощью (14); тогда для элементов преобразуемых ¿-го и у-го столбцов справедливы соотношения (вытекающие из (20)):

V X 2 /V \ 2

ak,i) + ( ak,j )

k = 1, 2,

(60)

Замечание 3. Классические ортогональные преобразования по виду носителей (множеств пар индексов (р, д), которые присваиваются ненулевым элементам) тех специальных форм, к которым приводятся матрицы, целесообразно классифицировать на два вида:

первый вид - множество пар индексов ^ определяющих носитель специальной формы, не разделяет множество пар индексов из Ш на два непустых подмножества; иначе - Ш - связное множество;

второй вид - множество пар индексов ^ определя-

ющих носитель специальной формы, разделяет множество пар индексов из CJ на подмножества; иначе

- CJ не является связным множеством.

Замечание 4. Носителем строки (столбца) заданной матрицы (полный носитель которой есть множество пар индексов (p, q), 1 < p < N, 1 < q < M) называется множество пар индексов вида (p, k), p = const, k = 1, 2,..., M (соответственно - множество пар индексов (k, q), q = const, k = 1, 2,..., N). Каждая пара строк (столбцов) в матрице с заданным полным носителем, которая преобразуется к специальной форме с носителем J (для ненулевых элементов), принадлежит к одному из трех типов:

тип 1: обе строки (оба столбца) имеют носители с парами индексов, принадлежащих J;

тип 2 - из двух строк (столбцов) одна из них (один из них) имеет носитель с парами индексов, принадлежащими J, соответственно вторая (второй) имеет носитель с парами индексов, принадлежащими CJ;

тип 3 - из двух строк (столбцов) ни одна из них (ни один из них) не имеет носителя, содержащего пары индексов, принадлежащих J (оба носителя принадлежат CJ).

Одновременно ясно, что:

1) в случае матриц размера (nxn) имеются пары строк (столбцов) только типа 1;

2) в случае матриц размера NxM, N>M, пары строк могут принадлежать к любому из трех типов, а пары столбцов принадлежат только типу 1;

3) в случае матриц размера NxM, N<M, пары столбцов могут принадлежать к любому из трех типов, а пары строк принадлежат только типу 1.

Замечание 5. Пусть задана пара строк (столбцов) с носителями 1-го или 2-го типов, т.е. пара строк (столбцов), у которых имеется множество пар индексов, принадлежащих J. При конструировании орто-тональных преобразований приведения заданной матрицы к заданному специальному виду (с носителем ненулевых элементов J) последние строятся в виде произведений элементарных матриц плоского вращения, определяемых парами строк (столбцов), у которых параметры вращений определяются из условий:

а)

?(0

в случае пар строк;

б)

(i)

2 + a j) 2 max

E(J) E(J) V

2 V 2

+ a(j) = max

E(J) E(J) Ф

(61)

(62)

1) Ясно, что в общем случае это не исходная матрица, а текущая — получаемая в процессе некоторого общего ортогонального преобразования.

2) См. предыдущее примечание.

- в случае пар столбцов.

При этом в (61) использованы соотношения (12) и (19), а в (62) использованы обозначения (14) и (20). Ясно, что при этом под А понимается не исходная,

^ V

а текущая матрица, соответственно под А и А - матрицы, получаемые после преобразования. При этом

2

2

364

СТРАХОВ, СТРАХОВ: АППРОКСИМАЦИОННЫй ПОДХОД К РЕШЕНИЮ ЗАДАЧ II

очевидно, что в силу “законов сохранения” (59) и (60) экстремальные задачи (61) и (62) должны быть приведены к следующему виду:

V"' -2

Е!.’

kEJi к£Л^

в случае пар строк, и к виду

V

Еа’

к£3*

,«+Е

к£Л*

?к,)

тах

Р

тах

Ф

(63)

(63’

В рассматриваемом случае 3), если

А = ТТ (р)ЛТ) (р) 1 < * < 3 < п ,

и матрица приводится к диагональной форме, целесообразно использовать постановку

22 агг - а))

тах

р

(65)

- в случае пар столбцов. При этом множества индексов ^, X,, X*, X* определены следующим образом:

1) Хг есть подмножество индексов к - таких, что (*, к) £ X, а (3, к) £ СХ;

2) X, есть подмножество индексов к - таких, что (3, к) £ X, а (*, к) £ СХ;

3) X* есть подмножество индексов к - таких, что (к, *) £ X, а (к,3) £ СХ;

4) X* есть подмножество индексов к - таких, что (к,3) £ X, а (к,*) £ Ш.

Замечание 6. Если пара строк (столбцов) из заданного носителя полной матрицы принадлежит типу 3, то при конструировании ортогонального преобразования исходной матрицы эта пара строк (столбцов) для построения элементарного ортогонального преобразования (с матрицей элементарного плоского вращения) не используется. При построении элементарных ортогональных преобразований целесообразно (хотя и не обязательно) применять такую стратегию:

а) вначале используются пары 2-го типа, если таковые имеются (они всегда имеются только в случае прямоугольных матриц);

б) далее используются пары 1-го типа (после того, как строки (столбцы) принадлежащие Ш, стали иметь нулевую эвклидову норму).

Замечание 7. Особым является случай квадратных симметричных положительно полуопределенных матриц. В этом случае, как и в рамках классических ортогональных преобразований, должны использоваться “симметричные” элементарные ортогональные преобразования (вида (57)); при этом очевидно, что носители X множеств пар индексов для матриц специальной формы всегда обладают симметрией (если (р, д) £ X, то и (д,р) £ X), поэтому в X всегда входят пары индексов вида (р,р), т.е. соответствующие диагональным элементам. Ясно, что в рассматриваемом случае матриц Л = Лт > 0 возможны случаи:

1) X = {(р,р)}^, т.е. матрица специального вида -диагональная;

2) X = {(р,р)}]_, т.е. матрица специального вида, хотя и содержит диагональ, но также имеет и другие элементы - в симметричных комбинациях.

Если же решается задача приведения матрицы к более общей, чем диагональная, специальной форме, то тогда в преобразовании (64) параметры матриц элементарного плоского вращения должны определяться из условия

4 - )■)+Еа2

г,к

а2

к,)

тах

Р

),к

(66)

здесь р > 0 - заданная константа, а суммирование осуществляется по таким парам индексов (*, к), (3, к),

(к, *), (к, 3), где к = *, к = 3, которые удовлетворяют следующим условиям:

1) если (*, к) £ X, то (3, к) £ Ш;

2) если (3, к) £ X, то (*,к) £ Ш;

3) если (к,*) £ X, то (к,3) £ ^;

4) если (к,3) £ X, то (к,*) £ Ш.

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

В случае прямоугольных матриц Л все процедуры приведения к специальному виду с помощью последовательности элементарных ортогональных преобразований, кроме процедуры приведения к стреловидным формам, целесообразно реализовывать в два этапа: первый этап - выметание матрицы на квадратный носитель (на носитель X, включающий пары индексов (р, д) вида: 1 < р < тт(Ж, М), 1 < д < тт(Ж, М)).

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

3) Ясно, что под А здесь понимается уже не исходная, текущая матрица.

2

р

2

к

а

элементарные ортогональные преобразования первого

где

Л = Ло ,

ТТ„„ (Рк )Лк-1, <67>

1 < *к < м, М +1 < 3к < N (68)

и параметры ортогональных преобразований находятся из условий:

,(гк)

= тах

Е рк

(69)

где а(гк) - суть *к-ая вектор-строка матрицы Лк. Эта последовательность ортогональных преобразований, реализующих “выметание” Л из носителя Ш в носитель X, осуществляется до тех пор, пока не выполнится (при некотором к) неравенство

Лк

2

Е(ОЛ)

Лк

2

Е(Л)

< е2

(70)

где е2 > 0 - заданная малая постоянная.

Далее ясно, что в случае прямоугольных матриц с М<М на первом этапе используются только элементарные ортогональные преобразования второго типа:

Л = Ло ,

Лк = Лк—1 Тгь,)ь (фк ) , (71)

где

1 < *к < N, N +1 < 3к < М (72)

и параметры ортогональных преобразований находятся из условия

1(^)

= тах

Е Фк

(73)

где а(гк) - *к-ый столбец Лк. Преобразования ведутся до того значения к, при котором выполнится неравенство (70).

Что касается второго этапа в рассматриваемых задачах, т.е. приведения полученных квадратных матриц к специальному виду, то здесь должны использоваться, в циклической последовательности, как преобразования строк, так и преобразования столбцов матрицы, с определением параметров ортогональных преобразований по формулам, приведенным в Замечании 5 из настоящего параграфа. Необходимо использовать последовательность ортогональных преобразований вида

Л = Л,

(о)

'Т т 'Т V

Л(1) = VI Л(о) = Л(о)

Л(1) = Л(о) и1 ,

- т V

Л(к) = V, Л(к-1)

V Л

Л(к) = Л(к) ик ,

в которой принимается

1 -1

(74)

Vk = П vkp)

И

Р=1

1 -1 П[ и

9=1

(9)

к

(75)

(76)

и каждая из матриц V(p) и И(д) формируется в форме произведения элементарных матриц плоского вращения:

1—р

Vkp) = П Тр,)

в = 1 1—р

рРк2

ик9) = П т, (фРк«0

(77)

(78)

где

3з = р + 1, р + 2,..., п . (79)

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

Наконец, очевидно, что если прямоугольная матрица Л преобразуется к матрице стреловидной формы - вида (48) при М>М и вида (51) при М<М, то тогда первый этап (преобразования прямоугольной матрицы к квадратной) исключается, и приведение осуществляется с помощью циклической последовательности преобразований (74), в которой принимается: N—1

V = П ^

Р=1

м—1

И = П и“ .

9=1

и каждая из ортогональных матриц vkp) и представляет из себя произведение элементарных матриц

(р)

к

(80)

типа:

к

2

2

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

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

Лг

2

E(CJ)

Ar

2

E(J)

(81)

/ О ,2 /О ,2

(ай) - (jv) = max ;

Р, ф

в случае 2) j' = j :

, ◦ \ 2 ,О \ 2

(ajjJ - = max ;

р,ф

в случае 3) j' = i :

/О \2 /О \2

(ай) - (о^,) = max ;

p, ф

в случае 4) j = i' :

/О \2 /О \2

Kj) - Kj') =max .

^, ф

(84)

(85)

(86)

Можно показать [см. Страхов, 1997е], что при использовании преобразований вращения Якоби-Страхова выполняются следующие соотношения:

- в случаях 1) и 2)

в случаях 3) и 4)

0

(88)

(89)

где п2 > 0 - априорно заданное малое число.

В заключение настоящего параграфа остановимся еще на одном специальном виде элементарных ортогональных преобразований III типа, которые именуются преобразованиями вращения Якоби-Страхова [см. Страхов, 1997е; Страхов, Тетерин, 1993] и которые имеют вид

О т

А = ТТ (р)ЛТ^ (ф) , (82)

при этом обязательно имеет место одно из следующих соотношений:

1) *' = *; 2) / = 3; 3) 3' = *; 4) 3 = г' . (83)

Соответствующие этим соотношениям элементарные ортогональные преобразования именуются преобразованиями вращения Якоби-Страхова 1-го, 2-го,

3-го и 4-го видов.

Параметры преобразований (82) с указанными соотношениями (83) для индексов определяются из следующих условий:

- в случае 1) *' = * :

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

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

Ортогональные преобразования подкласса 14). В рамках этого подкласса выделяются два ортогональных преобразования. Первое - типа I; оно используется в основном в случае матриц с М>М:

Л = Л о ,

Л к = Vт Л к—1 , к = 1, 2,... ; при этом матрицы Vk имеют структуру

N —1

(90)

Vk = П Vf” ,

p=1

(91)

N

vkР) = П Тр,г (рРк^ . (92)

г=р+1

При этом параметры матриц элементарного плоского вращения, фигурирующих в (92), выбираются

из условия максимума квадрата эвклидовой нормы Ортогональные преобразования подкласса !5). В

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

Если в (91) положить к = то, то придем к сингулярному разложению матрицы А,

гр ◦

Vх А = ◦ А W ,

Еа

£ А т М 4

т

0 N - М

4

Аі = иі ,

Еа =

V

Еа

v ^

где Еа = ^іадаі(А), и = и Ф есть ортогональ-

к=1

них используется следующая последовательность преобразований матрицы А:

А = Ао ,

V

А і = ^ Ао, Аі = А іИі

V

А 2 = Vх Аі, А2 = А 2И2

(98)

(93)

где аа = ¿*адаг(Л), V = П Vk, и W есть ортого-

к=1

нальная по столбцам (МхМ)-матрица.

Второе преобразование в подклассе ^) - типа II; оно используется в основном в случае матриц с N<М:

V

А = Ао ,

(94)

(95)

V V

Лк = Лк—1 Ик , к = 1, 2,..., ; при этом матрицы Ик имеют структуру

м—1

И = П и“ ,

?=1 м

икР) = П Тр,г (фРк^ . (96)

г=р+1

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

Если в (94) положить к = то, то снова придем к сингулярному разложению матрицы Л:

ЛИ = ФЕа ,

Л к ^т Л—1, Л = Л И , к = 1, 2,... .

Здесь Vk и Ик, как обычно, суть ортогональные по столбцам матрицы, размеров ^х^- и (МхМ)-, соответственно. При этом структура матриц Vk и Ик может быть различной, различными могут быть и способы назначения параметров тех элементарных матриц плоского вращения, произведениями которых представимы матрицы Vk и Ик.

В наиболее простом варианте, который используется для нахождения сингулярного разложения матрицы Л (в общем случае прямоугольной, при этом может быть как N>M, так и N<M, а также N=M) матрицы Vk имеют структуру (91)-(92), при этом параметры матриц элементарного плоского вращения в (92) снова определяются из условия максимума квадрата эвклидовой нормы р-ой вектор-строки, получаемой после преобразования, а матрицы Ик имеют структуру (95)-(96), при этом параметры матриц элементарного плоского вращения в (96) снова определяются из условия максимума квадрата эвклидовой нормы р-го вектора-столбца, получаемого после преобразования.

Как сказано выше, процесс (98) снова приводит (при к ^ то) к сингулярному разложению матрицы Л.

Ясно, что в рамках подкласса ^) возможны и другие, по структуре матриц Vk и Ик, ортогональные преобразования. Здесь особенно интересен случай матриц Л с N=M (т.е. квадратных), и особенно случай с Л = Лт > 0.

В этих случаях:

а) матрицы Vk и Ик выбираются непосредственно в форме произведений элементарных матриц плоского вращения, чаще всего вида

(97)

М

(к)

г=к+і

М

ик8) = П Тк,г (^)

ная по столбцам (МхМ)-матрица.

г=к + і

о

где

s = const, s = 1, 2, k = 1,2,..., M- 1

(100)

б) в случае матриц Л = Лт > 0 всегда в (99)-(100) принимается фГк) = рГк) (симметричные преобразования). При этом принципы задания параметров могут быть различными, но самый простой полностью аналогичен описанному выше.

Ортогональные преобразования подкласса ^). Здесь имеется всего два преобразования, оба они принадлежат типу II, различаются же лишь способом задания значений параметров матриц элементарного плоского вращения.

Итак, в преобразованиях подкласса ^) используется последовательность ортогональных преобразований (90), в которых ортогональные по столбцам матрицы берутся в форме:

M

p=2

k = 1, 2,

(101)

Ограничимся рассмотрением случая, когда ортогональные преобразования обеспечивают выметание первого подстолбца (возможны и варианты выметания первого столбца без нижней компоненты, либо последнего столбца без нижней или верхней компоненты и т.д.). В этом случае параметры вращения (cos и sin pPk)) могут выбираться одним из двух способов:

первый способ - из условия

1(k>p)

г(1)

- (2 + afc,p) (a11’p))

min

pPk)

(102)

^k,p ^ 0 ; второй способ - из условия ,(k,p)

,(k>P)

(1)

11

,(k,P)

11

min

pPk) ^p

(103)

здесь a

(k,p)

(1)

первый столбец текущей матрицы, получаемый после преобразования с матрицей элементарного плоского вращения Т1,р(рРк)), см. (101),

а(к,р) а11

первая (верхняя) компонента этого столбца.

§ 7. Теперь коротко об ортогональных преобразованиях, принадлежащих классу 2. Фактически здесь имеется всего два одинаковых по смыслу преобразования; в них фигурирует либо исходная система (4), либо система, получающаяся из (4) некоторым преобразованием (см. раздел 3, § 4 данной статьи). В

любом случае пусть имеется система

Лх = 7г; эта система преобразованием

B = VT Л, Vт г = рг приводится к одной из двух форм:

Вх = рг ;

1) Рг = С1в, 2) С’а ,

где

(105)

т

1 1

1

т

0 N1 - 1

1

т

0 N' - 1

1

т

1 1

1

(106)

(107)

и в общем случае N7 = N. Ясно, что в случае 1)

N'

У = ПТ1,Цpp) ,

p=1

(108)

а в случае 2)

N '-1

v= П

q=1

Tiq,N '

(109)

при этом параметры *р и выбираются по принципу малости модуля аннулируемого элемента (сначала аннулируется наименьший по модулю, далее наименьший по модулю среди оставшихся, и т.д.).

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

§ 8. Итак, остается лишь описать ортогональные преобразования 3-го класса, т.е. такие, в которых параметры матриц элементарного плоского вращения выбираются как по значениям элементов (текущих) матриц, так и по значениям (текущим) компонент вектора правой части.

Ясно, что в классе 3 могут использоваться ортогональные преобразования всех типов: I, II, III.

Ортогональные преобразования подкласса 31).

Эти преобразования относятся к типу II. В них используются матрицы Ик, которые представимы в форме

к

^шах

Ик = П Тгк,) (р)к) к=1

e =

p

2

2

E

2

(к)

в которых параметры р:- определяются из условий:

(111)

где а(к) есть либо вектор-столбец, либо парциальный вектор-столбец с номером 3 матрицы, получаемой после преобразования с элементарной матрицей Тгк,)(р)к)), а 7г к есть либо соответствующий вектор правой части, либо соответствующий парциальный вектор правой части.

Ортогональные преобразования подкласса 32).

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

А = Ао ,

V V

Ак = Ак-іик ,

(112)

- при само собой понятных обозначениях.

Ортогональные преобразования подкласса 33). Преобразования этого подкласса принадлежат типу I. В них используется последовательность преобразований

Л = Л о ,

в которых

А к = ^А й-і к = 1, 2,... ,

N

Vk = П V?’,

р=і

N

V?’ = П Тр:

:=Р+і

(118)

(119)

(120)

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

Первый способ: из условия

где

Ик

„(к)

М

П Ті.: (4

:=2

(к)

(113)

и параметры р:- находятся из условий

2

,(к,:)

(і)

,(к,:)

(і)

/»

2

тах

(к) Р: )

(114)

Ик = П иГк)

где

и(к)

М

П Тг.: М'

:=д+і

(к.г)

(115)

(116)

(к.г)

при этом параметры р: находятся из условий

а(к;Л/»

-,(к,:)

V)

/й

= тах

р

(к.г)

(117)

,(:)

— 1 + г

[к.р]

7к1

Второй способ: из условия

/(к.Р)

/:.й

р

тт

(к.р)

(121)

/(к.Р)

/:.й

Здесь а(к,)) - вектор-столбец (текущей) матрицы, получаемый после преобразования с элементарной матрицей Т1,) (р(к)), а 7г - вектор-столбец правой части.

Более общие преобразования рассматриваемого подкласса получаем в случае, если фиксируется набор из И. первых столбцов матрицы и матрицы Ик выбираются в форме

7(:)

[к.р]

= тах

р(к.р)

(122)

Здесь а(к) ] суть (текущая) у-ая вектор строка, полу-

(к.рЬ

чаемая после преобразования с матрицей Тр )- (р'

(к р)

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

Ортогональные преобразования подкласса 34).

В этих преобразованиях фиксируется номер одного столбца, обычно первого (возможен вариант, когда фиксируется последний столбец матрицы, но он рассматриваться не будет) текущей матрицы, и рассматривается последовательность преобразований II типа:

V

Л = Ло ,

V V

Аі = Аоиі

V V

Ак = Ак-іик к = 1, 2,... ,

и

2

2

2

2

Е

Е

2

2

2

Е

Е

в которых

M

Ufc = nTl-j j=2

(fc)

(124)

и параметры элементарных матриц плоского вращения выбираются таким образом, чтобы в конечном

V

итоге (у предельной матрицы Ато, но фактически, с

V

некоторой допустимой погрешностью, у матрицы А с некоторым не слишком большим k) выполнялись условия:

(то) (то)

«(1) ,а(.)

0, i = 2, 3, ...,M

(125)

(а(о),7^ =0, * = 2, 3,...,М. (126)

Добиться выполнения предельных соотношений (125)—(126) можно, выбирая параметры матриц Т1,) (р)к)) из условий:

р(к)Л(к) а(кЛ2 . (к)/ (к) 7 \2= • (127)

р) (va(l), а()^ +?) (Д)^ П) = П1111 , (127)

Р.

где

> 0

(fc)

(128)

суть априорно задаваемые веса.

Ясно, что в (125)-(127) а(к)), а) суть *-ый и у-ый векторы-столбцы (текущей) матрицы, получаемой после преобразования с матрицей элементарного плоского вращения Т1,)- (р)к)).

Существуют и другие способы определения параметров матриц элементарного плоского вращения, произведения которых определяют ортогональные матрицы ик в (123), основанные на прямом использовании условий (125) и (126) в некоторой последовательности, но они, по-видимому, менее эффективны по сравнению с описанным.

§ 9. В заключение раздела приведем некоторые общие соображения по проблеме устойчивости ортогональных преобразований. Эта проблема имеет два аспекта:

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

2) аспект устойчивости в проблеме нахождения параметров матриц элементарных плоских вращений.

Ясно, что на точность окончательного результата влияют оба типа погрешностей при вычислениях.

Что касается первого аспекта, то он изучен достаточно полно в случае гивенсовых вращений, см., например, [Воеводин, 1977; Воеводин, Кузнецов, 1984]; при этом большое значение имеет порядок выполнения операций - таким образом, чтобы уменьшить влияние погрешностей округления, связанных с неточностью определения параметров элементарных преобразований. Здесь следует подчеркнуть, что при реализации гивенсовых вращений по мере роста числа произведенных вращений относительная точность определения элементов матрицы, которые в дальнейшем должны быть аннулированы, заметно падает; это определяет накопление погрешности на множестве тех элементов, которые аннулируются в заключительной фазе преобразования. Именно в связи с указанным эффектом накопления погрешностей классические ортогональные преобразования “отказывают” в случае плотно заполненных матриц большой размерности, примеры чего уже имеются; они будут приведены в четвертой части работы.

Что касается новых ортогональных преобразований, то с позиций влияния погрешностей в нахождении параметров элементарных вращений, следует выделить два рода преобразований:

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

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

Очевидно, что влияние погрешностей округления при нахождении параметров элементарных ортогональных преобразований в случае преобразований 1го рода существенно выше, чем в случае преобразований 2-го рода; поэтому использование преобразований 2-го рода существенно более предпочтительно. Указанное соображение безусловно должен иметь в виду читатель при ознакомлении с содержанием двух последующих разделов настоящей статьи.

2. Новые конструктивные идеи нахождения устойчивых приближенных решений систем линейных алгебраических уравнений, являющиеся наиболее перспективными

§ 1. Как сказано в первой части работы [Страхов, Страхов, 1999], необходимо классифицировать не конкретные методы регуляризации систем линей-

1

ных алгебраических уравнений, а те конструктивные идеи, которые используются в конкретных методах. Три различных классификации 20 основных конструктивных идей были приведены в [Страхов, Страхов, 1999]. При этом к числу наиболее перспективных (по показателям точности - устойчивости -и быстродействия) были отнесены следующие конструктивные идеи4):

113) использования новых итерационных методов, основанных на априорной информации о существенной разделенности спектров полезного сигнала и помехи;

121) использования подхода, основанного на построении и решении последовательности систем малой размерности;

122) использования подхода авторегуляризации;

123) использования подхода фильтрации заданной системы линейных уравнений;

125) использования подхода регуляризации алгоритмов мультипликативных разложений квадратных матриц на треугольные сомножители;

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

135) использования итерационных процессов, основанных на ортогональных преобразованиях, изменяющих отношение эвклидовых норм парциальных векторов первого столбца.

Ниже (в §§ 2-8 данного раздела статьи) приводятся краткие описания сущности этих наиболее перспективных конструктивных идей.

§ 2. Понятие существенной разделенности спектров полезного сигнала и помехи дано в [Страхов, Страхов, 1999]. Итерационные методы, о которых пойдет речь, принадлежат к числу классических - в том смысле, что они предназначены для нахождения приближенных решений систем с симметричными положительно полуопределенными матрицами, и в них используются полиномы от матриц.

Итак, пусть дана система

Бад = рг = р + ¿р

в которой

Б = Бт > 0

Б1

(129)

(130)

4) Здесь используется номенклатура конструктивных идей, данная в [Страхов, Страхов, 1999].

ном случае может быть Б=А, если А = Ат > 0, и тогда т = х, а в других случаях: 1) Б = АтА (и тогда рг = Ат/г, т = х); 2) Б = ААт (и тогда т = г, X = Атт, рг = /г).

Итерационные процессы, о которых идет речь, определяются следующими соотношениями:

Р

(и) = рг - Бт(п)

р(и) = Р(„)(Б)рг ,

П

Р(п)(Б) = П П(к)(Б) ,

к=0

Р(0)(Б) = П(0)(Б) = Е ,

Р(и)(Б) = п(п)(Б)Р(„-1)(Б) п = 1, 2,...

(131)

(132)

(133)

(134)

(0) = 0

т(п) =д(„)(Б)рг ,

п = 1, 2,...

(135)

д(П)(Б)= (е - Р(П)(Б))Б

и) 1, 2,...

т(и) = т(и-1) + д(и)(Б)р() ?(и)(Б) = (Е - П(„)(Б^Б-1

(136)

(137)

(138)

В приведенных соотношениях Р(и)(Б), П(к)(Б), д(и)(Б), ^(к)(Б) суть алгебраические полиномы от матрицы Б, степеней Ж„, пк, N — 1, пк — 1, соответственно. Ясно, что если г, —то < г < +то, суть вещественная переменная и Р(„)(г), П(к)(г), д(п)(г), ад (г) суть соответствующие алгебраические полиномы от г, то должны выполняться следующие необходимые условия:

а) на полиномы П(к)(г) :

при этом ||Б|| = Лтах(Б) суть спектральная норма матрицы Б (ясно, что если ||Б|| > 1, то всегда можно перейти к ситуации (130), если найти (приближенно, но с избытком!) значение Лтах(Б) и перейти к матрице Б/Лтах(Б)). Ясно, что в (129)-(130) в част-

П(к)(0) = 1, к =1, 2,..., П(к)(£) < 1, 0 < £ < 1, к =1, 2,...

11т

к

0, 0 < г < 1 ;

п(к)(г)

б) на полиномы Р(и) (г) :

Р(п)(0) = 1, п =1, 2,...,

Р(и)(г)

< 1, 0 < г < 1, п = 1,2,... .

(139)

(140)

(141)

(142)

(143)

и

1

п

lim P(n)(t) =0, 0 < t — 1 .

(144)

Ясно, что (142) есть следствие (139), (143) есть следствие (140), (144) есть следствие (141) - в силу соотношений (132)-(134).

Кроме необходимых условий имеются еще и дополнительные желательные условия:

1) на полиномы П(к)(г) :

П(к) (1) =0, к = 1, 2, .. . ,

dn(k) (t)

dt

к = 1, 2,.

t=0

2) на полиномы P(n)(t) :

P(n)(1) =0, n =1, 2,...

dP (n) (t)

dt

0,

1, 2,... .

t=0

(145)

(146)

(147)

(148)

Здесь также ясно, что (147) есть следствие (145), а (148) есть следствие (146).

Кроме того, весьма желательно выполнение следующих условий на полиномы П(и)(г) :

1) если пк есть степень полинома, то все корни этого полинома располагаются на отрезке (0, 1); 2) если тк1) и ткПк) суть наименьший и наибольший корни полинома П(к)(г), то:

SUp n(fc)(t) — ^ < 1 .

rk1)<t<rk”fc)

(149)

Здесь ^к есть задаваемый параметр.

Очевидно, что при реализации всех указанных условий на полиномы будут иметь место соотношения:

(n+1)

<

(n)

1, 2,...

(150)

lim p(n) = p( oo) , (151)

n—— oo

где p( oo) - предельный вектор, который может быть и ненулевым (если Ker B не пусто).

Ясно, что в общем случае в процессе построения последовательных приближений могут находиться (хотя это и не гарантировано, но в большей части случаев это имеет место) пробные решения, для которых

¿2 • < min —

(n)

(152)

где и ¿тах - априорно заданные числа. Но в об-

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

1) использовать возможность задания полиномов

обеспечивающих различную скорость сходимости;

2) использовать конструкции нахождения так называемых “дробных приближений”.

Здесь нет возможности подробно рассматривать способы задания полиномов n(n)(t) и различные модификации описанного основного метода (“метода последовательного умножения полиномов”, предложенного В. Н. Страховым). Подробное рассмотрение метода последовательного умножения полиномов и его модификаций будет дано авторами в отдельной работе.

§ 3. Идея находить приближенные решения a;(fc), к = 1, 2,... , системы (4) как решения некоторых систем малой (последовательно повышающейся) размерности является, как это подчеркнуто в [Страхов, Страхов, 1999], исключительно важной и перспективной; эта идея предложена В. Н. Страховым, сначала в частной форме [Strakhov et al., 1995], и далее - в общей форме [Страхов, 1997d].

Данная конструктивная идея может быть использована в следующих модификациях:

1) основанная на выметании строк;

2) основанная на выметании столбцов;

3) основанная на выметании строк и столбцов;

4) основанная на приведении матриц с N>M к трапециевидной форме;

5) основанная на преобразовании систем с матрицами A = AT > 0, с выделением подсистем с преобладающей эвклидовой нормой.

Существуют и иные модификации рассматриваемой конструктивной идеи (например, основанных на выделении стреловидных подсистем с малой шириной стрелы, и т.д.), но указанные - важнейшие. Их краткое описание приводится ниже.

Модификация, основанная на выметании строк (N>M). Здесь используется последовательность ортогональных преобразований:

Ak

A = Ao ,

= VTAfc-1

(153)

k = 1, 2,...

в которой

N

Vk = П ,

к = 1, 2, . . .

(154)

p=k + 1

и параметры каждого элементарного преобразования (с матрицей Тк,р(ррк))) определяются под условием максимума квадрата эвклидовой нормы соответствующей (к-ой) строки, в которую “выметаются” остальные строки.

0

n

n

E

E

и

2

E

Системы малой размерности, формируемые по ходу преобразований, суть

Ак Акт г(к) = рг(к)

(155)

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

— (к) —(к) Л7-т (к-1)

торы рг получаются из векторов р£ = ^¿рг ,

рг0) = /г, к = 1, 2,... , отбрасыванием всех компонент, кроме к первых ( верхних), и приближенные решения ж(к) формируются по правилу:

углу). Системы (159) также решаются экономично с помощью разложения Холецкого, при этом уже выполненное разложение матрицы Бк эффективно используется при нахождении разложения Холецкого матрицы Бк+1.

Ясно, что векторы ж(к) - приближенные решения исходной системы (4) - формируются с помощью следующих соотношений:

(к) =

(к)

^(к)

(156)

(к)

(160)

Ак+1Ак+1 по-

Ясно, что матрицы систем (155) формируются по

правилу окаймления; именно, Бк+1 ‘ ‘

лучается из матрицы Бк = АкАдГ присоединением к последней еще одной строки и одного столбца, при этом элементы строки - скалярные произведения (к+1)-ой строки с предыдущими к строками - и с собой самой (элемент в правом нижнем углу). Ясно, что системы (155) экономично решаются с помощью разложения Холецкого, при этом уже выполненное разложение матрицы Бк эффективно используется при нахождении разложения Холецкого матрицы Бк+1.

Модификация, основанная на выметании столбцов (М<М). Здесь используется последовательность ортогональных преобразований

г(к) =и^И2(... (ИкХ(к))...)) .

Модификация, основанная на выметании строк и столбцов. Здесь используется последовательность ортогональных преобразований:

при этом

в которой

А = А0

А = А-1и к = 1, 2,...

N

и = П Тк,р (фРк)) р=к+1

(157)

(158)

V

А = А 0 ,

'Т V V ^

А1 = У]_ А0, А1 = А1И1 ,

-г т V V ^

А к = ук Ak-l, Ак = А к ик

к = 1, 2,... ;

N

Ук = П Тк+1,р(рРк))

р=к+2

М

ик = П Тк+1,Р(фРк)) .

р=к+2

(161)

(162)

и параметры каждого элементарного преобразования (с матрицей Тк,р(фРк))) определяются под условием максимума квадрата эвклидовой нормы соответствующего (к-го) столбца, в который “выметаются” остальные столбцы. Системы малой размерности, формируемые по ходу преобразований, суть

Определение параметров матриц элементарного вращения, фигурирующих в (162), осуществляется по тем же самым правилам, что в первой модификации (для сомножителей в представлениях матриц Ук) и во второй модификации (для сомножителей в представлениях матриц ик). В рамках рассматриваемой модификации используется последовательность систем линейных уравнений

Бкх(к) = АТ Ак х(к) = ргк) = АТ /г ,

(159)

(к)

(к)

(163)

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

Ясно, что матрицы Бк систем (159) формируются последовательно, по правилу окаймления - Бк+1 получается присоединением к Бк еще одной строки и одного столбца, при этом элементы строки - скалярные произведения (к+1)-го столбца с предыдущими столбцами и самим собой (элемент в правом нижнем

где Ак суть (к+1) х (к+1)-матрица, элементы кото-

V

рой суть элементы а(к) матриц Ак, 1 < * < к + 1,

□

1 < у < к +1, а (к+1)-векторы ф(к) имеют своими компонентами первые (верхние (к+1)) компоненты векторов ф(к) = Ут Ф(к 1), Ф(0) = /г. Ма-

□

трицы Ак+1 опять-таки формируются по принципу

0

г

окаймления матриц Ак. Нахождение решений систем (162) осуществляется с помощью разложения

□

матриц Ак в произведение двух треугольных с помощью компактной схемы Гаусса [см. Воеводин, 1977; Воеводин, Кузнецов, 1984; Уилкинсон, Райнш, 1976; Фаддеев, Фаддеева, 1963]. При этом ясно, что нахо-

□

ждение разложения матрицы Ак+1 осуществляется

□

эффективно, если известно разложение матрицы Ак.

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

) т к + 1 1

т

0 М - (к + 1)

1

з(к) = щ( и2(... (и^«)...;

Модификация, основанная на приведении матрицы к трапециевидной форме, М>М. В этой модификации используется последовательность ортогональных преобразований

А = А о ,

где

А і — ЛоИі, А і — У^А і

V „V

Ак — А к-іЩ, А к —Ук А к — 1, 2,... ,

м

и — П тк,р(фРк)), р=к

N

Ук — П Тк,р(рРк)) , р=к

б) у матриц Тк,р(рР>к)) параметры вращения выбираются из условия, что все элементы к-го подстолбца (т.е. расположенные в к-ом столбце ниже главной диагонали), последовательно аннулируются 5).

В итоге матрицы Ак, получаемые после к-го шага, имеют следующую (трапециевидную) структуру:

к

(М - к)

т

к К-к

1 Вк

т

N - к 0

1

(164)

(167)

где И.к суть верхняя треугольная матрица.

В данной модификации решается последовательность систем с верхними треугольными матрицами

Ик £(к) —

Рй

(к)

(168)

где И.к - блок матрицы Ак размера (кхк), гКк) - к-~(к) ,

вектор, рй - к-вектор, определяемый разложением

(к)

Рй ) —

-(к) Рй ) т к 1

V т

Р(к) N - к

Рй і

(о)

Рй ) — /й,

-/•(к) лтк /(к — і)

/й _ Ук /й .

После решения системы (168) находится вектор

(169)

(165)

У(к)

£(к) т к і

т

0 М - к

і

(170)

(166)

и параметры элементарных матриц плоского вращения в (166) выбираются следующим образом:

а) у матриц Тк,р(фРк)) возможно использование нескольких способов, два из которых важнейшие; первый из них - выметание парциальных столбцов (преобразуемой матрицы) с индексами р>к, в к-ый столбец, при этом выметаются все компоненты, начиная с к-го (к = 1, 2,... ); второй способ - из условия, что в р-ом столбце после преобразования с к-ым парциальным столбцом (образуемым компонентами, начиная с к-ой), будем иметь скалярное произведение с соответствующим парциальным столбцом текущего вектора правой части (/й(к-1); /й(к) = У^ /й(к-1), /(0) = /й) равное нулю;

и далее - искомое приближенное решение

3(к) =И1(И2(... (Ик ¿(к))...)) . (171)

Модификация для систем с симметричными положительно полуопределенными матрицами. Эта модификация отличается от описанной выше (основанной на выметании строк и столбцов) тем, что в ней используется симметричная схема ортогональных преобразований (III типа):

А —Ао і — Уі АоУі

Аі — УТ АоУ

Ак — Ук Ак-іУк

(172)

5) Здесь возможно использование более “аккуратного” алгоритма, в котором аннулирование элементов осуществляется в порядке их малости по абсолютной величине, начиная с наименьшего.

к

к — 1, 2,...

при этом (М—М—п)

N

Ук — П Тк+і,р (рРк))

(173)

р=к+2

параметры матриц Тк+і,р(Ррк)) выбираются из ус

и

ловия

к+і г) к+і г)

еМї+о +£№и:

і

І = і

(174)

(к)

к+і,к + і

2

(к)2

— тах

ррк)

(к)

А к ^(к) — флк)

(к) —

к) т к + 1 і

т

0 М - (к + 1)

і

(к)

Уі{У2{... {Ук^(к)}...}} .

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

Обозначим через Ш дополнение подмножества J пар индексов, определяющих носитель матриц специального вида, до полного множества пар индексов (1 < г < N, 1 < у < М). Ясно, что имеют место неравенства

здесь а ^' - элементы матриц Ак.

В данной модификации решается последователь-

□

ность систем с матрицами Ак, получающимися из матриц Ак в виде блоков в левом верхнем углу размеров (к+1)х(к+1):

А

к + і

2

Е(СЛ)

<

Ак

2

Е(СЛ)

Ак + і 2 — Ак

Е(Л)

2

Е(Л)

(177)

1, 2,...

где Ак - матрицы, получаемые в процессе преобразования, и

(175)

при само собой понятных обозначениях. Матрицы

□ □

Ак+і получаются из матриц Ак окаймлением; решение систем (175) находится методом Холецкого; при

□

этом разложение Ак+і эффективно находится, если

□

известно разложение Ак. Нахождение приближенного решения ж(к) по вектору ■£(к) с очевидностью основывается на соотношениях:

Ак

Ак

Е(Л) к——то 2

А

Е(СЛ) к—то

2

Е

0.

(178)

Основная идея авторегуляризации состоит в том, чтобы в процессе преобразования после каждого шага осуществлять приближенную замену матрицы Ак матрицей Вк, имеющей носителем подмножество пар индексов J, и решать системы

Вк£(к) — /йк) ;

(179)

(176)

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

§ 4. В данном параграфе будет описано содержание одной из важнейших конструктивных идей - идеи авторегуляризации [Страхов, 1997g; Страхов, Тетерин, 1991а, 1991Ь, 1991с, 1993].

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

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

Переход от матриц Ак к матрицам Вк может осуществляться несколькими способами, из которых три основных:

1) элементы а(к) матриц Ак, (г,у) € Ш, просто аннулируются;

2) элементы а(к) матриц Ак, (г, у) € Ш, аннулируются, но одновременно изменяются (увеличиваются по величине) диагональные элементы6) (все или только некоторые);

6) В первой части работы уже было сказано о том, что эле-

2

2

3) элементы а(к) матриц Ак, (г, у) € Ш, аннулируются, но при этом не только увеличиваются диагональные элементы, но изменяются (уменьшаются) и некоторые внедиагональные элементы Ак, принадлежащие носителю J.

Последний (третий) способ определяет комбинированный подход авторегуляризации и непараметрической регуляризации, [см. Страхов, 19971"; Страхов, Страхов, 1999] и следующий раздел настоящей части работы.

Практически наиболее важен второй способ авторегуляризации. В нем имеется несколько вариантов.

Первый вариант состоит в том, что если

§ 5. В данном параграфе рассматривается конструктивная идея фильтрации систем линейных алгебраических уравнений [Страхов, 1991с, 19971].

Имеется три типа фильтрации (т.е. удаления неинформативных компонентов системы7):

а) фильтрация строк (или уравнений); этот вид фильтрации представляется самым простым;

б) фильтрация столбцов матрицы (или фильтрация вектора преобразованных неизвестных);

в) одновременная фильтрация строк и столбцов матрицы.

Ясно, что фильтрация столбцов осуществляется с использованием ортогонального преобразования8)

Ак = Ак + ДАк ,

(180)

где Ак - матрица специальной формы с носителем J, ДАк - матрица с носителем Ш, то используется конструкция с

В

А

Вк = Ак + аЕ

при этом с очевидностью diag 6(к) =

diag а(к) + а

(181)

(182)

Е

(г,з)еол

г=сопБ1

а(?)2 + Е №'2

(г,з)есл

3 = сопб1

менты с парами равных индексов, т.е. определяющие главную диагональ системы, всегда принадлежат подмножеству

Л пар индексов, определяющих носитель матриц специальной формы.

А = АИ, г = Итж , (184)

после чего осуществляется переход к системе

(аБ + ВтВ) = ВТй = Рй , (185)

V

где матрица В определена блочным разложением А;

(186)

V

А =

V

ДА

N ,

4

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

О

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

Второй вариант состоит в том, что элементы Ак, принадлежащие Ш, некоторым образом “проектируются” на диагональ. Например, если диагональные

А (к)

элементы Ак суть а^ , то диагональные элементы Вк находятся по соотношениям:

где к ^ N; ясно, что в (185) Б = Бт > 0 есть достаточно хорошо обусловленная матрица размера (кхк), а > 0 - параметр регуляризации. При построении матрицы В (т.е. и задании значения к) используется неравенство

V 2

ДА

Е

В

< х2

(187)

где х2 > 0 - заданная малая постоянная. Как определяются приближенные решения ж(к), читателю должно быть совершенно ясно. Из сказанного однозначно следует, что ортогональные преобразования с матрицей И - это преобразования, основанные на выметании столбцов матрицы.

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

(188)

(183)

Этот вариант также применим в ситуациях, когда

О

Ак суть симметричные положительно полуопреде-ленные, либо треугольные матрицы.

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

8) Ясно, что фактически должна использоваться после-

V

довательность ортогональных преобразований А=Ао, Ак =

V

Ак—1^, к = 1, 2,... , где матрицы и представляются в виде произведений элементарных матриц плоского вращения, параметры которых определяются из условий максимумов эвклидовых норм соответствующих столбцов.

после чего осуществляется переход к системе ^аБ + ВВт^ а = й , где матрица В определена разложением

^ м ^

А:

(189)

т

В к

і

т

ДА 1

і

(190)

и к < шт(^ М), а Б = Бт > 0 - достаточно хорошо обусловленная матрица размера (кхк), а > 0 - параметр регуляризации. При построении матрицы В (здесь опять-таки надо иметь в виду, что фактически (188) описывается последовательностью преобразований А = Ао, А к = УТАк, где матрицы Ук представимы в форме элементарных матриц плоского вращения) с очевидностью используется неравенство, по смыслу аналогичное (187):

V 2

ДА

Е

В

— х2

При этом в (194)—(195) а(к)р) суть к-ая строка текущей матрицы (характеризуемой параметрами - индексами - элементарной матрицы плоского вращения), а /(АЙР) суть к-ая компонента текущего вектора правой части.

Оба способа гарантируют, что в уравнениях с номерами 1, 2,..., к, компоненты вектора правой части содержат существенно меньшие по относительной величине погрешности, чем в исходной системе. Здесь важно подчеркнуть, что останов процесса фильтрации уравнений системы продолжается до тех пор, пока не выполнятся одновременно два неравенства: (191) и следующее

N

£ (/'?)'

і=к + і

<8- ;

— ^Ш1П 1

(196)

(к)

здесь /> й) суть компоненты вектора правой части, полученного после преобразования с матрицей Ук. Отсюда следуют условия на выбор значений параметра ак .

Ясно, далее, что одновременная фильтрация строк и столбцов состоит в последовательных - попеременных - ортогональных преобразованиях строк и столбцов матрицы - с тем, чтобы получать матрицы размеров (кхк) в левом верхнем углу, по которым далее строятся первые и вторые трансформации Гаусса, и которые по выполнении определенных критериев регуляризуются; решения регуляризованных уравнений с определенными значениями а = аг и принимаются за приближенные решения ж(г) = жаг. Авторы считают, что читатели, после всего сказанного в настоящей, второй, части работы в состоянии выписать все необходимые соотношения самостоятельно.

(191)

где х2 - заданная постоянная; однако одного этого неравенства недостаточно, см. ниже.

Как определяются решения х(к) достаточно очевидно:

хак) — АЧа . (192)

Остается описать способ определения параметров матриц элементарного плоского вращения

N

Ук — П Т

р=к + і

(193)

(здесь Ук определены как сомножители У =

ктах находится по критерию (191)).

Имеется два основных способа.

Первый способ таков:

кп

^(к;р)

(к)

/ (к;р)

2

— тіп

Е ^Рк)

Второй способ таков:

-(к;р)

(к)

— тах

Ук,

к=і

§ 6. В данном параграфе будет рассмотрена одна из важнейших конструктивных идей, используемых в рамках новой теории регуляризации систем линейных алгебраических уравнений с приближенными данными - идея регуляризации треугольных разложений квадратных невырожденных матриц. Фактически речь идет о регуляризации двух разложений:

а) разложения Холецкого матриц А = Ат > 0;

б) компактного разложения Гаусса произвольных квадратных невырожденных (в общем случае даже несимметричных) матриц.

Напомним читателю, что если А = Ат > 0, то алгоритм Холецкого обеспечивает представление А в форме

А = ЬЬт , (197)

(195) где Ь - нижняя треугольная матрица с элементами ^,7, г — у < 0. Процедура построения матрицы Ь

(194)

2

Е

2

Е

описывается соотношениями:

ат 1

¿11

лА—. ¿ji = i—. j > 1

¿ii

i-i

— ¿ipupj.

p=l

¿ji =

т-i

— E¿t

p=l

¿ii =

i-i

— E l%

p=l

(198)

i = 2. 3..... n. j = i +1. i + 2..... n

а

¿ji =

i- i

ji — ¿ip¿jp

p=l ¿ii

j>i

Регуляризованный же алгоритм Холецкого (иначе регуляризация алгоритма Холецкого) состоит в ис- равенства): пользовании соотношений (в - параметр регуляри зации, 0 < в < 1):

Регуляризованный же алгоритм Гаусса (иначе -регуляризация алгоритма Гаусса) описывается соотношениями (в - параметр регуляризации, 0 < в < 1); здесь все ац > 0, что не уменьшает общности метода (если неравенство ац > 0 нарушается, то с помощью умножения А на диагональную матрицу Б с элементами ±1 можно добиться выполнения требуемого не-

¿11(в) = \/а11 + в . Т1(в) = і!-#- . j> 1.

Uli (в) = ац + в . uij(в) = (1 — в )а1т

¿11(в) ’

¿й(в)

i- i

+в — E¿2P(в). p=i

(199)

i-i

¿ji (в)

(1 — в )аг — 2_^ ¿ip(в )¿jp(в )

_______________р=________________

¿й(в)

¿,(в) = ЧвГ. j = 2. 3..... n .

иц(в ) i-i

uii (в) аїї + в ^ ^ ¿ip ( в) upi (в). i 2. 3.....n.

p=l

(203)

i-i

uij (в) (1 ^ ^¿ip (в)upj (в) .

p=l

j-i

Таким образом находятся элементы (в) нижней

треугольной матрицы Ь(в) и составляется регуляри-зованная система

(1 в)аг ^ >¿jp(в)upi(в)

¿ji (в)

p=l

Uii (в)

L(в ^т(в)в = й

(200)

которая и решается. Из сличения (198) и (199) видны простота и экономичность алгоритма.

Далее - пусть матрица А квадратная, размера (пхп), невырожденная, причем ведущие миноры матрицы А отличны от нуля. В этом случае компактная схема Гаусса для построения разложения А в произведение двух треугольных,

A = LU

(201)

где Ь - нижняя треугольная матрица специального вида - с единицами на главной диагонали, а И -верхняя треугольная матрица, описывается соотношениями:

uii = ац. ui,

а1,

j,i

аз'1 uii ’

j

2. 3............. n

i-i

^ ^ ¿ipupi. p=l

2. 3............. n

i = 2. 3..... n. j = i + 1. i + 2..... n .

Этот алгоритм позволяет находить элементы ¿j (в) и uj (в) верхней и нижней треугольной матриц в регуляризованной системе линейных алгебраических уравнений

L^U^xe = /й . (204)

которая экономично решается. Заметим, что у метода разложения Холецкого, а также у метода разложения Гаусса есть важное свойство сохранения профиля исходной матрицы А [см. Воеводин, 1977; Воеводин, Кузнецов, 1984]. Поэтому разложения Холец-кого и Гаусса очень эффективны в случае ленточных и стреловидных матриц. Соответствующие экономичные алгоритмы читатель легко распишет самостоятельно. Вопрос выбора значений параметра в, определяющих пробные решения, достаточно ясен.

В заключение сделаем следующее принципиальное замечание: регуляризованный алгоритм Холец-кого одновременно есть обычный алгоритм Холец-кого, примененный к матрице системы

а

ij

а

а

u

а

у которой

(Бл + вЕ) + (1 - в)(А - Бл)

0 < в < 1 ,

/Й ,

(205)

Рй

ь(і), /й

где Ба суть diag ац - диагональ матрицы А. Соответственно регуляризованный алгоритм Гаусса - это алгоритм, примененный к матрице системы (205), если А не обладает свойством А = Ат > 0.

Последнее замечание позволяет увидеть пути обобщения описанных алгоритмов. Именно, вместо системы (205) можно использовать более общие системы:

если А = Ат 0, то

т

1 1

і

т

0 М - 1

і

(213)

и матрица С — ВТВ имеет структуру

С —

Сіі 0

0 ДС

(214)

(Бл + а8) + (1. — (А — Бл)

Б — БТ > 0 ,

а,в

/й ,

(206)

В силу (213)—(214) нахождение вектора w сводится к решению одного уравнения с одним неизвестным

Сіі^і — (б(і), /й^

если же А = Ат, то системы можно брать в форме и к использованию соотношения

(БЛ + аЕ) + (1 — в)(А — Бл) ха,в — /й . (207)

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

Эта идея была в краткой форме описана в первой части работы [см. Страхов, Страхов, 1999], раздел 3, § 4). Она состоит в том, чтобы с помощью ортогонального преобразования 2-го типа,

и х

т

Юі 1

і

т

0 М - 1

і

(215)

(216)

где и - та же ортогональная матрица размера (М х М), что и в (208), (211).

Ясно, что уравнение (215) элементарным образом регуляризуется - переходом к уравнению

Сіі +

Сіі

Юі,

(і),

/й

(217)

В — Аи

(208)

где и - ортогональная по столбцам (МхМ)-матрица, обеспечить выполнение следующих требований к матрице С:

ь(і)> /й

0,

2, 3,

,М

Ь(і),Ь(і) 1—0, і — 2, 3, ...,М

(209)

(210)

где а > 0 - параметр регуляризации.

Данный метод достаточно эффективен, но его можно еще улучшить - по крайней мере в наиболее важном для практики случае, когда величина

2

/ неизвестна, но известны оценки этой вели-

Е

чины сверху и снизу. Суть улучшения - в сокращении времени счета на нахождение пробных решений, которое достигается переходом от “разовой” процедуры нахождения матрицы С = ВтВ со структурой (214), к некоторому итерационному процессу, основанному на нахождении последовательности матриц

Вк ,

где 6(7), у = 2, 3,...,М, суть столбцы матрицы В. Если теперь осуществить переход от системы

А —Во ,

Вк — Вк— іИк,

к — 1, 2,...

(218)

Вю — /й ,

ИТх

где Ик суть ортогональные (М х М)-матрицы, пред-(211) ставимые в виде произведений матриц элементарных плоских вращений

к ее первой трансформации Гаусса, то получим систему

Сю — ВТВю — Вт/й — рй , (212)

Ик

м

П т

р=2

і,М ррк)

(219)

е —

X

а

и

ю

B

Ясно, что если ввести последовательность матриц что дает: (р)

B(0) = B

k = Bk-i

B(p) — B(p-i) Ti (P(k))

Dk _ Dk -4,p [Pp ) ,

(k)

b(k)

b(i)

(220)

Следовательно, искомое приближенное решение ж(k) восстанавливается по соотношению:

и обозначить через b(k,p) первый столбец матрицы

Bkp), а через b(k,p) p-ый столбец этой матрицы, то

( (k) (k))

параметры I cos pp ', sin pp ' I матрицы элементарного плоского вращения T1 ,p ^ppk) j целесообразно находить из условия

Yk,p (^^ b(k;p)) +(l-7)(b(k;P), f) — mil) , (221)

Pp )

где Yk,p > 0 - заданный весовой множитель. Целесообразно принимать

?(k) — Ui{u2{... {ükw(k)}... }} .

(228)

Нетрудно показать, что имеет место предельное соотношение

г(к) ^ ж(то) , (229)

где ж(то) есть решение системы (211)—(212).

Но имеется и второй, в некотором смысле более эффективный, способ нахождения последовательности приближенных решений (к). Именно, пусть к системам (224) применяется первая трансформация Гаусса:

7k,i

b(k;p-1) b(k-1) b(i) ,b(i)

b(ííp-i,-b(ír,0 4‘(p-i’’/>)

2 •

(222)

Нетрудно показать, опираясь на соотношения

b(k)p) = b(k)p i) cos ppk) + 6(p^' sin p.

b(k-i) , b(p)

(k)

p

(223)

b(fc;p) - _b(fc;p-1) sin p(k) + b(k-1) cos p(k)

b(p) _ b(i) sin pP + b(p) cos pP ,

что условие (221) дает для величины

t -tan pPk), - 2 < pPk) < 2

алгебраическое уравнение четвертой степени; условием (221) определяется выбор корня уравнения.

Так вот, пусть b(k)) суть 1-ый столбец матрицы Bk,

b(k) b(k;M) B B(M)

b(1) _ b(1) , Bk _Bk .

Ясно, что приближенное решение w(k) системы

Bkw(k) _ /й (224)

можно найти так:

w(k) —

(k) wi т 1 4

т

0 N - 1

4

(k)b(k) / ;i b(i) - /й

min

(k)

i

СкЮ(к) — ВТВк ю(к) — вт/й — рйк),

к — 1, 2,.... (230)

Ясно, что при этом в первой строке (и первом столбце) матриц Ск — СТ > 0 ненулевыми в общем случае будут все элементы, а не только элемент с^. Равным образом у векторов рйк) ненулевыми в общем случае будут все компоненты, а не только первая (верхняя). Поэтому опять-таки можно искать приближенное решение Ю(к) системы (230) в форме (225) и находить скаляр ю(к) из условия наилучшего среднеквадратичного приближения

,,(k)_(k)

h c(i)

(k) РЙ )

min

(k)

і

(231)

w

(225)

(k)

где w1 есть решение задачи среднеквадратичного приближения:

w

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

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

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

к

1

2

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

Итак, пусть от системы (4) мы каким-то образом перешли к системе

Вг = рЙ , (232)

где В есть (^хМ)-матрица (в общем случае N = N), Рй — ^-вектор, рй = Р + ^Р, и пусть переход от приближенного решения 5 системы (232) к приближенному решению ж системы (4) однозначно определен. Допустим теперь, что в системе (232) вектор рй имеет только одну (первую или верхнюю) ненулевую компоненту р1 й :

Рй _

pm т 1 1

т

0 N — 1

i

(233)

Способы перехода от системы (4) к системе (232) будут более подробно рассмотрены в следующем разделе статьи.

Введем последовательность ортогональных преобразований системы (232):

B_Bo ,

Bk _ Bk — 1Uk ,

k _ 1,2,...

(234)

где Ик суть ортогональные по столбцам (М х М)-матрицы, представимые в форме

м

Ик = П Т1,р(РРк0 . (235)

р=2

Нетрудно понять, что и в данном случае можно использовать последовательность преобразований в форме (220), так что будут справедливы соотношения (223)-(224). Но при этом параметры преобразования будем находить либо по соотношению

(k;p)

(1)

- ь

(k;p)

1,1

b(1k1;p)

min

p(k)

p

(236)

(первый - основной - способ), либо по соотношению

(k;p)

(1)

— ( 2 + nk,i

(k;p)

1,1

min

ppk)

(237)

9) Заметим, что в общем случае можно рассматривать постановку с произвольным (по номеру) столбцом.

где nk p > 0 - заданная величина (второй способ).

В обоих случаях параметры ^cosppk), sinppk)j

определяются из условий (236) и (237) эффективно; при этом в случае (236) целесообразно вводить переменную t по (224).

В обоих случаях приближенные решения w(k) систем

Bkw(k) _ p¿ , (238)

где векторы p¿ определены по (233), w(k) опять-таки целесообразно искать в форме (225) и величины w(k) находить из условий (226), что дает для них соотношения (227). Соответствующие приближенные решения ж(k) ив данном случае находятся по соотношению (228).

Ясно, что отношение парциальных норм вектора первого столбца, фигурирующее в (236), с ростом k и p уменьшается; если при k ^ то оно стремится к нулю, то это означает, что решение задачи (238) стремится к решению одного уравнения с одним неизвестным. Но, конечно, в общем случае отношение, фигурирующее в (236), стремится лишь к некоторому (отличному от нуля) пределу; в этом случае существует лишь псевдорешение системы (4), тождественное тому, которое находится по методу наименьших квадратов. Но этого достаточно для того, чтобы среди приближенных решений, находимых по ходу описанного итерационного процесса, содержались пробные решения (быть может, для этого надо находить “дробные приближения” - решения соответствующих задач наилучшего среднеквадратичного приближения после каждого ортогонального преобразования с элементарной матрицей плоского вращения T1p^ppk)j). Авторы не сомневаются в том, что все аналитические детали (включая проблему “дробных приближений”) читатели легко рассмотрят самостоятельно. Им рекомендуется более обстоятельно рассмотреть соотношения конструктивных идей I33) и I35), рассмотренных в настоящем и предыдущем параграфах, и их связь с конструктивной идеей I31 ), рассматриваемой в следующем разделе настоящей статьи (второй части работы).

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

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

2

2

2

2

2

дыдущем разделе. Это прежде всего конструктивные идеи:

128) идея нахождения устойчивых приближенных решений систем большой (и сверхбольшой) размерности на основе нахождения устойчивых приближенных решений совокупности подсистем существенно меньшей размерности;

129) идея использования, в случае систем большой и сверхбольшой размерности, приближенных аналитических описаний компонент искомого вектора;

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

132) идея использования редукции задачи нахождения решения системы линейных алгебраических уравнений к последовательному решению двух линейных уравнений, каждое из которых содержит одну неизвестную;

133) идея использования специальных итерационных процессов, основанных на так называемых корреляционных преобразованиях.

Кроме того, обсуждаются еще две идеи, по мнению авторов еще менее перспективные:

111) идея использования постановок условных экстремальных задач, основанных на оценках предельной относительной погрешности приближенных решений;

127) идея использования приема ускорения сходимости: построения более быстро сходящейся последовательности приближенных решений.

В связи со сказанным о степени перспективности общих конструктивных идей (часть из которых кратко рассмотрена в первой части работы, часть - в предыдущем параграфе), то безусловно следует понимать, что приведенная оценка (классификация) идей по перспективности не только субъективна, но в существенном априорна - она основывается на неких общих соображениях и интуиции (на степени “привлекательности” идей), а не на значительном вычислительном опыте.

§ 2. Суть идеи нахождения устойчивого приближенного решения системы (4) на основе нахождения устойчивых приближенных решений подсистем очень проста. С помощью случайных выборок номеров уравнений формируются подсистемы

векторы /(д) вида

(4І)2< ||/ІЕ< (41)2 , (240)

где

¿(д) ) < тіп / —

/(д)|[ — (41)2,

(241)

¿(<г))2 = ^.¿2. (¿(<г) )2= ^.¿2

шіпу N тш5 \jmaxJ ^ тах

Принимается также, что числа Жд, хотя и могут быть разными, но желательно, чтобы имели место неравенства

' М"

N —

(242)

где V - заданное число (порядка 10^30), а р есть целая часть числа р.

Устойчивые приближенные решения подсистем (239) ищутся, исходя из условной экстремальной задачи

Их

шш,

х

А д X :

( )

(243)

где И. - заданная (М х М)-матрица-регуляризатор (невырожденная и достаточно хорошо обусловленная). Как показано в первой части работы, решение задачи (243) удовлетворяет системе линейных уравнений

ИтИж = А^А(д) , (244)

в которой Жд-вектор А( д) находится из решения линейной системы

в которой

ВдВ^

В д = АдИ"

( )

(245)

(246)

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

А х

( ), г ,

я = і, 2,...,д

(239)

( )

аЖд^Т + В дБ Л х(д)

Бт В

( )

(247)

где Ад суть матрицы размеров (Жд х М), /¿' суть

Жд-векторы, /¿д) = /(д) + ¿/(д); при этом, если известны константы и в неравенствах (1), то для подсистемы (239) вводятся априорные оценки на

где а > 0 - параметр регуляризации, Жд - матрица размера ^Жд х Жд^, обратимая и достаточно хорошо обусловленная (матрица-регуляризатор). Вопрос о нахождении устойчивого приближенного решения

V

г

і

системы (247) должен быть читателю ясен - здесь все должно делаться так, как если бы эта система, в совокупности с условиями (240)-(241), была исходной.

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

~(*0

Именно, вводится последовательность ж , к = 1, 2,... , приближенных решений системы (4), конструируемых по таким соотношениям:

;(1)

г(1)

= шт ;

е ^

здесь

где

Г(к)г

і(й-1) ^(й-1) , Р

Р

5(к х) - 5«

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

§ 3. Идея 129) использования параметрических представлений компонент вектора неизвестных пре-

( 1 м-1

дельно проста. Именно, пусть < суть

система линейно-независимых функций целочисленного аргумента і, 1 < і < М. Тогда приближенное представление компонент х^- вектора ж можно брать в форме

хі «Е й=1

сйРй(з’), 3 = 1, 2,..., М .

(253)

(248)

Ясно, что если подставить выражения ж^- в систему (4), которая предварительно записывается в форме

к = 2, 3,... .

Здесь ж(к) суть устойчивое приближенное решение к-ой подсистемы, а - числовой параметр, находимый из условия

м

53 а0')х' = /й ,

5=1

(254)

(249)

5(к) = /й - А г(к) = (1 - 5(к 1) + ^5(к), (250)

где а(^) суть векторы-столбцы размерности N, то получим систему

Б„= / , (255)

в которой Б^ суть матрица размера (Ж х V), а е(^ есть ^-вектор с компонентами е^, к = 1, 2,..., V . Матрица Б^ имеет векторы-столбцы

р-- = /й - АЖ(к) . (251)

Нетрудно понять, что из (249) и (250) следует такое выражение

м

Ь(Й) = Е і=1

аСЛрй(і), к = 1, 2,. .., V .

(256)

Ясно, что во многих случаях (практически всегда!) размерность вектора е(^, при которых впервые достигается неравенство

(252)

Рк

2

/й - В„С(")

2

4іп

(257)

Итак, конструкция построения устойчивых при-„ ~(*0

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

где е(^) есть соответствующее псевдорешение (решение) системы

ВТВе(") =БТ/й , (258)

будет удовлетворять неравенству

V < М (259)

(здесь знак ^ следует понимать так: в несколько раз меньше; обычно речь идет о том, что V меньше М в 5-10 раз).

Ясно также, что:

1) на практике необходимо использовать линейные аппроксимации компонент ж^- вида (253) с последовательно увеличивающимися значениями V (но не на единицу, а на некоторое число Дv, при Дv = 5 ^ 25);

2

Е

2) при достижении размерности V искомого вектора порядка 50-100 следует решать не системы (258), а соответствующие регуляризованные системы

аЕ + В,.В^ ) С

»

В,/й

(260)

неравенствами на

N

либо в форме

и = П Т1,р ы

р=2

N-1

И = ^ Тр,м (фР) ,

р=1

(262)

(263)

Ді)

,/й )=0,

і = 2, 3,

(265)

в случае (262), и

(7й),/й) =0, і = 1, 2,...,Ж - 1

в случае (263).

Ясно, что в системе

с использованием априорных условий на вектор соответствующей невязки, порожденных априорными

2

/ ;

Е

3) после нахождения векторов с(^) необходимо находить компоненты векторов 5 (^) по соотношениям (253).

Ах = ГАж = рй = Г/й ,

(267)

при выборе матрицы и по (262), (265) вектор правой части будет иметь вид

Рй =

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

Имеется два варианта использования идеи. Первый вариант может эффективно использоваться в случае N < М. Он основан на двух процедурах: а) выбор некоторой невырожденной матрицы В размера (Ж х N) и ортогонального преобразования этой матрицы

Г = ИТВ , (261)

где и суть ортогональная по столбцам М х М -матрица, представимая либо в форме

Р1,й т 1 1

т

0 N - 1

1

Р1,й = (7(1),/й) , (268)

а при выборе матрицы И по (263), (266) - вид:

Рй =

т

0 N - 1

1

т

р^й 1

1

Р^й = (7 ^ ),/й

(269)

Ясно, что если N = М и А=АТ, то при использовании данного варианта, вообще говоря, будет АТ = А. Но это часто не так нежелательно.

Заметим еще, что в некоторых случаях данный прием эффективен и в ситуации М>М; в этом случае просто целесообразно брать матрицу размера ( М х N); при этом соотношения (261)-(267) должны очевидным образом модифицироваться (что читателю предлагается сделать самостоятельно).

Второй вариант может использоваться как при N > М, так и при М<М; при этом в случае М=М и А = АТ необходимо использовать специальную модификацию, которая описывается ниже.

Условия эффективного использования данного второго варианта таковы:

N

N

1)

= 0, ]Т/,й = 0 , (270)

в зависимости от того, какая единственная компонента у вектора

Рй = Г/й (264)

должна быть ненулевой - первая (тогда и берется в виде (262)) или последняя (тогда и берется в виде (263)).

При этом параметры матриц плоского вращения выбираются из следующих условий - векторы-строки матрицы Г, 3 = 1, 2,..., N):

2) если Б(+) суть диагональная матрица с элементами

(271)

4+) = SІgn/iл

/ (+) ____

Л,й _

/і,й

0

то

(272)

(273)

В случае, когда выполняются условия (270), исходная система (4) заменяется системой

Ах = / й

(274)

и

которая получается, если к системе (4) дописывается (соответственно сверху, или снизу) еще одно уравнение вида

N

Е;

1

г=1

N

з ~~ N з '

г=1

(275)

Соответственно в случае, когда выполняются условия (271)-(273), исходная система (4) заменяется системой

А х = / , (276)

которая получается из исходной с помощью двух процедур:

во-первых, сначала осуществляется переход к системе

Б+Ах = Б+/ = /+ , (277)

где Б+ суть диагональная матрица с элементами по

(271);

во-вторых, к системе (277) дописывается еще одно уравнение (в зависимости от того, какая именно компонента в окончательной системе должна остаться ненулевой - первая или последняя); это последнее формулируется как среднее арифметическое из уравнений системы (277) (подобно тому, как это было в случае выполнения условий (270)).

В любом случае, после получения “расширенной” системы ((274) или (277)) эта последняя подвергается ортогональному преобразованию - путем умножения на ортогональную по столбцам матрицу Ут слева (на Ут умножается как матрица системы, так и вектор правой части). Матрица У берется либо в форме

N

У=ПТ1,зР(Рр) ,

р=1

(278)

N

У

Л +1(Рр)

р=1

(279)

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

Особо следует остановиться на случае систем (4) с симметричной положительно полуопределенной матрицей А. Для систем с подобными матрицами можно использовать и первый, и второй варианты, но с дополнительными условиями:

а) в варианте первом матрица Г должна быть симметричной положительно полуопределенной;

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

Покажем теперь, что если (п +1) х (п + 1) матрица А(пр) симметричная положительно полуопре-деленная, а (п + 1)-вектор рй правой части имеет только одну (последнюю) ненулевую компоненту,

А(пр)

г = Рй

= Ут

если в результирующем векторе правой части ненулевой должна остаться первая компонента, либо в форме

т

х п

4

т

£ 1

4

Рй =

PN +1,й

п

4

(280)

(281)

если в результирующем векторе правой части ненулевой должна остаться последняя (^ + 1)-ая) компонента.

Каждая из матриц элементарного плоского вращения в (278) и (279) реализует гивенсово вращение - обеспечивает аннулирование компоненты вектора правой части, находящегося в позиции, определяемой переменным индексом у матрицы элементарного плоского вращения (индексом в случае (278),

индексом гр в случае (279)). Способ формирования

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

Именно, осуществим стандартную регуляризацию по М. М. Лаврентьеву системы (280)-(281), т.е. перейдем к системе

(аБ + А(пр)) га =

Рй

(282)

где Б = Бт > 0 - достаточно хорошо обусловленная матрица, а > 0 - параметр регуляризации. Ясно, что нахождение приближенных решений исходной системы основано на нахождении решений системы (282) при разных а = ак, к = 1, 2,... . При этом при каждом матрица системы с помощью алгоритма

0

г

Холецкого преобразуется в произведение двух треугольных матриц (размера (п + 1) х (п + 1)):

S, ak) L (S, ak) рй •

(283)

LT(S; afc)

Рй

n+1,n+1

(«fc)

(284)

CBz = рй

(285)

= Fz

(286)

Рй =

т

0 n — 1

1

т

Р№,й 1

1

(287)

2) либо

B=C

T

Но если учесть, что вектор Рй имеет только одну ненулевую компоненту (равную рп+1,й), то нетрудно понять, что (283) редуцируется к системе с верхней треугольной матрицей

ясно, таким образом, что здесь рассматривается ситуация, когда система (285) получается из исходной с помощью ортогонального преобразования и одной из трансформаций Гаусса - первой или второй.

Имея в виду соотношения (288) и (289) (т.е. что выполняется одно из них), вводим ортогональное преобразование

где /п+1,п+1(а&) есть элемент матрицы Ь(Б; а&) в правом нижнем углу.

Итак, в случае Рй с одной (последней) ненулевой компонентой оказывается необходимым решать всего одну систему с верхней треугольной матрицей, а не две (решаемые последовательно) треугольные системы. Это повышает точность и устойчивость решения.

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

§ 5. Идея 132) использования редукции задачи нахождения решения системы линейных алгебраических уравнений с М неизвестными к последовательному решению двух уравнений с одним неизвестны каждое, достаточно важна и оригинальна, но, по-видимому, все же уступает по эффективности идее 1зз) - редукции к решению одного уравнения с одним неизвестным. Тем не менее проведем достаточно детальное рассмотрение этой идеи.

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

(290)

где Ь есть нижняя треугольная, а Ьт - верхняя треугольная матрица. Следовательно, система (285) оказывается редуцированной к системе

LLTz

Рй

(291)

Ясно, что система (291) должна быть регуляризо-вана, следующим образом:

L(per)L(peO,Tza = Рй ,

(292)

т (per) т (per),T

где а суть параметр регуляризации, а La и La

суть нижняя и верхняя треугольные матрицы. В

(per)

частности, переход от матриц L к матрицам La может осуществляться по соотношению

L(per) = L + aL ,

(293)

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

LLT + LLT

S

(294)

где Б = Бт > 0 - достаточно хорошо обусловленная матрица; при этом очевидно, что условие (294) есть следствие тождества

которая получена из исходной с помощью некоторых преобразований, при этом в общем случае

где F есть (m х п)-матрица, z и рй суть n-векторы, причем

L(per)L(Per),T = (l + aL)(LT + aLT) = LLT + a(LLT + LLT)+a2LLT .

(295)

и С и В суть матрицы размеров (пхш) и (тхп), соответственно, при этом выполняется одно из двух соотношений:

1) либо

С = Вт , (288)

Конструкция, определяемая соотношениями (293)-(295), выражает сущность конструктивной идеи I26) - использования специальной “треугольной регуляризации”, см. § 3 раздела 4 первой части работы

[Страхов, Страхов, 1999].

Ясно далее, что если /1рПг) (а) суть элемент ма-

(per)

трицы La , находящейся в нижнем правом углу, то решение системы (292) с произведением двух треугольных матриц в правой части, редуцируется к решению всего одной системы:

x

т (рег),т

р5

ри,5

а -^а

4риг)( а 1(р:гч«)

(296)

т

0 п — 1

1

т

1 1

1

И.

•и— 1

(«к =

ри,5

іПРиг)(а) ’

где Ки_і(а есть верхняя треугольная матрица, и

и— 1

^а = [ П Ти — к,и — &+і() ^а

т

1 1

1

т

0 п — 1

1

Поэтому система (298)-(299) редуцируется к решению всего одного уравнения с одним неизвестным:

откуда

,,(и—1) 1,1

»1,0

(«К.

ри,5

іПРиг)(а)

ри,5

1иРиГ) (а)г1д_1) (а)

т

»1,а 1

1

т

0 п — 1

1

Заметим теперь, что можно осуществить следующее ортогональное преобразование матрицы Ьарег),т:

Ь(рег)’т = Ио (а) ,

_ / ) (1) ^ ( ) (2) ( )

И1 (а] Р п,п-1 И0 (а) Т п—1,п(р1) ,

...................... (297)

( ) (1) ( ) (2) / ) Ик(а) = Р п— к + 1,п— к Ик—1(а) Т п —к,п—к + 1 (рк) ,

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

которое имеет следующий смысл: сначала с помощью матрицы перестановок Рп—к+1;П—к переставляются (п — к + 1)-ая и (п — к)-ая строки текущей матрицы, а затем с помощью матри(цы) элементарного плоского вращения Тп—кп—к+1(рк) аннулируется элемент получившейся матрицы, расположенный в позиции (п — к + 1, п — к) (т.е. единственный элемент, располагающийся левее главной диагонали текущей матрицы). Очевидно, что с помощью матриц перестановок Рп—к+1;П—к, к = 1, 2,... п — 1, действующих также на вектор е, система (296) приводится к виду

(302)

Все дальнейшее - переход от вектора »а к вектору га и далее к вектору ха, должно быть читателю совершенно очевидно, равно как и вопрос о построении пробных решений исходной системы (4).

§ 6. Теперь переходим к рассмотрению сущности конструктивной идеи І33) - использования специальных итерационных процессов, основанных на корреляционных преобразованиях.

Сущность этой идеи предельно проста - с помощью последовательности ортогональных преобразований

А = Ао ,

А1 = АоУ^

(303)

достичь выполнения двух условий:

1)

2)

'(то) 1(0 7

/5) =0,

'(то) (1),

і = 2, 3,

/5

'(то)

(1)

/5

= тах .

М ; (304) (305)

(298)

Здесь через а|з^), ] = 1, 2,..., М, обозначены столбцы предельной матрицы . Ясно, без длинных пояснений, что фактически речь идет о матрице Ак при некотором достаточно большом к.

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

(299)

м

V = п тЦрРк;10 т1,^рРк;20 ;

(306)

р=2

(300)

(301)

иными словами, реализуется последовательность преобразований

Ак— 1 = Ак0) ,

Акр) = Акр—1)Тт,р (рРк;2))Тт,р (рРк;1)), Р = 1, 2,..., М,

(307)

АкМ) = Ак .

Параметры матриц элементарных плоских вращений выбираются из следующих условий:

и

а

2

2

2

Е

Е

Є

1) матрицы ^Рр’1^ - из условия равенства нулю коэффициента корреляции - получаемой после преобразования р-го столбца текущей матрицы, и вектора правой части;

2) матрицы Т1р ^Ррк’2^ - максимума квадрата коэффициента корреляции первого столбца получаемой после преобразования текущей матрицы, с вектором правой части.

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

Далее - нетрудно понять, что термин “корреляционные” преобразования обусловлен именно способом определения параметров матриц элементарного плоского вращения.

Описанный “мультипликативный” итерационный процесс нетрудно преобразовать в “мультипликативно-аддитивный” - по схеме, описанной в конце § 2 первого раздела настоящей статьи, см. соотношения (31)-(33) и соответствующий текст в их “окрестности”.

§ 7. Переходим к описанию сущности конструктивной идеи постановки условных экстремальных задач, основанных на оценках предельной относительной погрешности приближенных решений. Эта идея очень важная, но она по-настоящему эффективна лишь для систем с квадратными (Ж = М) матрицами А. Тем не менее, поскольку условия симметричности и тем более положительной определенности матрицы А здесь роли не играют, важность идеи достаточно ясна. Только случай систем с квадратной матрицей А и рассматривается ниже.

Оценки предельной относительной погрешности для систем общего вида рассмотрены в первой части работы [Страхов, Страхов, 1999]. Если А - квадратная (п хп)-матрица, то принимается

( N У* — Ж/й Ир

;(Л, Ж; /й) = эир ------------------—-и-------Е

Л* = / ||х||е

<

6/Уе< 6

< \/2

где

Е-ЖЛ

л11еП

XV *Г

ж

6

1/2

п

(310)

1 - п II/«Не

Отсюда с очевидностью следует первая постановка - безусловной экстремальной задачи [Страхов, 1996а, 1996Ь; Страхов, Страхов, 1999]:

ЕЛ

XV* Г

Ж

= тіп . (311)

еж

Решение этой задачи известно и выражается простой формулой [см. Страхов, 1991Ь]:

Ж = Лт (^*)2е + ллт) = (^*)2Е + Атл) 1 Ат

(312)

х = Ж/й ,

(308)

где Ж есть (пхп)-матрица с элементами р^ (которые и являются подлежащими определению, см. ниже). Соответствующая оценка предельной относительной погрешности имеет вид [Страхов, Страхов, 1999]:

которая дает одновременно и обоснование значимости частного случая общего вариационного метода А. Н. Тихонова.

Соотношением (312) фактически задается значение параметра регуляризации а : а = (V*) . Ясно поэтому, что решение (308) с матрицей Ж в форме (312) во многих случаях будет неоптимальным. Необходимо более полным образом учитывать и наличие системы (4), и соответствующей априорной информации.

Более важное практическое значение имеют постановки условных экстремальных задач.

Первая постановка: к основному условию минимизации (311) добавляются линейные условия-равенства:

СЖ/й = р , (313)

где С есть (п х М)-матрица, п < N (обычно даже п ^ М), а р - заданный вектор). Условия (313) можно получить, подвергая исходную систему ортогональному преобразованию - в соответствии с идеями фильтрации уравнений и генерации последовательности систем малой размерности, см. предыдущие разделы настоящей статьи.

Вторая постановка: к основному условию минимизации (311) добавляются одно или несколько квадратичных условий типа равенства. В простейшем случае имеем одно условие

(Е — АЖ)/й

=62

(314)

2

Е

2

1

2

Е

где ¿2 - априорно заданное число. В более общем случае имеем совокупность из V квадратичных условий типа равенства:

Рд (Е — АЖ)/й

2

6

2

я = 1,2,...,д

(315)

стеме, а 62

если значение

6/

Л = У£АИт

(316)

где V и и суть ортогональные по столбцам (п х п)-матрицы, а Ид есть диагональная матрица с элементами А, то в этом случае матрицу Ж целесообразно брать в форме

ж = ихут

(317)

X = diagЖj, 1 < і < п .

При задании Ж по (316)-(317) получаем

(318)

(Е — ЖЛ) +(v*)

(Е — Х£а) + К)2

Е

ж

(319)

X

(Е — ^аХ)рй

= 62, р = Ут/й

(320)

а условия (314) - в условия

Рд (Е — ^аХ) Рй

62

(321)

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

Далее. Совершенно очевидно, что кроме задачи минимизации (311) можно рассматривать и более простую задачу минимизации

Ж

тіп

Ж

(322)

неизвестно, а известны только

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

Решения задач (311), (313) и (311), (314), либо (311), (315) достаточно сложно, и приводить его здесь не будем. Заметим лишь, что если известно сингулярное разложение матрицы А размера (п х п),

- при дополнительных условиях типа равенств: либо (313), либо (314), либо (315). При этом, если известно сингулярное разложение (316) матрицы А, то матрица Ж берется в форме (317) и при этом задача (322) оказывается эквивалентной задаче

X

ж2 = тіп

І=1

(323)

где ж* - компоненты диагональной матрицы X. Соответственно трансформируются и условия-равенства (линейные и квадратичные), фигурирующие в условных экстремальных задачах.

Далее - в первой части работы [Страхов, Страхов, 1999] была обсуждена конструкция регуляризации, основанная на переходе от матрицы А к матрице

В = ЛИ.

1

и далее к системе

В г = /й, г = И*

(324)

(325)

где У и и - те же, что в (316), а X есть диагональная матрица с элементами ж*,

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

X = Ж/й , (326)

то экстремальная задача перейдет в задачу

ЕВ

Ж

тіп

ж

где

||в||еП

1 — п

(327)

(328)

и, таким образом, вместо п2 неизвестных р^, 1 < г < п, 1 < ] < п, получаем всего п неизвестных. Постановка соответствующих условных экстремальных задач в этом случае очевидна; в особенности упрощаются условные экстремальные задачи с квадратичными условиями типа равенств. Именно, условие (314) переходит в условие

соответственно соотношения (313), (314), (315) примут вид

СЖ/й = р , (329)

2

(Е — ВЖ) /й

62

(330)

Р,(Е - ВЖ) / = ¿2, д =1, 2,..., д. (331)

е

На этом рассмотрение идеи постановки экстремальных задач нового типа можно закончить.

2

Е

Е

2

Е

2

2

Е

2

2

2

Е

Е

2

2

*

V

Е

Е

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

6/

/й

> 10-

Лх(к)

(к)

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

р

,(к;р)

/й — ЛХ(к,р) ,

а именно

1(к;р)

= тіп

Е т (к)

(340)

(341)

Это дает для величин параметров тРк) следующие выражения:

(332)

г(*0

х(к-1,р) х(к_1’Р^ Х(к’Р-1)

*(£: — 1,р) _ х(к,р-1)

(342)

Итак, пусть ж(к) суть последовательность приближенных решений системы (4), генерируемых с помощью какого-то метода (алгоритма); при этом конструкция этого метода может быть произвольной, но при этом имеем

(Лж, / - Аж) = 0 . (333)

Первый шаг в процедуре ускорения сходимости состоит во введении новых приближенных решений

жТк) = , (334)

где

(7й, Аж(к)Л)

. (335)

Как правило, достаточно использовать не очень большие значения р,

Р < Р

(343)

где Р = 8^12. Ясно, что вычислительные затраты на процесс построения приближений ж(к,р) невелик, но и эффект от данной схемы ускорения обычно не так

велик; можно ожидать, что при выполнении условия

2

(332) величины

будут на (5 ^ 15)% меньше

величин

(к)

. Однако нет сомнений в том, что по-

При этом получается (см. первую часть работы, § 4 раздела 3):

(АХк\/й - лжТк)) = 0 (336)

добного рода выигрыш будет полностью оправдан, если удастся уменьшить основной процесс построения ж(к) хотя бы на несколько шагов (по к).

4. О нахождении устойчивых приближенных решений систем линейных алгебраических уравнений с приближенно заданной правой частью в ситуации, когда оценки сверху и снизу для ||б/ || Е неизвестны (либо очень грубые)

г«

/й — АХ«

<

/й — ЛХ(к)

<

2

Е

(337)

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

г(к) = Х(к;0), к = 1, 2,... ;

вводится последовательность приближенных реше-

г(к;р)

ний х^^' по соотношениям

§ 1. Наиболее принципиальное отличие новой теории регуляризации систем линейных алгебраических уравнений от классической теории состоит в том, что

в ней, во-первых, рассматриваются постановки, в ко-

2

торых для

6/

вообще неизвестны оценки сверху

(338)

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

пользуются оценки

6/

сверху и снизу; б) приме-

(к;р) = (1 — 7рк)) х(к-1’р) + тРк)х(к’р-1) , (339)

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

2

'Р

2

2

Е

2

V

п

Е

Е

2

Е

и

2

2

Е

Е

2

Е

§ 2. Ясно, прежде всего, что при отсутствии апри-

2

N

орной информации о

¿/

необходимо ввести си-

1)

(к)

2

е .

= 2 ; е

/«

(344)

,(к)

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

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

Итак, предлагается использовать три типа показателей качества решения:

I тип - показатели качества, относящиеся к вектору р(к) = / - Лж(к) = / - /(к);

II тип - показатели качества, относящиеся к вектору /(к) = Лж(к);

III тип - показатели качества, характеризующие меру стабильности показателей качества первого и второго типа, или изменчивость этих показателей при переходе от ж(к) к ж(к+1).

Здесь ж(к) - к-ое приближенное решение, р(к) - невязка приближенного решения, /(к) - приближенная оценка вектора полезного сигнала /.

Основные гипотезы, лежащие в основе описываемой ниже алгоритмики, таковы.

I гипотеза - компоненты вектора помехи некорре-лированы и имеют нулевое математическое ожидание, при этом числа положительных и отрицательных величин равны.

II гипотеза - вектор полезного сигнала и помехи некоррелированы.

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

IV гипотеза - допустимые приближенные решения имеют показатели качества I и II типов, изменяющиеся достаточно плавно.

Указанные гипотезы прежде всего позволяют описать необходимую совокупность показателей качества I и II типов. По вектору невязки р(к) вычисляются следующие показатели качества решения:

2)

7к =

г=1

/«

(345)

3) г* - радиус корреляции вектора невязки р(к), определяемый как наименьшее из значений «, при которых выполняется неравенство

ч 2

№-8

Ер(к)р.(+»

/«

< X2

(346)

где х2 > 0 - заданное число (обычно принимается X = 0.3);

4)

0к = — к N

(347)

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

Ясно, что чем меньше значения всех четырех введенных показателей качества решения, тем последнее выше. При этом очевидно, что допустимыми (по данному показателю качества) можно считать только те приближенные решения, для которых выполняются неравенства:

7к < 7(пред) ,

Гк < г

(пред)

(348)

0к < 0(пред).

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

Далее - в соответствии с гипотезой II вводится показатель

^лж(к), / - лж(кЛ

кк = Л----------2------; (349)

Лж(к)

/«

при этом чем меньше величина кк, тем приближенное решение считается более приемлемым; одновременно вводится условие на допустимые решения:

кк < к(пред) . (350)

В соответствии же с гипотезой III вводится показатель качества

«к =

Е-

2

М к= ЫГ

(к)

1 м

м Е гЧ а(*),р к=1

(к)

2

2

е

2

е

е

е

2

е

е

где г2 (а(4),р(к)) есть квадрат коэффициента корреляции г-го вектор-столбца а(4) матрицы А с вектором невязки р(к) приближенного решения ж(к).

Ясно, что и в этом случае допустимые по данному показателю решения выделяются условием

(пред)

< 1.

(352)

(Нетрудно показать, что в случае решения ж(0), находимого по методу наименьших квадратов, для всех г: (а(4),/э) = 0, где р - невязка для решения ж(0)).

Наконец, в соответствии с гипотезой IV вводятся величины вида

Дрк =

Рк - Рк + 1 Рк + Рк + 1 , 2

(353)

где рк - любой из введенных выше параметров:

Рк = (Пк,7к,гк,0к, кк, «к) . (354)

В случае допустимых решений для всех введенных по (354) показателей качества должно быть

Дрк > 0 ,

(355)

т.е. показатели должны монотонно убывать.

Ясно, что всего введено 12 показателей качества (6 основных и 6 дополнительных - (353), характеризующих поведение показателей качества при переходе от к к к + 1), при этом 11 из них используются для характеристики допустимости приближенных решений (с этой целью не используется только показатель Пк).

В связи с введенными показателями сделаем два принципиальных замечания.

Замечание 1. Всюду выше предполагалось, что используется метод нахождения приближенных решений, в котором не используется прием точного учета характеристического уравнения метода наименьших квадратов, т.е. условия (Лж(к),/ - Аж(к)) = 0. Если в алгоритме нахождения приближенных решений ж(к) учет данного условия осуществляется, то показатель качества кк и соответствующий ему показатель Дрк подлежат исключению.

Замечание 2. На практике достаточно часто встречаются ситуации, в которых априорно известна,

2

¿/ снизу:

пусть и очень грубая, оценка величины

е

¿/

2

>

2

шт

(356)

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

Последовательно формируются значения к, к = 1, 2, . . . , и для каждого к автоматически формируются значения параметров счета, по которым, с помощью выбранного алгоритма, находятся приближенные решения ж(к). Для каждого приближенного решения ж(к), к = 1, 2,... , вычисляются значения шести основных показателей качества решения, а при к > 2 - еще и шести дополнительных показателей качества, характеризующих изменение основных показателей качества. Если для какого-то к впервые оказалось, что по семи (из 11) показателям качества данное приближенное решение является допустимым, то оно принимается в качестве допустимого по совокупности показателей качества.

Далее процесс нахождения последовательных приближений ж(к) продолжается, при этом в общем случае допустимые (по совокупности показателей качества) и недопустимые решения могут перемежаться. Однако, начиная с некоторого значения к, число недопустимых решений начнет преобладать. Вот тогда-то и должен сработать основной критерий останова процесса; он включается после того, как появилось заданное число V допустимых решений (обычно V невелико - от 3 до 5), и состоит в следующем: если после допустимого решения появилось подряд ^ недопустимых (следует принимать ^ от 5 до 10), то процесс построения приближенных решений прекращается, а по найденным допустимым решениям по принципу усреднения находится окончательное (искомое) устойчивое приближенное решение, для которого также вычисляются значения шести основных показателей качества, которые вместе с другими характеристиками и вектором ж, выдаются пользователю.

§ 4. Описанный алгоритм нахождения допустимых приближенных решений естественно назвать стационарным - в том смысле, что в нем отсутствуют процедуры адаптации к ходу вычислительного процесса. Но можно ввести процедуры нахождения допустимых решений нестационарные, или лучше сказать - адаптивные.

В адаптивных процедурах помимо указанных выше 12 показателей, используется еще несколько величин следующего вида

Совершенно ясно, что такая оценка позволяет иметь дополнительный (и по существу принципиально важный) критерий останова процесса построения приближенных решений ж(к).

^(к) =

ирж(к)

/«

р = 1, 2,..., Р ,

(357)

е

2

е

(обычно Р = 3^5), где Ир суть матрицы-регуляриза-торы, так что ^Рк) = шр (ж(к)), где шр(ж) суть функционалы, задающие отношения предпочтения на множестве приближенных решений системы (4).

Основная гипотеза состоит в том, что по мере

(к)

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

(к)

Поэтому значения величин шр используются только на стадии останова процесса, а именно следующим образом. Вначале допустимые решений выделяются как в стационарном процессе, с помощью 11 показателей качества (пяти основных и шести дополнительных). Но при этом для каждого выделяемого таким образом допустимого решения находятся и значения величин шРк), р = 1, 2,..., Р. После того, как найдено заданное число допустимых решений (обычно - от 6 до 10), по ним находятся средние ,(к) .

значения величин шр

1

Е4к)

(358)

(к) ___________

Шсг*4 < сршр

(359)

риал данной статьи позволяет утверждать, что существует не менее 15 действительно существенно различных методов) и по каждому из них с помощью описанных выше процессов (стационарных или адаптивных) находятся соответствующие устойчивые приближенные решения, которые обозначим Ж(г), г =

1, 2,..., Д, где г - номер метода, Д - общее их число.

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

Пр,<

Ж(р) — ж (?) Е

Ж(р) + Е ж(о0 Е

2

(360)

1 < р < Д, 1 < д < Д, р = д

которые позволяют оценить степень согласованности решений, находимых различными методами. Далее находится вектор

где

я

ж дгЖ(г)

Г=1

дг > 0, дг = 1

(361)

(362)

Далее вводятся коэффициенты ср, ср > 1 (обычно принимается ср = 102 ^103) и критерии допустимости приближенных решений:

суть весовые коэффициенты. Здесь возможно как использование простого осреднения

1

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

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

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

дг = Д ,

так и более тонких приемов, например

Аг

дг =

ЕА'

г=1

где

Аг =

1

Л - Лж(г)

Далее находятся векторы

2

Е

(г)

и величины

Дж(г) — ж — .

Ж(г), г

= 1,2,...,Д

Дж

(г)

= 1, 2,..., Д .

(363)

(364)

(365)

(366)

(367)

Если среди величин <гг имеются такие, которые сильно (не менее чем в 3 раза) превышают величину

1я

Д Е<

а Д^аг

г=1

Р

д

Е

то соответствующие им решения ж(г) исключаются, и процесс обработки повторяется заново.

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

§ 6. На первый взгляд кажется, что описанный подход к нахождению устойчивых приближенных решений в ситуации отсутствия оценок сверху и снизу

2

для ¿/ является излишне усложненным и гро-

Е

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

1) по допустимым решениям внутри каждого метода;

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

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

Заключение

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

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

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

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

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

Однако уже сейчас можно утверждать следующее: во-первых, основы новой теории регуляризации систем линейных алгебраических уравнений с приближенными данными безусловно созданы (в серии работ В. Н. Страхова, частично с соавторами, [см. Страхов, 1990а, 1990Ь, 1991а, 1991Ь, 1991с, 1992а, 1992Ь, 1995а, 1995Ь, 1996а, 1996Ь, 1996с, 1997а, 1997Ь, 1997с, 1997ё, 1997е, 1997!, 1997g, 1997И, 1998а, 1998Ь, 1998с, 1998ё, 1998е; Страхов, Страхов, 1998, 1900, 1900, 1900, 1900; Страхов, Тетерин, 1991а, 1991Ь, 1991с, 1993]);

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

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

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

1) разработка общей теории и общей методологии аппроксимационного подхода к решению задач гравиметрии и магнитометрии в рамках любых интерпретационных процессов;

2) разработка проблемы редуцирования конкретных задач гравиметрии и магнитометрии, в рамках

аппроксимационного подхода, к проблеме нахождения устойчивых приближенных решений систем линейных алгебраических уравнений, либо нахождения устойчивых приближенных решений последовательности подобных систем;

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

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

Приложение

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

Конструктивная идея редукции системы к одному уравнению с одним неизвестным настолько важна, что в данном Приложении приводятся еще два ее варианта:

а) для случая исходных систем с N < М (число уравнений в системе не превышает числа неизвестных);

б) для случая исходной системы с N > М.

§ 2. Если N < М, то система линейных уравнений

Лж = / = / + ¿/ (П.1)

может быть с помощью умножения на невырожденную ^ х N)-матрицу Г, векторы-строки 7(г), г = 1, 2,..., N -1, которой удовлетворяют условию

где

7(0,/й )=0, і = 1, 2, ...,Ж — 1

приведена к виду

Вж = (7(д ),/й

(П.2)

(П.3)

В = ГА,

т

0 N — 1

1

т

1 1

1

(П.4)

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

ь = ВИ,

= ит

(П.5)

где И есть ортогональная по столбцам (М х М)-матрица, а Ь есть нижняя треугольная матрица размера (Ж х N). С помощью преобразования (П.5) система (П.3) редуцируется к форме

7(« ),/й

(П.6)

т.е. фактически к решению одного уравнения с одним неизвестным:

1д,д ю« = (7(Д ),/й)

(П.7)

где > 0 - элемент матрицы N в правом нижнем углу Ь. Уравнение (П.7) эффективно регуляризуется - переходом к уравнению, содержащему параметр регуляризации а :

1«,Д + 1------ ) Ю«,а = ( 7(«), /й

¿Д,Д / У

откуда

1д,д (7(«), /й

Ю«,а =

(П.8)

(П.9)

Ясно далее, что приближенное решение исходной системы (П.1) находится с помощью соотношений

(П.10)

т

0 N — 1

1

т

ЮД,а 1

1

(П.11)

Эта конструкция достаточно экономична и устойчива. Можно даже утверждать, что она в некотором смысле наилучшая в случае недоопределенных систем.

Одновременно очевидно, что редукция исходной системы (1) к системе с вектором правой части всего с одной ненулевой компонентой - последней, может

ж

а

ж

а

быть построена с помощью и приема, основанного на введении дополнительного уравнения и далее -соответствующего ортогонального преобразования -способом, описанным в § 4 раздела 3 настоящей статьи. Этот способ будет даже несколько более экономичным, но по устойчивости, по-видимому, он будет уступать изложенному выше.

§ 3. Перейдем к рассмотрению случая N > М (т.е. N +1 > М). Здесь первый этап редукции состоит во введении дополнительного уравнения (которое объявляется первым) и последующего ортогонального преобразования, как это описано в § 4 раздела 3 статьи. В итоге получается система вида

Сж = е, (ПЛ2)

где С есть + 1) хМ)-матрица, - скаляр, и

т

1 1

1

0 т N

1

(П.13)

V:

1 0

0 ф

(П.15)

В первом варианте осуществляется переход к “усеченной” системе - путем отбрасывания элемента в позиции (М + 1, М), т.е. элемента км+1,м матрицы

о

Н. Обозначим матрицу усеченной системы через Н; ясно, что она имеет размер N х М; тогда системе (П.17) будет соответствовать (приближенно, но, как показывает анализ, с весьма высокой точностью) система:

Н ж = 0г ее, (п.18)

где

т

1 1

1

т

0 М — 1

1

(П.19)

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

НИ = И, И ж = ю

(П.20)

Дальнейшее достаточно просто - осуществляется ортогональное преобразование матрицы С

Н^ТС , (П.14)

где ^ +1) х ^ + 1)-матрица V имеет структуру

где И. - верхняя треугольная матрица, И - ортогональная по столбцам (МхМ)-матрица,

м-1

И 11 Тм—р, м—р+1 (Рр) р=1

(П.21)

и Ф есть ортогональная по столбцам (NхN)-матрица. Ясно, что

ФТе = е , (П.16)

и поэтому система (П.12) переходит в систему

Нж = 0Й е, (П.17)

где Н по определению есть (М + 1) х М -матрица типа Хессенберга, т.е. у нее ненулевыми являются (в общем виде) все элементы на диагонали, ниже главной, т.е. элементы с индексами (г,^), г -] = 1.

Далее возможно использовать два различных варианта редукции системы (П.17): в первом варианте обеспечивается приближенная (но высокоточная!) редукция к одному уравнению с одним неизвестным; во втором варианте обеспечивается приближенная редукция к переопределенной системе из двух уравнений с одной неизвестной (здесь точность приближенного решения, по-видимому, будет несколько выше).

и каждая из матриц элементарного плоского вращения, фигурирующая в правой части (П.21), реали-

е

зует гивенсово вращение, аннулирующее элемент Н с индексами (М -р, М -р + 1). Ясно, что при ортогональном преобразовании (П.20) система (П.18) переходит в систему

Ию = #й е

(П.22)

которая, в силу выражения для е, см. (П.19), приводит к одному уравнению с одним неизвестным:

Г11^1 =

(П.23)

Это уравнение, как мы уже хорошо знаем, элементарно регуляризуется (а > 0 - параметр регуляризации)

откуда

ю1,а — 2

Г11#й

г21 + а

(П.24)

(П.25)

После этого вектор жа восстанавливается по соот-

е

а

ношениям

Î

w1,a 1

1

Î

0 M - 1

1

(П.26)

Учитывая специфику матрицы С, см. (П.30), легко находим, что система (П.31) эквивалентна переопределенной системе вида

(П.33)

Г11 e¿

X W1 =

q1 0

= Uwa .

Далее нетрудно найти значение

ра f(5 Axa .

С очевидностью имеем 2 / ч 2

E

(П 27) Эта система снова регуляризуется таким образом:

(П.34)

2 a

, где Г11 + e¿

E Г11 X W1 a =

(П.28) q1 0

Pc

riiwi,a - e^j + (hM +1,

M XM,a

Є|а2

r21 + a

e2

r21 + a

+

h2 r2 e2

hM+1,M ' 11 e<5

r21 + a

2

X K2

+ KhM+1,M гп)

G = HU

Î

R M

1

Î

Q 1

1

(П.30)

где Q есть (1 х М)-матрица, т.е. попросту вектор-строка. Имеем систему

Gw = e¿ e ,

где

Î

1 1

1

Î

0 M

1

(П.31)

(П.32)

Используя процедуру среднеквадратичного приближения

a

Г11 +-----W1,a - I

Гц)

2

22 + wi,a qi = ,min

w1,a

(П.29) получаем

W1,

rn( r21 + a)e¿ (r21 + a)2+q2 r21

(П.35)

(П.36)

где K2 есть произведение квадратов синусов вида sin pp, фигурирующих в преобразованиях переменных. Ясно, что K2 при больших М(М > 100) исчезающе малая величина, вот почему рассмотренный первый вариант является вполне приемлемым.

Во втором варианте элемент +i,m не отбрасы-

вается, а снова осуществляется ортогональное преобразование с той же самой матрицей U, определенной по (П.21), причем все параметры преобразований (cos pp, sin ) берутся теми же самыми. В этом слу-

чае имеем

т.е. отличное от выражения для ад1а в первом варианте. Дальнейшее читателю должно быть предельно ясно - имеют место аналоги соотношений (П.26)-(П.29).

§ 4. Имеется еще и ряд других вариантов реализации конструктивной идеи редукции к одному уравнению с одним неизвестным, в частности - для случая симметричных положительно полуопределенных матриц. Рассмотрение их заняло бы слишком много места. Здесь же лишь обратим внимание читателя на тот принципиальный факт, что все идеи ^^-Гэб) по существу очень тесно связаны, и что это именно те идеи, которые могут (и, по мнению авторов, должны) обеспечить подлинный прогресс в случае систем линейных алгебраических уравнений большой и сверхбольшой размерности. При этом читатель должен в полной мере осознать и тот факт, что все конструктивные идеи в той или иной степени связаны

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

Литература

w

a

и

P

а

2

2

2

a

2

Воеводин В. В., О методе регуляризации, Журн. выч. математики и мат. физики, 9, (3), 673-675, 1969.

Воеводин В. В., Вычислительные основы линейной алгебры, 303 с., Наука, Москва, 1977.

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

Иванов В. К., Васин В. В., Танана В. П., Теория линейных некорректных задач и ее приложения, 206 с., Наука, Москва, 1978.

Лаврентьев М. М., О некоторых некорректных задачах математической физики, 91 с., Наука, Новосибирск,

1962.

Лисковец О. А., Вариационные методы решения неустойчивых задач, 343 с., Наука и техника, Минск, 1981.

Морозов В. А., Об устойчивых методах решения систем линейных алгебраических уравнений, Вычислительные методы линейной алгебры, с 57-62, Изд. СО АН СССР, Новосибирск, 1973.

Морозов В. А., Регулярные методы решения некорректно поставленных задач, 360 с., МГУ, Москва, 1974.

Морозов В. А., Регулярные методы решения некорректно поставленных задач, 239 с., Наука, Москва, 1987.

Страхов В. Н., Об устойчивых методах решения линейных задач геофизики, Ч. I, Постановки и основные конструктивные идеи, Изв. АН СССР, Физика Земли, (7), 3-27, 1990а.

Страхов В. Н., Об устойчивых методах решения линейных задач геофизики, Ч. II, Основные алгоритмы, Изв. АН СССР, Физика Земли, (8), 37-64, 1990Ь.

Страхов В. Н., Алгебраические методы в линейных задачах гравиметрии и магнитометрии, Основные идеи, Докл. АН СССР, 319, (1), 143-146, 1991а.

Страхов В. Н., Алгебраические методы в линейных задачах гравиметрии и магнитометрии, Постановки экстремальных задач, Докл. АН СССР, 319, (2), 342-346, 1991Ь.

Страхов В. Н., Метод фильтрации систем линейных алгебраических уравнений - основа для решения линейных задач гравиметрии и магнитометрии, Докл. АН СССР, 320, (3), 595-600, 1991с.

Страхов В. Н., О линейных некорректных задачах гравиметрии и магнитометрии, Условно-корректные задачи математической физики и анализа, с. 176-204, Наука, Новосибирск, 1992а.

Страхов В. Н., Алгоритмы редуцирования и трансформации аномалий силы тяжести, заданных на физической поверхности Земли, Интерпретация гравитационных и магнитных полей, с. 4-82, Наукова думка, Киев, 1992Ь.

Страхов В. Н., Основные направления развития теории и методологии интерпретации геофизических данных на рубеже ХХІ столетия, Ч. 2, Геофизика, (4), 10-20, 1995а.

Страхов В. Н., Геофизика и математика, Изв. РАН, Физика Земли, (12), 4-23, 1995Ь.

Страхов В. Н., Методологические проблемы теории и практики интерпретации данных в прикладной геофизике, Вопросы методологии интерпретации геофизических данных в прикладной геофизике, с. 4-20, ОИФЗ РАН, Москва, 1996а.

Страхов В. Н., Методологические особенности интерпретации данных гравиразведки и магниторазведки,

Вопросы методологии интерпретации геофизических данных в прикладной геофизике, с. 110-123, ОИФЗ РАН, Москва, 1996Ь.

Страхов В. Н., Развитие общей теории интерпретации геолого-геофизических данных, Основные достижения ОИФЗ РАН за 1992-1996 гг., Т.1, с. 48-53, ОИФЗ РАН, Москва, 1996с.

Страхов В. Н., Общая схема и основные итоги развития теории и практики интерпретации потенциальных полей в СССР и России в XX веке, Развитие гравиметрии и магнитометрии в XX веке, с. 98-120, ОИФЗ РАН, Москва, 1997а.

Страхов В. Н., Третья парадигма в теории и практике интерпретации потенциальных полей (гравитационных и магнитных аномалий), Ч. I, Вестн. ОГГГГН РАН, N0. (1(1)), 163-198, 1997Ь.

Страхов В. Н., Третья парадигма в теории и практике интерпретации потенциальных полей (гравитационных и магнитных аномалий), Ч. II, Вестн. ОГГГГН РАН, N0. (2(2)), 55-82, 1997с.

Страхов В. Н., Общая теория нахождения устойчивых приближенных решений систем линейных алгебраических уравнений с приближенно заданными правыми частями и матрицами, возникающих при решении задач геофизики, Вопросы теории и практики геологической интерпретации гравитационных, магнитных и электрических полей, с. 38-42, ОИФЗ РАН, Москва, 1997а.

Страхов В. Н., Математический аппарат, используемый при конструировании алгоритмов нахождений устойчивых приближенных решений систем линейных алгебраических уравнений, возникающих в задачах гравиметрии и магнитометрии, Вопросы теории и практики геологической интерпретации гравитационных, магнитных и электрических полей, с. 43-75, ОИФЗ РАН, Москва, 1997е.

Страхов В. Н., Экстремальные задачи, непараметрическая регуляризация и фильтрация в теории нахождения устойчивых приближенных решений систем линейных алгебраических уравнений с приближенно заданными правыми частями и матрицами, Вопросы теории и практики геологической интерпретации гравитационных, магнитных и электрических полей, с. 76-86, ОИФЗ РАН, Москва, 19971.

Страхов В. Н., Обобщенные ^Р-алгоритмы нахождения устойчивых приближенных решений систем линейных алгебраических уравнений с приближенно заданной правой частью, возникающие при решении линейных задач гравиметрии и магнитометрии, Вопросы теории и практики геологической интерпретации гравитационных, магнитных и электрических полей, с. 87-89, ОИФЗ РАН, Москва, 1997g.

Страхов В. Н., Линейные задачи гравиметрии и магнитометрии в постановках, адекватных геофизической практике, Вопросы теории и практики геологической интерпретации гравитационных, магнитных и электрических полей, с. 89-104, ОИФЗ РАН, Москва, 1997Ь.

Страхов В. Н., Современное состояние и перспективы развития теории интерпретации гравитационных и

магнитных аномалий, Вопросы теории и практики геологической интерпретации гравитационных, магнитных и электрических полей, с. 4-35, Воронеж, 1998а.

Страхов В. Н., Взгляд в будущее (перспективы развития теории интерпретации геофизических данных в начале ХХ1 века), III Междунар. конф. “Новые идеи в науках о Земле”, Избр. докл., с. 39-50, МГГА, Москва, 1998Ь.

Страхов В. Н., Г. А. Гамбурцев - один из основоположников гравиразведки, Григорий Александрович Гамбурцев, Воспоминания, очерки, статьи, с. 153-176, ОИФЗ РАН, Москва, 1998с.

Страхов В. Н., Третья парадигма в теории и практике интерпретации потенциальных полей (гравитационных и магнитных аномалий), Ч. III, Вестн. ОГГГГН РАН, 1, (3), 100-152, 1998а.

Страхов В. Н., Что делать? (О развитии гравиметрии и магнитометрии в России в начале XXI века), 24 с., ОИФЗ РАН, Москва, 1998е.

Страхов В. Н., Страхов А. В., О регуляризации метода наименьших квадратов, Геофиз. журн., (6), 18-38, 1998.

Страхов В. Н., Страхов А. В., О нахождении устойчивых приближенных решений систем линейных алгебраических уравнений большой и сверхбольшой размерности, возникающих в линейных задачах гравиметрии и магнитометрии, Геофиз. журн., 21, (1) 20-39, 1999.

Страхов В. Н., Страхов А. В., О решении линейных задач гравиметрии и магнитометрии, Геофиз. журн., 21, (1), 40-63, 1999.

Страхов В. Н., Страхов А. В., Аппроксимационный подход к решению задач гравиметрии и магнитометрии, I, Основная вычислительная проблема - регуляризация систем линейных алгебраических уравнений, Российский журнал наук о Земле, 1, 4, 271-299, 1999. (ИКЬ: http://eos.wdcb.ru/rjes/rje99020/rje99020.htm)

Страхов В. Н., Страхов А. В., Обобщения метода наименьших квадратов и регуляризованные алгоритмы нахождения устойчивых приближенных решений систем линейных алгебраических уравнений, возникающих при решении задач геофизики, I, Геофиз. журн., 21, (2), 325, 1999.

Страхов В. Н., Страхов А. В., Обобщения метода наименьших квадратов и регуляризованные алгоритмы нахождения устойчивых приближенных решений систем

линейных алгебраических уравнений, возникающих при решении задач геофизики, II, Геофиз. журн., 21, (3), 3-17, 1999.

Страхов В. Н., Тетерин Д. Е., Линейные трансформации гравитационных и магнитных аномалий в случае многоэлементных съемок при произвольных сетях наблюдений, Докл. АН СССР, 318, (3), 572-576, 1991a.

Страхов В. Н., Тетерин Д. Е., Метод авторегуляризации при решении задач линейных трансформаций гравитационных и магнитных аномалий, Докл. АН СССР, 318, (4), 867-871, 1991b.

Страхов В. Н, Тетерин Д. Е., О методе авторегуляризации для решения линейных задач гравиметрии и магнитометрии, Докл. АН СССР, 318, (4), 871-874, 1991c.

Страхов В. Н., Тетерин Д. Е., О методах сингулярных разложений матриц, No. 1915-93, Деп. в ВИНИТИ, 1993.

Танана В. П., Методы решения операторных уравнений, 155 c., Наука, Москва, 1981.

Тихонов А. Н., О решении некорректно поставленных задач и методе регуляризации, Докл. АН СССР, 151, (3), 501-504, 1963a.

Тихонов А. Н., О регуляризации некорректно поставленных задач, Докл. АН СССР, 153, (1), 49-52, 1963b.

Тихонов А. Н., О некорректных задачах линейной алгебры и устойчивом методе их решения, Докл. АН СССР, 163, (3), 591-594, 1965.

Тихонов А. Н., Арсенин В. Я., Методы решения некорректных задач, 286 с., Наука, Москва, 1979.

Тихонов А. Н., Гончарский А. В., Степанов В. В., Ягола А. Г., Регуляризующие алгоритмы и априорная информация, 198 с., Наука, Москва, 1985.

Уилкинсон Дж., Райнш К., Справочник алгоритмов на языке Алгол, Линейная алгебра, 248 с., Машиностроение, Москва, 1976.

Фаддеев Д. К., Фаддеева В. Н., Вычислительные методы линейной алгебры, Изд. 2-е, 734 с., Физматгиз, М.-Л.,

1963.

Strakhov V. N., Vorontsov S. V., Seki T., Linear inversions in helio-seismology: testing new regularisation techniques for linear algebraic equations, Cong’94 Helion- and Astero-Seismology, ASP Conf. Ser., Vol. 76, 1995.

(Поступила в редакцию 25 июня 1999.)