Алговики - Вклад [ru]

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

2025-11-04T19:13:19Z

BoldyrevVV MSU: add references

{{algorithm
| name = Матричный крестовой метод
| serial_complexity = <math>O((m + n) r^2)</math>
| pf_height = <math>O(r p)</math>
| pf_width = <math>O\left(\left\lceil \frac{m}{p} \right\rceil + \left\lceil \frac{n}{p} \right\rceil\right) r</math>
| input_data = <math>O((m + n) r)</math> (вычисляются лишь отдельные элементы матрицы)
| output_data = <math>O((m + n) r)</math> (факторы <math>U \in \mathbb{R}^{m \times r}</math>, <math>V \in \mathbb{R}^{r \times n}</math>)
}}

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

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

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

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

Ключевая идея основана на так называемом '''крестовом разложении'''<ref>S. A. Goreinov, E. E. Tyrtyshnikov, and N. L. Zamarashkin,“Pseudoskeleton Approximations of Matrices,”Dokl. Akad. Nauk343(2), 151–152 (1995) [Dokl. Math.52(1), 18–19 (1995)]</ref><ref>E. Tyrtyshnikov, “Incomplete Cross Approximation in theMosaic-Skeleton Method,” Computing64(4),367–380 (2000).</ref>.: ''Если'' известны <math>r</math> строк и <math>r</math> столбцов матрицы, пересечение которых образует невырожденную подматрицу, ''то'' вся матрица может быть восстановлена через интерполяцию по этим строкам и столбцам. На практике используется '''адаптивная схема''' (ACA и ACA+)<ref>Zhao, Kezhong & Vouvakis, M.N. & Lee, Jin-Fa. (2005). The Adaptive Cross Approximation Algorithm for Accelerated Method of Moments Computations of EMC Problems. Electromagnetic Compatibility, IEEE Transactions on. 47. 763 - 773. 10.1109/TEMC.2005.857898.</ref> с выбором элемента максимальной невязки на каждой итерации (аналог жадного алгоритма), что позволяет обойтись без поиска подматрицы максимального объёма.

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

В условиях современных суперкомпьютеров и больших размерностей (<math>m, n \sim 10^5</math> и выше) естественным шагом становится параллельная реализация, предложенная в работе Д. А. Желткова и Е. Е. Тыртышникова (2015)<ref>Желтков, Д.А. Параллельная реализация матричного крестового метода // Д.А. Желтков, Е.Е. Тыртышников // ВЫЧИСЛИТЕЛЬНЫЕ МЕТОДЫ И ПРОГРАММИРОВАНИЕ: НОВЫЕ ВЫЧИСЛИТЕЛЬНЫЕ ТЕХНОЛОГИИ. - 2015. № 3 16. - 369-375.</ref>. Она распределяет строки и столбцы матрицы между процессорами и использует коллективные операции для глобального выбора доминирующего элемента невязки.

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

В параллельной версии каждый процессор хранит часть строк и столбцов факторов <math>U</math> и <math>V</math>. На каждой итерации:

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

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

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

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

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

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

'''Алгоритм 2. Крестовый метод (<math>k = 1</math>)'''

''Инициализация:''
<math>I \gets {1,\dots,m}</math>, <math>J \gets {1,\dots,n}</math>, <math>r \gets 0</math>.
 ''Цикл:''
 1. <math>r \gets r + 1</math>;
 2. Найти строку с максимальной невязкой в первом доступном столбце:
 
:<math>i \gets \arg\max_{i \in I} |A_{i, J(1)} - [UV]{i, J(1)}|</math>;
 3. Найти столбец с максимальной невязкой в найденной строке:
 
:<math>j \gets \arg\max{j \in J} |A_{ij} - [UV]{ij}|</math>;
 4. Проверить критерий остановки:
 
:''Если''
:<math>|(A - UV){ij}| \cdot \sqrt{(m - r)(n - r)} < \varepsilon |UV|F</math>, ''то''
::<math>r \gets r - 1</math>,
::выйти из цикла;
 5. Обновить факторы:
:<math>u \gets \dfrac{A^{j} - [UV]^{j}}{\sqrt{\vert A_{ij}\vert}}</math>,
:<math>v \gets \dfrac{(A_{i} - [UV]_{i}) \cdot \sqrt{\vert A_{ij}\vert}}{A_{ij}}</math>;
 6. Расширить матрицы:
:<math>U \gets [Uu]</math>,
:<math>V \gets \begin{bmatrix} V \\ v \end{bmatrix}</math>;
 7. Исключить использованные индексы:
:<math>I \gets I \setminus \left\{ {i} \right\}</math>,
:<math>J \gets J \setminus \left\{ {j} \right\}</math>;
 Вернуть <math>\left( U, V, r \right)</math>.
 

 Здесь <math>A^{(j)}</math> обозначает <math>j</math>-й столбец матрицы <math>A</math>, а <math>A_{(i)}</math> — <math>i</math>-ю строку. Норма Фробениуса <math>|UV|_F</math> пересчитывается инкрементально за <math>O((m + n) r)</math> операций, что позволяет сохранить общую сложность алгоритма на уровне <math>O((m + n) r^2)</math>.

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

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

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

:<math>{b}_{m} \gets \frac{m}{p}</math>
:<math>{b}_{n} \gets \frac{n}{p}</math>

''Если''
:<math>rnk = 0</math>, ''то''
::<math>J_* \gets {1, \dots, n}</math>
''КонецЕсли''

:<math>I \gets { bm \cdot rnk + 1, \dots, bm \cdot (rnk + 1) }</math>
:<math>J \gets { bn \cdot rnk + 1, \dots, bn \cdot (rnk + 1) }</math>
:<math>r \gets 0</math>

''повторять''
:<math>r \gets r + 1</math>

''Eсли'' <math>rnk = 0</math>, ''то''
::<math>j \gets J^*(1)</math>
''КонецЕсли''
:<math>\texttt{Broadcast}(0, j, j)</math>

:<math>\texttt{Broadcast}\left(\left\lfloor \frac{j - 1}{p} \right\rfloor, V^j, v\right)</math>

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

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

''Eсли'' <math>rnk = 0</math>, ''то''
::<math>{i, column_maximum} \gets \max_{k \in [1; p]} (max_i[k])</math>
''КонецЕсли''

:<math>\texttt{Broadcast}(0, i, i)</math>

:<math>\texttt{Broadcast}\left(\left\lfloor \frac{i - 1}{p} \right\rfloor, U_i, u\right)</math>

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

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

''Eсли'' <math>rnk = 0</math>, ''то''
::<math>{j, maximum} \gets \max_{k \in [1; p]} (maxj[k])</math>
''КонецЕсли''

:<math>\texttt{Broadcast}(0, \left\{j, maximum\right\}, \left\{j, maximum\right\})</math>

:<math>\texttt{Broadcast}\left(\left\lfloor \frac{j - 1}{p} \right\rfloor, V^j, v\right)</math>

:<math>u_{upd} \gets \dfrac{(A^j - U v)[b_m \cdot rnk + 1 : b_m \cdot (rnk + 1)]}{\sqrt{ \vert maximum \vert }}</math>
:<math>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>U \gets [U u_{upd}]</math>
:<math>V \gets \begin{bmatrix} V \\ v_{upd} \end{bmatrix}</math>

:<math>I \gets I \setminus \left\{ i \right\}</math>
:<math>J \gets J \setminus \left\{ j \right\}</math>

''Eсли'' <math>rnk = 0</math>, ''то''
::<math>J_* \gets J_* \setminus \left\{ j \right\}</math>
''КонецЕсли''

''до тех пор, пока''
:<math>\operatorname{tol} \vert \vert UV \vert \vert _F < |maximum| \cdot \sqrt{(m - r)(n - r)}</math>

'''Вернуть''' <math>(U, V, r)</math>
 
 
Функция <math>\texttt{Broadcast}\left( i,a,b \right)</math> рассылает переменную a с процессора <math>i</math> в переменную <math>b</math> всех процессоров; функция <math>\texttt{Gather}\left( i,a,b \right)</math> собирает переменную <math>a</math> всех процессоров в переменную <math>b</math> (таким образом, <math>b</math> массив из переменных <math>a</math> процессора <math>i</math>.
 
Пересчёт нормы текущего приближения легко осуществляется параллельно в условиях реализованного в алгоритме распределения по процессорам матриц <math>U</math> и <math>V</math>. Каждый процессор может вычислять векторы
:<math>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>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>x = \sum_{rnk=0}^{p} x_{rnk}</math>, <math>y = \sum_{rnk=0}^{p} y_{rnk}</math>,
каждую с помощью одной операции коллективного взаимодействия процессоров. Вычисляем их скалярное произведение <math>(x, y)</math>. Очевидно, что <math>(x, y) = (U^T u, V v^T)</math>.

Слагаемые <math>(V v^T, U^T u)</math> и <math>||u||_2^2 ||v||_2^2</math>, требуемые для пересчёта нормы, вычисляются аналогично.

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

Каждый процессор, кроме первого, в данном алгоритме выполняет <math>O\left(\left(\left\lceil \frac{m}{p} \right\rceil + \left\lceil \frac{n}{p} \right\rceil\right) r^2\right)</math> операций и вычисляет <math>O\left(\left(\left\lceil \frac{m}{p} \right\rceil + \left\lceil \frac{n}{p} \right\rceil\right) r\right)</math> элементов матрицы. Первый процессор дополнительно выполняет <math>O(rp)</math> операций.

== Литература ==

<references />

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

2025-11-04T18:48:11Z

BoldyrevVV MSU: add comments

{{algorithm
| name = Матричный крестовой метод
| serial_complexity = <math>O((m + n) r^2)</math>
| pf_height = <math>O(r p)</math>
| pf_width = <math>O\left(\left\lceil \frac{m}{p} \right\rceil + \left\lceil \frac{n}{p} \right\rceil\right) r</math>
| input_data = <math>O((m + n) r)</math> (вычисляются лишь отдельные элементы матрицы)
| output_data = <math>O((m + n) r)</math> (факторы <math>U \in \mathbb{R}^{m \times r}</math>, <math>V \in \mathbb{R}^{r \times n}</math>)
}}

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

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

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

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

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

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

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

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

В параллельной версии каждый процессор хранит часть строк и столбцов факторов <math>U</math> и <math>V</math>. На каждой итерации:

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

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

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

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

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

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

'''Алгоритм 2. Крестовый метод (<math>k = 1</math>)'''

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

 Здесь <math>A^{(j)}</math> обозначает <math>j</math>-й столбец матрицы <math>A</math>, а <math>A_{(i)}</math> — <math>i</math>-ю строку. Норма Фробениуса <math>|UV|_F</math> пересчитывается инкрементально за <math>O((m + n) r)</math> операций, что позволяет сохранить общую сложность алгоритма на уровне <math>O((m + n) r^2)</math>.

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

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

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

:<math>{b}_{m} \gets \frac{m}{p}</math>
:<math>{b}_{n} \gets \frac{n}{p}</math>

:<math>\color{Red}\text{если}</math>
:<math>rnk = 0</math>, то
::<math>J_* \gets {1, \dots, n}</math>

:<math>I \gets { bm \cdot rnk + 1, \dots, bm \cdot (rnk + 1) }</math>
:<math>J \gets { bn \cdot rnk + 1, \dots, bn \cdot (rnk + 1) }</math>
:<math>r \gets 0</math>

''повторять''
:<math>r \gets r + 1</math>

'''Eсли''' <math>rnk = 0</math>, то
::<math>j \gets J^*(1)</math>
'''КонецЕсли'''
:<math>\texttt{Broadcast}(0, j, j)</math>

:<math>\texttt{Broadcast}\left(\left\lfloor \frac{j - 1}{p} \right\rfloor, V^j, v\right)</math>

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

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

'''Eсли''' <math>rnk = 0</math>, то
::<math>{i, column_maximum} \gets \max_{k \in [1; p]} (max_i[k])</math>
'''КонецЕсли'''

:<math>\texttt{Broadcast}(0, i, i)</math>

:<math>\texttt{Broadcast}\left(\left\lfloor \frac{i - 1}{p} \right\rfloor, U_i, u\right)</math>

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

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

'''Eсли''' <math>rnk = 0</math>, то
::<math>{j, maximum} \gets \max_{k \in [1; p]} (maxj[k])</math>
'''КонецЕсли'''

:<math>\texttt{Broadcast}(0, \left\{j, maximum\right\}, \left\{j, maximum\right\})</math>

:<math>\texttt{Broadcast}\left(\left\lfloor \frac{j - 1}{p} \right\rfloor, V^j, v\right)</math>

:<math>u_{upd} \gets \dfrac{(A^j - U v)[b_m \cdot rnk + 1 : b_m \cdot (rnk + 1)]}{\sqrt{ \vert maximum \vert }}</math>
:<math>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>U \gets [U u_{upd}]</math>
:<math>V \gets \begin{bmatrix} V \\ v_{upd} \end{bmatrix}</math>

:<math>I \gets I \setminus \left\{ i \right\}</math>
:<math>J \gets J \setminus \left\{ j \right\}</math>

'''Eсли''' <math>rnk = 0</math>, то
::<math>J_* \gets J_* \setminus \left\{ j \right\}</math>
'''КонецЕсли'''

''до тех пор, пока''
:<math>\operatorname{tol} \vert \vert UV \vert \vert _F < |maximum| \cdot \sqrt{(m - r)(n - r)}</math>

'''вернуть''' <math>(U, V, r)</math>
 
 
Функция <math>\texttt{Broadcast}\left( i,a,b \right)</math> рассылает переменную a спроцессора <math>i</math> в переменную <math>b</math> всех процессоров; функция <math>\texttt{Gather}\left( i,a,b \right)</math> собирает переменную <math>a</math> всех процессоров в переменную <math>b</math> (таким образом, <math>b</math> массив из переменных <math>a</math> процессора <math>i</math>.
 
Пересчёт нормы текущего приближения легко осуществляется параллельно в условиях реализованного в алгоритме распределения по процессорам матриц <math>U</math> и <math>V</math>. Каждый процессор может вычислять векторы
:<math>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>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>x = \sum_{rnk=0}^{p} x_{rnk}</math>, <math>y = \sum_{rnk=0}^{p} y_{rnk}</math>,
каждую с помощью одной операции коллективного взаимодействия процессоров. Вычисляем их скалярное произведение <math>(x, y)</math>. Очевидно, что <math>(x, y) = (U^T u, V v^T)</math>.

Слагаемые <math>(V v^T, U^T u)</math> и <math>||u||_2^2 ||v||_2^2</math>, требуемые для пересчёта нормы, вычисляются аналогично.

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

Каждый процессор, кроме первого, в данном алгоритме выполняет <math>O\left(\left(\left\lceil \frac{m}{p} \right\rceil + \left\lceil \frac{n}{p} \right\rceil\right) r^2\right)</math> операций и вычисляет <math>O\left(\left(\left\lceil \frac{m}{p} \right\rceil + \left\lceil \frac{n}{p} \right\rceil\right) r\right)</math> элементов матрицы. Первый процессор дополнительно выполняет <math>O(rp)</math> операций.

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

2025-11-04T18:37:30Z

BoldyrevVV MSU: upd parallel impl

{{algorithm
| name = Матричный крестовой метод
| serial_complexity = <math>O((m + n) r^2)</math>
| pf_height = <math>O(r p)</math>
| pf_width = <math>O\left(\left\lceil \frac{m}{p} \right\rceil + \left\lceil \frac{n}{p} \right\rceil\right) r</math>
| input_data = <math>O((m + n) r)</math> (вычисляются лишь отдельные элементы матрицы)
| output_data = <math>O((m + n) r)</math> (факторы <math>U \in \mathbb{R}^{m \times r}</math>, <math>V \in \mathbb{R}^{r \times n}</math>)
}}

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

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

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

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

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

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

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

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

В параллельной версии каждый процессор хранит часть строк и столбцов факторов <math>U</math> и <math>V</math>. На каждой итерации:

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

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

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

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

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

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

'''Алгоритм 2. Крестовый метод (<math>k = 1</math>)'''

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

 Здесь <math>A^{(j)}</math> обозначает <math>j</math>-й столбец матрицы <math>A</math>, а <math>A_{(i)}</math> — <math>i</math>-ю строку. Норма Фробениуса <math>|UV|_F</math> пересчитывается инкрементально за <math>O((m + n) r)</math> операций, что позволяет сохранить общую сложность алгоритма на уровне <math>O((m + n) r^2)</math>.

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

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

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

:<math>{b}_{m} \gets \frac{m}{p}</math>
:<math>{b}_{n} \gets \frac{n}{p}</math>

:<math>\color{Red}\text{если}</math>
:<math>rnk = 0</math>, то
::<math>J_* \gets {1, \dots, n}</math>

:<math>I \gets { bm \cdot rnk + 1, \dots, bm \cdot (rnk + 1) }</math>
:<math>J \gets { bn \cdot rnk + 1, \dots, bn \cdot (rnk + 1) }</math>
:<math>r \gets 0</math>

''повторять''
:<math>r \gets r + 1</math>

:'''если'''
:<math>rnk = 0</math>, то
::<math>j \gets J^*(1)</math>
:<math>\texttt{Broadcast}(0, j, j)</math>

:<math>\texttt{Broadcast}\left(\left\lfloor \frac{j - 1}{p} \right\rfloor, V^j, v\right)</math>

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

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

:'''если'''
:<math>rnk = 0</math>, то
::<math>{i, column_maximum} \gets \max_{k \in [1; p]} (max_i[k])</math>

:<math>\texttt{Broadcast}(0, i, i)</math>

:<math>\texttt{Broadcast}\left(\left\lfloor \frac{i - 1}{p} \right\rfloor, U_i, u\right)</math>

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

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

:'''если'''
:<math>rnk = 0</math>, то
::<math>{j, maximum} \gets \max_{k \in [1; p]} (maxj[k])</math>
:<math>\texttt{Broadcast}(0, \left\{j, maximum\right\}, \left\{j, maximum\right\})</math>

:<math>\texttt{Broadcast}\left(\left\lfloor \frac{j - 1}{p} \right\rfloor, V^j, v\right)</math>

:<math>u_{upd} \gets \dfrac{(A^j - U v)[b_m \cdot rnk + 1 : b_m \cdot (rnk + 1)]}{\sqrt{ \vert maximum \vert }}</math>
:<math>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>U \gets [U u_{upd}]</math>
:<math>V \gets \begin{bmatrix} V \\ v_{upd} \end{bmatrix}</math>

:<math>I \gets I \setminus \left\{ i \right\}</math>
:<math>J \gets J \setminus \left\{ j \right\}</math>

:'''если''' <math>rnk = 0</math>, то
::<math>J_* \gets J_* \setminus \left\{ j \right\}</math>

''до тех пор, пока''
:<math>\operatorname{tol} \vert \vert UV \vert \vert _F < |maximum| \cdot \sqrt{(m - r)(n - r)}</math>

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

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

2025-11-04T16:46:23Z

BoldyrevVV MSU: fix

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

2025-11-04T16:39:08Z

BoldyrevVV MSU: fix formula

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

2025-11-04T16:32:37Z

BoldyrevVV MSU: fix formula

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

2025-11-04T16:29:47Z

BoldyrevVV MSU: fix formula

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

2025-11-04T16:23:11Z

BoldyrevVV MSU: add description

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

2025-11-03T17:53:52Z

BoldyrevVV MSU: init

init