QR разложение прямоугольной матрицы является одним из базовых вычислительных алгоритмов для многих задач вычислительной математики. В частности, подобные вычисления возникают при решении СЛАУ, при решении задач наименьших квадратов и задач на собственные значения [1], и т.д. Возможность эффективно параллельным образом вычислять QR разложение разреженной матрицы в некоторых случаях означает возможность использования новых классов вычислительных алгоритмов, и поэтому подобные разработки представляют конкретный практический интерес.
В современной вычислительной практике существенное использование разреженности матрицы зачастую является практически единственным способом эффективного решения задач большой размерности. Относительно давно научились эффективно использовать разреженность вычислений при параллельном вычислении традиционных полных треугольных разложений типа LU разложения [2] (разложения Холецкого в симметричном случае) и при вычислении их неполных вариантов [3]. Что касается параллельного разреженного QR разложения, то эта тематика получила особое развитие относительно недавно, прежде всего в работах [4], [5] Тима Дэвиса с соавторами. В этих работах предложен мульти-фронтальный алгоритм построения QR разложения, в котором параллельно вычисляемые разреженные профильные блочно строчные QR разложения по дереву вычислений комбинируются группами для формирования объединяющих QR разложений. В данной работе рассматриваются идейно близкие алгоритмы параллельного вычисления разреженного QR разложения матрицы. Основные отличия предлагаемого в работе подхода состоят в следующем:
за счет введения начального нулевого блока вычисления проводятся по возможности с более разреженными Q факторами;
введено дополнительное строчное упорядочивание и биение, которое позволяет дополнительно уменьшить связность вычислений и заполнение Q факторов;
предложен алгоритм построения представления матрицы, удобного для параллельного вычисления разреженного QR разложения, на основе объемного биения расчетной сетки.
Данная работа является теоретической основой для практической работы [6], содержащей описание реализации алгоритма на гибридной параллельной MPI+threads+SIMD архитектуре и начальные результаты экспериментов на такой архитектуре для тестовых задач. Также работа является базовой для планируемой серии работ автора по новым параллельным итерационным алгоритмам решения СЛАУ и задач наименьших квадратов на основе композиции подпространств, порождаемых разреженными базисами.
Работа построена следующим образом. В Разделе 2 описывается блочное преобразование Хаусхолдера [7] и его применение для профильной разреженной QR факторизации прямоугольной матрицы, аналогичной используемым в работах [4], [5]. В Разделе 3 описывается алгоритм вычисления расширенного профильного разреженного QR разложения с более разреженным во многих практических случаях Q фактором для матрицы, расширенной квадратным нулевым начальным блоком. Проводится сравнение профильного и расширенного профильного разреженного QR разложения для некоторых матриц. В Разделе 4 описывается построение QR разложения матрицы в случае строчного биения матрицы. В Разделе 5 вводится понятие верхней квази-треугольной матрицы с разреженностью типа вложенных сечений и приводится параллельный алгоритм вычисления разреженного QR разложения для такой матрицы. В Разделе 6 описывается способ вычисления столбцевого и строчного упорядочиваний и биений, которые приводят матрицу, заданную на объемной сетке с известным объемным геометрических биением, к виду верхней квази-треугольной матрицы с разреженностью типа вложенных сечений.
Классическое преобразование Хаусхолдера -это унитарная матрица вида [1]:
Блочные преобразования Хаусхолдера (2.4) можно применить для вычисления разреженного QR разложения разреженной матрицы , составленной из мелко блочных столбцов. Для этого можно поступить например следующим образом. Для минимизации заполнения R-фактора QR разложения вычисляем структуру разреженности матрицы и вычисляем упорядочивание, минимизирующее заполнение матрицы с такой структурой при ее треугольной факторизации. Можно выбрать, например, упорядочивание вложенных сечений ND или упорядочивание RCM [2]. Это определяет упорядочивание мелко блочных столбцов матрицы . Далее упорядочивание мелко блочных строк можно выбрать последовательно нумеруя сначала элементы первого столбца , потом второго, и т.д. Переупорядоченная таким образом матрица будет иметь вид типа приведенного на Рисунке 1 слева. В качестве базового столбцевого упорядочивания для этой матрицы было выбрано упорядочивание ND.
Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org При построении так называемого профильного разреженного QR разложения таким образом упорядоченной мелко блочной матрицы будем действовать по аналогии со случаем плотной матрицы. По первому мелко блочному столбцу матрицы строим разреженное блочное преобразование Хаусхолдера с разреженностью столбца такой, чтобы обнулить мелкие блоки матрицы под первой блочной диагональю. Применяем транспонированное блочное преобразование Хаусхолдера ко второму мелко блочному столбцу, для полученного результата строим следующее разреженное блочное преобразование Хаусхолдера для обнуления элементов под второй мелко блочной диагональю, и т.д. На Рис. 1 справа показана мелко блочная разреженность Q-и R-факторов QR разложения для матрицы слева. Для этого теста заметно довольно сильное заполнение в Q-факторе внутри соответствующего профиля матрицы.
Профильное разреженное QR разложение матрицы, описанное выше, не является единственным вариантом построения QR разложения. Существует вариант, для которого более естественным образом выделяется его параллельная структура, и который во многих практически важных случаях приводит к меньшему заполнению в факторе Q. Это так называемое расширенное профильное QR разложение.
Схематически профильное и расширенное профильное QR изображены на Рисунке 2. Фактически расширенное профильное QR разложениеэто профильное QR разложение, примененное к матрице, расширенной сверху нулевым квадратным блоком. При этом структура разреженности Q фактора пополняется возможными дополнительными элементами на месте бывшего фактора R, и отсоединенными диагональными элементами, примыкающими к новому R фактору. Дополнительное заполнение блочных преобразований Хаусхолдера Q-фактора в над-диагональной части расширенного профильного QR разложения может быть значительным, как это уже было показано выше. Тем не менее, в частном случае матриц с числом строк значительно превышающем число столбцов этим дополнительным заполнением можно пренебречь на фоне возможного избыточного внутри профильного под-диагонального заполнения из-за взаимодействий преобразований Хаусхолдера на диагональном элементе. В контексте интересующих автора приложений число столбцов всегда будет суммарно в сотни и тысячи раз меньше числа строк, поэтому в последующей части изложения в качестве базового разложения будет использоваться только расширенное профильное QR разложение. Расширенное профильное QR разложение имеет еще одно естественное структурное преимущество. Если число строк в матрице меньше числа столбцов (что редко, но случается в приложениях из-за, например, дополнительного строчного биения матрицы), то в этом случае в профильном разложении нужно делать что-то искусственное, в то время как в расширенном разложении это не приводит к особым ситуациям.
Описанный в предыдущем разделе параллельный алгоритм вычисления QR разложения по блочным строкам можно сделать эффективным за счет использования дополнительной столбцевой разреженности матрицы.
Для начала рассмотрим двух-уровневую организацию вычислений для прямоугольной матрицы со структурой разреженности, изображенной на Рисунке 5 слева. Пусть число мелко блочных столбцов в матрицах , и равны соответственно , и соответственно, . Обобщая этот алгоритм на общий случай введем следующее определение. Определение 1. Прямоугольную мелко блочную матрицу с введенным на ней блочным строчным и столбцевым биениями будем называть верхней квази-треугольной -уровневой матрицей о труктурой разреженно ти типа вложенных ечений, если в терминах крупных блоков матрица является квадратной верхней треугольной и имеет структуру разреженности типа вложенных сечений, описываемой -уровневым деревом зависимостей вычислений.
В частности, матрица на Рисунке 5 слева в терминах Определения 1 является двух уровневой верхней квази-треугольной с двух уровневым бинарным деревом зависимостей вычислений. Аналогично, матрица на Рисунке 5 справа является трех уровневой верхней квази-треугольной с трех уровневым бинарным деревом зависимостей вычислений.
Как следует из предыдущего изложения, для эффективного вычисления QR разложения верхних квази-треугольных матриц с разреженностью типа вложенных сечений можно использовать следующий параллельный алгоритм:
1. Параллельно и независимо для каждой блочной строки матрицы строим расширенное профильное разреженное QR разложение на основе блочных преобразований Хаусхолдера; 2. Параллельно в порядке, определяемом деревом зависимостей вычислений, достраиваем QR разложения для объединяющих подматриц типа (5.4) для вычисления QR разложений соответствующих мелко блочных столбцевых окаймлений.
Существует множество способов, при помощи которых разреженную прямоугольную мелко блочную матрицу можно привести к виду многоуровневой верхней квази-треугольной матрицы типа вложенных сечений. Простейший из них например такой.
По структуре разреженности матрицы вычисляем упорядочивание вложенных сечений. По структуре переупорядоченной матрицы строим блочное биение и бинарное дерево нужной глубины, которое описывает разреженность этой матрицы в контексте вложенных сечений. Переупорядочиваем мелко блочные столбцы матрицы в соответствии с упорядочиванием вложенных сечений. Упорядочивание строк вводим последовательно нумеруя элементы первого мелко блочного столбца , затем второго и т.д. Столбцевое биение вводим в соответствии с биением вложенных сечений. Строчное биение вводим по мере того как заканчиваются столбцевые блоки. В результате получится переупорядоченная по столбцам и строкам матрица с введенным строчным и столбцевым биением типа изображенной на Рисунке 6. В некоторых случаях удобно ввести дополнительное строчное упорядочивание и биение для оптимизации затрат при обработке окаймлений матрицы. Для этого в каждой блочной строке последними ставим строки, в которых есть элементы в последнем окаймлении, и вводим новый строчный блок если такие строки есть; предпоследними ставим строки, которые не были перенесены в конец ранее и которые имеются в предпоследнем окаймлении, и ставим новый строчный блок если такие строки имеются, и т.д. Для тестовой задачи на Рисунке 6 получится новое строчное упорядочивание и строчное биение, изображенное на Рисунке 7. На рисунке видно заметно меньшее блочное заполнение в окаймлениях матрицы. Для задачи в таком дополнительном биении построить QR разложение можно параллельным алгоритмом, аналогичным описанному выше. Дополнительным является этап объединяющего строчного QR разложения в рамках одного узла дерева зависимостей, которое возникает по причине введения дополнительного строчного биения в каждой основной подматрице.
Описанный в предыдущем разделе основанный на явном вычислении упорядочивании ND способ построения упорядочиваний и биений, при которых мелко блочная матрица становится многоуровневой верхней квази-треугольной типа вложенных сечений, не является практически удобным для множества приложений, в которых матрица возникает как следствие аппроксимации уравнений на расчетной сетке. Для подобных задач необходима специальная методика построения соответствующих упорядочиваний и биений. В данном разделе предлагается подобная методика, которая основана на предположении о том, что столбцы и строки матрицы ассоциируются естественным образом с элементами расчетной сетки, и для этой сетки имеется ее геометрически локальное разбиение на достаточно малые объемыдомены.
Таким образом, предположим, что для расчетной области имеется разбиение на домены, и что строки и столбцы матрицы ассоциируются с расчетными доменами. Это означает, что возможно ввести нумерацию доменов, нумерацию строк и столбцов матрицы и строчное и столбцевое биения таким образом, чтобы в общем случае в прямоугольном блоке (1,1) были элементы связей первого домена между собой, в блоке (1,2) элементы связей первого домена со вторым, и т.д.
Для начала рассмотрим вопрос о построении дерева зависимостей вычислений. Для этого введем в рассмотрение две квадратных разреженных целочисленных матрицы размером , здесьполное число доменов. Первую матрицу будем называть матрицей вязей. В этой Упорядочим матрицу в соответствии с введенным биением и упорядоченным квази-бинарным деревом зависимостей следующим образом. Для этого введем соответствие между узлами дерева и введенными новыми строчным и столбцевым блочными биениями. Поместим каждый столбцевой блок в соответствии с его связями следующим образом. Внутренние столбцевые блоки домена поместим в соответствующий домену лист дерева. Столбцевые квази двух мерные границы всех доменов поместим в наименьший по номеру узел упорядоченного дерева зависимостей, в сыновьях которого имеются узлы дерева, содержащие соответствующие родительские для этой квази двух мерной границы 2 домена. Аналогично, столбцевые квази одномерные границы всех доменов поместим в наименьший по номеру узел упорядоченного дерева зависимостей, в сыновьях которого имеются узлы дерева, содержащие соответствующие родительские для этой квази одномерной границы 3 домена. Наконец, остальные столбцевые квази нуль-мерные границы всех доменов поместим в наименьший по номеру узел упорядоченного дерева зависимостей, в сыновьях которого имеются все узлы дерева, содержащие соответствующие родительские для этой квази-нульмерной границы >3 доменов. Строчные подблоки каждого домена разместим в тех узлах упорядоченного дерева зависимостей, в котором они встречаются в первый раз, если блочные столбцы заранее упорядочить в соответствии с описанным выше алгоритмом. Далее, все блочные столбцы в каждом узле упорядочим по их мерностичем больше мерность, тем раньше расположен столбцевой блок. После подобных преобразований очевидно матрица станетуровневой верхней квази-треугольной матрицей типа вложенных сечений, и ее структура разреженности будет иметь вид, подобный изображенному на Рисунке 7. Для уменьшения заполнения Q-и R-факторов QR разложений в столбцевых блоках узлов можно дополнительно использовать упорядочивание ND как было описано выше.
В работе представлен параллельный алгоритм вычисления QR разложения многоуровневой разреженной верхней квази-треугольной матрицы со структурой разреженности типа вложенных сечений. Базовым вычислительным примитивом алгоритма являются блочные разреженные преобразования Хаусхолдера. Предложены варианты алгоритмов, приводящих заданную матрицу к упомянутому виду, удобному для параллельного вычисления ее разреженного QR разложения, в том числе для задач аппроксимации уравнений на расчетной сетке. Результаты численных экспериментов с предложенным алгоритмом для тестовых задач на гибридной параллельной MPI+threads+SIMD архитектуре приведены в работе [6], также представленной в качестве доклада на конференцию.
Keywords: sparse rectangular matrix, upper quasi triangular matrix, nested dissection, volume partitioning, QR decomposition, Householder transformations, parallel algorithm, SLAE, least squares problem The paper considers algorithm for computing sparse QR decomposition of a specially ordered rectangular matrix. Decomposition is based on block sparse Householder transformations. For ordering computations the ND-type ordering for sparsity of ATA matrix can be used, here Aoriginal rectangular matrix. For mesh based problems the ordering can be constructed starting from appropriate volume partitioning of the computational mesh. Parallel computations are based on sparse QR decomposition for sets of rows with additional zero block at the beginning. The suggested algorithm is planned to be used as main computational kernel in the developed by the author parallel iterative algorithms for solving SLAEs and least squares problems. The corresponding algorithms will be based on composition of the subspaces represented by sparse bases.