Участник:Bashleon

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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