Log-Gaussian processes for AI-assisted TAS experiments

From a methodological perspective, we aim to discover an intensity function $i:X\to [0,\infty )$ on a rectangular set $X\subseteq R^n$, $n\in N$, with coordinates on a certain hyperplane in four-dimensional $Q$-$E$ space [sivia2011elementary]. As an example, Fig. 1a displays an intensity function defined on $X=[2.3,3.3]\times [2.5,5.5]\subseteq R^2$, where the two-dimensional hyperplane is spanned by the vector $(0,0,1)$ in $Q$ space with offset $(1,1,0)$ and energy transfer ($E$).

For directions in $Q$ space, we use relative lattice units (r.l.u.), while energy transfer $E$ is measured in micro-electron volts (meV). Note that, due to restrictions of an instrument, the intensity function might only be defined on a subset $X^{\ast }\subseteq X$ consisting of all measurement locations reachable. For more details on the TAS setting, we refer to the Supplementary Information (Sec. S1).

The intensity function $i$ is accessed by counting scattered neutrons on a detector device for a finite number of measurement locations $x\in X^{\ast }$ yielding noisy observations $I^{+}\left(x\right)\sim \text{Pois}(\lambda =i(x\left)\right)$. Note that detector counts are usually normalized by neutron counts on a monitor device, i. e., the corresponding unit of measurement is $\text{detector counts}/\left(M\text{monitor counts}\right)$, $M\in N$. Since Poisson distributions $\text{Pois}\left(\lambda \right)$ with a sufficiently large parameter $\lambda >0$ can be approximated by a normal distribution $N(\lambda ,\lambda )$, we assume that

where ${\eta }^{+}\sim N(0,1)$.

In the following, we repeatedly refer to “experiments” which are defined as a sequential collection of intensity observations.

Since demand for beam time at TAS is high but availability is limited, the goal of an approach is to perform the most informative experiments at the lowest possible cost. Beam time use can, for example, be optimized when excessive counting in uninformative regions like the background (Fig. 1b) is avoided but focused on regions of signal (Fig. 1c).

We quantify the benefit of an experiment $A$ by a benefit measure $\mu =\mu \left(A\right)\in R$ and its cost by a cost measure $c=c\left(A\right)\in R$. For now, it suffices to mention that cost is measured by experimental time (the time used for an experiment) and benefit is defined by reducing a relative weighted error between the target intensity function and a corresponding approximation constructed by the collected intensity observations. Note that, quantifying benefits this way, their computation is only possible in a synthetic setting with known intensity functions (as is our benchmark setting). Real neutron experiments do not meet this requirement and thus must be evaluated in a more qualitative way.

In our setting, an approach attempts to conduct an experiment $A$ with highest possible benefit using a given cost budget $C\ge 0$, i. e., it aims to maximize $\mu \left(A\right)$ while ensuring that $c\left(A\right)\le C$. The steps of a corresponding general experiment are given in Alg. 1.

Line 5 is most important and crucial for both cost and benefit of the experiment since it decides where to observe intensities, i. e., count neutrons, next.

From an algorithmic perspective, if we denote the current step of an experiment by $J\in N$, our approach implements the decision for the next measurement location by maximizing an objective function ${\varphi }_J:X^{\ast }\to R$. It balances an acquisition function ${\text{acq}}_J:X^{\ast }\to R$ and a cost function $c_J:X^{\ast }\to R$ that both depend on the current step $J$ and hence can change from step to step. The acquisition function ${\text{acq}}_J$ indicates the value of any $x\in X^{\ast }$ for improving the benefit of the experiment whereas the cost function $c_J$ quantifies the costs of moving the instrument axes from the current location $x_J\in X^{\ast }$ to $x$.

The metric $d:X^{\ast }\times X^{\ast }\to [0,\infty )$ used for our particular cost function

is formally specified in the Supplementary Information (Sec. S1, Eq. (10)). Our methodology concentrates on developing a useful acquisition function as the crucial component of our approach making it straightforward to find regions of signal.

We briefly describe here why log-Gaussian processes, our central methodological component, are suitable to identify regions of signal in a TAS experiment. Methodological details can be found in the Methods section.

Although the intensity function is not directly observable due to measurement noise (Eq. (1)), we aim to approximate it by the mean function of a log-Gaussian process

where $F$ is a Gaussian process. That is, after $J$ steps, we fit logarithmic intensity observations to $F$ yielding its posterior mean and variance function denoted by $m_J$ and ${\sigma }_J^2$, respectively. The acquisition function is then defined as the uncertainty of $I$ given by its posterior standard deviation, i. e.,

Observe the crucial detail that $m_J$ appears exponentially in this function which is the main reason why our approach is based on log-Gaussian processes. A posterior log-Gaussian process is thus able to find regions of signal just through maximizing its uncertainty. As illustration, regard Fig. 2 displaying the posterior of the Gaussian process $F$ together with logarithmic intensity observations (Fig. 2a) and the corresponding posterior of the log-Gaussian process $I$ (Fig. 2b).

Maximizing the acquisition function from Eq. (5) prioritizes regions with high intensities over regions with low intensities. This poses a problem when there are multiple signal regions with intensities of different magnitudes (Extended Data Fig. 6a). Indeed, measurement points are mainly placed in regions with higher intensities whereas regions with less signal are neglected (Extended Data Fig. 6b). In TAS, we are interested in each region of signal no matter of which intensity magnitude. We compensate for this potential problem by introducing an intensity threshold $\tau >0$ for observed intensities. That is, we truncate the observed intensities to a maximum of $\tau$ before fitting (Extended Data Fig. 6c). Consequently, measurement points get more evenly distributed among all signal regions (Extended Data Fig. 6d) since their placement is not biased due to large differences in their intensity values.

As another problem, intensity observations, in neutron experiments, contain background which is not part of the actual signal, i. e., even if there is no actual signal at a certain location, we might nonetheless observe a positive intensity there. If our approach does not compensate for regions of background, it might not recognize them as parasitic and hence consider them as regions of weak signal which potentially yields uninformative measurement points being placed there. Therefore, we subtract a background level $\gamma \in [0,\tau )$ from already threshold-adjusted intensity observations while ensuring a non-negative value.

Ultimately, the log-Gaussian process actually aims to approximate the modified intensity function

Observe that, by regarding adjusted intensity observations, the assumption of their noise being normally distributed is violated in general as their noise distribution is asymmetric. In particular, if $i\left(x_j\right)$ substantially exceeds the intensity threshold $\tau$, the distribution gets rather concentrated at $\tau$ and thus small in variance. We, however, do assume noise on adjusted intensities as if they were observed so without adjusting. This does not change the expected behaviour of getting a useful acquisition function but even ensures numerical stability since noise regularizes the computational problem of solving linear systems in GPR.

It remains to explain how we seek to compute suitable values for the background level $\gamma$ and the intensity threshold $\tau$ without knowing the intensity function. We estimate $\gamma =\gamma \left(A_0\right)$ and $\tau =\tau \left(A_0\right)$ by statistics of $J_0\in N$ initial measurement points

collected to initialize our approach. The arrangement of initial measurement locations (Extended Data Fig. 7) is specified later in the Methods section. While also the estimation of the background level $\gamma$ is described later, we define the estimation of the intensity threshold already at this point. This is necessary to understand the results of the benchmark, for which we set $\gamma =0$ manually due to its synthetic setting anyway. We define

where $m_{10}\ge \gamma \left(A_0\right)$ is the median of all intensity observations larger than their $9$-th decile and $\beta \in (0,1]$ is a parameter controlling the distinction between regions of strong and weak signals that needs to be set before starting an experiment. Note that this definition is, by the meaning of $m_{10}$, robust to outliers in $A_0$.

At SINQ (PSI), we investigated a sample of SnTe (tin telluride) in a real neutron experiment performed at the thermal TAS EIGER [stuhr2017thermal]. Our general aim is to reproduce known results from [li2014phonon, Fig. 1b] using our approach. Furthermore, we assess the robustness of the experimental result w.r.t. changes in the parameters for the background level and the intensity threshold (scenario 1) as well as challenge our approach using a coarser initialization grid on a modified domain $X$ with no initial measurement locations directly lying in a region of signal (scenario 2).

The software implementation of our approach communicates with the instrument control system NICOS (nicos-controls.org) which was configured on site. For the experimental setting at the instrument, we refer to the Supplementary Information (Sec. S2).

As mentioned, the benefit measure used for benchmarking, involving a known target intensity function, is not computable in this setting of a neutron experiment due to experimental artefacts like background and noise. We therefore evaluate the results of our approach in a more qualitative way for this experiment. Also, although the costs for moving the instrument axes contribute to the total experimental time, we do not consider them here for optimizing the objective function since they are approximately constant across the domain $X$.

For scenario 1, we adopt the setting from the original results [li2014phonon, Fig. 1b] which, in our context, means to investigate intensities on $X=[0,2]\times [2,12.4]$ along the vector $(1,1,0)$ in $Q$ space with offset $(0,0,3)$ and energy transfer. Initially, we performed measurements in a conventional mode for reference, i. e., we mapped $X$ with a grid of 11 columns containing 27 measurement points each (bottom row in Fig. 3).

These measurements were arranged in columns from bottom to top, i. e., from low to high energies, and from left to right, i. e., from small to large values for the coordinate along the vector in $Q$ space.

A direct comparison with our approach would put this mapping mode in disadvantage since its total experimental time could have been spent more efficiently. For this reason, the approach that we eventually take for comparison, the grid approach (top row in Fig. 3), takes the measurement points from the mapping mode but changes their order into four stages I-IV (Extended Data Fig. 8a-d). The first stage is from bottom to top and from left to right but only consists of every other grid point on both axes. The second stage is again from bottom to top but from right to left and also consists only of every other grid point on both axes. The third and fourth stage are analogous but consist of the remaining grid points that were skipped in the first two stages. With this order, the grid approach can observe intensities on each region of $X^{\ast }$ already after the first stage and hence can acquire more information in a faster way as with the conventional order.

We ran our approach in the default setting (Tab. 1), i. e., with automatically computed background level and intensity threshold, as well as with two alternative pairs of corresponding parameter values manually set in order to study the robustness of the results. In the default setting (row 2 in Fig. 3), the background level and intensity threshold were computed to ${\gamma }_0=58$ and ${\tau }_0=94.5$, respectively. To study the robustness of these results w.r.t. changes in the intensity threshold, the parameters for the first alternative (row 3 in Fig. 3) were set to ${\gamma }_1={\gamma }_0$ and ${\tau }_1=130$. The second alternative (row 4 in Fig. 3), for studying the robustness w.r.t. changes in the background level, was configured with ${\gamma }_2=30$ and ${\tau }_2={\tau }_0$. Eventually, the cost budget $C$ for each instance of our approach was determined by the total experimental time of the grid approach which was $\sim 11.5$ hours. Each respective initialization with the same grid of 61 measurement points took $\sim 2.9$ hours of experimental time. The results displayed in Fig. 3 show that, after initialization, our approach places measurement points mainly in regions of signal in all three settings. The grid approach, in contrast, has observed intensities at a considerable amount of uninformative measurement locations.

In scenario 2, we tried to challenge our approach with a more difficult initial setting. The initialization grid in scenario 1 (triangles in Fig. 3) indeed beneficially lies on an axis of symmetry of the intensity function. Also, partly as a consequence, some initial measurements points are already located in regions of signal. From a methodological perspective, it is interesting to investigate how our approach behaves after unfavorable initialization, i. e., if almost all initial intensity observations are lying in the background and hence almost no useful initial information is available. Hence, we reduced the number of initial measurement points from 61 to 25 and set $X=[0.2,2]\times [2,16]$ to break the axis of symmetry and measure larger energy transfers (Fig. 4a).

Estimating the background level and intensity threshold is not reasonable in this setting since the amount of initial information is too little. Therefore, we manually set $\gamma =50$ and $\tau =80$, respectively. The cost budget $C$ for our approach was not set explicitly here. In fact, we stopped the experiment manually after $\sim 7.3$ hours of total experimental time, with $\sim 1.1$ hours spent on initialization with the modified grid. The result is depicted in Fig. 4b. It shows that, even in this challenging situation, our approach keeps making reasonable decisions on measurement locations.

This section shows results of a benchmark on several two-dimensional (i. e., $n=2$) test case intensity functions (Extended Data Fig. 9) as a quantitative proof of concept. The benchmarking procedure quantifies the performance of an approach by how much benefit it is able to acquire for certain cost budgets in the context of predefined test cases. We briefly describe its setting in the following paragraph. For more details, we refer to the original work [teixeiraparente2022benchmarking].

A test case mainly includes an intensity function defined on a certain set $X$ and a synthetic TAS used for moving on $X^{\ast }$ with certain velocities and observing intensities. As mentioned, the cost measure is chosen to be the experimental time used, i. e., the sum of cumulative counting time and cumulative time for moving instrument axes. The benefit measure is defined by a relative weighted $L^2$ approximation error between the target intensity function and a linear interpolation $\hat{i}=\hat{i}\left(A\right)$ using observed intensities from experiments $A$. It encodes the fact that a TAS experimenter is more interested in regions of signal than in the background which suggests to use $i$ itself as a weighting function. However, an important constraint, which is controlled by a benchmark intensity threshold ${\tau }^{\ast }>0$ truncating weights to

as similarly described for the intensity threshold of our approach, is that separate regions of signal with different magnitudes of intensity might be equally interesting. Formally, we define

where $\parallel \cdot \parallel =\parallel \cdot {\parallel }_{L^2(X^{\ast },{\rho }_{i,{\tau }^{\ast }})}$ for the weighting function

For each test case, benefits are measured for several ascending cost budgets, called “milestone values”, to demonstrate the evolution of performance over time. We note that this synthetic setting includes neither background nor measurement noise as both are artefacts of real neutron experiments.

We run the benchmarking procedure with three approaches for comparison:

1. random approach,

2. grid approach,

3. our approach with different values for the intensity threshold parameter $\beta$.

The random approach places intensity observations uniformly at random in $X^{\ast }$. The grid approach is adopted from the section on the neutron experiment but using a square grid of dimension $P\times P$, $P\in N$.

For our approach, we set a zero background level, i. e.,  $\gamma =0$, manually since background is not included in this synthetic setting. Our approach is run with four variations of the intensity threshold parameter $\beta$ (Eq. (8)) in order to study corresponding sensitivities of the benchmark results.

The specific benchmarking procedure involves four milestone values according to the four stages (I-IV) of the grid approach, i. e., the $m$-th milestone value represents the experimental time needed to complete stage $m\in \{\text{I},\text{II},\text{III},\text{IV}\}$. Observe that they depend on the particular test case. The specific milestone values used are indicated in the Supplementary Information (Sec. S4). The number of columns/rows $P$ for the grid approach is chosen to be the maximum number such that the corresponding experimental time for performing an experiment in the described order does not exceed $9$ hours.

Since both the random and our approach contain stochastic elements, we perform a number of 100 repeated runs with different random seeds for each test case in order to see the variability of their results.

The results for each test case along with its corresponding intensity function are shown in Fig. 5.

Our approach has, on median average, performed significantly better than the random and grid approach in most test cases. Furthermore, its mostly thin and congruent shapes of variability (light color areas) demonstrate its reproducibility and its robustness w.r.t. changes in the intensity threshold parameter $\beta$. Examples of particular experiments performed by our approach are additionally provided in Extended Data Fig. 10 for each test case.

We generally intend to approximate a function of interest $f:X\to R$, which becomes the logarithm of the intensity function later. Using GPR for this, we have to assume that $f$ is a realization of a Gaussian process $F$.

Note that, for each $x\in X$, it particularly holds that

where ${\sigma }^2\left(x\right)\text{coloneqq}V\text{ar}\left(F\right(x\left)\right)={\kappa }_{\theta }(x,x)\ge 0$ is the (prior) variance function.

The kernel function ${\kappa }_{\theta }$ is an important component in GPR since it describes the correlation structure between random variables $F\left(x\right)$. Hence, it enables to include assumptions on realizations of $F$ and determines function properties like regularity, periodicity, symmetry, etc. In practice, it is crucial to choose a kernel function that matches with properties of the function of interest $f$.

We acquire knowledge on $f$ through noisy observations ${\hat{f}}_j$ at locations $x_j\in X^{\ast }$ (Eq. (2)) in our context. Therefore, for $j=1,\dots ,J$, we set

where ${\eta }_j\sim N(0,1)$ are i. i. d. random variables independent of $F\left(x_j\right)$ and $e\left(x\right)>0$ denotes the noise standard deviation at $x\in X$. Note that $\hat{F}\left(x_j\right)$ is a normally distributed random variable.

For clear notation, we define

and let

for any function $h:X\to R$.

After observations have been made, we are interested in the posterior, i. e., conditional, Gaussian process. It holds that

with the posterior mean function

and the posterior variance function

If necessary, the hyperparameters ${\theta }_J$ of the kernel function ${\kappa }_J\text{coloneqq}{\kappa }_{{\theta }_J}$ can be optimized using data ${\hat{f}}_J$. For this, we compute ${\theta }_J$ such that the logarithm of the so-called marginal likelihood, i. e.,

is maximized. A suitable optimizer is specified below. Note that the analytical expression for the integral in Eq. (22) is only feasible due to the Gaussian distributions involved. However, the computational cost of GPR is often hidden in this kernel optimization step since it requires solving linear systems and computing determinants [rasmussen2006gaussian, Sec. 2.3]. An appropriate criterion for stopping kernel optimizations that detects stagnant hyperparameters during an experiment is therefore provided below. Furthermore, observe that, for a fixed non-optimized kernel, ${\sigma }_J^2$ does not depend on observations ${\hat{f}}_J$ but only on locations $X_J$ they have been made at.

The posterior mean function $m_J:X\to R$, incorporating knowledge on $J$ noisy observations of the function of interest $f$, can now be used as an approximation to $f$, whereas the posterior variance function ${\sigma }_J^2:X\to [0,\infty )$ quantifies the corresponding uncertainties (Extended Data Fig. 11). Note that ${\sigma }_J^2\left(x\right)<{\sigma }^2\left(x\right)$ for each $x\in X$ since ${\left[{\kappa }_J\right(X_J,X_J)+\text{diag}(e\left(X_J\right){)}^2]}^{-1}$ is positive-definite. Since we have $m\equiv 0$ later, we can further simplify Eq. (20) to

The Gaussian process, which is fitted to logarithmic intensity observations in our approach, is exponentiated to the original linear scale in this section to become log-Gaussian. It is this transformation that directly yields a useful acquisition function.

Before we describe the details, we need to define a log-normal random variable and prove a technical detail that becomes important below.

The mean and variance of a log-normal random variable $Z$ are given by

Below, we look at the noise distribution of log-Gaussian processes in order to satisfy our assumption on intensity observations to contain normally distributed noise (Eq. (1)). The following result is fundamental for this and proved in the Supplementary Information (Sec. S5). It investigates the behaviour of log-normal random variables in the small variance limit.

Our result for single log-normal random variables gives pointwise convergence. However, convergence in distribution can also be proven for sums of log-normal random variables [dufresne2004log].

As mentioned, the exponentiation of a Gaussian process is log-Gaussian, by definition.

Using Eqs. (20) and (21), it immediately follows for the posterior log-Gaussian process that

In particular, its posterior mean function is

and its posterior variance function becomes

In the context of our methodology from the previous sections, we choose the function of interest to be the logarithm of the intensity function from the TAS setting (Supplementary Information, Sec. S1), i. e.,

If we define

relating to Eq. (28), i. e., $F=\log I$, Eq. (16) gives

for measurement locations $x_j\in X^{\ast }$. Note that, in contrast to the process $\hat{F}$ containing additive normally distributed noise (Eq. (16)), the noise of $\hat{I}$ is multiplicative and log-normal, i. e.,

However, referring to Eq. (1), the physical assumption on the noise of intensity observations $I^{+}$ is to be additive and normally distributed, i. e.,

where ${\eta }_j^{+}\sim N(0,1)$ and $e^{+}\left(x_j\right)=\sqrt{i\left(x_j\right)}$. Fortunately, the actual deviation of the two different noise distributions is negligible for sufficiently large intensities $i\left(x_j\right)$ as the following explanations demonstrate.

As a first step, let us determine $e\left(x_j\right)$ from Eq. (34) such that the variances of both noise distributions are equal, i. e.,

which yields

Note that the intensity term $i\left(x_j\right)$ in Eq. (38) is not known in practice but can be replaced by the corresponding observation ${\hat{i}}_j\approx i\left(x_j\right)$.

Since we aim to apply the small variance limit for log-normal random variables from above (Eq. (27)), we set

with

Observe that

Furthermore, note that the term $\mu$ also depends on the limit $i\left(x_j\right)\to \infty$ and hence on $\sigma \to 0^{+}$, in contrast to the setting above, but disappears in the calculations due to cancellations (see proof in Supplementary Information, Sec. S5). Using the small variance limit, we immediately get the following result (Extended Data Fig. 12).

If we set

the acquisition function ${\text{acq}}_J$ of our approach can now be defined as

which eventually gives Eq. (5) through Eq. (31).

As placing measurement points too close to each other has limited benefit for an efficient experiment, we mark areas around each measurement location $x_j\in X^{\ast }$ as “consumed” and ignore them as potential locations for new measurement points. Since the resolution function of a TAS yields ellipsoids in $Q$-$E$ space, it is natural to consider consumed areas in $X^{\ast }$ as ellipses. An ellipse with center point $x_j\in X^{\ast }$ is defined as

for a matrix-valued function

where $U\left(x_j\right)\in R^{n\times n}$ is a rotation matrix and with $r_k>0$. Then, the union of all ellipses at step $J$ is denoted by

It is included in the objective function ${\varphi }_J$ which is part of the final algorithm.

Incorporating all of the discussed methodological components, the steps for the final algorithm are listed in Alg. 2.