G Matthews

Michael Arbel

Danilo J Rezende

Arnaud Doucet
email: doucet@google.com>.

Continual Repeated

Alexander G D G Matthews
email: <alexmatthews@google.com>

Continual Repeated Annealed Flow Transport Monte Carlo

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

There are few algorithmic problems richer or more fundamental than the task of drawing samples from an unnormalized distribution and approximating its normalizing constant. In this paper we combine two important methods for this task specifically Sequential Monte Carlo (SMC) [START_REF] Del Moral | Sequential Monte Carlo samplers[END_REF] and variational inference with Normalizing flows (NFs) [START_REF] Jimenez | Variational inference with normalizing flows[END_REF].

The first of these components, SMC samplers, are an importance sampling based algorithm. They are a principled way to combine annealing, resampling and Markov chain Monte Carlo (MCMC) and as such enjoy continuing popularity [START_REF] Dai | An invitation to sequential Monte Carlo samplers[END_REF]. One potential disadvantage of SMC samplers is that at each successive temperature step in a SMC sampler there is importance sampling-albeit between adjacent distributions. For high dimensional target distributions this can lead either to estimators with high error or require many temperatures. Another potential disadvantage is that for finite particle numbers, SMC estimates of expectation w.r.t. the target distribution are biased. This can be mitigated by using the SMC sampler inside a Particle MCMC outer loop [START_REF] Andrieu | Particle Markov chain Monte Carlo methods[END_REF] at the cost of adding repeated SMC sampler calls. It is particularly desirable in this Particle MCMC context for the base SMC sampler to be fast and output low variance estimates of the normalizing constant.

The other component of our method is variational inference with NFs. Flow methods learn a differentiable invertible transformation (diffeomorphism) typically with a tractable Jacobian [START_REF] Papamakarios | Normalizing flows for probabilistic modeling and inference[END_REF]. The learnt mappings are both flexible and fast to evaluate but there are still some challenges around using flows. The tractability of training with reverse Kullback-Leibler divergence comes at the cost of its well known mode seeking behaviour. Further placing the full burden of modelling on a parameterized flow can lead to high memory usage in contrast to non-parametric sampling algorithms. Finally, as diffeomorphisms, flows preserve topological properties of their input. This property can represent a challenge whenever there is topological mismatch between the, typically, simple base distribution and the complex target [START_REF] Cornish | Relaxing bijectivity constraints with continuously indexed normalising flows[END_REF].

Methods combining annealed samplers with normalizing flows have the potential to fix the limitations of each component part. The use of flows can reduce the variance from importance sampling. The MCMC steps can reduce the topological and representational burden on the flows. Used correctly, annealing can reduce the mode seeking behaviour of the variational objective. Resampling can be used to focus computation on promising samples. The idea of combining fixed transformations with annealing goes back as far as [START_REF] Vaikuntanathan | Escorted free energy simulations[END_REF]. It is the training of the flows that still presents a challenge and the reason that such combined approaches have not met their full potential. Two recent papers target this problem [START_REF] Wu | Stochastic normalizing flows[END_REF][START_REF] Arbel | Annealed flow transport Monte Carlo[END_REF], but as we shall relate, neither is entirely satisfactory or as well suited to the task as the algorithm we propose.

CRAFT is a method for estimating an SMC sampler augmented with interleaved normalizing flows. Unlike many other approaches in this area that incorporate sampling ideas, it is able to cope with full MCMC steps with Metropolis accept/reject corrections and SMC resampling steps. The CRAFT objective is not a standard variational bound. Rather it uses a KL divergence for each transition between temperatures. We find that as well as outperforming other methods for estimating the flows in practice (Section 4.1 and Section 4.2), our learnt sampler is useful as a fast inner loop for a Particle MCMC sampler (Section 4.3).

Method

In this section we describe and motivate both the sampling method and how to learn the required flows. Section 2.1 describes the method for fixed flows, which is equivalent to an SMC sampler with added fixed normalizing flows. In Section 2.2, we discuss conventional evidence lower bound (ELBO) methods for learning the flows and their limitations. Finally, we present the CRAFT training method in Section 2.3, and also explain how we improve over the AFT method on which CRAFT builds, by solving what we call the sample replenishment problem.

Sampler for fixed flows

We start by describing the sampler for fixed normalizing flows. In this case, it corresponds to a generalization of the standard SMC sampler which adds deterministic normalizing flows to the usual steps in SMC. This generalization reduces to a standard SMC sampler when the normalizing flows are the identity. [START_REF] Arbel | Annealed flow transport Monte Carlo[END_REF] give a detailed history. The algorithm is described in Algorithm 3 which sequentially calls Algorithm 1.

We now describe the steps and derivation in more detail. We consider a sequence of distributions (π k (x))

K k=0 on R M where π k (x)= γ k (x)
Z k and γ k (x) can be evaluated pointwise. The final distribution π K (x) is the distribution with unknown normalizing constant Z K that we wish to approximate. For the initial distribution π 0 (x) we assume that we can draw exact samples and that the normalizing constant Z 0 = 1 without loss of generality. The intermediate distributions enable us to transition smoothly (or in physics parlance 'anneal') from the tractable π 0 (x) to our goal π K (x). While they will sometimes be of interest themselves, they will often only be used to help us construct approximations. The method returns N weighted particles

X i K , W i K N i=1
which can be used to provide an unbiased estimate of the normalizing constant Z K and consistent estimates of expectations under the target π K .

To build up intuition, consider a forward sampling process η producing a sequence

(X k ) K k=0 defined on R M ×(K+1)
such that the final sample X K approximates π K . The process starts with an initial sample X 0 drawn from some distribution π 0 that is successively transformed using an interleaved sequence of K NFs (T k ) K k=1 and Markov transition kernels

(K k ) K k=1 with invariant distributions (π k ) K k=1 : Y 1 =T 1 (X 0 ) X 1 ∼K 1 (Y 1 , •) • • • Y K =T K (X K-1 ) X K ∼K K (Y K , •) X 0 ∼π 0 X 1 X K-1 X K
Intuitively, the purpose of each flow T k is to 'transport' samples from a density π k-1 corresponding to an annealing temperature k-1 towards the next density π k at temperature k. Successive densities are usually similar, making the flow easier to learn compared to a flow directly transporting π 0 towards the target π K . Since it is generally hard to find a flow T k that perfectly transports between successive densities, the sampler employs two correction mechanisms: (1) the Markov transition kernel K k with invariant density π k further diffuses the samples towards π k and (2) Importance Sampling re-weights the samples to correct for any mismatch between π k and the density T # k π k-1 obtained by transporting π k-1 using the flow T k . Importance Sampling. Ultimately, we are interested in re-weighting the final sample X K to approximate the target π K . This could be theoretically achieved by evaluating the target to marginal ratio π K /η K , where η K is obtained by integrating the proposal distribution η over all previous variables (X k ) K-1 k=0 . To avoid computing such intractable integrals, we follow the standard approach in SMC and Annealed Importance Sampling (AIS) [START_REF] Radford | Annealed importance sampling[END_REF] of computing importance weights between the whole forward process η and an augmented target π admitting π K as marginal at time K. While there are multiple choices for the augmented target π, as discussed in [START_REF] Del Moral | Sequential Monte Carlo samplers[END_REF], our choice of π is made for tractability of the importance weights and corresponds to a formal backward process starting by sampling exactly from π K and then sequentially generating samples backwards in time using the reversal transformations of each flow T k and MCMC kernel K k . The reversal of each flow transformation is its inverse while the reversal of each forward Markov kernel is the Markov kernel Kk satisfying Kk (x , x). When choosing a reversible Markov kernel for the forward kernel K k , the reversal kernel is equal to the forward kernel, i.e. Kk =K k . For our choice of augmented target, the unnormalized importance weights take the form:

π k (x)K k (x, x )=π k (x )
w K (x 0:K-1 ) = K k=1 G k (x k-1 ), (1a)
G k (x k-1 ) := γ k (T k (x k-1 )) γ k-1 (x k-1 ) |∇T k (x k-1 )|. ( 1b
)
When the flow is the identity, equation (1b) reduces to the ratio of successive densities, which is the standard annealed importance sampling expression. For non-identity flows the T k dependent terms correct for the use of the flow. See Appendix A for a precise mathematical description of η, π and the derivation of the weights (1a)-(1b).

The algorithm is then implemented sequentially to maintain a set of N particles (X

(i) k , W i k ) N i=1
consisting in pairs of samples and corresponding normalized IS weights at each time k where the samples are obtained using the forward generating process η. The unnormalized IS weights w i k are computed recursively using

w i k =W i k-1 G k (X i k-1
), thus allowing to compute the normalized IS weights W i k by normalizing over the set of particles

W i k =w i k / N j=1 w j k . A consistent approximate expectation of a function f un- der the target π k is then given by N i=1 W i k f (X (i) k ). More- over, an unbiased estimate Z N k of the normalizing constant Z k is obtained sequentially from the previous one Z N k-1 us- ing the update Z N k =Z N k-1 ( N i=1 w i k ).
If the NFs transport perfectly between the temperatures, this algorithm enjoys the property that the normalizing constant estimate Z N K has a zero variance and the corresponding approximate expectations are unbiased with a variance corresponding to the one of the true target.

Resampling. Up to this point, we have described AIS with additional fixed normalizing flows. We now go a step further to incorporate resampling, a key ingredient of SMC methods that has proved to be very beneficial [START_REF] Arbel | Annealed flow transport Monte Carlo[END_REF][START_REF] Chopin | A sequential particle filter method for static models[END_REF][START_REF] Del Moral | Sequential Monte Carlo samplers[END_REF][START_REF] Hukushima | Population annealing and its application to a spin glass[END_REF]. Resampling consists in randomly selecting a particle X i k with probability W i k and repeating the operation N times to construct a set of N particles approximately distributed according to π k . This operation is equivalent to assigning a number of offspring N i k to each particle X i k and associating a uniform weight of 1/N to each offspring. The vector of all offspring (N i k ) N i=1 is then drawn from a multinomial distribution with weights

(W i k ) N i=1 under the constraint that N i=1 N i k =N . The ex- pectation under π k of a function f is then approximated using N i=1 N i k N f (X i k )
which is also an unbiased estimate of the weighted sum

N i=1 W i k f (X i k ).
Hence, resampling refocuses computational effort on promising particles (the ones falling in high density regions of π k ) whilst preserving the key analytic properties of the algorithm. However, to avoid the additional variance coming from the multinomial distribution, we only use it when the effective sample size

ESS N k :=( N i=1 (W i k ) 2
) -1 of the particles falls beyond some predefined proportion A∈[1/N, 1) (we use A = 0.3 in all our experiments) of the total number of samples N (Liu and Chen, 1995). Since resampling is followed by the Markov kernel step in each iteration of Algorithm 1 any degeneracy introduced during resampling can be reduced.

Algorithm 1 SMC-NF-step 1: Input:

Approximations (π N k-1 , Z N k-1 ) to (π k-1 , Z k-1 ), normalizing flow T k , unnormal- ized annealed targets γ k-1 and γ k and resampling threshold A ∈ [1/N, 1). 2: Output: Particles at iteration k: π N k = (X i k , W i k ) N i=1 , approximation Z N k to Z k . 3: Transport particles: Y i k = T k (X i k-1
). 4: Compute IS weights:

w i k ← W i k-1 G k (X i k-1 ) // unnormalized W i k ← w i k / N j=1 w j k // normalized 5: Estimate normalizing constant Z k : Z N k ← Z N k-1 ( N i=1 w i k ). 6: Compute effective sample size ESS N k . 7: if ESS N k ≤ N A then 8:
Resample N particles denoted abusively also Y i k according to the weights W i k , then set

W i k = 1 N . 9: end if 10: Sample X i k ∼ K k (Y i k , •). // MCMC 11: Return (π N k , Z N k ).
Algorithm 2 CRAFT-training 1: Input: Initial NFs {T k } 1:N , number of particles N , unnormalized annealed targets {γ k } K k=0 with γ 0 =π 0 and γ K =γ, resampling threshold A ∈ [1/N, 1). 2: Output: Learned flows T k and length J sequence of approximations (π N K , Z N K ) to (π K , Z K ). 3: for j = 1, . . . , J do 4: Sample X i 0 ∼π 0 and set W i 0 = 1 N and Z N 0 =1.

5:

for k = 1, . . . , K do 6: ĥ ← flow-grad T k , π N k-1 using eqn (8).

7:

π N k , Z N k ← SMC-NF-step π N k-1 , Z N k-1 , T k 8:
Update the flow T k using gradient ĥ.

9:

end for 10:

Yield (π N K , Z N K ) and continue for loop. 11: end for 12: Return learned flows {T k } K k=1 .

Algorithm 3 CRAFT-deployment

1: Input: Fixed/trained NFs {T k } K k=1 , number of parti- cles N , unnormalized annealed targets {γ k } K k=0 with γ 0 =π 0 , resampling threshold A ∈ [1/K, 1). 2: Output: Approximations (π N K , Z N K ) to (π, Z). 3: Sample X i 0 ∼ π 0 and set W i 0 = 1 N and Z N 0 = 1. 4: for k = 1, . . . , K do 5: π N k , Z N k ← SMC-NF-step π N k-1 , Z N k-1 , T k 6: end for 7: Return (π N K , Z N K ).

Limitations of ELBO based training

The idea of learning parameters of a forward generating Markov process η to approximate a target distribution π K and its normalizing constant Z K has already been explored in the literature, in particular in the context of VAEs training [START_REF] Salimans | Markov chain Monte Carlo and variational inference: Bridging the gap[END_REF][START_REF] Wu | Stochastic normalizing flows[END_REF][START_REF] Geffner | MCMC variational inference via uncorrected Hamiltonian annealing[END_REF][START_REF] Achille Thin | Monte Carlo variational auto-encoders[END_REF][START_REF] Zhang | Differentiable annealed importance sampling and the perils of gradient noise[END_REF]. In these references, the standard approach consisting of minimizing the KL divergence between the distribution η of the forward generating process and a suitable augmented target π admitting π K as marginal at time K is adopted:

KL[η||π] = log Z K -E η [log w K ] . (2)
However, minimizing this Kullback-Leibler divergence, equivalently maximizing the ELBO objective E η [log w K ], requires differentiating through the forward sampler which can be challenging when MCMC kernels are included as we now discuss in more detail.

Discontinuity of the forward sampler. The most widely used MCMC transition kernels, such as Hamiltonian Monte Carlo (HMC) or Metropolis-adjusted Langevin algorithm (MALA), rely on an Metropolis accept/reject mechanism to ensure invariance. However, these mechanisms are discontinuous functions of the input random variables of the forward sampler. To avoid this issue and enable use of the reparametrization trick, a popular solution is to use approximate MCMC kernels without Metropolis correction which have continuous densities and then use an approximate reverse distribution [START_REF] Salimans | Markov chain Monte Carlo and variational inference: Bridging the gap[END_REF][START_REF] Wu | Stochastic normalizing flows[END_REF][START_REF] Geffner | MCMC variational inference via uncorrected Hamiltonian annealing[END_REF][START_REF] Achille Thin | Monte Carlo variational auto-encoders[END_REF][START_REF] Zhang | Differentiable annealed importance sampling and the perils of gradient noise[END_REF]. The downside of this is the bias that accrues from not having the correction and the difficulty of approximating the reversal. This can mean many slow mixing proposals are required.

It is desirable to use in the forward generating process a combination of standard Metropolis-corrected MCMC samplers (e.g. HMC, MALA) and NFs as used by CRAFT, Stochastic Normalizing Flows (SNFs) [START_REF] Wu | Stochastic normalizing flows[END_REF] and AFT [START_REF] Arbel | Annealed flow transport Monte Carlo[END_REF]. In this case, the proposition below proven in Appendix B and extending the results of [START_REF] Achille Thin | Monte Carlo variational auto-encoders[END_REF] shows that the gradient of the ELBO in the flow parameters is the sum of two terms, one of which is a high variance score term.

Proposition 1. Let U be the set of continuous r.v.s (used to sample proposals) of distribution independent of any parameter and A the set of discrete r.v.s (used for accept/reject steps) during the forward generating process. Let p denote the joint distribution of A, U . The gradient of the ELBO E η [log w K ] in the flow parameters is:

E p [∇ log(w K • Φ) + log(w K • Φ)∇ log p] , (3)
where Φ is a differentiable re-parameterization of the trajectory X 0:K in terms of U for any A, i.e. X 0:K =Φ(U, A).

In particular, when the learned flows (T k ) K k=1 are optimal, i.e. T # k π k-1 = π k , it holds that:

E p [∇ log(w K • Φ)] = E p [log(w K • Φ)∇ log p] = 0.
The first term on the r.h.s. of the expression (3) corresponds to a reparameterization trick term. The second term, however, requires computing a high variance score ∇ log p. As pointed out by Thin et al. ( 2021), algorithms such as SNFs with Metropolis-adjusted MCMC kernels are omitting the second score term when optimizing the ELBO. While this term happens to vanish in the particular case where the flows are optimal, its contribution can, in general, be important especially when initializing the algorithm with suboptimal flows.

Discontinuity of the resampling steps. The literature has mostly focused on scenarios where the forward generating process η is a simple Markov process. SMC forward generating processes as used in CRAFT could be exploited to obtain a tighter ELBO and have been used in the context of state-space models [START_REF] Maddison | Filtering variational objectives[END_REF][START_REF] Tuan | Auto-encoding sequential Monte Carlo[END_REF][START_REF] Christian A Naesseth | Variational sequential Monte Carlo[END_REF]. However, the resampling steps of SMC correspond to sampling discrete distributions and lead to high variance gradient estimates. Omitting them can introduce significant bias [START_REF] Corenflos | Differentiable particle filtering via entropy-regularized optimal transport[END_REF]. These difficulties are further exacerbated in our context since, as discussed earlier, the MCMC steps are not differentiable either.

Emergence of annealing requires mixing. In [START_REF] Wu | Stochastic normalizing flows[END_REF] the ELBO corresponding to equation ( 2) has the form of the standard VI with NFs objective, with additional weight terms entering into w K due to the MCMC steps and annealing schedule. The additional log-weight corrections are

∆ log w SNF MC k = log γ k (x SNF k ) -log γ k (x SNF k ), (4)
where xSNF k is the value of the sample after the MCMC at temperature k has been applied and x SNF k is the value of the sample before it (see Appendix A). This is the only way the annealed densities enter into the Stochastic Normalizing Flow ELBO. Since this additional term is zero in the limit where the Markov kernel does not mix, it follows that the value of the SNF ELBO reduces to the VI with normalizing flows objective in this case. Consequently, for the training objective to be significantly different from the one without MCMC and annealing, this method requires the MCMC kernels to mix sufficiently, which may not be the case in challenging examples.

CRAFT training

Here we describe a method that neatly sidesteps the issues described with the standard ELBO. We take the learning objective to have the following form:

H = K k=1 KL[T # k π k-1 ||π k ] (5)
where the sum runs over each transition between temperatures. The minimum of the objective is zero and this is attained when each flow transports perfectly between successive temperatures. In general, breaking the annealing task down term-by-term means that we can reduce the mode seeking effect of the reverse KL divergence. Each KL term in the sum can be written as follows:

KL[T # k π k-1 ||π k ] = E π k-1 [D k ] + log Z k Z k-1 , (6a)
D k (x) := log γ k-1 (x) γ k (T k (x)) + log |∇ x T k (x)|. ( 6b
)
The expectation over π k-1 on the RHS of the equation is approximated using a streaming SMC approximation of π k-1 which will always be based on the corresponding learnt sampler for the current optimization step:

KL[T # k π k-1 ||π k ]-log Z k Z k-1 ≈ i W i k-1 D k (X i k-1
).

(7) The terms inside the summation require that we can evaluate the output of the normalizing flow T k (X) and the log determinant of the Jacobian log |∇ x T k (x)|. Similarly the gradients with respect to flow parameters θ k of flow T k are approximated:

∂H ∂θ k ≈ i W i k-1 ∂D k (X i k-1 ) ∂θ k . ( 8
)
When the flows transport perfectly, the gradient ∂H ∂θ k will have true value zero and there will be no bias in this particle approximation to it. More generally the gradient estimate will be biased since it is based on a normalized importance estimator. We will analyze this aspect further in Section 3.2, but for now we note that the normal guarantees from the SMC literature imply the gradient is consistent in the number of particles and the variance and bias are of order O(1/N ) [START_REF] Del | Mean Field Simulation for Monte Carlo Integration[END_REF].

The proposed algorithm is described in Algorithm 2. Each training loop iteration of CRAFT looks like a iteration of the test time sampler in Algorithm 3 with gradient estimation and flow parameter updates added. Note that the parameter update is applied after the flow transport step for that temperature so the first effect it has will be on the particles that come through in the next step. This allows us to view CRAFT as an optimization of the objective (5) and will allow stability analysis that we will detail later in the text. It also means that the normalizing constant estimators produced after each step are unbiased. Since all the constants produced in this way depend on the optimization sequence they are not independent. At test time when the parameters are fixed the independence is restored.

Relationship to AFT. We now describe the relationship between CRAFT and the AFT method of [START_REF] Arbel | Annealed flow transport Monte Carlo[END_REF]. Relating them is made more challenging by the fact there are two variants of the AFT algorithm described in [START_REF] Arbel | Annealed flow transport Monte Carlo[END_REF]. For clarity we include both in our supplement as Algorithms 5 and 7. We call Algorithm 5 simple AFTit is easier to describe and analyse but is not actually used in practice by [START_REF] Arbel | Annealed flow transport Monte Carlo[END_REF] for reasons we will describe.

Instead they use what we call practical AFT (Algorithm 7) a more complicated version that performs better in practice. Relative to CRAFT, both variants of AFT may be viewed as a one pass greedy optimization of the objective (5). The closest point of algorithmic comparison to CRAFT is simple AFT-we show the line difference to CRAFT in Algorithm 6. At each temperature step simple AFT optimizes the particle approximated loss to its numerical minimum, before then updating the particles using the new flow.

We now describe the sample replenishment problem in AFT. The learning of the flows is hampered by the fact that the sample complexity of estimation can be higher than the finite number of available particles. In simple AFT this manifests as over-fitting of the flows to the available particles and large biases in expectations. Regenerating new samples at each temperature separately would require computation quadratic in the number of temperatures, which is too slow. To mitigate the over-fitting in a manor consistent with the one pass paradigm of the paper, the authors added validation samples which were used for early stopping. Test samples were also added to remove finite sample residual bias. This significantly increased the complexity of the AFT method, changing simple AFT (Algorithm 5) to practical AFT (Algorithm 7). This gave substantial empirical improvements but the fundamental limitation of the sample replenishment problem remains, arising from the desire to use one temperature pass. Even when the particles are a perfect approximation to the desired distribution at each temperature, if the sample complexity of estimating the training gradients of the flows is high relative to the number of particles then the practical AFT method still struggles to train. By having multiple independent passes instead of one, CRAFT is able to fully address the sample replenishment problem in a different way. Since in CRAFT the parameter update is applied after flow transport at each temperature the analysis of the algorithm is simpler than that required for AFT, where sophisticated arguments must be invoked to ensure consistency of the particles. Later in Section 3.2 we will see that we can analyse CRAFT from the perspective of unbiased gradients. This interpretation does not apply to either variant of AFT.

From a practical perspective CRAFT performs better than AFT (Section 4.1). This is significant because [START_REF] Arbel | Annealed flow transport Monte Carlo[END_REF] already showed they could train samplers that performed well relative to strong baselines like SMC. CRAFT is easier to tune and more robust to hyperparameter misspecification than AFT. Conceptually it is much simpler, to the extent that it is easier to describe CRAFT as a standalone algorithm than in terms of the AFT algorithm on which it builds. CRAFT is closer to a traditional machine learning training paradigm with gradients applied to a single objective. This makes it easier to adopt and to scale through parallelism.

Recent work from [START_REF] Zimmermann | Nested variational inference[END_REF] investigates a framework that reduces to the AFT objective (5) for deterministic forward and backward transitions. Inspired by [START_REF] Arbel | Annealed flow transport Monte Carlo[END_REF] they investigate some toy examples with normalizing flows. There are good reasons to concentrate on the normalizing flow case as we have done. The reverse distribution is analytic and optimal for a flow but can be difficult to approximate otherwise [START_REF] Achille Thin | Monte Carlo variational auto-encoders[END_REF]. Further [START_REF] Zimmermann | Nested variational inference[END_REF] do not use Markov transition kernels in practice, which we found to be essential for good performance.

Using CRAFT within Particle MCMC

Whilst the SMC-NF estimator of the normalizing constant is unbiased, in general estimates of the expectations w.r.t. the target are asymptotically consistent but exhibit a O(1/N ) bias. In situations requiring high accuracy expectations it is desirable to have a method returning asymptotically unbiased estimates. However, increasing N for the bias to be negligible will often not be feasible for challenging applications because of memory requirements. MCMC algorithms have the advantage of providing consistent estimate of these expectations as compute time increases without having to store an increasing number of samples. Particle MCMC methods [START_REF] Andrieu | Particle Markov chain Monte Carlo methods[END_REF] provide a way to bring such benefits to SMC samplers. In particular, the so-called Particle Independent Metropolis-Hastings is an MCMC sampler that uses an SMC sampler with N particles and provides consistent estimates of expectations for any N as the number of iterations increases; see Appendix D for details. Since a single pass of CRAFT with fixed parameters can be thought of as an SMC sampler with additional deterministic transformations it can be used here. Computational trade-offs for Particle MCMC are analyzed by [START_REF] Michael K Pitt | On some properties of Markov chain Monte Carlo simulation methods based on the particle filter[END_REF].

Analysis of the training objective

Reformulating the CRAFT objective

The CRAFT objective (5) can be rewritten as a single KL divergence between certain product distributions:

H = KL K k=1 T # k π k-1 || K k=1 π k . ( 9
)
We have adopted the notation that the product symbol can be used to denote multi-variable product measures.

The proof above effectively applies the multi-variable generalization of the result for the KL-divergence between two product mea-

sures, KL[U ×V ||F ×G]=KL[U ||F ]+KL[V ||G].
Despite the non-standard form of the KL divergence, it is also possible to obtain a bound on the normalizing constant

log Z K ≥ - K k=1 E X∼π k-1 [D k ] (10)
This bound can be obtained starting from H ≥ 0, expanding each term using (6a) and noting that there is a telescoping cancellation between the intermediate normalizing constants. The bound is distinct from the one obtained if we take the logarithm of the unbiased estimate used in SMC samplers. In practice the bound is estimated using the particle estimate for each temperature and is effectively the training objective of CRAFT.

Re-thinking the KL criterion at the particle level

We discuss here a re-interpretation of the CRAFT training procedure in terms of KL-divergences between particle approximations and targets. An advantage of this interpretation is that it shows that we compute an unbiased gradient of a modified objective that is well-motivated even when the number of particles is moderate.

CRAFT training requires minimizing the objective (5).

As we do not have access to π k-1 , CRAFT approximates the intractable gradient of this objective using equation ( 8). We will compactly refer to the particle approximation as a weighted sum of delta masses

π N k-1 (dx)= N i=1 W i k-1 δ X i k-1 (dx).
While the use of such an approach could be of concern when N is moderate or at initialization, the following proposition shows that such concerns are ill-founded and that this approximate SMC based gradient is also the unbiased gradient of a welldefined and intuitive objective.

Proposition 2. Let πN k-1 (dx) = E[π N k-1 (dx)]
E X∼π N k-1 [∇ θ k D k (X)]. (11)
The l.h.s of this KL is the pushforward of the "average" SMC approximation to π k-1 through the flow T k .

To help understand this result and how it is derived, consider a simpler scenario where the random approximation π N k-1 of π k-1 has been obtained by a batch parallel importance sampling method. That is we sample

X i k-1 i.i.d. ∼ q k-1 for i=1, . . . , N , compute w(X i k-1 )=π k-1 (X i k-1 )/q k-1 (X i k-1 ). We then define W i k-1 ∝w(X i k-1 ) such that N i=1 W i k-1 =1 and finally re- turn π N k-1 (dx).
If we average over the random particle locations X 1:N k-1 , we obtain

πN k-1 (dx) = E X i k-1 i.i.d. ∼ q k-1 [π N k-1 (dx)]. (12)
Contrary to π N k-1 (dx) which is a discrete measure, the distribution πN k-1 (dx) admits a density and its Kullback-Leibler divergence with µ k = (T -1 k ) # π k is well-defined. It then follows directly from ( 12), diffeomorphism invariance, and iterated expectations that

KL[T # k πN k-1 ||π k ] = E X∼π N k-1 log πN k-1 (X) µ k (X) (13) = E X i k-1 i.i.d. ∼ q k-1 E X∼π N k-1 log πk-1 (X) µ k (X) .
A direct consequence of this identity is that an unbiased gradient of KL[T # k πN k-1 ||π k ] is indeed given by (11). To prove Proposition 2, we extend this argument to the scenario where πN k-1 (dx) has been obtained by using an SMC sampler where the particles are not independent because of resampling; see Appendix C. Taking a step back and returning to the original objective we note that once we start considering the marginal particle distribution πN k-1 instead of the target distribution π k-1 the flow parameters θ effect the KL divergences that follow them, since later particle distributions depend on them. This additional effect is not included in the CRAFT updates. Each flow takes whatever particle distribution it receives from the steps and is trained only to reduce the KL of the push-forward to the next target distribution π k , without reference to what follows.

Experiments

In this section we empirically investigate the performance of CRAFT. First, in Section 4.1 we give a case study demonstrating the empirical benefit of CRAFT relative to AFT, then in Section 4.2 we show that CRAFT outperforms Stochastic Normalizing flows in two challenging examples. We then show a compelling example use case for CRAFT as a learnt proposal for a particle MCMC sampler applied to lattice field theory. Code for the algorithms and examples can be found at https://github.com/deepmind/annealed_ flow_transport. Further experiments and details are included in the Appendix.

Empirically comparing CRAFT and AFT

Our first point of comparison for CRAFT is with AFT [START_REF] Arbel | Annealed flow transport Monte Carlo[END_REF]. We use what we call the practical AFT method in Section 2.3. To illustrate the benefit of solving the AFT sample replenishment problem we show the effect of varying the batch size in both algorithms. Both methods are given the same total budget of particles at train and test time. To highlight that the sample replenishment problem cannot be mitigated simply by using more MCMC we gave AFT 100 times more MCMC updates than CRAFT at train and test time. We use the 1024 dimensional log Gaussian Cox process (LGCP) example which is the most challenging from [START_REF] Arbel | Annealed flow transport Monte Carlo[END_REF]. As in that paper, we used a diagonal affine flow.

Figure 1 shows the results. We see that CRAFT outperforms AFT for all number of particles considered. This can be understood in terms of our better solution to the sample replenishment problem as described in Section 2.3. Since many of the most challenging modern sampling problems also have high memory usage this is also of considerable practical significance.

Comparing CRAFT and SNFs

We now compare CRAFT training with that of the corresponding SNFs. For the Markov transition kernel both methods use full Metropolis corrected HMC. We use the same normalizing flow family for both methods. As such this constitutes a direct comparison of the two approaches. We ran timed training experiments for both methods computing the corresponding unbiased estimates for the normalizing constant as it progressed. Note that at test time these estimates could be obtained faster by running the learnt sampler. We consider two different target densities following [START_REF] Arbel | Annealed flow transport Monte Carlo[END_REF] namely the 30 dimensional variational autoencoder latent space and the LGCP example. For the VAE example we used an affine inverse autoregressive flow [START_REF] Durk P Kingma | Improved variational inference with inverse autoregressive flow[END_REF]. For the LGCP we again used a diagonal affine flow.

Figure 2 shows the results. In both cases we find that CRAFT converges to a better value than SNFs. We attribute the worse final value for SNFs to the issues with the training discussed in Section 2.2.

CRAFT based Particle MCMC for lattice φ 4 theory

A classic proving ground for sampling algorithms is in the area of lattice field theory. Indeed, it was this application that motivated Hamiltonian/Hybrid Monte Carlo [START_REF] Duane | Hybrid Monte Carlo[END_REF]. State of the art quantum chromodynamics (QCD) calculations (Borsanyi et al., 2021) use lattice sampling run on multiple supercomputers. There has been a recent surge of interest in using normalizing flows to improve sampling in such models starting with [START_REF] Albergo | Flow-based generative models for Markov chain Monte Carlo in lattice field theory[END_REF], so it is natural to investigate the performance of CRAFT in this area.

Of the available lattice field theories we use the Euclidean φ 4 theory in two dimensions following [START_REF] Albergo | Flow-based generative models for Markov chain Monte Carlo in lattice field theory[END_REF]2021). We start from the continuum action:

S cont [φ] = ||∇φ(x)|| 2 2 + m 2 φ(x) 2 + λφ(x) 4 dx,
which we discretize on a 14 × 14 lattice, using lattice units, and periodic boundary conditions to obtain:

S latt (φ) = x {φ(x)ζ(x) + m 2 φ(x) 2 + λφ(x) 4 }, (14a) ζ(x) := µ [2φ(x) -φ(x + êµ ) -φ(x -êµ )] . (14b)
Here the sum over x runs over all lattice sites. The summation µ ∈ {x, y} runs over the dimensions of the lattice and êµ defines a single lattice step the in direction µ. λ ≥ 0 defines the coupling parameter and m 2 is the mass squared which counter-intuitively may be positive or negative. A probability distribution over the values of the field at the lattice sites is then defined using p(φ) = 1 Z exp {-S latt (φ)}. For positive λ and large φ the quartic term will dominate giving highly non-Gaussian tails. For negative m 2 , which will be our focus, the marginal distributions are bimodal.

As is typical of physical sciences applications, symmetry is crucial in this example. The target is invariant under any translation on the lattice. In Appendix H we show that the corresponding requirement on flows between successive temperatures is translation equivariance. We incorporate this in to the flow models under consideration using convolutional affine coupling layers [START_REF] Dinh | Density estimation using real NVP[END_REF] with periodic boundary conditions. Since high accuracy is required we are interested in methods that are asymptotically unbiased for general expectations. We take physically relevant expectations from (Albergo et al., 2019) namely the 'two point susceptibility' in the main text and the 'Ising energy density' in the supplement -both have similar results. We use Particle MCMC with three different independent proposal types. These are SMC, VI with normalizing flows (and no annealing) and CRAFT. We also compare against HMC that directly targets the target density. The results are shown in Figure 3. We see that HMC struggles with this target density because it is unable to mix between modes. Despite producing samples quickly, the VI proposal fails to cover all of the modes because of the use of reverse KL divergence. The SMC proposal does much better. As an annealing method it is less susceptible to multi-modality. Finally CRAFT converges the fastest, improving on SMC without flows and avoiding the mode seeking behaviour of pure VI. Note that by scaling arguments if we accelerate an MCMC method by a factor of R we expect that errors in expectations will be reduced by a factor of √ R on such error plots.

Conclusion

In this paper we have proposed and analyzed the CRAFT algorithm. At test time it is an SMC sampler with added flow transport steps. The training algorithm uses particle estimation to optimize a training objective that explicitly promotes transport between annealing temperatures. We have given a thorough analysis of the training including an insightful particle interpretation. We have shown that CRAFT compares well empirically and conceptually with existing methods to estimate flows in this context.

In terms of where such efforts fit into the bigger picture, there are many applications where high accuracy expecta-tions are required and algorithms like Particle MCMC become necessary due to the challenging nature of the target.

As such we believe the speed up we have achieved in the quantum field theory application is exciting both for that field and others.

Continual Repeated Annealed Flow Transport Monte Carlo: Supplementary Material

A. Extended proposal and target distributions

Assume the transport maps T l are fixed. We detail the extended proposal and target distributions used by Stochastic Normalizing Flows [START_REF] Wu | Stochastic normalizing flows[END_REF]. These were introduced early on in [START_REF] Vaikuntanathan | Escorted free energy simulations[END_REF]. The same target and proposal are also used by CRAFT and Annealed Flow Transport (Arbel et al., 2021) and the Gibbs flow algorithm [START_REF] Heng | Gibbs flow for approximate transport with applications to Bayesian computation[END_REF] when resampling is omitted. We include here these distributions for sake of completeness.

We sample

X 0 ∼ π 0 (•) at k = 0 then use Y k = T k (X k-1 ) followed by X k ∼ K k (Y k , •) at time k ≥ 1.
To simplify presentation, we do not use measure-theoretic notation. Hence, using the notation

M trans l (x, x ) = δ T l (x) (x ) and M mut l (x, x ) = K l (x,
x ), we obtain the following proposal at time k after the transport step

ηk (x 0:k-1 , y 1:k ) = π 0 (x 0 ) k-1 l=1 M trans l (x l-1 , y l )M mut l (y l , x l ) M trans k (x k-1 , y k ), (15)
and the target is

πk (x 0:k-1 , y 1:k ) = π k (y k )L trans k-1 (y k , x k-1 ) k-1 l=1 L mut l-1 (x l , y l )L trans l-1 (y l , x l-1 ) , (16)
where

L trans l-1 (x, x ) = δ T -1 l (x) (x )
and L mul l-1 (x, x ) = π l (x )M mul l (x , x)/π l (x). After the mutation step at time k, the proposal is

ηk (x 0:k , y 1:k ) = ηk (x 0:k-1 , y 1:k )M mut k (y k , x k ) = π 0 (x 0 ) k l=1 M trans l (x l-1 , y l )M mut l (y l , x l ) (17)
and the target is

πk (x 0:k , y 1:k ) = π k (x k ) k-1 l=0 L mut l (x l+1 , y l+1 )L trans l (y l+1 , x l ) . ( 18
)
Hence the incremental weight after a transport term at time k is of the form

πk (x 0:k-1 , y 1:k ) ηk (x 0:k-1 , y 1:k ) = πk-1 (x 0:k-1 , y 1:k-1 ) ηk-1 (x 0:k-1 , y 1:k-1 ) π k (y k )L trans k-1 (y k , x k-1 ) π k-1 (y k-1 )M trans k (x k-1 , y k ) incremental weight= Z k-1 Z k G k,T k (x k-1 ) , ( 19
)
while after the mutation step it is of the form

πk (x 0:k , y 1:k ) ηk (x 0:k , y 1:k ) = πk (x 0:k-1 , y 1:k ) ηk (x 0:k-1 , y 1:k ) π k (x k )L mul k-1 (x k , y k ) π k (y k )M mul k (y k , x k ) incremental weight=1 . ( 20
)
As shown in Appendix B.2 of [START_REF] Arbel | Annealed flow transport Monte Carlo[END_REF], we can also integrate out the random variables Y 1:k in these expressions as we can collapse the transport step and mutation step into one single Markov kernel

M col l (x l-1 , x l ) = M trans l (x l-1 , y l )M mut l (y l , x l )dy l = δ T l (x l-1 ) (y l )K l (y l , x l )dy l = K l (T l (x l-1 ), x l ). ( 21
)
Similarly we can collapse the backward kernels that have been used to defined the extended target distributions πk

L col l-1 (x l , x l-1 ) = L mut l-1 (x l , y l )L trans l-1 (y l , x l-1 )dy l = π l (T l (x l-1 ))|∇T l (x l-1 )|K l (T l (x l-1 ), x l ) π l (x l ) , ( 22
)
so the collapsed proposal and target can be written as

ηk (x 0:k ) = π 0 (x 0 ) k l=1 M col l (x l-1 , x l ), ( 23
) πk (x 0:k ) = π k (x k ) k-1 l=0 L col l (x l+1 , x l ). ( 24
)
For k = K, we write ηK = η and πK = π.

It follows that

Z K η(x 0:K ) π(x 0:K ) := w K (x 0:K-1 ) = γ K (x K ) γ 0 (x 0 ) K l=1 L col k-1 (x k , x k-1 ) M col k (x k-1 , x k ) (25) = γ K (x K )γ K-1 (x K-1 ) • • • γ 1 (x 1 ) γ 0 (x 0 )γ 1 (x 1 ) • • • γ K-1 (x K-1 ) K k=1 L col k-1 (x k , x k-1 ) M col k (x k-1 , x k ) = K k=1 γ k (x k )L col k-1 (x k , x k-1 ) γ k-1 (x k-1 )M col k (x k-1 , x k ) = K l=1 G k (x k-1 ), (26)
where, from ( 21) and ( 22), we obtain the so-called incremental weights

G k (x k-1 ) = γ k (x k )L col k-1 (x k , x k-1 ) γ k-1 (x k-1 )M col k (x k-1 , x k ) = γ k (T k (x k-1 )) γ k-1 (x k-1 ) |∇T k (x k-1 )| (27)
as given in (1b). In [START_REF] Wu | Stochastic normalizing flows[END_REF], the weight w K (x 0:K-1 ) is instead obtained by first computing recursively the product appearing on the RHS of ( 25), thus they consider formally1 in equation ( 12) of their paper a log-weight update at iteration k of the form

log L col k-1 (x k , x k-1 ) M col k (x k-1 , x k ) = log(γ k (T k (x k-1 )|∇T l (x l-1 )|) -log γ k (x k ) = log γ k (x SNF k ) -log γ k (x SNF k ) (28)
where xSNF k is the value of the sample after the MCMC at temperature k has been applied and x SNF k is the value of the sample before it as described in (4).

So we can rewrite the ELBO as

L = E ηK (x 0:K ) log w K (x 0:K-1 ) = E π0(x0)g(u 1:K )β(a 1:K |x0,u 1:K ) log w K (x 0:K-1 ) , ( 38
)
where we recall that the states x 0:K are deterministic functions of x 0 , u 1:K , a 1:K given by ( 36). Now when taking the gradient of the ELBO w.r.t. parameters of the NFs and/or stochastic kernels then we have

∇L = E π0(x0)g(u 1:K )β(a 1:K |x0,u 1:K ) ∇ log w K (x 0:K-1 ) + ∇ log β(a 1:K |x 0 , u 1:K ) • log w K (x 0:K-1 ) . ( 39
)
Now set u 0 =x 0 , u=u 0:K , u=a 1:K . Hence, the random variables (A, U ) admit a joint distribution p of the form p(a, u):=π 0 (u 0 )g(u 1:K )β(a 1:L |x 0 , u 1:K ). Moreover, define Φ the re-parametrization that maps (A, U ) to the samples

X 0:K , i.e. [Φ(A, U )] l = Φ l (U 0:l , A 1:l ) = X l .
Hence, we can express the gradient in Equation ( 39) as follows:

E p [∇ log(w K • Φ) + log(w K • Φ) ∇ log p] . (40)
The first term on the r.h.s. corresponds to the reparametrization trick while the second term is a score term, which in general can have high variance.

The gradient estimator of [START_REF] Wu | Stochastic normalizing flows[END_REF] effectively uses an unbiased estimator of the first reparameterization term but neglects the second term altogether.

In the particular case where the flows transport exactly between each temperature, i.e. T # l π l-1 =π l , the importance weight w K (x 0:K-1 ) has the constant value Z K . Hence, the second term in the expectation is a constant multiplied by a score function and consequently it vanishes:

E p [log(w K • g)∇ log p] ∝ E p [∇ log p] = ∇p = 0. (41)
Moreover, since the ELBO is maximized when using optimal flows, the gradient ∇L must vanish, which directly implies that the first term must also vanish. Therefore, despite the missing term in the gradient, the SNF learning rule has the correct fixed point in expectation.

C. Particle Interpretation of CRAFT

CRAFT is attempting to minimize w.r.t.

θ the KL divergence KL(T # k (θ)π k ||π k+1 ) = KL(π k ||(T -1 k (θ)) # π k+1 )
We have shown in Section 3.2 that we could re-interpret the criterion minimized by CRAFT at time k as the KL divergence between the expectation of the random measure approximating π k and the pullback of π k+1 by T k (θ). This was established in the case where the random measure is obtained using importance sampling. We show here that it is also valid for the case considered by CRAFT where it arises from an SMC-NF algorithm.

For sake of simplicity, we assume resampling is performed at any time step. Then in this case, the particles generated by the SMC-NF algorithm are distributed according to

qk (dx 1:N 0:k , dy 1:N 1:k ) = N i=1 π 0 (x i 0 ) k l=1 N i=1 δ T l (x i l-1 ) (dy i l ) N i=1 W a i l-1 l K l (y a i l-1 l , dx i l ) , (42)
where a i l-1 ∈ {1, ..., N } is the ancestral index of particle x l i . Using random particles arising from the SMC-NF algorithm, the distribution π k is then approximated by the random empirical measure

π N k (dx) = 1 N N i=1 δ X i k (dx) (43)
and its expectation w.r.t. ( 42) is denoted

πN k (dx) = E[π N k (dx)]. (44)
In this case, we have

KL(π N k ||T -1 k (θ) # π k+1 ) = E X∼π N k log πN k (X) (T -1 k (θ)) # π k+1 (X) = E (X 1:N 0:k ,Y 1:N 1:k )∼q k ,X∼π N k log πN k (X) (T -1 k (θ)) # π k+1 (X) = E (X 1:N 0:k ,Y 1:N 1:k )∼q k E X∼π N k log πk (X) (T -1 k (θ)) # π k+1 (X) . ( 45
)
So an unbiased gradient estimate w.r.t.

θ of -KL(π k ||T -1 k (θ) # π k+1 ) is given by -E X∼π N k ∇ θ log T -1 k (θ) # π k+1 , (46)
which is the gradient used by CRAFT to update the parameter of T k (θ).

D. Particle MCMC

We describe here the particle independent Metropolis-Hastings algorithm (PIMH) introduced in [START_REF] Andrieu | Particle Markov chain Monte Carlo methods[END_REF]. This is an MCMC algorithm using as an independent proposal an importance sampling or SMC-type algorithm such as an SMC combined with NFs learned by CRAFT (see Algorithm 3). It generates a Markov chain

((π N K (j), Z N K (j)) j≥1 where π N K (j)(dx) = N i=1 W i K (j)δ X i K (j) (dx) such that for any test function f 1 J J j=1 E X∼π N K (j) [f (X)] → E X∼π K [f (X)], where E X∼π N K (j) [f (X)] := N i=1 W i K (j)f (X i K (j)), (47)
as the number J of MCMC iterations goes to infinity. This is described in detail in Algorithm 4.

Algorithm 4 Particle Independent Metropolis-Hastings step P ((π N K (j),

Z N K (j)), (•, •)) 1: Input: Current approximations (π N K (j), Z N K (j)) to (π K , Z K ). 2: Output: New approximations (π N K (j + 1), Z N K (j + 1)) to (π K , Z K ). 3: Run an SMC algorithm to obtain (π ,N K , Z ,N K ) approximating (π K , Z K ) (such as Algorithm 3). 4: Set (π N K (j + 1), Z N K (j + 1)) = (π ,N K , Z ,N K ) with probability α(Z N K (j), Z ,N K ) := min 1, Z ,N K Z N K (j) . 5: Otherwise set (π N K (j + 1), Z N K (j + 1)) = (π N K (j), Z N K (j).
For fixed computational efforts, one could either run many MCMC iterations with N small or few MCMC iterations with N large. Under regularity conditions, it was shown in [START_REF] Michael K Pitt | On some properties of Markov chain Monte Carlo simulation methods based on the particle filter[END_REF] that N should be selected such that the variance of the estimate of the log-normalizing constant is close to 1 in order to minimize the asymptotic variance of the estimator (47).

E. Learning model parameters using CRAFT

We do not consider in the paper examples where π K depends on some parameters φ that one wants to learn. We explain here how this could be addressed.

Consider the following target

π φ K (x) = γ φ K (x) Z φ K . (48)
To estimate φ, we propose to learn φ by maximizing the log-normalizing constant/evidence arg min

φ log Z φ K . (49)

H. Incorporating symmetries

In applications to high dimensional target densities the target density often exhibits symmetries. Incorporating such symmetries in normalizing flows can give significant improvements in the performance and efficiency of normalizing flows.

Consider a case where we have a group G with elements g and corresponding group representation members D g . In this case we have a symmetric target density π T which is symmetric under the target density.

π T (D g x) = π T (x) ∀g ∈ G (51)
We also take a base density π B that is also symmetric under the group. We choose a normalizing flow that is equivariant under the action of the group so that:

T (D g x) = D g T (x) ∀g ∈ G (52)
When we push forward the base density we end up with a distribution T # π B that is invariant under the action of the group.

T # π B (D g x) = T # π B (x) ∀g ∈ G (53)
In standard normalizing flow modelling the symmetric pushforward distribution is then tuned to match the symmetric target distribution using reverse KL-divergence KL[T # π B ||π T ].

It is important to extend such a treatment of symmetries to the case of CRAFT. To this end we analyze the symmetry properties of a sequence of densities with a geometric annealing schedule. We consider an initial normalized density π 0 and unnormalized final density γ 1 . The unnormalized density of the geometric annealing schedule takes the form:

log γ β (x) = log(1 -β) log π 0 (x) + β log γ 1 (x) (54)
where the schedule time β ∈ [0, 1]. Clearly if both the initial density and final density are symmetric under the action of the group then so is γ β (x) and its normalized version π β .

In CRAFT we use a sum of KL divergences between the pushforward of an initial annealed density and the next target density KL[T #π β ||π β ]. In terms of symmetry there is a direct analogy to the case of standard variational inference with normalizing flows with which we started the discussion. In particular π β now takes the role of the symmetric base density and π β takes the role of the target. The symmetry requirement on the normalizing flow T is thus the same. It needs to be equivariant under the action of the group so that the push forward T #π β is symmetric in a manner that matches its target.

I. Additional experiment details

All experiments used a geometric (log-linear) annealing schedule. The initial distribution was always a standard multivariate normal. All experiments used HMC as the Markov kernel, which was tuned to get a reasonable acceptance rate based on preliminary runs of SMC. Normalizing flows were always initialized to the identity flow. Wherever a stochastic gradient optimizer was required we used the Adam optimizer [START_REF] Diederik | Adam: A method for stochastic optimization[END_REF].

In terms of software dependencies for our code we used Python, JAX [START_REF] Bradbury | JAX: composable transformations of Python+NumPy programs[END_REF], Optax [START_REF] Hessel | Optax: composable gradient transformation and optimisation[END_REF], Haiku [START_REF] Hennigan | Haiku: Sonnet for JAX[END_REF], and the TensorFlow probability JAX 'substrate' [START_REF] Dillon | Tensor-Flow Distributions[END_REF].

All experiments were carried out using a single Nvidia v100 GPU. We used sufficient CPUs and CPU RAM such that these were not the bottleneck.

I.1. Comparing CRAFT and AFT batch size performance

Here we give more details of the experiment described in Section 4.1. We used the original software for AFT which is available at https://github.com/deepmind/annealed_flow_transport. We used practical AFT (Algorithm 7), which is the version using early stopping on a set of validation particles and a 'hold out' test set.

We used the same MCMC method for both approaches, which was pre-tuned full Metropolis adjusted Hamiltonian Monte Carlo (HMC) with 10 leapfrog steps per iteration. For CRAFT we used 10 HMC steps per temperature and for AFT we used 1000 HMC steps per temperature meaning that AFT had 100 times more HMC steps per temperature. The HMC step sizes were linearly interpolated with 0 corresponding to the initial distribution and 1 corresponding to the final distribution, the interpolation points were [0, 0.25, 0.5.1] and the corresponding step sizes were [0.3, 0.3, 0.2, 0.2].

At test time AFT and CRAFT used the same number of particles for each experiment. To make it fair at train time, the total CRAFT particle budget was divided in to two halves for AFT, one half was used for the training particles and the other half was used for the validation particles. For each experiment the total number of train particles for either method was the same as the number used at test time and this is the number shown in Figure 1.

In terms of optimization both methods were well converged. AFT used an optimization step size of 1e-2, and 500 optimization iterations per temperature. CRAFT used 200 total optimization iterations, and a step size of 5e-2 which was reduced down to 1e-2 after 100 iterations.

I.2. Comparing CRAFT and Stochastic Normalizing Flows

To the best of knowledge, the HMC based Stochastic Normalizing Flow described in [START_REF] Wu | Stochastic normalizing flows[END_REF] was not implemented or experimented with in that work. We implemented it and found in preliminary experiments that it outperformed the random walk Metropolis Hastings used in the original work. It was also the most commensurate with our own HMC usage. Therefore we used the HMC SNF for the comparisons in this work.

In the spirit of the original SNF paper, we performed preliminary experiments with a SNF ELBO that incorporated resampling in the forward computation. For the reverse computation we neglected the resampling contribution just as is done for the score term arising from the discrete Metropolis-Hastings correction. In other words we proceeded as if the forward computation was compatible with the reparameterization gradient although this does not in fact correspond to an unbiased estimate of the gradient. These preliminary experiments indicated that this variant of SNF became unstable and challenging to train. The fact that CRAFT can cope with resampling is a benefit of the approach, and is expected from the form of the CRAFT training objective.

We observed that variational learning of MCMC step sizes is unstable in SNF, so we followed [START_REF] Wu | Stochastic normalizing flows[END_REF] who use a fixed step size during variational training for the majority of their experiments. Instead of a learnt step size, we tuned the step sizes separately in advance so that the corresponding MCMC would have a reasonable acceptance rate just as we do for CRAFT. Consequently we used the same Markov kernels for both CRAFT and SNFs.

For the SNF software we implemented them in JAX using similar modules to that of the CRAFT implementation. We verified this JAX implementation against the original publicly available SNF implementation at https://github.com/ noegroup/stochastic_normalizing_flows. A benefit of using the same basic modules as the CRAFT implementation is that the timed comparison is much less subject to unwanted differences arising from the different software libraries used. We confirmed using a Python execution profiler that both the CRAFT and SNF code was spending time on the expected core algorithmic computations and not on unexpected code paths. We observed some variability in the CRAFT and SNF step times.

The JAX SNF implementation exploits the connection between SNFs and Annealed Importance Sampling with normalizing flows. This equivalence can be easily seen from the form of the overall unnormalized weights in each case (Appendix A).

The gradient dynamics are unchanged by using this representation of the forward computation.

Relative to CRAFT we found that SNFs used significantly more memory in our experiments. This is due to the backward pass of the SNFs where gradients are passed through the whole forward computation whereas the gradients are local in CRAFT. To the benefit of SNFs, we reduced the batch size of both CRAFT and SNFs in the relevant experiments so that the SNF training would still fit on GPU and that the batch sizes would be comparable for the two algorithms.

The normalizing flow used for the VAE experiment was an Affine Inverse Autoregressive Flow. Similiar to [START_REF] Arbel | Annealed flow transport Monte Carlo[END_REF] we use a masked network [START_REF] Mathieu Germain | International Conference on Machine Learning[END_REF] to achieve the autoregressive dependency structure. The masked network had two hidden layers and unmasked would have 150 hidden units. The final weights and biases were initialized to give an identity flow. We used a leaky ReLU nonlinearity. Since the flows have high capacity in this case we only needed to use 3 temperatures in both cases. Since the cost of these inverse autoregressive flows is quadratic in the dimension they are prohibitively expensive for systems of higher dimensionality. We used 2 HMC steps per transition with 10 leapfrog steps in both cases. Both CRAFT and SNF had a batch size of 100.

For the LGCP SNF/CRAFT comparison we used 10 transitions in both cases. We used 1 HMC step per transition with 10 leapfrog steps in both cases. Both CRAFT and SNF had a batch size of 2000.

I.3. CRAFT based Particle MCMC for lattice φ 4 theory

Those readers who are less familiar with the background physics can gain intuition from the function form of equation ( 14). The terms over neighbouring lattice sites promote spatial correlation. The other terms all involve only a single lattice site.

The expectations we use to evaluate the algorithms are physically motivated observables described by [START_REF] Albergo | Flow-based generative models for Markov chain Monte Carlo in lattice field theory[END_REF]. In the main text we use the two point susceptibility. In Figure 5 we show that similar results show for the Ising energy density. The chains that were averaged to produce the expectations in 3 and 5 are shown in 6 and 7 respectively. The errors shown in the average value plots are computed by subtracting off the gold standard value.

We note that de Haan et al. ( 2021) incorporate additional symmetries of the target on top of the lattice translation symmetry we incorporate. This is at the expense of introducing numerical integration into the variational inference with normalizing flows. The convolutional affine coupling layers we use do not require such numerical integration, and are known to significantly outperform fully connected affine coupling layers in this context [START_REF] Michael S Albergo | Introduction to normalizing flows for lattice field theory[END_REF].

The parameters we consider for the theory are based on those from [START_REF] Albergo | Flow-based generative models for Markov chain Monte Carlo in lattice field theory[END_REF]. In particular we take the hardest parameter set E5 and make it more difficult by reducing m 2 from -4 to -4.75. As noted by [START_REF] Albergo | Flow-based generative models for Markov chain Monte Carlo in lattice field theory[END_REF], the reason this makes the problem more challenging is it increases the gap between the positive and negative modes. The resulting parameters we use are therefore the following: Lattice size 14 × 14; λ = 5.1 and m 2 = -4.75.

In these φ 4 experiments, upon preliminary investigation of computation time using a Python execution profiler, we found that time was sometimes being spent in parts of the code that might not be expected from algorithmic considerations alone.

In particular for fast samplers like VI, we found that storing results and computing expectations were an (unoptimized) bottleneck. To make timed comparisons fairer, we therefore performed separate timing estimation on code that ran only the core sampling code and not on computing expectations and storing results. These adjusted timing estimates are what is reported in our results.

The gold standard values of the error for each reported expectation are computed using SMC. We verified in cases of simple known expectations that SMC and CRAFT gave the correct answers whereas HMC and VI gave large errors, just as for the physically motivated observables reported.

The normalizing flow used a convolutional affine coupling layer [START_REF] Dinh | Density estimation using real NVP[END_REF] with periodic boundary conditions. We used alternating checkerboard masking patterns. The convolution kernels were of size 3 × 3. The conditioner of each flow had one hidden convolutional layer with 10 channels. The final weights and the biases of the convolution were initialized to give an identity flow. We used a ReLU nonlinearity.

SMC used 90 transitions. Since the masking pattern required two coupling layers to change each site we used two coupling layers per transition in CRAFT. CRAFT used 10 transitions. VI, which has no annealing, used 20 coupling layers. We confirmed there was no benefit to increasing the number of VI coupling layers. We used 10 HMC steps per annealing transition with 10 leapfrog steps for SMC and CRAFT. Direct HMC used 10 leapfrog iterations per step.

In this section we describe the two different variants of AFT from [START_REF] Arbel | Annealed flow transport Monte Carlo[END_REF] for clarity and completeness. We also describe the line difference between CRAFT and simple AFT. Note that the full source code for each method is also available at https://github.com/deepmind/annealed_flow_transport. Based on equation ( 7) we define:

L N k (T ) := KL[T # k π k-1 ||π k ] -log Z k Z k-1 ≈ i W i k-1 D k (X i k-1 ) (55)
which from equation ( 8) has approximate gradient given by i

W i k-1 ∂D k (X i k-1 ) ∂θ k . ( 56
)
Algorithm 5 Simple AFT: from Arbel, Matthews and Doucet 2021. This corresponds to Algorithm 1 of the prior paper.

1: Input: number of particles N , unnormalized annealed targets {γ k } K k=0 such that γ 0 = π 0 and γ K = γ, resampling threshold A ∈ [1/N, 1). 2: Output: Approximations π N K and Z N K of π and Z. 3: Sample X i 0 ∼ π 0 and set W i 0 = 1 N and Z N 0 = 1. 4: for k = 1, . . . , K do 5:

Solve T k ← argmin T ∈T L N k (T ) using e.g. SGD.

6:

π N k , Z N k ← SMC-NF-step π N k-1 , Z N k-1
, T k 7: end for Algorithm 6 CRAFT to simple AFT line difference: additions for CRAFT shown in green and removals from simple AFT shown in red 1: Input: number of particles N , unnormalized annealed targets {γ k } K k=0 such that γ 0 = π 0 and γ K = γ, resampling threshold A ∈ [1/N, 1). 2: Output: Length J sequence of approximations π N K and Z N K of π and Z. 3: for j = 1, . . . , J do 4: Sample X i 0 ∼ π 0 and set W i 0 = 1 N and Z N 0 = 1.

5:

for k = 1, . . . , K do 6:

Solve T k ← argmin T ∈T L N k (T ) using e.g. SGD.

7: ĥ ←flow-grad T k , π N k-1 using eqn (8).

8:

π N k , Z N k ← SMC-NF-step π N k-1 , Z N k-1 , T k 9:
Apply gradient update ĥ to flow T k 10:

end for 11:

Yield (π N K , Z N K ) and continue for loop. 12: end for

  denote the expectation of the random SMC approximation π N k of π k-1 w.r.t. to the law of the SMC-NF algorithm. An unbiased gradient of the objective KL[T # k πN k-1 ||π k ] is given by

Figure 1 .

 1 Figure 1. Comparison of normalizing constants between CRAFT and AFT for different numbers of particles. AFT is given 100 times more train/test MCMC steps than CRAFT. A gold standard value is shown in magenta.

Figure 2 .

 2 Figure 2. Timed comparison of normalizing constant estimates produced by CRAFT and SNF during training (higher is better). Three independent repeats are shown. A gold standard value is shown in magenta.

Figure 3 .

 3 Figure 3. Timed comparison of MCMC methods for the φ 4 example, based on fifteen repeats. CRAFT, SMC and VI serve as proposal mechanisms for Particle MCMC. HMC is applied directly to the target. Error is in estimating two point susceptibility, an example, physically relevant expectation. Note HMC never reaches the detailed level of error in the top row.

Figure 4 .

 4 Figure 4. Comparison of CRAFT, AFT, SMC and VI on the pines task from (Arbel et al., 2021) with a reduced number of HMC iterations.

Figure 5 .

 5 Figure 5. Timed comparison of MCMC methods for the φ 4 example, based on fifteen repeats. CRAFT, SMC and VI serve as proposal mechanisms for particle MCMC. HMC is applied directly to the target. Note HMC never reaches the detailed level of error in the top row.

Figure 6 .

 6 Figure 6. Five minutes of the raw chain used to compute the expectation averages for the SMC and CRAFT two point susceptibility. The estimated gold standard value is shown as the black dashed line. SMC and CRAFT are displayed because they are the two best algorithms. VI and HMC have much larger errors.

Figure 7 .

 7 Figure 7. Five minutes of the raw chain used to compute the expectation averages for the SMC and CRAFT Ising energy density. The estimated gold standard value is shown as the black dashed line. SMC and CRAFT are displayed because they are the two best algorithms. VI and HMC have much larger errors.

Continual

  Repeated Annealed Flow Transport Monte Carlo Algorithm 7 Practical AFT: from Arbel, Matthews and Doucet 2021. This corresponds to Algorithm 2 of the prior paper. 1: Input: number of training, test and validation particles N train , N test , N val , unnormalized annealed targets {γ k } K k=0 such that γ 0 = π 0 and γ K = γ, resampling thresholds A a ∈ [1/N a , 1) for a ∈ {train, test, val}, number of training iterations J. 2: Output: Approximations π Ntest K,test and Z Ntest K,test of π and Z. 3: for a ∈ {train, test, val} do 4:Sample X i,a 0 ∼ π 0 and set W i,a 0 ← 1 Na and Z N,a 0 ← 1. 5: end for 6: for k = 1, . . . , K do 7:Learn the flowT k ← EarlyStopLearnFlow(J, X Ntrain k-1,train , W Ntrain k-1,train , X Nval k-1,val , W Nval k-1,val ) 8:for a ∈ {train, test, val} do9: π Na k,a , Z Na k,a ← SMC-NF-step π Na k-1,a , Z Na k-1,a , T k 10:end for 11: end for Algorithm 8 Subroutine EarlyStopLearnFlow for practical AFT. 1: Input: Number of training iterations J, training and validation particles and weights X i,train kOuput: Estimated flow T k 3: Initialize flow to identity T k = I D . 4: Initialize list of flows T opt ← {T k }. 5: Initialize list of validation losses E ← Nval i=1 W i,val k-1 D k (X i,val k-1 ) 6: for j = 1, ..., J do using SGD step on L Ntrain k (T k ) with approx. gradient (56).

  flows T opt ← T opt ∪ {T k } flow with smallest validation error from the list of flows T opt .

Note that while G k (x k-1 ) corresponds to the Radon-Nikodym derivative between the measuresγ k (dx k )L col k-1 (x k , dx k-1 ) and γ k-1 (dx k-1 )M col k (x k-1 , dx k ), the ratio L col k-1 (x k , dx k-1 )/M col k (x k-1 , dx k )is not defined and taking its logarithm can only be interpreted as a shorthand notation for log G k (x k-1 ) -log(γ k (x k )/γ k-1 (x k-1 )).

Acknowledgements

The authors wish to thank Sébastien Racanière, Michalis Titsias, Phiala Shanahan, Jonas Köhler, Andy Ballard and the anonymous reviewers.

B. Gradient expression for ELBO based approaches

Assume we consider the ELBO objective as in [START_REF] Wu | Stochastic normalizing flows[END_REF], that is we are interested in maximizing w.r.t. to the NFs parameters and the stochastic kernel parameters L = E X 0:K ∼η K log w K (X 0:K-1 ) ,

where η(x 0:K ) = ηK (x 0:K ) is defined in (23) and, as shown in (26), we have

for πK (x 0:K ) defined in (24). It is thus trivial to check that w K (X 0:K-1 ) is an unbiased estimate of Z K when X 0:K ∼ ηK .

Thin et al. ( 2021) provide a gradient estimate of the ELBO for annealed importance sampling, which corresponds to the case where T l (x) = x for all l, when the MCMC kernels K l are of the Metropolis-Hastings type and we use a reparameterized proposal

We show here how one can directly extend these results to our settings. We will use the convention

By combining ( 21), ( 23), ( 31) and (32), we can rewrite the distribution ηk (x 0:K ) as

where, for a given realization of u 1:L , we draw sequentially the Bernoulli r.v. A l with

and, given x 0 , u 1:l , a 1:l , the state x l is given by

with Φ 0 (x 0 , u 1:0 , a 1:0 ) := x 0 . Finally we used the notation g(u

The gradient of this objective w.r.t. φ can be approximated as follows

where π N,φ K is an particle approximation of π φ K . This approximation can be obtained using a learned SMC/normalizing flow sampler obtained using CRAFT.

F. Further comparison of CRAFT and AFT

In this section we investigate the effect of reducing the number of MCMC iterations in CRAFT and AFT. As a benchmark we use the log Gaussian Cox process-with a large lattice discretization of 40 × 40 = 1600. To make the example harder we reduce the number of HMC iterations per temperature from 10 to 2. The rest of the experimental setup is kept the same. The goal here is to learn a fast flow augmented sampler at test time with reasonable training time. Since the diagonal affine flow used constitutes a negligible overhead the change corresponds to a 5 times reduction in test compute time relative to the previous work.

Figure 4 shows the results. We see that for 10 and 30 transitions that AFT struggles to learn a good flow. This is because the expectation of the objective required in the greedy training is not well estimated with fewer MCMC iterations. CRAFT has the same difficulty at the beginning of optimization but due to the repeated nature of the optimization it is able to recover.

Having seen the solution found by CRAFT one might reasonably ask if we could use different numbers of HMC iterations at train and test time for AFT-though this is not suggested in the original work of [START_REF] Arbel | Annealed flow transport Monte Carlo[END_REF]. After all, we know that with 10 iterations AFT performs well on the Cox process task, which suggests we could use that learnt flow with fewer HMC samples at test time. This is indeed the case. Training with extra MCMC iterations restores the behaviour of AFT to be similar to that of CRAFT in this case. However this is at the expense of introducing even more hyperparameters to AFT that require tuning. Relative to CRAFT, a user of AFT has to correctly tune the behaviour of three sets of particles (train/validation/test) and manage multiple separate optimization loops (one for each temperature). This extra modification of allowing different number of MCMC iterations is then a further complexity.

Note additionally that AFT is subject to the degradation in performance with batch size relative to CRAFT described in Section 4.1-we regard this case as the most direct demonstration of how CRAFT solves the sample replenishment problem in AFT.

G. Distributed and parallel computation

In this section we discuss the particular challenges we envisage scaling the algorithm to large scale compute infrastructure and potential remedies.

We have described the algorithm as a sequential algorithm but there is considerable scope for running it asynchronously in parallel. Each temperature can be considered to correspond to a node in a chain graph. Along the edges of the graph messages are passed. In the most direct implementation this corresponds to the particles and the weights. There is no need to run for a full pass along the graph before starting the next one. In some compute clusters communication locality is also important so it is helpful that the message passing is sparse and local.

An alternative to passing the particles and weights is to cycle the temperature at each node and pass gradient updates instead. The temperatures at each node can be offset by one so that the gradients are always received at the right time. Such a prescription is inspired by the distributed Algorithm 5 of [START_REF] Syed | Non-reversible parallel tempering: a scalable highly parallel MCMC scheme[END_REF]. For CRAFT, this method will reduce communication overhead between nodes in the case where the gradient updates require less memory than passing particles and weights. Although the gradient updates for neural networks can in general be large, incorporating symmetries as in Sections H and 4.3 or more generally adding domain knowledge can reduce network sizes substantially.

Another consideration is the case where only a few particles will fit onto a worker, so that it is desirable to parallelize a single iteration across multiple workers. The resampling step means that the algorithm is not batch parallel. Therefore resampling could introduce a significant communication overhead between the parallel workers. One possible remedy is only to resample within a worker.


Jinfeng Zhong
email: jinfeng.zhong@dauphine.euelsa.negre@lamsade.dauphine.fr

Elsa Negre

Amélioration des explications contrefactuelles pour les recommandations à l'aide de SHAP

Les explications dans les systèmes de recommandation aident les utilisateurs à mieux comprendre pourquoi de telles recommandations sont générées. Expliquer la recommandation est crucial pour renforcer la confiance et la satisfaction des utilisateurs. Étant donné que les systèmes de recommandation deviennent de plus en plus impénétrables, expliquer directement les recommandations devient parfois impossible. Les méthodes d'explication post-hoc qui n'élucident pas les mécanismes internes des systèmes de recommandation sont des approches populaires. Les méthodes d'explication post-hoc telles que SHAP génèrent des explications en construisant des modèles de substitution plus simples pour se rapprocher des modèles originaux. Cependant, l'application directe de telles méthodes pose plusieurs soucis : (1) Il se peut que les explications posthoc ne soient pas fidèles aux modèles de recommandation originaux puisque les mécanismes internes ne sont pas élucidés; (2) Les résultats retournés par des méthodes telles que SHAP ne sont pas triviales à comprendre pour la plupart des utilisateurs, car des connaissances mathématiques sont nécessaires. Dans ce travail, nous présentons une méthode d'explication des recommandations à l'aide de SHAP qui peut générer des explications contrefactuelles facilement compréhensibles avec une forte fidélité au modèle original.

Introduction

L'article résumé ici, intitulé "Shap-enhanced counterfactual explanations for recommendations" (Zhong et Negre, 2022a), a été accepté et présenté à la conférence internationale 37th ACM/SIGAPP Symposium on Applied Computing (SAC '22). Dans Zhong et Negre (2022a), nous nous intéressons à l'explicabilité dans les systèmes de recommandation (SRs). Plus spécifiquement, nous étudions des explications contrefactuelles à l'aide de SHAP [START_REF] Lundberg | A unified approach to interpreting model predictions[END_REF]. Les SRs sont de plus en plus déployés pour aider les utilisateurs à trouver des articles qui les intéressent. Par conséquent, les SRs sont devenus des outils d'aide à la décision importants et omniprésents. Cependant, les techniques actuelles des SRs deviennent de plus en plus impénétrables en raison de la complexité des modèles et de la grande dimensionnalité des données avec lesquelles les modèles fonctionnent. Des modèles "boîte noire" tels que les réseaux de neurones sont déployés pour générer des recommandations, ce qui soulève des inquiétudes concernant des propriétés indésirables telles que le manque de transparence et les préjugés d'illégalité [START_REF] Rudin | Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead[END_REF]. Dans certains domaines tels que la finance et la médecine où l'impact de mauvaises décisions ne peut être négligé, les utilisateurs s'attendent à ce que les SRs génèrent non seulement des recommandations de haute qualité, mais également des explications pour justifier de telles recommandations. Selon les modèles utilisés dans les SRs, diverses méthodes pour générer des explications ont été développées [START_REF] Zhang | Explainable recommendation : A survey and new perspectives[END_REF]. Cependant, l'application de ces méthodes est limitée à des techniques de recommandation spécifiques, ce qui limite leur adaptabilité, leur réutilisation et leur généricité. Cela signifie que chaque fois qu'une nouvelle technique de recommandation est développée, une méthode d'explication correspondante doit être développée. Ainsi, peut-on trouver une méthode générique pour expliquer pourquoi de telles recommandations sont générées ?

Les avancées récentes de l'Intelligence Artificielle (IA) explicable/interprétable éclairent ce sujet. Afin d'expliquer les sorties d'un système d'IA, deux stratégies sont possibles : (1) Adopter des modèles transparents, ce qui signifie que le fonctionnement des modèles soit compréhensible par les humains, e.g. les modèles linéaires simples et/ou les règles d'association ; (2) Fournir des explications post-hoc sans élucider précisément le fonctionnement des modèles. Les modèles "simples" sont faciles à comprendre, mais leurs applications sont limitées par une précision faible (compromis entre complexité et précision). LIME [START_REF] Ribeiro | why should i trust you ?" explaining the predictions of any classifier[END_REF] et SHAP [START_REF] Lundberg | A unified approach to interpreting model predictions[END_REF] sont, selon la littérature, des outils capables de fournir des explications post-hoc. Ils peuvent expliquer pourquoi une instance est proposée en estimant l'importance de chaque caractéristique de cette instance. Il existe des travaux qui expliquent les recommandations en calculant directement, grâce à LIME et/ou SHAP, l'importance de chacune de leurs caractéristiques. Les principales limites et préoccupations liées à l'application directe de telles méthodes incluent : (1) Il se peut que les explications post-hoc ne soient pas fidèles aux modèles de recommandations originaux et parfois ces explications sont contradictoires avec la recommandation ; (2) Pour certains utilisateurs, comprendre les résultats retournés par LIME et SHAP n'est pas trivial puisque des connaissances mathématiques plutôt complexes sont requises. Par exemple, les sorties de SHAP incluent la figure "force_plot" qui nécessite que les utilisateurs aient des connaissances en théorie des jeux. Par conséquent, LIME et SHAP sont de bons outils pour les développeurs de modèles pour visualiser les comportements des modèles et peuvent les aider à déboguer les modèles. Cependant, pour les utilisateurs, les importances/résultats retounés ne sont pas faciles à comprendre.

Comment utiliser ces méthodes d'explication post-hoc pour expliquer n'importe quel modèle de recommandation tout en évitant leurs limites ? Générer des explications contrefactuelles facilement compréhensibles et avec une forte fidélité au modèle de recommandations original est une solution potentielle. Contrairement à LIME et SHAP qui se rapprochent des modèles originaux, les explications contrefactuelles recherchent les changements minimaux nécessaires pour modifier les sorties des modèles originaux [START_REF] Wachter | Counterfactual explanations without opening the black box : Automated decisions and the gdpr[END_REF]. Par conséquent, les explications contrefactuelles ont une forte fidélité [START_REF] Kaffes | Model-agnostic counterfactual explanations of recommendations[END_REF]. Une explication contrefactuelle typique décrit une situation causale : "Si A ne s'est pas produit, alors B ne se serait pas produit" [START_REF] Molnar | Interpretable machine learning[END_REF]. Une méthode naïve pour rechercher des explications contrefactuelles consiste à énumérer tous les sous-ensembles de caractéristiques d'une instance [START_REF] Martens | Explaining data-driven document classifications[END_REF]. Une préoccupation majeure est que la complexité est exponentielle par rapport au nombre de caractéristiques utilisées dans le modèle original. Dans Zhong et Negre (2022a), nous montrons que SHAP peut être adapté pour générer des explications contrefactuelles sans énumérer tous les sous-ensembles de caractéristiques d'une instance. Plus de détails sur LIME et SHAP seront présentés dans la section 2.

LIME et SHAP

Dans cette section, nous comparons les deux méthodes et montrons pourquoi nous avons choisi SHAP.

Les similarités entre LIME et SHAP comprennent : (1) LIME et SHAP peuvent générer des explications post-hoc pour une seule instance à expliquer sans avoir à élucider les mécanismes internes des modèles. Par conséquent, LIME et SHAP peuvent être utilisés pour expliquer tous les modèles ; (2) LIME et SHAP fournissent une visualisation intuitive qui montre la contribution de chaque caractéristique, ce qui aide les développeurs à comprendre le fonctionnement de modèles complexes. Il convient de noter que pour les utilisateurs ordinaires, il n'est pas facile de comprendre les résultats obtenus par LIME et SHAP, car des connaissances mathématiques sont nécessaires.

Cependant, à chaque fois que LIME est appliqué, il sélectionne un nombre limité de caractéristiques et génère de nouveaux points de données. Par conséquent, il se peut que les modèles linéaires (par exemple, la régression de Ridge) générés soient différents (c'est-à-dire les caractéristiques utilisées et les coefficients renvoyés), ce qui soulève des inquiétudes quant à la stabilité de LIME [START_REF] Visani | Statistical stability indices for lime : obtaining reliable explanations for machine learning models[END_REF][START_REF] Zafar | Dlime : A deterministic local interpretable modelagnostic explanations approach for computer-aided diagnosis systems[END_REF]. SHAP est fondamentalement fondé sur la théorie des jeux donc SHAP assure un calcul équitable de la contribution de chaque caractéristique. Par conséquent, les propriétés souhaitées telles que la précision et la cohérence locales sont garanties [START_REF] Lundberg | A unified approach to interpreting model predictions[END_REF]. Finalement, nous avons choisi SHAP pour générer des explications post-hoc.

Nous invitons les lecteurs à se référer à notre article Zhong et Negre (2022a) ainsi que LIME [START_REF] Ribeiro | why should i trust you ?" explaining the predictions of any classifier[END_REF] et SHAP [START_REF] Lundberg | A unified approach to interpreting model predictions[END_REF] pour plus de détails sur la comparaison entre eux.

Notre proposition

Selon SHAP, le "jeu" est la prédiction et les "joueurs" sont les caractéristiques que le modèle de recommandation prend en compte. SHAP calcule la contribution de chaque caractéristique, dans la suite de cet article, nous notons la contribution de la caractéristique j comme ϕ j . Cependant, interpréter ϕ j comme une déclaration contrefactuelle n'est pas une manière judicieuse de résumer l'information puisque ϕ j est la contribution marginale moyenne de la caractéristique j [START_REF] Kumar | Problems with shapley-value-based explanations as feature importance measures[END_REF]. Bien qu'il ne soit pas approprié d'interpréter ϕ j comme une déclaration contrefactuelle, SHAP permet de guider la recherche d'explications contrefactuelles.

Nous proposons l'algorithme 1 pour générer des explications contrefactuelles. L'idée principale est d'appliquer SHAP pour sélectionner les caractéristiques importantes i.e. celles pour lesquelles ϕ j est positive. Ensuite, la valeur de la caractéristique qui est classée en tête est modifiée. Dans l'étape 13, x j est remplacé par x new j . Une autre liste de recommandations {i ′ 1 , i

  Generate counterfactual explanations using SHAP Input : X = {x 1 , x 2 , x 3 , . . . , x p }, information used to generate a list of recommendations ; f , recommender model ; i, a recommendation to be explained Output : Set of original features values E 1 to be changed, Set of counterfactual feature values E_counter 2 E_to_expand = {} % The indices of candidate features 3 E = {} 4 E_counter = {} 5 Stop = False % Indication of termination of algorithm 6 Generate recommendations {i 1 , i 2 , i 3 . . . , i max } for user u 7 Compute the contribution of each features in X using SHAP :rui = ϕ 0 + ϕ 1 + ϕ 2 + • • • + ϕ p 8 I p ⊂ X : {I 1 , I 2 , . . . },indices of features whose ϕ Ip >= 0 9 Sort I p according to ϕ Ip in a descending manner 10 for l ← 1 to length(I p ) do 11 E_to_expand = I p [: l] % The l first indices

	Algorithme 1 : 12 for j ∈ E_to_expand do
	13	Change x j to x new j
	14	Generate a new list of recommendations {i ′ 1 , i ′ 2 , i ′ 3 , . . . , i ′ max } given x new = {x 1 , x 2 , x new j , . . . , x p }
	15	if i / ∈ {i ′ 1 , i ′ 2 , i ′ 3 , . . . , i ′ max } then
	16	E = E ∪ x j
	17	E_counter = E_counter ∪ x new j
	18	Stop = True
	19	break
	20	if Stop = True then
	21	break
	22 if i ∈ {i ′ 1 , i ′ 2 , i ′ 3 , . . . , i ′ max } then
	23	E = {}
		′ 2 , i ′ 3 , . . . , i ′ max } est calculée étant donné x new = {x 1 , x 2 , x new j	. . . x p }. Si i n'est plus

24 return E, E_counter

quantifiée via le coefficient de chevauchement (Overlap Coefficient) (Vijaymeena et Kavitha, 2016) entre eux. Nous appliquons l'algorithme 1 pour expliquer un modèle de recommandation contextuelle que nous avons récemment proposé (Zhong et Negre, 2022b). Nous montrons que notre méthode est capable d'expliquer plus de 95% des recommandations, montrant ainsi une grande fidélité du modèle. Par rapport à d'autres méthodes, notre méthode peut réduire le temps de recherche car elle permet de réduire l'espace de recherche. De plus, les explications générées par notre méthode sont fidèles au modèle de recommandation original en identifiant les conditions contextuelles clés. Un autre avantage est que ces explications sont adaptées aux situations contextuelles des utilisateurs, ce qui pourrait être plus convaincant.

En résumé, les contributions de ce travail comprennent : (1) Nous commençons par résumer les travaux de l'état de l'art, nous montrons que les explications dans les SRs doivent être contrastives et sélectionnées ; (2) Nous comparons également LIME et SHAP, cela nous permet de préciser pourquoi nous avons choisi SHAP plutôt que LIME dans ce travail ; (3) Nous proposons ensuite un algorithme capable de générer des explications contrefactuelles dans les SRs à l'aide de SHAP. À notre connaissance, il s'agit du premier article qui combine la méthode de la valeur de Shapley et la méthode contrefactuelle pour générer des explications dans les SRs ; (4) Parallèlement, nous montrons que les explications contrefactuelles générées par notre algorithme sont contrastives et sélectionnées, deux propriétés privilégiées pour les explications en IA ; (5) Enfin, nous présentons une étude de cas pour montrer l'efficacité de notre algorithme. Nos expériences montrent que SHAP, sous certaines conditions, peut réduire le temps nécessaire à la recherche d'explications contrefactuelles.

dans la liste {i ′ 1 , i ′ 2 , i ′ 3 , . . . , i ′ max } alors x j est une explication contrefactuelle. Si i est encore dans la liste {i ′ 1 , i ′ 2 , i ′ 3 , . . . , i ′ max } après avoir essayé toutes les possibilités de x new j , alors les deux premières caractéristiques classées sont examinées, et ainsi de suite, jusqu'à ce que toutes les caractéristiques de I p soient explorées. Enfin, si l'ensemble des valeurs des caractéristiques originales E est vide, alors cela signifie que cette recommandation ne peut pas être expliquée par l'algorithme 1. Il est à noter qu'il est également possible d'expliquer une liste de recommandations avec notre algorithme 1. Dans ce cas, i devient une liste de recommandations {i 1 , i 2 , i 3 , . . . , i max } ; dans l'étape 15, la condition devient :

Summary

Explanations in recommender systems help users better understand why recommendations are generated, which is crucial for enhancing users' trust and satisfaction. As recommender systems become ever more inscrutable, directly explaining recommender systems sometimes becomes impossible. Post-hoc explanation methods that do not elucidate internal mechanisms of recommender systems are popular approaches. State-of-art post-hoc explanation methods such as SHAP can generate explanations by building simpler surrogate models to approximate the original models. However, directly applying such methods has several concerns: (1) Posthoc explanations may not be faithful to the original recommender systems since the internal mechanisms of recommender systems are not elucidated; (2) The outputs returned by methods such as SHAP are not trivial for plain users to understand since background mathematical knowledge is required. In this work, we present an explanation method enhanced by SHAP that can generate easily understandable counterfactual explanations with high fidelity.