Abstract Sequential experimental design to discover interventions that achieve a desired outcome is a key problem in various domains including science, engineering and public policy. When the space of possible interventions is large, making an exhaustive search infeasible, experimental design strategies are needed. In this context, encoding the causal relationships between the variables, and thus the effect of interventions on the system, is critical for identifying desirable interventions more efficiently. Here we develop a causal active learning strategy to identify interventions that are optimal, as measured by the discrepancy between the post-interventional mean of the distribution and a desired target mean. The approach employs a Bayesian update for the causal model and prioritizes interventions using a carefully designed, causally informed acquisition function. This acquisition function is evaluated in closed form, allowing for fast optimization. The resulting algorithms are theoretically grounded with information-theoretic bounds and provable consistency results for linear causal models with known causal graph. We apply our approach to both synthetic data and single-cell transcriptomic data from Perturb–CITE-sequencing experiments to identify optimal perturbations that induce a specific cell-state transition. The causally informed acquisition function generally outperforms existing criteria, allowing for optimal intervention design with fewer but carefully selected samples. Similar content being viewed by others Interpretation, extrapolation and perturbation of single cells Article 02 January 2026 Causal inference on human behaviour Article 23 August 2024 Systematic reconstruction of molecular pathway signatures using scalable single-cell perturbation screens Article 26 February 2025 Main An important problem across multiple disciplines, ranging from bioengineering to mechanical systems, operations research and environmental regulation, is the discovery of interventions on a system that can produce a desired outcome. With little prior knowledge of the outcome before performing the intervention, one can have a huge number of possible choices for the optimal design. In particular, the interventions in many applications are combinatorial, resulting in an exponential-size design space, making an exhaustive search infeasible. Examples include experimental design of genetic perturbations, such as those for cellular reprogramming in regenerative medicine1, optimal feedback control in mechanical systems2 as well as turbulent flows3, dynamic pricing strategies in customer networks4 and iterative intervention research for climate change adaption5. In this context, ‘active learning’ has been proposed as a machine-learning strategy to efficiently explore the search space6. Such methods sequentially and strategically acquire new interventions, with the goal of discovering an optimal or a close-to-optimal intervention using the fewest number of samples. Although it can be of interest to identify optimal interventions for ‘estimating’ particular quantities in the model (for example, ref. 7), in this work, we consider the optimality of an intervention with respect to ‘optimizing’ its effect. Specifically, this is done by successively (1) updating the model belief using samples acquired so far from different interventions and (2) selecting the next intervention to obtain samples from by constructing and optimizing an acquisition function, which prioritizes interventions that are more informative for the desired outcomes (Fig. 1). Fig. 1: Overview schematic of the active learning framework for optimal intervention design in causal models. Fig. 1: Overview schematic of the active learning framework for optimal intervention design in causal models. Full size image Iterative process of active learning for intervention design in causal models, in which the main design steps are to update the structural causal model using the obtained samples and to select the next intervention based on an acquisition function. Standard approaches towards this problem are correlation based. More precisely, the idea is to estimate associational relations between intervention and outcome to update the model belief and make decisions about which samples to acquire next. The two main approaches use either statistical theory, by minimizing the posterior variance of the outcome estimate8, or information theory, by maximizing the mutual information between its samples and the quantity of interest9. Although these methods are being widely applied, correlation-based approaches are not optimal when the underlying model is causal, as they do not take into account the structural information that can reduce the number of feasible models. Many systems that are relevant for applications are causal, in which an intervention can only have an impact on downstream variables. It is therefore of interest to develop methods that learn what is necessary about the underlying causal mechanisms to identify optimal interventions more quickly. These limitations of correlation-based approaches have been noted in the related bandit setting, in which the goal is to minimize the cumulative regret by selecting arms iteratively, and causal relations have been used to improve over standard regret bounds10,11. Recent works in the related field of Bayesian optimization have also considered exploiting causal structure12,13,14. We provide a detailed review of these works in Supplementary Information I. Our work extends these works in two ways: (1) rather than optimizing a single target node, we optimize the entire distribution mean, and (2) rather than considering discrete or finite interventions, we consider continuous-valued interventions. This is important for various applications, such as optimizing drug dosages or product prices. As a concrete example, consider cellular reprogramming in genomics, a problem of great interest for regenerative medicine1,15. The aim in this field is to reprogramme easily accessible cell types into a desired cell type via continuous-valued interventions such as the overexpression of particular transcription factors and genes. Because such interventions act on genes, which regulate each other through different pathways16, this problem can be formulated as optimal intervention design in a causal model represented by a directed network on genes. In particular, we model the underlying causal model in a Bayesian way using a structurally informed prior. We then construct an acquisition function based on this model and show how to efficiently evaluate and optimize it. The acquisition function both enjoys an information-theoretic bound and provably recovers the optimal intervention in the appropriate limit. We demonstrate experimentally that our algorithms outperform baselines on both synthetic data and the design of genetic perturbations in the context of single-cell gene expression data. Finally, we conclude with an outlook to future research directions and discuss other potential applications. Problem set-up The state of the system of interest is described by a p-dimensional random variable sampled from a structural causal model (SCM)17,18. Precisely, the ‘causal structure’ of the system is represented by a directed acyclic graph (DAG), and we assume that the joint distribution P of x factorizes with respect to the DAG; that is, , where pa(i) = {j ∈ [p]: j → i} denotes the parents of node i in the DAG. In this formulation, the conditional distributions are the ‘causal mechanisms’ that generate the variable xi from its parents xpa(i). Assuming a linear Gaussian model, then (1) where the real-valued coefficients Bik = 0 if k ∉ pa(i) and the exogenous noise variables with variance are mutually independent. An example of this model is given in Fig. 2. For simplicity, we assume that the system is centred to be mean zero, but intercepts can be easily added to equation (1) and the following approaches still apply. Fig. 2: Example causal model. Fig. 2: Example causal model. Full size image An example of a linear SCM with Gaussian noise on a four-node DAG, in which nodes and edges are labelled with variables and coefficients, respectively. An intervention, denoted by a vector , modifies the conditional distribution P(xi∣xpa(i)) into a new conditional distribution, Pa(xi∣xpa(i)). Every i ∈ [p] for which ai ≠ 0 is called an ‘intervention target’. In this work, we consider shift interventions19,20 (which are a special class of soft interventions21), in which the interventional distribution Pa(x) of the modified system under intervention is given by (2) where ai = 0 for every i that is not an intervention target. This can be written in matrix form as x = Bx + a + ϵ, or x = (I − B)−1(a + ϵ), where I is the pth order identity matrix and with Σ being the diagonal matrix with on its diagonal. These interventions can be used to model a broad class of genetic perturbations including clustered regularly interspaced short palindromic repeat (CRISPR) interference and activation22, and transcription factor overexpression23. Although we consider the case in which all variables can be intervened on, as for example in genetic perturbation experiments24, for other applications, it may be of interest to consider extensions in which only a subset of variables can be intervened on. We also note that an extension to hard interventions can be obtained by setting Bik to 0 for every i that is an intervention target and subsequently defining the post-interventional mean. Whether or not an intervention induces the desired outcome is decided through samples obtained in the intervened environment, that is, x ~ Pa. Since x is random, it is intuitive to use the average over multiple samples20,25, that is, the empirical estimate of the distribution mean, and compare this with the ‘target mean’ μ*, which is user specified and describes the desired outcome. Formally, we seek an intervention such that, after obtaining n samples , which we denote by x[n], the squared distance between the empirical mean and the target mean, that is, the quantity , is minimized. Denoting by ϵ(1),…, ϵ(n) (or short ϵ[n]) the exogenous noise vectors independently sampled from , then the squared distance between the empirical and target mean can be written as (3) using the matrix form. Here accounts for the finite number n of interventional samples, where n is user specified based on the available budget. However, as ϵ[n] does not depend on the choice of a, we can discard this term for the minimization of equation (3) with respect to a by considering its infinite-sample version. When the sample size goes to infinity, by the law of large numbers, almost surely, and thus, the optimal intervention a* achieves the minimum value of zero in equation (3) and has an explicit form a* = (I − B)μ*, which, as expected, depends on the unknown parameter B. In what follows, we assume that the DAG structure is given; that is, the sparsity pattern of B is known but not the edge weights Bik for Bik ≠ 0. This assumption is natural in many applications including fluid mechanics26 and optimal pricing4, in which there are pre-specified networks given by either known laws or prior information. For applications in which the network cannot be assumed to be known such as problems in biology, the common approach is to use existing data from various sources to learn the DAG first15,27. Further discussions on the implications of an unknown DAG structure are provided in Supplementary Information D. Let be the ‘current’ dataset, consisting of all samples obtained so far by performing interventions a(1),…, a(t); here, n1,…, nt denotes the number of samples obtained for each of the t interventions. To simplify notation, we will assume n1 = ⋯ = nt = n in the following, but all results still hold when using different sample sizes for each intervention. Given this dataset, the goal is to select the next intervention such that the resulting dataset contains as much information as possible about the underlying optimal intervention a*. The overall aim of iteratively picking interventions is to find the optimal intervention a* with a minimum number of samples. Design of the acquisition function To sequentially select the next best intervention, there are two important steps (Fig. 1): (1) updating the posterior of the edge weights in the causal model based on the samples in collected so far and (2) constructing an ‘acquisition function’, , such that the new dataset after adding samples from is most informative of the optimal intervention a* and can be evaluated and optimized efficiently. For the first step, we generalize the DAG–Wishart distribution28,29,30 to define a Bayesian model on the parameters B, which can be updated efficiently given a dataset (Methods). For the second step, we first characterize the uncertainty in estimating the optimality of an arbitrary intervention a. Recall that the optimality of a is given by the square distance in equation (3), which measures how close the intervention is to achieving the target mean. As B is unknown, we can only estimate this square distance based on the current collected samples in . The uncertainty of the estimation can be characterized by its variance . However, this quantity is typically hard to evaluate as it involves (I − B)−1, whose posterior does not have a closed form. We instead multiply the term inside the variance by I − B and characterize the following variance: (4) where we define (5) Interestingly, g(a) can be interpreted as a noisy version of the ‘optimality gap’ by noting that the gap between an arbitrary intervention a and the optimal intervention a* = (I − B)μ* can be written as , which is a version of g(a) without noise terms. Building upon equation (4), the next acquired intervention should be such that after adding its samples to to obtain , the uncertainty conditioned on is minimized. However, as the samples are unobserved before performing the intervention, we do not have access to yet. Therefore, when deciding which intervention to perform, we consider the a′-augmented dataset with ‘hypothetical samples’ of a′, which are n repetitions of the plug-in estimator defined as . This estimator is obtained through by using the maximum a posteriori estimate , which concentrates to B as grows31, and replacing ϵ with its mean 0. By denoting the feasible set of interventions by , a reasonable choice is to select a(t+1) based on integrating the uncertainty over all . Formally, we let (6) and define the causal integrated variance (CIV) acquisition function h as follows. Definition 1 The CIV acquisition function evaluated at a′ with current dataset is (7) where ν is a non-negative measure on . Intuitively, this acquisition function provides a one-step lookahead of the overall uncertainty after acquiring intervention a′. Minimizing it will prioritize interventions that are most informative towards estimating the optimal intervention. This acquisition function also automatically accounts for the causal model by using a posterior on B. Note that in this formulation, we can choose the measure ν. For example, a uniform measure treats each intervention equally and the resulting uncertainty captures how well we can estimate the entire landscape of the optimality gap g. In most cases, an overly concentrated measure (for example, a Dirac measure at a single point) is not preferred, as it can lead to the erroneous estimation of g(a) for most , which makes minimizing g hard. Inspired by a recent line of work on output-weighted acquisition functions32,33, we describe how to choose a non-uniform measure ν in the next section. We discuss how to optimize CIV to solve for a(t+1) in Methods, in which we show that the variance can be computed in closed form. Considering to be the unit hypersphere and ν to be the uniform measure on it, this then leads to an explicit formula for CIV, which enables fast gradient-based optimizers to be used. Making the acquisition function output weighted Although the use of a uniform measure ν places an equal weight on reducing the variance of estimating the optimality gap for all a in , as our goal is to identify the optimal intervention (that is, which minimizes the optimality gap), it is desirable to place more weight on interventions a in with a smaller optimality gap. Note that as the ambient dimension grows, the volume (and thus the probability) of interventions with the optimality gap under a certain threshold shrinks (Supplementary Information C). This motivates the following measure, which uses the inverse of the optimality gap probability to up-weight interventions a in that are closer to the optimal intervention: (8) where is the estimated optimal intervention, and the probability density function is the distribution on the optimality gap induced by the uniform distribution fa on . Figure 3a shows that, as desired, the proposed weighting puts more mass on interventions that are closer to b as the dimension increases. When p = 3, however, ν degenerates to a uniform weighting (Supplementary Information C). In this case, we can use the bimodal von Mises–Fisher distribution, which behaves similar to equation (8) in higher dimensions (Fig. 3b and Supplementary Information C). Note that the weighting proposed here is symmetric and also puts more mass on interventions that are closer to −b. The reason for this is that the optimal intervention can also be recovered by maximizing the optimality gap ∥a − a*∥2 which gives −a* (Supplementary Information C). Fig. 3: Illustration of output-weighted non-negative measure ν(a) on the space of all possible interventions . Fig. 3: Illustration of output-weighted non-negative measure ν(a) on the space of all possible interventions $${{{\mathcal{A}}}}$$ A . Full size image a, Values of da/dν(a) (up to constant multipliers) plotted against values of on the half sphere with . Solid lines correspond to in different dimensions; the dashed line is the non-degenerate replacement for p = 3 using the bimodal von Mises–Fisher distribution. b, Visualization of ν(a) in three dimensions using the bimodal von Mises–Fisher distribution. Here was randomly generated and the points on the sphere are coloured corresponding to the value of dν(a)/da, which is larger for directions that are more aligned with b. The corresponding causal output-weighted integrated variance (CIV-OW) acquisition function is given by (9) Methods for evaluating and optimizing CIV-OW are given in Supplementary Information C and Supplementary Figs. 1 and 2. Theoretical results We provide two interpretations of the introduced CIV acquisition function. First, taking an information-theoretic perspective, we show that the proposed uncertainty measure can be lower bounded by the ‘negative’ mutual information between the variables of interest and the newly acquired samples. Thus, minimizing the uncertainty corresponds to maximizing a lower bound to the mutual information, which means that CIV approximately prioritizes the most informative interventions. Second, taking a graphical perspective, we illustrate how the causal structure is utilized by the CIV acquisition function to identify an intervention that is asymptotically consistent with the optimal intervention. The details underlying these theoretical results are provided in Methods, and Fig. 4 illustrates our consistency results experimentally, showing that a gradient-based optimizer with initialization close to the optimal intervention a* converges to a* as t increases. Fig. 4: Convergence of the selected intervention a(t) to the optimal intervention a*. Fig. 4: Convergence of the selected intervention a(t) to the optimal intervention a*. Full size image a, DAG showing the targets of the optimal intervention a*. b, Minimizing the CIV acquisition function using a sample size of n = 1 at each time step yields interventions a(t) for t = 1,…, 50 that converge to a*. Applications Experiments on the synthetic dataset To create a generative model following the linear SCM in equation (1), we first generate a DAG with p nodes (Supplementary Information G and Supplementary Fig. 3). Then, we randomly draw edge weights, namely, Bik with k ∈ pa(i), from a uniform distribution bounded away from zero. Next, we generate a sparse set of intervention targets and a randomly sampled optimal intervention a* over these targets. Finally, we calculate the target mean μ* using a* and the ground-truth causal model. These steps construct a synthetic instance of a causal system and optimal intervention a*. A more detailed description of the above procedure is given in Methods and Supplementary Information G. We compare our two acquisition functions, CIV in equation (7) and CIV-OW in equation (9), against four relevant baselines. The ‘random’ baseline corresponds to a passive setting, in which each intervention is selected at random and no information from the collected samples is used. We also compare against three other active methods. The ‘greedy’ baseline selects the next intervention purely based on the current estimate of a* (which is given by , where B is estimated from ). The MaxV baseline seeks to select interventions that minimize the posterior variance of the estimate of the model parameters B. It uses as acquisition function a scalar version of this variance, namely, , where ∥⋅∥2 denotes the spectral norm of the covariance matrix of Bi,pa(i). As different rows of B (for example, Bi,pa(i) and Bj,pa(j) for i ≠ j) are independent (Definition 2), we use the maximum of the spectral norms over i ∈ [p]. Whereas our acquisition functions concentrate on estimating the model parameters that are relevant for a*, this baseline estimates the entire model. Finally, the CV baseline seeks to minimize the posterior variance of estimating a* and uses the spectral norm of its covariance matrix as acquisition function; that is, . This is in contrast to our approach that integrates the estimation uncertainty over the entire feasible set of interventions. To reduce evaluation noise effects, we run each method 20 times over 10 instances of a fixed DAG and optimal intervention a*. Figure 5 shows our results for 30-node DAGs with 10 intervention targets and sample size n = 1 per time step. Denoting the obtained samples at time step t with , we use the unbiased estimate of a* to obtain the current estimate of the target mean using the true model parameters, that is, . Figure 5a shows the decline in the relative distance between current and target mean across time steps. Figure 5b highlights the statistics of the last time step. As is apparent from these experiments, our proposed acquisition functions consistently outperform all baselines, with CIV-OW improving upon CIV by using an output-weighted measure. In Supplementary Information G, we examine the effect of varying different parameters including graph size, number of intervention targets and graph structure, demonstrating the robustness of our results (Extended Data Figs. 1–4). We also include three additional baselines from prior works9,12,34,35 (Supplementary Information G and I) on a ten-node DAG (Extended Data Fig. 5) and an experiment to examine the effect of misspecifying the underlying causal structure (Supplementary Fig. 4 and Extended Data Fig. 6). These results further demonstrate the effectiveness of our proposed method in terms of both accuracy and efficiency. Fig. 5: Comparison of different acquisition functions (random baseline, greedy baseline, MaxV baseline, CV baseline and our proposed CIV and CIV-OW) in a simulation study. Fig. 5: Comparison of different acquisition functions (random baseline, greedy baseline, MaxV baseline, CV baseline and our proposed CIV and CIV-OW) in a simulation study. Full size image a,b, Simulations on 10 instances of a 30-node DAG in which the optimal intervention a* has 10 targets. Each method is run 20 times and averaged. a, Relative distance between the target mean μ* and the best approximation across time step t. Lines denote the mean over the ten instances; shading corresponds to 1 s.d. b, Relative distance statistics of each method averaged over ten instances at the last time step in a. Small circles: relative distance statistics of the corresponding method on a specific instance. Whiskers: first quartile minus 1.5× the inter-quartile range to third quartile plus 1.5× the inter-quartile range. Experiments on the biological dataset We next study the performance of our method to identify the optimal intervention for inducing a desired cell-state change in human melanoma cells, thereby mimicking a cellular reprogramming experiment. For this, we use Perturb-CITE-sequencing(Perturb–CITE-seq) data from ref. 36 consisting of single-cell transcriptomic readouts for a large collection of patient-derived melanoma cells. Here an intervention is a genetic perturbation that targets one or multiple genes and drives the expression of these genes towards zero (more precisely, these are knockout interventions). Cell states are measured by the joint distribution of the expression of a collection of genes. Samples from these distributions correspond to gene expression vectors of individual melanoma cells. To avoid dealing with batch effects, we use only one of the screens from ref. 36 (namely, the control screen with no additional treatment), which contains 5,039 cells with no perturbation and 30,486 cells with interventions on subsets of 248 genes associated with immunotherapy resistance (Supplementary Fig. 5). Gene expression of each sampled cell is captured as a vector of log-transformed transcripts per million. For our optimal intervention design task, we focused on a particular functional context, namely, the p = 36 genes among the 248 interventional targets that are involved in interferon-γ signalling and immune response36 (Supplementary Information H). As multi-target interventional samples are extremely scarce in this dataset, with most such interventions having no more than one sample (Supplementary Fig. 6), we only consider single-target interventions. We use the observational samples, that is, the cells with no perturbation, to learn a DAG over these 36 genes (Methods, Supplementary Fig. 7 and Extended Data Fig. 7). We model each of the 36 single-target interventions as a downshift of the target gene by its observational mean; that is, knocking out target gene i ∈ [p] is modelled as a shift intervention a with ai being the negative of the observed mean of this gene, and all other entries being zero. This model assumes that the perturbations are effective at knocking out their target genes. This assumption is not always met by current technologies: we observe that some interventions do not downregulate their target genes in a statistically significant manner (Supplementary Information H). However, this is only known after performing the experiments, and thus, we use the idealized model. If another dataset or prescreen (usually on a smaller scale) is available to help determine the effectiveness of an intervention to target a specific gene, then we can easily replace the idealized model with this estimate. As we will show, the idealized model still extracts enough signal for our method to aid in finding good interventions. To evaluate our approach, we use the observational distribution as the source cell state and a particular single-node interventional distribution as the target cell state. The aim is to identity this single-node intervention or an intervention with similar effect using the least number of samples. For benchmarking purposes, we use a setting in which the optimal intervention is contained in the feasible set of interventions, but we note that this is not a requirement and the target cell state can be any desired distribution, the experimental design goal then being to identify an intervention that moves the distribution as close as possible to the target distribution. Figure 6a shows three representative examples of interventional distributions (targeting CDH19, SOX4 and HLA-C), visualized using a 2D projection along the most variable directions obtained using contrastive principal component analysis (PCA) (Methods). These examples highlight that the alterations induced by single-target interventions are subtle (Supplementary Fig. 8), with SOX4 being more similar to CDH19 than HLA-C as also corroborated by the squared distance between the interventional means (Fig. 6b). This observation is consistent with previous melanoma studies showing that HLA-C is associated with positive immune response37, whereas CDH19 and SOX4 are associated with metastasis and immune evasion38,39. Figure 6c shows the performance of our results with the optimal intervention targeting gene CDH19. More examples using other target genes are given in Supplementary Information H. Each method is run 50 times, in which each run starts with a warm-up set of 100 observational samples and then n = 10 interventional samples per time step. In Fig. 6c, we present the comparison between three acquisition functions: CIV and the two baselines random and greedy. More implementation details as well as the full results are given in Extended Data Figs. 8–10. Similar to what we observed in the synthetic data experiments, CIV outperforms the benchmarks in terms of distance to target mean across all time steps. This suggests that the proposed approach is beneficial for identifying perturbations to induce a desired cell-state change. Fig. 6: Results on perturbational single-cell gene expression dataset. Fig. 6: Results on perturbational single-cell gene expression dataset. Full size image a, Kernel density estimate plot of the two-dimensional projection of the interventional distributions obtained by targeting three representative genes: CDH19, SOX4 and HLA-C. Colors indicate the density of the interventional distribution targeting the corresponding gene. b, Squared distances between the target mean of the optimal intervention targeting CDH19 and other interventions. c, Comparison of different acquisition functions for identifying interventions that match the target mean presented as mean ± s.e.m. Discussion In this work, we developed an active learning framework for optimal intervention design in causal models. Our method has two main ingredients: (1) modelling and updating the edge weights in the causal model using a Bayesian approach (using the DAG–Bayesian linear regression (BLR) distribution) and (2) optimizing the next intervention from which to obtain samples using a class of causally aware acquisition functions (CIV acquisition functions). The DAG–BLR distribution respects the underlying causal structure and allows for efficient posterior updates. The proposed CIV class of acquisition functions prioritizes the most informative intervention with respect to identifying the optimal intervention for moving the system towards a desired mean by minimizing an uncertainty quantity weighted based on the directions of interest. Importantly, the designed CIV acquisition function allows for efficient optimization by having tractable closed-form evaluations. In addition, we showed that the introduced acquisition function is characterized by attractive theoretical properties, such as by mutual information bounds and consistency. Finally, we demonstrated the developed active learning framework on both synthetic data and a biological dataset. In both cases, the designed acquisition function outperforms empirical analogues, allowing for accurate predictions with fewer experiments. We made various assumptions that may be limiting for some applications and motivate future research directions. First, we focused on the setting with known causal structure. An important future research direction is to consider the case in which the causal structure is partially or entirely unknown. We discuss potential avenues to approach this problem in Supplementary Information D. Second, we considered linear SCMs with additive Gaussian noise. A potential extension to the nonlinear setting could be achieved using kernels with linearity over the feature space. Similar derivations to our work could be used to evaluate or approximate the CIV acquisition functions in this setting. Third, in our model, we assumed ‘causal sufficiency’17, which excludes the existence of latent confounders as well as the possibility of performing interventions only on some nodes. This is violated when some system variables are unobserved and thus cannot be intervened on or specified in the desired state. A possible approach is to use a more agnostic model between interventions and their effects as proposed, for example, in ref. 12; however, this will generally lead to weaker results and a loss of structural information. Although we discussed our work in the context of applications to cellular reprogramming, we envision our framework to be applicable broadly for sequential design problems arising in complex systems. In Supplementary Information J, we discuss several other applications and how they fit into our proposed framework. Methods Posterior update of edge weights The current dataset induces a belief on the model parameters B. We consider the Bayesian setting in which the belief corresponds to a distribution. To account for the known causal structure, we assume a generalization of the DAG–Wishart distribution28,29, which places a prior on B that respects the causal structure and allows for efficient posterior updates (Supplementary Information A). The prior is as follows. Definition 2 The DAG–BLR prior models B and the noise variances Σ jointly as , where each is (10) Here denotes the inverse-gamma distribution and are hyperparameters satisfying certain constraints specified in Supplementary Information A. This prior is consistent with the DAG structure as it sets Bik = 0 for all k ∉ pa(i). Although it was developed in the observational setting40, the following lemma shows that it can be extended to the interventional setting and is a conjugate prior for the model in equation (2). Thus, the posterior lies in the same family of distributions. Lemma 1 The posterior corresponding to the DAG–BLR prior satisfies , where is (11) The hyperparameters are specified in Supplementary Information A. The hyperparameters can be updated easily using the new samples obtained at step t; thus, this choice of a ‘conjugate’ prior allows for efficient posterior updates. We make one comment about the DAG–Wishart distribution. In previous literature28,29, this distribution has been extended to model the posterior beliefs of both the edge weights and the DAG structure. The posterior for the DAG structure is known as the Bayesian Gaussian equivalence score. Thus, although we here assume that the DAG is known or prefixed and we only update the posterior on the non-zero entries of B, our framework can be extended to the unknown DAG setting by placing a probability on the DAG structure using the Bayesian Gaussian equivalence score, as discussed in Supplementary Information D. Evaluation of CIV in closed form We discuss how to optimize CIV to solve for a(t+1). We start by providing a closed form of the variance . For ease of reading, we here provide the formula for the special case in which the exogenous noise variances Σ are known. The general formula, which is similar in flavour but more complicated, including the proof is given in Supplementary Information B. Proposition 1 Conditioning on Σ, we have (12) where only depends on the augmented dataset , whereas and the constant c do not depend on a′. We note that the resulting expression is a function of bi and vi, which can be interpreted as follows: Lemma 1 can be used to rewrite as the estimated shift of the optimal intervention at node i, and as the covariance matrix of Bi,pa(i) scaled in the direction of the target mean μ*. Next, we discuss how to integrate equation (12) over the space of possible interventions to evaluate . Here we consider to be the hypersphere and ν to be the uniform measure on the hypersphere (other types of feasible sets , for example, with sparsity constraints, can be considered, and similar derivations can be used to identify if a closed-form integration exists). This corresponds to fixing the magnitude of the intervention and only optimizing over its direction. This is suitable for various applications including the biological problem considered below, as the strength of the intervention is often designed separately or prefixed manually by the experimenter. We also note that for many linear problems, the uncertainty is decreased by making the magnitude of the selected point larger32. We show in Supplementary Information B that this is the case also for the problem considered in this paper. Considering the hypersphere allows us to avoid this ambiguity. The following proposition provides the resulting formula when Σ is known; see Supplementary Information B for the general case including the proof. Proposition 2 For and ν being a uniform measure, the CIV acquisition function evaluated at a′, conditioned on Σ, is (13) where c1 > 0 and c2 are constants that do not depend on a′. The gradient of this objective function can be calculated explicitly also in the general case when Σ is unknown, and thus, we can use gradient-based optimization methods with a ball constraint to solve for (see the section Implementations below for details). For other types of feasible sets , for example, with sparsity constraints or consisting of hard interventions, the appropriate optimization method needs to be adjusted accordingly. We also note that the objective function is not necessarily convex (Supplementary Information B), and gradient-based optimizers may therefore only find a local minimum. Mutual information bound To provide an information-theoretic interpretation of the CIV acquisition function, we use the ‘relative decay’ of the uncertainty measured by (14) where the expectation is taken with respect to a new sample x′ from a′ whose distribution is given by . To simplify the notation, we restrict the discussion in this section to the case of n = 1. Similar results can easily be derived for sample sizes n > 1. The following theorem shows that the mutual information between g(a) and a sample x′ from a′ can be lower bounded by this relative decay up to a multiplicative factor. The proof can be found in Supplementary Information E. Theorem 1 Conditioning on , the following inequality holds for any for the mutual information between g(a) and a sample x′ from Pa′: (15) where ρ2 is a constant that does not depend on a′. This result holds in general for any function g(a) and its posterior modelling on , but the constant ρ in the bound depends on the specific choice of g and its posterior modelling. For example, for the prominent BLR task, Theorem 1 holds with ρ2 = 1/2 (Supplementary Information E). Note that as does not depend on , minimizing the uncertainty is equivalent to maximizing its relative decay in equation (14). The uncertainty used by the CIV acquisition function (equation 7) is a computationally efficient estimator of in the relative decay (equation 14), where the expectation with respect to is replaced by a plug-in estimator . It follows from Theorem 1 that minimizing the integrated variance in equation (7) corresponds to selecting an intervention that maximizes a lower bound on the mutual information between the resulting samples and g(a). This perspective shows how the CIV acquisition function connects to a prominent line of previous works, known as Bayesian active learning by disagreement (BALD)7,9,41, in which the idea is to select interventions that maximize the information gain at each step. For example, if the quantity of interest is g(a), then BALD seeks to find a′ such that the information gain measured by the entropy decay (16) is maximized. Note that this term equals . As a consequence, BALD directly operates on the mutual information, which involves high-dimensional moments of the random variable g(a). This highlights an important drawback of BALD, namely, that it is generally computationally difficult to evaluate and often requires approximation techniques such as Markov chain Monte Carlo9. By contrast, our CIV acquisition function only depends on the first and the second moments of g(a), thereby allowing for direct computation and efficient optimization. Consistency of CIV for identifying the optimal intervention We now provide a graphical interpretation of the proposed CIV acquisition function. For simplicity, we restrict our discussion to the expression developed in Proposition 2 for the case of uniform measure and known variances. The general setting is discussed in Supplementary Information F. As c1,c2 in Proposition 2 are constants that do not depend on a′, the new intervention is obtained by minimizing ; with a slight abuse of notation, let (17) The term can be understood as the epistemic uncertainty of node i in the direction of interest, namely, as , the variance of Bi,pa(i) characterizes the epistemic uncertainty of estimating xi using its parent nodes; in vi, this uncertainty, presented as a covariance matrix, is transformed into a scalar value by scaling it in the direction of the target mean μ*. As the number of time steps t increases, one can show (Supplementary Information E) using the Bernstein–von Mises theorem31 that vi = O(1/t), which implies that the second-order terms are dominated by the first-order terms in equation (17). Note that their coefficients are larger for nodes i with larger noise variances and larger estimated optimal shift values . This scaling is intuitive, as it puts emphasis on reducing the variance of nodes in which there is still a high uncertainty (given by ) or that require high shift values (estimated by ). As we show in the following theorem, the combined effect of these terms on acquiring new interventions is that a* will become an approximate local minimizer of equation (17) as t increases as shown in the following theorem, in which we use the conventional asymptotic notation, o, O and Θ, with respect to the time step t. Theorem 2 For all , the CIV acquisition function either decays linearly with respect to t, that is, , or degenerates to the constant zero, that is, . The optimal intervention a* satisfies and . This theorem shows that as t increases, the gradient at a* vanishes to zero faster than the acquisition function and the Hessian becomes almost positive semi-definite. This suggests that gradient-based optimizers with an initialization close enough to a* converge to the optimal intervention a* as t → ∞. In experiments, we observe this to hold even for a moderate number of time steps t as illustrated for a ten-node DAG example in Fig. 4. The proof of Theorem 2 can be found in Supplementary Information F. Implementations Optimization of the acquisition functions was performed using a gradient-based solver, more precisely, sequential least squares programming with the nonlinear constraint ∥a(t)∥2 ≤ 1. As proven in Supplementary Information B, optimizing the proposed acquisition functions with this constraint always outputs a feasible a(t) such that ∥a(t)∥2 = 1. For other types of feasible sets, for example, a(t) ∈ [0,1]p, one can use other solvers such as limited-memory Broyden–Fletcher–Goldfarb–Shanno with box constraints. Sequential least squares programming was implemented using the SciPy package42 and initialized in two ways: (1) using the intervention a(t−1) that was selected in the previous time step and (2) using the current estimate of the optimal intervention. Of the two resulting solutions, we used the one with the better acquisition function value as a(t). For the experiments on synthetic data, we generated the DAGs using the NetworkX package43 and the CausalDAG package44. We used different graph types and sizes, as detailed in Supplementary Information G. The edge weights were uniformly sampled from [−1,−0.25] ∪ [0.25,1] to ensure that the parameters were bounded away from zero. The exogenous noise levels were set to 1, and the resulting linear Gaussian SCM was then rescaled with our implementation of the standardized model as described in Appendix F of ref. 45. For the DAG–BLR prior (equation 10), we set the edge-weight-related hyperparameters as follows: mi = 0 and Mi equal to the identity matrix for all i ∈ [p]. We assumed the variances to be known and used the known-variance formula. The transcriptomic dataset analysed in the biological application was processed using the Scanpy package46. The kernel density estimate plots shown in Fig. 6a share the same coordinate axes and were produced by projecting the high-dimensional dataset along the first two principal components obtained using contrastive PCA47, with the unperturbed samples as ‘background’ and all interventional samples as ‘foreground’. The purpose of contrastive PCA is to identify low-dimensional structures in a foreground dataset relative to a background dataset. The greedy sparsest permutation algorithm48 was used to learn the DAG structure from the unperturbed data. The parameters used in the greedy sparsest permutation are given in Supplementary Information H and Fig. 7. For the DAG–BLR prior, we set the hyperparameters as follows: mi = 0, Mi equal to the identity matrix, αi = 2 and βi = 0. These hyperparameters were then updated using posterior formula 1 based on 100 samples from the unperturbed data. Similar warm-up steps are used in the active learning literature7,32. As we usually have access to some observational samples in applications of optimal intervention design, we recommend using these to obtain dataset-informed hyperparameters of the DAG–BLR prior. After this warm-up step using unperturbed data, interventions were acquired based on the known-variance formula, in which we used the estimated mean of the noise variances as a proxy. The codebase for updating the model posterior and optimizing the causally aware acquisition functions proposed in this work can be obtained via refs. 49. We also provide a notebook to extract and process the Perturb–CITE-seq data from ref. 36. Our codebase can be used to replicate all the main results and figures as well as for other user-defined applications. All methods can be run efficiently on a CPU for a moderate number of variables.
