Участник: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] (включает перемножение матрицы размера [math]\displaystyle{ n \times d }[/math] на матрицу размера [math]\displaystyle{ d \times n }[/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{ VU^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
    translation_vec = Y_mean - s * X_mean * R

    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.