Year: 2019

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.

heading: 1 INTRODUCTION
text: 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.
Figure 1(a) shows Dr. May-Britt Moser, who together with Dr. Edvard Moser, won the 2014 Nobel Prize for Physiology or Medicine, for their discovery of the grid cells (Hafting et al. (2005); Fyhn et al. (2008); Yartsev et al. (2011); Killian et al. (2012); Jacobs et al. (2013); Doeller et al. (2010))
∗Equal contributions.
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.
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.

heading: 2 CONTRIBUTIONS AND RELATED WORK
text: 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)).

heading: 3 REPRESENTATIONAL MODEL OF GRID CELLS
text: Consider an agent navigating within a domainD = [0, 1]× [0, 1] (actually the shape does not matter, and it can be R2). We can discretize D into an N × N lattice. N = 40 in our experiments. Let x = (x1, x2) ∈ 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.

heading: 3.1 SUB-MODEL 1 ABOUT MOTION ALGEBRA: VECTOR-MATRIX MULTIPLICATION
text: 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:
Motion : xt +∆x −−−−−→ xt+1
↓ ↓ ↓ v(xt) M(∆x)× −−−−−→ v(xt+1)
(2)
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 M = I + M̃(∆x) such that each element of M̃(∆x) is a quadratic (or polynomial) function of ∆x = (∆x1,∆x2):
M̃ij(∆x) = β (1) ij ∆x1 + β (2) ij ∆x2 + β (11) ij ∆x 2 1 + β (22) ij ∆x 2 2 + β (12) ij ∆x1∆x2, (3)
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.

heading: 3.2 DISENTANGLED BLOCKS OR MODULES
text: 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.

heading: 3.3 SUB-MODEL 2 ABOUT LOCAL GEOMETRY: MAGNIFIED LOCAL ISOMETRY
text: 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,
〈v(k)(x), v(k)(x+ ∆x)〉 = d(1− αk|∆x|2), (4) 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.

heading: 3.4 ROTATION AND PROJECTION
text: 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 maxx〈v(k), v(k)(x)〉 where the maximization is within the local neighborhood.

heading: 3.5 ERROR CORRECTION
text: 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, s2Id), and Id 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.

heading: 3.6 HEXAGON GRID PATTERNS
text: 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) = (ei〈aj ,x〉, j = 1, 2, 3)>, and a1, a2, a3 are three 2D vectors so that the angle between ai and aj is 2π/3 for ∀i 6= j and |aj | = 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 ai and aj is 2π/3 for ∀i 6= j, patterns of v(x) over x have hexagon periodicity. Moreover, the scale of the patterns is controlled by the length of aj , 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.

heading: 3.7 SUB-MODEL 3 ABOUT GLOBAL GEOMETRY: ADJACENCY KERNEL
text: Because of the periodicity, each block (v(k)(x),∀x) does not form a global codebook of 2D positions, i.e., there can be x 6= y, but v(k)(x) = v(k)(y), i.e., v(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:
〈v(x), v(y)〉 = K∑ k=1 〈v(k)(x), v(k)(y)〉 = (Kd)f(|x− y|), (6)
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 ( −r2/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 maxv∈C〈u, v〉. The whole vector v is like aK-arm clock, with each v(k) being an arm rotating at a speed ωk|∆x| =√
2αk|∆x|.

heading: 3.8 LOCALIZATION AND HEAT MAP
text: (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 maxx〈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)×N2 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 N2-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 N2-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:
Localization : v V >× −−−−−→ h (heat map and decoding to δx)
δx V× −−−−−→ v(x) (encoding) (8)
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.

heading: 3.9 PATH INTEGRAL
text: 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 x0 and motion sequences {∆x1, ...,∆xT }, and the output is the prediction of one’s current position xT . We first encode the initial position x0 as v(x0). Then, by the motion model, the hidden vector v(xT ) at time T can be predicted as:
v(xT ) = ∏1
t=T M(∆xt)v(x0). (9)
We can then decode xT from v(xT ).

heading: 3.10 PATH PLANNING
text: Our representation system can plan direct path from the starting position x0 to the target position y by steepest ascent on the inner product 〈v(x), v(y)〉, i.e., let v0 = v(x0), the algorithm iterates
∆xt = arg max ∆x∈∆
〈v(y),M(∆x)vt−1〉, (10)
vt = M(∆xt)vt−1, (11)
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).

heading: 3.11 PLACE CELLS AND SCALE
text: 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(−r2/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.

heading: 3.12 GROUP REPRESENTATION IN MOTOR CORTEX
text: 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.

heading: 4 LEARNING REPRESENTATION
text: 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 L1 = Ex,∆x [ ‖v(x+ ∆x)−M(∆x)v(x)‖2 ] , (12)
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
L1,T = Ex,∆x1,...∆xT [ ‖v(x+ ∆x1 + ...+ ∆xT )−M(∆xT )...M(∆x1)v(x)‖2 ] , (13)
where (∆xt, 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 L2,k = Ex,∆x [ (〈v(k)(x), v(k)(x+ ∆x)〉 − (d(1− αk|∆x|2))2 ] , (14)
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 L2 = ∑K k=1 L2,k.
For sub-model (3) on global adjacency kernel, the loss function is L3 = Ex,y [ ((Kd)f(|x− y|)− 〈v(x), v(y)〉)2 ] , (15)
where both x and y are sampled uniformly from D. Let v(x) = (vi(x), i = 1, ...,Kd), we also impose a regularization loss L0 = ∑Kd i=1(Ex[vi(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 Ex[vi(x)] within a certain range is biologically plausible, because vi(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 Ex[vi(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 Ex[vi(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,∆xt, t = 1, ..., T )) for L1, (x,∆x) for L2, and (x, y) for L3 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.

heading: 5 EXPERIMENTS
text: 5.1 LEARNING SINGLE BLOCKS: HEXAGON PATTERNS AND METRICS
We first learn a single block with fixed αk by minimizing L1 + λ2L2,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. Figure 3(a) shows the learned results with 6 units and different values of α. The scale or metric of the lattice is controlled by α. The units within a block have patterns with similar scale and arrangement, yet different phases. Figure 3(b) shows the learned results with 100 units and α = 72, indicating that the grid-like patterns are stable and easy to learn even when the number of units is large.

heading: 5.2 LEARNING MULTIPLE HEXAGON BLOCKS AND METRICS
text: We learn multiple blocks by minimizing L1 + λ2L2 + λ3L3. Instead of manually assigning αk, we learn αk by gradient descent, simultaneously with v and M .
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 L1 +λ3L3 for simplicity.

heading: 5.3 PATH INTEGRAL
text: Figure 5(a) shows an example of path integral result (time duration T̃ = 40), where we use single step motion loss L1,T=1 in learning. Gaussian kernel with σ = 0.08 is used as the adjacency measure in L3. 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 L1,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).

heading: 5.4 PATH PLANNING
text: 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 ∆ = {∆xi = (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∆xi) = M
k(∆xi). 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 ∆xt at time t by
∆xt = arg max ∆x∈∆
[ 〈v(y),M(∆x)vt〉 − a〈v(z),M(∆x)vt〉b ] , (16)
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 {zi}mi=1. They can be included by summing over the kernels of every obstacle {a〈v(zi),M(∆x)vt〉b}mi=1 and choosing ∆xt at time t by ∆xt = arg max∆x∈∆ [ 〈v(y),M(∆x)vt〉 − ∑m i=1 a〈v(zi),M(∆x)vt〉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)).

heading: 6 DISCUSSION: ROTATIONIST-CONNECTIONIST MODEL?
text: In terms of general modeling methodology, a typical recurrent network is of the form vt = tanh(W (vt−1, δt)), where vt is the latent vector, δt is the input change or action. vt−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 vt = M(δt)vt−1, where M(δt) is a matrix that is non-linear in δt. M(δt) is a matrix representation of δt. vt 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

heading: ACKNOWLEDGEMENT
text: We thank the three reviewers for their insightful comments and suggestions. Part of the work was done while Ruiqi Gao was an intern at Hikvision Research Institute during the summer of 2018. She thanks Director Jane Chen for her help and guidance. We also thank Jiayu Wu for her help with experiments and Zilong Zheng for his help with visualization. The work is supported by DARPA XAI project N66001-17-2-4029; ARO project W911NF1810296; ONR MURI project N00014-161-2007; and a Hikvision gift to UCLA. We gratefully acknowledge the support of NVIDIA Corporation with the donation of the Titan Xp GPU used for this research.

heading: A PROOF FOR SECTION 3.6
text: Proof. (aj , j = 1, 2, 3) forms a tight frame in the 2D space, in that for any vector x in 2D,∑3 j=1〈x, aj〉2 ∝ |x|2.
Since C∗C = I , we have
〈v(x), v(y)〉 = v(x)∗v(y) (17) = e(x)∗C∗Ce(y) (18)
= 3∑ j=1 ei〈aj ,y−x〉. (19)
Then we have
RE(〈v(x), v(x+ ∆x)〉) = 3∑ j=1 cos(〈aj ,∆x〉) (20)
≈ 3∑ j=1 (1− 〈aj ,∆x〉2/2) (21) = 3− 3∑ j=1 〈aj ,∆x〉2/2 (22)
= 3(1− α|∆x|2), (23) where |aj | = 2 √ α.
The self motion from v(x) to v(x+ ∆x) is
v(x+ ∆x) = Ce(x+ ∆x) (24) = CD(∆x)e(x) (25) = CD(∆x)C∗v(x) (26) = M(∆x)v(x), (27)
where D(∆x) = diag(ei〈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
f(|x− y|) = 〈v(x), v(y)〉 (28)
≈ K∑ k=1 3∑ j=1 ei〈akj ,y−x〉 (29)
= K∑ k=1 〈v(k)(x), v(k)(y)〉. (30)

heading: B HEXAGON GRID PATTERNS AND METRICS

heading: B.1 SIMULATED INPUT DATA
text: 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 x0. At each time step t, a small motion ∆xt (≤ 3 grids in each direction) is randomly sampled with the restriction of not leading the agent outside the boundary, resulting in a simulated trajectory {x0 + ∑t i=1 ∆xi}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 L3, while consecutive position sequences (x,∆xt, t = 1, ..., T ) are randomly sampled as the input to the motion loss L1, with length specified by T (which is usually much smaller than the whole length of the trajectory T̃ ).

heading: B.2 LEARNED SINGLE BLOCK UNITS WITH DIFFERENT BLOCK SIZES
text: 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.
Figure 7 displays the response maps of the learned single block units with different block sizes. For block sizes 4 and 5, the learned maps show square lattice patterns. For block sizes greater than or equal to 6, the learned maps show hexagon lattice patterns.

heading: B.3 LEARNED MULTIPLE BLOCK UNITS AND METRICS

heading: B.3.1 WITH DIFFERENT BLOCK SIZES
text: Figure 8 displays the response maps of the multiple block units with different block sizes. The metrics of the multiple blocks are learned automatically.

heading: B.3.2 WITH DIFFERENT SHAPES OF AREA
text: Figure 9 displays the response maps of the multiple block units with different shapes of the area, such as circle and triangle.

heading: B.4 QUANTITATIVE ANALYSIS OF SPATIAL ACTIVITY
text: 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
αk Learned blocks
5.7
7.7
10.4
10.7
11.6
12.7
18.0
31.2
31.2
37.8
56.0
68.0
73.5
74.3
86.5
92.5
(a) Block size = 2
αk Learned blocks
4.6
6.7
10.6
11.0
13.9
23.4
31.3
38.9
39.3
56.7
64.5
72.3
75.4
80.5
84.4
93.0
(b) Block size = 3
αk Learned blocks
2.0
6.9
12.3
18.7
19.5
25.7
28.9
39.7
41.9
52.3
55.1
61.1
66.3
67.7
70.6
88.3
(c) Block size = 4
αk Learned blocks
1.8
9.1
12.4
18.5
28.6
32.8
41.9
48.6
49.9
58.0
64.5
73.0
74.9
77.8
81.7
97.1
(d) Block size = 5
remaining similar within each block. We show the histograms of the grid orientations and scales in Figure 10(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 figure 10(d). Interestingly,
the grid scale is nicely proportional to the learned 1/ √ αk.

heading: B.5 ABLATION STUDIES
text: We conduct ablation studies to assess various assumptions in our model.

heading: B.5.1 LOSS TERMS
text: We learn v and M with multiple blocks by minimizing L1 + λ2L2 + λ3L3, which consists of (1) a motion loss L1, (2) a local isometry loss L2, and (3) a global adjacency loss L3. 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 L3, or only the isometry loss L2. 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.

heading: B.5.2 ASSUMPTIONS OF MOTION MATRIX
text: 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 L2 or the global adjacency loss L3 in addition to the motion loss L1. 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.

heading: C MODELING EGOCENTRIC MOTION
text: The model can be generalized to handle egocentric motion that involves head direction.

heading: C.1 COUPLING TWO GRID SYSTEMS
text: The egocentric motion consists of angular velocity in the change of head direction and the spatial velocity along the current head direction. We couple two grid systems, one for head direction and the other for the spatial position. For notational simplicity, we put the “hat” notation on top of v and M notation to denote the vector and matrix for the head direction.
Specifically, the agent has a head direction θ, which may change over time. The agent moves along its head direction with scalar motion. We can discretize the range of θ, [0, 2π], into N̂ equally spaced values {θi, i = 1, ..., N̂}, and introduce a vector representation of self-direction. That is, the agent represents its self-direction by a d̂-dimensional hidden vector v̂(θ).

heading: C.1.1 MOTION SUB-MODEL
text: Suppose at a position x, the agent has a head direction θ. The self-motion is decomposed into (1) a scalar motion δ along the head direction θ and then (2) a head direction rotation ∆θ. The self-motion ∆x = (δ cos θ, δ sin θ). We assume
v̂(θ + ∆θ) = M̂(∆θ)v̂(θ), (31) v(x+ ∆x) = M(δ, v̂(θ))v(x), (32)
where M̂(∆θ) is a d̂ × d̂ matrix that depends on ∆θ, and M(∆θ, v̂(θ)) is a d × d matrix that depends on δ and v̂(θ). M̂(∆θ) is the matrix representation of the head direction rotation ∆θ, while M(δ, v̂(θ)) is the matrix representation of scalar motion δ along the direction θ.
We can model M(δ, v̂(θ)) by an attention (or selection) mechanism:
pi = 〈v̂(θ), v̂(θi〉b∑N̂ i=1〈v̂(θ), v̂(θi)〉b , (33)
M(δ, v̂(θ)) = ∑N̂
i=1 piM
(i)(δ), (34)
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 = (pi, 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.

heading: C.1.2 LOCALIZATION SUB-MODEL
text: 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.

heading: C.1.3 LOSS FUNCTION FOR LEARNING
text: For learning the system, we can first learn v̂ and M̂ (or the coefficients that parametrize M̂ ) by minimizing
Eθ1,θ2 [ (f̂(|θ1 − θ2|)− 〈v̂(θ1), v̂(θ2)〉)2 ] + λEθ,∆θ [ ‖v̂(θ + ∆θ)− M̂(∆θ)v̂(θ)‖2 ] , (35)
and then we learn v and M (or the coefficients that parametrize M ) as before.

heading: C.2 LEARNED UNITS FOR SELF-DIRECTION AND SELF-POSITION
text: 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.

heading: C.3 CLOCK AND TIMESTAMP
text: 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).

heading: D ERRORS

heading: D.1 ERROR CORRECTION
text: 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 vt, we decode it to xt = arg maxx〈vt, v(x)〉 and then re-encode it to the hidden vector v(xt). Actually the whole process can be accomplished by projection to the codebook sub-manifold without explicit decoding, by obtaining the vector: arg maxv∈C〈vt, 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.

heading: D.2 NOISY SELF-MOTION INPUT
text: 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 L1+λ2L2,k. A batch of 500,000 examples of (x, y) and (x,∆xt, t = 1, ..., T ) is sampled online at every iteration for training. Figure 14 shows the learned units.
Next we learn multiple blocks by minimizing L1 +λ3L3, and use the learned models to perform 3D path integral and path planning. For simplicity, we remove L2. 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 mn2 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.

heading: E.4 LEARNING IN 1D
text: 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,∆t2). 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.

Review summaries:
This paper proposes a simple and elegant approach to learning "grid-cell like" representations that uses a high-dimensional encoding of position, together with a matrix for propagating position that involves only local connections among the elements of the vector.  The vectors are also constrained to have their inner products reflect positional similarity.  The paper also shows how such a representation may be used for path planning.