Участник:Бугаков Юрий/Построение матрицы Адамара произвольного размера: различия между версиями

Материал из Алговики
Перейти к навигации Перейти к поиску
← Предыдущая правкаСледующая правка →
Строка 5: Строка 5:




Важную роль в алгебре и комбинаторике играют матрицы Адамара, которые впервые были введены в математический обиход в конце прошлого века одним из крупнейших французских математиков Жаком Адамаром (1865-1963). Их применение в науке и технике посвящены тысячи публикаций. Они предоставляют эффективные возможности для организации хранения, обработки и передачи информации.
Важную роль в алгебре и комбинаторике играют матрицы Адамара, которые впервые были введены в математический обиход в конце прошлого века одним из крупнейших французских математиков [https://ru.wikipedia.org/wiki/%D0%90%D0%B4%D0%B0%D0%BC%D0%B0%D1%80,_%D0%96%D0%B0%D0%BA Жаком Адамаром(1865-1963)]. Их применение в науке и технике посвящены тысячи публикаций. Они предоставляют эффективные возможности для организации хранения, обработки и передачи информации.


Квадратная матрица Н порядка ''m'' с элементами ±1 называется матрицей Адамара, если выполняется условие
Квадратная матрица Н порядка ''m'' с элементами ±1 называется матрицей Адамара, если выполняется условие

Версия от 14:15, 14 октября 2016

Свойства и структура алгоритма

Общее описание алгоритма

Важную роль в алгебре и комбинаторике играют матрицы Адамара, которые впервые были введены в математический обиход в конце прошлого века одним из крупнейших французских математиков Жаком Адамаром(1865-1963). Их применение в науке и технике посвящены тысячи публикаций. Они предоставляют эффективные возможности для организации хранения, обработки и передачи информации.

Квадратная матрица Н порядка m с элементами ±1 называется матрицей Адамара, если выполняется условие

[math]\displaystyle{ H_m\,H_m^T = m\,E_m }[/math]

Нетрудно заметить, что различные строки матрицы Адамара попарно ортогональны. Также, можно увидеть из определения, что m четно и

[math]\displaystyle{ 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]\displaystyle{ H_n = H_{1} \otimes H_{n-1} }[/math]
[math]\displaystyle{ \begin{align} H_1 = \frac{1}{\sqrt2} &\begin{pmatrix}\begin{array}{rr} 1 & 1\\ 1 & -1 \end{array}\end{pmatrix} \end{align} }[/math]

где [math]\displaystyle{ \otimes }[/math] представляет собой тензорное произведение.

Представим некоторые частные примеры матриц Адамара:

[math]\displaystyle{ H_0 = +1, }[/math]
[math]\displaystyle{ H_1 = \frac{1}{\sqrt2} \begin{pmatrix}\begin{array}{rr} 1 & 1\\ 1 & -1 \end{array}\end{pmatrix} }[/math]
[math]\displaystyle{ 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]\displaystyle{ 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]

Математическое описание алгоритма

Основой алгоритма является тот факт, что любую матрицу Адамара [math]\displaystyle{ H_n }[/math], размерностью [math]\displaystyle{ 2^n\times 2^n }[/math] можно вычислить поэлементно [ссылочка=)]:

[math]\displaystyle{ h_{k,l} = \frac{1}{2^\frac{n}{2}} (-1)^{\sum_j k_j l_j} }[/math],

где kj и lj коэффициенты двоичного разложения чисел k и l (1 или 0)

[math]\displaystyle{ 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]\displaystyle{ 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].

Вычислительное ядро алгоритма

Вычислительное ядро составляют вычисления элементов матрицы [math]\displaystyle{ h_{k,l} = \frac{1}{2^\frac{n}{2}} (-1)^{\sum_j k_j l_j} }[/math].

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

Макроструктура алгоритма

Основную часть метода составляют вычисления значений степени -1 для каждого элемента:

[math]\displaystyle{ \sum_{j=0}^n k_j l_j }[/math].

Как видно по формуле, получаемое значение равна количеству единиц в двоичной записи результата побитового умножения k и l. Существует эффективный метод вычисления количества единиц , сложность которого равна m операций вычитания и m операций побитового умножения, где m-искомое количество единиц.

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

Последовательность исполнения метода следующая:

1. Побитовое умножение [math]\displaystyle{ k\&l }[/math].

2. Вычисление количества единиц двоичной записи.

3. Определение знака элемента по степени.

Последовательная сложность алгоритма

Для построения матрицы Адамара порядка [math]\displaystyle{ N=2^n }[/math] в первом шаге потребуется [math]\displaystyle{ N^2 }[/math] побитовых умножений k&l. Во втором шаге, как говорилось выше, сложность метода вычисления количества единиц числа в двоичном представлении равна удвоенному количеству единиц. В наихудшем случае сложность равна [math]\displaystyle{ 2log_2N }[/math], при k=l=N-1, а в наилучшем 0. Для нахождения степени всех элементов матрицы потребуются [math]\displaystyle{ 2N^2[\frac{log_2N}{4}] }[/math] операций. В третьем шаге для определения знака потребуется [math]\displaystyle{ N^2 }[/math] операций.

Информационный граф

Опишем граф алгоритма[1][2][3] как аналитически, так и в виде рисунка.

Граф алгоритма состоит из трёх групп вершин, расположенных в целочисленных узлах трёх областей разной размерности.

Первая группа вершин расположена в одномерной области, соответствующая ей операция вычисляет функцию SQRT. Единственная координата каждой из вершин [math]\displaystyle{ i }[/math] меняется в диапазоне от [math]\displaystyle{ 1 }[/math] до [math]\displaystyle{ n }[/math], принимая все целочисленные значения.

Аргумент этой функции

  • при [math]\displaystyle{ i = 1 }[/math] — элемент входных данных, а именно [math]\displaystyle{ a_{11} }[/math];
  • при [math]\displaystyle{ i \gt 1 }[/math] — результат срабатывания операции, соответствующей вершине из третьей группы, с координатами [math]\displaystyle{ i - 1 }[/math], [math]\displaystyle{ i }[/math], [math]\displaystyle{ i - 1 }[/math].

Результат срабатывания операции является выходным данным [math]\displaystyle{ l_{ii} }[/math].

Вторая группа вершин расположена в двумерной области, соответствующая ей операция [math]\displaystyle{ a / b }[/math]. Естественно введённые координаты области таковы:

  • [math]\displaystyle{ i }[/math] — меняется в диапазоне от [math]\displaystyle{ 1 }[/math] до [math]\displaystyle{ n-1 }[/math], принимая все целочисленные значения;
  • [math]\displaystyle{ j }[/math] — меняется в диапазоне от [math]\displaystyle{ i+1 }[/math] до [math]\displaystyle{ n }[/math], принимая все целочисленные значения.

Аргументы операции следующие:

  • [math]\displaystyle{ a }[/math]:
    • при [math]\displaystyle{ i = 1 }[/math] — элементы входных данных, а именно [math]\displaystyle{ a_{j1} }[/math];
    • при [math]\displaystyle{ i \gt 1 }[/math] — результат срабатывания операции, соответствующей вершине из третьей группы, с координатами [math]\displaystyle{ i - 1, j, i - 1 }[/math];
  • [math]\displaystyle{ b }[/math] — результат срабатывания операции, соответствующей вершине из первой группы, с координатой [math]\displaystyle{ i }[/math].

Результат срабатывания операции является выходным данным [math]\displaystyle{ l_{ji} }[/math].

Третья группа вершин расположена в трёхмерной области, соответствующая ей операция [math]\displaystyle{ a - b * c }[/math]. Естественно введённые координаты области таковы:

  • [math]\displaystyle{ i }[/math] — меняется в диапазоне от [math]\displaystyle{ 2 }[/math] до [math]\displaystyle{ n }[/math], принимая все целочисленные значения;
  • [math]\displaystyle{ j }[/math] — меняется в диапазоне от [math]\displaystyle{ i }[/math] до [math]\displaystyle{ n }[/math], принимая все целочисленные значения;
  • [math]\displaystyle{ p }[/math] — меняется в диапазоне от [math]\displaystyle{ 1 }[/math] до [math]\displaystyle{ i - 1 }[/math], принимая все целочисленные значения.

Аргументы операции следующие:

  • [math]\displaystyle{ a }[/math]:
    • при [math]\displaystyle{ p = 1 }[/math] элемент входных данных [math]\displaystyle{ a_{ji} }[/math];
    • при [math]\displaystyle{ p \gt 1 }[/math] — результат срабатывания операции, соответствующей вершине из третьей группы, с координатами [math]\displaystyle{ i, j, p - 1 }[/math];
  • [math]\displaystyle{ b }[/math] — результат срабатывания операции, соответствующей вершине из второй группы, с координатами [math]\displaystyle{ p, i }[/math];
  • [math]\displaystyle{ c }[/math] — результат срабатывания операции, соответствующей вершине из второй группы, с координатами [math]\displaystyle{ p, j }[/math];

Результат срабатывания операции является промежуточным данным алгоритма.

Описанный граф можно посмотреть на рис.1 и рис.2, выполненных для случая [math]\displaystyle{ n = 4 }[/math]. Здесь вершины первой группы обозначены жёлтым цветом и буквосочетанием sq, вершины второй — зелёным цветом и знаком деления, третьей — красным цветом и буквой f. Вершины, соответствующие операциям, производящим выходные данные алгоритма, выполнены более крупно. Дублирующие друг друга дуги даны как одна. На рис.1 показан граф алгоритма согласно классическому определению , на рис.2 к графу алгоритма добавлены вершины , соответствующие входным (обозначены синим цветом) и выходным (обозначены розовым цветом) данным.

Рисунок 1. Граф алгоритма без отображения входных и выходных данных. SQ - вычисление квадратного корня, F - операция a-bc, Div - деление.
Рисунок 2. Граф алгоритма с отображением входных и выходных данных. SQ - вычисление квадратного корня, F - операция a-bc, Div - деление, In - входные данные, Out - результаты.

Ресурс параллелизма алгоритма

Для разложения матрицы порядка [math]\displaystyle{ n }[/math] методом Холецкого в параллельном варианте требуется последовательно выполнить следующие ярусы:

  • [math]\displaystyle{ n }[/math] ярусов с вычислением квадратного корня (единичные вычисления в каждом из ярусов),
  • [math]\displaystyle{ n - 1 }[/math] ярус делений (в каждом из ярусов линейное количество делений, в зависимости от яруса — от [math]\displaystyle{ 1 }[/math] до [math]\displaystyle{ n - 1 }[/math]),
  • по [math]\displaystyle{ n - 1 }[/math] ярусов умножений и сложений/вычитаний (в каждом из ярусов — квадратичное количество операций, от [math]\displaystyle{ 1 }[/math] до [math]\displaystyle{ \frac{n^2 - n}{2} }[/math].

Таким образом, в параллельном варианте, в отличие от последовательного, вычисления квадратных корней и делений будут определять довольно значительную долю требуемого времени. При реализации на конкретных архитектурах наличие в отдельных ярусах ЯПФ отдельных вычислений квадратных корней может породить и другие проблемы. Например, при реализации на ПЛИСах остальные вычисления (деления и тем более умножения и сложения/вычитания) могут быть конвейеризованы, что даёт экономию и по ресурсам на программируемых платах; вычисления же квадратных корней из-за их изолированности приведут к занятию ресурсов на платах, которые будут простаивать большую часть времени. Таким образом, общая экономия в 2 раза, из-за которой метод Холецкого предпочитают в случае симметричных задач тому же методу Гаусса, в параллельном случае уже имеет место вовсе не по всем ресурсам, и главное - не по требуемому времени.

При этом использование режима накопления требует совершения умножений и вычитаний в режиме двойной точности, а в параллельном варианте это означает, что практически все промежуточные вычисления для выполнения метода Холецкого в режиме накопления должны быть двойной точности. В отличие от последовательного варианта это означает увеличение требуемой памяти почти в 2 раза.

При классификации по высоте ЯПФ, таким образом, метод Холецкого относится к алгоритмам со сложностью [math]\displaystyle{ O(n) }[/math]. При классификации по ширине ЯПФ его сложность будет [math]\displaystyle{ O(n^2) }[/math].

Входные и выходные данные алгоритма

Входные данные: n - размерность матрицы [math]\displaystyle{ H_n }[/math].

Объём входных данных: 1 (число n).

Выходные данные: матрица Адамара [math]\displaystyle{ H_n }[/math] размерностью [math]\displaystyle{ 2^n\times 2^n }[/math].

Объём выходных данных: [math]\displaystyle{ 2^{2n} }[/math].

Свойства алгоритма

Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов, как хорошо видно, является квадратичным (отношение кубической к линейной).

При этом вычислительная мощность алгоритма, как отношение числа операций к суммарному объему входных и выходных данных – всего лишь линейна.

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

Дуги информационного графа, исходящие из вершин, соответствующих операциям квадратного корня и деления, образуют пучки т. н. рассылок линейной мощности (то есть степень исхода этих вершин и мощность работы с этими данными — линейная функция от порядка матрицы и координат этих вершин). При этом естественно наличие в этих пучках «длинных» дуг. Остальные дуги локальны.

Наиболее известной является компактная укладка графа — его проекция на треугольник матрицы, который перевычисляется укладываемыми операциями. При этом «длинные» дуги можно убрать, заменив более дальнюю пересылку комбинацией нескольких ближних (к соседям).

Эквивалентное возмущение [math]\displaystyle{ M }[/math] у метода Холецкого всего вдвое больше, чем возмущение [math]\displaystyle{ \delta A }[/math], вносимое в матрицу при вводе чисел в компьютер: [math]\displaystyle{ ||M||_{E} \leq 2||\delta A||_{E} }[/math]

Это явление обусловлено положительной определённостью матрицы. Среди всех используемых разложений матриц это наименьшее из эквивалентных возмущений.

ЧАСТЬ Программная реализация алгоритма

Особенности реализации последовательного алгоритма

В простейшем варианте генерацию матрицы Адамара для нашего алгоритма можно записать так: 

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) // 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;
    }
}

Локальность данных и вычислений

Возможные способы и особенности параллельной реализации алгоритма

Масштабируемость алгоритма и его реализации

Динамические характеристики и эффективность реализации алгоритма

Выводы для классов архитектур

Существующие реализации алгоритма

Примеры библиотек, содержащих функции генерации матриц Адамара:

  1. MATLAB
  2. WOLFRAM
  3. Python

Литература

  1. ↑ Воеводин В.В. Математические основы параллельных вычислений// М.: Изд. Моск. ун-та, 1991. 345 с.
  2. ↑ Воеводин В.В., Воеводин Вл.В. Параллельные вычисления. – СПб.: БХВ - Петербург, 2002. – 608 с.
  3. ↑ Фролов А.В.. Принципы построения и описание языка Сигма. Препринт ОВМ АН N 236. М.: ОВМ АН СССР, 1989.
Источник — https://xn--80aecvaldz.xn--80ahdrnmvk.xn--p1ai/w/ru/index.php?title=Участник:Бугаков_Юрий/Построение_матрицы_Адамара_произвольного_размера&oldid=7590