Участник:Смирнова Александра/Нахождение собственных чисел квадратной матрицы методом QR разложения (3)
![]() | Эта работа успешно выполнена Преподавателю: в основное пространство, в подстраницу Данное задание было проверено и зачтено. Проверено Frolov и VadimVV. |
Нахождение собственных чисел квадратной матрицы методом QR разложения | |
Последовательный алгоритм | |
Последовательная сложность | O(n^3)+N*O(n^2) |
Объём входных данных | n^2 |
Объём выходных данных | n |
Параллельный алгоритм | |
Высота ярусно-параллельной формы | O(n^2)+N*O(n) |
Ширина ярусно-параллельной формы | O(n^2) |
Основные авторы описания:
- Смирнова А.С. - описание теоретической части
- Киямова А. - описание программистской части
Содержательно все пункты обсуждались авторами вместе, и их вклад был равным.
Задача нахождения собственных значений и собственных векторов для матрицы A заключается в поиске таких чисел \lambda, которые удовлетворяют уравнению:
Ax=\lambda x, при этом, числа \lambda называются собственными значениями, а вектора x - собственными векторами[1].
Данная задача является одной из важнейших задач линейной алгебры. Собственные вектора и собственные значения применяются в различных областях науки: в аналитической геометрии, при решении систем интегральных уравнений, в математической физике. Однако не существует простых алгоритмов прямого вычисления собственных значений для матриц в общем случае, поэтому данная задача на практике решается численными методами. Одним из таких методов является QR-алгоритм.
Содержание
- 1 Свойства и структура алгоритма
- 1.1 Общее описание алгоритма
- 1.2 Математическое описание алгоритма
- 1.3 Вычислительное ядро алгоритма
- 1.4 Макроструктура алгоритма
- 1.5 Схема реализации последовательного алгоритма
- 1.6 Последовательная сложность алгоритма
- 1.7 Информационный граф
- 1.8 Ресурс параллелизма алгоритма
- 1.9 Входные и выходные данные алгоритма
- 1.10 Свойства алгоритма
- 2 Программная реализация алгоритма
- 2.1 Особенности реализации последовательного алгоритма
- 2.2 Локальность данных и вычислений
- 2.3 Возможные способы и особенности параллельной реализации алгоритма
- 2.4 Масштабируемость алгоритма и его реализации
- 2.5 Динамические характеристики и эффективность реализации алгоритма
- 2.6 Выводы для классов архитектур
- 2.7 Существующие реализации алгоритма
- 3 Литература
1 Свойства и структура алгоритма
1.1 Общее описание алгоритма
QR-алгоритм — это численный метод в линейной алгебре, предназначенный для решения полной проблемы собственных значений, то есть отыскания всех собственных чисел матрицы. При этом алгоритм позволяет найти и собственные вектора матрицы. Он был разработан в конце 1950-х годов независимо В. Н. Кублановской(Россия) и Дж. Фрэнсисом(Англия). Открытию QR-алгоритма предшествовал LR-алгоритм, который использовал LU-разложение вместо QR-разложения. В настоящее время LR-алгоритм используется очень редко ввиду своей меньшей эффективности, однако он был важным шагом на пути к открытию QR-алгоритма[2].
Суть базового QR-алгоритма заключается в итерационном приведении матрицы A к некоторой подобной ей матрице A_N при помощи QR-разложения. Матрица A_N является правой верхней треугольной матрицей, а значит ее диагональ содержит собственные значения. В силу подобия матриц A и A_N их наборы собственных значений совпадают. Таким образом задача поиска собственных значений матрицы A сводится к задаче выведения матрицы A_N и поиска собственных значений для нее, что является тривиальной задачей.
Однако базовый QR-алгоритм может обладать очень низкой скоростью сходимости, поэтому существует несколько способов ускорить его:
- Перед итерациями привести матрицу A к подобной ей матрице A_H, которая будет иметь форму Хессенберга. Данный шаг позволит ускорить процесс QR-разложения.
- Использовать QR-алгоритм со сдвигами. Это позволит уменьшить количество итераций алгоритма.
В дальнейшем, в данной статье под модифицированным алгоритмом будет пониматься алгоритм, использующий сдвиги и матрицу с формой Хессенберга. Под базовым алгоритмом будет пониматься алгоритм, не использующий данные приемы.
1.2 Математическое описание алгоритма
1.2.1 QR-разложение матрицы
Основой алгоритма является тот факт, что любую вещественную матрицу можно разложить на произведение двух матриц следующего вида:
- A=QR, где Q - ортогональная матрица (Q^T=Q^{-1}), R - верхняя треугольная матрица.
Данное разложение называется QR-разложением.
Есть несколько алгоритмов вычисления QR-разложения матрицы[3] [4]. Кратко опишем их в данной статье.
1.2.1.1 Метод Хаусхолдера (отражений) QR-разложения квадратной матрицы
- Основная статья: Метод Хаусхолдера (отражений) QR-разложения квадратной матрицы
Суть метода Хаусхолдера заключается в последовательном приведении матрицы A к верхней треугольной форме при помощи домножения ее на так называемые матрицы отражения. Получившаяся треугольная матрица будет искомой матрицей R, а матрица Q будет равна произведению сопряженных матриц отражения.
На i-ом шаге задача i-ой матрицы отражения заключается в обнулении всех поддиагональных элементов i-го столбца матрицы A (при этом столбцы левее i-го не изменяются). Таким образом, алгоритм состоит из n-1 шагов, на каждом из которых вычисляется очередная матрица отражения, после чего найденное отражение применяется к матрице, которая является результатом предыдущего шага.
Матрица отражения имеет вид E-\frac{1}{\gamma }vv^*, где вектор v вычисляется из текущего i-го столбца матрицы при помощи использования операции скалярного произведения векторов. Данное представление матриц отражения позволяет хранить их в виде одного вектора и сводит домножение матрицы отражения на текущую матрицу к арифметическим операциям над векторами (скалярное произведение и сложение векторов).
1.2.1.2 Метод Гивенса (вращений) QR-разложения квадратной матрицы
- Основная статья: Метод Гивенса (вращений) QR-разложения квадратной матрицы
Суть метода Гивенса заключается в последовательном приведении матрицы A к верхней треугольной форме при помощи домножения ее на так называемые матрицы вращения. Получившаяся треугольная матрица будет искомой матрицей R, а матрица Q будет равна произведению сопряженных матриц вращения.
На каждом шаге задачей матрицы вращения является обнуление одного поддиагонального элемента. Вначале обнуляются все поддиагональные элементы 1-го столбца, затем 2-го и так далее до (n-1)-го. Таким образом, алгоритм состоит из \frac{n(n-1)}{2} шагов на каждом из которых вычисляется очередная матрица вращения, после чего она применяется к матрице, которая является результатом предыдущего шага.
Матрицы вращения T_{ij} по ствоей структуре очень близки к единичным матрицам, за исключением четырех элементов: элементы с номерами ii и jj содержат некоторое число-параметр c, элементы с номерами ij и ji содержат числа-параметры -s и s соответственно. Вычисление параметров c и s происходит на каждом шаге в зависимости от текущей матрицы и состоит из простых численных арифметических операций. Умножение матрицы вращения на текущую матрицу может быть представлено как последовательное изменение элементов с номерами ik и kk для всех столбцов k находящихся правее столбца i. Каждое такое изменение по своей структуре эквивалентно операции перемножения двух комплексных чисел.
1.2.2 QR-алгоритм нахождения собственных чисел
Пусть матрица A - матрица, для которой мы хотим найти собственные значения. Тогда итерационный процесс строится следующим образом:
- A_0=A.
- Пусть имеется матрица A_k, тогда матрица A_{k+1} строится следующим образом:
- Строится QR-разложение: A_k=Q_kR_k.
- Вычисляется A_{k+1}=R_kQ_k.
Заметим, что 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.
Таким образом матрицы A_{k+1} и A_{k} подобны для \forall k, а значит, в силу транзитивности подобия, все матрицы A_{k} подобны матрице A и имеют одинаковый набор собственных значений.
1.2.3 Сходимость алгоритма
Предположим, что для \forall m выполнены следующие условия:
- 1. A=X\Lambda X^{-1}, где \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} .
- 2. \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 , где \{\lambda_1,...,\lambda_m\} = \lambda(\Lambda_1), \{\lambda_{m+1},...,\lambda_{m+r}\} = \lambda(\Lambda_2) .
- 3. Ведущая подматрица порядка m в X^{-1} невырождена.
Тогда при k \rightarrow \infty последовательность матриц A_k сходится к матрице с верхней треугольной формой[5].
Таким образом, на практике необходимо выполнять итерации до тех пор пока матрица A_k не станет треугольной (также можно продолжать выполнять их пока искомая матрица не будет найдена с некоторой заранее заданной точностью \varepsilon). Если итерации закончились на шаге N, то числа на диагонали матрицы A_N будут считаться собственными значениями матрицы A.
1.2.4 Вещественный вариант QR-алгоритма
Если вещественная матрица A имеет различные вещественные собственные значения, то, как было описано ранее, QR-алгоритм сходится к матрице с верхней треугольной формой, на диагонали которой находятся собственные значения. Однако вещественная матрица может иметь комплексные собственные значения. В данном случае алгоритм будет сходиться не к верхней треугольной матрице, а к блочной верхней треугольной матрице, которая на диагонали содержит блоки 1-го и 2-го порядка. Блоки 1-го порядка содержат различные вещественные собственные значения, блоки 2-го порядка соответствуют парам комплексных сопряженных собственных значений[6].
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}.
В дальнейшем, в данной статье, матрицы, имеющие вышеописанную форму, будут называться квазитреугольными.
1.2.5 Приведение матрицы к форме Хессенберга
- Основная статья: Метод Хаусхолдера (отражений) приведения матрицы к хессенберговой (почти треугольной) форме
Матрицей, имеющей форму Хессенберга, называется такая матрица, у которой все элементы, находящиеся ниже первой поддиагонали, равны нулю (a_{ij}=0 при i\lt j+1). Пример такой матрицы приведен ниже:
\begin{bmatrix} 1& 2& 3& 4\\ 2& 5& 6& 7\\ 0& 3& 8& 9\\ 0& 0& 4& 1 \end{bmatrix}
Любую матрицу A можно привести к подобной ей матрице A_H, имеющей форму Хессенберга (в силу подобия данные матрицы будут иметь одинаковые собственные значения). Наличие нулевых элементов в данной матрице позволяет ускорить процесс QR-разложения, причем данное ускорение будет иметь место на каждой итерации алгоритма, т.к. матрица с формой Хессенберга инвариантна относительно QR-итерации. Ускорение можно получить за счет использования модифицированного метода Гивенса(вращений) QR-разложения, который из-за изначального наличия нулевых элементов будет состоять не из \frac{n(n-1)}{2} шагов, а из n-1 шагов (будет происходить домножение только на те матрицы вращения T_{ij}, у которых i=j+1).
Одним из способов приведения матрицы к форме Хессенберга является Метод Хаусхолдера (отражений)[7]. Суть метода Хаусхолдера заключается в последовательном приведении матрицы A к форме Хессенберга при помощи домножения ее на так называемые матрицы отражения.
На i-ом шаге задача i-ой матрицы отражения заключается в обнулении поддиагональных элементов i-го столбца матрицы A (при этом столбцы левее i-го не изменяются). Таким образом, алгоритм состоит из n-2 шагов, на каждом из которых вычисляется очередная матрица отражения, после чего найденное отражение применяется к матрице, которая является результатом предыдущего шага.
Матрица отражения имеет вид E-\frac{1}{\gamma }ww^*, где вектор w вычисляется из текущего i-го столбца матрицы при помощи использования операции скалярного произведения векторов. Данное представление матриц отражения позволяет хранить их в виде одного вектора и сводит домножение матрицы отражения на текущую матрицу к арифметическим операциям над векторами (скалярное произведение и сложение векторов).
В основной статье "Метод Хаусхолдера (отражений) приведения матрицы к хессенберговой (почти треугольной) форме" данный алгоритм описан более подробно с той лишь разницей, что для QR-алгоритма нет необходимости обнулять наддиагональные элементы, как это сделано в вышеупомянутой статье.
1.2.6 QR-алгоритм со сдвигами
QR-алгоритм со сдвигами позволяет сократить количество итераций, необходимых для сходимости[8]. Пусть у нас есть матрица A_k, тогда процесс перехода к матрице A_{k+1} выглядит следующим образом:
- На каждом шаге подбирается число \nu_k и ищется следующее QR-разложение: A_k-\nu_kE=Q_kR_k.
- Вычисляется матрица A_{k+1}: A_{k+1} = R_kQ_k+\nu_kE.
При этом сохраняется свойство подобия матриц A_k и A_{k+1}:
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.
В качестве параметра \nu_k может быть взят последний диагональный элемент матрицы A_k.
1.3 Вычислительное ядро алгоритма
1.3.1 Базовый алгоритм
QR-алгоритм обладает двумя вычислительными ядрами, которые повторяются на каждой итерации:
- Вычисление QR-разложения матрицы: A_k=Q_kR_k. Существует несколько методов решения данной задачи:
- Метод Хаусхолдера (отражений): вычислительное ядро данного алгоритма состоит из операций скалярного произведения, необходимых для вычисления матрицы отражения, и из операций скалярного произведения, необходимых для пересчета матрицы на каждом шаге.
- Метод Гивенса (вращений): вычислительное ядро данного алгоритма состоит из численных арифметических операций, необходимых для вычисления параметров матрицы вращения, и из операций (эквивалентных перемножению комплексных чисел), необходимых для пересчета матрицы на каждом шаге.
- Перемножение двух плотных матриц: A_{k+1}=R_kQ_k.
1.3.2 Модифицированный алгоритм
Модифицированный QR-алгоритм обладает тремя вычислительными ядрами, первое вычисляется единожды, второе и третье повторяются на каждой итерации:
- Приведение изначальной матрицы к форме Хессенберга: вычислительное ядро данного алгоритма состоит из операций скалярного произведения, необходимых для вычисления матрицы отражения, и из операций скалярного произведения, необходимых для пересчета матрицы на каждом шаге.
- Вычисление QR-разложения матрицы: A_k-\nu_kE=Q_kR_k при помощи модифицированного метода Гивенса (вращений). Описание вычислительного ядра не отличается от приведенного в описании базового алгоритма, единственная разница заключается в значительно меньшем количестве матриц вращения.
- Перемножение двух плотных матриц: A_{k+1}-\nu_kE=R_kQ_k.
1.4 Макроструктура алгоритма
Как уже было описано ранее, базовый QR-алгоритм содержит в себе две макрооперации:
- Вычисление QR-разложения матрицы: A_k=Q_kR_k.
- Перемножение двух плотных матриц: A_{k+1}=R_kQ_k.
В случае модифицированного 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-разложения и перемножения матриц представлен в статьях, посвященных данным алгоритмам).
- QR-разложение матрицы.
- Метод Хаусхолдера (отражений) имеет сложность \frac{4}{3}n^3.
- Метод Гивенса (вращений) имеет сложность 2n^3.
- Перемножение двух плотных матриц имеет сложность n^3.
- Проверка матрицы на квазитреугольную форму состоит из набора сравнений элементов с номерами ij (i\gt j+1) с нулем (таких элементов \frac{(n-2)(n-1)}{2}). Поддиагональные элементы с номерами ij (i=j+1) должны быть проверены на соответствие необходимой квазитреугольной форме. Для этого достаточно для каждого такого элемента проверить следующее условие: a_{i+1,i}==0 \vee a_{i,i-1}==0 \wedge a_{i+2,i+1}==0 (либо поддиагональный элемент равен 0, либо, в противном случае, соседние поддиагональные элементы равны 0, чтобы текущий элемент соответствовал блоку 2-го порядка). Таких проверок поддиагональных элементов будет n-1. После данных проверок следует набор логических операций "И" между результатами всех сравнений (таких операций будет \frac{n(n-1)}{2}-1).
- Сравнение новой матрицы с предыдущей состоит из операций вычитания и сравнения для каждой пары ненулевых соответствующих элементов (такие пары имеют номера элементов ij (i \geq j+1), количество таких пар равно \frac{n(n+1)}{2}+(n-1)), а также из набора логических операций "И" между результатами сравнений (таких операций будет \frac{n(n+1)}{2}+(n-1)-1).
Итого, в сумме получаем O(n^3) - сложность алгоритма на каждой итерации. Если алгоритм остановился на итерации с номером N, то общая сложность алгоритма будет равна N*O(n^3).
1.6.2 Модифицированный алгоритм
Отличия от базового алгоритма заключаются в следующем:
- Метод Хаусхолдера (отражений) имеет последовательную сложность O(n^3).
- QR-разложение матрицы, имеющей форму Хессенберга будет иметь сложность порядка O(n^2)[8].
Стоит заметить, что вычисление сдвигов по сравнению с другими операциями не обладает высокой сложностью, однако позволяет значительно сократить количество необходимых итераций.
Итого, в сумме получаем O(n^2) - сложность алгоритма на каждой итерации. Сложность всего алгоритма, с учетом приведения матрицы к форме Хессенберга, будет равна O(n^3)+N*O(n^2)
1.7 Информационный граф
1.7.1 Базовый алгоритм
На рисунках 1 и 2 изображен информационный граф QR-алгоритма. Вершины данного графа обозначают следующее:
- QR - операция QR-разложения матрицы.
- *M - операция перемножения матриц.
- Triang - операция проверки матрицы на квазитреугольную форму.
- Dif - операция проверки различия элементов двух матриц не более чем на некоторое заданное число.
- Iteri - итерация алгоритма
Подробные графы операций QR-разложения (Метод Хаусхолдера (отражений), Метод Гивенса (вращений)) и перемножения матриц можно найти в статьях, посвященных этим алгоритмам. Далее рассмотрим информационные графы операций Triang (рис.3) и Dif (рис.4) на примере матрицы размера 5 \times 5. Графы для других матриц выглядят аналогичным образом.
Вершины V соответствуют операции сравнения с 0. Вершины F соответствуют проверке поддиагональных элементов на соответствие квазитреугольной форме, которая была описана в предыдущем разделе. Вершины & соответствуют логической операции "И".
Вершины -V соответствуют операции вычитания и сравнения с 0. Вершины & соответствуют логической операции "И".
1.7.2 Модифицированный алгоритм
На рисунках 5 и 6 изображен информационный граф модифицированного QR-алгоритма. Дополнительные вершины для этого графа обозначают следующее:
- Shift - операция сдвига
- H - операция приведения матрицы к форме Хессенберга
1.8 Ресурс параллелизма алгоритма
1.8.1 Базовый алгоритм
На макроуровне (который можно увидеть на информационном графе QR-алгоритма) алгоритм не обладает ресурсами параллелизма. Все макрооперации внутри итерации, а также сами итерации выполняются последовательно (за исключением операций Triag и Dif, которые могут выполняться параллельно). Основной ресурс параллелизма заложен отдельно в каждой из макроопераций. На каждой итерации алгоритм имеет следующие параллельные характеристики:
- QR-разложение матрицы (описание ресурсов параллелизма для алгоритмов QR-разложения можно найти в статьях, посвященных этим алгоритмам).
- Метод Хаусхолдера (отражений) имеет высоту ярусно-параллельной формы O(n^2) и ширину ярусно-параллельной формы O(n).
- Метод Гивенса (вращений) имеет высоту ярусно-параллельной формы 11n-16 и ширину ярусно-параллельной формы O(n^2).
- Перемножение двух плотных матриц имеет высоту ярусно-параллельной формы n и ширину ярусно-параллельной формы n^2.
- Проверка матрицы на выходе из итерации.
- Проверка матрицы на квазитреугольную форму состоит из одного яруса сравнений для каждого элемента и последующих ярусов, вычисляющих итоговый результат при помощи логической операции "И". Для вычисления итогового результата можно использовать метод сдваивания, который дает высоту ярусно-параллельной формы порядка логарифма от количества элементов, над которыми совершается операция. Таким образов высота ярусно-параллельной формы будет равна O(log_2n). Ширина ярусно-параллельной формы достигается на ярусе сравнений для каждого элемента и равна O(n^2).
- Сравнение новой матрицы с предыдущей состоит из одного яруса сравнений для каждого элемента и последующих ярусов, вычисляющих итоговый результат при помощи логической операции "И". Высота ярусно-параллельной формы будет равна O(log_2n). Ширина ярусно-параллельной формы достигается на ярусе сравнений для каждого элемента и равна O(n^2).
Таким образом, основной вклад в высоту ярусно-параллельной формы одной итерации вносит операция QR-разложения матрицы и она будет равна O(n), если использовать метод Гивенса. Ширина ярусно-параллельной формы будет равна O(n^2). Если алгоритм остановился на итерации с номером N, то параллельные характеристики для всего алгоритма будут равны N*O(n) для высоты и O(n^2) для ширины ярусно-параллельной формы.
1.8.2 Модифицированный алгоритм
Отличия от базового алгоритма заключаются в следующем:
- Приведение матрицы к форме Хессенберга имеет высоту ярусно-параллельной формы O(n^2) и ширину ярусно-параллельной формы O(n^2).
- QR-разложение матрицы модифицированным методом Гивенса (вращений) имеет высоту ярусно-параллельной формы O(n) и ширину ярусно-параллельной формы O(n).
- Операция сдвига (вычитание из диагональных элементов матрицы одного числа) имеет высоту ярусно-параллельной формы O(1) и ширину ярусно-параллельной формы O(n).
Таким образом, итоговые характеристики всего модифицированного алгоритма будут следующие: высота ярусно-параллельной формы O(n^2)+N*O(n) и ширина ярусно-параллельной формы O(n^2).
1.9 Входные и выходные данные алгоритма
Входные данные:
- квадратная вещественная плотная матрица A: A \in \R^{n \times n}.
Объем входных данных:
- n^2 (необходимо хранить все элементы матрицы).
Выходные данные:
- собственные значения матрицы A.
Объем выходных данных:
- n (квадратная матрица размера n \times n имеет ровно n собственных значений при этом некоторые из них могут быть комплексными).
1.10 Свойства алгоритма
- Cоотношение последовательной и параллельной сложности алгоритма квадратично, что дает довольно большой выигрыш.
- Вычислительная мощность, которая показывает, сколько операций приходится на единицу переданных данных, равна \frac{N*O(n^3)}{n^2+n}=N*O(n), а значит перемещение данных для их обработки не будет составлять большой проблемы.
- Алгоритм является недетерминированным, т.к. заранее неизвестно сколько итераций необходимо совершить до момента сходимости.
- Скорость сходимости алгоритма зависит от собственных значений. Чем ближе по модулю соседние собственные значения, тем меньше скорость сходимости. Этот недостаток призван решить QR-алгоритм со сдвигами.
2 Программная реализация алгоритма
2.1 Особенности реализации последовательного алгоритма
2.2 Локальность данных и вычислений
2.3 Возможные способы и особенности параллельной реализации алгоритма
2.4 Масштабируемость алгоритма и его реализации
Масштабируемость алгоритма (рис.7) исследовалась на основе реализации алгоритма в библиотеке ScaLAPACK v2.0.0: функция pdhseqr. Время работы было замерено для различных размеров матрицы, а также различного числа процессоров.
При увеличении числа процессоров время работы программы увеличивается при малом размере матрицы, распараллеливание алгоритма неэффективно. Однако при больших размерах матрицы (более 10^3\times10^3) время работы параллельной реализации значительно меньше последовательной.
Все вычисления были проведены на суперкомпьютере Regatta. Библиотека ScaLAPACK v2.0.0 была установлена вручную. Компиляция проводилась с использованием компилятора mpicxx с опциями -lscalapack и -llapack с указанием полного пути до используемых библиотек.
2.5 Динамические характеристики и эффективность реализации алгоритма
2.6 Выводы для классов архитектур
2.7 Существующие реализации алгоритма
2.7.1 Последовательные реализации
- LAPACK (Linear Algebra PACKage) — библиотека с открытым исходным кодом, содержащая методы для решения основных задач линейной алгебры. Написана на языке Fortran с использованием библиотеки BLAS и является развитием пакета LINPACK. Находится в открытом доступе в соответствии с модифицированной лицензией BSD, в том числе и для коммерческого использования.
- ALGLIB - это кросс-платформенная библиотека численного анализа, поддерживающая несколько языков программирования (C++, C#, Pascal, VBA) и несколько операционных систем (Windows, Linux, Solaris). ALGLIB является свободным программным обеспечением, которое использует двойное лицензирование: оно может быть использовано в соответствии с лицензией GPL (версии 2+), а для использования в коммерческих целях необходимо купить отдельную лицензию.
- Полный алгоритм нахождения собственных значений: подпрограмма RMatrixEVD.
- QR-разложение: подпрограмма rmatrixqr реализует QR-разложение для вещественных матриц, cmatrixqr – для комплексных матриц.
- Перемножение матриц: для перемножения матриц ALGLIB использует соответствующую реализацию библиотеки BLAS.
- Eigen – библиотека шаблонов для линейной алгебры, написанная на языке C++. Eigen – свободно распространяемое программное обеспечение. Начиная с версии 3.1.1, оно лицензируется MPL2, на ранние версии распространяется действие лицензии LGPL3+.
- Полный алгоритм нахождения собственных значений: модуль Eigenvalues module.
- QR-разложение: модуль QR module.
- Перемножение матриц: операторы * и *=.
2.7.2 Параллельные реализации
- ScaLAPACK (Scalable Linear Algebra PACKage) — библиотека с открытым исходным кодом, включающая в себя подмножество процедур LAPACK, переработанных для использования на MPP-компьютерах. ScaLAPACK разработана с использованием PBLAS и BLACS, и предназначена для вычислений на любом компьютере или кластере поддерживающим MPI или PVM. Библиотека в настоящее время написана на языке Fortran. Находится в открытом доступе в соответствии с модифицированной лицензией BSD, в том числе и для коммерческого использования.
- PLAPACK (Parallel Linear Algebra Package) — пакет функций LAPACK для параллельного решения задач линейной алгебры. Пакет функций PLAPACK является альтернативой библиотеке ScaLAPACK. Для осуществления межпроцессорных коммуникаций в PLAPACK использован интерфейс передачи сообщений MPI. При передаче сообщений в PLAPACK в основном используются коллективные операции, такие, как обобщенная передача данных от одного процесса всем процессам (MPI_Scatter), обобщенная передача данных от всех процессов одному процессу (MPI_Gather), широковещательная рассылка (MPI_Bcast) и другие. PLAPACK включает интерфейсы для языков Fortran и С.
3 Литература
- ↑ Ильин В.А., Ким Г.Д. "Линейная алгебра и аналитическая геометрия".
- ↑ Wikipedia: QR-algorithm
- ↑ Метод Хаусхолдера (отражений) QR-разложения квадратной матрицы
- ↑ Метод Гивенса (вращений) QR-разложения квадратной матрицы
- ↑ Тыртышников Е.Е. "Методы численного анализа" — М., Академия, 2007. - 320 c.
- ↑ R. Granat, Bo Kagstrom, D. Kressner "LAPACK Working Note #216: A novel parallel QR algorithm for hybrid distributed memory HPC systems".
- ↑ Метод Хаусхолдера (отражений) приведения матрицы к хессенберговой (почти треугольной) форме
- ↑ Перейти обратно: 8,0 8,1 Бахвалов Н.С., Жидков Н.П., Кобельков. Г.М. "Численный методы"— 6-е изд. — М. : БИНОМ. Лаборатория знаний, 2008. — 636 с.