supReferences Supplement

# Functional additive models on manifolds of planar shapes and forms

Almond Stöcker, Lisa Steyer, and Sonja Greven
School of Business and Economics, Humboldt-Universität zu Berlin

###### Abstract

The “shape” of a planar curve and/or landmark configuration is considered its
equivalence class under translation, rotation and scaling, its “form” its
equivalence class under translation and rotation while scale is preserved. We
extend generalized additive regression to models for such shapes/forms as
responses respecting the resulting quotient geometry by employing the squared
geodesic distance as loss function and a geodesic response function to map the
additive predictor to the shape/form space. For fitting the model, we propose
a Riemannian $L_{2}$-Boosting algorithm well suited for a potentially large
number of possibly parameter-intensive model terms, which also yields
automated model selection. We provide novel intuitively interpretable
visualizations for (even non-linear) covariate effects in the shape/form space
via suitable tensor-product factorization. The usefulness of the proposed
framework is illustrated in an analysis of 1) astragalus shapes of wild and
domesticated sheep and 2) cell forms generated in a biophysical model, as well
as 3) in a realistic simulation study with response shapes and forms motivated
from a dataset on bottle outlines.

Keywords: functional regression, boosting, shape analysis, tensor-product
model, visualization

##  1 Introduction

In many imaging data problems, the coordinate system of recorded objects is
arbitrary or explicitly not of interest. Statistical shape analysis (Dryden
and Mardia, 2016) addresses this point by identifying the ultimate object of
analysis as the shape of an observation, reflecting its geometric properties
invariant under translation, rotation and re-scaling, or as its form (or size-
and-shape) invariant under translation and rotation. This paper establishes a
flexible additive regression framework for modeling the shape or form of
planar (potentially irregularly sampled) curves and/or landmark configurations
in dependence on scalar covariates. A rich shape analysis literature has been
developed for 2D or 3D landmark configurations – presenting for instance
selected points of a bone or face – which are considered elements of Kendall’s
shape space (see, e.g. Dryden and Mardia, 2016). In many 2D scenarios,
however, observed points describe a curve reflecting the outline of an object
rather than dedicated landmarks (Adams et al., 2013). Considering outlines as
images of (parameterized) curves shows a direct link to functional data
analysis (FDA, Ramsay and Silverman, 2005) and, in this context, we speak of
functional shape/form data analysis. As in FDA, functional shape/form data can
be observed on a common and often dense grid (regular/dense design) or on
curve-specific often sparse grids (irregular/sparse design). While in the
regular case, analysis often simplifies by treating curve evaluations as
multivariate data, more general irregular designs gave rise to further
developments in sparse FDA (e.g. Yao et al., 2005; Greven and Scheipl, 2017),
explicitly considering irregular measurements instead of pre-smoothing curves.
To the best of our knowledge, we are the first to consider irregular/sparse
designs in the context of functional shape/form analysis.

Shapes and forms are examples of manifold data. Petersen and Müller (2019)
propose “Fréchet regression” for random elements in general metric spaces,
which requires estimation of a (potentially negatively) weighted Fréchet mean
for each covariate combination. Their implicit rather then explicit model
formulation renders model interpretation difficult. More explicit model
formulations have been developed for the special case of a Riemannian
geometry. Besides tangent space models (Kent et al., 2001), extrinsic models
(Lin et al., 2017) and models based on unwrapping (Jupp and Kent, 1987;
Mallasto and Feragen, 2018), a variety of manifold regression models have been
designed based on the intrinsic Riemannian geometry. Starting from geodesic
regression (Fletcher, 2013), which extends linear regression to curved spaces,
these include MANOVA (Huckemann et al., 2010), polynomial regression (Hinkle
et al., 2014), smoothing splines (Kume et al., 2007), regression along
geodesic paths with non-constant speed (Hong et al., 2014), or kernel
regression (Davis et al., 2010) and Kriging (Pigoli et al., 2016). However,
mostly only one metric covariate or categorical covariates are considered,
possibly in hierarchical model extensions for longitudinal data (Muralidharan
and Fletcher, 2012; Schiratti et al., 2017). By contrast, Zhu et al. (2009);
Shi et al. (2009); Kim et al. (2014) generalize geodesic regression to
regression with multiple covariates focusing on symmetric positive-definite
(SPD) matrix responses. Cornea et al. (2017) develop a general generalized
linear model (GLM) analogue regression framework for responses in a symmetric
manifold and apply it to shape analysis. Recently, Lin et al. (2020) proposed
a Lie group additive regression model for Riemannian manifolds focusing on SPD
matrices rather than shapes.

In FDA, there is a much wider range of developed regression methods (see
overviews in Morris, 2015; Greven and Scheipl, 2017). Among the most flexible
models are functional additive models (FAMs) for (univariate) functional
responses (in contrast to FAMs with functional covariates (Ferraty et al.,
2011)) with different strategies existing to model a) response functions and
b) smooth covariate effects. For a), basis expansions in spline (Brockhaus et
al., 2015), functional principal component (FPC) bases (Morris and Carroll,
2006) or both (Scheipl et al., 2015) are employed as well as wavelets (Meyer
et al., 2015), sometimes directly expanding functions to model on coefficients
and sometimes expanding only predictions while keeping the raw measurements.
Other approaches effectively evaluate curves on grids or apply pre-smoothing
techniques instead (e.g., Jeon and Park, 2020). For b), again penalized spline
basis approaches are employed (Scheipl et al., 2015; Brockhaus et al., 2015),
or local linear/polynomial (Müller and Yao, 2008; Jeon et al., 2022) or other
kernel-based approaches (Jeon and Park, 2020; Jeon et al., 2021). The
different approaches come with different theoretical and practical advantages,
but similiarities such as regarding asymptotic behavior are also known from
scalar nonparametric regression (Li and Ruppert, 2008). Advantages of the
fully basis expansion based approach summarized in Greven and Scheipl (2017)
include its appropriateness for sparse irregular functional data and its
modular extensibility to functional mixed models (Scheipl et al., 2015; Meyer
et al., 2015) and non-standard response distributions (Brockhaus et al., 2015;
Stöcker et al., 2021). For bivariate or multivariate functional responses,
which are closest to functional shapes/forms but without invariances, Rosen
and Thompson (2009); Zhu et al. (2012); Olsen et al. (2018) consider linear
fixed effects of scalar covariates, the latter also allowing for warping. Zhu
et al. (2017); Backenroth et al. (2018) consider one or more random effects
for one grouping variable, linear fixed effects and common dense grids for all
functions. Volkmann et al. (2021) combine the FAM model class of Greven and
Scheipl (2017) with multivariate FPC analysis (Happ and Greven, 2018) to model
multivariate (sparse) functional responses.

This paper establishes an interpretable FAM framework for modeling the shape
or form of planar (potentially irregularly sampled) curves and/or landmark
configurations in dependence on scalar covariates, extending $L_{2}$-Boosting
(Bühlmann and Yu, 2003; Brockhaus et al., 2015) to Riemannian manifolds for
model estimation. The three major contributions of our regression framework
are: 1\. We introduce additive regression with shapes/forms of planar curves
and/or landmarks as response, extending FAMs to non-linear response spaces or,
vice versa, extending GLM-type regression on manifolds for landmark shapes
both to functional shape manifolds and to include (non-linear) additive model
effects. 2. We propose a novel Riemannian $L_{2}$-Boosting algorithm for
estimating regression models for this type of manifold response, and 3. a
visualization technique based on tensor-product factorization yielding
intuitive interpretations even of multi-dimensional smooth covariate effects
for practitioners. Although related tensor-product model transformations based
on higher-order SVD have been used, i.a., in control engineering (Baranyi et
al., 2013), we are not aware of any comparable application for visualization
in FAMs or other statistical models for object data. Despite our focus on
shapes and forms, transfer of the model, Riemannian $L_{2}$-Boosting, and
factorized visualization to other Riemannian manifold responses is intended in
the generality of the formulation and the design of the provided R package
manifoldboost (developer version on github.com/Almond-S/manifoldboost). The
versatile applicability of the approach is illustrated in three different
scenarios: an analysis of the shape of sheep astragali (ankle bones)
represented by both regularly sampled curves and landmarks in dependence on
categorical “demographic” variables; an analysis of the effects of different
metric biophysical model parameters (including smooth interactions) on the
form of (irregularly sampled) cell outlines generated from a cellular Potts
model; and a simulation study with irregularly sampled functional shape and
form responses generated from a dataset of different bottle outlines and
including metric and categorical covariates.

In Section 2, we introduce the manifold geometry of irregular curves modulo
translation, rotation and potentially re-scaling, which underlies the
intrinsic additive regression model formulated in Section 3. The Riemannian
$L^{2}$-Boosting algorithm is introduced in Section 4. Section 5 analyzes
different data problems, modeling sheep bone shape responses (Section 5.1) and
cell outlines (Section 5.2). Section 5.3 summarizes the results of simulation
studies with functional shape and form responses. We conclude with a
discussion in Section 6.

##  2 Geometry of functional forms and shapes

Riemannian manifolds of planar shapes (and forms) are discussed in various
textbooks at different levels of generality, in finite (Dryden and Mardia,
2016; Kendall et al., 1999) or potentially infinite dimensions (Srivastava and
Klassen, 2016; Klingenberg, 1995). Starting from the Hilbert space
$\mathcal{Y}$ of curve representatives $y$ of a single shape or form
observation, we successively characterize its quotient space geometry under
translation, rotation and re-scaling including the respective tangent spaces.
Building on that, we introduce Riemannian exponential and logarithmic maps and
parallel transports needed for model formulation and fitting, and the sample
space of (irregularly observed) functional shapes/forms.

To make use of complex arithmetic, we identify the two-dimensional plane with
the complex numbers, $\mathds{R}^{2}\cong\mathds{C}$, and consider a planar
curve to be a function $y:\mathds{R}\supset\mathcal{T}\rightarrow\mathds{C}$,
element of a separable complex Hilbert space $\mathcal{Y}$ with a complex
inner product $\langle\cdot,\cdot\rangle$ and corresponding norm $\|\cdot\|$.
This allows simple scalar expressions for the group actions of translation
$\operatorname{Trl}=\\{y\overset{\operatorname{Trl}_{\gamma}}{\longmapsto}y+\gamma\,\\!\text{1}:\gamma\in\mathds{C}\\}$
with $\\!\text{1}\in\mathcal{Y}$ canonically given by
$\\!\text{1}:t\mapsto\frac{1}{\|t\mapsto 1\|}$ the real constant function of
unit norm; re-scaling
$\operatorname{Scl}=\\{y\overset{\operatorname{Scl}_{\lambda}}{\longmapsto}\lambda\cdot(y-\text{0}_{y})+\text{0}_{y}:\lambda\in\mathds{R}^{+}\\}$
around the centroid $\text{0}_{y}=\langle\\!\text{1}\,,y\rangle\\!\text{1}$
(which we consider more natural than using 0, the zero element of
$\mathcal{Y}$, mostly chosen in the literature); and rotation
$\operatorname{Rot}=\\{y\overset{\operatorname{Rot}_{u}}{\longmapsto}u\cdot(y-\text{0}_{y})+\text{0}_{y}:u\in\mathds{S}^{1}\\}$
around $\text{0}_{y}$ with
$\mathds{S}^{1}=\\{u\in\mathds{C}:|u|=1\\}=\\{\exp(\omega\sqrt{\text{-}1}):\omega\in\mathds{R}\\}$
reflecting counterclockwise rotations by $\omega$ radian measure.
Concatenation yields combined group actions $G$ as direct products, such as
the rigid motions
$G=\operatorname{Trl}\times\operatorname{Rot}=\\{\operatorname{Trl}_{\gamma}\circ\operatorname{Rot}_{u}:\gamma\in\mathds{C},u\in\mathds{S}^{1}\\}\cong\mathds{C}\times\mathds{S}^{1}$
(see Supplement S.1.1 for more details). The two real-valued component
functions of $y$ are identified with the real part
$\operatorname{Re}\\!\left(y\right):\mathcal{T}\rightarrow\mathds{R}$ and
imaginary part
$\operatorname{Im}\\!\left(y\right):\mathcal{T}\rightarrow\mathds{R}$ of
$y=\operatorname{Re}\\!\left(y\right)+\operatorname{Im}\\!\left(y\right)\sqrt{\text{-}1}$.
While the complex setup is used for convenience, the real part of
$\langle\cdot,\cdot\rangle$ constitutes an inner product
$\operatorname{Re}\\!\left(\langle
y_{1},y_{2}\rangle\right)=\langle\operatorname{Re}\\!\left(y_{1}\right),\operatorname{Re}\\!\left(y_{2}\right)\rangle+\langle\operatorname{Im}\\!\left(y_{1}\right),\operatorname{Im}\\!\left(y_{2}\right)\rangle$
for $y_{1},y_{2}\in\mathcal{Y}$ on the underlying real vector space of planar
curves. Typically $\operatorname{Re}\\!\left(y\right),\
\operatorname{Im}\\!\left(y\right)$ are assumed square-intregrable with
respect to a measure $\nu$ and we consider the canonical inner product
$\langle y_{1},y_{2}\rangle=\int{y}^{\dagger}_{1}y_{2}d\nu$ where
${y}^{\dagger}$ denotes the conjugate transpose of $y$, i.e.
${y}^{\dagger}(t)=\operatorname{Re}\\!\left(y\right)(t)-\operatorname{Im}\\!\left(y\right)(t)\sqrt{\text{-}1}$
is simply the complex conjugate, but for vectors
$\mathbf{y}\in\mathds{C}^{k}$, the vector ${\mathbf{y}}^{\dagger}$ is also
transposed. For curves, we typically assume $\nu$ to be the Lebesgue measure
on $\mathcal{T}=[0,1]$; for landmarks, a standard choice is the counting
measure on $\mathcal{T}=\\{1,\dots,k\\}$.

The ultimate response object is given by the _orbit_ $[y]_{G}=\\{g(y):g\in
G\\}$ (or short $[y]$) of $y\in\mathcal{Y}$, the equivalence class under the
respective combined group actions $G$: with
$G=\operatorname{Trl}\times\operatorname{Rot}\times\operatorname{Scl}$,
${[y]=}[y]_{\operatorname{Trl}\times\operatorname{Rot}\times\operatorname{Scl}}=\\{\lambda
u\,y+\gamma\,\\!\text{1}:\lambda\in\mathds{R}^{+},u\in\mathds{S}^{1},\gamma\in\mathds{C}\\}$
is referred to as the _shape_ of $y$ and, for
$G=\operatorname{Trl}\times\operatorname{Rot}$,
${[y]=}[y]_{\operatorname{Trl}\times\operatorname{Rot}}=\\{uy+\gamma\,\\!\text{1}:u\in\mathds{S}^{1},\gamma\in\mathds{C}\\}$
as its _form_ or _size-and-shape_.
$\mathcal{Y}_{/G}=\\{[y]_{G}:y\in\mathcal{Y}\\}$ denotes the quotient space of
$\mathcal{Y}$ with respect to $G$. The description of the Riemannian geometry
of $\mathcal{Y}_{/G}$ involves, in particular, a description of the tangent
spaces $T_{[y]}\mathcal{Y}_{/G}$ at points $[y]\in\mathcal{Y}_{/G}$, which can
be considered local vector space approximations to $\mathcal{Y}_{/G}$ in a
neighborhood of $[y]$. For a point $q$ in a manifold $\mathcal{M}$ the tangent
vectors $\beta\in T_{q}\mathcal{M}$ can, i.a., be thought of as gradients
$\dot{c}(0)$ of paths
$c:\mathds{R}\supset(-\delta,\delta)\rightarrow\mathcal{M}$ at $0$ where they
pass through $c(0)=q$. Besides their geometric meaning, they will also play an
important role in the regression model, as additive model effects are
formulated on tangent space level. Choosing suitable representatives
$\widetilde{y}^{G}\in[y]_{G}\subset\mathcal{Y}$ (or short $\widetilde{y}$) of
orbits $[y]_{G}$, we use an identification of tangent spaces with suitable
linear subspaces $T_{[y]_{G}}\mathcal{Y}_{/G}\subset\mathcal{Y}$.

Form geometry: Starting with translation as the simplest invariance, an orbit
$[y]_{\operatorname{Trl}}$ can be one-to-one identified with its centered
representative $\widetilde{y}^{\operatorname{Trl}}=y-\langle
y,\\!\text{1}\rangle\,\\!\text{1}$ yielding an identification
$\mathcal{Y}_{/\operatorname{Trl}}\cong\\{y\in\mathcal{Y}:\langle
y,\\!\text{1}\,\rangle=0\\}$ with a linear subspace of $\mathcal{Y}$. Hence,
also $T_{[y]}\mathcal{Y}_{/\operatorname{Trl}}=\\{y\in\mathcal{Y}:\langle
y,\\!\text{1}\,\rangle=0\\}$. For rotation, by contrast, we can only find
local identifications with Hilbert subspaces (i.e. charts) around reference
points $[p]_{\operatorname{Trl}\times\operatorname{Rot}}$ we refer to as
“poles”. Moreover, we restrict to
$y,p\in\mathcal{Y}^{*}=\mathcal{Y}\setminus[\text{0}\,]_{\operatorname{Trl}}$
eliminating constant functions as degenerate special cases in the translation
orbit of zero. For each $[y]_{\operatorname{Trl}\times\operatorname{Rot}}$ in
an open neighborhood around $[p]_{\operatorname{Trl}\times\operatorname{Rot}}$
which can be chosen with
$\langle\widetilde{y}^{\operatorname{Trl}},\widetilde{p}^{\operatorname{Trl}}\rangle\neq
0$, $y$ can be uniquely rotation aligned to $p$, yielding a one-to-one
identification of the form $[y]_{\operatorname{Trl}\times\operatorname{Rot}}$
with the aligned representative given by
$\widetilde{y}^{\operatorname{Trl}\times\operatorname{Rot}}=\frac{\langle\widetilde{y}^{\operatorname{Trl}},\widetilde{p}^{\operatorname{Trl}}\rangle}{|\langle\widetilde{y}^{\operatorname{Trl}},\widetilde{p}^{\operatorname{Trl}}\rangle|}\widetilde{y}^{\operatorname{Trl}}=\underset{y^{\prime}\in[y]_{\operatorname{Trl}\times\operatorname{Rot}}}{\operatorname{argmin}\;}\|y^{\prime}-p\|$
(compare Fig. 1). While
$\widetilde{y}^{\operatorname{Trl}\times\operatorname{Rot}}$ depends on $p$,
we omit this in the notation for simplicity. All
$\widetilde{y}^{\operatorname{Trl}}$ rotation aligned to
$\widetilde{p}^{\operatorname{Trl}}$ lie on the hyper-plane determined by
$\operatorname{Im}\\!\left(\langle\widetilde{y}^{\operatorname{Trl}},\widetilde{p}^{\operatorname{Trl}}\rangle\right)=0$
(Figure 1), which yields
$T_{[p]}\mathcal{Y}_{/\operatorname{Trl}+\operatorname{Rot}}^{*}=\\{y\in\mathcal{Y}:\langle
y,\\!\text{1}\rangle=0,\ \operatorname{Im}\\!\left(\langle
y,p\rangle\right)=0\\}$ with normal vectors
$\zeta^{(1)}=\\!\text{1},\zeta^{(2)}=\sqrt{\text{-}1}\,\\!\text{1},\zeta^{(3)}=\sqrt{\text{-}1}\,p$.
Note that, despite the use of complex arithmetic,
$T_{[p]}\mathcal{Y}_{/\operatorname{Trl}\times\operatorname{Rot}}^{*}$ is a
real vector space not closed under complex scalar multiplication. The geodesic
distance of $[y]_{\operatorname{Trl}\times\operatorname{Rot}}$ to the pole
$[p]_{\operatorname{Trl}\times\operatorname{Rot}}$ is given by
$d([y]_{\operatorname{Trl}\times\operatorname{Rot}},[p]_{\operatorname{Trl}\times\operatorname{Rot}})=\|\widetilde{y}^{\operatorname{Trl}\times\operatorname{Rot}}-\widetilde{p}^{\operatorname{Trl}}\|=\underset{y^{\prime}\in[y]_{\operatorname{Trl}\times\operatorname{Rot}},p^{\prime}\in[p]_{\operatorname{Trl}\times\operatorname{Rot}}}{\operatorname{argmin}}\|y^{\prime}-p^{\prime}\|$.
It reflects the length of the shortest path (i.e. the geodesic) between the
forms and the minimum distance between the orbits as sets.

Shape geometry: To account for scale invariance in shapes
$[y]_{\operatorname{Trl}\times\operatorname{Rot}\times\operatorname{Scl}}$,
they are identified with normalized representatives
$\widetilde{y}^{\operatorname{Trl}\times\operatorname{Rot}\times\operatorname{Scl}}=\frac{\widetilde{y}^{\operatorname{Trl}\times\operatorname{Rot}}}{\|\widetilde{y}^{\operatorname{Trl}\times\operatorname{Rot}}\|}$.
Motivated by the normalization, we borrow the well-known geometry of the
sphere $\mathds{S}=\\{y\in\mathcal{Y}:\|y\|=1\\}$, where
$T_{p}\mathds{S}=\\{y\in\mathcal{Y}:\operatorname{Re}\\!\left(\langle
y,p\rangle\right)=0\\}$ is the tangent space at a point $p\in\mathds{S}$ and
geodesics are great circles. Together with translation and rotation
invariance, the shape tangent space is then given by
$T_{[p]}\mathcal{Y}_{/\operatorname{Trl}\times\operatorname{Rot}\times\operatorname{Scl}}^{*}=T_{[p]}\mathcal{Y}_{/\operatorname{Trl}\times\operatorname{Rot}}^{*}\cap
T_{p}\mathds{S}=\\{y\in\mathcal{Y}:\langle y,\\!\text{1}\rangle=0,\ \langle
y,p\rangle=0\\}$ with normal vector $\zeta^{(4)}=p$ in addition to
$\zeta^{(1)},\zeta^{(2)},\zeta^{(3)}$ above. The geodesic distance
$d([p]_{\operatorname{Trl}\times\operatorname{Rot}\times\operatorname{Scl}},[y]_{\operatorname{Trl}\times\operatorname{Rot}\times\operatorname{Scl}})=\arccos|\langle\widetilde{y}^{\operatorname{Trl}\times\operatorname{Rot}\times\operatorname{Scl}},\widetilde{p}^{\operatorname{Trl}\times\operatorname{Rot}\times\operatorname{Scl}}\rangle|$
corresponds to the arc-length between the representatives. This distance is
often referred to as Procrustres distance in statistical shape analysis.

We may now define the maps needed for the regression model formulation. Let
$\widetilde{y}$ and $\widetilde{p}$ be shape/form representatives of $[y]$ and
$[p]$ rotation aligned to the shape/form pole representative $p$. Generalizing
straight lines to a Riemannian manifold $\mathcal{M}$, geodesics
$c:(-\delta,\delta)\rightarrow\mathcal{M}$ can be characterized by their
“intercept” $c(0)\in\mathcal{M}$ and “slope” $\dot{c}(0)\in
T_{c(0)}\mathcal{M}$. The _exponential map_
$\operatorname{Exp}_{q}:T_{q}\mathcal{M}\rightarrow\mathcal{M}$ at a point
$q\in\mathcal{M}$ is defined to map $\beta\mapsto c(1)$ for $c$ the geodesic
with $q=c(0)$ and $\beta=\dot{c}(0)$. It maps $\beta\in T_{q}\mathcal{M}$ to a
point $\operatorname{Exp}_{q}(\beta)\in\mathcal{M}$ located
$d(q,\operatorname{Exp}_{q}(\beta))=\|\beta\|$ apart of the pole $q$ in the
direction of $\beta$. On the form space
$\mathcal{Y}_{/\operatorname{Trl}\times\operatorname{Rot}}$, the exponential
map is simply given by
$\operatorname{Exp}_{[p]{{}_{\operatorname{Trl}\times\operatorname{Rot}}}}(\beta)=\left[\widetilde{p}{{}^{\operatorname{Trl}\times\operatorname{Rot}}}+\beta\right]{{}_{\operatorname{Trl}\times\operatorname{Rot}}}$.
On the shape space
$\mathcal{Y}_{/\operatorname{Trl}\times\operatorname{Rot}{\times}\operatorname{Scl}}$,
identification with exponential maps on the sphere yields
$\operatorname{Exp}_{[p]{{}_{G}}}(\beta)=\left[\cos(\|\beta\|)\widetilde{p}{{}^{G}}+\sin(\|\beta\|)\frac{\beta}{\|\beta\|}\right]{{}_{G}}$
with $G=\operatorname{Trl}\times\operatorname{Rot}\times\operatorname{Scl}$.
In an open neighborhood $\mathcal{U}$, $q\in\mathcal{U}\subset\mathcal{M}$,
$\operatorname{Exp}_{q}$ is invertible yielding the
$\operatorname{Log}_{q}:\mathcal{U}\rightarrow T_{q}\mathcal{M}$ map from the
manifold to the tangent space at $q$. For forms, it is given by
$\operatorname{Log}_{[p]_{\operatorname{Trl}\times\operatorname{Rot}}}([y]_{\operatorname{Trl}\times\operatorname{Rot}})=\widetilde{y}^{\operatorname{Trl}\times\operatorname{Rot}}-\widetilde{p}^{\operatorname{Trl}\times\operatorname{Rot}}$
and, for shapes, by
$\operatorname{Log}_{[p]_{G}}([y]_{G})=d([p]_{G},[y]_{G})\frac{\widetilde{y}^{G}-\langle\widetilde{p}^{G},\widetilde{y}^{G}\rangle\widetilde{p}^{G}}{\|\widetilde{y}^{G}-\langle\widetilde{p}^{G},\widetilde{y}^{G}\rangle\widetilde{p}^{G}\|}$
with $G=\operatorname{Trl}\times\operatorname{Rot}\times\operatorname{Scl}$.
Finally, $\operatorname{Transp}_{q,q^{\prime}}:T_{q}\mathcal{M}\rightarrow
T_{q^{\prime}}\mathcal{M}$ parallel transports tangent vectors
$\varepsilon\mapsto\varepsilon^{\prime}$ isometrically along a geodesic
$c(\tau)$ connecting $q$ and $q^{\prime}\in\mathcal{M}$ such that the slopes
${\operatorname{Transp}_{q,q^{\prime}}(}\dot{c}(q){)}{=}\dot{c}(q^{\prime})$
are identified and all angles are preserved. For shapes,
$\operatorname{Transp}_{[{y}]_{G},[p]_{G}}(\varepsilon)=\varepsilon-\langle\varepsilon,{\widetilde{p}}^{G}\rangle\frac{{\widetilde{y}}^{G}+{\widetilde{p}}^{G}}{1+\langle{\widetilde{y}}^{G},{\widetilde{p}}^{G}\rangle}$,
with $G=\operatorname{Trl}\times\operatorname{Rot}\times\operatorname{Scl}$,
takes the form of the parallel transport on a sphere replacing the real inner
product with its complex analogue. For forms, it changes only the
$\operatorname{Im}\\!\left(\langle\varepsilon,\widetilde{p}\rangle\right)$
coordinate orthogonal to the real ${\widetilde{y}}$-$\widetilde{p}$-plane as
in the shape case, while the remainder of $\varepsilon$ is left unchanged as
in a linear space. This yields
$\operatorname{Transp}_{[{y}]{{}_{G}},[p]{{}_{G}}}\left(\varepsilon\right)=\varepsilon-\operatorname{Im}\\!\left(\langle\widetilde{p}{{}^{G}}/\|\widetilde{p}{{}^{G}}\|,\varepsilon\rangle\right)\frac{{\widetilde{y}^{G}}/\|{\widetilde{y}^{G}}\|+\widetilde{p}{{}^{G}}/\|\widetilde{p}{{}^{G}}\|}{1+\langle{\widetilde{y}^{G}}/\|{\widetilde{y}^{G}}\|,\widetilde{p}{{}^{G}}/\|\widetilde{p}{{}^{G}}\|\rangle}\sqrt{\text{-}1}$,
with $G=\operatorname{Trl}\times\operatorname{Rot}$, for form tangent vectors.
While equivalent expressions for the parallel transport in the shape case can
be found, e.g., in Dryden and Mardia (2016); Huckemann et al. (2010), a
corresponding derivation for the form case is given in Supplement S.1.2
including a discussion of the quotient space geometry in differential
geometric terms.

Based on this understanding of the response space, we may now proceed to
consider a sample of curves $y_{1},\dots,y_{n}\in\mathcal{Y}$ representing
orbits $[y_{1}],\dots,[y_{n}]$ with respect to group actions $G$. In the
functional case, with the domain $\mathcal{T}=[0,1]$, these curves are usually
observed as evaluations
$\mathbf{y}_{i}=(y_{i}(t_{i1}),\dots,y_{i}(t_{ik_{i}}))^{\top}$ on a finite
grid $t_{i1}<\dots<t_{ik_{i}}\in\mathcal{T}$ which may differ between
observations. In contrast to the regular case with common grids, this more
general data structure is referred to as irregular functional shape/form data.
To handle this setting, we replace the original inner product
$\langle\cdot,\cdot\rangle$ on $\mathcal{Y}$ by individual $\langle
y_{i},y_{i}^{\prime}\rangle_{i}={\mathbf{y}}^{\dagger}_{i}\mathbf{W}_{i}\mathbf{y}^{\prime}_{i}$
providing inner products on the $k_{i}$-dimensional space
$\mathcal{Y}_{i}=\mathds{C}^{k_{i}}$ of evaluations
$\mathbf{y}_{i},\mathbf{y}_{i}^{\prime}$ on the same grid. The symmetric
positive-definite weight matrix $\mathbf{W}_{i}$ can be chosen to implement an
approximation to integration w.r.t. the original measure $\nu$ with a
numerical integration measure $\nu_{i}$ such as given by the trapezoidal rule.
Alternatively, $\mathbf{W}_{i}=\frac{1}{k_{i}}\mathbf{I}_{k_{i}}$ with
$k_{i}\times k_{i}$ identity matrix $\mathbf{I}_{k_{i}}$ presents a canonical
choice that is analog to the landmark case for $k_{i}\equiv k$. Moreover,
data-driven $\mathbf{W}_{i}$ could also be motivated from the covariance
structure estimated for (potentially sparse) $y_{1},\dots,y_{n}$ along the
lines of Yao et al. (2005); Stöcker et al. (2022). While this is beyond the
scope of this paper, potential procedures are sketched in Supplement S.7. With
the inner products given for $i=1,\dots,n$, the sample space naturally arises
as the Riemannian product
$\mathcal{Y}^{*}_{1/G}\times\dots\times\mathcal{Y}^{*}_{n/G}$ of the orbit
spaces, with the individual geometries constructed as described above.

Figure 1:  Left: Quotient space geometry: assuming $p$ and $y$ centered,
translation invariance is not further considered in the plot; given pole
representative $p$, we express $y=\frac{\operatorname{Re}\\!\left(\langle
p,y\rangle\right)}{\|p\|^{2}}p+\frac{\operatorname{Im}\\!\left(\langle
p,y\rangle\right)}{\|p\|^{2}}ip+(y-\frac{\langle
p,y\rangle}{\|p\|^{2}}p)\in\mathcal{Y}$ in its coordinates in $p$ and $ip$
direction, subsuming all orthogonal directions in the third dimension. In this
coordinate system, the rotation orbit $[y]_{\operatorname{Rot}}$ corresponds
to the dotted horizontal circle, and is identified with the aligned
$\widetilde{y}:=\widetilde{y}^{\operatorname{Rot}}$ in the half-plane of $p$;
$[y]_{\operatorname{Rot}\times\operatorname{Scl}}$ is identified with the unit
vector
$\widetilde{y}^{\operatorname{Rot}\times\operatorname{Scl}}=\frac{\widetilde{y}}{\|\widetilde{y}\|}$
projecting $\widetilde{y}$ onto the hemisphere depicted by the vertical
semicircle. Form and shape distances between $[p]$ and $[y]$ correspond to the
length of the geodesics $c(\tau)$ (thick lines) on the plane and sphere,
respectively. Right: Geodesic line $c(\tau)$ between $p=c(0)$ and
$p^{\prime}=c(1)$, Log-map projecting $y$ to $\varepsilon\in
T_{p}\mathcal{M}$, parallel transport $\text{Transp}_{pp^{\prime}}$ forwarding
$\varepsilon$ to $\varepsilon^{\prime}\in T_{p^{\prime}}\mathcal{M}$, and Exp-
map projecting $\varepsilon^{\prime}$ onto $\mathcal{M}$ visualized for a
sphere. Tangent spaces, identified with subspaces of the ambient space, are
depicted as gray planes above the respective poles. The parallel transport
preserves all angles between tangent vectors and identifies
$\dot{c}(0)\cong\dot{c}(1)$.

##  3 Additive Regression on Riemannian Manifolds

Consider a data scenario with $n$ observations of a random response covariate
tuple $\left(Y,\mathbf{X}\right)$, where the realizations of $Y$ are planar
curves $y_{i}:\mathcal{T}\rightarrow\mathds{C}$, $i=1,\dots,n$, belonging to a
Hilbert space $\mathcal{Y}$ defined as above and potentially irregularly
measured on individual grids $t_{i1}<\dots<t_{ik_{i}}\in\mathcal{T}$. The
response object $[Y]$ is the equivalence class of $Y$ with respect to
translation, rotation and possibly scale and the sample
$[y_{1}],\dots,[y_{n}]$ is equipped with the respective Riemannian manifold
geometry introduced in the previous section. For $i=1,\dots,n$, realizations
$\mathbf{x}_{i}\in\mathcal{X}$ of a covariate vector $\mathbf{X}$ in a
covariate space $\mathcal{X}$ are observed. $\mathbf{X}$ can contain several
categorical and/or metric covariates.

For regressing the mean of $[Y]$ on $\mathbf{X}=\mathbf{x}$, we model the
shape/form $[\mu]$ of $\mu\in\mathcal{Y}$ as

$[\mu]=\operatorname{Exp}_{[p]}\left(h(\mathbf{x})\right)=\operatorname{Exp}_{[p]}\left(\sum_{j=1}^{J}h_{j}(\mathbf{x})\right),$
(1)

with an additive predictor $h:\mathcal{X}\rightarrow
T_{[p]}\mathcal{Y}^{*}_{/G}$ acting in the tangent space at an “intercept”
$[p]\in\mathcal{Y}^{*}_{/G}$. Generalizing an additive model
“$Y=\mu+\epsilon=p+h(\mathbf{x})+\epsilon$” in a linear space, we implicitly
define $[\mu]$ as the conditional mean of $[Y]$ given $\mathbf{X}=\mathbf{x}$
by assuming zero-mean “residuals” $\epsilon$. In their definition, we follow
Cornea et al. (2017) but extend to the functional shape/form and additive
case. We assume local linearized residuals
$\varepsilon_{[\mu]}=\operatorname{Log}_{[\mu]}([Y])$ in
$T_{[\mu]}\mathcal{Y}^{*}_{/G}$ to have mean
$\mathds{E}\left(\varepsilon_{[\mu]}\right)=\text{0}$, which corresponds to
$\mathds{E}\left(\varepsilon_{[\mu]}(t)\right)=0$ for ($\nu$-almost) all
$t\in\mathcal{T}$. Here, we assume $[Y]$ is sufficiently close to $[\mu]$ with
probability 1 such that $\operatorname{Log}_{[\mu]}$ is well-defined, which is
the case whenever $\langle\widetilde{Y},\widetilde{\mu}\rangle\neq 0$ for
centered shape/form representatives $\widetilde{Y}$ and $\widetilde{\mu}$, an
un-restrictive and common assumption (compare also Cornea et al., 2017).
However, residuals $\varepsilon_{[\mu]}$ for different $[\mu]$ belong to
separate tangent spaces. To obtain a formulation in a common linear space
instead, local residuals are mapped to residuals
$\epsilon=\operatorname{Transp}_{[\mu],[p]}(\varepsilon_{[\mu]})$ by parallel
transporting them from $[\mu]$ to the common covariate independent pole $[p]$.
After this isometric mapping into $T_{[p]}\mathcal{Y}^{*}_{/G}$, we can
equivalently define the conditional mean $[\mu]$ via
$\mathds{E}\left(\epsilon\right)=\text{0}$ for the transported residuals
$\epsilon$.
$\operatorname{Exp}_{[p]}$ maps the additive predictor
$h(\mathbf{x})=\sum_{j=1}^{J}h_{j}(\mathbf{x})\in T_{[p]}\mathcal{Y}^{*}_{/G}$
to the response space. It is analogous to a response function in GLMs but
depends on $[p]$. While other response functions could be used, we restrict to
the exponential map here, such that the model contains a geodesic model
(Fletcher, 2013) – the direct generalization of simple linear regression – as
a special case for $h(\mathbf{x})=\beta x_{1}$ with a single covariate $x_{1}$
and tangent vector $\beta$. Typically, it is assumed that $h$ is centered such
that $\mathds{E}\left(h(\mathbf{X})\right)=\text{0}$, and the pole $[p]$ is
the overall mean of $[Y]$ defined, like the conditional mean, via residuals of
mean zero.

###  3.1 Tensor-product effect functions $h_{j}$

Scheipl et al. (2015) and other authors employ tensor-product (TP) bases for
functional additive model terms. This naturally extends to tangent space
effects, which we model as

$h_{j}(\mathbf{x})=\sum_{r,l}\theta_{j}^{(r,l)}\,b_{j}^{(l)}(\mathbf{x})\,\partial_{r}$

with the TP basis given by the pair-wise products of $m$ linearly independent
tangent vectors $\partial_{r}\in T_{[p]}\mathcal{Y}^{*}_{/G},\ r=1,\dots,m,$
and $m_{j}$ basis functions $b_{j}^{(l)}:\mathcal{X}\rightarrow\mathds{R},\
l=1,\dots,m_{j},$ for the $j$-th covariate effect depending on one or more
covariates. The real coefficients can be arranged as a matrix
$\\{\theta_{j}^{(r,l)}\\}_{r,l}=\boldsymbol{\Theta}_{j}\in\mathds{R}^{m\times
m_{j}}$. Also for infinite-dimensional $T_{[p]}\mathcal{Y}^{*}_{/G}$ and a
general non-linear dependence on $x$, a basis representation approach requires
truncation to finite dimensions $m$ and $m_{j}$ in practice. Choosing the
bases to capture the essential variability in the data, their size can be
extended with increasing data size and computational resources.

While, in principle, the basis $\\{\partial_{r}\\}_{r}$ could also vary across
effects $j=1,\dots,J$, we assume a common basis for notational simplicity,
which presents the typical choice. Due to the identification of
$T_{[p]}\mathcal{Y}^{*}_{/G}$ with a subspace of the function space
$\mathcal{Y}$, the $\\{\partial_{r}\\}_{r}$ may be specified using a function
basis commonly used in additive models: Let
$b_{0}^{(l)}:\mathcal{T}\rightarrow\mathds{R}$, $l=1,\dots,m_{0}$ be a basis
of real functions, say a B-spline basis (other typical bases used in the
literature include wavelet (Meyer et al., 2015) or FPC bases (Müller and Yao,
2008)). Then we construct the tangent space basis as
$\partial_{r}=\sum_{l=1}^{m_{0}}\left(z_{p}^{(l,r)}+z_{p}^{(m_{0}+l,r)}\sqrt{\text{-}1}\right)b_{0}^{(l)}$,
employing the same basis for the $1$\- and $\sqrt{\text{-}1}$-dimension before
transforming it with a basis transformation matrix
$\mathbf{Z}_{p}=\\{z_{p}^{(l,r)}\\}_{l,r}\in\mathds{R}^{2m_{0}\times m}$ with
$m<2m_{0}$ implementing the linear tangent space constraints (Section 2).
Practically, $\mathbf{Z}_{p}$ is obtained as null space basis matrix of the
matrix
$(\operatorname{Re}\\!\left(\mathbf{C}\right),\operatorname{Im}\\!\left(\mathbf{C}\right))$
with $\mathbf{C}=\\{\langle b_{0}^{(l)},\zeta^{(r)}\rangle\\}_{r,l}$ (or with
the empirical inner product on the product space of irregular curves instead)
constructed from the normal vectors $\zeta^{(r)}\in\mathcal{Y}$,
$r=1,\dots,2m_{0}-m,$ to $T_{[p]}\mathcal{Y}^{*}_{/G}$. For closed curves, we
additionally choose $\mathbf{Z}_{p}$ to enforce periodicity, i.e.
$\partial_{r}(t)=\partial_{r}(t+t_{0})$ for some $t_{0}\in\mathds{R}$ (compare
Hofner et al., 2016).

Given the tangent space basis, we may now modularly specify the usual additive
model basis functions $b_{j}^{(l)}:\mathcal{X}\rightarrow\mathds{R}$,
$l=1,\dots,m_{j}$, for the $j$-th covariate effect to obtain the full
functional additive model “tool box” offered by, e.g., Brockhaus et al.
(2015). A linear effect – linear in the tangent space – of the form
$h_{j}(\mathbf{x})=\beta z$ with a scalar (typically centered) covariate $z$
in $\mathbf{x}$ and $\beta\in T_{[p]}\mathcal{Y}^{*}_{/G}$ is simply
implemented by a single function $b_{j}^{(1)}(\mathbf{x})=z$. A smooth effect
of the generic form $h_{j}(\mathbf{x})(t)=f(z,t)$ can be implemented by
choosing, e.g., a B-spline basis (Asymptotic properties of penalized B-splines
and connections to kernel estimators are discussed, e.g., by Wood et al.
(2016); Li and Ruppert (2008)). For a categorical covariate effect of the form
$h_{j}(\mathbf{x}):\\{1,\dots,K\\}\rightarrow T_{[p]}\mathcal{Y}^{*}_{/G}$,
$\kappa\mapsto\beta_{\kappa}$, the basis
$\mathbf{b}_{j}(\mathbf{x}):\kappa\mapsto\mathbf{e}_{\kappa}\in\mathds{R}^{K-1}$
maps category $\kappa$ to a usual contrast vector $\mathbf{e}_{\kappa}$ just
as in standard linear models. Here, we typically use effect-encoding to obtain
centered effects. Moreover, TP interactions of the model terms described
above, as well as group-specific effects and smooth effects with additional
constraints (Hofner et al., 2016) can be specified in the model formula,
relying on the mboost framework introduced by Hothorn et al. (2010), which
also allows to define custom effect designs. For identification of an overall
mean intercept $[p]$, sum-to-zero constraints yielding
$\sum_{i=1}^{n}h_{j}(\mathbf{x}_{i})=\text{0}$ for observed covariates
$\mathbf{x}_{i}$ can be specified, and similar constraints can be used to
distinguish linear from non-linear effects and interactions from their
marginal effects (Kneib et al., 2009). Different quadratic penalties can be
specified for the coefficients $\boldsymbol{\Theta}_{j}$, allowing to
regularize high-dimensional effect bases and to balance effects of different
complexity in the model fit (cf. Section 4).

###  3.2 Tensor-product factorization

The multidimensional structure of the response objects makes it challenging to
graphically illustrate and interpret additive model terms, in particular when
it comes to non-linear (interaction) effects, or when effect sizes are
visually small. To solve this problem, we suggest to re-write estimated TP
effects ${\hat{h}}_{j}$ with estimated coefficient matrix
${\widehat{\boldsymbol{\Theta}}}_{j}$ as

${\hat{h}}_{j}(\mathbf{x})=\sum_{r=1}^{m_{j}^{\prime}}\xi^{(r)}_{j}{\hat{h}}_{j}^{(r)}(\mathbf{x})$

factorized into $m_{j}^{\prime}=\min(m_{j},m_{0})$ components consisting of
covariate effects ${\hat{h}}_{j}^{(r)}:\mathcal{X}\rightarrow\mathds{R},\
r=1,\dots,m_{j}^{\prime},$ in corresponding orthonormal directions
$\xi^{(r)}_{j}\in T_{[p]}\mathcal{Y}^{*}_{/G}$ with
$\langle\xi_{j}^{(r)},\xi_{j}^{(l)}\rangle=\mathds{1}({r=l})$, i.e. $1$ if
$r=l$ and 0 otherwise. Assuming
$\mathds{E}\left(b_{j}^{(l)}(\mathbf{X})^{2}\right)<\infty$,
$l=1,\dots,m_{j}$, for the underlying effect basis, the ${\hat{h}}_{j}^{(r)}$
are specified to achieve decreasing component variances
$v_{j}^{(1)}\geq\dots\geq v_{j}^{(m_{j}^{\prime})}\geq 0$ given by
$v_{j}^{(r)}=\mathds{E}\left({\hat{h}}_{j}^{(r)}(\mathbf{X})^{2}\right)$. In
practice, the expectation over the covariates $\mathbf{X}$ and the inner
product $\langle.,.\rangle$ are replaced by empirical analogs (compare
Supplement Corollary 3). Due to orthonormality of the $\xi_{j}^{(r)}$, the
component variances add up to the total predictor variance
$\sum_{r=1}^{m_{j}^{\prime}}v_{j}^{(r)}=v_{j}=\mathds{E}\left(\langle{\hat{h}}_{j}(\mathbf{X}),{\hat{h}}_{j}(\mathbf{X})\rangle\right)$.
Moreover, the TP factorization is optimally concentrated in the first
components in the sense that for any $l\leq m_{j}^{\prime}$ there is no
sequence of $\xi^{(r)}_{*}\in\mathcal{Y}$ and
${\hat{h}}_{*}^{(r)}:\mathcal{X}\rightarrow\mathds{R}$, such that
$\mathds{E}\left(\|{\hat{h}}_{j}(\mathbf{X})-\sum_{r=1}^{l}\xi^{(r)}_{*}{\hat{h}}_{*}^{(r)}(\mathbf{X})\|^{2}\right)<\mathds{E}\left(\|h_{j}(\mathbf{X})-\sum_{r=1}^{l}\xi_{j}^{(r)}{\hat{h}}^{(r)}_{j}(\mathbf{X})\|^{2}\right)$,
i.e. the series of the first $l$ components yields the best rank $l$
approximation of ${\hat{h}}_{j}$. The factorization relies on SVD of (a
transformed version of) the coefficient matrix
${\widehat{\boldsymbol{\Theta}}}_{j}$ and the fact that it is well-defined is
a variant of the Eckart-Young-Mirsky theorem (proof in Supplement S.2).
Particularly when large shares of the predictor variance are explained by the
first component(s), the decomposition facilitates graphical illustration and
interpretation: choosing a suitable constant $\tau\neq 0$, an effect direction
$\xi^{(r)}_{j}$ can be visualized by plotting the pole representative $p$
together with $\operatorname{Exp}_{p}(\tau\,\xi^{(r)}_{j})$ on the level of
curves, while accordingly re-scaled
$\frac{1}{\tau}{\hat{h}}_{j}^{(r)}(\mathbf{x})$ is displayed separately in a
standard scalar effect plot. Adjusting $\tau$ offers an important degree of
freedom for visualizing $\xi^{(r)}_{j}$ on an intuitively accessible scale
while faithfully depicting $\xi^{(r)}_{j}{\hat{h}}_{j}^{(r)}(\mathbf{x})$.
When based on the same $\tau$, different covariate effects can be compared
across the plots sharing the same scale. We suggest
$\tau=\max_{j}\sqrt{v_{j}}$, the maximum total predictor standard deviation of
an effect, as a good first choice.

Besides factorizing effects separately, it can also be helpful to apply TP
factorization to the joint additive predictor, yielding

$h(\mathbf{x})=\sum_{r=1}^{m^{\prime}}\xi^{(r)}{\hat{h}}^{(r)}(\mathbf{x})=\sum_{r=1}^{m^{\prime}}\xi^{(r)}\left({\hat{h}}_{1}^{(r)}(\mathbf{x})+\dots+{\hat{h}}_{J}^{(r)}(\mathbf{x})\right),\quad
m^{\prime}=\min(\sum_{j}m_{j},m),$

with again $\xi^{(r)}\in T_{[p]}\mathcal{Y}^{*}_{/G}$ orthonormal and the
corresponding variance concentration in the first components, but now
determined w.r.t. entire additive predictors
${\hat{h}}^{(r)}=\sum_{j=1}^{J}{\hat{h}}_{j}^{(r)}$ spanned by all covariate
basis functions in the predictor. In this representation, the first component
yields a geodesic additive model approximation where the predictor moves along
a geodesic line $c(\tau)=\operatorname{Exp}_{[p]}\left(\xi^{(1)}\tau\right)$
with the signed distance $\tau\in\mathds{R}$ from $[p]$, modeled by a scalar
additive predictor ${\hat{h}}^{(1)}(\mathbf{x})$ composed of covariate effects
analogous to the original model predictor. In Section 5, we illustrate its
potential in three different scenarios.

##  4 Component-wise Riemannian $L_{2}$-Boosting

Component-wise gradient boosting (e.g. Hothorn et al., 2010) is a step-wise
model fitting procedure accumulating predictors from smaller models, so called
base-learners, to built an ensemble predictor aiming at minimizing a mean loss
function. To this end, the base-learners are fit (via least squares) to the
negative gradient of the loss function in each step and the best fitting base-
learner is added to the current ensemble predictor. Due to its versatile
applicability, inherent model selection, and slow over-fitting behavior,
boosting has proven useful in various contexts (Mayr et al., 2014). Boosting
with respect to the least squares loss function
$\ell(y,\mu)=\frac{1}{2}(y-\mu)^{2}$, $y,\mu\in\mathds{R}$, is typically
referred to as $L_{2}$-Boosting and simplifies to repeated re-fitting of
residuals $\varepsilon=y-\mu=-\nabla_{\mu}\ell(y,\mu)$ corresponding to the
negative gradient of the loss function. For $L_{2}$-Boosting with a single
learner, Bühlmann and Yu (2003) show how fast bias decay and slow variance
increase over the boosting iterations suggest stopping the algorithm early
before approaching the ordinary (penalized) least squares estimator. Lutz and
Bühlmann (2006) prove consistency of component-wise $L^{2}$-Boosting in a
high-dimensional multivariate response linear regression setting and Stöcker
et al. (2021) illustrate in extensive simulation studies how stopping the
boosting algorithm early based on curve-wise cross-validation applies desired
regularization when fitting (even highly autocorrelated) functional responses
with parameter-intense additive model base-learners and, thus, leads to good
estimates even in challenging scenarios.
When generalizing to least squares on Riemannian manifolds with the loss
$\frac{1}{2}d^{2}([y],[\mu])$ given by the squared geodesic distance, the
negative gradient
$-\nabla_{[\mu]}\frac{1}{2}d^{2}([y],[\mu])=\operatorname{Log}_{[\mu]}([y])=\varepsilon_{[\mu]}$
(compare e.g. Pennec, 2006) corresponds to the local residuals
$\varepsilon_{[\mu]}$ defined in Section 3. This analogy to $L_{2}$-Boosting
motivates the presented generalization where local residuals are further
transported to residuals $\epsilon$ in a common linear space.

Consider the pole $[p]$ known and fixed for now. Assuming its existence, we
aim to minimize the population mean loss

$\sigma^{2}(h)=\mathds{E}\left(d^{2}\left([Y],\operatorname{Exp}_{[p]}\left(h(\mathbf{X})\right)\right)\right)$

with the point-wise minimizer
$h^{\star}(\mathbf{x})=\underset{h:\mathcal{X}\rightarrow
T_{[p]}\mathcal{Y}_{/G}^{*}}{\operatorname{argmin}\;}\mathds{E}\left(d^{2}\left([Y],\operatorname{Exp}_{[p]}\left(h(\mathbf{X})\right)\right)\mid\mathbf{X}=\mathbf{x}\right)$
minimizing the conditional expected squared distance. Fixing a covariate
constellation $\mathbf{x}\in\mathcal{X}$, the prediction
$[\mu]=\operatorname{Exp}_{[p]}\left(h^{\star}(\mathbf{x})\right)$ corresponds
to the Fréchet mean (Karcher, 1977) of $[Y]$ conditional on
$\mathbf{X}=\mathbf{x}$. In a finite-dimensional context, Pennec (2006) show
that $\mathds{E}\left(\varepsilon_{[\mu]}\right)=\text{0}$ for a Fréchet mean
$[\mu]$ if residuals $\varepsilon_{[\mu]}$ are uniquely defined with
probability one. This indicates the connection to our residual based model
formulation in Section 3. We fit the model by reducing the empirical mean loss
$\hat{\sigma}^{2}(h)=\frac{1}{n}\sum_{i=1}^{n}d_{i}^{2}\left([y_{i}],\operatorname{Exp}_{[p]}\left(h(\mathbf{x}_{i})\right)\right),$
where we replace the population mean by the sample mean and compute the
geodesic distances $d_{i}$ with respect to the inner products
$\langle\cdot,\cdot\rangle_{i}$ defined for the respective evaluations of
$y_{i}$.

A base-learner corresponds to a covariate effect
$h_{j}(\mathbf{x})=\sum_{r,l}\theta_{j}^{(r,l)}\,b_{j}^{(l)}(\mathbf{x})\,\partial_{r}$,
$\boldsymbol{\Theta}_{j}=\\{\theta_{j}^{(r,l)}\\}_{r,l}$, which is repeatedly
fit to the transported residuals $\epsilon_{1},\dots,\epsilon_{n}$ by
penalized least-squares (PLS) minimizing
$\sum_{i=1}^{n}\|\epsilon_{i}-h_{j}(\mathbf{x}_{i})\|_{i}^{2}+\lambda_{j}\operatorname{tr}{\left(\boldsymbol{\Theta}_{j}\mathbf{P}_{j}\boldsymbol{\Theta}_{j}^{\top}\right)}+\lambda\operatorname{tr}{\left(\boldsymbol{\Theta}^{\top}\mathbf{P}\boldsymbol{\Theta}\right)}$.
Via the penalty parameters $\lambda_{j},\lambda\geq 0$ the effective degrees
of freedom of the base-learners are controlled (Hofner et al., 2011) to
achieve a balanced “fair” base-learner selection despite the typically large
and varying number of coefficients involved in the TP effects. The symmetric
penalty matrices $\mathbf{P}_{j}\in\mathds{R}^{m_{j}\times m_{j}}$ and
$\mathbf{P}\in\mathds{R}^{m\times m}$ (imposing, e.g., a second-order
difference penalty for B-splines in either direction) can equivalently be
arranged as a $m_{j}m\times m_{j}m$ penalty matrix
$\mathbf{R}_{j}=\lambda_{j}(\mathbf{P}_{j}\otimes\mathbf{I}_{m})+\lambda(\mathbf{I}_{m_{j}}\otimes\mathbf{P})$
for the vectorized coefficients
$\operatorname{vec}{(\boldsymbol{\Theta}_{j})}=(\theta^{(1,1)}_{j},\dots,\theta^{(m,1)}_{j},\dots,\theta^{(m,m_{j})})^{\top}$,
where $\otimes$ denotes the Kronecker product. The standard PLS estimator is
then given by
$\operatorname{vec}{(\widehat{\boldsymbol{\Theta}}_{j})}=\left(\boldsymbol{\Psi}_{j}+\mathbf{R}_{j}\right)^{-1}\boldsymbol{\psi}_{j}$
with
$\boldsymbol{\Psi}_{j}=\sum_{i=1}^{n}\left\\{\operatorname{Re}\\!\left(\langle
b^{(l)}_{j}(\mathbf{x}_{i})\partial_{r},b^{(l^{\prime})}_{j}(\mathbf{x}_{i})\partial_{r^{\prime}}\rangle_{i}\right)\right\\}_{\begin{subarray}{c}(r,l)=(1,1),\dots,(m,1),\dots,(m,m_{j})\\\
(r^{\prime},l^{\prime})=(1,1),\dots,(m,1),\dots,(m,m_{j})\end{subarray}}\in\mathds{R}^{m\,m_{j}\times
m\,m_{j}}$ and
$\boldsymbol{\psi}_{j}=\sum_{i=1}^{n}\left\\{\operatorname{Re}\\!\left(\langle
b^{(l)}_{j}(\mathbf{x}_{i})\partial_{r},\epsilon_{i}\rangle_{i}\right)\right\\}_{(r,l)=(1,1),\dots,(m,1),\dots,(m,m_{j})}\in\mathds{R}^{m\,m_{j}}$.
In a regular design, using the functional linear array model (Brockhaus et
al., 2015) can save memory and computation time by avoiding construction of
the complete matrices. The basis construction of $\\{\partial_{r}\\}_{r}$ via
a transformation matrix $\mathbf{Z}_{p}$ (Section 3.1) is reflected in the
penalty by setting
$\mathbf{P}=\mathbf{Z}_{p}^{\top}(\mathbf{I}_{2}\otimes\mathbf{P}_{0})\mathbf{Z}_{p}$
with $\mathbf{P}_{0}$ the penalty matrix for the un-transformed basis
$\\{b_{0}^{(r)}\\}_{r}$.

In each iteration of the proposed Algorithm 1, the best-performing base-
learner is added to the current ensemble additive predictor $h(\mathbf{x})$
after multiplying it with a step-length parameter $\eta\in(0,1]$. Due to the
additive model structure this corresponds to a coefficient update of the
selected covariate effect. Accordingly, after repeated selection, the
effective degrees of freedom of a covariate effect, in general, exceed the
degrees specified for the base-learner. They are successively adjusted to the
data. To avoid over-fitting, the algorithm is typically stopped early before
reaching a minimum of the empirical mean loss. The stopping iteration is
determined, e.g., by re-sampling strategies such as bootstrapping or cross-
validation on the level of shapes/forms.

# Initialization:

Geometry : specify geometry (shape/form) and pole representative $p$

Hyper-parameters : Step-length $\eta\in(0,1]$, number of boosting iterations

Base-learners : $h_{j}(\mathbf{x})$ with penalty matrix $\mathbf{R}_{j}$ and

initial coefficient matrix $\boldsymbol{\Theta}_{j}=\boldsymbol{0}$

for  _$j=1$ to $J$_ do # Prepare penalized least-squares (PLS)

# set up $m\,m_{j}\times m\,m_{j}$ matrix:
$\boldsymbol{\Psi}_{j}\leftarrow\sum_{i=1}^{n}\left\\{\operatorname{Re}\\!\left(\langle
b^{(l)}_{j}(\mathbf{x}_{i})\partial_{r},b^{(l^{\prime})}_{j}(\mathbf{x}_{i})\partial_{r^{\prime}}\rangle_{i}\right)\right\\}_{\begin{subarray}{c}(r,l)=(1,1),\dots,(m,1),\dots,(m,m_{j})\\\
(r^{\prime},l^{\prime})=(1,1),\dots,(m,1),\dots,(m,m_{j})\end{subarray}}$

end

repeat # boosting steps

for  _$i=1,\dots,n$_ do # Compute current transported residuals

$[\mu_{i}]\leftarrow\operatorname{Exp}_{[p]}(h(\mathbf{x}_{i}))$

$\varepsilon_{[\mu_{i}]}\leftarrow\operatorname{Log}_{[\mu_{i}]}([y_{i}])$

$\epsilon_{i}\leftarrow\operatorname{Transp}_{[\mu_{i}],[p]}(\varepsilon_{[\mu_{i}]})$

end

for  _$j=1,\dots,J$_ do # PLS fit to residuals

# $m\,m_{j}$ vector:
$\boldsymbol{\psi}_{j}\leftarrow\sum_{i=1}^{n}\left\\{\operatorname{Re}\\!\left(\langle
b^{(l)}_{j}(\mathbf{x}_{i})\partial_{r},\epsilon_{i}\rangle_{i}\right)\right\\}_{(r,l)=(1,1),\dots,(m,1),\dots,(m,m_{j})}$

$\widehat{\boldsymbol{\Theta}}_{j}=\\{\hat{\theta}_{j}^{(r,l)}\\}_{r,l}\leftarrow$
Solve(_
$\left(\boldsymbol{\Psi}_{j}+\mathbf{R}_{j}\right)\operatorname{vec}(\boldsymbol{\Theta})=\boldsymbol{\psi}_{j}$
_)

end

$\hat{\jmath}\leftarrow\underset{j\in\\{1,\dots,J\\}}{\operatorname{argmin}\;}\sum_{i=1}^{n}\|\epsilon_{i}-\sum_{r,l}\hat{\theta}^{(r,l)}_{j}b_{j}^{(l)}(\mathbf{x})\partial_{r}\|_{i}^{2}$;
# Select base-learner

$\boldsymbol{\Theta}_{\hat{\jmath}}\leftarrow\boldsymbol{\Theta}_{\hat{\jmath}}+\eta\,\widehat{\boldsymbol{\Theta}}_{\hat{\jmath}}$;
# Update selected model coefficients

until  _Stopping criterion (e.g. minimal cross-validation error)_

Algorithm 1 Component-wise Riemannian $L^{2}$-Boosting

The pole $[p]$ is, in fact, usually not a priori available. Instead we
typically assume
$[p]=\underset{q\in\mathcal{Y}^{*}}{\operatorname{argmin}\;}\mathds{E}\left(d^{2}([Y],[q])\right)$
is the overall Fréchet mean, also often referred to as Riemannian center of
mass for Riemannian manifolds or as Procrustes mean in shape analysis (Dryden
and Mardia, 2016). Here, we estimate it as
$[p]=\operatorname{Exp}_{[p_{0}]}(h_{0})$ in a preceding Riemannian
$L^{2}$-Boosting routine. The constant effect $h_{0}\in
T_{[p_{0}]}\mathcal{Y}^{*}_{/G}$ in the intercept-only special case of our
model is estimated with Algorithm 1 based on a preliminary pole
$[p_{0}]\in\mathcal{Y}^{*}_{/G}$. For shapes and forms, a good candidate for
$p_{0}$ can be obtained as the standard functional mean of a reasonably well
aligned sample $y_{1},\dots,y_{n}\in\mathcal{Y}$ of representatives.

The proposed Riemannian $L_{2}$-Boosting algorithm is available in the R (R
Core Team, 2018) package manifoldboost (github.com/Almond-S/manifoldboost).
The implementation is based on the package FDboost (Brockhaus et al., 2020),
which is in turn based on the model-based boosting package mboost (Hothorn et
al., 2010).

##  5 Applications and Simulation

###  5.1 Shape differences in astragali of wild and domesticated sheep

In a geometric morphometric study, Pöllath et al. (2019) investigate shapes of
sheep astragali (ankle bones) to understand the influence of different living
conditions on the micromorphology of the skeleton. Based on a total of $n=163$
shapes recorded by Pöllath et al. (2019), we model the astragalus shape in
dependence on different variables, including domestication status
(wild/feral/domesticated), sex (female/male/NA), age
(juvenile/subadult/adult/NA), and mobility (confined/pastured/free) of the
animals as categorical covariates. The sample comprises sheep of four
different populations: Asiatic wild sheep (Field Museum, Chicago; Lay, 1967;
Zeder, 2006), feral Soay sheep (British Natural History Museum, London;
Clutton-Brock et al., 1990), and domestic sheep of the Karakul and Marsch
breed (Museum of Livestock Sciences, Halle (Saale); Schafberg and Wussow,
2010). Table S1 in Supplement S.3 shows the distribution of available
covariates within the populations. Each sheep astragalus shape, $i=1,\dots,n$,
is represented by a configuration composed of 11 selected landmarks in a
vector $\mathbf{y}_{i}^{\text{lm}}\in\mathds{C}^{11}$ and two vectors of
sliding semi-landmarks $\mathbf{y}_{i}^{\text{c1}}\in\mathds{C}^{14}$ and
$\mathbf{y}_{i}^{\text{c2}}\in\mathds{C}^{18}$ evaluated along two outline
curve segments, marked on a 2D image of the bone (dorsal view). Several
example configurations are displayed in Supplement Figure S1. In general, we
could separately specify smooth function bases for the outline segments
$y_{i}^{\text{c1}}$ and $y_{i}^{\text{c2}}$, respectively. Due to their
systematic recording, we assume, however, that not only landmarks but also
semi-landmarks are regularly observed on a fixed grid, and refrain from using
smooth function bases for simplicity. Accordingly, shape configurations can
directly be identified with their evaluation vectors
$\mathbf{y}_{i}=\big{(}\mathbf{y}_{i}^{\text{lm}\top},\mathbf{y}_{i}^{\text{c1}\top},\mathbf{y}_{i}^{\text{c2}\top}\big{)}^{\top}\in\mathds{C}^{43}=\mathcal{Y}$,
and the geometry of the response space
$\mathcal{Y}^{*}_{/\operatorname{Trl}\times\operatorname{Rot}\times\operatorname{Scl}}$
widely corresponds to the classic Kendall’s shape space geometry, with the
difference that, considering landmarks more descriptive than single semi-
landmarks, we choose a weighted inner product
$\langle\mathbf{y}_{i},\mathbf{y}_{i}^{\prime}\rangle=\mathbf{y}_{i}^{\dagger}\mathbf{W}\mathbf{y}_{i}^{\prime}$
with diagonal weight matrix $\mathbf{W}$ with diagonal
$\big{(}\boldsymbol{1}_{11}^{\top},\frac{3}{14}\boldsymbol{1}_{14}^{\top},\frac{3}{18}\boldsymbol{1}_{18}^{\top}\big{)}^{\top}$
assigning the weight of three landmarks to each outline segment. We model the
astragalus shapes
$[\mathbf{y}_{i}]\in\mathcal{Y}^{*}_{/\operatorname{Trl}\times\operatorname{Rot}\times\operatorname{Scl}}$
as

$\displaystyle[\boldsymbol{\mu}_{i}]$
$\displaystyle=\operatorname{Exp}_{[\mathbf{p}]}\left(\boldsymbol{\beta}_{\mathrm{status}_{i}}+\boldsymbol{\beta}_{\mathrm{{pop}}_{i}}+\boldsymbol{\beta}_{\mathrm{age}_{i}}+\boldsymbol{\beta}_{\mathrm{sex}_{i}}+\boldsymbol{\beta}_{\mathrm{mobility}_{i}}\right)$

with the pole $[\mathbf{p}]\in\mathcal{Y}^{*}_{/G}$ specified as overall mean
and the conditional mean
$[\boldsymbol{\mu}_{i}]\in\mathcal{Y}^{*}_{/\operatorname{Trl}\times\operatorname{Rot}\times\operatorname{Scl}}$
depending on the effect coded covariate effects
$x_{ij}\mapsto\boldsymbol{\beta}_{x_{ij}}\in
T_{[\mathbf{p}]}\mathcal{Y}^{*}_{/\operatorname{Trl}\times\operatorname{Rot}\times\operatorname{Scl}}$.
For identifiability, the population and mobility effects are centered around
the status effect, as we only have data on different populations/mobility
levels for domesticated sheep. All base-learners are regularized to one degree
of freedom by employing ridge penalties for the coefficients of the covariate
bases $\\{b_{j}^{(l)}\\}_{l}$ while the coefficients of the response basis
(the standard basis for $\mathds{C}^{43})$ are left un-penalized. With a step-
length of $\eta=0.1$, 10-fold shape-wise cross-validation suggests early
stopping after $89$ boosting iterations. Due to the regular design, we can
make use of the functional linear array model (Brockhaus et al., 2015) for
saving computation time and memory, which lead to 8 seconds of initial model
fit followed by 47 seconds of cross-validation. To interpret the categorical
covariate effects, we rely on TP factorization (Figure 2). The first component
of the status effect explains about 2/3 of the variance of the status effect
and over 50% of the cumulative effect variance in the model. In that main
direction, the effect of feral is not located between wild and domestic, as
might be naively expected. By contrast, the second component of the effect
seems to reflect the expected order and still explains a considerable amount
of variance. Similar to Pöllath et al. (2019), we find little influence of
age, sex and mobility on the astragalus shape. Yet, all covariates were
selected by the boosting algorithm.

Figure 2:  Left: Shares of different factorized covariate effects in the total
predictor variance. Right: Factorized effect plots showing the two components
of the status effect (rows): in the right column, the two first directions
$\boldsymbol{\xi}^{(1)}_{1},\boldsymbol{\xi}^{(2)}_{1}\in
T_{[\mathbf{p}]}\mathcal{Y}^{*}_{/\operatorname{Trl}+\operatorname{Rot}+\operatorname{Scl}}$
are visualized via line-segments originating at the overall mean shape (empty
circles) and ending in the shape resulting from moving 1 unit into the target
direction (solid circles; large: landmarks; small: semi-landmarks along the
outline); in the left column, the status effect in the respective direction is
depicted. As illustrated in the middle plot, an effect of 1 would correspond
to the full extend of the direction shown to the right.

Visually, differences in estimated mean shapes are rather small, which is, in
our experience, quite usual for shape data. With differences in size, rotation
and translation excluded by definition, only comparably small variance remains
in the observed shapes. Nonetheless, TP factorization provides accessible
visualization of the effect directions and allows to partially order the
effect levels in each direction.

###  5.2 Cellular Potts model parameter effects on cell form

The stochastic biophysical model proposed by Thüroff et al. (2019), a cellular
Potts model (CPM), simulates migration dynamics of cells (e.g. wound healing
or metastasis) in two dimensions. The progression of simulated cells is the
result of many consecutive local elementary events sampled with a Metropolis-
algorithm according to a Hamiltonian. Different parameters controlling the
Hamiltonian have to be calibrated to match real live cell properties
(Schaffer, 2021). Considering whole cells, parameter implications on the cell
form are not obvious. To provide additional insights, we model the cell form
in dependence on four CPM parameters considered particularly relevant: the
bulk stiffness $x_{i1}$, membrane stiffness $x_{i2}$, substrate adhesion
$x_{i3}$, and signaling radius $x_{i4}$ are subsumed in a vector
$\mathbf{x}_{i}$ of metric covariates for $i=1,\dots,n$. Corresponding sampled
cell outlines $y_{i}$ were provided by Sophia Schaffer in the context of
Schaffer (2021), who ran underlying CPM simulations and extracted outlines.
Deriving the intrinsic orientation of the cells from their movement
trajectories, we parameterize $y_{i}:[0,1]\rightarrow\mathds{C}$, clockwisely
relative to arc-length such that $y_{i}(0)=y_{i}(1)$ points into the movement
direction of the barycenter of the cell. With an average of
$k=\frac{1}{n}\sum_{i=1}^{n}k_{i}\approx 43$ samples per curve (after sub-
sampling preserving $95\%$ of their inherent variation, as described in
Volkmann et al., 2021, Supplement), the evaluation vectors
$\mathbf{y}_{i}\in\mathds{C}^{k_{i}}$ are equipped with an inner-product
implementing trapezoidal rule integration weights. Example cell outlines are
depicted in Supplement Figure S4. The results shown below are based on cell
samples obtained from 30 different CPM parameter configurations. For each
configuration, 33 out of 10.000 Monte-Carlo samples were extracted as
approximately independent. This yields a dataset of $n=990=30\times 33$ cell
outlines.

As positioning of the irregularly sampled cell outlines $y_{i}$,
$i=1,\dots,n$, in the coordinate system is arbitrary, we model the cell forms
$[y_{i}]\in\mathcal{Y}^{*}_{/\operatorname{Trl}+\operatorname{Rot}}$. Their
estimated overall form mean $[p]$ serves as pole in the additive model

$\displaystyle[\mu_{i}]=\operatorname{Exp}_{[p]}\big{(}h(\mathbf{x}_{i})\big{)}=\operatorname{Exp}_{[p]}\big{(}\sum_{j}\beta_{j}x_{ij}+\sum_{j}f_{j}(x_{ij})+\sum_{j\neq{\ddot{\jmath}}}f_{j{\ddot{\jmath}}}(x_{ij},x_{i{\ddot{\jmath}}})\big{)}$

where the conditional form mean $[\mu_{i}]$ is modeled in dependence on
tangent-space linear effects with coefficients $\beta_{j}\in
T_{[p]}\mathcal{Y}_{/\operatorname{Trl}\times\operatorname{Rot}}$ and non-
linear smooth effects $f_{j}$ for covariate $j=1,\dots,4$, as well as smooth
interaction effects $f_{j{\ddot{\jmath}}}$ for each pair of covariates
$j\neq{\ddot{\jmath}}$. All involved (effect) functions are modeled via a
cyclic cubic P-spline basis $\\{b_{0}^{(r)}\\}_{r}$ with $7$ (inner) knots and
a ridge penalty, and quadratic P-splines with $4$ knots for the covariates
$x_{ij}$ equipped with a second order difference penalty for the $f_{j}$ and
ridge penalties for interactions. Covariate effects are mean centered and
interaction effects $f_{j{\ddot{\jmath}}}(x_{j},x_{\ddot{\jmath}})$ are
centered around their marginal effects
$f_{j}(x_{j}),f_{\ddot{\jmath}}(x_{\ddot{\jmath}})$, which are in turn
centered around the linear effects $\beta_{j}x_{j}$ and
$\beta_{\ddot{\jmath}}x_{\ddot{\jmath}}$, respectively. Resulting predictor
terms involve 69 (linear effect) to 1173 (interaction) basis coefficients but
are penalized to a common degree of freedom of 2 to ensure a fair base-learner
selection. We fit the model with a step-size of $\eta=0.25$ and stop after
2000 boosting iterations observing no further meaningful risk reduction, since
no need for early-stopping is indicated by 10-fold form-wise cross-validation.
Due to the increased number of data points and coefficients, the irregular
design, and the increased number of iterations, the model fit takes
considerably longer than in Section 5.1, with about 50 initial minutes
followed by 8 hours of cross-validation. However, as usual in boosting, model
updates are large in the beginning and only marginal in later iterations, such
that fits after $1000$ or $500$ iterations would already yield very similar
results.

Observing that the most relevant components point into similar directions, we
jointly factorize the predictor as
${\hat{h}}(\mathbf{x}_{i})=\sum_{r}\xi^{(r)}{\hat{h}}^{(r)}(\mathbf{x}_{i})$
with TP factorization. The first component explains about 93% of the total
predictor variance (Supplement Fig. S3), indicating that, post-hoc, a good
share of the model can be reduced to the geodesic model
$[{\hat{\mu}}_{i}]=\operatorname{Exp}_{[p]}(\xi^{(1)}{\hat{h}}^{(1)}(\mathbf{x}_{i}))$
illustrated in Figure 3. A positive effect in the direction $\xi^{(1)}$ makes
cells larger and more keratocyte / croissant shaped, a negative effect –
pointing into the opposite direction – makes them smaller and more mesenchymal
shaped / elongated.

Figure 3: Center: the main direction $\xi^{(1)}$ of the model illustrated as
vectors pointing from the overall mean cell form $[p]$ (grey curve) to the
form $\operatorname{Exp}_{[p]}(\xi^{(1)})$ (blue dots), which are both
oriented as cells migrating rightwards. Left: Effects of the bulk stiffness
$x_{i1}$ into the direction $\xi^{(1)}$. A vertical line from 0, corresponding
to $[p]$, to 1, corresponding to the full extent of $\xi^{(1)}$, underlines
the connection between the plots and helps to visually asses the amount of
change for a given value of $x_{i1}$. Right: The overall effect of $x_{i1}$
and membrane stiffness $x_{i2}$, comprising linear, smooth and interaction
effects, as a 3D surface plot. The heat map plotted on the surface shows only
the interaction effect $f^{(1)}_{12}(x_{i1},x_{i2})$ illustrating deviations
from the marginal effects, which are of particular interest for CPM
calibration.

The bulk stiffness $x_{i1}$ turns out to present the most important driving
factor behind the cell form, explaining over 75% of the cumulative variance of
the effects (Supplement Fig. S2). Around 80% of its effect are explained by
the linear term reflecting gradual shrinkage at the side of the cells with
increasing bulk stiffness.

###  5.3 Realistic shape and form simulation studies

To evaluate the proposed approach, we conduct simulation studies for both form
and shape regression for irregular curves. We compare sample sizes
$n\in\\{54,162\\}$ and average grid sizes
$k=\frac{1}{n}\sum_{i=1}^{n}k_{i}\in\\{40,100\\}$ as well as an extreme case
with $k_{i}=3$ for each curve but $n=720$, i.e. where only random triangles
are observed (yet, with known parameterization over $[0,1]$). We additionally
investigate the influence of nuisance effects and compare different inner
product weights. While important results are summarized in the following,
comprehensive visualizations can be found in Supplement S.5.

Simulation design: We simulate models of the form
$[\mu]=\operatorname{Exp}_{[p]}\left(\beta_{\kappa}+f_{1}(z_{1})\right)$ with
overall mean $[p]$, a binary effect with levels $\kappa\in\\{0,1\\}$ and a
smooth effect of $z_{1}\in[-60,60]$. We choose a cyclic cubic B-spline basis
with $27$ knots for $T_{[p]}\mathcal{Y}^{*}_{/G}$, placing them irregularly at
1/27-quantiles of unit-speed parameterization time-points of the curves. Cubic
B-splines with $4$ regularly placed knots are used for covariates in smooth
effects. True models are based on the bot dataset from R package Momocs
(Bonhomme et al., 2014) comprising outlines of 20 beer ($\kappa=0$) and 20
whiskey ($\kappa=1$) bottles of different brands. A smooth effect is induced
by the 2D viewing transformations resulting from tilting the planar outlines
in a 3D coordinate system along their longitudinal axis by an angle of up to
60 degree towards the viewer ($z_{1}=60$) and away ($z_{1}=-60$) (i.e. in a
way not captured by 2D rotation invariance). Establishing ground truth models
based on a fit to the bottle data, we simulate new responses
$[y_{1}],\dots,[y_{n}]$ via residual re-sampling (Supplement S.5) to preserve
realistic autocorrelation. Subsequently, we randomly translate, rotate and
scale $y_{1},\dots,y_{n}\in\mathcal{Y}$ somewhat around the aligned form/shape
representatives to obtain realistic samples.
The implied residual variance
$\frac{1}{n}\sum_{i=1}^{n}\|\epsilon_{i}\|_{i}^{2}=\frac{1}{n}\sum_{i=1}^{n}d_{i}^{2}([y_{i}],[\mu_{i}])$
on simulated datasets ranges around 105% of the predictor variance
$\frac{1}{n}\sum_{i=1}^{n}\|h(\mathbf{x}_{i})\|_{i}^{2}=\frac{1}{n}\sum_{i=1}^{n}d_{i}^{2}([\mu_{i}],[p])$
in the form scenario and around 65% in the shape scenario. All simulations
were repeated 100 times, fitting models with the model terms specified above
and three additional nuisance effects: a linear effect $\beta z_{1}$
(orthogonal to $f_{1}(z_{1})$), an effect $f_{2}$ of the same structure as
$f_{1}$ but depending on an independently uniformly drawn variable $z_{2}$,
and a constant effect $h_{0}\in T_{[p]}\mathcal{Y}^{*}_{/G}$ to test centering
around $[p]$. Base-learners are regularized to 4 degrees of freedom (step-
length $\eta=0.1$). Early-stopping is based on 10-fold cross-validation.

Figure 4: Left: First (row 1) and second (row 2) main components of the smooth
effect $f_{1}(z_{1})$ in the form scenario obtained from TP factorization.
Normalized component directions are visualized as bottle outlines after
transporting them to the true pole (gray solid outline). Underlying truth
(orange solid lines / areas) are plotted together with five example estimates
for $n=162$ and $k=100$ (black solid lines) and the extremely sparse $k_{i}=3$
setting (gray dashed lines). Center: Conditional means for both bottle types
with fixed metric covariate $z_{1}=0$ in the shape scenario with $n=54$ and
$k=40$. Five example estimates (black solid outlines) are plotted in front of
the underlying truth (olive-green areas). Right: rMSE of shown example
estimates (jittered colored diamonds) contextualized with boxplots of rMSE
distributions observed in respective simulation scenarios.

Form scenario: In the form scenario, the smooth covariate effect $f_{1}$
offers a particularly clear interpretation. TP factorization decomposes the
true effect into its two relevant components, where the first (major)
component corresponds to the bare projection of the tilted outline in 3D into
the 2D image plane and the second to additional perspective transformations
(Fig. 4). For this effect, we observe a median relative mean squared error
$\operatorname{rMSE}(\hat{h}_{j})=\sum_{i=1}^{n}\|\hat{h}_{j}(\mathbf{x}_{i})-h_{j}(\mathbf{x}_{i})\|_{i}^{2}/\sum_{i=1}^{n}\|h(\mathbf{x}_{i})\|_{i}^{2}$
of about 3.7% of the total predictor variance for small data settings with
$n=54$ and ${k}=100$ (5.9% with $k=40$), which reduces to 1.5% for $n=162$
(for both $k=40$ and $k=100$). It is typical for functional data that, from a
certain point, adding more (highly correlated) evaluations per curve leads to
distinctly less improvement in the model fit than adding further observations
(compare, e.g., also Stöcker et al., 2021). In the extreme $k_{i}=3$ scenario,
we obtain an rMSE of around 15%, which is not surprisingly considerably higher
than for the moderate settings above. Even in this extreme setting (Fig. 4),
the effect directions are captured well, while the size of the effect is
underestimated. Rotation alignment based on only three points (which are
randomly distributed along the curves) might considerably differ from the full
curve alignment, and averaging over these sub-optimal alignments masks the
full extend of the effect. Still, results are very good given the sparsity of
information in this case. Having a simpler form, the binary effect
$\beta_{\kappa}$ is also estimated more accurately with an rMSE of around 1.5%
for $n=54$, $k=100$ (1.9% for $k=40$) and less than 0.8% for $n=162$ (for both
$k=40$ and $k=100$). The pole estimation accuracy varies on a similar scale.

Shape scenario: Qualitatively, the shape scenario shows a similar picture. For
$k=40$, we observe median rMSEs of 2.8% ($n=54$) and 2.2% ($n=162$) for
$f_{1}(z_{1})$, and 1.5% and 0.6% for the binary effect $\beta_{\kappa}$. For
$k=100$, accuracy is again slightly higher.

Nuisance effects and integration weights: Nuisance effects in the model where
generally rarely selected and, if selected at all, only lead to a marginal
loss in accuracy. The constant effect is only selected sometimes in the
extreme triangle scenarios, when pole estimation is difficult. We refer to
Brockhaus et al. (2017), who perform gradient boosting with functional
responses and a large number of covariate effects with stability selection,
for simulations with larger numbers of nuisance effects and further discussion
in a related context, as variable selection is not our main focus here.
Finally, simulations indicate that inner product weights implementing a
trapezoidal rule for numerical integration are slightly preferable for typical
grid sizes ($k=40,100$), whereas weights of $1/k_{i}$ equal over all grid
points within a curve gave slightly better results in the extreme $k_{i}=3$
settings.

All in all, the simulations show that Riemannian $L_{2}$-Boosting can
adequately fit both shape and form models in a realistic scenario and captures
effects reasonably well even for a comparably small number of sampled outlines
or evaluations per outline.

##  6 Discussion and Outlook

Compared to existing (landmark) shape regression models, the presented
approach extends linear predictors to more general additive predictors
including also, e.g., smooth nonlinear model terms and interactions, and
yields the first regression approach for functional shape as well as form
responses. Moreover, we propose novel visualizations based on TP factorization
that, similar to FPC analysis, enable a systematic decomposition of the
variability explained by an additive effect on tangent space level. Yielding
meaningful coordinates for model effects, its potential for visualization will
be useful also for FAMs in linear spaces and also beyond our model framework,
such as we exemplarily illustrate for the non-parametric approach of Jeon and
Park (2020) in Supplement S.8.

Instead of operating on the original evaluations
$\mathbf{y}_{i}\in\mathds{C}^{k_{i}}$ of response curves $y_{i}$ as in all
applications above, another frequently used approach expands $y_{i}$,
$i=1,\dots,n$, in a common basis first, before carrying out statistical
analysis on coefficient vectors (compare Ramsay and Silverman (2005); Morris
(2015) and Müller and Yao (2008) for smoothing spline, wavelet or FPC
representations in FDA or Bonhomme et al. (2014) in shape analysis).
Shape/form regression on the coefficients is, in fact, a special case of our
approach, where the inner product is evaluated on the coefficients instead of
evaluations (Supplement S.6).

The proposed model is motivated by geodesic regression. However, in the
multiple linear predictor, a linear effect of a single covariate does, in
general, not describe a geodesic for fixed non-zero values of other covariate
effects. Or put differently,
$\operatorname{Exp}_{[p]}\left(h_{1}+h_{2}\right)\neq\operatorname{Exp}_{\operatorname{Exp}_{[p]}\left(h_{1}\right)}\left(h_{2}\right)\neq\operatorname{Exp}_{\operatorname{Exp}_{[p]}\left(h_{2}\right)}\left(h_{1}\right)$
in general. Thus, hierarchical geodesic effects of the form
$\operatorname{Exp}_{\operatorname{Exp}_{[p]}\left(h_{1}\right)}\left(h_{2}\right)$,
relevant, i.a., in mixed models for hierarchical/longitudinal study designs
(Kim et al., 2017), present an interesting future extension of our model.
Moreover, an “elastic” extension based on the square-root-velocity framework
(Srivastava and Klassen, 2016) presents a promising direction for future
research, as do other manifold responses.

## Acknowledgement

We sincerely thank Nadja Pöllath for providing carefully recorded sheep
astragalus data and important insights and comments, and Sophia Schaffer for
running and discussing cell simulations and providing fully processed cell
outlines. Moreover, we gratefully acknowledge funding by grant GR 3793/3-1
from the German research foundation (DFG).

SUPPLEMENTARY MATERIAL

Supplementary material with further details is provided in an online
supplement.

## References

  * Adams et al. (2013) Adams, D., F. Rohlf, and D. Slice (2013).  A field comes of age: geometric morphometrics in the 21st century.  Hystrix, the Italian Journal of Mammalogy 24(1), 7–14.
  * Backenroth et al. (2018) Backenroth, D., J. Goldsmith, M. D. Harran, J. C. Cortes, J. W. Krakauer, and T. Kitago (2018).  Modeling motor learning using heteroscedastic functional principal components analysis.  Journal of the American Statistical Association 113(523), 1003–1015.
  * Baranyi et al. (2013) Baranyi, P., Y. Yam, and P. Várlaki (2013).  Tensor product model transformation in polytopic model-based control.  CRC press.
  * Bonhomme et al. (2014) Bonhomme, V., S. Picq, C. Gaucherel, and J. Claude (2014).  Momocs: Outline analysis using R.  Journal of Statistical Software 56(13), 1–24.
  * Brockhaus et al. (2017) Brockhaus, S., M. Melcher, F. Leisch, and S. Greven (2017).  Boosting flexible functional regression models with a high number of functional historical effects.  Statistics and Computing 27(4), 913–926.
  * Brockhaus et al. (2020) Brockhaus, S., D. Rügamer, and S. Greven (2020).  Boosting functional regression models with FDboost.  Journal of Statistical Software 94(10), 1–50.
  * Brockhaus et al. (2015) Brockhaus, S., F. Scheipl, and S. Greven (2015).  The Functional Linear Array Model.  Statistical Modelling 15(3), 279–300.
  * Bühlmann and Yu (2003) Bühlmann, P. and B. Yu (2003).  Boosting with the L2 loss: regression and classification.  Journal of the American Statistical Association 98(462), 324–339.
  * Clutton-Brock et al. (1990) Clutton-Brock, J., K. Dennis-Bryan, P. L. Armitage, and P. A. Jewell (1990).  Osteology of the Soay sheep.  Bull. Br. Mus. Nat. Hist. 56(1), 1–56.
  * Cornea et al. (2017) Cornea, E., H. Zhu, P. Kim, J. G. Ibrahim, and the Alzheimer’s Disease Neuroimaging Initiative (2017).  Regression models on Riemannian symmetric spaces.  Journal of the Royal Statistical Society: Series B 79(2), 463–482.
  * Davis et al. (2010) Davis, B. C., P. T. Fletcher, E. Bullitt, and S. Joshi (2010).  Population shape regression from random design data.  International journal of computer vision 90(2), 255–266.
  * Dryden and Mardia (2016) Dryden, I. L. and K. V. Mardia (2016).  Statistical Shape Analysis: With Applications in R.  John Wiley & Sons.
  * Ferraty et al. (2011) Ferraty, F., A. Goia, E. Salinelli, and P. Vieu (2011).  Recent advances on functional additive regression.  Recent Advances in Functional Data Analysis and Related Topics, 97–102.
  * Fletcher (2013) Fletcher, P. T. (2013).  Geodesic regression and the theory of least squares on Riemannian manifolds.  International Journal of Computer Vision 105(2), 171–185.
  * Greven and Scheipl (2017) Greven, S. and F. Scheipl (2017).  A general framework for functional regression modelling (with discussion and rejoinder).  Statistical Modelling 17(1-2), 1–35 and 100–115.
  * Happ and Greven (2018) Happ, C. and S. Greven (2018).  Multivariate functional principal component analysis for data observed on different (dimensional) domains.  Journal of the American Statistical Association 113(522), 649–659.
  * Hinkle et al. (2014) Hinkle, J., P. T. Fletcher, and S. Joshi (2014).  Intrinsic polynomials for regression on Riemannian manifolds.  Journal of Mathematical Imaging and Vision 50(1), 32–52.
  * Hofner et al. (2011) Hofner, B., T. Hothorn, T. Kneib, and M. Schmid (2011).  A framework for unbiased model selection based on boosting.  Journal of Computational and Graphical Statistics 20(4), 956–971.
  * Hofner et al. (2016) Hofner, B., T. Kneib, and T. Hothorn (2016).  A unified framework of constrained regression.  Statistics and Computing 26(1-2), 1–14.
  * Hong et al. (2014) Hong, Y., N. Singh, R. Kwitt, and M. Niethammer (2014).  Time-warped geodesic regression.  In International Conference on Medical Image Computing and Computer-Assisted Intervention, pp. 105–112. Springer.
  * Hothorn et al. (2010) Hothorn, T., P. Bühlmann, T. Kneib, M. Schmid, and B. Hofner (2010).  Model-based boosting 2.0.  Journal of Machine Learning Research 11, 2109–2113.
  * Huckemann et al. (2010) Huckemann, S., T. Hotz, and A. Munk (2010).  Intrinsic MANOVA for Riemannian manifolds with an application to Kendall’s space of planar shapes.  IEEE Transactions on Pattern Analysis and Machine Intelligence 32(4), 593–603.
  * Jeon et al. (2022) Jeon, J. M., Y. K. Lee, E. Mammen, and B. U. Park (2022).  Locally polynomial hilbertian additive regression.  Bernoulli 28(3), 2034–2066.
  * Jeon and Park (2020) Jeon, J. M. and B. U. Park (2020).  Additive regression with hilbertian responses.  The Annals of Statistics 48(5), 2671–2697.
  * Jeon et al. (2021) Jeon, J. M., B. U. Park, and I. Van Keilegom (2021).  Additive regression for non-euclidean responses and predictors.  The Annals of Statistics 49(5), 2611–2641.
  * Jupp and Kent (1987) Jupp, P. E. and J. T. Kent (1987).  Fitting smooth paths to spherical data.  Journal of the Royal Statistical Society: Series C (Applied Statistics) 36(1), 34–46.
  * Karcher (1977) Karcher, H. (1977).  Riemannian center of mass and mollifier smoothing.  Communications on Pure and Applied Mathematics 30(5), 509–541.
  * Kendall et al. (1999) Kendall, D. G., D. Barden, T. K. Carne, and H. Le (1999).  Shape and shape theory, Volume 500.  John Wiley & Sons, LTD.
  * Kent et al. (2001) Kent, J. T., K. V. Mardia, R. J. Morris, and R. G. Aykroyd (2001).  Functional models of growth for landmark data.  Proceedings in Functional and Spatial Data Analysis 109115.
  * Kim et al. (2014) Kim, H. J., N. Adluru, M. D. Collins, M. K. Chung, B. B. Bendlin, S. C. Johnson, R. J. Davidson, and V. Singh (2014).  Multivariate general linear models (mglm) on Riemannian manifolds with applications to statistical analysis of diffusion weighted images.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 2705–2712.
  * Kim et al. (2017) Kim, H. J., N. Adluru, H. Suri, B. C. Vemuri, S. C. Johnson, and V. Singh (2017).  Riemannian nonlinear mixed effects models: Analyzing longitudinal deformations in neuroimaging.  In 2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pp. 5777–5786.
  * Klingenberg (1995) Klingenberg, W. (1995).  Riemannian geometry.  de Gruyter.
  * Kneib et al. (2009) Kneib, T., T. Hothorn, and G. Tutz (2009).  Variable selection and model choice in geoadditive regression models.  Biometrics 65(2), 626–634.
  * Kume et al. (2007) Kume, A., I. L. Dryden, and H. Le (2007).  Shape-space smoothing splines for planar landmark data.  Biometrika 94(3), 513–528.
  * Lay (1967) Lay, D. M. (1967).  A study of the mammals of iran: resulting from the street expedition of 1962-63.  In Fieldiana: Zoology 54. Field Museum of Natural History.
  * Li and Ruppert (2008) Li, Y. and D. Ruppert (2008).  On the asymptotics of penalized splines.  Biometrika 95(2), 415–436.
  * Lin et al. (2017) Lin, L., B. St. Thomas, H. Zhu, and D. B. Dunson (2017).  Extrinsic local regression on manifold-valued data.  Journal of the American Statistical Association 112(519), 1261–1273.
  * Lin et al. (2020) Lin, Z., H.-G. Müller, and B. U. Park (2020).  Additive models for symmetric positive-definite matrices, Riemannian manifolds and Lie groups.  arXiv preprint arXiv:2009.08789.
  * Lutz and Bühlmann (2006) Lutz, R. W. and P. Bühlmann (2006).  Boosting for high-multivariate responses in high-dimensional linear regression.  Statistica Sinica, 471–494.
  * Mallasto and Feragen (2018) Mallasto, A. and A. Feragen (2018).  Wrapped gaussian process regression on riemannian manifolds.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 5580–5588.
  * Mayr et al. (2014) Mayr, A., H. Binder, O. Gefeller, and M. Schmid (2014).  The evolution of boosting algorithms.  Methods of information in medicine 53(06), 419–427.
  * Meyer et al. (2015) Meyer, M. J., B. A. Coull, F. Versace, P. Cinciripini, and J. S. Morris (2015).  Bayesian function-on-function regression for multilevel functional data.  Biometrics 71(3), 563–574.
  * Morris (2015) Morris, J. S. (2015).  Functional Regression.  Annual Review of Statistics and its Applications 2, 321–359.
  * Morris and Carroll (2006) Morris, J. S. and R. J. Carroll (2006).  Wavelet-based functional mixed models.  Journal of the Royal Statistical Society, Series B 68(2), 179–199.
  * Müller and Yao (2008) Müller, H.-G. and F. Yao (2008).  Functional additive models.  Journal of the American Statistical Association 103(484), 1534–1544.
  * Muralidharan and Fletcher (2012) Muralidharan, P. and P. T. Fletcher (2012).  Sasaki metrics for analysis of longitudinal data on manifolds.  In 2012 IEEE conference on computer vision and pattern recognition, pp. 1027–1034. IEEE.
  * Olsen et al. (2018) Olsen, N. L., B. Markussen, and L. L. Raket (2018).  Simultaneous inference for misaligned multivariate functional data.  Journal of the Royal Statistical Society: Series C 67(5), 1147–1176.
  * Pennec (2006) Pennec, X. (2006).  Intrinsic statistics on riemannian manifolds: Basic tools for geometric measurements.  Journal of Mathematical Imaging and Vision 25(1), 127–154.
  * Petersen and Müller (2019) Petersen, A. and H.-G. Müller (2019).  Fréchet regression for random objects with euclidean predictors.  The Annals of Statistics 47(2), 691–719.
  * Pigoli et al. (2016) Pigoli, D., A. Menafoglio, and P. Secchi (2016).  Kriging prediction for manifold-valued random fields.  Journal of Multivariate Analysis 145, 117–131.
  * Pöllath et al. (2019) Pöllath, N., R. Schafberg, and J. Peters (2019).  Astragalar morphology: Approaching the cultural trajectories of wild and domestic sheep applying geometric morphometrics.  Journal of Archaeological Science: Reports 23, 810–821.
  * R Core Team (2018) R Core Team (2018).  R: A Language and Environment for Statistical Computing.  Vienna, Austria: R Foundation for Statistical Computing.
  * Ramsay and Silverman (2005) Ramsay, J. O. and B. W. Silverman (2005).  Functional Data Analysis.  Springer New York.
  * Rosen and Thompson (2009) Rosen, O. and W. K. Thompson (2009).  A Bayesian regression model for multivariate functional data.  Computational statistics & data analysis 53(11), 3773–3786.
  * Schafberg and Wussow (2010) Schafberg, R. and J. Wussow (2010).  Julius Kühn. Das Lebenswerk eines agrarwissenschaftlichen Visionärs.  Züchtungskunde 82(6), 468–484.
  * Schaffer (2021) Schaffer, S. A. (2021).  Cytoskeletal dynamics in confined cell migration: experiment and modelling.  PhD thesis, LMU Munich.  DOI: 10.5282/edoc.28480.
  * Scheipl et al. (2015) Scheipl, F., A.-M. Staicu, and S. Greven (2015).  Functional additive mixed models.  Journal of Computational and Graphical Statistics 24(2), 477–501.
  * Schiratti et al. (2017) Schiratti, J.-B., S. Allassonnière, O. Colliot, and S. Durrleman (2017).  A bayesian mixed-effects model to learn trajectories of changes from repeated manifold-valued observations.  The Journal of Machine Learning Research 18(1), 4840–4872.
  * Shi et al. (2009) Shi, X., M. Styner, J. Lieberman, J. G. Ibrahim, W. Lin, and H. Zhu (2009).  Intrinsic regression models for manifold-valued data.  In International Conference on Medical Image Computing and Computer-Assisted Intervention, pp. 192–199. Springer.
  * Srivastava and Klassen (2016) Srivastava, A. and E. P. Klassen (2016).  Functional and Shape Data Analysis.  Springer-Verlag.
  * Stöcker et al. (2021) Stöcker, A., S. Brockhaus, S. A. Schaffer, B. v. Bronk, M. Opitz, and S. Greven (2021).  Boosting functional response models for location, scale and shape with an application to bacterial competition.  Statistical Modelling 21(5), 385–404.
  * Stöcker et al. (2022) Stöcker, A., M. Pfeuffer, L. Steyer, and S. Greven (2022).  Elastic full Procrustes analysis of plane curves via Hermitian covariance smoothing.
  * Thüroff et al. (2019) Thüroff, F., A. Goychuk, M. Reiter, and E. Frey (2019, dec).  Bridging the gap between single-cell migration and collective dynamics.  eLife 8, e46842.
  * Volkmann et al. (2021) Volkmann, A., A. Stöcker, F. Scheipl, and S. Greven (2021).  Multivariate functional additive mixed models.  Statistical Modelling.
  * Wood et al. (2016) Wood, S. N., N. Pya, and B. Säfken (2016).  Smoothing parameter and model selection for general smooth models.  Journal of the American Statistical Association 111(516), 1548–1563.
  * Yao et al. (2005) Yao, F., H. Müller, and J. Wang (2005).  Functional data analysis for sparse longitudinal data.  Journal of the American Statistical Association 100(470), 577–590.
  * Zeder (2006) Zeder, M. A. (2006).  Reconciling rates of long bone fusion and tooth eruption and wear in sheep (Ovis) and goat (Capra).  Recent advances in ageing and sexing animal bones 9, 87–118.
  * Zhu et al. (2009) Zhu, H., Y. Chen, J. G. Ibrahim, Y. Li, C. Hall, and W. Lin (2009).  Intrinsic regression models for positive-definite matrices with applications to diffusion tensor imaging.  Journal of the American Statistical Association 104(487), 1203–1212.
  * Zhu et al. (2012) Zhu, H., R. Li, and L. Kong (2012).  Multivariate varying coefficient model for functional responses.  Annals of statistics 40(5), 2634–2666.
  * Zhu et al. (2017) Zhu, H., J. S. Morris, F. Wei, and D. D. Cox (2017).  Multivariate functional response regression, with application to fluorescence spectroscopy in a cervical pre-cancer study.  Computational Statistics and Data Analysis 111, 88–101.

##  Appendix S Online Supplementary Material to

Functional additive models on manifolds of planar shapes and forms
by Almond Stöcker, Lisa Steyer and Sonja Greven

###  S.1 Geometry of functional forms and shapes

####  S.1.1 Translation, rotation and re-scaling as normal subgroups

We consider the following invariances of a response curve $y\in\mathcal{Y}$
with respect to the transformations $\mathcal{Y}\rightarrow\mathcal{Y}$ given
by the group actions of translation
$\operatorname{Trl}=\\{y\overset{\operatorname{Trl}_{\gamma}}{\longmapsto}y+\gamma\,\\!\text{1}:\gamma\in\mathds{C}\\}$
with some $\\!\text{1}\in\mathcal{Y}\setminus\\{\text{0}\\}$ (for curves
typically $\\!\text{1}:t\mapsto\frac{1}{\|t\mapsto 1\|}$ the real constant
function of unit norm), re-scaling
$\operatorname{Scl}=\\{y\overset{\operatorname{Scl}_{\lambda}}{\longmapsto}\lambda\cdot(y-\text{0}_{y})+\text{0}_{y}:\lambda\in\mathds{R}^{+}\\}$
around a reference point $\text{0}_{y}\in\mathds{C}$, and rotation
$\operatorname{Rot}=\\{y\overset{\operatorname{Rot}_{u}}{\longmapsto}u\cdot(y-\text{0}_{y})+\text{0}_{y}:u\in\mathds{S}^{1}\\}$
around $\text{0}_{y}$ with
$\mathds{S}^{1}=\\{u\in\mathds{C}:|u|=1\\}=\\{\exp(\omega\sqrt{\text{-}1}):\omega\in\mathds{R}\\}$
the circle group reflecting counterclockwise rotations by $\omega$ radian
measure. In the literature, the reference point is usually omitted setting
$\text{0}_{y}=\text{0}$, which can be done without loss of generality under
translation invariance (i.e. in particular for shapes/forms). However, keeping
other possible combinations of invariances in mind, we explicitly refer to an
individual reference point and suggest the centroid
$\text{0}_{y}=\langle{\\!\text{1}}\,,y\rangle\\!\text{1}$ or, more generally,
$\text{0}_{y}=a(y)\,\\!\text{1}$ for some linear functional
$a:\mathcal{Y}\rightarrow\mathds{R}$. Assuming
$\operatorname{Trl}_{\gamma}(0_{y})=0_{\operatorname{Trl}_{\gamma}(y)}$, as
for the centroid, the definition of re-scaling and rotation around
$\text{0}_{y}$ ensures that $\operatorname{Trl}_{\gamma}$,
$\operatorname{Scl}_{\lambda}$ and $\operatorname{Rot}_{u}$ commute – and that
$\operatorname{Trl}$, $\operatorname{Rot}$ and $\operatorname{Scl}$ present
normal subgroups of the combined group actions $\\{y\mapsto\lambda
uy+\gamma:\gamma\in\mathds{C},\lambda\in\mathds{R}^{+},u\in\mathds{S}^{1}\\}$
of shape invariances. Thus, the combined group actions can be written as the
direct product (or direct sum)
$\operatorname{Trl}\times\operatorname{Scl}\times\operatorname{Rot}=\\{\operatorname{Trl}_{\gamma}\circ\operatorname{Scl}_{\lambda}\circ\operatorname{Rot}_{u}:\gamma\in\mathds{C},\lambda\in\mathds{R}^{+},u\in\mathds{S}^{1}\\}\cong\mathds{C}\times\mathds{R}^{+}\times\mathds{S}^{1}$
and invariances with respect to $\operatorname{Trl}$, $\operatorname{Scl}$,
$\operatorname{Rot}$ can be modularly accounted for in arbitrary order.
$\operatorname{Trl}\times\operatorname{Rot}$, for instance, describe rigid
motions. The ultimate response object is then given by the _orbit_
$[y]_{G}=\\{g(y):g\in G\\}$ (or short $[y]$), i.e. the equivalence class with
respect to the direct sum $G$ generated by the chosen combination of
$\operatorname{Trl}$, $\operatorname{Scl}$ and $\operatorname{Rot}$.
$[y]_{\operatorname{Trl}\times\operatorname{Scl}\times\operatorname{Rot}}$ is
referred to as the _shape_ of $y$ and
$[y]_{\operatorname{Trl}\times\operatorname{Rot}}$ as its _form_ or _size-and-
shape_ [compare][]DrydenMardia2016ShapeAnalysisWithApplications; studying
$[y]_{\operatorname{Scl}}$ is closely related to _directional_ data analysis
MardiaJupp2009DirectionalStatistics where the direction of $y$ is analyzed
independent of its size $\|y\|$.

####  S.1.2 Parallel transport of form tangent vectors

To confirm that the parallel transport in the form space
$\mathcal{Y}^{*}_{/\operatorname{Trl}\times\operatorname{Rot}}$ can be carried
out via representatives in $\mathcal{Y}$ as described in the main manuscript,
we closely follow Huckemann2010intrinsicMANOVA in their derivation of shape
parallel transport. Necessary differential geometric notions and statements
are briefly introduced in the following before stating the main result in
Lemma 1. For a more profound introduction, we recommend [in particular, p. 21,
43, 93, 124, 337, 402]Lee2018RiemannianManifolds, as well as tu2011manifolds
for an illustrative introduction into some of the concepts, and [in
particular, p.103- 107]Klingenberg1995RiemannianGeometry for an introduction
in the light of potentially infinite dimensional manifolds.

The entire argument crucially relies on properties known for Riemannian
submersions between differentiable manifolds $\widetilde{\mathcal{M}}$ and
$\mathcal{M}$, which allow to relate the structure of $\mathcal{M}$ back to
$\widetilde{\mathcal{M}}$. A submersion is a smooth surjective function
$\Phi:\widetilde{\mathcal{M}}\rightarrow\mathcal{M}$, for which also the
differential $d\Phi:T_{\widetilde{q}}\widetilde{\mathcal{M}}\rightarrow
T_{q}\mathcal{M}$, $\widetilde{q}\in\widetilde{\mathcal{M}}$,
$q=\Phi{(\widetilde{q})}$, is surjective at each
$\widetilde{q}\in\widetilde{\mathcal{M}}$. For $q\in\mathcal{M}$, the
$\Phi^{-1}(\\{q\\})$ are submanifolds of $\mathcal{M}$, and
$T_{\widetilde{q}}\widetilde{\mathcal{M}}=T_{\widetilde{q}}\Phi^{-1}(\\{q\\})\oplus
H_{\widetilde{q}}\widetilde{\mathcal{M}}$ can be decomposed into the vertical
space $T_{\widetilde{q}}\Phi^{-1}(\\{q\\})=\ker(d\Phi)$ and its orthogonal
complement $H_{\widetilde{q}}\widetilde{\mathcal{M}}$, the horizontal space.
When restricted to the horizontal space,
$d\Phi\big{|}_{H_{\widetilde{q}}\widetilde{\mathcal{M}}}:H_{\widetilde{q}}\widetilde{\mathcal{M}}\rightarrow
T_{q}\mathcal{M}$ presents a linear isomorphism. A submersion $\Phi$ is called
Riemannian submersion if
$d\Phi\big{|}_{H_{\widetilde{q}}\widetilde{\mathcal{M}}}$ is also isometric.
It gives rise to an identification $T_{q}\mathcal{M}\cong
H_{\widetilde{q}}\widetilde{\mathcal{M}}$ of tangent spaces of $\mathcal{M}$
with horizontal spaces on $\widetilde{\mathcal{M}}$. Such an identification
underlies the presentation of the response geometry in Section 2 of the main
manuscript.
By construction, the quotient map
$\Phi:\mathcal{Y}^{*}\rightarrow\mathcal{Y}^{*}_{/\operatorname{Trl}\times\operatorname{Rot}},y\mapsto[y]$
presents a Riemannian submersion: Since
$[p]=\\{up+\gamma\\!\text{1}:u\in\mathds{S}^{1},\gamma\in\mathds{C}\\}$ embeds
$\mathds{S}^{1}\times\mathds{R}^{2}$ in $\mathcal{Y}$, and, since the tangent
spaces of $\mathds{S}^{1}$ and $\mathds{R}^{2}$ are well-known, the vertical
space is given by $T_{p}[p]\cong\\{\lambda
p\sqrt{\text{-}1}+\gamma\,\\!\text{1}:\lambda\in\mathds{R},\gamma\in\mathds{C}\\}\subset\mathcal{Y}$,
with orthogonal complement
$H_{p}\mathcal{Y}^{*}\cong\\{y\in\mathcal{Y}:\langle y,\\!\text{1}\rangle=0,\
\operatorname{Im}\\!\left(\langle y,p\rangle\right)=0\\}$ (see also Figure 1
in the main manuscript for an illustration). While $\Phi$ is obviously
surjective, surjectivity and isometry of $d\Phi|_{H_{p}\mathcal{Y}^{*}}$ can
be seen by expressing $\Phi$ in terms of the charts for
$\mathcal{Y}^{*}_{/\operatorname{Trl}\times\operatorname{Rot}}$: for a given
$p\in[p]\in\mathcal{Y}_{/\operatorname{Trl}\times\operatorname{Rot}}^{*}$, the
map
$\widetilde{(\cdot)}:[y]\mapsto\widetilde{y}^{\operatorname{Trl}\times\operatorname{Rot}}$
provides a chart $\mathcal{U}_{[p]}\rightarrow\mathcal{V}_{p}$, i.e. an
isomorphism from
$\mathcal{U}_{[p]}=\\{[y]\in\mathcal{Y}_{/\operatorname{Trl}\times\operatorname{Rot}}^{*}:\langle\widetilde{p}^{\operatorname{Trl}},\widetilde{y}^{\operatorname{Trl}}\rangle\neq
0\\}$ to $\mathcal{V}_{p}=\\{y\in\mathcal{Y}:\operatorname{Im}\\!\left(\langle
p,y\rangle\right)=0,\operatorname{Re}\\!\left(\langle
p,y\rangle\right)>0,\langle\\!\text{1},y\rangle=0\\}$ used to establish the
differential structure on
$\mathcal{Y}_{/\operatorname{Trl}\times\operatorname{Rot}}^{*}$. Expressed in
this chart,
$\widetilde{\Phi}(y)=\widetilde{(\cdot)}\circ\Phi(y)=\widetilde{y}^{\operatorname{Trl}\times\operatorname{Rot}}$
is the identity for all
$y\in\mathcal{V}_{p}\subset\Phi^{-1}\left(\mathcal{U}_{[p]}\right)$. Thus,
since $T_{p}\mathcal{V}_{p}=H_{p}\mathcal{Y}^{*}$, also
$d\widetilde{\Phi}\big{|}_{H_{p}\mathcal{Y}^{*}}$ is the identity, which is
obviously an isometric isomorphism. The latter carries over to
$d\Phi\big{|}_{H_{p}\mathcal{Y}^{*}}$ independent of the given chart.
The isometric isomorphism
$d\Phi\big{|}_{H_{p}\mathcal{Y}^{*}}:H_{p}\mathcal{Y}^{*}\rightarrow
T_{[p]}\mathcal{Y}_{/\operatorname{Trl}\times\operatorname{Rot}}^{*}$ yields
the identification
$T_{[p]}\mathcal{Y}_{/\operatorname{Trl}\times\operatorname{Rot}}^{*}\cong\\{y\in\mathcal{Y}:\langle
y,\\!\text{1}\rangle=0,\ \operatorname{Im}\\!\left(\langle
y,p\rangle\right)=0\\}$, which we rely on in the main manuscript. Unlike
there, we denote
$d\Phi\big{|}_{H_{y}\mathcal{Y}^{*}}^{-1}:\xi\mapsto\widetilde{\xi}$ also for
tangent vectors in the following, to make the identification of
$\xi=d\Phi(\widetilde{\xi})\in
T_{[y]}\mathcal{Y}_{/\operatorname{Trl}\times\operatorname{Rot}}^{*}$ with the
corresponding $\widetilde{\xi}\in H_{y}\mathcal{Y}^{*}$, usually referred to
as horizontal lift, explicit in the notation.

The covariant derivative (Levi-Civita connection)
$\nabla^{\mathcal{M}}_{V}W\in T\mathcal{M}$ of a vector-field $W\in
T\mathcal{M}$ along a vector-field $V\in T\mathcal{M}$ provides a derivative
of vector-fields in the tangent bundle
$T\mathcal{M}=\\{T_{q}\mathcal{M}:q\in\mathcal{M}\\}$ of a Riemannian manifold
$\mathcal{M}$. As a derivation in $W$ and a linear function in $V$,
$\nabla^{\mathcal{M}}$ fulfills a set of properties identifying it as unique
generalization of ordinary directional derivatives of the components of
$W:q\mapsto W_{q}\in T_{q}\mathcal{M}$ into the direction $V_{q}\in
T_{q}\mathcal{M}$. For a submanifold $\mathcal{M}$ of a linear space
$\mathcal{Y}$, $\nabla^{\mathcal{M}}_{V}W$ corresponds to the ordinary
directional derivative orthogonally projected into $T_{q}\mathcal{M}$. For the
linear case (with $\mathcal{M}=\mathcal{Y}$), the covariant derivative of a
vector field $W(\tau):=W_{c{(\tau)}}$ along a differentiable curve $c(\tau)$
is directly given as

$\nabla^{\mathcal{Y}}_{\dot{c}(t)}W(\tau)=\dot{W}(\tau)=\frac{d}{d\tau}W(\tau).$
(2)

In analogy to straight lines, geodesic curves $c(\tau)$ are characterized by

$\nabla^{\mathcal{M}}_{\dot{c}(\tau)}\dot{c}(\tau)=\text{0},$

i.e. curves with zero ‘second derivative’. More generally, a vector-field $W$
is called parallel along a curve $c(\tau)$ if

$\nabla^{\mathcal{M}}_{\dot{c}(\tau)}W(\tau)=\text{0}.$ (3)

According to that the parallel transport
$\operatorname{Transp}_{q,q^{\prime}}^{c}:T_{q}\mathcal{M}\rightarrow
T_{q^{\prime}}\mathcal{M}$ along a curve
$c:[\tau_{0},\tau_{1}]\rightarrow\mathcal{M}$ between $c(\tau_{0})=q,\
c(\tau_{1})=q^{\prime}\in\mathcal{M}$ is defined to map tangent vectors
$\varepsilon=W(\tau_{0})\mapsto\varepsilon^{\prime}=W(\tau_{1})$ for some
vector field $W$ parallel along $c$ (fulfilling Equation 3). If the curve $c$
is clear from context, we omit it in the notation. This is especially the case
in the following, where $c$ can be chosen as the unique geodesic between two
forms $[p]$ and $[p^{\prime}]$ with $\langle p,p^{\prime}\rangle\neq 0$,
yielding a canonical connection (in this case, $c$ corresponds to the line
between $p$ and the aligned $\tilde{p}^{\prime}$; for $\langle
p,p^{\prime}\rangle=0$, by contrast, it is easy to see that for each
$u\in\mathds{S}^{1}$ the line between $p$ and $up^{\prime}$ corresponds to a
different geodesic; the second case can, however, be neglected).

The possibility to effectively carry out the parallel transport between forms
$[p],[p^{\prime}]$ on suitable representatives
$p,p^{\prime}\in\mathcal{Y}^{*}$ stems from the following theorem and
subsequent Corollary [compare, e.g,][p.
103-105]Klingenberg1995RiemannianGeometry.

######  Theorem 1.

Let $\Phi:\widetilde{\mathcal{M}}\rightarrow\mathcal{M}$ be a Riemannian
submersion between manifolds $\widetilde{\mathcal{M}}$ and $\mathcal{M}$, and
$V,W\in T\mathcal{M}$ vector-fields. Then

$\nabla^{\widetilde{\mathcal{M}}}_{\widetilde{V}}\widetilde{W}=\widetilde{\left(\nabla^{\mathcal{M}}_{V}W\right)}+\frac{1}{2}[\widetilde{V},\widetilde{W}]^{\perp}$

where $\widetilde{Z}\in H\widetilde{\mathcal{M}}$ denotes the horizontal lift
of $Z\in T\mathcal{M}$ to the horizontal bundle
$H\widetilde{\mathcal{M}}=\\{H_{\widetilde{p}}\widetilde{\mathcal{M}}:\widetilde{p}\in\widetilde{\mathcal{M}}\\}$,
$Z=d\Phi\left(\widetilde{Z}\right)$, and
$[\widetilde{V},\widetilde{W}]^{\perp}$ is the the Lie bracket
$[\widetilde{V},\widetilde{W}]=\widetilde{V}\circ\widetilde{W}-\widetilde{W}\circ\widetilde{V}$
orthogonally projected $(\cdot)^{\perp}:T\widetilde{\mathcal{M}}\rightarrow
ker\left(d\Phi\right)$ to the vertical space.

######  Corollary 1.

Let $\Phi:\widetilde{\mathcal{M}}\rightarrow\mathcal{M}$ be a Riemannian
submersion and $c:(\tau_{0},\tau_{1})\rightarrow\mathcal{M}$ a smooth curve on
$\mathcal{M}$ with
$\widetilde{c}:(\tau_{0},\tau_{1})\rightarrow\widetilde{\mathcal{M}}$ its
horizontal lift, i.e., $\Phi\circ\widetilde{c}=c$,
$d\Phi\circ\dot{\widetilde{c}}=\dot{c}$ and $\dot{\widetilde{c}}(\tau)\in
H_{\widetilde{c}(\tau)}\widetilde{\mathcal{M}}$ horizontal (i.e.
$\dot{\widetilde{c}}=\widetilde{\dot{c}}$). Then

  1. i)

a vector-field $W=d\Phi\circ\widetilde{W}\in T\mathcal{M}$ along $c$ is
parallel if and only if

$\nabla_{\dot{\widetilde{c}}}^{\widetilde{\mathcal{M}}}\widetilde{W}=\frac{1}{2}[\dot{\widetilde{c}},\widetilde{W}]^{\perp}$

for the horizontal vector-field $\widetilde{W}\in T\widetilde{\mathcal{M}}$
along $\widetilde{c}$.

  2. ii)

$c$ is a geodesic if and only if $\widetilde{c}$ is a geodesic.

While i) yields the basis for confirming the parallel transport computation,
ii) is the underlying fact behind the identification of geodesics in form and
shape spaces with geodesics of suitably aligned representatives. Note that,
while Huckemann2010intrinsicMANOVA generally restrict their discussion to
finite dimensional manifolds, the theorem does in fact not have this
restriction. Based on these preparations, we can now verify the presented
parallel transport along the lines of Huckemann2010intrinsicMANOVA, but for
forms rather than shapes and explicitly based on a separable Hilbert space
$\mathcal{Y}$ rather than on $\mathds{C}^{k}$. Note that, while identifying
$H_{p}\mathcal{Y}^{*}\cong
T_{[p]}\mathcal{Y}_{/\operatorname{Trl}\times\operatorname{Rot}}^{*}$ and
$\varepsilon\cong d\Phi(\varepsilon)$ in the main manuscript, they are
distinguished here for clarity.

######  Lemma 1.

Let $p,p^{\prime}\in\mathcal{Y}^{*}$ with $\langle p,p^{\prime}\rangle\neq 0$
centered and mutually rotation aligned representatives of forms
$[p],[p^{\prime}]\in\mathcal{Y}^{*}_{/\operatorname{Trl}\times\operatorname{Rot}}$
(i.e. $p=\widetilde{p}^{{\operatorname{Trl}\times\operatorname{Rot}}}$ for
notational simplicity and $p^{\prime}$ accordingly), let $\varepsilon\in
T_{[p]}\mathcal{Y}^{*}_{/\operatorname{Trl}\times\operatorname{Rot}}$ with
horizontal lift $\tilde{\varepsilon}\in H_{p}\mathcal{Y}^{*}$, and let
$\Phi:y\mapsto[y]$ denote the quotient map. Then

$\operatorname{Transp}_{[p],[p^{\prime}]}\left(\varepsilon\right)=d\Phi\left(\tilde{\varepsilon}-\operatorname{Im}\\!\left(\langle
p^{\prime}/\|p^{\prime}\|,\tilde{\varepsilon}\rangle\right)\frac{p/\|p\|+p^{\prime}/\|p^{\prime}\|}{1+\langle
p/\|p\|,p^{\prime}/\|p^{\prime}\|\rangle}\sqrt{\text{-}1}\right)$ (4)

implements the form parallel transport via its horizontal lift.

###### Proof.

For $\langle p,p^{\prime}\rangle\neq 0$ aligned and centered, the unique unit-
speed geodesic (uniqueness can be seen using Corollary 1 ii)) between $[p]$
and $[p^{\prime}]$ is described by
$\tau\rightarrow[p+\tau\frac{p^{\prime}-p}{\|p^{\prime}-p\|}]$. Yet, to
simplify the argument, we choose a unit-angular speed parameterization
instead. It takes the form
$c(\tau):=[\widetilde{c}(\tau)]:=[\rho(\tau)\gamma(\tau)]$ with
$\gamma(\tau)=\cos(\tau)\beta+\sin(\tau)\beta^{\prime}$ where
$\beta=\frac{p}{\|p\|}$ and
$\beta^{\prime}=\frac{p^{\prime}-\langle\beta,p^{\prime}\rangle\beta}{\|p^{\prime}-\langle\beta,p^{\prime}\rangle\beta\|}=\frac{\frac{p^{\prime}}{\|p^{\prime}\|}-\langle\frac{p}{\|p\|},\frac{p^{\prime}}{\|p^{\prime}\|}\rangle\frac{p}{\|p\|}}{\|\frac{p^{\prime}}{\|p^{\prime}\|}-\langle\frac{p}{\|p\|},\frac{p^{\prime}}{\|p^{\prime}\|}\rangle\frac{p}{\|p\|}\|}$
form an orthonormal basis of the real plain containing the horizontal
geodesic. With $\widetilde{c}(0)=p$ and
$\widetilde{c}(\arccos\langle\frac{p}{\|p\|},\frac{p^{\prime}}{\|p^{\prime}\|}\rangle)=p^{\prime}$,
$\widetilde{c}(\tau)$ describes the line connecting $p$ and $p^{\prime}$ in
polar coordinates.
$[\gamma(\tau)]_{\operatorname{Trl}\times\operatorname{Rot}\times\operatorname{Scl}}$
corresponds to the shape geodesic between
$[p]_{\operatorname{Trl}\times\operatorname{Rot}\times\operatorname{Scl}}$ and
$[p^{\prime}]_{\operatorname{Trl}\times\operatorname{Rot}\times\operatorname{Scl}}$,
and $\rho(\tau)=\|\widetilde{c}(\tau)\|$ reflects the size of the geodesic
$c(\tau)$. An explicit definition of $\rho(\tau)$ is not needed.

Due to the alignment of $p$ and $p^{\prime}$, $\dot{\gamma}(\tau)$ and also

$\widetilde{W}(\tau):=\tilde{\varepsilon}+\operatorname{Im}\\!\left(\langle\beta^{\prime},\tilde{\varepsilon}\rangle\right)\left(\dot{\gamma}(\tau)-\beta^{\prime}\right)\sqrt{\text{-}1}$
(5)

are horizontal along $\widetilde{c}(\tau)$, i.e. $\widetilde{W}(\tau)\in
H_{\widetilde{c}(\tau)}\mathcal{Y}^{*}$ for each $\tau$, if
$\tilde{\varepsilon}$ is horizontal, i.e. if
$\operatorname{Im}\\!\left(\langle
p,\tilde{\varepsilon}\rangle\right)=\langle\\!\text{1},\tilde{\varepsilon}\rangle=0$.
More concretely, this holds as

$\displaystyle\operatorname{Im}\\!\left(\langle\widetilde{c}(\tau),\widetilde{W}(\tau)\rangle\right)$
$\displaystyle=\rho(\tau)\;\left(\operatorname{Im}\\!\left(\langle\gamma(\tau),\tilde{\varepsilon}\rangle\right)+\operatorname{Im}\\!\left(\langle\beta^{\prime},\tilde{\varepsilon}\rangle\right)\operatorname{Re}\\!\left(\langle\gamma(\tau),\dot{\gamma}(\tau)-\beta^{\prime}\rangle\right)\right)$
$\displaystyle\overset{\tilde{\varepsilon}\text{
horizontal}}{=}\rho(\tau)\;\left(\sin(\tau)\operatorname{Im}\\!\left(\langle\beta^{\prime},\tilde{\varepsilon}\rangle\right)+\operatorname{Im}\\!\left(\langle\beta^{\prime},\tilde{\varepsilon}\rangle\right)(0-\sin(\tau)\underbrace{\|\beta^{\prime}\|^{2}}_{=1})\right)=0$

and, obviously, also $\langle\\!\text{1},\widetilde{W}(\tau)\rangle=0$ as this
is the case for all involved vectors. Moreover, $\widetilde{W}$ is smooth and
$\tilde{\varepsilon}\mapsto\widetilde{W}\left(\arccos\langle\frac{p}{\|p\|},\frac{p^{\prime}}{\|p^{\prime}\|}\rangle\right)$
yields the transport formulated in Equation (4), which follows from basic
trigonometric relations. In detail, it follows from plugging

$\displaystyle\dot{\gamma}\left(\arccos\langle\frac{p}{\|p\|},\frac{p^{\prime}}{\|p^{\prime}\|}\rangle\right)-\beta^{\prime}$
$\displaystyle=\langle\frac{p}{\|p\|},\frac{p^{\prime}}{\|p^{\prime}\|}\rangle\
\beta^{\prime}-\sqrt{1-\langle\frac{p}{\|p\|},\frac{p^{\prime}}{\|p^{\prime}\|}\rangle^{2}}\
\beta-\beta^{\prime}$
$\displaystyle=\left(\langle\frac{p}{\|p\|},\frac{p^{\prime}}{\|p^{\prime}\|}\rangle-1\right)\overbrace{\frac{\frac{p^{\prime}}{\|p^{\prime}\|}-\langle\frac{p}{\|p\|},\frac{p^{\prime}}{\|p^{\prime}\|}\rangle\frac{p}{\|p\|}}{\sqrt{1-\langle\frac{p}{\|p\|},\frac{p^{\prime}}{\|p^{\prime}\|}\rangle^{2}}}}^{\beta^{\prime}}\
-\sqrt{1-\langle\frac{p}{\|p\|},\frac{p^{\prime}}{\|p^{\prime}\|}\rangle^{2}}\
\frac{p}{\|p\|}$
$\displaystyle=\frac{\left(\langle\frac{p}{\|p\|},\frac{p^{\prime}}{\|p^{\prime}\|}\rangle-1\right)\frac{p^{\prime}}{\|p^{\prime}\|}}{\sqrt{1-\langle\frac{p}{\|p\|},\frac{p^{\prime}}{\|p^{\prime}\|}\rangle^{2}}}$
$\displaystyle\qquad+\frac{-\langle\frac{p}{\|p\|},\frac{p^{\prime}}{\|p^{\prime}\|}\rangle^{2}\frac{p}{\|p\|}+\langle\frac{p}{\|p\|},\frac{p^{\prime}}{\|p^{\prime}\|}\rangle\frac{p}{\|p\|}-\left(1-\langle\frac{p}{\|p\|},\frac{p^{\prime}}{\|p^{\prime}\|}\rangle^{2}\right)\frac{p}{\|p\|}}{\sqrt{1-\langle\frac{p}{\|p\|},\frac{p^{\prime}}{\|p^{\prime}\|}\rangle^{2}}}$
$\displaystyle=\frac{-\left(1-\langle\frac{p}{\|p\|},\frac{p^{\prime}}{\|p^{\prime}\|}\rangle\right)\left(\frac{p^{\prime}}{\|p^{\prime}\|}+\frac{p}{\|p\|}\right)}{\sqrt{1-\langle\frac{p}{\|p\|},\frac{p^{\prime}}{\|p^{\prime}\|}\rangle^{2}}}$

and

$\displaystyle\operatorname{Im}\\!\left(\langle\beta^{\prime},\tilde{\varepsilon}\rangle\right)\overset{\text{$\tilde{\varepsilon}$
horizontal}}{=}\frac{\operatorname{Im}\\!\left(\langle
p^{\prime},\tilde{\varepsilon}\rangle\right)}{\sqrt{1-\langle\frac{p}{\|p\|},\frac{p^{\prime}}{\|p^{\prime}\|}\rangle^{2}}}$

into the definition of $W(\tau)=d\Phi(\widetilde{W}(\tau))$ using (5).

Hence, due to Corollary 1 i), we mainly need to show

$\nabla_{\dot{\widetilde{c}}}^{\mathcal{Y}^{*}}\widetilde{W}=\frac{1}{2}[\dot{\widetilde{c}},\widetilde{W}]^{\perp}$
(6)

where the left-hand side may directly be computed as

$\displaystyle\nabla_{\widetilde{\dot{c}}(\tau)}^{\mathcal{Y}^{*}}\widetilde{W}(\tau)$
$\displaystyle\overset{\eqref{covariant_derivative}}{=}\dot{\widetilde{W}}(\tau)\overset{\eqref{W}}{=}-\operatorname{Im}\\!\left(\langle\beta^{\prime},\tilde{\varepsilon}\rangle\right)\gamma(\tau)\sqrt{\text{-}1}$

since $\ddot{\gamma}(\tau)=-\gamma(\tau)$.

On the right-hand side, the orthogonal projection of a vector-field
$V(\tau):=V_{\widetilde{c}(\tau)}\in T_{\widetilde{c}(\tau)}\mathcal{Y}^{*}$
along $\widetilde{c}(\tau)$ into the vertical spaces (of which
$\\{\gamma(\tau),\\!\text{1},\sqrt{\text{-}1}\\!\text{1}\\}$ constitute an
orthonormal basis) is given by

$\displaystyle V^{\perp}(\tau)$
$\displaystyle=\frac{\operatorname{Re}\\!\left(\langle\sqrt{\text{-}1}\,\widetilde{c}(\tau),V(\tau)\rangle\right)}{\|\widetilde{c}(\tau)\|^{2}}\widetilde{c}(\tau)\sqrt{\text{-}1}+\langle\\!\text{1},V(\tau)\rangle\,\\!\text{1}$
$\displaystyle=\frac{\omega^{\operatorname{Rot}}(V(\tau))}{\rho(\tau)}\gamma(\tau)\sqrt{\text{-}1}+\omega^{\operatorname{Trl}}(V(\tau))\,\\!\text{1}$

with the 1-forms $\omega^{\operatorname{Rot}}$ and
$\omega^{\operatorname{Trl}}$ defined as

$\displaystyle\omega^{\operatorname{Rot}}(V_{p})$
$\displaystyle:=\operatorname{Re}\\!\left(\langle\sqrt{\text{-}1}\,p,V_{p}\rangle\right)=\operatorname{Re}\\!\left(-\sqrt{\text{-}1}\,\langle
p,V_{p}\rangle\right)$
$\displaystyle=\operatorname{Re}\\!\left(-\sqrt{\text{-}1}\,\left(\operatorname{Re}\\!\left(\langle
p,V_{p}\rangle\right)+\operatorname{Im}\\!\left(\langle
p,V_{p}\rangle\right)\sqrt{\text{-}1}\right)\right)$
$\displaystyle=\operatorname{Re}\\!\left(-\sqrt{\text{-}1}\,\operatorname{Re}\\!\left(\langle
p,V_{p}\rangle\right)+\operatorname{Im}\\!\left(\langle
p,V_{p}\rangle\right)\right)$ $\displaystyle=\operatorname{Im}\\!\left(\langle
p,V_{p}\rangle\right).$

and $\omega^{\operatorname{Trl}}(V_{p})=\langle\\!\text{1},V_{p}\rangle$ for
$p\in\mathcal{Y}^{*}$.

Thus, to confirm (6) and complete the proof, it remains to show
$\omega^{\operatorname{Rot}}([\dot{\widetilde{c}},\widetilde{W}])=-2\operatorname{Im}\\!\left(\langle\beta^{\prime},\tilde{\varepsilon}\rangle\right)\rho$
and $\omega^{\operatorname{Trl}}([\dot{\widetilde{c}},\widetilde{W}])=0$. For
this, we use some statements on the exterior derivative $d\omega$ of a 1-form
$\omega$ subsumed in the following auxiliary lemma (proven later):

######  Lemma 2.

Let $V,W$ be smooth vector-fields.

  1. i)

For any smooth 1-form $\omega$ it holds that
$\omega\left([V,W]\right)=V\left(\omega\left(W\right)\right)-W\left(\omega\left(V\right)\right)-d\omega\left(V,W\right)$.

  2. ii)

For $\omega^{\operatorname{Rot}}$ defined above,
$d\omega^{\operatorname{Rot}}(V,W)=2\,\operatorname{Im}\\!\left(\langle
V,W\rangle\right)$.

  3. iii)

For $\omega^{\operatorname{Trl}}$ defined above,
$d\omega^{\operatorname{Trl}}(V,W)=0$.

Using further that

$\displaystyle\omega^{\operatorname{Rot}}\left(\dot{\widetilde{c}}(\tau)\right)$
$\displaystyle=\operatorname{Im}\\!\left(\langle\gamma(\tau),\rho(\tau)\,\dot{\gamma}(\tau)\rangle\right)+\operatorname{Im}\\!\left(\langle\gamma(\tau),\dot{\rho}(\tau)\,\gamma(\tau)\rangle\right)=0$

and

$\displaystyle\omega^{\operatorname{Rot}}\left(\widetilde{W}(\tau)\right)$
$\displaystyle=\rho(\tau)\operatorname{Im}\\!\left(\langle\gamma(\tau),\tilde{\varepsilon}\rangle\right)+\rho(\tau)\operatorname{Im}\\!\left(\langle\beta^{\prime},\tilde{\varepsilon}\rangle\right)\underbrace{\langle\gamma(\tau),\dot{\gamma}(\tau)-\beta^{\prime}\rangle}_{\in\mathds{R}}$
$\displaystyle=\rho(\tau)\left(\sin(\tau)\operatorname{Im}\\!\left(\langle\beta^{\prime},\tilde{\varepsilon}\rangle\right)-\operatorname{Im}\\!\left(\langle\beta^{\prime},\tilde{\varepsilon}\rangle\right)\sin(\tau)\right)=0$

we then have

$\displaystyle\omega^{\operatorname{Rot}}([\dot{\widetilde{c}}(\tau),\widetilde{W}(\tau)])$
$\displaystyle=\underbrace{\dot{\widetilde{c}}(\tau)\left(\omega^{\operatorname{Rot}}(\widetilde{W}(\tau))\right)}_{=0}-\underbrace{\widetilde{W}(\tau)\left(\omega^{\operatorname{Rot}}(\dot{\widetilde{c}}(\tau))\right)}_{=0}-d\omega^{\operatorname{Rot}}\left(\dot{\widetilde{c}}(\tau),\widetilde{W}(\tau)\right)$
$\displaystyle=-2\,\big{(}\operatorname{Im}\\!\left(\langle\rho(\tau)\dot{\gamma}(\tau),\widetilde{W}(\tau)\rangle\right)+\underbrace{\operatorname{Im}\\!\left(\langle\dot{\rho}(\tau)\gamma(\tau),\widetilde{W}(\tau)\rangle\right)}_{=0\text{
($\widetilde{W}$ horizontal, $\rho$ and $\dot{\rho}$ real)}}\big{)}$
$\displaystyle=-2\rho(\tau)\,\operatorname{Im}\\!\left(\langle\cos(\tau)\beta^{\prime}-\sin(\tau)\beta,\tilde{\varepsilon}\rangle+\langle\dot{\gamma}(\tau),\operatorname{Im}\\!\left(\langle\beta^{\prime},\tilde{\varepsilon}\rangle\right)\left(\dot{\gamma}(\tau)-\beta^{\prime}\right)\sqrt{\text{-}1}\rangle\right)$
$\displaystyle=-2\rho(\tau)\,\big{(}\cos(\tau)\operatorname{Im}\\!\left(\langle\beta^{\prime},\tilde{\varepsilon}\rangle\right)+\operatorname{Im}\\!\left(\langle\beta^{\prime},\tilde{\varepsilon}\rangle\right)\underbrace{\langle\dot{\gamma}(\tau),\dot{\gamma}(\tau)-\beta^{\prime}\rangle}_{\in\mathds{R}}\big{)}$
$\displaystyle=-2\rho(\tau)\,\big{(}\cos(\tau)\operatorname{Im}\\!\left(\langle\beta^{\prime},\tilde{\varepsilon}\rangle\right)+\operatorname{Im}\\!\left(\langle\beta^{\prime},\tilde{\varepsilon}\rangle\right)-\operatorname{Im}\\!\left(\langle\beta^{\prime},\tilde{\varepsilon}\rangle\right)\cos(\tau)\big{)}$
$\displaystyle=-2\rho(\tau)\,\operatorname{Im}\\!\left(\langle\beta^{\prime},\tilde{\varepsilon}\rangle\right)$

and

$\displaystyle\omega^{\operatorname{Trl}}([\dot{\widetilde{c}}(\tau),\widetilde{W}(\tau)])$
$\displaystyle=\underset{=\langle\\!\text{1},\widetilde{W}\rangle=0}{\dot{\widetilde{c}}(\tau)\left(\underbrace{\omega^{\operatorname{Trl}}\left(\widetilde{W}(\tau)\right)}\right)}-\underset{=\langle\\!\text{1},\dot{\widetilde{c}}\rangle=0}{\widetilde{W}(\tau)\left(\underbrace{\omega^{\operatorname{Trl}}\left(\dot{\widetilde{c}}(\tau)\right)}\right)}-\underbrace{d\omega^{\operatorname{Trl}}\left(\dot{\widetilde{c}}(\tau),\widetilde{W}(\tau)\right)}_{=0}=0,$

where tangent vectors $\dot{\widetilde{c}}(\tau)$ and $\widetilde{W}(\tau)$
are interpreted as directional derivatives. These are the two equations that
remained to show. ∎

###### Proof of Lemma 2.

  1. i)

See, e.g., Lee2018RiemannianManifolds, Proposition B.12 on page 402. This is a
standard result. Note that based on an alternative (yet also common)
definition of the wedge product and, hence, the exterior derivative,
Huckemann2010intrinsicMANOVA and others write
$\omega\left([V,W]\right)=V\left(\omega\left(W\right)\right)-W\left(\omega\left(V\right)\right)-2\,d\omega\left(V,W\right)$
instead. In this case, we also have
$d\omega^{\operatorname{Rot}}(V,W)=\,\operatorname{Im}\\!\left(\langle
V,W\rangle\right)$ in ii) compensating for the different factor in the proof
of Lemma 1.

  2. ii)

Let $\\{e_{r}\\}_{r}$ be an orthonormal $\mathds{C}$-basis of $\mathcal{Y}$ (a
complete orthonormal system existing since $\mathcal{Y}$ is separable) and
$\\{\vartheta^{(r)}(y)\\}_{r}=\langle e_{r},y\rangle$ the corresponding dual
basis. The tangent vectors $\partial_{{\operatorname{Re}},r}\big{|}_{p}\cong
e_{r}$ and
$\partial_{{\operatorname{Im}},r}\big{|}_{p}\cong\sqrt{\text{-}1}\,e_{r}$,
$p\in\mathcal{Y}$ together form an $\mathds{R}$-basis of
$T_{p}\mathcal{Y}^{*}\cong\mathcal{Y}$. The dual 1-forms are given by
$d^{{\operatorname{Re}},r}(V_{p}):=V_{p}\left({\operatorname{Re}}\circ\vartheta^{(r)}\right)\cong{\operatorname{Re}}\circ\vartheta^{(r)}(V_{p})$
and
$d^{{\operatorname{Im}},r}(V_{p}):=V_{p}\left({\operatorname{Im}}\circ\vartheta^{(r)}\right)\cong{\operatorname{Im}}\circ\vartheta^{(r)}(V_{p})$
where we identify tangent vectors either with directional derivatives
$V_{p}(f)=\frac{d}{d\tau}\left(f\circ\operatorname{Exp}_{p}(\tau
V_{p})\right)\big{|}_{\tau=0}$ of functions
$f:\mathcal{M}\rightarrow\mathds{R}$ or with elements of $\mathcal{Y}$, and
the equality follows from $\mathcal{M}=\mathcal{Y}^{*}$, and
${\operatorname{Re}}\circ\vartheta^{(r)}$,
${\operatorname{Im}}\circ\vartheta^{(r)}$ linear. With this given, we have

$\displaystyle\omega^{\operatorname{Rot}}(V_{p})$
$\displaystyle=\operatorname{Im}\\!\left(\langle\sum_{r}\langle e_{r},p\rangle
e_{r},V_{p}\rangle\right)$ (7)
$\displaystyle=\sum_{r}\operatorname{Im}\\!\left(\langle p,e_{r}\rangle\langle
e_{r},V_{p}\rangle\right)$
$\displaystyle=\sum_{r}\operatorname{Re}\\!\left(\langle
e_{r},p\rangle\right)\operatorname{Im}\\!\left(\langle
e_{r},V_{p}\rangle\right)-\operatorname{Im}\\!\left(\langle
e_{r},p\rangle\right)\operatorname{Re}\\!\left(\langle
e_{r},V_{p}\rangle\right)$
$\displaystyle=\sum_{r}{\operatorname{Re}}\circ\vartheta^{(r)}(p)\,d^{{\operatorname{Im}},j}\left(V_{p}\right)-{\operatorname{Im}}\circ\vartheta^{(r)}(p)\,d^{{\operatorname{Re}},j}\left(V_{p}\right)$

and thus, expressing the exterior derivative in terms of wedge products

$\displaystyle d\omega^{\operatorname{Rot}}$
$\displaystyle=\sum_{r}\sum_{l}\partial_{{\operatorname{Re}},r}\left({\operatorname{Re}}\circ\vartheta^{(r)}\right)\
d^{{\operatorname{Re}},l}\wedge
d^{{\operatorname{Im}},r}+\partial_{{\operatorname{Im}},r}\left({\operatorname{Re}}\circ\vartheta^{(r)}\right)\
d^{{\operatorname{Im}},l}\wedge d^{{\operatorname{Im}},r}$
$\displaystyle\quad-\partial_{{\operatorname{Re}},r}\left({\operatorname{Im}}\circ\vartheta^{(r)}\right)\
d^{{\operatorname{Re}},l}\wedge
d^{{\operatorname{Re}},r}-\partial_{{\operatorname{Im}},r}\left({\operatorname{Im}}\circ\vartheta^{(r)}\right)\
d^{{\operatorname{Im}},l}\wedge d^{{\operatorname{Re}},r}$
$\displaystyle=\sum_{r}\sum_{l}d^{{\operatorname{Re}},l}\left(\partial_{{\operatorname{Re}},r}\right)\
d^{{\operatorname{Re}},l}\wedge
d^{{\operatorname{Im}},r}+d^{{\operatorname{Re}},l}\left(\partial_{{\operatorname{Im}},r}\right)d^{{\operatorname{Im}},l}\wedge
d^{{\operatorname{Im}},r}$
$\displaystyle\quad-d^{{\operatorname{Im}},l}\left(\partial_{{\operatorname{Re}},r}\right)d^{{\operatorname{Re}},l}\wedge
d^{{\operatorname{Re}},r}-d^{{\operatorname{Im}},l}\left(\partial_{{\operatorname{Im}},r}\right)d^{{\operatorname{Im}},l}\wedge
d^{{\operatorname{Re}},r}$
$\displaystyle=\sum_{r}d^{{\operatorname{Re}},r}\wedge
d^{{\operatorname{Im}},r}-d^{{\operatorname{Im}},r}\wedge
d^{{\operatorname{Re}},r}$
$\displaystyle=2\sum_{r}d^{{\operatorname{Re}},r}\wedge
d^{{\operatorname{Im}},r}$

which evaluates to

$\displaystyle d\omega^{\operatorname{Rot}}(V,W)$
$\displaystyle=2\sum_{r}\Big{(}d^{{\operatorname{Re}},r}\left(V\right)d^{{\operatorname{Im}},r}\left(W\right)-d^{{\operatorname{Im}},r}\left(V\right)d^{{\operatorname{Re}},r}\left(W\right)\Big{)}$
$\displaystyle=2\,\operatorname{Im}\\!\left(\langle V,W\rangle\right)$

where the last equation follows from a computation analogous to (7).

  3. iii)

By choosing w.l.o.g. $e_{1}=\\!\text{1}$, we obtain

$\displaystyle\omega^{\operatorname{Trl}}(V)=d^{{\operatorname{Re}},1}+\sqrt{\text{-}1}\,d^{{\operatorname{Im}},1}$

which immediately yields $d\langle\\!\text{1},\cdot\rangle=0$, since
$d\,d^{{\operatorname{Re}},1}=d\,d^{{\operatorname{Im}},1}=0$.

∎

###  S.2 Tensor-product factorization

The optimality of the proposed tensor-product factorization follows from the
Eckart-Young-Mirsky theorem (EYM) which can be found, e.g., in [page
139]Gentle2007MatrixAlgebra for matrices and, in more general terms, in [page
111]Hsing2015TheoreticalFoundations for Hilbert-Schmidt operators. In the
following, we present a tensor-product version of EYM designed for our needs.
The optimality of the tensor-product factorization is then illustrated in two
corollaries – first in a theoretical model setting and second for the
empirical decomposition on evaluations which can be practically conducted on
given data. Consider two real vector spaces $\mathcal{B}_{j}$,
$j\in\\{0,1\\}$, with positive semi-definite bilinear forms
$\langle\cdot,\cdot\rangle_{j}:\mathcal{B}_{j}\times\mathcal{B}_{j}\rightarrow\mathds{R}$
inducing semi-norms $\|\cdot\|_{j}$. Assuming $\mathcal{B}_{1}$ to be, in
fact, a function space of functions $f:\mathcal{X}\rightarrow\mathds{R}$ on
some set $\mathcal{X}$, the (vector space) tensor product
$\,\mathcal{B}_{1}\otimes\mathcal{B}_{0}\,$ of $\mathcal{B}_{0}$ and
$\mathcal{B}_{1}$ is the vector space spanned by all $f\otimes
y:\mathcal{X}\rightarrow\mathcal{B}_{0},\mathbf{x}\mapsto f(\mathbf{x})\,y$
with $f\in\mathcal{B}_{1}$ and $y\in\mathcal{B}_{0}$. By linear extension, a
symmetric positive semi-definite bilinear form on
$\,\mathcal{B}_{1}\otimes\mathcal{B}_{0}\,$ is defined by $\langle\,f\otimes
y\,,\,f^{\prime}\otimes
y^{\prime}\,\rangle_{\,\mathcal{B}_{1}\otimes\mathcal{B}_{0}\,}=\langle
f,f^{\prime}\rangle_{1}\,\langle y,y^{\prime}\rangle_{0}$ for all
$f,f^{\prime}\in\mathcal{B}_{1},y,y^{\prime}\in\mathcal{B}_{0}$. It induces a
semi-norm $\|\cdot\|_{\,\mathcal{B}_{1}\otimes\mathcal{B}_{0}\,}$ on the
tensor product space.

######  Theorem 2 (Eckart-Young-Mirsky for finite-dimensional tensor-
products).

Let $\mathcal{B}_{0},\mathcal{B}_{1}$ be semi-normed vector spaces as defined
above and
$h=\sum_{r=1}^{m_{0}}\sum_{l=1}^{m_{1}}\theta^{(r,l)}\,b_{1}^{(l)}\otimes
b_{0}^{(r)}\,\in\,\mathcal{B}_{1}\otimes\mathcal{B}_{0}\,$ expressed as a
finite linear-combination with
$b_{1}^{(1)},\dots,b_{1}^{(m_{1})}\in\mathcal{B}_{1}$,
$b_{0}^{(1)},\dots,b_{0}^{(m_{0})}\in\mathcal{B}_{0}$, and coefficient matrix
$\\{\theta^{(r,l)}\\}_{r,l}=\boldsymbol{\Theta}\in\mathds{R}^{m_{0}\times
m_{1}}$. Then we can optimally decompose
$h=\sum_{r=1}^{m}d^{(r)}\,\xi_{1}^{(r)}\otimes\xi_{0}^{(r)}\,$ with
$m=\min\\{m_{0},m_{1}\\}$, $d^{(1)}\geq\dots\geq d^{(m)}\geq 0$ and
$\langle\xi_{j}^{(r)},\xi_{j}^{(l)}\rangle_{j}=\mathds{1}(r=l)$ for
$\xi_{j}^{(r)}\in\mathcal{B}_{j}$ , in the sense that for any $L\leq m$

$\|h-\sum_{r=1}^{L}d^{(r)}\,\xi_{1}^{(r)}\otimes\xi_{0}^{(r)}\,\|_{\,\mathcal{B}_{1}\otimes\mathcal{B}_{0}\,}\
\leq\
\|h-\sum_{r=1}^{L}d^{(r)}_{\star}\,\xi_{1\star}^{(r)}\otimes\xi_{0\star}^{(r)}\,\|_{\,\mathcal{B}_{1}\otimes\mathcal{B}_{0}\,}$
(8)

for all $d^{(r)}_{\star}\in\mathds{R}$ and
$\xi_{j\star}^{(r)}\in\mathcal{B}_{j}$, $j\in\\{0,1\\}$, $r=1,\dots,L$.
Arranging $\mathbf{D}=\operatorname{diag}(d^{(1)},\dots,d^{(m)})$ and
expressing $\xi_{j}^{(r)}=\sum_{l}u_{j}^{(l,r)}b_{j}^{(l)},r=1,\dots,m,$ with
the coefficient matrices
$\\{u_{j}^{(l,r)}\\}_{l,r}=\mathbf{U}_{j}\in\mathds{R}^{m_{j}\times m}$,
$j\in\\{0,1\\}$, an optimal decomposition is obtained as follows:

  1. i)

If for $j\in\\{0,1\\}$ the Gram matrices $\mathbf{G}_{j}=\\{\langle
b_{j}^{(r)},b_{j}^{(l)}\rangle_{j}\\}_{r,l}$ are the identity
$\mathbf{G}_{j}=\mathbf{I}_{m_{j}}$, the matrices $\mathbf{D}$ and
$\mathbf{U}_{j}$, $j\in\\{0,1\\}$, are directly determined via SVD of the
coefficient matrix
$\boldsymbol{\Theta}=\mathbf{U}_{0}\mathbf{D}\mathbf{U}_{1}^{\top}$.

  2. ii)

In general, there are suitable matrices
$\mathbf{M}_{j}\in\mathds{R}^{\operatorname{rank}\mathbf{G}_{j}\times m_{j}}$,
$j\in\\{0,1\\}$, such that
$\boldsymbol{\Xi}=\mathbf{V}_{0}\mathbf{D}\mathbf{V}_{1}^{\top}$ is the SVD of
the matrix
$\boldsymbol{\Xi}=\mathbf{M}_{0}\boldsymbol{\Theta}\mathbf{M}_{1}^{\top}$ and
$\mathbf{U}_{j}=\mathbf{M}_{j}^{-}\mathbf{V}_{j}$ with generalized inverse
$\mathbf{M}_{j}^{-}=\mathbf{M}_{j}^{\top}(\mathbf{M}_{j}\mathbf{M}_{j}^{\top})^{-1}$.

    1. (a)

In general, a suitable matrix is given by
$\mathbf{M}_{j}=\sqrt{\mathbf{G}_{j}}^{\top}$ with
$\mathbf{G}_{j}=\sqrt{\mathbf{G}_{j}}\sqrt{\mathbf{G}_{j}}^{\top}$ a Cholesky
decomposition.

    2. (b)

If the $b_{j}^{(r)}$ can be identified with vectors
$\mathbf{b}_{j}^{(r)}\in\mathds{R}^{m_{j}^{\prime}}$ of some length
$m_{j}^{\prime}$, arranged as column vectors of a “design matrix”
$\mathbf{B}_{j}\in\mathds{R}^{m_{j}^{\prime}\times m_{j}}$, such that $\langle
b_{j}^{(r)},b_{j}^{(l)}\rangle_{j}=(\mathbf{b}_{j}^{(r)})^{\top}\mathbf{W}_{j}\mathbf{b}_{j}^{(l)}$,
with $r,l=1,\dots,m_{j}$, for a symmetric positive definite weight matrix
$\mathbf{W}_{j}$, we may equivalently set $\mathbf{M}_{j}=\mathbf{R}_{j}$
based on the QR-decomposition
$\sqrt{\mathbf{W}_{j}}^{\top}\mathbf{B}_{j}=\mathbf{Q}_{j}\mathbf{R}_{j}$. In
this case, design matrices $\mathbf{E}_{j}$ of vector representatives
$\boldsymbol{\xi}_{j}^{(r)}$ for the $\xi_{j}^{(r)}$ can, alternatively, be
obtained as
$\mathbf{E}_{j}=\sqrt{\mathbf{W}_{j}}^{-\top}\mathbf{Q}_{j}\mathbf{V}_{j}$
(where $\mathbf{W}_{j}$ is typically diagonal and, hence,
$\sqrt{\mathbf{W}_{j}}^{-\top}$ fast to compute).

###### Proof.

  1. i)

For $j\in\\{0,1\\}$, denote the column vectors of $\mathbf{U}_{j}$ by
$\mathbf{u}_{j}^{(r)}$, $r=1,\dots,m$, and consider the space of $m_{0}\times
m_{1}$ matrices equipped with the inner product
$\langle\boldsymbol{\Theta}_{1},\boldsymbol{\Theta}_{2}\rangle_{F}=\operatorname{tr}{\left(\boldsymbol{\Theta}_{1}^{\top}\boldsymbol{\Theta}_{2}\right)}$,