Участник:Bashleon

Вариант 4 Башлыков Леонид Александрович 616

Математическое описание алгоритма Ланцоща для вычисления собственных значений симметричной матрицы А требует 1) математического описания метода Ланцоща для построения крыловского подпространства и 2) математического описания т.н. процедуры Рэлея-Ритца.

1 Метод Ланцоща для построения крыловского подпространства

Пусть дано уравнение $A * x = b$, причём вектор $b$ известен и имеется способ вычисления $A * x$.

Вычислим $y_{1} = b$ , $y_{2} = A * b$, $y_{3} = A^2 * b$,...,$y_{n} = A^(n - 1) * b$, где $n$ - порядок матрицы $A$. Определим матрицу $K = [y_{1},...,y_{n}]$.

Тогда $A * K = [A * y_{1},...,A * y_{n - 1}, A * y_{n}] = [y_{2},...,y_{n},A^n * y_{1}]$. В предположении, что матрица $K$ - невырождена, можно вычислить вектор $c = -K^-1 * A^n * y_{1}$. Поскольку первые $n - 1$ столбцов в $A * K$ совпадают с последними $n - 1$ столбцами в $K$, то справедливо: $A * K = K * [e_{2},e_{3},...,e_{n},-c] = K * C$, то есть $K^-1 * A * K = C$, где $C$ - верхняя хессенбергова матрица.

Заменим $K$ ортогональной матрицей $Q$ так, что при любом $k$ линейные оболочки первых $k$ столбцов в $K$ и $Q$ - одно и то же подпространство. Это подпространство называется крыловским и точно определяется как $\kappa_{k}(A, b)$ - линейная оболочка векторов $b$,$A * b$,...,$A^(k - 1) * b$. Суть алгоритма Ланцоша заключается в вычислении стольких первых столбцов в $Q$, сколько необходимо для получения требуемого приближения к решению $A * x = b$ ($A * x = \lambda * x$).

Используем QR-разложение матрицы $K = Q * Q$
Тогда $K^-1 * A * K = (R^-1 * Q^T) * A * (Q * R) = C$, откуда $Q^T * A * Q = R * C * R^-1 = H$

Матрица $H$ является верхней хессенберговой в силу того, что верхней хессенберговой является матрица $C$, а матрицы $R$ и $R^-1$ - верхние треугольные.

Однако алгоритм Ланцоша подразумевает симметричность матрицы $A$. Положим $Q = [Q_{k}, Q_{u}]$, где $Q_{k} = [q_{1},...,q_{k}]$ и $Q_{u} = [q_{k+1},...,q_{n}]$. В итоге получается, что $H = Q^T * A * Q = [Q_{k},Q_{u}]^T * A * [Q_{k},Q_{u}]$ Из симметричности $A$, таким образом, следует симметричность и трехдиагональность матрицы $H = T$. Теперь можно вывести две основные формулы, используемые в методе Ланцоша:

• $A * Q_{j} = \beta_{j - 1} * q_{j - 1} + \alpha_{j} * q_{j} + \beta_{j} * q_{j + 1}$ - получается приравниванием столбцов $j$ в обеих частях равенства $A * Q = Q * T$.
• $q_{j}^T * A * q_{j} = \alpha_{j}$ - получается домножением обеих частей предыдущего соотношения на $q_{j}^T$ и учетом ортогональности $Q$.

2 Процедура Рэлея-Ритца

Пусть $Q = [Q_{k}, Q_{u}]$ - произвольная ортогональная матрица порядка $n$, причём $Q_{k}$ и $Q_{u}$ имеют размеры $n * k$ и $n * (n - k)$.

Пусть $T = Q^T * A * Q = [Q_{k},Q_{u}]^T * A * [Q_{k},Q_{u}]$. Обратим внимание, что в случае $k = 1$ выражение $T_{k} = Q_{k}^T * A * Q_{k}$ - отношение Рэлея $\rho(Q_{1},A)$, то есть число $(Q_{1}^T * A * Q_{1}) / (Q_{1}^T * Q_{1})$. Определим теперь суть процедуры Рэлея-Ритца:

• Процедура Рэлея-Ритца - интерпретация собственных значений матрицы $T = Q^T * A * Q$ как приближений к собственным значениям матрицы $A$. Эти приближения называются числами Ритца. Пусть $T_{k} = V * \lambda * V^T$ - спектральное разложение матрицы $T_{k}$. Столбцы матрицы $Q_{k} * V$ рассматриваются как приближения к соответствующим собственным векторам и называются векторами Ритца.

Существует несколько важных фактов о процедуре Рэлея-Ритца, характеризующих получаемые приближения собственных значений матрицы $A$.

• Минимум величины $||A * Q_{k} - Q_{k} * R||_{2}$ по всем симметричным $k * k$ матрицам $R$ достигается при $R = T_{k}$. При этом $||A * Q_{k} - Q_{k} * R||_{2} = ||T_{ku}||_{2}$. Пусть $T_{k} = V * \lambda * V^T$ - спектральное разложение матрицы $T_{k}$. Минимум величины $||A * P_{k} - P_{k} * D||_{2}$, когда $P_{k}$ пробегает множество $n * k$ матрицы с ортонормированными столбцами, таких, что $span(P_{k}) = span(Q_{k})$, а $D$ пробегает множество диагональных $k * k$ матриц также равен $||T_{ku}||_{2}$ и достигается для $P_{k} = Q_{k} * V$ и $D = \lambda$.

Пусть $T_{k}$, $T_{ku}$, $Q_{k}$ - те же самые матрицы. Пусть $T_{k} = V * \lambda * V^T$ - всё то же спектральное разложение, причём $V = [v_{1},...,v_{k}]$ - ортогональная матрица, а $\lambda = diag(\theta_{1},...,\theta_{k})$. Тогда:

• Найдутся $k$ необязательно наибольших собственных значений $\alpha_{1},...,\alpha_{k}$ матрицы $A$, такие, что $|\theta_{i} - \alpha_{i}| \lt = ||T_{ku}||_{2}$ для $i = 1,...,k$. Если $Q_{k}$ вычислена методом Ланцоша,то правая часть этого выражения равняется $\beta_{k}$, то есть единственному элементу блока $T_{ku}$, который может быть отличен от нуля. Указанный эдемент находится в правом верхнем углу блока.
• $||A * (Q_{k} * v_{i}) - (Q_{k} * v_{i}) * \theta_{i}||_{2} = ||T_{ku} * v_{i}||_{2}$. Таким образом, разность между числом Ритца $\theta_{j}$ и некоторым собственным значением $\alpha$ матрицы $A$ не превышает величины $||T_{ku} * v_{i}||_{2}$, которая может быть много меньше, чем $||T_{ku}||_{2}$. Если матрица $Q_{k}$ вычислена алгоритмом Ланцоша, то эта величина равна $\beta_{k} * |v_{i}(k)|$, где $v_{i}(k)$ - k-ая компонента вектора $v_{i}$. Таким образом, можно дещево вычислить норму невязки - левой части выражения, не выполняя ни одного умножения вектора на матрицы $Q_{k}$ или $A$.
• Не имея информации о спектре матрицы $T_{u}$, нельзя дать какую-либо содержательную оценку для погрешности в векторе Ритца $Q_{k} * v_{i}$. Если известно, что $\theta_{j}$ отделено расстоянием, не меньшим $g$, от прочих собственных значений матриц $T_{k}$, $T_{u}$ и матрица $Q_{k}$ вычислена алгоритмом Ланцоша, то угол $\theta$ между $Q_{k} * v_{i}$ и точным собственным вектором $A$ можно оценить формулой $\frac{1}{2} * \sin{2 * \theta} \lt = \frac{b_{k}}{g}$

В алгоритме Ланцоша из ортонормированных векторов строится матрица $Q_{k} = [q_{1},...,q_{k}]$ и в качестве приближенных собственных значений $A$ принимаются числа Ритца - собственные значения симметричной трёхдиагональной матрицы $T_{k} = Q_{k}^T * A * Q_{k}$.

3 Макроструктура - пункт 1.4

Метод Ланцоша соединяет метод Ланцоша для построения крыловского подпространства в процедурой Рэлея-Ритца, подразумевающей вычисление собственных значений симметричной трёхдиагональной матрицы. Первая часть алгоритма представляет собой строго последовательные итерации, на каждой из которых вначале строится очередной столбец матрицы $Q_{j}$ из первых $j$ ортонормированных векторов Ланцоша и матрица $T_{j} = Q^T_{j} * A * Q_{j}$ порядка $j$. Итерационный процесс завершается, когда $j = k$. Вторая часть алгоритма представляет собой поиск собственных значений и собственных векторов симметричной трёхдиагональной матрицы $T_{j}$, что реализуется отдельным алгоритмом, на практике используется QR-алгоритм.
Особенность методов крыловского подпространства состоит в предположении, что матрица $A$ доступна в виде "черного ящика", а именно подпрограммы, вычисляющей по заданному вектору $z$ произведение $y = A * z$ (а также, возможно, произведение $y = A^T * z$, если матрица $A$ несимметрична). Другими словами, прямой доступ к элементам матрицы и их изменение не используются, т.к.

• Умножение матрицы на вектор - наиболее дешевая нетривиальная операция, которую можно проделать с разреженной матрицей: в случае разреженной $A$, содержащей $m$ ненулевых элементов, для матрично-векторного умножения нужны $m$ скалярных умножений и не более $m$ сложений.
• $A$ может быть не представлена в виде матрицы явно, а доступна именно через подпрограмму для вычисления произведений $A * x$.

Таким образом, как QR-алгоритм поиска собственных значений, собственных векторов матрицы $T_{j}$ и оценки погрешности в них результирующей трехдиагональной матрицы, так и умножение матрицы на вектор внутри каждой итерации могут быть рассмотрены в качестве отдельных макроопераций. Проблемой является то, что именно эти манипуляции порождают основную вычислительную сложность алгоритма. Поскольку матрица, к которой применяется QR-алгоритм, имеет порядок $k$ и на практике бывает мала, сделаем допущение о плотности матрицы $A$. Тогда умножение этой матрицы на вектор внутри каждой итерации нужно будет реализовывать в рамках самого алгоритма. Еще одним важным замечанием является момент вычисления собственных значений матрицы $T_{j}$. Её собственные значения, полученные на итерации $p$, никак не используются на итерации $p + 1$, поэтому целесообразно вынести эту процедуру за основной цикл. Однако на практике в силу малых размеров матрицы $T_{j}$ вычисление соответствующих величин производится непосредственно в теле основного цикла, что позволяет сразу оценить достигнутую точность вычислений.