Blake Woodworth
email: blakewoodworth@gmail.com

Francis Bach
email: francis.bach@inria.fr

Alessandro Rudi
email: alessandro.rudi@inria.fr

Non-Convex Optimization with Certificates and Fast Rates Through Kernel Sums of Squares

We consider potentially non-convex optimization problems, for which optimal rates of approximation depend on the dimension of the parameter space and the smoothness of the function to be optimized. In this paper, we propose an algorithm that achieves close to optimal a priori computational guarantees, while also providing a posteriori certificates of optimality. Our general formulation builds on infinite-dimensional sums-of-squares and Fourier analysis, and is instantiated on the minimization of multivariate periodic functions.

Introduction

A well-designed optimization algorithm provides two important types of guarantees. First, it guarantees a priori that its output will achieve a certain degree of accuracy, with computational complexity that is hopefully adaptive to the specific properties of function to be optimized and possibly even optimal over a certain class of algorithms or functions to minimize. Second, it provides an a posteriori certificate, i.e., an explicit bound on the solution's accuracy that we can calculate once we run the algorithm. There are many examples of such well-designed optimization algorithms in the convex setting, which often use some form of convex duality (see, e.g., [START_REF] Nemirovski | Accuracy certificates for computational problems with convex structure[END_REF].

In this paper, our goal is to provide a well-designed algorithm for non-convex optimization,

c * = inf x∈X f (x), (1)
with X ⊆ R d and f potentially non-convex. In general, this task is extremely difficult, and in the worst case the computational cost must be exponential in the dimension d. However, it is known [START_REF] Novak | Deterministic and Stochastic Error Bounds in Numerical Analysis[END_REF] that in order for an algorithm to achieve the optimal computational complexity in solving Eq. ( 1), it must be adaptive to the degree of differentiability of f . That is, it should be able 1 to overcome the curse of dimensionality in terms of the approximation error ε, when the function is very smooth. More precisely, if f is m-times differentiable, then the computational complexity for finding c * with error ε should be C d ε -d/m , where C d is exponential in d in the worst case, but the dependence on the accuracy scales as only ε -d/m , which becomes quite mild once m approaches d. In this case, the curse of dimensionality is relegated just to the constant C d , making it possible to efficiently solve non-convex problems to high accuracy as long as d is relatively small, which has applications to tasks like hyperparameter tuning, industrial process optimization, and more.

Establishing well-designed algorithms for non-convex optimization is a difficult task. Many nonconvex optimization algorithms in the literature lack a priori guarantees, a posteriori guarantees, or actually any guarantees at all, and many methods used in practice are based on heuristics or can only guarantee convergence to a local rather than global minimum. Some other algorithms have good a posteriori guarantees, but weak a priori bounds; for example, methods based on polynomial sum of squares [START_REF] Lasserre | Global optimization with polynomials and the problem of moments[END_REF][START_REF] Parrilo | Semidefinite programming relaxations for semialgebraic problems[END_REF] are not adaptive, a priori, to the smoothness and are therefore subject to the curse of dimensionality in terms of the accuracy ε. The new family of algorithms based on kernel sum of squares [START_REF] Rudi | Finding global minima via kernel approximations[END_REF] achieves quasi-optimal a priori guarantees, but without any certificate a posteriori.

Our Contribution. In this paper, we provide a general strategy to derive algorithms that compute a lower bound ĉ of c * with strong guarantees both a priori and a posteriori (see Corollary 12). As a particular example, we consider optimizing smooth, periodic functions f on X = [0, 1] d and derive an algorithm that: (a priori) approximates c * with almost optimal error ε and computational complexity that scales well with ε, and (a posteriori) provide a certificate of accuracy, which is adaptive to the specific instance of the problem.

The a priori guarantee is useful since it shows that the proposed algorithm has nearly optimal complexity, which is adaptive to the smoothness of the function to be optimized. The a posteriori certificate is particularly useful because the accuracy of our estimate of c * depends on the specific instance at hand, and may be much better than the worst-case, exponential-in-d constant would suggest. Indeed, better-than-worst-case performance is frequently observed in practice, but we need a certificate of accuracy in order to know when we are so lucky.

2 Deriving Well-Designed Algorithms for Smooth Non-Convex Optimization

Our approach begins with rewriting the problem (1) as finding the highest lower bound on f :

c * = max c∈R c such that f (x) ≥ c ∀x ∈ X.
The inequality constraint f -c ≥ 0 can be converted to an equality constraint by introducing a non-negative function g ≥ 0:

c * = max c∈R, g≥0 c such that f (x) -c -g(x) = 0 ∀x ∈ X.
Finally, we can rewrite this again in a penalized form

c = max c∈R, g≥0 c -∥f -c -g∥ L ∞ (X) . (2)
It may not yet be obvious what this accomplishes, but the following Lemma indicates its promise:

Lemma 1. The problem (2) is a concave maximization problem with solution c = c * , and for any feasible (c, g) such that c ∈ R and g ≥ 0, c * is lower-bounded by c -∥f -c -g∥ L ∞ (X) .

Proof. The objective V (c, g)

:= c -∥f -c -g∥ L ∞ (X) is concave because it is a linear term c minus a convex term ∥ • ∥ L ∞ (X)
composed with a linear function of the optimization variables.

Since for all x ∈ X, f (x) -c -g(x) ≥ -∥f -c -g∥ L ∞ (X) , for any feasible (c, g), then ∀x ∈ X, f (x) ≥ V (c, g). Thus, by minimizing with respect to x, c * ≥ c. In order to show the other inequality, we notice that

(c * , f -c * ) is feasible and V (c * , f -c * ) = c * .
So, we have reduced the non-convex minimization problem (1) to a concave maximization problem (2), which is an improvement. However, whereas the original problem had a d-dimensional optimization variable, our new problem requires optimizing over the infinite-dimensional space of non-negative functions, and it is not clear how to do this. Furthermore, the quantity ∥f -c-g∥ L ∞ (X) is just as difficult to compute as inf x∈X f (x) would be. We will now describe several modifications leading to a more tractable optimization problem that maintains the same desirable properties.

A more tractable formulation. Specifically, our approach revolves around introducing a more tractable norm ∥•∥ W on functions from X to R, and restricting g ∈ G, for G a tractable subset of all non-negative functions. We will specify W and G later. To provide guarantees on the method we only need them to satisfy the following:

Definition 2 (Norms and models for non-negative functions).

1. Let ∥ • ∥ W be a norm on the space of real-valued functions over X, such that ∥ • ∥ L ∞ (X) ≤ ∥ • ∥ W . Denote by W the associated Banach space.

2. Let G be a convex subset of the set of non-negative functions, such that G is a closed subset of W. Now, by restricting the problem in Eq. ( 2) on G and considering the norm ∥ • ∥ W , we obtain what can be a more tractable formulation,

c = max c∈R, g∈G c -∥f -c -g∥ W . (3)
It is, of course, not the case that any choice of ∥•∥ W and G would make Eq. ( 3) easy to solve. However, we will later discuss examples where Eq. ( 3) is much easier to solve than Eq. ( 2). Regardless, we will see in the next Theorem, that this formulation comes with strong guarantees, expressing the error of the algorithm directly in terms of the approximation properties of the class of models G for non-negative functions.

Theorem 3 (Tightness of Eq. ( 3)). Suppose that f -c ∈ W for any c ∈ R, then

c * -q ≤ c ≤ c * , q = min g∈G ∥f -c * -g∥ W . Moreover, if there exists g * ∈ G satisfying ∥f -c * -g * ∥ W = 0, then c = c * . Finally, any c ∈ R and g ∈ G leads to a lower bound c -∥f -c -g∥ W on c * . Proof. By construction c ≤ c * , since c -∥f -c -g∥ W ≤ c -∥f -c -g∥ L ∞ (X) ≤ c * and G ⊆ {g | g ≥ 0}.
Now the problem is well defined, since it corresponds to a maximization on the closed subset R × G of a concave and continuous objective function (for the topology inherited from W). By setting c = c * and optimizing over g, the optimized objective is exactly c * -q with q as above, thus c * -q ≤ c. Moreover,we have q = 0, i.e., c = c * when there exists g * ∈ G satisfying ∥f -c * -g * ∥ W = 0.

In the Theorem above we see that c is lower bound of c * by construction. Moreover, the error c * -c is bounded by q, the approximation error of f -c * with respect to the class of models for non-negative functions G that we consider, measured in the norm W . So far, this all holds for any ∥•∥ W and G satisfying Definition 2; we continue by analyzing a specific choice.

The model for non-negative functions. We would like to use a class of models G that can approximate smooth non-negative functions using as few parameters as possible, while remaining tractable to optimize over. We consider the class of PSD models introduced by Marteau-Ferey et al.

(2020) and defined as g(x) = ϕ(x) * Aϕ(x), A ⪰ 0, for a suitable map ϕ : X → C n and A ∈ C n×n Hermitian positive semidefinite1 , where ϕ * denotes the conjugate transpose of ϕ. By the definition of positive semidefiniteness, g(x) ≥ 0 for any x ∈ X, and this is also a tractable class to optimize over since g is linear in the parameters A. The approximation properties of the model will depend on the choice of the feature map ϕ, and have been shown to give rates in the order of n -m/d for specific choices [START_REF] Rudi | PSD representations for effective probability models[END_REF].

Here, however, we need to study the approximation error with respect to our norm of choice W , and we will see that different feature maps than the ones considered by [START_REF] Rudi | Finding global minima via kernel approximations[END_REF] will lead to better rates.

The norms. We need norms that bound ∥ • ∥ L ∞ (X) from above as tightly as possible, but are also easy to compute. We consider from this viewpoint two norms:

1. The "F norm": L 1 norm of the Fourier transform, i.e.,

∥u∥ F = R d |û(ω)|dω,
2. The "S norm": The norm associated to a richer reproducing kernel Hilbert space such as the Sobolev space of exponent (d + 1)/2. Let S(ω) non negative and integrable, we define

∥u∥ 2 S = C 2 S R d |û(ω)| 2 S(ω) dω,
and C 2 S = R d S(ω)dω. For example for the Sobolev case we set S(ω) = (1 + ∥ω∥ 2 ) (d+1)/2 . Lemma 4. The norms above satisfy ∥

• ∥ L ∞ (R d ) ≤ ∥ • ∥ F ≤ ∥ • ∥ S ,
|u(x)| = û(ω)e 2πi⟨x, ω⟩ dω ≤ ∥û∥ L 1 (R d ) ∥e 2πi⟨x, •⟩ ∥ L ∞ (R d ) ≤ ∥û∥ L 1 (R d ) =: ∥u∥ F , ∀x ∈ R d from which we conclude that ∥•∥ F always upper bounds ∥•∥ L ∞ (R d )
. Analogously, note that, for any S ≥ 0 and S ∈ L 1 (R d ), and for any u such that ∥u∥ S is finite, we have, by Cauchy-Schwartz,

∥u∥ 2 F = ∥ u∥ 2 L 1 (R d ) = S 1/2 u S 1/2 2 L 1 (R d ) ≤ ∥S 1/2 ∥ 2 L 2 (R d ) u S 1/2 2 L 2 (R d ) = C 2 S |û(ω)| 2 S(ω) = ∥u∥ 2 S .
Finally, for u with finite F -norm, ∥u∥ F = ∥u∥ S when S = | u| for which S ≥ 0, S ∈ L 1 (R d ).

The F norm and the S norm are of interest when we can assume that we have access to the Fourier transform of the target function f . In particular the norm S can be also computed in closed form in specific scenarios. For example, consider the case when f is of the form

f = M j=1 β j h(x -x i ),
for some β 1 , . . . , β M ∈ R and x 1 , . . . , x M ∈ R d . This arises, e.g., for mixtures of Gaussians, and learning linear models or RBF networks. If we know the Fourier transform of h, then

∥f ∥ 2 S := M i,j=1 β i β j H(x i -x j ),
where H is the inverse Fourier transform of the function S(ω)| ĥ(ω)| 2 .

Perhaps more interestingly, Lemma 4 shows that the F norm is weaker than the norm S for any S, meaning that the F norm allows us to automatically adapt to certain structures in the function. For example, suppose that f (x) = h(P x) for some unknown

P ∈ R d ′ ×d with d ′ ≪ d and h : R d ′ → R.
Using the S norm for a certain S that depends on P would allow us to adapt to the low-dimensional structure and depend on d ′ rather than d, but this requires knowing P . On the other hand, the F norm is always weaker, so we can take advantage of the low-dimensional structure automatically.

PSD Models for Periodic Functions on [0, 1] d and Their Approximation Properties

The goal of the section is to provide a self-contained introduction to the approximation properties of PSD models. In particular, we consider the problem of approximating smooth 1-periodic functions (which corresponds to X = [0, 1] d ) using PSD models where ϕ is a subset of the Fourier basis. This setting, while already being of interest for practical applications, allows for an elementary proof which highlights the main conceptual steps of the derivation.

The main results of this section are Theorem 6 and, in particular, Theorem 8. With a more refined proof based on the same strategy, it is also possible to obtain results that hold for other scenarios beyond periodic functions on the torus and for more general maps ϕ. See, for example, [START_REF] Rudi | PSD representations for effective probability models[END_REF] for the approximation of non-periodic C m functions on subsets of R d via PSD models based on a finite dimensional feature map defined with respect to the Gaussian kernel, or [START_REF] Rudi | Finding global minima via kernel approximations[END_REF] for a feature map defined with respect to any kernel that satisfy some algebraic property, such as the Sobolev kernel.

We have seen in Theorem 3 that the optimization error of Eq. ( 3) depends on the approximation error of the function f with respect to the class of models for non-negative functions. So, in our analysis there will be three main ingredients: a class of models G t parametrized by its bandwidth t, which will depend on the space of functions associated with a feature map ϕ t ; the space of functions where f lives, which we denote H ρ ; and the norm that we use to measure the approximation error, in our case, ∥ • ∥ F .

We start by introducing G t , parametrized by a bandwidth t ∈ N. We associate each entry in

ϕ t (x) with an element of {k ∈ Z d : |k| ≤ t} where |k| = j |k j |, with n = #({k ∈ Z d : |k| ≤ t}) ≤ (2t+1) d , i.e., n = O(t d ). So for each |k| ≤ t, we define the feature map ϕ t : X → C n elementwise as (ϕ t (x)) k = e k (x)
, where e k is the k-th Fourier component, i.e., e k (x) = e -2πik ⊤ x . Consider the class of PSD models

G t = {g A,t | A ∈ C n×n , A ⪰ 0}, g A,t (x) = ϕ t (x) * Aϕ t (x). (4)
We are thus considering the feature map ϕ t associated with the classical band-limited space of functions. This choice is convenient for our analysis, but there are also many other choices of finite dimensional feature maps for PSD models that can have good approximation properties [START_REF] Rudi | PSD representations for effective probability models[END_REF][START_REF] Rudi | Finding global minima via kernel approximations[END_REF].

We consider continuous, 1-periodic functions f on R d , i.e., functions satisfying f (x + k) = f (x) for any x ∈ R d and k ∈ Z d . We note that these can therefore be identified with continuous, periodic functions on the torus [0, 1] d . We now introduce H ρ , where ρ ∈ ℓ 1 (Z d ) is a strictly positive summable sequence. The space is a separable Hilbert space of periodic functions defined as

H ρ = {f ∈ L 2 (X) | ∥f ∥ ρ < ∞}, where ∥f ∥ 2 ρ = k∈Z d | f k | 2 /ρ k < ∞ and f k = [0,1] d e k (x)f (x)dx is the Fourier series associated to f . One classical example of ρ is ρ k = (1 + ∥k∥ 2 ) -m , with m > d/2, corresponding to the Sobolev space H m
2,per of periodic functions whose derivatives up to order m are squared integrable [START_REF] Wahba | Spline Models for Observational Data[END_REF], or the space of periodic entire functions (of order 1), corresponding to ρ k = exp(-σ∥k∥), for some σ > 0.

We will soon present a Theorem showing that the PSD models G t can approximate functions of the form f = q j=1 u 2 j for q ∈ N and u j ∈ H ρ , using a small t depending on the decreasing quantity

R 2 t := |k|>t ρ k .
First, we start with a Lemma concerning the norm ∥ • ∥ F of the pointwise product of functions. This part of the proof is crucial and is handled differently in the other settings (e.g., [START_REF] Rudi | Finding global minima via kernel approximations[END_REF].

Lemma 5. Let f, g be 1-periodic functions on X with ∥f ∥ F , ∥g∥ F < ∞ and denote by f • g their pointwise product, i.e., (f

• g)(x) = f (x)g(x). Then ∥f • g∥ F ≤ ∥f ∥ F ∥g∥ F .
Proof. By the convolution property of Fourier series,

( f • g) k = j∈Z d f j g k-j
. By the Young inequality for the convolution of discrete sequences, we have for any two sequences u

, v ∈ ℓ 1 (Z d ), j,k∈Z d | u j v k-j | ≤ ( k∈Z d |u k |)( k∈Z d |v k |).
The result is obtained by applying this inequality on the Fourier series of f • g and noting that k∈Z d | f k | is exactly ∥f ∥ F , and the same for g. Now we are ready to state the first theorem that on the approximation error for the PSD models described above.

Theorem 6. Let f (x) = T j=1 u j (x) 2 for functions u j ∈ H ρ and T ∈ N, then min g∈Gt ∥f -g∥ F ≤ C ′ f R t .
where

C ′ f 2 = T j=1 ∥u j ∥ ρ ∥u j ∥ F .
Proof. Denote by u j,t the function, u j,t (x) = |k|≤t ( u j ) k e k (x) (a low-pass filtered version of u), and by v j,t the n-dimensional vector (v j,t ) k = u j,t , for any |k| ≤ t. Now, define Ā ∈ C n×n as

Ā = T j=1 v j,t v * j,t . Since, by construction, v * j,t ϕ(x) = |k|≤t ( u j ) k e k (x) = u j,t (x), then g Ā,t (x) = ϕ t (x) * Aϕ t (x) = T j=1 ϕ t (x) * (v j,t v * j,t )ϕ t (x) = T j=1 (v * j,t ϕ t (x)) 2 = T j=1 u j,t (x) 2 . Now note that u 2 j -u 2 j,t = (u j + u j,t ) • (u j -u j,t ), then, by using Lemma 5, ∥f -g Ā,t ∥ F = T j=1 (u 2 j -u 2 j,t ) F ≤ T j=1 (∥u j ∥ F + ∥u j,t ∥ F )∥u j -u j,t ∥ F .
We conclude noting that ∥u j,t ∥ F ≤ ∥u j ∥ F by construction and, by Cauchy-Schwartz,

∥u j -u j,t ∥ F = |k|>t | u j | = |k|>t √ ρ t | u j | √ ρ t ≤ R t ∥u j,t ∥ ρ ≤ R t ∥u j ∥ ρ .
Therefore,

min g∈Gt ∥f -g∥ F = min A∈C n×n ,A⪰0 ∥f -g A,t ∥ F ≤ ∥f -g Ā,t ∥ F ≤ R t C ′ f .
The theorem above controls the approximation error of the PSD models of bandwidth t when the target function can be written in terms of a sum of squares of functions belonging to an H ρ for a given ρ. In general it is not clear how to guarantee when a function f can be characterized as a sum of squares of functions in a given space. Luckily, in the case of m-times differentiable functions, there exists an easy geometrical characterization. We are going to use this fact, to specify the result above for the case when f is an m-times differentiable function. First, we need the following lemma, which is the adaptation to periodic functions of Theorem 2 of [START_REF] Rudi | Finding global minima via kernel approximations[END_REF] (more specifically, of Corollary 2, page 23) .

Lemma 7 [START_REF] Rudi | Finding global minima via kernel approximations[END_REF]). Let f be an m + 2-times differentiable non-negative periodic function. Assume that the minimizers of f in X are finitely many and with strictly positive Hessian.

Then, there exists

Q ∈ N and z 1 , . . . , z Q periodic m-times differentiable functions, such that f = Q j=1 z 2 j .
The proof of the lemma above is reported in Appendix A. Now we are ready to specify Theorem 6 in the case of an m-times differentiable function.

Theorem 8. Let f be an (m + d/2 + 2)-times differentiable periodic function, with m > 0 and let c * be its global minimum. Assume that the minimizers of f in X are finitely many and with strictly positive Hessian. Then, for any t ∈ N,

min g∈Gt ∥f -c * -g∥ F ≤ C f t -m ,
where the constant C f depends only on f, m, d.

The proof is self-contained and reported in Appendix B. It is obtained by first applying Lemma 7 on f , and Theorem 6 on the resulting characterization. To make this possible and to obtain a sharp rate, a crucial step is to show that the resulting functions belong to the space H ρ for a specific ρ satisfying ρ k ∝ |k| -2m-d , then deriving the bound on the associated residual R t .

The Resulting Problem and the Associated A Priori Guarantees

Now the problem Eq. ( 3), with the PSD models (4) and the F norm, takes the following form

c = max c∈R, A∈C n×n c -∥f -c -g A ∥ F such that A ⪰ 0, (5)
and, combining Theorem 3 and Theorem 8 gives the following a priori guarantee Corollary 9. Let f be an (m + d/2 + 2)-times differentiable, 1-periodic function with m > 0, and let c * be its global minimum. Also, let f have finitely many minimizers in X, which each have strictly positive Hessian. Then, for any t ∈ N

0 ≤ c * -c ≤ C f t -m .
Expressing t with respect to n, the dimension of the matrix A, we have t = O(n 1/d ). The bound above, then reads as

0 ≤ c * -c ≤ C ′ n -m/d .
This shows that the solution, c is always a lower bound of the global minimum c * and converges to c * with a rate depending on the dimension of the matrix A and the degree of differentiability of f . E.g., when m ≥ d, the error goes to zero as quick as n -1 . In the following section we see how to solve the optimization problem Eq. ( 5) in practice, by making use of the fact that ∥ • ∥ F , in the case of the torus, is a sum, which makes it easy to write Eq. ( 5) as a stochastic optimization objective.

Solving the Optimization Problem

We now describe the process of solving the optimization problem in Eq. ( 3) in the specific case of the F norm and a PSD model g A (x) = ϕ(x) * Aϕ(x) parametrized by positive semidefinite A ∈ C n×n . For now, we consider an arbitrary feature map ϕ, but we will also contextualize our results in the specific case of ϕ t , the map introduced in Section 2.1. A serious challenge to solving (5) is that computing ∥f -c -g A ∥ F or its subgradients exactly will typically be intractable because the F norm is the series:

∥f -c -g A ∥ F = k∈Z d | fk -c1 {k=0} -g A k |.
To circumvent this issue, we recast the problem as a stochastic optimization objective. In particular, we introduce a probability measure, π, supported on Z d and rewrite

∥f -c -g A ∥ F = k∈Z d π k | fk -c1 {k=0} -g A k | π k = E k∼π | fk -c1 {k=0} -g A k | π k .
Written this way, we can now attack our objective using any number of methods from the stochastic optimization arsenal, such as projected stochastic gradient ascent.

To see how π should be chosen, we first note that g A (x) = ⟨A, ϕ(x)ϕ(x) * ⟩, and use k) . Thus, our optimization problem now reads

M (k) = ϕϕ * k ∈ C n×n to denote the k-th Fourier component of ϕϕ * , so that g A k = A, M (
c = max c∈R, A∈C n×n c -E k∼π | fk -c1 {k=0} -A, M (k) | π k such that A ⪰ 0.
Noting that c only appears in two terms, we can also eliminate this variable by solving

max c c -| f0 -c -⟨A, M (0) ⟩| = f0 -⟨A, M (0) ⟩.
Putting this all together, we want to solve the stochastic concave maximization problem

c = max A⪰0 E k∼π [L k (A)] (6)
where

L k (A) =    1 π 0 f0 -A, M (0) k = 0 -1 π k fk -A, M (k) k ̸ = 0. (7) Algorithm 1 PROJECTED STOCHASTIC GRADIENT ASCENT Initialize A 0 = 0 for t = 0, 1, . . . , T -1 do Ãt+1 = A t + η∇L kt (A t ) for k t ∼ π A t+1 = min A ∥A -Ãt+1 ∥ F rob. s.t. A ⪰ 0, ∥A∥ F rob. ≤ R. Return: ĀT = 1 T T t=1 A t
Using projected stochastic gradient ascent yields the following error guarantee:

Theorem 10. Let R ≥ ∥A * ∥ Frob. be an upper bound the norm of a maximizing A * , and let ĀT be the output of Algorithm 1, with an optimally-chosen constant stepsize η and

π k ∝ ∥M (k) ∥ F rob. + (1 + d j=1 (2πk j ) d+1 ) -1 . Then ∥ ĀT ∥ Frob.
≤ R and for any δ ∈ (0, 1), with probability 1 -δ

E k∼π L k ( ĀT ) ≥ c - 20R log(2/δ) 1 + k∈Z d ∥M (k) ∥ F rob. √ T .
The proof, which we defer to Appendix C, simply requires proving that the functions L k are Lipschitz-continuous and then appealing to Proposition 2.2 from [START_REF] Nemirovski | Robust stochastic approximation approach to stochastic programming[END_REF]. This result bounds, a priori, the optimization error incurred in trying to estimate A * which realizes the maximum of (5). In Section 5, we combine this with Theorem 8 and the yet to be presented Theorem 11 to state our a priori guarantees.

A Posteriori Certification

Obviously, it is nice to know a priori that our estimate ĀT will be close to attaining the optimum of (5). However, with this estimate in hand, what we really want is to compute a lower bound on c * , so we need to actually evaluate E k∼π [L k ( ĀT )], which is non-trivial since π has infinite support.

Things are easier when f and ϕϕ * are band-limited, meaning that for some K, |k| > K implies fk = 0 and M (k) = 0. Specifically, we can choose π k ∝ ∥M (k) ∥ Frob. , which is only supported on {k : |k| ≤ K}, and then easily compute E k∼π [L k ( ĀT )] to obtain an exact lower bound on c * .

However, if one or both of f and ϕϕ * are not band-limited, then we are forced to estimate the value of an infinite sum. One approach is to draw samples k ∼ π and estimate the value using a sample average, and under suitable conditions on f and the matrices M (k) , this allows us to accurately estimate the value of the lower bound with high-probability. Alternatively, under stronger conditions on f and the matrices M (k) , we can compute L k for a finite set of k's and deterministically bound the contribution of the remaining, uncomputed terms. The following Theorem indicates the accuracy of these methods:

Theorem 11. Let f satisfy the conditions of Theorem 8 with m > d/2. Then for any K and k 1 , . . . , k K i.i.d.

∼ π, for any δ ∈ (0, 1) and A ⪰ 0, with probability 1 -δ,

c * ≥ c ≥ ĉ1-δ := 1 K K i=1 L k i (A) -Err 1-δ ≥ E k∼π [L k (A)] -2Err 1-δ Err 1-δ := ( √ d + 1∥f ∥ C d+1 (X) + ∥A∥ Frob. ) 2 log(2/δ) K 1 + k∈Z d ∥M (k) ∥ Frob.
,

where ∥f ∥ C d+1 (X) = max 1≤j≤d max 1≤q≤d+1 ∥ ∂ q ∂x q j f ∥ L ∞ (X) .
In addition, for any K, the following holds deterministically:

c * ≥ c ≥ ĉ1 := k:|k|≤K π k L k (A) -Err 1 ≥ E k∼π [L k (A)] -2Err 1 Err 1 := k:|k|>K | fk | + ∥A∥ Frob. ∥M (k) ∥ Frob.
The proof, which we defer to Appendix D, analyzes ĉ1-δ and ĉ1 separately. For the former, we first show that L k (A) is bounded for each k, and then apply Hoeffding's inequality. For the latter, we decompose the sum over k ∈ Z d into those k's with |k| ≤ K, and those k's with |k| > K, and then upper bound this second portion of the sum.

The Theorem shows that the sample average has additive error that decays with 1/ √ K with high probability. Furthermore, this lower bound is tractable given the parameter ζ and enough knowledge of our feature map for us to upper bound k∈Z d ∥M (k) ∥ Frob. . For the deterministic lower bound on c * , we need to have some control over how quickly ∥M (k) ∥ Frob. decays with increasing |k|, but if the feature map is chosen so that this decay is (eventually) rapid, then this lower bound can be tight. In the particular case of G t introduced in Section 2.1, we show in Appendix E that we can bound k∈Z d ∥M (k) ∥ Frob. ≤ n(8t) d , and for K ≥ 2t, k:|k|>K ∥M (k) ∥ Frob. = 0. Therefore, ĉ1-δ can provide a tight approximation of c using K ≫ n(8t) d samples, and ĉ1 can once the maximum bandwidth is set K ≥ 2t, as long as the Fourier coefficients fk decay sufficiently quickly.

With Theorem 11, we can use the solution returned by our optimization algorithm to compute a lower bound on c * , one that holds with high probability and one that holds deterministically. However, to actually compute a certificate of the accuracy of our lower bound, we also need an upper bound on c * . Getting some upper bound on c * is as easy as evaluating f (x) at any point x, although most x's will not be close to minimizing f , so this may not give us much information about c * . Of course, there are many better ways, and the bulk of the non-convex optimization literature is devoted to designing algorithms for computing approximate minimizers of f , i.e., upper bounds on c * . Upper bounds for c * are easier to produce, f (x 0 ) for any point x 0 ∈ X is a valid upper bound. We can use the point x 0 produced for example by [START_REF] Rudi | Finding global minima via kernel approximations[END_REF], that converges provably to a global minimizer with a rate that avoids the curse of dimensionality. In our experiments (in low dimensions), we simply compute f (x 1 ), . . . , f (x N ) for N random points and upper bound c * ≤ min i f (x i ), which allows for tight enough certificates.

A Priori and A Posteriori Guarantees

In the previous sections, we have described a method for estimating c * in the case of periodic functions on [0, 1] d , and all the pieces are in place to state our method's a priori and a posteriori guarantees. To summarize so far:

1. Theorem 3 shows that the solution of the relaxed problem (3), c, is a lower bound on c * which is tight up to the error of approximating f -c * with the class of non-negative functions, G. 2. In Theorem 8, we bound this approximation error for smooth, periodic functions with respect to G t , the class of PSD models defined in (4) with the band-limited kernel ϕ t . 3. But, we need to actually solve (3) defined using the F norm and G t . So, in Theorem 10, we bound the optimization error of the solution returned by projected stochastic gradient ascent. 4. However, our optimization algorithm returns the parameters ĀT of a PSD model, and to compute a lower bound on c * we need to actually evaluate the value of the objective at ĀT . So, finally, Theorem 11 bounds the estimation error when using ĀT to estimate lower bounds ĉ1-δ and ĉ1 on c * that holds with high probability and deterministically, respectively.

Therefore, our a priori guarantees amount to combining (Approximation Error) + (Optimization Error) + (Estimation Error). On the other hand, given any PSD model parameters, A, we can evaluate an a posteriori bound on the error by upper bounding c * ≤ f (x) for any x and lower bounding c * using Theorem 11. The following Corollary summarizes these guarantees:

Corollary 12. For the F norm and family of PSD models G t defined using ϕ t , under the conditions of Theorems 8, 10, and 11, let ĉ1-δ ( ĀT ) and ĉ1 ( ĀT ) be lower bound estimates defined in Theorem 11. Then for any δ ∈ (0, 1), we provide the following a priori guarantee with probability 1 -2δ:

c * ≥ ĉ1-δ ( ĀT ) ≥ c * -C f t -m -C d nt d 20R log(2/δ) √ T + 2( √ d + 1∥f ∥ C d+1 (X) + R) 2 log(2/δ) √ K ,
with C d = 8 d . At the same time, given any point x and parameters A for the PSD model, we guarantee a posteriori that f (x) ≥ c * ≥ ĉ1 (A) and f (x) ≥ c * ≥ ĉ1-δ (A) with probability 1 -δ.

The Corollary follows immediately by combining Theorems 8, 10, and 11. Since

t = O(n 1/d ), by choosing T = K = O(n 4+2m/d ), we have c * ≥ ĉ1-δ ( ĀT ) ≥ c * -C ′ n -m/d ,
when n is the dimension of the matrix ĀT . In this case the algorithm has a complexity that is

O(T n 3 + Kn 2 ) = O(n 7+2m/d ).
In particular, for the class of (m + d/2 + 2)-times differentiable functions with m > d/2, we achieve a bound c * ≥ ĉ1-δ ( ĀT ) ≥ c * -C ′ n -1 , with a computational cost of O(n 8 ). There is a lot of room for improvement in the constants of the exponents, but the considered algorithm shows that it is possible to obtain the global optimum of a function with both a posteriori guarantees and an a priori error rate that is adaptive to the degree of differentiability of the function to minimize and that avoids the curse of dimensionality for very smooth functions.

Empirical Evaluation

Finally, we apply our method to two simple non-convex optimization problems in one and two dimensions. The results are summarized in Figures 1 and2, and all of the details of the experiments are deferred to Appendix F, in which we describe a new feature map ϕ, and describe a more practical algorithm for solving (5) based on reparametrizing A = U U * [START_REF] Burer | A nonlinear programming algorithm for solving semidefinite programs via low-rank factorization[END_REF].

Discussion

Convex duality. Following [START_REF] Rudi | Finding global minima via kernel approximations[END_REF], we can provide a dual interpretation to the use of PSD models. Indeed, the minimimization problem we solve is

inf µ probability measure X f (x)dµ(x)
which we reformulate as

inf µ signed measure X f (x)dµ(x) such that X dµ(x) = 1 and X Φ(x)Φ(x) * dµ(x) ≽ 0.
Given that we expect the solution of the original problem to Dirac measures supported at global minimizers, we can add constraint that are satisfied by Diracs, such as,

X |dµ(x)| ⩽ 1 or Ω(µ) ⩽ 1,
for any norm Ω on signed measure that is larger than the total variation norm. The first constraint leads to a dual problem sup

c∈R, B≽0 c -f -c1 -ϕ(•) ⊤ Bϕ(•) ∞ ,
while the second one leads to

sup c∈R, B≽0 c -Ω * f -c1 -ϕ(•) ⊤ Bϕ(•) .
It turns out that the dual of the S norm and of F norm are domininating the total variation norm, and thus have a dual interpretation. Thus, our method for obtaining a posteriori certificates directly extends to optimization problems that are defined through probabilty measures and already tackled by kernel sum-of-squares, such as optimal transport [START_REF] Vacher | A dimension-free computational upper-bound for smooth optimal transport estimation[END_REF], or optimal control [START_REF] Berthier | Infinite-dimensional sumsof-squares for optimal control[END_REF].

Comparison to previous work on kernel sums-of-square. Compared to [START_REF] Rudi | Finding global minima via kernel approximations[END_REF], the subsampling is now done differently than before: the constraint X ϕ(x)ϕ(x) * dµ(x) ≽ 0 is replaced by the projection on the span of ϕ(x (1) , . . . , ϕ(x (n) ) being a positive semidefinite matrix. This is a relaxation in the dual, which still leads to a lower bound for the optimization problem.

Figure 1: For f : R → R as described in Appendix F, we use PSD models using feature maps of dimension n = 25 and n = 100. To the left, we see that for n = 25, the model g(x) + c does not approximate f well, so our a posteriori error guarantee is 0.24. But, when n = 100, the model approximates f much better, and our a posteriori guarantee is 0.01. To the right, we plot the absolute difference between fk and (g + c) k for small k's, which is what drives the difference in performance between n = 25 and n = 100.

This also suggests a candidate optimal solution when applied to the torus. Indeed, at optimality, we expect µ to be close to a Dirac at x 0 , and then (in 1D for simplicity), μ1 should be close to e -2iπx 0 , and we can read off a candidate minimizer as the argument of the first Fourier coefficients (we could imagine using more than one).

Figure 2: For f : R 2 → R as described in Appendix F, we plot the a posteriori error guarantee of our algorithm's estimate vs. n, i.e min j f (x j ) -ĉ which is the difference between the minimum value of f achieved on a random grid of points, which upper bounds c * , and our estimate of the minimum, ĉ1 , which lower bounds c * . The error is at most 10% of the function's range with n = 150, and can be made less than 1% with n = 750.

A Proof of Lemma 7

Proof. The proof of the adaptation of Corollary 2 of [START_REF] Rudi | Finding global minima via kernel approximations[END_REF] to the periodic setting is organized in four steps that are summarized in this paragraph. By translation f (x) = f (x -b) for a suitable vector b ∈ X, we show that all the global minima are contained in a closed set A strictly contained in an open set Ω strictly contained in the closed set [0, 1] d . Then Corollary 2 of [START_REF] Rudi | Finding global minima via kernel approximations[END_REF] shows that there exist q functions u 1 , . . . u q ∈ C m+2 (R d ) such that f (x) = q j=1 u j (x) 2 for any x ∈ Ω. Then we build two bump functions such that 0 ≤ η, ν ≤ 1 and η 2 + ν 2 = 1 on R d , moreover η = 0 on X \ Ω and 1 on A (so ν = 1 on X \ Ω and 0 on A).

Then we create a new function ṽ = f • ν and we show that its periodic extension ṽper satisfies

ṽper ∈ C m+2 (R d ).
We prove the same for ũj,per , the periodic extension of ũj = u j • ν. The result is obtained by noting that ũ1,per , . . . , ũq,per , ṽper are m-times differentiable periodic functions that satisfy f (x) = ṽper (x + b) 2 + j ũj,per (x + b) 2 for any x ∈ R d .

Step 1. Translating the functions, applying Corollary 2 of [START_REF] Rudi | Finding global minima via kernel approximations[END_REF]. Let s ∈ N and M = {x 1 , . . . , x s } be the set of minimizers on [0, 1) d . First, we work with a translated periodic

function f = f (x -b), where b = (τ /2, . . . , τ /2) ∈ R d and τ is the minimum distance of a point in M from [0, 1] d \ [0, 1) d (note that τ > 0 by construction). Let A = [2τ /3, 1 -2τ /3] d and Ω = ∪ x∈(τ /2,1-τ /2) d B τ /6 (x), where B τ /6 (x)
is the open ball of radius τ /6 centered in x. Note that A ⊂ Ω ⊂ X and that A, X are closed, while Ω is open. Now in the translated version f all the zeros are in A. By applying Corollary 2 of [START_REF] Rudi | Finding global minima via kernel approximations[END_REF] on f and Ω, we obtain that there exists q ∈ N and u 1 , . . . , u q ∈ C m (R d ) such that

f (x) = q j=1 u j (x) 2 , ∀x ∈ Ω.
Step 2. Building the bump functions. Let now α, β be two infinitely differentiable non-negative functions on R d such that α = 0 on R d \ Ω and is strictly positive on Ω, while β = 0 on A, strictly positive on

R d \ A. Since α 2 + β 2 > 0 on R d , the function √ • ∈ C ∞ ((0, ∞)) and α, β ∈ C ∞ (R d ), then η = α/ α 2 + β 2 and ν = β/ α 2 + β 2 are C ∞ (R d ) and satisfy η = 0 on R d \ Ω, η = 1 on A, analogously ν = 1 on R d \ Ω and 0 on A, moreover η 2 + ν 2 = 1 on R d .
Step 3. Construction of ṽper and ũj,per

. Let C = [τ, 1 -τ ] d . By construction, {x 1 , . . . , x s } ⊂ C and so f > 0 on the set X \ C. Then ( f ) 1/2 ∈ C m+2 (X \ C) since it is the composition of √ • ∈ C ∞ ((0, ∞)) and f that is f > 0 on X \ C. Then ṽ = ( f ) 1/2 • ν ∈ C m+2 (X), since ( f ) 1/2 is m + 2 times differentiable
on the set X \ C and ν is infinitely differentiable and 0 on A ⊃⊃ C. Denote by ṽper the periodic extension of ṽ, i.e. ṽper (x + k) = ṽ(x) for any x ∈ X and k ∈ Z d .

We now prove that ṽper satisfies ṽper ∈ C m+2 (R d ). First note that it is m + 2 times differentiable in the interior of each cube X + k for k ∈ Z d , since ṽ has this property on the interior of X. Moreover it has the same property also in a neighbourhood of the set S + k with S = [0, 1] d \ (0, 1] d and k ∈ Z d . Indeed, let B τ /6 (x + k) be the open ball of radius τ /6 around x + k, with x ∈ S and k ∈ Z d . On B τ /6 (x + k) the function ṽ is equal to ( f ) 1/2 , which in that region is m + 2 times differentiable, since we are in a translation of X \ Ω.

Define now ũj = u j • η and denote by ũj,per its periodic extension. Similary to the case of ṽper , since ũj is identically 0 on X \ Ω and it is m-times differentiable on X, we can prove that the periodic extension of ũj satisfies ũj,per ∈ C m (R d ).

Step 4. Conclusion. Note that R d = ∪ k∈Z d X + k . Then, by expanding the definitions, we have that for all x ∈ X, k ∈ Z d ,

ṽper (x + k) 2 + q j=1 ũj,per (x + k) 2 = ṽj (x) 2 + q j=1 ũj (x) 2 = (( f ) 1/2 (x)ν(x)) 2 + q j=1 (u j (x)η(x)) 2 = f (x)ν(x) 2 + η(x) 2 q j=1 u j (x) 2 =      f (x) x ∈ X \ Ω f (x)ν(x) 2 + η(x) 2 q j=1 u j (x) 2 x ∈ Ω \ A q j=1 u j (x) 2 x ∈ A = f (x),
where in the last two steps we use the fact that j u j (x) 2 = f (x) on x ∈ Ω and the fact that η 2 + ν 2 = 1 everywhere and, in particular, η = 0 on X \ Ω, and on A, while ν = 1 on X \ Ω and 0 on A. The proof is concluded by taking the translated version of by the vector b. I.e. z j (x) = ũj,per (x + b) for any x ∈ R d and j = 1, . . . , q and moreover z q+1 (x) = ṽper (x + b) for any x ∈ R d , and Q = q + 1.

B Proof of Theorem 8

Proof. First, let s = m + d/2. Note that, by applying Lemma 7 to f , we have that there exist Q ∈ N functions u 1 , . . . , u Q , that are periodic, s-times differentiable and which provide a new characterization of f as f = Q j=1 u 2 j . The desired result is obtained by applying Theorem 6 to this characterization of f . To apply Theorem 6, we need to find a suitable ρ ∈ ℓ 1 (Z d ) such that H ρ contains u 1 , . . . , u Q . In particular, we choose ρ k = (1 + d j=1 (2πk j ) s ) -2 . We prove now that u j ∈ H ρ and we characterize the resulting convergence rate.

Denote i =

√ -1 and by ∂ q j u the function ∂ q j u = ∂ q ∂x q j u, for all j ∈ {1, . . . , d} and q ∈ {0, . . . , s}, and u an s-times differentiable periodic function. Denote by ∥u∥ C s (X) the following norm ∥u∥ C s (X) = max 1≤j≤d max 0≤q≤s ∥∂ q j u∥ L ∞ (X) . With the notation we are using of the Fourier series we have ( ∂ q j u) k = (2πik j ) q u k (see, e.g., Prop. 3.1.2 of Grafakos, 2008). Now, by the Plancherel's identity (Prop 3.1.16 of [START_REF] Grafakos | Classical Fourier Analysis[END_REF] and the fact that X has volume equal to 1,

k∈Z d (2πk j ) 2q | u k | 2 = k∈Z d |(2iπk j ) q u k | 2 = X |∂ q j u| 2 dx ≤ ∥∂ q j u∥ 2 L ∞ (X) ≤ ∥u∥ 2 C s (X) < ∞,
where the norm ∥u∥ C s (X) is finite, since, for any u that is periodic and s-times differentiable, we have that ∂ q ∂x q j u, with q ≤ s, is also continuous and periodic, so uniformly bounded on X. Now, since ( d j=0 a j ) 2 ≤ c d d j=0 a 2 j for any a j ≥ 0, with c = d + 1, by expanding the definition of ρ k ,

k∈Z d | u k | 2 ρ k = k∈Z d u k + d j=1 (2πk j ) s u k 2 ≤ c k∈Z d | u k | 2 + d j=1 k∈Z d |(2iπk j ) s u| 2 ≤ c∥u∥ 2 C s (X) .
(8) This proves that the functions u 1 , . . . , u Q , that are periodic and s-times differentiable, belong to the space H ρ . Now we can apply Theorem 6, which gives min g∈Gt ∥f -g∥ F ≤ C ′ f R t , for any t ∈ N, where now

C ′2 f = c Q j=1 ∥u j ∥ F ∥u j ∥ C m (X)
and R t is bound as follows.

Since, ( d j=0 a s j ) 2 ≥ c s ( d j=0 a j ) 2s for any a j ≥ 0, with c s = (d + 1) -2(s-1) , (see, e.g., page 11 of [START_REF] Grafakos | Classical Fourier Analysis[END_REF], and the cardinality of the set of vectors in Z d summing up to a given number corresponds to # k |k| = r = r+d-1 d-1 ≤ C d r d-1 for any r, d ∈ N (e.g., page 52 of [START_REF] Brualdi | Introductory Combinatorics[END_REF], with C d = (2e) k-1 , we have

R 2 t = |k|>t ρ k ≤ |k|>t 1 1 + c s |k| 2s = r>t # k |k| = r 1 + c s r 2s ≤ r>t C d r d-1 1 + c s r 2s ≤ C d (s -d) c s t -(2s-d) ,
where, in the last step, we used the fact that

C d r d-1 1+csr 2s ≤ C d cs r -(2s-d+1) , moreover, r>t r -(2s-d+1) ≤ ∞ t x -(2s-d+1) dx = t -(2s-d) 2s-d . The final constant is then C 2 f = C ′2 f C d /((s -d)c s ).

C Proof of Theorem 10

Proof. When L k is concave and G-Lipschitz w.r.t. the Frobenius norm for all k, then the average of the iterates of projected stochastic gradient ascent with optimal constant stepsize applied to a problem of the form in Eq. ( 6) will have error bounded by (Nemirovski et al., 2009, Proposition 2.2)

c -E k∼π [L k ( ĀT )] ≤ a 0 • GR log(2/δ) √ T (9)
with probability at least 1 -δ. For our particular objective, it is easy to see that for all k ∈ Z d , sup

A,A ′ |L k (A) -L k (A ′ )| ∥A -A ′ ∥ F rob. ≤ ∥M (k) ∥ F rob. π k .
Therefore, with our choice of π k ∝ ∥M (k) ∥ F rob. + (1 + d j=1 (2πk j ) d+1 ) -1 , we can bound the parameter of Lipschitz continuity for all k by

G ≤ k∈Z d   ∥M (k) ∥ F rob. + 1 + d j=1 (2πk j ) d+1 -1   ≤ 1 + k∈Z d ∥M (k) ∥ F rob.
Plugging this into Eq. ( 9) completes the proof.

D Proof of Theorem 11

Proof. First, we prove the high-probability bound. We begin by arguing that L

k (A) is bounded. Let µ k = (1 + d j=1 (2πk j ) d+1 ) -1 so that π k ∝ ∥M (k) ∥ F rob. + µ k . For each k ̸ = 0, we have |L k (A)| = 1 π k fk -A, M (k) ≤ k ′ ∈Z d ∥M (k ′ ) ∥ F rob. + µ k ′ ∥M (k) ∥ F rob. + µ k | fk | + ∥A∥ F rob. ∥M (k) ∥ F rob. ≤ 1 + k ′ ∈Z d ∥M (k ′ ) ∥ F rob. | fk | µ k + ∥A∥ F rob.
. Now we need to bound | fk |/µ k . Note that, by applying Eq. ( 8), with u = f , s = d+1 and ρ k = µ 2 k , we have that for any k

∈ Z d , |f k | µ k ≤ k∈Z d |f k | 2 µ 2 k 1/2 ≤ √ d + 1∥f ∥ C d+1 (X) , where ∥f ∥ C d+1 (X) = max j=1,...,d max q=1,...,d+1 ∥ ∂ q ∂x q j f ∥ L ∞ (X)
. The result then follows by Hoeffding's inequality: for any δ ∈ (0, 1)

P E k∼π [L k (A)] - 1 K K i=1 L k i (A) ≥ ( √ d + 1∥f ∥ C d+1 (X) + ∥A∥ F rob. ) 2 log(2/δ) K 1 + k∈Z d ∥M (k) ∥ F rob. ≤ δ.
Rearranging and noting that E k∼π [L k (A)] ≤ c ≤ c * , completes the first half of the proof.

For the second set of bounds, we note that

c ≥ E k∼π [L k (A)] = k:|k|≤K π k L k (A) - k:|k|>K fk -A, M (k) ≥ k:|k|≤K π k L k (A) - k:|k|>K | fk | + ∥A∥ F rob. ∥M (k) ∥ F rob.
This completes the proof.

E Bound on ∥M (k) ∥ F rob. for ϕ t

The k 1 , k 2 entry of (ϕϕ * )(x) is equal to e k 1 (x)e k 2 (x) * , so

[M (k) ] k 1 ,k 2 = [0,1] d e k 1 (x)e k 2 (x) * e -2πix ⊤ k dx = [0,1] d e -2πix ⊤ (k+k 1 -k 2 ) dx = 1 {k=k 2 -k 1 } .
Therefore, the entries of M (k) are bounded by 1, and if |k| ≥ 2t, then M (k) = 0. Therefore, we can bound

∥M (k) ∥ F rob. ≤ n1 {|k|≤2t} . Since #{|k| ≤ 2t} ≤ (4t + 1) d ≤ 8 d t d when t ≥ 1, then k∈Z d ∥M (k) ∥ F rob. ≤ n8 2d t 2d .

F Experimental Details

Here, we describe how we applied our approach to two simple non-convex optimization problems to show its promise. As described, we can lower bound c * by solving the stochastic concave maximization problem Eq. ( 6) using an algorithm like projected stochastic gradient ascent. However, each iteration requires projecting the algorithm's iterate onto the PSD cone, which is computationally expensive.

Therefore, following [START_REF] Burer | A nonlinear programming algorithm for solving semidefinite programs via low-rank factorization[END_REF], we reparametrize A = U U * , which is always positive semidefinite, using new parameters U ∈ C n×n , yielding the unconstrained objective

c = max U ∈C n×n E k∼π [L k (U U * )]. (10)
Due to the non-linear reparametrization, the objective is no longer concave, but just as [START_REF] Burer | A nonlinear programming algorithm for solving semidefinite programs via low-rank factorization[END_REF] exhibit for linear SDPs, we find that stochastic gradient ascent on U succeeds for our problem when we optimize a smooth surrogate for our non-differentiable objective. Specifically, for k ̸ = 0, we replace

L k (A) = -1 π k fk -A, M (k) → Lk (A) = -1 π k (απ k ) 2 + fk -⟨A, M (k) ⟩ 2 ,
where α is small scalar. Choosing α larger makes the objective smoother, but makes Lk a worse approximation of L k . In our experiments, we tune α, along with the other hyperparametersincluding the stepsize, η and the number of iterations, T -with cross validation.

A random non-convex objective. We constructed a family of non-convex periodic functions on [0, 1] and [0, 1] 2 to test our algorithm. The functions are defined in terms of their Fourier series, with fk ∼ N(0, 1/(1 + |k|) 2 ) + i • N(0, 1/(1 + |k|) 2 ) for each k with |k| ≤ 15 in the 1D case and |k| ≤ 4 in the 2D case. We then adjust the Fourier components so that they satisfying the necessary property fk * = f * k . The value of f itself is then computed on a grid of points, x 1 , . . . , x N , and is rescaled by dividing by max i f (x i ) -min j f (x j ) so that f 's range is of order 1. Since ϕ i (x) = ϕ(x -x i ) and ϕ j (x) = ϕ(x -x j ), we need to know how to compute the k-th Fourier coefficient of x → φ(x -y)φ(x -z). We have: φ(x -y)φ(x -z) = n,m∈Z ρ |n|+|m| e 2iπn(x-y)+2iπm(x-z) = n,m∈Z ρ |n|+|m| e -2iπny-2iπmz e 2iπ(n+m)x .

For simplification and by symmetry, we can consider T (y, x) = φ(x -y)φ(x), so that ϕ(xy)ϕ(x -z) = T (y -z, x -z). Thus, this k-th Fourier coefficient is simply 1-ρ 2 e 2iπ(y-z) + e -2iπ(y+kz) -e -2iπ(ky+z) e -2iπz -e -2iπy + e -2iπky 1-ρ 2 e -2iπ(y-z) y ̸ = z ρ k e -2ikπz k + 1+ρ 2 1-ρ 2 y = z.

Therefore, we have that the i, jth entry of M (k) for k ≥ 0 is given by

[M (k) ] ij =    ρ k e -2ikπx j
1-ρ 2 e 2iπ(x i -x j ) + e -2iπ(x i +kx j ) -e -2iπ(kx i +x j ) e -2iπx j -e -2iπx i + e -2iπkx i 1-ρ 2 e -2iπ(x i -x j )

x i ̸ = x j ρ k e -2iπx i k + 1+ρ 2 1-ρ 2

x i = x j , and for k < 0, we have M (k) = M (-k) * . This allows us to compute [M (k) ] ij in the 1D case, which is the above function of k, x i , and x j ; denote this function h(k, x i , x j ). For the multidimensional case, since the feature map decomposes over coordinates, we simply have

[M (k) ] ij = d a=1 h(k, x i [a], x j [a]).
With this in hand, we can implement our algorithm.

Special cases and bounds. Now, we try to control ∥M (k) ∥ F rob. , which is needed to compute the a posteriori error guarantees.

First, we note that in the special case k = 0, we get:

[M (0) ] ij = 1 -ρ 4 1 + ρ 4 -2ρ 2 cos 2π(x i -x j )
.

Moreover, we have: T (0) k = ρ |k| |k| + 1 + ρ 2 1 -ρ 2 . We also have, since φ is always non-negative (see ( 11 Therefore, in 1D

|[M (k) ] ij | = |e -2iπkx j T (x i -x j ) k | ≤ | T (x i -x j ) k | ≤ T (0) k = ρ |k| |k| + 1 + ρ 2 1 -ρ 2
Therefore, we can upper bound in 1D

∥M (k) ∥ F rob. = n i,j=1 [M (k) ] 2 ij ≤ nρ |k| |k| + 1 + ρ 2 1 -ρ 2 .
In multiple dimensions, we can further bound

∥M (k) ∥ F rob. ≤ n d d i=1 ρ |k i | |k i | + 1 + ρ 2 1 -ρ 2 ≤ n d ρ |k| |k| d + 1 + ρ 2 1 -ρ 2 d = n 1 -ρ 2 1 + ρ 2 d ρ |k| 1 + |k| 1-ρ 2 1+ρ 2 d d ≤ n 1 -ρ 2 1 + ρ 2 d ρ |k| e |k| 1-ρ 2 1+ρ 2 = n 1 -ρ 2 1 + ρ 2 d ρe 1-ρ 2 1+ρ 2 |k| = ζ ρ|k| ,
where ζ is a constant independent of k and ρ < 1. Therefore, ∥M (k) ∥ F rob. decays exponentially quickly as |k| increases, which ensures that k∈Z d ∥M (k) ∥ F rob. is finite and not too large, and that k:|k|>K ∥M (k) ∥ F rob. goes to zero as K increases, which can allow for a tight a posteriori guarantee using Theorem 11.

T

  (y) k = n+m=k ρ |n|+|m| e -2iπny = n∈Z ρ |n|+|n-k| e -2iπny .Moreover, we will need to compute e -2ikπz T (y -z) k . We directly have T (y) -k = T (y) * k , so we can consider k ≥ 0 andT (y) k = 0 n=-∞ ρ k-2n e -2iπny +

  φ(x -y)dx = T (0) k .

We need complex numbers because of Fourier analysis, but this extends to symmetric real matrices.

Acknowledgements This work was supported by the French government under the management of the Agence Nationale de la Recherche as part of the "Investissements d'avenir" program, reference ANR-19-P3IA-0001 (PRAIRIE 3IA Institute). We also acknowledge support from the European Research Council (grants SEQUOIA 724063 and REAL 947908).

A different feature map The feature map that we use for our experiments has the form ϕ n,ρ (x) = φn,ρ (x[1]) • φn,ρ (x[2]) • • • • • φn,ρ (x[d]), where n ∈ N and ρ ∈ (0, 1) are hyperparameters to be chosen later, φn,ρ : R → C n , and • denotes the hadamard product, so the feature map decomposes over the coordinates of x. To define φn,ρ , we sample n points x 1 , . . . , x n uniformly at random from [0, 1] d , and set

The function φ is chosen so that its Fourier components φk = ρ |k| decay exponentially quickly with k. In Appendix F.1 below, we show how to compute the matrices M (k) that are needed to implement our algorithm, and we bound ∥M (k) ∥ F rob. , which is needed to compute the a posteriori guarantees. For this feature map, larger n allows for a more expressive, but more computationally expensive PSD model and Figures 1 and2 demonstrate the effect of n on our a posteriori accuracy guarantees in one and two dimensions, respectively. The parameter ρ, which we choose using cross-validation, clearly affects the Fourier components of the PSD model that we learn, with smaller ρ making them decay more quickly with k.

F.1 Analysis of the Feature Map

In what follows, we will drop the subscripts n and ρ and consider these hyperparameters to be fixed and arbitrary. We recall the definition

We further note that

In this Appendix, we show how to compute M (k) , the kth Fourier component of ϕϕ * , which is needed to implement our algorithm. Since ϕ

) decomposes across coordinates, this essentially boils down to computing the 1D version d times and multiplying across dimensions.

So, for now we focus on the case d = 1, and attempt to compute


Bruno Grouiez ¥

Virginie Zeninari
email: virginie.zeninari@univ-reims.fr

Lilian Joly

Bertrand Parvitte

Pulsed quantum cascade laser spectroscopy with intermediate-size pulses: Application to NH 3 in the 10 µm region

Keywords: 07, 07, Df Sensors; 07, 57, Ty Infrared spectrometers

come     Pulsed quantum-cascade-laser spectroscopy with intermediate-size pulses: application to NH3 in the 10 µm region B. Grouiez, Virginie Zeninari, Lilian Joly, Bertrand

Intermediate-size pulse quantum cascade laser spectrometry

Quantum-cascade lasers (QCLs) are ideal mid-IR sources for trace gas monitoring. They can be fabricated to operate at any of a very wide range of wavelengths from 3 to 30 µm where the absorption gas coefficients are stronger than in the near-infrared. Various techniques using continuous wave (cw) QCLs have been tested in the Reims laser group demonstrating very beautiful results from a spectroscopic point of view [START_REF] Joly | Spectroscopic study of the ν 1 band of SO 2 using a continuous wave DFB QCL at 9.1 µm[END_REF][START_REF] Parvitte | Preliminary results of heterodyne detection with quantum cascade lasers in the 9.1µm region[END_REF][START_REF] Joly | Water vapor isotope ratio measurements in air with a quantum cascade laser spectrometer[END_REF][START_REF] Grossel | New improvements in methane detection using a Helmholtz resonant photoacoustic laser sensor: a comparison between near-IR diode lasers and mid-IR quantum cascade lasers[END_REF][START_REF] Grossel | Photoacoustic detection of nitric oxide with a Helmholtz resonant quantum cascade laser sensor[END_REF][START_REF] Grossel | Optimization of a compact photoacoustic quantum cascade laser spectrometer for atmospheric flux measurements: application to the detection of methane and nitrous oxide[END_REF][START_REF] Joly | Development of a spectrometer using a cw DFB quantum cascade laser operating at room temperature for the simultaneous analysis of N 2 O and CH 4 in the Earth's atmosphere[END_REF][START_REF] Grossel | Quantum cascade laser spectroscopy of N 2 O in the 7.9 µm region for the in situ monitoring of the atmosphere[END_REF]. A comparison between the cw and the pulsed modes for a quantum cascade laser has been realized in [START_REF] Grouiez | Comparison of a Quantum Cascade Laser used in both cw and pulsed mode. Application to the study of SO 2 lines around 9 µm[END_REF]. Despite the good results obtained with cw devices, pulsed mode operation has long been the main way for using QCLs [START_REF] Faist | Quantum cascade lasers[END_REF] and remains the only available technique for several spectral ranges [START_REF] Devenson | High temperature operation of 3.3 µm quantum cascade lasers[END_REF].

Two different ways of pulsing the laser have been demonstrated.

The first way is used by most spectrometers based on pulsed QCLs [START_REF] Namjou | Sensitive absorption spectroscopy with a room temperature distributed feedback quantum cascade laser[END_REF][START_REF] Kosterev | Trace gas detection in ambient air with a thermoelectrically cooled QC-DFB laser[END_REF] and is based on a short-duration current pulse applied to the QCL. This current generates a pulse in the spectral domain, which can be tuned through the absorption feature of interest by use of a slowly varying subthreshold current ramp or by use of a variation of the temperature of the laser chip but at much lower speed, owing to the high thermal capacity of the laser submount and laser base. The main drawback associated with this method is related to the laser line broadening that reduces spectral resolution and sensibility. The "line broadening" that is often attributed to the laser linewidth is in fact related to the detector's bandwidth [START_REF] Grouiez | Alternative method for gas detection using pulsed quantum cascade laser spectrometers[END_REF].

The other method for pulsing QCL differs from the previous one. It uses the linear frequency downchirp that is induced when a long top-hat current pulse is applied to a QCL. The temporal profile of a generated pulse is recorded and used to infer the spectral nature of the pulse. The extent of the frequency window that is covered depends on the duration of the applied current pulse [START_REF] Normand | Characterisation of the spectral behaviour of pulsed quantum cascade lasers using a high resolution Fourier transform infrared spectrometer[END_REF]. The spectral profiles of molecular gases can be recorded in real time [START_REF] Grouiez | Comparison of a Quantum Cascade Laser used in both cw and pulsed mode. Application to the study of SO 2 lines around 9 µm[END_REF][START_REF] Normand | Fast, real-time spectrometer based on a pulsed quantum cascade laser[END_REF]. With this "long pulses" method, the resolution of the QCL spectrometer is limited not by the effective linewidth but by the temporal resolution of the detection system. Under some experimental conditions (low-pressure spectra, for example) the correlation between temporal and spectral profiles is not necessary straightforward. Particularly, if the duration of the interaction between the laser pulse and the molecular gas is short in comparison with the collisional lifetime, saturation and rapid-passage effects can occur as demonstrated in [START_REF] Mcculloch | Observation of saturation and rapid passage signals in the 10.25 micron spectrum of ethylene using a frequency chirped quantum cascade laser[END_REF][START_REF] Mcculloch | Highly sensitive detection of trace gases using the time-resolved frequency downchirp from pulsed quantum cascade lasers[END_REF].

Although this technique has numerous advantages when compared to the short-pulse technique, it suffers from several drawbacks: it requires much more specific instrumentation than the shortpulse technique, such as high-speed detectors and preamplifiers, a high-speed acquisition board, or a digital sampling oscilloscope. Moreover, the use of such long pulses is not recommended for most of the currently commercially available pulsed QCLs provided by Alpes Lasers [19]. Lasers allowing the use of long pulses seem to be available only for fewer wavelengths.

As none of those methods was completely satisfactory, a third way has been developed that uses intermediate length pulses, typically between 30 and 80 ns [START_REF] Grouiez | Alternative method for gas detection using pulsed quantum cascade laser spectrometers[END_REF]. This interval of pulse duration was chosen to reach the linear frequency variation area that follows the transient at the beginning of the pulse without damaging the QC laser. We expected to obtain a narrower linewidth than with short pulses. The spectral range that is scanned inside so short pulses is usually small and we use either a subthreshold ramp (duty cycle of 0.3% with a frequency of 100 Hz) or a substrate temperature variation (duty cycle of 1% and ramping rate of 0.06cm - 1 /min) to increase the spectral range. These experiments were performed with a setup described in Fig. 1.

Figure 1

The pulse generator (Berkeley Nucleonics Corporation (BNC 565)) delivers 70ns impulsions with a 1% duty cycle that are used in conjunction with a temperature or a subthreshold current ramp in order to obtain the laser wavelength variation. When using a temperature ramp, the laser temperature is scanned step by step, with typically 100 ms between successive steps.

The series of impulsions are sent (using the driver from Alpes Laser starter kit which in part composed by a DEI PCO-7120 card) to the QC laser placed inside a Laser Laboratory Housing from Alpes Lasers. The laser beam passes through the absorption path cell containing ammonia (L = 20 cm) and is collected by SAT1903N high-speed MCT detector. The detector's signal is then sent to a digital sampling oscilloscope (Lecroy Waverunner104 Xi) so that each optical pulse is recorded. A part of the laser beam is sent through a Germanium etalon (with a FSR of 0.048 cm -1 ) to a Kolmar detector in order to obtain the wavenumber calibration. Absorption features can be seen either inside each pulse (i.e. versus time) or at a constant delay inside successive pulses, versus subtreshold ramp/substrate temperature.

Each of spectra presented in this paper have been built using the intermediate size method described in our previous paper [START_REF] Grouiez | Alternative method for gas detection using pulsed quantum cascade laser spectrometers[END_REF]. Each pulse is recorded by the digital oscilloscope and a post-treatment permits one to rebuild the spectrum by selecting a constant delay inside pulses.

The baseline due to the power variation of the laser during the ramp is also taken into account during this post-treatment operation. Fig 2 presents the experimental data corresponding to 3000 successive measurements recorded during a substrate temperature ramp ranging from 17 to 27°C. Each measurement is in fact the average of 10 successive pulses and corresponds to one step of the temperature scanning. The total time for one measurement of about 7 µs is negligible when compared to the duration of one temperature step. Figure 2a) is a 3D representation of the 3000 measurements made between 17 and 27°C. The signal amplitude (measured optical power) is coded using the color scale given on the right of the figure, from blue (no signal) to red (maximum of signal). The pulse recorded at T=21°C is plotted in black.

Figure 2b corresponds to a top view of Fig2a

). Each measurement corresponds to one horizontal line on Fig. 2b). The horizontal scale represents the pulse length and the vertical scale gives the substrate temperature corresponding to each measurements. According to the color coding, the maximum of absorption appear as the thin blue lines.

Figure 2 Figure 2b) can be used to follow the moving of the absorption line inside the pulse during the substrate temperature ramp. This is materialized by the arrows on the Figure 2b). This moving is an important clue to understand the variation of the instantaneous wavenumber versus the substrate temperature. In the middle of the pulse (typically between 1 and 1.4 10 -7 s on Fig. 2b), the emitted wavenumber varies in a very regular manner: the frequency downchirp can be exactly compensated by a substrate temperature variation and the absorption line shifts nearly linearly. In both transients at the beginning and at the end of the pulse, the emitted wavenumber varies very differently at a given substrate temperature the same absorption feature can be seen twice in one pulse (see for example between 22 and 23°C and for a time delay about 0.8 10 -7 s). The arrows allow us to follow the value of the instantaneous chirp frequency which is directly linked to the variation in heating during the pulse. Once can see that the sense of this one change rapidly and is also more important in transients than in the rest of the pulse.

Figure 2 also allows us to compare the spectral range due to frequency downchirp and to the substrate temperature variation: at a given temperature (21°C for example) only one absorption feature (at 967.997761 cm -1 ) can be seen inside a pulse. The substrate temperature variation (for example for a delay inside the pulse of 1.2 10 -7 s) allows to scan several absorption features (one line at 967.997761 cm -1 and a doublet at 967.774745 and 967.738421 cm -1 ). The first alternative can be used for narrow and isolated spectral lines (i.e. at low pressure) whereas the second alternative should be preferred at higher pressures. Results on NH 3 in the 10 µm region will be presented in this paper: in the low pressure range ( pressures the spectrum due to the chirp of the laser into the pulse let oscillations appears after the absorption. This phenomenon is named rapid-passage effect and was observed with the intrapulse technique in [START_REF] Mcculloch | Observation of saturation and rapid passage signals in the 10.25 micron spectrum of ethylene using a frequency chirped quantum cascade laser[END_REF].

Figure 3

The oscillations observed during the pulse cannot be explained by the Beer-Lambert law. The lifetimes of the levels that are linked to the Einstein coefficients and that are the basics of this law cannot take into account the individual comportment of the atoms. In this case there is no more equilibrium in the interaction between radiation and matter. A more realistic representation can be obtained thanks to a semi-classical model as it was initially proposed by Duxbury team [START_REF] Mcculloch | Observation of saturation and rapid passage signals in the 10.25 micron spectrum of ethylene using a frequency chirped quantum cascade laser[END_REF][START_REF] Mcculloch | Highly sensitive detection of trace gases using the time-resolved frequency downchirp from pulsed quantum cascade lasers[END_REF][START_REF] Duxbury | Rapid passage induced population transfer and coherences in the 8 micron spectrum of nitrous oxide[END_REF]. In this model the electromagnetic field is represented in a classical manner via the Maxwell equations and the atomic part is described by the quantum mechanic. This description is known as the Maxwell-Bloch equations. The differential Bloch equations are obtained by combining the coherent effect of the wave with the effect of random phenomenon such as spontaneous emission and collisions. These equations can be presented by the following system [START_REF] Allen | Eberly Optical resonance and two-level atoms[END_REF]:

v δω u τ 1 dt u d 2 + - = (1) m v τ 1 u δω dt v d 1 2 Ω - - - = (2) 1 0 1 τ m m v Ω dt m d - + = (3)
In this system u and v represent the mean values of the two crossed variables of the Bloch system. They are linked to the electric dipolar moment because the stationary solutions of this system permit to determine the dielectric susceptibility of the medium. m corresponds to the inversion of population and m 0 represents the initial value of the inversion of population. δω represents the difference between the pulsation of the incident wave and the central pulsation of the transition and is related to the frequency downchirp tuning rate. τ 1 and τ 2 are linked to the lifetime of the levels. Finally Ω 1 represents the Rabi pulsation.

Analytical solutions of this system can be found when δω remains constant during the experiment. Due to frequency downchirp δω varies during each pulse and only numerical resolution of this system can be done. One important difference between the analytical solutions and the numerical solutions is that the oscillations in u , v and m do not have a constant frequency. Moreover, as it was demonstrated in [START_REF] Duxbury | Rapid passage induced population transfer and coherences in the 8 micron spectrum of nitrous oxide[END_REF], the resolution of this system do not allow to reproduce the observed signals and several effects must be included in the calculation:

• The velocity of the molecules affects the value of δω due to Doppler effect.

The system of coupled equation must be solved for several class of velocities (200). • The Rabi pulsation used in eq. 4a presents some recorded signals at several pressures and Fig. 4b presents the results of the calculations. At this moment we do not obtain full agreement between measurements and calculations. As only numerical solutions can be obtained, they cannot be used in a typical inversion process like least-squares fitting. However the overall agreement on the three pressure range is good.

The best solution to obtain spectra that are satisfactory from a spectroscopic point of view is to increase the pressure. In this case the lifetime τ 1 becomes very small and the second term of equation ( 3)

1 0 τ m m -
becomes larger than the first one 1 v Ω thus the oscillations disappear. Part 3 of this paper will present results obtained at higher pressures where rapidpassage effects do not exist anymore.

Mid-pressure (50 mbar range) spectra of ammonia

In this part either temperature or current ramps were added to the 70ns length duration pulses.

A spectrum corresponding to a 40 ns delay after the beginning of the pulse is presented in Fig. 5. This spectrum is very satisfactory from a spectroscopic point of view and can be directly compared to a calculation using Voigt profile and HITRAN parameters [START_REF] Rothman | The HITRAN 2004 molecular spectroscopic database[END_REF]. The fitted and experimental spectra are very similar. Numerous spectra were recorded at different pressures from 50 to 80 mbar. For all those spectra the collisional linewidth was fitted and compared to the HITRAN value. The fitting procedure has been described in detail in [START_REF] Zéninari | Laboratory spectroscopic calibration of infrared tunable laser spectrometers for the in situ sensing of the Earth and Martian atmospheres[END_REF]. The experimental and HITRAN collisional linewidths are plotted vs. pressure in Fig. 6. One can remark than the two lines are parallel however only the HITRAN curve passes through zero. This is the demonstration of the existence of an apparatus function. Note that the difference between the two techniques is that the temperature spectrum is longer to record (typ. 140 s) than the current one (typ. 10 ms). On the other side the current spectrum is less wide than the temperature one. For these pressures there is no problem but it may be a drawback for atmospheric pressure spectra.

Figure 7 Even if the experimental and HITRAN calculated spectra are similar, some differences remain. To confirm that these differences were due to an apparatus function, we calculate the HITRAN spectrum using Voigt profile and convoluted it by a Gaussian profile with Δν = 0.003 cm -1 . These calculations are also displayed in Fig 7a and7b. For these two figures one can clearly see that the experimental spectra and the new calculation are much well fitted. The small difference remaining in the two spectra (on the right part of the wing) is probably due to a remaining part of the rapid passage effect. The recorded and calculated spectra allows estimating the effective linewidth of the laser to be about 0.003 cm -1 , which is significantly shorter than typical values reported with shorter pulses (typically 0.015 cm -1 [START_REF] Namjou | Sensitive absorption spectroscopy with a room temperature distributed feedback quantum cascade laser[END_REF]). Thus the intermediate-size pulse technique is more satisfactory from a spectroscopic point of view.

Conclusion

Pulsed quantum-cascade-laser spectrometers are still used to detect atmospheric gases. We have previously demonstrated [START_REF] Grouiez | Alternative method for gas detection using pulsed quantum cascade laser spectrometers[END_REF] that neither the interpulse nor the intrapulse technique was totally satisfactory to use the pulsed QCLs. The interpulse (short-pulse) technique gives distorted spectra and the intrapulse (long-pulse) technique cannot be applied to all QCLs.

Thus we have developed the intermediate-size pulse technique for gas detection using pulsed QCL spectrometers. This technique was applied in this paper to record infrared spectra of ammonia in the 10 µm region. However NH 3 spectra recorded in the mbar range show typical oscillations after the absorption due to the rapid-passage effect. This phenomenon had been previously observed with the intrapulse technique [START_REF] Mcculloch | Observation of saturation and rapid passage signals in the 10.25 micron spectrum of ethylene using a frequency chirped quantum cascade laser[END_REF]. In this paper, the resolution of the Maxwell-Bloch equations has demonstrated a good agreement in comparison with experimental results at various pressures.

To realize correct spectroscopic measurements NH 3 spectra were recorded in the 50 mbar range where the oscillations disappear. Thus the intermediate-size pulse technique gives reliable measurements for NH 3 detection using a temperature or a current ramp. Moreover the typical apparatus function (0.003 cm -1 HWHM) is far lower from the typical apparatus function of the interpulse QCL spectrometers (0.015 cm -1 HWHM).

Captions

  1 mbar range -Part 2) and in the mid-pressure range (50 mbar range -Part 3). 2 Low pressure (mbar range) spectra of ammonia An example of recorded pulse of ammonia at low-pressure with the intermediate-size pulse technique and without added ramp is presented in Fig. 3. Fig. 3a features a 70 ns pulse at a laser temperature without absorption and Fig.3b features a recorded pulse at a laser temperature where ammonia absorbs. In both cases the first part of the pulse presents a transient and is followed by a regular frequency downchirp area. The last part of the pulse also presents a transient. The Fabry-Perot etalon transmission shown on Fig 3b allows us to determine the frequency downchirp tuning rate in this useful part of the pulse. At these low

  2 and 3 is related to the transition dipole moment. The arbitrary orientation of the molecule versus the laser field must also be included in the calculation. These two effects smooth the very fast oscillations observed in the solutions of Eqs 1, 2 and 3 and allow us to obtain a satisfactory agreement with observed signals. Some examples are presented on Fig 4.

Figure

Figure 5 Figure 6

 56 Figure 5Figure 6

Figure 1 :

 1 Figure 1: Schematics of the experimental set-up.

Figure 2 :

 2 Figure 2: (a) 3D view of successive experimental laser pulses (NH 3 /air spectrum, P=60 mbar, Concentration # 0.3%, Temperature # 296 K, Length = 20 cm). Successive laser pulses are presented versus the laser substrate temperature. The black trace corresponds to one pulse. (b) Top view of Fig2(a).

Figure 3 :

 3 Figure 3: (a) record of a pulse with the intermediate-size pulse technique and (b) observation of rapid-passage effect in a 70 ns pulse with 0.37 mbar of pure ammonia.

Figure 4 :

 4 Figure 4: (a) Recorded signals for different gas pressures. (b) Calculated signals for the same experimental conditions.

Figure 5 :

 5 Figure 5: Example of spectrum recorded with a temperature ramp of 10 K and comparison with a calculated spectrum using Voigt profile for a pressure = 71 mbar and concentration of NH 3 in air is 1.3 %.

Figure 6 :

 6 Figure 6: Fitted line width vs. pressure for the experimental spectra. The HITRAN values are also displayed. See text for more details.

Figure 7 :

 7 Figure 7: Transmission experimental spectra with (a) a temperature ramp (P = 58 mbar and 0.28 % of NH 3 in air) and (b) a current ramp (P = 51 mbar and 0.8 % of NH 3 in air). On both parts of the figure are also displayed the calculated spectra using Voigt profile with the HITRAN data and another calculation introducing a convolution of the HITRAN spectra with a Gaussian lineshape (see text for more details).

  Figure 1

Figure 6

 6 Figure 6

Acknowledgments

The authors are indebted to the Région Champagne-Ardenne for financially supporting of the PhD of Bruno GROUIEZ.


L Joly

V Zéninari

T Decarpenterie

J Cousin

B Grouiez

D Mammez

G Durry

M Carras

X Marcadet

B Parvitte

V Zeninari

Continuous-wave quantum cascade lasers absorption spectrometers for trace gas detection in the atmosphere

Keywords: Quantum cascade laser, Infrared gas detection; Atmospheric applications PACS Numbers: 07.57.Ty, 42.55.Px

Mid infra-red absorption spectrometry based on continuous-wave distributed feedback (DFB) quantum cascade laser (QCL) is more and more widely used for trace gas detection and pollution monitoring. The main advantages of this technique are high sensitivity, high selectivity and a potential for extreme compactness. Various examples of trace gas detection for atmospheric detection will be presented in this paper. Commercial QCLs available on the shelves were first implemented. A cryogenic QCL emitting at 6.7 µm was used to demonstrate the detection of water vapor and its isotopes. A room-temperature QCL was then used to simultaneously detect methane and nitrous oxide at 7.9 µm. Recently, we have developed a room-temperature top grating DFB QCL designed around 4.5 µm for the demonstration of N 2 O detection in the ppb range. Atmospheric applications of these spectrometers will be presented. The improvements of QCL performances make it now possible to develop instruments that are more and more compact and therefore compatible with in situ applications.

Absorption spectrometers for trace gas detection

Air pollution is currently a major environmental issue. Numerous chemical agents identified in urban and peri-urban atmospheres have a direct impact on the environment and/or on the human health as greenhouses gases (carbon dioxide, methane, nitrous oxide, ozone, CFCs …) and Toxic Industrial Chemicals (TIC) (i.e. : Ammonia, Sulfur dioxide, Sulfuric acid…). For either the detection or the control of chemical agents, a flexible sensor adaptable to the variety of pollutants has to be developed. The main characteristics of this sensor must be:

Very high sensitivity (in the ppb range); Capable to give fast measurements of pollutant densities (few seconds of acquisition); Portable device to make on site measurements (outdoor as well as indoor).

The currently available gas detection methods may be split into two groups, (i) direct methods, which monitor a physical parameter of the target gas, and (ii) indirect methods, which use a chemical reaction or indicator to assess the concentration of the gas being sensed. This division of methods may be further split into optical sensing techniques and non-optical techniques. Indirect methods are used in a variety of fashions depending on the specific need or application. They are usually based on low cost technologies while being able to reach sensitivity in the ppb range. The small size, low power consumption and relative durability of their designs lend themselves well to portable detection instruments. They are also used in a variety of fixed applications as well. However, they can be poisoned by exposure to adverse chemicals or conditions and can be sensitive to groups of gases rather than a specific gas in a group. They also may exhibit non-reversible behavior, which, in many cases, may be undesirable. In contrast, direct methods are generally very accurate and highly selective means of detecting gas and measuring their concentrations. For the most sensitive chemical analyses of gases, laboratory techniques such as gas chromatography (GC) or the combination of gas chromatography and mass spectrometry (GC/MS) are often applied. These techniques achieve both the highest sensitivity and high selectivity, but are complicated. Moreover, major disadvantages like long measurement time, the need for additional carrier gases, high instrument cost, high maintenance cost, limited detector lifetime and the need for frequently recalibration prevent these methods from being applied efficiently on a large scale and in various situations.

2

Direct spectroscopic optical methods are easier to implement and provide fast, robust, high sensitivity and selectivity since every molecule has a unique IR spectrum (fingerprint). Gas sensing based on Tunable Diode Laser Spectroscopy (TDLS) is an established technology that has progressed in parallel with the advances in near-infrared (NIR) semiconductor diode lasers in the 0.7-3 μm range [1][2][3]. TDLS has been the preferred technology for the purpose of gas detection in industrial applications (chemical emission monitoring, process control, high sensitivity trace detection). The robustness of TDLS makes it possible to manufacture sensors with long life time (20 years) that need no recalibration. However, NIR TDLS still shows some drawbacks which prevent to use it for in situ atmosphere detection : the use of NIR overtones much weaker than fundamental Mid Infra-Red (MIR) absorption lines limit their sensitivity and the size of the current high sensitivity systems does not make them transportable.

The use of QCL as a MIR source in TDL spectrometers permits to overcome these problems [4][5][6][7][8]. In recent years, the performance of QCLs has improved markedly and they are today the only semiconductor laser sources operating at room temperature in pulsed or continuous modes over the 3-12 µm wavelength range. This semiconductor technology has moreover a potential for extreme compactness, high lifetime, high reliability and low cost. QCLs were invented 15 years ago and have now reached a certain degree of maturity [START_REF] Faist | Quantum Cascade Laser[END_REF]. However, the business of QCLs is still in its infancy. There are today only two commercial companies in Europe which are really able to fabricate and deliver quantum cascade lasers (a Swiss SME called Alpes Laser founded in 1998 by one of the co-inventors of the device and Alcatel-Thales III-V Lab).

Single frequency operation of the laser is required for spectroscopy application. It is usually achieved by introducing a distributed feedback (DFB) structure into the QCL active region in order to favor a particular mode. Pulsed DFB QCLs are commercially available over the whole 4-10 µm range. Pulsed operation mode operation of those devices at Room Temperature (RT) on Thermal Electro Cooler (TEC) is now easily achieved. On the other hand, Continuous Wave (CW) RT operation lasers on TEC are more lab devices than product and their price is prohibitive. As the demand is going to rapidly increase, in particular due to the development of MIR gas sensing systems, the price should reduce in the near future.

The challenge of designing and fabricating single mode CW, RT QCL has been partially overcome using buried gratings. However, this technique is quite costly and heavy since it requires an overgrowth step. The use of surface metal grating to fabricate QCL DFB is a way to reduce cost [START_REF] Carras | Laser semi-conducteur à contre-réaction répartie forte[END_REF]. Indeed, one major advantage of surface metallic grating is that no overgrowth is necessary. Moreover, by this approach we can achieve the cm -1 precision criterion over the device emission wavelength range. This has allowed pulsed DFB QCL with both a high output power and a high SMSR (Side Mode Suppression Ratio) to be demonstrated. We showed that robustness of top grating design is such that the fabrication yield can be very high (close to 100%) [11]. By reducing the difficulty of fabricating such diodes, costs would be automatically scaled down. The combination of this approach with improved heat dissipation processing has been used in this work to proceed toward continuous waves operation at 4.5 µm.

In this paper we will present examples of realizations of direct absorption spectrometers based on quantum cascade lasers and applied to gas detection in the atmosphere: water vapor isotopes at 6.7 µm, methane and nitrous oxide at 7.9 µm and nitrous oxide at 4.5 µm. The schematic representation of such a spectrometer is presented in Fig. 1.

Figure 1: Schematic representation of a quantum cascade laser spectrometer

The QCL is driven by a temperature controller and a current source. The laser emission wavelength is scanned over the molecular transition by applying a triangular current ramp (typ. 100 Hz). The QCL beam is usually collected by a lens adapted to the laser wavelength and is separated into two parts by a beam splitter. The reflected beam is coupled with a Fabry-Pérot etalon used for relative frequency calibration. The second beam passes through a gas cell. The beams are focused onto two detectors adapted to the laser wavelength. Signals are further recorded and treated in order to obtain the atmospheric molecule concentrations.

All measurements are usually made at room temperature. To retrieve the gas concentration for a particular line we applied a nonlinear least-squares fit to the molecular transmission, using a Voigt profile to model the line shape. The fitting procedure was described for example in [12].

H O isotopes detection with a 6.7 µm cryogenic QCL

The first spectrometer presented in this paper was developed in order to measure H 2 O and its isotopologues. In the stratosphere, a constant increase in the amount of water vapor of roughly 1% per year is observed. Even if the stratosphere appears to be strongly dehydrated, a change in the weak content of H 2 O could have a strong effect on the radiative and chemical equilibrium of the stratosphere and further alter ozone recovery. The coupling processes that transport and mix humid tropospheric air masses into the stratosphere, essentially in tropical regions where strong convective events occur, are largely unknown [13][14][15]. The in situ measurement of H 2 O isotopologues can be helpful in addressing this issue, as water-vapor fractionation can be used for diagnosing the transport of an air mass through the upper troposphere (UT) and the lower stratosphere (LS) and associated dehydration mechanisms.

The developed spectrometer is represented in Fig. 2 [16]. The dimensions of this spectrometer are quite big: 160 cm x 60 cm x 60 cm. One reason is that the QCL from Alpes Lasers -Switzerland was housed inside a cryostat filled with liquid nitrogen every 8 h. The QCL works from 80 to 125 K, and the threshold current is 380 mA. The tunability range of the QCL is 1483.5 cm -1 at T=124 K and I=710 mA to 1489 cm -1 at T=87.4 K and I=400 mA. The average output power is #30 mW, the laser linewidth is approximately 3-4 MHz, and there are no mode hops over the tunability range. One scan of the laser emission wavelength with the temperature over selected water-vapor transitions goes from 1483 to 1487 cm -1 . [17]. The second reason of the big size of the spectrometer is that the cell is a multipass White-type cell of 48 m effective length with BaF 2 windows. The total air pressure was measured with an uncertainty of 0.5% by an MKS baratron manometer with a 1000 Torr full scale. The air sample pressures vary from 15 to 75 Torr. We studied eight samples. One example of a recorded spectrum with the fitted Voigt profile and the residual below is presented in Fig. 3.

The tunability of the QCL enables it to reach many water-vapor isotopologues transitions. Some of these lines are used for the retrieval of the isotopologues concentrations. The concentrations derived from the fitting procedure are compared with the calculations based on the saturated vapor pressure in air at the temperature measurement and the partial water-vapor pressure obtained from a hygrometer. The concentrations obtained from the fitting procedure and those of the hygrometer are in quite good agreement for H 2 18 O and H 2 16 O. The HDO value is quite different, as explained below. The experimental H 2 18 O / H 2 16 O isotopic ratio is 2014 parts in 10 6 whereas the standard value is 2006.

The experimental HDO / H 2 16 O isotopic ratio is 236 parts in 10 6 whereas the standard value is 156. Our value seems not to be reasonable. However, the White cell used for the static measurements was used previously for intensity measurements of HDO in the laboratory [18].

The sample gases for the study of HDO lines consist of evaporation from various mixtures of liquid-enriched H 2 O and liquid H 2 O. So our cell was probably polluted by the outgassing of HDO trapped in the optical cell used during the previous laboratory experiment.

Fig. 3: Example of a recorded spectrum for air and fitted spectrum from 1484.9 to 1485.2 cm -1 . The (experimental minus calculated) residual is shown in the lower panel. The current excursion can cover two water vapor isotopologues, H 2 18 O and H 2 16 O within 1 s.

For measurements performed directly in the field, the problem of contamination of outgassed or trapped water vapor will be eradicated by use of an optical multipass cell operated open to the atmosphere. Nevertheless, this experiment has demonstrated that quantum cascade lasers are ideal sources for these challenging developments. These results have been obtained with a liquid-nitrogen cooled source however the same results may be obtained nowadays with a room-temperature cw source as it will be demonstrated in the next paragraphs. [START_REF] Mosier | [END_REF].

A photograph of the developed spectrometer is depicted in Fig.

4 [21]. The RT QCL purchased from Alpes Lasers (Switzerland) emits in the 7.9 µm region on the ν 4 band of methane [17]. The dimensions of this spectrometer are always quite big: 120 cm x 80 cm x 40 cm. Since the availability of RT CW QCL, the size of the whole system is not anymore limited by the source. The current limitation comes mainly from the absorption multipass cell which is a new type of multiple-reflection cell (L = 140 m) presented in [22]. The QCL is mounted onto a copper substrate attached to a thermoelectric cooler and contained in a and CH 4 lines are isolated enough, and there is no interference with water vapor lines.

Pressure (MKS Baratron) and temperature (PT100) sensors permit measurements of these two important parameters for high resolution spectroscopy. The sample gas (outside the laboratory) is brought from atmospheric pressure to #40 mbar across the inlet orifice to reduce the pressure broadening of the absorption lines, increasing the distinctive nature of the absorbing species spectral signature. The system operates at a gas flow of 5 L/min. A dust filter is set at the inlet of the measurement head to protect the gas system and the mirrors of the cell from pollution. A WXT510 weather transmitter from Vaisala is used to measure thermodynamic conditions outside the laboratory. An example of a recorded spectrum is shown in Fig. 5. The line shape was modeled with a Voigt profile and the resulting residual is also displayed in Fig. 5. The residual features an M shape, which is typical when using the Voigt model approach. The fitting can be improved and the residual reduced by using profiles such as the Rautian-Sobel'man line shape functions. However the residual with the Voigt profile is < 1%. Moreover the use of more complicated profiles would increase the processing time with no significant improvement in the accuracy of the concentration retrieval [2]. The instrument has been used to record a series of ambient air measurements. Ambient air samples were collected from outside the laboratory building. The spectrometer was located in the laboratory. A filter was put at the tube entrance to prevent dust entering the optical cell.

The vegetation all around the laboratory was mainly a meadow with some oaks. The third spectrometer presented in this paper was developed in order to measure N 2 O in the atmosphere. It was particularly dedicated to N 2 O measurements for agronomical studies because agricultural soils contribute substantially to this budget. N 2 O emission occurs due to the use of mineral and organic fertilizers in agricultural activity, through the activity of microorganisms, and to the expirations of human beings [24]. The production of N 2 O by soils primarily occurs through two biological processes: (1) denitrification, the anaerobic microbial respiration using nitrate and nitrite as electrons acceptors [25] and (2) nitrification, the exothermic oxidation of ammonium to nitrite and nitrate [START_REF] Prosser | Nitrification[END_REF]. A better understanding of the soil-atmosphere exchange is necessary to control N 2 O emissions with new cropping systems.

N 2 O emissions are very variable with time and space [START_REF] Duxbury | [END_REF][28] and difficult to measure. Two main methods are available today [29]: (1) the static method that is based on the measurement of gas concentration within closed chambers that are put at the soil surface during a given period (typically 3 h) and ( 2) the micrometeorological method that is based on the measurement of the gradient of gas concentration in the atmosphere. Both methods require gas concentration measurements that are done either in the laboratory after atmospheric sampling or directly in the field. Available gas analyzers limit the possibility to measure the gas concentration on a wide scale and over a long period while modern atmospheric research on gas exchange between the biosphere and the atmosphere requires sensitive, reliable, and fast response chemical sensors. Therefore, techniques for fast and simultaneously sensitive trace gas measurements based on tunable diode laser absorption spectroscopy have been successfully applied to micrometeorological trace gas flux measurement techniques [30].

A photograph of the developed spectrometer is depicted in Fig. 6. This spectrometer is a faster and a more compact version of the previous described ones. The QCL especially developed in this work by Alcatel-Thales III-V Lab emits in the 4.5 µm region where a fundamental band of nitrous oxide is present. This wavelength region where nitrous oxide absorptions are stronger permits to reach the same kind of detection limit that the previous spectrometer (part 3 of this paper) but using a cell with a smaller path length. Thus the cell volume is also lower and may be refreshed more easily permitting fast measurements such as working at 6.7 µm was used to demonstrate the detection of atmospheric water vapor and its isotopes. A room-temperature 7.9 µm-QCL spectrometer was then used to simultaneously detect methane and nitrous oxide in air. These two spectrometers demonstrate that cw QCLs are ideal sources for trace gas detection in the atmosphere. Moreover QCLs are still in progress in terms of cw, RT, etc. Thus the spectrometer characteristics are also in progress and these spectrometers are more and more efficient and give reliable interesting measurements for various applications.

However the size of these spectrometers (typ. 120 cm x 60 cm x 40 m) mainly due to the multipass cell and to the laser source was still a brake to real in situ measurements. The third part of this paper presents a room-temperature QCL spectrometer which was used to demonstrate the detection of N 2 O around 4.5 µm. With this kind of spectrometer which is more compact by using low volume cell and RT laser, in situ atmospheric applications can now be realized. The preliminary results demonstrated in this paper show the capabilities of the instrument for agronomical applications. The next step will be to use this spectrometer in fields and to automate it.

Figure 2 :

 2 Figure 2: Three-dimensional view of the H 2 O isotopologues set-up

  hermetically sealed enclosure (Alpes Lasers LLH). Heat produced by the cooler is removed by circulating water through the base of the enclosure. The QCL operates between -30 °C and +30 °C and the maximum current is 550mA. The temperature is controlled and stabilized by a model ILX LightWave temperature controller and the current source is an ILX Lightwave model LDX-3232. The tunability range of this laser is from 1265 cm -1 (T QCL = -29 °C, I QCL = 550mA) to 1274 cm -1 (T QCL = -29 °C, I QCL = 290mA). The QCL power is #30mW. The synthetic spectrum calculated in this region shows the best spectral window to simultaneously detect N 2 O and CH 4 . This spectral range is from 1270.7 to 1271.2 cm -1 . In this region N 2 O

Figure 4 :

 4 Figure 4: Photograph of the CH 4 and N 2 O sensor set-up.

  Stability and laboratory tests of the instrument based on the Allan variance method have been realized[21]. The log-log plot of the Allan variance versus the integration time (τ) in seconds demonstrates than the Allan variance decreases and, hence, only the white noise is dominant up to 1280 s integration time. The variance minimum at τ Allan = 1280 s corresponds to a minimum σ Allan of 0.5 ppb; and for CH 4 the variance minimum at τ Allan = 1280 s corresponds to a minimum σ Allan of 3 ppb. These are the maximum precisions reachable, which satisfy the field measurement requirements asking for sensitivities to flux variability of 1 ppb for N 2 O and 10 ppb for CH 4 .

Figure 5 :

 5 Figure 5: Example of a recorded spectrum for air and fitted spectrum from 1270.65 to 1271.15 cm -1 . The (experimental minus calculated) residual is shown in the lower panel. The current excursion can cover three lines within 1 s: one CH 4 line around 1270.78 cm -1 , a small H 2 O line around 1270.95 cm -1 and one N 2 O line around 1271.08 cm -1 . The retrieved volume mixing ratios are 1.92 ppm of CH 4 and 320 ppb of N 2 O.

N 2 O 4 N 2 O

 242 and CH 4 concentrations have been monitored during a period of 1 day. The mean N 2 O volume mixing ratio was 319 ppb and the mean CH 4 volume mixing ratio was 1.93 ppm, which are typical background troposphere boundary layer values. Over this one-day period, no cycle for N 2 O nor CH 4 was observed. Moreover no correlation between concentrations and weather conditions at ground level were observed.The instrument was further used to measure concentrations of CH 4 and N 2 O when N fertilizer is added to the soil[21]. Added mineral N fertilizer is susceptible to gaseous loss via microbial denitrification, particularly if fertilization is accompanied by irrigation, which creates anaerobic microhabitats suitable for denitrifying bacteria. Denitrification can lead to the release of either N 2 or N 2 O as gaseous end products. Thus the study of the time response of N 2 O fluxes to a fertilizer addition in a soil sample is of high interest[23]. Each soil sample enriched with fertilizer was located in a closed chamber. The box containing the soil sample was then connected to the spectrometer. For this measurement 3.8 g of ammonium nitrate (33% N) were added to the soil sample (mass of 4 kg) that was previously saturated with water (this supply was equivalent to 300 kg of N.ha -1 ). Then N 2 O and CH 4 concentrations were monitored during a period of 9 h after the addition of a commercial fertilizer. The N 2 O volume mixing ratio varies from 320 ppb up to > 800 ppb after 9 h. The CH 4 volume mixing ratio remains constant at #1.93 ppm. These experiments have demonstrated the capabilities of the developed instrument for a long time online atmospheric accurate measurement. detection with a 4.5 µm room-temperature QCL from III-V lab

Eddy covariance measurements around 10 Figure 6 :

 106 Figure 6: Photograph of the N 2 O detector set-up.

Figure 7 :

 7 Figure 7: Retrieved N 2 O concentrations in air with the III-V lab QCL spectrometer

CH 4 and N 2 O detection with a 7.9 µm room-temperature QCL The

  second spectrometer presented in this paper was developed in order to measure CH 4 and N 2 O in the atmosphere. These two gases are well known to be involved in the global warming phenomenon. The global warming potentials of N 2 O and CH 4 at 100 years on the horizon are 296 and 23 times larger than that of CO 2 , respectively[START_REF]Climate Change 2007: the Physical Science Basis, Summary for Policymakers[END_REF]. Worldwide, concentrations of CH 4 have increased by 148% and those of N 2 O by 18% compared to the levels in preindustrial times. Globally, the soil-atmosphere exchange of greenhouse gases is thought to contribute roughly 30% and 70% to the annual emissions of CH 4 and N 2 O, respectively

Thus the final idea of the work is to implement a laser source onboard the SDLA [1-2] balloonborne tunable diode laser spectrometer to investigate the mechanisms that control H 2 O in the tropical tropopause layer.

3