Vector-valued Gaussian Processes on Riemannian Manifolds via Gauge Independant Projected Kernels
Year: 2021

Abstract:
Gaussian processes are machine learning models capable of learning unknown functions in a way that represents uncertainty, thereby facilitating construction of optimal decision-making systems. Motivated by a desire to deploy Gaussian processes in novel areas of science, a rapidly-growing line of research has focused on constructively extending these models to handle non-Euclidean domains, including Riemannian manifolds, such as spheres and tori. We propose techniques that generalize this class to model vector fields on Riemannian manifolds, which are important in a number of application areas in the physical sciences. To do so, we present a general recipe for constructing gauge independent kernels, which induce Gaussian vector fields, i.e. vector-valued Gaussian processes coherent with geometry, from scalar-valued Riemannian kernels. We extend standard Gaussian process training methods, such as variational inference, to this setting. This enables vector-valued Gaussian processes on Riemannian manifolds to be trained using standard methods and makes them accessible to machine learning practitioners.

heading: 1 Introduction
text: Gaussian processes are an effective model class for learning unknown functions. They are particularly attractive for use within data-efficient decision systems, including Bayesian optimization [3, 32, 39], model-based reinforcement learning [34, 7], and active learning [21]. In these settings, Gaussian processes can represent and propagate uncertainty, as well as encode inductive biases as prior information in order to drive data efficiency. A key aspect of prior information is the geometry of the domain on which the Gaussian process is defined, which often encodes key properties, such as symmetry. Following the growing deployment of Gaussian processes, a number of recent works have focused on how to define Gaussian processes on non-Euclidean domains in ways that reflect their geometric structure [2, 1].
In many applications, such as climate science, quantities of interest are vector-valued. For example, global wind velocity modeling must take into account both speed and direction, and is represented by a vector field. On geometric domains, the mathematical properties of vector fields can differ noticeably from their Euclidean counterparts: for instance, one can prove that every smooth vector field on a sphere must vanish in at least one point [26]. Behavior such as this simultaneously highlights the need to represent geometry correctly when modeling vector-valued data, and presents a number of non-trivial technical challenges in constructing models that are mathematically sound.
∗Equal contribution. Code: https://github.com/MJHutchinson/ExtrinsicGaugeIndependantVectorGPs. For a general implementation, see https://github.com/GPflow/GeometricKernels/.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
In particular, even the classical definition of a vector-valued Gaussian process—that is, a random function with multivariate Gaussian marginals at any finite set of points—already fails to be a fully satisfactory notion when considering smooth vector fields on a sphere. This is because tangent vectors at distinct points live within different tangent spaces, and it is not clear how to construct a cross-covariance between them that does not depend on a completely arbitrary choice of basis vectors within each space. Constructions that are independent of this choice of basis are called gauge independent, and recent work [45, 5, 15] in geometric machine learning has focused on satisfying this key property for convolutional neural networks that deal with non-Euclidean data.
Our contributions include the following. We (a) present a differential-geometric formalism for defining Gaussian vector fields on manifolds in a coordinate-free way, suitable for Gaussian process practitioners with minimal familiarity with differential geometry, (b) present a universal and fully constructive technique for defining prior Gaussian vector fields on Riemannian manifolds, which we term the projected kernel construction, and (c) discuss how to adapt key components in the computational Gaussian process toolkit, such as inducing point methods, to the vector field setting.
The structure of the paper is as follows. In Section 2, we define vector-valued Gaussian processes on smooth manifolds. We start by reviewing the multi-output Gaussian process set-up, which is typically used in machine learning. We then detail a differential-geometric formalism for defining vectorvalued Gaussian processes on smooth manifolds. In Section 3, we provide a concrete construction for these Gaussian processes on Riemannian manifolds and discuss how they can be trained using variational sparse approximations. Section 4 showcases Gaussian vector fields on two tasks, namely weather imputation from satellite observations and learning the dynamics of a mechanical system.

heading: 2 Vector-valued Gaussian Processes on Smooth Manifolds
text: A vector-valued Gaussian process (GP) is a random function f : X → Rd such that, for any finite set of points x ∈ Xn, the random variable f(x) ∈ Rn×d is jointly Gaussian. Every such GP is characterized by its mean function µ : X → Rd and matrix-valued covariance kernel k : X ×X → Rd×d, which is a positive-definite function in the matrix sense. These functions satisfy E(f(x)) = µ(x) and Cov(f(x),f(x′)) = k(x,x′) for any x,x′ ∈ X . Here, dependence between function values is encoded in the kernel’s variability along its input domain, and correlations between different dimensions of the vector-valued output are encoded in the matrix that the kernel outputs.
Consider a function f with y = f(x) + ε, where ε ∼ N(0, σ2I) and training data (x,y). Placing a GP prior f ∼ GP(0, k) on the unknown function results in a GP posterior, whose mean and covariance are given by
E(f | y) = K(·)x(Kxx + σ2I)−1y Cov(f | y) = K(·,·) −K(·)x(Kxx + σ2I)−1Kx(·). (1) Here, (·) denotes an arbitrary set of test locations, Kxx = k(x,x) is the kernel matrix, and K(·)x = k(x, ·) is the cross-covariance matrix between function values evaluated at the training and test inputs. The GP posterior can also be written as
(f | y)(·) = f(·) + K(·)x(Kxx + σ2I)−1(y − f(x)− ε), ε ∼ N(0, σ2I) (2) where f(·) is the prior GP, and equality holds in distribution [47, 48]. These expressions form the foundation upon which Gaussian-process-based methods in machine learning are built.
Recent works have studied techniques for working with the expressions (1) and (2) when the input domain X is a Riemannian manifold, focusing both on defining general classes of kernels [2], and on efficient computational techniques [47, 48]. In this setting, namely for f : X → R, defining kernels already presents technical challenges: the seemingly-obvious first choice one might consider, namely the geodesic squared exponential kernel, is ill-defined in general [13]. We build on these recent developments to model vector fields on manifolds using GPs. We do not consider manifold-valued generalizations of Gaussian processes, for instance f : R→ X: various constructions in this setting are instead studied by Stroock [40], Émery [11], Mallasto and Feragen [29], and Mallasto et al. [30]. To begin, we review what a vector field on a manifold actually is.

heading: 2.1 Vector Fields on Manifolds
text: Let X be a d-dimensional smooth manifold with TxX denoting its tangent space at x. Let TX = {(x, v) |x ∈ X, v ∈ TxX} be its tangent bundle, and let T ∗X = {(x, φ) |x ∈ X,φ ∈ T ∗xX}
Basis coefficients Frame Vector field
be its cotangent bundle—endow both spaces with the structure of smooth manifolds. Define the projection map projX : TX → X by projX(x, v) = x. A vector field on X is a map that assigns each point in X to a tangent vector attached to that point. More formally, a vector field is a cross-section, or simply a section, of the tangent bundle, which is a map f : X → TX , such that projX f(x) = x for all x.1 A vector field is called smooth if this map f is smooth.
To represent a vector field on a manifold numerically, one must choose a basis in each tangent space, which serves as a coordinate system for vectors in the tangent space. On many manifolds it is impossible to choose these basis vectors in a way that they vary smoothly in space.2 This can be handled by working with local coordinates, or with bases that are non-smooth. Any chosen set of basis vectors is arbitrary, so objects constructed using them should not depend on this choice. Constructions that satisfy this property are called gauge independent. This notion is illustrated in Figure 1, and will play a key role in the sequel.

heading: 2.2 Gaussian Vector Fields
text: Upon reflecting on the above considerations in the context of GPs, the first issue one encounters is that, for a random vector field f : X → TX , it is not clear what it means for finite-dimensional marginal distributions to be multivariate Gaussian given that f takes its values in a bundle rather than a vector space. The first step towards constructing Gaussian vector fields, therefore, involves adapting the notion of finite-dimensional marginal distributions appropriately.
Definition 1. Let X be a smooth manifold. We say that a random vector field f is Gaussian if for any finite set of locations x1, . . . , xn ∈ X , the random vector (f(x1), . . . , f(xn)) ∈ Tx1X⊕ . . .⊕TxnX is Gaussian (either in the sense of duality or in any basis: see Appendix A for details).
This definition is near-identical to the Euclidean case: the only difference is that finite-dimensional marginals are now supported in a direct sum of tangent spaces, instead of Rn×d. With this definition, the standard multi-output GP properties, such as conditioning, carry over, virtually unmodified. Definition 1 is a natural choice: if we embed our manifold into Euclidean space, the induced GP is a vector-valued GP as defined in the beginning of Section 2.
1A vector field is not the same as a map f̃ : X → Rd: an output value f(x) ∈ TX formally consists of both a copy of the input point x, and vector within the tangent space TxX at this point. This encodes the geometric structure of the underlying manifold. The algebraic requirement projX f(x) = x for all x ensures that the tangent vector chosen correctly corresponds to the point at which it is attached.
2If a smooth choice of basis vectors existed, it would define a smooth non-vanishing vector field. On the sphere, by the hairy ball theorem, all smooth vector fields vanish in at least one point, so no such bases exist.
Proposition 2. Let X be a manifold and emb : X → Rp be a smooth embedding. Let f be a Gaussian vector field (as defined in Definition 1), and let femb : emb(X)→ Rp be its pushforward along the embedding. Then femb is a vector-valued Gaussian process in the Euclidean sense.
All proofs in this work can be found in Appendix A. Having established the notion of a vector-valued Gaussian process on a smooth manifold, we proceed to deduce what mathematical objects play the role of a mean function and kernel, so that it is clear what ingredients are needed to construct and determine such a process.
The former is clear: the mean of a Gaussian vector field should be an ordinary vector field, and will determine the mean vector at all finite-dimensional marginals. The kernel, on the other hand, is less obvious: because distinct tangent vectors live in different tanget spaces, it is unclear whether or not a Gaussian vector field is characterized by an appropriate notion of a matrix-valued kernel, or by something else. Now, we define the right notion for kernel in this setting. Definition 3. We call a symmetric function k : T ∗X × T ∗X → R fiberwise bilinear if for all pairs of points x, x′ ∈ X
k(λαx + µβx, γx′) = λk(αx, γx′) + µk(βx, γx′) (3)
holds for any αx, βx ∈ T ∗xX , γx′ ∈ T ∗x′X and λ, µ ∈ R, and positive semi-definite if for any set of covectors αx1 , . . . , αxn ∈ T ∗X , we have ∑n i=1 ∑n j=1 k(αxi , αxj ) ≥ 0. We call a symmetric fiberwise bilinear positive semi-definite function a cross-covariance kernel.
This coordinate-free function should be viewed as analogous to ((x,v), (x′,v′)) 7→ vTKx,x′v′ in the Euclidean setting, where v,v′ multiply the matrix-valued kernel from both sides. Its coordinate representation, which more closely matches the Euclidean case, will be explored in the sequel. To show that this is indeed the right notion, we prove the following result. Theorem 4. The system of marginal distributions of a Gaussian vector field on a smooth manifold X is uniquely determined by a mean vector field µ : X → TX and a cross-covariance kernel k : T ∗X × T ∗X → R. Moreover, this correspondence is one-to-one.
By virtue of defining and characterizing all Gaussian vector fields, Theorem 4 assures us the definition of a kernel introduced is the correct mathematical notion. The constructions presented here are all intrinsic or, in other words, coordinate-free, and do not involve the use of bases. To understand how to perform numerical calculations with these kernels we proceed to study their coordinate representations with respect to a specific choice of basis.

heading: 2.3 Matrix-valued Kernels
text: In Section 2.2, we defined what a Gaussian vector field on a manifold is. However, by nature of the manifold setting, the resulting objects are more abstract than usual and do not describe how it can be represented numerically. We now develop a point of view suitable for this task.
To this end, we introduce a frame F on X , also known as a gauge in physical literature, which is a collection of (not necessarily smooth) vector fields e1, . . . , ed on X such that at each point x ∈ X , the set of vectors e1(x), . . . , ed(x) forms a basis of TxX . The frame allows us to express a vector field f on X as simply a vector-valued function f = (f1, . . . , fd) : X → Rd, such that f(x) = ∑d i=1 f
i(x)ei(x) for all x ∈ X . The corresponding coframe F ∗ is defined as a collection e1, . . . , ed of covector fields (one-forms) on X such that 〈 ei(x)|ej(x) 〉 = δij for all x ∈ X , where δij is the Kronecker delta. In the following proposition, we show that if f is a Gaussian vector field on X (in the sense of Definition 1), then the corresponding vector representation f expressed in a given frame is a vector-valued GP in the standard sense. Proposition 5. Let f be a Gaussian vector field defined on X with cross-covariance kernel k : T ∗X × T ∗X → R. Given a frame F = (e1, . . . , ed) on X , define f : X → Rd as above. Then f is a vector-valued GP in the usual sense with kernel KF : X ×X → Rd×d given by
KF (x, x ′) = k(e 1(x), e1(x′)) . . . k(e1(x), ed(x′)) ... . . .
... k(ed(x), e1(x′)) . . . k(ed(x), ed(x′))
 , (4)
where (ei), with raised indices, is the coframe corresponding to (ei). Conversely, given a vectorvalued GP f = (f1, . . . , fd) : X → Rd and a frame F = (e1, . . . , ed) on X , f(·) :=∑d i=1 f i(·)ei(·) defines a Gaussian vector field on X .
This result shows precisely how numerical representations of a Gaussian vector field depends on the choice of frame. While this representation is not invariant under this choice, it is equivariant, meaning that a transformation in the frame results in an appropriate transformation of the kernel. To make this notion precise, we introduce a matrix subgroup G = GL(d,R), called the gauge group, that acts on Rd by a standard matrix-vector multiplication. Given two frames F, F ′ on X , an abstract vector fx ∈ TxX has two vector representations fx,f ′ x in the respective frames. We say that F
′ is obtained from F by a gauge transformation with respect to a matrix field A : X → G ⊆ Rd×d, if
f ′x = A(x)fx (5)
holds for all x ∈ X , and we write F ′ = AF . Note that A(x) need not be smooth in x. We see that the gauge transformation is therefore just a linear change of basis of the frame F at each point for which one can identify vectors in TxX as elements in Rd. The corresponding gauge dependant matrix-valued kernels must also respect this transformation rule, a statement of gauge independence. Corollary 6. Let F be a frame on X and KF : X × X → Rd×d be the corresponding matrix representation (4) of a cross-covariance kernel k : T ∗X ×T ∗X → R. This satisfies the equivariance condition
KAF (x, x ′) = A(x)KF (x, x ′)A(x′)T , (6) where A : X → G ⊆ Rd×d is a gauge transformation. All cross-covariance kernels in the sense of Proposition 4 arise this way.
Hence, one way to define a Gaussian vector field on a manifold is to find a gauge independent kernel. In summary, we have described Gaussian vector fields in a coordinate-free differential-geometric language, and deduced enough properties to confirm the objects defined truly deserve to be called GPs. In doing so, we have both introduced the necessary formalism to the GP community, and obtained a recipe for defining kernels, through a simple condition atop standard matrix-valued kernels. To proceed towards practical machine learning methods, we therefore study techniques for constructing such kernels explicitly.

heading: 3 Model Construction and Bayesian Learning for Riemannian Manifolds
text: In Section 2, we introduced a notion of a Gaussian vector field. We now study how to use vector fields for machine learning purposes. This entails two primary issues: (a) how to construct practically useful kernels, and (b) once a kernel is constructed, how to train Gaussian processes.
To construct a Gaussian vector field prior, the preceding theory tells us that we need to specify a mean vector field and a cross-covariance kernel. From the definition, it is not at all obvious how to specify a natural kernel, and experience with the scalar-valued case—where the innocuous-looking geodesic squared exponential kernel is generally not positive semi-definite on most Riemannian manifolds [13]—suggests that the problem is delicate, i.e., simply guessing the kernel’s form is unlikely to succeed. Our goal, therefore, is to introduce a general construction for building wide classes of kernels from simple building blocks.
The same issues are present if we consider variational approximations to posterior GPs, such as the inducing point framework of Titsias [43]: these are formulated using matrix-vector expressions involving kernel matrices, and it is important for the approximate posterior covariance to be gauge independent in order to lead to a valid approximate process. We proceed to address these issues.

heading: 3.1 Projected Kernels
text: Here, we introduce a general technique for defining cross-covariance kernels k : T ∗X × T ∗X → R and for working with such functions numerically. Section 2 gives us a promising strategy to construct a suitable kernel—namely, it suffices to find a gauge independent matrix-valued kernel. At first glance, it is not obvious how to construct such a kernel in the manifold setting. On many manifolds, such as the sphere, owing to the hairy ball theorem, every frame must be discontinuous: therefore, constructing a continuous kernel in such a choice of frame appears difficult.
Scalar processes Embedded process Projected process
To both get around these obstacles, and aid numerical implementation, we propose to isometrically embed the manifold into Euclidean space.3 Doing so greatly simplifies these issues by virtue of making it possible to represent the manifold using a single global coordinate system. On the other hand, the main trade-off from this choice is that by its extrinsic nature, the construction can make theoretical analysis more difficult. To proceed, we need two ingredients.
1. An isometric embedding emb : X → Rd′ of the manifold. 2. A vector-valued Gaussian process f ′ : X → Rd′ in the standard sense.4
A simple choice which reflects the geometry of the manifold is to take f ′ to be d′ independent scalar-valued GPs on X .
By standard results in differential geometry, any smooth map φ : X → X ′ between two manifolds induces a corresponding linear map on the tangent spaces dxφ : TxX → Tφ(x)X ′, which can loosely be thought of as mapping φ to its first-order Taylor expansion at x. Thus, an embedding emb : X → Rd′ , induces a map dxemb : TxX → Temb(x)Rd ′ . Now fixing a frame F on X , each tangent space TxX can be identified with Rd, so without loss of generality, the map dxemb can be expressed simply as a position-dependent matrix PTx ∈ Rd
′×d. Taking the transpose, we obtain Px ∈ Rd×d ′ , which we call the projection matrix. The desired Gaussian vector field on X , with respect to F , is then constructed as f(x) = Pxf ′(x). This procedure is illustrated in Figure 2: there we see that to get a vector field on an R3-embedded sphere, we may take a vector-valued function on it and project its values to make vectors tangential to this sphere, thus obtaining a valid vector field. Since the projection operator preserves smoothness and since we can take a smooth vector-valued GP to begin with, it is clear that this approach may be used to build smooth vector fields.
We prove that (a) the resulting expression is, indeed, a kernel, and that (b) no expressivity is lost via the construction because all cross-covariance kernels arise this way.
Proposition 7. Let (X, g) be a Riemannian manifold, emb : X → Rd′ be an isometric embedding and F be a frame on X . We denote by P(·) : X → Rd×d ′ the associated projection matrix under F , and let f ′ : X → Rd′ be any vector-valued Gaussian process with matrix-valued kernel κ : X × X → Rd′×d′ . Then, the vector-valued function f = Pf ′ defines a Gaussian vector field f on X using the construction in Proposition 5, whose kernel under the frame F has matrix representation
KF (x, x ′) = Pxκ(x, x ′)PTx′ . (7)
3An embedding emb : X → Rd′ is called isometric if it preserves the metric tensor. By Nash’s Theorem [26], such embeddings exist for any d-dimensional manifold, with an embedded dimension d′ ≤ 2d+ 1.
4We emphasize again that f ′ is not a Gaussian vector field because it is not a random section. In particular, note that d′ > d for most embeddings.
Moreover, all cross-covariance kernels k : T ∗X×T ∗X → R arise this way. We call a kernel defined this way a projected kernel.
To construct these kernels we require scalar-valued kernels on manifolds to use as a basic building block. These are studied in the general Riemannian setting by Lindgren et al. [28] and Borovitskiy et al. [2]: relying on these kernels is the only reason we require the Riemannian structure. It is also possible to obtain such kernels using embeddings, following Lin et al. [27]. Similar techniques to those we consider are used by Freeden and Schreiner [14] to construct vector-valued zonal kernels on the sphere: in contrast, we work with arbitrary manifolds. The projection kernel idea is a very general way to build kernels for vector fields by combining scalar kernels, but effective scalar kernels, naturally, rely on Riemannian structure.5 Figure 3 shows random samples from Gaussian processes constructed with the described kernels.
The projected kernel construction both makes it easy to define cross-covariance kernels on general manifolds, and describes a straightforward way to implement them numerically by representing the embedded manifold in coordinates and calculating the resulting matrix-vector expressions. The constructed kernel depends on the embedding, but can be transformed appropriately if switching to a different embedding. Embeddings, in turn, are available for most manifolds of practical interest, and are obtained automatically for manifolds approximated numerically as meshes. Everything described is constructive and fully compatible with the modern automatic-differentiation-based machine learning toolkit, and most operations for constructing and/or sampling from specialized priors [44, 22, 23], including on spaces such as the sphere where specific analytic tools are available [6, 8, 10, 9]. With these kernels in hand, we thus proceed to study training methods.

heading: 3.2 Gauge Independent Variational Approximations
text: We now discuss variational inference for training GPs in the Riemannian vector field setting. Approximations, such as the inducing-point framework by Titsias [43] and Hensman et al. [17], approximate the posterior GP with another GP, termed the variational approximation. The latter is typically constructed by specifying a multivariate Gaussian at a set of test locations with a parameterized mean and kernel matrix. For example, Opper and Archambeau [31] consider N(m,S), where
m = K(·)z(Kzz + Σ) −1µ S = K(·,·) −K(·)z(Kzz + Σ)−1Kzz(Kzz + Σ)−1Kz(·). (8)
The variational parameters include a set of inducing locations z, a mean vector µ, and a blockdiagonal cross-covariance matrix Σ. Training proceeds by optimizing these parameters to minimize the Kullback–Leibler divergence of the variational distribution from the true posterior, typically using mini-batch stochastic gradient descent.
In the last decade, a wide and diverse range of inducing point approximations suited for many different settings have been proposed [43, 31, 25, 44, 49]. The vast majority of them employ coordinate-
5The structure of a smooth manifold is not rigid enough to define natural kernels. For instance, the smooth structure of the sphere is indistinguishable from the smooth structure of an ellipsoid or even of the dragon manifold from Borovitskiy et al. [2], but their Riemannian structures differ considerably.
dependent matrix-vector expressions. This raises the question, which of these constructions can be adapted to define valid variational approximations in the vector field setting?
To proceed, one can choose a frame and formulate a given variational approximation using matrices defined with respect to this frame. To ensure well-definedness, one must ensure that all these matrices, such as the kernel matrix and the variational parameter Σ in (8), are gauge independent. These considerations can be simplified adopting the pathwise view of GPs, and examining the random variables directly. For example, the variational approximation of Opper and Archambeau [31] shown previously in (8) can be reinterpreted pathwise as the GP
(f | y)(·) ≈ f(·) + K(·)z(Kzz + Σ)−1(µ− f(z)− ε) ε ∼ N(0,Σ) (9) where we view the matrices K(·)z,Kzz,Σ as linear operators between direct sums of tangent spaces: K(·)z : Tz1X ⊕ . . .⊕ TzmX → T(·)X and Kzz,Σ : Tz1X ⊕ . . .⊕ TzmX → Tz1X ⊕ . . .⊕ TzmX . By virtue of being defined at the level of vector fields using components that are all intrinsically valid, the posterior covariance of the resulting variational approximation is automatically gauge independent. Hence, checking gauge independence is then equivalent to deducing the domains and ranges of these operators from their coordinate representations, and checking if they are compatible. This applies to any variational family that can be constructed in the given manner.
The vast majority of inducing point constructions can be interpreted in this manner and thus extend readily to the Riemannian vector field setting by simply representing the necessary matrices in a chosen frame. In particular, the classical approach of Titsias [43] is gauge independent.

heading: 4 Illustrated Examples
text: Here, we showcase a number of examples that illustrate potential use cases of the models developed.

heading: 4.1 Dynamical Systems Modeling
text: Here, we show how Gaussian vector fields can be used to learn the equations of motion of a physical system—an important task in imitation learning, model-based reinforcement learning, and robotics. GPs are an attractive model class in this area owing to their ability to represent and propagate
uncertainty, which enables them to separate what is known about an environment from what is not, thereby driving data-efficient exploration.
For a prototype physical system, we consider an ideal pendulum, whose configuration space is the circle S1, representing the angle of the pendulum, with zero being at the bottom of the loop, and whose position-momentum state-space is the cylinder S1 × R. We consider conservative dynamics with additional friction applied at the pivot. Since this system is non-conservative, we cannot just learn the Hamiltonian of the system, but must learn the vector field over the state space that defines the dynamics of the system. The true dynamics of the system are given by the differential equations
H = p 2 2ml2 +mgl(1− cos(q)) dq dt = ∂H ∂p
dp dt = −∂H ∂q − b m p, (10)
whereH is the Hamiltonian of the system defining the conservative part of the dynamics, q and p are the position and momentum of the pendulum, m is the mass, l is the length, g is the gravitational field strength and b is a friction parameter. Experimental details can be found in Appendix B.
To learn this model, we initialise the system at two start points, and evolve the system using leapfrog integration. From these observations of position, we backward Euler integrate the momentum of the system, and from these position-momentum trajectories we estimate observations of the dynamics field. Using these observations, we condition a sparse GP. The result is an estimate of the system dynamics with suitable uncertainty estimates. In order to compute rollouts of these dynamics, we use pathwise sampling of this sparse GP [47, 48] for speed together with leapfrog integration.
Results can be seen in Figure 4. While the Euclidean GP performs reasonably well at the start of the rollouts, once the trajectory crosses the discontinuity caused by looping the angle back around to zero, the system starts to make incoherent predictions: this is due to the discontinuity arising from wrap-around condition of the angle. The manifold vector-valued GP does not have this issue as the learned and sampled dynamics fields are continuous throughout the state-space.

heading: 4.2 Weather Modeling
text: In this experiment, we show how vector-valued GPs on manifolds can be used in the context of meteorology, where geometric information often plays an important role in accurately modeling global weather fields [46, 38, 19]. Data assimilation in numerical weather forecasting refers to the practice of using observed data to update predictions of the atmosphere closer to the truth. Uncertainty plays a critical role here: it is not usually possible to observe the weather at all locations on the globe simultaneously, and taking into account observation uncertainty is crucial in numerical weather forecasting during the data assimilation step [24, 36]. In this section, we explore Gaussian processes as a tool for carrying out global interpolation of wind fields, while simultaneously performing uncertainty quantification, mirroring optimal interpolation techniques in data assimilation [20].
We consider a simplified setting, where the goal is to interpolate the wind velocity observed by the Aeolus satellite [37], which uses LiDAR sensors to measure wind velocity directly. To mimic this setting, we use an hour of the Aeolus satellite track during the period 2019/01/01 09:00-10:00 for the input locations and the wind velocity data (10m above ground) from the ERA5 atmospheric reanalysis data [18] interpolated at these locations, to simulate measurements taken from the Aeolus satellite. We subtract the weekly historical average wind velocity from the observations, before training the GP models, where the historical mean is computed from the hourly wind data (10 m above ground) from the WeatherBench dataset [35], available from 1979–2018. Further details can be found in Appendix B. We compare the results of a Matérn-3/2 manifold vector-valued GP regression model fitted on the wind anomaly observations along the Aeolus trajectory, with the results from a Euclidean Matérn-3/2 multi-output GP trained on the same data, except projected onto a latitude-longitude map.
Results are shown in Figure 5, where the benefits of using a manifold vector-valued GP become clear. When the satellite crosses the left/right boundary in the lat/lon projection, the outputs from the Euclidean vector-valued GP give rise to a spurious discontinuity in the uncertainty along the solid pink line. In addition, predictions become less certain in the Euclidean case as the satellite approaches the poles, which is simply an artifact of the distortion caused by projecting the spherical data onto the plane. By construction, the manifold vector-valued GP is able to avoid both of these issues, resulting in a more realistic prediction with much more uniform uncertainty along the satellite trajectory from pole to pole. In addition, the predictions from the manifold GP are more certain overall, due to the useful structural bias embedded in the kernel.

heading: 5 Conclusion
text: In this paper, we propose techniques that generalize Gaussian processes to model vector fields on Riemannian manifolds. This is done by first providing a well-defined notion of such processes on manifolds and then introducing an explicit method to construct them in a way that respects the underlying geometry. By virtue of satisfying the key condition of gauge independence, our construction is coordinate-free and thus meaningful on manifolds. In addition to this, we extend standard Gaussian process training methods, such as variational inference, to this setting, and verify that such methods are also compatible with gauge independence. This theoretical work gives practitioners additional tools for stochastic modeling of vector fields on manifolds. As such, its societal impact will be mainly determined by the applications that belong to the domain of future work. We demonstrate our techniques on a series of examples in modeling dynamical systems and weather science, and show that incorporating geometric structural bias into probabilistic modeling is beneficial in these settings to obtain coherent predictions and uncertainties.

heading: 6 Acknowledgements and Funding Disclosure
text: Michael is supported by the EPSRC Centre for Doctoral Training in Modern Statistics and Statistical Machine Learning (EP/S023151/1). Section 3 and Appendix A of this paper were solely financially supported by the RSF grant No21-11-00047 (https://rscf.ru/en/project/21-11-00047/), but was contributed to by all authors in equal proportion with the remaining work. Other sections were supported by the remaining funding. We also acknowledge support from Huawei Research. We thank Peter Dueben for useful suggestions in Section 4.2, including, in particular, making us aware of the Aeolus satellite. We are grateful to Nick Sharp for introducing us to the process used for making the three-dimensional figures appearing in this work.

Review summaries:
The paper proposes a natural extension of vector-valued GPs to data residing on Riemannian manifolds. The construction is straight-forward (replace matrix products with bilinear forms), which I consider a good thing. This also imply that it is relatively easy to apply the developed theory in practice. Toy demonstrations are provided to show that this is indeed the case. The paper is exceptionally well-written.

The paper suggests a construction of Covariance functions of multi-output Gaussian processes for vector fields on manifolds. The construction is independent of coordinate choices, but does depend on a Nash embedding of the manifold. In addition, the paper describes a clear mathematical framework of coordinate independent constructions of Gaussian processes. The approach is extended to variational approximations of such Gaussian processes, where the independence of the approach w.r.t. coordinate choices is again taken into consideration.

The authors propose a class of vector Gaussian processes (GP) on a class of Riemannian manifolds.  An important feature of the proposed GP is its ability to model vector fields on the Riemannian manifold. This is a broad application in physical sciences. One of the main novelties of the work is the development of gauge-equivariant kernels operators. Variation inference algorithms are extended for posterior inference.  Numerical examples demonstrate the practical performance of the proposed GP.

This paper proposed a way to define vector-valued Gaussian processes on Riemannian manifolds. A computational framework is provided to facilitate statistical inference using these models. The main contribution is an explicit methodology for characterizing vector-valued Gaussian processes through extrinsic kernels satisfying a gauge constraint, thus bypassing abstract mechanisms often employed to define stochastic processes on Riemannian manifolds. Two examples are provided to illustrate the power of the proposed models.