Slow Feature Analysis. Разбор метода и реализация на Python с нуля

В этой статье я хочу рассказать про метод обучения без учителя – “Анализ медленных признаков” (Slow Feature Analysis), далее SFA. Метод был разработан в 2002 году Лоренцом Вискоттом и Терренсом Сейновски.

SFA можно использовать для выделения стабильных сигналов на фоне шума, такие как отслеживание объектов на видео, трендов цен из финансовых данных, признаков износа по вибрациям оборудования.

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

Объяснение алгоритма и реализация его на Python

Пусть у нас есть многомерный временной сигнал длины T и количеством признаков N такой, что:

Где значения временного сигнала, а t = 1, 2, …, T.

Для этой задачи будет создан искусственный пример – три различных сигнала:

• Синус с белым шумом

• Синус с угасающим розовым шумом

• Синус с экстремальными значениями и броуновским блужданием.

Примечание: наличие розового шума, броуновского блуждания и выбросов сделано намеренно, чтобы подсветить границы применимости метода. Как будет показано далее на графиках, SFA успешно выделяет общий тренд, но его весовые коэффициенты для зашумленных каналов стремятся к нулю, эффективно отключая “сломанный сигнал”. Это показывает как силу, так и ограничение метода – он не восстанавливает сигнал из шума, а стремится игнорировать источник шума.

В примере будут использоваться исходные три столбца, но для более сложных задач можно создавать степенные, логарифмические и более сложные признаки и их композиции.

import numpy as np
import pandas as pd

np.random.seed(1)
pd.set_option('display.float_format', '{:.4f}'.format)

data_size = 500
# 4 периода синуса
x = np.linspace(0, 4 * np.pi, data_size) 

sin_wave = np.sin(x)
# Создание синуса с белым шумом
white_noise = np.random.normal(0, 0.2, data_size)
series_1 = sin_wave + white_noise

# Создание синуса угасающим розовым шумом
pink_noise = (lambda x: x / np.std(x))(
    np.real(
        np.fft.ifft(
            np.concatenate([
                [0],
                (1 / np.sqrt(np.abs(np.fft.fftfreq(data_size)[1:]) + 1e-10))
                * np.exp(2j * np.pi * np.random.random(data_size - 1))
            ])
        )
    )
)

decay_envelope = 2.25 * np.exp(-2.25 * x / x.max()) + 0.25

pink_noise = pink_noise * decay_envelope

series_2 = sin_wave + 0.4 * pink_noise

# Создание синуса с выбросами и броуновским движением
brownian = np.cumsum(np.random.normal(0, 0.5, data_size))
extreme_values = np.zeros(data_size)
outlier_indices = np.random.choice(data_size, size=int(0.5 * data_size), replace=False)
extreme_values[outlier_indices] = np.random.uniform(-1, 1, size=len(outlier_indices))
series_3 = sin_wave + extreme_values + 0.4 * brownian

df = pd.DataFrame({
    'series_1': series_1,
    'series_2': series_2,
    'series_3': series_3
})
  

Мы хотим найти такие функции , чтобы выходные сигналы изменялись наиболее медленно. Это достигается минимизацией следующего функционала:

Где производная по времени, а треугольные скобки ⟨⟩ означают усреднение по времени. Так как у нас дискретный случай, то производная по времени будет считаться следующим образом:

Для статьи был выбран самый простой способ, но при необходимости можно использовать и центральные разности:

Как и у многих оптимизационных задач здесь есть ограничения:

1. Нулевое среднее. Это нужно чтобы исключить неинформативные константные признаки.

2. Единичная дисперсия. Все признаки нужно привести к одному масштабу.

3. Некоррелированность. Это нужно чтобы каждый следующий признак не дублировал информацию.

В настоящий момент данные образуют следующее облако точек:

Чтобы наши данные удовлетворяли этим требованиям мы сделаем следующие действия:

Для центрирования следует вычесть из каждого элемента среднее:

df_centr = pd.DataFrame()

for column in df.columns:
    df_centr[column] = df[column] - df[column].mean()
  

Вычисляем ковариационную матрицу. Ковариационная матрица признаков определяется как математическое ожидание внешнего произведения центрированного вектора этих признаков, то получается следующее:

A = (df_centr.T @ df_centr) / (len(df_centr) - 1)
  

A = np.cov(df_centr.T)
  

Теперь нужно сделать сферизацию, ZCA whitening или, как это еще называют, отбеливание Махаланобиса. Для этого нужно выполнить сингулярное разложение ковариационной матрицы. Получится следующее: , важно заметить, что для наших целей элементы матрицы D будут в степени . Полученная формула будет такой:

Далее нам нужно умножить:

Матрица A_s выполняет преобразование, которое приводит данные к единичной ковариационной матрице , где I – единичная матрица.

# Выполняем сингулярное разложение матрицы A
U, D, Vt = np.linalg.svd(A)
# Создаём диагональную матрицу из обратных квадратных корней сингулярных значений
D_inv_sqrt = np.diag(D ** (-0.5))  
# Умножение U @ D_inv_sqrt @ U.T даёт матрицу, которая будет использоваться для сферирования данных
A_s = U @ D_inv_sqrt @ U.T
# Применяем сферирование к df
df_sphered = df_centr @ A_s
  

После отбеливания облако точек выглядит так:

Вычисляем производные сферированных данных:

Примечание: усреднение ⟨⟩ идет по t от 1 до T-1, так как для последнего элемента T нельзя вычислить по формуле .

df_derivative = df_sphered.diff().dropna().reset_index(drop=True)
  

Центрируем полученные производные:

for column in df_derivative.columns:
    df_derivative[column] = df_derivative[column] - df_derivative[column].mean()
  

Создаем матрицу ковариации B из полученных данных.

B = (df_derivative.T @ df_derivative) / (len(df_derivative) - 1)
  

B = np.cov(df_derivative.T)
  

Вспомним, что мы хотели минимизировать среднюю квадратичную скорость изменения проекции:

А есть ни что иное как матрица ковариации . Поэтому:

Примечание: для общности вывода мы сначала рассмотрим случай без сферизации, где ограничение записывается как . Затем подставим как частный случай.

Для поиска решений используем метод множителей Лагранжа.
Чтобы найти значение j-того признака запись формулы будет иметь следующий вид:

Берем производную по w и приравниваем к нулю:

Умножаем слева на , где k – индекс любого уже найденного решения (k < j).

У нас получается:

Из-за условия ортогональности и .

Так как при мы получаем . Из этого у нас следует то, что и мы можем переписать уравнение в следующем виде:

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

Так как ранее предварительно сферировали данные, их ковариационная матрица стала единичной – , поэтому получается , что можем записать еще более просто: .

Где это собственное значение, чем оно меньше, тем более медленный признак. Поэтому сортируем по возрастанию значения < <…<

В результате вычисления собственных значений матрицы получаем:

Здесь видно, что первая компонента изменяется практически в 6 раз медленнее третьей.

Признаки получатся из решения где самый медленный признак, следующий и так далее. Поэтому решение будет выглядеть как:

Где матрица собственных весов.

# Вычисляем собственные значения и собственные векторы марицы B
lambdas, w = np.linalg.eigh(B)
# Получаем индексы для сортировки собственных значений
sort_idx = np.argsort(lambdas)
# Сортируем собственные значения
lambdas_sorted = lambdas[sort_idx]
# Сортируем собственные векторы
w_sorted = w[:, sort_idx]
# Эвристика, которая нужна для воспроизводимости, так как иногда
 из-за особенностей вычисления график будет отражаться. 
# Например библиотека sksfa может вернуть перевёрнутый график 
# по сравнению с ручной реализацией. 
# Умножение на -1 исправит его положение.
for i in range(w_sorted.shape[1]):
    if np.sum(w_sorted[:, i]) < 0:
        w_sorted[:, i] *= -1
# Умножаем матрицу данных на матрицу отсортированных 
# собственных векторов, получая медленные признаки.
y = df_sphered @ w_sorted
  

Вот как выглядят графики медленных признаков:

Рассмотрим вклад исходных признаков в полученные решения:

w_original = A_s @ w_sorted

weights_df = pd.DataFrame(
    w_original,
    index=df_sphered.columns,
    columns=[f'Веса для медленного признака {i+1}' for i in range(w_original.shape[1])]
)

print("Вклад исходных признаков в каждую медленную компоненту:")
display(weights_df)
  

Как видно в таблице, модель присваивает в самой медленной компоненте наибольший вес первому, самому стабильному признаку. Мы можем изменять параметры и полученные результаты подтвердят то, что SFA склонен отказываться от наиболее шумных источников.

Реализация алгоритма с помощью библиотеки sksfa

В интернете есть реализации этого алгоритма. Я рассмотрю библиотеку sksfa для описанной задачи:

Код с использованием sksfa

from sksfa import SFA
import numpy as np
import pandas as pd

sfa = SFA(n_components=3)

y_sksfa = sfa.fit_transform(df.values)


w, b = sfa.affine_parameters()

# Чтобы график компонент не был перевернутым, 
# а веса с противоположным знаком, умножу w и y_sksfa на -1
w[0, :] *= -1
w[2, :] *= -1
y_sksfa[:,0] *= -1
y_sksfa[:,2] *= -1

weights_df = pd.DataFrame(
    w.T,
    index=df.columns,
    columns=[f'Веса для медленного признака {i+1}' for i in range(w.shape[1])]
)

print("Вклад исходных признаков в каждую медленную компоненту (sksfa):")
display(weights_df)
  

Собственные значения, полученные через sksfa, совпадают с ручным расчётом:

Сравнение показало, что ручная реализация SFA и библиотека sksfa работают идентично. Собственные значения тоже совпадают, а после приведения знаков весов и компонент к общему виду совпадают веса и сами медленные признаки. Так что вы можете использовать любое удобное решение в своих проектах.

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

Полезные ссылки:

1. Оригинальная статья: Wiskott, L. and Sejnowski, T. J. (2002). Slow feature analysis: Unsupervised learning of invariances. Neural computation, 14(4):715-770.

2. sksfa – implementation of Slow Feature Analysis for scikit-learn.

3. Код в Jupyter Notebook – юпитер ноутбук с кодом из статьи.