Алгоритм Ланцоша с выборочной ортогонализацией: различия между версиями
[непроверенная версия] | [непроверенная версия] |
Qmakar (обсуждение | вклад) |
Qmakar (обсуждение | вклад) |
||
Строка 89: | Строка 89: | ||
3) Находим первый вектор матрицы <math> Q: \; \; q_{1} = \frac{b_{j}}{\|b\|_2}</math> | 3) Находим первый вектор матрицы <math> Q: \; \; q_{1} = \frac{b_{j}}{\|b\|_2}</math> | ||
4) Начинается цикл, повторяющийся <math> k </math> раз по переменной <math> j </math> | 4) Начинается цикл, повторяющийся <math> k </math> раз по переменной <math> j </math> | ||
− | 4.1) Считаем произведение матрицы на вектор: <math> | + | 4.1) Считаем произведение матрицы на вектор: <math> z_d = \sum\limits_{y=1}^{n} a_{dy} q_{j_y}, \; d = 1,\,\dots\,, n, </math> где |
<math>z_x - </math> компоненты вектора Ланцоша <math>z. </math> | <math>z_x - </math> компоненты вектора Ланцоша <math>z. </math> | ||
− | 4.2) Скалярно умножаем <math>q_j^T</math> и <math> z: \; \; \alpha_j = \sum\limits_{ | + | 4.2) Скалярно умножаем <math>q_j^T</math> и <math> z: \; \; \alpha_j = \sum\limits_{d=1}^{n}q_{j_d} z_d</math> |
− | 4.3) <math>t</math> присвоить <math>j - 1. </math> В цикле от первого до высчитанного <math> t</math>-го собственного вектора провести | + | 4.3) Вычисляем линейную комбинацию векторов <math>z=z-\alpha_jq_j-\beta_{j-1}q_{j-1}, </math> |
+ | 4.4) <math>t</math> присвоить <math>j - 1. </math> В цикле по <math> i </math> от первого до высчитанного <math> t</math>-го собственного вектора провести | ||
выборочную ортогонализацию к сошедшимся векторам Ритца. То есть: | выборочную ортогонализацию к сошедшимся векторам Ритца. То есть: | ||
− | 4. | + | 4.4.1) Высчитать норму матрицы <math>\|T_{t}\| </math> |
− | 4. | + | 4.4.2) Взять <math>t</math>-ю координату вектора <math> v_i</math> и понять, выполняется ли равенство:<math>\beta_{t}|v_i(t)| \leqslant \sqrt{\varepsilon}\|T_{t}\| \, \, </math> |
+ | Если выполнено, то | ||
+ | 4.4.2.1) Найти вектор Ритца <math>y_{i,t_x} = \sum\limits_{d=1}^{t} q_{d_x} v_{i_d}, \; x = 1,\,\dots\,, n,</math>, где <math>y_{i,t_x}</math> - <math>x</math>-я компонента вектора Ритца <math>y</math>. | ||
+ | 4.4.2.2) Скалярно умножаем <math>y_{i,t}^T</math> и <math> z: \; \; \gamma = \sum\limits_{d=1}^{n}y_{i,t_d} z_d</math> | ||
+ | 4.4.2.3) Вычисляем линейную комбинацию векторов <math>z = z-\gamma\, y_{i,t} </math> | ||
=== Последовательная сложность алгоритма=== | === Последовательная сложность алгоритма=== |
Версия 00:40, 12 декабря 2016
Содержание
- 1 Свойства и структура алгоритма
- 1.1 Общее описание алгоритма
- 1.2 Математическое описание алгоритма
- 1.3 Вычислительное ядро алгоритма
- 1.4 Макроструктура алгоритма
- 1.5 Схема реализации последовательного алгоритма
- 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 Свойства и структура алгоритма
1.1 Общее описание алгоритма
Алгоритм Лáнцоша представляет собой мощный метод для нахождения нескольких собственных значений и собственных векторов симметричной матрицы А и для решения систем линейных уравнений. Алгоритм особенно эффективен, если матрица А разреженная и большого размера. Однако любая практическая реализация этого алгоритма страдает от ошибок округления, т.к. векторы Ланцоша теряют взаимную ортогональность. Для того чтобы поддерживать некоторый уровень ортогональности, появились методы полной переортогонализации и выборочной ортогонализации. В этой работе мы рассмотрим последний метод в качестве способа для поддержания ортогональности среди векторов Ланцоша. Он обладает почти столь же высокой точностью, как алгоритм с полной переортогонализацией, и почти столь же низкой стоимостью, как алгоритм без ортогонализации [1].
Дается вещественная симметричная матрица [math]A = A^T[/math],
- [math] 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] \, \; (1), [/math]
случайный вектор [math]b [/math], являющийся первым приближением собственного вектора матрицы, [math]k [/math] - количество собственных значений и собственных векторов, которые мы хотим найти, т.е. количество итераций.
На каждой итерации строится матрица [math]Q_j = [q_1, q_2, \dots, q_j][/math] размерности [math]n \times j[/math], состоящая из ортонормированных векторов Ланцоша [math]z[/math]. В качестве приближенных собственных значений берутся числа Ритца [math]\theta_i [/math], т.е. собственные значения симметричной трехдиагональной матрицы [math]T_j = Q^T_j A Q_j[/math] размерности [math]j \times j[/math].
- [math] 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]q_j [/math] теряют ортогональность вследствие приобретения больших компонент в направлениях векторов Ритца [math]y_{i,j} = Q_j v_i [/math], отвечающих сошедшимся числам Ритца [math] \theta_i [/math]. Поэтому чтобы построить [math]q_j [/math], предлагается на каждом шаге следить за оценками погрешностей [math]\beta_{t}|v_i(t)|, i = 1 \dots t, t = j - 1 [/math], где [math]v_i(t) [/math] - [math]t[/math]-я компонента собственного вектора [math]v_i [/math]. И когда какая-то оценка становится слишком малой, проводить ортогонализацию вектора Ланцоша [math]z [/math]. Величина [math]\beta_{t}|v_i(t)| [/math] считается малой, если она меньше, чем [math]\sqrt{\varepsilon}||T_{t}|| [/math], где [math]\varepsilon[/math] - доступная машинная точность чисел.
После следует вычисление собственных значений [math] \theta_j [/math] и собственных векторов [math]v_j [/math] полученной трехдиагональной матрицы [math]T_j[/math], для чего существует, например, метод "разделяй и властвуй"[2]
1.2 Математическое описание алгоритма
[math] \beta_0=0,q_0=0[/math] [math] q_{1} = \frac{b_{j}}{\|b\|_2}[/math], где [math] \|b\|_2 = \sqrt{\sum\limits_{j=1}^{n} b_j^2}[/math] [math] for\, j=1\,\, to\, \, k\, \, do:[/math] [math]z=Aq_j, [/math] [math]\alpha_j=q_j^Tz, [/math] [math]z=z-\alpha_jq_j-\beta_{j-1}q_{j-1}, [/math] [math]t = j - 1, [/math] [math]for\, i=1\,\, to\, \, t\, \, do: [/math] [math]if\, \beta_{t}|v_i(t)| \leqslant \sqrt{\varepsilon}\|T_{t}\| \, \, then:[/math] [math]z = z-(y^T_{i,t},z)y_{i,t} [/math], где [math]y_{i,t} = Q_{t}v_i[/math] [math]\beta_{j}=\|z\|_2 [/math] [math]q_{j+1}=z/\beta_{j}, [/math] Строим матрицу [math] T_j[/math] (2), вычисляем собственные значения [math] \theta_j [/math] и собственные векторы [math]v_j [/math] полученной матрицы [math]T_j[/math].
1.3 Вычислительное ядро алгоритма
Выделены следующие вычислительные ядра:
- Умножение матрицы на вектор, [math]z=Aq_j, [/math].
- Ортогонализация по отношению к сошедшимся векторам [math]z = z-(y^T_{i,t},z)y_{i,t} [/math] для [math]i = 1 \dots t[/math]
1.4 Макроструктура алгоритма
Были выделены следующие макрооперации:
1. Умножение матрицы на вектор,
[math]z=Aq_j, [/math].
2. Вычисление вектора [math]q_{j+1} [/math] с помощью линейной комбинации других векторов:
[math]\alpha_j=q_j^Tz, [/math]
[math]z=z-\alpha_jq_j-\beta_{j-1}q_{j-1}, [/math]
[math]q_{j+1}=z/\|z\|_2[/math]
3. Ортогонализация вектора Ланцоша с помощью скалярного произведения:
[math]z = z-(y^T_{i,t},z)y_{i,t} [/math]
4. Для вычисления собственных значений матрицы будет использоваться алгоритм "разделяй и властвуй".
1.5 Схема реализации последовательного алгоритма
[math] \beta_0=0,q_0=0[/math]
[math] q_{1} = \frac{b_{j}}{\|b\|_2}[/math], где [math] \|b\|_2 = \sqrt{\sum\limits_{j=1}^{n} b_j^2}[/math] [math] for\, j=1\,\, to\, \, k\, \, do:[/math] [math]z=Aq_j, [/math] [math]\alpha_j=q_j^Tz, [/math] [math]z=z-\alpha_jq_j-\beta_{j-1}q_{j-1}, [/math] [math]t = j - 1, [/math] [math]for\, i=1\,\, to\, \, t\, \, do: [/math] [math]if\, \beta_{t}|v_i(t)| \leqslant \sqrt{\varepsilon}\|T_{t}\| \, \, then:[/math] [math]z = z-(y^T_{i,t},z)y_{i,t} [/math], где [math]y_{i,t} = Q_{t}v_i[/math] [math]\beta_{j}=\|z\|_2 [/math] [math]q_{j+1}=z/\beta_{j}, [/math] Строим матрицу [math] T_j[/math] (2), вычисляем собственные значения [math] \theta_j [/math] и собственные векторы [math]v_j [/math] полученной матрицы [math]T_j[/math].
1) Инициализируем векторы [math] \beta_0=0,q_0=0[/math] 2) Считаем норму вектора [math] b: \; \; \|b\|_2 = \sqrt{\sum\limits_{j=1}^{n} b_j^2}[/math] 3) Находим первый вектор матрицы [math] Q: \; \; q_{1} = \frac{b_{j}}{\|b\|_2}[/math] 4) Начинается цикл, повторяющийся [math] k [/math] раз по переменной [math] j [/math] 4.1) Считаем произведение матрицы на вектор: [math] z_d = \sum\limits_{y=1}^{n} a_{dy} q_{j_y}, \; d = 1,\,\dots\,, n, [/math] где [math]z_x - [/math] компоненты вектора Ланцоша [math]z. [/math] 4.2) Скалярно умножаем [math]q_j^T[/math] и [math] z: \; \; \alpha_j = \sum\limits_{d=1}^{n}q_{j_d} z_d[/math] 4.3) Вычисляем линейную комбинацию векторов [math]z=z-\alpha_jq_j-\beta_{j-1}q_{j-1}, [/math] 4.4) [math]t[/math] присвоить [math]j - 1. [/math] В цикле по [math] i [/math] от первого до высчитанного [math] t[/math]-го собственного вектора провести выборочную ортогонализацию к сошедшимся векторам Ритца. То есть: 4.4.1) Высчитать норму матрицы [math]\|T_{t}\| [/math] 4.4.2) Взять [math]t[/math]-ю координату вектора [math] v_i[/math] и понять, выполняется ли равенство:[math]\beta_{t}|v_i(t)| \leqslant \sqrt{\varepsilon}\|T_{t}\| \, \, [/math] Если выполнено, то 4.4.2.1) Найти вектор Ритца [math]y_{i,t_x} = \sum\limits_{d=1}^{t} q_{d_x} v_{i_d}, \; x = 1,\,\dots\,, n,[/math], где [math]y_{i,t_x}[/math] - [math]x[/math]-я компонента вектора Ритца [math]y[/math]. 4.4.2.2) Скалярно умножаем [math]y_{i,t}^T[/math] и [math] z: \; \; \gamma = \sum\limits_{d=1}^{n}y_{i,t_d} z_d[/math] 4.4.2.3) Вычисляем линейную комбинацию векторов [math]z = z-\gamma\, y_{i,t} [/math]
1.6 Последовательная сложность алгоритма
1.7 Информационный граф
1.8 Ресурс параллелизма алгоритма
1.9 Входные и выходные данные алгоритма
Входные данные
Матрица [math] A [/math], начальный вектор [math] b [/math] и число собственных значений, которые мы хотим найти [math] k [/math].
Ввиду симметричности из матрицы [math] A [/math] достаточно хранить только ее диагональ и элементы, лежащие выше ее. Таким образом, объем входных данных:
[math]\frac{n (n + 1)}{2} + n + 1 = \frac{n^2+3n+2}{2}.[/math]
Выходные данные
Матрица собственных значений [math] \Lambda = diag(\lambda_1 \dots \lambda_k) [/math], матрица собственных векторов [math] V = [v_1 \dots v_k] [/math] размера [math] n \times k [/math]
Так, объем выходных данных: [math] k + kn [/math]
1.10 Свойства алгоритма
Фишечки, клевые свойства, сложность реализаций.
2 Программная реализация алгоритма
2.1 Особенности реализации последовательного алгоритма
2.2 Локальность данных и вычислений
2.3 Возможные способы и особенности параллельной реализации алгоритма
2.4 Масштабируемость алгоритма и его реализации
ЭТО
2.5 Динамические характеристики и эффективность реализации алгоритма
2.6 Выводы для классов архитектур
2.7 Существующие реализации алгоритма
Алгоритм Ланцоша реализован в различных пакетах, библиотеках и проектах
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