Участник:AleksLevin/Алгоритм Ланцоша вычисления собственных значений симметричной матрицы для точной арифметики (без переортогонализации): различия между версиями

Материал из Алговики
Перейти к навигации Перейти к поиску
Строка 65: Строка 65:
  
 
== Схема реализации последовательного алгоритма ==
 
== Схема реализации последовательного алгоритма ==
 +
С целью полностью не дублировать пункт <math style="vertical-align:0%;>\mathbf{1.2}</math> данной статьи, где представлен очень наглядный и доступный псевдокод алгоритма, постараемся лишь коротко просуммировать и обобщить весь перечень последовательных действий в рассматриваемом алгоритме.
  
 +
В качестве входных данных имеем: вещественная симметричная матрица <math style="vertical-align:0%;>A</math> размера <math style="vertical-align:-10%>n\times n</math> и известный вещественный вектор <math style="vertical-align:0%;>b</math> длины <math style="vertical-align:0%>n</math>.
 +
 +
Далее:
 +
* Инициализируем вещественную константу <math style="vertical-align:-20%;>\beta_0</math> и векторы <math style="vertical-align:-10%;>q_0</math> и <math style="vertical-align:-10%;>q_1</math>
 +
* Итерационно для всех <math style="vertical-align:-20%>j=\overline {1,N}</math>:
 +
** Вычисляем вспомогательный вектор <math style="vertical-align:0%>z</math>
 +
** Вычисляем, используя скалярное произведение векторов, значение <math style="vertical-align:-30%>\alpha_j</math> - элемент на главной диагонали трёхдиагональной матрицы <math style="vertical-align:-30%;>T_j</math>
 +
** Вычисляем значение <math style="vertical-align:-20%>\beta_j</math> - элементы непосредственно над и под главной диагональю той же матрицы <math style="vertical-align:-30%;>T_j</math>
 +
** Вычисляем с помощью нового значения вектора <math style="vertical-align:0%>z</math>  очередной столбец матрицы <math style="vertical-align:-20%>Q</math>, он же вектор Ланцоша : <math style="vertical-align:-20%>\,\,q_{j+1}</math>
 +
** Теперь, когда матрица <math style="vertical-align:-30%;>T_j</math> сформирована на очередном шаге, находим её собственные значения и вектора, например с помощью процедур Рэлея-Ритца (о используемых на практике методах нахождения подробнее было сказано ранее в математическом описании алгоритма)
 +
** Переходим на следующую итерацию: <math style="vertical-align:-20%>j</math>++
  
 
== Последовательная сложность алгоритма ==
 
== Последовательная сложность алгоритма ==

Версия 12:58, 11 ноября 2016

Основные авторы описания: Левин А.Д. (студент, кафедра вычислительных методов, 604 группа)

1 Свойства и структура алгоритмов

1.1 Общее описание алгоритма

Данный алгоритм появился в 1950 г. и носит имя венгерского физика и математика Корнелия Ланцоша (венг. Lánczos Kornél). Алгоритм Ланцоша относится к итерационным методам вычисления собственных значений для матриц столь больших порядков [math] n[/math], что к ним нельзя применить прямые методы из-за ограничений по времени и памяти.

Сам Ланцош указывал, что его метод предназначен для отыскания нескольких собственных векторов симметричных матриц, хотя к методу сразу было обращено внимание, как к способу приведения всей матрицы к трёхдиагональному виду. Двадцатью годами позже канадский математик Крис Пэж показал, что, несмотря на чувствительность к округлениям, алгоритм Ланцоша - эффективное средство вычисления некоторых [math]k[/math] собственных чисел и соответствующих им собственных векторов [1, c.276].

Алгоритм Ланцоша для вычисления собственных значений симметричной матрицы [math]A[/math] соединяет в себе метод Ланцоша для построения последовательности подпространств Крылова и ортонормированных векторов Ланцоша и процедуру Рэлея-Ритца получения оптимальных приближений собственных значений и соответствующих собственных векторов исходной матрицы [math]A[/math] [2, с.381].

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

1.2 Математическое описание алгоритма

Для лучшего понимания описания, данного в этом пункте статьи, рекомендуется ознакомиться с параграфом 6.6 Методы Крыловского подпространства [2, с.313]. Здесь же дано краткое описание всех переменных, математических операций и необходимый теоретический минимум.

Алгоритм Ланцоша для вычисления [math]k[/math] собственных значений и собственных векторов вещественной симметричной матрицы [math]A=A^T[/math] в точной арифметике [2, с.381]:

[math] \begin{align} q_1 = & b/ \|b\|_2,\; \beta_0 = 0,\; q_0 = 0\\ for \; & j = 1 \; to \; k \\ & z = A\,q_j\\ & \alpha_j = q^T_j z\\ & z = z - \alpha_j q_j - \beta_{j-1}q_{j-1}\\ & \beta_j = \|z\|_2\\ & If \; (\beta_j == 0) \; then\; stop\; the\; algorithm \\ & \; q_{j+1} = z / \beta_j \\ & compute\; eigenvalues\; and \;eigenvectors\;of \;matrix \;\;T_j= Q^T_j A Q\;\;and\;estimate \;the\; errors\\ end \; & for \end{align} [/math]


В продемонстрированном выше алгоритме [math]b[/math] - заданный вещественный вектор. Также полагается известным алгоритм вычисления произведения матрицы [math]A[/math] на вектор [math]x[/math]. Введём матрицу Крылова, определяемую следующим соотношением: [math]K_j = [b,Ab,A^2b,...,A^{j-1}b][/math].

Далее, на практике, матрица [math]K[/math] заменяется матрицей [math]Q[/math], такой, что при любом числе [math]k[/math] линейные оболочки первых [math]k[/math] столбцов в [math]K[/math] и [math]Q[/math] являются одним и тем же подпространством [2, c.315]. Тогда матрица [math]Q[/math], в отличие от матрицы [math]K[/math], хорошо обусловлена и легко обратима. В результате получаем матрицу [math]Q_j = [q_1, q_2, \dots, q_j][/math] размерности [math]n\times j[/math], столбцы которой ортогональны, являются базисом подпространства Крылова и называются векторами Ланцоша.

В алгоритме Ланцоша вычислению подлежит столько первых столбцов в матрице [math]Q_j[/math], сколько необходимо для получения требуемого приближения к решению [math]A\,x\,=b\,\,(A\,x=\lambda \, x)[/math].

Получение нулевого [math]\beta_j[/math] в итерационном процессе - событие, желательное в том смысле, что оно свидетельствует о том, что найдено некоторое подпространство, являющееся в точности инвариантным. На практике же нулевое и близкие к нулю значения получаются редко [3, стр. 429].

Затем на каждом шаге цикла формируется симметричную трёхдиагональную матрицу [math]T_j = Q^T_j A Q[/math], к которой применяется процесс Рэлея-Ритца для поиска её собственных значений (на практике для поиска этих собственных значений можно использовать любой из специальных методов, изложенных в [3, [math]\S[/math] 8.4]). Эти собственные значения, они же числа Ритца, и полагаются приближёнными собственными значениями исходной матрицы [math]A[/math].

1.3 Вычислительное ядро алгоритма

Вычислительное ядро рассматриваемого алгоритма (наиболее ресурсно-затратная часть алгоритма) - формирование на каждом шаге цикла промежуточного вектора [math]z[/math], вычисляемого по формуле: [math]z = A\,q_j[/math]

1.4 Макроструктура алгоритма

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

В первой части алгоритма, где итерационно строятся Крыловское подпространство и трёхдиагональная матрица [math]T_j[/math], можно выделить следующие макрооперации:

  • умножение матрицы на вектор (в свою очередь представляет собой умножение вектора на число и сложение векторов);
  • вычисление нормы (скалярное произведение векторов + вычисление квадратного корня);
  • скалярное произведение векторов;
  • сложение векторов и их умножение на вещественные числа

Умножение матрицы на вектор - весьма дешевая и нетривиальная операция, если предполагать, что исходная матрица [math]A[/math] разреженная. Именно она подлежит оптимизации с использованием техник распараллеливания, поскольку она является вычислительным ядром алгоритма, а весь упомянутый процесс построения матрицы [math]T_j[/math] выполняется последовательно.

Вторая часть алгоритма - поиск собственных значений матрицы [math]T_j[/math], в нашем случае с помощью процедуры Рэлея-Ритца, может рассматриваться как отдельная макрооперация.

1.5 Схема реализации последовательного алгоритма

С целью полностью не дублировать пункт [math]\mathbf{1.2}[/math] данной статьи, где представлен очень наглядный и доступный псевдокод алгоритма, постараемся лишь коротко просуммировать и обобщить весь перечень последовательных действий в рассматриваемом алгоритме.

В качестве входных данных имеем: вещественная симметричная матрица [math]A[/math] размера [math]n\times n[/math] и известный вещественный вектор [math]b[/math] длины [math]n[/math].

Далее:

  • Инициализируем вещественную константу [math]\beta_0[/math] и векторы [math]q_0[/math] и [math]q_1[/math]
  • Итерационно для всех [math]j=\overline {1,N}[/math]:
    • Вычисляем вспомогательный вектор [math]z[/math]
    • Вычисляем, используя скалярное произведение векторов, значение [math]\alpha_j[/math] - элемент на главной диагонали трёхдиагональной матрицы [math]T_j[/math]
    • Вычисляем значение [math]\beta_j[/math] - элементы непосредственно над и под главной диагональю той же матрицы [math]T_j[/math]
    • Вычисляем с помощью нового значения вектора [math]z[/math] очередной столбец матрицы [math]Q[/math], он же вектор Ланцоша : [math]\,\,q_{j+1}[/math]
    • Теперь, когда матрица [math]T_j[/math] сформирована на очередном шаге, находим её собственные значения и вектора, например с помощью процедур Рэлея-Ритца (о используемых на практике методах нахождения подробнее было сказано ранее в математическом описании алгоритма)
    • Переходим на следующую итерацию: [math]j[/math]++

1.6 Последовательная сложность алгоритма

1.7 Информационный граф

1.8 Ресурс параллелизма алгоритма

1.9 Входные и выходные данные алгоритма

1.10 Свойства алгоритма

2 Программная реализация алгоритма

2.1 Особенности реализации последовательного алгоритма

2.2 Локальность данных и вычислений

2.3 Возможные способы и особенности параллельной реализации алгоритма

2.4 Масштабируемость алгоритма и его реализации

2.5 Динамические характеристики и эффективность реализации алгоритма

2.6 Выводы для классов архитектур

2.7 Существующие реализации алгоритма

3 Литература

[1] Парлетт Б. - Симметричная проблема собственных значений. Численные методы //М.: Мир. - 1983. - С. 276-294

[2] James W. Demmel - Вычислительная линейная алгебра. Теория и приложения //Мир. - 2001. С. 381-391

[3] Голуб Дж., Ван Лоун Ч. - Матричные вычисления//М.: Мир. - 1999. - С. 426-436