# Functional Output Regression with Infimal Convolution:
Exploring the Huber and $\epsilon$-insensitive Losses

Alex Lambert  Dimitri Bouche  Zoltán Szabó  Florence d’Alché-Buc

###### Abstract

The focus of the paper is functional output regression (FOR) with convoluted
losses. While most existing work consider the square loss setting, we leverage
extensions of the Huber and the $\epsilon$-insensitive loss (induced by
infimal convolution) and propose a flexible framework capable of handling
various forms of outliers and sparsity in the FOR family. We derive
computationally tractable algorithms relying on duality to tackle the
resulting tasks in the context of vector-valued reproducing kernel Hilbert
spaces. The efficiency of the approach is demonstrated and contrasted with the
classical squared loss setting on both synthetic and real-world benchmarks.

kernel methods, functional data

problemequation problem

##  1 Introduction

Functional data analysis (FDA, Ramsay & Silverman 1997; Wang et al. 2016) has
attracted a growing attention in the field of machine learning and statistics,
with applications for instance in biomedical signal processing (Ullah & Finch,
2013), epidemiology monitoring and climate science (Ramsay & Silverman, 2007).
The key assumption is that we have access to densely-measured observations, in
which case functional data description becomes the most natural and adequate.
An important subfield of FDA is functional output regression (FOR) which
focuses on regression problems where the output variable is a function. There
are numerous ways to tackle the FOR problem family. The simplest approach is
to assume linear dependence between the inputs and the outputs (Morris, 2015).
However, in order to cope with more complex dependencies, various nonlinear
approaches have been designed. In nonparametric statistics, Ferraty et al.
(2011) proposed a Banach-valued Nadaraya-Watson estimator. The flexibility of
kernel methods (Steinwart & Christmann, 2008) and the richness of the
associated reproducing kernel Hilbert spaces (RKHSs; Micchelli et al. 2006)
have proven to be particularly useful in the area, with works involving tri-
variate regression problem (Reimherr et al., 2018), and approximated kernel
ridge regression (KRR) using orthonormal bases (Oliva et al., 2015). In the
operator-valued kernel (Pedrick, 1957; Carmeli et al., 2010) literature,
examples include function-valued KRR with double representer theorem (Lian,
2007), solvers based on the discretization of the loss function (Kadri et al.,
2010), purely functional methods relying on approximate inversion of integral
operators (Kadri et al., 2016) or techniques relying on finite-dimensional
coefficients of the functional outputs in a dictionary basis (Bouche et al.,
2021).

Most of these works employ the square loss which induces an estimate of the
conditional expectation of the functional outputs given the input data.
However defective sensors or malicious attacks can lead to erroneous or
contaminated measurements (Hubert et al., 2015), resulting in local or global
functional outliers. The square loss is expected to be badly affected in those
cases and considering alternative losses is a natural way to obtain reliable
and robust prediction systems. For scalar-valued outputs, the Huber loss
(Huber, 1964) and the $\epsilon$-insensitive loss (Lee et al., 2005) are
particularly popular and well-suited to construct outlier-robust estimators.
In the FDA setting, robustness has been investigated using Bayesian methods
(Zhu et al., 2011), trading the mean for the median (Cadre, 2001), using
bounded loss functions (Maronna & Yohai, 2013), or leveraging principal
component analysis (Kalogridis & Van Aelst, 2019).

In the operator-valued kernel literature, $\epsilon$-insensitive losses for
vector-valued regression have been proposed by Sangnier et al. (2017) for
finite-dimensional outputs. The use of convex optimization tools such as the
infimal convolution operator and parametric duality leads to efficient solvers
and provides sparse estimators. This idea is exploited by Laforgue et al.
(2020) where a generalization of this approach to infinite-dimensional outputs
encompassing both the Huber and $\epsilon$-insensitive losses is developed.

In this paper, we extend the families of losses considered by Laforgue et al.
(2020) by leveraging specific $p$-norms in functional spaces to handle various
forms of outliers (with Huber loss) and sparsity (with $\epsilon$-insensitive
losses). We study the properties of their Fenchel-Legendre conjugates, and
derive the associated dual optimization problems, which require suitable
representations and approximations adapted to each situation to be manageable
computationally. We propose tractable optimization algorithms for
$p\in\\{1,2\\}$ in the Huber loss scenario, and for $p\in\\{2,+\infty\\}$ with
the $\epsilon$-insensitive family. Finally, we provide an empirical study of
the proposed algorithms over synthetic and real functional datasets.

The paper is structured as follows. After introducing the general problem in
Section 2, we focus in Section 3 on a generalized family of Huber losses and
propose loss-specific tractable optimization schemes, before turning to the
family of $\epsilon$-insensitive losses in Section 4. We illustrate the
benefits of the approach on several benchmarks in Section 5. Proofs are
deferred to the supplement.

##  2 Problem Formulation

In this section, we introduce the general setting of FOR in the context of vv-
RKHSs, chosen for their modeling flexibility. To benefit from duality
principles, we focus on losses that can be expressed as infimal convolutions
in the functional output space.

##### Notations:

Let $\mathscr{X}$ be an input set, $\Theta\subset\mathbb{R}$ a compact set
endowed with a Borel probability measure $\mu$,
$\mathscr{Y}:=L^{2}[\Theta,\mu]$ the space of square $\mu$-integrable real-
valued functions. For $p\in[1,+\infty[$ and $f\in\mathscr{Y}$, let
$\left\lVert{f}\right\rVert_{p}=\left[\int_{\Theta}|f(\theta)|^{p}\mathrm{d}\mu(\theta)\right]^{\frac{1}{p}}\in[0,+\infty]$;
$\left\lVert{\cdot}\right\rVert_{\infty}$ refers to the essential supremum. In
both cases, the norm is allowed to take infinite value
($\left\lVert{\cdot}\right\rVert_{p}\colon\mathscr{Y}\to\mathbb{R}\cup\\{+\infty\\}$).111This
assumption is natural in convex optimization which is designed to handle
functions taking infinite values. Two numbers $p$ and $q\in[1,+\infty]$ are
said to be conjugate exponents if $\frac{1}{p}+\frac{1}{q}=1$, with the
classical $0=\frac{1}{\infty}$ convention. The ball in $\mathscr{Y}$ of radius
$\epsilon>0$ and center $0$ w.r.t. $\left\lVert{\cdot}\right\rVert_{p}$ is
denoted by $\mathcal{B}^{p}_{\epsilon}$. The space of bounded linear operators
over $\mathscr{Y}$ is $\mathscr{L}(\mathscr{Y})$. An operator-valued kernel
(OVK) $K\colon\mathscr{X}\times\mathscr{X}\to\mathscr{L}(\mathscr{Y})$ is a
mapping such that $\sum_{i=1}^{n}\sum_{j=1}^{n}\langle
K(x_{i},x_{j})y_{i},y_{j}\rangle_{\mathscr{Y}}\geq 0$ for all
$(x_{i},y_{i})_{i=1}^{n}\subset\mathscr{X}\times\mathscr{Y}$ and positive
integer $n$. An OVK $K$ gives rise to a space of functions from $\mathscr{X}$
to $\mathscr{Y}$ called vector-valued RKHS (vv-RKHS); it is defined as
$\mathscr{H}_{K}:=\overline{\operatorname{Span}}\\{K(\cdot,x)y\,:\,(x,y)\in\mathscr{X}\times\mathscr{Y}\\}$,
where $\operatorname{Span}(\cdot)$ denotes the linear hull of its argument,
$\bar{\cdot}$ stands for closure, and $K(\cdot,x)y$ is the function
$x^{\prime}\in\mathscr{X}\mapsto K(x^{\prime},x)y\in\mathscr{Y}$ while keeping
$x\in\mathscr{X}$ fixed. The _Fenchel-Legendre conjugate_ of a function
$f\colon\mathscr{Y}\to\left[-\infty,+\infty\right]$ is defined as
$f^{\star}(z):=\sup_{y\in\mathscr{Y}}~{}\langle z,y\rangle_{\mathscr{Y}}-f(y)$
where $z\in\mathscr{Y}$. Given a convex set $\mathcal{C}\subset\mathscr{Y}$,
$\iota_{\mathcal{C}}(\cdot)$ is its indicator function
($\iota_{\mathcal{C}}(x)=0$ if $x\in\mathcal{C}$, and
$\iota_{\mathcal{C}}(x)=+\infty$ otherwise), and
$\operatorname{Proj}_{\mathcal{C}}(\cdot)$ is the orthogonal projection on
$\mathcal{C}$ when $\mathcal{C}$ is also closed. Given two functions
$f,g\colon\mathscr{Y}\to\left]-\infty,+\infty\right]$, their _infimal
convolution_ is defined as
$f\operatorname{\square}g~{}(y)=\inf_{y^{\prime}\in\mathscr{Y}}f(y-y^{\prime})+g(y^{\prime})$
for $y\in\mathscr{Y}$ and the _proximal operator_ of $f$ (when $f$ is convex,
lower semi-continuous) is
$\operatorname{prox}_{f}(y)=\mathop{\mathrm{arg\,min}}_{y^{\prime}\in\mathscr{Y}}\frac{1}{2}\left\lVert{y-y^{\prime}}\right\rVert_{\mathscr{Y}}^{2}+f(y^{\prime})$
for all $y\in\mathscr{Y}$. For a positive integer $n$, let
$[n]=\\{1,\ldots,n\\}$. For $p\in[1,\infty[$, the $p$-norm of a vector
$\mathbf{v}\in\mathbb{R}^{m}$ is denoted by
$\left\lVert{\mathbf{v}}\right\rVert_{p}=\left(\sum_{j=1}^{m}|v_{j}|^{p}\right)^{\frac{1}{p}}$,
and $\left\lVert{\mathbf{v}}\right\rVert_{\infty}=\max_{j\in[m]}{|v_{j}|}$.
Given a matrix $\mathbf{A}\in\mathbb{R}^{n\times m}$ and $p,s\in[1,+\infty]$,
$\left\lVert{\mathbf{A}}\right\rVert_{p,s}$ is the $s$-norm of the $p$-norms
of the rows of $\mathbf{A}$. The positive part of $x\in\mathbb{R}$ is denoted
by $|x|_{+}=\max(x,0)$.

Next we introduce the FOR problem in vv-RKHSs. Recall that $\mathscr{X}$ is a
set and $\mathscr{Y}=L^{2}[\Theta,\mu]$, the latter capturing the functional
outputs. Assume that we have i.i.d. samples $(x_{i},y_{i})_{i\in[n]}$ from a
random variable $(\mathsf{X},\mathsf{Y})\in\mathscr{X}\times\mathscr{Y}$.
Given a proper, convex lower-semicontinuous loss function
$L\colon\mathscr{Y}\to\mathbb{R}$ and a regularization parameter $\lambda>0$,
we consider the _regularized empirical risk minimization_ problem

$\inf_{h\in\mathscr{H}_{K}}~{}\frac{1}{n}\sum_{i\in[n]}L(y_{i}-h(x_{i}))+\frac{\lambda}{2}\left\lVert{h}\right\rVert_{\mathscr{H}_{K}}^{2},$
(0)

where $K\colon\mathscr{X}\times\mathscr{X}\to\mathscr{L}(\mathscr{Y})$ is a
_decomposable_ OVK of the form $K=k_{\mathscr{X}}T_{k_{\Theta}}$. Here
$k_{\mathscr{X}}\colon\mathscr{X}\times\mathscr{X}\to\mathbb{R}$ and
$k_{\Theta}\colon\Theta\times\Theta\to\mathbb{R}$ are continuous real-valued
kernels, and $T_{k_{\Theta}}\in\mathscr{L}(\mathscr{Y})$ is the integral
operator associated to $k_{\Theta}$, defined for all $f\in\mathscr{Y}$ as
$(T_{k_{\Theta}}f)(\theta):=\int_{\Theta}k(\theta,\theta^{\prime})f(\theta)\mathrm{d}\mu(\theta^{\prime})$
where $\theta\in\Theta$. Similarly to the scalar case (Wahba, 1990), the
minimizer of (0) enjoys a representer theorem (Micchelli & Pontil, 2005) and
writes as

$\hat{h}=\frac{1}{\lambda
n}\sum_{i\in[n]}k_{\mathscr{X}}(\cdot,x_{i})T_{k_{\Theta}}\hat{\alpha}_{i}$

for some coefficients $\\{\hat{\alpha}_{i}\\}_{i\in[n]}\subset\mathscr{Y}$.
However, the functional nature of these parameters renders the problem
extremely challenging, with quite few existing solutions. Particularly, even
in the case of the square loss

$L(f)=\frac{1}{2}\left\lVert{f}\right\rVert_{\mathscr{Y}}^{2}=\frac{1}{2}\int_{\Theta}f(\theta)^{2}\mathrm{d}\mu(\theta),$

the value of $\\{\hat{\alpha}_{i}\\}_{i\in[n]}$ can not be computed in closed
form, and some level of approximation is required. For instance in Lian
(2007), $L$ is approximated as a discrete sum, allowing for a double
application of the representer theorem and yielding tractable models. In Kadri
et al. (2016), the integral operator is traded for a finite rank approximation
based on its eigendecomposition, providing a computable closed-form expression
for the coefficients. Aiming at robustness, Laforgue et al. (2020) propose a
Huber loss based on infimal convolution, yet limited to a narrow choice of
kernels $k_{\Theta}$. Moreover, the lack of flexibility in the definition of
the loss prevents the resulting estimators from being robust to a large
variety of outliers.

The goal of this work is (i) to widen the scope of the FOR problem ((0)) by
considering losses capable of handling different forms of outliers and
sparsity, and (ii) to design efficient optimization schemes for the resulting
tasks. The proposed two loss families are based on infimal convolution and can
be written in the form

$\displaystyle\frac{1}{2}\left\lVert{\cdot}\right\rVert_{\mathscr{Y}}^{2}\operatorname{\square}g$
(1)

for appropriately chosen functions $g\colon\mathscr{Y}\to]-\infty,+\infty]$.
The key property which allows one to handle these convoluted losses from an
optimization perspective is the fact that

$\left(\frac{1}{2}\left\lVert{\cdot}\right\rVert_{\mathscr{Y}}^{2}\operatorname{\square}g\right)^{\star}=\frac{1}{2}\left\lVert{\cdot}\right\rVert_{\mathscr{Y}}^{2}+g^{\star},$

as it makes the associated (0) amenable for dual approaches. The losses with
the proposed dedicated optimization schemes are detailed in Section 3 and
Section 4, respectively.

The starting point for working with convoluted losses in vv-RKHSs is the
following lemma.

######  Lemma 2.1 (Dualization for convoluted losses; Laforgue et al. 2020).

Let $L$ be a loss function defined as
$L=\frac{1}{2}\left\lVert{\cdot}\right\rVert_{\mathscr{Y}}^{2}\operatorname{\square}g$
for some $g\colon\mathscr{Y}\to]-\infty,+\infty]$. The solution of (0) is
given by

$\hat{h}=\frac{1}{\lambda
n}\sum_{i\in[n]}k_{\mathscr{X}}(\cdot,x_{i})T_{k_{\Theta}}\hat{\alpha}_{i},$
(2)

with $(\hat{\alpha}_{i})_{i\in[n]}\in\mathscr{Y}^{n}$ being the solution of
the dual task

$\displaystyle\inf_{(\alpha_{i})_{i\in[n]}\in\mathscr{Y}^{n}}$
$\displaystyle\sum_{i\in[n]}\left[\frac{1}{2}\left\lVert{\alpha_{i}}\right\rVert^{2}_{\mathscr{Y}}-\langle\alpha_{i},y_{i}\rangle_{\mathscr{Y}}+g^{\star}(\alpha_{i})\right]$
(2) $\displaystyle+\frac{1}{2\lambda
n}\sum_{i,j\in[n]}k_{\mathscr{X}}(x_{i},x_{j})\left\langle\alpha_{i},T_{k_{\Theta}}\alpha_{j}\right\rangle_{\mathscr{Y}}.$

Solving (2) in the general case for various $g$ and $k_{\Theta}$ raises
multiple challenges which have to be handled simultaneously. (2) is often
referred to as a _composite_ optimization problem, with a differentiable term
consisting of a quadratic part added to a non-differentiable term induced by
$g^{\star}$. The first challenge is to be able to compute the proximal
operator associated to $g^{\star}$. The second and third difficulties arise
from the fact that the dual variables are functions
($\alpha_{i}\in\mathscr{Y}$) and hence managing them computationally requires
specific care. Particularly, evaluation of
$\left\langle\alpha_{i},T_{k_{\Theta}}\alpha_{j}\right\rangle_{\mathscr{Y}}$
can be non-trivial. In addition, one has to design a finite-dimensional
description of the dual variables that is compatible with $T_{k_{\Theta}}$ and
the proximal operator of $g^{\star}$. The primary focus and technical
contribution of the paper is to handle these challenges, after which proximal
gradient descent optimization can be applied. We detail our proposed solution
in the next section.

##  3 Learning with $\mathscr{Y}$-Huber Losses

In this section, we propose a generalized Huber loss on $\mathscr{Y}$ based on
infimal convolution, followed by an efficient dual optimization approach to
solve the corresponding (0). This loss (as illustrated in Section 5) shows
robustness against different kind of outliers. Our proposed loss on
$\mathscr{Y}$ relies on functional $p$-norms where $p\in[1,+\infty]$.

######  Definition 3.1 ($\mathscr{Y}$-Huber loss).

Let $\kappa>0$ and $p\in[1,+\infty]$. We define the Huber loss with parameters
$(\kappa,p)$ as

$H_{\kappa}^{p}:=\frac{1}{2}\left\lVert{\cdot}\right\rVert_{\mathscr{Y}}^{2}\operatorname{\square}\kappa\left\lVert{\cdot}\right\rVert_{p}.$

Notice that in the specific case of $\mathscr{Y}=\mathbb{R}$, $H_{\kappa}^{p}$
reduces to the classical Huber loss on the real line for arbitrary $p$. Our
following result describes the behavior of $H_{\kappa}^{p}$.

######  Proposition 3.2.

Let $\kappa>0$, $p\in[1,+\infty]$, and $q$ the conjugate exponent of $p$. Then
for all $f\in\mathscr{Y}$,

$H_{\kappa}^{p}(f)=\frac{1}{2}\left\lVert{\operatorname{Proj}_{\mathcal{B}_{\kappa}^{q}}(f)}\right\rVert_{\mathscr{Y}}^{2}+\kappa\left\lVert{f-\operatorname{Proj}_{\mathcal{B}_{\kappa}^{q}}(f)}\right\rVert_{p}.$

Remark: For general $p$, the value of $H_{\kappa}^{p}(f)$ can not be computed
straighforwardly due to the complexity of the projection on
$\mathcal{B}_{\kappa}^{q}$. As we show however using a dual approach (0) is
still computationally manageable. For $p=2$, one gets back the loss
investigated by Laforgue et al. (2020).

The following proposition is a key result of this work that allows to leverage
$p$-norms as suitable candidates for $g$ in (1). It extends to $\mathscr{Y}$
the well-known finite-dimensional case that can e.g. be found in Bauschke et
al. (2011).

######  Proposition 3.3.

Let $p,q\in[1,+\infty]$ such that $\frac{1}{p}+\frac{1}{q}=1$. Then

$\left\lVert{\cdot}\right\rVert_{p}^{\star}=\iota_{\mathcal{B}^{q}_{1}}(\cdot).$

Our next result provides the dual of (0), and shows the impact of the
parameters $(p,\kappa)$.

######  Proposition 3.4 (Dual Huber).

Let $\kappa>0$, $p\in[1,+\infty]$, and $\frac{1}{p}+\frac{1}{q}=1$. The dual
of (0) with loss $H_{\kappa}^{p}$ writes

$\displaystyle\inf_{(\alpha_{i})_{i=1}^{n}\in\mathscr{Y}^{n}}$
$\displaystyle\sum_{i=1}^{n}\left[\frac{1}{2}\left\lVert{\alpha_{i}}\right\rVert_{\mathscr{Y}}^{2}-\langle\alpha_{i},y_{i}\rangle_{\mathscr{Y}}\right]$
(2) $\displaystyle+\frac{1}{2\lambda
n}\sum_{i=1}^{n}\sum_{j=1}^{n}k_{\mathscr{X}}(x_{i},x_{j})\left\langle\alpha_{i},T_{k_{\Theta}}\alpha_{j}\right\rangle_{\mathscr{Y}}$
$\displaystyle\mathrm{s.t}\,\left\lVert{\alpha_{i}}\right\rVert_{q}\leq\kappa,\quad
i\in[n].$

Remarks:

  * •

Influence of $\kappa$ and $p$: The difference between using the square loss
and the Huber loss $H_{\kappa}^{p}$ lies in the constraint on the $q$-norm of
the dual variables. The parameter $p$ influences the _shape_ of the ball via
the dual exponent $q$ ($\left\lVert{\alpha_{i}}\right\rVert_{q}$) defining the
admissible region for dual variables, and $\kappa$ determines its _size_
($\left\lVert{\alpha_{i}}\right\rVert_{q}\leq\kappa$). As $\kappa$ grows, the
constraint becomes void and we recover the solution of the classical ridge
regression problem. In Appendix C, we explore how different choices of $p$ can
affect the sensitivity of the loss to two different types of outliers.

  * •

Partially observed data: The observed data $(y_{i})_{i\in[n]}$ enter into (2)
only via their scalar product with the dual variables
$(\alpha_{i})_{i\in[n]}$. However in real life scenarii one never fully
observe the $y_{i}$ functions and these inner products are to be estimated.
One can instead assume access to a sampling at some locations
$(\theta_{j})_{j=1}^{m}$ which can be used to approximate the inner products
(see Section 3.1).

Let us now recall the challenges to be tackled to solve (2). Firstly, as
$\mathscr{Y}$ is infinite-dimensional no finite parameterization of the dual
variables can be assumed a priori. Secondly, even computing the different
terms of the objective function is non-trivial. Indeed, computing the
quadratic term corresponding to the regularization is not straightforward as
it involves the terms
$\left\langle\alpha_{i},T_{k_{\Theta}}\alpha_{j}\right\rangle_{\mathscr{Y}}$,
which require the knowledge of the action of the integral operator
$T_{k_{\Theta}}$. The third difficulty comes from handling the constraints.
Gradient-based optimization algorithms will require the projection of the dual
variables on the feasible set $\mathcal{B}_{\kappa}^{q}\subset\mathscr{Y}$,
which can be intractable to evaluate for some choices of $q$ depending on the
chosen representation.

The next proposition ensures the tractability of the projection step for
specific choices of $q$.

######  Proposition 3.5 (Projection on $\mathcal{B}_{\kappa}^{q}$).

Let $\kappa>0$. The projection on $\mathcal{B}_{\kappa}^{q}$ is tractable for
$q=2$ and $q=\infty$ and can be expressed for all
$(\alpha,\theta)\in\mathscr{Y}\times\Theta$ as

$\displaystyle\operatorname{Proj}_{\mathcal{B}_{\kappa}^{2}}(\alpha)$
$\displaystyle=\min{\left(1,\frac{\kappa}{\left\lVert{\alpha}\right\rVert_{2}}\right)}\alpha,$
(3)
$\displaystyle\left(\operatorname{Proj}_{\mathcal{B}_{\kappa}^{\infty}}(\alpha)\right)(\theta)$
$\displaystyle=\operatorname{sign}{(\alpha(\theta))}\min{(\kappa,|\alpha(\theta)|)}.$
(4)

The projection operator for $p=q=2$ in Equation 3 simply consists of a
multiplication by a scalar involving the $2$-norm of the dual variable. In the
$p=1$ case (i.e., $q=\infty$), the projection Equation 4 involves a
_pointwise_ projection. A suitable representation must guarantee the
feasibility of this projection, which requires a pointwise control over the
dual variables.

In order to solve (2), we propose to use two different representations. In
Section 3.1, we advocate representing the dual variables by linear splines and
approximating the action of $T_{k_{\Theta}}$ by Monte-Carlo (MC) sampling.
Splines allow pointwise control of the dual variable which is well-suited to
both $p=1$ and $p=2$. Our alternative approach (elaborated in Section 3.2)
relies on a finite-rank approximation of $T_{k_{\Theta}}$ using its
eigendecomposition. This method is applicable when $p=2$ with the
complementary advantage of performing dimensionality reduction.

Table 1: Correspondence between the quantities involved in (2) depending on the representation. | Linear splines | Eigenvectors of $T_{k_{\Theta}}$
---|---|---
$\sum_{i\in[n]}\frac{1}{2}\left\lVert{\alpha_{i}}\right\rVert_{\mathscr{Y}}^{2}$ | $\frac{1}{2m}\operatorname{Tr}\left(\mathbf{A}\mathbf{A}^{\top}\right)$ | $\operatorname{Tr}\left(\frac{1}{2}\mathbf{A}\mathbf{A}^{\top}\right)$
$\sum_{i\in[n]}\langle\alpha_{i},y_{i}\rangle_{\mathscr{Y}}$ | $\frac{1}{m}\operatorname{Tr}\left(\mathbf{A}\mathbf{Y}^{\top}\right)$ | $\operatorname{Tr}\left(\mathbf{A}\mathbf{R}^{\top}\right)$
$\sum_{i,j\in[n]}k_{\mathscr{X}}(x_{i},x_{j})\left\langle\alpha_{i},T_{k_{\Theta}}\alpha_{j}\right\rangle_{\mathscr{Y}}$ | $\frac{1}{m^{2}}\operatorname{Tr}\left(\mathbf{K}_{\mathscr{X}}\mathbf{A}\mathbf{K}_{\Theta}\mathbf{A}^{\top}\right)$ | $\operatorname{Tr}\left(\mathbf{K}_{\mathscr{X}}\mathbf{A}\mathbf{\Delta}\mathbf{A}^{\top}\right)$

###  3.1 The Linear Spline Based Approach

In this section we introduce a linear spline based representation for the dual
variables to tackle the challenges outlined.

Linear splines represent an easy-to-handle function class which provides
pointwise control over the dual variables as they are encoded by their
evaluations at some knots. While the class lacks smoothness, the dual
variables are smoothed out by $T_{k_{\Theta}}$ in the estimator expression
from Equation 2. Indeed, given that the kernel $k_{\Theta}$ is $2s$-times
continuously differentiable, the RKHS $\mathscr{H}_{k_{\Theta}}$ (where
$T_{k_{\Theta}}$ maps) consists of $s$-times continuously differentiable
functions (Zhou, 2008), a desirable property in many settings making linear
splines good candidates for modeling the dual variables.

A linear spline is a piecewise linear curve which can be encoded by a set of
ordered locations or anchor points $(\theta_{j})_{j\in[m]}\in\Theta^{m}$, and
by a vector of size $m$ corresponding to the evaluation of the spline at these
points. We choose the anchors to be distributed i.i.d. according to $\mu$; in
practice, we often take the locations to be those of the available sampling of
the observed data $(y_{i})_{i\in[n]}\in\mathscr{Y}^{n}$.

Fixing the anchors, the $n$ dual variables are encoded by the matrix of
evaluations
$\mathbf{A}=[\mathbf{a}_{i}]_{i\in[n]}=[\alpha_{i}(\theta_{j})]_{i\in[n],j\in[m]}\in\mathbb{R}^{n\times
m}$ with $\mathbf{a}_{i}$ being the $i^{th}$ row of $\mathbf{A}$. The action
of $T_{k_{\Theta}}$ on a function $\alpha\in\mathscr{Y}$ is then approximated
using MC approximation as

$T_{k_{\Theta}}\alpha\approx\frac{1}{m}\sum_{j\in[m]}k_{\Theta}(\cdot,\theta_{j})\alpha(\theta_{j}),$

resulting in the estimator

$\displaystyle h(x)(\theta)=\frac{1}{\lambda
nm}\sum_{i\in[n]}k_{\mathscr{X}}(x,x_{i})\sum_{j\in[m]}a_{ij}k_{\Theta}(\theta,\theta_{j}).$

Using this $\mathbf{A}$ parameterization of $h$, the different terms in (2)
are approximated as follows.

  * •

Squared norm of the dual variables: We approximate the squared norm of the
dual variables using MC sampling with locations $(\theta_{j})_{j\in[m]}$:

$\sum_{i\in[n]}\frac{1}{2}\left\lVert{\alpha_{i}}\right\rVert_{\mathscr{Y}}^{2}\approx\frac{1}{2m}\operatorname{Tr}\left(\mathbf{A}^{\top}\mathbf{A}\right).$

  * •

Scalar product with the data: We encode the evaluations of the observed
functions $(y_{i})_{i\in[n]}$ at locations $(\theta_{j})_{j\in[m]}$ in a
matrix
$\mathbf{Y}=[y_{i}(\theta_{j})]_{i\in[n],\,j\in[m]}\in\mathbb{R}^{n\times m}$,
and again use MC approximation

$\sum_{i\in[n]}\langle\alpha_{i},y_{i}\rangle_{\mathscr{Y}}\approx\frac{1}{m}\operatorname{Tr}\left(\mathbf{A}\mathbf{Y}^{\top}\right).$

  * •

Regularization term: Encoding the dual variables as linear splines hinders the
exact computation of the quadratic terms
$\langle\alpha_{i},T_{k_{\Theta}}\alpha_{j}\rangle_{\mathscr{Y}}$, which we
propose to approximate using a MC approximation of $T_{k_{\Theta}}$. Letting
$\mathbf{K}_{\mathscr{X}}\in\mathbb{R}^{n\times
n},\mathbf{K}_{\Theta}\in\mathbb{R}^{m\times m}$ be the Gram matrices
respectively associated to the data $(x_{i})_{i\in[n]}$ and kernel
$k_{\mathscr{X}}$, and to the locations $(\theta_{j})_{j\in[m]}$ and kernel
$k_{\Theta}$, the regularization term can be rephrased as $\frac{1}{\lambda
nm^{2}}\operatorname{Tr}\left(\mathbf{K}_{\mathscr{X}}\mathbf{A}\mathbf{K}_{\Theta}\mathbf{A}^{\top}\right)$.

  * •

Constraints: The constraints
$\left\lVert{\alpha_{i}}\right\rVert_{q}\leq\kappa$ can be handled similarly.
Particularly, let $\alpha_{i}\in\mathscr{S}_{m}$ with associated evaluations
vector $\mathbf{a}_{i}\in\mathbb{R}^{m}$ and $q\in[1,+\infty[$. We trade the
integral expression in the norm for an MC sum, resulting in the constraint
$\left\lVert{\mathbf{a}_{i}}\right\rVert_{q}\leq m^{\frac{1}{q}}\kappa$. This
expression also holds for $q=+\infty$, as
$\left\lVert{\alpha_{i}}\right\rVert_{\infty}\leq\kappa$ iff.
$\left\lVert{\mathbf{a}_{i}}\right\rVert_{\infty}\leq\kappa$.

Gathering the different terms (summarized in Table 1) yields the following
relaxation of (2):

$\displaystyle\inf_{\mathbf{A}\in\mathbb{R}^{n\times m}}$
$\displaystyle\operatorname{Tr}\left(\frac{1}{2}\mathbf{A}\mathbf{A}^{\top}-\mathbf{A}\mathbf{Y}^{\top}+\frac{1}{2\lambda
nm}\mathbf{K}_{\mathscr{X}}\mathbf{A}\mathbf{K}_{\Theta}\mathbf{A}^{\top}\right)$
(4) $\displaystyle{s.t.~{}}\left\lVert{\mathbf{A}}\right\rVert_{q,\infty}\leq
m^{\frac{1}{q}}\kappa.$

Remark: The decomposable assumption on the kernel $K$ plays a role in the
regularization. It has the effect of disentangling the action of both Gram
matrices $\mathbf{K}_{\mathscr{X}}$ and $\mathbf{K}_{\Theta}$.

We propose to tackle (4) using accelerated proximal gradient descent (APGD),
where the proximal step amounts to projecting the coefficients on the $q$-ball
of radius $m^{1/q}\kappa$. The technique is summarized in Algorithm 1. The
gradient stepsize $\gamma$ can be computed exactly from the parameters of the
problem. Indeed, for guaranteed convergence, one must set $\gamma<\frac{2}{C}$
where $C$ is the Lipschitz constant associated to the gradient of the
objective function; here $C\leq 1+\frac{1}{\lambda
n}\left\lVert{\mathbf{K_{\mathscr{X}}}}\right\rVert_{\mathrm{op}}\left\lVert{\mathbf{K_{\Theta}}}\right\rVert_{\mathrm{op}}.$
The initialization can either be the null matrix
$\mathbf{A}^{(0)}=\mathbf{0}_{\mathbb{R}^{n\times m}}$ or the solution of the
unconstrained optimization problem which can be obtained in closed form. This
solution (i) can dramatically reduce the number of iterations for small
$\epsilon$ or large $\kappa$, and (ii) can be computed in
$\mathcal{O}(n^{3}+m^{3})$ time exploiting the Kronecker structure inherited
from the separable kernel with a Sylvester equation solver (Sima, 1996;
Dinuzzo et al., 2011). Since the objective function in (4) is the sum of two
functions, one convex and differentiable with Lipschitz continuous gradient
(the quadratic form) and one convex and lower semi-continuous (the indicator
function of the constraint set), the optimal worst case complexity is
$\mathcal{O}\left(\frac{1}{T^{2}}\right)$ (Beck & Teboulle, 2009). The time
complexity per iteration is dominated by the computation of the matrix
$\mathbf{K}_{\mathscr{X}}\mathbf{V}\mathbf{K}_{\Theta}$ which is
$\mathcal{O}(n^{2}m+m^{2}n)$.

input : Gram matrices $\mathbf{K}_{\mathscr{X}},\mathbf{K}_{\Theta}$, data
matrix $\mathbf{Y}$, regularization parameter $\lambda$, loss parameters
$(\kappa,p)$ or $(\epsilon,p)$, gradient step $\gamma$

init : $\mathbf{A}^{(0)},\mathbf{A}^{(-1)}=\bm{0}\in\mathbb{R}^{n\times m}$

for  _epoch $t$ from $0$ to $T-1$_ do

// gradient step

$\mathbf{V}=\mathbf{A}^{(t)}+\frac{t-2}{t+1}\left(\mathbf{A}^{(t)}-\mathbf{A}^{(t-1)}\right)$

$\mathbf{U}=\mathbf{V}-\gamma\left(\mathbf{V}+\frac{1}{\lambda
nm}\mathbf{K}_{\mathscr{X}}\mathbf{V}\mathbf{K}_{\Theta}-\mathbf{Y}\right)$

// projection step

if  _$p=2$_ then

for  _row $i\in[n]$_ do

$\mathbf{a}_{i}^{(t+1)}=\min\left(\frac{\sqrt{m}\kappa}{\left\lVert{\mathbf{u}_{i}}\right\rVert_{2}},1\right)\mathbf{u}_{i}$
// if $H_{\kappa}^{2}$

$\mathbf{a}_{i}^{(t+1)}=\left|1-\frac{\gamma\epsilon}{\sqrt{m}\left\lVert{\mathbf{u}_{i}}\right\rVert_{2}}\right|_{+}\mathbf{u}_{i}$
// if $\ell_{\epsilon}^{2}$

else

for  _row $i\in[n]$_ do

for  _column $j\in[m]$_ do

$a_{ij}^{(t+1)}=\operatorname{sign}(u_{ij})\min\left(\kappa,\left|u_{ij}\right|\right)$
// if $H_{\kappa}^{1}$

$a_{ij}^{(t+1)}=\operatorname{sign}{(u_{ij})}\left||u_{ij}|-\frac{\gamma\epsilon}{m}\right|_{+}$
// if $\ell_{\epsilon}^{\infty}$

return _$\mathbf{A}^{(T)}$_

Algorithm 1 APGD with linear splines

###  3.2 The Eigendecomposition Approach

In this section, we propose an alternative finite-dimensional description of
the dual variables relying on an approximate eigendecomposition of
$T_{k_{\Theta}}$ when $p=2$. The rationale of this approach is to decrease the
number of parameters needed to represent the estimator by selecting directions
well-suited to $T_{k_{\Theta}}$, namely the dominant $r$ eigenvectors of
$T_{k_{\Theta}}$. As computing the eigendecomposition of $T_{k_{\Theta}}$ is
generally intractable, we propose an approximation detailed in the following.

Let us consider the problem of finding a continuous eigenvector $\psi$ of a
sampled version of the integral operator $T_{k_{\Theta}}$ with eigenvalue
$\delta>0$ (Hoegaerts et al., 2005):

$\frac{1}{m}\sum_{j\in[m]}k_{\Theta}(\theta,\theta_{j})\psi(\theta_{j})=\delta\psi(\theta)\text{
for }\forall\theta\in\Theta.$

By evaluating it at points $(\theta_{j})_{j\in[m]}$, one gets that $(\delta
m,\psi(\theta_{j})_{j\in[m]})$ form the eigensystem of the Gram matrix
$\mathbf{K}_{\Theta}\in\mathbb{R}^{m\times m}$. These can be computed using
for instance singular value decomposition, and by substitution one arrives at
an approximated eigenbasis of dimension at most $m$; this can be used as a
proxy instead of the true eigenvectors of $T_{k_{\Theta}}$. By using the first
$r$ of these vectors for some $r<m$ we are able to lower the size of the
parameterization of the model. We store in a diagonal matrix
$\bm{\Delta}\in\mathbb{R}^{r\times r}$ the first $r$ eigenvalues.

The problem is now parameterized by a matrix
$\mathbf{A}=[\mathbf{a}_{i}]_{i\in[n]}\in\mathbb{R}^{n\times r}$ with each row
$\mathbf{a}_{i}\in\mathbb{R}^{r}$ encoding the coefficients of the dual
variable $\alpha_{i}$ on the $(\psi_{j})_{j\in[r]}$ basis. The estimator then
reads as

$h(x)(\theta)=\frac{1}{\lambda
n}\sum_{i\in[n],\,j\in[r]}a_{ij}\delta_{j}k_{\mathscr{X}}(x,x_{i})\psi_{j}(\theta).$

We store in $\mathbf{R}=[R_{ij}]_{i\in[n],j\in[r]}\in\mathbb{R}^{n\times r}$
the scalar products between the observed data and the eigenbasis:
$R_{ij}=\langle y_{i},\psi_{j}\rangle_{\mathscr{Y}}$. The correspondence
between the different optimization terms are summarized in Table 1; the
optimization reduces to

$\displaystyle\inf_{\mathbf{A}\in\mathbb{R}^{n\times r}}~{}~{}$
$\displaystyle\operatorname{Tr}\left(\frac{1}{2}\mathbf{A}\mathbf{A}^{\top}-\mathbf{A}\mathbf{R}^{\top}+\frac{1}{2\lambda
n}\mathbf{K}_{\mathscr{X}}\mathbf{A}\mathbf{\Delta}\mathbf{A}^{\top}\right)$
$\displaystyle{s.t.~{}}\left\lVert{\mathbf{A}}\right\rVert_{2,\infty}\leq\kappa.$
(5)

Again, one can use APGD to solve this task; the resulting computations are
deferred to Algorithm 2 in the supplement.

Remark: For $\kappa$ large enough, the solution reduces to that of Kadri et
al. (2016). However, our approximation allows handling a wide range of kernels
in contrast to the analytical knowledge imposed for the eigensystem of
$T_{k_{\Theta}}$ by Kadri et al. (2016); Laforgue et al. (2020).

##  4 Learning with $\epsilon$-insensitive Losses

In this section, we propose a generalized $\epsilon$-insensitive version of
the square loss on $\mathscr{Y}$ involving infimal convolution, and derive
tractable dual optimization algorithm to solve (0). This loss induces sparsity
on the matrix of coefficients as illustrated in Section 5.

######  Definition 4.1 ($\epsilon$-insensitive loss).

Let $\epsilon>0$ and $p\in[1,+\infty]$. We define the $\epsilon$-insensitive
version of the square loss with parameters $(\epsilon,p)$ as

$\ell_{\epsilon}^{p}:=\frac{1}{2}\left\lVert{\cdot}\right\rVert_{\mathscr{Y}}^{2}\operatorname{\square}\iota_{\mathcal{B}^{p}_{\epsilon}}(\cdot).$

When $\mathscr{Y}=\mathbb{R}$, the loss reduces to the classical
$\epsilon$-insensitive version of the square loss regardless of $p$. The
following proposition (counterpart of Proposition 3.2) sheds light on the
effect of the infimal convolution on the square loss.

######  Proposition 4.2.

Let $\epsilon>0$ and $p\in[1,+\infty]$. Then for all $f\in\mathscr{Y}$,

$\ell_{\epsilon}^{p}(f)=\frac{1}{2}\left\lVert{f-\operatorname{Proj}_{\mathcal{B}_{\epsilon}^{p}}(f)}\right\rVert^{2}_{\mathscr{Y}}.$
(6)

Remark: Proposition 4.2 means that $\ell_{\epsilon}^{p}(f)=0$ when
$\left\lVert{f}\right\rVert_{p}\leq\epsilon$, i.e. small residuals do not
contribute to the risk. For general $p$, $\ell_{\epsilon}^{p}(f)$ is not
straightforward to compute due to the complexity of
$\operatorname{Proj}_{\mathcal{B}_{\epsilon}^{p}}(f)$. As we however use a
dual approach, (0) can still be tackled computationally.

The next result shows how to dualize (0) when the proposed
$\epsilon$-insensitive loss is used.

######  Proposition 4.3 (Dual $\epsilon$-insensitive).

Let $\epsilon\geq 0,p\in[1,+\infty]$, and $\frac{1}{p}+\frac{1}{q}=1$. The
dual of (0) writes as

$\displaystyle\inf_{(\alpha_{i})_{i\in[n]}\in\mathscr{Y}^{n}}$
$\displaystyle\sum_{i\in[n]}\left[\frac{1}{2}\left\lVert{\alpha_{i}}\right\rVert_{\mathscr{Y}}^{2}-\langle\alpha_{i},y_{i}\rangle_{\mathscr{Y}}+\epsilon\left\lVert{\alpha_{i}}\right\rVert_{q}\right]$
(6) $\displaystyle+\frac{1}{2\lambda
n}\sum_{i\in[n],\,j\in[n]}k_{\mathscr{X}}(x_{i},x_{j})\left\langle\alpha_{i},T_{k_{\Theta}}\alpha_{j}\right\rangle_{\mathscr{Y}}.$

Remark (influence of $\epsilon$ and $p$): Compared to the square loss,
$\ell_{\epsilon}^{p}$ induces an additional term
$\sum_{i\in[n]}\epsilon\left\lVert{\alpha_{i}}\right\rVert_{q}$ in the dual.
Setting $\epsilon=0$ recovers the square loss case.

The challenges involving the representation of the dual variables, and the
computability of the different terms composing (6) are similar to those evoked
in Section 3. We have however traded the constraints on the $q$-norms of the
dual variables against an additional non-smooth term. As for the Huber loss
family in Section 3, we address this convex non-smooth optimization problem
through the APGD algorithm. The proximal step involves the computation of
$\operatorname{prox}_{\gamma\epsilon\left\lVert{\cdot}\right\rVert_{q}}$ for a
suitable gradient stepsize $\gamma>0$, which is the focus of the next
proposition.

######  Proposition 4.4 (Proximal $q$-norm).

Let $\epsilon>0$. The proximal operator of
$\epsilon\left\lVert{\cdot}\right\rVert_{q}$ is computable for $q=1$ and
$q=2$, and given for all $(\alpha,\theta)\in\mathscr{Y}\times\Theta$ by

$\displaystyle\left(\operatorname{prox}_{\epsilon\left\lVert{\cdot}\right\rVert_{1}}(\alpha)\right)(\theta)$
$\displaystyle=\operatorname{sign}{(\alpha(\theta))}\left|\left|\alpha(\theta)\right|-\epsilon\right|_{+},$
(7)
$\displaystyle\operatorname{prox}_{\epsilon\left\lVert{\cdot}\right\rVert_{2}}(\alpha)$
$\displaystyle=\alpha\left|1-\frac{\epsilon}{\left\lVert{\alpha}\right\rVert_{\mathscr{Y}}}\right|_{+}.$
(8)

We recognize in Equation 7 a pointwise soft thresholding, and in Equation 8 an
analogous to the block soft thresholding, both are known to promote sparsity.

To solve (6) we rely on the two kinds of finite representations introduced
previously. In Section 4.1, we tackle the $p=2$ and $p=\infty$ case with the
linear splines based method from Section 3.1, before using dimensionality
reduction from Section 3.2 for the case $p=2$ in Section 4.2.

###  4.1 The Linear Spline Based Approach

Similarly to what was presented in Section 3.1, we use linear splines to
represent the dual variables $(\alpha_{i})_{i\in[n]}$ as they allow a
pointwise control over the dual variable and thus give rise to computable
proximal operators. Keeping the notations, the optimization boils down to

$\displaystyle\inf_{\mathbf{A}\in\mathbb{R}^{n\times m}}$
$\displaystyle\operatorname{Tr}\left(\frac{1}{2}\mathbf{A}\mathbf{A}^{\top}-\mathbf{A}\mathbf{Y}^{\top}+\frac{1}{2\lambda
nm}\mathbf{K}_{\mathscr{X}}\mathbf{A}\mathbf{K}_{\Theta}\mathbf{A}^{\top}\right)$
$\displaystyle+\frac{\epsilon}{m^{\frac{1}{q}}}\sum_{i\in[n]}\left\lVert{\mathbf{a}_{i}}\right\rVert_{q}.$

We use APGD to solve it with steps detailed in Algorithm 1. When $q=1$, the
proximal operator is the soft thresholding operator, akin to promote sparsity
in the dual coefficients.

###  4.2 The Eigendecomposition Approach

We mobilize the eigendecomposition technique from Section 3.2 to solve (6) in
the case $p=2$. Using the same notation as in Problem 3.2, we get the
following task:

$\displaystyle\inf_{\mathbf{A}\in\mathbb{R}^{n\times
r}}\operatorname{Tr}\left(\frac{1}{2}\mathbf{A}\mathbf{A}^{\top}-\mathbf{A}\mathbf{R}^{\top}+\frac{1}{2\lambda
n}\mathbf{K}_{\mathscr{X}}\mathbf{A}\mathbf{\Delta}\mathbf{A}^{\top}\right)+\epsilon\sum_{i\in[n]}\left\lVert{\mathbf{a}_{i}}\right\rVert_{2}.$

APGD is applied to tackle this problem; the details are deferred to Algorithm
2 in the supplement. Notice that the proximal operator in this case is the
block soft thresholding operator, known to promote structured row-wise
sparsity.

##  5 Numerical Experiments

In this section, we demonstrate the efficiency of the proposed convoluted
losses. The implementation is done in Python, and is available in the form of
an open source package at https://github.com/allambert/foreg.
The experiments are centered around two key directions:

  1. 1.

The first goal is to understand the accuracy-sparsity tradeoffs of the
$\epsilon$-insensitive loss as a function of the regularization $\lambda$ and
insensitivity parameter $\epsilon$.

  2. 2.

Our second aim is to quantify the robustness of the Huber losses w.r.t.
different forms of outliers with a particular focus on global versus local
ones. To gain further insight into the robustness w.r.t. corruption, we
designed 3 types of functional outliers with distinct characteristics.

Our proposed losses are investigated on 3 benchmarks: a synthetic one
associated to Gaussian processes, followed by two real-world ones arising in
the context of neuroimaging and speech analysis. We investigate both questions
on the synthetic data, and provide further insights for the first and the
second question on the neuroimaging and the speech dataset, respectively.

We now detail the 3 functional outlier types used in our experiments on
robustness to study the effect of local and global corruption. Local outliers
affect the functions only on small portions of $\Theta$ whereas global ones
contaminates them in their entirety. To corrupt the functions
$(y_{i})_{i\in[n]}$, we first draw a set $I\subset\\{0,\ldots,n\\}$ of size
$\lfloor\tau n\rfloor$ corresponding to the indices to contaminate;
$\tau\in[0,1]$ being the proportion of contaminated functions. Then, we
perform different kinds of corruption:
$\bullet$ Type 1: Let $\omega$ be the permutation defined for $j\in[|I|]$ as
$\omega(I_{j})=I_{j+1}$ if $j<|I|$ and $\omega(I_{|I|})=I_{1}$, then for $i\in
I$, the data point $(x_{i},y_{i})$ is replaced by $(x_{i},-y_{\omega(i)})$.
$\bullet$ Type 2: Given covariance parameters
$\boldsymbol{\sigma}\in\mathbb{R}^{r}$ and an intensity parameter $\zeta>0$,
we draw a Gaussian process $g_{c}\sim\mathcal{GP}(0,k_{\sigma_{c}})$ for
$c\in[r]$ where $k_{\sigma_{c}}$ is the Gaussian covariance function with
standard deviation $\sigma_{c}$. Then, for $i\in I$, we replace
$(x_{i},y_{i})$ with $\left(x_{i},\sum_{c\in[r]}a_{ic}g_{c}\right)$ where the
coefficients $a_{ic}$ are drawn i.i.d. from a uniform distribution
$\mathcal{U}([-0.5\zeta,0.5\zeta])$.
$\bullet$ Type 3: For each $i\in I$, a randomly chosen fraction $\xi\in[0,1]$
of the discrete observations for $y_{i}$ is replaced by random draws from a
uniform distribution $\mathcal{U}([-b_{\text{max}},b_{\text{max}}])$, where
$b_{\text{max}}:=\underset{i,j}{\max}|y_{i}(\theta_{j})|$. The corruptions of
Type 1 and 2 are global ones whereas Type 3 is a local one. In terms of the
characteristics of the different corruptions, for Type 1 the properties of the
outlier functions remain close to those of the non-outlier ones, whereas with
Type 2 they become completely different. Finally, for corrupted data in the
hyperparameter choice using cross-validation the mean was replaced with
median.

For the losses $H_{\kappa}^{1}$ and $\ell_{\epsilon}^{\infty}$ we solve the
problem based on the representation with linear splines (see Section 3.1 and
Section 4.1 respectively); this is the only possible approach. However for the
losses $H_{\kappa}^{2}$ and $\ell_{\epsilon}^{2}$ we exploit the
representation using a truncated basis of approximate eigenfunctions (see
Section 3.2 and Section 4.2 respectively), in doing so we reduce the
computational cost compared to the linear splines approach. Concerning
optimization, we deployed the APGD method (Beck & Teboulle, 2009) with
backtracking line search, and adaptive restart (O’Donoghue & Candès, 2015).
The initialization in APGD was carried out with the closed-form solution
available for the square loss using a Sylvester solver.

Table 2: MSEs and sparsity on the DTI dataset $\lambda$ | Metric | $\nicefrac{{1}}{{2}}\lVert\cdot\rVert_{\mathscr{Y}}^{2}$ | $H_{\kappa}^{2}$ | $H_{\kappa}^{1}$ | $\ell^{2}_{\epsilon}$ | $\ell_{\epsilon}^{\infty}$
---|---|---|---|---|---|---
10-5 |  MSE (10-1) |  2.5$\pm$0.19 |  2.21$\pm$0.31 |  2.21$\pm$0.31 |  2.41$\pm$0.26 |  2.5$\pm$0.23
| Sparsity | - | - | - |  27.4$\pm$17.2% |  85.9$\pm$10.7%
10-3 |  MSE (10-1) |  2.18$\pm$0.27 |  2.23$\pm$0.32 |  2.21$\pm$0.32 |  2.2$\pm$0.29 |  2.18$\pm$0.28
| Sparsity | - | - | - |  3.4$\pm$6.9% |  12.7$\pm$10.5%

Regarding the performance measure applied for evaluation, let
$((y_{i}(\theta_{ij}))_{j\in[m_{i}]})_{i\in[n]}$ be the set of observed
discretized functions and let
$(\widehat{y}_{i}(\theta_{ij}))_{j\in[m_{i}]})_{i\in[n]}$ be an estimated set
of discretized functions, where $(\theta_{ij})_{j\in[m_{i}]}$ denotes the
observation locations for $y_{i}$. We used the mean squared error defined as

$\text{MSE}:=\frac{1}{n}\sum_{i\in[n]}\sum_{j\in[m_{i}]}[y_{i}(\theta_{ij})-\widehat{y}_{i}(\theta_{ij})]^{2}.$

When $m_{i}=m$ for all $i$, we normalize it by $m$ and define
NMSE$:=\frac{1}{m}$ MSE.

###  5.1 Experiments on the synthetic dataset

The impact of the different losses are investigated in detail on a function-
to-function synthetic dataset whose construction is detailed in Section B.2.
The kernels $k_{\mathscr{X}}$ and $k_{\Theta}$ are chosen to be Gaussian and
the experiments are averaged over 20 draws with training and testing samples
of size 100.

Figure 1: Interaction between regularization $\lambda$ and insensitivity
$\epsilon$ for the loss $\ell_{\epsilon}^{2}$ (1st row) and
$\ell_{\epsilon}^{\infty}$ (2nd row).

$\epsilon$-insensitive loss: To study the interaction between $\lambda$ and
$\epsilon$ and the resulting sparsity-accuracy trade-offs, we added i.i.d.
Gaussian noise with standard deviation $0.5$ to the observations of the output
functions. The resulting MSE values are summarized in Fig. 1. For both the
$\ell_{\epsilon}^{2}$ and the $\ell_{\epsilon}^{\infty}$ loss, one can reduce
$\lambda$, increase $\epsilon$ and get a fair amount of sparsity while making
a small compromise in terms of accuracy.

Huber loss: We investigate the robustness of the Huber loss to different types
of outliers while selecting both $\lambda$ and $\kappa$ through robust cross
validation. The resulting MSE values are summarized in Fig. 2. As it can be
seen in the first row of the figure, the losses $H_{\kappa}^{1}$ and
$H_{\kappa}^{2}$ are significantly more robust to global outliers than the
square loss, both when the outliers’ intensity $\zeta$ and the proportion
$\tau$ of contaminated samples vary. The second row of the figure shows that
when dealing with local outliers, the closer one gets to the whole sample
being contaminated ($\tau=1$), the less robust $H_{\kappa}^{2}$ becomes. On
the other hand, $H_{\kappa}^{1}$ shows remarkable robustness as it can be
observed at the bottom right panel (in case of $\tau=1$). One can interpret
this phenomenon by noticing that the loss $H_{\kappa}^{2}$ can penalize less
big discrepancy between functions in the $\|.\|_{\mathscr{Y}}$ norm sense, but
if all samples are contaminated locally a little, the outliers are meddled in
the norm and so $H_{\kappa}^{2}$ becomes inefficient.

Figure 2: Robustness to different types of outliers.

###  5.2 Experiments on the DTI dataset

Table 3: MSEs on speech data VT | $\nicefrac{{1}}{{2}}\lVert\cdot\rVert_{\mathscr{Y}}^{2}$ | $H_{\kappa}^{2}$ | $H_{\kappa}^{1}$ | $\ell^{2}_{\epsilon}$ | $\ell_{\epsilon}^{\infty}$
---|---|---|---|---|---
LP |  6.58$\pm$0.62 |  6.59$\pm$0.62 |  6.59$\pm$0.64 |  6.58$\pm$0.62 |  6.58$\pm$0.62
LA |  4.65$\pm$0.55 |  4.65$\pm$0.55 |  4.66$\pm$0.55 |  4.64$\pm$0.55 |  4.64$\pm$0.55
TBCL |  4.26$\pm$0.46 |  4.26$\pm$0.46 |  4.27$\pm$0.46 |  4.26$\pm$0.46 |  4.26$\pm$0.46
TBCD |  4.67$\pm$0.37 |  4.68$\pm$0.38 |  4.7$\pm$0.38 |  4.67$\pm$0.38 |  4.67$\pm$0.38
VEL |  2.94$\pm$0.5 |  2.94$\pm$0.5 |  2.95$\pm$0.5 |  2.94$\pm$0.5 |  2.94$\pm$0.5
GLO |  7.25$\pm$0.65 |  7.26$\pm$0.65 |  7.25$\pm$0.64 |  7.25$\pm$0.65 |  7.25$\pm$0.65
TTCL |  3.76$\pm$0.21 |  3.76$\pm$0.21 |  3.74$\pm$0.2 |  3.73$\pm$0.21 |  3.73$\pm$0.21
TTCD |  5.93$\pm$0.34 |  5.94$\pm$0.34 |  5.93$\pm$0.35 |  5.92$\pm$0.34 |  5.92$\pm$0.34
Table 4: MSEs on contaminated speech data VT |  Type 1 Outliers ($\tau=0.1$) |  Type 3 outliers ($\tau=0.1$, $\xi=0.1$)
---|---|---
$\nicefrac{{1}}{{2}}\lVert\cdot\rVert_{\mathscr{Y}}^{2}$ | $H_{\kappa}^{2}$ | $H_{\kappa}^{1}$ | $\nicefrac{{1}}{{2}}\lVert\cdot\rVert_{\mathscr{Y}}^{2}$ | $H_{\kappa}^{2}$ | $H_{\kappa}^{1}$
LP |  9.4$\pm$0.75 |  9.4$\pm$0.66 |  9.19$\pm$0.79 |  7.53$\pm$0.58 |  7.62$\pm$0.59 |  7.0$\pm$ 0.59
LA |  5.72$\pm$0.76 |  5.63$\pm$0.71 |  5.52$\pm$0.69 |  5.06$\pm$0.6 |  5.11$\pm$0.6 |  5.09$\pm$0.55
TBCL |  6.71$\pm$0.96 |  6.14$\pm$0.97 |  5.98$\pm$0.93 |  5.06$\pm$0.51 |  5.16$\pm$0.48 |  4.72$\pm$ 0.54
TBCD |  5.8$\pm$0.41 |  5.86$\pm$0.44 |  5.83$\pm$0.44 |  5.18$\pm$0.4 |  5.26$\pm$0.41 |  5.08$\pm$ 0.4
VEL |  4.37$\pm$0.56 |  3.76$\pm$0.62 |  3.76$\pm$0.59 |  3.52$\pm$0.57 |  3.54$\pm$0.58 |  3.41$\pm$ 0.57
GLO |  9.61$\pm$0.87 |  9.51$\pm$0.86 |  9.53$\pm$0.84 |  7.94$\pm$0.61 |  8.02$\pm$0.61 |  7.76$\pm$ 0.61
TTCL |  15.06$\pm$2.22 |  9.51$\pm$0.63 |  9.48$\pm$0.6 |  5.89$\pm$0.43 |  5.91$\pm$0.45 |  6.62$\pm$0.66
TTCD |  8.15$\pm$0.48 |  7.96$\pm$0.49 |  8.02$\pm$0.51 |  6.63$\pm$0.44 |  6.74$\pm$0.42 |  6.36$\pm$ 0.39

In our next experiment we considered the DTI benchmark222This dataset was
collected at the Johns Hopkins University and the Kennedy-Krieger Institute..
The dataset contains a collection of fractional anisotropy profiles deduced
from diffusion tensor imaging scans, and we take the first scans of the
$n=100$ multiple sclerosis patients. The profiles are given along two tracts,
the corpus callosum and the right cortiospinal. The goal is to predict the
latter function from the former, which can be framed as a function-to-function
regression problem. When some functions admit missing observations, we fill in
the gaps by linear interpolation, and later use the MSE as metric. We use a
Gaussian kernel for $k_{\mathscr{X}}$ and a Laplacian one for $k_{\Theta}$ and
average over 10 runs with $n_{\text{train}}=70$ and $n_{\text{test}}=30$.

Similarly to our experiences gained on the synthetic dataset, a compromise can
be made between the two parameters $\lambda$ and $\epsilon$ to get increased
sparsity, as can be observed in Table 2. Moreover, we highlight that even for
optimal regularization with respect to the square loss $\lambda=10^{-3}$, one
gets a fair amount of sparsity while getting the same score with
$\ell_{\epsilon}^{\infty}$ and a very small difference with
$\ell_{\epsilon}^{2}$.

###  5.3 Speech data

In this section, we focus on a speech inversion problem (Mitra et al., 2009).
Particularly, our goal is to predict a vocal tract (VT) configuration that
likely produced a speech signal (Richmond, 2002). This benchmark encompasses
$n=413$ synthetically pronounced words to which 8 VT functions are associated:
LA, LP, TTCD, TTCL, TBCD, TBCL, VEL, GLO. This is then a time-
series–to–function regression problem. We predict the VT functions separately
in eight subproblems.

Since the words are of varying length, we use the MSE as metric and extend
symmetrically the signals to match the longest word for in training. We encode
the input sounds through 13 mel-frequency cepstral coefficients (MFCC) and
normalize the VT functions’ values to the range [-1, 1]. We average over $10$
train-test splits taking $n_{\text{train}}=250$ and $n_{\text{test}}=163$.
Finally we take an integral Gaussian kernel on the standardized MFCCs (see
Appendix B for further details) as $k_{\mathscr{X}}$ and a Laplace kernel as
$k_{\Theta}$.

We first compare all the losses on untainted data in Table 3. Then to evaluate
the robustness of the Huber losses, we ran experiments on contaminated data
with two configurations. In the first case, we added Type 1 (global) outliers
with $\tau=0.1$ and in the second one, we added Type 3 (local) outliers with
$\tau=0.1$ and $\xi=0.1$. The results are displayed in Table 4. In the
contaminated setting, one gets results similar to ones obtained on the
synthetic dataset. The loss $H_{\kappa}^{1}$ works especially well for local
outliers whereas the loss $H_{\kappa}^{2}$ is robust only to global outliers.

##  6 Conclusion

In this paper we introduced generalized families of loss functions based on
infimal convolution and p-norms for functional output regression. The
resulting optimization problems were handled using duality principles. Future
work could focus on extending these techniques to a wider choice of $p$ using
iterative techniques for the proximal steps.

## Acknowledgements

AL, DB, and FdB received funding from the _Télécom Paris Research and Testing
Chair on Data Science and Artificial Intelligence for Digitalized Industry and
Service_. AL also obtained additional funding from the _ERC Advanced Grant
E-DUALITY (787960)_.

## References

  * Bauschke et al. (2011) Bauschke, H. H., Combettes, P. L., et al.  _Convex analysis and monotone operator theory in Hilbert spaces_.  Springer, 2011.
  * Beck & Teboulle (2009) Beck, A. and Teboulle, M.  A fast iterative shrinkage-thresholding algorithm for linear inverse problems.  _SIAM Journal on Imaging Sciences_ , 2(1):183–202, 2009.
  * Bouche et al. (2021) Bouche, D., Clausel, M., Roueff, F., and d’Alché Buc, F.  Nonlinear functional output regression: A dictionary approach.  In _International Conference on Artificial Intelligence and Statistics (AISTATS)_ , pp. 235–243, 2021.
  * Cadre (2001) Cadre, B.  Convergent estimators for the l1-median of Banach valued random variable.  _Statistics: A Journal of Theoretical and Applied Statistics_ , 35(4):509–521, 2001.
  * Carmeli et al. (2010) Carmeli, C., De Vito, E., Toigo, A., and Umanitá, V.  Vector valued reproducing kernel Hilbert spaces and universality.  _Analysis and Applications_ , 8(01):19–61, 2010\.
  * Dinuzzo et al. (2011) Dinuzzo, F., Ong, C. S., Gehler, P., and Pillonetto, G.  Learning output kernels with block coordinate descent.  In _International Conference on Machine Learning (ICML)_ , pp. 49–56, 2011.
  * Ferraty et al. (2011) Ferraty, F., Laksaci, A., Tadj, A., Vieu, P., et al.  Kernel regression with functional response.  _Electronic Journal of Statistics_ , 5:159–171, 2011.
  * Hoegaerts et al. (2005) Hoegaerts, L., Suykens, J. A., Vandewalle, J., and De Moor, B.  Subset based least squares subspace regression in RKHS.  _Neurocomputing_ , 63:293–323, 2005.
  * Huber (1964) Huber, P. J.  Robust estimation of a location parameter.  _The Annals of Mathematical Statistics_ , pp. 73–101, 1964.
  * Hubert et al. (2015) Hubert, M., Rousseeuw, P. J., and Segaert, P.  Multivariate functional outlier detection.  _Statistical Methods and Applications_ , 24:177–202, 2015\.
  * Kadri et al. (2010) Kadri, H., Duflos, E., Preux, P., Canu, S., and Davy, M.  Nonlinear functional regression: a functional RKHS approach.  In _International Conference on Artificial Intelligence and Statistics (AISTATS)_ , pp. 374–380, 2010.
  * Kadri et al. (2016) Kadri, H., Duflos, E., Preux, P., Canu, S., Rakotomamonjy, A., and Audiffren, J.  Operator-valued kernels for learning from functional response data.  _Journal of Machine Learning Research_ , 17(20):1–54, 2016.
  * Kalogridis & Van Aelst (2019) Kalogridis, I. and Van Aelst, S.  Robust functional regression based on principal components.  _Journal of Multivariate Analysis_ , 173:393–415, 2019\.
  * Laforgue et al. (2020) Laforgue, P., Lambert, A., Brogat-Motte, L., and d’Alché Buc, F.  Duality in RKHSs with infinite dimensional outputs: Application to robust losses.  In _International Conference on Machine Learning (ICML)_ , pp. 5598–5607, 2020.
  * Lee et al. (2005) Lee, Y.-J., Hsieh, W.-F., and Huang, C.-M.  epsilon-SSVR: A smooth support vector machine for epsilon-insensitive regression.  _IEEE Transactions on Knowledge & Data Engineering_, (5):678–685, 2005.
  * Lian (2007) Lian, H.  Nonlinear functional models for functional responses in reproducing kernel Hilbert spaces.  _Canadian Journal of Statistics_ , 35(4):597–606, 2007.
  * Maronna & Yohai (2013) Maronna, R. A. and Yohai, V. J.  Robust functional linear regression based on splines.  _Computational Statistics & Data Analysis_, 65:46–55, 2013.
  * Micchelli & Pontil (2005) Micchelli, C. and Pontil, M.  On learning vector-valued functions.  _Neural Computation_ , 17:177–204, 2005.
  * Micchelli et al. (2006) Micchelli, C., Xu, Y., and Zhang, H.  Universal kernels.  _Journal of Machine Learning Research_ , 7:2651–2667, 2006\.
  * Mitra et al. (2009) Mitra, V., Ozbek, Y., Nam, H., Zhou, X., and Espy-Wilson, C. Y.  From acoustics to vocal tract time functions.  In _IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , pp. 4497–4500, 2009.
  * Moreau (1965) Moreau, J. J.  Proximité et dualité dans un espace hilbertien.  Technical report, 1965.  (https://hal.archives-ouvertes.fr/hal-01740635).
  * Morris (2015) Morris, J. S.  Functional regression.  _Annual Review of Statistics and Its Application_ , 2:321–359, 2015.
  * Oliva et al. (2015) Oliva, J., Neiswanger, W., Póczos, B., Xing, E., Trac, H., Ho, S., and Schneider, J.  Fast function to function regression.  In _International Conference on Artificial Intelligence and Statistics (AISTATS)_ , pp. 717–725, 2015.
  * O’Donoghue & Candès (2015) O’Donoghue, B. and Candès, E.  Adaptive restart for accelerated gradient schemes.  _Foundations of Computational Mathematics_ , 15:715–732, 2015.
  * Parikh & Boyd (2014) Parikh, N. and Boyd, S.  Proximal algorithms.  _Foundations and Trends in Optimization_ , 1(3):127–239, 2014.
  * Pedrick (1957) Pedrick, G.  Theory of reproducing kernels for Hilbert spaces of vector-valued functions.  Technical report, University of Kansas, Department of Mathematics, 1957\.
  * Ramsay & Silverman (1997) Ramsay, J. and Silverman, B.  Functional data analysis, 1997.
  * Ramsay & Silverman (2007) Ramsay, J. O. and Silverman, B. W.  _Applied functional data analysis: methods and case studies_.  Springer, 2007.
  * Reimherr et al. (2018) Reimherr, M., Sriperumbudur, B., Taoufik, B., et al.  Optimal prediction for additive function-on-function regression.  _Electronic Journal of Statistics_ , 12(2):4571–4601, 2018.
  * Richmond (2002) Richmond, K.  _Estimating Articulatory Parameters from the Acoustic Speech Signal_.  PhD thesis, The Center for Speech Technology Research, Edinburgh University, 2002.
  * Sangnier et al. (2017) Sangnier, M., Fercoq, O., and d’Alché-Buc, F.  Data sparse nonparametric regression with $\epsilon$-insensitive losses.  In _Asian Conference on Machine Learning (ACML)_ , pp. 192–207, 2017.
  * Sima (1996) Sima, V.  _Algorithms for Linear-Quadratic Optimization_.  Chapman and Hall/CRC, 1996.
  * Steinwart & Christmann (2008) Steinwart, I. and Christmann, A.  _Support vector machines_.  Springer Science & Business Media, 2008.
  * Ullah & Finch (2013) Ullah, S. and Finch, C. F.  Applications of functional data analysis: A systematic review.  _BMC medical research methodology_ , 13(1):1–12, 2013.
  * Wahba (1990) Wahba, G.  _Spline Models for Observational Data_.  SIAM, CBMS-NSF Regional Conference Series in Applied Mathematics, 1990\.
  * Wang et al. (2016) Wang, J.-L., Chiou, J.-M., and Müller, H.-G.  Functional data analysis.  _Annual Review of Statistics and Its Application_ , 3:257–295, 2016.
  * Zhou (2008) Zhou, D.-X.  Derivative reproducing properties for kernel methods in learning theory.  _Journal of Computational and Applied Mathematics_ , 220:456–463, 2008.
  * Zhu et al. (2011) Zhu, H., Brown, P. J., and Morris, J. S.  Robust, adaptive functional regression in functional mixed model framework.  _Journal of the American Statistical Association_ , 106(495):1167–1179, 2011.

The supplement is structured as follows. We present the proofs of our results
in Appendix A. Appendix B complements the main part of the paper by providing
additional algorithmic and experimental details. Finally, Appendix C includes
additional plots and insights about the loss functions.

##  Appendix A Proofs

In this section, we present the proofs of our results. In Section A.1 we
recall some definitions from convex optimization used throughout the proofs,
with focus on Fenchel-Legendre conjugation and proximal operators, followed by
the proofs themselves (Section A.2-A.8).

###  A.1 Reminder on Convex Optimization

Recall that $\mathscr{Y}:=L^{2}[\Theta,\mu]$ where $\Theta\subset\mathbb{R}$
is a compact set endowed with a probability measure $\mu$.

######  Definition A.1 (Proper, convex, lower semi-continuous functions).

We denote by $\Gamma_{0}(\mathscr{Y})$ the set of functions
$J\colon\mathscr{Y}\to\left]-\infty,+\infty\right]$ that are

  1. 1.

_proper_ :
$\operatorname{dom}(J):=\left\\{f\in\mathscr{Y}\,:\,J(f)<+\infty\right\\}\neq\varnothing$,

  2. 2.

_convex_ : $J(tf+(1-t)g)\leq tJ(f)+(1-t)J(g)$ for $\forall
f,g\in\mathscr{Y},\forall t\in[0,1]$, and

  3. 3.

_lower semicontinuous_ : $\varliminf_{g\to f}J(g)\geq J(f)$ for $\forall
f\in\mathscr{Y}$, where $\varliminf$ denotes limit inferior.

######  Definition A.2.

The _Fenchel-Legendre conjugate_ of a function
$J\colon\mathscr{Y}\to\left[-\infty,+\infty\right]$ is defined as

$\displaystyle J^{\star}(f)$ $\displaystyle:=\sup_{g\in\mathscr{Y}}~{}\langle
f,g\rangle_{\mathscr{Y}}-J(g),\quad f\in\mathscr{Y}.$

The Fenchel-Legendre conjugate of a function $J$ is always convex. It is also
involutive on $\Gamma_{0}(\mathscr{Y})$, meaning that $(J^{\star})^{\star}=J$
for any $J\in\Gamma_{0}(\mathscr{Y})$. We gather in Table 5 examples and
properties of Fenchel-Legendre conjugates.

We now introduce the infimal convolution operator following Bauschke et al.
(2011).

######  Definition A.3 (Infimal convolution).

The _infimal convolution_ of two functions
$L,J\colon\mathscr{Y}\to\left]-\infty,+\infty\right]$ is

$L\operatorname{\square}J\colon\left(\begin{aligned}
\mathscr{Y}&\to[-\infty,+\infty]\\\
f&\mapsto\inf_{g\in\mathscr{Y}}L(f-g)+J(g)\end{aligned}\right).$

One key property of the infimal convolution operator is that it behaves nicely
under Fenchel-Legendre conjugation, as it is detailed in the following
proposition.

######  Lemma A.4 (Bauschke et al. 2011, Proposition 13.24).

Let $L,J\colon\mathscr{Y}\to\left]-\infty,+\infty\right]$. Then

$\displaystyle(L\operatorname{\square}J)^{\star}=L^{\star}+J^{\star}.$

We now define the proximal operator, used as a replacement for the classical
gradient step in the presence of non-differentiable objective functions.

######  Definition A.5 (Proximal operator, Moreau 1965).

The _proximal operator_ (or proximal map) is defined as

$\displaystyle\operatorname{prox}_{J}(f)$
$\displaystyle:=\mathop{\mathrm{arg\,min}}_{g\in\mathscr{Y}}~{}J(g)+\frac{1}{2}\left\lVert{f-g}\right\rVert^{2}_{\mathscr{Y}},\text{
for }(J,f)\in\Gamma_{0}(\mathscr{Y})\times\mathscr{Y}.$ (9)

One advantage of working with functions in $\Gamma_{0}(\mathscr{Y})$ is that
the proximal operator is always well-defined. Its computation is doable for
various losses thanks to the following lemma.

######  Lemma A.6 (Moreau decomposition, Moreau 1965).

Let $J\in\Gamma_{0}(\mathscr{Y})$ and $\gamma>0$. Then

$\operatorname{Id}=\operatorname{prox}_{\gamma
J}(\cdot)+\gamma\operatorname{prox}_{J^{\star}/\gamma}(\cdot/\gamma),$ (10)

where $\operatorname{Id}$ stands for the identity operator.

Table 5: Properties of Fenchel-Legendre conjugate, for any $J,L\colon\mathscr{Y}\to\left[-\infty,+\infty\right]$ and $p,q\in[1,+\infty]$ s.t. $\frac{1}{p}+\frac{1}{q}=1$.  Function | Fenchel-Legendre conjugate
---|---
$\frac{1}{2}\left\lVert{\cdot}\right\rVert_{\mathscr{Y}}^{2}$ | $\frac{1}{2}\left\lVert{\cdot}\right\rVert_{\mathscr{Y}}^{2}$
$\left\lVert{\cdot}\right\rVert_{p}$ | $\iota_{\mathcal{B}^{q}_{1}}$
$\epsilon J$ |  $\epsilon J^{\star}(\frac{\cdot}{\epsilon})$ for all $\epsilon>0$
$J(\cdot-y)$ |  $J^{\star}+\langle\cdot,y\rangle_{\mathscr{Y}}$ for all $y\in\mathscr{Y}$
$L\operatorname{\square}J$ | $L^{\star}+J^{\star}$

We remind the reader that we want to solve

$\inf_{h\in\mathscr{H}_{K}}~{}\frac{1}{n}\sum_{i\in[n]}L(y_{i}-h(x_{i}))+\frac{\lambda}{2}\left\lVert{h}\right\rVert_{\mathscr{H}_{K}}^{2},$
(10)

where $K\colon\mathscr{X}\times\mathscr{X}\to\mathscr{L}(\mathscr{Y})$ is a
_decomposable_ OVK of the form $K=k_{\mathscr{X}}T_{k_{\Theta}}$. Here
$k_{\mathscr{X}}\colon\mathscr{X}\times\mathscr{X}\to\mathbb{R}$ and
$k_{\Theta}\colon\Theta\times\Theta\to\mathbb{R}$ are continuous real-valued
kernels, and $T_{k_{\Theta}}\in\mathscr{L}(\mathscr{Y})$ is the integral
operator associated to $k_{\Theta}$, defined for all $\alpha\in\mathscr{Y}$ by

$(T_{k_{\Theta}}\alpha)(\theta)=\int_{\Theta}k(\theta,\theta^{\prime})\alpha(\theta)\mathrm{d}\mu(\theta^{\prime})\text{
for all }\theta\in\Theta.$

We also remind the reader to the dual of (10) when the loss writes as an
infimal convolution.

######  Lemma A.7 (Dualization for convoluted losses, Laforgue et al. 2020).

Let $L$ be a loss function defined as
$L=\frac{1}{2}\left\lVert{\cdot}\right\rVert_{\mathscr{Y}}^{2}\operatorname{\square}g$
for some $g\colon\mathscr{Y}\to]-\infty,+\infty]$. Then the solution of (10)
is given by

$\hat{h}=\frac{1}{\lambda
n}\sum_{i\in[n]}k_{\mathscr{X}}(\cdot,x_{i})T_{k_{\Theta}}\hat{\alpha}_{i},$
(11)

with $(\hat{\alpha}_{i})_{i\in[n]}\in\mathscr{Y}^{n}$ being the solution of
the dual task

$\displaystyle\inf_{(\alpha_{i})_{i\in[n]}\in\mathscr{Y}^{n}}~{}\sum_{i\in[n]}\left[\frac{1}{2}\left\lVert{\alpha_{i}}\right\rVert^{2}_{\mathscr{Y}}-\langle\alpha_{i},y_{i}\rangle_{\mathscr{Y}}+g^{\star}(\alpha_{i})\right]+\frac{1}{2\lambda
n}\sum_{i,j\in[n]}k_{\mathscr{X}}(x_{i},x_{j})\left\langle\alpha_{i},T_{k_{\Theta}}\alpha_{j}\right\rangle_{\mathscr{Y}}.$
(11)

###  A.2 Proof of Proposition 3.3

Before going through the proof, let us recall Hölder’s inequality.

######  Lemma A.8 (Hölder’s inequality).

Let $p,q\in[1,+\infty]$ be conjugate exponents, in other words
$\frac{1}{p}+\frac{1}{q}=1$. Let $\Theta$ be a measurable space enriched with
probability measure $\mu$. Then for any $f,g:\Theta\rightarrow\mathbb{R}$
measurable functions one has

$\int_{\Theta}\left\lvert
f(\theta)g(\theta)\right\rvert\mathrm{d}\mu(\theta)\leq\left\lVert{f}\right\rVert_{p}\left\lVert{g}\right\rVert_{q}.$

Moreover, if $p\in]1,+\infty[$, $f\in L^{p}[\Theta,\mu]$ and $g\in
L^{q}[\Theta,\mu]$, then equality is attained if and only if $\left\lvert
f\right\rvert^{p}$ and $\left\lvert g\right\rvert^{p}$ are linearly dependent
in $L^{1}[\Theta,\mu]$.

We now introduce a lemma useful to the proof of Proposition 3.3.

######  Lemma A.9.

Let $p,q\in]1,+\infty[$ be conjugate exponents and $f\in\mathscr{Y}$ such that
$1<\left\lVert{f}\right\rVert_{q}<+\infty$. Then there exist $h\in\mathscr{Y}$
and $C>0$ such that

$\langle f,h\rangle_{\mathscr{Y}}-\left\lVert{h}\right\rVert_{p}\geq C.$

Moreover, one can choose $h$ such that whenever $f(\theta)=0$, $h(\theta)=0$
also holds.

##### Proof of Lemma A.9

Let $p,q\in]1,+\infty[$ be conjugate exponents and $f\in\mathscr{Y}$ such that
$1<\left\lVert{f}\right\rVert_{q}<+\infty$. We know that Hölder’s inequality
becomes an equality if and only if $\left\lvert f\right\rvert^{q}$ and
$\left\lvert g\right\rvert^{p}$ are linearly dependent in $L^{1}[\Theta,\mu]$.
To that end, let $g\colon\Theta\to\mathbb{R}$ be defined as

$\displaystyle g(\theta)=\operatorname{sign}(f(\theta))\left\lvert
f(\theta)\right\rvert^{\frac{q}{p}}\quad\text{where $\theta\in\Theta$}.$ (12)

It is to be noted that $g$ does not necessarily belong to $\mathscr{Y}$, yet
it belongs to $L^{p}[\Theta,\mu]$. By construction, we have

$\int_{\Theta}f(\theta)g(\theta)\mathrm{d}\mu(\theta)=\left\lVert{f}\right\rVert_{q}\left\lVert{g}\right\rVert_{p}.$
(13)

We consider a sequence $(g_{n})_{n\in\mathbb{N}}\in\mathscr{Y}^{\mathbb{N}}$
such that $g_{n}(\theta)=\operatorname{sign}{(g(\theta))}\min{(\left\lvert
g(\theta)\right\rvert,n)}$ with $(n,\theta)\in\mathbb{N}\times\Theta$. As
$\left\lvert g_{n}(\theta)\right\rvert\leq n$ for all
$(n,\theta)\in\mathbb{N}\times\Theta$ and $\mu$ is a probability measure, the
functions $g_{n}$ belong to $\mathscr{Y}$. Since (i)
$g_{n}(\theta)\xrightarrow{n\to\infty}g(\theta)$ for all $\theta\in\Theta$ and
(ii) $\left\lvert g_{n}(\theta)\right\rvert\leq\left\lvert
g(\theta)\right\rvert$ for any $n\in\mathbb{N}$ holds $\mu$-almost everywhere,
the dominated convergence theorem in $L^{p}[\Theta,\mu]$ ensures that
$\left\lVert{g-g_{n}}\right\rVert_{p}\xrightarrow{n\to\infty}0$. Consequently,
it holds that for all $n\in\mathbb{N}$,

$\displaystyle\left\lvert\int_{\Theta}f(\theta)g(\theta)\mathrm{d}\mu(\theta)-\int_{\Theta}f(\theta)g_{n}(\theta)\mathrm{d}\mu(\theta)\right\rvert$
$\displaystyle=\left\lvert\int_{\Theta}f(\theta)[g(\theta)-g_{n}(\theta)]\mathrm{d}\mu(\theta)\right\rvert\stackrel{{\scriptstyle(a)}}{{\leq}}\int_{\Theta}\left\lvert
f(\theta)\right\rvert\left\lvert
g(\theta)-g_{n}(\theta)\right\rvert\mathrm{d}\mu(\theta)$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}\left\lVert{f}\right\rVert_{q}\left\lVert{g-g_{n}}\right\rVert_{p}.$

In (a) we used that the absolute value of the integral can be upper bounded by
the integral of the absolute value, in (b) the Hölder’s inequality was
invoked. Thus by
$\left\lVert{g-g_{n}}\right\rVert_{p}\xrightarrow{n\to\infty}0$ and
$\left\lVert{f}\right\rVert_{q}<+\infty$, this means that $\langle
f,g_{n}\rangle_{\mathscr{Y}}\xrightarrow{n\to\infty}\int_{\Theta}f(\theta)g(\theta)\mathrm{d}\mu(\theta)\stackrel{{\scriptstyle\eqref{eq:Holder-
eq}}}{{=}}\left\lVert{f}\right\rVert_{q}\left\lVert{g}\right\rVert_{p}$, and
for all $\epsilon>0$, there exist $N\in\mathbb{N}$ such that for all $n\geq
N$, $\langle
f,g_{n}\rangle_{\mathscr{Y}}>(\left\lVert{f}\right\rVert_{q}-\epsilon)\left\lVert{g}\right\rVert_{p}$.
In particular for $\epsilon=\frac{\left\lVert{f}\right\rVert_{q}-1}{2}>0$, we
have $\langle
f,g_{N}\rangle_{\mathscr{Y}}>\frac{1+\left\lVert{f}\right\rVert_{q}}{2}\left\lVert{g}\right\rVert_{p}$.
Then,

$\displaystyle\langle
f,g_{N}\rangle_{\mathscr{Y}}-\left\lVert{g_{N}}\right\rVert_{p}$
$\displaystyle\stackrel{{\scriptstyle(c)}}{{\geq}}\langle
f,g_{N}\rangle_{\mathscr{Y}}-\left\lVert{g}\right\rVert_{p}\stackrel{{\scriptstyle(d)}}{{\geq}}\frac{1+\left\lVert{f}\right\rVert_{q}}{2}\left\lVert{g}\right\rVert_{p}-\left\lVert{g}\right\rVert_{p}\geq\underbrace{\frac{\left\lVert{f}\right\rVert_{q}-1}{2}\left\lVert{g}\right\rVert_{p}}_{>0}.$

In (c) we used that
$\left\lVert{g_{N}}\right\rVert_{p}\leq\left\lVert{g}\right\rVert_{p}$, (d) is
implied by $\langle
f,g_{N}\rangle_{\mathscr{Y}}>\frac{1+\left\lVert{f}\right\rVert_{q}}{2}\left\lVert{g}\right\rVert_{p}$.
Taking $h=g_{N}$ and
$C=\frac{\left\lVert{f}\right\rVert_{q}-1}{2}\left\lVert{g}\right\rVert_{p}$
yields the announced result, by noticing that (12) shows that $f(\theta)=0$
also implies $h(\theta)=g_{N}(\theta)=g(\theta)=0$. ∎

We are now ready to prove Proposition 3.3, which is the building block for
dualizing optimization problems resulting from the generalized Huber and
$\epsilon$-insensitive losses whose definition can be found respectively in
Definition 3.1 and Definition 4.1. The proposition is an extension of the
well-studied finite-dimensional case to the space $\mathscr{Y}$.

######  Proposition (3.3).

Let $p,q\in[1,+\infty]$ such that $\frac{1}{p}+\frac{1}{q}=1$. Then

$\left\lVert{\cdot}\right\rVert_{p}^{\star}=\iota_{\mathcal{B}^{q}_{1}}(\cdot).$
(14)

##### Proof

The proof is structured as follows. We first consider the case of $p=1$,
followed by $p\in]1,+\infty]$, and $p=+\infty$. The reasoning in all cases
rely heavily on Hölder’s inequality. Throughout the proof it is assumed that
$f\in\mathscr{Y}$.

##### Case $p=1$:

The reasoning goes as follows: we show that
$\left\lVert{f}\right\rVert_{\infty}\leq 1$ implies
$\left\lVert{\cdot}\right\rVert_{1}^{\star}(f)=0$, and
$\left\lVert{f}\right\rVert_{\infty}>1$ gives
$\left\lVert{\cdot}\right\rVert_{1}^{\star}(f)=+\infty$, which allow one to
conclude that
$\left\lVert{\cdot}\right\rVert_{1}^{\star}=\iota_{\mathcal{B}^{\infty}_{1}}(\cdot)$.

  * •

When $\left\lVert{f}\right\rVert_{\infty}\leq 1$: Exploiting Hölder’s
inequality, it holds that

$\langle
f,g\rangle_{\mathscr{Y}}\leq\left\lVert{f}\right\rVert_{\infty}\left\lVert{g}\right\rVert_{1}\quad\text{
for all }g\in\mathscr{Y}.$

Since $\left\lVert{f}\right\rVert_{\infty}\leq 1$, this implies that

$\langle f,g\rangle_{\mathscr{Y}}-\left\lVert{g}\right\rVert_{1}\leq
0\quad\text{ for all }g\in\mathscr{Y}.$

The supremum being attained for $g=0$, we conclude that
$\left\lVert{\cdot}\right\rVert_{1}^{\star}(f)=0$.

  * •

When $\left\lVert{f}\right\rVert_{\infty}>1$: Let
$A=\left\\{\theta\in\Theta\,:\,\left\lvert
f(\theta)\right\rvert>\frac{1+\left\lVert{f}\right\rVert_{\infty}}{2}\right\\}$.
By the definition of the essential supremum, $\mu(A)>0$. We define
$g\colon\Theta\to\mathbb{R}$ to be the function:
$g(\theta)=\operatorname{sign}{(f(\theta))}$ if $\theta\in A$ and $0$
otherwise. Since $g$ is bounded, $g\in\mathscr{Y}$. Denoting by $t>0$ a
running parameter, it holds that

$\displaystyle\langle
f,tg\rangle_{\mathscr{Y}}-\left\lVert{tg}\right\rVert_{1}$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}\langle
f,tg\rangle_{\mathscr{Y}}-t\mu(A)=t\int_{\Theta}f(\theta)g(\theta)\mathrm{d}\mu(\theta)-t\mu(A)\stackrel{{\scriptstyle(a)}}{{=}}t\int_{A}\left\lvert
f(\theta)\right\rvert\mathrm{d}\mu(\theta)-t\mu(A)$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\geq}}t\mu(A)\frac{1+\left\lVert{f}\right\rVert_{\infty}}{2}-t\mu(A)=t\underbrace{\mu(A)\frac{\left\lVert{f}\right\rVert_{\infty}-1}{2}}_{>0}\xrightarrow{t\to\infty}+\infty.$

In (a) we used the definition of $g$, (b) is implied by the fact that
$\left\lvert
f(\theta)\right\rvert>\frac{1+\left\lVert{f}\right\rVert_{\infty}}{2}$ for all
$\theta\in A$. Thus $\left\lVert{\cdot}\right\rVert_{1}^{\star}(f)=+\infty$,
which concludes the proof.

##### Case $p\in]1,+\infty[$:

The reasoning proceeds as follows: we show that (i)
$\left\lVert{f}\right\rVert_{q}\leq 1$ implies
$\left\lVert{\cdot}\right\rVert_{p}^{\star}(f)=0$, (ii)
$1<\left\lVert{f}\right\rVert_{q}<+\infty$ gives
$\left\lVert{\cdot}\right\rVert_{p}^{\star}(f)=+\infty$, and (iii)
$\left\lVert{f}\right\rVert_{q}=+\infty$ results in
$\left\lVert{\cdot}\right\rVert_{p}^{\star}(f)=+\infty$. This allows us to
conclude that
$\left\lVert{\cdot}\right\rVert_{p}^{\star}=\iota_{\mathcal{B}_{1}^{q}}(\cdot)$.

  * •

When $\left\lVert{f}\right\rVert_{q}\leq 1$: By Hölder’s inequality, it holds
that

$\langle
f,g\rangle_{\mathscr{Y}}\leq\left\lVert{f}\right\rVert_{q}\left\lVert{g}\right\rVert_{p}\text{
for all }g\in\mathscr{Y}.$

Exploiting $\left\lVert{f}\right\rVert_{q}\leq 1$, we get that

$\langle f,g\rangle_{\mathscr{Y}}-\left\lVert{g}\right\rVert_{p}\leq 0\text{
for all }g\in\mathscr{Y}.$

The supremum being reached for $g=0$; we conclude that
$\left\lVert{\cdot}\right\rVert_{p}^{\star}(f)=0$.

  * •

When $1<\left\lVert{f}\right\rVert_{q}<+\infty$: According to Lemma A.9, there
exist $g\in\mathscr{Y}$ and $C>0$ such that

$\langle f,g\rangle_{\mathscr{Y}}-\left\lVert{g}\right\rVert_{p}\geq C.$

Denoting by $t>0$ a running parameter, one arrives at

$\langle f,tg\rangle_{\mathscr{Y}}-\left\lVert{tg}\right\rVert_{p}\geq
tC\xrightarrow{t\to\infty}+\infty.$

This shows that $\left\lVert{\cdot}\right\rVert_{p}^{\star}(f)=+\infty$.

  * •

When $\left\lVert{f}\right\rVert_{q}=+\infty$: We consider the sequence of
functions $(f_{n})_{n\in\mathbb{N}}$ defined as $f_{n}(\theta)=f(\theta)$ if
$\left\lvert f(\theta)\right\rvert\leq n$ and $f_{n}(\theta)=0$ otherwise,
where $(n,\theta)\in\mathbb{N}\times\Theta$. Each $f_{n}$ is bounded, thus
belongs to $L^{q}[\Theta,\mu]$, and the monotone convergence theorem applied
to the functions $\left\lvert f_{n}\right\rvert^{q}$ states that
$\left\lVert{f_{n}}\right\rVert_{q}\xrightarrow{n\to\infty}\left\lVert{f}\right\rVert_{q}=+\infty$.
Thus, there exists $N\in\mathbb{N}$ such that
$\left\lVert{f_{N}}\right\rVert_{q}>1$. We can then apply Lemma A.9 to get
$g\in\mathscr{Y}$ and $C>0$ such that

$\langle f_{N},g\rangle_{\mathscr{Y}}-\left\lVert{g}\right\rVert_{p}\geq C.$

According to Lemma A.9, $g(\theta)=0$ whenever $f_{N}(\theta)=0$, which
ensures that

$\langle f,g\rangle_{\mathscr{Y}}=\langle f_{N},g\rangle_{\mathscr{Y}}.$

Taking a running parameter $t>0$, this means that

$\langle f_{N},tg\rangle_{\mathscr{Y}}-\left\lVert{tg}\right\rVert_{p}=\langle
f,tg\rangle_{\mathscr{Y}}-\left\lVert{tg}\right\rVert_{p}\geq
tC\xrightarrow{t\to\infty}+\infty,$

which shows that $\left\lVert{\cdot}\right\rVert^{\star}(f)=+\infty$.

##### Case $p=+\infty$:

The reasoning goes as follows: we show that
$\left\lVert{f}\right\rVert_{1}\leq 1$ implies
$\left\lVert{\cdot}\right\rVert_{\infty}^{\star}(f)=0$, and that
$\left\lVert{f}\right\rVert_{1}>1$ gives
$\left\lVert{\cdot}\right\rVert_{\infty}^{\star}(f)=+\infty$, which allows one
to conclude that
$\left\lVert{\cdot}\right\rVert_{\infty}^{\star}=\iota_{\mathcal{B}^{1}_{1}}(\cdot)$.

  * •

When $\left\lVert{f}\right\rVert_{1}\leq 1$: By applying Hölder’s inequality
we get that $\langle
f,g\rangle_{\mathscr{Y}}\leq\left\lVert{f}\right\rVert_{1}\left\lVert{g}\right\rVert_{\infty}\text{
for all }g\in\mathscr{Y}$. Using the condition that
$\left\lVert{f}\right\rVert_{1}\leq 1$, this means that $\langle
f,g\rangle_{\mathscr{Y}}-\left\lVert{g}\right\rVert_{\infty}\leq 0\text{ for
all }g\in\mathscr{Y}$. Since the supremum is reached for $g=0$, we get that
$\left\lVert{\cdot}\right\rVert_{\infty}^{\star}(f)=0$.

  * •

When $\left\lVert{f}\right\rVert_{1}>1$: Let
$g\colon\theta\mapsto\operatorname{sign}{(f(\theta))}$. Since $g$ is bounded
by $1$, it belongs to $\mathscr{Y}$, and $\langle
f,g\rangle_{\mathscr{Y}}=\left\lVert{f}\right\rVert_{1}$. Running a free
parameter $t>0$, this means that $\langle
f,tg\rangle_{\mathscr{Y}}-t\left\lVert{g}\right\rVert_{\infty}=t\underbrace{(\left\lVert{f}\right\rVert_{1}-1)}_{>0}\xrightarrow{t\to\infty}+\infty$
which implies that
$\left\lVert{\cdot}\right\rVert_{\infty}^{\star}(f)=+\infty$.

∎

###  A.3 Proof of Proposition 3.2

######  Proposition (3.2).

Let $\kappa>0$, $p\in[1,+\infty]$, and $q$ the dual exponent of $p$ (i.e.,
$\frac{1}{p}+\frac{1}{q}=1$). Then for all $f\in\mathscr{Y}$,

$H_{\kappa}^{p}(f)=\left\\{\begin{matrix}\begin{array}[]{cl}\frac{1}{2}\|f\|_{\mathscr{Y}}^{2}&\textrm{if
}\left\lVert{f}\right\rVert_{q}\leq\kappa\\\\[4.26773pt]
\frac{1}{2}\left\lVert{\operatorname{Proj}_{\mathcal{B}_{\kappa}^{q}}(f)}\right\rVert_{\mathscr{Y}}^{2}+\kappa\left\lVert{f-\operatorname{Proj}_{\mathcal{B}_{\kappa}^{q}}(f)}\right\rVert_{p}&\textrm{otherwise}.\end{array}\end{matrix}\right..$

##### Proof

Let us introduce the notation
$R(g)=\frac{1}{2}\left\lVert{f-g}\right\rVert_{\mathscr{Y}}^{2}+\kappa\left\lVert{g}\right\rVert_{p}$
where $f\in\mathscr{Y}$, $g\in\mathscr{Y}$. Then

$\displaystyle H_{\kappa}^{p}(f)$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}\inf_{g\in\mathscr{Y}}R(g)\stackrel{{\scriptstyle(b)}}{{=}}R(\operatorname{prox}_{\kappa\left\lVert{\cdot}\right\rVert_{p}}(f))\stackrel{{\scriptstyle(c)}}{{=}}\frac{1}{2}\left\lVert{\operatorname{Proj}_{\mathcal{B}_{\kappa}^{q}}(f)}\right\rVert_{\mathscr{Y}}^{2}+\kappa\left\lVert{f-\operatorname{Proj}_{\mathcal{B}_{\kappa}^{q}}(f)}\right\rVert_{p},$
(15)

where (a) follows from the definition of the infimal convolution, (b) is
implied by that of the proximal operator using that
$\kappa\left\lVert{\cdot}\right\rVert_{p}\in\Gamma_{0}(\mathscr{Y})$. (c) is a
consequence of the Moreau decomposition (Lemma A.6) as

$\displaystyle\operatorname{prox}_{\kappa\left\lVert{\cdot}\right\rVert_{p}}(f)$
$\displaystyle=f-\operatorname{prox}_{\left(\kappa\left\lVert{\cdot}\right\rVert_{p}\right)^{\star}}(f)\stackrel{{\scriptstyle(d)}}{{=}}f-\operatorname{prox}_{\iota_{\mathcal{B}_{\kappa}^{q}}}(f)\stackrel{{\scriptstyle(e)}}{{=}}f-\operatorname{Proj}_{\mathcal{B}_{\kappa}^{q}}(f),$
(16)

where in (d) and (e) we used that

$\displaystyle\left(\kappa\left\lVert{\cdot}\right\rVert_{p}\right)^{\star}$
$\displaystyle\stackrel{{\scriptstyle(f)}}{{=}}\iota_{\mathcal{B}_{\kappa}^{q}}(\cdot)\text{
with $\frac{1}{p}+\frac{1}{q}=1$},$ (17)
$\displaystyle\operatorname{prox}_{\iota_{\mathcal{B}_{\kappa}^{q}}}$
$\displaystyle\stackrel{{\scriptstyle(g)}}{{=}}\operatorname{Proj}_{\mathcal{B}_{\kappa}^{q}}.$
(18)

(f) follows from the facts listed in the 3rd and the 2nd line of Table 5:

$\displaystyle\left(\kappa\left\lVert{\cdot}\right\rVert_{p}\right)^{\star}$
$\displaystyle=\kappa\left(\left\lVert{\cdot}\right\rVert_{p}\right)^{\star}(\cdot/\kappa)=\kappa\iota_{\mathcal{B}_{1}^{q}}(\cdot/\kappa)=\iota_{\mathcal{B}_{\kappa}^{q}}(\cdot).$

(g) is implied by
$\iota_{\mathcal{B}_{\kappa}^{q}}=\iota_{\mathcal{B}_{1}^{q}}(\cdot/\kappa)$,
the precomposition rule of proximal operators
($\operatorname{prox}_{f(\alpha\cdot)}=\frac{1}{\alpha}\operatorname{prox}_{\alpha^{2}f}(\alpha\cdot)$
holding for any $\alpha>0$; see (2.2) in (Parikh & Boyd, 2014)), and
$\operatorname{prox}_{\iota_{\mathcal{B}_{1}^{q}}}=\operatorname{Proj}_{\mathcal{B}_{1}^{q}}$:

$\displaystyle\operatorname{prox}_{\iota_{\mathcal{B}_{\kappa}^{q}}}$
$\displaystyle=\operatorname{prox}_{\iota_{\mathcal{B}_{1}^{q}}(\cdot/\kappa)}=\kappa\operatorname{prox}_{\frac{1}{\kappa^{2}}\iota_{\mathcal{B}_{1}^{q}}}(\cdot/\kappa)=\kappa\operatorname{prox}_{\iota_{\mathcal{B}_{1}^{q}}}(\cdot/\kappa)=\kappa\operatorname{Proj}_{\mathcal{B}_{1}^{q}}(\cdot/\kappa)=\operatorname{Proj}_{\mathcal{B}_{\kappa}^{q}}.$

Finally we note that when
$f=\operatorname{Proj}_{\mathcal{B}_{\kappa}^{q}}(f)$ ($\Leftrightarrow
f\in\mathcal{B}_{\kappa}^{q}\Leftrightarrow\left\lVert{f}\right\rVert_{q}\leq\kappa$),
$\eqref{eq:Hkp}$ simplifies to $\frac{1}{2}\|f\|_{\mathscr{Y}}^{2}$. ∎

###  A.4 Proof of Proposition 3.4

######  Proposition (3.4).

Let $p,\kappa\in[1,+\infty]\times]0,+\infty[$ and $\frac{1}{p}+\frac{1}{q}=1$.
Then the dual of (10) with loss $H_{\kappa}^{p}$ writes as

$\displaystyle\inf_{(\alpha_{i})_{i\in[n]}\in\mathscr{Y}^{n}}$
$\displaystyle\sum_{i\in[n]}\frac{1}{2}\left\lVert{\alpha_{i}}\right\rVert_{\mathscr{Y}}^{2}-\langle\alpha_{i},y_{i}\rangle_{\mathscr{Y}}+\frac{1}{2\lambda
n}\sum_{i,j\in[n]}^{n}k_{\mathscr{X}}(x_{i},x_{j})\left\langle\alpha_{i},T_{k_{\Theta}}\alpha_{j}\right\rangle_{\mathscr{Y}}\text{
s.t. }\left\lVert{\alpha_{i}}\right\rVert_{q}\leq\kappa\text{ for }\forall
i\in[n].$ (18)

##### Proof

Applying Lemma A.7 and (17) give the result. ∎

###  A.5 Proof of Proposition 3.5

######  Proposition (3.5).

Let $\kappa>0$. The projection on $\mathcal{B}_{\kappa}^{q}$ is tractable for
$q=2$ and $q=\infty$ and can be expressed for all
$(\alpha,\theta)\in\mathscr{Y}\times\Theta$ as

$\displaystyle\operatorname{Proj}_{\mathcal{B}_{\kappa}^{2}}(\alpha)$
$\displaystyle=\min{\left(1,\frac{\kappa}{\left\lVert{\alpha}\right\rVert_{2}}\right)}\alpha,\quad\text{if
}\alpha\neq 0$ (19)
$\displaystyle\left(\operatorname{Proj}_{\mathcal{B}_{\kappa}^{\infty}}(\alpha)\right)(\theta)$
$\displaystyle=\operatorname{sign}{(\alpha(\theta))}\min{(\kappa,|\alpha(\theta)|)}.$
(20)

##### Proof

The projection on the $2$-ball of radius $\kappa$ is similar to the finite-
dimensional case $(\mathscr{Y}=\mathbb{R}^{d})$ for which Equation 19 is well-
known.

We now turn to the case of $q=\infty$. Let $\alpha\in\mathscr{Y}$. By
definition,

$\operatorname{Proj}_{\mathcal{B}_{\kappa}^{\infty}}(\alpha)=\mathop{\mathrm{arg\,min}}_{y\in\mathscr{Y}}\frac{1}{2}\int_{\Theta}\left[\alpha(\theta)-y(\theta)\right]^{2}\mathrm{d}\mu(\theta)+\iota_{\mathcal{B}_{\kappa}^{\infty}}(y).$
(21)

Since $\alpha\in\mathscr{Y}$, the function $g$ defined as

$g(\theta)=\operatorname{sign}{(\alpha(\theta))}\min{(\kappa,|\alpha(\theta)|)},\quad\theta\in\Theta,$

is measurable, it is in $\mathscr{Y}$, and one can verify easily that it is
the solution of Equation 21; it corresponds to taking the pointwise projection
of $\alpha$ on the segment $[-\kappa,\kappa]$. ∎

###  A.6 Proof of Proposition 4.2

######  Proposition (4.2).

Let $\epsilon>0$ and $p\in[1,+\infty]$. Then
$\ell_{\epsilon}^{p}(f)=\frac{1}{2}\left\lVert{f-\operatorname{Proj}_{\mathcal{B}_{\epsilon}^{p}}(f)}\right\rVert^{2}_{\mathscr{Y}}$
for all $f\in\mathscr{Y}$.

##### Proof

Let
$R(g)=\frac{1}{2}\left\lVert{f-g}\right\rVert_{\mathscr{Y}}^{2}+\iota_{\mathcal{B}_{\epsilon}^{p}}(g)$
where $f\in\mathscr{Y}$ and $g\in\mathscr{Y}$. Then

$\displaystyle\ell_{\epsilon}^{p}(f)$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}\inf_{g\in\mathscr{Y}}R(g)\stackrel{{\scriptstyle(b)}}{{=}}R\left(\operatorname{prox}_{\iota_{\mathcal{B}_{\epsilon}^{p}}}(f)\right)\stackrel{{\scriptstyle(c)}}{{=}}R\left(\operatorname{Proj}_{\mathcal{B}_{\epsilon}^{p}}(f)\right)\stackrel{{\scriptstyle(d)}}{{=}}\frac{1}{2}\left\lVert{f-\operatorname{Proj}_{\mathcal{B}_{\epsilon}^{p}}(f)}\right\rVert^{2}_{\mathscr{Y}},$

where (a) follows from the definition of the infimal convolution, (b) is
implied by that of the proximal operator and by
$\iota_{\mathcal{B}_{\epsilon}^{p}}\in\Gamma_{0}(\mathscr{Y})$, (c) is the
consequence of
$\operatorname{prox}_{\iota_{\mathcal{B}_{\epsilon}^{p}}}(f)=\operatorname{Proj}_{\mathcal{B}_{\epsilon}^{p}}(f)$
implied by (18), in (d) the definition of $R$ was applied. ∎

###  A.7 Proof of Proposition 4.3

######  Proposition (4.3).

Let $(p,\epsilon)\in[1,+\infty]\times]0,+\infty[$, and
$\frac{1}{p}+\frac{1}{q}=1$. Then, the dual of (10) writes as

$\displaystyle\inf_{(\alpha_{i})_{i\in[n]}\in\mathscr{Y}^{n}}$
$\displaystyle\sum_{i\in[n]}\left[\frac{1}{2}\left\lVert{\alpha_{i}}\right\rVert_{\mathscr{Y}}^{2}-\langle\alpha_{i},y_{i}\rangle_{\mathscr{Y}}+\epsilon\left\lVert{\alpha_{i}}\right\rVert_{q}\right]+\frac{1}{2\lambda
n}\sum_{i,j\in[n]}k_{\mathscr{X}}(x_{i},x_{j})\left\langle\alpha_{i},T_{k_{\Theta}}\alpha_{j}\right\rangle_{\mathscr{Y}}.$
(21)

##### Proof

Applying Lemma A.7 with

$\displaystyle\left(\iota_{\mathcal{B}_{\epsilon}^{p}}(\cdot)\right)^{\star}=\epsilon\left\lVert{\cdot}\right\rVert_{q}$
(22)

gives the result. (22) is the consequence of (17) and the involutive property
of the Fenchel-Legendre conjugate. ∎

###  A.8 Proof of Proposition 4.4

######  Proposition (4.4).

Let $\epsilon>0$. The proximal operator of
$\epsilon\left\lVert{\cdot}\right\rVert_{q}$ is computable for $q=1$ and
$q=2$, and given for all $(\alpha,\theta)\in\mathscr{Y}\times\Theta$ by

$\displaystyle\left(\operatorname{prox}_{\epsilon\left\lVert{\cdot}\right\rVert_{1}}(\alpha)\right)(\theta)$
$\displaystyle=\operatorname{sign}{(\alpha(\theta))}\left|\left|\alpha(\theta)\right|-\epsilon\right|_{+},$
(23)
$\displaystyle\operatorname{prox}_{\epsilon\left\lVert{\cdot}\right\rVert_{2}}(\alpha)$
$\displaystyle=\alpha\left|1-\frac{\epsilon}{\left\lVert{\alpha}\right\rVert_{\mathscr{Y}}}\right|_{+}\quad\text{if
}\alpha\neq 0.$ (24)

##### Proof

By (16) we know that

$\operatorname{prox}_{\epsilon\left\lVert{\cdot}\right\rVert_{q}}(f)=f-\operatorname{Proj}_{\mathcal{B}_{\epsilon}^{p}}(f).$
(25)

The projection operator is known from Proposition 3.5 in the case of $p=2$ and
$p=\infty$, which allows to express the proximal operator of the $q$-norm for
$q=2$ and $q=1$ and by substituting Equation 19 and Equation 20 into Equation
25. ∎

##  Appendix B Additional Details

In this section we present additional algorithmic details as well as
complement the numerical experiments presented in the main document.

###  B.1 Algorithmic Details

Algorithm 2 fully describes how to learn models with the representation
relying on the eigendecomposition of the integral operator developed in
Section 3.2 and Section 4.2.

input : Gram matrix $\mathbf{K}_{\mathscr{X}}$, matrix of eigenvalues
$\mathbf{\Delta}$, data scalar product matrix $\mathbf{R}$, regularization
parameter $\lambda$, loss parameters $(\kappa,2)$ or $(\epsilon,2)$, gradient
step $\gamma$

init : $\mathbf{A}^{(0)},\mathbf{A}^{(-1)}=\bm{0}\in\mathbb{R}^{n\times r}$

for  _epoch $t$ from $0$ to $T-1$_ do

// gradient step

$\mathbf{V}=\mathbf{A}^{(t)}+\frac{t-2}{t+1}\left(\mathbf{A}^{(t)}-\mathbf{A}^{(t-1)}\right)$

$\mathbf{U}=\mathbf{V}-\gamma\left(\mathbf{V}+\frac{1}{\lambda
n}\mathbf{K}_{\mathscr{X}}\mathbf{V}\mathbf{\Delta}-\mathbf{R}\right)$

// proximal step

for  _row $i\in[n]$_ do

$\mathbf{a}_{i}^{(t+1)}=\min\left(\frac{\kappa}{\left\lVert{\mathbf{u}_{i}}\right\rVert_{2}},1\right)\mathbf{u}_{i}$
// if $H_{\kappa}^{2}$

$\mathbf{a}_{i}^{(t+1)}=\left|1-\frac{\gamma\epsilon}{\left\lVert{\mathbf{u}_{i}}\right\rVert_{2}}\right|_{+}\mathbf{u}_{i}$
// if $\ell_{\epsilon}^{2}$

return _$\mathbf{A}^{(T)}$_

Algorithm 2 APGD with eigenbasis representation

###  B.2 Synthetic Data

Below we detail the generation process of the synthetic dataset (Section
B.2.1), we expose in full detail the parameters used in the experiments
(Section B.2.2; see Fig. 1 and 2 in the main paper), and we provide additional
illustration for the interaction between the Huber loss’ $\kappa$ and the
regularization parameter $\lambda$ (Section B.2.3).

####  B.2.1 Generation Process

Given covariance parameters
$(\boldsymbol{\sigma}^{\text{in}},\boldsymbol{\sigma}^{\text{out}})\in\mathbb{R}^{r}\times\mathbb{R}^{r}$
for $c\in[r]$ we draw and fix Gaussian processes
$g_{c}^{\text{in}}\sim\mathcal{GP}(0,k_{\sigma^{\text{in}}_{c}})$ and
$g_{c}^{\text{out}}\sim\mathcal{GP}(0,k_{\sigma^{\text{out}}_{c}})$.

Figure 3: Examples from the toy dataset and corresponding type 2 outliers.

We then generate $n$ samples as
$\left(\sum_{c\in[r]}u_{ic}g_{c}^{\text{in}},\sum_{c\in[r]}u_{ic}g_{c}^{\text{out}}\right)_{i\in[n]}$,
where the coefficients $u_{ic}$ are drawn i.i.d. according to a uniform
distribution $\mathcal{U}([-0.5,0.5])$. In the experiments, we take $r=4$ and
set
$\boldsymbol{\sigma}^{\text{in}}=\boldsymbol{\sigma}^{\text{out}}=(0.05,0.1,0.5,0.7)$.
We show input and output functions drawn in this manner in the first and
second row of Fig. 3. In the bottom row we display outliers of Type 2 with
$\boldsymbol{\sigma}=(0.01,0.05,1,4)$ and intensity $\zeta=2$. For the
contaminated indices $i$ in $I$ we add the corresponding outlier to the
function $y_{i}$

####  B.2.2 Experimental Details

We provide here the full details of the parameters used for the experiments on
the toy dataset. For all experiments, we fix the parameter $\rho^{\text{in}}$
of the input Gaussian kernel
$k_{\mathscr{X}}:(x_{1},x_{2})\longmapsto\exp\left(-\rho\|x_{0}-x_{1}\|_{\mathscr{X}}^{2}\right)$
to $\rho^{\text{in}}=0.01$ and that of the output Gaussian kernel to
$\rho^{\text{out}}=100$. Indeed, since we are only given discrete observations
for the input functions as well, we use the available observations to
approximate the norms in the above kernels. For the experiments on robustness
which results are displayed in Fig. 2 of the main paper, we select via cross-
validation the regularization parameter $\lambda$ and the $\kappa$ parameters
of the Huber loss, considering values in a geometric grid of size $10$ ranging
from $10^{-6}$ to $10^{-3}$ for $\lambda$ and values in a geometric grid of
size $25$ ranging from $10^{-3}$ to $10^{-1}$ for $\kappa$.

####  B.2.3 Additional Illustrations

To highlight the interaction between the regularization parameter $\lambda$
and the parameter $\kappa$ of the Huber loss, we plot the NMSE values for
various values of $\lambda$ and $\kappa$ using the toy dataset corrupted with
the two main types of contamination used in the main paper, Type 2 and Type 3
outliers. The results are displayed in Fig. 4 and confirm that by making
$\kappa$ and $\lambda$ vary, when the data are corrupted, we can always find a
configuration that is significantly more robust than the square loss. In
accordance with one’s expectation, when dealing with local outliers (Fig.
4(b)), the loss $H_{\kappa}^{1}$ is much more efficient than the loss
$H_{\kappa}^{2}$. However, when dealing with global outliers (Fig. 4(a)), the
two losses perform equally well.

(a) Type 2 outliers with $\tau=0.2$ and $\zeta=2$

(b) Type 3 outliers with $\tau=0.3$ and $\xi=0.3$

Figure 4: NMSE as a function of $\lambda$ and Huber losses’ $\kappa$ with two
types of outliers.

###  B.3 DTI Data

In this section we provide details regarding the experiments on the DTI
dataset. For this dataset, we use a Gaussian kernel as input kernel and a
Laplace kernel as output kernel, for the first we fix its parameter to
$\rho^{\text{in}}=1.25$, and for the second, defined as
$k_{\Theta}:(\theta_{1},\theta_{2})\longmapsto\exp(-\rho^{\text{out}}\|x_{0}-x_{1}\|_{\mathscr{X}})$,
we fix its parameter to $\rho^{\text{out}}=10$. We consider two values of
$\lambda$, the first one ($\lambda=10^{-5}$) is chosen too small for the
square loss to highlight the additional sparsity-inducing regularization
possibilities offered by the $\epsilon$-insensitive loss through the parameter
$\epsilon$, while the second one ($\lambda=10^{-3}$) corresponds to a near-
optimal value for the square loss. We do cross-validate the parameters of the
losses. For the loss $\ell_{\epsilon}^{2}$ we consider values of $\epsilon$ in
a geometric grid of size 50 ranging from $10^{-3}$ to $10^{-1}$, while for the
loss $\ell_{\epsilon}^{\infty}$, we search in a geometric grid of the same
size, however ranging from $10^{-3}$ to $10^{-0.5}$. For the Huber losses
$H_{\kappa}^{1}$ and $H_{\kappa}^{2}$, we search for $\kappa$ using a
geometric grid of size 50 ranging this time from $10^{-4}$ to $10^{-1}$.

###  B.4 Speech Data

This section is dedicated to additional details about the experiments carried
out on the speech benchmark.

Input kernel: As highlighted in the main paper, we encode the input sounds
through 13 mel-frequency cepstral coefficients (MFCC). To deal with this
particular input data type we used the following kernel. Let
$((x_{ij}^{(v)})_{v\in[13]})_{i\in[n],j\in[m]}$ be the transformed input data
where $v$ serves as an index for the MFCC number. The number of locations $m$
is the same for all $i\in[n]$ since we extend the signals to match the longest
one to be able to train the models. We then center and reduce each MFCC using
all samples and sampling locations to compute the mean and standard deviation;
let $((\tilde{x}_{ij}^{(v)})_{v\in[13]})_{i\in[n],j\in[m]}$ be the resulting
standardized input data. Denoting by $\tilde{x}_{i}$ the element
$((\tilde{x}_{ij}^{(v)})_{v\in[13]})_{j\in[m]}$, we then use the following
kernel:

$k_{\mathscr{X}}:(\tilde{x}_{0},\tilde{x}_{1})\longmapsto\frac{1}{m}\sum_{j\in[m]}\exp\left(-\rho^{\text{in}}\sum_{v\in[13]}\left(\tilde{x}_{0j}^{(v)}-\tilde{x}_{1j}^{(v)}\right)^{2}\right).$

Experimental details: For all the experiments (with or without corruption), we
select the parameter of the input kernel $\rho^{\text{in}}$, the
regularization parameter and the parameters of the losses using cross-
validation. We fix the parameter of the Laplace output kernel to
$\rho^{\text{out}}=10$. However, to reduce the computational burden, we
perform the selection of the parameter $\rho^{\text{in}}$ only for the square
loss, and then take the corresponding values for the other losses. For this
parameter values in a geometric grid of size $15$ ranging from $10^{-2}$ to
$10^{-0.5}$ are considered. For $\lambda$, the search space is a geometric
grid of size $10$ ranging from $10^{-10}$ to $10^{-6}$. Finally, for the
$\epsilon$-insensitive loss, values of $\epsilon$ in a geometric grid of size
$80$ ranging from $10^{-5}$ to $10^{-1}$ are considered, while for the Huber
losses we search for $\kappa$ in a geometric grid of size $100$ ranging from
$10^{-7}$ to $1$.

##  Appendix C Illustration of Loss Functions

Figure 5: Sensitivity of $H_{\kappa}^{p}$ to outliers for various $p\in[1,2]$

In this section we illustrate the differences between our proposed convoluted
losses in several ways. In Section C.1 we study empirically how the choice of
$p$ affects the sensitivity of the Huber loss $H_{\kappa}^{p}$ to different
kind of outliers. In Section C.2, we plot some of our proposed losses when
they are defined either on $\mathbb{R}$ or $\mathbb{R}^{2}$.

###  C.1 Discussion on the Choice of $p$ for $H_{\kappa}^{p}$

As it is highlighted in the main paper, solving Problem ((4)) for
$p\not\in\\{1,2\\}$ is unpractical since it involves the computation of $n$
projections on a $q$-ball at each APGD iteration. Performing such projection
is feasible (it is a convex optimization problem) but it has to be done in an
iterative way. In our case, to run APGD with such inner iterations turns out
to be too time consuming. However we still can approximately calculate the
losses $H_{\kappa}^{p}$ using Proposition 3.2 for any $p$ (computing the
involved projection iteratively). We thus propose to leverage this possibility
to study empirically the sensitivity of the Huber losses $H_{\kappa}^{p}$ to
global and local outliers, for different values of $p$.

The impact of the outliers on the solution of a regularized empirical risk
minimization problem is partly determined by the contribution of the outliers
to the data-fitting term relatively to the contribution of the normal
observations. In order to investigate this aspect, we study and define next a
quantity which we call Robustness Ratio.

Let $(e_{i})_{i\in[n]}\in(L^{2}[\Theta,\mu])^{n}$ be a set of functional
residuals and let $(\tilde{e}_{i})_{i\in[n]}$ be the same functional residuals
but contaminated with outliers. In practice, we have to choose a probability
distribution to draw the functions $(e_{i})_{i\in[n]}$ from, and an outlier
distribution to corrupt those. For the functions $(e_{i})_{i\in[n]}$ we use
our synthetic data generation process (see Section B.2.1), and for the
outliers, we consider the same type 2 and type 3 outliers as in the
experiments in Section 5.1 from the main paper. We then define the Robustness
Ratio as

Robustness Ratio $\displaystyle:=\inf_{\kappa\geq
0}\frac{1}{n}\sum_{i\in[n]}\frac{H_{\kappa}^{p}(\tilde{e}_{i})}{H_{\kappa}^{p}(e_{i})}.$

The best value of this quantity is $1$; it means that the loss is not affected
at all by the outliers, but it is indeed not possible to reach such value. In
practice, we restrain our study to $p\in[1,2]$. For each $p$ we reduce the
search for $\kappa$ to different empirical quantiles of the $q$-norms of the
uncorrupted functions $(e_{i})_{i\in[n]}$, where $q$ is the dual exponent of
$p$. It makes sense to do so since $\kappa$ corresponds to a $q$-norm
threshold which separates observations considered to be outliers from those
deemed normal (see Proposition 3.2). We consider the
$\\{0.5,0.6,0.7,0.8,0.9,0.95,0.99\\}$-th such empirical quantiles. For each
$p$, we compute the robustness ratio for $\kappa$ equal to each of those
quantiles, and then for each level of corruption, we select the value which
minimizes the ratio. This indeed corresponds to an ideal setting, since in
practice, we never have access to the uncorrupted data and we can never
optimize $\kappa$ in this way. Thus the robustness ratio reflects more of a
general robustness property of the loss in an optimal setting.

In accordance with one’s expectation, when the data is contaminated with
global outliers (left panel of Fig. 5), it is better to choose $p=2$ whereas
when the contamination is local (right panel of Fig. 5), $p=1$ is almost the
best choice; even though it seems that choosing $p$ slightly bigger than $1$
could be a tad better. Even though, we highlight that this analysis based on
the Robustness Ratio has its limits; indeed we do not take into account the
interplay between the data-fitting term and the regularization term which
takes place during optimization. This certainly explains why we found the
losses $H_{\kappa}^{1}$ and $H_{\kappa}^{2}$ to perform equally well in
practice whereas based only on the Robustness Ratio analysis (left panel of
Fig. 5) we would have said otherwise. The findings in presence of local
outliers (right panel of Fig . 5) are nevertheless coherent with what we
observed in practice for the losses $H_{\kappa}^{1}$ and $H_{\kappa}^{2}$ in
our experiments.

###  C.2 Loss Examples in 1d and 2d

In this section, we plot several of the proposed losses when they are defined
on $\mathbb{R}$ and $\mathbb{R}^{2}$. In Fig. 6, we compare the Huber (Fig.
6(b)) and the $\epsilon$-insensitive (Fig. 6(a)) losses with the square loss
when they are defined on $\mathbb{R}$.

Then in Fig. 7 we highlight the influence of $p$ on the shape of the
$\epsilon$-insensitive loss $\ell_{\epsilon}^{p}$ defined on $\mathbb{R}^{2}$.
We set $\epsilon=1$ and consider values of $p\in\\{1.01,1.5,2,3,5,+\infty\\}$.
We display $\ell_{\epsilon}^{1.01}$ in Fig. 7(a), $\ell_{\epsilon}^{1.5}$ in
Fig. 7(b), $\ell_{\epsilon}^{2}$ in Fig. 7(c), $\ell_{\epsilon}^{3}$ in Fig.
7(d), $\ell_{\epsilon}^{5}$ in Fig. 7(e) and $\ell_{\epsilon}^{\infty}$ in
Fig. 7(f).

Finally, in Fig. 8 we underline the influence that the parameter $p$ has on
our proposed Huber losses when it is defined on $\mathbb{R}^{2}$; we take
$\kappa=0.8$ and we display $H_{\kappa}^{2}$ in Fig. 8(a), $H_{\kappa}^{1.5}$
in Fig. 8(b), $H_{\kappa}^{1.25}$ in Fig. 8(c) and $H_{\kappa}^{1}$ in Fig.
8(d).

(a) $\epsilon$-insensitive ($\epsilon=1$) and square loss

(b) Huber loss ($\kappa=0.8$) and square loss

Figure 6: Illustrations of the different losses defined on $\mathbb{R}$.

(a) $\ell_{\epsilon}^{1.01}$ ($\epsilon=1$)

(b) $\ell_{\epsilon}^{1.5}$ ($\epsilon=1$)

(c) $\ell_{\epsilon}^{2}$ ($\epsilon=1$)

(d) $\ell_{\epsilon}^{3}$ ($\epsilon=1$)

(e) $\ell_{\epsilon}^{5}$ ($\epsilon=1$)

(f) $\ell_{\epsilon}^{\infty}$ ($\epsilon=1$)

Figure 7: Examples of the proposed $\epsilon$-insensitive losses defined on
$\mathbb{R}^{2}$ for different values of $p$.

(a) $H_{\kappa}^{2}$ ($\kappa=0.8$)

(b) $H_{\kappa}^{1.5}$ ($\kappa=0.8$)

(c) $H_{\kappa}^{1.25}$ ($\kappa=0.8$)

(d) $H_{\kappa}^{1}$ ($\kappa=0.8$)

Figure 8: Examples of the proposed Huber losses defined on $\mathbb{R}^{2}$
for different values of $p$.