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

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

7. Итерационные методы решения СЛАУ

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

Ax=b ( 7.49)
система линейных уравнений относительно неизвестного вектора x\in R^n с симметричной положительно определенной матрицей A размерности n \times n и правой частью вектором b\in R^n Точное решение системы (7.49) обозначим через x^*

Итерационный метод , предназначенный для решения системы (7.49), генерирует последовательность векторов x^{(s)} \in R^n, s=0,1,2,..., каждый из которых может рассматриваться исследователем как приближенное решение системы (7.49). Начальное приближение – вектор x^{(0)} \in R^n – задает исследователь. Итерационный метод называется сходящимся , если для любого начального приближения x^{(0)} \in R^n последовательность x^{(s)} \in R^n, s=0,1,2,..., сходится к точному решению x^*, т.е.

\lim_{s\to\infty}\| x^{(s)}-x^*\|=0 ( 7.50)
где через \|\| обозначена любая векторная норма.

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

\|x^{(s)}-x^{(s-1)}\| < \varepsilon_1 ( 7.51)
где x^* – приближение, полученное на итерации с номером s, x^{(s-1)} – приближение, полученное на предыдущей итерации. Параметр точности метода \varepsilon_1 задает исследователь. Если условие (7.51) выполнено, метод завершает работу и вектор x^{(s)} трактуется как приближенное решение системы (7.49). Величина \varepsilon_2 , определяемая как
\varepsilon_2 =\| x^{(s)}- x^{(s-1)}\| ( 7.52)
называется достигнутой точностью метода и сообщается исследователю.

Второй критерий определяется максимальным числом итераций N, на которое готов пойти исследователь. Если номер итерации, которую предстоит выполнить, превышает NN, итерация не выполняется, метод завершает работу, вектор x^{(N)} трактуется как приближенное решение системы и достигнутая точность метода \varepsilon_2 также сообщается исследователю.

При изучении сходимости итерационных методов и при их реализации в качестве векторной нормы \|\| обычно используют евклидову норму \|\|_2, \|\|_1, определяемую суммой модулей всех компонент вектора, норму \|\|_\infty,определяемую максимальным модулем компоненты, а также энергетическую норму \|\|_A, порождаемую симметричной, положительно определенной матрицей A. Так как в конечномерных нормированных пространствах указанные нормы эквивалентны, для доказательства сходимости метода или для его реализации можно пользоваться любой из перечисленных норм.

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

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

Метод простой итерации относится к классу явных стационарных одношаговых итерационных методов решения линейных систем вида с симметричной положительно определенной матрицей A. Метод определяется формулой Ax=b с симметричной положительно определенной матрицей A. Метод определяется формулой

\frac {x^{(s+1)}-x^{(s)}}{\tau}+Ax^{(s)}=b
, где x^{(s)} – текущее приближение, x^{(s+1)} – следующее приближение, \tau \ne 0 – фиксированное число, являющееся параметром метода.

Формула для вычисления нового приближения по предыдущему может быть записана в виде

x^{(s+1)}= -\tau(Ax^{(s)}-b)+x^{(s)}=-\tau \cdot r^{(s)} + x^{(s)}
, где r^{(s)}= Ax^{(s)}-bневязка приближения с номером s.

Далее нетрудно записать явные формулы для отыскания компонент нового вектора x^{(s+1)}:

x_i^{(s+1)}= -\tau\Bigl(\sum_{j=1}^n a_{ij}x_j^{(s)}-b_i\Bigr)+x_i^{(s)}.

Справедлива следующая теорема о сходимости: если матрица A симметрична и положительно определена и параметр метода \tau лежит в интервале (0,\lambda_{max}), где \lambda_{max} - максимальное собственное число матрицы A, метод сходится к точному решению системы (7.49) с любого начального приближения. Известно, что оптимальным значением параметра \tau является значение

\tau^*=\frac {2}{\lambda_{min}+\lambda_{max}}
, где \lambda_{min} и \lambda_{max} – минимальное и максимальное собственные числа матрицы A соответственно.

Напомним, что погрешностью приближения с номером s называют вектор z^{(s)}, определяемый как

z^{(s)}=x^*-x^{(s)}
, где x^* – точное решение системы Ax=b. Для метода простой итерации с оптимальным параметром справедлива следующая оценка:
\|z^{(s+1)}\|_2\le\left(\frac{\mu_A-1}{\mu_A+1}\right)^{s+1}\|z^{(0)}\|_2.

Здесь через \|\|_2 обозначена среднеквадратичная норма вектора, a через \mu_A – число обусловленности матрицы A, согласованное с векторной нормой \|\|_2 и в силу этого вычисляемое как \mu_A=\lambda_{max}/\lambda_{min}.

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

t_1=2n(n+1).

Таким образом, выполнение L итераций метода потребует

T_1=2n(n+1)L
операций.

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

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

Анализ последовательного алгоритма показывает, что основные затраты на s-й итерации – порядка O(n^2) операций – состоят в умножении матрицы A на вектор текущего приближения x^{(s)}. Как результат, при организации параллельных вычислений могут быть использованы известные методы параллельного умножения матрицы на вектор (например, параллельный алгоритм матрично-векторного умножения при ленточном горизонтальном разделении матрицы). Дополнительные вычисления, имеющие порядок сложности O(n), представляют собой операции сложения и умножение на скаляр для векторов. Организация таких вычислений также может быть выполнена в многопоточном режиме.

Несмотря на простоту распараллеливания, данный метод редко применяется, т.к. он обладает существенно меньшей скоростью сходимости по сравнению с иными методами, к рассмотрению которых мы переходим.

Метод верхней релаксации

Рассмотрим систему (7.49) с симметричной, положительно определенной матрицей A размера n\times n. Обозначим через D диагональную матрицу n\times n такую, что ее главная диагональ совпадает с главной диагональю матрицы A. Через L обозначим нижнюю треугольную матрицу n\times n такую, что ее ненулевые (поддиагональные) элементы также совпадают с элементами A, а главная диагональ является нулевой. Аналогично обозначим через R верхнюю треугольную матрицу n\times n, ненулевые (наддиагональные) элемен- ты которой совпадают с элементами A, а главная диагональ также является нулевой. В этом случае для A справедливо представление в виде

A=L+D+R ( 7.5)
Метод верхней релаксации является представителем стационарных одношаговых итерационных методов линейной алгебры и записывается в виде
\frac{(D+\omega L)(x^{(s+1)}-x^{(s)})}{\omega}+Ax^{(s)}=b ( 7.54)
Здесь x^{(s)} – приближение, полученное на итерации с номером s, x^{(s+1)} – следующее приближение, \omega – число (параметр метода), матрицы A, L, D и вектор b определены выше.

Необходимым условием сходимости метода релаксации с любого начального приближения x^{(0)} к точному решению задачи x^* является выполнение условия \omega \in (0,2). Если же матрица A симметрична и положительно определена, то выполнение данного условия является также и достаточным. При этом если \omega \in (01), то говорят о методе нижней релаксации , а при \omega \in (1,2) - о методе верхней релаксации , при \omega =1 метод релаксации будет совпадать с известным методом Зейделя .

Скорость сходимости метода верхней релаксации определяется выбором параметра \omega . Известно [4], что при решении некоторых классов разреженных систем уравнений, метод Зейделя требует O(n^2) итераций, а при надлежащем выборе итерационного параметра \omega метод будет сходиться за O(n) итераций.

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

\omega_{opt}\approx 2-O(h)
. где h – шаг сетки, на которой проводилась дискретизация.

В некоторых случаях можно более точно оценить оптимальный параметр. Известной оценкой [17] является выражение

\omega_{opt}=\frac {2}{1+\sqrt{1-\rho^2(D^{-1}(R+L))}
, где \rho – спектральный радиус матрицы, а R, D, L – из (7.53). Возникающий при этом вопрос об оценке спектрального радиуса может быть решен численно, например, при помощи степенного метода [2].

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

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

Получим формулы для отыскания x^{(s+1)} по предыдущему приближению x^{(s)} в явном виде.

(D+\omega L)(x^{(s+1)}-x^{(s)})+\omega Ax^{(s)}=\omega b, \\
Dx^{(s+1)}+\omega Lx^{(s+1)}-Dx^{(s)}-\omega Lx^{(s)}+\omega Ax^{(s)}=\omega b,\\
Dx^{(s+1)}=-\omega Lx^{(s+1)}+Dx^{(s)}-\omega (A-L)Dx^{(s)}+\omega b.

С учетом того, что A-L=R+D, получаем

Dx^{(s+1)}=-\omega Lx^{(s+1)} +(1-\omega)Dx^{(s)}-\omega Rx^{(s)}+\omega b.

Далее нетрудно записать явные формулы для отыскания компонент нового вектораx^{(s+1)}:

a_{ii}x_i^{(s+1)}=-\omega \sum_{j=1}^{i-1}a_{ii}x_i^{(s+1)}+(1-\omega)a_{ii}x_i^{(s)}-\omega \sum_{j=j+1}^{n}a_{ii}x_i^{(s)} +\omega b_i ( 7.55)

Как следует из формулы (7.55), при подсчете i-й компоненты нового приближения все компоненты, индекс которых меньше i, берутся из нового приближения x^{(s+1)}, а все компоненты, индекс которых больше либо равен i – из старого приближения x^{(s)} . Таким образом, после того, как i-я компонента нового приближения вычислена, i-я компонента старого приближения нигде использоваться не будет. Напротив, для подсчета следующих компонент вектора x^{(s+1)} компоненты с индексом, меньшим или равным i, будут использоваться "в новой версии". В силу этого обстоятельства для реализации метода достаточно хранить только одно (текущее) приближение x^{(s)} , а при расчете следующего приближенияx^{(s+1)} использовать формулу (7.55) для всех компонент по порядку и постепенно обновлять вектор x^{(s)} .

Организация параллельных вычислений

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

Анализ последовательного алгоритма показывает, что основные затраты на s -й итерации – порядка  O(n^2) операций – заключаются в операциях матричного умножения Lx^{(s+1)}и Rx^{(s+1)}. Однако напрямую распараллелить эти операции не представляется возможным, т.к. в методе верхней релаксации не только итерации осуществляются последовательно, но и вычисление компонент вектора очередного приближения x^{(s+1)}также осуществляется последовательно, начиная с первой.

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

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

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

 

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

Скачал архив и незнаю как ничать изучать материал. Видео не воспроизводится (скачено очень много кодеков, различных плееров -- никакого эффекта. Максимум видно часть изображения без звука). При старте 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.