Участник:Бугаков Юрий/Построение матрицы Адамара произвольного размера
Содержание
- 1 Свойства и структура алгоритма
- 1.1 Общее описание алгоритма
- 1.2 Математическое описание алгоритма
- 1.3 Вычислительное ядро алгоритма
- 1.4 Макроструктура алгоритма
- 1.5 Схема реализации последовательного алгоритма
- 1.6 Последовательная сложность алгоритма
- 1.7 Информационный граф
- 1.8 Ресурс параллелизма алгоритма
- 1.9 Входные и выходные данные алгоритма
- 1.10 Свойства алгоритма
- 2 ЧАСТЬ Программная реализация алгоритма
- 2.1 Особенности реализации последовательного алгоритма
- 2.2 Локальность данных и вычислений
- 2.3 Возможные способы и особенности параллельной реализации алгоритма
- 2.4 Масштабируемость алгоритма и его реализации
- 2.5 Динамические характеристики и эффективность реализации алгоритма
- 2.6 Выводы для классов архитектур
- 2.7 Существующие реализации алгоритма
- 3 Литература
1 Свойства и структура алгоритма
1.1 Общее описание алгоритма
Важную роль в алгебре и комбинаторике играют матрицы Адамара, которые впервые были введены в математический обиход в конце прошлого века одним из крупнейших французских математиков Жаком Адамаром (1865-1963). Их применение в науке и технике посвящены тысячи публикаций. Они предоставляют эффективные возможности для организации хранения, обработки и передачи информации.
Квадратная матрица Н порядка m с элементами ±1 называется матрицей Адамара, если выполняется условие
- [math] H_m\,H_m^T = m\,E_m [/math]
Нетрудно заметить, что различные строки матрицы Адамара попарно ортогональны. Также, можно увидеть из определения, что m четно и
- [math] H_{2m} = \begin{pmatrix} H_{m} & H_{m}\\ H_{m} & -H_{m}\end{pmatrix} [/math]
С матрицами Адамара связан ряд нерешенных проблем, одна из которых состоит в следующем. Мы уже видели, что порядок m матрицы Адамара при m ≥ 3 может быть лишь четным. Более того, при m ≥ 4 порядок обязан делиться на 4. До сих пор остается открытым вопрос: для любого ли m, кратного 4, существует матрица Адамара порядка m? Неизвестно, в частности, существует ли матрица Адамара порядка 268 (это наименьший порядок, кратный 4, для которого матрица Адамара еще не построена).
Часто матрицу Адамара нормализуют и рассматривают размерности степени 2. В дальнейшем будем рассматривать только такие случаи. Переопределим матрицу Адамара следующим образом:
Матрица Адамара Hn - это матрица размера 2n × 2n, для которой справедливо равенство
- [math]H_n = H_{1} \otimes H_{n-1}[/math]
- [math]\begin{align} H_1 = \frac{1}{\sqrt2} &\begin{pmatrix}\begin{array}{rr} 1 & 1\\ 1 & -1 \end{array}\end{pmatrix} \end{align}[/math]
где [math] \otimes [/math] представляет собой тензорное произведение.
Представим некоторые частные примеры матриц Адамара:
- [math] H_0 = +1, [/math]
- [math] H_1 = \frac{1}{\sqrt2} \begin{pmatrix}\begin{array}{rr} 1 & 1\\ 1 & -1 \end{array}\end{pmatrix} [/math]
- [math] H_2 = \frac{1}{2} \begin{pmatrix}\begin{array}{rrrr} 1 & 1 & 1 & 1\\ 1 & -1 & 1 & -1\\ 1 & 1 & -1 & -1\\ 1 & -1 & -1 & 1 \end{array}\end{pmatrix} [/math]
- [math] H_3 = \frac{1}{2^{\frac{3}{2}}} \begin{pmatrix}\begin{array}{rrrrrrrr} 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1\\ 1 & -1 & 1 & -1 & 1 & -1 & 1 & -1\\ 1 & 1 & -1 & -1 & 1 & 1 & -1 & -1\\ 1 & -1 & -1 & 1 & 1 & -1 & -1 & 1\\ 1 & 1 & 1 & -1 & -1 & -1 & -1 & -1\\ 1 & -1 & 1 & 1 & -1 & 1 & -1 & 1\\ 1 & 1 & -1 & 1 & -1 & -1 & 1 & 1\\ 1 & -1 & -1 & -1 & -1 & 1 & 1 & -1\\ \end{array}\end{pmatrix} [/math]
1.2 Математическое описание алгоритма
Основой алгоритма является тот факт, что любую матрицу Адамара [math]H_n[/math], размерностью [math]2^n\times 2^n[/math] можно вычислить поэлементно [ссылочка=)]:
- [math]h_{k,l} = \frac{1}{2^\frac{n}{2}} (-1)^{\sum_j k_j l_j}[/math],
где kj и lj коэффициенты двоичного разложения чисел k и l (1 или 0)
- [math]k = \sum^{n}_{i=0} {k_i 2^i} = k_{n} 2^{n} + k_{n-1} 2^{n-1} + \cdots + k_1 2 + k_0[/math]
и
- [math]l = \sum^{n}_{i=0} {l_i 2^i} = l_{n} 2^{n} + l_{n-1} 2^{n-1} + \cdots + l_1 2 + l_0[/math].
1.3 Вычислительное ядро алгоритма
Вычислительное ядро составляют вычисления элементов матрицы [math]h_{k,l} = \frac{1}{2^\frac{n}{2}} (-1)^{\sum_j k_j l_j}[/math].
Помимо этой формулы чаще всего используется метод тензорного произведения. В связи с трудностями распараллеливания и постоянного выделения памяти при каждом шаге рекурсивного вызова этого метода, мы решили выбрать первый способ вычисления
1.4 Макроструктура алгоритма
Основную часть метода составляют вычисления значений степени -1 для каждого элемента:
- [math]\sum_{j=0}^n k_j l_j[/math].
Как видно по формуле, получаемое значение равна количеству единиц в двоичной записи результата побитового умножения k и l. Существует эффективный метод вычисления количества единиц , сложность которого равна m операций вычитания и m операций побитового умножения, где m-искомое количество единиц.
1.5 Схема реализации последовательного алгоритма
Последовательность исполнения метода следующая:
1. Побитовое умножение [math]k\&l[/math].
2. Вычисление количества единиц двоичной записи.
3. Определение знака элемента по степени.
1.6 Последовательная сложность алгоритма
Для построения матрицы Адамара порядка [math]N=2^n[/math] в последовательном варианте вначале потребуется [math]N^2[/math] побитовых умножений k\$l. Как говорилось выше, сложность метода вычисления количества единиц в числа в двоичном предчтавлении
- [math]\frac{n(n-1)}{2}[/math] делений,
- [math]\frac{n^3-n}{6}[/math] сложений (вычитаний),
- [math]\frac{n^3-n}{6}[/math] умножений.
Умножения и сложения (вычитания) составляют основную часть алгоритма.
При этом использование режима накопления требует совершения умножений и вычитаний в режиме двойной точности (или использования функции вроде DPROD в Фортране), что ещё больше увеличивает долю умножений и сложений/вычитаний во времени, требуемом для выполнения метода Холецкого.
При классификации по последовательной сложности, таким образом, метод Холецкого относится к алгоритмам с кубической сложностью.
1.7 Информационный граф
Опишем граф алгоритма[1][2][3] как аналитически, так и в виде рисунка.
Граф алгоритма состоит из трёх групп вершин, расположенных в целочисленных узлах трёх областей разной размерности.
Первая группа вершин расположена в одномерной области, соответствующая ей операция вычисляет функцию SQRT. Единственная координата каждой из вершин [math]i[/math] меняется в диапазоне от [math]1[/math] до [math]n[/math], принимая все целочисленные значения.
Аргумент этой функции
- при [math]i = 1[/math] — элемент входных данных, а именно [math]a_{11}[/math];
- при [math]i \gt 1[/math] — результат срабатывания операции, соответствующей вершине из третьей группы, с координатами [math]i - 1[/math], [math]i[/math], [math]i - 1[/math].
Результат срабатывания операции является выходным данным [math]l_{ii}[/math].
Вторая группа вершин расположена в двумерной области, соответствующая ей операция [math]a / b[/math]. Естественно введённые координаты области таковы:
- [math]i[/math] — меняется в диапазоне от [math]1[/math] до [math]n-1[/math], принимая все целочисленные значения;
- [math]j[/math] — меняется в диапазоне от [math]i+1[/math] до [math]n[/math], принимая все целочисленные значения.
Аргументы операции следующие:
- [math]a[/math]:
- при [math]i = 1[/math] — элементы входных данных, а именно [math]a_{j1}[/math];
- при [math]i \gt 1[/math] — результат срабатывания операции, соответствующей вершине из третьей группы, с координатами [math]i - 1, j, i - 1[/math];
- [math]b[/math] — результат срабатывания операции, соответствующей вершине из первой группы, с координатой [math]i[/math].
Результат срабатывания операции является выходным данным [math]l_{ji}[/math].
Третья группа вершин расположена в трёхмерной области, соответствующая ей операция [math]a - b * c[/math]. Естественно введённые координаты области таковы:
- [math]i[/math] — меняется в диапазоне от [math]2[/math] до [math]n[/math], принимая все целочисленные значения;
- [math]j[/math] — меняется в диапазоне от [math]i[/math] до [math]n[/math], принимая все целочисленные значения;
- [math]p[/math] — меняется в диапазоне от [math]1[/math] до [math]i - 1[/math], принимая все целочисленные значения.
Аргументы операции следующие:
- [math]a[/math]:
- при [math]p = 1[/math] элемент входных данных [math]a_{ji}[/math];
- при [math]p \gt 1[/math] — результат срабатывания операции, соответствующей вершине из третьей группы, с координатами [math]i, j, p - 1[/math];
- [math]b[/math] — результат срабатывания операции, соответствующей вершине из второй группы, с координатами [math]p, i[/math];
- [math]c[/math] — результат срабатывания операции, соответствующей вершине из второй группы, с координатами [math]p, j[/math];
Результат срабатывания операции является промежуточным данным алгоритма.
Описанный граф можно посмотреть на рис.1 и рис.2, выполненных для случая [math]n = 4[/math]. Здесь вершины первой группы обозначены жёлтым цветом и буквосочетанием sq, вершины второй — зелёным цветом и знаком деления, третьей — красным цветом и буквой f. Вершины, соответствующие операциям, производящим выходные данные алгоритма, выполнены более крупно. Дублирующие друг друга дуги даны как одна. На рис.1 показан граф алгоритма согласно классическому определению , на рис.2 к графу алгоритма добавлены вершины , соответствующие входным (обозначены синим цветом) и выходным (обозначены розовым цветом) данным.
1.8 Ресурс параллелизма алгоритма
Для разложения матрицы порядка [math]n[/math] методом Холецкого в параллельном варианте требуется последовательно выполнить следующие ярусы:
- [math]n[/math] ярусов с вычислением квадратного корня (единичные вычисления в каждом из ярусов),
- [math]n - 1[/math] ярус делений (в каждом из ярусов линейное количество делений, в зависимости от яруса — от [math]1[/math] до [math]n - 1[/math]),
- по [math]n - 1[/math] ярусов умножений и сложений/вычитаний (в каждом из ярусов — квадратичное количество операций, от [math]1[/math] до [math]\frac{n^2 - n}{2}[/math].
Таким образом, в параллельном варианте, в отличие от последовательного, вычисления квадратных корней и делений будут определять довольно значительную долю требуемого времени. При реализации на конкретных архитектурах наличие в отдельных ярусах ЯПФ отдельных вычислений квадратных корней может породить и другие проблемы. Например, при реализации на ПЛИСах остальные вычисления (деления и тем более умножения и сложения/вычитания) могут быть конвейеризованы, что даёт экономию и по ресурсам на программируемых платах; вычисления же квадратных корней из-за их изолированности приведут к занятию ресурсов на платах, которые будут простаивать большую часть времени. Таким образом, общая экономия в 2 раза, из-за которой метод Холецкого предпочитают в случае симметричных задач тому же методу Гаусса, в параллельном случае уже имеет место вовсе не по всем ресурсам, и главное - не по требуемому времени.
При этом использование режима накопления требует совершения умножений и вычитаний в режиме двойной точности, а в параллельном варианте это означает, что практически все промежуточные вычисления для выполнения метода Холецкого в режиме накопления должны быть двойной точности. В отличие от последовательного варианта это означает увеличение требуемой памяти почти в 2 раза.
При классификации по высоте ЯПФ, таким образом, метод Холецкого относится к алгоритмам со сложностью [math]O(n)[/math]. При классификации по ширине ЯПФ его сложность будет [math]O(n^2)[/math].
1.9 Входные и выходные данные алгоритма
Входные данные: n - размерность матрицы [math]H_n[/math].
Объём входных данных: 1 (число n).
Выходные данные: матрица Адамара [math]H_n[/math] размерностью [math]2^n\times 2^n[/math].
Объём выходных данных: [math]2^{2n}[/math].
1.10 Свойства алгоритма
Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов, как хорошо видно, является квадратичным (отношение кубической к линейной).
При этом вычислительная мощность алгоритма, как отношение числа операций к суммарному объему входных и выходных данных – всего лишь линейна.
При этом алгоритм почти полностью детерминирован, это гарантируется теоремой о единственности разложения. Использование другого порядка выполнения ассоциативных операций может привести к накоплению ошибок округления, однако это влияние в используемых вариантах алгоритма не так велико, как, скажем, отказ от использования режима накопления.
Дуги информационного графа, исходящие из вершин, соответствующих операциям квадратного корня и деления, образуют пучки т. н. рассылок линейной мощности (то есть степень исхода этих вершин и мощность работы с этими данными — линейная функция от порядка матрицы и координат этих вершин). При этом естественно наличие в этих пучках «длинных» дуг. Остальные дуги локальны.
Наиболее известной является компактная укладка графа — его проекция на треугольник матрицы, который перевычисляется укладываемыми операциями. При этом «длинные» дуги можно убрать, заменив более дальнюю пересылку комбинацией нескольких ближних (к соседям).
Эквивалентное возмущение [math]M[/math] у метода Холецкого всего вдвое больше, чем возмущение [math]\delta A[/math], вносимое в матрицу при вводе чисел в компьютер: [math] ||M||_{E} \leq 2||\delta A||_{E} [/math]
Это явление обусловлено положительной определённостью матрицы. Среди всех используемых разложений матриц это наименьшее из эквивалентных возмущений.
2 ЧАСТЬ Программная реализация алгоритма
2.1 Особенности реализации последовательного алгоритма
В простейшем варианте генерацию матрицы Адамара можно записать так:
void adamar(int** H, int n) // N = 2^n - размерность матрицы
{
int i,j,N;
float h;
h = pow(1/2,n/2);
N = pow(n,2);
for(i=0;i<N;++i)
{
for(j=0;j<N;++j)
{
ij = i & j;
for (l=0;ij;++l) ij&=ij-1; // l - количество единиц бинарного представления ij
if(l%2==0) H[i][j] = h;
else H[i][j] = -h;
}
}
}
void adamar(int** H, int N, int i, int j, int h)
{
if(n==1) H[i][j] = h;
else
{
adamar(H, N/2, i, j, h);
adamar(H, N/2, i, N/2+j, h);
adamar(H, N/2, N/2+i, j, h);
adamar(H, N/2, N/2+i, N/2+j, -h);
}
}
void adamar(int** H, int n) // N = 2^n - размерность матрицы
{
int i,j,N;
float h;
h = pow(1/2,n/2);
N = pow(n,2);
k=2;
H[0][0] = h;
while(k<=N)
{
for (i=0; i<k; ++i)
{
for (j=0; j<k; ++j)
{
if (i<k/2 && j<k/2)
{
}
else
{
if(i>=k/2 && j>=k/2)
{
H[i][j] = -(H[i-k/2][j-k/2]);
}
else
{
if (i>j)
{
H[i][j] = H[i-k/2][j];
}
else
{
H[i][j] = H[i][j-k/2];
}
}
}
}
}
k=k*2;
}
}
2.2 Локальность данных и вычислений
2.3 Возможные способы и особенности параллельной реализации алгоритма
2.4 Масштабируемость алгоритма и его реализации
2.5 Динамические характеристики и эффективность реализации алгоритма
2.6 Выводы для классов архитектур
2.7 Существующие реализации алгоритма
3 Литература
- ↑ Воеводин В.В. Математические основы параллельных вычислений// М.: Изд. Моск. ун-та, 1991. 345 с.
- ↑ Воеводин В.В., Воеводин Вл.В. Параллельные вычисления. – СПб.: БХВ - Петербург, 2002. – 608 с.
- ↑ Фролов А.В.. Принципы построения и описание языка Сигма. Препринт ОВМ АН N 236. М.: ОВМ АН СССР, 1989.