Title: LEARNING GRID CELLS AS VECTOR REPRESENTA-TION OF SELF-POSITION COUPLED WITH MATRIX REPRESENTATION OF SELF-MOTION

Abstract: This paper proposes a representational model for grid cells. In this model, the 2D self-position of the agent is represented by a high-dimensional vector, and the 2D self-motion or displacement of the agent is represented by a matrix that transforms the vector. Each component of the vector is a unit or a cell. The model consists of the following three sub-models. (1) Vector-matrix multiplication. The movement from the current position to the next position is modeled by matrix-vector multiplication, i.e., the vector of the next position is obtained by multiplying the matrix of the motion to the vector of the current position. (2) Magnified local isometry. The angle between two nearby vectors equals the Euclidean distance between the two corresponding positions multiplied by a magnifying factor. (3) Global adjacency kernel. The inner product between two vectors measures the adjacency between the two corresponding positions, which is defined by a kernel function of the Euclidean distance between the two positions. Our representational model has explicit algebra and geometry. It can learn hexagon patterns of grid cells, and it is capable of error correction, path integral and path planning.

1. INTRODUCTION:

Imagine you are walking in your living room in the dark at night without any visual cues. Purely based on your self-motion, you know where you are while you are walking simply by integrating your self-motion. This is called path integral (Hafting et al. (2005) ; Fiete et al. (2008) ; McNaughton et al. (2006) ). You can also plan your path to the light switch or to the door. This is called path planning (Fiete et al. (2008) ; Erdem & Hasselmo (2012) ; Bush et al. (2015) ). You need to thank your grid cells for performing such navigation tasks. in 2005. Their thesis advisor, Dr. John O'keefe, shared the prize for his discovery of the place cells (O' Keefe (1979) ). Both the place cells and grid cells are used for navigation. The discoveries of these cells were made by recording the activities of the neurons of the rat when it moves within a square region. See Figure 1 (b) . Some neurons in the Hippocampus area are place cells. Each place cell fires when the rat moves to a particular location, and different place cells fire at different locations. The whole collection of place cells cover the whole square region. The discovery of grid cells was much more surprising and unexpected. The grid cells exist in the Entorhinal cortex. Each grid cell fires at multiple locations, and these locations form a regular hexagon grid. See Figure 1(c ). The grid cells have been discovered across mammalian species, including human. See Figure 1(d) .

In this paper, we propose a representational model to explain the hexagon patterns of the grid cells, and to explain how the grid cells perform path integral and path planning. We shall show that the grid cells are capable of error correction, which provides a justification for the grid cells.

v(x) v(x + ∆x) ×M (∆x) v(x) v(x + ∆x) ω|∆x| v(x) v(y) f (|x -y|)

(1) Vector-matrix multiplication

(2) Magnified local isometry (3) Global adjacency kernel Figure 2 : Grid cells form a high-dimensional vector representation of 2D self-position. Three submodels: (1) Local motion is modeled by vector-matrix multiplication.

(2) Angle between two nearby vectors magnifies the Euclidean distance.

(3) Inner product between any two vectors measures the adjacency which is a kernel function of the Euclidean distance.

Figure 2 illustrates our model. The 2D self-position of the agent is represented by a highdimensional vector, and the 2D self-motion or displacement of the agent is represented by a matrix that acts on the vector. Each component of the vector is a unit or a cell. In Figure 2 , x denotes the 2D self-position, v(x) is the high-dimensional vector representation of the 2D x. ∆x is the self-motion or one-step displacement. M (∆x) is the matrix representation of ∆x. The model consists of the following three sub-models. (1) Vector-matrix multiplication. The movement from the current position x to the next position x + ∆x is modeled by matrix-vector multiplication, i.e., the vector of the next position v(x + ∆x) is obtained by multiplying the matrix of the motion, M (∆x), to the vector of the current position v(x).

(2) Magnified local isometry. The angle between two nearby vectors equals the Euclidean distance |∆x| between the two corresponding positions multiplied by a magnifying factor. (3) Global adjacency kernel. The inner product between two vectors v(x), v(y) measures the adjacency between the two corresponding positions, which is defined by a kernel function f of the Euclidean distance |x -y| between the two positions x and y. One additional feature is that the whole vector v(x) is partitioned into multiple sub-vectors, and each sub-vector is driven by an associated sub-matrix. The whole system is like a multi-arm clock, with each arm rotating at a magnified speed and spanning a 2D sub-manifold on a high-dimensional sphere.

Our experiments show that sub-models (1) and ( 2) are sufficient for the emergence of the hexagon grid patterns of the grid cells. Sub-model (2) makes the vector representation robust to noises or errors due to the magnification of the distance. Sub-model (3) enables unique decoding of the position from the vector, because f is unimodal and peaked at 0, and it may serve as the link between the grid cells and place cells. Together with sub-model ( 1), sub-model (3) enables path integral and path planning, because the adjacency v(x), v(y) informs the Euclidean distance |x -y|. All the three sub-models can be implemented by one-layer neural networks.

2. CONTRIBUTIONS AND RELATED WORK:

The following are the contributions of our work.

(1) We propose a representational model for grid cells, where the self-position is represented by a vector and the self-motion is represented by a matrix that acts on the vector. (2) We show that our model can learn hexagon grid patterns.

(3) We show our model is capable of path integral, path planning, and error correction.

Many mathematical and computational models (Burak & Fiete (2009) ; Sreenivasan & Fiete (2011) ; Blair et al. (2007) ; de Almeida et al. (2009) ) have been proposed to explain the formation and function of grid cells. Compared to previous computational models on grid cells, our model only makes very generic assumptions about the algebra and geometry of the representational scheme, without assuming Fourier plane waves or clock arithmetics.

Recently, deep learning models (Banino et al. (2018) ; Cueva & Wei ( 2018)) have been proposed to learn the grid-like units for navigation. Our work was inspired by them. Compared with these models, our model has explicit algebra and geometry. In terms of algebra, our model has explicit matrix representation of self-motion, and the change of self-position is modeled by vector-matrix multiplication. In terms of geometry, our model assumes that the vector rotates while the agent moves, and our model assumes magnified local isometry and global adjacency kernel based on the angles between the vectors.

Expressing the adjacency kernel as the inner product between the vectors is related to the kernel trick (Cortes & Vapnik (1995) ), random Fourier basis (Ng et al. ( 2002)), and spectral clustering (Ng et al. (2002) ).

3. REPRESENTATIONAL MODEL OF GRID CELLS:

Consider an agent navigating within a domain D = [0, 1]×[0, 1] (actually the shape does not matter, and it can be R 2 ). We can discretize D into an N × N lattice. N = 40 in our experiments. Let x = (x 1 , x 2 ) ∈ D be the self-position of the agent. x is 2D. Suppose the agent wants to represent its self-position by a d-dimensional hidden vector v(x). We introduce the following three sub-models.

3.1. SUB-MODEL 1 ABOUT MOTION ALGEBRA: VECTOR-MATRIX MULTIPLICATION:

Suppose at a position x, the self-motion or one-step displacement is ∆x, so that the agent moves to

x + ∆x after one step. We assume that v(x + ∆x) = M (∆x)v(x),

(1) where M (∆x) is a d × d matrix that depends on ∆x. While v(x) is the vector representation of the self-position x, M (∆x) is the matrix representation of the self-motion ∆x. We can illustrate the motion model by the following diagram:

EQUATION

See Figure 2 (1). Both v(x) and M (∆x) are to be learned.

We can discretize ∆x, and learn a motion matrix M for each ∆x. We can also learn a parametric model for M . To this end, we can further parametrize

EQUATION

where the coefficients β are to be learned. The above may be considered a second-order Taylor expansion which is expected to be accurate for small ∆x ∈ ∆, where ∆ is the allowed collection of one-step displacements, e.g., ±3 grid points in each direction on the 40 × 40 grid.

The motion model can be considered a linear recurrent neural network (RNN). However, if we are to interpret M as the weight matrix, then the weight matrix is dynamic because it depends on the motion ∆x. One may implement it by discretizing ∆x, so that we have a finite set of ∆x, and thus a finite set of M (∆x). Then at each time step, the RNN switches between the finite set of motion matrices. This is like the gearing operation of a multi-speed bicycle.

3.2. DISENTANGLED BLOCKS OR MODULES:

For the sake of estimation accuracy and computational efficiency, we further assume that M (∆x) is block diagonal, i.e., we can divide v into K blocks of sub-vectors, v = (v (k) , k = 1, ..., K), and

v (k) (x + ∆x) = M (k) (∆x)v (k) (x)

. That is, the vector v consists of sub-vectors, each of which rotates in its own subspace, so that the dynamics of the sub-vectors are disentangled.

The assumption of disentangled blocks in our model is related to the modular organization of grid cells in neuroscience (Stensola et al. (2012) ), where a module refers to a region of grid cells where all cells share similar grid scales and orientations.

3.3 SUB-MODEL 2 ABOUT LOCAL GEOMETRY: MAGNIFIED LOCAL ISOMETRY

The above motion algebra alone is not sufficient for learning the vector matrix representations, because a trivial solution is that all the v(x) are the same, and the matrix M (∆x) is always identity. We need to properly displace the vectors v(x). So we shall model v(x), v(y) both locally and globally.

Let d be the dimensionality of v (k) , i.e., the number of grid cells within the k-th block. For the local geometry, we assume that for each block,

EQUATION

) for all x and ∆x such that α k |∆x| 2 ≤ c, i.e., ∆x ∈ ∆(α k ) = {∆x : α k |∆x| 2 ≤ c}. In our experiments, we take c = 1.5. α k can be either designed or learned.

Based on sub-model (4), v (k) (x) 2 = d for every x. The inner product on the left hand side is d cos(∆θ) where ∆θ is the angle between v (k) (x) and v (k) (x + ∆x). 1α k |∆x| 2 on the right hand side may be considered a second order Taylor expansion of a function f (r) such that f (0) = 1, f (0) = 0, i.e., 0 is the maximum, and

f (0) = -2α k . It is also an approximation to cos( √ 2α k |∆x|). Let ω k = √ 2α k , we have Magnified local isometry : ∆θ = ω k |∆x|,

(5) i.e., the angle between v(x) and v(x + ∆x) magnifies the distance |∆x| by a factor of ω k uniformly for all x. See Figure 2 (2).

The factor ω k defines the metric of block k.

3.4. ROTATION AND PROJECTION:

Since v (k) (x) is a constant for all x, M (k) (∆x) is an orthogonal matrix, and the self-motion is represented by a rotation in the d-dimensional space. ω k |∆x| is the angular speed of rotation of the sub-vector k. (v (k) (x), ∀x) forms a 2D sub-manifold on the sphere in the d-dimensional space. v (k) is like an arm of a clock except that the clock is not a 1D circle, but a 2D sub-manifold. This 2D sub-manifold becomes a local codebook for the 2D positions within a local neighborhood. For a vector v (k) , we can decode its position by projecting it onto the 2D sub-manifold, to get x = arg max x v (k) , v (k) (x) where the maximization is within the local neighborhood.

3.5. ERROR CORRECTION:

The neurons are intrinsically noisy and error prone. For ω k = √ 2α k 1, the magnification offers error correction because v (k) (x) and v (k) (x + ∆x) are far apart, which is resistant to noises or corruptions. That is, projection to the 2D sub-manifold codebook removes the noises.

Specifically, suppose we have a vector u = v (k) (x) + , where ∼ N(0, s 2 I d ), and I d is the d-dimensional identity matrix. We can decode x from u based on the codebook by maximizing u, v (k) 

(y) = v (k) (x), v (k) (y) + , v (k) (y) = d(1 -α k |y -x| 2 ) +

√ dsZ over y that is within a local neighborhood of x, where Z ∼ N(0, 1). The optimal y will be close to x because if y deviates from x by ∆x, then the first term will drop by dα k |∆x| 2 , which cannot be made up by the second term √ dsZ due to noises, unless the noise level s is extremely large. From the above analysis, we can also see that the larger d is, the more resistant the system is to noise, because the first term is of order d and the second term is of order √ d.

In addition to additive noises, the system is also resistant to multiplicative noises including dropout errors, i.e., multiplicative Bernoulli 0/1 errors. The dropout errors may occur due to noises, aging, or diseases. They may also be related to the asynchronous nature of the neuron activities in computing.

3.6. HEXAGON GRID PATTERNS:

While the magnified local isometry enables error correction, it also causes global periodicity. Because the angle between nearby v (k) (x) and v (k) (x + ∆x) is magnified, when we vary x, the vector v (k) (x) will rotate at a magnified speed ω k |∆x|, so that it quickly rotates back to itself, like an arm of a clock. Thus each unit of v (k) (x) is periodic with ω k determining the periodicity.

Our experiments show that sub-models ( 1) and ( 2) are sufficient for the emergence of hexagon grid patterns. In fact, we have the following analytical solution (see Appendix for a proof): Theorem 1. Let e(x) = (e i aj ,x , j = 1, 2, 3) , and a 1 , a 2 , a 3 are three 2D vectors so that the angle between a i and a j is 2π/3 for ∀i = j and |a j | = 2 √ α for ∀j. Let C be a random 3 × 3

complex matrix such that C * C = I. Then v(x) = Ce(x), M (∆x) =

Cdiag(e(∆x))C * satisfy equation 1 and equation 4 approximately for all x and small ∆x.

v(x) amounts to a 6-dimensional real vector. Since the angle between a i and a j is 2π/3 for ∀i = j, patterns of v(x) over x have hexagon periodicity. Moreover, the scale of the patterns is controlled by the length of a j , i.e., the scaling parameter α.

We want to emphasize that sub-models ( 1) and ( 2) are about local ∆x, where we do not make any assumptions about global patterns, such as Fourier basis. In contrast, the solution in the above theorem is global. That is, our model assumes much less than the solution in the theorem. Our experiments show that the hexagon patterns will emerge as long as the number of units is greater than or equal to 6.

3.7. SUB-MODEL 3 ABOUT GLOBAL GEOMETRY: ADJACENCY KERNEL:

Because of the periodicity, each block (v (k) (x), ∀x) does not form a global codebook of 2D positions, i.e., there can be k) does not encode x uniquely. We can combine multiple blocks to resolve the global ambiguity. Specifically, let v(x) = (v (k) (x), k = 1, ..., K) be the whole vector, we assume the following global adjacency sub-model for the whole vector:

EQUATION

EQUATION

where (v(x), M (∆x), α k , ∀x, ∆x, k) are to be learned.

Recall d is the number of grid cells in each block, and K is the number of blocks. f (r) is the adjacency kernel that decreases monotonically as the Euclidean distance r = |x -y| increases.

One example of f is the Gaussian kernel f (r) = exp -r 2 /2σ 2 . Another example is the exponential kernel f (r) = exp (-r/σ) . As a matter of normalization, we assume f (0) = 1, which is the maximum of f (r).

Since f (0) = 1, v(x) 2 = Kd for any x, and v(x), v(y) = (Kd) cos θ, where θ is the angle between v(x) and v(y), and we have Global adjacency : cos θ = f (|x -y|).

(7) The angle between any two vectors is always less than π/2. See Figure 2 (3).

By fitting the multiple sub-vectors together, we still retain the error correction capacity due to magnified local isometry, meanwhile we eliminate the ambiguity by letting v (k) (x), v (k) (y) for different k cancel each other out by destructive interference as y moves away from x, so that we obtain unique decoding of positions. Let C = {v(x), x ∈ D} be the codebook sub-manifold, error correction of a vector u is obtained by projection onto C: arg max v∈C u, v .

The whole vector v is like a K-arm clock, with each v (k) being an arm rotating at a speed ω k |∆x| = √ 2α k |∆x|.

3.8. LOCALIZATION AND HEAT MAP:

(v(x), ∀x) forms a global codebook for x. It is a 2D sub-manifold on the sphere in the (Kd)dimensional space. For a vector v, we can decode its position by its projection on the codebook manifold. Since f (r) is monotonically decreasing, h(x) = v, v(x) gives us the heat map to decode the position of v uniquely. Let the decoded position be x, then x = arg max x v, v(x) . We can obtain the one-hot representation δ x of x by non-maximum suppression on the heat map h(x).

Let V = (v(x), ∀x) be the (Kd) × N 2 matrix (recall the domain D = [0, 1] 2 is discretized into an N × N lattice)

, where each column is a v(x). We can write the heat map h(x) = v, v(x) as a N 2 -dimensional vector h = V v, which serves to decode the position x encoded by v. Conversely, for a one-hot representation of a position x, i.e., δ x , which is a one-hot N 2 -dimensional vector, we can encode it by v = V δ x . Both the encoder and decoder can be implemented by a linear neural network with connection weights V and V respectively, as illustrated by the following diagram:

EQUATION

Note that in decoding v → h(x) → δ x and encoding δ x → v, we do not represent or operate on the 2D coordinate x explicitly, i.e., x itself is never explicitly represented, although we may use the notation x in the description of the experiments.

3.9. PATH INTEGRAL:

Path integral (also referred to as dead-reckoning) is the task of inferring the self-position based on self-motion (e.g., imagine walking in a dark room). Specifically, the input to path integral is a previously determined initial position x 0 and motion sequences {∆x 1 , ..., ∆x T }, and the output is the prediction of one's current position x T . We first encode the initial position x 0 as v(x 0 ). Then, by the motion model, the hidden vector v(x T ) at time T can be predicted as:

EQUATION

We can then decode x T from v(x T ).

3.10. PATH PLANNING:

Our representation system can plan direct path from the starting position x 0 to the target position y by steepest ascent on the inner product v(x), v(y) , i.e., let v 0 = v(x 0 ), the algorithm iterates

EQUATION

EQUATION

)

where ∆ is the set of allowable displacements ∆x.

When a rat does path planning, even if it is not moving, its grid cells are expected to be active. This may be explained by the above algorithm. In path planning, the rat can also fantasize bigger step sizes that are beyond its physical capacity, by letting ∆ include physically impossible large steps.

In general, our representation scheme

(v(x), M (∆x), f (|x -y|)) mirrors (x, ∆x, |x -y|).

Thus the learned representation is capable of implementing existing path planning algorithms in robotics (Siegwart et al. (2011) ) even though our system does not have explicit coordinates in 2D (i.e., the 2D coordinates x are never explicitly represented by two neurons).

3.11. PLACE CELLS AND SCALE:

For each x, we may interpret v(x), v(y) /(Kd) = f (|x -y|) as a place cell whose response is f (|x -y|) if the agent is at y, or at least, we may consider f (|x -y|) as an internal input to the place cell based on self-motion, in addition to external visual cues.

For the Gaussian kernel f (r) = exp(-r 2 /2σ 2 ), the choice of σ determines the scale. In path integral, for accurate decoding of the position x from the vector v in the presence of noises, we want σ to be small so that f (r) drops to zero quickly. However, in path planning, for a vector v, we also want to know the Euclidean distance between its position x to a target position y, which may be far away from x. The distance is |x -y| = f -1 ( v, v(y) /(Kd)). For accurate estimation of long distance in the presence of noises, we need σ to be large, so that the slope of f -1 is not too big. Perhaps we need multiple f (r) with different σ, and for each f (r) we have

(Kd)f (r) = K k=1 γ k v (k) (x), v (k) (y)

, where different f (r) have different coefficients (γ k ) while sharing the same (v (k) (x)). We shall study this issue in future work.

3.12. GROUP REPRESENTATION IN MOTOR CORTEX:

Where does the self-motion ∆x come from? It comes from the movements of head and legs, i.e., each step of navigation involves a whole process of path integral and path planning of the movements of head and legs (as well as arms). In general, we can use the same system we have developed for the movements of head, legs, arms, fingers, etc. Their movements form groups of actions. In navigation, the movements belong to 2D Euclidean group. In body movements, the movements belong to various Lie groups. Our method can be used to learn the representational systems of these groups, and such representations may exist in motor cortex. We leave this problem to future investigation.

4. LEARNING REPRESENTATION:

The square domain D is discretized into a 40 × 40 lattice, and the agent is only allowed to move on the lattice. We can learn (v(x), ∀x) and (M (∆x), ∀∆x) (or the β coefficients that parametrize M ) by minimizing the following loss functions.

For sub-model (1) on vector-matrix multiplication, the loss function is

EQUATION

)

where x is sampled from the uniform distribution on D = [0, 1] 2 , and ∆x is sampled from the uniform distribution within a certain range ∆ 1 (3 grid points in each direction in our experiments).

The above motion loss is a single-step loss. It can be generalized to multi-step loss

EQUATION

where (∆x t , t = 1, ..., T ) is a sequence of T steps of displacements, i.e., a simulated trajectory.

For sub-model (2) on magnified local isometry, the loss function is

EQUATION

)

where for fixed α k , ∆x is sampled uniformly within the range ∆ 2 (α k ) = {∆x : α k |∆x| 2 ≤ c} (c = 1.5 in our experiments). We can define L 2 = K k=1 L 2,k . For sub-model (3) on global adjacency kernel, the loss function is

EQUATION

)

where both x and y are sampled uniformly from D.

Let v(x) = (v i (x), i = 1, ..., Kd), we also impose a regularization loss

L 0 = Kd i=1 (E x [v i (x) 2

] -1) 2 , to enforce uniform energy among the grid cell. It also helps to break the symmetry caused by the fact that the loss function is invariant under the transformation v(x) → Qv(x), ∀x, where Q is an arbitrary orthogonal matrix. This loss is not crucial though, so we will make it implicit for the rest of the paper.

Fixing the magnitude of E x [v i (x)] within a certain range is biologically plausible, because v i (x) is a single cell. For mathematical and computational convenience, we can also normalize v(x) so that v(x

) 2 = 1, v(x), v(y) = f (|x -y|), v (k) (x), v (k) (x + ∆x) = (1 -α k |∆x| 2 )/K, and E x [v i (x) 2 ] = 1/(Kd). When learning a single block, we can normalize v (k) (x) so that v (k) (x) 2 = 1, v (k) (x), v (k) (x + ∆x) = 1 -α k |∆x| 2 and E x [v i (x) 2 ] = 1/d.

The total loss function is a linear combination of the above losses, where the weights for combining the losses are chosen so that the weighted losses are of the same order of magnitude.

The loss function is minimized by Adam optimizer (Kingma & Ba ( 2014)) (lr = 0.03) for 6, 000 iterations. A batch of 30, 000 examples, i.e., (x, ∆x t , t = 1, ..., T )) for L 1 , (x, ∆x) for L 2 , and (x, y) for L 3 are sampled at each learning iteration as the input to the loss function. In the later stage of learning (≥ 4, 000 iterations), v(x) = 1 is enforced by projected gradient descent, i.e., normalizing each v(x) after the gradient descent step. We first learn a single block with fixed α k by minimizing L 1 + λ 2 L 2,k (we shall drop the subscript k in this subsection for simplicity). Figure 3 shows the learned units over the 40 × 40 lattice of x. We learn multiple blocks by minimizing L 1 + λ 2 L 2 + λ 3 L 3 . Instead of manually assigning α k , we learn α k by gradient descent, simultaneously with v and M .

5. EXPERIMENTS:

In Figure 4 (a), we show the learned units v(x) over the 40 × 40 lattice of x and the learned metrics α k . A Gaussian kernel with σ = 0.08 is used for the global adjacency measure f (|x -y|). Block size is set to 6 and each row shows the learned units belonging to the same block. The scales of the firing fields are controlled by the learned α k .

Figure 4 (b) illustrates the combination of multiple blocks. For the localization model, given a vector v, the heat map of a single block v (k) , v (k) (x) has periodic firing fields and cannot determine a location uniquely. However, ambiguity disappears by combing the heat maps of multiple blocks, which have firing fields of multiple scales and phases that add up to a Gaussian kernel v, v(x) . The Gaussian kernel informs the place cell, by which the location is determined uniquely. For a motion ∆x, every block rotates in its own subspace with motion matrix M (k) (∆x), resulting in phase shifting in the heat map of each single block.

In the subsequent experiments, we shall learn the grid cells by minimizing L 1 + λ 3 L 3 for simplicity. Figure 5 (a) shows an example of path integral result (time duration T = 40), where we use single step motion loss L 1,T =1 in learning. Gaussian kernel with σ = 0.08 is used as the adjacency measure in L 3 . We find single-step loss is sufficient for performing path integral. The mean square error remains small (∼ 1.2 grid) even after 400 steps of motions (figure 5(b)). The error is averaged over 1, 000 episodes. The motion loss can be generalized to multi-step L 1,T , as shown by equation 13. In Figure 5 (b), we show that multi-step loss can improve the performance slightly.

In Figure 5 (c) we compare the learned models with fixed number of units (96) but different block sizes. We also compare the performance of models using Gaussian kernel (σ = 0.08) and exponential kernel (σ = 0.3) as the adjacency measure in the localization model. The result shows that models with Gaussian kernel and block size ≥ 3, and with exponential kernel and block size ≥ 4 have performances comparable to the model learned without block-diagonal assumption (block size = 96).

For path planning, we assume continuous x and ∆x. First, we design the set of allowable motions ∆. For a small length r, we evenly divide [0, 2π] into n directions {θ i , i = 1, ..., n}, resulting in n candidate motions ∆ = {∆x i = (r cos(θ i ), r sin(θ i )), i = 1, ..., n}. These n small motions serve as motion basis. Larger motion can be further added to ∆ by estimating the motion matrix M (k∆x i ) = M k (∆x i ). The starting position and destination can also be any continuous values, where the encoding to the latent vector is approximated by bilinear interpolation of nearest neighbors on the lattice.

In the experiments, we choose r = 0.05 and n = 100, and add another set of motions with length 0.025 to enable accurate planning. The system is learned with exponential kernel (σ = 0.3) as global adjacency to encourage connection of long distance. Figure 6 (a) shows planning examples with six settings of motion ranges ∆. Including larger motions accelerates the planning process so that it finishes with less steps. We define one episode to be a success if the distance between the end position and the destination is less than 0.025. We achieve a success rate of > 99% over 1, 000 episodes for all the six settings. The learned system can also perform planning with obstacles, where the global adjacency between the agent's current position and an obstacle serves as a repulsive forces. Specifically, suppose z is an obstacle to avoid, the agent can choose the motion ∆x t at time t by

EQUATION

)

where a and b are the scaling and annealing parameters. Figure 6 (b) shows the planning result with a dot obstacle laid on the direct path between the starting position and destination, with tuning of a and b. We choose a = 0.5 and b = 6 in subsequent experiments.

Now suppose we have more complicated obstacles {z i } m i=1 . They can be included by summing over the kernels of every obstacle {a v(z i ), M (∆x)v t b } m i=1 and choosing ∆x t at time t by

∆x t = arg max ∆x∈∆ v(y), M (∆x)v t - m i=1 a v(z i ), M (∆x)v t b . Figure 6(c)

shows some examples, where the obstacles mimicking walls, large objects and simple mazes.

The above method is related to the potential field method in robotics (Siegwart et al. (2011) ).

6. DISCUSSION: ROTATIONIST-CONNECTIONIST MODEL?:

In terms of general modeling methodology, a typical recurrent network is of the form v t = tanh(W (v t-1 , δ t )), where v t is the latent vector, δ t is the input change or action. v t-1 and δ t are concatenated and linearly mixed by W , followed by coordinate-wise tanh non-linearity. We replace it by a model of the form v t = M (δ t )v t-1 , where M (δ t ) is a matrix that is non-linear in δ t . M (δ t ) is a matrix representation of δ t . v t can be interpreted as neuron activities. We can discretize the value of δ t into a finite set {δ}, and each M (δ) can be stored as synaptic connection weights that drive the neuron activities. In prediction, the input δ t activates M (δ t ). In planning or control, all the {M (δ)} are activated, and among them the optimal δ is chosen.

The matrix representation of the above model is inspired by the group representation theory, where the group elements are represented by matrices acting on the vectors (Fulton & Harris (2013) ). It underlies much of modern mathematics and holds the key to the quantum theory (Zee (2016) ). Perhaps it also underlies the visual and motor cortex, where neurons form rotating sub-vectors driven by matrices representing groups of transformations. One may call it a rotationist-connectionist model.

PROJECT PAGE:

http://www.stat.ucla.edu/ ˜ruiqigao/gridcell/main.html A PROOF FOR SECTION 3.6

Proof. (a j , j = 1, 2, 3) forms a tight frame in the 2D space, in that for any vector x in 2D,

EQUATION

)

Then we have

EQUATION

EQUATION

where

EQUATION

where D(∆x) = diag(e i aj ,∆x , j = 1, 2, 3).

If K > 1 and block size = 6, we can fit the multiple blocks together by a Fourier expansion of the kernel function

EQUATION

)

B HEXAGON GRID PATTERNS AND METRICS

B.1 SIMULATED INPUT DATA:

We obtain input data for learning the model by simulating agent trajectories with the number of steps equal to T in the square domain D. D is discretized into a 40 × 40 lattice and the agent is only allowed to move on the lattice. The agent starts at a random location x 0 . At each time step t, a small motion ∆x t (≤ 3 grids in each direction) is randomly sampled with the restriction of not leading the agent outside the boundary, resulting in a simulated trajectory {x 0 + t i=1 ∆x i } T t=1 . T is set to 1, 000 to obtain trajectories that are uniformly distributed over the whole area. Although the trajectories used for training are restricted to the lattice and with small motions, in Section 5.4 we show that the learned model can be easily generalized to handle continuous positions and large motions. In training, pairs of locations (x, y) are randomly sampled from each trajectory as the input to the adjacency loss L 3 , while consecutive position sequences (x, ∆x t , t = 1, ..., T ) are randomly sampled as the input to the motion loss L 1 , with length specified by T (which is usually much smaller than the whole length of the trajectory T ).

B.2 LEARNED SINGLE BLOCK UNITS WITH DIFFERENT BLOCK SIZES:

In (Blair et al. (2007) ), the grid cell response is modeled by three cosine gratings with different orientations and phases. In our model, we learn such patterns of different scales without inserting artificial assumptions. 

B.4 QUANTITATIVE ANALYSIS OF SPATIAL ACTIVITY:

We assess the spatial activity of the learned units quantitatively using measures adopted from the neuroscience literature. Specifically, we quantify the hexagonal regularity of the grid-like patterns using the gridness score (Langston et al. (2010) ; Sargolini et al. (2006) ). The measure is derived from the spatial autocorrelogram of each unit's response map. A unit is classified as a grid cell if its gridness score is larger than 0. For those units that are classified as grid cells, grid scale and orientation can be further derived from the autocorrelogram following Sargolini et al. (2006) . Figure 10 (a) summarizes the results. 76 out of 96 learned units are classified as grid cells. Most units with large learned α k are classified as grid cells, while those units with small α k are not due to the lack of a full period of hexagonal patterns. The grid scales and orientations vary among units, while 

B.5 ABLATION STUDIES:

We conduct ablation studies to assess various assumptions in our model.

B.5.1 LOSS TERMS:

We learn v and M with multiple blocks by minimizing L 1 + λ 2 L 2 + λ 3 L 3 , which consists of (1) a motion loss L 1 , (2) a local isometry loss L 2 , and (3) a global adjacency loss L 3 . Figure 11 shows the learned units when using only some of the components. Grid-like patterns do not emerge if using only the adjacency loss L 3 , or only the isometry loss L 2 . If using only the motion loss, the motion equation is approximately satisfied at the beginning of training, since both v and M are initialized from small values that are close to 0. The system cannot be learned without an extra term to push v(x) apart from each other. Figure 11 (c) shows the learned units using the adjacency loss and the motion loss, leaving out the isometry loss. Grid-like patterns still emerge (also some strip-like patterns), although less obvious than the ones learned using the full loss.

B.5.2 ASSUMPTIONS OF MOTION MATRIX:

In another ablation study, we drop the quadratic parametrization by β coefficients and the block diagonal assumption of the motion matrix. A separate motion matrix M (∆x) is learned for each displacement ∆x, and v(x) is assumed to be a single block. We use either the local isometry loss L 2 or the global adjacency loss L 3 in addition to the motion loss L 1 . The results are shown in Figure 12 . With the isometry loss, the setting is similar to the one of learning single block, except that the parametrization of M (∆x) is dropped. As shown in figure 12 (a), the learned units resemble the ones learned with parametrized M (∆x), but when the block size is large, the grid-like patterns are less obvious. With the adjacency loss, grid-like patterns do not emerge any more. We can model M (δ, v(θ)) by an attention (or selection) mechanism:

EQUATION

EQUATION

where M (i) (δ) is the matrix representation of scalar motion δ given the head direction θ i . The inner product v(θ), v(θ i ) , that informs the angular distance between θ and θ i , serves as the attention weight. b is an annealing (inverse temperature) parameter. If b → ∞, p = (p i , i = 1, ..., N ) becomes a one-hot vector for selection. We can further assume that M (∆θ) and M (i) (δ) are block diagonal, and learn a parametric model for each of them by the second order Taylor expansion at ∆θ and δ respectively.

C.1.2 LOCALIZATION SUB-MODEL:

For the localization sub-model, we define the adjacency measures of self-direction and self-position separately. Let f (|θ 1θ 2 |) be the adjacency measure between two directions θ 1 and θ 2 . We use von Mises kernel f (r) = exp((cos(r) -1)/σ 2 ), where f (0) = 1. For adjacency measure between self-positions x and y, we keep it the same as described in section 3.7. 

EQUATION

and then we learn v and M (or the coefficients that parametrize M ) as before.

C.2 LEARNED UNITS FOR SELF-DIRECTION AND SELF-POSITION:

Figure 13 shows a result of learning such an egocentric motion model by displaying the response curves of the learned units in the head direction system, v(θ), for θ ∈ [0, 2π], as well as the response maps of the learned multiple block units in the self-position system, v(x), for x ∈ [0, 1] 2 .

C.3 CLOCK AND TIMESTAMP:

We may re-purpose the head direction system as a clock, by interpreting θ ∈ [0, 2π] as the time on a clock, and v(θ) as a timestamp for events happening over time. This may be related to the recent neuroscience observations in Tsao et al. (2018) .

D ERRORS D.1 ERROR CORRECTION:

Unlike commonly used embedding in machine learning, here we embed a 2D position into a highdimensional space, and the embedding is a highly distributed representation or population code. The advantage of such a redundant code lies in its tolerance to errors. We show that the learned system is tolerant to various sources of errors. Specifically, in both path integral and path planning tasks, at every time step t, we randomly add (1) Gaussian noises or (2) dropout masks to the hidden units and see if the system can still perform the tasks well. We find that the decoding-encoding process (DE) is important for error correction. That is, at each time step t, given the noisy hidden vector v t , we decode it to x t = arg max x v t , v(x) and then re-encode it to the hidden vector v(x t ). Actually the whole process can be accomplished by projection to the codebook sub-manifold without explicit decoding, by obtaining the vector: arg max v∈C v t , v , where C = {v(x), x ∈ [0, 1] 2 }.

Table 1 shows the error correction results tested on path integral and path planning tasks. Each number is averaged over 1, 000 episodes. We compute the overall standard deviation of {v(x)} for all x and treat it as the reference standard deviation (s) for the Gaussian noise. For dropout noise, we set a percentage to drop at each time step. With the decoding-encoding process, the system is quite robust to the Gaussian noise and dropout error, and the system still works even if 70% units are silenced at each step.

D.2 NOISY SELF-MOTION INPUT:

Besides adding noises to the hidden units, we also experiment with adding noises to the self-motion ∆x, and compare the performance of path integral with the one performed in the original 2D coordinates. Specifically, at each time step, we add Gaussian noises to the self-motion ∆x. For path integral, we compute the mean square error between the predicted locations and ground truth locations. Besides, we also compute the predicted locations using the 2D coordinates with the noisy self-motions. Its mean square error serves as the reference error. Table 2 shows the result, indicating that the error of the learned system is close to the error in the original 2D coordinates, i.e., our system does not blow up the noise. 

E 3D ENVIRONMENT AND 1D TIME:

The system can be generalized to 3D environments. Specifically, we assume the agent navigates within a domain D = [0, 1] × [0, 1] × [0, 1], which is discretized into a 40 × 40 × 40 lattice. We learn a parametric model for motion matrix M in a residual form M = I + M (∆x), where M (∆x) is approximated by the second order Taylor expansion of ∆x.

We first learn a single block with fixed α k by minimizing L 1 +λ 2 L 2,k . A batch of 500,000 examples of (x, y) and (x, ∆x t , t = 1, ..., T ) is sampled online at every iteration for training. Figure 14 shows the learned units. Next we learn multiple blocks by minimizing L 1 + λ 3 L 3 , and use the learned models to perform 3D path integral and path planning. For simplicity, we remove L 2 . We use 8 blocks of units with block size 8 and exponential kernel (σ = 0.3) for path planning. A batch of 200,000 examples is sampled online at every iteration for training.

E.1 3D PATH INTEGRAL:

Figure 15 shows some results of 3D path integral with duration T = 30. Gaussian Kernel (σ = 0.08) is used as the adjacency measure.

E.2 3D PATH PLANNING:

Figure 16 shows some results of 3D simple path planning. Exponential kernel (σ = 0.3) is used as the adjacency measure. We design a set of allowable motions ∆: m lengths of radius r are used, and for each r, we evenly divide the inclination θ ∈ [0, π] and the arimuth α ∈ [0, 2π] into n directions, which results in a motion pool with mn 2 candidate motions ∆ = {(r sin(θ) cos(α), r sin(θ) sin(α), r cos(θ))}. We use m = 2, n = 90 in this experiment.

E.3 3D PATH PLANNING WITH OBSTACLES:

Figure 17 shows some examples of 3D path planning with a cuboid obstacle. a = 38 and b = 24 in equation 16. 

E.4 LEARNING IN 1D:

The system can also be applied to 1D. Inspired by Tsao et al. (2018) , the learned system in 1D may serve as a timestamp of events. We assume domain D = [0, 1] and discretize it into 100 time points. A parametric M is learned by a residual form M = I + M (∆t), where each element of M (∆t) is parametrized as a function of (∆t, ∆t 2 ). 16 blocks of hidden units with block size 6 are used.

Figure 18 visualizes the learned units over the 100 time points. The response wave of every unit shows strong periodicity of a specific scale or frequency. Within each block, the response waves of units have similar patterns, with different phases. 

:



Figures and Tables:

Figure 1: Place cells and grid cells. (a) The rat is moving within a square region. (b) The activity of a neuron is recorded. (c) When the rat moves around (the curve is the trajectory), each place cell fires at a particular location, but each grid cell fires at multiple locations that form a hexagon grid. (d) The place cells and grid cells exist in the brains of both rat and human. (Source of pictures: internet)

Figure 3: Learned units of a single block with fixed α. (a) Learned single block with 6 units. Every row shows the learned units with a given α. (b) Learned single block with 100 units and α = 72.

Figure 4: (a) Response maps of learned units of the vector representation and learned scaling parameters α k . Block size equals 6 and each row shows the units belonging to the same block. (b) Illustration of block-wise activities of the units (where the activities are rectified to be positive).

Figure 5: (a) Path integral prediction. The black line depicts the real path while red dotted line is the predicted path by the learned model. (b) Mean square error over time step. The error is average over 1, 000 episodes. The curves correspond to different numbers of steps used in the multi-step motion loss. (c) Mean square error performed by models with different block sizes and different kernel types. Error is measured by number of grids.

Figure 6: (a) Planning examples with different motion ranges. Red star represents the destination y and green dots represent the planned position {x0 + t i=1 ∆xi}. (b) Planning examples with a dot obstacle. Left figure shows the effect of changing scaling parameter a, while right figure shows the effect of changing annealing parameter b. (c) Planning examples with obstacles mimicking walls, large objects and simple mazes.

Figure 7: Response maps of learned single block units with different block sizes

each block. We show the histograms of the grid orientations and scales in Figure10(b) and 10(c) respectively. Moreover, we average the grid scales within each block and make the scatter plot of the averaged grid scales and the learned 1/ √ α k in figure10(d). Interestingly, the grid scale is nicely proportional to the learned 1/ √ α k .

Figure 8: Learned multiple block units and metrics with different block sizes

Figure 10: (a) Autocorrelograms of the learned units' response maps. Gridness scores are calculated based on the autocorrelograms. A unit is classified as a grid cell if the gridness score is larger than 0. The gridness score is shown in red color if a unit fails to be classified as a grid cell. For those units that are classified as grid cells, gridness score, scale and orientation are listed sequentially in black color. Orientation is computed using a camera-fixed reference line (0 • ) and in counterclockwise direction. (b) Histogram of grid orientations. (c) Histogram of grid scales. (d) Scatter plot of averaged grid scales within each block versus the corresponding learned 1/ √ α k .

(a) L3 only. (b) L2 only. (c) L1 + λ3L3, without L2.

Figure 11: Ablation study of the components in the training loss. (a) Learn the model using only the localization loss with global adjacency. (b) Learn the model using only the localization loss with local adjacency. (c) Using global adjacency and motion loss, leaving out local adjacency.

Figure 12: Learned units by dropping the parametrization and the block diagonal assumption of the motion matrix M (∆x).

Figure 13: (a) Response curves of learned units in the head direction model, v(θ). Each block shows the units belonging to the same sub-vector in the model. The horizontal axis represents the angle within a range of [0, 2π], while the vertical axis indicates the values of responses. (b) Response maps of the learned multiple block units of the self-position model, v(x), for x ∈ [0, 1] 2 .

Figure 14: Learned units of a single block with fixed α in 3D environment. Every row shows the learned units with a given α.

Figure 15: Examples of 3D path integral with duration T = 30.

Figure 16: Examples of 3D simple path planning, where the agent is capable of planning a direct trajectory.

Figure 17: Examples of 3D path planning with a cuboid obstacle.

Figure 18: Response curves of learned units in 1D. Each block shows the units belonging to the same subvector in the motion model. The horizontal axis represents the time points, while the vertical axis indicates the value of responses.

Error correction results on the vector representation. The performance of path integral is measured by mean square error between predicted locations and ground truth locations; while for path planning, the performance is measured by success rate. Experiments are conducted using several noise levels: Gaussian noise with different standard deviations in terms of the reference standard deviation s and dropout mask with different percentages. DE means implementing decoding-encoding process when performing the tasks.

Table content: <table><tr><td/><td/><td colspan="3">Path integral: MSE</td><td/><td/><td>Path planning: success rate</td></tr><tr><td>Noise type</td><td>1s</td><td>0.75s</td><td>0.5s</td><td colspan="2">0.25s 0.1s</td><td>1s</td><td colspan="2">0.75s 0.5s 0.25s 0.1s</td></tr><tr><td colspan="5">Gaussian (DE) 1.687 1.135 0.384 0.017</td><td>0</td><td colspan="3">0.928 0.959 0.961 0.977 0.985</td></tr><tr><td>Gaussian</td><td colspan="8">6.578 2.999 1.603 0.549 0.250 0.503 0.791 0.934 0.966 0.982</td></tr><tr><td>Noise type</td><td>70%</td><td>50%</td><td>30%</td><td>10%</td><td>5%</td><td colspan="2">70% 50% 30% 10%</td><td>5%</td></tr><tr><td colspan="9">Dropout (DE) 2.837 1.920 1.102 0.109 0.013 0.810 0.916 0.961 0.970 0.978</td></tr><tr><td>Dropout</td><td colspan="8">19.611 16.883 14.137 3.416 0.602 0.067 0.186 0.603 0.952 0.964</td></tr></table>

Peer review (iclr2019):
 Updated score from 6 to 7 after the authors addressed my comments below. Previous review: This paper builds upon the recent work on computational models of grid cells that rely on trainable (parametric) models such as recurrent neural networks [Banino et al, 2018; Cueva & Wei, 2018]. It focuses entirely on path integration in 2D and 3D, from velocity inputs only, and it relies on two sub-networks: the motion model (an RNN) and the localization model (a feed-forward network). The avowed goal of the paper is to build a very simple and linear model for grid cells. By linearly embedding the position x into a high-dimensional hidden vector v(x) (e.g., 96 elements), it can model motion using a linear model relying on matrix-vector multiplication: v(x + dx) = M(dx) v(x), where dx is the 2D or 3D displacement, v(.) is a vector and M(.) is a matrix. The embeddings v(.) are learnable and the paper assumes a square or cubic grid of N*N or N*N*N possible positions x (with N=40); these embeddings are also normalized to unit length and obey the kernel constraint that the dot-product between any two positions* vectors v(x) and v(y) is a Gaussian or an exponential function. The motion matrix is represented as block diagonal, where each block is a rotation of subvector v_k(x) into v_k(x + dx), where each block corresponds to a specific grid cell, and where the diagonal block is further expressed as a quadratic function of dx_1, dx_2, dx_3 elements of the displacement vector. The strengths of the paper are that: 1) The supervision of the localization subnetwork only depends on Euclidean proximity between two positions x and y and therefore uses relative positions, not absolute ones. Similarly, the path integration supervision of the motion model uses only relative displacements. 2) The resulting rate maps of the hidden units seem perfect; the model exhibits multi-scale grid behaviour. 3) The idea of using disentangled blocks, rather than injecting noise or using dropout and a softmax bottleneck as in [Banino et al, 2018], is interesting. 4) The model accumulates little path integration error over 1000 step-long episodes. The weakness of the paper is its simplicity: 1) The assumption that A(x, y) can be modeled by a Gaussian or exponential (Laplacian?) kernel is limiting, in particular for positions x and y that are far apart. 2) There is no discussion about what egocentric vs. allocentric referentials, and dx is assumed to be aligned with (x, y) axes (which are also the axes defining the bounding box of the area). 3) Unlike the other work on learning path integration using an RNN, the linear matrix model can only handle allocentric displacements dx_1, dx_2 (and optional dx_3 in 3D). 4) No consideration is given to non-square areas: would the network also exhibit grid-like behavior if the area was circular? 5) What happens if the quadratic parameterisation of block diagonals is dropped? 6) The paper did not use metrics accepted in the neuroscience community for computing a gridness score of the grid cells (although the grid cell nature is evident). There should however be metrics for quantifying how many units represent the different scales, offsets and orientations. 7) The authors did not consider (but mentioned) embedding locations from vision, and did not consider ambiguous position embeddings. The experiments about path planning are unconvincing. First of all, the algorithm requires to input absolute positions of every obstacle into equation (9) - (10), which assumes that there is perfect information about the map. Secondly, the search algorithm is greedy and it is not obvious how it would handle a complex maze with cul-de-sac. Saying that "there is no need for reinforcement learning or sophisticated optimal control" is very misleading: the problem here is simplified to the extreme, and fully observed, and any comparison with deep RL algorithms that can handle partial observations is just out of place. In summary, the authors have introduced an interesting and elegant model for grid cells that suffers from simplifications. The part on path planning should be cut and replaced with more analysis of the grid cells and an explanation of how the model would handle egocentric velocity.