Участник:Blizn/Хранение ненулевых элементов разреженной матрицы. Умножение разреженной матрицы на вектор.

Файл:Symbol wait.svg

Эта работа ждет рассмотрения преподавателем Дата последней правки страницы: 02.11.2016 Авторы этой статьи считают, что задание выполнено.

Умножение разреженной матрицы на вектор
Последовательный алгоритм
Последовательная сложность	[math]\displaystyle{ O(l) }[/math]
Объём входных данных	[math]\displaystyle{ 2l+m+n+1 }[/math]
Объём выходных данных	[math]\displaystyle{ n }[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы	[math]\displaystyle{ O(m) }[/math]
Ширина ярусно-параллельной формы	[math]\displaystyle{ O(n) }[/math]

Выполнила: И.В. Близнякова (611 группа).

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

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

Разрежённая матрица — это матрица с преимущественно нулевыми элементами. В противном случае, если бо́льшая часть элементов матрицы ненулевые, матрица считается плотной.

Среди специалистов нет единства в определении того, какое именно количество ненулевых элементов делает матрицу разрежённой. Разные авторы предлагают различные варианты. Для матрицы порядка n число ненулевых элементов:

• есть [math]\displaystyle{ O(n) }[/math]. Такое определение подходит разве что для теоретического анализа асимптотических свойств матричных алгоритмов;
• в каждой строке не превышает 10 в типичном случае;
• ограничено [math]\displaystyle{ n^{1+\gamma} }[/math], где [math]\displaystyle{ \gamma \lt 1 }[/math].
• таково, что для данного алгоритма и вычислительной системы имеет смысл извлекать выгоду из наличия в ней нулей.

Хранение разреженной матрицы

Формат RR(C)O

Рассмотрим сначала формат RR(C)O. Сокращенное название данного формата происходит от английского словосочетания "Row - wise Representation Complete and Ordered" (строчное представление, полное и упорядоченное). В данном формате вместо одного двумерного массива, используются три одномерных. Значения ненулевых элементов матрицы и соответствующие им столбцовые индексы хранятся в этом формате по строкам в двух массивах [math]\displaystyle{ AN }[/math] и [math]\displaystyle{ JA }[/math]. Массив указателей [math]\displaystyle{ IA }[/math], используется для ссылки на компоненты массивов [math]\displaystyle{ AN }[/math] и [math]\displaystyle{ JA }[/math], с которых начинается описание очередной строки. Последняя компонента массива [math]\displaystyle{ IA }[/math] содержит указатель первой свободной компоненты в массивах [math]\displaystyle{ AN }[/math] и [math]\displaystyle{ JA }[/math], т.е. равна числу ненулевых элементов матрицы, увеличенному на единицу. Здесь уместно привести пример.

Рассмотрим матрицу [math]\displaystyle{ A }[/math]:

[math]\displaystyle{ \begin{pmatrix} 0 & 0 & 0 & 2 & 0 \\ 1 & 0 & 3 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 6 & 0 & 0 & 0 \\ 0 & 0 & 4 & 0 & 0 \\ \end{pmatrix} }[/math],

тогда ее представление в формате RR(C)O будет иметь вид:

  IA = [ 1 2 4 4 5 6 ]
  JA = [ 4 1 3 2 3 ]
  AN = [ 2 1 3 6 4 ]

Т.е. массив [math]\displaystyle{ AN }[/math] содержит все не нулевые элементы исходной матрицы [math]\displaystyle{ A }[/math], массив [math]\displaystyle{ JA }[/math] номер столбца в котором находится соответствующий элемент из [math]\displaystyle{ AN }[/math] и наконец массив [math]\displaystyle{ IA }[/math] содержит номер с которого начинается описание элементов в массивах [math]\displaystyle{ JA }[/math] и [math]\displaystyle{ AN }[/math]. Таким образом информация об элементах 2-ой строки матрицы хранится в элементах с [math]\displaystyle{ IA[2] = 2 }[/math] по [math]\displaystyle{ IA[3] - 1 = 3 }[/math] включительно массивов [math]\displaystyle{ JA }[/math] и [math]\displaystyle{ AN }[/math]. Можно обратить внимание, что [math]\displaystyle{ IA[3] = IA[4] = 4 }[/math], а это означает, что 3-я строка матрицы [math]\displaystyle{ A }[/math] нулевая.

В общем случае описание [math]\displaystyle{ r }[/math]-й строки матрицы A хранится в компонентах с [math]\displaystyle{ IA[r] }[/math] до [math]\displaystyle{ IA[r + 1] - 1 }[/math] включительно массивов [math]\displaystyle{ AN }[/math] и [math]\displaystyle{ JA }[/math]. Если [math]\displaystyle{ IA[r + 1] = IA[r] }[/math], то это означает, что [math]\displaystyle{ r }[/math]-я строка нулевая. Количество элементов в массиве [math]\displaystyle{ IA }[/math] на единицу больше, чем число строк исходной матрицы, а количество элементов в массивах [math]\displaystyle{ JA }[/math] и [math]\displaystyle{ AN }[/math] равно числу ненулевых элементов исходной матрицы.

Данный способ представления называют полным, поскольку представлена вся матрица [math]\displaystyle{ A }[/math], упорядоченным, поскольку элементы каждой строки матрицы [math]\displaystyle{ A }[/math] хранятся в соответствии с возрастанием столбцовых индексов, и строчным, поскольку информация о матрице [math]\displaystyle{ A }[/math] указывается по строкам.

Массивы [math]\displaystyle{ IA }[/math] и [math]\displaystyle{ JA }[/math] представляют портрет (структуру) матрицы [math]\displaystyle{ A }[/math], задаваемый как множество списков смежности ассоциированного с [math]\displaystyle{ A }[/math] графа. Если алгоритм, реализующий какую-либо операцию над разреженными матрицами, разбит на этапы символической обработки, на котором определяется портрет результирующей матрицы, и численной обработки, на котором определяются значения элементов результирующей матрицы, то массивы [math]\displaystyle{ IA }[/math] и [math]\displaystyle{ JA }[/math] заполняются на первом этапе, а массив [math]\displaystyle{ AN }[/math] - на втором.

Формат RR(C)U

Рассмотрим теперь формат RR(C)U.

Сокращенное название данного формата происходит от английского словосочетания "Row - wise Representation Complete and Unordered" (строчное представление, полное, но неупорядоченное). Формат RR(C)U отличается от RR(C)O тем, что в данном случае соблюдается упорядоченность строк, но внутри каждой строки элементы исходных матриц могут храниться в произвольном порядке. Для матрицы [math]\displaystyle{ A }[/math] нашего примера вполне можно было бы использовать и строчное представление, полное, но неупорядоченное такое:

  IA = [ 1 2 4 4 5 6 ]
  JA = [ 4 3 1 2 3 ]
  AN = [ 2 1 3 6 4 ]

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

Замечания

Несколько замечаний по поводу рассмотренных форматов представления:

1. Очевидно, что представление матрицы в формате RR(C)O так же является и представлением в формате RR(C)U, но не наоборот.
2. Из представления матрицы в формате RR(C) нельзя получить информацию о точном количестве столбцов исходной матрицы.
3. Целесообразно (в вопросе экономии памяти) использовать представления RR(C) в случае, если матрица содержит значительное число нулевых элементов.

Умножение разреженной матрицы на вектор

Важным приложением этих алгоритмов является вычисление векторов Ланцоша, необходимое при итерационном решении линейных уравнений методом сопряженных градиентов, а также при вычислении собственных значений и собственных векторов матрицы. Достоинство этих процедур, с вычислительной точки зрения, состоит в том, что единственная требуемая матричная операция - это повторное умножение матрицы на последовательность заполненных векторов; сама матрица не меняется.

Мы рассмотрим умножение разреженной матрицы общего вида, хранимой в форме RR(C)U посредством массивов [math]\displaystyle{ IA }[/math], [math]\displaystyle{ JA }[/math], [math]\displaystyle{ AN }[/math] на заполненный вектор-столбец.

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

Исходные данные: разреженная матрица общего вида [math]\displaystyle{ A }[/math] с элементами [math]\displaystyle{ a_{ij} }[/math] ([math]\displaystyle{ i = 1,...,n }[/math] и [math]\displaystyle{ j = 1,...,m }[/math]). Заполненный вектор-столбец [math]\displaystyle{ b }[/math] с элементами [math]\displaystyle{ b_{j} }[/math] ([math]\displaystyle{ j =1,...,m }[/math]).

Вычисляемые данные: заполненный вектор-столбец [math]\displaystyle{ c }[/math] с элементами [math]\displaystyle{ c_{i} }[/math] ([math]\displaystyle{ i = 1,...,n }[/math]).

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

[math]\displaystyle{ c_{i} = \sum\limits_{k = 1}^{l_{i}} a_{i,j=j(k)}b_{j=j(k)} }[/math],

где [math]\displaystyle{ l_{i} }[/math] - количество ненулевых элементов строки [math]\displaystyle{ i }[/math] матрицы [math]\displaystyle{ A }[/math], [math]\displaystyle{ j(k) }[/math] - индекс [math]\displaystyle{ k }[/math]-го ненулевого элемента матрицы [math]\displaystyle{ A }[/math].

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

Вычислительное ядро последовательной версии умножения разреженной матрицы на вектор можно составить из множественных (всего их [math]\displaystyle{ n }[/math]) вычислений скалярных произведений строк матрицы:

[math]\displaystyle{ c_{i} = \sum\limits_{k = 1}^{l_{i}} a_{i,j=j(k)}b_{j=j(k)} }[/math].

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

Как записано и в описании ядра алгоритма, основную часть метода составляют множественные (всего [math]\displaystyle{ n }[/math]) вычисления сумм:

[math]\displaystyle{ c_{i} = \sum\limits_{k = 1}^{l_{i}} a_{i,j=j(k)}b_{j=j(k)} }[/math],

которые могут вычисляться в произвольном порядке.

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

Далее предполагаем, что разреженная матрица общего вида [math]\displaystyle{ A }[/math] хранится в форме RR(C)U посредством массивов [math]\displaystyle{ IA }[/math], [math]\displaystyle{ JA }[/math], [math]\displaystyle{ AN }[/math]. Последовательность исполнения метода следующая:

Выполнять для [math]\displaystyle{ i }[/math] от [math]\displaystyle{ 1 }[/math] до [math]\displaystyle{ n }[/math]

1. [math]\displaystyle{ c_{i} = 0 }[/math]
2. [math]\displaystyle{ IAA = IA[i] }[/math]
3. [math]\displaystyle{ IAB = IA[i + 1] - 1 }[/math]
4. [math]\displaystyle{ c_{i} = \sum\limits_{k = IAA}^{IAB} AN[k]b_{JA[k]} }[/math].

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

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

Для умножения разреженной матрицы общего вида, хранимой в форме RR(C)U, размером [math]\displaystyle{ n \times m }[/math] на заполненный вектор [math]\displaystyle{ m \times 1 }[/math] в последовательном (наиболее быстром) варианте требуется:

• [math]\displaystyle{ l }[/math] сложений,
• [math]\displaystyle{ l }[/math] умножений.

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

При классификации по последовательной сложности, таким образом, алгоритм умножения разреженной матрицы на вектор относится к алгоритмам [math]\displaystyle{ O(l) }[/math].

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

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

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

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

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

где [math]\displaystyle{ k = k(i) = IA[i+1]-IA[i] }[/math].

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

• [math]\displaystyle{ a }[/math] - элемент входных данных, а именно [math]\displaystyle{ AN[IA_{i}+j-1] }[/math].
• [math]\displaystyle{ b }[/math] - элемент входных данных, а именно [math]\displaystyle{ b[IA_{i}+j-1] }[/math].

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

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

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

где [math]\displaystyle{ k = k(i) = IA[i+1]-IA[i] }[/math].

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

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

Результат срабатывания операции:

• при [math]\displaystyle{ j \lt k - 1 }[/math] является промежуточным данным алгоритма;
• при [math]\displaystyle{ j = k - 1 }[/math] является выходным данным [math]\displaystyle{ c_{i} }[/math].

Описанный граф можно посмотреть на рис.1. Здесь вершины первой группы обозначены красным цветом и отмечены знаком умножения, вершины второй - зелёным цветом и знаком сложения. Вершины, соответствующие входным данным обозначены белым цветом и выходным - синим.

Рисунок 1. Граф алгоритма умножения разреженной матрицы на вектор

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

Для умножения разреженной матрицы общего вида, хранимой в форме RR(C)U, размерности [math]\displaystyle{ n \times m }[/math] на заполненный вектор [math]\displaystyle{ m \times 1 }[/math] в параллельном варианте требуется последовательно выполнить следующие ярусы:

• Не более чем [math]\displaystyle{ m }[/math] сложений и умножений ([math]\displaystyle{ n }[/math] вычислений в каждом из ярусов)

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

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

Входные данные: разреженная матрица общего вида [math]\displaystyle{ A }[/math] размерности [math]\displaystyle{ n \times m }[/math], хранимая в форме RR(C)U посредством массивов [math]\displaystyle{ IA }[/math], [math]\displaystyle{ JA }[/math], [math]\displaystyle{ AN }[/math]. Заполненный вектор-столбец [math]\displaystyle{ b }[/math] с элементами [math]\displaystyle{ b_{j} }[/math] размерности [math]\displaystyle{ m \times 1 }[/math].

Объём входных данных: [math]\displaystyle{ 2l+m+n+1 }[/math], где [math]\displaystyle{ l }[/math] - количество ненулевых элементов в матрице [math]\displaystyle{ A }[/math].

Выходные данные: заполненный вектор-столбец [math]\displaystyle{ c }[/math] с элементами [math]\displaystyle{ c_{i} }[/math] ([math]\displaystyle{ i = 1,...,n }[/math]).

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

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

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

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

Пусть [math]\displaystyle{ l }[/math] - количество ненулевых элементов матрицы [math]\displaystyle{ A }[/math] размерности [math]\displaystyle{ n \times m }[/math]. Пусть [math]\displaystyle{ x }[/math] - объем памяти, используемый для хранения значения элемента матрицы, [math]\displaystyle{ y }[/math] - объём памяти, используемый для хранения номера столбца или строки. В таком случае для хранения матрицы в стандартном представлении нам потребуется объем памяти, равный [math]\displaystyle{ x \cdot n \cdot m }[/math], для хранения в формате RR(C) - [math]\displaystyle{ y(n + 1) + (x + y)l }[/math]. Таким образом хранение в формате RR(C) не является эффективным для матриц, в которых [math]\displaystyle{ l \gt \frac{xnm - y(n+1)}{x+y} }[/math].

Программная реализация алгоритма

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

Существующие реализации:

SciPy [1],

MatLab [2],

Intel MKL [3].