# Geometric control of algebraic systems

Benoît Legat  Raphaël M. Jungers<EMAIL_ADDRESS>ICTEAM institute,
UCLouvain, Louvain-la-Neuve, Belgium. (e-mail: ).
<EMAIL_ADDRESS>ICTEAM institute, UCLouvain, Louvain-la-Neuve,
Belgium. (e-mail: ).

###### Abstract

In this paper, we present a geometric approach for computing the controlled
invariant set of a continuous-time control system. While the problem is well
studied for in the ellipsoidal case, this family is quite conservative for
constrained or switched linear systems. We reformulate the invariance of a set
as an inequality for its support function that is valid for any convex set.
This produces novel algebraic conditions for the invariance of sets with
polynomial or piecewise quadratic support function. We compare it with the
common algebraic approach for polynomial sublevel sets and show that it is
significantly more conservative than our method.

###### keywords:

Lyapunov methods; Control of constrained systems; Control of switched systems;
Convex optimization.

††thanks: RJ is a FNRS honorary Research Associate. This project has received
funding from the European Research Council (ERC) under the European Union’s
Horizon 2020 research and innovation programme under grant agreement No 864017
- L2C. RJ is also supported by the Walloon Region, the Innoviris Foundation,
and the FNRS (Chist-Era Druid-net).

##  1 Introduction

Computing111A preliminary version of this work appears in Legat (2020).
controlled invariant set is paramount in many applications (Blanchini and
Miani (2015)). Indeed, the existence of a controlled invariant set is
equivalent to the stabilizability of a control system (Sontag (1983)) and a
(possibly nonlinear) stabilizable state feedback can be deduced from the
controlled invariant set (Barmish (1985)).

The stabilizability of a linear time-invariant (LTI) control system is
equivalent to the stability of its uncontrollable subspace (which is readily
accessible in its Controllability Form) (Wonham, 1985, Section 2.4). Indeed,
the eigenvalues of its controllable subspace can be fixed to any value by a
proper choice of linear state feedback. The resulting controlled system is
stable hence an invariant ellipsoid can be determined by solving a system of
linear equations (Liapounoff (1907)). This set is also controlled invariant
for the control system. When a control system admits an ellipsoidal controlled
invariant set, it is said to be _quadratically stabilizable_. When there
exists a linear state feedback such that the resulting autonomous system
admits an ellipsoidal invariant set, it is said to be _quadratically
stabilizable via linear control_.

While the stabilizability of LTI control systems is equivalent to their
quadratic stabilizability via linear control, it is no longer the case for
_uncertain_ or _switched_ systems (Petersen (1985)). Furthermore, it is often
desirable for _constrained_ systems to find a controlled invariant set of
maximal volume (or which is maximal in some direction (Ahmadi and Gunluk
(2018))). For such problem, the method detailed above is not suitable as it
does not take any volume consideration but more importantly, the maximal
volume invariant set may not be an ellipsoid and may not be rendered stable
via a linear control. For this reason, the Linear Matrix Inequality (LMI) (4)
was devised to encapsulate the controlled invariance of an ellipsoid via
linear control (Boyd et al., 1994, Section 7.2.2) and the conservatism of the
choice of linear control was analysed (Sontag (1983)). As the linearity of the
control was found to be conservative for uncertain systems (Petersen (1985)),
the LMI (13) was found to encapsulate controlled invariance of an ellipsoid
via _any_ state-feedback Barmish (1985).

These LMIs have had a tremendous impact on control, but unfortunately the
approach is limited to ellipsoids due to its algebraic nature. We reinterpret
it in a geometric/behavioural framework, based on convex analysis, which
allows us to formulate a general condition for the controlled invariance of
arbitrary convex sets via any state-feedback in Theorem 7. While this
condition reduces to (13) for the special case of ellipsoids, it provides a
new methods for computing controlled invariant convex sets with convex
polynomial and piecewise quadratic support functions.

In Section 2, we review the classical LMI’s for the invariance and controlled
invariance of ellipsoids and discusses the challenges for its generalization
to sublevel sets of polynomials. In Section 3, we develop a generic condition
of control invariance for continuous-time systems using our geometric
approach. We particularize it for ellipsoids (resp. sets with polynomial and
piecewise quadratic support functions) in Section 3.1 (resp. Section 3.2 and
Section 3.3). We illustrate these new results with a numerical example in
Section 4.

#### Reproducibility

The code used to obtain the results is published on codeocean (Legat and
Jungers (2021)). The set programs are reformulated by SetProg (Legat et al.
(2019)) into a Sum-of-Squares program which is reformulated into a
semidefinite program by SumOfSquares (Weisser et al. (2019)) which is solved
by Mosek v8 (ApS (2017)) through MathOptInterface (Legat et al. (2020a)).

##  2 Algebraic approach

Computing an ellipsoidal invariant set for an autonomous system $\dot{x}=Ax$
where $A\in\mathbb{R}^{n\times n}$ consists in searching for an ellipsoidal
set

$\mathcal{E}_{P}=\\{\,x\in\mathbb{R}^{n}\mid x^{\top}Px\leq 1\,\\}$

that satisfies the _Nagumo condition_ (Blanchini and Miani, 2015, Theorem
4.7): $x^{\top}PAx\leq 0$ for any $x\in\mathbb{R}^{n}$. The Nagumo condition
for ellipsoids is equivalent to the LMI:

$A^{\top}P+PA\preceq 0.$ (1)

which allows to search for ellipsoidal invariant sets using semidefinite
programming.

Consider the continuous-time control linear system

$\dot{x}=Ax+Bu$ (2)

where $A\in\mathbb{R}^{n_{x}\times n_{x}}$ and $B\in\mathbb{R}^{n_{x}\times
n_{u}}$ with the following definition of invariance.

######  Definition 1 (Controlled invariant set)

A set $S$ is _controlled invariant_ for system (2) if for any state $x_{0}\in
S$, there exists a control $u(t)$ such that the trajectory with initial state
and control $u$ remains in $S$.

With the presence of the control $u$ in the system, the Nagumo condition
becomes:

$\forall x\in\mathbb{R}^{n_{x}},\exists
u\in\mathbb{R}^{n_{u}},x^{\top}P(Ax+Bu)\leq 0.$ (3)

The control term $u$, or more precisely the existential quantifier $\exists$
prevents us to transform this into an LMI directly.

Fixing the control to a linear state feedback $u(x)=Kx$ for some matrix $K$
allows to fallback to the case of autonomous system $\dot{x}=(A+BK)x$. Then,
the invariance condition can be formulated as the _Bilinear Matrix Inequality_
(BMI):

$A^{\top}P+PA+K^{\top}B^{\top}P+PBK\preceq 0$

While the matrix inequality is bilinear in $K$ and $P$, and BMI’s are NP-hard
to solve in general (Toker and Ozbay (1995)), a clever algebraic manipulation
allows to reformulate it as a Linear Matrix Inequality (LMI) in $Q:=P^{-1}$
and $Y:=KQ$, where the sought controlled invariant ellipsoid is given by
$\mathcal{E}_{P}$, see e.g. (Boyd et al., 1994, Section 7.2.1, Section 7.2.2)
and (Blanchini and Miani, 2015, Section 4.4.1) for more details. The linear
matrix inequality is

$QA^{\top}+AQ+Y^{\top}B^{\top}+BY\preceq 0.$ (4)

Because the algebraic manipulation which allows to reformulate the BMI into an
LMI is done at the level of matrices, it is not clear how this approach can be
generalized to other families of sets such as the ones considered in Section
3.2 and Section 3.3. Moreover, searching for ellipsoidal controlled invariant
sets may be rather restrictive and the conservativeness is amplified for the
class of hybrid systems.

One attempt to generalize it to controlled invariant sublevel sets of
polynomials of degree $d$ is developed in Prajna et al. (2004). While the
method allows to consider systems defined by polynomial equations, we show
below that for linear systems, it has significant restrictions for $d>2$ and
the case $d=2$ reduces to the ellipsoidal case given by (4). This suggests
that the methods is essentially restricted to compute invariant sublevel sets
of polynomials of a most twice the degree of the polynomial equations. On the
other hand, the condition we find in Corollary 10 is necessary and sufficient.
The remaining of this section particularize the approach for linear systems,
more details can be found in Prajna et al. (2004).

Let $x_{[d]}$ denote the vector of all monomials of degree $d$ with the
variables $x_{i}$. Consider the set $\\{\,x\mid x_{[d]}^{\top}Px_{[d]}\leq
1\,\\}$ for a symmetric positive definite matrix $P$ and a state feedback of
the form $u(x)=K_{1}(x)x+K_{2}(x)x_{[d]}$ where $K_{1}(x),K_{2}(x)$ are
matrices of the appropriate dimensions whose entries are polynomials in $x$.
The invariance of the set for the autonomous system
$\dot{x}=Ax+BK_{1}(x)x+BK_{2}(x)x_{[d]}$ is equivalent to the nonnegativity of
the polynomial

$x_{[d]}^{\top}PM(x)BK_{2}(x)x_{[d]}+x_{[d]}^{\top}PF(x)x$

where $F(x)=M(x)BK_{1}(x)+M(x)A$ and $M(x)$ is the jacobian of the
transformation $x\mapsto x_{[d]}$. This polynomial can be rewritten in the
matrix form:

$\begin{bmatrix}x_{[d]}\\\
x\end{bmatrix}^{\top}\begin{bmatrix}PM(x)BK_{2}(x)&PF(x)\\\
0&0\end{bmatrix}\begin{bmatrix}x_{[d]}\\\ x\end{bmatrix}.$ (5)

The following is therefore a sufficient condition for the invariance:

$\forall x\in\mathbb{R}^{n},\quad\begin{bmatrix}2PM(x)BK_{2}(x)&PF(x)\\\
F^{\top}(x)P&0\end{bmatrix}\preceq 0.$ (6)

After a similarity transformation with the block diagonal matrix
$\operatorname{BlockDiag}(P^{-1},I_{n_{x}})$ where
$I_{n_{x}}\in\mathbb{R}^{n_{x}\times n_{x}}$ is the identity matrix, the
condition is rewritten as:

$\forall x\in\mathbb{R}^{n},\quad\begin{bmatrix}2M(x)BY(x)&F(x)\\\
F^{\top}(x)&0\end{bmatrix}\preceq 0$ (7)

where $Y(x)=K_{2}(x)P^{-1}$.

By Proposition 14, for (7) to hold we need all entries of $F(x)$ to be zero
polynomials which is quite restictive. The example below discusses the
conservativeness of (5) compared to (6).

###### Example 2

Consider the autonomous system $\dot{x}=-x$ and the invariant set $\\{\,x\mid
x^{4}\leq 1\,\\}$. The condition (5) is satisfied:

$\begin{bmatrix}x^{2}\\\ x\end{bmatrix}^{\top}\begin{bmatrix}0&-2x\\\
0&0\end{bmatrix}\begin{bmatrix}x^{2}\\\ x\end{bmatrix}=-x^{4}$

while the matrix

$\begin{bmatrix}0&-2x\\\ -2x&0\end{bmatrix}$

of conditions (6) and (7) is indefinite for any nonzero $x$.

##  3 Geometric approach

In this section we derive a characterization of the controlled invariance of a
closed convex set under the form of an inequality for its support function. We
start by showing the equivalence of the notion of invariance with another
class of systems that directly models the geometric behaviours of the
trajectories of control systems with unconstrained input.

Consider the continuous-time _algebraic_ linear systems:

$E\dot{x}=Cx.$ (8)

with the following definition of invariance.

######  Definition 3 (Invariant set for an algebraic system)

A set $S$ is _invariant_ for system (8) if for any state $x_{0}\in S$, there
exists a trajectory of the system that remains in $S$.

Note the use of “there exists” instead of “for all” in the definition of
invariance as both versions exists; see (Legat et al., 2020b, Remark 4) for
more details.

###### Proposition 4

Given a subset $\mathcal{S}\subseteq\mathbb{R}^{n}$ and matrices
$A\in\mathbb{R}^{r\times n},B\in\mathbb{R}^{r\times m}$, the following holds:

$A\mathcal{S}+B\mathbb{R}^{m}=\pi_{\operatorname{Im}(B)^{\perp}}{}^{-1}\pi_{\operatorname{Im}(B)^{\perp}}{}A\mathcal{S}$

where $\pi_{\operatorname{Im}(B)^{\perp}}{}$ is any projection matrix onto the
orthogonal subspace of $\operatorname{Im}(B)$.

###### Proof 3.1

Given $x\in\mathcal{S}$ and $y\in\mathbb{R}^{r}$, we have $y\in
A\\{x\\}+B\mathbb{R}^{m}$ if and only if $y-Ax\in\operatorname{Im}(B)$ or
equivalently,
$\pi_{\operatorname{Im}(B)^{\perp}}y=\pi_{\operatorname{Im}(B)^{\perp}}Ax$.

###### Proposition 5

Consider a control system (2) and an arbitrary projection matrix
$\pi_{\operatorname{Im}(B)^{\perp}}{}$ onto the orthogonal subspace of
$\operatorname{Im}(B)$. A set $S$ is controlled invariant for the control
system (2) with $\mathcal{U}=\mathbb{R}^{n_{u}}$, as defined in Definition 1,
if and only if it is controlled invariant for the algebraic system

$\pi_{\operatorname{Im}(B)^{\perp}}{}\dot{x}=\pi_{\operatorname{Im}(B)^{\perp}}{}Ax,$

as defined in Definition 3.

###### Proof 3.2

By Proposition 4, there exists $u\in\mathbb{R}^{n_{u}}$ such that
$\dot{x}=Ax+Bu$ if and only if
$\pi_{\operatorname{Im}(B)^{\perp}}{}\dot{x}=\pi_{\operatorname{Im}(B)^{\perp}}{}Ax$.
As the input $u$ is unconstrained, the result follows.

The Nagumo condition for algebraic systems has the following form.

###### Proposition 6

A closed set $\mathcal{S}$ is invariant for system (8), as defined in
Definition 3, if and only if

$\forall x\in\partial\mathcal{S},\exists y\in T_{S}(x),Ey=Ax.$ (9)

See Section B for a brief review of the concepts of convex geometry needed for
the remaining of this section. The invariance condition (9) can be rewritten
in terms of _exposed faces_.

######  Theorem 7 (Controlled invariance of convex set)

A convex set $\mathcal{C}$ is invariant for system (8) with matrices
$C,E\in\mathbb{R}^{r\times n}$ if and only if

$\forall z\in\mathbb{R}^{r},\forall x\in F_{\mathcal{C}}(E^{\top}z),\langle
z,Cx\rangle\leq 0.$ (10)

###### Proof 3.3

As $\mathcal{C}$ is convex, $T_{\mathcal{C}}(x)$ is a convex cone. By
definition of the polar of a cone, $x\in ET_{\mathcal{C}}(x)$ if and only if
$\langle y,x\rangle\leq 0$ for all $y\in[ET_{\mathcal{C}}(x)]^{\circ}$. By
Proposition 20, $[ET_{\mathcal{C}}(x)]^{\circ}=E^{-\top}N_{\mathcal{C}}(x)$.
Therefore, the set $\mathcal{C}$ is invariant if and only if

$\forall x\in\partial\mathcal{C},\forall z\in
E^{-\top}N_{\mathcal{C}}(x),\langle z,Cx\rangle\leq 0.$ (11)

By Proposition 19, we have

$\\{\,(x,z)\in\partial\mathcal{C}\times\mathbb{R}^{r}\mid E^{\top}z\in
N_{\mathcal{C}}(x)\,\\}=\\\
\\{\,(x,z)\in\partial\mathcal{C}\times\mathbb{R}^{r}\mid x\in
F_{\mathcal{C}}(E^{\top}z)\,\\}.$

As we show in the remaining of this section, Theorem 7 allows to reformulate
the invariance as an inequality in terms of the support function. This allows
to combine the invariance constraint with other set constraints that can be
formulated in terms of support functions. Moreover, for an appropriate family
of sets, also called _template_ , the set program can be automatically
rewritten into a convex program combining all constraints using _set
programming_ Legat et al. (2019); Legat (2020). For this reason, we only focus
on the invariance constraint and do not detail how to formulate the complete
convex programs with the objective and all the constraints needed to obtain
the results of Section 4 as these problems are decoupled.

This allows to formulate the invariance of a convex set in terms of its
support function if it is differentiable. We generalize this result with a
relaxed notion of differentiability in Theorem 12.

###### Theorem 8

Consider a nonempty closed convex set $\mathcal{C}$ such that
$\delta^{*}(\cdot|\mathcal{C})$ is differentiable. Then $\mathcal{C}$ is
invariant for system (8) with matrices $C,E\in\mathbb{R}^{r\times n}$ if and
only if

$\forall z\in\mathbb{R}^{r},\langle
z,C\nabla\delta^{*}(E^{\top}z|\mathcal{C})\rangle\leq 0.$ (12)

###### Proof 3.4

By Proposition 21,
$F_{\mathcal{C}}(E^{\top}z)=\\{\nabla\delta^{*}(E^{\top}z|\mathcal{C})\\}$
hence (10) is equivalent to (12).

###  3.1 Ellipsoidal controlled invariant set

In this section, we particularize Theorem 8 to the case of ellipsoids. Since
the support function of an ellipsoid $\mathcal{E}_{P}$ is
$\delta^{*}(y|\mathcal{E}_{P})=\sqrt{y^{\top}P^{-1}y}$, we have the following
corollary of Theorem 8.

######  Corollary 9 (Barmish (1985))

Given a positive definite matrix $P$, the ellipsoid $\mathcal{E}_{P}$ is
controlled invariant for system (8) if and only if

$CP^{-1}E^{\top}+EP^{-1}C^{\top}\preceq 0.$ (13)

Observe that for the trivial case $\operatorname{Im}(B)=\mathbb{R}^{n}$ for
system (2), Proposition 5 produces a system (8) with $r=0$ hence the LMI (13)
will be trivially satisfied for any $P^{-1}$, which is expected.

In comparison to (4), for a system (8) with matrices
$C,E\in\mathbb{R}^{r\times n}$, the LMI (4) has size $n\times n$ while the LMI
(13) has only size $r\times r$. The characterization of controlled invariance
of ellipsoids using (13) can also be obtained by applying an elimination
procedure to reduce (4); see (Boyd et al., 1994, Equation (7.11)). However,
uncertain or switched system may need a nonlinear state feedback to be
quadratically stabilizabilized Petersen (1985). For such systems, (4) is
conservative since it assumes a linear feedback while (13) does not assume
anything about the feedback. It was shown in Barmish (1985) that if (13) is
satisfied then a stabilizing nonlinear continuous state feedback can be
deduced from the solution $P$. There is even a closed form for the feedback in
case of single input (Barmish, 1985, Eq. (15)).

###  3.2 Polynomial controlled invariant set

In this section, we derive the algebraic condition for the controlled
invariance of a set with polynomial support function. This template is
referred to as _polyset_ ; see (Legat, 2020, Section 1.5.3).

###### Corollary 10

Given a homogeneous222A polynomial is _homogeneous_ if all its monomials have
the same total degree nonnegative polynomial $p(x)$ of degree $2d$, the set
$\mathcal{C}$ defined by the support function
$\delta^{*}(y|\mathcal{C})=p(y)^{\frac{1}{2d}}$ is invariant for system (8)
with matrices $C,E\in\mathbb{R}^{r\times n}$ if and only if the polynomial

$z^{\top}C\nabla p(E^{\top}z)$ (14)

is nonpositive for all $z\in\mathbb{R}^{r}$.

###### Proof 3.5

We have

$\nabla\delta^{*}(y|\mathcal{C})=\frac{1}{p(y)^{1-\frac{1}{2d}}}\nabla p(y).$

If $p(y)$ is identically zero, this is trivially satisfied. Otherwise,
$p(y)^{1-\frac{1}{2d}}$ is nonnegative and is zero in an algebraic variety of
dimension $n-1$ at most. Therefore, (12) is equivalent to (14).

While verifying the nonnegativity of a polynomial is co-NP-hard, a sufficient
condition can be obtained via the standard Sum-of-Squares programming
framework; see Blekherman et al. (2012).

###  3.3 Piecewise semi-ellipsoidal controlled invariant set

Johansson and Rantzer (1998) study the computation of piecewise quadratic
Lyapunov functions for continuous-time autonomous piecewise affine systems.
Legat et al. (2021) present a convex programming approach to compute
_piecewise semi-ellipsoidal_ controlled invariant sets of discrete-time
control systems. In this section, we show that Theorem 8 provides the
corresponding condition for continuous-time.

We recall (Legat et al., 2021, Definition 2) below.

###### Definition 11

A _polyhedral conic partition_ of $\mathbb{R}^{n}$ is a set of $m$ polyhedral
cones $(\mathcal{P}_{i})_{i=1}^{m}$ with nonempty interior such that for all
$i\neq j$, $\dim(\mathcal{P}_{i}\cap\mathcal{P}_{j})<n$ and
$\cup_{i=1}^{m}\mathcal{P}_{i}=\mathbb{R}^{n}$.

Piecewise semi-ellipsoids have a support function of the form

$\delta^{*}(y|\mathcal{C})=\sqrt{y^{\top}Q_{i}y}\qquad y\in\mathcal{P}_{i}$
(15)

where $(\mathcal{P}_{i})_{i=1}^{m}$ is a polyhedral conic partition. The
support function additionally has to satisfy (Legat et al., 2021, (2) and (3))
to ensure its continuity and convexity.

###### Theorem 12

Consider a polyhedral conic partition $(\mathcal{P}_{i})_{i=1}^{m}$ and a
nonempty closed convex set $\mathcal{C}$ such that

$\delta^{*}(y|\mathcal{C})=f_{i}(y)\qquad y\in\mathcal{P}_{i}$

for differentiable functions $f_{i}:\mathcal{P}_{i}\to\mathbb{R}$. The set
$\mathcal{C}$ is invariant for system (8) with matrices
$C,E\in\mathbb{R}^{r\times n}$ if and only if, for all $i=1,\ldots,m$ and
$z\in\mathbb{R}^{r}$ such that $E^{\top}z\in\mathcal{P}_{i}$, we have

$\langle z,C\nabla f_{i}(E^{\top}z)\rangle\leq 0.$ (16)

###### Proof 3.6

Given $z\in\mathbb{R}^{r}$ such that $E^{\top}z$ is in the intersection of the
boundary of $\mathcal{C}$ and the interior of $\mathcal{P}_{i}$, the support
function is differentiable at $E^{\top}z$ hence, by Proposition 21,
$F_{\mathcal{C}}(E^{\top}z)=\\{\nabla f_{i}(E^{\top}z)\\}$. The condition (10)
is therefore reformulated as (16).

Given a subset $I$ of $\\{1,\ldots,m\\}$ and $z\in\mathbb{R}^{r}$ such that
$E^{\top}z$ is in the intersection of the boundary of $\mathcal{C}$ and
$\cap_{i\in I}\mathcal{P}_{i}$, $F_{\mathcal{C}}(E^{\top}z)$ is the convex
hull of $\nabla\delta^{*}(E^{\top}z|\mathcal{C})$ for each $i\in I$. For any
convex combination (i.e., nonnegative numbers summing to 1)
$(\lambda_{i})_{i\in I}$, (16) implies that

$\langle z,C\sum_{i\in I}\lambda_{i}\nabla f_{i}(E^{\top}z)\rangle=\sum_{i\in
I}\lambda_{i}\langle z,C\nabla f_{i}(E^{\top}z)\rangle\leq 0.$

###### Corollary 13

A piecewise semi-ellipsoid $\mathcal{C}$ is invariant for system (8) with
matrices $C,E\in\mathbb{R}^{r\times n}$ if and only if the quadratic form

$z^{\top}CP_{i}^{-1}E^{\top}z+z^{\top}EP_{i}^{-1}C^{\top}z$ (17)

is nonpositive for all $i=1,\ldots,m$ and $z\in\mathbb{R}^{r}$ such that
$E^{\top}z\in\mathcal{P}_{i}$.

The condition (17) amounts to verifying the positive semidefiniteness of a
quadratic form when restricted to a polyhedral cone. When this cone is the
positive orthant, this is called the _copositivity_ which is co-NP-complete to
decide (Murty and Kabadi (1987)). However, a sufficient LMI is given in (Legat
et al., 2021, Proposition 2) and a necessary and sufficient condition is given
by a hierarchy of Sum-of-Squares programs (Parrilo, 2000, Chapter 5).

##  4 Numerical example

In this section, we study a simple numerical example to illustrate our new
approach. Suppose we want to compute a controlled invariant set for the
control system

$\dot{x}=\begin{pmatrix}0&1\\\ 0&0\\\ \end{pmatrix}x+\begin{pmatrix}0\\\
1\end{pmatrix}u$

with the state constraint $x\in[-1,1]^{2}$ and input constraint $u\in[-1,1]$.
We represent the state set $[-1,1]^{2}$ and its polar in green in Figure 1 and
Figure 2.

The union of controlled invariant sets is controlled invariant. Moreover, by
linearity, the convex hull of the unions of controlled invariant sets is
controlled invariant. Therefore, there exists a _maximal_ controlled
invariant, i.e., a controlled invariant set in which all controlled invariant
sets are included, for any family that is closed under union (resp. convex
hull); it is the union (resp. convex hull) of all controlled invariant sets
included in $[-1,1]^{2}$.

For this simple planar system, the maximal controlled invariant set can be
obtained by hand, it is

$\\{\,(x_{1},x_{2})\in[-1,1]^{2}\mid x_{1}x_{2}\leq 0\text{ or }|x_{1}|\leq
1-x_{2}^{2}/2\,\\}$

Its polar is given by

$\displaystyle\,\\{\,x\mid x_{1}x_{2}\leq 0\text{ and }|x_{1}-x_{2}|\leq
1\,\\}$ $\displaystyle\cup$ $\displaystyle\,\\{\,x\mid x_{1}(x_{1}-x_{2})\leq
0\text{ and }|x_{1}/2+x_{2}|\leq 1\,\\}$ $\displaystyle\cup$
$\displaystyle\,\\{\,x\mid x_{2}(x_{2}-x_{1})\leq 0\text{ and
}(2x_{1}-\operatorname*{sign}(x_{1}))^{2}+2x_{2}^{2}\leq 1\,\\}$

We represent it in yellow in Figure 1 and Figure 2.

As Proposition 5 requires the input to be unconstrained, we will consider
projections onto the first two dimensions of controlled invariant sets of the
following control system:

$\dot{x}=\begin{pmatrix}0&1&0\\\ 0&0&1\\\
0&0&0\end{pmatrix}x+\begin{pmatrix}0\\\ 0\\\ 1\end{pmatrix}u.$

with the state constraint $x\in[-1,1]^{3}$ and no input constraint.

Following Proposition 5, we consider the algebraic system

$\begin{pmatrix}1&0&0\\\ 0&1&0\\\ \end{pmatrix}\dot{x}=\begin{pmatrix}0&1&0\\\
0&0&1\\\ \end{pmatrix}x$

with the state constraint $x\in[-1,1]^{3}$.

While the _maximal_ invariant set is well defined, it is not the case anymore
when we restrict the set to belong to the family of ellipsoids, polysets or
piecewise semi-ellipsoids for a fixed polyhedral conic partition as these
families are not invariant under union nor convex hull. The objective used to
determine which invariant set is selected depends on the particular
application. Let $\mathcal{D}$ be the convex hull of
$\\{(-1+\sqrt{3},-1+\sqrt{3}),(-1,1),(1-\sqrt{3},1-\sqrt{3}),(1,-1)\\}$. For
this example, we maximize $\gamma$ such that $\gamma\mathcal{D}$ is included
in the projection of the invariant set onto the first two dimensions. We
represent $\gamma\mathcal{D}$ in red in Figure 1 and Figure 2.

For the ellipsoidal template considered in Section 3.1, the optimal solution
is shown in Figure 1 as ellipsoids corresponds to polysets of degree 2. The
optimal objective value is $\gamma\approx 0.81$.

For the polyset template considered in Section 3.2, the optimal solution are
represented in Figure 1. The optimal objective value for degree 4 (resp. 6, 10
and 20) is $\gamma\approx 0.91$. (resp. $\gamma\approx 0.93$, $\gamma\approx
0.96$ and $\gamma\approx 0.98$).

For the piecewise semi-ellipsoidal template, we consider polyhedral conic
partitions made of the conic hull of each facet of the polytope with extreme
points

$(\cos(\alpha)\cos(\beta),\sin(\alpha)\cos(\beta),\sin(\beta))$ (18)

where $\alpha=0,2\pi/m_{1},4\pi/m_{1},\ldots,2(m_{1}-1)\pi/m_{1}$ and
$\beta=-\pi/2,\ldots,-2\pi/(m_{2}-1),-\pi/(m_{2}-1),0,\pi/(m_{2}-1),2\pi/(m_{2}-1),\ldots,\pi/2$.

The optimal objective value for $m=(4,3)$ (resp. $(8,5)$) is $\gamma\approx
0.89$ (resp. $\gamma\approx 0.92$). The corresponding optimal solution is
shown in Figure 2.

Figure 1: In blue are the solution for polysets of different degrees. The
degrees from top to bottom are respectively 2, 4, 6, 10 and 20. The green set
is the safe set $[-1,1]^{2}$, the yellow set is the maximal controlled
invariant set and the red set is $\gamma\mathcal{D}$. The sets are represented
in the primal space in left figures and in polar space in the right figures.

Figure 2: In blue are the solution for piecewise semi-ellipsoids for two
different polyhedral conic partitions. The partitions from top to bottom are
as described in (18) with $m=(4,3)$ (resp. $(8,5)$). The green set is the safe
set $[-1,1]^{2}$, the yellow set is the maximal controlled invariant set and
the red set is $\gamma\mathcal{D}$. The sets are represented in the primal
space in left figures and in polar space in the right figures.

##  5 Conclusion

We proved a condition for continuous-time controlled invariance of a convex
set based on its support function. We particularized the condition for three
templates: ellipsoids, polysets and piecewise semi-ellipsoids. In the
ellipsoidal case, it reduces to a known LMI, in the polyset case, it provides
a condition significantly less conservative than the existing one333Indeed, it
is equivalent to invariance by Corollary 10 and we showed in Section 2 that
the existing approach is quite conservative. and in the piecewise semi-
ellipsoidal case, it provides the first convex programming approach for
continuous-time controlled invariance to the best of our knowledge.

As future work, we aim to apply this framework to other families such as the
_piecewise polysets_ defined in Legat (2020). Moreover, instead of considering
a uniform discretization of the hypersphere as in (18), a more adaptive
methods could be considered. The sensitivity information provided by the dual
solution of the optimization program could for instance determine which pieces
of the partition should be refined.

Finally, as the discrete-time version of this work developed in Legat et al.
(2020b, 2021) also requires the set to be represented by its support function
for the optimization program to be convex, both these methods and the method
of this paper could be combined to compute controlled invariant sets for
hybrid automata using the condition of this paper for the invariance subject
to the dynamics of each mode and the condition of Legat et al. (2020b, 2021)
for the invariance subject to each reset map.

## References

  * Ahmadi and Gunluk (2018) Ahmadi, A.A. and Gunluk, O. (2018).  Robust-to-Dynamics Optimization.  _arXiv e-prints_ , arXiv:1805.03682.
  * ApS (2017) ApS, M. (2017).  MOSEK Optimization Suite Release 8.1.0.43.  URL: http://docs.mosek.com/8.1/intro.pdf.
  * Barmish (1985) Barmish, B.R. (1985).  Necessary and sufficient conditions for quadratic stabilizability of an uncertain system.  _Journal of Optimization theory and applications_ , 46(4), 399–408.
  * Blanchini and Miani (2015) Blanchini, F. and Miani, S. (2015).  _Set-theoretic methods in control_.  Springer, second edition.
  * Blekherman et al. (2012) Blekherman, G., Parrilo, P.A., and Thomas, R.R. (2012).  _Semidefinite Optimization and Convex Algebraic Geometry_.  Society for Industrial and Applied Mathematics, Philadelphia, PA.
  * Boyd et al. (1994) Boyd, S.P., El Ghaoui, L., Feron, E., and Balakrishnan, V. (1994).  _Linear matrix inequalities in system and control theory_ , volume 15.  SIAM.
  * Hiriart-Urruty and Lemaréchal (2012) Hiriart-Urruty, J.B. and Lemaréchal, C. (2012).  _Fundamentals of convex analysis_.  Springer Science & Business Media.
  * Johansson and Rantzer (1998) Johansson, M. and Rantzer, A. (1998).  Computation of piecewise quadratic lyapunov functions for hybrid systems.  _IEEE Transactions on Automatic Control_ , 43, 555–559.
  * Legat (2020) Legat, B. (2020).  _Set programming : theory and computation_.  Ph.D. thesis, UCLouvain.
  * Legat et al. (2020a) Legat, B., Dowson, O., Garcia, J.D., and Lubin, M. (2020a).  Mathoptinterface: a data structure for mathematical optimization problems.  Second round of reviews.
  * Legat et al. (2019) Legat, B., Jungers, R.M., Parrilo, P.A., and Tabuada, P. (2019).  Set Programming with JuMP.  In _The Third Annual JuMP-dev Workshop_.
  * Legat et al. (2021) Legat, B., Raković, S.V., and Jungers, R.M. (2021).  Piecewise semi-ellipsoidal control invariant sets.  _IEEE Control Systems Letters_ , 5(3), 755–760.  10.1109/LCSYS.2020.3005326.
  * Legat et al. (2020b) Legat, B., Tabuada, P., and Jungers, R.M. (2020b).  Sum-of-squares methods for controlled invariant sets with applications to model-predictive control.  _Nonlinear Analysis: Hybrid Systems_ , 36, 100858.  10.1016/j.nahs.2020.100858.
  * Legat and Jungers (2021) Legat, B. and Jungers, R.M. (2021).  Continuous-time controlled invariant sets, a geometric approach.  https://www.codeocean.com/.
  * Liapounoff (1907) Liapounoff, A. (1907).  Problème général de la stabilité du mouvement.  _Annales de la Faculté des sciences de Toulouse : Mathématiques_ , 9, 203–474.  URL http://eudml.org/doc/72801.
  * Murty and Kabadi (1987) Murty, K.G. and Kabadi, S.N. (1987).  Some np-complete problems in quadratic and nonlinear programming.  _Mathematical Programming: Series A and B_ , 39(2), 117–129.
  * Parrilo (2000) Parrilo, P.A. (2000).  _Structured semidefinite programs and semialgebraic geometry methods in robustness and optimization_.  Ph.D. thesis, Citeseer.
  * Petersen (1985) Petersen, I. (1985).  Quadratic stabilizability of uncertain linear systems: existence of a nonlinear stabilizing control does not imply existence of a linear stabilizing control.  _IEEE Transactions on Automatic Control_ , 30(3), 291–293.
  * Prajna et al. (2004) Prajna, S., Papachristodoulou, A., and Wu, F. (2004).  Nonlinear control synthesis by sum of squares optimization: A lyapunov-based approach.  In _5th Asian Control Conference_ , volume 1, 157–165. IEEE.
  * Rockafellar (2015) Rockafellar, R.T. (2015).  _Convex analysis_.  Princeton university press.
  * Schneider (2013) Schneider, R. (2013).  _Convex bodies: the Brunn–Minkowski theory_.  151\. Cambridge University Press.
  * Sontag (1983) Sontag, E.D. (1983).  A lyapunov-like characterization of asymptotic controllability.  _SIAM Journal on Control and Optimization_ , 21(3), 462–471.
  * Toker and Ozbay (1995) Toker, O. and Ozbay, H. (1995).  On the np-hardness of solving bilinear matrix inequalities and simultaneous stabilization with static output feedback.  In _Proceedings of the American Control Conference_ , volume 4, 2525–2526. IEEE.
  * Weisser et al. (2019) Weisser, T., Legat, B., Coey, C., Kapelevich, L., and Vielma, J.P. (2019).  Polynomial and moment optimization in julia and jump.  In _JuliaCon_.  URL https://pretalx.com/juliacon2019/talk/QZBKAU/.
  * Wonham (1985) Wonham, W.M. (1985).  Linear multivariable control: A geometric approach.  In _Applications of Mathematics_ , volume 10. Springer, third edition.  10.1007/978-1-4612-1082-5.

##  Appendix A Block matrices

We have the following result for block matrices.

###### Proposition 14

Consider a symmetric matrix $A\in\mathbb{R}^{n\times n}$ and a matrix
$B\in\mathbb{R}^{n\times m}$. If the matrix

$C=\begin{bmatrix}A&B\\\ B^{\top}&0\end{bmatrix}$

is positive semidefinite then $B$ is zero.

###### Proof A.1

If $C$ is positive semidefinite, then there exists an integer $r$ and matrices
$X\in\mathbb{R}^{r\times n},Y\in\mathbb{R}^{r\times m}$ such that

$C=\begin{bmatrix}X^{\top}\\\
Y^{\top}\end{bmatrix}\begin{bmatrix}X&Y\end{bmatrix}.$

From $Y^{\top}Y=0$, we deduce that $Y=0$ hence $B=X^{\top}Y=0$.

##  Appendix B Convex geometry

######  Definition 15 (Support function (Rockafellar, 2015, p. 28))

Consider a convex set $\mathcal{C}$. The _support function_ of $\mathcal{C}$
is defined as

$\delta^{*}(y|\mathcal{C})=\sup_{x\in\mathcal{C}}\langle y,x\rangle.$

######  Definition 16 (Polar set)

For any convex set $\mathcal{C}$ the polar of $\mathcal{C}$, denoted
$\mathcal{C}^{\circ}$, is defined as

$\mathcal{C}^{\circ}=\\{\,y\mid\delta^{*}(y|\mathcal{C})\leq 1\,\\}.$

We define the _tangent cone_ as follows (Blanchini and Miani, 2015, Definition
4.6).

######  Definition 17 (Tangent cone)

Given a closed convex set $\mathcal{S}$ and a distance function
$d(\mathcal{S},x)$, the _tangent cone_ to $\mathcal{S}$ at $x$ is defined as
follows:

$T_{\mathcal{S}}(x)=\left\\{\,y\mid\lim_{\tau\to 0}\frac{d(\mathcal{S},x+\tau
y)}{\tau}=0\,\right\\}$

where the distance is defined as

$d(\mathcal{S},x)=\inf_{y\in\mathcal{S}}\|x-y\|$

where $\|\cdot\|$ is a norm. The tangent cone is a convex cone and is
independent of the norm used.

For a convex set $\mathcal{C}$, the _normal cone_ is the polar of the tangent
cone $N_{\mathcal{C}}(x)=T_{\mathcal{C}}^{\circ}(x)$.

The _exposed face_ (also called the _support set_ , e.g., in (Schneider, 2013,
Section 1.7.1)) is defined as follows (Hiriart-Urruty and Lemaréchal, 2012,
Definition 3.1.3).

######  Definition 18 (Exposed face)

Consider a nonempty closed convex set $\mathcal{C}$. Given a vector $y\neq 0$,
the _exposed face_ of $\mathcal{C}$ associated to $y$ is

$F_{\mathcal{C}}(y)=\\{\,x\in\mathcal{C}\mid\langle
x,y\rangle=\delta^{*}(y|\mathcal{C})\,\\}.$

The exposed faces and normal cones are related by the following property
(Hiriart-Urruty and Lemaréchal, 2012, Proposition C.3.1.4).

###### Proposition 19

Consider a nonempty closed convex set $\mathcal{C}$. For any $x\in\mathcal{C}$
and nonzero vector $y$, $x\in F_{\mathcal{C}}(y)$ if and only if $y\in
N_{\mathcal{C}}(x)$.

Given a set $\mathcal{S}$ and a matrix $A$, let $A^{-\top}$ denote the
preimage $\\{\,x\mid A^{\top}x\in\mathcal{S}\,\\}$.

######  Proposition 20 ((Rockafellar, 2015, Corollary 16.3.2))

For any convex set $\mathcal{C}$ and linear map $A$,

$\displaystyle(A\mathcal{C})^{\circ}$
$\displaystyle=A^{-\top}\mathcal{C}^{\circ}.$

where $\mathcal{C}^{\circ}$ denotes the polar of the set $\mathcal{C}$.

When the support function is differentiable at a given point,
$F_{\mathcal{C}}$ is a singleton and may be directly obtained using the
following result:

######  Proposition 21 ((Rockafellar, 2015, Corollary 25.1.2))

Given a nonempty closed convex set $\mathcal{C}$, if
$\delta^{*}(y|\mathcal{C})$ is differentiable at $y$ then
$F_{\mathcal{C}}(y)=\\{\nabla\delta^{*}(y|\mathcal{C})\\}$.

In fact, for nonempty compact convex sets, the differentiability at $y$ is
even a necessary and sufficient conditions for the unicity of
$F_{\mathcal{C}}(y)$ (Schneider, 2013, Corollary 1.7.3).


changing linearly” means that (a) for a fixed $i$, there is some $N\geq 0$
such that for every $(e_{i},e_{i+1})$, the set
$\textnormal{NE}^{\textnormal{ext}}_{2}({\mathcal{G}}_{\vec{\alpha}}^{\texttt{c}},S_{1},S_{2})$
has exactly $N$ elements and (b) there are linear functions
$\mathpzc{t}_{2}^{n}:(e_{i},e_{i+1})\to\Pi_{2}$, $n=1,\dots,N$, such that for
every $\texttt{c}\in(e_{i},e_{i+1})$,
$\textnormal{NE}^{\textnormal{ext}}_{2}({\mathcal{G}}_{\vec{\alpha}}^{\texttt{c}},S_{1},S_{2})=\\{\mathpzc{t}_{2}^{n}(\texttt{c})\mid
n=1,\dots,N\\}.$

$\left[\begin{array}[]{ccccccccc|c}1&1&\cdots&1&&\dots&&0&\phantom{\shortminus}0&1\\\
&&&&&&&&&\\\
u_{1}(a_{1},b_{1})&u_{1}(a_{1},b_{2})&\cdots&u_{1}(a_{1},b_{m})&&\dots&&0&\shortminus
1&\alpha_{a_{1}}\\\
u_{1}(a_{2},b_{1})&u_{1}(a_{2},b_{2})&\cdots&u_{1}(a_{2},b_{m})&&\dots&&0&\shortminus
1&\alpha_{a_{2}}\\\
\vdots&\vdots&\vdots&\vdots&&&&\vdots&\phantom{\shortminus}\vdots&\vdots\\\
u_{1}(a_{n},b_{1})&u_{1}(a_{n},b_{2})&\cdots&u_{1}(a_{n},b_{m})&&\dots&&0&\shortminus
1&\alpha_{a_{n}}\\\
u_{1}(a^{\prime}_{1},b_{1})&u_{1}(a^{\prime}_{1},b_{2})&\cdots&u_{1}(a^{\prime}_{1},b_{m})&1&0&\dots&0&\shortminus
1&\alpha_{a^{\prime}_{1}}\\\
u_{1}(a^{\prime}_{2},b_{1})&u_{1}(a^{\prime}_{2},b_{2})&\cdots&u_{1}(a^{\prime}_{2},b_{m})&0&1&&0&\shortminus
1&\alpha_{a^{\prime}_{2}}\\\
\vdots&\vdots&\ddots&\vdots&\vdots&&\ddots&\vdots&\phantom{\shortminus}\vdots&\vdots\\\
u_{1}(a^{\prime}_{n^{\prime}},b_{1})&u_{1}(a^{\prime}_{n^{\prime}},b_{2})&\cdots&u_{1}(a^{\prime}_{n^{\prime}},b_{m})&0&0&\dots&1&\shortminus
1&\alpha_{a^{\prime}_{n^{\prime}}}\\\ \end{array}\right],$ Figure 9: The
matrix form $Ax_{\texttt{c}}^{\textbf{T}}=y^{\texttt{c}}$ of the LP
(B.14)-(B.20), where the columns are indexed by
$x_{\texttt{c}}=(\pi^{\texttt{c}}_{2}(b_{1}),\dots,\pi^{\texttt{c}}_{2}(b_{1}),w_{a^{\prime}_{1}},\dots,w_{a^{\prime}_{n^{\prime}}},\gamma)$.
(The numbers $m$, $n$, and $n^{\prime}$ stand for the size of $S_{2}$,
$S_{1}$, and $\mathcal{A}_{1}\setminus S_{1}$ respectively.) The additional
constraints are $\pi^{\texttt{c}}_{2}(b_{j})\geq 0$ and $w_{a}^{\prime}\geq
0$. Note that because P1’s utilities are adjusted independently of P2’s
actions, the adjustments $\alpha_{a}$ and $\alpha_{a^{\prime}}$ can be moved
to right-hand side of the equation.

###### Proof.

Let $({\mathcal{G}},\vec{\alpha})$ be a game with linearly adjustable payoffs
and suppose that c is such that there exists some NE in
$({\mathcal{G}}_{\vec{\alpha}}^{\texttt{c}},S_{1},S_{2})$.

As the first step, we rewrite the relevant part of the LP from earlier. Using
Lemma 32 (in combination with (1) from Lemma 33), we see that a policy
$\pi_{2}$ lies in
$\textnormal{NE}^{\textnormal{ext}}_{2}({\mathcal{G}}_{\vec{\alpha}}^{\texttt{c}},S_{1},S_{2})$
if and only if it is a basic feasible solution of the following LP:

$\displaystyle\sum\nolimits_{b\in S_{2}}\pi_{2}(b)=1$ (B.7) $\displaystyle
u_{1}(a,\pi_{2})-\texttt{c}\alpha_{a}=\gamma$ $\displaystyle\textnormal{ for
}a\in S_{1}$ (B.8) $\displaystyle
u_{1}({a^{\prime}},\pi_{2})-\texttt{c}\alpha_{a^{\prime}}\geq\gamma$
$\displaystyle\textnormal{ for }{a^{\prime}}\in a_{1}\setminus S_{1}$ (B.9)
where the variables satisfy (B.10) $\displaystyle\pi_{2}(b)=0$
$\displaystyle\textnormal{ for }b\in\mathcal{A}_{2}\setminus S_{2}$ (B.11)
$\displaystyle\pi_{2}(b)\geq 0$ $\displaystyle\textnormal{ for }a\in S_{2}$
(B.12) $\displaystyle\gamma\in\mathbb{R}.$ (B.13)

We turn all of the inequalities into equalities by introducing slack variables
$w_{a^{\prime}}\geq 0$, ${a^{\prime}}\in\mathcal{A}_{1}\setminus S_{1}$:

$\displaystyle\sum\nolimits_{b\in S_{2}}\pi_{2}(b)=1$ (B.14) $\displaystyle
u_{1}(a,\pi_{2})-\texttt{c}\alpha_{a}=\gamma$ $\displaystyle\textnormal{ for
}a\in S_{1}$ (B.15) $\displaystyle
u_{1}({a^{\prime}},\pi_{2})-\texttt{c}\alpha_{a}+w_{a^{\prime}}=\gamma$
$\displaystyle\textnormal{ for }{a^{\prime}}\in\mathcal{A}_{1}\setminus S_{1}$
(B.16) where the variables satisfy (B.17) $\displaystyle\pi_{2}(b)=0$
$\displaystyle\textnormal{ for }b\in\mathcal{A}_{2}\setminus S_{2}$ (B.18)
$\displaystyle\pi_{2}(b)\geq 0$ $\displaystyle\textnormal{ for }a\in S_{2}$
(B.19) $\displaystyle\gamma\in\mathbb{R}.$ (B.20)

Second, we rewrite the LP (B.14)-(B.20) in a matrix form. We denote the
relevant variables as
$x_{\texttt{c}}=(\pi^{\texttt{c}}_{2}(b_{1}),\dots,\pi^{\texttt{c}}_{2}(b_{n}),w_{{a^{\prime}}_{1}},\dots,w_{{a^{\prime}}_{n^{\prime}}},\gamma)$,
where $n\ :=\lvert\mathcal{A}_{1}\rvert$, $n^{\prime}\
:=\lvert\mathcal{A}_{1}\setminus S_{1}\rvert$. In this notation, there will be
some matrix $A$ and right-hand side $y^{\texttt{c}}=(1,0,\dots,0)$ for which
some $\pi_{2}$ is a solution of (B.14)-(B.20) if and only if it satisfies
$Ax_{\texttt{c}}^{\textbf{T}}=y^{\texttt{c}}$ and
$\pi^{\texttt{c}}_{2}(b_{j}),w_{a^{\prime}}\geq 0$. However, we can
additionally use the fact that the adjustment cost $\texttt{c}\alpha_{a}$ that
P1 pays does not depend on the action of P2. This allows us the matrix form
depicted in Figure 9, where all the $\texttt{c}\alpha_{a}$-s have been moved
to the right-hand side.

As the third step, we note that while the matrix from Figure 9 might have
linearly dependent rows, we can always replace it by a matrix whose rows are
linearly independent. To see this, note first that clearly no row
corresponding to one of the actions ${a^{\prime}}\in\mathcal{A}_{1}\setminus
S_{1}$ can be expressed as a linear combination of any other rows, because of
the $1$-s in the bottom-right corner of the matrix. Second, it is possible
that one of the rows corresponding to some $a\in S_{1}$ can be expressed as a
linear combination of the other rows corresponding to $S_{1}$. Suppose that
$a$ is such action and $\lambda_{i}\in\mathbb{R}$ are the corresponding
weights. There are two options: If
$\sum_{i}\lambda_{i}\alpha_{a_{i}}=\alpha_{a}$, the condition corresponding to
$a$ can be omitted, since it is already subsumed by the conditions
corresponding to $S_{1}\setminus\\{a\\}$. Conversely, if
$\sum_{i}\lambda_{i}\alpha_{a_{i}}\neq\alpha_{a}$, the system of equations
from Figure 9 will be unsolvable. However, this case is ruled out by our
assumption that c is such that $\textnormal{NE}({\mathcal{G}},S_{1},S_{2})$ is
non-empty. In summary: For the remainder of the proof, we can assume that the
rows of the matrix $A$ are linearly independent.

Fourth, we identify the basic solution of the LP given by Figure 9. For the
purpose of this step, denote the set of column-indices of $A$ as
$\mathcal{I}:=S_{2}\cup\left(\mathcal{A}_{1}\setminus
S_{1}\right)\cup\\{\gamma\\}$. For a “basis” $B\subseteq\mathcal{I}$, we use
$A_{B}$ to denote the sub-matrix of $A$ consisting of the columns indexed by
$B$. By $\mathcal{B}$, we denote the set of all $B$-s for which the sub-matrix
$A_{B}$ is regular. Finally, for $B\in\mathcal{B}$, we denote by
$x_{\texttt{c}}^{B}$ the basic solution corresponding to $B$ — i.e., the
solution of $Ax_{\texttt{c}}^{\textbf{T}}=y^{\texttt{c}}$ for which all the
variables indexed by $\mathcal{I}\setminus B$ are equal to $0$. By definition
of a BFS, the basic feasible solutions of the LP are precisely all the vectors
of the form $x_{\texttt{c}}^{B}$, $B\in\mathcal{B}$.

Fifth, we show that every basic (not necessarily feasible solution of the LP
given by Figure 9 changes linearly with c. To see this, note that each basic
(not necessarily feasible) solution $x_{\texttt{c}}^{B}$ can be written as the
vector $A_{B}^{-1}y^{\texttt{c}}$, extended by $0$-s at the indices
$\mathcal{I}\setminus B$. (Since $A_{B}$ is assumed to be regular, the inverse
exists.) Since the matrix $A$ does not depend on c and $y^{\texttt{c}}$ only
depends on c linearly, the mapping $\texttt{c}\mapsto x_{\texttt{c}}^{B}$ is
linear.

Finally, we conclude the proof. To do this, recall that a basic solution
$x_{\texttt{c}}^{B}$ is feasible if all of the variables
$\pi^{\texttt{c}}_{2}(b),w_{a}$ are non-negative. For every
$x_{\texttt{c}}^{B}$, the set of the values of c for which all of these
definitions are satisfied is going to be some (possibly empty or trivial)
closed interval $[e^{B}_{0},e^{B}_{1}]$. By taking the set $\\{e^{B}_{i}\mid
i=0,1,\,B\in\mathcal{B}\\}\cup\\{-\infty,\infty\\}$ and reordering it as an
increasing sequence, we obtain the desired breakpoint set $E$. This completes
the whole proof ∎

###### Proof of Lemma 30.

Let $({\mathcal{G}},\vec{\alpha})$ be a game with linearly adjustable payoffs.
To get the desired sequence of breakpoints, we take – for every pair
$S_{1}\subseteq\mathcal{A}_{1}$ and $S_{2}\subseteq\mathcal{A}_{2}$ – some set
$E(S_{1},S_{2})$ of breakpoints given by Proposition 35 and define
$E:=\bigcup\nolimits_{S_{1}\subseteq\mathcal{A}_{1},S_{2}\subseteq\mathcal{A}_{2}}E(S_{1},S_{2})$.
We then enumerate $E$ as a strictly increasing sequence $(e_{i})_{i}$ To prove
our result, let $\texttt{c}\in[e_{i},e_{i+1}]$, and
$\pi^{\texttt{c}_{0}}\in\textnormal{NE}({\mathcal{G}}_{\vec{\alpha}}^{\texttt{c}_{0}})$.
We finding a trajectory $\mathpzc{t}_{2}$ and that satisfies the conclusion of
the lemma.

By 31, there are some sets $S_{1}$, $S_{2}$ for which
$\pi^{\texttt{c}_{0}}\in\textnormal{NE}({\mathcal{G}}_{\vec{\alpha}}^{\texttt{c}_{0}},S_{1},S_{2})$.
By Lemma 33 and the subsequent observation, there are some basic feasible
solutions $\beta^{1},\dots,\beta^{N}$ (of the LP from Lemma 32) and convex
combination $\lambda_{1},\dots,\lambda_{N}$ such that
$\pi^{\texttt{c}_{0}}_{2}=\sum\nolimits_{i=1}^{N}\lambda_{i}\beta^{i}_{2}$. By
Proposition 35, there are some linear trajectories
$\mathpzc{t}^{1}_{2},\dots,\mathpzc{t}^{N}_{2}$ such that
$\mathpzc{t}_{2}^{i}(\texttt{c}_{0})=\beta^{i}_{2}$ and for every
$\texttt{c}\in(e_{i},e_{i+1})$,
$\mathpzc{t}_{2}^{i}(\texttt{c})\in\textnormal{NE}^{\textnormal{ext}}_{2}({\mathcal{G}}_{\vec{\alpha}}^{\texttt{c}},S_{1},S_{2})$.
Moreover, by Lemma 34, we have
$(\pi^{\texttt{c}_{0}}_{1},\mathpzc{t}_{2}^{i}(\texttt{c}))\in\textnormal{NE}({\mathcal{G}}_{\vec{\alpha}}^{\texttt{c}},S_{1},S_{2})$
for every $\texttt{c}\in(e_{i},e_{i+1})$. By the convexity of NE (Lemma 33),
$\mathpzc{t}_{2}:=\sum\nolimits_{i=1}^{N}\mathpzc{t}_{2}^{i}$ is a linear
trajectory for which $\mathpzc{t}_{2}(\texttt{c}_{0})=\pi^{\texttt{c}_{0}}$
and for every $\texttt{c}\in(e_{i},e_{i+1})$,
$(\pi^{\texttt{c}_{0}},\mathpzc{t}_{2}^{i}(\texttt{c}))\in\textnormal{NE}({\mathcal{G}}_{\vec{\alpha}}^{\texttt{c}},S_{1},S_{2})$.
Moreover, since the utilities in ${\mathcal{G}}_{\vec{\alpha}}^{\texttt{c}}$
depend continuously on c, this also implies that
$(\pi^{\texttt{c}_{0}},\mathpzc{t}_{2}^{i}(\texttt{c}))\in\textnormal{NE}({\mathcal{G}}_{\vec{\alpha}}^{\texttt{c}},S_{1},S_{2})$
for $\texttt{c}\in\\{e_{i},e_{i+1}\\}$. (By 31,) this concludes the whole
proof. ∎


# Advancing Set-Conditional Set Generation:
Diffusion Models for Fast Simulation of Reconstructed Particles

Dmitrii Kobylianskii<EMAIL_ADDRESS>Weizmann Institute of
Science, Israel  These authors contributed equally  Nathalie Soybelman
<EMAIL_ADDRESS>Weizmann Institute of Science, Israel  These
authors contributed equally  Nilotpal Kakati  Weizmann Institute of Science,
Israel  Etienne Dreyer  Weizmann Institute of Science, Israel  Benjamin
Nachman  Lawrence Berkeley National Laboratory, Berkeley, USA  Eilam Gross
Weizmann Institute of Science, Israel

###### Abstract

The computational intensity of detector simulation and event reconstruction
poses a significant difficulty for data analysis in collider experiments. This
challenge inspires the continued development of machine learning techniques to
serve as efficient surrogate models. We propose a fast emulation approach that
combines simulation and reconstruction. In other words, a neural network
generates a set of reconstructed objects conditioned on input particle sets.
To make this possible, we advance set-conditional set generation with
diffusion models. Using a realistic, generic, and public detector simulation
and reconstruction package (COCOA), we show how diffusion models can
accurately model the complex spectrum of reconstructed particles inside jets.

##  I Introduction

For experiments at the Large Hadron Collider (LHC), the demand for simulated
data has increased throughout recent data-taking periods and will surge by
roughly a factor of ten during the High-Luminosity runs [1]. The standard
simulation pipeline used at the general-purpose experiments ATLAS and CMS is
designed to closely mirror that of real data [2, 3]. Initially, proton-proton
collision events are generated, followed by parton shower, hadronization, and
secondary decays. The set of stable, or “truth”, particles then enter the
detector volume, and their material interactions are modeled on a microscopic
level using Geant4 [4, 5, 6]. The resulting energy deposits in the tracker and
calorimeter cells seed the formation of tracks and clusters, respectively.
Subsequently, a reconstruction algorithm refines these data into objects
suitable for statistical analysis. Depending on the algorithm, the
reconstructed “particle flow objects” (PFOs) either represent particles
directly or else consist of tracks and calorimeter clusters with kinematics
adjusted to maintain energy flow at the jet level.

From the perspective of computational efficiency, the main constraint in the
traditional simulation pipeline comes from the highly sophisticated, but
expensive model of particle-detector interactions provided by Geant4. In a bid
to remove this constraint, various machine learning approaches are currently
being explored [7, 8]. These rapid surrogate models aim to drastically reduce
computational time while otherwise performing similarly to the traditional
simulation pipeline.

The majority of research efforts in this direction have concentrated on fast
calorimeter simulation [7, 8, 9, 10]. In this approach, the objective is to
replace the slowest link in the simulation chain while preserving the other
steps, including reconstruction. Utilizing the features of an incoming
particle as a conditioning factor, a generative model is employed to produce
the detector response, which can then be aggregated for all particles within
an event [11]. Within the scope of the Fast Calorimeter Challenge [12, 13,
14], a wide array of methods have been developed and analyzed. Various
architectures, including those based on Generative Adversarial Networks (GANs)
[15, 16], normalizing flows [17, 18] and diffusion models [19, 20, 21, 22, 23]
have been introduced. These approaches represent the calorimeter’s response
through diverse formats such as images, point clouds, or graphs.

A complementary approach is to replace the entire simulation pipeline by
combining event generation, detector simulation, and event reconstruction.
While reconstruction has typically not been a computational bottleneck, this
is changing with improvements in fast detector simulations and with increased
combinatorics from higher instantaneous luminosities [24, 25]. While we are
not aware of a model demonstrated for entire events at the level of individual
PFOs, there are a number of proposals using high-level objects [26]. This
approach is highly efficient as it encompasses all simulation steps in a
single process, eliminating the need for additional processing. However, it is
worth noting that this method is process-dependent, necessitating optimization
and retraining of the network for each analysis. A related strategy entails
creating a generative model for individual jets, the most complex part of
entire events [27, 28, 29, 30, 31, 32, 33, 34]. In this scenario, the
constituents of the jet are generated conditioned on both the jet features and
the particle type. This approach is less process-dependent than learning
entire events, but it still requires the model to learn about the physics of
jet substructure and when moving to entire events, it necesitates event-
splitting and additional processing for other objects within the event.

Our goal is to develop a fast and automated end-to-end strategy that starts at
truth particles and ends with reconstructed particles. By only learning the
detector response and not the full event generation, we strive to be process
independent. Going all the way to PFOs addresses the computational bottleneck
from reconstruction and keeps the dimensionality of the problem manageable. In
the pursuit of a process-independent, full-event method, a novel strategy for
fast simulation known as FlashSim [35, 36] has been introduced111The model
presented in Ref. [37] has similar goals, but there is no matching between the
truth and reconstructed objects during training.. This method involves
utilizing stable particles to predict high-level objects such as jets, fat
jets, muons, electrons, and more. For each object, a separate network is
trained using distinct input and target variables. Extending the end-to-end
methodology to indiviudal truth and reconstructed particles requires going
beyond all examples mentioned above because we need to generate a variable-
length set of particles conditioned on another variable-length set of
particles. A prototype for this approach was presented in Ref. [38] using a
simplified version of the problem. This initial model solely considered
charged particles within a jet, with smeared tracks as targets. The
architecture incorporated graph neural networks with slot attention
mechanisms. Now, transitioning towards a more realistic approach that
encompasses both neutral and charged particles while targeting reconstructed
particles, we explore diffusion models, showcasing enhanced performance
compared to the baseline slot attention model. These models are studied using
a detailed and realistic detector simulation and reconstruction. We will show
that continuous normalizing flows [39] are the most effective approach of the
methods we studied. This forms the bassis of a companion paper [40] that shows
how an end-to-end fast simulation of the CMS experiment can outperform
classical baselines.

This paper is organized as follows. Section II introduces the dataset, with
both the detailed simulation and reconstruction steps. The neural network
approaches are described in Sec. III and trained in Sec. IV. Numerical results
are presented in Sec. V and the paper ends with conclusions and outlook in
Sec. VI.

##  II Dataset

The network takes as input a single jet of truth particles entering the
detector222We consider single jets here as the most complex objects; future
work will explore entire events.. For the generation task, the target is the
set of reconstructed particles. Instead of relying on a parametrized smearing
model (e.g. Delphes [41]), we obtain reconstructed particles using a realistic
detector simulation followed by a particle flow algorithm described below.
Each input or output reconstructed particle is represented by its momentum,
direction, and charge ($p_{\text{T}}$, $\eta$, $\phi$, $|q|$). Notably, charge
prediction is a new feature introduced in this study, as our previous work
utilized a toy model focusing solely on charged particles. Generally,
reconstructed particles are classified into five classes: charged hadrons,
neutral hadrons, photons, electrons, and muons. Given the significant class
imbalance within the dataset, we simplify the task by predicting only whether
the object is neutral or charged.

###  II.1 Truth Event Generation

We concentrate on a localized reconstruction of particles within a single jet.
For this purpose, we utilize Pythia8.3 [42] to generate a single quark with
momentum ranging between 10 GeV and 200 GeV, with the initial direction
randomly selected within the ranges $|\eta|<2.5$ and $-\pi<\phi<\pi$.
Subsequent to the parton shower and hadronization, only stable particles with
momentum exceeding 1 GeV and $|\eta|<3.0$ are retained. The set of particles
contained within the jets exhibits an average cardinality of $N=6.7$, with a
maximum of 25 particles and a minimum of 1 particle per truth event, as
measured on the whole dataset.

###  II.2 Detector simulation

The detector simulation uses COCOA [43], a Geant4-based configurable
calorimeter simulation toolkit. The detector comprises three electromagnetic
and three hadronic calorimeter layers, simulating ATLAS-like materials.
Geometrically, the coverage is divided into a barrel region ($|\eta|<1.5$) and
two endcaps ($1.5<|\eta|<3.0$). The cell depth follows a $1/\cosh{\eta}$
profile to maintain a constant effective interaction depth across $\eta$.
Moreover, the tracking region of the detector is subjected to a magnetic field
of 3.8 T. To avoid particles created upstream of the calorimeter, our detector
simulation is simplified by assuming that the tracker contains no material.

Tracking effects are emulated by first computing the trajectories of charged
particles within the magnetic field, followed by applying smearing to the
track parameters $q/p$, $\theta$, and $\phi$. Additionally, we discard tracks
originating far from the beamline, with a transverse radius exceeding $75$ mm
($250$ mm) in the barrel (endcap) region.

Following the simulation, a topological clustering algorithm is employed to
group cells based on their deposited energy, expected noise levels, and
proximity. These calorimeter clusters, along with the cells belonging to them
and the set of tracks, are the input for the reconstruction algorithm.

###  II.3 Particle Reconstruction

To reconstruct particles, we utilize a recently developed machine learning
algorithm called HGPflow [44]. The design principle of HGPflow is to learn the
energy assignment between the input set of tracks and calorimeter clusters and
an output set of particles using a hypergraph prediction network [45]. We
trained a version of HGPflow on a statistically independent dataset of 294k
single jets and evaluated this model to obtain a set of reconstructed
particles for each jet in our training dataset. These predictions serve as the
target for our fast simulation models.

###  II.4 Preprocessing

To mitigate contamination from detector noise potentially converted into
neutral particles during HGPflow reconstruction, we employ a filtering
process. Neutral reconstructed particles located outside a cone of $R=0.4$
around the truth-jet axis are removed since these typically correspond to
clusters of calorimeter noise. Subsequently, before inputting the object
features into the networks, we apply relative scaling. For each event
individually, we scale $\log p_{T}$, $\eta$, and $\phi$ based on the mean and
standard deviations of the truth particles. The resulting scaling quantities
of $\eta$ and $\log p_{T}$ are then added as global features to retain
knowledge of the absolute jet properties. Additionally, we remove events where
truth particles have a $|\phi|>2.8$ to avoid edge effects.

##  III Fast simulation algorithms

In this study, we evaluate the performance of two set-conditional set
generation approaches. Both methods rely on graph neural networks, as they can
naturally represent unordered, set-valued data through a graph structure.
Currently, we use fully connected graphs because we have single jets. This can
be readily generalized when considering full events where different event
regions are loosely or not at all connected. Edges within the graph allow
inter-particle relations to be encoded through message passing. Simplified
illustrations of the network architectures are depicted in Fig. 1.

Figure 1: Simplified sketches of the two proposed network architectures for
the set generation task. EDM and CFM extend the Slot Attention approach
presented in detail in [38] with an iterative denoising procedure.

###  III.1 Slot Attention (SA)

The slot-attention approach, initially proposed in Ref. [38], constitutes the
baseline model. In this approach, the generation task is divided into two
stages: cardinality and feature prediction. Cardinality is inferred from the
updated representations of truth particles following message passing.
Subsequently, based on the predicted cardinality, particles are initialized
with random noise and iteratively refined through a slot-attention block [46],
which attends to the updated truth particle representations. The introduction
of noise during initialization is crucial to prevent deterministic
predictions. A WGAN [47]-like setup is used to train the model whereby a
distribution-level distance is used to steer the generator. In our previous
and current work, we recognize(d) that the learned noise model exhibits room
for improvement, indicated by the limited precision of the predictions.

We explore diffusion models in order to address these limitations. Diffusion
models are state-of-the-art in generative AI and aim to directly learn
explicitly or implicity the probability density (or the gradient of the
density). Details of the EDM Diffusion (EDM) and Continuous Flow Matching
(CFM) approaches are detailed in the subsequent section. The CFM is based on
normalizing flows, but is still a variation on diffusion models; we will refer
to both EDM and CFM as ‘diffusion’ while EDM will refer specifically to the
score-based approach described below.

###  III.2 Diffusion Models

In this approach, we maintain the cardinality prediction method as originally
introduced in the slot-attention approach and solely modify the feature
prediction task. After noting that a single slot-attention block lacks the
expressiveness needed to transform complete noise into sufficiently accurate
features, we opt for a different approach. This involves replacing the feature
prediction component with a diffusion process or flow matching. Starting with
a fully noised set of particles, we use a transformer-based network to
iteratively predict a feature estimate in the case of diffusion or a vector
field in the case of the flow matching objective. This output is used to
update the noised set until we achieve the final, denoised reconstructed
particles.

Network and preconditioning
---
Skip scaling $c_{\text{skip}}(\sigma)$ | $\sigma_{\text{data}}^{2}/(\sigma^{2}+\sigma_{\text{data}}^{2})$
Output scaling $c_{\text{out}}(\sigma)$ | $\sigma\cdot\sigma_{\text{data}}/\sqrt{\sigma^{2}+\sigma_{\text{data}}^{2}}$
Input scaling $c_{\text{in}}(\sigma)$ | $1/\sqrt{\sigma^{2}+\sigma_{\text{data}}^{2}}$
Training
Loss weighting $\lambda(\sigma)$ | ($\sigma^{2}+\sigma_{\text{data}}^{2})/(\sigma\cdot\sigma_{\text{data}})^{2}$
Sampling
ODE solver | 4-th order PNDM [48] method
Time steps | $(\sigma_{\text{max}}^{1/\rho}+\frac{i}{N-1}(\sigma_{\text{min}}^{1/\rho}-\sigma_{\text{max}}^{{1/\rho}}))^{\rho}$
Parameters |  $\sigma_{\text{min}}=0.002$, $\sigma_{\text{max}}=80$, $\sigma_{\text{data}}=1.1$, $\rho=7$
Table 1: Parameters of the EDM diffusion model formulation.

####  III.2.1 EDM Formulation

We adopt score-based diffusion, a method where a diffusion process gradually
perturbs the original data $x$ with Gaussian noise while the neural network
learns the time-dependent score function $\nabla_{x}\log p(x;\sigma)$. This
function enables us to reverse the noise process, starting from pure noise
$x_{T}\sim\mathcal{N}(0,1)$ and iteratively denoising it to sample from the
original data distribution $x_{0}\sim p_{\text{data}}$.

We adapt the EDM strategy introduced in Ref. [49] for score-based denoising
diffusion models. During training, we sample noise rates $\sigma$ from the
log-normal distribution, defined by $\log(\sigma)\sim\mathcal{N}(-0.8,0.8)$.
The network receives scaled noisy data
$x=c_{in}(\sigma)(x_{0}+\sigma\epsilon)$ as input, where
$\epsilon\sim\mathcal{N}(0,1)$. To enhance network predictions, we combine the
neural network output $F_{\theta}(x;\frac{1}{4}\ln{\sigma},c)$ with a skip
connection
$D_{\theta}(x;\sigma,c)=c_{skip}(\sigma)x+c_{out}(\sigma)F_{\theta}(x;\frac{1}{4}\ln{\sigma},c)$,
where $c$ denotes contextual information. The weighted loss function minimized
during training is then defined as:

$\mathcal{L}=E_{\sigma,x_{0},\epsilon,c}\big{[}\lambda(\sigma)\left\|D_{\theta}(x;\sigma,c)-x_{0}\right\|^{2}\big{]}\,.$
(1)

We adopted specific design choices from Ref. [49], as detailed in Tab. 1.
During the inference process, we utilize the network predictions to compute
the score function $\nabla_{x}\log
p(x;\sigma,c)=(D_{\theta}(x;\sigma,c)-x)/\sigma$ and solve the diffusion
ordinary differential equation (ODE) using PNDM method from Liu et al. [48].

####  III.2.2 Flow Matching Formulation

Continuous normalizing flows (CNF) [39] are a type of generative model that
constructs a transformation from samples $x_{0}$ of a base distribution
$p_{0}$ to samples $x_{1}$ of a target distribution $p_{1}$ using an ordinary
differential equation (ODE):

$\mathrm{d}x=v_{\theta}(t,x)\mathrm{d}t,$ (2)

where the vector field $v_{\theta}$ is represented by a neural network and $t$
ranges from 0 to 1. Lipman et al. [50] recently proposed a Flow Matching (FM)
objective, which allows the vector field to be learned through a regression
task:

$\mathcal{L}_{\text{FM}}(\theta)=\mathbb{E}_{t,p_{t}(x)}||v_{\theta}(t,x)-u_{t}(x)||^{2}\,.$
(3)

The function $u_{t}$ defines the desired probability path between $p_{0}$ and
$p_{1}$, though it is generally not known in a closed form. This issue can be
addressed with the Conditional Flow Matching (CFM) objective [50, 51, 52],
optimizing which is equivalent to optimizing the FM objective. In this
formulation, $u_{t}$ and the probability path $p_{t}$ are built in a sample-
conditional manner:

$\mathcal{L}_{\text{CFM}}(\theta)=\mathbb{E}_{t,q(z),p_{t}(x|z)}||v_{\theta}(t,x)-u_{t}(x|z)||^{2},$
(4)

where $t\in[0,1]$, $z\sim q(z)$, and $x\sim p_{t}(x|z)$. Our approach uses the
modified CFM framework from [50]. By identifying the condition $z$ with the
single sample $x_{1}$ (in our context, a set of PFOs), we define:

$\displaystyle p_{t}(x|z)$
$\displaystyle=\mathcal{N}(x|tx_{1},(t\sigma-t+1)^{2})\,,$ (5) $\displaystyle
u_{t}(x|z)$ $\displaystyle=\frac{x_{1}-(1-\sigma)x}{1-(1-\sigma)t}\,,$

which outlines a linear trajectory between a standard normal distribution and
a Gaussian distribution centred at $x_{1}$ with a standard deviation of
$\sigma$, chosen as $10^{-4}$.

Following the strategy in Ref. [53], we sample $t$ from a power-law
distribution $p(t)\propto t^{\frac{1}{1+\alpha}},\ t\in[0,1]$. Setting
$\alpha=0$ results in a uniform distribution, while $\alpha>0$ gives more
weight to probability paths corresponding to larger $t$ values. Our
experiments showed that $\alpha=2$ achieves the best results.

####  III.2.3 Architecture Description

The EDM and CFM models share the same transformer-based architecture; network
hyperparameters are shown in Table 2. To improve sampling quality, we
incorporate self-conditioning on the network output from the previous
timestep, following a technique proposed in Ref. [54]. This involves
concatenating $x_{t}$ with the previously estimated $\tilde{x}/\tilde{v}$. For
training the network in this context, we set $\tilde{x}/\tilde{v}=0$ with a
probability of $p=0.5$ reverting to the model without self-conditioning. The
remaining time, first estimate $\tilde{x}/\tilde{v}$ and then utilize it for
the self-conditioning with detached gradients.

Figure 2: Schematic representation of training and generation processes of the
EDM and CFM models. For the training, noise is applied to the particle flow
objects. Using the updated truth information from the first network, the
denoised properties are obtained through a cross-attention block. For
generation, the network iteratively denoises the objects using the ODE solver.

  1. a)

Input set Encoding and Cardinality Prediction
This segment is similar to the slot attention approach, as depicted at the top
of Fig. 1. The set of truth particles characterized by $p_{T},\eta,\phi$, and
$|q|$ is enriched by a sinusoidal positional encoding, embedded using a
multilayer-perceptron (MLP) and subsequently updated through two transformer
encoder blocks. An overall embedding $G(T)$ of the truth event $T$ is
aggregated using the node-level representation. It is used to predict the
categorical distribution of the output set cardinality via an MLP.

  2. b)

Feature prediction
This network component takes as input the set of PFOs with noised features
$x_{t}$, combined with self-conditioning features $\tilde{x}/\tilde{v}$ and
the current timestep. These features are concatenated with the sinusoidal
positional encoding and embedded using an MLP. The embedding process also
incorporates the global representation of the truth event $G(T)$, event
scaling information, and the sinusoidal embedding of the timestep as
contextual information. Then, updated features are processed through three
multi-head cross-attention blocks. Each block employs two cross-attention
mechanisms to first update the PFOs using the truth particles and then update
the truth particles using the PFOs. The final representation is then processed
by an MLP together with the contextual information to produce $F_{\theta}$ or
$v_{\theta}$. In the diffusion model, the output is combined with a skip
connection to generate the denoised estimate $D_{\theta}$. A schematic
representation of this process is shown in Fig. 2.

  3. c)

Inference
Since our goal is to predict the cardinality of the generated set of
particles, we split the inference process into two steps. Initially, we
predict the cardinality of the generated particles from the set of truth
particles and initialize them with Gaussian noise. Subsequently, we employ the
PNDM method [48] in conjunction with our network to iteratively sample the
reconstructed particle features from the noise. Following generation, objects
with a $\Delta R>0.5$ with respect to the truth-jet are removed to stabilize
against outliers.

Hyperparameters
---
Batch size | 100
Optimizer | AdamW [55]
Weight decay | 0.01
Learning rate | $10^{-4}$
# of epochs | 150
# of time steps | 25
Gradient norm clip | 1.0
Max output particles | 51
# of transformer blocks | 2
# of cross attention blocks | 3
Trainable parameters | 1,742,519
Table 2: Network hyperparameters of the diffusion models.

##  IV Training and Loss

For the Slot Attention model, we adopt a familiar loss function outlined in
Ref. [38]. The central element involves the application of the Hungarian
algorithm [56] for particle matching. For each pair of input and output
particles from the two sets, we calculate the squared differences for the
particle features $\log p_{\text{T}}$ and $\eta$. Due to the periodicity of
$\phi$, we use the cosine loss $2(1-\cos\Delta\phi)$, which recovers the
squared difference when Taylor expanded to the second order. Additionally, we
include the Binary Cross Entropy (BCE) loss for the charge.

Hungarian matching determines the pairing configuration that minimizes the sum
of all pair losses per event. The total loss comprises this sum and the Cross-
Entropy (CE) loss for the cardinality prediction.

Originally, the set approach was chosen to ensure permutation invariance, as
no enforced ordering has a physical motivation. However, for the diffusion, we
found that due to the sampling, better performance is achieved if ordering is
applied so that every particle stays on its own trajectory. Consequently,
there is no need for Hungarian matching, and we use regular Mean Squared Error
(MSE) loss for the diffusion models. We also weigh the EDM model’s losses as
defined in Eq. 1.

The training, validation, and test datasets comprise 1,000,000, 20,000, and
300,000 truth events, respectively. The trainings were performed on an NVIDIA
RTX A5000.

Figure 3: Cardinality distributions for the generated PFOs and target sets for
different input set cardinalities.

##  V Results

Given the task of set generation, we can assess multiple aspects of network
performance. First, we analyze properties of the set as a whole, such as set
cardinality and inclusive particle distributions. Second, we examine the
properties of individual set constituents. Additionally, for our application,
we aim to study resolution modeling, which entails examining the variation in
generated particle properties given the same input values. These aspects are
discussed separately in the following subsections.

###  V.1 Set-based Performance

The cardinality prediction network follows the same principles and structure
for all used approaches; hence, similar results are expected. Minor deviations
may occur as the networks are trained simultaneously with the property
prediction network using combined losses. Figure 3 compares the generated and
target set cardinalities for different input set cardinalities. The
distributions for low cardinalities peak at the truth cardinality, with the
width of the distribution increasing as the number of truth particles rises.
Higher particle counts in the event lead to more instances of
misreconstructions and inefficiencies. For very high cardinalities, the peak
position is lower than the truth cardinality since dense environments are more
prone to reconstructing several neutral particles as one. Both networks
closely match the target distributions, even for high cardinalities where the
event count is low, with predictions aligning with targets within statistical
uncertainties. Next, we examine the marginal distributions of the set
constituents depicted in Fig. 4. These histograms encompass all particles from
all events in the test set. The diffusion models outperform the baseline model
in the $p_{T}$ distributions for low energies. Both models perform well in the
intermediate and high $p_{T}$ range. The Slot Attention model slightly
outperforms the diffusion models in the distribution of $\eta$. The bias comes
from the fact that low $p_{T}$ particles tend to have higher $|\eta|$ and vice
versa. However, the mismatch results in a deviation of a few percent, which is
not deemed very significant. The distributions of $\phi$ closely match the
target distribution within statistical errors for both networks.

Figure 4: Generated and target distributions of $p_{T}$, $\eta$ and $\phi$ for
all particles in the test set. Figure 5: Difference of truth and reconstructed
cardinality per event split by neutral and charged constituents.

The feature prediction task involves classifying particles as charged or
neutral. Given the initially predicted total cardinality, the network must
learn the number of charged and neutral constituents in the set. Since the
total cardinality prediction displays good agreement, we can separately
compare the neutral and charged contributions to analyze class prediction
performance. Figure 5 illustrates the difference between the number of truth
and reconstructed particles per event, separated for neutral and charged
constituents. The target distribution for charged particles is a delta
function since track presence ensures perfect classification. However, the
reconstruction of neutral particles is more prone to inefficiencies and fakes,
resulting in a wider distribution. Slot Attention demonstrates excellent
agreement with the target, while the diffusion models exhibit some
misclassification, as achieving binary prediction accuracy is challenging in
diffusion models.

Overall, we conclude that both models are able to reproduce the salient
properteis of the sets of PFOs.

Figure 6: Residual distributions of $p_{T}$, $\eta$ and $\phi$ for particles
matched between ground truth particles and target/fast simulated reconstructed
particles.

###  V.2 Constituent-based Performance

A more nuanced investigation involves a per-event comparison of set
constituents. While Fig. 4 illustrates that on average, the generated
particles have correct features, mismatches among certain features may occur
within an event, which cannot be captured in this plot. To explore this
further, we conduct Hungarian matching using the same configuration as during
training. However, here, we match the truth set with the output set. With the
matched pairs, we examine their differences in $p_{T}$, $\eta$, and $\phi$.

These residual distributions are depicted in Fig. 6. For angular quantities,
the diffusion models perform better than Slot Attention in matching the peak
of the distribution but exhibits wider tails originating from outliers,
particularly in $\eta$. Conversely, Slot Attention displays a narrower
distribution in $\eta$, suggesting that generated particles are not
sufficiently smeared with respect to the truth. This trend is also evident in
the right tail of the $p_{T}$ distribution. Here, the diffusion models
demonstrate better agreement. However, EDM deviates from the target
distribution in the peak area similarly to the SA, while CFM has a very good
match. Notably, this area is primarily influenced by low-$p_{T}$ particles,
which are challenging to model.

Figure 7: Residual distributions of $p_{T}$ for matched particles split for
different truth $p_{T}$.  Figure 8: Per event Hungarian cost for truth
matching normalized by number of matched particles in the event.

To gain deeper insights, we depict the residual $p_{T}$ distribution for
different bins of the matched truth particle $p_{T}$ in Fig. 7. Initially
focusing on the target distributions, we observe that low-energy particles
form a narrow distribution with sharply declining tails comprised of
mismatched or noisy particles. As the truth $p_{T}$ increases, the
distributions widen, likely driven by worsening track momentum resolution. All
tend to overestimate the residuals for low energies. For higher energies, we
observe excellent agreement for the diffusion models, while slot attention
predicts distributions that are too narrow, underestimating smearing effects
and predicting values too close to the truth.

We can also assess constituent-based performance at the event level. Similar
to the procedure in the training loss, we can extract the Hungarian cost for
each event from the truth and reconstructed particle features, which can
conveniently be used as a summary metric. It is depicted in Fig. 8.
Intuitively, this plot can be understood as follows: If the generated
distribution is shifted to the left of the target, i.e., the Hungarian cost is
too low, the predictions are too close to the truth and smearing effects are
not well modeled. Conversely, a shift to the right indicates excessive
smearing of truth quantities and possibly outliers. The double-peak structure
originates from events with and without high $p_{T}$ particles.

We can summarize previously made observations: Slot Attention exhibits a clear
left-shift in the bulk of the data, indicating overly precise predictions.
Both networks fail to match the left tail of the distribution, which consists
of events with very low energies. In the right peak, we observe a slight shift
to the right for the diffusion models. While the diffusion models shows room
for improvement on the edges of the distributions, they match the bulk of the
data very well and clearly outperform the baseline model. The CFM seems
somewhat superior to the EDM, in particular for the per-particle Hungarian
cost in Fig. 8.

Acknowledging that matching-based comparison is the easiest method for
evaluating the performance of constituent modeling, it’s important to note a
few flaws that can potentially introduce biases. Firstly, particles without a
match are not included in the comparison, and since the matching is performed
with the truth, a substantial number of particles can be missed, as evident
from Fig. 3. Additionally, in some cases, a generated event may have the same
number of particles as the truth, but one particle was not reconstructed,
while a particle from noise with very different features was created. This
scenario can occur in reality, but the matching would pair the unreconstructed
truth particle with the fake one, resulting in a high Hungarian cost. We also
note that when summing the costs of $p_{T}$, $\eta$, $\phi$ and class we weigh
them so every variable on average has approximately the same contribution.
Changing the balance would yield different matching results and can influence
the performance metrics.

Despite these flaws, matching-based comparison still serves as a useful
metric, as these effects will equally occur for both targets and predictions
as long as the cardinality is modeled correctly.

###  V.3 Resolution Performance

A crucial aspect for a good detector simulation is accurate resolution
modeling. The stochastic nature of particle-material interactions leads to
non-trivial smearing of truth particle features. Conducting detector
simulation and reconstruction multiple times on the same truth event results
in different outcomes each time. To evaluate how well our model reproduces
this variability, we introduce a replica dataset, similarly to the concept
introduced in Ref. [38]. We generate 10,000 unique truth events, and for each
truth event, we run detector simulation and reconstruction 100 times. After
performing truth matching for all replicas, each truth particle has a matched
distribution of particles. For each truth particle, we compute the mean and
standard deviation of the matched distribution as a proxy for the detector
resolution. This is summarized in bins of truth $p_{T}$ in Fig. 9.

Figure 9: Means and sigmas of the per-particle replica $p_{T}$-distributions
binned by the matched truth particle $p_{T}$. Figure 10: Event display with
replicas. The detector simulation is run 100 times for the displayed truth
event.

While the means are modeled comparably well by both models, Slot Attention
largely underestimates the standard deviations and therefore doesn’t
accurately model the resolution. However, the diffusion models match the
target resolution significantly better. EDM features a slight trend towards
underestimating the resolution.

We note that utilizing the replica dataset for training with an extended loss
function, as demonstrated in Ref. [38], is an option. In such a scenario, the
Slot Attention model would yield better performance in resolution modelling.
However, the diffusion models consistently outperform it. Since we found no
significant advantage for the latter, we opted for the simpler setup without
replicas.

For demonstrative purposes, we present an event display with replicas in Fig.
10. The showcased event was selected because it contains particles across a
broad $p_{T}$ spectrum without exceeding a count of particles that would
classify it as a rare event. The event display aligns with our previous
discussions. Slot Attention fails to accurately model the distribution and
does not replicate the low-energy noise. In contrast, the diffusion models
exhibit a much closer alignment with the target distribution but are prone to
a few noisy outliers.

##  VI Conclusion

In this study, we have advanced a previously introduced approach for fast
simulation of particles within dense environments such as jets. Our method
involves directly generating reconstructed particles from input truth
particles, effectively replacing the conventional processes of detector
simulation and reconstruction in a single step.

Central to our approach is the introduction of the diffusion models, which
leverage cutting-edge diffusion techniques on graph-valued data. We conducted
a comprehensive analysis of the networks’ performance, focusing on three key
aspects: overall set quantities, features of individual set constituents, and
resolution modeling, a significant challenge in detector simulation. Both
diffusion models outperform slot attention with respect to resolution
modelling without the need for a training dataset containing multiple detector
replicas of the same underlying truth event. Overall, the CFM model
outperformed the EDM, seen most clearly for the per-event Hungarian matching
cost (Fig. 8).

To accurately assess resolution performance, we introduced the replica
dataset, enabling us to evaluate generation outcomes multiple times per ground
truth. Our findings demonstrate that the diffusion models significantly
outperform the baseline model in resolution modelling while exhibiting
comparable performance in overall set constituent features.

Looking ahead, our ultimate objective is to develop a fast simulation package
that is universally applicable across various events. Future directions could
involve scaling up the model to process full events and, therefore, handle
higher cardinalities, although the performance implications of such scaling
remain uncertain, particularly for events with several hundred particles.
Alternatively, implementing event partitioning could address the cardinality
issue, but would require training on a diverse range of detector signatures,
including jets of varying energy ranges and isolated objects. Using the
results of this paper, Ref. [40] makes progress on some of these challenges
towards a full end-to-end approach.

## Data and Code Availability

The code for this paper can be found at
https://github.com/dkobylianskii/fastsim-advanced. The COCOA dataset can be
found at Zenodo at https://zenodo.org/records/11383444.

## Acknowledgement

We thank Yaron Lipman and Vinicius Mikuni for fruitful discussions. ED, EG,
NK, DK, and NS are supported by the BSF-NSF grant 2028 and the ISF Research
Center 494. BN is supported by the U.S. Department of Energy (DOE), Office of
Science under contract DE-AC02-05CH11231.

## References

  * Brüning  _et al._ [2022] O. Brüning, H. Gray, K. Klein, M. Lamont, M. Narain, R. Polifka, and L. Rossi, The scientific potential and technological challenges of the High-Luminosity Large Hadron Collider program, Rept. Prog. Phys. 85, 046201 (2022).
  * ATLAS Collaboration [2010] ATLAS Collaboration, The ATLAS Simulation Infrastructure, Eur. Phys. J. C 70, 823 (2010).
  * Hildreth  _et al._ [2015] M. Hildreth, V. N. Ivanchenko, D. J. Lange, and M. J. Kortelainen, CMS Full Simulation for Run-2, J. Phys. Conf. Ser. 664, 072022 (2015).
  * Agostinelli  _et al._ [2003] S. Agostinelli _et al._ (GEANT4), GEANT4: A simulation toolkit, Nucl. Instrum. Meth. A506, 250 (2003).
  * Allison  _et al._ [2006] J. Allison _et al._ , Geant4 developments and applications, IEEE Transactions on Nuclear Science 53, 270 (2006).
  * Allison  _et al._ [2016] J. Allison _et al._ , Recent developments in geant4, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 835, 186 (2016).
  * Paganini  _et al._ [2018] M. Paganini, L. de Oliveira, and B. Nachman, Accelerating Science with Generative Adversarial Networks: An Application to 3D Particle Showers in Multilayer Calorimeters, Phys. Rev. Lett. 120, 042003 (2018), arXiv:1705.02355 [hep-ex] .
  * Hashemi and Krause [2023] B. Hashemi and C. Krause, Deep Generative Models for Detector Signature Simulation: An Analytical Taxonomy (2023), arXiv:2312.09597 [physics.ins-det] .
  * Adelmann  _et al._ [2022] A. Adelmann, W. Hopkins, E. Kourlitis, M. Kagan, G. Kasieczka, C. Krause, D. Shih, V. Mikuni, B. Nachman, K. Pedro, and D. Winklehner, New directions for surrogate models and differentiable programming for high energy physics detector simulation (2022), arXiv:2203.08806 [hep-ph] .
  * Badger  _et al._ [2023] S. Badger _et al._ , Machine learning and LHC event generation, SciPost Phys. 14, 079 (2023), arXiv:2203.07460 [hep-ph] .
  * Aad  _et al._ [2022] G. Aad _et al._ (ATLAS), AtlFast3: the next generation of fast simulation in ATLAS, Comput. Softw. Big Sci. 6, 7 (2022).
  * Giannelli  _et al._ [2022] M. F. Giannelli, G. Kasieczka, C. Krause, B. Nachman, D. Salamani, D. Shih, and A. Zaborowska, Fast Calorimeter Simulation Challenge 2022 - Dataset 1, 10.5281/zenodo.6368338 (2022).
  * Faucci Giannelli  _et al._ [2022a] M. Faucci Giannelli, G. Kasieczka, C. Krause, B. Nachman, D. Salamani, D. Shih, and A. Zaborowska, Fast Calorimeter Simulation Challenge 2022 - Dataset 2, 10.5281/zenodo.6366271 (2022a).
  * Faucci Giannelli  _et al._ [2022b] M. Faucci Giannelli, G. Kasieczka, C. Krause, B. Nachman, D. Salamani, D. Shih, and A. Zaborowska, Fast Calorimeter Simulation Challenge 2022 - Dataset 3, 10.5281/zenodo.6366324 (2022b).
  * Giannelli and Zhang [2023] M. F. Giannelli and R. Zhang, CaloShowerGAN, a Generative Adversarial Networks model for fast calorimeter shower simulation (2023), arXiv:2309.06515 [physics.ins-det] .
  * Scham  _et al._ [2023] M. A. W. Scham, D. Krücker, B. Käch, and K. Borras, DeepTreeGAN: Fast Generation of High Dimensional Point Clouds (2023), arXiv:2311.12616 [hep-ex] .
  * Krause and Shih [2023] C. Krause and D. Shih, Fast and accurate simulations of calorimeter showers with normalizing flows, Phys. Rev. D 107, 113003 (2023).
  * Ernst  _et al._ [2023] F. Ernst, L. Favaro, C. Krause, T. Plehn, and D. Shih, Normalizing Flows for High-Dimensional Detector Simulations (2023), arXiv:2312.09290 [hep-ph] .
  * Mikuni and Nachman [2022] V. Mikuni and B. Nachman, Score-based generative models for calorimeter shower simulation, Phys. Rev. D 106, 092009 (2022), arXiv:2206.11898 [hep-ph] .
  * Amram and Pedro [2023] O. Amram and K. Pedro, Denoising diffusion models with geometry adaptation for high fidelity calorimeter simulation, Phys. Rev. D 108, 072014 (2023).
  * Buhmann  _et al._ [2024] E. Buhmann, F. Gaede, G. Kasieczka, A. Korol, W. Korcari, K. Krüger, and P. McKeown, CaloClouds II: ultra-fast geometry-independent highly-granular calorimeter simulation, Journal of Instrumentation 19 (04), P04020.
  * Mikuni and Nachman [2024a] V. Mikuni and B. Nachman, CaloScore v2: single-shot calorimeter shower simulation with diffusion models, Journal of Instrumentation 19 (02), P02001.
  * Kobylianskii  _et al._ [2024] D. Kobylianskii, N. Soybelman, E. Dreyer, and E. Gross, CaloGraph: Graph-based diffusion model for fast shower generation in calorimeters with irregular geometry (2024), arXiv:2402.11575 [hep-ex] .
  * CER [2022] _ATLAS Software and Computing HL-LHC Roadmap_, Tech. Rep. (CERN, Geneva, 2022).
  * Software and Computing [2022] C. O. Software and Computing,  _CMS Phase-2 Computing Model: Update Document_, Tech. Rep. (CERN, Geneva, 2022).
  * Hashemi  _et al._ [2019] B. Hashemi, N. Amin, K. Datta, D. Olivito, and M. Pierini, LHC analysis-specific datasets with Generative Adversarial Networks, (2019), arXiv:1901.05282 [hep-ex] .
  * Buhmann  _et al._ [2023a] E. Buhmann, G. Kasieczka, and J. Thaler, EPiC-GAN: Equivariant point cloud generation for particle jets, SciPost Phys. 15, 130 (2023a).
  * Mikuni  _et al._ [2023] V. Mikuni, B. Nachman, and M. Pettee, Fast point cloud generation with diffusion models in high energy physics, Phys. Rev. D 108, 036025 (2023).
  * Buhmann  _et al._ [2023b] E. Buhmann, C. Ewen, D. A. Faroughy, T. Golling, G. Kasieczka, M. Leigh, G. Quétant, J. A. Raine, D. Sengupta, and D. Shih, EPiC-ly Fast Particle Cloud Generation with Flow-Matching and Diffusion (2023b), arXiv:2310.00049 [hep-ph] .
  * Leigh  _et al._ [2024] M. Leigh, D. Sengupta, J. A. Raine, G. Quétant, and T. Golling, Faster diffusion model with improved quality for particle cloud generation, Phys. Rev. D 109, 012010 (2024).
  * Kansal  _et al._ [2021] R. Kansal, J. Duarte, H. Su, B. Orzari, T. Tomei, M. Pierini, M. Touranakou, D. Gunopulos, _et al._ , Particle cloud generation with message passing generative adversarial networks, Advances in Neural Information Processing Systems 34, 23858 (2021).
  * Vaselli  _et al._ [2024a] F. Vaselli, F. Cattafesta, P. Asenov, and A. Rizzi, End-to-end simulation of particle physics events with flow matching and generator oversampling (2024a), arXiv:2402.13684 [hep-ex] .
  * Birk  _et al._ [2024] J. Birk, A. Hallin, and G. Kasieczka, OmniJet-$\alpha$: The first cross-task foundation model for particle physics, (2024), arXiv:2403.05618 [hep-ph] .
  * Mikuni and Nachman [2024b] V. Mikuni and B. Nachman, OmniLearn: A Method to Simultaneously Facilitate All Jet Physics Tasks, (2024b), arXiv:2404.16091 [hep-ph] .
  * Vaselli  _et al._ [2023] F. Vaselli, A. Rizzi, F. Cattafesta, and G. Cicconofri (CMS),  _FlashSim prototype: an end-to-end fast simulation using Normalizing Flow_, Tech. Rep. (CERN, Geneva, 2023).
  * Vaselli  _et al._ [2024b] F. Vaselli, F. Cattafesta, P. Asenov, and A. Rizzi, End-to-end simulation of particle physics events with Flow Matching and generator Oversampling, (2024b), arXiv:2402.13684 [hep-ex] .
  * Howard  _et al._ [2022] J. N. Howard, S. Mandt, D. Whiteson, and Y. Yang, Learning to simulate high energy particle collisions from unlabeled data, Sci. Rep. 12, 7567 (2022), arXiv:2101.08944 [hep-ph] .
  * Soybelman  _et al._ [2023] N. Soybelman, N. Kakati, L. Heinrich, F. A. Di Bello, E. Dreyer, S. Ganguly, E. Gross, M. Kado, and J. Shlomi, Set-conditional set generation for particle physics, Mach. Learn. Sci. Tech. 4, 045036 (2023).
  * Chen  _et al._ [2019] R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. Duvenaud, Neural ordinary differential equations (2019), arXiv:1806.07366 [cs.LG] .
  * Dreyer  _et al._ [2024] E. Dreyer _et al._ , Parnassus: An Automated Approach to Accurate, Precise, and Fast Detector Simulation and Reconstruction, (2024), arXiv:2405.xxxxx [hep-ex] .
  * Selvaggi [2014] M. Selvaggi, DELPHES 3: A modular framework for fast-simulation of generic collider experiments, J. Phys. Conf. Ser. 523, 012033 (2014).
  * Sjostrand  _et al._ [2008] T. Sjostrand, S. Mrenna, and P. Z. Skands, A Brief Introduction to PYTHIA 8.1, Comput. Phys. Commun. 178, 852 (2008).
  * Charkin-Gorbulin  _et al._ [2023] A. Charkin-Gorbulin _et al._ , Configurable calorimeter simulation for AI applications, Mach. Learn. Sci. Tech. 4, 035042 (2023).
  * Di Bello  _et al._ [2023] F. A. Di Bello, E. Dreyer, S. Ganguly, E. Gross, L. Heinrich, A. Ivina, M. Kado, N. Kakati, L. Santi, J. Shlomi, and M. Tusoni, Reconstructing particles in jets using set transformer and hypergraph prediction networks, The European Physical Journal C 83, 10.1140/epjc/s10052-023-11677-7 (2023).
  * Zhang  _et al._ [2023] D. W. Zhang, G. J. Burghouts, and C. G. M. Snoek, Pruning edges and gradients to learn hypergraphs from larger sets (2023), arXiv:2106.13919 [cs.LG] .
  * Locatello  _et al._ [2020] F. Locatello, D. Weissenborn, T. Unterthiner, A. Mahendran, G. Heigold, J. Uszkoreit, A. Dosovitskiy, and T. Kipf, Object-centric learning with slot attention, CoRR abs/2006.15055 (2020), 2006.15055 .
  * Arjovsky  _et al._ [2017] M. Arjovsky, S. Chintala, and L. Bottou, Wasserstein gan (2017), arXiv:1701.07875 [stat.ML] .
  * Liu  _et al._ [2022a] L. Liu, Y. Ren, Z. Lin, and Z. Zhao, Pseudo Numerical Methods for Diffusion Models on Manifolds (2022a), arXiv:2202.09778 [cs, math, stat].
  * Karras  _et al._ [2022] T. Karras, M. Aittala, T. Aila, and S. Laine, Elucidating the design space of diffusion-based generative models (2022), arXiv:2206.00364 [cs.CV] .
  * Lipman  _et al._ [2023] Y. Lipman, R. T. Q. Chen, H. Ben-Hamu, M. Nickel, and M. Le, Flow Matching for Generative Modeling (2023), arXiv:2210.02747 [cs, stat].
  * Tong  _et al._ [2024] A. Tong, K. Fatras, N. Malkin, G. Huguet, Y. Zhang, J. Rector-Brooks, G. Wolf, and Y. Bengio, Improving and generalizing flow-based generative models with minibatch optimal transport (2024), arXiv:2302.00482 [cs.LG] .
  * Liu  _et al._ [2022b] X. Liu, C. Gong, and Q. Liu, Flow Straight and Fast: Learning to Generate and Transfer Data with Rectified Flow (2022b), arXiv:2209.03003 [cs].
  * Dax  _et al._ [2023] M. Dax, J. Wildberger, S. Buchholz, S. R. Green, J. H. Macke, and B. Schölkopf, Flow Matching for Scalable Simulation-Based Inference (2023), arXiv:2305.17161 [cs].
  * Chen  _et al._ [2023] T. Chen, R. Zhang, and G. Hinton, Analog bits: Generating discrete data using diffusion models with self-conditioning (2023), arXiv:2208.04202 [cs.CV] .
  * Loshchilov and Hutter [2019] I. Loshchilov and F. Hutter, Decoupled weight decay regularization (2019), arXiv:1711.05101 [cs.LG] .
  * Kuhn [1955] H. W. Kuhn, The Hungarian method for the assignment problem,  _Naval research logistics quarterly_ , Naval research logistics quarterly 2, 83 (1955).