# Optimal convex $M$-estimation via score matching

Oliver Y. Feng∗,‡, Yu-Chun Kao†, Min Xu† and Richard J. Samworth‡

∗Department of Mathematical Sciences, University of Bath
†Department of Statistics, Rutgers University
‡Statistical Laboratory, University of Cambridge

###### Abstract

In the context of linear regression, we construct a data-driven convex loss
function with respect to which empirical risk minimisation yields optimal
asymptotic variance in the downstream estimation of the regression
coefficients. Our semiparametric approach targets the best decreasing
approximation of the derivative of the log-density of the noise distribution.
At the population level, this fitting process is a nonparametric extension of
score matching, corresponding to a log-concave projection of the noise
distribution with respect to the Fisher divergence. The procedure is
computationally efficient, and we prove that our procedure attains the minimal
asymptotic covariance among all convex $M$-estimators. As an example of a non-
log-concave setting, for Cauchy errors, the optimal convex loss function is
Huber-like, and our procedure yields an asymptotic efficiency greater than
$0.87$ relative to the oracle maximum likelihood estimator of the regression
coefficients that uses knowledge of this error distribution; in this sense, we
obtain robustness without sacrificing much efficiency. Numerical experiments
confirm the practical merits of our proposal.

##  1 Introduction

In linear models, the Gauss–Markov theorem is the primary justification for
the use of ordinary least squares (OLS) in settings where the Gaussianity of
our error distribution may be in doubt. It states that, provided the errors
have a finite second moment, OLS attains the minimal covariance among all
linear unbiased estimators; recent papers on this topic include Hansen,
(2022), Pötscher and Preinerstorfer, (2022) and Lei and Wooldridge, (2022).
Nevertheless, it is now understood that biased, non-linear estimators can
achieve lower mean squared error than OLS when the noise distribution is
appreciably non-Gaussian (Stein, 1956b, ; Hoerl and Kennard,, 1970; Zou and
Yuan,, 2008; Dümbgen et al.,, 2011). However, it remains unclear how best to
fit linear models in a computationally efficient and adaptive fashion, i.e.
without knowledge of the error distribution.

Consider a linear model where $Y_{i}=X_{i}^{\top}\beta_{0}+\varepsilon_{i}$
for $i=1,\dotsc,n$. Recall that an $M$-estimator of
$\beta_{0}\in\mathbb{R}^{d}$ based on a loss function
$\ell\colon\mathbb{R}\to\mathbb{R}$ is defined as an empirical risk minimiser

$\hat{\beta}\in\operatorname*{argmin}_{\beta\in\mathbb{R}^{d}}\frac{1}{n}\sum_{i=1}^{n}\ell(Y_{i}-X_{i}^{\top}\beta),$
(1)

provided that this exists. If $\ell$ is differentiable on $\mathbb{R}$ with
negative derivative $\psi=-\ell^{\prime}$, then
$\hat{\beta}\equiv\hat{\beta}_{\psi}$ solves the corresponding estimating
equations

$\frac{1}{n}\sum_{i=1}^{n}X_{i}\psi(Y_{i}-X_{i}^{\top}\hat{\beta}_{\psi})=0$
(2)

and is referred to as a $Z$-estimator. We study a random design setting in
which $(X_{1},Y_{1}),\dotsc,(X_{n},Y_{n})$ are independent and identically
distributed, with $X_{1},\dotsc,X_{n}$ being $\mathbb{R}^{d}$-valued
covariates that are independent of real-valued errors
$\varepsilon_{1},\dotsc,\varepsilon_{n}$ with density $p_{0}$. Suppose further
that $\mathbb{E}\\{X_{1}\psi(\varepsilon_{1})\\}=0$. This means that
$\hat{\beta}_{\psi}$ is _Fisher consistent_ in the sense that the population
analogue of (2) is satisfied by the true parameter $\beta_{0}$, i.e.
$\mathbb{E}\\{X_{1}\psi(Y_{1}-X_{1}^{\top}\beta_{0})\\}=0$. Under suitable
regularity conditions, including $\psi$ being differentiable and
$\mathbb{E}(X_{1}X_{1}^{\top})\in\mathbb{R}^{d\times d}$ being invertible, we
have

$\sqrt{n}(\hat{\beta}_{\psi}-\beta_{0})\overset{d}{\to}N_{d}\bigl{(}0,V_{p_{0}}(\psi)\cdot\\{\mathbb{E}(X_{1}X_{1}^{\top})\\}^{-1}\bigr{)}\quad\text{as
}n\to\infty,\quad\text{where
}V_{p_{0}}(\psi):=\frac{\mathbb{E}\psi^{2}(\varepsilon_{1})}{\\{\mathbb{E}\psi^{\prime}(\varepsilon_{1})\\}^{2}}$
(3)

(e.g. van der Vaart,, 1998, Theorems 5.21, 5.23 and 5.41). Since the
covariates and errors are assumed to be independent, they contribute
separately to the limiting covariance above (a special case of the ‘sandwich’
formula (Huber,, 1967; Young and Shah,, 2023)): the matrix
$\mathbb{E}(X_{1}X_{1}^{\top})^{-1}$ depends only on the covariate
distribution, whereas the scalar $V_{p_{0}}(\psi)$ depends on the loss
function $\ell$ (through $\psi=-\ell^{\prime}$) and on the error distribution.

If the errors $\varepsilon_{1},\varepsilon_{2},\dotsc$ have a known absolutely
continuous density $p_{0}$ on $\mathbb{R}$, then we can define the maximum
likelihood estimator $\hat{\beta}^{\mathrm{MLE}}$ by taking $\ell=-\log p_{0}$
in (1). In this case, $\psi=-\ell^{\prime}$ is the _score function (for
location)_ 111The score is usually defined as a function of a parameter
$\theta\in\mathbb{R}$ as the derivative of the log-likelihood; the link with
our terminology comes from considering the location model
$\\{p_{0}(\cdot+\theta):\theta\in\mathbb{R}\\}$, and evaluating the score at
the origin. $\psi_{0}:=(p_{0}^{\prime}/p_{0})\mathbbm{1}_{\\{p_{0}>0\\}}$.
Under appropriate regularity conditions (e.g. van der Vaart,, 1998, Theorem
5.39), including that the _Fisher information (for location)_
$i(p_{0}):=\int_{\mathbb{R}}\psi_{0}^{2}\,p_{0}=\int_{\\{p_{0}>0\\}}(p_{0}^{\prime})^{2}/p_{0}$
is finite, we have

$\sqrt{n}\,(\hat{\beta}^{\mathrm{MLE}}-\beta_{0})\overset{d}{\to}N_{d}\biggl{(}0,\frac{\\{\mathbb{E}(X_{1}X_{1}^{\top})\\}^{-1}}{i(p_{0})}\biggr{)}$
(4)

as $n\to\infty$. The limiting covariance
$\\{\mathbb{E}(X_{1}X_{1}^{\top})\\}^{-1}/i(p_{0})$ constitutes an efficiency
lower bound; see Remark 15 below. In fact, it can be seen directly that
$1/i(p_{0})$ is the smallest possible value of the asymptotic variance factor
$V_{p_{0}}(\psi)$ in the limiting covariance of
$\sqrt{n}(\hat{\beta}_{\psi}-\beta_{0})$ in (3). Indeed, by the Cauchy–Schwarz
inequality,

$V_{p_{0}}(\psi)=\frac{\int_{\mathbb{R}}\psi^{2}p_{0}}{\bigl{(}\int_{\mathbb{R}}\psi^{\prime}p_{0}\bigr{)}^{2}}=\frac{\int_{\mathbb{R}}\psi^{2}p_{0}}{\bigl{(}\int_{\mathbb{R}}\psi
p_{0}^{\prime}\bigr{)}^{2}}\geq\frac{1}{\int_{\\{p_{0}>0\\}}(p_{0}^{\prime})^{2}/p_{0}}=\frac{1}{i(p_{0})}\in(0,\infty),$
(5)

whenever the integration by parts in the second step is justified, and
equality holds if and only if there exists $\lambda\neq 0$ such that
$\psi(\varepsilon_{1})=\lambda\psi_{0}(\varepsilon_{1})$ almost surely. This
leads to an equivalent variational definition of the Fisher information; see
Huber and Ronchetti, (2009, Theorem 4.2), which we restate as Proposition 36
in Section 6.4. Thus, when (4) holds, $\hat{\beta}^{\mathrm{MLE}}$ has minimal
asymptotic covariance among all $Z$-estimators $\hat{\beta}_{\psi}$ for which
(3) is valid, with the score function $\psi_{0}$ being the optimal choice of
$\psi$.

Our goal in this work is to choose $\psi$ in a data-driven manner, such that
the corresponding loss function $\ell$ in (1) is convex, and such that the
scale factor $V_{p_{0}}(\psi)$ in the asymptotic covariance (3) of the
downstream estimator of $\beta_{0}$ is minimised. Convexity is a particularly
convenient property for a loss function, since for the purpose of
$M$-estimation, it leads to more tractable theory and computation. Indeed, the
empirical risk in (1) becomes convex in $\beta$, so its local minimisers are
global minimisers. In particular, when $\ell$ is also differentiable,
$\hat{\beta}_{\psi}$ is a $Z$-estimator satisfying (2) if and only if it is an
$M$-estimator satisfying (1). The existence, uniqueness and
$\sqrt{n}$-consistency of $\hat{\beta}_{\psi}$ are then guaranteed under
milder conditions on $\ell$ than for generic loss functions (Yohai and
Maronna,, 1979; Maronna and Yohai,, 1981; Portnoy,, 1985; Mammen,, 1989;
Arcones,, 1998; He and Shao,, 2000). Furthermore, an important practical
advantage is that we can compute $\hat{\beta}_{\psi}$ efficiently using convex
optimisation algorithms with guaranteed convergence (Boyd and Vandenberghe,,
2004, Chapter 9).

In view of the discussion above, our first main contribution in Section 2 is
to determine the optimal population-level convex loss function in the sense
described in the previous paragraph. For a uniformly continuous error density
$p_{0}$, this amounts to finding

$\psi_{0}^{*}\in\operatorname*{argmin}_{\psi\in\Psi_{\downarrow}(p_{0})}V_{p_{0}}(\psi),$
(6)

where $\Psi_{\downarrow}(p_{0})$ denotes the set of decreasing, right-
continuous functions $\psi$ satisfying
$\int_{\mathbb{R}}\psi^{2}p_{0}<\infty$. We will actually define the ratio
$V_{p_{0}}(\psi)$ in a slightly more general way than in (5) to allow us to
handle non-differentiable functions $\psi$. This turns out to be convenient
because, for instance, the robust Huber loss $\ell_{K}$ given by

$\ell_{K}(z):=\begin{cases}z^{2}/2\quad&\text{if }|z|\leq K\\\
K|z|-K^{2}/2\quad&\text{if }|z|>K\end{cases}$ (7)

for $K\in(0,\infty)$ has a non-differentiable negative derivative
$\psi_{K}:=-\ell_{K}^{\prime}$ satisfying $\psi_{K}(z)=(-K)\vee(-z)\wedge K$
for $z\in\mathbb{R}$.

In Section 2.1, we show that minimising $V_{p_{0}}(\cdot)$ over
$\Psi_{\downarrow}(p_{0})$ is equivalent to minimising the _score matching_
objective

$D_{p_{0}}(\psi):=\mathbb{E}\bigl{\\{}\psi^{2}(\varepsilon_{1})+2\psi^{\prime}(\varepsilon_{1})\bigr{\\}},$
(8)

over $\psi\in\Psi_{\downarrow}(p_{0})$, provided that we take appropriate care
in defining this expression when $\psi$ is not absolutely continuous. This
observation allows us to obtain an explicit characterisation of the
‘projected’ score function $\psi_{0}^{*}$ in terms of $p_{0}$ and its
distribution function $F_{0}$. Indeed, to obtain $\psi_{0}^{*}$ at
$z\in\mathbb{R}$, we can first consider $p_{0}\circ F_{0}^{-1}$ (whose domain
is $[0,1]$), then compute the right derivative of its least concave majorant,
before finally applying the resulting function to $F_{0}(z)$. The negative
antiderivative $\ell_{0}^{*}$ of $\psi_{0}^{*}$ is then the optimal convex
loss function we seek. An important property is that
$\mathbb{E}\psi_{0}^{*}(Y_{1}-X_{1}^{\top}\beta_{0})=0$, which ensures that
$\ell_{0}^{*}$ correctly identifies the estimand $\beta_{0}$ on the population
level; equivalently, $\hat{\beta}_{\psi_{0}^{*}}$ is Fisher consistent.

Figure 1: Top row: Plots of the score function $\psi_{0}$ (green) and
projected score function $\psi_{0}^{*}$ (blue); Bottom row: their respective
negative antiderivatives, namely the negative log-density $-\log p_{0}$
(green) and optimal convex loss function $\ell_{0}^{*}$ (blue), for each of
the following non-log-concave distributions (from left to right): (a)
Student’s $t_{2}$; (b) symmetrised Pareto (66) with $\sigma=2$ and $\alpha=3$;
(c) Gaussian mixture $0.4N(-2,1)+0.6N(2,1)$.

Note that $\hat{\beta}^{\mathrm{MLE}}$ is a convex $M$-estimator if and only
if $\ell=-\log p_{0}$ is convex, i.e. $p_{0}$ is log-concave, in which case
$\psi_{0}^{*}=\psi_{0}$ by (5). We will be especially interested in error
densities $p_{0}$ that are not log-concave, for which the efficiency lower
bound in (5) cannot be achieved by a convex $M$-estimator corresponding to a
decreasing function $\psi$. We interpret the minimum ratio
$V_{p_{0}}(\psi_{0}^{*})$ as an analogue of the inverse Fisher information,
serving as the crucial part of the efficiency lower bound for convex
$M$-estimators. To reinforce the link with score matching, we will see in
Section 2.2 that the density proportional to $e^{-\ell_{0}^{*}}$ is the best
log-concave approximation to $p_{0}$ with respect to the Fisher divergence
defined formally in (20) below. This is typically different from the well-
studied log-concave projection with respect to Kullback–Leibler divergence,
and indeed the latter may yield considerably suboptimal covariance for the
resulting convex $M$-estimator; see Proposition 6. In concrete examples where
$p_{0}$ has heavy tails (e.g. a Cauchy density) or is multimodal (e.g. a
mixture density), we compute closed-form expressions for the projected score
function and the optimal convex loss function in Section 2.3. In particular,
$\ell_{0}^{*}$ turns out to be a robust Huber-like loss function in the Cauchy
case. More generally, when the errors are heavy-tailed in the sense that their
(two-sided) hazard function is bounded, Lemma 4 shows that the projected score
function $\psi_{0}^{*}$ is bounded, in which case the corresponding convex
loss $\ell_{0}^{*}$ grows at most linearly in the tails and hence is robust to
outliers. A major advantage of our framework over the use of the Huber loss is
that it does not require the choice of a transition point $K$ (see (7) above)
between quadratic and linear regimes (which in a regression context amounts to
a choice of scale for the error distribution). In fact, the antitonic score
projection, and hence the Fisher divergence projection, is affine equivariant
(Remark 8), which reflects the fact that we optimise $V_{p_{0}}(\cdot)$ in (6)
over a class $\Psi_{\downarrow}(p_{0})$ that is closed under multiplication by
non-negative scalars.

In Section 3, we turn our attention to a linear regression setting where the
error density $p_{0}$ is unknown. To ensure that $\beta_{0}$ is identifiable,
we assume either that $p_{0}$ is symmetric (Section 3.1) or that the model
contains an explicit intercept term (Section 3.2). We aim to construct a
semiparametric $M$-estimator of $\beta_{0}$ that achieves minimal covariance
among all convex $M$-estimators, but since $p_{0}$ and hence the optimal loss
function $\ell_{0}^{*}$ are unknown, we seek to estimate $\beta_{0}$ and
$\psi_{0}^{*}$ simultaneously. We first show that on the population level, the
pair $(\beta_{0},\psi_{0}^{*})$ solves a joint score matching problem subject
to a Fisher consistency constraint; see (28) and (29). This motivates an
alternating optimisation procedure where we start with an arbitrary
initialiser $\bar{\beta}_{n}$, and compute a kernel density estimate of the
error distribution based on the residuals. We can then apply the linear-time
Pool Adjacent Violators Algorithm (PAVA) to obtain the projected score
function of the density estimate, before minimising its negative
antiderivative using Newton optimisation techniques to yield an updated
estimator. This process could then be iterated to convergence, but if we
initialise with a $\sqrt{n}$-consistent pilot estimator $\bar{\beta}_{n}$,
then one iteration of the alternating algorithm above suffices for our
theoretical guarantees, and moreover it ensures that the procedure is
computationally efficient. We prove that a three-fold cross-fitting version of
our algorithm (with the different steps computed on different folds) yields an
estimator $\hat{\beta}_{n}$ that is $\sqrt{n}$-consistent and asymptotically
normal, with limiting covariance attaining our efficiency lower bound for
convex $M$-estimators. Consistent estimation of the information quantity
$1/V_{p_{0}}(\psi_{0}^{*})$ is straightforward using our nonparametric score
matching procedure, so combining this with our asymptotic distributional
result for $\hat{\beta}_{n}$, we can then perform inference for $\beta_{0}$
(Section 3.3).

Section 4 is devoted to a numerical study of the empirical performance and
computational efficiency of our antitonic score matching estimator. These
corroborate our theoretical findings: our proposed approach achieves smaller
estimation error (sometimes dramatically smaller) compared with alternatives
such as OLS, the least absolute deviation (LAD) estimator, a semiparametric
one-step estimator, and a semiparametric $M$-estimator based on the log-
concave MLE of the noise distribution. Moreover, the corresponding confidence
sets for $\beta_{0}$ are smaller, while retaining nominal coverage. Finally,
we perform a runtime analysis to show that the improved statistical
performance comes without sacrificing computational scalability.

The proofs of all results in Sections 2 and Section 3 are given in the
appendix in Sections 6.1 and 6.3 respectively. The appendix (Section 6) also
contains additional examples for Section 2 (Section 6.2) and auxiliary results
(Sections 6.4 and 6.5).

###  1.1 Related work

Score matching (Hyvärinen,, 2005; Lyu,, 2012) is an estimation method designed
for statistical models where the likelihood is only known up to a
normalisation constant (e.g. a partition function) that may be infeasible to
compute; see the recent tutorial by Song and Kingma, (2021) on ‘energy-based’
models. Instead of maximising an approximation to the likelihood, score
matching circumvents this issue altogether by estimating the derivative of a
log-density, i.e. the score function. More precisely, given a differentiable
density $p_{0}$ on $\mathbb{R}^{d}$ with score function $\psi_{0}:=(\nabla
p_{0}/p_{0})\mathbbm{1}_{\\{p_{0}>0\\}}$, the population version of the
procedure aims to minimise

$\mathbb{E}_{p_{0}}\bigl{(}\|{\psi(\varepsilon)-\psi_{0}(\varepsilon)}\|^{2}\bigr{)}$
(9)

over a suitable class $\Psi$ of differentiable functions
$\psi\equiv(\psi_{1},\dotsc,\psi_{d})\colon\mathbb{R}^{d}\to\mathbb{R}^{d}$,
where $\varepsilon\sim p_{0}$. Hyvärinen, (2005) used integration by parts to
show that it is equivalent to minimise

$D_{p_{0}}(\psi):=\mathbb{E}\bigl{\\{}\|{\psi(\varepsilon)}\|^{2}+2(\nabla\cdot\psi)(\varepsilon)\bigr{\\}}$
(10)

over $\psi\in\Psi$, where
$\nabla\cdot\psi:=\sum_{j=1}^{d}\partial\psi_{j}/\partial x_{j}$. The score
matching estimator based on data $\varepsilon_{1},\dotsc,\varepsilon_{n}$ in
$\mathbb{R}^{d}$ is then defined as a minimiser of the empirical analogue
$\hat{D}_{n}(\psi):=n^{-1}\sum_{i=1}^{n}\\{\|{\psi(\varepsilon_{i})}\|^{2}+2(\nabla\cdot\psi)(\varepsilon_{i})\\}$
over $\psi\in\Psi$; see also Cox, (1985). Such estimators are important in the
context of Langevin Monte Carlo (Parisi,, 1981; Roberts and Tweedie,, 1996;
Betancourt et al.,, 2017; Cheng et al.,, 2018) and diffusion models (Li et
al.,, 2023). The appearance of the score function in the underlying (reverse-
time) stochastic differential equations can be related to Tweedie’s formula,
which underpins empirical Bayes denoising (Efron,, 2011; Derenski et al.,,
2023).

Likelihood maximisation corresponds to distributional approximation with
respect to the Kullback–Leibler divergence; on the other hand, score matching
seeks to minimise the Fisher divergence (Johnson,, 2004, Section 1.3) from a
class of densities to the target $p_{0}$, in view of the equivalence between
the optimisation objectives (9) and (10); see (18) below. Sriperumbudur et
al., (2017) studied infinite-dimensional exponential families indexed by
reproducing kernel Hilbert spaces, and proposed and analysed a density
estimator that minimises a penalised empirical Fisher divergence. Koehler et
al., (2022) used isoperimetric inequalities to investigate the statistical
efficiency of score matching relative to maximum likelihood, thereby
quantifying the effect of eliminating normalisation factors. Lyu, (2012)
observed that Fisher divergence and Kullback–Leibler divergence are related by
an analogue of de Bruijn’s identity (Johnson,, 2004, Appendix C; Cover and
Thomas,, 2006, Section 17.7), which links Fisher information and Shannon
entropy. From an information-theoretic perspective, Johnson and Barron, (2004)
proved central limit theorems that establish convergence in Fisher divergence
to a limiting Gaussian distribution. Ley and Swan, (2013) extended Stein’s
method to derive information inequalities that bound a variety of integral
probability distances in terms of the Fisher divergence.

Score matching has been generalised in different directions and applied to a
variety of statistical problems including graphical modelling (e.g.
Hyvärinen,, 2007; Vincent,, 2011; Lyu,, 2012; Mardia et al.,, 2016; Song et
al.,, 2020; Yu et al.,, 2020, 2022; Lederer and Oesting,, 2023; Benton et
al.,, 2024), where it exhibits excellent empirical performance while being
computationally superior to full likelihood approaches. In particular, score-
based algorithms for generative modelling, via Langevin dynamics (Song and
Ermon,, 2019) and diffusion models (Song et al.,, 2021), have achieved
remarkable success in machine learning tasks such as the reconstruction,
inpainting and artificial generation of images; see e.g. Jolicoeur-Martineau
et al., (2020), De Bortoli et al., (2022) and many other references therein.
In these applications, score matching is applied to a class of functions
parametrised by the weights of a deep neural network. On the other hand,
different statistical considerations lead us to develop a nonparametric
extension of score matching in Section 2, which we use to construct data-
driven convex loss functions for efficient semiparametric estimation. We see
that it is by minimising the Fisher divergence instead of the Kullback–Leibler
divergence to the error distribution that one obtains a convex $M$-estimator
with minimal asymptotic variance.

The framework in Section 3.1 includes as a special case the classical location
model in which we observe $Y_{i}=\theta_{0}+\varepsilon_{i}$ for
$i=1,\dotsc,n$, where $\theta_{0}\in\mathbb{R}$ is the parameter of interest
and $\varepsilon_{1},\dotsc,\varepsilon_{n}$ are independent errors with an
unknown density $p_{0}$ that is symmetric about 0. Starting from the seminal
paper of Stein, 1956a , a series of works (e.g. van Eeden,, 1970; Stone,,
1975; Beran,, 1978; Bickel,, 1982; Schick,, 1986; Faraway,, 1992; Dalalyan et
al.,, 2006; Gupta et al.,, 2023) showed that adaptive, asymptotically
efficient estimators of $\theta_{0}$ can be constructed; see also Doss and
Wellner, (2019) and Laha, (2021) for approaches based on the further
assumption that $p_{0}$ is log-concave. Many of these traditional
semiparametric procedures have drawbacks that limit their practical utility.
In particular, the estimated likelihood may have multiple local optima and it
may be difficult to guarantee convergence of an optimisation algorithm to a
global maximum (van der Vaart,, 1998, Example 5.50). This is one of the
reasons why prior works often study a one-step estimator resulting from a
single iteration of Newton’s method (Bickel,, 1975; Jin,, 1990; Mammen and
Park,, 1997; Laha,, 2021), rather than full likelihood maximisation, though
finite-sample performance may remain poor and sensitive to tuning (see Section
4). By contrast, our focus is not on classical semiparametric adaptive
efficiency per se; instead, we directly study the theoretical properties of a
minimiser of the empirical risk with respect to an estimated loss function,
whose convexity ensures that the estimator can computed efficiently by
iterating gradient descent or Newton’s method to convergence.

Recently, Kao et al., (2023) constructed a location $M$-estimator that can
adaptively attain rates of convergence faster than $n^{-1/2}$ when the
symmetric error density is compactly supported and suitably irregular (e.g.
discontinuous at the boundary of its support). They considered
$\ell^{q}$-location estimators
$\hat{\theta}_{q}:=\operatorname*{argmin}_{\theta\in\mathbb{R}}\sum_{i=1}^{n}|Y_{i}-\theta|^{q}$
based on univariate observations $Y_{1},\dotsc,Y_{n}$, and used Lepski’s
method to select an exponent $\hat{q}\in[2,\infty)$ that minimises a proxy for
the asymptotic variance of $\hat{\theta}_{q}$. The resulting estimator
$\hat{\theta}_{\hat{q}}$ is shown to be minimax optimal up to poly-logarithmic
factors, and the procedure is extended to linear regression models with
unknown symmetric errors. By comparison with Kao et al., (2023), we study
‘regular’ regression models where the Fisher information is finite and minimax
rates faster than $n^{-1/2}$ are impossible to achieve. We aim to minimise the
asymptotic variance as an end in itself, over the entire nonparametric class
of convex loss functions rather than $\ell^{q}$-loss functions specifically.

We finally mention the connections between our work and robust statistics,
which deals with heavy-tailed noise distributions and data that may be
contaminated by random or adversarial outliers. As mentioned previously,
robust loss functions are designed to be tolerant to such data corruption;
examples include the Huber loss (7), a two-parameter family of loss functions
considered by Barron, (2019), and the antiderivatives of Catoni’s influence
functions (Catoni,, 2012). The Huber loss functions $\ell_{K}$ originally
arose as solutions to the following minimax asymptotic variance problem: for
every $\epsilon\in(0,1)$, there exists a unique $K\equiv K_{\epsilon}>0$ such
that

$\psi_{K}\equiv-\ell_{K}^{\prime}=\operatorname*{argmin}_{\psi\in\Psi}\sup_{P\in\mathcal{P}_{\epsilon}^{\mathrm{sym}}(\Phi)}\underbrace{\frac{\int_{\mathbb{R}}\psi^{2}\,dP}{\int_{\mathbb{R}}\psi^{\prime}\,dP}}_{=:\,V_{P}(\psi)},$

where $\Psi$ consists of all ‘sufficiently regular’
$\psi\colon\mathbb{R}\to\mathbb{R}$, and the symmetric
$\epsilon$-contamination neighbourhood
$\mathcal{P}_{\epsilon}^{\mathrm{sym}}(\Phi)$ contains all univariate
distributions of the form $P=(1-\epsilon)N(0,1)+\epsilon Q$ for some symmetric
distribution $Q$. The pioneering paper of Huber, (1964) also developed
variational theory for minimising $\sup_{P\in\mathcal{P}}V_{P}(\cdot)$ more
generally when $\mathcal{P}$ is a convex class of distributions, such as an
$\epsilon$-contamination or a Kolmogorov neighbourhood of a symmetric log-
concave density (Huber and Ronchetti,, 2009, Section 4.5). See Donoho and
Montanari, (2015) for a high-dimensional extension of this line of work. An
alternative to the Huber loss that seeks robustness without serious efficiency
loss relative to OLS is the composite quantile regression (CQR) estimator of
Zou and Yuan, (2008); in fact, our approach is always at least as efficient as
CQR (see Lemma 17). Other recent papers on robust convex $M$-estimation
include Chinot et al., (2020) and Brunel, (2023); see also the notes on robust
statistical learning theory by Lerasle, (2019).

More closely related to our optimisation problem (6) is the work of Hampel,
(1974) on optimal $B$-robust estimators, which have minimal asymptotic
variance subject to an upper bound on the gross error sensitivity (Hampel et
al.,, 2011, Section 2.4). In our linear regression setting with
$\varepsilon_{1}\sim p_{0}$, this amounts to

$\text{minimising }V_{p_{0}}(\psi)\text{ over all `regular' $\psi$ such that
}\int_{\mathbb{R}}\psi\,p_{0}=0\text{ and }\sup_{z\in\mathbb{R}}|\psi(z)|\leq
b$ (11)

for some suitable $b>0$ (van der Vaart,, 1998, Example 5.29; Hampel et al.,,
2011, p. 121 and Section 2.5d). In particular, when $p_{0}$ is a standard
Gaussian density, $\psi_{K}$ is again optimal for some $K\equiv K_{b}>0$ that
depends non-linearly on $b$. By contrast with (6) however, the Fisher
consistency condition $\mathbb{E}\psi(\varepsilon_{1})=0$ must be explicitly
included as a constraint in (11), and moreover the $L^{\infty}$ bound on
$\psi$ means that the set of feasible $\psi$ is not closed under non-negative
scalar multiplication. Consequently, the resulting optimal location
$M$-estimators are generally not scale invariant (Hampel et al.,, 2011, p.
105). In robust regression, adaptive selection of scale parameters is a non-
trivial problem (e.g. van der Vaart,, 1998, Section 5.4; Huber and Ronchetti,,
2009, Section 7.7; Loh,, 2021); see also Figure 4 below. Finally, we mention
that in a proportional asymptotic regime where $n/d\to\kappa\in(1,\infty)$ for
a sequence of linear models with log-concave errors and independent Gaussian
covariates, Bean et al., (2013) derived the (unpenalised) convex $M$-estimator
with minimal expected out-of-sample prediction error. In this setting, the
optimisation objective is no longer $V_{p_{0}}(\psi)$ but instead the solution
to a pair of non-linear equations involving the proximal operator of the
convex loss function (El Karoui et al.,, 2013).

###  1.2 Notation

Throughout this paper, we will adopt the convention $0/0:=0$ and write
$[n]:=\\{1,\dotsc,n\\}$ for $n\in\mathbb{N}$. For a function
$f\colon\mathbb{R}\to\mathbb{R}$, let
$\|{f}\|_{\infty}:=\sup_{z\in\mathbb{R}}|f(z)|$. Recall that $f$ is
_symmetric_ (i.e. even) if $f(z)=f(-z)$ for all $z\in\mathbb{R}$, and
_antisymmetric_ (i.e. odd) if $f(z)=-f(-z)$ for all $z\in\mathbb{R}$. For an
open set $U\subseteq\mathbb{R}$, we say that $f\colon U\to\mathbb{R}$ is
_locally absolutely continuous_ on $U$ if it is absolutely continuous on every
compact interval $I\subseteq U$. Equivalently, there exists a measurable
function $g\colon U\to\mathbb{R}$ such that for every compact subinterval
$I\subseteq U$, we have $\int_{I}|g|<\infty$ and
$f(z_{2})=f(z_{1})+\int_{z_{1}}^{z_{2}}g$ for all $z_{1},z_{2}\in I$. In this
case, $f$ is differentiable Lebesgue almost everywhere on $U$, with
$f^{\prime}=g$ almost everywhere.

Given a Borel probability measure $P$ on $\mathbb{R}$, we write $L^{2}(P)$ for
the set of all Lebesgue measurable functions $f$ on $\mathbb{R}$ such that
$\|{f}\|_{L^{2}(P)}:=\bigl{(}\int_{\mathbb{R}}f^{2}\,dP\bigr{)}^{1/2}<\infty$.
Denote by $\langle f,g\rangle_{L^{2}(P)}:=\int_{\mathbb{R}}fg\,dP$ the
$L^{2}(P)$-inner product of $f,g\in L^{2}(P)$. For measures $\mu,\nu$ on a
general measurable space $(\mathcal{X},\mathcal{A})$, we say that $\mu$ is
absolutely continuous with respect to $\nu$, and write $\mu\ll\nu$, if
$\mu(A)=0$ whenever $\nu(A)=0$ for $A\in\mathcal{A}$. The notation
$\mu\centernot{\ll}\nu$ indicates that $\mu$ is not absolutely continuous with
respect to $\nu$.

For a function $F\colon[0,1]\to\mathbb{R}$, we write $\hat{F}$ for its least
concave majorant on $[0,1]$. Denote by $F^{(\mathrm{L})}(u)$ and
$F^{(\mathrm{R})}(u)$ respectively the left and right derivatives of $F$ at
$u\in[0,1]$, whenever these are well-defined. Given an integrable function
$f\colon(0,1)\to\mathbb{R}$ with antiderivative $F\colon[0,1]\to\mathbb{R}$
given by $F(u):=\int_{0}^{u}f$, define
$\widehat{\mathcal{M}}_{\mathrm{R}}f\colon[0,1]\to[-\infty,\infty]$ by

$(\widehat{\mathcal{M}}_{\mathrm{R}}f)(u):=\begin{cases}\hat{F}^{(\mathrm{R})}(u)\;&\text{for
}u\in[0,1)\\\ \hat{F}^{(\mathrm{L})}(1)\;&\text{for }u=1,\end{cases}$

so that $(\widehat{\mathcal{M}}_{\mathrm{R}}f)(1)=\lim_{u\nearrow
1}(\widehat{\mathcal{M}}_{\mathrm{R}}f)(u)$ by Rockafellar, (1997, Theorem
24.1). Furthermore, define
$\widehat{\mathcal{M}}_{\mathrm{L}}f\colon[0,1]\to[-\infty,\infty]$ by
$(\widehat{\mathcal{M}}_{\mathrm{L}}f)(u):=(\widehat{\mathcal{M}}_{\mathrm{R}}g)(1-u)$
for $u\in[0,1]$, where $g(u):=f(1-u)$ for all such $u$.

##  2 The antitonic score projection

###  2.1 Construction and basic properties

The aim of this section is to define formally and solve the optimisation
problem (6) that yields the minimal asymptotic covariance of the regression
$M$-estimator in (3). Let $P_{0}$ be a probability measure on $\mathbb{R}$
with a uniformly continuous density $p_{0}$, which necessarily satisfies
$p_{0}(\pm\infty):=\lim_{z\to\pm\infty}p_{0}(z)=0$. Letting
$\operatorname{supp}p_{0}:=\\{z\in\mathbb{R}:p_{0}(z)>0\\}$, define
$\mathcal{S}_{0}\equiv\mathcal{S}(p_{0}):=\bigl{(}\inf(\operatorname{supp}p_{0}),\sup(\operatorname{supp}p_{0})\bigr{)}$,
which is the smallest open interval that contains $\operatorname{supp}p_{0}$.
We write $\Psi_{\downarrow}(p_{0})$ for the set of all $\psi\in L^{2}(P_{0})$
that are decreasing and right-continuous. Observe that
$\Psi_{\downarrow}(p_{0})$ is a convex cone, i.e.
$c_{1}\psi_{1}+c_{2}\psi_{2}\in\Psi_{\downarrow}(p_{0})$ whenever
$\psi_{1},\psi_{2}\in\Psi_{\downarrow}(p_{0})$ and $c_{1},c_{2}\geq 0$.
Moreover, every $\psi\in\Psi_{\downarrow}(p_{0})$ is necessarily finite-valued
on $\mathcal{S}_{0}$, so the corresponding Lebesgue–Stieltjes integral
$\int_{\mathcal{S}_{0}}p_{0}\,d\psi\in[-\infty,0]$ is well-defined.

For $\psi\in\Psi_{\downarrow}(p_{0})$ with
$\int_{\mathbb{R}}\psi^{2}\,dP_{0}>0$, let

$V_{p_{0}}(\psi):=\frac{\int_{\mathbb{R}}\psi^{2}\,dP_{0}}{\bigl{(}\int_{\mathcal{S}_{0}}p_{0}\,d\psi\bigr{)}^{2}}\in[0,\infty],$
(12)

where we have modified the denominator in (5) to extend the original
definition of the asymptotic variance factor to non-differentiable functions
in $\Psi_{\downarrow}(p_{0})$ such as $z\mapsto-\operatorname{sgn}(z)$. As a
first step towards minimising $V_{p_{0}}(\psi)$ over
$\psi\in\Psi_{\downarrow}(p_{0})$, note that
$V_{p_{0}}(c\psi)=V_{p_{0}}(\psi)$ for every $c>0$, so any minimiser is at
best unique up to a positive scalar. Ignoring unimportant edge cases where the
denominator in (12) is zero or infinity, our optimisation problem can
therefore be formulated as a constrained minimisation of the numerator in (12)
subject to the denominator being equal to 1. This motivates the definition of
the Lagrangian

$D_{p_{0}}(\psi,\lambda):=\int_{\mathbb{R}}\psi^{2}\,dP_{0}+2\lambda\int_{\mathcal{S}_{0}}p_{0}\,d\psi\in[-\infty,\infty)$
(13)

for $\psi\in\Psi_{\downarrow}(p_{0})$ and $\lambda\geq 0$. If $\psi$ is
locally absolutely continuous on $\mathcal{S}_{0}$ with derivative
$\psi^{\prime}$ Lebesgue almost everywhere, then

$D_{p_{0}}(\psi,\lambda)=\int_{\mathbb{R}}\psi^{2}\,dP_{0}+2\lambda\int_{\mathcal{S}_{0}}\psi^{\prime}p_{0}=\int_{\mathbb{R}}(\psi^{2}+2\lambda\psi^{\prime})\,dP_{0}=\mathbb{E}\bigl{(}\psi^{2}(\varepsilon_{1})+2\lambda\psi^{\prime}(\varepsilon_{1})\bigr{)}$
(14)

when $\varepsilon_{1}\sim P_{0}$, which we recognise as the score matching
objective (8) in the introduction when $\lambda=1$.

The formal link between $V_{p_{0}}(\cdot)$ and $D_{p_{0}}(\cdot\,,\lambda)$ is
that for $\psi\in\Psi_{\downarrow}(p_{0})$ with
$\int_{\mathbb{R}}\psi^{2}\,dP_{0}>0$, we have
$\int_{\mathcal{S}_{0}}p_{0}\,d\psi\leq 0$ and
$c\psi\in\Psi_{\downarrow}(p_{0})$ for all $c\geq 0$, so

$\inf_{c\geq 0}D_{p_{0}}(c\psi,\lambda)=\inf_{c\geq
0}\,\Bigl{(}c^{2}\int_{\mathbb{R}}\psi^{2}\,dP_{0}+2c\lambda\int_{\mathcal{S}_{0}}p_{0}\,d\psi\Bigr{)}=-\frac{\lambda^{2}\bigl{(}\int_{\mathcal{S}_{0}}p_{0}\,d\psi\bigr{)}^{2}}{\int_{\mathbb{R}}\psi^{2}\,dP_{0}}=-\frac{\lambda^{2}}{V_{p_{0}}(\psi)}$
(15)

for every $\lambda\geq 0$. Thus, minimising $V_{p_{0}}(\cdot)$ over
$\Psi_{\downarrow}(p_{0})$ is equivalent to minimising
$D_{p_{0}}(\cdot\,,\lambda)$ up to a scalar multiple, but
$D_{p_{0}}(\cdot,\lambda)$ is a convex function that is more tractable than
$V_{p_{0}}(\cdot)$. Both Cox, (1985, Proposition 1) and Theorem 2 below
indicate that $\lambda=1$ is the ‘canonical’ choice of Lagrange multiplier in
(13), so we define $D_{p_{0}}(\psi):=D_{p_{0}}(\psi,1)$ for
$\psi\in\Psi_{\downarrow}(p_{0})$.

By exploiting this connection with score matching together with ideas from
monotone function estimation, we prove in Theorem 2 below that the solution to
our asymptotic variance minimisation problem is the function $\psi_{0}^{*}$
that we construct explicitly in the following lemma.

######  Lemma 1.

Let $P_{0}$ be a distribution with a uniformly continuous density $p_{0}$ on
$\mathbb{R}$. Let $F_{0}\colon[-\infty,\infty]\to[0,1]$ be the corresponding
distribution function, and for $u\in[0,1]$, define

$F_{0}^{-1}(u):=\inf\\{z\in[-\infty,\infty]:F_{0}(z)\geq
u\\}\quad\text{and}\quad J_{0}(u):=(p_{0}\circ F_{0}^{-1})(u).$

Then both $J_{0}$ and its least concave majorant $\hat{J}_{0}$ on $[0,1]$ are
continuous, with $p_{0}=J_{0}\circ F_{0}$ on $\mathbb{R}$, and

$\psi_{0}^{*}:=\hat{J}_{0}^{(\mathrm{R})}\circ F_{0}$

is decreasing and right-continuous as a function from $\mathbb{R}$ to
$[-\infty,\infty]$, provided that we set
$\hat{J}_{0}^{(\mathrm{R})}(1):=\lim_{u\nearrow
1}\hat{J}_{0}^{(\mathrm{R})}(u)$. Moreover, $\psi_{0}^{*}(z)\in\mathbb{R}$ if
and only if $z\in\mathcal{S}_{0}$.

We refer to $J_{0}$ as the _density quantile function_ (Parzen,, 1979; Jones,,
1992). In the case where $p_{0}$ is a standard Cauchy density, Figure 2
presents a visualisation of $J_{0}$ and its least concave majorant
$\hat{J}_{0}$, as well as the corresponding score functions
$\psi_{0}=p_{0}^{\prime}/p_{0}$ and $\psi_{0}^{*}$.

Figure 2: _Left_ : The density quantile function $J_{0}$ and its least concave
majorant $\hat{J}_{0}$ for a standard Cauchy density. _Right_ : The
corresponding score functions $\psi_{0}$ and $\psi_{0}^{*}$.

######  Theorem 2.

In the setting of Lemma 1, the following statements hold.

  1. (a)

$\int_{\mathbb{R}}\psi_{0}^{*}\,dP_{0}=0$.

  2. (b)

Let $i^{*}(p_{0}):=\int_{\mathbb{R}}(\psi_{0}^{*})^{2}\,dP_{0}$. Then
$\inf_{\psi\in\Psi_{\downarrow}(p_{0})}D_{p_{0}}(\psi,\lambda)=-\lambda^{2}i^{*}(p_{0})$
for every $\lambda>0$.

  3. (c)

Suppose that $i^{*}(p_{0})<\infty$. Then for each $\lambda>0$, the function
$\lambda\psi_{0}^{*}$ is the unique minimiser of $D_{p_{0}}(\cdot\,,\lambda)$
over $\Psi_{\downarrow}(p_{0})$. Moreover, for every
$\psi\in\Psi_{\downarrow}(p_{0})$ such that
$\int_{\mathbb{R}}\psi^{2}\,dP_{0}>0$, we have

$V_{p_{0}}(\psi)\geq
V_{p_{0}}(\psi_{0}^{*})=\frac{1}{i^{*}(p_{0})}\in(0,\infty),$ (16)

with equality if and only if $\psi=\lambda\psi_{0}^{*}$ for some $\lambda>0$.

  4. (d)

Assume further that $p_{0}$ is absolutely continuous on $\mathbb{R}$ with
derivative $p_{0}^{\prime}$ Lebesgue almost everywhere, corresponding score
function222Our convention $0/0=0$ means that
$\psi_{0}=(p_{0}^{\prime}/p_{0})\mathbbm{1}_{\\{p_{0}>0\\}}$.
$\psi_{0}:=p_{0}^{\prime}/p_{0}$ and Fisher information
$i(p_{0}):=\int_{\mathbb{R}}\psi_{0}^{2}\,p_{0}$. Then

$\psi_{0}^{*}=\widehat{\mathcal{M}}_{\mathrm{R}}(\psi_{0}\circ
F_{0}^{-1})\circ F_{0}$ (17)

and $0<i^{*}(p_{0})\leq i(p_{0})$, with equality if and only if $p_{0}$ is
log-concave. In particular, if $i(p_{0})<\infty$, then the conclusions of (c)
hold.

Some remarks are in order here. As mentioned in the introduction, Theorem 2(a)
ensures the Fisher consistency of the regression $Z$-estimator
$\hat{\beta}_{\psi_{0}^{*}}$ defined in (2). This reflects the fact that
$\psi+c\in\Psi_{\downarrow}(p_{0})$ whenever $\psi\in\Psi_{\downarrow}(p_{0})$
and $c\in\mathbb{R}$; see (46) and the first-order stationarity condition
(47). Moreover, (53) in the proof of Theorem 2 _(d)_ shows that

$D_{p_{0}}(\psi,\lambda)=\int_{\mathbb{R}}(\psi-\lambda\psi_{0})^{2}\,dP_{0}-\int_{\mathbb{R}}(\lambda\psi_{0})^{2}\,dP_{0}=\|{\psi-\lambda\psi_{0}}\|_{L^{2}(P_{0})}^{2}-\lambda^{2}i(p_{0})$

for all $\psi\in\Psi_{\downarrow}(p_{0})$, so if $i(p_{0})<\infty$, then
Theorem 2 _(c)_ ensures that

$\lambda\psi_{0}^{*}\in\operatorname*{argmin}_{\psi\in\Psi_{\downarrow}(p_{0})}D_{p_{0}}(\psi,\lambda)=\operatorname*{argmin}_{\psi\in\Psi_{\downarrow}(p_{0})}\|{\psi-\lambda\psi_{0}}\|_{L^{2}(P_{0})}^{2}.$
(18)

Thus, in the terminology of Section 6.5, $\lambda\psi_{0}^{*}$ is a version of
the $L^{2}(P_{0})$-antitonic333Antitonic means decreasing, in contrast to
isotonic (increasing) (Groeneboom and Jongbloed,, 2014, Section 2.1).
projection of $\lambda\psi_{0}$ onto $\Psi_{\downarrow}(p_{0})$. Indeed, the
explicit representation (17) of $\psi_{0}^{*}$ as a ‘monotonisation’ of
$\psi_{0}$ (see the right panel of Figure 2) is consistent with that given in
Proposition 43 for a general $L^{2}(P)$-antitonic projection, where $P$ is a
univariate probability measure with a continuous distribution function.

When $i(p_{0})<\infty$, the inequality
$i^{*}(p_{0})=\int_{\mathbb{R}}(\psi_{0}^{*})^{2}\,dP_{0}\leq\int_{\mathbb{R}}\psi_{0}^{2}\,dP_{0}=i(p_{0})$
in Theorem 2 _(d)_ follows from the fact that the $L^{2}(P_{0})$-antitonic
projection onto the convex cone $\Psi_{\downarrow}(p_{0})$ is 1-Lipschitz with
respect to $\|{{\cdot}}\|_{L^{2}(P_{0})}$; see (115) in Lemma 44. A
statistical explanation of this information inequality arises from the fact
that the $1/i(p_{0})$ is the infimum of the asymptotic variance functional
$V_{p_{0}}(\psi)$ in (5) over all sufficiently regular
$\psi\colon\mathbb{R}\to\mathbb{R}$; see Huber and Ronchetti, (2009, Theorem
4.2), which we restate as Proposition 36. On the other hand, by (16),
$1/i^{*}(p_{0})$ is the minimum value of $V_{p_{0}}(\cdot)$ over the
restricted class $\Psi_{\downarrow}(p_{0})$, so in view of our discussion in
the introduction, it can be interpreted as an information lower bound for
convex $M$-estimators. When $i^{*}(p_{0})<\infty$, the ratio

$\mathrm{ARE}^{*}(p_{0}):=\frac{i^{*}(p_{0})}{i(p_{0})}$

therefore quantifies the price we pay in statistical efficiency for insisting
that our loss function be convex, and will be referred to as the antitonic
relative efficiency. By Theorem 2 _(d)_ , $\mathrm{ARE}^{*}(p_{0})\leq 1$ with
equality if and only if $p_{0}$ is log-concave, so we can regard
$1-\mathrm{ARE}^{*}(p_{0})$ as a measure of departure from log-concavity; see
Section 2.2 below. Example 11 shows that $\mathrm{ARE}^{*}(p_{0})\approx
0.878$ when $p_{0}$ is the Cauchy density, whereas Example 25 in the appendix
yields a density $p_{0}$ for which $\mathrm{ARE}^{*}(p_{0})=0$. More
generally, in Lemma 7 below, we provide a simple lower bound on
$\mathrm{ARE}^{*}(p_{0})$ that is reasonably tight for heavy-tailed densities
$p_{0}$.

When $p_{0}$ is only uniformly continuous (and not absolutely continuous), the
score function and Fisher information cannot be defined as above, but we
nevertheless refer to $\psi_{0}^{*}$ and $i^{*}(p_{0})$ as the _antitonic
projected score function_ (see Lemma 5 below) and _antitonic information (for
location)_ respectively.

######  Remark 3.

Since $F_{0}$ is continuous, we have $(F_{0}\circ F_{0}^{-1})(u)=u$ for all
$u\in(0,1)$. The concave function $\hat{J}_{0}$ is therefore absolutely
continuous on $[0,1]$ with derivative
$\hat{J}_{0}^{(\mathrm{R})}=\psi_{0}^{*}\circ F_{0}^{-1}$ Lebesgue almost
everywhere (Rockafellar,, 1997, Corollary 24.2.1), so

$i^{*}(p_{0})=\int_{\mathbb{R}}(\psi_{0}^{*})^{2}\,dP_{0}=\int_{0}^{1}(\psi_{0}^{*}\circ
F_{0}^{-1})^{2}=\int_{0}^{1}\bigl{(}\hat{J}_{0}^{(\mathrm{R})}\bigr{)}^{2}.$

In the setting of Theorem 2 _(d)_ above, a straightforward calculation (cf.
(52) in the proof of Theorem 2) shows that the density quantile function
$J_{0}=p_{0}\circ F_{0}^{-1}$ is absolutely continuous on $[0,1]$ with
derivative $J_{0}^{\prime}=\psi_{0}\circ F_{0}^{-1}$ Lebesgue almost
everywhere. Therefore, $\psi_{0}^{*}\circ
F_{0}^{-1}=\hat{J}_{0}^{(\mathrm{R})}=\widehat{\mathcal{M}}_{\mathrm{R}}(\psi_{0}\circ
F_{0}^{-1})$ almost everywhere and

$i(p_{0})=\int_{\mathbb{R}}\psi_{0}^{2}\,dP_{0}=\int_{0}^{1}(\psi_{0}\circ
F_{0}^{-1})^{2}=\int_{0}^{1}(J_{0}^{\prime})^{2}.$

Finally in this subsection, we define the _two-sided hazard function_
$h_{0}\colon\mathbb{R}\rightarrow[0,\infty)$ of $p_{0}$ by

$h_{0}(z):=\frac{p_{0}(z)}{F_{0}(z)\wedge\bigl{(}1-F_{0}(z)\bigr{)}}=\begin{cases}p_{0}(z)/F_{0}(z)&\text{if
}F_{0}(z)\leq 1/2\\\ p_{0}(z)/\bigl{(}1-F_{0}(z)\bigr{)}&\text{if
}F_{0}(z)>1/2,\end{cases}$ (19)

where in accordance with our convention $0/0=0$, we have $h_{0}(z)=0$ whenever
$F_{0}(z)\in\\{0,1\\}$. The following simple lemma provides a necessary and
sufficient condition on the two-sided hazard function for $\psi_{0}^{*}$ to be
appropriately bounded, which means that any negative antiderivative
$\ell_{0}^{*}$ grows at most linearly, and therefore ensures robustness.

######  Lemma 4.

In the setting of Lemma 1, define $z_{\min}:=\inf(\operatorname{supp}p_{0})$
and $z_{\max}:=\sup(\operatorname{supp}p_{0})$. Then

  1. (a)

$\lim_{z\to-\infty}\psi_{0}^{*}(z)<\infty$ if and only if $\limsup_{z\searrow
z_{\min}}h_{0}(z)<\infty$, in which case $z_{\min}=-\infty$;

  2. (b)

$\lim_{z\to\infty}\psi_{0}^{*}(z)>-\infty$ if and only if $\limsup_{z\nearrow
z_{\max}}h_{0}(z)<\infty$, in which case $z_{\max}=\infty$.

Recall that a Laplace density has a constant two-sided hazard function, as
well as a score function whose absolute value is constant. Roughly speaking,
the conditions on $h_{0}$ in Lemma 4 are satisfied by densities whose tails
are heavier than those of the Laplace density (Samworth and Johnson,, 2004),
for which it is particularly attractive to have bounded projected score
functions.

###  2.2 The log-concave Fisher divergence projection

Let $P_{0}$ and $P_{1}$ be Borel probability measures on $\mathbb{R}$ such
that $P_{0}\ll P_{1}$. Write $\operatorname{supp}P_{0}$ for the support of
$P_{0}$ (i.e. the smallest closed set $S$ satisfying $P_{0}(S)=1$), and
$\operatorname{Int}(\operatorname{supp}P_{0})$ for its interior. Suppose that
there exists a Radon–Nikodym derivative $dP_{0}/dP_{1}$ that is continuous on
$\operatorname{Int}(\operatorname{supp}P_{0})$, and also strictly positive and
differentiable on some subset
$E\subseteq\operatorname{Int}(\operatorname{supp}P_{0})$ such that
$P_{0}(E^{c})=0$.444This condition precludes $P_{0}$ from having any isolated
atoms, so in particular, $P_{0}$ cannot be a discrete measure. The Fisher
divergence (also known as the Fisher information distance555This is not to be
confused with the Fisher information (or Fisher–Rao) metric (Amari and
Nagaoka,, 2000, Chapter 2), a Riemannian metric on a manifold of probability
distributions.) from $P_{1}$ to $P_{0}$ is defined to be

$I(P_{0},P_{1}):=\int_{E}\biggl{(}\Bigl{(}\log\frac{dP_{0}}{dP_{1}}\Bigr{)}^{\prime}\biggr{)}^{2}\,dP_{0}.$
(20)

If $P_{0},P_{1}$ do not satisfy the assumptions above, then we define
$I(P_{0},P_{1}):=\infty$. In the case where $P_{0},P_{1}$ have Lebesgue
densities $p_{0},p_{1}$ respectively that are both locally absolutely
continuous on $\mathbb{R}$, we have

$I(p_{0},p_{1})\equiv
I(P_{0},P_{1})=\begin{cases}\displaystyle\int_{\\{p_{0}>0\\}}\biggl{(}\Bigl{(}\log\frac{p_{0}}{p_{1}}\Bigr{)}^{\prime}\biggr{)}^{2}\,p_{0}=\int_{\mathbb{R}}(\psi_{0}-\psi_{1})^{2}\,dP_{0}\;&\text{if
}\operatorname{supp}p_{0}\subseteq\operatorname{supp}p_{1}\\\
\infty\;&\text{otherwise},\end{cases}$

where we denote by $\psi_{j}:=(\log
p_{j})^{\prime}\mathbbm{1}_{\\{p_{j}>0\\}}=p_{j}^{\prime}/p_{j}$ the
corresponding score functions for $j\in\\{0,1\\}$. For further background on
the Fisher divergence, see Johnson, (2004, Definition 1.13), Yang et al.,
(2019, Section 2) and references therein.

The following lemma establishes the connection between the projected score
function and the Fisher divergence.

######  Lemma 5.

In the setting of Lemma 1, there is a unique continuous log-concave density
$p_{0}^{*}$ on $\mathbb{R}$ such that
$\operatorname{supp}p_{0}^{*}=\mathcal{S}_{0}$ and $\log p_{0}^{*}$ has right
derivative $\psi_{0}^{*}$ on $\mathcal{S}_{0}$. In particular,
$\psi_{0}^{*}=(\log p_{0}^{*})^{\prime}$ Lebesgue almost everywhere on
$\mathcal{S}_{0}$. Furthermore, if $p_{0}$ is itself log-concave, then
$p_{0}^{*}=p_{0}$.

When $p_{0}$ is absolutely continuous, (18) indicates that the log-concave
density $p_{0}^{*}$ (with score function
$\psi_{0}^{*}\mathbbm{1}_{\mathcal{S}_{0}}$) in Lemma 5 minimises $I(p_{0},p)$
over the class $\mathcal{P}_{\mathrm{LC}}$ of all univariate log-concave
densities $p$. Moreover, $p_{0}^{*}=p_{0}$ if and only if $p_{0}$ is log-
concave. Even when $p_{0}$ is only uniformly continuous, we refer to
$p_{0}^{*}$ as the log-concave Fisher divergence projection of $p_{0}$. In
contrast, the log-concave maximum likelihood projection $p_{0}^{\mathrm{ML}}$
of the distribution $P_{0}$ (Dümbgen et al.,, 2011; Barber and Samworth,,
2021) can be interpreted as a minimiser of Kullback–Leibler divergence rather
than Fisher divergence over the class of upper semi-continuous log-concave
densities. By Dümbgen et al., (2011, Theorem 2.2), $p_{0}^{\mathrm{ML}}$
exists and is unique if and only if $P_{0}$ is non-degenerate and has a finite
mean (but not necessarily a Lebesgue density). On the other hand, moment
conditions are not required for $p_{0}^{*}$ to exist and be unique, but
$p_{0}^{*}$ is only defined in Lemma 5 when $P_{0}$ has a uniformly continuous
density on $\mathbb{R}$. As we will discuss in Section 2.4, the non-existence
of the Fisher divergence projection for discrete measures $P_{0}$ has
consequences for our statistical methodology.

When $p_{0}$ is not log-concave, $p_{0}^{\mathrm{ML}}$ usually does not
coincide with $p_{0}^{*}$ even when both exist, and moreover the associated
regression $M$-estimators

$\hat{\beta}_{\psi_{0}^{\mathrm{ML}}}\in\operatorname*{argmax}_{\beta\in\mathbb{R}^{d}}\sum_{i=1}^{n}\log
p_{0}^{\mathrm{ML}}(Y_{i}-X_{i}^{\top}\beta)\quad\text{and}\quad\hat{\beta}_{\psi_{0}^{*}}\in\operatorname*{argmax}_{\beta\in\mathbb{R}^{d}}\sum_{i=1}^{n}\log
p_{0}^{*}(Y_{i}-X_{i}^{\top}\beta)$ (21)

are generally different; see Examples 12 and 26 below. In fact, the following
result shows that there exist error distributions $P_{0}$ for which the
asymptotic covariance of $\hat{\beta}_{\psi_{0}^{\mathrm{ML}}}$ is arbitrarily
large compared with that of the optimal convex $M$-estimator
$\hat{\beta}_{\psi_{0}^{*}}$, even when the latter is close to being
asymptotically efficient in the sense of (4).

######  Proposition 6.

For every $\epsilon\in(0,1)$, there exists a distribution $P_{0}$ with a
finite mean and an absolutely continuous density $p_{0}$ such that
$i(p_{0})<\infty$, and the log-concave maximum likelihood projection
$q_{0}\equiv p_{0}^{\mathrm{ML}}$ has corresponding score function
$\psi_{0}^{\mathrm{ML}}:=q_{0}^{(\mathrm{R})}/q_{0}\in\Psi_{\downarrow}(p_{0})$
satisfying

$\frac{V_{p_{0}}(\psi_{0}^{*})}{V_{p_{0}}(\psi_{0}^{\mathrm{ML}})}\leq\epsilon\quad\text{and}\quad\mathrm{ARE}^{*}(p_{0})\geq
1-\epsilon.$ (22)

Our next result provides a simple lower bound on the antitonic information.

######  Lemma 7.

Suppose that $p_{0}$ is an absolutely continuous density on $\mathbb{R}$ with
$i(p_{0})<\infty$. Then $p_{0}$ is bounded with $i^{*}(p_{0})\geq
4\|{p_{0}}\|_{\infty}^{2}$, so

$\mathrm{ARE}^{*}(p_{0})\geq\frac{4\|{p_{0}}\|_{\infty}^{2}}{i(p_{0})},$

with equality if and only if $p_{0}^{*}$ is a Laplace density, i.e. there
exist $\mu\in\mathbb{R}$ and $\sigma>0$ such that
$p_{0}^{*}(z)=(2\sigma)^{-1}\exp(-|z-\mu|/\sigma)$ for all $z\in\mathbb{R}$.

######  Remark 8.

A reassuring property of the antitonic projection is its affine equivariance:
if $p_{0}$ is a uniformly continuous density, then for $a>0$ and
$b\in\mathbb{R}$, the density $z\mapsto ap_{0}(az+b)=:p_{a,b}(z)$ has
antitonic projected score function and log-concave Fisher divergence
projection given by

$\psi_{a,b}^{*}(z):=a\psi_{0}^{*}(az+b)\quad\text{and}\quad
p_{a,b}^{*}(z):=ap_{0}^{*}(az+b)$

respectively for $z\in\mathbb{R}$. It follows that
$1/V_{p_{a,b}}(\psi_{a,b}^{*})=i^{*}(p_{a,b})=\int_{-\infty}^{\infty}(\psi_{a,b}^{*})^{2}\,p_{a,b}=a^{2}i^{*}(p_{0})$,
so because $\|{p_{a,b}}\|_{\infty}=a\|{p_{0}}\|_{\infty}$, both the antitonic
relative efficiency and the lower bound in Lemma 7 are affine invariant in the
sense that they remain unchanged if we replace $p_{0}$ with $p_{a,b}$.

Similarly, if $P_{0}$ has a finite mean, then by the affine equivariance of
the log-concave maximum likelihood projection (Dümbgen et al.,, 2011, Remark
2.4), $p_{a,b}^{\mathrm{ML}}(z)=ap_{0}^{\mathrm{ML}}(az+b)$ and hence
$\psi_{a,b}^{\mathrm{ML}}=a\psi_{0}^{\mathrm{ML}}(az+b)$ for $z\in\mathbb{R}$.
Thus,
$V_{p_{a,b}}(\psi_{a,b}^{\mathrm{ML}})=V_{p_{0}}(\psi_{0}^{\mathrm{ML}})/a^{2}$,
so the first ratio in (22) is also affine invariant. Consequently, for any
$C\in(0,\infty)$, $\mu\in\mathbb{R}$ and $\epsilon\in(0,1)$, there exists a
density $p_{0}$ satisfying (22) with $i(p_{0})=C$ and
$\int_{\mathbb{R}}zp_{0}(z)\,dz=\mu$.

The final result in this subsection relates properties of densities and their
log-concave Fisher divergence projections.

######  Proposition 9.

For a uniformly continuous density $p_{0}\colon\mathbb{R}\to\mathbb{R}$, the
log-concave Fisher divergence projection $p_{0}^{*}$ and its corresponding
distribution function $F_{0}^{*}\colon[-\infty,\infty]\to\mathbb{R}$ have the
following properties.

  1. (a)

Denote by $\mathcal{T}$ the set of $z\in\mathcal{S}_{0}$ such that
$\psi_{0}^{*}$ is non-constant on every open interval containing $z$. For
$z\in\mathcal{T}$, we have

$\frac{p_{0}^{*}(z)}{F_{0}^{*}(z)}\leq\frac{p_{0}(z)}{F_{0}(z)}\quad\text{and}\quad\frac{p_{0}^{*}(z)}{1-F_{0}^{*}(z)}\leq\frac{p_{0}(z)}{1-F_{0}(z)},$
(23)

whence $p_{0}^{*}(z)\leq p_{0}(z)$.

  2. (b)

$\|{p_{0}^{*}}\|_{\infty}\leq\|{p_{0}}\|_{\infty}$ and
$i(p_{0}^{*})=-\int_{\mathbb{R}}p_{0}^{*}\,d\psi_{0}^{*}\leq-\int_{\mathbb{R}}p_{0}\,d\psi_{0}^{*}=i^{*}(p_{0})$.

We can define the two-sided hazard function $h_{0}^{*}$ of the log-concave
Fisher divergence projection $p_{0}^{*}$ analogously to $h_{0}$ in (19). Since
$p_{0}^{*}$ is log-concave, its density quantile function
$J_{0}^{*}:=p_{0}^{*}\circ(F_{0}^{*})^{-1}$ has decreasing right derivative
$(J_{0}^{*})^{(\mathrm{R})}=\psi_{0}^{*}\circ(F_{0}^{*})^{-1}$ by (the proof
of) Lemma 5, so $J_{0}^{*}$ is concave on $[0,1]$. Thus,

$0=J_{0}^{*}(0)\leq
J_{0}^{*}(u)-u(J_{0}^{*})^{(\mathrm{R})}(u)\quad\text{and}\quad
0=J_{0}^{*}(1)\leq J_{0}^{*}(u)+(1-u)(J_{0}^{*})^{(\mathrm{R})}(u)$

for all $u\in[0,1]$, so for $z\in\mathcal{T}$, Lemma 1 and (23) in Proposition
9 _(a)_ imply that

$\displaystyle|\psi_{0}^{*}(z)|=\bigl{|}(J_{0}^{*})^{(\mathrm{R})}\bigl{(}F_{0}^{*}(z)\bigr{)}\bigr{|}\leq\frac{J_{0}^{*}\bigl{(}F_{0}^{*}(z)\bigr{)}}{F_{0}^{*}(z)\wedge\bigl{(}1-F_{0}^{*}(z)\bigr{)}}$
$\displaystyle=\frac{p_{0}^{*}(z)}{F_{0}^{*}(z)\wedge\bigl{(}1-F_{0}^{*}(z)\bigr{)}}=h_{0}^{*}(z)$
$\displaystyle\leq\frac{p_{0}(z)}{F_{0}(z)\wedge\bigl{(}1-F_{0}(z)\bigr{)}}=h_{0}(z).$

Proposition 9 _(b)_ provides inequalities on the supremum norm and antitonic
information of the log-concave Fisher divergence projection. In particular,
the Fisher information of the projected density is at most the antitonic
information of the original density.

###  2.3 Examples

######  Example 10.

Let $p_{0}$ be the $\mathrm{Beta}(a,b)$ density given by
$p_{0}(z)=z^{a-1}(1-z)^{b-1}\mathbbm{1}_{\\{z\in(0,1)\\}}/\mathrm{B}(a,b)$ for
$a,b>1$, where $\mathrm{B}$ denotes the beta function. Then $p_{0}$ is
uniformly continuous and log-concave on $\mathbb{R}$, so

$\psi_{0}^{*}(z)=\psi_{0}(z)=(\log
p_{0})^{\prime}(z)=\frac{a-1}{z}-\frac{b-1}{1-z}$

for all $z\in(0,1)$, while $\psi_{0}^{*}(z)=\infty$ for $z\leq 0$ and
$\psi_{0}^{*}(z)=-\infty$ for $z\geq 1$. We have
$i^{*}(p_{0})=i(p_{0})=\int_{0}^{1}\psi_{0}^{2}\,p_{0}<\infty$ if and only if
$a,b>2$, in which case the conclusions of Theorem 2(c, d) hold.

Figure 3: _Left_ : The negative log-density $\ell_{0}=-\log p_{0}$ and the
optimal convex loss function $\ell_{0}^{*}=-\log p_{0}^{*}$ when $p_{0}$ is
the standard Cauchy density. _Right_ : The corresponding densities $p_{0}$ and
$p_{0}^{*}$.

######  Example 11.

Let $p_{0}$ be the standard Cauchy density given by
$p_{0}(z)=1/\bigl{(}\pi(1+z^{2})\bigr{)}$ for $z\in\mathbb{R}$. Then $p_{0}$
is absolutely continuous on $\mathbb{R}$ with
$i(p_{0})=\int_{\mathbb{R}}(p_{0}^{\prime})^{2}/p_{0}=1/2$ and
$\lim_{z\to\pm\infty}h_{0}(z)=0$. We will derive an explicit expression for
$\psi_{0}^{*}$, which is necessarily bounded by Lemma 4. In contrast to the
previous example, $p_{0}$ is not log-concave, so $\psi_{0}^{*}$ does not
coincide with $\psi_{0}=p_{0}^{\prime}/p_{0}\colon z\mapsto-2z/(1+z^{2})$.
Indeed,

$\displaystyle F_{0}(z)$
$\displaystyle=\frac{1}{2}+\frac{\arctan(z)}{\pi}\;\;\text{for
}z\in\mathbb{R},\qquad
F_{0}^{-1}(u)=\tan\biggl{(}\pi\Bigl{(}u-\frac{1}{2}\Bigr{)}\biggr{)}=-\cot(\pi
u)\;\;\text{for }u\in(0,1),$ $\displaystyle J_{0}(u)$
$\displaystyle=\frac{1}{\pi\bigl{(}1+\cot^{2}(\pi
u)\bigr{)}}=\frac{\sin^{2}(\pi u)}{\pi}=\frac{1-\cos(2\pi
u)}{2\pi}\;\;\text{for }u\in[0,1].$

Let $\zeta_{0}\approx 2.33$ be the unique $\zeta\in(0,\pi)$ satisfying
$\zeta=\tan(\zeta/2)$, and define $u_{0}:=\zeta_{0}/(2\pi)\in(0,1/2)$. Then we
can verify that $\hat{J}_{0}$ is linear on $[0,u_{0}]$ and on $[1-u_{0},1]$,
with $\hat{J}_{0}(u)=J_{0}(u)$ for $u\in[u_{0},1-u_{0}]\cup\\{0,1\\}$; see the
left panel of Figure 2. It follows that

$\displaystyle\hat{J}_{0}^{(\mathrm{R})}(u)$
$\displaystyle=\begin{cases}\,\sin\zeta_{0}&\text{for }u\in[0,u_{0}]\\\
\,\sin(2\pi u)&\text{for }u\in[u_{0},1-u_{0}]\\\ \,-\sin\zeta_{0}&\text{for
}u\in[1-u_{0},1];\end{cases}$
$\displaystyle\psi_{0}^{*}(z)=\hat{J}_{0}^{(\mathrm{R})}\bigl{(}F_{0}(z)\bigr{)}$
$\displaystyle=\begin{cases}\,\sin\zeta_{0}=-2z_{0}/(1+z_{0}^{2})&\text{for
}z\in(-\infty,-z_{0}]\\\ \,-\sin(2\arctan
z)=-2z/(1+z^{2})=\psi_{0}(z)&\text{for }z\in[-z_{0},z_{0}]\\\
\,-\sin\zeta_{0}&\text{for }z\in[z_{0},\infty),\end{cases}$

where $z_{0}:=\cot(\pi u_{0})=\cot(\zeta_{0}/2)\approx 0.43\in(0,1)$ satisfies
$z_{0}\arctan(1/z_{0})=1/2$. Thus, $\psi_{0}^{*}(z)=\psi_{0}\bigl{(}(z\wedge
z_{0})\vee(-z_{0})\bigr{)}$ for $z\in\mathbb{R}$; see the right panel of
Figure 2. An antiderivative $\phi_{0}^{*}$ of $\psi_{0}^{*}$ is given by

$\phi_{0}^{*}(z):=\int_{0}^{z}\psi_{0}^{*}-\log\pi=\begin{cases}\,-\log\bigl{(}\pi(1+z^{2})\bigr{)}&\text{for
}z\in[-z_{0},z_{0}]\\\
\,-(|z|-z_{0})\sin\zeta_{0}-\log\bigl{(}\pi(1+z_{0}^{2})\bigr{)}&\text{for
}z\in\mathbb{R}\setminus[-z_{0},z_{0}].\end{cases}$

As illustrated in the left panel of Figure 3, $\ell_{0}^{*}:=-\phi_{0}^{*}$ is
a symmetric convex function that is approximately quadratic on
$[-z_{0},z_{0}]$ and linear outside this interval, so in this respect, it
resembles the Huber loss function (7). This is significant as far as
$M$-estimation is concerned, since Huber-like loss functions are designed
precisely to be robust to outliers, such as those that arise in regression
problems with heavy-tailed Cauchy errors. As discussed in the introduction,
$\ell_{0}^{*}$ is optimal in the sense that the resulting regression
$M$-estimator
$\hat{\beta}_{\psi_{0}^{*}}\in\operatorname*{argmin}_{\beta\in\mathbb{R}^{d}}\sum_{i=1}^{n}\ell_{0}^{*}(Y_{i}-X_{i}^{\top}\beta)$
has minimal asymptotic covariance (3) among all convex $M$-estimators. By
direct computation,

$\frac{1}{V_{p_{0}}(\psi_{0}^{*})}=i^{*}(p_{0})=\frac{1}{2}-\frac{2\zeta_{0}\cos(2\zeta_{0})-\sin(2\zeta_{0})}{4\pi}\approx
0.439,\;\;\text{so}\;\;\mathrm{ARE}^{*}(p_{0})=\frac{i^{*}(p_{0})}{i(p_{0})}\approx
0.878$

in this case, meaning that the restriction to convex loss functions results in
only a small loss of efficiency relative to the maximum likelihood estimator.
This may well be outweighed by the increased computational convenience of
optimising a convex empirical risk function as opposed to a Cauchy likelihood
function, which typically has several local extrema; see van der Vaart, (1998,
Example 5.50) for a discussion of the difficulties involved. Lemma 7 yields
the bound $\mathrm{ARE}^{*}(p_{0})\geq
8\|{p_{0}}\|_{\infty}^{2}=8/\pi^{2}\approx 0.811$.

Figure 4: Plot of the asymptotic relative efficiency $r(K)$ of the Huber
$M$-estimator $\hat{\beta}_{\psi_{K}}$ compared with the optimal convex
$M$-estimator.

For $K>0$, the Huber regression $M$-estimator
$\hat{\beta}_{\psi_{K}}\in\operatorname*{argmin}_{\beta\in\mathbb{R}^{d}}\sum_{i=1}^{n}\ell_{K}(Y_{i}-X_{i}^{\top}\beta)$
defined with respect to (7) has asymptotic relative efficiency

$r(K):=\frac{V_{p_{0}}(\psi_{0}^{*})}{V_{p_{0}}(\psi_{K})}=\frac{\pi\bigl{(}\pi
K^{2}+2K-2(1+K^{2})\arctan K\bigr{)}}{4i^{*}(p_{0})\arctan^{2}K}$

compared with the optimal convex $M$-estimator $\hat{\beta}_{\psi_{0}^{*}}$;
see Figure 4. The maximum value $\sup_{K>0}r(K)\approx 0.9998$ is attained at
$K^{*}\approx 0.394$. Moreover,

$\lim_{K\to 0}r(K)=\frac{4}{\pi^{2}i^{*}(p_{0})}\approx 0.922$

is the asymptotic relative efficiency
$V_{p_{0}}(\psi_{0}^{*})/V_{p_{0}}(\psi)$ of the least absolute deviation
(LAD) estimator
$\hat{\beta}_{\psi}\in\operatorname*{argmin}_{\beta\in\mathbb{R}^{d}}\sum_{i=1}^{n}|Y_{i}-X_{i}^{\top}\beta|$,
for which $\psi(\cdot):=-\operatorname{sgn}(\cdot)$ and
$V_{p_{0}}(\psi)=1/\bigl{(}4p_{0}(0)^{2}\bigr{)}=\pi^{2}/4$. On the other
hand, the Huber loss $\ell_{K}$ in (7) converges pointwise to the squared
error loss as $K\to\infty$, so
$\lim_{K\to\infty}V_{p_{0}}(\psi_{K})=\int_{\mathbb{R}}z^{2}\,p_{0}(z)\,dz=\infty$
and hence $\lim_{K\to\infty}r(K)=0$. We recall the difficulties of choosing
$K$, and its connection to the choice of scale, from the discussion in the
introduction.

The Cauchy density $p_{0}$ and its log-concave Fisher divergence projection
$p_{0}^{*}:=e^{\phi_{0}^{*}}$ are plotted in the right panel of Figure 3.
Since $p_{0}=p_{0}^{*}$ on $[-z_{0},z_{0}]$ and $\psi_{0}^{*}$ is constant on
$\mathbb{R}\setminus[-z_{0},z_{0}]$, it turns out that
$i^{*}(p_{0})=-\int_{\mathbb{R}}p_{0}\,d\psi_{0}^{*}=-\int_{\mathbb{R}}p_{0}^{*}\,d\psi_{0}^{*}=i(p_{0}^{*})$,
so both inequalities in Proposition 9(b) are in fact equalities in this
example.

In Example 26 in the appendix, we take $P_{0}$ to be a scaled $t_{2}$
distribution, which has a finite first moment (unlike the Cauchy distribution
in Example 11), and verify that $\hat{\beta}_{\psi_{0}^{\mathrm{ML}}}$ and
$\hat{\beta}_{\psi_{0}^{*}}$ in (21) are different convex $M$-estimators.
Example 27 features a symmetrised Pareto density with polynomially decaying
tails, where the optimal convex loss function $\ell_{0}^{*}$ is a scale
transformation of the robust absolute error loss $z\mapsto|z|$.

Moving on from heavy-tailed distributions, we now consider a density $p_{0}$
that fails to be log-concave because it is not unimodal.

######  Example 12.

For $\rho\in(0,1)$ and $\mu>0$, let $p_{0}$ be the two-component Laplace
mixture density given by

$p_{0}(z)=\frac{1-\rho}{2}e^{-|z+\mu|}+\frac{\rho}{2}e^{-|z-\mu|}$

for $z\in\mathbb{R}$. The corresponding score function $\psi_{0}=(\log
p_{0})^{\prime}=p_{0}^{\prime}/p_{0}$ satisfies

$\psi_{0}(z)=\begin{cases}1\;\;&\text{for }z<-\mu\\\ \dfrac{\rho
e^{z}-(1-\rho)e^{-z}}{\rho e^{z}+(1-\rho)e^{-z}}\;\;&\text{for
}z\in(-\mu,\mu)\\\ -1\;\;&\text{for }z>\mu,\end{cases}$

so $\psi_{0}$ is strictly increasing on $(-\mu,\mu)$ and constant on both
$(-\infty,\mu)$ and $(\mu,\infty)$. Since $J_{0}^{\prime}=\psi_{0}\circ
F_{0}^{-1}$ on $(0,1)\setminus\\{F_{0}(-\mu),F_{0}(\mu)\\}$, it follows that
$J_{0}$ is convex on $[F_{0}(-\mu),F_{0}(\mu)]$ while being linear on both
$[0,F_{0}(-\mu)]$ and $[F_{0}(\mu),1]$. Therefore, $\hat{J}_{0}=J_{0}$ on
$[0,F_{0}(-\mu)]\cup[F_{0}(\mu),1]$ and $\hat{J}_{0}$ is linear on
$[F_{0}(-\mu),F_{0}(\mu)]$ with
$J_{0}\bigl{(}F_{0}(\pm\mu)\bigr{)}=p_{0}(\pm\mu)$, so

$\psi_{0}^{*}(z)=\hat{J}_{0}^{(\mathrm{R})}\bigl{(}F_{0}(z)\bigr{)}=\begin{cases}1\;\;&\text{for
}z<-\mu\\\
\dfrac{p_{0}(\mu)-p_{0}(-\mu)}{F_{0}(\mu)-F_{0}(-\mu)}=2\rho-1\;\;&\text{for
}z\in[-\mu,\mu)\\\ -1\;\;&\text{for }z\geq\mu,\end{cases}$

and hence

$p_{0}^{*}(z)=\begin{cases}Ce^{z+\mu}\;\;&\text{for }z\leq-\mu\\\
Ce^{(2\rho-1)(z+\mu)}\;\;&\text{for }z\in[-\mu,\mu]\\\ Ce^{(2\rho-1)2\mu}\cdot
e^{-(z-\mu)}\;\;&\text{for }z\geq\mu,\end{cases}\qquad\text{where
}C:=\frac{1}{1+e^{(2\rho-1)2\mu}+\int_{0}^{2\mu}e^{(2\rho-1)z}\,dz}.$

By direct calculation, $p_{0}^{*}(-\mu)=C<p_{0}(-\mu)$ and
$p_{0}^{*}(\mu)=Ce^{(2\rho-1)2\mu}<p_{0}(\mu)$, so
$\|{p_{0}^{*}}\|_{\infty}=p_{0}^{*}(-\mu)\vee p_{0}^{*}(\mu)<p_{0}(-\mu)\vee
p_{0}(\mu)=\|{p_{0}}\|_{\infty}$ and

$i^{*}(p_{0})=-\int_{\mathbb{R}}p_{0}\,d\psi_{0}^{*}=2(1-\rho)\cdot
p_{0}(-\mu)+2(1+\rho)\cdot p_{0}(\mu)>2(1-\rho)\cdot
p_{0}^{*}(-\mu)+2(1+\rho)\cdot p_{0}^{*}(\mu)=i(p_{0}^{*}).$

We mention that Proposition 6 can be proved by slightly modifying the
construction of $p_{0}$ in Example 12.

###  2.4 Estimation of the projected score function

Theorem 2 motivates a general-purpose nonparametric procedure to estimate the
antitonic projected score function $\psi_{0}^{*}$ based on a sample
$\varepsilon_{1},\dotsc,\varepsilon_{n}\overset{\mathrm{iid}}{\sim}p_{0}$. We
will first explain why a naive approach fails before describing our
methodology. For a locally absolutely continuous function
$\psi\colon\mathbb{R}\to\mathbb{R}$ with derivative $\psi^{\prime}$, the
empirical analogue of
$D_{p_{0}}(\psi)=\mathbb{E}\bigl{(}\psi^{2}(\varepsilon_{1})+2\psi^{\prime}(\varepsilon_{1})\bigr{)}$
in (14) is

$\hat{D}_{n}(\psi)\equiv\hat{D}_{n}(\psi;\varepsilon_{1},\dotsc,\varepsilon_{n}):=\frac{1}{n}\sum_{i=1}^{n}\\{\psi^{2}(\varepsilon_{i})+2\psi^{\prime}(\varepsilon_{i})\\}.$
(24)

Recall from the introduction that score matching estimates the score function
$\psi_{0}=p_{0}^{\prime}/p_{0}$ (associated with a locally absolutely
continuous $p_{0}$) by an empirical risk minimiser
$\hat{\psi}_{n}\in\operatorname*{argmin}_{\psi\in\Psi}\hat{D}_{n}(\psi)$ over
an appropriate class of functions $\Psi$.

However, to obtain a monotone score estimate, we cannot minimise
$\psi\mapsto\hat{D}_{n}(\psi)$ directly over the class
$\Psi_{\downarrow}^{\mathrm{ac}}$ of all decreasing, locally absolutely
continuous $\psi\colon\mathbb{R}\to\mathbb{R}$. Indeed,
$\inf_{\psi\in\Psi_{\downarrow}^{\mathrm{ac}}}\hat{D}_{n}(\psi)=-\infty$, as
can be seen by constructing differentiable approximations to a decreasing step
function whose jumps are at the data points
$\varepsilon_{1},\dotsc,\varepsilon_{n}$. To circumvent this issue, we instead
propose the following estimation strategy.

Antitonic projected score estimation: Consider smoothing the empirical
distribution of $\varepsilon_{1},\dotsc,\varepsilon_{n}$, for example by
convolving it with an absolutely continuous kernel
$K\colon\mathbb{R}\to\mathbb{R}$ to obtain a kernel density estimator
$z\mapsto\tilde{p}_{n}(z):=n^{-1}\sum_{i=1}^{n}K_{h}(z-\varepsilon_{i})$,
where $h>0$ is a suitable bandwidth and $K_{h}(\cdot):=h^{-1}K(\cdot/h)$. We
can then define the smoothed empirical score matching objective

$\tilde{D}_{n}(\psi):=D_{\tilde{p}_{n}}(\psi)=\int_{-\infty}^{\infty}\psi^{2}\,\tilde{p}_{n}+2\int_{\mathcal{S}_{0}}\tilde{p}_{n}\,d\psi$

for $\psi\in\Psi_{\downarrow}(\tilde{p}_{n})$, which approximates the
population expectation in the definition of $D_{p_{0}}(\psi)$. Then by Theorem
2,

$\hat{\psi}_{n}:=\hat{J}_{n}^{(\mathrm{R})}\circ\tilde{F}_{n}\in\operatorname*{argmin}_{\psi\in\Psi_{\downarrow}(\tilde{p}_{n})}\tilde{D}_{n}(\psi),$
(25)

where $\tilde{F}_{n}$ denotes the distribution function corresponding to
$\tilde{p}_{n}$, and $J_{n}:=\tilde{p}_{n}\circ\tilde{F}_{n}^{-1}$. This is
reminiscent of maximum smoothed likelihood estimation of a density or
distribution function (e.g. Eggermont and LaRiccia,, 2000; Groeneboom and
Jongbloed,, 2014, Sections 8.2 and 8.5). By evaluating the antiderivative of
$J_{n}$ on a suitably fine grid, we may obtain a piecewise affine
approximation to $\hat{J}_{n}$ (and hence a piecewise constant approximation
to its derivative) using PAVA, whose space and time complexities scale
linearly with the size of the grid (Samworth and Shah,, 2024, Section 10.3.1).

Figure 5: Kernel-based estimates of the projected score function and optimal
convex loss function based on a sample of size $n=2000$ from the Cauchy
distribution.

More generally, we can use $\varepsilon_{1},\dotsc,\varepsilon_{n}$ to
construct a generic (not necessarily monotone) score estimator
$\tilde{\psi}_{n}$ and an estimate $\hat{F}_{n}$ of the distribution function
$F_{0}$ corresponding to the density $p_{0}$. By analogy with the explicit
representation (17) of $\psi_{0}^{*}$, we then define the decreasing score
estimate

$\hat{\psi}_{n}:=\widehat{\mathcal{M}}_{\mathrm{R}}(\tilde{\psi}_{n}\circ\hat{F}_{n}^{-1})\circ\hat{F}_{n}.$
(26)

As explained above, (an approximation to) $\hat{\psi}_{n}$ can be computed
efficiently using isotonic regression algorithms. In particular, if
$\hat{F}_{n}$ is taken to be the empirical distribution function of
$\varepsilon_{1},\dotsc,\varepsilon_{n}$, then by Proposition 43,
$\hat{\psi}_{n}^{\mathrm{L}}:=\widehat{\mathcal{M}}_{\mathrm{L}}(\tilde{\psi}_{n}\circ\hat{F}_{n}^{-1})\circ\hat{F}_{n}$
is an antitonic least squares estimator based on
$\bigl{\\{}\bigl{(}\varepsilon_{i},\tilde{\psi}_{n}(\varepsilon_{i})\bigr{)}:i\in[n]\bigr{\\}}$.
Our decreasing score estimate can be taken to be either
$\hat{\psi}_{n}^{\mathrm{L}}$ or the closely related $\hat{\psi}_{n}$: for
every $z\in\mathbb{R}$, we have
$\hat{\psi}_{n}(z)=\hat{\psi}_{n}^{\mathrm{L}}\bigl{(}\varepsilon(z)\bigr{)}$,
where $\varepsilon(z)$ is the smallest element of
$\\{\varepsilon_{1},\dotsc,\varepsilon_{n}\\}$ that is either strictly greater
than $z$ or equal to $\max_{i\in[n]}\varepsilon_{i}$.

Figure 6: Kernel-based estimates of the projected score function and optimal
convex loss function based on a sample of size $n=10^{4}$ from the Gaussian
mixture distribution $0.4N(-2,1)+0.6N(2,1)$.

The transformation (26) may be applied to any appropriate initial score
estimator $\tilde{\psi}_{n}$. Since a misspecified parametric method may
introduce significant error at the outset, we seek a nonparametric estimator.
For instance, we may take $\tilde{\psi}_{n}$ to be a ratio of kernel density
estimates of $p_{0}^{\prime}$ and $p_{0}$, which may be truncated for
theoretical and practical convenience to avoid instability in low-density
regions; see (31) in Section 3.1. Observe that (25) is a special case of (26)
with $\tilde{\psi}_{n}=\tilde{p}_{n}^{\prime}/\tilde{p}_{n}$ and
$\hat{F}_{n}=\tilde{F}_{n}$.

##  3 Semiparametric $M$-estimation via antitonic score matching

Suppose that we observe independent and identically distributed pairs
$(X_{1},Y_{1}),\dotsc,(X_{n},Y_{n})$ satisfying

$Y_{i}=X_{i}^{\top}\beta_{0}+\varepsilon_{i}\quad\text{for }i\in[n],$ (27)

where $X_{1},\dotsc,X_{n}$ are $\mathbb{R}^{d}$-valued covariates that are
independent of errors $\varepsilon_{1},\dotsc,\varepsilon_{n}$ with an unknown
absolutely continuous (Lebesgue) density $p_{0}$ on $\mathbb{R}$. We do not
necessarily insist that $\mathbb{E}(\varepsilon_{1})=0$, because we also have
in mind settings where we might want to assume that a quantile (e.g. the
median) of $\varepsilon_{1}$ is zero. In fact, we do not even assume that
$\mathbb{E}(\varepsilon_{1})$ is well-defined, which allows us to consider
Cauchy errors, for instance, as a running example. A price to pay for this
generality is that if an intercept term is present in (27), i.e. $X_{id}=1$
for each $i\in[n]$, then the corresponding component $\beta_{0d}$ of
$\beta_{0}$ is unidentifiable. On a related point, a necessary condition for
$\beta_{0}$ to be identifiable is that $\mathbb{E}(X_{1}X_{1}^{\top})$ is
positive definite. Indeed, if this matrix is singular, then there exists
$v\in\mathbb{R}^{d}\setminus\\{0\\}$ such that $X_{i}^{\top}v=0$ for all $i$
almost surely; in that case, the joint distribution of our observed data is
unchanged if we replace $\beta_{0}$ with $\beta_{0}+v$. This identifiability
issue is discussed in greater detail following the statement of Proposition
13.

For $\beta\in\mathbb{R}^{d}$, define $q_{\beta}\colon\mathbb{R}\to\mathbb{R}$
by

$q_{\beta}(z):=\mathbb{E}p_{0}\bigl{(}z-X_{1}^{\top}(\beta_{0}-\beta)\bigr{)},$

so that $q_{\beta}$ is the density of
$Y_{1}-X_{1}^{\top}\beta=\varepsilon_{1}+X_{1}^{\top}(\beta_{0}-\beta)$.

Our approach to estimating $\beta_{0}$ with data-driven convex loss functions
is motivated by the following population-level optimisation problem. We seek
to minimise the augmented score matching objective

$\displaystyle
Q(\beta,\psi):=D_{q_{\beta}}(\psi)=\mathbb{E}\biggl{(}\int_{-\infty}^{\infty}\psi^{2}\bigl{(}z+X_{1}^{\top}(\beta_{0}-\beta)\bigr{)}\,p_{0}(z)\,dz+2\int_{\mathbb{R}}p_{0}\bigl{(}z-X_{1}^{\top}(\beta_{0}-\beta)\bigr{)}\,d\psi(z)\biggr{)}$
(28)

jointly over $\beta\in\mathbb{R}^{d}$ and
$\psi\in\Psi_{\downarrow}(q_{\beta})$ satisfying

$\mathbb{E}\bigl{(}X_{1}\psi(Y_{1}-X_{1}^{\top}\beta)\bigr{)}=0.$ (29)

Recall from (14) that if $\psi\in\Psi_{\downarrow}(q_{\beta})$ is locally
absolutely continuous on $\mathbb{R}$, then

$Q(\beta,\psi)=\mathbb{E}\psi^{2}(Y_{1}-X_{1}^{\top}\beta)+2\mathbb{E}\psi^{\prime}(Y_{1}-X_{1}^{\top}\beta).$

The constraint (29) is the population analogue of the estimating equations (2)
based on $\psi$. It forces $\beta$ to be a minimiser of the convex population
risk function $b\mapsto\mathbb{E}\ell(Y_{1}-X_{1}^{\top}b)=:L_{\psi}(b)$ over
$\mathcal{D}_{\psi}$, where $\ell$ is any negative antiderivative of $\psi$,
and $\mathcal{D}_{\psi}$ is the convex set of $b\in\mathbb{R}^{d}$ such that
$L_{\psi}(b)$ is well-defined in $[-\infty,\infty]$.666Indeed, $\ell$ is
convex with $\psi=-\ell^{(\mathrm{R})}$, so for
$\beta,\beta^{\prime}\in\mathbb{R}^{d}$, we have
$\ell(Y_{1}-X_{1}^{\top}\beta^{\prime})\geq\ell(Y_{1}-X_{1}^{\top}\beta)+(\beta-\beta^{\prime})^{\top}X_{1}\psi(Y_{1}-X_{1}^{\top}\beta).$
If $\beta$ satisfies (29), then taking expectations in the display above shows
that $L_{\psi}(\beta^{\prime})\geq L_{\psi}(\beta)$ for all
$\beta^{\prime}\in\mathcal{D}_{\psi}$..

The following proposition characterises the global joint minimiser of our
constrained score matching optimisation problem. For $j\in[d]$, write
$\mathsf{e}_{j}$ for the $j$th standard basis vector in $\mathbb{R}^{d}$.

######  Proposition 13.

For an absolutely continuous density $p_{0}$ with
$\mathcal{S}(p_{0})=\mathbb{R}$, let $\psi^{*}_{0}$ be the projected score
function defined in (17). Assume that
$i^{*}(p_{0})=\int_{\mathbb{R}}(\psi_{0}^{*})^{2}\,p_{0}<\infty$ and that
$\mathbb{E}(X_{1}X_{1}^{\top})\in\mathbb{R}^{d\times d}$ is positive definite.
Let

$\Gamma:=\bigl{\\{}(\beta,\psi):\beta\in\mathbb{R}^{d},\,\psi\in\Psi_{\downarrow}(q_{\beta}),\,\mathbb{E}\bigl{(}X_{1}\psi(Y_{1}-X_{1}^{\top}\beta)\bigr{)}=0\bigr{\\}}.$

  1. (a)

Suppose that $p_{0}$ is symmetric. Denote by
$\Psi_{\downarrow}^{\mathrm{anti}}(q_{\beta})$ the set of all
$\psi\in\Psi_{\downarrow}(q_{\beta})$ such that
$\psi(-z)=-\lim_{z^{\prime}\nearrow z}\psi(z^{\prime})$, and let
$\Gamma^{\mathrm{anti}}:=\\{(\beta,\psi)\in\Gamma:\psi\in\Psi_{\downarrow}^{\mathrm{anti}}(q_{\beta})\\}$.
Then

$\displaystyle(\beta_{0},\psi_{0}^{*})=\operatorname*{argmin}_{(\beta,\psi)\in\Gamma^{\mathrm{anti}}}Q(\beta,\psi).$

  2. (b)

If instead $X_{1d}=1$ almost surely, then

$\displaystyle\bigl{(}\beta_{0}+c\mathsf{e}_{d},\psi_{0}^{*}(\cdot+c)\bigr{)}\in\operatorname*{argmin}_{(\beta,\psi)\in\Gamma}\,Q(\beta,\psi)$

for all $c\in\mathbb{R}$. Moreover, all minimisers are of this form.

Thus, when we restrict attention to (right-continuous versions of)
antisymmetric functions $\psi$ in the setting of Proposition 13 _(a)_ , the
unique minimiser of our population-level optimisation problem is given by the
pair consisting of the vector $\beta_{0}$ of true regression coefficients, and
the antitonic score projection $\psi_{0}^{*}$ of the error density $p_{0}$.
This motivates our statistical methodology below, where we seek to minimise
sample versions of this population-level objective. From the proof of
Proposition 13 _(a)_ , we see that for any
$(\beta,\psi)\in\Gamma^{\mathrm{anti}}$, we also have
$(\beta_{0},\psi)\in\Gamma^{\mathrm{anti}}$; moreover,
$Q(\beta_{0},\psi_{0}^{*})=\inf_{\psi\in\Psi_{\downarrow}^{\mathrm{anti}}(q_{\beta_{0}})}Q(\beta_{0},\psi)=-1/i^{*}(p_{0})<0$.
On the other hand, if $\beta\neq\beta_{0}$, and if
$\psi\in\Psi_{\downarrow}^{\mathrm{anti}}(q_{\beta})$ is such that
$(\beta,\psi)\in\Gamma^{\mathrm{anti}}$, then $Q(\beta,\psi)\geq 0$. In other
words, feasible pairs $(\beta,\psi)$ with $\beta\neq\beta_{0}$ are well-
separated from the optimal solution in terms of their objective function
values.

In Proposition 13 _(b)_ , where the errors need not have mean zero and where
we no longer restrict attention to antisymmetric $\psi$, the pair
$(\beta_{0},\psi_{0}^{*})$ still minimises our objective up to appropriate
translations to account for the lack of identifiability of the intercept term.
Indeed, the joint distribution of $(X_{1},Y_{1})$ does not change if we
replace $\beta_{0}$ and $\varepsilon_{1}$ with $\beta_{0}+c\mathsf{e}_{d}$ and
$\varepsilon_{1}-c$ respectively, for any $c\in\mathbb{R}$. This explains why
the minimiser in _(b)_ is not unique; observe that $\psi_{0}^{*}(\cdot+c)$ is
the projected score corresponding to the density $p_{0}(\cdot+c)$ of
$\varepsilon_{1}-c$. Nevertheless, _(b)_ indicates that translations of the
intercept term are the only source of non-uniqueness, so that the first $d-1$
components of $\beta_{0}$ can still be recovered by solving the constrained
optimisation problem above on the population level.

Since the objective function $(\beta,\psi)\mapsto Q(\beta,\psi)$ is not
jointly convex, we could seek to minimise $Q(\beta,\psi)$ subject to the
constraint (29) by alternating the following two steps:

  1. I.

For a fixed $\beta$, minimise the (convex) score matching objective
$\psi\mapsto Q(\beta,\psi)=D_{q_{\beta}}(\psi)$ based on the density
$q_{\beta}$ of $Y_{1}-X_{1}^{\top}\beta$, without imposing (29). This
constraint must be omitted here as otherwise $\beta$ will not be updated in
Step II below.

  2. II.

For a fixed decreasing and right-continuous $\psi$, minimise the convex
function $\beta\mapsto\mathbb{E}\ell(Y_{1}-X_{1}^{\top}\beta)$, where $\ell$
is a negative antiderivative of $\psi$.

In an empirical version of this alternating minimisation algorithm based on
$(X_{1},Y_{1}),\dotsc,(X_{n},Y_{n})$, we can estimate $\psi_{0}^{*}$ in Step I
by minimising a sample analogue of $D_{q_{\beta}}(\psi)$ over $\psi$. As
discussed in Section 2.4, the unknown density of $Y_{1}-X_{1}^{\top}\beta$ can
be approximated by smoothing the empirical distribution of the residuals
$(Y_{i}-X_{i}^{\top}\beta)_{i=1}^{n}$ from the current estimate of
$\beta_{0}$. Step II then involves finding an $M$-estimator (1) of $\beta_{0}$
based on the convex loss function induced by the current estimate of
$\psi_{0}^{*}$. In practice, we can apply any suitable convex optimisation
algorithm such as gradient descent, with Newton’s method being a faster
alternative when the score estimate is differentiable. Steps I and II can then
be iterated to convergence.

The following two subsections focus in turn on the linear models in parts
_(a)_ and _(b)_ of Proposition 13. We will analyse a specific version of the
above procedure that is initialised with a pilot estimator $\bar{\beta}_{n}$
of $\beta_{0}$. Provided that $(\bar{\beta}_{n})$ is $\sqrt{n}$-consistent, we
show that a single iteration of Steps I and II yields a semiparametric convex
$M$-estimator of $\beta_{0}$ that achieves ‘antitonic efficiency’ as
$n\to\infty$.

###  3.1 Linear regression with symmetric errors

Under the assumption that $p_{0}$ is symmetric, we first approximate
$\psi_{0}^{*}$ via antitonic projection of kernel-based score estimators. We
do not observe the errors $\varepsilon_{1},\ldots,\varepsilon_{n}$ directly,
so in view of the discussion above, in the algorithm below we use the
residuals from the pilot regression estimator to construct our initial score
estimators. Assume throughout that $n\geq 3$.

  1. 1.

Sample splitting: Partition the observations into three folds indexed by
disjoint $I_{1},I_{2},I_{3}\subseteq[n]$ such that $|I_{1}|=|I_{2}|=\lfloor
n/3\rfloor$ and $|I_{3}|=n-|I_{1}|-|I_{2}|$ respectively. For notational
convenience, let $I_{j+3}:=I_{j}$ for $j\in\\{1,2\\}$.

  2. 2.

Pilot estimators: Let $\psi\colon\mathbb{R}\to\mathbb{R}$ be differentiable,
antisymmetric and strictly decreasing with
$0<\lim_{z\to-\infty}\psi(z)<\infty$. For $j\in\\{1,2,3\\}$, let
$\bar{\beta}_{n}^{(j)}$ be an initial $Z$-estimator satisfying the estimating
equations

$\sum_{i\in
I_{j}}X_{i}\psi\bigl{(}Y_{i}-X_{i}^{\top}\bar{\beta}_{n}^{(j)}\bigr{)}=0.$
(30)

  3. 3.

Antitonic projected score estimation: For $j\in\\{1,2,3\\}$ and $i\in
I_{j+1}$, define out-of-sample residuals
$\hat{\varepsilon}_{i}:=Y_{i}-X_{i}^{\top}\bar{\beta}_{n}^{(j)}$. Letting
$K\colon\mathbb{R}\to[0,\infty)$ be a differentiable kernel and $h\equiv
h_{n}>0$ be a bandwidth that will be specified below, define a kernel density
estimator $\tilde{p}_{n,j}$ of $p_{0}$ by

$\tilde{p}_{n,j}(z):=\frac{1}{|I_{j+1}|}\sum_{i\in
I_{j+1}}K_{h}(z-\hat{\varepsilon}_{i})$

for $z\in\mathbb{R}$, where $K_{h}(\cdot)=h^{-1}K(\cdot/h)$. In addition, let
$\tilde{S}_{n,j}:=\bigl{\\{}z\in\mathbb{R}:|\tilde{p}_{n,j}^{\prime}(z)|\leq\alpha_{n},\,\tilde{p}_{n,j}(z)\geq\gamma_{n}\bigr{\\}}$,
where $\alpha_{n}\in(0,\infty]$ and $\gamma_{n}\in(0,\infty)$ are suitably
chosen truncation parameters, and define
$\tilde{\psi}_{n,j}\colon\mathbb{R}\to\mathbb{R}$ by

$\tilde{\psi}_{n,j}(z):=\frac{\tilde{p}_{n,j}^{\prime}(z)}{\tilde{p}_{n,j}(z)}\mathbbm{1}_{\\{z\in\tilde{S}_{n,j}\\}}.$
(31)

Writing $\tilde{F}_{n,j}$ for the distribution function corresponding to
$\tilde{p}_{n,j}$, let
$\hat{\psi}_{n,j}:=\widehat{\mathcal{M}}_{\mathrm{R}}(\tilde{\psi}_{n,j}\circ\tilde{F}_{n,j}^{-1})\circ\tilde{F}_{n,j}$
be an antitonic projected score estimate, in accordance with (26). Finally,
define an estimator
$\hat{\psi}_{n,j}^{\mathrm{anti}}\in\Psi_{\downarrow}^{\mathrm{anti}}(p_{0})$
of $\psi_{0}^{*}$ by

$\hat{\psi}_{n,j}^{\mathrm{anti}}(z):=\frac{\hat{\psi}_{n,j}(z)-\hat{\psi}_{n,j}(-z)}{2}$

for $z\in\mathbb{R}$.

  4. 4.

Plug-in cross-fitted convex $M$-estimator: For $j\in\\{1,2,3\\}$, let
$\hat{\ell}_{n,j}^{\mathrm{sym}}\colon\mathbb{R}\to\mathbb{R}$ be the induced
convex loss function given by
$\hat{\ell}_{n,j}^{\mathrm{sym}}(z):=-\int_{0}^{z}\hat{\psi}_{n,j}^{\mathrm{anti}}$,
and define

$\hat{\beta}_{n}^{(j)}\in\operatorname*{argmin}_{\beta\in\mathbb{R}^{d}}\sum_{i\in
I_{j+2}}\hat{\ell}_{n,j}^{\mathrm{sym}}(Y_{i}-X_{i}^{\top}\beta)$ (32)

to be a corresponding $M$-estimator of $\beta_{0}$. Finally, let

$\hat{\beta}_{n}^{\dagger}:=\frac{\hat{\beta}_{n}^{(1)}+\hat{\beta}_{n}^{(2)}+\hat{\beta}_{n}^{(3)}}{3}.$
(33)

######  Remark.

For each $j\in\\{1,2,3\\}$, we claim that $\bar{\beta}_{n}^{(j)}$ and
$\hat{\beta}_{n}^{(j)}$ always exist in Steps 2 and 4 respectively. Indeed,
either $\hat{\psi}_{n,j}^{\mathrm{anti}}\equiv 0$, or
$\lim_{z\to-\infty}\hat{\psi}_{n,j}^{\mathrm{anti}}(z)>0>\lim_{z\to\infty}\hat{\psi}_{n,j}^{\mathrm{anti}}(z)$.
In the former case, $\hat{\ell}_{n,j}^{\mathrm{sym}}\equiv 0$ and any
$\beta\in\mathbb{R}^{d}$ minimises the objective function in (32). On the
other hand, in the latter case, $\hat{\ell}_{n,j}^{\mathrm{sym}}$ is a finite,
convex function on $\mathbb{R}$ that is _coercive_ in the sense that
$\hat{\ell}_{n,j}^{\mathrm{sym}}(z)\to\infty$ as $|z|\to\infty$. Thus,
$\theta\equiv(\theta_{1},\dotsc,\theta_{n})\mapsto\sum_{i\in
I_{j+2}}\hat{\ell}_{n,j}^{\mathrm{sym}}(Y_{i}-\theta_{i})=:\mathcal{L}_{n,j}(\theta)$
is also convex and coercive on the column space of the design matrix
$X:=(X_{1}\;\cdots\;X_{n})^{\top}\in\mathbb{R}^{n\times d}$, and hence
$\beta\mapsto\mathcal{L}_{n,j}(X\beta)$ attains its minimum on
$\mathbb{R}^{d}$. The existence of $\bar{\beta}_{n}^{(j)}$ is guaranteed by
similar reasoning. If $\hat{\psi}_{n,j}^{\mathrm{anti}}$ is also continuous,
then any minimiser $\hat{\beta}_{n}^{(j)}$ satisfies the estimating equations
$\sum_{i\in
I_{j+2}}X_{i}\hat{\psi}_{n,j}^{\mathrm{anti}}(Y_{i}-X_{i}^{\top}\hat{\beta}_{n}^{(j)})=0$.

In (32), $\hat{\beta}_{n}^{(j)}$ is certainly not unique if $X$ has does not
have full column rank, but this happens with asymptotically vanishing
probability as $n\rightarrow\infty$ if $\mathbb{E}(X_{1}X_{1}^{\top})$ is
invertible. On the other hand, $\hat{\beta}_{n}^{(j)}$ is unique if $X$ does
have full column rank and $\hat{\psi}_{n,j}^{\mathrm{anti}}$ is strictly
decreasing, in which case $\mathcal{L}_{n,j}$ is strictly convex on the column
space of $X$. In practice, our antitonic score estimators may have constant
pieces. Nevertheless, we emphasise that the conclusion of Theorem 14 below
applies to all sequences of minimisers $(\hat{\beta}_{n}^{(j)})$, so it is
unaffected by non-uniqueness issues.

In the above procedure, we use the observations indexed by $I_{1},I_{2},I_{3}$
to construct the pilot estimator $\bar{\beta}_{n}^{(1)}$, antitonic score
estimate $\hat{\psi}_{n,1}^{\mathrm{anti}}$ and semiparametric $M$-estimator
$\hat{\beta}_{n}^{(1)}$ respectively. Since each fold (specifically the last
one) contains only about one-third of all the data, sample splitting reduces
the efficiency of $\hat{\beta}_{n}^{(1)}$. We remedy this by cross-fitting
(Chernozhukov et al.,, 2018, Definition 3.1), which involves cyclically
permuting the folds to obtain $\hat{\beta}_{n}^{(2)},\hat{\beta}_{n}^{(3)}$
analogously to $\hat{\beta}_{n}^{(1)}$, and then averaging these three
estimators. This reduces the limiting covariance of $\hat{\beta}_{n}^{(1)}$ by
a factor of three in the theory below, where we show that
$\hat{\beta}_{n}^{(1)},\hat{\beta}_{n}^{(2)},\hat{\beta}_{n}^{(3)}$ are
‘asymptotically independent’ in a precise sense.

A different version of cross-fitting (Chernozhukov et al.,, 2018, Definition
3.2) instead averages the empirical risk functions across all three folds, and
outputs a single estimator

$\hat{\beta}_{n}^{\ddagger}\in\operatorname*{argmin}_{\beta\in\mathbb{R}^{d}}\sum_{j=1}^{3}\sum_{i\in
I_{j+2}}\hat{\ell}_{n,j}^{\mathrm{sym}}(Y_{i}-X_{i}^{\top}\beta),$ (34)

whose existence is similarly guaranteed by the convexity of
$\hat{\ell}_{n,j}^{\mathrm{sym}}$ for $j\in\\{1,2,3\\}$.

For a sequence of regression models (27) indexed by $n\in\mathbb{N}$, we make
the following assumptions on the parameters in our procedure.

  1. (A1)

$\mathbb{E}(X_{1}X_{1}^{\top})\in\mathbb{R}^{d\times d}$ is positive definite
and $\max_{i\in[n]}\|{X_{i}}\|\,\alpha_{n}/\gamma_{n}=o_{p}(n^{1/2})$.

  2. (A2)

$\alpha_{n}\to\infty$, $\gamma_{n},h_{n}\to 0$,
$nh_{n}^{3}\gamma_{n}^{2}\to\infty$ and for every fixed $a>0$, we have
$(h_{n}\vee n^{-a})(\alpha_{n}/\gamma_{n})^{2}\to 0$ as $n\to\infty$.

  3. (A3)

The kernel $K$ is non-negative, twice continuously differentiable and
supported on $[-1,1]$.

  4. (A4)

$p_{0}$ is an absolutely continuous density on $\mathbb{R}$ such that
$i(p_{0})<\infty$ and $\int_{\mathbb{R}}|z|^{\delta}\,p_{0}(z)\,dz<\infty$ for
some $\delta>0$.

  5. (A5)

There exists $t_{0}>0$ such that for $t\in\\{-t_{0},t_{0}\\}$, the antitonic
projected score function $\psi_{0}^{*}$ satisfies
$\int_{\mathbb{R}}(\psi_{0}^{*})^{2}(z+t)\,p_{0}(z)\,dz<\infty$.

The conditions in (A2) on the truncation parameters $\alpha_{n},\gamma_{n}$
and bandwidth $h_{n}$ are mild; for instance, we may take
$\alpha_{n}=\gamma_{n}^{-1}=\log n$ and $h_{n}=n^{-b}$ for some $b\in(0,1/3)$.
Our procedure does not require knowledge of the exponent $\delta>0$ in (A4).
Recall from Lemma 1 that $\psi_{0}^{*}(z)$ is finite if and only if
$z\in\mathcal{S}_{0}=\bigl{(}\inf(\operatorname{supp}p_{0}),\sup(\operatorname{supp}p_{0})\bigr{)}$,
so if (A5) holds, then $\mathcal{S}_{0}=\mathbb{R}$ and $\psi_{0}^{*}$ is
necessarily finite-valued on $\mathbb{R}$. The conditions (A4)–(A5) are
satisfied by a variety of commonly-encountered densities $p_{0}$ with
different tail behaviours, ranging from all $t_{\nu}$ densities with $\nu>0$
degrees of freedom (including the Cauchy density as a special case $\nu=1$) to
lighter-tailed Weibull, Laplace, Gaussian and Gumbel densities.

Under the assumptions above, we first establish the $L^{2}(P_{0})$-consistency
of the initial estimates $\tilde{\psi}_{n,j}$ of the score function
$\psi_{0}$, from which it follows that the antitonic functions
$\hat{\psi}_{n,j}$ consistently estimate the population-level projected score
$\psi_{0}^{*}$ in $L^{2}(P_{0})$; see Lemmas 29 and 30 in Section 6.3. This
enables us to prove that our semiparametric convex $M$-estimators
$\hat{\beta}_{n}^{\dagger},\hat{\beta}_{n}^{\ddagger}$ are
$\sqrt{n}$-consistent and have the same limiting Gaussian distribution as the
‘oracle’ convex $M$-estimator
$\hat{\beta}_{\psi_{0}^{*}}:=\operatorname*{argmin}_{\beta\in\mathbb{R}^{d}}n^{-1}\sum_{i=1}^{n}\ell_{0}^{*}(Y_{i}-X_{i}^{\top}\beta)$,
where $\ell_{0}^{*}$ denotes an optimal convex loss function with right
derivative $\psi_{0}^{*}$.

######  Theorem 14.

Suppose that  _(A1)–(A5)_ hold for the linear model (27) with symmetric error
density $p_{0}$. Then for any sequence of estimators $(\hat{\beta}_{n})$ with
$\hat{\beta}_{n}\in\\{\hat{\beta}_{n}^{\dagger},\hat{\beta}_{n}^{\ddagger}\\}$
for each $n$, we have

$\sqrt{n}(\hat{\beta}_{n}-\beta_{0})\overset{d}{\to}N_{d}\biggl{(}0,\frac{\\{\mathbb{E}(X_{1}X_{1}^{\top})\\}^{-1}}{i^{*}(p_{0})}\biggr{)}$

as $n\to\infty$.

######  Remark 15.

Since $X_{1}$ and $\varepsilon_{1}$ are independent, $(X_{1},Y_{1})$ has joint
density $(x,y)\mapsto p_{0}(y-x^{\top}\beta_{0})$ with respect to the product
measure $P_{X}\otimes\mathrm{Leb}$ on $\mathbb{R}^{d}\times\mathbb{R}$, where
we write $P_{X}$ for the distribution of $X_{1}$, and $\mathrm{Leb}$ for
Lebesgue measure on $\mathbb{R}$. Therefore, the score function
$\dot{\ell}_{\beta_{0}}\colon\mathbb{R}^{d}\times\mathbb{R}\to\mathbb{R}^{d}$
is given by

$\dot{\ell}_{\beta_{0}}(x,y):=x\psi_{0}(y-x^{\top}\beta_{0}),$

where $\psi_{0}=p_{0}^{\prime}/p_{0}$. If
$i(p_{0})=\mathbb{E}\bigl{(}\psi_{0}(\varepsilon_{1})^{2}\bigr{)}$ is finite,
then because $\mathbb{E}\psi_{0}(\varepsilon_{1})=0$, the Fisher information
matrix at $\beta_{0}\in\mathbb{R}^{d}$ is

$I(\beta_{0}):=\operatorname{Cov}_{\beta_{0}}\dot{\ell}_{\beta_{0}}(X_{1},Y_{1})=\operatorname{Cov}\\{X_{1}\psi_{0}(\varepsilon_{1})\\}=\mathbb{E}\bigl{\\{}X_{1}X_{1}^{\top}\psi_{0}(\varepsilon_{1})^{2}\bigr{\\}}=\mathbb{E}(X_{1}X_{1}^{\top})\,i(p_{0})\in\mathbb{R}^{d\times
d}.$

Since $i(p_{0})\geq i^{*}(p_{0})>0$ by Theorem 2(d), $I(\beta_{0})$ is
positive definite if and only if $\mathbb{E}(X_{1}X_{1}^{\top})$ is positive
definite. In this case, the convolution and local asymptotic minimax theorems
(van der Vaart,, 1998, Chapter 8) indicate that $\hat{\beta}^{\mathrm{MLE}}$
in (4) has the ‘optimal’ limiting distribution

$N_{d}\bigl{(}0,I(\beta_{0})^{-1}\bigr{)}=N_{d}\biggl{(}0,\frac{\\{\mathbb{E}(X_{1}X_{1}^{\top})\\}^{-1}}{i(p_{0})}\biggr{)}$

among all (regular) sequences of estimators of $\beta_{0}$.

By analogy with the previous display, the limiting covariance in Theorem 14
can be written as the inverse $I^{*}(\beta_{0})^{-1}$ of the antitonic
information matrix
$I^{*}(\beta_{0}):=\mathbb{E}(X_{1}X_{1}^{\top})\,i^{*}(p_{0})$. By Theorem 2,
$1/i^{*}(p_{0})=V_{p_{0}}(\psi_{0}^{*})=\min_{\psi\in\Psi_{\downarrow}(p_{0})}V_{p_{0}}(\psi)$,
so we see from (3) that $I^{*}(\beta_{0})^{-1}$ is the smallest attainable
limiting covariance among all convex $M$-estimators $\hat{\beta}_{\psi}$ based
on a fixed $\psi\in\Psi_{\downarrow}(p_{0})$. We can therefore interpret
$I^{*}(\beta_{0})^{-1}$ as an antitonic efficiency lower bound.

###  3.2 Linear regression with an intercept term

For $d\geq 2$, now consider the linear model

$Y_{i}=\mu_{0}+\tilde{X}_{i}^{\top}\theta_{0}+\varepsilon_{i}\quad\text{for
}i\in[n],$ (35)

where $\mu_{0}$ is an explicit intercept term, so that
$\beta_{0}=(\theta_{0},\mu_{0})$ and $X_{i}=(\tilde{X}_{i}^{\top},1)^{\top}$
in (27) for $i\in[n]$. As mentioned previously, a necessary condition for
$\theta_{0}$ to be identifiable is that

$\mathbb{E}(X_{1}X_{1}^{\top})=\begin{pmatrix}\mathbb{E}(\tilde{X}_{1}\tilde{X}_{1}^{\top})&\\!\\!\\!\mathbb{E}(\tilde{X}_{1})\\\
\mathbb{E}(\tilde{X}_{1})^{\top}&\\!\\!\\!1\end{pmatrix}$

is positive definite, which is equivalent to
$\operatorname{Cov}(\tilde{X}_{1})$, i.e. the Schur complement of $1$ in
$\mathbb{E}(X_{1}X_{1}^{\top})$, being positive definite.

For the functional $\phi\colon(\theta,\mu)\mapsto\theta$ of interest, whose
derivative is represented by
$D\phi(\theta,\mu)=(I_{d-1}\;\;0)\in\mathbb{R}^{(d-1)\times d}$, the
information lower bound is the top-left $(d-1)\times(d-1)$ submatrix of
$I(\beta_{0})^{-1}$, namely

$D\phi(\beta_{0})\;I(\beta_{0})^{-1}\;D\phi(\beta_{0})^{\top}=\frac{\operatorname{Cov}(\tilde{X}_{1})^{-1}}{i(p_{0})}\in\mathbb{R}^{(d-1)\times(d-1)};$
(36)

see van der Vaart, (1998, Chapters 8 and 25.3). In Section 6.3.2, we verify
that the efficient score function (van der Vaart,, 1998, Chapter 25.4)
$\tilde{\ell}_{\beta_{0}}\colon\mathbb{R}^{d}\times\mathbb{R}\to\mathbb{R}^{d-1}$
is given by

$\tilde{\ell}_{\beta_{0}}(x,y):=\bigl{(}\tilde{x}-\mathbb{E}(\tilde{X}_{1})\bigr{)}\cdot\psi_{0}(y-x^{\top}\beta_{0}),$
(37)

where $x=(x_{1},\ldots,x_{d})$ and $\tilde{x}=(x_{1},\ldots,x_{d-1})$.

We will construct an adaptive convex $M$-estimator of $\beta_{0}$ that
asymptotically achieves antitonic efficiency. Similarly to Section 3.1, we
employ three-fold cross-fitting with the convention $I_{j+3}=I_{j}$ for
$j\in\\{1,2\\}$. We follow Steps 1 and 2 of the previous procedure to obtain
pilot estimators
$\bar{\beta}_{n}^{(j)}\equiv(\bar{\theta}_{n}^{(j)},\bar{\mu}_{n}^{(j)})$ of
$\beta_{0}=(\theta_{0},\mu_{0})$ satisfying $\sum_{i\in
I_{j}}X_{i}\psi\bigl{(}Y_{i}-X_{i}^{\top}\bar{\beta}_{n}^{(j)}\bigr{)}=0$ for
$j\in\\{1,2,3\\}$, but make the following modifications to subsequent steps.

  1. 3′.

Antitonic projected score estimation: For $j\in\\{1,2,3\\}$ and $i\in
I_{j+1}$, use the out-of-sample residuals
$\hat{\varepsilon}_{i}:=Y_{i}-X_{i}^{\top}\bar{\beta}_{n}^{(j)}$ to construct
the initial kernel-based score estimator $\tilde{\psi}_{n,j}$ and its
antitonic projection
$\hat{\psi}_{n,j}:=\widehat{\mathcal{M}}_{\mathrm{R}}(\tilde{\psi}_{n,j}\circ\tilde{F}_{n,j}^{-1})\circ\tilde{F}_{n,j}$
as before. Since $p_{0}$ is not symmetric in general, we estimate
$\psi_{0}^{*}$ by $\hat{\psi}_{n,j}$ instead of
$\hat{\psi}_{n,j}^{\mathrm{anti}}$.

  2. 4′.

Plug-in cross-fitted convex $M$-estimator: For $j\in\\{1,2,3\\}$, let
$\hat{\ell}_{n,j}\colon\mathbb{R}\to\mathbb{R}$ be the induced convex loss
function given by $\hat{\ell}_{n,j}(z):=-\int_{0}^{z}\hat{\psi}_{n,j}$, and
define

$\displaystyle\hat{\theta}_{n}^{(j)}\in\operatorname*{argmin}_{\theta\in\mathbb{R}^{d-1}}\sum_{i\in
I_{j+2}}\hat{\ell}_{n,j}\bigl{(}Y_{i}-\bar{\mu}_{n}^{(j)}-\bar{X}_{n,j}^{\top}\bar{\theta}_{n}^{(j)}-(\tilde{X}_{i}-\bar{X}_{n,j})^{\top}\theta\bigr{)},$
(38)

where $\bar{X}_{n,j}:=|I_{j+2}|^{-1}\sum_{i\in I_{j+2}}\tilde{X}_{i}$.
Finally, let
$\hat{\theta}_{n}^{\dagger}:=\bigl{(}\hat{\theta}_{n}^{(1)}+\hat{\theta}_{n}^{(2)}+\hat{\theta}_{n}^{(3)}\bigr{)}/3$.
Alternatively, define

$\hat{\theta}_{n}^{\ddagger}\in\operatorname*{argmin}_{\theta\in\mathbb{R}^{d-1}}\sum_{j=1}^{3}\sum_{i\in
I_{j+2}}\hat{\ell}_{n,j}\bigl{(}Y_{i}-\bar{\mu}_{n}^{(j)}-\bar{X}_{n,j}^{\top}\bar{\theta}_{n}^{(j)}-(\tilde{X}_{i}-\bar{X}_{n,j})^{\top}\theta\bigr{)}.$

For the same reasons as in Section 3.1, there exists a minimiser
$\hat{\theta}_{n}^{(j)}$ in (38) if

$\lim_{z\to-\infty}\hat{\psi}_{n,j}(z)>0>\lim_{z\to\infty}\hat{\psi}_{n,j}(z),$
(39)

and $\hat{\theta}_{n}^{(j)}$ is unique if $\hat{\psi}_{n,j}$ is strictly
decreasing and the design matrix has full column rank. If $\hat{\psi}_{n,j}$
is continuous, then any minimiser $\hat{\theta}_{n}^{(j)}$ satisfies

$\sum_{i\in
I_{j+2}}(\tilde{X}_{i}-\bar{X}_{n,j})\cdot\hat{\psi}_{n,j}\bigl{(}Y_{i}-\bar{\mu}_{n}^{(j)}-\bar{X}_{n,j}^{\top}\bar{\theta}_{n}^{(j)}-(\tilde{X}_{i}-\bar{X}_{n,j})^{\top}\hat{\theta}_{n}^{(j)}\bigr{)}=0.$

This is a variant of the efficient score equations (van der Vaart,, 1998,
Chapter 25.8) based on (37) and suitable estimators of the nuisance parameters
$\mu_{0},\psi_{0}^{*}$ respectively. Since we do not assume knowledge of the
population mean $\mathbb{E}(\tilde{X}_{1})\in\mathbb{R}^{d-1}$, we replace it
with its sample analogue $\bar{X}_{n,j}$ in each of the three folds.

Under the regularity conditions in Section 3.1, it turns out that the
projected score functions $\hat{\psi}_{n,j}$ are also
$L^{2}(P_{0})$-consistent estimators of $\psi_{0}^{*}$ in this setting. It
follows that with probability tending to 1 as $n\to\infty$, (39) holds and
hence $\hat{\theta}_{n}^{\dagger}$ and $\hat{\theta}_{n}^{\ddagger}$ exist. We
then adapt the proof strategy for Theorem 14 above to establish a similar
antitonic efficiency result.

######  Theorem 16.

Suppose that  _(A1)–(A5)_ hold in the linear model (35). Then for any sequence
of estimators $(\hat{\theta}_{n})$ such that
$\hat{\theta}_{n}\in\\{\hat{\theta}_{n}^{\dagger},\hat{\theta}_{n}^{\ddagger}\\}$
for each $n$, we have

$\sqrt{n}(\hat{\theta}_{n}-\theta_{0})\overset{d}{\to}N_{d-1}\biggl{(}0,\frac{\operatorname{Cov}(\tilde{X}_{1})^{-1}}{i^{*}(p_{0})}\biggr{)}$

as $n\to\infty$.

Zou and Yuan, (2008) introduced a robust alternative to least squares called
the composite quantile regression estimator, which remains
$\sqrt{n}$-consistent when the error variance is infinite but also has
asymptotic efficiency at least 70% that of OLS under suitable conditions on
the error density $p_{0}$. This is achieved by borrowing strength across
several quantiles of the conditional distribution of $Y_{1}$ given $X_{1}$,
rather than targeting just a single quantile (e.g. the conditional median)
using the quantile loss $\ell_{\tau}\colon\mathbb{R}\to\mathbb{R}$ given by
$\ell_{\tau}(z):=(\tau-\mathbbm{1}_{\\{z<0\\}})z$ for $\tau\in(0,1)$. More
precisely, for $K\in(0,\infty)$, let $\tau_{k}:=k/(K+1)$ for $k\in[K]$ and
define

$(\hat{\theta}_{n}^{\mathrm{CQ},K},\hat{\mu}_{n,1},\dotsc,\hat{\mu}_{n,K})\in\operatorname*{argmin}_{(\theta,\mu_{1},\dotsc,\mu_{K})}\,\sum_{k=1}^{K}\sum_{i=1}^{n}\ell_{\tau_{k}}(Y_{i}-\mu_{k}-\tilde{X}_{i}^{\top}\theta),$
(40)

where the $\operatorname*{argmin}$ is taken over all
$\theta\in\mathbb{R}^{d-1}$ and $\mu_{1},\dotsc,\mu_{K}\in\mathbb{R}$. Under
assumption (2) of Zou and Yuan, (2008) in our random design setting, it
follows from their Theorems 2.1 and 3.1 that

$\sqrt{n}(\hat{\theta}_{n}^{\mathrm{CQ},K}-\theta_{0})\overset{d}{\to}N_{d-1}\bigl{(}0,V_{p_{0},\mathrm{CQ},K}\cdot\operatorname{Cov}(\tilde{X}_{1})^{-1}\bigr{)},$

where, writing $J_{0}=p_{0}\circ F_{0}^{-1}$ for the density quantile function
of the errors, we have

$V_{p_{0},\mathrm{CQ},K}:=\frac{\sum_{k,k^{\prime}=1}^{K}\tau_{k\wedge
k^{\prime}}(1-\tau_{k\vee
k^{\prime}})}{\bigl{(}\sum_{k=1}^{K}J_{0}(\tau_{k})\bigr{)}^{2}}\to\frac{1}{12\,\bigl{(}\int_{0}^{1}J_{0}\bigr{)}^{2}}=:V_{p_{0},\mathrm{CQ}}$

as $K\to\infty$. For every $p_{0}$ satisfying their assumptions, Zou and Yuan,
(2008, Theorem 3.1) established that the CQR estimator (in the notional limit
$K\to\infty$) has asymptotic relative efficiency

$\frac{V_{p_{0}}(\psi_{\mathrm{OLS}})}{V_{p_{0},\mathrm{CQ}}}=12\biggl{(}\int_{0}^{1}J_{0}\biggr{)}^{2}\int_{\mathbb{R}}z^{2}\,p_{0}(z)\,dz>\frac{6}{e\pi}\approx
0.703$

relative to OLS, where $\psi_{\mathrm{OLS}}(z):=-z$ for $z\in\mathbb{R}$.
Nevertheless, Theorem 16 and the following result imply that the asymptotic
covariance of the ‘limiting’ CQR estimator is always at least that of the
semiparametric convex $M$-estimator $\hat{\theta}_{n}$ in our framework.
Moreover, the former estimator can have arbitrarily low efficiency relative to
the latter, even when $p_{0}$ is log-concave.

######  Lemma 17.

For every uniformly continuous density $p_{0}$, we have
$V_{p_{0},\mathrm{CQ}}\geq 1/i^{*}(p_{0})$, with equality if and only if
either $i^{*}(p_{0})=\infty$ or $p_{0}$ is a logistic density of the form

$p_{0}(z)=\frac{\lambda e^{-\lambda(z-\mu)}}{(1+e^{-\lambda(z-\mu)})^{2}}$

for $z\in\mathbb{R}$, where $\mu\in\mathbb{R}$ and $\lambda>0$. Moreover,

$\inf_{p_{0}\in\mathcal{P}_{\mathrm{LC}}:i(p_{0})<\infty}\,\frac{1}{i^{*}(p_{0})\,V_{p_{0},\mathrm{CQ}}}=0.$

###  3.3 Inference

To perform asymptotically valid inference for $\beta_{0}$ based on Theorems 14
and 16, we require a consistent estimator of the antitonic information
$i^{*}(p_{0})$ when $p_{0}$ is unknown.

######  Lemma 18.

Suppose that  _(A1)–(A5)_ hold for a linear model (27) in which either $p_{0}$
is symmetric or $X_{1d}=1$ almost surely. For $j\in\\{1,2,3\\}$, letting
$\breve{\varepsilon}_{i}:=Y_{i}-X_{i}^{\top}\bar{\beta}_{n}^{(j)}$ for $i\in
I_{j+2}$, we have

$\hat{\imath}_{n}:=\frac{1}{n}\sum_{j=1}^{3}\sum_{i\in
I_{j+2}}\hat{\psi}_{n,j}(\breve{\varepsilon}_{i})^{2}\overset{p}{\to}i^{*}(p_{0})$

as $n\to\infty$.

The antitonic information matrix $I^{*}(\beta_{0})$ can therefore be estimated
consistently by the observed antitonic information matrix

$\hat{I}_{n}:=\frac{\hat{\imath}_{n}}{n}\sum_{i=1}^{n}X_{i}X_{i}^{\top}=\frac{\hat{\imath}_{n}}{n}X^{\top}X,$

where $X=(X_{1}\;\cdots\;X_{n})\in\mathbb{R}^{n\times d}$ has full column rank
with probability tending to 1 under (A1). It follows from Theorem 14 and Lemma
18 that

$\sqrt{n}\hat{I}_{n}^{1/2}(\hat{\beta}_{n}-\beta_{0})\overset{d}{\to}N_{d}(0,I_{d})$

when $p_{0}$ is symmetric. Thus, writing $\hat{\beta}_{n,j}$ for the $j$th
component of $\hat{\beta}_{n}$ and $\hat{v}_{j}:=(\hat{I}_{n}^{-1})_{jj}$ for
the $j$th diagonal entry of $\hat{I}_{n}^{-1}$, we have that

$\biggl{[}\hat{\beta}_{n,j}-z_{\alpha/2}\sqrt{\frac{\hat{v}_{j}}{n}}\,,\,\hat{\beta}_{n,j}+z_{\alpha/2}\sqrt{\frac{\hat{v}_{j}}{n}}\biggl{]}$

is an asymptotic $(1-\alpha)$-level confidence interval for the $j$th
component of $\beta_{0}$, where $z_{\alpha/2}$ denotes the
$(1-\alpha/2)$-quantile of the standard normal distribution. Moreover,

$\bigl{\\{}b\in\mathbb{R}^{d}:n(\hat{\beta}_{n}-b)^{\top}\hat{I}_{n}(\hat{\beta}_{n}-b)\leq\chi_{d}^{2}(\alpha)\bigr{\\}}$

is an asymptotic $(1-\alpha)$-confidence ellipsoid for $\beta_{0}$, where
$\chi_{d}^{2}(\alpha)$ denotes the $(1-\alpha)$-quantile of the $\chi_{d}^{2}$
distribution. Similarly, when $X_{1d}=1$ almost surely, define
$\widetilde{I}_{n}:=n^{-1}\sum_{i=1}^{n}\hat{\imath}_{n}(\tilde{X}_{i}-\bar{X}_{n})(\tilde{X}_{i}-\bar{X}_{n})^{\top}$,
where $\bar{X}_{n}:=n^{-1}\sum_{i=1}^{n}\tilde{X}_{i}$, and let
$\tilde{v}_{j}:=(\widetilde{I}_{n}^{-1})_{jj}$ for $j\in[d-1]$. Then by
Theorem 16 and Lemma 18,

$\biggl{[}\hat{\theta}_{n,j}-z_{\alpha/2}\sqrt{\frac{\tilde{v}_{j}}{n}}\,,\,\hat{\theta}_{n,j}+z_{\alpha/2}\sqrt{\frac{\tilde{v}_{j}}{n}}\biggl{]}$

is an asymptotic $(1-\alpha)$-level confidence interval for the $j$th
component of $\theta_{0}$ for each $j\in[d-1]$, and

$\displaystyle\bigl{\\{}v\in\mathbb{R}^{d-1}:n(\hat{\theta}_{n}-v)^{\top}\widetilde{I}_{n}(\hat{\theta}_{n}-v)\leq\chi_{d-1}^{2}(\alpha)\bigr{\\}}$
(41)

is an asymptotic $(1-\alpha)$-confidence ellipsoid for $\theta_{0}$. Lemma 18
also ensures that standard linear model diagnostics, either based on
heuristics such as Cook’s distances (Cook,, 1977), or formal goodness-of-fit
tests (Janková et al.,, 2020), can be applied.

##  4 Numerical experiments

In our numerical experiments, we generate covariates
$\tilde{X}_{1},\dotsc,\tilde{X}_{n}\overset{\mathrm{iid}}{\sim}N_{d-1}(\mathbf{1}_{d-1},I_{d-1})$,
where $\mathbf{1}_{d-1}$ denotes a $(d-1)$-dimensional all-ones vector, and
responses $Y_{1},\dotsc,Y_{n}$ according to the linear model (35) with
$\|\theta_{0}\|_{2}=3$. We generate independent errors
$\varepsilon_{1},\ldots,\varepsilon_{n}\overset{\mathrm{iid}}{\sim}P_{0}$ for
the following choices of $P_{0}$:

  1. (i)

$P_{0}$ is standard Gaussian.

  2. (ii)

$P_{0}$ is standard Cauchy.

  3. (iii)

Gaussian scale mixture: $P_{0}=\frac{1}{2}N(0,1)+\frac{1}{2}N(0,16)$.

  4. (iv)

Gaussian location mixture:
$P_{0}=\frac{1}{2}N\bigl{(}-\frac{3}{2},\frac{1}{100}\bigr{)}+\frac{1}{2}N\bigl{(}\frac{3}{2},\frac{1}{100}\bigr{)}$.
This distribution is similar to the one constructed in the proof of
Proposition 6 to show that log-concave maximum likelihood estimation of the
error distribution may result in arbitrarily large efficiency loss.

  5. (v)

Smoothed uniform: $P_{0}$ is the distribution of $U+\frac{1}{10}Z$, where
$U\sim\mathrm{Unif}[-1,1]$ and $Z\sim N(0,1)$ are independent.

  6. (vi)

Smoothed exponential: $P_{0}$ is the distribution of
$W-1+\frac{\sqrt{3}}{10}Z$, where $W\sim\operatorname{Exp}(1)$ and $Z\sim
N(0,1)$ are independent. We choose the standard deviation
$\frac{\sqrt{3}}{10}$ for the Gaussian component so that the ratio of the
variances of the non-Gaussian and the Gaussian components is the same as that
in the smoothed uniform setting.

In cases (i), (v) and (vi) $P_{0}$ is log-concave, while for the other
settings it is not.

| Oracle | ASM | Alt | LCMLE | 1S | LAD | OLS
---|---|---|---|---|---|---|---
Standard Gaussian | 8.34 | 8.88 | 8.86 | 10.73 | 9.72 | 13.11 | 8.34
Standard Cauchy | 20.07 | $\bf{20.61}$ | 20.76 | 23.85 | 22.13 | 21.36 | $2.03\times 10^{6}$
Gaussian scale mixture | 31.36 | $\bf{32.01}$ | 32.33 | 34.89 | 36.18 | 34.97 | 72.34
Gaussian location mixture | 0.16 | 0.17 | $\bf{0.16}$ | 1.51 | 18.07 | 319.65 | 18.51
Smoothed uniform | 1.02 | 1.29 | $\bf{1.15}$ | 1.52 | 2.07 | 7.92 | 2.78
Smoothed exponential | 1.78 | 2.13 | $\bf{2.02}$ | 2.54 | 3.36 | 8.27 | 8.65
Table 1: Squared estimation error ($\times 10^{3}$) for different estimators,
with $n=600$ and $d=6$.

Figure 7: Plots of the average squared error loss ($\times 10^{3}$) of
different estimators for noise distributions (i)–(vi), with $n=600$ and $d=6$.
In each plot, the red dashed line indicates the corresponding value for the
oracle convex $M$-estimator, and we omit the estimators that have very large
estimation error (see Table 1 for full details).

We compared the performance of two versions of our procedure with an oracle
approach and four existing methods. The first variant of our procedure, which
we refer to as ASM (antitonic score matching) in all of the plots, is as
described in Section 3, except that we do not perform sample splitting, cross-
fitting or truncation of the initial score estimates. These devices are
convenient for theoretical analysis but not essential in practice. More
precisely, ASM first constructs a pilot estimator
$(\bar{\theta}_{n},\bar{\mu}_{n})$, and then uses the vector of residuals
$\bigl{(}\hat{\varepsilon}_{1},\ldots,\hat{\varepsilon}_{n}\bigr{)}$, where
$\hat{\varepsilon}_{i}:=Y_{i}-\bar{\mu}_{n}-\tilde{X}_{i}^{\top}\bar{\theta}_{n}$,
to obtain an initial kernel-based score estimator $\tilde{\psi}_{n}$, formed
using a Gaussian kernel and the default Silverman’s choice of bandwidth
(Silverman,, 1986, p. 48). Following Step 3′ in Section 3.2, we estimate the
antitonic projected score and the corresponding convex loss function by
$\hat{\psi}_{n}:=\widehat{\mathcal{M}}_{\mathrm{R}}(\tilde{\psi}_{n}\circ\tilde{F}_{n}^{-1})\circ\tilde{F}_{n}$
and $\hat{\ell}_{n}$ respectively, where $\tilde{F}_{n}$ denotes the
distribution function associated with the kernel density estimate. Finally, we
use Newton’s algorithm to compute our semiparametric estimator

$\hat{\theta}_{n}^{\mathrm{ASM}}\in\operatorname*{argmin}_{\theta\in\mathbb{R}^{d-1}}\sum_{i=1}^{n}\hat{\ell}_{n}\bigl{(}Y_{i}-\bar{\mu}_{n}-\bar{X}_{n}^{\top}\bar{\theta}_{n}-(\tilde{X}_{i}-\bar{X}_{n})^{\top}\theta\bigr{)},$

where $\bar{X}_{n}:=n^{-1}\sum_{i=1}^{n}\tilde{X}_{i}$.

In the second version of our procedure, which we refer to as Alt in our plots,
we implement the empirical analogue of the alternating optimisation procedure
described at the beginning of Section 3. We start with an uninformative
initialiser
$(\hat{\theta}_{n}^{(0)},\hat{\mu}_{n}^{(0)})=(0,0)\in\mathbb{R}^{d-1}\times\mathbb{R}$,
and then alternate between the following steps for $t\in\mathbb{N}$:

  1. I.

Compute residuals
$\hat{\varepsilon}_{i}^{(t-1)}:=Y_{i}-\hat{\mu}_{n}^{(t-1)}-\tilde{X}_{i}^{\top}\hat{\theta}_{n}^{(t-1)}$
for $i\in[n]$ and hence estimate the antitonic projected score
$\hat{\psi}_{n}^{(t-1)}$ and the corresponding convex loss
$\hat{\ell}_{n}^{(t-1)}$ as for ASM.

  2. II.

Update
$(\hat{\theta}_{n}^{(t)},\hat{\mu}_{n}^{(t)})\in\operatorname*{argmin}_{(\theta,\mu)\in\mathbb{R}^{d-1}\times\mathbb{R}}\sum_{i=1}^{n}\hat{\ell}_{n}^{(t-1)}(Y_{i}-\mu-\tilde{X}_{i}^{\top}\theta)$.

We iterate these steps until convergence of the empirical score matching
objective

$\hat{D}_{n}(\hat{\psi}_{n}^{(t)};\hat{\varepsilon}_{1}^{(t)},\dotsc,\hat{\varepsilon}_{n}^{(t)})=\frac{1}{n}\sum_{i=1}^{n}\bigl{\\{}\hat{\psi}_{n}^{(t)}(\hat{\varepsilon}_{i}^{(t)})+2(\hat{\psi}_{n}^{(t)})^{\prime}(\hat{\varepsilon}_{i}^{(t)})\bigr{\\}}$

defined in (24). We note that the iterates $\hat{\mu}_{n}^{(t)}$ and
$\hat{\psi}_{n}^{(t)}$ are not guaranteed to converge due to non-
identifiability of the intercept term $\mu_{0}$ in (35), but
$\hat{\theta}_{n}^{(t)}$ and the score matching objective values did indeed
converge in all of our experiments.

The alternative approaches that we consider are as follows:

  * •

Oracle: The $M$-estimator $\hat{\theta}_{n}^{\mathrm{oracle}}$, where

$\begin{pmatrix}\hat{\theta}_{n}^{\mathrm{oracle}}\\\
\hat{\mu}_{n}^{\mathrm{oracle}}\end{pmatrix}\in\operatorname*{argmin}_{(\theta,\mu)\in\mathbb{R}^{d-1}\times\mathbb{R}}\;\sum_{i=1}^{n}\ell_{0}^{*}(Y_{i}-\mu-\tilde{X}_{i}^{\top}\theta)$

is defined with respect to the optimal convex loss function $\ell_{0}^{*}$.
Although $\hat{\theta}_{n}^{\mathrm{oracle}}$ asymptotically attains antitonic
efficiency in the sense of Theorem 16, it is not a valid estimator in our
semiparametric framework since it requires knowledge of $p_{0}$.

  * •

LAD: The least absolute deviation estimator $\hat{\theta}_{n}^{\mathrm{LAD}}$,
where

$\begin{pmatrix}\hat{\theta}_{n}^{\mathrm{LAD}}\\\
\hat{\mu}_{n}^{\mathrm{LAD}}\end{pmatrix}\in\operatorname*{argmin}_{(\theta,\mu)\in\mathbb{R}^{d-1}\times\mathbb{R}}\;\sum_{i=1}^{n}|Y_{i}-\mu-\tilde{X}_{i}^{\top}\theta|.$

  * •

OLS: The ordinary least squares estimator $\hat{\theta}_{n}^{\mathrm{OLS}}$,
where

$\begin{pmatrix}\hat{\theta}_{n}^{\mathrm{OLS}}\\\
\hat{\mu}_{n}^{\mathrm{OLS}}\end{pmatrix}\in\operatorname*{argmin}_{(\theta,\mu)\in\mathbb{R}^{d-1}\times\mathbb{R}}\;\sum_{i=1}^{n}(Y_{i}-\mu-\tilde{X}_{i}^{\top}\theta)^{2}.$

  * •

1S: The semiparametric one-step method, where we start with a pilot estimator
$(\bar{\theta}_{n},\bar{\mu}_{n})$, compute a (not necessarily decreasing)
nonparametric score estimate, and then update $\bar{\theta}_{n}$ with a single
Newton step instead of solving the estimating equations exactly. Our
implementation follows van der Vaart, (1998, Chapter 25.8): we split the data
into two folds of equal size indexed by $I_{1}$ and $I_{2}$, and then use the
residuals
$\hat{\varepsilon}_{i}:=Y_{i}-\tilde{X}_{i}^{\top}\bar{\theta}_{n}-\bar{\mu}_{n}$
for $i\in I_{1}$ and $i\in I_{2}$ separately to obtain kernel density
estimates $\hat{p}_{n,1},\hat{p}_{n,2}$ of $p_{0}$ (constructed as for ASM).
Defining the score estimates
$\hat{\psi}_{n,j}:=\hat{p}_{n,j}^{\prime}/\hat{p}_{n,j}$ for $j\in\\{1,2\\}$,
we output the cross-fitted estimator $\hat{\theta}_{n}^{\mathrm{1S}}$, where

$\displaystyle\begin{pmatrix}\hat{\theta}_{n}^{\mathrm{1S}}\\\
\hat{\mu}_{n}^{\mathrm{1S}}\end{pmatrix}:=\begin{pmatrix}\bar{\theta}_{n}\\\
\bar{\mu}_{n}\end{pmatrix}-$ $\displaystyle\biggl{(}\sum_{i\in
I_{1}}\hat{\psi}_{n,2}(\hat{\varepsilon}_{i})^{2}X_{i}X_{i}^{\top}+\sum_{i\in
I_{2}}\hat{\psi}_{n,1}(\hat{\varepsilon}_{i})^{2}X_{i}X_{i}^{\top}\biggr{)}^{-1}$
$\displaystyle\hskip 142.26378pt\biggl{(}\sum_{i\in
I_{1}}\hat{\psi}_{n,2}(\hat{\varepsilon}_{i})X_{i}+\sum_{i\in
I_{2}}\hat{\psi}_{n,1}(\hat{\varepsilon}_{i})X_{i}\biggr{)}.$

  * •

LCMLE: We estimate the error density $p_{0}$ using the log-concave maximum
likelihood estimator (Cule et al.,, 2010; Dümbgen et al.,, 2011, 2013). More
precisely, again writing $\mathcal{P}_{\mathrm{LC}}$ for the set of univariate
log-concave densities, start with a pilot estimator
$(\hat{\theta}_{n}^{(0)},\hat{\mu}_{n}^{(0)})$ and alternate the following two
steps for $t\in\mathbb{N}$:

$\displaystyle\hat{p}_{n}^{(t)}$
$\displaystyle\in\operatorname*{argmax}_{p\in\mathcal{P}_{\mathrm{LC}}}\;\sum_{i=1}^{n}\,\log
p(Y_{i}-\hat{\mu}_{n}^{(t-1)}-\tilde{X}_{i}^{\top}\hat{\theta}_{n}^{(t-1)})$
$\displaystyle\begin{pmatrix}\hat{\theta}_{n}^{(t)}\\\
\hat{\mu}_{n}^{(t)}\end{pmatrix}$
$\displaystyle\in\operatorname*{argmax}_{(\theta,\mu)\in\mathbb{R}^{d-1}\times\mathbb{R}}\;\sum_{i=1}^{n}\,\log\hat{p}_{n}^{(t)}(Y_{i}-\mu-\tilde{X}_{i}^{\top}\theta).$

Where required, we took the pilot estimator to be the least absolute deviation
estimator for all methods except in the Gaussian location mixture setting
(iv). Here, we chose
$(\hat{\theta}_{n}^{\mathrm{OLS}},\hat{\mu}_{n}^{\mathrm{OLS}})$ to be the
pilot estimator since $p_{0}(0)$ is close to 0 and hence
$V_{p_{0}}(\psi)=1/\bigl{(}4p_{0}(0)^{2}\bigr{)}$ is very large for
$\hat{\theta}_{n}^{\mathrm{LAD}}$, where
$\psi(\cdot)=-\operatorname{sgn}(\cdot)$.

$n$ | 600 | 1200 | 2400
---|---|---|---
ASM | 0.26 | 0.61 | 1.93
Alt | 4.35 | 9.99 | 27.85
LCMLE | 0.64 | 2.21 | 8.34
1S | 0.10 | 0.15 | 0.29
Standard Gaussian

$n$ | 600 | 1200 | 2400
---|---|---|---
ASM | 0.27 | 0.61 | 1.81
Alt | 5.03 | 10.52 | 28.64
LCMLE | 0.71 | 1.36 | 2.72
1S | 0.11 | 0.15 | 0.29
Gaussian scale mixture

$n$ | 600 | 1200 | 2400
---|---|---|---
ASM | 0.21 | 0.58 | 1.96
Alt | 1.84 | 3.72 | 11.36
LCMLE | 3.15 | 7.09 | 9.29
1S | 0.12 | 0.15 | 0.29
Gaussian location mixture

Table 2: Mean execution time in seconds for different semiparametric
estimators when $d=100$.

Figure 8: Plot of $n\times(\text{average squared loss})$ for different
estimators, with $d=100$ and sample sizes $n\in\\{600,1200,2400\\}$
corresponding to the yellow, green, and pink bars respectively. The black,
blue and red dashed lines indicate the values of $n\times(\text{average
squared loss})$ for the oracle convex $M$-estimator when $n=600,1200,2400$
respectively.

Figure 9: Comparison of the mean squared errors of different convex
$M$-estimators for sample sizes $n\in\\{50,100,150,200,250,300\\}$ and $d=6$.

Figure 10: Log-log plots of the average squared error loss against sample size
for three convex $M$-estimators. In each experiment, we set
$n\in\\{800,1600,3200,6400\\}$ and $d=2$.

Figure 11: Q-Q plots of the standardised errors of estimates of two components
of $\theta_{0}$. The blue and green vertical lines mark the 90% and 95%
theoretical quantiles respectively. In each experiment, we set $n=600$ and
$d=4$, and perform 8000 repetitions.

(a) Standard normal

(b) Gaussian mixture

Figure 12: The projection of the $95\%$ confidence ellipsoid of
$\hat{\theta}^{\mathrm{ASM}}$ onto the first two dimensions (blue) and the
projection of the $95\%$ confidence ellipsoid of
$\hat{\theta}_{n}^{\mathrm{OLS}}$ onto the first two dimensions (black
dashed). The blue dot, black circle and red star correspond to
$\hat{\theta}^{\mathrm{ASM}}$, $\hat{\theta}_{n}^{\mathrm{OLS}}$, and
$\theta_{0}$ respectively.

| $\theta_{0,1}$ | $\theta_{0,2}$ | $\theta_{0,3}$
---|---|---|---
95% CI | 0.965 | 0.965 | 0.964
90% CI | 0.924 | 0.925 | 0.924
Standard Cauchy

| $\theta_{0,1}$ | $\theta_{0,2}$ | $\theta_{0,3}$
---|---|---|---
95% CI | 0.952 | 0.955 | 0.955
90% CI | 0.905 | 0.908 | 0.908
Gaussian mixture

| $\theta_{0,1}$ | $\theta_{0,2}$ | $\theta_{0,3}$
---|---|---|---
95% CI | 0.950 | 0.954 | 0.952
90% CI | 0.903 | 0.902 | 0.901
Standard Gaussian

Table 3: Empirical coverage of the confidence intervals for each coordinate of
$\theta_{0}$. In each experiment, we set $n=600$ and $d=4$, and perform 8000
repetitions.

| Cauchy | Mixture | Gaussian
---|---|---|---
$i^{*}(p_{0})$ | 0.44 | 0.46 | 1.0
RMSE($\hat{\jmath}_{n}$) | 0.01 | 0.05 | 0.1
(a) Root mean squared errors of $\hat{\jmath}_{n}$.

| Cauchy | Mixture | Gaussian
---|---|---|---
$\dfrac{\mathrm{Vol}(\hat{C}_{n}^{\mathrm{ASM}})}{\mathrm{Vol}(\hat{C}_{n}^{\mathrm{OLS}})}$ | 0.001 | 0.49 | 1.14
(b) Mean ratio of volumes of $\hat{C}_{n}^{\mathrm{ASM}}$ and
$\hat{C}_{n}^{\mathrm{OLS}}$.

Table 4: Accuracy of estimates of the antitonic information, and comparison of
the volumes of the two confidence ellipsoids. In each experiment, we set
$n=600$ and $d=4$, and perform 8000 repetitions.

###  4.1 Estimation accuracy

In all of the experiments in this subsection, we set $d=6$, $\mu_{0}=2$ and
drew $\theta_{0}$ uniformly at random from the centred Euclidean sphere in
$\mathbb{R}^{d-1}$ of radius 3. For each of the estimators above, the average
squared Euclidean norm errors $\|{\hat{\theta}_{n}-\theta_{0}}\|^{2}$ were
computed over 200 repetitions. In the first experiment, we consider each error
distribution in turn and compare the average squared estimation error when
$n=600$. The results are presented in Table 1 and Figure 7. We observe that
our proposed procedures ASM and Alt have the lowest estimation error except in
the case of standard Gaussian error, where OLS coincides with the oracle
convex loss estimator and has a slightly lower estimation error. It is
interesting that the one-step estimator (1S) can perform very poorly in finite
samples, and in particular, may barely improve on its initialiser. In all
settings considered, ASM and Alt have comparable error.

Next, we investigate the running time and estimation accuracy of the
semiparametric estimators ASM, Alt, LCMLE and 1S for moderately large sample
sizes $n\in\\{600,1200,2400\\}$ and dimension $d=100$. We consider in turn the
error distributions (i), (iii) and (iv); see Table 2 and Figure 8. We see that
ASM in particular is competitive in terms of its running time, and that the
performance of our procedures does not deteriorate relative to its competitors
for this larger choice of $d$.

Our third numerical experiment demonstrates that the estimation error of ASM
is comparable to that of the oracle convex $M$-estimator even for small sample
sizes. For $n\in\\{50,100,150,200,300\\}$ and $d=6$, and error distributions
(i), (ii) and (iv), we record in Figure 9 the average estimation error of
Oracle, ASM, and LCMLE. With the exception of the smallest sample size $n=50$
in the Gaussian location mixture setting, the estimation error of ASM tracks
that of the oracle very closely; on the other hand, LCMLE is somewhat
suboptimal, particularly in the Gaussian location mixture setting.

Finally in this subsection, we verify empirically the suboptimality of LCMLE
relative to ASM for larger sample sizes. We let $d=2$ and
$n\in\\{800,1600,3200,6400\\}$, and consider error distribution (iv) above as
well as the $t_{2}$ distribution (Example 26). We include the $t_{2}$
distribution in place of the standard Cauchy because when $P_{0}$ does not
have a finite first moment, its population level log-concave maximum
likelihood projection does not exist (and indeed we found empirically that
LCMLE was highly unstable in the Cauchy setting). Figure 10 displays the
average squared error loss of the oracle convex $M$-estimator, ASM, and LCMLE.
We see that the relative efficiency of ASM with respect to Oracle approaches 1
as $n$ increases, while LCMLE has an efficiency gap that does not vanish even
for large $n$, in agreement with the asymptotic calculation in Example 26.

###  4.2 Inference

In this subsection, we assess the finite-sample performance of the inferential
procedures for $\theta_{0}$ described in Section 3.3, but do not perform
sample splitting or cross-fitting. Here, we instead estimate the antitonic
information $i^{*}(p_{0})$ by

$\hat{\jmath}_{n}:=\frac{1}{n}\sum_{i=1}^{n}\hat{\psi}_{n}(\hat{\varepsilon}_{i})^{2},$

rather than $\hat{\imath}_{n}$, where
$\hat{\varepsilon}_{i}:=Y_{i}-\bar{\mu}_{n}-\tilde{X}_{i}^{\top}\bar{\theta}_{n}$
for $i\in[n]$ is the $i$th residual of the pilot estimator
$(\bar{\theta}_{n},\bar{\mu}_{n})$. We then compute
$\tilde{\mathcal{J}}_{n}:=\frac{\tilde{\jmath}_{n}}{n}\sum_{i=1}^{n}(\tilde{X}_{i}-\bar{X}_{n})(\tilde{X}_{i}-\bar{X}_{n})^{\top}$.

We first verify that the standardised error vector
$\sqrt{n}\tilde{\mathcal{J}}_{n}^{1/2}(\hat{\theta}_{n}^{\mathrm{ASM}}-\theta_{0})$
has a distribution very close to $N_{d-1}(0,I_{d-1})$ under three different
noise distributions: standard Gaussian, standard Cauchy, and the Gaussian
mixture $P_{0}=\frac{2}{3}N(0,1)+\frac{1}{3}N\bigl{(}\frac{1}{2},9\bigr{)}$.
Setting $n=600$ and $d=4$, we present the Q-Q plot for 8000 repetitions in
Figure 11. Moreover, Table 4(a) demonstrates that the average squared
estimation error $\bigl{(}\hat{\jmath}_{n}-i^{*}(p_{0})\bigr{)}^{2}$ is small.

Next, we show that the empirical coverage of the confidence ellipsoid

$\hat{C}_{n}^{\mathrm{ASM}}:=\bigl{\\{}v\in\mathbb{R}^{d-1}:n(\hat{\theta}_{n}^{\mathrm{ASM}}-v)^{\top}\tilde{\mathcal{J}}_{n}(\hat{\theta}_{n}^{\mathrm{ASM}}-v)\leq\chi^{2}_{d-1}(\alpha)\bigr{\\}}$

is close to the nominal level $1-\alpha$ when $n=600$. Indeed, for the
standard Gaussian and Gaussian mixture noise distributions above, Table 3
reveals that we maintain nominal coverage while for standard Cauchy errors,
our confidence set is slightly conservative. We also compare the volume of
$\hat{C}_{n}^{\mathrm{ASM}}$ with that of the OLS confidence ellipsoid

$\hat{C}_{n}^{\mathrm{OLS}}:=\bigl{\\{}v\in\mathbb{R}^{d-1}:n(\hat{\theta}_{n}^{\mathrm{OLS}}-v)^{\top}\hat{I}_{n}^{\mathrm{OLS}}(\hat{\theta}_{n}^{\mathrm{OLS}}-v)\leq\chi^{2}_{d-1}(\alpha)\bigr{\\}},$

where
$\hat{I}_{n}^{\mathrm{OLS}}:=n^{-1}\sum_{i=1}^{n}\hat{\sigma}^{-2}_{n}(\tilde{X}_{i}-\bar{X}_{n})(\tilde{X}_{i}-\bar{X}_{n})^{\top}$
and
$\hat{\sigma}_{n}^{2}:=n^{-1}\sum_{i=1}^{n}(Y_{i}-\hat{\mu}^{\mathrm{OLS}}_{n}-\tilde{X}_{i}^{\top}\hat{\theta}_{n}^{\mathrm{OLS}})^{2}$.
We plot examples of $\hat{C}_{n}^{\mathrm{ASM}}$ and
$\hat{C}_{n}^{\mathrm{OLS}}$ in Figure 12, and compute the average ratio of
their volumes in Table 4(b). The ASM ellipsoid is appreciably smaller than the
OLS ellipsoid except when the noise distribution is standard Gaussian, in
which case the OLS ellipsoid is slightly smaller.

##  5 Discussion

Despite the Gauss–Markov theorem, one of the messages of this paper is that
the success of ordinary least squares is relatively closely tied to Gaussian
or near-Gaussian error distributions. Our antitonic score matching approach
represents a middle ground that frees the practitioner from the Gaussian
straitjacket while retaining the convenience and stability of working with
convex loss functions. The Fisher divergence projection framework brings
together previously disparate ideas on shape-constrained estimation, score
matching, information theory and classical robust statistics. Given the
prevalence of procedures in statistics and machine learning that are
constructed as optimisers of pre-specified loss functions, we look forward to
seeing how related insights may lead to more flexible, data-driven approaches
that combine robustness and efficiency.

Acknowledgements: The authors thank Cun-Hui Zhang for helpful discussions. The
research of YCK and MX was supported by National Science Foundation grants
DMS-2311299 and DMS-2113671; OYF and RJS were supported by Engineering and
Physical Sciences Research Council Programme Grant EP/N031938/1, while RJS was
also supported by European Research Council Advanced Grant 101019498.

## References

  * Amari and Nagaoka, (2000) Amari, S.-i. and Nagaoka, H. (2000).  Methods of Information Geometry, volume 191.  American Mathematical Society.
  * Arcones, (1998) Arcones, M. A. (1998).  Asymptotic theory for $M$-estimators over a convex kernel.  Econometric Theory, 14(4):387–422.
  * Barber and Samworth, (2021) Barber, R. F. and Samworth, R. J. (2021).  Local continuity of log-concave projection, with applications to estimation under model misspecification.  Bernoulli, 27(4):2437–2472.
  * Barron, (2019) Barron, J. T. (2019).  A general and adaptive robust loss function.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 4331–4339.
  * Bean et al., (2013) Bean, D., Bickel, P. J., El Karoui, N., and Yu, B. (2013).  Optimal $M$-estimation in high-dimensional regression.  Proceedings of the National Academy of Sciences, 110(36):14563–14568.
  * Benton et al., (2024) Benton, J., Shi, Y., De Bortoli, V., Deligiannidis, G., and Doucet, A. (2024).  From denoising diffusions to denoising Markov models.  Journal of the Royal Statistical Society, Series B: Statistical Methodology, to appear.
  * Beran, (1978) Beran, R. (1978).  An efficient and robust adaptive estimator of location.  The Annals of Statistics, 6(2):292–313.
  * Betancourt et al., (2017) Betancourt, M., Byrne, S., Livingstone, S., and Girolami, M. (2017).  The geometric foundations of Hamiltonian Monte Carlo.  Bernoulli, 23(4A):2257–2298.
  * Bickel, (1975) Bickel, P. J. (1975).  One-step Huber estimates in the linear model.  Journal of the American Statistical Association, 70(350):428–434.
  * Bickel, (1982) Bickel, P. J. (1982).  On adaptive estimation.  The Annals of Statistics, 10(3):647–671.
  * Bobkov, (1996) Bobkov, S. G. (1996).  Extremal properties of half-spaces for log-concave distributions.  The Annals of Probability, 24(1):35–48.
  * Boyd and Vandenberghe, (2004) Boyd, S. P. and Vandenberghe, L. (2004).  Convex Optimization.  Cambridge University Press.
  * Brunel, (2023) Brunel, V.-E. (2023).  Geodesically convex $M$-estimation in metric spaces.  In The Thirty Sixth Annual Conference on Learning Theory, pages 2188–2210. PMLR.
  * Catoni, (2012) Catoni, O. (2012).  Challenging the empirical mean and empirical variance: a deviation study.  Annales de l’Institut Henri Poincaré – Probabilités et Statistiques, 48(4):1148–1185.
  * Chen and Samworth, (2013) Chen, Y. and Samworth, R. J. (2013).  Smoothed log-concave maximum likelihood estimation with applications.  Statistica Sinica, 23(3):1373–1398.
  * Cheng et al., (2018) Cheng, X., Chatterji, N. S., Bartlett, P. L., and Jordan, M. I. (2018).  Underdamped Langevin MCMC: a non-asymptotic analysis.  In Conference on Learning Theory, pages 300–323. PMLR.
  * Chernozhukov et al., (2018) Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W., and Robins, J. (2018).  Double/debiased machine learning for treatment and structural parameters.  The Econometrics Journal, 21(1):C1–C68.
  * Chinot et al., (2020) Chinot, G., Lecué, G., and Lerasle, M. (2020).  Robust statistical learning with Lipschitz and convex loss functions.  Probability Theory and Related Fields, 176(3):897–940.
  * Cook, (1977) Cook, R. D. (1977).  Detection of influential observation in linear regression.  Technometrics, 19(1):15–18.
  * Cover and Thomas, (2006) Cover, T. M. and Thomas, J. A. (2006).  Elements of Information Theory.  John Wiley & Sons, 2nd edition.
  * Cox, (1985) Cox, D. D. (1985).  A penalty method for nonparametric estimation of the logarithmic derivative of a density function.  Annals of the Institute of Statistical Mathematics, 37(2):271–288.
  * Cule et al., (2010) Cule, M., Samworth, R., and Stewart, M. (2010).  Maximum likelihood estimation of a multi-dimensional log-concave density.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 72(5):545–607.
  * Dalalyan et al., (2006) Dalalyan, A. S., Golubev, G. K., and Tsybakov, A. B. (2006).  Penalized maximum likelihood and semiparametric second-order efficiency.  The Annals of Statistics, 34(1):169–201.
  * De Bortoli et al., (2022) De Bortoli, V., Mathieu, E., Hutchinson, M., Thornton, J., Teh, Y. W., and Doucet, A. (2022).  Riemannian score-based generative modelling.  Advances in Neural Information Processing Systems, 35:2406–2422.
  * Derenski et al., (2023) Derenski, J., Fan, Y., James, G., and Xu, M. (2023).  An empirical Bayes shrinkage method for functional data.  Submitted.
  * Donoho and Montanari, (2015) Donoho, D. L. and Montanari, A. (2015).  Variance breakdown of Huber $M$-estimators: $n/p\in(1,\infty)$.  arXiv preprint arXiv:1503.02106.
  * Doss and Wellner, (2019) Doss, C. R. and Wellner, J. A. (2019).  Univariate log-concave density estimation with symmetry or modal constraints.  Electronic Journal of Statistics, 13(2):2391–2461.
  * Dümbgen et al., (2011) Dümbgen, L., Samworth, R., and Schuhmacher, D. (2011).  Approximation by log-concave distributions, with applications to regression.  The Annals of Statistics, 39(2):702–730.
  * Dümbgen et al., (2013) Dümbgen, L., Samworth, R. J., and Schuhmacher, D. (2013).  Stochastic search for semiparametric linear regression models.  In From Probability to Statistics and Back: High-Dimensional Models and Processes–A Festschrift in Honor of Jon A. Wellner, volume 9, pages 78–91. Institute of Mathematical Statistics.