Julia Linhart
email: julia.linhart@inria.fr

Alexandre Gramfort
email: alexandre.gramfort@inria.fr

Pedro L C Rodrigues
email: pedro.rodrigues@inria.fr

L-C2ST: Local Diagnostics for Posterior Approximations in Simulation-Based Inference

Many recent works in simulation-based inference (SBI) rely on deep generative models to approximate complex, high-dimensional posterior distributions. However, evaluating whether or not these approximations can be trusted remains a challenge. Most approaches evaluate the posterior estimator only in expectation over the observation space. This limits their interpretability and is not sufficient to identify for which observations the approximation can be trusted or should be improved. Building upon the well-known classifier two-sample test (C2ST), we introduce ℓ-C2ST, a new method that allows for a local evaluation of the posterior estimator at any given observation. It offers theoretically grounded and easy to interpret -e.g. graphical -diagnostics, and unlike C2ST, does not require access to samples from the true posterior. In the case of normalizing flow-based posterior estimators, ℓ-C2ST can be specialized to offer better statistical power, while being computationally more efficient. On standard SBI benchmarks, ℓ-C2ST provides comparable results to C2ST and outperforms alternative local approaches such as coverage tests based on highest predictive density (HPD). We further highlight the importance of local evaluation and the benefit of interpretability of ℓ-C2ST on a challenging application from computational neuroscience.

Introduction

Expressive simulators are at the core of modern experimental science, enabling the exploration of rare or challenging-to-measure events in complex systems across various fields such as population genetics [START_REF] Tavaré | Inferring coalescence times from DNA sequence data[END_REF], astrophysics [START_REF] Dax | Group equivariant neural posterior estimation[END_REF], cosmology [START_REF] Lucas Makinen | Lossless, scalable implicit likelihood inference for cosmological fields[END_REF], and neuroscience [START_REF] Sanz Leon | The virtual brain: a simulator of primate brain network dynamics[END_REF][START_REF] Gonçalves | Training deep neural density estimators to identify mechanistic models of neural dynamics[END_REF][START_REF] Bittner | Interrogating theoretical models of neural computation with emergent property inference[END_REF][START_REF] Hashemi | Simulation-based inference for whole-brain network modeling of epilepsy using deep neural density estimators[END_REF]. These simulators implicitly encode the intractable likelihood function p(x | θ) of the underlying mechanistic model, where θ represents a set of relevant parameters and x ∼ Simulator(θ) is the corresponding realistic observation. The main objective is to infer the parameters associated with a given observation using the simulator's posterior distribution p(θ | x) [START_REF] Cranmer | The frontier of simulation-based inference[END_REF]. However, classical methods for sampling posterior distributions, such as MCMC [START_REF] Robert | Monte Carlo statistical methods[END_REF] and variational inference [START_REF] Rezende | Variational inference with normalizing flows[END_REF], rely on the explicit evaluation of the model-likelihood, which is not possible when working with most modern simulators.

Simulation-based inference (SBI) [START_REF] Cranmer | The frontier of simulation-based inference[END_REF] addresses this problem by estimating the posterior distribution on simulated data from the joint distribution. This can be done after choosing a prior distribution p(θ) over the parameter space and using the identity p(θ, x) = p(x | θ)p(θ). In light of recent developments in the literature on deep generative models, different families of algorithms have been proposed to approximate posterior distributions in SBI [START_REF] Cranmer | The frontier of simulation-based inference[END_REF]. Certain works use normalizing flows [START_REF] Papamakarios | Normalizing flows for probabilistic modeling and inference[END_REF] to directly learn the posterior density function (neural posterior estimation, NPE [START_REF] Greenberg | Automatic posterior transformation for likelihood-free inference[END_REF]) or aim for the likelihood (neural likelihood estimation, NLE [START_REF] Papamakarios | Sequential neural likelihood: Fast likelihood-free inference with autoregressive flows[END_REF]). Other approaches reframe the problem in terms of a classification task and aim for likelihood ratios (neural ratio estimation, NRE [START_REF] Hermans | Likelihood-free MCMC with amortized approximate ratio estimators[END_REF]). However, appropriate validation remains a challenge for all these paradigms, and principled statistical approaches are still needed before SBI can become a trustworthy technology for experimental science. This topic has been the goal of many recent studies. For instance, certain proposals aim at improving the posterior estimation by preventing over-confidence [START_REF] Hermans | A crisis in simulation-based inference? beware, your posterior approximations can be unfaithful[END_REF] or addressing model misspecification [START_REF] Frazier | Model misspecification in approximate Bayesian computation: consequences and diagnostics[END_REF] to ensure conservative [START_REF] Delaunoy | Towards reliable simulation-based inference with balanced neural ratio estimation[END_REF] and more robust posterior estimators [START_REF] Ward | Robust neural posterior estimation and statistical model criticism[END_REF][START_REF] Lemos | Robust simulation-based inference in cosmology with bayesian neural networks[END_REF]. Another approach is the development of a SBI benchmark [START_REF] Lueckmann | Benchmarking simulation-based inference[END_REF] for comparing and validating different algorithms on many standard tasks. While various validation metrics exist, Lueckmann et al. [START_REF] Lueckmann | Benchmarking simulation-based inference[END_REF] show that, overall, classifier two sample tests (C2ST) [START_REF] Lopez-Paz | Revisiting classifier two-sample tests[END_REF] are currently the most powerful and flexible approach. Based on standard methods for binary classification, they can scale to high-dimensions as well as handle non-Euclidean data spaces [START_REF] Kim | Global and local two-sample tests via regression[END_REF][START_REF] Pandeva | E-valuating classifier two-sample tests[END_REF]. Typical use-cases include tests for statistical independence and the evaluation of sample quality for generative models [START_REF] Lopez-Paz | Revisiting classifier two-sample tests[END_REF]. Implicitly, C2ST is used in algorithms such as noise contrastive estimation [START_REF] Michael | Noise-contrastive estimation of unnormalized statistical models, with applications to natural image statistics[END_REF] and generative adversarial networks [START_REF] Goodfellow | Generative adversarial networks[END_REF], or to estimate likelihood-to-evidence ratios [START_REF] Hermans | Likelihood-free MCMC with amortized approximate ratio estimators[END_REF]. To be applied in SBI settings, however, C2ST requires access to samples from the true target posterior distribution, which renders it useless in practice. Simulationbased calibration (SBC) [START_REF] Talts | Validating bayesian inference algorithms with simulation-based calibration[END_REF] bypasses this issue by only requiring samples from the joint distribution. Implemented in standard packages of the field (sbi [START_REF] Tejero-Cantero | sbi: A toolkit for simulationbased inference[END_REF], Stan [START_REF] Carpenter | Stan: A probabilistic programming language[END_REF]), it has become the go-to validation method for SBI [START_REF] Lueckmann | Benchmarking simulation-based inference[END_REF][START_REF] Hermans | A crisis in simulation-based inference? beware, your posterior approximations can be unfaithful[END_REF] and has been further studied in recent works [START_REF] Modrák | Simulation-based calibration checking for bayesian computation: The choice of test quantities shapes sensitivity[END_REF][START_REF] Dalmasso | Confidence sets and hypothesis testing in a likelihood-free inference setting[END_REF][START_REF] Gelman | Bayesian workflow[END_REF]. Coverage tests based on the highest predictive density (HPD) as used in [START_REF] Hermans | A crisis in simulation-based inference? beware, your posterior approximations can be unfaithful[END_REF][START_REF] Delaunoy | Towards reliable simulation-based inference with balanced neural ratio estimation[END_REF], can be seen as a variant of SBC that are particularly well adapted to multivariate data distribution.

Nevertheless, a big limitation of current SBI diagnostics remains: they only evaluate the quality of the posterior approximation globally (in expectation) over the observation space and fail to give any insight of its local behavior. This hinders interpretability and can lead to false conclusions on the validity of the estimator [START_REF] Zhao | Diagnostics for conditional density models and bayesian inference algorithms[END_REF][START_REF] Linhart | Validation diagnostics for SBI algorithms based on normalizing flows[END_REF]. There have been attempts to make existing methods local, such as local-HPD [START_REF] Zhao | Diagnostics for conditional density models and bayesian inference algorithms[END_REF] or local-multi-PIT [START_REF] Linhart | Validation diagnostics for SBI algorithms based on normalizing flows[END_REF], but they depend on many hyper-parameters and are computationally too expensive to be used in practice. In this work, we present ℓ-C2ST, a new local validation procedure based on C2ST that can be used to evaluate the quality of SBI posterior approximations for any given observation, without using any data from the target posterior distribution. ℓ-C2ST comes with necessary, and sufficient, conditions for the local validity of multivariate posteriors and is particularly computationally efficient when applied to validate NPE with normalizing flows, as often done in SBI literature [START_REF] Dax | Group equivariant neural posterior estimation[END_REF][START_REF] Gonçalves | Training deep neural density estimators to identify mechanistic models of neural dynamics[END_REF][START_REF] Lemos | Robust simulation-based inference in cosmology with bayesian neural networks[END_REF][START_REF] Ward | Robust neural posterior estimation and statistical model criticism[END_REF][START_REF] Dax | Real-time gravitational wave science with neural posterior estimation[END_REF][START_REF] Vasist | Neural posterior estimation for exoplanetary atmospheric retrieval[END_REF][START_REF] Jallais | Inverting brain grey matter models with likelihood-free inference: a tool for trustable cytoarchitecture measurements[END_REF]. Furthermore, ℓ-C2ST offers graphical tools for analysing the inconsistency of posterior approximations, showing in which regions of the observation space the estimator should be improved and how to act upon, e.g. signs of positive / negative bias, signs of over / under dispersion, etc.

In what follows, we first introduce the SBI framework and review the basics of C2ST. Then, we detail the ℓ-C2ST method and prove asymptotic theoretical guarantees. Finally, we report empirical results on two SBI benchmark examples to analyze the performance of ℓ-C2ST and a non-trivial neuroscience use-case that showcases the need of a local validation method.

Validating posterior approximations with classifiers

Consider a model with parameters θ ∈ R m and observations x ∈ R d obtained via a simulator. In what follows, we will always assume the typical simulation-based inference setting, meaning that the likelihood function p(x | θ) of the model cannot be easily evaluated. Given a prior distribution p(θ), it is possible to generate samples from the joint pdf p(θ, x) as per:

Θ n ∼ p(θ) ⇒ X n = Simulator(Θ n ) ∼ p(x | Θ n ) ⇒ (Θ n , X n ) ∼ p(θ, x) . (1)
Let N s be a fixed simulation budget and

{(Θ n , X n )} Ns n=1 = D train ∪ D cal with D train ∩ D cal = ∅.
The data from D train are used to train an amortized2 approximation q(θ | x) ≈ p(θ | x), e.g. via NPE [START_REF] Greenberg | Automatic posterior transformation for likelihood-free inference[END_REF], and those from D cal to diagnose its local consistency [START_REF] Zhao | Diagnostics for conditional density models and bayesian inference algorithms[END_REF].

Definition 1 (Local consistency).

A conditional density estimator q is said to be locally consistent at x o with the true posterior density p if, and only if, the following null hypothesis holds:

H 0 (x o ) : q(θ | x o ) = p(θ | x o ), ∀θ ∈ R m . (2)
We can reformulate H 0 (x o ) as a binary classification problem by partitioning the parameter space into two balanced classes: one for samples from the approximation (C = 0) and one for samples from the true posterior (C = 1), as in

Θ | (C = 0) ∼ q(θ | x o ) vs. Θ | (C = 1) ∼ p(θ | x o ) , (3)
for which the optimal Bayes classifier [START_REF] Hastie | The elements of statistical learning[END_REF] is

f ⋆ xo (θ) = argmax 1 -d ⋆ xo (θ), d ⋆ xo (θ) with d ⋆ xo (θ) = P(C = 1 | Θ = θ; x o ) = 1 -P(C = 0 | Θ = θ; x o ) = p(θ|xo) p(θ|xo)+q(θ|xo) . (4)
It is a standard result [START_REF] Lopez-Paz | Revisiting classifier two-sample tests[END_REF][START_REF] Kim | Global and local two-sample tests via regression[END_REF] to relate (2) with (3) as per

H 0 (x o ) holds ⇐⇒ d ⋆ xo (θ) = P(C = 1 | Θ = θ; x o ) = 1 2 ∀θ
(5) When the classes are non-separable, the optimal Bayes classifier will be unable to make a decision and we can assume that it behaves as a Bernoulli random variable [START_REF] Lopez-Paz | Revisiting classifier two-sample tests[END_REF].

Classifier Two-Sample Test (C2ST)

The original version of C2ST [START_REF] Lopez-Paz | Revisiting classifier two-sample tests[END_REF] uses [START_REF] Dalmasso | Confidence sets and hypothesis testing in a likelihood-free inference setting[END_REF] to define a test statistic for H 0 (x o ) based on the accuracy of a classifier f xo trained on a dataset defined as

Θ q n ∼ q(θ | x o ) C=0 and Θ p n ∼ p(θ | x o ) C=1 and D = {(Θ q n , 0)} N n=1 ∪ {(Θ p n , 1)} N n=1 . (6)
The classifier accuracy is then empirically estimated over 2N v samples (N v samples in each class) from a held-out validation dataset D v generated in the same way as D:

tAcc (f xo ) = 1 2N v 2Nv n=1 I f xo (Θ q n ) = 0 + I f xo (Θ p n ) = 1 . (7)
Theorem 1 (Local consistency and classification accuracy). If f xo is Bayes optimal 3 and N v → ∞, then tAcc (f xo ) = 1/2 is a necessary and sufficient condition for the local consistency of q at x o .

See Appendix A.1 for a proof. The intuition is that, under the null hypothesis H 0 (x o ), it is impossible for the optimal classifier to distinguish between the two data classes, and its accuracy will remain at chance-level [START_REF] Lopez-Paz | Revisiting classifier two-sample tests[END_REF]. In the context of SBI, C2ST has been used to benchmark a variety of different procedures on toy examples where the true posterior is known and can be sampled [START_REF] Lueckmann | Benchmarking simulation-based inference[END_REF]. This is why we call this procedure an oracle C2ST, since it uses information that is not available in practice.

Regression C2ST. Kim et al. [START_REF] Kim | Global and local two-sample tests via regression[END_REF] argues that the usual C2ST based on the classifier's accuracy may lack statistical power because of the "binarization" of the posterior class probabilities. They propose to instead use probabilistic classifiers (e.g. logistic regression) of the form

f xo (θ) = I d xo (θ) > 1 2 where d xo (θ) = P(C = 1 | θ; x o ) (8
) and define the test statistics in terms of the predicted class probabilities d xo instead of the predicted class labels. The test statistic is then the mean squared distance between the estimated class posterior probability and one half:

tMSE (f xo ) = 1 N v Nv n=1 d xo (Θ q n ) - 1 2 2 + 1 N v Nv n=1 d xo (Θ p n ) - 1 2 2 (9)
Theorem 2 (Local consistency and regression). If f xo is Bayes optimal and N v → ∞, then tMSE (f xo ) = 0 is a necessary and sufficient condition for the local consistency of q at x o .

See Appendix A.2 for a proof. The numerical illustrations in Kim et al. [START_REF] Kim | Global and local two-sample tests via regression[END_REF] give empirical evidence that Regression C2ST has superior statistical power as compared to its accuracy-based counterpart, particularly for high-dimensional data spaces. Furthermore, it offers tools for interpretation and visualization: evaluating the predicted class probabilities d xo (θ) for any θ ∈ R m informs the regions where the classifier is more (or less) confident about its choice [START_REF] Lopez-Paz | Revisiting classifier two-sample tests[END_REF][START_REF] Kim | Global and local two-sample tests via regression[END_REF].

ℓ-C2ST: Local Classifier Two-Sample Tests

The oracle C2ST framework is not applicable in practical SBI settings, since it requires access to samples from the true posterior distribution to (1) train a classifier and (2) evaluate its performance in discriminating data from q and p. This section presents a new method called local C2ST (ℓ-C2ST) capable of evaluating the local consistency of a posterior approximation requiring data only from the joint pdf p(θ, x) which can be easily sampled as per [START_REF] Bittner | Interrogating theoretical models of neural computation with emergent property inference[END_REF].

(1) Train the classifier. We define a modified version of the classification framework (3) with:

(Θ, X) | (C = 0) ∼ q(θ | x)p(x) vs. (Θ, X) | (C = 1) ∼ p(θ, x) . (10)
The optimal Bayes classifier is now

f ⋆ (θ, x) = argmax {1 -d ⋆ (θ, x), d ⋆ (θ, x)} with d ⋆ (θ, x) = p(θ, x) p(θ, x) + q(θ | x)p(x) = p(θ | x) p(θ | x) + q(θ | x) = d ⋆ x (θ) , (11)
where one can notice the direct relation with the Bayes classifier for (3). Therefore, using data sampled as in [START_REF] Dheur | A large-scale study of probabilistic calibration in neural network regression[END_REF], it is possible to train a classifier f (θ, x) and write

f xo (θ) = f (θ, x o ) for each x o .
See Algorithm 1 for details on the implementation of this procedure.

(2) Evaluate the classifier. Define a new test statistic that evaluates the MSE-statistic for a classifier f and its associated predicted probabilities d using data samples from only the class associated to the posterior approximation (C = 0):

tMSE0 (f, x o ) = 1 N v Nv n=1 d(Θ q n , x o ) - 1 2 2 with Θ q n ∼ q(θ | x o ) . ( 12
)
Theorem 3 (Local consistency and single class evaluation). If f is Bayes optimal and N v → ∞, then tMSE0 (f, x o ) = 0 is a necessary and sufficient condition for the local consistency of q at x o . Proof. Let d be an estimator of P(C = 1 | Θ, X) and f = I(d > 0.5). Suppose that f = f ⋆ is Bayes optimal and let x o be a fixed observation. We have that

lim Nv→∞ tMSE0 (f ⋆ , x o ) = d ⋆ xo (θ) - 1 2 2 q(θ | x o )dθ .
Because of the squared term in the integral and q being a p.d.f., we have that

lim Nv→∞ tMSE0 (f ⋆ , x o ) = 0 ⇐⇒ d ⋆ xo (θ) = P(C = 1 | θ; x o ) = 1 2 .
This new statistical test can thus be used to assess the local consistency of posterior approximation q without using any sample from the true posterior distribution p, but only from the joint pdf. Furthermore, it is amortized, so a single classifier is trained for (10) that can then be used for any choice of conditioning observation x o . This is not the case in the usual oracle C2ST framework. QDA-C2ST (power)

t MSE t MSE0 2
Figure 1: Results for the C2ST framework when p = N (0, I 2 ) and q = N (0, σ 2 I 2 ). Left panel portrays the test statistics for the optimal Bayes classifier and the right panel shows the test's empirical power with QDA. Single-class accuracy test ( tAcc0 ) fails to detect when p ̸ = q but tMSE0 behaves correctly.

Figure 1 illustrates the behavior of different test statistics to discriminate samples from two bivariate Normal distributions whose covariance mismatch is controlled by a single scaling parameter σ. Note that the optimal Bayes classifier for this setting can be obtained via quadratic discriminant analysis (QDA) [START_REF] Hastie | The elements of statistical learning[END_REF]. The results clearly show that even though t MSE0 exploits only half of the dataset (i.e. samples from class C = 0) it is capable of detecting when p and q are different (σ ̸ = 1). The plot also includes the results for a one-class test statistic based on accuracy values ( tAcc0 ) which, as opposed to tMSE0 , has no guarantees for being a necessary and sufficient condition for local consistency. Not surprisingly, it fails to reject the null hypothesis for various choices of σ.

The assumptions of Theorem 3 are never met in practice: datasets are finite and one rarely knows which type of classifier is optimal for a given problem. Therefore, the values of tMSE0 in the null hypothesis (p = q) tend to fluctuate around one-half, and it is essential to determine a threshold for deciding whether or not H 0 (x o ) should be rejected. In ℓ-C2ST, these threshold values are obtained via a permutation procedure [START_REF] Efron | Institute of mathematical statistics monographs: Computer age statistical inference: Algorithms, evidence, and data science series number 5[END_REF] described in Algorithm 1. This yields N H estimates of the test statistic under the null hypothesis and can be used to calculate p-values for any given α significance level as described in Algorithm 2. These estimates can also be used to form graphical summaries known as PP-plots, which display the empirical CDF of the probability predictions versus the nominal probability level. These plots show how the predicted class probability d(x o ) deviates from its theoretical value under the null hypothesis (i.e. one half) as well as (1 -α) confidence regions; see Algorithm B.1 available in the appendix for more details and Figure 4 for an example.

/ for n = 1, . . . , N cal do Θ q n ∼ q(θ | X n ) W 2n = (Θ q n , X n ); C 2n = 0 /* Sample from q(θ | x)p(x) */ W 2n+1 = (Θ n , X n ); C 2n+1 = 1 /* Sample from p(θ, x) */ D ← {W n , C n } 2N cal n=1 /*
); p-value p(x o ) Generate N v samples Θ q n ∼ q(θ | x o ) with predicted probabilities d(Θ q n , x o ) and d h (Θ q n , x o ) /* Compute test statistics */ tMSE0 (x o ) ← 1 Nv n d(Θ q n , x o ) -1 2 2 for h = 1, . . . , N H do th (x o ) ← 1 Nv n d h (Θ q n , x o ) -1 2 2 /* Compute p-value */ p(x o ) ← 1 N H h I th (x o ) > tMSE0 (x o ) return tMSE0 (x o ), p(x o )

The case of normalizing flows

The ℓ-C2ST framework can be further improved when the posterior approximation q is a conditional normalizing flow [START_REF] Papamakarios | Normalizing flows for probabilistic modeling and inference[END_REF], which we denote q ϕ . Given a Gaussian base distribution u(z) = N (0, I m ) and a bijective transform T ϕ (•; x) with Jacobian J T ϕ (•; x) we have

q ϕ (θ | x) = u(z)| det J T ϕ (z; x)| -1 , θ = T ϕ (z; x) ∈ R m . (13)
In other words, normalizing flows (NF) are invertible neural networks that define a map between a latent space where data follows a Gaussian distribution and the parameter space containing complex posterior distributions. This allows for both efficient sampling and density evaluation:

Z ∼ N (0, I m ) ⇒ Θ q = T ϕ (Z; x) ∼ q ϕ (θ | x) , (14)
q ϕ (θ | x) = u(T -1 ϕ (θ; x))| det J T -1 ϕ (θ; x)| . (15)
Our main observation is that the inverse transform T -1 ϕ can also be used to characterize the local consistency of the conditional normalizing flow in its latent space, yielding a much simpler and computationally less expensive statistical test for posterior local consistency. Theorem 4 (Local consistency and normalizing flows). Given a posterior approximation q ϕ based on a normalizing flow, its local consistency at x o can be characterized as follows:

p(θ | x o ) = q ϕ (θ | x o ) ⇐⇒ p(T -1 ϕ (θ; x o ) | x o ) = u(z) , ∀θ ∈ R m . ( 16
) Proof. Let Θ ∼ p(θ | x o ). Following (14), we have that Θ ∼ q ϕ (θ | x o ) if, and only if, Θ = T ϕ (Z; x o ) with Z ∼ N (0, I m ).
Applying the inverse transformation of the flow gives us

T -1 ϕ (Θ; x o ) = T -1 ϕ (T ϕ (Z; x o ); x o ) = Z ∼ N (0, I m )
, which concludes the proof. Based on Theorem 4 we propose a modified version of our statistical test named ℓ-C2ST-NF. The new null hypothesis associated with the consistency of the posterior approximation q ϕ at x o is

H NF 0 (x o ) : p(T -1 ϕ (θ; x o ) | x o ) = N (0, I m ) , (17)
which leads to a new binary classification framework

(Z, X) | (C = 0) ∼ N (0, I m )p(x) vs. (Z, X) | (C = 1) ∼ p(T -1 ϕ (θ; x), x) . (18)
Algorithm 3 describes how to sample data from each class and train a classifier to discriminate them.

The classifier is then evaluated on N v samples Z n ∼ N (0, I m ) which are independent of x o .

A remarkable feature of ℓ-C2ST-NF is that calculating the test statistics under the null hypothesis is considerably faster than for ℓ-C2ST. In fact, for each null trial h = 1, . . . , N H we use the dataset D cal only for recovering the samples X n and then independently sample new data Z n ∼ N (0, I m ). As such, it is possible to pre-train the classifiers without relying on a permutation procedure (cf. Algorithm 4), and to re-use them to quickly compute the validation diagnostics for any choice of x o or new posterior estimator q ϕ of the given inference task. It is worth mentioning that this is not possible with the usual ℓ-C2ST, as it depends on q and it would require new simulations for each trial.

Algorithm 3: ℓ-C2ST-NF -training the classifier on the joint distribution Input: NF posterior estimator q ϕ ; calibration data

D cal = {Θ n , X n } N cal n=1 ; classifier f Output: estimate d of the class probabilities /* Construct classification training set */ for n in 1, . . . , N cal do Z n ∼ N (0, I m ); Z q n = T -1 ϕ (Θ n ; X n ) /* inverse NF-transformation */ W 2n = (Z n , X n ); C 2n = 0 W 2n+1 = (Z q n , X n ); C 2n+1 = 1 D ← {W n , C n } 2N cal n=1 /*
Z n ∼ N (0, I m ) for n = 1, . . . , 2N cal D ← {(Z 2n , X n ), 0} N cal n=1 ∪ {(Z 2n+1 , X n ), 1} N cal n=1 /*

Experiments

All experiments were implemented with Python and the sbi package [START_REF] Tejero-Cantero | sbi: A toolkit for simulationbased inference[END_REF] combined with Py-Torch [START_REF] Paszke | PyTorch: An Imperative Style, High-Performance Deep Learning Library[END_REF] and nflows [START_REF] Durkan | nflows: normalizing flows in PyTorch[END_REF] for neural posterior estimation 4 . Classifiers on the C2ST framework use the MLPClassifier from scikit-learn [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF] with the same parameters as those used in sbibm [START_REF] Lueckmann | Benchmarking simulation-based inference[END_REF].

Two benchmark examples for SBI

We illustrate ℓ-C2ST on two examples: Two Moons and SLCP. These models have been widely used in previous works from SBI literature [START_REF] Greenberg | Automatic posterior transformation for likelihood-free inference[END_REF][START_REF] Papamakarios | Sequential neural likelihood: Fast likelihood-free inference with autoregressive flows[END_REF] and are part of the SBI benchmark [START_REF] Lueckmann | Benchmarking simulation-based inference[END_REF]. They both represent difficult inference tasks with locally structured multimodal true posteriors in, respectively, low (θ ∈ R 2 , x ∈ R 2 ) and high (θ ∈ R5 , x ∈ R 8 ) dimensions. See [START_REF] Lueckmann | Benchmarking simulation-based inference[END_REF] for more details. To demonstrate the benefits of ℓ-C2ST-NF, all experiments use neural spline flows [START_REF] Durkan | Neural spline flows[END_REF] trained under the amortized paradigm for neural posterior estimation (NPE) [START_REF] Papamakarios | Normalizing flows for probabilistic modeling and inference[END_REF]. We use implementations from the sbibm package [START_REF] Lueckmann | Benchmarking simulation-based inference[END_REF] to ensure uniform and consistent experimental setups. Samples from the true posterior distributions for both examples are obtained via MCMC and used to compare ℓ-C2ST(NF) to the oracle-C2ST framework. We include results for local-HPD implemented using the code repository of the authors of [START_REF] Zhao | Diagnostics for conditional density models and bayesian inference algorithms[END_REF] with default hyper-parameters and applied to HPD. 5 First, we evaluate the local consistency of the posterior estimators of each task over ten different observations x o with varying N train and fixed N cal = 10 4 . The first column of Figure 2 displays the MSE statistics for the oracle and ℓ-C2ST frameworks. As expected, we observe a decrease in the test statistics as N train increases: more training data usually means better approximations. For Two Moons the statistic of ℓ-C2ST decreases at the same rate as oracle-C2ST, with notably accurate results for ℓ-C2ST-NF. However, in SLCP both ℓ-C2ST statistics decrease much faster than the oracle. This is possibly due to the higher dimension of the observation space in the latter case, which impacts the training-procedure of ℓ-C2ST on the joint distribution.

We proceed with an empirical analysis based on 50 random test runs for each validation method and computing their rejection-rates to the nominal significance level of α = 0.05. Column 4 in Figure 2 confirms that the false positive rates (or type I errors) for all tests are controlled at desired level. Column 2 of Figure 2 portrays the true positive rates (TPR), i.e. rejecting H 0 (x o ) when p ̸ = q, of each test as a function of N train . Both ℓ-C2ST strategies decrease with N train as in Column 1, with higher TPR for ℓ-C2ST-NF in both tasks. The oracle has maximum TPR and rejects the local consistency of all posterior estimates across all observations at least 90% of the time. Note that SLCP is designed to be a difficult model to estimate 6 , meaning that higher values of TPR are expected ( t ̸ → 0 in Column 1). In Two Moons, the decreasing rejection rate can be seen as normal as it reflects the convergence of the posterior approximator ( t → 0 in Column 1): as N train increases, the task of differentiating the estimator from the true posterior becomes increasingly difficult.

We fix N train = 10 3 (which yields inconsistent q ϕ in both examples) and investigate in Column 3 of Figure 2 how many calibration samples are needed to get a maximal TPR in each validation method. SLCP is expected to be an easy classification task, since the posterior estimator is very far from the true posterior (large values of t in Column 1). We observe similar performance for ℓ-C2ST-NF and ℓ-C2ST, with slightly faster convergence for the latter. Both methods perform better than local-HPD, that never reaches maximum TPR. Two Moons represents a harder discrimination task, as q ϕ is already pretty close to the reference posterior (see Column 1). Here, ℓ-C2ST-NF attains maximum power at N cal = 2000 and outperforms all other methods. Surprisingly, the regression-based oracle-C2ST performs comparably to local-HPD, converging to TPR = 1 at N cal = 5000.

Jansen-Rit Neural Mass Model (JRNMM)

We increase the complexity of our examples and consider the well known Jansen & Rit neural mass model (JRNMM) [START_REF] Jansen | Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns[END_REF]. This is a neuroscience model which takes parameters θ = (C, µ, σ, g) ∈ R 4 as input and generates time series x ∈ R 1024 with properties similar to brain signals obtained in neurophysiology. Each parameter has a physiologically meaningful interpretation, but they are not relevant for the purposes of this section; the interested reader is referred to [START_REF] Luiz | HNPE: Leveraging Global Parameters for Neural Posterior Estimation[END_REF] for more details.

The approximation q ϕ of the model's posterior distribution is a conditioned masked autoregressive flow (MAF) [START_REF] Papamakarios | Masked autoregressive flow for density estimation[END_REF] with 10 layers. We follow the same experimental setting from [START_REF] Luiz | HNPE: Leveraging Global Parameters for Neural Posterior Estimation[END_REF], with a uniform prior over physiologically-relevant parameter values and a simulated dataset from the joint distribution including N train = 50 000 training samples for the posterior estimator and N cal = 10 000 samples to compute the validation diagnostics. An evaluation set of size N v0 = 10 000 is used for ℓ-C2ST-NF.

We first investigate the global consistency of our approximation, which informs whether q ϕ is consistent (or not) on average with the model's true posterior distribution. We use standard tools for this task such as simulation-based calibration (SBC) [START_REF] Talts | Validating bayesian inference algorithms with simulation-based calibration[END_REF] and coverage tests based on highest predictive density (HPD) [START_REF] Zhao | Diagnostics for conditional density models and bayesian inference algorithms[END_REF]. The results are shown in the left panel of Figure 3. We observe that the empirical cdf of the global HPD rank-statistic deviates from the identity function (black dashed line), indicating that the approximator presents signs of global inconsistency. We also note that the marginal SBC-ranks are unable to detect any inconsistencies in q ϕ .

We use ℓ-C2ST-NF to study the local consistency of q ϕ on a set of nine observations x (i) o defined as7

x (i) o = JRNMM(θ (i) o ) with θ (i) o = (135, 220, 2000, g (i) o ) and g (i) o ∈ [-20, +20] . (19)
The right panel of Figure 3 shows that the test statistics of ℓ-C2ST-NF vary in a U-shape, attaining higher values as g o deviates from zero and at the borders of the prior. The plot is overlaid with the 95% confidence region, illustrating how much the test statistics deviate from the local null hypothesis. We demonstrate the interpretability of the results for ℓ-C2ST-NF with a focus on the behavior of q ϕ when conditioned on an observation for which g o = 10. The local PP-plot in Figure 4 summarises the test result: the predicted class probabilities deviate from 0.5, outside of the 95%-CR, thus rejecting the null hypothesis of local consistency at g o = 10. The rest of Figure 4 displays 1D and 2D histograms of samples Θ q ∼ q ϕ (θ | x o ) within the prior region, obtained by applying the learned T ϕ to samples Z ∼ N (0, I 4 ). The color of each histogram bin is mapped to the intensity of the corresponding predicted probability in ℓ-C2ST-NF and informs the regions where the classifier is more (resp. less) confident about its choice of predicting class 0.8 This relates to regions in the parameter space where the posterior approximation has too much (resp. not enough mass) w.r.t. to the true posterior: q ϕ > p (resp. q ϕ < p). We observe that the ground-truth parameters are often outside of the red regions, indicating positive bias for µ and σ and negative bias for g in the 1D marginal. It also shows that the posterior is over-dispersed in all 2D marginals. See Appendix D for results on all observations x (i) o .

Discussion

We have presented ℓ-C2ST, an extension to the C2ST framework tailored for SBI settings which requires only samples from the joint distribution and is amortized along conditioning observations. Strikingly, empirical results show that, while ℓ-C2ST does not have access to samples from the true posterior distribution, it is actually competitive with the oracle-C2ST approach that does. This comes at the price of training a binary classifier on a potentially large dataset to ensure the correct calibration of the predicted probabilities. Should this be not the case, some additional calibration step for the classifier can be considered [START_REF] Dheur | A large-scale study of probabilistic calibration in neural network regression[END_REF].

Notably, ℓ-C2ST allows for a local analysis of the consistency of posterior approximations and is more sensible, precise, and computationally efficient than its concurrent method, local-HPD. Appendix F.4 provides a detailed discussion of these statements, based on results obtained for additional benchmark examples. When exploiting properties of normalizing flows, ℓ-C2ST can be further improved as demonstrated by encouraging results on difficult posterior estimation tasks (see Table 2 in Appendix F). We further analyze the benefits of this -NF version in Appendix F.3. ℓ-C2ST provides necessary, and sufficient, conditions for posterior consistency, features that are not shared by other standard methods in the literature (e.g. SBC). When applied to a widely used model from computational neuroscience, the local diagnostic proposed by ℓ-C2ST offered interesting and useful insights on the failure modes of the SBI approach (e.g. poor estimates on the border of the prior intervals), hence demonstrating its potential practical relevance for works leveraging simulators for scientific discoveries.

Limitations and Perspectives

Training a classifier with finite data. The proposed validation method leverages classifiers to learn global and local data structures and shows great potential in diagnosing conditional density estimators. However, it's validity is only theoretically guaranteed by the optimality of the classifier when N v → ∞. In practice, this can never perfectly be ensured. Figure 6 in Appendix F.2 shows that depending on the dataset, ℓ-C2ST can be more or less accurate w.r.t. the true C2ST. Therefore, one should always be concerned about false conclusions due to a far-from-optimal classifier and make sure that the classifier is "good enough" before using it as a diagnostic tool, e.g. via cross-validation. Note, however, that the MSE test statistic for ℓ-C2ST is defined by the predicted class probabilities and not the accuracy of the classifier, thus one should also check how well the classifier is calibrated.

Why Binary Classification? An alternative to binary classification would have been to train a second posterior estimate in order to assess the consistency of the first one. Indeed, one could ask whether training a classifier is inherently easier than obtaining a good variational posterior, the response to which is non-trivial. Nevertheless, we believe that adding diversity into the validation pipeline with two different ML approaches might be preferable. Furthermore, building our method on the C2ST framework was mainly motivated by the popularity and robustness of binary classification: it is easy to understand and has a much richer and stable literature than deep generative models. As such, we believe that choosing a validation based on a binary classifier has the potential of attracting the interest of scientists across various fields, rather than solely appealing to the probabilistic machine learning community.

Possible extensions and improvements. Future work could focus on leveraging additional information of q while training the classifier as in [START_REF] Yao | Discriminative calibration[END_REF]. For example, by using more samples from the posterior estimator q or its explicit likelihood function (which is accessible when q is a normalizing flow). On a different note, split-sample conformal inference could be used to speed up the p-value calculations (avoiding the time-consuming permutation step in Algorithm 1).

In summary, our article shows that ℓ-C2ST is theoretically valid and works on several datasets, sometimes even outperforming local-HPD, which to our knowlegde is the only other existing local diagnostic. Despite facing some difficulties for certain examples (just like for other methods as well), an important feature of ℓ-C2ST is that one can directly leverage from improvements in binary classification to adapt and enhance it for any given dataset and task. This makes ℓ-C2ST a competitive alternative to other validation approches, with great potential of becoming the go-to validation diagnostic for SBI practitioners. The integral of p -q over {p = q} is zero, which implies that q<p p(θ) -q(θ) dθ = B=q≤p p(θ) -q(θ) dθ > 0 meaning that the region {θ ∈ S, p(θ) < q(θ)} is is of strictly positive measure, which concludes the proof.

A.2 Proof of Theorem 2

Theorem 2 (Local consistency and regression). If f xo is Bayes optimal and N v → ∞, then tMSE (f xo ) = 0 is a necessary and sufficient condition for the local consistency of q at x o .

Proof. Let d xo be an estimator of

P(C = 1 | Θ; x o ) defined on S xo such that f xo = I dx o > 1 2 .
As this region contains all the data points

Θ q n ∼ q(θ | x o ) (C n = 0) and Θ p n ∼ p(θ | x o ) (C n = 1), the mean squared error tMSE (f xo ) over the dataset D = {(Θ n , C n )} 2N
n=1 is well defined [START_REF] Dey | Calibrated predictive distributions via diagnostics for conditional coverage[END_REF] and

tMSE (f xo ) -→ Nv→∞ t MSE (f xo ) = 1 2 d xo (θ) - 1 2 2 q(θ | x o ) + p(θ | x o ) dθ .
This integral is zero if, and only if d xo (θ) -1 2 2 = 0 for every θ ∈ S xo (all terms are non-negative and q(θ

| x o ) + p(θ | x o ) > 0) 9 , which is equivalent to d xo (θ) = 1 2 for every θ ∈ S xo . Assuming f xo = f ⋆ xo is optimal, we have that d xo (θ) = d * xo (θ) = P(C = 1 | θ;
x o ) and we conclude the proof using the result from equation ( 5) (and knowing that outside of S xo , p = q = 0):

t MSE (f ⋆ xo ) = 0 ⇐⇒ d ⋆ xo (θ) = P(C = 1 | θ; x o ) = 1 2 , ∀θ ∈ S xo ⇐⇒ (5)
H 0 (x o ) holds .

A.3 Proof of Theorem 3

Theorem 3 (Local consistency and single class evaluation). If f is Bayes optimal and N v → ∞, then tMSE0 (f, x o ) = 0 is a necessary and sufficient condition for the local consistency of q at x o .

Proof. Let d be an estimator of P(C = 1 | Θ; X) and f = I d> 1 2 the associated classifier, both defined on S = {w = (θ, x) ∈ R m × R d , q(θ, x) + p(θ, x) > 0} . Suppose that f = f ⋆ is optimal and let x o ∈ R d be a fixed observation. As explained in section 3, we have that

d ⋆ (θ, x o ) = P(C = 1 | θ; x o ) = d ⋆ xo (θ) and f ⋆ (θ, x o ) = f ⋆ xo (θ), ∀θ ∈ S xo . Consider the support S q,xo = {θ ∈ R m , q(θ | x o ) > 0} ⊂ S xo containing all data points Θ q n ∼ q(θ | x o ) from our single-class validation set D v0 = {(Θ q n , 0)} Nv0 n=1
. Therefore our single-class test statistic tMSE0 is well defined [START_REF] Durkan | nflows: normalizing flows in PyTorch[END_REF] and

tMSE0 (f ⋆ , x o ) -→ Nv0→∞ t MSE0 (f ⋆ , x o ) = d ⋆ xo (θ) - 1 2 2 q(θ | x o )dθ
With the same arguments as in the proof of Theorem 2 in A.2, we get that

t MSE0 (f ⋆ , x o ) = 0 ⇐⇒ d ⋆ xo (θ) = P(C = 1 | θ; x o ) = 1 2 , ∀θ ∈ S q,xo ⇐⇒ H 0 (x o ) holds ,
where the second equivalence is true because S q,xo = S xo for p(.

| x o ) = q(. | x o ). Therefore, t MSE0 (f ⋆ , x o ) = 0 is a necessary and sufficient condition for H 0 (x o ).
N.B. Single-class accuracy (Acc 0 ) does not provide a sufficient condition for H 0 (x o ).

Proof. Following the proof of Theorem 1 in A.1, we have that

H 0 (x o ) holds ⇐⇒ Acc(f ⋆ xo ) = 1 2 ⇐⇒ P(f ⋆ xo (Θ q ) = 0) + P(f ⋆ xo (Θ p ) = 1) = 1 . This means that Acc 0 (f ⋆ , x o ) = P(f ⋆ xo (Θ q ) = 0) = 1 2 can only be a sufficient condition for H 0 (x o ) if P(f ⋆ xo (Θ q ) = 0) = P(f ⋆ xo (Θ p ) = 1
), which is not generally true. In conclusion, evaluating Acc 0 (f ⋆ , x o ) only provides a necessary condition for the local null hypothesis (see (⇒) in A.1).

B Algorithms B.1 Generating PP-plots with ℓ-C2ST

Algorithm 5: ℓ-C2ST -local PP-plots for any x o Input: evaluation data set D v0 ; an observation x o ; estimate d of the class probabilities; estimates {d 1 , . . . , d N H } under the null; grid G of PP-levels in (0,1); significance level α Output: empirical CDF-values { F (l; x o )} l∈G of predicted class-0 probabilities;

(

1 -α)-confidence bands {L l (x o ), U l (x o )} l∈G Predict class-0 probabilities {d 0 (v, x 0 ) = 1 -d(v, x o )} v∈Dv0 /* d is an estimate of class 1 */ for l in G do /* Compute empirical CDFs at l */ F (l; x o ) ← 1 Nv0 I d0(v,x0)≤l for h = 1, . . . , N H do Fh (l; x o ) ← 1 Nv0 I d0(v,x0)≤l /* Compute confidence bands at l */ L l (x o ), U l (x o ) ← q α 2 ({ Fh (l; x o )} N H h=1 ), q 1-α 2 ({ Fh (l; x o )} N H h=1 ) /* quantiles */ return { F (l; x o )} l∈G , L l (x o ), U l (x o )} l∈G

C Run-times for each validation procedure

To complete the benchmark experiments from Section 4.1, we analyze the run-times of each validation method to compute the test statistic 10 for an NPE of the SLCP-task, obtained for different values of N train . Results are displayed in Table 1 and computed on average over all given observations x o and for N cal -values that ensure high test power. We here focus solely on the SLCP-task, as the higher dimensional observation space allows to illustrate the differences between local methods (trained on the joint data space, (θ, x) ∈ R 5+8 ) and the oracle (trained on the parameter space only, θ ∈ R 5 ).

As expected, the computation time increases with the sample size N cal (left vs. right part of Table 1). While the oracle C2ST has close to constant run-time for fixed N cal , local methods become faster with increasing N train . We observe that ℓ-C2ST(-NF) has comparable run-times to the oracle C2ST: the amortization cost is negligible, in particular for difficult tasks involving "good" posterior estimators (high N train ). However, this is not the case for local-HPD, which is the most expensive method. Indeed, it involves (1) the costly computation of the HPD statistics on the joint and (2) the training of several (default is n c = 11) classifiers, both of which increase with the sample size and dimension of the data space.

N cal =

D Graphical diagnostics (ℓ-C2ST-NF) for the JRNMM posterior estimator

The following figures present additional results for the interpretability of ℓ-C2ST applied to the JRNMM neural posterior estimator. They complete

E On the cross-entropy loss for ℓ-C2ST

This section aims at facilitating the understanding of ℓ-C2ST by proving the result of equation [START_REF] Durkan | Neural spline flows[END_REF].

As detailed in section 3, the first step in ℓ-C2ST consists in training a classifier to distinguish between N cal data points (Θ n , X n ) and (Θ q n , X n ) from the joint distributions p(θ, x) and q(θ, x) respectively. Here, the same conditioning observations {X n } N cal n=1 are used to construct the data samples for each class (see Algorithm 1). We show that this does not affect the objective function and convergence of the classifier.

The theoretical cross-entropy loss function to distinguish between data (Θ, X) from class C = 1 and class C = 0 is defined by

l CE (d) := - 1 2 E (Θ,X)|C=1 [log (d(Θ, X))] - 1 2 E (Θ,X)|C=0 [log (1 -d(Θ, X))] . (20)
Note that the we have equal marginals X | (C = 1) ∼ p(x) = q(x) = X | (C = 0). 11 This allows us to take the expectation over X and approximate [START_REF] Hashemi | Simulation-based inference for whole-brain network modeling of epilepsy using deep neural density estimators[END_REF] via Monte-Carlo for only one set of conditioning observations {X n } N cal n=1 , and with data points Θ n and Θ q n respectively associated to class C = 1 and C = 0 for a given X n :

l CE (d) = E X - 1 2 E Θ|X,C=1 [log (d(Θ, X))] - 1 2 E Θ|X,C=0 [log (1 -d(Θ, X))] ≈ - 1 2N cal N cal n=1 log (d(Θ n , X n )) + log (1 -d(Θ q n , X n )) . (21)
Therefore, we can train a classifier to minimize (20) using data from the joint distributions with same conditioning observations. The obtained estimate d = arg min{l CE (d)} of the class probabilities is defined for every

(θ, x) ∈ R m × R d by d(θ, x) ≈ p(θ,x)
p(θ,x)+q(θ,x) . As p(x) = q(x), we recover the class probabilities of the optimal Bayes classifier on the conditional data space for any given x ∈ R d :

d ⋆ (θ, x) = p(θ | x) p(θ | x) + q(θ | x o ) = d ⋆ x (θ) . (22)
For an example, see works related to neural ratio estimation (NRE) [START_REF] Hermans | Likelihood-free MCMC with amortized approximate ratio estimators[END_REF][START_REF] Delaunoy | Towards reliable simulation-based inference with balanced neural ratio estimation[END_REF]: these algorithms implicitly use a classifier to distinguish between the joint and marginal distributions p(θ, x) and p(θ)p(x). Like in our case, both classes are modeled using the same observations X n obtained via the simulator.

F Additional Experiments on several benchmark examples

The results on the two benchmark examples in Figure 2 give first intuitions about the validity and behaviour of ℓ-C2ST(-NF), our proposal, w.r.t. the oracle C2ST and the alternative local-HPD methods. In this section, Figure 5 extends Figure 2 with results on additional benchmark examples and more detailed results on the correlation between the test statistics of ℓ-C2ST(-NF) and the oracle C2ST for different conditioning observations are shown in Figure 6 and Table 3. We refer the reader to Table 2 for a description of all benchmark examples (data dimensionality, posterior structure, challenges) and a summary of the main results comparing local methods.

While investigating the scalability of the algorithms to high dimensional data spaces, these additional experiments help to further analyze how well ℓ-C2ST(-NF) captures the true C2ST, when it outperforms local-HPD, and better understand the benefit of the -NF version. These points are further detailed in the following subsections. The goal is to create a first guideline for when our proposal can and should be used, while raising awareness of its limitiations (i.e. when it should not be used or at least be adpated for improved performance).

F.1 Scalability to high dimensions

First of all, note that in the specific case of SBI, the dimension of the parameter space is typically of order 10 0 to 10 1 and m ≈ 10 2 is already often considered as high dimensional. The observation space, however, can be of higher dimension (e.g. d ≈ 10 3 for time-series), but summary statistics are often used to reduce the dimension of the observations to the order of d ≈ 10 Gaussian Linear Uniform tasks were not included in this analysis, as here, the difficulty of the classification task has more impact on statistical power than the data dimensionality 12 .

Interestingly, we observe in the Bernoulli GLM Raw task, that ℓ-C2ST-NF scales well to the highdimensional observation space (faster convergence maximum TPR compared to Bernoulli GLM), while the normal ℓ-C2ST and local-HPD significantly lose in statistical power. It should also be noted that local-HPD performs significantly worse in medium dimensions (cf. Bernoulli GLM or even SLCP) than in low dimensions (Gaussian Mixture and Two Moons), though this could be because of the complex posterior structure.

F.2 Accuracy of ℓ-C2ST(-NF) w.r.t. the true C2ST Figures 2 and 5 compare our method to the oracle C2ST, but only in terms of statistical power, as the local analysis is limited to the averaged results over 10 different reference observations x o . We here provide a more detailed local analysis that examines how the ℓ-C2ST(-NF) test statistics correlate with those from the oracle C2ST, by plotting them against each other. The results obtained for the above mentioned 10 reference observations are shown in Figure 6a. To allow for more robust conclusions, we also show results obtained for a total of 100 different reference observations (cf. Figure 6b), as well as quantitative results on the statistical significance of the correlation indices in Table 3.

Overall, the scattered points are not too far from the diagonal, which indicates that there is some correlation between the test statistics for ℓ-C2ST(-NF) and those from the oracle C2ST. This correlation becomes weaker when N train becomes larger, since the test statistics in these cases tend to zero and can start to be confused with noise. This observation is consistent with the results in Table 3, showing the p-values of the Pearson test, a standard tests for the statistical significance of the correlation indices between the scores.

Another general trend is that the scattered points tend to be below the diagonal, indicating that ℓ-C2ST(-NF) is less sensible to local inconsistencies than the oracle C2ST. This behaviour was expected, as ℓ-C2ST is trained on the joint and thus less precise. Interestingly this doesn't apply to

F.3 Benefit of the -NF version

The results of all benchmark experiments indicate that the -NF version of ℓ-C2ST works better when the (true) posterior distribution of the model is "more complicated" than a Gaussian distribution (see Table 2 at the beginning of Appendix F). This is for example the case for the Two Moons and SLCP tasks: the posterior distributions are globally multi-modal and locally structured. We observe in Column 3 of Figure 2 in the main paper that the -NF version requires less samples (i.e. lower N cal ) to reach maximum power/TPR. This is also the case for the additional Bernoulli GLM task (see Column 3 of Figure 5 in Appendix F.1). In contrast, for Gaussian Mixture and Gaussian Linear Uniform, where the posterior is Gaussian, the normal ℓ-C2ST is as powerful or even better than its -NF counterpart.

Finally, we refer the reader to the analysis in Section F.1 to point out an interesting observation: for the Bernoulli GLM, the -NF version scales much better to high dimensional observation spaces than the normal ℓ-C2ST. Note that this does not allow us to make any general conclusions, but it might be worth further investigating this result.

F.4 Advantages of ℓ-C2ST w.r.t. local-HPD First of all, it is important to mention that having uniform HPD-values is not a sufficient condition for asserting the null hypothesis of consistency (see end of Section 3.3 in [START_REF] Zhao | Diagnostics for conditional density models and bayesian inference algorithms[END_REF]). This is a clear disadvantage compared to our proposal, which provides a necessary and sufficient proxy for inspecting local posterior consistency. Furthermore, the HPD methodology summarizes the whole information concerning θ into a single scalar, while in ℓ-C2ST we handle the θ-vector in its multivariate form. In medium-high θ-dimensions (m > 2 as in Bernoulli GLM) or for complex posterior distributions (SLCP), such summarized information might discard too much information and not be enough to satisfactorily assess the consistency of the posterior estimator. Indeed, the tasks where local-HPD has similar statistical power to ℓ-C2ST are either in low dimensions (Two Moons) and / or have a Gaussian posterior (Gaussian Mixture / Gaussian Linear Uniform). See Table 2 for a summary of those results w.r.t. data dimensionality and posterior structure. For a detailed analysis of the scalability to high dimensional data spaces see Appendix F.1.

Finally, local-HPD in its naive implementation is much less efficient than ℓ-C2ST (see Appendix C). Note however that a new "amortized" version of local-HPD has recently been proposed [START_REF] Dey | Calibrated predictive distributions via diagnostics for conditional coverage[END_REF] and could be interesting to look at.

Algorithm 1 :

 1 ℓ-C2ST -training the classifier on data from the joint distribution Input: posterior estimator q; calibration data D cal = {Θ n , X n } N cal n=1 ; classifier f ; number of samples N H from the distribution under the null hypothesis Output: estimate d of the class probabilities; estimates {d 1 , . . . , d N H } under the null hypothesis /* Construct classification training set *

Figure 2 :

 2 Figure2: Results on two examples from the SBI benchmark: SLCP and Two Moons. We compare ℓ-C2ST and ℓ-C2ST-NF (dashed) to the oracle C2ST and local-HPD. Columns 1 and 2 display the test statistic and empirical power as a function of N train , while Columns 3 and 4 show the empirical power and type I error for varying N cal . The ℓ-C2ST(-NF) statistics are comparable to the oracle, as their decreasing behaviour reflects the convergence of NPE to the true posterior for large training datasets. We also note that ℓ-C2ST-NF is uniformly better than ℓ-C2ST (i.e. higher power for all N train and increases faster with N cal ), and both reach maximum statistical power with smaller N cal than local-HPD. All Type I errors are controlled at α = 0.05. Experiments were performed over 10 different observations x o (mean and std) and Columns 2-4 used additional 50 random test runs.

Figure 3 :

 3 Figure 3: Results for global and local tests on JRNMM. Left: PP-plots for the marginal SBC and global HPD rank statistics. Right: Test statistics for ℓ-C2ST-NF on observations with varying g o . SBC fails to detect any inconsistency of q ϕ , while HPD only provides a global assessment, unlike ℓ-C2ST which locally explains the inconsistencies in q ϕ .

Figure 4 :

 4 Figure 4: Graphical diagnostics of ℓ-C2ST-NF for JRNMM. Top right panel displays the empirical CDF of the classifier (blue) overlaid with the theoretical CDF of the null hypothesis (step function at 0.5) and 95% confidence region of estimated classifiers under the null displayed in gray. The pairplot displays histograms of samples from q ϕ within the prior region and dashed lines indicate values of θ o used to generate the conditioning observation x o . The predicted probabilities are mapped on the colors of the bins in the histogram. Blue-green (resp. orange-red) regions indicate low (resp. high) predicted probabilities of the classifier. Yellow regions correspond to chance level, thus q ϕ ≈ p.

AppendicesA 16 A. 1 16 A. 2 17 A. 3 17 B 19 B. 1 19 C 20 D 21 E 23 F 24 F. 1 24 F. 2 25 F. 3 27 F. 4

 1611621731719119202123241242253274 Proofs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Proof of Theorem 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Proof of Theorem 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Proof of Theorem 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Generating PP-plots with ℓ-C2ST . . . . . . . . . . . . . . . . . . . . . . . . . . . Run-times for each validation procedure . . . . . . . . . . . . . . . . . . . . . . . . Graphical diagnostics (ℓ-C2ST-NF) for the JRNMM posterior estimator . . . . . . . . On the cross-entropy loss for ℓ-C2ST . . . . . . . . . . . . . . . . . . . . . . . . . . . Additional Experiments on several benchmark examples . . . . . . . . . . . . . . . Scalability to high dimensions . . . . . . . . . . . . . . . . . . . . . . . . . . . . Accuracy of ℓ-C2ST(-NF) w.r.t. the true C2ST . . . . . . . . . . . . . . . . . . . . Benefit of the -NF version . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Advantages of ℓ-C2ST w.r.t. local-HPD . . . . . . . . . . . . . . . . . . . . . . . 27

Figure 4 .

 4 Results are shown for all given observations associated to ground-truth gain values g o = -20, -15, -10, -5, 0, 5, 10, 15, 20.In each Figure, the top right panel displays the empirical CDF of the classifier (blue) overlaid with the theoretical CDF of the null hypothesis (step function at 0.5) and 95% confidence region of estimated classifiers under the null displayed in gray. The pairplot displays histograms of samples from the posterior estimator q ϕ within the prior region and dashed lines indicate values of θ o used to generate the conditioning observation x o . The predicted probabilities are mapped on the colors of the bins in the histogram. Blue-green (resp. orange-red) regions indicate low (resp. high) predicted probabilities of the classifier. Yellow regions correspond to chance level, thus q ϕ ≈ p.

1 .

 1 In Section 4 we analyze the results obtained for rather low-dimensional benchmark examples (Two-Moons:m = 2, d = 2, SLCP: m = 5, d = 8). As an extension of Figure 2, Figure 5 includes additional benchmark examples with low and medium dimensionality: Gaussian Mixture (m = 2, d = 2) Gaussian Linear Uniform (m = 10, d = 10) and Bernoulli GLM (m = 10, d = 10). Furthermore, the Bernoulli GLM Raw task allows us to analyze how our method scales to high dimensional observation spaces only (without parameter-space / task variability): it considers raw observation data with d = 100, as opposed to sufficient summary statistics in the Bernoulli GLM task.Column 3 in Figure5shows that ℓ-C2ST requires more data to converge to the oracle C2ST (at maximum power TPR = 1) as the data dimensionality increases: N cal ≈ 2000 for Two Moons and SLCP, but N cal ≈ 5000 for the Bernouilli GLM task. Note that the Gaussian Mixture and

Figure 5 :

 5 Figure 5: Results on additional sbibm benchmark examples using the same experimental setup as for Figure 2 in Section 4 (50 test runs, mean / std over 10 different reference observations).

  Algorithm 4: ℓ-C2ST-NF -precompute the null distribution for any estimatorInput: calibration data D cal = {Θ n , X n } N cal n=1 ;classifier f ; number of null samples N H Output: estimates {d 1 , . . . , d N H } of the class probabilities under the null for h in 1, . . . N H do /* Construct classification training set */ Sample

	Get estimate d of the class probabilities	*/
	Train the classifier f on D
	d ← f probability
	return d

  Get estimate d of the class probabilities */ Train the classifier f on D d h ← f probability return {d 1 , . . . , d N H }

Table 1 :

 1 Run-time (in seconds) to compute the test statistic for the SLCP task (mean over observations). C2ST has close to constant run-time for fixed N cal . Local methods become faster with increasing N train and ℓ-C2ST(-NF) stays comparable to the oracle C2ST, even for N cal = 10 000. While the amortization cost of ℓ-C2ST is not an issue, local-HPD is always at least 30 times slower.

				5 000			N cal = 10 000
	N train	10 2	10 3	10 4	10 5	10 2	10 3	10 4	10 5
	oracle C2ST	5.47	4.52	5.95	7.56	16.36	18.03	23.3	15.33
	ℓ-C2ST	5.92	5.09	1.78	1.81	27.62	34.06	17.9	3.65
	ℓ-C2ST-NF	6.98	6.84	6.38	1.72	43.99	25.01	18.56	7.62
	local-HPD	282.19 282.85 279.38 282.5 956.91 938.21 682.92 530.45

Table 2 :

 2 Description of benchmark examples and summary of main results for local methods.

		Dimension	Posterior	Challenge	Better local
		(θ, x)	structure		method
	SLCP	low	4 symmetrical complex posterior	ℓ-C2ST-NF
		(5, 8)	modes
	Two Moons	low	bi-modal,	global and local	local-HPD /
		(2, 2)	crescent shape	structure	ℓ-C2ST-NF
	Gaussian Mixture	low	2D Gaussian	large vs. small	all similar
		(2, 2)		s.t.d. in GMM
	Gaussian Linear	medium	multivariate	dimensionality	local-HPD /
	Uniform	(10, 10)	Gaussian	scaling	ℓ-C2ST
	Bernoulli GLM	medium	unimodal,	dimensionality	ℓ-C2ST(-NF)
		(10, 10)	concave	scaling
	Bernoulli GLM Raw	high	unimodal,	raw observations	ℓ-C2ST-NF
		(10, 100)	concave	(no summary stats)

i.e. the approximation q(θ | x) is close to p(θ | x) on average for all values of x ∈ R d , so we can quickly generate samples from the posterior for any choice of conditioning observation without redoing any training.

i.e. it is the classifier with lowest possible classification error for the dataset D.

Code is available at https://github.com/JuliaLinhart/lc2st.

The average run-times for each validation method are provided in Appendix C.

SLCP = simple likelihood, complex posterior

Note that the uniform prior for g is defined in [-20, +20] when training q ϕ with NPE.

Specifically, we compute the average predicted probability of class 0 for data points Z ∼ N (0, I4) corresponding to samples T ϕ (Z; xo) ∼ q ϕ (θ | xo) within each histogram bin.

Note that this integral can also be zero if Sx o is of measure zero. But as mentioned at the beginning of this appendix, this has generally no practical implications.

Note that in order to compute to compute p-values, we need to compute the test statistic NH times under the null hypothesis. The number of classifiers we need to train depends on how many we need to compute the test statistic (NH for ℓ-C2ST vs. NH × nc for local-HPD). In summary, if ℓ-C2ST is more efficient in computing a single test statistic, it will also be more efficient to compute NH test statistics.

The joint distributions p(θ, x) and q(θ, x) are both modeled using samples X ∼ p(x | Θ) obtained from the prior Θ ∼ p(θ), which implies that the marginals p(x) and q(x) are both defined by p(x | θ)p(θ)dθ.

Local methods a TPR < 1 at Ntrain = 1000 (seeColumn 2), which means that the classification task is harder and requires more data: local methods never reach maximum TPR = 1, and even the oracle C2ST (with MSE test statistic) needs N cal = 2000, at least four times more than for the other tasks.

Acknowledgments

Julia Linhart is recipient of the Pierre-Aguilar Scholarship and thankful for the funding of the Capital Fund Management (CFM). Alexandre Gramfort was supported by the ANR BrAIN (ANR-20-CHIA0016) grant while in his role at Université Paris-Saclay, Inria.

A Proofs

In what follows, we assume that it is sufficient for the null hypothesis H 0 to hold on any set C ⊆ R m of strictly positive measure, rather than requiring it to hold for all points θ ∈ R m . In practice, it generally has no practical implications if the posterior estimator is inconsistent with the true posterior (q ̸ = p) on a set of measure zero, since those sets don't have any real statistical significance.

For the proofs of Theorems 1 and 2, we will consider a classifier f xo defined for a fixed observation

A.1 Proof of Theorem 1 Theorem 1 (Local consistency and classification accuracy). If f xo is Bayes optimal and N v → ∞, then tAcc (f xo ) = 1/2 is a necessary and sufficient condition for the local consistency of q at x o .

Proof. As S xo contains all data points

n=1 is well defined [START_REF] Dax | Group equivariant neural posterior estimation[END_REF] and

Let's show that if f xo is Bayes optimal, then H 0 (x o ) holds ⇐⇒ Acc(f xo ) = 1 2 . (⇒): Suppose H 0 (x o ) holds. The optimality of f xo implies that f xo (Θ q ) and f xo (Θ p ) are Bernoulli random variables B( 12 ) (see interpretation of equation ( 5)), and so Acc(

Suppose that H 0 (x o ) does not hold, there exists a set

of strictly positive measure (w.r.t. p or q, which ever is non zero on that set). We can decompose C into the direct sum of

Either A or B is necessarily of strictly positive measure. Let's say A (see Lemma 1 for the symmetric case).

A is exactly the region where Prob(C = 1 | Θ = θ) > 1 2 and thus where f xo (θ) = 1; B is the region where f xo (θ) = 0. We therefore get that: A is of strictly positive measure), so the integral in the last equality is strictly positive and we get Acc(f x ) > 1 2 , which concludes the proof.

Lemma 1. Let p and q be two probability density functions defined on a space S. If there exists a set A ⊆ S of strictly positive measure such that ∀θ ∈ A, q(θ) < p(θ), then there exists a set B ⊆ S of strictly positive measure such that ∀θ ′ ∈ B, q(θ ′ ) > p(θ ′ ).

Proof. We know that S p = S q = 1 or equivalently, using

By grouping the integrals over A on one side and the ones over B on the other, we get:

which is non-negative because A is assumed to be of strictly positive measure and q -p > 0 everywhere in A.

the Gaussian Mixture task. This could be due to a big variability in the local consistency of q: while trained on the joint ℓ-C2ST(-NF) could overfit on the "bad" observations, resulting in higher test statistics for observations where the true C2ST statistic would be small.

Finally, across all benchmark examples we observe results that are consistent with the ones from Figure 5: higher correlation means higher statistical power.


Zaineb Amor

Pierre-Antoine Comby

Caroline Le Ster

Alexandre Vignaud

Philippe Ciuciu

NON-CARTESIAN NON-FOURIER FMRI IMAGING FOR HIGH-RESOLUTION RETINOTOPIC MAPPING AT 7 TESLA

Keywords: fMRI, non-Cartesian, Compressed Sensing, ultra-high field MRI, non-Fourier imaging, 7T

Functional MRI (fMRI) is an invaluable tool for neuroscience, but it necessitates a trade-off between spatial and temporal resolution to maintain a reasonable temporal signal-to-noise ratio (tSNR). Non-Cartesian acquisition schemes are more efficient sampling strategies as compared to Cartesian ones and were proposed as a way to push the limits of spatiotemporal resolution further. 3D-SPARKLING is a novel non-Cartesian scheme recently evaluated for fMRI applications. However, 3D-SPARKLING is highly sensitive to B0 field imperfections. Such imperfections are detrimental to fMRI applications, especially at ultra-high magnetic fields. In this work, we collect measurements of the static and dynamic B0 field perturbations concurrently with 3D-SPARKLING fMRI data acquisition and retrospectively correct these perturbations during image reconstruction. The advantages of this image reconstruction strategy are assessed on the statistical sensitivity to the BOLD contrast during a retinotopic mapping fMRI experiment. Importantly, a gain of 44% (resp., 159%) additionally activated voxels was quantified when adopting this brain activity-enhanced image reconstruction technique at a type-I statistical control level of 0.001 without multiple comparisons correction (resp., 0.05 with false discovery rate correction). Additionally, significantly improved retinotopic maps were retrieved on the cortical surface.

INTRODUCTION

Functional MRI is today one of the most widespread neuroimaging techniques thanks to its non-invasiveness and high spatial resolution compared to other modalities such as electroencephalography (EEG) or magnetoencephalography (MEG). In contrast, fMRI data has a low temporal resolution on average (around 1s) as the data acquisition, performed in k-space, is relatively slow at least in the Cartesian imaging setting, where the Nyquist criterion must be met. Accelerated Cartesian acquisition schemes such as 2D-SMS EPI [START_REF] Moeller | Multiband multislice GE-EPI at 7 Tesla, with 16-fold acceleration using partial parallel imaging with application to high spatial and temporal whole-brain fMRI[END_REF] or 3D-EPI [START_REF] Poser | Three dimensional echo-planar imaging at 7 Tesla[END_REF], which sample the k-space using readouts that span parallel planes, combined with parallel imaging [START_REF] Klaas P Pruessmann | Encoding and Reconstruction in Parallel MRI[END_REF][START_REF] Deshmane | Parallel MR Imaging[END_REF], enable undersampling without aliasing artifacts in the image domain by exploiting the data redundancy from the multiple receiver channels during volume-wise inverse Fourier transform reconstruction [START_REF] K P Pruessmann | SENSE: Sensitivity encoding for fast MRI[END_REF][START_REF] Griswold | Generalized Autocalibrating Partially Parallel Acquisitions (GRAPPA)[END_REF]. Non-Cartesian fMRI methods such as spiral readouts [START_REF] Kasper | Advances in spiral fMRI: A highresolution study with single-shot acquisition[END_REF] or TURBINE [START_REF] Graedel | Ultrahigh Resolution fMRI at 7T Using Radial-Cartesian TURBINE Sampling[END_REF] have been popularized in recent years due to their ability to sample the k-space more efficiently. 3D-SPARKLING [START_REF] G R Chaithya | Optimizing full 3D SPARKLING trajectories for highresolution Magnetic Resonance imaging[END_REF] is a new encoding scheme that was assessed for fMRI and proved to be competitive with state-of-the-art 3D-EPI to detect evoked brain activity during a visual task [START_REF] Amor | Prospects of non-cartesian 3d-sparkling encoding for functional mri: A preliminary case study for retinotopic mapping[END_REF]. A comparison with 2D-SMS EPI would be more complicated as 2D and 3D acquisitions generate different contrasts.

As an accelerated and massively under-sampled non-Cartesian acquisition method, 3D-SPARKLING requires a nonlinear Compressed Sensing (CS) image reconstruction method involving a nonuniform Fourier operator and sparsity-promoting priors. Proximal gradient methods such as the Fast Iterative Shrinkage-Tresholding Algorithm (FISTA) [START_REF] Beck | A fast iterative shrinkagethresholding algorithm for linear inverse problems[END_REF] or the Promixal Optimized Gradient Method (POGM) [START_REF] Kim | Adaptive restart of the optimized gradient method for convex optimization[END_REF] can be used to solve the associated optimization problem. Furthermore, unlike Cartesian methods where ∆B0 inhomogeneities result mainly in geometric distortions and can be corrected by means of preprocessing steps such as the TOPUP approach [START_REF] Smith | Advances in functional and structural MR image analysis and implementation as FSL[END_REF], the impact of B0 imperfections (static inhomogeneities in space and temporal fluctuations over time) is more complex in 3D-SPARKLING acquisitions and result in a superposition of blurring, signal loss and geometric distortions. To correct for both static and dynamic imperfections, the model of the collected signal in k-space during each sampling trajectory must be extended from the classic (non-uniform) Fourier transform to a pseudo-Fourier operator.

In this work, we retrospectively perform static and up-to-thefirst-order dynamic B0 perturbations correction on resting-state and task-based fMRI data collected during a retinotopic mapping experiment. This correction is achieved by means first of collecting additional external measurements of the field fluctuations using a field Camera [START_REF] De Zanche | NMR Probes for Measuring Magnetic Fields and Field Dynamics in MR Systems[END_REF][START_REF] Barmet | A transmit/receive system for magnetic fieldmonitoring of in vivo mri[END_REF](Skope Clip-on camera) and second through the definition of a non-Fourier forward operator that encodes an extended signal model taking B0 field imperfections [START_REF] B P Sutton | Fast, iterative image reconstruction for mri in the presence of field inhomogeneities[END_REF] into account. A study has been conducted over three different healthy volunteers, and the impact of B0 imperfections corrections is assessed for image quality, tSNR, and sensitivity to the BOLD contrast.

THEORY

Extended signal model

3D-SPARKLING is a segmented k-space encoding scheme, meaning that the k-space data is collected over N shot consecutive readouts (also called shots) separated by different excitation RF pulses. All readouts have the same duration denoted T obs , and each one of them samples the center of the k-space at the echo time denoted T E. Furthermore, in the case of parallel imaging (multiple receive coils), all coils collect simultaneously their own data in the k-space domain convolved with their own sensitivity profile. During the s-th readout, the NMR signal µ s ℓ (t) collected by the ℓ-th coil for each t ∈ [T E -T obs 2 , T E + T obs 2 ] can be modeled by the idealized signal model in Eq. (1):

µ s ℓ (t) = FOV xℓ (r)e -2ıπks(t)•r dr , (1)
where xℓ (r) = ζ ℓ (r)x(r), x(r) and ζ ℓ (r) respectively encode the image and sensitivity map associated with the ℓ-th coil at the spatial position r. Note that ks(t) (in m -1 ) is the prescribed (theoretical) kspace position at time t (in s) relative to the s-th readout. Therefore, a non-uniform Fourier transform binds µ s ℓ (t) to xℓ (r). Eq. (1) supposes a homogeneous B0 and accurate application of the prescribed trajectories (k=[ks] 1≤s≤N shot ) by the MR system. Unfortunately, in the case of highly under-sampled 3D-SPARKLING, strong degradation in the image quality occurs due to B0 imperfections arising from the patient (susceptibility changes at the tissue-air interfaces, physiological noise) and the system (eddy currents, temperature drifts, system instabilities) [START_REF] Vannesjo | Gradient system characterization by impulse response measurements with a dynamic field camera[END_REF][START_REF] Kasper | Monitoring, analysis, and correction of megnetic field fluctuations in Echo Planar Imaging time series[END_REF][START_REF] Bollmann | Analysis and correction of field fluctuations in fMRI data using field monitoring[END_REF]. One way to overcome this issue is to take such imperfections in the signal model into account during MR image reconstruction: Eq. ( 1) then becomes

µ s ℓ (t) = e -2iπt∆B dyn 0,s FOV xℓ (r)e -2ıπ[∆B stat 0 (r)t+ ks(t)•r] dr (2)
where ∆B stat 0 (r) (in Hz) denotes the static inhomogeneities of the B0 field in space and ∆B dyn 0,s (in Hz) and ks = ks +δks (in m -1 ) denote respectively its zeroth order dynamic fluctuation and the measured trajectory (deviated from the prescribed one due to first order fluctuation δks). ∆B dyn 0,s is slowly varying and considered constant during a shot. In Eq. ( 2), the term ∆B stat 0 (r)t depends on the image domain making the integral dependent both on the image and k-space domains which is not compatible with the usual Fourier transform model.

Linear approximation of the non-Fourier model

According to Eq. ( 2), the discretized adjoint operator can be written as follows:

xℓ (rn) = t e 2ıπt∆B dyn 0,s µ s ℓ (t)e 2ıπ[∆B stat 0 (r)t+ ks(t)•r] (3)
The mixed term e 2ıπ∆B stat 0 (rm)t = P bm,pcp,n is split in a P -rank linear decomposition using a SVD (as implemented in [START_REF] Daval-Frérot | Iterative static field map estimation for offresonance correction in non-cartesian susceptibility weighted imaging[END_REF] following the methods in [START_REF] B P Sutton | Fast, iterative image reconstruction for mri in the presence of field inhomogeneities[END_REF][START_REF] Fessler | Toeplitz-based iterative image reconstruction for MRI with correction for magnetic field inhomogeneity[END_REF]) and we then recover a linear combination of P (non-uniform) Fourier transforms as follows:

xℓ (rn) = P p=1 cp,n T obs tm=0 bm,p e 2ıπt∆B dyn 0,s µ s ℓ (tm) µ s ℓ (tm)
e 2ıπ ks(tm)•rn . (4)

The higher the number of interpolators P in the sum, the more accurate the approximation, however, at the expense of additional computing time.

Since the term related to ∆B dyn 0,s is outside of the integral in Eq. ( 2), the zeroth order dynamic fluctuations can be corrected by simply demodulating each µ s ℓ (t) by the corresponding e 2iπt∆B dyn 0,s to obtain µ s ℓ (t). As Eq. ( 4) holds for all frequencies ks(tm) and locations rn across all the N shot readouts, we can summarize the perturbed acquisition in Eq. ( 4), as a linear combination of adjoint non-uniform Fourier transforms F Ω , yielding a coil-specific image xℓ from the measured frequencies at locations Ω and associated corrected values ( µ ℓ ):

xℓ = P p=1 cp ⊙ F * Ω (bp ⊙ µ ℓ ) = F * P, Ω ( μℓ ), (5)
where ⊙ denotes the element-wise product.

Sparse reconstruction

In this work, we perform a CS-based reconstruction, where each volume is reconstructed independently from the others using Lchannel k-space data. Therefore, for each complex-valued volume (xj ∈ C N , j = 1, .., Nt), we consider the complex-valued k-space data y j,ℓ = [ µ 1 ℓ , . . . , µ

N shot ℓ
] + n j,ℓ , associated with the set of sampled frequencies points Ωj = Ω = [ k1, . . . , kN shot ] relative to the j-th volume. Future work will address time or scan-varying under-sampling patterns Ωj. Here n j,ℓ is an additive zero-mean white Gaussian noise with variance σ2 ℓ and y j,ℓ ∈ C M . N and M are, respectively, the number of voxels in each volume and that of the k-space measurements collected by each coil. Then the reconstruction problem consists in minimizing the following objective function

xj = arg min x∈C N 1 2 L ℓ=1 σ -2 ℓ F P, Ω S ℓ xj -ỹj,ℓ 2 2 + λg(Ψ(xj)) (6)
where S ℓ denotes the sensitivity profile of the ℓ-th coil and F P, Ω the linearized over P interpolators off-resonance corrected Fourier operator. Ψ is the sparsifying basis (e.g. a tight frame), g is typically chosen as the l1-norm, and λ is a regularization parameter.

Although more advanced low-rank + sparse regularization could be envisaged easily (see [START_REF] Andrii Y Petrov | Improving temporal resolution in fMRI using a 3D spiral acquisition and low rank plus sparse (L+ S) reconstruction[END_REF] for instance), this requires correcting the number of degrees of freedom during subsequent statistical analysis. Due to this, every volume xj for j = 1, . . . , Nt, is reconstructed independently by solving Eq. ( 6). Practically, in our experiments, the sparsifying basis (Ψ) we used was the sym8 wavelet basis and the other parameters were set to L = 32, P = 30, λ = 10 -8 .

MATERIALS AND METHODS

Multi-modal data acquisition

Functional MRI data was collected at 7T (7T Magnetom investigational MRI, Siemens Healthineers, Erlangen, Germany) using a 1Tx-32Rx head coil (Nova Medical, Willmington, CO, USA) from 3 healthy volunteers during resting-state and a retinotopic mapping experiment 1 at a spatial resolution of 1mm3 , a temporal resolution (T R vol ) of 2.4s, T R shot /T E = 50/20ms and a 3D fieldof-view (FOV) of (192,192,128) mm 3 . Concurrently and for each acquired MR volume, 16 NMR probes (a field camera) were used to monitor and record the zeroth order field fluctuations ∆B dyn 0 = [∆B dyn 0,1 , . . . , ∆B dyn 0,N shot ] and measure the trajectories played by the MR system k = [ k1, . . . , kN shot ]. Additionally, an external ∆B stat 0 map as well as external sensitivity maps were acquired using a threeecho 3D gradient recalled echo (3D GRE) sequence.

Functional MRI data reconstruction and preprocessing

The fMRI volumes were reconstructed by solving the minimization problem in Eq. ( 6) using 15 iterations of the POGM algorithm for each volume. This method is implemented in the pysap-mri [START_REF] Gueddari | PySAP-MRI: A python package for mr image reconstruction[END_REF] plugin 2 of the pySAP package [START_REF] Farrens | PySAP: Python Sparse Data Analysis Package for multidisciplinary image processing[END_REF]. Motion correction and co-registration of the functional and the T1weighted anatomical MR images were applied using SPM12 3 . The Fig. 1. Comparison of the (A) mean images and (B) tSNR maps yielded by the resting-state fMRI scans sequence collected in the three volunteers (one raw per subject) reconstructed with (A, right col.) and without (A, left col.) correcting off-resonance artifacts due to B0 inhomogeneities. Fully corrected means that static and up-to-the-first-order dynamic terms were taken into account during image reconstruction as in Eq. ( 6). The orange and green arrows depict examples of recovered signal and better-reconstructed details, respectively. The improved image quality and tSNR maps are systematic across the three volunteers. segmentation of the T1-weighted MRI scan into the cortical surface meshes was performed using FreeSurfer 7 .

Statistical analysis of the task-based fMRI data

Retinotopic mapping fMRI data was analyzed for each participant separately, using a first-level general linear model (GLM) where two paradigm-related regressors (parametric, continuous and sinusoidal) were used to capture the BOLD fluctuations elicited by the stimulus presentation. The design matrix also encompasses 6 motion regressors, a polynomial drift and a baseline. A Fisher test over the taskrelated regressors was used to estimate the global effect of interest after thresholding the F-statistic maps over the entire brain for two different strategies: (i) p < 0.001 without correcting for multiple comparisons (ii) p < 0.05 with false discovery rate (FDR) control. The thresholded F-statistic maps were further used to delineate regions of interest (ROIs) in which the retinotopic mapping was computed. The statistical analysis was implemented using Nilearn [START_REF] Abraham | Machine learning for neuroimaging with scikit-learn[END_REF].

RESULTS

Improved image quality and increased tSNR

Fig. 1-(A) compares, subject-wise, the mean images derived from the resting-state sequence of fMRI volumes reconstructed with and without correcting artifacts due to static and dynamic B0 perturba-Table 1. Gain in % of median tSNR in corrected images relative to the native tSNR (uncorrected images) computed over the brain mask for the three different volunteers. The highest gains (in bold font) are retrieved when static and up-to-the-first-order field terms are jointly corrected.

Gain in % of median tSNR Terms corrected

V#1 V#2 V#3 Average ∆B stat 0 +5 +6 +4 +5 ∆B dyn 0 +20 +11 +22 +18 ∆B dyn 0 & δk +26 +13 +23 +21 ∆B stat 0 & ∆B dyn 0 +28 +18 +29 +25 ∆B stat 0 & ∆B dyn 0 & δk +34 +20 +29 +28
tions. We note that the image quality is systematically improved with the proposed correction across the three volunteers. Indeed, performing B0 perturbations correction improves the overall T * 2 contrast. Additionally, the lost signal is recovered, and anatomical details are sharper as depicted, respectively, with orange and green arrows. Fig. 1-(B) shows the boost in tSNR associated with B0 perturbations correction, notably in the anterior and posterior cortex and, along the edges of the brain suggesting that subtle head movementinduced field fluctuations related to breathing were compensated.

Tab. 1 reports the relative gain in % of median tSNR computed over the brain mask when correcting the different field terms during image reconstruction. Although the increase is systematic across volunteers (V#1-3), the relative gain for V#2 is lower than in other volunteers when dynamic fluctuations are corrected. This is likely because this volunteer has less intense breathing therefore, the native fMRI data is less affected by signal fluctuations. Otherwise, the relative gain reaches a plateau around 30% of gain at maximum.

Increased sensitivity to the BOLD effect and enhanced retinotopic maps

Tab. 2 summarizes the systematic boost in the number of activated voxels and maximum z-score values when static and up-to-the-firstorder dynamic B0 fluctuations imperfections were corrected. The number of activated voxels extracted using the thresholding alternative (i) (resp, (ii)) is, on average, 43.3% ± 17.2% (resp, 159,3% ± 38.6%) larger. The reported figures are consistent between the first and third participants. V#2, however, reveals fewer activated voxels, especially after FDR correction. This is likely due to larger (strong) head movement amplitudes. Additionally, Fig. 2 shows an example (V#1) of the projection of the retinotopic maps on the cortical surface. We notice clearly that the two visual hemifields project onto the contra-lateral hemispheres in the occipital cortex, a well-known mirroring feature of the primary visual cortex. However, a larger spatial extent of the retinotopic maps is associated with B0 field imperfections correction. Furthermore, the zoomed-in region of interest (ROI) in the left hemisphere illustrates a sharp gain in sensitivity along the cortical surface, suggesting that the activated voxels recovered when performing B0 perturbations correction are highly relevant.

DISCUSSION AND CONCLUSION

The systematic gain in tSNR on average associated with B0 imperfections correction we observed is generally consistent with previous studies even though the reported figures vary [START_REF] Vannesjo | Retrospective correction of physiological field fluctuations in high-field brain MRI using concurrent field monitoring[END_REF][START_REF] Schwarz | Correction of physiological field fluctuations in high-and lowresolution 3D-EPI acquisitions at 7 Tesla[END_REF]. Such differences arise from the experimental conditions (2D vs 3D acquisition sequence, 3T vs 7T B0 magnetic field, number of scans, volunteer's movement, etc.). The increase in sensitivity associated with B0 field imperfections correction we observed, notably at a more conservative statistical threshold is likely due to the gain in tSNR. Furthermore, the recovery of lost signal plays a major role in retrieving activated voxels. The enhanced image quality is also expected to yield better spatial specificity as blurring is extremely reduced.

In this work, we consider that the static and dynamic field terms evolve independently from each other since such an approximation is easy to implement and remains accurate enough. Nevertheless, it does not reflect the physics of the experiment faithfully. In fact, a truthful model would consider static and dynamic B0 imperfections Fig. 2. The projection of the retinotopic maps collected in V#1 and associated with data reconstructed without and with static and upto-the-first order field terms correction on the cortical surface. These maps were computed for an activation mask corresponding to a pvalue of 0.001 and without correcting for multiple comparisons correction (i.e. alternative (i)).

as evolving jointly: A ∆B0 map could be estimated for each volume in this case analogously to what is proposed in [START_REF] Dymerska | A method for the dynamic correction of B0-related distortions in single-echo EPI at 7 T[END_REF][START_REF] Haskell | FieldMapNet MRI: Learning-based mapping from single echo time BOLD fMRI data to fieldmaps with model-based reconstruction[END_REF]. Such a solution should also account for strong head movements to be complete. Furthermore, we performed typical sequential volume-wise reconstructions where the consecutive blocks of k-space data are considered independent from one another. Despite being easy to implement and efficient, such a strategy neglects the temporal structure of the fMRI scans. Global strategies that consider the entire scan data to leverage the temporal structure of the fMRI signal during reconstruction have also been proposed and can be further combined with B0 perturbations correction: Typically, low/fixed rank plus sparse or subspace-constrained reconstructions [START_REF] Andrii Y Petrov | Improving temporal resolution in fMRI using a 3D spiral acquisition and low rank plus sparse (L+ S) reconstruction[END_REF] can be carried out in the same framework, which could also improve the current temporal resolution.

In this work, we demonstrated a systematic and significant benefit in image quality, tSNR, and statistical sensitivity of brain activity detection and localization when correcting static and dynamic B0 inhomogeneities. We can further improve the image reconstruction process by considering a global 4D ∆B0 map or more advanced priors and even soon deep learning self-supervised strategies, which are still not available for 3D fMRI data at this time.

COMPLIANCE WITH ETHICAL STANDARDS

The experimental protocol was approved by the national ethics committee under the protocol identifier CPP 100048 issued by the national Comité de Protection des Personnes (CPP Sud Méditerranée 4 number 180913, IDRCB:2018-A011761-53). All participants gave their written informed consent.

Table 2 .

 2 Increase in the number of activated voxels detected and the maximum z-score values extracted from task-based fMRI scans with/out correcting B0 imperfections. The highest values (in bold font) are obtained when static and up-to-the-first-order B0 terms -Full Correction-are used. V#2 reveals the lowest statistical significance.

	Volunteer	p < 0.001 uncorrected	p < 0.05 FDR corrected	Z-score (max)
	#1	6456 9506 3253 7722	8.02	9.95
	#2	4367 7405 1059 2204	6.75	7.59
	#3	8529 10 823 2503 7745 10.83 12.08
	Average	6450 9244 2271 5890	8.54	9.87

https://github.com/hbp-brain-charting/public_ protocols

https://github.com/CEA_COSMIC/pysap-mri

https://www.fil.ion.ucl.ac.uk/spm/doc/biblio/