# First equilibrium reconstruction for ITER with the code NICE

B. Faugeras11footnotetext: Corresponding author. J. Blum  C. Boulbe

###### Abstract

In this short paper we present the first application of the IMAS compatible
code NICE to equilibrium reconstrution for ITER geometry. The inverse problem
is formulated as a least square problem and the numerical methods implemented
in NICE in order to solve it are presented. The results of a numerical
experiment are shown: a reference equilibrium is computed from which a set of
synthetic magnetic measurements are extracted. Then these measurements are
used successfully to reconstruct the equilibrium of the plasma.

##  1 Introduction

Numerical reconstruction of the plasma equilibrium in a tokamak is an
important and long standing subject in fusion plasma science [34, 39, 27, 37,
24, 2]. The resolution of this inverse problem consists in the computation of
the poloidal flux function and of the plasma boundary as well as the
identification of two non-linear source term functions known as $p^{\prime}$
and $ff^{\prime}$ in the Grad-Shafranov equation [18, 35, 26]. It is needed on
the one hand for a posteriori analysis of experimental equilibrium
configurations and on the other hand for real time control of the plasma
during a discharge. The basic set of measurements needed and used are magnetic
probes and flux loops which provide values of the poloidal magnetic field and
flux at several points surrounding the vacuum vessel and the plasma. All free
boundary reconstruction codes (e.g. [25, 6, 28, 40, 3, 4, 10, 9, 30])
primarily use these magnetic measurements which proved to be sufficient to
identify correctly the plasma boundary and the averaged plasma current density
profile [4]. The goal of this paper is to present a first test of the
adaptation of the numerical tools developed by the authors, namely the code
NICE [9], to the foreseen ITER configuration. NICE stands for "Newton direct
and Inverse Computation for Equilibrium".

Next Section 2 is devoted to the formulation of the direct model and the
inverse problem. In Section 3 we discuss the numerical methods which we have
developed for their resolution and finally in Section 4 a first test numerical
experiment is presented for ITER configuration.

##  2 Inverse problem formulation

###  2.1 Free-boundary plasma equilibrium

The equations which govern the equilibrium of a plasma in the presence of a
magnetic field in a tokamak are on the one hand Maxwell’s equations satisfied
in the whole of space (including the plasma):

$\displaystyle\nabla\cdot\mathbf{B}=0,$
$\displaystyle\nabla\times(\frac{\mathbf{B}}{\mu})=\mathbf{j},$ (2.1)

and on the other hand the equilibrium equation for the plasma itself

$\nabla p=\mathbf{j}\times\mathbf{B},$ (2.2)

where $\mathbf{B}$ is the magnetic field, $\mu$ is the magnetic permeability,
$p$ is the kinetic pressure and $\mathbf{j}$ is the current density. We refer
to standard text books (e.g. [15, 2, 38, 17, 23]) and to [20] for details of
the derivation and only state the needed equations in what follows which is a
summary of what can be found in [9].

Introducing a cylindrical coordinate system
$({\boldsymbol{e}}_{r},{\boldsymbol{e}}_{\phi},{\boldsymbol{e}}_{z})$ ($r=0$
is the major axis of the tokamak torus) and assuming axial symmetry equations
(2.1) and (2.2) reduce to the following equation for the magnetic poloidal
flux $\psi(r,z)$ in the poloidal plane
$\Omega_{\infty}=(0,\infty)\times(-\infty,\infty)$:

$-\Delta^{*}\psi=j_{\phi},$ (2.3)

where $j_{\phi}$ is the toroidal component of $\mathbf{j}$, and the second
order elliptic differential operator $\Delta^{*}$ is defined by

$\Delta^{*}\psi:=\partial_{r}\left(\frac{1}{\mu_{0}r}\partial_{r}\psi\right)+\partial_{z}\left(\frac{1}{\mu_{0}r}\partial_{z}\psi\right):=\nabla\cdot\left(\frac{1}{\mu_{0}r}\nabla\psi\right).$
(2.4)

Here $\nabla$ is the 2D operator in the $(r,z)$-plane and $\mu_{0}$ is the
magnetic permeability of vacuum (in this work we consider only air-transformer
tokamaks such as ITER).

The magnetic field can be decomposed in poloidal and toroidal components

$\displaystyle\mathbf{B}=\mathbf{B}_{p}+\mathbf{B}_{\phi},$
$\displaystyle\mathbf{B}_{p}=\frac{1}{r}[\nabla\psi\times{\boldsymbol{e}}_{\phi}],\
$
$\displaystyle{\boldsymbol{B}}_{\phi}=B_{\phi}\mathbf{e}_{\phi}=\frac{f}{r}{\boldsymbol{e}}_{\phi},$
(2.5)

where $f$ is the diamagnetic function. Equation (2.5) shows that the magnetic
surfaces are generated by the rotation of the iso-flux lines around the axis
of the torus.

The toroidal component of the current density $j_{\phi}$ is zero everywhere
outside the plasma domain and the poloidal field coils (and possibly the
passive structures). The different sub-domains of the poloidal plane of a
schematic tokamak (see Fig. 1) as well as the corresponding expression for
$j_{\phi}$ are described below:

rz$\Omega_{\mathrm{c}_{i}}$$\Omega_{\mathrm{ps}}$$\Omega_{\mathrm{p}}$$\Omega_{\mathrm{L}}$$\partial\Omega_{\mathrm{L}}$$0$
(a)  (b)

Figure 1: Left (a): schematic representation of the poloidal plane of a
tokamak. $\Omega_{\mathrm{p}}$ is the plasma domain, $\Omega_{\mathrm{L}}$ is
the limiter domain accessible to the plasma, $\Omega_{\mathrm{c}_{i}}$
represent poloidal field coils and the central solenoid (corresponding to the
IMAS pf$\\_$active IDS), $\Omega_{\mathrm{ps}}$ the passive structures. Right
(b): example of a plasma whose boundary is defined by the contact with limiter
(left) or by the presence of an X-point (right).

  * $\bullet$

$\Omega_{\mathrm{L}}$ is the domain accessible to the plasma. Its boundary is
the limiter $\partial\Omega_{\mathrm{L}}$.

  * $\bullet$

$\Omega_{\mathrm{p}}$ is the plasma domain where equations (2.2) and (2.1)
imply that $p$ and $f$ are constant on each magnetic surface i.e. $p=p(\psi)$
and $f=f(\psi)$. One then deduces the so-called Grad-Shafranov equilibrium
equation in the plasma [18, 35, 26]

$-\Delta^{*}\psi=rp^{\prime}(\psi)+\frac{1}{\mu_{0}r}(ff^{\prime})(\psi).$
(2.6)

The right-hand side of (2.6) is the toroidal component $j_{\phi}$ of the
current density in the plasma.

The plasma domain is unknown, i.e.
$\Omega_{\mathrm{p}}=\Omega_{\mathrm{p}}(\psi)$, and this is a free boundary
problem. This domain is defined by its boundary which is the outermost closed
$\psi$ iso-contour contained within the limiter $\Omega_{\mathrm{L}}$. The
plasma can either be limited if this iso-contour is tangent to the limiter
$\partial\Omega_{\mathrm{L}}$ or defined by the presence of an X-point (see
Fig. 1). Functions $p^{\prime}$ and $ff^{\prime}$ are zero outside
$\Omega_{\mathrm{p}}$.

The current density is non-linear in $\psi$ due to the non-linear functions
$p^{\prime}$ and $ff^{\prime}$ and the definition of the plasma domain
$\Omega_{\mathrm{p}}(\psi)$. While $\Omega_{\mathrm{p}}(\psi)$ is fully
determined for a given $\psi$, the two functions $p^{\prime}$ and
$ff^{\prime}$ are not determined in this modelization. It is the goal of the
inverse equilibrium reconstruction problem to determine them. For now let us
consider that we are given two functions $\mathcal{A}(\psi_{\mathrm{N}})$ and
$\mathcal{B}(\psi_{\mathrm{N}})$ such that

$j_{\phi}=\lambda(\displaystyle\frac{r}{r_{0}}\mathcal{A}(\psi_{\mathrm{N}})+\displaystyle\frac{r_{0}}{r}\mathcal{B}(\psi_{\mathrm{N}})).$
(2.7)

Here $r_{0}$ is the major radius of the tokamak vacuum chamber and $\lambda$
is a scaling coefficient. The normalized poloidal flux
$\psi_{\mathrm{N}}(r,z)$ is

$\psi_{\mathrm{N}}(r,z)=\frac{\psi(r,z)-\psi_{\mathrm{a}}(\psi)}{\psi_{\mathrm{b}}(\psi)-\psi_{\mathrm{a}}(\psi)}\,.$
(2.8)

with $\psi_{\mathrm{a}}$ and $\psi_{\mathrm{b}}$ being the flux values at the
magnetic axis and at the boundary of the plasma.

  * $\bullet$

Each domain $\Omega_{\mathrm{c}_{i}}$, $i=1,\dots N_{C}$ represent one of the
$N_{C}$ coils carrying currents that is to say poloidal field coils as well as
the central solenoid (this corresponds to the IMAS pf$\\_$active IDS). The
expression of the current density in the $i$-th coil is

$j_{\phi}=\displaystyle\frac{I_{i}}{S_{i}},$ (2.9)

where $S_{i}$ is the section area of the coil and $I_{i}$ is a given measured
current.

  * $\bullet$

$\Omega_{\mathrm{ps}}$ represents passive structures where the induced current
density is assumed to be $0$ in this work but can be considered to be measured
and given in the same form as Eq. (2.9)

To sum up, given functions $\mathcal{A}$ and $\mathcal{B}$, and currents
$\mathbf{I}=\\{I_{i}\\}_{i=1}^{N_{C}}$ in the coils, the free-boundary
equilibrium equation for $\psi(r,z)$ on $\Omega_{\infty}$ is the following
non-linear boundary value problem

$\left\\{\begin{array}[]{l}-\Delta^{*}\psi=\begin{cases}\lambda(\displaystyle\frac{r}{r_{0}}\mathcal{A}(\psi_{\mathrm{N}})+\displaystyle\frac{r_{0}}{r}\mathcal{B}(\psi_{\mathrm{N}}))&\text{in
}\Omega_{\mathrm{p}}(\psi),\\\ \displaystyle\frac{I_{i}}{S_{i}}&\text{in
}\Omega_{\mathrm{c}_{i}},\\\ 0&\text{elsewhere},\end{cases}\\\
\psi(0,z,t)=0,\\\ \lim\limits_{\|(r,z)\|\to+\infty}\psi(r,z)=0.\\\
\end{array}\right.$ (2.10)

This formulation on an infinite domain is not used directly in computations
where we use finite elements on a truncated bounded domain. The infinite
domain is reduced to a semi circular computational domain by an uncoupling
procedure [1, 16]. We chose a semi-circle $\Gamma$ of radius $\rho_{\Gamma}$
surrounding the coil domains $\Omega_{\mathrm{c}_{i}}$ and define the
computation domain $\Omega$ having boundary
$\partial\Omega=\Gamma\cup\Gamma_{0}$, where
$\Gamma_{0}=\\{(0,z),\,z\in[-\rho_{\Gamma},\rho_{\Gamma}]\\}$.

The weak formulation of the equilibrium problem on which the finite element
method relies uses a function space $V$ defined in [20] and can be written as:
Given function $\mathcal{A}$ and $\mathcal{B}$, and currents $\mathbf{I}$,
find $\psi\in V$ such that for all $\xi\in V$

$\mathsf{a}(\psi,\xi)+\mathsf{c}(\psi,\xi)-\,\mathsf{J}_{\mathrm{p}}(\psi,\xi;\mathcal{A},\mathcal{B})-\mathsf{\ell}(\mathbf{I},\xi)=0,$
(2.11)

where

$\displaystyle\mathsf{a}(\psi,\xi)$
$\displaystyle:=\int_{\Omega}\frac{1}{\mu_{0}r}\nabla\psi\cdot\nabla\xi\,drdz,$
(2.12)
$\displaystyle\mathsf{J}_{\mathrm{p}}(\psi,\xi;\mathcal{A},\mathcal{B})$
$\displaystyle:=\int_{\Omega_{\mathrm{p}}(\psi)}\lambda\left(\frac{r}{r_{0}}\mathcal{A}(\psi_{\mathrm{N}})+\frac{r_{0}}{r}\mathcal{B}(\psi_{\mathrm{N}})\right)\xi\,drdz,$
$\displaystyle\mathsf{\ell}(\mathbf{I},\xi)$
$\displaystyle:=\sum_{i=1}^{N_{C}}\frac{I_{i}}{S_{i}}\int_{\Omega_{\mathrm{c}_{i}}}\xi\,drdz,$

and the bilinear form $\mathsf{c}:V\times V\to\mathbb{R}$ is accounting for
the boundary conditions at infinity. We refer to [20] for its precise
expression and to [19, Chapter 2.4] for the details on its the derivation.
Alternative approaches for the incorporation of boundary conditions at
infinity were more recently presented in [12].

###  2.2 The inverse reconstruction problem

Magnetics constitute the basic set of experimental measurements used in
equilibrium reconstruction for the identification of functions $\mathcal{A}$
and $\mathcal{B}$. They consist in measurements of projections of the poloidal
magnetic field, $\mathbf{B}_{p}\cdot\mathbf{d}$ at several locations around
the vacuum vessel of the tokamak (the unit vector $\mathbf{d}$ varies with
each B-probe) and of flux loops measurements, noted $F(\psi)$, at several
locations too.

At this point we have defined a direct model given by the equilibrium equation
(2.11), control variables $\mathcal{A}$ and $\mathcal{B}$, and measurements to
which are attached experimental errors represented by the standard deviations
$\sigma$s in Eq. (2.14) below. The identification problem can now be
formulated as a constrained minimization problem for the following cost
function:

$J(\psi,\mathcal{A},\mathcal{B}):=J_{obs}(\psi)+R(\mathcal{A})+R(\mathcal{B}),$
(2.13)

where the least-square misfit term is

$\begin{array}[]{l}J_{obs}(\psi):=\displaystyle\sum_{i=1}^{N_{B}}\frac{1}{2\sigma_{Bi}^{2}}((\mathbf{B}_{p}(r_{i},z_{i})\cdot\mathbf{d}_{i})-B^{i}_{p,obs})^{2}+\displaystyle\sum_{i=1}^{N_{F}}\frac{1}{2\sigma_{Fi}^{2}}(F(\psi)_{i}-F^{i}_{obs})^{2}\\\
\end{array}$ (2.14)

and the regularization term $R(\mathcal{A})$ is defined as

$R(\mathcal{A}):=\displaystyle\frac{\varepsilon_{\mathcal{A}}}{2}\int_{0}^{1}(\mathcal{A}^{\prime\prime}(x))^{2}dx+\frac{\alpha_{\mathcal{A}}}{2}|\mathcal{A}(1)|^{2}.$
(2.15)

The $R(\mathcal{B})$ term is defined similarly. Parameters $\varepsilon$
enable to tune the smoothness of the identified functions whereas parameters
$\alpha$ tune the penalization to zero of their value on the plasma boundary.
The inverse problem consists in the minimization of $J$ (2.13) under the
constraint of the model equation (2.11).

##  3 Numerical methods

###  3.1 Discretization of the direct model

Equilibrium equation (2.11) is discretized using a P1 finite element method
based on a triangular mesh. From now on let us also assume that functions
$\mathcal{A}$ and $\mathcal{B}$ are decomposed in a basis of functions
${\phi_{i}}$ defined on $[0,1]$. We use cubic spline functions in this work
and

$\displaystyle\mathcal{A}(x)=\sum_{i=1}^{N_{\mathcal{A}}}u_{\mathcal{A}i}\phi_{i}(x),$
$\displaystyle\mathcal{B}(x)=\sum_{i=1}^{N_{\mathcal{B}}}u_{\mathcal{B}i}\phi_{i}(x).$
(3.1)

Let us denote $\mathbf{u}=(\mathbf{u}_{\mathcal{A}},\mathbf{u}_{\mathcal{B}})$
of size $N_{u}$ the vector of degrees of freedom of $\mathcal{A}$ and
$\mathcal{B}$ in the decomposition basis. Classically approximating $\psi$ by
$\psi_{h}=\displaystyle\sum_{i=1}^{N}\psi_{i}\lambda_{i}(r,z)$ on the finite
element approximation space as well as the operators of (2.11) and taking all
basis elements $\lambda_{i}$ as test functions leads to the following non-
linear system of $N$ equations:

$({\boldsymbol{A}}+{\boldsymbol{C}})\boldsymbol{\psi}-\mathbf{J}_{p}(\boldsymbol{\psi},\mathbf{u})-{\boldsymbol{L}}\mathbf{I}=0$
(3.2)

where $\boldsymbol{\psi}$ denotes the vector of finite element coefficients
$\\{\psi_{i}\\}_{i=1}^{N}$ and other notations are obvious.

###  3.2 The discrete identification problem

Using the discrete variables of the preceding section, cost function (2.13)
can also be discretized leading to the following expression

$\begin{array}[]{l}J(\boldsymbol{\psi},\mathbf{u}):=\displaystyle\frac{1}{2}||{\boldsymbol{H}}\boldsymbol{\psi}-\mathbf{m})||^{2}+\displaystyle\frac{1}{2}||{\boldsymbol{R}}\mathbf{u}||^{2}.\end{array}$
(3.3)

In order to lighten notations the $\displaystyle\frac{1}{\sigma}$ terms have
been dropped and are assumed to be included in the observation operator
${\boldsymbol{H}}$ and in the measurements $\mathbf{m}$. The last term
involving matrice ${\boldsymbol{R}}$ is the discretization of the
regularization terms in which we have gathered the contributions from
functions $\mathcal{A}$ and $\mathcal{B}$.

The discrete identification problem can now be stated as

$\underset{\boldsymbol{\psi},\mathbf{u}}{\min}J(\boldsymbol{\psi},\mathbf{u})$
(3.4)

subject to the constraint of the non-linear model

$({\boldsymbol{A}}+{\boldsymbol{C}})\boldsymbol{\psi}-\mathbf{J}_{p}(\boldsymbol{\psi},\mathbf{u})-{\boldsymbol{L}}\mathbf{I}=0$
(3.5)

This problem is solved thanks to a quasi-SQP algorithm with reduced Hessian
(QSQP). SQP methods are well documented [31, 21] and for fusion application we
refer to [5, Appendix A] and [9, 11]. An SQP method can be seen as a Newton
method to solve the non-linear system given by the fisrt order optimality
condition for the Lagrangian of the PDE-constrained optimization problem.

The QSQP method we use is the following 2 steps iterative algorithm:

  1. 1.

control variable update step

${\boldsymbol{M}}(\mathbf{u}^{n+1}-\mathbf{u}^{n})=-\mathbf{h}$ (3.6)

  2. 2.

state variable update step

$\boldsymbol{\psi}^{n+1}=\boldsymbol{\psi}^{n}+\delta\boldsymbol{\psi}+{\boldsymbol{S}}(\mathbf{u}^{n+1}-\mathbf{u}^{n})$
(3.7)

where

$\delta\boldsymbol{\psi}=-[({\boldsymbol{A}}+{\boldsymbol{C}})\boldsymbol{\psi}^{n}-D_{\boldsymbol{\psi}}\mathbf{J}_{p}(\boldsymbol{\psi}^{n},\mathbf{u}^{n})]^{-1}[({\boldsymbol{A}}+{\boldsymbol{C}})\boldsymbol{\psi}^{n}-\mathbf{J}_{p}(\boldsymbol{\psi}^{n},\mathbf{u}^{n})-{\boldsymbol{L}}\mathbf{I}],$
(3.8)
${\boldsymbol{S}}=-[({\boldsymbol{A}}+{\boldsymbol{C}})\boldsymbol{\psi}^{n}-D_{\boldsymbol{\psi}}\mathbf{J}_{p}(\boldsymbol{\psi}^{n},\mathbf{u}^{n})]^{-1}[-D_{\mathbf{u}}\mathbf{J}_{p}(\boldsymbol{\psi}^{n},\mathbf{u}^{n})],$
(3.9)

${\boldsymbol{M}}={\boldsymbol{R}}^{T}{\boldsymbol{R}}+{\boldsymbol{S}}^{T}{\boldsymbol{H}}^{T}{\boldsymbol{H}}{\boldsymbol{S}}$
(3.10)

and

$\mathbf{h}={\boldsymbol{R}}^{T}{\boldsymbol{R}}\mathbf{u}^{n}+{\boldsymbol{S}}^{T}{\boldsymbol{H}}^{T}({\boldsymbol{H}}\boldsymbol{\psi}^{n}-\mathbf{m})+{\boldsymbol{S}}^{T}{\boldsymbol{H}}^{T}{\boldsymbol{H}}\delta\boldsymbol{\psi}$
(3.11)

At each iteration this algorithm demands the resolution of $N_{u}+1$ linear
systems (3.8)-(3.9) of size $N$ involving the same matrix with different righ-
hand sides which can be done very efficiently and of one smaller linear system
of size $N_{u}$ in (3.6).

The performance of the QSQP method used for the resolution of the
identification problem relies on the accuracy of the derivative terms
$D_{\psi}\mathbf{J}_{p}$, $D_{\mathbf{u}}\mathbf{J}_{p}$. In this work we have
implemented the exact derivatives of the fully discretized operators. This
essential but very technical work is not further detailed here and we refer to
[20, section 3.2 and 3.3] for details.

##  4 Numerical experiment

The numerical methods presented in the previous sections are implemented in
the code NICE [9] with which the following numerical experiment is conducted.
ITER machine description providing the points defining the limiter contour,
the poloidal field coils and a description of the magnetic sensors (196
poloidal magnetic field probes and 22 differential flux loops) are read from
the wall, pf$\\_$active and magnetics Interface Data Structure (IDS) from the
ITER Integrated Modelling and Analysis Suite (IMAS) [22]. From this machine
description NICE builds the triangular mesh used with the finite elements
computations.

We generate a reference equilibrium by solving the so-called inverse static
equilibrium problem [20, 9] that is to say finding the currents in the coils
giving a desired prescribed plasma boundary. For this reference computation
the unknown function are given analytically as
$\mathcal{A}(x)=(1-x^{1.5})^{0.9}$ and $\mathcal{B}(x)=(1-x^{0.9})^{1.5}$. The
scaling factor $\lambda$ in the current density (2.7) is computed such that
the total plasma current is $I_{p}=12$ $[MA]$. The vacuum toroidal field is
$B_{0}=5.3$ $[T]$ and $r_{0}=6.2$ $[m]$. Synthetic magnetic measurements are
computed from this reference equilibrium.

Then in a second step these measurements are plugged in cost function (3.3)
and the optimization problem is solved using the iterative QSQP algorithm
presented above. The initial guess for this resolution consists in a given
circular plasma domain in which the flux $\psi$ is a constant and outside of
which it is $0$, as well as affine functions
$\mathcal{A}(x)=\mathcal{B}(x)=1-x$. Convergence is assumed when the relative
residue for the vector of unknowns $\mathbf{X}=(\boldsymbol{\psi},\mathbf{u})$
satisfies
${||\mathbf{X}^{k+1}-\mathbf{X}^{k}||}/{||\mathbf{X}^{k}||}<10^{-12}$.

The finite element mesh is composed of $N=11856$ nodes among which $11772$
correspond to free values of $\psi$ (the remaining correspond to the imposed
boundary condition $\psi=0$ on the axis $r=0$). Each function to be identified
is decomposed in $11$ cubic splines defined on $[0,1]$ with knots at
$0,0.1,...,1$. Therefore $\mathbf{X}$ is a vector of size $11772+2\times
11=11794$.

The regularization parameters $\varepsilon$ are tuned to their lowest value,
typically $10^{-2}$, avoiding oscillations in the reconstructed profiles or
non convergence of the code. The penalization parameters are set to
$\alpha=10^{12}$.

Starting from the initial guess described above, which is far from the
reference solution we want to recover, the algorithm needs 16 iterations to
converge to a relative residue of $4\times 10^{-14}$. In this test
configuration one iteration takes about $800$ ms on a laptop with an Intel CPU
at $2.90$ GHz. This computation time is highly dependent on the mesh size and
is limited by the performance of the linear solver. The reconstructed flux map
is shown on Figure 2. The plasma boundary is perfectly recovered. The fit to
measurements is excellent with a root mean square error of $10^{-4}$ on
magnetic probes and flux loops.

Figure 2: Reconstructed poloidal flux map ($\psi$ in [Wb]). The plasma
boundary with an X-point and the magnetic axis are shown in red.

Figure 3 shows the reference and identified $p^{\prime}$ and $ff^{\prime}$
profiles. A typical feature of the equilibrium reconstruction inverse problem
using only magnetic measurements appears: the computed error bars on the
reconstructed profiles increase when approaching the magnetic axis. However
even if the $p^{\prime}$ and $ff^{\prime}$ profiles are not perfectly
identified for low $\psi_{\mathrm{N}}$ values, the flux surface averaged
current density profile $jtor=<j_{\phi}/r>/<1/r>$ and the safety factor $q$
are very well recovered.

Figure 3: Top row: reference and reconstructed $p^{\prime}$ and $ff^{\prime}$
profiles as function of the normalized flux $\psi_{\mathrm{N}}$ with computed
error bars. Bottom row: same for the average toroidal current density $jtor$
and safety factor $q$ profiles.

##  5 Conclusion

The equilibrium code NICE is ready to use ITER data. A first equilibrium
reconstruction numerical exercice using synthetic magnetic measurements has
been successfully conducted. The code is fully IMAS compatible, can read and
write IDS. It is also ready to use other measurements such as pressure
measurements, motional Stark effect (MSE) as well as interferometry and
polarimetry. The polarimetry Stokes vector modelization can also be used [13].
These internal measurements will most likely be needed to improve the accuracy
of the reconstruction of the $p^{\prime}$ and $ff^{\prime}$ profiles for ITER
discharges. The code NICE is already routinely used at WEST and has been
tested on different tokamaks and validated against other codes within the
EUROfusion program [29, 14, 8, 36, 33]. NICE is a robust code which exhibits
excellent convergence properties thanks to the use of Newton and SQP methods.
The code can be used for equilibrium reconstruction but also for direct static
or evolutive simulations. It can also be used to solve the inverse problems
consisting in finding the currents or voltages in the poloidal field coils
which enable to have a desired plasma shape [9]. Finally the code NICE can be
used with high order finite elements providing a smooth representation of the
magnetic flux and field in the plasma [7, 32].

## Acknowledgment

We would like to thank Masanari Hosokawa, Simon Pinches and Mireille Schneider
for their help in providing the latest ITER machine description. We would also
like to thank the anonymous reviewer for the numerous remarks and constructive
criticism from which the paper has benefitted.

## References

  * [1] R. Albanese, J. Blum, and O. Barbieri.  On the solution of the magnetic flux equation in an infinite domain.  In EPS. 8th Europhysics Conference on Computing in Plasma Physics (1986), pages 41–44, 1986.
  * [2] J. Blum.  Numerical Simulation and Optimal Control in Plasma Physics with Applications to Tokamaks.  Series in Modern Applied Mathematics. Wiley Gauthier-Villars, Paris, 1989\.
  * [3] J. Blum, C. Boulbe, and B. Faugeras.  Real-time plasma equilibrium reconstruction in a tokamak.  In Journal of Physics: Conference Series. Proceedings of the 6th International Conference on Inverse Problems in Engineering: Theory and Practice, volume 135, page 012019. IOP Publishing, 2008.
  * [4] J. Blum, C. Boulbe, and B. Faugeras.  Reconstruction of the equilibrium of the plasma in a tokamak and identification of the current density profile in real time.  Journal of Computational Physics, 231(3):960 – 980, 2012.
  * [5] J. Blum, H Heumann, E. Nardon, and X. Song.  Automating the design of tokamak experiment scenarios.  J. Computational Physics, 394:594–614, 2019.
  * [6] J. Blum, E. Lazzaro, J. O’Rourke, B. Keegan, and Y. Stefan.  Problems and methods of self-consistent reconstruction of tokamak equilibrium profiles from magnetic and polarimetric measurements.  Nuclear Fusion, 30(8):1475, 1990.
  * [7] A. Elarif, B. Faugeras, and F. Rapetti.  Tokamak free-boundary plasma equilibrium computation using finite elements of class C0 and C1 within a mortar element approach.  J. Computational Physics, 439:110388, 2021.  .
  * [8] G. L. Falchetto, P. Strand, R. Coelho, D. Coster, J. Ferreira, T. Jonsson, D Yadikin, R. Dumont, B. Faugeras, J. Hollocombe, P. Huynh, J. Joly, D. Kalupin, E. Lerche, J. Morales, M. Poradzinski, P. Siren, E. Tholerus, D. Van Eester, J. Varje, W. Zwingmann, JET Contributors, and EUROfuion-IM Team.  Multi-machine analysis of eu experiments using the eurofusion integrated modelling (eu-im) framework.  In European Physical Society, editor, 46th EPS Conference on Plasma Physics, volume 43C of Europhysics Conference Abstracts (ECA), Milan, Italy, 2019.
  * [9] B. Faugeras.  An overview of the numerical methods for tokamak plasma equilibrium computation implemented in the NICE code.  Fusion Eng. Design, 160:112020, 2020.  https://hal.archives-ouvertes.fr/hal-02955053.
  * [10] B. Faugeras, J. Blum, C. Boulbe, P. Moreau, and E. Nardon.  2D interpolation and extrapolation of discrete magnetic measurements with toroidal harmonics for equilibrium reconstruction in a Tokamak.  Plasma Phys. Control Fusion, 56:114010, 2014.
  * [11] B. Faugeras, J. Blum, H. Heumann, and C. Boulbe.  Optimal control of a coupled partial and ordinary differential equations system for the assimilation of polarimetry stokes vector measurements in tokamak free-boundary equilibrium reconstruction with application to ITER.  Comput. Phys. Comm., 217(Supplement C):43 – 57, 2017.
  * [12] B. Faugeras and H Heumann.  FEM-BEM coupling methods for tokamak plasma axisymmetric free-boundary equilibrium computations in unbounded domains.  J. Computational Physics, 343(Supplement C):201 – 216, 2017.
  * [13] B. Faugeras, F. Orsitto, and JET Contributors.  Equilibrium reconstruction at JET using Stokes model for polarimetry.  Nuclear Fusion, 58(10):106032, 2018.
  * [14] L. Fleury, H. Ancher, J.F. Artaud, B. Faugeras, M. Geynet, F. Hollocombe, F. Imbeaux, P. Maini, and J. Morales.  WEST plasma reconstruction chain and IMAS related tools.  In SOFT 2020, Croatia, 09 2020.
  * [15] J. P. Freidberg.  Ideal Magnetohydrodynamics.  Plenum US, 1987.
  * [16] G.N. Gatica and G.C. Hsiao.  The uncoupling of boundary integral and finite element methods for nonlinear boundary value problems.  J. Math. Anal. Appl., 189(2):442–461, 1995.
  * [17] Johan Peter Goedbloed and Stefaan Poedts.  Principles of magnetohydrodynamics: with applications to laboratory and astrophysical plasmas.  Cambridge university press, 2004.
  * [18] H. Grad and H. Rubin.  Hydromagnetic equilibria and force-free fields.  Proceedings of the 2nd UN Conf. on the Peaceful Uses of Atomic Energy, 31:190, 1958.
  * [19] Virginie Grandgirard.  Modélisation de l’équilibre d’un plasma de tokamak.  PhD thesis, Université de Franche-Comté, 1999.
  * [20] H. Heumann, J. Blum, C. Boulbe, B. Faugeras, G. Selig, J.-M. Ané, S. Brémond, V. Grangirard, P. Hertout, and E. Nardon.  Quasi-static free-boundary equilibrium of toroidal plasma with CEDRES++: computational methods and applications.  J. Plasma Physics, 2015.
  * [21] M. Hinze, R. Pinnau, M. Ulbrich, and S. Ulbrich.  Optimization with PDE constraints, volume 23 of Mathematical Modelling: Theory and Applications.  Springer, New York, 2009.
  * [22] F. Imbeaux, S.D. Pinches, J.B. Lister, Y. Buravand, T. Casper, B. Duval, B. Guillerminet, M. Hosokawa, W. Houlberg, P. Huynh, S.H. Kim, G. Manduchi, M. Owsiak, B. Palak, M. Plociennik, G. Rouault, O. Sauter, and P. Strand.  Design and first applications of the ITER integrated modelling & analysis suite.  Nuclear Fusion, 55(12):123006, 2015.
  * [23] S.C. Jardin.  Computational methods in plasma physics.  Boca Raton, FL : CRC Press/Taylor & Francis, 2010.
  * [24] L.L. Lao.  Separation of $\beta_{p}$ and $l_{i}$ in tokamaks of non-circular cross-section.  Nuclear Fusion, 25(11):1421, 1985.
  * [25] L.L. Lao, J.R. Ferron, R.J. Geoebner, W. Howl, H.E. St. John, E.J. Strait, and T.S. Taylor.  Equilibrium analysis of current profiles in Tokamaks.  Nuclear Fusion, 30(6):1035, 1990.
  * [26] R. Lüst and A. Schlüter.  Axialsymmetrische magnetohydrodynamische Gleichgewichtskonfigurationen.  Z. Naturforsch. A, 12:850–854, 1957.
  * [27] J.L. Luxon and B.B. Brown.  Magnetic analysis of non-circular cross-section tokamaks.  Nuclear Fusion, 22(6):813–821, 1982.
  * [28] P.J. Mc Carthy, P. Martin, and W. Schneider.  The CLISTE Interpretive Equilibrium Code.  Technical Report IPP Report 5/85, Max-Planck-Institut fur Plasmaphysik, 1999.
  * [29] A. Merle, R. Coelho, F. Carpanese, S. Dixon, M. Dunne, B. Faugeras, L. Fleury, J. Hollocombe, F. Imbeaux, L. Kogan, M. Romanelli, O. Sauter, W. Zwingmann, ASDEX-Upgrade team, JET contributors, MAST team, TCV team, and EUROfusion-IM team.  Equilibrium reconstruction of discharges from eurofusion tokamaks using the wpcd scientific workflows.  In European Physical Society, editor, 47th EPS Conference on Plasma Physics, volume 45A of Europhysics Conference Abstracts (ECA), page P2.1036, Milan, Italy, June 2021.
  * [30] J.-M. Moret, B.P. Duval, H.B. Le, S. Coda, F. Felici, and H. Reimerdes.  Tokamak equilibrium reconstruction code LIUQE and its real time implementation.  Fusion Eng. Design, 91(0):1–15, 2015.
  * [31] J. Nocedal and S. J. Wright.  Numerical optimization.  Springer Series in Operations Research and Financial Engineering. Springer, New York, second edition, 2006.
  * [32] F. Rapetti, B. Faugeras, and C. Boulbe.  High-order finite elements in tokamak free-boundary plasma equilibrium computations.  In ICOSAHOM 2020(2021) International Conference on Spectral and High Order Methods, Vienna, Austria, July 2021.
  * [33] R. Santos, R. Coelho, P. Rodrigues, B. Faugeras, H. Fernandes, B. B. Carvalho, D. Corona, H. Figueiredo, and H. Alves.  Plasma boundary reconstruction in ISTTOK using magnetic diagnostic data.  Journal of Instrumentation, 14(09):C09019–C09019, 2019.
  * [34] V D Shafranov.  Determination of the parameters $\beta_{p}$ and $l_{i}$ in a tokamak for arbitrary shape of plasma pinch cross-section.  Plasma Physics, 13(9):757, 1971.
  * [35] V.D. Shafranov.  On magnetohydrodynamical equilibrium configurations.  Soviet Journal of Experimental and Theoretical Physics, 6:545, 1958\.
  * [36] X. Song, X.M. Song, B. Li, J. Zhou, E. Nardon, H. Heumann, B. Faugeras, J.X. Li, Sh Wang, SH.Y. Liang, J.Z. Zhang, T.F Sun, W.B Li, Zh.H. Huang, L. Liu, Z.C. Yang, H.X. Wang, X.Q. Ji, W.L. Zhong, and HL-2M Team.  Plasma initiation and preliminary magnetic control in the HL-2M tokamak.  Nuclear Fusion, 61(8):086010, 2021.  https://doi.org/10.1088/1741-4326/ac09fc.
  * [37] D.W. Swain and G.H. Neilson.  An efficient technique for magnetic analysis for non-circular, high-beta tokamak equilibria.  Nuclear Fusion, 22(8):1015–1030, 1982.
  * [38] J. Wesson.  Tokamaks.  The International Series of Monographs in Physics. Oxford University Press, 2004.
  * [39] L.E. Zakharov and V.D. Shafranov.  Equilibrium of a toroidal plasma with noncircular cross-section.  Sov. Phys. Tech. Phys., 18(2):151–156, 1973.
  * [40] W. Zwingmann.  Equilibrium analysis of steady state tokamak discharges.  Nuclear Fusion, 43:842–850, 2003.


Weyl anomaly. The question that remains is whether any such polynomial, both
$\delta$-invariant and vanishing on the constraint surface, actually exists.

## Acknowledgments

C.I. is deeply grateful to M. Fröb for many patient and very helpful lessons
about FieldsX.

This work is supported in part by the Italian Istituto Nazionale di Fisica
Nucleare and by Research Projects, F.R.A. 2022 of the Università di Genova.
The work of A.W. is supported by the UKRI Frontier Research Grant,
underwriting the ERC Advanced Grant “Generalized Symmetries in Quantum Field
Theory and Quantum Gravity”.

##  Appendix A  $d=4$, $\mathcal{N}=1$ Lie superconformal algebra

In this appendix we review our conventions for the $d=4$, $\mathcal{N}=1$
superconformal algebra. The bosonic and fermionic generators, the
corresponding gauge fields and BRST ghosts are listed in Table 1.

Bosonic Symmetry | Generator | Gauge field |  | Ghost
---|---|---|---|---
Local Lorentz | $J_{ab}$ | spin connection | $\omega_{\mu}{}^{ab}$ | $\Omega^{ab}$
Weyl | $W$ | dilaton | $b_{\mu}$ | $\sigma$
$U(1)_{R}$ chiral R-symmetry | R | R-gauge field | $a_{\mu}$ | $\alpha$
Diffeomorphisms | $P_{a}$ | vielbein | $e_{\mu}{}^{a}$ | $\xi^{\mu}$
Special conformal | $K_{a}$ | conformal vielbein | $f_{\mu}{}^{a}$ | $\theta^{a}$
Fermionic Symmetry | Generator | Gauge field |  | Ghost
Supersymmetry | $Q_{\alpha}$ | gravitino | $\psi_{\mu}{}^{\alpha}$ | $\zeta^{\alpha}$
Conformal supersymmetry | $S_{\alpha}$ | conformal gravitino | $\tilde{\psi}_{\mu}{}^{\alpha}$ | $\eta^{\alpha}$
Table 1:  $\mathfrak{su}(2,2|1)$ symmetries and generators, with their
associated gauge fields and BRST ghosts.

The (anti)-commutation relations defining the $d=4,\;\mathcal{N}=1$
superconformal algebra are: 202020We take spinor contractions in the
$\searrow$ direction. Hence
$\lambda^{\alpha}\chi_{\alpha}=-\lambda_{\alpha}\chi^{\alpha}$. E.g.
$\zeta^{\alpha}(\Gamma^{a})_{\alpha\beta}\zeta^{\beta}=-\zeta^{\alpha}(\Gamma^{a})_{\alpha}{}^{\beta}\zeta_{\beta}=-\overline{\zeta}\,\Gamma^{a}\zeta$.

$[J_{ab},J_{cd}]=\eta_{ac}\,J_{db}-\eta_{bc}\,J_{da}+\eta_{bd}\,J_{ca}-\eta_{ad}\,J_{cb},\\\
$ $\displaystyle[J_{bc},P_{a}]$
$\displaystyle=\eta_{ac}\,P_{b}-\eta_{ab}\,P_{c},\qquad$
$\displaystyle[J_{bc},K_{a}]$
$\displaystyle=\eta_{ac}\,K_{b}-\eta_{ab}\,K_{c},$
$\displaystyle[P_{a},P_{b}]$ $\displaystyle=0,$ $\displaystyle[K_{a},K_{b}]$
$\displaystyle=0,$ $\displaystyle[W,P_{a}]$ $\displaystyle=P_{a},$
$\displaystyle[W,K_{a}]$ $\displaystyle=-K_{a},$
$[P_{a},K_{b}]=2\,(\eta_{ab}\,W+J_{ab}),\\\ $
$\displaystyle[J_{ab},Q_{\alpha}]$
$\displaystyle=\tfrac{1}{2}\,(\Gamma_{ab})_{\alpha}{}^{\beta}Q_{\beta},$
$\displaystyle[J_{ab},S_{\alpha}]$
$\displaystyle=\tfrac{1}{2}\,(\Gamma_{ab})_{\alpha}{}^{\beta}\,S_{\beta},$
$\displaystyle[W,Q]$ $\displaystyle=\tfrac{1}{2}\,Q,$ $\displaystyle[W,S]$
$\displaystyle=-\tfrac{1}{2}\,S,$ $\displaystyle[R,Q_{\alpha}]$
$\displaystyle=-\tfrac{3}{2}\,i\,(\Gamma_{5})_{\alpha}{}^{\beta}\,Q_{\beta},\qquad$
$\displaystyle[R,S_{\alpha}]$
$\displaystyle=\tfrac{3}{2}\,i\,(\Gamma_{5})_{\alpha}{}^{\beta}S_{\beta},$
$\displaystyle[P_{a},Q]$ $\displaystyle=0,$ $\displaystyle[K_{a},S]$
$\displaystyle=0,$ $\displaystyle[P_{a},S_{\alpha}]$
$\displaystyle=i\,(\Gamma_{a})_{\alpha}{}^{\beta}\,Q_{\beta},$
$\displaystyle[K_{a},Q_{\alpha}]$
$\displaystyle=-i\,(\Gamma_{a})_{\alpha}{}^{\beta}\,S_{\beta},$
$\displaystyle\\{Q_{\alpha},Q_{\beta}\\}$
$\displaystyle=-2\,(\Gamma^{a})_{\alpha\beta}\,P_{a},$
$\displaystyle\\{S_{\alpha},S_{\beta}\\}$
$\displaystyle=2\,(\Gamma^{a})_{\alpha\beta}\,K_{a},$
$\\{Q_{\alpha},S_{\beta}\\}=2\,i\,W\delta_{\alpha\beta}+2\,i\,(\Gamma^{ab})_{\alpha\beta}\,\tfrac{1}{2}\,J_{ab}+2\,(\Gamma_{5})_{\alpha\beta}\,R.$
(210)

If we collectively denote such generators by $\\{T_{i}\\}$ with $1\leq i\leq
24$, $T_{i}$ is bosonic for $1\leq i\leq 16$ and fermionic for $17\leq i\leq
24$. The _grading_ $|i|$ of $T_{i}$ is defined to be:

$|i|=\begin{cases}0\;\;(\text{mod}\;2),\hskip 14.22636pt\text{if $T_{i}$ is
bosonic}\hskip 14.22636pt(1\leq i\leq 16),\\\ 1\;\;(\text{mod}\;2),\hskip
14.22636pt\text{if $T_{i}$ is fermionic}\hskip 14.22636pt(17\leq i\leq
24).\end{cases}$ (211)

Ghosts are fields which have _opposite_ statistics, (i.e. $\mathbb{Z}_{2}$
gradings) with respect to the generator $T_{i}$ to which they correspond:

$|g^{i}|=|i|+1.$ (212)

The superLie bracket is written as

$[T_{j},T_{k}]=f_{j}{}^{i}{}_{k}\,T_{i},$ (213)

where $[\cdot,\cdot]$ denotes the commutator or anti-commutator

$[T_{i},T_{j}]\equiv T_{i}\,T_{j}-(-)^{|i||j|}T_{j}\,T_{i},$ (214)

Hence

$f_{j}{}^{i}{}_{k}=-(-)^{|j||k|}f_{k}{}^{i}{}_{j}$ (215)

The super-Jacobi equation is equivalent to the statement that
$f_{j}{}^{i}{}_{k}$ is a super-invariant tensor:

$(-)^{|i||k|}f_{i}{}^{l}{}_{j}\,f_{l}{}^{m}{}_{k}+(-)^{|j||i|}f_{j}{}^{l}{}_{k}\,f_{l}{}^{m}{}_{i}+(-)^{|k||j|}f_{k}{}^{l}{}_{i}\,f_{l}{}^{m}{}_{j}=0.$
(216)

Both ghosts $g^{i}$ and generators $T_{i}$ are graded, so that $g\equiv
g^{i}\,T_{i}$ is odd. Hence

$[g,g]=\tilde{f}_{j}{}^{i}{}_{k}\,g^{j}\,g^{k}\,T_{i},$ (217)

where the “Grassmann envelope structure constants” $\tilde{f}_{j}{}^{i}{}_{k}$
are related to the structure constants $f_{j}{}^{i}{}_{k}$ as follows

$\tilde{f}_{j}{}^{i}{}_{k}=(-)^{|j|(|k|+1)}f_{j}{}^{i}{}_{k}.$ (218)

Hence

$\tilde{f}_{k}{}^{i}{}_{j}=(-)^{(1+|j|)(1+|k|)}\,\tilde{f}_{j}{}^{i}{}_{k}.$
(219)

The structure constants of $\mathfrak{su}(2,2|1)$ are:

$f^{[ef]}_{[ab][cd]}=\eta_{[c[a}\,\delta^{[e}_{d]}\,\delta^{f]}_{b]},\\\ $
$\displaystyle f_{[bc],a}^{d}$ $\displaystyle=\eta_{a[c}\,\delta^{d}_{b]},$
$\displaystyle f_{[bc],\tilde{a}}^{\tilde{d}}$
$\displaystyle=\eta_{\tilde{a}[c}\,\delta^{\tilde{d}}_{b]},$ $\displaystyle
f^{b}_{W,a}$ $\displaystyle=\delta^{a}_{b},$ $\displaystyle
f^{\tilde{b}}_{W,\tilde{a}}$ $\displaystyle=-\delta^{\tilde{b}}_{\tilde{a}},$
$\displaystyle f^{[cd]}_{a,\tilde{b}}$
$\displaystyle=\delta^{[c}_{a}\delta^{d]}_{\tilde{b}},$ $\displaystyle
f^{W}_{a,\tilde{b}}$ $\displaystyle=2\,\eta_{a\tilde{b}},$ $\displaystyle
f^{\beta}_{[ab],\alpha}$
$\displaystyle=\tfrac{1}{2}\,(\Gamma_{ab})_{\alpha}{}^{\beta},$ $\displaystyle
f^{\tilde{\beta}}_{[ab],\tilde{\alpha}}$
$\displaystyle=\tfrac{1}{2}\,(\Gamma_{ab})_{\tilde{\alpha}}{}^{\tilde{\beta}},$
$\displaystyle f^{\alpha}_{W,\beta}$
$\displaystyle=\tfrac{1}{2}\,\delta^{\alpha}_{\beta},$ $\displaystyle
f^{\tilde{\alpha}}_{W,\tilde{\beta}}$
$\displaystyle=-\tfrac{1}{2}\,\delta^{\tilde{\alpha}}_{\tilde{\beta}},$
$\displaystyle f_{R,\alpha}^{\beta}$
$\displaystyle=-\tfrac{3}{2}\,i\,(\Gamma_{5})_{\alpha}{}^{\beta},\quad$
$\displaystyle f_{R,\tilde{\alpha}}^{\tilde{\beta}}$
$\displaystyle=\tfrac{3}{2}\,i\,(\Gamma_{5})_{\tilde{\alpha}}{}^{\tilde{\beta}},$
$\displaystyle f^{\beta}_{a,\tilde{\alpha}}$
$\displaystyle=i\,(\Gamma_{a})_{\tilde{\alpha}}{}^{\beta},$ $\displaystyle
f^{\tilde{\beta}}_{\tilde{a},{\alpha}}$
$\displaystyle=-i\,(\Gamma_{a})_{\alpha}{}^{\tilde{\beta}},$ $\displaystyle
f^{a}_{\alpha\beta}$ $\displaystyle=-2\,(\Gamma^{a})_{\alpha\beta},$
$\displaystyle f^{\tilde{a}}_{\tilde{\alpha}\tilde{\beta}}$
$\displaystyle=2\,(\Gamma^{\tilde{a}})_{\tilde{\alpha}\tilde{\beta}},$
$\displaystyle
f^{W}_{\alpha\tilde{\beta}}=2\,i\,\delta_{\alpha\tilde{\beta}},$
$\displaystyle
f^{ab}_{\alpha\tilde{\beta}}=2\,i\,(\Gamma^{ab})_{\alpha\tilde{\beta}},$
$\displaystyle
f^{R}_{\alpha\tilde{\beta}}=2\,(\Gamma_{5})_{\alpha\tilde{\beta}}.$ (220)

The mass dimensions of ghosts and gauge fields in $\mathcal{N}=1\;d=4$
conformal supergravity are fixed from their BRST transformations taking into
account that $s$ is dimensionless. For the standard bosonic YM symmetries for
which $s\,c=-c^{2}$, one has

$[\Omega^{ab}]=[\sigma]=[\alpha]=0.$ (221)

From the BRST transformation rules for the diffeormorphism ghost $\xi\mu$ we
obtain

$[\gamma^{\mu}]=[\xi^{\mu}]=-1.$ (222)

From the BRST transfomations for th supersymmetry ghosts we deduce

$[\zeta]=-\tfrac{1}{2},\qquad[\eta]=\tfrac{1}{2},\qquad[\theta]=1.$ (223)

For tensorial connections and curvatures we have

$[h_{\mu}]=[g]+1,\qquad[H_{\mu\nu}]=[h_{\mu}]+1.$ (224)

Note that that the physical fields $e_{\mu}{}^{a}$, $a_{\mu}$, $b_{\mu}$ and
$\psi_{\mu}{}^{\alpha}$ have canonical mass dimensions, $0,1,1,\tfrac{1}{2}$
respectively.212121Restoring the gravitational constant, which we put to 1,
graviton and gravitino would get the familiar mass dimensions, $1$ and
$\tfrac{3}{2}$. Instead the composite fields
$f_{\mu}{}^{a},\omega_{\mu}{}^{ab},\tilde{\psi}_{\mu}{}^{\alpha}$ have non-
canonical higher mass dimensions $2,1,\tfrac{3}{2}$. Table 2 lists the $W$ and
$R$ charges of the fields of the theory, that we took into account in Section
3 to construct the possible forms for the $\lambda_{0}^{i}$’s.

$\bm{h}$ | $W\text{-charge}$ | $R\text{-charge}$
---|---|---
$\quad e^{a}$ | $\quad\;\;\;1$ | $\quad\;\;\;0$
$\quad\bm{\omega}^{ab}$ | $\quad\;\;\;0$ | $\quad\;\;\;0$
$\quad\bm{b}$ | $\quad\;\;\;0$ | $\quad\;\;\;0$
$\quad\bm{a}$ | $\quad\;\;\;0$ | $\quad\;\;\;0$
$\quad\bm{f}^{a}$ | $\quad-1$ | $\quad\;\;\;0$
$\quad\bm{\psi}^{\alpha}$ | $\quad\;\;\;\tfrac{1}{2}$ | $\quad-\tfrac{3}{2}$
$\quad\bm{\tilde{\psi}}^{\alpha}$ | $\quad-\tfrac{1}{2}$ | $\quad\;\;\;\tfrac{3}{2}$
Table 2:  Weyl and R-charges of the connections.

##  Appendix B Supertrace for Lie superalgebras

In this appendix we review a few general properties of supertraces for Lie
superalgebras, relevant for superconformal anomalies. For more details see
KAC19778 , Berezin1987a .

If $V=V_{\bar{0}}\oplus V_{\bar{1}}$ is a $\mathbb{Z}_{2}$ graded vector
space, a _homogeneous basis_ of $V=V_{\bar{0}}\oplus V_{\bar{1}}$, where
$m=\dim{V_{\bar{0}}}$ and $n=\dim{V_{\bar{1}}}$ is of the form

$\\{e^{(b)}_{1},...e^{(b)}_{m},\;e^{(f)}_{m+1},...e^{(f)}_{m+n}\\},$ (225)

where the superscripts b and f stand for “bosonic” (even) and “fermionic”
(odd) respectively.

We can define linear _representations_ of Lie superalgebras, by associating to
each algebra generator an element of $End(V)=l(V)_{\bar{0}}\oplus
l(V)_{\bar{1}}$, which are matrices in a given basis, and by defining the
super-Lie bracket in terms of the graded-commutator of matrices.

The matrices representing bosonic $B\in l(V)_{\bar{0}}$ and fermionic $F\in
l(V)_{\bar{1}}$ operators respectively have the form

$B=\begin{pmatrix}B_{bb}&\vline&0\\\ \hline\cr
0&\vline&B_{ff}\end{pmatrix},\qquad F=\begin{pmatrix}0&\vline&F_{bf}\\\
\hline\cr F_{fb}&\vline&0\end{pmatrix}.$ (226)

Let $A\in l(V)$ be a generic operator. In block-diagonal form it is written as

$A=\begin{pmatrix}\alpha&\vline&\gamma\\\
\hline\cr\delta&\vline&\beta\end{pmatrix}.$ (227)

The _supertrace_ of $A$ is defined as

$\text{str}(A)=\tr(\alpha)-\tr(\beta).$ (228)

The supertrace is independent of the choice of (homogeneous) basis. It has has
the following properties

  * •

Consistency: $\text{str}(BF)=0=\text{str}(FB)$ where $B\in l(V)_{\bar{0}}$ and
$F\in l(V)_{\bar{1}}$.

  * •

Supersymmetry: $\text{str}(TA)=(-)^{|T||A|}\text{str}(AT)$ $\forall\,T,A\in
l(V)$.

  * •

Invariance: $\text{str}([T,A])=0$ $\forall\;T,A\in l(V)$ where $[\cdot,\cdot]$
denotes the super-Lie bracket.

If $\\{T_{i}\\}_{i\in I}$ are a linear representation $R$ of a super-Lie-
algebra we can define a superinvariant tensor as follows

$K^{R}_{ijk}=2\,\text{str}_{R}(T_{i}T_{j}T_{k})=\text{str}_{R}([T_{i},T_{j}]T_{k})+\text{str}_{R}(\\{T_{i},T_{j}\\}T_{k})=C(R)f_{ijk}+d_{ijk}(R).$
(229)

$K^{R}_{ijk}$ satisfies the following equation, consequence of the properties
of the supertrace:222222
$0=\text{str}([T_{l},T_{i}T_{j}T_{k}])=\text{str}([T_{l},T_{i}]T_{j}T_{k})+(-)^{|l||i|}\text{str}(T_{i}[T_{l},T_{j}]T_{k})+(-)^{|l|(|i|+|j|)}\text{str}(T_{i}T_{j}[T_{l},T_{k}]).$

$f_{m}{}^{l}{}_{i}\,K^{R}_{ljk}+(-)^{|l||i|}\,f_{m}{}^{l}{}_{j}\,K^{R}_{ilk}+(-)^{|l|(|i|+|j|)}\,f_{m}{}^{l}{}_{k}\,K^{R}_{ijl}=0.$
(230)

$C(R)$ is the _index_ of the representation, such that

$\text{str}_{R}(T_{i}T_{j})=C(R)g_{ij}$ (231)

and $g_{ij}$ is the Lie superalgebra Cartan-Killing metric, defined
as232323For $\mathfrak{su}(2,2|1)$ this metric is non-degenerate.

$g_{ij}=\text{str}_{\mathrm{super\text{-}adj}}(T_{i}T_{j}).$ (232)

$f_{ijk}$ are therefore related to the structure constants $f_{i}{}^{l}{}_{j}$
as follows:

$f_{ijk}=f_{i}{}^{l}{}_{j}\,g_{lk}.$ (233)

$f_{ijk}$ does not depend on the representation $R$ because the Cartan-Killing
metric is the _unique_ rank-two (super)-symmetric invariant tensor for a given
Lie-superalgebra. $f_{ijk}$ is completely “anti-symmetric” in the graded
sense, that is

$f_{jik}=-(-)^{|i||j|}f_{ijk},\qquad f_{ikj}=-(-)^{|k||j|}f_{ijk}.$ (234)

$d_{ijk}(R)$ is instead completely “symmetric” in the graded sense

$d_{jik}(R)=(-)^{|i||j|}d_{ijk}(R),\qquad d_{ikj}(R)=(-)^{|k||j|}d_{ijk}(R).$
(235)

The tensor $d_{ijk}(R)$ could in principle, for a generic superalgebra, depend
on the representation. However we computed solutions to the invariance
equation (230) for $d=4\,,\,\mathcal{N}=1$ superconformal algebra and we
obtained two linearly independent solutions: one which coincides (up to a
multiplicative constant) with the lowered structure constants $f_{ijk}$ and
another solution with precisely the symmetry properties of $d_{ijk}(R)$. Hence
there is a _unique_ rank three invariant tensor of $\mathfrak{su}(2,2|1)$ with
its symmetry properties, up to a multiplicative constant. Therefore,

$d_{ijk}(R)=2\,A(R)\,d_{ijk},$ (236)

where $d_{ijk}$ is independent of the representation $R$.

$d_{ijk}$ is related to the tensor $\tilde{d}_{ijk}$ which defines the
invariant polynomial $P_{3}(\bm{H})$, (see eq. (127)), as follows

$\\!\\!\text{str}_{R}(\bm{H}^{3})=A(R)\,(-)^{|i||j|+|i||k|+|j||k|}\,d_{ijk}\,\bm{H}^{i}\,\bm{H}^{j}\,\bm{H}^{k}\equiv
A(R)\,\tilde{d}_{ijk}\bm{H}^{i}\bm{H}^{j}\bm{H}^{k}=A(R)\,P_{3}(\bm{H}),$
(237)

that is

$\tilde{d}_{ijk}=(-)^{|i||j|+|i||k|+|j||k|}\,d_{ijk}.$ (238)

The sign factor relating $d_{ijk}$ to $\tilde{d}_{ijk}$, which is caused by
the fact that both the curvatures and generators are graded, is invariant
under exchange and cyclic permutation of its indices. Therefore does not
change the symmetry properties (235). It is important to keep in mind that the
“invariance” equation satisfied by the tensor $\tilde{d}_{ijk}$ — which
ensures BRST invariance of $P_{3}(\bm{H})$ — is different, although
equivalent, from eq. (230) valid for ${d}_{ijk}$:

$f_{m}{}^{l}{}_{i}\,\tilde{d}_{ljk}+(-)^{|j||m|}f_{m}{}^{l}{}_{j}\,\tilde{d}_{ilk}+(-)^{|m||k|+|m||j|}\,f_{m}{}^{l}{}_{k}\,\tilde{d}_{ijl}=0.$
(239)

The coefficient $A(R)$ in equations (236) is the _anomaly coefficient_ : it
describes the contribution to the superconformal anomaly of matter in a
representation $R$ of the superconformal algebra.

## References

  * (1) S.L. Adler, _Axial vector vertex in spinor electrodynamics_ , _Phys. Rev._ 177 (1969) 2426.
  * (2) J.S. Bell and R. Jackiw, _A PCAC puzzle: $\pi^{0}\to\gamma\gamma$ in the $\sigma$ model_, _Nuovo Cim. A_ 60 (1969) 47.
  * (3) C. Becchi, A. Rouet and R. Stora, _Renormalization of Gauge Theories_ , _Annals Phys._ 98 (1976) 287.
  * (4) R. Stora, _Continuum gauge theories_ , in _New developments in quantum field theory and statistical mechanics, Cargese 1976, Ed. NATO ASI Ser B vol. 26_ , P.M. M. Levy, ed., pp. 201–224, Plenum Press (1977), DOI.
  * (5) F. Langouche, T. Schucker and R. Stora, _Gravitational Anomalies of the Adler-bardeen Type_ , _Phys. Lett. B_ 145 (1984) 342.
  * (6) M.B. Green and J.H. Schwarz, _Anomaly Cancellation in Supersymmetric D=10 Gauge Theory and Superstring Theory_ , _Phys. Lett. B_ 149 (1984) 117.
  * (7) J.M. Maldacena, _The Large N limit of superconformal field theories and supergravity_ , _Adv. Theor. Math. Phys._ 2 (1998) 231 [hep-th/9711200].
  * (8) E. Witten, _Anti-de Sitter space and holography_ , _Adv. Theor. Math. Phys._ 2 (1998) 253 [hep-th/9802150].
  * (9) O. Piguet and K. Sibold, _Algebraic characterization of superconformal anomalies_ , _Nucl. Phys. B_ 197 (1982) 247.
  * (10) L. Bonora, P. Cotta-Ramusino and C. Reina, _Conformal Anomaly and Cohomology_ , _Phys. Lett. B_ 126 (1983) 305.
  * (11) L. Bonora, P. Pasti and M. Tonin, _Gravitational and Weyl Anomalies_ , _Phys. Lett. B_ 149 (1984) 346.
  * (12) L. Bonora, P. Pasti and M. Tonin, _Cohomologies and Anomalies in Supersymmetric Theories_ , _Nucl. Phys. B_ 252 (1985) 458.
  * (13) L. Bonora, P. Pasti and M. Bregola, _Weyl Cocycles_ , _Class. Quant. Grav._ 3 (1986) 635.
  * (14) L. Bonora, P. Pasti and M. Tonin, _The Chiral Anomaly in Supersymmetric Gauge Theories Coupled to Supergravity_ , _Phys. Lett. B_ 167 (1986) 191.
  * (15) S. Deser and A. Schwimmer, _Geometric classification of conformal anomalies in arbitrary dimensions_ , _Phys. Lett. B_ 309 (1993) 279 [hep-th/9302047].
  * (16) F. Brandt, _Local BRST cohomology in minimal D = 4, N=1 supergravity_ , _Annals Phys._ 259 (1997) 357 [hep-th/9609192].
  * (17) F. Brandt, _Extended BRST cohomology, consistent deformations and anomalies of four-dimensional supersymmetric gauge theories_ , _JHEP_ 04 (2003) 035 [hep-th/0212070].
  * (18) N. Boulanger, _Algebraic Classification of Weyl Anomalies in Arbitrary Dimensions_ , _Phys. Rev. Lett._ 98 (2007) 261302 [0706.0340].
  * (19) L. Bonora and S. Giaccari, _Weyl transformations and trace anomalies in N=1, D=4 supergravities_ , _JHEP_ 08 (2013) 116 [1305.7116].
  * (20) R. Stora, _Algebraic structure and topological origin of anomalies_ , in _Progress in Gauge Field Theory_ , G. ’t Hooft, A. Jaffe, H. Lehmann, P.K. Mitter, I.M. Singer and R. Stora, eds., (Boston, MA), pp. 543–562, Springer US (1984), DOI.
  * (21) B. Zumino, _Chiral Anomalies and Differential Geometry: Lectures given at Les Houches, August 1983_ , in _Les Houches Summer School on Theoretical Physics: Relativity, Groups and Topology_ , pp. 1291–1322, 10, 1983\.
  * (22) J. Manes, R. Stora and B. Zumino, _Algebraic Study of Chiral Anomalies_ , _Commun. Math. Phys._ 102 (1985) 157.
  * (23) W.A. Bardeen and B. Zumino, _Consistent and Covariant Anomalies in Gauge and Gravitational Theories_ , _Nucl. Phys. B_ 244 (1984) 421.
  * (24) F. Brandt, N. Dragon and M. Kreuzer, _All consistent Yang-Mills anomalies_ , _Phys. Lett. B_ 231 (1989) 263.
  * (25) M. Kaku, P.K. Townsend and P. van Nieuwenhuizen, _Superconformal Unified Field Theory_ , _Phys. Rev. Lett._ 39 (1977) 1109.
  * (26) M. Kaku, P.K. Townsend and P. van Nieuwenhuizen, _Properties of Conformal Supergravity_ , _Phys. Rev. D_ 17 (1978) 3179.
  * (27) P. van Nieuwenhuizen, _Supergravity as a Yang–Mills Theory_ , in _50 years of Yang-Mills theory_ , G. ’t Hooft, ed., pp. 433–456 (2005), DOI [hep-th/0408137].
  * (28) P. Van Nieuwenhuizen, _Supergravity_ , _Phys. Rept._ 68 (1981) 189.
  * (29) E.S. Fradkin and A.A. Tseytlin, _Conformal Supergravity_ , _Phys. Rept._ 119 (1985) 233.
  * (30) L. Baulieu and M.P. Bellon, _$p$ Forms and Supergravity: Gauge Symmetries in Curved Space_, _Nucl. Phys._ B266 (1986) 75.
  * (31) L. Baulieu, M.P. Bellon and S. Ouvry, _Off-shell BRS invariant construction of N=1 D=4 conformal supergravity_ , _Phys. Lett. B_ 171 (1986) 396.
  * (32) C. Imbimbo and D. Rosa, _The topological structure of supergravity: an application to supersymmetric localization_ , _JHEP_ 05 (2018) 112 [1801.04940].
  * (33) M.B. Fröb, C. Imbimbo and N. Risso, _Deformations of supergravity and supersymmetry anomalies_ , _JHEP_ 12 (2021) 009 [2107.03401].
  * (34) M.B. Fröb, _FieldsX – An extension package for the xAct tensor computer algebra suite to include fermions, gauge fields and BRST cohomology_ , 2008.12422.
  * (35) I. Papadimitriou, _Supersymmetry anomalies in $\mathcal{N}=1$ conformal supergravity_, _JHEP_ 04 (2019) 040 [1902.06717].
  * (36) V. Kac, _Lie superalgebras_ , _Advances in Mathematics_ 26 (1977) 8.
  * (37) F.A. Berezin, _Lie superalgebras_ , in _Introduction to Superanalysis_ , A.A. Kirillov, ed., (Dordrecht), pp. 231–244, Springer Netherlands (1987), DOI.


# Efficient and Accurate Learning of Mixtures of Plackett-Luce Models

Duc Nguyen1 Anderson Y. Zhang2

###### Abstract

Mixture models of Plackett-Luce (PL) – one of the most fundamental ranking
models – are an active research area of both theoretical and practical
significance. Most previously proposed parameter estimation algorithms
instantiate the EM algorithm, often with random initialization. However, such
an initialization scheme may not yield a good initial estimate and the
algorithms require multiple restarts, incurring a large time complexity. As
for the EM procedure, while the E-step can be performed efficiently,
maximizing the log-likelihood in the M-step is difficult due to the
combinatorial nature of the PL likelihood function (Gormley and Murphy 2008).
Therefore, previous authors favor algorithms that maximize surrogate
likelihood functions (Zhao et al. 2018, 2020). However, the final estimate may
deviate from the true maximum likelihood estimate as a consequence. In this
paper, we address these known limitations. We propose an initialization
algorithm that can provide a provably accurate initial estimate and an EM
algorithm that maximizes the true log-likelihood function efficiently.
Experiments on both synthetic and real datasets show that our algorithm is
competitive in terms of accuracy and speed to baseline algorithms, especially
on datasets with a large number of items.

##

###  1 Introduction

Learning to rank is an active area of research with wide-ranging applications
in recommendation systems, information retrieval, crowdsourcing and the social
sciences. The Plackett-Luce (PL) model (Plackett 1975; Luce 1959) is one of
the most fundamental ranking models. In a universe of $n$ items, the PL model
posits that item $i$ has a latent _utility_ $\theta^{*}_{i}\in\mathbb{R}$. The
probability of observing a full ranking $\pi$ given by the user (most
preferred item first) is given as

$\mathbb{P}^{PL}(\pi=[\pi_{1},..,\pi_{n}]\,\lvert\theta^{*})=\prod_{i=1}^{n-1}\frac{\mathrm{exp}\left(\theta^{*}_{\pi_{i}}\right)}{\sum_{j=i}^{n}\mathrm{exp}\left(\theta^{*}_{\pi_{j}}\right)}\,.$
(1)

The maximum likelihood estimate (MLE) can be obtained using iterative
algorithms such as the Minorize-Maximize (MM) algorithm of Hunter (2004) and
enjoys favorable theoretical properties (Hajek, Oh, and Xu 2014). In recent
years, an algorithm known as Luce spectral ranking (LSR) (Maystre and
Grossglauser 2015) has become the method of choice for maximum likelihood
inference for PL models. LSR outputs the MLE just like MM but is often much
faster.

The PL model is closely connected to the Bradley-Terry-Luce (BTL) model (Luce
1959) for _pairwise comparisons_. For two items $i\neq j$, the probability
that $i$ is ranked ahead of $j$ _in a ranking_ is equal to the probability
that $i$ beats $j$ in a _pairwise comparison_ under the BTL model. That is,

$\mathbb{P}^{PL}(\pi(i)<\pi(j))=\mathbb{P}^{BTL}_{ij}=\frac{1}{1+\mathrm{exp}\left(-(\theta^{*}_{i}-\theta^{*}_{j})\right)}\,,$
(2)

where $\pi(i)$ is the position of item $i$ in ranking $\pi$.

The classical PL model assumes that there is a universal preference ordering
of the items according to their utilities. However, in practice, there might
be multiple subpopulations of users with different preference profiles which
cannot be fully captured by a single PL model. In such settings, a mixture of
PL models is a more appropriate modeling assumption.

Problem Descriptions. Consider a mixture model with $K$ components and $n$
items for some constant $K$. Let
$\beta^{*}=[\beta^{*}_{1},\ldots,\beta^{*}_{K}]^{\top}$,
${\beta^{*}}^{\top}\mathbf{1}=1$ denote the mixing distribution. For component
$k\in[K]$ (where $[a]$ denotes $[1,\ldots,a]$), the utility parameters for the
items are

$\theta^{*k}=[\theta^{*k}_{1},\ldots,\theta^{*k}_{n}]^{\top}\in\mathbb{R}^{n}\,.$

Let $\bm{\theta}^{*}=[\theta^{*1},\ldots,\theta^{*K}]\in\mathbb{R}^{n\times
K}$ denote the concatenation of the $K$ sets of parameters. A _ranking
dataset_ $\Pi$ is a collection of full rankings.

Consider the following generative model for a ranking dataset of size $m$. For
$l\in[m]$, let ${z^{*}_{l}}\in[K]$ denote the mixture component membership
where $\mathbb{P}({z^{*}_{l}}=k)=\beta^{*}_{k}$. Then a permutation $\pi_{l}$
is drawn from the PL distribution parametrized by $\theta^{{*z_{l}}}$. That
is,

$\mathbb{P}^{PL}(\pi_{l}=[\pi_{l,1},\ldots,\pi_{l,n}]\,\lvert\,{z^{*}_{l}},\bm{\theta}^{*})=\prod_{i=1}^{n-1}\frac{\mathrm{exp}\left(\theta^{*z^{*}_{l}}_{\pi_{l,i}}\right)}{\sum_{j=i}^{n}\mathrm{exp}\left(\theta^{*z_{l}^{*}}_{\pi_{l,j}}\right)}\,,$
(3)

where $\pi_{l,i}$ denote the $i$-th item in permutation $\pi_{l}$. The reader
may recognize two identifiability issues here. The first is parameter
translation. For each component, the distributions parametrized by
$\theta^{*k}$ and $\theta^{*k}+c\cdot\mathbf{1}_{n}$ are the same for any
$c\in\mathbb{R}$. The second is mixture components (columns of
$\bm{\theta}^{*}$) relabeling. To account for these issues, we consider the
following error metric.

$\text{dist}(\bm{\theta},\bm{\theta}^{*}):=\min_{R\in\mathcal{O}^{K\times
K}}\lVert N(\bm{\theta})R-N(\bm{\theta}^{*})\rVert_{F}\,,$ (4)

where $\mathcal{O}^{K\times K}$ is the set of all permutation matrices (Strang
et al. 1993, Chapter 2) of size $K\times K$ and $N$ is the normalization
operator (i.e.,
$N(\bm{\theta})_{\cdot\,k}=\theta^{k}-\frac{1}{n}(\theta^{k})^{\top}\mathbf{1}_{n}$).

Prior Works. Generalizing the PL model to mixtures adds a layer of complexity
to the inference problem. In general, the likelihood function is non-convex in
the model parameters. Most previously proposed algorithms instantiate the EM
algorithm (Dempster, Laird, and Rubin 1977). As a general recipe, an EM
algorithm is initialized with some parameter $\bm{\theta}^{(0)}$ (e.g., using
random initialization). It then repeats the following two steps for
$t=1,2,\ldots$ until convergence.

The E-step computes the posterior class probability conditioned on the current
estimate:

$q_{l}^{k}=\mathbb{P}({z^{*}_{l}}=k\,\lvert\,\pi_{l},\bm{\theta}^{(t-1)})\propto\beta_{k}\cdot\mathbb{P}^{PL}(\pi_{l}\,\lvert{z^{*}_{l}}=k,\bm{\theta}^{(t-1)})\,$
(5)

for $l\in[m],k\in[K]$ where $\mathbb{P}^{PL}$ is given in Equation (3) and
$\beta$ the prior class probability. Thanks to the closed form of the PL
likelihood function, the E-step can be done efficiently. The M-step obtains
the next estimate $\bm{\theta}^{(t)}$ by maximizing the joint log-likelihood
function which decomposes into $K$ _weighted log-likelihood functions_.
Namely,

$\bm{\theta}^{(t)}=\arg\max_{\bm{\theta}}\sum_{k=1}^{K}\bigg{(}\sum_{l=1}^{m}q_{l}^{k}\log\mathbb{P}^{PL}(\pi_{l},{z^{*}_{l}}=k\,\lvert\,\bm{\theta})\bigg{)}\,.$
(6)

Due to the combinatorial nature of the PL likelihood function, the derivative
of the log likelihoood function has a complicated form. As a result,
maximizing the (weighted) log-likelihood via gradient-based algorithms quickly
becomes inefficient as $n$ grows.

The first practical approach towards solving the M-step uses the Minorize-
Maximize algorithm of Hunter (2004), yielding the EMM algorithm of Gormley and
Murphy (2008). While guaranteed to solve the M-step, it has been observed that
the MM subroutine converges slowly even for datasets with a moderate number of
items (e.g., Figure 2). Motivated by practical concerns, researchers have
developed _pseudo-likelihood_ estimators that optimize, instead of the true
log-likelihood function, _alternative objective functions_. Two such
algorithms are the Generalized Method of Moments (GMM) of Azari Soufiani et
al. (2013) and Composite Marginal Likelihood (CML) of Zhao and Xia (2018). It
has been observed experimentally that GMM is considerably faster than MM and
CML is even faster than GMM with comparable accuracy. Besides maximum
likelihood (ML) inference methods, previous authors have also proposed
Bayesian inference algorithms (Guiver and Snelson 2009; Mollica and Tardella
2017). In this paper, we focus primarily on ML algorithms but include
additional experiments with Bayesian methods in the supplementary materials.

Using GMM and CML to solve the M-step gives us the EM-GMM algorithm (Zhao,
Villamil, and Xia 2018) and the EM-CML algorithm (Zhao, Liu, and Xia 2020),
respectively. The only non-EM algorithm for learning PL mixtures that we are
aware of is a GMM-based algorithm proposed in Zhao, Piech, and Xia (2016);
Zhao and Xia (2019). However, the construction of the algorithm is quite ad-
hoc and the authors did not show extension of the algorithm to more than 2
mixture components. In addition, previous authors primarily restrict their
experiments to datasets with a small number of items such as the SUSHI
datasets (Kamishima 2003) with $n=10$. It is unknown how the previous methods
perform when $n$ is large. Recent works have also studied PL mixtures learning
with features and partial rankings (Tkachenko and Lauw 2016; Liu et al. 2019).
While we include possible extensions of our algorithm in the supplementary
materials, our main focus in this paper is an _improved algorithm for the
classical setting_.

Our Contributions. We propose a new EM algorithm for learning mixtures of PL
models that

  * •

Has a provably accurate initialization procedure with a finite sample error
guarantee, the first of its kind in the literature;

  * •

Efficiently maximizes the weighted log-likelihood function in the M-step
without using a surrogate likelihood or objective function, thus returning the
true maximum likelihood estimate;

  * •

Performs competitively with the previously proposed algorithms in terms of
accuracy and speed, and is scalable to datasets with $n\geq 100$.

###  2 The Spectral EM Algorithm

In this section, we present our algorithmic contributions. Section 2.1
describes the spectral initialization algorithm and Section 2.2 describes the
EM refinement procedure.

####  2.1 Spectral Initialization

The initialization for our algorithm is delegated to spectral clustering
(Algorithm 1) and a least squares minimization algorithm (Algorithm 2). To
apply spectral clustering, we first embed each ranking $\pi_{l}$ into a
‘pairwise vector’ – $X_{l}\in\\{0,1\\}^{{n\choose 2}}$ where each entry
corresponds to a pair of items. As an overload of notation, we use
$d=(d_{1},d_{2})$ where $d_{1}<d_{2}$ to denote the entry corresponding to the
pair $(d_{1},d_{2})$. Define

$X_{l,d}(\pi)=\begin{cases}1&\text{ if }\pi_{l}(d_{1})<\pi_{l}(d_{2})\\\
0&\text{ otherwise }\end{cases}\,.$ (7)

Let $X\in\mathbb{R}^{m\times{n\choose 2}}$ denote the concatenation of the
embeddings of $m$ rankings in dataset $\Pi$. Given a target number of
components $K$, Algorithm 1 can then be applied to the rows of $X$ to obtain
$K$ clusters, $\\{\hat{C}^{k}\\}_{k=1}^{K}\subseteq[m]$.

For each cluster of rankings $\hat{C}^{k}$, we estimate the preference
probability for a pair $(i,j)$ as

$\hat{P}^{k}_{ij}=\frac{1}{\lvert\hat{C}^{k}\rvert}\sum_{l\in\hat{C}^{k}}\mathbf{1}[\pi_{l}(i)<\pi_{l}(j)]\,.$
(8)

From the preference probability estimates for all pairs, Algorithm 2 recovers
the utility parameter $\hat{\theta}^{k}$. It applies the logit function on the
pairwise probabilities and solves a constrained least squares minimization
problem, which can be efficiently done using off-the-shelf solvers (Virtanen
et al. 2020). Algorithm 3 summarizes the spectral initialization algorithm.

Algorithm 1 Spectral Clustering with Adaptive Dimension Reduction

Input:  Dataset $\Pi=\\{\pi_{1},\ldots,\pi_{m}\\}$, number of mixture
components $K$ and threshold $T$.
Output:  $K$ clusters of rankings.

1: Embed the rankings as the rows of a matrix $X\in\\{0,1\\}^{m\times{n\choose
2}}$ according to Equation (7).

2: Perform SVD: $X=USV^{\top}$, where the singular values are arranged from
largest to smallest.

3: Let $\hat{r}$ be largest index in $[K]$ such that the difference between
the successive singular values is greater than $T$, i.e.,
$\hat{r}=\max\\{a\in[K]\,:\,S_{aa}-S_{(a+1)(a+1)}\geq T\\}\,.$

4: Run k-means on the rows of $XV_{1:\hat{r}}$ with $K$ clusters:
$\big{(}\hat{z},\\{\hat{c}_{k}\\}_{k=1}^{K}\big{)}=\underset{\begin{subarray}{c}z\in\\{1,\ldots,K\\}^{m}\\\
\\{c_{k}\\}\in\mathbb{R}^{\hat{r}}\end{subarray}}{\arg\min}\sum_{l=1}^{m}\lVert
V_{1:\hat{r}}^{\top}X_{l}-c_{z_{l}}\rVert_{2}^{2}\,.$

5: Return clusters $\hat{C}^{k}=\\{l\in[m]\,:\,\hat{z}_{l}=k\\}$ for
$k\in[K]$.

Algorithm 2 Least Squares Parameter Estimation

Input: Pairwise preference matrix $\hat{P}\in\mathbb{R}^{n\times n}$.
Output: Normalized parameter estimate $\hat{\theta}$.

1: Solve the least squares optimization problem
$\hat{\theta}=\arg\max_{\theta:\theta^{\top}\mathbf{1}=0}\sum_{i\neq
j}(\hat{\phi}_{ij}-(\theta_{i}-\theta_{j}))^{2}\,,$ where
$\hat{\phi}_{ij}=\ln(\hat{P}_{ij}/(1-\hat{P}_{ij}))\,.$

Algorithm 3 Spectral Initialization

Input:  Dataset $\Pi=\\{\pi_{1},\ldots,\pi_{m}\\}$, number of mixture
components $K$.
Output:  Parameter estimates for $K$ mixture components
$\hat{\bm{\theta}}=[\hat{\theta}^{1},\ldots,\hat{\theta}^{K}]\in\mathbb{R}^{n\times
K}$.

1: Run Algorithm 1 on $\Pi$ with $T=\sqrt{n}\sqrt{m+n}\sqrt{\log n}$ to obtain
$K$ clusters $\hat{C}^{1},\ldots,\hat{C}^{K}$.

2: Estimate the pairwise preference probabilities $\hat{P}^{k}_{ij}$ per
Equation (8) for each cluster.

3: Run Algorithm 2 on $\\{\hat{P}^{k}\\}_{k=1}^{K}$ and return the parameter
estimates for $K$ mixture components.

Remarks. The application of spectral clustering to mixtures of PL models has
also appeared in a manuscript by Shah and Song (2018). There, the authors
apply the classical spectral clustering algorithm – clustering the rows of
$XV_{1:K}$ – and their analysis requires a spectral gap condition which is
hard to verify. We use spectral clustering with adaptive dimension reduction
and our analysis does not require any spectral gap condition (Zhang and Zhou
2022). Furthermore, we focus on _parameter estimation_ while Shah and Song
(2018) only focus on clustering, resulting in different theoretical
guarantees. The choice of threshold $T$ in Algorithm 3 is to satisfy a mild
technical condition in the analysis of spectral clustering. In our
experiments, the performance of the EM algorithm does not seem to critically
depend on this threshold.

Intuition behind Algorithm 2. Recall the connection between the PL model and
the BTL model in Equation (2). Suppose we observe a large sample drawn from a
single PL distribution. Then $\hat{P}_{ij}\approx
P_{ij}^{*}=1/\big{(}\mathrm{exp}\left(-(\theta_{i}^{*}-\theta_{j}^{*})\right)\big{)}$
and
$\hat{\phi}_{ij}=\ln(\hat{P}_{ij}/1-\hat{P}_{ij})\approx\theta^{*}_{i}-\theta^{*}_{j}$.
Solving the least squares optimization problem recovers
$\hat{\theta}\approx\theta^{*}$. In the mixture setting, if the estimates
$\hat{P}^{k}$’s are accurate, we obtain good parameter estimates (e.g.,
Theorem 3.1). Rajkumar and Agarwal (2016) apply a similar idea in their
algorithm for _ranking from comparisons of $O(n\log n$) pairs under a single
BTL model_. They first apply the logit function on the pairwise preference
probabilities, followed by a low rank matrix completion algorithm (Keshavan,
Montanari, and Oh 2009). Their algorithm _produces a ranking_. On the other
hand, our goal is mixture learning and the resulting theoretical analysis is
different.

####  2.2 Iterative Refinement via EM

The Weighted LSR Algorithm. As noted before, we wish to maximize the weighted
log-likelihood (6) efficiently. Towards this goal, we generalize the Luce
spectral ranking (LSR) algorithm (Maystre and Grossglauser 2015) to
incorporates sample weights. The original LSR algorithm produces the MLE. Our
generalized algorithm outputs the weighted MLE (see Theorem 3.2).

The intuition behind LSR is an interpretation of the _PL ranking generative
process as a sequence of choices_ (Plackett 1975). Given a ranking $\pi_{l}$,
define its _choice breaking_ as

$\mathcal{B}_{\pi_{l}}=\\{(\pi_{l,1},\\{\pi_{l,1},\ldots,\pi_{l,n}\\},l),\ldots,(\pi_{l,n-1},\\{\pi_{l,n-1},\pi_{l,n}\\},l)\\}\,.$

Each tuple $(i,A,l)\in\mathcal{B}(\pi_{l})$ is a _choice enumeration_ of the
ranking $\pi_{l}$. Given a ranking dataset $\Pi=\\{\pi_{1},\ldots,\pi_{m}\\}$,
define the _choice breaking_ of $\Pi$ as the union of all ranking-level choice
breakings:

$\mathcal{B}_{\Pi}=\mathcal{B}_{\pi_{1}}\cup\ldots\cup\mathcal{B}_{\pi_{m}}\,.$
(9)

Note that $\lvert\mathcal{B}_{\Pi}\rvert=m(n-1)$. When the dataset $\Pi$ is
clear from context, we simply use $\mathcal{B}$ to denote the dataset-level
choice breaking.

_We now introduce sample weights_. Firstly, define the ‘weight’ of a choice
breaking $\mathcal{B}$ with weight vector $w$ and parameter
$\theta\in\mathbb{R}^{n}$ as

$\gamma(\mathcal{B},w,\theta)=w^{\top}\bigg{(}\frac{1}{\sum_{j\in
A}e^{\theta_{j}}}\bigg{)}_{(i,A,l)\in\mathcal{B}}\,,$ (10)

where $w\in\mathbb{R}^{m(n-1)}_{+}$ is an arbitrary weight vector;
$\bigg{(}\frac{1}{\sum_{j\in
A}e^{\theta_{j}}}\bigg{)}_{(i,A,l)\in\mathcal{B}}$ is also vector in
$\mathbb{R}^{m(n-1)}_{+}$ where each entry corresponds to a choice enumeration
$(i,A,l)$.

The reader may recognize that the weight vector $w$ has the same size as the
choice breaking while sample weights are often given at the ranking level –
each ranking $\pi_{l}$ is assigned a weight $q_{l}$ for $l\in[m]$ as in (6).
Given sample weights $q=(q_{1},\ldots,q_{m})$, one simply sets

$w=[\underbrace{q_{1},\ldots,q_{1}}_{n-1\text{
terms}},\underbrace{q_{2},\ldots,q_{2}}_{n-1\text{
terms}},\ldots,\underbrace{q_{m},\ldots,q_{m}}_{n-1\text{ terms}}]^{\top}\,.$
(11)

Given a choice breaking $\mathcal{B}$ and items $i,j$, define the set of
choice enumerations where $i$ ‘beats’ $j$ as

$\mathcal{B}_{i\succ j}=\\{(i,A,l)\in\mathcal{B}\,:\,j\in A\\}\,.$

As a shorthand notation, for a weight vector $w$ corresponding to choice
breaking $\mathcal{B}$, define $w_{j\succ i}$ as the sub-vector of $w$
corresponding to $\mathcal{B}_{i\succ j}$.

Similarly to the original LSR algorithm, we construct a Markov chain (MC) and
recover PL parameters from its stationary distribution. This MC has $n$
states. Given choice breaking $\mathcal{B}$, weight vector $w$ and parameter
$\theta$, the pairwise transition probabilities of $M$ are given as

$M_{ij}=\begin{cases}\frac{1}{d}\cdot\gamma(\mathcal{B}_{j\succ i},{w_{j\succ
i}},\theta)&\text{if }i\neq j\\\ 1-\frac{1}{d}\cdot\sum_{k\neq
i}\gamma(\mathcal{B}_{k\succ i},w_{k\succ i},\theta)&\text{if
}i=j\end{cases}\,,$ (12)

where $d$ is a sufficiently large normalization constant such that $M$ does
not contain any negative entries. Intuitively, $M_{ij}$ is proportional to the
sum of the weights of all choice enumerations where $j$ ‘beats’ $i$.

Algorithm 4 summarizes the weighted LSR algorithm. It repeatedly constructs a
Markov chain based on the current estimate, computes its stationary
distribution and recovers the next estimate until convergence. When sample
weights are not given, the weighted LSR algorithm reduces to the original LSR
algorithm.

Algorithm 4 Weighted Luce Spectral Ranking

Input:  Dataset $\Pi=\\{\pi_{1},\ldots,\pi_{m}\\}$, (optional) weight vector
$q\in\mathbb{R}^{m}_{+}$ and (optional) initial estimate
$\hat{\theta}^{(0)}\in\mathbb{R}^{n}$.
Output:  Normalized estimate of the item parameters
$\hat{\theta}\in\mathbb{R}^{n}$.

1: Obtain choice breaking $\mathcal{B}$ from $\Pi$ per Equation (9).

2: If the weight vector $q$ is not given, set $q=\mathbf{1}_{m}$.

3: Construct $w$ from $q$ per Equation (11).

4: If the initial estimate is not given, set
$\hat{\theta}^{(0)}=\mathbf{0}_{n}$.

5: For $t=1,\ldots$ until convergence 5.1: Construct a Markov chain $M$ with
pairwise transition probability per Equation (12) from choice breaking
$\mathcal{B}$, weight vector $w$ and parameter $\hat{\theta}^{(t-1)}$. 5.2:
Compute the stationary distribution of $M$ (e.g., via power iteration), $p$
and return the normalized estimate
$\hat{\theta}^{(t)}=\log(p)-\big{(}\frac{1}{n}\sum_{i=1}^{n}\log(p)\big{)}\cdot\mathbf{1}_{n}$

The EM-LSR Algorithm. In the E-step, we compute the posterior class
probabilities $q^{k}\in\mathbb{R}^{m},k\in[K]$. The M step consists of $K$
maximization problem as shown in Equation (6). These can be solved in parallel
by running Algorithm 4 on $\Pi$ using $q^{k}$ as sample weights for $k\in[K]$.
Algorithm 5 summarizes the overall algorithm.

Algorithm 5 Spectral EM (EM-LSR)

Input:  Dataset $\Pi=\\{\pi_{1},\ldots,\pi_{m}\\}$, number of components $K$,
prior distribution $\beta$, (optional) initial estimate
$\hat{\bm{\theta}}^{(0)}\in\mathbb{R}^{n\times K}$.
Output:  Normalized estimate
$\hat{\bm{\theta}}=[\hat{\theta}^{1},\ldots,\hat{\theta}^{K}]$.

1: If $\hat{\bm{\theta}}^{(0)}$ is not given, run Algorithm 3 on $\Pi$ with
$K$ mixture components and set $\hat{\bm{\theta}}^{(0)}$ to the output.

2: For $t=1,2,\ldots$ until convergence 2.1: E-step – Compute the class
posterior probabilities
$q^{k}_{l}=p(z^{*}_{l}=k|\pi_{l},\hat{\bm{\theta}}^{(t-1)})$ for
$l\in[m],k\in[K]$. 2.2: M-step – Estimate $\hat{\theta}^{k(t)}$ by running
Algorithm 4 on $\Pi$ with sample weight vector
$q^{k}=[q^{k}_{1},\ldots,q^{k}_{m}]$ and initial estimate
$\hat{\theta}^{k(t-1)}$ for $k\in[K]$.

In another EM-based approach for learning PL mixtures, Liu et al. (2019) use
the _unweighted LSR_ algorithm. There, the E-step remains the same. The key
differences lie in initialization (they use random initialization) and in the
M-step. Our algorithm maximizes the weighted log-likelihood via weighted LSR
and is therefore an exact EM algorithm. On the other hand, Liu et al. use the
posterior class probabilities to perform a random clustering of the rankings
and then run unweighted LSR on each cluster, making their algorithm an
_inexact_ EM algorithm. From additional experiments in the supplementary
materials, one can observe that the stochastic M-step actually leads to _worse
estimates without a significant reduction in inference time_.

###  3 Theoretical Analysis

In this section, we study the theoretical properties of EM-LSR. Section 3.1
presents the finite sample error guarantee for the spectral initialization
algorithm. Section 3.2 focuses on the analysis of the M-step.

####  3.1 Spectral Initialization

Central to the analysis of the spectral initialization algorithm is the
accuracy of spectral clustering (Algorithm 1). Our analysis starts from the
fact that, under the pairwise representation in Equation (7), the PL
distribution exhibits _sub-gaussian characteristics_ (Vershynin 2018; Shah and
Song 2018). The detailed descriptions of these characteristics are not
immediately important to our discussions so we refer the interested reader to
the supplementary materials. However, we emphasize that these characteristics
also appear in a broad class of ranking models known as _random utility
models_ (RUMs) that subsume the PL model. The spectral clustering algorithm is
model-agnostic. It can be applied to mixtures of sub-gaussian distributions
and enjoys high clustering accuracy if the signal-to-noise ratio (SNR) is
high. We also show how, by changing the mapping function used in Algorithm 2,
we can perform parameter estimation for a general RUM, not just PL. Thanks to
this flexibility, Algorithm 3 can be a useful tool for learning mixtures of
general RUMs.

We now consider an expressive generative model for mixtures of $K$ PLs where
Algorithm 3 produces a provably accurate estimate. The generative model
assumes that for all mixture components, only the utilities of the first $L$
items are different while the those of the remaining $n-L$ items are the same.
This model reflects the phenomenon where users from different sub-populations
differ in their preference among a few items while the remaining items are
essentially interchangeable. Intuitively, one would expect that when $L$ is
small, so is the difference between the subpopulations and it is harder to
separate the rankings into the correct clusters. On the other hand, when $L$
is large, the difference among the subpopulations is large and it is easier to
separate the clusters. The following theorem captures this intuition.

######  Theorem 3.1.

Consider a mixture of $K$ Plackett-Luce models with uniform mixing
probabilities. Suppose that $\theta^{*k}_{i}=0\,\forall i\in[L+1:n]$ and
$\theta^{*k}_{1:L}\sim\mathbb{N}(0,I_{L})$ for $k\in[K]$. Fix a constant
$\alpha>0$. There exist constants $c,c_{1},C_{1},C_{2},D$ such that if $m\geq
c\max\\{K^{4},Kn\\}$ then the output $\hat{\bm{\theta}}$ of Algorithm 3
satisfies the following. If $L\geq c_{1}\,\mathrm{exp}\left(C_{1}\sqrt{\log
n}\right)$, then

$\text{dist}(\hat{\bm{\theta}},\bm{\theta}^{*})=O\bigg{(}\mathrm{exp}\left(D\sqrt{\log
n}\right)\bigg{(}\sqrt{\frac{K^{2}n\log
n}{m}}+\frac{\sqrt{Kn}}{e^{L^{0.99}}}\bigg{)}\bigg{)}\,$

with probability
$1-O(\frac{K}{n^{8}})-O(K^{2}n^{2}\mathrm{exp}\left(-L^{0.99}\right))$. If
$L\geq C_{2}n^{\alpha}$ and assuming that $n=\omega(\log m)$, then

$\text{dist}(\hat{\bm{\theta}},\bm{\theta}^{*})=O\bigg{(}\mathrm{exp}\left(D\sqrt{\log
n}\right)\,\sqrt{\frac{K^{2}n\log n}{m}}\,\bigg{)}\,$

with probability
$1-O(\frac{K}{n^{8}})-O(K^{2}n^{2}\mathrm{exp}\left(-n^{\alpha}\right))$.

The first error bound is a sum of two terms. The first is the estimation error
incurred by Algorithm 2 which diminishes with increasing $m$. The second comes
from the clustering error incurred by Algorithm 1 and is controlled by the SNR
of the generative model. One can also check that $\mathrm{exp}\left(\sqrt{\log
n}\right)=o(n^{\alpha})$ for any $\alpha>0$ and $\mathrm{exp}\left(\sqrt{\log
n}\right)=\omega(\log n)$. When $L\approx\mathrm{exp}\left(O(\sqrt{\log
n})\right)$ (low SNR), there is significant clustering error and the second
term scales approximately as
$O(\frac{\sqrt{n}}{e^{L}})=O\big{(}\frac{1}{\text{poly}(n)}\big{)}$. Hence,
Algorithm 3 converges to within a small radius around $\bm{\theta}^{*}$ given
a sufficiently large $m$. However, when $L$ is polynomial in $n$ (high SNR),
estimation error dominates clustering error, giving us the second error bound
which diminishes with sample size $m$. In this regime, the spectral
initialization algorithm works well as a _standalone mixture learning
algorithm_. Note that this guarantee holds even for a small $\alpha>0$, when
the fraction of ‘informative’ items diminishes: $L/n=o(1)$. Our proposed
generative model is new and could be a useful analysis framework for future
works. To the best of our knowledge, the finite sample error bounds are also
the _first of their kind in the literature_.

####  3.2 Iterative Refinement via EM

Accuracy of the M-step. The following theorem generalizes Theorem 1 of Maystre
and Grossglauser (2015).

######  Theorem 3.2.

The output of weighted LSR (Algorithm 4) is the maximum weighted log-
likelihood estimate:

$\theta_{q}^{\text{MLE}}:=\arg\max_{\theta}\sum_{l=1}^{m}\,q_{l}\cdot\log\mathbb{P}^{\text{PL}}(\pi_{l},z_{l}\,\lvert\,\theta)\,.$

As noted before, the EMM algorithm is an alternative approach that exactly
solves the M-step using the (weighted) MM algorithm. In other words, _assuming
perfect numerical precision and the same initialization_ , EMM and EM-LSR will
produce the same final estimate. However, our EM-LSR algorithm is often much
faster than EMM (e.g., Figure 1).

Convergence of EM. It is well known that the EM algorithm converges to a
stationary point (Wu 1983). There is, unfortunately, no guarantee how close
such a point is to the global optimum. However, assuming correct model
specification and that the initial estimate falls within a neighborhood around
$\theta^{*}$ which satisfies certain high SNR conditions, the EM algorithm
will converge to $\theta^{*}$ (Wang et al. 2015; Wu et al. 2016; Balakrishnan,
Wainwright, and Yu 2017). The area around $\theta^{*}$ where this desirable
behaviour occurs is referred to as the _basin of attraction_. We leave the
detailed characterization of the basin of attraction as a subject of future
studies.

True Likelihood versus Surrogate Likelihood. For two other commonly used EM
algorithms in the literature – EM-CML and EM-GMM – previous authors use random
initialization. On the other hand, ours uses spectral initialization. However,
initialization is not the only differentiating characteristic of our
algorithm. In fact, our algorithm, EM-CML and EM-GMM are _fundamentally
different EM-based algorithms_. To see why, one needs to inspect the objective
function of the M-step. Suppose that all three algorithms are initialized at
some $\hat{\bm{\theta}}^{(0)}$. Let $\\{q^{k}_{l}\\}_{l\in[m]}^{k\in[K]}$
denote the posterior class probabilities conditioned on
$\hat{\bm{\theta}}^{(0)}$ per Equation (5).

In the first iteration, EM-LSR and EMM maximize the weighted log-likelihood.

$\hat{\bm{\theta}}_{\text{LSR}}^{(1)}=\arg\max_{\bm{\theta}}\sum_{l=1}^{m}\sum_{k=1}^{K}\bigg{[}q^{k}_{l}\cdot\log\mathbb{P}^{PL}(\pi_{l},z_{l}\,\lvert\,\theta^{k})\bigg{]}\,.$

On the other hand, EM-CML maximizes the _composite (surrogate) marginal
likelihood_. $\hat{\bm{\theta}}_{\text{CML}}^{(1)}=\arg\max_{\bm{\theta}}$

$\sum_{l=1}^{m}\sum_{k=1}^{K}\bigg{[}\sum_{\begin{subarray}{c}i,j:\\\
\pi_{l}(i)<\pi_{l}(j)\end{subarray}}q^{k}_{l}\log\bigg{(}\frac{1}{1+\mathrm{exp}\left(-(\theta^{k}_{i}-\theta^{k}_{j})\right)}\bigg{)}\bigg{]}\,.$

Lastly, EM-GMM _minimizes_ the following function.

$\hat{\bm{\theta}}_{\text{GMM}}^{(1)}=\arg\min_{\bm{\theta}}\sum_{k=1}^{K}\sum_{i\neq
j}\bigg{(}\hat{F}^{k}_{ij}-\frac{1}{1+\mathrm{exp}\left(-(\theta^{k}_{i}-\theta^{k}_{j})\right)}\bigg{)}^{2}\,,$

where
$\hat{F}^{k}_{ij}=\frac{\sum_{l=1}^{m}\mathbf{1}[\pi_{l}(i)<\pi_{l}(j)]\,q_{l}^{k}}{\sum_{l=1}^{m}q_{l}^{k}}$.

One can see that the objective functions are different and so are their
solutions. Hence, even if we initialize all three algorithms with the same
estimate, their trajectories will be different in general. While EM-LSR and
EMM converges to the true MLE when initialized within the basin of attraction,
this _may not be true for EM-GMM and EM-CML_. This difference is supported by
our experiments, where even with the same initialization, the algorithms
produce different final estimates.

###  4 Experiments

We compare our spectral EM algorithm to the following baselines: EMM, EM-GMM
and EM-CML.

(a) For a small number of mixture components, all methods are quite accurate.
As the theory implies, EMM and EM-LSR produce similar estimates.

(b) EM-LSR is comparatively efficient (figure shows total inference time).

(c) For a moderate number of mixture components, EM-LSR is the most accurate
method (EMM not shown due to timeout).

(d) EM-LSR and EM-CML are the only two methods that are efficient for a
moderate number of mixture components.

Figure 1: $\ell_{2}$ error and inference time on synthetic datasets. EM-LSR
(in blue) is competitive in terms of accuracy and speed to the baseline
algorithms.

Synthetic Datasets. We simulate data from the generative model as described in
Theorem 3.1. Specifically, we set $n=100$ and $L=5$ while varying the number
of mixture components $K$ for different experiments. Figure 1 shows estimation
error and total inference time against the sample size $m$, averaged over 25
trials. Experimentally, spectral initialization consistently gives better
initial estimates than both random initialization and GMM initialization
(Zhao, Piech, and Xia 2016). To keep a fair comparison, we use spectral
initialization for all algorithms. When $K$ is small (e.g., Figures 1(a) and
1(b)) all four methods are quite accurate. When the number of mixture
components are moderate (e.g., Figures 1(c) and 1(d)), the advantages that EM-
LSR enjoys over the other methods become more apparent. While EMM becomes too
inefficient for practical purposes, EM-LSR remains relatively efficient and
produces more accurate estimates than both EM-CML and EM-GMM.

Real Datasets. We include commonly used datasets in previous works such as
APA, Irish Elections (West, North, Meath) and SUSHI all with $n<15$. We
partition all the rankings with a 80-20 training-testing split; and the train
rankings into 80% for inference and 20% for validation. $K$ is chosen using
Bayesian Information Criterion (Gelman, Hwang, and Vehtari 2014) on the
validation set and the log-likelihood of the final model is evaluated using
the test set. For these datasets, EM-LSR and EMM are the most accurate while
EM-CML is the fastest, especially on datasets with a large $m$ such as the
Irish election datasets. We have a possible explanation for the relative speed
between EM-LSR and EM-CML. The bottle neck in these EM algorithms is the
M-step. The most time-consuming procedure in the M-step of EM-LSR is
constructing the Markov chain in Algorithm 4 with time complexity $O(mn^{2})$.
For EM-CML, it is solving a constrained concave maximization problem via SLSQP
(Virtanen et al. 2020) and may scale at least as $\Omega(n^{3})$111SLSQP
solves a sequence of quadratic optimization problems with $n$ variables. Each
solves a linear system with $n$ variables and $n$ equations and generally
takes $O(n^{3})$ (Strang et al. 1993).. Therefore, EM-CML tends to be faster
for datasets with a small $n$ and a large $m$. However, its inference time
could grow significantly with $n$.

Indeed, the setting where EM-LSR outperforms the baselines is when $n$ is
large. We perform additional experiments on the ML-10M movie ratings datasets
(Harper and Konstan 2015). To generate rankings, we first run a low rank
matrix completion algorithm (Zitnik and Zupan 2012) on the user-item rating
matrix to fill in the missing entries. We then select $n$ movies from the set
of all movies and the rankings are obtained from the completed matrix. Figure
2 shows the performance of the four methods on two versions of the ML-10M
datasets with $n=25$ and $n=100$ given increasing training data up to 14k. In
the supplementary materials, we also include additional experiments,
strategies to extend EM-LSR to handle partial rankings with ties and
comparisons to a Bayesian method (Mollica and Tardella 2017).

Dataset | Test log-likelihood | Inference time (s)
---|---|---
| EM-LSR | EM-CML | EM-GMM | EMM | EM-LSR | EM-CML | EM-GMM | EMM
APA ($n=5$) | -4.619 | -4.656 | -4.617 | -4.614 | 598 | 33.3 | 2.2K | 9.24K
West ($n=9$) | -11.9 | -12.008 | -11.904 | -11.896 | 810 | 199 | 5.75K | 25.8K
Sushi ($n=10$) | -13.64 | -14.0 | -13.773 | -13.766 | 746 | 24.6 | 489 | 1.22K
North ($n=12$) | -18.67 | -18.923 | -18.744 | -18.711 | 1.51K | 120 | 3.09K | 14K
Meath ($n=14$) | -23.645 | -23.885 | -23.69 | -23.633 | 1.48K | 497 | 29.9K | 69.1K
ML-10M ($n=25$) | -49.191 | -50.095 | -49.766 | -49.186 | 3.71K | 2.47K | 25.7K | 63K
ML-10M ($n=50$) | -130.499 | -132.209 | -132.143 | NA | 5.77K | 6.8K | 125K | NA
ML-10M ($n=100$) | -325.873 | -329.629 | -331.53 | NA | 11.5K | 27.2K | 492K | NA
ML-10M ($n=150$) | -550.462 | -557.923 | NA | NA | 14.7K | 62.1K | NA | NA
ML-10M ($n=200$) | -787.244 | -799.036 | NA | NA | 24.5K | 81.3K | NA | NA

Table 1: Test log-likelihood and inference time on real datasets. ‘NA’ denotes
not available due to timeout.

(a) With a small number of items, EMM and EM-LSR are more accurate.

(b) However, EM-LSR is also competitive in terms of efficiency.

(c) For a larger set of items, EM-LSR is the most accurate method (EMM not
shown due to timeout).

(d) EM-LSR is comparatively scalable for larger datasets.

Figure 2: Test log-likelihood and inference time on ML-10M datasets. For
larger datasets, EM-LSR (in blue) is more accurate while being competitive in
speed to the baseline algorithms.

###  5 Conclusion

We have proposed an accurate and efficient algorithm for learning a mixture of
Plackett-Luce models. For future works, we would like to consider other
initialization methods such as the method of moments or tensor decomposition.
Detailed characterization of the basin of attraction within which the EM
algorithm converges to the true parameter is also a challenging open question.
On a more practical note, incorporating the representation power of deep
neural networks into our algorithm will further increase its utility for large
scale recommendation systems applications.

###  6 Acknowledgements

The authors thank the anonymous reviewers for their thoughtful suggestions and
comments. The authors are supported by NSF Grant DMS-2112099. A.Z.
acknowledges financial support from the Alfred H. Williams Faculty Scholar
award. Any opinions expressed in this paper are those of the authors and do
not necessarily reflect the views of the National Science Foundation.

### References

  * Azari Soufiani et al. (2013) Azari Soufiani, H.; Chen, W.; Parkes, D. C.; and Xia, L. 2013.  Generalized method-of-moments for rank aggregation.  _Advances in Neural Information Processing Systems_ , 26.
  * Balakrishnan, Wainwright, and Yu (2017) Balakrishnan, S.; Wainwright, M. J.; and Yu, B. 2017.  Statistical guarantees for the EM algorithm: From population to sample-based analysis.  _The Annals of Statistics_ , 45(1): 77–120.
  * Caron and Doucet (2012) Caron, F.; and Doucet, A. 2012.  Efficient Bayesian inference for generalized Bradley–Terry models.  _Journal of Computational and Graphical Statistics_ , 21(1): 174–196.
  * Caron, Teh, and Murphy (2014) Caron, F.; Teh, Y. W.; and Murphy, T. B. 2014.  Bayesian nonparametric Plackett–Luce models for the analysis of preferences for college degree programmes.  _The Annals of Applied Statistics_ , 8(2): 1145–1181.
  * Dempster, Laird, and Rubin (1977) Dempster, A. P.; Laird, N. M.; and Rubin, D. B. 1977.  Maximum likelihood from incomplete data via the EM algorithm.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 39(1): 1–22.
  * Diaconis (1989) Diaconis, P. 1989.  A generalization of spectral analysis with application to ranked data.  _The Annals of Statistics_ , 949–979.
  * Gelman, Hwang, and Vehtari (2014) Gelman, A.; Hwang, J.; and Vehtari, A. 2014.  Understanding predictive information criteria for Bayesian models.  _Statistics and computing_ , 24(6): 997–1016.
  * Gormley and Murphy (2008) Gormley, I. C.; and Murphy, T. B. 2008.  Exploring voting blocs within the Irish electorate: A mixture modeling approach.  _Journal of the American Statistical Association_ , 103(483): 1014–1027.
  * Guiver and Snelson (2009) Guiver, J.; and Snelson, E. 2009.  Bayesian inference for Plackett-Luce ranking models.  In _proceedings of the 26th annual international conference on machine learning_ , 377–384.
  * Hajek, Oh, and Xu (2014) Hajek, B.; Oh, S.; and Xu, J. 2014.  Minimax-optimal inference from partial rankings.  _Advances in Neural Information Processing Systems_ , 27.
  * Harper and Konstan (2015) Harper, F. M.; and Konstan, J. A. 2015.  The movielens datasets: History and context.  _Acm transactions on interactive intelligent systems (tiis)_ , 5(4): 1–19.
  * Hunter (2004) Hunter, D. R. 2004.  MM algorithms for generalized Bradley-Terry models.  _The annals of statistics_ , 32(1): 384–406.
  * Jiao and Vert (2015) Jiao, Y.; and Vert, J.-P. 2015.  The Kendall and Mallows kernels for permutations.  In _International Conference on Machine Learning_ , 1935–1944. PMLR.
  * Kamishima (2003) Kamishima, T. 2003.  Nantonac collaborative filtering: recommendation based on order responses.  In _Proceedings of the ninth ACM SIGKDD international conference on Knowledge discovery and data mining_ , 583–588.
  * Keshavan, Montanari, and Oh (2009) Keshavan, R.; Montanari, A.; and Oh, S. 2009.  Matrix completion from noisy entries.  _Advances in neural information processing systems_ , 22.
  * Lei and Rinaldo (2015) Lei, J.; and Rinaldo, A. 2015.  Consistency of spectral clustering in stochastic block models.  _The Annals of Statistics_ , 43(1): 215–237.
  * Liu et al. (2019) Liu, A.; Zhao, Z.; Liao, C.; Lu, P.; and Xia, L. 2019.  Learning plackett-luce mixtures from partial preferences.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 33, 4328–4335.
  * Lomeli et al. (2018) Lomeli, M.; Rowland, M.; Gretton, A.; and Ghahramani, Z. 2018.  Antithetic and Monte Carlo kernel estimators for partial rankings.  _arXiv preprint arXiv:1807.00400_.
  * Luce (1959) Luce, R. D. 1959.  _Individual choice behavior: A theoretical analysis_.  Courier Corporation.
  * Mattei and Walsh (2013) Mattei, N.; and Walsh, T. 2013.  Preflib: A library for preferences http://www. preflib. org.  In _International Conference on Algorithmic Decision Theory_ , 259–270. Springer.
  * Maystre (2015) Maystre, L. 2015.  Choix.  https://github.com/lucasmaystre/choix.
  * Maystre and Grossglauser (2015) Maystre, L.; and Grossglauser, M. 2015.  Fast and accurate inference of Plackett–Luce models.  _Advances in neural information processing systems_ , 28.
  * Mollica and Tardella (2017) Mollica, C.; and Tardella, L. 2017.  Bayesian Plackett–Luce mixture models for partially ranked data.  _Psychometrika_ , 82(2): 442–458.
  * Peng and Wei (2007) Peng, J.; and Wei, Y. 2007.  Approximating k-means-type clustering via semidefinite programming.  _SIAM journal on optimization_ , 18(1): 186–205.
  * Plackett (1975) Plackett, R. L. 1975.  The analysis of permutations.  _Journal of the Royal Statistical Society: Series C (Applied Statistics)_ , 24(2): 193–202.
  * Rajkumar and Agarwal (2016) Rajkumar, A.; and Agarwal, S. 2016.  When can we rank well from comparisons of O (n$\backslash$log (n)) non-actively chosen pairs?  In _Conference on Learning Theory_ , 1376–1401. PMLR.
  * Shah and Song (2018) Shah, D.; and Song, D. 2018.  Learning RUMs: Reducing Mixture to Single Component via PCA.  _arXiv preprint arXiv:1812.11917_.
  * Strang et al. (1993) Strang, G.; Strang, G.; Strang, G.; and Strang, G. 1993.  _Introduction to linear algebra_ , volume 3.  Wellesley-Cambridge Press Wellesley, MA.
  * Thurstone (1927) Thurstone, L. L. 1927.  A law of comparative judgment.  In _Scaling_ , 81–92. Routledge.
  * Tkachenko and Lauw (2016) Tkachenko, M.; and Lauw, H. W. 2016.  Plackett-luce regression mixture model for heterogeneous rankings.  In _Proceedings of the 25th ACM International on Conference on Information and Knowledge Management_ , 237–246.
  * Vershynin (2018) Vershynin, R. 2018.  _High-dimensional probability: An introduction with applications in data science_ , volume 47.  Cambridge university press.
  * Virtanen et al. (2020) Virtanen, P.; Gommers, R.; Oliphant, T. E.; Haberland, M.; Reddy, T.; Cournapeau, D.; Burovski, E.; Peterson, P.; Weckesser, W.; Bright, J.; van der Walt, S. J.; Brett, M.; Wilson, J.; Millman, K. J.; Mayorov, N.; Nelson, A. R. J.; Jones, E.; Kern, R.; Larson, E.; Carey, C. J.; Polat, İ.; Feng, Y.; Moore, E. W.; VanderPlas, J.; Laxalde, D.; Perktold, J.; Cimrman, R.; Henriksen, I.; Quintero, E. A.; Harris, C. R.; Archibald, A. M.; Ribeiro, A. H.; Pedregosa, F.; van Mulbregt, P.; and SciPy 1.0 Contributors. 2020.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  _Nature Methods_ , 17: 261–272.
  * Wang et al. (2015) Wang, Z.; Gu, Q.; Ning, Y.; and Liu, H. 2015.  High dimensional em algorithm: Statistical optimization and asymptotic normality.  _Advances in neural information processing systems_ , 28.
  * Wu et al. (2016) Wu, C.; Yang, C.; Zhao, H.; and Zhu, J. 2016.  On the convergence of the em algorithm: A data-adaptive analysis.  _arXiv preprint arXiv:1611.00519_.
  * Wu (1983) Wu, C. J. 1983.  On the convergence properties of the EM algorithm.  _The Annals of statistics_ , 95–103.
  * Zhang and Zhou (2022) Zhang, A. Y.; and Zhou, H. H. 2022.  Leave-one-out Singular Subspace Perturbation Analysis for Spectral Clustering.  _arXiv preprint arXiv:2205.14855_.
  * Zhao, Liu, and Xia (2020) Zhao, Z.; Liu, A.; and Xia, L. 2020.  Learning Mixtures of Plackett-Luce Models with Features from Top-$l$ Orders.  _arXiv preprint arXiv:2006.03869_.
  * Zhao, Piech, and Xia (2016) Zhao, Z.; Piech, P.; and Xia, L. 2016.  Learning mixtures of Plackett-Luce models.  In _International Conference on Machine Learning_ , 2906–2914. PMLR.
  * Zhao, Villamil, and Xia (2018) Zhao, Z.; Villamil, T.; and Xia, L. 2018.  Learning mixtures of random utility models.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 32.
  * Zhao and Xia (2018) Zhao, Z.; and Xia, L. 2018.  Composite marginal likelihood methods for random utility models.  In _International Conference on Machine Learning_ , 5922–5931. PMLR.
  * Zhao and Xia (2019) Zhao, Z.; and Xia, L. 2019.  Learning mixtures of plackett-luce models from structured partial orders.  _Advances in Neural Information Processing Systems_ , 32.
  * Zitnik and Zupan (2012) Zitnik, M.; and Zupan, B. 2012.  Nimfa: A Python Library for Nonnegative Matrix Factorization.  _Journal of Machine Learning Research_ , 13: 849–853.