# Kohler-Jobin meets Ehrhard: the sharp lower bound for the Gaussian principal
frequency while the Gaussian torsional rigidity is fixed, via rearrangements

Orli Herscovici, Galyna V. Livshyts  School of Mathematics, Georgia Institute
of Technology, Atlanta, GA<EMAIL_ADDRESS>School of Mathematics,
Georgia Institute of Technology, Atlanta, GA<EMAIL_ADDRESS>
###### Abstract.

In this note, we provide an adaptation of the Kohler-Jobin rearrangement
technique to the setting of the Gauss space. As a result, we prove the
Gaussian analogue of the Kohler-Jobin’s resolution of a conjecture of Pólya-
Szegö: when the Gaussian torsional rigidity of a (convex) domain is fixed, the
Gaussian principal frequency is minimized for the half-space. At the core of
this rearrangement technique is the idea of considering a “modified” torsional
rigidity, with respect to a given function, and rearranging its layers to
half-spaces, in a particular way; the Rayleigh quotient decreases with this
procedure.

We emphasize that the analogy of the Gaussian case with the Lebesgue case is
not to be expected here, as in addition to some soft symmetrization ideas, the
argument relies on the properties of some special functions; the fact that
this analogy does hold is somewhat of a miracle.

##  1\. Introduction

Consider a log-concave measure $\mu$ on $\mathbb{R}^{n}$ with density
$e^{-V}$, for some convex function $V$, and its associated Laplacian

$L\cdot=\Delta\cdot-\langle\nabla\cdot,\nabla V\rangle.$

Let $K\subset\mathbb{R}^{n}$ be a simply-connected domain. Let
$W^{1,2}(K,\mu)$ stand for the weighted Sobolev space and $\nabla$ for the
weak gradient. Letting $W_{0}^{1,2}(K,\mu)=W^{1,2}(K,\mu)\cap\\{u:u|_{\partial
K}=0\\}$ (where the boundary value is understood in the sense of the trace),
define the $\mu$-torsional rigidity as

$T_{\mu}(K)=\sup_{u\in
W_{0}^{1,2}(K,\mu)}\frac{(\int_{K}u\,d\mu)^{2}}{\int_{K}|\nabla u|^{2}d\mu}.$

See section 2 for the details, or, for example, Pólya and Szegö [18].

The $\mu$-principal frequency of a domain $K$ is defined to be

$\Lambda_{\mu}(K):=\inf_{u\in W_{0}^{1,2}(K,\mu)}\frac{\int_{K}|\nabla
u|^{2}d\mu}{\int_{K}u^{2}d\mu}.$

Note that the torsional rigidity is monotone increasing while the principle
frequency is monotone decreasing, i.e. whenever $K\subset L,$ we have
$T_{\mu}(K)\leq T_{\mu}(L)$ and $\Lambda_{\mu}(K)\geq\Lambda_{\mu}(L)$.

In the case when $\mu$ is the Lebesgue measure and $L=\Delta,$ these
quantities have been studied extensively, and are intimately tied with the
subject of isoperimetric inequalities. See, e.g. Kawohl [10], Pólya and Szegö
[18], Burchard [4], Lieb and Loss [16], Kesavan [11], or Vazquez [20]. In
particular, the Faber-Krahn inequality [9, 14, 15] states that the (Lebesgue)
principal frequency of a domain $K$ of a fixed Lebesgue measure is _minimized_
when $K$ is an Euclidean ball. The result of Saint-Venant (see e.g. [18])
states that, conversely, the torsional rigidity of a domain $K$ of a fixed
Lebesgue measure is _maximized_ when $K$ is an Euclidean ball.

The easiest way to prove these results is via rearrangements. For a set $K$ in
$\mathbb{R}^{n}$, denote by $K^{*}$ the Euclidean ball of the same Lebesgue
measure as $K$. Recall that the Schwartz rearrangement of a non-negative
function $u:K\rightarrow\mathbb{R}$ is the function
$u^{*}:K^{*}\rightarrow\mathbb{R}$ whose level sets $\\{u^{*}\geq t\\}$ are
all Euclidean balls, and such that $|\\{u\geq t\\}|=|\\{u^{*}\geq t\\}|$
(where $|\cdot|$ stands for the Lebesgue measure.) The Pólya-Szegö principle
[18] (which is a consequence of the isoperimetric inequality) implies that
$\int_{K}|\nabla u|^{2}dx\geq\int_{K^{*}}|\nabla u^{*}|^{2}dx$, while the
definition of the symmetrization yields that
$\int_{K}u^{2}dx=\int_{K^{*}}(u^{*})^{2}dx$, and, for a non-sign-changing
function $u,$ $\int_{K}u\,dx=\int_{K^{*}}u^{*}dx$. Therefore, the Faber-Krahn
and the Saint-Venant results follow (together with some additional information
that the extremal function for the torsional rigidity is non-negative, as can
be deduced via the maximal principle).

In the case of the Gaussian measure $\gamma$ (which is the measure with
density $\frac{1}{(\sqrt{2\pi})^{n}}\,e^{-\frac{x^{2}}{2}}$), the analogue of
the Schwartz rearrangement was developed by Ehrhard [7, 8]. The Euclidean
balls are replaced with the Gaussian isoperimetric regions, which are half-
spaces (see Sudakov and Tsirelson [19], Borel [5]). In the Gaussian world,
$K^{*}$ is the half-space of the same Gaussian measure as the domain $K$, and
$u^{*}$ is the function whose level sets are half-spaces, and such that
$\gamma(\\{u\geq t\\})=\gamma(\\{u^{*}\geq t\\}).$ The Pólya-Szegö principle
is replaced with the analogous Ehrhard principle, which yields
$\int_{K}|\nabla u|^{2}d\gamma\geq\int_{K^{*}}|\nabla u^{*}|^{2}d\gamma.$ As a
result, whenever the Gaussian measure of the domain $K$ is fixed, the Gaussian
principal frequency $\Lambda_{\gamma}(K)$ is minimized when $K$ is a half-
space (see Carlen and Kerce [6]), and the Gaussian torsional rigidity
$T_{\gamma}(K)$ is maximized when $K$ is a half-space (see e.g. Livshyts
[17]).

In the Lebesgue world, Pólya and Szegö asked another natural question: if for
a set $K,$ not its measure, but its torsional rigidity is fixed, then is the
principal frequency still minimized on the Euclidean ball? This question was
answered in the affirmative by Kohler-Jobin [12, 13] back in the 1970s. For
further generalizations and applications, as well as a nice exposition of the
topic, see Brasco [2].

In this paper, we develop the Gaussian analogue of the Kohler-Jobin
rearrangement, and show

######  Theorem 1.1.

For any convex domain $K\subset\mathbb{R}^{n}$, letting $H$ be the half-space
such that $T_{\gamma}(K)=T_{\gamma}(H)$, we have
$\Lambda_{\gamma}(K)\geq\Lambda_{\gamma}(H).$

We would like to emphasize that the analogy of our work and the works of
Kohler-Jobin [12, 13] and Brasco [2] is not to be expected, like in the case
of the inequalities of Faber-Krahn and Saint-Venant! Our proof relies heavily
on the particular properties of certain special functions, and not just on
soft properties of rearrangements. It is somewhat miraculous that Theorem 1.1
is true!

This paper is organized as follows. In Section 2 we discuss some
preliminaries. In Section 3 we discuss the Gaussian modified torsional
rigidity. In Section 4 we define the Gaussian version of the Kohler-Jobin
rearrangement, study its properties, and prove Theorem 1.1.

Acknowledgement. The authors are supported by the NSF CAREER DMS-1753260.

##  2\. Preliminaries

###  2.1. General

Fix a measurable connected set $K$ in $\mathbb{R}^{n}$. Let $\gamma$ denote
the Gaussian measure, that is the measure with density
$\frac{1}{(\sqrt{2\pi})^{n}}\,e^{-\frac{x^{2}}{2}}.$ Let $W^{1,2}(K,\gamma)$
be the Sobolev space of square integrable functions with respect to $\gamma,$
whose (weak) gradient is also square integrable, and let
$W^{1,2}_{0}(K,\gamma)$ denote the space of all $v(x)\in W^{1,2}(K,\gamma)$
vanishing on $\partial K$.

For a point $x\in\partial K,$ we denote by $n_{x}$ the unit outer normal
vector to the boundary of $K$ at $x$; if $K$ is convex then $n_{x}$ is defined
almost everywhere on the boundary of $K$.

Recall the co-area formula (see e.g. Brothers and Zeimer [3] or [21]) which
states that

$\int_{\mathbb{R}^{n}}|\nabla
u|f\,dx=\int_{0}^{\infty}\int_{\\{u=s\\}}f\,d\mathcal{H}^{n-1}ds,$

where $f\in W^{1,p}(\mathbb{R}^{n})$, $1\leq p\leq\infty$, is a compactly
supported function, $u\in C^{1}(\mathbb{R}^{n})$ is a nonnegative Borel
function, and $\mathcal{H}^{n-1}$ denotes the $n-1$ dimensional Hausdorff
measure. By replacing $f(x)$ with $f(x)e^{-\frac{x^{2}}{2}}$, we get

(1) $\int_{\mathbb{R}^{n}}|\nabla
u|f\,d\gamma=\int_{0}^{\infty}\int_{\\{u=s\\}}f\,d\gamma_{\partial\\{u=s\\}}ds,$

where by $d\gamma_{\partial M}$, for a surface $M,$ we denote the measure
$e^{-\frac{x^{2}}{2}}d\mathcal{H}^{n-1},$ i.e. the boundary measure with the
Gaussian weight.

Recall the Ornstein-Uhlenbeck operator [1] defined as

(2) $Lu=\Delta u-\langle x,\nabla u\rangle,$

which verifies the integration by parts identity for the standard Gaussian
measure:

$\int_{K}vLu\,d\gamma=-\int_{K}\langle\nabla u,\nabla
v\rangle\,d\gamma+\int_{\partial K}v\langle\nabla
u,n_{x}\rangle\,d\gamma_{\partial K}.$

Recall also the following classical existence result, see e.g. [17].

######  Theorem 2.1 (Dirichlet boundary condition).

Let $K$ be a Lipschitz domain. Let $F\in L^{2}(K,\gamma)$. Then there exists a
unique function $u\in W^{1,2}(K,\gamma)$ which is a solution of

$\begin{cases}Lu=F&\,on\,\,\,K,\\\ u=0&\,on\,\,\,\partial K.\end{cases}$

Moreover, if $\partial K$ is $C^{2}$ then $u\in C^{2}(int(K))\cap
C^{1}(\bar{K}).$

Furthermore, if $F\leq 0$ then $u\geq 0,$ and the level sets of $u$ are
nested.

Note that the “furthermore” part of the theorem follows from the maximal
principle.

###  2.2. Equivalent definitions of the Gaussian torsional rigidity

######  Definition 2.2.

For a convex domain $K\subset\mathbb{R}^{n}$, define the Gaussian torsional
rigidity by

(3) $T_{\gamma}(K)=\sup_{v\in
W^{1,2}_{0}(K,\gamma)}\frac{(\int_{K}v\,d\gamma)^{2}}{\int_{K}|\nabla
v|^{2}d\gamma}.$

This object (with a different normalization!) was considered in [17]. Let us
also state

######  Definition 2.3.

For a convex domain $K\subset\mathbb{R}^{n}$, define the Gaussian torsion
function of $K$ to be the unique function $v\in W^{1,2}(K)$ satisfying the
differential equation

$\begin{cases}Lv=-1&\,on\,\,\,K,\\\ v=0&\,on\,\,\,\partial K.\end{cases}$

We note that by the maximal principle, the Gaussian torsion function is
necessarily non-negative.

The Gaussian torsional rigidity admits several equivalent definitions.

######  Proposition 2.4.

The following are equivalent:

  1. (1)

$T_{\gamma}(K)=\inf_{u\in W^{1,2}(K,\gamma):\,\,Lu=-1}\int_{K}|\nabla
u|^{2}d\gamma;$

  2. (2)

$T_{\gamma}(K)$ is the Gaussian torsional rigidity, i.e.
$T_{\gamma}(K)=\sup_{\begin{subarray}{c}v\in
W^{1,2}_{0}(K,\gamma)\end{subarray}}\frac{(\int_{K}v\,d\gamma)^{2}}{\int_{K}|\nabla
v|^{2}d\gamma};$

  3. (3)

$T_{\gamma}(K)=\sup_{\begin{subarray}{c}v\in
W^{1,2}_{0}(K,\gamma)\end{subarray}}\left({-}\int_{K}|\nabla
v|^{2}d\gamma{+2}\int_{K}v\,d\gamma\right);$

  4. (4)

$T_{\gamma}(K)=\int_{K}vd\gamma=\int_{K}|\nabla
v|^{2}d\gamma={-}\int_{K}|\nabla v|^{2}d\gamma{+2}\int_{K}v\,d\gamma,$ where
$v$ is the Gaussian torsion function of $K.$

###### Proof.

For any $v\in W^{1,2}_{0}(K,\gamma)$ and any $u\in W^{1,2}(K,\gamma)$ with
$Lu=-1$, we have, by Cauchy’s inequality,

$\int_{K}|\nabla u|^{2}d\gamma\geq\frac{\left(\int_{K}\langle\nabla u,\nabla
v\rangle\,d\gamma\right)^{2}}{\int_{K}|\nabla
v|^{2}d\gamma}=\frac{\left(\int_{K}v\,d\gamma\right)^{2}}{\int_{K}|\nabla
v|^{2}d\gamma},$

where the last identity is obtained integrating by parts. Therefore,

$\inf_{\begin{subarray}{c}u\in W^{1,2}(K,\gamma):\\\
Lu=-1\end{subarray}}\int_{K}|\nabla
u|^{2}d\gamma\geq\sup_{\begin{subarray}{c}v\in
W^{1,2}_{0}(K,\gamma)\end{subarray}}\frac{(\int_{K}v\,d\gamma)^{2}}{\int_{K}|\nabla
v|^{2}d\gamma}.$

In order to see that the equality is attained, consider (on both sides) the
function $v\in W^{1,2}_{0}(K,\gamma)$ satisfying the differential equation

$\begin{cases}Lv=-1&\,on\,\,\,K,\\\ v=0&\,on\,\,\,\partial K.\end{cases}$

Thus (1) and (2) are equivalent.

In order to see that (1) and (3) are equivalent, use, in place of Cauchy’s
inequality, the fact that for any $v\in W^{1,2}_{0}(K,\gamma)$ and any $u\in
W^{1,2}(K,\gamma)$ with $Lu=-1$,

$\int_{K}|\nabla u|^{2}d\gamma\geq\int_{K}\left(2\langle\nabla u,\nabla
v\rangle-|\nabla v|^{2}\right)d\gamma,$

and argue in the same manner as before.

The equivalence of (4) to the rest of the definitions follows automatically by
the choice of the extremizing function. Note that with this choice of $v,$ we
have, integrating by parts, $\int_{K}v\,d\gamma=\int_{K}|\nabla
v|^{2}d\gamma.$ ∎

We shall record the following classical fact:

######  Lemma 2.5.

[torsional rigidity is monotone] For a pair of convex sets $K$ and $L,$
$T_{\gamma}(K)\geq T_{\gamma}(L)$ whenever $L\subset K.$

###### Proof.

Let $v$ be the torsion function on $L$ (i.e. the function $v$ maximizing the
righthand side of the (3)), and extend it by zero to the whole of $K$; let us
call the resulting function $\tilde{v}.$ Note that the result is still a
Sobolev function since $v|_{\partial L}=0.$ We have

$T_{\gamma}(L)=\frac{(\int_{L}v\,d\gamma)^{2}}{\int_{L}|\nabla
v|^{2}d\gamma}=\frac{(\int_{K}\tilde{v}\,d\gamma)^{2}}{\int_{K}|\nabla\tilde{v}|^{2}d\gamma}\leq
T_{\gamma}(K),$

where in the last line, the definition of the torsional rigidity was used
again. ∎

Let us also recall the Gaussian analogue of the Saint-Venant theorem (see e.g.
Proposition 5.6 in Livshyts [17]):

######  Lemma 2.6.

Of all the measurable sets $K$ with Gaussian measure $a\in[0,1]$, the
torsional rigidity is maximizied for the half-space
$\\{x_{1}\leq\Phi^{-1}(a)\\}.$

Here the function $\Phi^{-1}(t)$ is the inverse function of the standard
Gaussian distribution
$\Phi(t)=\int\limits_{-\infty}^{t}\frac{1}{\sqrt{2\pi}}\exp(-\frac{s^{2}}{2})\,ds$.

###  2.3. Gaussian torsional rigidity of a half-space

Let us denote the right half-space

$H_{s}=\\{x_{1}\geq s\\}.$

Define the function $T:\mathbb{R}\rightarrow\mathbb{R}$ to be the Gaussian
torsional rigidity of the half-space:

$T(s)=T_{\gamma}(H_{s}).$

We shall show

######  Lemma 2.7.

$T^{\prime}(s)=-\sqrt{2\pi}\,e^{\frac{s^{2}}{2}}\gamma(H_{s})^{2}.$

###### Proof.

Note that the torsion function $v_{s}:H_{s}\rightarrow\mathbb{R}$ is given by

$v_{s}(x)=V(x_{1})-V(s),$

where $V(t)$ is such a function that

$V^{\prime}(t)=e^{\frac{t^{2}}{2}}\int_{t}^{\infty}e^{-\frac{\tau^{2}}{2}}d\tau.$

Indeed, $v_{s}\in W^{1,2}(H_{s})$, and one may check directly that
$Lv_{s}=-1$, and clearly $v_{s}(x)=0$ when $x\in\partial H_{s}.$ Therefore, by
(4) of Proposition 2.4,

$T(s)=T_{\gamma}(H_{s})=\int_{H_{s}}|\nabla
v_{s}|^{2}d\gamma=\frac{1}{\sqrt{2\pi}}\int_{s}^{\infty}V^{\prime}(t)^{2}\,e^{-\frac{t^{2}}{2}}dt,$

where in the last passage we used Fubini’s theorem. Thus

$T^{\prime}(s)=-\frac{1}{\sqrt{2\pi}}\,V^{\prime}(s)^{2}\,e^{-\frac{s^{2}}{2}}=-\sqrt{2\pi}\,e^{\frac{s^{2}}{2}}\left(\frac{1}{\sqrt{2\pi}}\int_{s}^{\infty}e^{-\frac{\tau^{2}}{2}}d\tau\right)^{2},$

which finishes the proof. ∎

###  2.4. Equivalent definitions of the Gaussian principal frequency

######  Definition 2.8.

For a convex domain $K\subset\mathbb{R}^{n}$, define the Gaussian principal
frequency $\Lambda_{\gamma}(K)$ by

$\Lambda_{\gamma}(K)=\inf_{\begin{subarray}{c}v\in
W^{1,2}_{0}(K,\gamma)\end{subarray}}\frac{\int_{K}|\nabla
v|^{2}d\gamma}{\int_{K}v^{2}d\gamma}.$

We outline the following classical fact (see e.g. Pólya and Szegö for the
Lebesgue version [18].)

######  Proposition 2.9.

$\Lambda_{\gamma}(K)$ is the first eigenvalue of $L$ with the Dirichlet
boundary condition on $K$.

###### Proof.

We aim to show that $\Lambda_{\gamma}(K)$ is the smallest positive number
$\Lambda_{\gamma}$ such that there exists a function $u\in W_{0}^{1,2}(K)$
with

(4) $Lu=-\Lambda_{\gamma}u.$

Let $v\in W_{0}^{1,2}(K)$ be a nonnegative function satisfying the
differential equation (4) in the interior of the domain $K$ and $v|_{\partial
K}=0$.

First of all we assume that $f=qv$ for a smooth function $q$ defined in
$K\cup\partial K$. Then

$\displaystyle\int_{K}\left(|\nabla
f|^{2}-\Lambda_{\gamma}f^{2}\right)d\gamma$
$\displaystyle=\int_{K}\left(\left<\nabla(qv),\nabla(qv)\right>-\Lambda_{\gamma}q^{2}v^{2}\right)d\gamma$
$\displaystyle=\int_{K}\left(q^{2}|\nabla v|^{2}+2vq\left<\nabla q,\nabla
v\right>+v^{2}|\nabla q|^{2}+q^{2}vLv\right)d\gamma$
$\displaystyle=\int_{K}|\nabla q|^{2}v^{2}d\gamma+F,$

where

(5) $F=\int_{K}\left(q^{2}|\nabla v|^{2}+2vq\left<\nabla q,\nabla
v\right>+q^{2}vLv\right)d\gamma$

Integrating by parts and following the definition (2), we get

$\displaystyle\int_{K}2vq\left<\nabla q,\nabla v\right>d\gamma$
$\displaystyle=-\frac{1}{2}\int_{K}\left<\nabla q^{2},\nabla
v^{2}\right>d\gamma$ $\displaystyle=-\frac{1}{2}\int_{K}q^{2}Lv^{2}d\gamma$
$\displaystyle=-\frac{1}{2}\int_{K}q^{2}(2vLv+2|\nabla v|^{2})d\gamma.$

Substituting it into (5), we obtain

$F=\int_{K}\left(q^{2}|\nabla v|^{2}-q^{2}vLv-q^{2}|\nabla
v|^{2}+q^{2}vLv\right)d\gamma=0,$

which implies

$\int_{K}|\nabla
f|^{2}d\gamma=\Lambda_{\gamma}\int_{K}f^{2}d\gamma+\int_{K}|\nabla
q|^{2}v^{2}d\gamma,$

and therefore

$\frac{\int_{K}|\nabla
f|^{2}d\gamma}{\int_{K}f^{2}d\gamma}\geq\Lambda_{\gamma}.$

The equality is attained if $q$ is a constant function, and the statement
follows by approximation. ∎

We shall record the following classical fact.

######  Lemma 2.10.

[principal frequency is monotone decreasing] For a pair of convex sets $K$ and
$L,$ $\Lambda_{\gamma}(K)\geq\Lambda_{\gamma}(L)$ whenever $K\subset L.$

The proof of this Lemma is completely analogous to the proof of Lemma 2.5,
except we consider infimum in place of supremum, and argue that infimum over a
larger class is smaller than the infimum over a smaller class. We leave the
details to the reader.

##  3\. Modified torsional rigidity

Let us use the notation

(6) $T_{\gamma}[v]=-\int_{K}|\nabla v|^{2}d\gamma+2\int_{K}v\,d\gamma,\quad
v\in W^{1,2}_{0}(K,\gamma).$

It was shown in the Proposition 2.4 that

$T_{\gamma}(K)=\sup_{v\in W^{1,2}_{0}(K,\gamma)}T_{\gamma}[v].$

Consider a non-negative function $u$ on $K$ which vanishes at the boundary and
let $u_{\max}$ be the maximum value of $u$ on K.

######  Definition 3.1.

Define $Cl(u)$ to be the class of functions $v(x)\in W^{1,2}_{0}(K,\gamma)$
which can be presented as $v=\varphi\circ u(x)$, for some real-valued function
$\varphi(x)$ defined on the interval $[0,u_{\max}]$ and satisfying
$\varphi(0)=0$.

The _Gaussian modified torsional rigidity_ $T^{mod}_{\gamma}(K;u(x))$ with
respect to a reference function $u\in W_{0}^{1,2}(K,\gamma)$ is given by

(7) $T^{mod}_{\gamma}(K;u)=\max_{v\in Cl(u)}T_{\gamma}[v].$

Note that $T_{\gamma}[v]=T_{\gamma}[\varphi\circ u(x)]$, and thus
$T^{mod}_{\gamma}(K;u),$ depends only on the collection of the level sets of
$u(x)$.

Fix a non-negative $u(x)\in W_{0}^{1,2}(K,\gamma)$ whose level sets are
nested. Let us define

$\displaystyle K_{t}$ $\displaystyle=\\{x\in K\,|\,u(x)>t\,\\},$
$\displaystyle\partial K_{t}$ $\displaystyle=\\{x\in K\,|\,u(x)=t\,\\},$

and

$\displaystyle\gamma(K_{t})$ $\displaystyle=\int_{K_{t}}d\gamma,$ (8)
$\displaystyle\ell(t)$ $\displaystyle=\int_{\partial K_{t}}|\nabla
u|\,d\gamma_{\partial K_{t}}.$

We shall show the following analogue of the [12, Lemma 1], [13, Lemma 1] (see
also a wonderful exposition by Brasco [2]), in which we are going to find
exactly the maximizing function for the modified torsional rigidity, given an
arbitrary reference function $u(x)\in W_{0}^{1,2}(K,\gamma)$.

######  Proposition 3.2.

Fix a non-negative reference function $u\in W_{0}^{1,2}(K,\gamma)$ with nested
level sets. The Gaussian modified torsional rigidity of $K$ with respect to
$u$ is given by

$T^{mod}_{\gamma}(K;u(x))=\int_{0}^{u_{\max}}\frac{\gamma(K_{t})^{2}}{\ell(t)}\,dt.$

###### Proof.

Suppose $v=\varphi\circ u,$ for some function
$\varphi:\mathbb{R}\rightarrow\mathbb{R}$ such that $\varphi(0)=0.$ Suppose
without loss of generality that $\varphi\in C^{2}(\mathbb{R})$ (otherwise we
may approximate it by a smooth function.)

Consider the first summand from the functional $T_{\gamma}[v]$ defined in (6):

(9) $\int_{K}|\nabla
v|^{2}d\gamma=\int_{K}|\nabla\varphi(u(x))|^{2}d\gamma=\int_{K}(\varphi^{\prime})^{2}|\nabla
u|^{2}d\gamma.$

By (9) and by the co-area formula (1),

(10) $\int_{K}|\nabla
v|^{2}d\gamma=\int_{0}^{u_{\max}}\left|\frac{d\varphi}{dt}\right|^{2}\int_{\partial
K_{t}}|\nabla u|\,d\gamma_{\partial
K_{t}}\,dt=\int_{0}^{u_{\max}}\left|\frac{d\varphi}{dt}\right|^{2}\ell({t})\,dt.$

Next, again using the co-area formula (1), we get

(11) $\frac{d\gamma(K_{t})}{dt}=\frac{d}{dt}\left(\int_{K_{t}}\frac{1}{|\nabla
u|}|\nabla u|\,d\gamma\right)=$
$\frac{d}{dt}\left(\int_{t}^{\infty}\int_{\partial K_{s}}\frac{1}{|\nabla
u|}\,d\gamma_{\partial K_{s}}ds\right)=-\int_{\partial K_{t}}\frac{1}{|\nabla
u|}\,d\gamma_{\partial K_{t}}.$

Therefore, combining (11) with (yet again!) the co-area formula (1), we see

$\displaystyle\int_{K}v(x)\,d\gamma$
$\displaystyle=\int_{K}\varphi(u(x))\,d\gamma=\int_{0}^{u_{\max}}\varphi(t)\int_{\partial
K_{t}}\frac{1}{|\nabla u|}d\gamma_{\partial K_{t}}dt$ (12)
$\displaystyle=\int_{0}^{u_{\max}}\varphi(t)\left(-\frac{d\gamma(K_{t})}{dt}\right)dt=\int_{0}^{u_{\max}}\gamma(K_{t})\left(\frac{d\varphi}{dt}\right)dt,$

where the last equality follows by the one-dimensional integration by parts,
together with the facts that $\gamma(K_{u_{max}})=0$ and $\varphi(0)=0$.

Combining (10) and (12), we write the functional $T_{\gamma}[v]$ (defined in
(6)) with $v=\varphi\circ u$ as follows:

$T_{\gamma}[\varphi(u)]=-\int_{0}^{u_{\max}}\left|\frac{d\varphi}{dt}\right|^{2}\ell(t)\,dt+2\int_{0}^{u_{\max}}\gamma(K_{t})\,\frac{d\varphi}{dt}\,dt.$

Since $\ell(t)>0$ for $t<u_{\max}$, we can rewrite the last equation as

$\displaystyle T_{\gamma}[v]$
$\displaystyle=\int_{0}^{u_{\max}}\ell(t)\left[2\,\frac{d\varphi}{dt}\frac{\gamma(K_{t})}{\ell(t)}-\left|\frac{d\varphi}{dt}\right|^{2}\right]dt$
$\displaystyle=\int_{0}^{u_{\max}}\ell(t)\left[-\left(\frac{d\varphi}{dt}-\frac{\gamma(K_{t})}{\ell(t)}\right)^{2}+\left(\frac{\gamma(K_{t})}{\ell(t)}\right)^{2}\right]dt$
(13)
$\displaystyle\leq\int_{0}^{u_{\max}}\frac{\gamma^{2}(K_{t})}{\ell(t)}\,dt.$

Considering the function
$\varphi_{0}(t)=\int_{0}^{t}\frac{\gamma(K_{s})}{\ell(s)}\,ds,$ we have
$d\varphi_{0}/dt=\gamma(K_{t})/\ell(t)$, and the equality is attained in (13);
therefore,

$T_{\gamma}(K;u)=\sup_{\varphi}\,T_{\gamma}[\varphi(u)]=T_{\gamma}[\varphi_{0}(u)]=\int_{0}^{u_{\max}}\frac{\gamma^{2}(K_{t})}{\ell(t)}\,dt,$

which completes the proof. ∎

######  Corollary 3.3.

The functional $T_{\gamma}(K_{t};u(x)-t)$ is the Gaussian modified torsional
rigidity of a level set $K_{t}$ with respect to a function $u(x)-t$, and it is
given by

$T_{\gamma}(K_{t};u(x)-t)=\int_{t}^{u_{\max}}\frac{\gamma^{2}(K_{s})}{\ell(s)}\,ds.$

##  4\. The Kohler-Jobin-style symmetrization in the Gauss space, and the
proof of Theorem 1.1

Let $K$ be a domain in $\mathbb{R}^{n}$. Let $u\in W_{0}^{1,2}(K)$ be a non-
negative function.

In this section, we define the analogue of the Kohler-Jobin symmetrization in
the Gauss space. The general idea is to correspond the level sets of $u$ to
half-spaces whose torsional rigidity is the same as the modified torsional
rigidity of the level sets of $u.$ Now, these half-spaces are level sets of
$u^{\dagger}$, but composed with a specific function. Let us describe the
details below.

In view of Corollary 3.3, define the distribution function of the modified
torsional rigidity of $K$ with respect to $u$ by

(14)
$D(t)=T_{\gamma}(K_{t};u(x)-t)=\int_{t}^{u_{\max}}\frac{\gamma^{2}(K_{s})}{\ell(s)}\,ds.$

By Corollary 3.3,

(15)
$(D^{-1})^{\prime}(\tau)=\frac{1}{D^{\prime}(D^{-1}(\tau))}=-\frac{\ell(D^{-1}(\tau))}{\gamma^{2}(K_{D^{-1}(\tau)})}.$

Next, recall from the sub-section 2.3 the notation $H_{s}=\\{x_{1}\geq s\\}$
and

$T(s)=T_{\gamma}(H_{s}).$

Let $K^{\dagger}$ be the half-space of the form $\\{x_{1}\geq s\\}$ such that

$t_{0}:=T^{mod}_{\gamma}(K;u)=T_{\gamma}(K^{\dagger}).$

We shall define the _Ehrhrard-Kohler-Jobin rearrangement_ of $u$ to be the
non-negative function $u^{\dagger}\in W_{0}^{1,2}(K^{\dagger})$, non-
decreasing in $x_{1}$, and constant in $x_{i}$ for $i\neq 1,$ given by

$u^{\dagger}(x)=f\circ T(x_{1}),$

where $f:\mathbb{R}\rightarrow\mathbb{R}$ is the (non-increasing) function
given by

(16)
$\begin{cases}f^{\prime}(\tau)=\frac{(D^{-1})^{\prime}(\tau)\gamma(K_{D^{-1}(\tau)})}{\gamma(H_{T^{-1}(\tau)})},\\\
f(t_{0})=0.\end{cases}$

In other words, $u^{\dagger}=f(\tau)$ on the boundary of the right half-space
whose torsional rigidity is $\tau$, and so we have, for all $\tau\geq 0,$

(17) $T_{\gamma}(K_{D^{-1}(\tau)})\geq
T^{mod}_{\gamma}(K_{D^{-1}(\tau)};u-D^{-1}(\tau))=T_{\gamma}(H_{T^{-1}(\tau)}).$

(Just to clarify, the last identity is a tautology $\tau=\tau$). In view of
(17) together with Lemma 2.6 (the Gaussian analogue of Saint-Venant theorem),
we conclude that

(18) $\gamma(K_{D^{-1}(\tau)})\geq\gamma(H_{T^{-1}(\tau)}).$

Note also that the condition $f(t_{0})=0$ means that $u^{\dagger}=0$ on the
boundary of $K^{\dagger}$, just as intended.

We are ready to prove

######  Theorem 4.1.

The function $u^{\dagger}$ has the following properties:

  1. (1)

$\int_{K}|\nabla u|^{2}\,d\gamma=\int_{K^{\dagger}}|\nabla
u^{\dagger}|^{2}\,d\gamma;$

  2. (2)

$\int_{K}u\,d\gamma=\int_{K^{\dagger}}u^{\dagger}\,d\gamma;$

  3. (3)

For any convex non-decreasing function
$F:\mathbb{R}^{+}\rightarrow\mathbb{R}^{+}$ we have

$\int_{K}F(u)\,d\gamma\leq\int_{K^{\dagger}}F(u^{\dagger})\,d\gamma.$

###### Proof.

Part (1). By the co-area formula (1), combined with the definition (8) of the
function $l(s),$ we have

(19) $\int_{K}|\nabla u|^{2}d\gamma=\int_{0}^{u_{max}}l(s)\,ds.$

Applying change of variables $\tau=D(s),$ with the function $D$ defined in
(14), we write

(20)
$\int\limits_{0}^{u_{max}}l(s)\,ds=-\int\limits_{0}^{t_{0}}(D^{-1})^{\prime}(\tau)\cdot
l\circ
D^{-1}(\tau)d\tau=\int\limits_{0}^{t_{0}}\left((D^{-1})^{\prime}(\tau)\right)^{2}\gamma(K_{D^{-1}(\tau)})^{2}d\tau,$

where in the last passage we used (15). Combining (19) and (20), we get

(21) $\int_{K}|\nabla
u|^{2}d\gamma=\int_{0}^{t_{0}}\left((D^{-1})^{\prime}(\tau)\right)^{2}\gamma(K_{D^{-1}(\tau)})^{2}d\tau.$

Next, for any $\tau\in\mathbb{R},$ in view of the definition of $u^{\dagger},$
we have

(22) $\int\limits_{\\{x_{1}=T^{-1}(\tau)\\}}|\nabla
u^{\dagger}|d\gamma_{\partial
H_{T^{-1}(\tau)}}=\frac{1}{\sqrt{2\pi}}\,e^{-\frac{T^{-1}(\tau)^{2}}{2}}f^{\prime}(\tau)T^{\prime}(T^{-1}(\tau))=-f^{\prime}(\tau)\gamma(H_{T^{-1}(\tau)})^{2},$

where in the last passage we used the fact that, by Lemma 2.7,
$T^{\prime}(s)=-\sqrt{2\pi}e^{\frac{s^{2}}{2}}\gamma(H_{s})^{2}$.

Thus we write, by the change of variables $s=f(\tau),$ and again by the co-
area formula:

(23) $\int_{K^{\dagger}}|\nabla
u^{\dagger}|^{2}d\gamma=\int_{0}^{u_{\max}^{\dagger}}\int_{\\{u^{\dagger}=s\\}}|\nabla
u^{\dagger}|d\gamma_{\partial\\{u^{\dagger}=s\\}}ds=$
$-\int_{0}^{t_{0}}f^{\prime}(\tau)\int_{\\{u^{\dagger}=f(\tau)\\}}|\nabla
u^{\dagger}|d\gamma_{\partial\\{u^{\dagger}=f(\tau)\\}}d\tau=\int_{0}^{t_{0}}f^{\prime}(\tau)^{2}\gamma(H_{T^{-1}(\tau)})^{2}d\tau,$

where in the last passage we used (22) together with the fact that
$u^{\dagger}(x)=f(\tau)$ whenever $x_{1}=T^{-1}(\tau).$ Finally, by (21), (23)
and (16), we conclude part (1).

Part (2). We write the layer-cake representation and apply the change of
variables $\tau=D(t)$ and $\tau=f^{-1}(t)$:

$\int_{K}u\,d\gamma=\int_{0}^{u_{max}}\gamma(\\{u\geq
t\\})\,dt=-\int_{0}^{t_{0}}\gamma(K_{D^{-1}(\tau)})(D^{-1})^{\prime}(\tau)\,d\tau=$
$-\int_{0}^{t_{0}}\gamma(H_{T^{-1}(\tau)})f^{\prime}(\tau)\,d\tau=\int_{K^{\dagger}}u^{\dagger}\,d\gamma,$

where in the last passage we used the definition of $f$ and $u^{\dagger}$
(16).

Part (3). Similarly to part (2), we use the layer cake formula together with
the additional change of variables $s=F^{-1}(t)$ and write

$\int\limits_{K}F(u)\,d\gamma=\int\limits_{F(0)}^{F(u_{max})}\gamma(\\{F(u)\geq
t\\})\,dt=-\int\limits_{0}^{t_{0}}F^{\prime}(D^{-1}(\tau))\gamma(K_{D^{-1}(\tau)})(D^{-1})^{\prime}(\tau)\,d\tau$

and

$\int\limits_{K^{\dagger}}F(u^{\dagger})\,d\gamma=\int\limits_{F(0)}^{F(u^{\dagger}_{max})}\gamma(\\{F(u^{\dagger})\geq
t\\})dt=-\int\limits_{0}^{t_{0}}F^{\prime}(f(\tau))\gamma(H_{T^{-1}(\tau)})f^{\prime}(\tau)d\tau.$

In view of the fact that $F$ is convex, and therefore $F^{\prime}$ is non-
decreasing, and in view of (16), in order to conclude part (3), we are only
left to show that for all $\tau\in[0,t_{0}],$

(24) $D^{-1}(\tau)\leq f(\tau).$

Note that by (16) and (18), in view of the fact that both $D$ and $f$ are non-
increasing, we have

$-(D^{-1})^{\prime}(\tau)\leq-f^{\prime}(\tau).$

We conclude that

$D^{-1}(\tau)=-\int_{\tau}^{t_{0}}(D^{-1})^{\prime}(\tau)\,d\tau\leq-\int_{\tau}^{t_{0}}f^{\prime}(\tau)\,d\tau=f(\tau),$

and (24) follows. ∎

######  Remark 4.2.

Note that there is not much room in this construction, and one really has to
be very precise when constructing the symmetrization. Indeed, if, say,
$\int_{K}|\nabla u|^{2}d\gamma>\int_{K^{\dagger}}|\nabla
u^{\dagger}|^{2}d\gamma$ (rather than equals), but
$\int_{K}u\,d\gamma=\int_{K^{\dagger}}u^{\dagger}\,d\gamma,$ then
$T_{\gamma}^{mod}(K;u)<T_{\gamma}(K^{\dagger})$, in view of the definition of
the modified torsional rigidity (7) as well as Proposition 2.4.

As a corollary, using $F(t)=t^{2},$ we immediately get

######  Corollary 4.3.

For any $u\in W_{0}^{1,2}(K),$ the Rayleigh quotient does not increase under
the $\dagger-$symmetrization. In other words,

$\frac{\int_{K}|\nabla
u|^{2}d\gamma}{\int_{K}u^{2}d\gamma}\geq\frac{\int_{K^{\dagger}}|\nabla
u^{\dagger}|^{2}d\gamma}{\int_{K^{\dagger}}(u^{\dagger})^{2}d\gamma}.$

This immediately implies

######  Corollary 4.4.

For any convex set $K$, taking $K^{\dagger}$ to be the $\dagger$-rearrangement
of $K$ with respect to the first Dirichlet eigenfunction of $L$, we have
$\Lambda_{\gamma}(K)\geq\Lambda_{\gamma}(K^{\dagger})$.

###### Proof.

Use the Corollary 4.3 with the function $u$ which minimizes the Rayleigh
quotient (i.e. with the first Dirichlet eigenfunction of $L$), to see that

$\Lambda_{\gamma}(K)=\frac{\int_{K}|\nabla
u|^{2}d\gamma}{\int_{K}u^{2}d\gamma}\geq\frac{\int_{K^{\dagger}}|\nabla
u^{\dagger}|^{2}d\gamma}{\int_{K^{\dagger}}(u^{\dagger})^{2}d\gamma}\geq\Lambda_{\gamma}(K^{\dagger}).$

∎

Finally, we are ready to outline the

Proof of Theorem 1.1. By the definition of the $\dagger-$rearrangement and the
modified torsional rigidity, $T_{\gamma}(K)\geq
T_{\gamma}(K;u)=T_{\gamma}(K^{\dagger}),$ where $u$ is taken to be the first
Dirichlet eigenfunction of $L$ on $K.$ Next, by Corollary 4.4,
$\Lambda_{\gamma}(K)\geq\Lambda_{\gamma}(K^{\dagger}).$

Therefore, we have shown that for any convex set $K$ there exists a half-space
$K^{\dagger}$ such that simultaneously,
$\Lambda_{\gamma}(K)\geq\Lambda_{\gamma}(K^{\dagger})$ and $T_{\gamma}(K)\geq
T_{\gamma}(K^{\dagger})$.

Consider $H$ to be the half-space such that $T_{\gamma}(K)=T_{\gamma}(H)$.
Then, by Lemma 2.5, $K^{\dagger}\subset H$, and thus by Lemma 2.10,
$\Lambda_{\gamma}(K^{\dagger})\geq\Lambda_{\gamma}(H)$. Hence
$\Lambda_{\gamma}(K)\geq\Lambda_{\gamma}(H)$. This finishes the proof.
$\square$

######  Remark 4.5.

Note that a more general result also follows from Theorem 4.1: given a convex
function $F:\mathbb{R}^{+}\rightarrow\mathbb{R},$ define
$\Lambda_{\gamma}^{F}(K)$ to be the minimal value of the following energy:

$\Lambda_{\gamma}^{F}(K)=\inf_{u\in W_{0}^{1,2}(K)}\frac{\int_{K}|\nabla
u|^{2}d\gamma}{\int_{K}F(u)\,d\gamma}.$

For a domain $K$, let $\tilde{K}$ be the half-space with the same Gaussian
torsional rigidity as $K$. Then
$\Lambda^{F}_{\gamma}(K)\geq\Lambda^{F}_{\gamma}(\tilde{K}).$

## References

  * [1] V. I. Bogachev, _Gaussian Measures_ , Math. Surveys Monogr. 62, Amer. Math. Soc., 1998.
  * [2] L. Brasco, _On torsional rigidity and principal frequencies: an invitation to the Kohler-Jobin rearrangement technique._ ESAIM Control Optim. Calc. Var. 20 (2014), no. 2, 315–338.
  * [3] J.E. Brothers and W.P. Ziemer, _Minimal rearrangements of Sobolev functions._ Journal fur die reine und angenwandte Mathematik 384 (1988), 153–179.
  * [4] A. Burchard, _A short course on rearrangement inequalities_ , lecture notes, 2009.
  * [5] C. Borell, _The Brunn-Minkowski inequality in Gauss spaces,_ Invent. Math. 30 (1975), 207–216.
  * [6] E.A. Carlen and C. Kerce, _On the cases of equality in Bobkov’s inequality and Gaussian rearrangement._ Calc. Var. 13 (2001), 1–18.
  * [7] A. Ehrhard, _Symetrisation dans l’espace de Gauss,_ Math. Scand. 53 (1983), 281–301.
  * [8] A. Ehrhard, _Elements extremaux pour les inegalites de Brunn-Minkowski Gaussiennes,_ Ann. Inst. H. Poincaré Probab. Statist. 22 (1986), no. 2, 149–168.
  * [9] G. Faber, _Beweis, dass unter allen homogenen Membranen von gleicher Fläche und gleicher Spannung die kreisförmige den tiefsten Grundton gibt_ Sitzungsber. Bayer. Akad. Wiss. München, Math.-Phys. Kl., 1923, 169–172.
  * [10] B. Kawohl, _Rearrangements and Convexity of Level Sets in PDE,_ Lecture Notes in Mathematics 1150, Springer, Berlin, 1985.
  * [11] S. Kesavan, _Symmetrization and Applications,_ Series in Analysis, Volume 3, World Scientific, 2006.
  * [12] M.-T. Kohler-Jobin, _Symmetrization with equal Dirichlet integrals._ SIAM J. Math. Anal. 13 (1982), no. 1, 153–161.
  * [13] M.-T. Kohler-Jobin, _Une méthode de comparaison isopérimétrique de fonctionnelles de domaines de la physique mathématique. I. Une démonstration de la conjecture isopérimétrique $P\lambda^{2}\geq\pi j^{4}_{0}/2$ de Pólya et Szegö._ (in French) Z. Angew. Math. Phys. 29 (1978), no. 5, 757–766.
  * [14] E. Krahn, _Über eine von Rayleigh formulierte Minimaleigenschaft des Kreises_ Math. Ann., 94 (1925), 97–100.
  * [15] E. Krahn, _Über Minimaleigenschaften der Kugel in drei und mehr Dimensionen_ Acta Comm. Univ. Tartu (Dorpat), A9 (1926), 1–44 (English transl.: Ü. Lumiste and J. Peetre (eds.), Edgar Krahn, 1894-1961, A Centenary Volume, IOS Press, 1994, Chap. 6, 139–174).
  * [16] E. H Lieb, M. Loss, _Analysis, Second Edition, Graduate studies in mathematics,_ vol. 14, 2001, American Mathematical Society, Providence.
  * [17] G.V. Livshyts, _On a conjectural symmetric version of Ehrhard’s inequality._ arXiv:2103.11433.
  * [18] G. Pólya, G. Szegö, _Isoperimetric Inequalities in Mathematical Physics._ Annals of Mathematics Studies, No. 27 Princeton University Press, Princeton, N. J., 1951.
  * [19] V. N. Sudakov and B. S. Tsirel’son, _Extremal properties of half-spaces for spherically invariant measures._ Problems in the theory of probability distributions, II. Zap. Nauch. Leningrad Otdel. Mat. Inst. Steklov 41 (1974), 14–24 (in Russian).
  * [20] J.L. Vazquez, _The Porous Medium Equation, Mathematical Theory_ , Oxford University Press, 2007.
  * [21] https://en.wikipedia.org/wiki/Coarea_formula?wprov=sfti1


11institutetext: Department of Biomedical Engineering, Linköping University,
Linköping, Sweden 22institutetext: Center for Medical Image Science and
Visualization, Linköping University, Linköping, Sweden 33institutetext:
Division of Statistics and Machine learning, Department of Computer and
Information Science, Linköping University, Linköping, Sweden 33email:
{marco.domenico.cirillo, david.abramian<EMAIL_ADDRESS>
# Vox2Vox: 3D-GAN for Brain Tumour Segmentation

Marco Domenico Cirillo 1122 David Abramian 1122 Anders Eklund 112233

###### Abstract

Gliomas are the most common primary brain malignancies, with different degrees
of aggressiveness, variable prognosis and various heterogeneous histological
sub-regions, i.e., peritumoral edema, necrotic core, enhancing and non-
enhancing tumour core. Although brain tumours can easily be detected using
multi-modal MRI, accurate tumor segmentation is a challenging task. Hence,
using the data provided by the BraTS Challenge 2020, we propose a 3D volume-
to-volume Generative Adversarial Network for segmentation of brain tumours.
The model, called Vox2Vox, generates realistic segmentation outputs from
multi-channel 3D MR images, segmenting the whole, core and enhancing tumor
with mean values of 87.20%, 81.14%, and 78.67% as dice scores and 6.44mm,
24.36mm, and 18.95mm for Hausdorff distance 95 percentile for the BraTS
testing set after ensembling 10 Vox2Vox models obtained with a 10-fold cross-
validation.

###### Keywords:

MRI Vox2Vox Generative Adversarial Networks deep learning artificial
intelligence 3D image segmentation.

##  1 Introduction

Gliomas are the most frequent intrinsic tumours of the central nervous system.
Based on the presence or absence of marked mitotic activity, necrosis and
ﬂorid microvascular proliferation, a malignancy grade, WHO grade II, III or IV
[33], is assigned. Gliomas with WHO grade II are also called low grade gliomas
(LGG), whereas gliomas with higher WHO grade are called high grade gliomas
(HGG). Although both these brain tumour types can easily be detected, they
have a diffuse, infiltrative way of growing in the brain, and they exhibit
peritumoural edema, such as an increase in water content in the area
surrounding the tumour. This makes it arduous to define the tumour border by
visual assessment, both in analysis and also during surgery [5].

For this reason, researchers recently started resorting to powerful
techniques, able to segment complex objects and, in this way, guide the
surgeons during the operation with a suitable accuracy. Indeed, machine
learning [27] and especially deep learning [14, 16, 19, 22, 25, 31] can
provide state-of-the-art segmentation results.

###  1.1 Related Works

Nowadays generative adversarial networks (GANs) [10] are gaining popularity in
computer vision, since they can learn to synthesise virtually any type of
image. Specifically, GANs can be used for style transfer [9], image synthesis
from noise [17], image to image translation [15], and also image segmentation
[29]. GANs have become especially popular in medical imaging [34] since
medical imaging datasets are much smaller compared to general computer vision
datasets such as ImageNet. Additionally, in medical imaging it is common to
collect several image modalities for each subject before proceeding with the
analysis, and, when this is not possible, CycleGAN introduced in [35] can be
used to synthesize the missing modalities.

GANs have also been used for medical image segmentation. Indeed, Han et al. in
[11] proposed a GAN to segment multiple spinal structures in MRIs; Li et al.
in [21] developed a novel transfer-learning framework using a GAN for robust
segmentation of different human epithelial type 2 (HEp-2) cells; Dong et al.
in [7] implemented a U-Net style GAN for accurate and timely organs-at-risk
(OARs) segmentation; Nema et al. in [26] designed a 2D GAN, called RescueNet,
to segment brain tumours from MR images; etc. Anyhow, Yi et al. [34] provide a
complete and recent review of GANs applied in medicine.

Hence, inspired by these works and especially by the Pix2Pix GAN [15], which
can generate an image of type A from a paired image of type B, the aim of this
project is to do 3D image segmentation using 3D Pix2Pix GAN, named Vox2Vox, to
segment brain gliomas. While a normal convolutional neural network, such as
U-Net [28], performs the segmentation pixel by pixel, or voxel by voxel,
through maximizing a segmentation metric or metrics (i.e. dice score,
intersection over union, etc), a GAN will also punish segmentation results
that do not look realistic. Our hypothesis is that this can result in better
segmentations.

##  2 Method

###  2.1 Data

The MR images used for this project are the Multimodal Brain tumour
Segmentation Challenge (BraTS) 2020 training ones [1, 2, 3, 4, 23]. The BraTS
2020 training dataset contains MR volumes of shape $240\times 240\times 155$
from 369 patients, and for each patient four types of MR images were
collected: native (T1), post-contrast T1-weighted (T1Gd), T2-weighted (T2),
and T2 Fluid Attenuated Inversion Recovery (FLAIR). The BraTS 2020 validation
dataset contains MR volumes from 125 patients. The images were acquired from
19 different institutions with different clinical protocols. The training set
was segmented manually, by one to four raters, following the same annotation
protocol, and their annotations were approved by experienced neuro-
radiologists, whereas no segmentation was provided for the validation set.
Moreover, all data were co-registered to the same anatomical template,
interpolated to the same resolution (1 mm3) and skull-stripped. Figure 1 shows
an example of one training T1 MR image overlapped with its true segmentation
in the theree different planes.

Figure 1: From left to right there is a T1 MR image in the sagittal, coronal
and transverse plane overlapped with its true segmentation. Peritumoural edema
(ED), necrotic and non-enhancing tumour core (NCR/NET), and GD-enhancing
tumour (ET) are highlighted in yellow, red and cyan respectively.

###  2.2 Image Pre-processing and Augmentation

For each MR image intensity normalization is done per channel, whereas the
background voxels are fixed to 0. On the other hand, the grey-scale ground-
truths are transformed into categorical, so each target has four channels, as
the number of the classes to segment: background, peritumoural edema (ED),
necrotic and non-enhancing tumour core (NCR/NET), and GD-enhancing tumour (ET)
labeled with 0, 1, 2, 4 respectively. BraTS, moreover, released 125 and 166
additional volumes as validation and test sets, respectively.

Since these volumes are memory demanding, patch augmentation is applied to
extract one sub-volume of $128\times 128\times 128$ from each original volume.
In this way, only 23.5% of the whole training set is used in every training
epoch. Moreover, in order to prevent the networks from overfitting and
memorizing the exact details of the training images, random 3D flipping,
random 3D rotations between 0°and 30°, power-law gamma intensity
transformation (gain and gamma randomly chosen between [0.8-1.2]), elastic
deformation with square deformation grid with displacements sampled from from
a normal distribution with standard deviation 5 voxels [28], or a combination
of these with probability 0.5 are applied as image augmentation techniques.

###  2.3 Model Architecture

The Vox2Vox model, as the Pix2Pix one [15], consists of a generator and a
discriminator. The generator, illustrated by Figure 2, is built with U-Net and
Res-Net [13] architecture style, see Figure 2; whereas the discriminator is
build with PatchGAN [15] architecture style, see Figure 3. Following are the
details of both model’s architecture111In the model descriptions we use I, E,
B, D, and O to refer to Input(s), Encoder, Bottleneck, Decoder, and Output
respectively.. The generator consists of:

  * I:

a 3D image with 4 channels: T1, T2, T1Gd, and T2 FLAIR;

  * E:

four down-sampling blocks, each of them made by 3D convolutions using kernel
size 4, stride 2 and same padding, followed by instance normalization [32] and
Leaky ReLU activation function. The number of filters used at the first 3D
convolution is 64 and at each down-sampling the number is doubled;

  * B:

four residual blocks, each made by 3D convolutions using kernel size 4, stride
1 and same padding, followed by instance normalization and Leaky ReLU
activation function. Every convolution-normalization-activation output is
concatenated with the previous one;

  * D:

three up-sampling blocks, each of them made by 3D transpose convolutions using
kernel size 4 and stride 2, followed by instance normalization and ReLU
activation. Each 3D convolution input is concatenated with the respective
encoder output layer;

  * O:

a 3D transpose convolution using 4 filters (as the number of the classes to
segment), kernel size 4 and stride 2, followed by softmax activation function.
The output has shape $128\times 128\times 128\times 4$ and constitutes the
segmentation prediction for each class.

On the other hand, the discriminator consists of:

  * I:

the 3D image with 4 channels and its segmentation ground-truth or the
generator’s segmentation prediction;

  * E:

same as the generator;

  * O:

one 3D convolution using 1 filter, kernel size 4, stride 1 and same padding.
The output has shape $8\times 8\times 8\times 1$ and constitutes the quality
of the segmentation prediction created by the generator.

Figure 2: The generator model. Figure 3: The discriminator model.

All 3D convolution and 3D transpose convolution layers use kernel size 4 (as
in [15]), the He et al. weight initialization method[12], and same padding.
Moreover, all Leaky ReLU layers have slope coefficient 0.3.

In order to reduce GPU memory consumption, we decided to used a single
convolution layer in the encoder and decoder blocks, while a standard U-Net
uses two convolutions layers per block. As our convolutions use stride 2, it
means that the encoder/decoder will directly downsample/upsample the input
volume a factor 2. In future work we will investigate if using two
convolutions layers per blocks improves segmentation performance.

###  2.4 Losses

Since Vox2Vox contains two models, the generator and the discriminator, two
loss functions are used. The discriminator loss, $L_{D}$, is the sum of the
$L_{2}$ error of the discriminator output, $D(\cdot,\cdot)$, between the
original image $x$ and the respective ground-truth $y$ with a tensor of ones,
and the $L_{2}$ error of the discriminator output between the original image
and the respective segmentation prediction $\hat{y}$ given by the generator
with a tensor of zeros, i.e.:

$L_{D}=L_{2}\left[D(x,y),\mathbf{1}\right]+L_{2}\left[D(x,\hat{y}),\mathbf{0}\right]\,,$
(1)

whereas, the generator loss, $L_{G}$, is the sum of the $L_{2}$ error of the
discriminator output between the original image and the respective
segmentation prediction given by the generator with a tensor of ones, and the
generalized dice loss [24, 30], $GDL(\cdot,\cdot)$, between the ground-truth
and the generator’s output multiplied by the scalar weight coefficient
$\alpha\geq 0$, i.e.:

$L_{G}=\,L_{2}\left[D(x,\hat{y}),\mathbf{1}\right]+\alpha\,GDL\left(y,\hat{y}\right)\,.$
(2)

By looking at Equation 2, it is easy to conclude that if $\alpha=0$: Vox2Vox
is a pure GAN and it minimizes only the unsupervised loss given by the
discriminator; whereas if $\alpha\rightarrow\infty$: Vox2Vox ignores the
discriminator, and behaves as a 3D U-Net with the architecture shown in Figure
2.

###  2.5 Optimization and Regularization

Both the generator and the discriminator are trained using the Adam optimizer
[18] with the parameters: $\lambda=2\cdot 10^{-4}$, $\beta_{1}=0.5$, and
$\beta_{2}=0.999$. Dropout regularization with a dropout probability of 0.2 is
used after each 3D convolutional operation in the generator’s bottleneck (see
Figure 2). Moreover, as Yi et al. reported in [34], the discriminator loss
helps the generator to guarantee the spatial consistency in the final results,
behaving as a shape regularizer. In other words, the discriminator takes care
that the generated brain segmentation looks realistic (i.e. like manual
segmentations). In the end, we expect that Vox2Vox performs better with a
trade-off $\alpha$ which does not disregard completely the discriminator loss
and, at the same time, does not disregard the generator either.

###  2.6 Model Ensembling and Post-processing

Model ensembling is a technique that combines the outputs of several models in
order to obtain more robust predictions [16]. Hence, once the Vox2Vox model is
built, $M$-fold cross-validation over the training data can be done, which
results in $M$ models. Instead of performing the ensembling independently for
each voxel, we propose a neighborhood averaging ensembling, by training a CNN
that combines the outputs of the $M$ Vox2Vox models. This CNN, called
Ensembler, is simply made by just one 3D convolution, which takes as input a
tensor with shape $(128\times 128\times 128\times 4M)$, and returns a tensor
with shape $(128\times 128\times 128\times 4)$, using stride 1, kernel size 3
and softmax activation function. Since the input of our Ensembler are softmax
outputs from each Vox2Vox model, the probabilities were zero-centered by
subtracting 0.5 prior to training. In the end, the Ensembler is trained over
100 epochs, minimizing the generalized dice loss, and early stopping with a
patience of 10 epoch (for the validation loss) is used.

Furthermore, since the $M$ models are trained to detect 4 classes, it is
normal that they sometimes detect a class which is not present for a subject
and the ensembling follows such mis-segmentation. Therefore, post-processing
can be useful in order to reduce mis-segmentation. Since the mis-segmentation
normally results in small false positives, a cluster size threshold can be
used to remove clusters smaller than a specific volume $V$.

##  3 Results

The Vox2Vox model is implemented using Python 3.7, Tensorflow 2.1 and its
Keras library. The model is trained and validated on sub-volumes of size
$128\times 128\times 128$ from 369 and 125 subjects respectively, using batch
size 4, over 200 epochs on a computer equipped with 128 GB RAM and an Nvidia
GeForce RTX 2080 Ti graphics card with 11 GB of memory. Once the training is
completed, the 166 test volumes are cropped in order to have shape $160\times
192\times 128$. In this way, the testing set can be given as input to the
fully convolutional Vox2Vox, because each axis is now divisible by $2^{4}=16$,
where 4 is the generator’s and discriminator’s depth.

Table 1 reports the dice and the Hausdorff distance 95 percentile scores for
all the classes of interest varying the $\alpha$ parameter in Equation 2 for
the training dataset. Note that the classes are reset as: whole tumour (WT =
ET $\cup$ ED $\cup$ NCR/NET), tumour core (TC = ET $\cup$ NCR/NET) and
enhancing tumour (ET). It also clearly shows that the best trade-off for the
$\alpha$ parameter is 5, but at the same time that it cannot detect the
enhancing tumor (ET) class properly over the whole training dataset. Anyway,
it also shows that the discriminator helps to achieve good results, because
the metrics obtained with $\alpha=5$ are better than when the model only
considers the generator loss (high values of $\alpha$) and also when the model
is a pure GAN ($\alpha=0$).

Anyway, it has to be pointed out that the Vox2Vox GAN is trained to detected
the classes presented in the data and not the combination of them.

Table 1: Mean dice score and Hausdorff distance 95 percentile for the different brain tumour areas over the training set. The metrics here are obtained training the model with sub-volumes of $128\times 128\times 128$ voxels, for different values of $\alpha=0,1,3,5,10,25,50,100,200,250$. | Dice score [%] | Hausdorff distance 95 [mm]
---|---|---
$\alpha$ |  WT |  TC |  ET |  WT |  TC |  ET
0 |  58.96 |  23.57 |  37.63 |  77.74 |  104.35 |  73.84
1 |  86.36 |  72.48 |  62.50 |  13.12 |  17.26 |  40.50
3 |  92.98 |  90.82 |  79.12 |  4.49 |  4.06 |  31.32
5 |  93.21 |  91.70 |  79.31 |  3.80 |  3.18 |  30.91
10 |  87.50 |  80.66 |  70.55 |  12.09 |  9.92 |  34.40
25 |  84.35 |  83.01 |  71.04 |  9.62 |  8.75 |  32.53
50 |  87.87 |  76.79 |  61.54 |  9.55 |  8.33 |  31.57
100 |  91.81 |  89.22 |  77.31 |  5.28 |  5.03 |  31.63
200 |  88.72 |  90.16 |  77.62 |  6.33 |  3.38 |  31.19
250 |  86.65 |  85.67 |  78.35 |  9.32 |  8.90 |  31.87

Moreover, it is evident that the ET class is the most problematic one to
detect. Indeed, there are just 27 subjects (7.3%) in the training dataset that
do not contain any ET voxels. So, the post-processing should focus on that
class in the end, see section 2.6.
With the alpha parameter set to 5, 10-fold cross-validation ($M=10$) over the
training set is applied, using all the image augmentation techniques listed in
section 2.2 are applied to avoid overfitting. Every epoch takes approximately
20 minutes to complete using a Keras.utils.Sequence generator. Successively,
once the 10 Vox2Vox models are trained, their outputs are combined by the
Ensembler CNN explained previously in section 2.6.

As post-processing, a threshold $th=1000$ voxels ($V=1$cm3) is set for the ET
class: if the final segmentation has a number of ET voxels fewer than $th$
ones, the ET voxels are converted into the NCR/NET class. Table 2 reports the
metrics calculated on the training, validation and test sets for each class by
the CBICA Image Processing Portal222https://ipp.cbica.upenn.edu, the name of
this group for the challenge is IMT_AE. with our proposed ensembling and post-
processing.

Table 2: Mean, standard deviation and median dice score and Hausdorff distance 95 percentile for the three different brain tumour classes over the training, validation, and test set. The predictions are calculated ensembling 10 models trained after a 10-fold cross-validation and post-processed with a threshold $th=1000$ voxels for the ET class. The values reported here were calculated by the CBICA Image Processing Portal. |  | Dice score [%] | Hausdorff distance 95 [mm]
---|---|---|---
dataset |  |  WT |  TC |  ET |  WT |  TC |  ET
|  Mean |  91.63 |  89.25 |  79.56 |  3.66 |  3.52 |  30.04
training |  StdDev |  6.36 |  11.49 |  24.74 |  3.87 |  4.57 |  96.67
|  Median |  93.39 |  92.50 |  87.16 |  2.44 |  2.23 |  1.73
|  Mean |  89.26 |  79.19 |  75.04 |  6.39 |  14.07 |  36.00
validation |  StdDev |  8.28 |  22.30 |  28.68 |  11.44 |  47.56 |  105.28
|  Median |  91.75 |  88.13 |  85.87 |  3.0 |  3.74 |  2.23
|  Mean |  87.20 |  81.14 |  78.67 |  6.44 |  24.36 |  18.95
test |  StdDev |  14.07 |  25.47 |  20.17 |  11.09 |  79.78 |  74.91
|  Median |  91.42 |  90.59 |  83.65 |  3.16 |  3.0 |  2.23

The training metrics for the WT and TC class decreased compared to those
reported in Table 1, probably due to the image augmentation techniques
introduced during the training, but just slightly; on the other hand, the ET
ones slightly increased, probably thanks to the ensembling and the post-
processing. The metrics obtained for both sets are high, but there are still
some bad predictions that compromise the mean values, which also explains the
large standard deviations. For this reason, the median values for each class
reported in Table 2 may be more representative of the typical segmentation
performance of the model.

##  4 Conclusions

Table 2 establishes that ensembling multiple Vox2Vox models generates high
quality segmentation outputs that looks realistic thanks to their
discriminators’ support, achieving median values of: 93.40%, 92.49%, and
86.48% dice scores and 2.44mm, 2.23mm, and 1.73mm Hausdorff distance 95
percentile over the training dataset; 91.75%, 88.13%, and 83.14% and 3.0mm,
3.74mm, and 2.44mm over the validation dataset; and 91.42%, 90.59%, and 83.65%
and 3.16mm, 3.0mm and 2.23mm over the testing dataset for whole tumour, core
tumour and enhancing tumour respectively.

Moreover, as future work, the Vox2Vox can be improved for the following BraTS
challenges in many ways, i.e.: training the model to optimize the BraTS
metrics and not those provided in the data, optimizing the PatchGAN
architecture, and ensembling Vox2Vox trained with different augmentation
techniques, following the suggestions reported in [6].

In the end, the Vox2Vox model can be used not only for image segmentation but
also for further image augmentation. Indeed, Vox2Vox could be combined with a
3D noise to image GAN [8, 20] which can synthesize realistic segmentation
outputs, that are then translated to realistic MR volumes. The combination of
these two GANs might result in a really fast batch generation of MR images
with their targets.

## Acknowledgement

This study was supported by LiU Cancer, VINNOVA Analytic Imaging Diagnostics
Arena (AIDA), and the ITEA3 / VINNOVA funded project Intelligence based
iMprovement of Personalized treatment And Clinical workflow supporT (IMPACT).
Funding was also provided by the Center for Industrial Information Technology
(CENIIT) at Linköping University.

## References

  * [1] Bakas, S., Akbari, H., Sotiras, A., Bilello, M., Rozycki, M., Kirby, J., Freymann, J., Farahani, K., Davatzikos, C.: Segmentation labels and radiomic features for the pre-operative scans of the TCGA-GBM collection. The Cancer Imaging Archive (2017)
  * [2] Bakas, S., Akbari, H., Sotiras, A., Bilello, M., Rozycki, M., Kirby, J., Freymann, J., Farahani, K., Davatzikos, C.: Segmentation labels and radiomic features for the pre-operative scans of the TCGA-LGG collection. The Cancer Imaging Archive (2017)
  * [3] Bakas, S., Akbari, H., Sotiras, A., Bilello, M., Rozycki, M., Kirby, J.S., Freymann, J.B., Farahani, K., Davatzikos, C.: Advancing the cancer genome atlas glioma MRI collections with expert segmentation labels and radiomic features. Scientific data 4, 170117 (2017)
  * [4] Bakas, S., Reyes, M., Jakab, A., Bauer, S., Rempfler, M., Crimi, A., Shinohara, R.T., Berger, C., Ha, S.M., Rozycki, M., et al.: Identifying the best machine learning algorithms for brain tumor segmentation, progression assessment, and overall survival prediction in the BRATS challenge. arXiv preprint arXiv:1811.02629 (2018)
  * [5] Blystad, I.: Clinical Applications of Synthetic MRI of the Brain, vol. 1600. Linköping University Electronic Press (2017)
  * [6] Cirillo, M.D., Abramian, D., Eklund, A.: What is the best data augmentation approach for brain tumor segmentation using 3d u-net? arXiv preprint arXiv:2010.13372 (2020)
  * [7] Dong, X., Lei, Y., Wang, T., Thomas, M., Tang, L., Curran, W.J., Liu, T., Yang, X.: Automatic multiorgan segmentation in thorax CT images using U-net-GAN. Medical physics 46(5), 2157–2168 (2019)
  * [8] Eklund, A.: Feeding the zombies: Synthesizing brain volumes using a 3D progressive growing GAN. arXiv preprint arXiv:1912.05357 (2019)
  * [9] Gatys, L.A., Ecker, A.S., Bethge, M.: Image style transfer using convolutional neural networks. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 2414–2423 (2016)
  * [10] Goodfellow, I., Pouget-Abadie, J., Mirza, M., Xu, B., Warde-Farley, D., Ozair, S., Courville, A., Bengio, Y.: Generative adversarial nets. In: Advances in neural information processing systems. pp. 2672–2680 (2014)
  * [11] Han, Z., Wei, B., Mercado, A., Leung, S., Li, S.: Spine-GAN: Semantic segmentation of multiple spinal structures. Medical image analysis 50, 23–35 (2018)
  * [12] He, K., Zhang, X., Ren, S., Sun, J.: Delving deep into rectifiers: Surpassing human-level performance on imagenet classification. In: Proceedings of the IEEE international conference on computer vision. pp. 1026–1034 (2015)
  * [13] He, K., Zhang, X., Ren, S., Sun, J.: Deep residual learning for image recognition. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 770–778 (2016)
  * [14] Isensee, F., Kickingereder, P., Wick, W., Bendszus, M., Maier-Hein, K.H.: No new-net. In: International MICCAI Brainlesion Workshop. pp. 234–244. Springer (2018)
  * [15] Isola, P., Zhu, J.Y., Zhou, T., Efros, A.A.: Image–to–image translation with conditional adversarial networks. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 1125–1134 (2017)
  * [16] Kamnitsas, K., Bai, W., Ferrante, E., McDonagh, S., Sinclair, M., Pawlowski, N., Rajchl, M., Lee, M., Kainz, B., Rueckert, D., et al.: Ensembles of multiple models and architectures for robust brain tumour segmentation. In: International MICCAI Brainlesion Workshop. pp. 450–462. Springer (2017)
  * [17] Karras, T., Aila, T., Laine, S., Lehtinen, J.: Progressive growing of GANs for improved quality, stability, and variation. ICLR (2018)
  * [18] Kingma, D.P., Ba, J.: Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980 (2014)
  * [19] Kong, X., Sun, G., Wu, Q., Liu, J., Lin, F.: Hybrid pyramid U-Net model for brain tumor segmentation. In: International Conference on Intelligent Information Processing. pp. 346–355. Springer (2018)
  * [20] Kwon, G., Han, C., Kim, D.s.: Generation of 3D brain MRI using auto-encoding generative adversarial networks. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 118–126 (2019)
  * [21] Li, Y., Shen, L.: cC-GAN: A robust transfer-learning framework for HEp-2 specimen image segmentation. IEEE Access 6, 14048–14058 (2018)
  * [22] McKinley, R., Meier, R., Wiest, R.: Ensembles of densely-connected CNNs with label-uncertainty for brain tumor segmentation. In: International MICCAI Brainlesion Workshop. pp. 456–465. Springer (2018)
  * [23] Menze, B.H., Jakab, A., Bauer, S., Kalpathy-Cramer, J., Farahani, K., Kirby, J., Burren, Y., Porz, N., Slotboom, J., Wiest, R., et al.: The multimodal brain tumor image segmentation benchmark (BRATS). IEEE transactions on medical imaging 34(10), 1993–2024 (2014)
  * [24] Milletari, F., Navab, N., Ahmadi, S.A.: V-net: Fully convolutional neural networks for volumetric medical image segmentation. In: 2016 Fourth International Conference on 3D Vision (3DV). pp. 565–571. IEEE (2016)
  * [25] Myronenko, A.: 3D MRI brain tumor segmentation using autoencoder regularization. In: International MICCAI Brainlesion Workshop. pp. 311–320. Springer (2018)
  * [26] Nema, S., Dudhane, A., Murala, S., Naidu, S.: RescueNet: An unpaired GAN for brain tumor segmentation. Biomedical Signal Processing and Control 55, 101641 (2020)
  * [27] Polly, F., Shil, S., Hossain, M., Ayman, A., Jang, Y.: Detection and classification of HGG and LGG brain tumor using machine learning. In: 2018 International Conference on Information Networking (ICOIN). pp. 813–817. IEEE (2018)
  * [28] Ronneberger, O., Fischer, P., Brox, T.: U-net: Convolutional networks for biomedical image segmentation. In: International Conference on Medical image computing and computer-assisted intervention. pp. 234–241. Springer (2015)
  * [29] Sato, M., Hotta, K., Imanishi, A., Matsuda, M., Terai, K.: Segmentation of Cell Membrane and Nucleus by Improving Pix2pix. In: BIOSIGNALS. pp. 216–220 (2018)
  * [30] Sudre, C.H., Li, W., Vercauteren, T., Ourselin, S., Cardoso, M.J.: Generalised dice overlap as a deep learning loss function for highly unbalanced segmentations. In: Deep learning in medical image analysis and multimodal learning for clinical decision support, pp. 240–248. Springer (2017)
  * [31] Topol, E.J.: High-performance medicine: the convergence of human and artificial intelligence. Nature medicine 25(1), 44–56 (2019)
  * [32] Ulyanov, D., Vedaldi, A., Lempitsky, V.: Instance normalization: The missing ingredient for fast stylization. arXiv preprint arXiv:1607.08022 (2016)
  * [33] Wesseling, P., Capper, D.: Who 2016 classification of gliomas. Neuropathology and applied neurobiology 44(2), 139–150 (2018)
  * [34] Yi, X., Walia, E., Babyn, P.: Generative adversarial network in medical imaging: A review. Medical image analysis p. 101552 (2019)
  * [35] Zhu, J.Y., Park, T., Isola, P., Efros, A.A.: Unpaired image-to-image translation using cycle-consistent adversarial networks. In: Proceedings of the IEEE international conference on computer vision. pp. 2223–2232 (2017)


# A tutorial on MDL hypothesis testing for graph analysis

Peter Bloem
Knowledge Representation and Reasoning Group
VU University Amsterdam
Steven de Rooij
Mathematical Institute
University of Leiden

###### Abstract

This document provides a tutorial description of the use of the MDL principle
in complex graph analysis. We give a brief summary of the preliminary
subjects, and describe the basic principle, using the example of analysing the
size of the largest clique in a graph. We also provide a discussion of how to
interpret the results of such an analysis, making note of several common
pitfalls.

##  1 Introduction

When analysing graph structure, it can be difficult to determine whether
patterns found are due to chance, or due to structural aspects of the process
that generated the data. Hypothesis tests are often used to support such
analyses. These allow us to make statistical inferences about which null
models are responsible for the data, and they can be used as a heuristic in
searching for meaningful patterns.

The minimum description length (MDL) principle [6, 4] allows us to build such
hypothesis tests, based on efficient descriptions of the data. Broadly: we
translate the regularity we are interested in into a _code_ for the data, and
if this code describes the data more efficiently than a code corresponding to
the null model, by a sufficient margin, we may reject the null model.

This is a _frequentist_ approach to MDL, based on hypothesis testing. Bayesian
approaches to MDL for model _selection_ rather than model _rejection_ are more
common, but for the purposes of pattern analysis, a hypothesis testing
approach provides a more natural fit with existing literature.

We provide a brief illustration of this principle based on the running example
of analysing the size of the largest clique in a graph. We illustrate how a
_code_ can be constructed to efficiently represent graphs with large cliques,
and how the description length of the data under this code can be compared to
the description length under a code corresponding to a _null model_ to show
that the null model is highly unlikely to have generated the data.

##  2 Preliminaries

We require some basic notions from graph theory and from information theory.
We describe these briefly, to define our notation, and to provide a concise
introduction. For a more detailed introduction, please follow the references
provided.

#### Graphs

A graph $G$ of size $n$ is a tuple $(V_{G},E_{G})$ containing a set of _nodes_
(or _vertices_) $V_{G}$ and a set of _links_ (or _edges_) $E_{G}$. For
convenience in defining probability distributions on graphs, we take $V_{G}$
to be the set of the first $n$ natural numbers. $E_{G}$ contains pairs of
elements from $V_{G}$. For the dimensions of the graph, we use the functions
$n(G)=|V_{G}|$ and $m(G)=|E_{G}|$. If a graph $G$ is _directed_ , the pairs in
$E_{G}$ are ordered, if it is _undirected_ , they are unordered. A
_multigraph_ has the same definition as a graph, but with $E_{G}$ a multiset,
i.e. the same link can occur more than once.

We will limit ourselves here to datasets that are (directed or undirected)
_simple graphs_ : i.e. no link connects a node to itself and the same link
cannot occur more than once. Probability distributions on the set of simple
graphs are usually the most complex, so that we can trust that a method that
works for simple graphs is easily translated to other types.

Two graphs $G$ and $H$ are _isomorphic_ if there exists a bijection
$f:N_{G}\to N_{H}$ on the nodes of $G$ such that two nodes $a$ and $b$ are
adjacent in $G$ if and only if $f(a)$ and $f(b)$ are adjacent in $H$. If two
graphs $G$ and $H$ are isomorphic, we say that they belong to the same
_isomorphism class_.

#### Codes and MDL

The intuition behind MDL is that compressing data and learning from it are
much the same process: in both cases, we are analyzing data to find patterns
that are characteristic for the data, or patterns that represent the data
well. This is not just an intuition, there exists a very precise
correspondence between optimizing for probability and optimizing for
description length. We will detail the basic principle below, and give an
intuitive example of how we apply the method in the next section. For more
details, we refer the reader to [4].

Let ${\mathbb{B}}$ be the set of all finite-length binary strings. We use
$|b|$ to represent the length of $b\in{\mathbb{B}}$. Let
$\log(x)=\log_{2}(x)$. A _code_ for a set of objects $\cal X$ (usually graphs,
in our case) is an injective function $f:{\cal X}\to{\mathbb{B}}$, mapping
objects to binary code words. All codes in this paper are _prefix-free_ : no
code word is the prefix of another. We will denote a _codelength function_
with the letter $L$, ie. $L(x)=|f(x)|$. It is common practice to compute $L$
directly, without explicitly computing the codewords. In fact, we will adopt
the convention of referring to $L$ itself as a code.

A well known result in information theory is the association between codes and
probability distributions, implied by the _Kraft inequality_ : for each
probability distribution $p$ on $\cal X$, there exists a prefix-free code $L$
such that for all $x\in\cal X$: $-\log p(x)\leq L(x)<-\log p(x)+1$. Inversely,
for every prefix-free code $L$ for $\cal X$, there exists a probability
distribution $p$ such that for all $x\in\cal X$: $p(x)=2^{-L(x)}$. For proofs,
see [4, Section 3.2.1] or [3, Theorem 5.2.1]. To explain the intuition, note
that we can easily transform a code $L$ into a sampling algorithm for $p$ by
feeding the decoding function random bits until it produces an output. To
transform a probability distribution to a code, techniques like arithmetic
coding [7] can be used.

As explained in [4, page 96], the fact that $-\log p^{*}(x)$ is real-valued
and $L^{*}(x)$ is integer-valued can be safely ignored and we may _identify_
codes with probability distributions, allowing codes to take non-integer
values.

When we need to encode a single choice from a finite set $S$ of options, we
can use the code with length $\log|S|$, corresponding to a uniform
distribution on $S$. In some cases, we can allow codes with multiple codewords
for a single object. If a code $L$ has multiple codewords for some object $x$,
we may indicate the choice for a particular codeword by a parameter $a$ as
$L(x;a)$.

##  3 Model Rejection by Codelength

Figure 1: A simple example of graph analysis by MDL: finding cliques. (a) The
data. (b) A representation of the data that exploits the existence of a large
clique. We first store a set indicating which nodes belong to the clique. We
then store the graph, omitting those links that belong to the clique. (c) An
illustration of the principle of a bound. For a complete code on graphs, using
$L^{\text{uniform}}$, we need to first store the size of the graph using some
arbitrary code on the integers $L^{\mathbb{N}}$. To make sure that it does not
matter which code we choose for $L^{\mathbb{N}}$, we require our code
($L^{\text{clique}}$) to beat the bound $B(G)$ (by at least $k$ bits).

In this section we provide an intuitive discussion of the principle of MDL
hypothesis testing for graph analysis. We will illustrate the idea on a simple
problem: determining whether a single large clique in a graph is an unusual
feature, or to be expected under a particular null model. For certain null
models, such as the Erdős-Renyi (ER) model, the expected size of the largest
clique is well understood [2]. For other models, less so. Without resorting to
sampling a large number of graphs from the model, how can we determine that a
particular size of clique is unusually large under a given null model? We will
define two codes, $L^{\text{null}}$ and $L^{\text{clique}}$, and show that if
$L^{\text{clique}}$ compresses the data better than $L^{\text{null}}$ by
enough bits, it is very unlikely that the data came from the null model, and
we cannot dismiss the clique as a consequence of random chance.

Consider an undirected graph $G$ containing a large clique: all nodes in some
subset $V_{C}\subseteq V_{G}$ are directly connected to one another. We can
describe the graph by first describing the set $V_{C}$, and then describing
$G$ in a canonical manner. Since every node in $V_{C}$ is connected to every
other node in the clique, we can omit these links from the second part of our
description, shortening the total description length, if $V_{C}$ is large
enough. This gives us a code $L^{\text{clique}}$ with short codelengths for
graphs containing a single large clique, and—by the correspondence mentioned
in the preliminaries—also a probability distribution $p^{\text{clique}}$ with
high probabilities for graphs with a large clique. The principle is
illustrated in Figure 1.

It is important to be aware that this is a code with _multiple codewords for
the same graph_. To convert this to a probability distribution, we sum the
probabilities of all codewords:
$p^{\text{clique}}(G)=\sum_{c}2^{-L^{\text{clique}}(G;c)}$, where each $c$
identifies a clique in $G$.111We can convert $p^{\text{clique}}$ back to a
another code by taking the negative logarithm:
$L^{\text{clique*}}(G)=-\log(G)$. $\bar{L}^{\text{clique*}}$ is a superior
code (it always lowerbounds $L^{\text{clique}}$), but inefficient to compute.
See also Section 4.1 As we shall see, we can use such codes without having to
compute this sum, _so long as they are not used as the null model_. For null
models we require a code with a single codeword per graph.

Of course, there is no guarantee that, of all the distributions with a bias
towards large cliques, $p^{\text{clique}}$ is the precise distribution from
which we sampled our data. Luckily, it does not need to be: so long as we can
show that our clique-based model encodes the data more efficiently than the
null model, the presence of the clique disproves the hypothesis that the data
came from the null model. We do not know where it came from, but it did not
come from the null model.

Under the hypothesis that the null model was the source of the data, we can
show that the probability that any other model compresses the data better by
$k$ bits or more, decays exponentially in $k$. This is known as the _no-
hypercompression inequality_ [4, p103]. Let $p^{\text{null}}(x)$ be any
probability distribution, with $L^{\text{null}}(x)=-\log p^{\text{null}}(x)$
and let $L(x)$ be _any_ code, then we have

$p^{\text{null}}\left[L^{\text{null}}(x)-L(x)\geq k\right]\leq
2^{-k}\mbox{\,.}$ (1)

To provide some intuition for why this should be so, imagine that
$L^{\text{null}}$ and $L$ are both codes for uniform distributions:
$L^{\text{null}}$ provides each graph in $X^{\text{null}}$ with a codeword of
constant length $n$, and $L$ provides each graph in $X$ with a codeword of
constant length $n-k$. Then, $L^{\text{null}}$ provides enough codewords for
$2^{n}$ graphs and, $L$ provides enough codewords for $2^{n-k}$ graphs. Thus,
if a graph can be represented by both (and thus compressed by $k$ bits), its
probability under the null model is at most $2^{-k}$. This is the basic
principle allowing us to use compression for statistical purposes: _there are
exponentially fewer short codes than long ones, so compressible data is rare_.

To give a sense of scale; under the null model, the probability that
$L^{\text{clique}}$ will compress the data better than the null model by 10
bits or more is less than one in one-thousand. For twenty bits, we get one in
a million, for thirty bits, one in a billion, and so on.

We can interpret this procedure as a hypothesis test: the difference in
compression $\Delta$ between the null model and the alternative model is a
statistic [4, Example 14.2]. The no-hypercompression inequality gives us a
bound on the probability $p^{\text{null}}(\Delta\geq k)$. To reject the null
model with significance level $\alpha$, we must find some code on the set of
all graphs and show that it compresses the data better than the null model by
$k$ bits, with $2^{-k}\leq\alpha$. Any code will do, as long as it was chosen
before seeing the data. The better we design our code, the greater the power
of the test.

Note that $\Delta$ is also the logarithm of the likelihood ratio between the
null model and $L$, so we can see this as a likelihood ratio test. We can also
interpret the difference in codelength between two models $p_{a}$ and $p_{b}$
as the logarithm of the _Bayes factor_ $p_{a}(x)/p_{b}(x)$ [4, Section
14.2.3].

Now, while our test _rejects_ the null model, it does not necessarily
_confirm_ anything about the pattern we used to compress the data. We would
like to make it as likely as possible that it was the pattern we are
interested in (e.g. the clique) that allowed us to reject the null model, and
not some other aspect of the alternative model. This is a complex problem and
it is difficult to provide any guarantees. If the alternative model stores,
say, the linkset of a graph more efficiently than the null model does, it may
well be that the alternative model beats the null models on all graphs in our
domain, regardless of the pattern used.

To rule out most such artifacts, we must re-use the same model we used for the
null model within the alternative model to store every part of the graph
except the pattern. There is usually an intuitive way to do this: for
instance, in the alternative model for the clique, we store the graph in two
parts: first the nodes belonging to the clique, and then the full graph, minus
any links in the clique. If we use, say, the Erdős-Renyi model as our null
model, we also use it to store the graph in this second part of the code. In
this way, we know that any compression achieved must be due to the
clique.222The lack of guarantees may be a cause for skepticism. Consider,
however, that a hypothesis test by itself is _never_ proof that the pattern
found is meaningful, only that the null hypothesis is incorrect. See also
Section 4.1.

###  3.1 Rejecting Multiple Null Models

A final benefit of this method is that we can reject multiple null models with
a single test. In many situations we will have a function $B(G)$ that
lowerbounds any code in some set $\cal L$. If our alternative model provides a
codelength below $B(G)-k_{\alpha}$ with $k_{\alpha}$ the number of bits
required for our chosen $\alpha$, we can reject all of $\cal L$.

As an example, Let ${\mathcal{G}}_{n}$ be the set of all undirected graphs of
size $n$. For our null model, we define a uniform code on such graphs:
$L^{\text{uniform}}_{n}=\log|{\mathcal{G}}_{n}|$. This null model captures the
idea that the size of the graph is the only informative statistic: given the
size, all graphs are equally likely. However, it is _parametrized_. It is
currently not a code on _all_ graphs, just those of size $n$. To turn it into
a code that can represent all graphs, we need to encode the parameter $n$ as
well, with some code $L^{\mathbb{N}}$ over the natural numbers

$L^{\text{complete}}(G)=L^{\mathbb{N}}(n(G))+L^{\text{uniform}}_{n(G)}(G)\mbox{\,.}$

This is called _two-part coding_ : we encode the parameters of a model first,
and then the data given the parameters. For some parametrized model
$L_{\theta}$, we can choose any code for $\theta$ to make it complete. We will
call the set of all such complete codes the _two-part codes on $L_{\theta}$_.
Note that we can simply concatenate the two codewords, since all codes are
prefix-free.

Which code we choose for the parameter is arbitrary. We may be able to reject
the uniform code for one choice of $L^{\mathbb{N}}$, or several, but how can
we prove that $L^{\text{complete}}$ will be rejected whatever $L^{\mathbb{N}}$
we choose? Instead of choosing an arbitrary code for the size, we can use the
_bound_ $B(G)=L^{\text{uniform}}_{n(G)}(G)$ as our null model. This is not a
code, but it _is_ a lower bound for any two-part code on
$L^{\text{uniform}}_{n}$. If $L^{\text{clique}}(G)$ is shorter than $B(G)$, it
is also shorter than $L^{\text{complete}}(G)$ whatever the choice of
$L^{\mathbb{N}}$.333 In probabilistic terms, the code on the parameter
corresponds to a prior on the parameter. The two-part codes correspond to
maximum likelihood posterior probabilities:
$p(\hat{\theta})p(x\mid\hat{\theta})$. Our bound corresponds to the maximal
likelihood of the data: $p(x\mid\hat{\theta})$. This shows us that the bound
applies not only to the two-part codes, but also to the full Bayesian mixture:
$\sum_{\theta}p(x\mid\theta)p(\theta)\leq\sum_{\theta}p(x\mid\hat{\theta})p(\theta)=p(x\mid\hat{\theta})$

Note that when we store the rest of the graph within $L^{\text{clique}}$ we
_cannot_ use $B(G)$ in place of $L^{\text{complete}}(G)$. We want a
_conservative_ hypothesis test: the probability of incorrectly rejecting a
null model may be lower than $\alpha$ but never higher. By this principle,
bounds chosen in place of either model should always decrease $\Delta$. The
code corresponding to the null model must always be lowerbounded, and the
optimal code for the alternative model must always be upperbounded.

This is also the reason that we allow the alternative code to have multiple
codewords for one graph. Replacing such a code with the complete equivalent
model, found by summing over all codewords, would increase $\Delta$,
strengthening the hypothesis test. By evaluating only one codeword, we are
effectively using an upper bound on this complete model, trading off power for
efficiency of computation.

##  4 Discussion

We will conclude by discussing several subtleties n the interpretation of this
method.

###  4.1 Multiple Testing and Significance as Pattern Evidence

In some cases, we may test multiple patterns for the same graph. For instance,
we may test several cliques found in the same data, to strengthen our attempt
to reject the null hypothesis. At first glance this may seem like multiple
testing, raising the question of whether we should adjust our rejection
region. Does the probability that one or more of the cliques we find to be
significantly large is a “false positive” increase with the number of
subgraphs checked?

This question highlights an important issue in the use of hypothesis testing
in pattern analysis: the question about the pattern we are trying to answer—is
$C$ a meaningful clique?—is not the question we are testing for. Strictly
interpreted, the hypothesis test investigates the question _whether the full
data came from the null model_. Thus, while we test multiple cliques, we only
test them on a _single_ null-hypothesis. To perform a valid single null-
hypothesis test, we can simply select the best compressing clique out of all
the cliques tested, and use the log-factor of that clique as the significance
with which we reject the hypothesis that the null model generated the data.

This is a valid code: we can check as many patterns as we like before choosing
the one we will use to compress the data with. Since it is a valid code, it is
a valid, _single_ hypothesis test, and there is no need to correct for
multiple testing.444If we test multiple _null models_ on a single dataset,
multiple testing _should_ be taken in to account.

To make this more precise, assume that we have some data $G$. The null model
assigns codelength $L^{\text{null}}(G)$ and the clique code assigns
$L^{\text{clique}}(G;C)$ for a particular clique $C$. The clique code gives us
many different ways to describe this graph (one code word for each clique). To
remove this inefficiency, we can define a code that sums the probability mass
of all these codewords, so that it has only a single codeword per graph:

$\displaystyle p^{\text{clique}}(G;C)$
$\displaystyle=\text{exp}_{2}\left[{-L^{\text{clique}}(G;C)}\right]$
$\displaystyle p^{\text{clique*}}(G)$
$\displaystyle=\sum_{C}p^{\text{clique}}(G;C)$ $\displaystyle
L^{\text{clique*}}(G)$ $\displaystyle=-\log p^{\text{clique*}}(G)$

The codeword that $L^{\text{clique*}}$ provides for $G$ is smaller than any
codeword provided by $L^{\text{clique}}$.

We can now use $L^{\text{clique*}}$ to perform a _single_ test: if
$L^{\text{clique*}}$ compresses better than $L^{\text{null}}$, then by the no-
hypercompression inequality (1), we can reject the null model. We cannot
compute $L^{\text{clique*}}$, for any reasonable size of graph, but we _can_
upperbound it. If we can find only _one_ codeword that allows
$L^{\text{clique}}$ to compress better than $L^{\text{null}}$ by the required
margin, we have proved that $L^{\text{clique*}}$ does so as well.555As usual,
if we do not find such a codeword, we do not know whether or not the null
hypothesis can be rejected. This allows us to test as many cliques as we like,
and we may use exhaustive search or sampling to find these: the real test is
whether $L^{\text{clique*}}$ compresses the data better than the null model.

Of course, in general practice, we are not actually interested in whether or
not the null model is true. What we are investigating is whether a given
_pattern_ is “true”. We are using the framework of hypothesis testing merely
as a _heuristic_ to guide unsupervised pattern mining. Does the general
intuition of multiple hypothesis testing still apply? Depending on the type of
pattern under investigation, we may not be able to exclude artifacts similar
to false positives: every subgraph of the largest clique may prove a
“positive” pattern, even though it is not the the pattern we are ultimately
looking for.

However, increasing the number of patterns checked does not raise the amount
of such false positives. This is because we analyze multiple patternss on a
_single dataset_. The process generating the data is stochastic, but once the
data is known, the log-factors of all patterns are determined. Checking more
patterns will get us closer to the true number of positives, but it will not
put us in danger of an increased number of _false_ positives.

###  4.2 Limitations and Extensions

Whether the approach described here is useful for a particular type of
pattern, depends on its properties. One illustrative example is an _anti-
motif_ : a subgraph that occurs less than expected under the null model. In
principle, any deviation from the null-model can be used to compress the data.
However, computing such a code efficiently for anti-motifs may be a
challenging task. For other patterns (hub nodes, motifs, $k$-regularity),
suitable codes are easier to devise.

Note also that the constraints on the null model are stronger than those on
the alternative model. In our code, we technically use the length of one of
many codewords, rather than the true probability (which would require a sum
over all codewords). In the null model, such approximations are _not_
appropriate, and each graph must have a single codeword. For some models (such
as the configuration model [5]), this makes the codelength so expensive to
compute that we may gain little over traditional methods of hypothesis
testing.

Many null models can easily be used in this framework. For instance, for the
$G(n,m)$ and $G(n,p)$ models can the codelength can easily be computed. For
other models, like the preferential attachment model [1], we can easily come
up with a code with multiple codewords for each graph, but a code with a
single codeword per graph is more complex to compute efficiently.

#### Acknowledgements

This publication was supported by the Dutch national program COMMIT, by the
Netherlands eScience center, and by the Amsterdam Academic Alliance Data
Science (AAA-DS) Program Award to the UvA and VU Universities.

## References

  * [1] Réka Albert and Albert-László Barabási.  Statistical mechanics of complex networks.  Reviews of modern physics, 74(1):47, 2002.
  * [2] Béla Bollobás and Paul Erdős.  Cliques in random graphs.  In Mathematical Proceedings of the Cambridge Philosophical Society, volume 80, pages 419–427. Cambridge University Press, 1976.
  * [3] Thomas M. Cover and Joy A. Thomas.  Elements of information theory (2. ed.).  Wiley, 2006.
  * [4] P.D. Grünwald.  The minimum description length principle.  The MIT Press, 2007.
  * [5] Mark Newman.  Networks: an introduction.  Oxford University Press, 2010.
  * [6] Jorma Rissanen.  Modeling by shortest data description.  Automatica, 14(5):465–471, 1978.
  * [7] Jorma Rissanen and Glen G Langdon.  Arithmetic coding.  IBM Journal of research and development, 23(2):149–162, 1979.