Computationally Efficient Kalman Filter Approaches for Fitting Smoothing Splines

Abstract

Smoothing spline models have shown to be effective in various fields (e.g., engineering and biomedical sciences) for understanding complex signals from noisy data. As nonparametric models, smoothing spline ANOVA (Analysis Of variance) models do not fix the structure of the regression function, leading to more flexible model estimates (e.g., linear or nonlinear estimates). The functional ANOVA decomposition of the regression function estimates offers interpretable results that describe the relationship between the outcome variable, and the main and interaction effects of different covariates/predictors. However, smoothing spline ANOVA (SS-ANOVA) models suffer from high computational costs, with a computational complexity of ON3 for N observations. Various numerical approaches can address this problem. In this chapter, we focus on the introduction to a state space representation of SS-ANOVA models. The estimation algorithms based on the Kalman filter are implemented within the SS-ANOVA framework using the state space representation, reducing the computational costs significantly.

1. Introduction

Smoothing spline ANOVA (SS-ANOVA) has been widely used in various applications [1, 2, 3]. The representer theorem enables an exact solution of regression function in SS-ANOVA models by minimizing a regularized function in a finite-dimensional space, even though the problem resides in a infinite-dimensional space. While SS-ANOVA models have strong theoretical properties, the estimation algorithms used to fit these models are computational intensive, with a computational complexity of $O\left(N^3\right)$ for datasets with $N$ observations. Numerous approaches have been developed to reduce the heavy computational costs of SS-ANOVA [4, 5, 6, 7]. For example, Kim and Gu (2004) proposed to select a $q\ll N$ basis functions from $N$ ones and reduced the computational complexity to $O\left({\mathit{Nq}}^2\right)$. Sun et al. (2021) synergistically combined asymptotic results with the smoothing parameter estimates based on randomly selected samples with sizes of $\stackrel{\sim }{N}\ll N$ to reduce the computational complexity of selecting smoothing parameters to $O\left({\stackrel{\sim }{N}}^3\right)$.

In this book chapter, we focus on the estimation approaches based on the Kalman filter. The Kalman filter was originally created to solve linear filtering and prediction problems used to generate simulations for the Apollo 11 project [8]. More recently, it has been implemented in a variety of engineering and biomedical fields [9, 10]. The Kalman filter is naturally used to fit state space models, methods that use recursive calculations on each observation entered one at a time and resulting in calculations on more accurate unknown variables after each iteration. The Kalman filter updates the state of the dynamic system given a new observation based on the state after the previous observation and the information gained from the new observation. This memory-less property and its simple recursive formulas make Kalman filter approaches computationally efficient, making them a useful tool for big data analytics. SS-ANOVA models can be reformulated to a state space representation, allowing computationally efficient Kalman filter-based model fitting and reducing computational costs to $O\left(N\right)$ for estimating univariate smoothing spline models [11]. An extension to the bivariate setting also significantly reduces the computational costs of SS-ANOVA models [12, 13].

Section 2 of this chapter will provide the theoretical background of SS-ANOVA models. Section 3 provides a brief background on state space models. The state space representation of SS-ANOVA models can be found in Section 4, along with a simulation study under the univariate setting in Section 5. Section 6 concludes this chapter.

2. Smoothing spline ANOVA models

We assume the data $\left(y_i,{\mathbf{x}}_i\right)$ and $i=1,2,\dots N$ are independent and identically distributed where $y_i\in Y\in \mathbb{R}$ is the outcome/response variable and ${\mathbf{x}}_i\in \mathcal{X}\in {\mathbb{R}}^d$ represents the covariates/predictors. A nonparametric model can then be written by

$$y_i=f\left({\mathbf{x}}_i\right)+e_i,$$E1

where $f$ is a function of covariates and $e_i\sim \mathcal{N}\left(0,{\sigma }^2\right)$ represents the random errors. For this nonparametric model, the structure of $f$ is not fixed and can be estimated by minimizing a penalized least squares score,

$$\frac{1}{N}\sum _{i=1}^N{\left(y_i-f\left({\mathbf{x}}_i\right)\right)}^2+\mathit{\lambda J}\left(f\right),$$E2

where the first term measures the goodness of fit of $f$, and the smoothing parameter $\lambda$ controls the trade-off between the goodness of fit and the roughness of $f$ measured by $J\left(f\right)$ [2, 3, 14]. SS-ANOVA models can also handle responses from exponential families and/or correlated responses. For readers who are interested in these topics, more examples of the model estimation and implementation exist (e.g., [2, 14]). The nonparametric estimation allows $f$ to vary in a high-dimensional (possibly infinite) space leading to more flexible results that can balance the bias-variance trade-off [2, 15]. The functional analysis of variance (ANOVA) is applied to the regression function $f$ to improve the interpretability of model estimates by decomposing the function into main and interaction effects of covariates. These main and interaction effects can be estimated in the corresponding subspaces of the reproducing kernel Hilbert space (RKHS), which is introduced in the next section.

2.1 ANOVA

2.1.1 Classical ANOVA

Classical ANOVA can be used to help to understand the decomposition of regression function in (1). We use a one-way classical ANOVA model as an example. The outcome $y_i$ can be modeled by $y_{\mathit{ij}}={\mu }_i+e_{\mathit{ij}}$, where ${\mu }_i$ is the mean treatment levels with $i=1,2,\cdots ,K_1$ and $j=1,2,\cdots ,K_2$. The terms in this model can be rewritten as

$$y_{\mathit{ij}}=\mu +{\delta }_i+e_{\mathit{ij}},$$E3

where $\mu$ is the overall mean effect and ${\delta }_i$ is the treatment effect. Side conditions are added to ensure the uniqueness of this decomposition. Now consider the univariate nonparametric function in (1). The regression function can be written as

$$f\left(x\right)=\mathit{Af}+\left(I-A\right)f=f_0+f_1$$E4

where $A$ is an averaging operator that averages the effect of $x$ and $I$ is the identity operator. We also need to add some side conditions for this decomposition to ensure the uniqueness of the decomposition of the regression function.

2.1.2 Functional ANOVA

The multivariate function $f\left(x_{⟨1⟩},x_{⟨2⟩},\dots ,x_{⟨d⟩}\right)$ on a $d$-dimensional product domain $\mathcal{X}={\prod }_{j=1}^d{\mathcal{X}}_j\in {\mathcal{R}}^d$ can be decomposed similarly to the classical ANOVA in the RKHS. The construction of the RKHS on ${\mathrm{\Pi }}_{j=1}^d{\mathcal{X}}_j$ is by taking the tensor product over the marginal domains ${\mathcal{X}}_j$. We need the following theorem to construct the tensor-product space.

Theorem 1.1 If $R_1\left(x_{⟨1⟩},{\stackrel{\sim }{x}}_{⟨1⟩}\right)$ is nonnegative definite on ${\mathcal{X}}_1$ and $R_2\left(x_{⟨2⟩},{\stackrel{\sim }{x}}_{⟨2⟩}\right)$ is nonnegative definite on ${\mathcal{X}}_2$, then $R\left(x_{⟨1⟩},x_{⟨2⟩}\right)=R_1\left(x_{⟨1⟩},{\stackrel{\sim }{x}}_{⟨1⟩}\right)R_2\left(x_{⟨2⟩},{\stackrel{\sim }{x}}_{⟨2⟩}\right)$ is nonnegative definite on $\mathcal{X}={\mathcal{X}}_1\times {\mathcal{X}}_2$.

Theorem 1.1 implies that the RKHS $\mathcal{H}$ on ${\mathrm{\Pi }}_{j=1}^d{\mathcal{X}}_j$ has the reproducing kernel $R={\mathrm{\Pi }}_{j=1}^d{R}_{⟨j⟩}$, where $R_{⟨j⟩}$ is the reproducing kernel for ${\mathcal{H}}_{⟨\gamma ⟩}$ on ${\mathcal{X}}_{⟨\gamma ⟩}$. Additionally, the Hilbert space ${\mathcal{H}}_{⟨j⟩}$ can be decomposed into ${\mathcal{H}}_{⟨j⟩}={\mathcal{H}}_{⟨j⟩\left(0\right)}\oplus {\mathcal{H}}_{⟨j⟩\left(1\right)}$, where ${\mathcal{H}}_{⟨j⟩\left(0\right)}$ is the null space and ${\mathcal{H}}_{⟨j⟩\left(1\right)}$ is the orthogonal complement to ${\mathcal{H}}_{⟨j⟩\left(0\right)}$. Then $\mathcal{H}={\otimes }_{j=1}^d{\mathcal{H}}_{⟨j⟩}$ can be decomposed as

$$\begin{array}{c}\mathcal{H}={\otimes }_{j=1}^d\left({\mathcal{H}}_{⟨j⟩\left(0\right)}\oplus {\mathcal{H}}_{⟨j⟩\left(1\right)}\right)\\ ={\oplus }_{\mathcal{S}}\left\{\left({\otimes }_{j\in \mathcal{S}}{\mathcal{H}}_{⟨j⟩\left(1\right)}\right)\otimes \left({\otimes }_{j\cancel{\in }\mathcal{S}}{\mathcal{H}}_{⟨j⟩\left(0\right)}\right)\right\}\\ ={\oplus }_{\mathcal{S}}{\mathcal{H}}_{\mathcal{S}}\end{array}$$E5

where $\mathcal{S}$ denotes all of the subsets of $\left\{1,\dots ,d\right\}$. The term ${\mathcal{H}}_{\mathcal{S}}$ has the reproducing kernel $R_{\mathcal{S}}\propto {\mathrm{\Pi }}_{j\in \mathcal{S}}{R}_{⟨j⟩\left(1\right)}$.

In general, the inner product in $\mathcal{H}$ can be specified as

$$J\left(f,g\right)=\sum _{j=1}^B{\theta }_j^{-1}{⟨f_j,g_j⟩}_j$$E6

where ${\theta }_j\ge 0$ are tuning parameters and $B$ is the number of smoothing parameters. The roughness penalty in (2) can be written in the form of (6). Then the reproducing kernel associated with (6) can be written as

$$R=\sum _{j=1}^B{\theta }_j{R}_j,$$E7

where $R_j$ is the reproducing kernel for the corresponding tensor product RKHS in (5).

2.2 An example of RKHS on ${\left[0,1\right]}^2$

We will use one example of an RKHS on ${\left[0,1\right]}^2$ to demonstrate the decomposition of RKHS. More examples of discrete and/or continuous domains can be found in Gu (2013) [2]. We consider the following tensor sum decomposition on $\left[0,1\right]$,

$$\begin{array}{l}{\mathcal{H}}_{⟨1⟩}=\left\{f:{\int }_0^1{\left(f^{\left(2\right)}\left(x\right)\right)}^2\mathit{dx}<\infty \right\}\\ =\left\{f:f\propto 1\right\}\oplus \left\{f:f\propto k_1\left(x\right)\right\}\\ \oplus \left\{f:{\int }_0^1\mathit{fdx}={\int }_0^1{f}^{\left(1\right)}\mathit{dx}={\int }_0^1{f}^{\left(2\right)}\mathit{dx}=0,f^{\left(2\right)}\in {\mathrm{ℒ}}_2\left[0,1\right]\right\}\\ ={\mathcal{H}}_{⟨1⟩\left(00\right)}\oplus {\mathcal{H}}_{⟨1⟩\left(01\right)}\oplus {\mathcal{H}}_{⟨1⟩\left(1\right)},\end{array}$$E8

where ${\mathcal{H}}_{⟨1⟩01}\oplus {\mathcal{H}}_{⟨1⟩1}$ forms the contrast in the one-way ANOVA decomposition, and the function $k_r$ is a scaled Bernoulli polynomial function with $k_r\left(x\right)=B_r\left(x\right)/r!$ [2, 16]. The RKHS has three reproducing kernels $R_{⟨1⟩00}\left(x,\stackrel{\sim }{x}\right)=1$, $R_{⟨1⟩01}\left(x,\stackrel{\sim }{x}\right)=k_1\left(x\right)k_1\left(\stackrel{\sim }{x}\right)$, and $R_{⟨1⟩1}\left(x,\stackrel{\sim }{x}\right)=k_2\left(x\right)k_2\left(\stackrel{\sim }{x}\right)-k_4\left(|x-\stackrel{\sim }{x}|\right)$.

Now consider the RKHS $\mathcal{H}$ on $\left[0,1\right]\times \left[0,1\right]$. Here $\mathcal{H}$ can be the tensor product spaces of ${\mathcal{H}}_{⟨1⟩}$ on $\left[0,1\right]$ and ${\mathcal{H}}_{⟨2⟩}$ on $\left[0,1\right]$. Based on the tensor sum decomposition in (8), we have

$$\begin{array}{l}\mathcal{H}={\mathcal{H}}_{⟨1⟩}\otimes {\mathcal{H}}_{⟨2⟩}\\ =\left({\mathcal{H}}_{⟨1⟩\left(00\right)}\oplus {\mathcal{H}}_{⟨1⟩\left(01\right)}\oplus {\mathcal{H}}_{⟨1⟩\left(1\right)}\right)\otimes \left({\mathcal{H}}_{⟨2⟩\left(00\right)}\oplus {\mathcal{H}}_{⟨2⟩\left(01\right)}\oplus {\mathcal{H}}_{⟨2⟩\left(1\right)}\right)\\ =\left({\mathcal{H}}_{⟨1⟩\left(00\right)}\otimes {\mathcal{H}}_{⟨2⟩\left(00\right)}\right)\oplus \left({\mathcal{H}}_{⟨1⟩\left(00\right)}\otimes {\mathcal{H}}_{⟨2⟩\left(01\right)}\right)\oplus \left({\mathcal{H}}_{⟨1⟩\left(01\right)}\otimes {\mathcal{H}}_{⟨2⟩\left(00\right)}\right)\\ \oplus \left({\mathcal{H}}_{⟨1⟩\left(01\right)}\otimes {\mathcal{H}}_{⟨2⟩\left(01\right)}\right)\oplus \left({\mathcal{H}}_{⟨1⟩\left(00\right)}\otimes {\mathcal{H}}_{⟨2⟩\left(1\right)}\right)\oplus \left({\mathcal{H}}_{⟨1⟩\left(1\right)}\otimes {\mathcal{H}}_{⟨2⟩\left(00\right)}\right)\\ \oplus \left({\mathcal{H}}_{⟨1⟩\left(01\right)}\otimes {\mathcal{H}}_{⟨2⟩\left(1\right)}\right)\oplus \left({\mathcal{H}}_{⟨1⟩\left(1\right)}\otimes {\mathcal{H}}_{⟨2⟩\left(01\right)}\right)\oplus \left({\mathcal{H}}_{⟨1⟩\left(1\right)}\otimes {\mathcal{H}}_{⟨2⟩\left(1\right)}\right),\end{array}$$E9

where the first four terms in (9) are in the null space ${\mathcal{H}}_{\left(\mathbf{0}\right)}$ and the remaining five terms are in the orthogonal complement (i.e., ${\mathcal{H}}_{\left(\mathbf{1}\right)}$). The reproducing kernel for the cubic spline (m=2) for each subspace can be found in Table 1.

Subspace	Reproducing Kernel
${\mathrm{\mathscr{H}}}_{⟨1⟩\left(00\right)}\otimes {\mathrm{\mathscr{H}}}_{⟨2⟩\left(00\right)}$	$1$
${\mathrm{\mathscr{H}}}_{⟨1⟩\left(01\right)}\otimes {\mathrm{\mathscr{H}}}_{⟨2⟩\left(00\right)}$	$k_1\left(x_{⟨1⟩}\right)k_1\left({\stackrel{\sim }{x}}_{⟨1⟩}\right)$
${\mathrm{\mathscr{H}}}_{⟨1⟩\left(01\right)}\otimes {\mathrm{\mathscr{H}}}_{⟨2⟩\left(01\right)}$	$k_1\left(x_{⟨1⟩}\right)k_1\left({\stackrel{\sim }{x}}_{⟨1⟩}\right)k_1\left(x_{⟨2⟩}\right)k_1\left({\stackrel{\sim }{x}}_{⟨2⟩}\right)$
${\mathrm{\mathscr{H}}}_{⟨1⟩\left(1\right)}\otimes {\mathrm{\mathscr{H}}}_{⟨2⟩\left(00\right)}$	$k_2\left(x_{⟨1⟩}\right)k_2\left({\stackrel{\sim }{x}}_{⟨1⟩}\right)-k_4\left(|x_{⟨1⟩}-{\stackrel{\sim }{x}}_{⟨1⟩}|\right)$
${\mathrm{\mathscr{H}}}_{⟨1⟩\left(1\right)}\otimes {\mathrm{\mathscr{H}}}_{⟨2⟩\left(01\right)}$	$\left[k_2\left(x_{⟨1⟩}\right)k_2\left({\stackrel{\sim }{x}}_{⟨1⟩}\right)-k_4\left(|x_{⟨1⟩}-{\stackrel{\sim }{x}}_{⟨1⟩}|\right)\right]k_1\left(x_{⟨2⟩}\right)k_1\left({\stackrel{\sim }{x}}_{⟨2⟩}\right)$
${\mathrm{\mathscr{H}}}_{⟨1⟩\left(1\right)}\otimes {\mathrm{\mathscr{H}}}_{⟨2⟩\left(1\right)}$	$\left[k_2\left(x_{⟨1⟩}\right)k_2\left({\stackrel{\sim }{x}}_{⟨1⟩}\right)-k_4\left(|x_{⟨1⟩}-{\stackrel{\sim }{x}}_{⟨1⟩}\right)\right]\left[k_2\left(x_{⟨2⟩}\right)k_2\left({\stackrel{\sim }{x}}_{⟨2⟩}\right)-k_4\left(|x_{⟨2⟩}-{\stackrel{\sim }{x}}_{⟨2⟩}|\right)\right]$

Table 1.

Subspaces and their corresponding reproducing kernels for the RKHS, $\mathrm{\mathscr{H}}$, on $\left[0,1\right]\times \left[0,1\right]$.

2.3 Estimation

The estimation algorithm of SS-ANOVA models relies on the following representer theorem.

Theorem 1.2 (Representer Theorem) There exist coefficient vectors $\mathit{\xi }={\left({\xi }_1,\dots ,{\xi }_M\right)}^{\prime }\in {\mathcal{R}}^M$ and $\mathbf{c}={\left(c_1,\dots ,c_N\right)}^{\prime }\in {\mathcal{R}}^N$ such that the minimizer of (2) in $\mathcal{H}={\mathcal{H}}_0\oplus {\mathcal{H}}_1$ has the following representation:

$$f\left(\mathbf{x}\right)=\sum _{m=1}^M{\xi }_m{\varphi }_m\left(\mathbf{x}\right)+\sum _{i=1}^N{c}_{i}R\left({\mathbf{x}}_i,\mathbf{x}\right),$$E10

where $\left\{{\varphi }_m,m=1,\dots ,M\right\}$ are the basis functions of the null space ${\mathcal{H}}_0$ and $R\left(\cdot ,\cdot \right)$ is the reproducing kernel of ${\mathcal{H}}_1$.

Taking into consideration model (1), the function $f$ can be estimated by minimizing (2). Using the representer theorem, the function $f$ can be written as

$$\mathbf{f}=\mathbf{S}\mathit{\xi }+\mathbf{Qc}$$E11

where $\mathbf{f}={\left(f\left({\mathbf{x}}_1\right),\dots ,f\left({\mathbf{x}}_N\right)\right)}^{\prime }$, $\mathbf{S}$ is a $N\times M$ matrix (e.g., $M=2$ for cubic spline) where the ${\left(i,j\right)}^{\mathit{th}}$ entry is ${\varphi }_j\left({\mathbf{x}}_i\right)$ and $\mathbf{Q}$ is a $N\times N$ matrix where the ${\left(i,j\right)}^{\mathit{th}}$ entry is $R\left({\mathbf{x}}_i,{\mathbf{x}}_j\right)$ of the form (7). Plugging (11) into (2), the penalized least squares can be written as

$$\frac{1}{N}{\left(\mathbf{y}-\mathbf{S}\mathit{\xi }-\mathbf{Qc}\right)}^{\prime }\left(\mathbf{y}-\mathbf{S}\mathit{\xi }-\mathbf{Qc}\right)+\lambda {\mathbf{c}}^{\prime }\mathbf{Qc}.$$E12

We differentiate (12) with respect to $\mathit{\xi }$ and $\mathbf{c}$ to obtain the linear system

$$\begin{array}{c}\left(\mathbf{Q}+\mathit{N\lambda }\mathbf{I}\right)\mathbf{c}+\mathbf{S}\mathit{\xi }=\mathbf{y},\\ {\mathbf{S}}^{\prime }\mathbf{c}=0.\end{array}$$E13

For given smoothing parameters, solving for $\mathbf{c}$ and $\mathit{\xi }$ provides the estimation for $f$. The selection of smoothing parameters for SS-ANOVA models is introduced below.

2.4 Selection of smoothing parameters

The smoothing parameters $\lambda /\mathit{\theta }$ balance the trade-off between the goodness of fit for $f$ and the roughness of $f$ in (2). Choosing the optimal smoothing parameters is data specific and should be performed prior to nonparametric regression analysis. Several smoothing parameter selection methods have been developed for the SS-ANOVA models [17]. Generalized cross-validation (GCV) is one of the most popular methods for selecting the optimal smoothing parameters $\lambda /\mathit{\theta }$ [18, 19].

To avoid overparameterization, let $\mathit{\lambda }={\left(\lambda /{\theta }_1,\dots \lambda /{\theta }_B\right)}^{\prime }$. The GCV score is defined as

$$V\left(\mathit{\lambda }\right)=\frac{N^{-1}{\mathbf{y}}^{\prime }{\left(I-A\left(\mathit{\lambda }\right)\right)}^2\mathbf{y}}{{\left[N^{-1}\mathit{tr}\left(I-\mathit{\alpha A}\left(\mathit{\lambda }\right)\right)\right]}^2}$$E14

where $A\left(\mathit{\lambda }\right)$ is symmetric matrix similar to the hat matrix in linear regression, and $\mathit{tr}\left(\cdot \right)$ represents trace. The parameter $\alpha \ge 1$ is a fudge factor [4]. When $\alpha =1$ it is the original GCV score. Larger $\alpha$’s yield smoother estimates. By default, we set $\alpha =1.4$. Then optimal smoothing parameters $\mathit{\lambda }$ can be chosen by minimizing the GCV score (14) using Newton-Raphson methods.

2.5 Computational complexity

In this section, we will discuss the computation complexity for calculating $\mathbf{c}$ and $\mathit{\xi }$ from (12). One requires $N^3/3+O\left(N^2\right)$ operations to obtain estimates of $\mathbf{c}$ and $\mathit{\xi }$ for the fixed smoothing parameters. In practice, the optimization of the smoothing parameter is also needed, which requires operations of $4{\mathit{BN}}^3/3+O\left(N^2\right)$, where $B$ is the number of smoothing parameters. Therefore, to minimize (12), the estimation algorithms have a computational complexity of $O\left(N^3\right)$. The following sections will discuss how the Kalman filter can be used to fit SS-ANOVA models, which reduces the computation complexity to $O\left(N\right)$.

3. State space models

State space methodology was traditionally used to study dynamic problems (e.g., space tracking settings) because the procedure allows for “real-time” updating as data are collected [20]. In this chapter, we use the linear Gaussian state space model as an example to introduce concepts of state space models using the Kalman filter approach. More applications of state space models can be found in a study by Douc, Moulines and Stoffer (2014) and Durbin and Koopman (2001) [21, 22]. A state space model consists of two equations: state equation and measurement equation. The state equation describes the dynamics of the state variables:

$${\mathbf{z}}_{t+1}={\mathbf{G}}_t{\mathbf{z}}_t+{\mathtt{\Psi }}_t{\mathit{\eta }}_t,{\mathit{\eta }}_t\stackrel{\mathrm{iid}}{\sim }\mathcal{N}\left(0,{\mathtt{\Delta }}_t\right),$$E15

where ${\mathbf{z}}_t$ is the $h\times 1$ state vector, ${\mathit{\eta }}_t$ is the $g\times 1$ disturbance vector with zero mean and a covariance matrix ${\mathtt{\Delta }}_t$, and ${\mathbf{G}}_t$ and ${\mathtt{\Psi }}_t$ are fixed design matrices of dimensions $h\times h$ and $h\times g$, respectively. The measurement equation shows the relationship between the observed variable and the unobserved state variable:

$${\mathbf{o}}_t={\mathtt{\Phi }}_t{\mathbf{z}}_t+{\mathit{\epsilon }}_t,{\mathit{\epsilon }}_t\stackrel{\mathrm{iid}}{\sim }\mathcal{N}\left(0,{\mathit{V}}_t\right),$$E16

where ${\mathbf{o}}_t$ is a $n\times 1$ vector with $n$ observations, ${\mathit{\epsilon }}_t$ is a $n\times 1$ disturbance vector with zero mean and a covariance matrix ${\mathbf{V}}_t$, and ${\mathtt{\Phi }}_t$ is a fixed design matrix of dimension $n\times h$. The initial state vector ${\mathbf{z}}_0$ is assumed to be normally distributed with mean ${\mu }_0$ and covariance matrix $\mathbf{P}$. The two vectors ${\mathit{\eta }}_t$ and ${\mathit{\epsilon }}_t$ are assumed to be mutually uncorrelated, that is,

$$\left(\begin{array}{c}{\mathit{\eta }}_t\\ {\mathit{\epsilon }}_t\end{array}\right)\stackrel{\mathrm{iid}}{\sim }\mathcal{N}\left[\mathbf{0},\left(\begin{array}{cc}{\mathtt{\Delta }}_t& \mathbf{0}\\ \mathbf{0}& {\mathbf{V}}_t\end{array}\right)\right].$$E17

These two vectors are also uncorrelated to the initial state vector ${\mathbf{z}}_0$.

The Kalman filter utilizes a set of recursive equations to estimate ${\mathbf{z}}_t$, given the observations ${\mathbf{O}}_t=\left\{{\mathbf{o}}_1,\dots ,{\mathbf{o}}_t\right\}$ at time $t$ and its error variance matrix ${\mathbf{P}}_t$ [13]. Define

$$\begin{array}{l}{\widehat{\mathbf{z}}}_t=\mathbb{E}\left[{\mathbf{z}}_t|{\mathbf{O}}_t\right],\\ {\mathbf{P}}_t=\mathbb{E}\left[\left({\mathbf{z}}_t-{\widehat{\mathbf{z}}}_t\right)\right({\left({\mathbf{z}}_t-{\widehat{\mathbf{z}}}_t\right)}^{\prime }\mid {\mathbf{O}}_t],\end{array}$$E18

where ${\widehat{\mathbf{z}}}_t$ is the Kalman filter estimation of ${\mathbf{z}}_t$, with ${\mathbf{z}}_0=\mathbb{E}\left[{\mathbf{z}}_0\right]={\mu }_0$, and ${\mathbf{P}}_0=\mathbf{P}$. From the state and measurement Eqs. (15) and (16), the estimated ${\mathbf{z}}_{t\mid t-1}$ and the covariance matrix given ${\mathbf{z}}_{t-1}$, and ${\mathbf{P}}_{t-1}$ become

$$\begin{array}{c}{\widehat{\mathbf{z}}}_{t\mid t-1}=\mathbb{E}\left[{\mathbf{z}}_t|{\mathbf{O}}_{t-1}\right]={\mathbf{G}}_t{\mathbf{z}}_{t-1},\\ {\mathbf{P}}_{t\mid t-1}=\mathbb{E}\left[\left({\mathbf{z}}_t-{\widehat{\mathbf{z}}}_{t\mid t-1}\right)\right({\left({\mathbf{z}}_t-{\widehat{\mathbf{z}}}_{t\mid t-1}\right)}^{\prime }\mid {\mathbf{O}}_{t-1}]={\mathbf{G}}_t{\mathbf{P}}_{t-1}{\mathbf{G}}_t^{\prime }+{\mathtt{\Psi }}_t{\mathtt{\Delta }}_t{\mathtt{\Psi }}_t^{\prime },\end{array}$$E19

and

$${\mathbf{o}}_{t\mid t-1}={\mathtt{\Phi }}_t{\widehat{\mathbf{z}}}_{t\mid t-1}.$$E20

Then the prediction error vector ${\mathbf{v}}_t$ is

$${\mathbf{v}}_t={\mathbf{o}}_t-{\mathbf{o}}_{t\mid t-1}={\mathbf{o}}_t-{\mathtt{\Phi }}_t{\widehat{\mathbf{z}}}_{t\mid t-1}={\mathtt{\Phi }}_t\left({\mathbf{z}}_t-{\widehat{\mathbf{z}}}_{t\mid t-1}\right)+{\mathit{\epsilon }}_t,$$E21

with the covariance matrix

$$\mathbb{E}\left[{\mathbf{v}}_t{\mathbf{v}}_t^{\prime }\right]={\mathtt{\Lambda }}_t={\mathtt{\Phi }}_t{\mathbf{P}}_{t\mid t-1}{\mathtt{\Phi }}_t^{\prime }+{\mathbf{V}}_t.$$E22

Using the facts of the joint distribution of ${\mathbf{z}}_t$ and ${\mathbf{v}}_t$, the Kalman filter estimator ${\widehat{\mathbf{z}}}_t=\mathbb{E}\left({\mathbf{z}}_t|{\mathbf{O}}_t\right)$ at time $t$ and its covariance matrix can be updated using

$$\begin{array}{c}{\widehat{\mathbf{z}}}_t={\widehat{\mathbf{z}}}_{t\mid t-1}+{\mathbf{P}}_{t\mid t-1}{\mathtt{\Phi }}_t^{\prime }{\mathtt{\Lambda }}_t^{-1}\left({\mathbf{o}}_t-{\mathtt{\Phi }}_t{\widehat{\mathbf{z}}}_{t\mid t-1}\right)\\ ={\mathbf{G}}_t{\mathbf{z}}_{t-1}+{\mathbf{P}}_{t\mid t-1}{\mathtt{\Phi }}_t^{\prime }{\mathtt{\Lambda }}_t^{-1}{\mathbf{v}}_t,\\ {\mathbf{P}}_t={\mathbf{P}}_{t\mid t-1}-{\mathbf{P}}_{t\mid t-1}{\mathtt{\Phi }}_t^{\prime }{\mathtt{\Lambda }}_t^{-1}{\mathtt{\Phi }}_t{\mathbf{P}}_{t\mid t-1}.\end{array}$$E23

We further define Kalman gain as

$${\mathbf{K}}_t={\mathbf{P}}_{t\mid t-1}{\mathtt{\Phi }}_t^{\prime }{\mathtt{\Lambda }}_t^{-1}.$$E24

Then the filtered estimate of ${\mathbf{z}}_t$ and its covariance matrix ${\mathbf{P}}_t$ is

$$\begin{array}{l}{\widehat{\mathbf{z}}}_t={\widehat{\mathbf{z}}}_{t\mid t-1}+{\mathbf{K}}_t{\mathbf{v}}_t,\\ {\mathbf{P}}_t=\left(\mathbf{I}-{\mathbf{K}}_t{\mathtt{\Phi }}_t\right){\mathbf{P}}_{t\mid t-1}.\end{array}$$E25

4. State space representation of SS-ANOVA models

Due to the Markov structure of SS-ANOVA models after reparameterization, the SS-ANOVA models can be represented by the state space models, allowing for efficient estimation by algorithms based on the Kalman filter. Wecker and Ansley (1983) showed the state space representation for univariate smoothing spline models [11]. Such an approach reduces the computational complexity of smoothing splines from $O\left(N^3\right)$ to $O\left(N\right)$. Based on the fast algorithm for the multivariate Kalman filter, Qin and Guo (2006) extended the univariate case to multivariate SS-ANOVA models [12]. The two-dimensional procedure was implemented, with computational complexity of $O\left(n_1{n}_2^3\right)$ for data of size $N=n_1{n}_2$. The extension to higher dimensions was also discussed. In this chapter, we focus on the univariate setting to demonstrate the procedure to derive the state space representation of smoothing splines.

4.1 Univariate setting

We can use the state space formulation to represent model (1) ($d=1$) and apply the Kalman filter algorithm to estimate the model parameters of smoothing spline models efficiently. Based on the pioneered work in Wahba (1978) [23], the univariate function $f\left(x\right)$ can be written in the following form

$$f\left(x\right)=\sum _{\nu =0}^{m-1}{\alpha }_{\nu }\frac{{\left(x-x_l\right)}^{\nu }}{\nu !}+\sqrt{\lambda }\sigma {\int }_{x_l}^x\frac{{\left(x-h\right)}^{m-1}}{\left(m-1\right)!}\mathit{dW}\left(h\right),$$E26

where the covariate $x\in \left[x_l,x_u\right]$ and $W\left(h\right)$ is a Wiener process with the unit dispersion parameter. When all $\alpha$’s have the diffuse prior distribution, the conditional expectation of $f\left(x\right)$ given all data is the function minimizing (2) with the smoothing parameter $1/\lambda$. The model (26) can be rewritten as

$$y_i={\mathbf{1}}^{\prime }\mathbf{U}\left(x_i\right)+e_i,i=1,\cdots ,N,$$E27

where ${\mathbf{1}}^{\prime }=\left[1,0,\cdots ,0\right]$ and $\mathbf{U}\left(x\right)={\left[U^{\left(m\right)}\left(x\right),\cdots ,U^{\left(1\right)}\left(x\right)\right]}^{\prime }$ is the $m$-dimensional stochastic process. In particular, we define

$$U^{\left(j\right)}\left(x\right)=\sum _{\nu =0}^{j-1}{U}^{\left(m-\nu \right)}\left(x_l\right)\frac{{\left(x-x_l\right)}^{\nu }}{\nu !}+\sqrt{\lambda }\sigma {\int }_{x_l}^x\frac{{\left(x-h\right)}^{m-1}}{\left(m-1\right)!}\mathit{dW}\left(h\right),j=m,\cdots ,1.$$E28

The vector $\mathbf{U}$ contains $U^{\left(m\right)}\left(x\right)$ and its first $\left(m-1\right)$ derivatives. Let ${\alpha }_{\nu }=U^{\left(m-\nu \right)}\left(x_l\right)$, and we can easily verify the model (27).

The state space formulation relies on the Markov structure of $\mathbf{U}\left(x\right)$, which is demonstrated below. We define a $m\times m$ matrix, ${\mathtt{\Gamma }}_m\left(x_b,x_a\right)$, for any $x_b$ and $x_a$ within the interval $\left[x_l,x_u\right]$ as

$${\mathtt{\Gamma }}_m\left(x_b,x_a\right)=\left(\begin{array}{cccc}1& \left(x_b-x_a\right)& \dots & \frac{{\left(x_b-x_a\right)}^{m-1}}{\left(m-1\right)!}\\ & 1& \dots & \frac{{\left(x_b-x_a\right)}^{m-2}}{\left(m-2\right)!}\\ & & & .\\ & & & .\\ & & & 1\end{array}\right).$$E29

We can easily verify

$${\mathtt{\Gamma }}_m\left(x_c,x_a\right)={\mathtt{\Gamma }}_m\left(x_c,x_b\right){\mathtt{\Gamma }}_m\left(x_b,x_a\right).$$E30

To show the Markov structure of $\mathbf{U}\left(x\right)$, we also define a $m\times 1$ random vector $\omega \left(x_b,x_a\right)={\left[{\omega }^{\left(1\right)}\left(x_b,x_a\right),\cdots ,{\omega }^{\left(m\right)}\left(x_b,x_a\right)\right]}^{\prime }$, where

$${\omega }^{\left(\nu \right)}\left(x_b,x_a\right)=\sqrt{\lambda }\sigma {\int }_{x_a}^{x_b}\frac{{\left(x_b-h\right)}^{\nu -1}}{\left(\nu -1\right)!}\mathit{dW}\left(h\right),\nu =1,\cdots ,m.$$E31

For any $\nu =1,\cdots ,m$, we have

$${\omega }^{\left(\nu \right)}\left(x_c,x_a\right)={\omega }^{\left(\nu \right)}\left(x_c,x_b\right)+\sum _{j=0}^{\nu -1}\frac{{\left(x_c-x_b\right)}^j}{j!}{\omega }^{\left(\nu -j\right)}\left(x_b,x_a\right).$$E32

Thus we have

$$\omega \left(x_b,x_a\right)={\mathtt{\Gamma }}_m\left(x_c,x_b\right)\mathit{\omega }\left(x_b,x_a\right)+\mathit{\omega }\left(x_c,x_b\right).$$E33

We now apply (30) and (33) to obtain the Markov structure of $\mathbf{U}$

$$\begin{array}{c}\mathbf{U}\left(x_b\right)={\mathtt{\Gamma }}_m\left(x_b,x_l\right)\mathbf{U}\left(x_l\right)+\omega \left(x_b,x_l\right)\\ ={\mathtt{\Gamma }}_m\left(x_b,x_a\right){\mathtt{\Gamma }}_m\left(x_a,x_l\right)\mathbf{U}\left(x_l\right)+{\mathtt{\Gamma }}_m\left(x_b,x_a\right)\mathit{\omega }\left(x_a,x_l\right)+\mathit{\omega }\left(x_b,x_a\right)\\ ={\mathtt{\Gamma }}_m\left(x_b,x_a\right)\mathbf{U}\left(x_a\right)+\mathit{\omega }\left(x_b,x_a\right).\end{array}$$E34

The Markov structure is the key to the state space representation of SS-ANOVA models. For $x_l=x_1\le \cdots \le x_n=x_u$, we have

$$y_i=\sum _{\nu =0}^{m-1}{\alpha }_{\nu }\frac{{\left(x_i-x_l\right)}^{\nu }}{\nu !}+{\mathbf{1}}^{\prime }\mathtt{\Omega }\left(x_i\right)+e_i,$$E35

as the measurement equation, where $\mathrm{\Omega }\left(x_i\right)={\left[{\mathrm{\Omega }}^{\left(m\right)}\left(x_i\right),\cdots ,{\mathrm{\Omega }}^{\left(1\right)}\left(x_i\right)\right]}^{\prime }$ and ${\mathrm{\Omega }}^{\left(\nu \right)}\left(x_i\right)={\omega }^{\left(\nu \right)}\left(x_i,x_l\right)$ for $\nu =1,\cdots ,m$ and $i=1,\cdots ,N$. From (33), we have the state equation

$$\mathtt{\Omega }\left(x_i\right)={\mathtt{\Gamma }}_m\left(x_i,x_{i-1}\right)\mathtt{\Omega }\left(x_{i-1}\right)+\mathit{\omega }\left(x_i,x_{i-1}\right),$$E36

for $i=1,\cdots ,N$. Given the parameters $\lambda$ and $\sigma$, the Kalman filtering and smoothing algorithms can be implemented to perform estimation for this state space representation.

5. Simulation studies

We used simulated data to compare the estimates of smoothing splines with those based on state space representation under the univariate setting. The following model was used to simulate $N=\mathrm{1,000}$ observations.

$$y_i=7sin\left({\pi }^{\ast }{x}_i\right)+e_i,$$E37

where $x_i=t_i/100$, $t_i=1,\cdots \mathrm{,1,000}$, and $e_i\sim \mathcal{N}\left(0,1\right)$. To apply the univariate SS-ANOVA model to the simulated data, we used the ssanova function in the gss package (version 2.2-3) [24]. The GCV algorithm was used to select the smoothing parameter of SS-ANOVA models. To implement the Kalman filtering algorithm to fit the smoothing splines, we used Eq. (35) as the measurement equation and Eq. (36) as the state equation. The parameters $\lambda$ and $\sigma$ were set to $0.01$ and $1$, respectively. We fitted the cubic spline (i.e., $m=2$) in the simulation studies. In the state Eq. (36), we have

$${\mathtt{\Gamma }}_2\left(x_i,x_{i-1}\right)=\left(\begin{array}{cc}1& \left(x_i-x_{i-1}\right)\\ 0& 1\end{array}\right),$$E38

for $i=2,\cdots \mathrm{,1,000}$. The ${\mathit{\nu \nu }}^{\text{'}}$th element of the variance matrix of $\mathit{\omega }\left(x_i,x_{i-1}\right)$ is

$${\mathit{\lambda \sigma }}^2\frac{{\left(x_i-x_{i-1}\right)}^{\nu +{\nu }^{\prime }-1}}{\left(\nu +{\nu }^{\prime }-1\right)\left(\nu -1\right)!\left({\nu }^{\prime }-1\right)!},$$E39

where $\nu ,{\nu }^{\prime }=1,\cdots ,m$. Given the above information, we utilized the fkf function in the FKF package (version 0.2.3) to implement the Kalman filtering and smoothing algorithms for the SS-ANOVA model. The details of iteration procedures are available in the help document of fkf, which is similar to the procedures described in Section 3. Figure 1 shows the similarity between the SS-ANOVA model fit with the generic algorithm from the gss package and the model fit based on the state space representations in (35) and (36).

Figure 1.

Comparison between the SS-ANOVA model fit and the model fit with the Kalman filter given $\lambda =0.01$ on simulated data.

6. Conclusions

In this chapter, we have introduced the theoretical foundation (e.g., representer theorem) and estimation algorithms of SS-ANOVA models. Given tensor product operations, the SS-ANOVA models can handle the multivariate data and study the main and interaction effects in the corresponding subspaces via functional ANOVA. The estimation algorithms of SS-ANOVA models need $O\left(N^3\right)$ operations, which might be prohibitive computationally for analyzing super large data. Utilizing the Markov structure of SS-ANOVA models, the Kalman filter can be used to fit SS-ANOVA models when reparameterized into a state space formulation [11]. This state space representation reduces the computational complexity from $O\left(N^3\right)$ to $O\left(N\right)$ for the univariate case, allowing SS-ANOVA models to be applicable to big data applications. Additional research has been done to extend this representation to the multidimensional setting [12]. For the two-dimensional data with the dimensions of $n_1$ and $n_2$, the SS-ANOVA models can be fitted with the computational complexity of $O\left(n_1{n}_2^3\right)$, where $N=n_1{n}_2$. Furthermore, we provided a simulated example to compare estimates from the state space representation and the estimates from the SS-ANOVA model for the univariate case.

Acknowledgments

Sun’s research was supported by the 18th Mile TRIF Funding from the University of Arizona and the NIH grant U19AG065169. LaFleur’s research was supported by the NIH grant U19AG065169.

Conflict of interest

The authors declare no conflict of interest.