Уровень алгоритма

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

Материал из Алговики
Перейти к навигации Перейти к поиску
 
(не показано 85 промежуточных версий 5 участников)
Строка 1: Строка 1:
−
{{Assignment}}
+
{{Assignment|VadimVV|Frolov}}
−
 
 
 
{{algorithm
 
{{algorithm
 
| name              = Нахождение собственных чисел квадратной матрицы методом QR разложения
 
| name              = Нахождение собственных чисел квадратной матрицы методом QR разложения
−
| serial_complexity = <math>N*O(n^3)</math>
+
| serial_complexity = <math>O(n^3)+N*O(n^2)</math>
−
| pf_height        = <math>N*O(n)</math>
+
| pf_height        = <math>O(n^2)+N*O(n)</math>
 
| pf_width          = <math>O(n^2)</math>
 
| pf_width          = <math>O(n^2)</math>
 
| input_data        = <math> n^2 </math>
 
| input_data        = <math> n^2 </math>
Строка 10: Строка 9:
 
}}
 
}}
  
−
Основные авторы описания: [[Участник:Смирнова_Александра|Смирнова А.С.]], [[Участник:Александра_Киямова|Киямова А.]]
+
Основные авторы описания:  
 +
*[[Участник:Смирнова_Александра|Смирнова А.С.]] - описание теоретической части
 +
*[[Участник:Александра_Киямова|Киямова А.]] - описание программистской части
 +
Содержательно все пункты обсуждались авторами вместе, и их вклад был равным.
 +
 
  
 
----
 
----
Строка 16: Строка 19:
 
Задача нахождения собственных значений и собственных векторов для матрицы <math>A</math> заключается в поиске таких чисел <math>\lambda</math>, которые удовлетворяют уравнению:
 
Задача нахождения собственных значений и собственных векторов для матрицы <math>A</math> заключается в поиске таких чисел <math>\lambda</math>, которые удовлетворяют уравнению:
  
−
<math style="vertical-align:0%">Ax=\lambda x</math>, при этом, числа <math>\lambda</math> называются собственными значениями, а вектора <math>x</math> - собственными векторами <ref>Ильин В.А., Ким Г.Д. "Линейная алгебра и аналитическая геометрия"</ref>
+
<math style="vertical-align:0%">Ax=\lambda x</math>, при этом, числа <math>\lambda</math> называются собственными значениями, а вектора <math>x</math> - собственными векторами<ref>Ильин В.А., Ким Г.Д. "Линейная алгебра и аналитическая геометрия".</ref>.
  
 
Данная задача является одной из важнейших задач линейной алгебры. Собственные вектора и собственные значения применяются в различных областях науки: в аналитической геометрии, при решении систем интегральных уравнений, в математической физике. Однако не существует простых алгоритмов прямого вычисления собственных значений для матриц в общем случае, поэтому данная задача на практике решается численными методами. Одним из таких методов является QR-алгоритм.
 
Данная задача является одной из важнейших задач линейной алгебры. Собственные вектора и собственные значения применяются в различных областях науки: в аналитической геометрии, при решении систем интегральных уравнений, в математической физике. Однако не существует простых алгоритмов прямого вычисления собственных значений для матриц в общем случае, поэтому данная задача на практике решается численными методами. Одним из таких методов является QR-алгоритм.
  
−
== ЧАСТЬ. Свойства и структура алгоритма ==
+
== Свойства и структура алгоритма ==
  
 
=== Общее описание алгоритма ===
 
=== Общее описание алгоритма ===
Строка 30: Строка 33:
 
Суть базового 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_N</math> при помощи QR-разложения. Матрица <math>A_N</math> является правой верхней треугольной матрицей, а значит ее диагональ содержит собственные значения. В силу подобия матриц <math>A</math> и <math>A_N</math> их наборы собственных значений совпадают. Таким образом задача поиска собственных значений матрицы <math>A</math> сводится к задаче выведения матрицы <math>A_N</math> и поиска собственных значений для нее, что является тривиальной задачей.
  
−
Однако базовый QR-алгоритм может обладать очень низкой скоростью сходимости, поэтому существуют несколько способов ускорить его:
+
Однако базовый QR-алгоритм может обладать очень низкой скоростью сходимости, поэтому существует несколько способов ускорить его:
−
* Перед итерациями привести матрицу <math>A</math> к подобной ей матрице <math>A_H</math>, которая будет иметь форму Хессенберга. Данный шаг позволит ускорить процесс QR-разложения
+
* Перед итерациями привести матрицу <math>A</math> к подобной ей матрице <math>A_H</math>, которая будет иметь форму Хессенберга. Данный шаг позволит ускорить процесс QR-разложения.
−
* Использовать QR-алгоритм со сдвигами
+
* Использовать QR-алгоритм со сдвигами. Это позволит уменьшить количество итераций алгоритма.
−
 
+
В дальнейшем, в данной статье под модифицированным алгоритмом будет пониматься алгоритм, использующий сдвиги и матрицу с формой Хессенберга. Под базовым алгоритмом будет пониматься алгоритм, не использующий данные приемы.
−
В данной статье будет рассмотрен только базовый QR-алгоритм.
 
  
 
=== Математическое описание алгоритма ===
 
=== Математическое описание алгоритма ===
Строка 44: Строка 46:
 
Основой алгоритма является тот факт, что любую вещественную матрицу можно разложить на произведение двух матриц следующего вида:
 
Основой алгоритма является тот факт, что любую вещественную матрицу можно разложить на произведение двух матриц следующего вида:
  
−
::<math>A=QR</math>, где <math>Q</math> - ортогональная матрица (<math>Q^T=Q^{-1}</math>), <math>R</math> - верхняя треугольная матрица
+
::<math>A=QR</math>, где <math>Q</math> - ортогональная матрица (<math>Q^T=Q^{-1}</math>), <math>R</math> - верхняя треугольная матрица.
  
 
Данное разложение называется QR-разложением.
 
Данное разложение называется QR-разложением.
Строка 58: Строка 60:
 
Суть метода Хаусхолдера заключается в последовательном приведении матрицы <math style="vertical-align:0%">A</math> к верхней треугольной форме при помощи домножения ее на так называемые матрицы отражения. Получившаяся треугольная матрица будет искомой матрицей <math style="vertical-align:0%">R</math>, а матрица <math>Q</math> будет равна произведению сопряженных матриц отражения.
 
Суть метода Хаусхолдера заключается в последовательном приведении матрицы <math style="vertical-align:0%">A</math> к верхней треугольной форме при помощи домножения ее на так называемые матрицы отражения. Получившаяся треугольная матрица будет искомой матрицей <math style="vertical-align:0%">R</math>, а матрица <math>Q</math> будет равна произведению сопряженных матриц отражения.
  
−
На <math style="vertical-align:0%">i</math>-ом шаге задача <math style="vertical-align:0%">i</math>-ой матрицы отражения заключается в обнулении всех поддиагональных элементов <math style="vertical-align:0%">i</math>-го столбца матрицы <math style="vertical-align:0%">A</math> (при этом столбцы левее <math style="vertical-align:0%">i</math>-го не изменяются). Таким образом, алгоритм состоит из <math style="vertical-align:0%">n</math> шагов, на каждом из которых вычисляется очередная матрица отражения, после чего найденное отражение применяется к матрице, которая является результатом предыдущего шага.
+
На <math style="vertical-align:0%">i</math>-ом шаге задача <math style="vertical-align:0%">i</math>-ой матрицы отражения заключается в обнулении всех поддиагональных элементов <math style="vertical-align:0%">i</math>-го столбца матрицы <math style="vertical-align:0%">A</math> (при этом столбцы левее <math style="vertical-align:0%">i</math>-го не изменяются). Таким образом, алгоритм состоит из <math style="vertical-align:0%">n-1</math> шагов, на каждом из которых вычисляется очередная матрица отражения, после чего найденное отражение применяется к матрице, которая является результатом предыдущего шага.
  
 
Матрица отражения имеет вид <math>E-\frac{1}{\gamma }vv^*</math>, где вектор <math>v</math> вычисляется из текущего <math style="vertical-align:0%">i</math>-го столбца матрицы при помощи использования операции скалярного произведения векторов. Данное представление матриц отражения позволяет хранить их в виде одного вектора и сводит домножение матрицы отражения на текущую матрицу к арифметическим операциям над векторами (скалярное произведение и сложение векторов).
 
Матрица отражения имеет вид <math>E-\frac{1}{\gamma }vv^*</math>, где вектор <math>v</math> вычисляется из текущего <math style="vertical-align:0%">i</math>-го столбца матрицы при помощи использования операции скалярного произведения векторов. Данное представление матриц отражения позволяет хранить их в виде одного вектора и сводит домножение матрицы отражения на текущую матрицу к арифметическим операциям над векторами (скалярное произведение и сложение векторов).
Строка 75: Строка 77:
  
 
Пусть матрица <math>A</math> - матрица, для которой мы хотим найти собственные значения. Тогда итерационный процесс строится следующим образом:
 
Пусть матрица <math>A</math> - матрица, для которой мы хотим найти собственные значения. Тогда итерационный процесс строится следующим образом:
−
#<math>A_0=A</math>
+
#<math>A_0=A</math>.
 
#Пусть имеется матрица <math>A_k</math>, тогда матрица <math>A_{k+1}</math> строится следующим образом:
 
#Пусть имеется матрица <math>A_k</math>, тогда матрица <math>A_{k+1}</math> строится следующим образом:
−
::::*Строится QR-разложение: <math>A_k=Q_kR_k</math>
+
::::*Строится QR-разложение: <math>A_k=Q_kR_k</math>.
−
::::*Вычисляется <math>A_{k+1}=R_kQ_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}=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> и имеют одинаковый набор собственных значений.
 
Таким образом матрицы <math>A_{k+1}</math> и <math>A_{k}</math> подобны для <math>\forall k</math>, а значит, в силу транзитивности подобия, все матрицы <math>A_{k}</math> подобны матрице <math>A</math> и имеют одинаковый набор собственных значений.
Строка 90: Строка 92:
 
\Lambda_1 & 0\\  
 
\Lambda_1 & 0\\  
 
0 & \Lambda_2  
 
0 & \Lambda_2  
−
\end{bmatrix}, \Lambda_1\in\mathbb{C}^{m\times m},\Lambda_2\in\mathbb{C}^{r\times r} </math>
+
\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 | >\left | \lambda_{m+1}  \right | \geq ...\geq  \left | \lambda_{m+r}  \right | > 0 </math>, где <math>\{\lambda_1,...,\lambda_m\} = \lambda(\Lambda_1), \{\lambda_{m+1},...,\lambda_{m+r}\} = \lambda(\Lambda_2) </math>
+
::2. <math>\left | \lambda_1  \right | \geq ...\geq \left | \lambda_m  \right | >\left | \lambda_{m+1}  \right | \geq ...\geq  \left | \lambda_{m+r}  \right | > 0 </math>, где <math>\{\lambda_1,...,\lambda_m\} = \lambda(\Lambda_1), \{\lambda_{m+1},...,\lambda_{m+r}\} = \lambda(\Lambda_2) </math>.
  
  
−
::3. Ведущая подматрица порядка <math style="vertical-align:0%">m</math> в <math style="vertical-align:0%">X^{-1}</math> невырождена
+
::3. Ведущая подматрица порядка <math style="vertical-align:0%">m</math> в <math style="vertical-align:0%">X^{-1}</math> невырождена.
  
  
 
Тогда при <math style="vertical-align:0%"> k \rightarrow \infty </math> последовательность матриц <math>A_k</math> сходится к матрице с верхней треугольной формой<ref>Тыртышников Е.Е. "Методы численного анализа" — М., Академия, 2007. - 320 c.</ref>.
 
Тогда при <math style="vertical-align:0%"> k \rightarrow \infty </math> последовательность матриц <math>A_k</math> сходится к матрице с верхней треугольной формой<ref>Тыртышников Е.Е. "Методы численного анализа" — М., Академия, 2007. - 320 c.</ref>.
  
−
Таким образом, на практике необходимо выполнять итерации до тех пор пока матрица <math>A_k</math> не станет треугольной (также можно продолжать выполнять их пока искомая матрица не будет найдена с некоторой заранее заданной точностью <math>\varepsilon</math>). Если итерации закончились на шаге <math>N</math>, то числа на диагонали матрицы <math>A_N</math> будут считаться собственными значениями матрицы <math>A</math>
+
Таким образом, на практике необходимо выполнять итерации до тех пор пока матрица <math>A_k</math> не станет треугольной (также можно продолжать выполнять их пока искомая матрица не будет найдена с некоторой заранее заданной точностью <math>\varepsilon</math>). Если итерации закончились на шаге <math>N</math>, то числа на диагонали матрицы <math>A_N</math> будут считаться собственными значениями матрицы <math>A</math>.
  
 
==== Вещественный вариант QR-алгоритма ====
 
==== Вещественный вариант QR-алгоритма ====
  
−
Если вещественная матрица <math>A</math> имеет различные вещественные собственные значения, то, как было описано ранее, QR-алгоритм сходится к матрице с верхней треугольной формой, на диагонали которой находятся собственные значения. Однако вещественная матрица может иметь комплексные собственные значения. В данном случае алгоритм будет сходиться не к верхней треугольной матрице, а к блочной верхней треугольной матрице, которая на диагонали содержит блоки 1-го и 2-го порядка. Блоки 1-го порядка содержат различные вещественные собственные значения, блоки 2-го порядка соответствуют парам комплексных сопряженных собственных значений<ref>R. Granat, Bo Kagstrom, D. Kressner "LAPACK Working Note #216: A novel parallel QR algorithm for hybrid distributed memory HPC systems"</ref>.
+
Если вещественная матрица <math>A</math> имеет различные вещественные собственные значения, то, как было описано ранее, QR-алгоритм сходится к матрице с верхней треугольной формой, на диагонали которой находятся собственные значения. Однако вещественная матрица может иметь комплексные собственные значения. В данном случае алгоритм будет сходиться не к верхней треугольной матрице, а к блочной верхней треугольной матрице, которая на диагонали содержит блоки 1-го и 2-го порядка. Блоки 1-го порядка содержат различные вещественные собственные значения, блоки 2-го порядка соответствуют парам комплексных сопряженных собственных значений<ref>R. Granat, Bo Kagstrom, D. Kressner "LAPACK Working Note #216: A novel parallel QR algorithm for hybrid distributed memory HPC systems".</ref>.
  
 
<math>A_N=   
 
<math>A_N=   
Строка 118: Строка 120:
 
  \vdots&  \ddots&  \ddots&  \ddots&  \ddots&  \ddots&  \ddots&  \ddots& \bullet\\  
 
  \vdots&  \ddots&  \ddots&  \ddots&  \ddots&  \ddots&  \ddots&  \ddots& \bullet\\  
 
  0&  \cdots&  \cdots&  \cdots&  \cdots&  \cdots&  0&  0& \blacksquare
 
  0&  \cdots&  \cdots&  \cdots&  \cdots&  \cdots&  0&  0& \blacksquare
 +
\end{bmatrix}</math>.
 +
 +
В дальнейшем, в данной статье, матрицы, имеющие вышеописанную форму, будут называться квазитреугольными.
 +
 +
==== Приведение матрицы к форме Хессенберга ====
 +
 +
:''Основная статья:'' [[Метод Хаусхолдера (отражений) приведения матрицы к хессенберговой (почти треугольной) форме]]
 +
 +
Матрицей, имеющей форму Хессенберга, называется такая матрица, у которой все элементы, находящиеся ниже первой поддиагонали, равны нулю (<math>a_{ij}=0</math> при <math>i<j+1</math>). Пример такой матрицы приведен ниже:
 +
 +
<math>\begin{bmatrix}
 +
1&  2&  3& 4\\
 +
2&  5&  6& 7\\
 +
0&  3&  8& 9\\
 +
0&  0&  4& 1
 
\end{bmatrix}</math>
 
\end{bmatrix}</math>
  
−
В дальнейшем, в данной статье, матрицы, имеющие вышеописанную форму, будут называться квазитреугольными.
+
Любую матрицу <math>A</math> можно привести к подобной ей матрице <math>A_H</math>, имеющей форму Хессенберга (в силу подобия данные матрицы будут иметь одинаковые собственные значения). Наличие нулевых элементов в данной матрице позволяет ускорить процесс QR-разложения, причем данное ускорение будет иметь место на каждой итерации алгоритма, т.к. матрица с формой Хессенберга инвариантна относительно QR-итерации. Ускорение можно получить за счет использования модифицированного метода Гивенса(вращений) QR-разложения, который из-за изначального наличия нулевых элементов будет состоять не из <math>\frac{n(n-1)}{2}</math> шагов, а из <math>n-1</math> шагов (будет происходить домножение только на те матрицы вращения <math>T_{ij}</math>, у которых <math>i=j+1</math>).
 +
<!-- Одним из способов приведения матрицы к форме Хессенберга является преобразование Хаусхолдера<ref>[https://en.wikipedia.org/wiki/Householder_transformation#Tridiagonalization%7C Wikipedia: Householder transformation]</ref>. Данный алгоритм состоит из <math>n-1</math> итераций. Пусть имеется матрица <math>A^{(k)}</math>, тогда для перехода к матрице <math>A^{(k+1)}</math> проделывается следующее:
 +
 
 +
* Вычисление параметра <math>\alpha</math>: <math>\alpha=-sgn(a_{k+1,k}^{(k)})\sqrt{\sum_{j=k+1}^{n}(a_{j,k}^{(k)})^2}</math>
 +
 
 +
* Вычисление параметра <math>r</math>: <math>r=\sqrt{\frac{1}{2}(\alpha^2-a_{k+1,k}^{(k)}\alpha )}</math>
 +
 
 +
* Вычисление вектора <math>v^{(k)}=\begin{bmatrix}v_1\\ v_2\\ \vdots \\ v_n\end{bmatrix}</math>, где <math>v_{1}^{(k)}=v_{2}^{(k)}=...=v_{k}^{(k)}=0, v_{k+1}^{(k)}=\frac{a_{k+1,k}^{(k)}-\alpha}{2r}</math>, а остальные элементы <math>v_{j}^{(k)}=\frac{a_{jk}^{(k)}}{2r}</math>.
 +
 
 +
* Вычисление матрицы <math>P^{(k)}</math>: <math>P^{(k)}=E-2v^{(k)}(v^{(k)})^T</math>
 +
 
 +
* Вычисление матрицы <math>A^{(k+1)}</math>: <math>A^{(k+1)}=P^{(k)}A^{(k)}(P^{(k)})^T</math>. Важно отметить, что форма, которой обладают матрицы <math>P^{(k)}</math>, позволяет свести данную операцию перемножения матриц к менее затратным операциям перемножения матрицы на вектор. -->
 +
 
 +
Одним из способов приведения матрицы к форме Хессенберга является Метод Хаусхолдера (отражений)<ref>[[Метод_Хаусхолдера_(отражений)_приведения_матрицы_к_хессенберговой_(почти_треугольной)_форме|Метод Хаусхолдера (отражений) приведения матрицы к хессенберговой (почти треугольной) форме]]</ref>. Суть метода Хаусхолдера заключается в последовательном приведении матрицы <math style="vertical-align:0%">A</math> к форме Хессенберга при помощи домножения ее на так называемые матрицы отражения.
 +
 
 +
На <math style="vertical-align:0%">i</math>-ом шаге задача <math style="vertical-align:0%">i</math>-ой матрицы отражения заключается в обнулении поддиагональных элементов <math style="vertical-align:0%">i</math>-го столбца матрицы <math style="vertical-align:0%">A</math> (при этом столбцы левее <math style="vertical-align:0%">i</math>-го не изменяются). Таким образом, алгоритм состоит из <math style="vertical-align:0%">n-2</math> шагов, на каждом из которых вычисляется очередная матрица отражения, после чего найденное отражение применяется к матрице, которая является результатом предыдущего шага.
 +
 
 +
Матрица отражения имеет вид <math>E-\frac{1}{\gamma }ww^*</math>, где вектор <math>w</math> вычисляется из текущего <math style="vertical-align:0%">i</math>-го столбца матрицы при помощи использования операции скалярного произведения векторов. Данное представление матриц отражения позволяет хранить их в виде одного вектора и сводит домножение матрицы отражения на текущую матрицу к арифметическим операциям над векторами (скалярное произведение и сложение векторов).
 +
 
 +
В основной статье "[[Метод Хаусхолдера (отражений) приведения матрицы к хессенберговой (почти треугольной) форме]]" данный алгоритм описан более подробно с той лишь разницей, что для QR-алгоритма нет необходимости обнулять наддиагональные элементы, как это сделано в вышеупомянутой статье.
 +
 
 +
==== QR-алгоритм со сдвигами ====
 +
 
 +
QR-алгоритм со сдвигами позволяет сократить количество итераций, необходимых для сходимости<ref name="bahvalov">Бахвалов Н.С., Жидков Н.П., Кобельков. Г.М. "Численный методы"— 6-е изд. — М. : БИНОМ. Лаборатория знаний, 2008. — 636 с.</ref>. Пусть у нас есть матрица <math>A_k</math>, тогда процесс перехода к матрице <math>A_{k+1}</math> выглядит следующим образом:
 +
 
 +
* На каждом шаге подбирается число <math>\nu_k</math> и ищется следующее QR-разложение: <math>A_k-\nu_kE=Q_kR_k</math>.
 +
* Вычисляется матрица <math>A_{k+1}</math>: <math>A_{k+1} = R_kQ_k+\nu_kE</math>.
 +
 
 +
При этом сохраняется свойство подобия матриц <math>A_k</math> и <math>A_{k+1}</math>:
 +
 
 +
<math>A_{k+1} = R_kQ_k+\nu_kE=Q_{k}^{T}Q_kR_kQ_k+\nu_kE=Q_{k}^{T}(A_k-\nu_kE)Q_k+\nu_kE=Q_{k}^{T}A_kQ_k-Q_{k}^{T}(\nu_kE)Q_k+\nu_kE=Q_{k}^{T}A_kQ_k-\nu_kE+\nu_kE=Q_{k}^{T}A_kQ_k</math>.
 +
 
 +
В качестве параметра <math>\nu_k</math> может быть взят последний диагональный элемент матрицы <math>A_k</math>.
 +
<!--===== Подбор параметра <math>\nu_k</math>=====
 +
 
 +
//TODO-->
  
 
=== Вычислительное ядро алгоритма ===
 
=== Вычислительное ядро алгоритма ===
 
----
 
----
  
 +
==== Базовый алгоритм ====
 
QR-алгоритм обладает двумя вычислительными ядрами, которые повторяются на каждой итерации:
 
QR-алгоритм обладает двумя вычислительными ядрами, которые повторяются на каждой итерации:
  
−
1. Вычисление QR-разложения матрицы: <math>A_k=Q_kR_k</math>. Существует несколько методов решения данной задачи:
+
# Вычисление QR-разложения матрицы: <math>A_k=Q_kR_k</math>. Существует несколько методов решения данной задачи:
−
* [[Метод_Хаусхолдера_(отражений)_QR-разложения_квадратной_матрицы,_вещественный_вариант|Метод Хаусхолдера (отражений)]]
+
#* [[Метод_Хаусхолдера_(отражений)_QR-разложения_квадратной_матрицы,_вещественный_вариант|Метод Хаусхолдера (отражений)]]: вычислительное ядро данного алгоритма состоит из операций скалярного произведения, необходимых для вычисления матрицы отражения, и из операций скалярного произведения, необходимых для пересчета матрицы на каждом шаге.
−
::Вычислительное ядро данного алгоритма состоит из операций скалярного произведения, необходимых для вычисления матрицы отражения, и из опреаций скалярного произведения, необходимых для пересчета матрицы на каждом шаге.
+
#* [[Метод_Гивенса_(вращений)_QR-разложения_квадратной_матрицы_(вещественный_вариант)|Метод Гивенса (вращений)]]: вычислительное ядро данного алгоритма состоит из численных арифметических операций, необходимых для вычисления параметров матрицы вращения, и из операций (эквивалентных перемножению комплексных чисел), необходимых для пересчета матрицы на каждом шаге.
−
* [[Метод_Гивенса_(вращений)_QR-разложения_квадратной_матрицы_(вещественный_вариант)|Метод Гивенса (вращений)]]
+
#[[Перемножение_плотных_неособенных_матриц_(последовательный_вещественный_вариант)|Перемножение]] двух плотных матриц: <math>A_{k+1}=R_kQ_k</math>.
−
::Вычислительное ядро данного алгоритма состоит из численных арифметических операций, необходимых для вычисления параметров матрицы вращения, и из операций (эквивалентных перемножению комплексных чисел), необходимых для пересчета матрицы на каждом шаге.
 
  
−
2. [[Перемножение_плотных_неособенных_матриц_(последовательный_вещественный_вариант)|Перемножение]] двух плотных матриц: <math>A_{k+1}=R_kQ_k</math>
+
==== Модифицированный алгоритм ====
 +
Модифицированный QR-алгоритм обладает тремя вычислительными ядрами, первое вычисляется единожды, второе и третье повторяются на каждой итерации:
 +
# Приведение изначальной матрицы к форме Хессенберга: вычислительное ядро данного алгоритма состоит из операций скалярного произведения, необходимых для вычисления матрицы отражения, и из операций скалярного произведения, необходимых для пересчета матрицы на каждом шаге.
 +
# Вычисление QR-разложения матрицы: <math>A_k-\nu_kE=Q_kR_k</math> при помощи модифицированного метода Гивенса (вращений). Описание вычислительного ядра не отличается от приведенного в описании базового алгоритма, единственная разница заключается в значительно меньшем количестве матриц вращения.
 +
# [[Перемножение_плотных_неособенных_матриц_(последовательный_вещественный_вариант)|Перемножение]] двух плотных матриц: <math>A_{k+1}-\nu_kE=R_kQ_k</math>.
  
 
=== Макроструктура алгоритма ===
 
=== Макроструктура алгоритма ===
 
----
 
----
  
−
Как уже было описано ранее, QR-алгоритм содержит в себе две макрооперации:
+
Как уже было описано ранее, базовый QR-алгоритм содержит в себе две макрооперации:
  
−
1.Вычисление QR-разложения матрицы: <math>A_k=Q_kR_k</math>.
+
#Вычисление QR-разложения матрицы: <math>A_k=Q_kR_k</math>.
 +
#Перемножение двух плотных матриц: <math>A_{k+1}=R_kQ_k</math>.
  
−
2.Перемножение двух плотных матриц: <math>A_{k+1}=R_kQ_k</math>
+
В случае модифицированного QR-алгоритма появляется еще макрооперация приведения матрицы к форме Хессенберга.
  
 
=== Схема реализации последовательного алгоритма ===
 
=== Схема реализации последовательного алгоритма ===
 
----
 
----
 +
 +
==== Базовый алгоритм ====
  
 
Опишем необходимый для реализации цикл при помощи псевдокода:
 
Опишем необходимый для реализации цикл при помощи псевдокода:
  
−
  A - исходная матрица
+
  A - исходная матрица.
−
  curA - текущая матрица, на основе которой будет вычисляться QR-разложение на каждом шаге
+
  curA - текущая матрица, на основе которой будет вычисляться QR-разложение на каждом шаге.
−
  nextA - новая матрица, полученная после перемножения матриц R и Q  
+
  nextA - новая матрица, полученная после перемножения матриц R и Q.
 
  triangular - функция, проверяющая, имеет ли матрица квазитреугольную форму.
 
  triangular - функция, проверяющая, имеет ли матрица квазитреугольную форму.
 
  difference - функция, проверяющая, что элементы двух матриц, стоящие на одинаковых местах, различаются не более чем на eps (данная функция нужна чтобы проверять не только сходимость матрицы к квазитреугольной форме, но и сходимость ее ненулевых элементов).
 
  difference - функция, проверяющая, что элементы двух матриц, стоящие на одинаковых местах, различаются не более чем на eps (данная функция нужна чтобы проверять не только сходимость матрицы к квазитреугольной форме, но и сходимость ее ненулевых элементов).
Строка 162: Строка 221:
 
     findQRdecomposition(curA,Q,R)
 
     findQRdecomposition(curA,Q,R)
 
     nextA = R*Q
 
     nextA = R*Q
 +
}
 +
 +
==== Модифицированный алгоритм ====
 +
 +
Опишем необходимый для реализации цикл при помощи псевдокода:
 +
 +
A - исходная матрица.
 +
curA - текущая матрица, на основе которой будет совершаться сдвиг и вычисляться QR-разложение на каждом шаге.
 +
nextA - новая матрица, полученная после перемножения матриц R и Q.
 +
nu - параметр сдвига
 +
shiftA - текущая матрица со сдвигом.
 +
triangular - функция, проверяющая, имеет ли матрица квазитреугольную форму.
 +
difference - функция, проверяющая, что элементы двух матриц, стоящие на одинаковых местах, различаются не более чем на eps (данная функция нужна чтобы проверять не только сходимость матрицы к квазитреугольной форме, но и сходимость ее ненулевых элементов).
 +
 +
curA = A
 +
nextA = A
 +
'''while''' ( '''not''' (triangular(nextA) & difference(curA,nextA,eps)) )
 +
{
 +
    curA = nextA
 +
    nu = computeNu(curA)
 +
    shiftA = curA - nu*E
 +
    findQRdecomposition(shiftA,Q,R)
 +
    nextA = R*Q + nu*E
 
  }
 
  }
  
 
=== Последовательная сложность алгоритма ===
 
=== Последовательная сложность алгоритма ===
 
----
 
----
−
Подсчитаем сложность одной итерации QR-алгоритма (расчет сложностей для QR-разложения и перемножения матриц представлен в статьях, посвященных данным алгоритмам)
 
  
−
# QR-разложение матрицы
+
==== Базовый алгоритм ====
−
#* [[Метод_Хаусхолдера_(отражений)_QR-разложения_квадратной_матрицы,_вещественный_вариант|Метод Хаусхолдера (отражений)]] имеет сложность  <math>\frac{4}{3}n^3</math>
+
 
−
#* [[Метод_Гивенса_(вращений)_QR-разложения_квадратной_матрицы_(вещественный_вариант)|Метод Гивенса (вращений)]] имеет сложность  <math style="vertical-align:0%">2n^3</math>
+
Подсчитаем сложность одной итерации QR-алгоритма (расчет сложностей для QR-разложения и перемножения матриц представлен в статьях, посвященных данным алгоритмам).
−
# [[Перемножение_плотных_неособенных_матриц_(последовательный_вещественный_вариант)|Перемножение двух плотных матриц]] имеет сложность  <math style="vertical-align:0%">n^3</math>
+
 
−
# Проверка матрицы на квазитреугольную форму состоит из набора сравнений элементов с номерами <math>ij</math>  (<math>i>j+1</math>) с нулем (таких элементов <math>\frac{(n-2)(n-1)}{2}</math>). Поддиагональные элементы с номерами <math>ij</math>  (<math>i=j+1</math>) должны быть проверены на соответствие необходимой квазитреугольной форме. Для этого достаточно для каждого такого элемента проверить следующее условие: <math>a_{i+1,i}==0 \vee a_{i,i-1}==0  \wedge a_{i+2,i+1}==0</math> (либо поддиагональный элемент равен 0, либо, в противном случае, соседние поддиагональные элементы равны 0, чтобы текущий элемент соответствовал блоку 2-го порядка). Таких проверок поддиагональных элементов будет <math>n-1</math>. После данных проверок следует набор логических операций "И" между результатами всех сравнений (таких операций будет <math>\frac{n(n-1)}{2}-1</math>)
+
# QR-разложение матрицы.
−
# Сравнение новой матрицы с предыдущей состоит из операций вычитания и сравнения для каждой пары ненулевых соответствующих элементов (такие пары имеют номера элементов <math>ij</math> (<math>i \geq j+1</math>), количество таких пар равно <math>\frac{n(n+1)}{2}+(n-1)</math>), а также из набора логических операций "И" между результатами сравнений (таких операций будет <math>\frac{n(n+1)}{2}+(n-1)-1</math>)
+
#* [[Метод_Хаусхолдера_(отражений)_QR-разложения_квадратной_матрицы,_вещественный_вариант|Метод Хаусхолдера (отражений)]] имеет сложность  <math>\frac{4}{3}n^3</math>.
 +
#* [[Метод_Гивенса_(вращений)_QR-разложения_квадратной_матрицы_(вещественный_вариант)|Метод Гивенса (вращений)]] имеет сложность  <math style="vertical-align:0%">2n^3</math>.
 +
# [[Перемножение_плотных_неособенных_матриц_(последовательный_вещественный_вариант)|Перемножение двух плотных матриц]] имеет сложность  <math style="vertical-align:0%">n^3</math>.
 +
# Проверка матрицы на квазитреугольную форму состоит из набора сравнений элементов с номерами <math>ij</math>  (<math>i>j+1</math>) с нулем (таких элементов <math>\frac{(n-2)(n-1)}{2}</math>). Поддиагональные элементы с номерами <math>ij</math>  (<math>i=j+1</math>) должны быть проверены на соответствие необходимой квазитреугольной форме. Для этого достаточно для каждого такого элемента проверить следующее условие: <math>a_{i+1,i}==0 \vee a_{i,i-1}==0  \wedge a_{i+2,i+1}==0</math> (либо поддиагональный элемент равен 0, либо, в противном случае, соседние поддиагональные элементы равны 0, чтобы текущий элемент соответствовал блоку 2-го порядка). Таких проверок поддиагональных элементов будет <math>n-1</math>. После данных проверок следует набор логических операций "И" между результатами всех сравнений (таких операций будет <math>\frac{n(n-1)}{2}-1</math>).
 +
# Сравнение новой матрицы с предыдущей состоит из операций вычитания и сравнения для каждой пары ненулевых соответствующих элементов (такие пары имеют номера элементов <math>ij</math> (<math>i \geq j+1</math>), количество таких пар равно <math>\frac{n(n+1)}{2}+(n-1)</math>), а также из набора логических операций "И" между результатами сравнений (таких операций будет <math>\frac{n(n+1)}{2}+(n-1)-1</math>).
 +
 
 +
Итого, в сумме получаем <math>O(n^3)</math> - сложность алгоритма на каждой итерации. Если алгоритм остановился на итерации с номером <math>N</math>, то общая сложность алгоритма будет равна <math>N*O(n^3)</math>.
 +
 
 +
==== Модифицированный алгоритм ====
 +
 
 +
Отличия от базового алгоритма заключаются в следующем:
 +
# [[Метод Хаусхолдера (отражений) приведения матрицы к хессенберговой (почти треугольной) форме|Метод Хаусхолдера (отражений)]] имеет последовательную сложность <math>O(n^3)</math>.
 +
# QR-разложение матрицы, имеющей форму Хессенберга будет иметь сложность порядка <math>O(n^2)</math><ref name="bahvalov"/>.
 +
 
 +
Стоит заметить, что вычисление сдвигов по сравнению с другими операциями не обладает высокой сложностью, однако позволяет значительно сократить количество необходимых итераций.
  
−
Итого, в сумме получаем <math>O(n^3)</math> - сложность алгоритма на каждой итерации. Если алгоритм остановился на итерации с номером <math>N</math>, то общая сложность алгоритма будет равна <math>N*O(n^3)</math>
+
Итого, в сумме получаем <math>O(n^2)</math> - сложность алгоритма на каждой итерации. Сложность всего алгоритма, с учетом приведения матрицы к форме Хессенберга, будет равна <math>O(n^3)+N*O(n^2)</math>
  
 
=== Информационный граф ===
 
=== Информационный граф ===
 
----
 
----
 +
==== Базовый алгоритм ====
 +
На рисунках 1 и 2 изображен информационный граф QR-алгоритма. Вершины данного графа обозначают следующее:
 +
* '''QR''' - операция QR-разложения матрицы.
 +
* ''' *<sub>M</sub>''' - операция перемножения матриц.
 +
* '''Triang''' - операция проверки матрицы на квазитреугольную форму.
 +
* '''Dif''' - операция проверки различия элементов двух матриц не более чем на некоторое заданное число.
 +
* '''Iter<sub>i</sub>''' - итерация алгоритма
  
−
На рисунке ниже изображен информационный граф QR-алгоритма. Вершины данного графа обозначают следующее:
+
[[Файл:QR_Alg_Base_All.png|center|frame|Рисунок 1. Информационный граф базового QR-алгоритма]]
−
* '''QR''' - операция QR-разложения матрицы
+
[[Файл:QR_Alg_Base_Iter.png|center|frame|Рисунок 2. Информационный граф итерации базового QR-алгоритма]]
−
* ''' *<sub>M</sub>''' - операция перемножения матриц
 
−
* '''Triang''' - операция проверки матрицы на квазитреугольную форму
 
−
* '''Dif''' - операция проверки различия элементов двух матриц не более чем на некоторое заданное число
 
  
−
[[Файл:QR_Alg_Graf.png|center|frame|Информационный граф QR-алгоритма]]
+
Подробные графы операций QR-разложения ([[Метод_Хаусхолдера_(отражений)_QR-разложения_квадратной_матрицы,_вещественный_вариант|Метод Хаусхолдера (отражений)]], [[Метод_Гивенса_(вращений)_QR-разложения_квадратной_матрицы_(вещественный_вариант)|Метод Гивенса (вращений)]]) и [[Перемножение_плотных_неособенных_матриц_(последовательный_вещественный_вариант)|перемножения]] матриц можно найти в статьях, посвященных этим алгоритмам. Далее рассмотрим информационные графы операций '''Triang''' (рис.3) и '''Dif''' (рис.4) на примере матрицы размера <math style="vertical-align:0%">5 \times 5</math>. Графы для других матриц выглядят аналогичным образом.
  
−
Подробные графы операций QR-разложения ([[Метод_Хаусхолдера_(отражений)_QR-разложения_квадратной_матрицы,_вещественный_вариант|Метод Хаусхолдера (отражений)]], [[Метод_Гивенса_(вращений)_QR-разложения_квадратной_матрицы_(вещественный_вариант)|Метод Гивенса (вращений)]]) и [[Перемножение_плотных_неособенных_матриц_(последовательный_вещественный_вариант)|перемножения]] матриц можно найти в статьях, посвященных этим алгоритмам. Далее рассмотрим информационные графы операций '''Triang''' и '''Dif''' на примере матрицы размера <math style="vertical-align:0%">5 \times 5</math>. Графы для других матриц выглядят аналогичным образом.
+
[[Файл:QR_Alg_Graf_Triang.png|center|frame|Рисунок 3. Информационный граф операции '''Triang''', которая проверяет соответствие квазитреугольной форме]]
−
 
 
−
[[Файл:QR_Alg_Graf_Triang.png|center|frame|Информационный граф операции '''Triang''', которая проверяет соответствие квазитреугольной форме]]
 
  
 
Вершины '''V''' соответствуют операции сравнения с 0. Вершины '''F''' соответствуют проверке поддиагональных элементов на соответствие квазитреугольной форме, которая была описана в предыдущем разделе. Вершины '''&''' соответствуют логической операции "И".
 
Вершины '''V''' соответствуют операции сравнения с 0. Вершины '''F''' соответствуют проверке поддиагональных элементов на соответствие квазитреугольной форме, которая была описана в предыдущем разделе. Вершины '''&''' соответствуют логической операции "И".
  
−
[[Файл:QR_Alg_Graf_Dif2.png|center|frame|Информационный граф операции '''Dif''', которая проверяет сходимость ненулевых элементов]]
+
[[Файл:QR_Alg_Graf_Dif2.png|center|frame|Рисунок 4. Информационный граф операции '''Dif''', которая проверяет сходимость ненулевых элементов]]
  
 
Вершины '''-V''' соответствуют операции вычитания и сравнения с 0. Вершины '''&''' соответствуют логической операции "И".
 
Вершины '''-V''' соответствуют операции вычитания и сравнения с 0. Вершины '''&''' соответствуют логической операции "И".
 +
 +
==== Модифицированный алгоритм ====
 +
 +
На рисунках 5 и 6 изображен информационный граф модифицированного QR-алгоритма. Дополнительные вершины для этого графа обозначают следующее:
 +
* '''Shift''' - операция сдвига
 +
* '''H''' - операция приведения матрицы к форме Хессенберга
 +
 +
[[Файл:QR_Alg_Mod_All.png|center|frame|Рисунок 5. Информационный граф модифицированного QR-алгоритма]]
 +
[[Файл:QR_Alg_Mod_Iter.png|center|frame|Рисунок 6. Информационный граф итерации модифицированного QR-алгоритма]]
  
 
=== Ресурс параллелизма алгоритма ===
 
=== Ресурс параллелизма алгоритма ===
 
----
 
----
 +
==== Базовый алгоритм ====
  
 
На макроуровне (который можно увидеть на информационном графе QR-алгоритма) алгоритм не обладает ресурсами параллелизма. Все макрооперации внутри итерации, а также сами итерации выполняются последовательно (за исключением операций '''Triag''' и '''Dif''', которые могут выполняться параллельно). Основной ресурс параллелизма заложен отдельно в каждой из макроопераций. На каждой итерации алгоритм имеет следующие параллельные характеристики:
 
На макроуровне (который можно увидеть на информационном графе QR-алгоритма) алгоритм не обладает ресурсами параллелизма. Все макрооперации внутри итерации, а также сами итерации выполняются последовательно (за исключением операций '''Triag''' и '''Dif''', которые могут выполняться параллельно). Основной ресурс параллелизма заложен отдельно в каждой из макроопераций. На каждой итерации алгоритм имеет следующие параллельные характеристики:
  
−
# QR-разложение матрицы (описание ресурсов параллелизма для алгоритмов QR-разложения можно найти в статьях, посвященных этим алгоритмам)
+
# QR-разложение матрицы (описание ресурсов параллелизма для алгоритмов QR-разложения можно найти в статьях, посвященных этим алгоритмам).
−
#* [[Метод_Хаусхолдера_(отражений)_QR-разложения_квадратной_матрицы,_вещественный_вариант|Метод Хаусхолдера (отражений)]] имеет высоту ярусно-параллельной формы <math>O(n^2)</math> и ширину ярусно-параллельной формы <math>O(n)</math>
+
#* [[Метод_Хаусхолдера_(отражений)_QR-разложения_квадратной_матрицы,_вещественный_вариант|Метод Хаусхолдера (отражений)]] имеет высоту ярусно-параллельной формы <math>O(n^2)</math> и ширину ярусно-параллельной формы <math>O(n)</math>.
−
#* [[Метод_Гивенса_(вращений)_QR-разложения_квадратной_матрицы_(вещественный_вариант)|Метод Гивенса (вращений)]] имеет высоту ярусно-параллельной формы <math>11n-16</math> и ширину ярусно-параллельной формы <math>O(n^2)</math>
+
#* [[Метод_Гивенса_(вращений)_QR-разложения_квадратной_матрицы_(вещественный_вариант)|Метод Гивенса (вращений)]] имеет высоту ярусно-параллельной формы <math>11n-16</math> и ширину ярусно-параллельной формы <math>O(n^2)</math>.
−
# [[Перемножение_плотных_неособенных_матриц_(последовательный_вещественный_вариант)|Перемножение двух плотных матриц]] имеет высоту ярусно-параллельной формы <math>n</math> и ширину ярусно-параллельной формы <math style="vertical-align:0%">n^2</math>
+
# [[Перемножение_плотных_неособенных_матриц_(последовательный_вещественный_вариант)|Перемножение двух плотных матриц]] имеет высоту ярусно-параллельной формы <math>n</math> и ширину ярусно-параллельной формы <math style="vertical-align:0%">n^2</math>.
−
# Проверка матрицы на выходе из итерации
+
# Проверка матрицы на выходе из итерации.
−
#* Проверка матрицы на квазитреугольную форму состоит из одного яруса сравнений для каждого элемента и последующих ярусов, вычисляющих итоговый результат при помощи логической операции "И". Для вычисления итогового результата можно использовать метод сдваивания, который дает высоту ярусно-параллельной формы порядка логарифма от количества элементов, над которыми совершается операция. Таким образов высота ярусно-параллельной формы будет равна <math>O(log_2n)</math>. Ширина ярусно-параллельной формы достигается на ярусе ссравнений для каждого элемента и равна <math>O(n^2)</math>.
+
#* Проверка матрицы на квазитреугольную форму состоит из одного яруса сравнений для каждого элемента и последующих ярусов, вычисляющих итоговый результат при помощи логической операции "И". Для вычисления итогового результата можно использовать метод сдваивания, который дает высоту ярусно-параллельной формы порядка логарифма от количества элементов, над которыми совершается операция. Таким образов высота ярусно-параллельной формы будет равна <math>O(log_2n)</math>. Ширина ярусно-параллельной формы достигается на ярусе сравнений для каждого элемента и равна <math>O(n^2)</math>.
 
#* Сравнение новой матрицы с предыдущей состоит из одного яруса сравнений для каждого элемента и последующих ярусов, вычисляющих итоговый результат при помощи логической операции "И". Высота ярусно-параллельной формы будет равна <math>O(log_2n)</math>. Ширина ярусно-параллельной формы достигается на ярусе сравнений для каждого элемента и равна <math>O(n^2)</math>.
 
#* Сравнение новой матрицы с предыдущей состоит из одного яруса сравнений для каждого элемента и последующих ярусов, вычисляющих итоговый результат при помощи логической операции "И". Высота ярусно-параллельной формы будет равна <math>O(log_2n)</math>. Ширина ярусно-параллельной формы достигается на ярусе сравнений для каждого элемента и равна <math>O(n^2)</math>.
  
 
Таким образом, основной вклад в высоту ярусно-параллельной формы одной итерации вносит операция QR-разложения матрицы и она будет равна <math>O(n)</math>, если использовать метод Гивенса. Ширина ярусно-параллельной формы будет равна <math>O(n^2)</math>. Если алгоритм остановился на итерации с номером <math>N</math>, то параллельные характеристики для всего алгоритма будут равны <math>N*O(n)</math> для высоты и <math>O(n^2)</math> для ширины ярусно-параллельной формы.
 
Таким образом, основной вклад в высоту ярусно-параллельной формы одной итерации вносит операция QR-разложения матрицы и она будет равна <math>O(n)</math>, если использовать метод Гивенса. Ширина ярусно-параллельной формы будет равна <math>O(n^2)</math>. Если алгоритм остановился на итерации с номером <math>N</math>, то параллельные характеристики для всего алгоритма будут равны <math>N*O(n)</math> для высоты и <math>O(n^2)</math> для ширины ярусно-параллельной формы.
 +
 +
==== Модифицированный алгоритм ====
 +
 +
Отличия от базового алгоритма заключаются в следующем:
 +
# Приведение матрицы к форме Хессенберга имеет высоту ярусно-параллельной формы <math>O(n^2)</math> и ширину ярусно-параллельной формы <math>O(n^2)</math>.
 +
# QR-разложение матрицы модифицированным методом Гивенса (вращений) имеет высоту ярусно-параллельной формы <math>O(n)</math> и ширину ярусно-параллельной формы <math>O(n)</math>.
 +
# Операция сдвига (вычитание из диагональных элементов матрицы одного числа) имеет высоту ярусно-параллельной формы <math>O(1)</math> и ширину ярусно-параллельной формы <math>O(n)</math>.
 +
 +
Таким образом, итоговые характеристики всего модифицированного алгоритма будут следующие: высота ярусно-параллельной формы <math>O(n^2)+N*O(n)</math> и ширина ярусно-параллельной формы <math>O(n^2)</math>.
  
 
=== Входные и выходные данные алгоритма ===
 
=== Входные и выходные данные алгоритма ===
Строка 217: Строка 333:
  
 
'''Входные данные''':  
 
'''Входные данные''':  
−
:::квадратная вещественная плотная матрица <math>A</math>: <math>A \in \R^{n \times n}</math>
+
:::квадратная вещественная плотная матрица <math>A</math>: <math>A \in \R^{n \times n}</math>.
  
 
'''Объем входных данных''':  
 
'''Объем входных данных''':  
−
:::<math>n^2</math> (необходимо хранить все элементы матрицы)
+
:::<math>n^2</math> (необходимо хранить все элементы матрицы).
  
 
'''Выходные данные''':  
 
'''Выходные данные''':  
−
:::собственные значения матрицы <math>A</math>
+
:::собственные значения матрицы <math>A</math>.
  
 
'''Объем выходных данных''':  
 
'''Объем выходных данных''':  
−
:::<math>n</math> (квадратная матрица размера <math>n \times n</math> имеет ровно <math>n</math> собственных значений при этом некоторые из них могут быть комплексными)
+
:::<math>n</math> (квадратная матрица размера <math>n \times n</math> имеет ровно <math>n</math> собственных значений при этом некоторые из них могут быть комплексными).
  
 
=== Свойства алгоритма ===
 
=== Свойства алгоритма ===
Строка 235: Строка 351:
 
*Скорость сходимости алгоритма зависит от собственных значений. Чем ближе по модулю соседние собственные значения, тем меньше скорость сходимости. Этот недостаток призван решить QR-алгоритм со сдвигами.
 
*Скорость сходимости алгоритма зависит от собственных значений. Чем ближе по модулю соседние собственные значения, тем меньше скорость сходимости. Этот недостаток призван решить QR-алгоритм со сдвигами.
  
−
== ЧАСТЬ. Программная реализация алгоритма ==
+
== Программная реализация алгоритма ==
  
 
=== Особенности реализации последовательного алгоритма ===
 
=== Особенности реализации последовательного алгоритма ===
Строка 248: Строка 364:
 
=== Масштабируемость алгоритма и его реализации ===
 
=== Масштабируемость алгоритма и его реализации ===
 
----
 
----
 +
Масштабируемость алгоритма (рис.7) исследовалась на основе реализации алгоритма в библиотеке [http://www.netlib.org/scalapack/#_scalapack_version_2_0_0 ScaLAPACK v2.0.0]: функция [http://www.netlib.org/scalapack/explore-html/d3/d3a/pdhseqr_8f.html#a2dbf17fffb167fe38d2a0544bb9fd5db pdhseqr]. Время работы было замерено для различных размеров матрицы, а также различного числа процессоров.
 +
 +
[[Файл:QRPerfomance.png|center|frame|Рисунок 7. Масштабируемость QR-алгоритма]]
 +
 +
При увеличении числа процессоров время работы программы увеличивается при малом размере матрицы, распараллеливание алгоритма неэффективно. Однако при больших размерах матрицы (более <math style="vertical-align:0%"> 10^3\times10^3</math>) время работы параллельной реализации значительно меньше последовательной.
 +
 +
[https://algowiki-project.org/pool/File:QR-Kiyamova-Smirnova.docx Код программы]
 +
 +
Все вычисления были проведены на суперкомпьютере [http://hpc.cs.msu.su/regatta Regatta]. Библиотека [http://www.netlib.org/scalapack/#_scalapack_version_2_0_0 ScaLAPACK v2.0.0] была установлена вручную.
 +
Компиляция проводилась с использованием компилятора mpicxx с опциями -lscalapack и -llapack с указанием полного пути до используемых библиотек.
  
 
=== Динамические характеристики и эффективность реализации алгоритма ===
 
=== Динамические характеристики и эффективность реализации алгоритма ===
Строка 261: Строка 387:
  
 
# LAPACK ([http://netlib.org/lapack/ Linear Algebra PACKage]) — библиотека с открытым исходным кодом, содержащая методы для решения основных задач линейной алгебры. Написана на языке Fortran с использованием библиотеки [http://www.netlib.org/blas/ BLAS] и является развитием пакета LINPACK. Находится в открытом доступе в соответствии с [http://www.netlib.org/lapack/LICENSE.txt модифицированной лицензией BSD], в том числе и для коммерческого использования.  
 
# LAPACK ([http://netlib.org/lapack/ Linear Algebra PACKage]) — библиотека с открытым исходным кодом, содержащая методы для решения основных задач линейной алгебры. Написана на языке Fortran с использованием библиотеки [http://www.netlib.org/blas/ BLAS] и является развитием пакета LINPACK. Находится в открытом доступе в соответствии с [http://www.netlib.org/lapack/LICENSE.txt модифицированной лицензией BSD], в том числе и для коммерческого использования.  
−
#* Полный алгоритм нахождения собственных значений: функция [http://www.netlib.org/lapack/explore-html/d8/d66/dhseqr_8f_a273e6c6bdb731ed8d6ba3cb060f3f739.html#a273e6c6bdb731ed8d6ba3cb060f3f739 dhseqr]  
+
#* Полный алгоритм нахождения собственных значений: функция [http://www.netlib.org/lapack/explore-html/d8/d66/dhseqr_8f_a273e6c6bdb731ed8d6ba3cb060f3f739.html#a273e6c6bdb731ed8d6ba3cb060f3f739 dhseqr].
−
#* QR-разложение: функция [http://www.netlib.org/lapack/explore-html/dd/d9a/group__double_g_ecomputational_ga3766ea903391b5cf9008132f7440ec7b.html#ga3766ea903391b5cf9008132f7440ec7b dgeqrf]  
+
#* QR-разложение: функция [http://www.netlib.org/lapack/explore-html/dd/d9a/group__double_g_ecomputational_ga3766ea903391b5cf9008132f7440ec7b.html#ga3766ea903391b5cf9008132f7440ec7b dgeqrf].
−
#* Перемножение матриц: функция [http://www.netlib.org/lapack/explore-html/d1/d54/group__double__blas__level3_gaeda3cbd99c8fb834a60a6412878226e1.html#gaeda3cbd99c8fb834a60a6412878226e1 dgemm]  
+
#* Перемножение матриц: функция [http://www.netlib.org/lapack/explore-html/d1/d54/group__double__blas__level3_gaeda3cbd99c8fb834a60a6412878226e1.html#gaeda3cbd99c8fb834a60a6412878226e1 dgemm].
 
# [http://alglib.sources.ru/ ALGLIB] - это кросс-платформенная библиотека численного анализа, поддерживающая несколько языков программирования (C++, C#, Pascal, VBA) и несколько операционных систем (Windows, Linux, Solaris). ALGLIB является свободным программным обеспечением, которое использует двойное лицензирование: оно может быть использовано в соответствии с лицензией GPL (версии 2+), а для использования в коммерческих целях необходимо купить отдельную лицензию.
 
# [http://alglib.sources.ru/ ALGLIB] - это кросс-платформенная библиотека численного анализа, поддерживающая несколько языков программирования (C++, C#, Pascal, VBA) и несколько операционных систем (Windows, Linux, Solaris). ALGLIB является свободным программным обеспечением, которое использует двойное лицензирование: оно может быть использовано в соответствии с лицензией GPL (версии 2+), а для использования в коммерческих целях необходимо купить отдельную лицензию.
−
#* Полный алгоритм нахождения собственных значений: подпрограмма [http://alglib.sources.ru/eigen/nonsymmetric/nonsymmetricevd.php RMatrixEVD]  
+
#* Полный алгоритм нахождения собственных значений: подпрограмма [http://alglib.sources.ru/eigen/nonsymmetric/nonsymmetricevd.php RMatrixEVD].
−
#* QR-разложение: подпрограмма [http://alglib.sources.ru/matrixops/qr.php  rmatrixqr] реализует QR-разложение для вещественных матриц, [http://alglib.sources.ru/matrixops/qr.php cmatrixqr] – для комплексных матриц  
+
#* QR-разложение: подпрограмма [http://alglib.sources.ru/matrixops/qr.php  rmatrixqr] реализует QR-разложение для вещественных матриц, [http://alglib.sources.ru/matrixops/qr.php cmatrixqr] – для комплексных матриц.
−
#* Перемножение матриц: для перемножения матриц ALGLIB использует соответствующую реализацию библиотеки [http://www.netlib.org/blas/ BLAS]
+
#* Перемножение матриц: для перемножения матриц ALGLIB использует соответствующую реализацию библиотеки [http://www.netlib.org/blas/ BLAS].
 
# [http://eigen.tuxfamily.org/index.php?title=Main_Page Eigen] – библиотека шаблонов для линейной алгебры, написанная на языке C++. Eigen – свободно распространяемое программное обеспечение. Начиная с версии 3.1.1, оно лицензируется MPL2, на ранние версии распространяется действие лицензии LGPL3+.
 
# [http://eigen.tuxfamily.org/index.php?title=Main_Page Eigen] – библиотека шаблонов для линейной алгебры, написанная на языке C++. Eigen – свободно распространяемое программное обеспечение. Начиная с версии 3.1.1, оно лицензируется MPL2, на ранние версии распространяется действие лицензии LGPL3+.
−
#* Полный алгоритм нахождения собственных значений: модуль [http://eigen.tuxfamily.org/dox/group__Eigenvalues__Module.html Eigenvalues module]  
+
#* Полный алгоритм нахождения собственных значений: модуль [http://eigen.tuxfamily.org/dox/group__Eigenvalues__Module.html Eigenvalues module].
−
#* QR-разложение: модуль [http://eigen.tuxfamily.org/dox-devel/group__QR__Module.html  QR module]
+
#* QR-разложение: модуль [http://eigen.tuxfamily.org/dox-devel/group__QR__Module.html  QR module].
−
#* Перемножение матриц: [http://eigen.tuxfamily.org/dox/group__TutorialMatrixArithmetic.html#title5  операторы * и *=]
+
#* Перемножение матриц: [http://eigen.tuxfamily.org/dox/group__TutorialMatrixArithmetic.html#title5  операторы * и *=].
  
 
==== Параллельные реализации ====
 
==== Параллельные реализации ====
 
# ScaLAPACK ([http://www.netlib.org/scalapack/ Scalable Linear Algebra PACKage]) — библиотека с открытым исходным кодом, включающая в себя подмножество процедур [http://netlib.org/lapack/ LAPACK], переработанных для использования на MPP-компьютерах. ScaLAPACK разработана с использованием PBLAS и BLACS, и предназначена для вычислений на любом компьютере или кластере поддерживающим MPI или PVM. Библиотека в настоящее время написана на языке Fortran. Находится в открытом доступе в соответствии с [http://www.netlib.org/scalapack/LICENSE модифицированной лицензией BSD], в том числе и для коммерческого использования.
 
# ScaLAPACK ([http://www.netlib.org/scalapack/ Scalable Linear Algebra PACKage]) — библиотека с открытым исходным кодом, включающая в себя подмножество процедур [http://netlib.org/lapack/ LAPACK], переработанных для использования на MPP-компьютерах. ScaLAPACK разработана с использованием PBLAS и BLACS, и предназначена для вычислений на любом компьютере или кластере поддерживающим MPI или PVM. Библиотека в настоящее время написана на языке Fortran. Находится в открытом доступе в соответствии с [http://www.netlib.org/scalapack/LICENSE модифицированной лицензией BSD], в том числе и для коммерческого использования.
−
#* Полный алгоритм нахождения собственных значений: функция [http://www.netlib.org/scalapack/explore-html/d3/d3a/pdhseqr_8f.html#a2dbf17fffb167fe38d2a0544bb9fd5db pdhseqr]
+
#* Полный алгоритм нахождения собственных значений: функция [http://www.netlib.org/scalapack/explore-html/d3/d3a/pdhseqr_8f.html#a2dbf17fffb167fe38d2a0544bb9fd5db pdhseqr].
−
#* QR-разложение: функция [http://www.netlib.org/scalapack/explore-html/da/d85/pdgeqrf_8f.html#ad48b15b5976e51cbc40a0e1ab2720fa5 pdgeqrf]
+
#* QR-разложение: функция [http://www.netlib.org/scalapack/explore-html/da/d85/pdgeqrf_8f.html#ad48b15b5976e51cbc40a0e1ab2720fa5 pdgeqrf].
−
#* Перемножение матриц: функция [http://www.netlib.org/scalapack/explore-html/d6/da2/pdgemm___8c.html#a2959b3974ab5e83564517aa1b99d2433 pdgemm]
+
#* Перемножение матриц: функция [http://www.netlib.org/scalapack/explore-html/d6/da2/pdgemm___8c.html#a2959b3974ab5e83564517aa1b99d2433 pdgemm].
 
# PLAPACK ([http://www.cs.utexas.edu/users/plapack/ Parallel Linear Algebra Package]) — пакет функций LAPACK для параллельного решения задач линейной алгебры. Пакет функций PLAPACK является альтернативой библиотеке [http://www.netlib.org/scalapack/ ScaLAPACK]. Для осуществления межпроцессорных коммуникаций в PLAPACK использован интерфейс передачи сообщений MPI. При передаче сообщений в PLAPACK в основном используются коллективные операции, такие, как обобщенная передача данных от одного процесса всем процессам  (MPI_Scatter), обобщенная передача данных от всех процессов одному процессу (MPI_Gather), широковещательная рассылка (MPI_Bcast) и другие. PLAPACK включает интерфейсы для языков Fortran и С.  
 
# PLAPACK ([http://www.cs.utexas.edu/users/plapack/ Parallel Linear Algebra Package]) — пакет функций LAPACK для параллельного решения задач линейной алгебры. Пакет функций PLAPACK является альтернативой библиотеке [http://www.netlib.org/scalapack/ ScaLAPACK]. Для осуществления межпроцессорных коммуникаций в PLAPACK использован интерфейс передачи сообщений MPI. При передаче сообщений в PLAPACK в основном используются коллективные операции, такие, как обобщенная передача данных от одного процесса всем процессам  (MPI_Scatter), обобщенная передача данных от всех процессов одному процессу (MPI_Gather), широковещательная рассылка (MPI_Bcast) и другие. PLAPACK включает интерфейсы для языков Fortran и С.  
−
#* QR-разложение: функция [http://www.cs.utexas.edu/users/plapack/C_Manual/#PLA_QR PLA_QR]
+
#* QR-разложение: функция [http://www.cs.utexas.edu/users/plapack/C_Manual/#PLA_QR PLA_QR].
−
#* Перемножение матриц: функция [http://www.cs.utexas.edu/users/plapack/C_Manual/#PLA_Gemm PLA_Gemm]
+
#* Перемножение матриц: функция [http://www.cs.utexas.edu/users/plapack/C_Manual/#PLA_Gemm PLA_Gemm].
  
 
== Литература ==
 
== Литература ==

Текущая версия на 11:16, 10 декабря 2016

Symbol confirmed.svgЭта работа успешно выполнена
Преподавателю: в основное пространство, в подстраницу
Данное задание было проверено и зачтено.
Проверено Frolov и VadimVV.



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


Основные авторы описания:

Содержательно все пункты обсуждались авторами вместе, и их вклад был равным.


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

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

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

Содержание

1 Свойства и структура алгоритма

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

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

Суть базового 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-алгоритм со сдвигами. Это позволит уменьшить количество итераций алгоритма.

В дальнейшем, в данной статье под модифицированным алгоритмом будет пониматься алгоритм, использующий сдвиги и матрицу с формой Хессенберга. Под базовым алгоритмом будет пониматься алгоритм, не использующий данные приемы.

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

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

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

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

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

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

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]-го не изменяются). Таким образом, алгоритм состоит из [math]n-1[/math] шагов, на каждом из которых вычисляется очередная матрица отражения, после чего найденное отражение применяется к матрице, которая является результатом предыдущего шага.

Матрица отражения имеет вид [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] сходится к матрице с верхней треугольной формой[5].

Таким образом, на практике необходимо выполнять итерации до тех пор пока матрица [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-го порядка соответствуют парам комплексных сопряженных собственных значений[6].

[math]A_N= \begin{bmatrix} \blacksquare& \bullet& \bullet& \cdots& \cdots& \cdots& \cdots& \cdots& \bullet\\ 0& \blacksquare& \blacksquare& \bullet& \ddots& \ddots& \ddots& \ddots& \vdots\\ 0& \blacksquare& \blacksquare& \bullet& \bullet& \ddots& \ddots& \ddots& \vdots\\ \vdots& 0& 0& \blacksquare& \bullet& \bullet& \ddots& \ddots& \vdots\\ \vdots& \ddots& 0& 0& \blacksquare& \bullet& \bullet& \ddots& \vdots\\ \vdots& \ddots& \ddots& 0& 0& \blacksquare& \blacksquare& \ddots& \vdots\\ \vdots& \ddots& \ddots& \ddots& 0& \blacksquare& \blacksquare& \ddots& \bullet\\ \vdots& \ddots& \ddots& \ddots& \ddots& \ddots& \ddots& \ddots& \bullet\\ 0& \cdots& \cdots& \cdots& \cdots& \cdots& 0& 0& \blacksquare \end{bmatrix}[/math].

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

1.2.5 Приведение матрицы к форме Хессенберга

Основная статья: Метод Хаусхолдера (отражений) приведения матрицы к хессенберговой (почти треугольной) форме

Матрицей, имеющей форму Хессенберга, называется такая матрица, у которой все элементы, находящиеся ниже первой поддиагонали, равны нулю ([math]a_{ij}=0[/math] при [math]i\lt j+1[/math]). Пример такой матрицы приведен ниже:

[math]\begin{bmatrix} 1& 2& 3& 4\\ 2& 5& 6& 7\\ 0& 3& 8& 9\\ 0& 0& 4& 1 \end{bmatrix}[/math]

Любую матрицу [math]A[/math] можно привести к подобной ей матрице [math]A_H[/math], имеющей форму Хессенберга (в силу подобия данные матрицы будут иметь одинаковые собственные значения). Наличие нулевых элементов в данной матрице позволяет ускорить процесс QR-разложения, причем данное ускорение будет иметь место на каждой итерации алгоритма, т.к. матрица с формой Хессенберга инвариантна относительно QR-итерации. Ускорение можно получить за счет использования модифицированного метода Гивенса(вращений) QR-разложения, который из-за изначального наличия нулевых элементов будет состоять не из [math]\frac{n(n-1)}{2}[/math] шагов, а из [math]n-1[/math] шагов (будет происходить домножение только на те матрицы вращения [math]T_{ij}[/math], у которых [math]i=j+1[/math]).

Одним из способов приведения матрицы к форме Хессенберга является Метод Хаусхолдера (отражений)[7]. Суть метода Хаусхолдера заключается в последовательном приведении матрицы [math]A[/math] к форме Хессенберга при помощи домножения ее на так называемые матрицы отражения.

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

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

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

1.2.6 QR-алгоритм со сдвигами

QR-алгоритм со сдвигами позволяет сократить количество итераций, необходимых для сходимости[8]. Пусть у нас есть матрица [math]A_k[/math], тогда процесс перехода к матрице [math]A_{k+1}[/math] выглядит следующим образом:

  • На каждом шаге подбирается число [math]\nu_k[/math] и ищется следующее QR-разложение: [math]A_k-\nu_kE=Q_kR_k[/math].
  • Вычисляется матрица [math]A_{k+1}[/math]: [math]A_{k+1} = R_kQ_k+\nu_kE[/math].

При этом сохраняется свойство подобия матриц [math]A_k[/math] и [math]A_{k+1}[/math]:

[math]A_{k+1} = R_kQ_k+\nu_kE=Q_{k}^{T}Q_kR_kQ_k+\nu_kE=Q_{k}^{T}(A_k-\nu_kE)Q_k+\nu_kE=Q_{k}^{T}A_kQ_k-Q_{k}^{T}(\nu_kE)Q_k+\nu_kE=Q_{k}^{T}A_kQ_k-\nu_kE+\nu_kE=Q_{k}^{T}A_kQ_k[/math].

В качестве параметра [math]\nu_k[/math] может быть взят последний диагональный элемент матрицы [math]A_k[/math].

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

1.3.1 Базовый алгоритм

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

  1. Вычисление QR-разложения матрицы: [math]A_k=Q_kR_k[/math]. Существует несколько методов решения данной задачи:
    • Метод Хаусхолдера (отражений): вычислительное ядро данного алгоритма состоит из операций скалярного произведения, необходимых для вычисления матрицы отражения, и из операций скалярного произведения, необходимых для пересчета матрицы на каждом шаге.
    • Метод Гивенса (вращений): вычислительное ядро данного алгоритма состоит из численных арифметических операций, необходимых для вычисления параметров матрицы вращения, и из операций (эквивалентных перемножению комплексных чисел), необходимых для пересчета матрицы на каждом шаге.
  2. Перемножение двух плотных матриц: [math]A_{k+1}=R_kQ_k[/math].

1.3.2 Модифицированный алгоритм

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

  1. Приведение изначальной матрицы к форме Хессенберга: вычислительное ядро данного алгоритма состоит из операций скалярного произведения, необходимых для вычисления матрицы отражения, и из операций скалярного произведения, необходимых для пересчета матрицы на каждом шаге.
  2. Вычисление QR-разложения матрицы: [math]A_k-\nu_kE=Q_kR_k[/math] при помощи модифицированного метода Гивенса (вращений). Описание вычислительного ядра не отличается от приведенного в описании базового алгоритма, единственная разница заключается в значительно меньшем количестве матриц вращения.
  3. Перемножение двух плотных матриц: [math]A_{k+1}-\nu_kE=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].

В случае модифицированного QR-алгоритма появляется еще макрооперация приведения матрицы к форме Хессенберга.

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

1.5.1 Базовый алгоритм

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

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.5.2 Модифицированный алгоритм

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

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

curA = A
nextA = A
while ( not (triangular(nextA) & difference(curA,nextA,eps)) )
{
    curA = nextA
    nu = computeNu(curA)
    shiftA = curA - nu*E
    findQRdecomposition(shiftA,Q,R)
    nextA = R*Q + nu*E
}

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

1.6.1 Базовый алгоритм

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

  1. QR-разложение матрицы.
  2. Перемножение двух плотных матриц имеет сложность [math]n^3[/math].
  3. Проверка матрицы на квазитреугольную форму состоит из набора сравнений элементов с номерами [math]ij[/math] ([math]i\gt j+1[/math]) с нулем (таких элементов [math]\frac{(n-2)(n-1)}{2}[/math]). Поддиагональные элементы с номерами [math]ij[/math] ([math]i=j+1[/math]) должны быть проверены на соответствие необходимой квазитреугольной форме. Для этого достаточно для каждого такого элемента проверить следующее условие: [math]a_{i+1,i}==0 \vee a_{i,i-1}==0 \wedge a_{i+2,i+1}==0[/math] (либо поддиагональный элемент равен 0, либо, в противном случае, соседние поддиагональные элементы равны 0, чтобы текущий элемент соответствовал блоку 2-го порядка). Таких проверок поддиагональных элементов будет [math]n-1[/math]. После данных проверок следует набор логических операций "И" между результатами всех сравнений (таких операций будет [math]\frac{n(n-1)}{2}-1[/math]).
  4. Сравнение новой матрицы с предыдущей состоит из операций вычитания и сравнения для каждой пары ненулевых соответствующих элементов (такие пары имеют номера элементов [math]ij[/math] ([math]i \geq j+1[/math]), количество таких пар равно [math]\frac{n(n+1)}{2}+(n-1)[/math]), а также из набора логических операций "И" между результатами сравнений (таких операций будет [math]\frac{n(n+1)}{2}+(n-1)-1[/math]).

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

1.6.2 Модифицированный алгоритм

Отличия от базового алгоритма заключаются в следующем:

  1. Метод Хаусхолдера (отражений) имеет последовательную сложность [math]O(n^3)[/math].
  2. QR-разложение матрицы, имеющей форму Хессенберга будет иметь сложность порядка [math]O(n^2)[/math][8].

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

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

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

1.7.1 Базовый алгоритм

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

  • QR - операция QR-разложения матрицы.
  • *M - операция перемножения матриц.
  • Triang - операция проверки матрицы на квазитреугольную форму.
  • Dif - операция проверки различия элементов двух матриц не более чем на некоторое заданное число.
  • Iteri - итерация алгоритма
Рисунок 1. Информационный граф базового QR-алгоритма
Рисунок 2. Информационный граф итерации базового QR-алгоритма

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

Рисунок 3. Информационный граф операции Triang, которая проверяет соответствие квазитреугольной форме

Вершины V соответствуют операции сравнения с 0. Вершины F соответствуют проверке поддиагональных элементов на соответствие квазитреугольной форме, которая была описана в предыдущем разделе. Вершины & соответствуют логической операции "И".

Рисунок 4. Информационный граф операции Dif, которая проверяет сходимость ненулевых элементов

Вершины -V соответствуют операции вычитания и сравнения с 0. Вершины & соответствуют логической операции "И".

1.7.2 Модифицированный алгоритм

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

  • Shift - операция сдвига
  • H - операция приведения матрицы к форме Хессенберга
Рисунок 5. Информационный граф модифицированного QR-алгоритма
Рисунок 6. Информационный граф итерации модифицированного QR-алгоритма

1.8 Ресурс параллелизма алгоритма

1.8.1 Базовый алгоритм

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

  1. QR-разложение матрицы (описание ресурсов параллелизма для алгоритмов QR-разложения можно найти в статьях, посвященных этим алгоритмам).
  2. Перемножение двух плотных матриц имеет высоту ярусно-параллельной формы [math]n[/math] и ширину ярусно-параллельной формы [math]n^2[/math].
  3. Проверка матрицы на выходе из итерации.
    • Проверка матрицы на квазитреугольную форму состоит из одного яруса сравнений для каждого элемента и последующих ярусов, вычисляющих итоговый результат при помощи логической операции "И". Для вычисления итогового результата можно использовать метод сдваивания, который дает высоту ярусно-параллельной формы порядка логарифма от количества элементов, над которыми совершается операция. Таким образов высота ярусно-параллельной формы будет равна [math]O(log_2n)[/math]. Ширина ярусно-параллельной формы достигается на ярусе сравнений для каждого элемента и равна [math]O(n^2)[/math].
    • Сравнение новой матрицы с предыдущей состоит из одного яруса сравнений для каждого элемента и последующих ярусов, вычисляющих итоговый результат при помощи логической операции "И". Высота ярусно-параллельной формы будет равна [math]O(log_2n)[/math]. Ширина ярусно-параллельной формы достигается на ярусе сравнений для каждого элемента и равна [math]O(n^2)[/math].

Таким образом, основной вклад в высоту ярусно-параллельной формы одной итерации вносит операция QR-разложения матрицы и она будет равна [math]O(n)[/math], если использовать метод Гивенса. Ширина ярусно-параллельной формы будет равна [math]O(n^2)[/math]. Если алгоритм остановился на итерации с номером [math]N[/math], то параллельные характеристики для всего алгоритма будут равны [math]N*O(n)[/math] для высоты и [math]O(n^2)[/math] для ширины ярусно-параллельной формы.

1.8.2 Модифицированный алгоритм

Отличия от базового алгоритма заключаются в следующем:

  1. Приведение матрицы к форме Хессенберга имеет высоту ярусно-параллельной формы [math]O(n^2)[/math] и ширину ярусно-параллельной формы [math]O(n^2)[/math].
  2. QR-разложение матрицы модифицированным методом Гивенса (вращений) имеет высоту ярусно-параллельной формы [math]O(n)[/math] и ширину ярусно-параллельной формы [math]O(n)[/math].
  3. Операция сдвига (вычитание из диагональных элементов матрицы одного числа) имеет высоту ярусно-параллельной формы [math]O(1)[/math] и ширину ярусно-параллельной формы [math]O(n)[/math].

Таким образом, итоговые характеристики всего модифицированного алгоритма будут следующие: высота ярусно-параллельной формы [math]O(n^2)+N*O(n)[/math] и ширина ярусно-параллельной формы [math]O(n^2)[/math].

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] собственных значений при этом некоторые из них могут быть комплексными).

1.10 Свойства алгоритма

  • Cоотношение последовательной и параллельной сложности алгоритма квадратично, что дает довольно большой выигрыш.
  • Вычислительная мощность, которая показывает, сколько операций приходится на единицу переданных данных, равна [math]\frac{N*O(n^3)}{n^2+n}=N*O(n)[/math], а значит перемещение данных для их обработки не будет составлять большой проблемы.
  • Алгоритм является недетерминированным, т.к. заранее неизвестно сколько итераций необходимо совершить до момента сходимости.
  • Скорость сходимости алгоритма зависит от собственных значений. Чем ближе по модулю соседние собственные значения, тем меньше скорость сходимости. Этот недостаток призван решить QR-алгоритм со сдвигами.

2 Программная реализация алгоритма

2.1 Особенности реализации последовательного алгоритма

2.2 Локальность данных и вычислений

2.3 Возможные способы и особенности параллельной реализации алгоритма

2.4 Масштабируемость алгоритма и его реализации

Масштабируемость алгоритма (рис.7) исследовалась на основе реализации алгоритма в библиотеке ScaLAPACK v2.0.0: функция pdhseqr. Время работы было замерено для различных размеров матрицы, а также различного числа процессоров.

Рисунок 7. Масштабируемость QR-алгоритма

При увеличении числа процессоров время работы программы увеличивается при малом размере матрицы, распараллеливание алгоритма неэффективно. Однако при больших размерах матрицы (более [math] 10^3\times10^3[/math]) время работы параллельной реализации значительно меньше последовательной.

Код программы

Все вычисления были проведены на суперкомпьютере Regatta. Библиотека ScaLAPACK v2.0.0 была установлена вручную. Компиляция проводилась с использованием компилятора mpicxx с опциями -lscalapack и -llapack с указанием полного пути до используемых библиотек.

2.5 Динамические характеристики и эффективность реализации алгоритма

2.6 Выводы для классов архитектур

2.7 Существующие реализации алгоритма

2.7.1 Последовательные реализации

  1. LAPACK (Linear Algebra PACKage) — библиотека с открытым исходным кодом, содержащая методы для решения основных задач линейной алгебры. Написана на языке Fortran с использованием библиотеки BLAS и является развитием пакета LINPACK. Находится в открытом доступе в соответствии с модифицированной лицензией BSD, в том числе и для коммерческого использования.
    • Полный алгоритм нахождения собственных значений: функция dhseqr.
    • QR-разложение: функция dgeqrf.
    • Перемножение матриц: функция dgemm.
  2. ALGLIB - это кросс-платформенная библиотека численного анализа, поддерживающая несколько языков программирования (C++, C#, Pascal, VBA) и несколько операционных систем (Windows, Linux, Solaris). ALGLIB является свободным программным обеспечением, которое использует двойное лицензирование: оно может быть использовано в соответствии с лицензией GPL (версии 2+), а для использования в коммерческих целях необходимо купить отдельную лицензию.
    • Полный алгоритм нахождения собственных значений: подпрограмма RMatrixEVD.
    • QR-разложение: подпрограмма rmatrixqr реализует QR-разложение для вещественных матриц, cmatrixqr – для комплексных матриц.
    • Перемножение матриц: для перемножения матриц ALGLIB использует соответствующую реализацию библиотеки BLAS.
  3. Eigen – библиотека шаблонов для линейной алгебры, написанная на языке C++. Eigen – свободно распространяемое программное обеспечение. Начиная с версии 3.1.1, оно лицензируется MPL2, на ранние версии распространяется действие лицензии LGPL3+.

2.7.2 Параллельные реализации

  1. ScaLAPACK (Scalable Linear Algebra PACKage) — библиотека с открытым исходным кодом, включающая в себя подмножество процедур LAPACK, переработанных для использования на MPP-компьютерах. ScaLAPACK разработана с использованием PBLAS и BLACS, и предназначена для вычислений на любом компьютере или кластере поддерживающим MPI или PVM. Библиотека в настоящее время написана на языке Fortran. Находится в открытом доступе в соответствии с модифицированной лицензией BSD, в том числе и для коммерческого использования.
    • Полный алгоритм нахождения собственных значений: функция pdhseqr.
    • QR-разложение: функция pdgeqrf.
    • Перемножение матриц: функция pdgemm.
  2. PLAPACK (Parallel Linear Algebra Package) — пакет функций LAPACK для параллельного решения задач линейной алгебры. Пакет функций PLAPACK является альтернативой библиотеке ScaLAPACK. Для осуществления межпроцессорных коммуникаций в PLAPACK использован интерфейс передачи сообщений MPI. При передаче сообщений в PLAPACK в основном используются коллективные операции, такие, как обобщенная передача данных от одного процесса всем процессам (MPI_Scatter), обобщенная передача данных от всех процессов одному процессу (MPI_Gather), широковещательная рассылка (MPI_Bcast) и другие. PLAPACK включает интерфейсы для языков Fortran и С.
    • QR-разложение: функция PLA_QR.
    • Перемножение матриц: функция PLA_Gemm.

3 Литература

  1. ↑ Ильин В.А., Ким Г.Д. "Линейная алгебра и аналитическая геометрия".
  2. ↑ Wikipedia: QR-algorithm
  3. ↑ Метод Хаусхолдера (отражений) QR-разложения квадратной матрицы
  4. ↑ Метод Гивенса (вращений) QR-разложения квадратной матрицы
  5. ↑ Тыртышников Е.Е. "Методы численного анализа" — М., Академия, 2007. - 320 c.
  6. ↑ R. Granat, Bo Kagstrom, D. Kressner "LAPACK Working Note #216: A novel parallel QR algorithm for hybrid distributed memory HPC systems".
  7. ↑ Метод Хаусхолдера (отражений) приведения матрицы к хессенберговой (почти треугольной) форме
  8. ↑ 8,0 8,1 Бахвалов Н.С., Жидков Н.П., Кобельков. Г.М. "Численный методы"— 6-е изд. — М. : БИНОМ. Лаборатория знаний, 2008. — 636 с.
Источник — https://algowiki-project.org/w/ru/index.php?title=Участник:Смирнова_Александра/Нахождение_собственных_чисел_квадратной_матрицы_методом_QR_разложения_(3)&oldid=20510