Slow feature analysis consists of three main steps
Once the proper connection is made to PCA, this method becomes easy to understand and implement.
Polynomial expansion is often applied as a preprocessing step for SFA. All papers published on the topic avoid to write down explicitely the intended meaning behind nonlinear expansion. One reason for this might be that due to the wide use of the SFA-tk, rSFA or MDP packages for computing SFA which hide much of the details. Berkes (2005) is more precise on this question. For a random vector of size 3 for example, its corresponding polynomial expansion of degree 2 becomes:
\begin{equation} \mathbf{h}(\mathbf{x}) = [x_1^2, x_1x_2, x_1x_3, x_2^2, x_2x_3, x_3^2, x_1, x_2, x_3] \end{equation}That is, $\mathbf{h}$ consists of all the monomials up to degree 2. Using the built-in product
function, we can generate all possible n-tuple for which the sum is less than or equal to the specified degree.
from itertools import product
monomials = lambda n, d: [t for t in product(range(d+1), repeat=n) if sum(t)<=d and sum(t) > 0]
expand = lambda X, d: np.apply_along_axis(lambda x, mn: np.prod(np.power(x, mn), axis=1), 1, X, monomials(X.shape[1], d))
X = np.array([[2, 3, 4], [5, 6, 7]])
Y = expand(X, 2)
print sorted(Y[0])
print sorted(Y[1])
[2, 3, 4, 4, 6, 8, 9, 12, 16] [5, 6, 7, 25, 30, 35, 36, 42, 49]
We can compare the above output with the Matlab expansion.m
function from sfa_tk
import pymatbridge as pymat
from pymatbridge import Matlab
pymat.load_ipython_extension(get_ipython())
mlab = Matlab()
mlab.start()
mlab.run_code('expansion([], [2 3 4; 5 6 7])')
Starting MATLAB on http://localhost:65268 visit http://localhost:65268/exit.m to shut down same ....MATLAB started and connected! Starting MATLAB on http://localhost:65284 visit http://localhost:65284/exit.m to shut down same ...MATLAB started and connected!
{u'content': {u'code': u'expansion([], [2 3 4; 5 6 7])', u'datadir': u'/private/tmp/MatlabData/', u'figures': [], u'stdout': u'\nans =\n\n 2 3 4 4 6 8 9 12 16\n 5 6 7 25 30 35 36 42 49\n\n'}, u'success': u'true'}
In order to satisfy the constraints in the optimization problem, we need to subtract the mean and transform the data so that it has unit variance: this can be done using PCA. Let's first generate some synthetic data to see the effect of whitening.
import numpy as np
from matplotlib import pyplot as plt
mu1 = np.array([0.0, 0.0])
sigma1 = np.array([[1.0, 0.5], [0.5, 1.0]])
X = np.random.multivariate_normal(mu1, sigma1, 1000)
plt.plot(X[:,0], X[:,1], '.')
[<matplotlib.lines.Line2D at 0x106abd550>]
We will then subtract the mean of $X$ and take its SVD decomposition \begin{equation} \mathbf{X} = \mathbf{U} \Sigma \mathbf{V} \end{equation}
The columns of $U$ are called left-singular vectors of $X$. They correspond to the eigenvectors of $XX^\top$ (covariance). We then apply the whitening transformation using dot product \begin{equation} \mathbf{X}_W = \mathbf{X}(\mathbf{V}\Sigma^{-1/2})^\top \end{equation}
def whiten(X):
U, S, V = np.linalg.svd(X, full_matrices=False)
U *= sqrt(X.shape[0])
return U
X = whiten(X)
plt.plot(X[:,0], X[:,1], '.')
[<matplotlib.lines.Line2D at 0x106af8690>]
We compute the temporal derivative of the sphered signal using a discrete approximation.
derivative = lambda X: X[1:, :]-X[:-1, :]
Using the temporal derivative of the whiten signal, we can now obtain a basis in terms of its slowly varying components. It suffices here to once more to take the SVD decomposition. Numpy returns the singular values in the standard descending order. Therefore the first singular value is associated with the direction in which the variance in the temporal manifold is the highest. Since slow feature analysis aims at recovering the slow components, we want to pick the left singular vectors associated with the smallest singular values first. The slow features are the projection of the whiten signal onto the singular components of least variance.
t = np.linspace(0, 2*np.pi, 5000)
x1 = np.sin(t)+2*np.power(np.cos(11.*t), 2)
x2 = np.cos(11.*t)
X = np.column_stack((x1, x2))
plt.subplot(1, 2, 1)
plt.plot(X[:,1], X[:,0], '.')
Y = whiten(expand(X, 2))
U, S, V = np.linalg.svd(derivative(Y), full_matrices=False)
plt.subplot(1, 2, 2)
Z = Y[:-1].dot(V[-1,:])
plt.plot(Z)
[<matplotlib.lines.Line2D at 0x1020a10d0>]
%%matlab
path(path, '/Users/pierrelucbacon/workspace/sfa_tk/sfa');
path(path, '/Users/pierrelucbacon/workspace/sfa_tk/lcov');
T = 5000;
t = linspace(0, 2*pi, T);
x1 = sin(t)+cos(11*t).^2;
x2 = cos(11*t);
x = [x1; x2]';
% plot the input signal
subplot(1,3,1); plot(x2,x1);
set(gca, 'Xlim', [-1.2, 1.2], 'YLim', [-1.2,2.2]);
title('input signal x(t)');
xlabel('x_2(t)'); ylabel('x_1(t)');
[y, hdl] = sfa2(x);
subplot(1,3,2); plot(t, y(:,1));
set(gca, 'Xlim', [0,2*pi], 'PlotBoxAspectRatio', [1,1,1])
title('output of the slowest varying function');
xlabel('t'); ylabel('y_1(t)');
init preprocessing close preprocessing whitening and dimensionality reduction {Warning: PACK can only be used from the MATLAB command line.} > In lcov_clear at 11 In sfa2_step at 41 In sfa_step at 34 In sfa2 at 27 In web_eval at 44 In webserver at 241 init expansion step close expansion step perform slow feature analysis {Warning: PACK can only be used from the MATLAB command line.} > In lcov_clear at 11 In sfa2_step at 96 In sfa_step at 34 In sfa2 at 29 In web_eval at 44 In webserver at 241 {Warning: PACK can only be used from the MATLAB command line.} > In lcov_clear at 11 In sfa2_step at 98 In sfa_step at 34 In sfa2 at 29 In web_eval at 44 In webserver at 241 SFA2 closed