# Multi-objective Bayesian optimisation for design of Pareto-optimal current
drive profiles in STEP

Theodore Brown, Stephen Marsden, Vignesh Gopakumar, Alexander Terenin, Hong
Ge, and Francis Casson Manuscript received 4th October 2023. This work was
supported by the UKAEA/STFC FARSCAPE programme.Theodore Brown was with the
Department of Engineering, University of Cambridge, and is now with the
Department of Electronic and Electrical Engineering, University College
London. He is also with the United Kingdom Atomic Energy Authority (e-mail:
theo.brown@ukaea.uk).Stephen Marsden, Vignesh Gopakumar, and Francis Casson
are with the United Kingdom Atomic Energy Authority.Alexander Terenin was with
the Department of Engineering, University of Cambridge, and is now with the
Center for Data Science for Enterprise and Society, Cornell University.Hong Ge
is with the Department of Engineering, University of Cambridge.

###### Abstract

The safety factor profile is a key property in determining the stability of
tokamak plasmas. To improve the safety factor profile in the United Kingdom’s
proposed Spherical Tokamak for Energy Production (STEP), we apply _multi-
objective Bayesian optimisation_ —a machine learning technique in which an
uncertainty-aware predictive model guides the optimisation process based on
the observed data—to design electron-cyclotron heating profiles. The resulting
procedure generates sets of solutions that represent optimal tradeoffs between
objectives, enabling decision makers to understand the compromises made during
the design process. The solutions from our method score as well as those
generated in previous work by a genetic algorithm; however, our method leads
to a greater degree of solution diversity and interpretability, providing more
information to tokamak designers without compromising performance. Our results
suggest that the region of reversed central safety factor in STEP can be
reduced at the cost of letting rational safety factor surfaces move slightly
inwards, while keeping the minimum safety factor above two.

###### Index Terms:

Fusion reactor design, optimisation methods, Gaussian processes, machine
learning

##  I Introduction

The starting point in the development of high-performance tokamak scenarios is
the selection of the desired plasma properties during the main phase of
operation, known as the flat-top operating point (FTOP). The chosen flat-top
scenario determines not only the energy-generating performance, but also the
controllability and stability of the plasma [1].

The design space for tokamak devices is very large, involving the interaction
of many hundreds of parameters and many different physical regimes:
consequently, it is very difficult to find designs that meet the criteria and
constraints for a successful device. Exploring the space is a particular
challenge for the UK’s flagship Spherical Tokamak for Energy Production
(STEP), as the novelty of spherical tokamak technology and power-generating
plasmas introduces significant uncertainty into the design process. In such a
vast and changing space, it is important to be able to quantify the tradeoffs
between competing objectives and quickly observe the effects of reformulating
or adjusting the FTOP specification.

FTOP candidates for STEP are simulated using JETTO, a plasma modelling code
that iteratively solves coupled magnetic equilibrium and plasma transport
equations [2]. Currently, JETTO takes around 3 hours to simulate a new FTOP on
state-of-the-art high performance computing, which involves simulating for
several resistive timescales and converging to a steady state with simple
feedback controllers.111Benchmarks performed on a single 2.20GHz Intel Xeon
Platinum 8276 CPU with 3380 MiB of RAM, hosted by the Cambridge Service for
Data-driven Discovery. As the main cost is evolving the temporal rather than
the spatial grid, this particular problem has minimal benefit from further
parallelisation. As scenario optimisation requires many hundreds of JETTO
runs, the computational cost involved in iterating a design is significant.

Therefore, we ask: _how can we better resolve the tradeoffs between different
design objectives?_ We demonstrate that multi-objective analysis provides a
robust and interpretable approach to design optimisation, enabling the STEP
team to better understand and formulate the FTOP specification.

Secondly, we ask: _how can we accelerate the design iteration process?_ We
demonstrate that using a probabilistic method allows a diverse set of optimal
solutions to be generated without any loss of performance, reducing the
computational cost involved in iterating a design compared to previous
methods.

Many plasma properties of interest vary from the core to the edge of the
plasma as 1-D functions of radius, and are known as _profiles_. We demonstrate
how multi-objective optimisation can be applied to design electron-cyclotron
resonance heating (ECRH) deposition profiles, with the goal of shaping safety
factor ($q$) profiles. We do so via _Bayesian optimisation_ , a generic
technique that can be applied to any tokamak design problem.

##  II Design optimisation in fusion research

Multi-objective Bayesian optimisation has previously been applied to the
design of toroidal field coils for future fusion reactors, demonstrating
improved performance compared to multi-objective genetic algorithms [3].

The optimisation of STEP current drive profiles using a scalar genetic
algorithm was presented in [4]. Our paper extends this work by applying multi-
objective Bayesian optimisation to the same problem, and providing a
comparison in terms of performance and quality of solutions.

For a discussion of the plasma scenarios used in our simulations, see [5, 6].
A comprehensive review of the STEP scenarios will be presented in a
forthcoming paper.

##  III Multi-objective Bayesian optimisation of safety factor profiles for
STEP

We now present the safety factor optimisation problem, and formulate it as a
multi-objective Bayesian optimisation task.

###  III-A Safety factor properties and ECRH parameterisation

TABLE I: Objective functions for the safety factor profile ($q$-profile) ID |  Desired property | Untransformed metric | Target
---|---|---|---
1 |  Minimise height of on-axis reversed shear | $\|q_{0}-q_{\mathrm{min}}\|$ | $0$
2 |  Minimise width of on-axis reversed shear | $\hat{\rho}(q=q_{\mathrm{min}})$ | $0$
3 |  $q$ increasing | $\frac{1}{N}\sum_{i=1}^{N}\mathbf{1}\Big{(}\Big{[}\frac{dq}{d\hat{\rho}}\Big{]}_{i}>0\Big{)}$ | $1$
4 |  $q_{\mathrm{min}}$ within safe range | $q_{\mathrm{min}}$ | $[2.2,2.5]$
5 |  Maximise radius of $q=3$ surface | $\hat{\rho}(q=3)$ | $1$
6 |  Maximise radius of $q=4$ surface | $\hat{\rho}(q=4)$ | $1$

(a)

(b)

Figure 1: Example random ECRH profiles (i.e., before any optimisation).
‘Piecewise linear’ consists of 8 linear segments. ‘Bézier’ is a single order 7
Bézier curve. Segment ends and Bézier control points are shown with markers.

The shape of the $q$-profile plays an important role in improving the plasma
confinement and minimising the impact of magnetohydrodynamic instabilities.

The planned scenarios for STEP exhibit a monotonically increasing $q$-profile
with $q(0)>2$ [5]. Monotonic $q$ profiles have been identified as contributing
to increased stable normalised beta, $\beta_{N}$, [7] the avoidance of
internal transport barriers [8] and reduced fast particle instabilities [9].
Maintaining a monotonic $q$ also means that the plamsa is operating in a
regime where magnetohydrodynamic instabilities are generally better understood
[10]. Keeping $q(0)>2$ avoids the formation of the $(2,1)$ neoclassical
tearing mode (NTM), allowing stable operation at high $\beta_{N}$ [11].
Additional NTMs occur at rational values of $q$ [1], and have been shown to
exhibit improved stability if the rational $q$ surfaces are nearer the edge of
the plasma [12] (although NTMs can also be destabilised by high magnetic
shear, which is likely to occur at large minor radii).

The exact mathematical representations of the objective functions have a
significant impact on the solutions generated. To quantify monotonicity, we
compute the gradient of the $q$-profile at $N$ radial locations, and then find
the fraction of locations for which the gradient is positive. In addition,
STEP FTOP candidates often have a region of reversed $q$ at the plasma centre,
due to the vanishing of the bootstrap current on axis: we penalise this
undesirable behaviour with two objectives that quantify the radial width and
magnitude of reversed shear. We set the target range for $q_{\text{min}}$ to
be $2.2-2.5$. We transform each of these desired $q$ properties using a custom
function, which we call SoftHat, that rescales them to be maximisation
objectives in $[0,1]$. The function is detailed in Appendix A. The objectives
are shown in Table I, prior to transformation.

Note that our choice of objectives is not intended to be definitive: instead,
the understanding obtained from optimisation against these heuristics will be
used to narrow future high-resolution evaluations of the scenario space.

At present, we focus on purely theoretical design of ECRH power density
profiles (hereafter referred to as ‘QECE’). We do not seek to design a
launcher configuration or incoporate the physics of current drive. Instead, we
assume that any arbitrary profile can be achieved by the ECRH launchers, and
define general function parameterisations that can represent a wide variety of
possible QECE profiles.

The first parameterisation we consider is a piecewise linear function with 12
parameters, introduced in [4]. This function was designed to produce a limited
range of plausible QECE profiles, with an on-axis peak and optional off-axis
peak.

The second parameterisation is a Bézier curve, which is widely used in
engineering design to represent complicated smooth shapes [13]. Bézier curves
are a useful choice as they are bounded by the convex hull defined by the
control points; hence, the curve can be constrained to be positive, which is
difficult to achieve under other smooth parameterisations (such as cubic
splines). The Bézier curve we use has 8 parameters.

The resulting power density profile is scaled such that the total ECRH power
is a predefined constant value. Example profiles from each parameterisation
are shown in Figure 1. See Appendix B for detailed mathematical definitions.

###  III-B Multi-objective optimisation

The ECRH-$q$ optimisation problem is an example of a multi-objective
optimisation (MOO) problem. A common approach to MOO problems is to cast the
multi-objective (vector) objective function to a single-objective (scalar) one
by computing a weighted sum of the components. However, using a weighting
function introduces preconceptions about the desired balance between the
objectives, and the optimisation routine will only produce solutions that
align with the pre-determined tradeoff. Moreover, tuning the weights is often
a resource-intensive process, and commonly relies on human intuition.

Instead, taking a ‘true’ multi-objective (vector) approach allows us to gain
understanding of how the objectives interact and conflict, and facilitates
reasoned decision-making after observing a wide variety of possible solutions
to the tradeoff. To perform true multi-objective analysis, MOO algorithms seek
to find a representative set of Pareto optimal solutions. More precisely, a
solution $s$ is _Pareto-optimal_ if there is no other solution that is (i) at
least as good as $s$ in all objectives and (ii) better than $s$ in at least
one objective. The set of all Pareto-optimal solutions is known as the Pareto
front.

###  III-C Bayesian optimisation

Bayesian optimisation (BO) is a principled approach to performing gradient-
free global optimisation of costly black-box functions. For a tutorial on BO,
see [14]; we provide a brief summary of the method below and in Figure 2.

Global optimisation tasks involve a compromise between exploitation and
exploration: the optimiser must decide whether to evaluate trial solutions
that are close to its current best guess of the optimal solution
(‘exploitation’) or try solutions from an as-yet unobserved region
(‘exploration’). Many global optimisation algorithms tackle this dilemma using
stochastic exploration; for example, in a genetic algorithm, random mutations
are introduced in every generation. Bayesian optimisation instead works by
building a probabilistic model of the mapping from inputs $s$ to objectives
$\phi(s)$, and using the model’s uncertainty predictions to resolve the
explore-exploit tradeoff.

Firstly, a probabilistic model is fit to all of the previously observed input-
objective pairs, $\mathcal{D}=\left\\{(s_{i},\phi(s_{i})\right\\}_{i=1}^{n}$,
by maximising the marginal likelihood.222An initial set of observations is
generated by selecting inputs using random or pseudo-random sampling, and
observing the corresponding outputs. The model used is a Gaussian process
(GP). A GP is an extension of the Gaussian distribution to infinite
dimensions, and can be understood as representing a probability distribution
over functions [15].

###### Definition 1

Gaussian process. A Gaussian process over functions
$f:\mathcal{X}\to\mathbb{R}$ is defined by a mean function
${m:\mathcal{X}\to\mathbb{R}}$ and a positive semi-definite covariance
function ${k:\mathcal{X}\times\mathcal{X}\to\mathbb{R}}$ such that
$k(\cdot,\cdot)\geq 0$. A GP is denoted by

$f\sim\mathcal{GP}\left(m,k\right).$ (1)

Then, the values of $f$ at a finite set of evaluation points
$\mathbf{x}=\left[x_{1},\dots,x_{N}\right]\subset\mathcal{X}$ are distributed
according to

$f(\mathbf{x})\sim\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma}),$ (2)

where the vector $\boldsymbol{\mu}=\left[m(x_{1}),m(x_{2}),\dots\right]^{T}$
and the matrix $\boldsymbol{\Sigma}$ has elements
$\left[\boldsymbol{\Sigma}\right]_{ij}=k(x_{i},x_{j})\ \forall
i,j\in\\{1,\dots,N\\}$.

In our example problem, for any candidate ECRH profile $s$, the fitted GP
provides a predictive probability distribution for each objective value
$\phi(s)$ given in Table I.

To select the next points to trial, BO uses the predictive model in
conjunction with an _acquisition function_ ,
$\alpha:\mathcal{S}\to\mathbb{R}$, which estimates the ‘usefulness’ of
observing the true objective value of a previously unseen candidate
$s\in\mathcal{S}$. Gradient ascent is performed on $\alpha$, using only the GP
predictions, to find the candidate solution that maximises the acquisition
function. The objective value of the proposed candidate is then evaluated (in
our case, by running JETTO to generate a $q$-profile, and computing the
objective values from the $q$-profile properties), and the model is updated to
account for the new observation.

The choice of acquisition function determines the approach to resolving the
explore-exploit tradeoff. In scalar optimisation tasks, a common acquisition
function is expected improvement, which provides useful intuition for the
multi-objective setting.

###### Definition 2

Expected improvement. Given a set of observations
$\mathcal{D}=\left\\{(s_{i},\boldsymbol{\phi}(s_{i})\right\\}_{i=1}^{n}$ with
the highest observed objective denoted by
$\boldsymbol{\phi}^{*}=\max_{i}\boldsymbol{\phi}(s_{i})$, and a probabilistic
model with conditional density
$p(\cdot|\cdot,\mathcal{D}):\mathbb{R}\times\mathcal{S}\to\mathbb{R}$, the
expected improvement is given by

$\alpha_{\mathrm{EI}}(s)=\int\max(\phi-\phi^{*},0)p(\phi|s,\mathcal{D})d\phi.$
(3)

The multi-objective analogue of this is Expected Hypervolume Improvement
(EHVI), which replaces the one-dimensional improvement with the increase in
multi-dimensional volume. It can be readily extended to support _batch
optimisation_ , where a set of candidates are generated concurrently, and
optimisation where the objectives are corrupted with noise. The resulting
function, qNEHVI, is a current state-of-the-art acquisition functions for
multi-objective Bayesian optimisation [16].

Figure 2: Schematic of Bayesian optimisation for the ECRH-$q$ optimisation
problem.

###  III-D Experimental setup

For our experiments, we fixed a budget of 8 evaluation steps, each with 32
parallel JETTO evaluations. This corresponds to approximately 24 hours of
wall-clock time, and is the same budget used for the scalar genetic algorithm
method in [4].

We use BoTorch [17] for Bayesian optimisation. We use a GP model with a
Matérn-$\frac{5}{2}$ kernel and independent outputs, using the default BoTorch
hyperparameters. These are standard choices in Bayesian optimisation, which we
selected for their simplicity. Note that a Matérn-$\frac{5}{2}$ kernel models
less-smooth functions than, say, a squared exponential kernel, leading to a
model that is slightly more robust. We generate the batches of $32$ candidate
solutions by performing joint optimisation of the qNEHVI acquisition function;
among the acquisition optimisation approaches we tried, this performed best.

We generate initial candidates using Sobol sampling. Our budget of $8\times
32$ evaluations was split between Sobol initialisation and Bayesian
optimisation steps. We found that the choice of ECRH parameterisation affected
the ideal number of Sobol candidates to evaluate: for the piecewise linear
ECRH profile we used $96$ Sobol samples, and for the Bézier profile we used
$128$. This discrepancy may reflect the greater diversity of profiles that can
be generated using the Bézier parameterisation. In our case, external
constraints limited the number of parallel observations to 32: in practice,
the $N$ Sobol samples can be observed completely in parallel, leading to
$1\times N$ Sobol evaluations and $7\times 32$ BayesOpt evaluations.

Our code is available on GitHub,333https://github.com/theo-brown/jetto-mobo
including a flexible Python package for using Bayesian optimisation with
JETTO.

(a) Piecewise linear

(b) Bézier

Figure 3: Solutions to the ECRH-$q$ optimisation problem generated by multi-
objective BO, compared to the solution from a scalar genetic algorithm (GA)
with piecewise linear ECRH (from [4]). The GA’s solution is still on the
Pareto front; however, multi-objective BO identifies many other optimal
solutions, as well as potential gains from different weighting of the
objectives.

##  IV Results

Our results are shown in Figure 3, alongside the original genetic algorithm
(GA) solution from [4]. We highlight a few key conclusions:

  1. (1)

Under both parameterisations, the GA solution is in the Pareto optimal set:
that is, no solution is as good as the GA in all objectives and better in at
least one objective.

  2. (2)

Significant improvements could be made in the central reversed shear and
monotonicity criteria, if the $q=3$ surface is allowed to move closer to the
centre. For example, Solution 1 in Figure 3 and Solution 4 in Figure 3
strongly outperform the previous solution in objectives 1-3, in exchange for a
modest loss of performance in objective 5.

  3. (3)

To completely flatten the reversed shear at the core, a very sharp on-axis
ECRH peak is required—see Figure 3, Solution 3. Such a steep peak may not be
physically realisable. Our Bézier profile was constrained to have zero
gradient at the origin, which meant it was unlikely to generate profiles with
sufficiently sharp peaks.

##  V Discussion and further work

We now discuss some of the challenges encountered, areas for improvement, and
overall takeaways.

_ECRH parameterisations._ One of the challenges in designing an ECRH
parameterisation is to ensure sufficient flexibility without making the search
space too large. For example, we found that an unconstrained piecewise linear
function was overly expressive, which meant that it was difficult to find
‘good’ regions of parameter space. It would be illuminating to perform studies
with ECRH parameterisations that better reflect realistic heating
configurations, such as Gaussian bumps representing the deposition from EC
beam launchers; this would allow optimisation of machine-relevant scenarios.
From an algorithmic perspective, it would also be interesting to explore
BayesOpt directly on an infinite-dimensional space of functions, rather than
on a finite-dimensional space such as the one defined by Bézier curves.

_Objective formulation._ The solutions generated by multi-objective Bayesian
optimisation are often tightly grouped around a few objective values, which
makes it difficult to differentiate between solutions. Further work in shaping
the objective functions, whether by adjusting the metrics, targets, or
transformations, could provide improvements in interpretability and encourage
greater diversity in the Pareto set.

_Objectives and constraints._ In general, multi-objective optimisation
algorithms tend to work best with a relatively small number of objectives, and
can exhibit significant decay in performance in so-called _many_ -objective
problems. Strategies for tackling many-objective problems include performing
dimensionality reduction to find a smaller set of uncorrelated objectives,
guiding the search algorithm by introducing objectives sequentially, and
scalarisation methods [18, 19]. In our example case, no decay in performance
is immediately evident, as the solutions are at least as good as those
generated previously. However, it is clear that some of our objectives are
correlated. For example, as $q$ is increasing everywhere to the right of
$q_{min}$, objectives 2 and 3 ($q_{\mathrm{min}}$ close to centre and $q$
increasing) are measuring almost identical quantities. Similarly, the
locations of the $q=3$ and $q=4$ surfaces are closely related (objectives 5
and 6). Moreover, many of the ‘objectives’ that we selected could be
reformulated as constraints: for example, the $2<q_{\mathrm{min}}<3$ objective
should perhaps be treated as a constraint, as any profiles with $q<2$ will not
be considered for STEP. Extensions such as these are a promising avenue for
future work, and could enable optimisation of more complex plasma properties.

_Acquisition functions._ Although qNEHVI has been shown to perform well in
multi-objective tasks, it has not been empirically demonstrated to be the best
choice in many-objective problems. Moreover, qNEHVI slows down at large batch
sizes, as the joint optimisation over samples becomes infeasible. Further work
is needed to explore whether alternative options such as variants of Thompson
sampling [20] deliver improved performance at scale for tokamak design
problems.

##  VI Conclusion

Using multi-objective Bayesian optimisation, we find a diverse range of ECRH
heating profiles that represent optimal tradeoffs between the desired
properties of the $q$-profile for STEP, without any increase in compute
requirements. With the same number of evaluations, our method generates
solutions that are as good as those generated with a genetic algorithm
operating on a scalarised objective. However, our solutions are widely spaced
along the Pareto front and hence represent diverse tradeoffs: in contrast, the
genetic algorithm’s solutions are clustered around a specific tradeoff,
determined by the scalar weighting of the objectives (see [4]). Although it is
possible to generate multi-objective sets of solutions using a genetic
algorithm, previous work on similar problems has demonstrated that the results
are lower-quality and dominated by those generated using Bayesian optimisation
[3].

Our set of solutions also highlights that slackening the requirement for the
$q=3$ surface to be towards the edge allows for dramatic improvements in
minimising the reversed shear at the centre of the plasma. We hope our work
encourages fusion scientists and engineers to adopt multi-objective Bayesian
optimisation to facilitate informed decision making when designing future
devices.

## References

  * [1] M. L. Walker, P. De Vries, F. Felici, and E. Schuster, “Introduction to tokamak plasma control,” in _American Control Conference_ , 2020, pp. 2901–2918.
  * [2] G. Cenacchi and A. Taroni, “JETTO: a free-boundary plasma transport code,” European Nuclear Energy Authority, Tech. Rep., 1988.
  * [3] T. Nunn, V. Gopakumar, and S. Kahn, “Shaping of magnetic field coils in fusion reactors using Bayesian optimisation,” in _NeurIPS Workshop on Gaussian Processes, Spatiotemporal Modeling, and Decision-making Systems_ , 2022\.
  * [4] S. Marsden, F. Casson, B. Patel, E. Tholerus, and T. Wilson, “Using genetic algorithms to optimise current drive in STEP,” in _Europhysics Conference Abstracts_ , vol. 46A, 2022, p. P2b.123.
  * [5] H. Meyer, “The physics of the preferred plasma scenario for STEP,” in _Europhysics Conference Abstracts_ , vol. 46A, 2022, p. O2.107.
  * [6] D. Kennedy, M. Giacomin, F. J. Casson, D. Dickinson, W. A. Hornsby, B. S. Patel, and C. M. Roach, “Electromagnetic gyrokinetic instabilities in the Spherical Tokamak for Energy Production (STEP) part I: linear physics and sensitivity,” 2023. [Online]. Available: arxiv.org/abs/2307.01670
  * [7] J. E. Menard, M. Bell, R. Bell, D. Gates, S. Kaye, B. LeBlanc, S. Sabbagh, E. Fredrickson, S. Jardin, R. Maingi _et al._ , “Unified ideal stability limits for advanced tokamak and spherical torus plasmas,” Princeton Plasma Physics Laboratory, United States, Tech. Rep., 2003.
  * [8] J. Connor, T. Fukuda, X. Garbet, C. Gormezano, V. Mukhovatov, M. Wakatani _et al._ , “A review of internal transport barrier physics for steady-state operation of tokamaks,” _Nuclear Fusion_ , vol. 44, no. 4, p. R1, 2004.
  * [9] R. Buttery, B. Covele, J. Ferron, A. Garofalo, C. Holcomb, T. Leonard, J. Park, T. Petrie, C. Petty, G. Staebler _et al._ , “DIII-D research to prepare for steady state advanced tokamak power plants,” _Journal of Fusion Energy_ , vol. 38, no. 1, pp. 72–111, 2019.
  * [10] S. Mahajan and R. Hazeltine, “Linear theory of electromagnetic instabilities in plasmas with a hollow q-profile,” _Nuclear Fusion_ , vol. 22, no. 9, p. 1191, 1982.
  * [11] J. Menard, M. Bell, R. Bell, D. Gates, S. Kaye, B. LeBlanc, R. Maingi, S. Sabbagh, V. Soukhanovskii, D. Stutman _et al._ , “Aspect ratio scaling of ideal no-wall stability limits in high bootstrap fraction tokamak plasmas,” _Physics of Plasmas_ , vol. 11, no. 2, pp. 639–646, 2004.
  * [12] P. Maget, H. Lütjens, R. Coelho, B. Alper, M. Brix, P. Buratti, R. Buttery, E. De la Luna, N. Hawkes, G. Huysmans _et al._ , “Modelling of (2,1) NTM threshold in JET advanced scenarios,” _Nuclear Fusion_ , vol. 50, no. 4, p. 045004, 2010.
  * [13] H. N. Fitter, A. B. Pandey, D. D. Patel, and J. M. Mistry, “A review on approaches for handling Bézier curves in CAD for manufacturing,” _Procedia Engineering_ , vol. 97, pp. 1155–1166, 2014.
  * [14] R. Garnett, _Bayesian Optimization_. Cambridge University Press, 2023.
  * [15] C. K. Williams and C. E. Rasmussen, _Gaussian Processes for Machine Learning_. MIT Press, 2006.
  * [16] S. Daulton, M. Balandat, and E. Bakshy, “Parallel Bayesian optimization of multiple noisy objectives with Expected Hypervolume Improvement,” in _Advances in Neural Information Processing Systems_ , 2021.
  * [17] M. Balandat, B. Karrer, D. Jiang, S. Daulton, B. Letham, A. G. Wilson, and E. Bakshy, “BoTorch: a framework for efficient Monte-Carlo Bayesian optimization,” in _Advances in Neural Information Processing Systems_ , 2020\.
  * [18] P. J. Fleming, R. C. Purshouse, and R. J. Lygoe, “Many-objective optimization: an engineering design perspective,” in _International Conference on Evolutionary Multi-criterion Optimization_ , 2005, pp. 14–32.
  * [19] R. Kasimbeyli, Z. K. Ozturk, N. Kasimbeyli, G. D. Yalcin, and B. I. Erdem, “Comparison of some scalarization methods in multiobjective optimization,” _Bulletin of the Malaysian Mathematical Sciences Society_ , vol. 42, p. 1875–1905, 2019.
  * [20] K. Kandasamy, A. Krishnamurthy, J. Schneider, and B. Póczos, “Parallelised Bayesian optimisation via Thompson sampling,” in _Artificial Intelligence and Statistics_ , 2018.

(a)

(b)

(c)

Figure 4: Demonstration of the SoftHat transformation function. The function
first increases as a squared exponential, passing through
$(x_{\textsc{lower}},y_{\textsc{lower}})$, then is uniformly $1$ between start
and end, before decreasing as a squared exponential, passing through
$(x_{\textsc{upper}},y_{\textsc{upper}})$. Setting $y_{\textsc{lower}}$ or
$y_{\textsc{upper}}$ to $0$ recovers a square pulse (‘top-hat’) function.

##  Appendix A SoftHat transformation

We defined a transformation to scale the objective values to $[0,1]$, where
$1$ is the ‘best’ value. Our key requirements were:

  * •

Easy to interpret

  * •

Smooth

  * •

Can represent maximisation and minimisation objectives

  * •

Can represent objectives where the desired value is within a given range

Our transformation is:

$\displaystyle\textsc{SoftHat}\left(x,x_{\textsc{lower}},y_{\textsc{lower}},\textsc{start},\textsc{end},x_{\textsc{upper}},y_{\textsc{upper}}\right)$
$\displaystyle=\begin{cases}\exp\left(-k_{\textsc{lower}}\left(x-\textsc{start}\right)^{2}\right)&x<\textsc{start}\\\
1&\textsc{start}\leq x\leq\textsc{end}\\\
\exp\left(-k_{\textsc{upper}}\left(x-\textsc{end}\right)^{2}\right)&x>\textsc{end}\mathrm{,}\\\
\end{cases}$

where

$k_{\textsc{lower}}=-\frac{\log(y_{\textsc{lower}})}{\left(x_{\textsc{lower}}-\textsc{start}\right)^{2}}$
$k_{\textsc{upper}}=-\frac{\log(y_{\textsc{upper}})}{\left(x_{\textsc{upper}}-\textsc{end}\right)^{2}}\mathrm{.}$

This choice of the coefficients sets the rate of decay such that the function
passes through $(x_{\textsc{lower}},y_{\textsc{lower}})$ and
$(x_{\textsc{upper}},y_{\textsc{upper}})$.

TABLE II: Parameter bounds for the piecewise linear parameterisation. Parameter | Label | Lower | Upper
---|---|---|---
$p_{0}$ | Normalised on-axis power | $0$ | $1$
$p_{1}$ | Normalised radius of $b$ | $0$ | $0.1$
$p_{2}$ | Normalised power at $b$ as a fraction of on-axis power | $0$ | $1$
$p_{3}$ | Normalised power fraction halfway between $b$ and $d$ | $0$ | $1$
$p_{4}$ | Normalised radius of $d$ | $0.1$ | $0.3$
$p_{5}$ | Normalised power at $d$ | $0$ | $1$
$p_{6}$ | Normalised power fraction $\frac{1}{3}$ between $d$ and $g$ | $0$ | $1$
$p_{7}$ | Normalised power fraction $\frac{2}{3}$ between $d$ and $g$ | $0$ | $1$
$p_{8}$ | Normalised radial distance between $d$ and $g$ | $0.1$ | $0.8$
$p_{9}$ | Normalised power at $g$ | $0$ | $1$
$p_{10}$ | Normalised power halfway between $g$ and turnoff ($i$) | $0$ | $1$
$p_{11}$ | Normalised radial distance between $g$ and turnoff | $0$ | $0.5$

TABLE III: Relationship between parameters and node
locations for the piecewise
linear parameterisation. Node | $x$ | $y$
---|---|---
$a$ | $0$ | $p_{0}$
$b$ | $p_{1}$ | $p_{2}a_{y}$
$c$ | $\frac{b_{x}+d_{x}}{2}$ | $p_{3}b_{y}+(1-p_{3})d_{y}$
$d$ | $p_{4}$ | $p_{5}$
$e$ | $\frac{2}{3}d_{x}+\frac{1}{3}g_{x}$ | $p_{7}f_{y}+(1-p_{7})d_{y}$
$f$ | $\frac{1}{3}d_{x}+\frac{2}{3}g_{x}$ | $p_{6}g_{y}+(1-p_{6})d_{y}$
$g$ | $d_{x}+p_{8}$ | $p_{9}$
$h$ | $\frac{g_{x}+i_{x}}{2}$ | $p_{10}g_{y}$
$i$ | $g_{x}+p_{11}$ | $0$

##  Appendix B ECRH power density profile parameterisations

###  B-A Piecewise linear

The piecewise linear function parameterisation was designed to permit a
variety of shapes; in particular, it can create a) a sharp peak on axis to
prevent the formation of a current hole, and either b) a trough next to the
on-axis peak and create an off-axis peak, or c) monotonically decreasing
profiles with no trough. The function is parameterised by 12 parameters (Table
III), defining 9 nodes (Table III). In Table III, components $x$ and $y$ are
the node’s position in normalised radius and normalised ECRH power density
respectively. A straight line is interpolated between each of the adjacent
pairs of nodes to form the normalised ECRH power density profile.

This formulation ensures that the nodes remain in the intended order. In
addition, defining the nodes in relation to one another means that adjusting
the radial position of the off-axis peak does not affect its overall shape.

###  B-B Bézier

The 8-dimensional parameter vector $\left(p_{0},\dots,p_{7}\right)$ is
combined with two smoothing hyperparameters $\delta_{0}$ and $\delta_{1}$,
which adjust the gradient of the curve at the end points, to produce the
control point vector
$\left((0,p_{0}),(\delta_{0},p_{0}),(p_{1},p_{2}),\dots,(p_{5},p_{6}),(p_{7}-\delta_{1},0),(p_{7},0)\right)$.
Note that this constrains the profile to finish at 0.


# AI Fairness via Domain Adaption
Application to AMD Detection

Neil Joshi  Johns Hopkins University Applied Physics Laboratory11100 Johns
Hopkins RoadLaurelMDUSA<EMAIL_ADDRESS>and  Phil Burlina  Johns Hopkins
University Applied Physics Laboratory11100 Johns Hopkins RoadLaurelMDUSA
<EMAIL_ADDRESS>
(2020)

###### Abstract.

While deep learning (DL) approaches are reaching human-level performance for
many tasks, including for diagnostics AI, the focus is now on challenges
possibly affecting DL deployment, including AI privacy, domain generalization,
and fairness. This last challenge is addressed in this study. Here we look at
a novel method for ensuring AI fairness with respect to protected or sensitive
factors. This method uses domain adaptation via training set enhancement to
tackle bias-causing training data imbalance. More specifically, it uses
generative models that allow the generation of more synthetic training samples
for underrepresented populations. This paper applies this method to the use
case of detection of age related macular degeneration (AMD). Our experiments
show that starting with an originally biased AMD diagnostics model the method
has the ability to improve fairness.

††copyright: none††journalyear: 2020††journal: PRE

##  1\. Introduction

Deep learning (DL) has now widely and successfully been applied to tasks
ranging from recognizing faces, detecting people, prescreening diseases or
performing lesion segmentation (pekala2019deep, ; burlina2019assessment, ),
and has outperformed past machine learning approaches (burlina2011automatic,
). However, new challenges have emerged and are potentially impacting future
deployment of AI, including: privacy issues (shokri2017membership, );
challenges that entail the existence of adversarial machine learning or data
poisoning to defeat DL models (carlini2017adversarial, ); the lack of
availability of sufficient training data (burlina2020low, ); or AI fairness
which may lead to lack of parity in prediction performance
(burlina2020addressing, ), the last concern is our focus here. With regard to
AI fairness one important cause of bias is data imbalance which is
specifically addressed herein.

Addressing AI fairness has recently gotten significant attention
(mehrabi2019survey, ; hardt2016, ; bolukbasi2016man, ; kinyanjui2020fairness,
; quadrianto2019discovering, ). Work in addressing AI bias falls roughly into
the following basic taxonomy and the following types of approaches: some
methods alter the model, some alter the data, and some perform a form of post-
processing via re-calibration.

Regarding approaches that alter the model one path that is now intensely
investigated looks at adversarial two player methods that tackle the possible
conditional dependence of predictions made by DL models (e.g. predicting AMD)
on protected attributes, and use an adversarial technique
(goodfellow2014generative, ) to promote independence and fairness (see
(ganin2016domain, ; beutel2017data, ) or studies like (wadsworth2018achieving,
) and (zhang2018mitigating, )). (beutel2017data, ) in particular uses an
adversarial player that predicts protected factors from the feature embeddings
(as does (wadsworth2018achieving, )) with the adversarial network trying to
minimize a loss function related to making an accurate attribute prediction,
while the main network is both trying to estimate the labels for the task at
hand (for example in our case AMD vs. no AMD) while disallowing the second
adversarial network from guessing correctly the protected factors
(zhang2018mitigating, ). A related approach is used in (song2019learning, )
via employing a loss function that uses information theoretic measures.
(quadrianto2019discovering, ) uses a method that employs a related adversarial
approach but makes changes to the image (instead of making changes to the
embedding representation) and applies this to masking of visual features in
faces.

One possible drawback of this first type of approaches is that a) adversarial
two player systems can yield challenges for training the networks and b) these
methods do not address inherent issues of bias emanating from data imbalance.
This motivates the type of approach used here which instead uses training data
augmentation and domain adaptation.

To address this important cause of bias, our approach proceeds via
synthesizing more data for underrepresented populations and does this via a
generative method which has the unique ability to perform fine control of
protected attributes, performing a sort of domain adaptation. This strategy
has yet to be explored to a large extent compared to the adversarial approach
to bias mitigation, and is a promising area of research that is focused on
specific aspects of generative models, which are briefly reviewed next.

Generative models have the ability to generate new data (karras2019style, )
and therefore the potential of addressing the data imbalance challenges in
bias. Generative methods –broadly speaking– learn to sample from the
underlying training data distribution so as to generate new samples who are
statistically indistinguishable from the training distribution. Those models
broadly fall into several categories that encompass: generative adversarial
networks (GANs) (karras2019style, ; grover2019bias, ), autoencoders,
variational autoencoders (VAEs) (kingma2013auto, ; louizos2015variational, ;
zhao2017infovae, ), generative autoregressive models, and invertible flow-
based latent vector models.

Generative methods have evolved and have culminated recently in approaches
leading to GANs that achieve realism on high resolution images (a former
limitation of GANs). This includes approaches such as BigGAN (brock2018large,
), which relies on SAGAN (self-attention GANs). (brock2018large, ) exploits
larger batch sizes which appear to improve performance. That study noted also
that larger networks has a comparable positive effect. Another positive effect
results from usage of the truncation method which consists of, for the
generator, sampling from a standard normal distribution in training, while
sampling instead from a truncated normal distribution at inference/generation
time, where samples above a certain threshold are re-sampled. Truncating with
a lower threshold allows to trade-off between higher fidelity and lower
diversity. While this model has resulted in many improvements to yield high-
resolution image synthesis, it did not aim for alteration of the image akin to
style transfer.

As an alternative, StyleGAN (karras2019style, ) used a multi-scale design that
was successful at addressing the generation of high resolution images (512x512
or 1024x1024), and allowed for stylistic image alterations. In that approach a
style vector W is employed to affect factors of the image (at the low scale,
coarse factors like skin tone, and at higher scales, fine factors like hair).
As such this method is able to perform stylistic mixing.

However all these methods have limitations for addressing the problem of data
imbalance we consider here because – while they generate realistic images –
these generative actions are essentially uncontrolled alterations to the
image. Even when the method allows style mixing (karras2019style, ) this
mixing does not allow specific control of a given factor in the latent space
$Z$ or style space $W$, without alteration of other concomitant (entangled)
factors. In sum, those methods, and derivative approaches including the method
in (burlina2019assessment, ) fail to generate images specifically for a
missing or under represented protected factor (eg retinal images for a
specific subpopulation, say African Americans with AMD lesions).

This motivates the search for other methods that allow fine control of
individual semantic images factors (e.g. images of dark skin individuals with
Lyme disease (burlina2020ai, )). While work has emerged that starts to address
that problem (paul2020unsupervised, ), the challenge is that when control of
attributes is achieved it may come at the cost of entanglement among latent
factors that controls those attributes. While disentanglement was linked to
fairness empirically in (Locatello2019OnTF, ), methods that achieve this
disentanglement are yet to be fully fleshed out. These limitations motivate
our approach here, which uses latent space manipulation to control specific
attributes while keeping other attributes unchanged.

In summary: fairness is achieved here via generation of more training data and
this generation is done in a way that finely controls attributes; This method
is able to perform debiasing when applied to the specific use case of age
related macular degeneration studied here.

##  2\. Methods

The general strategy is as follows: our method performs generation of
additional data with fine control of attributes. It does this for attribute
conditions and populations less or not represented originally in the training
dataset to aim to reduce dataset imbalance. For the specific use case of AMD
the goal is to generate more fundus images for underrepresented populations
(e.g. dark skin individuals with AMD). This explains why we desire to generate
data with fine control of attributes, e.g. change markers for ethnicity or
race as well as image markers for the disease, while leaving the rest of the
image and image markers unaffected.

These images are generated as follows: starting with healthy African American
individuals our goal is to generate new images that include lesions but keep
other image characteristics unchanged such as vasculature and possibly other
image markers like those described in (poplin2018prediction, ). The approach
uses StyleGAN as a basis. StyleGAN includes a fully connected network that
takes a latent vector $Z$ and transforms it into a 512-length intermediate
latent vector $W=\\{w_{i}\\}$ for all scales $i$ that controls a multiscale
generation process via an adaptive instance normalization (AdaIN) operation
and allows stylistic mixing of images (hence the so called ’style’ vector).
This style vector $W$ was originally used in StyleGAN to affect some factors
of the image at the low scale, coarse factors like skin tone, and at the
higher scale, fine factors like hair, by doing essentially style mixing.
Direct style mixing cannot be used for our purposes which require controlled
data generation since style mixing – along with some useful characteristics
changes added to the image (e.g. image lesions) – may impact other
characteristics and other image markers and changes that are undesirable (e.g.
changes in the vasculature).

Our process works as follows: it first trains StyleGAN on the available set of
training data, which is biased. It then uses a generated starter image sample
in latent space that corresponds to a retina of an African American
individual; we denote this by the tuple (image. latent vector) ($I$,
$\\{w_{i}\\}$), for all scales i in StyleGAN. Then it moves this sample along
a trajectory in latent space that takes its initial position to a location in
that space that entails it having more disease lesions (AMD), turning the
initial $\\{w_{i}\\}$ into $\\{w_{i}^{\prime}\\}$. This new latent space
vector $\\{w_{i}^{\prime}\\}$ sample is then turned into an actual image
$I^{\prime}$ via the generative model. This process is repeated for as many
samples that are needed to reestablish balance in the dataset. Once a balanced
dataset is generated a new AMD diagnostics model is then trained with the
augmented data and performance is compared to the originally biased model.

To achieve this trajectory in latent space a gradient descent direction of
motion is used. This gradient descent is used to maximize a loss function that
is obtained as follows: A simple classifier (a neural network) is trained on
the original images to distinguish AMD vs no AMD affected retinas. This
classifier is then used to classify synthetic samples of retinas and those
samples representations in latent space are used to train a new classifier
(AMD vs. no AMD) that now operates in latent space. This second classifier’s
output loss function, evaluated for a desired label that corresponds to having
AMD, is then used as loss function in the gradient decent process. The net
effect of the gradient decent is to ’impart’ the disease (AMD) to the original
image with minimal alteration of the rest of the image. Hence a fundus image
from a African American individual that is likely healthy with regard to AMD
is transformed into a realistic fundus retinal image of an African American
individual that is AMD-affected.

##  3\. Experiments

To measure the effect of the debiasing on an originally biased classifier for
AMD we look at several metrics principally: accuracy, sensitivity,
specificity, F1, kappa score, PPV and NPV, along with average precision and
ROC AUC.

Figure 1. Examples of African American fundus images, starting with healthy (left), transformed to affected (right), through our domain adaptation approach. Metrics | Baseline | New Method
---|---|---
| Method | Domain Adaptation
Accuracy | 71.75 (5.03) | 82.79 (4.22)
Sensitivity | 59.09 (7.77) | 74.03 (6.93)
Specificity | 84.42 (5.73) | 91.56 (4.39)
PPV | 79.13 (7.43) | 89.76 (5.27)
NPV | 67.36 (6.62) | 77.90 (6.04)
Weighted Kappa | 0.4351 | 0.6558
F1 | 0.6766 | 0.8114
Average Precision | 0.8469 (0.0402) | 0.9284 (0.0288)
ROCAUC | 0.8049 (0.0443) | 0.9169 (0.0308)
Test Set Subset Analysis: |  |
Accuracy (Caucasians) | 80.52 (6.26) | 85.71 (5.53)
Accuracy (African Americans) | 62.99 (7.63) | 79.87 (6.33)
Larger Leftover Set Analysis (614 images): |  |
Accuracy (Leftover dataset of African Americans/AMD) | 42.51 (3.91) | 54.23 (3.94)
Table 1. Results for debiasing comparing results for the baseline model and the synthetic debiased model using latent space manipulation showing performance improvement as well as debiasing for the proposed approach. Test data | Healthy | Affected (AMD)
---|---|---
Caucasians | 1843 | 1843
African Americans | 3686 | 0
Table 2. Characteristic table for training dataset for baseline AMD diagnostic model. Test data | Healthy | Affected (AMD)
---|---|---
Caucasians | 77 | 77
African Americans | 77 | 77
Table 3. Characteristic table for testing dataset. This dataset is used both
for the baseline as well as the new (domain adapted) diagnotics model.

We use a retinal disease use case focused on AMD that is based on using the
AREDS dataset. AREDS is a dataset (burlina2011automatic, ;
burlina2019assessment, ) made available by request to the NIH. It contains
individuals’ retinas with labels corresponding to a 4-scale AMD severity. It
also includes information on ethnicity. Individuals self declared as Caucasian
or African American are used in this study.

We built a model that predicts no AMD (severity 1 and 2) vs. referable AMD
(severity 3 and 4). We employ a use case of data imbalance where the training
partition was such that no data was made available to the baseline model for
African American individuals that were affected by AMD, all other types of
data were made available. The training dataset is described in the
characteristic table 2 for training. Note that an equal number of healthy and
diseased retinal fundi were used in training so as not to incur artificial
bias in the disease prediction overall. Only a reduced number of images were
left over for testing considering we had to equally balance those partitions
resulting in Table 3. A left over dataset of retinal images from African
American individuals with AMD was used for additional testing (about 614
images). Results are shown in 1 comparing the original baseline AMD
diagnostics model and the model. It is apparent looking at these results that
the baseline model trained on the originally unbalanced dataset exhibits bias
(Accuracy on Caucasian is 80.52% (6.26%) vs African American individuals,
62.99% (7.63%)) while the new model trained on the augmented dataset has
accuracies of 85.71% (5.53%) for Caucasian individuals vs 79.87% (6.33%) for
African American individuals.

##  4\. Discussion

The proposed method that generates more data via latent space manipulation and
is able to address domain adaptation shows performance improvements in regard
to most metrics (e.g. the baseline system has a ROC AUC of 0.8049 (0.0443) vs
an ROC AUC of 0.9169 (0.0308) for the new diagnotic model using augmented data
via latent space manipulation).

What is noticeable however is that the new diagnostic model is able to improve
overall performance as well as reduce the performance gap between populations.

Looking at the 95% confidence interval allows us to conclude positively that
the performance improvement in metrics such as accuracy and ROC AUC between
baseline and new system. Looking at the overall system performance also
indicates similar improvements when considering the other performance metrics.

Similarly, looking at accuracy and noting the difference between the accuracy
for Caucasian and African American individuals shows significance in terms of
the improvement made, from baseline model to new model; it shows that there
was definitely a bias in the original model; it also shows that this bias (gap
in accuracy between performance of African American individuals and Caucasian
individuals) became within confidence interval for the new diagnostic model,
showing that the new model was more fair. The metric measuring accuracy on a
leftover dataset of African American individuals with AMD also demonstrated
similar improvement in performance in the new model.

Future work will entail the development of novel metrics to better understand
the tradeoffs between accuracy and fairness which are not discussed here.

##  5\. Conclusion

We study a method for debiasing that generates more data for underrepresented
populations and addresses domain adaptation as a means of tackling bias. We
show that we are able to use this method to address AI fairness in a use case
that entails developing a DL based detection system for age related macular
degeneration. The experiments show that this approach is indeed able to
improve fairness.

## References

  * (1) Alex Beutel, Jilin Chen, Zhe Zhao, and Ed H Chi.  Data decisions and theoretical implications when adversarially learning fair representations.  arXiv preprint arXiv:1707.00075, 2017.
  * (2) Tolga Bolukbasi, Kai-Wei Chang, James Y Zou, Venkatesh Saligrama, and Adam T Kalai.  Man is to computer programmer as woman is to homemaker? debiasing word embeddings.  In Advances in neural information processing systems, pages 4349–4357, 2016.
  * (3) Andrew Brock, Jeff Donahue, and Karen Simonyan.  Large scale gan training for high fidelity natural image synthesis.  arXiv preprint arXiv:1809.11096, 2018.
  * (4) Philippe Burlina, David E Freund, Bénédicte Dupas, and Neil Bressler.  Automatic screening of age-related macular degeneration and retinal abnormalities.  In 2011 Annual International Conference of the IEEE Engineering in Medicine and Biology Society, pages 3962–3966. IEEE, 2011.
  * (5) Philippe Burlina, Neil Joshi, William Paul, Katia D Pacheco, and Neil M Bressler.  Addressing artificial intelligence bias in retinal disease diagnostics.  arXiv preprint arXiv:2004.13515 and Translational Vision Science and Technology (accepted), 2020.
  * (6) Philippe Burlina, William Paul, Philip Mathew, Neil Joshi, Katia D Pacheco, and Neil M Bressler.  Low-shot deep learning of diabetic retinopathy with potential applications to address artificial intelligence bias in retinal diagnostics and rare ophthalmic diseases.  JAMA ophthalmology, 2020.
  * (7) Philippe M Burlina, Neil Joshi, Katia D Pacheco, TY Alvin Liu, and Neil M Bressler.  Assessment of deep generative models for high-resolution synthetic retinal image generation of age-related macular degeneration.  JAMA ophthalmology, 137(3):258–264, 2019.
  * (8) Philippe M Burlina, Neil J Joshi, Phil A Mathew, William Paul, Alison W Rebman, and John N Aucott.  Ai-based detection of erythema migrans and disambiguation against other skin lesions.  Computers in Biology and Medicine, 125:103977, 2020.
  * (9) Nicholas Carlini and David Wagner.  Adversarial examples are not easily detected: Bypassing ten detection methods.  In Proceedings of the 10th ACM Workshop on Artificial Intelligence and Security, pages 3–14, 2017.
  * (10) Yaroslav Ganin, Evgeniya Ustinova, Hana Ajakan, Pascal Germain, Hugo Larochelle, François Laviolette, Mario Marchand, and Victor Lempitsky.  Domain-adversarial training of neural networks.  The Journal of Machine Learning Research, 17(1):2096–2030, 2016\.
  * (11) Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative adversarial nets.  In Advances in neural information processing systems, pages 2672–2680, 2014.
  * (12) Aditya Grover, Jiaming Song, Ashish Kapoor, Kenneth Tran, Alekh Agarwal, Eric J Horvitz, and Stefano Ermon.  Bias correction of learned generative models using likelihood-free importance weighting.  In Advances in Neural Information Processing Systems, pages 11058–11070, 2019.
  * (13) Moritz Hardt, Eric Price, and Nati Srebro.  Equality of opportunity in supervised learning.  In Advances in neural information processing systems, pages 3315–3323, 2016.
  * (14) Tero Karras, Samuli Laine, and Timo Aila.  A style-based generator architecture for generative adversarial networks.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 4401–4410, 2019.
  * (15) Diederik P Kingma and Max Welling.  Auto-encoding variational bayes.  arXiv preprint arXiv:1312.6114, 2013.
  * (16) Newton M Kinyanjui, Timothy Odonga, Celia Cintas, Noel CF Codella, Rameswar Panda, Prasanna Sattigeri, and Kush R Varshney.  Fairness of classifiers across skin tones in dermatology.  In International Conference on Medical Image Computing and Computer-Assisted Intervention, pages 320–329. Springer, 2020.
  * (17) Francesco Locatello, Gabriele Abbati, Tom Rainforth, S. Bauer, B. Schölkopf, and Olivier Bachem.  On the fairness of disentangled representations.  In NeurIPS, 2019.
  * (18) Christos Louizos, Kevin Swersky, Yujia Li, Max Welling, and Richard Zemel.  The variational fair autoencoder.  arXiv preprint arXiv:1511.00830, 2015.
  * (19) Ninareh Mehrabi, Fred Morstatter, Nripsuta Saxena, Kristina Lerman, and Aram Galstyan.  A survey on bias and fairness in machine learning.  arXiv preprint arXiv:1908.09635, 2019.
  * (20) William Paul, I-Jeng Wang, Fady Alajaji, and Philippe Burlina.  Unsupervised semantic attribute discovery and control in generative models.  arXiv preprint arXiv:2002.11169 and Neural Computation (accepted), 2020.
  * (21) Mike Pekala, Neil Joshi, TY Alvin Liu, Neil M Bressler, D Cabrera DeBuc, and Philippe Burlina.  Deep learning based retinal oct segmentation.  Computers in Biology and Medicine, 114:103445, 2019.
  * (22) Ryan Poplin, Avinash V Varadarajan, Katy Blumer, Yun Liu, Michael V McConnell, Greg S Corrado, Lily Peng, and Dale R Webster.  Prediction of cardiovascular risk factors from retinal fundus photographs via deep learning.  Nature Biomedical Engineering, 2(3):158, 2018.
  * (23) Novi Quadrianto, Viktoriia Sharmanska, and Oliver Thomas.  Discovering fair representations in the data domain.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 8227–8236, 2019.
  * (24) Reza Shokri, Marco Stronati, Congzheng Song, and Vitaly Shmatikov.  Membership inference attacks against machine learning models.  In 2017 IEEE Symposium on Security and Privacy (SP), pages 3–18. IEEE, 2017.
  * (25) Jiaming Song, Pratyusha Kalluri, Aditya Grover, Shengjia Zhao, and Stefano Ermon.  Learning controllable fair representations.  In The 22nd International Conference on Artificial Intelligence and Statistics, pages 2164–2173, 2019.
  * (26) Christina Wadsworth, Francesca Vera, and Chris Piech.  Achieving fairness through adversarial learning: an application to recidivism prediction.  In Conference on Fairness, Accountability, and Transparency in Machine Learning (FATML), 2018.
  * (27) Brian Hu Zhang, Blake Lemoine, and Margaret Mitchell.  Mitigating unwanted biases with adversarial learning.  In Proceedings of the 2018 AAAI/ACM Conference on AI, Ethics, and Society, pages 335–340, 2018.
  * (28) Shengjia Zhao, Jiaming Song, and Stefano Ermon.  Infovae: Information maximizing variational autoencoders.  arXiv preprint arXiv:1706.02262, 2017.


# Worldsheet from worldline

Umut Gürsoy  Guim Planella Planas  Institute for Theoretical Physics, Utrecht
University, 3584 CE Utrecht, The Netherlands

###### Abstract

We take a step toward a “microscopic” derivation of gauge-string duality. In
particular, using mathematical techniques of Strebel differentials and
discrete exterior calculus, we obtain bosonic string worldsheet action for
string embedded in d+1 dimensional asymptotically AdS space from multi-loop
Feynman graphs of quantum field theory of scalar matrices in d-dimensions in
the limit of diverging loop number. Our work is building on the program
started by ’t Hooft in 1974 ’t Hooft (1974), this time including the
holographic dimension which we show to emerge from the continuum of Schwinger
parameters of Feynman diagrams.

##  I Introduction

Despite being widely accepted and having passed many non-trivial tests, gauge-
gravity duality Maldacena (1998); Gubser  _et al._ (1998); Witten (1998) lacks
a satisfactory derivation as an equivalence between quantum gauge theory and
string theory. In particular, apart from certain low dimensional examples and
integrable class of quantum field theories, e.g. Kazakov  _et al._ (1985);
Ginsparg and Moore (1993); Eberhardt  _et al._ (2019); Gaberdiel and Gopakumar
(2013, 2011), we are unaware of any “microscopic” derivation of the string
world sheet directly from the gauge theory path integral. This derivation
would entail, in a concrete sense, completing the program started by ‘t
Hooft’s seminal paper ’t Hooft (1974), that is, to derive the string path
integral — this time including the holographic dimension — from the ribboned
Feynman graphs of gauge theories.

A promising starting point for this endeavor is Gopakumar’s observation
Gopakumar (2004a, b, 2005). See Aharony  _et al._ (2021); Alday  _et al._
(2023); Douglas  _et al._ (2011); Lee (2014); Heemskerk  _et al._ (2009);
Razamat (2010); Ooguri and Vafa (2002); Eberhardt  _et al._ (2019); Berkovits
and Vafa (2008) for an incomplete list of other examples and alternative
approaches. In particular, Gopakumar (2005) proposed to use Strebel
differentials to map the skeleton graphs of n-point functions of gauge
invariant operators in free CFTs to the worldsheet of a putative string
theory, see also Gopakumar and Mazenc (2022); Gaberdiel  _et al._ (2021).
Furthermore, in Domingo Gallegos  _et al._ (2022); Gallegos  _et al._ (2022) a
curious similarity between (discretized) open and closed strings and a
graphical representation of one- and two-trees that appear in Symanzik
polynomial representation of high-loop Feynman graphs, see e.g. Lam and Lebrun
(1969); Bogner and Weinzierl (2010); Weinzierl (2022); Dai and Siegel (2007);
Gallegos  _et al._ (2022), were observed and it was proposed that embedding of
the sting in the holographic dimension arises from a continuum limit of the
Schwinger parameters of graphs with large number of loops.

We merge these observations in this work and explicitly construct the string
worldsheet action that arises from the continuum limit of Feynman graphs of a
quantum field theory in d-dimensions. We show that the target space of the
string that correspond to the QFTs we consider in this work is asymptotically
AdSd+1.

##  II Generating function

Consider expansion of a QFT generating function in Feynman diagrams. We
restrict ourselves to QFTs of N$\times$N matrices that are massless scalars
with a generic polynomial type interaction. Contribution of a graph $\Gamma$
is given by

$\mathcal{Z}_{\Gamma}\propto\int\prod_{i}d^{d}X_{i}\prod_{e_{ij}\in\Gamma}G(X_{i},X_{j})\prod_{v_{i}\in\mathcal{S}_{\Gamma}}J(X_{i})$
(1)

with the proportionality constant including all factors of the couplings and
symmetry factors, $e_{ij}$ and $v_{i}$ denoting the edges and the vertices and
$J$ denoting the sources. We assume, following ’t Hooft (1974), that Feynman
graphs are put in ribboned form and grouped as discretized string worldsheets
with genus $g$. We first rewrite (1) in a form that allows us to read off the
corresponding (discretized) worldsheet action.

A key role is played by the Schwinger parametrization of the propagators
$G(X_{i}-X_{j})$ with Schwinger parameters $\sigma_{ij}$. If there are $m$
edges between a pair of vertices then we can merge them in an effective
propagator Gopakumar (2004b) written in terms of a single Schwinger parameter
as

$\displaystyle G(X_{i},X_{j})^{m}=$
$\displaystyle\frac{i^{m}\Gamma\left(\frac{d}{2}-1\right)^{m}}{2^{\frac{md}{2}}\Gamma\left(m\left(\frac{d}{2}-1\right)\right)}\cdot$
$\displaystyle\cdot\int_{0}^{\infty}d\sigma_{ij}\,\sigma_{ij}^{m\left(\frac{d}{2}-1\right)-1}e^{-\frac{\sigma_{ij}(X_{i}-X_{j})^{2}}{4}}\,.$
(2)

We assume vanishing masses throughout this paper—see Gallegos  _et al._ (2022)
for generalization of this expression to massive scalars. This procedure
constructs the so-called “skeleton graph” $\bar{\Gamma}$ of $\Gamma$ Gopakumar
(2004b) with all duplicated edges fused into a single one and each skeleton
edge assigned a length given by the corresponding $\sigma_{ij}$. The partition
function is then expressed as

$\mathcal{Z}[J]=\sum_{\bar{\Gamma}}\int
d[X]d[\sigma]\mathcal{P}_{\bar{\Gamma}}e^{-\frac{1}{4}\sum\sigma_{ij}(X_{i}-X_{j})^{2}}$
(3)

where $[X]$ and $[\sigma]$ denote the collection of all vertex positions and
graph metrics of the skeleton graph and $\mathcal{P}_{\bar{\Gamma}}$ is a
polynomial in $\sigma_{ij}$ which depends on the graph structure, coupling
constants and sources.

To determine it, we introduce auxiliary variables $\xi_{i}$ at the vertices
$v_{i}$ and represent each edge connected to a vertex by $i\xi_{i}$ so that
for each propagator between $v_{i}$ and $v_{j}$ we must include a factor of
$-\xi_{i}\xi_{j}$. It is straightforward to show that

$\displaystyle\mathcal{P}_{\bar{\Gamma}}=\frac{1}{n!}\left(i\prod_{i}U\left(-i\partial_{\xi_{i}}\right)\right)\prod_{\bar{e}_{ij}\in\bar{\Gamma}}\sum_{m}\tfrac{(-1)^{m}\Gamma\left(\frac{d}{2}-1\right)^{m}}{2^{\frac{md}{2}}\Gamma\left(m\left(\frac{d}{2}-1\right)\right)}\cdot$
$\displaystyle\left.\cdot\sigma_{ij}^{m\left(\frac{d}{2}-1\right)-1}\frac{(\Tr(\xi_{i}^{\dagger}\xi_{j}+\xi_{j}^{\dagger}\xi_{i}))^{m}}{2(m!)}\right|_{\xi=0}\,$
(4)

gives a non-vanishing contribution only when the number of edges connected to
the vertex $v_{i}$ is precisely equal to the number of derivatives and only
when the coordination number of the vertices equal to the one on some diagram
$\Gamma$ whose skeleton diagram is $\bar{\Gamma}$. The symmetry factors of the
graph are correctly accounted for by the factors of $1/n!$ and $1/m!$ in this
expression. The latter arises from swapping of propagators between the same
pair of vertices and the former from swapping of vertices. This gives the
correct factor as long as we sum over all nonequivalent permutations of
vertices. This sum will be implemented below as part of the integral over the
moduli space, which should now be the moduli space of Riemann surfaces with
$n$ marked points counting separately all nonequivalent permutations of these
points.

We recognize the exponential in (3) as a kinetic term

$S_{K}=\frac{1}{4}\sum_{ij}\sigma_{ij}(X_{i}-X_{j})^{2}\,,$ (5)

and rewrite $\mathcal{P}_{\bar{\Gamma}}$ above as an interaction potential

$S_{I}[\xi,\theta]=-\sum\log(iU(\theta_{i}))-i\xi_{i}\theta_{i}+h(\sigma_{ij},\xi_{i}\xi_{j})\,,$
(6)

where the function $h$ is

$\displaystyle h(\sigma_{ij},\xi_{i}\xi_{j})$
$\displaystyle=\log(\sum_{m}\tfrac{\Gamma\left(\frac{d}{2}-1\right)^{m}(-1)^{m}}{2^{\frac{md}{2}}\Gamma\left(m\left(\frac{d}{2}-1\right)\right)}\sigma_{ij}^{m\left(\frac{d}{2}-1\right)-1}(\xi_{i}\xi_{j})^{m})$
$\displaystyle=-\log\sigma_{ij}+f\left(\sigma_{ij}^{\frac{d}{2}-1}\xi_{i}\xi_{j}\right)\,.$
(7)

To obtain the expression (6) we followed the trick

$\left(\cdots\right)\bigg{|}_{\xi=0}=\int_{-\infty}^{\infty}d\xi\int_{-\infty}^{\infty}d\theta\,e^{-i\xi\theta}\left(\cdots\right)\,,$

and performed integration by parts inside the path integral.

##  III Triangulation and continuum

We view the sum of the kinetic and potential terms in (5) and (6) as
discretized worldsheet action of the dual string theory and take the continuum
limit to obtain the continuous worldsheet theory. This should be done with
care employing the technology of Strebel differentials and Discrete Exterior
calculus which we shortly review below, and in more detail in the
Supplementary material.

###  III.1 Strebel differentials

Theory of Strebel differentials Mulase and Penkava (1998) provides an
isomorphism between the space of metric ribboned graphs and the moduli space
on punctured Riemann surfaces with a positive real number assigned to every
puncture. To every propagator with length $\sigma$ corresponds a vertical
strip in the complex plane of width $\sigma$ Gopakumar (2005); Gopakumar and
Mazenc (2022). Then, these strips are glued into a manifold as depicted in
figure 1 and detailed in the Supplementary material. The vertices of $\Gamma$
become punctures in the manifold and a positive real number $a_{i}$ associated
to the $i$-th puncture is given by the sum of all the lengths of edges
incident to the corresponding vertex $a_{i}=\sum_{j}\sigma_{ij}$.

Figure 1: Gluing of the strips at (a) a face and (b) a vertex. In (c) we show
an example of the Riemann surface obtained from a simple graph of genus $0$
with $4$ vertices corresponding to a sphere with $4$ infinite cylinders
attached.

This construction naturally provides a metric for the Riemann surface given by
the flat metric on each strip, corresponding to a particular gauge fixing of
Weyl and diffeomorphism symmetries. This metric exhibits a double pole around
the punctures (vertices of $\Gamma$)

$g_{Str}=\frac{a_{i}^{2}}{4\pi^{2}}\frac{dud\bar{u}}{\left|u\right|^{2}}\,.$
(8)

###  III.2 Discrete Exterior Calculus

Theory of Strebel differentials provide the worldsheet as a Riemann surface
but it does not determine the worldsheet action from (5) and (6). For this,
one first needs to extend the map between double poles of the Riemann surface
(8) and location of vertices in space time, $X_{i}^{\mu}$, to entire functions
$X^{\mu}(u,\bar{u})$ on the Riemann surface. We can address this issue through
the machinery of Discrete Exterior Calculus used to solve geometric problems
on manifolds through numerical methods where the specified data is necessarily
finite Desbrun  _et al._ (2005).

Consider a Riemann surface $\Sigma$ with a triangulation, see below, in terms
of points $v_{i}$, intervals $e_{ij}$ and triangles $f_{ijk}$ such that the
edges of $f_{ijk}$ are $e_{ij}$, $e_{jk}$ and $e_{ki}$ and the endpoints of
$e_{ij}$ are $v_{i}$ and $v_{j}$. One discretizes data on the manifold by
integrating over the relevant simplexes of our triangulation. That is, one
integrates a function over the points $v_{i}$, a 1-form over the intervals
$e_{ij}$ and a 2-form over the triangles $f_{ijk}$. This allows us to use
Stokes’ theorem without issue. For instance, take the $1$-form $df$ specified
over the edges $e_{ij}$, then one has

$(df)_{e_{ij}}=\int_{e_{ij}}df=\int_{\partial e_{ij}}f=(f)_{j}-(f)_{i}\,.$ (9)

This defines for us a discrete exterior derivative which is unambiguous in the
discretization lengths by virtue of Stokes’ theorem.

Another crucial operation is the Hodge dual. We first define the dual
$\star\sigma^{k}$ of a $k$-simplex $\sigma^{k}$ by the convex span of the
circumcenters of all the simplexes that contain $\sigma^{k}$ and $\sigma^{k}$
itself, see Figure 2.

Figure 2: Depiction of part of the triangulation (a) with its dual shown in
(b). Cells painted in the same color are dual to each other.

The hodge dual of a $k$-form $\alpha$ is then defined by

$(\star\alpha)_{\star\sigma^{k}}=\frac{\left|\star\sigma^{k}\right|}{\left|\sigma^{k}\right|}(\alpha)_{\sigma^{k}}$
(10)

where $\left|\cdot\right|$ stands for the volume in the corresponding
dimension 111Note that $\star\alpha$ is not discretized in the same
triangulation as $\alpha$, but instead in a dual version of it.. The last
operation we need is the wedge product Desbrun  _et al._ (2005)

$(\alpha\wedge\beta)_{\sigma^{k}\cup\star\sigma^{k}}=\frac{\left|\sigma^{k}\cup\star\sigma^{k}\right|}{\left|\sigma^{k}\right|\left|\star\sigma^{k}\right|}(\alpha)_{\sigma^{k}}(\beta)_{\star\sigma^{k}}$
(11)

where $\sigma^{k}\cup\star\sigma^{k}$ is the smallest convex set containing
both $\sigma^{k}$ and $\star\sigma^{k}$. If we decrease the size of our
triangulation we can expect all these operations to converge Mohamed  _et al._
(2018); Schulz and Tsogtgerel (2020) since surfaces are locally flat and thus
well approximated by subsets of $\mathbb{R}^{2}$.

###  III.3 Triangulation

Finally, we need to provide a triangulation of the Riemann surface—which we
obtained using Strebel differentials above—on which discrete exterior calculus
can operate. This is tricky, as the punctures are located at the endpoints of
semi-infinite cylinders, see Fig. 1. We can, however, remove these cylinders
by a Weyl transformation with a factor of $\epsilon^{2}\left|u\right|^{2}$ in
the region $0<\left|u\right|<1$ on the Strebel metric (8). This leaves a flat
punctured disk which can be completed by restoring the marked point. Disks
corresponding to different poles are identified at the horizontal curves at
$\left|u\right|=1$ defining a continuous Weyl transformation on the entire
manifold. The resulting manifold is a collection of sections of flat disks
with opening angles $\frac{2\pi\sigma_{ij}}{a_{i}}$ and radii
$\frac{a_{i}\epsilon}{2\pi}$ identified along their curved edges, see Fig. 3.
The factor of $\epsilon$ enables us to adjust the volume $V$ of the manifold,
i.e. to keep it fixed even when we introduce more poles or change the $a_{i}$.
One finds

$V=\sum_{i}\frac{\epsilon^{2}a_{i}^{2}}{4\pi}\qquad\Rightarrow\qquad\epsilon=\sqrt{\frac{4\pi
V}{\sum a_{i}^{2}}}\sim\frac{1}{\sqrt{n}}\,,$ (12)

What we achieved so far is not quite a triangulation because sections of discs
are not simplexes, but we can split each section further into small triangles
that approximate the discs arbitrarily well. Note that at strong coupling,
where diagrams with $n\gg 1$ dominate the partition function, the mesh
describing our manifold becomes finer and finer.

We will need the following geometric data, which we can easily deduce from
Fig. 3. Denoting the poles by indices $i$ and $j$ and the auxiliary vertices
introduced to make the discs into triangles by an index $a$, we obtain
$\left|e_{ia}\right|=a_{i}\epsilon/2\pi$ for length of the auxiliary edges and
$\left|\star e_{ia}\right|=\epsilon\delta\sigma/2$ for their duals. Their
convex span has area $\left|e_{ia}\cup\star
e_{ia}\right|=a_{i}\delta\sigma\epsilon^{2}/8\pi$. The area of the dual of a
pole is $\left|\star v_{i}\right|=a_{i}^{2}\epsilon^{2}/16\pi$ and for the
auxiliary vertices lying between disk $i$ and disk $j$ we find $\left|\star
v_{a}\right|=3(a_{i}+a_{j})\delta\sigma\epsilon^{2}/16\pi$.

Figure 3: After we perform a Weyl transformation the infinite cylinders become
flat disks split into sections that are identified along its curved edge. We
triangulate the manifold by splitting this sections into small triangles.

###  III.4 Continuum

We are ready to take the continuum limit of the discrete action (5) and (6).
Discrete exterior derivative of $X^{\mu}$ is given by

$(\differential
X^{\mu})_{e_{ia}}=\frac{a_{i}}{a_{i}+a_{j}}\left(X_{j}-X_{i}\right)$ (13)

where the vertex $v_{a}$ lies between the disks $i$ and $j$ and we chose to
linearly interpolate the known values of $X^{\mu}$ to the auxiliary vertices
on the edges of the disk. Following the definitions above one finds

$(adX^{\mu}\wedge\star dX_{\mu})_{e_{ia}\cup\star e_{ia}+e_{ja}\cup\star
e_{ja}}=\frac{\pi\delta\sigma}{4}(X_{j}-X_{i})^{2}$ (14)

Summing over the entire manifold yields

$S_{K}=\frac{1}{\pi}\int
d^{2}z\sqrt{\textrm{det}g}\,a\,\nabla_{\alpha}X^{\mu}\nabla_{\beta}X^{\nu}g^{\alpha\beta}\eta_{\mu\nu}\,,$
(15)

up to corrections of order $\epsilon\sim 1/\sqrt{n}$.

The continuum limit of the potential term (6) is trickier, as we need to find
a function $\sigma$ such that

$(\sigma)_{v_{a}}=\sigma_{ij}$ (16)

so that we can express it in terms of the data on the manifold. Because the
Strebel differential is unique for given puncture locations and residues
$a_{i}$, this function should be expressible in them. To find it, we compute
the second functional derivative of $S_{K}$ directly from the continuum
expression and using the discretized form $\ref{kinetic})$. Equating the two
yields,

$a(z)-n\sigma(z)\int_{S_{z}}\star f_{p}=0$ (17)

where $\sigma(z)$ assigns to every point the $\sigma_{ij}$ of the strip it
belongs to and $S_{z}$ refers to the strip that contains $z$. The function
$f_{p}$ is defined through its Hodge dual by $\star
f_{p}=\frac{1}{n}\sum_{i}\delta(z-z_{i})$ giving the density of poles
normalized to $1$. To evaluate this integral we assume that $f_{p}$ has some
finite width and is distributed preserving the symmetry of the disk around
$v_{i}$, then we obtain

$\int_{S_{z}}\star
f_{p}=\frac{\sigma_{ij}}{na_{i}}+\frac{\sigma_{ij}}{na_{j}}=\left(\frac{2\sigma}{na}\right)_{v_{a}}\\!\\!\Rightarrow\quad\sigma(z)=\frac{a(z)}{\sqrt{2}}\,.$
(18)

The other quantity in (6) for which we need a continuum approximation for is
$\xi_{i}\xi_{j}$. It is straightforward to show that

$\xi_{i}\xi_{i}=\left(\xi^{2}+\frac{8\xi}{a\pi}\frac{\nabla\xi\cdot\nabla
a}{nf_{p}}-\frac{4}{\pi}\frac{(\nabla\xi)^{2}}{nf_{p}}\right)_{v_{a}}+\mathcal{O}\left(\frac{1}{n^{\frac{3}{2}}}\right)\,.$
(19)

Using the definition of $f_{p}$ to turn (6) into an integral we can express
its continuum limit as

$\displaystyle
S_{I}=-n\int\left(\log(iU(\theta))-i\xi\theta+h\left(\frac{a}{\sqrt{2}},\xi^{2}\right)\right)\star
f_{p}$
$\displaystyle\phantom{=}-\int\partial_{\xi^{2}}h\left(\frac{a}{\sqrt{2}},\xi^{2}\right)\left(\frac{8\xi}{a\pi}\nabla\xi\cdot\nabla
a-\frac{4}{\pi}(\nabla\xi)^{2}\right)\star 1$ (20)

up to $\epsilon\sim 1/\sqrt{n}$ corrections. For more details of this
computation we refer to Gursoy and Planella-Planas .

Putting everything together we find the following expression for the
generating function

$\mathcal{Z}[J]=\sum_{n}\frac{1}{n!}\int\mathcal{D}X^{\mu}\mathcal{D}\xi\mathcal{D}\theta\mathcal{D}a\mathcal{D}[g]\left|\frac{\delta\sigma}{\delta(a,g)}\right|e^{-S_{K}-S_{I}}$
(21)

where $S_{K}$ and $S_{I}$ are given by (15) and (III.4). The path integrals
are understood to be discretized on the triangulation described above and
$\mathcal{D}[g]$ refers to an integral over the moduli space of the Riemann
surfaces with marked points.

Leaving the Jacobian $\left|\delta\sigma/\delta(a,g)\right|$ aside for a
moment, we now consider the large $n$ limit of (21). Clearly the path
integrals over $\xi$ and $\theta$ will be dominated by the classical saddle of
(III.4). These can be further be solved in powers of $a$ around $a\to\infty$
which is also a natural limit to take, as it corresponds to the UV limit of
the QFT, which we expect to be a CFT—see Discussion below. Restricting to the
case of real fields with a single trace potential given by
$U=\Tr\left(\sum\frac{\lambda_{k}}{k!}\Phi^{k}+J^{T}\Phi\right)$, we find the
following saddles for $\xi$ and $\theta$

$\displaystyle\xi$
$\displaystyle=\frac{\alpha}{\sqrt{N}}\left(\frac{a}{\sqrt{2}}\right)^{-\frac{d-2}{4}}+\mathcal{O}\left(a^{-\frac{d}{4}}\right)$
(22) $\displaystyle\theta$
$\displaystyle=-\frac{ik\sqrt{N}}{\alpha}\left(\frac{a}{\sqrt{2}}\right)^{\frac{d-2}{4}}+\mathcal{O}\left(a^{\frac{d}{4}-1}\right)$
(23)

where $k$ is the order of the highest interaction in $U$ and $\alpha$ is a
numerical constant that satisfies

$2\alpha^{2}f^{\prime}(\alpha^{2})=k$ (24)

which has solutions for real $\alpha$ 222This holds only for $d\geq 3$
otherwise the terms that dominate the contribution from the potential are no
longer the highest order interactions.. More general situations can be
considered just as easily with similar results. Using this solution in (III.4)
we obtain

$S_{I}=C_{0}+C_{1}\int
d^{2}z\sqrt{\textrm{det}g}\frac{\nabla_{\alpha}a\nabla^{\alpha}a}{a^{2}}$ (25)

where the constant $C_{0}$ is undetermined and
$C_{1}=\left(1/4\pi+k(2d^{2}+7d-22)/16\pi\right)$ is not universal but depends
on our choice to discretize in (13). We further assume that $a(z)=1$ at all
the zeroes. Eqs. (25) and (15) constitute our final result for the worldsheet
action in the large $n$ and $a\to\infty$ limits. We find see that together
with an AdSd+1 sigma model
$ds^{2}=du^{2}+e^{-2u}\eta_{\mu\nu}dX^{\mu}dX^{\nu}$ with $a(z)=\exp(-2u(z))$
near the boundary. Deviations from this limit relate geometry toward the
interior to the RG flow of the QFT. We leave study of such deviations to
future work.

Consider finally the Jacobian in (21). Luckily, this was already obtained by
Kontsevich Kontsevich (1992), see also Dijkgraaf (1992); Mukhi (2005)

$\prod
d\sigma_{ij}=\left(\prod\frac{da_{j}}{4}\right)\times\frac{\Omega^{3g-3+n}}{(3g-3+n)!}2^{-5g+5}$
(26)

with $\Omega$ given by

$\Omega=\sum a_{i}^{2}\omega_{i}$ (27)

and $\omega_{i}$ being the first Chern class of the manifold evaluated at the
$i$-th marked point. While this contribution does not seem to affect our
conclusion above, it is desirable to provide its physical interpretation, see
the Discussion below.

##  IV Discussion

We note that the same $\omega_{i}$ in (27) appear in topological $2$d gravity
Gopakumar (2005). It turns out that the physical operators of this theory are
operators $\mathcal{O}_{n}$ indexed by the cohomology classes of the moduli
space $\mathcal{M}_{g,n}$. This suggest that one may replace insertions of
$\omega_{i}$ by insertions of $\mathcal{O}_{n}$. If that is the case, the non-
trivial measure can be written as a collection of topological terms while
simultaneously introducing a ghost sector to the theory as would be expected
from the fact that we are only integrating over the moduli space of metrics in
our derivation.

At first sight it is surprising that we obtain AdS target space from a generic
class of large N QFT in the large loop, large $a$ limit. There are two hidden
assumptions behind our result however. First, we ignored issues of
renormalization and assumed that the field mass continues to vanish in the
full quantum theory. One can show Gursoy and Planella-Planas  that issues of
renormalization do not alter our results, see also Gallegos  _et al._ (2022),
as long as the QFT is sufficiently close to a fixed point. Second, such a
fixed point is indeed found, precisely when the large loop limit is well
defined, see for example Zamolodchikov (1980). Indeed, Zamolodchikov
considered the so-called “fishnet” theories and showed that large vertex limit
exists at a specific coupling where the theory becomes a CFT for
$(d,k)=(3,6)$, $(4,4)$ and $(6,3)$, precisely the combinations argued to
possess a continuum limit in Gallegos  _et al._ (2022). In Basso and Zhong
(2019) the $(4,4)$ Fishnet theory has been shown to lead to string propagation
in AdS5 using integrability techniques. Our derivation provides an alternative
without resorting to integrability.

What happens beyond the CFT limit is of utmost importance. This can be
obtained by going beyond the large $a$ limit in the generic expressions (15),
(III.4). On the other hand, it is straightforward, yet cumbersome, to include
fermions and gauge fields in our derivation that is essentially based on the
Schwinger parametrization of the generating function, see for example Lam and
Lebrun (1969). Further generalizing our method to include renormalization
properly, one could, in principle be able to derive the worldsheet action dual
to QCD!

Finally, we note that our computation is valid for arbitrary genus. In
particular, (21) holds even though the Strebel metric varies yet remaining
unique. It is tempting to ask whether we can extend our analysis to non-
perturbative quantum gravity. In principle, our starting point, eqs. (5) and
(6) remains valid for arbitrarily small $N$. See Gallegos  _et al._ (2022) for
a discussion on this issue.

There are various other exciting questions to explore in the future, for
instance string propagation on blackholes through extending our analysis to
thermal field theory and dynamical gravity through, perhaps, generalizing to
ensemble of QFTs.

##  V Acknowledgments

We are grateful to Rajesh Gopakumar, Edward Mazenc and Claire Moran for
discussions. This work was supported, in part, by the Netherlands Organisation
for Scientific Research (NWO) under the VICI grant VI.C.202.104.

## References

  * ’t Hooft (1974) G. ’t Hooft, Nucl. Phys. B 72, 461 (1974).
  * Maldacena (1998) J. M. Maldacena, Adv. Theor. Math. Phys. 2, 231 (1998), arXiv:hep-th/9711200 .
  * Gubser  _et al._ (1998) S. S. Gubser, I. R. Klebanov, and A. M. Polyakov, Phys. Lett. B 428, 105 (1998), arXiv:hep-th/9802109 .
  * Witten (1998) E. Witten, Adv. Theor. Math. Phys. 2, 253 (1998), arXiv:hep-th/9802150 .
  * Kazakov  _et al._ (1985) V. A. Kazakov, A. A. Migdal, and I. K. Kostov, Phys. Lett. B 157, 295 (1985).
  * Ginsparg and Moore (1993) P. H. Ginsparg and G. W. Moore, in  _Theoretical Advanced Study Institute (TASI 92): From Black Holes and Strings to Particles_ (1993) pp. 277–469, arXiv:hep-th/9304011 .
  * Eberhardt  _et al._ (2019) L. Eberhardt, M. R. Gaberdiel, and R. Gopakumar, JHEP 04, 103 (2019), arXiv:1812.01007 [hep-th] .
  * Gaberdiel and Gopakumar (2013) M. R. Gaberdiel and R. Gopakumar, J. Phys. A 46, 214002 (2013), arXiv:1207.6697 [hep-th] .
  * Gaberdiel and Gopakumar (2011) M. R. Gaberdiel and R. Gopakumar, Phys. Rev. D 83, 066007 (2011), arXiv:1011.2986 [hep-th] .
  * Gopakumar (2004a) R. Gopakumar, Phys. Rev. D 70, 025009 (2004a), arXiv:hep-th/0308184 .
  * Gopakumar (2004b) R. Gopakumar, Phys. Rev. D 70, 025010 (2004b), arXiv:hep-th/0402063 .
  * Gopakumar (2005) R. Gopakumar, Phys. Rev. D 72, 066008 (2005), arXiv:hep-th/0504229 .
  * Aharony  _et al._ (2021) O. Aharony, S. M. Chester, and E. Y. Urbach, JHEP 03, 208 (2021), arXiv:2011.06328 [hep-th] .
  * Alday  _et al._ (2023) L. F. Alday, T. Hansen, and J. A. Silva, Phys. Rev. Lett. 131, 161603 (2023), arXiv:2305.03593 [hep-th] .
  * Douglas  _et al._ (2011) M. R. Douglas, L. Mazzucato, and S. S. Razamat, Phys. Rev. D 83, 071701 (2011), arXiv:1011.4926 [hep-th] .
  * Lee (2014) S.-S. Lee, JHEP 01, 076 (2014), arXiv:1305.3908 [hep-th] .
  * Heemskerk  _et al._ (2009) I. Heemskerk, J. Penedones, J. Polchinski, and J. Sully, JHEP 10, 079 (2009), arXiv:0907.0151 [hep-th] .
  * Razamat (2010) S. S. Razamat, JHEP 03, 049 (2010), arXiv:0911.0658 [hep-th] .
  * Ooguri and Vafa (2002) H. Ooguri and C. Vafa, Nucl. Phys. B 641, 3 (2002), arXiv:hep-th/0205297 .
  * Berkovits and Vafa (2008) N. Berkovits and C. Vafa, JHEP 03, 031 (2008), arXiv:0711.1799 [hep-th] .
  * Gopakumar and Mazenc (2022) R. Gopakumar and E. A. Mazenc, (2022), arXiv:2212.05999 [hep-th] .
  * Gaberdiel  _et al._ (2021) M. R. Gaberdiel, R. Gopakumar, B. Knighton, and P. Maity, JHEP 05, 073 (2021), arXiv:2011.10038 [hep-th] .
  * Domingo Gallegos  _et al._ (2022) A. Domingo Gallegos, U. Gürsoy, and N. Zinnato, J. Phys. Conf. Ser. 2191, 012012 (2022), arXiv:2202.07767 [hep-th] .
  * Gallegos  _et al._ (2022) D. Gallegos, U. Gursoy, and N. Zinnato, (2022), arXiv:2211.16514 [hep-th] .
  * Lam and Lebrun (1969) C. S. Lam and J. P. Lebrun, Nuovo Cim. A 59, 397 (1969).
  * Bogner and Weinzierl (2010) C. Bogner and S. Weinzierl, Int. J. Mod. Phys. A 25, 2585 (2010), arXiv:1002.3458 [hep-ph] .
  * Weinzierl (2022) S. Weinzierl,  _Feynman Integrals_ (2022) arXiv:2201.03593 [hep-th] .
  * Dai and Siegel (2007) P. Dai and W. Siegel, Nucl. Phys. B 770, 107 (2007), arXiv:hep-th/0608062 .
  * Mulase and Penkava (1998) M. Mulase and M. Penkava, “Ribbon graphs, quadratic differentials on riemann surfaces, and algebraic curves defined over $\bar{q}$,” (1998), arXiv:math-ph/9811024 [math-ph] .
  * Desbrun  _et al._ (2005) M. Desbrun, A. N. Hirani, M. Leok, and J. E. Marsden, “Discrete exterior calculus,” (2005), arXiv:math/0508341 [math.DG] .
  * Mohamed  _et al._ (2018) M. S. Mohamed, A. N. Hirani, and R. Samtaney, International Journal for Computational Methods in Engineering Science and Mechanics 19, 194 (2018).
  * Schulz and Tsogtgerel (2020) E. Schulz and G. Tsogtgerel, Discrete & Computational Geometry 63, 346 (2020).
  * (33) U. Gursoy and G. Planella-Planas, “to appear,” .
  * Kontsevich (1992) M. Kontsevich, Commun. Math. Phys. 147, 1 (1992).
  * Dijkgraaf (1992) R. Dijkgraaf, “Intersection theory, integrable hierarchies and topological field theory,” (1992), arXiv:hep-th/9201003 [hep-th] .
  * Mukhi (2005) S. Mukhi, “Topological matrix models, liouville matrix model and c=1 string theory,” (2005), arXiv:hep-th/0310287 [hep-th] .
  * Zamolodchikov (1980) A. B. Zamolodchikov, Phys. Lett. B 97, 63 (1980).
  * Basso and Zhong (2019) B. Basso and D.-l. Zhong, JHEP 01, 002 (2019), arXiv:1806.04105 [hep-th] .

##  Appendix A Supplementary material: Strebel Differentials

Let $\Sigma$ be a Riemann surface of genus $g$ with $n$ marked points
$\left\\{p_{1},\dots,p_{n}\right\\}$ equipped with a complex structure. The
basic objects we are interested in are quadratic meromorphic forms, that is,
objects $\chi\in\Omega^{(1,0)}_{1}(\Sigma)\otimes\Omega^{(1,0)}_{1}(\Sigma)$
which in holomorphic coordinates locally take the form

$\chi=\phi(z)\differential z^{2}$ (28)

with $\phi(z)$ being a holomorphic function that transforms as
$\phi(w)=\left(\frac{\differential w}{\differential z}\right)^{2}\phi(z)$
under holomorphic changes of coordinates. Using these differentials one can
define a special set of curves $\gamma$ such that
$\chi(\dot{\gamma},\dot{\gamma})>0$. We call such curves horizontal since for
$\Sigma=\mathcal{C}$ with the usual metric and choosing $\chi=\differential
z^{2}$ they correspond to lines of constant imaginary part. We similarly
denote by vertical curves those for which $\chi(\dot{\gamma},\dot{\gamma})<0$.
At regular points of $\phi$ the vertical and horizontal curves form a grid
which specify a set of coordinates for which $\chi=\differential z^{2}$. Near
zeroes or poles of $\phi$ this is no longer true and we depict the local
structure of vertical and horizontal lines in Figure 4. Note that from a
quadratic differential one can define a metric for the manifold given by

$g_{Str}=|\phi(z)|\differential z\differential\bar{z}$ (29)

In the special coordinates defined by vertical and horizontal curves this
metric reduces to patches of the Euclidean metric in $\mathcal{C}$ along with
curvature localized at the zeroes and poles of $\phi$.

Figure 4: Structure of a Strebel differential near a double pole (a) and a
single zero (b). The darker lines correspond to vertical curves and the
lighter ones to horizontal curves.

The main result of the theory of Strebel differentials is given by the
following theorem

###### Theorem 1

Given a Riemann surface $\Sigma$ of genus $g$ with $n$ marked points
$\left\\{p_{1},\dots,p_{n}\right\\}$ and a set of $n$ positive real numbers
$\left\\{a_{1},\dots,a_{n}\right\\}$ one can find a unique quadratic
meromorphic form $\chi$ such that

  * •

$\chi$ locally takes the form $\phi(z)\differential z^{2}$ with $\phi(z)$
holomorphic on $M/\left\\{p_{1},\dots,p_{n}\right\\}$

  * •

$\phi(z)$ has a double pole at every $p_{i}$

  * •

The union of all non-compact horizontal curves is a set of measure $0$

  * •

Every compact horizontal curve $\gamma$ is a simple loop around some $p_{i}$
and it satisfies

$a_{i}=\int_{\gamma}\sqrt{\phi(z)}\differential z$ (30)

Such a quadratic form is called a Strebel differential and the set of non-
compact horizontal curves is called its critical graph. These corresponds to
horizontal lines connecting two zeroes such that the critical graph can be
regarded as the embedding of some abstract graph in our manifold.

In light of theorem 1 one can deduce that Strebel differentials parametrize
the space $\mathcal{M}_{g,n}\times\mathbb{R}^{n}_{+}$ where
$\mathcal{M}_{g,n}$ is the moduli space of Riemann surfaces of genus $g$ with
$n$ punctures and one can take $g_{Str}$ as a representative of the classes of
metrics up to Weyl and diffeomorphism transformations for some fixed
collection of residues at every marked point.

From a metric ribboned graph $\Gamma$ one can construct a Riemann surface with
its corresponding Strebel differential such that its critical graph is given
by $\Gamma^{*}$, that is, the dual graph of $\Gamma$ obtained by replacing
vertices by faces and vice-versa. The construction goes as follows.

Let $\Gamma$ be an arbitrary metric ribboned graph with all faces bounded by
at least three distinct edges. For every edge $e_{i}$ of $\Gamma$ with length
$\sigma_{i}$ let us take an infinite strip in $\mathbb{C}$ oriented along the
imaginary axis and of width $\sigma_{i}$. On these strips $\differential
z_{i}^{2}$ defines a holomorphic quadratic form. Take any vertex $v_{i}$.
Since $\Gamma$ is ribboned we have an ordering of the edges incident to this
vertex, then we can glue the corresponding strips into a cylinder according to
this ordering. This can be done by mapping all these strips into a single disk
according to the map

$u_{i}(z_{k})=\exp\left(2\pi
i\frac{z_{k}+\sum_{j<k}\sigma_{j}}{\sum_{j}\sigma_{j}}\right)$ (31)

Pulling back the quadratic differentials on the strips to the $u$-disk we
obtain

$\differential
z_{k}^{2}=-\frac{(\sum\sigma_{j})^{2}}{4\pi^{2}}\frac{\differential
u_{i}^{2}}{u^{2}_{i}}$ (32)

so that the form develops a second order pole at the points $u_{i}=0$. Here
$k$ runs over all the vertices incident to $v_{i}$ according to the ordering
induced by the ribbon structure of $G$. We also do something similar for the
faces. Take any face $f_{i}$ with $m\geq 3$ edges bounding it. We now glue the
strips of the edges forming the loop into a disk by mapping them with

$w_{i}(z_{k})=e^{\frac{2\pi ik}{m}}z_{k}^{\frac{2}{m}}$ (33)

Now the quadratic differentials are mapped to

$\differential z_{k}^{2}=\frac{m^{2}}{4}w_{i}^{m-2}\differential w_{i}^{2}$
(34)

so that they develop a zero of order $m-2$ at the points $w_{i}=0$. The
collection of the strips with coordinates $z_{i}$ for the edges, the disks
with coordinates $u_{i}$ for the vertices and the disks with coordinates
$w_{i}$ for the faces define an atlas with holomorphic transition functions
given by $u_{i}(z_{k})$ and $w_{i}(z_{k})$. Thus they completely characterize
a Riemann surface with the same genus as the graph. Moreover the form $\chi$
with local expression on the strips given by $\differential z_{i}^{2}$ is the
unique Strebel differential with marked points given by $u_{i}=0$ and residues
$\sum_{e_{j}\in v_{i}}\sigma_{j}\equiv a_{i}$. In Figure 1 we give a drawing
of how the strips are glued at vertices and faces and an example of the
resulting Riemann surface from a simple graph.

There is no obstacle to following the same construction for diagrams that have
two-edge loops, however the strips corresponding to the two edges get glued
together into a single strip without any zero and the information that they
were originally separate is lost.

This construction is actually an isomorphism between the set of skeleton
ribboned graphs with lengths assigned to each edge with
$\mathcal{M}_{g,n}\times\mathbb{R}_{+}^{n}$ with the positive numbers $a_{i}$
being the coordinates on $\mathbb{R}_{+}^{n}$. Note additionally that we can
construct a triangulation of our Riemann surface with $0$-cells given by the
zeroes and the poles, $1$-cells given by the horizontal curves connecting
zeroes to zeroes and the vertical curves connecting zeroes to poles and
$2$-cells given by half-strips split at $\imaginary(z_{i})=0$.

Finally in the main text we have used some properties of Strebel differentials
for large number of punctures that we now prove. In particular, directly using
the solution for the saddles (22-23) on equation (III.4) a term of the form
$n\int\log a\star f_{p}$ appears that contributes to the kinetic term for
$a(z)$ by using the following identity

$-4\pi\left(nf_{p}-\frac{m}{2}f_{q}\right)=\nabla^{2}\log a$ (35)

where $m$ is the number of zeroes of the Strebel differential (including
multiplicities) and $f_{q}$ is the density of zeroes defined similarly to
$f_{p}$ replacing poles by zeroes. Note that the term with $f_{q}$ vanishes
when integrated against $loga$ due to the assumptions that $a(z)=1$ at the
zeroes. To show that this is true when we have a large number of punctures we
note that the holomorphicity of the Strebel differential implies

$\sqrt{\textrm{det}g_{0}}\;\nabla_{0}^{2}\log{\left|\phi(z)\right|}=-4\pi\sum_{p\in\text{poles}}\delta^{(2)}(z-p)+2\pi\sum_{q\in\text{zeroes}}\delta^{(2)}(z-q)$
(36)

This is all expressed in the metric given by $g_{0}=\differential
z\differential\bar{z}$ where $z$ are the coordinates that make the Strebel
metric take the form of equation (29). We would like to use instead the metric
$g$ where our triangulation has been defined which is related to $g_{0}$ and
$g_{Str}$ by a Weyl transformation. In particular we have

$\sqrt{\textrm{det}g}=e^{2\omega}\left|\phi(z)\right|\sqrt{\textrm{det}g_{0}}=e^{2\omega}\sqrt{\textrm{det}g_{Str}}$
(37)

where $\omega$ is the parameter of some Weyl transformation that we need to
find. Using the transformation properties of the Ricci scalar under Weyl
transformations we can find

$\nabla^{2}\omega=-\frac{R}{2}+2\pi\left(nf_{p}(z)-\frac{mf_{q}(z)}{2}\right)$
(38)

Which can be directly solved using the Green’s function to find

$\left|\phi(x,y)\right|=\left|C\right|\exp(-4\pi\int\left(nf_{p}(z^{\prime})-\frac{m}{2}f_{q}(z^{\prime})\right)G(z^{\prime};z)\differential^{2}z^{\prime})$
(39)

for some complex constant $C$. This must have the right residue at each of the
poles and imposing such a condition is the way to find $f_{q}$. Then, if we
want to compute the residue at some pole $p_{i}$ we can restrict the integral
that computes $\log\left|\phi\right|$ to integrate only over
$\left|z-p_{i}\right|>\epsilon$ and take the limit $\epsilon\to 0$ this
excludes precisely the contribution that gives the pole of interest and
nothing more. Demanding that the absolute value of $\phi$ reduces to the
residue when we remove the contribution from a pole we obtain the relation

$\lim_{\epsilon\to
0}|C|\exp(-4\pi\int_{\left|z^{\prime}-p_{i}\right|>\epsilon}\left(f_{p}(z^{\prime})-\frac{1}{2}f_{q}(z^{\prime})\right)G(z^{\prime};p_{i})\differential^{2}z^{\prime})=\frac{a(p_{i})}{2\pi}$
(40)

Assuming that the triangulation is dense enough for $f_{p}$ and $f_{q}$ to be
sufficiently continuous we directly obtain equation (35) by taking a logarithm
and a Laplacian.