Что такое разреженная матрица

Разреженные матрицы

Разреженные матрицы обеспечивают эффективное устройство хранения данных double или данных logical , которые имеют большой процент нулей. В то время как полный (или плотный ) матрицы хранят каждый элемент в памяти независимо от значения, разреженные матрицы хранят только ненулевые элементы и их индексы строки. Поэтому использование разреженных матриц может значительно уменьшать объем памяти, требуемый для хранения данных.

Весь MATLAB ® встроенная арифметика, логическая, и операции индексации, может быть применен к разреженным матрицам, или к смесям разреженных и полных матриц. Операции на разреженных матрицах возвращают разреженные матрицы, и операции на полных матрицах возвращают полные матрицы. Для получения дополнительной информации смотрите Вычислительные Преимущества Разреженных матриц Построения и Разреженных матриц.

Функции

Создание

spalloc	Выделите место для разреженной матрицы
spdiags	Извлеките и создайте разреженную полосу и диагональные матрицы
speye	Разреженная единичная матрица
sprand	Разреженная равномерно распределенная случайная матрица
sprandn	Разреженная нормально распределенная случайная матрица
sprandsym	Разреженная симметричная случайная матрица
sparse	Создайте разреженную матрицу
spconvert	Импортируйте из внешнего формата разреженной матрицы

Манипуляция

issparse	Определите, разреженно ли введенный
nnz	Количество ненулевых элементов матрицы
nonzeros	Ненулевые элементы матрицы
nzmax	Сумма устройства хранения данных выделяется для ненулевых элементов матрицы
spfun	Примените функцию к ненулевым элементам разреженной матрицы
spones	Замените ненулевые элементы разреженной матрицы на единицы
spparms	Установите параметры для стандартных программ разреженной матрицы
spy	Визуализируйте шаблон разреженности
find	Найдите индексы и значения ненулевых элементов
full	Преобразуйте разреженную матрицу в полное устройство хранения данных

Переупорядочение алгоритмов

dissect	Вложенная перестановка рассечения
amd	Аппроксимируйте минимальную перестановку степени
colamd	Столбец аппроксимированная минимальная перестановка степени
colperm	Разреженная перестановка столбца на основе ненулевого количества
dmperm	Разложение Дулмаге-Мендельсона
randperm	Случайная перестановка
symamd	Симметричная аппроксимированная минимальная перестановка степени
symrcm	Разреженное упорядоченное расположение обратного алгоритма Катхилла-Макки

Итеративные методы и предварительные формирователи

pcg	Предобусловленный метод сопряженных градиентов
minres	Метод минимальных невязок
symmlq	Симметричный метод LQ
gmres	Обобщенный метод минимальных невязок (с перезапусками)
bicg	Бисопряженный метод градиентов
bicgstab	Бисопряженные градиенты стабилизировали метод
bicgstabl	Бисопряженные градиенты стабилизированный (l) метод
cgs	Методы сопряженных градиентов придали методу квадратную форму
qmr	Метод квази-минимальных невязок
tfqmr	Метод квази-минимальных невязок без транспонирования
lsqr	Метод LSQR
equilibrate	Матрица, масштабирующаяся для улучшенного создания условий
ichol	Неполная факторизация Холесского
ilu	Неполная LU-факторизация

Умножение разреженных матриц

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

Понятие разреженной матрицы можно определить многими способами, суть которых состоит в том, что в разреженной матрице «много» нулевых элементов. Обычно говорят, что матрица разрежена, если она содержит $O\left(n\right)$ отличных от нуля элементов. В противном случае матрица считается плотной. Типичным случаем разреженности является ограниченность числа ненулевых элементов в одной строке от 1 до k, где $k\ll n$ . Задачи линейной алгебры с разреженными матрицами возникают во многих областях, например, при решении дифференциальных уравнений в частных производных, при решении многомерных задач локальной оптимизации. В п. 7.3 мы уже познакомились с методами решения СЛАУ с трехдиагональной матрицей, которая, являясь ленточной, относится также и к классу разреженных матриц.

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

Многие алгоритмы, тривиальные для случая плотных матриц, в разреженном случае требуют более тщательного подхода. Во многих алгоритмах обработки разреженных матриц можно выделить два этапа: символический и численный. На символическом этапе формируется портрет результирующей матрицы (т.е. определяются места ненулевых элементов в структуре матрицы); на численном этапе определяются значения ненулевых элементов результирующей матрицы. В качестве примера типовых операций с разреженными матрицами в данной главе будет рассмотрен алгоритм умножения разреженной матрицы на плотный вектор, а также круг вопросов, возникающих при использовании метода Холецкого для решения систем линейных уравнений с разреженной матрицей.

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

Существуют различные форматы хранения разреженных матриц. Одни предназначены для хранения матриц специального вида (например, ленточных), другие обеспечивают работу с матрицами общего вида. Ниже рассмотрим некоторые весьма распространенные способы представления разреженных матриц, информацию о других способах можно найти, например, в [10].

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

массив ненулевых элементов матрицы A (обозначим его как values);
массив номеров строк матрицы A, соответствующих элементам массива values (обозначим его как rows);
массив номеров столбцов матрицы A, соответствующих элементам массива values (обозначим его как cols);

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

В качестве примера рассмотрим разреженную матрицу

$A=\begin1 & -1 & 0 & -3 & 0\\ -2 & 5 & 0 & 0 & 0\\ 0 & 0 & 4 & 6 & 4\\ -4 & 0 & 2 & 7 & 0\\ 0 & 8 & 0 & 0 & -5 \end$

( 7.40)

которая может быть представлена в координатном формате как

values=(1, -1, -3, -2, 5, 4, 6, 4, -4, 2, 7, 8, -5);

rows=(1, 1, 1, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5);

cols=(1, 2, 4, 1, 2, 3, 4, 5, 1, 3, 4, 2, 5).

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

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

В соответствии с рассматриваемой схемой для хранения матрицы A требуется три одномерных массива:

массив ненулевых элементов матрицы A, в котором они перечислены по строкам от первой до последней (обозначим его опять как values);
массив номеров столбцов для соответствующих элементов массива values (обозначим его как cols);
массив указателей позиций, с которых начинается описание очередной строки (обозначим его pointer). Описание k-й строки хранится в позициях с pointer[k]-й по (pointer[k+1]–1)-ю массивов values и cols. Если pointer[k]=pointer[k+1], то k-я строка пустая. Если матрица A состоит из n строк, то длина массива pointer будет n+1.

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

Для примера рассмотрим представление матрицы (7.40) в разреженном строчном формате:

values=(1, -1, -3, -2, 5, 4, 6, 4, -4, 2, 7, 8, -5);

cols=(1, 2, 4, 1, 2, 3, 4, 5, 1, 3, 4, 2, 5);

pointer=(1, 4, 6, 9, 12, 14).

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

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

После рассмотрения строчного формата хранения очевидным является и разреженный столбцовый формат. В этом случае ненулевые элементы матрицы A перечисляются в порядке их появления в столбцах матрицы, а не в строках. Все ненулевые элементы хранятся по столбцам в массиве values; индексы строк ненулевых элементов – в массиве rows; элементы массива pointer указывают на позиции, с которых начинается описание очередного столбца.

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

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

$A=\begin 1 & -1 & 0 & -3 & 0\\ -1 & 5 & 0 & 0 & 0\\ 0 & 0 & 4 & 6 & 4\\ -3 & 0 & 6 & 7 & 0\\ 0 & 0 & 4 & 0 & -5 \end» /> может быть представлена в разреженном строчном формате как values=(1, -1, -3, 5, 4, 6, 4, 7, -5); cols=(1, 2, 4, 2, 3, 4, 5, 4, 5); pointer=(1, 4, 5, 8, 9, 10). Если большинство диагональных элементов заданной симметричной матрицы отличны от нуля (например, у симметричной положительно определенной матрицы все диагональные элементы положительны), то они могут храниться в отдельном массиве BD, а разреженным форматом представляется только верхний треугольник В. Завершая обзор форматов хранения разреженных матриц, можно еще раз отметить, что выбор способа хранения матрицы полностью определяется алгоритмами, которые планируется в дальнейшем использовать для обработки данной матрицы. <h6>Базовые алгоритмы обработки разреженных матриц</h6> Реализация матричных операций является тривиальной в случае плотной или ленточной матрицы, но не столь очевидна для разреженных матриц, хранящихся в одном из экономичных форматов (здесь и далее мы будем рассматривать разреженный строчный формат хранения матрицы). <h6>Умножение матрицы на вектор</h6> Рассмотрим операцию умножения разреженной матрицы на плотный вектор. Результатом этой операции будет заполненный вектор, а не разреженная матрица. Поэтому в алгоритмах умножения вычисления выполняются прямо, без символического этапа. Среди алгоритмов, в которых встречается данная операция, можно отметить итерационные методы решения систем линейных уравнений. Достоинство данных методов, с вычислительной точки зрения, состоит в том, что единственная требуемая матричная операция, это повторное умножение матрицы на последовательность заполненных векторов; сама матрица не меняется. Итак, рассмотрим умножение разреженной матрицы общего вида, хранимой в строчном формате посредством массивов values, cols, pointer на заполненный векторстолбец b, хранимый в одномерном массиве. Результатом будет новый заполненный вектор <img decoding=$

также размещаемый в одномерном массиве.

Пусть n – число строк матрицы. Порядок, в котором накапливаются скалярные произведения, определяется порядком хранения элементов матрицы. Для каждой ее строки i мы находим с помощью массива индексов pointer значения первой pointer[i] и последней pointer[i+1]–1 позиций, занимаемых элементами строки i в массивах values и cols. Затем, чтобы вычислить скалярное произведение строки i и вектора b, мы просто просматриваем values и cols на отрезке от pointer[i] до pointer[i+1]–1; каждое значение, хранимое в cols[pointer[j]], есть столбцовый индекс и используется для извлечения из массива b элемента, который должен быть умножен на соответствующее число из массива values. Результат каждого умножения прибавляется к c[i].

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

Как правило, количество элементов в одной строке разреженной матрицы размера $n\times n$ ограничено некоторой константой k, $k\ll n$ . Значит, в процессе умножения такой матрицы на вектор количество вычислительных операций для получения скалярного произведения примерно одинаково для всех базовых подзадач. При использовании систем с общей памятью число потоков p будет меньше числа базовых подзадач n, и мы можем объединить базовые подзадачи таким образом, чтобы каждый поток выполнял несколько таких задач, соответствующих непрерывной последовательности строк матрицы А. В этом случае по окончании вычислений каждая базовая подзадача определяет набор элементов результирующего вектора с. Распределение подзадач между потоками может быть выполнено произвольным образом.

Транспонирование матрицы

A^T

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

A^T[i,j]=A^T[j,i]

и простейшая реализация данной операции в «плотном» случае оказывается эффективной. Однако в случае разреженной матрицы ситуация не столь очевидна, и простейшая реализация может существенно ухудшить показатели эффективности.

Будем формировать результирующую матрицу построчно. Для этого можно брать столбцы исходной матрицы и создавать из них строки результирующей матрицы. Но операция выделения из CRS-матрицы столбца №i является трудоемкой, т.к. данные в векторе values хранятся по строкам и для выборки данных по столбцу нужно просмотреть всю матрицу, что приводит к квадратичной (от числа ненулевых элементов) трудоемкости алгоритма. Необходимо другое решение. Подробно проблема транспонирования разреженной матрицы обсуждается в книге [10], здесь же рассмотрим основные идеи описанного в [10] алгоритма.

Сформируем N одномерных векторов для хранения целых чисел (IntVectors), а также N векторов для хранения вещественных чисел (RealVectors). N в данном случае соответствует числу столбцов исходной матрицы.
В цикле просмотрим все строки исходной матрицы, для каждой строки – все ее элементы. Пусть текущий элемент находится в строке i, столбце j, его значение равно v. Тогда добавим числа i и v в j-ые вектора для хранения целых и вещественных чисел (соответственно). Тем самым в векторах мы сформируем строки транспонированной матрицы.
Последовательно скопируем данные из векторов в CRS структуру транспонированной матрицы (сols и values), попутно формируя массив pointer.

Рассмотрим пример (7.17).

При обходе исходной матрицы A формируются вектора IntVectors и RealVectors (число 3 расположено в строке 0 и столбце 1, следовательно, в IntVectors[1] добавляется 0, в RealVectors[1] добавляется 3 и т.д.). Далее вектора последовательно формируют структуру : обратите внимание на порядок следования элементов в массиве values матрицы и его соответствие порядку элементов в массиве RealVectors (аналогично cols и IntVectors). Для формирования pointer достаточно подсчитать количество элементов в каждом из N векторов. pointer[0] всегда равно нулю, pointer[i]=pointer[i-1] + «Количество элементов в векторе i-1».

А^Т

Таким образом, алгоритм транспонирует матрицу за линейное время, что значительно лучше исходного тривиального алгоритма. Недостатком же алгоритма в том виде, как он изложен, является использование дополнительной памяти. В книге [10] приводится описание алгоритма, лишенного этого недостатка. Основная идея состоит в использовании структур данных матрицы для промежуточных результатов вычислений.

Транспонирование разреженной матрицы

Рис. 7.17. Транспонирование разреженной матрицы

Матричное умножение

Приступим к обсуждению алгоритмов умножения разреженных матриц. Прежде всего рассмотрим алгоритм умножения непосредственно по определению.

Умножение матриц

Рис. 7.18. Умножение матриц

Определение предполагает, что элемент в строке i и столбце j матрицы C вычисляется как скалярное произведение i-ой строки матрицы A и j-го столбца матрицы B. Какие особенности вносит разреженное представление?

N^2

Во-первых, используемая структура данных, построенная на основе формата CRS, предполагает хранение только ненулевых элементов, что усложняет программирование вычисления скалярного произведения, но одновременно уменьшает количество арифметических операций. При вычислении скалярных произведений нет необходимости умножать нули и накапливать полученный нуль в частичную сумму, что положительно влияет на сокращение времени счета. Пусть, например, в первом и втором векторах находится по 1% ненулевых элементов, при этом только десятая часть этих элементов расположена на соответствующих друг другу позициях. В этом случае расчет с использованием информации о структуре векторов может использовать в 1000 раз меньшее число умножений и сложений. Учитывая, что таких пар векторов , получается существенное сокращение объема вычислений. К сожалению, не все так просто. Учет структуры векторов тоже требует машинного времени. Необходимо выполнить сопоставление номеров ненулевых элементов с целью обнаружения пар значений, которые необходимо перемножить и накопить в частичную сумму.

Во-вторых, необходимо научиться выделять вектора в матрицах A и B. В соответствии с определением, речь идет о строках матрицы A и столбцах матрицы B. Выделить строку матрицы в формате CRS не представляет труда: i-я строка может быть легко найдена, так как ссылки на первый элемент poiner[i] и последний элемент (pointer[i+1]-1) известны, что позволяет получить доступ к значениям элементов и номерам столбцов, хранящихся в массивах values и cols соответственно. Таким образом, проход по строке выполняется за время, пропорциональное числу ненулевых элементов в указанной строке, а проход по всем строкам – за время, пропорциональное NZ, где NZ – количество ненулевых элементов в матрице.

Проблема возникает с выделением столбца. Чтобы найти элементы столбца j необходимо просмотреть массив cols ( $\sim NZ$ операций) и выделить все элементы, у которых в соответствующей ячейке массива cols записано число j. Если это нужно проделать для каждого столбца, необходимо $\sim NZ\cdot N$ операций, что является неэффективным.

B^T

Возможное, но не единственное решение проблемы состоит в транспонировании матрицы B: после вычисления появится возможность эффективно работать со столбцами исходной матрицы B, а сам алгоритм транспонирования работает достаточно быстро. Другой вариант – для каждой CRSматрицы, которая может понадобиться в столбцовом представлении, дополнительно хранить транспонированный портрет. Сэкономив время на транспонировании, придется смириться с расходами времени на поддержание дополнительного портрета в актуальном состоянии. Есть и другие варианты решения проблемы, никак не связанные с транспонированием матрицы B (см., например, эффективный алгоритм Густавсона [11]). В дальнейшем мы будем использовать подход, основанный на транспонировании матрицы B.

B^T

В-третьих, необходимо научиться записывать посчитанные элементы в матрицу C. Учитывая, что C хранится в формате CRS, важно избежать перепаковок. Для этого нужно обеспечить пополнение матрицы C ненулевыми элементами последовательно, по строкам – слева направо, сверху вниз. Это означает, что часто используемый в вычислениях на бумаге метод «возьмем первый столбец матрицы B, запишем его над матрицей A и перемножим на все строки…» не подходит. Нужно делать наоборот – брать первую строку матрицы A и умножать ее по очереди на все столбцы матрицы B (строки матрицы ). В этом случае обеспечивается последовательное пополнение матрицы С, позволяющее дописывать элементы в массивы values и cols, а также формировать массив pointer.

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

B^T

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

Еще один непростой момент: в процессе вычисления скалярного произведения на бумаге (в точной арифметике) может получиться нуль, причем не только в том случае, когда при сопоставлении векторов соответствующих друг другу пар не обнаружилось, но и просто как естественный результат. В арифметике с плавающей точкой нуль может получиться еще и в связи с ограничениями на представление чисел и погрешностью вычислений (см. стандарт IEEE-754). Нули, получившиеся в процессе вычислений, можно как хранить в матрице C (в векторе values в соответствующей позиции), так и не хранить; оба подхода имеют право на существование.

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

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

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

Встать на начало обоих векторов $\left(ks =\cdots , ls =\cdots \right)$
Сравнить текущие элементы . Если значения совпадают, накопить в сумму произведение A.Value[ks] * B.Value[ls] и увеличить оба индекса, в противном случае – увеличить один из индексов, в зависимости от того, какое значение больше (например, $A.Col[ks] >B.Col[ls] \rightarrow ls++» />).</li> </ol> Шаг 2 выполняется до тех пор, пока не кончатся элементы хотя бы в одном из векторов. Данный алгоритм имеет линейную трудоемкость, но он не очень хорошо соответствует архитектуре компьютера, т.к. содержит очень много ветвлений. Существуют более эффективные схемы: например, в книге [10] предлагается следующий алгоритм вычисления скалярного произведения разреженных векторов. <ol> <li>Создадим дополнительный целочисленный массив X длины N. Инициализируем его числом -1. Обнулим вещественную переменную S.</li> <li>Просмотрим все ненулевые элементы первого вектора <img decoding=$ . Пусть такой элемент с порядковым номером i расположен в столбце с номером . В этом случае запишем i в j-ю ячейку массива X.
Просмотрим в цикле все ненулевые элементы второго вектора . Пусть элемент с порядковым номером k расположен в столбце с номером . Проверим значение . Если оно равно -1, в первом векторе нет соответствующего элемента, т.е. умножение выполнять не нужно. Иначе умножаем и накапливаем сумму в S.

Данный алгоритм также имеет линейную (относительно числа ненулевых элементов в векторе) трудоемкость, и не содержит избыточных ветвлений.

В заключение отметим, что рассмотренные нами базовые алгоритмы демонстрируют типичные трудности, которые возникают при реализации очевидных «плотных» алгоритмов в «разреженном» случае.

Метод Холецкого для разреженных матриц

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

( 7.41)

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

В качестве примера рассмотрим систему уравнений

$\begin4 &1 & 2& 0.5& 2 \\ 1 & 0.5& 0& 0& 0\\ 2 &0 &3 & 0& 0 \\ 0.5 & 0& 0& 0.625& 0\\ 2 & 0& 0& 0& 16\\ \end \begin x_1\\ x_2\\ x_3\\ x_4\\ x_5\\ \end= \begin 17\\ 3\\ 7\\ 6\\ 12 \end$

( 7.42)

Используя расчетные формулы метода (7.13) и (7.14), можно вычислить фактор Холецкого для данной задачи

$L=\begin 2 &0 & 0& 0& 0 \\ 0.5 & 0.5& 0& 0& 0\\ 1 &-1 &1 & 0& 0 \\ 0.25 & -0.25& -0.5& 0.5& 0\\ 1 & -1& -2& -3& 1\\ \end» /> а затем, выполнив обратный ход алгоритма по формулам (7.15), (7.16), найти искомое решение <img decoding=$

Этот пример иллюстрирует наиболее важный факт, относящийся к применению метода Холецкого для разреженных матриц: матрица обычно претерпевает заполнение. Это значит, что в матрице L появляются ненулевые элементы в позициях, где в нижней треугольной части A стояли нули.

$y_i=x_<i+1> Предположим теперь, что мы перенумеровали переменные в соответствии с правилом , i=1, 2,\cdots , 4, y_5=x_1″ /> и переупорядочили уравнения так, чтобы первое стало последним, а остальные сдвинулись бы на единицу вверх. Тогда мы получим эквивалентную систему уравнений <table border=$ $\begin0.5 &0 & 0& 0 & 1 \\ 0 & 3& 0& 0& 2\\ 0 &0 &0.625 & 0& 0.5 \\ 0 & 0& 0& 16& 2\\ 1 & 2& 0.5& 2& 4\\ \end \begin y_1\\ y_2\\ y_3\\ y_4\\ y_5\\ \end= \begin 3\\ 7\\ 6\\ 12\\ 17 \end$ ( 7.43)

Очевидно, что эта перенумерация переменных и переупорядочение уравнений равносильны симметричной перестановке строк и столбцов матрицы A, причем та же перестановка применяется и к вектору b. Решение новой системы есть переупорядоченный вектор x. Указанную перенумерацию легко записать в виде произведения матицы A на матрицу перестановки P. Напомним, что матрица P называется матрицей перестановки, если в каждой строке и столбце матрицы находится лишь один единичный элемент. В нашем примере

$P=\begin0 &1 & 0& 0& 0 \\ 0 & 0& 1& 0& 0\\ 0 &0 &0 & 1& 0 \\ 0 & 0& 0& 0& 1\\ 1 & 0& 0& 0& 0\\ \end$

( 7.44)

а матрицу системы (7.43) можно получить как

$\overline A=PAP^T$

При этом связь с исходной постановкой задачи выражается соотношением

$\left(PAP^T\right)\left(Px\right)=Pb$

10^<–3> Ниже приведен фактор Холецкого для новой системы уравнений (с точностью до » />) <img decoding=

0.707 &0 & 0& 0& 0 \\ 0 & 1.732& 0& 0& 0\\ 0 &0 &0.79 & 0& 0 \\ 0 & 0& 0& 4& 0\\ 1.414 & 1.154& 0.632& 0.5& 0.129\\ \end» />

Выполнив обратный ход, можно получить переупорядоченный вектор решения

y^T=[2, 1, 8, 0.5, 2]

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

Задача нахождения перестановки, минимизирующей заполнение при вычислении фактора Холецкого, является NPтрудной (т.к. вообще говоря, минимум должен вычисляться по всем n! перестановкам). Поэтому на практике используют эвристические алгоритмы, которые хотя и не гарантируют, что находят оптимальную перестановку, но, как правило, дают приемлемый результат.

Таким образом, при решении разреженной системы методом Холецкого можно выделить следующие этапы:

Переупорядочивание – вычисление матрицы перестановки P;
Символическое разложение – построение портрета матрицы L и выделение памяти под хранение ненулевых элементов;
Численное разложение – вычисление значений матрицы L и размещение их в выделенной памяти.
Обратный ход – решение двух треугольных систем уравнений.

Этапы 1 и 2 специфичны именно для разреженных матриц, тогда как этапы 3 и 4 выполняются для любых задач.

Ниже мы рассмотрим два известных алгоритма определения матрицы перестановки – методы минимальной степени и вложенных сечений.

Метод минимальной степени

Для формулирования метода минимальной степени нам потребуются некоторые сведения из теории графов.

Рассмотрим неориентированный граф G. Формально граф G можно рассматривать как упорядоченную пару двух множеств

$G=\left(V, E \right)$

где $V=<1,\cdots. n>» /> – множество вершин, <img decoding=$ -множество связей между вершинами (ребер). Так как рассматриваемый граф – неориентированный, то множество E удовлетворяет свойству

$\left(e_1,e_2\right)\in E \leftrightarrow \left(e_2,e_1\right)\in E$

$x\in V, y\in V$

Пара вершин называется смежной, если существует ребро, их соединяющее, т.е.

Для любой вершины $x\in V$ можно определить множество $Adj\left(x\right)$ смежных ей вершин как

$Adj\left(x\right)=\lbrace y \in V :\left(x,y\right)\in E\rbrace$

$Adj\left(x\right)$

Мощность множества называется степенью вершины

$deq\left(x\right)=\left|Adj\left(x\right)\right|$

Граф, каждые две вершины которого соединены ребром, называется полным.

Подграф исходного графа – граф, содержащий некое подмножество вершин данного графа и некое подмножество инцидентных им рёбер.

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

Одной из форм задания графа $G=\left(V, E\right)$ является матрица смежности $A=\left(a_<ij>\right)» /> размера <img decoding=$ такая, что

$a_<ij> =\begin\neq 0,\left(i,j\right)\in E,\\ 0,\left(i,j\right)\notin E\end» /> В нашем случае матрицу A системы линейных уравнений (7.41) можно рассмотреть как матрицу смежности соответствующего графа G(A), за исключением диагональных элементов, которым будут соответствовать петли. На рис. 7.19 приведен граф, соответствующий матрице из примера (7.42). <img decoding=$

Рис. 7.19. Пример графа матрицы

После того, как введены все необходимые определения, можно перейти к описанию алгоритма минимальной степени.

Пусть исходной матрице A из (7.41) соответствует граф G(A). Алгоритм минимальной степени строит последовательность графов исключения G_i , каждый их которых получен из предыдущего удалением вершины с минимальной степенью и созданием клики между всеми вершинами, которые были смежными с удаленной. В случае, когда вершин с минимальной степенью несколько, выбирается любая. Алгоритм продолжается до тех пор, пока в очередном графе есть вершины. По мере удаления вершин, их номер записывается в перестановку $\pi$ , по которой впоследствии строится матрица перестановки P.

$\pi$

В качестве примера работы алгоритма на рис. 7.20 приведена последовательность графов исключения и соответствующая перестановка , порождаемая алгоритмом минимальной степени для задачи (7.42).

Последовательность графов исключения

Рис. 7.20. Последовательность графов исключения

$\pi$

Полученная перестановка = соответствует матрице

$P=\begin 0 &1 & 0& 0& 0 \\ 0 & 0& 1& 0& 0\\ 0 &0 &0 & 1& 0 \\ 1 & 0& 0& 0& 0\\ 0 & 0& 0& 0& 1\\ \end» /> <img decoding=$

которая отличается от уже рассмотренной нами перестановки (7.44). Проверим ожидаемые свойства переупорядоченной матрицы: новая матрица системы уравнений, получаемая по формуле , будет

$\overline A=PAP^T=\begin0.5 &0 & 0& 1& 0 \\ 0 & 3& 0& 2& 0\\ 0 &0 &0.625 & 0.5& 0 \\ 1 & 2& 0.5& 4& 2\\ 0 & 0& 0& 2& 16\\ \end$

( 7.45)

10^<–3> Фактор Холецкого, соответствующий данной матрице, будет (с точностью до » />) <img decoding=

0.707 &0 & 0& 0& 0 \\ 0 & 1.732& 0& 0& 0\\ 0 &0 &0.79 & 0& 0 \\ 1.414 & 1.154& 0.632& 0.516& 0\\ 0 & 0& 0& 3.873& 1\\ \end» />

$\overline A$

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

Разреженные матрицы

Матрицы с математическим описанием местоположения элементов

для преобразования индексов матрицы в индекс вектора;
для получения значения элемента матрицы из ее упакованного представления по двум индексам (строка, столбец);
для записи значения элемента матрицы в ее упакованное представление по двум индексам.

Матрицы со случайным расположением элементов

0	6	0	9	0
2	0	7	8	4
10	0	0	0	0
0	12	0	0	0
0	0	3	0	5

Один из основных способов хранения подобных разреженных матриц заключается в запоминании ненулевых элементов в одномерном массиве записей с идентификацией каждого элемента массива индексами строки и столбца матрицы (см. Рисунок 6). Такой способ хранения называется последовательным представлением разреженных матриц. Рисунок 6. Последовательное представление разреженных матриц Доступ к элементу матрицы A с индексами i и j выполняется выборкой индекса i из поля Row, индекса j из поля Colum и значения элемента из поля Value. Следует отметить, что элементы матрицы обязательно запоминаются в порядке возрастания номеров строк для ускорения поиска. Такое представление матрицы A сокращает используемый объем памяти. Для больших матриц экономия памяти является очень актуальной проблемой. Однако последовательное представление разреженных матриц имеет определенные недостатки. Так включение и исключение новых элементов матрицы вызывает необходимость перемещения большого числа существующих элементов. Если включение новых элементов и их исключение осуществляется часто, то можно использовать описываемый ниже метод связанных структур. Метод связанных структур переводит статическую структуру матрицы в динамическую. Эта динамическая структура реализована в виде циклических списков. Для такого представления разреженных матриц требуется следующая структура элемента списка: type PElement = ^TypeElement; TypeElement = record Left: PElement; Up: PElement; Value: TypeData; Row: integer; Colum: integer; end; Рисунок 7 показывает связанную структуру для разреженной матрицы A.Рисунок 7. Представление разреженных матриц в виде связанных структур Циклический список представляет отдельную строку или столбец. Список столбца может содержать общие элементы с одним или более списком строки. Для того чтобы обеспечить использование более эффективных алгоритмов включения и исключения элементов, все списки строк и столбцов имеют головные элементы. Головной элемент каждого списка строки содержит ноль в поле Colum. Аналогично, каждый головной элемент в списке столбца имеет ноль в поле Row. Строка или столбец, содержащие только нулевые элементы, представлены головными вершинами, у которых поле Left или Up указывает само на себя. Может показаться странным, что указатели в этой многосвязной структуре направлены вверх и влево, вследствие чего при сканировании циклического списка элементы матрицы встречаются в порядке убывания номеров строк и столбцов. Такой метод представления используется для упрощения включения новых элементов в структуру. Предполагается, что новые элементы, которые должны быть добавлены к матрице, обычно располагаются в порядке убывания индексов строк и индексов столбцов. Если это так, то новый элемент всегда добавляется после головного и не требуется никакого просмотра списка.

Стек

Стек – это структура данных, в которой новый элемент всегда записывается в ее начало (вершину) и очередной читаемый элемент также всегда выбирается из ее начала. Здесь используется принцип «последним пришел – первым вышел» (LIFO: Last Input – First Output). Стек можно реализовывать как статическую структуру данных в виде одномерного массива, а можно как динамическую структуру – в виде линейного списка. При реализации стека в виде статического массива необходимо резервировать массив, длина которого равна максимально возможной глубине стека, что приводит к неэффективному использованию памяти. Одновременно, работать с такой реализацией проще и быстрее. При такой реализации дно стека будет располагаться в первом элементе массива, а рост стека будет осуществляться в сторону увеличения индексов. Одновременно, необходимо отдельно хранить значение индекса элемента массива, являющегося вершиной стека. Можно обойтись без отдельного хранения индекса, если в качестве вершины стека всегда использовать первый элемент массива, но в этом случае, при записи или чтении из стека, необходимо будет осуществлять сдвиг всех остальных элементов, что приводит к дополнительным затратам вычислительных ресурсов. Стек как динамическую структуру данных легко организовать на основе линейного списка. Поскольку работа всегда идет с заголовком стека, то есть не требуется осуществлять просмотр элементов, удалению и вставку элементов в середину или конец списка, то достаточно использовать экономичный по памяти линейный однонаправленный список. Для такого списка достаточно хранить указатель вершины стека, который указывает на первый элемент списка. Если стек пуст, то списка не существует и указатель принимает значение nil. Рисунок 8. Стек и его организация Поскольку стек, по своей сути, является структурой с изменяемым количеством элементов, то основное внимание уделим динамической реализации стека. Как уже говорилось выше, для такой реализации целесообразно использовать линейный однонаправленный список. Поэтому, при описании динамической реализации будем использовать определения и операции, приведенные в п. 2.2.6.1. Описание элементов стека аналогично описанию элементов линейного однонаправленного списка, где DataType является типом элементов стека. Поэтому здесь приводить его не будем. Основные операции, производимые со стеком:

записать (положить в стек);
прочитать (снять со стека);
очистить стек;
проверка пустоты стека.

Реализацию этих операций приведем в виде соответствующих процедур, которые, в свою очередь, используют процедуры операций с линейным однонаправленным списком. procedure PushStack(NewElem: TypeData; var ptrStack: PElement); begin InsFirst_LineSingleList(NewElem, ptrStack); end; procedure PopStack(var NewElem: TypeData, ptrStack: PElement); begin if ptrStack <> nil then begin NewElem := ptrStack^.Data; Del_LineSingleList(ptrStack, ptrStack); end; end; procedure ClearStack(var ptrStack: PElement); begin while ptrStack <> nil do Del_LineSingleList(ptrStack, ptrStack); end; function EmptyStack(var ptrStack: PElement): boolean; begin if ptrStack = nil then EmptyStack := true else EmptyStack := false; end;

Разреженные матрицы: как ученые ускорили машинное обучение на GPU

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

/ фото alantankenghoe CC

Трудности тренировки крупных нейронных сетей на GPU

Графические процессоры (GPU) лучше подходят для машинного обучения, чем центральные процессоры (CPU). Технические особенности помогают GPU выполнять одновременно множество матричных операций, которые используются для обучения нейронных сетей.

Чтобы добиться схожего результата на центральном процессоре, придется выстроить инфраструктуру из нескольких кластеров CPU, что очень дорого. Система Google для тренировки нейросетей на CPU стоила порядка 5 млрд долларов. Сегодня ученые из Стэнфорда построили систему с аналогичной вычислительной мощностью на GPU всего за 33 тыс. долларов.

Однако здесь есть трудности: использовать весь потенциал GPU на ресурсоемких задачах не так просто. Для обработки данные должны храниться в памяти GPU, однако её объем невелик, что затрудняет тренировку крупных моделей. Например, модель VGG-16 требует около 14 ГБ, в то время как объем памяти Nvidia Titan X составляет 12 ГБ. И эту карту Nvidia позиционирует как один из самых мощных GPU для глубокого обучения.

Как верно заметил EvilGenius18 в комментариях, 7 декабря компания Nvidia представила новую карту Titan V на архитектуре Volta. Она обладает вычислительной мощностью 110 TFLOPS на задачах глубокого обучения, что в 9 раз больше, чем у предшественницы.

При этом для эффективной тренировки крупных моделей нейронных сетей используют различные подходы. Один из них — обработка данных на графическом процессоре последовательными партиями, когда CPU выступает временным контейнером. Минус такого подхода — расходование ресурсов на перенос данных.

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

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

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

Решение от OpenAI

Команда OpenAI разработала программное обеспечение, которое моделирует работу крошечных ядер, способных взаимодействовать с такими матрицами. Ядра опробовали на обучении сетей, анализирующих обзоры на сайтах Amazon и IMDB. Как сообщает команда, уровень ошибок в работе со сводом данных IMDB был снижен с 5,91% до 5,01%.

Ядра реализованы с использованием CUDA, программно-аппаратной архитектуры параллельных вычислений от Nvidia. Но модель OpenAI пока доступна только для TensorFlow. Скотт Грей (Scott Gray), член команды Open AI, сказал, что решение может быть распространено на другие архитектуры, кроме Google TPU2. Компания Nvidia уже знает о работе OpenAI и готова оптимизировать свои системы.

Альтернативные проекты

Концепция разреженных матриц получила свое воплощение в компиляторе с открытым исходным кодом под названием Taco. О проекте, над которым работает команда ученых из Массачусетского технологического института в партнерстве с Adobe Research, стало известно в ноябре. Разработчики искали способ автоматизировать процесс обработки чисел в разреженных матрицах. И использовали для этого тензоры.

О своих разработках в области машинного обучения в декабре отчиталась и компания IBM. Решение ИТ-гиганта — DuHL — предлагает новый метод переноса данных с CPU на GPU. Основная задача технологии — определить, какая информация наиболее важна для алгоритма обучения, и передать её сети в правильном порядке. Исследования показали, что новый подход на основе DuHL в 10 раз быстрее по сравнению с классическим методом последовательной передачи данных между процессорами. Следующая цель компании — предложить DuHL как услугу в облаке.

Но в IBM не первыми придумали переносить GPU-вычисления в облако. Подобные проекты, работающие в том числе по модели IaaS, уже известны. Изначально услугу vGPU предоставляла компания Nvidia. Сейчас этим занимаются и AMD, и Intel.

Об OpenAI

OpenAI — это некоммерческая исследовательская организация, основанная главой Tesla Илоном Маском. Она ставит своей задачей продвижение и развитие искусственного интеллекта на благо человечества. Организация плотно сотрудничает с другими учреждениями и исследователями, предоставляя открытый доступ к своим разработкам.

P.S. Еще несколько материалов из нашего корпоративного блога:

Как IaaS-облака помогают в проектах по сервисной разработке: кейс Azoft
BIM-технологии: обзор решени для информационного моделирования
Балансировка нагрузки в облаке IaaS: зачем нужна и как работает

Что такое разреженная матрица

Разреженные матрицы

Функции

Создание

Манипуляция

Переупорядочение алгоритмов

Итеративные методы и предварительные формирователи

Умножение разреженных матриц

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

Транспонирование матрицы

Матричное умножение

Метод Холецкого для разреженных матриц

Метод минимальной степени

Разреженные матрицы

Матрицы с математическим описанием местоположения элементов

Матрицы со случайным расположением элементов

Разреженные матрицы: как ученые ускорили машинное обучение на GPU

Трудности тренировки крупных нейронных сетей на GPU

Решение от OpenAI

Альтернативные проекты

Об OpenAI

Добавить комментарий Отменить ответ