Решение систем линейных уравнений с разреженными матрицами

Решение систем линейных уравнений с разреженными матрицами

Во многих областях науки и техники, при решении различных задач, возникает необходимость расчета систем линейных уравнений, которые иногда являются сильно разреженными. Поскольку обычно такие системы уравнений имеют большой порядок, то найти их решение (в сколь-нибудь обозримые сроки) возможно только с применением вычислительной техники. Для решения подобных задач можно использовать прямые методы (Гаусса, Крамера), и приближенные (Зейделя, градиентного спуска и др.). Прямые методы дают точное решение (ошибка может быть только из-за выхода величин при преобразованиях за пределы точности типов для используемой платформы), но обычно требуют большее количество итераций. Приближенные методы обычно работают быстрее, но могут требовать особых начальных условий или какой-либо подготовки исходных данных, и решение не всегда сходится. В данной статье рассматривается реализация метода Гаусса, модифицированная для решения разреженных систем линейных уравнений. Данный метод не является «математическим», а скорее «алгоритмическим» и показывает, как путем простых приемов можно ускорить расчет в сотни раз.

Условия применения

Рассматриваемый метод подходит для любых невырожденных (имеющих одно решение) систем линейных уравнений. В дальнейшем под словом матрица будет пониматься матрица коэффициентов линейных уравнений. Единственное условие — метод имеет смысл только для сильно разреженных матриц больших размерностей N>1000, где не менее двух третей коэффициентов равны 0. Для обычных матриц метод не даст выигрыша ни в скорости, ни в экономии памяти.

Описание метода

Метод Гаусса заключается в обнулении нижней (или верхней) полуматрицы a[i,k], где строка i=2..N, столбец k=1..i-1. Используя известное правило, что решение системы не изменится, если к одному из ее уравнений добавить другое, умноженное на произвольный коэффициент, строки матрицы поочередно, начиная с первой, домножаются на специальные коэффициенты и суммируются со всеми ниже лежащими строками. Коэффициенты подбираются таким образом, чтобы каждая строка обнуляла соответствующий ей (i=k) ниже лежащий столбец (только его поддиагональную часть). Так первая строка обнуляет первый столбец A[2..N,1], вторая – второй A[3..N,2] и т.д. до предпоследней строки. Коэффициенты рассчитываются так: например для того чтобы обнулить первый столбец нужно первую строку домножить на k=-A[2,1]/A[1,1] и прибавить ко второй, затем на k=-A[3,1]/A[1,1] и прибавить к третей и т.д. После операция повторяется со второй строкой (она домножается на соответствующие коэффициенты и складывается со всеми нижележащими строками) и т.д. до предпоследней. В результате получается матрица, под главной диагональю которой остаются одни нули. Теперь в последнем уравнении остался только один коэффициент – соответственно можно легко вычислить значение последней неизвестной: X[N]=B[N]/A[N,N], где B – вектор свободных членов уравнений, который кстати, также участвует в преобразованиях при обнулении матрицы A. В предпоследнем уравнении осталось два коэффициента, и последняя неизвестная уже известна, можно вычислить предпоследнюю: X[N-1]=(B[N-1]-A[N-1,N]*X[N])/A[N-1,N-1]. Таким образом, используя обратный ход, можно вычислить все неизвестные. В процессе преобразования матрицы не должны обнуляться коэффициенты главной диагонали. Если такое происходит, то матрица не была факторизована или неправильно составлены уравнения. Или же данная система имеет не одно решение. Более подробно о методе Гаусса можно узнать здесь.

Приведенный выше алгоритм легко реализуем. Количество итераций известно заранее и зависит только от порядка матрицы. Если его запрограммировать в классическом стиле, то он легко будет справляться с небольшими матрицами N

1000. Однако матрицы иногда бывают порядка 100 000 и более. На расчет такой матрицы классическим методом Гаусса уйдут годы, поскольку количество итераций будет огромным. Если данную матрицу хранить в обычном двумерном массиве, то при использовании типа double она будет занимать в памяти 100 000 х 100 000 х 8 = 800 000 000 байт = 800 Мб что совсем не мало. А если учесть что размер матрицы возрастает пропорционально квадрату порядка матрицы, то для матрицы порядка 400 000 это будет в 16 раз больше, т.е. 12.8 Гб!

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

• Более компактное хранение информации.
• Уменьшение количества итераций.

Для разреженных матриц возможно и то, и другое:

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

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

• для главной диагонали это одна координата d, которая равна соответствующему индексу элемента главной диагонали
• для верхней и нижней полуматрицы координаты d,j, где j – смещение от главной диагонали вправо или вниз (как показано на рисунке)

Диаграмма классов

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

Реализация метода

Все исходные коды приведены на C#, и представляют полный набор классов, необходимых для функционирования метода. Код написан максимально просто для понимания принципа работы метода, и поэтому не претендует на максимальную производительность и не является оптимальным. О необходимых доработках сказано в заключении. Код снабжен комментариями и не нуждается в дополнительном разъяснении.

Реализация интерфейса строк и столбцов:

Реализация интерфейса матрицы:

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

На данном этапе нужно решить в каком формате хранить значения ячеек матрицы. Для очень сильно разреженных матриц, когда в каждой строке/столбце матрицы ненулевых значений совсем мало, данные лучше хранить в обычных массивах, так как поиск в небольшом массиве работает быстро. Если же матрица не сильно разрежена, то лучше использовать хэш-таблицу. Ниже приведены реализации интерфейса IRow для обоих случаев:

Реализация интерфейса строки c использованием массива для хранения значений

IРеализация интерфейса строки c использованием хэш-таблицы для хранения значений:

Полная реализация класса «матрица»:

И, наконец, реализация метода Гаусса:

Заключение

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

• Как минимум нужно сделать обработку исключительных ситуаций, например таких как деление на ноль или выход индекса за пределы массива
• Реализовать расчетный алгоритм в отдельном потоке команд. Если вызывать метод calculate из основного UI потока, это остановит прием сообщений окружения («повесит программу») до момента завершения расчета
• Если необходимо рассчитать несколько систем линейных уравнений, которые отличаются только правой частью (свободными членами), то естественно нет необходимости каждый раз обнулять нижнюю полуматрицу. Для этого нужно доработать класс Gauss таким образом, чтобы он, обнуляя матрицу, сохранял историю изменений. После этого для расчета неизвестных нужно будет преобразовать только вектор свободных членов, передаваемых в метод calculate в соответствии с историей. Историю можно хранить в отдельной матрице, или в нижней полуматрице основной матрицы, которая все равно не используется при обратном ходе
• Если данный алгоритм планируется реализовать на платформе, где возможна прямая работа с указателями (Delphi, C++), то лучше по возможности использовать их, это позволит в некоторых местах выполнять меньшее количество преобразований
• Для некоторых матриц необходима доработка алгоритма факторизации, для исключения попадания нулей на главную диагональ
• В реализациях интерфейса IRow, при создании экземпляров массивов и хэш-таблиц, в явном виде не задавалось свойство capacity. На практике, чтобы оптимально использовать быстродействие и экономить память, нужно подобрать компромисное значение capacity. Так с увеличением capacity увеличивается быстродействие и расход памяти, с уменьшением – наоборот. Оптимальное значение capacity зависит от структуры матрицы и не превышает максимального количества ненулевых элементов в строке или столбце матрицы

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

Литература

1. Эндрю Троелсен. С# и платформа .NET. Библиотека программиста. — СПб.: Питер, 2004.

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

ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ ОБРАЗОВАТЕЛЬНОЕ

УЧРЕЖДЕНИЕ ВЫСШЕГО ПРОФЕССИОНАЛЬНОГO

“ЮЖНЫЙ ФЕДЕРАЛЬНЫЙ УНИВЕРСИТЕТ”

РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ

АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ ДЛЯ БОЛЬШИХ РАЗРЕЖЕННЫХ МАТРИЦ

В.А.Еремеев. Решение систем линейных алгебраических уравнений для больших разреженных матриц: Учебно-методическое пособие. – г.Ростов-на-Дону, 2008 г. – 39 с.

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

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

2. векторные и матричные нормы;

3. итерационные методы;

4. методы подпространств Крылова.

Пособие предназначено для преподавания дисциплин в рамках магистерской образовательной программы “Математическое моделирование и компьютерная механика”.

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

Пособие подготовлено в рамках гранта ЮФУ K-07-T-41.

1 Задачи линейной алгебры с разреженными матрицами на примере дискретизации уравнения Пуассона 4

1.1 Решение уравнения −x 00 = f (t ) . . . . . . . . . . . . . . . 4

1.2 Решение уравнения Пуассона . . . . . . . . . . . . . . . . 6

2 Векторные и матричные нормы 11

2.1 Векторные нормы . . . . . . . . . . . . . . . . . . . . . . 11

2.2 Матричные нормы . . . . . . . . . . . . . . . . . . . . . . 13

2.3 Связь векторных и матричных норм . . . . . . . . . . . 15

3 Итерационные методы 19

3.1 Определения и условия сходимости итерационных методов 19

3.2 Метод простой итерации . . . . . . . . . . . . . . . . . . . 23

3.3 Метод Якоби . . . . . . . . . . . . . . . . . . . . . . . . . 24

3.4 Метод Гаусса-Зейделя . . . . . . . . . . . . . . . . . . . . 25

3.5 Метод последовательной верхней релаксации (SOR) . . . 26

3.6 Метод симметричной последовательной верхней релак-

сации (SSOR) . . . . . . . . . . . . . . . . . . . . . . . . . 27

4 Методы подпространств Крылова 30

4.1 Инвариантные подпространства . . . . . . . . . . . . . . 30

4.2 Степенная последовательность и подпространства Кры-

4.3 Итерационные методы подпространств Крылова . . . . . 33

Дополнительная литература 39

1 Задачи линейной алгебры с разреженными матрицами на примере дискретизации уравнения Пуассона

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

1.1 Решение уравнения −x 00 = f (t )

Рассмотрим простейшую краевую задачу

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

так, чтобы t ₀ = 0, t_{n +1} = 1. Используем обозначения

Для дискретизации второй производной воспользуемся центральной конечной разностью

,

которая аппроксимирует x 00 (t_i ) c точностью O (h 2 ). Тогда исходная краевая задача сводится к системе линейных алгебраических уравнений (СЛАУ) относительно

,

.

Эту СЛАУ можно записать также в матричном виде

Обратим внимание на следующие свойства матрицы A .

• A – разреженная матрица, она трехдиагональная;

• A – симметричная и, более того, положительно определенная. Это свойство она унаследовала от свойств оператора исходной краевой задачи;

• Структура A достаточно простая, ее элементы вычисляются по простым формулам. Следовательно, для ее хранения в памяти компьютера нет необходимости использовать типы данных, предназначенные для хранения заполненных матриц.

Обсудим размерность полученной СЛАУ. Предположим, что требуемая точность аппроксимации равна ² = 10 −6 . Поскольку ² ∼ h 2 , то

−3 , т.к. h ∼ √² . Т.е. требуемый шаг сетки нужно выбрать порядка 10 мы имеем матрицу размерности 10 3 × 10 3 . На практике, размерность должна быть выше, чтобы удовлетворить большей точности.

Рис. 1: Конечно-разностная сетка. Естественная нумерация узлов.

1.2 Решение уравнения Пуассона Рассмотрим более сложную краевую задачу

где u = u (x,y ) – неизвестная функция, f (x,y ) – известная, заданные на единичном квадрате (x,y ) ∈ Ω = [0, 1]×[0, 1]>, Γ = ∂ Ω – граница Ω. Введем на Ω равномерную сетку, так что координаты узлов даются формулами x_i = hi, y_i = hi, i = 0. N + 1, h = 1/ (N + 1).

Для дискретизации вторых производных в операторе Лапласа ∆ воспользуемся формулами для центральной конечной разности

,

Таким образом, исходная краевая задача сводится к СЛАУ

−u i −1,j + 4u i,j − u i +1,j − u i,j −1 − u i,j +1 = h 2f i,j (1) относительно u_i,j . Для того, чтобы свести ее к стандартному виду A x = b, необходимо преобразовать u_i,j к выражению с одним индексом, т.е. перенумеровать каким-то образом. Выбор нумерации влияет, вообще говоря, на структуру матрицы A .

Рассмотрим случай N = 3. Соответствующая сетка показана на рис. 1, где использована естественная нумерация внутренних узлов. Полые кружки соответствуют граничным узлам, в которых известно значение функции, а сплошные – внутренним.

Если преобразовать конечно-разностное уравнение (1) в матричное уравнение A x = b, вводя вектор неизвестных по правилу u 1, 1 = x 1, u 1, 2 = x 2, u 1, 3 = x 3, u 2, 1 = x 4. u 3, 3 = x 9,

Видно, что A является симметричной ленточной разреженной матрицей с диагональным преобладанием.

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

Обсудим размерность полученной СЛАУ. Предположим как и ранее, что требуемая точность аппроксимации равна ² = 10 −6 . Поскольку ² ∼ h 2 , то требуемый шаг сетки нужно выбрать порядка 10 −3 , т.к.

√ 3 × 10 3 = 10 6 . Таh ∼ ² . Число узлов здесь оказывается равным 10 кому же числу равно число неизвестных n . Таким образом, мы имеем матрицу размерности 10 6 × 10 6 .

Нетрудно догадаться, что если рассмотреть, не квадрат, а куб, то число неизвестных будет примерно равно n = 10 9 . В расчетной практике встречаются задачи, где нужно определять не одну функцию, а несколько. Например, в гидродинамике при учете теплопереноса имеем четыре скалярных неизвестных (три компоненты поля скорости и поле температуры). Если рассматривать конечно-разностную аппроксимацию этой задачи с точностью ² = 10 −6 , то получим размерность n = 4 · 10 9 .

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

Рассмотрим другую нумерацию узлов, показанную на рис. 2. Здесь использовано так называемое красно-черное разбиение (или чернобелое). Вначале нумеруются узлы, сумма индексов которых – четное число, а потом узлы, сумма индексов которых дает нечетное число. Таким образом, вектор переменных вводится по правилу вектор неизвестных по правилу

– это красные неизвестные, а

Рис. 2: Конечно-разностная сетка. Красно-черное разбиение.

Соответствующая красно-черному разбиению матрица дается формулой

Решение разреженных СЛАУ больших размерностей средствами ManagedCuda в .NET

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

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

Постановка задачи

Рассматривается классическая задача решения СЛАУ:
Ax=b
Здесь матрица A размером nxn, x — вектор неизвестных размером n, b — вектор известных свободных членов размером n.

Автор статьи занимается разработкой специализированного программно-вычислительного комплекса (ПВК) для моделирования и оптимизации нестационарных течений газа в системах газоснабжения. В решаемых мною задачах A является положительно определенным якобианом с диагональным преобладанием некоторой системы нелинейных уравнений. A получается в результате матричных преобразований матрицы инциденций графа системы газоснабжения и других матриц. В моих практических расчетах A обычно заполнена на 6-10%, остальное нули. Размер же ее зависит от размера конкретной системы газоснабжения и варьируется n от 10 до 1000.

Наш ПВК разрабатывается под .NET 4.0 на C#. Расчетный модуль и вся математика также разрабатывается на C#. Изначально для решения СЛАУ я написал собственную, доморощенную, библиотеку, не использующую технологию разреженных матриц. Для решения СЛАУ применял метод LU декомпозиции. И до поры меня все устраивало, пока я не начал заниматься задачей оптимизацией нестационарных режимов течения газа, где приходится методом динамического программирования осуществлять большое количество переборов значений управляющих параметров и, соответственно, много раз решать СЛАУ. Стандартный профилировщик Visual Studio показал, что в процессе выполнения программы около 40% затрат приходится на решение СЛАУ.

Именно в этот момент я понял, что пора что-то менять.

Математическая библиотека Math.Net Numerics

Проведя анализ имеющихся математических библиотек под .Net я решил остановиться на библиотеке Math.Net Numerics. Подробно ознакомиться с ней вы можете здесь.

Меня заинтересовали ее ключевые возможности:

• Реализована технология разреженные матриц и векторов;
• Имеются встроенные все необходимые вектор-матричные операции;
• Присутствуют встроенные решатели;
• Интуитивно понятный API.

Приведу листинг с примером решения СЛАУ средствами Math.Net Numerics:

Как видно, все очень элегантно выглядит, но практически сразу меня настигло разочарование — встроенный решатель Math.Net Numerics работал значительно медленнее моего. Меня это совершенно не устроило, однако претензии к реализованным в библиотеке вектор-матричным операциям не столь существенные. Поэтому полностью отказываться от Math.Net Numerics я не стал, оставив код, где производится работа с векторами и матрицами.

Но тут мне вспомнился весьма удачный опыт коллег по аспирантуре, которые для решения задач подземной гидродинамики использовали графические процессоры. В распоряжении у меня имеется ноутбук с видеокартой nVidia GeeForce GT 540M с 2 GB памяти и поддержкой технологии CUDA. Было принято решение попробовать эту технологию на деле, о чем теперь не жалею.

Библиотека ManagedCuda

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

Меня заинтересовала библиотека cuSPARSE. При первом знакомстве с ней я наткнулся на проблемы:

• Непосредственная работа с библиотекой возможна только через C/C++. Конечно, проблема решится, если использовать качественный враппер, который очень не хотелось писать. По результатам поиска был найден CUDA-враппер для .Net под названием Cudafy.
• cuSPARCE позволяет решать СЛАУ с треугольными матрицами, а в моем случае матрицы таковыми не являются. Однако cuSPARSE содержит методы факторизации матриц, приводящие их к треугольной форме (метод LU факторизации, разложение Холецкого). Однако в API Cudafy отсутствовали соответствующие методы, поэтому от Cudafy пришлось отказаться.

После продолжения поиска я открыл для себя библиотеку ManagedCuda, которая является высокоуровневой оберткой над CUDA API. Все её возможности я перечислять не буду — их можно найти на официальном сайте, но остановлюсь на том, как можно, используя Math.Net Numerics и ManagedCuda, элегантно и эффективно решать СЛАУ.

Идея заключается в использовании SparseMatrix из Math.Net Numerics для формирования и хранения разреженной матрицы в формате CSR, который принимают на вход cuSPARSE и ManagedCuda. Далее приводится листинг соответствующей программы:

Вычислительный эксперимент: анализ временных затрат на решению СЛАУ различными технологиями

Дабы не утомлять читателя сухим текстом и листингами программ, рассмотрим результаты расчетов СЛАУ размерностей от 50 до 500 с помощью различных технологий:

• Доморощенный решатель, написанный автором статьи. Назовем его Mani.Net.
• Стандартный решатель библиотеки Math.Net Numerics.
• Метод LU декомпозиции библиотеки ManagedCuda.

Матрица A заполняется на 10% случайным образом.

Рисунок наглядно демонстрирует, почему мне пришлось отказаться от Math.Net Numerics для решения СЛАУ.
Отметим, при размерности матрицы равной 500 моему собственному решателю (Mani.Net) потребовалось 1170 мс, Math.Net Numerics — 12968 мс, ManagedCuda — 70 мс.

Заключение

Следует ожидать в комментариях замечания, мол де не на всех машинах есть GPU nVidia с поддержкой CUDA. Действительно, это так. Поэтому наше приложение настроено на две конфигурации компиляции: с собственной библиотекой решения СЛАУ и с ManagedCuda.

Насчет Mani.Net отмечу, что это не реклама моей собственной библиотеки. Найти ее нигде невозможно и никому её я передавать не буду. Нет, я не жадный. Мне стыдно за код.

Спасибо за прочтение статьи! Буду рад узнать о вашем мнении и замечаниях в комментариях.