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

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

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

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

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

[math]\displaystyle{ H_n\,H_n^T = n\,E_n }[/math]

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

Матрицу Адамара можно определить эквивалентным образом:

[math]\displaystyle{ H_{n} = H_1\otimes H_{n-1} }[/math]

[math]\displaystyle{ \begin{align} H_1 = &\begin{pmatrix}\begin{array}{rr} 1 & 1\\ 1 & -1 \end{array}\end{pmatrix} \end{align} }[/math]

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

[math]\displaystyle{ H_{n} = \begin{pmatrix} H_{n-1} & H_{n-1}\\ H_{n-1} & -H_{n-1}\end{pmatrix} }[/math]

[math]\displaystyle{ \begin{align} H_1 = &\begin{pmatrix}\begin{array}{rr} 1 & 1\\ 1 & -1 \end{array}\end{pmatrix} \end{align} }[/math]

С матрицами Адамара связан ряд нерешенных проблем, одна из которых состоит в следующем. Мы уже видели, что порядок n матрицы Адамара при n ≥ 3 может быть лишь четным. Более того, при n ≥ 4 порядок обязан делиться на 4. До сих пор остается открытым вопрос: для любого ли n, кратного 4, существует матрица Адамара порядка n? Неизвестно, в частности, существует ли матрица Адамара порядка 268 (это наименьший порядок, кратный 4, для которого матрица Адамара еще не построена).

Часто размерность матрицы Адамара считают равной степени 2 и нормализируют её. Определение принимает следующий вид:

Матрицей Адамара H_n эта матрица размера 2ⁿ × 2ⁿ, для которой справедливо равенство

[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_n }[/math]. Для этого представим порядковые номера k и l элемента [math]\displaystyle{ h_{kl} }[/math] в виде двоичного разложения

[math]\displaystyle{ k = \sum^{m-1}_{i=0} {k_i 2^i} = k_{m-1} 2^{m-1} + k_{m-2} 2^{m-2} + \cdots + k_1 2 + k_0 }[/math]

и

[math]\displaystyle{ l = \sum^{m-1}_{i=0} {l_i 2^i} = l_{m-1} 2^{m-1} + l_{m-2} 2^{m-2} + \cdots + l_1 2 + l_0 }[/math]

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

Тогда каждый элемент матрицы Адамара [math]\displaystyle{ H_n }[/math] имеет вид

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

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

[math]\displaystyle{ \begin{align} H_0 = &+1\\ H_1 = \frac{1}{\sqrt2} &\begin{pmatrix}\begin{array}{rr} 1 & 1\\ 1 & -1 \end{array}\end{pmatrix} \end{align} }[/math]

[math]\displaystyle{ \begin{align} 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}\\ 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}\\ \end{align} }[/math]

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

Исходные данные: положительно определённая симметрическая матрица [math]\displaystyle{ A }[/math] (элементы [math]\displaystyle{ a_{ij} }[/math]).

Вычисляемые данные: нижняя треугольная матрица [math]\displaystyle{ L }[/math] (элементы [math]\displaystyle{ l_{ij} }[/math]).

Формулы метода:

[math]\displaystyle{ \begin{align} l_{11} & = \sqrt{a_{11}}, \\ l_{j1} & = \frac{a_{j1}}{l_{11}}, \quad j \in [2, n], \\ l_{ii} & = \sqrt{a_{ii} - \sum_{p = 1}^{i - 1} l_{ip}^2}, \quad i \in [2, n], \\ l_{ji} & = \left (a_{ji} - \sum_{p = 1}^{i - 1} l_{ip} l_{jp} \right ) / l_{ii}, \quad i \in [2, n - 1], j \in [i + 1, n]. \end{align} }[/math]

Существует также блочная версия метода, однако в данном описании разобран только точечный метод.

В ряде реализаций деление на диагональный элемент выполняется в два этапа: вычисление [math]\displaystyle{ \frac{1}{l_{ii}} }[/math] и затем умножение на него всех (видоизменённых) [math]\displaystyle{ a_{ji} }[/math] . Здесь мы этот вариант алгоритма не рассматриваем. Заметим только, что он имеет худшие параллельные характеристики, чем представленный.

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

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

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

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

Как записано выше основную часть метода составляют множественные вычисления элементов

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

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

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

1. [math]\displaystyle{ l_{11}= \sqrt{a_{11}} }[/math]

2. [math]\displaystyle{ l_{j1}= \frac{a_{j1}}{l_{11}} }[/math] (при [math]\displaystyle{ j }[/math] от [math]\displaystyle{ 2 }[/math] до [math]\displaystyle{ n }[/math]).

Далее для всех [math]\displaystyle{ i }[/math] от [math]\displaystyle{ 2 }[/math] до [math]\displaystyle{ n }[/math] по нарастанию выполняются

3. [math]\displaystyle{ l_{ii} = \sqrt{a_{ii} - \sum_{p = 1}^{i - 1} l_{ip}^2} }[/math] и

4. (кроме [math]\displaystyle{ i = n }[/math]): [math]\displaystyle{ l_{ji} = \left (a_{ji} - \sum_{p = 1}^{i - 1} l_{ip} l_{jp} \right ) / l_{ii} }[/math] (для всех [math]\displaystyle{ j }[/math] от [math]\displaystyle{ i + 1 }[/math] до [math]\displaystyle{ n }[/math]).

После этого (если [math]\displaystyle{ i \lt n }[/math]) происходит переход к шагу 3 с бо́льшим [math]\displaystyle{ i }[/math].

Особо отметим, что вычисления сумм вида [math]\displaystyle{ a_{ji} - \sum_{p = 1}^{i - 1} l_{ip} l_{jp} }[/math] в обеих формулах производят в режиме накопления вычитанием из [math]\displaystyle{ a_{ji} }[/math] произведений [math]\displaystyle{ l_{ip} l_{jp} }[/math] для [math]\displaystyle{ p }[/math] от [math]\displaystyle{ 1 }[/math] до [math]\displaystyle{ i - 1 }[/math], c нарастанием [math]\displaystyle{ p }[/math].

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

Для разложения матрицы порядка n методом Холецкого в последовательном (наиболее быстром) варианте требуется:

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

Умножения и сложения (вычитания) составляют основную часть алгоритма.

При этом использование режима накопления требует совершения умножений и вычитаний в режиме двойной точности (или использования функции вроде DPROD в Фортране), что ещё больше увеличивает долю умножений и сложений/вычитаний во времени, требуемом для выполнения метода Холецкого.

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

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

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

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

Первая группа вершин расположена в одномерной области, соответствующая ей операция вычисляет функцию 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].

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

Входные данные: плотная матрица [math]\displaystyle{ A }[/math] (элементы [math]\displaystyle{ a_{ij} }[/math]). Дополнительные ограничения:

• [math]\displaystyle{ A }[/math] – симметрическая матрица, т. е. [math]\displaystyle{ a_{ij}= a_{ji}, i, j = 1, \ldots, n }[/math].
• [math]\displaystyle{ A }[/math] – положительно определённая матрица, т. е. для любых ненулевых векторов [math]\displaystyle{ \vec{x} }[/math] выполняется [math]\displaystyle{ \vec{x}^T A \vec{x} \gt 0 }[/math].

Объём входных данных: [math]\displaystyle{ \frac{n (n + 1)}{2} }[/math] (в силу симметричности достаточно хранить только диагональ и над/поддиагональные элементы). В разных реализациях эта экономия хранения может быть выполнена разным образом. Например, в библиотеке, реализованной в НИВЦ МГУ, матрица A хранилась в одномерном массиве длины [math]\displaystyle{ \frac{n (n + 1)}{2} }[/math] по строкам своего нижнего треугольника.

Выходные данные: нижняя треугольная матрица [math]\displaystyle{ L }[/math] (элементы [math]\displaystyle{ l_{ij} }[/math]).

Объём выходных данных: [math]\displaystyle{ \frac{n (n + 1)}{2} }[/math] (в силу треугольности достаточно хранить только ненулевые элементы). В разных реализациях эта экономия хранения может быть выполнена разным образом. Например, в той же библиотеке, созданной в НИВЦ МГУ, матрица [math]\displaystyle{ L }[/math] хранилась в одномерном массиве длины [math]\displaystyle{ \frac{n (n + 1)}{2} }[/math] по строкам своей нижней части.

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

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

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

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

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

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

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

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