# A Type of Nonlinear Fréchet Regressions

Lu Lin and Ze Chen
Zhongtai Securities Institute for Financial Studies, Shandong University,
Jinan, China
The corresponding author. Email<EMAIL_ADDRESS>Lin’s research was
supported by National Natural Science Foundation (NNSF) of China (No.11971265)
and the National Statistical Science Research Project (2022LD03). Chen’s
research was supported by the Postdoctoral Fellowship Program of CPSF
(GZC20231478).

###### Abstract

The existing Fréchet regression is actually defined within a linear framework,
since the weight function in the Fréchet objective function is linearly
defined, and the resulting Fréchet regression function is identified to be a
linear model when the random object belongs to a Hilbert space. Even for
nonparametric and semiparametric Fréchet regressions, which are usually
nonlinear, the existing methods handle them by local linear (or local
polynomial) technique, and the resulting Fréchet regressions are (locally)
linear as well. We in this paper introduce a type of nonlinear Fréchet
regressions. Such a framework can be utilized to fit the essentially nonlinear
models in a general metric space and uniquely identify the nonlinear structure
in a Hilbert space. Particularly, its generalized linear form can return to
the standard linear Fréchet regression through a special choice of the weight
function. Moreover, the generalized linear form possesses methodological and
computational simplicity because the Euclidean variable and the metric space
element are completely separable. The favorable theoretical properties (e.g.
the estimation consistency and presentation theorem) of the nonlinear Fréchet
regressions are established systemically. The comprehensive simulation studies
and a human mortality data analysis demonstrate that the new strategy is
significantly better than the competitors.

Keywords: Metric space, Hilbert space, Fréchet regression, nonlinearity,
generalized linearity.

##  1 Introduction

###  1.1 Related issues and our goals

When data types are out of the Euclidean framework and lack a general
algebraic structure, the methods of metric space and Hilbert space in
probability and statistics offer alternative tools for analyzing and modeling
these abstract contexts. These methodologies shift the starting point away
from the Euclidean sample space, and replace it with a metric space or a
Hilbert space according to the features of the abstract and complex elements.
The early work of metric space methods in probability and statistics can date
back at least to the Fréchet mean in Fréchet (1948), which is defined for
random elements in a metric space as the minimizer of an expected loss defined
within that space. Similarly, the Hilbert space methods employ notions such as
inner product, orthogonality, and projection within linear subspaces, as
alternative approaches to statistical inference on abstract and complex
elements (see Small and Mcleish (1994) for the systematic introduction).

In recent years, the increasing complexity and abstraction of data types
driven by application demands, the advancements in measurement techniques and
improvements in data collection and storage, have sparked significant interest
in statistical methodologies for analyzing random objects in non-Euclidean
spaces (see, e.g., Marron and Alonso (2014), Wang and Marron (2007) and
Faraway (2014)). As a prominent issue in non-Euclidean space, the Fréchet
regression was introduced by Petersen and Müller (2019), in which the response
is a random object defined in a metric space. Specifically, the Fréchet
regression function is defined as the minimizer of an expected weighted loss
defined within a metric space. The variable selection and model averaging for
Fréchet regression were further proposed by Tucker et al. (2023) and Yan et
al. (2023), respectively. Moreover, the Fréchet regression was recently
extended into semiparametric models by Ghosal and Meiring (2023) and
Bhattacharjee and Müller (2023).

The existing Fréchet regressions, however, are actually defined within a
linear framework. This is because the weight function in the Fréchet
regression objective function is linearly defined, and the resulting Fréchet
regression function is only suitable for the linear model. Even for
nonparametric and semiparametric Fréchet regressions, which are often
nonlinear, current methods handle them by local linear (or local polynomial)
methods or approximately linear techniques, the resulting Fréchet regressions
are locally or approximately linear as well (see, e.g., Petersen and Müller
(2019), Schötz (2022), Ghosal and Meiring (2023) and Bhattacharjee and Müller
(2023)).

To the best of our knowledge, the pursuit of nonlinear Fréchet regression
remains an open question. Constructing the framework of nonlinear Fréchet
regression may be highly nontrivial. The main difficulties are as follows. In
addition to that the framework should possess the common statistical
properties, it should exhibit adaptability in the following aspects:

  1. (1)

it can correctly and even explicitly identify the underlying nonlinear
regression function when the underlying regression is definitely nonlinear;

  2. (2)

it should encompass a wide range of nonlinear functions, and yet return to a
standard linear framework when the underlying regression is indeed linear;

  3. (3)

it is easily implementable and computationally efficient especially for the
case when the estimation procedure for the unknown parameters in a nonlinear
function cannot be separated from the metric space.

The challenges mainly come from the feature that the random object in a metric
does not live in a linear space, thereby rendering common mathematical
operations such as addition and multiplication inapplicable to these objects.
This feature results in a major problem for applying parameter based methods
to define the nonlinear Fréchet regression and estimate the unknown
parameters. These challenges will be successively investigated in the
following sections.

Besides, like the common regression in Euclidean space, the Fréchet regression
is originally defined by a conditional expected loss (Petersen and Müller,
2019). The classic approach to tackling the conditional expectation problem
involves transforming it into a system of unconditional expectations. Then,
the unconditional expectation methods, such as the optimally weighted
generalized method of moments, can be employed to realize statistical
inference. This transformation is mathematically convenient, as the objective
function constructed from unconditional moments can be analytically computed.
Early work on this approach can be found in Hansen (1982), and present
variants are discussed in Bennett and Kallus (2023), among many others.
Particularly, when the unconditional moments are defined in a product space of
Euclidean and metric spaces, it is expected that the variables in Euclidean
space and the objects in metric space can be separated into different factors
such that the procedures of modeling and analyzing can be implemented
separately in the respective spaces.

In this paper, our goal is to establish the framework of nonlinear Fréchet
regression. The key is to elaborately design the related Fréchet objective
functions in the following form: they are expected weighted losses, in which
the weights are a type of nonlinear functions of covariates defined in an
Euclidean space, and the loss is the distance defined in a metric space.
Consequently, the Fréchet regression functions can be defined as the minimizer
of these Fréchet objective functions. Such a structural nonlinear framework is
essentially different from that of linear Fréchet regression. The new
proposals can be utilized to fit true nonlinear models, and their separable
structures can be easily analyzed and computed. Furthermore, the underlying
nonlinear models can be explicitly identified by the proposed Fréchet
regressions under a Hilbert space. Particularly, the generalized linear
Fréchet regression can return to the standard linear Fréchet regression
through a special choice of the weight functions. Theoretically, some
favorable properties, such as the estimation consistency and the convergence
rate, can be established systematically.

###  1.2 Problem setup

Let $\Omega$ be a metric space, equipped with a metric $d$. The Fréchet
regression function of an object response $Y\in\Omega$ on a covariate $X\in
R^{p}$ is defined as the minimizer of the conditional expectation of a loss
defined in the metric space (Petersen and Müller, 2019). Formally, the Fréchet
regression function is originally written as

$m_{\oplus}(x)=\operatorname*{arg\,min}_{\omega\in\Omega}E[d^{2}(Y,\omega)|X=x].$
(1.1)

However, it is difficult to solve the conditional expectation problem,
especially for the issue in a metric space. Just like the classical strategies
for the problems in Euclidean spaces, an approach to the conditional
expectation in a metric space is to reduce it to an unconditional expectation
with some weight factors (see, e.g., Hansen (1982) and Bennett and Kallus
(2023). Then, in order to easily implement the Fréchet regression, with
$\mu=E[X]$ and $\Sigma=Var[X]$, the existing literature (Petersen and Müller,
2019) utilizes the weight $s(X,x)=1+(X-\mu)^{T}\Sigma^{-1}(x-\mu)$ to
reconstruct the Fréchet regression function as the minimizer of the
unconditionally expected weighted loss; namely, the Fréchet regression
function is equivalently rewritten as

$m_{\oplus}(x)=\operatorname*{arg\,min}_{\omega\in\Omega}E[s(X,x)d^{2}(Y,\omega)].$
(1.2)

In the above, the weight function $s(X,x)$ is in a parameter framework and
exhibits linearity in $x$. When the Fréchet regression is redefined by the
above unconditional expectation, the resulting model is indeed linear in $x$
with the structure similar to the model (1.3) under a Hilbert space (see
Theorem 5 of Petersen and Müller (2019)). For emphasizing the linearity and
differentiating it from the nonlinearity studied in this paper, here we call
$m_{\oplus}(x)$ the linear Fréchet regression function, denoted by LFR for
short, and refer to $E[s(X,x)d^{2}(Y,\omega)]$ as the linear Fréchet objective
function. By the unconditional expectation expression, the LFR function can be
computed analytically for the case where the Euclidean variable $X$ and the
metric space object $Y$ are separated into different factors $s(X,x)$ and
$d^{2}(Y,\omega)$.

The LFR function defined in (1.2) always passes through the central point
$(\mu,\omega_{\oplus})$, since $s(\cdot,\mu)\equiv 1$ implies that
$m_{\oplus}(\mu)=\omega_{\oplus}$, where
$\omega_{\oplus}=\operatorname*{arg\,min}_{\omega\in\Omega}E[d^{2}(Y,\omega)]$
is the Fréchet mean of $Y$ (Fréchet, 1948). Particularly, consider the common
linear regression model in an Euclidean space defined by

$Y=\beta^{(0)}+\beta^{T}(X-\mu)+\varepsilon,$ (1.3)

where the response $Y\in R$, $\mu=E[X]\in R^{p}$, $\beta^{(0)}\in R$ and
$\beta\in R^{p}$. Under the above linear model defined in an Euclidean space,
the LFR function defined in (1.1) happens to be the minimizer of the
unconditional weighted expectation given in (1.2) with the linear weight
$s(X,x)$.

  1. However, it is intuitive that if the underlying model is actually nonlinear, the regression function should NOT always pass through the central point $(\mu,\omega_{\oplus})$. A typical example is that if $m_{\oplus}(x)=m(x)$ is a strictly convex function, then it always holds that $\omega_{\oplus}>m(\mu)$, which is against the rule of passing through the central point $(\mu,\omega_{\oplus})$. Furthermore, a linear Fréchet objective function uniquely identifies a linear regression function in a Hilbert space. Thus, if the underlying model is nonlinear, the weight in the objective should NOT always be of the linear framework.

The above discussions show that under the nonlinear case, the framework of
Fréchet regression should have essential differences from those in the linear
case. These observations motivate us to develop new methodology and theory for
constructing nonlinear Fréchet regression.

###  1.3 Article structure

The remainder of the paper is then organized in the following way. In Section
2, the nonlinear models in Euclidean space are first reexamined to motivate
the development of the nonlinear model in a metric space. Then, the definition
of nonlinear Fréchet regression is proposed, and particularly, a generalized
linear Fréchet regression is introduced. The related representation theorems
are investigated. For illustration, some important examples are discussed. In
Section 3, the estimation methods and the algorithms are presented, and the
asymptotic properties are established. The main simulation studies and the
main result of a human mortality data analysis are presented respectively in
Section 4 and Section 5. Section 6 concludes the paper and discusses the
possible future works. The Supplementary Materials include the following
materials: the regularity conditions and the proofs for lemmas and theorems,
some implementation details, and some additional simulation results and the
details of the real data analysis.

##  2 Nonlinear Fréchet regression modeling

Here, our goal is to establish the framework of nonlinear Fréchet regression
and investigate its basis properties.

###  2.1 A type of nonlinear Euclidean regressions

It is well-known that nonlinear regression is widely used in various fields. A
famous example of nonlinear regression is the Gompertiz model defined by

$Y=\theta_{1}\exp{[-\exp(\theta_{2}-\theta_{3}^{X})]}+\varepsilon.$ (2.1)

This is a sigmoidal growth model and commonly appears in some application
fields such as biology and economics. For more examples of nonlinear
regression models, particularly in biostatistical contexts, see Ratkowsky
(1983), Seber and Wild (2003) and Wei (1998), among others. Generally, a
nonlinear regression model has the following unified form:

$Y=f(X,\beta)+\varepsilon$ (2.2)

for a given function $f$ with covariate $X\in R^{p}$ and unknown parameter
$\beta\in R^{q}$. Particularly, when $p=q$, a special structure for the
regression function is $f(x,\beta)=g(\beta^{T}(x-\mu))$, where $g$ is a known
link function. In this case, the resulting model is a generalized linear
model, that is,

$Y=g(\beta^{T}(X-\mu))+\varepsilon.$ (2.3)

Under an exponential family distribution, the link function $g$ has a
canonical or natural choice. It is known that generalized linear regression is
one of the most extensive models in the field of statistics (see, e.g.,
McCulloch and Searle (2001) and Härdle et al. (2004)).

Denote $X=(X^{(1)},\cdots,X^{(p)})^{T}$ and
$\mu=(\mu^{(1)},\cdots,\mu^{(p)})^{T}=E[X]$, and suppose that $\Sigma=Var[X]$
exists and is positive. To motivate the methodological development, we
consider a type of nonlinear Euclidean regression of $Y\in R$ on $X\in R^{p}$
as

$m^{\cal
N}(x,\beta)=E[Y|X=x]=E[Y]+\sum_{j=1}^{p}\sigma^{(j)}f_{j}\left(x,\beta\right),$
(2.4)

where $\beta\in\mathscr{B}$ is an unknown parameter vector, the parameter
space $\mathscr{B}$ is a subset of $R^{q}$, $f_{j}$’s are known regression
functions, and $\sigma^{(j)}=E[Y(X^{(j)}-\mu^{(j)})]$. The above model (2.4)
is a special case of general regression (2.2), but it is a big class of
nonlinear regressions duo to the arbitrariness of the choices of $f_{j}$, and
it contains the famous Gompertiz model (2.1) as its special case. Also the
model can be derived from exponential families nonlinear models (see, e.g.,
Wei (1998)). The first condition for nonlinear modeling is the assumption of
$\sigma^{(j)}\neq 0$ for some $j$. This condition is easy to satisfy;
otherwise, it can be replaced by
$\sigma^{(j)}_{t}=E[Y(t(X^{(j)})-E[t(X^{(j)})])]\neq 0$ for some function $t$
and some index $j$. The condition $\sigma^{(j)}_{t}\neq 0$ is always satisfied
if there exist some (nonlinear) correlations between the response $Y\in R$ and
covariates $X^{(j)}\in R$. For the nonlinear Euclidean regression, the
nonlinear weight function is then defined by

$s^{\cal N}(X,x,\beta)=1+\sum_{j=1}^{p}(X^{(j)}-\mu^{(j)})f_{j}(x,\beta).$
(2.5)

The weight depends only on the Euclidean variables and parameters, and the
known functions $f_{j}$, which is free of $\sigma^{(j)}$ although the
underlying regression depends on $\sigma^{(j)}$. Particularly, under the
generalized linear case, the Euclidean regression function is designed as

$m^{\cal
N}(x,\beta)=E[Y|X=x]=E[Y]+\sum_{j=1}^{p}\sigma^{(j)}f_{j}\left(\beta^{T}(x-\mu)\right),$
(2.6)

and the corresponding generalized linear weight function is defined by

$s^{\cal
N}(X,x,\beta)=1+\sum_{j=1}^{p}(X^{(j)}-\mu^{(j)})f_{j}(\beta^{T}(x-\mu)).$
(2.7)

As aforementioned, the model (2.6) can be derived from the exponential
families generalized linear models.

Since $E[s^{\cal N}(X,x,\beta)]=1$ for all $x$ and $\beta$, the regression
function $m^{\cal N}(x,\beta)$ defined in (2.4) and (2.6) can be recast as

$\begin{split}m^{\cal N}(x,\beta)=E\left[s^{\cal
N}(X,x,\beta)Y\right],\end{split}$ (2.8)

which is the weighted expectation of $Y$ with the weight function $s^{\cal
N}(X,x,\beta)$ as in (2.5) or (2.7). Then, under Euclidean space, the
nonlinear regression function $m^{\cal N}(x,\beta)$ in (2.4) or (2.6) is
indeed the solution to the following optimization problem:

$m^{\cal N}(x,\beta)=\operatorname*{arg\,min}_{y\in R}E[s^{\cal
N}(X,x,\beta)d_{E}^{2}(Y,y)]$ (2.9)

for given $\beta$, where $d_{E}$ is the standard Euclidean metric.

###  2.2 Nonlinear Fréchet regression

The above alternative formulation (2.9) of the nonlinear Euclidean regression
function provides the key to defining the nonlinear Fréchet regression
function on an arbitrary metric space $(\Omega,d)$. By simply replacing the
Euclidean metric $d_{E}$ with a more general metric $d$ that is suitable for
responses in metric space $\Omega$, the nonlinear Fréchet regression function
is formally defined as

$m_{\oplus}^{\cal
N}(x,\beta)=\operatorname*{arg\,min}_{\omega\in\Omega}E[s^{\cal
N}(X,x,\beta)d^{2}(Y,\omega)]$ (2.10)

for given $\beta$ with the nonlinear weight function $s^{\cal N}(X,x,\beta)$
defined by (2.5) or (2.7). We call $m_{\oplus}^{\cal N}(x,\beta)$ the
nonlinear Fréchet regression function, denoted by NLFR for brevity, because
the weight function $s^{\cal N}(X,x,\beta)$ is in a nonlinear framework. A
further reason is that it will be proved later that the NLFR can uniquely
identify an underlying nonlinear regression function in a Hilbert space. The
pivotal features of the above definition are as follows:

  1. (1)

The loss (or the distance) $d^{2}(Y,\omega)$ is defined on the given metric
space;

  2. (2)

The nonlinearity of the regression function $m_{\oplus}^{\cal N}(x,\beta)$ is
completely determined by the nonlinearity of the weight function $s^{\cal
N}(X,x,\beta)$ that is only related to Euclidean variables;

  3. (3)

The Fréchet objective function is defined by an unconditional expectation, and
the Euclidean variable $X$ and the metric space object $Y$ are separated into
different factors $s^{\cal N}(X,x,\beta)$ and $d^{2}(Y,\omega)$.

In addition, it is necessary to achieve the adaptiveness to the linear case as
in (1.2). In other words, with some special choices of the functions
$f_{j}(\cdot)$ in the generalized linear case (2.7), the resulting model can
be reduced to an LFR. To this end, let $\beta=\Sigma^{-1}\sigma$ with
$\sigma=(\sigma^{(1)},\cdots,\sigma^{(p)})^{T}$, and rewrite $f_{j}$ as

$f_{j}\left(\sigma^{T}\Sigma^{-1}(x-\mu)\right)=f_{j}\Big{(}\sum_{i=1}^{p}\sigma^{(i)}\sum_{k=1}^{p}v^{(ik)}(x^{(k)}-\mu^{(k)})\Big{)}$

with $(v^{(ik)})_{i,k=1}^{p}=\Sigma^{-1}$. If each link function
$f_{j}(\cdot)$ is chosen as the following linear function:

$f_{j}(u)=\frac{u-\sum_{i\neq
j,i=1}^{p}\sigma^{(i)}\sum_{k=1}^{p}v^{(ik)}(x^{(k)}-\mu^{(k)})}{\sigma^{(j)}},$
(2.11)

then, $f_{j}\left(\sigma^{T}\Sigma^{-1}(x-\mu)\right)$ has the following
linear representation:

$f_{j}\left(\sigma^{T}\Sigma^{-1}(x-\mu)\right)=\sum_{k=1}^{p}v^{(jk)}(x^{(k)}-\mu^{(k)}).$

We then have the following finding:

  1. Based on the choice of $f_{j}(\cdot)$ in (2.11), the generalized linear weight $s^{\cal N}(X,x,\beta)$ given in (2.7) reduces to the linear weight $s(X,x)$. Consequently, the resulting model returns to the standard LFR defined in Petersen and Müller (2019).

Next, we study the basic properties of the NLFR and the representation
theorems. It is clear that the newly defined NLFR has the following favorable
properties:

  1. (1)

When the NLFR reduces to an LFR (i.e., all $f_{j},j=1,\cdots,p$, are given in
(2.11)), the regression curve always passes through the central point
$(\mu,\omega_{\oplus})$.

  2. (2)

When the NLFR function is essentially nonlinear (i.e., at least one of
$f_{j}(\cdot),j=1,\cdots,p$, is different from the function in (2.11)), the
regression curve does not passes through the central point
$(\mu,\omega_{\oplus})$ except for the case of $f_{j}(0)=0$ for
$j=1,\cdots,p$.

  3. (3)

The nonlinearity of the Fréchet regression is completely determined by the
nonlinearity of the Euclidean weight function $s^{\cal N}(X,x,\beta)$.

The above properties (1) and (2) ensure that the NLFR is adaptive to “the
center rule” according to whether the Fréchet regression is linear or not. The
property (3) shows how to check the nonlinearity.

Moreover, another issue is about how to explicitly identify the NLFR function
$m_{\oplus}^{\cal N}(x,\beta)$ when the random object $Y$ belongs to a Hilbert
space. The following representation theorem states the details.

######  Theorem 2.1.

Suppose that $\Omega$ is a Hilbert space equipped with an inner product
$\langle\cdot,\cdot\rangle$ and the corresponding norm
$\|\cdot\|^{2}_{\Omega}$. If $E\|Y\|^{2}_{\Omega}<\infty$, then, there exist
unique elements $\beta^{(0)}\in\Omega$ and $\sigma^{(j)}\in\Omega$ such that
$E\langle Y,\omega\rangle=\langle\beta^{(0)},\omega\rangle$ and
$E\langle(X^{(j)}-\mu^{(j)})Y,\omega\rangle=\langle\sigma^{(j)},\omega\rangle$
for all $\omega\in\Omega$. With these elements $\beta^{(0)}$ and
$\sigma^{(j)}$, the NLFR function $m_{\oplus}^{\cal N}(x,\beta)$ defined in
(2.10) has the following explicit representation:

$m^{\cal
N}_{\oplus}(x,\beta)=\beta^{(0)}+\sum_{j=1}^{p}\sigma^{(j)}f_{j}\left(x,\beta^{T}\right)\mbox{
or }m^{\cal
N}_{\oplus}(x,\beta)=\beta^{(0)}+\sum_{j=1}^{p}\sigma^{(j)}f_{j}\left(\beta^{T}(x-\mu)\right)$
(2.12)

for given $\beta$ according the weight function given by (2.5) or (2.7).

The representations in the theorem provide a convenient way for theoretical
analysis and practical application. Besides, similar to the Euclidean case as
in (2.8), the NLFR function can be expressed as a weighted expectation of $Y$.
The following corollary gives the details.

######  Corollary 2.1.

Under the conditions of Theorem 2.1, in a Hilbert space and for given $\beta$,
the NLFR function $m_{\oplus}^{\cal N}(x,\beta)$ defined in (2.10) has the
following representation:

$m^{\cal N}_{\oplus}(x,\beta)=E[s^{\cal N}(X,x,\beta)Y].$ (2.13)

The corollary provides a new perspective to the regression function: under a
Hilbert space, the NLFR function $m^{\cal N}_{\oplus}(x,\beta)$ can be
regarded as a weighted expectation of the response $Y$ with nonlinear weight
$s^{\cal N}(X,x,\beta)$. Thus, it could be expected that we can develop
various Fréchet regressions by weighted expectation with various weights. We
will continue to study this interesting issue in the future.

Next, we suggest the following examples that often appear in statistics to
illustrate the conditions of Theorem 2.1.

Example 1. Let $\Omega$ be the set of quantile functions $G^{-1}$ of the
probability distributions $G$ on $R$ such that $\int_{R}x^{2}dG(x)<\infty$,
equipped with the Wasserstein metric $d=d_{W}$. Under Wasserstein metric, the
Wasserstein distance is given by
$d^{2}_{W}(G_{1},G_{2})=\int_{0}^{1}(G^{-1}_{1}(t)-G^{-1}_{2}(t))^{2}dt$ for
two quantile functions $G_{1}^{-1}$ and $G_{2}^{-1}$ in $\Omega$. With these
conditions, $\Omega$ is a Hilbert space. In this case, the resulting Fréchet
regression function has the representations as in (2.12) and (2.13).

Example 2. Suppose that $\Omega$ is the set of correlation matrices of a fixed
dimension $r$, equipped with the Frobenius metric $d_{F}$. In this case,
$\Omega$ is a Hilbert space, and the resulting Fréchet regression function has
the representations as in (2.12) and (2.13) as well.

Example 3. It is known that the tangent vector field is often used to
formalize the notion of the Riemannian manifold. Let $\Omega$ be the tangent
space of a bounded Riemannian manifold of dimension $r$, and $d$ be defined by
the inner product in the tangent space. In the Hilbert space $\Omega$, the
resulting Fréchet regression function has the representations as in (2.12) and
(2.13) as well.

These examples and the corresponding representations for the NLFR functions
provide a basis for the simulation studies and real data analysis given in
Sections 4 and 5.

###  2.3 A separable representation for generalized linear form

In the previous subsection, the definition and the representation theorem for
the NLFR $m^{\cal N}_{\oplus}(x,\beta)$ are all based on the condition: the
parameter vector $\beta$ is given. Then, the estimation method proposed in the
next section depends on a profile method, which needs an iterative estimation
method to implement the estimation procedure. In this subsection, we introduce
a separable representation for generalized linear form to avoid the
theoretical and computational complexity.

For the generalized linear regression functions $f_{j}(\beta^{T}(x-\mu))$ in
(2.7), we consider a particular choice of $\beta$ as
$\beta=c_{h}\Sigma^{-1}\sigma_{h}$, where $h(y):y\in\Omega\mapsto R$ is a
given measurable function, and the constant $c_{h}\neq 0$ depends on $h$. The
choices of $h$ and $c_{h}$ will be given later. To demonstrate this structural
assumption, we provide the following motivating examples.

Example 4. Under linear Euclidean regression (1.3), the regression function is
$m(x)=E[Y]+\beta^{T}(x-\mu)$. It can be seen that the least square
representation of $\beta$ is $\beta=\Sigma^{-1}\sigma$ with
$\sigma=(\sigma^{(1)},\cdots,\sigma^{(p)})^{T}$ and
$\sigma^{(j)}=E[Y(X^{(j)}-\mu^{(j)})]$. Then the regression function can be
rewritten as $m(x)=E[Y]+\sigma^{T}\Sigma^{-1}(x-\mu)$. In this simple case,
$\beta=c_{h}\Sigma^{-1}\sigma_{h}$ with $c_{h}=1$ and $\sigma_{h}=\sigma$.

Example 5. Under an Euclidean space, if regression is expressed by standard
orthogonal basis $\\{\varphi_{1}(x),\cdots,\varphi_{p}(x)\\}$, i.e.,
$m(x)=E[Y]+\sum_{j=1}^{p}\theta^{(j)}\varphi_{j}(x)$, then it can be seen that
$\theta^{(j)}=\sigma^{(j)}$ and $\mu^{(j)}=0$ if $X^{(j)}$ is replaced by
$\varphi_{j}(X)$. In this case, we also have
$\beta=c_{h}\Sigma^{-1}\sigma_{h}$ with $c_{h}=1$ and $\sigma_{h}=\sigma$.

Example 6. Consider the following generalized linear Euclidean regression
function: $m^{\cal N}(x)=E[Y]+\sigma^{(j)}g(\beta^{T}(x-\mu))$ with
$\sigma^{(j)}=E[Y(X^{(j)}-\mu^{(j)})]\neq 0$ for a fixed $j$, where $g(\cdot)$
is a given nonlinear link function. The existing literature (see, e.g., Li and
Duan (1989) and Wang et al. (2012)) shows that if $E[X|\beta^{T}X]$ is a
linear function of $\beta^{T}X$, then $\beta=c_{h}\Sigma^{-1}\sigma_{h}$ for a
constant $c_{h}$ that depends on a function $h$.

Motivated by these observations, we introduce a separable generalized linear
weight function as

$s^{\cal
N}(X,x,c_{h})=1+\sum_{j=1}^{p}(X^{(j)}-\mu^{(j)})f_{j}\left(c_{h}\sigma_{h}^{T}\Sigma^{-1}(x-\mu)\right)$
(2.14)

for a nonzero constant $c_{h}$. Here the separability means that $\beta$ has
the separable expression $\beta=c_{n}\Sigma^{-1}\sigma_{h}$, and the
estimation procedures for parameter vector $\beta$ and the regression function
are separable (for the details see the next section). In this case, the
estimators of $\sigma_{h}^{T}$ and $\Sigma^{-1}$ can be constructed by the
data only from Euclidean space. This alternative method applies to situation
where there is some linear trend in the regression function (Wang, et al.,
2012), for example, $f_{j}\left(\beta^{T}(x-\mu)\right)=(\beta^{T}x+c)^{2}$
for a constant $c\neq 0$, which contains a linear term $2c\beta^{T}x$,
although the regression function is really nonlinear; otherwise it is possible
that $c_{h}=0$, for example,
$f_{j}\left(\beta^{T}(x-\mu)\right)=(\beta^{T}x)^{2}$, without any linear
information.

Then, for a given function $h$ (i.e. given $c_{h}$), the separable nonlinear
Fréchet regression function is formally defined by

$m_{\oplus}^{\cal
N}(x,c_{h})=\operatorname*{arg\,min}_{\omega\in\Omega}E[s^{\cal
N}(X,x,c_{h})d^{2}(Y,\omega)]$ (2.15)

with the separable weight function $s^{\cal N}(X,x,c_{h})$ given in (2.14). We
denote the separable nonlinear Fréchet regression by SNLFR for short.

##  3 Estimations and asymptotic properties

###  3.1 Estimation for NLFR and its asymptotic properties

Here, we first introduce the estimation method for the NLFR function
$m_{\oplus}^{\cal N}(x,\beta)$ defined in (2.10). Let
$(X_{i},Y_{i}),i=1,\cdots,n$, be independent and identically (i.i.d.)
observations of $(X,Y)$, and $\widehat{\mu}$ and $\widehat{\Sigma}$ be the
empirical estimators of $\mu$ and $\Sigma$, respectively. For example, the
estimators can be chosen as $\widehat{\mu}=\frac{1}{n}\sum_{i=1}^{n}X_{i}$ and
$\widehat{\Sigma}=\frac{1}{n}\sum_{i=1}^{n}(X_{i}-\widehat{\mu})(X_{i}-\widehat{\mu})^{T}$.
Similar to the notations used in the previous section, the estimated weight
function is denoted by $\widehat{s}^{\cal
N}(X,x,\beta)=1+\sum_{j=1}^{p}(X^{(j)}-\widehat{\mu}^{(j)})f_{j}(x,\beta)$ or
$\widehat{s}^{\cal
N}(X,x,\beta)=1+\sum_{j=1}^{p}(X^{(j)}-\widehat{\mu}^{(j)})f_{j}(\beta^{T}(x-\mu))$.
Then, given $\beta$, the estimator of the NLFR function can be obtained by

$\widehat{m}_{\oplus}^{\cal
N}(x,\beta)=\operatorname*{arg\,min}_{\omega\in\Omega}\frac{1}{n}\sum_{i=1}^{n}\widehat{s}^{\cal
N}(X_{i},x,\beta)d^{2}(Y_{i},\omega).$ (3.1)

Similar to the representation in Theorem 2.1, in a Hilbert space and given
$\beta$, the estimator of NLFR in (3.1) can be explicitly expressed as

$\widehat{m}^{\cal
N}_{\oplus}(x,\beta)=\widehat{\beta}^{(0)}+\sum_{j=1}^{p}\widehat{\sigma}^{(j)}f_{j}(x,\beta)\mbox{
or }\widehat{m}^{\cal
N}_{\oplus}(x,\beta)=\widehat{\beta}^{(0)}+\sum_{j=1}^{p}\widehat{\sigma}^{(j)}f_{j}(\beta^{T}(x-\widehat{\mu})),$
(3.2)

where $\widehat{\beta}^{(0)}=\frac{1}{n}\sum_{i=1}^{n}Y_{i}$ and
$\widehat{\sigma}^{(j)}=\frac{1}{n}\sum_{i=1}^{n}Y_{i}(X_{i}^{(j)}-\widehat{\mu}^{(j)})$.
The representation ensures that the estimator $\widehat{m}^{\cal
N}_{\oplus}(x,\beta)$ is identified as a nonlinear model as well, and its
structure matches up to the theoretical version given in Theorem 2.1. Also,
the estimator can be constructed by a weighted sum of $Y_{i}$ with the
estimated nonlinear weight $\widehat{s}^{\cal N}(X_{i},x,\beta)$, that is,

$\widehat{m}^{\cal
N}_{\oplus}(x,\beta)=\frac{1}{n}\sum_{i=1}^{n}\widehat{s}^{\cal
N}(X_{i},x,\beta)Y_{i},$ (3.3)

which also matches up to the theoretical version given in Corollary 2.1. The
above representations are useful for theoretical analysis and numerical
calculation.

In the estimation procedure for the NLFR estimator $\widehat{m}^{\cal
N}_{\oplus}(x,\beta)$ in (3.1), the parameter vector $\beta$ cannot be
estimated separately from the estimation procedure. Motivated by the profile
likelihood in semiparametric model (see, e.g., Severini and Wong (1992)), a
profile estimation procedure is then designed as follows:

  1. (1)

Given $\beta$, an initial estimator of $m_{\oplus}^{\cal N}(x,\beta)$ is given
as the solution to the optimization problem in (3.1).

  2. (2)

Based on the obtained estimator $\widehat{m}_{\oplus}^{\cal N}(x,\beta)$, the
estimator of $\beta$ can be attained as the solution to the following
optimization problem:

$\widehat{\beta}_{\widehat{m}}=\operatorname*{arg\,min}_{\beta\in\mathscr{B}}\frac{1}{n}\sum_{i=1}^{n}d^{2}(Y_{i},\widehat{m}_{\oplus}^{\cal
N}(X_{i},\beta)).$

  3. (3)

Final estimator of the regression function is $\widehat{m}_{\oplus}^{\cal
N}(x,\widehat{\beta}_{\widehat{m}})$.

Particularly, under the situation where $Y$ is defined in a Hilbert space, the
profile algorithm can be designed as follows.

  1. Step

1\. Construct the estimators $\widehat{\mu}$, $\widehat{\beta}^{(0)}$ and
$\widehat{\sigma}^{(j)}$ by the aforementioned methods.

  2. Step

2\. Using the representation $\widehat{m}_{\oplus}^{\cal N}(X_{i},\beta)$
given in (3.2) or (3.3), construct the estimator of parameter $\beta$ by

$\widehat{\beta}=\operatorname*{arg\,min}_{\beta\in\mathscr{B}}\frac{1}{n}\sum_{i=1}^{n}d^{2}(Y_{i},\widehat{m}_{\oplus}^{\cal
N}(X_{i},\beta)).$ (3.4)

  3. Step

3\. With the estimator $\widehat{\beta}$ obtained in Step 2, compute the
estimator of the regression function by the following explicit representation
$\widehat{m}^{\cal
N}_{\oplus}(x,\widehat{\beta})=\widehat{\beta}^{(0)}+\sum_{j=1}^{p}\widehat{\sigma}^{(j)}f_{j}(x,\widehat{\beta})$
or $\widehat{m}^{\cal
N}_{\oplus}(x,\widehat{\beta})=\widehat{\beta}^{(0)}+\sum_{j=1}^{p}\widehat{\sigma}^{(j)}f_{j}(\widehat{\beta}^{T}(x-\widehat{\mu}))$.

Note that the solution of parameter estimator $\widehat{\beta}$ in (3.4) does
not have an exact explicit representation. Then, the numerical iterative
solution is required in Step 2.

For the estimation consistency, we only present the conclusions for the case
where the regression functions are of the form of $f_{j}(x,\beta)$. For the
case of $f_{j}(\beta^{T}(x-\mu))$, the asymptotic theory is similar; the
detail is omitted here. The regularity conditions are given in the
Supplementary Materials. These regularity conditions are mainly related to,
for example, the existence and uniqueness of the Fréchet regression functions
and their estimators, and that the expected loss $E[s^{\cal
N}(X,x,\beta)d^{2}(Y,\omega)]$ is significantly minimized at $m_{\oplus}^{\cal
N}(x,\beta)$. Moreover, the existing literature (Petersen and Müller, 2019)
has proven that the models in Examples 1-3 satisfy these regularity
conditions.

######  Theorem 3.1.

Under the regularity conditions C1-C4 given in the Supplementary Materials,

(i) if all the regression functions $f_{j}(x,\beta)$ are uniformly continuous
for $x\in R^{p}$ and $\beta\in\mathscr{B}$, then it holds that

$d_{E}(\widehat{\beta},\beta_{*})=o_{p}(1)\mbox{ and
}d(\widehat{m}_{\oplus}^{\cal N}(x,\widehat{\beta}),m_{\oplus}^{\cal
N}(x,\beta_{*}))=o_{p}(1)$ (3.5)

for any fixed $x\in R^{p}$, where
$\beta_{*}=\operatorname*{arg\,min}_{\beta\in\mathscr{B}}E[d^{2}(Y,m_{\oplus}^{\cal
N}(X,\beta))]$.

(ii) if without the uniform continuity, suppose that all the regression
functions $f_{j}(x,\beta)$ are continuous for $x\in R^{p}$ and
$\beta\in\mathscr{B}$, and the boundedness of $\|x\|_{E}\leq B$,
$\mbox{diam}(\mathscr{B})\leq B$ and $\|X_{i}\|_{E}\leq B$ holds for a
constant $B>0$ and uniformly for all $i=1,\cdots,n$, then, the above
consistency in (3.5) is satisfied as well.

Furthermore, the following theorem presents the convergence rate.

######  Theorem 3.2.

(i) If all the regression functions $f_{j}(x,\beta)$ are uniformly continuous
for $x\in R^{p}$ and $\beta\in\mathscr{B}$, under the regularity conditions
C1-C6 in the Supplementary Materials, for any fixed $x\in R^{p}$, it holds
that

$d(\widehat{m}_{\oplus}^{\cal N}(x,\widehat{\beta}),m_{\oplus}^{\cal
N}(x,\beta_{*}))=O_{p}(n^{-1/(2(\gamma-1))}),$

where the constant $\gamma>1$ is defined in the regularity condition C6.

(ii) If all the regression functions $f_{j}(x,\beta)$ are uniformly continuous
for $x\in R^{p}$ and $\beta\in\mathscr{B}$, under the regularity conditions
C1-C8, it holds that

$\sup_{x\in R^{p}}d(\widehat{m}_{\oplus}^{\cal
N}(x,\widehat{\beta}),m_{\oplus}^{\cal
N}(x,\beta_{*}))=O_{p}(n^{-1/(2(\alpha^{\prime}-1))})$

for some constant $\alpha^{\prime}>\alpha$, where $\alpha>1$ is defined in the
regularity condition C8.

If without the condition of uniform continuity, we also can get the same
convergence rate for the case of $\|x\|_{E}\leq B$ and
$\mbox{diam}(\mathscr{B})\leq B$ (the details are omitted here). The
parameters $\gamma$ and $\alpha$ defined in the regularity conditions depend
on the model conditions. In some common cases, such as Examples 1-3, the
parameters can be definitely determined as $\gamma=\alpha=2$ (Petersen and
Müller, 2019). This implies the convergence rates of parameter estimation. It
is not surprising to attain such convergence rates, because the models are
actually in an underlying parametric framework.

###  3.2 Estimation for SNLFR and its asymptotic properties

Let
$\widehat{\sigma}_{h}=\frac{1}{n}\sum_{i=1}^{n}h(Y_{i})(X_{i}-\widehat{\mu})$.
Given $h$ and $c_{h}$, the estimator of the weight function can be expressed
as $\widehat{s}^{\cal
N}(X,x,c_{h})=1+\sum_{j=1}^{p}(X^{(j)}-\widehat{\mu}_{t}^{(j)})f_{j}\left(c_{h}\widehat{\sigma}_{h}^{T}\widehat{\Sigma}^{-1}(x-\widehat{\mu})\right)$.
Then, for given $h$ and $c_{h}$, we obtain the estimator of the SNLFR function
$m_{\oplus}^{\cal N}(x,c_{h})$ as

$\widehat{m}_{\oplus}^{\cal
N}(x,c_{h})=\operatorname*{arg\,min}_{\omega\in\Omega}\frac{1}{n}\sum_{i=1}^{n}\widehat{s}^{\cal
N}(X_{i},x,c_{h})d^{2}(Y_{i},\omega).$ (3.6)

The main difference from the estimation procedure of the NLFR proposed in the
previous subsection is that here we does not need the profile method to
estimate the SNLFR function $\widehat{m}_{\oplus}^{\cal N}(x,c_{h})$, because
$c_{h}$ is an external parameter and independent of the model, and it only
depends on the function $h$. Besides, similar to Theorem 2.1, we have the
following representation: under a Hilbert space with the conditions of Theorem
2.1, it holds that

$\widehat{m}^{\cal
N}_{\oplus}(x,c_{h})=\widehat{\beta}^{(0)}+\sum_{j=1}^{p}\widehat{\sigma}^{(j)}f_{j}\left(c_{h}\widehat{\sigma}^{T}_{h}\widehat{\Sigma}^{-1}(x-\widehat{\mu})\right)=\frac{1}{n}\sum_{i=1}^{n}\widehat{s}^{\cal
N}(X_{i},x,c_{h})Y_{i}.$ (3.7)

The above representations match up to the theoretical versions given in
Theorem 2.1 and Corollary 2.1.

Theoretically, a special condition given in the following theorems is the
linearity condition on $X$: $E[X|\beta^{T}X]$ is a linear function of
$\beta^{T}X$. This linearity condition is satisfied when $X$ has an elliptical
distribution (see, e.g., Li (1991), Cook (1998) and Cook and Ni (2005)). With
these conditions, the following theorem presents the estimation consistency:

######  Theorem 3.3.

Suppose that $E[X|\beta^{T}X]$ is a linear function of $\beta^{T}X$, and all
the link functions $f_{j}(u)$ are continuous for $u\in R$.

(i) Under the regularity conditions C1’-C3’ in the Supplementary Materials,
for given $h$ and $c_{h}$, and any fixed $x\in R^{p}$, it holds that
$d(\widehat{m}_{\oplus}^{\cal N}(x,c_{h}),m_{\oplus}^{\cal
N}(x,c_{h}))=o_{p}(1).$

(ii) Under the regularity conditions C1’, C2’ and C4’ in the Supplementary
Materials, for given $h$ and $c_{h}$, it holds that

$\sup_{\|x\|_{E}\leq B}d(\widehat{m}_{\oplus}^{\cal
N}(x,c_{h}),m_{\oplus}^{\cal N}(x,c_{h}))=o_{p}(1)$

for some constant $B>0$. In addition to the above conditions, if all the link
functions $f_{j}(u)$ are uniformly continuous for $u\in R$, then,

$\sup_{x\in R^{p}}d(\widehat{m}_{\oplus}^{\cal N}(x,c_{h}),m_{\oplus}^{\cal
N}(x,c_{h}))=o_{p}(1).$

Furthermore, the following theorem presents the convergence rate.

######  Theorem 3.4.

Suppose that $E[X|\beta^{T}X]$ is a linear function of $\beta^{T}X$, and all
the link functions $f_{j}(u)$ are continuous for $u\in R$.

(i) Under the regularity conditions C1’-C6’ in the Supplementary Materials,
for given $h$ and $c_{h}$, and any fixed $x\in R^{p}$, it holds that

$d(\widehat{m}_{\oplus}^{\cal N}(x,c_{h}),m_{\oplus}^{\cal
N}(x,c_{h}))=O_{p}(n^{-1/(2(\gamma-1))}),$

where the constant $\gamma>1$ is defined in the regularity condition C6’.

(ii) Furthermore, under the regularity conditions C1’-C8’, it holds that

$\sup_{\|x\|_{E}\leq B}d(\widehat{m}_{\oplus}^{\cal
N}(x,c_{h}),m_{\oplus}^{\cal
N}(x,c_{h}))=O_{p}(n^{-1/(2(\alpha^{\prime}-1))})$

for some constants $\alpha^{\prime}>\alpha$ and $B>0$, where $\alpha>1$ is
defined in the regularity condition C8’. In addition to the above conditions,
if all the link functions $f_{j}(u)$ are uniformly continuous for $u\in R$,
then,

$\sup_{x\in R^{p}}d(\widehat{m}_{\oplus}^{\cal N}(x,c_{h}),m_{\oplus}^{\cal
N}(x,c_{h}))=O_{p}(n^{-1/(2(\alpha^{\prime}-1))}).$

Before ending this subsection, we discuss how to choose the function
$h(Y):Y\in\Omega\mapsto R$ and how to identify the constant $c_{h}$. It can be
seen from the procedure of estimating the SNLFR function that we need the
function $h(Y)$ only for independently identifying the parameter vector
$\beta$ in the link functions $f_{j}(\beta^{T}(x-\mu))$ with the special
choice of $\beta$ as $\beta=c_{h}\Sigma^{-1}\sigma_{h}$. That is to say that
after identifying $\beta=c_{h}\Sigma^{-1}\sigma_{h}$, the estimation procedure
for $\beta$ is mostly independent of the metric space data, making the
estimation procedure easy to implement.

As shown in Example 6, the constant $c_{h}$ depends only on the chosen
function $h$, implying that we should first choose $h$ and then determine the
constant $c_{h}$. As stated in the existing literature (Wang et al., 2012), a
common choice of $h(\cdot)$ is $h(y)=G(y)$, the distribution function of $Y$,
if it is well-defined in a metric space. When the distribution $G(y)$ is
unknown in practice, we may use the empirical distribution function
$\widehat{G}(y)$ of $Y$. If without the distribution in a metric space, $h(Y)$
can be chosen as a one-to-one correspondence between the metric space $\Omega$
and a subset of $R$. Then, for a given $h$, an optimal choice of the constant
$c_{h}$ is defined by

$c^{*}_{h}=\operatorname*{arg\,min}_{c_{h}\in
R}\frac{1}{n}\sum_{i=1}^{n}d^{2}(Y_{i},\widehat{m}_{\oplus}^{\cal
N}(X_{i},c_{h})).$ (3.8)

By the representation Theorem 2.1, the SNLFR function
$\widehat{m}_{\oplus}^{\cal N}(x,c_{h})$ is uniquely identified when $h$ is
given. It guarantees that the above optimal choice of $c^{*}_{h}$ is unique as
well under a Hilbert space. The ideal optimality (3.8) leads to the following
empirical version:

$c^{e}_{h}=\operatorname*{arg\,min}_{c_{h}\in\mathscr{C}}\frac{1}{n}\sum_{i=1}^{n}d^{2}(Y_{i},\widehat{m}_{\oplus}^{\cal
N}(X_{i},c_{h}))$ (3.9)

for a given finite subset $\mathscr{C}\subset R$.

More generally, the optimal strategy of simultaneous choices for both
$h(\cdot)$ and $c_{h}$ will be given in the Supplementary Materials.

##  4 Numerical analyses

The main simulation studies are reported in this section, and some additional
numerical results will be given in the Supplementary Materials. In the
simulation experiments, two types of responses are considered, the first one
is the one-dimensional probability distributions and the second one is
symmetric positive definitive matrices. For a comprehensive comparison, we
consider the two competitors, the linear Fréchet regression (LFR), also known
as global Fréchet regression, and the local Fréchet regression (LoFR)
(Petersen and Müller, 2019). The computations for the LFR and the LoFR can be
carried out by R package frechet (Chen et al., 2020). Following Petersen and
Müller (2019), a grid of bandwidths $\\{0.05,0.1,\ldots,0.3\\}$ together with
the Gaussian kernel are used for constructing LoFR, and the optimal bandwidth
is selected by 10-fold cross-validation.

For given the independently generated testing data
$\\{\widetilde{X}_{i},\widetilde{Y}_{i}\\}_{i=1}^{\widetilde{n}}$ and a
specific metric $d(\cdot,\cdot)$, we adopt two types of mean squared errors
$\mathrm{MSE}_{Y}=\frac{1}{\widetilde{n}}\sum_{i=1}^{\widetilde{n}}d^{2}(\widehat{m}_{\oplus}(\widetilde{X}_{i}),\widetilde{Y}_{i})$
and
$\mathrm{MSE}_{m}=\frac{1}{\widetilde{n}}\sum_{i=1}^{\widetilde{n}}d^{2}(\widehat{m}_{\oplus}(\widetilde{X}_{i}),m_{\oplus}(\widetilde{X}_{i}))$
to evaluate the performance of each method, where
$\widehat{m}_{\oplus}(\cdot)$ is an estimated Fréchet regression function
based respectively on the SNLFR, the NLFR, the LFR and the LoFR. In the
following simulations, all the results are obtained by computing the averaged
values of $\mathrm{MSE}_{Y}$ and $\mathrm{MSE}_{m}$ with 100 replications. In
addition, the quality of the parameter estimator of $\beta$ in the NLFR is
measured by the averaged squared error across replications defined by
$\mathrm{ASE}_{\beta}=\frac{1}{100}\sum_{r=1}^{100}\|\widehat{\beta}^{[r]}-\beta\|_{2}^{2},$
where $\widehat{\beta}^{[r]}$ is the parameter estimator in the $r$th
replication.

###  4.1 Fréchet regression for probobility distributions

Let $\Omega$ be the metric space of probability distributions on $R$ with
finite second-order moments, equipped with the Wasserstein metric distance
$d_{W}^{2}(\cdot,\cdot)$ (see Example 1 for a detailed introduction). In this
subsection, the responses $Y$ represent the distribution function with the
quantile function $Q(Y)$. Following Petersen and Müller (2019), for the sake
of notational simplicity, we also denote the quantile function corresponding
to $Y$ as $Y$. To achieve a computationally manageable implementation of the
Wasserstein metric distance, the discrete approximation of
$d_{W}^{2}(G_{1},G_{2})$ is computed as
$m^{-1}\sum_{i=1}^{m}(G_{1}^{-1}(t_{i})-G_{2}^{-1}(t_{i}))^{2}$ based on an
equally spaced grid $\\{t_{1},\ldots,t_{m}\\}$ on the interval $[0,1]$ with
$m=20$ (Tucker et al., 2023).

The regression function is given by
$m_{\oplus}(x)=E(Y(\cdot)|X=x)=U_{0}+\alpha^{\mathrm{T}}g(x)+(V_{0}+\alpha^{\mathrm{T}}g(x))(\Phi^{-1}(\cdot)+1),$
where $\Phi(\cdot)$ is the standard normal distribution function,
$\alpha=(\alpha_{1},\ldots,\alpha_{p})^{\mathrm{T}}\in R^{p}$, $U_{0},V_{0}\in
R$ and $g(x)=(g_{1}(x,\beta),\ldots,g_{p}(x,\beta))^{\mathrm{T}}$ with the
known component functions $g_{1},\ldots,g_{p}$. For the quantile function, the
random response $Y$ is generated by adding noise as follows:
$Y(\cdot)=U+V(\Phi^{-1}(\cdot)+1)$ with $U\mid X\sim
N(U_{0}+\alpha^{\mathrm{T}}g(X),v_{1})$ and $V\mid
X\sim\operatorname{Gamma}((V_{0}+\alpha^{\mathrm{T}}g(X))^{2}/v_{2},v_{2}/(V_{0}+\alpha^{\mathrm{T}}g(X)))$.

We next provide some implementation details in the SNLFR and the NLFR, such as
the determination of the function $h(Y)$ and the link function
$f(x)=(f_{1}(x,\beta),\ldots,f_{p}(x,\beta))^{T}$. To better reflect the
information of $h(Y)$, we choose
$h(Y)=U-E(U)+(V-E(V))(\int_{0}^{1}\Phi^{-1}(t)dt+1)=U+V-E(U+V)$ for the SNLFR.
Note that based on an equally spaced grid $\\{t_{1},\ldots,t_{m}\\}$ on
$[0,1]$, a discrete representation of the response $Y$ is
$(Y(t_{1}),\ldots,Y(t_{m}))^{\mathrm{T}}$. Thus, the function $h(Y)$ can be
calculated by $Y^{(m)}-E[Y^{(m)}]$ with
$Y^{(m)}=m^{-1}\sum_{i=1}^{m}Y(t_{i})$. Here, the expectation $E[Y^{(m)}]$ can
be approximately computed by its empirical estimator given a random sample.
The optimal choice of $c_{h}$ is further derived via (3.9) in Section 3.2,
where the optim function in the R program can be used with the option “Brent”.
Moreover, we discuss how to choose the link function $f(x)$. Let
$\sigma=(\sigma^{(1)},\ldots,$ $\sigma^{(p)})^{\mathrm{T}}$. By the
representation Theorem 2.1, it is readily shown that
$\beta^{(0)}=U_{0}+\alpha^{\mathrm{T}}E[g(X)]+(V_{0}+\alpha^{\mathrm{T}}E[g(X)])(\Phi^{-1}(\cdot)+1)$
and $\sigma=\mathrm{Cov}(X,g(X))(2\alpha+\alpha\Phi^{-1}(\cdot))$, where
$\mathrm{Cov}(X,g(X))=E[(X-\mu)(g(X)-E[g(X)])^{\mathrm{T}}]$. According to
$m_{\oplus}(x)=\beta^{(0)}+\sigma^{\mathrm{T}}f(x)$, the link function $f(x)$
can be determined by $\mathrm{Cov}(g(X),X)f(x)=g(x)-E[g(X)].$

In the following, we compare different Fréchet regression procedures under
various models and the settings of sample sizes and dimensions.

Model 1.1. (Linear Fréchet Regression) Consider
$(X^{\prime(1)},\ldots,X^{\prime(p)})\sim N_{p}(0,\Sigma)$ with $\Sigma$
$=(0.5^{|i-j|})_{p\times p}$, and then set $X^{(j)}=2\Phi(X^{\prime(j)})-1$
for $j=1,\ldots,p$. Clearly, $X^{(j)}\sim\mathcal{U}(-1,1)$ for each $j$. Let
$\alpha=(1,0,\ldots,0)^{\mathrm{T}}$,
$g(x)=(\beta^{\mathrm{T}}(x-\mu),0,\ldots,0)^{\mathrm{T}}$ and
$f(x)=(f_{1}(x,\beta),0,\ldots,0)^{\mathrm{T}}$. Consider $p\in\\{2,5\\}$ and
$n\in\\{100,200,500\\}$. For $p=2$, we choose $\beta=(1,-0.5)^{\mathrm{T}}$
and $V_{0}=2$. For $p=5$, we choose $\beta=(1,-0.5,2,1.5,-1)^{\mathrm{T}}$ and
$V_{0}=6.5$, the additional parameters are set as $U_{0}=0$, $v_{1}=1$,
$v_{2}=0.5$ and $\widetilde{n}=500$.

Model 1.2. (Generalized linear Fréchet Regression) Let
$(X^{(1)},\ldots,X^{(p)})\sim N_{p}(0,\Sigma)$ with
$\Sigma=(0.5^{|i-j|})_{p\times p}$ and
$g(x)=((\beta^{\mathrm{T}}(x-\mu)+1)^{2},0,\ldots,0)^{\mathrm{T}}$. The
additional parameter $V_{0}$ is set as 0.5. The other settings are the same as
in Model 1.1.

Model 1.3. (Nonlinear Fréchet Regression) Let
$g(x)=(g_{1}(x,\beta),0,\ldots,0)^{\mathrm{T}}$ with
$g_{1}(x,\beta)=\sum_{j=1}^{p_{1}}(\beta^{(j)}(x^{(j)}-\mu^{(j)})+1)^{2}+\sum_{j=p_{1}+1}^{p}\exp{(\beta^{(j)}(x^{(j)}-\mu^{(j)}))}$,
where $\beta^{(j)}$ is the $j$th element in $\beta$. Consider the case of
$p_{1}=1$ for $p=2$ and $p_{1}=3$ for $p=5$. The other settings are the same
as in Model 1.2.

For Models 1.1–1.3 with $p=2$, the averaged values of $\mathrm{MSE}_{Y}$ and
$\mathrm{MSE}_{m}$ with standard errors in parentheses over 100 repetitions
are reported on Table 1. It can be seen from Table 1 that the SNLFR, the NLFR
and the LFR yield similar $\mathrm{MSE}_{Y}$ and $\mathrm{MSE}_{m}$ values for
Model 1.1, and moreover, these $\mathrm{MSE}_{Y}$ and $\mathrm{MSE}_{m}$
values are smaller than those of the LoFR, especially when the sample size is
small. For Model 1.2 with a generalized linear structure, the NLFR has the
best performance with the SNLFR coming in a close second. The performance of
LFR is the worst for Model 1.2 since the LFR is only suitable for linear
Fŕechet regression. For the nonlinear model Model 1.3, the NLFR still
maintains the best performance and significantly outperforms the LFR and the
LoFR under both the $\mathrm{MSE}_{Y}$ and $\mathrm{MSE}_{m}$ measures.

Table 1:  The averaged $\mathrm{MSE}_{Y}$ and $\mathrm{MSE}_{m}$ of various
methods and the associated standard errors (in parenthesis) for Models 1.1–1.3
with $p=2$.

|  |  | $\mathrm{MSE}_{Y}$ | $\mathrm{MSE}_{m}$
---|---|---|---|---
|  |  | SNLFR | NLFR | LFR | LoFR | SNLFR | NLFR | LFR | LoFR
| Model 1.1 | $n=100$ | 1.914 | 1.914 | 1.916 | 2.063 | 0.054 | 0.054 | 0.055 | 0.206
|  | (0.011) | (0.012) | (0.011) | (0.015) | (0.004) | (0.004) | (0.004) | (0.010)
| $n=200$ | 1.912 | 1.911 | 1.912 | 1.981 | 0.025 | 0.025 | 0.025 | 0.096
|  | (0.011) | (0.010) | (0.010) | (0.011) | (0.002) | (0.002) | (0.002) | (0.003)
| $n=500$ | 1.896 | 1.896 | 1.896 | 1.923 | 0.011 | 0.011 | 0.011 | 0.039
|  | (0.011) | (0.011) | (0.011) | (0.011) | (0.001) | (0.001) | (0.001) | (0.001)
| Model 1.2 | $n=100$ | 2.127 | 1.979 | 7.424 | 5.143 | 0.234 | 0.085 | 5.516 | 3.244
|  | (0.023) | (0.015) | (0.099) | (0.586) | (0.022) | (0.005) | (0.092) | (0.583)
| $n=200$ | 1.982 | 1.928 | 7.101 | 3.144 | 0.091 | 0.039 | 5.248 | 1.256
|  | (0.016) | (0.010) | (0.088) | (0.133) | (0.010) | (0.002) | (0.088) | (0.132)
| $n=500$ | 1.918 | 1.898 | 7.148 | 2.509 | 0.036 | 0.015 | 5.262 | 0.682
|  | (0.012) | (0.011) | (0.092) | (0.146) | (0.003) | (0.001) | (0.084) | (0.147)
| Model 1.3 | $n=100$ | – | 1.984 | 13.051 | 7.087 | – | 0.102 | 11.222 | 5.227
|  | – | (0.012) | (0.172) | (1.753) | – | (0.006) | (0.172) | (1.763)
| $n=200$ | – | 1.917 | 12.529 | 4.399 | – | 0.043 | 10.549 | 2.513
|  | – | (0.011) | (0.171) | (0.943) | – | (0.002) | (0.161) | (0.936)
| $n=500$ | – | 1.882 | 11.887 | 2.237 | – | 0.015 | 10.012 | 0.367
|  | – | (0.010) | (0.165) | (0.038) | – | (0.001) | (0.163) | (0.037)

The simulation results in the case of $p=5$ are summarized in Table 2. It can
be observed that although the $\mathrm{MSE}_{Y}$ and $\mathrm{MSE}_{m}$ values
of NLFR are larger relative to the case of $p=2$, they remain the smallest for
Model 1.1–1.3 among all competitors. It is noteworthy that LFR shows quite
poor performance compared to NLFR for Models 1.2–1.3, indicating that it is
not suitable to fit generalized linear and nonlinear structures. What needs to
be explained is that SNLFR produces larger $\mathrm{MSE}_{Y}$ and
$\mathrm{MSE}_{m}$ values than those of the NLFR for Model 1.2. It is because
the function $h(Y)$ causes a loss of the information of $Y$ on covariates,
especially when the sample size is small and the dimensionality of covariates
is high. In addition, as the increase of sample size, the $\mathrm{MSE}_{Y}$
and $\mathrm{MSE}_{m}$ values of the SNLFR remarkably decrease, which shows
that the SNLFR is also effective in dealing with the generalized linear
structures. Considering the curse of dimensionality for the LoFR, the package
frechet only handles cases where the dimension of $X$ is not larger than 2.
For the case of $p=5$, the simulation results of the LoFR are not shown in
Table 2.

Table 2:  The averaged $\mathrm{MSE}_{Y}$ and $\mathrm{MSE}_{m}$ of various
methods and the associated standard errors (in parenthesis) for Models 1.1–1.3
with $p=5$.

|  |  | $\mathrm{MSE}_{Y}$ | $\mathrm{MSE}_{m}$
---|---|---|---|---
|  |  | SNLFR | NLFR | LFR | LoFR | SNLFR | NLFR | LFR | LoFR
| Model 1.1 | $n=100$ | 1.996 | 1.986 | 1.986 | – | 0.125 | 0.119 | 0.119 | –
|  | (0.013) | (0.012) | (0.013) | – | (0.007) | (0.006) | (0.006) | –
| $n=200$ | 1.921 | 1.917 | 1.917 | – | 0.061 | 0.056 | 0.056 | –
|  | (0.010) | (0.010) | (0.010) | – | (0.003) | (0.002) | (0.002) | –
| $n=500$ | 1.894 | 1.892 | 1.892 | – | 0.024 | 0.023 | 0.023 | –
|  | (0.010) | (0.010) | (0.010) | – | (0.001) | (0.001) | (0.001) | –
| Model 1.2 | $n=100$ | 288.349 | 2.248 | 740.977 | – | 286.633 | 0.368 | 739.117 | –
|  | (26.216) | (0.086) | (14.123) | – | (26.212) | (0.083) | (14.182) | –
| $n=200$ | 141.654 | 2.047 | 697.442 | – | 139.878 | 0.171 | 695.946 | –
|  | (11.447) | (0.045) | (10.926) | – | (11.471) | (0.044) | (10.883) | –
| $n=500$ | 56.771 | 1.917 | 694.413 | – | 54.915 | 0.049 | 692.475 | –
|  | (5.012) | (0.012) | (12.014) | – | (5.005) | (0.003) | (12.002) | –
| Model 1.3 | $n=100$ | – | 2.913 | 346.129 | – | – | 1.023 | 344.159 | –
|  | – | (0.148) | (29.707) | – | – | (0.149) | (29.661) | –
| $n=200$ | – | 2.594 | 306.084 | – | – | 0.680 | 304.674 | –
|  | – | (0.103) | (5.95) | – | – | (0.099) | (5.938) | –
| $n=500$ | – | 2.148 | 319.947 | – | – | 0.279 | 318.232 | –
|  | – | (0.034) | (8.611) | – | – | (0.034) | (8.661) | –

Additionally, we report the performance of the estimator $\widehat{\beta}$ in
the NLFR for Models 1.1–1.3 in terms of $\mathrm{ASE}_{\beta}$ in Figure 1. It
can be seen that the $\mathrm{ASE}_{\beta}$ values of $\widehat{\beta}$
decrease and approach 0 with the increase of sample size $n$ in both cases of
$p=2$ and $p=5$, which confirms the consistency conclusion of Theorem 3.1.
When $p$ is larger, it is expected that a larger sample size is required for
the convergence of $\widehat{\beta}$.

Figure 1: The root $\mathrm{ASE}_{\beta}$ curves of the estimators
$\widehat{\beta}$ in NLFR for Models 1.1–1.3.

###  4.2 Fréchet regression for symmetric positive definite matrices

In this subsection, suppose that $\Omega$ is the set of symmetric positive
definite (SPD) matrices of a fixed dimension $m$ equipped with the Frobenius
metric $d_{F}$. Here, in evaluating the performance of various methods, we
simultaneously consider the Cholesky decomposition metric $d_{C}$ for a
comprehensive comparison. The Cholesky decomposition metric $d_{C}$ is defined
as follows: For $P_{1},P_{2}\in\Omega$, we can derive
$P_{1}=(P_{1}^{1/2})^{\mathrm{T}}P_{1}^{1/2}$ and
$P_{2}=(P_{2}^{1/2})^{\mathrm{T}}P_{2}^{1/2}$ based on Cholesky decomposition
with upper triangle matrices $P_{1}^{1/2}$ and $P_{2}^{1/2}$. Then the
Cholesky decomposition distance between $P_{1}$ and $P_{2}$ is defined as
$d_{C}(P_{1},P_{2})=\|P_{1}^{1/2}-P_{2}^{1/2}\|_{F}$, where $\|\cdot\|_{F}$ is
the Frobenius norm.

The regression function is specified by
$m_{\oplus}(x)=E(Y|X=x)=(U_{0}+\alpha^{\mathrm{T}}g(x))I_{m}+(V_{0}+\alpha^{\mathrm{T}}g(x))P_{m},$
where $U_{0}+\alpha^{\mathrm{T}}g(x)>0$, $V_{0}+\alpha^{\mathrm{T}}g(x)>0$,
$I_{m}$ denotes the $m\times m$ identity matrix and $P_{m}$ is a given
$m\times m$ positive semi-definite matrices. In our simulation, $P$ can be
given for each replication by the following procedure: (1) let $Z$ be a given
$m\times m$ matrix with independent $N(0,1)$ elements; (2) compute the
orthonormal basis vectors $z_{1},\ldots,z_{m^{\prime}}$ of $Z$ with
$m^{\prime}=\mathrm{rank}(Z)$; (3) obtain an $m\times m$ matrix $Q$ with
$Q=(z_{1},\ldots,z_{m^{\prime}})$ for $m^{\prime}=m$ and
$Q=(z_{1},\ldots,z_{m^{\prime}},0,\ldots)$ for $m^{\prime}<m$; (4)
$P=Q^{\mathrm{T}}DQ$ with $D=\mathrm{diag}\\{1,\ldots,m\\}$. For symmetric
positive definite matrices, the random response $Y$ is generated by adding
noise as follows: $Y=U^{2}I_{m}+VP_{m}$ with $U\mid X\sim
N((U_{0}+\alpha^{\mathrm{T}}g(X)-v_{1})^{1/2},v_{1})$ and $V\mid
X\sim\operatorname{Gamma}((V_{0}+\alpha^{\mathrm{T}}g(X))^{2}/v_{2},v_{2}/(V_{0}+\alpha^{\mathrm{T}}g(X)))$.
The transformation function $h(Y)$ in the SNLFR is chosen as
$h(Y)=\mathrm{trace}(Y-E(Y))$, where $\mathrm{trace}(\cdot)$ denotes the trace
of a matrix. In this way, we hope that $h(Y)$ retains the information of
$U^{2}$ and $V$ as much as possible. In simulation, $E(Y)$ is approximately
calculated by its empirical estimator based on a given sample. Further, an
optimal $c_{h}$ can be selected based on (2.12). In addition, by the
representation theorems, we have
$\beta^{(0)}=(U_{0}+\alpha^{\mathrm{T}}E[g(X)])I_{m}+(V_{0}+\alpha^{\mathrm{T}}E[g(X)])P_{m}$
and $\sigma=\mathrm{Cov}(X,g(X))\alpha(I_{m}+P_{m})$. Using
$m_{\oplus}(x)=\beta^{(0)}+\sigma^{\mathrm{T}}f(x)$, the link function $f(x)$
is determined by $\mathrm{Cov}(g(X),X)f(X)=g(X)-E[g(X)]$.

We next give Models 2.1-2.3 with the responses are the SPD matrices to assess
the performance of each method under different settings.

Model 2.1. (Linear Fréchet Regression) Let
$(X^{\prime(1)},\ldots,X^{\prime(p)})\sim N_{p}(0,\Sigma)$ with $\Sigma$
$=(0.5^{|i-j|})_{p\times p}$, and then set $X^{(j)}=2\Phi(X^{\prime(j)})-1$
for $j=1,\ldots,p$. Set $\alpha=(1,0,\ldots,0)^{\mathrm{T}}$,
$g(x)=(\beta^{\mathrm{T}}(x-\mu),0,\ldots,0)^{\mathrm{T}}$ and
$f(x)=(f_{1}(x,\beta),0,\ldots,0)^{\mathrm{T}}$. Let $p\in\\{2,5\\}$ and
$n\in\\{100,200\\}$. For $p=2$, $\beta=(1,-0.5)^{\mathrm{T}}$, $U_{0}=3$ and
$V_{0}=2$. For $p=5$, $\beta=(1,-0.5,2,1.5,-1)^{\mathrm{T}}$, $U_{0}=8$ and
$V_{0}=6.5$. The additional parameters are set as $m=3$, $v_{1}=1$,
$v_{2}=0.5$ and $\widetilde{n}=500$.

Model 2.2. (Generalized linear Fréchet Regression) Let
$(X^{(1)},\ldots,X^{(p)})\sim N_{p}(0,\Sigma)$ with
$\Sigma=(0.5^{|i-j|})_{p\times p}$ and
$g(x)=((\beta^{\mathrm{T}}(x-\mu)+1)^{2},0,\ldots,0)^{\mathrm{T}}$. The
additional parameters are set as $U_{0}=1.5$ and $V_{0}=0.5$. The other
settings are the same as in Model 2.1.

Model 2.3. (Nonlinear Fréchet Regression) Let
$g(x)=(g_{1}(x,\beta),0,\ldots,0)^{\mathrm{T}}$ with
$g_{1}(x,\beta)=\sum_{j=1}^{p_{1}}(\beta^{(j)}(x^{(j)}-\mu^{(j)})+1)^{2}+\sum_{j=p_{1}+1}^{p}\exp{(\beta^{(j)}(x^{(j)}-\mu^{(j)}))}$.
Set $p_{1}=1$ for $p=2$ and $p_{1}=3$ for $p=5$. The other settings are the
same as in Model 2.2.

To save space, the simulation results of various methods for Models 2.1–2.3
are recorded in Tables 3–4 (with Cholesky decomposition metric) and Figures
S1–S2 (with Frobenius metric) in Section S.1 of the Supplementary Materials.
For linear fréchet regression Model 2.1, under the Cholesky decomposition and
Frobenius metric, the NLFR and SNLFR perform with almost indistinguishable
differences and they are both superior to the LFR and the LoFR, especially for
the small sample size. Not surprisingly, we can observe from Tables 3–4 and
Figures S1–S2 that the NLFR is the best performer for the generalized linear
and nonlinear models, and its superiority is more apparent than other methods
when $p$ is larger. The SNLFR has a satisfactory performance when $p$ is
smaller and the difference between the NLFR and the SNLFR is insignificant
when the sample size is large for $p=5$.

Table 3:  The averaged $\mathrm{MSE}_{Y}$ and $\mathrm{MSE}_{m}$ of various
methods and the associated standard errors (in parenthesis) for Models 2.1–2.3
with $p=2$ under the Cholesky decomposition metric $d_{C}$.

|  |  | $\mathrm{MSE}_{Y}$ | $\mathrm{MSE}_{m}$
---|---|---|---|---
|  |  | SNLFR | NLFR | LFR | LoFR | SNLFR | NLFR | LFR | LoFR
| Model 2.1 | $n=100$ | 2.295 | 2.295 | 2.312 | 2.409 | 0.975 | 0.975 | 1.034 | 1.133
|  | (0.044) | (0.044) | (0.046) | (0.047) | (0.041) | (0.041) | (0.043) | (0.044)
| $n=200$ | 2.245 | 2.245 | 2.267 | 2.315 | 0.920 | 0.921 | 0.986 | 1.037
|  | (0.050) | (0.050) | (0.052) | (0.052) | (0.046) | (0.046) | (0.048) | (0.048)
| $n=500$ | 2.198 | 2.198 | 2.223 | 2.238 | 0.897 | 0.897 | 0.965 | 0.981
|  | (0.046) | (0.046) | (0.048) | (0.048) | (0.042) | (0.042) | (0.044) | (0.044)
| Model 2.2 | $n=100$ | 2.639 | 2.612 | 3.121 | 6.478 | 1.204 | 1.187 | 1.731 | 5.079
|  | (0.060) | (0.060) | (0.063) | (1.083) | (0.058) | (0.058) | (0.059) | (1.081)
| $n=200$ | 2.566 | 2.558 | 3.091 | 3.101 | 1.151 | 1.147 | 1.715 | 1.736
|  | (0.049) | (0.048) | (0.049) | (0.079) | (0.047) | (0.047) | (0.048) | (0.079)
| $n=500$ | 2.508 | 2.495 | 3.001 | 2.822 | 1.101 | 1.093 | 1.646 | 1.467
|  | (0.056) | (0.056) | (0.057) | (0.082) | (0.053) | (0.053) | (0.055) | (0.083)
| Model 2.3 | $n=100$ | – | 3.108 | 4.001 | 4.381 | – | 1.847 | 2.731 | 3.114
|  | – | (0.088) | (0.089) | (0.173) | – | (0.084) | (0.085) | (0.171)
| $n=200$ | – | 3.043 | 3.929 | 4.198 | – | 1.764 | 2.650 | 2.952
|  | – | (0.096) | (0.099) | (0.312) | – | (0.093) | (0.095) | (0.304)
| $n=500$ | – | 3.004 | 3.846 | 3.334 | – | 1.746 | 2.599 | 2.090
|  | — | (0.088) | (0.088) | (0.100) | – | (0.085) | (0.085) | (0.096)

Table 4:  The averaged $\mathrm{MSE}_{Y}$ and $\mathrm{MSE}_{m}$ of various
methods and the associated standard errors (in parenthesis) for Models 2.1–2.3
with $p=5$ under the Cholesky decomposition metric $d_{C}$.

|  |  | $\mathrm{MSE}_{Y}$ | $\mathrm{MSE}_{m}$
---|---|---|---|---
|  |  | GLFR | NLFR | LFR | LoFR | GLFR | NLFR | LFR | LoFR
| Model 2.1 | $n=100$ | 4.924 | 4.924 | 5.008 | – | 3.684 | 3.684 | 3.780 | –
|  | (0.166) | (0.166) | (0.170) | – | (0.160) | (0.160) | (0.164) | –
| $n=200$ | 4.840 | 4.849 | 4.913 | – | 3.593 | 3.600 | 3.679 | –
|  | (0.164) | (0.164) | (0.166) | – | (0.159) | (0.160) | (0.164) | –
| $n=500$ | 4.317 | 4.317 | 4.386 | – | 3.110 | 3.110 | 3.194 | –
|  | (0.165) | (0.165) | (0.169) | – | (0.161) | (0.161) | (0.165) | –
| Model 2.2 | $n=100$ | 17.489 | 8.138 | 28.766 | – | 15.780 | 6.806 | 26.874 | –
|  | (0.840) | (0.546) | (0.352) | – | (0.814) | (0.537) | (0.340) | –
| $n=200$ | 12.006 | 6.204 | 27.694 | – | 10.507 | 4.893 | 25.872 | –
|  | (0.721) | (0.293) | (0.296) | – | (0.708) | (0.289) | (0.289) | –
| $n=500$ | 9.150 | 7.058 | 27.590 | – | 7.730 | 5.733 | 25.694 | –
|  | (0.318) | (0.289) | (0.307) | – | (0.309) | (0.283) | (0.294) | –
| Model 2.3 | $n=100$ | – | 6.267 | 15.457 | – | – | 4.996 | 13.994 | –
|  | – | (0.222) | (0.230) | – | – | (0.217) | (0.229) | –
| $n=200$ | – | 6.078 | 15.811 | – | – | 4.815 | 14.333 | –
|  | – | (0.223) | (0.239) | – | – | (0.220) | (0.239) | –
| $n=500$ | – | 5.780 | 15.245 | – | – | 4.510 | 13.776 | –
|  | – | (0.219) | (0.234) | – | – | (0.215) | (0.227) | –

Furthermore, the performance of the estimator $\widehat{\beta}$ concerning the
$\mathrm{ASE}_{\beta}$ measure for $p=2$ and $p=5$ is exhibited in Figure S3
in Section S.2 of the Supplementary Materials. It can be seen that when the
sample size changes from 100 to 500, the $\mathrm{ASE}_{\beta}$ values
gradually become smaller and closer to 0 for both $p=2$ and $p=5$.

To sum up, these comprehensive numerical results show that the proposed SNLFR
and NLFR methods produce results that are significantly superior to commonly
competitive methods for the Fréchet regression models with generalized linear
and nonlinear structures. Also, these results can well verify the theoretical
conclusions about consistency given in the previous sections.

##  5 Real data analysis

In this section, we apply the proposed two approaches to the human mortality
dataset across countries. The age-at-death distributions are considered as
random object responses of interest, and the objective is to explore the
relationship between the age-at-death distributions and country-specific
covariates. The data are obtained from United Nations Databases
(http://data.un.org/) and UN World Population Prospects 2019 Databases
(https://population.un.org/wpp/Download) in the form of life tables for
$n=162$ countries in the period 2015-2020. For the sake of space, the detailed
analysis is presented in Section S.2 of the Supplementary Materials. From
these detailed analysis we can see that NLFR outperforms remarkably the
baseline LFR and LoFR methods in all scenarios. Also, the SNLFR performs
slightly inferior to the LoFR in the case of the small sample size of the
training set, but it will gradually outperform the LoFR with the increase of
sample size of the training set. Overall, the NLFR procedure is more
recommendable in practice.

##  6 Conclusions and future works

As shown in the introduction, the existing Fréchet regression methods identify
a model only in a linear framework. Even in the existing nonparametric and
semiparametric Fréchet regressions, the models are handled by locally or
approximately linear techniques, and the resulting models are locally or
approximately linear in nature. We in the previous sections proposed two types
of Fréchet regressions: generalized linear and nonlinear Fréchet regressions.
These frameworks can be utilized to fit the essentially nonlinear models.
Particularly, the SNLFR can return to the standard LFR when the underlying
model is indeed linear. Structurally, the nonlinear Fréchet objective
functions are defined by an unconditional expectation, and the related
Euclidean variable and the metric space object are separated into different
factors. In particular, under a Hilbert space, the NLFR and SNLFR have
explicit representations. These favorable structures make the estimation
procedures easy for theoretical analysis and numeral computation. Moreover,
the estimation consistency was established systematically in the previous
sections. Some comprehensive simulation studies and real data analysis were
given demonstrating that the estimation methods are easy to use, and the
performances of the estimations are significantly better than the competitors.

Although in a Hilbert space, the new methods are computationally simple, its
calculation burden is relatively heavy because the profile estimation is
involved when $Y$ is in a general metric space. Moreover, unlike the LFR, the
variable selection in SNLFR cannot be transformed into a ridge regression
framework. Thus, it is difficult to implement the variable selection in SNLFR
with high-dimensional covariates, and it is also difficult to extend the
methods to more general nonlinear models. These are interesting issues and are
worth further study in the future.

## References

  * Bennett and Kallus (2023) Bennett, A. and Kallus, N. (2023).  “The variational method of moments.”  _Journal of the Royal Statistical Society, Series B_ , 85, 810–841.
  * Bhattacharjee and Müller (2023) Bhattacharjee, S. and Müller, H. (2023).  “Single index Fréchet regression.”  _The Annals of Statistics_ , 51, 1770–1798.
  * Chen et al. (2020) Chen, Y., Gajardo, A., Fan, J., Zhong, Q., Dubey, P., Han, K., Bhattacharjee, S., and Müller, H. (2020).  “frechet: Statistical analysis for random objects and non-Euclidean data.”  _R package version 0.2.0._
  * Cook (1998) Cook, R. (1998).  “Regression graphics: Ideas for studying regressions through graphics.”  _John Wiley, New York_.
  * Cook and Ni (2005) Cook, R. and Ni, L. (2005).  “Sufficient dimension reduction via inverse regression: a minimum discrepancy approach.”  _Journal of the American Statistical Association_ , 100, 410–428.
  * Faraway (2014) Faraway, J.J. (2014).  “Regression for non-Euclidean data using distance matrices.”  _Journal of Applied Statistics_ , 41, 2342–2357.
  * Fréchet (1948) Fréchet, M. (1948).  “Les éléments aléatoires de nature quelconque dans un espace distancié.”  _Annales de l’institut Henri Poincaré_ , 10, 215–310.
  * Ghosal and Meiring (2023) Ghosal, A. and Meiring, W. (2023).  “Fréchet single index models for object response regression.”  _Electronic Journal of Statistics_ , 17, 1074–1112.
  * Hansen (1982) Hansen, L.P. (1982).  “Large sample properties of generalized method of moments estimators.”  _Econometrica_ , 50, 1029–1054.
  * Härdle et al. (2004) Härdle, W., Müller, M., Sperlich, S., and Werwatz, A. (2004).  “Nonparametric and semiparametric models.”  _Springer, Berlin Heidelberg_.
  * Li (1991) Li, K.C. (1991).  “Sliced inverse regression for dimension reduction.”  _Journal of the American Statistical Association_ , 86, 316–327.
  * Li and Duan (1989) Li, K.C. and Duan, N. (1989).  “Regression analysis under like violation.”  _The Annals of Statistics_ , 17, 1009–1052.
  * Marron and Alonso (2014) Marron, J.S. and Alonso, A.M. (2014).  “Overview of object oriented data analysis.”  _Biometrical Journal_ , 56, 732–753.
  * McCulloch and Searle (2001) McCulloch, C.E. and Searle, S.R. (2001).  “Generalized, linear, and mixed models.”  _John Wile & Sons, Inc_.
  * Petersen and Müller (2019) Petersen, A. and Müller, H. (2019).  “Fréchet regression for random objects with euclidean predictors.”  _The Annals of Statistics_ , 47, 691–719.
  * Schötz (2022) Schötz, C. (2022).  “Nonparametric regression in nonstandard spaces.”  _Electronic Journal of Statistics_ , 16, 4679–4741.
  * Severini and Wong (1992) Severini, T.A. and Wong, W.H. (1992).  “Profile likelihood and conditionally parametric models.”  _The Annals of Statistics_ , 20, 1768–1802.
  * Small and Mcleish (1994) Small, C.G. and Mcleish, D.F. (1994).  “Hilbert space methods in probability and statistical inference.”  _A Wiley-Interscience Publication, John Wiley & Sons, Inc._
  * Tucker et al. (2023) Tucker, D.C., Wu, Y., and Müller, H. (2023).  “Variable selection for global Fréchet regression.”  _Journal of the American Statistical Association_ , 118, 1023–1037.
  * Wang and Marron (2007) Wang, H. and Marron, J.S. (2007).  “Object oriented data analysis: Sets of trees.”  _The Annals of Statistics_ , 35, 1849–1873.
  * Wang et al. (2012) Wang, T., Xu, P., and Zhu, L. (2012).  “Non-convex penalized estimation in high-dimensional models with single-index structure.”  _Journal of Multivariate Analysis_ , 109, 221–235.
  * Wei (1998) Wei, B.C. (1998).  “Exponential families nonlinear model.”  _Springer-Verlag Singapore Pte. Ltd._
  * Yan et al. (2023) Yan, X., Zhang, X., and Zhao, P. (2023).  “Frequentist model averaging for global Fréchet regression.”  _arXiv preprint arXiv: 2309.01691_.