# Integrated Model Selection and Scalability in Functional Data Analysis Through Bayesian Learning
## 2 Formulation  
### Generative Model  
In conventional fPCA, the data are assumed to be samples of functions that are elements of an appropriately smooth function space [1]. Using this assumption, data samples acquired at discrete points are typically interpolated to the continuum using tools such as splines, Fourier basis functions, or wavelets. In our work, we assume that the functions \(y_{i}:R^{M}\mapsto R\) are outcomes of an \(M\)-dimensional stochastic process. As in classicalfPCA, we assume that \(y_{i}\) is in a class of functions that can be approximated through a truncated, finite expansion, which is a weighted summation of \(K\) kernel functions \(\{\phi_{k}\}_{k=1}^{K}\):  
\[y_{i}(x)=\sum_{k=1}^{K}w_{ik}\phi_{k}(x), \tag{1}\]  
with \(w_{ik}\) being random variables, \(\phi_{k}(x)=\mathcal{K}(x,\mathbb{X}_{k})\), \(\mathcal{K}\) is the kernel function, and \(\mathbb{X}_{k}\) is the \(k\)-th location.  
Thus, \(y_{i}(x)\)s are realizations of a finite-dimensional stochastic process. Conventionally, the covariance operator of these functions is discretized, and the leading eigenfunctions form the estimated principal component loadings, following the Karhunen-Loeve theorem [1]. By contrast, we establish a flexible Bayesian framework of fPCA following the form of probabilistic PCA [26], where the principal subspace is identifiable up to an arbitrary rotation and does not enforce the orthogonality of the loadings. Nevertheless, it is straightforward to recover classical eigenfunctions from the final covariance estimation over an arbitrary grid in this Bayesian framework.  
The observed data are \(P\) independent, noisy samples of the functions \(\{y_{i}\}_{i=1}^{P}\) at index \(\{X_{i}\in R^{N_{i}\times M}\}_{i=1}^{P}\), where \(N_{i}\) is the number of measured samples for the \(i\)th function \(y_{i}\) and \(X_{in}\in R^{M}\) is the location of the \(n\)th measurement in the domain of the sample. The observations are \(\{Y_{i}\}_{i=1}^{P}\), where \(Y_{in}=y_{i}(X_{in})+E_{in}\), where \(E_{in}\) is white Gaussian noise of variance \(\sigma^{2}\).  
We also assume that the functions span a low-dimensional subspace of dimension \(J<<K\). We model this stochastically by assuming that the weights, \(w_{i}\in R^{K}\), are given by \(w_{ik}=\sum_{j=1}^{J}Z_{ij}W_{jk}+\tilde{Z}_{k}\), where \(W\in R^{J\times K}\) are the principal component loading coefficients and \(Z_{i}\in R^{J}\) are standard normal variables [26]. This model is therefore  
\[Y_{in}=\sum_{k=1}^{K}\left(\left(\sum_{j=1}^{J}\Bigl{(}Z_{ij}W_{jk}\Bigr{)}+ \tilde{Z}_{k}\right)\phi_{k}(X_{in})\right)+E_{in}=(Z_{i}W+Z)\Phi_{i\cdot n}+ E_{in}, \tag{2}\]  
where \(\Phi_{i\cdot n}=[\phi_{1}(X_{in}),\ldots,\phi_{K}(X_{in})]^{T}\) are the evaluations of the basis functions at the \(n\)-th index of the \(i\)-th sample function.  
The choice of the kernel family usually benefits from knowledge of the dataset's characteristics, such as the periodicity or domain geometry. Our framework is flexible across various kernel families, but we employ Gaussian kernels for both one-dimensional and multidimensional data by default, with the initial length-scale selection carried out through cross-validation, which will be refined through our sparse prior described below. To avoid disproportionately favoring larger length scales, we normalize each kernel's scaling coefficient using its square integral over the observational domain.