Stratified Splitting for Efficient Monte Carlo Integration

TL;DR

This paper introduces Stratified Splitting, a new Sequential Monte Carlo method that efficiently estimates high-dimensional integrals, including rare-event probabilities, with unbiasedness and polynomial complexity.

Contribution

The paper presents a novel Stratified Splitting technique that improves Monte Carlo integration efficiency and handles indicator functions, with proven polynomial complexity.

Findings

Accurate estimates for high-dimensional integrals.

Handles indicator functions in rare-event problems.

Achieves polynomial complexity in certain variants.

Abstract

The efficient evaluation of high-dimensional integrals is of importance in both theoretical and practical fields of science, such as data science, statistical physics, and machine learning. However, exact computation methods suffer from the curse of dimensionality. However, due to the curse of dimensionality, deterministic numerical methods are inefficient in high-dimensional settings. Consequentially, for many practical problems, one must resort to Monte Carlo estimation. In this paper, we introduce a novel Sequential Monte Carlo technique called Stratified Splitting. The method provides unbiased estimates and can handle various integrand types including indicator functions, which are used in rare-event probability estimation problems. Moreover, we demonstrate that a variant of the algorithm can achieve polynomial complexity. The results of our numerical experiments suggest that the Stratified Splitting method is capable of delivering accurate results for a variety of integration problems.

Tables

Table 1. Table 1 : The SSA results for the Normal copula credit risk model with 21 21 21 factors and 1 , 000 1 000 1,000 obligors.

$v$	$\hat{c}$	RE	$v$	$\hat{c}$	RE
10000	$1.68\text{×}{10}^{04}$	0.67 %	34000	$3.80\text{×}{10}^{04}$	0.05 %
14000	$2.09\text{×}{10}^{04}$	0.51 %	38000	$4.11\text{×}{10}^{04}$	0.05 %
18000	$2.46\text{×}{10}^{04}$	0.21 %	40000	$4.26\text{×}{10}^{04}$	0.08 %
22000	$2.82\text{×}{10}^{04}$	0.19 %	44000	$4.56\text{×}{10}^{04}$	0.08 %
24000	$2.99\text{×}{10}^{04}$	0.23 %	48000	$4.86\text{×}{10}^{04}$	0.02 %
28000	$3.32\text{×}{10}^{04}$	0.21 %	50000	$5.01\text{×}{10}^{04}$	0.02 %
30000	$3.48\text{×}{10}^{04}$	0.14 %

Equations

$${\mathbb{E}}_f\left[\phi (\mathbf{X})\right]=\sum _{\mathcal{X}}\phi (\mathbf{x})f(\mathbf{x})\text{or}{\mathbb{E}}_f\left[\phi (\mathbf{X})\right]={\int }_{\mathcal{X}}\phi (\mathbf{x})f(\mathbf{x})\text{d}\mathbf{x},$$

$$\mathcal{X}=\bigcup _{t=1}^n{\mathcal{Z}}_t,{\mathcal{Z}}_{t_1}\cap {\mathcal{Z}}_{t_2}=\varnothing \mathrm{for}1\le t_1<t_2\le n.$$

$$f_t(\mathbf{x})$$

$$g_t(\mathbf{x})=f\left(\mathbf{x}\mid \mathbf{x}\in {\mathcal{Z}}_t\right)=\frac{f(\mathbf{x})\text{mathbbm}1\{\mathbf{x}\in {\mathcal{Z}}_t\}}{{\mathbb{P}}_f(\mathbf{X}\in {\mathcal{Z}}_t)}\text{for}t=1,\dots ,n.$$

$$z\stackrel{\mathrm{def}}{=}{\mathbb{E}}_f\left[\phi \left(\mathbf{X}\right)\right]=\sum _{t=1}^n{\mathbb{E}}_f\left[\phi \left(\mathbf{X}\right)\mid \mathbf{X}\in {\mathcal{Z}}_t\right]{\mathbb{P}}_f\left(\mathbf{X}\in {\mathcal{Z}}_t\right)=\sum _{t=1}^n{z}_t.$$

$${\hat{\Phi }}_t=\frac{1}{|{\mathcal{Z}}_t|}\sum _{\mathbf{Z}\in {\mathcal{Z}}_t}\phi (\mathbf{Z})\text{for all}t=1,\dots ,n.$$

$$p_t$$

$$=\prod _{i=0}^{t-1}{r}_i-\prod _{i=0}^t{r}_i=(1-r_t)\prod _{i=0}^{t-1}{r}_i.$$

$$\hat{Z}=R^{-1}\sum _{j=1}^R{\hat{Z}}^{(j)}.$$

$$\mathrm{RE}=\frac{\sqrt{\mathrm{Var}\left(\hat{Z}\right)}}{\mathbb{E}\left[\hat{Z}\right]\sqrt{R}}.$$

$$\mathbb{P}\left(z(1-\epsilon )\le \hat{Z}\le z(1+\epsilon )\right)\ge 1-\delta .$$

$${\mu }_t(A)={\int }_A{f}_t(u)\lambda (\mathrm{d}u),$$

$${\kappa }_t^{\tau }(A\mid \mathbf{x})=\mathbb{P}\left(X_t^{(\tau )}\in A\mid X_t^{(0)}=\mathbf{x}\right),$$

$$\Vert {\kappa }_t^{\tau }(\cdot \mid \mathbf{x})-{\mu }_t{\Vert }_{\mathrm{TV}}=\underset{A}{\sup }\left|{\kappa }_t^{\tau }(A\mid \mathbf{x})-{\mu }_t(A)\right|.$$

$${\tau }_{\mathrm{mix}}(\epsilon ,\mathbf{x})=\min \{\tau :\Vert {\kappa }_t^{\tau }(\cdot \mid \mathbf{x})-{\mu }_t{\Vert }_{\mathrm{TV}}\le \epsilon \}.$$

$$S(\mathbf{w},\mathbf{x})=\sum _{i=1}^k{w}_i{x}_i={\mathbf{w}}^{\top }\mathbf{x}.$$

$$\mathbb{E}\left[S(\mathbf{w},\mathbf{X})\mid S(\mathbf{w},\mathbf{X})\le \gamma \right],$$

$$\mathbb{P}(S(\mathbf{X})\le \gamma )=\mathbb{E}\left[\text{mathbbm}1\{S(\mathbf{X})\le \gamma \}\right].$$

$${\mathcal{X}}_b=\left\{\mathbf{x}\in \{0,1{\}}^k:\sum _{i=1}^k{w}_i{x}_i\le b\right\}\text{for }b\in \mathbb{R},$$

$$\Vert \hat{\pi }-\pi {\Vert }_{\mathrm{TV}}\le \frac{\epsilon a}{4(b-a)},\text{and}m\ge \frac{(b-a^2)\ln (2/\delta )}{2(\epsilon /4{)}^2{a}^2},$$

$$\mathbb{P}\left((1-\epsilon ){\mathbb{E}}_{\pi }[X]\le \frac{1}{m}\sum _{i=1}^m{X}_i\le (1+\epsilon ){\mathbb{E}}_{\pi }[X]\right)\ge 1-\delta .$$

$$m=\frac{(\overline{w}-\underline{w}{)}^2\ln (2/\delta )}{2(\epsilon /4{)}^2{\underline{w}}^2}$$

$${\Vert {\hat{\mu }}_{\gamma }-{\mu }_{\gamma }\Vert }_{\mathrm{TV}}\le \frac{\epsilon \underline{w}}{4(\overline{w}-\underline{w})}.$$

$$m=\frac{(\overline{w}-\underline{w}{)}^2\ln (2/\delta )}{(\epsilon /4{)}^2{\underline{w}}^2}\underset{\overline{w}=\mathcal{O}(\mathcal{P}(k))\underline{w}}{\underbrace{=}}\mathcal{O}(\mathcal{P}(k))\frac{\ln (2/\delta )}{{\epsilon }^2}.$$

$$n=\lfloor \frac{\left({\sum }_{i=1}^k{w}_i\right)-\gamma }{\underline{w}}\rfloor ,$$

$$\mathbb{E}\left[\text{mathbbm}1\{S(\mathbf{X})\le \gamma \}\right]=\mathbb{E}\left[\text{mathbbm}1\{S(\mathbf{X})\le \gamma +n\underline{w}\}\mid S(\mathbf{X})\le \gamma +n\underline{w}\right]\mathbb{P}\left(S(\mathbf{X})\le \gamma +n\underline{w}\right)$$

$$+\mathbb{E}\left[\text{mathbbm}1\{S(\mathbf{X})\le \gamma \}\mid \gamma +n\underline{w}<S(\mathbf{X})\le \sum _{i=1}^k{w}_i\right]\mathbb{P}\left(\gamma +n\underline{w}<S(\mathbf{X})\le \sum _{i=1}^k{w}_i\right)$$

$$=\underset{(\ast )}{\underbrace{\mathbb{E}\left[\text{mathbbm}1\{S(\mathbf{X})\le \gamma \}\mid S(\mathbf{X})\le \gamma +n\underline{w}\right]}}\frac{2^k-1}{2^k}+\left(\sum _{i=1}^k{w}_i\right)\frac{1}{2^k},$$

$$\hat{Z}=\prod _{t=0}^n{\hat{R}}_t.$$

$$\Vert {\hat{\mu }}_t-{\mu }_t{\Vert }_{\mathrm{TV}}\le \epsilon {\underline{r}}_t/32n,\text{and}|{\mathcal{X}}_t|\ge 3072n^2\ln (4n^2/\delta )/{\epsilon }^2{\underline{r}}_t^2,$$

$$\lfloor \frac{\left({\sum }_{i=1}^k{w}_i\right)-\tau }{\underline{w}}\rfloor \le \frac{k\overline{w}}{\underline{w}}\underset{\overline{w}=\mathcal{O}(\mathcal{P}(k))\underline{w}}{\underbrace{=}}\mathcal{O}(\mathcal{P}(k)).$$

$$L=\sum _{i=1}^k{l}_i{X}_i,$$

Full text

1 Introduction

We consider the evaluation of expectations and integrals of the form

[TABLE]

where $\mathbf{X}\sim f$ is a random variable taking values in a set ${\mathscr{X}}\subseteq\mathbb{R}^{d}$, $f$ is a probability density function (pdf) with respect to the Lebesgue or counting measure, and $\varphi:\mathscr{X}\rightarrow\mathbb{R}$ is a real-valued function.

The evaluation of such high-dimensional integrals is of critical importance in many scientific areas, including statistical inference [13], rare-event probability estimation [1], machine learning [41, 25], and cryptography [28]. An important application is the calculation of the normalizing constant of a probability distribution, such as the marginal likelihood (model evidence). However, often obtaining even a reasonably accurate estimate of $\mathbb{E}_{f}\left[\varphi(\mathbf{X})\right]$ can be hard [39]. Deterministic computation methods that use Fubini’s theorem [12] and quadrature rules or extrapolations [11] suffer from the curse of dimensionality, with the number of required function evaluations growing exponentially with the dimension. Many other deterministic and randomized methods have been proposed to estimate high-dimensional integrals. Examples include Bayesian quadrature, sparse grids, and various Monte Carlo, quasi-Monte Carlo, Nested Sampling, and Markov Chain Monte Carlo (MCMC) algorithms [35, 30, 32, 19, 42, 26]. There are also procedures based on the sequential Monte Carlo (SMC) approach [8] that provide consistent and unbiased estimators that have asymptotic normality. However, many popular methods are not unbiased and some are not known to be consistent. For example, Nested Sampling produces biased estimates and even its consistency, when Markov Chain Monte Carlo (MCMC) is used, remains an open problem  [7]. An alternative criterion that one might use to describe the efficiency of algorithms is computational complexity, which often used in a computer science context. Here, an algorithm is considered to be efficient if it scales polynomially (rather than exponentially) in the size of the problem.

In this paper, we propose a novel SMC approach for reliable and fast estimation of high-dimensional integrals. Our method extends the Generalized Splitting (GS) algorithm of [4], to allow the estimation of quite general integrals. In addition, our algorithm is specifically designed to perform efficient sampling in regions of ${\mathscr{X}}$ where $f$ takes small values and $\varphi$ takes large values. In particular, we present a way of implementing stratification for variance reduction in the absence of knowing the strata probabilities. A benefit of the proposed Stratified Splitting algorithm (SSA) is that it provides an unbiased estimator of $\mathbb{E}_{f}\left[\varphi(\mathbf{X})\right]$, and that it can be analyzed in a non-asymptotic setting. In particular, we prove that a simplified version of SSA can provide polynomial complexity, under certain conditions. We give a specific example of a $\#$P complete problem where polynomial efficiency can be achieved, providing a polynomial bound on the number of samples required to achieve a predefined error bound.

The SSA uses a stratified sampling scheme (see, e.g., [39], Chapter 5), defining a partition of the state space into strata, and using the law of total probability to deliver an estimator of the value of the integral. To do so, one needs to obtain a sample population from each strata and know the exact probability of each such strata. Under the classical stratified sampling framework, it is assumed that the former is easy to achieve and the latter is known in advance. However, such favorable scenarios are rarely seen in practice. In particular, obtaining samples from within a stratum and estimating the associated probability that a sample will be within this stratum is hard in general [23]. To resolve this issue, the SSA incorporates a multi-level splitting mechanism [24, 4, 38, 9] and uses an appropriate MCMC method to sample from conditional densities associated with a particular stratum.

The rest of the paper is organized as follows. In Section 2 we introduce the SSA, explain its correspondence to a generic multi-level sampling framework, and prove that the SSA delivers an unbiased estimator of the expectation of interest. In Section 3, we provide a rigorous complexity analysis of the approximation error of the proposed method, under a simplified setting. In Section 4, we introduce a challenging estimation problem called the weighted component model, and demonstrate that SSA can provide an arbitrary level of precision for this problem with polynomial complexity in the size of the problem. In Section 5, we report our numerical findings on various test cases that typify classes of problems for which the SSA is of practical interest. Finally, in Section 6 we summarize the results and discuss possible directions for future research. Detailed proofs are given in the appendix.

2 Stratified splitting algorithm

2.1 Generic multilevel splitting framework

We begin by considering a very generic multilevel splitting framework, similar to [14]. Let $\mathbf{X}\sim f$ be a random variable taking values in a set ${\mathscr{X}}$. Consider a decreasing sequence of sets ${\mathscr{X}}={\mathscr{X}}_{0}\supseteq\cdots\supseteq{\mathscr{X}}_{n}=\emptyset$ and define ${\mathscr{Z}}_{t}={\mathscr{X}}_{t-1}\setminus{\mathscr{X}}_{t}$, for $t=1,\ldots,n$. Note that ${\mathscr{X}}_{t-1}=\bigcup_{i=t}^{n}{{\mathscr{Z}}_{i}}$ and that $\{{\mathscr{Z}}_{t}\}$ yields a partition of ${\mathscr{X}}$; that is,

[TABLE]

We can define a sequence of conditional pdfs

[TABLE]

where $\mathbbm{1}$ denotes the indicator function. Also, define

[TABLE]

Our main objective is to sample from the pdfs $f_{t}$ and $g_{t}$ in (2) and (3), respectively. To do so, we first formulate a generic multilevel splitting framework, given in Algorithm 1.

Note that the samples in $\{{\cal X}_{t}\}_{1\leq t\leq n}$ and $\{{\cal Z}_{t}\}_{1\leq t\leq n}$ are distributed according to $f_{t}$ and ${g}_{t}$, respectively, and these samples can be used to handle several tasks. In particular, the $\{{\cal X}_{t}\}_{1\leq t\leq n}$ sets allow one to handle the general non-linear Bayesian filtering problem [17]. Moreover, by tracking the cardinalities of the sets $\{{\cal X}_{t}\}_{1\leq t\leq n}$ and $\{{\cal Z}_{t}\}_{1\leq t\leq n}$, one is able to tackle hard rare-event probability estimation problems, such as delivering estimates of $\mathbb{P}_{f}\left(\mathbf{X}\in\mathscr{Z}_{n}\right)$ [4, 26, 40]. Finally, it was recently shown by [43] that Algorithm 1 can be used as a powerful variance minimization technique for any general SMC procedure. In light of the above, we propose taking further advantage of the sets $\{{\cal X}_{t}\}_{1\leq t\leq n}$ and $\{{\cal Z}_{t}\}_{1\leq t\leq n}$, to obtain an estimation method suitable for general integration problems.

2.2 The SSA set-up

Following the above multilevel splitting framework, it is convenient to construct the sequence of sets $\{{\mathscr{X}}_{t}\}_{0\leq t\leq n}$ by using a performance function $S:{\mathscr{X}}\rightarrow\mathbb{R}$, in such a way that $\{{\mathscr{X}}_{t}\}_{0\leq t\leq n}$ can be written as super level-sets of $S$ for chosen levels $\gamma_{0},\ldots,\gamma_{n}$, where $\gamma_{0}$ and $\gamma_{n}$ are equal to $\inf_{\mathbf{x}\in{\mathscr{X}}}S(\mathbf{x})$ and $\sup_{\mathbf{x}\in{\mathscr{X}}}S(\mathbf{x})$, respectively. In particular, ${{\mathscr{X}}_{t}=\left\{\mathbf{x}\in{\mathscr{X}}:S(\mathbf{x})\geq\gamma_{t}\right\}}$ for $t=0,\ldots,n$. The partition $\{{\mathscr{Z}}_{t}\}_{1\leq t\leq n}$, and the densities $\{f_{t}\}_{1\leq t\leq n}$ and $\{{g}_{t}\}_{1\leq t\leq n}$, are defined as before via (1), (2), and (3), respectively. Similarly, one can define a sequence of sub level-sets of $S$; in this paper we use the latter for some cases and whenever appropriate.

Letting $z_{t}\stackrel{{\scriptstyle\tiny\mathrm{def}}}{{=}}\mathbb{E}_{f}\left[\varphi\left(\mathbf{X}\right)\mid\mathbf{X}\in{\mathscr{Z}}_{t}\right]\mathbb{P}_{f}\left(\mathbf{X}\in{\mathscr{Z}}_{t}\right)$ for $t=1,\ldots,n$, and combining (1) with the law of total probability, we arrive at

[TABLE]

The SSA proceeds with the construction of estimators $\widehat{Z}_{t}$ for $z_{t}$ for $t=1,\ldots,n$ and, as soon as these are available, we can use (4) to deliver the SSA estimator for $z$, namely $\widehat{Z}=\sum_{t=1}^{n}{\widehat{Z}_{t}}$.

For $1\leq t\leq n$, let $\varphi_{t}\stackrel{{\scriptstyle\tiny\mathrm{def}}}{{=}}\mathbb{E}_{f}\left[\varphi\left(\mathbf{X}\right)\mid\mathbf{X}\in{\mathscr{Z}}_{t}\right]$, $p_{t}\stackrel{{\scriptstyle\tiny\mathrm{def}}}{{=}}\mathbb{P}_{f}\left(\mathbf{X}\in{\mathscr{Z}}_{t}\right)$, and let $\widehat{\Phi}_{t}$ and $\widehat{P}_{t}$ be estimators of $\varphi_{t}$ and $p_{t}$, respectively. We define $\widehat{Z}_{t}=\widehat{\Phi}_{t}\,\widehat{P}_{t}$, and recall that, under the multilevel splitting framework, we obtain the sets $\{{\cal X}_{t}\}_{1\leq t\leq n}$, and $\{{\cal Z}_{t}\}_{1\leq t\leq n}$. These sets are sufficient to obtain unbiased estimators $\{\widehat{\Phi}_{t}\}_{1\leq t\leq n}$ and $\{\widehat{P}_{t}\}_{1\leq t\leq n}$, in the following way.

We define $\widehat{\Phi}_{t}$ to be the (unbiased) Crude Monte Carlo (CMC) estimator of $\varphi_{t}$, that is,

[TABLE] 2. 2.

The estimator $\widehat{P}_{t}$ is defined similar to the one used in the Generalized Splitting (GS) algorithm of [4]. In particular, the GS product estimator is defined as follows. Define the level entrance probabilities $r_{0}\stackrel{{\scriptstyle\tiny\mathrm{def}}}{{=}}1$, $r_{t}\stackrel{{\scriptstyle\tiny\mathrm{def}}}{{=}}\mathbb{P}_{f}\left(\mathbf{X}\in\mathscr{X}_{t}\mid\mathbf{X}\in\mathscr{X}_{t-1}\right)$ for $t=1,\ldots,n$, and note that $\mathbb{P}_{f}\left(\mathbf{X}\in\mathscr{X}_{t}\right)=\prod_{i=0}^{t}{r_{i}}$. Then, for $t=1,\ldots,n$, it holds that

[TABLE]

This suggests the estimator $\widehat{P}_{t}=(1-\widehat{R}_{t})\prod_{i=0}^{t-1}\widehat{R}_{i}$ for $p_{t}$, where $\widehat{R}_{0}\stackrel{{\scriptstyle\tiny\mathrm{def}}}{{=}}1$ and $\widehat{R}_{t}=\frac{|{\cal Y}_{t}|}{|{\cal X}_{t}|}=\frac{|{\cal X}_{t}\setminus{\cal Z}_{t}|}{|{\cal X}_{t}|}$ for all $t=1,\ldots,n$.

In practice, obtaining the $\{{\cal X}_{t}\}_{1\leq t\leq n}$ and $\{{\cal Z}_{t}\}_{1\leq t\leq n}$ sets requires sampling from the conditional pdfs in (2) and (3). However, for many specific applications, designing such a procedure can be extremely challenging. Nevertheless, we can use the $\mathcal{Y}_{t}={\cal X}_{t}\setminus{\cal Z}_{t}$ set from iteration $t$ to sample $\mathcal{X}_{t+1}$ from $f_{t+1}$ for each $t=1,\ldots,n-1$, via MCMC. In particular, the particles from the $\mathcal{Y}_{t}$ set can be “split”, in order to construct the desired set $\mathcal{X}_{t+1}$ for the next iteration, using a Markov transition kernel $\kappa_{t+1}\left(\cdot\mid\cdot\right)$ whose stationary pdf is $f_{t+1}$, for each $t=1,\ldots,n-1$. Algorithm 2 summarizes the general procedure for the SSA. The particular splitting step described in Algorithm 2 is a popular choice [4, 3, 37], especially for hard problems with unknown convergence behavior of the corresponding Markov chain.

Remark 1** (Independent SSA)**

Obviously, the SSA above produces dependent samples in sets ${\cal X}_{t}$ and ${\cal Z}_{t}$. It will be convenient to also consider a version of SSA where the samples are independent. This can for example be achieved by simulating multiple independent runs of SSA, until $N$ realizations $\mathbf{X}\in{\cal X}_{t}$ have been produced. Of course this Independent SSA (ISSA) variant is much less efficient than the original algorithm, but the point is that the computational effort (the run time complexity) remains polynomial — $\mathcal{O}(C^{2}n)$ versus $\mathcal{O}(Cn)$, where $C$ is the complexity of SSA for each level. Hence, from a complexity point of view we may as well consider the independent SSA variant. This is what will be done in Theorem 2.

Theorem 1** (Unbiased estimator)**

Algorithm 2 outputs an unbiased estimator; that is, it holds that $\mathbb{E}\left[\widehat{Z}\right]=\mathbb{E}_{f}\left[\varphi\left(\mathbf{X}\right)\right]=z.$

Proof 1

See Appendix A.

We next proceed with a clarification for a few remaining practical issues regarding the SSA.

Determining the SSA levels. It is often difficult to make an educated guess how to set the values of the level thresholds. However, the SSA requires the values of $\{\gamma_{t}\}_{1\leq t\leq n}$ to be known in advance, in order to ensure that the estimator is unbiased. To resolve this issue, we perform a single pilot run of Algorithm 2 using a so-called rarity parameter $0<\rho<1$. In particular, given samples from an ${\cal X}_{t}$ set, we take the $\rho\,|{\cal X}_{t}|$ performance quantile as the value of the corresponding level $\gamma_{t}$, and form the next level set. Such a pilot run helps to establish a set of threshold values adapted to the specific problem. After the completion of the pilot run we simply continue with a regular execution of Algorithm 2 using the level threshold values observed in the pilot run.

Controlling the SSA error. A common practice when working with a Monte Carlo algorithm that outputs an unbiased estimator, is to run it for $R$ independent replications to obtain $\widehat{Z}^{(1)},\ldots,\widehat{Z}^{(R)}$, and report the average value. Thus, for a final estimator, we take

[TABLE]

To measure the quality of the SSA output, we use the estimator’s relative error (RE), which is defined by

[TABLE]

As the variance and expectation of the estimator are not known explicitly, we report an estimate of the relative error by estimating both terms from the result of the $R$ runs.

3 Complexity

In this section, we present a theoretical complexity analysis of the Independent SSA variant, discussed in Remark 1. Our focus is on establishing time-complexity results, and thus our style of analysis is similar to that used for approximate counting algorithms. For an extensive overview, we refer to [29, Chapter 10]. We begin with a definition of a randomized algorithm’s efficiency.

Definition 1** ([29])**

A randomized algorithm gives an $(\varepsilon,\delta)$-approximation for the value $z$ if the output $\widehat{Z}$ of the algorithm satisfies

[TABLE]

With the above definition in mind, we now aim to specify the sufficient conditions for the SSA to provide an $(\varepsilon,\delta)$-approximation to $z$. A key component in our analysis is to construct a Markov chain $\left\{X^{(m)}_{t},\,m\geq 0\right\}$ with stationary pdf $f_{t}$ (defined in (2)), for all $1\leq t\leq n$, and to consider the speed of convergence of the distribution of $X^{(m)}_{t}$ as $m$ increases. Let $\mu_{t}$ be the probability distribution corresponding to $f_{t}$, so

[TABLE]

for all Borel sets $A$, where $\lambda$ is some base measure, such as the Lebesgue or counting measure. To proceed, we have

[TABLE]

for the $\tau$-step transition law of the Markov chain. Consider the total variation distance between $\kappa_{t}^{\tau}(\cdot\mid\mathbf{x})$ and $\mu_{t}$, defined as:

[TABLE]

An essential ingredient of our analysis is the so-called mixing time (see [36] and [27] for an extensive overview), which is defined as

[TABLE]

Let $\widehat{\mu}_{t}=\kappa_{t}^{\tau}(\cdot\mid\mathbf{x}\,)$ be the SSA sampling distribution at steps $1\leq t\leq n$ where, for simplicity, we suppress $\mathbf{x}$ in the notation of $\widehat{\mu}_{t}$.

Finally, similar to $\mu_{t}$ and $\widehat{\mu}_{t}$, let $\nu_{t}$ be the probability distribution corresponding to the pdf ${g}_{t}$ (defined in (3)), and let $\widehat{\nu}_{t}$ be the SSA sampling distribution, for all $1\leq t\leq n$.

Theorem 2 details the main efficiency result for the ISSA. We reiterate that from a complexity point of view the independence setting does not impose a theoretical limitation, as the run time complexity remains polynomial. An advantage is that by using ISSA we can engage powerful concentration inequalities [6, 20].

Theorem 2** (Complexity of the ISSA)**

Let $\varphi$ be a strictly positive real-valued function, $a_{t}=\min_{\mathbf{x}\in\mathscr{Z}_{t}}{\{\varphi(\mathbf{x})\}}$, $b_{t}=\max_{\mathbf{x}\in\mathscr{Z}_{t}}{\{\varphi(\mathbf{x})\}}$, and $\underline{r}_{t}=\min{\{r_{t},1-r_{t}\}}$ for $1\leq t\leq n$. Let $\widehat{\mu}_{t}$ and $\widehat{\nu}_{t}$ be the ISSA sampling distributions at steps $1\leq t\leq n$, for ${\cal X}_{t}$ and ${\cal Z}_{t},respectively$. Then, the ISSA gives an $(\varepsilon,\delta)$-approximation to $z=\mathbb{E}_{f}[\varphi(\mathbf{X})]$, provided that for all $1\leq t\leq n$ the following holds.

$\displaystyle\|\widehat{\mu}_{t}-\mu_{t}\|_{\mathrm{TV}}\leq\frac{\varepsilon\,\underline{r}_{t}}{32n}\quad\text{and}\quad|{\cal X}_{t}|\geq\frac{3072\,n^{2}\ln(4n^{2}/\delta)}{\varepsilon^{2}\underline{r}_{t}^{2}}.$** 2. 2.

$\|\widehat{\nu}_{t}-\nu_{t}\|_{\mathrm{TV}}\leq\frac{\varepsilon\,a_{t}}{16(b_{t}-a_{t})}\quad\text{and}\quad|{\cal Z}_{t}|\geq\frac{128(b_{t}-a_{t})^{2}\ln(4n/\delta)}{\varepsilon^{2}a_{t}^{2}}$.

Proof 2

See Appendix A.

In some cases, the distributions of the states in $\mathscr{X}_{t}$ and $\mathscr{Z}_{t}$ generated by Markov chain defined by the kernel $\kappa_{t}^{\tau}$, approach the target distributions $\mu_{t}$ and $\nu_{t}$ very fast. This occurs for example when there exists a polynomial in $n$ (denoted by $\mathcal{P}(n)$), such that the mixing time [27] is bounded by $\mathcal{O}(\mathcal{P}(n))$, $(b_{t}-a_{t})^{2}/a^{2}=\mathcal{O}(\mathcal{P}(n))$, and $\underline{r}_{t}=\mathcal{O}(1/\mathcal{P}(n))$ for all $1\leq t\leq n$. In this case, the ISSA becomes a fully polynomial randomized approximation scheme (FPRAS) [29]. In particular, the ISSA results in a desired $(\varepsilon,\delta)$-approximation to $z=\mathbb{E}_{f}[\varphi(\mathbf{X})]$ with running time bounded by a polynomial in $n$, $\varepsilon^{-1}$, and $\ln(\delta^{-1})$. Finally, it is important to note that an FPRAS algorithm for such problems is essentially the best result one can hope to achieve [22].

We next illustrate the use of Theorem 2 with an example of a difficult problem for which the ISSA provides an FPRAS.

4 FPRAS for the weighted component model

Consider a system of $k$ components. Each component $i$ generates a specific amount of benefit, which is given by a positive real number $w_{i}$, $i=1,\ldots,k$. In addition, each component can be operational or not.

Let $\mathbf{w}=(w_{1},\ldots,w_{k})^{\top}$ be the column vector of component weights (benefits), and $\mathbf{x}=(x_{1},\ldots,x_{k})^{\top}$ be a binary column vector, where $x_{i}$ indicates the $i$th component’s operational status for $1\leq i\leq k$. That is, if the component $i$ is operational $x_{i}=1$, and $x_{i}=0$ if it is not. Under this setting, we define the system performance as

[TABLE]

We further assume that all elements are independent of each other, and that each element is operational with probability $1/2$ at any given time. For the above system definition, we might be interested in the following questions.

Conditional expectation estimation. Given a minimal threshold performance $\gamma\leq\sum_{i=1}^{k}w_{i}$, what is the expected system performance? That is to say, we are interested in the calculation of

[TABLE]

where $\mathbf{X}$ is a $k$-dimensional binary vector generated uniformly at random from the $\{0,1\}^{k}$ set. This setting appears (in a more general form), in a portfolio credit risk analysis [16], and will be discussed in Section 5. 2. 2.

Tail probability estimation [1]. Given the minimal threshold performance $\gamma$, what is the probability that the overall system performance is smaller than $\gamma$? In other words, we are interested in calculating

[TABLE]

The above problems are both difficult, since a uniform generation of $\mathbf{X}\in\{0,1\}^{k}$, such that $\mathbf{w}^{\top}\mathbf{X}\leq\gamma$, corresponds to the knapsack problem, which belongs to #P complexity class  [44, 31]. In this section, we show how one can construct an FPRAS for both problems under the mild condition that the difference between the minimal and the maximal weight in the $\mathbf{w}$ vector is not large. This section’s main result is summarized next.

Proposition 1

Given a weighted component model with $k$ weights, $\mathbf{w}=(w_{1},\ldots,w_{k})$ and a threshold $\gamma$, let $\underline{w}=\min\{\mathbf{w}\}$, $\overline{w}=\max\{\mathbf{w}\}$. Then, provided that $\overline{w}=\mathcal{O}\left(\mathcal{P}(k)\right)\underline{w}$, there exists an FPRAS for the estimation of both (5) and (6).

Prior to stating the proof of Proposition 1, define

[TABLE]

and let $\mu_{b}$ be the uniform distribution on the $\mathscr{X}_{b}$ set. [31] introduce an MCMC algorithm that is capable of sampling from the $\mathscr{X}_{b}$ set almost uniformly at random. In particular, this algorithm can sample $\mathbf{X}\sim\widehat{\mu}_{b}$, such that $\left\|\widehat{\mu}_{b}-\mu_{b}\right\|_{\mathrm{TV}}\leq\varepsilon$. Moreover, the authors show that their Markov chain mixes rapidly, and in particular, that its mixing time is polynomial in $k$ and is given by $\tau_{\mathrm{mix}}(\varepsilon)=\mathcal{O}\left(k^{9/2+\varepsilon}\right)$. Consequentially, the sampling from $\widehat{\mu}_{b}$ can be performed in $\mathcal{O}(\mathcal{P}(k))$ time for any $\varepsilon>0$.

The proof of Proposition 1 depends on the following technical lemma.

Lemma 1

Let $X\sim\widehat{\pi}$ be a strictly positive univariate random variable such that $a\leq X\leq b$, and let $X_{1},\ldots,X_{m}$ be its independent realizations. Then, provided that

[TABLE]

it holds that:

[TABLE]

Proof 3

See Appendix A.

We next proceed with the proof of Proposition 1 which is divided into two parts. The first for the conditional expectation estimation and the second for the tail probability evaluation.

Proof 4** (Proposition 1: FPRAS for (5))**

With the powerful result of [31] in hand, one can achieve a straightforward development of an FPRAS for the conditional expectation estimation problem. The proof follows immediately from Lemma 1. In particular, all we need to do in order to achieve an $(\varepsilon,\delta)$ approximation to (5) is to generate

[TABLE]

samples from $\widehat{\mu}_{\gamma}$, such that

[TABLE]

Recall that the mixing time is polynomial in $k$, and note that the number of samples $m$ is also polynomial in $k$, thus the proof is complete, since

[TABLE]

Proof 5** (Proposition 1: FPRAS for (6))**

In order to put this problem into the setting of Theorem 2 and achieve an FPRAS, a careful definition of the corresponding level sets is essential. In particular, the number of levels should be polynomial in $k$, and the level entrance probabilities $\{r_{t}\}_{1\leq t\leq n}$, should not be too small. Fix

[TABLE]

to be the number of levels, and set $\gamma_{t}=\gamma+(n-t)\,\underline{w}$ for $t=0,\ldots,n$. For general $\gamma$ it holds that

[TABLE]

where the last equality follows from the fact that there is only one vector $\mathbf{x}=(1,1,\ldots,1)$ for which $\gamma+n\,\underline{w}<S(\mathbf{x})\leq\sum_{i=1}^{k}{w_{i}}$. That is, it is sufficient to develop an efficient approximation to $(*)$ only, since the rest are constants.

We continue by defining the sets ${\mathscr{X}}={\mathscr{X}}_{\gamma_{0}}\supseteq\cdots\supseteq{\mathscr{X}}_{\gamma_{n}}$ via (7), and by noting that for this particular problem, our aim is to find $\mathbb{P}(\mathbf{X}\in{\mathscr{X}}_{\gamma_{n}})$, so the corresponding estimator simplifies to (see Section 2.2 (2)),

[TABLE]

In order to show that the algorithm provides an FPRAS, we will need to justify only condition (1) of Theorem 2, which is sufficient in our case because we are dealing with an indicator integrand. Recall that the formal requirement is

[TABLE]

where $\mu_{t}$ is the uniform distribution on ${\mathscr{X}}_{\gamma_{t}}$ for $t=0,\ldots,n$, and each sample in ${\cal X}_{t}$ is distributed according to $\widehat{\mu}_{t}$. Finally, the FPRAS result is established by noting that the following holds.

From Lemma 6 in the Appendix, we have that $\underline{r}_{t}\geq\frac{1}{k+1}$ for $1\leq t\leq n$. 2. 2.

The sampling from $\widehat{\mu}_{t}$ can be performed in polynomial (in $k$) time **[31]**. 3. 3.

The number of levels $n$ (and thus the required sample size $\{|{\cal X}_{t}|\}_{1\leq t\leq n}$) is polynomial in $k$, since

[TABLE]

Unfortunately, for many problems, an analytical result such as the one obtained in this section is not always possible to achieve. The aim of the following numerical section is to demonstrate that the SSA is capable of handling hard problems in the absence of theoretical performance.

5 Numerical experiments

5.1 Portfolio credit risk

We consider a portfolio credit risk setting [16]. Given a portfolio of $k$ assets, the portfolio loss $L$ is the random variable

[TABLE]

where $l_{i}$ is the risk of asset $i\in\{1,\ldots,k\}$, and $X_{i}$ is an indicator random variable that models the default of asset $i$. Under this setting (and similar to Section 4), one is generally interested in the following.

Conditional Value at Risk (CVaR). Given a threshold (value at risk) $v$, calculate the conditional value at risk $c=\mathbb{E}[L\mid L\geq v]$. 2. 2.

Tail probability estimation. Given the value at risk, calculate the tail probability $\mathbb{P}(L\geq v)=\mathbb{E}\left[\mathbbm{1}{\{L\geq v\}}\right]$.

The SSA can be applied to both problems as follows. For tail probability estimation, we simply set $\varphi{}(\mathbf{x})=\mathbbm{1}{\{\sum_{i=1}^{k}{l_{i}\;x_{i}}\geq v\}}$. For conditional expectation, we set $\varphi{}(\mathbf{x})=\sum_{i=1}^{k}{l_{i}\;x_{i}}$.

Note that the tail probability estimation (for which the integrand is the indicator function), is a special case of a general integration. Recall that the GS algorithm of [4] works on indicator integrands, and thus GS is a special case of the SSA. Consequently, in this section we will investigate the more interesting (and more general) scenario of estimating an expectation conditional on a rare event.

As our working example, we consider a credit risk in a Normal Copula model and, in particular, a $21$ factor model from [16] with $1,000$ obligors. Thus, the integrals of interest are $1021$ dimensional.

The SSA setting is similar to the weighted component model from Section 4. We define a $k$-dimensional binary vector $\mathbf{x}=(x_{1},\ldots,x_{k})$, for which $x_{i}$ stands for the $i$th asset default ($x_{i}=1$ for default, and [math] otherwise). We take the performance function $S(\mathbf{x})$ to be the loss function (8). Then, the level sets are defined naturally by $\mathscr{X}_{t}=\{\mathbf{x}\;:\;S(\mathbf{x})\geq\gamma_{t}\}$, (see also Section 2.2). In our experiment, we set $\gamma_{0}=0$ and $\gamma_{n}=1+\sum_{i=1}^{k}{l_{i}}$. In order to determine the remaining levels $\gamma_{1},\ldots,\gamma_{n-1}$, we execute a pilot run of Algorithm 2 with $N=1,000$ and $\rho=0.1$. As an MCMC sampler, we use a Hit-and-Run algorithm [26, Chapter 10, Algorithm 10.10], taking a new sample after $50$ transitions.

It is important to note that despite the existence of several algorithms for estimating $c$, the SSA has an interesting feature, that (to the best of our knowledge) is not present in other methods. Namely, one is able to obtain an estimator for several CVaRs via a single SSA run. To see this, consider the estimation of $c_{1},\ldots,c_{s}$ for $s\geq 1$. Suppose that $v_{1}\leq\cdots\leq v_{s}$ and note that it will be sufficient to add these values to the $\{\gamma_{t}\}$ (as additional levels), and retain $s$ copies of $\widehat{P}_{t}$ and $\widehat{Z}_{t}$. In particular, during the SSA execution, we will need to closely follow the $\gamma$ levels, and as soon as we encounter a certain $v_{j}$ for $1\leq j\leq s$, we will start to update the corresponding values of $\widehat{P}_{t}^{(j)}$ and $\widehat{Z}_{t}^{(j)}$, in order to allow the corresponding estimation of $c_{j}=\mathbb{E}[L\mid L\geq v_{j}]$. Despite that such a procedure introduces a dependence between the obtained estimators, they still remain unbiased.

To test the above setting, we perform the experiment with a view to estimate $\{c_{j}\}_{1\leq j\leq 13}$ using the following values at risk:

[TABLE]

The execution of the SSA pilot run with the addition of the desired VaRs (levels) from (9) (marked in bold), yields the following level values of $(\gamma_{0},\ldots,\gamma_{21})$:

[TABLE]

Table 1 summarizes the results obtained by executing $1$ pilot and $5$ regular independent runs of the SSA. For each run, we use the parameter set that was specified for the pilot run ($N=1,000$ and burn-in of $50$). The overall execution time (for all these $R=1+5=6$ independent runs) is 454 seconds. The SSA is very accurate. In particular, we obtain an RE of less than $1\%$ for each $\widehat{c}$ while employing a very modest effort.

The obtained result is especially appealing, since the corresponding estimation problem falls into rare-event setting [15]. That is, a CMC estimator will not be applicable in this case.

5.2 Self-avoiding walks

In this section, we consider random walks of length $n$ on the two-dimensional lattice of integers, starting from the origin. In particular, we are interested in estimating the following quantities:

$c_{n}$: the number of SAWs of length $n$, 2. 2.

$\Delta_{n}$: the expected distance of the final SAW coordinate to the origin.

To put these SAW problems into the SSA framework, define the set of directions, $\mathscr{X}=\{\mathrm{Left,Right,Up,Down}\}^{n}$, and let $f$ be the uniform pdf on $\mathscr{X}$. Let $\xi(\mathbf{x})$ denote the final coordinate of the random walk represented by the directions vector $\mathbf{x}$. We have $c_{n}=\mathbb{E}_{f}\left[\mathbbm{1}{\{\mathbf{X}\text{ is SAW}\}}\right]$ and \Delta_{n}=\mathbb{E}_{f}\left[\|\xi(\mathbf{X})\|\,\,\big{|}\,\mathbbm{1}{\{\mathbf{X}\text{ is SAW}\}}\right].

Next, we let $\mathscr{X}_{t}\subseteq\mathscr{X}$ be the set of all directions vectors that yield a valid self-avoiding walk of length at least $t$, for $0\leq t\leq n$. In addition, we define $\mathscr{Z}_{t}$ to be the set of all directions vectors that yield a self-avoiding walk of length (exactly) $t$, for $1\leq t\leq n$. The above gives the required partition of $\mathscr{X}$. Moreover, the simulation from $f_{t}(\mathbf{x})=f\left(\mathbf{x}\mid\mathbf{x}\in{\mathscr{X}}_{t-1}\right)$, reduces to the uniform selection of the SAW’s direction at time $1\leq t\leq n$.

Our experimental setting for SAWs of lengths $n$ is as follows. We set the sample size of the SSA to be $N_{t}=1000$ for all $t=1,\ldots,n$. In this experiment, we are interested in both the probability that $\mathbf{X}$ lies in $\mathscr{Z}_{n}$, and the expected distance of $\mathbf{X}\in\mathscr{Z}_{n}$ (uniformly selected) to the origin. These give us the required estimators of $c_{n}$ and $\Delta_{n}$, respectively. The leftmost plot of Fig. 1 summarizes a percent error (PE), which is defined by

[TABLE]

where $\widehat{c}_{n}$ stands for the SSA’s estimator of $c_{n}$.

In order to explore the convergence of the SSA estimates to the true quantity of interest, the SSA was executed for a sufficient number of times to obtain $3\%$ and $1\%$ relative error (RE) [39], respectively. The exact $c_{n}$ values for $n=1,\ldots,71$ were taken from [18, 21]); naturally, when we allow a smaller RE, that is, when we increase $R$, the estimator converges to the true value $c_{n}$, as can be observed in leftmost plot of Fig. 1. In addition, the rightmost plot of Fig. 1 shows that regardless of the RE, the required number of independent SSA runs $(R)$ divided by SAW’s length $(n)$, is growing linearly with $n$.

To further investigate the SSA convergence, we consider the two properties of SAW’s. In particular, the following holds [33, 10, 2, 34]

$\mu=\lim_{n\to\infty}c_{n}^{\frac{1}{n}}\in[\underline{\mu},\overline{\mu}]=[2.62002,2.679192495]$. 2. 2.

$\underline{\Delta}_{n}=n^{\frac{1}{4}}\leq\Delta_{n}\approx n^{3/4}=\widetilde{\Delta}_{n}\ll n=\overline{\Delta}_{n}$.

Fig. 2 summarizes our results compared to these bounds. In particular, we run the SSA to achieve the $3\%$ RE (for $\hat{c}_{n}$) for $1\leq n\leq 200$. It can be clearly observed, that the estimator $\widehat{c}_{n}^{1/n}$ converges toward the $[\underline{\mu},\overline{\mu}]$ interval as $n$ grows. It is interesting to note that, at least for small $n$, $\Delta_{n}$ seems to grow at a rate smaller than the suggested $3/4$.

6 Discussion

In this paper we described a general procedure for multi-dimensional integration, the SSA, and applied it to various problems from different research domains. We showed that this method belongs to a very general class of SMC algorithms and developed its theoretical foundation. The proposed SSA is relatively easy to implement and our numerical study indicates that the SSA yields good performance in practice. However, it is important to note that generally speaking, the efficiency of the SSA and similar sequential algorithms is heavily dependent on the mixing time of the corresponding Markov chains that are used for sampling. A rigorous analysis of the mixing time for different problems is thus of great interest. Finally, based on our numerical study, it will be interesting to apply the SSA to other practical problems.

7 Acknowledgements

This work was supported by the Australian Research Council Centre of Excellence for Mathematical & Statistical Frontiers, under grant number CE140100049.

Appendix

Appendix A Technical arguments

Proof of Theorem 1. Recall that ${\mathscr{Z}}_{1},\ldots,{\mathscr{Z}}_{n}$ is a partition of the set ${\mathscr{X}}$, so from the law of total probability we have

[TABLE]

By the linearity of expectation, and since $\widehat{Z}=\sum_{t=1}^{n}{\widehat{Z}_{t}}$, it will be sufficient to show that for all $t=1,\ldots,n$, it holds that

[TABLE]

To see this, we need the following.

Although that the samples in the $\mathcal{Z}_{t}$ set for $1\leq t\leq n$ are not independent due to MCMC and splitting usage, they still have the same distribution; that is, for all $t$, it holds that

[TABLE] 2. 2.

From the unbiasedness of multilevel splitting [26, 4], it holds for all $1\leq t\leq n$ that

[TABLE]

Combining (10) and (11) with a conditioning on the cardinalities of the $\mathcal{Z}_{t}$ sets, we complete the proof with:

[TABLE]

$\Box$

Proof of Lemma 1. Recall that

[TABLE]

for any function $\varphi:\mathbb{R}\rightarrow[a,b]$, (Proposition 3 in [36]). Hence,

[TABLE]

Combining this with the fact that $X\geq a$, we arrive at

[TABLE]

Next, since

[TABLE]

we can apply the [20] inequality, to obtain

[TABLE]

for

[TABLE]

Finally, we complete the proof by combining (12) and (13), to obtain:

[TABLE]

$\Box$

Proof of Theorem 2. The proof of this theorem consists of the following steps.

In Lemma 2, we prove that an existence of an $\left(\varepsilon,\frac{\delta}{n}\right)$-approximation to $\{z_{t}\}_{1\leq t\leq n}$ implies an existence of an $\left(\varepsilon,\delta\right)$-approximation to $z=\sum_{t=1}^{n}{z_{t}}$. 2. 2.

In Lemma 3, we prove that an existence of an $\left(\frac{\varepsilon}{4},\frac{\delta}{2n}\right)$-approximation to $\{\varphi_{t}\}_{1\leq t\leq n}$ and $\{p_{t}\}_{1\leq t\leq n}$ implies an $\left(\varepsilon,\frac{\delta}{n}\right)$-approximation existence to $\{z_{t}\}_{1\leq t\leq n}$. 3. 3.

In Lemmas 1, 4 and 5, we provide the required $\left(\frac{\varepsilon}{4},\frac{\delta}{2n}\right)$-approximations to $\varphi_{t}$ and $p_{t}$ for $1\leq t\leq n$.

Lemma 2

Suppose that for all $t=1,\ldots,n$, an $\left(\varepsilon,\frac{\delta}{n}\right)$-approximation to $z_{t}$ exists. Then,

[TABLE]

Proof 6

From the assumption of the existence of the $\left(\varepsilon,\frac{\delta}{n}\right)$-approximation to $z_{t}$ for each $1\leq t\leq n$, we have

[TABLE]

*By using Boole’s inequality (union bound), we arrive at *

[TABLE]

that is, it holds for all $t=1,\ldots,n$, that

[TABLE]

and hence,

[TABLE]

$\Box$

Lemma 3

Suppose that for all $t=1,\ldots,n$, there exists an $\left(\frac{\varepsilon}{4},\frac{\delta}{2n}\right)$-approximation to $\varphi_{t}$ and $p_{t}$. Then,

[TABLE]

Proof 7

By assuming an existence of $\left(\frac{\varepsilon}{4},\frac{\delta}{2n}\right)$-approximation to $\varphi_{t}$ and $p_{t}$, namely:

[TABLE]

and combining it with the union bound, we arrive at

[TABLE]

where $(*)$ follows from the fact that for any $0<|\varepsilon|<1$ and $n\in\mathbb{N}$ we have

[TABLE]

To see that (14) holds, note that by using exponential inequalities from [5], we have that $\left(1-\frac{\varepsilon/2}{n}\right)^{n}\geq 1-\frac{\varepsilon}{2}\geq 1-\varepsilon$. In addition, it holds that $|\text{e}^{\varepsilon}-1|<7\varepsilon/4$, and hence:

[TABLE]

To complete the proof of Theorem 2, we need to provide $\left(\frac{\varepsilon}{4},\frac{\delta}{2n}\right)$-approximations to both $\varphi_{t}$ and $p_{t}$. However, at this stage we have to take into account a specific splitting strategy, since the SSA sample size bounds depend on the latter. Here we examine the independent setting, for which the samples in each $\{{\cal X}_{t}\}_{1\leq t\leq n}$ set are independent. That is, we use multiple runs of the SSA at each stage $(t=1,\ldots,n)$ of the algorithm execution. See Remark 1 for further details.

Remark 2** (Lemma 1 for binary random variables)**

For a binary random variable $X\in\{0,1\}$, with a known lower bound on its mean, Lemma 1 can be strengthened via the usage of [6] bound instead of the [20] inequality. In particular, the following holds.

Let $X\sim\widehat{\pi}(x)$ be a binary random variable and let $X_{1},\ldots,X_{m}$ be its independent realizations. Then, provided that $\mathbb{E}_{\widehat{\pi}}[X]\geq\mathbb{E}_{\widehat{\pi}}^{\prime}[X]$,

[TABLE]

it holds that:

[TABLE]

The corresponding proof is almost identical to the one presented in Lemma 1. The major difference is the bound on the sample size in (13), which is achieved via the Chernoff bound from [29, Theorem 10.1] instead of Hoeffding’s inequality.

Lemma 4

Suppose that $a_{t}=\min_{\mathbf{x}\in\mathscr{Z}_{t}}{\{\varphi(\mathbf{x})\}}$, $b_{t}=\max_{\mathbf{x}\in\mathscr{Z}_{t}}{\{\varphi(\mathbf{x})\}}$ for all $t=1,\ldots,n$. Then, provided that the samples in the ${\cal Z}_{t}$ set are independent, and are distributed according to $\widehat{\nu}_{t}$ such that

[TABLE]

then $\widehat{\Phi}_{t}={|{\cal Z}_{t}|^{-1}}\sum_{\mathbf{X}\in{\cal Z}_{t}}{\varphi(\mathbf{X})}$ is an $\left(\frac{\varepsilon}{4},\frac{\delta}{2n}\right)$-approximation to $\varphi_{t}$.

Proof 8

The proof is an immediate consequence of Lemma 1. In particular, note that

[TABLE]

and that

[TABLE]

Lemma 5

Suppose that the samples in the ${\cal X}_{t}$ set are independent, and are distributed according to $\widehat{\mu}_{t}$, such that

[TABLE]

where $\underline{r}_{t}=\min{\{r_{t},1-r_{t}\}}$ for $1\leq t\leq n$. Then, $\widehat{P}_{t}$ is an $\left(\frac{\varepsilon}{4},\frac{\delta}{2n}\right)$-approximation to $p_{t}$.

Proof 9

Recall that $\widehat{P}_{t}=\left(1-\widehat{R}_{t}\right)\prod_{j=0}^{t-1}{\widehat{R}_{j}}$ for $t=1,\ldots,n$. Again, by combining the union bound with (14), we conclude that the desired approximation to $p_{t}$ can be obtained by deriving the $\left(\frac{\varepsilon}{8n},\frac{\delta}{2n^{2}}\right)$-approximations for each $r_{t}$ and $1-r_{t}$. In this case, the probability that for all $t={1,\ldots,n}$, $\widehat{R}_{t}/r_{t}$ satisfies $1-\varepsilon/8n\leq{\widehat{R}_{t}}/r_{t}\leq 1+\varepsilon/8n$ is at least $1-\delta/2n^{2}$. The same holds for $({1-\widehat{R}_{t}})/(1-r_{t})$, and thus, we arrive at:

[TABLE]

The bounds for each $\widehat{R}_{t}$ and $(1-\widehat{R}_{t})$ are easily achieved via Remark 2. In particular, it is not very hard to verify that in order to get an $\left(\frac{\varepsilon}{8n},\frac{\delta}{2n^{2}}\right)$-approximation, it is sufficient to take

[TABLE]

and

[TABLE]

Lemma 6

Suppose without loss of generality that $\mathbf{w}=(w_{1},\ldots,w_{k})$ satisfies $w_{1}\leq w_{2}\leq\cdots\leq w_{k}$, that is $\underline{w}=w_{1}$. Then, for $\mathscr{X}_{b}$ and $\mathscr{X}_{b-w_{1}}$ sets defined via (7), it holds that:

[TABLE]

Proof 10

For any $b\in\mathbb{R}$ and $\mathbf{x}=(x_{1},\ldots,x_{k})$, define a partition of $\mathscr{X}_{b}$ via

[TABLE]

Then, the following holds.

For any $\mathbf{x}\in\mathscr{X}_{b}^{(w_{1})}$, replace $x_{1}=1$ with $x_{1}=0$, and note that the resulting vector is in $\mathscr{X}_{b-w_{1}}^{(-w_{1})}$ set, since its performance is at most $b-w_{1}$, that is $\left|\mathscr{X}_{b}^{(w_{1})}\right|\leq\left|\mathscr{X}_{b-w_{1}}^{(-w_{1})}\right|$. Similarly, for any $\mathbf{x}\in\mathscr{X}_{b-w_{1}}^{(-w_{1})}$, setting $x_{1}=1$ instead of $x_{1}=0$, results in a vector which belongs to the $\mathscr{X}_{b}^{(w_{1})}$ set. That is:

[TABLE] 2. 2.

For any $\mathbf{x}\in\mathscr{X}_{b}^{(w_{1})}$, replace $x_{1}=1$ with $x_{1}=0$ and note that the resulting vector is now in the $\mathscr{X}_{b}^{(-w_{1})}$ set, that is $|\mathscr{X}_{b}^{(w_{1})}|\leq|\mathscr{X}_{b}^{(-w_{1})}|$. In addition, for any $\mathbf{x}\in\mathscr{X}_{b}^{(-w_{1})}$, there are at most $k-1$ possibilities to replace $\mathbf{x}$’s non-zero entry with zero and set $x_{1}=1$, such that the result will be in the $\mathscr{X}_{b}^{(w_{1})}$ set. That is, $\left|\mathscr{X}_{b}^{(-w_{1})}\right|\leq(k-1)\left|\mathscr{X}_{b}^{(w_{1})}\right|+1$, (where $+1$ stands for the vector of zeros), and we arrive at

[TABLE]

Combining (15) and (16), we complete the proof by noting that

[TABLE]

$\Box$