Участник:Bashleon

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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