АЛГОРИТМ СИНГУЛЯРНОГО РАЗЛОЖЕНИЯ МАТРИЦ Текст научной статьи по специальности «Компьютерные и информационные науки»

УДК 621.3

АЛГОРИТМ СИНГУЛЯРНОГО РАЗЛОЖЕНИЯ МАТРИЦ

Сафонова Анастасия Владимировна

аспирант кафедры радиотехнических систем. ФГБОУ ВПО «Рязанский государственный радиотехнический университет».

E-mail: anastasia-altair@mail.ru. Адрес: 390005, г. Рязань, ул. Гагарина, д. 59/1. Аннотация: Предложен алгоритм сингулярного разложения матриц, имеющий быструю сходимость. Представлена схема работы алгоритма на устройствах с параллельной обработкой информации. Получен график зависимости недиагональной нормы от числа требуемых итераций. Проведено сравнение предлагаемого алгоритма с алгоритмом BLV с точки зрения числа повторений, требующихся для сокращения недиагональной нормы. Показано преимущество предложенного алгоритма перед алгоритмом BLV при работе с матрицами размерностью больше 32.

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

Введение

Большинство современных алгоритмов сверхразрешения требуют разложения кросс-спектральной матрицы по собственным векторам или сингулярного разложения матрицы полученных данных [1]. Кроме того, сингулярное разложение матриц является важным компонентом многих других алгоритмов обработки сигнала и используется в алгоритмах обработки изображений [2], а также для оценки каналов в MIMO системах [3-4]. Таким образом, для практического применения требуется проведение сингулярного разложения в реальном времени [5]. Для антенных решеток (АР), имеющих малые размеры, удовлетворить данному требованию несложно. Однако, если используются большие матрицы, то достижение быстрой сходимости является трудной задачей

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

На данный момент разработано множество алгоритмов сингулярного разложения матриц

[7]. Один из наиболее известных методов, обеспечивающий сходимость матриц при сингулярном разложении, - метод Якоби [8]. Этот метод может быть использован только для симметричных матриц и имеет квадратичную

сходимость (сходимость порядка 2). Кроме того, доказано, что он является более точным, чем методы на QR основе [9].

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

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

Метод Якоби использует последовательность плоских вращений для диагонализации матрицы Л. Для сингулярного разложения используется двустороннее плоское вращение Якоби [10].

Л;+1 = JíГ Л, J г. (1)

Такое преобразование называется вращением, а матрица J - матрицей вращения Якоби. Преобразование (1) обладает следующими свойствами[11]:

1. Свойство симметрии: ац = а^ ;

2. Элементы стоящие на главной диагонали, остаются неизменны. Изменяются только элементы, стоящие в ь ой и j- ой строке и в j- м и ь м столбце;

3. Собственные векторы и собственные значения сохраняются.

Матрица Якоби - это единичная матрица, четыре элемента которой с индексами (р,р),

(р,ц), (ц,р), (ц,ц) можно заменить следующими значениями: 3рр = СО8(0), 3рц = sm(в),

Зцр = ~sin(0), Лцц = cos(в). Косинус и синус

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

pp

0

qq.

ci si

~si ci J p L aqp

a

a

pp pq

aq

— s

qq _|L r

q

aij a(i+1)j

У

a(i+1)j

ai(j +1) a(i +1)j

a»■ ai(j+1) a(i+1)j

cl sl —sl cl

r

—sr

где С1 = cos(вl), сг = cos(вr), Sl = sm(в^), sr = sin(вr), а'рр, аЦц - диагональные элементы, полученные после диагонализации. Замечено, что левостороннее вращение Якоби влияет только на элементы в столбцах р и ц, а правостороннее вращение Якоби - только на элементы в строках р и ц. Таким образом, в матрице размером N N/2 подзадач могут решаться параллельно.

Форсит и Хенриси расширили метод Якоби для работы с обычными матрицами [12]. Позже цикличная версия их подхода была реализована в виде Вгеп^ик-Уап Loan(BLV) алгоритма для систолических антенных решеток [13]. Алгоритм BLV использует двустороннюю трансформацию, что доказывает сохранение квадратичной сходимости метода Якоби [14]. Хестен предложил одностороннюю вариацию метода Якоби, но представление сингулярных векторов при этом было не таким точным, как в двустороннем случае [15].

Преимуществом алгоритма ВЦУ для небольших АР является простая структура и невысокие вычислительные затраты, но для сходимости матриц больших размеров требуется огромное число итераций [16]. В предлагаемом методе проблема медленной сходимости реша-

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

Предлагаемый алгоритм заключается в следующем:

1. Проводим преобразование диагональных элементов матрицы.

1.1 Находим N/2 больших элементов диагонали матрицы (хI;у{), 7 = 1,...,N / 2 , таких, что: Х+1 * xl,...,7, х+1 * уь...7, у+1 * xl,...,7,

Уг+1 * Л —7, X * У .

1.2. Формируем подматрицы Адиаг размерностью 2 х 2 :

а хх Ах

а а

_ УЛ УУ.

1.3. К полученным матрицам Адиаг применяем сингулярное разложение и находим левые и правые параметры вращения в[ и вг соответственно.

1.4. Применяем двустороннее плоское вращение Якоби:

a

диаг

T

T

c

s

r

r

0

c

r

X

s

г

X

c

г

Л диаг = Т1 Лдиаг Jг

2. Проводим преобразование столбцов матрицы.

2.1.Формируем подматрицу Лст размерностью 2 х 2 . Если хI < У1, то

А. с х Лс 2

Лс

Л =

Лст

Ч ¡У,

Л

С2гУг

7 = 1,...,( N /2 -1),

Ск * X, у , к = 1,...,( N - 2),

или

Л =

ст

Л

С2 мУ

Л

С2 У

Л Л

С2 ]У1

2.2. Применяем правостороннее вращение к подматрице столбцов:

л' = Л т

^ст ^ст^г

2.3. Заменяем подматрицу в исходной матрице.

3. Проводим преобразование строк матрицы.

3.1. Формируем подматрицу Астр размерностью 2 х 2 . Если х7 < у7, то

" Л Л

х/2 ¡-1 х

Л1? ~ Л,

Л

стр

У/2,

УГ,

или

Л

стр

ЛЛ

УГ УГ ]

ЛЛ

ХГ2 ¡-1 х?2

где Гк * х7,у7, к = 1,...,(N - 2).

3.2. Применяем левостороннее вращение к подматрице строк:

л' = т,л

стр " I стр

3.3. Заменяем подматрицу в исходной матрице.

3.4. Выводим полученную информацию и ждем новые входные данные.

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

Входные данные

НБЭ — Л ФПДЭ

ЭП

■ !—!■

ПДЭ

-1

ФПНДЭ

^ ПСТ |,

ЭП

[_

ПСТ

ф ПСТ

114

N/2

ФПДЭ

ЭП

ПДЭ

"1-

ПСТ

ФПНДЭ

ЭП

ПСТ

ПСТ

МК

-I

ПСТР

ПСТР

Ф N/2

ПСТР

ПСТР

Выходные данные

Контроллер

г

Рис. 1. Схема реализации предложенного алгоритма

1

1

N/2-1

N/2-1

N/2-1

N/2-1

но велика [17], поэтому алгоритм рассчитан на работу именно на таких микросхемах.

Реализация предложенного алгоритма на

устройствах с параллельной обработкой информации

Схема реализации предлагаемого алгоритма представлена на рис. 1.

На рис. 1 приняты следующие обозначения: НБЭ - блок нахождения большего элемента, ФПДЭ - блок формирования подматрицы диагональных элементов, ФПНДЭ - блок формирования подматрицы недиагональных элементов, ЭП - элемент памяти, ПДЭ - блок преобразования матрицы диагональных элементов, ПСТ - блок преобразования столбцов матрицы, ПСТР - блок преобразования строк матрицы, МК - матричный коммутатор.

Закрашенные стрелки показывают связь параметров вращения, а белые - связь матричных элементов. Так как предложенный метод является итерационным, необходимо продумать критерий остановки вычислений. Заранее рассчитывается допустимое отклонение, и, если полученный результат удовлетворяет исходному условию, начинается следующая итерация. Эта функция включена в контроллер. Контроллер также распределяет полученную информацию между входными и выходными банками памяти, если итерации еще не завершены. Предложенная схема также предназначена для увеличения пропускной способности, так что контроллер должен быть разработан соответствующим образом. Блок НБЭ находит N/2 больших элементов, как это указано в шаге 1.1 алгоритма.

Основными модулями являются блоки формирования подматриц диагональных и недиагональных элементов, а также блоки, производящие преобразование столбцов и строк матриц и ПДЭ. Блок ФПДЭ проводит формирование подматрицы Адиаг . Блок ФПНДЭ отвечает за формирование подматриц Лст и Астр, описанных в п.2.1 и п. 3.1 предложенного алгоритма. Блок ПДЭ рассчитывает левые и правые параметры вращения и передает их в модули преобразования строк и столбцов. Операции, производимые в данном блоке, соответствуют пп. 1.3 -1.4 алгоритма. Более подробно принцип работы ПДЭ описан в [18]. В блоках ПСТ и ПСТР рассчитанные параметры

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

Результаты моделирования

В результате моделирования было проведено сравнение предложенного алгоритма с алгоритмом ВЦУ с точки зрения числа повторений, требующихся для сокращения недиагональной нормы (суммы квадратов всех недиагональных элементов) с заданной точностью (в данном эксперименте точность - 10-15). Под числом повторений подразумевается количество итераций, деленное на N где N - размерность матрицы N х N. Моделирование проводилось для матриц с размерностями от 4 х 4 до 128х128.

На рис. 2 показана зависимость числа повторений от размерности используемой матрицы для ВЦУ и для предложенного алгоритма.

—в— BLV

• Предложенный алгоритм

JS

1

Размерность матрицы (N)

Рис. 2. Зависимость числа повторений от размерности используемой матрицы

Было получено, что число повторений, необходимых для реализации предложенного алгоритма, при использовании матриц больших размеров примерно одинаково. Это связано с тем, что при использовании матриц с N > 32 число итераций практически пропорционально размерности матрицы. Поэтому число повторений сокращается с увеличением размера матрицы, хотя данное свойство является преимуществом только при работе с большими матрицами.

На рис. 3 представлен график зависимости недиагональной нормы от числа итераций для предложенного алгоритма (ПА) и алгоритма ВЦУ.

... -----BLV N=128 t » ПА N=128 _.~~BLV N=64

§ \ 4. 1 \ \ ч Ъ \ 4 \ -М-ПА N=64 _ _ BLV N=8 — tt\ N=8

\ 4 \ \ \ \ \ 4 4, 4.

i \ \\ у 4 \ V 4

i

О 200 400 600 800 1000

Число итераций

Рис. 3. Зависимость недиагональной нормы от _числа итераций_

По полученным данным можно сделать вывод, что при использовании матриц малой размерности (N=8) оба алгоритма обладают примерно одинаковой сходимостью. Однако, при работе с большими матрицами (N=64, N=128) число требуемых итераций при использовании предложенного алгоритма сокращается примерно в 2 раза по сравнению с применением BLV, что является существенным при реализации алгоритмов.

Литература

1. Кошелев В.И., Сафонова А.В. Модифицированный Propagator метод оценки направления прихода радиосигнала // Вестник Рязанского государственного радиотехнического университета. -2014.- №47. - С. 53-58.

2. M. Rahmati, M. S. Sadri, and M. A. Naeini, "FPGA based singular value decomposition for image processing applications," in Application-Specific Systems, Architectures and Processors, 2008. ASAP 2008. International Conference on. Pp. 185-190.

3. Матричные разложения для формирования и обработки сигналов в MIMO системах связи в каналах с МСИ / Ю. Б. Нечаев, А. А. Малютин // Радиотехника. - 2012. - № 8. - С. 40-44 .

4. Y.L. Chen, C.Z. Zhan, T.J. Jheng, and A.Y.Wu, "Reconfigurable adaptive singular value decomposition engine design for high-throughput MIMO-OFDM systems," Very Large Scale Integration (VLSI) Systems, IEEE Transactions on, vol. 21, no. 4, pp. 747-760, 2013.

5.Сафонова А.В. Эффективность алгоритма оце-

Поступила 25 декабря 2015 г.

нивания угловых координат источника радиосигнала при различных методах обработки входных реализаций // Радиотехнические и телекоммуникационные системы. - 2015. - №2. - С. 54-60.

6. Вержбицкий В. М. Численные методы: Линейная алгебра и нелинейные уравнения - М.: Высшая школа, 2000. - 266 с.

7. Van der Vorst H.A., van Dooren P. Parallel algorithms for numerical linear algebra - Amsterdam: North-Holland, 1990. - 320 pp.

8. Основы матричных вычислений / Д.Уоткинс; Пер. с англ. - М.: БИНОМ. Лаборатория знаний, 2006. - 664 с.

9. J. Demmel and V. K., "Jacobi's method is more accurate than QR," SIAM Journal on Matrix Analysis and Applications, vol. 13, no. 4, pp. 1204-1245, 1992.

10. Ширапов Д.Ш. Численные методы линейной алгебры: Учебное пособие. - Улан-Удэ: Изд-во ВСГТУ, 2003. - 96 с.

11. Парлетт Б. Симметричная проблема собственных значений. Численные методы: Пер. с англ. - М.:Мир, 1983, 384 с.

12. G. Forsythe and P. Henrici, The cyclic jacobi method for computing the principal values of a complex matrix, Trans. Amer. Math. Soc., 94, 1958, pp. 1-23.

13. F. T. L. R. P. Brent and C. V. Loan, "Computation of the singular value decomposition using mesh-connected processors," Numerische Mathematik, 1985, pp. 479-491.

14. F. T. Luk and P. Haesun, "A proof of convergence for two parallel Jacobi SVD algorithms", Computers, IEEE Transactions on, vol. 38, no. 6, 1989, pp. 806-811.

15. M. R. HESTENES, Inversion of matrices by biorthogonalization and related results, J. Soc. Indust. Appl. Math., vol.6, 1958, pp. 51-90.

16. N. D. Hemkumar and J. R. Cavallaro, "A systolic VLSI architecture for complex SVD," in Circuits and Systems, 1992. ISCAS '92. Proceedings., 1992 IEEE International Symposium on, vol. 3, pp. 1061-1064.

17. Мелехин В.Ф., Павловский Е.Г. Вычислительные машины: учебник для студ. учреждений высш. проф. образования - М.: Издательский центр Академия, 2013. - 384 с.

18. A. Ahmedsaid, A. Amira, and A. Bouridane, "Improved SVD systolic array and implementation on FPGA," in Field-Programmable Technology (FPT), 2003. Proceedings. 2003 IEEE International Conference on, pp. 35-42.

English

Algorithm of matrix singular value decomposition

Anastasia Vladimirovna Safonova - Graduate student Radio Engineering Systems Department "The Ryazan state radio engineering university". E-mail: anastasia-altair@mail.ru.

Address: 390005, Ryazan, Gagarin str., 59/1.

Abstract: Matrix singular value decomposition is an important component of signal processing algorithms and is used in majority known super-resolution algorithms, in image processing algorithms, and also for channels assessment in MIMO systems. Applied application of these algorithms requires implementation of singular value decomposition in real time. This requirement is easy to meet for small-size antenna arrays. However it is a difficult task to achieve fast convergence if big matrixes are used. The article introduces algorithm of matrix singular value decomposition having fast convergence and also its operation scheme in devices with parallel information processing is presented. As a result of modeling, comparison was made between proposed algorithm and BLV algorithm from the point of view of iterations number required to reduce off-diagonal norm with specified accuracy. Dependency diagram of off-diagonal norm on required iterations number is obtained. Dependency of iterations number on applied matrix size for BLV and for proposed algorithm is shown. Advantage of proposed algorithm over BLV algorithm is manifested when applying matrixes of size more than 32. It was proved that iterations number required for proposed algorithm implementation when using big sized matrixes is approximately the same. It is due to the fact that iterations number is almost proportional to matrix size when using such matrixes. Thus, use of proposed algorithm with big sized matrixes application enables to reduce twice required iterations number compared to BLV algorithm and that is tangible benefit in algorithms implementation.

Key words: singular value decomposition, off-diagonal norm, BLV algorithm, SVD fast convergence.

References

1. Koshelev V. I., Safonova A.V. Modified Propagator estimation method of incoming radio signal direction. -Vestnik Ryazanskogo gosudarstvennogo radiotekhnicheskogo universiteta. - 2014. - No. 47. - P. 53-58.

2. M. Rahmati, M. S. Sadri, and M. A. Naeini FPGA based singular value decomposition for image processing applications in Application-Specific Systems, Architectures and Processors, 2008. ASAP 2008. International Conference on, pp. 185-190.

3. Matrix decompositions for generating and processing signals in MIMO communication systems with MSI channels. - B. Nechayev, A. A. Malyutin// Radiotekhnika. - 2012. - No. 8. - P. 40-44.

4. Y.L. Chen, C.Z. Zhan, T.J. Jheng, and A. Y.Wu Reconfigurable adaptive singular value decomposition engine design for high-throughput MIMO-OFDM systems, Very Large Scale Integration (VLSI) Systems, IEEE Transactions on, vol. 21, no. 4, pp. 747-760, 2013.

5. Safonova A.V. Algorithm efficiency for estimation of source radio signal angular coordinates with different methods of processing of input implementations. - Radiotekhnicheskiye i telekommunikatsionnye sistemy. - 2015. -No. 2. - P. 54-60.

6. Verzhbitsky V. M. Numerical methods: Linear algebra and nonlinear equations - M.: Vysshaya shkola, 2000 -266 p.

7. Van der Vorst H.A., van Dooren P. Parallel algorithms for numerical linear algebra - Amsterdam: North-Holland, 1990 - 320 pp.

8. Fundamentals of matrix computations. - D. Watkins; Transl. from English - M.: BINOM. Laboratory of knowledge, 2006. - 664 p.

9. J. Demmel and V. K. Jacobi's method is more accurate than QR, SIAM Journal on Matrix Analysis and Applications, vol. 13, no. 4, pp. 1204-1245, 1992.

10. Shirapov D.Sh. Numerical methods of linear algebra: Manual. - Ulan-Ude: Publ. h. of VSGTU, 2003. - 96 p.

11. Parlett B. Symmetrical Problem of Eigen Values. Numerical methods: Transl. from Engl. - M.: Mir, 1983, 384 pages.

12. G. Forsythe and P. Henrici The cyclic jacobi method for computing the principal values of a complex matrix, Trans. Amer. Math. Soc., 94, 1958, pp. 1-23.

13. F.T.L.R.P. Brent and C.V. Loan Computation of the singular value decomposition using mesh-connected processors Numerische Mathematik, 1985, pp. 479-491.

14. F.T. Luk and P. Haesun A proof of convergence for two parallel Jacobi SVD algorithms Computers, IEEE Transactions on, vol. 38, no. 6, 1989, pp. 806-811.

15. M.R. HESTENES, Inversion of matrices by biorthogonalization and related results, J. Soc. Indust. Appl. Math., vol.6, 1958, pp. 51-90.

16. N.D. Hemkumar and J.R. Cavallaro A systolic VLSI architecture for complex SVD in Circuits and Systems, 1992. ISCAS '92. Proceedings, 1992 IEEE International Symposium on, vol. 3, pp. 1061-1064.

17. Melekhin V.F., Pavlovsky E.G. Computers: Textbook for students of Higher Professional Training Institutions - M.: Publishing center Akademiya, 2013. - 384 p.

18. A. Ahmedsaid, A. Amira, and A. Bouridane Improved SVD systolic array and implementation on FPGA in Field-Programmable Technology (FPT), 2003. Proceedings. 2003 IEEE International Conference on, pp. 35-42.