Binary Regression and Classification with Covariates in Metric Spaces

TL;DR

This paper introduces a novel regression and classification framework for binary data with covariates in general metric spaces, providing theoretical bounds and demonstrating practical effectiveness through simulations and fMRI data analysis.

Contribution

It develops the first regression model and minimax bounds for binary responses with covariates in general metric spaces, including Riemannian manifolds.

Findings

Proposed estimator achieves optimal error bounds in key metric spaces.

Classifier's excess risk is bounded and shown to be optimal in Riemannian manifolds.

Numerical experiments confirm the practical utility of the methods.

Abstract

Inspired by logistic regression, we introduce a regression model for data tuples consisting of a binary response and a set of covariates residing in a metric space without vector structures. Based on the proposed model we also develop a binary classifier for metric-space valued data. We propose a maximum likelihood estimator for the metric-space valued regression coefficient in the model, and provide upper bounds on the estimation error under various metric entropy conditions that quantify complexity of the underlying metric space. Matching lower bounds are derived for the important metric spaces commonly seen in statistics, establishing optimality of the proposed estimator in such spaces. Similarly, an upper bound on the excess risk of the developed classifier is provided for general metric spaces. A finer upper bound and a matching lower bound, and thus optimality of the proposed classifier, are established for Riemannian manifolds. To the best of our knowledge, the proposed regression model and the above minimax bounds are the first of their kind for analyzing a binary response with covariates residing in general metric spaces. We also investigate the numerical performance of the proposed estimator and classifier via simulation studies, and illustrate their practical merits via an application to task-related fMRI data.

Tables

Table 1. Table 1: Empirical estimation performance of the proposed method in Case 1 with μ ∗ = I 3 superscript 𝜇 subscript 𝐼 3 \mu^{*}=I_{3} and β ∗ = 3 ​ Σ A ​ R superscript 𝛽 3 subscript Σ 𝐴 𝑅 \beta^{*}=3\Sigma_{AR} . The standard deviations based on the 500 500 500 independent Monte Carlo replicates are given in the parentheses. The value of σ X 2 superscript subscript 𝜎 𝑋 2 \sigma_{X}^{2} is approximated numerically.

${\sigma }_X^2$	$n$	$d({\mu }^{\ast },{\hat{\mu }}_n)$	$d({\beta }^{\ast },{\hat{\beta }}_n)$	RMSE
2.248	100	0.145(0.047)	1.137(0.442)	0.155(0.062)
	500	0.063(0.020)	0.438(0.138)	0.027(0.008)
4.496	100	0.205(0.066)	0.921(0.403)	0.173(0.079)
	500	0.089(0.028)	0.348(0.112)	0.030(0.009)

Table 2. Table 2: Empirical classification performance of classifiers when σ X 2 = 2.248 superscript subscript 𝜎 𝑋 2 2.248 \sigma_{X}^{2}=2.248 . The standard deviations based on the 500 500 500 independent Monte Carlo replicates are in the parentheses.

Model	$n$	Method	Accuracy	AUC	Sensitivity	Specificity	EER
Case 1 (${\mu }^{\ast }=I_3$) (${\beta }^{\ast }=3{\Sigma }_{AR}$)	100	SPDNet	0.609(0.112)	0.657(0.132)	0.603(0.186)	0.615(0.194)	0.040(0.111)
		$k$-NN	0.583(0.118)	0.596(0.100)	0.571(0.227)	0.589(0.234)	0.067(0.133)
		kSVM	0.584(0.121)	0.592(0.101)	0.568(0.287)	0.595(0.283)	0.066(0.135)
		KDC	0.561(0.117)	0.571(0.093)	0.583(0.311)	0.538(0.321)	0.088(0.138)
		Proposed	0.622(0.114)	0.685(0.118)	0.612(0.169)	0.634(0.163)	0.028(0.106)
	500	SPDNet	0.638(0.048)	0.694(0.054)	0.624(0.086)	0.653(0.090)	0.022(0.037)
		$k$-NN	0.621(0.052)	0.621(0.050)	0.624(0.115)	0.619(0.102)	0.039(0.048)
		kSVM	0.642(0.048)	0.642(0.047)	0.647(0.107)	0.637(0.103)	0.018(0.038)
		KDC	0.627(0.054)	0.628(0.049)	0.639(0.192)	0.615(0.186)	0.033(0.052)
		Proposed	0.651(0.049)	0.709(0.053)	0.652(0.073)	0.649(0.071)	0.009(0.030)
Case 2	100	SPDNet	0.636(0.115)	0.670(0.135)	0.631(0.191)	0.637(0.200)	0.147(0.126)
		$k$-NN	0.614(0.115)	0.618(0.107)	0.606(0.228)	0.611(0.218)	0.168(0.132)
		kSVM	0.634(0.116)	0.631(0.107)	0.618(0.266)	0.638(0.254)	0.149(0.133)
		KDC	0.600(0.115)	0.600(0.100)	0.608(0.280)	0.579(0.284)	0.183(0.133)
		Proposed	0.664(0.115)	0.693(0.129)	0.659(0.169)	0.672(0.170)	0.118(0.125)
	500	SPDNet	0.689(0.048)	0.706(0.055)	0.676(0.084)	0.703(0.086)	0.094(0.053)
		$k$-NN	0.680(0.050)	0.681(0.050)	0.685(0.089)	0.676(0.091)	0.102(0.055)
		kSVM	0.717(0.046)	0.717(0.046)	0.724(0.080)	0.710(0.082)	0.066(0.048)
		KDC	0.665(0.065)	0.665(0.064)	0.678(0.145)	0.652(0.151)	0.118(0.067)
		Proposed	0.723(0.047)	0.720(0.054)	0.726(0.069)	0.720(0.070)	0.060(0.048)
Case 3	100	SPDNet	0.707(0.105)	0.724(0.131)	0.422(0.216)	0.852(0.110)	0.060(0.099)
		$k$-NN	0.673(0.111)	0.585(0.105)	0.867(0.136)	0.300(0.240)	0.093(0.118)
		kSVM	0.673(0.106)	0.550(0.091)	0.933(0.102)	0.166(0.228)	0.094(0.114)
		KDC	0.655(0.112)	0.534(0.077)	0.929(0.147)	0.129(0.232)	0.112(0.122)
		Proposed	0.638(0.109)	0.727(0.123)	0.599(0.141)	0.717(0.189)	0.128(0.120)
	500	SPDNet	0.726(0.047)	0.762(0.054)	0.454(0.100)	0.867(0.051)	0.048(0.044)
		$k$-NN	0.704(0.047)	0.618(0.051)	0.890(0.062)	0.346(0.121)	0.069(0.048)
		kSVM	0.717(0.048)	0.633(0.056)	0.901(0.057)	0.365(0.128)	0.056(0.046)
		KDC	0.660(0.048)	0.527(0.047)	0.948(0.105)	0.105(0.182)	0.113(0.054)
		Proposed	0.667(0.051)	0.755(0.054)	0.627(0.064)	0.746(0.080)	0.106(0.052)

Table 3. Table 3: Comparison of classifiers.

Method	Accuracy	AUC	Sensitivity	Specificity
SPDNet	0.8667	0.9205	0.8613	0.8690
$k$-NN	0.8333	0.8349	0.8401	0.8297
kSVM	0.8436	0.8444	0.8667	0.8221
KDC	0.7744	0.7811	0.8313	0.7309
Logistic-BOLD	0.7641	0.8120	0.8490	0.6984
Proposed	0.8718	0.9060	0.9158	0.8376

Equations

$$\text{log}\bigg{(}\frac{P_{X}}{1-P_{X}}\bigg{)}=d(\mu^{*},X)d(\mu^{*},\beta^{*})\cos(\angle_{\mu^{*}}(X,\beta^{*})),$$

$${\bar{\angle }}_p(q,r)=\arccos \frac{d^2(p,q)+d^2(p,r)-d^2(q,r)}{2d(p,q)d(p,r)}.$$

$${\angle }_p(\gamma ,\eta )=\underset{r\to 0}{\lim }\underset{0<t,s<r}{\sup }{\bar{\angle }}_p(\gamma (t),\eta (s)).$$

$$({\text{Hess}}_{p}f)(\xi )={\nabla }_{\xi }\text{grad }f$$

$$\text{logit}(P_X)=⟨X-{\mu }^{\ast },{\beta }^{\ast }-{\mu }^{\ast }{⟩}_E=d_E({\mu }^{\ast },X)d_E({\mu }^{\ast },{\beta }^{\ast })\cos ({\angle }_{{\mu }^{\ast }}(X,{\beta }^{\ast })),$$

$${\mu }^{\ast }=\underset{\mu \in \mathcal{M}}{\mathrm{argmin}}\mathbb{E}[d^2(X,\mu )],$$

$$\text{logit}(P_X)=h({\beta }^{\ast };X,{\mu }^{\ast }):=d({\mu }^{\ast },X)d({\mu }^{\ast },{\beta }^{\ast })\cos ({\angle }_{{\mu }^{\ast }}(X,{\beta }^{\ast })),$$

$${\hat{\mu }}_n=\underset{\mu \in \mathcal{M}}{\mathrm{argmin}}\frac{1}{n}\sum _{i=1}^n{d}^2(X_i,\mu ).$$

$${\beta }^{\ast }=\underset{\beta \in \mathcal{M}}{\mathrm{argmax}}L(\beta )\text{with}L(\beta ):=\mathbb{E}\left[Yh(\beta ;X,{\mu }^{\ast })-\log (1+e^{h(\beta ;X,{\mu }^{\ast })})\right].$$

$$L_n(\beta )$$

$${\hat{\beta }}_n\in \underset{\beta \in \mathcal{M}}{\mathrm{argmax}}{L}_n(\beta ).$$

$$|d(p,q)\cos ({\angle }_p(u,q))-d(p,r)\cos ({\angle }_p(u,r))|\le C_U{d}^{{\alpha }_U}(q,r).$$

$$\mathbb{E}\big{[}|h(\beta;X,\mu^{*})-h(\beta^{*};X,\mu^{*})|^{2}\big{]}\geq C_{L}\lambda_{X}d^{\alpha_{L}}(\beta,\beta^{*}).$$

$$v_{L,n}={\left(\frac{D}{n}\right)}^{\frac{1}{2({\alpha }_L-{\alpha }_U)}}$$

$$v_{L,n}=\begin{cases}\big{(}\frac{1}{n}\big{)}^{\frac{1}{2(\alpha_{L}-\alpha_{U})}}&if~{}D<2\alpha_{U},\\ \big{(}\frac{\log^{2}n}{n}\big{)}^{\frac{1}{2(\alpha_{L}-\alpha_{U})}}&if~{}D=2\alpha_{U},\\ {\big{(}\frac{1}{n}\big{)}^{\frac{\alpha_{U}}{D(\alpha_{L}-\alpha_{U})}}}&if~{}D>2\alpha_{U}\end{cases}$$

$$\mathbb{E}\big{[}d^{2}(\beta^{*},\hat{\beta}_{n})\big{]}\leq cv_{L,n}^{2}.$$

$$\underset{{\hat{\beta }}_n}{\inf }\underset{\beta \in \mathcal{M}}{\sup }{\mathbb{E}}_{\beta }[d^2(\beta ,{\hat{\beta }}_n)]\gtrsim \{\begin{array}{ll}{\left(\frac{D}{{\sigma }_X^{2}n}\right)}^{\frac{1}{{\alpha }_U}}& \text{ when }N\left(ϵ,B({\beta }^{\ast },\delta ),d\right)\asymp {\left(\frac{\delta }{ϵ}\right)}^D,\\ {\left(\frac{1}{{\sigma }_X^{2}n}\right)}^{\frac{1}{{\alpha }_U}}& \text{ when }\log N\left(ϵ,B({\beta }^{\ast },\delta ),d\right)\asymp {\left(\frac{\delta }{ϵ}\right)}^D.\end{array}$$

$$\mathbb{E}[d^2(\beta ,{\hat{\beta }}_n)]\asymp D/n.$$

$$v_{G,n}={\left(\frac{D}{n}\right)}^{\frac{1}{2({\alpha }_L-{\alpha }_U+{\gamma }_0)}}$$

$$v_{G,n}=\begin{cases}\big{(}\frac{1}{n}\big{)}^{\frac{1}{2(\alpha_{L}-\alpha_{U}+D/2)}}&if~{}D<2\alpha_{U},\\ \big{(}\frac{\log^{2}n}{n}\big{)}^{\frac{1}{2\alpha_{L}}}&if~{}D=2\alpha_{U},\\ {\big{(}\frac{1}{n}\big{)}^{\frac{\alpha_{U}}{D\alpha_{L}}}}&if~{}D>2\alpha_{U}\end{cases}$$

$$\mathbb{E}\big{[}d^{2}(\beta^{*},\hat{\beta}_{n})\big{]}\leq cv_{G,n}^{2}.$$

$$\underset{{\hat{\beta }}_n}{\inf }\underset{\beta \in \mathcal{M}}{\sup }{\mathbb{E}}_{\beta }[d^2(\beta ,{\hat{\beta }}_n)]\gtrsim \{\begin{array}{ll}{\left(\frac{D}{{\sigma }_X^{2}n}\right)}^{\frac{1}{{\alpha }_U}}& \text{ when }N\left(ϵ,\mathcal{M},d\right)\asymp {\left(\frac{1}{ϵ}\right)}^D,\\ {\left(\frac{1}{{\sigma }_X^{2}n}\right)}^{\frac{1}{{\alpha }_U+D/2}}& \text{ when }\log N\left(ϵ,\mathcal{M},d\right)\asymp {\left(\frac{1}{ϵ}\right)}^D.\end{array}$$

$$\hat{s}(X)=\{\begin{array}{rl}1& \text{if}{\hat{P}}_X\ge 1/2\text{or, equivalently, }h({\hat{\beta }}_n;X,{\hat{\mu }}_n)\ge 0,\\ 0& \text{otherwise}.\end{array}$$

$$s^{\ast }(X)=\{\begin{array}{rl}1& \text{if}{P}_X\ge 1/2\text{or, equivalently, }h({\beta }^{\ast };X,{\mu }^{\ast })\ge 0,\\ 0& \text{otherwise}.\end{array}$$

$$\mathcal{E}(s^{*},\hat{s})\leq c\left(\mathbb{E}\left[d^{2}(\mu^{*},\hat{\mu}_{n})\right]+\mathbb{E}\big{[}d^{2}(\beta^{*},\hat{\beta}_{n})\big{]}\right)^{1/2}.$$

$$\mathcal{E}(s^{\ast },\hat{s})\le c\sqrt{D/n}.$$

$$\underset{\hat{s}}{\inf }\underset{s^{\ast }\in \mathcal{C}}{\sup }\mathcal{E}(s^{\ast },\hat{s})\gtrsim \sqrt{D/n},$$

$$\mathrm{RMSE}=\sqrt{\frac{1}{|{\mathcal{D}}_{\mathrm{test}}|}\sum _{i\in {\mathcal{D}}_{\mathrm{test}}}{\left(\text{logit}({\hat{P}}_{X_i})-\text{logit}(P_{X_i})\right)}^2}.$$

$$\mathrm{EER}=\frac{1}{|{\mathcal{D}}_{\mathrm{test}}|}\sum _{i\in {\mathcal{D}}_{\mathrm{test}}}\{\mathcal{L}(\hat{s}(X_i))-\mathcal{L}(s^{\ast }(X_i))\},$$

Full text

Keywords: Alexandrov geometry, excess risk, Hadarmard space, logistic regression, manifold, minimax optimality.

MSC2020 subject classification: 62R20, 62J12, 62H30.

1 Introduction

Object data sampled from general metric spaces emerge recently in many domains, posing fundamental challenges to traditional statistical methodologies geared to vector structures. Examples of such data include phylogenetic trees arising from evolutionary biology (Billera et al., 2001), distribution functions sampled from the nonlinear Wasserstein geometry for modeling differentiating cells in developmental biology (Schiebinger et al., 2019), symmetric positive-definite matrices from a nonlinear Riemannian manifold for representing diffusion tensors and brain connectivities in neuroscience (e.g., Arsigny et al., 2007; Dryden et al., 2009; Yuan et al., 2012), graph Laplacians for studying networks (Dubey and Müller, 2022), shapes of tumors (Bharath and Kurtek, 2020) in medical sciences, and more in Marron and Dryden (2021). Although there is a body of literature on regression analysis for this type of data (e.g., Hein, 2009; Faraway, 2014; Petersen and Müller, 2019; Qiu et al., 2022; Ying and Yu, 2022; Tucker et al., 2023; Zhang et al., 2023, among many others), very few studies explore parametric regression models specifically designed for binary response variables with metric-space valued covariates.

In the analysis of a binary response variable linked to a set of covariates, the crucial concepts of odds and log odds emerge across various fields, with logistic regression serving as a widely used parametric model to examine the impact of covariates on the odds. Motivated by this background, this paper presents a new model designed to analyze a binary response affected by covariates that reside in a metric space $\mathcal{M}$ without a vector structure. Specifically, for a binary variable $Y$ and a covariate $X\in\mathcal{M}$, we propose the model

[TABLE]

where $P_{X}:=\mathbb{P}(Y=1|X)$ is success probability, $d$ is the distance on $\mathcal{M}$, $\mu^{*}\in\mathcal{M}$ is the (Fréchet) mean of $X$, $\beta^{*}$ is the coefficient object to be estimated, and $\angle_{\mu^{*}}(X,\beta^{*})$ is the angle between the geodesics emanating from $\mu^{*}$ to $X$ and $\beta^{*}$. The left-hand side of (1) represents the log odds. The right-hand side of (1), referred to as the Alexandrov inner product between the two geodesics, reduces to $(X-\mu^{*})^{\top}(\beta^{*}-\mu^{*})$ when $\mathcal{M}=\mathbb{R}^{D}$, thus mimicking and extending the linear covariate effect on the log odds. Similarly, $d(\mu^{*},X)\cos(\angle_{\mu^{*}}(X,\beta^{*}))$ mimics the projection of the geodesic connecting $\mu^{*}$ to $X$ onto the geodesic connecting $\mu^{*}$ to $\beta^{*}$; see Figure 1 for an illustration. Intuitively, the parameter $\beta^{*}$ points to a “direction” in which the log odds changes at the fastest rate, and the influence of $X$ on the log odds is determined by its projection onto this direction. This direction could therefore offer insights into the association between the response and the covariates, highlighting the practical merit of the proposed model (1); see Section 7 for a demonstration. In Section 3, we provide further details about (1) and develop an efficient estimator $\hat{\beta}_{n}$ for $\beta^{*}$ based on the principle of maximum likelihood. To the best of our knowledge, this represents the first parametric model designed for the odds and its association with covariates from general metric spaces.

The proposed model (1) also induces a classifier for data residing in general metric spaces. Via modeling the log odds in a parametric form, this classifier is less prone to the curse of dimensionality, and thus is preferable when interpretable results are desired or when the sample size is small relative to the complexity of $\mathcal{M}$; see Section 7 for a numeric demonstration. Meanwhile we should remark that, when the association between log odds and covariates is highly nonlinear, dedicated classifiers (e.g., Gottlieb et al., 2016; Qiao et al., 2019; Huang and Van Gool, 2017, among others), which are generally nonparametric, may produce more accurate classification results at the cost of interpretability; see Section 6 for a numeric illustration.

In addition to the aforementioned methodological contributions, we establish upper and matching lower bounds on the estimation error $d^{2}(\hat{\beta}_{n},\beta^{*})$ under various metric entropy conditions that characterize complexity of a metric space. This demonstrates the optimality of the proposed estimator for metric spaces commonly encountered in statistics. For instance, under certain regularity conditions, the proposed estimator $\hat{\beta}_{n}$ achieves the optimal rate $\mathbb{E}[d^{2}(\hat{\beta}_{n},\beta^{*})]\asymp n^{-1/\alpha}$, where $\alpha$ is linked to the geometric characteristics of $\mathcal{M}$; see Section 4 for details. To the best of our knowledge, this is the first optimality result for regression of logistic type like (1) on general metric spaces.

Another major theoretical contribution is the establishment of upper and corresponding lower bounds on the excess risk of the classifier induced by the proposed mode (1). Roughly speaking, we show that, when the covariate resides in a $D$-dimensional Riemannian manifold, the excess risk decays at the optimal rate $\sqrt{D/n}$; see Section 5 for details. To our knowledge, this represents the first optimality result in classification for manifold-valued objects. Notably, our results encompass the optimal rate in Abramovich and Grinshtein (2019) for Euclidean logistic regression as a special case. Additionally, we develop an upper bound on the excess risk for general metric spaces. This upper bound, presented in Theorem 5 and expressed in terms of the estimation error, characterizes the influence of the geometric features and complexity of the underlying metric space on the excess risk.

Theoretical analysis of the proposed estimator and classifier is challenged by the absence of a vector structure in a general metric space. For instance, the regression model (1) is defined geometrically, and our theoretical analysis is significantly complicated by the nonlinear metric geometry of the underlying space, including but not limited to the space curvature. This contrasts with the Euclidean case where the presence of rich vector and algebraic structures makes the analysis much more tractable. In addition, the likelihood functions (see Section 3), combining geometric constructions with the nonlinear logistic function, introduces extra complexity into our theoretical studies; see Claim LABEL:lemma:maximal-ineq-expectation, Corollary LABEL:cor:fmse and Lemma LABEL:lem:mid-seg-comparison-model-space of Supplementary Material for examples.

The rest of the paper is organized as follows. Section 2 provides a brief introduction on metric and Riemannian geometry. In Section 3, we introduce the proposed model (1) in details and develop an estimator for the model parameter. We study theoretical properties of the estimator in Section 4 and the excess risk of classification in Section 5. Section 6 is devoted to assessing the empirical performance of the proposed method via simulation studies. An application to task-related functional magnetic resonance imaging (fMRI) data is presented in Section 7. We conclude the article in Section 8. All proofs are relegated to Supplementary Material.

2 Background on Metric and Riemannian Geometry

Below we briefly review some geometric concepts that are essential to our developments, while refer readers to Burago et al. (2001) and Lee (2018) for a comprehensive treatment on metric and Riemannian geometry; see also Holopainen (2009) and Section S1 of Shao et al. (2022) for self-contained notes on basic concepts of metric and Riemannian geometry.

Geodesic. For a metric space $(\mathcal{M},d)$ and two points $p,q\in\mathcal{M}$, a geodesic connecting $p$ and $q$ is a map $\gamma$ from a closed interval $[0,T]\subset\mathbb{R}$ to $\mathcal{M}$ with $T=d(p,q)$ such that, $\gamma(0)=p$, $\gamma(T)=q$ and $d(\gamma(t),\gamma(t^{\prime}))=|t-t^{\prime}|$ for all $t,t^{\prime}\in[0,T]$. The image of $\gamma$ on $[0,T]$ is called a geodesic segment. Intuitively, a geodesic segment is a path/curve with a shortest length (with respect to the underlying metric) that connects two given points. A metric space $(\mathcal{M},d)$ is a geodesic space if every two points in $\mathcal{M}$ can be connected by a geodesic, and is a uniquely geodesic space if this geodesic is unique. We use $\gamma_{p}^{q}$ to denote the geodesic connecting $p$ and $q$ in a uniquely geodesic space.

Curvature. Curvature quantifies the deviation of a general metric space from being flat. A standard way to define curvature is to compare geodesic triangles on the underlying space with geodesic triangles on the following reference spaces $(M_{\kappa}^{2},d_{\kappa})$:

• •

When $\kappa=0$, $M_{\kappa}^{2}=\mathbb{R}^{2}$ with the standard Euclidean distance;

• •

When $\kappa<0$, $M_{\kappa}^{2}$ is the hyperbolic space $\mathbb{H}^{2}=\{(x,y,z)\in\mathbb{R}^{3}:x^{2}+y^{2}-z^{2}=-1\text{ and }z>0\}$ with the hyperbolic distance function $d_{\kappa}(p,q)=\cosh^{-1}(z_{p}z_{q}-x_{p}x_{q}-y_{p}y_{q})/\sqrt{-\kappa}$, where $p=(x_{p},y_{p},z_{p})$ and $q=(x_{q},y_{q},z_{p})$;

• •

When $\kappa>0$, $M_{\kappa}^{2}$ is the sphere $\mathcal{S}^{2}=\{(x,y,z)\in\mathbb{R}^{3}:x^{2}+y^{2}+z^{2}=1\}$ with the angular distance function $d_{\kappa}(p,q)=\cos^{-1}(x_{p}x_{q}+y_{p}y_{q}+z_{p}z_{q})/\sqrt{\kappa}$.

A geodesic triangle $\triangle(p,q,r)$ with vertices $p,q,r$ in a uniquely geodesic space $\mathcal{M}$ is formed by three geodesic segments, termed sides, that connect $p$ to $q$, $p$ to $r$ and $q$ to $r$, respectively. The sum $d(p,q)+d(p,r)+d(r,q)$ is called the perimeter of $\triangle(p,q,r)$. A comparison triangle of $\triangle(p,q,r)$ in the reference space $M_{\kappa}^{2}$ is a geodesic triangle on $M_{\kappa}^{2}$ consisting of vertices $\bar{p},\bar{q},\bar{r}\in M_{\kappa}^{2}$ and sides $\gamma_{\bar{p}}^{\bar{q}},\gamma_{\bar{p}}^{\bar{r}},\gamma_{\bar{q}}^{\bar{r}}\subset M_{\kappa}^{2}$ such that $d_{\kappa}(\bar{p},\bar{q})=d(p,q)$, $d_{\kappa}(\bar{p},\bar{r})=d(p,r)$ and $d_{\kappa}(\bar{q},\bar{r})=d(q,r)$. In addition, $\bar{x}\in\gamma_{\bar{p}}^{\bar{q}}$ is a comparison point of $x\in\gamma_{p}^{q}$ if $d_{\kappa}(\bar{p},\bar{x})=d(p,x)$. Let $D_{\kappa}=\pi/\sqrt{\kappa}$ if $\kappa>0$ and $D_{\kappa}=\infty$ otherwise. We say the (global) curvature of $\mathcal{M}$ is lower (upper, respectively) bounded by $\kappa$ if every geodesic triangle with perimeter less than $2D_{\kappa}$ satisfies the following property: There exists a comparison triangle $\triangle(\bar{p},\bar{q},\bar{r})$ in $M_{\kappa}^{2}$ such that $d(x,y)\geq d_{\kappa}(\bar{x},\bar{y})$ ($d(x,y)\leq d_{\kappa}(\bar{x},\bar{y})$, respectively) for all $x\in\gamma_{p}^{q}$ and $y\in\gamma_{p}^{r}$ and their comparison points $\bar{x}$ and $\bar{y}$ on $\triangle(\bar{p},\bar{q},\bar{r})$. Metric spaces with lower or upper bounded curvature are called Alexandrov spaces, and complete geodesic spaces with curvature upper bounded by [math] are called Hadamard spaces.

Angle. For a metric space $(\mathcal{M},d)$ and three distinct points $p,q,r\in\mathcal{M}$, the comparison angle between $q$ and $r$ at $p$, denoted by $\bar{\angle}_{p}(q,r)$, is defined by

[TABLE]

Geometrically, for $p,q,r\in\mathcal{M}$ that form a geodesic triangle $\triangle(p,q,r)$, the comparison angle $\bar{\angle}_{p}(q,r)$ is the Euclidean angle at $\bar{p}$ of the comparison triangle $\triangle(\bar{p},\bar{q},\bar{r})$ in $\mathbb{R}^{2}$. The Alexandrov angle between two geodesics $\gamma$ and $\eta$ emanating from $p$ in a uniquely geodesic space, denoted by $\angle_{p}(\gamma,\eta)$, is defined by

[TABLE]

In a curved space, the comparison angle $\bar{\angle}_{p}(\gamma(t),\eta(s))$ varies with $t$ and $s$ even though the directions of $\gamma$ and $\eta$ remain unchanged. In contrast, the Alexandrov angle $\angle_{p}(\gamma,\eta)$ intrinsically depends only on the directions of $\gamma$ and $\eta$ (i.e., invariant to $t$ and $s$), but not on their lengths. For three distinct points $p,q,r\in\mathcal{M}$, we define the angle $\angle_{p}(q,r)=\angle_{p}(\gamma_{p}^{q},\gamma_{p}^{r})$.

Riemannian Manifolds. A differentiable manifold is a topological manifold with a maximal differentiable atlas that induces a differential structure on the manifold. When all transition maps in the atlas are $C^{\infty}$ differentiable, the manifold is a smooth manifold. For each point $p$ on a smooth manifold $\mathcal{M}$, there is a linear space $T_{p}\mathcal{M}$ attached to $p$, called the tangent space at $p$. The elements in $T_{p}\mathcal{M}$ are called tangent vectors. A Riemannian manifold is a smooth manifold $\mathcal{M}$ with a Riemannian metric that defines a positive-definite inner product $\langle\cdot,\cdot\rangle_{p}$ on the tangent space $T_{p}\mathcal{M}$ for each $p\in\mathcal{M}$ and smoothly varying with $p$. An affine connection on a smooth manifold provides a way to compare tangent vectors at distinct points on the manifold via the notion of parallel transport. Each Riemannian manifold has a canonical affine connection, the Levi–Civita connection, which is torsion-free and compatible with the Riemannian metric. For a smooth function $f:\mathcal{M}\rightarrow\mathbb{R}$, the Hessian of $f$ at $p\in\mathcal{M}$ is defined by

[TABLE]

for all $\xi\in T_{p}\mathcal{M}$, where $\nabla$ is the Levi–Civita connection and grad is the gradient operator (Page 45, Lee, 2018). It can be seen that is a linear operator on the tangent space $T_{p}\mathcal{M}$.

3 Binary Regression in Metric Spaces

3.1 Model

To motivate our model for a binary response with a covariate residing in a metric space, we observe that the Euclidean logistic regression model with the covariate $X=(X_{1},\ldots,X_{D})\in\mathbb{R}^{D}$, when appropriately reparameterized, can be represented as follows:

[TABLE]

where $P_{X}=\mathbb{P}(Y=1\mid X)$ is the success probability, $\text{logit}(t)=\log(\frac{t}{1-t})$ is the logit function for $t\in[0,1]$, ${\beta}^{*}=({\beta}_{1}^{*},\ldots,{\beta}_{D}^{*})$ is the true coefficient vector, $\langle\cdot,\cdot\rangle_{E}$ is the standard Euclidean inner product, $\mu^{*}=\mathbb{E}X$ is the mean, $d_{E}$ is the standard distance on $\mathbb{R}^{D}$, and $\angle_{\mu^{*}}(X,\beta^{*})$ is the Euclidean angle between $X$ and $\beta^{*}$ at $\mu^{*}$. Without loss of generality, in the above we have assumed that the intercept is zero; an extension to the general case is straightforward (see also Remark 1).

The logistic regression model in the form of (4) inspires us to model the log odds influenced by a metric-space valued covariate $X$, as follows. Let $(\mathcal{M},d)$ be a uniquely geodesic metric space and $X\in\mathcal{M}$, where $d$ is the distance on $\mathcal{M}$. First, the Euclidean mean is replaced with its popular extension in general metric spaces, i.e., the Fréchet mean

[TABLE]

which is well studied in the literature and is reduced to the Euclidean mean $\mathbb{E}X$ when $\mathcal{M}=\mathbb{R}^{D}$. In addition, the angle $\angle_{\mu^{*}}(X,\beta^{*})$ in (4) is naturally generalized to general metric spaces by the Alexandrov angle between the geodesics $\gamma_{\mu^{*}}^{X}$ and $\gamma_{\mu^{*}}^{\beta^{*}}$ emanating from $\mu^{*}$. Then we propose to model the effect of the metric-space valued covariate $X\in\mathcal{M}$ on the log odds by

[TABLE]

where $\beta^{*}\in\mathcal{M}$ is a population parameter to be estimated and $\angle_{\mu^{*}}(X,\beta^{*})=\angle_{\mu^{*}}(\gamma_{\mu^{*}}^{X},\gamma_{\mu^{*}}^{\beta^{*}})$ is the Alexandrov angle.

The right-hand side of (6), termed the Alexandrov inner product (Lin and Müller, 2021), is interpreted as the conditional log odds of $Y=1$ and $Y=0$ given $X$. In addition, $d(\mu^{*},X)\cos(\angle_{\mu^{*}}(X,\beta^{*}))$ in (6) is interpreted as the projection of the geodesic $\gamma_{\mu^{*}}^{X}$ onto the geodesic $\gamma_{\mu^{*}}^{\beta^{*}}$, and thus $d(\mu^{*},\beta^{*})$ represents the expected change in log odds based on a one-unit change of this projection. Alternatively, the projection $d(\mu^{*},\beta^{*})\cos(\angle_{\mu^{*}}(X,\beta^{*}))$ represents the expected change in log odds based on a one-unit change of $X$ along the geodesic $\gamma_{\mu^{*}}^{X}$ emanating from $\mu^{*}$ to $X$; in Section LABEL:sec:odds-ratio of Supplementary Material, we also show that $d(\mu^{*},\beta^{*})\cos(\angle_{\mu^{*}}(X,\beta^{*}))$ represents a log odds ratio. These indicate that the geodesic $\gamma_{\mu^{*}}^{\beta^{*}}$ identifies the direction in which the log odds change most rapidly, thereby encapsulating critical information of the covariates for distinguishing the two classes of the response variable; see the applications below for a demonstration.

The proposed model (6) has at least two potential applications. First, based on the relationship between the log odds and the geodesic $\gamma_{\mu^{*}}^{\beta^{*}}$ in the above, practitioners may use (6) as a regression model to gain insights into distinctive patterns in covariates between the two classes via examining an estimate of the geodesic $\gamma_{\mu^{*}}^{\beta^{*}}$. For an illustration, take the example of investigating how the brain functional connectivity is associated with different tasks related to certain human functions (Salvo et al., 2021), such as the motor and language processing tasks. To study this question of scientific interest, one may focus on 8 brain regions associated with these two tasks and utilize a set of $8\times 8$ SPD (symmetric positive-definite) matrices derived from fMRI images, each matrix representing connectivity among these regions when a subject performs a task. Applying model (6) and a maximum likelihood estimator (given below) to analyze this data provides the geodesic estimate $\gamma_{\hat{\mu}}^{\hat{\beta}}(t)$, as shown in Figure 2. As $t$ increases, it is observed that Region 1 (the superior temporal gyrus) becomes more active, while its connectivity with Regions 5–8 (encompassing the primary motor and visual cortex) diminishes. Thus, practitioners may infer that there is a reduced functional interaction between the superior temporal gyrus and the primary motor and visual cortex during the language processing task. Further demonstration of this application is available in Section 7.

Secondly, the classifier derived from (6), which classifies $X$ into the class $Y=1$ if $P_{X}\geq 0.5$ and into $Y=0$ otherwise, allows practitioners to classify a metric-space valued covariate $X$ based on an estimate of the conditional success probability $P_{X}=\mathbb{P}(Y=1\mid X=x)=\{1+\exp(-h(\beta^{*};x,\mu^{*}))\}^{-1}$; see Section 5 for more details. For instance, in the context of analyzing functional connectivity and task-specific brain functions, the classifier can be used to predict the task being executed, utilizing the observed $8\times 8$ SPD matrix that depicts connectivity among the 8 brain regions. As outlined in Section 7, our findings suggest that this model-based classifier can surpass nonparametric methods like the $k$-nearest-neighbors algorithm and neural networks in performance. This could be attributed to the model’s interpretability and its efficiency in handling the relatively high dimensionality inherent to the data of $8\times 8$ SPD matrices, which have an intrinsic dimension of 36, where nonparametric classifiers might struggle due to the curse of dimensionality.

Remark 1*.*

The model (6) can naturally accommodate both Euclidean and metric-space valued covariates. For instance, suppose $X_{1}\in\mathbb{R}^{D}$ and $X_{2}\in\mathcal{M}$. To include both $X_{1}$ and $X_{2}$ into the proposed model, we can treat $X=(X_{1},X_{2})$ jointly as a covariate on the product metric space $\mathbb{R}^{D}\times\mathcal{M}$, and then apply the model (6) to $X$. For instance, with $X_{1}\equiv 1$, we can include an intercept term in the model (6). A similar strategy can be used to accommodate multiple covariates potentially residing in different metric spaces.

Remark 2*.*

The proposed model (6) encompasses the Euclidean logistic model as a special case when $\mathcal{M}=\mathbb{R}^{D}$. In light of the resemblance between our model (6) and the Euclidean logistic model (4), the former may be referred to as the metric-space logistic regression model, although it does not possess all properties of the latter due to the lack of Euclidean geometry and algebraic structures in (6). Nonetheless, we opt to use the term “binary regression” in this paper rather than “logistic regression” to prevent potential misconceptions.

The strength of the proposed model (6) stems from its resistance to the curse of dimensionality (when compared to nonparametric methods) and its ability to help practitioners investigate the relationship between the log odds and metric-space valued covariates, as demonstrated in the above and in Section 7. Alternative models that are inspired by other features of Euclidean logistic regression are beyond the scope of this paper and left for future studies.

3.2 Maximum Likelihood Estimator

To estimate the parameters $\mu^{*}$ and $\beta^{*}$, let $(X_{i},Y_{i})$, $i=1,\ldots,n$, be independent and identically distributed (i.i.d.) observations of $(X,Y)\sim F$ for a joint distribution $F$ on $\mathcal{M}\times\{0,1\}$. We first estimate the Fréchet mean $\mu^{*}$ by the sample Fréchet mean

[TABLE]

For $\beta^{*}$, according to the proposition below, we observe that it is the maximizer of the population log-likelihood function $L(\beta)$, i.e.,

[TABLE]

Proposition 1**.**

If $(X,Y)$ is a random pair sampled from the model (6) with the parameters $\beta^{*}$ and $\mu^{*}$, then $\beta^{*}\in\operatorname*{argmax}_{\beta\in\mathcal{M}}L(\beta)$. Moreover, $\beta^{*}=\operatorname*{argmax}_{\beta\in\mathcal{M}}L(\beta)$ when $L(\beta)$ has a unique maximizer.

This motivates us to estimate $\beta^{*}$ by maximizing the empirical version of $L(\beta)$,

[TABLE]

where $h(\beta;X_{i},\hat{\mu}_{n})=d(\hat{\mu}_{n},X_{i})d(\hat{\mu}_{n},\beta)\cos(\angle_{\hat{\mu}_{n}}(X_{i},\beta))$. The maximum likelihood estimator (MLE) of $\beta^{*}$ is then

[TABLE]

Existence and uniqueness of $\hat{\beta}_{n}$, as well as $\mu^{*}$, $\beta^{*}$ and $\hat{\mu}_{n}$, are discussed in Section 4.1.

Optimization methods for (9) in both general metric spaces and Riemannian manifolds, including how to compute the Alexandrov inner product in these spaces, are detailed in Section LABEL:sec:RGD of Supplementary Material.

4 Theoretical Properties for Estimation

4.1 Assumptions

To study the theoretical properties of $\hat{\beta}_{n}$ in (9), we require the following two geometric assumptions on the metric space $\mathcal{M}$.

Assumption 1*.*

$\mathcal{M}$ is a complete uniquely geodesic metric space.

Assumption 2*.*

There exist constants $C_{U}>0$ and $\alpha_{U}\in(0,1]$ such that, for all $p,q,r,u\in\mathcal{M}$,

[TABLE]

Assumption 1 is standard in metric-space data analysis. With this assumption, any pair of points in $\mathcal{M}$ can be connected by a unique geodesic. Assumption 2 is a regularity condition on the curvature of $\mathcal{M}$, which plays a similar role of the weak quadruple inequality in Schötz (2019). These two assumptions are satisfied with $\alpha_{U}=1$ for Hadamard spaces and some subspaces of positively curved Alexandrov spaces, according to Propositions LABEL:prop:had-lip and LABEL:prop:alex-lip in Section LABEL:sec:proposition-lip-metric of Supplementary Material.

Example*.*

Among many others, below are some concrete examples of metric spaces satisfying both Assumption 1 and Assumption 2 (with $\alpha_{U}=1$).

1. (a)

The space of SPD matrices, endowed with either of the Riemannian metrics of Moakher (2005); Arsigny et al. (2007); Lin (2019), for modeling diffusion tensors and brain functional connectivity in brain sciences. 2. (b)

The positive quadrant $\{(x_{0},\ldots,x_{D})\in\mathcal{S}^{D}:x_{j}\geq 0\text{ for }j=0,\ldots,D\}$ of the unit sphere $\mathcal{S}^{D}$, endowed with the angular distance, for modeling compositional data of $D+1$ components. 3. (c)

The tree space (Billera et al., 2001) in evolutionary biology. 4. (d)

The Wasserstein space of probability distributions defined on a compact domain of $\mathbb{R}$, endowed with the optimal transport geometry, for analyzing human mortality data (Chen and Müller, 2021).

We also require the following two distributional assumptions on $(X,Y)$. Their variants are commonly adopted for metric-space data analysis, and their Euclidean counterparts are commonly seen in the literature of $M$-estimation; see the discussions below.

Assumption 3*.*

The distribution of $X\in\mathcal{M}$ has a bounded support, and the objects $\mu^{*}$ and $\hat{\mu}_{n}$ exist and are unique, the later almost surely. Moreover, $\hat{\beta}_{n}$ exists almost surely.

Assumption 4*.*

There exist constants $\lambda_{X}>0$, $C_{L}>0$ and $\alpha_{L}>1+{2/m}$ for some constant $m\geq 4$ with $\mathbb{E}[d^{m}(\mu^{*},X)]<\infty$, such that, for any $\beta\in\mathcal{M}$,

[TABLE]

Existence and uniqueness of the Fréchet means $\mu^{*}$ and $\hat{\mu}_{n}$, often presumed in data analysis within general metric spaces, are influenced by the distribution of $X$. Nonetheless, in certain metric spaces, the requirement of a finite second moment of $X$ is adequate. For example, for any Hadamard space $\mathcal{M}$, both $\mu^{*}$ and $\hat{\mu}_{n}$ exist and are unique when $\mathbb{E}[d^{2}(X,u)]<\infty$ for some $u\in\mathcal{M}$ (Proposition 4.3, Sturm, 2003), regardless of the distribution of $X$. For other metric spaces, a convexity condition on the distance function can guarantee existence and uniqueness of both $\mu^{*}$ and $\hat{\mu}_{n}$ when data reside in a geodesically convex region (Lin and Müller, 2021). See also Bhattacharya and Patrangenaru (2003); Afsari (2011) for sufficient conditions on $X$ when $\mathcal{M}$ is a Riemannian manifold.

For $\beta^{*}$, Proposition 1 establishes its existence. In Section LABEL:subsec:uniqueness-MLE of Supplementary Material, Proposition LABEL:prop:unique-MLE further establishes the uniqueness of $\beta^{*}$ under additional Assumption 4, with two concrete examples provided in the same section. For $\hat{\beta}_{n}$, our theoretical analysis does not require uniqueness of $\hat{\beta}_{n}$, since the theoretical results in the sequel hold for any measurable choice of the maximizers $\hat{\beta}_{n}$ in (9); such measurability is commonly adopted in the M-estimation theory (e.g., Pages 285–286 of van der Vaart and Wellner, 1996) and object data analysis (e.g., Theorem 2.3 of Bhattacharya and Patrangenaru, 2003).

The existence of the MLE $\hat{\beta}_{n}$ appears to be relatively subtle. When $\mathcal{M}$ is complete and totally bounded, $\hat{\beta}_{n}$ exists (Proposition LABEL:prop:exist-MLE in Supplementary Materials). For unbounded metric spaces, specific distributional conditions for the pair $(X,Y)$ are required. This necessity applies even within Euclidean spaces $\mathcal{M}=\mathbb{R}^{D}$, where conditions such as the overlap condition (Albert and Anderson, 1984) and the proper convexity condition (Silvapulle, 1981) are pertinent. To illustrate this type of distributional conditions outside of Euclidean realms, in Section LABEL:subsec:ex-as-exist of Supplementary Material we construct an analogous condition through the example of the SPD matrix manifold. A detailed investigation into the distributional conditions on $(X,Y)$ in the context of metric spaces, whether sufficient or necessary for $\hat{\beta}_{n}$’s existence, is a subject deserving of its own comprehensive study in future research.

Assumption 4 is common in the $M$-estimation theory. It regulates the behavior of $L_{n}-L$ near $\beta^{*}$ and is related to the rate of convergence (see, Ahidar-Coutrix et al., 2020, Section 3). In addition, if $\mathcal{M}=\mathbb{R}^{D}$, then Assumption 4 is reduced to the usual condition that $\mathbb{E}[(X-\mu^{*})(X-\mu^{*})^{\top}]$ is positive-definite with $\lambda_{X}$ being the smallest eigenvalue of $\mathbb{E}[(X-\mu^{*})(X-\mu^{*})^{\top}]$. Therefore, Assumption 4 generalizes this basic distributional condition from Euclidean space to general metric spaces. Examples of $X$ satisfying Assumption 4 are given in Section LABEL:subsec:ex-as3-sphere of Supplementary Material.

Remark 3*.*

Although certain conditions of Assumption 3 may be satisfied for certain metric spaces according to the above discussions, we shall highlight that both Assumption 3 and Assumption 4 generally relate to the distribution of $(X,Y)$ and are not solely dependent on the geometry of $\mathcal{M}$, even in the special case of Euclidean space. Hence, it is generally not possible to specify which spaces meet these two assumptions without considering the underlying distributions.

Remark 4*.*

We do not assume the space $\mathcal{M}$ to be bounded. In this general setting, the bounded support or similar conditions on $X$ are typically assumed to establish the consistency and convergence rates of the MLE $\hat{\beta}_{n}$ even for the Euclidean logistic regression model (4) (Fahrmeir and Kaufmann, 1985, 1986; Yin et al., 2006; van de Geer, 2008; Meier et al., 2008; Bach, 2010; Liang and Du, 2012; Negahban et al., 2012; Ning and Liu, 2017; Guo et al., 2021). This boundedness condition on the data support, as one might expect, could potentially be relaxed to allow the data space to expand at a certain rate as the sample size grows; a similar relaxation in Euclidean spaces has been explored, for example, by Fahrmeir and Kaufmann (1986). Alternatively, it may be replaced with certain conditions (e.g., sub-Gaussianity) on the tail probability of the distribution of $X$. Since such relaxations involve much heavier technicalities without providing additional insights, we have chosen to omit them in this paper to maintain a more focused presentation.

Remark 5*.*

The constraint $m\geq 4$ in Assumption 4 is imposed to obtain the rate of convergence for the expected squared estimation error \mathbb{E}\big{[}d^{2}(\beta^{*},\hat{\beta}_{n})\big{]}; for the expected estimation error \mathbb{E}\big{[}d(\beta^{*},\hat{\beta}_{n})\big{]}, $m\geq 2$ is sufficient. Especially, if $\mathcal{M}$ is bounded, the constraint $\alpha_{L}>1+2/m$ in Assumption 4 can be reduced to $\alpha_{L}>1$, since now $\mathbb{E}[d^{m}(\mu^{*},X)]<\infty$ for any arbitrarily large but fixed $m$. In addition, $\alpha_{L}>1$ is sufficient for establishing the convergence rate of $d(\beta^{*},\hat{\beta}_{n})$ in probability.

In the following subsections we present convergence rates of $\hat{\beta}_{n}$ under two types of metric entropy conditions that are widely used in the literature to quantify space complexity.

Notation. For two non-negative sequences $\{a_{n}\}$ and $\{b_{n}\}$, we write $a_{n}\lesssim b_{n}$ (respectively, $a_{n}\gtrsim b_{n}$) if there is a constant $c>0$ not depending on $n$, such that $a_{n}\leq cb_{n}$ (respectively, $a_{n}\geq cb_{n}$) for all sufficiently large $n$. We write $a_{n}\asymp b_{n}$ if and only if both $a_{n}\lesssim b_{n}$ and $a_{n}\gtrsim b_{n}$ are true.

4.2 Rates of Convergence and Optimality with Local Metric Entropy Condition

The first type of the condition on metric entropy, termed the (local) metric entropy condition in this paper, concerns the metric entropy in a local ball of $\mathcal{M}$. Let $B(\beta^{*},\delta)\subset\mathcal{M}$ be the (closed) ball of the radius $\delta>0$ centered at $\beta^{*}$ and $N\left(\epsilon,B(\beta^{*},\delta),d\right)$ be its covering number (Wainwright, 2019, Section 5.1) using balls of the radius $\epsilon>0$. Metric entropy then refers to $\log N\left(\epsilon,B(\beta^{*},\delta),d\right)$. We consider the following two metric entropy conditions of different scales.

Assumption 5*.*

There is a constant $D>0$, such that one of the following holds:

1. (a)

(Log-polynomial metric entropy) $N\left(\epsilon,B(\beta^{*},\delta),d\right)\lesssim\left(\frac{\delta}{\epsilon}\right)^{D}$, 2. (b)

(Polynomial metric entropy) $\log N\left(\epsilon,B(\beta^{*},\delta),d\right)\lesssim\left(\frac{\delta}{\epsilon}\right)^{D}$.

Assumption 5(a), in which the metric entropy $\log N\left(\epsilon,B(\beta^{*},\delta),d\right)$ is of the order $D\log(\delta/\epsilon)$, is often satisfied by finite-dimensional spaces, such as the spaces (a)–(c) in Example Example. More generally, Assumption 5(a) holds for $D$-dimensional Riemannian manifolds of bounded curvature and metric spaces that admit an Ahlfors–David regular measure (see Example 2.3, Ahidar-Coutrix et al., 2020).

Assumption 5(b), when compared to 5(a), allows the covering number to grow in a much faster rate and thus accommodates spaces of much higher complexity, such as the Wasserstein space in Example Example(d). In fact, for this Wasserstein space, one can show $\log N\left(\epsilon,B(\beta^{*},\delta),d_{W}\right)\asymp{\delta}/{\epsilon}$, where $d_{W}$ is the $2$-Wasserstein distance, based on the metric entropy of the set of all non-decreasingly bounded functions on $[0,1]$ (Yang, 1999) and by following the arguments in the proof of Proposition 1 in Petersen and Müller (2019). Thus, this Wasserstein space satisfies Assumption 5(b) with $D=1$, but not Assumption 5(a).

Define

[TABLE]

when $\mathcal{M}$ satisfies Assumption 5(a), and

[TABLE]

when $\mathcal{M}$ satisfies Assumption 5(b). With this definition, the following theorems provide upper and lower bounds on the rate of convergence for the estimator $\hat{\beta}_{n}$ under the metric entropy assumption.

Theorem 1**.**

Suppose that Assumptions 1–5 hold. Then, there is a constant $c>0$ depending on the constants in the assumptions, such that for any measurable sequence $\hat{\beta}_{n}$ indexed by $n$ and determined by (9),

[TABLE]

Theorem 2**.**

Suppose $\sigma_{X}^{2}=\mathbb{E}[d^{2}(\mu^{*},X)]<\infty$ and $\mathcal{M}$ satisfies Assumptions 1 and 2. If $(X_{i},Y_{i}),i=1,\ldots,n$, are i.i.d. sampled from the model (6), then

[TABLE]

In the above, the infimum is taken over all estimators $\hat{\beta}_{n}$, and $\mathbb{E}_{\beta}$ is the expectation with respect to the distribution induced by the model (6) with $\beta^{*}=\beta$.

According to Theorem 1 and Theorem 2, for the commonly encountered scenario where $\alpha_{L}=2\alpha_{U}$, the proposed estimator $\hat{\beta}_{n}$ is rate-optimal under Assumption 5(a) of the log-polynomial metric entropy condition. It is also optimal under Assumption 5(b) of the polynomial metric entropy condition with $D\leq 2\alpha_{U}$ (up to a logarithmic factor when $D=2\alpha_{U}$). To the best of our knowledge, this is the first optimality result in parametric binary regression with covariates residing in a general metric space.

Among nonlinear metric spaces, Riemannian manifolds, which satisfy Assumption 5(a), are most commonly seen in statistics and have a wide range of applications (e.g., Pennec et al., 2006; Dryden et al., 2009). As an application of Theorems 1 and 2, the following corollary demonstrates that, for certain Riemannian manifolds, the proposed estimator is rate-optimal in estimation. Recall that a Hadamard manifold is a simply connected and complete Riemannian manifold that is also a Hadamard space. Let $h(\beta^{*})=h(\beta^{*};X,\mu^{*})$ be viewed as a function of $\beta^{*}$ and let $\Lambda_{\min}(\cdot)$ denote the smallest eigenvalue of a linear operator. Recall the Hessian defined in (3).

Corollary 1**.**

For a $D$-dimensional Hadamard manifold $\mathcal{M}$ of bounded sectional curvature, if the distribution of $X$ has a bounded support, the smallest eigenvalue $\lambda_{X}$ of $\mathbb{E}[\text{Hess}_{\beta^{*}}h]$ satisfies $\lambda_{X}>0$, and $\hat{\beta}_{n}$ exists almost surely, then

[TABLE]

The same conclusion holds for a complete and simply connected $D$-dimensional Riemannian manifold $\mathcal{M}$ with non-negative sectional curvature upper bounded by $\kappa>0$ and of diameter no larger than $D_{\kappa}/2$ such that $\inf_{p,q\in\mathcal{M}}\Lambda_{\min}(\text{Hess}_{q}d_{p}^{2})>0$, where $d_{p}^{2}(\cdot)=d^{2}(\cdot,p)$.

A proof of this corollary follows directly from verifying Assumptions 1–5, which is presented in the proof of Theorem 6. Examples of bounded and unbounded Riemannian manifolds that satisfy the requirement for the Hessian of $h$ are provided in Section LABEL:subsec:ex-as3-sphere of Supplementary Material.

4.3 Rates of Convergence and Optimality with Global Metric Entropy Assumption

In the case that $\mathcal{M}$ is totally bounded, the global metric entropy condition (Yang and Barron, 1999), concerning the minimal number of balls to cover the entire metric space $\mathcal{M}$, is another widely adopted condition on metric entropy.

Assumption 6*.*

There is a constant $D>0$, such that one of the following holds:

1. (a)

(Log-polynomial global metric entropy) $N\left(\epsilon,\mathcal{M},d\right)\lesssim\left(\frac{1}{\epsilon}\right)^{D}$, 2. (b)

(Polynomial global metric entropy) $\log N\left(\epsilon,\mathcal{M},d\right)\lesssim\left(\frac{1}{\epsilon}\right)^{D}$.

Assumption 6 can be deduced from Assumption 5 when $\mathcal{M}$ is bounded, and hence is weaker than the latter. This may make Assumption 6 relatively easier to check. All finite-dimensional compact Riemannian manifolds satisfy Assumption 6(a). A collection of $D$-dimensional Riemannian manifolds embedded in a common compact subspace of $\mathbb{R}^{D}$, when endowed with the Hausdorff distance, satisfies Assumption 6(b) (Theorem 9, Genovese et al., 2012). Assumption 6 is widely adopted in the literature of classification and/or metric-space data analysis (Ahidar-Coutrix et al., 2020; Yang, 1999; Genovese et al., 2012; Lin and Müller, 2021, among others). For instance, Ahidar-Coutrix et al. (2020) derived upper bounds on the rate of convergence of the generalized Fréchet mean in metric spaces under Assumption 6(b).

Let $\gamma_{0}=W_{0}(1/2)\approx 0.35173$, where $W_{0}$ is the Lambert W function. With Assumption 6, define

[TABLE]

when $\mathcal{M}$ satisfies Assumption 6(a), and

[TABLE]

when $\mathcal{M}$ satisfies Assumption 6(b). The following theorems provide upper and lower bounds on \mathbb{E}\big{[}d^{2}(\beta^{*},\hat{\beta}_{n})\big{]} under each condition of Assumption 6.

Theorem 3**.**

Suppose that Assumptions 1–4 and 6 hold. Then, there is a constant $c>0$ depending on the constants in the assumptions, for any measurable sequence $\hat{\beta}_{n}$ indexed by $n$ and determined by (9), such that

[TABLE]

Theorem 4**.**

Suppose $\sigma_{X}^{2}=\mathbb{E}[d^{2}(\mu^{*},X)]<\infty$ and $\mathcal{M}$ satisfies Assumptions 1 and 2. If $(X_{i},Y_{i}),i=1,\ldots,n$, are i.i.d. sampled from the model (6), then

[TABLE]

In the above, the infimum is taken over all estimators $\hat{\beta}_{n}$ and $\mathbb{E}_{\beta}$ is the expectation with respect to the distribution induced by the model (6) with $\beta^{*}=\beta$.

Theorems 3 and 4 show that, if $\alpha_{L}=2\alpha_{U}$, then the estimator $\hat{\beta}_{n}$ is rate-optimal under the polynomial global metric entropy condition with $D\leq 2\alpha_{U}$ (up to a logarithmic factor when $D=2\alpha_{U}$). Under the log-polynomial global metric entropy condition, there is a small gap $\gamma_{0}$ between the upper and lower bounds. Similar gaps are also observed in other estimation problems (Yang and Barron, 1999; Ma and Wu, 2015). Closing such a gap is challenging especially for nonlinear metric spaces and left for future study.

5 Logistic Classifiers

Plugging the estimators $\hat{\beta}_{n}$ and $\hat{\mu}_{n}$ into the model (6), we can estimate the success probability $\mathbb{P}(Y=1\mid X=x)$ and obtain a classifier by classifying the instance $x\in\mathcal{M}$ accordingly. Specifically, with $\hat{P}_{X}=\{1+\exp(-h(\hat{\beta}_{n};X,\hat{\mu}_{n}))\}^{-1}$ denoting the estimated success probability, we propose the classifier

[TABLE]

Of theoretical interest is the rate of the excess risk of $\hat{s}$ relative to the Bayes classifier. For a random pair $(X,Y)\in\mathcal{M}\times\{0,1\}$ sampled from the model (6), the Bayes classifier $s^{*}$, minimizing the error probability, is

[TABLE]

The excess risk of $\hat{s}$ is defined by $\mathcal{E}(s^{*},\hat{s})=\mathbb{E}\mathcal{L}(\hat{s})-\mathbb{E}\mathcal{L}(s^{*})$, where $\mathcal{L}(s)=\mathbb{I}_{\{s(X)\neq Y\}}$ signifies whether an error is made by the classifier $s$ and $\mathbb{I}_{\{\cdot\}}$ is the indicator function. The excess risk, commonly investigated in the literature of classification (Massart and Nédélec, 2006; Audibert and Tsybakov, 2007), quantifies the classification accuracy of $\hat{s}$ (relative to the Bayes classifier $s^{\ast}$), since $\mathbb{P}\{\hat{s}(X)=Y\}=1-\mathbb{E}\mathcal{L}(\hat{s})=(1-\mathbb{E}\mathcal{L}(s^{*}))-\mathcal{E}(\hat{s},s^{*})$, where $1-\mathbb{E}\mathcal{L}(s^{*})$ is the classification accuracy of the Bayes classifier $s^{*}$.

We first provide an upper bound on the excess risk in terms of the estimation errors of $\hat{\beta}_{n}$ and $\hat{\mu}_{n}$ for Alexandrov spaces. This bound requires some regularity on the Alexandrov inner product $h(\beta^{*};x,\mu)$ defined in (6). We say $h(\beta^{*};x,\mu)$ is regular near $\mu^{*}$ if there exist constants $C_{1},C_{2}>0$ and arbitrarily small but fixed constants $0<r_{1}\leq r_{2}$ such that, for all $\mu\in B(\mu^{*},r_{2})$, the following two Lipschitz continuity conditions are satisfied:

$|h(\beta^{*};x,\mu)-h(\beta^{*};x,\mu^{*})|\leq C_{2}d(\mu,\mu^{*})$ for all $x\in B(\mu^{*},r_{2})$, and

• 2)

$|d(\beta^{*},\mu)\cos(\angle_{\mu}(x,\beta^{*}))-d(\beta^{*},\mu^{*})\cos(\angle_{\mu^{*}}(x,\beta^{*}))|\leq C_{1}d(\mu,\mu^{*})$ for all $x$ with $d(x,\mu^{*})=r_{1}$.

This mild regularity condition, guaranteeing adequate (Lipschitz) continuity of the Alexandrov inner product $h(\beta^{*};x,\mu)$ near $\mu^{*}$, holds for all $\mu^{*}$ in Riemannian manifolds (e.g., the SPD matrix manifold and the unit sphere in Example Example) since $h(\beta^{*};x,\mu)$ is smooth in both $x$ and $\mu$, and for the Wasserstein space in Example Example(d). It also generally holds for other metric spaces when $\mu^{*}$ is a “nonsingular” point, e.g., any interior point of an orthant in Example Example(c) of the Billera–Holmes–Vogtmann tree space (Billera et al., 2001).

Remark 6*.*

The constants $C_{1}$ and $C_{2}$ in the above typically depend on the constants $r_{1}$ and $r_{2}$, respectively. For instance, in the special case of Euclidean space (i.e., when $\mathcal{M}=\mathbb{R}^{D}$), $C_{1}$ diverges as $r_{1}$ approaches 0, while $C_{2}$ grows as $r_{2}$ approaches $\infty$. Therefore, the above regularity of $h(\beta^{*};x,\mu)$ can only be assumed locally near $\mu^{*}$ (but not globally in the entire space $\mathcal{M}$); this significantly intensifies the complexity of the proofs (e.g., see Corollary LABEL:cor:fmse and Lemma LABEL:lem:mid-seg-comparison-model-space of Supplementary Material).

The following theorem upper bounds the excess risk in terms of the estimation error of $\hat{\beta}_{n}$ and $\hat{\mu}_{n}$, where we recall $D_{k}=\pi/\sqrt{\kappa}$ if $\kappa>0$ and $D_{\kappa}=\infty$ otherwise.

Theorem 5**.**

Suppose $\mathcal{M}$ is a complete Alexandrov space with curvature upper bounded by $\kappa$ and of diameter no larger than $D_{\kappa}/2$. Assume that the Alexandrov inner product $h(\beta^{*};x,\mu)$ is regular near $\mu^{*}$. Then, under Assumptions 3, there exists a constant $c>0$ such that

[TABLE]

The assumed space in the above theorem includes Hadamard spaces with $\kappa=0$ ($D_{\kappa}=\infty$ in this case) and positively curved spaces such as spheres with $\kappa>0$. The convergence rate of $\hat{\mu}_{n}$ established in Schötz (2019); Ahidar-Coutrix et al. (2020) and the rate of $\hat{\beta}_{n}$ in Theorems 1 and 3 imply $\mathbb{E}\left[d^{2}(\mu^{*},\hat{\mu}_{n})\right]\lesssim\mathbb{E}[d^{2}(\beta^{*},\hat{\beta}_{n})]$. Consequently, the above theorem, together with Theorems 1 and 3, provides upper bounds on the rate of the excess risk under the log-polynomial and polynomial metric entropy conditions.

For Riemannian manifolds, we can obtain a finer upper bound on the excess risk by exploiting the Riemannian structure. Let $h(\beta^{*})=h(\beta^{*};X,\mu^{*})$ be viewed as a function of $\beta^{*}$ and let $\Lambda_{\min}(\cdot)$ denote the smallest eigenvalue of a linear operator. Recall the Hessian defined in (3).

Theorem 6**.**

For a $D$-dimensional Hadamard manifold of bounded sectional curvature, suppose that the distribution of $X$ has a bounded support, the smallest eigenvalue $\lambda_{X}$ of $\mathbb{E}[\text{Hess}_{\beta^{*}}h]$ satisfies $\lambda_{X}>0$, and $\hat{\beta}_{n}$ exists almost surely. Then, there is a constant $c>0$ depending on $\lambda_{X}$ and $\sigma_{X}^{2}=\mathbb{E}[d^{2}(\mu^{*},X)]$, such that

[TABLE]

The same conclusion holds for a complete and simply connected $D$-dimensional Riemannian manifold $\mathcal{M}$ with non-negative sectional curvature upper bounded by $\kappa>0$ and of diameter no larger than $D_{\kappa}/2$ such that $\inf_{p,q\in\mathcal{M}}\Lambda_{\min}(\text{Hess}_{q}d_{p}^{2})>0$, where $d_{p}^{2}(\cdot)=d^{2}(\cdot,p)$.

Below we provide a matching lower bound for the class $\mathcal{C}=\{s(X)=\mathbb{I}_{\{h(\beta;X,\mu^{*})\geq 0\}}:\beta\in\mathcal{M}\}$ for Riemannian manifolds, establishing the asymptotic minimax optimality of the proposed classifier in terms of both $n$ and the dimension $D$; a similar class of classifiers was considered by Abramovich and Grinshtein (2019) in the Euclidean setting. To the best of our knowledge, this is the first optimality result for classifiers on covariates without a vector structure.

Theorem 7**.**

For a $D$-dimensional Riemannian manifold $(\mathcal{M},d)$, if $(X_{i},Y_{i}),i=1,\ldots,n$, are i.i.d. drawn from the model (6), then

[TABLE]

where the infimum is taken over all classifiers.

Remark 7*.*

For general metric spaces, Gottlieb et al. (2016) provides a lower bound of the rate $n^{-1}$ (without a matching upper bound) under a restrictive condition that the Bayes classifier has zero error probability. This lower bound, when applied to Riemannian manifolds, is sub-optimal, as it requires the restrictive condition and does not include the dimension $D$. Finding matching lower bounds and establishing optimality for the excess risk in general metric spaces are much more challenging and left for future exploration.

6 Simulation Studies

We illustrate the numeric performance of the proposed regression from the perspectives of both estimation and classification by using the manifold $\mathcal{M}={\mathcal{S}_{3}^{+}}$ of $3\times 3$ symmetric positive-definite (SPD) matrices endowed with the Log-Cholesky metric (Lin, 2019). For estimation, we consider

• •

Case 1: data are sampled from the proposed model (6).

For the Fréchet mean $\mu^{*}$, we consider two settings, namely, the $3\times 3$ identity matrix $I_{3}$ and the matrix $\Sigma_{AR}=(0.5^{|i-j|})_{1\leq i,j\leq 3}$. For the coefficient $\beta^{*}$, we also consider two settings, namely, $3I_{3}$ and $3\Sigma_{AR}$. For classification, we consider the following additional data models:

• •

Case 2: $\text{logit}(P_{X})=2\sin(\pi h(3\Sigma_{AR};X,I_{3}))+h(3\Sigma_{AR};X,I_{3})$,

• •

Case 3: $\text{logit}(P_{X})=(\log 3)\sin(\pi\vec{X}_{1})+\frac{\sqrt{3}}{2}\vec{X}_{2}^{2}+\frac{\sqrt{3}}{4}\exp(\vec{X}_{3})+\frac{3}{4}\vec{X}_{5}+2(\log\frac{3}{2})(\vec{X}_{4}+\vec{X}_{6})$,

where $P_{X}$ is the conditional success probability, $h$ is defined in (6) and $\vec{X}=(\vec{X}_{1},\ldots,\vec{X}_{6})$ is the vector representation of $X\in{\mathcal{S}_{3}^{+}}$, i.e., the $6$-dimensional vector obtained by stacking elements in the lower triangular part of $X$. The model in Case 2 includes a nonlinear term of $h$, while the model in Case 3 is nonparametric and highly nonlinear in the covariate; both are designed to challenge the proposed classifier.

For each of the above models, we synthesize data under different settings. Specifically, we consider two sample sizes, $n=100$ and $n=500$. The data are synthesized, as follows. For each $i=1,\ldots,n$, we draw the covariate $X_{i}=\mathrm{Exp}_{\mu^{*}}S_{i}$ with ${\vec{S_{i}}}\sim N(0,rI_{6})$, where $r=1,4$ determines the variance $\sigma_{X}^{2}=\mathbb{E}[d^{2}(\mu^{*},X_{i})]$ and $\mathrm{Exp}_{\mu^{*}}$ denotes the Riemannian exponential map at $\mu^{*}$. Finally, we compute $P_{X_{i}}$ according to the models and sample $Y_{i}$ via Bernoulli trial with the success probability $P_{X_{i}}$. For each model and each setting, we perform $500$ independent Monte Carlo replicates. In each replicate, we split the dataset into a training set $\mathcal{D}_{\mathrm{train}}$ (80%) and a test set $\mathcal{D}_{\mathrm{test}}$ (20%).

Estimation quality. As mentioned in the above, we assess the quality of the proposed estimators $\hat{\mu}_{n}$ and $\hat{\beta}_{n}$ computed from $\mathcal{D}_{\mathrm{train}}$ in Case 1, and quantify estimation quality by $d(\mu^{*},\hat{\mu}_{n})$ and $d(\beta^{*},\hat{\beta}_{n})$, respectively. In addition, based on $\hat{\mu}_{n}$ and $\hat{\beta}_{n}$, for each data point $i\in\mathcal{D}_{\mathrm{test}}$ in the test set, we compute the estimated success probability $\hat{P}_{X_{i}}$ which is of interest in practice. To quantify the estimation quality of the success probability, we consider the root-mean-squared error (RMSE)

[TABLE]

Tables 1 summarizes the results of the simulation studies with $\mu^{*}=I_{3}$ and $\beta^{*}=3\Sigma_{AR}$. We see that, the estimation errors and RMSE decrease as the sample size grows, as indicated by our theoretical results in Section 4. We also observe that the error $d(\beta^{*},\hat{\beta}_{n})$ decreases as $\sigma_{X}^{2}$ grows. Intuitively, this is because, from the perspective of estimating $\beta^{*}$, a larger $\sigma_{X}^{2}$ indicates a larger signal-to-noise ratio, making estimation of the coefficient $\beta^{*}$ relatively easier; see Example 8.4 in Duchi (2019) for a similar observation in the Euclidean space. This phenomenon is also dictated by Theorems 2 and 4, where the lower bounds decrease with respect to $\sigma_{X}^{2}$. In contrast, the error $d(\mu^{*},\hat{\mu}_{n})$ increases as $\sigma_{X}^{2}$ increases. This is because, like estimating the mean of a Gaussian sample in the Euclidean space, the error of estimating the Fréchet mean $\mu^{*}$ is scaled up according to the variance $\sigma_{X}^{2}$ of $X$. Similar results are observed for the other settings presented in Section LABEL:sec:further-numerical of Supplementary Material.

Classification performance. We compare the proposed classifier with the SPDNet (Huang and Van Gool, 2017), $k$-NN (Chaudhuri and Dasgupta, 2014), kernel support vector machine (kSVM, Jayasumana et al., 2013) and kernel density classifier (Section 6.6.2, Hastie et al., 2009). The details of these methods on SPD matrices are given in Section LABEL:sec:further-numerical of Supplementary Material. For evaluation, we consider all models in the above, and measure classification performance by accuracy, area under the receiver operating characteristic curve (AUC), sensitivity, specificity and empirical excess risk (EER) defined by

[TABLE]

where $\hat{s}$ is a classifier, $s^{*}$ is the Bayes classifier, and $\mathcal{L}(\cdot)$ is defined in Section 5.

Table 2 summarizes the results with $\sigma_{X}^{2}=2.248$ (corresponding to $r=1$) and with $\mu^{*}=I_{3}$, $\beta^{*}=3\Sigma_{AR}$ for Case 1; similar results for other settings are presented in Section LABEL:sec:further-numerical of Supplementary Material. First, in Case 1 where the data are sampled from the proposed model, the EER of the proposed classifier decreases as the sample size grows, which numerically illustrates the theoretical results in Section 5. Also, as expected, compared with other methods, the proposed classifier has superior performance in Case 1. For Case 2, although the data model contains a nonlinear term of the Alexandrov inner product, the proposed classifier still outperforms the other methods. In Case 3 that strongly favors nonparametric methods, the proposed classifier can produce reasonable but suboptimal results relative to the other classifiers. In summary, there is no universally best method, and the proposed classifier has excellent performance when the data model is linear or moderately nonlinear while the other classifiers generally produce better results when the data model is highly nonlinear and the manifold has relatively low dimensions.

7 Data Application

We apply the proposed method to the task-based fMRI data from the Human Connectome Project (Van Essen et al., 2013). In the project, each subject was asked to complete several tasks that are related to certain human functions. For instance, in the motor task, subjects were asked to tap their fingers, squeeze their toes, or move their tongue with the given visual cues, while in the language processing task, subjects were asked to answer related questions after brief auditory stories, or complete addition and subtraction math problems which are presented auditorily. During each task, the temporal blood-oxygen-level dependent (BOLD) signals from each brain region were recorded. Based on the BOLD signals, functional connectivity of the brain, represented by covariance matrices of the BOLD signals, can be derived. Of interest is the effect of functional connectivity on the motor and language processing tasks.

We focus on $201$ subjects performing only the motor task and $193$ subject performing only the language processing task; data111Available at https://www.humanconnectome.org/study/hcp-young-adult/document/500-subjects-data-release. collected from these subjects form a dataset of the sample size $n=394$. We consider $m=8$ brain regions related to motor skills and language processing, and apply a preprocessing pipeline described in Section LABEL:sec:preprocess-application of Supplementary Material to obtain the observations $(X_{i},Y_{i})\in{\mathcal{S}_{8}^{+}}\times\{0,1\}$ for $i=1,\ldots,n$, where $X_{i}$ is the functional connectivity, $Y_{i}=0$ for the subject completing the motor task and $Y_{i}=1$ for the subject completing the language processing task. As the Euclidean distance on ${\mathcal{S}_{8}^{+}}$ severely suffers from the “swelling effect” that artificially inflates the variation in statistical analysis of diffusion tensors and functional connectivities (Arsigny et al., 2007), to avoid this undesirable effect, we adopt the Log-Cholesky distance (Lin, 2019) on ${\mathcal{S}_{8}^{+}}$ specifically designed for eliminating the swelling effect; the Log-Euclidean (Arsigny et al., 2007) and affine-invariant (Moakher, 2005) distances yield similar results that are thus not presented.

We fit the proposed model (6) by using the above data and depict the estimated parameters in Figure 2, where we plot the geodesic $\gamma_{\hat{\mu}_{n}}^{\hat{\beta}_{n}}$ emanating from $\hat{\mu}_{n}$ in both directions. To assess the covariate effect, we test the null hypothesis $H_{0}:\beta^{*}=\mu^{*}$, which corresponds to no different effect of functional connectivity on the two tasks, by using the permutation test outlined in Section LABEL:subsec:perm-test of Supplementary Material. The empirical $p$-value is less than 0.001, strongly suggesting that functional connectivity has significantly different effect on motor skills and language processing.

The geodesic $\gamma_{\hat{\mu}_{n}}^{\hat{\beta}_{n}}$ in Figure 2 is normal to the empirical decision boundary between the two tasks, with the positive direction (larger odds) extending into the language processing task and the negative direction (smaller odds) extending into the motor task. From the figure we also find that, for any $t$, if the functional connectivity is changed from $\gamma_{\hat{\mu}_{n}}^{\hat{\beta}_{n}}(t)$ to $\gamma_{\hat{\mu}_{n}}^{\hat{\beta}_{n}}(t+\delta)$ with $\delta=0.00375$ (i.e., the functional connectivity advances along the direction $\gamma_{\hat{\mu}_{n}}^{\hat{\beta}_{n}}$ by $\delta$ unit), then the change in the log odds is $1.486$. Moreover, we observe that the odds becomes much larger than 1 as Region 1 (the superior temporal gyrus) becomes more active, which suggests that Region 1 is substantially more involved in auditory and language processing; the finding matches the previous study on this brain region (Bigler et al., 2007).

In addition, we observe that the correlation between Region 1 and other regions exhibits different patterns in the two tasks. For example, the correlation between Region 1 and Region 4 (which is part of the inferior frontal gyrus related to speech and language processing) in the direction of larger odds is visibly stronger than that in the opposite direction, suggesting a stronger functional connectivity between Region 1 and Region 4 in the language processing task. In contrast, the correlation between Region 1 and Regions 5–8 in the direction of larger odds is weaker than that in the opposite direction. This may indicate that, in the language processing task which does not invoke motor and visual processing, the functional connectivity between Region 1 and Regions 5–8 which include the primary motor cortex and visual cortex may be suppressed to certain extent.

The in-sample ROC (receiver operating characteristic) curve in Figure 3 shows that the proposed model is adequate to classify the two tasks based on the functional connectivity. In addition, as in Section 6, we compare the proposed classifier with SPDNet, $k$-NN, kernel SVM (kSVM) and kernel density classifier (KDC). The $10$-fold cross-validated accuracy, AUC, sensitivity and specificity are presented in Table 3. From the table, it can be seen that, even compared with SPDNet that is based on neural networks, the proposed classifier has slightly higher accuracy, and outperforms the other methods, perhaps because some nonparametric classifiers start to suffer from the curse of dimensionality as the manifold ${\mathcal{S}_{8}^{+}}$ has 36 dimensions. Moreover, the proposed classifier is more interpretable. For example, as discussed in the above, the geodesic $\gamma_{\hat{\mu}_{n}}^{\hat{\beta}_{n}}$ in Figure 2, which is the direction normal to the empirical decision boundary between the two tasks, reveals interesting patterns of functional connectivity in these two tasks.

The aforementioned different patterns of inter-connectivities (the off-diagonal entries of each $X_{i}\in{\mathcal{S}_{8}^{+}}$) also play a significant role in classifying the two tasks. To illustrate this, we fit a Euclidean logistic regression model (denoted by Logistic-BOLD in Table 3) that includes only the BOLD signals (diagonal entries of each $X_{i}\in{\mathcal{S}_{8}^{+}}$) as covariates. Compared with the model (6), this naïve logistic model, which does not consider inter-connectivities, has substantially lower accuracy and other deteriorated performance measures.

In summary, the proposed model can uncover interesting and interpretable effects of functional connectivity on the motor and language processing tasks, while offering near state-of-the-art performance in classifying the two tasks.

8 Concluding Remarks

The proposed model (6), which incorporates covariates located in a metric space without a vector structure, offers an intuitive approach to modeling the effect of metric-space valued covariates on the log odds, and the model parameters enable practitioners to gain insights into the fundamental distinctions between the two classes under consideration. While we have established the first set of optimality results for regression with a binary response and covariates on general metric spaces, several open problems remain for future studies. For example, optimality in the case of $D>2\alpha_{U}$ remains largely unknown. In addition, under the log-polynomial global metric entropy condition, there is a small gap between the upper and lower bounds as presented in Theorems 3 and 4, respectively. Closing this gap is challenging and, to the best of our knowledge, remains unresolved even in the Euclidean setting. Furthermore, finding the minimax lower bound on the excess risk for general metric spaces awaits future exploration.

Acknowledgments

Z. Lin’s research is partly supported by the NUS startup grant A-0004816-01-00. Data were provided in part by the Human Connectome Project, WU-Minn Consortium (PI: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University.