Участник: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] столбцов матрицы, пересечение которых образует невырожденную подматрицу, то вся матрица может быть восстановлена через интерполяцию по этим строкам и столбцам. На практике используется адаптивная схема (ACA и ACA+) с выбором элемента максимальной невязки на каждой итерации (аналог жадного алгоритма), что позволяет обойтись без поиска подматрицы максимального объёма.

Важнейшее преимущество метода — малая вычислительная стоимость: для построения приближения ранга [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].
Повторять:

1. [math]\displaystyle{ r \gets r + 1 }[/math];
   
2. Найти подмножества индексов [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]);
   
3. Выбрать столбцы: [math]\displaystyle{ C_r \gets A[:, J_*] }[/math];
   
4. Выбрать строки: [math]\displaystyle{ R_r \gets A[I_*, :] }[/math];
   
5. Выделить ядро: [math]\displaystyle{ \hat{A_i} \gets A[:, I_* J_*] }[/math];
   
6. Обновить текущую невязку: [math]\displaystyle{ A \gets A - C_r \hat{A}_{r}^{-1} R_r }[/math];
   
7. Исключить использованные индексы: [math]\displaystyle{ I \gets I \setminus I_* }[/math], [math]\displaystyle{ J \gets J \setminus J_* }[/math];
   

До выполнения критерия остановки.

В результате работы этого алгоритма матрица [math]\displaystyle{ \sum_{i=1}^{r} C_i\hat{A}_{i}^{-1} R_i }[/math] будет приближать матрицу [math]\displaystyle{ A }[/math]. В рассматриваемой работе используется частный случай этой схемы при [math]\displaystyle{ k = 1 }[/math], то есть на каждой итерации выбирается один элемент с максимальной по модулю невязкой. В этом случае подматрица [math]\displaystyle{ \hat{A}_r }[/math] — скаляр, а обновление факторов происходит по явным формулам.

С учeтом данного выбора алгоритм выглядит следующим образом.

Алгоритм 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].
Цикл:
1. [math]\displaystyle{ r \gets r + 1 }[/math];
2. Найти строку с максимальной невязкой в первом доступном столбце:

[math]\displaystyle{ i \gets \arg\max_{i \in I} |A_{i, J(1)} - [UV]{i, J(1)}| }[/math];

3. Найти столбец с максимальной невязкой в найденной строке:

[math]\displaystyle{ j \gets \arg\max{j \in J} |A_{ij} - [UV]{ij}| }[/math];

4. Проверить критерий остановки:

Если

[math]\displaystyle{ |(A - UV){ij}| \cdot \sqrt{(m - r)(n - r)} \lt \varepsilon |UV|F }[/math],
то

[math]\displaystyle{ r \gets r - 1 }[/math],

выйти из цикла;

5. Обновить факторы:

[math]\displaystyle{ u \gets \dfrac{A^{j} - [UV]^{j}}{\sqrt{\vert A_{ij}\vert}} }[/math],

[math]\displaystyle{ v \gets \dfrac{(A_{i} - [UV]_{i}) \cdot \sqrt{\vert A_{ij}\vert}}{A_{ij}} }[/math];

6. Расширить матрицы:

[math]\displaystyle{ U \gets [Uu] }[/math],

[math]\displaystyle{ V \gets \begin{bmatrix} V \\ v \end{bmatrix} }[/math];

7. Исключить использованные индексы:

[math]\displaystyle{ I \gets I \setminus \left\{ {i} \right\} }[/math],

[math]\displaystyle{ J \gets J \setminus \left\{ {j} \right\} }[/math];

Вернуть [math]\displaystyle{ \left( U, V, r \right) }[/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] — локально на каждом процессоре для его части строк/столбцов.

Входные параметры: [math]\displaystyle{ m }[/math] — число строк матрицы, [math]\displaystyle{ n }[/math] — число столбцов, [math]\displaystyle{ A(i,j) }[/math] — функция вычисления элемента матрицы, [math]\displaystyle{ p }[/math] — число процессоров, [math]\displaystyle{ rnk }[/math] — номер текущего процессора (от 0 до [math]\displaystyle{ p-1 }[/math]).

Алгоритм 3. Параллельный крестовый метод

[math]\displaystyle{ {b}_{m} \gets \frac{m}{p} }[/math]

[math]\displaystyle{ {b}_{n} \gets \frac{n}{p} }[/math]

Если

[math]\displaystyle{ rnk = 0 }[/math], то

[math]\displaystyle{ J_* \gets {1, \dots, n} }[/math]

КонецЕсли

[math]\displaystyle{ I \gets { bm \cdot rnk + 1, \dots, bm \cdot (rnk + 1) } }[/math]

[math]\displaystyle{ J \gets { bn \cdot rnk + 1, \dots, bn \cdot (rnk + 1) } }[/math]

[math]\displaystyle{ r \gets 0 }[/math]

повторять

[math]\displaystyle{ r \gets r + 1 }[/math]

Eсли [math]\displaystyle{ rnk = 0 }[/math], то

[math]\displaystyle{ j \gets J^*(1) }[/math]

КонецЕсли

[math]\displaystyle{ \texttt{Broadcast}(0, j, j) }[/math]

[math]\displaystyle{ \texttt{Broadcast}\left(\left\lfloor \frac{j - 1}{p} \right\rfloor, V^j, v\right) }[/math]

[math]\displaystyle{ i \gets \arg\max_{i \in I} \left| A_{ij} - [U v]_i \right| }[/math]

[math]\displaystyle{ \texttt{Gather}\left( 0, \left\{ i, \vert A_{ij} - [U v]_i \vert \right\}, max_i \right) }[/math]

Eсли [math]\displaystyle{ rnk = 0 }[/math], то

[math]\displaystyle{ {i, column_maximum} \gets \max_{k \in [1; p]} (max_i[k]) }[/math]

КонецЕсли

[math]\displaystyle{ \texttt{Broadcast}(0, i, i) }[/math]

[math]\displaystyle{ \texttt{Broadcast}\left(\left\lfloor \frac{i - 1}{p} \right\rfloor, U_i, u\right) }[/math]

[math]\displaystyle{ j \gets \arg\max_{j \in J} \left| A_{ij} - [u V]_j \right| }[/math]

[math]\displaystyle{ \texttt{Gather}(0, \left\{j, \vert A_{ij} - [u V]_j \vert \right\}, max_j) }[/math]

Eсли [math]\displaystyle{ rnk = 0 }[/math], то

[math]\displaystyle{ {j, maximum} \gets \max_{k \in [1; p]} (maxj[k]) }[/math]

КонецЕсли

[math]\displaystyle{ \texttt{Broadcast}(0, \left\{j, maximum\right\}, \left\{j, maximum\right\}) }[/math]

[math]\displaystyle{ \texttt{Broadcast}\left(\left\lfloor \frac{j - 1}{p} \right\rfloor, V^j, v\right) }[/math]

[math]\displaystyle{ u_{upd} \gets \dfrac{(A^j - U v)[b_m \cdot rnk + 1 : b_m \cdot (rnk + 1)]}{\sqrt{ \vert maximum \vert }} }[/math]

[math]\displaystyle{ v_{upd} \gets \dfrac{(A_i - u V)[b_n \cdot rnk + 1 : b_n \cdot (rnk + 1)] \cdot \sqrt{ \vert maximum \vert }}{maximum} }[/math]

[math]\displaystyle{ U \gets [U u_{upd}] }[/math]

[math]\displaystyle{ V \gets \begin{bmatrix} V \\ v_{upd} \end{bmatrix} }[/math]

[math]\displaystyle{ I \gets I \setminus \left\{ i \right\} }[/math]

[math]\displaystyle{ J \gets J \setminus \left\{ j \right\} }[/math]

Eсли [math]\displaystyle{ rnk = 0 }[/math], то

[math]\displaystyle{ J_* \gets J_* \setminus \left\{ j \right\} }[/math]

КонецЕсли

до тех пор, пока

[math]\displaystyle{ \operatorname{tol} \vert \vert UV \vert \vert _F \lt |maximum| \cdot \sqrt{(m - r)(n - r)} }[/math]

Вернуть [math]\displaystyle{ (U, V, r) }[/math]

Функция [math]\displaystyle{ \texttt{Broadcast}\left( i,a,b \right) }[/math] рассылает переменную a с процессора [math]\displaystyle{ i }[/math] в переменную [math]\displaystyle{ b }[/math] всех процессоров; функция [math]\displaystyle{ \texttt{Gather}\left( i,a,b \right) }[/math] собирает переменную [math]\displaystyle{ a }[/math] всех процессоров в переменную [math]\displaystyle{ b }[/math] (таким образом, [math]\displaystyle{ b }[/math] массив из переменных [math]\displaystyle{ a }[/math] процессора [math]\displaystyle{ i }[/math].

Пересчёт нормы текущего приближения легко осуществляется параллельно в условиях реализованного в алгоритме распределения по процессорам матриц [math]\displaystyle{ U }[/math] и [math]\displaystyle{ V }[/math]. Каждый процессор может вычислять векторы

[math]\displaystyle{ x_{rnk} = u^T[b_m \cdot rnk + 1 : b_m \cdot (rnk + 1)] \cdot U[b_m \cdot rnk + 1 : b_m \cdot (rnk + 1), :] }[/math],

[math]\displaystyle{ y_{rnk} = v^T[b_n \cdot rnk + 1 : b_n \cdot (rnk + 1)] \cdot V[:, b_n \cdot rnk + 1 : b_n \cdot (rnk + 1)] }[/math].

Находим суммы [math]\displaystyle{ x = \sum_{rnk=0}^{p} x_{rnk} }[/math], [math]\displaystyle{ y = \sum_{rnk=0}^{p} y_{rnk} }[/math], каждую с помощью одной операции коллективного взаимодействия процессоров. Вычисляем их скалярное произведение [math]\displaystyle{ (x, y) }[/math]. Очевидно, что [math]\displaystyle{ (x, y) = (U^T u, V v^T) }[/math].

Слагаемые [math]\displaystyle{ (V v^T, U^T u) }[/math] и [math]\displaystyle{ ||u||_2^2 ||v||_2^2 }[/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{ O(rp) }[/math] операций.