Реализация алгоритма Гивенса в ПЛИС
А. М. Сергиенко, д-р техн. наук,
Т. М. Лесик,
Национальный техн. ун-т Украины «КПИ»
(Украина, 02056, Киев, пр. Победы, 37, тел.(044) 4549337, Е-mail: aser@comsys.kpi.ua )
Вычисление матрицы Гивенса в ПЛИС
A new method of Givens rotation coefficient calculation is proposed. The method consists in finding Pythagorean triplets, which form the numerator and denominator of the sine and cosine coefficients. Such coefficient representation provide the computation error minimization due to the fact that the coefficients have precise values.
Ключевые слова: FPGA, ПЛИС, матрица Гивенса.
Введение.
Вычисление матрицы Гивенса часто применяется при решении различных задач линейной алгебры, таких как, QR-разложение матриц, нахождение собственных чисел и векторов, среднеквадратичная интерполяция [1]. При QR-разложении методом Гивенса матрица A последовательно умножается на матрицы вращения Гивенса Pi,j, в результате чего аннулируются поддиагональные элементы матрицы A = [ai,j]. Например, если необходимо аннулировать элемент a1,2 = х, стоящий под элементом a1,1 = y, то после преобразования получим: A’ = P1,2*A, где P1,2 – слабозаполненная матрица, в диагонали которой стоит ортогональная матрица P = ((с,-s)Т,(s,c)Т) и единичная матрица, причем
c = y/z , s = x/z , z =√ x2 + y2 (1)
Здесь c = cos ϕ, s = sin ϕ, где ϕ – угол поворота матрицы Гивенса. Вычисление коэффициентов (1) – сложная задача при реализации алгоритма Гивенса в специализированных процессорах и требует повышенной точности [2]. Она решается итерационными алгоритмами методом Ньютона [3], Волдера [3,4] или по алгоритму вычисления корня из суммы квадратов [5]. Итерационный алгоритм Волдера основан на элементарных поворотах вектора (y,x) и широко используется в специализированных вычислителях на базе программируемых логических интегральных схем (ПЛИС) [6,7,8]. Однако для получения п-разрядных результатов необходимо выполнить п итераций элементарных поворотов. Кроме того, требуются дополнительные итерации для коррекции увеличения длины вектора (y,x).
В докладе предлагается способ вычисления матрицы Гивенса, который обходится без извлечения квадратного корня и позволяет получить повышенную точность. Способ основан на свойствах пифагоровых троек.
Пифагоровы тройки и вычисления в компъютерах.
Решение задачи Пифагора состоит в нахождении всех прямоугольных треугольников с целочисленными сторонами, т.е. в решении диофантового уравнения:
x2 + y2 = z2 (2)
Формулы (1) и (2) показывают, что для некоторых углов ϕ′ поворачивающие коэффициенты равны рациональным дробям:
sin ϕ′ = x/z; cos ϕ′ = y/z (3)
Задача Пифагора может быть решена рядом способов [9]. В результате, коэффициенты sinϕ′, cosϕ′ представимы в компьютере точными значениями рациональных дробей (3).
В [10,11] показано, что представление поворачивающих коэффициентов с помощью пифагоровых троек (x,y,z) позволяет существенно уменьшить погрешность вычисления скользящего дискретного преобразования Фурье и делает эффективным реализацию этого алгоритма в ПЛИС. Эти коэффициенты вычисляются заранее и хранятся в постоянной памяти.
Пифагоровы тройки в матрице Гивенса.
Умножение на матрицу Гивенса означает вращение двумерного вектора:
aʼi,j = (ai,j⋅у + ai,k⋅х)/z ; aʼi,k = (ai,k⋅у – ai,j⋅х)/z (4)
где вектор (ai,j,ai,k) составлен из пар элементов матрицы А. Применение поворачивающих коэффициентов (3) в матрице Гивенса позволяет сохранить длину вектора (aʼi,j,aʼi,k) после поворота в силу равенства (2). Поэтому операция поворота (4) выполняется с минимальной погрешностью.
Эта погрешность может быть минимизирована при использовании арифметики целых чисел или рациональных дробей достаточной розрядности. В случае целочисленного представления данных деление на z относится ко всем элементам пар строк матрицы А. Поэтому оно может не выполняться. Тогда число 1/z следует запомнить как масштабный коэффициент, который учитывается при формировании окончательного результата.
Вычисление поворачивающих коэффициентов.
Поиск пифагоровых троек с требуемыми свойствами представляет решение диофантова уравнения, которое ищется с помощью переборного алгоритма [9]. Поэтому такой алгоритм из-за своей сложности не годится для формирования матрицы Гивенса. Тем не менее, для малых углов ϕ можно предложить простой поиск пифагоровых троек.
Пусть ai,j > 0, ai,k > 0 , ai,j >> ai,k. Тогда
m = ]ai,j / ai,k[;
x = 4m ; y = 4m2 – 1; z = 4m2 + 1 (5)
где ]*[ – число, округленное до целого, у/х ≈ m.
Для ai,j < р⋅ai,k целесообразно использовать табличный поиск коэффициентов, где p равно нескольким сотням или тысячам. Такая таблица имеет р входов. Умножение на табличные коэффициенты позволяет после поворота (4) получить вектор (aʼi,j,aʼi,k), в котором aʼi,j > p⋅|aʼi,k| . Если этого недостаточно для обнуления ai,k, то для вектора (aʼi,j,aʼi,k) следует найти коэффициенты по формулам (5). В результате получим две пифагоровы тройки: (x1,y1,z1), (x2,y2,z2). Результирующая тройка равна
у = y1y2 – x1x2 ; х = x1y2 + y1x2 ; z = z1z2 (6)
Матрица Гивенса, сформированная такой тройкой, обеспечивает уменьшение второго элемента вектора (ai,j,ai,k) более чем в р2 раз. Такое уточнение коэффициентов можно повторить, причем после каждой итерации уточнения ошибка обнуления ai,k уменьшается более чем в р раз.
Пусть a1,1 = 100000, a1,2 = 57735, т.е. вектор (a1,1,a1,2) имеет угол 30°. Отношению a1,2/a1,1 ≈ 74/128 соответствует табличное значение пифагоровой тройки (120,209,241). После вычислений по формулам (4) получаем aʼ1,1 = 27828200/241 ; aʼ1,2 = 66615/241 .
Согласно (5), m = 417; вторая тройка равна (1668,695555,695557). Результирующая тройка – (83815212,145170835,167629237). Синус и косинус, которые ей
отвечают по (3), имеют погрешность 3,5⋅10–6 , что меньше относительной погрешности представления исходных данных. При использовании коэффициентов этой тройки как коэффициентов Гивенса (1) после поворота вектора (a1,1,a1,2) получим a1,2 = 0,5.
Недостаток представления матрицы Гивенса с помощью пифагоровой тройки состоит в том, что целые коэффициенты x, y, z могут иметь большую разрядность, что затрудняет их практическое использование. Тем не менее, их можно усечь до n старших разрядов, где n – разрядность мантиссы вычислителя. Так, после деления на 103 коэффициентов тройки, найденной в приведенном примере, погрешность вычислений увеличивается незначительно. В этом случае можно также упростить поиск коэффициентов по (5), т.е. x = 1; y = z = m. Но тогда это не будет пифагоровой тройкой.
Вычисление матрицы Гивенса в ПЛИС.
Было разработано экспериментальное устройство для вычисления коэффициентов x, y, z для 18-разрядных целых чисел ai,j, ai,k, которые обеспечивают аннулирование числа ai,k. Устройство работает в конвейерном режиме с глубиной конвейера, равной десяти. Оно занимает в ПЛИС Xilinx Virtex4 447 CLB slices, 3 BlockRAM, 4 блока умножения DSP48 и работает на тактовой частоте до 145 МГц.
Выводы.
Предложенный способ формирования матрицы Гивенса позволяет выполнять вычисления без использования извлечения квадратного корня. В отличие от алгоритма Волдера, способ имеет в log2р раз меньшее число итераций, где р – объем таблицы начальных приближений коэффициентов. Способ имеет эффективную реализацию в ПЛИС.
1. Ортега Дж. Введение в параллельные и векторные методы решения линейных систем.-М.: Мир, 1991.-366 с.
2. Brisebarre N., Jodes M., Kornerup P. Augmented precision square roots, 2-D norms, and discussion on correctly rounding √x2 + y2 // Proc. of the IEEE Symp. on Computer Arithmetic. Tuebingen, Germany. –2011.
3. Сергиенко А.М. Отображение алгоритма QR-разложения Гивенса в многопроцессорную систему // Электронное моделирование. −2004. –Т.26, №5. − С. 43−53.
4. Ahmed H.M., Desolme J.M. and Morf M. Higly Concurrent computing structures for matrix arithmetic and signal processing // Computer,-1982.-V. 15.-p. 65-82.
5. Moler C., Morrison D. Replacing Square Roots by Pythagorean Sums // IBM J. Res. Develop. –1983. –V27. – N6. – P. 577-581.
6. Sergyienko А., Maslennikov О. Implementation of Givens QR Decomposition in FPGA // Lecture Notes in Computer Science. − Berlin: Springer. −2002. −V. 2328. −Р. 453−459.
7. Cavallaro J.R., Luk, F.T. CORDIC arithmetic for an SVD processor. // Journal of Parallel and Distributed Computing. –1998. –N5. –P. 271-290.
8. Liu Y., Bouganis С.-S., Cheung P. Y. K. Hardware Efficient Architectures for Eigenvalue Computation // Proc. DATE’06. – 2006. – P.953-958.
9. Koshy Т. Elementary Number Theory with Applications. Academic Press. –2007. –773 p.
10.Сергієнко А.М., Мєлковська В.М. Спосіб демодуляції сигналів з багаточастотною модуляцією //Вісник НТУУ “КПІ”: “Інформатика, управління та обчислювальна техніка”: зб. наук. праць. –2008, –№48. –С. 82−84.
11. Сергиенко А.М., Лесик Т.М. Вычисление линейных рекуррентных последовательностей в ПЛИС // Тр. междунар. конф. –Киев. –2011. –С. 135-138.