# Robust risk management via multi-marginal optimal transport 111The work of
BP was partially supported by National Sciences and Engineering Research
Council of Canada Discovery Grant number 04658-2018. He is also pleased to
acknowledge the generous hospitality of the Institut de Mathématiques d’Orsay,
Université Paris-Sud during his stay in November of 2019 as a missionaire
scientifique invité, when this work was partially completed. L.N. was
supported by a public grant as part of the Investissement d’avenir project,
reference ANR-11-LABX-0056-LMH, LabEx LMH, from the CNRS PEPS JCJC (2022) and
H-CODE.

Hamza Ennaji Normandie Université, ENSICAEN, UNICAEN, CNRS, GREYC, France.
<EMAIL_ADDRESS>Quentin Mérigot Université Paris-Saclay, CNRS,
Laboratoire de mathématiques d’Orsay, 91405, Orsay, France.
<EMAIL_ADDRESS>Luca Nenna Université Paris-Saclay,
CNRS, Laboratoire de mathématiques d’Orsay, 91405, Orsay, France.
<EMAIL_ADDRESS>Brendan Pass Department of Mathematical
and Statistical Sciences, 632 CAB, University of Alberta, Edmonton, Alberta,
Canada, T6G 2G1<EMAIL_ADDRESS>
###### Abstract

We study the problem of maximizing a spectral risk measure of a given output
function which depends on several underlying variables, whose individual
distributions are known but whose joint distribution is not. We establish and
exploit an equivalence between this problem and a multi-marginal optimal
transport problem. We use this reformulation to establish explicit, closed
form solutions when the underlying variables are one dimensional, for a large
class of output functions. For higher dimensional underlying variables, we
identify conditions on the output function and marginal distributions under
which solutions concentrate on graphs over the first variable and are unique,
and, for general output functions, we find upper bounds on the dimension of
the support of the solution. We also establish a stability result on the
maximal value and maximizing joint distributions when the output function,
marginal distributions and spectral function are perturbed; in addition, when
the variables one dimensional, we show that the optimal value exhibits
Lipschitz dependence on the marginal distributions for a certain class of
output functions. Finally, we show that the equivalence to a multi-marginal
optimal transport problem extends to maximal correlation measures of multi-
dimensional risks; in this setting, we again establish conditions under which
the solution concentrates on a graph over the first marginal.

###### Contents

  1. 1 Introduction
  2. 2 Preliminaries and problem formulation
    1. 2.1 Maximal spectral risk measures over couplings of given marginals
  3. 3 Equivalence between maximizing spectral risk measures and multi-marginal optimal transport
  4. 4 Solutions for one-dimensional assets and compatible payouts
    1. 4.1 Sensitivity analysis and maximal river flow
  5. 5 Higher dimensional assets
  6. 6 Stability
  7. 7 Multidimensional measures of risk

##  1 Introduction

In a variety of problems in operations research, a variable of interest
$b=b(x_{1},x_{2},\dotsc,x_{d})$ depends on several underlying random
variables, whose individual distributions are known (or can be estimated) but
whose joint distribution is not. A natural example arises in finance, when one
considers the payout of a derivative depending on several underlying assets.
An estimate of the distribution of the asset values themselves can often be
inferred from the prices of vanilla call and put options for a wide range of
strike prices; since these options are widely traded, their prices are readily
available. However, estimating the joint distribution would require prices of
a wide range of derivatives with payouts depending on all the variables, which
are typically much scarcer. Other examples arise when the variables represent
parameters in a physical system whose individual distributions can be
estimated empirically or through modeling (or a combination of both) but whose
dependence structure cannot; an example of this flavour, originating in [18],
in which the output $b$ represents the simulated height of a river at risk of
flooding, and the underlying variables include various design parameters and
climate dependent factors can be found in Section 4.1 below. A third example
comes from election projections in political science, where $b$ might
represent the outcome of an election and the $x_{i}$ vote shares in different
regions, whose individual distributions can be modelled from polling data.

Methods used in risk management to evaluate the resulting aggregate level of
risk naturally depend on the distribution of the output variable $b$, and
therefore, in turn, on the joint distribution of the $x_{i}$. A natural
problem is therefore to determine bounds, or worst-case scenarios, on these
risk measures; that is, to maximize a given risk measure over all possible
joint distributions of the $x_{i}$ with known marginal distributions. Problems
such as this have received extensive attention within the financial risk
management community; see, for example, [32]. In the simplest of these
problems, the underlying variables $x_{i}$ are typically real valued and the
output variable $b$ is often assumed to have a particular structure (in many
cases, it is a weighted sum of the $x_{i}$, reflecting the value of a
portfolio built out of underlying assets with values $x_{i}$, or a function of
this weighted sum). In these cases, explicit solutions for the maximizing
couplings of the $x_{i}$ can sometimes be obtained (see, for example, [29, 30,
36]).

For more general output functions $b$, and possibly multi-dimensional
underlying variables $x_{i}$, much less is known about the dependence
structure of the maximizing joint distributions. A recent paper of Ghossoub-
Hall-Saunders [15] studies this more general setting systematically, allowing
the underlying variables to take values in very general spaces and the output
$b$ to take a very general form, and focuses on _spectral risk measures_ (see
Definition 10 below). They observe that the resulting problem is a
generalization of optimal transport (see Section 2 for a brief overview of the
optimal transport problem); in fact, in the simplest case, when the spectral
function is identically equal to one and the number $d$ of underlying
variables is $2$, the problem is exactly a classical optimal transport problem
with surplus function $b$. For more general spectral risk measures, the
problem has the _same constraint set_ as the optimal transport problem, but a
more general, _non-linear objective functional_ ; they adapt the duality
theory of optimal transport to this setting and establish results on the
stability with respect to perturbations of $b$ and the marginal distributions.
Although the analysis in [15] focused on the $d=2$ case, they note that their
results can be extended to the $d\geq 3$ setting, in which case maximizing
spectral risks becomes a generalization of the multi-marginal optimal
transport problem (see Section 2 for a brief overview of this problem). As in
the two marginal case, one maximizes a concave objective function over the set
of couplings of the given marginals – when the spectral function is constant,
one obtains multi-marginal optimal transport, a linear maximization.

The main purpose of this paper is to establish and exploit a simpler, but
equivalent, reformulation of this problem. Specifically, we show that for
_any_ spectral risk measure, the maximization can in fact be formulated as a
traditional multi-marginal optimal transport problem with $d+1$ marginals: the
given marginals distributions of the $x_{i}$ as well as another distribution
arising from the particular form of the spectral function. Although this
formulation slightly increases the underlying dimension, it makes the
maximization problem linear and much more tractable – indeed, the results and
techniques, both theoretical and computational, in the substantial literature
on multi-marginal optimal transport become directly applicable.

Though the structure of solutions to multi-marginal optimal transport problems
is in general a notoriously delicate issue (see [26] for an overview) there
are many important cases when the problem is well understood. As a consequence
of our result, when the underlying variables $x_{i}$ are all one dimensional,
we derive an explicit characterization of solutions for a substantial class of
output functions $b$, through a careful refinement of the existing theory of
multi-marginal optimal transport on one dimensional ambient spaces. For $b$
falling outside this class, explicit solutions are likely generally
unattainable; however, our formulation of the problem can potentially
facilitate the use of a very broad range of computational methods for optimal
transport to approximate solutions numerically (see Section 8.1 in [2] and the
references therein).

For underlying variables in higher dimensional spaces, explicit solutions are
generally not possible. However, we show that known techniques from multi-
marginal optimal transport can be adapted to identify conditions under which
the solution concentrates on a graph over the first variable; solutions having
this structure, commonly referred to as _Monge solutions_ in the optimal
transport literature, are analogous to the well known comonotone couplings
which often maximize spectral risk measures in the one dimensional case, since
they essentially assert that knowing the value of the first variable
completely determines the values of the others. On the other hand, for many
output functions $b$, conditions ensuring this structure of a solution do not
hold; in this case, we demonstrate that available techniques can still be used
bound the dimension of the support of the optimizer. In addition to being of
theoretical interest, this fact is potentially important in future work on the
computation of solutions.222In general, the joint distributions of interest
are probability measures on an $n\cdot d$ dimensional space, where $d$ is the
number of variables, and $n$ the dimension of each variable (assuming, for
simplicity of exposition, that all variables are of the same dimension). When
maximizers concentrate on graphs over $x_{1}$, they are essentially much
simpler, $n$-dimensional objects. More generally, it is possible to prove that
optimizers concentrate on $m$ dimensional subsets of the $n\cdot d$
dimensional ambient space, where $m$ satisfies $n\leq m\leq n\cdot d$ and is
determined by the output function $b$; this has consequences for the covering
number of the support, which will potentially play an important role in future
computation.

We go on to show that our formulation, combined with standard stability
results for optimal transport problems, implies that the maximal value of the
spectral risk measure, and the maximizing joint distributions, are stable with
respect to perturbations of the marginals, output function and spectral
function. When the underlying variables are one dimensional and the output
function satisfies additional hypotheses, we further leverage the connection
to optimal transport to show that the maximal value’s dependence on the
marginals is in fact Lipschitz.

We take particular note of the special but important case when the spectral
risk measure in question is the _expected shortfall_ ; in this case, the
problem further reduces to a _multi-marginal partial transport problem_ on the
$d$ original distributions. In this case, in higher dimensions, results in [6,
12] and [20] immediately apply, implying uniqueness and graphical structure of
solutions for a certain class of output functions $b$. Furthermore, the fact,
exposed by our work here, that optimal partial transport can be reformulated
as a multi-marginal optimal transport problem does not seem to have been
observed before, and may be of independent interest in the optimal transport
community.

We also extend these ideas to the setting where the output function $b$ is
multi-variate valued, using the maximal correlation risk measures developed in
[11]. We are again able to identify conditions under which the solution
concentrates on a graph.

The paper is organized as follows: in Section 2, we first introduce the
traditional (two marginal) optimal transport problem and its multi-marginal
generalization, and establish certain preliminary results which we will need.
We then introduce spectral risk measures, and the problem of maximizing them
over joint distributions with fixed marginals, which is the main object of
interest in this paper (see (13) below). In Section 3, we establish the
equivalence between this maximization problem and multi-marginal optimal
transport. In Section 4, we establish explicit solutions for the maximizing
joint distributions when the underlying marginal variables are one dimensional
and in the fifth section, we study how these result generalize to the setting
where the underlying variables are higher dimensional. Section 6 is focused on
the stability of the problem with respect to variations in the marginals,
output function and spectral risk measure, while the final section is reserved
for the extension of some of our results to maximal correlation risk measures
on multi-variate output functions.

##  2 Preliminaries and problem formulation

#### Assumption and notations

Let $\mu_{i}\in\mathcal{P}(X_{i})$ be probability measures on Polish spaces
$X_{i}$, for $i=0,1,\dotsc,d$. Given a mapping $T:X_{i}\rightarrow X_{j}$, we
say that $T$ pushes $\mu_{i}$ forward to $\mu_{j}$, and write

$T_{\\#}\mu_{i}=\mu_{j},\;\text{if}\;\mu_{i}(T^{-1}(A))=\mu_{j}(A)\;\text{for
all Borel}\;A\subseteq X_{j}.$

We will let
$\Gamma(\mu_{0},\mu_{1},\dotsc,\mu_{d})\subset\mathcal{P}(X_{0}\times
X_{1}\times\dotsm\times X_{d})$ denote the set of probability measures on
$X_{0}\times X_{1}\times\dotsm\times X_{d}$ whose marginals are the $\mu_{i}$;
that is, $\mu_{i}=\Big{(}(x_{0},x_{1},x_{2},\dotsc,x_{d})\mapsto
x_{i}\Big{)}_{\\#}\pi$ for each $i$ and each
$\pi\in\Gamma(\mu_{0},\mu_{1},\dotsc,\mu_{d})$.

######  Definition 1 ($(1,\dotsc,d)$-marginal).

For $\pi\in\mathcal{P}(X_{0}\times X_{1}\times\dotsm\times X_{d})$, we call
the _$(1,\dotsc,d)$ -marginal_ of $\pi$, the projection of $\pi$ onto
$X_{1}\times\dotsm\times X_{d}$; that is,
$\Big{(}(x_{0},x_{1},x_{2},\dotsc,x_{d})\mapsto(x_{1},x_{2},\dotsc,x_{d})\Big{)}_{\\#}\pi$.
Note that if $\pi\in\Gamma(\mu_{0},\mu_{1},\dotsc,\mu_{d})$ then clearly its
$(1,\dotsc,d)$-marginal is in $\Gamma(\mu_{1},\dotsc,\mu_{d})$.

As usual we define the cumulative distribution function (c.d.f) of a
probability measure $\mu\in\mathcal{P}(\mathbb{R})$ as
$F_{\mu}(x)=\mu((-\infty,x])$. The generalized or pseudo-inverse of $F_{\mu}$
is defined as

$F_{\mu}^{-1}(m)=\inf\\{x\in\mathbb{R}\;|\;F_{\mu}(x)\geq m\\}.$

We recall that $F_{\mu}^{-1}(m)$ is the value of the $m$th quantile.

#### Optimal Transport problem

Given two probability distributions $\mu_{0}$ and $\mu_{1}$ on Polish spaces
$X_{0}$ and $X_{1}$, respectively, and a surplus function $s:X_{0}\times
X_{1}\to\mathbb{R}$, the optimal transport problem (in its Monge formulation)
consists in maximising

$\int_{X_{0}}s(x,T(x))\mathrm{d}\mu_{0}(x)$ (1)

under the constraint that $T_{\\#}\mu_{0}=\mu_{1}$ (namely $\mu_{1}$ is the
image measure of $\mu_{0}$ through the map $T$). This is a delicate problem
since the mass conservation constraint $T_{\\#}\mu_{0}=\mu_{1}$ is highly
nonlinear (and the feasible set may even be empty, for instance if $\mu_{0}$
is a Dirac mass and $\mu_{1}$ is not). For this reason, it is common to study
a relaxed formulation of (1) which allows mass splitting; that is,

$\max_{\pi\in\Gamma(\mu,\nu)}\int_{X_{0}\times X_{1}}s(x,y)\mathrm{d}\pi(x,y)$
(2)

where, as above, $\Gamma(\mu_{0},\mu_{1})$ consists of all probability
measures on $X_{0}\times X_{1}$ having $\mu_{0}$ and $\mu_{1}$ as marginals.
Note that this is a (infinite dimensional) linear programming problem and that
there exists solutions under very mild assumptions on $s,\mu_{0}$ and
$\mu_{1}$. A solution $\pi$ of (2) is called an optimal transport plan; in
particular, if an optimal plan of (2) has a deterministic form
$\pi=(\mathrm{id},T)_{\sharp}\mu_{0}$ (which means that no splitting of mass
occurs and $\pi$ is concentrated on the graph of $T$) then $T$ is an optimal
transport map, i.e., a solution to (1). It is therefore sometimes called a
_Monge solution_ of (2). The linear problem (2) also has a convenient dual
formulation

$\min_{\begin{subarray}{c}u,v\\\ u(x)+v(y)\geq
s(x,y)\end{subarray}}\int_{X_{0}}u(x)\mathrm{d}\mu_{0}(x)+\int_{X_{1}}v(y)\mathrm{d}\mu_{1}(y)$
(3)

where $u(x)$ and $v(y)$ are the so-called Kantorovich potentials. Optimal
Transport theory for two marginals has developed very rapidly in the 25 last
years; there are well known conditions on $s$, $\mu_{0}$ and $\mu_{1}$ which
guarantee that there is a unique deterministic optimal plan and we refer to
the textbooks by Santambrogio [33] and Villani [34, 35], for a detailed
exposition.

#### Multi-marginal optimal transport

In this paper we consider a generalization of (1) to the case in which more
than two marginals are involved: given probability measures
$\mu_{0},\mu_{1},\dotsc,\mu_{d}$ on spaces $X_{0},X_{1},\dotsc,X_{d}$, and a
surplus function $s(x_{0},x_{1},\dotsc,x_{d})$, the multi-marginal optimal
transport problem is to maximize

$\int_{X_{0}\times X_{1}\times\dotsm\times
X_{d}}s(x_{0},x_{1},\dotsc,x_{d})\mathrm{d}\pi(x_{0},x_{1},\dotsc,x_{d})$ (4)

over $\pi\in\Gamma(\mu_{0},\mu_{1},\dotsc,\mu_{d})$. This problem has
attracted increasing attention in recent years since it arises naturally in
many different settings such as economics [10], quantum chemistry [5, 9],
fluid mechanics [4], etc. Our results in Section 3 will bring forth a new
application to maximizing spectral risk measures.

It is well known that under mild conditions (for example, continuity of $s$
and compact support of the $\mu_{i}$ is more than sufficient) (4) admits a
solution; see, for example, [26]. Although the structure of solutions to (4)
can generally depend delicately on $s$, there is a growing theory, and an
important class for which solutions can be derived essentially explicitly,
which we describe below. In addition, there are now powerful numerical tools
to compute solutions for problems falling outside this class (for a more
detailed discussion, we refer to Section 8.1 in [2] and the references
therein).

The most tractable of these problems occurs when the underlying asset spaces
are one dimensional, and the mixed second derivatives of $s$ interact in a
certain way.

######  Definition 2 (Compatibility).

Suppose each $X_{i}\subset\mathbb{R}$ is a bounded real interval. Assume that
$s\in\mathcal{C}^{2}(X_{0}\times X_{1}\times\dotsm\times X_{d})$. We say that
$s$ is strictly compatible if for each three distinct indices
$i,j,k\in\\{0,1,\dotsc,d\\}$ and each $(x_{0},\ldots,x_{d})\in
X_{0}\times\dotsm\times X_{d}$ we have

$\frac{\partial^{2}s}{\partial x_{i}\partial
x_{j}}\Bigg{[}\frac{\partial^{2}s}{\partial x_{k}\partial
x_{j}}\Bigg{]}\frac{\partial^{2}s}{\partial x_{k}\partial
x_{i}}(x_{0},\dotsc,x_{d})>0.$ (5)

We will say that $s$ is weakly compatible (or simply compatible) if for each
$i\neq j$, we have either $\frac{\partial^{2}s}{\partial x_{i}\partial
x_{j}}\geq 0$ throughout $X_{0}\times X_{1}\times\dotsm\times X_{d}$ or
$\frac{\partial^{2}s}{\partial x_{i}\partial x_{j}}\leq 0$ throughout
$X_{0}\times X_{1}\times\dotsm\times X_{d}$, and

$\frac{\partial^{2}s}{\partial x_{i}\partial
x_{j}}\Bigg{[}\frac{\partial^{2}s}{\partial x_{k}\partial
x_{j}}\Bigg{]}\frac{\partial^{2}s}{\partial x_{k}\partial
x_{i}}(x_{0},\dotsc,x_{d})\geq 0,$ (6)

for each distinct $i,j,k\in\\{0,1,\dotsc,d\\}$ and each
$(x_{0},\ldots,x_{d})\in X_{0}\times\dotsm\times X_{d}$.

The fact that the mixed partials $\frac{\partial^{2}s}{\partial x_{i}\partial
x_{j}}$ do not change signs under the compatibility condition allows us to
partition the set $\\{0,1,2,\dotsc,d\\}=S_{+}\cup S_{-}$ of indices into
disjoint subsets $S_{+}$ and $S_{-}$ such that $0\in S_{+}$ and for each
$i\neq j$, $\frac{\partial^{2}s}{\partial x_{i}\partial x_{j}}\geq 0$
throughout $X_{0}\times X_{1}\times\dotsm\times X_{d}$ if either both $i$ and
$j$ are in $S_{+}$ or if both are in $S_{-}$, and
$\frac{\partial^{2}s}{\partial x_{i}\partial x_{j}}\leq 0$ throughout
$X_{0}\times X_{1}\times\dotsm\times X_{d}$ otherwise. The inequalities are
strict if strict compatibility holds . A well known special case of
compatibility is captured by the following definition.

######  Definition 3 (Supermodularity).

Suppose each $X_{i}\subset\mathbb{R}$ is a connected real interval. Assume
that $s\in\mathcal{C}^{2}(X_{0}\times X_{1}\times\dotsm\times X_{d})$. We say
that $s$ is supermodular if $\frac{\partial^{2}s}{\partial x_{i}\partial
x_{j}}\geq 0$ for all $i\neq j$; we say $s$ is strictly supermodular if
$\frac{\partial^{2}s}{\partial x_{i}\partial x_{j}}>0$ for all $i\neq j$

######  Remark 4.

Notice that if $s$ is supermodular, then it is clearly compatible. On the
other hand, if $s$ is compatible, define
$\tilde{s}(\tilde{x}_{0},\tilde{x}_{1},\dotsc,x_{d})=s(x_{0},x_{1},\dotsc,x_{d})$,
where $\tilde{x}_{0}=x_{0}$ and for each $i=1,2,\dotsc,d$ we set
$\tilde{x}_{i}=x_{i}$ if $i\in S_{+}$ and $\tilde{x}_{i}=-x_{i}$ if $i\in
S_{-}$. Then it is straightforward to show that $\tilde{s}$ is supermodular.
Therefore, compatibility is simply supermodularity up to a change of
variables.

Both supermodularity and compatibility may be extended to non $C^{2}$
functions, by considering partial differences; for simplicity of exposition,
we stick with the version relying on second derivatives here.

We then introduce the following special coupling in
$\Gamma(\mu_{0},\mu_{1},\dotsc,\mu_{d})$.

######  Definition 5.

For a compatible $s$, we define the _$s$ \- comonotone coupling_ by:

$\pi=(G_{0},G_{1},\dotsc,G_{d})_{\\#}\mathrm{Leb}_{[0,1]},$ (7)

where $G_{0}=F^{-1}_{\mu_{0}}$ and for each $i=1,2,\ldots,d$

$G_{i}(m)=\begin{cases}F^{-1}_{\mu_{i}}(m)\;&\hbox{if}\;i\in S_{+},\\\
F^{-1}_{\mu_{i}}(1-m)\;&\hbox{if}\;i\in S_{-}.\end{cases}$ (8) 333There may be
multiple ways to define the $G_{i}$ (if, for example, for some $i$,
$\frac{\partial^{2}s}{\partial x_{i}\partial x_{j}}=0$ everywhere for all
$j$); in this case, the results below apply to all resulting such $s$ \-
comonotone couplings. Note, however, that the $G_{i}$ are uniquely defined
when compatibility is strict.

######  Remark 6.

If $s$ is supermodular, we can take $S_{-}$ to be empty and
$S_{+}=\\{0,1,\dotsc,n\\}$; the $s$-comonotone coupling is then exactly the
well known classical comonotone coupling.

The following result can be found in [7] for strictly supermodular functions
as well as in [22] for compatible functions (where its proof appears together
with the formulation of the compatibility conditions and the observation that
it is equivalent to strict supermodularity after changing coordinates).

######  Theorem 7.

Suppose that $s$ is strictly compatible. Then the $s$-comonotone coupling (7)
is the unique optimizer in (4).

######  Remark 8.

A classical but important case occurs when $d=2$ and
$s(x_{0},x_{1})=x_{0}x_{1}$. It is well known that there is only one measure
in $\Gamma(\mu_{0},\mu_{1})$ with monotone increasing support. This will be
used in the proof of Lemma 16 below.

For $d\geq 3$, the $s$-comonotone $\pi$ is characterized by the fact that its
twofold marginals,
$\pi_{ij}=((x_{0},x_{1},\ldots,x_{d})\mapsto(x_{i},x_{j}))_{\\#}\pi$ are the
monotone increasing (respectively decreasing) couplings if
$\frac{\partial^{2}s}{\partial x_{i}\partial x_{j}}>0$ (respectively
$\frac{\partial^{2}s}{\partial x_{i}\partial x_{j}}<0$); compatibility ensures
existence of a $\pi$ with this property.

Though the following result does not seem to be available in the literature,
it is a straightforward consequence of Theorem 7 and the well known stability
of optimal transport with respect to perturbations of the cost function. It
asserts that the optimality of the measure constructed in Theorem 7 still
holds if the strong compatibility assumption is relaxed to weak compatibility,
although the uniqueness assertion may not.

######  Proposition 9.

Suppose that $s$ is weakly compatible. Then the $s$-comonotone coupling (7) is
optimal in (4).

###### Proof.

For $\epsilon>0$, define

$\begin{split}s_{\epsilon}(x_{0},x_{1},\dotsc,x_{d})&=s(x_{0},x_{1},\dotsc,x_{d})+\epsilon\sum_{i,j\in
S_{+}}(x_{i}-x_{j})^{2}+\epsilon\sum_{i,j\in S_{-}}(x_{i}-x_{j})^{2}\\\
&-\epsilon\sum_{i\in S_{+},j\in S_{-}}(x_{i}-x_{j})^{2}\end{split}$

notice now that $s_{\epsilon}$ is strictly compatible and so Theorem 7 implies
that $\pi=(G_{0},G_{1},\dotsc,G_{d})_{\\#}\mathrm{Leb}_{[0,1]}$ is optimal in
(4) for each $\epsilon$. Therefore, for any other
$\tilde{\pi}\in\Gamma(\mu_{0},\mu_{1},\dotsc,\mu_{d})$ we have

$\int s_{\epsilon}\mathrm{d}\tilde{\pi}\leq\int s_{\epsilon}\mathrm{d}\pi$

Since the $s_{\epsilon}$ converge uniformly to $s$ as $\epsilon\rightarrow 0$,
we can take the limit in the inequality above and obtain

$\int s\mathrm{d}\tilde{\pi}\leq\int s\mathrm{d}\pi.$

As $\tilde{\pi}\in\Gamma(\mu_{0},\mu_{1},\dotsc,\mu_{d})$ was arbitrary, this
yields the desired result. ∎

We also note that, as for the two marginal problem, (4) admits a very useful
dual formulation:

$\begin{split}\inf\bigg{\\{}\sum_{i=0}^{d}\int_{X_{i}}u_{i}(x_{i})\mathrm{d}\mu_{i}(x_{i})\;&|\;u_{i}\in\mathcal{C}_{b}(X_{i}),\forall
i=0,\dotsc,d,\\\ &\sum_{i=0}^{d}u_{i}(x_{i})\geq
s(x_{0},\dotsc,x_{d})\bigg{\\}}.\end{split}$ (9)

As with the primal problem (4), it is well known that under mild conditions a
solution $(u_{0},u_{1},\dotsc,u_{d})$ to (9) exists; see, for example, [14, 8,
23]. In particular, given optimal solutions $\pi$ and $(u_{0},\dotsc,u_{d})$
to (4) and (9), respectively, the following optimality condition holds

$\sum_{i=0}^{d}u_{i}(x_{i})=s(x_{0},\dotsc,x_{d}),\;\pi-\text{a.e.}$ (10)

moreover, by Proposition 9 for a compatible $s$ (10) rewrites

$\sum_{i=0}^{d}u_{i}(G_{i}(m))=s(G_{0}(m),\dotsc,G_{d}(m)),\;\mathrm{Leb}-\text{a.e..}$

###  2.1 Maximal spectral risk measures over couplings of given marginals

Let the real line describe a certain level of risk (e.g., the level of
radiation in the nuclear power plant) and $\mu$ be a probability measure on
$\mathbb{R}$ which can be interpreted as the distribution of risk. We will
consider the following form of quantifier of the risk associated with $\mu$.

######  Definition 10 (spectral risk measure).

A functional $\mathcal{R}:\mathcal{P}(\mathbb{R})\rightarrow\mathbb{R}$ is
called a _spectral risk measure_ if it takes the form
$\mathcal{R}=\mathcal{R}_{\alpha}$, where

$\mathcal{R}_{\alpha}(\mu):=\int_{0}^{1}F_{\mu}^{-1}(m)\alpha(m)\mathrm{d}m.$
(11)

######  Remark 11.

Note that by taking $\alpha$ to be defined and real valued on $[0,1]$, we are
tacitly assuming that it is bounded, since
$\alpha(0)\leq\alpha(m)\leq\alpha(1)$ by monotonicity; the same assumption was
made in [15]. This is mostly for technical convenience; we expect that most of
our results can be extended to the more general case
$\alpha:[0,1)\rightarrow\mathbb{R}_{+}$, where possibly $\lim_{m\rightarrow
1}\alpha(m)=\infty$, under suitable hypothesis (for instance, decay conditions
on $\alpha_{\\#}\mathrm{Leb}_{[0,1]}$).

For a given $\alpha$, we will often refer to the spectral risk measure
$\mathcal{R}_{\alpha}$ as the _$\alpha$ -risk_, and $\alpha$ as the the
_spectral function._

Spectral risk measures are pervasive in risk management and insurance. Indeed,
it is well known that any coherent risk measure which is additive for
comontonic random variables must be spectral (see for instance [13]).

We begin with the following variational characterization of spectral risk
measures. This result is well known to experts (see, for example, [31] and
[11], in which this fact is the basis for a multi-variate extension). We
include a brief proof here for the convenience of the reader

######  Lemma 12.

For any spectral risk measure and any probability measure $\mu$ with
$\int_{\mathbb{R}}x\mathrm{d}\mu(x)>-\infty$,

$\mathcal{R}_{\alpha}(\mu)=\max_{\pi\in\Gamma(\alpha_{\sharp}\mathrm{Leb}_{[0,1]},\mu)}\int_{\mathbb{R}\times\mathbb{R}}xy\mathrm{d}\pi(x,y),$
(12)

where we have denoted $\Gamma(\alpha_{\sharp}\mathrm{Leb}_{[0,1]},\mu)$ the
space of probability measures on $\mathbb{R}^{2}$ with marginals
$\alpha_{\sharp}\mathrm{Leb}_{[0,1]}$ and $\mu$.

###### Proof.

Let us re-write the right hand side of (12) as the following optimal transport
problem between $\mu_{0}=\alpha_{\sharp}\mathrm{Leb}_{[0,1]}$ and
$\mu_{1}=\mu$

$\max_{\gamma\in\Gamma(\mu_{0},\mu_{1})}\int_{[0,1]\times\mathbb{R}}s(x,y)\mathrm{d}\pi(x,y),$

where $s(x,y)=xy$. Then, since the surplus is strictly supermodular there
exists a unique optimal $\pi$ of the form
$\pi=(G_{0},G_{1})_{\sharp}\mathrm{Leb}_{[0,1]}$, where
$G_{0}(m)=F^{-1}_{\mu_{0}}(m)=F^{-1}_{\alpha_{\sharp}\mathrm{Leb}_{[0,1]}}(m)=\alpha(m)$
and $G_{1}(m)=F^{-1}_{\mu}(m)$, by Theorem 7. The desired result follows. ∎

######  Corollary 13.

$\mathcal{R}_{\alpha}$ is concave on $\mathcal{P}(\mathbb{R})$.

###### Proof.

Let $\mu_{0},\mu_{1}\in\mathcal{P}(\mathbb{R})$ and let $\pi_{0},\pi_{1}$ be
optimal for the problem (12). Then, with $\mu_{t}=(1-t)\mu_{0}+t\mu_{1}$ and
$\pi_{t}=(1-t)\gamma_{0}+\gamma_{1}$ one has
$\pi_{t}\in\Gamma(\alpha_{\sharp}\mathrm{Leb}_{[0,1]},\mu_{t})$ so that

$\begin{split}\mathcal{R}_{\alpha}(\mu_{t})&\geq\int
xy\mathrm{d}\pi_{t}(x,y)\\\ &=(1-t)\int xy\mathrm{d}\pi_{0}(x,y)+t\int
xy\mathrm{d}\pi_{1}(x,y)\\\
&=(1-t)\mathcal{R}_{\alpha}(\mu_{0})+t\mathcal{R}_{\alpha}(\mu_{1}).\end{split}$

∎

Return now to the case in which there are $d$ parameters which enter into the
estimation of risk through the output function $b$, and we wish to evaluate
the worst case scenario for the $\alpha$-risk
$\mathcal{R}_{\alpha}(b_{\\#}\gamma)$ of $b_{\\#}\gamma$ among couplings
$\gamma\in\Gamma(\mu_{1},\mu_{2},\dotsc,\mu_{d})$ of the marginal
distributions $\mu_{i}$. That is, we want to maximize:

$\max_{\gamma\in\Gamma(\mu_{1},\mu_{2},\dotsc,\mu_{d})}\mathcal{R}_{\alpha}(b_{\\#}\gamma).$
(13)

######  Remark 14 (Maximizing Expected Shortfall is optimal partial
transport).

A special case of particular importance in risk management applications occurs
when $\alpha=\alpha_{m_{0}}:=\frac{1}{m_{0}}\mathbf{1}_{[1-m_{0},1]}$, in
which case (11) is also know as the _Expected Shortfall_ , or _Conditional
Value at Risk_.

In this setting, the maximization problem (13) is actually equivalent to a
well known variant of the optimal transport problem; indeed, it can be
reformulated into

$\max_{\gamma\in\Gamma_{m_{0}}(\mu_{1},\dotsc,\mu_{d})}\frac{1}{m_{0}}\int
b(x_{1},\dotsc,x_{d})\mathrm{d}\gamma(x_{1},\dotsc,x_{d})$

where $\Gamma_{m_{0}}(\mu_{1},\dotsc,\mu_{d})$ denotes the set of non-negative
measures $\gamma$ on $X_{1}\times\dotsm\times X_{d}$ with total mass $m_{0}$,
such that its $i-$marginal $\gamma_{i}$ is dominated by $\mu_{i}$ for each
$i$; that is
$\int_{X_{i}}\phi\mathrm{d}\gamma_{i}\leq\int_{X_{i}}\phi\mathrm{d}\mu_{i}$
for all non-negative test functions $\phi\in\mathcal{C}^{0}(X_{i})$.

This is known as the _optimal partial transport problem_ when $d=2$ [6, 12],
and the _multi-marginal optimal partial transport problem_ [20] when $d>2$.

As we will show below, the problem is in fact equivalent to an ordinary multi-
marginal optimal transport problem with an additional marginal; see Theorem 18
and Remark 19 below.

######  Proposition 15.

Let $b:X_{1}\times X_{2}\times\dotsm\times X_{d}\rightarrow\mathbb{R}$ be an
upper semi-continuous function bounded from below. Then the map
$\gamma\mapsto\mathcal{R}_{\alpha}(b_{\\#}\gamma)$ is concave.

###### Proof.

Using (12) one gets

$\displaystyle\mathcal{R}_{\alpha}(b_{\\#}\gamma)$
$\displaystyle=\max_{\sigma\in\Gamma(\alpha_{\sharp}\mathrm{Leb}_{[0,1]},b_{\\#}\gamma)}\int
x_{0}y\mathrm{d}\sigma(x_{0},y)$
$\displaystyle=\max_{\pi\in\Gamma(\alpha_{\sharp}\mathrm{Leb}_{[0,1]},\gamma)}\int
x_{0}b(z)\mathrm{d}\pi(x_{0},z),$

where $z=(x_{1},\dotsc,x_{d})$. Then concavity of
$\gamma\mapsto\mathcal{R}_{\alpha}(b_{\\#}\gamma)$ on
$\mathcal{P}(\mathbb{R}^{d})$ follows as in Corollary 13. ∎

##  3 Equivalence between maximizing spectral risk measures and multi-marginal
optimal transport

Our first main contribution is to show that the spectral risk maximization
problem (13) is equivalent to the multi-marginal optimal transport problem (4)
with $X_{0}=[\alpha(0),\alpha(1)]\subseteq\mathbb{R}$,
$\mu_{0}=\alpha_{\sharp}\mathrm{Leb}_{[0,1]}$, the other $X_{i}$ and $\mu_{i}$
representing the domains and distributions of the underlying variables,
respectively, and

$s(x_{0},x_{1},\dotsc,x_{d})=x_{0}b(x_{1},\dotsc,x_{d}).$ (14)

######  Lemma 16.

Suppose that $\pi\in\Gamma(\mu_{0},\mu_{1},\dotsc,\mu_{d})$ and let $\gamma$
be the $(1,...,d)$ -marginal of $\pi$. Then, for $s$ given by (14),

$\int_{\times_{i=0}^{d}X_{i}}s(x_{0},x_{1},\dotsc,x_{d})\mathrm{d}\pi(x_{0},x_{1},\dotsc,x_{d})\leq\mathcal{R}_{\alpha}(b_{\\#}\gamma)$

Furthermore, we have equality if and only if the support of

$\tau_{\pi}=\Big{(}(x_{0},x_{1},x_{2},\dotsc,x_{d})\mapsto(x_{0},b(x_{1},x_{2},\dotsc,x_{d}))\Big{)}_{\\#}\pi\in\mathcal{P}(\mathbb{R}^{2})$

is monotone increasing.

###### Proof.

We first note that

$\begin{split}\int_{\times_{i=0}^{d}X_{i}}s(x_{0},x_{1},\dotsc,x_{d})\mathrm{d}\pi(x_{0},\dotsc,x_{d})&=\int_{\times_{i=0}^{d}X_{i}}x_{0}b(x_{1},\dotsc,x_{d})\mathrm{d}\pi(x_{0},\dotsc,x_{d})\\\
&=\int_{\mathbb{R}^{2}}zx_{0}\mathrm{d}\tau_{\pi}(x_{0},z)\\\
&\leq\mathcal{R}_{\alpha}(b_{\\#}(\gamma)),\end{split}$

by Lemma 12. Furthermore, the uniqueness in Theorem 7 implies that the
inequality is in fact an equality if and only if the coupling $\tau_{\pi}$ has
monotone increasing support.

∎

######  Lemma 17.

Given any measure $\gamma\in\Gamma(\mu_{1},\dotsc,\mu_{d})$, there exists a
$\pi\in\Gamma(\mu_{0},\mu_{1},\dotsc,\mu_{d})$ whose $(1,\dotsc,d)$-marginal
is $\gamma$, such that
$\Big{(}(x_{0},x_{1},\dotsc,x_{d})\mapsto(x_{0},b(x_{1},x_{2},\dotsc,x_{d})\Big{)}_{\\#}\pi$
has monotone increasing support.

###### Proof.

Let $\tau=(F_{\mu_{0}}^{-1},F_{b_{\\#}\gamma}^{-1})_{\\#}\mathrm{Leb}_{[0,1]}$
be the comonotonic coupling of $\mu_{0}$ and $b_{\\#}\gamma$. We disintegrate
$\tau(x_{0},y)=\tau^{y}(x_{0})\otimes(b_{\\#}\gamma)(y)$ with respect to its
second marginal, so that for any measurable function $g(x_{0},y)$ we have

$\int_{\mathbb{R}^{2}}g(x_{0},y)\mathrm{d}\tau(x_{0},y)=\int_{\mathbb{R}}\int_{\mathbb{R}}g(x_{0},y)\mathrm{d}\tau^{y}(x_{0})\mathrm{d}(b_{\\#}\gamma)(y),$

and
$\gamma(x_{1},...,x_{d})=\gamma^{y}(x_{1},...,x_{d})\otimes(b_{\\#}\gamma)(y)$
with respect to $b$, so that for any measurable $h(x_{1},...,x_{d})$ we have

$\int_{X_{1}\times\dotsm\times
X_{d}}h(x_{1},\dotsc,x_{d})\mathrm{d}\gamma(x_{1},\dotsc,x_{d})=\int_{\mathbb{R}}\int_{b^{-1}(y)}h(x_{1},\dotsc,x_{d})\mathrm{d}\gamma^{y}(x_{1},\dotsc,x_{d})\mathrm{d}(b_{\\#}\gamma).$

Any measure
$\pi(x_{0},x_{1}...,x_{d})=\pi^{y}(x_{0},x_{1},...,x_{d})\otimes(b_{\\#}\gamma)(y)$
on $X_{0}\times(X_{1}\times\dotsm\times X_{d})$, whose conditional
probabilities
$\pi^{y}\in\Gamma(\tau^{y},\gamma^{y})\subseteq\mathcal{P}(X_{0}\times
X_{1}\times...\times X_{d})$ after disintegrating with respect to the mapping
$(x_{0},\dotsc,x_{d})\mapsto(x_{0},b(x_{1},...,x_{d}))$ are couplings of
$\tau^{y}$ and $\gamma^{y}$ for $b_{\\#}\gamma$ almost every $y$ then
satisfies the requirement. ∎

We are now ready to establish the equivalence between the maximal
$\alpha-$risk problem (13) and the multi-marginal optimal transport problem
(4).

######  Theorem 18.

A probability measure $\pi$ in $\Gamma(\mu_{0},\mu_{1},\dotsc,\mu_{d})$ is
optimal for (4) with cost function (14) if and only if its
$(1,\dotsc,d)$-marginal is optimal in (13), and
$\tau_{\pi}=\Big{(}(x_{0},x_{1},x_{2},\dotsc,x_{d})\mapsto(x_{0},b(x_{1},x_{2},\dotsc,x_{d}))\Big{)}_{\\#}\pi$
has monotone increasing support.

###### Proof.

First, suppose that $\pi$ is optimal in (4) and let $\gamma$ be its
$(1,...,d)$-marginal. Then for any other
$\tilde{\gamma}\in\Gamma(,\mu_{1},...,\mu_{d})$, construct
$\tilde{\pi}\in\Gamma(\mu_{0},\mu_{1},...,\mu_{d})$ as in Lemma 17. We then
have, by Lemma 16 and optimality of $\pi$

$\displaystyle\mathcal{R}_{\alpha}(b_{\\#}\tilde{\gamma})$ $\displaystyle=$
$\displaystyle\int_{X_{0}\times X_{1}\times\dotsm\times
X_{d}}s(x_{0},...,x_{d})\mathrm{d}\tilde{\pi}(x_{0},\dotsc,x_{d})$
$\displaystyle\leq$ $\displaystyle\int_{X_{0}\times X_{1}\times\dotsm\times
X_{d}}s(x_{0},\dotsc,x_{d})\mathrm{d}\pi(x_{0},\dotsc,x_{d})$
$\displaystyle\leq$ $\displaystyle\mathcal{R}_{\alpha}(b_{\\#}\gamma).$

This establishes that $\gamma$ is optimal in (13). Furthermore, we must have
equality throughout if we choose $\tilde{\gamma}=\gamma$. In particular, this
implies that $\int_{X_{0}\times X_{1}\times\dotsm\times
X_{d}}s(x_{0},\dotsc,x_{d})\mathrm{d}\pi(x_{0},\dotsc,x_{d})=\mathcal{R}_{\alpha}(b_{\\#}\gamma)$,
which, by the second assertion in Lemma 16, implies that $\tau_{\pi}$ has
monotone increasing support.

Conversely, assume that $\pi\in\Gamma(\mu_{0},\mu_{1},\dotsc,\mu_{d})$, such
that its $(1,\dotsc,d$)-marginal $\gamma$ is optimal in (13) and $\pi$ couples
$\mu_{0}$ and $b_{\\#}\gamma$ monotonically. Let now $\tilde{\pi}$ be an
optimizer in (4) and $\tilde{\gamma}$ its $(1,\dotsc,d)$-marginal. We then
have, using Lemma 16,

$\displaystyle\mathcal{R}_{\alpha}(b_{\\#}\gamma)$ $\displaystyle=$
$\displaystyle\int_{X_{0}\times X_{1}\times\dotsm\times
X_{d}}s(x_{0},x_{1},\dotsc,x_{d})\mathrm{d}\pi(x_{0},x_{1},\dotsc,x_{d})$
$\displaystyle\leq$ $\displaystyle\int_{X_{0}\times X_{1}\times\dotsm\times
X_{d}}s(x_{0},x_{1},\dotsc,x_{d})\mathrm{d}\tilde{\pi}(x_{0},x_{1},\dotsc,x_{d})$
$\displaystyle\leq$ $\displaystyle\mathcal{R}_{\alpha}(b_{\\#}\tilde{\gamma})$

Optimality of $\gamma$ in (13) then implies that we must actually have
equality; in particular $\pi$ is optimal in (4) as desired. ∎

######  Remark 19.

Returning to the Expected Shortfall Case,
$\alpha=\alpha_{m_{0}}=\frac{1}{m_{0}}\mathbf{1}_{[1-m_{0},1]}$, in view of
Remark 14, Theorem 18 shows that the optimal partial transport problem for the
surplus $b(x_{1},\dotsc,x_{d})$ can be transformed into a multi-marginal
transport problem for the surplus
$s(x_{0},\dotsc,x_{d})=x_{0}b(x_{1},\dotsc,x_{d})$ and additional marginal
$\mu_{0}=\alpha_{\sharp}\mathrm{Leb}_{[0,1]}$. To the best of our knowledge,
this equivalence between these two well studied mathematical problems has not
been observed before in the optimal transport literature.

In addition, this perspective, together with results in [6, 12, 20] allows us
to immediately identify conditions under which the active (that is, the part
that couples to $\alpha>0$) part of the optimal $\gamma$ in (13) is uniquely
determined and concentrates on a graph over $x_{1}$. Furthermore, algorithms
to compute the solution are readily available [3, 17].

These questions will be revisited for more general spectral risk measures
later on.

We close this section with the following immediate consequence of Theorem 18
and a standard existence result in optimal transport theory.

######  Corollary 20.

Suppose that $b$ is bounded above and upper-semicontinuous on $X_{1}\times
X_{2}\times\dotsm\times X_{d}$. Then there exists a maximizing $\gamma$ in
(13).

###### Proof.

Note that the boundedness of $\alpha$ implies that $s$ given by (14) is
bounded above and upper-semicontinuous on $X_{0}\times X_{1}\times
X_{2}\times\dotsm\times X_{d}$. The existence of an optimal $\pi$ in (23) is
then well known; see for example, Theorem 1.7 in [33]. The $(1,\dotsc,d)-$
marginal of $\pi$ then maximizes (13) by Theorem 18. ∎

##  4 Solutions for one-dimensional assets and compatible payouts

We now turn our attention to the structure of maximizers in (13) when the
underlying variables are one dimensional.

Assume that each $x_{i}\in\mathbb{R}$ and each $\mu_{i}$ is supported on an
interval, $X_{i}=[\underline{x}_{i},\overline{x}_{i}]$, for all
$i=1,\dotsc,d$. We note that in our setting, the first marginal, $\mu_{0}$ is
supported always on the interval
$[\underline{x}_{0},\overline{x}_{0}]:=[\alpha(0),\alpha(1)]$ with
$\alpha(0)\geq 0$.

######  Lemma 21.

Suppose that $b$ is compatible and monotone increasing in each $x_{i}\in
S_{+}$ and monotone decreasing for each $x_{i}\in S_{-}$. Then the
$s$-comonotone coupling (7) is optimal for (4) with surplus function given by
(14). The maximal value is

$\int_{0}^{1}\alpha(m)b(G_{1}(m),G_{2}(m),\dotsc,G_{d}(m))\mathrm{d}m.$ (15)

Furthermore, if in addition the monotonicity of $b$ with respect to each
argument is strict, the solution is unique on the support of $\alpha$; that
is, any other solution $\overline{\pi}$ coincides with
$\pi=(G_{\alpha},G_{1},\dotsc,G_{d})_{\\#}\mathrm{Leb}_{[0,1]}$ on
$(0,\alpha(1)]\times[\underline{x}_{1},\overline{x}_{1}]\times\dotsm\times[\underline{x}_{d},\overline{x}_{d}]$.

Clearly full uniqueness of the optimal $\pi$ cannot hold if $\alpha=0$ on a
set of positive measure, or, equivalently, $\mu_{0}(\\{0\\})>0$; in this case,
we can rearrange the part of the $(1,...,d)-$ marginal $\gamma$ of $\pi$ that
couples to $\alpha=0$ in any way without affecting the value of $\int
s\mathrm{d}\pi$. The preceding Lemma identifies conditions under which this is
the only source of non-uniqueness, so that the optimal $\pi$ is uniquely
determined on
$(0,\alpha(1)]\times[\underline{x}_{1},\overline{x}_{1}]\times\dotsm\times[\underline{x}_{d},\overline{x}_{d}]$.

Note also that we are able to obtain uniqueness off the set where $\alpha=0$
despite the fact that the surplus function may not be strongly compatible.
This is done by exploiting an idea developed in [27] (see in particular
Proposition 4.1 there); namely that one only needs strongly compatible
interactions between certain key pairs of variables (rather than all of them).

###### Proof.

The conditions on $b$ make the function $s$ compatible, and so Proposition 9
implies the optimality of the $s$-comonotone $\pi$. Note that since $\alpha$
is in fact the quantile function for the first marginal
$\mu_{0}=\alpha_{\sharp}\mathrm{Leb}_{[0,1]}$, (15) is exactly $\int
s\mathrm{d}\pi$.

Now, to prove the last assertion, we will take advantage of the fact that
compatibility is equivalent to supermodularity up to a change of variables
(recall Remark 4). After such a change of variables, our assumptions become
that $b$ is supermodular, and strictly increasing in each argument, and each
$G_{i}=F^{-1}_{\mu_{i}}$ is the quantile map of $\mu_{i}$.

It will suffice to prove the following claim: for any optimizer
$\overline{\pi}$, if $(x_{0},x_{1},\dotsc,x_{d})$ and
$(\tilde{x}_{0},\tilde{x}_{1},\dotsc,\tilde{x}_{d})$ are both in the support
of $\overline{\pi}$ and $x_{0}>\tilde{x}_{0}\geq 0$, then
$x_{i}\geq\tilde{x}_{i}$ for each $i$. To see that this is sufficient, note
that applying the claim for $\tilde{x}_{0}=0$ implies that the lowest $a$
portion of each mass, $\\{x_{i}<G_{i}(a)\\}$ must pair with $x_{0}=0$, where
$a=\mu_{0}(\\{0\\})$. Applying the claim with $\tilde{x}_{0}>0$ then implies
that the support of $\overline{\pi}$ on
$(0,\alpha(1)]\times[\underline{x}_{1},\overline{x}_{1}]\times\dotsm\times[\underline{x}_{d},\overline{x}_{d}]$
must be monotone increasing, which immediately implies the desired result.

To see the claim, we apply the $s$-monotonicity property found, in, for
example [25, 7]. We define $x^{+}_{i}=\max\\{x_{i},\tilde{x}_{i}\\}$ and
$x^{-}_{i}=\min\\{x_{i},\tilde{x}_{i}\\}$. We then have

$s(x^{+}_{0},x^{+}_{1},\dotsc,x^{+}_{d})+s(x^{-}_{0},x^{-}_{1},\dotsc,x^{-}_{d})\leq
s(x_{0},x_{1},\dotsc,x_{d})+s(\tilde{x}_{0},\tilde{x}_{1},\dotsc,\tilde{x}_{d}).$
(16)

Note that $x_{0}^{-}=\tilde{x}_{0}$. If $\tilde{x}_{0}=0$, then
$s(x^{-}_{0},x^{-}_{1},\dotsc,x^{-}_{d})=s(\tilde{x}_{0},\tilde{x}_{1},\dotsc,\tilde{x}_{d})=0$,
so that (16) becomes (after dividing by $x_{0}$ )

$b(x^{+}_{1},\dotsc,x^{+}_{d})\leq b(x_{1},\dotsc,x_{d})$

As each $x^{+}_{i}\geq x_{i}$ and $b$ is strictly monotone in each coordinate,
this can only happen if each $x_{i}^{+}=x_{i}$, which is equivalent to the
claim.

On the other hand, if $\tilde{x}_{0}>0$, we follow the approach in the proof
of Proposition 4.1 in [27]. A straightforward calculation yields that

$\displaystyle
s(x^{+}_{0},x^{+}_{1},\dotsc,x^{+}_{d})+s(x^{-}_{0},x^{-}_{1},\dotsc,x^{-}_{d})-s(x_{0},x_{1},\dotsc,x_{d})-s(\tilde{x}_{0},\tilde{x}_{1},\dotsc,\tilde{x}_{d})$
$\displaystyle=\int_{0}^{1}\int_{0}^{1}\sum_{i\neq
j=0}^{d}\dfrac{\partial^{2}s(x(\theta,\phi))}{\partial x_{i}\partial
x_{j}}(x_{i}^{+}-x_{i})(x_{j}-x_{j}^{-})\mathrm{d}\theta\mathrm{d}\phi$ (17)

where

$x(\theta,\phi)=\phi[\theta
x^{+}+(1-\theta)x]+(1-\phi)[\theta\tilde{x}+(1-\theta)x^{-}].$

Now, by supermodularity of $s$, each
$\dfrac{\partial^{2}s(x(\theta,\phi))}{\partial x_{i}\partial x_{j}}\geq 0$,
and as each $x_{i}^{+}-x_{i}\geq 0$ and $x_{j}-x_{j}^{-}\geq 0$, we have that
(17) is greater than or equal to $0$; by (16), it must then be equal to $0$,
which then implies that each individual term
$\dfrac{\partial^{2}s(x(\theta,\phi))}{\partial x_{i}\partial
x_{j}}(x_{i}^{+}-x_{i})(x_{j}-x_{j}^{-})$ in the sum must be $0$. Now, strict
monotonicity of $b$ with respect to each $x_{i}$ implies that when we take
$j=0$, we have $\dfrac{\partial^{2}s(x(\theta,\phi))}{\partial x_{i}\partial
x_{0}}=\dfrac{\partial b(x(\theta,\phi))}{\partial x_{i}}>0$, which, as
$x_{0}>\tilde{x}_{0}=x_{0}^{-}$ by assumption, implies that each
$x_{i}^{+}=x_{i}$. This completes the proof of the claim, and, consequently,
the Lemma.

∎

Thanks to Theorem 18, the preceding result then easily yields the following
characterization of solutions to (13).

######  Theorem 22.

Assume that $b$ is compatible, monotone increasing in each $x_{i}\in S_{+}$
and monotone decreasing in each $x_{i}\in S_{-}$. Then
$(G_{1},\dotsc,G_{d})_{\\#}\mathrm{Leb}_{[0,1]}$ maximizes (13) and the
maximal value is given by (15). If in addition the monotonicity is strict,
then letting $a=\mu_{0}\\{0\\}$ any other maximizer $\gamma$ must couple the
regions defined by $G_{i}([0,a])$, and must coincide with the $b$-comonotone
coupling $(G_{1},\dotsc,G_{d})_{\\#}\mathrm{Leb}_{[0,1]}$ on
$G_{1}([a,1])\times G_{2}([a,1])\times\dotsm\times G_{d}([a,1])$.

Note that in the special case when $b$ is supermodular and increasing with
respect to each argument, this asserts optimality of the comonotone coupling
Note that each region $G_{i}([0,a])$ is of mass $a$; it is an interval of the
form $[\underline{x_{i}},G_{i}(a)]$ if $i\in S_{+}$ (corresponding to the
lowest fraction of the mass) and of the form $[G_{i}(a),\overline{x_{i}}]$ if
$i\in S_{-}$ (corresponding to the highest fraction of the mass).

######  Example 23.

Consider $b(x_{1},\dotsc,x_{d})=x_{1}+\dotsb+x_{d}$. This corresponds to the
value of a portfolio composed of underlying assets $x_{1},\dotsc,x_{d}$. It is
clearly weakly supermodular and therefore compatible, and strictly monotone
increasing in each coordinate, so the $(1,...,d)$-marginal of the $\pi$
defined by (7) is optimal by Theorem 22.

Explicitly, the optimal value in (13) is:

$\mathcal{R}_{\alpha}(b_{\\#}\gamma)=\int_{0}^{1}\alpha(m)(F_{\mu_{1}}^{-1}(m)+\dotsb+F_{\mu_{d}}^{-1}(m))\mathrm{d}m=\sum_{i=1}^{d}\mathcal{R}_{\alpha}(\mu_{i}),$

which corresponds to coupling “good” and “bad” events (i.e., with $x_{i}$
small or large, respectively) together.

Note that this special case could alternatively be deduced by noting that
spectral risk measures are sub-additive for all couplings, and additive for
the comonotone coupling; see, for instance, [29] for the expected shortfall
case (arguments for other $\alpha$ are similar).

######  Example 24.

If $b(x_{1},\dotsc,x_{d})=-\sum_{1\leq i<j\leq d}(x_{i}-x_{j})^{2}$ and if
$\alpha\equiv 1$, problem (13) is equivalent to the computation of Wasserstein
barycenters [1], while with $\alpha=\alpha_{m_{0}}$ we get the partial
Wasserstein barycenter problem [20]. In both cases, since $b$ is supermodular
the answer can be calculated explicitly by Theorem 22.

######  Example 25.

For $\alpha=1$, it is well known that optimal transport for the surplus
function in the preceding example is equivalent to the same problem with
$b(x_{1},\dotsc,x_{d})=|\sum_{i=1}^{d}x_{i}|^{2}$ (see, for example, [14]).
Since the mean $M:=\int_{X_{1}\times\dotsm\times
X_{d}}\sum_{i=1}^{d}x_{i}\mathrm{d}\gamma(x_{1},\dotsc,x_{d})=\sum_{i=1}^{d}\int_{X_{i}}x_{i}d\mu_{i}(x_{i})$
of the sum $z=\sum_{i=1}^{d}x_{i}$ is completely determined by the marginals
$\mu_{i}$ for any $\gamma\in\Gamma(\mu_{1},\dotsc,\mu_{d})$, (4) is equivalent
to maximizing the variance of the sum $z$ over
$\gamma\in\Gamma(\mu_{1},\dotsc,\mu_{d})$. Interpreting this sum as the losses
of a portfolio comprised of assets with losses $x_{i}$, we note the variance
is sometimes used in this setting as a measure of risk. One drawback to this
approach is its symmetry; variance punishes good outcomes (when the sum
$\sum_{i=1}^{d}x_{i}$ is below the mean) as well as bad ones (when the sum
$\sum_{i=1}^{d}x_{i}$ is above the mean).

Choosing $\alpha=\alpha_{m_{0}}$ in (11) on the other hand measures only
downside risk beyond a particular level. Since the quadratic surplus is
clearly supermodular, Theorem 22 implies that the worse case scenario can be
computed explicitly.

More generally, surplus functions of the form
$b(x_{1},\dotsc,x_{d})=H(\sum_{i=1}^{d}x_{i})$ for a convex function $H$ are
supermodular so Theorem 22 applies in this setting as well. Note that the case
$H(z)=\max(z-K,0)$ reflects the payout of a basket call option on the
underlying assets.

Our last example requires more extensive explanation and so we devote a
separate subsection to it.

###  4.1 Sensitivity analysis and maximal river flow

In this section we briefly describe a recurring example from [18], used
throughout that paper to illustrate issues in sensitivity analysis. In that
setting, one wants to understand the influence of the dependence structure
(among other factors) between several contributing inputs on an output
behavior.

We consider a simple model which involves the height of river at risk of
flooding and compares it to the height of a dyke which protects industrial
facilities. The maximal annual overflow $S$ of a river is modelled by

$S=Z_{\nu}+\Bigg{(}\dfrac{Q}{BK_{s}\sqrt{\frac{Z_{m}-Z_{\nu}}{L}}}\Bigg{)}^{0.6}-H_{d}-C_{b},$

where $\Bigg{(}\dfrac{Q}{BK_{s}\sqrt{\frac{Z_{m}-Z_{\nu}}{L}}}\Bigg{)}^{0.6}$
is the maximal annual height of the river and the variables
$Q,K_{s},Z_{\nu},Z_{m},H_{d},C_{b},L,B$ are physical quantities whose values
are modelled as random variables due to their variation in time and space,
measurement inaccuracies, or uncertainty of their true values. Their
individual distributions are modelled in [18] (see Table 1 on p.4), and are
all absolutely continuous with respect to Lebesgue measure. The $\alpha$-risk
(11), where the $x_{i}$ are the variables
$Q,K_{s},Z_{\nu},Z_{m},H_{d},C_{b},L,B$ and $b=S$ is the overflow, quantifies
the risk of overflow. In [18], the variables were assumed to be independent,
although other dependence structures are certainly possible; (13) asks what is
the maximal risk over all possible dependence structures. Notice then that the
surplus function $b$ is compatible and satisfies the strict monotonicity with
respect to each variable required in Theorem 22. The $(1,...,d)$ marginal
$\gamma$ of the $s$-comonotone solution $\pi$ defined in (7) is optimal, and
the unique optimizer on the support of $\alpha$, $\\{\alpha>0\\}$. In this
case, the $G_{i}$ corresponding to $Z_{\nu},Q$ and $L$ are monotone
increasing, while the $G_{i}$ corresponding to the other variables are
decreasing.

They can be retrieved explicitly by computing the quantile functions
associated to each marginals as described above.

##  5 Higher dimensional assets

There is an analogous theory of multi-marginal optimal transport when the
underlying variables lie in more general spaces. Although in these cases it is
generally not possible to derive explicit solutions as we did above, it is
possible to prove analogous structural properties of optimal couplings,
namely, that solutions are of Monge type (that is, concentrated on graphs over
$x_{1}$), for certain surplus functions $s$ (see [14] for an early result in
this direction, for a particular $s$, and [19] and [23] for general,
sufficient conditions on $s$).

We present one such result below, when each $X_{i}\subseteq\mathbb{R}^{n}$, to
illustrate how the theory can be adapted to the present setting. For
simplicity, we restrict our attention to the $d=2$ case. The conditions we
impose on $b$ are closely related to the local differential condition found in
[23]; however, the proof in [23] was restricted to the case when each marginal
was supported on a space of the same dimension, whereas here our first
marginal $\mu_{0}$, corresponding to $x_{0}=\alpha$, has one dimensional
support while the other marginals are supported in $\mathbb{R}^{n}$. We modify
that proof to fit the present setting here. Similar results can be established
for larger $d$, using similar arguments. However, the conditions imposed on
$b$ become more stringent for larger $d$, and more complicated to state.

In what follows, we will assume that $b$ is $\mathcal{C}^{2}$;
$D_{x_{i}}b(x_{1},x_{2})=(\frac{\partial b}{\partial x_{i}^{1}},\frac{\partial
b}{\partial x_{i}^{2}},...,\frac{\partial b}{\partial x_{i}^{n}})$ represents
the gradient of the function
$b:\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}$ with respect to
the variable $x_{i}\in\mathbb{R}^{n}$, for $i=1,2$. Similarly,
$D^{2}_{x_{1}x_{2}}b(x_{1},x_{2})=\big{(}\frac{\partial^{2}b}{\partial
x_{1}^{r}\partial x_{2}^{l}}\big{)}_{rl}$ is the $n$ by $n$ matrix of mixed
second order derivatives (each entry is a derivative with respect to one
coordinate from each of $x_{1}$ and $x_{2}$).

######  Theorem 26.

Suppose that $d=2$, that the domains $X_{1},X_{2}\subseteq\mathbb{R}^{n}$ are
compact and that $\mu_{1}$ is absolutely continuous with respect to Lebesgue
measure. Assume that $x_{2}\mapsto D_{x_{1}}b(x_{1},x_{2})$ is injective for
each fixed $x_{1}\in X_{1}$, and that for each $(x_{1},x_{2})\in X_{1}\times
X_{2}$, $\det(D^{2}_{x_{1}x_{2}}b(x_{1},x_{2}))\neq 0$ and

$D_{x_{2}}b(x_{1},x_{2})\cdot[D^{2}_{x_{1}x_{2}}b(x_{1},x_{2})]^{-1}D_{x_{1}}b(x_{1},x_{2})>0.$
(18)

Then the part of the solution to (13) away from $\alpha=0$ concentrates on the
graph of a function over $x_{1}$. Furthermore, if $|\\{\alpha=0\\}|=0$, the
solution is unique.

######  Remark 27.

The injectivity of $x_{2}\mapsto D_{x_{1}}b(x_{1},x_{2})$ is known as the
_twist condition_ in the optimal transport literature; it is well known that,
together with the absolute continuity of $\mu_{1}$, it guarantees the Monge
structure of the solution to the two marginal optimal transport problem with
surplus $b$ (see, for example, [33]). The invertibility of
$D^{2}_{x_{1}x_{2}}b(x_{1},x_{2})$ is frequently referred to as _non-
degeneracy_ , and can be seen as a linearized version of the twist. Here we
require additional hypotheses as we are dealing with the more sophisticated 3
marginal problem with surplus $s(\alpha,x_{1},x_{2})=\alpha b(x_{1},x_{2})$,
which is equivalent to the spectral risk maximization (13) by Theorem 18.

###### Proof.

By Theorem 18, we can consider the multi-marginal optimal transport problem
(4) rather than (13). Let $\pi$ solve (4) and $u_{0},u_{1},u_{2}$ the optimal
solutions to the dual problem (9). Then $u_{1}$ is differentiable $\mu_{1}$
almost everywhere, by a standard argument, originally found in [21]. Fix an
$x_{1}$ where $u_{1}$ is differentiable. To show that the solution
concentrates on a graph, we need to prove that there is only one $x_{2}$ and
one non-zero $\alpha$ such that $(\alpha,x_{1},x_{2})$ is in the support of
$\pi$.

At such points, the inequality constraint in (9), together with the equality
(10) on the support of the optimizer and the envelope theorem implies,
whenever $(\alpha,x_{1},x_{2})$ is in the support of the optimal $\pi$

$Du_{1}(x_{1})=\alpha D_{x_{1}}b(x_{1},x_{2}).$

Note that for fixed $x_{1}$, the twist condition implies that for $\alpha\neq
0$, this uniquely defines $x_{2}=x_{2}(\alpha)$ as a function of $\alpha$
(that is, for each $\alpha$, there is at most one $x_{2}$ satisfying this
equation). Furthermore, we can differentiate implicitly with respect to
$\alpha$ to obtain,

$0=D_{x_{1}}b(x_{1},x_{2}(\alpha))+D^{2}_{x_{1}x_{2}}b(x_{1},x_{2}(\alpha))D_{\alpha}x_{2}(\alpha),$

or

$D_{\alpha}x_{2}(\alpha)=-[D^{2}_{x_{1}x_{2}}b(x_{1},x_{2}(\alpha))]^{-1}D_{x_{1}}b(x_{1},x_{2}(\alpha)).$
(19)

Now, it is well known (see for example [14] or [23]) that the $u_{i}$ can be
taken to be $s$-conjugate. Explicitly, this means that $u_{0}$, the potential
corresponding to $x_{0}=\alpha$, can be written as

$u_{0}(\alpha)=\sup_{x_{1},x_{2}}[\alpha
b(x_{1},x_{2})-u_{1}(x_{1})-u_{2}(x_{2})],$

and is therefore convex, as a supremum of affine functions. If
$(\alpha,x_{1},x_{2})$ is in the support of the optimal $\pi$, by the
optimality condition presented in section 2, we must have that $x_{1},x_{2}$
attains the supremum above, so that $b(x_{1},x_{2})$ lies in the
subdifferential of $u_{0}$ at $\alpha$:

$b(x_{1},x_{2})\in\partial u_{0}(\alpha).$

But we also have $x_{2}=x_{2}(\alpha)$, so that

$b(x_{1},x_{2}(\alpha))\in\partial u_{0}(\alpha).$ (20)

We claim that this can hold for at most one value of $\alpha$. Indeed,
differentiating $b(x_{1},x_{2}(\alpha))$ with respect to $\alpha$ and using
(19) yields

$\begin{split}D_{\alpha}b(x_{1},x_{2}(\alpha))&=D_{x_{2}}b(x_{1},x_{2}(\alpha))\cdot
D_{\alpha}x_{2}(\alpha)\\\
&=-D_{x_{2}}b(x_{1},x_{2}(\alpha))\cdot[D^{2}_{x_{1}x_{2}}b(x_{1},x_{2}(\alpha))]^{-1}D_{x_{1}}b(x_{1},x_{2}(\alpha)).\end{split}$

Under the given assumption, this quantity is negative, and therefore the left
hand side in (20) is a strictly decreasing function of $\alpha$; it can
therefore intersect the subdifferential $\partial u_{0}(\alpha)$, an
increasing set valued function of $\alpha$, at most once. The points
$(\alpha,x_{2}(\alpha))$, where $\alpha$ is the point at which this
intersection occurs, are then the only points such that $(\alpha,x_{1},x_{2})$
are in the support of $\pi$, which establishes that the solution $\gamma$ to
(13) concentrates on a graph.

Uniqueness then follows by a very standard argument; if $\pi_{0}$ and
$\pi_{1}$ are both solutions, by linearity of the functional
$\gamma_{1/2}=\frac{1}{2}[\gamma_{0}+\gamma_{1}]$ is also a solution. The
argument above implies that the supports of $\gamma_{0}$ and $\gamma_{1}$
concentrate on graphs $T_{0}$ and $T_{1}$ over $x_{1}$. The support of
$\gamma_{1/2}$ then concentrates on the union of the graphs of $T_{0}$ and
$T_{1}$. However, since $\gamma_{1/2}$ is also a solution, the argument above
implies that this set should concentrate on a single graph. This is not
possible unless $T_{0}=T_{1}$, $\mu_{1}$ almost everywhere, in which case
$\gamma_{0}$ and $\gamma_{1}$ coincide, establishing the uniqueness assertion.
∎

There are of course output functions $b$ for which the hypotheses in the
preceding result do not hold, and in these cases solutions are often not of
Monge type. Alternatively, one can estimate the dimension of the support of
the optimal $\pi$ in (4) using the signature of the off diagonal part of the
Hessian of $s(\alpha,x_{1},...,x_{d})=\alpha b(x_{1},...,x_{d})$, as in [22]
and [24]. We illustrate this below with the following result in the $d=2$
case, which is an immediate consequence of Theorem 3.1.3 and Lemma 3.3.2 in
[22].

######  Proposition 28.

Assume $X_{1},X_{2}\subseteq\mathbb{R}^{n}$ and choose any point
$(\alpha,x_{1},x_{2})\in(0,\infty)\times X_{1}\times X_{2}$ where
$D^{2}_{x_{1}x_{2}}b(x_{1},x_{2})$ is non-singular. Let $\pi$ be optimal in
(4). Then in a neighbourhood of $(\alpha,x_{1},x_{2})$ the support of $\pi$ is
contained in a Lipschitz submanifold of dimension at most $n+1$. Furthermore,
if (18) holds, the dimension is at most $n$.

For larger $d$, if each $X_{i}\subseteq\mathbb{R}^{n_{i}}$ lies in a space of
dimension $n_{i}$, the results in [22] and [24] imply that the support of
$\pi$ will again be contained in a Lipschitz submanifold of dimension $m$,
where $m$ depends on the mixed second derivatives of $b$ in an intricate way.
Generically, if each $n_{i}=n$ is the same, $m$ will lie between $n$ and
$(d-1)n+1$. Of course, the dimension of the support of the solution $\gamma$
to (13) will be no larger than $m$, as by Theorem 18 $\gamma$ is the
$(1,...,d)$-marginal of $\pi$.

##  6 Stability

Among other results, in [15] the authors show stability of the optimal value
$\mathcal{R}_{\alpha}(b_{\\#}\gamma)$ and measure $\gamma$ with respect to
weak convergence of the marginals and $L^{1}$ convergence of the payoff
function $b$. Below, we show that similar results (under slightly different
assumptions – see Remark 30) can be easily deduced by combining Theorem 18
with optimal transport theory.

In addition, we further leverage the connection to optimal transport to
improve the stability with respect to the marginals when the underlying assets
are one dimensional and $b$ is compatible; in this case, we show that under
appropriate conditions, $\mathcal{R}_{\alpha}$ in fact exhibits Lipschitz
dependence on the marginals.

######  Proposition 29.

Suppose that for each $i$ $\\{\mu_{i}^{k}\\}$ is a sequence of probability
measures in $\mathcal{P}(X_{i})$ converging weakly to $\mu_{i}$ as
$k\rightarrow\infty$ and $\\{b^{k}\\}$ is a sequence of continuous output
functions on $X_{1}\times X_{2}\times\dotsm\times X_{d}$ converging uniformly
to a continuous $b$, such that $b$ is bounded above, that is $b\leq C$, and
bounded below by a sum of integrable functions, that is,
$b(x_{1},...,x_{d})\geq\sum_{i=1}^{d}u_{i}(x_{i})$ with each $u_{i}\in
L^{1}(\mu_{i})$. Let $\alpha^{k}$ be a sequence of spectral functions,
uniformly bounded $\alpha^{k}\leq M$, such that the associated spectral
measure $\alpha^{k}_{\\#}\mathrm{Leb}_{[0,1]}$ converges weakly to
$\alpha_{\\#}\mathrm{Leb}_{[0,1]}$, where $\alpha$ is bounded. Then

$\max_{\gamma\in\Gamma(\mu^{k}_{1},\mu^{k}_{2},\dotsc,\mu^{k}_{d})}\mathcal{R}_{\alpha^{k}}(b^{k}_{\\#}\gamma)\rightarrow\max_{\gamma\in\Gamma(\mu_{1},\mu_{2},\dotsc,\mu_{d})}\mathcal{R}_{\alpha}(b_{\\#}\gamma).$

Furthermore, if $\gamma^{k}$ maximizes
$\mathcal{R}_{\alpha^{k}}(b^{k}_{\\#}\gamma)$ over
$\Gamma(\mu^{k}_{1},\mu^{k}_{2},\dotsc,\mu^{k}_{d})$, then the weak limit
$\gamma$ of any weakly convergent subsequence of the $\gamma^{k}$ maximizes
$\mathcal{R}_{\alpha}(b_{\\#}\gamma)$ over
$\Gamma(\mu_{1},\mu_{2},\dotsc,\mu_{d})$.

######  Remark 30.

This result should be compared with Proposition 7 in [15]; in that setting,
stronger assumptions (uniform Holder continuity) on the $b^{k}$ are required,
although a weaker notion of convergence is imposed on them. In addition, we
also obtain stability with respect to variations in the spectral function
$\alpha$.

###### Proof.

The proof follows almost immediately from Theorem 18 and a corresponding
stability result in optimal transport; namely, [35, Theorem 5.20]. Note that
as we can take $X_{0}=[0,M]$ to be bounded, uniform convergence of the $b^{k}$
implies uniform convergence of the corresponding $s^{k}$ defined by (14).

The only issue to note is that that [35, Theorem 5.20] applies only to the two
marginal problem. However, exactly the same argument as is used there applies
to the multi-marginal case to show that the multi-marginal version of the
$s^{k}$ -cyclic monotonicity property satisfied by the support of each
$\pi^{k}$ implies that the support of the limit measure $\pi$ is
$s$-cyclically monotone. This then implies that $\pi$ is optimal in (4) by
[16, Theorem 12].

∎

We now specialize to the case when the $x_{i}$ are one dimensional. Our main
result in this direction requires the following lemma. It is in fact a special
case of [28, Corollary 11]; we provide the short proof for the convenience of
the reader.

######  Lemma 31.

If $\alpha\leq M$, then
$\mu\in\mathcal{P}(\mathbb{R})\mapsto\mathcal{R}_{\alpha}(\mu)$ is
$M$-Lipschitz for the $p$-Wasserstein distance,

$\left|\mathcal{R}_{\alpha}(\mu)-\mathcal{R}_{\alpha}(\nu)\right|\leq
M\mathrm{W}_{p}(\mu,\nu).$ (21)

###### Proof.

On $\mathbb{R}$, the $p$-Wasserstein distance is given by the $\mathrm{L}^{p}$
distance between the inverse cumulative distribution functions
$\mathrm{W}_{p}(\mu,\nu)=\|F_{\mu}^{-1}-F_{\nu}^{-1}\|_{\mathrm{L}^{p}([0,1])}.$
The statement follows:

$\displaystyle\mathcal{R}_{\alpha}(\mu)-\mathcal{R}_{\alpha}(\nu)$
$\displaystyle=\int(F_{\mu}^{-1}(m)-F_{\nu}^{-1}(m))\alpha(m)\mathrm{d}m$
$\displaystyle\leq\|F_{\mu}^{-1}-F_{\nu}^{-1}\|_{\mathrm{L}^{p}([0,1])}\|\alpha\|_{\infty}.$

∎

######  Proposition 32.

Assume that each $X_{i}\subset\mathbb{R}$ and that the output function $b$ is
weakly compatible, monotone increasing in each $x_{i}\in S_{+}$ and monotone
decreasing in each $x_{i}\in S_{-}$. If in addition $b$ is $K$-Lipschitz with
respect to $\|\cdot\|_{p}$ on $\mathbb{R}^{d}$ and $\alpha$ is bounded by $M$,
then:

$\left|\max_{\gamma\in\Gamma(\mu_{1},\dotsc,\mu_{d})}\mathcal{R}_{\alpha}(b_{\\#}\gamma)-\max_{\gamma\in\Gamma(\nu_{1},\dotsc,\nu_{d})}\mathcal{R}_{\alpha}(b_{\\#}\gamma)\right|\leq
MK\left(\sum_{i}\mathrm{W}_{p}(\mu_{i},\nu_{i})\right)^{1/p}$ (22)

###### Proof.

By Theorem 22, the $(1,\dotsc,d)-$ marginals
$\gamma=(G_{1},\dotsc,G_{d})_{\\#}\mathrm{Leb}_{[0,1]}$ and
$\tilde{\gamma}=(\tilde{G}_{1},\dotsc,\tilde{G}_{d})_{\\#}\mathrm{Leb}_{[0,1]}$
of the $s$-monotone couplings $\pi$ and $\tilde{\pi}$ of the $\mu_{i}$ and
$\nu_{i}$, respectively, achieve the suprema of $\mathcal{R}_{\alpha}$.
Therefore, we get that
$(G_{1},\dotsc,G_{d},\tilde{G}_{1},\dotsc,\tilde{G}_{d})_{\\#}\mathrm{Leb}_{[0,1]}\in\Gamma(\gamma,\tilde{\gamma})\subset\mathcal{P}(\mathbb{R}^{2d})$,
and its projection onto each $X_{i}\times X_{i}$ is
$(G_{i},\tilde{G}_{i})_{\\#}\mathrm{Leb}_{[0,1]}$, which is the optimal
coupling between $\mu_{i}$ and $\nu_{i}$ for the cost function
$|x_{i}-\tilde{x}_{i}|^{p}$ in the definition of the $p-$Wasserstein distance.
Therefore, setting $G=(G_{1},\dotsc,G_{d})$ and
$\tilde{G}=(\tilde{G}_{1},\dotsc,\tilde{G}_{d})$,

$\displaystyle\mathrm{W}_{p}^{p}(b_{\\#}\gamma,b_{\\#}\tilde{\gamma})$
$\displaystyle\leq\int_{0}^{1}(b(G(t))-b(\tilde{G}(t)))^{p}\mathrm{d}t$
$\displaystyle\leq K^{p}\int\|G(t)-\tilde{G}(t)\|_{p}^{p}\mathrm{d}t$
$\displaystyle=K^{p}\int\sum_{i}|G_{i}(t)-\tilde{G}_{i}(t)|^{p}\mathrm{d}t$
$\displaystyle=K^{p}\sum_{i}\mathrm{W}_{p}(\mu_{i},\nu_{i})^{p}$

One concludes using Lemma 31 on the Lipschitz continuity of
$\mathcal{R}_{\alpha}$. ∎

##  7 Multidimensional measures of risk

We consider the framework in [31] in which risk is measured in a multi-
dimensional way. Instead of a single output variable, we now have several,
depending on the same underlying random variables; this corresponds to a
vector valued output function $b:X_{1}\times\dotsm\times
X_{d}\to\mathbb{R}^{n}$. A natural form of multi-variate risk measures on the
distribution $b_{\\#}\gamma$ of output variables is then the maximal
correlation measure from [31], which is defined by

$\mathcal{R}_{\nu}(b_{\\#}\gamma)=\max_{\eta\in\Gamma(b_{\\#}\gamma,\nu)}\int_{\mathbb{R}^{n}\times\mathbb{R}^{n}}z\cdot
y\;\mathrm{d}\eta$

for some probability measure $\nu\in\mathcal{P}(\mathbb{R}^{n})$. We note that
it was proven in [11] that any strongly coherent multi-variate risk measure
takes this form. We consider the problem of maximizing
$\mathcal{R}_{\nu}(b_{\\#}\gamma)$ over all
$\gamma\in\Gamma(\mu_{1},\dotsc,\mu_{d})$, where the $\mu_{i}$ as before
represent the distributions of the underlying variables. Exactly as in Theorem
18, one can show that this problem is equivalent to the multi-marginal problem

$\max_{\pi\in\Gamma(\nu,\mu_{1},\dotsc,\mu_{d})}\int
b(x_{1},\dotsc,x_{d})\cdot y\;\mathrm{d}\pi.$ (23)

This problem is more challenging than the case of a scalar valued $b$;
nonetheless, we are able to obtain the following results in particular cases.

######  Proposition 33.

Suppose that the underlying variables $x_{i}$ are one dimensional. Assume that
$\nu$ is concentrated on a smooth curve, that is
$\nu=f_{\\#}\mathrm{Leb}_{[0,1]}$ with $f:[0,1]\to\mathbb{R}^{d}$, where each
component of $f$ is positive and monotone increasing, and each component
$b_{j}$ of $b=(b_{1},b_{2},\dotsc,b_{n})$ is supermodular and monotone
increasing in each $x_{i}$. Then
$\Big{(}(m,x_{1},\dotsc,x_{d})\mapsto(f(m),x_{1},\dotsc,x_{d})\Big{)}_{\\#}\pi$
is optimal in (23),where
$\pi=(\mathrm{id},F_{\mu_{1}}^{-1},F_{\mu_{2}}^{-1},\dotsc,F_{\mu_{d}}^{-1})_{\\#}\mathrm{Leb}_{[0,1]}$
is the comonotone coupling of $\mu_{0},\mu_{1},\dotsc,\mu_{d}$, where
$\mu_{0}=\mathrm{Leb}_{[0,1]}$.

###### Proof.

We rewrite the problem (23) above as

$\max_{\pi\in\Gamma(\nu,\mu_{1},\dotsc,\mu_{d})}\int
b(x_{1},\dotsc,x_{d})\cdot
y\;\mathrm{d}\pi=\max_{\pi\in\Gamma(\mathrm{Leb}_{[0,1]},\mu_{1},\dotsc,\mu_{d})}\int
b(x_{1},\dotsc,x_{d})\cdot f(m)\;\mathrm{d}\pi.$

The objective function $b(x_{1},\dotsc,x_{d})\cdot
f(m)=\sum_{j=1}^{d}b_{j}(x_{1},\dotsc,x_{d})f_{j}(m)$ is now a supermodular
function of the one dimensional variables $m,x_{1},x_{2},...,x_{d}$, and so
Proposition 9 implies the desired result. ∎

For more general, diffuse $\nu$, we are able to prove that the solution is of
Monge form and unique, provided that $b$ is invertible.

######  Proposition 34.

Assume that $\nu$ is absolutely continuous with respect to $n$ dimensional
Lebesgue measure and $b:X_{1}\times\dotsm\times X_{d}\to\mathbb{R}^{n}$ is
invertible. Then there exists a unique solution to (23). Furthermore, it
concentrates on a graph over $y$.

Note that the invertibility assumption on $b$ implies that the sum of the
dimensions of the $X_{i}$ must be less than or equal to the dimension $n$ of
$\nu$

###### Proof.

Noting that the gradient of the surplus function in (23) with respect to $y$
is $b$, the invertibility condition on $b$ is the famous twist condition in
optimal transport, and so this result follows by a minor refinement of a
standard argument which we produce below.

Letting $u_{0}(y),u_{1}(x_{1}),\dotsc,u_{n}(x_{n})$ be the optimal solutions
to the dual problem, we note that as $\sum_{i=0}^{d}u_{i}(x_{i})\geq
b(x_{1},\dotsc,x_{d})\cdot y$, with equality by (10) in the support of any
optimzer $\pi$, the envelope theorem implies,

$Du_{0}(y)=D_{y}(y\cdot b(x_{1},\dotsc,x_{d}))=b(x_{1},\dotsc,x_{d})$ (24)

wherever $u_{0}$ is differentiable, at points $(y,x_{1},\dotsc,x_{d})$ in the
support of $\pi$. A now standard argument, found in [21], implies that $u_{0}$
is Lipschitz and therefore differentiable Lebesgue almost everywhere by
Rademacher’s theorem; the differentiability holds $\nu$ almost everywhere by
the absolute continuity of $\nu$. Therefore, (24) holds on a set of full $\pi$
measure. Invertibility of $b$ then means that this equation becomes
$(x_{1},\dotsc,x_{d})=b^{-1}(Du_{0}(y))$. Setting $T(y)=b^{-1}(Du_{0}(y))$,
this means that $\pi$ concentrates on the graph of
$T:\mathbb{R}^{n}\rightarrow\mathbb{R}^{d}$, as desired.

Uniqueness then follows exactly as in the proof of Theorem 26. ∎

## References

  * [1] M. Agueh and G. Carlier.  Barycenters in the Wasserstein space.  SIAM J. on Mathematical Analysis, 43(2):904–924, 2011.
  * [2] Jean-David Benamou.  Optimal transportation, modelling and numerical simulation.  Acta Numer., 30:249–325, 2021.
  * [3] Jean-David Benamou, Guillaume Carlier, Marco Cuturi, Luca Nenna, and Gabriel Peyré.  Iterative Bregman projections for regularized transportation problems.  SIAM J. Sci. Comput., 37(2):A1111–A1138, 2015.
  * [4] Yann Brenier.  The least action principle and the related concept of generalized flows for incompressible perfect fluids.  Journal of the American Mathematical Society, 2(2):225–255, 1989\.
  * [5] Giuseppe Buttazzo, Luigi De Pascale, and Paola Gori-Giorgi.  Optimal-transport formulation of electronic density-functional theory.  Physical Review A, 85(6):062502, 2012.
  * [6] Luis A. Caffarelli and Robert J. McCann.  Free boundaries in optimal transport and Monge-Ampère obstacle problems.  Ann. of Math. (2), 171(2):673–730, 2010.
  * [7] Guillaume Carlier.  On a class of multidimensional optimal transportation problems.  Journal of convex analysis, 10(2):517–530, 2003.
  * [8] Guillaume Carlier and Bruno Nazaret.  Optimal transportation for the determinant.  ESAIM Control Optim. Calc. Var., 14(4):678–698, 2008.
  * [9] Codina Cotar, Gero Friesecke, and Claudia Klüppelberg.  Density functional theory and optimal transportation with coulomb cost.  Communications on Pure and Applied Mathematics, 66(4):548–599, 2013\.
  * [10] I. Ekeland.  An optimal matching problem.  ESAIM Control. Optim. Calc. Var., 11:57–71, 2005.
  * [11] Ivar Ekeland, Alfred Galichon, and Marc Henry.  Comonotonic measures of multivariate risks.  Mathematical Finance: An International Journal of Mathematics, Statistics and Financial Economics, 22(1):109–132, 2012.
  * [12] Alessio Figalli.  The optimal partial transport problem.  Arch. Ration. Mech. Anal., 195(2):533–560, 2010.
  * [13] Hans Föllmer and Alexander Schied.  Stochastic finance, volume 27 of De Gruyter Studies in Mathematics.  Walter de Gruyter & Co., Berlin, 2002.  An introduction in discrete time.
  * [14] Wilfrid Gangbo and Andrzej Swiech.  Optimal maps for the multidimensional monge-kantorovich problem.  Communications on pure and applied mathematics, 51(1):23–45, 1998\.
  * [15] Mario Ghossoub, Jesse Hall, and David Saunders.  Maximum spectral measures of risk with given risk factor marginal distributions.  To appear in Mathematics of Operations Research. Preprint available at arXiv:2010.14673, 2020.
  * [16] Claus Griessler.  $C$-cyclical monotonicity as a sufficient criterion for optimality in the multimarginal Monge-Kantorovich problem.  Proc. Amer. Math. Soc., 146(11):4735–4740, 2018.
  * [17] Noureddine Igbida and Van Thanh Nguyen.  Augmented Lagrangian method for optimal partial transportation.  IMA J. Numer. Anal., 38(1):156–183, 2018.
  * [18] Bertrand Iooss and Paul Lemaître.  A review on global sensitivity analysis methods.  Uncertainty management in simulation-optimization of complex systems, pages 101–122, 2015.
  * [19] Young-Heon Kim and Brendan Pass.  A general condition for monge solutions in the multi-marginal optimal transport problem.  SIAM Journal on Mathematical Analysis, 46(2):1538–1550, 2014.
  * [20] Jun Kitagawa and Brendan Pass.  The multi-marginal optimal partial transport problem.  arXiv preprint arXiv:1401.7255, 2014.
  * [21] Robert J. McCann.  Polar factorization of maps on Riemannian manifolds.  Geom. Funct. Anal., 11(3):589–608, 2001.
  * [22] Brendan Pass.  Structural results on optimal transportation plans.  PhD thesis, University of Toronto, 2011.
  * [23] Brendan Pass.  Uniqueness and monge solutions in the multimarginal optimal transportation problem.  SIAM Journal on Mathematical Analysis, 43(6):2758–2775, 2011.
  * [24] Brendan Pass.  On the local structure of optimal measures in the multi-marginal optimal transportation problem.  Calculus of Variations and Partial Differential Equations, 43(3-4):529–536, 2012.
  * [25] Brendan Pass.  Regularity of optimal transportation between spaces with different dimensions.  Math. Res. Lett., 19(2):291–307, 2012.
  * [26] Brendan Pass.  Multi-marginal optimal transport: theory and applications.  arXiv preprint arXiv:1406.0026, 2014.
  * [27] Brendan Pass and Adolfo Vargas-Jiménez.  Monge solutions and uniqueness in multi-marginal optimal transport: weaker conditions on the cost, stronger conditions on the marginals.  arXiv preprint arXiv:2202.06783, 2022.
  * [28] Alois Pichler.  Evaluations of risk measures for different probability measures.  SIAM J. Optim., 23(1):530–551, 2013.
  * [29] Giovanni Puccetti.  Sharp bounds on the expected shortfall for a sum of dependent random variables.  Statistics & Probability Letters, 83(4):1227–1232, 2013.
  * [30] Giovanni Puccetti and Ludger Rüschendorf.  Sharp bounds for sums of dependent risks.  J. Appl. Probab., 50(1):42–53, 2013.
  * [31] Ludger Rüschendorf.  Law invariant convex risk measures for portfolio vectors.  Statist. Decisions, 24(1):97–108, 2006.
  * [32] Ludger Rüschendorf.  Mathematical risk analysis.  Springer Series in Operations Research and Financial Engineering. Springer, Heidelberg, 2013.  Dependence, risk bounds, optimal allocations and portfolios.
  * [33] Filippo Santambrogio.  Optimal transport for applied mathematicians.  Progress in Nonlinear Differential Equations and their Applications, 87\. Birkhäuser/Springer, Cham, 2015.  Calculus of variations, PDEs, and modeling.
  * [34] C. Villani.  Topics in Optimal Transportation, volume 58 of Graduate Studies in Mathematics.  American Mathematical Society, Providence, 2003.
  * [35] C. Villani.  Optimal Transport: Old and New, volume 338 of Grundlehren der mathematischen Wissenschaften.  Springer, 2009.
  * [36] Bin Wang and Ruodu Wang.  The complete mixability and convex minimization problems with monotone marginal densities.  J. Multivariate Anal., 102(10):1344–1360, 2011.


11institutetext: Zama 111 11email<EMAIL_ADDRESS>http://zama.ai

# Deep Neural Networks for Encrypted Inference with TFHE

Andrei Stoian  Jordan Frery  Roman Bredehoft  Luis Montero  Celia Kherfallah
Benoit Chevallier-Mames

###### Abstract

Fully homomorphic encryption (FHE) is an encryption method that allows to
perform computation on encrypted data, without decryption. FHE preserves the
privacy of the users of online services that handle sensitive data, such as
health data, biometrics, credit scores and other personal information. A
common way to provide a valuable service on such data is through machine
learning and, at this time, Neural Networks are the dominant machine learning
model for unstructured data.

In this work we show how to construct Deep Neural Networks (DNN) that are
compatible with the constraints of TFHE, an FHE scheme that allows arbitrary
depth computation circuits. We discuss the constraints and show the
architecture of DNNs for two computer vision tasks. We benchmark the
architectures using the Concrete stack222https://github.com/zama-ai/concrete-
ml, an open-source implementation of TFHE.

##  1 Introduction

Neural Networks (NNs) are machine learning (ML) models that have driven the
recent expansion of the field of Artificial Intelligence (AI). Their
performance on unstructured data such as images, sound and text is unmatched
by other ML techniques. Moreover, deep NNs obviate the need for complex
feature engineering and process raw data directly, making them easier to
deploy in production. Applications of NNs include image classification, face
recognition, voice assistants, and search engines, tools which today are a
staple of the user experience online. Deployment of such models in SaaS
applications raises a security risk: they are a target of malevolent entities
that seek to steal the sensitive user data these models process.

Privacy-preserving technologies, such as multi-party computing (MPC) and fully
homomorphic encryption (FHE), provide a solution to the risk of data leaks,
eliminating it by design. Notably, FHE encrypts user data and allows a third
party to process the data in its encrypted form, without needing to decrypt
it. Only the data owner can decrypt the result of the computation. Thus, an
attacker can only steal encrypted data they can not decrypt.

In this work we show how to build neural networks that are FHE compatible,
while minimizing the cryptography knowledge needed by the machine learning
practitioner. We based our work on the Concrete Library [7] which uses TFHE
[6], works over integers, provides a fast _programmable_ bootstrapping
mechanism, and performs exact computation.

##  2 Related work

Several alternative approaches exist for neural network inference over
encrypted data. All use NNs with integer weights and activations and many of
them rely on "leveled" fully homomorphic encryption schemes that do not use
bootstrapping, such as CKKS [5] and YASHE [3].

CryptoNets [9] uses YASHE which supports the computation of polynomials of
encrypted values. CryptoNets are NNs quantized to integers (of 5-10 bits) with
activation functions expressed as low-degree polynomials. CryptoNets achieve
99% accuracy on MNIST using a three layer network with an inference time of
570 seconds/image.

FHE-DiNN [4], a TFHE based approach, quantizes inputs, intermediate values and
weights to binary values. In this case, the training is done with hardSigmoid
activation which is swapped for the sign function in inference. However,
binary NNs are hard to train and do not perform well in many ML tasks such as
object detection and speech processing.

Another TFHE approach, SHE [11], uses bit series representation of encrypted
values and boolean gates. They run NNs that fit within a maximum
multiplicative depth budget and, by avoiding expensive multi-bit PBSs, they
achieve inference of a ShuffleNet on ImageNet with a latency of 18 000
seconds/image. They rely on logarithmic quantization of weights which allows
to reduce multiplicative depth for the convolution layers by using bit-shifts.
Sums, relu and maxpool are computed using boolean gates.

Leveled approaches such as SHE and CryptoNets are limited by the maximum
multiplicative depth budget, which, in turn, limits the supported network
types and their depth. Moreover, some schemes such as CKKS are approximate by
design, as the noise corrupts some of the message bits.

In this work we propose an approach to train arbitrary NNs which can have any
depth, number of neurons and activation functions. Furthermore, our approach
performs exact computation in FHE: the noise of the encryption scheme does not
corrupt the values that are processed. Thus results in FHE are the same as in
the clear - there is no degradation of accuracy when moving to encrypted
inference - which is a major advantage when putting models in production.

##  3 Neural Network Training for Encrypted Inference

Training NNs is usually done in floating point, but most FHE schemes,
including TFHE, only support integers. Consequently, quantization must be
used, and two main approaches exist:

  1. 1.

Post-training quantization is commonly used [9, 11], but, in this mode, NNs
lose accuracy when the quantization bit-width is lower than 7-8 bits. With
per-channel quantization, or logarithmic quantization, which are more complex
to implement, as few as 4 bits were used for weights and activations without
loss of accuracy [14].

  2. 2.

Quantization-aware Training (QAT), used in this work and in [4], is an
approach that adds quantizers to network activations and weights during
training. QAT enables extreme quantization with less than 4 bit weights and
activations.

To support arbitrarily deep NNs and any activation function, we make use of
the programmable boostrapping mechanism [8] (PBS) of TFHE. PBS reduces the
noise in accumulators of ciphertext leveled operations (addition,
multiplication with clear constants) but also allows to apply a lookup-table
(TLU) on its input ciphertext.

The TFHE PBS mechanism has a rather high computational cost, and this cost
depends on the number of bits of the encrypted value to be boostrapped. It is
convenient to keep the accumulator size low, in order to speed up the PBS
computation. However, reducing accumulator bit-width has a negative impact on
network prediction performance, so a compromise needs to be found.

We describe here a QAT strategy that can process all the intermediate
encrypted values as integers. In this way, training an FHE compatible network
becomes purely a machine learning problem and no cryptography knowledge is
needed by the practitioner. To build a TFHE compatible NN, the constraints on
the network architecture are the following:

  * •

All layers that sum or multiply two encrypted values, such as convolution conv
and fully-connected fc, must have quantized inputs. This is easily achieved
using QAT frameworks.

  * •

The bit-width of the accumulators of layers such as conv, fc must be bounded.
To achieve this, we use pruning.

To control the accumulator bit-width while keeping the training dynamics
stable, we use $L^{1}$-norm unstructured pruning. Figure 1 shows the impact of
pruning on the accumulator size for two quantization modes: narrow and wide
range.

While the inputs of conv and fc layers need to be quantized, it is possible to
use floating point layers for all univariate operations such as batch
normalization, quantization, and activations.

In our FHE compatible NNs the outputs of a conv or fc are processed by a
sequence of univariate operations that ends with quantization. This sequence
of functions takes integers and has integer outputs, but the intermediary
computations in these operations can use float parameters. Thus, batch
normalization, activation functions, neuron biases and any other univariate
transformation of conv or fc outputs does not need quantization. Figure 1
shows the architecture of the network during training and inference.

Figure 1: Left: accumulator size while varying the number of active neurons
during pruning for a 3-layer fully-connected network with 2 bit weights and
activations. Two quantization modes are shown: Narrow range uses values
$[-2^{b-1}+1,2^{b-1}-1]$, while Wide range uses $[-2^{b-1},2^{b-1}-1]$. Right:
the structure of a 2 layer convolutional network in training and during
inference. Univariate layers are fused to table-lookups, implemented with PBS.

##  4 Neural Network Inference using TFHE

Inference of our FHE NNs is based on quantized implementations of NN operators
that add or multiply together encrypted values. Convolutional, fully connected
and average pooling layers use the quantized formulation from [10]. Since
uniform quantization is used, we can define a quantized value $r$ as
$r=S(q-Z)$ where $S$ is the quantization scale, $Z$ is the quantization zero-
point and $q$ is the integer representation of the value. Next, the fully
connected layer, with inputs $x$, weights $w$ and outputs $o$, with per-tensor
quantization parameters $(S_{x},Z_{x})$, $(S_{w},Z_{w})$ can be written as:

$S_{o}(q_{o}^{k}-Z_{o})=\sum_{i=0}^{N}S_{x}(q_{x}^{i}-Z_{x})S_{w}(q_{w}^{(i,k)}-Z_{w})+b^{k}$
(1)

where $k$ is the index of a neuron in the layer and $N$ is the number of
connections of the neuron and $b^{k}$ is the bias of the $k$-th neuron. A
convolutional layer can be expressed by extending the sum to the height, width
and channel dimensions. Equation 2 can be re-written to separate integer and
floating point computations (note that zero-points $Z_{x},Z_{o},Z_{w}$ are
integers).

$q_{o}^{k}=b^{k}+Z_{o}+\frac{S_{x}S_{w}}{S_{o}}\sum_{i=0}^{N}(q_{x}^{i}-Z_{x})(q_{w}^{(i,k)}-Z_{w})$
(2)

Therefore we can separate the equation in a floating point univariate function
$f$ and a sum over products of encrypted inputs and clear weights:

$q_{o}^{k}=f(\Sigma)~{}~{}\textrm{where}~{}~{}f(q)=b^{k}+Z_{o}+\frac{S_{x}S_{w}}{S_{o}}q~{}~{}\textrm{and}~{}~{}\Sigma=\sum_{i=0}^{N}(q_{x}^{i}-Z_{x})(q_{w}^{(i,k)}-Z_{w})$
(3)

The univariate function $f$ in eq. 3 takes integer inputs. We compose this
function with the batch-normalization, and, finally, with the quantization
function $Q(x)=floor\left(\frac{x}{S_{x}}\right)+Z_{x}$. Thus $f$ becomes a
function defined on $\mathbb{Z}$, with values in $\mathbb{Z}$ and can be
implemented as a lookup table with a PBS in FHE, without any loss of
precision.

The complete NN computation can now be expressed over integers using the
following operations: multiplication of an encrypted value and a clear
constant, sums of encrypted integer values, table lookup of encrypted integer
values. In our implementation of TFHE, Concrete, we encode integers in two
different ways: integers up to 8 bits are encoded into a single ciphertext,
and integers between 9-16 bits are encoded with a CRT representation into
several ciphertexts as described in [2]. This contrasts to previous works,
such as [11], that encode each bit of an integer as an individual ciphertext
and use boolean gates to build arithmetic circuits.

An automated optimization process [2] determines the cryptographic parameters
of the circuit, based on several factors: (1) the _circuit bit-width_ ,
defined as the minimum bit-width necessary to encode the largest integer value
obtained anywhere in the NN’s integer-based evaluation, (2) the maximum 2-norm
of the integer weight tensors of the layers, and (3), the desired probability
of error of the PBS. The optimization process determines the cryptosystem
parameters (LWE dimension, polynomial size, GLWE dimension, etc.) to ensure a
fast execution, the target probability of failure and the security level
(using the lattice-estimator [1]). We set the PBS error probability
sufficiently low to ensure full correctness of the results, i.e. the results
in the clear are always the same as those in FHE, up to a user-defined error-
rate, e.g. $10^{-6}$, for one full NN inference.

##  5 Experimental Results

The networks were implemented in PyTorch with Brevitas [13] and converted to
FHE with Concrete-ML [12]. We ran experiments on two datasets with several
neural network architectures, in two quantization modes (see Figure. 1, left).
The test machine had an Intel i7-11800H CPU with 8 cores and we used 16
threads for the experiments.

Table 1: Experimental results obtained with Brevitas and Concrete-ML Network |  Quant. bits |  Active Neurons |  Narrow range |  Data-set |  Circuit bit-width |  Accuracy |  Inference time (s)
---|---|---|---|---|---|---|---
3-layer FCNN333Three FC layers with 192, 192 and 10 neurons |  2/2 |  150 |  Yes |  MNIST |  6 |  92.2% |  31
3-layer FCNN |  2/2 |  90 |  No |  MNIST |  7 |  96.5% |  77
3-layer FCNN |  2/2 |  190 |  No |  MNIST |  8 |  97.1% |  300
LeNet |  2/2 |  190 |  No |  MNIST |  8 |  97.6% |  2780
6-layer CNN444Four conv layers with 8,8,16 and 16 filters followed by a FC layer with 120 neurons and a final FC layer for classification |  2/2 |  190 |  No |  MNIST |  8 |  98.7% |  5072
VGG-9555Six conv layers with: 64, 64, 128, 128, 256, 256 filters, followed by two 512 neuron FC layers and a final FC layer for classification |  2/2666Inputs are quantized in 8 bits, but all other activations use 2 bits |  all |  Yes |  CIFAR10 |  13 |  87.5% |  18000777Estimated time for a 8 core machine, using 16 threads

##  6 Conclusion

Our approach to encrypted inference for Neural Networks shows several
advantages over other methods. First, we believe our method is easier to use
than other works, since the problem of making an FHE compatible network
becomes strictly an ML problem and no cryptography knowledge is needed.
Second, the computations in FHE are correct with respect to the computations
in the clear and, using TFHE, noise does not corrupt the encrypted values.
Thus, once a network is trained incorporating the quantization constraints,
the accuracy that is measured on clear data will be the same as that on
encrypted data. Finally, our approach, using PBS, shows competitive accuracies
in FHE and allows to convert arbitrary depth networks using any activation
function to FHE. Networks up to 9 layers were shown, but deeper NNs can easily
be implemented.

Preliminary code for the MNIST classifier is
available888https://github.com/zama-ai/concrete-
ml/tree/release/0.5.x/use_case_examples and code for the CIFAR10 classifier
will be released soon.

Many possible strategies can be employed to improve upon this work, in order
to support larger models, such as ResNet, on larger data-sets like ImageNet.
For example, a better pruning strategy could decrease the PBS count, per-
channel quantization can improve accuracy, and faster step functions in FHE
could improve the overall speed.

## References

  * [1] Albrecht, M.R., Player, R., Scott, S.: On the concrete hardness of learning with errors. Cryptology ePrint Archive, Paper 2015/046 (2015), https://eprint.iacr.org/2015/046, https://eprint.iacr.org/2015/046
  * [2] Bergerat, L., Boudi, A., Bourgerie, Q., Chillotti, I., Ligier, D., Orfila, J.B., Tap, S.: Parameter optimization & larger precision for (t)fhe. Cryptology ePrint Archive, Paper 2022/704 (2022), https://eprint.iacr.org/2022/704, https://eprint.iacr.org/2022/704
  * [3] Bos, J.W., Lauter, K., Loftus, J., Naehrig, M.: Improved security for a ring-based fully homomorphic encryption scheme. In: Stam, M. (ed.) Cryptography and Coding. pp. 45–64. Springer Berlin Heidelberg, Berlin, Heidelberg (2013)
  * [4] Bourse, F., Minelli, M., Minihold, M., Paillier, P.: Fast homomorphic evaluation of deep discretized neural networks. p. 483–512. Springer-Verlag, Berlin, Heidelberg (2018). https://doi.org/10.1007/978-3-319-96878-0_17
  * [5] Cheon, J.H., Kim, A., Kim, M., Song, Y.: Homomorphic encryption for arithmetic of approximate numbers. In: Takagi, T., Peyrin, T. (eds.) Advances in Cryptology – ASIACRYPT 2017. pp. 409–437. Springer International Publishing, Cham (2017)
  * [6] Chillotti, I., Gama, N., Georgieva, M., Izabachène, M.: Tfhe: Fast fully homomorphic encryption over the torus. Journal of Cryptology 33, 34–91 (2019)
  * [7] Chillotti, I., Joye, M., Ligier, D., Orfila, J.B., Tap, S.: Concrete: Concrete operates on ciphertexts rapidly by extending tfhe. In: WAHC 2020–8th Workshop on Encrypted Computing & Applied Homomorphic Cryptography. vol. 15 (2020)
  * [8] Chillotti, I., Joye, M., Paillier, P.: Programmable bootstrapping enables efficient homomorphic inference of deep neural networks. IACR Cryptol. ePrint Arch. 2021, 91 (2021)
  * [9] Gilad-Bachrach, R., Dowlin, N., Laine, K., Lauter, K., Naehrig, M., Wernsing, J.: Cryptonets: Applying neural networks to encrypted data with high throughput and accuracy. In: Balcan, M.F., Weinberger, K.Q. (eds.) Proceedings of The 33rd International Conference on Machine Learning. Proceedings of Machine Learning Research, vol. 48, pp. 201–210. PMLR, New York, New York, USA (20–22 Jun 2016), https://proceedings.mlr.press/v48/gilad-bachrach16.html
  * [10] Jacob, B., Kligys, S., Chen, B., Zhu, M., Tang, M., Howard, A., Adam, H., Kalenichenko, D.: Quantization and training of neural networks for efficient integer-arithmetic-only inference. In: 2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 2704–2713 (2018). https://doi.org/10.1109/CVPR.2018.00286
  * [11] Lou, Q., Jiang, L.: She: A fast and accurate deep neural network for encrypted data. In: Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., Garnett, R. (eds.) Advances in Neural Information Processing Systems. vol. 32. Curran Associates, Inc. (2019), https://proceedings.neurips.cc/paper/2019/file/56a3107cad6611c8337ee36d178ca129-Paper.pdf
  * [12] Meyre, A., Chevallier-Mames, B., Frery, J., Stoian, A., Bredehoft, R., Montero, L., Kherfallah, C.: Concrete-ml (2022), https://github.com/zama-ai/concrete-ml
  * [13] Pappalardo, A.: Xilinx/brevitas (2021). https://doi.org/10.5281/zenodo.3333552, https://doi.org/10.5281/zenodo.3333552
  * [14] Yvinec, E., Dapogny, A., Cord, M., Bailly, K.: SPIQ: data-free per-channel static input quantization. CoRR abs/2203.14642 (2022). https://doi.org/10.48550/arXiv.2203.14642, https://doi.org/10.48550/arXiv.2203.14642