Другие методы

Сдвиги и степенной метод

Напоминаем, что такое степенной метод для вычисления собственых значений.

\[ x_{k+1}:=A x_k, \quad x_{k+1}:=\frac{x_{k+1}}{\left\|x_{k+1}\right\|} . \]

• Он сходится к собственному вектору, соответствующему максимальному по модулю собственному значению.

• Сходимость может быть очень медленной.

• Используем сдвиги: если мы преобразуем матрицу

\[ A:=A-\lambda_k I, \]

то соответствующее собственное значение уменьшится (а нам нужно максимальное по модулю). Это не то, что нам было нужно!

Обратная итерация и итерация Релея

• Для того чтобы из малого собственного значения сделать большое, нам нужно обратить матрицу, и это приводит нас к методу обратной итерации

\[ x_{k+1}=(A-\lambda I)^{-1} x_k, \]

где \(\lambda\) - сдвиг, который близок к собственному значению, которое мы хотим найти. Аналогично степенному методу сходимость линейная.

• Для ускорения сходимости можно использовать итерацию Релея, которая задаётся с помощью адаптивного выбора параметра сдвига:

\[\begin{split} \begin{aligned} x_{k+1} &=\left(A-\lambda_k I\right)^{-1} x_k, \\ \lambda_k &=\frac{\left(A x_k, x_k\right)}{\left(x_k, x_k\right)} \end{aligned} \end{split}\]

В симметричном случае \(A=A^*\) сходимость локально кубическая, и локально квадратичная иначе.

Сингулярные значения и собственные значения

• Сингулярное разложение имеет вид

\[ A=U \Sigma V^* \]

и существует для любой матрицы.

• Его также можно считать способом приведения данной матрицы к диагональному виду с помощью двух унитарных преобразований:

\[ \Sigma=U^* A V . \]

• С помощью двусторонних преобразований Хаусхолдера мы можем привести любую матрицу к бидиагональной форме \(B\).

Неявный QR алгоритм (со сдвигами) вычисляет собственные значения. Но мы не можем применить его напрямую к бидиагональной матрице, поскольку она может быть недиагонализуема в общем случае. Однако задачу вычисления сингулярного разложения можно свести к симметричной задаче на собственные значения двумя способами:

1. Работать с трёхдиагональной матрицей

\[ T=B^* B \]

2. Работать с расширенной матрицей

\[\begin{split} T=\left[\begin{array}{cc} 0 & B \\ B^* & 0 \end{array}\right] \end{split}\]

Случай 1 практически реализуем, если не формировать матрицу \(T\) явно!

Таким образом, задача вычисления сингулярных чисел может быть сведена к задаче вычисления собственных чисел симметричной трёхдиагональной матрицы.

Метод «разделяй и властвуй»

Пусть у нас есть трёхдиагональная матрица и мы разделили её на блоки:

\[\begin{split} T=\left[\begin{array}{cc} T_1^{\prime} & B \\ B^{\top} & T_2^{\prime} \end{array}\right] \end{split}\]

Можем записать матрицу \(T\) в виде

\[\begin{split} T=\left[\begin{array}{cc} T_1 & 0 \\ 0 & T_2 \end{array}\right]+b_m v v^* \end{split}\]

где \(v v^*\) - матрица ранга \(1, v=(0, \ldots, 0,1,1,0, \ldots, 0)^T\).

Пусть мы уже разложили матрицы \(T_1\) и \(T_2\) :

\[ T_1=Q_1 \Lambda_1 Q_1^*, \quad T_2=Q_2 \Lambda_2 Q_2^* \]

Тогда (проверьте!),

\[\begin{split} \left[\begin{array}{cc} Q_1^* & 0 \\ 0 & Q_2^* \end{array}\right] T\left[\begin{array}{cc} Q_1 & 0 \\ 0 & Q_2 \end{array}\right]=D+\rho u u^*, \quad D=\left[\begin{array}{cc} \Lambda_1 & 0 \\ 0 & \Lambda_2 \end{array}\right] \end{split}\]

то есть мы свели задачу к задаче вычисления собственных значений у матрицы вида диагональная матрица плюс матрица малого ранга. Решение задачи изложено в теории в семинарских задачах 10 - 13.

Метод бисекции

• Абсолютно другой подход основан на идеи бисекции

• Дана матрица \(A\), инерция которой определяется как триплет \((\nu, \zeta, \pi)\), где \(\nu\) - число отрицательных, \(\zeta\) - число нулевых, и \(\pi\) - число положительных собственных значений.

• Если \(X\) невырождена, тогда \(\operatorname{Inertia}(A)=\operatorname{Inertia}\left(X^* A X\right)\)

Бисекция с помощью метода Гаусса

• Для данного z мы можем запустить метод Гаусса и получить разложение:

\[ A-z I=L D L^* \]

а инерция для диагональной матрицы вычисляется просто.

• С помощью инерции мы можем легко посчитать число собственных значений в заданном интервале.

• Пример: если \(\operatorname{Inertia}(A)=(5,0,2)\) и после сдвига \(\operatorname{Inertia}(A-z I)=(4,0,3), z \in[a, b]\) тогда это значит, что \(\lambda(A) \in[a, z]\).

Метод Якоби

• Вспомним что такое вращения Гивенса (Якоби): на плоскости им соответствуют ортогональные матрицы вида

\[\begin{split} \left(\begin{array}{cc} \cos \phi & \sin \phi \\ -\sin \phi & \cos \phi \end{array}\right) \end{split}\]

а в \(n\)-мерном пространстве мы выбираем два индекса і и ји вращаем относительно соответствующих элементов \(n\)-мерного вектора.

• Идея метода Якоби состоит в минимизации суммы квадратов недиагональных элементов:

\[ \Gamma(A)=\operatorname{off}\left(U^{*} A U\right), \quad \text { off }^{2}(X)=\sum_{i \neq j}\left|X_{i j}\right|^{2}=\|X\|_{F}^{2}-\sum_{i=1}^{n} x_{i i}^{2} \]

с помощью последовательных вращений Якоби для их зануления.

• Когда элементы выбраны, их легко занулить.

• Главный вопрос: в каком порядке нужно проводить зануление?

• Если мы всегда зануляем максимальный недиагональный элемент, метод имеет глобально линейную сходимость и локально квадратичную.

• На практике используется циклический порядок \(\text { (то есть, }(1,2),(1,3), \ldots,(2,3), \ldots) \)

Метод Якоби был первым численным методом для вычисления собственных значений, предложен в 1846. Преимущества:

• Большая константа в оценке сложности

• Очень точный (высокая относительная точность для малых собственных значений по сравнению с другими методами)

• Хорошая возможность параллелизации