Участник:BoldyrevVV MSU/Параллельная реализация матричного крестового метода

Матричный крестовой метод
Последовательный алгоритм
Последовательная сложность	[math]\displaystyle{ O((m + n) r^2) }[/math]
Объём входных данных	[math]\displaystyle{ O((m + n) r) }[/math] (вычисляются лишь отдельные элементы матрицы)
Объём выходных данных	[math]\displaystyle{ O((m + n) r) }[/math] (факторы [math]\displaystyle{ U \in \mathbb{R}^{m \times r} }[/math], [math]\displaystyle{ V \in \mathbb{R}^{r \times n} }[/math])
Параллельный алгоритм
Высота ярусно-параллельной формы	[math]\displaystyle{ O(r p) }[/math]
Ширина ярусно-параллельной формы	[math]\displaystyle{ O\left(\left\lceil \frac{m}{p} \right\rceil + \left\lceil \frac{n}{p} \right\rceil\right) r }[/math]

Основные авторы описания: Д.А. Желтков, Е.Е. Тыртышников

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

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

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

Ключевая идея основана на так называемом крестовом разложении: если известны [math]\displaystyle{ r }[/math] строк и [math]\displaystyle{ r }[/math] столбцов матрицы, пересечение которых образует невырожденную подматрицу, то вся матрица может быть восстановлена через интерполяцию по этим строкам и столбцам. На практике используется адаптивная схема с выбором элемента максимальной невязки на каждой итерации (аналог жадного алгоритма), что позволяет обойтись без поиска подматрицы максимального объёма.

Важнейшее преимущество метода — малая вычислительная стоимость: для построения приближения ранга [math]\displaystyle{ r }[/math] требуется [math]\displaystyle{ O((m + n) r) }[/math] вычислений элементов матрицы и [math]\displaystyle{ O((m + n) r^2) }[/math] арифметических операций.

В условиях современных суперкомпьютеров и больших размерностей ([math]\displaystyle{ m, n \sim 10^5 }[/math] и выше) естественным шагом становится параллельная реализация, предложенная в работе Желткова и Тыртышникова (2015). Она распределяет строки и столбцы матрицы между процессорами и использует коллективные операции для глобального выбора доминирующего элемента невязки.

Особенности параллельной реализации

В параллельной версии каждый процессор хранит часть строк и столбцов факторов [math]\displaystyle{ U }[/math] и [math]\displaystyle{ V }[/math]. На каждой итерации:

процессоры совместно находят глобальный максимум невязки с помощью операций [math]\displaystyle{ \texttt{Gather} }[/math] и [math]\displaystyle{ \texttt{Broadcast} }[/math]; соответствующие строки и столбцы обновляются локально; норма Фробениуса текущего приближения пересчитывается распределённо, без полного сбора данных. Эта стратегия позволяет сохранить низкую коммуникационную сложность: каждый процессор выполняет [math]\displaystyle{ O(r) }[/math] обменов, а длина каждого сообщения не превышает [math]\displaystyle{ r }[/math]. При этом трудоёмкость на один процессор снижается до [math]\displaystyle{ O\left(\left(\left\lceil \frac{m}{p} \right\rceil + \left\lceil \frac{n}{p} \right\rceil\right) r^2\right) }[/math], а количество вычисляемых элементов — до [math]\displaystyle{ O\left(\left(\left\lceil \frac{m}{p} \right\rceil + \left\lceil \frac{n}{p} \right\rceil\right) r\right) }[/math].

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

Пусть задана матрица [math]\displaystyle{ A \in \mathbb{R}^{m \times n} }[/math]. Цель — построить её малоранговое приближение вида [math]\displaystyle{ A \approx UV }[/math], где [math]\displaystyle{ U \in \mathbb{R}^{m \times r} }[/math], [math]\displaystyle{ V \in \mathbb{R}^{r \times n} }[/math], а ранг [math]\displaystyle{ r \ll \min(m,n) }[/math].

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

Алгоритм 1. Схема крестового метода

Инициализация:  [math]\displaystyle{ I \gets {1,\dots,m} }[/math], [math]\displaystyle{ J \gets {1,\dots,n} }[/math], [math]\displaystyle{ r \gets 0 }[/math].
Повторять:
 а. [math]\displaystyle{ r \gets r + 1 }[/math];
 б. Найти подмножества индексов [math]\displaystyle{ I_* \subset I }[/math], [math]\displaystyle{ J_* \subset J }[/math] (размера [math]\displaystyle{ k }[/math]), такие что подматрица [math]\displaystyle{ A[I_*, J_*] }[/math] имеет достаточно большой объём (процедура [math]\displaystyle{ \texttt{maxvol} }[/math]);
 в. Выбрать столбцы: [math]\displaystyle{ C_r \gets A[:, J_*] }[/math];
 г. Выбрать строки: [math]\displaystyle{ R_r \gets A[I_*, :] }[/math];
 д. Выделить ядро: [math]\displaystyle{ \hat{A}r \gets A[I, J_] }[/math];
 е. Обновить текущую невязку: [math]\displaystyle{ A \gets A - C_r \hat{A}r^{-1} R_r }[/math];
 ж. Исключить использованные индексы: [math]\displaystyle{ I \gets I \setminus I }[/math], [math]\displaystyle{ J \gets J \setminus J_ }[/math];
До выполнения критерия остановки.
В рассматриваемой работе используется частный случай этой схемы при [math]\displaystyle{ k = 1 }[/math], то есть на каждой итерации выбирается один элемент с максимальной по модулю невязкой. В этом случае подматрица [math]\displaystyle{ \hat{A}_r }[/math] — скаляр, а обновление факторов происходит по явным формулам.

Алгоритм 2. Крестовый метод ([math]\displaystyle{ k = 1 }[/math])

Инициализация:  [math]\displaystyle{ I \gets {1,\dots,m} }[/math], [math]\displaystyle{ J \gets {1,\dots,n} }[/math], [math]\displaystyle{ r \gets 0 }[/math], [math]\displaystyle{ U \gets [] }[/math], [math]\displaystyle{ V \gets [] }[/math]. Цикл:  
а. [math]\displaystyle{ r \gets r + 1 }[/math];  
б. Найти строку с максимальной невязкой в первом доступном столбце:   [math]\displaystyle{ i \gets \arg\max_{i \in I} |A_{i, J(1)} - [UV]{i, J(1)}| }[/math];  
в. Найти столбец с максимальной невязкой в найденной строке:   [math]\displaystyle{ j \gets \arg\max{j \in J} |A_{ij} - [UV]{ij}| }[/math];  
г. Проверить критерий остановки:   если [math]\displaystyle{ |(A - UV){ij}| \cdot \sqrt{(m - r)(n - r)} \lt \varepsilon |UV|F }[/math], то    [math]\displaystyle{ r \gets r - 1 }[/math], выйти из цикла;  
д. Обновить факторы:   [math]\displaystyle{ u \gets \dfrac{A^{(j)} - [UV]^{(j)}}{\sqrt{|(A - UV){ij}|}} }[/math],   [math]\displaystyle{ v \gets \dfrac{(A_{(i)} - [UV]{(i)}) \cdot \sqrt{|(A - UV){ij}|}}{(A - UV)_{ij}} }[/math];  
е. Расширить матрицы:   [math]\displaystyle{ U \gets [U ; u] }[/math], [math]\displaystyle{ V \gets \begin{bmatrix} V \ v \end{bmatrix} }[/math];  
ж. Исключить использованные индексы:   [math]\displaystyle{ I \gets I \setminus {i} }[/math], [math]\displaystyle{ J \gets J \setminus {j} }[/math];
Вернуть [math]\displaystyle{ U, V, r }[/math].
Здесь [math]\displaystyle{ A^{(j)} }[/math] обозначает [math]\displaystyle{ j }[/math]-й столбец матрицы [math]\displaystyle{ A }[/math], а [math]\displaystyle{ A_{(i)} }[/math] — [math]\displaystyle{ i }[/math]-ю строку. Норма Фробениуса [math]\displaystyle{ |UV|_F }[/math] пересчитывается инкрементально за [math]\displaystyle{ O((m + n) r) }[/math] операций, что позволяет сохранить общую сложность алгоритма на уровне [math]\displaystyle{ O((m + n) r^2) }[/math].

В параллельной реализации выбор индексов [math]\displaystyle{ i }[/math] и [math]\displaystyle{ j }[/math] выполняется совместно всеми процессорами с помощью коллективных операций (gather/broadcast), а обновление [math]\displaystyle{ U }[/math] и [math]\displaystyle{ V }[/math] — локально на каждом процессоре для его части строк/столбцов.