Wrapped Gaussian Process Functional Regression Model for Batch Data on Riemannian Manifolds
Jinzhao Liu Chao Liu Jian Qing Shi Tom Nye




DasAtom: A Divide-and-Shuttle Atom Approach to Quantum Circuit Transformation
    Yunqi Huang, Dingchao Gao, Shenggang Ying, and Sanjiang Li^*
        
Yunqi Huang and Sanjiang Li are with Centre for Quantum Software and Information (QSI), Faculty of Engineering and Information Technology, University of Technology Sydney, NSW 2007, Australia. Dingchao Gao and Shenggang Ying are with Institute of Software, Chinese Academy of Sciences, and University of Chinese Academy of Sciences, Beijing, China
Corresponding author (E-mail: Sanjiang.Li@uts.edu.au)

    September 5, 2024
============================================================================================================================================================================================================================================================================================================================================================================================================================================================================================



§ ABSTRACT 
Regression is an essential and fundamental methodology in statistical analysis. The majority of the literature focuses on linear and nonlinear regression in the context of the Euclidean space. However, regression models in non-Euclidean spaces deserve more attention due to collection of increasing volumes of manifold-valued data. In this context, this paper proposes a concurrent functional regression model for batch data on s by estimating both mean structure and covariance structure simultaneously. The response variable is assumed to follow a wrapped Gaussian process distribution. Nonlinear relationships between manifold-valued response variables and multiple Euclidean covariates can be captured by this model in which the covariates can be functional and/or scalar. The performance of our model has been tested on both simulated data and real data, showing it is an effective and efficient tool in conducting functional data regression on Riemannian manifolds.

 


§ INTRODUCTION

Regression models are ubiquitous and powerful tools in data analysis to reveal the relationship between independent variables and dependent variables. Most well known regression models are formulated under the assumption that the variables all lie in Euclidean vector spaces. However, statistical models of manifold-valued data are gaining popularity in various fields of scientific analysis, such as computational social science, medical imaging analysis and computer vision <cit.>. Together with the significant increase in manifold-valued data observed by current instrumentation, generalizations of Euclidean regression models to manifold-valued settings are needed. If the manifold structure is not accounted for then the reliability of model predictions in the target spaces is compromised. Consider, for instance, Figure <ref>, where a regression model is trained on manifold-valued data (represented by the black curves on a S2). The predictions of this model (illustrated by the blue curve) deviate from the manifold. However, when the manifold structure is incorporated into the model, the predictions (depicted by the red curve) align with the target manifold. This underscores the necessity of employing non-Euclidean structure for modelling manifold-valued data. An additional example demonstrating this principle in Kendall's shape space is provided in the Appendix.



In recent years, functional data have garnered increasing attention due to the continuous or intermittent recording of data at discrete time points. These data often exhibit non-linearity, and a common assumption is that they lie on a nonlinear manifold. For instance, image data, which can be influenced by random domain shifts, are known to reside in such manifolds . However, analyzing functional data poses challenges, particularly when dealing with spaces lacking global or local linear structures. Traditional methods like functional principal component analysis become infeasible in these cases. One specific difficulty arises from defining a suitable metric for data residing on a Riemannian manifold. Researchers have explored alternative approaches to address this issue. For example, <cit.> propose using point-wise Fréchet means for non-Euclidean time-varying random objects in general metric spaces. Additionally, they introduce point-wise distance trajectories between individual time courses and the estimated Fréchet mean trajectory. This innovative approach allows for a representation of time-varying random objects using functional data. In summary, the study of functional data on nonlinear manifolds presents exciting opportunities for advancing statistical methodologies and understanding complex data structures.

According to the manifold hypothesis, functional data can be mapped onto a low-dimensional nonlinear manifold, which is also known as manifold learning <cit.>. 
Gaussian processes (GPs) have proven to be powerful tools for modeling complex data structures. While GPs are commonly applied in Euclidean spaces, recent research has extended their applicability to Riemannian manifolds. <cit.> utilize the spectral theory of the Laplace-Beltrami operator to compute Matérn kernels on Riemannian manifolds. These kernels, widely used in physical sciences, capture both smoothness and spatial correlations. The spectral approach allows for efficient computation of Matérn covariance functions, enabling GPs to model data residing on curved surfaces. <cit.> introduce a vector-valued GP equipped with a matrix-valued kernel. This novel framework enables modeling of vector fields on Riemannian manifolds. Applications include geodesic flow modeling, fluid dynamics, and other scenarios where vector-valued data arise. Recently, <cit.> and <cit.> extend stationary GPs to compact Lie groups and non-compact symmetric spaces. These spaces play a critical role in spatiotemporal modeling. Those models define a metric based on geodesic distances between predictors, accommodating both functional and scalar predictors. In addition, <cit.> propose a regression model for functional data on manifolds. While existing approaches assume Euclidean response variables, certain problems involve responses residing directly on Riemannian manifolds. For instance, predicting flight trajectories based on flight time. In contrast to previous models, we propose a regression framework that maps Euclidean-valued predictors to manifold-valued functional responses.

A number of statistical approaches to analysis of manifold-valued data exist in the literature, especially in the context of principal component analysis (PCA). For example, <cit.> introduced principal geodesic analysis which is a generalisation of PCA on s. It replaces the first principal component with a principal geodesic constrained to pass through the intrinsic mean, with calculations performed approximately in Euclidean tangent space. <cit.> proposed a PCA technique for functional data on 2-dimensional Riemannian manifolds in which the authors adopt a regularization method by a smoothing penalty coherent with geodesic distances. Subsequently, <cit.> formulated Riemannian functional principal component analysis by mapping data to tangent spaces via the Riemannian logarithm map and then performing multivariate functional principal component analysis within the tangent space. 

In addition to generalizations of PCA to manifold data, generalizations of linear regression have also been studied. Similar to manifold-based PCA, the regression line is typically replaced with a geodesic within the manifold, which gives rise to more challenging optimisation problems in order to fit the models than in the standard Euclidean setting. Geodesic regression <cit.> is a generalization of linear regression to a manifold-valued setting in which an univariate independent variable is considered in ℝ and dependent variables lie in a . The author derived a gradient descent algorithm for model fitting via derivatives of the exponential map and Jacobi fields, because the least square method has no analytical solution under this setting. <cit.> then extended geodesic regression to the multivariate setting in which the independent variables are in ℝ^n and the dependent variable is still manifold-valued. In addition, the authors proposed a variational gradient descent method based on parallel transport which is more convenient for high-dimensional independent variables. In the context of shape analysis, <cit.> provided an analytic approach to geodesic regression in Kendall's shape space. Multiple linear regression has also been generalised by <cit.> for complex random objects in metric spaces with independent variables in ℝ^n. Using a least-squares approach, the authors derived asymptotic rates of convergence.

However, in many applications, non-Euclidean data cannot always be assumed to follow form of linear relationship. <cit.> introduced a kernel-based nonlinear regression model for both manifold-valued independent variables and manifold-valued dependent variables. <cit.> extended polynomial regression to s by introducing a class of curves on s which generalize geodesics and are analogs of polynomial curves in Euclidean space. <cit.> developed a regression model for a manifold-valued response variable in Riemannian symmetric spaces and covariates in Euclidean spaces with applications in medical imaging analysis. Moreover, <cit.> introduced an additive linear model for manifold-valued data which is based on the exponential map. In particular, they transform the manifold-valued data onto tangent spaces and then estimate the parameters of the additive linear model by a generalised least squares method.

Gaussian process regression (GPR) is a powerful non-linear and non-parametric model widely used for learning probability distributions over unknown functions. Operating within a Bayesian framework, GPR assumes that both the prior and likelihood follow normal distributions, allowing for the derivation of a posterior distribution that is also Gaussian. 
This ash been extended to solve the problems involving non-Gaussian data <cit.> and the use of other prior processes <cit.>. In addition, researchers have extended GPR to accommodate multi-output scenarios. Notable contributions include the work of <cit.>, <cit.>, and <cit.>, who have explored multi-output Gaussian process regression. Over the past decade, substantial progress has been made in developing Gaussian process regression models specifically tailored for manifolds, thereby expanding the model’s scope and applicability. One intriguing approach is the wrapped Gaussian process regression proposed by <cit.>. In this framework, the Riemannian manifold is linearized via the logrithm map, projecting data points onto a tangent space. This results in a non-parametric regression model with a probabilistic foundation. Importantly, the computational cost of wrapped Gaussian process regression remains relatively low, as it involves additional calculations only for the exponential map and logarithm map of each manifold-valued data point. The key innovation lies in defining a Gaussian distribution directly on the Riemannian manifold, leveraging insights from real-valued multivariate Gaussian distributions and the exponential map. This novel approach opens up avenues for Gaussian process regression on Riemannian manifolds. Within the Bayesian framework, researchers have derived manifold-valued posterior distributions, albeit under certain assumptions (e.g., the requirement of infinite injectivity radius). 

To the best of our knowledge, there is little current literature about regression models for functional batch data on s within a probabilistic framework. In light of this, we provide a nonlinear non-parametric regression model with uncertainty for batch data on a smooth . Called the wrapped Gaussian process functional regression model (WGPFR), it models the relationship of a functional manifold-valued response and mixed (scalar and functional) Euclidean-valued predictors by representing the mean structure and the covariance structure as different model terms. Specifically, a functional regression model is used for the mean structure with a functional manifold-valued response and scalar covariates, while a wrapped Gaussian process has been used for the covariance structure with functional covariates based on <cit.> who extend Gaussian process regression on s. In this way, the proposed WGPFR model extends the research field of non-parametric regression on manifolds.

The rest of this paper is organized as follows. In section <ref>, some basic concepts of functional data, s, Gaussian process regression and wrapped Gaussian process regression are reviewed. In section <ref>, we propose our model together with a method for inference via an efficient iterative algorithm.  Numerical experiments and real data analysis are reported in Section <ref>. Finally, we draw conclusions and discuss the results in Section <ref>.





§ BACKGROUND





 §.§ Functional Data


Functional data analysis is a field of study wherein each observation is a random function, manifesting as a curve, surface, or other variations over a continuum. This analysis is particularly pertinent to data recorded over identical intervals, with consistent frequency and numerous repeated observations. While data are frequently modelled by parametric models incorporating randomness in the realms of statistics or machine learning, functional data analysis contemplates a smooth process observed at discrete time points, offering greater flexibility and reduced parametric constraints, see e.g. <cit.> and <cit.>.

Furthermore, functional data can be procured from a diverse range of domains. For instance, in finance, volatility surfaces serve as a source of functional data. In the field of biomedical research, time-varying physiological signals yield functional data. Functional Magnetic Resonance Imaging (FMRI) is yet another domain that generates functional data. These examples underscore the broad applicability and interdisciplinary nature of functional data analysis.



 §.§ Preliminaries for Riemannian Manifolds


In this section, we review a few basic concepts of s, wrapped Gaussian distributions, wrapped Gaussian processes on manifolds and then set up basic notations. More detail on basic Riemannian geometry can be found in standard texts, such as that by <cit.>. 



  §.§.§ Concepts and Notation
 


Given a d-dimensional smooth differentiable manifold  and a tangent space T_p, for p∈, a Riemannian metric g_p:T_p× T_p→ℝ on  is a family of positive definite inner products which vary smoothly with p∈. Equipped with this Riemannian metric, we call the pair (, g) a Riemannian manifold.

The tangent bundle  of   is defined as a set  𝒯ℳ=∪_p∈ (p× T_p), where p is a point in  and v∈ T_p is a tangent vector at p. If γ is a smooth curve in ℳ, then its length is defined as the integral of ∂γ/∂ t where the norm is computed via the Riemannian metric at γ(t). For any pair (p,v)∈𝒯ℳ with v sufficiently small, there is a unique curve γ:[0,1]→ℳ of minimal length between two points  p∈ and q∈, with initial conditions  γ(0)=p, γ'(0)=v and γ(1)=q. Such a curve is called a geodesic. The  exponential map (p,v):× T_p→ at p∈ is defined as (p, v)=γ(1) for sufficiently small v∈ T_p. 
Consider the ball of largest radius around the origin in T_p on which  (p, v) is defined and let V(p)⊂ denote the image of (p,·) on this ball. 
Then the exponential map has an inverse on V(p), called the  Riemannian logarithm map at p, (p,·): V(p) → T_p. For any point q∈ V(p), the geodesic distance of p∈ and q∈ is given by d_ℳ(p,q)=(p,q). 

Figure <ref> shows the concepts above. For example, p and q are two points on a  ℳ, the tangent space at p is denoted as T_p(ℳ) and γ(t) refers to a geodesic between p and q. The tangent vector v is from p to q at T_p(ℳ) and p moves towarding v with distance v is q, which represents (p, q) and (p, v) respectively.



The exponential map and its inverse are analogs of vector space operations in the following way. Suppose p_ and q_ refer to two different points on  linked by a geodesic γ, and suppose the tangent vector from the former to the latter is denoted by pq_=(dγ/dt)_p. In addition, suppose p_E and q_E refer to two Euclidean data points and the vector from the former to the latter is given by pq_E. Table <ref> shows the relationship of operations between s and Euclidean spaces.
 


If (, d_) is a complete metric space that every Cauchy sequence of points in  has a limit that is also in , then the Hopf-Rinow theorem shows that every pair of points p,q∈ℳ is joined by at least one geodesic, a condition known as geodesic completeness, and, equivalently, (p,v) is defined for all (p,v)∈𝒯ℳ. However, given two points on ℳ, the geodesic between the points may not be unique, even if ℳ is geodesically complete. 
The injectivity radius at p∈ℳ, denoted p, is defined to be the radius of the largest ball at the origin of T_pℳ on which (p,·) is a diffeomorphism. It follows that the geodesic ball at p of radius p is contained within V(p) under which the Wrapped Gaussian process regression can be derived to a analytical form. In the following, we assume the data are on a Riemannian manifold with infinite injectivity radius, such as spheres and Kendall's shape spaces.

Some examples of s are shown in Appendix, including explicit formulae for the exponential and logarithm maps in certain cases.  



 §.§ Functional Data on Riemannian Manifolds


To study functional data located on s, <cit.> considered a d-dimensional complete Riemannian submanifold ^# of a Euclidean space ^d_0, with a geodesic distance d_ and a probability space (Ω, , P) with sample space Ω, σ-algebra  and probability P. The sample space of all -valued continuous functions on a compact interval ⊂ is denoted by ={ y: → | y ∈𝒞()}. Thus, we define -valued random functions to be functions y(t, ω), y: ×Ω→, satisfying y(·, ω) ∈. Let ={v: → ℝ^d |  ∫_v(t)^Tv(t)dt < ∞} refers to an ambient L^2 Hilbert space of d-dimensional square integrable functions on 𝒯. The inner product is defined as ⟨ v,u⟩=∫_v(t)^Tu(t)dt and the norm is given by v=⟨ v,v⟩^ for v,u∈.





 §.§ Gaussian Process Regression and Wrapped Gaussian Process Regression




  §.§.§ Euclidean Gaussian Process Regression


As described in <cit.>, a general Gaussian process regression model is defined as 

    y=f(x)+ϵ

 where  x∈ℝ^Q, y ∈ℝ, ϵ∼(0,σ^2).
  The regression function f is assumed to follow a Gaussian process prior, that is 

    f(·)∼ GP(μ(·),k(·,·;θ))

where μ(·)=E[f(·)] is the mean function of Gaussian process, and k(·,·;θ) is the covariance function with hyper-parameter θ.

Given training data = {x_i, y_i,i=1,…,n}, we have a multivariate Gaussian distribution that is

    (f(x_1),...,f(x_n))∼(μ(X),K(X,X))

where X refers to (x_1,...,x_n), μ(X) =  (μ(x_1),...,μ(x_n)) refers to a mean vector, and K(X,X) refers to a n× n covariance matrix in which the entry of the (i,j)-th element is k(x_i,x_j;θ).

To obtain the predictive distribution at a new input x^*, it is convenient to take advantage of the properties of Gaussian distribution:

    [      y; f(x^*) ]∼[ [   μ(X); μ(x^*) ]
    		,
    		[ K(X,X)+σ^2I  K(x^*,X)^T;    K(x^*,X)  K(x^*,x^*) ] ]

where K(x^*,X) is a 1× n covariance matrix, K(x^*,x^*) is scalar and their entries are similar to K(X,X). Fortunately, the predictive distribution has an analytical form which is

    p(f(x^*)|𝒟,x^*)∼(   K(x^*,X)^T(K(X,X)+σ^2I)^-1(y-μ(X))+μ(x^*),
       K(x^*,x^*)-K(x^*,X)^T(K(X,X)+σ^2I)^-1K(x^*,X))

For more details of Gaussian process regression, see <cit.>.

The kernel function in the covariance function can be chosen from a parametric family, such as the radial basis function (RBF) kernel:

    K(x,x';θ)=λ_1^2exp(-x-x'^2/2λ_2)

with the hyper-parameter θ={λ_1,λ_2} chosen by maximizing the marginal likelihood. Specifically, the hyper-parameters can be estimated by conjugate gradient descent algorithm <cit.> or Markov Chain Monte Carlo <cit.>.




  §.§.§ Wrapped Gaussian Distribution


In this subsection, we first review the definitions of the wrapped Gaussian distribution and wrapped Gaussian process on a Riemannian manifold ℳ, and then review wrapped Gaussian process regression (WGPR)  which is introduced by <cit.>.

Suppose  denotes a d-dimensional  and X is a random variable on . For μ∈ and a symmetric positive definite matrix K∈ℝ^d× d, we can define a wrapped Gaussian distribution as follows

    X =(μ,v), v∼(0,K)

where μ is called basepoint and v is a vector belongs to the tangent space at  μ which follows a multi-variate Gaussian distribution with zero mean and covariance matrix K (see Section 3.1 in <cit.>). 
The wrapped Gaussian distribution is formally denoted as

    X∼_(μ,K).


<cit.> describe a manifold-valued Gaussian process constructed by the wrapped Gaussian distribution, and the following material is based on their work. A collection (f(X_1),...,f(X_n)) of random points on a manifold  index over a set Ω is a wrapped Gaussian process, if every subcollection follows a jointly wrapped Gaussian distribution on .
It is denoted as 

    f(·)∼𝒢𝒫_(μ(·),k(·,·;θ))

where  μ(·)=E[f(·)] is the mean function on a Riemannian manifold, and k(·,·;θ) is a covariance function with hyper-parameter θ on a tangent space at  μ(·).

Consider the wrapped Gaussian process regression model which relates the manifold-valued response p∈ to the vector-valued predictor x∈ℝ^Q though a link function f:ℝ^Q→ which we assume to be a wrapped Gaussian process.  
Analogously to Gaussian process regression, the authors also derive inferences for wrapped Gaussian process regression in Bayesian framework. Specifically,
given the training data 𝒟_={(x_i,p_i)|x_i∈ℝ^Q,p_i∈, i=1,⋯,n}, the prior for f(·) is assumed to be a wrapped Gaussian process, which is f(·)∼𝒢𝒫_(μ(·),K(·,·;θ)). Then the joint distribution between the training outputs p=(p_1,⋯,p_n) and the test outputs p^* at x^* is given by

    ([ p^*;   p ]) ∼𝒩_(([ μ^*;   μ ]),([ K^* *   K^*;  K^*⊤     K ]))

where μ^*=μ(x^*), μ = (μ(x_1),...,μ(x_n)), K=K(x,x),
K^*=K(x,x^*),     K^* *=K(x^*, x^*) with x=(x_1,⋯,x_n).  

Therefore, by Theorem 1 in <cit.>, the predictive distribution for new data given the training data can be obtained by

    p^*|p∼Exp(μ^*,v^*),

where

    v^*∼𝒩_ℳ(m^*, Σ^*), m^*=K^* K^-1(μ,p), Σ^*=K^* *-K^* K^-1 K^*⊤.


Similar to the comparison of operators between Riemannian manifolds and Euclidean spaces, Table <ref> shows the analogs of Gaussian processes between these two spaces.







§ MAIN MODEL AND INFERENCE



Gaussian process regression is a powerful non-parametric tool in statistics and machine learning. In Section <ref>, we adapt GPR to functional data on manifolds  by proposing the WGPFR model for batch data on . Inference of the WGPFR model is presented in Section <ref>, where the statistical properties of the model parameters are studied and an efficient algorithm for inference is proposed. 



 §.§ Wrapped Gaussian Process Functional Regression Model


The data consist of discrete observations of M random functional curves on , so that the m-th curve is denoted as y_m(t), t ∈, m=1, …, M. We assume that all curves are observed at a set of times t_i, i=1,…,N. Throughout, i will index time and m will index curves. 
The set of observations is therefore y_mi=y_m(t_i), i=1,…,N, m=1,…,M. Associated with the m-th curve, we assume there is a real vector-valued functional covariate xm(t)∈ℝ^Q. We assume, as for the functional responses y_m(t), that  xm(t) is also observed at the times t_i, so that we are given points xm=xm(t_i). Finally, we assume there is a vector-valued covariate um∈ℝ^P observed for each curve y_m. Thus, the data for the m-th curve comprises y_mi, xm, t_i, i=1,…, N together with the vector um. These data, for fixed m, are called a batch. 

Next we define our hierarchical model for the data. It has the following overall structure, and details of each part are filled out in what follows. At the top level of the hierarchy, there is a mean curve μ_*(t)∈. There are two stages to obtain each curve y_m from μ_*. First there is a perturbation of μ_* to obtain an unobserved curve μ_m. We assume this perturbation depends on the vector covariate um and a functional model parameter β(t)∈ℝ^P, but not on the functional covariates (see Equation (<ref>)). Secondly, we assume y_m is obtained from μ_m via a vector field τ_m along μ_m which depends on the functional covariate xm(t) (see Equation (<ref>)). Figure <ref> shows a schematic of the overall model. 





  §.§.§ Definition of Mean Structure


It is difficult to infer non-linear relationships non-parametrically between a functional response variable and multi-dimensional functional covariates even when both of them are observed in a space for real-valued functions <cit.>. Therefore, we assume that the mean structure μ_m(·) depends on the scalar covariates u_m only. Specifically we assume 

    μ_m(t) = (μ_*(t), u_m^Tβ(t))

where μ_*(t)∈ and β(t)∈ T_μ_*(t) play  a similar role to the intercept term and slope term of the function-on-scalar linear regression model in Euclidean space. 
These are model parameters which we estimate. Furthermore, equation (<ref>) can be seen as a functional version of geodesic regression proposed in <cit.>, in which a geodesic regression model for manifold-valued response and scalar-valued predictor is proposed in the form y=(p,xv) with (p,v)∈𝒯ℳ where x∈ℝ is the predictor variable.




  §.§.§ Definition of Covariance Structure


Conditional on μ_m, the m-th curve y_m is modelled as

    y_m(t)=(μ_m(t),τ_m(t)), for m=1,...,M

where τ_m(t) ∈ T_μ_m(t) is a vector field along μ_m. When the data are observed with random noise, <cit.> defines the geodesic model as y=((p,xv),ϵ). Similarly, the proposed WGPFR model with error term can be defined as  

    y_m(t)=(μ_m(t),τ_m(t)+ϵ_m(t)),    for m=1,⋯,M;

where ϵ_m(t) are independent measurement errors following a multivariate normal distribution with mean zero and covariance matrix σ^2𝐈. For functional data, in addition to the measurement error term, the most important part we are interested is the underlying covariance structure contained in term τ_m(t), which is Equation (<ref>). In what follows, we explain how τ_m(t) is modelled as depending on the functional covariate xm(t).  

We assume that the vector field τ_m(t)∈ T_μ_m(t) can be represented in local coordinates as a d-dimensional vector function τ_m(t)=(τ_m,1(t),...,τ_m,d(t))∈^d. 
We further assume this vector-valued function follows a Gaussian process with mean zero and an underlying covariance structure k_m(·,·;θ), which is

    τ_m(t) ∼ GP(0, k_m(·,·;θ)).

Specifically, given two observation times t_i and t_mj from the m-th batch of data, we assume that the covariance of the two random vectors τ_m(t_i) and τ_m(t_mj) a function of the observed covariate values xm(t_i) and xm(t_mj).

Under the framework of the wrapped Gaussian process in <cit.>, we can represent the above Gaussian process as the wrapped Gaussian process for manifold functional data. That is, 

    y_m(t)∼𝒢𝒫_(μ_m(t),k_m(·,·;θ))

where μ_m(t)∈ and k_m(·,·;θ)∈ T_μ_m(t) are defined as above. 

To summarize, the WGPFR model studies the non-linear relationship between the manifold-valued functional response and mixed (functional and vector) Euclidean-valued predictors in two parts, as shown in Figure <ref>. The first part is the mean structure which is defined through a functional regression to model the manifold-valued functional mean with the vector-valued predictor. The second part is called the covariance structure, and for which the inference is based on a Bayesian method with the prior assumed to be Gaussian process.

Compared to the geodesic regression in <cit.>, the WGPFR additionally incorporates the information of dependence between data modelled by covariance structure which is the main focus of interest for functional data analysis. Compared to the WGPR model in  <cit.>, this paper further models the mean part which relates the functional mean on the manifold to the scalar predictor on Euclidean space, while the mean part in <cit.> is not related to scalar predictors. Note that the scalar predictor plays an important role in the prediction of the model, especially for the batch data we consider in this paper. The proposed WGPFR model is more complicated but much flexible than existing models and induces a more complex inference procedure which will be illustrated in detail in the next subsection.




 §.§ Inference for WGPFR





  §.§.§ Estimation of Mean Structure



Since the mean structure model (<ref>) can be seen as a functional version of the geodesic regression model <cit.>, an intuitive inference method consists of optimizing μ_*(t) and β(t) simultaneously for each fixed t. However, this would ignore information about the correlation of the function μ_* and β(t) at different times. On the other hand, if the number of observed times t_i is large, and we account fully for correlations between all different times, then the computational cost will be very high, unless we use some efficient methods, such as <cit.>. Due to the dependent structure, we cannot conduct optimization at each fixed t separately, otherwise the correlation information across different t will be lost. It can also be seen that inference of the functional mean structure model is indeed challenging.  

In this paper we borrow an idea of <cit.> to estimate μ_*(t) and β(t).  Note that μ_*(t) plays the role of an intercept in the linear regression in Euclidean space, which is usually defined as the mean of the response. Similarly, for functional data y_m(t), m=1,...,M on a , we can define the intrinsic population Fréchet mean function under the assumption of existence and uniqueness:

    μ_0(t)=p(t)∈arg min  E[d_(y(t),p(t))^2]

where d_(·,·) denotes geodesic distance. When t is fixed, μ_0(t) refers to a point on a ; when t is variable, μ_0(t) refers to a curve on a . Additionally, since each y_m(t) ∈ is continuous, μ_0(t) is also continuous (see equation (1) in <cit.>).

The intrinsic Fréchet mean function μ_0(t) can be a good choice for μ_*(t) since it is a natural mean of manifold-valued functional data and it is reasonable to assume y_m(t), m=1,...,M are confined to stay within the radius of injectivity at  μ_0(t) for all t∈𝒯. Then the inverse exponential map can be  well-defined. We estimate μ_*(t) at each time t via the sample Fréchet mean of the points y_m(t). This requires that the times t_i are the same across all curves m=1,…,M. Standard gradient descent methods <cit.> can be used to obtain the Fréchet sample mean at each time t=t_i to obtain estimates μ̂_0(t_i) of μ_*(t_i), i=1,…,N.  
Denote the logarithm process data as

    V_m(t)=(μ_0(t),y_m(t)),
 
where V_m(t) is a d-dimensional tangent vector of y_m(t) at μ_0(t) for a fixed t. 
Since we assume  is a d dimensional closed Riemannian submanifold of a Euclidean space,  it is reasonable to identify the tangent space T_μ_0(t)  as a subspace of the Euclidean space ℝ^d_0 with d_0≥ d. Then V_m(t) can be seen as ℝ^d_0 valued square integrable functions, which can be represented by a set of basis functions. Similarly, the slope function β(t) can also be represented by this set of basis functions. In this way, we can not only easily estimate the intercept term μ_0(t) but also the slope term β(t) can be obtained by least square methods on ℝ^d_0. This reduces the computation complexity compared with the optimization method in <cit.>. Next we illustrate the estimation procedure in detail.
 
For the estimation of the intrinsic Fréchet mean μ_0(t), the sample Fréchet mean is calculated by minimizing the following loss function for each t∈𝒯:

    μ̂_0(t) = p(t)∈arg min 1/M∑_m=1^Md_(y_m(t), p(t))^2.

It is natural to use the gradient descent algorithm introduced by <cit.> to estimate the sample Fréchet mean for each fixed t. Using <cit.>, the gradient is given by 

    ∇ d_(p,y_m(t))^2 = -2(y_m(t),p).

An implementation of the gradient descent algorithm to solve equation (<ref>) is given in Algorithm <ref>. In practice, we choose the step size l=1/N where N is the number of samples and ϵ=10^-8 which is relatively small.


 
For the estimation of  β(t)=(β_1(t),⋯,β_p(t)), the first step is to find a reasonable basis to represent β(t) in the tangent space. Under the framework of <cit.>, suppose we have an arbitrary system of K orthonormal basis functions

    Φ_K={ϕ_k∈ℍ | ϕ_k(t)∈ T_μ_0(t),⟨ϕ_k,ϕ_l ⟩=δ_kl,k,l=1,⋯,K},

where ℍ={ν:𝒯→ℝ^d,∫_𝒯ν(t)^Tν(t)dt<∞} is a L^2 Hilbert space of ℝ^d with inner product ⟨ν,u ⟩=∫_𝒯ν(t)^Tu(t)dt} and norm ‖ v‖=⟨ v,v ⟩^1/2, δ_kl=1 if k=l and 0 otherwise. Note that the values of each ϕ_k(t) at each time t∈𝒯 is restricted to the tangent space T_μ_0(t), which are identified with ℝ^d_0. Define the K-dimensional linear subspace of ℍ:

    ℳ_K(Φ_K):={ x(t)=∑_k=1^K a_kϕ_k(t)  for  t ∈𝒯| a_k∈ℝ}.

Then the slope functions  β_j(t)∈ T_μ_0(t), j=1,⋯,p  can be approximated by a linear span on ℳ_K(Φ_K) with expansion coefficients  b_jk: 

    β_j(t)≈∑_k=1^K b_jkϕ_k(t).

Then the estimation of β_j(t) can be transformed to the estimation of b_jk, k=1,2,⋯,K, j=1,2,⋯,p.

Defining V_m(t)=(μ_0(t),y_m(t)), V̂_m(t)=(μ̂_0(t),y_m(t)) and
W_m(t) =(μ_0(t),μ_m(t)), then the multiple linear function-on-scalar regression model becomes 

    V_m(t)=W_m(t)+e_m(t)≈∑_j=1^p∑_k=1^Ku_mjb_jkϕ_k(t)+e_m(t)

where e_m(t)=V_m(t)-W_m(t) which are assumed to be independent of the covariate u_m for m=1,⋯,M.

Since V_m(t) can not be obtained directly, we use V̂_m(t) as the response of the multiple linear regression model, so that the parameters of the regression model are estimated by minimizing the following loss function with observed data:

    L(b_jk)=∑_m=1^M∑_i=1^nV̂_m(t_i)-∑_j=1^p∑_k=1^Ku_mjb_jkϕ_k(t_i)^2.

The b_jk coefficients are calculated by standard least squares methods.

To be more specific, let
V̂=(V̂_im)_i=1,⋯,nd_0;m=1,⋯,M denote a nd_0× M matrix such that the m-th column denotes the m-th batch, and the i-th row denotes the vectored element (V̂_m1(t_1),⋯,V̂_m1(t_n),⋯,V̂_md_0(t_1),⋯,V̂_md_0(t_n)). Similarly, let 
Φ=(ϕ_ik)_i=1,⋯,nd_0;k=1,⋯,K denote a nd_0× K matrix; let U denote a M× p matrix for which each row u_m indicates a data point; and let B denote the p× K coefficient matrix with element b_jk for j=1,⋯,p, k=1,⋯,K.  Then the above loss function can be rewritten as

    L(B)=vec(V̂)-(Φ⊗ U)vec(B)^2.
 
The least squares estimate of B is obtained by

    vec(B̂)=((Φ⊗ U)^⊤(Φ⊗ U))^-1(Φ⊗ U)^⊤ vec(V̂)


Consistency of the estimator of the sample Fréchet mean was proved in <cit.>. The following theorem establishes consistency of the estimator of the regression coefficients; the proof is left to the appendix.


    
    Under the conditions (C1)-(C5) in the Appendix, the coefficient  vec(B̂) is a consistent estimator with probability tending to 1 as n →∞ in the sense that
    
    vec(B̂)-vec(B)=o_p(1).

 

It follows that the functional slope coefficients can be estimated by  β̂(t)=(β̂_1(t),⋯,β̂_p(t)) where

    β̂_j(t)= ∑_k=1^Kb̂_jkϕ_k(t)

for j=1,⋯,p.
The estimated mean structure for the m-th curve is given by

    μ̂_m(t)=(μ̂_0(t),u_mβ̂(t)).

The function μ̂_m(t) is expected to approximate μ_m(t). However we omit the rigorous proof of the consistency of μ̂_m(t) on the manifold due to some technique difficulties.



  §.§.§ Estimation of Covariance Structure



From the section above, we obtain an estimate of the mean structure μ_m(t) which is a continuous function on . In this section, we will focus on inference of the covariance structure τ_m(t)∈ T_μ_m(t). We assume that the tangent spaces T_μ_m(t) for t∈𝒯 can be identified with ^d_0 via some smooth local basis of orthonormal vector fields along μ_m(t). We have mentioned above that the covariance structure can be related to another functional covariate x_m(t), which is a ℝ^q-valued function. If we denote τ_m(t)=τ_m(x_m(t))=(τ_m1(t),⋯,τ_md_0(t)), then the correlation of different components in τ_m(t) could be estimated via a cross-covariance function model, such as the convolved Gaussian process <cit.>. However, if there are n observations τ_md(t_i), i=1,⋯,n, then the size of a covariance matrix in a Gaussian process of τ_md(t), d=1,⋯,d_0 is n× n while the size of a cross-covariance matrix in a convolved Gaussian process of τ_m(t) is n^d_0× n^d_0, which is computationally expensive. As a result, we will consider different dimensions independently. Specifically, we assume

    τ_md(t)∼ GP(0,K_md(x_m(t),x_m(t');θ_md)),

where K_md(·,·;θ_md) denotes a covariance kernel depending on hyper-parameter θ_md.

Given the value of τ_md(t_i), i=1,⋯,n, denote τ_md=(τ_md(t_m1),...,τ_md(t_mn)), then for any new input t_m^*, the conditional distribution of τ_md(t_m^*) is

    τ_md(t_m^*)|_md∼(μ_md^*, Σ_md^*), μ_md^* = k_md^*TK_md^-1τ_md, Σ_md^* = k_md^** - k_md^*TK_md^-1k_md^*

where k_md^*=(K_md(t_m1,t_m^*),...,K_md(t_mn,t_m^*)) and k_md^** = K_md(t_m^*,t_m^*).

In Gaussian process regression, the hyper-parameters for the d-th dimension, θ_md, can be estimated by maximizing the sum of marginal likelihood functions for each batch, i.e. maximizing

    (2π)^-n/2(K_md)^-e^-τ_md^TK_md^-1τ_md.


Under some regularity conditions, the estimator θ̂_md converges to θ_md <cit.>. With the estimated hyper-parameters, the estimated covariance structure of the d-th coordinate for the m-th batch is given by

    τ̂_md(t)|τ_md∼ GP(k_md^T(t)K_md^-1τ_md, K_md(t,t)-k_md(t)^TK_md^-1k_md(t)|θ̂_md))

where k_md(t)=(K_md(t,x_m(t_m1)),⋯,K_md(t,x_m(t_mn))). Under some mild conditions, the above estimator of the covariance structure is information consistent, as the following theorem states.


    
    Suppose (1) the underlying true vector field, τ_m(t) along μ_m(t) can be represented using a local basis by real-valued functions τ_md, d=1,…,d_0, and that each function has a Gaussian process prior with mean zero and bounded covariance function K_md(·,·,θ_md) for any covariate x_m(t), (2) K_md(·,·,θ_md) is continuous in θ_md, further assume the estimator θ̂_md converges to θ_md as n→∞ almost surely, and (3) the mean structure μ_m(t) is known.
    
    Then τ̂_md(·) is information consistent to the true τ_md(·), which means an estimator that guarantees convergence to the true parameter value as the sample size grows, if the reproducing kernel Hilbert space norm τ_md_k is bounded and the expected regret term  E_X(log|I+σ^-2C_n n|)=o(n),  where C_n n=(k(x_i, x_j ; θ)) is the covariance matrix over x_i, i=1, …, n, and σ^2 is the variance of the measurement error ϵ_m(t).	


The proof is given in Appendix. Note that the above theorem holds under the condition that μ_m(t) is known, but in practice, we only have an estimate of μ_m(t). Furthermore, the term τ_md in the posterior distribution can not be observed directly, and it is approximated by τ̃_md=(τ̃_md(t_m1),⋯,τ̃_md(t_mn)) which is the realization of the d-th element of  τ̃_m(t)=(τ̃_m1(t),⋯,τ̃_md_0(t))=(μ̂_m(t),y_m(t)).



  §.§.§ Update Mean Structure and Covariance Structure



After obtaining the estimated mean structure and covariance structure, we are able to make  predictions with given new inputs. In order to improve the performance of our model, we introduce an algorithm which can update the estimated mean structure and covariance structure iteratively.

The loss function of the m-th curve at time point t_i with the estimated mean structure and estimated covariance structure is given as

    E=∑_m=1^M∑_i=1^nd_(((μ̂_0(t_i),μ̂_m(t)), τ̂_m(t_i)), y_m(t_i))^2

where μ̂_m(t)=∑_j=1^p∑_k=1^Ku_mjb̂_jkϕ_k(t_i).

Given the above estimated covariance structure τ̂_m(t_i), the mean structure can be updated by a gradient descent algorithm where the gradient is

    ∇_μ_m(t_i)E = -d_μ_m(t_i)(μ_m(t_i),τ̂_m(t_i))^†(y_m(t_i),(μ_m(t_i),τ̂_m(t_i))),

where † is an adjoint with respect to the Riemannian inner product, which plays a similar role to parallel transport. In practice, a variational method for gradient descent <cit.> can be used as a substitute for gradient in Equation (<ref>). Thereafter, we update the mean structure from μ̂_m(t) to μ̂_m^(1)(t) (here the superscript ^(1) means the 1-st iteration). The updated coefficients b_jk^(1) can be estimated by minimizing the loss function

    L(b_jk)= (μ̂_0(t_i), μ̂^(1)_m(t_i))-∑_j=1^p	∑_k=1^Ku_m,jb_jkϕ_k(t)^2,

which is a linear least squares problem.

Given the updated mean structure μ̂_m^(1)(t), we can re-calculate the covariance structure τ̂_m^(1)(t). The only difference from Section <ref> is that the estimated mean structure is replaced by the updated mean structure. This two-way updating procedure is then repeated iteratively, stopping when the difference of the updated mean or covariance structure smaller than some given threshold. The algorithm of the updating procedure is specified in Algorithm <ref>.







§ NUMERICAL EXPERIMENTS



In this section we demonstrate the WGPFR model on two Riemannian manifolds: S^2 and Kendall's shape space. As previously, we suppose there are M curves on a Riemannian manifold; in what follows we will simulate data points on these random curves in different ways to form different scenarios of interpolation and extrapolation. For the extrapolation problem, we use our model to predict the last 15 data points (given all preceding data points) to test the long-term extrapolating performance, and the last 5 data points to test the short-term extrapolation performance. Moreover, since the Algorithm <ref> is calculated in Euclidean space, we use a optimizer from SciPy which is powerful and widely-used for minimization problems.



 §.§ Regression analysis on S^2



  §.§.§ Simulation scenario on S^2



Suppose observation times t are equally spaced in the interval [0, 1] with N points, where m∈{1,...,m_1,m_1+1,...,M} and M=m_1+m_2. To test the performance of our model with different number of observed curves and data points, we considered 30 and 60 batches and 20, 40 and 60 data points on each batch respectively, which is M ∈{30, 40} and n∈{20,40,60 }. The simulated data are shown in Figure <ref> and in the remainder of this subsection we explain how the data were generated.



As introduced previously, the WGPFR model consists of three parts: the intrinsic Fréchet mean function, the mean structure and the covariance structure. In order to simulate data, we first define the Fréchet mean function μ_0(t) on S^2 by

    μ_0(t) = (p, (sin(tπ/2)^2, sin(tπ)^3,0))

where p=(0,0,1) refers to a point on a S^2.

Suppose there are two batches, and that we generate two mean structures based on batch-specific covariates u_m and Fourier basis functions β(t), defined by

    μ_m(t) = (μ_0(t), u_s^Tβ(t)), m=1,...,m_1+m_2

where u_m=(1,0) for m=1,...,m_1 and u_m=(1,1) for m=m_1+1,...,m_1+m_2 refer to the covariates for two scalar batches respectively and β(t)=(β_1(t),β_2(t)) refers to tangent vector-valued basis functions. Specifically, we define

    β_1(t)    = R_t(∑_i=1^20 w_1,i(ϕ_i(t), ϕ_i(t+1/2), ϕ_i(t+2/2))), 
    β_2(t)    = R_t(∑_i=1^20 w_2,i(ϕ_i(t), ϕ_i(t+1/2), ϕ_i(t+2/2))), i=1,...,20

where

    w_1,i   = i/120,
    
        w_2,i   = -1/2√(sin(i/60)), i=1,...,20,
    ϕ_i(t)    = 1/√(2),  for i=0,
    ϕ_i(t)    = sin(2π it),  for i = 1,3,5,...,19,
    ϕ_i(t)    = cos(2π it), for i = 2,4,6,...,20,

and R_t is a rotation matrix transferring the vector ∑_i=1^20 w_1,i(ϕ_i(t), ϕ_i(t+1/2), ϕ_i(t+2/2))) or 
∑_i=1^20 w_2,i(ϕ_i(t), ϕ_i(t+1/2), ϕ_i(t+2/2))) to T_μ_0(t), this ensures that the vectors β_1(t),β_2(t) are tangents at μ_0(t). For example, given a point p and its tangent vector 𝐯, the rotation matrix is defined by

    R_t(𝐯)=𝐯-⟨ v,p⟩/‖𝐯-⟨ v,p⟩‖‖𝐯‖.

We use the following covariance function

    k(t_i, t_j) = v_0exp(-1/2w_0 (t_i-t_j)^2)+a_0+a_1t_it_j+σ^2δ_ij

which is a combination of a squared exponential kernel (often used for stationary Gaussian processes) and a linear kernel (often used for non-stationary Gaussian process), this can be used for most of Gaussian process regression <cit.>. Since S^2 can be embedded into ^3, for each dimension of tangent vector of covariance structure, the hyper-parameters (v_0,w_0,a_0,a_1,σ) are given as θ_1=(0.012, 3, 0.01, 0.01, 0.02), θ_2=(0.017, 3.1, 0.011, 0.012, 0.015) and θ_3=(0.015, 3.2, 0.012, 0.013, 0.01). Given these three Gaussian processes with zero mean and kernel <ref>, we can generate a tangent vector for the covariance structure. However, in practice, the tangent vector may not in the target tangent space because of numerical rounding. For example, if the ideal tangent vector is (1,1,0), our generated tangent vector might be (0.99,1,0.01), which will result in the exponential map and inverse exponential map incalculable (defined in Equation (<ref>) and Equation (<ref>) respectively). Therefore, the tangent vector must be projected into the correct tangent space, which is denoted as v_3(t)|θ_1,θ_2,θ_3. The subscript 3 in v_3 refers to 3-dimensional. As a consequence, the formula to generate manifold-valued data for numerical experiments is written as

    y_s,m(t)=(μ_s(t),v_3(t)|θ_1,θ_2,θ_3), s=1,2, m=1,...,M_s.




  §.§.§ Model assessment on S^2


In this section, we delete some generated data points on a randomly selected curve, say y_r(t), in different ways to form training data sets and then calculate the predictions for these deleted data points (i.e. the test data). The performance of WGPFR model is assessed by comparing the root mean square error between predictions and their real data by Euclidean distance, since a sphere can be embedded in ℝ^3 easily.

Specifically, we selected 15 data points on y_r(t) uniformly at random as our test data set, and all remaining data points were used as the training data set, which is a typical interpolation problem. This scenario is denoted as Type 1 prediction. As to Type 2 prediction, the short-term forecasting, the last 5 data points on y_r(t) are considered as test data set which is a typical extrapolation problem. Analogously to the short-term forecasting, we also test the performance of WGPFR for long-term forecasting, the Type 3 prediction. In Type 3 prediction, we choose the last 15 data points as test data. 

For comparison, in each scenario, the same training data set is used for several other models, such as functional linear regression on Riemannian manifolds (FLRM) and wrapped Gaussian process Fréchet mean regression (WGFmR). Specifically, FLRM is the mean structure (<ref>) without covariance structure; WGPFmR consists of mean structure and covariance structure in which the mean structure is the sample Fréchet mean point for all training data. In addition, the WGPFmR model does not have the updating part. We compare the performance of these three models not only to show a significant improvement of predictive accuracy from concurrent model (consider mean structure and covariance structure simultaneously), but also from the intervatively updating algorithm <ref>. Moreover, we can use some models, such as Gaussian process functional regression without manifold structure to fit the data and make predictions. However, such models are meaningless since the inferences cannot be guaranteed on the right space which is manifold. Thus, it might be not suitable as a baseline and we only consider comparison with manifold structure.

We replicate each simulation study 100 times. Thus, we test the performance of our model on thousands of test data points. The numerical results reported in Table <ref> are the average of root-mean-square-error (RMSE) in every single replication. Using the embedding S^2 ⊆^3, it is reasonable to use the Euclidean norm between points as a distance function (chordal metric), which provides a method to calculate the RMSE. 

In Type 1 prediction, Type 2 prediction and Type 3 prediction, WGPFR always provides the best prediction. The predictive accuracy of FLRM is better than that of WGPFmR in most scenarios, because FLRM only learns the mean structure of the training data which might be more useful for trend forecasting. In extrapolation, the long-term forecasting and short-term forecasting, since test data are distant from the training data, the output of GRP is usually inaccurate and then the mean structure mainly determines the accurate of prediction. Because we used the Fréchet mean point as the mean function of WGFmR, the prediction of WGFmR is very poor. 

Comparing long-term and short-term forecasting, since the test data of the former is closer to the training data than that of the latter, the predictions of short-term forecasting are more accurate than that of long-term forecasting. In addition, from the table we can see that when the numbers of curves are fixed, the RMSE between prediction and real data decreases with the increasing number of points; when the number of points is fixed, the RMSE also decreases with the increasing number of curves. 



Since Gaussian process regression is a Bayesian model, we also compare the log-predictive likelihood of the covariance structure which provides the randomness of this model. <cit.> introduce a calculation of log pointwise predictive density in practice. Thus we compute this index by summing the log pointwise predictive density in each dimension and the result is shown in Table <ref>. 

Moreover, discrete Fréchet distance is a metric to measure the similarity between two curves which considers the locations and orderings between predictions and real data. This method is widely used in GPS data, especially for comparing trajectories of animals, vehicles or other moving objects <cit.>. The results show that our model provides small predictive error under this metric. In other words, our model is effective even the measurement is different from its objective function. The numerical results are shown in Table <ref>.

We also show the RMSE of β in Table <ref>. We can see that, when the number of curves is fixed, the RMSE between estimated β and real β decreases with the increasing number of data points; however, when the number of data points is fixed, the RMSE between estimated β and real β is almost the same (0.1779 compared to 0.1773 and 0.0212 compared to 0.0214).






 §.§ Regression Analysis on Kendall's Shape Space


We firstly generated μ_0(t), for t∈ (0,1), which describes a varying shape. Specifically, the shape is a circle at the beginning and becomes a square in the end. The varying shapes are shown in Figure <ref> and they are generated based on an elliptic equation after removing scale, rotation and translation.



We generated data corresponding to 80 landmarks, given by an element of ℂ^2× 80(See the Appendix for more details about Kendall's shape space)

and the mean shape function μ_0(t) is defined as

    μ_0,z(t) =
        
        (z/20-1,√((z/20-1)^2+(1-(z/20-1)^2)(1-t/2)^2)i),    for  z=1,...,40 
    
        (1-z/20,-√((1-z/20)^2+(1-(z/20-1)^2)(1-t/2)^2)i),    for  z=41,...,80


Then, we use the scalar covariates to generate tangent vectors from the mean function μ_0(t) to mean structures μ_1(t) and μ_2(t) by

    μ_1(t)    =(μ_0(t), ∑_p=1^Pu_1,psin(t)^3sin(r(t))), 
    μ_2(t)    =(μ_0(t), ∑_p=1^Pu_2,psin(t)^3sin(r(t))cos(r(t))),

where P refers to the number of elements in 𝐮_1, u_1,1=0, u_1,2=0, u_2,1=1, u_2,2=2 and r(t) refers to a ranking function of t. For example, when t=t_3, r(t)=3. Thus, we obtain 3 curves in a batch and each curve is consist of 10 points (curves). 

The generated μ_0(t) and the mean structures μ_1(t) and μ_2(t) are shown in Figures <ref> and <ref> respectively.



We estimate the mean structure and covariance structure as discussed in Section <ref> and Section <ref>. The training data consist of all data in batch one together with the first and the second curve in batch two; the difference between these types of prediction is the handling of the third curve in batch two. Specifically, for Type 1 prediction, 8 time points are randomly selected to form the training data and the remaining 2 points are used as the test data; with regard to Type 2 prediction, we add the first half into the training data and the remaining for the test data; as for the Type 3 prediction, the last 3 time points are considered as the test data and the remaining are added into the training data. Each numerical experiment is repeated for 100 times. In addition, to compare the performances among different regression models, we also tested the WGPFR, approximate wrapped Gaussian process functional regression model, functional linear regression model on  (the mean structure) and wrapped Gaussian process Fréchet mean regression model. 

The RMSE are given in Table <ref>. Similar to the simulation study on a S^2, the WGPFR model achieves the best prediction results in both interpolation and extrapolation. Moreover, the time consuming of simulation study on S^2 is less than 10 minutes while that of shape space is about one hour, since the later has higher dimension. Specifically, there are 3 sets of estimated hyper-parameters for a wrapped Gaussian process on a S^2 while 160 sets of estimated hyper-parameters for a wrapped Gaussian process on a shape space.



Thus we also show a log pointwise predictive density and the results are in Table <ref>.






 §.§ Flight Trajectory Data



In this section, we test our model on a real data set. The earth is roughly a sphere and can be modelled as a copy of S^2. Certain data sets, for example hurricane trajectories, are therefore considered as a manifold-valued random curve <cit.>. Here we study flight trajectory data, shown in Figure <ref> in which the red curves represent flights from Shanghai to London by British Airways and the black curves represent flight trajectories of Eastern China Airlines between the same destinations (the data were downloaded from www.variflight.com on 2020). Therefore, these sets of trajectories can naturally be split into two batches and the model with common mean structure can be used. 



The original data includes time, height, speed, longitude and latitude of each flight at discrete time points. We select the position of each airplane as the response variable, which can be transformed onto S^2 using longitude and latitude; in addition, the company and time are regarded as non-functional and functional covariates respectively. Before we train the model, it is necessary to pre-process the raw data. In this step, we firstly set the takeoff time of each flight to be 0 and the landing time to be 1, excluding taxi time. 25 trajectories of each company were selected in which the number of observed data points in every flight is greater than 600. In order to obtain smooth manifold-valued curves from the data, some kernel smoothing functions with small bandwidth were applied to the longitude and latitude of the training data. For computational reasons, we choose every 6-th data points of each smoothed trajectory as training data (100 data points in each curve totally).

To model the mean structure for flight trajectory data, we use company as the batch-specific covariate. In practice, for Eastern China Airlines, the covariate is defined as 0; and for British Airways, the covariate is defined as 1. Estimation of the mean structure and covariance structure were described in Section <ref> and <ref>, respectively. The parameters of the basis functions in mean structure and hyper-parameters in the covariance structure were be updated iteratively, as described in Section <ref>. In addition, the initial values of hyper-parameters in convariance structure are drawn from a standard Gaussian distribution independently. Afterwards, the predictions of the WGPFR model has been compared to FLRM and WGFmR for the same flight trajectory data.

The overall prediction of WGPFR model should be better than the other models, since the FLRM model only learns from mean structure and ignores the dependent error. In addition, the mean function of WGFmR is a point on S^2 and the prediction should have significant error when data are far away from the Fréchet mean. This is verified by numerical results in Table <ref>. 

The training data and test data are generated in the following way. We randomly selected another flight trajectory of British Airways which had >600 time observations. After the same pre-processing steps, 50 data points were added into the training data and the remaining 50 data points are used as the test data. In order to test the performance of interpolation and extrapolation for the real data set, we form the Type 1 prediction by randomly choosing these 50 data points and form the forecasting by selecting the first 50 points in the trajectory as training data. In addition, we test the capability for short-term forecasting and long-term forecasting by supposing different test data. Specifically, the 50-th to 60-th points and the 50-th to final data points of the flight trajectory can be selected respectively to form two scenarios of prediction which are denoted as short-term and long-term, respectively.

The predictive accuracy of WGPFR compared to FLR and WGFmR on the flight trajectory data is show in Table <ref>, which shows the mean of  root-mean-squared-error of 20 repeated simulations to reduce the effects of random seeds. We can see that the WGPFR model outperforms FLR and WGFmR for interpolation. For the short-term prediction, the RMSE is much smaller than that of the long-term prediction. As mentioned previously, the reason is that GPR provides little extra information when the test data are distant from the training data. The prediction of WGFmR is less accurate since the mean structure (Fréchet mean) is only a manifold-valued point and the test data set are not close to that point.







§ CONCLUSION AND DISCUSSION


In this paper we studied a novel functional regression model within a probabilistic framework for nonlinear manifold-valued response variables and real-valued covariates. Specifically, a wrapped Gaussian process functional regression (WGPFR) model is proposed to model the regression relationship into mean structure and covariance structure. For the mean structure, we proposed a functional version of geodesic regression model for the manifold-valued functional response and scalar covariates, where the inference of the mean structure is based on the intrinsic Fréchet mean function and the traditional functional linear regression  model is used on the tangent vectors of the manifold-valued data at the intrinsic Fréchet mean function to estimate the scalar covariates in an ambient Euclidean space. A wrapped Gaussian process prior is proposed to model the covariance structure in the tangent space at the mean function of each batch, where inference is conducted in a Bayesian way.Furthermore, an iterative approach based on a gradient descent algorithm or variational gradient descent algorithm are also applied to update mean structure and covariance structure efficiently. The mean structure captures the relation between functional dependent variable and batch-based scalar independent variable. Meanwhile, the covariance structure models the nonlinear concurrent relationship between functional output and multiple functional covariates. This idea in WGPFR model avoids the curse of dimensionality by using multiple functional covariates as input for functional regression, which promotes the flexibility of this method. 
 
Future research endeavors could encompass various extensions and enhancements to the proposed WGPFR model. Initially, the WGPFR model, as delineated in this paper, is predicated on the assumption of Riemannian manifolds with infinite injectivity radius. The potential exists to adapt this model to other manifolds that may not conform to this assumption.

Secondly, the current paper presumes the independence of Gaussian processes across different dimensions. Future investigations could explore methodologies as suggested by  <cit.> and <cit.>, where Gaussian processes exhibit dependence. However, such an approach would necessitate the development of additional computational efficiencies.

Thirdly, the proposed model could be applied to other intriguing real data sets. For instance, the performance of our model could be evaluated on medical imaging data. Repeated measurements of functional magnetic resonance imaging (fMRI), for example, could be construed as data residing in a specific manifold. Furthermore, the recurrently scanned shapes of the corpus callosum from diverse patients could be modeled to predict the likelihood of Alzheimer’s disease based on certain relative factors <cit.>. In addition, consideration could be given to non-Gaussian data as <cit.> and alternative methods to define mean and covariance structures.

The convergence of parameter estimation is substantiated in Theorem <ref> and <ref>. However, the assurance of convergence for the iterative optimization Algorithm <ref> remains an open question for future exploration. <cit.> demonstrated that, given suitable initial parameters, the estimation error of conditional maximization problems is confined within the bounds of statistical error and optimization error with a high degree of probability, potentially laying the groundwork for convergence. Nonetheless, it is crucial to acknowledge the distinction between graphical spaces and Riemannian manifolds. The approach proposed by <cit.>, which establishes that the convergence rate of GPR is optimal under certain preconditions and is upper bounded based on the smoothness of the correlation function in Euclidean space, could potentially be extrapolated to Riemannian manifolds to facilitate the convergence of covariance structure estimation. This proposition, however, requires rigorous validation.



§ ACKNOWLEDGMENTS

JQS’s work is supported by funds of National Key R&D Program of China (2023YFA1011400), the National Natural Science Foundation of China (No. 12271239) and Shenzhen Fundamental Research Program (No. 20220111). Chao Liu is supported in part by China Postdoctoral Science Foundation, No.2021M691443, No.2021TQ0141 and SUSTC Presidential Postdoctoral Fellowship.










§ APPENDIX



 §.§ Comparison of A Model with or without A Manifold Structure

In Kendall's shape space, the prediction of a regression model without a manifold structure loses the shape framework while a model with manifold structure still keep it.







 §.§ Examples of Riemannian Manifolds


Sphere S^n:
The two-dimensional sphere S^2={x∈ℝ^3:x=1} is a simple example of a Riemannian manifold. 
It is easy to show that the tangent vectors at a point p∈ S^2 are the vectors which are orthogonal to p. 
A Riemannian metric is inherited from the ambient Euclidean metric on ℝ^3, and it exactly the Euclidean inner product between tangent vectors. 
It is easy to show, using this metric, that the shortest curve passing between two points p,q∈ S^2 is a the great circle, which is contained within the intersection of S^2 and a plane containing p, q and the origin. 
The geodesic is unique exactly when one such plane exists, or in other words, when p and q are not antipodal. 
The formula of geodesic distance is given by

    d_(p, q)=cos^-1(p^T q).

The formula of exponential map is

    (p, v)=cos(‖ v‖)p + sin(‖ v‖)v/‖ v‖, for all  v∈ T_p={v:p^Tv=0}.

The logarithm map is defined for all pairs p,q which are not antipodal:

    (p,q)=u/‖ u‖d_(p,q),  where u=q-(p^Tq)p, for all p+q≠ 0.





































These concepts can be to the hyper-sphere S^n, where n>2. 



Kendall's Shape Space:
As a well-developed probability space, Kendall's shape space provides a very useful space for illustrate new theory and applications in statistics and machine learning. 
In 2-dimensional shape analysis <cit.>, a landmark is a point that can be unambiguously identified on all objects in a population. Landmarking is the process of identifying and labelling a set of k landmarks on each object in an analysis.
This process generates a 2k dimensional (real) vector, which can also be represented in a k-dimensional complex space ℂ^k.
The shape of an object is the equivalence class of the object when translation, rotation and scaling of the object are removed. 
Translation is removed by moving the centroid of each object to be origin point. The landmark vectors then lie in a complex subspace V={(z_i)∈ℂ^k|∑_i=1^k z_i=0}, which is a copy of ℂ^k-1. Two configurations in ℂ^k-1 are equivalent if they are related by scaling and rotation. The set of equivalence classes can be show to be equal to complex projective space Cℙ^k-2 and is known as Kendall's shape space <cit.>. The following formulae of exp map, log map on Kendall's shape space are from <cit.>.





Analogously to the sphere, we write down formulae which specify the Riemannian geometry on Kendall's shape space. Suppose p and q are two data points in a general Cℙ^k-2, the geodesic distance can be computed by

    d_(p, q) = arccos|p^*q|/pq

Note that this expression is invariant under multiplication of p and q by arbitrary non-zero complex constants. 
The inverse exponential map is given by

    (p,q)=θ·(q-π_pq)/ q-π_pq, θ=arccos|⟨ p,q⟩|/pq, π_pq=p·⟨ p,q⟩/p^2


The exponential map is defined by

    (p,v)=cosθ· p + psinθ/θv, θ=(p,q)


The tangent space T_p of a point p has the same dimensionality as the Kendall's shape space.


 §.§ Proofs


In order to prove the theorems in Section 3, we need the following assumptions for the Riemannian manifold  and  the manifold-valued data y(t).

	
  * (C1) ℳ is a closed Riemannian submanifold of a Euclidean space ℝ^d_0, with geodesic distance d_ℳ induced by the Euclidean metric.
	
  * (C2) Trajectories y(t) are continuous for t ∈𝒯 almost surely.
	
  * (C3) For all t ∈𝒯, the Fréchet mean μ_0(t) and sample  Fréchet mean μ̂_0(t) exist and are unique  almost surely.
	
  * (C4) Almost surely, trajectories y(t) lie in a compact set S_t⊂ B_ℳ(μ_0(t), r) for t ∈𝒯, where B_ℳ(μ_0(t), r) ⊂ℳ is an open ball centered at μ_0(t) with radius r<inf _t ∈𝒯inj_μ_0(t).
	
  * (C5) For any ϵ>0,
	
    inf _t ∈𝒯inf _p: d_ℳ(p, μ_0(t))>ϵ M(p, t)-M(μ_0(t), t)>0 .





 §.§ Proof of Theorem 1

 
Proof. Recall the notations and the multiple functional regression model V_m(t)=W_m(t)+e_m(t), denote the vector of the realizations at  t_i, m=1,·,M, i=1,⋯,n of V_m(t),W_m(t) and e_m(t) as V, W and E, respectively. Then we have

    vec(B̂)-vec(B)   =	((Φ⊗ U)^⊤ (Φ⊗ U))^-1(Φ⊗ U)^⊤ vec(V̂)-vec(B)
       =((Φ⊗ U)^⊤ (Φ⊗ U))^-1(Φ⊗ U)^⊤(vec(W+E)+vec(V̂-V))-vec(B)
       =((Φ⊗ U)^⊤ (Φ⊗ U))^-1(Φ⊗ U)^⊤((Φ⊗ U)vec(B)+vec(E)+vec(V̂-V))-vec(B)
       =((Φ⊗ U)^⊤ (Φ⊗ U))^-1(Φ⊗ U)^⊤ vec(E)
       +((Φ⊗ U)^⊤ (Φ⊗ U))^-1(Φ⊗ U)^⊤ vec(V̂-V)
       =̇A_1+A_2

For the first term, as e_m(t) are assumed to be independent to the covariate u_m for m=1,⋯,M, so that E(A_1)=0, Var(A_1)=O_p(1/M). Therefore A_1=o(1) as M→∞.

For the second term, note that 

    V̂_m(t)-V_m(t)=(μ̂_0(t),y_m(t))-(μ_0(t),y_m(t))≤ d_(μ̂_0(t),μ_0(t))=O_p(1/√(M))

where the last term has been proofed in <cit.> under the above conditions (C1)-(C5).
So the second term also has A_2=o(1) as  M→∞.  ▪
 

 


 §.§ Proof of Theorem 2

 
Before proofing the Theorem 2, we first introduce some notation to simplify the proof.
Recall that τ_md=(τ_md(t_m1),⋯,τ_md(t_mn)), for simplicity, we omit the subscript md and denote the n observation at the m-th batch and the d-th dimension as z_n=(z_1,⋯,z_n)=̇(τ_md(t_m1),⋯,τ_md(t_mn)), and denote the corresponding covariates 
X_n=(x_1,⋯,x_n) with x_i∈𝒳 which are independently drawn from a distribution 𝒰(x). Denote the observed data 𝒟_n={(x_i,z_i),i=1,⋯,n}. 
 
 
Let τ_0(·) be the true underlying function. Assume that the underlying process τ(·) ∼GP(0, K(·, · ; θ)) where all the subscript md are omitted. Denote

    p_g p(z_n)    =∫_ℱ p(z_1, …, z_n|τ(X_n)) d p_n(τ) 
    
     	p_0(z_n)    =p(z_1, …, z_n|τ_0(X_n))

then p_g p(z_n) is the Bayesian predictive distribution of z_n based on the GPR model. Note that p_n(τ) depends on n since the hyperparameters of τ(·) is estimated from the data.
 
Proof. 	
It suffices to show 
    1/n E_X_n(D[p_0(z_n),   p_g p(z_n)]) → 0    as  n →∞,

Note that

    D[p_0(z_n),  p_g p(z_n)]    =∫_𝒵^n p_0(z_1, ⋯, z_n) logp_0(z_1, ⋯, z_n)/  p_g p(z_1, ⋯, z_n) d z_1⋯ d z_n
       =∫_𝒵^n p_0(z_1, ⋯, z_n)[-log    p_g p(z_1, …, z_n)+log p_0(z_1, …, z_n)] d z_1⋯ d z_n .

It suffices to give an upper bound for the term -log   p_g p(z_1, …, z_n)+log p_0(z_1, …, z_n). 


Let ℋ be the Reproducing Kernel Hilbert Space (RKHS) associated with the
covariance function k(·, · ; θ), and ℋ_n the span of {k(·, x_i ; θ)}, i.e. ℋ_n={f(·): f(x)=
∑_i=1^nα_i k(x, x_i ; θ), for any .α_i∈ℝ}. 
Assume the true underlying function τ_0∈ℋ_n,
then τ_0(·) can be expressed as

    τ_0(·)=∑_i=1^nα_i k(·, x_i ; θ) ≜ K(·) α

where K(·)=(k(·, x_1 ; θ), …, k(·, x_n ; θ)) and α=(α_1, …, α_n)^T . By the properties of
RKHS,τ_0_k^2=α^TC_n nα, and (τ_0(x_1), …, τ_0(x_n))^T=C_n nα, where C_n n=(k(x_i, x_j ; θ)) is
the covariance matrix over x_i, i=1, …, n.

Let P and P̅ be any two measures on ℱ, then it yields by Fenchel-Legendre duality relationship that, for any functional g(·) on ℱ,

    E_P̅[g(τ)] ≤log E_P[e^g(τ)]+D[P̅, P]

Let g(τ)=log p(z_1, …, z_n|τ) for any z_1, …, z_n in 𝒵 and τ∈ℱ, let P be the measure induced by G P(0, k(·, · ; θ̂_n)), hence its finite dimensional distribution at z_1, …, z_n is p̃(z_1, …, z_n)=N(0, Ĉ_n n), and

    E_P[e^g(τ)]=E_P[p(z_1, …, z_n|τ)]=∫_ℱ p(z_1, ⋯, z_n|τ) d p_n(τ)=p_g p(z_n)

where Ĉ_n n is defined in the same way as C_n n but with θ being replaced by its estimator
θ̂_n.

Let P̅ be the posterior distribution of τ(·) on ℱ which has a prior distribution G P(0, k(·, · ; θ))
and normal likelihood ∏_i=1^n N(ẑ_i ; τ(x_i), σ^2), where

    ẑ_n ≜([ ẑ_1;   ⋮; ẑ_n ])=(C_n n+σ^2I) α

and σ^2 is a constant to be specified. In other words, we assume a model z=τ(x)+η
with η∼ N(0, σ^2) and τ(·) ∼ G P(0, k(·, · ; θ)), and ẑ_n   is a set of observations at x_1, …, x_n . Thus, P̅(τ)=p(τ|ẑ_n, X_n) is a probability measure on ℱ. Therefore,
by Gaussian process regression, the posterior of (τ_1, …, τ_n) ≜(τ(x_1), …, τ(x_n)) is

    p̅(τ_1, ⋯, τ_n)    ≜ p(τ_1, ⋯, τ_n|ẑ, X_n) 
       =N(C_n n(C_n n+σ^2I)^-1ẑ, C_n n(C_n n+σ^2I)^-1σ^2) 
       =N(C_n nα, C_n n(C_n n+σ^2I)^-1σ^2) 
       =N(C_n nα, C_n n B^-1)

where B=I+σ^-2C_n n.

It follows that

D[P̅, P]    =∫_ℱlog(d P̅/d P) d P̅
   =∫_R^np̅(τ_1, …, τ_n) logp̅(τ_1, …, τ_n)/p̃(τ_1, …, τ_n) d τ_1… d τ_n
   =1/2[log|Ĉ_n n|-log|C_n n|+log |B|+tr(Ĉ_n n^-1C_n n B^-1)+(C_n nα)^TĈ_n n^-1(C_n nα)-n] 
   =1/2[-log|Ĉ_n n^-1C_n n|+log |B|+tr(Ĉ_n n^-1C_n n B^-1)+τ_0_k^2.
   .+α^TC_n n(Ĉ_n n^-1C_n n-I) α-n]

On the other hand,

    E_P̅[g(τ)]=E_P̅[log p(z_1, …, z_n|τ)]=∑_i=1^n E_P̅[log p(z_i|τ(x_i))]

By Taylor's expansion, expanding log p(z_i|τ(x_i)) to the second order at τ_0(x_i) yields

    log p(z_i|τ(x_i))=   log p(z_i|τ_0(x_i))+.d[log p(z_i|τ(x_i))]/d τ(x_i)|_τ(x_i)=τ_0(x_i)(τ(x_i)-τ_0(x_i)) 
       +.1/2d^2[log p(z_i|τ(x_i))]/[d τ(x_i)]^2|_τ(x_i)=τ̃(x_i)(τ(x_i)-τ_0(x_i))^2

where τ̃(x_i)=τ_0(x_i)+λ(τ(x_i)-τ_0(x_i)) for some 0 ≤λ≤ 1.
For Gaussian probability density function, 
it follows that

    E_P̅[log p(z_i|τ(x_i))]= log p(z_i|τ_0(x_i))-σ^2/2Var[τ(x_i)]

so that

    E_P̅[log p(z_1, …, z_n|τ)]=log p_0(z_1, …, z_n)-σ^2/2tr(Var[τ(X_n)])=log p_0(z_1, …, z_n)-σ^2/2tr(C_n n B^-1)

Therefore,

-log p_g p(z_1, …, z_n)+log p_0(z_1, …, z_n) 
≤   -log E_P[e^g(τ)]+E_P̅[g(τ)]+σ^2/2tr(C_n n B^-1) 
≤    D[P̅, P]+σ^2/2tr(C_n n B^-1) 
=   1/2τ_0_k^2+1/2[-log|Ĉ_n n^-1C_n n|+log |B|+tr(Ĉ_n n^-1C_n n B^-1+σ^2 C_n n B^-1).
   .+α^TC_n n(Ĉ_n n^-1C_n n-I) α-n].

Since the covariance function is continuous in θ and θ̂_n→θ we have Ĉ_n n^-1C_n n-I→ 0
as n →∞. Therefore there exist some positive constants K and ϵ such that

    -log|Ĉ_n n^-1C_n n|<K
      α^TC_n n(Ĉ_n n^-1C_n n-I) α<K 
    tr((Ĉ_n n^-1C_n n-I) B^-1)<K

Thus

    -log p_g p(z_1, …, z_n)+log p_0(z_1, …, z_n) <1/2τ_0_k^2+1/2log |B|+3/2K

for any τ_0(·) ∈ℋ_n.
Taking infimum over τ_0 and applying Representer Theorem (see Lemma 2 in <cit.>) we obtain

    -log p_g p(z_1, …, z_n)+log p_0(z_1, …, z_n) ≤1/2τ_0_k^2+1/2log|I+σ^-2C_n n|+3/2K

for all τ_0(·) ∈ℋ. 

Therefore, we obtain that

    1/n E_X_n(D[p_0(z_n),  p_g p(z_n)]) ≤1/2 nτ_0_k^2+1/2 n E_X_n(log|I+σ^-2C_n n|)+3/2K/n→ 0

as n →∞.  ▪


0.2in





