Участник:ArtyomKhakimov/Алгоритм Ланцоша для арифметики с плавающей точкой с выборочной ортогонализацией

Авторы: Хакимов А. С.

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

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

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

На вход алгоритму подаётся [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] \, \;[/math]

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

Матрица [math]Q_j = [q_1, q_2, \dots, q_j][/math] размерности [math]n \times j[/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], для чего существует, например, метод "разделяй и властвуй"

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

[math] q_{1} = b_{j}/\|b\|_2, \beta_0=0, q_0=0[/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]/* Провести выборочную ортогонализацию по отношению,  [/math]
    [math]   к сошедшимся векторам Ритца */,  [/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].