Алгоритм Ланцоша с выборочной ортогонализацией

Материал из Алговики
Перейти к навигации Перейти к поиску

Авторы: Беслаев М. ?, Тевдорадзе С.З. 614. Авторы внесли равнозначный вклад в каждый из разделов данной статьи.

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

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

Алгоритм Лáнцоша представляет собой мощный метод для нахождения нескольких собственных значений и собственных векторов симметричной матрицы А и для решения систем линейных уравнений. Алгоритм особенно эффективен, если матрица А разреженная и большого размера. Однако любая практическая реализация этого алгоритма страдает от ошибок округления, т.к. векторы Ланцоша теряют взаимную ортогональность. Для того чтобы поддерживать некоторый уровень ортогональности, появились методы полной переортогонализации и выборочной ортогонализации. В этой работе мы рассмотрим последний метод в качестве способа для поддержания ортогональности среди векторов Ланцоша. Он обладает почти столь же высокой точностью, как алгоритм с полной переортогонализацией, и почти столь же низкой стоимостью, как алгоритм без ортогонализации [1].


Дается вещественная симметричная матрица [math]\displaystyle{ A = A^T }[/math],

[math]\displaystyle{ A = \begin{pmatrix} a_{11} & a_{12} & a_{13} & \cdots & a_{1\ n-1} & a_{1\ n} \\ a_{12} & a_{22} & a_{23} & \cdots & a_{2\ n-1} & a_{2\ n} \\ a_{13} & a_{23} & a_{33} & \cdots & a_{3\ n-1} & a_{3\ n} \\ \vdots & \vdots & \ddots & \ddots & \ddots & \vdots \\ a_{1\ n-1} & \cdots & \cdots & a_{n-2\ n-1} & a_{n-1\ n-1} & a_{n-1\ n} \\ a_{1\ n} & \cdots & \cdots & a_{n-2\ n} & a_{n-1\ n} & a_{n\ n} \\ \end{pmatrix} }[/math] [math]\displaystyle{ \, \; (1), }[/math]

случайный вектор [math]\displaystyle{ b }[/math], являющийся первым приближением собственного вектора матрицы, [math]\displaystyle{ k }[/math] - количество собственных значений и собственных векторов, которые мы хотим найти, т.е. количество итераций.


На каждой итерации строится матрица [math]\displaystyle{ Q_j = [q_1, q_2, \dots, q_j] }[/math] размерности [math]\displaystyle{ n \times j }[/math], состоящая из ортонормированных векторов Ланцоша [math]\displaystyle{ z }[/math]. В качестве приближенных собственных значений берутся числа Ритца [math]\displaystyle{ \theta_i }[/math], т.е. собственные значения симметричной трехдиагональной матрицы [math]\displaystyle{ T_j = Q^T_j A Q_j }[/math] размерности [math]\displaystyle{ j \times j }[/math].

[math]\displaystyle{ T_j = \begin{pmatrix} \alpha_1 & \beta_1 \\ \beta_1 & \alpha_2 & \beta_2 \\ & \beta_2 & \ddots & \ddots \\ & & \ddots & \ddots & \beta_{j-1} \\ & & & \beta_{j-1} & \alpha_j \end{pmatrix}\; (2). }[/math]

Однако, векторы [math]\displaystyle{ q_j }[/math] теряют ортогональность вследствие приобретения больших компонент в направлениях векторов Ритца [math]\displaystyle{ y_{i,j} = Q_j v_i }[/math], отвечающих сошедшимся числам Ритца [math]\displaystyle{ \theta_i }[/math]. Поэтому чтобы построить [math]\displaystyle{ q_j }[/math], предлагается на каждом шаге следить за оценками погрешностей [math]\displaystyle{ \beta_{t}|v_i(t)|, i = 1 \dots t, t = j - 1 }[/math], где [math]\displaystyle{ v_i(t) }[/math] - [math]\displaystyle{ t }[/math]-я компонента собственного вектора [math]\displaystyle{ v_i }[/math]. И когда какая-то оценка становится слишком малой, проводить ортогонализацию вектора Ланцоша [math]\displaystyle{ z }[/math]. Величина [math]\displaystyle{ \beta_{t}|v_i(t)| }[/math] считается малой, если она меньше, чем [math]\displaystyle{ \sqrt{\varepsilon}||T_{t}|| }[/math], где [math]\displaystyle{ \varepsilon }[/math] - доступная машинная точность чисел.


После следует вычисление собственных значений [math]\displaystyle{ \theta_j }[/math] и собственных векторов [math]\displaystyle{ v_j }[/math] полученной трехдиагональной матрицы [math]\displaystyle{ T_j }[/math], для чего существует, например, метод "разделяй и властвуй"[2]

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

[math]\displaystyle{  \beta_0=0,q_0=0 }[/math]
[math]\displaystyle{  q_{1} = \frac{b_{j}}{\|b\|_2} }[/math], где [math]\displaystyle{  \|b\|_2 = \sqrt{\sum\limits_{j=1}^{n} b_j^2} }[/math]
[math]\displaystyle{  for\, j=1\,\, to\, \, k\, \, do: }[/math]
    [math]\displaystyle{ z=Aq_j,   }[/math]
    [math]\displaystyle{ \alpha_j=q_j^Tz,  }[/math]
    [math]\displaystyle{ z=z-\alpha_jq_j-\beta_{j-1}q_{j-1},   }[/math]
    [math]\displaystyle{ t = j - 1,   }[/math]
    [math]\displaystyle{ for\, i=1\,\, to\, \, t\, \, do:  }[/math]
        [math]\displaystyle{ if\,  \beta_{t}|v_i(t)| \leqslant \sqrt{\varepsilon}\|T_{t}\| \, \,  then: }[/math]
            [math]\displaystyle{ z = z-(y^T_{i,t},z)y_{i,t}  }[/math], где [math]\displaystyle{ y_{i,t} = Q_{t}v_i }[/math] 
    [math]\displaystyle{ \beta_{j}=\|z\|_2  }[/math]
    [math]\displaystyle{ q_{j+1}=z/\beta_{j},  }[/math]
    Строим матрицу  [math]\displaystyle{  T_j }[/math] (2), вычисляем собственные значения  [math]\displaystyle{  \theta_j  }[/math] и собственные векторы [math]\displaystyle{ v_j  }[/math] полученной матрицы [math]\displaystyle{ T_j }[/math].

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

Выделены следующие вычислительные ядра:

  • Умножение матрицы на вектор, [math]\displaystyle{ z=Aq_j, }[/math].


  • Ортогонализация по отношению к сошедшимся векторам [math]\displaystyle{ z = z-(y^T_{i,t},z)y_{i,t} }[/math] для [math]\displaystyle{ i = 1 \dots t }[/math]


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

Были выделены следующие макрооперации:

1. Умножение матрицы на вектор,

[math]\displaystyle{ z=Aq_j, }[/math].

2. Вычисление вектора [math]\displaystyle{ q_{j+1} }[/math] с помощью линейной комбинации других векторов:

[math]\displaystyle{ \alpha_j=q_j^Tz, }[/math]

[math]\displaystyle{ z=z-\alpha_jq_j-\beta_{j-1}q_{j-1}, }[/math]

[math]\displaystyle{ q_{j+1}=z/\|z\|_2 }[/math]

3. Ортогонализация вектора Ланцоша с помощью скалярного произведения:

[math]\displaystyle{ z = z-(y^T_{i,t},z)y_{i,t} }[/math]

4. Для вычисления собственных значений матрицы будет использоваться алгоритм "разделяй и властвуй".

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

[math]\displaystyle{ 1) }[/math] Инициализируются векторы [math]\displaystyle{  \beta_0=0,q_0=0, }[/math]
[math]\displaystyle{ 2) }[/math] Считается норма вектора [math]\displaystyle{  b: \; \;  \|b\|_2 = \sqrt{\sum\limits_{j=1}^{n} b_j^2}, }[/math]
[math]\displaystyle{ 3) }[/math] Находится первый вектор матрицы [math]\displaystyle{  Q: \; \; q_{1} = \frac{b_{j}}{\|b\|_2}, }[/math]
[math]\displaystyle{ 4) }[/math] Начинается цикл, повторяющийся [math]\displaystyle{  k  }[/math] раз по переменной [math]\displaystyle{  j:  }[/math]
   [math]\displaystyle{ 4.1) }[/math] Считается произведение матрицы на вектор: [math]\displaystyle{  z_d = \sum\limits_{y=1}^{n} a_{dy} q_{j_y}, \; d = 1,\,\dots\,, n,  }[/math] где 
   [math]\displaystyle{ z_d -  }[/math] компоненты вектора Ланцоша [math]\displaystyle{ z,  }[/math]
   [math]\displaystyle{ 4.2) }[/math] Скалярно умножается [math]\displaystyle{ q_j^T }[/math]  и [math]\displaystyle{  z: \; \; \alpha_j = \sum\limits_{d=1}^{n}q_{j_d} z_d, }[/math]
   [math]\displaystyle{ 4.3) }[/math] Вычисляется линейная комбинация векторов [math]\displaystyle{ z=z-\alpha_jq_j-\beta_{j-1}q_{j-1},   }[/math]
   [math]\displaystyle{ 4.4) }[/math] [math]\displaystyle{ t }[/math] Присваивается [math]\displaystyle{ j - 1.   }[/math] В цикле по [math]\displaystyle{  i  }[/math] от первого до посчитанного [math]\displaystyle{  t }[/math]-го собственного вектора проводится
   выборочная ортогонализация к сошедшимся векторам Ритца. То есть:
      [math]\displaystyle{ 4.4.1) }[/math] Считается норма матрицы [math]\displaystyle{ \|T_{t}\|_2,  }[/math] [math]\displaystyle{  \| T \|_2 =  \sup\limits_{\| x \|_2 = 1} \| T x \|_2 = \sup\limits_{(x, x) = 1} \sqrt{(Tx, Tx)} }[/math], подчиненная векторной норме [math]\displaystyle{  \| x \|_2 = \sqrt{\sum_{i = 1}^n |x_i|^2}  }[/math]. Т.е. 
      [math]\displaystyle{  \|    T_{t} \|_2 =  \max\limits_{i=1 \dots t} \theta_i }[/math], где [math]\displaystyle{  \theta_i }[/math] - собственные значения матрицы [math]\displaystyle{  T, }[/math]
      [math]\displaystyle{ 4.4.2) }[/math] Берется [math]\displaystyle{ t }[/math]-ю координата вектора [math]\displaystyle{  v_i }[/math]  и проверяется, выполняется ли равенство:[math]\displaystyle{ \beta_{t}|v_i(t)| \leqslant \sqrt{\varepsilon}\|T_{t}\| \, \, , }[/math]
      Если выполняется, то
          [math]\displaystyle{ 4.4.2.1) }[/math] Ищется вектор Ритца [math]\displaystyle{ y_{i,t_d} = \sum\limits_{x=1}^{t} q_{x_d} v_{i_x}, \; d = 1,\,\dots\,, n, }[/math], где  [math]\displaystyle{ y_{i,t_d} }[/math] -  [math]\displaystyle{ d }[/math]-я компонента вектора Ритца [math]\displaystyle{ y, }[/math]
          [math]\displaystyle{ 4.4.2.2) }[/math] Скалярно умножается [math]\displaystyle{ y_{i,t}^T }[/math]  и [math]\displaystyle{  z: \; \; \gamma = \sum\limits_{d=1}^{n}y_{i,t_d} z_d, }[/math]
          [math]\displaystyle{ 4.4.2.3) }[/math] Вычисляется линейная комбинация векторов [math]\displaystyle{ z = z-\gamma\,  y_{i,t},  }[/math]
          [math]\displaystyle{ 4.4.2.4) }[/math] Увеличивается [math]\displaystyle{ i  }[/math]. Производится возврат к пункту [math]\displaystyle{ 4.4.2, }[/math]
   [math]\displaystyle{ 4.5) }[/math] Ищется норма вектора [math]\displaystyle{  z: \; \; \beta_{j}=\|z\|_2,  }[/math]
   [math]\displaystyle{ 4.6) }[/math] Ищется вектор Ланцоша [math]\displaystyle{ q_{j+1}=z/\beta_{j},  }[/math]
   [math]\displaystyle{ 4.7) }[/math] Вычисляется новое собственное значение [math]\displaystyle{  \theta_j  }[/math] и собственный вектор [math]\displaystyle{  v_j  }[/math] для полученной матрицы [math]\displaystyle{  T_j, }[/math]
   [math]\displaystyle{ 4.8) }[/math] Увеличивается [math]\displaystyle{  j  }[/math], производится возврат к шагу [math]\displaystyle{ 4.1. }[/math]

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

Рассмотрим последовательную сложность алгоритма. 

[math]\displaystyle{ 1) }[/math] Инициализируются векторы,
[math]\displaystyle{ 2) }[/math] Считается норма вектора: [math]\displaystyle{  n  }[/math] сложений и умножений, одно вычисление корня,
[math]\displaystyle{ 3) }[/math] Находится первый вектор матрицы: [math]\displaystyle{  n  }[/math] делений,
[math]\displaystyle{ 4) }[/math] Начинается цикл, повторяющийся [math]\displaystyle{  k  }[/math] раз по переменной : [math]\displaystyle{    j  }[/math]:
   [math]\displaystyle{ 4.1) }[/math] Считается произведение матрицы на вектор: [math]\displaystyle{    n  }[/math] сложений и умножений,
   [math]\displaystyle{ 4.2) }[/math] Скалярно умножается два вектора : [math]\displaystyle{    n  }[/math] сложений и умножений,
   [math]\displaystyle{ 4.3) }[/math] Вычисляется линейная комбинация векторов : [math]\displaystyle{    2n  }[/math] вычитаний и умножений,
   [math]\displaystyle{ 4.4) }[/math] Повторяется цикл : [math]\displaystyle{    j-1  }[/math] раз. На каждом из этапов:
      [math]\displaystyle{ 4.4.1) }[/math] Считается норма матрицы: : [math]\displaystyle{    j-1  }[/math] сравнений,
      [math]\displaystyle{ 4.4.2) }[/math] Проверяется равенство::  2 умножения, одно сравнение
          [math]\displaystyle{ 4.4.2.1) }[/math] Ищется вектор Ритца : [math]\displaystyle{    n*(j-1)  }[/math] сложений и умножений,
          [math]\displaystyle{ 4.4.2.2) }[/math] Скалярно умножаются два вектора:  [math]\displaystyle{    n  }[/math] сложений и умножений,
          [math]\displaystyle{ 4.4.2.3) }[/math] Вычисляется линейная комбинация векторов : [math]\displaystyle{    n  }[/math] умножений и вычитаний,
          [math]\displaystyle{ 4.4.2.4) }[/math] Увеличивается [math]\displaystyle{ i  }[/math]. Одно сложение.
     Таким образом, т.к. внутренний цикл повторяется [math]\displaystyle{    j-1  }[/math] раз, а [math]\displaystyle{    j  }[/math] в свою очередь пробегает от [math]\displaystyle{ 1 \dots k  }[/math], значит, полное количество операций в пункте [math]\displaystyle{ 4.4) }[/math], которое будет занимать весь внутренний цикл будет: сравнений [math]\displaystyle{  \frac{(k-1)k(k+1)}{3} + \frac{k(k-1)}{2}  }[/math], умножений: [math]\displaystyle{  2\frac{k(k-1)}{2}  + \frac{(k-1)k(k+1)}{3} + \frac{(2n)k(k-1)}{2}  }[/math], сложений [math]\displaystyle{ \frac{(k-1)k(k+1)}{3} + \frac{(n+1)k(k-1)}{2} }[/math] и вычитаний [math]\displaystyle{  \frac{(n)k(k-1)}{2} }[/math]
   [math]\displaystyle{ 4.5) }[/math] Ищется норма вектора: [math]\displaystyle{  n  }[/math] сложений и умножений
   [math]\displaystyle{ 4.6) }[/math] Ищется вектор Ланцоша [math]\displaystyle{  n  }[/math] делений
   [math]\displaystyle{ 4.7) }[/math] Увеличивается [math]\displaystyle{  j  }[/math], одно сложение. 
   
   Так, пункты  [math]\displaystyle{ 4.1 - 4.3, 4.5 - 4.7 }[/math] занимают [math]\displaystyle{  k*(n+n+n+1) }[/math] сложений, [math]\displaystyle{ k(n+n+2n+n) }[/math] умножений, [math]\displaystyle{ 2nk }[/math] вычитаний и [math]\displaystyle{ n }[/math] делений.



Подводя итог, общее количество операций умножения и деления [math]\displaystyle{ n + n + 2\frac{k(k-1)}{2}  + \frac{(k-1)k(k+1)}{3} + \frac{(2n)k(k-1)}{2} +  k(n+n+2n+n) + n }[/math], т.е. [math]\displaystyle{ O( k^3+nk^2+k^2+kn+n ) }[/math]
общее количество сложений и вычитаний: [math]\displaystyle{  n + \frac{(k-1)k(k+1)}{3} + \frac{(n+1)k(k-1)}{2}+ \frac{(n)k(k-1)}{2} + k*(n+n+n+1) + 2nk }[/math], т.е.  [math]\displaystyle{ O( k^3+nk^2+kn+n ) }[/math]
И [math]\displaystyle{  \frac{(k-1)k(k+1)}{3} + \frac{k(k-1)}{2}  }[/math] сравнений ( [math]\displaystyle{ O( k^3+k^2) }[/math]).


рис. 1: Информационный граф

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

На рисунке (рис. 1) изображен информационный граф. Ниже приведено его описание:

  • Вершины 1 соответствуют операции [math]\displaystyle{ Aq_j, }[/math].
  • Вершины 2 соответствуют операции скалярного произведения [math]\displaystyle{ q_j^Tz }[/math].
  • Вершины 3 соответствуют операции [math]\displaystyle{ z-\sum\nolimits_{i=1}^{j-1}(z^Tq_i)q_i }[/math].
  • Вершины 4 соответствуют операции [math]\displaystyle{ z-\alpha_jq_j-\beta_{j-1}q_{j-1} }[/math].
  • Вершины 5 соответствуют операции [math]\displaystyle{ ||z|| }[/math].
  • Вершины 6 соответствуют операции [math]\displaystyle{ z/\beta_j }[/math].

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

Хотя алгоритм является итерационным, возможно распараллелить внутри каждой итерации умножение матрицы на вектор и процесс переортогонализации Грамма-Шмидта.

Процесс умножения матриц matvec можно распараллелить несколькими способами[1].

Как видно из графа, для выполнения j-й операции, необходимо выполнить следующие ярусы:

  1. [math]\displaystyle{ n }[/math] ярусов сложений с [math]\displaystyle{ n }[/math] операциями умножения в каждом (вычисление [math]\displaystyle{ z }[/math]);
  2. [math]\displaystyle{ n }[/math] ярусов сложений с [math]\displaystyle{ 1 }[/math] операцией умножения в каждом (вычисление [math]\displaystyle{ \alpha_j }[/math]);
  3. [math]\displaystyle{ j }[/math] ярусов сложений с [math]\displaystyle{ n }[/math] операциями умножения в каждом, [math]\displaystyle{ n }[/math] ярусов сложений с [math]\displaystyle{ j }[/math] операциями умножения в каждом, [math]\displaystyle{ 1 }[/math] ярус вычитаний размера [math]\displaystyle{ n }[/math] (первая переортогонализация);
  4. [math]\displaystyle{ j }[/math] ярусов сложений с [math]\displaystyle{ n }[/math] операциями умножения в каждом, [math]\displaystyle{ n }[/math] ярусов сложений с [math]\displaystyle{ j }[/math] операциями умножения в каждом, [math]\displaystyle{ 1 }[/math] ярус вычитаний размера [math]\displaystyle{ n }[/math] (вторая переортогонализация);
  5. [math]\displaystyle{ n }[/math] ярусов сложений с [math]\displaystyle{ 1 }[/math] операцией умножения в каждом (вычисление [math]\displaystyle{ \beta_j }[/math]);
  6. [math]\displaystyle{ 1 }[/math] ярус делений размера [math]\displaystyle{ n }[/math] (вычисление [math]\displaystyle{ q_{j + 1} }[/math]).

Если учесть весь ресурс параллелизма, можно получить мультипликативную вычислительную сложность всего алгоритма: [math]\displaystyle{ O(n^2k+nk^2)/p) }[/math].

Вычислительная мощность данного алгоритма, как отношение числа операций к суммарному объему входных и выходных, данных равняется [math]\displaystyle{ O(k) }[/math].

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

Входные данные

Матрица [math]\displaystyle{ A }[/math], начальный вектор [math]\displaystyle{ b }[/math] и число собственных значений, которые мы хотим найти [math]\displaystyle{ k }[/math].

Ввиду симметричности из матрицы [math]\displaystyle{ A }[/math] достаточно хранить только ее диагональ и элементы, лежащие выше ее. Таким образом, объем входных данных:

[math]\displaystyle{ \frac{n (n + 1)}{2} + n + 1 = \frac{n^2+3n+2}{2}. }[/math]

Выходные данные

Матрица собственных значений [math]\displaystyle{ \Lambda = diag(\theta_1 \dots \theta_k) }[/math], матрица собственных векторов [math]\displaystyle{ V = [v_1 \dots v_k] }[/math] размера [math]\displaystyle{ n \times k }[/math]

Так, объем выходных данных: [math]\displaystyle{ k + kn }[/math]

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

Если [math]\displaystyle{ A }[/math] эрмитова матрица, то алгоритм Ланцоша и Bi-Lanczos сходятся к одинаковым трехдиагональным матрицам Ритца[2].

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

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

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

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

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

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

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

ЭТО

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

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

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

Алгоритм Ланцоша реализован в различных пакетах, библиотеках и проектах

NAG Library http://www.nag.com/content/nag-library C, C++, Fortran, C#, MATLAB, R

ARPACK https://people.sc.fsu.edu/~jburkardt/m_src/arpack/arpack.html MATLAB

GrapLab https://turi.com/products/create/open_source.html C++

Gaussian Belief Propagation Matlab Package http://www.cs.cmu.edu/~bickson/gabp/ MATLAB

The IETL Project http://www.comp-phys.org/software/ietl/ C++

LANSO/PLANSO http://web.cs.ucdavis.edu/~bai/ET/lanczos_methods/overview_PLANSO.html Fortran

Julia Math https://github.com/JuliaMath/IterativeSolvers.jl Julia

Литература

  1. ↑ "The Lanczos Algorithm With Partial Reorthogonalization", Horst D. Simon, Mathematics of computation, volume 42, number 165, pages 115-142 (1984)
  2. ↑ Дж. Деммель «Вычислительная линейная алгебра», c. 232, алгоритм 5.2
Источник — https://xn--80aecvaldz.xn--80ahdrnmvk.xn--p1ai/w/ru/index.php?title=Алгоритм_Ланцоша_с_выборочной_ортогонализацией&oldid=21413