Экономичный метод многократного решения расширенных регуляризованных нормальных систем уравнений Текст научной статьи по специальности «Математика»

Вычислительная математика

УДК 519.654

ЭКОНОМИЧНЫЙ МЕТОД МНОГОКРАТНОГО РЕШЕНИЯ РАСШИРЕННЫХ РЕГУЛЯРИЗОВАННЫХ НОРМАЛЬНЫХ СИСТЕМ УРАВНЕНИЙ

В. В. Долишний, А. И. Жданов

Самарский государственный аэрокосмический университет им. ак. С. П. Королёва (национальный исследовательский университет),

443086, Россия, Самара, Московское ш., 34.

E-mails: VDolishniy@mail.ru, ZhdanovAleksan@yandex.ru

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

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

Введение. Многие прикладные задачи регрессионного анализа и математической физики сводятся к решению системы линейных алгебраических уравнений (СЛАУ) вида

Ах = Ь, АеГх", 6 = & + £eR"\ т^п, (1)

в которой точное значение вектора правой части b неизвестно, а известно его значение & = & + £, где £ — возмущение.

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

Требуется найти из системы (1) вектор х такой, чтобы норма разности \\х—жЦг, где ж — вектор решения невозмущённой системы, была минимальной.

В реальных задачах характерной особенностью матрицы системы (1), как правило, является плохая обусловленность (спектральное число обусловленности [1] К2{А) l/emach, где £mach ~ машинное эпсилон) и достаточно большая размерность (m, п > 104). Матрица А может иметь неполный ранг (или неполный численный ранг) и, возможно, т^> п.

Александр Иванович Жданов (д.ф.-м.н., проф.), зав. кафедрой, каф. прикладной математики. Долишний Василий Владимирович, аспирант, каф. прикладной математики.

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

где а > 0 — параметр регуляризации, согласованный с возмущениями в исходных данных.

Эффективность решения задачи регуляризации (2)) зависит от выбора параметра а. Для выбора параметра регуляризации а известны различные подходы [2]. При этом часто возникает необходимость решения системы (2) при различных значениях параметра а.

Известно, что наиболее эффективным численным алгоритмом решения системы (2) при различных значениях параметра регуляризации является алгоритм, основанный на сингулярном разложении матрицы системы А [3]. Однако решение задачи большой размерности с его помощью вызывает значительные трудности при реализации на ЭВМ. Например, использование программ из пакета МАТЬАВ при размерности матрицы т,п ~ 104 проблематично, так как требует значительных объёмов оперативной памяти. Однако сингулярное разложение — самое надёжное в смысле точности вычислений.

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

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

1. Метод расширенных регуляризованных нормальных систем. В [4] был

предложен подход, основанный на замене системы (2) эквивалентной расширенной регуляризованной нормальной системой (РРНС)

где ш = л/а; Im-> In ~ единичные матрицы размерности т х т и п х п соответственно; z = rw-1; г = Ъ — Ах — вектор невязки системы (1). Обозначим

Преимуществом решения РРНС (3) вместо системы (2) является то, что

(2)

(3)

К2(Аш)

шах

К2(АтА + W2/)) 2

где кг(^4) —спектральное число обусловленности матрицы А; ¿тах(^4)> ¿тт(^4) — максимальное и минимальное сингулярные числа матрицы А, соответственно [4].

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

2. Решение расширенной регуляризованной нормальной системы при помощи двухдиагонализации матрицы. Представим матрицу А в виде

где В € М"1*™ — нижняя двухдиагональная матрица, II € Ктхт и V € Мгахга — 0рТ0г0нальные матрицы.

Разложение (4) можно получить, используя, например, модификацию процесса двухдиагонализации, описанного в [5].

Очевидно, что разложение (4) существует для произвольной матрицы А. Лемма. Система (3) эквивалентна системе

где IV = иТг] V = Утх, с1 = IIт1); В, II, V удовлетворяют (4).

Доказательство. Система (3) является совместной и определённой [4].

( иТ о \ / шіт ивуТ \ ( и о \ { ит О V О Ут ) V УВтит -ш1п у ^ 0 ^ У V 0 Ут ) \ х ) ~

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

Перейдём к решению расширенной регуляризованной нормальной системы (5). Обозначим

А = иВУт

(4)

. Легко проверить, что матрица Н — ортогональная. Умножая систему (3) слева на Нт, учитывая разложение (4) и то, что ННТ = I, имеем

или после преобразований

В

Вт —ш1п

(6)

а ы = (й — Ву)ш 1

Введём следующую перестановку натуральных чисел от 1 до т + п:

1, т + 1,2, т + 2,3, т + 3,..., п, т + п, п + 1, п + 2,..., т — 1, т. (8)

В случае, когда т = п, ряд (8) обрывается на числе 2п.

Пусть Р — матрица перестановок, соответствующая перестановке (8). Обозначим _

Тш = РТВШР Є м(™+")х (">+"). (9)

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

Учитывая (7), (9), перепишем систему (5) в виде

где у

ТшУ = /,

“')-7=рТ(о

Рассмотрим два случая: т > п и т = п.

Если т > п, то матрица Тш будет выглядеть следующим образом:

(10)

= РТ

Тш =

0(тп—п—1)х(2п+1) где ~ трёхдиагональная матрица вида

/0(2га+1) X (т—п— 1) \

) Є ^(т+га)х(т+га)) (11)

^1щп—п— 1 )

( и &1Д

¿»і,і -и £>2,1

&2Д Ш £>2,2

V

є

!га+1)х(2га+1)

(12)

£>га+1,га Ш /

В случае т = п матрице В будет соответствовать матрица Тш

\

в 2га X 2га

/ ^ бід

£>і,і £>2,1

£>2,1 и ¿»2,2

V

є

(13)

-ш }

Для простоты и единообразия описания дальнейшего процесса решения в случае тп = п дополним систему (10) фиктивным уравнением

£>га+1,га + Шу2п+1 — ¡2п+1,

(14)

в котором положим Ъп+1>п = ¡2п+1 = 0. Тогда можно считать, что матрица при т = п имеет вид

=

о

Тш

1x2 га

О2

:гах 1

ш

!га+1)х(2га+1)

(15)

Учитывая (11), (14), (15), систему уравнений (10) можно переписать в виде

Яшу{1) = /(1),

Шу(2) = /(2)

и решение (10) сводится к решению системы

яшу{1) = 1{1),

где

у(1) =

У\

У2п+1

.„(2) =

У2п+2

Уп-\-п

/(1) =

( /1 \

\ І2п+1 /

/(2) =

(16

^ /2га+2 ^

\ /п+т /

Непосредственное решение системы (16) как системы с трёхдиагональной квадратной матрицей может привести к большим погрешностями, так как главная диагональ матрицы как следует из (12), (13) (15), содержит параметр регуляризации со, значение которого может быть мало.

Опишем метод получения устойчивого решения системы (16). Перепишем

(16) в виде

/ ш ¿»1,1 Ь\л -Ш Ь2л

\ ( VI \

У2

V

Ьп+1,п Ш / V У2п+1 /

( /1 \

/2

V І2п+1 )

Проведём исключение из этой системы переменных с чётными индексами. Для этого из строк с номерами 2г — 1 и 2г + 1, 1 ^ г ^ п вычтем строку с номером 2г, умноженную на —Ьігі/со и на —Ъ^\,г/со соответственно. Умножая полученные уравнения на со, получим систему

Сши =

(/2г Ьі іУ2і— 1 іУ2і-\-\ Л .

У2І = ^---------------!------—------------!--------) , 1 < * < П,

(17)

где

=

/ С1 С1 с; с2 с'2

и =

сп сп+1

( У 1 \

Уз

V У2га+1 /

/ 51 \

92

V 9п+і /

Сі = д + и2, сп+1 = Ь2+1га + ш2, Сі = б2^! + + 002, 2 < і < щ

с'і = Ьг,Д+м, 1 ^ г ^ Щ 91 = /і^ + /2&1,Ъ 9п-\-1 = /2га+1^ + ¡2пЬп+1,п,

9і = І2І-1Ш + /гг—2¿»г,г— 1 + /гЛ.г, 2 ^ І ^ П.

Определим числа

i

Vi = Yl(k, l^i^n + 1,

k=1

где

Cl = 1; a = \ 'lU’-il1, 0 < \bk-i,k-i\ < \bk,k~i\, 2<k^n + 1_

1, \bk-l,k-l\ ^ |6fc,fc_i|,

Пусть D = diag(r?i, r?2, • • •, i?ra+i) € R(ra+1)x(ra+1) - диагональная матрица. Рассмотрим следующую СЛАУ:

Сшф = д, (18)

где Сш = СшО, ф = D~1u.

Теорема. Для системы (18) метод прогонки будет корректен и устойчив при любом ш > 0.

Доказательство. Докажем корректность. Прогонка будет корректной, если знаменатели Aj прогоночных коэффициентов ¿¿,1^г^п + 1, не обращаются в нуль. Для матрицы Сш имеем

с (b1Ab2,l)r]2 b1Ab2,l А _ 2 2. n.

1 /72 I o\ г2 I 2 5 ^1 л UJ 0,

Clm (¿1 ! + Ш2)г] 1 bf ! + LO2

Г c'2fj3 (62,263,2)%

О 2 = — ~

C2V2 + c[r] i<5i + Щ>2 + Ш2)г) 2 + (61,162,1)^1^1

(62,263,2)% 62,263,2

(¿2,1 + ¿2,2 + w2)% - ^61,162,1^^2^) т bh + 62,2 + W2 -

Д2 = б| I + 6о 2 + W2 - .о1’1,2’1» ^ 6о 2 + W2 > 0.

■Сз)

я, +

Далее, по индукции, пусть для всех k, 1 ^ к ^ г, выполняется

fc6fc A i

г- bk,kbk+l,k л д ^ , 2 I 2 ^ п

—---------7------Cfc+1) Ак ^ Ък к + U) >0.

Для к = г + 1 имеем

^+1^+2 _ (6j+l,i+l6i+2,i+l)i?i+2

¿г+1 — —

Ci+iVi+i + с'г?А (б2,, * + б2,, i+1 + uj2)r]i+1 + (6M6i+M)r?A

ъ2 ъ2

л 1,2 , 1,2 I 2 г,г г+1,г ^

Аг+1 — &І+М + г+1,г+1 + Ш-------------------------^

н

ь2 ь2

\ 1,2 , 7,2 , 2 г,г г+М х і2 , 2 ^ п

^ &І+М + г+1,г+1 + ^ — , 2 , 2 ^ ^г+1,г+1 + Ш

'-'А 4 I ^

Таким образом,

Л — ^г+1,г+1^г+2,г-|-1 л » > /,2 , 2 ^ г,

Ог+1 — д Сг+2) ¿\+1 Рг+1,г+1 + ^ >0

Аг+1 и, следовательно,

с- Ьг,Л+ 1,г Л ^ ,2 , 2 ^ п

¿г — Т Сг+1) А» ^ 6^ + Сс> >0

^■г

выполняется для всех г = 1,2,...,п + 1.

Корректность прогонки для системы (18) доказана.

Докажем устойчивость. Прогонка будет устойчивой, если для всех г = 1, 2,... , п + 1, выполняется условие |<^| < 1.

Учитывая полученные ранее неравенства, имеем

N - « ^N<*1 < «Г-* < !• ¡ = 1,2,...,«; |<5„+.1=0.

Аг ъ2^ + со2 Ъ^ + со2

Таким образом, для системы (18) метод прогонки будет корректным и устойчивым для любого со > 0. □

Отметим, что иногда требуется найти не всё решение ^ ~ ^ системы

(5), а только вектор невязки гю. В этом случае в соответствующей системе

(17) нет необходимости находить неизвестные ?/2г> % = 1,2,..., гг, поскольку, как следует из (8), (10), эти компоненты вектора у соответствуют вектору V системы (5).

3. Вычислительная сложность алгоритма. Проведём анализ вычислительных затрат предложенного алгоритма. Наиболее трудоёмким является приведение матрицы исходной системы (1) к двухдиагональному виду (получение разложения (4)). Известно два варианта этой процедуры: двухдиагонализа-ция Хаусхолдера и Д-двухдиагонализация [5]. Первый даёт лучший результат при т ^ 5/3 п, второй —при т ^ 5/3 п. В случае, когда явное формирование матриц и и V не требуется, как в предложенном алгоритме, двухдиагонали-зация Хаусхолдера требует 4тп2 — 4п3/3 + 0(п2) флопов, затраты на выполнение Е-двухдиагонализация составляют 2тп2 + 2п3 + 0(п2) флопов [5].

Следует заметить, что решение системы (2) получается из решения системы (6) при помощи ортогонального преобразования, которое сохраняет евклидову норму вектора. Это обстоятельство часто позволяет использовать его, не вычисляя регуляризованного решения исходной системы. В этом случае количество арифметических операций, необходимых на одну итерацию решения системы (6), составляет 17п. Если всё-таки необходимо при каждом параметре регуляризации вычисления решение исходной системы (2), то

Время вычисления параметра регуляризации методом перекрёстной значимости

~~ ~____размерность метод — _____ 512 х 512 1024 х 1024 1536 х 1536 2048 х 2048

на основе РРНС, сек 11,19 66,45 201,05 446,64

на основе SVD, сек 15,06 112,66 382,77 918,58

каждая итерация потребует 17п + п2 операций. Все остальные вычисления, связанные с реализацией алгоритма, требуют порядка 0(п2) операций и не вносят существенного вклада в общие вычислительные затраты при m,n 1.

Итого вычислительная сложность предложенного алгоритма составляет Атп2 — 4п3/3 + т(Пп[+п2]) + О(п2) (при m ^ 5/3 п) или 2тп2 + 2п3 + + т{Пп[+п2}) + О(п2) (при m ^ 5/3 п) флопов, где г —количество вычисляемых решений.

Оценим необходимый для реализации алгоритма объём оперативной памяти. Разложение (4) может быть записано на место исходной матрицы при условии, что матрицы U и V хранятся в факторизованном виде [5]. Их явное формирование для реализации предложенного алгоритма, как упоминалось ранее, не требуется. Более того, матрица U нужна лишь для получения правой части системы (5). Таким образом, предложенный алгоритм не требует дополнительной оперативной памяти и вычисления по нему могут быть реализованы на месте матрицы исходной системы.

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

Одним из подходов к вычислению параметра регуляризации, не требующим априорной информации, является метод перекрёстной значимости [6].

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

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

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

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Беклемишев Д. В. Дополнительные главы линейной алгебры. М.: Наука, 1983. 336 с. [Beklemishev D. V. Additional Chapters of Linear Algebra. Moscow: Nauka, 1983. 336 pp.]

2. Морозов В. А. Регулярные методы решения некорректно поставленных задач. М.: Наука, 1987. 240 с. [Morozov V. A. Regular Methods for Solving Ill-Posed Problems. Moscow: Nauka, 1987. 240 pp.]

3. Chung J., Nagy J.G., O’Leary D.P. A weighted-GCV method for Lanczos-hybrid regularization// Electron. Trans. Numer. Anal., 2008. Vol. 28, no. 12. Pp. 149-167.

4. Жданов А. И. Об одном численно устойчивом алгоритме решения систем линейных алгебраических уравнений неполного ранга// Вестн. Сам. гос. техн. ун-та. Сер. Физ,-

мат. науки, 2008. №1(16). С. 149-153. [Zhdanov A. I. About one numerical stable algorithm for solving system linear algebraic equations of defect rank // Vestn. Samar. Gos. Tekhn. Univ. Ser. Fiz.-Mat. Nauki, 2008. no. 1(16). Pp. 149-153].

5. Голуб Дж., Ван Лоун Ч. Матричные вычисления. М.: Мир, 1999. 548 с. [Golub G., Van Loan С. Matrix Computations. Moscow: Mir, 1999. 548 pp.]

6. Bjdrck A. Least Squares Methods / In: Handbook of Numerical Analysis. I: Solution of Equations in Rn. Part 1.; eds. P. G. Ciarlet and J. L. Lions. Amsterdam: Elsevier/North Holland. Pp. 466-647.

Поступила в редакцию 19/X/2011; в окончательном варианте — 27/XI/2011.

MSC: 15А06

ECONOMY METHOD FOR MULTIPLE SOLVING OF AUGMENTED REGULARIZED NORMAL SYSTEM OF EQUATIONS

V. V. Dolishniy, A. I. Zhdanov

S. P. Korolyov Samara State Aerospace University

(National Research University),

34, Moskovskoe sh., Samara, 443086, Russia.

E-mails: VDolishniy@mail .ru, ZhdanovAleksan@yandex.ru

In this work, a new numerical algorithm for solving augmented regularized normal system of equations with several regularization parameters is proposed. The computational costs and the required amount of RAM of the proposed algorithm are analyzed. The proposed numerical algorithm is compared with the algorithm based on the singular value decomposition.

Key words: ill-posed problems, augmented regularized normal system of equations, numerical algorithm.

Original article submitted 19/X/2011; revision submitted 27/XI/2011.

Vasiliy V Dolishniy, Postgraduate Student, Dept, of Applied Mathematics. Alexander I. Zhdanov (Dr. Sci. (Phys. & Math)), Head of Dept., Dept, of Applied Mathematics.