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

Авторы страницы: Сычева Е.А и Хахалин А.С.

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

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

Разреженная матрица — матрица с большим количеством нулевых элементов.

Среди специалистов нет единства в определении того, какое именно количество ненулевых элементов делает матрицу разреженной. Матрицу порядка [math]n[/math] называют разреженной, если число ее ненулевых элементов :

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

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

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

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

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

1.1.1 Диагональная схема хранения ленточных матриц

С ленточными матрицами связана простейшая и наиболее широко применяемая стратегия использования нулевых элементов матрицы. Матрицу [math]A[/math] называют ленточной, если все ее ненулевые элементы заключены внутри ленты, образованной диагоналями, параллельными главной диагонали. Таким образом, [math]a_{ij} = 0[/math], если [math]|i - j| \gt b[/math], и [math]a_{k, k-b} \ne 0[/math], либо [math]a_{k, k+b} \ne 0[/math] хотя бы для одного значения [math]k[/math]. Здесь [math]b[/math] — полуширина, а [math]2b + 1[/math] — ширина ленты. Лентой матрицы [math]A[/math] называется множество элементов, для которых [math]|i - j| \leq b[/math]. Другими словами, для всякой строки [math]i[/math] ленте принадлежат все элементы со столбцовыми индексами от [math]i - b[/math] до [math]i + b[/math], т. е. [math]2b + 1[/math] элементов.

Если матрица симметрична, то достаточно хранить ее полу­ленту. Верхняя полулента состоит из элементов, находящихся в верхней части ленты, т. е. таких, что [math]0 \lt j - i \leq \beta [/math]; нижняя полулента состоит из элементов нижней части ленты, т. е. таких, что [math]0 \lt i - j \leq \beta[/math]; в обоих случаях в строке [math]\beta[/math] элементов.

Диагональная схема хранения симметричной ленточной ма­трицы [math]A[/math] порядка [math]n[/math] и полуширины ленты [math]\beta[/math] представляет собой массив [math]AN (I, J)[/math] размером [math]n \times (\beta + 1)[/math]. Главная диагональ хранится в последнем столбце, а нижние кодиагонали — в остальных столбцах со сдвигом на одну позицию вниз при каждом смещении влево.

Для несимметричной матрицы [math]A[/math] необходим массив размером [math]n \times(2\beta + 1)[/math]; нижняя полулента и главная диагональ хранятся прежним образом, а верхние кодиагонали — в правой части массива со сдвигом на одну позицию вверх при каждом смещении вправо. Диагональная схема удобна, если [math]\beta \lt \lt n [/math]. Она является схемой прямого доступа в том смысле, что имеется простое взаимно однозначное соответствие между положением элемента в матрице [math]A[/math] и его положением в массиве: [math]a_{ij}[/math] хранится компонентой [math]AN (i, j - i + \beta + 1)[/math].

1.1.2 Профильная схема хранения симметричных матриц

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

Для каждой строки [math]i[/math] симметричной матрицы [math]A[/math] положим:

[math]\beta_{i} = i - j_{min} (i)[/math]

где [math]j_{min} (i)[/math] — минимальный столбцовый индекс строки [math]i[/math], для которого [math]a_{ij} \ne 0[/math]. Таким образом, первый ненулевой элемент строки [math]i[/math] находится на [math]\beta_{i}[/math]- позиций левее главной диагонали.

Оболочка матрицы [math]A[/math] — это множество элементов [math]a_{ij}[/math], для которых [math]0 \lt i - j \lt \beta_{i}[/math] В строке [math]i[/math] оболочке принадлежат все элементы со столбцовыми индексами от [math]j_{min} (i)[/math] до [math]i - 1[/math], всего [math]\beta_{i}[/math] элементов. Диагональные элементы не входят в оболочку. Про­филь матрицы [math]A[/math] определяется как число элементов в оболочке:

[math]profile (A) = \sum\limits_{i}\beta_{i}[/math]

При использовании профильной схемы все элементы оболочки, упорядоченные по строкам, хранятся, включая нули, в одномер­ном массиве [math]AN[/math]. Диагональный элемент данной строки помещается в ее конец. Длина массива [math]AN[/math] равна сумме профиля и порядка [math]A[/math]. Кроме того, необходим массив указателей [math]IA[/math]; элементы этого массива - указатели расположения диагональ­ных элементов в [math]AN[/math]. Так, при [math]i \gt 1[/math] элементы строки [math]i[/math] находятся в позициях от [math]IA (i - 1) + 1[/math] до [math]IA (i)[/math]. Единственный элемент [math]a_{11}[/math] 1-й строки хранится в [math]AN (1)[/math].

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

1.1.3 Связные схемы разреженного хранения

Рассмотрим разреженную матрицу [math]A[/math]. Ее ненулевые элементы, которых очень мало в сравнении с числом нулей, рассеяны по всей матрице, образуя портрет матрицы или шаблон нулей—ненулей. Опишем схему хранения предложенную Кнутом. Ненулевые элементы хранятся в компактной форме и в произвольном порядке в одномерном массиве [math]AN[/math]. Информация о положении элементов может храниться двумя дополнительными параллельными одномерными массивами [math]I[/math] и [math]J[/math]; здесь для каждого ненулевого элемента содержатся его строчный и столбцовый индексы. Таким образом, для каждого [math]a_{ij}\ne0[/math] в памяти находится тройка [math](a_{ij},i,j)[/math]. Далее, чтобы можно было легко отыскать элементы произвольной строки или столбца, необходимы ещё пара указателей для каждой тройки, а также указатели входа для строк и столбцов, сообщающие начало каждого строчного или столбцового списка. Пусть [math]NR[/math] ("next nonzero element in the same row — "следующий ненулевой элемент той же строки") — массив, хранящий строчные указатели, а [math]NC[/math] ("next nonzero element in the same column — "следующий ненулевой элемент того же столбца") — массив столбцовых указателей. Пять массивов [math]AN,I,J,NR,NC[/math] имеют общую длину, и их одноименные позиции соответствуют друг другу. Пусть [math]JR[/math] и [math]JC[/math] — массивы, содержащие указатели входа для строк и столбцов, расположенные в соответсвии с порядком строк и столбцов матрицы. Тогда, например, для просмотра [math]i[/math] строки необходимо

1. определить входной индекс [math]j=JR[i][/math]
2. определить элемент [math]AN[j][/math]
3. найти индекс следующего элемента [math]j=NR[j][/math]
4. если [math]j=0[/math], то мы обошли все элемент [math]i[/math]-ой строки, иначе, перейти в п.2

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

1.1.4 Разреженный строчный формат

Значения ненулевых элементов и соответсвующие столбцовые индексы хранятся в этой схеме по строкам в двух массивах — [math]AN[/math] и [math]JA[/math] соответсвенно. Дополнительно также используется массив указателей [math]IA[/math], отмечающих позиции массивов [math]AN[/math] и [math]JA[/math], с которых начинается описание очередной строки. Дополнительная компонента в [math]IA[/math] содержит указатель первой свободной позиции в [math]AN[/math] и [math]JA[/math]. В общем случае описание [math]r[/math]-й строки [math]A[/math] хранится в позициях с [math]IA (r)[/math] до [math]IA (r + 1) - 1[/math] массивов [math]JA[/math] и [math]AN[/math], за исключением равенства [math]IA (r + 1) = IA (r)[/math], означающего, что [math]r[/math]-я строка пуста. Если матрица [math]A[/math] имеет [math]m[/math] строк, то [math]IA[/math] содержит [math]m + 1[/math] позиций.

Схема предъявляет минимальные требования к памяти и в то же время оказывается очень удобной для нескольких важных операций над разряженными матрицами: сложения, умножения, перестановок строк и столбцов, транспонирования, решения линейных систем с разреженными матрицами коэффициентов как прямым, так и итерационными методами.

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

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

Для вычислений с разреженными матрицами создано несколько библиотек:

• SparseLib++ (C++)
• uBLAS (C++, в составе Boost)
• SPARSPAK (Fortran)
• CSparse (C)

и другие.

3 Литература

<references \>