{'fraction_non_alphanumeric': 0.04706437332172146, 'fraction_numerical': 0.011593729588504246, 'mean_word_length': 3.981652205350687, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 24, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Geometrical Study of Matching Pursuit Parametrization

This paper studies the effect of discretizing the parametrization of a dictionary used for Matching Pursuit decompositions of signals. Our approach relies on viewing the continuously parametrized dictionary as an embedded manifold in the signal space on which the tools of differential (Riemannian) geometry can be applied. The main contribution of this paper is twofold. First, we prove that if a discrete dictionary reaches a minimal density criterion, then the corresponding discrete MP (dMP) is equivalent in terms of convergence to a weakened hypothetical continuous MP. Interestingly, the corresponding weakness factor depends on a density measure of the discrete dictionary. Second, we show that the insertion of a simple geometric gradient ascent optimization on the atom dMP selection maintains the previous comparison but with a weakness factor at least two times closer to unity than without optimization. Finally, we present numerical experiments confirming our theoretical predictions for decomposition of signals and images on regular discretizations of dictionary parametrizations.


Introduction
There has been a large effort in the last decade to develop analysis techniques that decompose non-stationary signals into elementary components, called atoms, that characterize their salient features [1][2][3][4][5]. In particular, the matching pursuit (MP) algorithm has been extensively studied [2,[6][7][8][9][10][11] to expand a signal over a redundant dictionary of elementary atoms, based on a greedy process that selects the elementary function that best matches the residual signal at each iteration. Hence, MP progressively isolates the structures of the signal that are coherent with respect to the chosen dictionary, and provides an adaptive signal representation in which the more significant coefficients are first extracted. The progressive nature of MP is a key issue for adaptive and scalable communication applications [12,13].
A majority of works that have considered MP for practical signal approximation and compression define the dictionary based on the discretization of a parametrized prototype function, typically a scaled/modulated Gaussian function or its second derivative [6,14,15]. An orthogonal 1-D or 2-D wavelet basis is also a trivial example of such a discretization even if in that case MP is not required to find signal coefficients; a simple wavelet decomposition is computationally more efficient. Works that do not directly rely on a prototype function either approximate such a parametrized dictionary based on computationally efficient cascades of filters [16][17][18], or attempt to adapt a set of parametrized dictionary elements to a set of training signal samples based on vector quantization techniques [19,20]. Thus, most earlier works define their dictionary by discretizing, directly or indirectly, the parameters of a prototype function.
The key question is then: how should the continuous parameter space be discretized ? A fine discretization results in a large dictionary which approximates signals efficiently with few atoms, but costs both in terms of computational complexity and atom index entropy coding. Previous works have studied this trade-off empirically [6,15]. In contrast, our paper focuses on this question in a formal way. It provides a first attempt to quantify analytically how the MP convergence is affected by the discretization of the continuous space of dictionary function parameters.
Our compass to reach this objective is the natural geometry of the continuous dictionary. This dictionary can be seen as a parametric (Riemannian) manifold on which the tools of differential geometry can be applied. This geometrical approach, of increasing interest in the signal processing literature, is inspired by the works [21,22] on Image Appearance Manifolds, and is also closely linked to manifolds of parametric probability density function associated to the Fisher information metric [23]. Some preliminary hints were also provided in a Riemannian study of generalized correlation of signals with probing functions [24].
The outcome of our study is twofold. On the one hand, we analyze how the rate of convergence of the continuous MP (cMP) is affected by the discretization of the prototype function parameters. We demonstrate that the MP using that discretized dictionary (dMP) converges like a weak continuous MP, i.e. a MP algorithm where the coefficient of the selected atom at each iteration overtakes only a percentage (the weakness factor ) of the largest atom magnitude. We describe then how this weakness factor decreases as the so-called density radius 1 of the discretization increases. This observation is demonstrated experimentally on images and randomly generated 1-D signals.
On the other hand, to improve the rate of convergence of discrete MP without resorting to a finer but computationally heavier discretization, we propose to exploit a geometric gradient ascent method. This allows to converge to a set of locally optimal continuous parameters, starting from the best set of parameters identified by a coarse but computationally light discrete MP. Each atom of the MP expansion is then defined in two steps. The first step selects the discrete set of parameters that maximizes the inner product between the corresponding dictionary function and the residual signal. The second step implements a (manifold 2 ) gradient ascent method to compute the prototype function parameters that maximize the inner product function over the continuous parameter space. As a main analytical result, we demonstrate that this geometrically optimized discrete MP (gMP) is again equivalent to a continuous MP, but with a weakness factor that is two times closer to unity than for the non-optimized dMP. Our experiments confirm that the proposed gradient ascent procedure significantly increases the rate of convergence of MP, compared to the non-optimized discrete MP. At an equivalent convergence rate, the optimization allows reduction of the discretization density by an order of magnitude, resulting in significant computational gains.
The paper is organized as follows. In Section 2, we introduce the notions of parametric dictionary in the context of signal decomposition in an abstract Hilbert space. This dictionary is then envisioned as a Hilbert manifold, and we describe how its geometrical structure influences its parametrization using the tools of differential geometry. Section 3 surveys the definition of (weak) continuous MP providing a theoretical optimal rate of convergence for further comparisons with other greedy decompositions. A "discretization autopsy" of this algorithm is performed in Section 4 and a resulting theorem explaining the dependences of the dMP convergence relatively to this sampling is proved. A simple but illustrative example of a 1-D dictionary, the wavelet (affine) dictionary, is then given. The optimization scheme announced above is developed in Section 5. After a review of gradient ascent optimization evolving on manifolds, the geometrically optimized MP is introduced and its theoretical rate of convergence analyzed in a second theorem. Finally, in Section 6, experiments are performed for 1-D and 2-D signal decompositions using dMP and gMP on various regular discretizations of dictionary parametrizations. We provide links to previous related works in Section 7 and conclude with possible extensions in Section 8.

Dictionary, Parametrization and Differential Geometry
Our object of interest throughout this paper is a general real "signal", i.e. a real function f taking value on a measure space X. More precisely, we assume f in the set of finite energy signals, i.e. f ∈ L 2 (X, dµ) = {u : X → R : u 2 = X |u(x)| 2 dµ(x) < ∞}, for a certain integral measure dµ(x). Of course, the natural comparison of two functions u and v in L 2 (X, dµ) is realized through the scalar product This very general framework can be specialized to 1-D signal or image decomposition where X is given respectively by R or R 2 , but also to more special spaces like the two dimensional sphere S 2 [25] or the hyperboloid [26]. In the sequel, we will write simply L 2 (X) = L 2 (X, dµ).
In the following sections, we will decompose f over a highly redundant parametric dictionary of real atoms. These are obtained from smooth transformations of a real mother function g ∈ L 2 (X) of unit norm. Formally, each atom is a function g λ (x) = [U (λ)g](x) ∈ L 2 (X), for a certain isometric operator U parametrized by elements λ ∈ Λ and such that g λ = g = 1. The parametrization set Λ is a continuous space where each λ ∈ Λ corresponds to P continuous components λ = {λ i } 0≤i≤P −1 of different nature. For instance, in the case of 1-D signal or image analysis, g may be transformed by translation, modulation, rotation, or (anisotropic) dilation operations, each associated to one component λ i of λ. Our dictionary is then the set dict(g, U, Θ) g λ (x) = [U (λ)g](x) : λ ∈ Θ , for a certain subset Θ ⊆ Λ. In the rest of the paper, we write dict(Θ) = dict(g, U, Θ), assuming g and U implicitly given by the context. For the case Θ = Λ, we write D = dict(Λ).
We assume that g is twice differentiable over X and that the functions g λ (x) are twice differentiable on each of the P components of λ. In the following, we write ∂ i for the partial derivative with respect to λ i , i.e. ∂ ∂λ i , of any element (e.g. g λ (x), g λ , u , . . . ) depending on λ, and ∂ ij = ∂ i ∂ j . From the smoothness of U and g, we have ∂ ij = ∂ ji on quantities built from these two ingredients.
Let us now analyze the geometrical structure of Λ. Rather than an artificial Euclidean distance we use a distance introduced by the dictionary D itself seen as a P -dimensional parametric submanifold of L 2 (X) (or a Hilbert manifold 4 [27]). The dictionary distance d D is thus the distance in the embedding space From this embedding, we can define an intrinsic distance in D, namely the geodesic distance. This later has been used in a similar context in the work of Grimes and Donoho [22] and we follow here their approach. For our two points λ a , λ b , assume that we have a smooth curve γ : [0, 1] → Λ with γ(t) = γ 0 (t), · · · , γ P −1 (t) , such that γ(0) = λ a and γ(1) = λ b . The length L(γ) of this curve in D is thus given by L(γ) The geodesic distance between λ a and λ b in Λ is the length of shortest path between these two points, i.e.
where γ(λ a → λ b ) is any differentiable curve γ(t) linking λ a to λ b for t equals to 0 and 1 respectively. We denote by γ λaλ b the optimal geodesic curve joining λ a and λ b on the manifold D, i.e. such that L(γ λaλ b ) = d G (λ a , λ b ), and we assume henceforth that it is always possible to define this curve between two points of Λ. Note that by construction, In the language of differential geometry, the parameter space Λ is a Riemannian manifold M = (Λ, G ij ) with metric G ij (λ) = ∂ i g λ , ∂ j g λ . Indeed, for any differentiable curve γ : t ∈ [−δ, δ] → γ(t) ∈ Λ with δ > 0 and γ(0) = λ, we have withu(t) = d dt u(t), and where Einstein's summation convention is used for simplicity 6 . The vector ξ i =γ i (0) is by definition a vector in the tangent space T λ Λ of Λ in λ. The meaning of relation (2) is that the metric G ij (λ) allows the definitions of a scalar product and a norm in each T λ Λ. The norm of a vector ξ ∈ T λ Λ is therefore noted |ξ| 2 = |ξ| 2 λ ξ i ξ j G ij (λ), with the correspondence d dt g γ(t) | t=0 = |γ|. For the consistency of further Riemannian geometry developments, we assume that our dictionary D is non-degenerate, i.e. that it induces a positive definite metric G ij . Appendix A provides additional details.
We conclude this section with the arc length (or curvilinear ) parametrization "s" [28] of a curve γ(s). It is such that . From its definition, the curvilinear parameter s is the one which measures at each point γ(s) the length of the segment of curve already travelled on γ from γ(0). Therefore, in this parametrization, λ a = γ λaλ b (0) and λ b = γ λaλ b (d G (λ a , λ b )).

Matching Pursuit in Continuous Dictionary
Let us assume that we want to decompose a function f ∈ L 2 (X) into simpler elements (atoms) coming from a dictionary dict(Θ), given a possibly uncountable and infinite subset Θ ⊆ Λ. Our general aim is thus to find a set of coefficients {c m } such that f (x) is equal or well approximated by f app (x) = m c m g λm (x) with a finite set of atoms {g λm } ⊂ dict(Θ).
Formally, for a given weakness factor α ∈ (0, 1], a General Weak(α) Matching Pursuit decomposition of f [2,29], written MP(Θ, α), in the dictionary dict(Θ) is performed through the following greedy 7 algorithm : The quantity R m+1 f is the residual of f at iteration m + 1. Since it is orthogonal to atom g λ m+1 , Notice that the selection rule (3c) concerns the square of the real scalar product g λ , R m f . Matching Pursuit atom selection is typically defined over the absolute value | g λ , R m f |. However, we prefer this equivalent quadratic formulation first to avoid the abrupt behavior of the absolute value when the scalar product crosses zero, and second for consistency with the quadratic optimization framework to be explained in Section 5. Finally, to allow the non-weak case where α = 1, we assume that a maximizer g u ∈ dict(Θ) of g, u 2 always exists for any u ∈ L 2 (X).
If Θ is uncountable, our general Matching Pursuit algorithm is named continuous Matching pursuit. In particular, for Θ = Λ, we write cMP(α) = MP(Λ, α). The rate of convergence (or convergence) of the cMP(α), characterized by the rate of decay of R m f with m, can be assessed in certain particular cases. For instance, if there exists a Hilbert space S ⊆ L 2 (X) containing D = dict(Λ) such that then the cMP(α) converges inside S. In fact, the convergence is exponential [30] We name β = β(S, D) the greedy factor since it charaterizes the MP convergence (greediness). The existence of the greedy factor β is obvious for instance for finite dimensional space [30], i.e. f ∈ C N , with finite dictionary (finite number of atoms).
For a finite dictionary in an infinite dimensional space, as L 2 (X), the existence of β is not guaranteed over the whole space. However, there exists on the space of functions given by linear combination of dictionary elements, the number of terms being restricted by the dictionary (cumulative) coherence [29].
In the case of an infinite dictionary in an infinite dimension space where the greedy factor vanishes, cMP(α) convergence is characterized differently on the subspace of linear combination of countable subsets of dictionary elements. This question is addressed separately in a companion Technical Report [31] to this article. We now consider only the case where a non-zero greedy factor exists to characterize the rate of convergence of MP using continuous and discrete dictionaries.

Discretization effects of Continuous Dictionary
The greedy algorithm cMP(α) using the dictionary D is obviously numerically unachievable because of the intrinsic continuity of its main ingredient, namely the parameter space Λ. Any computer implementation needs at least to discretize the parametrization of the dictionary, more or less densely, leading to a countable set Λ d ⊂ Λ. This new parameter space leads naturally to the definition of a countable subdictionary D d = dict(Λ d ). Henceforth, elements of Λ d are labelled with roman letters, e.g. k, to distinguish them from the continuous greek-labelized elements of Λ, e.g λ.
For a weakness factor α ∈ (0, 1], the discrete Weak(α d ) Matching Pursuit algorithm, or dMP(α), of a function f ∈ L 2 (X) over D d is naturally defined as dMP(α) = MP(Λ d , α). The replacement of Λ by Λ d in the MP algorithm (3) leads obviously to the following question that we address in the next section.

Discretization Autopsy
By working with D d instead of D, the atoms selected at each iteration of dMP(α) are of course less optimal than those available in the continuous framework. Answering Question 1 requires a quantitative measure of the induced loss in the MP coefficients. More concretely, defining the score function S u (λ) = g λ , u 2 for some u ∈ L 2 (X), we must analyze the difference between a maximum of S u computed over Λ and that obtained from Λ d . This function u will be next identified with the residue of dMP(α) at any iteration to characterize the global change in convergence.
We propose to found our analysis on the geometric tools described in Section 2.
An umbilical point for which ∂ i S u (λ a ) = 0 for all i, is obviously critical in any direction. An umbilical point corresponds geometrically either to maxima, minima or saddlepoints of S u relatively to Λ. where Proof. Let us define the twice differentiable function ψ(s) using the Cauchy-Schwarz (CS) inequality in L 2 (X) in the last equation. The result follows from the fact that d ds g γ(s) = 1. The previous Lemma is particularly important since it bounds the loss in coefficient value when we decide to choose S u (λ b ) instead of the optimal S u (λ a ) in function of the geodesic distance d G (λ a , λ b ) between the two parameters. To obtain a more satisfactory control of this difference, we need however a new property of the dictionary.
We start by defining the principal curvature in the point λ ∈ Λ as where γ ξ is the unique geodesic in M starting from λ = γ ξ (0) and with γ ′ ξ (0) = ξ, for a direction ξ of unit norm in T λ Λ.
If K does not exist (not bounded K λ ), by extension, D is said to be of zero condition number.
The notion of condition number has been introduced by Niyogi et al. [32] to bound the local curvature of an embedded manifold 8 in its ambient space, and to characterize its self-avoidance. Essentially, it is the inverse of the maximum radius of a sphere that, when placed tangent to the manifold at any point, intersects the manifold only at that point [33,34]. Our quantity K −1 is then by construction a similar notion for the dictionary D seen as a manifold in L 2 (X). However, it does not actually prevent manifold self-crossing on large distance due to the locality of our differential analysis 9 .
where G ij = G ij (λ) is the inverse 10 of G ij . 8 In their work, the condition number, named there τ −1 , of a manifold M ′ measures the maximal "thickness" τ of the normal bundle, the union of all the orthogonal complement of every tangent plane at every point of the manifold.
9 A careful study of local self-avoidance of well-conditioned dictionary would have to be considered but this is beyond the scope of this paper. 10 Using Einstein convention, this means The proof is given in Appendix B since it uses some elements of differential geometry not essential in the core of this paper. The interested reader will find also there a slightly lower bound than the bound presented in (8), exploiting covariant derivatives, Laplace-Beltrami operator and scalar curvature of M [28]. We can state now the following corollary of Proposition 1.

Corollary 1. In the conditions of Proposition 1, if D has a non-zero condition number
Therefore, in the dMP(α) decomposition of f based on D d , even if at each iteration the exact position of the continuous optimal atom of D is not known, we are now able to estimate the convergence rate of this MP provided we introduce a new quantity characterizing the set Λ d .
We say that Λ d covers Λ with a radius ρ d .
This radius characterizes the density of Λ d inside Λ. Given any λ in Λ, one is guaranteed that there exists an element k of Λ d close to λ, i.e. within a geodesic distance ρ d .
Theorem 1. Given a Hilbert space S ⊆ L 2 (X) with a non zero greedy factor β, and a dictionary then, for functions belonging to S, a dMP(α) algorithm using D d = dict(Λ d ) is bounded by the exponential convergence rate of a cMP(α ′ ) using D with a weakness parameter given by Proof. Notice first that since f ∈ S and D d ⊂ D ⊂ S, R m f ∈ S for all iteration m of dMP. Let us take the (m + 1) th step of dMP(α) and write u = R m f . We have of course Denote by gλ the atom of D that best represents R m f , i.e. S u (λ) = sup λ∈Λ S u (λ). Ifk is the closest element ofλ in Λ d , we have d G (λ,k) ≤ ρ d from the covering property of Λ d , and the Proposition 1 tells us that, which is the exponential convergence rate of the Weak(α) Matching Pursuit in D when β exists [29,30].
The previous proposition has an interesting interpretation : a weak Matching Pursuit decomposition in a discrete dictionary corresponds, in terms of rate of convergence, to a weaker Matching Pursuit in the continuous dictionary from which the discrete one is extracted.
About the hypotheses of the proposition, notice first that the existence of a greedy factor inside S concerns the continuous dictionary D and not the discrete one D d . Consequently, this condition is certainly easier to fulfill from the high redundancy of D. Second, the density requirement, ρ d < β/ √ 1 + K, is just sufficient since the Proposition 1 does not state that it achieves the best bound for the control of |S u (λ a ) − S u (λ b )| when λ a is critical. It is interesting to note that this inequality relates ρ d , a quantity that characterizes the discretization Λ d in Λ, to β and K, which depend only on the dictionary. In particular, β represents the density of D inside S ⊂ L 2 (X), and K depends on the shape of the atoms through the curvature of the dictionary.
Finally note that as β < 1 (from definition (4)) and K > 1 (Prop. 2), the density radius must at least satisfy ρ d < 1 √ 2 to guarantee that our analysis is valid.

A Simple Example of Discretization
Let us work on the line with L 2 (X) = L 2 (R, dt), and check if the hypothesis of the previous theorem can be assessed in the simple case of an affine (wavelet-like) dictionary.
We select a symmetric and real mother function g ∈ L 2 (R) well localized around the origin, e.g. a Gaussian or a Mexican Hat, normalized such that g = 1. The parameter set Λ is related to the affine group, the group of translations and dilations G aff . We identify λ = (λ 0 = b, λ 1 = a), where b ∈ R and a > 0 are the translation and the dilation parameters respectively. The dictionary D is defined from the atoms g Our atoms are nothing but the wavelets of a Continuous Wavelet Transform if g is admissible [35], and U is actually the representation of the affine group on L 2 (R) [36].
In the technical report [31], we prove that the associated metric is given by G ij (λ) = a −2 W , where W is a constant 2 × 2 diagonal matrix depending only of the mother function g and its first and second derivatives. Since G ij (λ) = a 2 W −1 , K can be bounded by a constant also associated to g and its first and second order time derivatives.
Finally, given the τ -adic parameter discretization with τ > 1 and a 0 , b 0 > 0, the density radius ρ d of Λ d is shown to be bounded by ρ d ≤ Ca −1 0 b 0 + D ln τ , with C and D depending only of the norms of g and its first derivative. This bound has two interesting properties. First, as for the grid Λ d , it is invariant under the change (b 0 , a 0 ) → (2b 0 , 2a 0 ). Second, it is multiplied by 2 n if we realize a "zoom" of factor 2 n in our τ -adic grid, in other words, if (b 0 , τ ) → (2 n b 0 , τ 2 n ). By the same argument, the true density radius has also to respect these rules. Therefore, we conjecture that Unfortunately, even for this simple affine dictionary, the existence of β = β(S, D) is non trivial to prove. However, if the greedy factor exists, the control of τ , a 0 and b 0 over ρ d tells us that it is possible to satisfy the density requirement for convenient values of these parameters.

Optimization of Discrete Matching Pursuits
The previous section has shown that under a few assumptions a dMP is equivalent, in terms of rate of convergence, to a weaker cMP in the continuous dictionary from which the discrete one has been sampled.
Question 2. Can we improve the rate of convergence of a dMP, not with an obvious increasing of the dictionary sampling, but by taking advantage of the dictionary geometry ?
Our approach is to introduce an optimization of the discrete dMP scheme. In short, at each iteration, we propose to use the atoms of D d as the seeds of an iterative optimization, such as the basic gradient descent/ascent, respecting the geometry of the manifold M = (Λ, G ij ).
Under the same density hypothesis of Theorem 1, we show that in the worst case and if the number of optimization steps is large enough, an optimized discrete MP is again equivalent to a continuous dMP, but with a weakness factor two times closer to unity than for the non-optimized discrete MP.
In this section, we first introduce the basic gradient descent/ascent on a manifold. Next, we show how this optimization can be introduced in the Matching Pursuit scheme to defined the geometrically optimized MP (gMP). Finally, the rate of convergence of this method is analyzed.

Gradient Ascent on Riemannian Manifolds
Given a function u ∈ L 2 (X) and S u (λ) = g λ , u 2 , we wish to find the parameter that maximizes S u , i.e. λ * = arg max λ∈Λ S u (λ) (P.1) Equivalently, by introducing h u,λ = g λ , u g λ , we can decide to find λ * by the minimization If we are not afraid to get stuck on local maxima (P.1) or minima (P.2) of these two non-convex problems, we can solve them by using well known optimization techniques such as gradient descent/ascent, or Newton or Newton-Gauss optimizations. We present here a basic gradient ascent of the Problem (P.1) that respect the geometry of M = (Λ, G ij ) [37]. This method increases iteratively the value of S u by following a path in Λ, composed of geodesic segments, driven by the gradient of S u .
Given a sequence of step size t r > 0, the gradient ascent of S u starting from λ 0 ∈ Λ is defined by the following induction [38] : where ξ r (λ 0 ) = |∇S u (φ r (λ 0 ))| −1 ∇S u (φ r (λ 0 )) is the gradient direction obtained from the gradient ∇ i S u = G ij ∂ j S u , and γ(s, λ 0 , ξ 0 ) is the geodesic starting at λ 0 = γ(0, λ 0 , ξ 0 ) with the unit velocity ξ 0 = ∂ ∂s γ(0, λ 0 , ξ 0 ). Notice that ∇ i is the natural notion of gradient on a Riemannian manifold. Indeed, as for the Euclidean case, with Practically, in our gradient ascent, we use the linear first order approximation of γ, i.e. φ r+1 (λ) = φ r (λ) + t r ξ r (λ), valid for small value of t r (error in O(t 2 r )). This is actually an optimization method since ∂ i S u (φ r (λ)) ξ i r = |∂S u (φ r (λ))| > 0 and S u (φ r+1 (λ)) = S u (φ r (λ)) + t r |∂S u (φ r (λ))| + O(t 2 r ) ≥ S u (φ r (λ)), for a convenient step size t r > 0. At each step of this gradient ascent, the value t r is chosen so that S u is increased. This can be done for instance by a line search algorithm [39]. From the positive definiteness of G ij and G ij , a fixed point More sophisticated algorithms such as Newton or Newton-Gauss can be developed to solve the Problem (P.2) on a Riemannian manifolds [38,40] even if, unlike to the flat case, a direct definition of the Hessian does not exist on differentiable manifolds. However, we will not use them here as our aim is to prove that a dMP driven by the very basic optimization above provides already a better rate of convergence than the non-optimized dMP.

Optimized Discrete Matching Pursuit Algorithm
Let us optimize each step of a discrete MP using the gradient ascent of the previous section.

Definition
Given sequence of positive integers κ m and a weakness factor 0 < α ≤ 1, the geometrically optimized discrete matching pursuit (gMP(α)) is defined by Notice that the best atom g ν m+1 is selected in the set Φ m {g φκ m (k) : k ∈ Λ d } ⊂ D. Elements of Φ m are determined by applying the optimization function φ r : Λ d → Λ of our gradient ascent defined in (11) on elements of Λ d . In consequence, Φ m depends on R m f and is thus different at each iteration m.

Rate of convergence
The following theorem characterizes the rate of convergence of the optimized Matching Pursuit defined in (12).
Theorem 2. Given the notations and the conditions of Theorem 1, there exists a sequence of positive integers κ m such that, the gMP(α) decomposition of functions in S ⊂ L 2 (X) optimized κ m steps at each iteration m, is bounded by the same rate of convergence as a cMP(α ′′ ) using the corresponding continuous dictionary D with α ′′ = α(1 − 1 2 β −2 ρ d (1 + K)) 1/2 ≤ α. In other words, for α = 1, a gMP is equivalent to a cMP with a weakness factor two times closer to unity than the one reached by a dMP in the same conditions. Before proving this result, let us introduce some new lemmata.
Proof. Without loss of generality, assume that S u (λ a ) ≥ S u (λ b ). If this is not the case, we can switch the labels a and b. Let us define λ(θ) = γ(θd) with θ ∈ [0, 1] on the geodesic γ = γ λaλ b . We have λ a = λ(0) and λ b = λ(1). Using the Corollary 1, the two following inequalities hold : . In other words, the critical nature of λ a and λ b divides by two the bound on the decreasing of S u between them compared to the situation where only one of these points is critical.
Lemma 2. Given a function u ∈ L 2 (X), assume that S u (λ) has a global maximum at λ M , i.e. ∂ i S u (λ M ) = 0 for all i, and write T k = {φ r (k) : r ∈ N} the trajectory of the gradient ascent described in (11) starting from a point k ∈ Λ d . There exists a λ ′ ∈ T k that can be reached in a finite number of optimization steps, such that For the sake of clarity, the proof of this technical Lemma is placed in Appendix C. The main idea is to find a point in the trajectory T k that is closer to λ M than k, and that is also critical in the direction of λ M so that Lemma 1 can be applied. Let us now enter in the proof of the previous proposition.
Proof of Theorem 2. In our gMP(α) decomposition of a function f ∈ S ⊂ L 2 (X) defined before, given the iteration m + 1 where u = R m f is analyzed, denote byλ the parameter of the atom in D maximizing S u , i.e. S u (λ) = sup λ∈Λ S u (λ).
In Theorem 2, even if the sequence of optimization steps κ m is proved to exist, it is actually unknown. One practical way to overcome this problem is to observe how the ratio |∇Su| Su decreases at each optimization steps, and to stop the procedure once this value falls below a predefined threshold. This follows from the idea that the closer to a local maximum S u (φ r (k)) is, the smaller must be the optimization step. As it is often the case in optimization problems, an upper bound on the number of optimization steps can be fixed jointly to this threshold test.

Experiments
In this section, dMP and gMP decompositions of 1-D and 2-D signals are studied experimentally in different situations. These will imply different classes of signals and different discretization of parametrization of various densities.
Prior to these experiments, some remarks have to be made about dMP and gMP implementations. First, for both algorithms, as described in Equations (3) and (12), a full-search has to be performed in D d = dict(Λ d ) to compute all the squared scalar products S u of the current residue u = R m f , with atoms g k . We decide thus to reduce the computational complexity of this full-search with the help of the Fast Fourier Transform (FFT). One component (for 1-D signals) or two components (for 2-D signals) of the parametrization correspond indeed to a regular grid of atoms positions, which makes S u a discrete correlation relatively to these parameters. Moreover, as described in detail in [41,42], we apply the fast boundary renormalization of atoms, where atoms of D truncated by the limit of the signal remain valid atoms, i.e. of unit norm, and features that suddenly terminate at the signal boundary are correctly caught in the procedure. Notice that all our dMP and gMP experiments are performed in the non-weak case, i.e. α = 1.
Second, for the Gradient-Ascent optimization, we realize some simplifications to the initial formulation : the best discrete atom only is optimized at each MP iteration and κ m = κ > 0 for all m, with κ typically equal to 5 or 10. Even if these two restrictions are not optimal compared to the method described in the theoretical results, the gain of the optimization in the quality of signals reconstructions is already impressive. We also set all the step sizes to t r = χ > 0, with χ = 0.1 in all our experiments. Then, at each optimization step r, we adaptively decrease the step parameter t r by dividing it by 2 if the ascent condition is not met, i.e. if S u (φ r+1 (k)) < S u (φ r (k)). If after 10 divisions, the ascent condition still does not hold, the optimization process is simply stopped.
Finally, let us mention that our algorithms are written in MATLAB c and are consequently not truly optimized. The different computation times that we provide through this section allow us only to compare various schemes, as for dMP and gMP decomposition of the same signal. All of our experiments were realized on a Pentium 1.73 GHz with 1Gb of memory.

One Dimensional Analysis
This section analyzes the benefit obtained from gMP, and from an increased density of the discrete dictionary, when decomposing some specific classes of randomly generated 1-D signals. In our experiments, each 1-D signal is of unit norm and has N = 2 13 samples. Each signal consists of the sum of 100 random bursts, each burst being a rectangular or Gaussian window, depending on the class of the signal. The position and magnitude of each burst is selected at random, according to a uniform distribution. The duration of the rectangular window and the standard deviation of the Gaussian function are selected uniformly within the range [ 1 2 L, 3 4 L], for L = 2 8 . The mother function of the dictionary is the Mexican Hat function g(t) ∝ (1 − t 2 ) e −t 2 /2 . Its scale and translation parameters are sampled as defined in Section 4.2, following the τ -adic discretization Λ d = {(nb 0 τ j , a 0 τ j ) : j, n ∈ Z}, with a 0 = 1. We work in the non-weak case, i.e. α = 1, for dMP and gMP, and we set κ = 10 for gMP. Figures 1(a) and 1(b) analyze how the energy R m f 2 of the residual decreases with the number m of MP iterations for the random Gaussian and rectangular signals, respectively. Notice that only a small number of iterations are studied (twelve) since our analysis aims at analyzing the behaviour of dMP and gMP on one class of signals. However the current residual R m f belongs only approximately to the considered class on small m when not many atoms have been substracted to f = R 0 f . Results presented are averaged over 20 trials. In each graph, two distinct discretizations of the Mexican Hat parameters are considered to provide two discrete dictionaries, with one (b 0 = 1, log 2 τ = 0.25) being two times denser than the other (b 0 = 2, log 2 τ = 0.5), according to the behavior 11 of the density radius ρ d analyzed in Section 4.2. Both discrete and geometrically optimized MP are studied for each dictionary. We observe that gMP significantly outperforms dMP, and that an increased density of the dictionary also speeds up the MP convergence. By comparing Figure 1(a) and 1(b), we also observe that the residual energy decreases much faster for Gaussian signals than for rectangular ones, which unsurprisingly reveals that the Mexican Hat dictionary is better suited to represent Gaussian structures. Figures 1(c)-1(f) further analyze the impact of the discretization of the dictionary parameters on MP convergence. In these figures, we introduce the notion of normalized atom energy (NAE) to measure the convergence rate of a particular dictionary dealing with a specific class of signals at a specific MP iteration step. Formally, the NAE denotes the expected value of the best squared atom coefficient computed on a normalized signal when this one is randomly generated within a specific class of signals. Mathematically, NAE = E g λ * , u u 2 , where u is a sample signal of the class and the g λ * the associated best atom for a fixed greedy algorithm (dMP or gMP). We show the dependence of NAE on the discretization for the 1 st and 30 th iteration 12 for both rectangular and Gaussian signals. Results are averaged over 500 trials.
By considering the dMP and gMP curves in Figures 1(c)-1(f), we first observe that the NAE is significantly higher for gMP than for dMP, which confirms the advantage of using gradient ascent optimization to refine the parameters of the atoms extracted by dMP. Note that the NAE for a Gaussian random signal (Fig. 1(c)-1(d)) is nearly one order of magnitude higher than for a rectangular one (Fig.1(e)-1(f)). This confirms that the Mexican Hat dictionary better matches the Gaussian structures than the rectangular ones. In all cases, the NAE sharply decreases with the iteration index, which is not a surprise as the coherence between the signal and the dictionary decreases as MP expansion progresses.
To better understand the penalty induced by the discretization of the continuous dictionary, we now analyze how the rate of convergence for a particular class of signals behaves compared to the reference provided by a signal composed of a single Mexican Hat function. For that purpose, an additional curve, denoted dMP a , has been plotted in each graph. This curve is expected to provide an upper bound to the penalty induced by a sparser dictionary. Specifically, dMP a plots the energy captured during the 1 st step of the dMP expansion of a random (scale and position) Mexican Hat function, as a function of the discretization parameter log 2 τ . As the Mexican Hat is the generative function of the dictionary, the 1 st step of the MP expansion would capture the entire function energy if the entire continuous dictionary were used, but is particularly penalized by a discretization of the dictionary. In each graph of Figures 1(c)-1(f), to compare dMP a with dMP, the dMP a curve obtained with pure atoms (i.e. unit coefficients) is scaled to correspond to atoms whose energy is set to the NAE expected from the expansion of the corresponding class of signals with a continuous dictionary. In practice, the NAE expected with a continuous dictionary is estimated based on the NAE computed with gMP and the densest dictionary (log 2 τ = 0.25).  Residual energy as a function of the MP iteration. dMP (b 0 , log 2 τ ) and gMP (b 0 , log 2 τ ) refer to discrete and optimized MP, computed on a discretization Λ d = {(nb 0 τ j , a 0 τ j ) : j, n ∈ Z} of the continuous Mexican Hat dictionary. (c)-(f) Normalized atom energy (NAE) as a function of the log 2 τ discretization parameter. b 0 is set to one in all cases. dMP and gMP respectively refer to discrete and optimized MP. dMP a provides a lower bound to the decrease of NAE with log 2 τ , and is formally described in the text.
The approximation is reasonable as we observe that gMP saturates for small log 2 τ values, i.e. for large densities. We first observe that both the dMP and the dMP a curves nearly coincide in Figure 1(c). Hence, the MP expansion of a Gaussian signal is penalized as much as the one of a Mexican Hat function by a reduction of the discrete dictionary density. We then observe that the penalty induced by a reduction of density decreases as the coherence between signal and dictionary structures drops. This is for example the case when the signal to represent is intrinsically sharper than the dictionary structures ( Fig. 1(e)-1(f)), or because the coherent structures have been extracted during the initial MP steps (Fig. 1(d)). This last observation is of practical importance because it reveals that using a coarsely discretized dictionary incurs a greater penalty during the first few iterations of the MP expansion than during the subsequent ones. For compression applications, it might thus be advantageous to progressively decrease the density of the dictionary along the expansion process, the cost associated to the definition of the atom indices decreasing with the density of the dictionary 13 . Hence, it might be more efficientin a rate-distortion sense -to use a dense but expensive dictionary during the first MP iterations, so as to avoid penalizing the MP convergence rate, but a sparser and cheaper during subsequent steps, so as to save bits. We plan to investigate this question in details in a future publication.

Two Dimensional Analysis
This section analyzes experimentally the effect of discretizing a dictionary on the Matching Pursuit decomposition of images, i.e. with the Hilbert space L 2 (R 2 ).

Parametrization and Dictionary
We use the same dictionary as in [41]. Its mother function g is defined by a separable product of two 1-D behaviors : a Mexican Hat wavelet in the xdirection, and a Gaussian in the y-direction, i.e. g(x) = ( 4 3π ) 1/2 (1 − x 2 ) exp(− 1 2 |x| 2 ), where x = (x, y) ∈ R 2 and g = 1 [30]. Notice that g is infinitely differentiable.
The dictionary is defined by the translations, rotations, and anisotropic dilations of g. Mathematically, these transformations are represented by operators T b , R θ , and D a , respectively. These are given by = (a 1 , a 2 ), a 1 , a 2 ∈ R * + , while r θ is the usual 2 × 2 rotation matrix r θ and d a = diag(a 1 , a 2 ).
In other words, we have a parametrization of P = 5 dimensions and Λ = {λ = (λ 0 , . . . , Obviously, the dictionary D is complete in L 2 (R 2 ). Indeed, translations, rotations and isotropic dilations alone are already enough to constitute a wavelet basis of L 2 (X) since g is an admissible wavelet [35,43]. Finally, as requested in the previous section, from the smoothness of g and of the transformations U above, the atoms g λ of our dictionary D are twice differentiable on each component λ i .
Spatial Sampling For all our experiments, images are discretized on a Cartesian regular grid of pixels, i.e. an image f takes its values on the grid X = [0, N x ) × [0, N y ) ∩ Z 2 , with N x and N y the "x" and "y" sizes of the grid. We work in the continuous approximation, that is we assume  that the grid X is fine enough to guarantee that the scalar products ·, · and norms · are well estimated from their discrete counterparts. This holds of course for band-limited functions on L(R 2 ). In consequence, in order to respect this continuous approximations and to have dictionary atoms smaller than the image size, the mother function g of our dictionary D must be dilated in a particular range of scales so that g λ is essentially band-limited, i.e. a 1 , a 2 ∈ [a m , a M ]. According to the definition of g above, we set experimentally a m = 0.7 and a M = min(N x , N y ).
Discrete Parameter Space We decide to sample regularly Λ so that to have N pix = N x N y positions b, J 2 scales a 1 and a 2 selected logarithmically in the range [a m , a M ], and K orientations evenly spaced in [0, π) , with J, K ∈ N. At the end, we obtain the discretized parameter set Λ d = Λ d (N pix , J, K) = (b, θ n , a 1j , a 2j ′ ), b ∈ X , n ∈ [0, K − 1], j, j ′ ∈ [0, J − 1] , and the corresponding dictionary D d (N pix , J, K) = dict (Λ d (N pix , J, K)). The number of atoms in the dictionary is simply |D d | = J 2 K N pix .  Results We start our experiment by decomposing the venerable image of Barbara. 300 atoms were selected by dMP and gMP for various J and K. Results are presented in Table 1. In these tests, the best quality obtained for dMP corresponds obviously to the finest grid, i.e. J = 5 and K = 8 (26.63 dB, Fig.2(a)), with a computational time (CT) of 4634s. With 10 optimization steps (κ = 10), the gMP for the coarsest parametrization (J = 3 and K = 4) is equivalent to the best dMP result with a PSNR of 26.68 dB and a CT of only 950s, i.e. almost five time faster. This is also far better than the dMP on the same grid (24.30 dB). The visual inspection of the dMP image (J = 5, K = 8, Fig.2(a)) and the gMP image (J = 3, K = 4, κ = 10, 2(b)) is also instructive. Most of the features of the gMP results are well represented (e.g. Barbara's mouth, eyes, nose, hair, ...). However, the regular pattern of the chair in the background of the picture, which needs a lot of similar atoms, is poorly drawn. This can be explained by the fact that this highly directional structure has to be represented by a lot of similarly oriented and scaled atoms with similar amplitude. The fine grid of dMP has therefore more chance to correctly fit these atoms, while the gMP on its coarse grid is deviated in its optimization process to more prominent structure with higher amplitudes. Notice finally, the best optimized result (PSNR 28.09 dB) is obtained for κ = 10 on the grid associated to J = 5 and K = 8 orientations.
For our second experiment, we compare dMP and gMP (κ = 10) 300 atoms approximation of well known 128×128 pixels pictures, namely Lena, Baboon, Cameraman, GoldHill, and Peppers, on the same parametrization grid (J = 4, K = 8). For a computational time slightly higher (5%) than the dMP decomposition, we reach in all cases a significantly higher PSNR with gMP than with dMP, i.e. the dB gain is within the range [0.87, 2.03].

Related Works
A similar approach to our geometric analysis of MP atom selection rule has been proposed in [24]. In that paper, a dictionary of (L 2 -normalized) wavelets is seen as a manifold associate to a Riemannian metric. However, the authors restrict their work to wavelet parametrization inherited from Lie group (such as the affine group). They also work only on the L 2 (dictionary) distance between dictionary atoms and do not introduce intrinsic geodesic distance. They define a discretization of the parametrization Λ such that, in our notations, G ij ∆λ i ∆λ j < ǫ, with ∆λ(k) the local width of the cell localized on k ∈ Λ d . There is however no analysis of the effect of this discretization on the MP rate of convergence.
In [14], the author uses a 4-dimensional Gaussian chirp dictionary to analyze 1-D signals with MP algorithm. He develops a fast procedure to find the best atom of this dictionary in the representation of the current MP residual by applying a two-step search. First, by setting the chirp rate parameter to zero, the best common Gabor atom is found with full search procedure taking advantage of the FFT algorithm. Next, a ridge theorem proves that starting from this Gabor atom, the best Gaussian chirp atom can be approximated with a controlled error. The whole method is similar to the development of our optimized matching pursuit since we start also from a discrete parametrization to find a better atom in the continuous one. However, our approach is more general since we are not restricted to a specific dictionary. We use the intrinsic geometry of any smooth dictionary manifold to perform a optimization driven by a geometric gradient ascent.

Conclusions
In this paper, we have adopted a geometrical framework to study the effect of dictionary discretization on the rate of convergence associated to MP. In a first step, we have derived an upper bound for this rate using geometrical quantities inherited from the dictionary seen as a manifold, such as the geodesic distance, the condition number of the dictionary, and the covering property of the discrete set of atoms in the continuous dictionary. We have also shown in a second step how a simple optimization of the parameters selected by the discrete dictionary, can lead theoretically and experimentally to important gain in the approximation of (general) signals.
In a future study, it could be interesting to see how our methods extend to other greedy algorithms, like the Orthogonal Matching Pursuit (OMP) [44]. However, this extension has to be performed carefully since we need to characterized the convergence of continuous OMP, as it is here for the one of MP induced by the existence of a greedy factor.
Our work paves the way for future extensions and advances in two practical fields. As explained in our 1-D experiments, a first idea could be to analyze carefully the benefit -in a rate-distortion sense -of using a dense but expensive dictionary during the first MP iterations, so as to avoid penalizing the MP convergence rate, but a sparser and cheaper dictionary during subsequent steps, so as to save bits. We plan to investigate this question in details in a future publication.
Another idea is to analyze the behaviors of gMP in the Compressive Sensing (CS) formalism, that is after random projection of the signal and atoms. Matching Pursuit is already used currently as a retrieval algorithm of CS of sparse signals [45][46][47]. However, recent results [48] suggests also that for manifold of bounded condition number, their geometrical structure (metric, distances) is essentially preserved after random projection of their points in a smaller space than the ambient one. If a natural definition of random projection in our continuous formalism can be formulated, a natural question is thus to check if the gradient ascent technique survives after random projection of the residual and the atoms on the same subspace. This could lead to dramatic computation time reduction, up to controlled errors that could be even attenuated by the greedy iterative procedure.

Acknowledgements
LJ wishes to thank Prof. Richard Baraniuk and his team at Rice University (Houston, TX, USA) for the helpful discussions about general "manifolds processing" and Compressive Sensing. LJ is also very grateful to R. Baraniuk for having accepted and funded him during a short postdoctoral stay at Rice University. LJ and CDV are funded by the Belgian FRS-FNRS. We would like to thank Dr. David Kenric Hammond (LTS2/EPFL, Switzerland) for his careful proofreading and the referees for valuable comments on this paper.
A Complements on the Geometry of (Λ, G ij ) In this short appendix, we provide some additional information on the geometrical concepts developed in Section 2. First, as explained in that section, the parameter space Λ of the dictionary D = dict(Λ) is linked to a Riemannian manifold M = (Λ, G ij ) with a structure inherited from the dictionary D ⊂ L 2 (X). From the geodesic definition (1) and the metric relation (2), we see that the curve γ λaλ b (t) ∈ Λ is thus also a geodesic in M. In other words, it is defined only from the metric G ij and not anymore from the full behavior of atoms of D ⊂ L 2 (X). In [31], we explain also that M is in fact an immersed manifold [28] in the Hilbert manifold D ⊂ L 2 (X), and G ij is the associated pullback metric. All the geometric quantities of the Riemannian analysis of M, such as Christoffel's symbols, covariant derivatives, curvature tensors, etc. can be defined. This is actually done in the following appendices of this paper.
Second, some important designations can be introduced. The metric G ij (λ) is a (covariant) tensor of rank-2, i.e. described by two subscript indices, on M. This means that G ij satisfies a specific transformation under changes of coordinates in T λ Λ such that the values of the bilinear form 14 G λ (ξ, ζ) ξ i ζ j G ij (λ) that it induces are unmodified 15 . A function f : Λ → R is a scalar field on M, or rank-0 tensor. A vector field ζ i (λ) on this manifold, which associates to each point λ a vector in the tangent plane T λ Λ, is a function ζ : Λ → T λ Λ ≃ R P also named (contravariant) rank-1 tensor, i.e. with one superscript. More generally, a rank-(m, n) tensor is a quantity T i 1 ··· im j 1 ··· jn (λ) m-times contravariant and n-times covariant such that G i 1 k 1 · · · G imkm ξ k 1 1 · · · ξ km m T i 1 ··· im j 1 ··· jn (λ) ζ j 1 1 · · · ζ jn n is invariant under change of coordinates in T λ Λ for any vectors {ξ 1 , · · · , ξ m , ζ 1 , · · · , ζ n } in this space.
We need now some elements of differential geometry. Since γ is a geodesic in M, it respects the second order differential equation γ ′′ k + Γ k ij γ ′ i γ ′ j = 0, where the values Γ k ij = 1 2 G lk ∂ j G li + ∂ i G jl − ∂ l G ij are the Christoffel's symbols [28] derived from the metric G ij . Therefore, we get where ∇ i g γ = ∂ i g γ and ∇ ij g γ = ∇ i ∇ j g γ = ∂ ij g γ − ∂ k g γ Γ k ij are by definition the first order i and the second order ij covariant derivatives of g γ respectively [28]. In addition, we can easily compute that for M = (Λ, G ij ), Γ k ij = G kl ∂ ij g λ , ∂ l g λ .
Using (16) and the expression (17) of the Christoffel's symbols above, we have d 2 ds 2 g γ(s) 2 = ∇ ij g γ γ ′ i γ ′ j 2 ≤ ∇ ij g γ , ∇ kl g γ G ik G jl , where we used |γ ′ | = 1 and the Cauchy-Schwarz (CS) inequality expressed in the Einstein's summation notation on rank-2 tensors. This latter states that, for the tensors A ij = ∇ ij g γ and B ij = γ ′ i γ ′ j , |A ij B ij | 2 ≤ |A ij A kl G ki G lj | |B ij B kl G ki G lj |, the equality holding if the two tensors are multiple of each other. We prove in [31] the general explanation for rank-n tensor as a simple consequence of the positive-definiteness of G ij . Therefore, taking γ = γ ξ , and since γ ξ (0) = λ, In the companion Technical Paper [31], we prove that this inequality is also equivalent to where R is the scalar curvature of M, i.e. the quantity R = R ijkl G ik G jl contracted from the curvature tensor R iklm = 1 2 (∂ kl G im + ∂ im G kl − ∂ km G il − ∂ il G km ) + G np (Γ n kl Γ p im − Γ n km Γ p il ), and ∆g λ = G ij ∇ i ∇ j g λ is the Laplace-Beltrami operator applied on g λ . The curvature R requires only the knowledge of G ij (λ) (and its derivatives), implying just one step of scalar products computations, i.e. integrations in L 2 (X).
The reader who does not want to deal with differential geometry can however get rid of the covariant derivatives of Equation (18) by replacing them by usual derivatives. This provides however a weaker bound. Indeed, using the expression (17) of the Christoffel's symbols, some easy calculation provides 0 ≤ ∇ ij g λ , ∇ kl g λ G ik G jl = ∂ ij g λ , ∂ kl g λ G ik G jl − a ijk a lmn G il G jm G kn , with a ijk = ∂ ij g λ , ∂ k g λ .
(ii) If ζ r decreases monotically for r > 0 : Since ζ r ≥ 0, the limit lim r→∞ ζ r exists and converges to ζ ∞ < ζ 0 . However, it is not guaranteed that the sequence {φ r (k)} converges to a point of Λ. Fortunately, since for all r > 0, φ r (k) remains in the finite volume V 0 = {λ ∈ Λ : d M (λ) ≤ d M (k)}, this sequence is bounded in the finite dimensional space Λ. Therefore, from the Bolzano-Weierstrass theorem on the metric space (Λ, d G (·, ·)), we can find a convergent subsequence {r i ∈ N : r i+1 > r i } such that lim i→∞ φ r i (k) = k ∞ ∈ V 0 . On this point, we will have ∇ i S u (k ∞ ) = 0 for all i. So, k ∞ is an umbilical point and, from Lemma 1, From now on, we abuse notation and write φ r i (k) = φ i (k). Since ζ 2 ∞ = d G (λ M , k ∞ ) 2 < d G (λ M , k) 2 = ζ 2 0 , we can find a δ > 0 such that d G (λ M , k ∞ ) 2 + δ < d G (λ M , k) 2 . Therefore, because lim i→∞ S u (φ i (k)) = S u (k ∞ ) by continuity of S u , and since S u (φ i (k)) increases monotically with i, there exists a i ′ > 0 such that S u (k ∞ ) − S(φ i ′ (k)) ≤ 1 2 u 2 δ (1 + K). With . This gives the result and concludes the proof.