Участник: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].