# Adaptive Reduced Multilevel Splitting

Frédéric Cérou, Patrick Héas, Mathias Rousset

► To cite this version:

Frédéric Cérou, Patrick Héas, Mathias Rousset. Adaptive Reduced Multilevel Splitting. Statistics and Computing, 2025, 36 (1), pp.13. 10.1007/s11222-025-10724-5. hal-04364553v2

HAL Id: hal-04364553

https://inria.hal.science/hal-04364553v2

Submitted on 17 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Adaptive Reduced Multilevel Splitting

Frédéric Cérou<fnref n="1" /> · Patrick Héas<fnref n="1" /> · Mathias Rousset<fnref n="1" />

Received: 24 October 2024 / Accepted: 5 September 2025
© The Author(s), under exclusive licence to Springer Science+Business Media, LLC, part of Springer Nature 2025

## Abstract

This paper considers the classical problem of sampling with Monte Carlo methods a target rare event distribution defined by a score function that is very expensive to compute. We assume we can build an approximate surrogate score using evaluations of the true score with error bounds. This work proposes a fully adaptive algorithm to sequentially sample surrogate rare event distributions with increasing target levels. An essential contribution consists in sampling at each iteration the surrogate rare event at a critical level corresponding to a specific cost. This cost is related to importance sampling of the target. The critical level is calculated solely from the reduced score and its error bound. From a practical point of view, sampling the proposal sequence is performed by extending the framework of the popular adaptive multilevel splitting algorithm to the use of score approximations. Numerical experiments evaluate the proposed importance sampling algorithm in terms of computational complexity versus squared error. In particular, we investigate the performance of the algorithm when simulating rare events related to the solution of a parametric PDE, which is approximated by a reduced basis.

**Keywords** Rare event simulation · reduced modeling · importance sampling · adaptive multilevel splitting · relative entropy

# 1 Introduction

## 1.1 Context and objective

Let $\pi(dx)$ denote a given, easily simulable, reference probability distribution on a state space $\mathcal{X}$. A random variable with distribution $\pi$ may, for instance, describe unknown physical parameters, uncertainty on initial/boundary conditions, or additional noise modeling either a chaotic behavior or the interaction with an environment. We also consider a real-valued score function $S^*: \mathcal{X} \to \mathbb{R}$ which is assumed to be of the form

$$S^*(x) = \text{score}(\Psi^*(x)),$$

where score is an easily computable function of physical interest, and $\Psi^*(x)$ is the outcome of a presumably exact numerical computation of a complex physical system

parametrized by $x \in \mathcal{X}$. With the previous ingredients, one can define the following family of conditional probability distributions on $\mathcal{X}$:

$$\eta_l^* \stackrel{\text{def}}{=} \frac{1}{p_l^*} \mathbf{1}_{S^* > l} \pi \quad l \in \mathbb{R} \qquad (1.1)$$

defined by conditioning the distribution $\pi$ on the event $\{S^* > l\}$ which may be interpreted as a certain failure region/event. The quantity $p_l^* \stackrel{\text{def}}{=} \int_{\mathcal{X}} \mathbf{1}_{S^* > l} d\pi = \pi(\{S^* > l\})$ is the associated normalization, here the probability of the considered (possibly rare) event $\{S^* > l\}$, for $l \in \mathbb{R}$. $\eta_l^*$ is the probability distribution of the system conditioned to the failure event.

We consider in this paper the classical problem of designing an algorithmic Monte Carlo procedure that samples according to $\eta_{l_{\max}}^*$ – for one specific large $l_{\max} \in \mathbb{R}$ – and estimate the associated rare event probability $p_{l_{\max}}^*$. The final event of interest is thus $\{S^* > l_{\max}\}$ while the intermediate sets are $\{S^* > l\}$ with $l < l_{\max}$.

In this work, we will focus on scenarios with two types of difficulty. A first difficulty arises when $l_{\max}$ is large, generating a very small (yet strictly positive) probability $p_{l_{\max}}^* > 0$. This rare event context confronts the practitioner with a first computational challenge: the majority of the probability mass

<footnote>Frédéric Cérou
frederic.cerou@inria.fr</footnote>

<footnote>Patrick Héas
patrick.heas@inria.fr</footnote>

<footnote>Mathias Rousset
mathias.rousset@inria.fr</footnote>

<footnote n="1">¹ Inria and Irmar, University of Rennes, Rennes, France</footnote>
of the conditional distribution $\eta_{l_{\max}}^*$ is concentrated in specific areas of $\mathcal{X}$, which are usually unknown. A direct sampling with a Monte Carlo simulation of $\pi$ requires a sample size of order at least $1/p_{l_{\max}}^*$, which is often infeasible. In that context, a popular efficient strategy to simulate a sample according to $\eta_{l_{\max}}^*$ and to estimate the associated normalizing constant is to resort to ‘Sequential Monte Carlo’ (SMC) methods. Such methods start with a Monte Carlo sample of clones (also called below ‘particles’) of the system distributed according to $\pi$ and then sequentially samples $\eta_l^*$ for increasing values of $l$. To do so, those methods resort to a combination of Importance Sampling (IS), re-sampling (selection) of weighted samples, and mutations of particles based on suitable Markov Chain Monte Carlo transitions. In this work we will focus on a specific instance called *Adaptive Multilevel Splitting* (AMS, see Cérou et al. (2012); Cérou and Guyader (2016)) which moreover chooses adaptively the sequence of levels $l$ where re-sampling is performed.

A second difficulty appears when the numerical evaluation of $\Psi^*$ is extremely expensive, so that the number of evaluations required in, say, a usual SMC or AMS approach becomes prohibitive. In order to circumvent this issue, we assume that the physical model $\Psi^*$ can be approximated by a reduced model denoted $\Psi$, the cost of evaluating $\Psi$ being small as compared to one evaluation of the true score with $\Psi^*$. The reduced model $\Psi$ we consider is assumed to come with two key features. To explain these, let us denote from now on the approximate score function by

$$S(x) = \text{score}(\Psi(x)) \in \mathbb{R}.$$

The first feature is an *a posteriori* error quantification associated with $S$ for every input parameter $x \in \mathcal{X}$. This error quantification comes in the form of an explicit error bound denoted $E(x) = \text{error}(\Psi(x)) \in \mathbb{R}$, and satisfying the pointwise estimate:

$$|S(x) - S^*(x)| \le E(x), \quad \forall x \in \mathcal{X}. \qquad (1.2)$$

Condition (1.2) is a sufficient condition that guarantees that the support of the distributions sampled by the ARMS algorithm studied here indeed contain the target rare event $\{S^* > l_{\max}\}$. This condition may not be necessary and can be weakened (see numerical Example #2b in Section 4) into an approximate bound at the price of theoretical guarantees). The second key feature assumes that we are able to update the reduced model, through a procedure of the following form

$$\{(X_1, \Psi^*(X_1)), \dots, (X_k, \Psi^*(X_k))\} \\ \xrightarrow{\text{reduced model}} \Psi = \Psi^{(X_1, \dots, X_k)},$$

where $(X_1, \dots, X_k)$ is any sequence of $k$ states in $\mathcal{X}$ – called here a sample of snapshots – together with the solution of

the associated true model $\Psi^*$ evaluated for each of those snapshots. The output is the reduced model $\Psi$. To improve the readability of the document, we will use a superscript ($k$) to specify the dependency on $(X_1, \dots, X_k)$, or even omit this dependency if the context is unambiguous,

$$(S^{(X_1, \dots, X_k)}, E^{(X_1, \dots, X_k)}) = (S^{(k)}, E^{(k)}) = (S, E).$$

## 1.2 Previous works

*Reduced models for rare events* The idea of using a simpler surrogate model in Monte Carlo simulation for target distributions is not new and there is a large amount of literature on the subject. Most references are concerned with either kriging on the one hand, and reduced basis methods on the other hand. For kriging approaches Echard et al. (2011); Moustapha et al. (2022); Razaaly et al. (2020); Zhang and Taflanidis (2018); Balesdent et al. (2013); Ling et al. (2019); Bect et al. (2017), the surrogate score function $S: \mathcal{X} \to \mathbb{R}$ is a random field, distributed according to the posterior of a prior Gaussian field conditioned by fitting the exact values of the evaluated snapshots. The estimate (1.2) is however not satisfied almost surely, and one has to resort to averaged error estimates. For reduced basis approaches Chen and Quarteroni (2013); Heinkenschloss et al. (2020); Gallimard (2019); Héas (2020), the surrogate score function is deterministic, and constructed by projection of the function $\Psi^*$ on a subspace where it is efficiently approximated and characterized with precise error bounds of the form (1.2). This model reduction framework is known to be very efficient for the numerical approximation of problems involving the repeated solution of parametric Partial Differential Equations (PDEs) Benjamin et al. (2018). Besides these two main approaches, other families of reduced models Cannamela et al. (2008); Elferson et al. (2022); Wagner et al. (2020); Peherstorfer et al. (2018) or combination of them Elferson et al. (2022); Wagner et al. (2020); Peherstorfer et al. (2018) have been proposed for rare event simulation.

Arguably, the most sensitive feature of any of these reduced model approaches is the specific collection of snapshots $(X_k, \Psi^*(X_k))$, $k \ge 1$ with which the reduced model is updated along the iterations. A natural approach consists in sampling sequentially and adaptively the collection of snapshots, the reduced model being updated after each new evaluation of the true model at iteration $k$. The updated reduced model (and the associated error estimates) can then be used to pick the next snapshot in presumably interesting areas. Sampling the snapshot $X_k$ rely on two ingredients: i) a probability distribution used to generate candidates, and ii) a weighting function, also called learning function, used to pick the presumably best candidate.

The idea of ingredient ii) can be found in the seminal work of Echard et al. (2011) where selecting snapshots is called
active learning. In the latter approach, snapshots are selected by optimizing a criterion (the learning function) that will most likely reduce the uncertainty of the surrogate model. An exhaustive set of references and a very quantitative review of methods related to active learning and kriging is provided in Moustapha et al. (2022). Some other references, as Heinkenschloss et al. (2020), pick the snapshots $(X_k)_{k \ge 1}$ uniformly among candidates using a decreasing sequence of subsets that contain the target rare event, constructed using the iteratively updated error estimates.

However, as we will argue in the present work, ingredient i) is also of high importance. Indeed, a key general phenomenon in a rare-event and high-dimensional setting is the following: the true target distribution $\eta_l^*$ for $l$ large is extremely concentrated around very specific but unknown points in $\mathcal{X}$, so that, unless the pool of candidates among which snapshots $(X_1, \dots, X_k)$ are picked is partly also located in such regions, the information given by the reduced model may eventually lack the accuracy required in order to describe the rare event. As a consequence, unless one is able to sample an extremely large amount of candidates, it is critical to select such candidates in areas that carefully track the true rare event. To make the picture clearer, consider at iteration $k$ the following surrogates of the target conditional distributions (1.1) built using the reduced score:

$$ \eta_l \stackrel{\text{def}}{=} \frac{1}{p_l} \mathbf{I}_{S>l} d\pi \quad l \in \mathbb{R}, \qquad (1.3) $$

where the normalization $p_l \stackrel{\text{def}}{=} \int_{\mathcal{X}} \mathbf{I}_{S>l} d\pi$, represents an approximation of the probability of the rare event. Then, for large $l$, the reduced rare event is also rare, hence very concentrated, and thus usually too far away from the true target rare event. If we want to guarantee sufficient mass in the region of the target rare event, we propose to specify a maximum possible level $l$ below which the true conditional distribution $\eta_l^*$ and the reduced one $\eta_l$ do match $\eta_l^* \simeq \eta_l$. This delicate issue is only rarely addressed in the literature. Most references avoid it and generate pools of candidates for the choice of snapshot relying on the reference distribution $\pi$. To the best of our knowledge, the only related pieces of work deal with the reverse problem of selecting (in a hierarchy of pre-computed surrogates) Héas (2020) or refining Bect et al. (2017) the reduced score $S$, so that it is appropriate for a certain level $l$ according to some error quantification.

The framework of the present work is to use the following guidelines to design an efficient and adaptive algorithm that synchronizes: i) the sequential Monte Carlo sampling of rare events using a reduced model (that is sampling the flow $l \mapsto \eta_l$), and ii) the sequential update of the reduced model with evaluations of the true model.

Sampling strategy Strategies for Monte Carlo sampling of the reduced conditional distributions $\eta_l$ can be grouped in two broad categories: on the one hand importance sampling (IS), usually through a parametric approach, and on the other hand splitting methods, i.e., sequential Monte Carlo Del Moral et al. (2006) methods, which are non-parametric, and in particular AMS Cérou et al. (2012) also called subset simulation Siu-Kui and Beck (2001).

Many approaches to IS of the reduced distribution $\eta_l$ have been proposed in recent years. In particular, cross-entropy methods Boer et al. (2005) have been proposed relying on kriging Balesdent et al. (2013) or reduced basis approximations Gallimard (2019); Héas (2020). In order to increase the effectiveness of the cross-entropy method, it has also been suggested to approximate the parameter variable $X$ by constraining its support in well-chosen low dimensional subspaces of $\mathcal{X}$ Uribe et al. (2021). Other recent contributions aim to build specific adaptive proposals for IS Razaaly et al. (2020); Zhang and Taflanidis (2018) or for stratified sampling Cannamela et al. (2008) (that can be interpreted as a form of IS). However, the effectiveness of IS approaches depends heavily on the choice of parametric family of proposals, for which we don't necessarily have any prior knowledge.

On the other hand, splitting methods have also been used in conjunction with (adaptively constructed) kriging-based surrogate models Ling et al. (2019); Bect et al. (2017). We refer again to the review on methods related to active learning and kriging provided in Moustapha et al. (2022). However, in these works, the splitting method samples an approximation $\eta_l$ of the target probability (1.1), which may be biased. We note nevertheless that in Elfverson et al. (2022) subset simulation is performed using nested multiresolution sets (provided by an accuracy hierarchy of the surrogate model), yielding an estimate with an associated a posteriori error.

Adapting splitting methods to reduced scores can be optionally supplemented with a bridging procedure. The goal of this bridge is to sample a new proposal $\eta_l'$ defined with some updated reduced score $S'$ and level $l'$, starting from the samples of the current proposal $\eta_l$. In general, this sampling is not trivial. For example, the nesting underlying a AMS simulation is broken, since in general updating the reduced score and the level will imply $\{S' > l'\} \not\subset \{S > l\}$. There is very little work in the literature on this bridging problem. Let us mention the work Wagner et al. (2020), which proposes a bridging procedure with a constant $l' = l$ level. Specifically, the authors use the sequential Monte Carlo method with soft levels and tempering instead of targeting distributions defined by conditioning on surrogate rare events. The authors mention that such an approach outperforms the multilevel subset simulation Ullmann and Papaioannou (2015) algorithm they previously proposed, which dealt with the nesting problem by computing correction terms and performing Markov chain Monte Carlo simulations.
## 1.3 Contributions

In this work, we propose a sequential algorithm called Adaptive Reduced Multilevel Splitting (ARMS) relying on reduced scores and error bounds. The algorithm provides an answer to the questions and limitations listed in Section 1.2.

*Summary of the ARMS algorithm.* We consider the family of proposal distributions $(\eta_l)_{l \in \mathbb{R}}$ defined in (1.3) (built using the reduced score). The goal is to use them in order to sample the final target conditional distribution $\eta_{l_{\max}}^*$. The proposed algorithm is sketched below.

**Algorithm 0 (Sketch of ARMS)**

At each iteration $k$, the algorithm performs the following tasks:

▷ Proposals sampling and search of the maximal admissible level.

1. The algorithm simulates the proposal $\eta_l^{(k)}$ for increasing values of $l$ with a Sequential Monte Carlo procedure as detailed in Section 3.1. This is done from an initial level, denoted by $l = l_b^{(k)}$, up to a maximal admissible level, denoted by $l^{(k)}$, as defined in the item below.

2. The maximal admissible level $l^{(k)}$ is defined by an admissibility constraint ensuring that the proposal $\eta_{l^{(k)}}^{(k)}$ does sufficiently match the true conditional target distribution $\eta_{l^{(k)}}^*$. The maximal admissible level $l^{(k)}$ reflects the accuracy of the current reduced model $S^{(k)}$. The matching is defined using a cross-entropy criterion and a surrogate target constructed with the error bound $E^{(k)}$ on the reduced score. The criterion and surrogate are detailed in Section 2.1 and 2.2 and their empirical approximations are given in Section 3.2.

▷ Snapshot sampling & reduced score update

3. The algorithm then samples a new snapshot of the solution of the true model $(X^{(k)}, \Psi^*(X^{(k)}))$, where $X^{(k)}$ is sampled according to a weighted version of the distribution $\eta_{l^{(k)}}^{(k)}$. The weights (or learning function) may favor states $x$ with large reduced score error $E^{(k)}(x)$. Details are given in Section 2.3.

4. It updates the reduced score $S^{(k)} \to S^{(k+1)}$ and the associated error $E^{(k)} \to E^{(k+1)}$.

▷ Bridging

5. It computes (using records of previous iterations) the next initial level $l_b^{(k+1)}$. The nesting (inclusion) of the support of $\eta_{l_b^{(k+1)}}^{(k+1)}$ in the support of a recorded past proposal must be guaranteed. The definition of $l_b^{(k+1)}$ is given in Section 2.4.

6. It bridges with a Sequential Monte Carlo simulation a record of a proposal (computed at a previous iterations) to the next proposal $\eta_{l_b^{(k+1)}}^{(k+1)}$. Details are provided in Section 3.3.

The Monte Carlo aspect of the algorithm is implemented using an adaptive Sequential Monte Carlo strategy called Adaptive Multilevel Splitting (AMS), similar to the approach referred to as Nested Sampling. As in any Sequential Monte Carlo simulation methodology, the flow $\eta_l^{(k)}$ for increasing values of $l$ is approximated by the empirical distribution of a sample (population) of particles (clones).

The proposed method can then yield importance sampling estimators, see details in Section 2.6. Consider $H$ snapshots sampled starting from a random iteration index $K_0$ (defined in Section 2.6) without re-weighting: the snapshots being sampled according to the proposal

$$X^{k+1} \sim \eta_{l^{(k)}}^{(k)},$$

it is then natural to consider the following estimator:

$$\hat{\pi}^H \stackrel{\text{def}}{=} \frac{1}{H} \sum_{k'=K_0+1}^{K_0+H} \frac{d\pi}{d\eta_{l^{(k'-1)}}^{(k'-1)}} (X^{(k')}) \delta_{X^{(k')}}, \quad (1.4)$$

which by construction, is an *unbiased* empirical estimator of the reference distribution $\pi$ but *restricted on the intersection of the support of the proposals* (see Section 2.6 for a proof). Note that unbiasedness is true up to the adaptivity features in sampling the proposal $\eta_l$, which is done in practice using AMS, possibly introducing a bias expected to be of order $1/N$ where $N$ is the number of Monte Carlo clones. Note also that the necessity to restrict the estimator $\hat{\pi}^H$ on the common support of the proposals requires to check in principle throughout the algorithm that these supports do contain the target rare event $\{S^* \ge l_{\max}\}$, using the error bound (1.2). An IS estimator e.g. of the rare event probability $p_l^*$ is then

$$\hat{p}_l = \hat{\pi}^H (\mathbf{1}_{S^* > l}).$$

A idealized (that is without the Monte Carlo aspect detailed) pseudo-code is given in Section 2.7.

*First methodological contribution: Item 2 in Algorithm 0.* Our first methodological contribution consists in performing at each iteration $k$ snapshot sampling (using the proposal $\eta_{l^{(k)}}^{(k)}$) at a maximal admissible level $l^{(k)}$. We use a specific criterion to adaptively compute the maximum admissible level. At each iteration $k$, we will use the information given by the error quantification (1.2) ensure that the proposal $\eta_{l^{(k)}}^{(k)}$ remains a reasonable importance distribution for sampling the true target $\eta_{l^{(k)}}^*$. The evaluation of the efficiency of an importance proposal distribution will be done using a quantification of the logarithmic cost of importance sampling. This quantification is defined by the relative entropy or Kullback-Leibler divergence between probabilities defined
as

$$ 
\begin{aligned}
\ln \text{cost} &\approx \text{Ent}(\eta_{l^{(k)}}^* | \eta_{l^{(k)}}^{(k)}) \stackrel{\text{def}}{=} \int \ln \frac{d\eta_{l^{(k)}}^*}{d\eta_{l^{(k)}}^{(k)}} d\eta_{l^{(k)}}^* \\
&= -\ln \int \mathbf{1}_{\{S^* > l^{(k)}\}} d\eta_{l^{(k)}}^{(k)},
\end{aligned}
$$

when $\eta_{l^{(k)}}^*$ is dominated by $\eta_{l^{(k)}}^{(k)}$ and extended to $+\infty$ otherwise. In the present rare event context, this is simply minus the log probability of the current true rare event $\{S^* > l^{(k)}\}$ under the proposal distribution $\eta_{l^{(k)}}^{(k)}$. This choice is discussed and justified in Section 2.1 and Appendix A. The latter cannot be evaluated exactly, but we can leverage the existence of error estimates in order to obtain an approximation of it denoted $\tilde{\eta}_{l^{(k)}}^{(k)}$. Relying on this cost, we are able to trigger snapshot sampling when the reduced modeling becomes critically inaccurate, the associated (maximum) level being $l^{(k)}$. Snapshot sampling is then performed using the distribution $\eta_{l^{(k)}}^{(k)}$, which we may reweight by favoring states with higher error quantification, as mentioned above.

**Second methodological contribution: Items 3 and 4 in Algorithm 0.** First remark that a simple way to proceed after the update of the reduced model is simply to perform a fresh AMS simulation starting from the initial reference distribution $\pi$. in order to significantly reduce the number of evaluations of the reduced model, our second methodological contribution is to propose instead a *bridging* procedure. At iteration $k$, this procedure searches in the simulation history for a well-chosen previous proposal denoted $\eta_{l^{(k_b)}}^{(k_b)}$ with $k_b \le k$, so that the updated proposal denoted $\eta_{l^{(k+1)}}^{(k+1)}$, with $l^{(k+1)}$ a well-designed level, is dominated and sufficiently close to the former. The latter proposal is then simulated starting from the former one using an AMS framework.

**Numerical results.** In numerical experiments, we study the proposed algorithm, called *Adaptive Reduced Multilevel Splitting* (ARMS) and show the empirical convergence of the ARMS algorithm on a toy model and a rare-event problem based on a realistic elliptic PDE and a *reduced basis* approach for model reduction. Our study confirms that the proposed relative entropy criterion plays a crucial role in ARMS and prevent the simulation to be taken into spurious regions of the state-space remote from the region of interest, resulting in a large variance of the IS estimator. Moreover, by evaluating the computational complexity in relation to the squared error of the estimate, we show empirically that the bridging procedure reduces in a certain regime the simulation cost while achieving a comparable squared error.

The paper is organized as follows. The main methodological concepts underlying the ARMS algorithm are presented in Section 2. We end this section by introducing an idealized

version of the ARMS algorithm, which assumes the exact evaluation of expectations and exact sampling of distributions. In Section 3, we define ARMS, which consists in a practical implementation of this idealized algorithm using AMS simulations. In Section 4, we numerically evaluate the proposed algorithm for various rare event problems. A final section deals with conclusions and perspectives. The appendix collects some technical proofs and details.

# 2 Key Elements and Idealized Algorithm

## 2.1 Importance sampling cost

Let $l_{\max}$ denote a large level of interest, and let us assume for the sake of the argument within the present section, that one is able to simulate a random variable $X$ distributed according to the proposal distribution $X \sim \eta_l = \frac{1}{p_l} \mathbf{1}_{S^* > l} d\pi$, the normalization $p_l$ being known. *Importance sampling* of the target rare event $\{S^* > l\}$ can be performed in our context if the distribution $\eta_l$ dominates the measure $\mathbf{1}_{S^* > l} d\pi$, and simply amounts to remark that one can estimate the averages

$$ p_l \mathbb{E}[\mathbf{1}_{S^*(X)>l} \varphi(X)] = \pi(\varphi \mathbf{1}_{S^*>l}), $$

where $\varphi$ denotes a generic test function. For instance taking $\varphi = \mathbf{1}$, the above yields an unbiased estimator of the normalization probability $p_l^*$.

The domination requirement, denoted $\eta_l \gg \eta_l^*$, holds true by definition if and only if one has the inclusion $\{S^* > l\} \subseteq \{S > l\}$. This is not guaranteed here as $S$ is an approximation of the true score $S^*$, and since $S^*$ is not computed for each sample states, one cannot check it in a systematic fashion. Moreover, even if the above domination condition holds, importance sampling becomes quickly unrealistic if the target distribution diverges too widely from the proposal, a typical phenomenon in high dimension. As discussed in Appendix A, it has been thoroughly argued in Chatterjee and Diaconis (2018) (see also Cérou et al. (2022)) that the cost of importance sampling – in terms of the required sample size $K$– is quite generically and roughly given by the exponential of the relative entropy between $\eta_l^*$ and $\eta_l$:

$$ K \approx e^{\text{Ent}(\eta_l^* | \eta_l)}. \qquad (2.1) $$

This fact will justify the choice of various heuristics that will be found in this paper. We will call $\text{Ent}(\eta | \mu)$ the *log cost of importance sampling* of the target $\eta$ by the proposal $\mu$. In the present rare event case, note that under the condition that $\eta_l \gg \eta_l^*$, relative entropy is simply the logarithm of small
probability ratios:

$$ \mathrm{Ent}(\eta_l^* | \eta_l) = \ln \frac{p_l}{p_l^*} = -\ln \int \mathbf{1}_{S^{*>l}} d\eta_l. $$

The domination between measures is implied by finite relative entropy, and we will use throughout the paper the condition $\mathrm{Ent}(\eta | \mu) < +\infty$ in order to check the domination condition $\mu \gg \eta$ for any measures $\mu$ and $\eta$.

## 2.2 Maximal admissible level

In this section, we describe how to choose the target level $l^{(k)}$, which is suited to the current score function $S^{(k)}$ at iteration $k$ of the algorithm $-k$ being the number of sampled snapshots of the true model which have been used to enrich the initial score approximation $S^{(0)}$. The idea is to tune $l^{(k)}$ with a criterion evaluating that surpassing this level should require, in some sense, a better surrogate model, and thus new snapshots. The adjusted proposal $\eta_{l^{(k)}}^{(k)}$ will then be used to sample a new snapshot, and obtain an updated score function $S^{(k+1)}$ for the next iteration $k + 1$. In short, we are going to develop a heuristic criterion on levels to detect when a new snapshot is needed.

The general insight enabling to define $l^{(k)}$ is the following. Consider the following problem: one is given a family of rare event target distributions $(\eta_l^*)_{l \in \mathbb{R}}$ in the form of (1.1), a given computational budget, and a family of importance sampling proposal distributions $(\eta_l^{(k)})_{l \in \mathbb{R}}$ in the form of (1.3). One wants to perform some sequential algorithm to sample the distributions $\eta_l$ along an increasing sequence of levels $l$ as high as possible, and then use importance sampling. However, we are given limited resources preventing sampling for levels that are unreasonably high. As an idealized criterion formalizing admissible levels, we propose the importance sampling cost (2.1) of the true target. The (ideal) maximal admissible (hence critical) level is thus given by the smallest level for which this cost matches a certain given numerical threshold $c$, that is, $\delta > 0$ being a very small numerical level increment:

$$ \text{Find first } l \text{ : } \mathrm{Ent}(\eta_{l+\delta}^* | \eta_{l+\delta}^{(k)}) > c. $$

By doing so, we ensure that the snapshot sampled with respect to $\eta_{l^{(k)}}^{(k)}$ will not diverge too much from the target distribution $\eta_{l^{(k)}}^*$ while increasing levels as much as possible. As an example, if, as it is the case in the present context, $\eta_l^{(k)}$ has been constructed using a reduced model based on $k$ snapshots $(X_1, \dots, X_k)$, the maximal admissible level $l^{(k)}$ will typically be of order $\max(S^*(X_1), \dots, S^*(X_k))$, as the quality of the approximation $S^* - S^{(k)}$ will deteriorate in areas with higher scores.

In a practical context, the log cost $\mathrm{Ent}(\eta_l^* | \eta_l^{(k)})$ has to be approximated using some prior information on $\eta_l^*$, for instance using some form of error quantification. We propose in this work the following method, which is by no means unique. We consider a pessimistic (respectively an optimistic) surrogate of $\eta_l^*$, denoted by $\tilde{\eta}_l^{(k)}$ (respectively by $\hat{\eta}_l^{(k)}$), constructed with the reduced score $S^{(k)}$ and its error quantification $E^{(k)}$, and typically associated with a much lower (respectively higher) rare event probability than $\eta_l^*$. In the present work, we will choose (2.3)-(2.4) below, which satisfy:

$$ \tilde{\eta}_l^{(k)} \ll \eta_l^* \ll \hat{\eta}_l^{(k)}. $$

The maximal admissible level will be estimated using the pessimistic $\tilde{\eta}_l$ as a surrogate for the true target. Given an initial level $l_b$ and a very small level increment $\delta > 0$, this yields the following rigorous definition of the maximal admissible level:

$$ \begin{aligned} l^{(k)}(c, l_b) &= \min\{l \ge l_b : \mathrm{Ent}(\tilde{\eta}_{l+\delta}^{(k)} | \eta_{l+\delta}^{(k)}) > c \text{ or} \\ & \qquad \mathrm{Ent}(\hat{\eta}_{l_{\max}}^{(k)} | \eta_{l+\delta}^{(k)}) = +\infty\}. \end{aligned} \quad (2.2) $$

For convenience, the dependency on $c$ and $l_b$ of $l^{(k)}$ will be removed if the context makes it clear. Note that the additional constraint $\mathrm{Ent}(\tilde{\eta}_{l_{\max}}^{(k)} | \eta_l^{(k)}) < +\infty$ in (2.2) ensures that the support of $\eta_l^{(k)}$ for $l < l^{(k)}(c, l_b)$ contains the support of the final target $\eta_{l_{\max}}^*$.

In the present rare event case, the importance distribution is (1.3), the point-wise error bound is (1.2) and we propose the pessimistic and optimistic surrogate of the target (1.1)<fnref n="1" />

$$ \tilde{\eta}_l^{(k)} \stackrel{\text{def}}{=} \frac{1}{\tilde{p}_l^{(k)}} \mathbf{1}_{S^{(k)} - E^{(k)} > l} \pi \quad \text{with} \quad \tilde{p}_l^{(k)} \stackrel{\text{def}}{=} \pi(\mathbf{1}_{S^{(k)} - E^{(k)} > l}), \quad (2.3) $$

and

$$ \hat{\eta}_l^{(k)} \stackrel{\text{def}}{=} \frac{1}{\hat{p}_l^{(k)}} \mathbf{1}_{S^{(k)} + E^{(k)} > l} \pi \quad \text{with} \quad \hat{p}_l^{(k)} \stackrel{\text{def}}{=} \pi(\mathbf{1}_{S^{(k)} + E^{(k)} > l}). \quad (2.4) $$

**Remark 2.1** Using the surrogate (2.3), the particularization of the worst-case log cost in (2.2) is $\mathrm{Ent}(\tilde{\eta}_l^{(k)} | \eta_l^{(k)}) = \ln \frac{p_l^{(k)}}{\tilde{p}_l^{(k)}}$, that is to say the logarithm of the ratio of the approximate probability $\pi(\{S^{(k)} > l\})$ over the smaller pessimistic probability $\pi(\{S^{(k)} - E^{(k)} > l\})$.

<footnote n="1">¹ $\mu(f)$ denotes $\int f(x)\mu(dx)$ so that $\mu f/\mu(f)$ denotes the normalized measure.</footnote>
**Remark 2.2** It can easily be checked that the critical level $l^{(k)}(c, l_b)$ tends to infinity when the error $E^{(k)}$ vanishes pointwise.

## 2.3 Snapshot sampling and learning function

Given the importance distribution $\eta_{l^{(k)}}^{(k)}$, we now discuss the procedure used to sample a new snapshot. A possible choice used in the present work is the penalization of the importance distribution as follows:

$$X \sim \frac{1}{\eta_{l^{(k)}}^{(k)} (e^{\tau^{(k)} E^{(k)}}) e^{\tau^{(k)} E^{(k)}} \eta_{l^{(k)}}^{(k)}},$$

with $\tau^{(k)} \ge 0$. The above penalization is similar to the so-called learning function in the literature on 'active learning' Moustapha et al. (2022). In this work, we simply choose an exponential function of the error $c$ with a penalty parameter $\tau^{(k)}$ adaptively chosen at iteration $k$ of the algorithm. In particular, $\tau^{(k)}$ will be chosen to be high in the beginning of the algorithm and then be reduced to 0 at some point. The idea here is that $\tau^{(k)}$ enables to find a compromise between two objectives: first i) the snapshots are used to refine the reduced model $S^{(k)}$, and ii) the snapshots are also used in estimation by importance sampling. A large $\tau^{(k)}$ is associated with the objective i). Indeed, favoring candidates with high errors in a $\eta_{l^{(k)}}^{(k)}$-distributed sample enables a faster learning of the true score $S^*$. A small or vanishing $\tau^{(k)}$ is associated with objective ii). Indeed, the distribution $\eta_{l^{(k)}}^{(k)}$ is constructed to be a robust importance sampling distribution for the target $\eta_{l^{(k)}}^*$ so that states with smaller error are needed for estimation, ensuring a limited variance for importance sampling. In other words, keeping the values of $\tau^{(k)}$ too large will lead to a non-robust importance sampling estimate, while setting them too low from the start will slow the learning of the true score and incur an additional cost. In practice, there is a wide range of values that allow for reasonably accurate learning of the score while ensuring importance sampling robustness.

In the numerical experiments of Section 4, we will consider a very large $\tau^{(k)}$ in a first phase of the algorithm (until a given number of snapshots are taken in the rare event $\{S^* > l_{\max}\}$ of interest), and then, because estimation is then possible, we will take $\tau^{(k)} = 0$.

No attempt is made to optimize the learning function and the choice of the penalty parameter; those studies are left for future work, taking into account the studies on active learning methods as found in the survey Moustapha et al. (2022).

## 2.4 Bridging to distributions with updated score approximation

We will now focus on what happens after the update of the reduced model that follows the sampling of a new snapshot. We are given a certain level $l^{(k)}$ and an associated importance distribution $\eta_{l^{(k)}}^{(k)}$. The problem is to find a new initial level which we will denote $l_b(k+1)$ and then sample the new importance distribution $\eta_{l_b(k+1)}^{(k+1)}$ associated with the updated score approximation $S^{(k+1)}$.

We emphasize that this bridging methodology is only optional. It is perfectly possible (and we will provide some numerical comparisons) to set $l_b(k+1) = -\infty$ which amounts to start at the next iteration with an importance distribution $\eta_{l_b(k+1)}^{(k+1)} = \pi$ given by the reference distribution.

This initial level $l_b(k+1)$ will then be used in the subsequent $k+1$ iteration to determine (with definition (2.2)) the critical level at iteration $k+1$: $l^{(k+1)}(c, l_b(k+1))$ for $S^{(k+1)}$. This will complete the main loop of the algorithm.

However, if one tries to sample $\eta_{l_b(k+1)}^{(k+1)}$ starting from $\eta_{l^{(k)}}^{(k)}$, using some form of importance or sequential sampling, there is an issue with inclusion. For instance, simply setting the initial level $l_b(k+1)$ to $l^{(k)}$, the inclusion $\{S^{(k+1)} > l^{(k)}\} \subseteq \{S^{(k)} > l^{(k)}\}$ will generally not be valid, making the task infeasible.

In order to circumvent this issue, we propose to keep in memory the importance distributions $\eta_{l^{(k')}}^{(k')}$ for $k' \le k$. We then set a quantile $\theta = \frac{M}{N} \in (0, 1)$ associated with an effective sample size $N - M \in \mathbb{N}$, and find an index

Find: $k_b \le k$

and a level

Find: $l_b(k+1)$

verifying:

$$\text{Ent} \left( \eta_{l_b(k+1)}^{(k+1)} | \eta_{l^{(k)}}^{(k_b)} \right) \le -\ln(1-\theta). \qquad (2.5)$$

This condition ensures that a relatively straightforward procedure enables to sample the distribution $\eta_{l_b(k+1)}^{(k+1)}$ starting from samples of the previous importance distribution $\eta_{l^{(k_b)}}^{(k_b)}$. In particular, the inclusion condition $\{S^{(k+1)} \ge l_b(k+1)\} \subseteq \{S^{(k_b)} > l^{(k_b)}\}$ will hold true.

In addition to the condition (2.5), we also want to ensure that the two conditions in (2.2) used to define the maximal admissible level $l^{(k+1)}$ hold for this initial level $l_b(k+1)$. This is necessary to ensure formally that $l_b(k+1) \le l^{(k+1)}$.
Explicitly:

$$ \mathrm{Ent}(\tilde{\eta}_{l_b(k+1)}^{(k+1)} | \eta_{l_b(k+1)}^{(k+1)}) \le c, \qquad (2.6) $$

$$ \mathrm{Ent}(\hat{\eta}_{l_{\max}}^{(k+1)} | \eta_{l_b(k+1)}^{(k+1)}) < +\infty. \qquad (2.7) $$

A pair $(k_b, l_b(k+1)) \in [0, k] \times \mathbb{R}$ satisfying (2.5), (2.6) and (2.7) determine two consecutive proposal distributions $\eta_{l_b(k_b)}^{(k_b)}$ and $\eta_{l_b(k+1)}^{(k+1)}$ “bridging” the score approximations $S^{(k_b)}$ and $S^{(k+1)}$. This pair can be determined according to various strategies.

* A trivial feasible pair is to set $k_b = 0$ such that $l^{(k_b)} = l_b(k+1) = -\infty$. As explained in Section 3, this will account in practice to restart the AMS simulation with the new score $S^{(k+1)}$ from the beginning until it reaches the critical level $l^{(k+1)}(c, -\infty)$. However, this option can lead to a significant computational overhead if the complexity of the reduced score evaluation is not negligible.

* A non-trivial feasible pair of interest is given by the largest (lexicographic) pair $(k_b, l_b(k+1)) \in [0, k] \times \mathbb{R}$ satisfying the three conditions (2.5), (2.6) and (2.7). This means first selecting the importance sampling distribution $\eta_{l_b(k_b)}^{(k_b)}$, among the already computed distributions, from the most recent to the oldest, and then to search the largest level $l_b(k+1)$ such that the three conditions hold. We will see in Section 3 that in practice, this solution will lead to an important computation saving, but will require nevertheless to save the history of the previous AMS simulations.

## 2.5 Stopping criteria

The maximum level $l_{\max}$ is setting the range of the computation of interest so that the level $l^{(k)}(c, l_b)$ defined in (2.2) is constrained to stay below $l_{\max}$. The practical definition of $l^{(k)}(c, l_b)$ is thus modified as the minimum between $l_{\max}$ and the critical level given by (2.2). The algorithm can thus continue as long as one wishes, for instance until enough snapshot samples are obtained.

In addition, in order to avoid unnecessarily enriching an already accurate reduced score and obtain computational savings, it is possible to stop updating the reduced score if the worst-case log cost at the critical level $l^{(k)}$ is almost zero and if this critical level is $l_{\max}$. More precisely we define the update stopping index<fnref n="2" /> as the first index $k$ such that

$$ \mathrm{Ent}(\tilde{\eta}_{l^{(k)}}^{(k)} | \eta_{l^{(k)}}^{(k)}) < \epsilon \quad \& \quad l^{(k)} \simeq l_{\max}, \qquad (2.8) $$

where $\epsilon \ge 0$ is a given threshold parameter close to the machine precision. Stopping the main ARMS loop with this criteria must not however be implemented when using a non-trivial bridging procedure. We will indeed see in Section 3 that by doing so, the normalization constant is no longer updated resulting in a sustained variance while new snapshots are sampled. Yet this sustained variance always occurs when using a non-trivial bridging procedure. We will see in Section 3 that using the latter stop-update rule with a non-trivial bridge usually largely compensate additional variance with computational (and memory) savings.

## 2.6 The estimator

We are now ready to discuss the estimation of the probability $p_{l_{\max}}^*$. It is more convenient to work with the un-normalized version of the target distribution $\eta_{l_{\max}}^*$ defined in (1.1), which will be denoted from now on by $\gamma_{l_{\max}}^*$ and defined by

$$ \gamma_{l_{\max}}^*(\varphi) \stackrel{\mathrm{def}}{=} \pi(\varphi \mathbf{1}_{S^* > l_{\max}}), $$

for the level of interest $l_{\max} \in \mathbb{R}$ and some test function $\varphi$. Note that $p_{l_{\max}}^* = \gamma_{l_{\max}}^*(\mathbf{1})$.

The importance sampling estimator $\hat{\gamma}_{l_{\max}}^H$ of $\gamma_{l_{\max}}^*$ is defined using a random iteration index which we will denote by $K_0$. Only the snapshots sampled at and after iteration $K_0+1$ will be used to build our importance sampling estimator. More precisely, $K_0$ must be of optional type, in the sense that $K_0$ is defined as the first iteration index $k$ for which a certain property depending on the snapshots $X_1, \dots, X_k$ holds true.

We propose the following option: $K_0$ is the first iteration index $k$ such that the score of the past snapshots have reached the level $l_{\max}$ (the rare event) $j_0$ times, for some prescribed $j_0 \ge 1$. At iteration $K_0 + H$ (for some $H > 0$) of the algorithm, the estimator of $\gamma_{l_{\max}}^*$ is then given by:

$$ \hat{\gamma}_{l_{\max}}^H \stackrel{\mathrm{def}}{=} \frac{1}{H} \sum_{k'=K_0+1}^{K_0+H} f_{l_{\max}}^{(k')}(X_{k'+1}) \delta_{X_{k'+1}}, \qquad (2.9) $$

where, if we choose for $k' \ge K_0$ to set the learning parameter as $\tau = 0$, the IS density ratio can be simplified as

$$ f_{l_{\max}}^{(k')}(x) \stackrel{\mathrm{def}}{=} \mathbf{1}_{S^* > l_{\max}}(x) \frac{d\pi}{d\eta_{l^{(k')}}^{(k')}}(x). $$

Setting a learning parameter $\tau > 0$ in the early stages of the algorithm is crucial as the proposal law has a low probability mass over the set of rare events, which implies that snapshots taken so far are still insufficient to learn an accurate reduced model in this rare event regime. However, once $j_0$ snapshots have already been sampled from the set of rare events, the

<footnote n="2">² We call a stopping index any exit criterion of the algorithm that is independent of the future snapshots.</footnote>
choice of $\tau = 0$ implies that the snapshot $X_{k'+1}$ is drawn according to the proposition $d\eta_{l(k')}^{(k')}$, thus favoring the robustness of the IS estimator rather than the accuracy of an already refined reduced model. Note that : i) the entropic criterion in (2.2) and (2.7) ensures that $\eta_{l_{\max}}^* \ll \eta_{l(k)}^{(k)}$, and ii) $f_{l_{\max}}^{(k)}(x)$ can be evaluated point-wise as long as we have an estimate of the normalizations constant of the proposal $\eta_{l(k)}^{(k)}$.

As shown in Appendix B with a martingale argument, because of the optional property of $K_0$, this estimator is unbiased

$$ \mathbb{E} \hat{\gamma}_{l_{\max}}^{H} = \gamma_{l_{\max}}^{*}, $$

and its variance decreases in $O(H^{-1})$. Note that samples $X_k$ in (2.9) are assumed to be drawn exactly from the sequence of importance distributions $\eta_{l(k-1)}^{(k-1)}$. In practice, this assumption will not hold exactly because of the adaptive features used in the sampling of the proposals by the Monte Carlo routine, yielding a small bias of order the inverse of the sample size. This issue is discussed in Section 3.

## 2.7 The idealized importance sampling algorithm

Based on the concepts introduced in Sections 2.1 to 2.6, we can now describe in this section the algorithm that will lead to an estimation of $\gamma_{l_{\max}}^*$ based on importance sampling.

### Algorithm 1 (Idealized ARMS)

**Require:** target level $l_{\max}$, initial reduced score $S^{(1)}$ and error $E^{(1)}$ functions, budget of snapshots $K$, worst-case log cost threshold $c$, reduced score precision $\epsilon$, trigger parameter $j_0$, quantile $\theta \in (0, 1)$ for the bridge between proposals, boolean restart variable for restarting or bridging, learning intensity $\tau$

▷ (initialization)

1: Set $l^{(0)} = l_b(1) = -\infty$, $H = 0$, $j = 0$, $k = 0$, $p_l^{(1)} = 1$

2: while $K > 0$ do

3: $k \leftarrow k + 1$

▷ (searching the maximal admissible level, see Section 2.2)

4: Simulate the flow of proposal distributions $\eta_l^{(k)}$ for $l \le l^{(k)}$. The maximal level $l^{(k)}$ must be a level $l \in [l_b(k), l_{\max}]$ satisfying both of the following two conditions:

$$ \text{Ent}(\tilde{\eta}_l^{(k)} | \eta_l^{(k)}) \le c $$

$$ \text{Ent}(\hat{\eta}_{l_{\max}}^{(k)} | \eta_l^{(k)}) < +\infty. $$

$l^{(k)}$ is defined as the smallest $l \in [l_b(k), l_{\max}]$ (in practice $l^{(k)} = l - \delta$ for some $\delta > 0$ very small) violating one of the above two conditions; if search fails, then $l^{(k)} = l_{\max}$.

5: Compute 3 the associated normalisation $p_l^{(k)}$.

▷ (snapshot sampling & score update, see Section 2.3)

6: if $j > j_0$ then draw $X \sim \eta_{l(k)}^{(k)}$, else $X \sim \frac{e^{\tau E^{(k)}}}{\eta_{l(k)}^{(k)}(e^{\tau E^{(k)}})} \eta_{l(k)}^{(k)}$ end if

7: Compute the snapshot $S^*(X) = \text{score}(\Psi^*(X))$

8: if $\text{Ent}(\tilde{\eta}_{l(k)}^{(k)} | \eta_{l(k)}^{(k)}) > \epsilon$ or $l^{(k)} < l_{\max}$ then

9: Update reduced model: build $(S^{(k+1)}, E^{(k+1)})$ by enriching the previous surrogate model $(S^{(k)}, E^{(k)})$ with $(X, S^*(X))$.

10: else

11: Don't update the reduced model: set $(S^{(k+1)}, E^{(k+1)}) \leftarrow (S^{(k)}, E^{(k)})$.

12: end if

13: if $S^*(X) > l_{\max}$ then $j \leftarrow j + 1$ end if
▷ (importance sampling estimation, see Section 2.6)

14: if $j = j_0$ then $K_0 = k$ end if

15: if $k > K_0$ then

16: $H = k - K_0$

17:

$$ \hat{\gamma}_{l_{\max}}^{H} \leftarrow \frac{1}{H} \left( (H-1)\hat{\gamma}_{l_{\max}}^{H-1} + p_{l(k)}^{(k)} \mathbf{1}_{S^*(X)>l_{\max}} \delta X \right) $$

18: end if

19: Set $K \leftarrow K - 1$

20: if $(S^{(k+1)}, E^{(k+1)}) \neq (S^{(k)}, E^{(k)})$ then
▷ (bridging, see Section 2.4)

21: if restart then

22: Find the largest (lexicographic) pair $(k', l') \in [0, k] \times \mathbb{R}$ satisfying the following three conditions:

$$ \text{Ent}(\eta_{l'}^{(k+1)} | \eta_{l(k')}^{(k')}) \le -\ln(1 - \theta) $$

$$ \text{Ent}(\tilde{\eta}_{l'}^{(k+1)} | \eta_{l'}^{(k+1)}) \le c, $$

$$ \text{Ent}(\hat{\eta}_{l_{\max}}^{(k+1)} | \eta_{l'}^{(k+1)}) < +\infty, $$

and denote the solution pair by $(k_b(k+1), l_b(k+1)) \leftarrow (k', l')$.

23: Simulate the proposal $\eta_{l_b(k+1)}^{(k+1)}$ using $\eta_{l_b(k+1)}^{(k+1)}$; then compute the associated normalization $p_{l_b(k+1)}^{(k+1)}$.

24: else

25: Set $l_b(k+1) = -\infty$

26: Set $p^{(k+1)} = 1$ and reinitialize $\eta_{l_b(k+1)}^{(k+1)}$ with $\pi$

27: end if

28: else

29: Set $l_b(k+1) = l^{(k)}$

30: Set $p^{(k+1)} = p^{(k)}$ and $\eta_{l_b(k+1)}^{(k+1)} = \eta_{l(k)}^{(k)}$

31: end if

32: end while

33: Return Estimate $\hat{\gamma}_{l_{\max}}^{H}$

This algorithm is an idealized concept and, in practice, it is not directly computable. Indeed it assumes exact evaluation of expectations and exact sampling according to the proposal and the surrogate distributions which are defined up to a normalization constant.

## 3 Practical Implementation: AMS Simulations

In this section, we propose a practical implementation of the idealized algorithm introduced in Section 2.7 using an AMS methodology. The latter will simulate the flow of proposals for growing level parameters $l$ by providing Monte Carlo samples with approximately correct distributions for large sample sizes.

<footnote n="3">3 With a sequential Monte Carlo approach for instance, see next sections.</footnote>
## 3.1 AMS

We wish to draw a sample of size $N$ according to the importance sampling distribution $\eta_{l_i}^{(k)}$ in the form (1.3), where, for ease of presentation, we temporarily assume that the target level $l^{(k)}$ is known. We suppose we start with a given sample of size $N$ approximately distributed according to the law $\eta_{l_b}^{(k)}$ with $l_b < l^{(k)}$, as well as with a given estimator of the associated normalization $p_{l_b}^{(k)}$. The AMS algorithm Cérou et al. (2012) will be used here to generate a finite sequence of $N$-samples approximately distributed according to a sequence of nested distributions $\{\eta_{l_i}^{(k)}\}_{i \ge 1}$, as well as estimators of the associated normalizations $\{p_{l_i}^{(k)}\}_{i \ge 1}$, i.e., such that $l_b < l_1 < \dots < l_m = l^{(k)}$, where $m$ is the number of intermediate levels. At iteration $i \ge 0$, the $N$-sample of particles (clones) will be denoted

$$ (\Xi_1^{(i)}, \dots, \Xi_N^{(i)}) $$

and we omit in notation the dependence in $(k)$. The associated empirical distribution will be denoted

$$ \eta_{l_i}^{N,(k)} := \frac{1}{N} \sum_{n=1}^{N} \delta_{\Xi_n^{(i)}} \approx \eta_{l_i}^{(k)}, $$

and the estimator of the normalization constant $p_{l_m}^{(k)}$ will be denoted $p_{l_m}^{N,(k)}$. As a consequence, AMS will yield an approximation of the reference distribution $\pi$ restricted to the current reduced score $S^{(k)}$ as follows:

$$ p_{l_m}^{N,(k)} \eta_{l_m}^{N,(k)} (\varphi \mathbf{1}_{S^{(k)} > l^{(k)}}) \approx \pi (\varphi \mathbf{1}_{S^{(k)} > l^{(k)}}), $$

for $\varphi$ any given test function.

Let us describe briefly the AMS algorithm. At the $i$-th iteration, the algorithm performs successively: i) a weighting and selection step, and ii) a mutation step. The weighting step $i$) assigns weights to the empirical distribution $\eta_{l_i}^{N,(k)}$ so that it approximates $\eta_{l_{i+1}}^{(k)}$, up to a normalization constant given by an empirical approximation of the normalization ratios. More precisely, the empirical distribution at level $l_{i+1}$ and the associated normalizations are obtained by weighting the current $N$-sample:

$$ \frac{p_{l_{i+1}}^{N,(k)}}{p_{l_i}^{N,(k)}} \eta_{l_{i+1}}^{N,(k)} := \frac{1}{N} \sum_{j=1}^{N} \mathbf{1}_{S^{(k)} > l_{i+1}} \delta_{\Xi_j^{(i)}}, $$

where the equality is well-defined because $l_i < l_{i+1}$. This empirical distribution of course requires the determination (or a choice) of $l_{i+1}$. AMS performs an adaptive choice of $l_{i+1}$ as follows. AMS starts by sorting the $N$ particles according

to their scores in ascending order. AMS then considers an adaptive target distribution defined by a fixed number $M < N$ of particles. The latter are assigned a weight equal to zero, and the level $l_{i+1}$ is defined by the $M$-th smallest particle score

$$ l_{i+1} = S^{(k)}(\Xi_{(M)}^{(i)}), $$

in which $(1), \dots, (N)$ is the usual shorthand notation for increasing order statistics (based on scores, here). In other words, this corresponds to set $l_{i+1}$ as the highest level $l$ in a set given by a maximal relative entropy

$$ \begin{aligned} & l_{i+1} \\ := & \sup \left\{ l \in [l_i, l_{\max}] : \begin{array}{l} \text{Ent}(\eta_l^{(k)} | \eta_{l_i}^{(k)}) \approx -\ln \eta_{l_i}^{N,(k)}(\mathbf{1}_{S^{(k)} > l}) \\ < -\ln \left( \frac{N-M}{N} \right) \end{array} \right\}, \end{aligned} $$

with the quantile (proportion of killed particles) satisfies $M = \lfloor \theta N \rfloor$.

In the selection step, the first $M$ particles with zero-weights are removed and surviving particles are duplicated. The $M$ resampled particles are chosen uniformly among the surviving particles. The AMS approximation of the normalization constant thus follows the recursion<fnref n="4" />

$$ p_{l_{i+1}}^{N,(k)} = p_{l_i}^{N,(k)} \eta_{l_i}^{N,(k)} (\mathbf{1}_{S^{(k)} > l_{i+1}}) = p_{l_i}^{N,(k)} (1 - M/N), \quad (3.1) $$

yielding

$$ p_{l_i}^{N,(k)} = p_{l_m}^{N,(k)} = p_{l_b}^{N,(k)} (1 - M/N)^m. $$

The mutation step ii) then consists of a Markovian transition that leaves $\eta_{l_{i+1}}^{(k)}$ invariant and is applied individually and independently a prescribed number of times to the $M$ particles after resampling, the goal being to bring diversity in the $N$-sample to counterbalance the effect of selection. In practice, Markovian transitions can be performed using a Metropolis-Hastings algorithm or other techniques studied in the Markov Chain Monte Carlo (MCMC) literature. We mention that the Markovian kernel is often chosen to be adaptive in order to maintain a given Metropolis acceptance rate Garthwaite et al. (2016).

<footnote n="4">⁴ To ease the presentation we have assumed here that particles have different scores. However, in the case $\exists \mathcal{J} \subset \mathbb{N}$ such that $S^{(k)}(\Xi_{(M)}^{(i)}) = S^{(k)}(\Xi_{(M+j)}^{(i)})$ for $j \in \mathcal{J}$, a good practice is to kill in addition the set of particles $\{S^{(k)}(\Xi_{(M+j)}^{(i)})\}_{j \in \mathcal{J}}$, and locally change the value of $M$ in (3.1) to $M + \text{card}(\mathcal{J})$ accordingly.</footnote>
Note that AMS typically involves many evaluations of the score function underlying the sequence of targeted distributions. This paper focuses precisely on cases where the true score $S^*$, which is too expensive to evaluate many times, is substituted by a reduced score $S^{(k)}$, which is much cheaper to evaluate.

## 3.2 Reaching the critical level with AMS

Let us now describe the main Monte Carlo sampling routine of the algorithm. We recall that after $k$ iterations, one is given with a current reduced score $S^{(k)}$, an associated point-wise error estimate $E^{(k)}$, an initial level $l_b(k)$ and a sample of size $N$ approximately distributed according to $\eta_{l_b(k)}^{(k)}$, that we denote using a sample-size superscript: $\eta_{l_b(k)}^{N,(k)}$. One want to perform an AMS simulation targeting the distribution $\eta_{l_b(k)}^{(k)}$, where we recall that $l^{(k)}$ is the maximal admissible level above $l_b(k)$. The latter is formally defined by (2.2) (for $N = +\infty$). When using sample-based empirical distributions, one will need to adapt the latter conditions.

The AMS simulation will proceed as explained in Section 3.1, creating a sequence $(l_i)_{i \ge 0}$ of increasing levels starting from $l_b = l_b(k)$, creating a sequence of empirical distributions $\eta_{l_i}^{N,(k)}$ for $i \ge 0$. We then check for each proposed level increase from $l_i$ to $l_{i+1}$ that the entropic conditions (2.2) using the empirical distribution $\eta_{l_i}^{N,(k)}$ to approximate the latter. More precisely, we remark that

$$ \begin{aligned} \operatorname{Ent}(\tilde{\eta}_{l_{i+1}}^{(k)} | \eta_{l_{i+1}}^{(k)}) &= \ln \frac{\eta_{l_i}^{(k)}(\mathbf{1}_{S^{(k)} > l_{i+1}})}{\eta_{l_i}^{(k)}(\mathbf{1}_{S^{(k)} - E^{(k)} > l_{i+1}})} \\ &\approx \ln \frac{\eta_{l_i}^{N,(k)}(\mathbf{1}_{S^{(k)} > l_{i+1}})}{\eta_{l_i}^{N,(k)}(\mathbf{1}_{S^{(k)} - E^{(k)} > l_{i+1}})}, \end{aligned} \quad (3.2) $$

$$ \begin{aligned} \operatorname{Ent}(\hat{\eta}_{l_{\max}}^{(k)} | \eta_{l_{i+1}}^{(k)}) &= \ln \frac{\eta_{l_i}^{(k)}(\mathbf{1}_{S^{(k)} > l_{i+1}})}{\eta_{l_i}^{(k)}(\mathbf{1}_{S^{(k)} + E^{(k)} > l_{\max}})} \\ &\approx \ln \frac{\eta_{l_i}^{N,(k)}(\mathbf{1}_{S^{(k)} > l_{i+1}})}{\eta_{l_i}^{N,(k)}(\mathbf{1}_{S^{(k)} + E^{(k)} > l_{\max}})}, \end{aligned} \quad (3.3) $$

so that (2.2) is re-written as:

$$ \begin{aligned} l^{(k)} = l_m := & \min\{l_i, i \ge 0 : \ln \frac{\eta_{l_i}^{N,(k)}(\mathbf{1}_{S^{(k)} > l_{i+1}})}{\eta_{l_i}^{N,(k)}(\mathbf{1}_{S^{(k)} - E^{(k)} > l_{i+1}})} \\ & \ge c \text{ or } \ln \frac{\eta_{l_i}^{N,(k)}(\mathbf{1}_{S^{(k)} > l_{i+1}})}{\eta_{l_i}^{N,(k)}(\mathbf{1}_{S^{(k)} + E^{(k)} > l_{\max}})} = +\infty \}. \end{aligned} $$

The above means that when the AMS simulation proposes to increment the level $l_i$ to $l_{i+1}$ and that at least one of these

conditions is violated, then the critical level set to $l^{(k)} = l_i$ is reached.

**Remark 3.1** If the logarithm of the ratio in (3.3) is finite when $N \to \infty$ for all intermediate levels $l_j$ from $-\infty$ to $l_i$ then we can verify recursively that, at a formal level,

$$ \begin{aligned} & \lim_{N \to \infty} \ln \frac{\eta_{l_i}^{N,(k)}(\mathbf{1}_{S^{(k)} > l_{i+1}})}{\eta_{l_i}^{N,(k)}(\mathbf{1}_{S^{(k)} + E^{(k)} > l_{\max}})} < +\infty \\ & \Leftrightarrow \operatorname{Ent}(\hat{\eta}_{l_{\max}}^{(k)} | \eta_{l_{i+1}}^{(k)}) < +\infty, \end{aligned} $$

which precisely means that the rare event of interest is a subset of the support of the current proposal distribution $\eta_{l_{i+1}}^{(k)}$.

The procedure to reach the critical level in practice is summarized by the following algorithm.

**Algorithm 2 (AMS simulation up to the critical level)**

**Require:** level of interest $l_{\max}$, initial distribution $\eta_{l_b(k)}^{N,(k)}$ and normalization $p_{l_b(k)}^{N,(k)}$, related reduced score $S^{(k)}$ and error $E^{(k)}$, quantile $\theta \in [1/N, 1[$, worst-case log cost threshold $c$.

1:

2: Set $i = 0, c_{try} = 0$

3: **while** $l_i < l_{\max}$ & $c_{try} \le c$ & $d_{try} < +\infty$ **do**

4: Order scores $S^{(k)}$ relative to particles of $\eta_{l_i}^{N,(k)}$

5: Set $l_{i+1} = \min(l_M, l_{\max})$ with $l_M$ the M-th smallest score

6: Compute

$$ c_{try} = \ln \frac{\eta_{l_i}^{N,(k)}(\mathbf{1}_{S^{(k)} > l_{i+1}})}{\eta_{l_i}^{N,(k)}(\mathbf{1}_{S^{(k)} - E^{(k)} > l_{i+1}})} \approx \operatorname{Ent}(\tilde{\eta}_{l_{i+1}}^{(k)} | \eta_{l_{i+1}}^{(k)}) $$

as well as

$$ d_{try} = \ln \frac{\eta_{l_i}^{N,(k)}(\mathbf{1}_{S^{(k)} > l_{i+1}})}{\eta_{l_i}^{N,(k)}(\mathbf{1}_{S^{(k)} + E^{(k)} > l_{i+1}})} \approx \operatorname{Ent}(\hat{\eta}_{l_{\max}}^{(k)} | \eta_{l_{i+1}}^{(k)}) $$

as defined in (3.2) and (3.3)

7: **if** $l_{i+1} < l_{\max}$ & $c_{try} \le c$ & $d_{try} < +\infty$ **then**

8: Kill the $M \ge \lfloor \theta N \rfloor$ particles with the lowest score or with score equal to $l_{i+1}$

9: Update $p_{l_{i+1}}^{N,(k)}$ with recursion (3.1)

10: Perform $M$ duplications (uniformly among the survivors)

11: Mutate the duplicates using a MCMC leaving target $\eta_{l_{i+1}}^{(k)}$ invariant.

12: Denote the resulting sample $\eta_{l_{i+1}}^{N,(k)}$.

13: **else**

14: Set $l^{(k)} = l_i$

15: **end if**

16: Set $i \leftarrow i + 1$

17: **end while**

18: **Return** Empirical measure $\eta_{l_{i(k)}}^{N,(k)}$ and normalizing constant $p_{l_{i(k)}}^{N,(k)}$

The following remark is useful in practice, as it provides some prior information on the worst-case log cost threshold parameter $c$.
**Remark 3.2** (Worst-case log cost: a lower bound in practice) Consider, in Algorithm 2 above, the empirical distribution $\eta_{l_i}^{N,(k)}$. By construction of AMS, $l_{i+1}$ is the M-th smallest score among the N scores. Set $c_{\min}(M, N) \stackrel{\text{def}}{=} \ln \frac{N-M}{N-M-1}$. As long as $E^{(k)}$ is non zero, at least one particle of the sample, the M-th particle, lies outside the level set obtained by subtracting the error to the score. $c_{\min}(M, N)$ is thus the lower bound of the entropy between the empirical measures with and without the error term, that is

$$c_{\min}(M, N) \le \ln \frac{\eta_{l_i}^{N,(k)}(\mathbf{I}_{S^{(k)}>l_{i+1}})}{\eta_{l_i}^{N,(k)}(\mathbf{I}_{S^{(k)}-E^{(k)}>l_{i+1}})}.$$

Setting the numerical threshold to the minimal possible value $c = c_{\min}(M, N)$ can be considered a safe strategy, leading to a slower but safer convergence of the algorithm.

### 3.3 Bridging with an AMS step

As exposed in Section 2.4, the bridging procedure at the end of iteration $k$ consist in selecting a past index $= k_b := k_b(k) \le k$ and a new initial level $l_b(k+1)$ which will be used in the AMS simulation with the new score approximation $S^{(k+1)}$ starting from a past (appropriately recorded) empirical distribution $\eta_{l_i}^{N,(k_b)}$. As defined previously, bridging consists in choosing, using a lexicographic search, a feasible pair $(k_b, l_b(k+1)) \in [0, k] \times \mathbb{R}$, where feasibility means satisfying the quantile condition (2.5) and the two entropic conditions (2.6) and (2.7). It is now necessary to define empirical counterparts to the latter conditions. As previously, the two entropic conditions are approximated by

$$\text{Ent}(\tilde{\eta}_{l'}^{(k+1)} | \eta_{l'}^{(k+1)}) \approx \ln \frac{\eta_{l}^{N,(k')} (\mathbf{I}_{S^{(k+1)}>l'})}{\eta_{l}^{N,(k')} (\mathbf{I}_{S^{(k+1)}-E^{(k)}>l'})} \le c, \quad (3.4)$$

$$\text{Ent}(\hat{\eta}_{l_{\max}}^{(k+1)} | \eta_{l'}^{(k+1)}) \approx \ln \frac{\eta_{l}^{N,(k')} (\mathbf{I}_{S^{(k+1)}>l'})}{\eta_{l}^{N,(k')} (\mathbf{I}_{S^{(k+1)}+E^{(k)}>l_{\max}})} < +\infty, \quad (3.5)$$

where in the above $(k', l')$ denote feasible pairs that are candidates for the solution $(k_b, l_b(k+1))$. Note that the well-posedness of (3.5) can be justified by Remark 3.1. Concerning the quantile condition, we remark that, given some $k'$ and $l'$, and assuming that the domination relation stands

$$\eta_{l}^{(k')} \gg \eta_{l'}^{(k+1)} \Leftrightarrow \{S^{(k+1)} > l'\} \subset \{S^{(k')} > l^{(k')}\}, \quad (3.6)$$

an empirical approximation of the quantile condition (2.5) is

$$\eta_{l}^{N,(k')} (\mathbf{I}_{S^{(k+1)} \le l'}) \le \theta. \quad (3.7)$$

In order to verify (3.6), we propose to browse the history of past AMS simulations and check the inclusion of the support of $\eta_{l'}^{(k+1)}$ in the one of $\eta_{l}^{(k')}$ using all the recorded N-samples targeting the past proposal distributions $\eta_{l}^{N,(\tilde{k})}$ with $\tilde{k} < k_b$, whose samples are scattered between the reference distribution and the current distribution. This check must be carried out for all $\tilde{k} < k_b$ and is given by

$$\eta_{l(\tilde{k})}^{N,(\tilde{k})} (\mathbf{I}_{S^{(k')} > l^{(k')}} \mathbf{I}_{S^{(k+1)} > l'} + \mathbf{I}_{S^{(k+1)} \le l'}) = 1. \quad (3.8)$$

To sum up, bridging consist in selecting some feasible pair, i.e., some pair $(k_b, l_b(k+1))$ satisfying (3.4), (3.5), (3.7) and (3.8). Note that once the pair has been determined, the normalization estimate can be updated as

$$p_{l_b(k+1)}^{N,(k+1)} = p_{l_b(k_b)}^{N,(k_b)} \eta_{l^{(k_b)}}^{N,(k_b)} (\mathbf{I}_{S^{(k+1)}>l_b(k+1)}). \quad (3.9)$$

Then the $\tilde{M}(k_b)$ particles (at max given by $\tilde{M}(k_b) = \lfloor \theta N \rfloor$ when (3.7) is saturated) with updated scores below or equal to $l_b(k+1)$ are replaced by duplicating surviving particles and performing mutation targeting $\eta_{l_b(k+1)}^{(k+1)}$ (in an AMS line of thinking). The next AMS simulation with $S^{(k+1)}$ can then be started from the level $l_b(k+1)$. The following iteration then starts up to the critical level $l^{(k+1)}(c, l_b(k+1))$, an so on.

As mentioned in Section 2.4, among the feasible set, we focus on either the trivial pair $(k_b, l_b(k+1)) = (0, -\infty)$, or the largest (lexicographic) pair $(k_b, l_b(k+1)) \in [0, k] \times \mathbb{R}$. A simple way to solve this maximisation problem is to begin with the highest possible index $k' := k$ and then to decrease $k'$ by a unit if no feasible level $l'$ has been found for this index. To search for the largest feasible level $l'$ at a given index $k'$, one can proceed as follows. Begin by setting $l'$ to the highest level satisfying (3.7), then lower this level until (3.4) and (3.5) hold and finally perform for all $\tilde{k} < k_b$ the check (3.8).

The non-trivial bridging procedure is summarized by the following algorithm.

**Algorithm 3 (Bridging)**

**Require:** level of interest $l_{\max}$, current index $k$, for all $k' \le k$, particles of distribution $\eta_{l^{(k')}}^{N,(k')}$, normalization $p_{l^{(k')}}^{N,(k')}$ and related reduced score $S^{(k')}$ and error $E^{(k')}$, reduced score $S^{(k+1)}$, the quantile $\theta \in [1/N, 1[$, worst-case log cost threshold $c$.

1:

2: Set $k' = k$ and feasible = false
    ▷ (find largest feasible pair)

3: while !feasible & $k' > 0$ do

4: Compute $S^{k+1}$ for particles in $\eta_{l^{(k')}}^{(k')}$ and sort their scores

5: Set $\tilde{M} = \lfloor N\theta \rfloor$

6: while !feasible & $\tilde{M} > 0$ do

7: Set $l' = \min(l_{\tilde{M}}, l_{\max})$ with $l_{\tilde{M}}$ the $\tilde{M}$-th smallest score
8: Compute (3.4): $c_{try} = \ln \frac{\eta_{l(k')}^{N,(k')} (1_{S(k+1)>l'})}{\eta_{l(k')}^{N,(k')} (1_{S(k+1)-E(k+1)>l'})} \approx \text{Ent}(\tilde{\eta}_{l'}^{(k+1)} | \eta_{l'}^{(k+1)})$

9: Compute (3.5): $d_{try} = \ln \frac{\eta_{l(k')}^{N,(k')} (1_{S(k+1)>l'})}{\eta_{l(k')}^{N,(k')} (1_{S(k+1)+E(k+1)>l_{max}})} \approx \text{Ent}(\hat{\eta}_{l_{max}}^{(k+1)} | \eta_{l'}^{(k+1)})$

10: if $c_{try} \le c \& d_{try} < +\infty$ then

11: Set feasible = true

12: for $\tilde{k} \in [0, k' - 1]$ do

13: Compute $S^{(k+1)}$ and $S^{(k)}$ for particles in $\eta_{l(k)}^{(\tilde{k})}$

14: if $\eta_{l(k)}^{N,(\tilde{k})} (1_{S^{(k)>l(k)}} \mathbf{1}_{S^{(k+1)>l'}} + \mathbf{1}_{S^{(k+1)}<l'}) \ne 1$

15: Set feasible = false

16: end if

17: end for

18: Set $\tilde{M} \leftarrow \tilde{M} - 1$

19: end while

20: Set $k' \leftarrow k' - 1$

21: end while

22: if feasible then

▷ ($k' > 0$: perform an AMS step)

23: Set the largest pair as $(k_b, l_b(k+1)) = (k', l')$

24: Compute $p_{l_b(k+1)}^{N,(k+1)}$ using (3.9)

25: Kill the $\tilde{M}$ particles in $\eta_{l(k_b)}^{N,(k_b)}$ of smallest score $S^{(k+1)}$ or with score equal to $l_b(k+1)$

26: Perform duplications (uniformly among the survivors) to replace the killed particles

27: Simulate $\eta_{l_b(k+1)}^{N,(k+1)}$ by mutating the duplicates to target $\eta_{l_b(k+1)}^{(k+1)}$

28: else

▷ ($k' = 0$: initialize a fresh AMS)

29: Set $l_b(k+1) = -\infty$

30: Set $p_{l_b(k+1)}^{N,(k+1)} = 1$

31: Simulate $\eta_{l_b(k+1)}^{N,(k+1)}$ (by sampling the reference distribution $\pi$)

32: end if

33: Return Empirical measure $\eta_{l_b(k+1)}^{N,(k+1)}$ and normalizing constant $p_{l_b(k+1)}^{N,(k+1)}$

## 3.4 AMS-based importance sampling estimator

Using the sequence of empirical distributions $\{\eta_{l(k)}^{N,(k)}\}_{k=1}^K$ generated by AMS simulations, we can write a closed-form approximation of the importance sampling estimator $\hat{\gamma}_{l_{max}}^H$ of $\gamma_{l_{max}}^*$ defined by (2.9). The AMS approximation of the estimator $\hat{\gamma}_{l_{max}}^H$ is defined by

$$ \hat{\gamma}_{l_{max}}^{H,N}(\varphi) \stackrel{\text{def}}{=} \frac{1}{H} \sum_{k=K_0+1}^{K_0+H} p_{l(k)}^{N,(k)} \varphi(X_{k+1}) \mathbf{1}_{S^*>l_{max}}(X_{k+1}), \quad (3.10) $$

In particular, the AMS approximation of the estimator $\hat{\gamma}_{l_{max}}^H(\mathbf{1})$ of $p^* = \gamma_{l_{max}}^*(\mathbf{1})$ is

$$ \hat{p}_{IS} \stackrel{\text{def}}{=} \hat{\gamma}_{l_{max}}^{H,N}(\mathbf{1}), \quad (3.11) $$

AMS also provides as by-products a sequence of approximate estimators of $\pi(\{S^{(k)} \ge l_{max}\})$ given by

$$ \hat{p}_{AMS} \stackrel{\text{def}}{=} p_{l(k)}^{N,(k)} \eta_{l(k)}^{N,(k)} (\mathbf{1}_{S^{(k)}>l_{max}}), \quad (3.12) $$

This estimator can be an increasingly accurate approximation of $p^*$, but as our numerical experiments show, it can also be highly biased, even in the case where the reduced score $S^{(k)}$ converges uniformly to $S^*$ as $k$ increases.

## 3.5 The ARMS algorithm

The complete methodology proposed in this work is summarized in the complete algorithm presented below named adaptive reduced multilevel splitting (ARMS). ARMS is a practical implementation of the idealized importance sampling estimator (Algorithm 1), built on the AMS simulation framework described above.

### Algorithm 4 (ARMS)

**Require:** target level $l_{max}$, initial reduced score $S^{(1)}$ and error $E^{(1)}$, budget of snapshots $K$, the quantile $\theta \in (0, 1)$ proposing the bridge level and level increases in AMS, trigger parameter $j_0$, number of particles $N$, worst-case log cost threshold $c$, reduced score precision $\epsilon$, boolean restart variable for restarting or bridging, learning intensity $\tau$.

▷ (initialization)

1: Set $l^{(0)} = l_b(1) = -\infty$, $H = 0$, $j = 0$, $k = 0$, $p_{l^{(0)}}^{N,(0)} = 1$

2: Set $\eta_{l_b(1)}^{N,(1)}$ by drawing a $N$-sample of the reference distribution $\pi$

3: while $K > 0$ do

4: $k \leftarrow k + 1$

5: ▷ (AMS simulation to critical level)

6: Use Algorithm 2 (see Section 3.2) from $l_b = l_b(k)$ to:

• estimate the critical level $l^{(k)} \in [l_b(k), l_{max}]$

• simulate $\eta_{l(k)}^{N,(k)}$ starting from $\eta_{l_b(k)}^{N,(k)}$

• compute the associated normalisation $p_{l(k)}^{N,(k)}$.

▷ (snapshot sampling & score update)

6: if $j > j_0$ then draw $X \sim \eta_{l(k)}^{N,(k)}$, else $X \sim \frac{\eta_{l(k)}^{N,(k)} e^{\tau E^{(k)}}}{\eta_{l(k)}^{N,(k)} (e^{\tau E^{(k)}})}$ end if

7: Compute the snapshot $S^*(X)$

8: if $l^{(k)} < l_{max}$ or Ent($\hat{\gamma}_{l(k)}^{(k)} | \eta_{l(k)}^{(k)}$) $\simeq$ $-\ln \eta_{l(k)}^{(k)} (\mathbf{1}_{S^{(k)}-E^{(k)}>l^{(k)}}) > \epsilon$ then

9: Build $(S^{(k+1)}, E^{(k+1)})$ by enriching $(S^{(k)}, E^{(k)})$ with $(X, S^*(X))$.

10: else

11: Set $(S^{(k+1)}, E^{(k+1)}) \leftarrow (S^{(k)}, E^{(k)})$.

12: end if

13: if $S^*(X) > l_{max}$ then $j \leftarrow j + 1$ end if

▷ (importance sampling)

14: if $j = j_0$ then $K_0 = k$ end if

15: if $k > K_0$ then

16: $H = k - K_0$

17: Update estimate (see Section 3.4)

$$ \hat{\gamma}_{l_{\max}}^{H,N} \leftarrow \frac{1}{H} \left( (H-1) \hat{\gamma}_{l_{\max}}^{H-1,N} + p_{l(k)}^{N,(k)} \mathbf{1}_{S^*(X)>l_{\max}} \delta_X \right) $$

18: end if

19: Set $K \leftarrow K - 1$

20: if $(S^{(k+1)}, E^{(k+1)}) \ne (S^{(k)}, E^{(k)})$ then
21: **if restart then**

▷ (AMS step for bridging)

22: Use Algorithm 3 (see Section 3.3) to:

• find the largest feasible pair $(k_b, l_b(k+1))$

• simulate $\eta_{l_b(k+1)}^{N,(k+1)}$ from $\eta_{l_b(k)}^{N,(k_b)}$

• compute the associated normalisation $p_{l_b(k+1)}^{N,(k+1)}$.

23: **else**

▷ (or AMS reinitialization)

24: Set $l_b(k+1) = -\infty$

25: Set $p_{l_b(k+1)}^{N,(k+1)} = 1$

26: Set $\eta_{l_b(k+1)}^{N,(k+1)}$ by drawing a N-sample of the reference distribution $\pi$

27: **end if**

28: **else**

29: Set $l_b(k+1) = l^{(k)}$

30: Set $p_{l_b(k+1)}^{N,(k+1)} = p_{l_b(k)}^{N,(k)}$, $\eta_{l_b(k+1)}^{N,(k+1)} = \eta_{l_b(k)}^{N,(k)}$

31: **end if**

32: **end while**

33: **Return** Estimate $\hat{\gamma}_{\max}^{H,N}$

We comment hereafter the ARMS algorithm. The simulation starts from an i.i.d. sample with the reference distribution. ARMS then iterates through four blocks until the computing budget allocated to snapshots is exhausted.

The first block (lines 4 to 5) consists in increasing the level until reaching the maximal admissible level presented in Section 2.2 and implemented in practice by an AMS routine in Section 3.2. In short, it simulates the distribution flow $l \mapsto \eta_l^{N,(k)}$ with $l \in [l_b(k), l^{(k)}]$, where $l^{(k)}$ is the critical level where IS of the true target $\eta_l^*$ by the current reduced model $S^{(k)}$ is expected to be too expensive.

The second block (lines 6 to 14) samples a new sample according to the proposal $\eta_l^{N,(k)}$, penalized by some learning function in the early iterations. It then computes the related snapshot, updates the reduced score and finally updates the learning function. This block follows the guideline presented in Section 2.3.

ing procedure searches in past simulations to find a certain pair $(k_b, l_b(k+1))$ which enables to simulate from $\eta_{l_b(k)}^{N,(k_b)}$ towards $\eta_{l_b(k+1)}^{N,(k+1)}$ in a single AMS step. Otherwise (in the case restart = true), a fresh AMS simulation is initialized at each score approximation update. We mention that the bridging routine is avoided in the case where the stopping criterion (2.8) is reached. Indeed, in this case, no bridging is needed, since the reduced score is not updated anymore and the next initial level can simply be set to $l^{(k)}$.

## 4 Numerical Simulations

We first consider a toy example given by an explicit one dimensional true score. Reduced modeling is performed by cubic spline interpolation. The rare event simulation is in this case one-dimensional but nonetheless challenging, as the score presents several maxima and moreover splines may produce highly oscillating approximations. We then focus our attention on a more realistic scenario in which the rare event is defined through the high-dimensional solution of an elliptic PDE with random inputs. The PDE solution is approximated by a reduced basis approach which yields the reduced models. All our numerical simulations are carried out using an implementation of the ARMS algorithm in the Julia language, available at https://gitlab.inria.fr/pheas/arms.

### 4.1 Examples

We will consider the problem of estimation of $p^* = \pi(\mathbf{1}_{S^* > l_{\max}})$ with $S^*(X) = \text{score}(\Psi^*(X))$, for various models $\Psi^*$ and reduced models $\Psi$, various dimensions $q = \text{dim}(X)$ of the parameter space, various model dimensions $h = \text{dim}(\Psi^*(X))$ and different score functions.

<table><thead><tr><th></th><th>Ψ<sup>*</sup></th><th>Ψ</th><th>dim(X)</th><th>dim(Ψ<sup>*</sup>(X))</th><th>score</th><th>p<sup>*</sup></th><th>l<sub>max</sub></th></tr></thead><tbody><tr><td>Example #1</td><td>analytical</td><td>spline</td><td>1</td><td>1</td><td>|·|</td><td>2.18e-08</td><td>90</td></tr><tr><td>Example #2a</td><td>PDE</td><td>RB</td><td>4</td><td>5101</td><td>L<sup>1</sup></td><td>6.31e-04</td><td>0.5</td></tr><tr><td>Example #2b</td><td>PDE</td><td>RB</td><td>4</td><td>5101</td><td>L<sup>∞</sup></td><td>1.47e-04</td><td>3.0</td></tr><tr><td>Example #2c</td><td>PDE</td><td>RB</td><td>25</td><td>5101</td><td>L<sup>1</sup></td><td>2.48e-14</td><td>0.5</td></tr></tbody></table>

The third block (lines 15 to 19) is the estimation phase beginning when the algorithm has reached the rare event in some sense. This block performs importance sampling with the proposal $\eta_l^{N,(k)}$ and updates the IS estimate (3.10).

The fourth and final block (lines 20 to 31) is the restart or bridge routine presented in Section 2.4 and implemented in practice by an AMS step in Section 3.3. In short, if the restart option is not triggered (restart = false), then the bridg-

We detail hereafter these examples.

#### 4.1.1 Example #1: spline approximation of a closed-form score

In a first numerical experiment, we design a one-dimensional analytic toy model. Cubic spline data interpolation is used to build the surrogate $S(x)$ from the snapshots and the error
<figure><img src="image_3.png" /><figcaption>Fig. 1 Illustration of the algorithm adaptive behavior (Example #1). The four plots illustrate the algorithm at iteration k = 1, 20, 40, 60. They display in blue the score function S* (the rare event being</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>S*(X) ≥ 90), in red the reduced score approximation S at iteration k, in green the snapshots (X₁, ..., Xₖ) used to build the reduced score</figcaption></figure>

bound is set to two times the absolute value of the true error $E(x) = 2|S(x) - S^*(x)|$.

Using spline interpolation is not recommended in practice, as it easily creates wide oscillations that can artificially create spurious high score areas that mislead the sampling routines. It is nonetheless an interesting benchmark for difficult situations where the reduced modeling produces very sub-optimal approximations in some areas.

The score function, depicted in Figure 1, is defined for $x \in \mathcal{X} = \mathbb{R}_+$ by $S^*(x) = |\Psi^*(x)|$ with

$$ \Psi^*(x) = \begin{cases} l_{\max} & \text{if } x \le 1/l_{\max}, \\ 1/x + f(x) & \text{else} \end{cases}, $$

where $f(x) = \alpha(\mathbf{1}_{x_0 \le x < x_1} \sin(x - x_0)^2 + \mathbf{1}_{x_1 \le x} - (\sin(x_1 - x_0)^2 - \beta(x - x_1)))$, given parameters $x_1 > x_0 > 1/l_{\max}$ and $\alpha, \beta > 0$. The reference distribution $\pi$ is a log-normal distribution such that $\ln(x) \sim \mathcal{N}(1.5, 1.5)$. Note that although here the rare event is $\{S^* \ge l_{\max}\}$ rather than $\{S^* > l_{\max}\}$, it does not change the behavior of the AMS or ARMS algorithms. The parameter are set to $x_0 = 0.5$,

$x_1 = 5, \alpha = 15, \beta = 0.1$ and $l_{\max} = 90$, corresponding to the rare event probability (given in closed-form) and given in the table above.

Figure 1 shows a typical evolution of the spline-based score approximation as well as the selected snapshots over the iterations of the algorithm.

### 4.1.2 Example #2: reduced basis approximation of a PDE score

In our second numerical experiment, the score function is defined as $S^*(x) = \|\Psi^*(x)\|_p$ with $p = 1$ or $p = \infty$, $\Psi^*(x)$ being the high-fidelity numerical approximation of the solution of a bi-dimensional PDE parametrized by a q-dimensional vector $x = (x_1 \cdots x_q)^\top \in \mathcal{X} = \mathbb{R}_+^q$. The chosen PDE corresponds to the well-studied thermal block problem, which models heat diffusion in an heterogeneous media. The reference distribution for $x$ will be a log-normal distribution with independent components with $\ln(x_i) \sim \mathcal{N}(1.5, 1.5)$, for $i = 1, \dots, q$. The probability of interest will correspond to the occurrence of the temperature exceeding the given thresh-
old $l_{\max}$ either in average (Examples #2a and #2c) or locally (Example #2b) over the spatial domain $\Omega = [0, 1]^2$. PDE solutions will be approximated using the popular reduced basis (RB) method.

We begin by describing the parametric PDE and the high-fidelity numerical approximation of its solution. The temperature function $\Psi$ over $\Omega$ is ruled by the elliptic problem

$$ \begin{cases} -\nabla \cdot (\kappa(\cdot; x)\nabla\Psi) = 1 & \text{in } \Omega \\ \Psi = 0 & \text{on } \partial\Omega \end{cases} $$

where the function $\kappa(\cdot; x) : \Omega \to \mathbb{R}_+$ represents a parametric diffusion coefficient function. This PDE models the stationary temperature field $\Psi$ of an heterogeneous medium which is uniformly heated in the domain and thermostatted at the boundary.

The high-fidelity approximation $\Psi^*(x)$ of the solution of this diffusion problem entails the solution of a variational problem. Precisely, let $\mathcal{V}$ denote the Hilbert space $\mathcal{H}_0^1(\Omega)$ and $\|v\|_{\mathcal{V}} = \|\nabla v\|_{L^2}$ the norm induced by the inner product defined over $\mathcal{V}$. For any $x \in \mathcal{X}$, we define the solution $\Psi^*(x) \in \mathcal{V}$ of the variational problem $a(\Psi(x), v; x) = \int_{\Omega} v d\Omega$, $\forall v \in \mathcal{V}$, where $a : \mathcal{V} \times \mathcal{V} \to \mathbb{R}$ is the bilinear form $a(u, v; x) = \int_{\Omega} \kappa(\cdot; x) \nabla u \cdot \nabla v d\Omega$. We consider the piecewise constant diffusion coefficient $\kappa(\cdot; x) = \sum_{i=1}^q \mathbf{1}_{\Omega_i}(\cdot) x_i$ with the partition $\Omega = \cup_{i=1}^q \Omega_i$ and $\Omega_i$'s being non-overlapping squares. Setting the parameter space $\mathcal{X} = (\mathbb{R}_+^q)^q$ implies that the bilinear form is strongly coercive, and the Lax-Milgram lemma ensures the variational problem has a unique solution.

Consider next a finite $h$-dimensional subspace $\mathcal{V}_h \subset \mathcal{V}$. The high-fidelity approximation $\Psi^*(x) \in \mathcal{V}_h$ of the solution $\Psi^*(x)$ is defined for any $x \in \mathcal{X}$ by

$$ \Psi^*(x) = \arg \min_{v \in \mathcal{V}_h} \|\Psi^*(x) - v\|^2, $$

with the energy norm induced by the bilinear form $\|v\| = \sqrt{a(v, v; x)}$. The high-fidelity approximation $\Psi^*(x)$ is obtained in practice by solving a large system of $h$ linear equations. Benjamin et al. (2018). For the function score defined as the $L^p$-norm, the rare event $\{S^*(X) = \|\Psi^*(X)\|_{L^p} > l_{\max}\}$ corresponds to reaching a critical temperature $l_{\max}$, e.g. a melting temperature, on average in the domain ($p=1$) or in some point of the domain ($p=\infty$).

We are interested in the construction of an inexpensive approximation to the set of high-fidelity solutions

$$ \mathcal{M} = \{\Psi^*(x) \in \mathcal{V}_h : x \in \mathcal{X}\}, $$

by evaluating only some of its elements. By selecting a well-chosen set of $K$ snapshots $\{\Psi^*(x_i)\}_{i=1}^K$, the RB framework approximates $\mathcal{M}$ in a well-chosen linear subspace $\mathcal{V}_K \subset \mathcal{V}_h$ of much lower dimension $K \ll h$. Using this framework yields $i)$ an approximation $\Psi$ of the high-fidelity solution $\Psi^*$ which is used to compute the approximation $S(x) = \|\Psi(x)\|_{L^p}$ of the score $S^*(x)$, ii) a so-called a posteriori error bound on $\|\Psi^*(x) - \Psi(x)\|_{\mathcal{V}_h}$ denoted by $\Delta\Psi(x)$. The high-fidelity solution, the RB approximation and the a posteriori error bounds are computed via the *pymor* Python® toolbox Milk et al. (2016) available at https://github.com/pymor/pymor.

To use the ARMS algorithm, we need to design an error bound $E(x)$ majorizing the absolute error on the score $|S^*(x) - S(x)|$, relying on the a posteriori error bound on $\|\Psi^*(x) - \Psi(x)\|_{\mathcal{V}_h}$. In the case of $p=1$, the error bound satisfies

$$ \begin{aligned} |S^*(x) - S(x)| &= |\|\Psi^*(x)\|_{L^1} - \|\Psi(x)\|_{L^1}| \\ &\le \|\Psi^*(x) - \Psi(x)\|_{L^1} \le \|\Psi^*(x) - \Psi(x)\|_{L^2} \\ &\le \pi^{-1}\|\Psi^*(x) - \Psi(x)\|_{\mathcal{V}_h} \\ &\le \pi^{-1}\Delta\Psi(x) = E(x), \end{aligned} $$

where we have used the triangle inequality, the Cauchy-Schwarz inequality, the Poincaré inequality with the constant $\pi^{-1}$ on the unit square domain and the a posteriori error bound. As the finite element basis is an interpolating one, we will assume $h$ sufficiently large so that we can compute the $L^1$-norm as the average of the finite element coefficients. In the case of $p = \infty$, except in the one-dimensional case, one can not in general rigorously upper bound $\|\Psi^*(x) - \Psi(x)\|_{L^\infty}$ by $\|\Psi^*(x) - \Psi(x)\|_{\mathcal{V}_h}$. We will nevertheless use for this experiment the rough approximation $|S^*(x) - S(x)| \approx \|\Psi^*(x) - \Psi(x)\|_{\mathcal{V}_h} \le \Delta\Psi(x) = E(x)$, possibly violating assumption (1.2).

Let us comment on two important features of the RB approximation in regards to our rare event simulation framework. First, we emphasize that our algorithm enriches the RB approximation space using a generalization of the popular weak greedy algorithm. Indeed, the latter selects at each iteration the snapshot in the solution set $\mathcal{M}$ which possess the worst a posteriori error bound for the current RB approximation. This setup corresponding in our algorithm to the degenerated case where the learning parameter is set to $\tau = +\infty$. At the other end, the setup where $\tau = 0$ will correspond to sampling according to the proposal, independently of the approximation error. In the general setup where $0 < \tau < \infty$, snapshot selection will balance these two cases, as the RB will be enriched with a sample drawn according to a weighted version of the proposal distribution, with weights proportional to the function $e^{\tau E^{(k)}}$.
Second, as in standard RB approaches, we here extensively rely on the affine parametric dependence of the RB approximation by splitting the assembly of the reduced matrices and vectors in an *off-line* and *on-line* phases. The former, which is performed once for all at each RB update, entails the computation of all the *h*-dependent and *x*-independent structures. In the case of the finite element discretization of our elliptic problem, computing scalar products, matrix-vector products or the high-fidelity solution has a complexity of $O(h)$, but with a very large constant (thanks to the generalized minimal residual method Baker et al. (2005) used to solve sparse linear system in pymor). It thus requires off-line $O(hqK^2)$ operations for computing reduced matrices and vectors used to assemble the RB system, and $O(hq^2K^2)$ operations for computing those used to build the error estimator. In the latter, on-line phase, for a given parameter $x \in \mathcal{X}$, we solve the RB system in $O(K^3 + qK^2)$ and provide the error bound in $O(q^2K^2)$, with a cost independent on $h$. For further details, we refer the reader to the introductory book on RB Benjamin et al. (2018).

As there is no analytic value for the rare event probabilities sought, we perform standard AMS calculations with a score based on this high-fidelity approximation $\Psi^*(x)$ and using $N = 1e3$ particles. Averaging 10 runs, we obtain, after intensive computation, an empirical Monte Carlo approximation for $p^* = \pi(\mathbf{1}_{S^*>l_{\max}})$ which is presented in the table above.

Figure 2 and 3 show typical evolutions along the ARMS algorithm iterations for Examples #2a and #2c. We observe that in the first iterations of the algorithm, the RB approximation of the high-fidelity solution of the PDE is too rough to characterize the rare event of interest. While the number of iterations increases, RB gains in accuracy, enabling importance sampling of the rare event.

In the context of Example #2a, a typical rare event sampled by ARMS corresponds to a diffusion coefficient one order of magnitude lower than other diffusion coefficients. To provide a physical interpretation of this rare event, note that each of the 4 subdomains touches the zero boundary condition, so that no coefficient configuration can isolate a subdomain. Consequently, increasing the average temperature can only be achieved by slowing down the heat (contributed by the homogeneous source) lost by the system, which corresponds to configurations with low diffusion coefficients. As these configurations lie in the tail of the reference distribution, such events are rare, and of these events, the least rare corresponds to a decrease in one of the 4 parameters.

In the context of Example #2c, a typical rare event sampled by ARMS corresponds to setting the value of most diffusion coefficients touching the zero boundary condition, one order of magnitude lower than the value of other diffusion coefficients. The physical interpretation of this rare event is clear: isolating the domain boundary slows the loss of heat con-

tributed by the source, and therefore induces an increase in the average temperature within the domain. Note that such a configuration is extremely rare (on the order of $1e-14$) according to the reference distribution, as 11 coefficients out of 25 must be sampled in the tail of the reference distribution.

## 4.2 ARMS setup & evaluation criteria

We now detail the configuration of the ARMS algorithm. The mutation Markov Chain Monte Carlo transition is a Metropolis accept/reject applied on an Ornstein-Uhlenbeck proposal, with an adaptive variance keeping the average acceptance rate over particles in a reasonable range (see e.g. Garthwaite et al. (2016)).

In the early iterations of our algorithm ($k \le K_0$), snapshots are sampled according to the penalized importance distribution where the learning parameter is set to $\tau = +\infty$ (which corresponds to the weak greedy algorithm setting in the RB framework); then for $k > K_0$ (once the sampled snapshots have hit $j_0$ times the rare event set), the snapshots are sampled according to the current proposals (which is equivalent to sample the penalized importance distributions with $\tau = 0$).

For each experiment, we give in the table below the ARMS parameters: the number of particles $N$, the number of particles killed at each selection step $M$, the snapshot budget $K$, the number of Markov transitions at each mutation step $T$, the number of independent ARMS runs for each problem $R$, the number of hits before the estimator is triggered $j_0$ and the entropic threshold $c$. The stopping parameter is set either close to machine precision ($\epsilon > 0$) or such that the reduced score updates never stop ($\epsilon = 0$).

<table><thead><tr><th></th><th>N</th><th>θ = M/N</th><th>K</th><th>T</th><th>R</th><th>j<sub>0</sub></th><th>c</th></tr></thead><tbody><tr><td>Example #1</td><td>[2.0e2, 8.0e4]</td><td>0.3</td><td>300</td><td>1e2</td><td>100</td><td>5</td><td>5.0e-2</td></tr><tr><td>Example #2a</td><td>[2.5e2, 4.0e2]</td><td>0.1</td><td>300</td><td>3e1</td><td>50</td><td>2</td><td>3.0</td></tr><tr><td>Example #2b</td><td>[1.0e2, 4.0e3]</td><td>0.3</td><td>200</td><td>3e1</td><td>100</td><td>5</td><td>5.0e-2</td></tr><tr><td>Example #2c</td><td>[1.0e2, 2.0e2]</td><td>0.1</td><td>300</td><td>3e1</td><td>50</td><td>2</td><td>3.0</td></tr></tbody></table>

**Remark 4.1** (Range for the entropic threshold $c$) $c$ is taken in the range [0.05, 3] which corresponds to the probabilities of the current surrogate rare event under the proposal given by [0.05, 0.95], arguably a reasonable range of tuning.

The initial surrogate is built using the 10 snapshots drawn according to the reference density for Example #1, while we build an initial 5-dimensional RB approximation for Examples #2a, #2b and #2c.

We plot the median and quantiles related to the IS estimate $\hat{p}_{IS}$ and to the AMS estimate $\hat{p}_{AMS}$ as a function of the number of snapshots $k$. ARMS performance is quantified by
<figure><img src="image_5.png" /><figcaption>Fig. 2 Illustration of the algorithm adaptive behavior (Example #2a). The figure displays the high-fidelity solution of the PDE (left), its reduced basis approximation (middle left) and the related absolute error (middle right), evaluated at the snapshot $X^{k+1}$ (right) sampled at iteration $k$ by our ARMS algorithm.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 3 Illustration of the algorithm adaptive behavior (Example #2c). The figure displays the high-fidelity solution of the PDE (left), its reduced basis approximation (middle left) and the related absolute error (middle right), evaluated at the snapshot $X^{k+1}$ (right) sampled at iteration $k$ by our ARMS algorithm.</figcaption></figure>

evaluating the expected cost as a function of the relative root expected square error. More specifically, after $k$ snapshots, the expected relative square error is approximated over runs

by the empirical average

$$ \frac{\mathbb{E}(\hat{p} - p^*)^2}{p^{*2}} \approx \frac{1}{R} \sum_{r=1}^{R} \frac{(\hat{p}_r - p^*)^2}{p^{*2}}, \quad (4.1) $$
<figure><img src="image_7.png" /><figcaption>Fig. 4 Influence of the entropic criterion threshold c (Example #1).</figcaption></figure>

Upper and lower left plots: median (in white) and quantiles (in shades of blue-grey) related to the IS estimate (3.11) as a function of the number of snapshots (equal to the number of iteration k of the algorithm) for N = 1e3, a reduced model gain of ×25, and 3 different values of

the worst log cost threshold c. The red line denotes the true rare event probability. Lower right plot: expected relative square error (4.1) as a function of the expected cost (4.2) (grey curves) for the IS estimate (3.11). Comparison to the performance of standard AMS with S* and N ∈ [1e2, 1e4] (blue curve).

where $\hat{p}_r$ denotes, at the r-th run and after k snapshots, either AMS estimate $\hat{p}_{AMS}$ given in (3.12), or the importance sampling estimate $\hat{p}_{IS}$ given in (3.11) (with the convention that $\hat{p}_{IS} = 0$ if $H = 0$). After k snapshots, the expected cost (in units of the cost of evaluating S*) is approximated over the runs by the empirical average

$$ \text{cost}(k) \approx k + \frac{1}{R} \sum_{r=1, k'=1}^{R, k} N_r^{(k')} \text{gain}^{(k')} \quad (4.2) $$

with $\text{gain}^{(k')}$ representing the ratio of the cost of evaluating the reduced score $S^{(k')}$ to the cost of the true score $S^*$, and where $N_r^{(k')}$ represents the number of evaluations of the reduced score at the k'-th iteration and r-th execution of the ARMS algorithm.

Please note that for the first example, $\text{gain}^{(k')}$ will be set artificially to emulate realistic cases in which the cost of

evaluating the surrogate is indeed cheaper than that of the true model (see Figure 5).

We will compare the performance of ARMS with the empirical average of the estimates obtained by performing R standard AMS runs with the true score $S^*$.

## 4.3 Results

We begin by analyzing the results obtained with the toy model of Example #1. Figure 4 evaluates the influence of the threshold c on the entropic criterion of Section 2.2. The three quantiles show the distribution of the IS estimate as a function of the number of snapshots for different values of the entropy criterion threshold c. As expected, we note (upper left) that a too large threshold produces over 40% of erroneous IS estimates, associated with zero probability regardless of the number of reduced score updates. Although the remaining 60% of estimates nevertheless become accurate after approx-
imately 60 updates of the reduced score, this demonstrates that increasing the critical level too quickly tends to concentrate particles around unimportant local minima. This results in sampling snapshots and refining the reduced score exclusively in regions of state space unrelated to the rare event of interest. To avoid this undesirable scenario, it is possible to set the threshold to its minimum value $c_{min}$ (upper right), which has the disadvantage of delaying the convergence of the IS estimate. Indeed, the plot shows convergence of all estimates towards a non-zero value close to $p^*$, but with a considerable delay, since in this case about a hundred more updates of the reduced score are required. A well-tuned value between these extrema (lower left) achieves convergence after about 70 snapshots, with a characterization of the distribution of IS estimates comparable to the safe distribution obtained with the minimum value.

This trend is confirmed by plotting (lower right) the relative root expected square error (4.1) as a function of the expected cost (4.2). In the case of a not-too-high threshold $c_{min}$ or $28c_{min}$, given a certain precision (i.e., after sampling enough snapshots), we observe a gain in expected cost around one order of magnitude. Clearly, a value of $42c_{min}$ is too high and do not achieve comparable accuracy, regardless of the number of snapshots. As we will analyze later, we point out that the examples with PDE-based scores and the $L^1$-norm (see the results for Examples #2a and #2b) were more robust, and we were able to use $c = 3$, a reasonable maximum value.

Figure 5 illustrates how the gain in terms of expected cost over accuracy can be dramatic, as the cost of evaluating $S^{(k)}$ becomes negligible compared to $S^*$.

Figure 6 shows that the AMS estimate (3.12) significantly underestimates the true probability, demonstrating that relying solely on reduced scores $S^{(k)}$ (here spline approximations) of the true score $S^*$ induces a significant bias. Indeed, in our specific example, the non-differentiability of $S^*$ at the rare event boundary leads to high-frequency oscillations (which tend to decrease with the number of interpolation points, but never disappear completely) in the spline model. The consequence is that a significant proportion of particles in the rare event set $\{x \in \mathcal{X} : S^*(x) \ge l_{max}\}$ are to be excluded from the set $\{x \in \mathcal{X} : S^{(k)}(x) \ge l_{max}\}$, leading to a negative bias. On the contrary, we observe that the IS estimate (3.11) is not biased, independently of the problem of convergence of $S^{(k)}$ to $S^*$.

The evaluation above focused on the ARMS algorithm in which the simulation is restarted after each reduced score update and in which the stop criterion introduced in Section 2.5 is never triggered ($\epsilon = 0$). Figure 7 extends the analysis by evaluating the influence of bridging the simulations between two iterations, as presented in Section 3.3. The entropic criterion is set to $c = 28 \times c_{min}$. As shown in

<figure><img src="image_8.png" /><figcaption>Fig. 5 Influence of the reduced score gain (Example #1). Expected cost (4.2) as a function of the relative expected squared error (4.1) for the IS estimate (3.11), for different values of the ratio between the cost of evaluating $S^{(k)}$ and $S^*$ (grey curves) and for $N = 1e3$ and the threshold $c = 28 \times c_{min}$. Comparison to the performance of standard AMS with $S^*$ and $N \in [1e2, 1e4]$ (blue curve).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 6 Bias of the reduced score based estimate (Example #1). Median (in white) and quantiles (in shades of blue-grey) related to the AMS estimate (3.12) as a function of the number of snapshots; the red line denotes the true rare event probability. Number of samples $N = 1e3$, entropic criterion set to $c = 28 \times e_{min}$, no bridging, never stopping the reduced score update ($\epsilon = 0$) and gain between surrogate and exact is assumed to be 25.</figcaption></figure>

Figure 7, the gain in computational cost brought about by bridging is offset by an increase in variance. One can remark that, when using the bridging procedure, and for a given sample size, we cannot expect to increase the accuracy for the IS estimate (3.11) indefinitely, contrary to when the simulation is restarted after each reduced score update. Indeed, the estimate depends crucially on the variance of the normalization constant $p_{l^{(k)}}^{N,(k)}$ which generates a residual error (upper right box of Figure 7). In the case of bridging, this residual error will indeed propagate across the iterations according to (3.9) instead of being re-estimated by restarting a fresh simulation at each reduced score update. Nevertheless, the plot of expected squared error versus predicted cost in the
<figure><img src="image_10.png" /><figcaption>**Fig. 7** Influence of bridging (Example #1). Median (in white) and quantiles (in shades of blue-grey) related to the IS estimate (3.11) (above) and of the cost (4.2) (above), as a function of the number of snapshots when bridging simulation or not; the red line denotes the true rare event probability. Number of samples *N* = 1e3, entropic criterion set to *c* = 28 × *e*<sub>min</sub>,</figcaption></figure>

stopping updates of the reduced score after convergence only in the case of bridging and gain between surrogate and exact is assumed to be 25. Note that, although imperceptible in the log-log plot at bottom right, the cost increases strictly with the number of snapshots.

upper left-hand panel of Figure 8 shows that, in the case of small sample sizes yielding a relative square error above 20%, some saving of order of 2 is achieved by a bridging procedure compared to restarting the simulation, while both procedures are equivalent for larger sample sizes.

On the other hand, we notice that enabling the stopping criterion to stop the main ARMS loop without bridging is a bad idea as it induces a spread around the mean of the distribution of the IS estimate as illustrated in Figure 9. Indeed, similarly to the above explanation, the IS estimate will be affected by the apparent bias of the estimates $p_{l^{(k)}}^{N,(k)}$ in (3.11), all identical for indices $k$ above some random time (determined by the stopping criterion). The phenomenon does not appear when the stopping criterion is off (Figure 7), since we restart the simulation and calculate a new normalization constant at each iteration. The stopping criterion is, in fact, simply used to marginally reduce the cost in the bridging procedure, without changing the outcome (not shown in plot).

We then move on to Example #2, the PDE-based scores and their RB approximations. Overall, the trends of the algorithm for the RB approximation of a PDE score are very similar to those of the previous example. We thus omit the precise evaluation made for the previous example, and jump directly to the analysis of performances in Figure 8. In the case of Examples #2a and #2b, we observe that about 10 to 20 more snapshots (i.e., ARMS iterations) are needed to reach convergence, with a median estimate around $p^*$. In the case of Example #2c, about 170 to 200 additional snapshots are needed to reach a good accuracy of the reduced model and converge towards $p^*$. Indeed, in the latter situation, the parameter space has a much higher dimensionality (25 dimensions versus 4), and learning the reduced model for parameters in the (initially unknown) region where rare events occur is a difficult task performed progressively over the course of ARMS iterations.
<figure><img src="image_11.png" /><figcaption>**Fig. 8 Performance of ARMS (Examples #1 & 2).** Relative expected square error (4.1) as a function of the expected cost (4.2) of the IS estimate (3.11) at the final iteration (once the snapshots budget *K* is</figcaption></figure>

exhausted), for a varying sample size *N*. Curves can be compared to standard ARMS estimation using the score $S^*$, for a varying sample size *N*.

In all examples, Figure 8 shows that, in terms of the “cost × square error” criterion, the proposed ARMS method enables important savings as compared to the use of the true model alone. Savings are consistently around a factor 10 in all cases. As mentioned in Section 4.1.2, in the case of Example #2b, ARMS uses error estimates that are not hard bounds for the $L^\infty$ score. Nevertheless, it is surprising that ARMS manages to accurately estimate the desired probability and achieve performances comparable to those obtained with the $L^1$ score.

In terms of a 'CPU time × square error' criterion, the savings picture is only similar for low-dimensional parameter spaces of RB models. As shown in the table below (for $K = 200$ and $h = 5101$), for Example #2a (where $q = 4$), the gain in cost for reaching a square error of $1e-1$ is similar to the gain in CPU time. However, for Example #2c (where $q = 25$), although there is a cost saving, there is no CPU time saving.

<table><thead><tr><th></th><th>CPU</th><th></th><th>CPU×(4.1)</th><th>(4.2))×(4.1)</th></tr></thead><tbody><tr><td></td><td>S(x)</td><td>S(x) &amp; E(x) v.s.</td><td>ARMS</td><td>ARMS</td></tr><tr><td></td><td>v.s.</td><td>S*(x)</td><td>v.s.</td><td>v.s.</td></tr><tr><td></td><td>S*(x)</td><td></td><td>AMS</td><td>AMS</td></tr><tr><td>#2a (q=4)</td><td>×63</td><td>×29</td><td>×17</td><td>×18</td></tr><tr><td>#2c (q=25)</td><td>×12</td><td>×3</td><td>×1</td><td>×3</td></tr></tbody></table>

In particular, the table shows that varying the parameter space dimension *q* of the RB model from 4 to 25 slows the speed of solving an RB system (scaling in $O(K^3+qK^2)$) and of estimating the error (scaling in $O(q^2K^2)$). This explains in part the ARMS performance drop. The latter is however mainly due to the dramatic change in the balance of CPU time usage by its routines: we observed in our experiments that the cost of an RB update (scaling in $O(hq^2K^2)$) is almost negligible for $q = 4$ compared to the cost of the simulations
<figure><img src="image_12.png" /><figcaption>**Fig. 9 An ineffective approach: using the stopping criterion without bridging (Example #1).** Median (in white) and quantiles (in shades of blue-grey) related to the IS estimate (3.11) as a function of the number of snapshots when stopping the reduced score updates ($\epsilon > 0$) without bridging; the red line denotes the true rare event probability. Number of samples $N = 1e3$, entropic criterion set to $c = 28 \times c_{min}$ and gain between surrogate and exact is assumed to be 25.</figcaption></figure>

to the critical levels and the bridging routines (both scaling in $O(M(q^2K^2 + K^3))$), but requires more than three quarter of the computational burden for $q = 25$. Let us remark that the core of the ARMS algorithm relies on an AMS methodology for the tempering and bridging routines, which is known to be stable in high dimensions Beskos et al. (2014). The drop in computation saving in high dimensions is here inherent to the RB formalism and not the ART algorithm itself. More precisely, the reduced basis assembly which is intended to be performed once for all off-line is done in the ARMS algorithm on-line at each of its iterations. The quadratic computational cost of this procedure prohibits a high dimension $q$. We should nevertheless mention that there are possible remedies to avoid this slowdown. Indeed, to avoid excessive cost of RB updates not affordable on-line, recent works challenge the RB off-line/on-line paradigm. In particular, the adaptive reduced modeling scheme of Peherstorfer and Willcox (2015) proposes inexpensive on-line updates to the RB basis while maintaining its dimension. Other popular works on active subspace offer ways to circumvent the problem by learning a low-dimensional representation of the parameter space Constantine et al. (2014).

# 5 Conclusions and Perspectives

In addition to demonstrating the unbiasedness of our estimator in an idealized setting, we evaluate in our numerical experiments the empirical convergence of the ARMS algorithm on a toy model and a rare-event problem based on a more realistic PDE. Our study confirms that entropy criteria play a crucial role in ARMS. In particular, as expected,

too high a value of the logarithm of the cost in the entropy criterion $i$) can take ARMS into spurious regions of the state-space remote from the region of interest, resulting in a large variance of the IS estimator. Too small a value gives a more desirable variance but a slower rate of convergence. We show empirically that the approximation procedure using the entropy criterion $iii)$ significantly reduces simulation cost (around a factor 10) while achieving a comparable squared error. Similarly, in both experiments, we observe a squared error gain of the order of one order of magnitude for a given simulation cost. Our numerical experiments show that bridging is a good option for small computational budgets, with a saving of around 2 compared with restarting the simulation at each iteration. However, for larger computational budgets, both procedures yield comparable cost-error trade-offs.

The prospects arising from this work are numerous. These include answering open theoretical questions, including the consistency of ARMS, which is proven in the present paper only in a formal sense. No straightforward answers seem available, as we are confronted with a non-standard framework in which interacting particles in AMS simulations have distributions with varying supports which depend on random approximations of the reduced score. Future research directions also include extending ARMS to smooth proposal distributions, with the advantage of removing nesting and domination constraints and the possibility of being more flexible in the choice of error estimates that need not have strict error bounds.

# Appendix A. Sample Size for Importance Sampling

## A.1 Generic quantification rules

Importance sampling of some generic distribution $\eta \propto \gamma$ known up to a normalizing constant $p = \gamma(1)$ using a proposal distribution $\mu$ is based on the following identity

$$ \gamma(\varphi) = p \eta(\varphi) = \mathbb{E}_{\mu}[f(X)\varphi(X)] \quad \forall \varphi: \mathcal{X} \to \mathbb{R} $$

where $f(x) = \frac{d\gamma}{d\mu}(x)$ is computable in closed-form.

Unfortunately, importance sampling becomes quickly useless if the target distribution diverges too widely from the proposal, a typical phenomenon in high dimension. This can be readily seen by computing the relative variance of the importance sampling estimator for $\varphi = 1$:

$$ \mathbb{V}\mathrm{ar}_{\mu}(f(X))/p^2 = \eta\left(\frac{d\eta}{d\mu}\right) - 1 = e^{\mathrm{Ent}_2(\eta|\mu)} - 1, \quad (A.1) $$
where we have denoted $Ent_2(\eta | \mu) = \ln \eta(\frac{d\eta}{d\mu})$ the order 2 Rényi entropy (also known as $\chi^2$-divergence). The latter scales linearly in $d$ if $\mu$ and $\eta$ are a $d$-fold product measure, generating a variance exponential with $d$.

As is thoroughly argued in Chatterjee and Diaconis (2018) (see also Cérou et al. (2022)), the variance $(e^{Ent_2(\eta|\mu)} - 1)$ is in general a pessimistic quantification of importance sampling, whose cost is better estimated using a $\mathbb{L}^1$-error, or a success probability within a tolerance precision. For a target $\eta$ with proposal $\mu$, it has been proven in Chatterjee and Diaconis (2018) that the cost of importance sampling – in terms of sample size $K$– is roughly given (up to a reasonable tail condition on the log-density $\ln \frac{d\eta}{d\mu}$), by the exponential of the relative entropy between $\eta$ and $\mu$: $K \sim e^{Ent(\eta|\mu)}$.

**Remark A.1** Jensen inequality

$$Ent_2(\eta | \mu) = \ln \eta(\frac{\eta}{\mu}) \geq \eta(\ln \frac{\eta}{\mu}) = \mathrm{Ent}(\eta | \mu),$$

shows that indeed, variance is a more pessimistic estimator of the required sample size for importance sampling.

## A.2 Particular case of interest

We will particularize the above discussion to the family of distributions studied in our paper: distributions of propositions of the form $\mu \propto \mathbf{I}_B d\pi$ and a target $\eta \propto \mathbf{I}_A d\pi$ such that $A \subseteq B$.

In accordance with the importance sampling identity, the normalization $p$ is estimated without bias using the random variable $f(X) = \pi(B) \mathbf{I}_A(X)$ where $X \sim \mu$.

In this remarkable case, it turns out that the relative variance and the exponential of the relative entropy reduces to the same quantity. In other words, we reach the limit of Jensen's inequality in the Remark A.1

$$Ent_2(\eta | \mu) = \mathrm{Ent}(\eta | \mu) = \ln \frac{\pi(B)}{\pi(A)}.$$

Therefore, in our study, we will refer indifferently to the exponential of the relative entropy $e^{Ent(\eta|\mu)}$ or to the relative variance $e^{Ent_2(\eta|\mu)}$ to designate the importance sampling cost of the target $\eta$ by the proposal $\mu$.

# Appendix B. Consistency of the Estimator

This section analyses the idealized algorithm obtained with $N = +\infty$ many particles (clones) in ARMS;

For simplicity (but the general case is similar), we study the estimator (2.9) in the case of a test function $\varphi = \mathbf{I}$.

We thus consider the rare event probability estimator: $\hat{p}_l^H = \hat{\gamma}_l^H(\mathbf{I})$.

**Proposition B.1** Assume that for each $k$ the support of $\eta_{l(k)}^{(k)}$ contains the target rare event $\{S^* > l_{\max}\}$. The estimator $\hat{p}_l^H$ is unbiased, i.e., $\mathbb{E} \hat{p}_l^H = p_l^*$, and its relative variance is given by

$$\mathbb{V}_{\mathrm{ar}}[\hat{p}_l^H]/p_l^{*2} = \frac{1}{H^2} \sum_{k=K_0+1}^{K_0+H} \left( \frac{\mathbb{E}(p_{l(k)}^{(k)})}{p_l^*} - 1 \right), \quad (A.2)$$

where the random variable $p_{l(k)}^{(k)}$ is a function of the past snapshots.

**Proof** We first show that for any test function $\varphi$ the random sequence $H \mapsto M_H^\varphi \stackrel{\mathrm{def}}{=} H(\hat{\gamma}_l^H(\varphi) - \gamma_l^*(\varphi))$, is a martingale with respect to the filtration generated by the snapshots. Indeed, defining this filtration by $\mathcal{F}_H = \sigma(X_1, \dots, X_{K_0+1+H})$, it follows from the definition (2.9) that

$$\begin{aligned} \mathbb{E}[M_H^\varphi | \mathcal{F}_{H-1}] &= \mathbb{E}[\sum_{k=K_0+1}^{K_0+H} f_l^{(k)}(X_{k+1})\varphi(X_{k+1}) \\ &\quad - H\gamma_l^*(\varphi) | \mathcal{F}_{H-1}], \\ &= M_{H-1}^\varphi + \mathbb{E}[f_l^{(K_0+H)}(X_{K_0+H+1}) \\ &\quad \varphi(X_{K_0+H+1}) - \gamma_l^*(\varphi) | \mathcal{F}_{H-1}], \\ &= M_{H-1}^\varphi + \pi(\varphi | \mathbf{I}_{S^{*>l}}) - \gamma_l^*(\varphi) = M_{H-1}^\varphi. \end{aligned}$$

Now, since $M_H^\varphi$ is a martingale, the estimator is unbiased $\mathbb{E} \hat{\gamma}_l^H(\varphi) = \gamma_l^*(\varphi)$, and its variance is given by its quadratic variations given by the sum over $k$ of the local variances. In particular, using (A.1), the variance of the rare event probability estimator is then

$$\mathbb{V}_{\mathrm{ar}}[\hat{p}_l^H]/p_l^{*2} = \frac{1}{H^2} \sum_{k=K_0+1}^{K_0+H} \left( \mathbb{E}(e^{Ent_2(\eta_l^*|\eta_{l(k)}^{(k)})}) - 1 \right),$$

which for $\eta_l^* \propto \mathbf{I}_{S^{*>l}} d\pi$ and $\eta_{l(k)}^{(k)} \propto \mathbf{I}_{S^{(k)>l(k)}} d\pi$ yields (A.2). $\square$

**Author Contributions** F.C., P.H and M.R wrote the main text. P.H. performed the numerical simulations

**Data Availability** No datasets were generated or analysed during the current study.

## Declarations

**Competing interests** The authors declare no competing interests.
## References

Baker, A.H., Jessup, E.R., Manteuffel, T.: A technique for accelerating the convergence of restarted gmres. SIAM J. Matrix Anal. Appl. 26(4), 962–984 (2005)

Balesdent, M., Morio, J., Marzat, J.: Kriging-based adaptive importance sampling algorithms for rare event estimation. Struct. Saf. 44, 1–10 (2013)

Bect, J., Li, L., Vazquez, E.: Bayesian subset simulation. SIAM/ASA Journal on Uncertainty Quantification 5(1), 762–786 (2017)

Beskos, A., Crisan, D., Jasra, A.: On the stability of sequential Monte Carlo methods in high dimensions. Ann. Appl. Probab. 24(4), 1396–1445 (2014)

Boer, P.-T.D., Kroese, D.P., Mannor, S., Rubinstein, R.Y.: A tutorial on the cross-entropy method. Ann. Oper. Res. 134(1), 19–67 (2005)

Cannamela, C., Garnier, J., Iooss, B.: Controlled stratification for quantile estimation. The Annals of Applied Statistics 2(4), 1554–1580 (2008)

F. Cérou, P. Del Moral, T. Furon, and A. Guyader. Sequential Monte Carlo for rare event estimation. *Stat. Comput.*, 22(3):795–808, 2012

Cérou, F., Guyader, A.: Fluctuation analysis of adaptive multilevel splitting. Annals of Applied Probability 26(6), 3319–3380 (2016)

Cérou, F., Del Moral, P., Furon, T., Guyader, A.: Sequential Monte Carlo for rare event estimation. Stat. Comput. 22(3), 795–808 (2012)

Chatterjee, S., Diaconis, P.: The sample size required in importance sampling. Ann. Appl. Probab. 28(2), 1099–1135 (2018)

Chen, P., Quarteroni, A.: Accurate and efficient evaluation of failure probability for partial different equations with random input data. Comput. Methods Appl. Mech. Eng. 267, 233–260 (2013)

Constantine, P.G., Dow, E., Wang, Q.: Active subspace methods in theory and practice: applications to kriging surfaces. SIAM J. Sci. Comput. 36(4), A1500–A1524 (2014)

Del Moral, P., Doucet, A., Jasra, A.: Sequential monte carlo samplers. J. R. Stat. Soc. Ser. B Stat Methodol. 68(3), 411–436 (2006)

Echard, B., Gayton, N., Lemaire, M.: AK-MCS: an active learning reliability method combining kriging and monte carlo simulation. Struct. Saf. 33(2), 145–154 (2011)

Benjamin Echard, Nicolas Gayton, and Maurice Lemaire. AK-MCS: an active learning reliability method combining kriging and monte carlo simulation. *Structural Safety*, 33(2):145–154, 2011

Gallimard, L.: Adaptive reduced basis strategy for rare-event simulations. Int. J. Numer. Meth. Eng. 120(3), 283–302 (2019)

Garthwaite, P.H., Fan, Y., Sisson, S.A.: Adaptive optimal scaling of metropolis-hastings algorithms using the robbins-monro process. Communications in Statistics-Theory and Methods 45(17), 5098–5111 (2016)

Héas, P.: Selecting reduced models in the cross-entropy method. SIAM/ASA Journal on Uncertainty Quantification 8(2), 511–538 (2020)

Heinkenschloss, M., Kramer, B., Takhtaganov, T.: Adaptive reduced-order model construction for conditional value-at-risk estimation. SIAM/ASA Journal on Uncertainty Quantification 8(2), 668–692 (2020)

Ling, C., Zhenzhou, L., Feng, K., Zhang, X.: A coupled subset simulation and active learning kriging reliability analysis method for rare failure events. Struct. Multidiscip. Optim. 60(6), 2325–2341 (2019)

Milk, R., Rave, S., Schindler, F.: pymor-generic algorithms and interfaces for model order reduction. SIAM J. Sci. Comput. 38(5), S194–S216 (2016)

Moustapha, M., Marelli, S., Sudret, B.: Active learning for structural reliability: Survey, general framework and benchmark. Struct. Saf. 96, 102174 (2022)

Peherstorfer, B., Willcox, K.: Online adaptive model reduction for nonlinear systems via low-rank updates. SIAM J. Sci. Comput. 37(4), A2123–A2150 (2015)

Peherstorfer, B., Willcox, K., Gunzburger, M.: Survey of multifidelity methods in uncertainty propagation, inference, and optimization. SIAM Rev. 60(3), 550–591 (2018)

Peherstorfer, B., Willcox, K., Gunzburger, M.: Survey of multifidelity methods in uncertainty propagation, inference, and optimization. SIAM Rev. 60(3), 550–591 (2018)

Razaaly, N., Crommelin, D., Congedo, P.M.: Efficient estimation of extreme quantiles using adaptive kriging and importance sampling. Int. J. Numer. Meth. Eng. 121(9), 2086–2105 (2020)

Siu-Kui, A., Beck, J.L.: Estimation of small failure probabilities in high dimensions by subset simulation. Probab. Eng. Mech. 16(4), 263–277 (2001)

Ullmann, E., Papaioannou, I.: Multilevel estimation of rare events. SIAM/ASA Journal on Uncertainty Quantification 3(1), 922–953 (2015)

Uribe, F., Papaioannou, I., Marzouk, Y.M., Straub, D.: Cross-entropy-based importance sampling with failure-informed dimension reduction for rare event simulation. SIAM/ASA Journal on Uncertainty Quantification 9(2), 818–847 (2021)

Wagner, F., Latz, J., Papaioannou, I., Ullmann, E.: Multilevel sequential importance sampling for rare event estimation. SIAM J. Sci. Comput. 42(4), A2062–A2087 (2020)

Zhang, J., Taflanidis, A.A.: Adaptive kriging stochastic sampling and density approximation and its application to rare-event estimation. ASCE-ASME Journal of Risk and Uncertainty in Engineering Systems, Part A: Civil Engineering 4(3), 04018021 (2018)

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.