# Leading the Pack: N-player Opponent Shaping

Alexandra Souly1 Timon Willi1 Akbir Khan2 Robert Kirk1
Chris Lu2 Edward Grefenstette1 Tim Rocktäschel1
1University College London 2University of Oxford
<EMAIL_ADDRESS>
###### Abstract

Reinforcement learning solutions have great success in the 2-player general
sum setting. In this setting, the paradigm of Opponent Shaping (OS), in which
agents account for the learning of their co-players, has led to agents which
are able to avoid collectively bad outcomes, whilst also maximizing their
reward. These methods have currently been limited to 2-player game. However,
the real world involves interactions with many more agents, with interactions
on both local and global scales. In this paper, we extend Opponent Shaping
(OS) methods to environments involving multiple co-players and multiple
shaping agents. We evaluate on over 4 different environments, varying the
number of players from 3 to 5, and demonstrate that model-based OS methods
converge to equilibrium with better global welfare than naive learning.
However, we find that when playing with a large number of co-players, OS
methods’ relative performance reduces, suggesting that in the limit OS methods
may not perform well. Finally, we explore scenarios where more than one OS
method is present, noticing that within games requiring a majority of
cooperating agents, OS methods converge to outcomes with poor global welfare.

##  1 Introduction

Multi agent systems are increasingly deployed to the real world, in which
interactions may be either cooperative, competitive or both. These can be
considered general-sum games where interactions combine opportunities for
coordination with self-interested motivations, where each player is trying to
maximise their own reward. Real-world examples of such games include global
pollution, deforestation, over-fishing and arms races. In these settings,
agents are incentivised to exploit the system and/or other players for their
own gains by defecting, but cooperation would result in greater social
welfare.

General sum games have multiple equilibria, many of which do not produce high
social welfare. Naive learning algorithms that work well in zero-sum settings
do not perform well in general-sum settings: they learn the Nash-equilibrium
strategies by best-responding to the opponent’s past behaviour, but in these
settings many Nash-equilibria often coincides with worst-case outcomes for all
participants [Letcher et al., 2019].

Opponent shaping (OS) methods are a family of approaches that take into
account the evolving social dynamics instead of simply best-responding to the
opponent’s behaviour [Foerster et al., 2018a, Lu et al., 2022, Zhao et al.,
2022, Willi et al., 2022, Khan et al., 2023]. They act to influence other
participant’s learning dynamics to result in a favorable equilibrium. These
methods have been successfully used in 2-player general sum settings to
achieve high rewards, however they do not always achieve pro-social outcomes,
and can sometimes successfully exploit their opponent.

The study of OS methods has so far been limited to only two interacting
agents, whilst the real world has many more agents. The interactions of more
players give rise to more complicated and unstable social dynamics, making
cooperation harder to achieve. In particular, in diverse populations
equilibrium selection requires modelling not only a single co-player’s beliefs
but also their beliefs of other co-players.

Figure 1: Normalised global welfare with a single shaping agent. We report the
converged (normalised) global welfare for training a shaping agent with unseen
co-players across games and player variants. Across most games shaper
outperforms LOLA and as number of players increase, both OS methods fail to
find better welfare maximising equilibria than the naive baseline. Results are
reported over 5 seeds with mean and standard deviation.

In this work, we introduce a natural extension for OS methods to more-
than-2-player setting. We then evaluate this extension over a total of 12
varying games. We find that OS methods consistently help achieve higher global
welfare than baseline methods. While OS-methods are only maximising their own
reward, they do so while effectively avoiding collective welfare minimising
outcomes, creating better social dynamics and thus also improving their co-
players’ rewards. We also find that between model-based and model-free OS
methods, model-free methods perform better in games with more players. In some
instances where OS methods are unable to exploit their opponents, they are
also willing to take lower rewards in order to avoid a collectively bad
outcome. This is promising for scenarios that involve some more capable self-
interested agents in a population of less capable agents.

##  2 Background

Model-based OS methods such as LOLA Foerster et al. [2018a], COLA Willi et al.
[2022], POLA Zhao et al. [2022], SOS Letcher et al. [2019] and Meta-MAPG Kim
et al. [2021] assume access to the other agents parameters, and as such,
anticipate their co-player’s updates as they learn. For example, LOLA
conditions its update on a co-player’s next set of parameters (effectively
thinking one step ahead). COLA, POLA and SOS are further adaptations and
extensions of LOLA, addressing different shortcomings of the original method.
Meta-MAPG also differentiates through the opponent’s current parameters via an
additional peer learning objective.

In comparison, model-free OS methods [Lu et al., 2022, Khan et al., 2023, Fung
et al., 2023, Lu et al., 2023] frame opponent-shaping as a meta-reinforcement
learning problem, where the goal is to maximise reward over a series of games
against a co-player that is learning. In these approaches a shaping agent
trains over $T$ trials consisting of $E$ episodes. Each trial consists of the
following: the shaping agent plays the game with the co-player for $E$
episodes, with the co-player updating after every episode according to its
update rule. After every episode, the shaping agent updates its hidden state
to add the played games into its memory, whilst the opponent takes a full
gradient learning step. After a trial the shaping agent updates its policy
parameters. Importantly, the shaping agent does not update during an episode
unlike the co-player.

##  3 Methods

We expand shaper Khan et al. [2023] to a multi-player setting and generalise
the algorithm to arbitrary numbers of shaping-agents and co-players. Co-
players can be any type of other agents, including shaper, other meta-agents
or naive learners.

We extend opponent-shaping to a generic number of shaping agents and co-
players naturally. Similar to single shaper single co-player shaping, shaping
agents play the game with their co-players for $E$ episodes but with all co-
players simultaneously updating after an episode. After every episode, all
shaping agents simultaneously update their respective hidden states. Shaping
agents learn to collectively shape their naive co-players.

We update shapers using evolutionary strategies [Salimans et al., 2017], which
are computationally better suited to the long timescale of trials and are
efficiently parallelisable. Given a game with $N$ shapers and $K$ co-players,
at the start of a trial we sample a population of groups of shapers. These
groups then each play against copies of the same co-players. At the end of the
trial, we sample a new generation of groups based on the parameters of the
most successful groups. Importantly, there is no mixture between the groups.
shapers do not play with members of other groups, only with the shapers in
their own group, therefore the whole group evolves together. We present a
multiplayer extension in Appendix B.

##  4 Experiments

For all games, one episode consists of 100 iterated cooperate-defect choices.
One trial consists of 1000 episodes. We evaluate on the following general sum
games:

Iterated Prisoner’s Dilemma (IPD) is a game in which agents must commit to
cooperating (staying silent) or defecting (implicating your comrades) [Yao and
Darwen, 2001]. For each player, assuming that the actions of others players is
fixed, defecting always results in a better payoff than cooperating (that is
to say, defection is the dominant choice). All players choosing cooperation is
global welfare maximising, but unstable with naive learners.

Snowdrift is a game in which players choose to shovel (cooperate) or ignore
the snow (defect) [Wettergren, 2021]. The shovelling players share a fixed
cost. If a player doesn’t shovel, they don’t incur a cost but still receive
the benefits of the snow is shovelled. If no one shovels, all players receive
a low payoff. This game differs from IPD in that mutual defection is not the
best choice - if no one else cooperates, it is still better for the player to
shovel alone than if everyone defects and the road stays blocked.

Tragedy of the Commons (ToC) is an N-player social dilemma. We use a general
formulation in which there is a social benefit B such that if enough players
pay a cost C, everyone gets the benefit. The social benefit can be thought of
as the shared resource, and the cost as the opportunity loss of not exploiting
the resource. This game encourages free-riding, as long as enough of the co-
players cooperate to prevent the tragedy.

Stag Hunt N players choose either to hunt a stag (cooperate) or to hunt a hare
(defect) [Pacheco et al., 2008]. A certain number of hunters are required to
successfully hunt a stag (in our case more than half the total population),
the hunt is unsuccessful otherwise. If there is a successful hunt, the payoff
is related to the number of cooperating hunters (more hunters implies bigger
game). The stag hunters pay the cost of the hunt but everyone equally shares
the benefits if the hunt is successful. This setup also encourages free-riding
similarly to ToC.

Figure 2: Normalised global welfare with multiple shaping agents. In IPD, the
optimal amount of shapers is around a third of the population, too many
decreases global welfare. In comparison in Snowdrift - a game where only a
single agent is required to cooperate, more shaping agents help coordinate
which naive learner is exploited. In ToC and Stag Hunt, where cooperation
requires a majority of players, increasing the number of shapers did not
increase welfare.  Figure 3: Normalised shaper versus naive co-player welfare.
We report the (normalised) converged scores attained by shapers and their
naive co-players. In all games, shapers either score similarly to their naive
co-players, or heavily exploit them, particularly in IPD and Snowdrift.

Methods \- We consider three different methods for shaping: naive learners, a
model-based OS method (LOLA) and the model-free OS method shaper. Naive
learners (NLs) are agents that don’t account for the learning or adaptation of
other agents in the environment. Updates are myopic, at the end of single
episodes, and make no modelling assumptions of other players. Our naive
learners are parameterized by a recurrent neural network and keep their hidden
states between episodes. Their networks are updated using proximal policy
optimisation (PPO) Schulman et al. [2017] and trained until convergence. For
our model-free OS method we use LOLA-PG - a policy gradient approximation to
estimate the co-player’s gradient update. We use the Infinitely Differentiable
Monte Carlo Estimator [Foerster et al., 2018b, DICE] to approximate the
gradients, and train until convergence. For our model free OS method we use
shaper, a meta-agent parameterised by a LSTM network, specifically a GRU with
hidden state of size 16. Models are trained until convergence.

##  5 Results, Discussion & Conclusion

Single Shaping Agent We firstly conduct experiments in which we evaluate a
single shaper with $N-1$ naive learners across the many environments. In line
with prior work [Lu et al., 2022] we find that shaper outperforms LOLA across
most tasks in improving global welfare (see Figure 1). We also find that
current OS methods provide diminishing results when applied on games with more
than 2 players. In all 3-player settings, shaper exploits its co-players, but
has varying success with more co-players (Figure 3). In particular, shaper has
difficulties on Staghunt and ToC, in scenarios where a majority of players is
required for cooperation to be rewarded.

We note that in ToC the difficulty of cooperation is not monotonic in number
of players. In the 4-player game, 3 are required to coordinate to get welfare
maximising outcome, whilst in the 5-player variant the required number is also
3. Hence, with randomly searching agents it is easier in the 5-player variant
to discover welfare maximising outcomes than the 4-player one. (50% vs 31%).

Many Shaping Agents With many shaping agents, we explore the stability of OS
methods when interacting with one another. In games where cooperation is
related to only a small percentage of a population cooperating, more shapers
improves global welfare (see Figure 2). In Snowdrift, only a single agent must
cooperate, therefore increasing the number of shapers increases performance as
exploiting the naive learners becomes easier. Conversely, in IPD, where
incentives for cooperation are lower, we find multiple shapers reduce global
welfare when the majority of the population become shapers, due to overly
exploiting the remaining naive learners (Figure 3). In games such as ToC and
Staghunt, cooperation is only beneficial if enough players coordinate. Two
shapers successfully coordinate in the 3 player game (where only one player
can defect), with one shaper cooperating and the other defecting. However they
do not succeed at coordination when more players or shapers are present,
suggesting a limitation to the paradigm.

In conclusion, while OS methods improve global welfare over naive learning in
multi-agent games, they might do so by exploitation of their co-players. When
the game dynamics do not allow for this, they are willing to take lower
rewards by cooperating. This is promising for improving the dynamics of
finite-resource games like ToC by including more capable agents in the
population. The method, however, faces limitations as the number of agents and
shapers increases, especially in games requiring majority cooperation. Future
work could extend this to the STORM environments [Willi et al., 2023,
Rutherford et al., 2023] with partial observability and temporally-extended
action spaces.

## References

  * Foerster et al. [2018a] Jakob N. Foerster, Richard Y. Chen, Maruan Al-Shedivat, Shimon Whiteson, Pieter Abbeel, and Igor Mordatch.  Learning with opponent-learning awareness, 2018a.
  * Foerster et al. [2018b] Jakob N. Foerster, Gregory Farquhar, Maruan Al-Shedivat, Tim Rocktäschel, Eric P. Xing, and Shimon Whiteson.  Dice: The infinitely differentiable monte carlo estimator.  In Jennifer G. Dy and Andreas Krause, editors, _Proc. of ICML_ , volume 80 of _Proceedings of Machine Learning Research_ , pages 1524–1533. PMLR, 2018b.  URL http://proceedings.mlr.press/v80/foerster18a.html.
  * Fung et al. [2023] Kitty Fung, Qizhen Zhang, Chris Lu, Timon Willi, and Jakob Nicolaus Foerster.  Analyzing the sample complexity of model-free opponent shaping.  In _ICML Workshop on New Frontiers in Learning, Control, and Dynamical Systems_ , 2023.  URL https://openreview.net/forum?id=Dm2fbPpU6v.
  * Khan et al. [2023] Akbir Khan, Timon Willi, Newton Kwan, Andrea Tacchetti, Chris Lu, Edward Grefenstette, Tim Rocktäschel, and Jakob Foerster.  Scaling opponent shaping to high dimensional games.  _arXiv preprint arXiv:2312.12568_ , 2023.
  * Kim et al. [2021] Dong-Ki Kim, Miao Liu, Matthew Riemer, Chuangchuang Sun, Marwa Abdulhai, Golnaz Habibi, Sebastian Lopez-Cot, Gerald Tesauro, and Jonathan P. How.  A policy gradient algorithm for learning to learn in multiagent reinforcement learning.  In Marina Meila and Tong Zhang, editors, _Proc. of ICML_ , volume 139 of _Proceedings of Machine Learning Research_ , pages 5541–5550. PMLR, 2021.  URL http://proceedings.mlr.press/v139/kim21g.html.
  * Kuhn [2019] Steven Kuhn.  Prisoner’s Dilemma.  In Edward N. Zalta, editor, _The Stanford Encyclopedia of Philosophy_. Metaphysics Research Lab, Stanford University, Winter 2019 edition, 2019.
  * Letcher et al. [2019] Alistair Letcher, Jakob N. Foerster, David Balduzzi, Tim Rocktäschel, and Shimon Whiteson.  Stable opponent shaping in differentiable games.  In _Proc. of ICLR_. OpenReview.net, 2019.  URL https://openreview.net/forum?id=SyGjjsC5tQ.
  * Lu et al. [2023] Chris Lu, Timon Willi, Alistair Letcher, and Jakob Nicolaus Foerster.  Adversarial cheap talk.  In _International Conference on Machine Learning_ , pages 22917–22941. PMLR, 2023.
  * Lu et al. [2022] Christopher Lu, Timon Willi, Christian A. Schröder de Witt, and Jakob N. Foerster.  Model-free opponent shaping.  In Kamalika Chaudhuri, Stefanie Jegelka, Le Song, Csaba Szepesvári, Gang Niu, and Sivan Sabato, editors, _International Conference on Machine Learning, ICML 2022, 17-23 July 2022, Baltimore, Maryland, USA_ , volume 162 of _Proceedings of Machine Learning Research_ , pages 14398–14411. PMLR, 2022.  URL https://proceedings.mlr.press/v162/lu22d.html.
  * Pacheco et al. [2008] Jorge Pacheco, Francisco Santos, and Max Souza.  Evolutionary dynamics of collective action in n-person stag hunt dilemmas.  _Proceedings. Biological sciences / The Royal Society_ , 276:315–21, 2008.  doi: 10.1098/rspb.2008.1126.
  * Rutherford et al. [2023] Alexander Rutherford, Benjamin Ellis, Matteo Gallici, Jonathan Cook, Andrei Lupu, Gardar Ingvarsson, Timon Willi, Akbir Khan, Christian Schroeder de Witt, Alexandra Souly, Saptarashmi Bandyopadhyay, Mikayel Samvelyan, Minqi Jiang, Robert Tjarko Lange, Shimon Whiteson, Bruno Lacerda, Nick Hawes, Tim Rocktaschel, Chris Lu, and Jakob Nicolaus Foerster.  Jaxmarl: Multi-agent rl environments in jax, 2023.
  * Salimans et al. [2017] Tim Salimans, Jonathan Ho, Xi Chen, Szymon Sidor, and Ilya Sutskever.  Evolution strategies as a scalable alternative to reinforcement learning, 2017.
  * Schulman et al. [2017] John Schulman, Filip Wolski, Prafulla Dhariwal, Alec Radford, and Oleg Klimov.  Proximal policy optimization algorithms, 2017.
  * Wettergren [2021] Thomas A. Wettergren.  Replicator dynamics of an n-player snowdrift game with delayed payoffs.  _Applied Mathematics and Computation_ , 404:126204, 2021.  ISSN 0096-3003.  doi: https://doi.org/10.1016/j.amc.2021.126204.  URL https://www.sciencedirect.com/science/article/pii/S0096300321002940.
  * Willi et al. [2022] Timon Willi, Alistair Letcher, Johannes Treutlein, and Jakob N. Foerster.  COLA: consistent learning with opponent-learning awareness.  In Kamalika Chaudhuri, Stefanie Jegelka, Le Song, Csaba Szepesvári, Gang Niu, and Sivan Sabato, editors, _International Conference on Machine Learning, ICML 2022, 17-23 July 2022, Baltimore, Maryland, USA_ , volume 162 of _Proceedings of Machine Learning Research_ , pages 23804–23831. PMLR, 2022.  URL https://proceedings.mlr.press/v162/willi22a.html.
  * Willi et al. [2023] Timon Willi, Akbir Khan, Newton Kwan, Mikayel Samvelyan, Chris Lu, and Jakob Foerster.  Pax: Multi-agent learning in jax.  https://github.com/ucl-dark/pax, 2023.
  * Yao and Darwen [2001] Xin Yao and Paul Darwen.  An experimental study of n-person iterated prisoner’s dilemma games.  _Informatica_ , 18, 2001.  doi: 10.1007/3-540-60154-6_50.
  * Zhao et al. [2022] Stephen Zhao, Chris Lu, Roger B Grosse, and Jakob Foerster.  Proximal learning with opponent-learning awareness.  _Advances in Neural Information Processing Systems_ , 35:26324–26336, 2022.

##  Appendix A Additional Environment Details

For each game, we present payoff matrices for their work.

IPD Global welfare increases linearly with the number of cooperators, but
agents are incentivised to defect for individually higher payoffs. Therefore
this game presents a free-riding problem. The payoff matrix we use is
structured as in Table 1: defection results in 1 more point than cooperation
given the actions of the co-players, and for each co-player cooperating our
reward increases by 2.

No. of C among remaining n-1 players
---
|  | 0 | $1$ | $2$ | … | $n-1$
Player A | C | $0$ | $2$ | $4$ | … | $2(n-1)$
D | $1$ | $3$ | $5$ | … | $2n-1$
Table 1: N-player IPD payoff where n is the total number of players, C is
cooperation action and D is defection action.

StagHunt

The payoffs are calculated with the following equations, where N is the
population size, and the reward is more than the individual cost of the hunt:

ifnum coop $\displaystyle\geq\lceil\frac{N}{2}\rceil:$ C payoff
$\displaystyle=\frac{\text{num coop}*\text{reward}}{N}-\text{hunt cost}$ D
payoff $\displaystyle=\frac{\text{num coop}*\text{reward}}{N}$ otherwise: C
payoff $\displaystyle=-\text{hunt cost}$ D payoff $\displaystyle=0$

In our experiments we set the cost of the hunt to 3 and the reward for a
successful hunt to 6.

Tragedy of the Commons

The payoff matrix we use is presented in Table 2, as per Kuhn [2019]. In our
experiments, we set the benefit to 5 and the cost to 3.

| more than T choose C | T or fewer choose C
---|---|---
C | benefit - cost | -cost
D | benefit | 0
Table 2: Tragedy of the Commons payoff where T is the threshold for resource
exhaustion, C is cooperation action and D is defection action

##  Appendix B Multi-Shaping Method

Following algorithm:

Algorithm 1 Training multiple shaper agents $\\{S\\}$ against co-players
$\\{-S\\}$

.

0: Stochastic game $\mathcal{M}$

1: Initialise all shaper parameters $\phi_{i}$ for $i\in\\{S\\}$

2: for $t=0$ to $T$ do

3: Initialise trial reward $\bar{J}=0$

4: Initialise co-players $\phi_{i}$ for $i\in\\{-S\\}$

5: Initialise all agent hidden state $h_{i}=\textbf{0}$ for
$i\in\\{S\\}\bigcup\\{-S\\}$

6: for $e=0$ to $E$ do

7: Roll out trajectories $J_{i},h^{\prime}_{i}$ =
$\mathcal{M}(\phi_{\\{S\\}},\phi_{-\\{S\\}},h_{\\{S\\}},h_{-\\{S\\}})$ with
$i\in\\{S\\}\bigcup\\{-S\\}$

8: Update co-player’s parameters $\phi_{-i}$ with $i\in\\{-S\\}$

9: Update all hidden states $h_{i}\leftarrow{}h^{\prime}_{i}$ with
$i\in\\{S\\}\bigcup\\{-S\\}$

10: Update trial reward $\bar{J}\leftarrow\bar{J}+J_{i}$

11: end for

12: Update shaper parameters $\phi_{i}$ with respect to $\bar{J}$ with
$i\in\\{S\\}$

13: end for

##  Appendix C Hyperparameters

Hyperparameter | Value
---|---
Number of Actor Hidden Layers | 1
Number of Critic Hidden Layers | 1
Torso GRU Size | [16]
Length of Trial | 1000
Length of Episode | 100
Number of Generations | until convergence
Number of parallel environments | 2
Number of parallel opponents | 10
Population Size | 100
OpenES sigma init | 0.04
OpenES sigma decay | 0.999
OpenES sigma limit | 0.01
OpenES init min | 0.0
OpenES init max | 0.0
OpenES clip min | -1e10
OpenES clip max | 1e10
OpenES lrate init | 0.01
OpenES lrate decay | 0.9999
OpenES lrate limit | 0.001
OpenES beta 1 | 0.99
OpenES beta 2 | 0.999
OpenES eps | 1e-8
Table 3: Hyperparameters for SHAPER

Hyperparameter | Value
---|---
Number of Minibatches | 10
Number of Epochs | 4
Gamma | 0.96
GAE Lambda | 0.95
PPO clipping epsilon | 0.2
Value Coefficient | 0.5
Clip Value | True
Max Gradient Norm | 0.5
Anneal Entropy | False
Entropy Coefficient Start | 0.1
Entropy Coefficient Horizon | 400000000
Entropy Coefficient End | 0.01
LR Scheduling | False
Learning Rate | 0.0003
ADAM Epsilon | 1e-5
With CNN | False
Table 4: Hyperparameters for PPO


# Behavior of Kreiss bounded $C_{0}$-semigroups on a Hilbert space

L. Arnold  Jana i Jedrzeja Śniadeckich 8, 00-656 Warszawa, Poland
<EMAIL_ADDRESS>
###### Abstract.

Let $0<\alpha\leq 1$. We prove that a $\alpha$-Kreiss bounded $C_{0}$
semigroup $(T_{t})_{t\geq 0}$ on a Hilbert space has asymptotics
$\left\|T_{t}\right\|=\mathcal{O}\big{(}t^{\alpha}/\sqrt{log(t)}\big{)}$.
Then, we give an application to perturbed wave equation.

###### Key words and phrases:

$C_{0}$-semigroup, Kreiss condition, perturbed wave equation

###### 2020 Mathematics Subject Classification:

Primary 47D06 ; Secondary 34D05, 35B45

Let $-A$ be the generator of a $C_{0}$-semigroup $(T_{t})_{t\geq 0}$ on a
complex Banach space $X$. For $\alpha>0$, we say that $(T_{t})_{t\geq 0}$ is
$\alpha$-Kreiss bounded (or simply Kreiss bounded when $\alpha=1$) if
$\sigma(A)\subset\overline{\mathbb{C}_{+}}=\\{\lambda\in\mathbb{C}:\,\operatorname{Re}(\lambda)\geq
0\\}$ and

(1)
$\left\|R(\lambda,A)\right\|\leq\frac{C}{(-\operatorname{Re}(\lambda))^{\alpha}},\quad\forall\operatorname{Re}(\lambda)<0.$

According to [5, (2.6)] (see also [8, Theorem 1.1.(1)]) , if $(T_{t})_{t\geq
0}$ is an $\alpha$-Kreiss bounded $C_{0}$-semigroup on a Hilbert space for
$\alpha>0$, then $\left\|T_{t}\right\|=O(t^{\alpha})$. Let us fix
$\gamma\in(0,1)$. In [3], the authors give an example (see [3, Example 4.4.])
of Kreiss bounded $C_{0}$-semigroup for which there exists a constant $C>0$
such that $\left\|T_{t}\right\|\geq Ct^{\gamma}$ for $t\geq 1$. A natural
question (listed as an open question in [7, Appendix B]) is whether the
estimate $\left\|T_{t}\right\|=O(t^{\alpha})$, for $\alpha$-Kreiss bounded
$C_{0}$-semigroup $(T_{t})_{t\geq 0}$ on a Hilbert space, is sharp. In the
following we answer this question for $\alpha\in(0,1]$.

######  Theorem 1.

Let $0<\alpha\leq 1$ and let $-A$ be the generator of $\alpha$-Kreiss bounded
$C_{0}$-semigroup $(T_{t})_{t\geq 0}$ on a Hilbert space. Then

(2)
$\left\|T_{t}\right\|=\underset{t\rightarrow\infty}{\mathcal{O}}\Big{(}t^{\alpha}/\sqrt{log(t)}\Big{)}.$

The proof of the Theorem 1 is based on techniques given in [2]. We first state
and prove two lemmas. It is known (see [4] or [9]) that a $C_{0}$-semigroup
$(T_{t})_{t\geq 0}$ on a Hilbert space, with negative generator $A$, is
bounded if and only if

$\left\|R(-r+i\cdot,A)x\right\|^{2}_{L^{2}(\mathbb{R},H)}\leq\frac{C}{r}\left\|x\right\|^{2},\quad
r>0$

and

$\left\|R(-r+i\cdot,A^{*})x^{*}\right\|^{2}_{L^{2}(\mathbb{R},H)}\leq\frac{C}{r}\left\|x^{*}\right\|^{2},\quad
r>0.$

The following lemma gives an analogous necessary condition for a
$C_{0}$-semigroup to be $\alpha$-Kreiss bounded, for $\alpha>0$.

######  Lemma 2.

Let $\alpha>0$ and $(T_{t})_{t\geq 0}$ be a $C_{0}$-semigroup on a Hilbert
space whose negative generator $A$. Assume $(T_{t})_{t\geq 0}$ is
$\alpha$-Kreiss bounded then

(3)
$\left\|R(-r+i\cdot,A)x\right\|^{2}_{L^{2}(\mathbb{R},H)}\leq\frac{C(1+r^{\alpha})^{2}}{r^{2\alpha}}\left\|x\right\|^{2},\quad
r>0$

and

(4)
$\left\|R(-r+i\cdot,A^{*})x^{*}\right\|^{2}_{L^{2}(\mathbb{R},H)}\leq\frac{C(1+r^{\alpha})^{2}}{r^{2\alpha}}\left\|x^{*}\right\|,\quad
r>0.$

###### Proof.

Let $x\in H$. First, we by [6, formula (7.1) section 1.7], we have

$R(-r+i\beta,A)x=-\int_{0}^{\infty}e^{i\beta s}e^{-rs}T_{s}xds,\quad
r>0,\beta\in\mathbb{R}.$

This means that for $r>0$,

(5) $R(-r+i\beta,A)x=-\mathcal{F}^{-1}(s\mapsto e^{-rs}T_{s}x)(\beta),$

where $\mathcal{F}:L^{2}(\mathbb{R},H)\mapsto L^{2}(\mathbb{R},H)$ is the
Fourier-Plancherel operator. By Fourier-Plancherel Theorem

(6) $\left\|\mathcal{F}^{-1}(s\mapsto
e^{-rs}T_{s}x)(\cdot)\right\|_{L^{2}(\mathbb{R},H)}=\frac{1}{\sqrt{2\pi}}\left\|e^{-r\cdot}T_{(\cdot)}x\right\|_{L^{2}(\mathbb{R}_{+},H)}.$

The Gearhart-Prüss Theorem (see [1, Theorem 5.2.1.]) says that for each $r>0$,
$\underset{t>0}{\sup}\left\|e^{-rt}T_{t}\right\|<\infty$. This implies that
there exists $K>0$ such that

$\underset{a>1}{\sup}\left\|e^{-a\cdot}T_{(\cdot)}x\right\|_{L^{2}(\mathbb{R}_{+},H)}\leq
K\left\|x\right\|$

and then

(7)
$\underset{a>1}{\sup}\left\|R(-a+i\beta,A)x\right\|_{L^{2}(\mathbb{R}_{+},H)}\leq\frac{K}{\sqrt{2\pi}}\left\|x\right\|.$

Now, let $r>0$ and $a=r+1>1$. By the resolvent identity, we have, for
$\beta\in\mathbb{R}$,

$\displaystyle R(-r+i\beta,A)x$
$\displaystyle=\big{(}I+(r-a)R(-r+i\beta,A)\big{)}R(-a+i\beta,A)x$
$\displaystyle=\big{(}I-R(-r+i\beta,A)\big{)}R(-a+i\beta,A)x.$

Using (1) and (7):

$\displaystyle\left\|R(-r+i\cdot,A)x\right\|^{2}_{L^{2}(\mathbb{R},H)}$
$\displaystyle\leq(1+\sup_{\beta\in\mathbb{R}}\left\|R(-r+i\beta,A)\right\|)^{2}\underset{a^{\prime}<-1}{\sup}\left\|R(-a^{\prime}+i\cdot,A)x\right\|^{2}_{L^{2}(\mathbb{R}_{+},H)}$
$\displaystyle\leq 2\pi
K^{2}\Big{(}1+\frac{C}{r^{\alpha}}\Big{)}^{2}\left\|x\right\|^{2}$
$\displaystyle\leq
C^{\prime}\frac{(1+r^{\alpha})^{2}}{r^{2\alpha}}\left\|x\right\|^{2}.$

In the same way we obtain (4).

∎

The next lemma says that if $(T_{t})_{t\geq 0}$ and $(T^{*}_{t})_{t\geq 0}$
satisfy a kind of Cesàro condition then $\left\|T_{t}\right\|$ satisfy (2).

######  Lemma 3.

Let $\alpha\in(0,1)$ and $(T_{t})_{t\geq 0}$ be a $C_{0}$-semigroup on a
Banach space $X$. Assume that there exists $C>0$ such that for each $t>1$,

$\int_{0}^{t}\left\|T_{s}x\right\|^{2}ds\leq
Ct^{2\alpha}\left\|x\right\|^{2},\quad x\in X,$

and

$\int_{0}^{t}\left\|T^{*}_{s}x^{*}\right\|^{2}ds\leq
Ct^{2\alpha}\left\|x^{*}\right\|^{2},\quad x^{*}\in X^{*}.$

Then, $\left\|T_{t}\right\|=\mathcal{O}(t^{\alpha}/\sqrt{log(t)})$.

###### Proof.

Let $t>2$ and $1<P<Q<t$. We have

$\displaystyle(Q-P)^{2}|\langle T_{t}x,x^{*}\rangle|^{2}$
$\displaystyle=\Big{|}\int_{P}^{Q}\langle T_{t-s}x,T_{s}^{*}x^{*}\rangle
ds\Big{|}^{2}$
$\displaystyle\leq\Big{(}\int_{P}^{Q}\left\|T_{t-s}x\right\|^{2}ds\Big{)}\Big{(}\int_{P}^{Q}\left\|T_{s}^{*}x^{*}\right\|^{2}ds\Big{)}$
$\displaystyle\leq\Big{(}\int_{P}^{Q}\left\|T_{t-s}x\right\|^{2}ds\Big{)}\Big{(}\int_{0}^{Q}\left\|T_{s}^{*}x^{*}\right\|^{2}ds\Big{)}$
$\displaystyle\leq
CQ^{2\alpha}\left\|x^{*}\right\|^{2}\Big{(}\int_{t-Q}^{t-P}\left\|T_{s}x\right\|^{2}ds\Big{)}$

Finally taking the supremum over $\\{\left\|x^{*}\right\|=1\\}$ one obtains,
since $\alpha\in(0,1]$,

(8) $\displaystyle\frac{(Q-P)^{2}}{Q^{2}}\left\|T_{t}x\right\|^{2}$
$\displaystyle\leq
CQ^{2(\alpha-1)}\int_{t-Q}^{t-P}\left\|T_{s}x\right\|^{2}ds$
$\displaystyle\leq C\int_{t-Q}^{t-P}\left\|T_{s}x\right\|^{2}ds$

Now we set $L:=\left\lfloor log(t)/log(2)\right\rfloor$. For $0\leq l\leq
L-1$, $Q=2^{l+1}$ and $P=2^{l}$ we have

$\left\|T_{t}x\right\|^{2}\leq
4C\int_{t-2^{l+1}}^{t-2^{l}}\left\|T_{s}x\right\|^{2}ds$

Therefore

$L\left\|T_{t}x\right\|^{2}=\sum_{l=0}^{L-1}\left\|T_{t}x\right\|^{2}\leq\sum_{l=0}^{L-1}4C\int_{t-2^{l+1}}^{t-2^{l}}\left\|T_{s}x\right\|^{p}ds\leq
4C\int_{0}^{t}\left\|T_{s}x\right\|^{2}ds\leq 4C^{2}t^{2\alpha}$

Hence for $t>2$,

$\left\|T_{t}x\right\|\leq\frac{2Ct^{\alpha}}{\sqrt{L}}\left\|x\right\|\leq\frac{C^{\prime}t^{\alpha}}{\sqrt{log(t)}}\left\|x\right\|.$

∎

We are now ready to prove Theorem 1.

###### Proof of Theorem 1.

Let $x\in H$ and $r>0$. By (5), (6) and Lemma 2,

$\left\|e^{-r\cdot}T_{(\cdot)}x\right\|_{L^{2}(\mathbb{R}_{+},H)}^{2}=\left\|R(-r+i\cdot,A)x\right\|^{2}_{L^{2}(\mathbb{R},H)}\leq\frac{C(1+r^{\alpha})^{2}}{r^{2\alpha}}\left\|x\right\|^{2}.$

Furthermore

$\left\|e^{-r\cdot}T_{(\cdot)}x\right\|^{2}_{L^{2}(\mathbb{R}_{+},H)}=\int_{\mathbb{R}}e^{-2rs}\left\|T_{s}x\right\|^{2}ds\geq
e^{-2rt}\int_{0}^{t}\left\|T_{s}x\right\|^{2}ds.$

Taking $t=\frac{1}{r}$ we obtain

$\int_{0}^{t}\left\|T_{s}x\right\|^{2}ds\leq
Ce^{2}t^{2\alpha}\big{(}1+\frac{1}{t^{\alpha}}\big{)}^{2}\left\|x\right\|^{2}.$

Hence, for $t>1$,

$\int_{0}^{t}\left\|T_{s}x\right\|^{2}ds\leq
C^{\prime}t^{2\alpha}\left\|x\right\|^{2},$

where $C^{\prime}=4e^{2}C$. Similarly, since $(T_{t}^{*})_{t\geq 0}$ is Kreiss
bounded $C_{0}$-semigroup on a Hilbert space, for $t>1$ and $x^{*}\in H$, we
have

$\int_{0}^{t}\left\|T^{*}_{s}x^{*}\right\|^{2}ds\leq
C^{\prime}t^{2\alpha}\left\|x^{*}\right\|^{2}.$

Finally, Lemma 3 allows us to conclude that
$\|T_{t}\|=\mathcal{O}(t^{\alpha}/\sqrt{log(t)})$.

∎

######  Remark 4.

  1. (1)

When $\alpha>1$, (8) with $Q=2$ and $P=1$ gives

$\left\|T_{t}x\right\|^{2}\leq
C2^{\alpha}\int_{t-2}^{t-1}\left\|T_{s}x\right\|^{2}ds\leq
C^{\prime}\int_{0}^{t}\left\|T_{s}x\right\|^{2}ds\leq
C^{\prime\prime}t^{2\alpha}\left\|x\right\|^{2}.$

Then $\left\|T_{t}\right\|=O(t^{\alpha})$, which gives another proof of [8,
Theorem 1.1.(1) with $g(s)=s^{\alpha}$ for $\alpha>1$].

  2. (2)

Let us notice that if $(T_{t})_{t\geq 0}$ satisfies the $\alpha$-Kreiss
bounded condition only for $\operatorname{Re}(\lambda)\in(0,1)$ that is

$\left\|R(\lambda,A)\right\|\leq\frac{C}{-\operatorname{Re}(\lambda)^{\alpha}},\quad\forall\operatorname{Re}(\lambda)\in(-1,0),$

then the conclusion
$\left\|T_{t}\right\|=\mathcal{O}(t^{\alpha}/\sqrt{log(t)})$ remains true.

Let give an application of the Theorem 1. Let $W^{2,s}(\mathbb{T}^{2})$ be the
second order Sobolev space equipped with the standard norm. Let $A$ be the
operator with domain $D(A)=W^{2,2}(\mathbb{T}^{2})\times
W^{2,1}(\mathbb{T}^{2})$ defined by

$A=\begin{pmatrix}0&-1\\\ -\Delta-M\frac{\partial}{\partial
x}&0\end{pmatrix},$

where $\Delta$ is the Laplacian with $D(\Delta)=W^{2,2}(\mathbb{T}^{2})$ and
$M:L^{2}(\mathbb{T}^{2})\rightarrow L^{2}(\mathbb{T}^{2})$ is the
multiplication operator given by $M(h)(x,y)=e^{iy}h(x,y)$. Then, it is known
(see [8, section 4] for details) that $-A$ generates a $C_{0}$-group
$(T_{t})_{t\in\mathbb{R}}$ on the Hilbert space
$H=W^{1,2}(\mathbb{T}^{2})\times L^{2}(\mathbb{T}^{2})$, satisfying
$\left\|T_{t}\right\|=\underset{t\rightarrow\pm\infty}{\mathcal{O}}(|t|e^{|t|/2})$.
With Theorem 1, we are able to improve this estimate:

######  Proposition 5.

Let $(T_{t})_{t\in\mathbb{R}}$ be the $C_{0}$-group on $H$ generated by $-A$.
Then

(9)
$\left\|T_{t}\right\|=\underset{t\rightarrow\pm\infty}{\mathcal{O}}\Bigg{(}\frac{|t|}{\sqrt{log(|t|)}}e^{|t|/2}\Bigg{)}.$

###### Proof.

According to [8, Lemma 4.4.], for $\lambda\in\mathbb{C}_{+}$, with
$\operatorname{Re}(\lambda)\in(0,1)$,

$\left\|R\Big{(}\lambda,A+\frac{1}{2}\Big{)}\right\|\leq\frac{C}{-\operatorname{Re}(\lambda)}\quad{and}\quad\left\|R\Big{(}\lambda,-A+\frac{1}{2}\Big{)}\right\|\leq\frac{C}{-\operatorname{Re}(\lambda)}.$

Then, according to Remark 4 (2), the semigroup $(e^{-1/2t}T_{t})_{t\geq 0}$
(resp. $(e^{-1/2t}T_{-t})_{t\geq 0}$ ) which is generated by
$-(A+\frac{1}{2})$ (resp. $A+\frac{1}{2}$), satisfies
$\left\|e^{-1/2t}T_{t}\right\|=\mathcal{O}(t/\sqrt{log(t)})$ (resp.
$\left\|e^{-1/2t}T_{-t}\right\|=\mathcal{O}(t/\sqrt{log(t)})$). This yields
the estimate (9).

∎

Acknowledgements. Author was supported by the ERC grant Rigidity of groups and
higher index theory under the European Union’s Horizon 2020 research and
innovation program (grant agreement no. 677120-INDEX).

## References

  * [1] A. Arendt, C. Batty, M. Hierber, and F. Neubrander.  Vector-valued Laplace Transforms and Cauchy Problems : Second Edition.  Monograph in Mathematics 96, Springer Basel AG, 2011.
  * [2] G Cohen, C Cuny, T Eisner, and M Lin.  Resolvent conditions and growth of powers of operators.  J. Math. Anal. Appl., 487, 2020.
  * [3] T Eisner and H Zwart.  Continuous-time kreiss resolvent condition on infinite-dimensional spaces.  Math. Comput., 75:1971–1985, 2006.
  * [4] A.M. Gomilko.  Conditions on the generator of a uniformly bounded $\textsc{C}_{0}$-semigroup.  Funct. anal. Appl. 33, pages 294–296, 1999.
  * [5] B Helffer and J Sjoestrand.  From resolvent bounds to semigroup bounds.  https://arxiv.org/pdf/1001.4171v1.pdf, 2010.
  * [6] A. Pazy.  Semigroups of Linear Operators and Applications to PDEs.  Springer, 1992.
  * [7] J Rozendaal.  Operator-valued $(l^{p},l^{q})$ fourier multipliers and stability theory for evolution equations.  https://arxiv.org/pdf/2204.02223.pdf, 2022.
  * [8] J Rozendaal and M Veraar.  Sharp growth rates for semigroups using resolvent bounds.  Journal of Evolution Equations, 18:1721–1744, 2018.
  * [9] D.-H. Shi and Feng D.X.  Characteristic conditions of the generation of $\textsc{C}_{0}$-semigroups in a hilbert space.  J. Math. Anal. Appl. 247 : 356-376, 2000.


Further author information: Christiaan G.A. Viviers: E-mail:
<EMAIL_ADDRESS>Telephone: +31 (0)6 206 60171

# Probabilistic 3D Segmentation for Aleatoric Uncertainty Quantification in
full 3D Medical Data

Christiaan G.A. Viviers  M.M. Amaan Valiuddin  Eindhoven University of
Technology, 5612 AP, Eindhoven, the Netherlands;  Peter H.N. de With
Eindhoven University of Technology, 5612 AP, Eindhoven, the Netherlands;  Fons
van der Sommen  Eindhoven University of Technology, 5612 AP, Eindhoven, the
Netherlands;

###### Abstract

Uncertainty quantification in medical images has become an essential addition
to segmentation models for practical application in the real world. Although
there are valuable developments in accurate uncertainty quantification methods
using 2D images and slices of 3D volumes, in clinical practice, the complete
3D volumes (such as CT and MRI scans) are used to evaluate and plan the
medical procedure. As a result, the existing 2D methods miss the rich 3D
spatial information when resolving the uncertainty. A popular approach for
quantifying the ambiguity in the data is to learn a distribution over the
possible hypotheses. In recent work, this ambiguity has been modeled to be
strictly Gaussian. Normalizing Flows (NFs) are capable of modelling more
complex distributions and thus, better fit the embedding space of the data. To
this end, we have developed a 3D probabilistic segmentation framework
augmented with NFs, to enable capturing the distributions of various
complexity. To test the proposed approach, we evaluate the model on the LIDC-
IDRI dataset for lung nodule segmentation and quantify the aleatoric
uncertainty introduced by the multi-annotator setting and inherent ambiguity
in the CT data. Following this approach, we are the first to present a 3D
Squared Generalized Energy Distance ($D^{2}_{\text{GED}}$) of 0.401 and a high
0.468 Hungarian-matched 3D IoU. The obtained results reveal the value in
capturing the 3D uncertainty, using a flexible posterior distribution
augmented with a Normalizing Flow. Finally, we present the aleatoric
uncertainty in a visual manner with the aim to provide clinicians with
additional insight into data ambiguity and facilitating more informed
decision-making. Our code is publicly available at:
https://github.com/cviviers/prob_3D_segmentation

###### keywords:

Aleatoric Uncertainty, Segmentation, Normalizing Flows, Volumetric Medical
Data

##  1 Introduction

With the broad acceptance of deep learning-based computer-aided diagnosis and
computer-aided detection (CAD) methods, increasingly more requirements are
being posed for their successful deployment. These CAD methods typically
assist clinicians with decision-making about potentially critical medical
procedures or the planning thereof. While most research has focused on
maximizing a specific accuracy metric, in practice, a high accuracy along with
a strong indication of potential uncertainty or ambiguity in the model output
is extremely valuable and highly relevant. Deep learning-based semantic
segmentation methods using convolutional neural networks have successfully
been adopted as CAD methods for a wide range of medical imaging modalities.
While research has been conducted towards quantifying the types of uncertainty
occurring when using a segmentation model, most of this work is limited to the
quantification of the uncertainty in two-dimensional slices or images, where
the latter often originate from a 3D volume such as in CT and MRI. However,
these approaches fail to exploit the rich 3D features that may help in
resolving ambiguities in the volume. In this research, we use the full 3D
volume as input to derive a reliable estimate of the aleatoric uncertainty in
the segmentation output. This approach potentially enables a better 3D
visualization of this uncertainty in clinical practice.

Two types of uncertainty are typically prevalent in deep learning-based image
analysis methods: aleatoric and epistemic uncertainty [1]. Aleatoric
uncertainty is an estimate of the intrinsic, irreducible ambiguity in data. It
is usually associated with inherent noise in the data and its acquisition
process. Epistemic uncertainty is the uncertainty about the model, either as a
result of the architecture or the true parameter values of the model due to
limited knowledge, e.g. a finite training set size. In this work, we employ
the LIDC-IDRI lung CT dataset [2], which makes use of multiple ground-truth
annotations per lung nodule. As described in earlier work of Valiuddin et al.
[3], the epistemic uncertainty – i.e. preferences, experiences and knowledge –
of the annotators manifests into aleatoric uncertainty when providing
annotations as ground-truth data. The different annotations per nodule adds
ambiguity during training of a segmentation network. During the annotation
process, the radiologist typically annotates on a single 2D plane of the 3D
volume. However, whilst annotating, full access to the other two views of the
CT scan are typically available on the same screen. This allows the annotator
to correct the annotation if it does not align with the other two views and as
a result, the annotator creates a true 3D annotation. In the LIDC-IDRI
dataset, as described in Section 3.3, annotators were allowed multiple rounds
of annotation, thereby potentially increasing the quality of the annotation by
using the full 3D information available.

In recent research, various methods have been proposed to quantify the
uncertainty arising in segmentation models or resulting from images. One
increasingly popular approach is the Probabilistic U-Net [4], as proposed by
Kohl et al.. They propose combining a 2D U-Net with a conditional variational
autoencoder (VAE) capable of learning a distribution over the possible
annotations and ultimately construct a generative segmentation model. The
Probabilistic U-Net provides compelling results in resolving the ambiguity in
an image. More recently, various improvements to this model have been proposed
[5, 3, 6], paving the way as one of the leading uncertainty quantification
methods. Selvan et al. [5] and Valiuddin et al. [3] propose adding a
Normalizing Flow to the posterior network of the Probabilistic U-Net. This
allows the model to move away from modelling the ambiguity as strictly axis-
aligned Gaussian and, instead, allows for a learned posterior distribution of
varying complexity.

This work provides the following contributions. First, we propose a 3D
probabilistic framework, which builds upon the research from both Kohl et al.
and Valiuddin et al. that exploits the full-3D spatial information to resolve
the uncertainty in the original CT volumes. Second, it is shown that more
diverse segmentations are obtained when the posterior distribution is enhanced
by a Normalizing Flow. This finding suggests that such modeling enables
capturing the uncertainty more accurately. Third, we test the proposed
method’s ability to capture uncertainty on the LIDC-IDRI lung nodule datatset
and are the first to present results in the 3D version of the
$D^{2}_{\text{GED}}$ metric, and show that a high segmentation accuracy is
obtained using a Hungarian-matched 3D IoU.

In this work, we present a 3D probabilistic segmentation model that exploits
the full 3D-spatial information to more accurately quantify the aleatoric
uncertainty, while maintaining 3D consistency. The proposed model is equipped
with a Normalizing Flow to eliminate the strictly Gaussian latent space that
is currently enforced to further improve the model’s ability to quantify the
uncertainty. More specifically, the model consists of a 3D U-Net and a 3D
conditional VAE enhanced with the Normalizing Flow, to generate a diverse set
of plausible segmentations.

##  2 Related Work

The Probabilistic U-Net [4] is capable of generating a diverse set of valid
segmentation hypotheses. This model was proposed by Kohl et al. as a method
for capturing the ambiguity in an image. Shi Hu et al. [6] showed how this
ambiguity can be interpreted as uncertainty. They propose adding variational
dropout and an additional inter-grader variability term to the training
objective of the Probabilistic U-Net. This change allows the model to capture
a combination of the epistemic and aleatoric uncertainty. Selvan et al. [5]
improve the quantification of aleatoric uncertainty by adding a Normalizing
Flow to the posterior network of the VAE [7] in the Probabilistic U-Net.
Valiuddin et al. [3] indicate that an additional metric, the Hungarian-matched
IoU, is necessary to effectively evaluate the performance of these uncertainty
quantification methods.

The aforementioned studies have all been conducted with a focus on 2D images
and slices from the original 3D volume, which is a step away from the domain
where the true uncertainty resides. In addition, these methods all heavily
focus on a specially crafted subset of the LIDC-IDRI dataset. The image
patches used to train and test the models exactly contain the regions where
the ambiguity resides, and nothing more. In practice, at test time, this is
not the case. In this work we use a fixed 3D volume surrounding the lung
nodules and consider all the image patches, with or without uncertainty.

##  3 Methods

###  3.1 Model Architecture

In this research, we extend the Probabilistic U-Net to the 3D domain and
address a key limitation by augmenting the posterior network with a
Normalizing Flow (NF). The network consists of a 3D U-Net, a Prior network,
Posterior network enhanced with an NF, and a Feature Combination network. By
combining the 3D spatial features extracted by the U-Net with samples taken
from a latent distribution encapsulating the solution space, a set of diverse,
but plausible segmentations can be generated. The standard deviation across
these predictions can be interpreted as the aleatoric uncertainty. The U-Net
[8] and 3D U-Net [9] have shown time and again their ability to segment
structures of interest at state-of-the-art performance. While we employ the 3D
U-Net to obtain the relevant 3D spatial information, this approach is generic
and allows for any other segmentation network to be used. A deep CNN
conditioned on the input CT scan is used to model a low-dimensional axis-
aligned Gaussian latent space, representing the segmentation variants (Prior
distribution). Another CNN-based axis-aligned Gaussian encoder (Posterior
network) that is conditioned on both the query CT scan and a ground-truth
segmentation, is utilized during training to model a posterior low-dimensional
latent space. Valiuddin et al. point out the shortcomings in modelling the
posterior distribution to be strictly Gaussian. As such, we augment our
posterior network with either a 2-step planar or radial flow [10], to
potentially increase the complexity of the captured posterior distribution,
and thereby provide more meaningful updates to the prior network during
training. Figure 1 portrays a detailed diagram of the proposed network
architecture.

Figure 1: Diagram of the 3D Probabilistic U-Net with an augmented flow
posterior. The bottom network depicts the 3D U-Net, the Prior and Posterior
network are shown at the top and a Feature combination network at the right
combines samples taken from the captured distributions. The diagram
additionally depicts both the training and testing configuration.

During training, the Probabilistic 3D U-Net makes use of the Posterior
network, Prior network, U-Net and the feature combination layers. Samples are
taken from the image-label conditional distribution captured by the Posterior
network and combined with the features extracted from the U-Net through the
Feature Combination network. The loss is then computed using Equation (1). The
Prior network follows the Posterior network during training, as enforced by
the KL-divergence, and thus learns to capture this image-label conditional
distribution from the image alone. At test time, the Posterior network is
discarded and samples are taken from the Prior network instead. It should be
noted that there is only one forward pass through the U-Net (for image feature
extraction) and the Prior network (to capture the image-conditional
distribution). However, multiple passes through the feature combination
network are made, in order to combine a new sample from the Prior distribution
with the image features.

###  3.2 Loss Function & Evaluation Criteria

In line with previous work on conditional variational autoencoders, our
training objective consists of minimizing the variational lower bound [11].
This entails minimizing a cross-entropy difference (in our case) between the
ground-truth segmentation (y) and a prediction (s), minimizing the Kullback-
Leibler ($\operatorname{KL}$) divergence between the posterior distribution
($p_{\phi}$) and the prior distribution ($p_{\theta}$), and finally, a
correction term for the density transformation through the Normalizing Flow
[12]. Given a query image ($\mathbf{x}$) and a posterior sample
($\mathbf{z}$), $p_{\psi}$ (a U-Net and Feature combination network) generates
a plausible segmentation (s). This loss term can formally be specified by

$\displaystyle\mathcal{L}(\mathbf{y},\mathbf{x},\theta,\phi,\psi)$
$\displaystyle=-\,\mathbb{E}_{p_{\phi}(\mathbf{z}|\mathbf{y},\mathbf{x})}[\operatorname{log}p_{\psi}(\mathbf{y}|\mathbf{z},\mathbf{x})]+\beta\cdot\left(\operatorname{KL}\left(\,p_{\phi}(\mathbf{z}_{0}|\mathbf{y},\mathbf{x})||p_{\theta}(\mathbf{z}|\mathbf{x})\,\right)-\mathbb{E}_{p_{\phi}(\mathbf{z}_{0}|\mathbf{y},\mathbf{x})}\left[\sum_{i=1}^{K}\log\left(\left|\operatorname{det}\frac{df_{i}}{d\mathbf{z}_{i-1}}\right|\right)\right]\right).$
(1)

These losses are combined and weighed using hyperparameter $\beta$ [13, 14].
The original research provides a detailed derivation of the elbo loss [11],
how it is used in the context of Probabilistic U-Net [4] and the NF-likelihood
objective [10, 3, 12].

The metric $D^{2}_{\text{GED}}$ or Squared Generalized Energy Distance has
become the de-facto metric in the context of uncertainty quantification and
the quantification of the distance between distributions of segmentations.
This metric is defined as

$\displaystyle
D^{2}_{\text{GED}}(P_{\text{GT}},P_{\text{Out}})=2\mathbb{E}\left[d(\mathbb{S},\mathbb{Y})\right]-\mathbb{E}\left[d(\mathbb{S},\mathbb{S}^{\prime})\right]-\mathbb{E}\left[d(\mathbb{Y},\mathbb{Y}^{\prime})\right],$
(2)

where $d$ is a distance measure, $1-\text{IoU}(x,y)$, in our case. The
parameters $\mathbb{S}$ and $\mathbb{S}^{\prime}$ are independent samples from
the predicted distribution $P_{\text{Out}}$. The parameters $\mathbb{Y}$ and
$\mathbb{Y}^{\prime}$ are the 4 samples from the ground-truth distribution
$P_{\text{GT}}$. In addition to the $D^{2}_{\text{GED}}$, we also report the
Hungarian-matched IoU. This compensates for a shortcoming in the
$D^{2}_{\text{GED}}$ that when the predictions are relatively poor, the metric
rewards sample diversity by definition. We duplicate the ground-truth set (4
annotations) to match the desired sample number when computing the Hungarian-
matched IoU. This measure calculates the distance between two discrete
distributions by determining an optimal coupling between the ground-truth and
prediction set subject to the IoU metric.

###  3.3 Dataset & Data Preparation

To evaluate the proposed method’s ability to resolve the ambiguity in the
data, we use the popular LIDC-IDRI dataset. This dataset contains the lung CT
scans from 1,010 patients with manual lesion annotations from up to 4 experts.
In total, there are 1,018 CT scans potentially containing multiple lung
nodules of different levels of malignancy. In this work, we have used the
annotations from a second reading, in which the radiologists were presented an
anonymized version of the annotations from other experts and were allowed to
make adjustments to their own annotations. Contrary to previous work, we use
every nodule in the dataset if it has been annotated by at least one
radiologist (potentially missed by three), regardless of the shape or severity
of the nodule. We pre-process the CT scans by clustering all nodule
annotations for a scan, by computing a distance measure between the
annotations. If an annotation is within one voxel spacing of that particular
CT scan from another annotation, it is grouped to belong to the same nodule.
The scan is resampled to a 0.5 mm along the $x$ and $y$-dimensions and 1 mm
along the $z$-dimension. This is followed by cropping the CT scan and
resulting annotations based on the center of mass of the first annotator’s
mask with a dimension of 96$\times$180$\times$180 voxels in the
$z,x,y$-dimensions. Finally, if the nodule does not have at least 4
annotations, the ground-truth (GT) masks are filled with empty annotations.
This addition is made to be consistent with previous work on this dataset [4,
3] and to capture the difficulty in detecting a nodule. This results in a
total of 2,651 3D patches, each containing a nodule and four annotations. An
example of the nodule in the CT scan and the four ground-truth annotations are
depicted in Figure 2.

###  3.4 Experiments

To compare the proposed approach against prior work, we conduct six
experiments. We train the (1) original Probabilistic 2D U-Net and the (2)
Radial NF-augmented Probabilistic 2D U-Net on 2D axial slices of the 3D
volume. In practice, we filter the slices based on the presence of at least
one positive annotation from any of the raters and use them for training, to
avoid a heavily imbalanced training set. The (3) 3D U-Net, (4) Probabilistic
3D U-Net and an (5-6) NF-augmented (Radial and Planar) Probabilistic 3D U-Net
are then trained on the 3D patches. In contrast to prior work where the 3D
lesion was sliced and split into 2D images, where some 2D slices potentially
land in the training set and some in the validation/test set, we conduct our
experiments on a per-lesion basis. This avoids any potential model bias caused
by the splitting and makes the proposed approach more clinically relevant,
since we can present the uncertainty for each lesion.

Figure 2: Example nodule in a slice from the CT scan and the four ground-truth
annotations.

In our implementation, we split the nodule data in a 70/15/15 training &
validation/test split. During training, we randomly sample one of the four
annotations to be used as ground-truth segmentation and crop the CT volume and
label to 64$\times$128$\times$128 voxels. In line with previous work, for the
3D Probabilistic U-Net, the dimensionality of the latent space is set to
$L=6$. The proposed framework is implemented in PyTorch and extends on the
work conducted by Wolny et al. [15]. We train using a batch size of 32 in the
2D case and 4 in the 3D case. An Adam optimizer with an initial learning rate
of $1\times 10^{-4}$ and a weight decay of $1\times 10^{-5}$ is used. We
reduce the learning rate by a factor of 0.2 if the validation loss does not
decrease after 20 epochs. The $\beta$ parameter is controlled using a cosine
cyclical annealing strategy as descibed by Fu et al. [13]. In all our 3D
Probabilistic U-Net experiments, we use the same hyperparameters and a
hardware configuration with an RTX 3090Ti GPU (available from Nvidia Inc.
Santa Clara, CA, USA). Training to completion takes about 2 days on the
average. For performance evaluations, we report results using a more readily
available RTX 2080Ti GPU.

##  4 Results

The results of our experiments are shown in Figure 3, Figure 4, Figure 5 and
Table 1. In Figure 3, example predictions from all the models used in our
experiments are showcased for qualitative evaluation. Here, $\mu$ GT refers to
the mean segmentation of the four raters and $\mu$ Pred is the mean of the
predictions. This mean prediction is the segmentation recommended by the
Probabilistic 3D U-Net. Additionally, the figure depicts the variation in the
segmentations. More specifically, the standard deviation of the ground-truth
labels ($\sigma$ GT) and the logits (after sigmoid activation) resulting from
the model predictions ($\sigma$ Pred) are depicted. In the figure it can be
observed that this deviation across the predictions can be interpreted as the
uncertainty. We scale the uncertainty heat map visualization to the the
maximum standard deviation of a particular model’s predictions. Additionally,
the figure depicts a rather conservative segmentation of part of the lesion
from the deterministic 3D U-Net segments, while the other models are capable
of producing a more accurate segmentation.

Figure 3: Example predictions for the same data slice with a nodule from the
3D U-Net, 2D & 3D Prob.U-Net in the test set. This same slice is used to get a
comparative sense of model performance.

Figure 4 depicts the predictions for the 2D and 3D Prob.U-Net for 2 slices
from a nodule in the test set. It can be seen that the 2D model misses the
nodule in Slice 34, while its 3D counterpart correctly detects it. Although
the 2D model has some uncertainty about the presence of the nodule, it is
rather low.

Figure 4: Zoomed example 3D & 2D Prob.U-Net predictions for 2 slices from a
nodule in the test set.

In Figure 5 multiple consecutive slices are depicted of a CT scan from our
test set and Prob. 3D U-Net predictions for a nodule. Slice 25 displays some
uncertainty from the model about the presence of a lesion although no rater
indicated its existence yet. In the next slice the lesion is clearly
delineated by the raters and the model captures the uncertainty in a similar
fashion as the disagreement between them. Slices 27-31 and 33-35 are not
shown, since the model correctly segments and captures the uncertainty in
comparison to the raters. Slice 38 reveals the large lesion as shown by the
annotations from the raters, but it rapidly disappears towards Slice 39.
However, the model still segments the lesion in Slice 39 and expresses high
uncertainty.

Figure 5: Example 3D Prob.U-Net predictions for a multiple slices from a
nodule in the test set.

We quantitatively compare the proposed approach with the 2D counterparts,
aiming to resolve the ambiguity in the LIDC-IDRI dataset in Table 1. We
compute the 2D $D^{2}_{\text{GED}}$ and 2D Hungarian IoU on a per-slice basis
and take the average across all the slices of the lesion, ignoring slices with
empty ground-truths and predictions. The 3D $D^{2}_{\text{GED}}$ and 3D
Hungarian IoU are immediately computed on a per-case level and then averaged
across the test set. In the case of the 2D models, during the forward pass of
a single 2D slice, an image conditional prior distribution is computed. We
then draw 16 samples from this distribution. For the next slice in the series
of the lesion, a completely new prior distribution is presented at inference
time. As such, the uncertainty captured by this prior distribution is
inconsistent over the individual slices of the lesion and it is not possible
to reconstruct a consistent and true 3D segmentation using this approach
(certainly not with 3D metrics such as 3D $D^{2}_{\text{GED}}$ and 3D IoU
calculation). For the 3D Probabilistic U-Net, we report the 2D Hungarian-
matched IoU and 2D $D^{2}_{\text{GED}}$ (2D IoU distance averaged along the
$z$-axis) and the 3D Hungarian-matched IoU and $D^{2}_{\text{GED}}$.

Model | 2D$\downarrow$ $D^{2}_{\text{GED}}$ | 2D$\uparrow$ IoU | 3D$\downarrow$ $D^{2}_{\text{GED}}$ | 3D$\uparrow$ IoU
---|---|---|---|---
Kohl et al. | 0.445 | 0.473 | N/A | N/A
Valiuddin et al. | 0.441 | 0.481 | N/A | N/A
3D U-Net | 1.283 | 0.332 | 1.263 | 0.383
3D Prob.U-Net | 0.427 | 0.510 | 0.422 | 0.457
\+ Planar Flow | 0.417 | 0.511 | 0.393 | 0.465
\+ Radial Flow | 0.429 | 0.520 | 0.401 | 0.468
Table 1: Evaluations on the LIDC-IDRI test set (15%) of the different methods
on the $D^{2}_{\text{GED}}$ and Hungarian IoU metric based on 16 samples.

We compare the inference time of the 3D U-Net and the Prob. 2D and 3D U-Net
per nodule volume ($64\times 128\times 128$ voxels). We do not include
additional results on the models with NF-augmented posteriors networks, since
the 2-step low-dimensional bijective transformation has a negligible
computational time footprint in comparison to the rest of the network. Table 2
showcases the computation time per operation for the different models and with
different batch sizes (BS). It can be seen that the Prob. 3D U-net has a
shorter inference time for a volume of this size, compared to its 2D
counterpart.

Operation | 3D U-Net(BS 1) | 2D PU-Net(BS 1) | 2D PU-Net(BS 64) | 3D PU-Net(BS 1)
---|---|---|---|---
Forward pass | 2.34 (ms) | 5.75 (ms) $\times 64$ | 397.27 (ms) | 124.31 (ms)
Sample + F-comb ($\times 1$) | N/A | 0.51 (ms) $\times 64\times 16$ | 157.09 (ms)$\times 16$ | 8.44 (ms)$\times 16$
Total | 2.34 (ms) | 892.29 (ms) | 2910.51 (ms) | 259.35 (ms)
Table 2: Inference time (per operation) of the 3D U-Net and 16 samples from
the Probabilistic 2D and 3D U-Net per nodule ($64\times 128\times 128$
voxels). BS is the Batch Size.

##  5 Discussion

This research extends the Probabilistic U-Net to the 3D domain to utilize the
rich 3D spatial information when resolving the uncertainty. We introduce the
Probabilistic 3D U-Net and employ recent improvements in the 2D Probabilistic
U-Net, by adding either a planar or radial flow to the posterior network. This
augmentation with NFs enables capturing distributions of various complexity,
thereby relaxing the strictly axis-aligned Gaussian constraint previously
employed. To test the model’s ability to capture the aleatoric uncertainty, we
use the LIDC-IDRI dataset for benchmark tests.

Section 4 displays the results of the conducted experiments. For qualitative
evaluations, Figure 3 showcases example predictions from all the models used
in this research. It can be noted that all of the models perform well on this
clearly defined lesions, except for the 3D U-Net. The 3D U-Net segments the
lesion in a conservative manner, possibly due to seeing many empty ground-
truth labels and not being able to capture this ambiguity in a meaningful way.
The Prob. 3D U-Net with planar flow expresses more uncertainty about various
parts of the lesion. Figure 4 highlights example predictions of the Prob. 3D
U-Net where information from the complete 3D volume is used to detect, segment
and resolve the uncertainty about the lesion in axial Slice 34 of the CT scan.
The same nodule is missed by the Prob. 2D U-Net, due to a lack of information
from prior slices. These results are also reflected in the 2D
$D^{2}_{\text{GED}}$ and Hungarian IoU, as shown in Table 1, with the 3D
models outperforming the 2D models. Interestingly, in Figure 5, the 3D spatial
awareness of the model is showcased through the uncertainty expressed in Slice
25. A rather large lesion is coming up (iterating through the CT slices in an
ascending order) and the model expresses uncertainty about the exact starting
position, since the raters segmented no lesion followed by a large lesion in
consecutive slices. The same phenomenon can be seen moving from Slice 38 to
39, although here the model incorrectly (according to the raters) segments the
lesion while it is still partially visible.

By the performance evaluations presented in Table 2, it can be observed that
it is the most efficient to present the uncertainty with a 3D Prob. U-Net. For
a volume of 64$\times$128$\times$128 voxels, the 64 2D slices can be passed
through the Prob. 2D U-Net in a large batch, but it scales poorly in
comparison to a single-slice forward pass. The forward pass of the Prob. 2D
U-Net with a batch size of 64 takes 397.27 ms to compute. Drawing 16 samples
from the Prior distribution and combining it with the 2D U-Net features
through the Feature combination networks takes 2513.44 ms (157.09
ms$\times$16). In total this approach will take 2910.51 ms to compute compared
to the approximate 10$\times$ speed improvement that is required for
calculating the uncertainty with the Prob. 3D U-Net (259.35 ms). It should be
noted that significant computational time drawbacks occur when using the Prob.
3D U-Net in comparison to the the standard 3D U-Net (2.34 ms for inference),
although this is not a realistic alternative since no uncertainty can be
expressed.

##  6 Conclusion

In CAD methods, it is important to provide clinicians with an accurate measure
of uncertainty when they evaluate and plan their procedures. Accurately
capturing and presenting segmentation uncertainty will increase clinician
confidence in model predictions and facilitates better informed decision-
making. Existing CT-based segmentation methods aim to do so by quantifying the
uncertainty from 2D image slices, whereas the true uncertainty resides in the
full 3D CT or MRI volume. We propose a novel 3D probabilistic segmentation
model that is capable of resolving and presenting the aleatoric uncertainty in
3D volumes through diverse and plausible nodule segmentations. The model
consists of a Deep 3D U-Net and a 3D conditional VAE that is augmented with an
Normalizing Flow (NF) in the posterior network. NFs allow for more flexible
distribution modelling and as such, we alleviate the strictly Gaussian
posterior distribution that was previously enforced. We test our approach on
the LIDC-IDRI lung nodule CT dataset. This is, to the best of our knowledge,
the first approach that presents the 3D Squared Generalized Energy Distance
($D^{2}_{\text{GED}}$) and 3D Hungarian-matched IoU for lung nodule
segmentation and uncertainty prediction. We quantify the uncertainty
prediction performance and achieve a 0.401 3D $D^{2}_{\text{GED}}$ and a
Hungarian-matched 3D IoU of 0.468 with the radial Prob. 3D U-Net. In addition,
since the model uses the full native 3D volumes, it is a step closer to the
practical application of accurately segmenting and presenting uncertainty in
3D CT data. Finally, we present the aleatoric uncertainty, computed as the
standard deviation across the model predictions, in a visual manner. This
enables an interpretable expression of the uncertainty and is potentially
providing clinicians additional insight into data ambiguity and allowing for
more informed decision-making.

## References

  * [1] Kendall, A. and Gal, Y., “What uncertainties do we need in bayesian deep learning for computer vision?,” CoRR (2017).
  * [2] Armato III, S., Mclennan, G., Bidaut, L., McNitt-Gray, M., Meyer, C., Reeves, A., Zhao, B., Aberle, D., Henschke, C., Hoffman, E., Kazerooni, E., Macmahon, H., Beek, E., Yankelevitz, D., Biancardi, A., Bland, P., Brown, M., Engelmann, R., Laderach, G., and Clarke, L., “The lung image database consortium (lidc) and image database resource initiative (idri): A completed reference database of lung nodules on ct scans,” Medical Physics 38, 915–931 (01 2011).
  * [3] Valiuddin, M., Viviers, C. G., van Sloun, R. J., and Sommen, F. v. d., “Improving aleatoric uncertainty quantification in multi-annotated medical image segmentation with normalizing flows,” 75–88, Springer (2021).
  * [4] Kohl, S. A., Romera-Paredes, B., Meyer, C., De Fauw, J., Ledsam, J. R., Maier-Hein, K. H., Eslami, S., Rezende, D. J., and Ronneberger, O., “A probabilistic u-net for segmentation of ambiguous images,” arXiv preprint arXiv:1806.05034 (2018).
  * [5] Selvan, R., Faye, F., Middleton, J., and Pai, A., “Uncertainty quantification in medical image segmentation with normalizing flows,” (2020).
  * [6] Hu, S., Worrall, D., Knegt, S., Veeling, B., Huisman, H., and Welling, M., “Supervised uncertainty quantification for segmentation with multiple annotations,” (2019).
  * [7] Rezende, D. and Mohamed, S., “Variational inference with normalizing flows,” in [International conference on machine learning ], 1530–1538, PMLR (2015).
  * [8] Ronneberger, O., Fischer, P., and Brox, T., “U-net: Convolutional networks for biomedical image segmentation,” in [Medical Image Computing and Computer-Assisted Intervention – MICCAI 2015 ], Navab, N., Hornegger, J., Wells, W. M., and Frangi, A. F., eds., 234–241, Springer International Publishing, Cham (2015).
  * [9] Çiçek, Ö., Abdulkadir, A., Lienkamp, S. S., Brox, T., and Ronneberger, O., “3d u-net: Learning dense volumetric segmentation from sparse annotation,” in [Medical Image Computing and Computer-Assisted Intervention – MICCAI 2016 ], Ourselin, S., Joskowicz, L., Sabuncu, M. R., Unal, G., and Wells, W., eds., 424–432, Springer International Publishing, Cham (2016).
  * [10] Kobyzev, I., Prince, S., and Brubaker, M., “Normalizing flows: An introduction and review of current methods,” IEEE Transactions on Pattern Analysis and Machine Intelligence , 1–1 (2020).
  * [11] Kingma, D. P. and Welling, M., “Auto-encoding variational bayes,” (2014).
  * [12] Van Den Berg, R., Hasenclever, L., Tomczak, J., and Welling, M., “Sylvester normalizing flows for variational inference,” in [34th Conference on Uncertainty in Artificial Intelligence 2018, UAI 2018 ], Globerson, A., Globerson, A., and Silva, R., eds., 34th Conference on Uncertainty in Artificial Intelligence 2018, UAI 2018, 393–402, Association For Uncertainty in Artificial Intelligence (AUAI) (Jan. 2018).  34th Conference on Uncertainty in Artificial Intelligence 2018, UAI 2018 ; Conference date: 06-08-2018 Through 10-08-2018.
  * [13] Fu, H., Li, C., Liu, X., Gao, J., Celikyilmaz, A., and Carin, L., “Cyclical annealing schedule: A simple approach to mitigating kl vanishing,” arXiv preprint arXiv:1903.10145 (2019).
  * [14] Alemi, A., Poole, B., Fischer, I., Dillon, J., Saurous, R. A., and Murphy, K., “Fixing a broken elbo,” in [International conference on machine learning ], 159–168, PMLR (2018).
  * [15] Wolny, A., Cerrone, L., Vijayan, A., Tofanelli, R., Barro, A. V., Louveaux, M., Wenzl, C., Strauss, S., Wilson-Sánchez, D., Lymbouridou, R., Steigleder, S. S., Pape, C., Bailoni, A., Duran-Nebreda, S., Bassel, G. W., Lohmann, J. U., Tsiantis, M., Hamprecht, F. A., Schneitz, K., Maizel, A., and Kreshuk, A., “Accurate and versatile 3d segmentation of plant tissues at cellular resolution,” eLife 9, e57613 (jul 2020).