Ankush Agarwal
email: ankush.agarwal@polytechnique.edu

Stefano De Marco
email: stefano.demarco@polytechnique.edu

Emmanuel Gobet
email: emmanuel.gobet@polytechnique.edu

Gang Liu
email: gang.liu1988@gmail.com

Study of new rare event simulation schemes and their application to extreme scenario generation

Keywords: Extreme scenario, rare event, Markov chains, ergodic properties, interacting particle systems AMS subject classifications (2010): 65C05, 65C40, 37M25, 65C35

Introduction

Our aim is to estimate quantities such as P(X ∈ A), E(ϕ(X)|X ∈ A) and E(ϕ(X)1 X∈A ) where X is a random variable which takes value in a measurable state space X and A ⊂ X is a small subset. Another problem of interest is to sample from the conditional distribution X|X ∈ A, i.e. the generation of extreme scenarios. These quantities readily appear in the context of reliability for various industrial fields such as communication systems, aircraft safety, nuclear reactors safety, and so forth. Furthermore, they are also related to financial risks, for example model, credit and actuarial risk. Notably, the problem of extreme scenario generation is known as stress-testing in the field of financial risk management. The current methodology is to introduce artificial market specific and idiosyncratic shocks in the model and then to compute the default probabilities/the expected loss. These handpicked scenarios are sometimes simplified and can misrepresent risks in different ways [START_REF] Breuer | Systematic stress tests with entropic plausibility constraints[END_REF] and, in a dynamic setting, can fail to capture the interdependence between various stochastic factors in the market (see Section 4.2 in [ECB08]). In order to study the stability of a financial -say, banking -system, it is essential to generate extreme scenarios given the bank/the system default. In this work (Section 3), we present some examples where this is achieved.

In the case of a rare event, where the probability P(X ∈ A) is extremely small, simple Monte Carlo method is inefficient. Alternatively, we can use Importance Sampling methods which can be very efficient but require information about the model and the rare event. The schemes based on Interacting Particle System (IPS) or splitting techniques or ergodic transformations are also well developed. These approaches can be combined. For a detailed description, see [START_REF] Gobet | Rare event simulation using reversible shaking transformations[END_REF]. We elaborate here on the splitting and shaking transformation approach designed therein. To implement this approach, we first create a nested sequence of more and more rare events:

X = A 0 ⊃ • • • ⊃ A k ⊃ • • • ⊃ A n = A.
Then, we can write the quantities above via the following products:

P(X ∈ A) = n k=1 P(X ∈ A k |X ∈ A k-1 ),
(1.1)

E(ϕ(X)1 X∈A ) = E(ϕ(X)|X ∈ A n ) n k=1 P(X ∈ A k |X ∈ A k-1 ). (1.2)
Thus, instead of estimating the rare event probability directly, we estimate each conditional probability (and expectation) on the right hand side. Two different schemes, Parallel-One-Path (POP) and Interacting Particle System (IPS), using reversible shaking transformation, are proposed in [START_REF] Gobet | Rare event simulation using reversible shaking transformations[END_REF] and further studied in [START_REF] Agarwal | Rare event simulation related to financial risks: efficient estimation and sensitivity analysis[END_REF]. POP is based on the time average of a Markov chain while IPS is based on the space average of a large Markovian particle system. As is shown in the above references, the shaking transformation is the tool which unifies IPS and POP methods. A similar approach to POP has been developed in the engineering sciences community (see [START_REF] Au | Estimation of small failure probabilities in high dimensions by subset simulation[END_REF]), and is called in that context subset simulation approach3 . We emphasize that the shaking transformation designed in [START_REF] Gobet | Rare event simulation using reversible shaking transformations[END_REF] avoids using explicit transition kernels and is presented with a simple parametrization. It facilitates quick implementation and makes it easier to tune shaking forces. Our framework can also handle general situations including stochastic processes driven by Brownian motion or Poisson process. Furthermore, we have developed an adaptive version of POP method and proved its consistency in [START_REF] Agarwal | Rare event simulation related to financial risks: efficient estimation and sensitivity analysis[END_REF].

In this paper, our aim is threefold:

• We provide efficient variants of POP and IPS methods. The IPS algorithm is revisited by incorporating extra resampling to enhance the independence property between particles and by reducing the size of the particle system to keep the same computational cost. The adaptive POP algorithm is modified to preserve the benefits of both adaptivity and parallelization.

• When dealing with stochastic processes, we test the robustness of the two algorithms w.r.t. the dimension of the discretization of the underlying noise (i.e. the number of Brownian increments) in order to address the issue of the curse of dimensionality. As a difference with usual MCMC algorithms where the dimension effect may be significant (see for example [START_REF] Pillai | Optimal scaling and diffusion limits for the Langevin algorithm in high dimensions[END_REF]), we do not observe any significant decrease of speed of convergence.

• We illustrate the application of these methods for extreme (rare-event) scenario generation.

2. Parallel-One-Path (POP) and Interacting Particle System (IPS) Method

Shaking transformation

POP and IPS methods are designed to estimate each term on the r.h.s. of (1.1) using a random transformation called shaking transformation (a local distribution-preserving perturbation). For a given random variable X ∈ X , K(•) is called a shaking transformation for X if (X, K(X)) d = (K(X), X). We will consider the case where K(X) = K(X, Y ), where K(•, •) is a deterministic function and Y is another random variable independent of X. Then, we can build a X -valued Markov chain as

X i+1 = K(X i , Y i ), X 0 d = X,
where Y i are independent copies of Y . This Markov chain has the distribution of X as its stationary distribution. Thus, we expect an ergodic theorem to hold, i.e. for any measurable Ψ : X → R, setting

Z i := Ψ(X i ) and Z = Ψ(X) we have 1 N N i=1 Z i → N →+∞ E (Z) a.s. (2.1)
as soon as Z is integrable. The convergence rate of the above approximation depends on X and K. In the Gaussian case, if X = (X 1 , . . . , X j , . . . ) is an infinite-dimensional centered Gaussian vector, we define

K(X, X ) j := ρ j X j + 1 -ρ 2 j X j (2.2)
given an independent copy X of X and the constants ρ j ∈ [-1, 1]. Then, we show the following:

Theorem 1 (Theorem 2 [ADGL15]). Assume Z = Ψ(X) is square integrable. If |ρ| := sup j |ρ j | < 1, then E 1 N N i=1 Z i -E (Z) 2 ≤ Var (Z) N 1 + |ρ| 1 -|ρ| . (2.
3)

The proof for the above general case can be conducted with the help of the generalized Gebelein inequality. In the following, we provide an elementary proof in one dimension (X d = X d = N (0, 1)) using the properties of Hermite polynomials. Recall that the j-th Hermite polynomial is defined by H 0 := 1, H j (x) := (-1) j j! e

x 2 2 d j dx j (e -x 2 2 ), j ≥ 1, and it satisfies the following properties:

H j (x) = H j-1 (x), (j + 1)H j+1 (x) = xH j (x) -H j-1 (x). Lemma 1. For ρ ∈ [-1, 1], we have E H j (ρx + 1 -ρ 2 X ) = ρ j H j (x).
Proof. Obviously, the result holds for j = 0 and j = 1. Now suppose it holds for 1, 2, . . . , j. We prove it for j + 1.

(j + 1)E H j+1 (ρx + 1 -ρ 2 X ) = E (ρx + 1 -ρ 2 X )H j (ρx + 1 -ρ 2 X ) -H j-1 (ρx + 1 -ρ 2 X ) = ρxρ j H j (x) + E 1 -ρ 2 X H j (ρx + 1 -ρ 2 X ) -ρ j-1 H j-1 (x).
It is easy to prove that for a smooth function h with polynomial growth, we have

E (h (X )) = E (X h(X )) . Then take h(x ) := H j (ρx + 1 -ρ 2 x ) and write E 1 -ρ 2 X H j (ρx + 1 -ρ 2 X ) = (1 -ρ 2 )E H j (ρx + 1 -ρ 2 X ) = (1 -ρ 2 )ρ j-1 H j-1 (x),
which allows to conclude.

Proof of Theorem 1 (dimension 1). We consider the case X d = N (0, 1). We denote the single shaking parameter by ρ and define

T ρ Ψ(X) := E (Ψ(K(X, X ))) = E Ψ(ρX + 1 -ρ 2 X )
where the expectation is taken only w.r.t. X .

Since Hermite polynomials form a complete orthonormal system in L 2 (X), there exist

(a j ∈ R, j ≥ 1) such that Z = Ψ(X) = E (Ψ(X))+ j≥1 a j H j (X).
Then, the above lemma gives

T ρ (Ψ(X)) = E(Ψ(X)) + j≥1 ρ j a j H j (X); E T ρ (Ψ(X)) -E (Ψ(X)) 2 ≤ j≥1 ρ 2j a 2 j ≤ |ρ| 2 Var (Ψ(X)).
It readily follows that, for i < i ,

Cov (Z i , Z i ) = E Z i -E (Z) T ρ i -i (Ψ(X i )) -E (Z) ≤ |ρ| i -i Var (Z) .
This allows us to get the upper bound (2.3), since the l.h.s. of (2.3) is equal to 1

N 2 1≤i,i ≤N Cov (Z i , Z i ).
Using the shaking transformation, we define a second random transformation, called shaking with rejection, which allows to formulate POP and IPS methods. For a fixed k ∈ {0, . . . , n}, we define

M K k (X) := K(X)1 K(X)∈A k + X1 K(X) / ∈A k . (2.4)
As earlier, we write

M K k (X) = M K k (X, Y ) when K(X) = K(X, Y
) and Y is another random variable independent of X.

POP method

POP method is designed to estimate each conditional probability P(X ∈ A k |X ∈ A k-1 ) separately. The key observation for this technique is that the conditional distribution X | X ∈ A k is invariant with respect to the shaking with rejection

M K k (•). Proposition 2 ([GL15, Proposition 2.2]). Let k ∈ {0, 1, • • • , n}. For any bounded measurable ϕ, we have E ϕ(M K k (X)) | X ∈ A k = E (ϕ(X) | X ∈ A k ) . (2.5)
This proposition enables us to define a Markov chain which has X|X ∈ A k as its invariant measure. Definition 1. For each k = 0, . . . , n, given a starting point X k,0 , define

X k,i := M K k (X k,i-1 ) = M K k (X k,i-1 , Y k,i-1 ) for i ≥ 1 (2.6) where (Y k,i ) i≥0 is a sequence of independent copies of Y , independent of X k,0 .
Therefore, if the above defined Markov chain is ergodic, from Birkhoff ergodic Theorem, we have the following approximation:

E (ϕ(X)|X ∈ A k ) ≈ 1 N N -1 i=0 ϕ(X k,i ), as N → +∞. (2.7) Setting ϕ ≡ 1 A k+1 for k ≤ n-1 yields an approximation of P (X ∈ A k+1 |X ∈ A k )
and the product of all the estimators provides an estimate of the rare event probability P (X ∈ A) thanks to (1.1). Notice that POP method is especially suitable for parallel implementation and it is very economical in terms of memory consumption.

Application to stress testing. Another application of the above approach is to generate samples for X | X ∈ A by using the last level Markov chain (X n,i ) i≥0 . We know that for an aperiodic ergodic Markov chain, its marginal distribution will converge to its unique stationary distribution. Thus, if we run the Markov chain X n,• with enough iterations, its marginal distribution will be close to the conditional distribution of X | X ∈ A. In other words, for N large, (X n,N , X n,2N , . . . , X n,lN ) will constitute a sample of size l of (almost independent) random variables with distribution (approximately equal to) X | X ∈ A. This can be used to build extreme scenarios, for example, in stress testing. As recalled in the introduction, this is an important question in risk assessment. This application is illustrated later in Section 3.

Adaptive POP. A suitable choice of nested subsets (A k ) 0≤k≤n is essential for an efficient application of POP method. If any of the conditional probabilities

P(X ∈ A k | X ∈ A k-1
) is small, we again fall in the case of rare event.

Ideally, we would wish to have all the conditional probabilities to be of the same order but this is not possible without some a priori knowledge. In order to overcome this difficulty, an adaptive version of POP method is developed in [START_REF] Agarwal | Rare event simulation related to financial risks: efficient estimation and sensitivity analysis[END_REF] together with a proof of its convergence. This adaptive version finds good nested subsets using the quantile estimators of each intermediate level Markov chain. However the subsets (A k ) are determined sequentially which rules out a parallel implementation. In subsection 2.4.1 we propose a variant of this adaptive POP method which allows to recover a fully parallel version of the algorithm.

IPS method

With the aim of estimating the rare event probability P (X ∈ A), we define a Markov chain as follows:

Definition 2. (X i ) 0≤i≤n-1 is a X -valued Markov chain, such that X 0 d = X, X i := M K i (X i-1 ) = M K i (X i-1 , Y i-1 ) for 1 ≤ i ≤ n -1, (2.8)
where (Y i ) 0≤i≤n-2 is a sequence of independent copies of Y , independent of X 0 .

With the above definition we obtain the following result. Theorem 3 ([GL15]). For any bounded measurable function ϕ : X → R, we have

E (ϕ(X)1 X∈A ) = E ϕ(X n-1 ) n-1 i=0 1 A i+1 (X i ) .
(2.9)

The interpretation of rare event probability as an expectation related to a Markov chain, i.e. in a Feynman-Kac measure way, enables the use of IPS method. The steps of this algorithm are described below in Algorithm 1 with J = 1. The idea behind this procedure is to imitate the evolutionary selection and mutation process via the interaction of a large number of particles (see, for example, [Del04]). As seen from Algorithm 1, each generation is used to approximate the conditional distribution X | X ∈ A k .

Variants of POP and IPS methods

POP method with level refinement

As shown through the numerical examples in [START_REF] Gobet | Rare event simulation using reversible shaking transformations[END_REF], POP method is efficient when provided with good choices of nested subsets, which is possible thanks to the aforementioned level adaptive version. However, the advantage of parallelization is lost. Here, we propose a variant with level refinement which allows to recover the possibility of parallelization.

To begin, we fix a threshold value q (for example q = 0.05) for each conditional probability under which the estimation using one Markov chain is considered insufficient, and we arbitrarily choose some nested subset A k (such as an equi-distant partition of the entire space according to some criteria function). Then, we run one Markov chain as defined in Definition 1 at each level in parallel to estimate all the conditional probabilities P(X ∈ A k | X ∈ A k-1 ). Next, we check if the estimator p (N ) k is larger than q and accept the estimate for each level for which it is true. Otherwise, if p (N ) k 0 < q for some k 0 , we change the algorithm to the adaptive scheme as in [START_REF] Agarwal | Rare event simulation related to financial risks: efficient estimation and sensitivity analysis[END_REF] to estimate P(X ∈ A k 0 | X ∈ A k 0 -1 ). Notice that the Markov chain already used for level k 0 need not be simulated again. Another scheme a way to reduce this dependency between generations. As previously seen, the shaking with rejection transformation is invariant with respect to the conditional distribution of X. Hence, if we apply several iterations of the transformation to obtain the next particles generation, the distribution of the system will be less influenced by the previous state. In order to keep the same computational cost, we reduce the size of the particle system. Thus, we run the proposed version of IPS algorithm with particle size M J and the transformation M K k (X) applied J times at each time step. Regarding the convergence analysis, we can show that for a given J, the convergence still holds as M J goes to infinity. As we see in the numerical experiments, this variant of IPS method has a better performance, compared to the standard algorithm (Subsection 2.3) without extra resampling (J = 1).

Numerical experiments

In this section, we provide new and complementary numerical tests on two examples chosen from [START_REF] Agarwal | Rare event simulation related to financial risks: efficient estimation and sensitivity analysis[END_REF]. For completeness, each example is recalled briefly.

Credit portfolio default probability

We first consider an example of a credit portfolio based on the asset values of N 0 different firms (S 1 (t), . . . , S N 0 (t) : 0 ≤ t ≤ T ). Each S i follows a Heston model with the same CIR volatility (σ(t) : 0 ≤ t ≤ T ) and with the same parameters, given by

N 0 S i (0) r ρ W σ(0) κ σ γ ρ σ T
125 90 0.06 0.10 0.4 3.5 0.4 0.7 -0.06 1

The parameter r stands for the drift of S i , γ for the volatility of volatility, σ for the long-term volatility, κ for its mean-reverting rate; the correlation between the Brownian motions driving S i and S j is ρ W and the correlation between each asset and its volatility is ρ σ .

Our aim is to calculate the probability that at least L defaults occur before T Table 3.1: Empirical standard deviation of IPS estimators of default probability for L = 100 and B = 36 based on 1000 algorithm macro-runs, with M = 10 4 and particle system size equal to M = M J . The last column is the empirical standard deviation of POP method using n = M = 10 4 iterations at each level.

P(L) = P N 0 i=1 1 {τ i (B)≤T } > L , 0 < L < N 0 , (3.1) ×10 -6 J = 1 J = 2 J = 3 J = 4 J = 5 J = 6 J = 7 J = 8 J = 9 J = 10 POP ρ = 0.
where τ i (B) is the default time for the firm i with default level B:

τ i (B) := inf t ≥ 0 : S i (t) ≤ B .
Different IPS-based methods have been proposed in [START_REF] Carmona | Interacting particle systems for the computation of rare credit portfolio losses[END_REF] and [START_REF] Agarwal | Rare event simulation related to financial risks: efficient estimation and sensitivity analysis[END_REF] to compute P(L) using a dynamic and a static point of view, respectively. We will continue with the static point view, i.e. applying shaking transformation on the path space. In order to express P(L) in the form of (1.1), we need to create a cascade of decreasing sets (A k ) 1≤k≤n . We set

A k := N 0 i=1 1 {τ i (S i (0)-(S i (0)-B i ) k n )≤T } > L , 1 ≤ k ≤ n.
We estimate P(L) for different values of L with N step = 50 time steps in the Euler discretization scheme for (S 1 , . . . , S N 0 , σ) as explained in [START_REF] Agarwal | Rare event simulation related to financial risks: efficient estimation and sensitivity analysis[END_REF].

The random variable X modeling the rare event is Gaussian, associated to the N step increments of the N 0 +1 Brownian motions and we use the Gaussian shaker (2.2) with ρ j = Cst.

For L = 100 and B = 36, a 99% confidence interval [4.92, 5.13] × 10 -6 is given in [START_REF] Agarwal | Rare event simulation related to financial risks: efficient estimation and sensitivity analysis[END_REF] as reference. In Table 3.1, we report the results for IPS estimators with different number J of shaking transformation iterations applied at each time for fixed n = 5 levels. As we can see from Table 3.1, the IPS method with fewer particles but extra resampling at each step have better performance compared to the case without resampling. Heuristically, a good choice is J = 4. However, even after extra resampling, POP method yields smaller standard deviation.

Typically, in classical IPS method such as in [START_REF] Carmona | Interacting particle systems for the computation of rare credit portfolio losses[END_REF], the choice of appropriate discretization step (∆ t ) in the Euler scheme is quite delicate. On one hand, small ∆ t reduces the bias but on the other hand, small ∆ t adds statistical errors because the number of selection/mutation steps increases in the IPS method used there. In our version of IPS and POP methods, the underlying random variable X lies in a large-dimensional space X = R (N 0 +1)Nstep . Therefore it is important to assess how the large dimension affects the statistical errors of the methods (like in MCMC sampler, for example, see [START_REF] Pillai | Optimal scaling and diffusion limits for the Langevin algorithm in high dimensions[END_REF]). The important point is that we use reversible transformation directly in the path space. This suggests that our methods are less sensitive to time-discretization. We investigate this problem in Figure 3.1, where we study the impact of N step in POP method with n = 5 levels. We report the numerical results only for the POP method, as the qualitative phenomenon for IPS method is the same. The graphs show that the choice of N step , when large, has no significant impact on the optimal value of the shaking parameter ρ, which corresponds to the minimum standard deviation of the conditional probability estimator at a given level. Additionally, both the standard deviation and the rejection rate are quite insensitive to N step . These are important advantages of the methods studied in this work. This allows to decouple the problem of bias reduction and the control of statistical convergence.

Stress-testing. In Figure 3.2, we exhibit extreme scenarios for the asset price of firm 1 and its volatility, in the situation of various defaults (level 1 and last level, n = 5, L = 100) with two different values of B. Each scenario is obtained by using 10 4 iterations of shaking with rejection (with ρ = 0.9). For this, we have used the first and last-level Markov chain X 1,. and X 5,. , respectively, to get samples from the distributions X | X ∈ A 1 and X | X ∈ A respectively, as explained in Subsection 2.2. We believe that such tools can be efficiently exploited by regulators and risk managers for a better risk assessment.

Model misspecification in financial engineering and extreme scenarios

We consider the Profit&Loss (PL) when a trader uses a Black-Scholes (BS) model with volatility σ -> 0 in order to hedge a European call option, while the true dynamics of the underlying S is given by a path-dependent volatility model with volatility switching between σ -and σ + > σ -. We propose a discrete-time version of the model based on a monitoring period ∆ t (say 1 week) and monitoring dates t i = i∆ t : if the underlying spot price drops below the average of previous four monitored prices, the level of volatility becomes σ + , otherwise it sticks to σ -(see [START_REF] Agarwal | Rare event simulation related to financial risks: efficient estimation and sensitivity analysis[END_REF] for more details). This model corresponds to the usual empirical observation that the underlying volatility is higher when price falls. The call option maturity is T , and [0, T ] is the trading period under consideration. For our numerical experiment, we take T = 1, ∆ t = T /50 and we assume that the trader's rebalancing period is 5∆ t . At times t 5∆tj , 0 ≤ j < 10, the trader holds δ j assets where δ j is given from the BS-Delta with volatility σ -and spot S 5∆tj . Hence at the maturity the trader's PL is given by P L trader := E trader (S T -K strk ) + + 9 j=0 δ j (S 5∆t(j+1) -S 5∆tj ) -(S T -K strk ) + .

Because of the mismatch between the realized model and the hedging model, the trader may incur large losses. Large loss probability P (A) with A = {P L trader ≤ L} has been estimated using POP and IPS methods in [START_REF] Agarwal | Rare event simulation related to financial risks: efficient estimation and sensitivity analysis[END_REF].

Here we are more interested to know what are the typical scenarios which generate large losses. We set S 0 = K strk = 10, σ -= 0.2, σ + = 0.27 and take L = -2.4. Applying the principle of subsection 2.2 on stress-testing, we can get samples from the distribution X | X ∈ A; 5 typical scenarios are reported in Figure 3.3. Each scenario is obtained using 10 4 iterations of shaking with rejection (with ρ = 0.9). As intuitively expected, typical extreme scenarios exhibit large fluctuations of S.

Figure 3 . 1 :

 31 Figure 3.1: The impact of the number of discretization times (N step ) on the optimal shaking parameter (ρ) at each level of POP method.

Figure 3 . 2 :

 32 Figure 3.2: Sample paths for the asset price of firm 1 at Level 1 and Level 5 in the POP method and the respective volatility sample paths.

Figure 3 . 3 :

 33 Figure 3.3: Typical paths of the underlying stock price which lead to large hedging loss

We would like to thank Robert Scheichl for directing our attention to this reference.

The author research is part of the Chair Financial The author research is part of the Finance for Energy Market Research Centre.

is to put new nested subsets between A k 0 and A k 0 -1 and run POP method again in parallel. Finally, the product of all the estimators provides an estimate of the rare event probability.

IPS method with extra resampling

Initialization:

For given integers J and M , set

∈ A i+1 } and independently of everything else; end Mutation step:

i,j-1 ) where

Algorithm 1: IPS method with extra resampling and reduced size for computing P (X ∈ A)

An intrinsic feature of IPS method is that different generations of particles are correlated with each other. Thus, if the empirical measure is inaccurate at the first generation, it is likely to be inaccurate at the following generation, which amplifies the variance of the final estimator. Here, we propose


Member, IEEE Pierre Chainais
email: pierre.chainais@ec-lille.fr.

Aymeric Leray
email: aymeric.leray@u-bourgogne.fr

Pierre Pfennig

Statistical performance analysis of a fast super-resolution technique using noisy translations

Keywords: high-resolution imaging, reconstruction algorithms, super-resolution, performance evaluation, error analysis I

The registration process is a key step for superresolution (SR) reconstruction. More and more devices permit to overcome this bottleneck by using a controlled positioning system, e.g. sensor shifting using a piezoelectric stage. This makes possible to acquire multiple images of the same scene at different controlled positions. Then a fast SR algorithm [1] can be used for efficient SR reconstruction. In this case, the optimal use of r 2 images for a resolution enhancement factor r is generally not enough to obtain satisfying results due to the random inaccuracy of the positioning system. Thus we propose to take several images around each reference position. We study the error produced by the SR algorithm due to spatial uncertainty as a function of the number of images per position. We obtain a lower bound on the number of images that is necessary to ensure a given error upper bound with probability higher than some desired confidence level. Such results give precious hints to the design of SR systems.

However, the positioning system (or any registration method) only approximately reaches the targeted positions with some small random error. Based on a statistical performance analysis, we study the influence of this error on the quality of the SR images reconstructed with a simple and fast SR algorithm [START_REF] Elad | A fast super-resolution reconstruction algorithm for pure translational motion and common space-invariant blur[END_REF] which assumes that displacements are exactly known. We also study the importance of using several acquisitions of the same targeted positions to compensate for postioning errors in order to optimize the number of images required to ensure a given quality of the SR image. While the chosen SR method is a priori not as efficient as state of the art methods [START_REF] Yang | Super resolution imaging, ch. Image superresolution: historical overview and future challenges[END_REF], the theoretical analysis of its statistical performance is possible, which would not likely be the case for other methods. Therefore, in addition to its rapidity, this method would come with theoretical guarantees on the quality of reconstruction. Moreover the adopted methodology paves the way to the analysis of more sophisticated SR methods, which is of great importance to give hints on the optimal coconception of integrated SR imaging sytems.

Over the last 30 years, several works have dealt with mathematical analysis of SR algorithms, e.g. [START_REF] Ng | Analysis of displacement errors in high-resolution image reconstruction with multisensors[END_REF]- [START_REF] Lin | Fundamental limits of reconstruction-based superresolution algorithms under local translation[END_REF]. The works described in [START_REF] Ng | Analysis of displacement errors in high-resolution image reconstruction with multisensors[END_REF]- [START_REF] Bose | Performance analysis of the tls algorithm for image reconstruction from a sequence of undersampled noisy and blurred frames[END_REF] essentially study the convergence of iterative methods for SR (e.g., conjugate gradient) including registration and deconvolution steps. They show that the reconstruction error decreases as the inverse of the number of LR images. In [START_REF] Baker | Limits on super-resolution and how to break them[END_REF], the difficulty of the inverse problem is characterized by the conditioning number of a matrix defined from the direct model which is proportional to r 2 s 2 (s = width of sensor pixels). When translations are uniformly distributed in (0, r) 2 , this conditioning number tends to 1 and a direct inversion is possible with high probability when a large number of images is used [START_REF] Traonmilin | On the Amount of Regularization for Super-Resolution Interpolation[END_REF]. In [START_REF] Champagnat | Statistical performance modeling for superresolution: a discrete data-continuous reconstruction framework[END_REF], the analysis was performed in the Fourier domain and showed that the mean square error decreases as the number of images increases when random translations are used. Ref. [START_REF] Robinson | Statistical performance analysis of superresolution[END_REF] quantifies the limitations of SR methods by computing Cramer-Rao lower bounds, also working in the Fourier domain. In the most favourable case where translations are known (no registration is needed), this bound is proportional to r/n if n is the number of images. All these works back to the 1980s [START_REF] Tsai | Multiframe image restoration and registration[END_REF] explain what makes SR difficult and how far more images can make it simpler. However, they have only expressed limited quantitative prediction beyond the qualitative 1/n behaviour of the reconstruction error. Our purpose is a detailed quantitative statistical error analysis of the simple Shift & Add method described in [START_REF] Elad | A fast super-resolution reconstruction algorithm for pure translational motion and common space-invariant blur[END_REF]. We obtain a lower bound on the number of images that is necessary to achieve a given error bound with high probability. The control of errors is crucial to produce nice looking results but also to ensure reliable scientific observations. The present study is performed in the Fourier domain. The error at each frequency component is quantitatively evaluated. The use of Hoeffding's inequality permits to compute upper bounds and confidence intervals of practical use are obtained.

A preliminary work was presented at ICASSP 2014 in [START_REF] Chainais | Quantitative control of the error bounds of a fast super-resolution technique for microscopy and astronomy[END_REF] with less general results because the assumptions were more restrictive (special uniform distribution). In this work, we use a more general and realistic assumption of bounded error on displacements. Furthermore, the potential presence of bias is taken into account and all mathematical proofs are given. The present results are tighter thanks to the use of Hoeffding's concentration inequality in place of the lose Bienaymé-Cebycev inequality. This article includes a numerical study and more detailed illustrations ; all useful Matlab codes are available.

Section II presents the setting and the model. Section III presents our main theoretical results which predict the required number of image acquisitions at each position to ensure some given confidence level in the reconstructed image. Sections III-A & III-B present the most technical aspects; proofs are in Appendix. Section III-C sums up our main theoretical results. Section IV presents numerical results. Section V discusses our contributions and some prospects.

II. A FAST AND CONTROLLED SUPER-RESOLUTION TECHNIQUE

A. The super-resolution problem

For a given SR factor r, the most common linear formulation of the general SR problem in the pixel domain is [START_REF] Elad | A fast super-resolution reconstruction algorithm for pure translational motion and common space-invariant blur[END_REF]:

Y k = D k H k F k Y HR + n k k = 1, .., K, (1)
where

Y HR = argmin Y K k=1 Y k -D k H k F k Y 2 2 . (2)
Other formulations based on the L1-norm or adding some regularization have also been proposed [START_REF] Farsiu | Advances and challenges in super-resolution[END_REF]. We focus on the method described in [START_REF] Elad | A fast super-resolution reconstruction algorithm for pure translational motion and common space-invariant blur[END_REF]: its simplicity makes it possible to quantitatively analyze its performances. Such a guarantee may be crucial for scientific imaging or the design of devices.

B. Super-resolution algorithm

Several usual assumptions are used in [START_REF] Elad | A fast super-resolution reconstruction algorithm for pure translational motion and common space-invariant blur[END_REF]. The PSF of the acquisition system is known and spatially homogeneous so that ∀k, H k = H. Decimation is the same for all images so that ∀k, D k = D in (1) & (2). We will also assume that the r 2 possible translated images at integer multiples (k, ℓ) ∈ (0, r -1) 2 of the HR scale are available to provide an optimal setting for SR [START_REF] Robinson | Statistical performance analysis of superresolution[END_REF]. Then the solution to the leastsquare error SR problem (2) consists of two steps. A blurred image Z = HY HR can be estimated by [START_REF] Elad | A fast super-resolution reconstruction algorithm for pure translational motion and common space-invariant blur[END_REF]:

(k ′ ) = [F HR X](k ′ ) Discrete Fourier Transform of HR image X Y (k) = [F LR Y ](k)
Z := HY HR = r 2 k=1 F t k D t Y k (3)
The operation in (3) is equivalent to a simple interlacing of LR images, see Fig. 1. Then the final HR image results from the deconvolution of Z, which can be done using any algorithm such as Wiener or Lucy [START_REF] Bovik | The Essential Guide to Image Processing[END_REF]. Such an approach separates the problem of SR into two steps of fusion (estimating Z) and deconvolution (deblurring to estimate Y HR ). This work focuses on the performance analysis of the fusion step only.

Recall that high frequency terms at some k ′ are preserved if and only if the PSF H(k ′ ) is not zero. Some prior information might be used to reconstruct missing frequencies [START_REF] Robinson | Statistical performance analysis of superresolution[END_REF]. This algorithm requires one idealized assumption: displacements (matrices F k ) are assumed to be exact integer multiples of HR pixels. In practice, this is only approximately true due to the finite precision of the positioning system. Our purpose is to study the influence of this approximation. One solution would be to carry out accurate sub-pixel registration. This would remain insufficient since state of the art techniques cannot ensure a precision much better than 0.1-0.01 pixel [START_REF] Foroosh | Extension of phase correlation to subpixel registration[END_REF], [START_REF] Robinson | Optimal registration of aliased images using variable projection with applications to super-resolution[END_REF]. Another possibility is to take n d ≥ 1 images for each required position so that the true Z will be replaced by the estimate:

X = Ẑ = d 1 n d n d j=1 (F d ) t D t Y de(j) (4)
Y d is the image of a scene Y translated by d where d = (d x , d y ) denotes the targeted displacement vector; d e (j) = d + b dj the real experimental displacement; b dj is the noise on the platform position. Note that in general (F d ) t F de(j) = I rN . One can hope to compensate from displacement inaccuracies by using multiple acquisitions at the same targeted position with some random error b dj around the expected value d. A realistic assumption [START_REF] Yamahata | Subnanometer translation of microelectromechanical systems measured by discrete fourier analysis of ccd images[END_REF], [START_REF] Lin | Fundamental limits of reconstruction-based superresolution algorithms under local translation[END_REF], [START_REF] Ben-Ezra | Video super-resolution using controlled subpixel detector shifts[END_REF] is that the position error is bounded by ǫ > 0 in LR pixel units or ǫ r = ǫr in HR pixel units. For a given targeted position, the positioning system will be reset between each acquisition so that positions are randomly distributed around the average position (which may be biased due to miscalibration). This averaging process is expected to enhance the SR quality. Fig. 2 illustrates typical results from this approach applied to a detail of Barbara for r = 2, ǫ = 0.1. The error on reconstructed high frequency components are compared for n d = 1 and n d = 32 images/position. A SNR gain of about 15dB is observed when using 32 images (note for later use that 10log 10 (32) = 15). Our aim is to reconstruct probably approximately correct (PAC) images by quantifying the number of images that should be taken per reference position to respect some given upper relative error bound of p (e.g. 0.10) with probability (confidence) higher than P (e.g. 0.90).

C. Aliasing effects and notations

To detail the effect of aliasing, we consider the relation between the estimated blurred HR image X defined by (4) and the LR images Y de(j) in the Fourier domain, see Fig. 1(b). For some integer n, the interval (-n : n) denotes the set of integers between -n and n (Matlab notations). When using the Discrete Fourier Transform (DFT), we denote by k the LR frequencies in D LR = (-N/2 : N/2 -1) 2 and k ′ the HR frequencies in D HR = (-rN/2 : rN/2 -1) 2 . Given some HR frequency k ′ , we need to deal with corresponding aliased terms in the LR image. The integer vector γ ∈ (-r : r) 2 is such that k = k ′ -γN ∈ D LR . We denote by α the integer vectors such that k + αN ∈ D HR . Sums d are over all the r 2 ideal displacements d ∈ (0 : r -1) 2 and sums over α are sums over all possible HR frequencies k α = k + αN (up to rN/2). The DFT of image Z is Z. To alleviate formulas, we introduce the normalized frequencies:

     q γ = 2π rN k ′ = 2π rN (k + γN ) = q + γ 2π r q α = 2π rN (k + αN ) = q + α 2π r (5)
where α, γ ∈ Z 2 . Note that q α ∈ (-π, π) 2 so that q α 1 ≤ 2π and q α 2 ≤ √ 2π. Back to (4), note that when D is the decimation operator, D t is an upsampling operation (inserting zeros between samples) that produces aliasing. If F HR is the HR DFT, for k ′ ∈ D HR :

[F HR D t Y de(j) ](k ′ ) = Y de(j) (k = k ′ -γN ) (6)
Taking phase shifts due to translations of (-d) associated to (F d ) t into account in the DFT of (4) yields:

X(k ′ ) = d 1 n d n d j=1 Y de(j) (k ′ -γN ) e 2iπ rN d•k ′ (7)
Since each observation is a decimated version of the blurred translated scene, one has in the spatial domain:

Y de(j) = DHF de(j) Y HR = DZ de(j) (8)
In the Fourier domain:

Y de(j) (k) = [F LR DF t HR Z de(j) ](k) (9)
and thanks to usual properties of the sum of roots of unity (see Appendix E):

Y de(j) (k) = 1 r 2 α Z(k α ) e -2iπ rN kα•de(j) (10)
where we have used the fact that the homogeneous blur operator (convolution) is diagonal in Fourier domain. One can explicitly see in [START_REF] Champagnat | Statistical performance modeling for superresolution: a discrete data-continuous reconstruction framework[END_REF] how the information at high frequencies k α = k + αN from the HR image is aliased at low frequency k in each LR image Y de(j) . By separating the desired main contribution at k ′ = k + γN and aliasing terms at k + αN for α = γ, one gets by reporting (10) in [START_REF] Bose | Performance analysis of the tls algorithm for image reconstruction from a sequence of undersampled noisy and blurred frames[END_REF]:

X(k ′ ) = Z(k ′ )G γ (k ′ ) + B(k ′ ) (11)
B(k ′ ) = α =γ Z(k + αN )G α (k ′ ) (12)
where

G α (k ′ ) = 1 r 2 n d d,j e -i 2π r (α-γ)d e iqα•b dj (13)
(except when k ′ x or k ′ y is equal to -rN/2). In the ideal case where b dj = 0 translations are exact multiples of HR pixels and one retrieves X = Z = H Y HR since G γ = 1 and G α = 0 for α = γ. The first term in [START_REF] Robinson | Statistical performance analysis of superresolution[END_REF] is the main approximation term, which should be as close as possible to Z(k ′ ). The second term B(k ′ ) in [START_REF] Robinson | Statistical performance analysis of superresolution[END_REF] is the aliasing term and should be as small as possible compared to the approximation term. Our purpose is to establish conditions for which X is a good approximation of Z within quantitative probabilistic bounds.

III. BOUNDS ON RECONSTRUCTION ERRORS

This section proves concentration inequalities that guarantee PAC SR. In this study, we make the general and realistic assumption that position errors are bounded so that b dj ∈ (-ǫ r , ǫ r ) 2 HR pixel units. We do not assume that IE[b dj ] = 0: the positioning system might be biased. In section III-A & III-B we deal with the coefficient G γ of Z in the main approximation term of [START_REF] Robinson | Statistical performance analysis of superresolution[END_REF]

A. Bound on the approximation term

G γ (k ′ ) Since one expects that 1 r 2 IEG γ ≃ 1, we start from |G γ (k ′ ) -1| ≤ (14) |G γ (k ′ ) -IE[G γ (k ′ )]| + |IE[G γ (k ′ )] -1| Noting that IE[G γ (k ′ )] = IE[e iqγ •b dj ], the Taylor development of the complex exponential function yields 1 IE[e iqγ •b d j ] -1 ≤ |q γ .IE[b dj ]| + IE (q γ .b dj ) 2 2 ≤ q γ 2 IE[b dj ] 2 + q γ 2 1 ǫ 2 r 2 B1 (15)
since b dj ∈ (-ǫ r , ǫ r ) 2 . Then we deal with the first term in ( 14) by introducing:

B G = 1 r 2 n d d,j (e iqγ •b d j -IEe iqγ •b dj ) (16)
To obtain concentration inequalities on |B G |, our approach goes in 3 steps: i) bound the real and imaginary parts thanks to properties of their power series expansions, ii) prove concentration inequalities by using Hoeffding's inequality for the sum of differences between random variables and their expectations, iii) bound |B G | by using Lemma 1 below to combine bounds on the real and imaginary parts. Lemma 1: (see proof in Appendix A) Let x 1 and x 2 two random variables in R. Let a 1 , a 2 > 0 and P 1 , P 2 ∈ (0, 1)

such that P (|x i | ≥ a i ) ≤ P i , i = 1, 2. Then P ( x 2 1 + x 2 2 ≥ a 2 1 + a 2 2 ) ≤ P 1 + P 2 (17)
P (|x 1 | + |x 2 | ≥ a 1 + a 2 ) ≤ P 1 + P 2 (18)
Let us recall Hoeffding's inequality. Hoeffding's inequality [START_REF] Boucheron | Concentration inequalities[END_REF].

Let {X i , 1 ≤ i ≤ n} a set of independent random 1 See Lemma 1 p. 512 in Feller (vol. 2) [20] on the Taylor development of exp(it) for t > 0. variables distributed over finite intervals [a i , b i ]. Let S = n i=1 (X i -IE[X i ]). For all t > 0, P (|S| ≥ t) ≤ 2 exp - 2t 2 n i=1 (b i -a i ) 2 (19)
This permits to prove that, see Appendix B:

P      |B G | ≥ √ 2 δ γ + q γ 3 1 ǫ 3 r 3 B2      ≤ 4e -c 2 n d /8 (20)
We obtain the final concentration inequality for the main approximation term by combining ( 15) and ( 20) and going back to [START_REF] Chainais | Quantitative control of the error bounds of a fast super-resolution technique for microscopy and astronomy[END_REF]:

P (|G γ (k ′ ) -1| ≥ B 1 + B 2 ) ≤ 4e -c 2 n d /8 (21)
for ǫ r ≤ 1/πr, where B 1 and B 2 are defined in ( 15) & [START_REF] Feller | An Introduction to Probability Theory and Its Applications[END_REF]. Let p ∈ (0, 1) the maximum relative error constraint, e.g., p = 0.1, and P 1 ∈ (0, 1) such that 1 -P 1 is the corresponding concentration probability. For sufficiently large p, one can define ∀k ′ ∈ D HR or q γ ∈ 2π rN D HR the adequate maximum coefficient c(q γ ) > 0 such that, neglecting the cubic term,

√ 2c(q γ ) q γ 1 ǫ+ q γ 2 IE[b d,j ] 2 d + q γ 2 1 ǫ 2 r 2 ≤ p (22)
c(q γ ) is a decreasing function of q γ 1 , which is minimum for maximal frequencies such that q γ 1 = 2π. For p large enough, one can define

c 1 (p) = min qγ c(q γ ) = c(π, π) (23) = 1 2 √ 2πǫ p - √ 2π IE[b d,j ] 2 d -2π 2 ǫ 2 r 2 Then (22) with c(q γ ) replaced by c 1 (p) is true for all q γ ∈ 2π rN D HR . When the averaged bias IE[e iqγ •b dj ] d is zero or remains negligible (≪ p/ √ 2π), c 1 (p) ≃ p -2π 2 ǫ 2 r 2 2 √ 2πǫ (24)
If ǫ ≤ 1/πr and c 1 (p) is well defined, (21) becomes:

P (|G γ (k ′ ) -1| ≥ p) ≤ 4 exp - c 1 (p) 2 n d 8 ( 25
)
∀k ′ ∈ D HR . Then the relative error remains bounded by p with probability larger than some P 1 ∈ (0, 1) if

n d ≥ 8 c 1 (p) 2 log 4 1 -P 1 (26)
The larger c 1 (p), the smaller the lower bound. This bound does not depend on the image content. In practice, it tells that, for n d large enough, the main approximation term in [START_REF] Robinson | Statistical performance analysis of superresolution[END_REF] is less than 100p% away from the targeted Z(k ′ ) with probability larger than P 1 . In ideal experimental conditions, with no bias and ǫr ≤ p/2π 2 ,

n d ≥ 8πǫ p -2π 2 ǫ 2 r 2 2 log 4 1 -P 1 ( 27
)
For instance, see Tab. II, for ǫ = 0.01, r = 2, p = 0.1 and P 1 = 0.90 (error ≤ 10% with ≥ 90% confidence level) this bound is n d ≥ 28. The concentration level (1 -P 1 ) can be very tight due to the logarithmic dependence of n d on (1 -P 1 ). At the same error level p = 0.1, the criterion becomes n d ≥ 45 for P 1 = 0.99. In contrast, a much larger n d ≥ 5.3 10 4 is necessary to guarantee an accuracy of 1% (p = 0.01) at P 1 = 0.90 confidence level. In summary, confidence is cheap while accuracy is expensive. Note that the position accuracy ǫ should essentially decrease proportionally to p as a finer reconstruction is desired. Moreover, given a desired SR factor r and a position accuracy ǫ, the relative error p is lower bounded by 2π2 ǫ 2 r 2 . For r = 2 and ǫ = 0.01, the smallest relative error p that can be guaranteed is p best = 0.008.

B. Bound on the aliasing terms (G

α , α = γ)
The ideal situation in [START_REF] Robinson | Statistical performance analysis of superresolution[END_REF] occurs when the translations d are exactly the r 2 possible multiples of HR pixels. Due to properties of complex roots of unity, all the aliasing terms G α (k ′ ) in ( 11) cancel for α = γ. Our aim is to bound the contribution of aliasing error terms when translations are noisy due to approximate control only. The adopted strategy is similar to that of previous section, see proof in Appendix C. We also use the properties of roots of unity and a standard assumption on the spectral content of the target image. We start from [START_REF] Lin | Fundamental limits of reconstruction-based superresolution algorithms under local translation[END_REF]:

G α (k ′ ) = 1 r 2 n d d,j e -i 2π r (α-γ)d e iqα•b dj (28)
Let

θ αd = 2π r (α -γ)d, d ∈ (0 : r -1) 2 (29)
Note that the set of the e iθ αd matches the set of products of complex roots of unity, see eq. ( 95)-(98) in Appendix E. The sum over translations d involves the sum of roots of unity, which is zero, in the computation of the aliasing term. In Appendix C, assuming that the variations of the bias IE[q α • b dj ] around IE[q α • b dj ] d for fixed d are negligible, we prove the following concentration inequalities. For αγ / ∈ {0, r/2} 2 :

P |G α (k ′ )| ≥ √ 2δ ′ α ≤ 4e -c 2 n d (30)
For αγ ∈ {0, r/2} 2 , (79) in Appendix C gives a deterministic bound on the real part. Moreover sin(θ αd ) = 0 in (80) so that one gets from [START_REF] Robinson | Optimal registration of aliased images using variable projection with applications to super-resolution[END_REF] in Lemma 1:

P (|G α (k ′ )| ≥ δ ′ α ) ≤ 2e -c 2 n d (31)
which is even tighter than (30). In the special case r = 2, all αγ are in {0, r/2} 2 = {0, 1} 2 so that we need (31) only and tighter bounds are obtained. We aim at taking into account the contribution of all terms Z α G α (k ′ ) for α = γ in [START_REF] Tsai | Multiframe image restoration and registration[END_REF]. Let assume that they are independent. This is at least approximately true for two main reasons. First one can show that the G α (k ′ ) are uncorrelated, see (111) in Appendix G and second the Z α carry information about very distinct frequencies in the image. Then we can use Lemma 2 (see proof in Appendix A): Lemma 2: Let x i , i = 1, ..., n, n independent random variables. Let a i > 0 and P i ∈ (0, 1) i = 1, ..., n, such that ∀i, P (|x i | ≥ a i ) ≤ P i . Then

P i |x i | ≤ i a i ≥ n i=1 (1 -P i ) (32)
Applying Lemma 2 to the set of (r 2 -1) possible α = γ from (30) yields a probabilistic bound on the relative aliasing error when Z γ = 0 2 :

P   α =γ Z α Z γ G α (k ′ ) ≤ √ 2 α =γ Z α Z γ δ ′ α   ≥ 1 -4e -c 2 n d r 2 -1 (33)
Given some desired relative error p ∈ (0, 1) and lower probability P 2 ∈ (0, 1), one needs to find whether there exists

c = c 2 (p) > 0 such that ∀k ′ ∈ D HR √ 2 α =γ Z α Z γ (c q α 1 ǫ + f (q α , ǫ r )) δ ′ α ≤ p, (34)
A necessary condition appears as

p > p 0 (ǫ, r, Z) = √ 2 α =γ Z α Z γ f (q α , ǫ r ) (35)
Then one can define

1 -4e -c2(p) 2 n d r 2 -1 ≥ P 2 , (37)
that is

n min d = 1 c 2 (p) 2 log   4 1 -P 1 r 2 -1 2   (38) (39)
In the special case r = 2, (31) yields the even tighter bound:

n min d = 1 c 2 ( √ 2p) 2 log 2 1 -P 1 3 2 . ( 40
)
One obtains a bound on the aliasing error relative to | Z(q γ )|:

P   α =γ Z α Z γ G α (k ′ ) ≤ p   ≥ P 2 (41)
This relative error provides a good estimate of the relative error on the HR image before deconvolution. It permits to evaluate the contribution of aliasing errors to the reconstructed blurred HR image Z. This necessitates the knowledge of the true HR image : one can also use the reconstructed image a posteriori to indicate which frequencies are most suspected to contribute to aliasing effects. Each specific image has a specific Fourier spectrum so that special aliasing effects may appear and make SR difficult, at least for a small set of frequencies for which the sum of aliasing terms in (41) may be particularly large. To propose a generic a priori estimate of the order of magnitude of this aliasing error, we need to make some assumptions on the content of images. It is well accepted that natural images often exhibit a power law energy spectrum ∝ 1/ k ′ 2(1+η)

2
where usually |η| ≪ 1 [22]- [START_REF] Chainais | Infinitely divisible cascades to model the statistics of natural images[END_REF]. Then

Z α Z γ = | H(q α )| | H(q γ )| q γ 2 q α 2 1+η ( 42
)
Therefore the strongest constraints are due to high frequencies (large k ′ or q γ ). Note the dependence on the blur kernel which acts as a low-pass filter: the presence of H in (42) will have adverse effects. Searching for lower-bounds, forthcoming computations consider the most favourable case when H = 1.

See section IV for a numerical illustration of the effect of a realistic Gaussian blur kernel. An approximate computation in App. D shows that the highest frequencies define c * 2 (p) as

c * 2 (p) = p -p * 0 (ǫ, r) a * (ǫ, H) (43)
where

p * 0 (ǫ, r) ≃ b 0 √ 2 η π 2 ǫ 2 r 2 (r 2 -1) (44)
a * (ǫ, H) = √ 2 α =γ | H(q α )| | H(q γ )| q γ 2 q α 2 1+η q α 1 ǫ ≃ a 0 2 1+η/2 ǫ(r 2 -1) (if H = 1) (45)
where the factor (r 2 -1) corresponds to the number of aliasing terms; the coefficient b 0 ≃ 2/3 for r = 2 and b 0 ≃ 1.2 for r ≥ 3 and it is almost independent of the size N of the image for N ≥ 32; a 0 ≃ 0.63 for r = 2 and a 0 ≃ 1.3 for r ≥ 3 (see Appendix D). In the general case, (44) & (45) interestingly permit to make explicit the dependence on r, ǫ and η. Thus, for a power-law spectrum image, the required minimum number n min d of images/position is:

n min d = a * 2 (ǫ, H) (p -p * 0 (ǫ, r)) 2 log   4 1 -P 1 r 2 -1 2   (46)
One observes that p 0 /ǫ 2 r 2 essentially depends on r as soon as ǫ is small enough. Figure 3 illustrates numerical orders of magnitude of reachable (p, ǫ) such that p > p * 0 (ǫ, r) for given r under the assumption of a power law spectrum. Pairs of acceptable parameters (p, ǫ) for which guaranteed error bounds exist are at the bottom right of each curve. Typical values can be evaluated numerically. For instance assuming η = 0, to guarantee an error smaller than 10%, r = 2, p = 0.1 ⇒ ǫ ≤ 0.036 or r = 6, p = 0.1 ⇒ ǫ ≤ 0.0035. Observe that ǫ should rapidly decrease as r becomes larger when some given error level p with high probability is desired. Note the logarithmic dependence on (1 -P

1 r 2 -1 2
) which permits to choose P 2 close to 1 without increasing n min d a lot. By using our results in the other way, one can also deduce a map of confidence intervals p(q) for fixed n d . In practice, the acquisition protocole may impose some fixed n d . Then one can set the value of c 2 (p) in (34) and compute a map of confidence intervals p(q) in the Fourier domain, taking into account the spectrum of the true HR image. Since it is not known, the Fourier transform may be replaced by its estimate. This procedure helps identifying which frequencies are more likely to contribute to aliasing errors.

C. Main results

The analysis of the estimate X of the blurred image Z = HY HR by the proposed algorithm gives in the spectral domain, see ( 11) & ( 12):

X(k ′ ) = Z(k ′ )G γ (k ′ ) + B(k ′ ) (47) B(k ′ ) = α =γ Z(k + αN )G α (k ′ ) (48)
Theorem 3 below gathers the necessary assumptions on the acquisition system (r, ǫ, IE[b dj ]), the scenes (spectrum exponent η in (42)) and the desired confidence level (p 1 & P 1 , p 2 & P 2 ) to obtain two fundamental concentration inequalities for the approximation and the aliasing terms respectively. Theorem 3: Acquisition system -Let r be the SR factor. Let 0 < ǫ < 1/πr be the maximum error of the positioning system (in LR pixel units). Assume bounded errors b dj on positions within (-ǫr, ǫr) in both x and y directions with a possible constant bias IE[b dj ] (in HR pixel units). Assume that n d images are taken for each one of the r necessary reference positions corresponding to d ∈ (0, r -1) 2 HR pixel units. Confidence intervals -Let p 1 ∈ (0, 1), resp. p 2 ∈ (0, 1) be the desired maximum relative error on the main approximation term, resp. the sum of aliasing terms, of the reconstructed image (p 1 & p 2 will generally be close to 0). Let P 1 ∈ (0, 1) be the desired level of confidence in the relative error p 1 due to the main approximation term. Let P 2 ∈ (0, 1) be the level of confidence in the relative error p 2 due to the aliasing term (P 1 and P 2 will be close to 1). Technical assumptions -Assume that one can define c 1 > 0 and c 2 > 0 by (dependences are omitted)

c 1 = 1 2 √ 2πǫ p 1 - √ 2π IE[b d,j ] 2 d -2π 2 ǫ 2 r 2 (49) c 2 (p 2 ) = inf qγ sup c {c : L γ (c) ≤ p 2 } (50)
where function f is defined by (77) and

L γ (c) = √ 2 α =γ Z α Z γ (c q α 1 ǫ + f (q α , ǫ r )) (51) If n d ≥ 8 c 1 (p 1 ) 2 log 4 1 -P 1 ( 52
)
then the following probabilistic inequality holds:

P ∀k ′ ∈ D HR , G γ (k ′ , n d ) r 2 -1 ≤ p 1 ≥ P 1 (53)
If

n d ≥              1 c 2 ( √ 2p) 2 log 2 1 -P 1 3 2 if r = 2, 1 c 2 (p) 2 log   4 1 -P 1 r 2 -1 2   if r ≥ 3. (54)
then the following concentration inequality holds:

P ∀k ′ ∈ D HR , B(k ′ , n d ) Z(k ′ ) ≤ p 2 ≥ P 2 (55)
Let us comment on Theorem 3. In ideal experimental conditions, with no positioning bias and ǫr ≤ p/2π 2 ,

c 1 (p) ≃ p -2π 2 ǫ 2 r 2 2 √ 2πǫ (56)
The quantity c 2 (p 2 ) can be computed numerically for some given specific image. A necessary condition to the existence of c 2 (p 2 ) is

p 2 > p 0 (ǫ, r, Z) = √ 2 sup qγ α =γ Z α Z γ f (q α , ǫ r ) (57)
In the most favourable case when H = 1 (no blur) and the image has a power law Fourier spectrum ∝ k ′ -2(1+η)

2

, (44) permits to estimate p 0 (ǫ, r, Z). Then c 2 (p 2 ) can be computed from (43) which is easy to use and gives quantitative indications about n d .

Corollary 4: Under the assumptions of Theorem 3 and denoting c 1 = c 1 (p 1 ) and c 2 = c 2 (p 2 ), if a sufficient number n d of images per position is used, one has the following concentration inequality which guarantees a small relative error with high probability:

P ∀k ′ ∈ D HR , X(k ′ ) -Z(k ′ ) Z(k ′ ) ≤ p 1 + p 2 ≥ P 2 -(1 -P 1 ) ≥ 1 -4e -c 2 2 n d r 2 -1 -4e -c 2 1 n d 8 ( 58
)
Proof : this is a direct consequence of Lemma 1 p. 4 applied to the sum of the approximation term |G γ /r 2 -1| and the aliasing term |B/ Z|. Corollary 4 gives a probabilistic bound to the total relative error on each frequency component of the reconstructed blurred image Z using the algorithm from [START_REF] Elad | A fast super-resolution reconstruction algorithm for pure translational motion and common space-invariant blur[END_REF] before the deconvolution step. Note that the bound in probability in (58) tends to 1 exponentially fast when n d → ∞. In practice, one can guarantee a global relative error ≤ 10% with probability ≥ 0.90 by choosing (p i , P i ) = (0.05, 0.95), i=1,2. This result provides a precise quantitative analysis of the reconstruction error. One limitation of the present study is that Z(k ′ ) = H(k ′ )Y HR (k ′ ) is the blurred super-resolved image resulting from the fusion of LR images. However the deconvolution step is common to every acquisition system and remains a limitation of any SR approach. Of course, the most favourable situation is when H(k ′ ) is close to 1, corresponding to a Dirac PSF in the spatial domain. Then Corollary 4 gives a good indication of the quality of high resolution imaging by using multiple acquisitions per positions.

In summary, we propose a detailed analysis of the reconstruction error of a fast method in the Fourier domain. It provides an a priori estimate of the number of images/position necessary to guarantee a given quality of reconstruction of each frequency (Fourier mode) with high probability. Based on Monte Carlo simulations, it also allows to estimate a posteriori a map of confidence levels in the frequency domain. Section IV will show numerically that these bounds are tight. We have worked on the intermediate reconstructed image Z before the deconvolution step that is common to most SR methods. Theorem 3 can be used based on the generic assumption of a power-law spectrum that is usual for natural images or more specifically for one specific image.

D. What about the SNR ?

We have demonstrated theoretical bounds to control the quality of the super resolved image in the Fourier domain. However this result deals with each frequency separately. Now we aim at identifying the dependence of the SNR between the reconstructed image and the ground truth. Again, this SNR deals with Z not Y HR and it measures the quality of the fusion step and does not consider the posterior deconvolution effects.

We consider the mean square error :

X -Z 2 2 = k X(k ′ ) -Z(k ′ ) Z(k ′ ) 2 α 2 k ′ •| Z(k ′ )| 2 (59)
and compare it to the energy of the original HR image. The | Z(k ′ )| are considered as fixed (the ground truth) while the α k ′ are random variables here (relative error estimates). Now we show that IEα 2 k ′ is of the order of 1/n d for all k ′ so that SN R ∝ log n d . From (58) in Corollary 4,

P (|α k ′ | ≤ p) ≥ 1 -4r 2 e -c 2 n d ( 60
)
where c 2 (p) = min(c 2 2 (p), c 2 1 (p)/8). Note from ( 43) & (56) that the typical order of magnitude of c 1 (p) and c 2 (p) is p/ǫ so that we can consider that there exists λ > 0 such that

c 2 ≥ λp 2 /ǫ 2 . Then IE[α 2 k ′ ] ≤ |α k ′ |≤p α 2 k ′ p(α k ′ )dα k ′ + |α k ′ |≥p α 2 k ′ p(α k ′ )dα k ′ ≤ p 2 + 2 ∞ n=1 (n+1)p np α 2 k ′ p(α k ′ )dα k ′ ≤ p 2 + 4r 2 ∞ n=1 e -λn 2 p 2 n d /ǫ 2 (n + 1) 2 p 2 ≤ p 2 (1 + e -λp 2 n d /ǫ 2 K(n d )) (61)
where K(n d ) is finite, decreasing with n d and independent of k ′ . Choosing p 2 = 1/n d , one obtains for all k ′ ∈ D HR ,

IE[α 2 k ′ ] ≤ 1 n d (1 + e -λ/ǫ 2 K(n d )) (62)
and consequently taking the expectation of (59),

IE X -Z 2 2 = k ′ IE[α 2 k ′ ]| Z(k ′ )| 2 ≤ 1 n d (1 + e -λ/ǫ 2 K(n d )) Z 2 2 (63)
Finally, using Parceval's equality we get:

SNR( X, Z) ≥ 10 log 10 n d + K ( 64
)
where K is a constant depending on the energy of the original image. As a function of the number of images per position n d , the SNR is improved with a magnitude of 10dB/decade. We can compare this result with the weak Cramer-Rao lower bound on the reconstruction error T weak ∝ 1/ √ K + 1 where K + 1 is the number of images in [START_REF] Robinson | Statistical performance analysis of superresolution[END_REF] : at best, the SNR can grow as log(number of images) as predicted by (64). This indicates that the proposed method is efficient at the best expected level [START_REF] Ng | Analysis of displacement errors in high-resolution image reconstruction with multisensors[END_REF]- [START_REF] Tsai | Multiframe image restoration and registration[END_REF]. Fig. 4 shows SNR computed for high frequencies only (the reconstructed HR part of the spectrum). Results were computed from 100 Monte-Carlo simulations with uniform distribution of position errors with ǫ = 0.01 for 11 images (Lena, Barbara, Boat...). This global indication that the SNR is ∝ log n d completes previous detailed bounds for each Fourier component.

IV. NUMERICAL RESULTS

To illuminate the complex interplay between the many parameters involved, we study the problem from various viewpoints. Section IV-A studies the lower-bound on the number n d of images per position to guarantee a given maximum errror level. Section IV-B compares our theoretical results to numerical estimates of probabilities from Monte-Carlo simulations. Section IV-C studies the connection between results in the Fourier domain and in the spatial domain. Section IV-D shows how the presence of noise and the nature of the blur operator influence the results. Monte-Carlo simulations use 100 realizations of the acquisition procedure assuming a uniform distribution of position errors in (-ǫ, ǫ). When no image is specified, the power law spectrum assumption is used.

A. How many images to ensure some maximum error level ? Fig. 5 shows the dependence of the required number of images n d (k ′ ) on the frequency k ′ to guarantee that aliasing contribution is less than p 1 = 0.05 with probability P 1 ≥ 0.95 when r = 4 and ǫ = 0.001 for an image with a power law spectrum. As expected, the recovery of high frequencies requires more LR images. The results are nearly independent of the size N of images for N ≥ 32. A similar picture (not shown) stands for the approximation term. In general (not always) the control of aliasing effects is the most constraining.

Tab. II gathers the constraints for various values of r and ǫ for images with a power law spectrum (η = 0 here). Numbers are computed from ( 52) & (54) in Theorem 3 for parameters (p i , P i ) = (0.05,0.95), i = 1, 2. This choice of equidistribution of error is certainly not optimal but of practical use with respect to Corollay 4 garantying an error level ≤ p 1 +p 2 = 0.10 with probability larger than P 2 -(1 -P 1 ) = 0.90. The larger r, the larger the need for multiple images. The smaller the positioning uncertainty ǫ, the smaller the lower bound on n d .

As an example, we consider a setting where the sensor has LR pixels of width ≃ 1 µm. The random bias on the positioning system can be reasonably expected to be between 1 and 10 nm corresponding to ǫ ≃ 0.001 -0.01 LR pixel. The acquition rate of images is usually of the order of 10 im./s (e.g. in a DSLR). In practice, r 2 displacements are used so that a minimum acquisition time of about r 2 × n d × 0.1s is necessary. For r = 2 and ǫ = 0.01, relative errors ≤ 5% on the restored image can be guaranteed with probability ≥ 0.95 C. How are the errors localized in the spatial domain ?

We use Monte-Carlo simulations to study the localization of errors in the spatial domain. By selecting the less reliable frequency components of an image where the aliasing error is ≥ 10% with probability ≥ 0.1, one can reconstruct the corresponding spatial counterpart to localize and quantify their contribution. For r = 2, ǫ = 0.01, these unreliable components weight for a SNR of -27.2dB. The present theoretical analysis permits such a selection of frequencies as well. For a given number n d of images/position one can reconstruct the spatial counterpart of the less reliable frequencies where the aliasing error is expected to be ≥ 10% with probability ≥ 0.1 according to Theorem 3. Fig. 8 shows such a picture for Barbara for r = 2, ǫ = 0.01 and n d = 256 im./pos., to be compared with the minimum n d = 157 in Tab. II. As expected, the spoiled regions are the most textured ones as well as some contours (better seen on screen). Remember that the analysis focused on the modulus of Fourier spectra while phases carry the location information. Maximum gray levels are about 4 and the standard deviation is of 0.67 (to compare with 255 in 8 bits). These "non reliable" components then weight for a SNR of -25.8 dB w.r.t. superresolved frequencies only. At least on this example, our theoretical predictions both qualitatively and quantitatively agree with Monte Carlo results. Our analysis not only gives indications to choose n d but also produces a detailed map of the error distribution both in the Fourier domain and in the spatial domain.

D. How do noise and PSF influence performances ?

The influence of noise and PSF are two important questions. The problem of noise is not the most critical: averaging numerous images attenuates additive noise. The present approach considers additive contributions of numerous images affected by independent realizations of noise: this naturally tends to increase the SNR. This is easily checked experimentally and not illustrated here for sake of briefness. If the noise in LR images was too strong to be compensated by simple averaging, the utility of SR would be questionable since the main concern would first be to access reliably denoised information at low resolution, giving up hopes for high resolution. Here we assume that LR images are of sufficient quality. The question of the PSF is a much bigger concern since it is involved in the error analysis. Of course frequencies where H(k ′ ) = 0 are lost and we already mentionned that the present analysis is not valid for these frequencies. Moreover the structure of aliasing is influenced by the PSF in an important manner, see (42). All the experiments above considered the ideal situation of a Dirac PSF where H(k ′ ) = 1 ∀k ′ . The lines 'PSF(0.5)' in Tab. II show how the lower bounds of n d are modified in presence of a Gaussian PSF of width 0.5. As expected it dramatically influences the estimates, e.g. for (r, ǫ) = (2, 0.001) as the bound becomes 18 in place of 1. The control of the PSF is a real stake in the design of a SR system: the present study permits to quantitatively evaluate its influence.

V. CONCLUSION

We have presented a theoretical analysis of a cheap and fast SR technique which takes benefit from any accurate controlled positioning system, e.g. piezoelectric actuators for sensor shifting, now currently available on many optical systems. Such an approach comes with some constraints. It requires a static scene captured using a static camera in good lighting conditions to avoid a high level of noise. It may also suffer from a lack of depth of field or an inhomogeneity of translations between images due to parallax for instance. However the statistical analysis of the algorithm proposed in [START_REF] Elad | A fast super-resolution reconstruction algorithm for pure translational motion and common space-invariant blur[END_REF] produces error confidence intervals as a function of the number of available images. This is made possible by the simplicity of the algorithm itself and by exploiting the averaging effect of LR images taken at positions that are randomly distributed around the same reference position. This approach is cheap and realistic to enhance the resolution of many devices. Even not state of the art, theoretical guarantees are a strong advantage of the approach when the reliability of the restored image is at stake, e.g. in scientific imaging (biology, astronomy...). This analysis considers a zero-mean noise which gets attenuated in the HR image reconstruction by fusing many images. The resulting probabilistic upper bounds are a good complement to the Cramer-Rao lower bounds in [START_REF] Robinson | Statistical performance analysis of superresolution[END_REF] and are nearly tight since the order of magnitudes are similar. Numerical experiments illustrate our results in both the Fourier and spatial domains as well as the effect of the PSF. A strong aspect of this work is in its predictions for practical implementation. Such results also give precious hints on the design of SR systems. Future works may investigate similar probabilistic bounds for more sophisticated SR algorithms where some reconstruction priors are used [START_REF] Farsiu | Fast and robust multiframe super resolution[END_REF]- [START_REF] Yang | Image super-resolution via sparse representation[END_REF]. [START_REF] Foroosh | Extension of phase correlation to subpixel registration[END_REF], see figure 9(a). [START_REF] Robinson | Optimal registration of aliased images using variable projection with applications to super-resolution[END_REF], see fig. 9(b) where the grey lozenge represents the region |x

APPENDIX A. Proofs of Lemma 1 & 2 Proof of Lemma 1: x 2 1 + x 2 2 ≥ a 2 1 + a 2 2 ⇒ |x 1 | ≥ a 1 or |x 2 | ≥ a 2 proves
|x 1 | + |x 2 | ≥ a 1 + a 2 ⇒ |x 1 | ≥ a 1 or |x 2 | ≥ a 2 proves
1 | + |x 2 | ≤ a 1 + a 2 . Proof of Lemma 2: ∀i, |x i | ≤ a i ⇒ i |x i | ≤ i a i so that P ( i |x i | ≤ a i ) ≥ P ({∀i, |x i | ≤ a i }). Since the x i are independent, P ({∀i, |x i | ≤ a i }) = i P (|x i | ≤ a i ). Noting that ∀i, P (|x i | ≤ a i ) ≥ (1 -P i ) concludes the proof. QED.

B. Proof of concentration inequality (20)

As far as the real part of B G in ( 16) is concerned:

Re (B G ) = 1 r 2 n d d,j cos (q γ .b dj ) -IE[cos (q γ .b dj )] (65)
Since |q α • b dj | ≤ q α 1 ǫ r , we apply Hoeffding's inequality to (78) and (81) for δ α = c q α 1 ǫ, and for α-γ / ∈ {0, r/2} 2 :

P (|Re (G α (k ′ ))| ≥ δ ′ α ) ≤ 2e - 2r 2 n 2 d c 2 4 d,j sin 2 θ αd (83) P (|Im (G α (k ′ ))| ≥ δ ′ α ) ≤ 2e - 2r 2 n 2 d c 2 4 d,j cos 2 θ αd (84)
where δ ′ α = c q α 1 ǫ + f (q α , ǫ r ). Using (97) & (98) in App. E, Lemma 1 yields inequalities (30) & (31).

D. Computing c 2 (p) in (36)

Here we estimate c 2 (p) in (36) under assumptions of Theorem 3. If one neglects the effect of blur, we aim at computing the maximum value of c 2 (p) such that for all q γ , a c 2 (p) + p 0 ≤ p.

(85)

after little reorganization of (34) where we use

p 0 ≃ √ 2 2 α =γ |Y (q α )| |Y (q γ )| q α 2 1 ǫ 2 r (86) a = √ 2 α =γ |Y (q α )| |Y (q γ )| q α 1 ǫ. (87)
as soon as ǫ r ≪ 1 so that cubic terms can be neglected. We first study (86). We focus on the highest frequencies only, typically q γ = (π -2π rN , π -2π rN ). As a consequence, note that q γ 1+η 2 ≃ ( √ 2π) 1+η . Then, one needs to detail: (89)

For r = 2, computations are easy and only 2 terms both equal to 1 appear in F (r, N ). For r ≥ 3, one can observe that q α 1 ∼ q α 2 (norms are equivalent) so that when η = 0 one expects that F (r, N ) ∝ (r 2 -1), the number of terms in β =(0,0) . This is due to the fact that q α 1 α =γ ≃ π for large r. As a result, one obtains in good approximation that : q α 1 ǫ (91)

p 0 ≃ b 0 √ 2 η π 2 ǫ 2 r 2 (r 2 -1) (90
Using that q α 1 ∼ q α 2 (within constant factors), one expects that when η = 0, a ∝ 2 1+η/2 (r 2 -1)ǫ (92)

Numerical estimates for values 2 ≤ r ≤ 8 show that a = a 0 × 2 1+η/2 ǫ(r 2 -1)

where a 0 varies with η around a typical value of 1.3 for η = 0, e.g. a 0 ≃ 0.95 if η = -0.2 and a 0 ≃ 1.85 if η = 0.2 for all r ≥ 3. For r = 2, one finds a 0 ≃ 0.63, resp. 1.14 and 3.04 when η = -0.2, resp. 0 and 0.2.

E. Properties of complex roots of unity

θ αd = 2π r (α -γ)d = 2π r δd (94)
where α and γ are integers in (0, r -1) 2 . The set of the θ αd matches the set of products of complex roots of unity so that:

F. Expectations IE[G α ]

Taking the expectation of ( 13) with respect to b dj yields: One remarks that

IEG α (k ′ ) =
β jℓ = χ((α 2 -α 1 )N ) if j = ℓ, χ(k ′ α1 )χ(-k ′ α2 ) if j = ℓ, (108)
IE[G α1 (k ′ )G * α2 (k ′ )] = δ α1α2 r 2 n d χ((α 2 -α 1 )N ) -χ(k ′ α1 )χ(-k ′ α2 ) = δ α1α2 r 2 n d 1 -|χ(k ′ α1 )| 2 (111)
so that the G αi , α i = γ, are uncorrelated. QED.

Fig. 1 .

 1 Fig. 1. (a) Spatial domain: black disks and thick grid are the original LR sampling grid, the thin grid is the target HR grid. Other symbols are positions of 3 translated LR images of 1/2 LR pixel (r = 2); (b) Fourier domain: the inner (red) square contains LR frequencies (-N/2, N/2) 2, the outer square is for HR frequencies (-rN/2, rN/2) 2 . Arrows represent aliasing, see[START_REF] Champagnat | Statistical performance modeling for superresolution: a discrete data-continuous reconstruction framework[END_REF].

Fig. 2 .

 2 Fig. 2. (a) Barbara, (b) Results from Algorithm [1] with r = 2, ǫ = 0.1, n d = 32 im./pos. ; HF reconstruction error (zoom on screen), (c) n d = 32 ⇒ SNR HF =25.0dB, (d) n d = 1 ⇒ SNR HF =10.0dB.

  and then turn to the contribution of the aliasing term B(k ′ ). The reader interested in our main results only can directly move to sections III-C & III-D. Proofs are in Appendices B & C.

c 2 (

 2 p) = inf qγ sup c {c : q γ obeys ineq.(34)} (36) If c 2 (p) > 0 is well defined, then there exists a minimum number of images per position n d such that

Fig. 3 .

 3 Fig. 3. Pairs of parameters (p, ǫ) for which SR with guaranteed error bounds is feasible are at the bottom right of the curve for each SR factor r indicated on the right margin, see (44).

Fig. 4 .

 4 Fig. 4. SNR for high frequencies only is proportional to log 10 n d . Results from 100 Monte-Carlo simulations with uniform distribution of positions with ǫ = 0.01 for 11 images (Lena, Barbara, Boat...).

Fig. 5 .

 5 Fig. 5. Minimum number n d (k ′ ) to guarantee an aliasing error ≤ 5% with probability ≥ 0.95 for all k ′ ; ǫ = 0.001, r = 4.

Fig. 8 .

 8 Fig. 8. (l.) Barbara, (r.) contribution of the less reliable frequencies.

v rN -2β/r 2 1 v rN -2β/r 1+η 2 F 2 - 1 )

 1221 (r,N )(88) where v rN = (1 -2/rN, 1 -2/rN ). The sum F (r, N ) can be computed numerically. It weakly depends on N so thatF (r, N ) ≃ if r = 2, 1.2(r 2 -1)if r ≥ 3.

) where b 0 = 2 / 3

 23 if r = 2 or b 0 ≃ 1.2 if r ≥ 3.Now let study coefficient (87) along the same lines.

2 2 ee i2πδ 1 -

 221 (θ αd ) = r 2 if αγ ∈ {0, r/2} 2 , r/2 otherwise. (97) d sin 2 (θ αd ) = 0 if αγ ∈ {0, r/2} 2 , r/2 otherwise. (98)Properties (95) and (96) come from the observation that d∈(0,r-1) iθ αd = in the r.h.s. is zero since α = γ and for any integer 1 ≤ δ ≤ r -1, e i2πδ/r = 0 (100)Now we prove (97) and (98). To this aim we need:cos 2 (θ αd ) = 1 + cos(2θ αd ) 2 (101) sin 2 (θ αd ) = 1 -cos(2θ αd ) 2 αd . For 0 ≤ δ ≤ r -1, r if δ ∈ {0, r/2}, 1 -e i4πδ1 -e i4πδ/r = 0 otherwise, (103) so that using (99) againd e i2θ αd = r 2 if αγ ∈ {0, r/2} 2 , 0 otherwise. (104)Taking the real part yields d cos(2θ αd ). The sum of (101) & (102) over d ∈ (0, r -1) 2 yield (97) & (98).

de

  -i 2π r (α-γ)•d IE e -i 2π rN k ′ α •b dj (105) Then let χ(k ′ ) = IE e -i 2π rN k ′ •b dj the characteristicfunction of the distribution of b dj . It results from properties of roots of unity above that d e -i 2π r (α-γ)d = 0 when α = γ, r 2 when α = γ (106)so that denoting Kronecker's symbol by δ γα :IEG α (k ′ ) = δ γα χ(k ′ ) (107) G. The G α are uncorrelatedThe correlation between G α1 and G α2 for α i = γ is:IE[G α1 G * α2 ] 2π rN (α1-γ)dN e +i 2π rN (α2-γ)d ′ N × n d j,ℓ=1 IE e +i 2π rN (k ′ +α1N )•b dj e -i 2π rN (k ′ +α2N )•b dℓ β jℓ

  α1 )χ(-k ′ α2 ) = n d χ((α 2α 1 )N ) -χ(k ′ α1 )χ(-k ′ α2 ) (110)As a consequence one finally gets:

  HR , Y k and n k are rearranged in lexicographic ordered vectors. The least squares optimization problem can be formulated as:

Y HR is the (desired) high resolution image to estimate from the K LR images {Y k , 1 ≤ k ≤ K}. We assume the unknown HR image Y HR is a periodic bandlimited image sampled above the Nyquist rate. Each image Y k is a LR observation of the same underlying scene translated by F k . The blur matrices H k model the point spread function (PSF) of the acquisition system and matrices D k are the decimation operator by a factor r. If Y HR is of size r 2 N 2 ×1 and Y k of size N 2 × 1, matrices F k and H k are of size (rN ) 2 × (rN ) 2 while D k are N 2 × (rN ) 2 ; n k is the noise, generally assumed to be Gaussian white noise so that E(n k n t k ) = σ 2 I. Images Y

  2 , the outer square is for HR frequencies (-rN/2, rN/2) 2 . Arrows represent aliasing, see[START_REF] Champagnat | Statistical performance modeling for superresolution: a discrete data-continuous reconstruction framework[END_REF].

	r	super-resolution factor, typically r = 2, 3...
	Y HR Y k	High Resolution (HR) image ∈ R (rN ) 2 Low Resolution (LR) image ∈ R N 2
	n k	noise in low resolution image Y k
	F k	translation operator on HR images
	H k D k	convolution blur operator on HR images decimation operator N 2 × (rN ) 2
	d	target position of one LR image
	Y

d HR image Y translated by d = (dx, dy) b dj error on displacement de(j) = d + b dj n d number of LR images around position d ǫ maximum positioning error in LR pixel units η exponent of the spectrum of natural images k, k ′ spatial frequency vectors resp. at LR and HR X

Note that one should first check that every term in the products are positive to ensure that the inequality above be relevant, which will be guaranteed by the final criterion.