Участник:OscarFox/Обобщенный прокрустов анализ

Автор: Фостенко О.А.

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

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

Обобщенный прокрустов анализ (Generalized Procrustes Analysis, GPA) - это метод, предназначенный для сравнения и выравнивания нескольких наборов из одного и того же числа точек (форм) в пространстве одной размерности. Цель GPA - найти оптимальные преобразования, включающие вращения, масштабирование и сдвиги, для наилучшего согласования форм друг с другом.

Формулировка похожей задачи встречается в работе B. F. Green 1952 года, аналитическое решение задачи выравнивания двух множеств точек друг относительно друга с учетом только вращений (ортогональная задача Прокруста) было предложено P. H. Schönemann в 1966 году, в 1970 году им же в соавторстве с R. M. Carroll предложено аналитическое решение для модификации, учитывающей, помимо вращений, смещения и масштабирование (классический прокрустов анализ), а непосредственно алгоритм GPA был разработан и опубликован J. C. Gower в 1975 году.

GPA - это итеративный алгоритм: он последовательно согласовывает все формы относительно некоторой «средней формы», обновляя её до тех пор, пока изменения не станут незначительными. К «средней форме» набор из точек приводится с помощью классического прокрустова анализа. Иллюстрация выравнивания, выполняемого с помощью как классического, так и обобщенного прокрустова анализа, приведена на рисунке 1.

Рис. 1. Иллюстрация прокрустова анализа на примере форм, отвечающих трехмерным человеческим позам

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

Описание данных и формул

Входные данные:

• [math]\displaystyle{ \{X_1, X_2,\dots,X_m\} }[/math] - набор входных форм, то есть матриц одного размера: [math]\displaystyle{ X_i \in \mathbb{R}^{n \times d}, \; i=1,\dots,m }[/math]. Каждая матрица [math]\displaystyle{ X_i }[/math] представляет собой набор из [math]\displaystyle{ n }[/math] точек в [math]\displaystyle{ d }[/math]-мерном пространстве. Точки с одинаковыми индексами (т.е. строки матриц с одинаковыми номерами) предполагаются соответствующими друг другу.

Выходные данные:

• Выровненные формы [math]\displaystyle{ \{ X_1', X_2',\dots,X_m'\} }[/math], где [math]\displaystyle{ X_i'\in \mathbb{R}^{n \times d}, \; i=1,\dots,m }[/math].
• Средняя форма (консенсус) [math]\displaystyle{ Y\in \mathbb{R}^{n \times d} }[/math].
• Опционально, коэффициенты масштаба [math]\displaystyle{ s_i\in \mathbb{R} }[/math].

Вспомогательные понятия:

• [math]\displaystyle{ CPA(X, Y): \mathbb{R}^{n \times d} \times \mathbb{R}^{n \times d} \rightarrow \mathbb{R}^{n \times d} }[/math] - функция, выполняющая классический прокрустов анализ, вращая, смещая и масштабируя форму [math]\displaystyle{ X }[/math] так, чтобы она наилучшим образом подходила под форму [math]\displaystyle{ Y }[/math]. Возвращает матрицу из [math]\displaystyle{ \mathbb{R}^{n \times d} }[/math], отвечающую преобразованной форме [math]\displaystyle{ X }[/math].

• [math]\displaystyle{ \| X \|_F = \sqrt{tr(X X^T)} }[/math] - норма Фробениуса матрицы [math]\displaystyle{ X }[/math].

Описание алгоритма (GPA)

Пусть [math]\displaystyle{ X_i' = s_i X_i R_i + \mathbf{1}_n t_i }[/math]. Здесь [math]\displaystyle{ s_i\in \mathbb{R} }[/math] есть константа масштаба, [math]\displaystyle{ R_i\in \mathbb{R}^{d \times d} }[/math] - ортогональная матрица, отвечающая вращению, [math]\displaystyle{ \mathbf{1}_n\in \mathbb{R}^{n \times 1} }[/math] - вектор-столбец из единиц, а [math]\displaystyle{ t_i\in \mathbb{R}^{1 \times d} }[/math] - вектор переноса.

В процессе работы GPA итеративно минимизирует следующую функцию:

[math]\displaystyle{ S = \displaystyle \sum_{i\lt j}^{m} \| X_i' - X_j' \|_F^2 }[/math]

по всевозможным смещениям [math]\displaystyle{ t_i }[/math], вращениям [math]\displaystyle{ R_i }[/math] и константам масштаба [math]\displaystyle{ s_i }[/math] для всех [math]\displaystyle{ i=1,\dots,m }[/math], с условием, что [math]\displaystyle{ \displaystyle \sum_{i=1}^{m} s_i^2 \| X_i \|_F^2 = \sum_{i=1}^{m} \| X_i \|_F^2 }[/math] (чтобы избежать зануления всех [math]\displaystyle{ s_i }[/math]).

Минимизация [math]\displaystyle{ S }[/math] сводится к минимизации [math]\displaystyle{ S_r = \displaystyle \sum_{i=1}^{m} tr(X_i' X_i'^T) - m \cdot tr(Y Y^T) }[/math], где [math]\displaystyle{ Y }[/math] - средняя форма (консенсус).

Алгоритм (GPA)

Шаг 1. Инициализация.

Задаем начальные значения выравниваемых форм, центрируя их:

Вычисляем смещение для центрирования форм: [math]\displaystyle{ \bar{X_i} = \frac{1}{n} \mathbf{1}_n^T X_i }[/math]

Центрируем формы: [math]\displaystyle{ X_i' := X_i - \mathbf{1}_n \bar{X_i}, \; i = 1,\dots,m. }[/math]

Далее нормируем масштаб форм (в совокупности):

[math]\displaystyle{ \displaystyle X_i' := \lambda X_i', \; \lambda = \sqrt{\frac{m}{ \sum_{i=1}^{m} \| X_i' \|_F^2}}, \; i=1,\dots,m. }[/math]

Инициализируем [math]\displaystyle{ s_i := \lambda, \; i=1,\dots,m }[/math].

Задаем начальную среднюю форму, выбирая одну из центрированных форм входного набора (например, первую): [math]\displaystyle{ Y := X_1' }[/math].

Шаг 2. Инициализация функции ошибки.

Выполнить подгонку каждой формы (без учета первой - так как выбрали ее за изначальный консенсус) [math]\displaystyle{ X_i' }[/math] к [math]\displaystyle{ Y }[/math] и оценить функцию ошибки:

[math]\displaystyle{ \forall i=2,\dots,m \; X_i' := CPA(X_i', Y) }[/math];

Пересчитать консенсус: [math]\displaystyle{ \displaystyle Y := \frac{1}{m} \sum_{i=1}^{m} X_i' }[/math];

посчитать начальное значение функции потерь: [math]\displaystyle{ S_r := m(1 - tr(Y Y^T)) }[/math].

Заметим, что на данном шаге в рамках классического прокрустова анализа ([math]\displaystyle{ CPA(X, Y) }[/math]) вычисляется только матрица поворота, но не смещение и масштаб (так как формы уже центрированы, а масштаб вычисляется далее).

Шаг 3. Выравнивание форм относительно средней.

Для всех текущих форм [math]\displaystyle{ X_i', \; i = 1,\dots,m }[/math] выполняется их выравнивание относительно средней формы (консенсуса). Заметим, что здесь также выполняется только поворот, но не смещение и масштабирование:

[math]\displaystyle{ X_i' }[/math] присваивается новое значение: [math]\displaystyle{ X_i' := CPA(X_i', Y) = X_i' R_i, \; i=1,\dots,m }[/math].

Выполняется вычисление нового консенсуса: [math]\displaystyle{ Y^* := \frac{1}{m} \sum_{i=1}^{m} X_i' }[/math];

также выполняется пересчет функции ошибки: [math]\displaystyle{ S_r^* := S_r - m \cdot tr(Y^* (Y^*)^T - Y Y^T) }[/math];

и устанавливаем [math]\displaystyle{ S_r^{**} := S_r^* }[/math].

Шаг 4.

Если масштабирование не требуется, перейти на шаг 6, присвоив [math]\displaystyle{ Y := Y^* }[/math], иначе перейти на следующий шаг.

Шаг 5. Масштабирование.

Для каждого [math]\displaystyle{ i=1,\dots,m }[/math]

выполняется подсчет [math]\displaystyle{ s_i^* := s_i \cdot \sqrt{ \frac{tr(X_i' Y^T)}{tr(X_i' X_i'^T) \cdot tr(Y Y^T)}} }[/math];

далее присваиваем [math]\displaystyle{ X_i' := \frac{s_i^*}{s_i} X_i' }[/math];

устанавливаем [math]\displaystyle{ s_i := s_i^* }[/math];

вычисляем новый консенсус: [math]\displaystyle{ Y^{**} := \frac{1}{m} \sum_{i=1}^{m} X_i' }[/math];

пересчитываем функцию ошибки [math]\displaystyle{ S_r^{**} := S_r^* - m \cdot tr(Y^{**}(Y^{**})^T - Y^* (Y^*)^T) = S_r - m \cdot tr(Y^{**}(Y^{**})^T - YY^T) }[/math];

присваиваем [math]\displaystyle{ Y := Y^{**} }[/math]

Шаг 6. Проверка критерия останова.

Если [math]\displaystyle{ S_r - S_r^{**} \gt \varepsilon }[/math], где [math]\displaystyle{ \varepsilon }[/math] - установленный заранее параметр точности, то установить [math]\displaystyle{ S_r := S_r^{**} }[/math] и перейти на шаг 3, иначе перейти на следующий шаг.

Шаг 7. Завершение итераций.

Вернуть результат: [math]\displaystyle{ Y, \; X_i' }[/math], а также опционально [math]\displaystyle{ s_i, \; i=1,\dots,m }[/math].

Описание подзадачи (CPA)

Задача классического прокрустова анализа (Classical Procrustes Analysis) без масштабирования и сдвигов заключается в нахождении матрицы вращения [math]\displaystyle{ R }[/math], минимизирующей [math]\displaystyle{ \| AR - B \|^2_F }[/math] для некоторых матриц [math]\displaystyle{ A }[/math] и [math]\displaystyle{ B }[/math] одного размера. Такая постановка также называется ортогональной задачей Прокруста (Orthogonal Procrustes Problem, OPA).

У этой задачи есть аналитическое решение. Такая оптимальная матрица поворота находится следующим образом:

Сначала вычисляется произведение [math]\displaystyle{ M = A^T B }[/math].

Далее для нахождения [math]\displaystyle{ R }[/math] требуется выполнить сингулярное разложение [math]\displaystyle{ M: M = U \Sigma V^T }[/math].

Тогда искомая матрица поворота [math]\displaystyle{ R }[/math] выражается как [math]\displaystyle{ R = VU^T }[/math].

Если [math]\displaystyle{ det(R) \lt 0 }[/math], то преобразование выполняется с отражением. Чтобы исключить такие случаи, достаточно домножить последний столбец матрицы [math]\displaystyle{ V }[/math] на [math]\displaystyle{ -1 }[/math]. Тогда для полученного [math]\displaystyle{ R }[/math] будет выполнено ограничение [math]\displaystyle{ det(R)=1 }[/math] при сохранении оптимальности.

Практические аспекты

В данном разделе алгоритм GPA был описан в общем виде, с учетом как его разновидности с масштабированием, так и без него. Важно заметить, что при реализации GPA с масштабированием на практике можно не пересчитывать консенсус и не вычислять функцию ошибки на шаге 3 (далее они все равно будут пересчитаны на шаге 5). Именно такой - оптимизированный - вариант алгоритма GPA с масштабированием будет описан далее.

Вычислительное ядро алгоритма

Вычислительное ядро каждой итерации GPA зависит от размерности выравниваемых матриц [math]\displaystyle{ X_i \in \mathbb{R}^{n \times d} }[/math]. Если [math]\displaystyle{ d \gt n }[/math], то основную сложность представляют сингулярное разложение и конструирование матрицы вращений в рамках классического прокрустова анализа для каждой формы: [math]\displaystyle{ CPA(X_i', Y) = X_i' R_i }[/math] (последовательная сложность - [math]\displaystyle{ O(d^3) }[/math]). Иначе, основные вычисления приходятся на матричные умножения видов [math]\displaystyle{ X^TY }[/math] и [math]\displaystyle{ XR }[/math], где [math]\displaystyle{ X, Y \in \mathbb{R}^{n \times d}, R \in \mathbb{R}^{d \times d} }[/math]. Эти матричные умножения имеют место как в рамках классического прокрустова преобразования, так и при обновлении форм (последовательная сложность - [math]\displaystyle{ O(nd^2) }[/math]).

Макроструктура алгоритма

В данном разделе используются те же обозначения, что и в разделе "Математическое описание алгоритма".

В качестве основных макроопераций можно выделить:

• Операция, сочетающая вычисление матрицы вращения с помощью классического прокрустова анализа [math]\displaystyle{ CPA(X, Y) }[/math] и обновление формы через умножение на полученную матрицу: на каждой итерации GPA (после нулевой) выполняется [math]\displaystyle{ m }[/math] таких операций (где [math]\displaystyle{ m }[/math] - число форм во входных данных). На нулевой же итерации производится [math]\displaystyle{ m-1 }[/math] операций этого типа. Сложность операции для каждой формы - [math]\displaystyle{ O(d^3 + nd^2) }[/math], где за вклад в сложность порядка [math]\displaystyle{ O(d^3) }[/math] отвечают сингулярное разложение в [math]\displaystyle{ CPA(X, Y) }[/math] и вычисление матрицы поворота, а за сложность [math]\displaystyle{ O(nd^2) }[/math] отвечают матричные умножения в рамках алгоритма (см. раздел "Вычислительное ядро алгоритма").
• Вычисление консенсусной формы (сложение матриц, отвечающих [math]\displaystyle{ m }[/math] формам, и умножение этой суммы на скаляр): 1 операция на 1 итерацию GPA. Сложность - [math]\displaystyle{ O(mnd) }[/math].
• Масштабирование форм (вычисление скаляров для каждой формы): [math]\displaystyle{ m }[/math] операций на 1 итерацию GPA, если используется вариант GPA с учетом масштаба (иначе таких операций нет) (на нулевой итерации скаляры принимаются равными одной и той же константе). Сложность - [math]\displaystyle{ O(nd) }[/math].
• Пересчет функции ошибки: 1 операция на 1 итерацию GPA. Сложность - [math]\displaystyle{ O(nd) }[/math] (требуется вычисление нормы консенсуса).

А также макрооперации для нулевой итерации в частности:

• Вычисление первоначальной нормировочной константы: последовательная сложность [math]\displaystyle{ O(mnd) }[/math].
• Первоначальные нормировка и центровка текущих форм: сложность [math]\displaystyle{ O(nd) }[/math] для нормировки каждой формы.
• Инициализация показателей масштаба: сложность [math]\displaystyle{ O(1) }[/math] для каждого при известной константе нормировки.

Представляется разумным отдельно описать макрооперации более низкого уровня, составляющие классический прокрустов анализ [math]\displaystyle{ CPA(X, Y) }[/math] и обновление форм [math]\displaystyle{ X_i' }[/math]:

• Сингулярное разложение матрицы из пространства [math]\displaystyle{ \mathbb{R}^{d \times d} }[/math].
• Перемножение матриц вида [math]\displaystyle{ X^TY }[/math], где [math]\displaystyle{ X, Y\in \mathbb{R}^{n \times d} }[/math].
• Перемножение матриц вида [math]\displaystyle{ XR }[/math], где [math]\displaystyle{ X \in \mathbb{R}^{n \times d}, R \in \mathbb{R}^{d \times d} }[/math].
• Перемножение матриц вида [math]\displaystyle{ UV^T }[/math], где [math]\displaystyle{ U, V \in \mathbb{R}^{d \times d} }[/math].

Схема реализации последовательного алгоритма

Приведем схему реализации с использованием псевдокода.

Схема реализации процедуры классического прокрустова преобразования

function CPA(X, Y, scale=True, translation=True):
    # Compute column-wise mean if translation is allowed
    if translation:
        X_mean = mean_columns(X)
        Y_mean = mean_columns(Y)
        X0 = X - X_mean
        Y0 = Y - Y_mean
    else:
        X0 = X
        Y0 = Y
        X_mean = 0
        Y_mean = 0

    # Compute rotation using SVD
    U, Σ, Vt = SVD(transpose(X0) * Y0)
    R = U * Vt
    if determinant(R) < 0:
        Vt[last_row] *= -1
        R = U * Vt

    # Compute scale if needed
    s = 1
    if scale:
        s = sum_elements(elementwise_multiply(X0 * R, Y0)) / Frobenius_norm_squared(X0)

    # Align X
    aligned_X = s * X0 * R + Y_mean
    if translation:
        translation_vec = Y_mean - s * X_mean * R
    else:
        translation_vec = 0

    return aligned_X, R, s, translation_vec

Схема реализации GPA

function GPA(X_list, scale=True, tolerance=1e-4):
    m = length(X_list)

    # Center matrices and normalize
    for i in 0..m-1:
        X_list[i] = X_list[i] - mean_columns(X_list[i])

    var_lambda = sqrt(m / sum(Frobenius_norm_squared(X) for X in X_list))
    for i in 0..m-1:
        X_list[i] *= var_lambda

    scaling_factors = [var_lambda for i in 0..m-1]

    # Initial consensus
    Y_consensus = X_list[0]
    for i in 1..m-1:
        X_list[i] = CPA(X_list[i], Y_consensus, scale=False, translation=False)[0]

    Y_consensus = average(X_list)
    Sr = m * (1 - Frobenius_norm_squared(Y_consensus))

    # Iterative alignment
    loop:
        for i in 0..m-1:
            X_list[i] = CPA(X_list[i], Y_consensus, scale=False, translation=False)[0]

        if scale:
            for i in 0..m-1:
                change_factor = sqrt(sum_elements(elementwise_multiply(X_list[i], Y_consensus))) /
                                     (Frobenius_norm_squared(X_list[i]) * Frobenius_norm_squared(Y_consensus)))
                s_i_star = scaling_factors[i] * change_factor
                X_list[i] *= s_i_star / scaling_factors[i]
                scaling_factors[i] = s_i_star

        Y_consensus_new = average(X_list)
        Sr_new = Sr - m * (Frobenius_norm_squared(Y_consensus_new) - Frobenius_norm_squared(Y_consensus))

        if Sr - Sr_new < tolerance:
            break

        Y_consensus = Y_consensus_new
        Sr = Sr_new

    return Y_consensus, X_list, scaling_factors

Последовательная сложность алгоритма

Как было отмечено ранее в разделе "Макроструктура алгоритма", на каждой итерации GPA выполняется [math]\displaystyle{ m }[/math] (или [math]\displaystyle{ m-1 }[/math] на нулевой итерации) операций [math]\displaystyle{ CPA(X, Y) }[/math], причем на них приходится основная сложность. Сложность каждой из операций типа [math]\displaystyle{ CPA(X, Y) }[/math] равна [math]\displaystyle{ O(d^3 + nd^2) }[/math]. Так как на одной итерации GPA выполняется последовательное приведение сразу [math]\displaystyle{ m }[/math] (или [math]\displaystyle{ m - 1 }[/math] на нулевой итерации) форм к консенсусу с помощью классического прокрустова преобразования, то последовательная сложность одной итерации GPA целиком есть [math]\displaystyle{ O(m(d^3 + nd^2)) }[/math]. Число итераций [math]\displaystyle{ k }[/math] заранее неизвестно, так как алгоритм итеративный. Однако, с учетом этого числа сложность всего GPA составит [math]\displaystyle{ O(km(d^3 + nd^2)) }[/math].

Сложность по памяти составит:

• [math]\displaystyle{ O(mnd) }[/math] для хранения входных и текущих матриц
• [math]\displaystyle{ O(nd + d^2) }[/math] для хранения матрицы поворота для одной формы, временных матриц и векторов, и учитывая сингулярное разложение

Так как любые временные матрицы размера [math]\displaystyle{ d \times d }[/math] хранятся лишь на протяжении обработки одной из форм, а формы обрабатываются последовательно, то суммарная сложность последовательного GPA по памяти равна [math]\displaystyle{ O(mnd + d^2) }[/math].

Информационный граф

GPA имеет свойство быстро сходиться (всего за несколько итераций), особенно при достаточно большом значении допустимой ошибки [math]\displaystyle{ \varepsilon }[/math]. Поэтому для наглядности далее приведен как граф для нулевой итерации, так и граф для последующих итераций GPA. Также отдельно рассматривается информационный граф для процедуры классического прокрустова анализа без сдвигов и масштабирования - именно такой вариант используется в процедуре GPA.

Рассмотрим информационные графы начальной и последующих итераций GPA (Рис.2 и Рис.3):

Первая группа вершин (ICtr) отвечает первоначальной центровке форм.

Вторая группа вершин (λ) отвечает вычислению константы нормализации на начальной итерации.

Третья группа вершин (s_init) отвечает инициализации показателей масштаба.

Четвертая группа вершин (Norm) отвечает нормировке форм с учетом вычисленной константы нормализации.

Пятая группа вершин (CPA) отвечает процедуре классического прокрустова преобразования.

Шестая группа вершин (Y) отвечает вычислению консенсуса на основе текущих форм.

Седьмая группа вершин (Loss) отвечает пересчету функции ошибки.

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

Девятая группа вершин (Halt) отвечает проверке критерия останова и возврату из процедуры GPA при его выполнении.

Рис. 2. Информационный граф начальной итерации GPA

Рис. 3. Информационный граф последующих итераций GPA

Рассмотрим информационный граф процедуры CPA (Рис.4):

Первая группа вершин (MatMul_dnd) отвечает перемножению двух матриц размерностей [math]\displaystyle{ d \times n }[/math] и [math]\displaystyle{ n \times d }[/math] соответственно.

Вторая группа вершин (SVD) отвечает сингулярному разложению матрицы размерности [math]\displaystyle{ d \times d }[/math].

Третья группа вершин (MatMul_ddd_detchk) отвечает перемножению двух матриц размерностей [math]\displaystyle{ d \times d }[/math] каждая и коррекции результата в случае наличия отражений.

Четвертая группа вершин (align) отвечает приведению выравниваемой формы, то есть умножению двух матриц размерностей [math]\displaystyle{ n \times d }[/math] и [math]\displaystyle{ d \times d }[/math] соответственно.

Рис. 4. Информационный граф процедуры CPA

Также для наглядности в рассмотренные графы добавлены следующие вершины:

in_X - вход, матрица размера [math]\displaystyle{ n \times d }[/math], отвечающая одной текущей форме;

in_L - вход, вещественное число, отвечающее значению функции ошибки;

in_Y - вход, матрица размера [math]\displaystyle{ n \times d }[/math], отвечающая консенсусной форме;

in_s - вход, вещественное число, отвечающее одному показателю масштаба;

Аналогичным образом помечены выходы.

Ресурс параллелизма алгоритма

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

Сложность SVD заслуживает отдельного внимания, так как есть различные алгоритмы и на практике

Сначала проанализируем параллельную сложность CPA. Отметим, что матричные умножения и высчитывание определителя матрицы поворота в CPA вносят вклад [math]\displaystyle{ O(\log d + \log n + d) = O(\log n + d) }[/math] в высоту ЯПФ и [math]\displaystyle{ O(d^2 + nd) }[/math] в ширину. Если принять сложность SVD (по высоте ЯПФ) за [math]\displaystyle{ O(d^2) }[/math], тогда суммарный вклад CPA в высоту ЯПФ - [math]\displaystyle{ O(\log n + d^2) }[/math]. Таким образом, SVD является определяющим фактором высоты ЯПФ в процедуре CPA, а матричные умножения определяют ширину.

Далее, рассмотрим нулевую итерацию GPA. Центровка каждой формы отдельно вносит вклад [math]\displaystyle{ O(\log n) }[/math] в высоту ЯПФ и [math]\displaystyle{ O(nd) }[/math] в ширину. Итого, будучи распараллеленной, центровка дает вклад [math]\displaystyle{ O(mnd) }[/math] в ширину (высота не меняется при учете параллельных вычислений для различных форм независимо). Вычисление изначального показателя нормировки занимает [math]\displaystyle{ O(\log(mnd)) }[/math] по высоте ЯПФ и [math]\displaystyle{ O(mnd) }[/math] по ширине. Далее нормировка производится за [math]\displaystyle{ O(1) }[/math] по высоте и [math]\displaystyle{ O(mnd) }[/math] по ширине ЯПФ. Вклад инициализации показателей масштаба пренебрежимо мал относительно упомянутых операций. Далее на нулевой итерации выполняется приведение форм с помощью CPA. Для [math]\displaystyle{ m }[/math] форм параллельное приведение их к консенсусу дает вклад [math]\displaystyle{ O(m(d^2 + nd)) }[/math] по ширине ЯПФ, а по высоте - столько же, сколько вносит 1 отдельная процедура CPA ([math]\displaystyle{ O(\log n + d^2) }[/math]). Вычисление консенсуса происходит за [math]\displaystyle{ O(\log m) }[/math] шагов при одновременном выполнении [math]\displaystyle{ O(mnd) }[/math] операций. Вычисление функции потерь не увеличивает высоту и ширину ЯПФ. Итого, ширина ЯПФ для нулевой итерации GPA составляет [math]\displaystyle{ O(mnd + md^2) }[/math], а высота ЯПФ равна [math]\displaystyle{ O(\log(mnd) + d^2) }[/math].

Для последующих итераций GPA также имеет место выполнение [math]\displaystyle{ m }[/math] процедур CPA одновременно (параллельно), что вносит такой же вклад в сложность, как и выполнение данных процедур на нулевой итерации. Далее, вычисление показателей масштаба происходит за [math]\displaystyle{ O(\log(nd)) }[/math] по высоте и [math]\displaystyle{ O(mnd) }[/math] по ширине ЯПФ (параллельно для [math]\displaystyle{ m }[/math] форм). Выполняемое далее обновление форм не вносит дополнительной сложности (по асимптотике). Далее, как и на нулевой итерации, происходит вычисление консенуса и пересчет функции ошибки, что вносит вклад [math]\displaystyle{ O(\log m) }[/math] в высоту ЯПФ и [math]\displaystyle{ O(mnd) }[/math] в ширину. Проверка критерия останова также несущественна по сложности. Итого, последующие итерации вносят вклад [math]\displaystyle{ O(\log(nd) + \log(m) + d^2) = O(\log(mnd) + d^2) }[/math] в высоту ЯПФ и [math]\displaystyle{ O(mnd + md^2) }[/math] в ширину.

Итого, если число итераций GPA равно [math]\displaystyle{ k }[/math] (также считаем, что оно не является константой и может неограниченно расти), то можно ожидать сложность в терминах ЯПФ следующего порядка: [math]\displaystyle{ O(k(\log(mnd) + d^2)) }[/math] по высоте и [math]\displaystyle{ O(mnd + md^2) }[/math] по ширине. Высота ЯПФ в значительной степени определяется параллельной сложностью (высотой ЯПФ) сингулярного разложения. Также можно отметить компромисс между высотой и шириной ЯПФ, на который может влиять степень распараллеливания вычисления сумм.

Входные и выходные данные алгоритма

Входные данные:

• [math]\displaystyle{ (X_1, X_2,\dots,X_m) }[/math] - набор входных форм, то есть вещественных матриц одного размера [math]\displaystyle{ n \times d }[/math];
• флаг, отвечающий за включение опции вычисления факторов масштаба (булева переменная);
• [math]\displaystyle{ \varepsilon }[/math] - параметр точности.

Объем входных данных: [math]\displaystyle{ mnd + 2 }[/math].

Выходные данные:

• [math]\displaystyle{ Y }[/math] - форма консенсуса (матрица размерности [math]\displaystyle{ n \times d }[/math]);
• [math]\displaystyle{ (X'_1, X'_2,\dots,X'_m) }[/math] - набор текущих форм, то есть вещественных матриц одного размера [math]\displaystyle{ n \times d }[/math], выровненных относительно консенсуса;
• [math]\displaystyle{ (s_1, s_2,\dots,s_m) }[/math] - факторы масштаба (вещественные числа).

Объем выходных данных: [math]\displaystyle{ mnd + nd + m }[/math].

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

Важно заметить, что на практике в задачах, решаемых с помощью этого алгоритма, исследователи часто имеют дело с облаками точек в пространствах размерности 2 (например, в специальных задачах NLP и задачах CV в медицине) или 3 (также в задачах CV). Поэтому в применении к приложениям может быть важна скорее эффективная реализация перемножения матриц, нежели тщательная оптимизация SVD (действительно, в пространствах размерности 2 и 3 в процедуре CPA будет производиться разложение лишь матриц размерностей [math]\displaystyle{ 2 \times 2 }[/math] и [math]\displaystyle{ 3 \times 3 }[/math] соответственно. Также в таком случае решающее значение в оценке сложности алгоритма будет иметь исключительно то, насколько распараллеливаются матричные умножение и сложение.

Обычно число точек [math]\displaystyle{ n }[/math] сильно больше числа форм [math]\displaystyle{ m }[/math] и размерности пространства [math]\displaystyle{ d }[/math].

Также, как было упомянуто ранее, число итераций [math]\displaystyle{ k }[/math] обычно мало по сравнению с остальными параметрами.

Если принять [math]\displaystyle{ k, d }[/math] константами, то сложность последовательной реализации есть [math]\displaystyle{ O(mn) }[/math], в то время как теоретическая сложность (высота ЯПФ) параллельной реализации будет [math]\displaystyle{ O(log(mn)) }[/math]. Очевиден выигрыш по времени.

Если [math]\displaystyle{ k, d }[/math] фиксированы, то процедура GPA выполняет константное число операций ([math]\displaystyle{ O(1) }[/math]) на единицу входных (и выходных) данных.

Программная реализация алгоритма

Масштабируемость алгоритма и его реализации

Исследование масштабируемости параллельной реализации алгоритма проводилось на суперкомпьютере «Ломоносов-2». Для этого была исследована авторская реализация алгоритма с использованием средств MPI и OpenMP. В ходе исследования было произведено множество запусков программы с различными значениями параметра, задающего размерность форм и их число, а также различным числом процессов.

Далее представлены два графика зависимости времени выполнения от объема входных данных и числа процессов в различных масштабах. Важно, что здесь для наглядности представления положено [math]\displaystyle{ m = n = d }[/math], т.е. число входных форм совпадает с числом строк и столбцов в этих формах. Соответственно, под размерностью матрицы понимается число [math]\displaystyle{ n=d }[/math]. Количество итераций равно двум (не включая нулевую).

• Размерность входа: от 2 включительно до 202 включительно с шагом 10.
• Число процессов: [1, 2, 4, 8, 16, 32, 64, 92].

Рис.5. График зависимости времени выполнения нескольких алгоритма от объема входных данных и числа процессов (увеличенный масштаб)

Рис.6. График зависимости времени выполнения нескольких алгоритма от объема входных данных и числа процессов

Из графиков видно, что при большом объеме входных данных имеется выигрыш по времени выполнения при росте числа процессоров, в то время как при небольших матрицах и небольшом числе форм выделение дополнительных процессоров лишь увеличивает накладные расходы, вследствие чего время выполнения растет. Также очевидна закономерность, при которой время выполнения в среднем растет с увеличением объема входных данных при неизменном числе процессоров. В условиях исследуемого алгоритма, конфигурации вычислительных узлов суперкомпьютера (многоядерные процессоры с подержкой гипертрединга) и реализации (используется библиотека Eigen с установленным ограничением на распараллеливание с помощью OpenMP в 3 нити - это значение было получено опытным путем как оптимальное) число процессоров можно смоделировать числом порожденных процессов MPI.

Эффективность распараллеливания есть величина [math]\displaystyle{ E = S/p }[/math], где [math]\displaystyle{ S = \frac{T_1}{T_p} }[/math] - ускорение, [math]\displaystyle{ p }[/math] - число используемых процессоров, и [math]\displaystyle{ T_i }[/math] - время работы программы на [math]\displaystyle{ i }[/math] процессорах.

Рис.7. График зависимости эффективности распараллеливания от объема входных данных и числа процессов

Рис.8. График зависимости ускорения от объема входных данных и числа процессов

По графикам можно наблюдать сверхлинейное ускорение в случае относительно небольшого числа процессоров. Это может быть связано с активным использованием кэша. Так как GPA в каждом процессе работает только с локальным набором форм (может быть небольшим по объему), то формы могут быть сохранены в кэше, а не быть загруженными из оперативной памяти и затем быть снова выгруженными, причем по несколько раз (что имеет место в последовательной реализации).

Максимальная эффективность составила [math]\displaystyle{ 611.8\% }[/math] в случае [math]\displaystyle{ n=202 }[/math], и при числе процессов [math]\displaystyle{ p=4 }[/math].

Минимальная эффективность составила [math]\displaystyle{ 0.0022\% }[/math] в случае [math]\displaystyle{ n=12 }[/math], и при числе процессов [math]\displaystyle{ p=92 }[/math].

Исходный код, описание параметров окружения и данные в числовом виде доступны в GitHub-репозитории.