Department of Physics University of Cambridge  Dr. Alpha Lee 0.25 Doctor of
Philosophy Wolfson College August 2022 LaTeX

# Applications of Gaussian Processes at Extreme Lengthscales: From Molecules
to Black Holes

Ryan-Rhys Griffiths

###### Abstract

In many areas of the observational and experimental sciences data is scarce.
Observation in high-energy astrophysics is disrupted by celestial occlusions
and limited telescope time while laboratory experiments in synthetic chemistry
and materials science are both time and cost-intensive. On the other hand,
knowledge about the data-generation mechanism is often available in the
experimental sciences, such as the measurement error of a piece of laboratory
apparatus.
Both characteristics make Gaussian processes (gps) ideal candidates for
fitting such datasets. gps can make predictions with consideration of
uncertainty, for example in the virtual screening of molecules and materials,
and can also make inferences about incomplete data such as the latent emission
signature from a black hole accretion disc. Furthermore, gps are currently the
workhorse model for Bayesian optimisation, a methodology foreseen to be a
vehicle for guiding laboratory experiments in scientific discovery campaigns.
The first contribution of this thesis is to use gp modelling to reason about
the latent emission signature from the Seyfert galaxy Markarian 335, and by
extension, to reason about the applicability of various theoretical models of
black hole accretion discs. The second contribution is to deliver on the
promised applications of gps in scientific data modelling by leveraging them
to discover novel and performant molecules. The third contribution is to
extend the gp framework to operate on molecular and chemical reaction
representations and to provide an open-source software library to enable the
framework to be used by scientists. The fourth contribution is to extend
current gp and Bayesian optimisation methodology by introducing a Bayesian
optimisation scheme capable of modelling aleatoric uncertainty, and hence
theoretically capable of identifying molecules and materials that are robust
to industrial scale fabrication processes.

###### keywords:

LaTeX PhD Thesis Engineering University of Cambridge

This thesis is the result of my own work and includes nothing which is the
outcome of work done in collaboration except as declared in the Preface and
specified in the text. I further state that no substantial part of my thesis
has already been submitted, or, is being concurrently submitted for any such
degree, diploma or other qualification at the University of Cambridge or any
other University or similar institution except as declared in the Preface and
specified in the text. It does not exceed the prescribed word limit for the
relevant Degree Committee

###### Acknowledgements.

I would like to thank all those who contributed in some part, indirect or
otherwise, to my research productivity over the past years. I would like to
thank Alpha Lee, my supervisor, firstly for giving me an opportunity to return
from the working world to pursue research at a time when it seemed like all
doors had been closed, and secondly, for striking a balance between giving me
the freedom to explore new topics and offering excellent guidance and support
when needed.
In rough chronological order, I would like to thank Philippe Schwaller who
introduced me to Alpha following a chance encounter on the streets of
Cambridge in April 2018. My life may have been very different save for that
meeting! Additionally, I would like to thank Philippe for his ongoing
collaboration and sharing his expertise in all things involving sequence data
and chemical reactions.
From my time at Prowler.io (now Secondmind Labs) from 2017-2018, I would like
to thank Alexis Boukouvalas for acting as a fantastic mentor and supporting me
in all my endeavours. I learned a great deal about both machine learning and
software engineering during our pair programming sessions, especially when
developing the code for adaptive sensor placement (2019_Grant). I would like
to thank James Hensman, Richard Turner and Carl Rasmussen for giving lectures
on Gaussian processes which sparked my interest in the topic. I would also
like to thank Adithya Devraj for the numerous interesting conversations about
machine learning and the philosophy of science.
I would like to thank my colleagues from Cambridge Spark for keeping me up to
speed with machine learning in industry from 2017-2022. In particular, Raoul
Gabriel-Urma, Petar Velickovic, Tim Hillel, Patrick Short, Catalina Cangea,
Sahan Bulathwela, Ilyes Khemakhen, Chris Davis, Fred Hallgren and Kevin
Lemagnen. Acting as a mentor for the Schmidt Data for Science Residency
program was a highlight where I had the opportunity to learn about areas
ranging from synthetic biology to geophysics and climate modelling.
I would like to thank my colleagues from the Lee group at TCM from 2018-2022,
namely, Philip Verpoort, Alex Aldrick, Penelope Jones, Yunwei Zhang, William
McCorkindale, Felix Faber, Alwin Bucher, Rhys Goodall, Janosh Riebesell, Rokas
Elijosius, David Kovacs and Emma King-Smith for social interactions and
academic discussions when I was not absent due to the global pandemic or
internships.
I would like to thank Bingqing Cheng, whom I first met in 2019 for involving
me in her work on the ASAP library (2020_Cheng), and for taking the time to
introduce me to a broad network of researchers applying machine learning to
problems in physics and materials science.
At the Institute of Astronomy I would like to thank Jiachen Jiang for
introducing me to the world of high-energy astrophysics in the summer of 2019
and for providing excellent guidance on the contents of Chapter 3, namely
modelling the multiwavelength variability of Mrk-335 using Gaussian processes
(2021_Mrk). I would also like to thank Douglas Buisson, Dan Wilkins and Luigi
Gallo for their feedback as well as Andy Fabian and Christopher Reynolds for
being kind enough to attend an astrophysics talk given by a PhD student
(myself) with no background in astrophysics!
At the Computer Lab I would like to thank Ben Day, Simon Mathis and Arian
Jamasb, whom I first met in 2020, for discussions on graphs, molecules,
proteins and antibodies. In particular, I would like to thank Arian for our
almost daily Slack conversations and for answering my endless lists of
questions!
I would like to thank my colleagues at Huawei Noah’s Ark Lab whom I began
working with in October 2020. I would like to thank Haitham Bou-Ammar for his
mentorship as well as Rasul Tutunov, Vincent Moens, Alexander Cowen-Rivers,
Alexander Maraval, Antoine Grosnit and Hang Ren whom I learned a great deal
from during our joint work (2020_Grosnit; 2020_Rivers; 2021_Grosnit).
I would like to thank Anthony Bourached, George Cann, Gregory Kell and David
Stork for their collaboration applying machine learning to artwork starting in
late 2020 (2021_Bourached_art; 2021_Cann; 2021_Stork; 2022_Kell). In
particular David has been an excellent mentor on the subject of research
practices.
I would like to thank Miguel Garcia-Ortegon, Vidhi Lalchand, James Wilson and
Luke Corcoran for informative discussions on heteroscedastic Bayesian
optimisation (2021_Griffiths), the topic of Chapter 6. I would like to thank
Ajmal Aziz and Edward Kosasih at the Institute for Manufacturing for
introducing me to supply chain logistics in the summer of 2021, and in
particular to applications of graph neural networks for supply chain problems
(2021_Aziz).
I would like to thank Jian Tang for supervising me at MILA from January 2022
as well as Bojana Rankovic, Sang Truong, Leo Klarner, Aditya Ravuri, Yuanqi
Du, Julius Schwartz, Austin Tripp, Alex Chan, Jacob Moss, Felix Opolka and
Chengzhi Guo for their contributions to the GAUCHE library.
I would like to thank Jake Greenfield for providing his photoswitch expertise
and in particular, for helping out in a tight spot by rediscovering an old
batch of lost molecules during a laboratory cleanup after the new batch had
been mislaid by courier following a 3-month journey through customs.
I would like to thank Henry Moss, whom I met virtually during the pandemic in
the summer of 2020 and with whom I began developing a Gaussian process library
for chemistry in the form of FlowMO (2020_flowmo). The evolved version,
GAUCHE, comprises the contents of Chapter 4. I also appreciate the daily Slack
discussions about all things to do with Bayesian optimisation and Gaussian
processes.
I would like to thank David Ginsbourger for hosting myself and Henry Moss in
Bern in April 2022 to discuss extensions of the ideas comprising Chapters 4
and 5 with Athénaïs Gautier and Anna Broccard. I would like to thank Ekansh
Verma and Souradip Chakraborty for involving me in their work on invariances
in Bayesian optimisation (2021_Verma) in the summer of 2022. I would like to
thank my long-time friend S. F. for rigorously inspecting the notation of the
final thesis in July 2022. I would like to thank Victor Prokhorov for always
providing interesting food for thought during our many conversations in
Cambridge.
On a personal note, I would like to thank Leandro Charanga and Monika
Jankauskaite for teaching me how to dance, Subhankar, Thomas C and Dan for
their advice on ethical dilemmas, Thomas M for discussions on mathematics and
Bachata, Teja for trying to teach me some gymnastics and my parents for their
ongoing support.

###### Table of Contents

  1. 1 Introduction
    1. 1.1 Motivation
    2. 1.2 Overview and Contributions
    3. 1.3 List of Publications
    4. 1.4 List of Software
  2. 2 Background
    1. 2.1 Gaussian Processes
      1. 2.1.1 Kernels
      2. 2.1.2 Predictions
      3. 2.1.3 Training
      4. 2.1.4 Bayesian Model Selection
    2. 2.2 Bayesian Optimisation
      1. 2.2.1 Black-Box Optimisation
      2. 2.2.2 Solution Methods
      3. 2.2.3 The Bayesian Optimisation Algorithm
      4. 2.2.4 The Surrogate Model
      5. 2.2.5 The Acquisition Function
  3. 3 Modelling Black Hole Signals with Gaussian Processes
    1. 3.1 Background on High-Energy Astrophysics
      1. 3.1.1 Black Holes
      2. 3.1.2 Active Galactic Nuclei
      3. 3.1.3 Seyfert Galaxies
      4. 3.1.4 Space Observatories
      5. 3.1.5 Accretion Disc Models
      6. 3.1.6 Markarian 335
    2. 3.2 Preface
    3. 3.3 Introduction
    4. 3.4 Modelling Markarian 335
      1. 3.4.1 Identifying the Flux Distribution
      2. 3.4.2 Noise Considerations
      3. 3.4.3 Lightcurve Simulations
      4. 3.4.4 Modelling Markarian 335 with Gaussian Processes
    5. 3.5 Structure Function Analysis
    6. 3.6 Lag and Coherence
    7. 3.7 Conclusions
  4. 4 Modelling Molecules with Gaussian Processes
    1. 4.1 Background on Molecular Machine Learning
      1. 4.1.1 Molecular Property Prediction
      2. 4.1.2 Chemical Reaction Optimisation
    2. 4.2 Preface
    3. 4.3 Introduction
    4. 4.4 Related Work
    5. 4.5 Molecular Kernels
      1. 4.5.1 Fingerprint Kernels
      2. 4.5.2 String Kernels
      3. 4.5.3 Graph Kernels
    6. 4.6 Experiments
      1. 4.6.1 Regression
      2. 4.6.2 Uncertainty Quantification (UQ)
      3. 4.6.3 Chemical Reaction Yield Prediction
      4. 4.6.4 Bayesian Optimisation
    7. 4.7 Conclusions
  5. 5 Molecular Discovery with Gaussian Processes
    1. 5.1 Preface
    2. 5.2 Introduction
    3. 5.3 Dataset Curation
    4. 5.4 Machine Learning Prediction Pipeline
      1. 5.4.1 Multioutput Gaussian Processes (MOGPs)
    5. 5.5 TD-DFT Performance Comparison
    6. 5.6 Human Performance Benchmark
    7. 5.7 Screening for Novel Photoswitches using the MOGP
      1. 5.7.1 Screening Criteria
      2. 5.7.2 Lead Candidates
    8. 5.8 Conclusions
  6. 6 Modelling Experimental Noise with Gaussian Processes
    1. 6.1 Preface
    2. 6.2 Introduction
    3. 6.3 Related Work
    4. 6.4 Heteroscedastic Bayesian Optimisation
      1. 6.4.1 The Most Likely Heteroscedastic Gaussian Process
      2. 6.4.2 Bayesian Optimisation with Aleatoric Noise Penalisation
      3. 6.4.3 Heteroscedastic Acquisition Functions
    5. 6.5 Experiments on Robustness to Aleatoric Uncertainty
      1. 6.5.1 Implementation
      2. 6.5.2 Heteroscedastic Sin Wave Function
      3. 6.5.3 Heteroscedastic Branin-Hoo Function
      4. 6.5.4 Soil Phosphorus Fraction Optimisation
      5. 6.5.5 Molecular Hydration Free Energy Optimisation
      6. 6.5.6 Heteroscedastic Acquisition Function Hyperparameters
      7. 6.5.7 Summary of Robustness Experiments
    6. 6.6 Ablation Study
      1. 6.6.1 Noise Properties
      2. 6.6.2 Initialisation Set Size
      3. 6.6.3 Summary of Ablation Experiments
    7. 6.7 Conclusions
  7. 7 Conclusion
    1. 7.1 Summary of Contributions
    2. 7.2 Future Work
      1. 7.2.1 GPs in Astrophysics
      2. 7.2.2 GPs in Chemistry
      3. 7.2.3 GPs in Scientific Experiments
  8. A Modelling Black Hole Signals with Gaussian Processes
    1. A.1 Additional Graphical Tests for Identifying the Flux Distribution
    2. A.2 Spectral Properties of the Examined Kernels
  9. B Modelling Experimental Noise with Gaussian Processes
    1. B.1 Heteroscedasticity of the Soil Phosphorus Fraction Dataset
    2. B.2 Additional Ablation Experiments
      1. B.2.1 Goldstein-Price Function
      2. B.2.2 Branin-Hoo Function
    3. B.3 Performance Impact of the Kernel Choice
  10. C Molecular Discovery with Gaussian Processes
    1. C.1 Sources of Experimental Data
    2. C.2 Dataset Visualisations
    3. C.3 Further Experiments
      1. C.3.1 Property Prediction
      2. C.3.2 Prediction Error as a Guide to Representation Selection
      3. C.3.3 Impact of Dataset Choice
      4. C.3.4 Human Performance Benchmark
      5. C.3.5 Confidence-Error Curves
      6. C.3.6 TD-DFT Benchmark
    4. C.4 Further Screening Details
      1. C.4.1 UV-Vis Absorption Spectroscopy
      2. C.4.2 Photoswitching
    5. C.5 Novelty of Screened Candidates relative to The Photoswitch Dataset

##  Chapter 1 Introduction

###  1.1 Motivation

The past decade has seen deep learning models achieve breakthroughs in
computer vision (2012_Krizhevsky), speech recognition (2013_Graves), and
natural language processing (2017_Vaswani). In fact, progress on developing
deep learning architectures has proceeded so rapidly that, as of 2021, machine
learning pioneer Andrew Ng has voiced the opinion that research into improving
deep architectures has plateaued, at least in the traditional domains of
vision, speech and language. Ng is now calling for a shift in focus towards
data-centric AI, arguing that the dataset, as opposed to the model, is now the
performance bottleneck in many real-world problems (2021_Ng).

In the natural sciences, however, model development is by no means a solved
problem. Large scientific datasets have existed for some time, such as those
generated by the Large Hadron Collider at CERN (2013_Cern), or the Chemical
Universe Database, GDB-17 (2012_Ruddigkeit), which enumerates 166 billion
small molecules. Developing effective models for scientific data is still an
active and fast-moving field of research however (2022_Kalinin). In contrast
to artificial data such as images, speech, and text, scientific data can often
be inexorably tied to causal paradigms, entailing challenges for purely data-
driven approaches seeking to achieve strong out-of-distribution (OOD)
performance. Lines of inquiry in this direction include incorporating
invariances due to symmetries into deep learning models for proteins and
molecules (2021_Jumper; 2020_Hermann), as well as causal mechanisms for
problems in physics (2021_Scholkopf).

A further challenge for building performant machine learning models for
scientific applications stems from the availability of data. While large
datasets in the sciences have undoubtedly been a key driver of research, there
are also areas of scientific discovery which will always be limited to small
data. Examples include molecular design, where one wishes to predict the
properties of a new class of molecule for which few experimental measurements
exist, as well as high-energy astrophysics, where one wishes to draw
inferences from astronomical time series with short observation periods. In
the past years researchers have achieved success in porting breakthroughs in
deep learning to large scientific datasets (2019_Bolgar; 2020_Chithrananda;
2022_White). Deep learning models, however, are known to struggle in small
data regimes to the extent that leading deep learning expert Yoshua Bengio
previously voiced a preference for a model called a Gaussian process (gp) for
small datasets (2011_Bengio). As such, leveraging them directly for small data
scientific discovery could prove to be difficult.

gps have received comparatively less attention relative to deep learning over
the past decade due to a variety of factors including a higher barrier to
entry in terms of the mathematical background required to use them, fewer
open-source software implementations, and perhaps most importantly, concerns
over the ability of gps to carry out representation learning, a stance summed
up in the following prescient quote from 2003_MacKay which foreshadows some of
the challenges currently encountered in supervised deep learning for the
sciences.

> "According to the hype of 1987, neural networks were meant to be intelligent
> models that discovered features and patterns in data. Gaussian processes in
> contrast are simply smoothing devices. How can Gaussian processes possibly
> replace neural networks? Were neural networks over-hyped, or have we
> underestimated the power of smoothing methods? I think both these
> propositions are true. The success of Gaussian processes shows that many
> real-world data modelling problems are perfectly well-solved by sensible
> smoothing methods. The most interesting problems, the task of feature
> discovery for example are not ones that Gaussian processes will solve. But
> maybe multilayer perceptrons can’t solve them either. Perhaps a fresh start
> is needed, approaching the problem of machine learning from a paradigm
> different from the supervised feedforward mapping."

One of the motivations for focussing on gps in this thesis, is the proposition
that many scientific discovery problems are instances of the real-world
problems described by MacKay. Furthermore, gps are more than just smoothing
devices. In addition to admitting exact Bayesian inference which can be used
to perform plausible reasoning (2003_Jaynes) over scientific hypotheses, gps
are also a longstanding workhorse of Bayesian optimisation (bo) and active
learning (2012_Settles), two methodologies that have already shown promise in
accelerating scientific discovery (2020_Pyzer; 2021_Shields). The goals of
this thesis are twofold: First, to showcase some of the use-cases for gps in
modelling scientific data and second, to extend current gp methodology and
software implementations to enable their application to scientific problems.
Specifically, the problem domains considered are:

  1. 1.

High-Energy Astrophysics \- It is challenging to test theories in high-energy
astrophysics due to the inability to perform physical experiments at the far
reaches of the universe. As such, the analysis of observational data is
important to guide the development of theory. It is shown how gp modelling can
play a role in performing inference over the structure of black hole accretion
discs and hence inform the development of future accretion disk theories.

  2. 2.

Photoswitch Chemistry \- In synthetic chemistry, new areas of chemical space
are constantly being explored and often little experimental data exists to
guide exploration. It is shown how gp modelling can be used for molecular
property prediction to prioritise the synthesis of novel molecules. We
validate the modelling approach with laboratory experiments, discovering new
and performant photoswitch molecules.

  3. 3.

Methodology/Software \- From a methodological standpoint a novel bo algorithm
is introduced that identifies and penalises input-dependent (heteroscedasatic)
measurement noise, an important consideration for the discovery of robust
materials suitable for industrial scale manufacturing. From a software
standpoint, an open-source gp library for chemistry is introduced, providing
implementations of bespoke kernels designed for common molecular and chemical
reaction representations.

Fig. 1.1: A pictorial overview of the thesis.

###  1.2 Overview and Contributions

A pictorial overview of the chapters of this thesis is available in Figure
1.1. The detailed summary and contributions of each chapter are as follows:

###### Chapter 2

The requisite background is provided on gps and bo, the machine learning
methodologies used across chapters of this thesis.

###### Chapter 3

A self-contained background is provided on the elements of high-energy
astrophysics required to understand our findings. The gapped lightcurves of
the Seyfert galaxy Markarian 335 (Mrk 335) are interpolated using gp modelling
with the intention of inferring the structure of the black hole accretion disk
through cross-correlation analysis. In a simulation study, Bayesian model
selection through the marginal likelihood is investigated as a means of
evaluating the most appropriate choice of gp kernel. Following gp modelling of
the observational data, it is found that the distance between the UV and X-ray
emission regions of Mrk 335 predicted by the Shakura-Sunyaev accretion disk
model is shorter than the light travel time measured using gp-based inference.
Tentative evidence is obtained for a short lag feature in the coherence and
lag spectra which could indicate the presence of an extended UV emission
region on the accretion disk where reverberation happens.

###### Chapter 4

A self-contained background is provided on the elements of molecular machine
learning required to understand the findings presented. GAUCHE is introduced,
a software library for Gaussian processes in chemistry, tackling the problem
of extending the gp framework to molecular representations such as graphs,
strings and bit vectors. By designing bespoke molecular kernels, the door is
opened to uncertainty quantification and bo directly on molecules and chemical
reactions.

###### Chapter 5

A small dataset of experimentally-determined properties for $405$ photoswitch
molecules is used in conjunction with the machinery made available in GAUCHE
to train a multioutput gp with a Tanimoto kernel to screen a large virtual
library of $7,265$ photoswitches, identifying 11 performant candidates
validated through laboratory experiment. Additionally, a predictive
performance comparison is conducted between the multioutput gp model and a
cohort of trained human photoswitch chemists with the gp model outperforming
the human experts. From a benchmark comparison against other machine learning
models, it is concluded that the curated dataset, as opposed to the choice of
model, is the key determinant of performance.

###### Chapter 6

A novel method for performing bo is introduced that is robust to experimental
measurement noise featuring a heteroscedastic gp surrogate model. From an
extensive empirical study, it is concluded that a moderately-sized
initialisation set is required for the model to be able to distinguish
heteroscedastic noise from intrinsic function variability. The chapter
concludes with recommendations on how future research might enable the
approach to be scaled to high-dimensional datasets.

###### Chapter 7

The thesis contributions are reviewed and discussed in the broader context of
identifying and enabling further applications of gps in the natural sciences.

###  1.3 List of Publications

What follows is the full list of publications co-authored during the PhD
process starting in October 2018. J1 (2021_Griffiths), J2 (2021_Mrk), J3
(2020_Griffiths) and W7 (2022_Gauche) comprise the thesis. In J1 and J2, all
coauthors acted in advisory roles, fine-tuning ideas and the final
manuscripts. I conducted all experiments, mathematical derivations and
implemented all code contributions. In J3, Aditya Raymond Thawani curated the
training dataset, designed and recruited participants for the human
performance comparison study and specified the set of performance criteria.
Jake Greenfield performed the spectral characterisation of the discovered
molecules in the Fuchter group laboratory at Imperial College London. I
conducted all machine learning experiments and implemented all code
contributions with the exception of the results in Table C.1 and Table C.2 of
Appendix C.3.1, where Penelope Jones, William McCorkindale, Arian Jamasb and
Henry Moss obtained results for the attentive neural process (ANP), smooth
overlap of atomic positions (SOAP) kernel, graph neural network (GNN) and
string kernel models respectively.

In W7, I ran all experiments excluding the Buchwald-Hartwig reaction
optimisation experiments which were run by Bojana Rankovic and the Weisfehler-
Lehman (WL) graph kernel table entries which were run by Aditya Ravuri. The
remaining co-authored articles are not included in the thesis for ease of
exposition.

J4 (2020_Griffiths) is a paper resulting from the continuation of my work from
the MPhil in Machine learning at the University of Cambridge in 2017. J5
(2020_Cheng) was work principally led by Dr. Bingqing Cheng. J6 (2020_Rivers)
and J7 (2020_Grosnit) were articles written during an internship at the Huawei
Noah’s Ark Lab. J8 (2020_Zagar) is a continuation of my work during an MSci at
Imperial College London in 2016. J9 (2022_Bourached) was principally led by
Anthony Bourached. C1 (2019_Grant) is a continuation of work undertaken whilst
a machine learning researcher at Secondmind Labs prior to commencement of the
PhD. C2-C5 (2022_Kell; 2021_Stork; 2021_Cann; 2021_Bourached_art) are articles
published in a domain unrelated to the topic of the thesis. The (unpublished)
workshop contributions, W1-W3 (2020_flowmo) are early versions of J3 and W7.
W4 (2018_Griffiths) is unrelated to the topic of the thesis although a figure
from this paper is used as Figure 4.2. W5 (2021_Aziz) was principally led by
Ajmal Aziz and so is not included in the thesis. W6 is a condensed version of
J8. P1 (2021_Grosnit) resulted from work undertaken whilst at Huawei Noah’s
Ark Lab and so is not included in the thesis though the subject matter is
related. P2 (bourached2021hierarchical) was principally led by Anthony
Bourached and so is not included in the thesis. P3 (2023_Frieder) and P4
(2022_Rankovic) was work led by Simon Frieder and Bojana Rankovic respectively
and is not included in the thesis.

Refereed Journal Papers

  1. [J1]

Griffiths RR, Aldrick A, Garcia-Ortegon M, Lalchand V, Lee, AA. Achieving
Robustness to Aleatoric Uncertainty with Heteroscedastic Bayesian
Optimisation. Machine Learning: Science and Technology. 2021.

  2. [J2]

Griffiths RR, Jiang J, Buisson D, Wilkins D, Gallo L, Ingram, A, Lee AA, Grupe
D, Kara M, Parker ML, Alston W, Bourached A, Cann G, Young A, Komossa S.
Modelling the Multiwavelength Variability of Mrk-335 using Gaussian Processes.
The Astrophysical Journal. 2021.

  3. [J3]

Griffiths RR, Greenfield JL, Thawani AR, Jamasb A, Moss HB, Bourached A, Jones
P, McCorkindale W, Aldrick AA, Fuchter, MJ, Lee AA. Data-Driven Discovery of
Molecular Photoswitches with Multioutput Gaussian Processes. Chemical Science.
2022.

  4. [J4]

Griffiths RR, Hernández-Lobato JM. Constrained Bayesian Optimization for
Automatic Chemical Design using Variational Autoencoders. Chemical Science.
2020.

  5. [J5]

Cheng B, Griffiths RR, Wengert S, Kunkel C, Stenczel T, Zhu B, Deringer VL,
Bernstein N, Margraf JT, Reuter K, Csanyi G. Mapping Datasets of Molecules and
Materials. Accounts of Chemical Research. 2020.

  6. [J6]

Cowen-Rivers A, Lyu W, Tutunov R, Wang Z, Grosnit A, Griffiths RR, Hao J, Wang
J, Bou-Ammar H. HEBO: Pushing the Limits of Sample-Efficient Hyper-parameter
Optimisation. Journal of Artificial Intelligence Research, 2022.

  7. [J7]

Grosnit A, Cowen-Rivers A, Tutunov R, Griffiths RR, Wang J, Bou-Ammar H. Are
We Forgetting About Compositional Optimisers in Bayesian Optimisation. Journal
of Machine Learning Research. 2021.

  8. [J8]

Zagar C, Griffiths RR, Podgornik R, Kornyshev AA. On the Voltage-Controlled
Self-Assembly of NP Arrays at Electrochemical Solid/Liquid Interfaces. Journal
of Electroanalytical Chemistry. 2020.

  9. [J9]

Bourached A, Griffiths RR, Gray R, Jha A, Nachev P. Generative Model-Enhanced
Human Motion Prediction. Applied AI Letters. 2021.

Refereed Conference Papers

  1. [C1]

Grant J, Boukouvalas A, Griffiths RR, Leslie D, Vaikili S, Munoz de Cote E.
Adaptive Sensor Placement for Continuous Spaces. International Conference on
Machine Learning. 2019.

  2. [C2]

Kell G, Griffiths RR, Bourached A, Stork D. Extracting Associations and
Meanings of Objects Depicted in Artworks through Bi-Modal Deep Networks,
Electronic Imaging 2022.

  3. [C3]

Stork D, Bourached A, Cann G, Griffiths RR. Computational Identification of
Significant Actors in Paintings through Symbols and Attributes, Electronic
Imaging, 2021.

  4. [C4]

Cann G, Bourached A, Griffiths RR, Stork D. Resolution Enhancement in the
Recovery of Underdrawings Via Style Transfer by Generative Adversarial Deep
Neural Networks, Electronic Imaging, 2021.

  5. [C5]

Bourached A, Cann G, Griffiths RR, Stork D. Recovery of Underdrawings and
Ghost-Paintings via Style Transfer by Deep Convolutional Neural Networks: A
Digital Tool for Art Scholars, Electronic Imaging, 2021.

Refereed Workshop Papers

  1. [W1]

Griffiths RR*, Moss H*. Gaussian Process Molecular Machine Learning with
FlowMO. NeurIPS Workshop on Machine Learning for Molecules. 2020 (Contributed
Talk - top 5%, * joint first authorship).

  2. [W2]

Griffiths RR, Jones P, McCorkindale W, Aldrick AA, Jamasb A, Day B.
Benchmarking Scalable Active Learning Strategies on Molecules. ICLR Workshop
on Fundamental Science in the Era of AI. 2020.

  3. [W3]

Griffiths RR, Thawani AR, Elijosius R. Enhancing the Diversity of Molecular
Machine Learning Benchmarks: An Open-Source Dataset for Molecular
Photoswitches. ICLR Workshop on Fundamental Science in the Era of AI. 2020.

  4. [W4]

Griffiths RR, Schwaller P, Lee AA. Dataset Bias in the Natural Sciences: A
Case Study in Chemical Reaction Prediction and Synthesis Design. NeurIPS
Workshop on Critiquing and Correcting Trends in Machine Learning. 2018.

  5. [W5]

Aziz A, Kosasih EE, Griffiths RR, Brintrup A. Data Considerations in Graph
Representation Learning for Supply Chain Networks. ICML Workshop on Machine
Learning for Data: Automated Creation, Privacy, Bias. 2021

  6. [W6]

Bourached A, Griffiths RR, Gray R, Jha A, Nachev P. Generative Model-Enhanced
Human Motion Prediction. NeurIPS Workshop on Interpretable Inductive Biases
and Physically-Structured Learning. 2020.

  7. [W7]

Griffiths RR, Klarner L, Moss Henry B., Ravuri A, Rankovic B, Truong S, Du Y,
Jamasb A, Schwartz J, Tripp A, Kell G, Bourached A, Chan A, Moss J, Guo C, Lee
AA, Schwaller P, Tang J, GAUCHE: A Library for Gaussian Processes in
Chemistry. ICML Workshop on AI4Science. 2022.

Preprints

  1. [P1]

Griffiths RR*, Grosnit A*, Tutunov R*, Maraval AM*, Cowen-Rivers A, Yang L,
Lin Z, Lyu W, Chen Z, Wang J, Peters J, Bou-Ammar H. High-Dimensional Bayesian
Optimisation with Variational Autoencoders and Deep Metric Learning. arXiv.
2021. (* joint first authorship)

  2. [P2]

Bourached A, Gray R, Griffiths RR, Jha A, Nachev P. Hierarchical Graph-
Convolutional Variational Autoencoding for Generative Modelling of Human
Motion. arXiv. 2021.

  3. [P3]

Frieder S, Pinchetti, L, Griffiths RR, Salvatori, T, Lukasiewicz, T, Petersen,
PC, Chevalier, A and Berner, J, 2023. Mathematical capabilities of ChatGPT..
arXiv. 2023.

  4. [P4]

Ranković, B, Griffiths, RR, Moss, HB and Schwaller, P. Bayesian optimisation
for additive screening and yield improvements in chemical reactions–beyond
one-hot encodings. ChemRxiv, 2022.

PhD Thesis

  1. [T1]

Griffiths RR, Applications of Gaussian Processes at Extreme Lengthscales: From
Molecules to Black Holes. University of Cambridge. 2022.

###  1.4 List of Software

The following list details the open-source software contributed to over the
duration of the PhD process:

  1. [S1]

Constrained Bayesian optimisation for automatic chemical design: Ryan-Rhys
Griffiths (2018). Code to reproduce the experiments from 2020_Griffiths.

Available at: https://github.com/Ryan-Rhys/Constrained-Bayesian-Optimisation-
for-Automatic-Chemical-Design

  2. [S2]

Mapping materials and molecules: Bingqing Cheng, Ryan-Rhys Griffiths, Tamas
Stenczel, Bonan Zhu, Felix Faber (2020). A software library containing
automatic selection tools for materials and molecules (2020_Cheng).

Available at: https://github.com/BingqingCheng/ASAP

  3. [S3]

Achieving robustness to aleatoric uncertainty with heteroscedastic Bayesian
optimisation: Ryan-Rhys Griffiths (2019). Code to reproduce the experiments
from 2021_Griffiths.

Available at: https://github.com/Ryan-Rhys/Heteroscedastic-BO

  4. [S4]

The photoswitch dataset: Ryan-Rhys Griffiths, Aditya Raymond Thawani, Arian
Jamasb, William McCorkindale, Penelope Jones (2020). Code to reproduce the
experiments from Chapter 5.

Available at: https://github.com/Ryan-Rhys/The-Photoswitch-Dataset

  5. [S5]

Modelling the multiwavelength variability of Mrk-335: Ryan-Rhys Griffiths
(2021). Code to reproduce the experiments from 2021_Mrk.

Available at: https://github.com/Ryan-Rhys/Mrk_335

  6. [S6]

An empirical study of assumptions in Bayesian optimisation: Alexander I.
Cowen-Rivers, Wenlong Lyu, Rasul Tutunov, Zhi Wang, Antoine Grosnit, Ryan-Rhys
Griffiths, Alexandre Max Maraval, Hao Jianye, Jun Wang, Jan Peters, Haitham
Bou-Ammar (2021). Code to reproduce the experiments from 2020_Rivers.

Available at: https://github.com/huawei-noah/HEBO/tree/master/HEBO

  7. [S7]

High-dimensional Bayesian optimisation with variational autoencoders and deep
metric learning: Antoine Grosnit, Rasul Tutunov, Alexandre Max Maraval, Ryan-
Rhys Griffiths, Alexander I. Cowen-Rivers, Lin Yang, Lin Zhu, Wenlong Lyu,
Zhitang Chen, Jun Wang, Jan Peters, Haitham Bou-Ammar. Code to reproduce the
experiments from 2021_Grosnit.

Available at: https://github.com/huawei-noah/HEBO/tree/master/T-LBO

  8. [S8]

Are we forgetting about compositional optimisers in Bayesian optimisation?:
Antoine Grosnit, Alexander I. Cowen-Rivers, Rasul Tutunov, Ryan-Rhys
Griffiths, Jun Wang, Haitham Bou-Ammar. Code to reproduce the experiments from
2020_Grosnit.

Available at: https://github.com/huawei-noah/HEBO/tree/master/T-LBO

  9. [S9]

FlowMO: Ryan-Rhys Griffiths and Henry Moss (2020). A GPflow library for
training Gaussian processes on molecular data (2020_Moss).

Available at: https://github.com/Ryan-Rhys/FlowMO

  10. [S10]

GAUCHE: Ryan-Rhys Griffiths, Leo Klarner, Henry Moss, Aditya Ravuri, Sang
Truong, Arian Jamasb, Austin Tripp, Bojana Rankovic, Philippe Schwaller
(2022). A software library for Gaussian processes in chemistry.

Available at https://github.com/leojklarner/gauche

  11. [S11]

Extracting associations and meanings of objects depicted in artworks through
bi-modal deep networks: Gregory Kell, Ryan-Rhys Griffiths (2021). Code to
reproduce the experiments from 2022_Kell.

Available at: https://github.com/gck25/fine_art_asssociations_meanings

##  Chapter 2 Background

In this chapter the requisite background is provided on Gaussian processes
(Chapters 3, 4, 5 and 6) and Bayesian optimisation (Chapters 4 and 6).

###  2.1 Gaussian Processes

In the context of machine learning, a Gaussian process (gp) is a Bayesian
nonparametric model for functions. gps are attractive models when limited data
is available, a setting common to many areas of the natural sciences, with
even notable deep learning experts voicing a preference for gps in the small
data regime (2011_Bengio). Furthermore, gps possess several important
properties for the applications in this thesis:

  1. 1.

Bayesian optimisation: gps have few hyperparameters that need to be determined
by hand which lends itself well to the repeated surrogate model hyperparameter
optimisation required by Bayesian optimisation.

  2. 2.

Astronomical time series: For astronomical time series, where noise processes
are often well understood, it is possible to incorporate this knowledge into
the design of the gp model.

  3. 3.

Molecules: gps maintain uncertainty estimates over molecular property values
through exact Bayesian inference. Uncertainty estimates are particularly
important when prioritising molecules for screening experiments.

A Gaussian process (gp) may be defined as a collection of random variables,
any finite subset of which have a joint Gaussian distribution
(2006_Rasmussen). In the cases considered in this thesis, the random variables
represent the value of the function $f(\mathbf{x})$ at location $\mathbf{x}$.
A stochastic process $f$ that follows a gp is written as

$f(\mathbf{x})\sim\mathcal{GP}\big{(}m(\mathbf{x}),k(\mathbf{x},\mathbf{x^{\prime}})\big{)}.$
(2.1)

The inputs to the gp may be scalars (e.g. time points in Chapter 3) or vectors
(e.g. molecular representations in Chapters 4 and 5). In the current
presentation we assume vector inputs $\mathbf{x}\in\mathbb{R}^{d}$ and we seek
to perform Bayesian inference over the latent function $f$ that represents the
mapping between the inputs $\\{\mathbf{x_{1}},\dotsc,\mathbf{x_{N}}\\}$ and
their function values $\\{f(\mathbf{x_{1}}),\dotsc,f(\mathbf{x_{N}})\\}$. The
gp is characterised by a mean function,

$m(\mathbf{x})=\mathbb{E}[f(\mathbf{x})],$ (2.2)

and a covariance function

$k(\mathbf{x},\mathbf{x^{\prime}})=\mathbb{E}[(f(\mathbf{x})-m(\mathbf{x}))(f(\mathbf{x^{\prime}})-m(\mathbf{x^{\prime}}))].$
(2.3)

In the absence of prior information on trends in the data, the mean function
is typically set to zero following standardisation of the outputs.
Standardisation, in this case refers to the common practice of subtracting the
mean and dividing by the standard deviation of the data when fitting the gp in
order to facilitate the identification of appropriate hyperparameters
(2008_Murray). The standardisation is reversed once the fitting procedure is
complete in order to obtain predictions on the original scale of the data.
$m(\mathbf{x})\equiv\mathbf{0}$ will be assumed henceforth for the sake of the
current presentation. The covariance function computes the pairwise covariance
between two random variables (function values). In the gp literature, the
covariance function is commonly referred to as the kernel. Informally, the
kernel is responsible for determining the properties of the functions which
the gp is capable of fitting e.g. smoothness and periodicity. The inductive
bias created by the choice of kernel is an important consideration in gp
modelling.

(a) SQE small lengthscale

(b) SQE large lengthscale

Fig. 2.1: gps with small and large lengthscales.

####  2.1.1 Kernels

The most widely-known kernel is the squared exponential (SQE) or radial basis
function (RBF) kernel,

$k_{\text{SQE}}(\mathbf{x},\mathbf{x^{\prime}})=\sigma_{f}^{2}\exp\bigg{(}\frac{-\left\lVert\mathbf{x}-\mathbf{x^{\prime}}\right\rVert^{2}}{2\ell^{2}}\bigg{)},$
(2.4)

where $\left\lVert\cdot\right\rVert$ is the Euclidean norm, $\sigma_{f}^{2}$
is the signal amplitude hyperparameter (vertical lengthscale) and $\ell$ is
the (horizontal) lengthscale hyperparameter. Although Equation 2.4 is written
with a single lengthscale shared across dimensions, for multidimensional input
spaces it is possible to optimise a lengthscale per dimension. We will adopt
the notation of $\theta$ to represent the set of kernel hyperparameters. An
illustration of gps with different lengthscales is given in Figure 2.1. It has
been argued by 2012_Stein that the smoothness assumptions of the SQE kernel
are unrealistic for many physical processes. As such, kernels such as the
Matérn,

$k_{\text{Mat\'{e}rn}}(\mathbf{x},\mathbf{x^{\prime}})=\frac{2^{1-\nu}}{\Gamma(\nu)}\bigg{(}\frac{\sqrt{2\nu}-\left\lVert\mathbf{x}-\mathbf{x^{\prime}}\right\rVert}{\ell}\bigg{)}^{\nu}K_{\nu}\bigg{(}\frac{\sqrt{2\nu}-\left\lVert\mathbf{x}-\mathbf{x^{\prime}}\right\rVert}{\ell}\bigg{)},$
(2.5)

are more commonly seen in the machine learning literature. Here $K_{\nu}$ is a
modified Bessel function of the second kind, $\Gamma$ is the gamma function
and $\nu$ is a non-negative hyperparameter of the kernel which is typically
taken to be either $\frac{3}{2}$ or $\frac{5}{2}$ (2006_Rasmussen). The
lengthscale hyperparameter $\ell$ can be thought of loosely as a decay
coefficient for the covariance between inputs as they become increasingly far
apart in the input space; the further apart the inputs are, the less
correlated they will be. The rational quadratic (RQ) kernel is defined as

$k_{\text{RQ}}(\mathbf{x},\mathbf{x^{\prime}})=\bigg{(}1+\frac{\left\lVert\mathbf{x}-\mathbf{x^{\prime}}\right\rVert^{2}}{2\alpha\ell^{2}}\bigg{)}^{-\alpha},$
(2.6)

where $\alpha,\ell>0$. The RQ kernel can be viewed as a scale mixture of SQE
kernels with different characteristic lengthscales. A comparison of the
functions drawn from gps with SQE and Matérn $\frac{5}{2}$ kernels is given in
Figure 2.2. The aforementioned kernels are defined over continuous input
spaces and are used in Chapters 3 and 6. For discrete input spaces such as
molecular representations it is necessary to define bespoke kernels which will
be introduced in Chapters 4 and 5.

(a) SQE

(b) Matérn $5/2$

Fig. 2.2: A comparison of the SQE and Matérn $\frac{5}{2}$ kernels.

####  2.1.2 Predictions

To obtain the predictive equations of gp regression, a mean function
$m(\mathbf{x})\equiv\mathbf{0}$ and kernel $k$ are specified and a gp prior
$p$ is placed over $f$,

$p(f(\mathbf{x})|\theta)=\mathcal{GP}\big{(}\mathbf{0},K_{\theta}(X,X^{\prime})\big{)}.$
(2.7)

The notation $K_{\theta}(X,X^{\prime})$ denotes a kernel matrix with entries
$[K]_{ij}=k(\mathbf{x}_{i},\mathbf{x}_{j})$ and the subscript notation is
chosen to indicate the dependence on the set of hyperparameters $\theta$ (e.g.
the signal variance $\sigma_{f}$ and lengthscale $\ell$ in Equation 2.4). We
suppress the explicit dependence on $\theta$ in the subsequent notation. It is
also necessary to specify a likelihood function

$p(y_{i}|f(\mathbf{x}_{i})),$ (2.8)

which depends on $f(\mathbf{x}_{i})$ only and is typically taken to be
Gaussian i.e.
$p(y_{i}|f(\mathbf{x}_{i}))=\mathcal{N}(y_{i}|f(\mathbf{x}_{i}),\sigma_{y}^{2})$.
The noise level $\sigma_{y}^{2}$ is most frequently assumed to be
homoscedastic, i.e. constant across the input domain. In Chapter 6,
heteroscedastic (input-dependent) noise is considered by introducing a
dependence $\sigma_{y}^{2}(\mathbf{x})$. The interpretation of
$\mathbf{y}_{i}$ is a noise-corrupted observation of the latent function
$f(\mathbf{x}_{i})$. Once data $\\{X,\mathbf{y}\\}$ has been observed, where
$X=\\{\mathbf{x}_{i}\\}_{i=1}^{N}$ and $\mathbf{y}=\\{y_{i}\\}_{i=1}^{N}$, the
joint prior distribution over the observations $\mathbf{y}$ and the predicted
function values $\mathbf{f}_{*}$ at test locations $X_{*}$ may be written

$\begin{bmatrix}\mathbf{y}\\\ \mathbf{f_{*}}\\\
\end{bmatrix}\sim\mathcal{N}\bigg{(}0,\begin{bmatrix}K(X,X)+\sigma_{y}^{2}I&K(X,X_{*})\>\\\
K(X_{*},X)\phantom{+\>\>\sigma{y}^{2}}&K(X_{*},X_{*})\end{bmatrix}\bigg{)},$
(2.9)

where $\mathcal{N}$ is the multivariate Gaussian probability density function
and $I\sigma_{y}^{2}$ represents the variance of iid Gaussian noise on the
observation vector $\mathbf{y}$. The joint prior in Equation 2.9 may be
conditioned on the observations through

$p(\mathbf{f_{*}}|\mathbf{y})=\frac{p(\mathbf{f_{*}},\mathbf{y})}{p(\mathbf{y})},$
(2.10)

which enforces that the joint prior agrees with the observations $\mathbf{y}$.
The posterior predictive distribution is then

$p(\mathbf{f_{*}}|X,\mathbf{y},X_{*})=\mathcal{N}\big{(}\mathbf{\bar{f}_{*}},\text{cov}(\mathbf{f_{*}})\big{)},$
(2.11)

with predictive mean at test locations $X_{*}$,

$\mathbf{\bar{f_{*}}}=K(X_{*},X)[K(X,X)+\sigma_{y}^{2}I]^{-1}\mathbf{y},$
(2.12)

and predictive uncertainty

$\text{cov}(\mathbf{f_{*}})=K(X_{*},X_{*})-K(X_{*},X)[K(X,X)+\sigma_{y}^{2}I]^{-1}K(X,X_{*}).$
(2.13)

Analysing the form of this expression one may notice that the first term
$K(X_{*},X_{*})$ in the expression for the predictive uncertainty
$\text{cov}(\mathbf{f}*)$ may be viewed as the prior uncertainty and the
second term $K(X_{*},X)[K(X,X)+\sigma_{y}^{2}I]^{-1}K(X,X_{*})$ can be thought
of as a subtractive factor that accounts for the reduction in uncertainty when
observing the data points $\mathbf{y}$. An illustration is given in Figure 2.3
of the posterior predictive distribution updates following data observation.

Fig. 2.3: An illustration of the gp posterior update on fitting 1, 2 and 4
data points (blue). The posterior distribution encodes the distribution over
possible functions that may explain the data.

####  2.1.3 Training

An important objective function for training gps is the log marginal
likelihood or evidence (1992_MacKay),

$\displaystyle\log p(\mathbf{y}|X,\theta)=$
$\displaystyle\underbrace{-\frac{1}{2}\mathbf{y}^{\top}(K_{\theta}(X,X)+\sigma_{y}^{2}I)^{-1}\mathbf{y}}_{\text{encourages
fit with data}}$ (2.14)
$\displaystyle\underbrace{-\frac{1}{2}\log|K_{\theta}(X,X)+\sigma_{y}^{2}I|}_{\text{controls
model capacity}}-\frac{N}{2}\log(2\pi).$

$N$ is the number of observations and $\theta$ again represents the set of
kernel hyperparameters to be optimised under the objective. The two terms in
the expression for the log marginal likelihood embody Occam’s Razor
(2001_Rasmussen) in their preference for selecting the simplest models that
explain the data well as illustrated in Figure 2.4. The first term in Equation
2.14 penalises functions that do not fit the data adequately whereas the
second term acts as a regulariser, disfavouring overly complex models. The
negative log marginal likelihood (NLML) is the gp training objective for all
experiments performed in this thesis.

Fig. 2.4: An illustration of the Bayesian Occam’s razor effect introduced by
2003_MacKay. If models are interpreted as probability distributions over
observations $y$, the x-axis may be viewed as the space of possible datasets.
The simple model (blue) explains datasets inside $S$ well, but not outside.
The complex model (red) explains datasets outside $S$ better, but worse inside
$S$. The probability density in $S$ must be lower so as to explain the
datasets outside $S$. Thus, if the dataset one wishes to model lies inside
$S$, Occam’s razor assigns preference to the simpler model.

####  2.1.4 Bayesian Model Selection

One desirable property of gps, and Bayesian models in general, is the ability
to carry out hierarchical modelling (1992_MacKay_Hierarchical). The three
tiers of the modelling hierarchy are:

  1. 1.

Model Parameters

  2. 2.

Model Hyperparameters

  3. 3.

Model Structures

In the case of the nonparametric gp framework, model parameters do not have
the same meaning as in parametric Bayesian models and are instead obtained
from the posterior distribution over functions. Model hyperparameters consist
of parameters of the kernel function such as signal amplitudes and
lengthscales as well as the likelihood noise. At the level of model
structures, the fit achieved by different kernels can be quantitatively
assessed by comparing the values of the optimised NLML objective permitting
Bayesian model selection, a procedure that is used in Chapter 3. The next
discussion point to be considered is an important application of gps in
mathematical optimisation.

###  2.2 Bayesian Optimisation

Bayesian optimisation (bo) (1962_Kushner; 1964_Kushner; 1975_Mockus;
1975_Zhilinskas; 1978_Mockus) is a data-efficient methodology for solving
black-box optimisation problems.

####  2.2.1 Black-Box Optimisation

In many problems in science and engineering we are interested in solving
global optimisation problems of the form

$\mathbf{x}^{\star}=\arg\max_{\mathbf{x}\in\mathcal{X}}f(\mathbf{x}),$ (2.15)

where $f:\mathcal{X}\rightarrow\mathbb{R}$ is a function over an input domain
$\mathcal{X}$ which is typically a compact subset of $\mathbb{R}^{d}$ (Chapter
6) but may also be non-numeric in the case of molecular representations such
as graphs and strings (Chapter 4). Equation 2.15 is also a black-box
optimisation problem in the sense that it possesses the following properties:

  1. 1.

Black-Box Objective: We do not have the analytic form of $f$ nor do we have
access to its gradients. We can, however, evaluate $f$ pointwise anywhere in
the input domain $\mathcal{X}$.

  2. 2.

Expensive Evaluations: Choosing an input $\mathbf{x}$ and evaluating
$f(\mathbf{x})$ takes a very long time or incurs a large financial cost.

  3. 3.

Noise: The evaluation of a given $\mathbf{x}$ is a noisy process. In addition,
this noise may vary across $\mathcal{X}$, making the underlying process
heteroscedastic.

A motivating example is molecular property optimisation where the input domain
$\mathcal{X}$ is a set of molecular graphs
$\\{\mathbf{x}:\mathbf{x}\in\mathcal{X}\\}$ and the black-box function
$f(\mathbf{x})$ is the property of the molecule to be optimised. $f$ maps a
molecule to its property, but its analytic form is unknown and so instead $f$
must be queried through experiment by synthesising a molecule and measuring
the value of its property under $f$. This is a time-consuming and financially
expensive process. In addition, the measurement process using laboratory
equipment is typically noisy.

####  2.2.2 Solution Methods

In the absence of an analytic form for the function to be optimised,
strategies for solving black-box optimisation problems tend to proceed by
sequentially evaluating the black-box function until the global optimum is
found or the evaluation budget is exhausted. Such strategies may be
represented by the abstract blueprint of sequential optimisation outlined in
Algorithm 1.

Algorithm 1 Sequential Optimisation

input: initial dataset $\mathcal{D}$ $\triangleright$ may be empty

repeat

$\mathbf{x}\leftarrow\text{Policy}(\mathcal{D})$ $\triangleright$ select the
next input

$y\leftarrow\text{Evaluate}(\mathbf{x})$ $\triangleright$ evaluate the black-
box at the chosen input

$\triangleright$ update the dataset

until termination condition reached $\triangleright$ e.g. evaluation budget
exhausted

return $\mathcal{D}$

Sequential optimisation algorithms differ in their choice of policy, or in
other words, how they make use of the dataset of evaluations $\mathcal{D}$.
Strategies may be non-adaptive in the sense that they ignore $\mathcal{D}$
completely, or they may be adaptive in the sense that they use the information
about the black-box function stored within $\mathcal{D}$ to inform the
selection of the next input $\mathbf{x}$ (2022_Garnett). Some of the most
relevant solution methods for black-box optimisation include:

Grid Search: Perhaps the most well-known strategy for black-box optimisation
problems, such as machine learning hyperparameter tuning, is grid search. Grid
search is a deterministic, non-adaptive strategy where the policy consists of
an exhaustive search through the input domain $\mathcal{X}$ by manually
specifying a subset of inputs to query. Typically the manually-specified
inputs are evenly spaced throughout the input domain and hence assume the form
of a grid. Grid search suffers from the curse of dimensionality (1957_Bellman)
since the number of inputs to evaluate grows exponentially as a function of
the dimensionality of x. Grid search is still a popular strategy in practice,
however, due to its ease of implementation and the fact that it is
embarrassingly parallel in so far as evaluations tend to be independent of
each other.

Random Search: This stochastic, non-adaptive strategy consists of draws from a
uniform density over the input domain $\mathcal{X}$. It has been demonstrated
empirically that in high dimensions, random search can often outperform grid
search due to its robustness to non-informative dimensions of the input space
(2012_Bergstra). Random search is used as a baseline strategy in Chapters 4
and 6.

Bayesian Optimisation: A third solution method is an adaptive strategy where
the policy is derived from Bayesian decision theory (2005_De_Groot;
1985_Berger; 2007_Robert) and formalises the approach to decision-making under
uncertainty with respect to the unknown objective function. bo, which is the
principal subject of Chapter 6 and plays a major role in Chapter 4, has
recently achieved notable and widely-publicised success as a component of
AlphaGo (2018_Yutian) as well as across applications including chemical
reaction optimisation (2021_Shields), robotics (2016_Calandra), and machine
learning hyperparameter optimisation (2021_Turner; 2020_Rivers). bo will be
the focus from hereon in.

Fig. 2.5: An illustration of the Bayesian optimisation algorithm.

####  2.2.3 The Bayesian Optimisation Algorithm

The bo algorithm, illustrated in Algorithm 2, implements the policy from
Algorithm 1 through the use of two components:

Algorithm 2 Bayesian Optimisation

input: initial dataset $\mathcal{D}$ $\triangleright$ may be empty

repeat

choose $\mathbf{x}$ by optimising $\alpha$, the acquisition function

$\mathbf{x}\leftarrow\operatorname*{arg\,max}_{\mathbf{x}\in\mathcal{X}}\alpha(\mathbf{x};\mathcal{D})$

$y\leftarrow\text{Evaluate}(\mathbf{x})$ $\triangleright$ evaluate the black-
box at the chosen input

$\triangleright$ update the dataset and surrogate model

until termination condition reached $\triangleright$ e.g. evaluation budget
exhausted

return $\mathcal{D}$

Surrogate Model: A flexible probabilistic model that captures the prior belief
about the behaviour of the black-box objective $f(\mathbf{x})$. A
probabilistic model is necessary to ensure that uncertainty in the values of
the black-box objective is maintained across the design space. This
uncertainty measure is then used to inform the data collection policy known as
the acquisition function. When a new data point is collected, the surrogate
model is updated by means of re-training.

Acquisition Function: The acquisition function
$\alpha(\mathbf{x},\mathcal{D})$ determines the next input on a given
iteration of bo by leveraging the uncertainty estimates of the surrogate model
to trade off exploration and exploitation. It is beneficial to explore regions
of the design space where the value of the objective is unknown, yet with a
finite budget of function evaluations it is desirable to exploit the knowledge
acquired to locate an input close to the global optimum of the function. From
a computational standpoint, the acquisition function should be cheaper to
evaluate relative to the black-box function. It should also be easy to
optimise (2018_Wilson; 2020_Grosnit; 2020_Schweidtmann).

The pseudocode for bo in Algorithm 2 does not represent a single instantiation
of an algorithm but rather a class of algorithms reflecting the broad range of
choices available for both the surrogate model and the acquisition function.
The set of criteria for choosing the surrogate model and the acquisition
function will now be discussed.

####  2.2.4 The Surrogate Model

The desiderata for the surrogate model in bo are often related to the quality
of the posterior distribution and the scalability of the model. In an
idealised scenario, viewing Bayesian inference as an optimal calculus for
dealing with incomplete information (2010_Turner; 1961_Cox; 2003_Jaynes;
2003_MacKay), one would obtain uncertainty estimates using full Bayesian
inference over the surrogate model posterior. Full Bayesian inference is
computationally demanding however and can be infeasible if the bo problem
features a large dataset or a long horizon of function evaluations. To date,
gps have been the model of choice for bo on small datasets due to the ability
to perform full Bayesian inference. gp surrogates have the following strengths
and weaknesses from the point of view of bo:

###### Strengths

  1. 1.

Full Bayesian inference, admits a closed-form posterior predictive
distribution via exact inference. In contrast, approximate inference methods
run the risk of degrading the quality of the uncertainty estimates
(2020_Foong). The importance of uncertainty estimate quality in obtaining
strong empirical performance is regularly emphasised in the bo literature
(2015_Shahriari; 2022_Garnett).

  2. 2.

We can perform Bayesian model selection at the hyperparameter level meaning
that we are more robust to overfitting. This is facilitated by an analytic
form for the marginal likelihood.

  3. 3.

Few of the gps hyperparameters needs to be determined by hand for example
through hyperparameter search routines. This makes gps well-suited to problems
such as bo in which running hyperparameter search per iteration of the bo loop
is not practically feasible (2003_MacKay).

###### Weaknesses

  1. 1.

Common choices of gp kernels are stationary kernels, meaning they cannot
accurately model situations in which the complexity of the objective function
varies in different regions of the input space. While non-stationary kernels,
warping functions (2020_Rivers; 2019_Balandat), deep gps (2013_Damianou;
2021_Hebbal), and normalising flows (2020_Maronas) are potential solutions,
they introduce additional complexity into the bo algorithm.

  2. 2.

The gp marginal distribution is not heavy-tailed. If outlier detection is a
concern for example, one may wish to employ a heavier-tailed distribution such
as the student T-process of 2014_Shah which has shown some success as a
surrogate for bo (2017_Cantin).

  3. 3.

The observation model assumes homoscedastic Gaussian noise. While
modifications to the standard gp framework exist to capture more complex noise
distributions (2021_Griffiths; 2021_Makarova), they likely require more data
in order to operate effectively.

  4. 4.

The most frequently cited downside of the gp framework for bo is the
computational complexity of performing full Bayesian inference. Computing the
inverse of the covariance matrix
$[K(\mathbf{X},\mathbf{X})+I\sigma^{2}_{y}]^{-1}$ is $O(N^{3})$ in the number
of data points $N$. This covariance matrix appears in the expression for the
marginal likelihood in addition to the predictive mean and covariance. A
mitigating factor is that, for a fixed set of kernel hyperparameters, the
Cholesky decomposition of this matrix may be computed once and stored,
yielding a complexity of $O(N^{2})$ for future predictions. In bo, however,
the kernel hyperparameters are recomputed each time a new data point is
collected. The $O(N^{3})$ complexity cannot be avoided in this instance.
Scalable surrogate model alternatives such as deep neural networks (DNNs)
(2015_Snoek; 2016_Springenberg; 2018_Perrone; 2021_White), sparse gps
(2018_Design; 2020_Griffiths), and transformers (2022_Maraval) have been
trialled but face challenges in terms of the quality of the model uncertainty
estimates.

  5. 5.

gps often struggle to model functions in high-dimensional, continuous input
spaces. In as little as $10$ input dimensions, the predictive capabilities of
gps can be impaired because the covariance function stipulates that inputs
separated by more than a few lengthscales are negligibly correlated
(2014_Garnett). As such, the majority of the input domain $\mathcal{X}$ may be
uncorrelated with the observed data making prediction challenging. Some
popular approaches in high-dimensional spaces include embedding methods such
as variational autoencoders (VAEs) which seek to learn a low-dimensional
embedding of the input data (2018_Design; 2020_Griffiths; 2021_Grosnit;
2021_Verma; 2022_Hie; 2022_Maus).

In the bo problems considered in this thesis, however, many of the
aforementioned limitations of gps do not apply. The scientific datasets lie in
the small data regime due to factors such as the expense of collecting
laboratory measurements of synthesised molecules or the limited observational
history of celestial objects and so the scalability of the surrogate model is
not an issue. Similarly, the only high-dimensional input space considered is
that of molecular fingerprints in which each input dimension is binary and so
the problem of extrapolation in high-dimensional, continuous input spaces is
avoided. The only exception is the case of the attempt to model
heteroscedastic noise distributions in Chapter 6. In this case a bespoke
heteroscedastic gp surrogate and acquisition function is devised.

####  2.2.5 The Acquisition Function

A sequential optimisation algorithm such as that defined in Algorithm 1
requires a policy or acquisition function $\alpha:\mathcal{X}\to\mathbb{R}$ to
provide a score for each potential observation.

###### Evaluation of Policies:

The ideal performance metric for a bo scheme would quantify how close the set
of queried inputs were to the global optimum of the black-box function. Regret
is one such metric for quantifying optimisation performance that is commonly
used in the analysis of optimisation algorithms. While there are many
formulations of regret in different contexts, the central idea is to compare
the values of the objective function visited during optimisation with the
value of the global optimum. The larger the gap between these values is, the
more regret is retrospectively incurred. The instantaneous regret is defined
as $r_{n}(\mathbf{x})=f^{*}-f(\mathbf{x}_{n})$, where $f^{*}$ is the global
optimum of the black-box and $f(\mathbf{x}_{n})$ is the value of the function
queried at the input $\mathbf{x}_{n}$ at iteration $n$. Two derived forms of
regret are typically used in theoretical analysis of optimisation algorithm
performance:

  1. 1.

Simple Regret \- gives the instantaneous regret at the final iteration of bo
as $r_{\tau}=f^{*}-\max f(\mathbf{x}_{\tau})$, where $\mathbf{x}_{\tau}$
represent the set of inputs queried at the terminal iteration $\tau$. This
metric has the advantage of not punishing the algorithm for explorative
queries early in the search procedure.

  2. 2.

Cumulative Regret \- is defined as $R_{N}=\frac{1}{N}\sum_{n=0}^{N}r_{n}$,
where $r_{n}$ is the instantaneous regret at iteration $n$. Thus the
cumulative regret is an average over all queries. The simple regret can be
obtained by taking the last term, $r_{N}$, in the expression for the
cumulative regret as the performance metric and setting $N=\tau$, where $\tau$
is the terminal iteration.

###### Designing an Optimal Policy:

In terms of designing an optimal policy, however, the regret metric cannot be
used directly since the global optimum $f^{*}$ is unknown. In this instance a
concept from Bayesian decision theory known as the utility
$u(a,\psi,\mathcal{D})$ can be applied, where $a$ represents the action i.e. a
choice of query location $\mathbf{x}$, $\psi$ represents the uncertain
elements in the optimisation problem e.g. the objective function values, and
$\mathcal{D}$ is the dataset of input/observation pairs collected so far.
Maximising the expected utility of the data returned by the optimisation
algorithm at a given iteration, however, requires consideration of the entire
remainder of the optimisation query budget. Under such long time horizons the
optimal policy becomes prohibitive to compute. Some attempts have been made to
approximate an idealised look-ahead policy (2009_Garnett; 2010_Ginsbourger;
2016_Osborne) but in practice most bo policies take the form of acquisition
functions; myopic heuristics that attempt to trade off exploration and
exploitation. When linked with the probabilistic surrogate model, this
translates to greedily selecting queries which have high (for maximisation
problems) predictive mean (exploitation) and high predictive variance
(exploration).

###### Classes of Acquisition Functions

While there exist a broad range of acquisition functions (2015_Shahriari;
2020_Grosnit), a large subset of commonly-used acquisitions can be divided
into three classes:

  1. 1.

Optimistic Acquisition Functions \- An example of this type of acquisition
function is the Upper Confidence Bound (UCB) (2010_Srinivas). In the bandits
literature these methods are described by the term optimism in the face of
uncertainty because they assign higher values to actions with high
uncertainty.

  2. 2.

Improvement-Based Acquisition Functions \- Are defined relative to some
incumbent target, typically taken to be the best queried value found so far in
the optimisation. Examples of this class include Probability of Improvement
(PI) (1964_Kushner) and Expected Improvement (EI) (1971_Saltenis; 1978_Mockus;
1998_Jones). The EI acquisition will be used and extended in Chapter 6.

  3. 3.

Information-Based Acquisition Functions \- In these methods, the posterior
over the unknown optimiser $\mathbf{x}_{*}$, induced implicitly by the
posterior distribution over objective functions, is used as a means of
selecting queries. Instances of this class of acquisition function include
Thompson Sampling (TS) (1933_Thompson), Entropy Search (ES) (2012_Hennig),
Predictive Entropy Search (PES) (2014_Lobato), General-Purpose Information-
Based Bayesian Optimisation (GIBBON) (moss2021gibbon), and the Informational
Approach to Global Optimization (IAGO) (2009_Villemonteix).

Additionally, ensembles of acquisition functions known as portfolios are
popular in practice and may perform better than any individual acquisition
function (2011_Hoffman; 2013_Hoffman; 2014_Shahriari; 2020_Rivers).

##  Chapter 3 Modelling Black Hole Signals with Gaussian Processes

Status: Published as Griffiths, RR., Jiang, J., Buisson, DJ., Wilkins, D.,
Gallo, LC., Ingram, A., Grupe, D., Kara, E., Parker, ML., Alston, W.,
Bourached, A. Cann, G., Young, A., Komossa, S., Modeling the Multiwavelength
Variability of Mrk 335 Using Gaussian Processes. The Astrophysical Journal,
2021.

###  3.1 Background on High-Energy Astrophysics

The chapter begins with a self-contained background on high-energy
astrophysics to aid in contextualising the findings.

####  3.1.1 Black Holes

John Michell was the first to posit the existence of black holes
(1784_Michell), describing them as dark stars due to the fact that no light
could escape from them. At the time, however, his work was largely ignored due
to the absence of a theory of gravity describing the behaviour of light in a
strong gravitational field. Following the introduction of the Einstein field
equations from General Relativity (1916_Einstein), Schwarzchild was the first
to calculate the radius of a black hole in the Schwarzschild metric
(1916_Schwarzschild). The Schwarzchild radius is

$r_{s}=\frac{2GM}{c^{2}},$ (3.1)

where $G$ is the gravitational constant, $M$ is the mass of the object and $c$
is the speed of light. Black holes are characterised according to their mass
and spin. When the mass of a black hole exceeds $10^{5}M_{\odot}$, it is
termed a supermassive black hole (SMBH), where $M_{\odot}$ is the solar mass
unit, approximately equal to the mass of the Sun.

####  3.1.2 Active Galactic Nuclei

The term Active Galactic Nucleus (AGN) was coined by Viktor Ambartsumian in
the early 1950s (1997_Victor). Ambartsumian argued that the nuclei of galaxies
were subject to explosions which caused large amounts of mass to be expelled,
and that for these explosions to occur, galactic nuclei must contain unknown
bodies of huge mass. Moreover, AGN were observed to be highly luminous with
unusual spectral properties, indicating that their power source could not be
ordinary stars. In 1964, some insight on the nature of AGN was offered by
Salpeter and Zeldovich (1964_Salpeter; 1964_Zeldovich), who proposed accretion
of gas onto a SMBH as the mechanism responsible for the power source of a
powerful class of AGN known as quasars. 1969_Lynden later paid testament to
the importance of the black hole accretion disc model by remarking that,

> With different values of the black hole mass and accretion rate these discs
> are capable of providing an explanation for a large fraction of the
> incredible phenomena of high-energy astrophysics.

Lynden-Bell’s statement is supported by the fact that AGN are one of the most
persistent luminous sources of electromagnetic radiation in the universe and
as such, may be leveraged to discover distant objects. Furthermore, the
evolution of AGN in cosmic time may be used to inform theoretical models of
the cosmos. It is estimated that one fifth of research astronomers work on
AGNs (1997_Peterson).

The observed properties of AGNs depend on the mass of the central SMBH, the
extent that the nucleus is obscured by dust, the orientation of the accretion
disc, the rate of gas accretion, as well as the presence or absence of
outflows of ionised matter along the axis of rotation known as jets. Some
subclasses of AGN include quasars, the most powerful form of AGN, blazars,
which contain a jet pointed toward the Earth, and Seyfert galaxies which are
characterised by broad emission lines in the optical band. It is the last of
these categories of AGN that is the subject of this chapter and will be
discussed next.

####  3.1.3 Seyfert Galaxies

In 1943, Carl Seyfert systematically studied a collection of bright AGN
possessing broad emission lines in the optical band (1943_Seyfert). The
eponymous Seyfert galaxies are further subdivided into Seyfert 1 (Sy1) and
Seyfert 2 (Sy2) galaxies based on their emission line range,
$1000-20,000\text{kms}^{-1}$ and $300-1000\text{kms}^{-1}$ respectively. The
orientation-based unified model of is one of the most popular means of
describing Seyfert galaxies and is based on the idea that classes of AGN are
physically similar but are viewed at different orientations (1993_Antonucci;
1995_Urry). Some features of the model include:

  * •

In the narrow line region there is ionised, low-velocity and low-density gas
extending to $100-1000$ parsecs (pc).

  * •

In the broad line region there are high-density, dust-free gas clouds located
at a distance of $0.01-1$ pc from the SMBH moving at Keplerian velocities.

  * •

There is an antisymmetric dusty structure known as a torus located at a
distance of $0.1=10$ pc from the SMBH.

  * •

There is a sub-pc accretion disc located around the SMBH which may be
optically thick or optically thin depending on the disc’s state.

  * •

There is an outflowing radio jet pointed in the general direction of the
accretion disc.

Fig. 3.1: The orientation-based unified model of AGN. Reprinted with
permission from 2019_Jiang.

A schematic for the orientation-based unified model is provided in Figure 3.1.
In Sy2 galaxies, the narrow line region is viewable by an edge-on observer due
to the fact that it is more extended relative to the broad line region. Both
the broad line region and the accretion disc are obscured by the middle plane
of the torus. In Sy1 galaxies, the observer is closer to the torus axis and
has an unobscured line of sight towards the nuclear region of the AGN. While
the orientation-based unified model can explain the spectral diversity of many
AGN, it has recently been challenged. For example, a different torus shape is
suggested by interferometry observations in the mid-infrared band of Sy1
galaxies which demonstrate that the majority of infrared emission originates
from dust in the polar region as opposed to the disc plane (2013_Honig).

####  3.1.4 Space Observatories

High energy X-rays function as the primary tool for examining the innermost
regions of AGN as they originate from the area closest to the central SMBH and
can easily penetrate absorbing materials along the line of sight. Optical/UV
emission is also useful in characterising the behaviour of AGN, for example,
in verifying reprocessing models of X-ray emission by computing lags between
X-ray and optical/UV lightcurves, where a lightcurve is a graph of the light
intensity of a celestial object or region with respect to time.

All astronomical data in this thesis originates from the Neil Gehrels Swift
observatory which was first launched by NASA in 2004 to detect and study
Gamma-Ray Bursts (GRBs) using the Burst Alert Telescope (BAT)
(2005_Barthelmy). While originally designed for the study of GRBs, Swift now
also functions as a multiwavelength observatory containing an X-ray Telescope
(XRT) (2005_Burrows) and a UV/Optical Telescope (UVOT) (2005_Roming). Swift is
also used to conduct long-term all sky surveys.

####  3.1.5 Accretion Disc Models

Theoretical models to explain accretion discs differ based on the physical
processes considered. Four representative examples are the Polish doughnut
(thick disc), Shakura-Sunyaev (thin disc), slim disc, and the advection-
dominated accretion flow (ADAF) models (2011_Abramowicz). These theoretical
models are not mutually exclusive in the sense that different aspects of real
physical systems may be best described by different models.

Polish Doughnut (Thick Disc) Model: The Polish Doughnut introduced by
Paczynski and collaborators in the 1970s and 80s (1980_Jaroszynski;
1980_Paczynski; 1981_Paczynski; 1982_Paczynski) is the minimal analytic
accretion disc model in so far as it only considers gravity and assumes a
perfect fluid. The Polish Doughnut Model is predicated on a general method for
constructing perfect fluid equilibria of matter orbiting an uncharged,
rotating black hole known as a Kerr black hole (1963_Kerr).

Thin Disc Models: The majority of analytic accretion disc models assume a
stationary and axially-symmetric state for matter undergoing accretion onto
the black hole with all physical quantities depending only on two spatial
coordinates, $r$, the radial distance from the disc centre, and $z$, the
vertical distance from the equatorial plane of symmetry. Unlike the Polish
Doughnut Model, which assumes vertically thick discs, in thin disc models,
$\frac{z}{r}\ll 1$ applies at all points within the matter distribution. In
1973, Shakura and Sunyaev introduced the canonical thin disc model
(1973_Shakura) by specifying additional physically reasonable assumptions that
allowed them to construct a set of algebraic equations from the standard set
of thin disc equations. The relativistic extension of the Shakura-Sunyaev
model was later proposed by 1973_Novikov.

Slim Disc Models: Slim discs are characterised by $\frac{z}{r}\leq 1$. Thin
disc models such as the Shakura-Sunyaev and Novikov–Thorne models assume that
viscous heating is balanced locally by radiative cooling i.e. the accretion
process is radially efficient and as such, all viscosity-generated heat is
radiated away. Although the assumption is valid if the accretion rate is
small, at a luminosity $L=0.3L_{\text{Edd}}$ 111$L_{\text{Edd}}$ is the
Eddington limit, the maximum achievable luminosity of a body subject to the
balance of an outward radiative force and an inward gravitational force. the
radial velocity is large and the disc is sufficiently thick to permit
advection to function as a cooling mechanism. At the highest luminosities,
thin disc models no longer apply as the cooling effect of advection becomes
comparable to radiative cooling. The standard thin disc model equations become
a two-dimensional system of ordinary differential equations with a critical
point for the slim disc case. These equations were first solved by
1988_Abramowicz and extended to a fully relativistic treatment by
1998_Beloborodov.

ADAF Models: Advection-dominated accretion flow models, introduced first in a
series of papers (1995_Narayan; 1994_Narayan; 1995_Abramowicz;
1996_Abramowicz; 1998_Gammie), assume that almost all viscously dissipated
energy is not radiated but advected into the black hole and applies when the
luminosity and mass accretion rate are low. As such, ADAF discs are typically
far less luminous than thin discs. Fully relativistic solutions to such discs
have been obtained numerically (1997_Abramowicz; 1997_Beloborodov). Further
information on ADAFs is available in 1997_Narayan.

####  3.1.6 Markarian 335

The accretion disc of Markarian 335 (Mrk 335) is the focus of study in this
chapter. Mrk 335 is a Sy1 galaxy located 324 million light-years from Earth in
the constellation of Pegasus. The central SMBH of Mrk 335 is notable for the
spinning rate of its corona at ca. 20% the speed of light. Relativistic
blurring of the reflection of the accretion disc has been used to infer the
geometry of the corona (2015_Wilkins_drive). By using gps to interpolate the
unevenly-sampled lightcurves of Mrk 335 and performing a cross-correlation
analysis, some insight into the structure of the accretion disc may be
obtained, and subsequently used to inform future developments in accretion
disc theories. Of the aforementioned disc theories, the Shakura-Sunyaev thin
disc model is the most relevant for Mrk 335 as its predictions for the extent
of UV emission match that from observation. The distance between the UV and
X-ray emission regions however is shorter than the light travel time measured
using gp-based inference on the observational data. The main contributions of
this chapter are now introduced.

###  3.2 Preface

The optical and UV variability of the majority of AGN may be related to the
reprocessing of rapidly-changing X-ray emission from a more compact region
near the central black hole. Such a reprocessing model would be characterised
by lags between X-ray and optical/UV emission due to differences in light
travel time. Observationally, however, such lag features have been difficult
to detect due to gaps in the lightcurves introduced through factors such as
source visibility or limited telescope time. In this chapter, gp regression is
employed to interpolate the gaps in the Swift X-ray and UV lightcurves of the
narrow-line Seyfert 1 galaxy Mrk 335. In a simulation study of five commonly-
employed analytic gp kernels, it is concluded that the Matérn $\frac{1}{2}$
and rational quadratic kernels yield the most well-specified models for the
X-ray and UVW2 bands of Mrk 335. In analysing the structure functions of the
gp lightcurves, a broken power law is obtained with a break point at 125 days
in the UVW2 band. In the X-ray band, the structure function of the gp
lightcurve is consistent with a power law in the case of the RQ kernel, whilst
a broken power law with a break point at 66 days is obtained from the Matérn
$\frac{1}{2}$ kernel. The subsequent cross-correlation analysis is consistent
with previous studies and furthermore, shows tentative evidence for a broad
X-ray-UV lag feature of up to 30 days in the lag-frequency spectrum. The
significance of the lag depends on the choice of gp kernel.

###  3.3 Introduction

AGN show strong and variable emission across multiple wavelengths. The UV
emission from an AGN is believed to be dominated by thermal emission from an
accretion disc close to the central SMBH (pringle81). The variability of
optical and UV AGN 222AGN with an UV and optical luminosity change of more
than 1 magnitude such as changing-look AGN, are not discussed in this chapter
cf. jiang21 for details. emission is stochastic and described by random
Gaussian fluctuations (welsh11; gezari13; zhu16; sanchez18; smith18; xin20)
with the autocorrelation functions of such fluctuations adhering to the
‘damped random walk’ model. The X-ray emission from an AGN is often found to
show faster variability relative to emission at longer wavelengths
(mushotzky93; gaskell03) and originates from a more compact region (morgan08;
chartas17).

The relationship between the UV and X-ray emission has been well studied. For
instance, correlations between the variability in two energy bands has been
seen in some individual sources (shemmer01; buisson17) while others do not
show significant evidence for similar correlation (smith07; 2018_Buisson). In
sources where correlation is found, lags that are related to the light travel
time between two emission regions are frequently observed. These lags are
often found to be on timescales of days and are longer than those predicted by
classical disc theories (1973_Shakura). Such lag amplitudes indicate a disc of
size a few times larger than expected (edelson00; shappee14; troyer16;
buisson17). Alternatively, some modified models have been proposed for the
underestimation of lags by the classical thin disc model, e.g. disc turbulence
(cai20), additional varying FUV illumination (gardner17), a tilted or
inhomogeneous inner disc (dexter11; starkey17) or an extended coronal region
(kammoun20). Much shorter lags, e.g. hundreds of seconds, in agreement with
the Shakura-Sunyaev model (1973_Shakura) have been rarely observed by
comparison e.g. in NGC-4395 (mchardy16).

The Neil Gehrels Swift Observatory has been monitoring the X-ray sky in the
past decade in tandem with simultaneous pointings in the optical and UV band.
In this work, we focus on the X-ray and UVW2 ($\lambda=$212 nm) lightcurves of
the narrow-line Seyfert 1 galaxy (NLS1) Mrk 335 obtained by XRT and UVOT, the
soft X-ray and UV/optical telescopes on Swift. Mrk 335 was one of the
brightest X-ray sources prior to 2007, before its flux diminished by
$10-50\times$ its original brightness (grupe07). The X-ray brightness has not
recovered since. During this low X-ray flux period, the UV brightness remains
relatively unchanged rendering Mrk 335 X-ray weak (tripathi20). The behavior
has been explained as a possible collapse of the X-ray corona (2013_Gallo_New;
2015_Gallo_New; 2014_Parker_New) and/or increased absorption in the X-ray
emitting region (grupe12; 2013_Longinotti; 2019_Longinotti; 2019_Parker).

Mrk 335 has been continuously monitored since 2007 making it one of the best-
studied AGN with Swift. Previous studies from the Swift monitoring program can
be found in grupe07; grupe12; gallo18; tripathi20; 2020_Komossa. The X-rays
are constantly fluctuating and regularly display large amplitude flaring
(2015_Wilkins). The UV are significantly variable, but at a much smaller
amplitude than the X-rays. gallo18 found tentative evidence for lags of
$\approx 20$ days based on cross-correlation analyses, suggesting a potential
reprocessing mechanism of the more variable X-ray emission in the UV emitter
of this source. One challenge faced by the Swift monitoring program is that
the lightcurves are not continuously sampled and hence standard Fourier
techniques cannot be applied. This uneven sampling of the lightcurves is
imposed by limited telescope time.

In the context of cross-correlation analysis, methods have been developed to
address the problem of unevenly-sampled lightcurves. In 2000_Reynolds, the
method of 1992_Press is extended to interpolate the lightcurve gaps using a
model of the covariance function, or equivalently the power spectrum, of the
lightcurve. In 1998_Bond; 2010_Miller; 2013_Zoghbi a maximum likelihood
approach is taken to fit models of the lightcurve power spectra which accounts
for the correlation between the lightcurves. In this paper we focus on a
relatively new approach to tackle unevenly-sampled lightcurves.

Gaussian processes (gps) confer a Bayesian nonparametric framework to model
general time series data (2013_Roberts; 2015_Tobar) and have proven effective
in tasks such as periodicity detection (2016_Durrande) and spectral density
estimation (2018_Tobar). More broadly gps have recently demonstrated modelling
success across a wide range of spatial and temporal application domains
including robotics (2011_Deisenroth; 2020_Greeff), Bayesian optimisation
(2015_Shahriari; 2020_Grosnit; 2020_Rivers) as well as areas of the natural
sciences such as molecular machine learning (2021_Nigam; 2020_Griffiths;
2020_flowmo; 2021_Griffiths; 2021_Hase) and genetics (2020_Moss). In the
context of astrophysics there is a recent trend favouring nonparametric models
such as gps due to the flexiblity afforded when specifying the underlying data
modelling assumptions. Applications have arisen in lightcurve modelling
(2021_Luger1; 2021_Luger_next), continuous-time autoregressive moving average
(CARMA) processes (2021_Yu), modelling stellar activity signals in radial
velocity data (2015_Rajpaul), lightcurve detrending (2016_Aigrain), learning
imbalances for variable star classification (2020_Lyon), inferring stellar
rotation periods (2018_Angus), estimating the dayside temperatures of hot
Jupiters (2019_Pass), exoplanet detection (2017_Jones; 2017_Czekala;
2020_Gordon; 2020_Langellier), spectral modelling (2012_Gibson; 2018_Nikolov;
2020_Diamond), as well as blazar variability studies (2015_Karamanavis;
2017_Karamanavis; 2020_Covino; 2020_Yang).

It has recently been demonstrated in lightcurve simulations by 2019_Wilkins
that a gp framework can compute time lags associated with X-ray reverberation
from the accretion disc that are longer and observed at lower frequencies than
can be measured by applying standard Fourier transform techniques to the
longest available continuous segments. It is for this principal reason that
gps are employed for the timing analysis in this chapter. Further desirable
facets of gps include the fact that, unlike parametric models, they do not
make strong assumptions about the shape of the underlying light curve
(2012_Wang_light). Additionally, Bayesian model selection may be performed at
the level of the covariance function or kernel allowing the quantitative
comparison of different models of the lightcurve power spectrum. Finally in
the cross-correlation analysis, a weaker modelling assumption is made than in
2013_Zoghbi in treating the X-ray and UV lightcurves as being independent
(2019_Wilkins).

The remainder of this chapter is outlined as follows: In Section 3.4
procedures used to fit gps to the X-ray and UVW2 bands are described,
including aspects such as identification of the flux distribution,
consideration of measurement noise as well as a simulation study to determine
the appropriate kernels. In Section 3.5 the structure functions of the gp-
interpolated lightcurves are compared with the observational structure
functions from gallo18. In Section 3.6 a cross-correlation analysis of the
X-ray and UVW2 bands is presented using the gp-interpolated lightcurves.
Finally, in Section 3.7 concluding remarks are provided about the discrepancy
between the observational and gp-derived structure functions as well as the
implications of the cross-correlation analysis, namely that the broad lag
features suggest an extended emission region of the disc in Mrk 335 during the
reverberation process. All code for reproducing the analysis is available at
https://github.com/Ryan-Rhys/Mrk_335.

###  3.4 Modelling Markarian 335

This Chapter considers the Swift X-ray and UVW2 lightcurves in time bins of
one day. The reader is referred to gallo18 for details of the data reduction
processes. The observational measurements used in this work run from
$54327-58626$ modified Julian days and comprise $509$ data points for the
X-ray band and $498$ data points for the UVW2 band. The latest UVOT
sensitivity calibration file (‘swusenscorr20041120v006.fits’) was considered
so as to account for the sensitivity loss with time in the UVW2 band333The
most up-to-date calibration files:
https://heasarc.gsfc.nasa.gov/docs/heasarc/caldb/swift. Only UVW2 data
collected by UVOT is considered because the UVW2 filter was most frequently
used in the archival observations..

####  3.4.1 Identifying the Flux Distribution

In order to assess the applicability of gps in modelling the flux distribution
of the X-ray and UVW2 bands of Mrk 335, a series of graphical distribution
tests were performed to determine the sample distribution. The histograms of
the log count rates for the X-ray, and flux for the UV bands, of Mrk 335 are
shown in Figure 3.2. The histograms show that the distribution of the UVW2
flux is approximately Gaussian-distributed whereas the X-ray count rate
distribution appears to be log-Gaussian distributed in line with the general
observation of 2005_Uttley that fluxes from accreting black holes tend to
follow log-Gaussian distributions. Further graphical distribution tests based
on probability-probability (PP) plots and empirical cumulative distribution
functions (ECDFs) are provided in Appendix A.1.

Furthermore, following 2019_Wilkins a Kolmogorov-Smirnov test for goodness-of-
fit was performed, where the null hypothesis is that the sample was drawn from
a Gaussian distribution. For the UVW2 flux values a p-value of $0.164$ was
obtained. For the raw X-ray count rates a p-value of $1.017e^{-20}$ was
obtained, and a p-value of $0.028$ for the log-transformed X-ray count rates.
As such, the null hypothesis that either UVW2 flux or log-transformed X-ray
count rates are drawn from a Gaussian distribution cannot be rejected at the
$1\%$ level of significance. The null hypothesis may however be rejected in
the case of the raw X-ray count rates, providing evidence that the raw X-ray
count rates should be log-transformed in order to be well-modelled by a
Gaussian distribution. As such, the raw X-ray count rates were log-transformed
and the UVW2 flux values were left unchanged.

(a) X-Ray Log Count Rates

(b) UVW2 Flux

Fig. 3.2: Histograms of the observed Swift X-ray log count rate and UVW2 flux
overlaid with Gaussian kernel density estimates. The raw UVW2 flux values have
been scaled by $1e^{14}$.

####  3.4.2 Noise Considerations

As noted by 2019_Wilkins fitting a gp to the logarithm of the count rate is
appropriate only in the limit of a large signal-to-noise ratio. In the case of
Mrk 335, the Poisson (shot) noise intrinsic to the photon detectors used to
obtain the flux measurements is over an order of magnitude smaller than the
flux measurement itself. As such the choice of the log-Gaussian process would
appear to be justified.

####  3.4.3 Lightcurve Simulations

Fig. 3.3: Residual plot. The normalised RSS metric is the sum of squared
residuals divided by the total number of discretised points (4390) comprising
the simulated lightcurve. A residual in this case represents the difference
between the Gaussian process predictive mean and the ground truth value of the
simulated lightcurve.

A simulation study was undertaken to quantitatively assess the abilities of
different kernels to interpolate gapped simulated lightcurves. Observational
power spectral densities (PSDs) of AGN are well-described by (broken) power
laws (2004_Mchardy). As such, the simulations employed a power law PSD with
index fit to the observational data. The goals with the study are twofold:
Firstly, although one cannot be sure of the true PSD for the observational
data, it is hoped that the simulations may afford a good proxy for identifying
performant kernels based on the fact that AGN typically exhibit power law-like
PSDs and secondly, it is desirable to asses the correlation between a kernel’s
ability to reconstruct the full simulated lightcurve and its marginal
likelihood value for the gapped data on which it is trained. If there is a
correlation, the marginal likelihood may be used as a metric for identifying
the appropriate kernel on the observational data.

One thousand simulated light curves with gaps were generated for the Mrk 335
X-ray and UV bands using the method of 1987_Davies, first applied in
astrophysics by 1995_Timmer. For each lightcurve there is access to the ground
truth functional form of the lightcurve before the introduction of gaps.
Computationally, the ground truth lightcurve was evaluated on a fine, discrete
grid of $4390$ time points whereas the gapped lightcurves were evaluated on a
coarser, unevenly-spaced grid of $498$ time points for the UV simulations and
$509$ time points for the X-ray simulations in line with the number of
observational data points. How well each gp kernel performs in recovering the
ground truth lightcurve was then quantified by measuring the normalised
residual sum of squared errors,

$\text{RSS}=\frac{1}{N}\sum_{i=1}^{N}(f(t_{i})-y_{i})^{2},$ (3.2)

where $f(t_{i})$ is the gp prediction at grid point $t_{i}$ and $y_{i}$ is the
true simulated count rate value. The RSS values were averaged over the one
thousand simulated lightcurves. An illustration of the RSS metric is provided
in Figure 3.3. In addition, the averaged negative log marginal likelihood
(NLML) values are computed for each kernel. Kernel hyperparameters were
selected via optimisation of the NLML using the SciPy optimiser of GPflow
(GPflow). The jitter level was fixed at 0.001, a small positive number to
ensure numerical stability. The output values (flux or the logarithm of the
count rate) were standardised according to their empirical mean and standard
deviation. A constant mean function set to the empirical mean of the data
following standardisation was employed.

The results of the simulation study are reported in Table 3.1. The NLML values
show correlation with RSS, thus providing evidence that NLML is an appropriate
metric for determining the gp kernel for the real observational data (for
which the ground truth lightcurve is of course not available). A paired t-test
was conducted to determine whether the RSS results were significant in terms
of identifying the best kernel. For the X-ray simulations, a t-statistic of
$9$ was obtained corresponding to a two-sided p-value of $5^{-20}$. For the
UVW2 simulations, a t-statistic of $-22$ was obtained corresponding to a two-
sided p-value of $9^{-85}$. As such, the null hypothesis that the performance
discrepancy between kernels on the RSS metric is due to chance variation
across $1000$ simulations, may be rejected at the $1\%$ level of significance.
Further rationalisation for why the top two performing kernels in the
simulation study are the Matérn $\frac{1}{2}$ and RQ kernels is offered in
Appendix A.2 (plotted by Douglas Buisson).

Table 3.1: Performance comparison of kernels based on the NLML on the simulated gapped X-ray and UV lightcurves and normalised residual sum of squared errors (RSS) on the ground truth simulated lightcurves. The mean NLML and RSS across 1000 simulations are reported with the standard error. UVW2 RSS values have an exponent of $-30$. Kernel | NLML | RSS
---|---|---
X-Ray |  |
Matérn$\frac{1}{2}$ | $\textbf{180.2}\pm\>\textbf{3.8}$ | $0.121\pm\>0.002$
Matérn$\frac{3}{2}$ | $420.7\pm\>3.3$ | $0.309\pm\>0.003$
Matérn$\frac{5}{2}$ | $523.5\pm\>2.9$ | $0.374\pm\>0.003$
Rational Quadratic | $\textbf{184.2}\pm\>\textbf{3.6}$ | $\textbf{0.117}\pm\>\textbf{0.002}$
Squared Exponential | $632.1\pm\>1.5$ | $0.554\pm\>0.004$
UVW2 |  |
Matérn$\frac{1}{2}$ | $\textbf{-399.0}\pm\>\textbf{5.2}$ | $\textbf{2.9}\pm\>\textbf{0.08}$
Matérn$\frac{3}{2}$ | $-298.3\pm\>6.0$ | $7.9\pm\>0.25$
Matérn$\frac{5}{2}$ | $-219.6\pm\>6.5$ | $17.0\pm\>0.41$
Rational Quadratic | $-349.2\pm\>5.4$ | $3.4\pm\>0.09$
Squared Exponential | $-65.0\pm\>7.4$ | $32.8\pm\>0.55$

####  3.4.4 Modelling Markarian 335 with Gaussian Processes

The fits to the observational data for the UVW2 and X-ray bands are shown in
Figure 3.4 and Figure 3.5 respectively. In an analogous fashion to the
simulation experiments five stationary kernels were evaluated: Matérn
$\frac{1}{2}$, Matérn $\frac{3}{2}$, Matérn $\frac{5}{2}$, rational quadratic
and squared exponential. The two kernels, rational quadratic and Matérn
$\frac{1}{2}$, which performed best in the simulation study in their abilities
to model power law-like PSDs are displayed. These kernels also have the most
favourable values under the NLML metric for the observational data. A constant
mean function set to the empirical mean of the data following standardisation
was again used. All kernel hyperparameters were optimised under the marginal
likelihood save for the noise level which was fixed to a constant value in the
standardised space. This constant noise value is computed by dividing the mean
output value in the standardised space by the mean signal-to-noise ratio in
the original space.

(a) UV Band | Matérn $\frac{1}{2}$ | Mean.

(b) UV Band | Matérn $\frac{1}{2}$ | Sample.

(c) UV Band | Rational Quadratic | Mean.

(d) UV Band | Rational Quadratic | Sample

Fig. 3.4: gp lightcurves for the UVW2 band. The shaded regions denote the gp
95% confidence interval. Both the gp mean and a sample from the gp posterior
are shown in separate plots. The insets are included to highlight the
variability of the fit.

(a) X-ray Band | Matérn $\frac{1}{2}$ | Mean.

(b) X-ray Band | Matérn $\frac{1}{2}$ | Sample.

(c) X-ray Band | Rational Quadratic | Mean.

(d) X-ray Band | Rational Quadratic | Sample

Fig. 3.5: gp lightcurves for the X-ray band. The shaded regions denote the gp
95% confidence interval. Both the gp mean and a sample from the gp posterior
are shown in separate plots. The insets are included to highlight the
variability of the fit.

###  3.5 Structure Function Analysis

Ideally one would like to examine the PSD of the gp fits to the observational
data. The PSD characterises the distribution of power over frequencies of a
given emission band and properties of the PSD can be linked to underlying
physical processes in the accretion disc. Computation of the PSD, while
possible, can be complicated by the uneven sampling of the observational data,
leading previous studies to instead perform a structure function analysis on
the Mrk 335 data (gallo18). While it is possible to extract the PSD from the
learned kernel (2019_Wilkins), in this work a structure function analysis of
the gp lightcurves was performed in order to compare directly against the
results of gallo18. The method described in 1985_Simonetti; 1992_Hughes;
1996_di_Clemente; 2001_Collier; gallo18 was followed. The binned structure
function is defined as:

$\text{SF}(\tau)=\frac{1}{N(\tau)}\sum_{i}\>[f(t_{i})-f(t_{i}+\tau)]^{2},$
(3.3)

where $\tau=t_{j}-t_{i}$ is the distance between pairs of points $i$ and $j$
such that $t_{j}>t_{i}$. The structure function is binned according to $\tau$
where the centres of each bin are given by $\tau_{i}=(i-\frac{1}{2})\delta$.
In this instance, $\delta$ is the structure function resolution. The same
$\delta$ as in gallo18 was used, namely 5.3 days for the structure function
computation over both the X-ray and UVW2 bands. $f(t_{i})$ gives the count
rate value at time point $t_{i}$ and $N(\tau)$ is the number of structure
function pairs in each bin $i$ with centre $\tau_{i}$. Accounting for
measurement noise by subtracting twice the mean noise variance from each
structure function bin, as performed in gallo18 was found to have negligible
effect on the gp structure functions and so was ignored. As in gallo18, the
structure function values were normalised by the global lightcurve variance.

The gp structure functions for the interpolated lightcurves are shown in
Figure 3.6. The $1\sigma$ gp error bars were obtained by computing the
structure function over 50 samples from the gp posterior. Each sample gives
rise to highly similar structure functions and so the errors are not visible
on the plot. The structure functions computed from the observational data, 509
and 498 data points for the X-ray and UV bands of Mrk 335 respectively, are
included for reference. In contrast to the gp structure function errorbars, in
the case of the observational data the error bars are computed as
$\frac{\sigma_{i}}{\sqrt{(}\frac{N_{i}}{2})}$ where $\sigma_{i}$ is the noise
standard deviation in bin $i$ and $N_{i}$ is the number of pairs in bin $i$.

The gp structure functions are compared against the observational structure
functions in Figure 3.6. In addition, the broken power law fits to the gp
structure functions are plotted, the parameters of which are given in Table
3.2. In the UVW2 band, both gp kernels yield structure functions possessing a
consistent break point at ca. 125 days. In the X-ray band the Matérn
$\frac{1}{2}$ kernel yields a break point at 66 days whereas the rational
quadratic kernel fit yields an unbroken power law. Given the discrepancy
between gp kernels, definite evidence for a break in the X-ray power law is
not found.

Of particular interest is whether the dip in the X-ray structure function is
an expected feature of the latent lightcurve or a measurement artefact arising
from uneven sampling. In order to assess the potential for the dip to be a
sampling artefact, simulations were performed using the Timmer and König
algorithm from Subsection 3.4.3. In this case structure functions of gapped
lightcurves were computed and compared against structure functions derived
from the ground truth lightcurves with no gaps. One representative simulation
is depicted in Figure 3.7. In this instance a similar dip to that found in the
observational data is observed in the X-ray band simulation. This highlights
the possiblity that the dip seen in the observational X-ray structure function
is a sampling artefact arising from gaps in the lightcurve.

(a) Matérn $\frac{1}{2}$ UVW2

(b) Rational Quadratic UVW2

(c) Matérn $\frac{1}{2}$ X-ray

(d) Rational Quadratic X-ray

Fig. 3.6: Comparison of observational and gp structure functions. The gp structure functions are consistent with those calculated from the observational data in the non-noise dominated regions. The dip at ca. $200$ days in the observational X-ray structure function is potentially a sampling artefact as demonstrated by simulation in Figure 3.7. Table 3.2: Parameters for the broken power law fits to the gp structure functions. $\alpha_{1}$ and $\alpha_{2}$ are the indices for the power law before and after the break point $\tau_{char}$. The break point $\tau_{char}$ is reported in days. Errors were computed using 200 bootstrap samples of the data points corresponding to the gp structure functions. The X-ray rational quadratic structure function was fit using a power law and as such only has a single index as a parameter. The Astropy library (astropy_1; astropy_2) was used to compute the (broken) power law fits using the simplex algorithm and least squares statistic for optimisation with the gp structure function uncertainties used as weights in the fitting. Waveband | Kernel | $\alpha_{1}$ | $\alpha_{2}$ | $\tau_{char}$
---|---|---|---|---
UVW2 | Matérn $\frac{1}{2}$ | $-0.72\pm\>0.03$ | $-0.26\pm\>0.01$ | $127\pm\>8$
UVW2 | Rational Quadratic | $-0.62\pm\>0.01$ | $-0.28\pm\>0.01$ | $125\pm\>5$
X-ray | Matérn $\frac{1}{2}$ | $-0.29\pm\>0.03$ | $-0.07\pm\>0.004$ | $66\pm\>8$
X-ray | Rational Quadratic | $-0.21\pm\>0.002$ | N/A | N/A

(a) Observational

(b) gp | Matérn $\frac{1}{2}$

Fig. 3.7: Structure function simulations. Pseudo-observational lightcurves are
produced by introducing gaps into the simulated lightcurves. The structure
function for the gapped lightcurve is shown in red in (a) whereas the
structure function of the gp interpolation is shown in red in (b). Both
structure functions are compared against the ground truth structure function
obtained from the full simulated ground truth lightcurve (no gaps). The dips
at $\tau=200\>\text{days}$ and $\tau=400\>\text{days}$ in the structure
function derived from the gapped observational simulation in 3.7(a) are
artefacts of the uneven sampling.

###  3.6 Lag and Coherence

Fig. 3.8: The coherence and lag spectra for Mrk 335, calculated by using 1000
pairs of gp lightcurve samples fit to the observed lightcurves. The error bars
are the standard errors of the corresponding measurements for the 1000
samples. Different panels are for different kernels. Positive lags imply that
the X-ray band leads the UVW2 band. Spectra plotted by Jiachen Jiang.

In this section, the coherence between the UVW2 and X-ray emission from Mrk
335 is computed in search of evidence of lag features in the Fourier frequency
domain. The coherence and lag spectra were estimated from one thousand pairs
of UVW2 and X-ray gp lightcurve samples drawn from the gp posterior for each
kernel. The lags are defined as the phase lags divided by the corresponding
Fourier frequency. A similar approach has been used in other disciplines (e.g.
fabian09; kara13). Both Matérn $\frac{1}{2}$ and rational quadratic kernels
are considered. The results are shown in Figure 3.8. These spectra were
plotted by Jiachen Jiang. Positive lags imply that the X-ray variability leads
the UVW2 variability. The error bars in the figure are the standard errors of
the corresponding measurements for the one thousand samples.

The coherence between the UVW2 and X-ray emission decreases with frequency,
suggesting more coherent variability at lower frequency. Positive lag features
are shown at the low frequencies in the range $f=$ 0.005–0.025 d-1. The
absolute value of the lag at $f=0.0039\pm 0.0014$ d-1 is estimated to be
$19\pm 22$ days for the Matern $\frac{1}{2}$ kernel applied to both
lightcurves and $29\pm 19$ days for the rational quadratic kernel, however
both measurements are consistent with zero lag in the $2\sigma$ uncertainty
range.

Tentative evidence of a shorter time lag at a higher frequency of $f=0.018\pm
0.006$ d-1 is also found. The longer lag feature at a lower frequency would
correspond to a more extended emission region while the shorter lag feature at
a higher frequency would correspond to a more compact region. This could be
explained by the presence of an extended UV emission region on the disc where
reverberation happens.

Given that the lags are consistent with zero lag within $2\sigma$ uncertainty
ranges, it is concluded that only tentative evidence for a broad lag feature
is found by applying gps to the UVW2 and X-ray lightcurves of Mrk 335.
Previous attempts to identify lags between two wavelengths of Mrk 335 based on
cross-correlation analysis in the time domain suggests similar results (e.g.
gallo18).

###  3.7 Conclusions

Following the interpolation of the unevenly-sampled lightcurves of Mrk 335
using gps, tentative evidence for broad lag features is found in the Fourier
frequency domain. The magnitude of the lags is consistent with previous cross-
correlation analyses. In addition, the broad lag features might suggest an
extended emission region e.g. of the disc in Mrk 335 during the reverberation
processes. If the corona is compact within 5 $R_{\rm g}$ in Mrk 335
(2015_Wilkins), our data suggest a possibly wide range of UVW2 emission radii.

The structure functions computed from the gp-interpolated lightcurves are
consistent with those derived from the observational data and furthermore,
illicit potential insights into the properties of the latent lightcurves. In
particular, it is shown through a simulation study that it is possible that
dips in the X-ray structure function may be produced by sampling artefacts
arising from gaps in the lightcurve. In contrast, the gp structure functions
show no dips. While this is not proof that the dip in the observational X-ray
structure function is due to a sampling artefact, it does allude to the
possibility. The UVW2 gp structure functions do not exhibit strong dependence
on the choice of kernel with both Matérn $\frac{1}{2}$ and rational quadratic
showing up a broken power law with breaks at 139 and 155 days respectively.
The X-ray structure functions however do show up differences between kernels
with the rational quadratic kernel predicting a power law and the Matérn 1/2
kernel predicting a broken power law.

From the gp modelling perspective, the ability to carry out Bayesian model
selection affords a quantitative means of comparing analytic kernels under the
marginal likelihood. It may be possible to incorporate further flexibility
into the fitting procedure by making use of more sophisticated methods of
kernel design (2014_Duvenaud) to allow the assessment of fits of sums and
products of analytic kernels or by leveraging advances in transforming gp
priors via Deep gps (2013_Damianou) or normalising flows (2020_Maronas). Such
approaches could be validated using simulation studies. Additionally,
modelling the cross-correlation using multioutput gps (2020_de_Wolff) may be
an interesting avenue for comparison against the approach taken here. Lastly,
Bayesian spectral density estimation (2018_Tobar) may afford further
flexibility through nonparametric modelling of the PSD in addition to
nonparametric modelling of the lightcurve in the time domain. These
improvements in Bayesian modelling machinery may help to minimise model
misspecification and as such, enable more robust inferences to be made about
the functional forms of the latent lightcurves.

##  Chapter 4 Modelling Molecules with Gaussian Processes

Status: Accepted as Griffiths, RR., Klarner, L., Moss, HB., Ravuri, A.,
Truong, S., Rankovic, B., Du, Y., Jamasb, A., Schwartz, J., Tripp, A., Kell,
G. Bourached, A., Chan, A., Moss, J., Chengzhi, G, Lee AA, Schwaller P., Tang,
J. GAUCHE: A library for Gaussian processes in chemistry. ICML Workshop on
AI4Science, 2022.

###  4.1 Background on Molecular Machine Learning

The chapter begins with a self-contained background on molecular machine
learning required to contextualise the findings of this chapter.

Fig. 4.1: The design-make-test cycle of molecular discovery.

Although application domains of machine learning in the molecular sciences are
constantly expanding (2020_Struble), an important subset of applications can
be taxonimised according to the role they play in enhancing the design-make-
test cycle (2012_Plowright) of molecular discovery campaigns, illustrated in
Figure 4.1. Molecule generation (2018_Design; 2020_Griffiths; 2017_Grammar;
2020_Jin; 2022_Du; 2022_Gao) is concerned with designing novel molecules using
generative models such as variational autoencoders (VAEs) (2014_Kingma) and
generative adversarial networks (GANs) (2014_Goodfellow). Chemical reaction
prediction (2019_Schwaller), reaction planning (2018_Coley_new), and synthesis
design (2022_Schwaller), illustrated in Figure 4.2, are focussed on improving
the throughput of the make stage of the design-make-test cycle by using
machine learning to suggest synthetic pathways to target molecules.

The molecular machine learning tasks considered in this thesis, are molecular
property prediction and chemical reaction optimisation. Molecular property
prediction is concerned with the design phase of the design-make-test cycle in
so far as it allows molecules to be prioritised for laboratory synthesis.
Chemical reaction optimisation, on the other hand, is concerned with the make
phase as it is a means of improving the yield of chemical reactions. Both
applications will be described in detail next.

Fig. 4.2: A, B and C are starting materials and P is the product of the
chemical reaction. Reaction planning involves finding a set of reagents to
transform a starting material into a product. Reaction prediction involves
predicting the product given a set of reactants. Synthesis design involves
working backwards from the product towards a set of reactants and reagents.
Machine learning-based solutions to all of these tasks would yield a blueprint
for a chemist to follow in synthesising a novel molecule.

####  4.1.1 Molecular Property Prediction

The laboratory synthesis of a novel molecule is a highly time-intensive
process. As such, there has been a great deal of interest in computational
approaches to prioritising molecules drawn from vast molecular databases in a
process known as high-throughput virtual screening (2015_Knapp). In this
fashion, theoretical techniques can be used to winnow the database down to a
few promising candidates for experimental chemists to follow up on. Before the
advent of machine learning, the dominant approach was to use first principles
quantum chemical calculations to compute molecular properties. Below, one such
first principles method is reviewed, density functional theory (DFT), which is
compared against machine learning model performance in Chapter 5, before
discussing machine learning approaches to molecular property prediction.

##### Density Functional Theory

DFT is a method of modelling the electronic structure of many-body systems
(brazdova_atomistic_2013), and has been applied across problems in physics,
chemistry, biology, and materials science (becke_perspective_2014). DFT is an
_ab initio_ , or first principles computational method because physical
constants are the only inputs to calculations based on the postulates of
quantum mechanics (leach_molecular_2001). Since the inception of DFT in
1964-1965, Kohn-Sham DFT (KS-DFT) has been one of the most frequently applied
electronic structure methods (becke_perspective_2014).

KS-DFT makes use of the Hohenberg-Kohn theorems
(hohenberg_inhomogeneous_1964), a trial electron density, and a self-
consistency scheme. KS-DFT executes a computational loop by starting with a
trial density, solving the Kohn-Sham equations, and obtaining the single
electron wavefunctions for the trial density; in the next step, an electron
density may be computed. If the computed density is consistent i.e. within a
tolerance threshold of the trial density, the theoretical ground state density
has been identified. If the densities are not consistent, however, the
computed density is taken as the new trial density, and the iterative loop
continues to be executed until the tolerance threshold is met. The accuracy of
DFT calculations, with exchange and correlation functionals, can be very high,
yet may exhibit significant fluctuations with the choice of functional,
pseudopotential, basis sets and cutoff energy (howard_assessing_2015).
Furthermore, these quantities are not always trivial to optimise.

##### Time-Dependent Density Functional Theory

Time-Dependent Density Functional Theory (TD-DFT) is a time-dependent analogue
of DFT based on the Runge-Gross (RG) theorem in place of the Hohenberg-Kohn
theorems (runge_density-functional_1984). The RG theorem states that a unique
delineation exists between the time-dependent electron density and the time-
dependent external potential. As such, a computational, time-dependent Kohn-
Sham system may be implemented (van_leeuwen_causality_1998) in a similar
fashion to KS-DFT. When TD-DFT has been used together with a linear response
theory (ullrich_time-dependent_2012), it has enjoyed success in the
calculation of electromagnetic spectra of medium and large molecules
(casida_progress_2012; burke_time-dependent_2005). A relevant application of
TD-DFT in this thesis is the computation of the
$\pi-\pi{{}^{*}}/\emph{n}-\pi{{}^{*}}$ electronic transitions wavelengths for
photoswitch molecules in Chapter 5.

##### Machine Learning Approaches

In contrast to first principles methods such as DFT, machine learning
approaches seek to carry out data-driven prediction. A key issue in data-
driven molecular property prediction is how best to featurise molecules. This
problem is commonly referred to as choosing a molecular representation. While
a great many base representations of molecules exist (2022_Wigh), some of the
most popular featurisations include graph-based, string-based and fingerprint
representations. The field of molecular representation learning is concerned
with learning representations on top of these base representations e.g.
(2015_Duvenaud) typically via deep learning. Below, a brief review is provided
of the molecular representations used in this thesis.

###### Graphs:

Molecules may be represented as an undirected, labeled graph
$\mathcal{G}=(\mathcal{V},\mathcal{E})$ where vertices,
$\mathcal{V}=\\{v_{1},\ldots,v_{N}\\}$, represent the atoms of an $N$-atom
molecule and edges, $\mathcal{E}\subset\mathcal{V}\times\mathcal{V}$,
represent covalent bonds between these atoms. Additional information may be
incorporated in the form of vertex and edge labels
$\mathcal{L}:\mathcal{V}\times\mathcal{E}\to\Sigma_{V}\times\Sigma_{E}$.
Common label spaces including attributes such as atom types (i.e. hydrogen,
carbon) as vertex labels and bond orders (i.e. single, double) as edge labels.

###### Strings:

The Simplified Molecular-Input Line-Entry System (SMILES) is a text-based
representation of molecules (1987_Anderson; 1988_Weininger), examples of which
are given in Figure 4.3. Self-Referencing Embedded Strings (SELFIES)
(2020_Krenn) is an alternative string representation to SMILES such that a
bijective mapping exists between a SELFIES string and a molecule.

Fig. 4.3: SMILES strings for structurally similar molecules. Similarity is
encoded in the string through common contiguous subsequences (black). Local
differences are highlighted in red. Note the molecules are chosen solely for
the purposes of illustrating the SMILES syntax.

###### Fingerprints:

Molecular fingerprints were first introduced for chemical database
substructure searching (1993_Christie), but were later repurposed for
similarity searching (1990_Johnson), clustering (1997_McGregor) and
classification (2017_Breiman). Extended Connectivity FingerPrints (ECFP)
(2010_Rogers) were introduced as part of the Pipeline project (2006_Hassan)
with the explicit goal of capturing features relevant for molecular property
prediction (2004_Xia). ECFP fingerprints operate by assigning initial numeric
identifiers to each atom in a molecule. These identifiers are subsequently
updated in an iterative fashion based on the identifiers of their neighbours.
The number of iterations corresponds to half the diameter of the fingerprint
and the naming convention reflects this. For example, ECFP6 fingerprints have
a diameter of 6, meaning that 3 iterations of atom identifier reassignment are
performed. Each level of iteration appends substructural features of
increasing non-locality to an array and the array is then hashed to a bit
vector reflecting the presence of absence of those substructures in the
molecule.

For property prediction applications a radius of 3 or 4 is recommended. A
radius of 3 is used for all experiments in the thesis. Fragment descriptors
are also used, which are count vectors, each component of which indicates the
number of a certain functional group present in a molecule. For example row 1
of the count vector could be an integer representing the number of aliphatic
hydroxl groups present in the molecule. Both fingerprint and fragment features
computed using RDKit are made use of (rdkit), as well as the concatenation of
the fingerprint and fragment feature vectors, a representation termed
fragprints which has shown strong empirical performance. Example
representations $\mathbf{x_{f}}$, for fingerprints and $\mathbf{x_{fr}}$, for
fragments might be

$\displaystyle\mathbf{x_{f}}$ $\displaystyle=\begin{bmatrix}1\\\ 0\\\
\vdots\\\ 1\end{bmatrix},\>\>\>\mathbf{x_{fr}}=\begin{bmatrix}3\\\ 0\\\
\vdots\\\ 2\end{bmatrix}.$

##### Bayesian Optimisation and Active Learning for Molecules

bo and active learning hold particular promise for accelerating high-
throughput virtual screening efforts for molecules (2017_Knapp; 2020_Pyzer).
The idea in this instance is that data-efficient molecular property
optimisation can be performed in a bo/active learning loop featuring
laboratory synthesis of the suggested molecules.

####  4.1.2 Chemical Reaction Optimisation

The yield of a chemical reaction is an important consideration for large-scale
production of a desired molecule. The percentage yield, of a chemical reaction
may be defined as

$\text{Percentage Yield}=\frac{\text{Actual Yield}}{\text{Theoretical
Yield}}.$ (4.1)

The Theoretical Yield assumes a flawless chemical reaction in which all
starting materials are converted to the desired product. In practice, chemical
reactions are not perfectly efficient due to factors such as reverse
reactions, in which the reactants and products exist in a state of chemical
equilibrium, as well as competing chemical reactions that form unwanted side
products. Factors such as the temperature of the reaction, the concentration
of reactant species, the choice of solvent as well as the presence of
reagents, molecular species that enhance the reaction without contributing
atoms to the product, are all determinants of the efficiency of the reaction
and hence the Actual Yield. Optimising such reaction parameters has recently
been tackled by machine learning approaches such as bo (2021_Shields). Common
reaction representations are detailed next.

##### Chemical Reaction Representations

A chemical reaction comprises reactants and reagents that transform into one
or more products together with reaction parameters such as temperature and
concentration. The reactants and reagents are instances of molecular species
which play different roles in the reaction. By means of illustration, the
high-throughput experiments by ahneman2018predicting on Buchwald-Hartwig
reactions feature a reaction design space consisting of 15 aryl and heteroaryl
halides, 4 Buchwald ligands, 3 bases, and 23 isoxazole additives. Below,
various means of featurising the reactant and reagent components of a chemical
reaction are introduced.

###### Concatenated molecular representations:

If the number of reactant and reagent categories is constant, the molecular
representations discussed above can be used to encode reactants and reagents.
The vectors for the individual reaction components may then be concatenated to
build the reaction representation (ahneman2018predicting;
sandfort2020structure). An additional concatenated representation, is the one-
hot-encoding (OHE) of the reaction categories where bits indicate the presence
or absence of a particular reactant/reagent. In the Buchwald-Hartwig example
above, the OHE would describe which of the aryl halides, Buchwald ligands,
bases, and additives are used in the reaction, resulting in a 44-dimensional
bit vector (chuang2018comment).

###### Differential reaction fingerprints:

Inspired by the hand-engineered difference reaction fingerprints by
schneider2015development, probst2022reaction recently introduced the
differential reaction fingerprint (DRFP). This reaction fingerprint is
constructed by taking the symmetric difference of the sets containing the
molecular substructures on both sides of the reaction arrow. The size of the
reaction bit vector generated by DRFP is independent of the number of reaction
components.

###### Data-driven reaction fingerprints:

schwaller2021mapping described data-driven reaction fingerprints using
Transformer models such as BERT (devlin2018bert), trained in a supervised or
an unsupervised fashion on reaction SMILES. The Transformer models can then be
fine-tuned on the task of interest to learn more specific reaction
representations (schwaller2021prediction). These representations are
designated using the acronym RXNFP. As with the DRFP, the size of data-driven
reaction fingerprints is also independent of the number of reaction
components.

In the next section, the focus of this chapter is introduced.

###  4.2 Preface

This chapter introduces GAUCHE, a library for GAUssian processes in CHEmistry.
gps have long been a cornerstone of probabilistic machine learning, affording
particular advantages for uncertainty quantification (UQ) and bo. Extending
gps to chemical representations however is nontrivial, necessitating kernels
defined over structured inputs such as graphs, strings and bit vectors. By
defining such kernels in GAUCHE, the door is opened to powerful tools for UQ
and bo in chemistry. Motivated by scenarios frequently encountered in
experimental chemistry, applications for GAUCHE are showcased in molecule
discovery and chemical reaction optimisation.

###  4.3 Introduction

Early-stage scientific discovery is typically characterised by the small data
regime due to the limited availability of high-quality experimental data
(2018_Zhang). Much of the novelty of discovery relies on the fact that there
is a lot of knowledge to gain in the small data regime. By contrast, in the
big data regime, discovery offers diminishing returns as much of the knowledge
about the space of interest has already been acquired. As such, machine
learning methodologies that facilitate search in small data regimes such as bo
(2018_Bombarelli; 2020_Griffiths; 2021_Shields) and active learning (AL)
(2019_Zhang; 2021_Jablonka) have great potential to expedite the rate at which
performant molecules, molecular materials, chemical reactions and proteins are
discovered.

To date in molecular machine learning, BNNs have been the surrogate of choice
to produce the uncertainty estimates that underpin bo and AL (2019_Ryu;
2019_Zhang; 2020_Hwang; 2020_Scalia). For small datasets, however, DNNs are
often not the model of choice. Notably, certain deep learning experts have
voiced a preference for gps in the small data regime 2011_Bengio. Furthermore,
for bo, gps possess particularly advantageous properties; first, they admit
exact as opposed to approximate Bayesian inference and second, few of their
parameters need to be determined by hand. In the words of Sir David MacKay
2003_MacKay,

> "Gaussian processes are useful tools for automated tasks where fine tuning
> for each problem is not possible. We do not appear to sacrifice any
> performance for this simplicity.”

The iterative model refitting required in bo makes it a prime example of such
an automated task. Although BNN surrogates have been trialled for bo
(2015_Snoek; 2016_Springenberg), gps remain the model of choice as evidenced
by the results of the recent NeurIPS Black-Box Optimisation Competition
2021_Turner.

Training gps on molecular inputs is non-trivial however. Canonical
applications of gps assume continuous input spaces of low and fixed
dimensionality. The most popular molecular input representations are
SMILES/SELFIES strings (1987_Anderson; 1988_Weininger; 2020_Krenn),
fingerprints (2010_Rogers; probst2018probabilistic; capecchi2020one) and
graphs (2015_Duvenaud; 2016_Kearnes). Each of these input representations
poses problems for gps. SMILES strings have variable length, fingerprints are
high-dimensional and sparse bit vectors, while graphs are also a form of non-
continuous input. To construct a gp framework over molecules, GAUCHE provides
GPU-based implementations of kernels that operate on molecular inputs,
including string, fingerprint and graph kernels. Furthermore, GAUCHE includes
support for protein and chemical reaction representations and interfaces with
the GPyTorch (2018_Gardner) and BoTorch (2019_Balandat) libraries to
facilitate usage for advanced probabilistic modelling and bo. The detailed
contributions of this chapter may be summarised as:

  1. 1.

The introduction of a gp framework for molecules and chemical reactions.

  2. 2.

The provision of an open-source, GPU-enabled library building on GPyTorch
(2018_Gardner), BoTorch (2019_Balandat), and RDKit rdkit.

  3. 3.

The use of black box graph kernels, from GraKel, (2020_GraKel), is extended to
gp regression via a GPyTorch interface, along with a limited set of graph
kernels implemented in native GPyTorch to enable optimisation of the graph
kernel hyperparameters under the marginal likelihood.

  4. 4.

Benchmark experiments are conducted, evaluating the utility of the gp
framework on regression, UQ and bo tasks.

GAUCHE is made available at https://github.com/leojklarner/gauche and includes
tutorials to guide users through the tasks considered in this paper.

###  4.4 Related Work

General-purpose gp and bo libraries do not specifically cater for molecular
representations. Likewise, general-purpose molecular machine learning
libraries do not specifically consider gps and bo. Here, existing libraries
are reviewed, highlighting the niche GAUCHE fills in bridging the gp and
molecular machine learning communities. To date, there has been little work on