Просмотр кода страницы Участник:Илья Карандеев(ВТМ, 403)/QR-Факторизация методом Хаусхолдера




== Свойства и структура алгоритма ==
== Программная реализация алгоритма ==


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

==== Реализация без распараллеливания ====
Ниже приведена моя реализация QR-факторизации методом Хаусхолдера. Она написана с помощью библиотеки cblas - 
оболочки замечательной фортрановской библиотеки BLAS для си. Это для одного процесса,
нераспараллеленная программа. Ее ядро (то что отнимает больше всего времени при счете) - функции
из cblas на 52-58 строке. После факторизации производится проверка ||Q * R - A|| / ||A||. Q - ортогональная матрица, R -
верхняя треугольная. Их вычисление цель реализации алгоритма. Ну и Q * R должно примерно равняться A, это следует из теории.
Вернее, без погрешности вычисления, в точности Q * R = A.  Вычисленный порядок погрешности на персональном компьютере
примерно 10^(-14), что есть довольно хорошо.


[[Файл:Code1.png|слева]]
[[Файл:Code2.png|слева]]

<br><br><br><br><br><br><br><br><br>
<br><br><br><br><br><br><br><br><br>
<br><br><br><br><br><br><br><br><br>
<br><br><br><br><br><br><br><br><br>
<br><br><br><br><br><br><br><br><br>
<br><br><br><br><br><br><br><br><br>
<br><br><br><br><br><br><br><br><br>
<br><br><br><br><br><br><br><br><br>
<br><br><br><br><br><br><br><br><br>
<br><br><br><br><br><br><br><br><br>
<br><br><br><br><br><br><br><br><br>
<br><br><br><br><br><br><br><br><br>

==== Распараллеленная реализация ====
А это распараллеленная реализация. Распараллеливание осуществляется при помощи библиотекки open MPI. Данные у каждого процесса 
свои. Они пересылаются друг другу, что и следует из названия библиотеки: "Message passing interface".  В этой реализации проверка
||Q*R - A|| / ||A|| выполняется только при условии, что размер матрицы кратен числу процессов. Написать проверку в этом случае
сложнее из-за того, что нужно пересылать размер отдела матрицы, которую обрабатывает процесс (а он, вообще говоря, может отличаться
у разных процессво, пусть и не более чем на единицу).
<br><br>
В изначальной версии программы ядром (тем, что занимает бОльшую часть времени при исполнении) являлась, как ни странно,
пересылка между процессами. Это узкое место было исправлено: рассылка одного процесса всем остальным стала исполняться
не при помощи функции MPI_Isend, а при помощи специальной функции - MPI_Bcast. Начальная версия программы, хоть и работала,
но не только не приносила выигрыша в скорости, но даже замедляла распараллеленную версию относительно исходной, однопроцессной. 

[[Файл:Code3.png|Слева]]
[[Файл:Code4.png|Слева]]
[[Файл:Code5.png|Слева]]
[[Файл:Code6.png|Слева]]



=== Локальность данных и вычислений ===
=== Возможные способы и особенности параллельной реализации алгоритма ===
=== Масштабируемость алгоритма и его реализации ===
<br><br>

Вот исследование масштабируемости параллельной реализации QR-факторизации квадртной матрицы методом Хаусхолдера. 
Оно проводилось на суперкомпьютере "Ломоносов-2" суперкомпьютерного комплекса МГУ. В первом эксперименте,
эксперименте на сильную масштабируемость, бралась матрица размером 14000 на 14000. Количество процессов в
эксперименте рассматривалось таким: 1, 2, 4, 8. 14, 28, 56, 112, 224, 448, 700. Этот выбор обусловлен тем,
что в одном узле на суперкомпьютере "Ломоносов-2" 14 ядер, а доступное количество узлов автору было равным
50. То есть максимум можно было задействовать 700 ядер.


[[Файл:My2.png.png|справа|Обычная шкала]]
[[Файл:My1.png|слева|Логарифмическая шкала]]

<br><br><br><br><br><br><br><br><br>
<br><br><br><br><br><br><br><br><br>
<br><br><br><br><br><br><br><br><br>


Следующая серия экспериментов на масштабируемость была проведена на 1, 2, 4, 14, 28 и 42 ядрах. Для матриц
от 200 до 8500, выбранных более-менее произвольно. Вот график в разных проекциях, полученный по результатам экспериментальных
запусков программы на  суперкомпьютере "Ломоносов-2".

[[Файл:My.3.png|слева]]
[[Файл:My4.png|слева]]
[[Файл:Figure_3.png|слева]]

<br><br><br><br><br><br><br><br><br>
<br><br><br><br><br><br><br><br><br>
<br><br><br><br><br><br><br><br><br>
<br><br><br><br><br><br><br><br><br>
<br><br><br><br><br><br><br><br><br>
<br><br><br><br><br><br><br><br><br>
<br><br><br><br><br><br><br><br><br>
<br><br><br><br><br><br><br><br><br>
<br><br><br><br><br><br><br><br><br>
<br><br><br><br><br><br><br><br><br>
<br><br><br><br><br><br><br><br><br>
<br><br><br><br><br><br><br><br><br>
<br><br><br><br><br><br><br><br><br>
<br><br><br><br><br><br><br><br><br>
<br><br><br><br><br><br><br><br><br>
<br><br><br><br><br><br><br><br><br>

В целом, все работает объяснимо: больше размерность матрицы → больше время счета. Больше процессов → меньше время счета.  
Можно заметить, что на малых размерах матрицы при большом количестве процессов начинает падать эффективность. Возможно, это из-за накладных расходов из-за обмена сообщениями между процессами.

<br><br><br><br>
=== Динамические характеристики и эффективность реализации алгоритма ===
=== Выводы для классов архитектур ===
=== Существующие реализации алгоритма ===
Большинство пакетов: LINPACK, LAPACK, SCALAPACK и прочие - используют для QR-разложения матриц именно метод Хаусхолдера, правда, в различных модификациях (обычно с использованием BLAS). Существует большая подборка исследовательских работ по блочным версиям. 

Можно найти и много пользовательских реализаций:

На языке Python:
https://www.quantstart.com/articles/QR-Decomposition-with-Python-and-NumPy/
<br>
https://gist.github.com/Hsankesara/cd35edb30825df19f182a6ecf96e126e

На языке R:
https://rpubs.com/aaronsc32/qr-decomposition-householder - R

Для MatLab:
https://www.mathworks.com/matlabcentral/answers/169648-qr-factorization-using-householder-transformations

== Литература ==
В. Б. Андреев: "Численные методы"