Abstract
This paper addresses the problem of spline interpolations on \(\mathcal {P}\), the space of probability density functions when only a few observations \(p_i \in \mathcal {P}\) are available. Given a finite set of \(n+1\) distinct time instants \(t_i\) and corresponding data points \(p_i \in \mathcal {P}\), we consider the general problem of estimating a spline as a special regularized function \(\gamma \) on \(\mathcal {P}\) with \(\gamma (t_i)=p_i\). In particular, we focus on estimating missing data using smooth temporal splines to overcome the discrete nature of observations. In addition to generalizing splines on \(\mathcal {P}\) with minimal squared-norm of the acceleration, we give numerical schemes for solving \(C^1\) and \(C^2\) splines from data points \(p_i \in \mathcal {P}\). The two solutions are then shown to be computationally efficient, geometrically simpler, extensible, and can be transposed to other spaces and applications.
Similar content being viewed by others
Data Availability
All data generated or analyzed during this study are included in this published article.
References
Ay N, Jost J, Le H, Schwachhöfer L (2017) Information geometry. Springer, Cham
Bardsley J, Cui T, Marzouk Y, Wang Z (2020) Scalable optimization-based sampling on function space. SIAM J Sci Comput 42(2):A1317–A1347
Bogfjellmo G, Modin K, Verdier O (2018) Numerical algorithm for \(c^{2}\)-splines on symmetric spaces. SIAM J Numer Anal 56(4):2623–2647
Celledoni E, Eidnes S, Owren B, Ringholm T (2018) Dissipative numerical schemes on riemannian manifolds with applications to gradient flows. SIAM J Sci Comput 40(6):A3789–A3806
Cencov N (1982) Statistical decision rules and optimal inference. Translations of Mathematical Monographs, 53
Crouch P, Kun G, Silva Leite F (1999) The de Casteljau algorithm on lie groups and spheres. J Dyn Control Syst 5:397–429
de Casteljau P (1959) Outillages méthodes de calcul. Technical Report, André Citroën Automobiles, Paris
Friedrich T (1991) Die fisher-information und symplektische strukturen. Math Nachr 153:273–296
Kahl K, Rottmann M (2018) Least angle regression coarsening in bootstrap algebraic multigrid. SIAM J Sci Comput 40(6):A3928–A3954
Kwang-Rae K, Ian L, Huiling L, Severn K (2020) Smoothing splines on Riemannian manifolds with applications to 3d shape space. CoRR, abs/1801.04978
Lancaster P, Lancaster D, Šalkauskas K (1986) Curve and surface fitting: an introduction. Academic Press, Computational mathematics and applications
Lin L, StThomas B, Zhu H, Dunson D (2017) Extrinsic local regression on manifold-valued data. J Am Stat Assoc 112(519):1261–1273
Machado L, Silva Leite F (2006) Fitting smooth paths on Riemannian manifolds. Int J Appl Math Stat 06(4):25–53
Petersen A, Müller H (2019) Fréchet regression for random objects with Euclidean predictors. Ann Stat 47(2):691–719
Popiel T, Noakes L (2006) \(c^{2}\) spherical Bézier splines. Comput Aided Geom Des 23:261–275
Popiel T, Noakes L (2007) Bézier curves and c2 interpolation in riemannian manifolds. J Approx Theory 148(2):111–127
Rao C (1945) Information and accuracy attainable in the estimation of statistical parameters. Bull Calcutta Math Soc 37:81–89
Sabine L, Michael W (2019) Iterative solution of saddle-point systems from radial basis function (rbf) interpolation. SIAM J Sci Comput 41(3):A1706–A1732
Samir C, Absil P-A, Srivastava A, Klassen E (2012) A gradient-descent method for curve fitting on Riemannian manifolds. Found Comput Math 12:49–73
Shingel T (2009) Interpolation in special orthogonal groups. IMA J Numer Anal 29(3):731–745
Srivastava A, Klassen E (2016) Functional and shape data analysis. Springer Series in Statistics
Wallner J, Nava Yazdani E, Grohs P (2007) Smoothness properties of lie group subdivision schemes. Multiscale Model Simul 6(2):493–505
Acknowledgements
The authors would like to thank the handling editor Jose Alberto Cuminato and two referees for their detailed comments.
Author information
Authors and Affiliations
Corresponding author
Ethics declarations
Conflict of interest
The authors declare that there is no conflict of interest.
Additional information
Communicated by Jose Alberto Cuminato.
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Appendices
Appendix 1. \(C^{2}\) Bèzier spline on \(\mathbb {R}^{m}\)
Let us consider the Euclidean case \(\mathbb {R}^{m}\). Given a list of \((n+1)\) interpolation points \(p_{0},...,p_{n}\) and the control points \({\widehat{b}}_{i}^{+}\) and \({\widehat{b}}_{i}^{-}\) in the right and in the left of \(p_{i}\), \(i=0,...,n\). The \(C^{1}\) differentiability condition at knots \(p_{i}\) allows us to express control points \({\widehat{b}}_{i}^{+}\) in terms of \({\widehat{b}}_{i}^{-}\) as:
Hence, the task now is reduced to search only control points \({\widehat{b}}_{i}^{-}\), for \(i=1,...,n-1\), that generate the \(C^{1}\) Bézier spline \(\beta \) in \(\mathbb {R}^{m}\). Replacing the new optimization variables in problem (2) gives (24), which is merely the problem of minimization of the mean square acceleration of the Bézier curve \(\beta \) in the Euclidean space \(\mathbb {R}^{m}\). The optimal solution \(Y=[{\widehat{b}}_{1}^{-},...,{\widehat{b}}_{n-1}^{-}]^{T} \in \mathbb {R}^{ (n-1) \times m}\) of that problem is the unique solution of a tridiagonal linear system
where A is a tridiagonal sparse square matrix of size \((n-1)\times (n-1)\) with a dominant diagonal, C a matrix of size \((n-1)\times (n+1)\) and P the matrix of \(p_i\)’s of size \((n+1)\times m\) given by:
We may now write the \(C^{2}\) differentiability condition. It is obvious that with this \(C^{2}\) condition the position of the control points \({\widehat{b}}_{i}^{-}\) and \({\widehat{b}}_{i}^{+}\) that generate the curve \(\beta \) will be modified. Therefore, it is more convenient to use another notation. Let us denote by \(b^{-}_{i}\) and \(b^{+}_{i}\) the new control points on the left and on the right hand side of the interpolation point \(p_{i}\), for \(i=1,...,n-1\). Computing the acceleration of \(\beta \) on respective intervals and taking into account that \(\beta \) is \(C^{1}\), we shall replace \(b_{1}^{+}\) by (36), \(b_{i}^{+}\) by (37), and \(b_{n-1}^{+}\) by (38). We deduce that:
We see at once that points that will be modified by the additional \(C^{2}\) condition are \({\widehat{b}}_{i}^{-}\) and hence \({\widehat{b}}_{i}^{+}\), for \(i=2,...,n-1\). The point \({\widehat{b}}_{1}^{-}\) remains invariant and consequently it will be the case for \({\widehat{b}}_{1}^{+}\). According to the \(C^{1}\) differentiability condition ensured at the first step, one can take \(b_{1}^{-}={\widehat{b}}_{1}^{-}\), with \({\widehat{b}}_{1}^{-}\) is the first row of the matrix Y obtained as a solution of the optimization problem (24). However, the endpoint \(p_{n}\) is affected as we can deduce from Eq. (50). Nevertheless, it follows that giving the control point \(b_{1}^{-}\) allows us to find all the other control points including \(b_{2}^{-}\) with Eq. (48) and hence \(b_{2}^{+}\) with (37), then \(b_{i+1}^{-}\) for \(i=2,...,n-2\) with (49) and therefore \(b_{i}^{+}\), for \(i=3,...,n-2\) with (37) and \(b_{n-1}^{+}\) with (38).
Appendix 2. Proof of Theorem 3
We now prove Theorem 3. The proof is based on the following two results given in Popiel and Noakes (2007).
Lemma 1
Let \(\psi _{1} \in \mathcal {M}\).
-
(i)
\((d\varphi _{\psi _{1}})_{\psi _{2}}^{-1} = (d\varphi _{\psi _{1}})_{\varphi _{\psi _{1}}(\psi _{2})}, \; \text {for all} \; \psi _{2} \in \mathcal {M} \).
-
(ii)
\((d\varphi _{\psi _{1}})_{ \text {Exp}_{\psi _{1}}(H)} \circ (d \text {Exp}_{\psi _{1}})_{H} = - (d \text {Exp}_{\psi _{1}})_{-H} \; \text {for all} \; H \in T_{\psi _{1}} \mathcal {M}\).
Theorem 5
Let \(t \longrightarrow \sigma _{j}(t,V_{0},...,V_{j}) \) be the Bézier curve of order j on \(\mathcal {M}\) with a number of control points \(V_{j}\) for \(i=0,...,j\). Then, \(\sigma _{j}(t; V_{0},...,V_{j})\) satisfies:
-
i)
\(\frac{D}{\textrm{d}t}|_{t=0} \; {\dot{\sigma }}_{j}(t; V_{0},...,V_{j})= j(j-1)\varOmega _{0}\), where
$$\begin{aligned} {\varOmega _{0}:=} \left\{ \begin{array}{ll} {\dot{\alpha }}(0,V_{1},V_{2}), &{} if\,V_{0}=V_{1} \\ \ (d\text {Exp}_{V_{0}})^{-1}_{{\dot{\alpha }}(0,V_{0},V_{1})}\left( {\dot{\alpha }} (0,V_{1},V_{2})-{\dot{\alpha }}(1,V_{0},V_{1}) \right) , &{} if\,V_{0}\ne V_{1} \ \end{array}\right. \end{aligned}$$ -
ii)
\(\frac{D}{\textrm{d}t}|_{t=1} \; {\dot{\sigma }}_{j}(t; V_{0},...,V_{j})= j(j-1)\varOmega _{j}\), where
$$\begin{aligned} {\varOmega _{j}:=} \left\{ \begin{array}{ll} -{\dot{\alpha }}(0,V_{j-2},V_{j-1}), &{} if\,V_{j-1}=V_{j} \ \\ (d\text {Exp}_{V_{j}})^{-1}_{-{\dot{\alpha }}(1,V_{j-1},V_{j})}\left( {\dot{\alpha }}(0,V_{j-1},V_{j}) -{\dot{\alpha }}(1,V_{j-2},V_{j-1}) \right) , &{} if\,V_{j-1}\ne V_{j} \ \end{array}\right. \end{aligned}$$
We will exploit a modified form of the Theorem 5 to obtain the proof of the Theorem 3.
Proof of Theorem 3
Part i) follows from theorem 2. We now prove ii). The Bézier curve \(\sigma \) is \(C^{2}\) on \(\mathcal {M}\) if and only if it satisfies the \(C^{2}\) differentiability condition at joint points \(\tilde{\psi _{i}}\), for \(i=1,...,n-1\). At the point \(\tilde{\psi _{1}}\), this means:
Applying Theorem. 5 yields: \(\sigma \) is \(C^{2}\) on \(\tilde{\psi _{1}}\) if and only if \(\varOmega _{2}-3\varOmega _{0}=0\) with:
Since \(\beta _{1}\) is a \(C^{1}\) Bézier curve on \(T_{\psi _{1}}\mathcal {M}\), we get that \({\dot{\alpha }}(1,\chi _{1}^{-},\tilde{\psi _{1}}) = {\dot{\alpha }}(0,\tilde{\psi _{1}},\chi _{1}^{+})\). By Lemma 1, we have
It follows that
Hence, \(\varOmega _{2}-3\varOmega _{0}=0\) if and only if
Nevertheless \( \varphi _{\tilde{\psi _{1}}} \left( \alpha (t,\tilde{\psi _{1}},\chi _{1}^{+}) \right) = \alpha (1-t,\chi _{1}^{-},\tilde{\psi _{1}}), \; \forall t\in [0,1] \). Differentiate this identity with respect to t, we obtain
Accordingly, Eqn. ( ) becomes
Now, Lemma. 1 shows that
It follows that \((\textrm{d}\varphi _{\tilde{\psi _{1}}})_{\chi _{1}^{-}}^{-1} \left( {\dot{\alpha }}(0,\chi _{1}^{+},\chi _{2}^{-})\right) = \frac{1}{3} \left( {\dot{\alpha }}(1,\tilde{\psi _{0}},\chi _{1}^{-}) -4 {\dot{\alpha }}(0,\chi _{1}^{-}, \tilde{\psi _{1}}) \right) \).
Consequently, with the exponential map at the point \(\chi _{1}^{+}\), we get
The proof of Part iii) follows in much the same way as Part ii). \(\square \)
Rights and permissions
Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.
About this article
Cite this article
Adouani, I., Samir, C. Numerical algorithms for spline interpolation on space of probability density functions. Comp. Appl. Math. 42, 127 (2023). https://doi.org/10.1007/s40314-023-02262-5
Received:
Revised:
Accepted:
Published:
DOI: https://doi.org/10.1007/s40314-023-02262-5