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

Итерационный метод решения системы линейных уравнений с использованием вейвлет-фильтров

В.А. Есаулов

Южно-Российский государственный политехнический университет (НПИ) имени М.И.

Платова

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

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

Введение

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

Построение итерационного метода решения СЛАУ с использованием

вейвлет-фильтров

Термин "вейвлет" (wavelet) в переводе с английского означает "маленькая (короткая) волна". Теория вейвлетов дает более гибкую технику

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

Одним из наиболее известных алгоритмов кратномасштабного анализа является алгоритм Малла [3]. В этом алгоритме два фильтра - сглаживающий А и детализирующий Э, рекурсивным образом применяются для получения данных для всех доступных масштабов. Как правило, используются фильтры с конечным импульсным откликом, в которых отсчеты анализируемого сигнала, попавшие в небольшое "окно", умножаются на заданный набор коэффициентов, полученные значения суммируются, и окно сдвигается для расчета следующего значения на выходе. Графическая схема алгоритма Малла приведена на рис. 1.

Рис. 1. Схема алгоритма Малла вейвлет-анализа сигнала Кроме ДВП известно пакетное вейвлет-преобразование [6], которое отличается от ДВП повторной фильтрацией детализирующих коэффициентов. Пакетное вейвлет-преобразование позволяет лучше контролировать процесс разделения спектра исходного сигнала на части, но

значительно увеличивает вычислительную сложность. Кроме того, вейвлет-пакетное разложение содержит избыточное число вейвлет-коэффициентов, которое можно сократить, если организовать поиск «наилучшего дерева». Рассмотрим использование вейвлетов в целях решения систем СЛАУ. Пусть дана система

Ах=Ь (1)

где А - матрица системы; Ь - столбец свободных членов; х - столбец неизвестных.

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

ЖАх = ЖЬ (2)

Соотношение (2) можно рассматривать как условие равенства вейвлет-образов правой и левой частей СЛАУ (1).

Учитывая свойства матрицы Ж, можно перепишем (2) как: (ЖАЖТ ))х = А и = ЖЬ = Ь (3)

где и = Жх - вейвлет-образ решения (1).

А = ЖАЖТ - предобусловленная главная матрица системы (1). где Ь = ЖЬ - вейвлет-образ правой части (1).

Матрица А, во многих случаях является сильно разреженной и имеет иерархическую ленточную структуру, что позволяет применять к (3) известные эффективные алгоритмы [7, 8].

Общая схема организации итерационного процесса решения СЛАУ имеет вид [9]:

= ^ + Т Р - АХк )= Хк + тЛ (4)

где тк - шаг итерационного процесса; гк - вектор невязки;

Хк, Хк+1 - величины, описывающие приближения решения СЛАУ. Из (4) вытекает, что как решение Хк, так и невязку гк можно рассматривать как временные вектор-функции. Это означает, что к ним можно применить вейвлет-преобразование и тем самым выявить закономерности их изменения во времени.

Пусть Ь и Н - аппроксимирующий и детализирующий вейвлет-фильтры порядка п , действующие на последовательность я = .

следующим образом:

п-1 п-1

ЬЗк = Е агЯк+г , Н}!к = Е вк+г (5)

'к

г=0 г=0

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

Зк = Ьзк + Нзк (6) Формула (6) отражает первый уровень декомпозиции сигнала я в алгоритме Малла. Тогда для т-го уровня декомпозиции получим соотношение:

—-1

^ = Ь-Зк +Е НЬзк , - > 1 (7)

]=0

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

Из (4) следует, что от характера изменения невязки гк будет зависеть

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

Гк+1 =(Р -Ткл)гк =Тгк (8)

Тогда изменение невязки на т-ом шаге может быть выражено непосредственно из (8) следующим образом:

гк+т =Ж -ткА)тгк =гтгк (9)

С учетом (8), (9) действие фильтров на гк будет иметь вид:

=ЕаГЧ , НГк =1дг ггк (10)

1=0 1=0 1=0 1=0

Перепишем гк в (4) с учетом эквивалентного представления (7)

хк+1 = хк + Ткгк = хк +Тк

т-1 ^

г т

Ьт + Н ^ V гк , т > 1 (11)

;=0 )

Предположим, что вклад одного из вейвлет-фильтров в вейвлет-разложения мал настолько, что им можно пренебречь. Исключив его из рассмотрения в (11), получим аппроксимацию вида:

Хк+1 = Хк +Ткгк х Хк + ТкТт гк , т > 1 (12)

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

Задачу о аппроксимации (11) можно рассматривать как задачу адаптивной фильтрации. В работе в качестве одного из вариантов ее решения был предложен подход [10]. Согласно ему, исключаемая составляющая определялась по уровню отношения ее энергии к энергии невязки гк. Если

это отношение было менее, чем заданная 8Е, то составляющая исключалась из дерева разложения.

Из (10) видно, что коэффициенты {а }П=01, {Д}-! определяют сглаживающие и детализирующие свойства фильтра. Из вида (10) можно сказать, что V и Н прогнозируют значения невязки согласно (9) на число шагов в пределах своего порядка. Это позволяет предположить возможность улучшения качества решения СЛАУ. Для конкретной СЛАУ точность ее решения и скорость сходимости процесса (12) будет определяться типом вейвлет-фильтра и уровнем разложения т.

Шаг тк для (12) можно рассчитывать несколькими способами. Одним из наиболее эффективных является расчет по методу минимальных невязок

[9]. В соответствии с ним шаг спуска должен доставлять минимум норме

к+1 •

1К+1Ц = \\Ь - лХк+11 = 1К -ТкЛТ-гк\ ^ тт (13)

Задача (13) в общем виде является нелинейной, причем ее вид задается как порядком фильтров (10), так и величиной уровня разложения -.

Решение (13) можно упростить, задавая значение тк в тт в (9), (10) постоянным, либо заменяя ее на известное на предыдущей значение шага тк -1. Тогда для оптимальной величины шага получим:

Тк =(Т4), (14)

\\ЛТ-Гк\\

Таким образом, формулы (12), (14) задают способ решения СЛАУ, основанный на многоуровневой вейвлет-аппроксимации невязки. Он отличается от изложенного в (1)-(3) тем, что объектом вейвлет-преобразования является невязка, а не элементы СЛАУ.

Вычислительный эксперимент

Целями проводимого эксперимента было исследование погрешности решения получаемого по методике (12), (14), в зависимости от типа вейвлет-фильтра и уровня вейвлет-декомпозиции. Максимальное число итераций задавалось как 3000. Критерием его останова являлось отношение ||Хк+1 - Хк||2 <8 , где 8 - допустимая погрешность. Величина погрешности

задавалась как 8 = 10-8. Пороговое значение 8Е задавалось как 8Е = 5 • 10-2 и принималось одинаковым для всех уровней разложения. Если значение погрешности составляло \Хк+1 - Хк|| 2 <8Е, то значение 8Е заменялось на

8Е =|\Хк+1 - Хк|| 2 .

Также рассматривалась СЛАУ, полученная предобуславливанием (1) по методу наименьших квадратов [9]:

ATAx = ATb (15)

Инструментальным средством реализации предложенного метода являлась среда Matlab. В качестве тестовой задачи бралась СЛАУ, полученная C. B. Shaw в [10]. Порядок системы задавался равным 64, а ее

число обусловленности составило cond(A)« 8.1 • 1018.

В табл. 1 приведены данные изменения относительной погрешности решения от уровня декомпозиции для задачи [10] c использованием вейвлета Хаара.

Таблица № 1.

Погрешность решения задачи [10] для случая вейвлета Хаара

Число уровней m 1 2 4 6 8 10

Погрешность решения (1), % 6.1 7.1 11.41 11.33 11.45 11.4

Погрешность решения (15), % 4.73 8.8 4.7 8.4 8.7 8.9

Помимо вейвлета Хаара, для сравнительного анализа использовались симметричные вейвлеты, используемые в широком спектре приложений [6].

В табл. 2 приведены данные изменения относительной погрешности решения от уровня декомпозиции для задачи [10] с использованием симметричного вейвлета второго порядка sym2.

Таблица № 2.

. Погрешность решения задачи [10] для случая вейвлета sym2

Число уровней m 1 2 4 6 8 10

Погрешность решения (1), % 6.7 6.4 7.2 6.23 11.04 11.05

Погрешность решения (15), % 4.71 8.1 9.2 12.8 15 4.69

В табл. 3 приведены данные изменения относительной погрешности решения от уровня декомпозиции для задачи [10] с использованием симметричного вейвлета четвертого порядка sy—4.

Таблица № 3.

. Погрешность решения задачи [10] для случая вейвлета sy—4

Число уровней т 1 2 4 6 8 10

Погрешность решения (1), % 5.2 6.59 7.1 7.17 6.56 10.4

Погрешность решения (15), % 7.8 15.3 16.62 16.5 16.8 14.2

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

В табл. 4 приведены данные изменения относительной погрешности решения от уровня декомпозиции для задачи [10] с использованием вейвлета bгor1.1.

Таблица № 4.

. Погрешность решения задачи [10] для случая вейвлета bгor1.1

Число уровней т 1 2 4 6 8 10

Погрешность решения (1), % 7.01 9.4 11.3 11.37 11.4 11.2

Погрешность решения (15), % 4.7 6.03 7.5 8.39 8.7 8.9

В табл. 4 приведены данные изменения относительной погрешности решения от уровня декомпозиции для задачи [10] с использованием вейвлета bгor1.5.

Таблица № 5.

. Погрешность решения задачи [10] для случая вейвлета bгor1.5

Число уровней т 1 2 4 8 10 12 14

Погрешность решения (1), % 5.5 5.69 6.6 7.8 9.3 10.8 11.1

Погрешность решения (15), % 4.83 16.1 15.4 9.5 9.1 9.6 7.54

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

Заключение

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

Литература

1. Шарый С.П. Курс вычислительных методов. Учеб. пособие. -Новосибирск: Новосиб. гос. ун-т. , 2014 г. , 279 с.

2. Годунов С.К. Решение систем линейных уравнений. - Новосибирск: Наука, 1980. - 178 с.

3. Добеши И. Десять лекций по вейвлетам.- Ижевск: НИЦ «Регулярная и хаотическая динамика», 2001.- 464 с.

4. Болдырев С.В. Применение гибридных самоорганизующихся нейронных сетей и быстрого дискретного вейвлет-преобразования для построения систем классификации сигналов // Инженерный вестник Дона, 2012. № 2 URL: ivdon.ru/ru/magazine/archive/n2y2012/773.

5. Астапов К.А. Анализ градиента для нейронных сетей с вейвлет-разложением целевого вектора // Инженерный вестник Дона, 2009. № 1 URL: ivdon.ru/ru/magazine/archive/n1y2009/104.

6. Матвеев Ю.Н., Симончик К.К., Тропченко А.Ю., Хитров М.В. Цифровая обработка сигналов. Учебное пособие по дисциплине «Цифровая обработка сигналов». - СПб: СПбНИУ ИТМО, 2013. - 166 с.

7. David L. Gines LU Factorization of Non-Standard Forms and Direct Multiresolution Solver / David L. Gines, G. Beylkin, J. Dunn. - Appl. Comput. Harmon. Anal. - 1998. - № 5(2). - pp. 156-201.

8. Юдин М.Н., Фарков Ю.А., Филатов Д.М. Введение в вейвлет-анализ: учеб.-практическое пособие. Моск. геологоразв. акад. М.,2001. - 72 с.

9. Турчак, Л.И., Плотников П.В. Основы численных методов: Учебное пособие. - 2-е изд., перераб. и доп. / Л.И. Турчак, П.В. Плотников/ - М.: физматлит, 2003. - 304 c.

10. A. Pikrakis, T. Giannakopoulos, and S. Theodoridis, A speech/music discriminator of radio recordings based on dynamic programming and bayesian networks. Multimedia, IEEE Transactions on, 10(5), (2008). pp. 846-857.

11. C. B. Shaw, Jr., Improvements of the resolution of an instrument by numerical solution of an integral equation, J. Math. Anal. Appl., 37 (1972), pp. 83112.

References

1. Sharyj S.P. Kurs vychislitel'nyh metodov. Ucheb. Posobie [The course of computing methods. Tutorial.]. Novosibirsk: Novosib. gos. un-t. , 2014 g., 279 p.

2. Godunov S.K. Reshenie sistem linejnyh uravnenij [Solving systems of linear equations]. Novosibirsk: Nauka, 1980. 178 p.

3. Dobeshi I. Desjat' lekcij po vejvletam [Ten Lections on Wavelets]. Izhevsk: NIC «Reguljarnaja i haoticheskaja dinamika», 2001. 464 p.

4. Boldyrev S.V. Inzenernyj vestnik Dona (Rus), 2012. № 2 URL: ivdon.ru/ru/magazine/archive/n2y2012/773.

5. Astapov K.A. Inzenernyj vestnik Dona (Rus), 2009. № 1 URL: ivdon.ru/ru/magazine/archive/n1y2009/104.

6. Matveev Ju.N., Simonchik K.K., Tropchenko A.Ju., Hitrov M.V. Cifrovaja obrabotka signalov. Uchebnoe posobie po discipline «Cifrovaja obrabotka signalov» [Digital signal processing. Tutorial on the discipline «Digital signal processing»]. SPb: SPbNIU ITMO, 2013. 166 p.

7. David L. Gines LU Factorization of Non-Standard Forms and Direct Multiresolution Solver. David L. Gines, G. Beylkin, J. Dunn. Appl. Comput. Harmon. Anal. 1998. № 5(2). pp. 156-201.

8. Judin M.N., Farkov Ju.A., Filatov D.M. Vvedenie v vejvlet-analiz: ucheb.-prakticheskoe posobie [An introduction to wavelet analysis: a training manual]. Mosk. geologorazv. akad. M., 2001. 72 p.

9. Turchak, L.I., Plotnikov P.V. Osnovy chislennyh metodov: Uchebnoe posobie [Basics of numerical methods: Textbook]. 2-e izd., pererab. i dop. L.I. Turchak, P.V. Plotnikov. M.: Fizmatlit, 2003. 304 p.

10. A. Pikrakis, T. Giannakopoulos, and S. Theodoridis, A speech/music discriminator of radio recordings based on dynamic programming and bayesian networks. Multimedia, IEEE Transactions on, 10(5), (2008). pp. 846-857.

11. C. B. Shaw, Jr., Improvements of the resolution of an instrument by numerical solution of an integral equation, J. Math. Anal. Appl., 37 (1972), pp. 83112.