CC BY
53
ПРЯМОЙ МЕТОД РЕШЕНИЯ СИСТЕМЫ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ НА ОСНОВЕ ВЕЙВЛЕТ-РАЗЛОЖЕНИЯ
О. В. Филиппов
Рассматривается прямой метод решения СЛАУ, в основе которого лежит нестандартное представление в вейвлет-базисе матрицы системы. Дается краткое введение в вейвлеты и описывается нестандартное представление оператора в вейвлет-базисе. Обсуждается, почему прямой метод при использовании вейвлетов становится эффективным, даже когда рассматриваются плотнозаполненные матрицы.
Введение
Существует много доступных методов решения больших сильноразреженных систем (CPG) линейных уравнений. Наиболее популярными методами решения для таких систем являются итерационные методы.
Прямые методы решения СРС могут также использоваться для решения таких систем. Их эффективность зависит от числа дополнительных ненулевых членов, полученных в ходе LU-факторизации матрицы системы. Напомним, что некоторые разреженные системы не имеют разреженных LU-факторов. В работе [2] предложен метод построения разреженного представления матриц, которые изначально плотно-заполнены.
Далее будет показано, что, если использовать нестандартное представление в вейвлет-базисе матрицы системы, соответствующие LU-факторы будут также сильноразреженными, это делает эффективным прямой метод решения СЛАУ.
Перечислим основные определения и свойства вейвлетов.
Вейвлет-функция (// строится по масштабирующей функции <р, которая обладает двумя свойствами:
а) система функций {(р(х-к),к е Zj ортогональна в L2(R),
б) функция является решением функционального уравнения
<p(x) = \[l'^ihkip(2x-k). (1)
к
Тогда для вейвлет-функции ц/ верно представление
i//(x) = j2Yi(-lfh]„k<p(2x-k). (2)
к
Наибольший интерес представляют компактно-сосредоточенные функции (р и ц/ . В этом случае фильтр Ик конечен.
Пусть К0 =span{cp(x-k),k&Z) = {f £l}{R),f(x) = YjCk<P{x-k)’Y,ck-00} и
Pj к{х) = 2/^(p^2Jх-к^. Через Vj обозначим линейную оболочку span^.^2(p{ll -к е zj.
Из соотношения (1) следует, что
...c:V_, cV0cV, с...
Определим подпространство W как ортогональное дополнение V] в VJ+], т.е.
VJ+1 =Vj®)Wj. Если \]v] = L2 (R) и Pi V} = {0}, то доказывается в [1], что
j j
Система функций |ц/; к (х)| и |(р0 к (д:)| является ортогональным базисом в пространстве
]} (Я). В частности, Уп - У0 Ф0<у<п~і ^ ■
Тогда система функций
>>о
где
уо7 (*)=IX* (*+/)и (х)=2Х* (х+1),
/ І
является ортогональным базисом в І? (0,1). Определим проектор Р на подпространство
и <2. на подпространстве УУ
PJ■.L2{R)^VJ
где 0] - /у, -Р]. Для линейного оператора Т: і2 (0,1) -> і2 (О.і) оператор Т] =Р]ТР] будет называться дискретизацией оператора Т масштаба у. В базисе г/>' к (х,у) с оператором Т(] можно
связать матрицу Т0 размера 2і х 27.
Введем ортогональную матрицу
Ч-.
(3)
где V' , и <2;_| - матричное представление проекторов Р , и 2 х соответственно. Пусть
■ 'Рх и у} -гР]у - дискретизации элементов хиуиз I? (0,1). Для того, чтобы свести
Т]Х}=У)
к масштабуу-1, применим (3) и (4)
№-№-№н*>Н™ну,)-
Получим
Л}_ Ъ-х' 4-і” Ч~
Си ТН. .5П. Л.
где с1;=Я}у, Э]='Р]У, = Ъ-Iх в то время
(4)
(5)
(6)
матричные
операторы
А]-\ : 1-^-1 -4;_] = Я- -1^2,-:
= бу- -1^-1
Л -> У]-\ <7-1 = ^у-
' I ТУ-1 = ^-
Матрицы Л}_\, Б/ . и С/содержат уточняющую информацию об операторе Т ^. На следующем шаге матрица Т;_х извлекается из (6), и к ней применяется ортогональная матрица У\^_2 > как и в (5)' Эта процедура продолжается на всех масштабах / = п -1,...,0. Результатом является нестандартная форма (Ш-форма) исходной матрицы Т0. Для у0 =3 Ш-форма матрицы Т0 показана на рис. 1. Здесь и 3,$ связаны соотношениями
^ =СД,у = 7о- !»•••,0,
и на последнем шаге
50 = + % •
Для того, чтобы конвертировать ^ ^ в вейвлет-представление |{^у}0< ] ,«01,
следует воспользоваться алгоритмом:
1. ^,=0,^ = 0.
2. = О, + 5Н), 5, = Т} (5;_, 4- 1Н).
3. Тогда = й? + , а на последнем шаге з0 - ,уГ) + 70.
Пусть система функций у/к (?) ортогональна и
СО
*(0='Е с*Мг)»
к--со
где коэффициенты ск вычисляются по функции
с* = |/(0^*(0Л-
Если ц/к ([) имеют нулевые моменты
0 = ту/к (1)Л,0<т<п-1,
тогда функция / хорошо аппроксимируется многочленом <п-1. Коэффициенты разложения будут малыми или нулевыми. Поэтому вейвлеты, имеющие нулевые моменты, могут использоваться для «сжатия» гладких данных.
Если Т является интегральным оператором и его ядро - гладкая функция вне диагонали, то матрица, полученная «дискретизацией», имеет много малых элементов вне диагонали. Такие матрицы после «сжатия» имеют разреженное представление в вейвлет-базисе. Поэтому N8-форма матрицы после «сжатия» становится сильно разреженной.
А в2
Л 3
<Л
А) А)
Со Ъ
й2
52 б2
4
~*0 *0
Рис. 1. Матричное уравнение в нестандартном представлении.
Алгоритм
Рассмотрим алгоритм решения СЛУ методом нестандартного представления в вейвлет-базисе. Его общий вид можно увидеть на рис. 2.
1. Исходное уравнение имеет вид: А*х = Ь, где А - матрица коэфициентов, Ь - вектор свободного члена, ах- вектор неизвестных
____________________________________________<2_____________________________________________
2. Приведение матрицы коэффициентов Л в нестандартную форму (М8-Рогт) и представление исходной матрицы ввиде произведения Ь*и, где I - это нижнетреугольная нестандартная форма, а II- верхнетреугольная нестандартная форма
_________________________________________________________________________
3. Приведение вектора свободного члена в нестандартную форму (ТЧБ-Рогт). В результате получаем нестандартную форму Ьо исходного вектора Ъ
______________________________________________________________
4. В ходе этих преобразований получаем уравнение Ь*и*х0 = Ьо, где Хо - это нестандартная форма искомого вектора
__________________________________а_________________________________
5. Решение уравнения Ь*у = Ьй методом нестандартного прямого прохода
______________________________________<2________________________________________
6. Решение уравнения и*х0 = у методом нестандартного обратного прохода
_______________________ъ±_________________________
7. Приведение вектора х0 к стандартной форме
Рис. 2. Алгоритм
Прокомментируем шаг 2. Выше уже описывался способ построения нестандартной формы матрицы, но для данного алгоритма предлагается совместить его с Ш-раэложением (более подробное описание данного подхода можно увидеть в работе [2]). Результирующая процедура получится более эффективной, чем если бы нестандартное и Ш-разложения применялись независимо.
В начале расчитаем вейвлет-разложение оператора Тп для «первого» масштаба:
Тп = Л-1 + + СпА + Тп-\ >
где
Л-1 =(2п-\Тп()п-\,
вп-1= Оп-\ТпРп-\>
Си-1 = Рп-\Т„(2п-\’
Тпл = Рп-\ГпРп-\ •
Используя стандартное Ш-разложение, получим:
Л-1 = А1-Л-!’
и зададим Ап_х = 1пА и Ап. Имея Ап_х и Ап_х, можно получить и Сп_х решая
Л-1^и-1 = Вп~] >
Сп-\А1-\=Сп-\
с помощью стандартных методов прямого и обратного прохода.
Для продолжения на следующем масштабе необходимо вейвлет-разложение оператора Тп_х и
СП_ХВП_Х. Эти процедуры могут быть скомбинированы, если вначале вычесть СП_ХВП_Х из Тп_х, а затем рассчитать
ги_1-сяЧ4ч=(Л-2-Д,-2)+(5я-2-Д,-2)+(с1|_2-си_2)+(2;_2-ги_2)) (7)
где
An-2~A„-2=e„-2(r„-i-C„-l£n-l)i2n-2>
B/t-2 ~Bn-2 -Qn^i^n-l ~Cn-\Bn-\)Pn-2>
Cn-2 ~СП-2 =Рп-г{Тп-\ ~Cn-lBn-l)Qn-2’
Pn-2 ~Tn-2 ~ P/i-2 {P>i-\ ~Cn-\B„-i)Pn-2-
Используя стандартное LU-разложение, получим
Л-2 ~ Л-2 = Ai-2^«-2>
где Ап_2 -Ln_2 и Ап_2 -Un_2 ■ Имея Д,_2 И Л_2, МОЖНО получить Вп_2 и С п_2 решая
Л-2 Д1-2 ~ -®И-2 ~ -®и-2 >
Сп-гЛ-2 = С„_2 -Сп_2,
с помощью стандартных методов прямого и обратного прохода.
Для продолжения необходимо рассчитать вейвлет-разложение Гд_2, СП__2ВП_2 и дополнительного члена Тп_2 из равенства (7).
^1-2 ~ Ai-2-^л-2-Тп-2 =(Л-3 ~Л-3.) + (-®и-3 " Д1-з) + (Сп-з -Си_з) + (Гп_з -7^_3),
где
Л-З ~ Л-З - бл-З (^я-2 " ^п-2^п-2 ~ Тп-2 )би-3>
■®и-3 ~-®п-3 = бя-3 (^в-2 ~Сп-гВп-2 — ^и-2)^п-3’
С„_з - С„_з = Рй_з (гп_2 - СП_2ВП_2 -Тп_2 ) gn_3,
тп.3-тп_3=рп_3(тя_2-сп_2вп.2-Т„_2)р„_у
Процесс продолжается до масштаба я, где рассчитываются члены Г0 и f0 .
Программа
На основе вышеприведенного алгоритма нами была разработана программа1 для решения СЛАУ. В качестве языка программирования был выбран Fortran90, для работы с сильно разреженными матрицами была использована библиотека SPARSKIT2.
Основную работу выполняет процедура solve UsingNSForm, на вход которой подается матрица коэффициентов, вектор свободного члена, векторы коэффициентов материнской и вейвлет-
функций, количество вейвлет-преобразований, которое необходимо будет проделать над матри-
цей коэффициентов в процессе получения нестандартной формы, и переменная, в которой будет возвращен вектор неизвестных. Стоит отметить, что все матрицы и векторы передаются в сжатом формате (так называемый, формат CSR).
Работу процедуры solve Us ingNSForm можно разделить на 5 частей.
1. Расчет нестандартного представления для матрицы коэффициентов. Это работу проделывает процедура createNSForm2Dim. Ей на вход подается матрица, которую необходимо преобразовать, векторы коэффициентов материнской и вейвлет-функций, количество преобразований и переменные, в которых будут возвращены верхнетреугольное
1 В момент написания данной статьи программа проходила процесс регистрации в ФГНУ «Государствен-ном координационном центре информационных технологий» Федерального агентства по образованию.
и нижнетреугольное нестандартные представления исходной матрицы. Данная процедура выполняется, четко следуя шагу 2 вышеприведенного алгоритма.
2. Расчет нестандартного представления для вектора свободного члена. Для этого используется процедура createNSFormlDim. Ей на вход подается вектор свободного члена, векторы коэффициентов материнской и вейвлет-функций, количество преобразований и переменная, в которой будет возвращен результат - нестандартная форма вектора. В процессе выполнения данной процедуры исходный вектор претерпевает заданное количество вейвлет-преобразований, причем получающиеся векторы запоминаются в переменную-результат.
3. Решение первого уравнения (см. шаг 5 алгоритма) в нестандартной форме методом прямого прохода. Для этого используется процедура тРот’ап^ЗиЬмИиНоп. На вход подается нижнетреугольная нестандартная форма, вектор свободного члена в нестандартном представлении, векторы коэффициентов материнской и вейвлет-функций и переменная, в которой вернется результат решения.
4. Решение второго уравнения (см. шаг 6 алгоритма) в нестандартной форме методом обратного прохода. Эту работу выполняет процедура тВасЪлюгйЗиЪьШиНоп. На вход подается верхнетреугольная нестандартная форма, результат, полученный на предыдущем шаге, векторы коэффициентов материнской и вейвлет-функций, переменная, в которой вернется результат решения.
5. Обратное вейвлет-преобразование для получения вектора решения в стандартной форме. Для этого используется процедура т\'ег$еЛУа\'е1е1Тгап81 Оип. на вход которой подается вектор материнских и вейвлет-коэффициентов первого масштаба искомого вектора, векторы коэффициентов материнской и вейвлет-функций, переменная, в которой вернется результирующий вектор.
Для проведения всех стандартных действий над матрицами и векторами используются процедуры и функции из библиотеки SPAR.SK.IT2. Она содержит процедуры для преобразования матриц из стандартного формата в «сжатый» и обратно, решения СЛАУ методами прямого и обратного прохода, Ш-разложения матрицы, умножения и сложения матриц и так далее.
Численные результаты
В этом разделе представлены численные результаты обработки нескольких примеров для демонстрации быстродействия и точности алгоритма и программы. Для тестирования был выбран компьютер с процессором 1п1е1(!<) РепПит(К) 4 СРи З.ОООН/. с 2048МЬ оперативной памяти.
Пример 1. Для первого примера была выбрана матрица вида:
Аи =
і*),
где N-2п ~ размерность матрицы, /,/ = 1.,2" и С = і/2" . Правая часть матричного уравнения подбиралась так, чтобы существовало известное точно решение.
Данная СЛАУ была решена при помощи нестандартного представления и при помощи библиотеки Ьараск. Из этого пакета была выбрана процедура ОЯСЕЗУ, реализующая один из наиболее часто используемых итерационных алгоритмов решения СЛАУ. Более подробно об этом алгоритме можно узнать из документации, расположенной по адресу http://www.netlib.org/lapack.
Результаты можно увидеть в табл. 1.
Таблица 1
N 4>(Ж) 12(Ж) Т Ьараск Д, (Ьараск) Ь2(Ьараск)
29 0,68 1,74 • 10~4 4,59 • 10“4 0,78 О і VI 1,36-10~14
2ш 3,35 4,25-10-4 1,19-10—3 6,18 3,33-Ю"15 2,48-Ю-14
2П 19,32 6,53-10-4 2,25 -10“3 48,18 9,05-10-15 8,39 -10-14
В первом столбце этой таблицы находится размерность обрабатываемой матрицы, во втором время в секундах работы алгоритма на основе нестандартного представления, в третьем и четвер-
том ошибка данного алгоритма, в пятом время работы в секундах процедуры из библиотеки Ьа-раск, в двух оставшихся - ошибка этой процедуры.
Пример 2. Во втором примере использовалась матрица вида:
'N-
), і* ]
/ 1 У
1, і = /
размерность матрицы и /,у' = Правая часть данного матричного уравнения
где N = 2"
подбиралась так, чтобы существовало известное точно решение.
Результаты можно увидеть в табл. 2, построенной аналогично табл. 1.
Таблица 2
N Tns 4с (NS) L2(NS) Т мLapack Lj., (Lapack) Ь2 (Ьараск)
29 0,73 1,46 -10~4 1,0Ы0“3 0,78 5,8Ы0“13 3,73-Ю”12
210 1,83 1,58-10-3 3,54-10-2 6,15 1,54-Ю-12 1,24-10""
2" 8,14 1,47-10~3 2,8Ы0~2 48,59 4,17 ■ 10-12 1 О Os
Здесь стоит отметить, что, хотя алгоритм на основе нестандартного представления работает и быстрее, чем процедура из библиотеки Ьараск, но точность полученного решения у него существенно ниже. Снижение точности происходит вследствие трешолдинга, производимого в момент вейвлет-преобразования и Ш-разложения. Суть трешолдинга состоит в том, что величины меньше некоего заданного порога считаются равными нулю. Все предыдущие примеры обрабатывались с порогом равным 1(Г5. В табл. 3 можно увидеть, как будет меняться точность и время решения, если уменьшать этот порог.
___________________________________Таблица 3
Порог Tns LX(NS) L,(NS)
10~8 43,55 1,51-Ю-8 7,75-10-8
10"9 44,38 1,29-10-9 2,13-Ю-9
10~п 48,01 6,78-10~12 7,84-10~П
Таким образом, видно, что для порога равного 1(ГП , ошибка решения становится одного порядка с ошибкой процедуры Ьараск. Что касается времени выполнения, то оно также становится сравнимым. Однако, стоит взглянуть что произойдет на следующей размерности матрицы (табл. 4).
_______________________________________________________________________Таблица 4
N TNS ^(NS) l2(ns) т Lapack Ьх (Ьараск■) Ь2(Ьараск)
212 362,87 1,31-Ю"10 5,01-10~ІО 387,36 1,46-10~п 2,03-Ю~10
Время работы алгоритма на основе нестандартного представления сопоставимо и даже несколько меньше, чем время работы итерационного алгоритма при сравнимой ошибке. Стоит отметить, что для вейвлет-преобразований, подбирая вейвлет-функцию, можно добиться сущест-венно-большего числа нулей в результате, а значит, увеличить общую скорость. Для данных примеров использовалась вейвлет-функция D4 (Добеши 4 порядка).
Итак, можно сделать вывод о том, что алгоритм решения СЛАУ на основе нестандартного представления работает сравнимо по времени выполнения и по ошибке результата с итерационными алгоритмами.
Литература
1. Добеши, И. Десять лекций по вейвлетам / И. Добеши. - Ижевск: НИЦ «Регулярная и хаотическая динамика». - 2001. - 464 с.
2. 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). - P. 156-201.
Поступила в редакцию 9 июня 2007 г.
