PFAx: Predictable Feature Analysis to Perform Control

arXiv:1712.00634v1 [cs.LG] 2 Dec 2017

Stefan Richthofer∗, Laurenz Wiskott†
Institut für Neuroinformatik,
Ruhr-Universität Bochum, Germany
December 5, 2017

Abstract
Predictable Feature Analysis (PFA) [Richthofer and Wiskott, 2015] is an algorithm that
performs dimensionality reduction on high dimensional input signal. It extracts those subsignals that are most predictable according to a certain prediction model. We refer to these
extracted signals as predictable features.
In this work we extend the notion of PFA to take supplementary information into account for improving its predictions. Such information can be a multidimensional signal like
the main input to PFA, but is regarded external. That means it won’t participate in the
feature extraction – no features get extracted or composed of it. Features will be exclusively
extracted from the main input such that they are most predictable based on themselves and
the supplementary information. We refer to this enhanced PFA as PFAx (PFA extended).
Even more important than improving prediction quality is to observe the effect of supplementary information on feature selection. PFAx transparently provides insight how the
supplementary information adds to prediction quality and whether it is valuable at all. Finally we show how to invert that relation and can generate the supplementary information
such that it would yield a certain desired outcome of the main signal.
We apply this to a setting inspired by reinforcement learning and let the algorithm learn
how to control an agent in an environment. With this method it is feasible to locally optimize
the agent’s state, i.e. reach a certain goal that is near enough. We are preparing a follow-up
paper that extends this method such that also global optimization is feasible.
∗ Electronic
† Electronic

address: stefan.richthofer@ini.rub.de; Corresponding author
address: laurenz.wiskott@ini.rub.de

1

1

Introduction

The reinforcement learning (RL) setting consists of an agent in an environment and usually involves
three signals – perception, action and reward. The agent has access to the perception and reward
signal and can generate the action signal as an output, which in turn yields consequences for future
perception and reward. The agent’s goal is to maximize reward over a certain time period.
In this work we focus on the action/perception cycle (figure 1) of RL scenarios with continuous
state space, perception signal and action signal. We – however – do not consider an arbitrary
reward signal for now, but just aim for reaching specific states indicated by their corresponding perception. A core issue of reinforcement learning is to represent the agent’s state in a way

action

agent

environment

perception
Figure 1: Perception/control cycle
that is of tractable complexity and at the same time allows to perform the desired task. In
[Richthofer and Wiskott, 2015] we have presented Predictable Feature Analysis (PFA) as an algorithm to focus on features of the environment that are intrinsically predictable. The idea is that
predictability is a crucial property for any task – without it, one cannot estimate consequences
of possible actions and would have to act randomly. There are – however – different notions of
predictability, e.g. differing in the choice of model, whether a model is used at all, what to predict
and what a prediction can be based on. These factors can crucially make a difference in how useful
predictability really is for an RL scenario. The original PFA approach selects features that are
most suitable to predict themselves.
We go one step further now and propose an extension to PFA that enables more fine-grained
control of what a prediction can be based on. This allows to incorporate the action or control
signal of RL settings and thus to aim for features that are not only predictable, but that can
be well manipulated. To accomplish this, we extend PFA to PFAx by enabling it to take any
supplementary information into account for prediction. Note that such supplementary information
will not be used for data extraction – PFAx will compose features exclusively from its main input
such that they are well predictable based on themselves and given the supplementary information.
While we have a concrete benefit for RL in mind, this enriches the approach and opens up various
new applications more generally.
Beneath improving prediction quality, this enhancement enables us to study the effect of supplementary information on feature selection. PFAx provides a coefficient matrix that transparently
shows how valuable the supplementary signal is for prediction of the extracted features. It turns
out that the prediction rule is invertible such that we can generate an action signal that is most
likely to result in a specific desired manipulation given the agent’s current state. This way we can
perform goal directed navigation in an RL environment.
2

1.1

Related work

A related SFA-inspired approach for RL is Contingent Feature Analysis (CFA) [Sprague, 2014].
The main difference to the PFAx based approach presented here is that CFA assumes a discrete
action set and uses a different selection criterion. PFAx assumes a continuous signal for control
– think of pressing buttons (discrete control) compared to turning a steering wheel or moving
levers (continuous control). CFA requires a no-operation (NOP) among the possible actions and
then selects features to behave as uncorrelated as possible to NOP induced behavior. Given that
continuous control is discretely sampled in digital processing, a comparison of PFAx and CFA
would be an interesting exercise in the future, but is not in the scope of this work.
There are a number of approaches with the same goal as described in
[Richthofer and Wiskott, 2015].
Forecastable Component Analysis (ForeCA) [Goerg, 2013], is a method based on the same
paradigm as PFA, but using a model-independent approach.
There are also differences
in scalability and we investigated a side-by-side comparison of PFA, ForeCA and SFA in
[Richthofer and Wiskott, 2015] to some extend.
Graph-based Predictable Feature Analysis (GPFA) [Weghenkel et al., 2017] analyses the distribution of data samples given previous samples. Low variance in that distribution is used to detect
predictability. The work draws links to graph embedding and an information-theoretic measure
of predictive information. It compares GPFA with SFA, PFA and ForeCA under various criteria,
including prediction quality and runtime performance.
Predictive Projections [Sprague, 2009] is an approach inspired by metric learning, more specifically by Neighborhood Components Analysis (NCA) [Goldberger et al., 2004]. It uses conjugate
gradient descent to find good projections for accurate prediction of future states in Markov Decision
Processes. While PFA is motivated by a comparable goal, it uses a different problem formulation
and optimization technique.
In [Box and Tiao, 1977] a comparable approach to PFA is formulated and PFA can be seen
as a generalization of the method presented there. We will refer to that method again in a later
section.

2

Extracting predictable features

We start with a comprehension of the PFA algorithm [Richthofer and Wiskott, 2015] and then
extend the method and notation to incorporate supplementary information. Given an input-signal
x(t) with n components, PFA’s objective is to find r most predictable output components, referred
to as “predictable features”.
The extraction itself is linear but can be enhanced by a non-linear expansion h. Expanding the input by monomials up to a desired degree, the extraction function essentially becomes a
polynomial of that degree1 . Note that by the Stone-Weierstrass theorem this technique can approximate any continuous function. This further covers regulated functions (i.e. piecewise continuous)
as these can be approximated by continuous functions. However, raising the expansion to high
degree might require a lot of training data and cause high computational costs. Applying PFA
1 For high-degree expansion, Legendre or Bernstein polynomials are preferable bases compared to monomials
because of better numerical stability.

3

in a hierarchical manner like is done with SFA in [Franzius et al., 2007a, Franzius et al., 2007b,
Schönfeld and Wiskott, 2015] can help to keep computation tractable.
As a prerequisite for PFA we shortly recall the SFA algorithm.

2.1

Recall SFA

Instead of predictability, SFA optimizes for slow variation. Extraction is performed by linear
transformation and projection. The extraction matrix is optimized over a finite training-phase
Ωt consisting of equidistant time points. Like mentioned above, a non-linear expansion h can be
applied to the signal. To avoid trivial constant solutions, output is constrained to have unit variance
and zero mean. Further more, the output signals must beP
pairwise uncorrelated to avoid redundant
output components. Mean is defined as hs(t)it := |Ω1t | t∈Ωt s(t) (average over training phase).
The initial step of SFA is to sphere the expanded signal over the training-phase, i.e. shift its mean
to zero and normalize the covariance-matrix to identity:
z̃(t)
z(t)

:=
:=

h(x(t)) − hh(x(t))it
− 21

with S := z̃z̃T

Sz̃(t)

(make mean-free)

(1)

(normalize covariance)

(2)

SFA then becomes the following optimization problem:
For i ∈ {1, . . . , r}
minimize
n
ai ∈R

subject to

aTi żżT ai
aTi hzi

=

0

(zero mean)

aTi zzT ai

=

1

(unit variance)

aTi zzT aj

=

0 ∀ j < i (pairwise uncorrelated)

(3)

Because of sphering it holds that hzi = 0 and zzT = I, transforming the constraints to
aTi aj = δij

(4)

For later equations a matrix notation of this constraint will be handy.
Ar := (a1 , . . . , ar ) ∈ Rn × r denoting the extraction matrix, (4) is equal to
∃A ∈ O(n) :

Ar

=

AIr

With
(5)

O(n) ⊂ Rn×n denotes the space of orthogonal transformations, i.e. AAT = I and Ir ∈ Rn×r
denotes the reduced identity matrix consisting of the first r Euclidean unit vectors as columns.
Choosing ai as eigenvectors of żżT , corresponding to the eigenvalues in ascending order, yields
Ar solving (3) globally. [Wiskott et al., 2011] describes this procedure in detail. In the following
m := ATr z denotes the extracted signal.

2.2

Modeling the PFA-problem

To measure predictability, we focus on linear, auto-regressive prediction as our model – it is
successfully used to model various time-related problems. In this model, a signal is predictable if it
4

can be approximated by a linear combination of p recent values. Formally, this yields the problem
of finding vectors a and b such that
aT z(t)

!

≈
=

b1 aT z(t − 1) + . . . + bp aT z(t − p)

(6)

T

a histz,p (t) b

(7)

with hist defined as the signal’s history of p time-steps:
:=

histz,p,∆ (t)

p
X

z(t − i∆)eTi

with ei ∈ Rp ,

(e1 , . . . , ep ) = Ip,p .

(8)

i=1

Ip,p denotes the p-dimensional identity and ei denotes the i-th p-dimensional Euclidean unit vector.
∆ defaults to 1: histz,p := histz,p,1 .

PFA

z(t)

aT z(t)

Figure 2: Illustration of PFA. Components are selected to be well
predictable based on themselves.
Like in SFA, we optimize parameters over Ωt and also adopt the constraints to avoid trivial
or redundant solutions. The first steps of PFA are indeed equal to those in SFA, i.e. non-linear
expansion and sphering. Where possible, we adopt notation from 2.1. Commonly, (6) is extended
to multiple dimensions as follows:
m(t)

!

≈

B1 m(t − 1) + . . . + Bp m(t − p)

with Bi ∈ Rn×n , diagonal

(9)

In [Richthofer and Wiskott, 2015] we explain why this restriction to diagonal Bi is not suitable for
PFA. We generalize it to:
m(t)

!

≈
=

B1 m(t − 1) + . . . + Bp m(t − p)
B vec(histz,p (t))

with Bi ∈ Rn×n

with B := (B1 , . . . , Bp ) ∈ R

n×np

(10)
(11)

In that formulation, each extracted component’s prediction can utilize all other extracted components. (9) and (10) are equal for n = 1. For (10) we can initially fit our data in full dimension
and search for the best-fitted components afterwards. This would not be possible for (9), because
fitting quality of each component is not invariant under the transformation used for extraction.
Nevertheless, we mention strategies to solve (9) in the appendix, section A.2.
Analytically, we obtain the following regression formula for an optimal B, given an orthogonal
extraction matrix Ar ∈ Rn×r :

−1
Bz (Ar ) := ATr zζ T Ar Ar T ζζ T Ar
∈ Rr×rp
(12)
| {z }
| {z }
n×np

np×np

5

In (12) we used ζ(t) := vec(histz,p (t)) ∈ Rnp and the following shortcut notation defined for any
matrix M ∈ Rn×m :
!
M

:=

M

Ip,p ⊗ M

0

..

=

∈ Rnp×mp

.
0
M
| {z }

(13)

p times M

Equation (12) can be derived as follows. For a given Ar the optimal B must solve
kATr z − BAr T ζk2

minimize
r×rp
B∈R

=: f (B)

(14)

We expand the value function f (B)
f (B)

=

ATr zzT Ar − ATr zζ T Ar BT − BAr T ζzT Ar + BAr T ζζ T Ar BT

(15)

and set its matrix derivative to zero
∂
f (B)
∂B

=

−2ATr zζ T Ar + 2BAr T ζζ T Ar

!

=

0

(16)

Solving (16) for B yields (12).
We define W to denote the optimal prediction matrix for the sphered signal z without any
extraction or transformation applied yet, i.e. set Ar = I in (12):
W

:=

Bz (I)

=

zζ T

ζζ T

−1

∈ Rn×np

(17)

It can happen that ζζ T is not (cleanly) invertible due to very small or zero-valued eigenvalues.
We regard it best practice to project away the eigenspaces corresponding to eigenvalues below a
critical threshold. These indicate redundancies in the signal and should not be used for prediction.
In an eigenvalue decomposition of ζζ T replace eigenvalues below the threshold by 0 and invert
the others. Use the resulting matrix as a proxy for ζζ T
Bz (A) = AT WA ∈ Rn×np .
If z is sphered, PFA yields the following problem:
minimize
A∈O(n)

−1

. For r = n and A ∈ O(n), we have

kATr z − Bz (Ar )Ar T ζk2

(18)

Projection composed with inversion renders (18) intractable by every method known to us (not
counting local or probabilistic methods). We propose the following tractable relaxation:
minimize
A∈O(n)

kATr z − ITr Bz (A)AT ζk2

=

kATr (z − Wζ)k2

(19)

Informally, (19) aims for components predictable based on the entire input, while (18) aims for
(0)
components predictable based on themselves. Let A∗r denote a global optimum of (18) and Ar a
global optimum of (19).
(19) is globally solvable by writing it as


T
minimize
Tr ATr (z − Wζ) (z − Wζ) Ar
(20)
A∈O(n)

6

z(t)

m(t) = ATr z(t)

PFA

Figure 3: Illustration of relaxation (19). Components are selected to
be well predictable based on the original input signal rather than on
themselves.
T

and choosing A such that it diagonalizes (z − Wζ) (z − Wζ)
and sorts the r smallest eigenvalues to the upper left. This can be described as performing PCA on the residuals of the
least squares fit. By some calculus, this can be shown to be equal to the method proposed in
(0)
[Box and Tiao, 1977]. To use Ar with (18), the prediction model must be refitted to the extracted
(0)
output by calculating Bz (Ar ) as defined in (12). In [Richthofer and Wiskott, 2015], we show that
the relaxation gap is related to the absolute prediction error of the optimal solution of (18). If
that error is zero, no relaxation gap exists at all. If the error is significant, the solution obtained
(0)
as Ar can suffer from overfitting being sub-optimal for (18). In [Richthofer and Wiskott, 2015],
we propose the following heuristic method to overcome this overfitting.

2.3

Avoiding overfitting

We propose the heuristics that signals well predictable in terms of (18) yield a lower errorpropagation to subsequent predictions than signals that are well predictable in terms of (19) but
not in terms of (18). The intuition is that in the second case, prediction is partly based on noisy
data – thus subsequent predictions inherit a higher error.
We define
−1
V :=
ζ(t + 1)ζ T (t) t ζζ T
∈ Rnp×np
(21)
V predicts ζ(t + 1) from ζ(t) like W predicts z(t) from ζ(t). Note that (21) is equal to (17)
with z(t) replaced by ζ(t + 1). The topmost n entries of ζ(t + 1) equal those of z(t) while the
others are down-shifted components of ζ(t). Indeed this shift operation is represented in V as we
can write it in terms of W:
W
V

=

I

0

The upper part consisting of W is responsible for predicting z(t), i.e. the topmost part of
ζ(t + 1) while the lower part performs the shift operation. Now we can perform iterated prediction
as follows:
z(t) ≈ WVi ζ(t − i)
(22)

7

m(t) = AT z(t)

PFA

z(t)

Figure 4: Illustration of iterated prediction. Components are selected to be
well predictable based on the original input signal and over several iterations
of prediction. This filters out false positives, i.e. components that are well
predictable based on poorly predictable components. These become expensive
as the error would propagate.
Taking k steps into account for prediction yields the following optimization problem:
k
X

minimize
A∈O(n)

kATr (z − WVi ζ(t − i)k2

(23)

t

i=0

It is globally solvable similarly to (19):
minimize

Tr

A∈O(n)



ATr

k
X


T
z − WVi ζ(t − i) z − WVi ζ(t − i)

A
t r



(24)

i=0

Solve (24) by diagonalizing

Pk

i=0


T
z − WVi ζ(t − i) z − WVi ζ(t − i)
(k)
Ar

t

and sorting the low-

denote the global solution of (23). How to optimally
est r eigenvalues to the upper left. Let
choose k for a certain problem is currently an open question, but we know from experiments that
in most cases increasing k (up to some unspecified value) improves the prediction error. Increasing
k too far, however, can lower the quality again, so it is up to experiments how to choose k. This
effect is illustrated in section 4.1, figure 8.

2.4

Taking supplementary information into account for prediction

The main idea in this work is to enable PFA to take arbitrary supplementary information into
account, yielding PFAx. This information shall not be used to extract data from it, but can
serve as a helper to improve prediction. Even more important is that PFAx can be used to find
subsignals for which the given supplementary information is most relevant. In the following, we
let u(t) denote this supplementary information.
To perform the desired analysis, we first need to extend the fitting formula (12) to take supplementary information into account. Our proposed prediction scheme is as follows:
m(t)

!

≈

B1 m(t − 1) + . . . + Bp m(t − p)

+

U1 u(t − 1) + . . . + Uq u(t − q)

=

B vec(histz,p (t)) + U vec(histu,q (t))
{z
}
|
{z
}
|
= ζ(t)

=: µ(t)

8

with Bi ∈ Rr×r
with Ui ∈ R

r×nu

(25)
(26)
(27)

z(t)

m(t)

PFAx

main input

u(t)

supplementary

Figure 5: Illustration of prediction using supplementary information. Components are selected to be well predictable if supplementary information is taken
into account.
For a given extraction matrix Ar the following formulas provide optimal values for B ∈ Rr×rp
and U ∈ Rr×nu q :

−1


(28)
zζ T − U(Ar ) µζ T Ar Ar T ζζ T Ar
B(Ar ) :=


−1
zµT − B(Ar ) Ar T ζµT
U(Ar ) :=
µµT
(29)
For convenience, since B(Ar ) and U(Ar ) are only defined implicitly, we provide an explicit formula
to obtain B(Ar ):




 −1
−1
−1
B(Ar ) :=
zζ T − zµT µµT
µζ T Ar Ar T ζζ T − ζµT µµT
µζ T Ar
(30)
If the arising inversions are not computable due to near-zero-eigenvalues, proceed like explained
in our comment on inverting ζζ T . In this section we define W := B(I) in analogy to (17),
but based on (30) rather than (12). Now we can adopt (20) to find an optimal Ar in terms of the
relaxation:


T
minimize
Tr ATr (z − Wζ − U(I)µ) (z − Wζ − U(I)µ) Ar
(31)
A∈O(n)

Just like (20), (34) can be solved by diagonalizing the central matrix and sorting the smallest
eigenvalues to the upper left. Applying the technique from section 2.3 to the optimization problem
of this section is a bit involved. We define


−1
:=
V
ζ(t + 1)ζ T − Inp,n U(I) µζ T
ζζ T
(32)
ẑ(i) (t)

:=

WVi ζ(t − i) + ITnp,n

i
X

Vj Inp,n U(I)µ(t − j)

(33)

j=0

Note that W and V simplify to their earlier definitions for U = 0. Also the predictor ẑ(i) simplifies
to the predictor defined in (22) for U = 0. For i = 0 it even boils down to Wζ, which was used as
predictor in (19). Consequently, we can also write the optimization problem in terms of ẑ(i) :
!
k
k
X
X


T
minimize
kATr (z − ẑ(i) )k2
= Tr ATr
z − ẑ(i) z − ẑ(i)
Ar
(34)
A∈O(n)

i=0

i=0

9

It can be solved by the familiar procedure:
Choose A such that it diagonalizes


Pk
(i)
(i) T
z − ẑ
z − ẑ
and sort the lowest r eigenvalues to the upper left.
i=0

3

Using extracted features to perform control

The idea here is that z represents some kind of an agent’s perception, while the supplementary
information u represents the agent’s control commands. The extracted predictable features are a
compact representation of perception aspects that are influenced by the control commands. Now
we consider the setting that after the training phase Ar , B(Ar ) and U(Ar ) are given and we want
to reach a goal position m∗ in feature space. We want to calculate what control must be chosen to
reach a desired position in feature space. For this purpose, we minimize the least square distance
between predicted features and goal features:
km∗ − B(Ar )ATr ζ(t + 1) − U(Ar )µ(t + 1)k2

minimize
n
u(t) ∈ R

u

=

km∗ − B(Ar )ATr ζ(t + 1) −

q
X

(35)


Uj (Ar )u(t − j + 1) −U1 u(t)k2

(36)

j=2

{z

|
=

=: u∗ ∈ Rr

}

ku∗ − U1 u(t)k2
|{z}

(37)

r×nu

∗
T
−1 T ∗
Solving this problem is straight forward by choosing u(t) := U−1
U1 u ,
1 u (or u(t) := (U1 U1 )
if U1 is not quadratic or not invertible). Note that for this solution squaring would not have been
necessary, i.e. it would also minimize ku∗ − U1 u(t)k. However, as soon as there are constraints on
u, the squared distance is much friendlier for optimization.

z(t)
control

m(t)

PFAx

perception

u(t)

Figure 6: Illustration of controlling predictable features. The influence of the
control signal onto prediction is inverted to calculate a control command that
would likely yield a desired outcome.
We will usually have constraints on u(t), reflecting the agent’s limited performance. A typical
case for future applications are linear constraints, which would yield an efficiently solvable quadratic
optimization problem (QP). In the following sections we will model a robot driving with constant
speed, so we need to deal with a normalized-length-constraint:
minimize
nu
u(t) ∈ R
ku(t)k = c

ku∗ − U1 u(t)k2

10

(38)

This problem is equivalent to the inhomogeneous eigenvalue problem
UT1 U1 u(t)
ku(t)k

= λu(t) + UT1 u∗

(39)

=

(40)

c

In [Mattheij and Söderlind, 1987] such problems are approached. For convenience, we provide a
method from there in A.3.

4

Experiments and Applications

In this section we apply the principles developed so far to some simulated problem settings. We
show how this is suitable to generate motor control for letting an agent reach a certain goal position
in a simple environment, based on different sensors. Once the environment gets complexer it turns
out that the principle is mainly suitable for local navigation and we show experiments that illustrate
this limitation. In a follow-up paper we will show how this principle can be extended to perform
global navigation reliably.

4.1

Using place cells for navigation

In this section we apply the developed techniques to the following setting: We have a virtual rat
– an agent – on a table and assume there is already a mechanism that obtains place cells from
vision. This is a reasonable assumption, because in [Franzius et al., 2007a, Franzius et al., 2007b,
Schönfeld and Wiskott, 2015] SFA has been shown to be capable of such a preprocessing. To find
out whether it is possible to use PFAx for goal-driven navigation, we place a certain number of 2DGaussians on the table to model place cells. That means, each component of x(t) is the activation
of one place cell, i.e. a 2D-Gaussian centered on a random position, evaluated for the rat’s position
at time t. In a training phase consisting of a random walk the rat explores the environment and
PFAx extracts the most predictable signals from place cell input x (z after sphering). Within
this section we always use 50 place cells, i.e. dim(x) = 50. To constrain the complexity of the
experiment we let the rat walk with constant speed. So the only control information is direction.
This is provided as a direction vector and fed into PFAx as supplementary information u (i.e. the
movement deltas are provided to PFAx).
After the training phase we choose a goal position and use the extraction matrix obtained by
PFAx to calculate corresponding goal features. Then we iterate the optimization from section 3,
where c is the rat’s constant speed. We actually apply the optimized u(t)-command after each
step, and increase t. Finally we observe the resulting path and modify the setup in various ways
to measure the robustness of this method.
Based on a random walk, PFAx can only find predictability w.r.t. to the provided movement
delta information. Since this incorporates just one time step, p = 1 is an exhaustively sufficient
value for all experiments. A higher value for p would only make sense if the movement simulation
of the agent would incorporate higher moments like acceleration. Adding such higher moments to
the simulation is postponed to future work.
Our first set of experiments is shown in figure 7. The yellow crosshairs always indicate the goal
and the colormap displays Euclidean distance of each point to the goal in r-dimensional feature
space. The navigation path is rendered in white on top of the colormap. To indicate the starting
point, some final steps of the random walk training phase are shown and the algorithmic navigation
11

is clearly distinguishable from the random walk by yielding a straight or slightly curved path. We
conclude from figure 7 that this kind of navigation is robust for a value of r that properly matches
the problem complexity, in this setting ideally r = 2. Larger values of r add irrelevant information
to the feature space and thus distract the navigation algorithm. However, for moderate values of r
this can be compensated by increasing the value of k – like shown in figure 8 or by increasing the
training phase – like shown in figure 9. Detecting a good value for r is subject of ongoing research.

1000 steps

10000 steps

1.0

1.0

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

r=2

0.0

0.0
0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

1.0

0.8

0.6

r=4
0.4

0.2

0.0

Figure 7: Compare r ∈ {2, 4}, |Ωt | ∈ {1000, 10000}, 50 place cells, k = 0
k=0

k=10

k=20

k=60

1.0

1.0

1.0

1.0

0.8

0.8

0.8

0.8

0.6

0.6

0.6

0.6

0.4

0.4

0.4

0.4

0.2

0.2

0.2

0.2

r=5

0.0

0.0

0.0

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

0.0
0.0

0.2

0.4

0.6

0.8

1.0

Figure 8: Varying k, |Ωt | = 8000, r = 5, 50 place cells
12

0.0

0.2

0.4

0.6

0.8

1.0

4000 steps

10000 steps

1.0

1.0

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

20000 steps

r=6

0.0

0.0

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

1.0

0.8

0.6

r=8
0.4

0.2

0.0

1.0

0.8

0.6

r=14
0.4

0.2

0.0

1.0

1.0

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

r=25

0.0

0.0
0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

Figure 9: More values of r and |Ωt |, 50 place cells, k = 0
13

0.4

0.6

0.8

1.0

4.2

Place-cell-based navigation with obstacle

To examine the capabilities of global navigation we add an obstacle to the environment. It turns
out that there exists no parameter set that would enable the Gaussian-modeled place-cell-based
navigation routine to measure the obstacle. Examining the color inside the obstacle we see that it
is not represented in feature space at all. Figure 10 illustrates an exemplary run that demonstrates
how the algorithm would want to go through the obstacle if collision detection is turned off. Note
that during training phase collision with the obstacle is prohibited.
1.0

0.8

0.6

0.4

0.2

0.0
0.0

0.2

0.4

0.6

0.8

1.0

Figure 10: r = 3, |Ωt | = 8000, 50 place cells, k = 0, collision detection turned off
We conclude that naively using Gaussians to model place cells is inherently unsuitable for
measuring obstacles, because in this model the place cells just shine through the walls of the
obstacle. So we refine the simulation to use a sensor that is sensitive to obstacle walls and cannot
look through them. Constructing place cells that do not shine through walls would not be feasible
with plain Gaussians. Raw vision input would lead to a high dimensional input, but we prefer
to keep the setting low dimensional for now and focus on the navigation task. We therefore
propose a sensor that is invariant under head direction, given that SFA has been shown to be
capable of finding head direction invariant features from vision input, see [Franzius et al., 2007a,
Franzius et al., 2007b, Schönfeld and Wiskott, 2015]. So, with SFA as a potential preprocessing
step in mind, these are reasonable simplifications.

4.3

Using a wall sensor for navigation

We introduce a virtual sensor that measures – for its current location – the visible fraction of each
wall segment in an overall 360° field of view. Figure 11 shows how the full field of view is split up
into sections occupied by each wall segment. In the example the sensor consists of 8 components,
one for each wall segment existing in the environment. The components belonging to the left and
right walls of the obstacle are invisible from the location illustrated in figure 11 and would emit a
value of 0 in that case.
Using a 360° field of view the sensor is by construction invariant under head direction. As an
initial experiment we show that navigation without obstacle is feasible. With 4 surrounding wall
14

2.0

1.5

1.0

0.5

0.0
0.0

0.5

1.5

1.0

2.0

Figure 11: Sensor that measures the fraction each wall segment occupies in an 360° field of view
segments the environment in figure 12 yields dim(x) = 4. Contrary to the place cells scenario we
apply a quadratic expansion on the sensor to compensate its low dimensionality. Expansions with
higher than quadratic degree provided no significant improvement for this setting. So we focus on
quadratic expansion for all wall sensor experiments in this paper. Further, we use a fixed-length,
exhaustive training phase of |Ωt | = 40000 steps for all such experiments.
1.0

0.8

0.6

0.4

0.2

0.0
0.0

0.2

0.4

0.6

0.8

1.0

Figure 12: Using the wall sensor for navigation without obstacle
Given that the navigation without obstacle works reliably with all sorts of parameter sets, we
do not provide experiments beyond the exemplary illustration in figure 12 and instead focus on
the setting with obstacle.

15

4.4

Wall-sensor-based navigation with obstacle

Using the wall sensor it turns out that an obstacle is now clearly represented in feature space – note
the different color inside the obstacle. Depending on r and k the feature space surrounding the
obstacle can guide the agent around. To show that the feature space itself is repelling the agent
from the obstacle, navigation is performed without collision detection. During training phase,
however, collision with any walls is prohibited.
In figure 13 we show a setting with a starting point rather far away from the obstacle. It
suggests that with the right parameters it is actually possible to avoid the obstacle.
With a starting point near the obstacle, like presented in figure 15, we observe somewhat equal
results – with the right values for r and k a proper navigation can be achieved, but the choice is
not entirely stable. While a higher value of k mostly improves the navigation, choosing it too high
can break it.

16

r=3

r=2
1.0

1.0

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

r=4

k=0

0.0

0.0
0.0

0.2

0.4

0.6

0.8

1.0

1.0

1.0

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

k=60

0.0

0.0
0.0

0.2

0.4

0.6

0.8

1.0
1.0

0.8

0.6

k=80
0.4

0.2

0.0

1.0

0.8

0.6

k=100
0.4

0.2

0.0

1.0

1.0

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

k=120

0.0

0.0
0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

Figure 13: Using the wall sensor for navigation with far obstacle
17

0.8

1.0

r=3

r=2
1.0

1.0

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

r=4

k=0

0.0

0.0
0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

1.0

0.8

0.6

k=60
0.4

0.2

0.0

1.0

1.0

0.8

0.8

0.6

0.6

0.4

0.4

k=80

0.2

0.2

0.0

0.0
0.0

0.2

0.4

0.6

0.8

1.0

1.0

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

1.0

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

k=100

0.0

0.0
0.0

0.2

0.4

0.6

0.8

1.0

1.0

1.0

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

k=120

0.0

0.0
0.0

0.2

0.4

0.6

0.8

1.0

1.0

1.0

1.0

0.8

0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

k=160

0.2

0.2

0.2

0.0

0.0

0.0

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

Figure 14: Using the wall sensor for navigation with near obstacle

18

Finally with a goal position closely behind the obstacle we present an example where the agent is
repelled from the obstacle, but can hardly find its way around it (figure 15). The sufficient value
k = 200 is rather unstable and already with k = 220 the agent would pass through the obstacle.
Running with collision detection, the wall itself would still guide the agent around the obstacle,
but even then the navigation is not very stable.

k=180

k=200

1.0

1.0

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

r=3

0.0

0.0
0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

Figure 15: Using wall sensor for navigation with a goal position closely behind the obstacle

5

Conclusion

We conclude that the presented approach can be used to learn the relation between a control
command signal and an agent’s state in an environment. Especially it can find a low dimensional
feature space that can efficiently represent the agent’s state in a way that allows goal driven
manipulation. We can generate the proper control signal to optimize the agent’s state locally and
approximately perform a gradient descent within feature space. This is, however, hardly sufficient
to perform global navigation, e.g. navigate around an obstacle. While some promising results
suggest that even global optimization is feasible, they highly depend on the choice of parameters r
and k. A theory for handling these parameters more systematically is subject to current research.
We are preparing a follow-up paper with an extension to this algorithm that enables stable global
navigation.

Acknowledgments
This work has been funded by a grant from the German Research Foundation (Deutsche
Forschungsgemeinschaft, DFG) to L. Wiskott (SFB 874, TP B3) and supported by the German
Federal Ministry of Education and Research within the National Network Computational Neuroscience – Bernstein Fokus: “Learning behavioral models: From human experiment to technical
assistance”, grant FKZ 01GQ0951.

19

References
[Box and Tiao, 1977] Box, G. E. P. and Tiao, G. C. (1977). A canonical analysis of multiple time
series. Biometrika, 64(2):pp. 355–365.
[Franzius et al., 2007a] Franzius, M., Sprekeler, H., and Wiskott, L. (2007a). Slowness and sparseness lead to place-, head direction-, and spatial-view cells. In Proc. 3rd Annual Computational
Cognitive Neuroscience Conference, Nov. 1–2, San Diego, USA, pages III–8.
[Franzius et al., 2007b] Franzius, M., Sprekeler, H., and Wiskott, L. (2007b). Slowness and sparseness lead to place, head-direction, and spatial-view cells. PLoS Computational Biology, 3(8):e166.
[Goerg, 2013] Goerg, G. (2013). Forecastable component analysis. In Dasgupta, S. and Mcallester,
D., editors, Proceedings of the 30th International Conference on Machine Learning (ICML-13),
volume 28, pages 64–72. JMLR Workshop and Conference Proceedings.
[Goldberger et al., 2004] Goldberger, J., Roweis, S. T., Hinton, G. E., and Salakhutdinov, R.
(2004). Neighbourhood components analysis. In Advances in Neural Information Processing
Systems 17 [Neural Information Processing Systems, NIPS 2004, December 13-18, 2004, Vancouver, British Columbia, Canada], pages 513–520.
[Mattheij and Söderlind, 1987] Mattheij, R. and Söderlind, G. (1987). On inhomogeneous eigenvalue problems. i. Linear Algebra and its Applications, 88-89(Supplement C):507 – 531.
[Richthofer and Wiskott, 2015] Richthofer, S. and Wiskott, L. (2015). Predictable feature analysis.
In 14th IEEE International Conference on Machine Learning and Applications, ICMLA 2015,
Miami, FL, USA, December 9-11, 2015, pages 190–196.
[Schönfeld and Wiskott, 2015] Schönfeld, F. and Wiskott, L. (2015). Modeling place field activity
with hierarchical slow feature analysis. Front Comput Neurosci, 9:51. 26052279[pmid].
[Sprague, 2009] Sprague, N. (2009). Predictive projections. In IJCAI 2009, Proceedings of the
21st International Joint Conference on Artificial Intelligence, Pasadena, California, USA, July
11-17, 2009, pages 1223–1229.
[Sprague, 2014] Sprague, N. (2014). Contingent features for reinforcement learning. In Artificial
Neural Networks and Machine Learning - ICANN 2014 - 24th International Conference on
Artificial Neural Networks, Hamburg, Germany, September 15-19, 2014. Proceedings, pages 347–
354.
[Weghenkel et al., 2017] Weghenkel, B., Fischer, A., and Wiskott, L. (2017). Graph-based predictable feature analysis. Mach. Learn., 106(9-10):1359–1380.
[Wiskott et al., 2011] Wiskott, L., Berkes, P., Franzius, M., Sprekeler, H., and Wilbert, N. (2011).
Slow feature analysis. Scholarpedia, 6(4):5282.

20

A
A.1

Appendix
Notation overview

This section gives an overview of the notation used in this paper.
x(t)
denotes the raw input signal.
u(t)
denotes the supplementary information signal.
:= {t0 , . . . , tk } denotes a discrete time sequence (considered as
Ωt
equidistant with step size normalized to 1). We usually refer to
Ωt as the
Ptraining phase.
1
hs(t)it∈S := |S|
t∈S s(t) denotes the average of some signal s over a
finite set S. For S = Ωt we just write hs(t)it or even hsi, if it
is obvious, what unbound variable is targeted.
h(x)
denotes the expansion function and usually consists of a set of
monomials of low degree.
z(t)
denotes h(x(t)) after sphering.
m(t)
denotes the optimized output signal (m for model ).
:= dim(h(x)) denotes the number of components to be analyzed
n
(after expansion).
:= dim(u) denotes the number of components in u
nu
r
denotes the number of extracted components (“features”).
A, a
denotes the matrix (or vector if r = 1) holding the linear composition of the output-signal. We set m(t) = AT z(t).
ai
denotes the i’th column of A, so we can write mi (t) = aTi z(t).
O(n)
⊂ Rn×n denotes the orthogonal group of dimension n, i.e. ∀ A ∈
O(n) : AAT = AT A = I
p
denotes the number of recent signal-values involved in the prediction. We also call it the prediction-order.
Is,r
denotes the s × r identity matrix (s counting rows, r counting
columns). For s = r this is a usual square identity, while in the
non-square case it consists of a square identity block in the top
or left area, filled up with zeroes to fit the given shape.
:= In,r
Ir
:= AIr
Ar
We frequently use the p-step time-history of a signal z, which we formalize by the following
function:
histz,p,∆ (t)

:=

p
X

z(t − i∆)eTi

with ei ∈ Rp

(41)

i=1

histz,p (t)

:=

histz,p,1 (t)

(42)

Here ei denotes the i-th p-dimensional euclidean unit vector, which is 1 at position i and 0 everywhere else.
Further more we sometimes use the Kronecker product ⊗ and the vec-operator defined as
follows:
21

For matrices A ∈ Rm×n and B ∈ Rk×l and with aij denoting the entries, ai the columns of A:


a11 B · · · a1n B

..  ∈ Rmk×nl
..
A ⊗ B :=  ...
(43)
.
. 
am1 B · · · amn B
 
a1
 .. 
(44)
vec(A) :=  .  ∈ Rmn
an
Additionally, we sometimes make use of the following shortcut:


A
0


..
A := Ip,p ⊗ A = 

.
0
|

{z

(45)

A
}

p times A

A.2

Extracting predictable single components

In section 2.2 we initially stated a prediction model that always scopes on single components. This
idea was not suitable for PFA because it contradicts the orthogonal agnosticity criterion. In this
section we propose a strategy to extract well predictable single components even though. We begin
by recalling our initial notion of linear auto regressive predictability:
aT z(t)

!

b1 aT z(t − 1) + . . . + bp aT z(t − p)

≈

T

=

a histz,p (t) b

(46)
(47)

It is possible to write this for multiple dimensions by constraining the coefficient-matrices to be
diagonal:
m(t)

!

≈

B1 m(t − 1) + . . . + Bp m(t − p)

with Bi ∈ Rn×n , diagonal

(48)

This model is not orthogonal agnostic, so a different approach than in section 2.2 is needed. To
minimize the least-squares-error of (46), the following optimization problem needs to be solved:
minimize
n
p

D

subject to

aT hzi

=

0

(zero mean)

aT zzT a

=

1

(unit variance)

a∈R , b∈R

2 E
aT (z − histz,p b)
(49)

Via analytic optimization it is straight forward to find the optimal a, if b is fixed and vice versa:
If b is fixed, choose a as the eigenvector corresponding to the smallest eigenvalue in
zzT − zbT histz,p − histTz,p bzT + histz,p bbT histTz,p

22

(50)

If a is fixed, choose b as
bT :=

histTz,p aaT histz,p

zT aaT histz,p

−1

By inserting (51) into (49) one could obtain a problem written in a only:

minimize
aT (z − histz,p zT aaT histz,p histTz,p aaT histz,p
n

(51)

−1

a∈R

subject to

2 
)
(52)

aT hzi

=

0

(zero mean)

aT zzT a

=

1

(unit variance)

Problem (52) is not efficiently globally solvable by any method known to us, which is mainly due
to the occurrence of a in a matrix-term under an inversion-symbol. However a possible strategy is
to approximate the solution by choosing an initial value for a or b and applying (50) and (51) in
turns until a stable state is reached.
As a reasonable initial value for this procedure we choose b such that it is the best predictor
of z on average, in absence of any a:
z(t)

!

≈

b1 z(t − 1) + . . . + bp z(t − p)

=

histz,p (t)b

(53)

To minimize the error of (53) on average over all components of z, we propose the following
least-squares optimization:
minimize
p
b∈R

(z − histz,p b)T (z − histz,p b)

(54)

The solution of this problem is
b :=

histTz,p histz,p

zT histz,p

−1

(55)

Solution (55) does not change, if we replace z by AT z with any orthogonal, full ranked A. However,
one quickly finds examples, where the procedure stabilizes in sub-optimal states. Though one can
partly overcome this issue by estimating better starting points, the method still has unknown
success-probability.
Probably a better possibility is to solve (49) with PFA as described in section 2 for r = 1.
After extracting one component either way, one can project z to the signal space uncorrelated
(i.e. orthogonal) to the extracted component. The extraction- and projection-procedure can be
repeated until any desired number of components is extracted.

A.3

Solving inhomogeneous eigenvalue problems

An inhomogeneous eigenvalue problem like (39) in general form can be stated as
Av

= λv + b

(56)

kvk

= c

(57)

with matrix A, vectors v and b, scalar c of appropriate dimensions.
While
[Mattheij and Söderlind, 1987] mainly focuses on direct numerical approaches, they also point out
23

a method to write (56) as an ordinary eigenvalue problem of higher dimension. Since inhomogeneous eigenvalues are different from homogeneous ones for b 6= 0, we can assume that (A − λI) is
invertible. Thus we can calculate v as
v

=

(A − λI)−1 b

(58)

and only have to find the inhomogeneous eigenvalues λ. To do so, we insert (58) into (57) and
obtain
bT (A − λI)−T (A − λI)−1 b = c
(59)
which implies that the block matrix

(A − λI)(A − λI)T
bT


b
c

(60)

can’t have full rank ((59) is the Schur complement of (60)). Consequently its determinant must
be zero. With the block matrix determinant formula


A B
det
= det(D) det(A − BD−1 C)
(61)
C D
we can state this as


bbT 
det (A − λI)(A − λI)T −
c

=

0

(62)

T

So ((A − λI)(A − λI)T − bbc ) must have at least one zero-valued eigenvalue and there exists a
corresponding eigenvector:
∃ x 6= 0 :


bbT 
λ2 x + λ(−AT − A)x + AAT −
x
c

=

0





(63)

(63) can be written as

∃ x 6= 0 :

T

bb
c

0
− AAT

I
T
A +A

 

x
λx


=

λ

x
λx

(64)

Finally we can 
use any ordinary eigenvalue algorithm to obtain λ as an eigenvalue of

0
I
.
bbT
−AAT AT +A
c

24

