Участник:Смирнова Александра/Нахождение собственных чисел квадратной матрицы методом QR разложения (3)

Нахождение собственных чисел квадратной матрицы методом QR разложения
Последовательный алгоритм
Последовательная сложность	[math]N*O(n^3)[/math]
Объём входных данных	[math] n^2 [/math]
Объём выходных данных	[math] n [/math]
Параллельный алгоритм
Высота ярусно-параллельной формы	-
Ширина ярусно-параллельной формы	-

Основные авторы описания: Смирнова А.С., Киямова А.

---

Задача нахождения собственных значений и собственных векторов для матрицы [math]A[/math] заключается в поиске таких чисел [math]\lambda[/math], которые удовлетворяют уравнению:

[math]Ax=\lambda x[/math], при этом, числа [math]\lambda[/math] называются собственными значениями, а вектора [math]x[/math] - собственными векторами

Данная задача является одной из важнейших задач линейной алгебры. Собственные вектора и собственные значения применяются в различных областях науки: в аналитической геометрии, при решении систем интегральных уравнений, в математической физике. Однако не существует простых алгоритмов прямого вычисления собственных значений для матриц в общем случае, поэтому данная задача на практике решается численными методами. Одним из таких методов является QR-алгоритм.

1 ЧАСТЬ. Свойства и структура алгоритма

1.1 Общее описание алгоритма

---

QR-алгоритм — это численный метод в линейной алгебре, предназначенный для решения полной проблемы собственных значений, то есть отыскания всех собственных чисел и собственных векторов матрицы. Был разработан в конце 1950-х годов независимо В. Н. Кублановской(Россия) и Дж. Фрэнсисом(Англия). Открытию QR-алгоритма предшествовал LR-алгоритм, который использовал LU-разложение вместо QR-разложения. В настоящее время LR-алгоритм используется очень редко ввиду своей меньшей эффективности, однако он был важным шагом на пути к открытию QR-алгоритма.

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

Однако базовый QR-алгоритм может обладать очень низкой скоростью сходимости, поэтому существуют несколько способов ускорить его:

• Перед итерациями привести матрицу [math]A[/math] к подобной ей матрице [math]A_H[/math], которая будет иметь форму Хессенберга. Данный шаг позволит ускорить процесс QR-разложения
• Использовать QR-алгоритм со сдвигами

В данной статье будет рассмотрен только базовый QR-алгоритм.

1.2 Математическое описание алгоритма

---

1.2.1 QR-разложение матрицы

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

[math]A=QR[/math], где [math]Q[/math] - ортогональная матрица ([math]Q^T=Q^{-1}[/math]), [math]R[/math] - верхняя треугольная матрица

Данное разложение называется QR-разложением.

Есть несколько алгоритмов вычисления QR-разложения матрицы. Кратко опишем их в данной статье.

1.2.1.1 Метод Хаусхолдера (отражений) QR-разложения квадратной матрицы

Основная статья: Метод Хаусхолдера (отражений) QR-разложения квадратной матрицы

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

На [math]i[/math]-ом шаге задача [math]i[/math]-ой матрицы отражения заключается в обнулении всех поддиагональных элементов [math]i[/math]-го столбца матрицы [math]A[/math] (при этом столбцы левее [math]i[/math]-го не изменяются). Таким образом, алгоритм состоит из n шагов, на каждом из которых вычисляется очередная матрица отражения, после чего найденное отражение применяется к матрице, которая является результатом предыдущего шага.

Матрица отражения имеет вид [math]E-\frac{1}{\gamma }vv^*[/math], где вектор [math]v[/math] вычисляется из текущего [math]i[/math]-го столбца матрицы при помощи использования операции скалярного произведения векторов. Данное представление матриц отражения позволяет хранить их в виде одного вектора и сводит домножение матрицы отражения на текущую матрицу к арифметическим операциям над векторами (скалярное произведение и сложение векторов).

1.2.1.2 Метод Гивенса (вращений) QR-разложения квадратной матрицы

Основная статья: Метод Гивенса (вращений) QR-разложения квадратной матрицы

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

На каждом шаге задачей матрицы вращения является обнуление одного поддиагонального элемента. Вначале обнуляются все поддиагональные элементы 1-го столбца, затем 2-го и так далее до [math](n-1)[/math]-го. Таким образом, алгоритм состоит из [math]\frac{n(n-1)}{2}[/math] шагов на каждом из которых вычисляется очередная матрица вращения, после чего она применяется к матрице, которая является результатом предыдущего шага.

Матрицы вращения [math]T_{ij}[/math] по ствоей структуре очень близки к единичным матрицам, за исключением четырех элементов: элементы с номерами [math]ii[/math] и [math]jj[/math] содержат некоторое число-параметр [math]c[/math], элементы с номерами [math]ij[/math] и [math]ji[/math] содержат числа-параметры [math]-s[/math] и [math]s[/math] соответственно. Вычисление параметров [math]c[/math] и [math]s[/math] происходит на каждом шаге в зависимости от текущей матрицы и состоит из простых численных арифметических операций. Умножение матрицы вращения на текущую матрицу может быть представлено как последовательное изменение элементов с номерами [math]ik[/math] и [math]kk[/math] для всех столбцов [math]k[/math] находящихся правее столбца [math]i[/math]. Каждое такое изменение по своей структуре эквивалентно операции перемножения двух комплексных чисел.

1.2.2 QR-алгоритм нахождения собственных чисел

Пусть матрица [math]A[/math] - матрица, для которой мы хотим найти собственные значения. Тогда итерационный процесс строится следующим образом:

1. [math]A_0=A[/math]
2. Пусть имеется матрица [math]A_k[/math], тогда матрица [math]A_{k+1}[/math] строится следующим образом:

• Строится QR-разложение: [math]A_k=Q_kR_k[/math]
• Вычисляется [math]A_{k+1}=R_kQ_k[/math]

Заметим, что [math]A_{k+1}=R_kQ_k={Q_{k}^{-1}}Q_kR_kQ_k={Q_{k}^{-1}}A_kQ_k={Q_{k}^{T}}A_kQ_k[/math]

Таким образом матрицы [math]A_{k+1}[/math] и [math]A_{k}[/math] подобны для [math]\forall k[/math], а значит, в силу транзитивности подобия, все матрицы [math]A_{k}[/math] подобны матрице [math]A[/math] и имеют одинаковый набор собственных значений.

1.2.3 Сходимость алгоритма

Предположим, что для [math]\forall m[/math] выполнены следующие условия:

1. [math]A=X\Lambda X^{-1}[/math], где [math]\Lambda =\begin{bmatrix} \Lambda_1 & 0\\ 0 & \Lambda_2 \end{bmatrix}, \Lambda_1\in\mathbb{C}^{m\times m},\Lambda_2\in\mathbb{C}^{r\times r} [/math]

2. [math]\left | \lambda_1 \right | \geq ...\geq \left | \lambda_m \right | \gt \left | \lambda_{m+1} \right | \geq ...\geq \left | \lambda_{m+r} \right | \gt 0 [/math], где [math]\{\lambda_1,...,\lambda_m\} = \lambda(\Lambda_1), \{\lambda_{m+1},...,\lambda_{m+r}\} = \lambda(\Lambda_2) [/math]

3. Ведущая подматрица порядка [math]m[/math] в [math]X^{-1}[/math] невырождена

Тогда при [math] k \rightarrow \infty [/math] последовательность матриц [math]A_k[/math] сходится к матрице с верхней треугольной формой.

Таким образом, на практике необходимо выполнять итерации до тех пор пока матрица [math]A_k[/math] не станет треугольной (также можно продолжать выполнять их пока искомая матрица не будет найдена с некоторой заранее заданной точностью [math]\varepsilon[/math]). Если итерации закончились на шаге [math]N[/math], то числа на диагонали матрицы [math]A_N[/math] будут считаться собственными значениями матрицы [math]A[/math]

1.2.4 Вещественный вариант QR-алгоритма

Если вещественная матрица [math]A[/math] имеет различные вещественные собственные значения, то, как было описано ранее, QR-алгоритм сходится к матрице с верхней треугольной формой, на диагонали которой находятся собственные значения. Однако вещественная матрица может иметь комплексные собственные значения. В данном случае алгоритм будет сходиться не к верхней треугольной матрице, а к блочной верхней треугольной матрице, которая на диагонали содержит блоки 1-го и 2-го порядка. Блоки 1-го порядка содержат различные вещественные собственные значения, блоки 2-го порядка соответствуют парам комплексных сопряженных собственных значений.

1.3 Вычислительное ядро алгоритма

---

QR-алгоритм обладает двумя вычислительными ядрами, которые повторяются на каждой итерации:

1. Вычисление QR-разложения матрицы: [math]A_k=Q_kR_k[/math]. Существует несколько методов решения данной задачи:

• Метод Хаусхолдера (отражений)

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

• Метод Гивенса (вращений)

Вычислительное ядро данного алгоритма состоит из численных арифметических операций, необходимых для вычисления параметров матрицы вращения, и из операций (эквивалентных перемножению комплексных чисел), необходимых для пересчета матрицы на каждом шаге.

2. Перемножение двух плотных матриц: [math]A_{k+1}=R_kQ_k[/math]

1.4 Макроструктура алгоритма

---

Как уже было показано ранее, QR-алгоритм содержит в себе две макрооперации:

1.Вычисление QR-разложения матрицы: [math]A_k=Q_kR_k[/math].

2.Перемножение двух плотных матриц: [math]A_{k+1}=R_kQ_k[/math]

1.5 Схема реализации последовательного алгоритма

---

Опишем необходимый для реализации цикл при помощи псевдокода:

A - исходная матрица
curA - текущая матрица, на основе которой будет вычисляться QR-разложение на каждом шаге
nextA - новая матрица, полученная после перемножения матриц R и Q 
triangular - функция, проверяющая, имеет ли матрица верхнюю треугольную форму.
difference - функция, проверяющая, что элементы двух матриц, стоящие на одинаковых местах, различаются не более чем на eps (данная функция нужна чтобы проверять не только сходимость матрицы к верхней треугольной форме, но и сходимость ее ненулевых элементов).

curA = A
nextA = A
while ( not (triangular(nextA) & difference(curA,nextA,eps)) )
{
    curA = nextA
    findQRdecomposition(curA,Q,R)
    nextA = R*Q
}

1.6 Последовательная сложность алгоритма

---

Подсчитаем сложность одной итерации QR-алгоритма (расчет сложностей для QR-разложения и перемножения матриц представлен в статьях, посвященных данным алгоритмам)

1. QR-разложение матрицы
   • Метод Хаусхолдера (отражений) имеет сложность [math]\frac{4}{3}n^3[/math]
   • Метод Гивенса (вращений) имеет сложность [math]2n^3[/math]
2. Перемножение двух плотных матриц имеет сложность [math]n^3[/math]
3. Проверка матрицы на верхнюю треугольную форму состоит из набора сравнений каждого поддиагонального элемента с нулем (таких элементов [math]\frac{n(n-1)}{2}[/math]) и набора логических операций "И" между результатами сравнений (таких операций будет [math]\frac{n(n-1)}{2}-1[/math])
4. Сравнение новой матрицы с предыдущей состоит из операций вычитания и сравнения для каждой пары ненулевых соответствующих элементов (таких пар будет [math]\frac{n(n+1)}{2}[/math]), а также из набора логических операций "И" между результатами сравнений (таких операций будет [math]\frac{n(n+1)}{2}-1[/math])

Итого, в сумме получаем [math]O(n^3)[/math] - сложность алгоритма на каждой итерации. Если алгоритм остановился на итерации с номером [math]N[/math], то общая сложность алгоритма будет равна [math]N*O(n^3)[/math]

1.7 Информационный граф

---

На рисунке ниже изображен информационный граф QR-алгоритма. Вершины данного графа обозначают следующее:

• QR - операция QR-разложения матрицы
• *_M - операция перемножения матриц
• Triang - операция проверки матрицы на верхнюю треугольную форму
• Dif - операция проверки различия элементов двух матриц не более чем на некоторое заданное число

Информационный граф QR-алгоритма

Подробные графы операций QR-разложения (Метод Хаусхолдера (отражений), Метод Гивенса (вращений)) и перемножения матриц можно найти в статьях, посвященных этим алгоритмам. Далее рассмотрим информационные графы операций Triang и Dif.

---

1.9 Входные и выходные данные алгоритма

---

Входные данные:

квадратная вещественная плотная матрица [math]A[/math]: [math]A \in \R^{n \times n}[/math]

Объем входных данных:

[math]n^2[/math] (необходимо хранить все элементы матрицы)

Выходные данные:

собственные значения матрицы [math]A[/math]

Объем выходных данных:

[math]n[/math] (квадратная матрица размера [math]n \times n[/math] имеет ровно [math]n[/math] собственных значений при этом некоторые из них могут быть комплексными)

---

---

---