extfieldsofstudy:
Mathematics
Computer Science

text:
Coupling effect and pole assignment in trajectory regulation of multi-agent systems

This paper revisits a well studied leader-following consensus problem of linear multi-agent systems, while aiming at follower nodes' transient performance. Conventionally, when not all follower nodes have access to the leader's state information, distributed observers are designed to estimate the leader's state, and the observers are coupled via communication network. Then each follower node only needs to track its observer's state independently, without interacting with its neighbors. This paper deliberately introduces certain coupling effect among follower nodes, such that the follower nodes tend to converge to each other cooperatively on the way they converge to the leader. Moreover, by suitably designing the control law, the poles of follower nodes can be assigned as desired, and thus transient tracking performance can also be adjusted.

Introduction
Cooperative control of multi-agent systems (MASs) has been extensively investigated for the past two decades. It remains to gain increasing attention in the control community for its widespread applications in the areas of microgrids, unmanned aerial vehicles (UAVs), unmanned ground vehicles (UGVs), social networks, and so on. For a comprehensive literature review, readers are referred to some survey papers [4,6,8] and references therein.
Among various problem formulations of multi-agent systems, a fundamental one is the cooperative tracking control problem, also known as leader-following consensus problem or trajectory regulation problem. In this scenario, there is one leader node and a group of follower nodes, and all the follower nodes are driven to track the trajectory of the leader node [9,14]. A salient feature of this problem is that only part of the follower nodes can or need to have access to the leader's information due to limited communication capability, and thus a distributed control algorithm is required.
It is also worth mentioning that almost all the existing works on distributed control of multi-agent systems focus on steady state collective behaviors, without considering transient performance of the group. An exception is a line of research called prescribed performance control [2,13] or funnel control [10] of multi-agent systems. Both prescribed performance control and funnel control provide simple control laws that can restrict the profiles of the synchronization error within the prescribed error bounds. These control laws use neither the system dynamics information nor the graph topology information, and can be applied to linear systems for sure. However, their control laws consist of certain time functions, generating the prescribed performance bound profiles. This will lead to non-autonomous nonlinear closed-loop systems, whose performance depends on the initial time. Moreover, the initial values must be restricted within the prescribed bounds. Also, it is well known that both transient and steady state performance of a linear system relies on the location of its poles. Then an interesting question is naturally raised that whether we can design a distributed control law such that the poles of all follower agents can be assigned as desired. This will provides an insight into the transient and steady state performance of multi-agent systems. To the best of our knowledge, this problem is still open.
Motivated by the above-mentioned statements, this paper aims to propose a novel distributed control law, which not only solves the trajectory regulation problem, but also considers the coupling effect between agents as well as desired pole assignment. More specifically, we intend to deliberately introduce coupling effect for cooperative regulation, while avoiding over-coupling that may cause different issues, or even violation of system stability. Technically, compared with [15], this paper gives a quantitative analysis of this trade-off in terms of stability, dominant pole assignment, and fully system decoupling and hence pole assignment. It is noted that a relevant work [12] also introduces coupling effect between follower nodes, but the design relies on the solution to an linear matrix inequality and each follower node must design its own coupling effect, and pole assignment is not considered.
The rest of this paper is organized as follows. The problem is formulated in Section 2. Rigorous analysis of stability and pole assignment is provided in Section 3. A numerical example in Section 4 illustrates the effectiveness of the proposed algorithm. Section 5 concludes the paper. Technical lemmas are summarized in Appendix.
Notations: The sets of real and complex numbers are denoted by R and C, respectively. For a matrix A, A T represents it transpose and A * its conjugate transpose. For a symmetric real matrix A, A > 0 (A ≥ 0) means A is positive definite (positive semi-definite). The determinant of a matrix A is denoted as det(A). The real part of a complex number m is denoted as Re [m]. The set N is defined as {1, 2, . . . , N } and 1 N is the N -dimensional column vector whose elements are 1. The Kronecker product is denoted by the operator ⊗. For vectors x 1 , · · · , x n , col(x 1 , · · · , x n ) = [x T 1 , · · · , x T n ] T represents the stacked vector.

Problem Formulation
The paper is concerned with control of a group of linear homogenous agents of the dynamics described bẏ where x i (t) ∈ R n is the state vector of the i-th agent and u i (t) ∈ R m the control input vector. The control task is to regulate every state trajectory x i (t) to the desired trajectory x 0 (t) ∈ R n described by the leader nodė that is, lim The task becomes complicated when x 0 (t) is unaccessible for some agents. Nevertheless, such a task has been well accomplished by using the framework of consensus of observers and trajectory regulation. More specifically, an observer is established for each agent as followṡ where e = col(e 1 , · · · , e N ) with e i =x 0i − x 0 and the consensus protocol function γ i with γ i (0) = 0 is designed such that lim t→∞ e(t) = 0. Note that the design of γ i has been well studied in literature, and is not an issue considered in this paper. See the following remark for one choice of γ i .

Remark 1
For homogenous MASs, the design of consensus tracking protocol has been well studied [11,?], e.g., in a directed graph G represented by the network adjacent weightᾱ ij ≥ 0, i, j ∈ N , with µ sufficiently large [11] (or µ can be designed in a decoupled manner µ = cK as in [14]) and g i = 0 for i ∈ S R , g i = 0 for i / ∈ S R , where S R ⊂ N is the set of nodes that have direct access to the reference signal x 0 (t). In this setting, if every follower node i ∈ N is reachable in G from the leader node, then one has lim t→∞ e(t) = 0.
The time variable (t) is ignored in (4) and the sequel for conciseness. Next, let η = col(η 1 , · · · , η N ) with A trajectory regulator with a feedback matrix F ∈ R m×n and a scalar gain As the closed-loop system iṡ it suffices to pick F and d i such that A−d i BF is Hurwitz. From above, the task (3) is accomplished by combining lim t→∞ e(t) = 0 and lim t→∞ η(t) = 0. Now, the main focus of this paper is on the design of regulator for lim t→∞ η(t) = 0. It is obvious that the simple regulator (6) is designed separately for each individual agent. Its simplicity has independent interest in many scenarios. However, in some other scenarios, researchers are also interested in more sophisticated interactions among agents when they obey the trajectory regulation protocol. For instance, individual regulation does not work when the agents intend to converge to each other before to the desired trajectory. For this purpose, we deliberately introduce coupling effect a ij among agents in the physical layer (see Fig. 1(b)), which results in where a ij ≥ 0, i, j ∈ N represents the coupling weight between agent i and j in the physical layer, with a ii = 0. Note thatᾱ ij denotes the coupling among observers in the cyber layer (see Fig. 1). It is worth mentioning that the leader model (2) has the role in the design of the observer network in (5), but it has no direct influence on the deliberate coupling network in (7). In other words, the main results given in this paper, based on the analysis of (7), may apply not only to a leader-following topology, but also to more general leaderless topologies.
With the above development, this paper aims to analyze the coupling effect of A, in comparison with D. First of all, it should be noted that (i) The diagonal entries of M are represented by D and the off-diagonal entries by −A. If d i = N j=1 a ij , i ∈ N , the matrix M becomes a Laplacian. But in general, the matrix M is adjusted as a diagonally dominant matrix for our purpose with a more significant D. (ii) The graph for the adjacency matrix A is not necessarily connected for the purpose of consensus. Especially, the consensus can be achieved by the controller with A = 0, i.e., (6).
In this paper, on one hand, we intend to introduce more significant coupling effect a ij on the physical layer for cooperative regulation. On the other hand, over-coupling may cause different issues even violation of system stability. So, the main objective of this paper is to conduct more specific analysis of the coupling effect of A relative to D from the following three aspects, which apply different upper boundary conditions on a ij . These conditions can be regarded as design criteria in selecting the strength of A for (8) or (9).
(1) The condition on A relative to D for stability of (8) or (9).
(2) A less significant A and a more significant D allow dominant pole assignment of (8) or (9). (3) The condition on A relative to D for a fully decoupled system matrix and hence pole assignment.

Main Results
The three aspects listed in the aforementioned main objective are addressed in this section in order. In the subsequent presentation, we will present the results on (9).

Stability Condition
Each agent uses its own feedback through d i and others via a ij . Intuitively, it needs more significant feedback from itself for stability. To explicitly characterize the relative significance, we introduce the condition under which M + M T is positive definite. This property plays an important role in proving the stability of (9).   (10), all the Geršgorin discs are strictly shifted to the right, not going through the origin any more, illustrated by the dotted circle.
Note that the stability of the system (9) can equivalently represented by a feedback structure illustrated in Fig. 3 in frequency domain. In other words, (9) can be rewritten as follows, with v = 0, Here v is an auxiliary input introduced for the convenience of presentation later. Now, the main result is stated below. Theorem 1 Under Assumption 1 and the condition (10), there exists F such that the system (9) is stable.
Proof: A constructive proof is given as follows. First of all, Lemma 1 shows that M + M T > 0 under the condition (10). There exists a scalar k > 0 such that M + M T > 2kI. Then, Z + Z T > 0 for Z = (M − kI)/k. That is, M = k(Z + I). Now, the system (11) can be rewritten asξ Under Assumption 1, there exists a solution X C > 0 to the control algebraic Riccati equation with Q C > 0. Let F = −kB T X C . By Lemma 3, the transfer function matrix T F (s) = F (sI −A+kBF ) −1 kB is positive real, so is T F (s) = I N ⊗T F (s). In other words, It is easy to see that the transfer function from v to y is whose poles are determined by det(I + kT F (s)Z) = 0. One has det[I + kT F (s)Z] = 0, ∀ Re[s] ≥ 0, due to (12) and Z + Z T > 0, using Lemma 4. It means that all the poles of (13) have negative real parts. The stability of (9) is thus proved.
Remark 3 A uniform selection of d i , i.e., can be used to satisfy the condition (10). It simplifies the protocol (7) with less parameters. Condition (10) requires sufficiently large d i . However, high gain of a controller may lead to some practical issues. To prevent the high gain issue, in practice, we can choose d i to the desired level and then scale a ij down to match the condition (10).

Pole Assignment
In this subsection, we aim to show that the system poles can be assigned when D is sufficiently large and A sufficiently small. In particular, the transient characteristics of the regulation behavior can be adjusted by pole assignment. We will use the inverse optimal linear quadratic regulator (LQR) technique to assign poles.
Here we assume rank(B) = m and the matrices A, B take the following form where A 22 and B 2 have the same dimension. Define the LQR probleṁ with Q = Q T ≥ 0 and R = R T > 0 for the convenience of presentation. Now, the main result is stated in the following theorem.
Theorem 2 Consider the system (9) of the structure (14) under Assumption 1. Let where σ 1 , · · · , σ m > 0 and V, Σ, K are selected such that w = −F z is a solution to the LQR problem (15). In particular, the eigenvalues of A 11 −A 12 K has the specified stable eigenvalues {s 1 , · · · , s n−m }, none of which is an eigenvalue of A 11 . If the system (9) has the following eigenvalue distribution: (i) there are (n−m)N eigenvalues λ j i of the form λ j i → s j , j = 1, · · · , n − m, i ∈ N ; and (ii) there are mN eigenvalues λ j+n−m i of the form Proof: First of all, it is noted that the system (9) can be rewritten asξ For the eigenvalues of a matrix continuously depend on its parameter variation, the eigenvalues λ j i , j = 1, · · · , n, i ∈ N , approach those ofĀ −BDF as A → 0. Furthermore, the eigenvalues ofĀ −BDF are those of A − d i BF .

Remark 4
The theorem shows that the eigenvalues of the closed-loop system approach arbitrarily specified stable poles s j for j = 1, . . . , n − m and other stable poles −d i σ j when all d i are sufficiently large and all a ij are sufficiently small. In particular, as all d i are sufficiently large, Re[s j ] −d i σ j . In other words, all specified s j are the dominant poles, enforcing the plant the specified transient characteristics.

Fully Decoupling Condition
In the previous subsection, it is proved that the dominant poles of the closed-loop system can be placed as desired such that the transient characteristics can be satisfied, provided that A → 0. It obviously contradicts to the main motivation of adding coupling effect to regulation, if we simply remove the coupling by letting A = 0. In this subsection, we will further study an explicit condition on the size of a ij under which the system can be fully decoupled and the pole assignment technique can be applied without the assumption of A → 0.
We first give a lemma for diagonalization of the coupling matrix M, followed by the main theorem.
where r i = N j=1 a ij , then M is diagonalizable. Proof: The diagonal entry of the i-th row of M is d i and the sum of the absolute values of the off-diagonal entries in the i-th row is r i = N j=1 a ij . As a result, one can define a Geršgorin disc centred at d i with radius r i , denoted D(d i , r i ). Under the condition (17), the Geršgorin discs are D(d κi , r κi ) in the order of i = 1, · · · , N , from left to right, and they do not intersect (see Fig. 4). Therefore, by Geršgorin theorem, all eigenvalues of M are distinct positive real numbers. Thus M is diagonalizable. Proof: By Lemma 2, there exists a nonsingular matrix T such that P = T MT −1 , P = diag(p 1 , · · · , p N ) With ζ = (T ⊗ I m )ξ, the system (9)   In the simulation, let F = [−0.3660, 0.9306] as in Theorem 2, and pick D = diag{2, 1.95, 2.5, 2.9, 3.4}ρ and A = [εa ij ] with a 12 = a 23 = a 34 = a 45 = a 51 = 1 and a ij = 0 otherwise, for two parameters ρ, ε > 0. The results are demonstrated in terms of x i = ξ i + x 0 with x 0 (t) governed by (2) and ξ(t) by (9), where the specific behavior of e(t) approaching zero is ignored.
x 51 ] T as the vector consisting of the first elements of the five agent states. The signal φ = max{|ξ 11 |, · · · , |ξ 51 |} represents the regulation error and ψ = max{ξ 11 , · · · , ξ 51 } − min{ξ 11 , · · · , ξ 51 } the difference among the agents. In the first case, we select ρ = 0.12 and ε = 0. As shown in Fig. 5, the agents reach consensus while their transient processes are independent and do not influence each other. In other words, the agents do not demonstrate a cooperative behavior. The profiles of the agent difference and the regulation error are shown in Fig. 6.
In the second case, we select ρ = 0.2 and ε = 0.2 as com-parison. All the agent states also converge a consensus on the same reference signal. With ε = 0.2, the agents have cooperation before achieving consensus. In the simulation, we pick the parameters such that the regulation error achieves zero with the same transient performance.
As shown in the bottom graph of Fig. 6, the error reduces from 3 to 0.3 (by 90%) at t = 15 for both two cases. It is obvious that, the difference among agents in the case with ε = 0.2 is significantly less than that in the case with ε = 0. More specifically, it takes t = 4.5 for the difference among agents to reduce below 0.5 in the former case while it takes t = 7.0 in the latter case, as shown in the top graph of Fig. 6.

Conclusion
This paper has studied a consensus tracking problem of linear multi-agent systems. Compared with the conventional observer based control approach, where each observer estimates the leader's information via neighborhood communication, this paper features itself in two aspects: first, coupling effect between follower nodes in the physical layer is deliberately introduced to take into account the cooperative behavior between all follower nodes before they converge to the leader node; second, dominant poles of follower nodes can be adjusted to obtain a desired transient performance.
The following lemma is adopted from [5] with slight modification (cf. Theorem 4.1 and Proposition 4.1 of [5]). where σ 1 , · · · , σ m > 0 and V, Σ, K are selected such that w = −F z is a solution to the LQR problem. In particular, the eigenvalues of A 11 −A 12 K has the specified stable eigenvalues {s 1 , · · · , s n−m }, none of which is an eigenvalue of A 11 . If d → ∞, then matrix A − dBF has the following eigenvalue distribution: (1) there are n − m eigenvalues λ j of the form λ j → s j , j = 1, · · · , n − m; and (2) there are m eigenvalues λ j+n−m of the form λ j+n−m → −dσ j , j = 1, · · · , m.


X-MAN: Explaining multiple sources of anomalies in video

Our objective is to detect anomalies in video while also automatically explaining the reason behind the detector's response. In a practical sense, explainability is crucial for this task as the required response to an anomaly depends on its nature and severity. However, most leading methods (based on deep neural networks) are not interpretable and hide the decision making process in uninterpretable feature representations. In an effort to tackle this problem we make the following contributions: (1) we show how to build interpretable feature representations suitable for detecting anomalies with state of the art performance, (2) we propose an interpretable probabilistic anomaly detector which can describe the reason behind it's response using high level concepts, (3) we are the first to directly consider object interactions for anomaly detection and (4) we propose a new task of explaining anomalies and release a large dataset for evaluating methods on this task. Our method competes well with the state of the art on public datasets while also providing anomaly explanation based on objects and their interactions.

Introduction
All over the world there are an ever growing number of surveillance monitoring systems, in streets, in shops and our homes. In the UK, it is estimated there are 500,000 cameras in Greater London alone with an estimate of at least 1 camera per 14 people in the whole of the country [17]. Usually these systems are used for post event analysis or online event detection, in automatic settings object detectors (e.g. people, faces) are typically used. A rather more difficult task is that of detecting anomalies, whereby the scene can be complex with many moving objects and the anomaly itself can only be defined by the environmental context. We also believe anomaly detection should be interpretable, if the model truly understands the scene then it should be able to explain it's response. We address this issue here by considering interpretable object interactions in video, producing an explainable anomaly detection model and building a  Figure 1. Overview. Every human-object pair in a single input frame of a video is encoded into a high-level and interpretable feature vector (HOI vector) based on their interactions, object category and person box size. Our method can describe anomalous events in a video by comparing this feature vector to normal ones, highlighting both where in the scene and also where in the feature vector the anomaly occurs. The above example illustrates the detected objects and the detected anomalous human-object interaction: "holding a bicycle" with the most salient feature for the anomaly being the "bicycle" object. Note that bicycle holding is anomalous here as previous data suggests it is unusual to observe bicycles from this camera. dataset for evaluating the quality of explanations.
Being able to both detect and explain anomalies has an important application in many types of video monitoring settings. For example, human operators needs to be able to decide on the appropriate response, i.e. simply staying more vigilant if someone is walking in the wrong direction or intervening in person when there is a fight. This is challenging, because (a) the most alarming anomalies involve complex interactions between objects, e.g. people pushing each other next to a road, (b) interpreting the anomaly re-  Figure 2. HOI vector. All pairs of detected objects (from an R-CNN network) are encoded into (1) a spatial and semantic encoding using a graph convolutional network (GCN), (2) human appearance feature and (3) object appearance feature. Two classifiers ingests these features to produce interaction probabilities and object class scores. These probabilities, together with bounding box size, are concatenated to form the HOI vector.
quires high level understanding of the video context, e.g. is anomalous motion in the scene due to a person running, or a tree blowing in the wind? and (c) the model also has to explain it's response.
In an effort to address all three of these issues, under the task of explainable anomaly detection, we make the following contributions: (1) we show how to build interpretable feature representations suitable for detecting anomalies with state of the art performance, (2) we propose an interpretable probabilistic anomaly detector which can describe the reason behind it's response using high level concepts, (3) we are the first to directly consider object interactions for anomaly detection and (4) we propose a new task of explaining anomalies and release an accompanying dataset, X-MAN (eXplanations of Multiple sources of ANomalies), the first dataset for evaluating explanations in anomaly detection. Our system (illustrated in Fig. 1) is capable of (i) outputting a high level description of the anomalous event e.g. holding a bicycle, (ii) recovering the most similar normal event e.g. holding a bag and (iii) is also able to identify the most salient features (from a deep neural network) causing the anomaly, e.g. the held bicycle object.

Anomaly detection
A common approach to anomaly detection with deep learning is to predict future 'normal' video frames from sparse feature representations [3,15,26], sometimes augmented with memory modules [19], and/or optical-flow images [13,20,21]. Anomalies are detected based on the assumption these models will find it difficult to generate abnormal frames. An alternative approach is using or learning these feature representations with out of distribution detection algorithms [4,22,6], as can be done for anomaly detection [18,24,27] or more recently and successfully, applying a one-vs-all cluster classification [7,8]. However, these methods do not explain the high level content of the scene, the 'typical' normal events or explain why an event is abnormal in an interpretable way. Our work is novel in that it uses traditional probabilistic modelling but leverages a deep network's ability to extract high level information about the scene to allow an event to be (1) interpretable and (2) explainable based on a saliency mapping over a learnt feature representation.

Anomaly explanation
To the best of our knowledge [5] is the only work tackling the task of describing the detected anomalies. Here, a sparse feature representation of the scene is simultaneously used as input to an anomaly classifier and object, action and attribute detector. The shortcoming of this approach is neglecting the interaction between objects. In many applications the interactions between objects are the crucial source of anomaly, for example it may not be abnormal for a person to jump, but it is abnormal for a person to jump over the gates of the underground. Moreover, this approach is still not fully explainable, i.e. input features for detecting anomalies cannot be inspected. Our method uses higher level reasoning than [5] (both high level object categories and human-object interactions) while also justifying the algorithms decision by returning interpretable salient features. Due to the lack of labelled data, there are also difficulties in forming a benchmark on the task of anomaly explanation. To this end, we build such a dataset and show how to make improvements in explaining anomalies by using human-object interactions.

Method
In Sec. 3.1 we outline how we encode the frames of a video into objects and human-object interactions. In Sec. 3.2 we detail the anomaly detection method, and in Sec. 3.3 how we interpret and explain the decisions from the anomaly detector.  The estimated means are interpreted as normal events, later used for explanation. Bottom: Test frame detected as anomalous, showing at least one HOI vector has low probability under the GMM. High-level information in the vector is inspected to interpret the anomaly. Saliency mapping highlights the most salient features by comparison to the closest normal HOI vector.

Encoding the video frame
Every frame of a video is encoded into a set of humanobject interaction vectors (HOI vectors), see Fig. 2. HOI vectors are interpretable as they are formed from high level interaction and object features. For every human-object pair in the scene a HOI vector is recovered and consists of object class score (one-hot encoding with object class confidence from the object detector in place of the "1"), the probabilities of each of the interaction classes and person bounding box dimensions. We use Faster R-CNN [23] for object detection, followed by Dual Relation Graph [1] for Human-Object Interaction recognition.

Anomaly detection
We formalize the anomaly detection as a probabilistic outlier detection problem. The scene is represented by a statistical model of the HOI vectors. An outlier/anomaly is declared when the probability of a HOI vector under this model is below a threshold.
All HOI vectors from normal (not anomalous) frames are first reduced in dimensionality using PCA. We train a Gaussian Mixture Model (GMM) from these representations of 'normal' as in Fig. 3(top). This model allows for a probabilistic modelling of the environment. During inference we declare a video frame anomalous if any of it's HOI vectors falls below a threshold probabilty under the GMM Fig. 3(bottom). Using our proposed saliency mapping, an interpretable explanation for the model can also be derived, this is detailed further below.
Temporal information In order to use make use of the temporal information provided by a video we smooth the per-frame anomaly scores with a Gaussian filter.

Anomaly explanation
Our method allows for explanation of abnormalities on three levels: an explanation of the tested event, the closest normal event and the salient features responsible for the deviation. We detail the methods used to explain each of these below: Tested event. In our framework, the vector describing an event consists of high-level information: object class, interaction class probabilities and bounding box size. Any event can be explained by a direct inspection of the corresponding HOI vector v, as shown in Fig.3. PCA is appled to reduce the dimensionality of the HOI vector: Normal event. We fit a GMM to the distribution of all PCA reduced HOI vectors from all normal video frames in the training set using M mixture components. The set of means for modes m ∈ {1, . . . , M }, denoted {µ m } represents our set of normal events.
Deviation from normal. The deviation ∆ m of an anomaly x (PCA-reduced HOI vector) from the closest normal event µ m can be measured as a probability of being generated by mixture component m ∈ {1, . . . , M }: Where Σ m is the covariance matrix of k and the quadratic term inside the exponential is a weighted distance of the anomaly from the closest normal event.
Saliency mapping. Expressing the exponential component of equation 1 using the original HOI vectors v: allows us to decompose the contribution of the summation (deviation from normal) across each element of the HOI vector, where i and j index the element. Each term in the summation consists of a difference between a pair of elements from the original HOI vector and normal event HOI vector y = W −1 µ. These distances are weighted by (W Σ −1 k W T ) ij and visualised as a 2D heatmap. 'Hot spots' being the feature with most contribution to the anomaly. The 2D location of the hot spots can explain which part of the HOI vector is responsible for the anomaly detection, e.g. if the hot spots lie in the upper left corner then the interaction class is responsible, yet if the hot spots spread along from the diagonal in the middle part of the vector, then the anomaly is due to an unexpected object , see Fig.4 for examples. A 3x3 heatmap can also be formed which better visualises the contribution from object, interaction and bounding box size by taking the sum within each component, see Fig.4.
Quantitatively, we use the values on the diagonal of the heatmap as scores to indicate the classes responsible for the HOI being flagged as anomalous. To obtain per-frame explanation scores, we keep the maximum score of each class across all HOI vectors in a frame. The explanation score for location being anomalous is taken by summing the scores corresponding to the bounding box dimensions. The sum of scores across all classes is not normalised, hence allowing for multiple correct explanations (a frame can contain more than anomalous event).

X-MAN dataset
We collect a new dataset for evaluating anomaly explanations: X-MAN (eXplanations of Multiple sources of ANomalies). The dataset supplements existing public anomaly detection datasets and provides labels describing the reason behind an anomalous event. Together, the existing public anomaly detection datasets and X-MAN, can be used as a benchmark for anomalous event localisation, i.e. detecting and categorising anomalies in video, therefore providing a measure of joint anomaly detection and explanation.

Labelling anomaly explanations
The labelling task was conducted by human operators using a custom-built labelling tool (see Fig. 5). The tool was used to label each anomalous frame in existing public datasets with explanations describing each of the anomalous events occurring in the frame on both a coarse level (e.g. action) and a fine-grained level (e.g. running). Anomalous events are seen only in testing, hence this dataset can be used for evaluation only.
Coarse labels. There are 3 possible classes of coarse labels: "object", "action" and "location". They capture the high-level reason behind an event being categorised as anomalous. Only these 3 classes of coarse labels are considered, because they cover all anomalies in the public anomaly datasets and any obvious real-life examples (e.g. guns or fighting).
Fine-grained labels. Fine-grained labels explain the anomalous events in detail, supplementing the coarse labels in a hierarchical manner. In the case of a coarse "object" label, the fine-grained label would specify which object made the frame anomalous, e.g. bicycle. For the coarse "action" label, fine-grained labels specify both the activity and the interaction object if one exists: examples include "run", "ride,motorcycle", "fight/hit,person". In case of actions with more than one object involved, e.g. person hitting another person with a bag, the interaction is decomposed into more than one interaction, each involving one object, e.g. "hit,person", "hit,bag (with)". Finally, for the coarse "location" label, fine-grained labels specify what object is in an unusual location, e.g. a human walking in a prohibited space or an unattended bag on the ground.
Fine-grained object and action labels use categories defined in MS COCO [12] and AVA Kinetics [10], respectively.
Both of these datasets capture a wide range of classes. However, it is expected that some anomaly categories will not be included, because anomalies are associated with rare events, hence they are unlikely to all be included in standard datasets.

Dataset statistics and metrics
Statistics. Our newly created dataset consists of 22,722 manually labelled frames in ShanghaiTech (17,362), Avenue (3,712) and UCSD Ped2 (1,648). Each frame contains between 1 and 5 explanation labels, each label being a different reason why the frame is anomalous (many  First 29 elements of the vector correspond to the interaction class, the following 81 to the object category, the last 2 to the box size. Red indicates which pairs of elements contributed the most to the weighted distance from the normal event. We also summarise the results in a 3x3 heatmap (right column) show the salient category. Examples of cases when the interaction class was anomalous (top) and when the object was anomalous (bottom).  Figure 5. Labelling tool. GUI displays the anomalous frames in the video and provides tools to label the explanations. The GUI has 5 distinct areas: (1) highlights the abnormal region frame, (2) allows for switching frames and reapplying labels from the previous frame, (3) provides interface to add anomalous object and action categories which are show in (4), and finally (5) creates corresponding buttons that are used to label the frame with explanations.
frames contain multiple anomalous events, e.g. one person running and one riding a bike). In total, there are 40,618 labels across all frames. The majority of anomalies (22,640) are due to actions, followed by anomalous objects (14,828). The remaining anomalies are due to an anomalous location. There are 42 anomalous actions (Fig. 6) and 13 anomalous objects. While the majority of actions can be found in AVA Kinetics, training videos with many action-object combinations in Fig. 6 would not be found (because they are anomalous), and hence rare, making the X-MAN dataset challenging and complementary to existing ones. Evaluation metric. The task of explaining anomalies is similar to that of action recognition: in each anomalous frame the system has to recognise the anomalous event, there are multiple possible events in each frame and the recognition has to be accurate temporally. Hence, similar to action recognition, we propose to use mean Average Precision (mAP) as the evaluation metric for the anomaly explanation task. The mean is taken across the different explanation classes in order to weight rare explanation classes equally to the common ones.

Experiments
We evaluate our approach for both detecting anomalies and explaining them. For anomaly explanation we utilise our new dataset (described above), for anomaly detection we compare to state of the art on existing public datasets described below:

UCSD.
A standard benchmark for anomaly detection. The training data contains only normal events, while testing data contains some abnormal events. 19600 frames captured using two different cameras: UCSD Ped1 and UCSD Ped2 which contains 16 training and 12 testing videos. Normal events include pedestrians walking, while abnormal events include trucks, cyclists and skateboarders. Following [5] we evaluate on Ped2 only as Ped1 is very low resolution.
Avenue. This dataset contains contains 16 training and 21 testing videos. All captured from the same scene, a total of 30,652 (15,328 training, 15,324 testing) frames. This is a challenging dataset because it includes a variety of events such as "running", "throwing bag", "pushing bike" and "wrong direction". We train from the videos in Avenue that contain normal events. This dataset focuses on dynamic events e.g. walking in an uncommon area in the scene and regards abnormal static events as normal e.g. standing in the same uncommon area.
Avenue17. Was originally proposed by [5] to better evaluate methods which can also detect static abnormal events. It consists of a subset of 17 videos from Avenue which exclude abnormal static events.
ShanghaiTech. The largest and the most complex anomaly detection dataset. Contains 274,515 normal training frames in 330 videos and 42,883 testing frames in 107 videos with 130 complex abnormal events, e.g. fighting, jumping over a fence. Videos are from 13 scenes with varying lighting and camera angles.
Evaluation metric. All test video frames from all datasets are marked as either containing or not containing an anomaly. Measuring the true and false positive rates against this ground truth, we use the standard metric of evaluating abnormal event detection: the area under the ROC curve (AUC).

Implementation details
We use an implementation of the object detector and DRG network from the authors of [1], and the provided pre-trained models from COCO [11] and V-COCO [2] datasets. The networks are trained to detect 29 interactions and 81 different objects including humans, humanhuman interactions are allowed. Each HOI vector is therefore 29 + 81 + 2 = 112 (the bounding box width and height of the human detection are included in the vector, because they carry information about proximity of the person to the camera, hence implicitly providing information about location). We apply PCA dimensionality reduction. The output size varies between 19 and 23, depending on the dataset, and is chosen to capture over 99 % of variance. The number of mixture components M for GMM fitting is chosen based on the Bayesian information criterion [25] (BIC) and estimating the 'elbow' of the BIC vs M curve. This varies between 2 and 7 depending upon the dataset.

Results
Quantitative anomaly detection state-of-the-art (SOTA) comparison. Table 1 summarises the AUC on all datasets. Our method performs competitively against all other methods. We believe our scores are slightly lower than the SOTA because it uses simpler clustering methods. However, it is the only method providing full explanations of the anomaly detection decisions. In a practical setting this is advantageous, allowing human operators of monitor systems to decide on an appropriate responses. Alternatively, it might allow for grading of the alert level raised by the anomaly, i.e. anomaly due to a gun is more alarming than an anomaly due to a person jumping.
Method UCSD Ped2 Avenue Avenue17 ShanghaiTech Kim et al. [9] 59.0 --Mahadevan et al. [16] 82. 9 --Lu et al. [14] -80.6 80.3 Hasan et al. [3] 90.0 70.2 76.9 60.9 Luo et al. [15] 92.2 81.7 -68.0 Liu et al. [13] 95.4 85.1 -72.8 Park et al. [19] 97.0 88.5 -70.5 Hinami et al. [5] 90.   Qualitative anomaly detection and explanation analysis. Fig.7. shows examples of correctly detected anomalies from our method all 3 datasets. The interpreted anomaly is shown as an object and interaction pair (in orange and blue respectively). Using our saliency mapping method, we show in the red box the high level feature causing the event to be detected as anomalous. Examples of detected anomalies due to unexpected objects, interactions and box sizes are demonstrated. In Fig.8. examples of false positive failures of our anomaly detector are shown. Note how using high level features and the saliency mapping we can explain these failure cases quite easily i.e. incorrect object detection. Please also see our supplementary video showing the X-MAN dataset and anomaly detections and explanations from our system.
Quantitative anomaly explanation analysis. Tab. 2 shows the mAP achieved by our system on the anomaly explanation task. We evaluate our method against a baseline, where we remove the interaction part of our HOI vector, keeping only the object classes and bounding box size. The dataset contains classes that follow mostly the COCO and AVA datasets, but HOI methods, including DRG (which we use here), are trained on the V-COCO dataset, hence there is a mismatch of class labels. To mitigate this effect, we evaluate our method on the subset of classes in the X-MAN dataset that overlap with the classes in V-COCO and COCO datasets. In Tab. 2 it can be seen that across all 3 datasets using interactions improves the explanation performance of our system, confirming our intuition that interactions are a crucial source of anomalies.

Conclusion
This paper addresses the problem of explainable anomaly detection. For this task we release a large new dataset for evaluation and propose an interpretable probabilistic anomaly detection framework which can reason about high level video content. Our model automatically generates three levels of explanation for it's response: explaining the anomalous event, describing the closest normal event and highlighting the salient features. The qualitative analysis shows the potential of our explainability module and how it is able to indicate which of the interpretable features led to an anomaly. We also believe this is the first method for detecting anomalies using object interactions. Our model performs competitively against the state-of-theart on 3 standard datasets. Quantitative evaluations on our newly gathered dataset show that using object interactions improves the anomaly explanations. Future work will focus on improving the interaction-object classifier on combinations not seen during the training phase, we believe this would help generalise the method further and increase performance of the anomaly detector.


extfieldsofstudy:
Computer Science
Medicine

text:
Assigning Confidence to Molecular Property Prediction

Introduction: Computational modeling has rapidly advanced over the last decades, especially to predict molecular properties for chemistry, material science and drug design. Recently, machine learning techniques have emerged as a powerful and cost-effective strategy to learn from existing datasets and perform predictions on unseen molecules. Accordingly, the explosive rise of data-driven techniques raises an important question: What confidence can be assigned to molecular property predictions and what techniques can be used for that purpose? Areas covered: In this work, we discuss popular strategies for predicting molecular properties relevant to drug design, their corresponding uncertainty sources and methods to quantify uncertainty and confidence. First, our considerations for assessing confidence begin with dataset bias and size, data-driven property prediction and feature design. Next, we discuss property simulation via molecular docking, and free-energy simulations of binding affinity in detail. Lastly, we investigate how these uncertainties propagate to generative models, as they are usually coupled with property predictors. Expert opinion: Computational techniques are paramount to reduce the prohibitive cost and timing of brute-force experimentation when exploring the enormous chemical space. We believe that assessing uncertainty in property prediction models is essential whenever closed-loop drug design campaigns relying on high-throughput virtual screening are deployed. Accordingly, considering sources of uncertainty leads to better-informed experimental validations, more reliable predictions and to more realistic expectations of the entire workflow. Overall, this increases confidence in the predictions and designs and, ultimately, accelerates drug design.

I. INTRODUCTION
While data-driven modeling has made large advances in image processing and speech recognition, groundbreaking contributions in drug discovery are harder to find. Often, only limited amounts of reliable in vitro and in vivo data suitable for supervised learning tasks are available. Additionally, uncertainties of experimental data, for instance regarding effectiveness and toxicity, are common and significant as the experimental design is complex, prone to noise and the outcome is dependent on many parameters, in particular compound dosage. Hence, these uncertainties are propagated to property prediction workflows when data-driven models are trained on that data.
While it is evident that machine learning (ML) prediction accuracy and data quality are intrinsically related, there are more subtle aspects that are sometimes neglected.
Dataset size, composition, coverage and training-test split have considerable influence on the final model performance but are insufficient to assess uncertainty alone without suitable models and methodologies for that purpose. For instance, using random or scaffold-based train/test splits does not lead to reliable measures of model performance and generalizability. Other common problems include high dataset bias, small sample size or low chemical diversity. However, some of these issues currently do not have a straightforward remedy as there is no commonly accepted way to describe chemical subspaces rendering it difficult to uncover all inherent biases.
Nevertheless, when the underlying data is not investigated appropriately model performance and generalizability tend to be overestimated leading to a significant number of false predictions and, ultimately, low confidence in the property prediction workflows.
In this article, we discuss methodologies and datasets for chemical properties important in drug design. In particular, we focus on the main sources of uncertainties for both simulation-based and data-driven property prediction. In doing so, we discuss uncertainties inherent in datasets, outputs of data-driven models, input features, and simulation of binding affinities. In addition, we dis-cuss the importance of uncertainty in generative models, especially for property-based molecule design. Finally, we close with our personal opinion on the most important aspects of uncertainty and confidence focusing on important problems and future avenues that will lead to higher predictive ability and models that naturally take uncertainty into account.

II. DATASET UNCERTAINTY
One of the most common challenges for data-driven molecular property prediction is dataset size and composition. For many important properties in drug discovery, especially pharmacologic properties of molecules such as absorption, distribution, metabolism, excretion and toxicity (ADMET), only a limited amount of high-quality data is available and it is usually only available for certain classes of molecules, which inherently introduces biases. To provide an overview of the amount and type of data available, Table 1 lists popular datasets for molecular property prediction. Accordingly, in this section, we discuss dataset characteristics that are important to consider for data-driven property prediction tasks along with methods and procedures to minimize their impact on property prediction performance.
In the framework of supervised learning, models learn relationships between input and output from the training data and, subsequently, predict properties for unlabelled input data, and the corresponding prediction ability is termed generalization. When predictors perform well on the training set, but performance on the unlabelled data is poor, the model is overfitted. This is a very common problem in ML, especially when the training set size is small (Figure 1-A). Learning theory states that collecting more data improves generalization if the data is independent and identically distributed [20]. However, real data is not a uniform sample of chemical space. Typically, molecules in a dataset are collected under specific criteria such as the number of atoms, the constituent elements, the similarity to known molecules or the availability of synthetic procedures introducing bias (Figure 1-B). In addition, there is an inherent bias in both industry and academia to publish only successful experiments which can be problematic for assembling training data as negative results are as important as positive ones for data-driven modeling.
First, the uncertainty of model predictions for unlabelled molecules is largely dependent on the difference of the structure and property distributions to the training data. To quantify this difference, in the field of quantitative structure-activity relationships (QSARs) [21], the concept of applicability domain (AD) is widely applied. The AD of a QSAR model is defined as "the response and chemical structure space in which the model makes predictions with a given reliability" [22] and makes the inherent bias of a model comprehensible (Figure 1-C). There are several methods to estimate ADs [23]. A common approach is to consider molecules with descriptors within a certain distance to the mean of the training data to be inside the AD. Notably, predictions for molecules inside the AD are considered reliable, outside the AD reliability cannot be guaranteed.
Secondly, a significant bias in the training data can cause models only to learn the inherent bias rather than physically meaningful relationships (Figure 1-B). When users are unaware of this bias, it inadvertently leads to overly optimistic conclusions concerning model performance. Recently [10], it was pointed out that there is significant hidden bias in the widely used dataset for structure-based virtual screening, the Directory of Useful Decoys: Enhanced (DUD-E) [24]. DUD-E includes two types of molecules: ligands, small molecules that bind to a receptor, and decoys, molecules with similar physical properties but dissimilar chemical structures compared to the ligands. Using DUD-E as a dataset for training, the performance of receptor-ligand and ligand-only ML models were compared to determine whether the receptor-ligand models also learned from the actual interactions or whether they just learned the inherent ligand bias. Ideally, the performance of these two types of models should be different as the binding energies are a function of both ligand and receptor structures, and the receptor structure can generally not be inferred from the ligand alone. However, the performance was found to be equivalent, the receptor-ligand models did not learn from the receptor structure, and it was concluded that "models may only learn the inherent bias in the dataset rather than physically meaningful features" [10]. This calls for general methods to test what property predictors actually learn as it would naturally lead to more powerful, generalizable models and to more confidence in the models. To minimize hidden biases, the importance of dataset diversity for the generalizability of ML models has been pointed out before [25].
Finally, small datasets pose a particular challenge as they are inherently biased due to their limited size and, additionally, overfitting is to be expected. An important strategy to reduce both bias and overfitting in ML models when data is sparse is transfer learning (Figure 1-D) [23,26]. When using that approach, models are first pre-trained on large, sometimes unlabelled, datasets, and subsequently trained on the property of interest. Importantly, transfer learning has been shown to improve the prediction performance for models trained on small datasets significantly [27,28]. For instance, recently [27], a language-based ML property prediction model was pretrained first via self-supervised learning on one million unlabeled molecules curated from ChEMBL. This pretraining improved the prediction performance after subsequent training on small datasets of molecular activity consistently, regardless of the specific property prediction FIG. 1. Common problems with datasets for data-driven molecular property prediction and strategies to improve confidence in the models. A) Small training datasets hamper generalization. B) Biased training datasets lead to models only learning the inherent bias rather than meaningful relationships. C) The applicability domain of a model defines the molecules for which the property predictions are expected to be reliable. D) Transfer learning via pre-training on large (unlabelled) datasets improves prediction accuracy for property prediction based on small training datasets.
task and the training dataset size. Notably, the smallest dataset only contained 642 samples [27].

III. UNCERTAINTY IN THE OUTPUTS
Supervised ML largely relies on an inductive inference approach to derive general rules from a finite number of training examples. Therefore, predictive models obtained using inductive inference are never formally correct and inherently uncertain.
There are often additional sources of uncertainty that stem from noise or imprecision on target measurements. For instance, molecular properties like binding affinity are often measured with an uncertainty, e.g. 1 µM ± 10%. The latter type is referred to as aleatoric uncertainty, while the former is referred to as epistemic or systematic uncertainty and can in theory be reduced in light of more training observations or increased knowledge about the correct model. Both types of uncertainty are important, although they are often not treated separately and the formal distinction between the two remains an active area of research [29]. Many ML methods are available for representing uncertainty. They can be categorized with respect to the way uncertainty is represented and whether or not they allow for differentiation between the two types (aleatoric and epistemic). Here we divide the discussion of methods derived from frequentist statistics and ones based on Bayesian inference.
A widely used class of frequentist methods are ensemble strategies which typically use a large set of distinct predictors as opposed to a single model. Ensembles of tree-based models such as random forests [30] are commonly used to produce frequentist probability estimates based on the relative output frequency of the ensemble members [31,32]. Ensemble learning can introduce variance through multiple ensemble models with different parameter initializations (ensembling) or through the use of randomly sampled training sets (bagging or boosting, Figure 2-a). Ensembles of neural networks can also provide computationally efficient, Database of commercially available compounds together with very simple estimated molecular properties for virtual screening.

billion
Inherently biased by currently synthesizable chemical space. Consequently, the molecular shapes have been shown to be highly biased against sphere-like molecules.

thousand
Biased towards small molecules only containing the elements C, H, N, O and F.
PubChemQC [3,4] Geometries and electronic properties of molecules with short string representations taken from PubChem.

million
Biased towards small molecules that have been reported in the literature before.
Tox21 [5] Toxicologic properties of molecules with respect to 12 different assays 13 thousand Biased towards environmental compounds and approved drugs.

ToxCast [6]
High-throughput screening and computational data for the toxicology of molecules from industry, consumer products and the food industry based on cell assays.

thousand
Biased towards molecules used in industry, consumer products and the food industry.
ClinTox [7] Drugs and drug candidates that made it to clinical trials and were either approved or failed.

thousand
Biased towards drugs that made it to clinical trials.
SIDER [8] Recorded adverse drug reactions of marketed drugs. 1.4 thousand Biased towards marketed drugs.
ChEMBL [9] Bioactive small molecules and their activities extracted from the literature, from clinical trials and from other databases.

million
Biased towards compounds for which bioactivity was published in the scientific literature.

DUD-E [10]
Ligand binding affinities against 102 distinct target proteins with both strong and weak binders.

thousand
Biased towards molecules that have been synthesized and evaluated for binding affinity.
AqSolDB [11] Aqueous solubility data of organic molecules taken from 9 different datasets. 10 thousand Biased towards organic molecules with relatively high aqueous solubility.
Olfaction Prediction Challenge [11] Olfactory perception of organic molecules at different concentrations.

thousand
Biased towards small and volatile organic molecules. Results biased by familiarity of smells.
FreeSolv [12] Experimental and computed hydration free energies of small and neutral molecules. 0.6 thousand Biased towards small and neutral molecules that have been studied in the literature both computationally and experimentally for hydration free energies.
ESOL [13] Experimental aqueous solubility combining datasets for small molecules from the literature, for medium-sized molecules used as pesticides and larger proprietary compounds from the pharmaceutical industry.

thousand
The sub-groups each have a different bias as they each have different application domains.

thousand
Biased towards molecules with distribution coefficients between -10 and 10.

million
Biased towards molecules of interest and molecules that are synthesizable.

thousand
Biased towards complexes with available crystal structures.
BBBP [19] The blood-brain penetration partition coefficient for molecules collected from the literature.

thousand
Biased towards molecules studied in the literature for blood-brain penetration.
readily parallelizable uncertainty estimates [33,34]. However, care must be taken using these methods as the uncertainties are often biased, poorly calibrated, and are only concerned with aleatoric uncertainty. Heteroscedastic aleatoric uncertainty can also be modeled using a single neural network trained with frequentist maximum likelihood inference that outputs a probability distribution ( Figure 2-b) [35]. Conformal prediction (CP) has recently received significant attention, in which one ML model is trained to predict a property and another to predict the uncertainty (Figure 2 d) [36][37][38][39]. CP methods are based on a rigorous mathematical framework and allow users to select intuitive confidence levels for their predictions, e.g. selecting a confidence level of 0.9 means that at most 10% of predictions will be outside the predicted range.
Methods inspired by Bayesian inference seek to update a prior distribution p(θ) on the space of all possible models in light of training data, D, yielding a posterior distribution p(θ|D). Given a molecular representation for prediction, x * , the predictive posterior distribution p(y * |x * , D) is in theory obtained by averaging all possible models weighted by their respective posteriors. Bayesian methods have experienced a recent resurgence in drug design applications due to increasing computational Gaussian processes (GPs) [40] allow for Bayesian inference in a non-parametric way, where the prior is determined by a mean and kernel function. The prior may be updated with observed data and for regression with Gaussian noise, the predictive posterior distribution is Gaussian with mean µ and variance σ 2 . The predictive variance σ 2 represents the total uncertainty but can be decomposed into the variance of the error term σ 2 which represents aleatoric uncertainty. Epistemic uncertainty is the difference between σ 2 − σ 2 and is determined by the hyperparameters of the kernel function such as the characteristic length scale. GPs are commonly used tools for property prediction in the low data regime (i.e. below 1000 data points) due to their known robustness to overfitting [40][41][42]. Hie et al. [43] used a GP trained on less than 100 compounds to screen a large library for candidates with nanomolar affinity for diverse kinases and whole-cell growth inhibition of Mycobacterium tuberculosis [43]. The authors show how uncertainty estimation can enable machine-guided discovery. However, the traditional implementation of GPs incurs cubic scaling with training set size, rendering them less desirable when training data is abundant.
Bayesian neural networks (BNNs) place probability distributions on neural network parameters w and the analytically intractable posterior distribution p(w|D) is typically estimated using variational Bayesian inference [44]. The predictive posterior distribution is approximated using Monte Carlo sampling, which intuitively involves averaging the predictions of multiple parameter configurations sampled from the posterior (Figure 2-c).
Replacing point-estimated network parameters with distributions allows for the quantification of epistemic uncertainty as well as heteroscedastic aleatoric uncertainty. Recently, such networks have been combined with state-of-the-art representation learning to give uncertain predictions of the physiological properties of small molecules [45].
For classification tasks, confidence calibration seeks to produce probability estimates which correctly reproduce the true correctness likelihood. In other words, given N predictions each with confidence 0.9 about a label, we would expect 90% of the predictions to give the correct label. Scalar valued confidence statistics such as expected calibration error and maximum calibration error have also been proposed. Many calibration methods such as isotonic regression or temperature scaling are based on a learned post-processing step in which a calibrated probability q i is derived from the model's output probability p i * [46]. Confidence calibration schemes can be applied as a scalable post-processing step to any predictive learning algorithm, and their ideas can be extended to regression tasks.
Confidence or uncertainty estimation in molecular property prediction is crucial for tasks that involve data-driven decision making where valuable resources or human risk is at stake. For example, BNNs have been employed to give intuitive probabilistic statements about the severity of the risk of drug-induced liver injury [47,48]. Predicting clinical outcomes from pre-clinical data and using uncertainty estimates to quantify and communicate risks can preserve time, resources, and patient well-being [49]. Predictive uncertainty is central to the field of active learning [50], which has been known in the context of drug discovery for almost two decades. Active learning strategies balance exploitation of the knowledge of a predictive model to retrieve active candidate compounds with the exploration of regions in compound space where the model is highly speculative [45,51,52]. Explorative behaviour of an active learner is paramount when structural novelty is prioritized. Elucidation of error sources in training data remains a critical application for uncertain predictive models. Ryu et al. modeled aleatoric and epistemic uncertainty separately using Bayesian graph convolutional neural networks for a variety of molecular property prediction tasks including classification of bioactivity and toxicity. Analysis of predictive aleatoric uncertainty allowed for the identification of poor quality training data for some molecular candidates [53].
In the near future, we believe predictive uncertainty will become an indispensable component of molecular property prediction.
Researchers will continue to systematically benchmark uncertainty quantification methods on various datasets [54]. We anticipate software packages that facilitate uncertain predictions will continue to appear [55].

IV. UNCERTAINTY IN THE INPUT FEATURES
Contrary to uncertainty in the properties being predicted, imprecision in the input features is a source of uncertainty that is rarely considered. Property prediction in drug design often relies on representations encoding a specific molecule, such as SMILES strings [58], or on three-dimensional structural data. However, these static descriptions of molecular systems are not representative of their dynamic nature in solution. In this section, we discuss how input-related uncertainty might arise when predicting chemical properties from molecular structures, and how some of the latest ML models might be used to capture this uncertainty in the confidence of the predictions.
One form of input uncertainty is due to the presence of multiple protonation and tautomeric states in solution. String and graph representations of molecules are static and unable to capture these multiple states and their probabilities in solution. SMILES strings representing two different tautomers are likely to return different property predictions when they in fact should be the same given the behaviour of these molecules in solution. This represents a challenge to many compound property predictions and cheminformatics tools and may be seen as a form of input uncertainty, given that the chemical species in solution can be modeled only approximately. Similarly, protonation states depend on the chemical environment of the molecule, there may not be one dominant state in solution and, finally, these may change upon interaction with a binding partner. All these aspects contribute to creating uncertainty about the input representation used in ML prediction tasks.
Three-dimensional structural data is often used to predict properties like drug binding affinity, selectivity, and mutational effects on protein stability and drug resistance. This data is typically derived by X-ray crystallography and, increasingly, cryo-electron microscopy (cryo-EM). Despite providing abundant information on protein-ligand interactions, which can be used to drive drug design, these structures are only a single snapshot of all possible protein-ligand conformations in solution. The protonated form of histidine has a pKa of approximately 6.0, which means that it is generally expected to be found in its neutral form at physiological pH. Yet, because of the local environment in protein binding pockets, it is often unclear which of the two is the major species, and they may both co-exist in non-negligible proportions. On the right-hand side, some of the tautomeric states of warfarin [57] are shown.
As such, this three-dimensional data is a single datapoint of a much broader input probability distribution. Some structure-based modeling approaches, like ensemble docking, try to take this uncertainty into account by sampling multiple conformations with molecular dynamics simulations. In such a way, ensemble docking tries to capture the distribution/uncertainty of protein-ligand structures. By generating an "ensemble" of drug target confirmations, the docking results are less likely to be overfitted to a specific protein conformation. However, fast scoring functions for ligand binding affinities are typically trained on databases of static structures, such as the PDBbind database, and do not capture this type of uncertainty in their predictions. In addition to the statistical uncertainty caused by conformational variability, differences between protein structures in crystallo and in vivo due to packing effects and cryo temperatures may result in systematic biases.
Ideally, one would like uncertainty about multiple possible chemical and conformational states in solution to be captured by the ML model and reflected in the confidence of the property predictions. Currently, only a handful of specialized ML approaches are able to capture uncertainty in the input features. Gaussian process latent variable models, together with other similar kernelbased approaches [59,60], are examples of such models [61,62]. Recent work has also started to investigate the effect of molecular ensembles in property predictions [63] with ML models. However, we are unaware of widely adopted supervised ML models that can either take into account known input uncertainty, or that can infer it directly from data. We believe this to be an underexplored research direction with much potential for practical applications. Such models could become an indispensable component of molecular property prediction, being able to capture or infer feature uncertainty in structural and chemical data.

V. UNCERTAINTY IN THE BINDING AFFINITIES
A key property of interest in drug discovery is the binding affinity of a ligand to a target receptor. Of the many physics-based models that have been developed and refined to predict binding affinity, here we focus on two widely-used techniques that utilize the three-dimensional structures of ligand and receptor, and a molecular-mechanics energy function to compute ligand affinities: molecular docking and free energy perturbation (FEP) methods.
Large uncertainties in docking estimates can be attributed to conceptual deficiencies in the scoring function, for example, a lack of proper accounting of solvation thermodynamics, which may require incorporation of explicit water molecules against an ensemble of receptor poses to better represent the physical system [73]. Docking scores are also made uncertain by implementation choices made for algorithmic efficiency. For example, flexible ligand or receptor degrees of freedom are usually restricted to enable efficient conformational searching. Recent studies using deep learning neural-network-based scoring functions show promise in reducing prediction uncertainties [72,74]. Regardless of these inaccuracies, docking approaches are likely to remain an essential component of virtual screening, as they can reduce and enrich the chemical space of potential binders, and provide initial bound-pose conformations for seeding molecular dynamics (MD) simulations.
In contrast to docking, FEP methods are much more expensive, requiring all-atom MD simulations to perform statistical sampling of ligand-receptor configurations [75][76][77][78]. The benefit is a best-in-class prediction of ligand binding affinity, achieving accuracies within chemical accuracy, i.e. 1.0 kcal/mol [79,80]. Assigning uncertainties to such predictions, however, requires consideration of the many error sources these methods introduce. These can be categorized as arising from both sampling and scoring inaccuracies.

A. Uncertainty in FEP Estimates due to Sampling Error
To understand uncertainty from statistical sampling error, consider that using MD simulation naïvely to predict ligand affinity would require equilibrium sampling of multiple reversible binding transitions which typically exceed the timescales of state-of-the-art MD, even on special-purpose high-performance computers [81]. FEP and related methods circumvent this problem by performing simulations in multiple thermodynamic ensembles for a series of alchemical intermediates and analyzing the resulting samples with a statistical free energy estimator. Absolute binding free energy (ABFE) is computed along alchemical paths that decouple all non-bonded interactions in solution of the ligand, including the binding to a receptor. Relative binding free energy (RBFE) is computed along with an alchemical transformation from one molecular topology to another.
Statistical free energy estimators that utilize multiensemble sampling include exponential averaging (EXP) [82], thermodynamic integration (TI) [83], the weighted-histogram method (WHAM) [84], the Bennett acceptance ratio (BAR) method [85], the multistate Bennett acceptance ratio (MBAR) method [86], and newer dynamical estimators like TRAM [87] and DHAM [88]. Uncertainties from statistical sampling can influence the confidence of calculated Gibbs free energies in complicated ways. Hence, regardless of which free estimator is used, a non-parametric bootstrap is preferred to estimate uncertainty in Gibbs free energies [76]. Similarly, uncertainties can be calculated over block-averages of a long and equilibrated simulation trajectory, or can alternatively be estimated across multiple independent simulations, in some cases initiated with different ligand starting poses. For FEP, which involves computing free energy differences between the ligand-bound state and the free ligand in solution, overall uncertainties are often obtained by summing individual uncertainties for each ensemble in quadrature [71,76].
Different free energy estimators have unique bias and variance due to finite sampling, but they all depend on having sufficient thermodynamic overlap between neighbouring ensembles along the thermodynamic path. For this reason, relative free energy estimation is often more accurate for pairs of similar molecules, since the alchemical transformation involves a smaller number of atoms, and hence inherently has a better overlap. Achieving good thermodynamic overlap is the main consideration guiding choices about the number and type of alchemical intermediates, which may require position or angular restraints on the ligand to focus sampling on particular ligand poses [71].
Most free-energy estimators assume that samples are statistically independent and drawn from the equilibrium distribution. Sufficient equilibration can be notoriously slow to achieve for some systems, e.g. ligands with "hidden barriers" between conformational states, a problem which Hamiltonian replica exchange and other enhanced sampling methods have been designed to address [16].
To avoid statistically correlated input, which results in artificially low uncertainty estimates, recommended strategies discard initial unequilibrated samples, and subsample snapshots at a frequency (2τ + 1), where τ is a computed correlation time [89]. An algorithm to automatically detect simulation equilibration in this way has been developed [90].
Related to FEP is a class of non-equilibrium work (NEW) free energy methods, based on the Crooks fluctuation theorem [91], which drives a system between alchemical endpoints and analyzes the distribution of non-equilibrium work values [77]. These methods can achieve comparable accuracies to FEP with fewer intermediates [92][93][94], yielding uncertainties that depend on the variance of the work distributions, in turn depending on the timescale of the alchemical schedule. Sampling uncertainty can also be introduced by artifacts from the implementation of MD. For example, uncertainties on thermodynamic and kinetic properties can depend greatly on solvation artifacts that depend on simulation box size [95]. Furthermore, surprisingly, variation in computed energies across different simulation packages [96], integrators, and thermostats may be non-negligible [97].

B. Uncertainty in FEP Estimates due to Scoring Error
The workhorse of FEP is a molecular mechanics energy function, called a force field, which has usually been designed to approximate a quantum mechanical potential energy surface, and tuned to reproduce experimental observables. Iterative cycles of force field development, parameterization and assessment are thus vital for reducing the uncertainty introduced by force field inaccuracy [94,98,99]. With new machine learning advances have come new ways to improve force fields for the prediction of ligand binding free energies. One product of the recently-launched Open Force Field Initiative is a new generation of force fields that eschews atom-types in favour of direct chemical perception, utilizing the SMIRKS-native Open Force Field (SMIRNOFF) scheme for parameter assignment [100]. Parameterization of Open Force Field v1.0.0, code-named Parsley [38], was made possible by the ForceBalance algorithm, a Bayesian algorithm for automatic and reproducible parameterization [101][102][103]. Polarizable force fields such as AMOEBA continue to make gains in accuracy for predictions of Gibbs free energy of binding [104,105], and there are examples of neural network-based potentials, which aim to reproduce molecular properties with DFT accuracy at force field cost [106,107].
Finally, uncertainty can also result from the assumptions of a selected model. For example, most fixed-charge force fields require the pre-selection of particular protonation states for both ligand and ionizable side chains of a receptor. Depending on the predicted pK a of acidic or basic groups, which may differ significantly between bound and unbound states, multiple simulations may be required to explore these possibilities. Overall, there are several factors to consider when designing a study to compute binding affinities for drug-like molecules. New methodologies and best practices can often be found through blind assessment challenges such as SAMPL [108] (Statistical Assessment of Modeling of Proteins and Ligands), which are created annually to assess predictions of binding affinities among other properties important to drug discovery.

VI. UNCERTAINTY IN GENERATIVE MODELING
Generative models can produce synthetic data samples, some of them by learning from reference data. For molecular design, generative models are used to tackle the inverse design problem [109], where novel structures with desired properties are generated. Molecular property estimators are often used in conjunction with these models [110]. In this section, we discuss the importance of property prediction uncertainty in the context of generative modeling. In particular, we consider variational autoencoders (VAEs) [111], generative adversarial networks (GANs) [112], reinforcement learning (RL) [113], and genetic algorithms (GAs) [114].
VAEs possess the ability to take discrete structures such as molecules and convert them into high-dimensional latent representations by way of neural networks (referred to as encoders). Simultaneously, a mapping from the latent representation back to the discrete structure is trained using a separate neural network (referred to as a decoder). During training, VAEs learn to encode and decode correct molecules and simultaneously shape the latent space into a gaussian distribution for sampling. For molecular design, Gómez-Bombarelli et al. [115] introduced an additional neural network that predicts molecular properties from the latent space representation. During the training of this network, gradient optimization steps are taken for both the parameters of the encoder network and the property prediction network. This shapes the latent space based on the predicted property.
The formulation of VAEs as generative models brings about three sources of uncertainty for the associated property prediction model: (1) Decoding points from a latent space does not always lead to valid molecules. However, for every point in the latent space one can predict molecular properties via the prediction model. Hence, these regions can be viewed as uncertain spaces where predictions have to be disregarded due to the lack of a physical manifestation. (2) The encoding of structures to the latent space is stochastic. Kingma et al. [111] introduced a reparameterization step that ensures differentiability through the latent space for training the overall architecture and imposition of a gaussian distribution. As a result, identical molecules are typically encoded at a range of points in the latent space. While these encoded molecules should have the exact same predicted property, the range of positions tends to give rise to different predictions. Hence, there are multiple property predictions associated with any given molecule representing the corresponding uncertainty. (3) Discrete objects within VAEs, in particular molecular string representations and adjacency matrices, are usually decoded from the latent space as a time-series. These discrete objects are generated one at a time, by sampling over the distribution of all user-defined actions, i.e. units such as the allowed characters in a string. Typically, the unit with the largest probability is selected. However, multiple units can have identical probabilities. Consequently, decoding the same point from the latent space multiple times can result in distinct molecules. Additionally, as the same latent space point represents all these different molecules simultaneously, they will have the exact same predicted property.
GANs [112,116,117] are generative models with joint training of two competing networks, a generator and a discriminator. Like VAEs, they rely on a gaussian latent space to draw samples from and decode them using a neural network (referred to as generator network). The generated samples are compared to references from a dataset by another neural network (referred to as discriminator). The generator and discriminator compete with each other, the generator tries to produce samples that appear more realistic to the discriminator while the discriminator tries to improve its reliability to distinguish generated from reference samples. Similar to VAEs, in the traditional formulation of GANs, property predictors can be associated with the latent space. The sources of uncertainty are essentially identical to VAEs.
RL is a field where agents learn to navigate environments to maximize rewards. In molecular design, RL agents build chemical compounds [118][119][120], and rewards are provided for achieving desired molecular properties. Actions constitute making changes to chemical structures, while rewards are calculated on-the-fly based on the resulting structures. In contrast, GAs [121][122][123] are generative models that generate structures without gradient-based optimization. The algorithm is initiated with a set of molecules, i.e. the initial population, that are ranked by calculating their fitness. For molecular design, the fitness is typically assessed based on a set of target properties. The fittest molecules are carried over to the next generation, the unfit molecules are replaced with mutations and crossovers of the best. For both RL and GAs, to offset the often prohibitive computational cost of obtaining accurate property predictions, cheap approximations such as ML models are utilized to estimate rewards and fitness, respectively. By using cheap property predictors, in high-dimensional spaces, the generation of structures can be skewed to uncertain regions of the chemical space. For instance, certain structures can be adversarial examples or outliers for property predictors, where the model can be overly confident or uncertain about its prediction. Hence, the resulting molecule modifications are uncertain due to the uncertainty in the rewards. Notably, Thiede et al. [124] demonstrated that the uncertainty from prediction networks can be used directly in the reward function to guide the generation of structures outside of the confidence region of the prediction network for curiositydriven molecular design. Interestingly, Shen et al. [125] propose PASITHEA, where a property prediction neural network is directly used for generating molecules with deep dreaming. An interesting direction would be the incorporation of uncertainty within property predictors for structure generation via the deep dreaming algorithm.
Generative models coupled with property predictors are yet to achieve wide adoption due to their complex setup, their high computational costs and the need to adapt them to specific problems. When used for properties that are inherently complicated and hard to predict, the role of uncertainty and prediction confidence will be very important. In addition, in models with latent spaces such as VAEs and GANs, the choice of molecular representation influences both the validity of the latent space and thereby the prediction confidence. For instance, the use of VAEs and GANs was demonstrated in conjunction with the SMILES [58] representation, where the probability to produce invalid strings is high. Consequently, both VAEs and GANs trained on SMILES can produce invalid molecules. Importantly, this problem has recently been solved by the SELFIES representation [126], which guarantees molecular validity by design. While SELF-IES increases model confidence by generating only valid molecules, it does not remove uncertainty from inherent model stochasticities such as probabilistic encoding and decoding. Furthermore, confidence in molecular property predictors can be increased by extending datasets used for training and validation of property predictors systematically. Recently, STONED [127], a simple generative model producing molecules at high speed via systematic navigation within SELFIES, has been introduced. Accordingly, based on an initial dataset, it can be used to diversify the molecule space systematically leading to better property prediction models.

VII. CONCLUSION
In this review, we inspected common sources of uncertainties encountered in data-driven drug design workflows, some of which are well-studied in the literature, and some of which are typically neglected. Hidden biases in data can lead to overly confident property predictions that fail to manifest in the real world. Data-driven prediction models that explicitly incorporate uncertainty convey insights into the confidence of the model itself and provide a handle to assess the prediction quality. Accounting for uncertainty in the input features is underappreciated in ML research but is inherent to complex systems such as the chemistry of the human body. In contrast, error estimation in the prediction of binding free energies is one of the central methods in the field to assess the quality of simulation workflows that guides the development of efficient sampling techniques. Uncertainties in generative models are another underappreciated topic that needs to be addressed to improve the success of property-driven molecular design approaches. Overall, it is evident that assigning confidence in molecular prop-erty predictions is intimately tied to understanding the associated sources of error, and we believe that future progress in this field will require the development and application of even more methods accounting for uncertainties.

VIII. EXPERT OPINION
Assigning confidence to predictions is directly connected to uncertainty identification and quantification. When training ML models for predicting complex molecular properties, users need to be aware of potential biases of the training data when assessing generalizability. The recent use of DUD-E to train data-driven models serves as a cautionary tale for models learning to exploit biases in the data rather than learning the underlying physics of intermolecular interactions [24]. Most importantly, the smaller the training dataset size the larger the danger of hidden bias and the higher the importance of accounting for them. In that regard, transfer learning, especially in extremely sparse data regimes, requires more research to establish it as the method of choice for molecular property prediction.
Moreover, not only do small datasets lead to significant biases but they also lead to higher epistemic uncertainty because the model does not have enough information to learn the true relationship between input features and output labels [29]. However, expanding datasets, especially for properties relevant in drug design, can be prohibitively expensive and very time-consuming, in particular in the case of experimental data. Furthermore, not every new datum added to an existing dataset provides the same amount of additional information. Accordingly, developing systematic methods to identify data to be added to the existing dataset to maximize the improvement in prediction accuracy and the reduction in epistemic uncertainty is the subject of ongoing research.
Additionally, the aleatoric uncertainty [29] that is inherent in the data used to train ML models is at least as important as the epistemic one. One prominent example from the field of drug discovery discussed in this review is the use of estimated and simulated binding affinities to assess the interaction strength between ligands and receptors. Notably, binding affinities estimated via molecular docking are only of limited accuracy, but their computation efficiency allows assembling large datasets. The main question is whether the estimated property, regardless of how efficient it can be computed, conveys any information about the real relationship between input features and model output. When noise dominates over valuable information, the largest dataset will not provide useful information for ML models to be learned. In summary, it is important to build datasets and construct models that explicitly and consistently account for both aleatoric and epistemic uncertainties. While uncertainty estimates are relatively straightforward to obtain in general, they are only rarely benchmarked or evaluated during model training and future work needs to establish robust procedures for that purpose.
Another important and often underappreciated source of uncertainty, especially in complex systems, are the input features. This is particularly true for chemistry where molecular properties are generally conformer-dependent, and most molecules adopt much more than just one conformation under ambient conditions. Data-driven models, in contrast, are largely trained on 2D molecular representations that neglect conformations entirely. Accordingly, incorporating 3D structures and appropriate features in ML models is an active field of research [63]. The main challenge is to offset the added cost of generating representative 3D conformations with higher prediction accuracy in the final models. Currently, 2D molecular representations are still the state of the art in the field but we believe 3D features will become the standard in the near future.
One topic of paramount importance in drug design is simulating binding affinity between ligands and receptors. For that purpose, docking is a popular technique that provides crude estimates of interaction geometries and strengths using extremely efficient scoring approaches. However, this efficiency comes at the price of accuracy and, hence, cannot match the quality of FEP simulations, which are several orders of magnitude more computationally demanding. Importantly, even among FEP approaches there is an inevitable trade-off between cost and accuracy as the extent of sampling conducted directly controls the quality of the results. Accordingly, for extremely accurate results [71], FEP methods need to capture multiple binding transitions and this can increase the simulation cost by orders of magnitude making it less practical. Finally, a large portion of the uncertainties from both docking approaches and MD simulations stem from the scoring functions and the force fields used, respectively. While systematic optimization procedures for force field parameters are being developed and implemented, the development of force fields that not only provide energy but also uncertainty estimates has received less attention and we believe that future work in that direction is necessary to deliver more accurate models.
Finally, the generative models for molecular design introduced recently rely on efficient molecular property predictors to generate molecules with both desirable structure and function [110]. However, an explicit account and treatment of uncertainties, unfortunately, is largely neglected in the field. While accounting for property prediction uncertainties seems relatively straightforward with current methodologies, dealing with uncertainties in the molecular structure itself, and modifications thereof remains unexplored territory in comparison. Accord-ingly, we believe that future research in that direction is paramount and has the potential to lead to significant advances in computer-driven molecular design in general and drug design in particular.