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

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

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

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

Дается вещественная симметричная матрица $A = A^T$,

$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}$ $\, \; (1),$

случайный вектор $b$, являющийся первым приближением собственного вектора матрицы, $k$ - количество собственных значений и собственных векторов, которые мы хотим найти, т.е. количество итераций.

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

$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).$

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

После следует вычисление собственных значений $\theta_j$ и собственных векторов $v_j$ полученной трехдиагональной матрицы $T_j$, для чего существует, например, метод "разделяй и властвуй"

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 Вычислительное ядро алгоритма

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

• Умножение матрицы на вектор, $z=Aq_j,$.

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

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

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

ЭТО

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

И ЭТО