Нижегородский государственный университет им. Н.И.Лобачевского
Опубликован: 02.10.2012 | Доступ: свободный | Студентов: 1633 / 164 | Длительность: 17:47:00
Специальности: Программист
Лекция 7:

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

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

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

Прямые методы решения СЛАУ

Методы решения систем линейных алгебраических уравнений (СЛАУ) относятся к численным методам алгебры. При формальном подходе решение подобных задач не встречает затруднений: решение системы можно найти, раскрыв определители в формуле Крамера. Однако при непосредственном раскрытии определителей решение системы с n неизвестными требует O\left(n!\right) арифметических операций; уже при n порядка 20 такое число операций недоступно для современных компьютеров. При сколько-нибудь больших n применение методов с таким порядком числа операций будет невозможно и в обозримом будущем. Другой причиной, по которой этот классический способ неприменим даже при малых n, является сильное влияние на окончательный результат округлений при вычислениях.

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

Например, в 80-х годах прошлого века точные методы применялись для решения систем до порядка 104, итерационные – до порядка 107, в 90-х – до порядков 105 и 108 соответственно. Современные суперкомпьютеры способны использовать точные методы при решении еще больших систем.

Мы будем рассматривать систему из n линейных алгебраических уравнений вида


\begin{array}{cccc}
a_{11}x_1 & +a_{12}x_2 & +\ldots +a_{1n}x_n & =b_1\\
a_{21}x_1 & +a_{21}x_2 & +\ldots +a_{2n}x_n & =b_2\\
\ldots\\
a_{n1}x_1 & +a_{n2}x_2 & +\ldots +a_{nn}x_n & =b_n\\
\end{array}
( 7.1)

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

Ax=b ( 7.2)

где A=\left(a_{ij}\right) есть вещественная матрица размера n\times n; b и x - вектора из n элементов.

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

Метод исключения Гаусса

В первую очередь рассмотрим алгоритмыы, предназначенные для решения системы

Ax=b ( 7.3)

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

Метод Гаусса основывается на возможности выполнения преобразований линейных уравнений, которые не меняют при этом решение рассматриваемой системы (такие преобразования носят наименование эквивалентных). К числу таких преобразований относятся:

  • умножение любого из уравнений на ненулевую константу,
  • перестановка уравнений,
  • прибавление к уравнению любого другого уравнения системы.

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

Последовательный алгоритм

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

На итерации i1\leq i<n, метода производится исключение неизвестной i для всех уравнений с номерами k, больших i(т.е.i<k\leq n) Для этого из этих уравнений осуществляется вычитание строки i, умноженной на константу \frac{a_{ki}}{a_{ii}} с тем, чтобы результирующий коэффициент при неизвестной x_{i} в строках оказался нулевым – все необходимые вычисления могут быть определены при помощи соотношений:

a'_{kj} = a_{kj}-\mu_{ki}a_{ij},\\ b'_{k}=b_{k}-\mu_{ki}b_{i},\\ 
i\leq j\leq n, i<k\leq n, 1\leq i<n, ( 7.4)

где \mu_{ki}=\frac{a_{ki}}{a_{ii}} - множители Гаусса.

В итоге приходим к системе Ux=c с верхней треугольной матрицей

U =
\left( \begin{array}{cccc}
u_{11} & u_{12} & \ldots &u_{1n}\\
0 & u_{22} & \ldots &u_{2n}\\
 &  & \ldots &\\
0 & 0 & \ldots &u_{nn}\\
\end{array} \right)

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

a_{ii}
имеет малое значение, то деление и умножение строк на этот элемент может приводить к накоплению вычислительной погрешности и вычислительной неустойчивости алгоритма.

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

y=\max_{i\leq k\leq n}\left|a_{ki}\right|

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

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

x_n
, после этого из предпоследнего уравнения становится возможным определение переменной
x_{n-1}
и т.д. В общем виде выполняемые вычисления при обратном ходе метода Гаусса могут быть представлены при помощи соотношений:

x_i=\frac{\left(c_i-\sum\limits^{n}_{j=i+1} u_{ij}x_{j}\right)}{u_{ii}}, i=n,\ldots ,1

Оценим трудоемкость метода Гаусса. При выполнении прямого хода число операций составит

\sum\limits^{n-1}_{i=1}2\left(n-i\right)^2=\frac{n\left(n-1\right)\left(2n-1\right)}{3}=\frac {2}{3}n^3+O\left(n^2\right)

Для выполнения обратного хода потребуется

\sum\limits^{n-1}_{i=1}\left(n-i\right)=\frac{n\left(n-1\right)}{2}=O\left(n^2\right)

Таким образом, общее время выполнения метода Гаусса при больших n можно оценить как

T_1=\frac{2}{3}n^3\tau

где \tau время выполнения одной операции.

Параллельный алгоритм

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

Для выполнения прямого хода метода Гаусса необходимо осуществить \left(n–1\right) итерацию по исключению неизвестных для преобразования матрицы коэффициентов A к верхнему треугольному виду. Выполнение итерации i, 1\leq i\leq n, прямого хода метода Гаусса включает ряд последовательных действий. Прежде всего, в самом начале итерации необходимо выбрать ведущую строку, которая при использовании метода главных элементов определяется поиском строки с наибольшим по абсолютной величине значением среди элементов столбца i, соответствующего исключаемой переменной x_i. Зная ведущую строку, подзадачи выполняют вычитание строк, обеспечивая тем самым исключение соответствующей неизвестной x_i.

При выполнении обратного хода метода Гаусса подзадачи выполняют необходимые вычисления для нахождения значения неизвестных. Как только какая-либо подзадача i, 1\leq i\leq n, определяет значение своей переменной x_i, это значение должно быть использовано всеми подзадачам с номерами k, k<i: подзадачи подставляют полученное значение новой неизвестной и выполняют корректировку значений для элементов вектора b.

Выделенные базовые подзадачи характеризуются одинаковой вычислительной трудоемкостью. Однако размер матрицы, описывающей систему линейных уравнений, является существенно большим, чем число потоков в программе (т.е.,p\ll n), и базовые подзадачи можно укрупнить, объединив в рамках одной подзадачи несколько строк матрицы. При этом применение последовательной схемы разделения данных для параллельного решения систем линейных уравнений приведет к неравномерной вычислительной нагрузке между потоками: по мере исключения (на прямом ходе) или определения (на обратном ходе) неизвестных в методе Гаусса для большей части потоков все необходимые вычисления будут завершены и они окажутся простаивающими. Возможное решение проблемы балансировки вычислений может состоять в использовании ленточной циклической схемы для распределения данных между укрупненными подзадачами. В этом случае матрица A делится на наборы (полосы) строк вида (см. рис. 7.1).

A=\left(A_1, A_2,\ldots ,A_p\right)^T\\
A_i=\left(a_{i_1},a_{i_2},\ldots ,a_{i_k}\right)\\
i_j=i+jp, 1\leq j\leq k, k=\frac{n}{p}
Ленточная схема

Рис. 7.1. Ленточная схема

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

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

  • поиск ведущей строки,
  • вычитание ведущей строки из всех строк, подлежащих обработке,
  • выполнение обратного хода.

Оценим трудоемкость рассмотренного параллельного варианта метода Гаусса. Пусть n есть порядок решаемой системы линейных уравнений, а p, p<n, обозначает число потоков. При разработке параллельного алгоритма все вычислительные операции, выполняемые алгоритмом Гаусса, были распределены между потоками параллельной программы. Следовательно, время, необходимое для выполнения вычислений на этапе прямого хода, можно оценить как

T_p=\frac{T_1}{p}

Подставив выражение T_1 получим, что время выполнения вычислений для параллельного варианта метода Гаусса описывается выражением:

T_p=\frac{2n^3\tau}{3p}

Теперь можно оценить величину накладных расходов, обусловленных организацией и закрытием параллельных секций. Пусть \delta – время, необходимое на организацию и закрытие параллельной секции. Параллельная секция создается при каждом выборе ведущей строки, при выполнении вычитания ведущей строки из остальных строк линейной системы, подлежащих обработке, а также при выполнении каждой итерации обратного хода метода Гаусса. Таким образом, общее число параллельных секций составляет 3\left(n–1\right).

Сводя воедино все полученные оценки можно заключить, что время выполнения параллельного метода Гаусса описывается соотношением

T_p=\frac{2n^3\tau}{3p}+3\left(n-1\right)\delta ( 7.5)
Связь метода Гаусса и LU-разложения

LU-разложение - представление матрицы A в виде

A=LU ( 7.6)

где L - нижняя треугольная матрица с диагональными элементами, равными единице, а U - верхняя треугольная матрица с ненулевыми диагональными элементами. LU-разложение также называют LU-факторизацией. Известно [4], что LU-разложение существует и единственно, если главные миноры матрицы A отличны от нуля.

Алгоритм LU-разложения тесно связан с методом исключения Гаусса. В самом деле, пусть мы решаем систему уравнений вида (7.2). Непосредственно проверяется, что преобразования k-го шага метода Гаусса равносильны домножению системы (7.2) слева на матрицу


M^{\left(k\right)}=\begin{bmatrix}
1&\quad&\quad&\quad&\quad&\quad&\quad\\
\quad&\ddots&\quad&\quad&\quad&\quad&\quad\\
\quad&\quad&1&\quad&\quad&\quad&\quad\\
\quad&\quad&-\mu_{k+1,k}&1&\quad&\quad&\quad\\
\quad&\quad&-\mu_{k+2,k}&\quad&1&\quad&\quad\\
\quad&\quad&\quad&\quad&\quad&\ddots&\quad\\
\quad&\quad&-\mu_{n,k}&\quad&\quad&\quad&1\\
\end{bmatrix}

где \mu_{ik}- множители Гаусса из (7.4). Как было рассмотрено в п. 7.1.1, прямой ход метода Гаусса преобразует исходную систему уравнений к виду

Ux=c

с верхней треугольной матрицей U. Зная матрицы M^{\left(i\right)}, можно записать матрицу U и вектор c как

\begin{array}{cc}
U=M^{\left(n-1\right)}M^{\left(n-2\right)}\ldots M^{\left(1\right)}A\\
c=M^{\left(n-1\right)}M^{\left(n-2\right)}\ldots M^{\left(1\right)}}b
\end{array}

Обозначим L^{-1}=M^{\left(n-1\right)}M^{\left(n-2\right)}\ldots M^{\left(1\right)} Можно непосредственно проверить, что

L=\begin{bmatrix}
1&\quad&\quad&\quad&\quad\\
\mu_{21}&1&\quad&\quad&\quad\\
\vdots&\vdots&\ddots&\quad&\quad\\
\mu_{n-1,1}&\mu_{n-1,2}&\ldots&1&\quad\\
\mu_{n,1}&\mu_{n,2}&\ldots&\mu_{n,n-1}&1
\end{bmatrix}

Отсюда получаем A=LU.

Таким образом, матрицу L можно получить как нижнюю треугольную матрицу коэффициентов Гаусса, а матрицу U - как верхнюю треугольную матрицу, получаемую в результате работы метода Гаусса. При этом очевидно, что трудоемкость получения LU-факторизации будет такой же-\frac{2}{3n^3}+O\left(n^2\right).

Рассмотренный нами алгоритм LU-факторизации реализован с помощью исключения по столбцу. Следует отметить, что можно сформулировать аналогичный алгоритм, основанный на исключении по строке. В самом деле, основная идея алгоритма с помощью исключения по столбцу заключается в том, что на i-й итерации ведущая строка с подходящими множителями вычитается из строк, лежащих ниже, чтобы занулить все элементы матрицы, расположенные в i-м столбце ниже диагонали. Между тем возможно и другое: на каждой i-й итерации можно вычитать из i-й строки все строки, расположенные выше, умноженные на подходящие коэффициенты, так, чтобы занулить все элементы i-й строки левее диагонали. При этом элементы матрицы L ниже главной диагонали и элементы матрицы U на главной диагонали и выше нее можно вычислять на месте матрицы А. Как и в случае исключения по столбцу, приведенная схема требует проведения \frac{2}{3n^3}+O\left(n^2\right) операций

Рассмотрим теперь еще один способ LU-факторизаци, называемый компактной схемой. Пусть матрица A\left(n\times n\right) допускает LU-разложение (7.6), где

L=\begin{bmatrix}
1&0&\ldots&0\\
l_{21}&1&\ldots&0\\
.&.&\ldots&.\\
l_{n1}&l_{n2}&\ldots&1\\
\end{bmatrix}
U=\begin{bmatrix}
u_{11}&u_{12}&\ldots&u_{1n}\\
0&u_{22}&\ldots&u_{2n}\\
.&.&\ldots&.\\
0&0&\ldots&u_{nn}\\
\end{bmatrix}

т.е.l_{ii}=1, l_{ij}=0 при i<j а u_{ij}=0 при j<i. Из соотношения (7.6) следует, что

a_{ij}=\sum\limits_{k=1}^{n}l_{ik}u_{kj},i,j=\overline{1,n}

Преобразуем эту сумму двумя способами:

\sum\limits_{k=1}^{n}l_{ik}u_{kj}=\sum\limits_{k=1}^{i-1}l_{ik}u_{kj}+l_{ii}u_{ij}+\sum\limits_{k=i+1}^{n}l_{ik}u_{kj}=\sum\limits_{k=1}^{i-1}l_{ik}u_{kj}+l_{ii}u_{ij}
\sum\limits_{k=1}^{n}l_{ik}u_{kj}=\sum\limits_{k=1}^{i-1}l_{ik}u_{kj}+l_{ij}u_{jj}+\sum\limits_{k=j+1}^{n}l_{ik}u_{kj}=\sum\limits_{k=1}^{j-1}l_{ik}u_{kj}+l_{ij}u_{jj}

Отсюда находим

l_{11}=1, u_{11}=a_{11}
\begin{array}{cc}
u_{ij}=a_{ij}-\sum\limits_{k=1}^{i-1}l_{ik}u_{kj}, i\leq j,i=\overline{1,j-1},j=\overline{2,n},\\
l_{ij}=\frac{\left(a_{ij}-\sum\limits_{k=1}^{i-1}l_{ik}u_{kj}\right)}{u_{jj}},i>j,j=\overline{1,n-1},i=\overline{2,n}
\end{array}

Оценка числа операций данного алгоритма LUфакторизаци также составляет

\frac{2}{3}n^3+O\left(n^2\right)

Если разложение (7.6) получено, то решение системы (7.2) сводится к последовательному решению двух систем уравнений с треугольными матрицами (обратный ход)

Ly=b, Ux=y ( 7.7)

Обратный ход требует O\left(n^2\right) операций.

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

Дмитрий Остапенко
Дмитрий Остапенко

поддерживаю выше заданые вопросы

 

Павел Каширин
Павел Каширин

Скачал архив и незнаю как ничать изучать материал. Видео не воспроизводится (скачено очень много кодеков, различных плееров -- никакого эффекта. Максимум видно часть изображения без звука). При старте ReplayMeeting и Start в браузерах google chrome, ie возникает script error с невнятным описанием. В firefox ситуация еще интереснее. Выводится: 

Meet Now: Кукаева Светлана Александровна. 

Meeting Start Time: 09.10.2012, 16:58:04
Meeting Stop Time: 09.10.2012, 18:45:18
Recording Duration:01:47:14

Downloading...

Your Web browser is not configured to play Windows Media audio/video files.

Make sure the features are enabled and available.

 

Олег Страхов
Олег Страхов
Россия
Дмитрий Кифель
Дмитрий Кифель
Казахстан, Темиртау