Active Multi-Information Source Bayesian Quadrature
Alexandra Gessner Javier Gonzalez Maren Mahsereci

Tightness-aware Evaluation Protocol for Scene Text Detection
    Yuliang Liu, Lianwen Jin^*, Zecheng Xie, Canjie Luo, Shuaitao Zhang, Lele Xie

College of Electronic Information Engineering, South China University of Technology 

liu.yuliang@mail.scut.edu.cn; lianwen.jin@gmail.com











    December 30, 2023
=======================================================================================================================================================================================================================================




Bayesian quadrature () is a sample-efficient probabilistic numerical method to solve integrals of expensive-to-evaluate black-box functions, yet so far, active  learning schemes focus merely on the integrand itself as information source, and do not allow for information transfer from cheaper, related functions.
Here, we set the scene for active learning in  when multiple related information sources of variable cost (in input and source) are accessible.
This setting arises for example when evaluating the integrand requires a complex simulation to be run that can be approximated by simulating at lower levels of sophistication and at lesser expense.
We construct meaningful cost-sensitive multi-source acquisition rates as an extension to common utility functions from vanilla  (), and discuss pitfalls that arise from blindly generalizing.
In proof-of-concept experiments we scrutinize the behavior of our generalized acquisition functions.
On an epidemiological model, we demonstrate that active multi-source  () allocates budget more efficiently than  for learning the integral to a good accuracy.

















§ INTRODUCTION




Integrals of expensive-to-evaluate functions arise in many scientific and industrial applications, for example when expected values need to be computed and each evaluation of the integrand requires the run of a complex computer simulation where an input is only known by its distribution e.g., in meteorology, astrophysics, fluid dynamics, biology, operations research, et cetera.
This complex simulation could be a Monte Carlo simulation, a finite-element or finite-volume simulation, or a stochastic model. 
Within reasonable budget, integration using Monte Carlo may not be feasible and alternative numerical integration schemes are needed that require fewer function evaluations.


Bayesian quadrature ()—a means of constructing posterior measures over the unknown value of the integral <cit.>—mitigates a high sample demand by encoding known or assumed structure of the integrand such as smoothness or regularity, usually via a Gaussian process ().
With its increased `data'[See <cit.> and <cit.> for a discussion on `data' in numerical solvers.] efficiency,  is a natural choice when function evaluations are precious <cit.>.
In the past,  has been applied in reinforcement learning <cit.>, filtering <cit.>,
and has been extended to probabilistic integrals <cit.>.


A complementary approach to sample efficiency is to make use of related, cheaper secondary information sources.
The task of finding approximations to computationally demanding numerical models is an area of active research all on its own (see ).
Secondary information sources of reduced cost and quality include numerical models that are run at a lower resolution (a coarser grid in a fluid dynamics application), model simplifications by neglecting details or by using an approximate model that is easier to solve numerically, and analytic approximations.
A primary source could be an elaborate Earth system model to simulate anthropogenic climate change. There exist a plethora of such models and secondary sources might parametrize important effects like albedo or neglect detailed land surface processes or ocean biogeochemistry <cit.>.


Multi-source modeling is a statistical technique for harvesting information from related functions by constructing correlated surrogates over multiple sources.
When the information sources are hierarchical in that they are ordered from most to least informative, this concept is known as multi-fidelity modeling <cit.>.

The notion of multi-source models is more overarching and includes settings in which sources do not exhibit an easily identifiable order, if any.
Each of the sources has its own cost function that quantifies the cost of evaluating the source at a certain input.
An input-dependent cost might arise when the simulation run to query the integrand needs to be refined for certain values of the input to ensure numerical stability. 
A linear instance of a multi-source model is a multi-output  aka. co-kriging <cit.>.
 with multi-output s to integrate several related functions has been studied by <cit.>, who impose properties on data—which they assume given—to prove theoretical guarantees and consistency of the Bayes estimator.


 leverages active learning schemes similar to Bayesian optimization <cit.> or experimental design <cit.>.
Through its argmax, an acquisition function identifies optimal future locations to query the integrand according to a user-defined metric.
Metrics of interest in  are information gain on the value of the integral or its predictive variance.
By optimizing the target per cost, active multi-source  () trades off improvement on the target (the integral) and ressources spent.
In Bayesian optimization, this setting has been explored by <cit.>.

We summarize our contributions:


  * We lay the foundations for active  for the task of integrating an expensive function that comes with cheaper approximations. 
We assign cost to function evaluations and generalize  acquisition functions to acquisition rates that trade off improvement on the integral against cost.

  * We find that some rates induce sane, others pathological acquisition policies. Pathologies were not present in the common  acquisition schemes that all give rise to the same degenerate policy, regardless of the acquisition's value. Cost-adapted rate policies do depend on these values and are thus intricately tied to the meaning of the acquisition function that encodes progress on the quadrature task. Simply put, all considered (even pathological) multi-source acquisition policies collapse onto a single policy for , as a corner case of .

  * We conduct proof-of-concept experiments which show that  improves upon  in that it spends less budget on learning the integral to an acceptable precision.
















§ MODEL



We wish to estimate the integral over the information source of interest (the primary source), w.l.o.g. indexed by 1, f_1:Ω↦, ↦ f_1() and integrated against the probability measure π on Ω⊆^D, 

    =: ∫_Ω f_1 () π()

in presence of L-1 not necessarily ordered or orderable secondary information sources f_2, …, f_L, with f_l:Ω↦.

Each source l∈Ł = {1,…, L} comes with an input-dependent cost c_l() which must be invested to query f_l at location .
For ease of interpretation and numerical stability we set c:Ł×Ω↦ [δ, 1] and 0<δ≤ 1.
This is equivalent to assuming there exists a c_min>0 and a c_max<∞ s.t. c_min≤ c_l() ≤ c_max and then normalizing w.r.t. c_max, δ = c_min/c_max.
In other words, no query takes an infinite amount of resources, nor does any evaluation come for free.
Normalization is not required and in practice, neither c_max nor c_min need to be known.

In this section we review the tools for building a statistical model that allows us to harvest information from both the primary and the secondary sources for learning the integral  of Eq. (<ref>), before turning to the decision theoretic problem of how to actively select locations and sources to query next in Section <ref>.




 §.§ VANILLA BQ



Let f:Ω↦, ↦ f() be a function and π a probability measure on Ω⊆^D that has an intractable integral f = ∫_Ω f() π().
In vanilla Bayesian quadrature (), we express our epistemic uncertainty about the value of f through a random variable .
The distribution over  is obtained by integrating a Gaussian process () prior that is placed over the integrand f, i.e. f ∼(m, k), where m: Ω↦, ↦ m() denotes the prior mean function and k: Ω×Ω↦, (, ')↦ k(, ') the covariance function or kernel. 
Observations come in form of potentially noisy function evaluations y = f() + with ∼ (0, σ^2).
Noise free evaluations are usually assumed in , but this might not be true for a black-box integrand.
Let  denote the matrix of N input locations =[_1 … _N] and =f() + the set of corresponding observations, summarized in ={, } (see  for an introduction to  inference).  
With the closure property of s, the posterior over  when conditioning on  is a univariate Gaussian distribution with posterior mean [] = m_ and variance [] = ∫_Ω∫_Ω k_π() π(') k_ that are integrals over the 's posterior mean m_ () and covariance k_ (, '). 
These expressions are detailed below for the general multi-source case that  is a subset of and further derivations can be found in <cit.>.
So as not to replace an intractable integral by another intractable integral, the kernel k(,') is chosen to be integrable against π().




 §.§ MULTI-SOURCE MODELS



We consider linear multi-source models, which can equally be phrased as multi-output Gaussian processes <cit.> over the vector-valued function = [f_1, …, f_L], : Ω↦^L. 
Non-linear models for multi-source modeling exist and have been considered by <cit.>. 
They do however come with the additional technical difficulty that the model may not be integrable analytically—a sensible pre-requisite for —and are thus another beast altogether.
The notation mimics the single-output case, that is, ∼(, ), where  is an L× L matrix-valued covariance function.
More precisely, the covariance between two sources f_l and f_l' at inputs  and ' is [f_l(), f_l' (')] = k_ll' (, ').
The kernel k_ll'(,') encodes not only characteristics of the individual sources (smoothness), but crucially the correlation between them. 
In the multi-source setting, observations come in source-location-evaluation triplets (l, , y_l) with y_l = f_l () + _l and source-dependent observation noise _l ∼(0,σ_l^2) as usually only one element of  is being observed (see supplementary material () for alternative representation as linear observations).


The dataset = {, , _} contains N data triplets from evaluating elements of  at N locations =[_1 … _N] with corresponding sources =[l_1 … l_N] and observations _ =  [f_l_1(_1) + _l_1 … f_l_N(_N) + _l_N]. 
The  posterior over  has mean and covariance

    m_l |  ()    = m_l() + _l (, ) _() (_ - _()),
    
      k_l l' |  (, ')    = k_ll' (,') - _l (, ) _()_ l (, '),

with the kernel Gram matrix _() = _(,) + _∈^N× N and _ =  (σ_l_1^2, …, σ_l_N^2).
A summary of the notation used can be found in Table <ref> in the supplementary material.




 §.§ MULTI-SOURCE BQ



The multi-source model of Section <ref> can be integrated and gives rise to a quadrature rule similar to  (cf. sec. <ref>).
Let  denote the random variable representing the integral of interest  of Eq. (<ref>).
The posterior over  given data triplets  is a univariate Gaussian with mean and variance

    [Z ]    = m_1 + _1  (·,)_() (_ - _()),
     [Z ]     = k_11 - _1 (·,)_()_ 1 (, ·),

where _1  (·,) = ∫_Ω_1  (, ) π() is the kernel mean and k_11=∫_Ω∫_Ω k_11 (, ') π() π(') the initial error, both of source 1.
Just as in , the kernel is required to be integrable analytically.

We choose an intrinsic coregionalization model () <cit.> with kernel

    k_ll' (,') = _ll'κ(,'),

where ∈^L× L is a positive definite matrix. 
Eq. (<ref>) is a simple extension of a standard kernel κ(, ') to the multi-source case which factors the correlation between the sources and input locations. 
If κ(, ') is integrable analytically, k_ll' (,') will be, too, and thus retains the favorable property of a -kernel. 
A typical choice for κ is the squared-exponential, aka. rbf kernel κ(,') = exp (--'_2^2/2λ^2) with no dependence on the sources l and l'. 
This model can easily be extended to a linear model of coregionalization (lmc) without challenging integrability of k. 
This would untie the lenthscales between sources, but would also introduce L-1 additional generally unknown kernel parameters.
The simpler icm is also used by <cit.> to establish convergence rates for a multi-output  rule.
















§ ACTIVE LEARNING





Active learning describes the automation of the decision-making about prospective actions to be taken by an algorithm in order to achieve a certain learning objective.
A heuristic measure of improvement towards the specified goal (here: learning the value of an integral) is defined through a utility function.
It transfers the decision-theoretic problem to an optimization problem, but usually an unfeasible one.
Therefore, the utility is commonly approximated by an acquisition function.
Optimizing the acquisition function induces an acquisition policy that pins down what action to take next.
To obtain a sequence of actions, the considered method (here: ) is placed in a loop where it is iteratively fed with N_⋆ new observations (_⋆, _⋆, __⋆) in the general multi-step look-ahead (non-myopic) approach.
A myopic approximation is to instead optimize for a single new observation triplet (l_⋆, _⋆, y_l_⋆) at a time.
Besides feasibility, the lack of exact model knowledge motivates a loop in which the model is repeatedly updated with new observations.

In Section <ref>, we recapitulate several utilities that are commonly used for .
All these utilities give rise to the same acquisition policy in the absence of cost and are thus not greatly differentiated between in the literature.
Intriguingly, the policies do not coincide for  if cost is accounted for in the acquisition functions, as will be shown and discussed in Section <ref>.





 §.§ MULTI-SOURCE BQ ACQUISITIONS



In the absence of any notion of evaluation cost (or if all sources come at the same cost), the utility functions from  generalize straightforwardly to the multi-source case.
The  case can be recovered by setting the number of sources to one.



  §.§.§ Mutual Information



From an information theoretic perspective, new source-location pairs (_⋆, _⋆) can be chosen such that they jointly maximize the mutual information () I [; __⋆] between the integral  and a set of new but yet unobserved observations __⋆ with y_l_⋆^i = f_l_⋆^i (_⋆^i) + ϵ_l^i_⋆.
In terms of the individual and joint differential entropies over  and __⋆, I [; __⋆] = H[] + H[__⋆] - H[, __⋆].
Sections <ref> and <ref> imply that both  and __⋆ are normally distributed and so is their joint.
The differential entropy of a multivariate normal distribution with covariance matrix ∈^M× M is H = M/2log (2π e) + 1/2log ||.
Since there is no explicit dependence on the value of __⋆, we (sloppily) express the mutual information as a function of the new source-location pairs (_⋆, _⋆),

    I [; _⋆, _⋆]= - 1/2log(1 - ρ^2_1 _⋆ |  (_⋆) ),

where we introduce the scalar correlation

    ρ^2_1 _⋆ |  (_⋆) := _1 _⋆ |  (·, _⋆) __⋆ |  (_⋆) __⋆ 1 |  (_⋆, ·)/[],

∈ [0,1], with the noise-corrected posterior covariance matrix __⋆ |  (_⋆) = __⋆_⋆ | (_⋆, _⋆) + __⋆∈^N_⋆× N_⋆.
In the one-step look-ahead case (N_⋆=1), 

    ρ_1 l_⋆ |  (_⋆) = k_1 l_⋆ |  (·, _⋆)/√(v_l_⋆ |  (_⋆)  [])

is the bivariate correlation between  and y_l_⋆.






  §.§.§ Variance-Based Acquisitions


Variance-based approaches attempt to select (_⋆, _⋆) such that the variance on  shrinks maximally. 
As , the integral variance reduction () normalized by the current integral variance [] can be written in terms of correlation ρ as

    Δ[; _⋆, _⋆]/[]   = [] -  [∪ (_⋆, _⋆, __⋆)]/[]
        = ρ^2_1 _⋆  |  (_⋆).


Eq. (<ref>) is a monotonic transformation of Eq. (<ref>) and therefore, both utility functions share the same global maximizer _⋆.
In fact, any monotonic transformation of Eq. (<ref>), whether interpretable or not, gives rise to the same acquisition policy.
This is because the policy only depends on the locations, but not the value of the utility function's global maximum.
Hence, in  it is equivalent to consider maximal shrinkage of the variance of the integral, minimization of the integral's standard deviation, or maximal increase of the integral's precision (), to name a few—they all lead to the same active learning scheme and have thus not been greatly distinguished between in previous work on active .




 §.§ COST-SENSITIVE ACQUISITIONS





When there is a location and/or source-dependent cost associated to evaluating the information sources (cf. Section <ref>), the utility function should trade off the improvement made on the integral against the budget spent for function evaluations.
This is achieved by considering the ratio of a cost-insensitive  utility and the cost function c__⋆ (_⋆) = ∑_i=1^N_⋆ c_l_i (_i).
Such a ratio can be interpreted as an acquisition rate and bears the units of the utility function divided by units of cost.
The notion of a rate becomes clearer when considering for example the mutual information utility Eq. (<ref>) with cost measured in terms of evaluation time: the unit is bits/second, a rate of information gain.

This construction has an important consequence: Modification of the  utility function (the numerator), even by a monotonic transformation, changes the maximizer of the cost-adapted acquisition rate and hence, also the acquisition policy.
In other words, the degeneracy of  acquisition functions in terms of the policy they induce in the absence of cost is lifted when evaluation cost is included, firstly, because the argmax of each acquisition is shifted differently with cost, and, secondly, because acquisition values from different sources are discriminated against each other now.
As will be discussed below, not all monotonic transformations yield a sensible acquisition policy; indeed, some display pathological behavior.


The adapted non-myopic acquisition rates for the  utilities mutual information ( Eq. (<ref>)) and integral variance reduction ( Eq. (<ref>)) are 

    α__⋆^ (_⋆)    -log(1 - ρ^2_1 _⋆ |  (_⋆) )/c__⋆ (_⋆)
    α__⋆^ (_⋆)    ρ^2_1 _⋆ |  (_⋆)/c__⋆(_⋆),

where we have dropped the factor 1/2 in  as an arbitrary scaling factor.
It is evident that these acquisition rates do no longer share their maximizer; yet they still induce a meaningful acquisition scheme.
Both  and  have the property to be zero at ρ^2=0 and thus never select points _⋆ that are uncorrelated with the integral , no matter the cost, locations that have already been observed exactly (with σ^2=0).

Such points do not update the posterior of the integral  when conditioned on.
In  these locations are the minimizers of all acquisition functions and thus excluded no matter their value.
This is not ensured for the cost-adapted acquisition rates and therefore, they additionally require the numerator to be zero at ρ^2=0.
Hence, not every monotonic transformation of the  utility produces a sane acquisition policy in the presence of cost.

Consider for example the valid transformation ρ^2 ↦ρ^2 - 1, which is -1 at ρ=0.
Maximizing this utility function corresponds to maximizing the negative integral variance, minimizing the integral variance, which is very commonly done in .
Since ρ^2∈ [0,1], ρ^2 -1 is negative everywhere and gets larger (takes a smaller negative value) with larger cost.
Hence when maximized, this acquisition would favor expensive evaluations.

More subtle is the misbehavior of the integral precision () which is positive everywhere and has the desired behavior of favoring low-cost evaluations.
In terms of the squared correlation ρ^2∈[0, 1] from Eq. (<ref>) (with simplified notation for convenience), the numerator of the  acquisition rate can be written as (1 - ρ^2).
This expression is non-zero at ρ^2 = 0 and therefore, it does not exclude points of zero correlation when they come at sufficiently cheap cost, and in experiments we observe it getting stuck re-evaluating at the location of minimum cost ad infinitum.
We conjecture that this is because  only encodes an absolute scale of the integral variance but does not quantify any “improvement” on the integral value.





Figure <ref> illustrates the augmentation of utility functions from  with multiple information sources and cost.
Figure <ref> displays the behavior of a few acquisitions, , , and .
The right plot shows these acquisitions as used in  in terms of the squared correlation ρ^2∈[0, 1] (Eq. (<ref>)) in the absence of cost.
All acquisitions are strictly monotonically increasing functions of ρ^2.
Among the sane acquisition rates that are zero at ρ^2=0, the differences in the corresponding policy can also be understood from the functional dependence on ρ.
 diverges at perfect correlation ρ^2→ 1.
Therefore, and since the cost c lies in [δ, 1],  will always take a `perfect step' to learn the integral exactly, it will always select the points _⋆ with correlation ρ^2 (_⋆) = 1, if the step is available and no matter the cost.
, however, is finite at ρ^2=1 and trades off cost against correlation even if the perfect _⋆ with ρ^2(_⋆)=1 exists.
These interpretations are reinforced by the left three plots of Figure <ref>, in which we plot , , and  versus a univariate x_⋆ for the synthetic choice ρ^2(x_⋆)=0.95sin^2(10 x_⋆), x_⋆∈[0, 0.2] and a myopic step (N_⋆=1).
In the pure  situation, the locations of all their maxima coincide, but as soon as a non-constant cost c(x_⋆) is applied, the shapes of the acquisition functions become relevant which discriminates their _⋆ and lifts the degeneracy in policies.
 tends more towards higher correlation than , the maximizer of which moves further towards locations of lower cost.
While  and  act differently, they are both sensible choices for acquisition functions in . 
In fact for low to mid-ranged values of ρ^2  0.5 where  is approximately a linear function of ρ^2 they roughly coincide.

The choice of acquisition ultimately depends on the application and the user, who may choose which measures of improvement on the integral and cost to trade off.
















§ EXPERIMENTS









The key practical applications for  is solving integrals of expensive-to-evaluate black-box functions that are accompanied by cheaper approximations, potentially in a setting where a finite budget is available.
Typical applications are models of complex nonlinear systems that need to be tackled computationally.
With evaluations being precious, the goal is to get a decent estimate of the integral with as little budget as possible, rather than caring about floating-point precision.
In the experiments, we focus on the rear vertices of the acquisition cube Figure <ref>, multiple sources with source and input-dependent or only source-dependent cost, and separate them into two main experiments:

  
  * A synthetic multi-source setting with cost that varies in source and location for the purpose of exploring and demonstrating the behavior of the acquisition functions derived in Section <ref>.
  
  * An epidemiological model of the spread of a disease with uncertain input, in which two sources correspond to simulations that differ in cost as well as quality of the quantity of interest.

Additionally, we present a bivariate experiment with three sources in the  Section <ref>.
We take a myopic approach to all scenarios in that we optimize the acquisition for a single source-location pair a time.
The implementation of the -model uses  <cit.> in .





 §.§ MULTI-SOURCE, VARIABLE COST



We initially consider a synthetic two-source setting with univariate input.
The cost functions depend on both source and location.
The experiment's purpose is to demonstrate our findings from Section <ref> and convey intuition about the behavior of the novel acquisition functions.
The sources we consider have been suggested by <cit.> with the primary source f_1(x) = (6x - 2)^2  sin(12x - 4) and the secondary source f_2(x) = 1/2 f_1(x) + 10x for x∈[0,1].
The cost functions both take the form of a scaled and shifted logistic function in a way that the cost lies in (0,1] (cf. Figure <ref> in the ).
The costs of both sources converge to the same value close to x_⋆=0; for larger x_⋆, f_2 is two orders of magnitude cheaper than f_1.

Figure <ref> shows snapshots of three consecutive query decisions taken by the  multi-source acquisition.
The  model (depicted in the left column) has been initialized with 3 datapoints in the primary and 5 in the secondary source and merely the noise variance was constrained to 10^-2.
The  acquisition given the current state of the  is shown on the right—the top left frame is shown for , , and  in Figure <ref> in the  to emphasize the pathology of  and to highlight the subtle difference between  and  in practice.
The acquisition function is optimized using the l-bfgs-b optimizer in .
We observe that  does not query f_2 where the source costs are almost identical for x_⋆ 0.2 (see Figure <ref> in supp. mat.).
This is because the two sources are not perfectly correlated and evaluating f_1 always conveys more information about  than f_2.
The fact that c_2 decreases with increasing x_⋆ is nicely represented in the increasing height of the maxima of the dashed acquisition function for the secondary source in the top left frame of Figure <ref>.





For assessing the performance of , we compare against  and a percentile estimator (pe) that both operate on the primary source.
The latter is obtained by separating the domain into intervals that contain the same probability mass and summing up the function values at these nodes.
For the uniform integration measure used here, this is equivalent to a right Riemann sum.
We assume that  inference comes at negligible cost as compared to the function evaluations and thus consider cost to be incurred purely by querying the information sources.


To render the integration problem more difficult, we modify the Forrester functions to vary on a smaller length scale by adding a sinusoidal term and adapting some parameters, s.t.
f_1(x) = (6x - 2)^2  sin(12x - 4)  - (2-x)^2 sin(36x) and f_2(x) = 3/4 f_1(x) + 16(x - 1/2) + 10 which we integrate from 0 to 1 against a uniform measure (cf. Figure <ref>, top left).
To avoid over- or underfitting, we set a conservative gamma prior on the lengthscale with a mode at a small fraction of the domain [0,1] for both  and , and assume zero observation noise.
With six[Due to the construction of =  +  (η⃗)] more hyperparameters than ,  is more prone to over-/underfitting, and we further set a prior on the coregionalization matrix  (cf. Section <ref>) with parameters estimated from the initial three data points using empirical Bayes.
This is to avoid initial over- or under-estimation of the correlation between sources, which would either cause the active scheme to select only f_2 or only f_1, respectively.
Compared to the previous experiment, the cost is changed to have a minimum, but still composed of a sum of logistic functions and normalized to be in (0,1] (Figure <ref>, top right).
The effect of these cost functions on the final state is depicted in the , Figure <ref>.
Furthermore, this setting reveals the pathology of the  acquisition (cf. Section <ref>) that everlastingly re-evaluates the secondary source at the location of minimal cost.
The convergence behavior of the well-behaved acquisition functions  and  are displayed in Figure <ref> (bottom) in comparison to  and pe.
The hyperparameters of the  are optimized after every newly acquired node, both for  and .
Figure <ref> shows the superior performance of both  methods in arriving close the true integral with little budget.
The vertical jumps in the  methods occur when f_2 is evaluated at cheaper cost.



 §.§ A SIMULATION OF INFECTIONS





We now consider multi-source models in which sources come with input-independent cost, a.k.a. multi-fidelity models (bottom rear  vertex in Figure <ref>).
We choose an epidemiological model in which evaluating the primary source requires running numerous stochastic simulations and the secondary source solves a system of ordinary differential equations.
Epidemiological models deal with simulating the propagation of an infectious disease through a population.
The sir model forms the base for many compartmental models and assumes a population of fixed size N where at any point in time, each individuum is in one of three states—susceptible, infected, and recovered (sir)—with sizes N_S, N_I, and N_R <cit.>.
The dynamics are determined by stochastic discrete-time events of individuals changing infection state, for which Poisson processes (exponentially distributed inter-event times) are commonly assumed (see ).
In the thermodynamic limit where N is large, the average dynamics is governed by a system of odes that does not admit a generic analytic solution.
There are two parameters in the sir model: the infection rate a, and the recovery rate b.
Model details and experiment setup can be found in Section <ref> ().

For the  experiment, we assume that we know b, but we are uncertain about a. 
We are interested in the expected maximum number of simultaneously infected individuals _a[max_t N_I(t)] and the time this maximum occurs _a[_t N_I(t)], which might be relevant for vaccination planning.
Querying the primary source f_1 for the quantities of interest as a function of a requires numerous realizations of a stochastic four-compartments epidemic model (an extension to the sir model) using the Gillespie algorithm <cit.>.
For each trajectory, the maximum value and time are computed and henceforth averaged over.
In our implementation, each query of f_1 takes ∼ 16 s on a laptop's cpu.
The secondary source f_2 solves the system of s for given a and computes the maximum value and time for the resulting function N_I(t), which takes about 8·10^-3s to evaluate.
As in previous experiments, we set a gamma prior on the lengthscale, a prior on the coregionalization matrix , and the noise variance to zero as in Section <ref>.
Both  and  are given the same initial value of f_1, and  additionally gets the value of f_2 at the same location, as well as one more random datum from f_2.
This is justified since  needs to learn more hyperparameters than  and secondary source evaluations are very cheap.
Otherwise, if the initial evaluations of f_2 were further apart than the prior lengthscale from the locations of the initial primary datum, virtually zero correlation would be inferred between the sources, and the primary source would be evaluated until a sampled location roughly coincides with locations where the secondary sources have been evaluated.



Figure <ref> shows the relative error of the  estimator against normalized cost as compared to  and pe for _a[max_t N_I(t)] (left) and _a[_t N_I(t)] (right).
The horizontal dashed line shows f_2, the integral of the secondary source with one evolution of a Monte Carlo estimator of f_2.
This illustrates that simply using the secondary source for the integral estimate might be computationally cheap, but results in an unknown bias.

In the left plot,  achieves a good estimate with one additional evaluation of f_1 only, while  takes another six evaluations.
Again, the vertical jumps for  are caused by evaluations of f_2.
The initial high confidence on the integral is caused by the choice of prior on the output scale from the initial data, which is located in the tail of the gamma prior on a.
Figure <ref> displays the order in which  evaluates primary and secondary source. 





§ DISCUSSION



The multi-source model presented in Section <ref> can be extended in various ways to increase its expressiveness by using a more elaborate kernel (one lengthscale per source), or by encoding knowledge about the functions to be integrated, a probabilistic integrand. 
Other applications might come with the complication that the cost function c is unknown a priori and needs to be learned during the active -loop from measurements of the amount of resource required during the queries. 
A simple example was presented in Section <ref> where the cost was parameterized by a constant, 

estimated during the initial observations. 
A probabilistic (in contrast to parametric) model upon the cost would induce an acquisition function which is not only conditioned on the uncertain model predictions but also on the uncertain cost predictions. 
Furthermore, as in other active learning schemes, non-myopic steps for acquiring multiple observations __⋆ at once might be beneficial especially when the multi-source model is already known, and does not benefit from being re-fitted to new data;
or when multiple evaluations of sources come at lower cost than evaluating sequentially.

On the experimental side, more elaborate applications of  in areas of active research are reserved for future work.




 §.§ CONCLUSION



We have placed multi-source  in a loop and thus enabled active learning to infer the integral of a primary source while including information from cheaper secondary sources.
We discovered that utilities that yield redundant acquisition policies in  give rise to various policies, some desirable and others pathological, when evaluation cost is accounted for.
Our experiments illustrate that with the sensible acquisition functions, the  algorithm allocates budget to information retrieval more efficiently than traditional methods do for solving expensive integrals.




  §.§.§ Acknowledgements

We thank Andrei Paleyes for software-related support, as well as Philipp Hennig, Motonobu Kanagawa, and Aaron Klein for useful discussions.
AG acknowledges partial funding by the European Research Council through ERC StG Action 757275 / PANAMA and support by the International Max Planck Research School for Intelligent Systems (IMPRS-IS).













[title=References]











                                                                                             


[


empty
SUPPLEMENTARY MATERIAL 
  

  Active Multi-Information Source Bayesian Quadrature
 0.2in plus 1fil minus 0.1in


]




§ MULTI-SOURCE MODELS AND MULTI-OUTPUT GPS


We have seen in Section <ref> that linear multi-source models can be phrased in terms of multi-output s.
Typically, the goal of multi-output s is to model a vector-valued function and observations come as a vector y⃗ =  () +, where y⃗∈^L.
In multi-source models, we wish to observe only elements of , y_l = f_l() + _l.
These observations can be written as projections of the vector-valued observations,

    y_l = _ly⃗

where _l denotes a vector with a 1 in the l^th coordinate and zero elsewhere.
Let ∈^NL denote the vector of N stacked vector-valued noisy observations [y⃗_1, …, y⃗_N].
Then the corresponding N observations of elements = [l_1 … l_N] is

    _ = [ _l_1    ⋯    0;    ⋮    ⋱    ⋮;    0    ⋯ _l_N ] =: ,

where  is a sparse NL× N matrix.
Note the delicate notational difference between the N observations of single elements of , _∈^N, and a single evaluation of the vector-valued function y⃗∈^L.
The covariance matrix between all of the observations is

    [_, _] =
       (
     [ (_1, _1)        ⋯ (_1, _N);        ⋮        ⋱        ⋮; (_N, _1)        ⋯ (_N, _N) ]_ (,) + ⊗_N× N)

where =  (σ_1^2,…,σ_L^2) ∈^L× L and _N× N is an N× N matrix with every element a 1.
Also, (, ) ∈^NL× NL.
With the following mappings, we arrive at the multi-source notation introduced in Section <ref>;

    _   = _(, );
     (⊗_N× N)    = _;
       = _;   etc.

Hence, the notational detour over vector-valued observations  is not required and evaluations of individual sources can be incorporated easily in the multi-source model.
From the mappings Eq. (<ref>) follow the posterior mean and covariance of the multi-source model Eq. (<ref>).








§ ADDITIONAL PLOTS FOR SECTION <REF>


Section <ref> showed two examples to demonstrate the behavior of our derived acquisition functions.
All relevant details and cross-references are in the captions.

















§ DETAILS FOR THE INFECTIONS MODEL





 §.§ THE SIR MODEL AND EXTENSIONS

When the population size is large, the sir (susceptible, infected, recovered) model can be described by the following system of ordinary differential equations,

    N_S/ t   = -a N_S N_I/N,
     N_I/ t   = a N_S N_I/N - b N_I,
     N_R/ t   = b N_I,

in which a is the rate of infection and b the rate of recovery.
It is the most basic of a series of compartmental epidemiological models.
Various extensions exist to accommodate additional effects vital dynamics, immunity, incubation time (cf. ).
Some of these extensions serve as a general model refinement, others are relevant to specific diseases.

Statistical properties, however, are not captured by the description through odes and call for a stochastic model.
The Gillespie algorithm <cit.> enables discrete and stochastic simulations in which every trajectory is an exact sample of the solution of the `master equation' that defines a probability distribution over solutions to a stochastic equation.
In the sir model, the rate constants are time-independent and thus, the underlying process is Markovian in which the event times are Poisson distributed.
Here, an event denotes the transition of one individual from one compartment to another (N_I → N_R).



 §.§ EXPERIMENTAL SETUP




For the  experiment, we assume that we know the recovery rate b, but we are uncertain about the infection rate a. 
Therefore, we rescale the odes and place a shifted gamma prior on a/b that starts at a/b=1 and has shape and scale parameters 5 and 4 respectively.
With this prior we encode our belief that the infection rate is significantly larger than the recovery rate so an offset of the epidemic is very likely.
Also, we set the population size to N=100 to be well below the thermodynamic limit and set one individual to be infected initially.
We are interested in the expected maximum number _a[max_t N_I(t)] of simultaneously infected individuals and the time this maximum occurs _a[_t N_I(t)], which might be relevant for vaccination planning.
Querying the primary source f_1 for the quantities of interest as a function of a requires numerous realizations of a stochastic four-compartments epidemic model using the Gillespie algorithm <cit.>; in addition to the base model (sir), we include the state `exposed', in which individuals are infected but not yet infectious.
The modified system of odes that also account for assumed known incubation time γ are

    N_S/ t   = -a N_S N_I/N,
     N_E/ t   = a N_S N_I/N - γ N_E,
     N_I/ t   = γ N_E - b N_I,
     N_R/ t   = b N_I,

where we set γ = 10b.
We absorb the prior on a/b in the black-box function for all methods.

Figure <ref> shows the sir and seir models (Eq. (<ref>) and (<ref>), respectively) with 10 stochastic trajectories (thin lines). 
The solid lines indicate the mean of 100 of these stochastic realizations, and the dashed lines show the solution of the odes, in both cases for the sir model.
We also use the sir model for solving the odes even though the stochastic model simulates the seir model.
The purpose of this is to mimic a case where secondary sources are simplified simulations in that minor components are deprecated.
In the stochastic case, there is not always an outbreak of the disease, the initially infected individuum recovers before infecting someone else.
This causes the average N_R to level off significantly below 1.
For the integrals, only outbreaks are taken into account.
The corresponding integrands for the quantities of interest are shown in Figure <ref>.







§ BIVARIATE LINEAR COMBINATION OF GAUSSIANS


We construct an integrand (primary source) f_1 in the 2D-domain [-3, 3]^2 as a linear combination of K=20 normalized Gaussian basis functions 

    Φ_k^1() = (2 π |_k^1|)^-1/2e^-1/2(-_k^1) (_k^1)^-1(-_k^1),
 
i.e., f_1() = ∑_k=1^Kz_k^1Φ^1_k().
 For this, we sample K=20 means uniformly _k^1∼Uniform [-3, 3]^2 in the 2D domain. 
We then sample corresponding covariance matices _k^1 according to _̆k^1∼(0, ), κ_k^1∼Uniform [0, 1]^2, and _k^1 := (κ_k^1) + _̆k^1(_̆k^1). 
The scalar weights z_k^1 are sampled from a standard Gaussian z_k^1∼(0, 1) and can be negative. 
Thus f_1 is not a probability density function but rather a linear combination of Gaussians with varying location, shape, and weight. 
We then construct secondary sources f_2 and f_3 consecutively by adding uniform noise to the means, and additive uniform noise to the diagonal of the covariance matrices. 
Thus, with each additional source, each of the K means get randomly but consecutively shifted up and right, and the basis functions Φ_k^i(), i=2, 3 randomly become wider and flatter.
Additionally we consecutively add Gaussian random noise to the weights z_k which ensures that the true integrals of the secondary sources differ from the integral of the primary source.
All sources are depicted in Figure <ref>; the primary source f_1 on the left, and secondary sources f_2 and f_3 in the middle and right respectively. 

The cost for evaluating the primary source is 1 everywhere, the cost of evaluating f_2 and f_3 are 5% of the primary cost each. 

The priors on the kernel lengthscale and coregionalization matrix  are set analogously to the other experiments already described in Section <ref>. 
 is initialized with one evaluation of the primary source and two evaluations each of the secondary sources which amounts to a total initial cost of 1.2 (initial evaluations shown as red dots in Figure <ref>). 
Vanilla- is initialized with three evaluations which are needed to get an initial guess for its hyperparameters (initial cost=3). 
The result is shown in Figure <ref> which plots relative error of the integral estimate versus the budget spent as well as two standard deviations of the relative error as returned by the model. 
It is apparent that  finds a good solution faster than vanilla-. 





Figure <ref> illustrates the sequence of sources choses by .
Secondary source f_2 is chosen more often than secondary source f_3 at equal evolution cost of 0.05. This is intuitive since f_2, by construction, provides more information about f_1 than f_3, but both secondary sources shrink the budget equally when queried. The percentage of number of evaluations for each source after spending a total budget of 50 is 15%, 57%, 28% for sources f_1, f_2, f_3 respectively. 











