Introduction

Data in the form of symmetric positive-definite matrices arise in many areas, including computer vision (Rathi et al., 2007;Caseiro et al., 2012), signal processing (Arnaudon et al., 2013;Hua et al., 2017), medical imaging (Fillard et al., 2007;Dryden et al., 2009) and neuroscience (Friston, 2011), among other fields and applications. For instance, they are used to model brain functional connectivity that is often characterized by covariance matrices of blood-oxygen-leveldependent signals (Huettel et al., 2008). In diffusion tensor imaging analysis (Le Bihan, 1991), a 3 × 3 symmetric positive matrix that is computed for each voxel describes the dominant shape of local diffusion of water molecules. While the analysis of symmetric positive-definite matrices as responses in a regression model is our primary emphasis in this paper, our results more generally extend to responses with a Lie group structure that include data on the torus.

The space S + of symmetric positive-definite matrices is a nonlinear metric space and, depending on the metric, forms a Riemannian manifold. Various metrics have been studied (Pigoli et al., 2014); an important criterion for choosing a metric is to avoid the swelling effect. This refers to the phenomenon that the determinant of the Fréchet mean of a set of symmetric positive-definite matrices may be substantially larger than that of any of the constituent matrices. The swelling effect becomes evident in the geodesics connecting two elements of S + (Arsigny et al., 2007), negatively affects specifically the Frobenius metric and various other metrics, and is problematic for many applications. These include diffusion tensor imaging, where a diffusion tensor is represented by a symmetric positive-definite matrix and the determinant quantifies diffusion of the tensor. The swelling effect will then lead to distortions when quantifying diffusion.

The abundance of S + -valued data in many areas stands in contrast with the relative sparsity of work on their statistical analysis, in particular regarding regression with S + -valued responses. Existing work includes Riemannian frameworks to analyse diffusion tensor images with a focus on averages and modes of variation (Pennec et al., 2006;Fletcher & Joshib, 2007) and various versions of nonparametric regression such as spline regression (Barmpoutis et al., 2007), local constant regression (Davis et al., 2010), intrinsic local linear regression (Yuan et al., 2012), wavelet regression (Chau & von Sachs, 2019) and Fréchet regression (Petersen et al., 2019). Various metric, manifold and Lie group structures have been proposed, for example, the trace metric (Lang, 1999), affine-invariant metric (Moakher, 2005;Pennec et al., 2006;Fletcher & Joshib, 2007), log-Euclidean metric (Arsigny et al., 2007), log-Cholesky metric (Lin, 2019), scalingrotation distance (Jung et al., 2015) and Procrustes distance (Dryden et al., 2009). As S + is a Riemannian manifold and more generally a metric space, regression techniques developed for general Riemannian manifolds (e.g., Pelletier, 2006;Shi et al., 2009;Davis et al., 2010;Steinke et al., 2010;Fletcher, 2013;Hinkle et al., 2014;Cornea et al., 2017, among many others) and metric spaces (Hein, 2009;Chen & Müller, 2022;Lin & Müller, 2021) also apply to S + .

Additive regression originating with Stone (1985) is known to be an efficient way of avoiding the well-known curse of the dimensionality problem that one faces in nonparametric regression when the dimension of the covariate vector increases, but so far has been by and large limited to the case of real-valued and functional responses. Examples for additive regression approaches for real-valued responses include the original work on smooth back fitting (Mammen et al., 1999) and extensions to generalized additive models (Yu et al., 2008), additive quantile models (Lee et al., 2010), generalized varying coefficient models (Lee et al., 2012), errors in variables (Han & Park, 2018) and functional or distributional responses (Scheipl et al., 2015;Park et al., 2018;Han et al., 2020;Jeon & Park, 2020).

In this paper, we develop an additive regression model for responses residing in a Lie group that includes the space S + as the primary example, with an extension to general Riemannian manifolds. This paper contains three major contributions. First, to the best of our knowledge, this is the first paper to study additive regression for S + -valued responses, counteracting the curse of dimensionality while maintaining a high degree of flexibility. Previous studies for modelling S + -valued responses focused on unstructured nonparametric regression such as local constant/polynomial regression that are well known to be subject to the curse of dimensionality when there are many predictors.

Second, by focusing on Abelian Lie group-valued responses that include data on the torus (Eltzner et al., 2018), as well as S + -valued responses, we propose a novel intrinsic group additive regression model that directly exploits the Abelian group structure of the responses. This sets our work apart, as previously only the general manifold structure of responses was considered in the few existing nonadditive regression approaches. To the best of our knowledge, this is the first work to connect additive regression with Lie groups. This connection extends the reach of additive models in a very natural way.

Third, we show that this group additive model can be transformed into an additive model on tangent spaces by utilizing the Riemannian logarithmic map, which paves the way for extending the additive model to the case of responses that lie on more general manifolds. While throughout we showcase the proposed approaches for the space S + of symmetric positive-definite matrices with suitable metrics, our results are by no means limited to this type of response, and are applicable to additive regression modelling for a much larger class of manifold-valued responses.

Preliminaries on differential geometry

We compile here some basic notions for Riemannian manifolds and Lie groups, referring readers to § S.1 of Shao et al. (2022) for a self-contained note on basic concepts of Riemannian geometry and to the text by Lee (2018) for a comprehensive treatment. In the following, L denotes a simply connected smooth manifold modelled on a D-dimensional Hilbert space.

The tangent space at y ∈ L, denoted by T y L, is a linear space consisting of velocity vectors α (0), where α : (-1, 1) → L represents a differentiable curve passing through y, i.e., α(0) = y. Each tangent space T y L is endowed with an inner product g y that varies smoothly with y and is a D-dimensional Hilbert space with the induced norm denoted by • y . The inner products {g y : y ∈ L} are collectively denoted by g, referred to as the Riemannian metric of L that also induces a distance d on L.

A geodesic γ is a constant-speed curve defined on [0, ∞) such that, for each t ∈ [0, ∞), the segment γ ([t, t + ]) is the shortest path connecting γ (t) and γ (t + ) for all sufficiently small > 0. The Riemannian exponential map Exp y at y ∈ L is a function mapping T y L into L and defined by Exp y (u) = γ (1) with γ (0) = y and γ (0) = u ∈ T y L. Conversely, γ y,u (t) = Exp y (tu) is a geodesic starting at y and with direction u.

For a tangent vector u ∈ T y L, the cut time c u is the positive number such that γ y,u ([0, c u ]) is a shortest path connecting γ y,u (0) and γ y,u (c u ), but γ y,u ([0, c u + ]) is not a shortest path for any > 0. It turns out that the exponential map Exp y (tu) for u ∈ T y L is invertible before the cut time c u . Formally, the inverse of Exp y , denoted by Log y and called the Riemannian logarithmic map at y, can be defined by Log y z = tu for z ∈ E y := {Exp y (tu

Let C ∞ (L) denote the collection of smooth real-valued functions defined on L. For a smooth function f ∈ C ∞ (L) and a tangent vector v ∈ T y L, the covariant derivative of f at y along the direction v, denoted by ∇ v f , is defined by

where γ :

A smooth vector field U is a smooth function defined on L such that U (y) ∈ T y L for all y ∈ L. The covariant derivative measures how fast a map changes along a direction and is defined for smooth vector fields as follows. Let (L) denote the collection of smooth vector fields on L. A connection is a map ∇ : (L)× (L) → (L), with (V , U ) → ∇ V U , that satisfies the following properties:

For f ∈ C ∞ (L) and a smooth vector field U , fU denotes a smooth vector field defined by ( fU )(y) = f (y)U (y) for all y ∈ L, and ∇ v f is a smooth real-valued function defined by

The value of ∇ V U at y depends on V only through its value at y (Proposition 4.5 of Lee, 2018), which makes the expression

, then we say that the connection ∇ is torsion-free. For a Riemannian manifold, there exists a unique connection that is both torsionfree and compatible with the Riemannian metric. It is called the Levi-Civita connection with Levi-Civita covariant derivative as its induced covariant derivative.

Curvature quantifies the degree of deviation from being flat. Define the map

The value of R(U , V , W ) at y depends only on the values of U , V , W at y, and therefore we can write R (u, v, w) for tangent vectors u, v, w at the same point. The sectional curvature at y ∈ L is a real-valued function on T y L × T y L defined for u, v ∈ T y L by

A Hadamard manifold is a complete and simply connected Riemannian manifold that has everywhere nonpositive sectional curvature and is thus a Hadamard space.

Given a curve γ (t) on L, t ∈ I for a real interval I , a vector field U along γ is a smooth map defined on I such that U (t) ∈ T γ (t) L. We say that U is parallel along γ if ∇ γ (t) U = 0 for all t ∈ I . In this paper, we primarily focus on parallel vector fields along geodesics. Let γ : [0, 1] → L be a geodesic connecting y and z, and let U be a parallel vector field along γ such that U (0) = u and U (1) = v. Then we say that v is the parallel transport of u along γ , denoted by τ y,z u = v. Parallel transport can be used as an intrinsic mechanism to compare tangent vectors residing at different points, e.g., via parallelly transporting the tangent vectors to the tangent space at a fixed point on L, where they can be easily compared.

A Lie algebra is a vector space g endowed with an alternating binary operation [•, •] : g×g → g satisfying the following axioms:

When L is a group such that the group operation '⊕' and inverse ι : y → y -1 are smooth, (L, ⊕) is a Lie group. We say that a vector field U on a Lie group (L, ⊕) is left invariant if U (y ⊕ z) = (D z L y ){U (z)} for all y, z ∈ L, where L y : z → y ⊕ z is the left translation induced by y and D z L y is the differential of L y at z. Right-invariant vector fields are defined in a similar fashion via right translations. One can show that any left-invariant vector field U is fully specified by its value U (e) at the group identity element e. Consequently, the collection of left-invariant vector fields, denoted by h, can be identified with the tangent space T e L. In addition, the space h gives rise to a Lie algebra, as follows. For two smooth vector fields U , V on the Lie group L, the Lie bracket [U , V ] is the vector field determined by

It can be shown that the vector space h endowed with the Lie bracket is a Lie algebra.

)v} for all y, z ∈ L and u, v ∈ T z L. Right-invariant metrics can be defined in a similar fashion. A metric is bi-invariant if it is both left invariant and right invariant.

The Lie exponential map, denoted by exp, that maps g into L is defined by exp(u) = γ (1), where γ : R → L is the unique one-parameter subgroup such that γ (0) = u ∈ g. Its inverse, if it exists, is called the Lie logarithmic map and denoted by log. When g is bi-invariant, exp = Exp e , that is, the Riemannian exponential map at the identity element coincides with the Lie exponential map.

For random elements Y ∈ L, for a Lie group L, the Fréchet function is F(y) = Ed 2 (y, Y ), where d is the Riemannian distance function induced by the metric on L. If L is a Hadamard manifold and F(y) < ∞ for some y ∈ L, and hence F(y) < ∞ for all y ∈ L according to the triangle inequality, the minimizer of F(y) exists and is unique (Sturm, 2003). It is known as the Fréchet mean and is denoted by E ⊕ Y .

Additive models for Lie groups

In this section we develop an additive regression model for Lie group-valued responses. The main motivating example of Lie groups is the space S + m of m × m symmetric positive matrices due to their ubiquity and practical relevance.

Example 1 (S + m with the log-Euclidean metric). The space S + m is a smooth submanifold of R m×m . Its tangent spaces are S m , the collection of m×m symmetric matrices. The matrix logarithm is a smooth bijective map between S + m and S m , and thus can be used to transfer the canonical Euclidean structure of S m to S + m . Recall that, for a symmetric matrix S, exp(S) = I m + ∞ j=1 S j /j! is a symmetric positive-definite matrix. The inverse of exp, denoted by log, exists and is called the matrix logarithmic map. Both exp and log are smooth maps between S + m and S m . The operation '⊕'defined as P 1 ⊕P 2 = exp{log(P 1 )+log(P 2 )} for P 1 , P 2 ∈ S + m turns S + m into an Abelian group. The canonical Riemannian metric on S m is trace(S 1 S 2 ) for S 1 , S 2 ∈ S m and can be transferred to a Riemannian metric on S + m given by g

, where D P log denotes the differential of the log at P. It turns out that g is a bi-invariant metric on (S + m , ⊕) that is isomorphic to the group of S m with the usual matrix addition as the group operation (Proposition 3.4 of Arsigny et al., 2007). This metric, termed the log-Euclidean metric and specifically designed to eliminate the swelling effect (Arsigny et al., 2007), turns S + m into a Hadamard manifold and a bi-invariant Abelian Lie group. This metric is invariant to permutation, but is not affine invariant.

Example 2 (S + m with the log-Cholesky metric). The log-Cholesky metric (Lin, 2019) utilizes the Cholesky decomposition to transfer a Lie group structure on the space of lower triangular matrices of positive diagonals to the space of symmetric positive-definitive matrices. Specifically, let lt(m) be the space of m×m lower triangular matrices and lt + (m) ⊂ lt(m) the subspace, such that L ∈ lt + (m) if all diagonal elements of L are positive. One can show that lt + (m) is a smooth submanifold of lt(m) and its tangent spaces are identified with lt(m). For a fixed L ∈ lt + (m), we define a Riemannian metric g on lt

where A ij denotes the element of A in the ith row and jth column. It is an Abelian Lie group with the operation defined by

is the diagonal part of L, that is, a diagonal matrix whose diagonals are equal to the respective diagonals of L. One can show that g is a bi-invariant metric for the Lie group lt + (m) with the group operation . It is well known that a symmetric positive-definite matrix P is associated with a unique matrix L in lt + (m) such that LL T = P. We refer to L as the Cholesky factor of P. For U , V ∈ T P S + m = S m , we define the metric

, where (S) 1/2 = L(S)+D(S)/2 for a matrix S. We also turn S + m into an Abelian Lie group with the operator ⊕ such that

where L 1 and L 2 are the Cholesky factors of P 1 and P 2 , respectively. The metric g, also specifically designed to eliminate the swelling effect (Lin, 2019), is a bi-invariant metric of the Lie group (S + m , ⊕) and turns S + m into a Hadamard manifold. In our experience, the log-Cholesky metric is computationally more efficient, but not invariant to permutation.

The torus serves as an example for a relevant Lie group that is not related to the space S + m and shows that our approach is not limited to symmetric positive matrices.

It is an Abelian Lie group when the group operation is the standard multiplication of complex numbers. In addition, S is a Riemannian submanifold of R 2 . A k-torus T k is the direct product of k copies of S, and is therefore also an Abelian Lie group. Its Lie algebra is g =

all u, v ∈ g. The Riemannian metric on T k is the product Riemannian metric of S, and this metric is bi-invariant. Data on tori occur in the statistical analysis of wind directions and RNA data (Eltzner et al., 2018).

Given scalar variables X 1 ∈ X 1 , . . . , X q ∈ X q , which are predictors that are paired with a Lie group-valued response Y and where X j ⊂ R, j = 1, . . . q, are compact domains, the proposed Lie group additive model is as follows, where we make use of the Lie group operation ⊕:

(1)

the group identity element e. For identifiability, we require that E ⊕ {w k (X k )} = e, k = 1, . . . , q. In addition, we assume existence and uniqueness of μ; see § 4 for more discussions. This model generalizes the additive model for Euclidean responses to L- valued responses. It includes the effect of the additive component functions on the mean response and of noise on the responses, neither of which can be additively modelled in the absence of a linear structure of the responses. The task at hand is to estimate the unknown parameter μ and the component functions w 1 , . . . , w q , given a sample of independently and identically distributed observations of size n. To overcome the challenge of the absence of a linear structure in L, the following result proves essential.

Proposition 1. If (L, ⊕) is an Abelian Lie group endowed with a bi-invariant metric g that turns L into a Hadamard manifold, then ( 1) is equivalent to

(2)

In (2), we transform the left-hand side of (1) by the Riemannian log map Log μ , so that ELog μ Y = 0, and the right-hand side of (1) by the Lie log map log, whereupon w 1 (X 1 ) ⊕ • • • ⊕ w q (X q ) is decomposed into additive components in the vector space T μ L, so that the standard smooth back-fitting algorithm can be adopted to estimate w 1 , . . . , w q . Specifically, let f k (X k ) = τ e,μ logw k (X k ) and ε = τ e,μ logζ . Then, according to Proposition 1, one may rewrite model ( 1 1). These considerations motivate us to estimate the component functions w k through estimation of the f k , as follows.

Step 1. Compute the sample Fréchet mean μ. Closed-form expressions of μ are available for special cases of L, including L = S + m with the log-Cholesky or log-Euclidean metric; see the Supplementary Material for details. Alternative numerical algorithms (e.g., Yang, 2007) are also available.

Step 2. Compute Log μY i . There are closed-form expressions available for L = S + m with the log-Cholesky or log-Euclidean metric; see the Supplementary Material. Numerical methods as described in § 5.3 of A. Brun's 2007 PhD thesis from the University of Linköping are also available.

Step 3. Solve the system of integral equations

where K h j is a kernel function with X j K h j (u, v) du = 1 for all v ∈ X j ; see Jeon & Park (2020).

Note that n -1 n i=1 Log μY i = 0 since μ is the sample Fréchet mean.

Step 4. Finally, estimate w k (x k ) by ŵk (x k ) = exp{τ μ,e fk (x k )}.

Step 3 is a multivariate version of the standard smooth back-fitting system of equations (Mammen et al., 1999). Since the tangent space T μL is also a Hilbert space, the above smooth back-fitting system of equations can be interpreted from a Bochner integral perspective; see Jeon & Park (2020), where also the empirical selection of bandwidths h k is discussed.

Remark 1. By Theorem 3.6 of Bröcker & tom Dieck (1985), any Abelian Lie group is isomorphic to T k × R s , where T k denotes a k-torus defined in Example 3. However, our model, estimation and theory do not rest on this isomorphism, as it does not lead to a natural extension to general Riemannian manifolds. Instead, in the above development we connect the additive model (1) on Lie groups to an additive model (2) on the tangent space at the Fréchet mean via Proposition 1. This allows an immediate extension to general Riemannian manifolds. Since neither our estimation method nor the theory presented in § 5 relies on the isomorphism to T k ×R s , our methods and theory cover both Lie groups and general Riemannian manifolds.

Extension to Riemannian manifolds

When S + m is endowed with the affine-invariant metric (Moakher, 2005;Pennec et al., 2006;Fletcher & Joshib, 2007), it is generally not an Abelian Lie group with a bi-invariant metric, and then model (1) depends on the order of operations and is thus not additive. Specifically, Proposition 1 ceases to hold. However, model (2) remains additive in all cases, suggesting a natural extension to accommodate other metrics and general Riemannian manifolds. For a Riemannian manifold M, which is not necessarily an Abelian Lie group, consider

where μ = E ⊕ Y , ε ∈ T μ M is centred and of finite variance, and f 1 , . . . , f q : R → T μ M are unknown functions to be estimated. Model (3) includes (1) as a special case by setting f k (x) = τ e,μ logw k (x) and ε = τ e,μ logζ according to Proposition 1, transforms the response into the tangent space T μ M and is thus applicable to general Riemannian manifolds, at the expense of interpretability in the original space, as per the group operation in model (1).

For general Riemannian manifolds that might feature positive sectional curvature, the Fréchet mean may not exist and therefore additional conditions are required for model (3). Specifically, if (M, g) is a general Riemannian manifold and Y is a random element on M, we assume the following.

Assumption 1. The minimizer of the Fréchet function F(•) = Ed 2 (•, Y ) exists and is unique. This is automatically satisfied when M is a Hadamard manifold, such as the space S + m with either of the log-Cholesky and log-Euclidean metrics. For other manifolds, we refer the reader to Bhattacharya & Patrangenaru (2003) and Afsari (2011) for conditions that imply Assumption 1.

For a nonempty subset A ⊂ M, let d(y, A) = inf {d(y, z) : z ∈ A} denote the distance between y and the set A. For a positive real number , set A = {y : d(y, A) < } and A -= M\(M\A) . When A = ∅, set A = ∅. The following assumption is only needed for the case where M is not a Hadamard manifold.

Assumption 2. It holds that Pr{Y ∈ E -0 μ } = 1 for some 0 > 0, where E μ is defined in § 2.

If Assumptions 1 and 2 are satisfied, the proposed manifold additive model (3) remains well defined, and the first three steps of the estimation method described in the previous section are still valid and can be employed to estimate f 1 , . . . , f q , with L replaced by M.

Theory

We first establish convergence rates and asymptotic normality of the estimators for the mean and the component functions for general manifolds in the manifold additive model (3), and then provide additional details for the space S + m endowed with either the log-Cholesky metric or the log-Euclidean metric. We consider a manifold M that satisfies one of the following conditions. Condition 1. Manifold M is a finite-dimensional Hadamard manifold that has sectional curvature bounded from below by c 0 0. Condition 2. Manifold M is a complete compact Riemannian manifold.

An example of a manifold satisfying Condition 1 is S + m endowed with the log-Cholesky metric, log-Euclidean metric or affine-invariant metric, while the unit sphere that is used to model compositional data (Dai & Müller, 2018) serves as an example of a manifold that satisfies Condition 2.

To establish the convergence rate of μ, we also make the following assumptions.

Assumption 3. Manifold M satisfies at least one of Conditions 1 and 2.

Assumption 4. For some constant c 2 > 0, F(y) -F(μ) c 2 d 2 (y, μ) when d(y, μ) is sufficiently small. Assumption 5. For some constant c 3 > 0, for all y, z ∈ M, the linear operator H y,z :

Assumption 4 is satisfied for Hadamard manifolds with c 2 = 1 according to Lemma S.7 of Lin & Müller (2021), and the CAT(0) inequality that holds for Hadamard manifolds (Chapter II.1 of Bridson & Häfliger, 1999). The assumption also holds for some manifolds of positive curvature when data concentrate on a small region; see Example 4 of Lin & Müller (2021). The operator H y,z in the technical Assumption 5 is the Hessian of the squared distance function d; see also equation (5.4) of Kendall & Le (2011). Assumption 5 is superfluous if the manifold M is compact and is satisfied by manifolds of zero curvature. It can also be replaced with a uniform moment condition on the operator norm of H z,Y over all z in a small local neighbourhood of μ.

We then obtain a parametric convergence rate for the estimate μ of the Fréchet mean μ.

Proposition 2. Assume that Assumptions 1, 3 and 4 hold and that Y is of second order. Then d( μ, μ) = O P (n -1/2 ).

To obtain convergence rates of the estimated component functions, we require some additional conditions that are standard in the literature on additive regression. Condition 3. The kernel function K is positive, symmetric, Lipschitz continuous and supported on [-1, 1] with K(x) dx = 1.

Condition 4. The bandwidths h k satisfy n 1/5 h k → α k > 0.

Condition 5. The joint density p of X 1 , . . . , X q is bounded away from zero and infinity on X ≡ X 1 × • • • × X q . The densities p kj are continuously differentiable for 1 j | = k q.

Condition 6. The additive functions f k are twice continuously (Fréchet) differentiable.

Without loss of generality, assume that X k = [0, 1] for all k and let

The moment condition on ε in the following theorem is required to control the effect of the error of μ as an estimator of μ on the discrepancies of Log μY i from Log μ Y i after parallel transport; see the Supplementary Material. It is a mild requirement and is satisfied for example when the manifold is compact or ε μ follows a subexponential distribution.

Theorem 1. Under Assumptions 1-5 and Conditions 3-6, if E ε α μ < ∞ for some α 10 and E( ε 2 μ | X j = •) are bounded on X j , respectively, for 1 j q, it holds that

where τ μ,μ is the parallel transport operator along geodesics.

The following corollary is an immediate consequence of Theorem 1.

Corollary 1. Under the conditions of Theorem 1, for S + m endowed with either the log-Cholesky metric or the log-Euclidean metric,

To derive the asymptotic distribution of fk , we define

where the α k are the constants in Condition 4. Let ( 1 , . . . , q ) be a solution of the system of equations

satisfying the constraints

We assume that the following conditions hold.

Condition 8. The ∂p/∂x k , k = 1, . . . , q, exist and are bounded on X = q k=1 X k .

Condition 7 is superfluous if the random noise ε is independent of the predictors X 1 , . . . , X q . Let N μ (x) be the product measure N {θ 1 (x 1 ), 1 (x 1 )} × • • • × N {θ q (x q ), q (x q )} on (T μ M) q , where N (θ, ) denotes a Gaussian measure on T μ M with mean vector θ and covariance operator . For a set A, let Int(A) denote the interior of A.

Theorem 2. Suppose that Assumptions 1-5 and Conditions 3-8 hold, that E ε α μ < ∞ for some α > 10 and that there exists

For M = S + m equipped with either the log-Cholesky metric or the log-Euclidean metric, the above asymptotic normality can be formulated on the Lie algebra g. To this end, assume that SPD 1 , . . . , SPD q and SPD 1 , . . . , SPD q are defined by ( 4)-( 9), with C k (x) and f k replaced by E{logζ ⊗ logζ | X k = x} and ψ k := logw k , respectively. Also, let c

The following corollary is an immediate consequence of Theorem 2, by noting that the manifold S + m when equipped with the log-Cholesky metric or the log-Euclidean metric satisfies Assumptions 1-4 when the second moment of the random noise ζ is finite.

Corollary 2. Assume that Conditions 3-8 hold and that E logζ α μ < ∞ for some α > 10. Furthermore, assume that there exists α > 5/2 such that E( logζ α e | X k = •) are bounded on X k for all 1 k q. For S + m endowed with either the log-Cholesky metric or the log-Euclidean metric, for x = (x 1 , . . .

These results generalize the work of Jeon & Park (2020). One of the major technical challenges that is addressed in Lemmas S1 and S2 in the Supplementary Material is to uniformly quantify the discrepancy between log μY i and log μ Y i due to μ | = μ, as the asymptotic behaviour of this discrepancy plays an important role in the theoretical analysis.

Simulations

To illustrate the numerical performance of the proposed manifold additive model estimators, we conducted simulations for L = S + m endowed with the log-Cholesky and log-Euclidean metrics, respectively. We consider two matrix dimensions, m = 3 and m = 10. We set X k = [0, 1] for k = 1, . . . , q. The predictors X 1 , . . . , X q are set to X 1 = B 2/3,2/3 (T 1 ) and X k = (T k ) for k = 2, . . . , q, where B 2/3,2/3 denotes the cumulative distribution function of the beta distribution with both shape parameters equal to 2/3, denotes the cumulative distribution function of the standard Gaussian distribution, and (T 1 , . . . , T q ) is sampled from the centred q-dimensional Gaussian distribution with covariance matrix whose (j, k)th entry is 1 if j = k and 1/5 if j | = k. Consequently, X 1 , . . . , X q are correlated and X 1 has a nonuniform distribution on [0, 1]. The mean μ is an m × m matrix whose (j, k)th entry is 2 -|j-k| . We then generate the response variable Y by Y = μ ⊕ w(X 1 , . . . , X q ) ⊕ ζ , where w(X 1 , . . . , X q ) = expτ μ,e f (X 1 , . . . , X q ) with three settings for f :

, where f k is defined as in setting I, while f 12 (x 1 , x 2 ) is an m × m matrix whose (j, l)th entry is g(x 1 ; j, l, q)g(x 2 , j, l, q);

The random noise ζ is generated according to logζ = p j=1 Z j v j , where p = m(m + 1)/2 is the dimension of T e S + m , Z 1 , . . . , Z p are independently sampled from N (0, σ 2 ) and v 1 , . . . , v p form an orthonormal basis of the tangent space T e S + m . The signal-to-noise ratio is measured by snr = E logw(X 1 , . . . , X q ) 2 e /E logζ 2 e . We tweak the value of the parameter σ 2 to cover two settings for the signal-to-noise ratio, namely, snr = 2 and snr = 4. Model I is additive, while models II and III are not additive. In particular, model III has no additive components and thus represents the most challenging scenario for the proposed additive regression. We consider q = 3 and q = 4 to probe the effect of the dimensionality of the predictor vector and study sample sizes n = 50, 100, 200.

The quality of the estimation is measured by the root-mean-square error

As a comparison method for the proposed manifold additive model, we also implemented the intrinsic local polynomial regression proposed in Yuan et al. (2012), which is a fully nonparametric approach. In addition, we implemented the following baseline method proposed by an anonymous reviewer, which provides a simple and straightforward approach based on Cholesky decomposition: each symmetric positive-definite matrix is represented by its Cholesky factor, and then a standard multivariate additive model is applied for the Cholesky factor. This simple method may be subject to the swelling effect as it does not use a swelling-free geometry on S + m ; see Example 1 of Lin (2019) for an illustration. Each simulation setting was repeated 100 times. The Monte Carlo prediction root-mean-square error and its standard error are shown in Table 1 for the log-Cholesky metric and in Table 2 for the log-Euclidean metric for m = 3; the results for m = 10 are similar and can be found in the Supplementary Material, where we also graphically illustrate the estimation quality of the proposed method for snr = 4, q = 3, m = 3 and the log-Euclidean metric. When the model is correctly specified as in setting I, the proposed model outperforms the other two methods by a significant margin. When the underlying model is not fully additive, but contains some additive components, such as the model in setting II, our proposed approach still outperforms the other two methods. When the true model has no additive components, such as the model in setting III, all methods have comparable performance when data are sampled from the log-Cholesky metric, while our method and the Cholesky-based method exhibit advantages for the log-Euclidean metric. Also, our method enjoys smaller prediction root mean squared error, compared to the baseline Cholesky method in settings I and II, while both have similar prediction root-mean-square error in setting III. As the true model is unknown in reality and our method is competitive even when there are no additive components, our proposed manifold additive model is overall preferable in applications. In terms of computational efficiency, we found that the log-Cholesky metric is computationally faster than the log-Euclidean metric. addition, the proposed method for the log-Cholesky metric and the baseline method have comparable computational efficiency, while intrinsic local polynomial regression is fastest in computation; see the Supplementary Material.

Application to diffusion tensor imaging

We applied the proposed additive model to diffusion tensors obtained from the Alzheimer's Disease Neuroimaging Initiative, ADNI. Diffusion tensors are 3 × 3 symmetric positive-definite matrices that characterize diffusion of water molecules in tissues, and convey rich information about brain tissues with important applications in tractography. They are utilized to aid in the diagnosis of brain-related diseases. In statistical modelling, diffusion tensors are typically considered to be random elements in the space S + 3 (Fillard et al., 2005;Arsigny et al., 2006;Lenglet et al., 2006;Pennec, 2006;Fletcher & Joshib, 2007;Dryden et al., 2009;Zhu et al., 2009;Zhou et al., 2016;Pennec, 2020). A traditional Euclidean framework for diffusion tensors suffers from significant swelling effects that undesirably inflate the diffusion tensors (Arsigny et al., 2007) and impede their interpretation. In our analysis we use diffusion tensors as responses under the log-Euclidean metric, which is designed to eliminate the swelling effect and relate these responses to several covariates.

We focus on the hippocampus, which plays a central role inAlzheimer's disease (Lindberg et al., 2012). In the ADNI study, brain images and assessment of memory, executive functioning and language ability were obtained for participating subjects. For each raw diffusion tensor image, a standard pre-processing protocol including denoising, eddy current and motion correction, skull stripping, bias correction and normalization was applied and then diffusion tensors for each hippocampal voxel were extracted, followed by computing their log-Euclidean mean. For each raw image, this resulted in an average diffusion tensor representing the typical hippocampal diffusion of the corresponding subject at the time of visit. To study the relation between the average hippocampal diffusion tensor and memory, executive functioning and language ability of the subject, we utilized the neuropsychological summary scores available from ADNI (Gibbons et al., 2012). We only included data from the first visit of subjects who were diagnosed as having either early mild cognitive impairment, mild cognitive impairment, late mild cognitive impairment or Alzheimer's disease. Records with missing values were excluded, under the assumption that missingness occurs at random. This resulted in 220 data tuples of the form (Y , X 1 , X 2 , X 3 ), where Y is the average diffusion tensor, which serves as response, while the predictors X 1 , X 2 , X 3 are standardized scores for memory, executive functioning and language ability, respectively.A subset of the data is presented in Fig. 1.

The estimated component functions ŵ1 (x 1 ), ŵ2 (x 2 ), ŵ3 (x 3 ) and their individual effect on the diffusion tensors are depicted in Fig. 2. For interpretation, we denote the coordinate system of X 1 X 2 X 3 Y -0.13 -0.85 0.20 -0.08 1.41 0.01 0.50 0.13 -0.89 -0.28 1.33 1.29 -0.64 -0.45 -0.15 -0.55 -1.67 -0.73 -1.33 1.13 0.34 0.54 1.60 -0.07 0.48 1.31 1.36 -1.02 -0.09 -0.22 0.33 -1.17 -0.97 -0.88 0.79 0.25 0.60 0.19 0.11 -0.34 0.62 0.99 -0.33 0.86 0.08 -6.93 x 10 -5 -4.86 x 10 -5 -4.

16 x 10 -5 -3.82 x 10 -5 -2.95 x 10 -5 -2.03 x 10 -5 -4.82 x 10 -6 -1.55 x 10 -6 6.07 x 10 -9 1.11 x 10 -6 2.06 x 10 -6 3.82 x 10 -6 6.64 x 10 -6 8.49 x 10 -6 9.53 x 10 -6 1.39 x 10 -5 8.16 x 10 -5 8.50 x 10 -4 9.45 x 10 -4 1.01 x 10 -3 1.04 x 10 -3 1.06 x 10 -3 1.11 x 10 -3 1.18 x 10 -3 1.57 x 10 -3 Fig. 1. A random sample of 15 observations from the data. colour key for ŵ1 , ŵ 2 and ŵ3 -0.017 345 -0.007 392 -0.005 043 -0.003 576 -0.002 424 -0.001 524 -0.000 842 -0.000 245 0.000 139 0.000 528 0.000 995 0.001 408 0.001 934 0.003 368 0.010 022 0.019 368 0.072 371 0.946 717 0.965 864 0.984 860 1.004 969 1.032 556 1.068 471 1.123 330 1.501 364 X1 ŵ1 -2.32 -1.95 -1.59 -1.22 -0.86 -0.49 -0.13 0.24 0.60 0.97 1.33 1.70 2.06 2.43 2.79 ⊕ μ = μ ⊕ ŵ1 X2 ŵ2 -2.92 -2.52 -2.12 -1.73 -1.33 -0.93 -0.54 -0.14 0.26 0.65 1.05 1.45 1.84 2.24 2.64 ⊕ μ = μ ⊕ ŵ2 X3 ŵ3 -4.42 -3.94 -3.47 -2.99 -2.52 -2.04 -1.57 -1.09 -0.62 -0.14 0.33 0.81 1.28 1.76 2.23 R 3 that was adopted to record the diffusion tensors by {e 1 , e 2 , e 3 }, so that the matrices presented in Figs. 1 and 2 indicate the coefficients of the corresponding diffusion tensors in this coordinate system. We find that the component functions have distinct effects on the outcome. Considering, for example, how ŵ3 , which encodes the fit for language ability as predictor and acts on the average diffusion tensor μ by the group operation ⊕, affects the diffusion along the directions e 1 and e 2 , we find that they become increasingly negatively correlated as the standardized language ability drops, indicating nonlinear changes in the diffusion that are associated with changing language ability. As a second example, consider ŵ1 , which encodes the effect of memory. A negative correlation between e 1 and e 2 at high memory levels changes into a positive correlation as memory levels drop in a nonlinear fashion, suggesting distinctive relationships between diffusion patterns and performance scores. This leads to the question of whether the covariates X k for k = 1, 2, 3 are significantly related to the spatially averaged hippocampal diffusion tensor and motivates testing the global null hypothesis H 0 : w k (x k ) = I 3 for all x k , (10)

where I 3 denotes the 3 × 3 identity matrix that is also the group identity of the Lie group (S + m , ⊕) in the log-Euclidean framework, that is, P ⊕ I 3 = P for all P ∈ S + 3 . Corollary 2 can be employed for testing this hypothesis: for a set H k of values of X 1 , a test based on the asymptotic normality of Corollary 2 can be implemented to obtain the p-value for testing the local null hypothesis H 0 : w k (x k ) = I 3 for each x k ∈ H k , followed by adjustment for multiple comparisons, for example, by the Benjamini-Hochberg method. Here, a natural choice of H k is the set of the observed values for X k in the data. The global null hypothesis (10) is then rejected if at least one adjusted p-value is less than the nominal level α. Implementing this approach and applying it to the ADNI data leads to rejecting the null hypothesis (10) at the level α = 0.05 for all k = 1, 2, 3, with the minimal corrected p-values 1.028 × 10 -9 , 1.80 × 10 -5 and < 10 -10 , respectively. This suggests that there are indeed associations between the spatially averaged hippocampal diffusion tensor and memory, executive functioning and language ability.

Discussion

There are at least three potential extensions of the proposed methods and theory. First, in the data application we consider the spatially averaged hippocampal diffusion tensor. Although significantly reducing the data noise level, the average may conceal some spatial structure of interest within the hippocampus. One way to address this problem is to view all hippocampal diffusion tensors derived from an image as an S + 3 -valued function Y : s → Y (s) ∈ S + 3 with s ranging over all hippocampal voxels. This functional perspective enables one to borrow information from neighbouring voxels to counter the high-level data noise. However, it is rather challenging to develop an additive model for S + 3 -valued functions (see Dubey & Müller, 2020). Second, we focus only on the first visit of each subject, and thus do not utilize all available data. To analyse the data of repeated visits that by default are correlated, the theory needs to be extended to account for such correlation, which seems nontrivial, especially when the number of visits may grow with the sample size. Third, when testing hypothesis (10), we perform multiple local tests by using the pointwise asymptotic normality of Corollary 2 and then make corrections for multiple comparisons. This approach, although sufficient for our data application, will be too conservative in general. An alternative method is to develop an asymptotic normality result for the random process { ŵk (x) : x ∈ X k }, which may lead to more powerful one-step tests. These extensions are left for future studies.

data, but did not participate in the analysis or writing of this report. A complete listing of ADNI investigators and funding for ADNI can be found at http://adni.loni.usc.edu/ wp-content/uploads/how_to_apply/ADNI_Acknowledgement_List.