Apostolos Chalkis

Ioannis Z Emiris

Vissarion Fisikopoulos

Panagiotis Repouskos

Elias Tsigaridas

Efficient Sampling from Feasible Sets of SDPs and Volume Approximation

Keywords: spectahedra, semidefinite-programming, sampling, random walk, polynomial eigenvalue problem, volume approximation, optimization

We present algorithmic, complexity, and implementation results on the problem of sampling points from a spectrahedron, that is the feasible region of a semidefinite program.

Our main tool is geometric random walks. We analyze the arithmetic and bit complexity of certain primitive geometric operations that are based on the algebraic properties of spectrahedra and the polynomial eigenvalue problem. This study leads to the implementation of a broad collection of random walks for sampling from spectrahedra that experimentally show faster mixing times than methods currently employed either in theoretical studies or in applications, including the popular family of Hit-and-Run walks. The different random walks offer a variety of advantages, thus allowing us to efficiently sample from general probability distributions, for example the family of log-concave distributions which arise in numerous applications. We focus on two major applications of independent interest: (i) approximate the volume of a spectrahedron, and (ii) compute the expectation of functions coming from robust optimal control.

We exploit efficient linear algebra algorithms and implementations to address the aforementioned computations in very high dimension. In particular, we provide a C++ open source implementation of our methods that scales efficiently, for the first time, up to dimension 200. We illustrate its efficiency on various data sets.

Introduction

Spectrahedra are probably the most well studied shapes after polyhedra. We can represent polyhedra as the intersection of the positive orthant with an affine subspace. Spectrahedra generalize polyhedra, in the sense that they are the intersection of the cone of positive semidefinite matrices -i.e., symmetric matrices with non-negative eigenvalues-with an affine space. In other words, a spectrahedron S ⊂ R n is the feasible set of a linear matrix inequality. That is, let

F(x) = A 0 + x 1 A 1 + • • • + x n A n , (1)
where A i are symmetric matrices in R m×m , then S = {x ∈ R n | F(x) 0}, where denotes positive semidefiniteness. We assume throughout that S is bounded of dimension n. Spectrahedra Figure 1: Left, a 3D elliptope (image from https://en.wikipedia.org/wiki/Spectrahedron); an elliptope -a special case of spectrahedron-is the set of all real, square symmetric matrices whose diagonal entries are all equal to one and whose eigenvalues are all non-negative. Right, a spectrahedron from [START_REF] Breiding | Random spectrahedra[END_REF]; it is called cubic spectrahedron.

are convex sets (Figures 2 and4) and every polytope is a spectrahedron, but not the opposite. They are the feasible regions of semidefinite programs [START_REF] Ramana | Some geometric results in semidefinite programming[END_REF] in the way that polyhedra are feasible regions of linear programs.

Efficient methods for sampling points in spectrahedra are crucial for many applications, such as volume approximation [START_REF] Cousins | Bypassing KLS: Gaussian cooling and an O * (n 3 ) volume algorithm[END_REF], integration [START_REF] Lovasz | Fast algorithms for logconcave functions: Sampling, rounding, integration and optimization[END_REF], semidefinite optimization [START_REF] Lovasz | Fast algorithms for logconcave functions: Sampling, rounding, integration and optimization[END_REF][START_REF] Kalai | Simulated annealing for convex optimization[END_REF], and applications in robust control analysis [START_REF] Calafiore | Robust convex programs: Randomized solutions and applications in control[END_REF][START_REF] Calafiore | Random walks for probabilistic robustness[END_REF][START_REF] Tempo | Randomized algorithms for analysis and control of uncertain systems: with applications[END_REF]. To sample in the interior or on the boundary of S , we employ geometric random walks [START_REF] Vempala | Geometric random walks: A survey[END_REF]. A geometric random walk on S starts at some interior point and at each step moves to a "neighboring" point that we choose according to some distribution, depending only on the current point; thus it is a special category of Markov chains. For example, in the so-called ball walk, we move to a point p that we choose uniformly at random in a ball of fixed radius δ, if p ∈ S . The complexity of a random walk depends on its mixing time -the number of steps required to bound the distance between the current and the stationary distribution-and the complexity of the basic geometric operations performed at each step of the walk; we call the latter per-step complexity.

The majority of geometric random walks are defined for general convex bodies and are based on an oracle; usually the membership oracle. There are also a few walks, e.g., Vaidya walk [START_REF] Chen | A sampling algorithm based on the volumetric barrier[END_REF] and the sub-linear ball walk [START_REF] Mangoubi | Faster polytope rounding, sampling, and volume computation via a sub-linear ball walk[END_REF], specialized for polytopes. Most results on their analysis focus on convergence and mixing time, while they define abstractly the operations they perform at each step and they enclose them in the corresponding oracle. That is why the complexity bounds involve the number of oracle calls.

To specialize a random walk for a family or representation of convex bodies one has to come up with efficient algorithms for the basic geometric operations to realize the (various) oracles. These operations should exploit the underlying geometric and algebraic properties and are of independent interest. They depend on efficient (numerical) linear algebra computations. More importantly, they dominate the per-step complexity and are hence crucial both for the overall complexity to sample a point from the target distribution, as well as for an efficient implementation.

The study of basic geometric operations to sample from non-linear convex objects finds its roots in non-linear computational geometry. During the last two decades, there are combined efforts [START_REF] Boissonnat | Effective computational geometry for curves and surfaces[END_REF][START_REF] Emiris | Nonlinear computational geometry[END_REF] to develop efficient algorithms for the basic operations (predicates) that are behind classical geometric algorithms, like convex hull, arrangements, Voronoi diagrams, to go beyond points and lines and handle curved objects. For this, one exploits efficiently structure and symmetry in linear algebra computations, and develops novel algebraic tools.

To our knowledge, only the Random Directions Hit and Run (W-HnR) random walk [START_REF] Smith | Efficient monte carlo procedures for generating points uniformly distributed over bounded regions[END_REF] has been studied for spectrahedra [START_REF] Calafiore | Random walks for probabilistic robustness[END_REF]. To exploit the various other walks, like Ball walk [START_REF] Vempala | Geometric random walks: A survey[END_REF], Billiard walk (W-Billard) [START_REF] Gryazina | Random sampling: Billiard walk algorithm[END_REF], and Hamiltonian Monte Carlo with boundary reflections (W-HMC-r) [START_REF] Afshar | Reflection, Refraction, and Hamiltonian Monte Carlo[END_REF], one needs to provide geometric operations. For example, we need to compute the reflection of a curve at the boundary, and the intersection point of a curve with the boundary (of a spectrahedron).

We should mention that there is a gap [START_REF] Cousins | A practical volume algorithm[END_REF][START_REF] Betancourt | A conceptual introduction to Hamiltonian Monte Carlo[END_REF] between the theoretical worst case bounds for the mixing times and the practical performance of the random walk algorithms. Furthermore, there are random walks without known theoretical mixing times, such as Coordinate Directions Hit and Run (W-CHnR), W-Billard and W-HMC-r. To study them experimentally, it is imperative to provide an efficient realization of the corresponding oracles.

Previous Work. Sampling convex sets via random walks has attracted a lot of interest in the last decades. Most of the works assume that the convex sets are polytopes; [START_REF] Mangoubi | Faster polytope rounding, sampling, and volume computation via a sub-linear ball walk[END_REF] provides an overview of the state-of-the-art. Random walks on spectrahedra are studied in [START_REF] Polyak | The D-decomposition technique for linear matrix inequalities[END_REF][START_REF] Dabbene | A randomized cutting plane method with probabilistic geometric convergence[END_REF], where they exploit the W-HnR and the computation of the intersection of the walk with the boundary reduces to a generalized eigenvalue problem.

The W-Billard [START_REF] Gryazina | Random sampling: Billiard walk algorithm[END_REF] is a general way of sampling in convex or non-convex shapes from the uniform distribution. A mathematical billiard consists of a domain D and a point-mass, moving freely in D [START_REF] Tabachnikov | Geometry and billiards. Student mathematical library[END_REF]. When this point-mass hits the boundary, it performs a specular reflection, albeit without losing velocity. An application of billiards is the study of optical properties of conics [START_REF] Tabachnikov | Geometry and billiards. Student mathematical library[END_REF]Section 4].

If the trajectory is not a line, but rather a parametric curve, then the intersection operation reduces to solving a polynomial eigenvalue problem (PEP); W-HMC-r requires this operation. PEP is a well-studied problem in computational mathematics, e.g., [START_REF] Tisseur | Backward error and condition of polynomial eigenvalue problems[END_REF], and it appears in many applications. There are important results both for the perturbation analysis of PEP [START_REF] Tisseur | Backward error and condition of polynomial eigenvalue problems[END_REF][START_REF] Berhanu | The polynomial eigenvalue problem[END_REF][START_REF] Dedieu | Perturbation theory for homogeneous polynomial eigenvalue problems[END_REF], as well as for the condition-based analysis of algorithms for the real and complex versions of PEP, if we assume random inputs [START_REF] Armentano | The polynomial eigenvalue problem is well conditioned for random inputs[END_REF][START_REF] Beltrán | The real polynomial eigenvalue problem is well conditioned on the average[END_REF].

For the closely related problem of volume computation, there is also an extensive bibliography [START_REF] Beler | Exact Volume Computation for Polytopes: A Practical Study[END_REF]. The bulk of the theoretical studies are either for general convex bodies [START_REF] Dyer | A random polynomial-time algorithm for approximating the volume of convex bodies[END_REF][START_REF] Cousins | Bypassing KLS: Gaussian cooling and an O * (n 3 ) volume algorithm[END_REF] or polytopes [START_REF] Lee | Convergence Rate of Riemannian Hamiltonian Monte Carlo and Faster Polytope Volume Computation[END_REF]. Practical algorithms and implementations exist only for polytopes [START_REF] Emiris | Practical polytope volume approximation[END_REF][START_REF] Cousins | A practical volume algorithm[END_REF]. Nevertheless, there are notable exceptions that consider algorithms for computing the volume of compact (basic) semi-algebraic sets. For example, in [START_REF] Lairez | Computing the volume of compact semi-algebraic sets[END_REF] they exploit the periods of rational integrals. In the same setting, [START_REF] Korda | Convergence rates of moment-sum-of-squares hierarchies for volume approximation of semialgebraic sets[END_REF][START_REF] Henrion | Approximate volume and integration for basic semialgebraic sets[END_REF] introduce numerical approximation schemes for volume computations, which rely on the moment-based algorithms and semi-definite programming.

Finally, sampling from a multivariate distribution is a central problem in numerous applications. For example, it is useful in robust control analysis [START_REF] Calafiore | Random walks for probabilistic robustness[END_REF][START_REF] Calafiore | Robust convex programs: Randomized solutions and applications in control[END_REF][START_REF] Tempo | Randomized algorithms for analysis and control of uncertain systems: with applications[END_REF] to overcome the worst case hardness as well as in integration [START_REF] Lovasz | Fast algorithms for logconcave functions: Sampling, rounding, integration and optimization[END_REF] and convex optimization [START_REF] Dabbene | A randomized cutting plane method with probabilistic geometric convergence[END_REF][START_REF] Kalai | Simulated annealing for convex optimization[END_REF].

Our contribution. We present a framework of basic geometric operations for computations with spectrahedra. In particular, we analyze the arithmetic and bit complexity of the three fundamental operations, membership, intersection, reflection, by reducing them to linear algebra computations. Based on this framework, we support a rich class of geometric random walks, which in turn we employ in order to build efficient methods for sampling points from spectahedra, under various probability distributions. We apply these tools to approximate the volume of spectrahedra, as well as to integrate over spectrahedral domains. This extends the limits of the state-of-the-art methods that sample from spectrahedra, which actually involve only the W-HnR [START_REF] Dabbene | A randomized cutting plane method with probabilistic geometric convergence[END_REF][START_REF] Calafiore | Random walks for probabilistic robustness[END_REF]. We offer an efficient C++ implementation of our algorithms as a development branch of the package volesti 1 , an open source library for high dimensional sampling and volume computation. While the implementation is in C++, there is also an R interface, for easier access to its functionality. Our implementation is based on state-of-the-art algorithms in numerical linear algebra to address computation in high dimension. Our software makes use of powerful C++ libraries such as Eigen [25] and Spectra [START_REF] Spectra | A header-only C++ library for large scale eigenvalue problems[END_REF].

We demonstrate the efficiency of our approach and implementation on problems from robust control and optimization. First, as a special case of integration, we approximate the volume of spectrahedra up to dimension 100; this is, to the best of our knowledge, the first time such computations for non-linear objects are performed in high dimensions. Then, we approximate the expected value of a function f : R n → [0, 1], whose argument is a random variable having uniform distribution over a spectrahedron of dimension 200.

Finally, we sample from the Boltzmann distribution using W-HMC-r; this exploits a random walk in a spectrahedron that employs a polynomial trajectory of degree two. Sampling using W-HMC-r from truncated distributions is a classical problem in computational statistics [START_REF] Chevallier | Hamiltonian Monte Carlo with boundary reflections, and application to polytope volume calculations[END_REF][START_REF] Afshar | Reflection, Refraction, and Hamiltonian Monte Carlo[END_REF]; alas, existing approaches handle either special distributions or special cases of constraints [START_REF] Pakman | Exact Hamiltonian Monte Carlo for Truncated Multivariate Gaussians[END_REF][START_REF] Lan | Spherical Hamiltonian Monte Carlo for Constrained Target Distributions[END_REF]. We equip W-HMC-r with geometric operations to handle log-concave densities truncated by linear matrix inequalities (LMI) constraints. A combination of Boltzmann distribution with a simulated annealing technique [START_REF] Kalai | Simulated annealing for convex optimization[END_REF] could lead to a practical solver for semidefinite programs (SDP).

A very short version of this paper has appeared as a poster in [START_REF] Chalkis | Sampling the feasible sets of SDPs and volume approximation[END_REF].

Paper organization. First, we introduce our notation. In Section 2 we introduce the basic geometric operations used to efficiently implement membership and boundary oracles. In Section 3, we develop the different types of random walks. Finally, Section 4 presents our implementation, our applications, and various experiments.

Notation. We denote by O, respectively O B , the arithmetic, respectively bit, complexity and we use O, respectively O B , to ignore (poly-)logarithmic factors. The bitsize of a univariate polynomial A ∈ Z[x] is the maximum bitsize of its coefficients. We use bold letters for matrices, A, and vectors, v; we denote by A i, j , respectively v i , their elements; A is the transpose and A * the adjoint of A.

If x = (x 1 , . . . , x n ), then F(x) = A 0 + n i=1 x i A i , see (1)
. For two points x and y, we denote the line through them by (x, y) and their segment as [x, y]. For a spectrahedron S , let its interior be S • and its boundary ∂S . We represent a probability distribution π with a probability density function π(x). When π is truncated to S the support of π(x) is S . If π is log-concave, then π(x) ∝ e -α f (x) , where f : R d → R a convex function. Finally, let B n be the n-dimensional unit ball and denote by ∂B n its boundary.

Basic geometric operations

Our toolbox for computations with spectrahedra and implementing random walks, consists of three basic geometric operations: membership, intersection, and reflection.

For a spectrahedron S , membership decides whether a point lies inside S , intersection computes the intersection of an algebraic curved trajectory C with the boundary ∂S , and reflection computes the reflection of an algebraic curved trajectory when it hits ∂S . The last two operations are required because random walks can move along non-linear trajectories inside convex bodies. For those studied in this paper, the trajectories are parametric polynomial curves, of various degrees. Computation with these curves reduces to solving the polynomial eigenvalue problem (PEP).

Analysis of PEP

To estimate the complexity of intersection we need to bound the complexity of PEP. The Polynomial Eigenvalue Problem (PEP) consists in computing λ ∈ R and x ∈ R m that satisfy the (matrix) equation

P(λ) x = 0 ⇔ (B d λ d + • • • + B 1 λ + B 0 )x = 0 , (2)
where P(λ) is a univariate matrix polynomial whose coefficients are matrices B i ∈ R m×m . We further assume that B d and B 0 are invertible. In general, there are δ = m d values of λ. We refer the reader to [START_REF] Tisseur | Backward error and condition of polynomial eigenvalue problems[END_REF] for a thorough exposition of PEP. One approach for solving PEP is to linearalize the problem and to express the λ's as the eigenvalues of a larger matrix. For this, we transform Equation ( 2) into a linear pencil of dimension δ. Following [START_REF] Berhanu | The polynomial eigenvalue problem[END_REF], the Companion Linearization consists in solving the generalized eigenvalue problem C 0 -λC 1 , where

C 0 =                     B d 0 • • • 0 0 I m . . . . . . . . . . . . . . . 0 0 • • • 0 I m                     and C 1 =                  B d-1 B d-2 • • • B 0 -I m 0 • • • 0 . . . . . . . . . . . . 0 • • • -I m 0                 
, and I m denotes the m × m identity matrix. The eigenvectors x and z are related as follows:

z = [1, λ, . . . , λ d-2 , λ d-1 ] ⊗ x.
To obtain an exact algorithm for PEP we exploit the assumption that B d is invertible so as to transform the problem to the following classical eigenvalue problem (λI d -C 2 )z = 0, where

C 2 =                  B d-1 B -1 d B d-2 B -1 d • • • B 0 B -1 d -I m 0 • • • 0 . . . . . . . . . . . . 0 • • • -I m 0                  .
The eigenvectors are roots of the characteristic polynomial of C 2 . Now the problem is to compute the eigenvalues of C ∈ R δ×δ . From a complexity point of view, the best algorithm relies on computing the roots of the characteristic polynomial [START_REF] Schönhage | The fundamental theorem of algebra in terms of computational complexity[END_REF], which is the method we follow in order to derive complexity bounds. Of course, in practice other methods are faster or more stable.

Lemma 2.1. Consider a PEP of degree d, involving matrices of dimension m × m, with integer elements of bitsize at most τ, see Equation [START_REF] Armentano | The polynomial eigenvalue problem is well conditioned for random inputs[END_REF]. There is a randomized algorithm for computing the eigenvalues and the eigenvectors of PEP up to precision = 2 -L , in O B (δ ω+3 L), where δ = md and L = Ω(δ 3 τ). The arithmetic complexity is O(δ 2.697 + δ lg(1/ )).

Proof. We can compute the characteristic polynomial of an N×N matrix M in O B (N 2.697+1 lg M ) using a randomized algorithm, see [START_REF] Kaltofen | On the complexity of computing determinants[END_REF] and references therein. Here M denotes the largest entry in absolute value. In our case, the elements of C 2 have bitsize O(δτ) and its characteristic polynomial is of degree d and coefficient bitsize O B (δ 2 τ). We compute it in O B (δ 2.697+1 δτ) = O B (δ 4.697 τ) and isolate all its real roots in O B (δ 5 + δ 4 τ) [START_REF] Pan | Univariate polynomials: nearly optimal algorithms for numerical factorization and root-finding[END_REF]; they correspond to the real eigenvalues of PEP. We can decrease the width of the isolating interval by a factor of = 2 -L for all the roots in O B (δ 3 τ + δL) [START_REF] Pan | Nearly optimal refinement of real roots of a univariate polynomial[END_REF]. Thus, the overall complexity is O B (δ 5 + δ 4.697 τ + δL).

It remains to compute the corresponding eigenvectors. For each eigenvalue λ we may compute the corresponding eigenvector z by Gaussian elimination and back substitution to the (augmented) matrix [λI δ -C 2 | 0]. This requires O(δ ω ) arithmetic operations. However, as λ is a root of the characteristic polynomial, one has to operate on algebraic numbers, which is highly nontrivial, and one needs to bound the number of bits needed to compute the elements of z correctly and to recover x.

Hence, we employ separation bounds for polynomial system adapted to eigenvector computation [START_REF] Emiris | Separation bounds for polynomial systems[END_REF]. One needs, as in the case of eigenvalues, O B (δ 4 + δ 3 τ) bits to isolate the coordinates of the eigenvectors. To decrease the width of the corresponding isolating intervals by a factor of = 2 -L , the number of bits becomes O B (δ 4 + δ 3 τ + L). Thus, we compute the eigenvectors in

O B (δ ω (δ 4 + δ 3 τ + L)) = O B (δ ω+4 + δ ω+3 τ + δ ω L).
For the arithmetic complexity we proceed as follows: We compute the characteristic polynomial in O(δ 2.697 ), we approximate its roots up to in O(δ lg(1/ )). Finally, we compute the eigenvectors with O(δ ω ) arithmetic operations. So the overall cost is O(δ 2.697 + δ lg(1/ )).

membership

The operation membership(F, p) decides whether a point p lies in the interior of a spectrahedron S = {x ∈ R n | F(x) 0}. For this, first, we construct the matrix F( p). If it is positive definite, then p ∈ S • . If it is positive semidefinite, then p ∈ ∂S , otherwise p ∈ R n \ S . The pseudo-code appears in Algorithm 1.

Algorithm 1: membership(F, p)

Input : An LMI F(x) 0 representing a spectrahedron S and a point p ∈ R n . Output: true if p ∈ S , false otherwise.

O B ((nm 2 + m 3.697 )(τ + σ)).
Proof. We construct F( p) in O(nm 2 ). Then, with O(m 2.697 ) operations we compute its characteristic polynomial [START_REF] Kaltofen | On the complexity of computing determinants[END_REF] and in O(m) we decide whether it has negative roots, for example by solving [START_REF] Pan | Univariate polynomials: nearly optimal algorithms for numerical factorization and root-finding[END_REF] or using fast subresultant algorithms [START_REF] Lecerf | On the complexity of the Lickteig-Roy subresultant algorithm[END_REF][START_REF] Lickteig | Sylvester-habicht sequences and fast cauchy index computation[END_REF]. For the bit complexity, the construction costs O B (nm 2 (τ + σ)) and the computation of the characteristic polynomial costs O B (m 2.697+1 (τ + σ)), using a randomized algorithm [START_REF] Kaltofen | On the complexity of computing determinants[END_REF]. One may test for negative roots, and thus eigenvalues, in O B (m 2 n(τ + σ)) [START_REF] Lickteig | Sylvester-habicht sequences and fast cauchy index computation[END_REF].

intersection

Consider a parametric polynomial curve C such that it has a non-empty intersection with a spectrahedron S . Throughout, we consider only the real trace of C. Assume that the value of the parameter t = 0 corresponds to a point p 0 , that lies in C ∩ S • . Furthermore, assume that the segment of C on which p 0 lies, intersects the boundary of S transversally at two points, say p - and p + . The operation intersection computes the parameters, t -and t + , corresponding to these two points. Figure 2 illustrates this discussion and the pseudo-code of intersection appears in Algorithm 2.

To prove correctness and estimate the complexity we proceed as follows: As before (see also Equation 1), S is the feasible region of linear matrix inequalities (LMI) F(x) 0. Consider the real trace of a polynomial curve C, with parameterization

Φ : R → R n t → Φ(t) := (p 1 (t), . . . , p n (t)), (3)
where p i (t) = d i j=0 p i, j t j are univariate polynomials in t of degree

d i , for i ∈ [m]. Also let d = max i∈[m] {d i }.
If the coefficients of the polynomials are integers, then we further assume that the maximum coefficient's bitsize is bounded by τ.

As t varies over the real line, there may be several disjoint intervals, for which the corresponding segment of C lies in S • . We aim to compute the endpoints, t -and t + , of a maximal such interval containing t = 0. Let p 0 = Φ(0); by assumption it holds F(Φ(0)) = F( p 0 ) 0.

The input of intersection (Algorithm 2) is F, the LMI representation of S , and Φ(t), the polynomial parameterization of C. Its crux is a routine, pep, that solves a polynomial eigenvalue problem. The following lemma exploits this relation.

Lemma 2.3 (pep and S ∩ C). Consider the spectrahedron S

= {x ∈ R n | F(x)
0}. Let Φ : R → R n be a parameterization of a polynomial curve C, such that Φ(0) ∈ S • . Let [t -, t + ] be the maximal interval containing 0, such that the corresponding part of C lies in S . Then, t -, respectively t + , is the maximum negative, respectively minimum positive, polynomial eigenvalue of F(Φ(t))x = 0, where

F(Φ(t)) = B 0 + tB 1 + • • • + t d B d .
Algorithm 2: intersection(F, Φ(t))

Input : An LMI F(x) 0 for a spectrahedron S and a parameterization Φ(t) of a polynomial curve C Require:

Φ(0) ∈ S • Output : t -, t + s.t. Φ(t -), Φ(t + ) ∈ ∂S 1 T := {t 1 ≤ t 2 ≤ • • • ≤ t } ← pep(F(Φ(t))); 2 t -← max{t ∈ T | t < 0}; // max negative polynomial eigenvalue 3 t + ← min{t ∈ T | t > 0}; // min positive polynomial eigenvalue 4 return t -, t + ;
Proof. The composition of F(x) and Φ(t) gives

F(Φ(t)) = A 0 + p 1 (t) A 1 + • • • + p n (t) A n . (4)
We rewrite Expression (4) by grouping the coefficients for each t i , i ∈ [d], then

F(Φ(t)) = B 0 + tB 1 + • • • + t d B d , (5)
where

B k = n j=0 p j,k A j , for 0 ≤ k ≤ d.
We use the convention that p j,k = 0, when k > d j . For t = 0, it holds, by assumption, that F(Φ(0)) = B 0 0: point Φ(0) lies in the interior of S . Actually, for any x ∈ S • it holds F(x) 0. On the other hand, if x ∈ ∂S , then F(x) 0. Our goal is to compute the maximal interval [t -, t + ] that contains 0 and for every t in it, we have F(Φ(t)) 0.

Starting from point Φ(0), by varying t, we move on the trajectory that C defines (in both directions) until we hit the boundary of S . When we hit ∂S , matrix F(Φ(t)) is not strictly definite anymore. Thus, its determinant vanishes.

Consider the function θ : R → R, where θ(t) = det F(Φ(t)) is a univariate polynomial in t. If a point Φ(t) is on the boundary of the spectrahedron, then θ(t) = 0. We opt to compute t -and t + , such that t -≤ 0 ≤ t + and θ(t -) = θ(t + ) = 0. At t = 0, θ(0) > 0 and the graph of θ is above the t-axis. So C intersects the boundary when the graph of θ touches the t-axis for the first time at t 1 ≤ 0 ≤ t 2 . It follows that t -= t 1 and t + = t 2 are the maximum negative and minimum positive roots of θ, or equivalently the corresponding polynomial eigenvalues of F(Φ(t)).

Lemma 2.4. Algorithm 2, intersection(F, Φ(t)), uses O((md) 2.697 + md lg L) arithmetic operations in order to compute the intersection, up to precision = 2 -L , of an LMI, F, consisting of n matrices of dimension m × m with a parametric curve, Φ(t), of degree d.

Proof. We have to construct PEP and solve it. Since Φ(t) has degree d, then

F(Φ(t)) = B 0 + tB 1 + • • • + t d B d is a PEP of degree d.
This construction costs O(dnm 2 ) operations. The solving phase, using Lemma 2.1, requires O((md) 2.697 + md lg L) arithmetic operations and dominates the complexity bound of the algorithm.

reflection

The reflection operation (Algorithm 3) takes as input an LMI representation, F, of a spectrahedron S and a polynomial curve C, given by a parameterization Φ. Assume that t = 0 corresponds to a point Φ(0) ∈ S • ∩ C. Starting from t = 0, we increase t along the positive real Algorithm 3: reflection (F, Φ(t))

Input : An LMI F(x) 0 for a spectrahedron S and a parameterization Φ(t) of a polynomial curve C. Require: (i) Φ(0) ∈ S • (ii) C intersects ∂S transversally at a smooth point. Output : t + such that Φ(t + ) ∈ ∂S and the direction of the reflection, s + , at this point.

1 t -, t + ← intersection (F, Φ(t)); 2 w ← ∇ det F(Φ(t + )); 3 w ← w w ; // Normalize w 4 s + ← dΦ dt (t + ) -2 ∇ dΦ dt (t + ), w w; 5 return t + , s + ;
semi-axis. As t changes, we move along the curve C through Φ(t), until we hit the boundary of S at the point p + := Φ(t + ) ∈ ∂S , for some t + > 0. Then, a specular reflection occurs at this point with direction s + ; this is the reflected direction. We output t + and s + . Figure 2 depicts the procedure.

The boundary of S , ∂S , with respect to the Euclidean topology, is a subset of the real algebraic set {x ∈ R n | det(F(x)) = 0}. The latter is a real hypersurface defined by one (determinantal) equation. For any x ∈ ∂S we have rank(F(x)) ≤ m -1. We assume that p + = Φ(t + ) is such that rank(F( p + )) = m -1. The normal direction at a point p ∈ ∂S , is the gradient of det F( p).

We compute the reflected direction using the following formula

s + = u -2 |w| 2 u, w w, ( 6
)
where w is the normalized gradient vector at the point Φ(t + ) and u = dΦ dt (t + ) is the direction of the trajectory at this point. We illustrate the various vectors in Figure 2.

m matrix F(x) is m -1. Then ∇ det(F(x)) = c • (v A 1 v, • • • , v A n v), (7)
where c = µ(F(x)) |v| 2 , µ(F(x)) is the product of the nonzero eigenvalues of F(x), and v is a non-trivial vector in the kernel of F(x). If rank(F(x)) ≤ m -2, then the gradient is the zero.

Proof. From Lemma Appendix A.2:

∂ det F(x) ∂x k = Trace (F(x) * A k ) . (8)
If rank(F(x)) ≤ -2, then F(x) * is the zero matrix. If we assume that rank(F(x)) = m -1, then using Lemma Appendix A.3:

Trace (F(x) * A k ) = Trace µ(F(x)) vu u v A k = µ(F(x)) u v • Trace vu A k = µ(F(x)) u v • u A k v.
However, since F(x) is symmetric, we can choose v = u, so:

µ(F(x)) u v • u A k v = µ(F(x)) |v| 2 • v A k v.
The algorithm reflection exploits Lemma 2.5. Nevertheless, it is not necessary to perform all then computations that the lemma indicates. For example, because we will normalize the resulting vector and we do not need its actual direction (internal or external), we can omit the computation of c. Moreover, the nonzero vector v, which satisfies F( p)v = 0, corresponds to the eigenvector w.r.t. the eigenvalue t + from the pep (Lemma 2.3). This is true because p = Φ(t + ) ∈ ∂S and thus det F(Φ(t + )) = 0.

At this point we should note that we compute the eigenvalues of pep up to some precision. Since the matrix-vector multiplication is backward stable, a small perturbation on v does not affect the computation of each coordinate of ∇ det(F(x)) [55, p. 104]. We quantify the accuracy of the computed ∇ det(F(x)) using floating point arithmetic as follows:

Lemma 2.6. The relative error of each coordinate of the gradient given in Lemma 2.5 when we compute it using floating point arithmetic with machine epsilon M is O( M σ max (A i ) ), for i ∈ [n], where σ max is the largest singular value of A i .

Proof. Let A ∈ R m×m be a symmetric matrix and consider the map f : v → v T Av. The relative condition number of f as defined in [55, p. 90] is

k(v) = ||J(v)|| || f (v)||/||v|| = 2 ||Av|| v T A i v = 2 σ max (A) σ 2 max (A) = 2 σ max (A)
,

where J(•) is the Jacobian of f . According to Theorem 15.1 in [55, p. 111], since the matrixvector multiplication is backward stable, the relative error of each coordinate in the gradient computation of Lemma 2.5 is O( M σ max (A i ) ), for 1 ≤ i ≤ n. Lemma 2.7. Let S be a spectrahedron represented by an LMI, F(x), consisting of n matrices of dimension m × m. Also let C be a parametric curve with parameterization Φ(t), involving polynomials of degree at most d. Algorithm 3, reflection(F, Φ(t)), computes the intersection, up to precision = 2 -L , of S with C, and the reflection of C at ∂S , by performing O((md) 2.697 + md lg L + dnm 2 ) arithmetic operations.

Proof. By inspecting Algorithm 3 we notice that the complexity of the algorithm depends on the construction of ∇ det(F(x)) and the call to intersection.

To compute ∇ det(F(x)) we just need to compute (v

A 1 v, • • • , v A n v).
If we have already computed v, then this computation requires O(nm 2 ) operations. The computation of the derivative of Φ(t) is straightforward, as Φ is a univariate polynomial. Taking into account the complexity of intersection, the total complexity for reflection is O((md) 2.697 + md lg L + dnm 2 + nm 2 ) = O((md) 2.697 + md lg L + dnm 2 ).

+ m lg L + nm 2 ) W-CHnR O(m 2.697 + m lg L + m 2 ) W-Billard O(ρ(m 2.697 + m lg L + nm 2 )) W-HMC-r O(ρ((dm) 2.697 + md lg L + dnm 2 ))
Table 1: The per-step complexity of the random walks in Section 3.

An example in 2D

Consider a spectrahedron S in the plane (Figure 3), given by an LMI F(x) = A 0 +x 1 A 1 +x 2 A 2 . The matrices A i , 0 ≤ i ≤ 2, are in the appendix.

Starting from the point p 0 = (-1, 1) , we walk along the line L with parameterization: Φ(t) = p 0 + tu, where u = (1.3, 0.8) . Then, intersection finds the intersection of S with L, by solving the degree one PEP, (B 0 + tB 1 )x = 0, where B 0 = F( p 0 ) and B 1 = u 1 A 1 + u 2 A 2 . Acquiring t -= -0.8 and t + = 0.5, we obtain the intersection point p 1 , which corresponds to p 0 + t + u = (-0.3, 1.4) .

To compute the direction of the trajectory, immediately after we reflect on the boundary of S at p 1 , reflection computes

w = ∇ det F(Φ(t + )) |∇ det F(Φ(t + ))| = (v A 1 v, v A 2 v) = (-0.2, -1) , ( 9
)
where v is the eigenvector of (B 0 + tB 1 )x = 0, with eigenvalue t + . The reflected direction is u = u -2 u, w w = (0.8, -1.3) .

Random walks

Using the basic geometric operations of Section 2, we implement and analyze three random walks for spectrahedra: Hit and Run (W-HnR), its variant Coordinate Directions Hit and Run (W-CHnR), Billiard Walk (W-Billard), and Hamiltonian Monte Carlo with reflections (W-HMC-r).

In Table 1, we present the per-step arithmetic complexity for each random walk.

Hit and Run

W-HnR (Algorithm 4) is a random walk that samples from any probability distribution π truncated to a convex body K; in our case a spectrahedron S . We should mention that there exist bounds for its mixing time only when π is log-concave distribution, for example the uniform distribution, which is O(n 3 ).

At the i-th step, W-HnR chooses uniformly at random a (direction of a) line , passing from its current position p i . Let p 1 and p 2 be the intersection points of with S . Let π be the restriction of π on the segment [ p 2 , p 2 ]. Then, we choose p i+1 from [ p 1 , p 2 ] w.r.t. the distribution π .

Algorithm 4: Hit-and-Run Walk (W-HnR)

Input : LMI F(x) 0 for a spectrahedron S & a point p i . Require:

p i ∈ S Output :
The point p i+1 of the (i + 1)-th step of the walk.

1 BO (F, interior point p i ) v ← R U(∂B n ); // choose direction 2 Φ(t) := p i + tv; // define trajectory 3 t -, t + ← intersection (F, Φ(t)); 4 p i+1 ← R [ p i + t -v, p i + t + v] w.r.t. π ; 5 return p i+1 ; Lemma 3.1. The per-step complexity of W-HnR is O(m 2.697 + m lg 1/ + nm 2 ).
Proof. The per-step complexity of W-HnR is dominated by the intersection, which requires O(nm 2 ) operations for the construction of the pep and O(m 2.697 +m lg 1/ ) for solving it; in the case where we want to approximate the intersection point up to a factor or = 2 -L (Lemma 2.4).

There is also a variation of W-HnR, the coordinate directions Hit and Run (W-CHnR) [START_REF] Smith | Efficient monte carlo procedures for generating points uniformly distributed over bounded regions[END_REF]. This walk chooses the direction vector randomly and uniformly among the vector basis {e i , i ∈ [n]}. In W-CHnR, for every step aside the first, the construction of the pep takes O(m 2 ) operations, and the complexity does not depend on the dimension n. The reason for this improvement is that to build the pep, F( p i + te j ) = F( p i ) + t A j , we can obtain the value of F( p i ) from F( p i ) = F( p i-1 ) + t A k , assuming that at the previous step we have chosed e k as direction. There is no theoretical mixing time for W-CHnR.

Billiard walk

W-Billard [START_REF] Polyak | Billiard walk -a new sampling algorithm for control and optimization[END_REF], Algorithm 5, samples a convex body K under the uniform distribution; no theoretical results for its mixing time exist. At i-th step, being at position p i , it chooses uniformly a direction vector v and a number , where = -τ ln η and η ∼ U(0, 1). Then, it moves at the direction of v for distance . If during the movement, it hits the boundary without having covered the required distance , then it continues on a reflected trajectory. If the number of reflections exceeds a bound ρ, it stays at p i . In [START_REF] Polyak | Billiard walk -a new sampling algorithm for control and optimization[END_REF] they experimentally conclude that W-Billard mixes faster when τ ≈ diam(K), where diam(K) is the diameter of K.

Lemma 3.2. The per-step complexity of W-Billard is O(ρ(m 2.697 + m lg L + nm 2 )), where ρ is the number of reflections. where κ is the condition number of ∇ 2 f [START_REF] Lee | Algorithmic theory of odes and sampling from well-conditioned logconcave densities[END_REF]. If we truncate π by considering its restriction in a convex body, then we can use boundary reflections (W-HMC-r), as in Algorithm 6, to ensure that the random walk converges to the target distribution [START_REF] Chevallier | Hamiltonian Monte Carlo with boundary reflections, and application to polytope volume calculations[END_REF]; however, in this case the mixing time is unclear.

The Hamiltonian dynamics behind HMC operate on a n-dimensional position vector p and a n-dimensional momentum v. So the full state space has 2n dimensions. The system is described by a function of p and v known as the Hamiltonian,

H( p, v) = U( p) + K(v) = f ( p) + 1 2 |v| 2 .
To sample from π, one has to solve the following system of Ordinary Differential Equations (ODE):

d p dt = ∂H( p, v) ∂v dv dt = - ∂H( p, v) ∂p ⇒            d p(t) dt = v(t) dv(t) dt = -α∇ f (p) . (10)
If π(x) is a log-concave density, then we can approximate the solution of the ODE with a low degree polynomial trajectory [START_REF] Lee | Algorithmic theory of odes and sampling from well-conditioned logconcave densities[END_REF], using the collocation method. A degree d = O(1/ log( )) suffices to obtain a polynomial trajectory with error O( ), for a fixed time interval, while we perform just O(1) evaluations of ∇ f (x). HMC at the i-th step uniformly samples a step from a proper interval to move on the trajectory implied by ODE [START_REF] Calafiore | Random walks for probabilistic robustness[END_REF], choses v randomly from N(0, I n ), and updates p using the ODE in [START_REF] Calafiore | Random walks for probabilistic robustness[END_REF], for t ∈ [0, ]. When π is truncated in a convex body, then W-HMC-r fixes an upper bound ρ on the number of reflections and reflects a polynomial trajectory as we describe in Section 2.4.

Each step of W-HMC-r, when π(x) is a log-concave density truncated by S , costs O(ρ((dm) 2.697 + md lg L + dnm 2 )), if we approximate the intersection points up to a factor 2 -L , where d is the degree of the polynomial that approximates the solution of the ODE [START_REF] Calafiore | Random walks for probabilistic robustness[END_REF].

Lemma 3.3. The per-step complexity of W-HMC-r is O(ρ((dm) 2.697 + md lg L + dnm 2 )), where ρ is the number of reflections.

Applications and experiments

This section demonstrates and compares the algorithms of Section 3 and the efficiency of our software on three applications that rely on sampling from spectrahedra.

We call walk length the number of the intermediate points that a random walk visits before producing a single sample. The longer the walk length of a random walk is, the smaller the distance of the current distribution to the stationary (target) distribution becomes. Typically we choose a sufficiently large length for the first sample, this procedure is often called "burning".

Our code is parameterized by the floating point precision of the computations. We use Eigen [25] for basic linear algebra operations, such as Cholesky decomposition and matrix multiplication. For eigenvalue computations, we employ Spectra [START_REF] Spectra | A header-only C++ library for large scale eigenvalue problems[END_REF], which is based on Eigen and offers crucial optimizations. First, it solves generalized eigenvalue problems of special structure; that is (B 0 -λB 1 )v = 0, when B 0 is positive semidefinite and B 1 symmetric. This operation is encountered when W-Billard or W-HnR call intersection. Second, it offers directly the computation of the largest eigenvalue which corresponds to t + after a simple transformation. Finally, Spectra provides ∼ 20x speedup over the default eigenvalue computation by Eigen. To the best of our knoledge, our software is the first that can sample efficiently from spectrahedra and estimates volumes up to a few hundred dimensions. It is accessible on github. 2For our experiments, we generate random spectrahedra following [START_REF] Dabbene | A randomized cutting plane method with probabilistic geometric convergence[END_REF]. In particular, to construct the LMI of Equation ( 1) we set A 0 to be positive semidefinite, i.e., A 0 = ZZ T + I m , where we pick the elements of Z ∈ R m×m uniformly at random from [0, 1]. Then, for A i , i = 1, . . . , n we set,

A i = Q 0 0 -Q , Q = Q + Q T ,
where we pick the elements of Q ∈ R (m/2)×(m/2) uniformly at random from [-1, 1]. We performed all the experiments on PC with Intel Core i7-6700 3.40GHz × 8 CPU and 32GB RAM.

Volume computation

We use the geometric operations (Section 2) and the random walks (Section 3) to compute the volume vol(S ) of spectrahedron S . Our implementation approximates vol(S ) within relative error 0.1 with high probability in a few minutes, for dimension n = 100.

A typical randomized algorithm for volume approximation exploits a Multiphase Monte Carlo (MMC) technique, which reduces volume approximation of convex body S to computing a telescoping product of ratios of integrals over S . In particular, for any sequence of functions { f 0 , . . . , f k }, where f i : R n → R, we have:

vol(S ) = S 1dx = S f k (x)dx S f k-1 (x)dx S f k (x)dx • • • S 1dx S f 0 (x)dx . ( 11
)
Notice that P f i-1 (x)dx

P f i (x)dx = P f i-1 (x) f i (x) f i (x) P f i (x)dx dx.
To estimate each ratio of integrals, we sample N points from a distribution proportional to f i and, we use the unbiased estimator 1 N N j=1 f i-1 (x j ) f i (x j ) . To exploit Equation [START_REF] Calafiore | Robust convex programs: Randomized solutions and applications in control[END_REF] we have to (i) fix the sequence such that k is as small as possible, (ii) select f i 's such that we can compute efficiently each integral ratio, and (iii) compute P f k (x)dx. The best theoretical result of [START_REF] Cousins | Bypassing KLS: Gaussian cooling and an O * (n 3 ) volume algorithm[END_REF] fixes a sequence of spherical Gaussians { f 0 , . . . , f k } with the mode being in S , parameterized by the variance. The overall complexity is O(n 3 ) membership calls. The implementation in [START_REF] Cousins | A practical volume algorithm[END_REF] is based on this algorithm but handles only convex polytopes in Hrepresentation as it requires the facets of the polytope and an inscribed ball to fix the sequence of Gaussians. Both the radius of the inscribed ball and the number of facets strongly influence the performance of the algorithm. So, it cannot handle efficiently the case of convex bodies without a facet description, e.g., zonotopes [START_REF] Cousins | A practical volume algorithm[END_REF], as it results a big sequence of ratios that spoil practical efficiency.

Our approach is to consider the f i 's as a sequence of indicator functions of concentric balls centered in S , as in [START_REF] Emiris | Practical polytope volume approximation[END_REF]. In particular, let f k and f 0 be the indicator functions of rB n and RB n respectively, while rB n ⊆ S ⊆ RB n and S i = (2 (k-i)/n rB n ) ∩ S for i = 0, . . . , k. Thus, it suffices to compute vol(rB n ) and apply the following:

vol(S ) = vol(S k ) vol(S k-1 ) vol(S k ) • • • vol(S 0 ) vol(S 1 ) , k = n lg(R/r) . (12)
Furthermore, we employ the annealing schedule from [START_REF] Chalkis | Practical volume estimation by a new annealing schedule for cooling convex bodies[END_REF] to minimize k, without computing neither an enclosed ball rB n nor an enclosing ball RB n of S . We do so by probabilistically bounding each ratio of Equation ( 12) in an interval [r, r + δ], which is given as input. To approximate each ratio of volumes, we sample uniformly distributed points from S i and count points in S i-1 . We follow the experimental results of Section 4.2 and use W-Billard which mixes faster than W-HnR. Table 2 reports the average volume, runtime, number of points generated for each S -n-m over 10 trials. We also compute a 95% confidence interval for the volume. Notice that for all cases the extreme values of each interval imply an error ≤ 0.1, which was the requested error. For n = 40 just a few seconds suffice to approximate the volume and for n = 100 our implementation takes a few minutes.

Expected value of a function

Randomized algorithms are commonly used for problems in robust control analysis to overcome the (worst case) hardness, especially in probabilistic robustness [START_REF] Blondel | A survey of computational complexity results in systems and control[END_REF][START_REF] Khargonekar | Randomized algorithms for robust control analysis and synthesis have polynomial complexity[END_REF][START_REF] Tempo | Randomized algorithms for analysis and control of uncertain systems: with applications[END_REF][START_REF] Calafiore | Robust convex programs: Randomized solutions and applications in control[END_REF]. A central problem is to approximate the integral of a function over a spectrahedron, e.g. [START_REF] Calafiore | Random walks for probabilistic robustness[END_REF][START_REF] Ray | A monte carlo approach to the analysis of control system robustness[END_REF] and thus uniform sampling is of particular interest. To put our experiments into perspective, we present experiment up to n = 200, while in [START_REF] Calafiore | Random walks for probabilistic robustness[END_REF] and [START_REF] Calafiore | Robust convex programs: Randomized solutions and applications in control[END_REF] they use only W-HnR for experiments in n ≤ 10. Our goal is to compute the expected value of a function f : R n → [0, 1], with respect to the measure given by the uniform distribution π over S , i.e., I = S f (x)π(x)dx. A standard approach is the Monte Carlo method, which suggests to sample N independent samples from π. Then,

ÊN [ f ] = 1 N N i=1 f (x i )
is an unbiased estimator for I. We employ the random walks of Section 3 to sample uniformly distributed points from S (i.e., W-HnR, W-CHnR, and W-Billard) and we experimentally compare their efficiency. It turns out that W-Billard mixes much faster and results to better accuracy (see Figures 5 & 6). This observation agrees with the experiments on the rate of convergence for W-HnR and W-Billard in [START_REF] Polyak | Billiard walk -a new sampling algorithm for control and optimization[END_REF]. To come to a decisive conclusion we need to perform a more detailed practical study on the mixing time of these random walks; we leave this study as future work.

The variance of an estimator is a crucial as it bounds the approximation error. Using Chebyshev's inequality and [START_REF] Lovász | Random walks and an O * (n 5 ) volume algorithm for convex bodies[END_REF], we have

Prob[| ÊN [ f ] -E[ f ]| ≤ ] ≤ var( ÊN [ f ]) 2 ≤ 4M N 2 , ( 13
)
where M is the mixing time of the random walk one uses to sample " close" to the uniform distribution from S . Thus, for fixed N and , the smaller the mixing time of the random walk is, the smaller the variance of estimator ÊN [ f ] and hence, the better the approximation. We estimate I when f := 1(S ∩ H), where H is the union of two half-spaces

H := {x | cx ≤ b 1 or cx ≥ b 2 },
where b 2 > b 1 and 1(•) is the indicator function. Note that I = vol(S ∩ H)/vol(S ). To estimate it we sample approximate uniformly distributed points from S and we count the number of points that lie in S ∩ H. & 6 illustrates these values, while the walk length increases. Notice that the st.d. is much smaller and the approximation more stable when W-Billard is used compared to both W-HnR and W-CHnR. As W-Billard mixes faster, we report in Table 3 the average time our software needs to sample N = 200 points for various walk lengths for W-Billard in n = 100, 200. The average time to generate a point is ≈ 0.3 and ≈ 7.2 milliseconds respectively.

Sampling from non-uniform distributions

The random walks of Section 3 open a promising avenue for approximating the optimal solution of a semidefinite program, that is min c, x , subject to x ∈ S .

(

) 14
We parameterize the optimization algorithm in [START_REF] Kalai | Simulated annealing for convex optimization[END_REF] with the choice of random walk and demonstrate that its efficiency relies heavily on the sampling method. We perform experiments with W-HMC-r and W-HnR, as both can sample from the distribution the algorithm requires. Deterministic approximations to the optimal solutions of these tests, were acquired via the SDPA library [START_REF] Yamashita | Implementation and evaluation of sdpa 6.0[END_REF]. The strategy to approximate the optimal solution x * of Equation ( 14), is based on sampling from the Boltzmann distribution, i.e., π(x) ∝ e -cx/T , truncated to S . The scalar T , is called temperature. As the temperature T diminishes, the mass of π tends to concentrate around its mode, which is x * . Thus, one could obtain a uniform point using the algorithm in [START_REF] Lovász | Simulated annealing in convex bodies and an O * (n 4 ) volume algorithm[END_REF], and then use it as a starting point to sample from π 0 ∝ e -cx/T 0 , where T 0 = R and S ⊆ RB n . Then, the cooling schedule T i+1 = T i (1 -1/ √ n) guarantees that a sample from π i yields a good starting point for π i+1 . After O( √ n) steps the temperature will be low enough, to sample a point within distance from x * with high probability. In [START_REF] Kalai | Simulated annealing for convex optimization[END_REF], they use only W-HnR. We also employ W-HMC-r. To sample from Boltzmann distributions with W-HMC-r, at each step, starting from p i and with momenta v i , the ODE of Equation ( 10

Its solution is the polynomial p(t) = -c 2T t 2 + v i t + p i , which is a parametric representation of a polynomial curve, see Equation [START_REF] Beltrán | The real polynomial eigenvalue problem is well conditioned on the average[END_REF].

In Table 4 we follow the cooling schedule described, after setting T 0 ≈ R and sampling the first uniform point with W-Billard. We give the optimal solution as input and we stop dropping T when an error ≤ 0.05 is achieved. Even in the case when the walk length is set equal to one, W-HMC-r still converges to the optimal solution. To the best of our knowledge, this is the first time that a randomized algorithm, which is based on random walks, is functional even when the walk length is set to one. On the other hand, we set the walk length of W-HnR O( √ n) or O(n) in our experiments. Notice that for the smaller walk length, its runtime decreases, but the method becomes unstable, as it sometimes fails to converge. For both cases its runtime is worse than that of W-HMC-r.

Lemma Appendix A.3 (Adjoint Matrix of A). Let A be a m × m matrix of rank r(A) = m -1. Then A * = µ(A) vu u v where µ( A) is the product of the m -1 non-zero eigenvalues of A, and x and y satisfy Av = A u = 0 (see chapter 3.2 in [START_REF] Magnus | Matrix Differential Calculus with Applications in Statistics and Econometrics (Revised Edition[END_REF]).

Appendix B. Matrices of the Example

The spectrahedron was randomly generated as in [START_REF] Dabbene | A randomized cutting plane method with probabilistic geometric convergence[END_REF]. Due to space considerations, the entries of the matrices are rounded to the first decimal.

A 0 =                         
                         (B.1) A 1 =                          0.5 -0.
                         (B.2) A 2 =                          2
                         (B.3)

Figure 2 :

 2 Figure 2: A spectrahedron S described by F(x) and a parameterized curve Φ. The point p 0 = Φ(0) lies in the interior of S , and the points p + = Φ(t + ) and p -= Φ(t -) on the boundary. Vector w is the surface normal of ∂S at p + and u is the direction of Φ at time t = t + .

Lemma 2 . 5 (

 25 Gradient of det F(x)). Assume that x ∈ ∂S and the rank of the m ×

Figure 3 :

 3 Figure 3: The i-th step of the W-Billard [Algorithm 5] (left) and of the W-HMC-r [Algorithm 6] (right) random walks.

Figure 4 :

 4 Figure 4: Samples from the uniform distribution with W-Billard (left) and from the Boltzmann distribution π(x) ∝ e -cx/T , where T = 1, c = [-0.09, 1] T , with W-HMC-r (right). The volume of this spectrahedron is 10.23.

4 S 3 Table 2 :

 432 -100-100 (1.21 ± 0.10)e-51 52962.7 362.(Volume of spectrahedra) For each S -n-m we run ν = 10 experiments; m is the matrix dimension in LMI and n the ambient dimension. µ stands for the average volume, s for the standard deviation; we give a confidence interval with level of confidence α = 0.05; t α,ν-1 is the critical value of student's distribution with ν -1 degrees of freedom. Points denotes the average number of points generated and T ime the average runtime in seconds. For all the above we set the error parameter e = 0.1.

Figure 5 :

 5 Figure 5: The standard deviation of ÊN [ f ] over M = 20 trials, estimating 2 functions with E[ f 1 ] = 0.0993 and E[ f 2 ] = 0.588. For each walk length we sample N = 200 points and we repeat M = 20 times.

Figure 6 :Table 3 :

 63 Figure 6: The mean value of the estimator ÊN [ f ] over M = 20 trials, estimating two functions with E[ f 1 ] = 0.0993 and E[ f 2 ] = 0.588. For each walk length we sample N = 200 points and we repeat M = 20 times. walk length 1 5 10 20 30 40 50 S -100-100 1.4 3.2 7.7 9.5 16.1 21.4 28.2 S -200-200 16.2 75.7 148 303 443 584 722 Table 3: Average time in sec to sample 200 points with W-Billard from 10 random spectrahedra S -n-m; n for the dimension that S -n-m lies; m for the dimension of the matrix in LMI.

Figure 7 :

 7 Figure 7: Sample a point from π(x) ∝ e -cx/T i and update the objective current best in each iteration, with T 0 ≈ diam(S ) and T i = T i-1 (1 -1/ √ n), i = 1, . . . 70. The walk length equals to one for W-HMC-r and 500 + 4n 2 = 10 500 for W-HnR.

) becomes d 2

 2 dt 2 p(t) = -c T , d dt p(0) = v i , p(0) = p i .

  -1 16.7 12.3 5.4 18.7

	16.7 3.7 12.3 8.7	5.1 10.4
	3.7	9.4	2.3	4	-2.3 -1
	12.3 2.3 26.8 18.7 7.1 16.7
	8.7	4	18.7 20	3.7 12.3
	5.1 -2.3 7.1 3.7	6.1	5.4
	10.4

https://github.com/GeomScale/volume_approximation/tree/sample_spectrahedra

https://github.com/GeomScale/volume_approximation/tree/sample_spectrahedra

Acknowledgements

ET is partially supported by ANR JCJC GALOP (ANR-17-CE40-0009), the PGMO grant ALMA, and the PHC GRAPE.

Algorithm 5: Billiard Walk (W-Billard)

Input : An LMI F(x) 0 for a spectrahedron S , a point p i , the diameter τ of S and a bound ρ on the number of reflections. Require: p i ∈ S Output : The point p i+1 of the (i + 1)-th step of the walk.

1 ← -τ ln η ; η ← R U((0, 1)); // choose length

Proof. The per-step complexity of W-Billard is dominated by the reflection, which requires O(m 2.697 + m lg L + nm 2 ) arithmetic operations (Lemma 2.7), when we want to approximate the intersection point up to a factor of 2 -L . Since we allow at most ρ reflections per step, the total complexity becomes O(ρ(m 2.697 + m lg L + nm 2 )).

Hamiltonian Monte Carlo with Reflections

Algorithm 6: HMC w reflection (W-HMC-r)

Input : An LMI F(x) 0 representing a spectrahedron S , a point p i , the diameter τ of S and a bound ρ to the number of reflections. Require: p i ∈ S Output : The point p i+1 of the (i + 1)-th step of the walk.

Compute trajectory Φ(t) from ODE (10);

Hamiltonian Monte Carlo (HMC), can be used to sample from any probability distribution π. Our focus lies again on the log-concave distributions, that is π(x) ∝ e -α f (x) . We exploit the setting in [START_REF] Lee | Algorithmic theory of odes and sampling from well-conditioned logconcave densities[END_REF] as they approximate the Hamiltonian trajectory with a polynomial curve. In this setting, if we assume that f is a strongly convex function, then the mixing time of HMC is O(k 1.5 log(n/ )), The average #iteration / runtime / failures over 10 generated S -n-m, to achieve relative error ≤ 0.05. The walk length is one for W-HMC-r and W 1 = 4 √ n and W 2 = 4n for W-HnR. With "failures" we count the number of times the method fails to converge. Also m is the dimension of the matrix in LMI and n is the dimension that S -n-m lies.

Appendix A. Additional proofs

To prove lemma 2.5 we will need the following lemmas.

Lemma Appendix A.1 (Partial Derivative of Determinant). Let A be a symmetric m × m matrix. Then

where c i j the cofactor of A i j .

Proof. From Laplace expansion:

Notice that c 1 j , • • • , c m j are independent of A i j , so we have

Proof. The function det F(x) is the composition of det A and A = F(x), so from Lemma Appendix A.1 and the chain rule:

where A k i j the i j-th element of matrix A k


Zuoying Wen

Xiaofeng Tang
email: tangxf@aiofm.ac.cn

Christa Fittschen

Cuihong Zhang

Tao Wang

Chengcheng Wang

Xuejun Gu

Weijun Zhang
email: wjzhang@aiofm.ac.cn

Online analysis of gas-phase radical reactions using vacuum ultraviolet lamp photoionization and time-of-flight mass spectrometry

A home-made vacuum ultraviolet (VUV) photoionization time-of-flight mass spectrometer has been developed and coupled to an atmospheric simulation chamber operated at atmospheric pressure and to a fast flow tube at low pressure (1-10 Torr). Gas sampling from the chamber is realized directly via a capillary effusive beam, and sampling from the flow tube is via a continuous molecular beam inlet. Both devices are connected simultaneously to the ionization chamber and can be switched in-between within minutes to study gas-phase radical reactions of atmospheric interest in a large range of reaction conditions and reaction times (from milliseconds in the flow tube to hours in the simulation chamber). A cage-shaped photoionization source combined with a commercial 10.6 eV krypton lamp has been developed to provide a high ion collection efficiency along the long light path in the cage. This way, a multiplexed detection with high sensitivity down to the sub-ppbv concentration range, e.g. a limit of detection (LOD) of 0.3 ppbv for benzene and 1.3 ppbv for the methyl radical, is obtained. The performance and suitability of the setup is illustrated by the study of the chlorine-initiated oxidation reaction of toluene in the atmospheric simulation chamber and in the fast flow tube. Stable products and reactive intermediates have been well determined and their reaction dynamics are discussed.

I. INTRODUCTION

Mass spectrometry as a versatile and sensitive analysis tool has the power to provide qualitative and quantitative information and has already been successfully applied in a wide range of fields. In comparison to hard ionization techniques, soft ionization techniques such as matrix-assisted laser desorption ionization (MALDI), electrospray ionization (ESI), chemical ionization (CI) and photoionization can produce high molecular ion yields with minimal fragmentation and have had a large influence on the development of mass spectrometry. [START_REF] Gross | Mass spectrometry[END_REF][START_REF] Jirasko | [END_REF] Among them, single-photon ionization with vacuum ultraviolet (VUV) radiations has the ability to ionize compounds including both polar and nonpolar species without matrix effects and has the potential to be a universal soft ionization method. The coupling of VUV photoionization with mass spectrometry has attracted a great deal of attention and provided valuable data at an unprecedented level of detail. [3][4][5] VUV photoionization mass spectrometry has made considerable progresses benefiting from the improvement of light sources including VUV lamps, lasers and synchrotron radiation. Synchrotron radiation has the merit of high photon flux, wide energy range and easy tunability, and is considered as an ideal single-photon ionization source. The description of VUV photoionization endstations at synchrotron facilities and their diverse applications can be found in the literature. [6][7][8][9] But the problem of synchrotron facilities is that it cannot be moved outside for field experiments and that the available beamtime for each user is very limited. For VUV lasers, frequency tripling of the third harmonic emitted from a Nd:YAG laser at 355 nm and four-wave mixing of UV pulses in a rare gas cell are often utilized to generate VUV photons. [10][11][12] But lasers are difficult to tune in the VUV range and the setups are sensitive to the fluctuation of laser pulse energy. 3 The combination of VUV lamp photoionization and mass spectrometer for use in kinetics was first developed by Bayes et al., and further developed by Gutman and Fockenberg et al. [13][14][15] But, unlike synchrotron radiation and lasers, VUV lamps usually have a low photon flux with a narrow spectral range. However, besides the low price, VUV lamps have the advantage of being robust with a small volume and of being almost maintenance free, well suited for mass spectrometer with a compact design and mobile experiments. In particular with an appropriate design of the ionization source, the disadvantage of VUV lamps, e.g. the low photon flux, can be compensated and then a detection with a high sensitivity has been approached. The limits of detection (LODs) of some typical VUV-lamp photoionization time-of-flight (TOF) mass spectrometers with different lamps and pressures inside the ionization source are listed in Table 1, and are found to be even better than those using lasers or synchrotron radiations. 7,8,16,17 The configuration of the ionization source is important for the detection sensitivity and the recent improvements were mainly benefiting from the pressure increase inside the ionization source. For example, after increasing the ionization source pressure up to 700 Pa, Wang et al obtained a LOD of 15 pptv (parts per trillion by volume) for aliphatic and aromatic hydrocarbons. 19 Employing a RF-powered VUV lamp with a photon flux ~2 × 10 14 photons s -1 and increasing the pressure inside the ionization source up to 1300 Pa, Shu et al. achieved LODs ~ 0.30-0.69 pptv for volatile aldehydes. 21 But the ionization mechanism of Shu's setup is complex: photoionization induces other ionization mechanisms, leading to adduct formation such as acetone and acetic acid. VUV-lamp photoionization mass spectrometers have been mainly applied on the detection of stable species like volatile organic compounds (VOCs), and seldom on reactive intermediates. Also up to now we only find very few reports on the LOD of reactive intermediates like radicals with VUV photoionization mass spectrometry in the literature. 7,17 One possible reason is that the available photoionization cross sections of radicals are very scarce [22][23][24][25] and so it is difficult to quantify them with photoionization mass spectrometry, especially for VUV lamps often with multiple output photon energies.

Another reason may be that within the above mentioned high pressure environment the concentration of reactive radicals rapidly declines during transfer into the ionization source.

In this manuscript, we present a VUV-lamp photoionization orthogonal acceleration reflectron time-of-flight mass spectrometer (VUVPI-TOFMS) for online analysis of radical reactions of atmospheric interest. A novel cage-shaped photoionization source combined with a commercial 10.6 eV krypton lamp has been developed. A focusing electric field 26 is then formed to allow all the ions produced along the long light path of the krypton lamp to be extracted, making the achieved LODs of the VUVPI-TOFMS comparable to the above setups with medium pressure ionization source in Table 1. Here, the pressure inside the cage-shaped photoionization source is low at ~10 -2 Pa, and both stable species and reactive intermediates can thus be transferred into the photoionization region to be ionized and analyzed.

In addition, two reactors, an atmospheric simulation chamber operated at atmospheric pressure and a fast flow tube at low pressure (1-10 Torr), are coupled simultaneously to the ionization chamber of the VUVPI-TOFMS in order to initiate radical reactions at different pressures and to study the time-resolved evolution of the gas mixture over different reaction times, in the minutes to hours range in the atmospheric simulation chamber and in the millisecond range in the fast flow tube. As far as we known, this is the first setup which allows the simultaneous connection of both, an atmospheric simulation chamber and a fast flow tube, independently to the same mass spectrometer in order to investigate atmosphere radical reactions at different pressures and reaction times, switching from one reactor to the other is possible within minutes. A multiplexed detection with a high sensitivity down to sub-ppbv concentration range has been obtained. As representative examples, the chlorine-initiated oxidation reaction of toluene has been studied in the atmospheric simulation chamber at atmospheric pressure and in the fast flow tube at 2 Torr, and preliminary results are also discussed.

II. EXPERIMENTAL

A schematic diagram of the setup is shown in

A. Chemical reactors

The configuration of the atmospheric simulation chamber has been introduced in our previous publications and is mainly made of a 150 L Teflon bag surrounded by several UV lamps emitting photons with wavelength at around 254 nm. 27,28 Clean air from a pure air generator (AADCO 737-series) was used as bath gas in the Teflon bag. The Teflon bag was operated at atmospheric pressure and presently was connected to the VUVPI-TOFMS directly via a fused silica capillary (Part 160-1010, Agilent Technologies, US) with a length of 20 cm and 0.1 mm inner diameter. As shown in Fig. 1, the capillary acted as an atmospheric pressure inlet and its tip penetrated directly into the photoionization region. The pressure of the ionization chamber increased from 10 -5 to 10 -2 Pa with the capillary sampling on and then an effusive beam was formed. Stable species could be transferred from the Teflon bag to the photoionization source of the VUVPI-TOFMS through the capillary.

The fast flow tube is similar to the previous design 29,30 and was composed of a 45 cm long Pyrex main tube with 20/16 mm outer/inner diameter, and a 50 cm long coaxial movable injector with 6/4 mm outer/inner diameter. Fluorine or chlorine atoms utilized as radical initiators were generated from diluted F 2 or Cl The atmospheric simulation chamber and the fast flow tube were operated independently and the sampling from one or the other reactor could be exchanged in minutes, just manually opening or closing the valves on their individual gas circuits shown in Fig. 1. The distinctive characteristics of the atmospheric simulation chamber is that it can be used for studying radical reactions at atmospheric pressure by following the concentration of stable species over reaction times from minutes to several hours, in contrary to the fast flow tube where the reaction time is much shorter (several 10 msec) but reactive intermediates like radicals can be detected through the molecular beam sampling.

B. Cage-shaped photoionization source

The species sampling from the capillary effusive beam or the continuous molecular beam cross the VUV light at a right angle in the photoionization region and species having an ionization energy below the VUV photon energy can be ionized. A commercial krypton lamp (PKS106, Heraeus, Germany) with two atomic resonant lines at hν = 10.0 and 10.6 eV was employed as the VUV light source. The krypton lamp was sealed with a MgF 2 window and its light entered into the photoionization region directly without focusing. The VUV light beam of the krypton lamp has a larger size than those of lasers and synchrotron radiations and its photon flux is ~1×10 11 photons s -1 . Presently to utilize the photon flux as complete as possible and thus to promote the photoionization efficiency, the krypton lamp was installed directly inside the ionization chamber, right above the photoionization region within the vacuum.

In addition, in contrary to our previous design, 30 a cage-shaped photoionization source with a cage electrode and a plate electrode has been developed to enhance the detection sensitivity, as shown in Fig. 2(a). The cage electrode has an inner diameter of 18 mm with four well-distributed holes (6 mm diameter) on the circle and one hole on the head to let the VUV light and gas in and out. The working pressure inside the photoionization region was low at about 10 -2 Pa and the molecular mean free path is calculated to be ~100 cm, much larger than the dimensions of the photoionization region, allowing both stable species and reactive radicals to be transferred and to survive in the photoionization region. The plate electrode with a 6 mm inner diameter was located coaxially on the right of the cage electrode and was insulated with polyether-ether-ketone (PEEK) columns.

C. Orthogonal acceleration reflectron TOF mass analyzer

A home-made orthogonal acceleration reflectron TOF mass analyzer 30 was employed to analyze ion masses and its detailed configuration of electrodes is presented in Fig. 3.

At the downstream of the cage-shaped photoionization source, an Einzel lens with three 2.

III. RESULTS AND DISCUSSION

A. Focusing effect of the cage-shaped photoionization source

In this section, we focus on the focusing effect of the cage-shaped photoionization source and its performance. The electric field and the ion trajectories inside the cageshaped photoionization source were simulated and are presented as blue dotted lines and red solid lines in Fig. 2(a). The voltages of the cage electrode and the plate electrode were optimized and set at 15 and 0 V, respectively. This way, a focusing electric field 26 was formed and all the ions produced along the long light path of the VUV lamp could be extracted and passed through the hole of the plate electrode to be analyzed and detected.

To compare the ion focusing effect of the cage-shaped photoionization source, the electric field and the ion trajectories inside the traditional ionization source with two parallel plate electrodes 8,30 were also simulated and are shown in Fig. 2(b). In this case near parallel electric field was formed between the two plate electrodes and only the small fraction of the ions within the large volume, located at the center of the photoionization region, was extracted through the second plate electrode to be detected.

The ion collection efficiencies of the focusing electric field (Fig. 2a, the cage-shaped photoionization source) and the parallel electric field (Fig. 2b) are calculated and listed in Table 3. In the ionization chamber the ions can be produced along the long light path of the krypton lamp and so the photoionization region should take a shape of a cigar. In the calculation the length of the photoionization region is fixed at 18 mm, i.e. the inner diameter of the cage electrode, whereas its diameter is varied from 1 to 3 mm, as the exact diameter of the VUV light is not known. The calculated results show that the ion collection efficiency is slightly dependent on the volume of the ionization region, whereas the ratio of the ion collection efficiencies of the focusing electric field (E cage ) to the parallel electric field (E para ) is almost independent. The ion collection efficiency of the focusing electric field is about three times higher than that of the parallel electric field, corresponding to the ratio of the inner diameter of the cage electrode to the one of the plate electrode. In addition, to experimentally demonstrate the focusing effect of the cage-shaped photoionization source, we have also performed VUV lamp photoionization experiments with the above two kinds of photoionization sources, the cage-shaped photoionization source and the parallel photoionization source, by using the capillary effusive beam inlet and the VUVPI-TOFMS. The photoionization mass spectra of a gas mixture of benzene, toluene and p-xylene were measured with the two photoionization sources respectively, whereas the other experimental conditions were the same, and are presented in Fig. 4. In the mass spectra, the intense peaks at m/z = 78, 92 and 106 are assigned as the ions of benzene, toluene and p-xylene, and the accompanied m/z = 79, 93 and 107 peaks can be assigned as their individual 13 C isotopic species. The ionization energies of benzene, toluene and p-xylene locate at 9.244, 8.828 and 8.445 eV, respectively, close to the photon energies of the krypton lamp. 31 Benzene, toluene and p-xylene had been softly ionized with the VUV photons and no fragment ions were observed in the mass spectra.

From the width of the ion peak of benzene, the mass resolution of the VUVPI-TOFMS is measured to be M/M ~ 1500 (FWHM, the full width at half maximum). In addition, the ion intensities of benzene, toluene and p-xylene in the mass spectrum obtained with the cage-shaped photoionization source have been increased for 2.6, 2.3 and 2.1 times, respectively, compared to those with the parallel photoionization source, agreeing well with the above simulated results.

B. LODs of stable species

A gas mixture of benzene, toluene and p-xylene with calibrated concentrations of 10 ppmv (Nanjing Special Gases Co. Ltd., China) was used to characterize the detection sensitivity of the VUVPI-TOFMS. The calibrated gas was dynamically diluted with helium through mass flow controllers and their photoionization mass spectra were measured by sampling through the capillary effusive beam. As shown in Fig. 5 The LOD of the VUVPI-TOFMS can be calculated with a signal-to-noise ratio of 3 from the photoionization mass spectrum via the equation E1, [START_REF] Williams | Twenty-Fourth Symposium (International) on Combustion[END_REF][START_REF] Heger | [END_REF] LOD = 3σc / (S -B) (E1)

where σ is the amplitude variance of the signal between peaks, c is the species concentration, S is the peak height and B is the mean value of the noise. In the mass spectrum of Fig. 1, the LOD of this VUVPI-TOFMS is lower than those of the setups using synchrotron radiations, lasers and some lamps as photoionization light source. In addition, although the pressure inside the cage-shaped photoionization source is low at only 10 -2 Pa, this VUVPI-TOFMS benefits from the ion focusing electric field which makes its LOD comparable to the setups with medium pressure ionization sources. 18

C. LOD of radicals

In comparison to stable species, radicals must be produced from reactions online and normally their concentrations are not known in advance. In photoionization mass spectra, the ion signal I is proportional to the specie's concentration N and its photoionization cross section σ, following the equation E2,

I = C σ N m β(P,T) (E2)
where C is the instrument factor, m is the ion mass and β is the mass discrimination factor. 23 But because the VUV lamp has multiple output photon energies and thus the photoionization cross sections are not known, it is difficult to quantify the absolute concentrations of radicals by using the VUV lamp photoionization with the equation E2.

As we suggested before, the radical concentration can be determined by measuring their time resolved concentration during the self-reaction in the fast flow tube, and a LOD of ~3 ppbv has been reported for the methyl radical (CH 3 ). 30 Very recently, to demonstrate the availability of the self-reaction method to measure the radical concentrations, we carried out synchrotron experiments at Hefei Light Source with the fast flow tube and the VUV photoionization mass spectrometer. Ethyl radicals (C 2 H 5 )

were produced through the reaction of fluorine atoms with ethane. The concentrations of ethyl radicals were measured at several fixed synchrotron photon energies with two different methods, the self-reaction method with the available self-reaction rate constant of C 2 H 5 and the photoionization cross section method with the calibrated concentrations of NO gas, together with the available photoionization cross sections of C 2 H 5 and NO.

The concentrations of C 2 H 5 measured with the two different methods agree well and their relative errors are within 10-15%. This work will be the subject of a forthcoming paper.

Presently the methyl radical self-reaction has been chosen again to measure the LOD of the new VUVPI-TOFMS with the cage-shaped photoionization source. The methyl radicals were produced in the fast flow tube through the H-abstraction reaction from methane by fluorine atoms, which were formed in the microwave discharge of F 2 diluted in helium. A typical VUV lamp photoionization mass spectrum was measured at a pressure of 2 Torr inside the fast flow tube and is presented in Fig. 6(a). The peak at m/z = 15 is assigned as CH 3 with its adiabatic ionization energy at 9.80 eV. 31,34 Without other reactant gases, CH 3 mainly decreased through self-reaction in the flow tube which has already been investigated in detail before. 35,36 The self-reaction of CH 3 mainly produces the product C 2 H 6 with its ionization energy located at 11.52 eV, above the photon energies of the VUV lamp, and so the molecule was not observed in the mass spectrum. 31 Briefly, the self-reaction of CH 3 is a second-order reaction and its inverse concentration (1/[CH 3 ]) changes linearly with reaction time, t, following the equation E4,

1 [CH 3 ] = 1 [CH 3 ] 0 + 2k𝑡 (E4)
where [CH 3 ] 0 is the initial concentration of CH 3 and k is the rate constant at 5.52×10 -11 cm 3 molecule -1 s -1 . 35,36 As shown in Fig. 6(b), the inverse intensities of CH 3 in the photoionization mass spectra with different initial concentrations of CH 3 change linearly with reaction time and agree well with the kinetics of a second-order reaction, also demonstrating that methyl radicals majorly performed self-reaction in the flow tube.

As shown in the equation E2, the ion signals of CH 3 in the photoionization mass spectra are proportional to their concentrations, I CH3 ∝ f[CH 3 ], and the scaling factor, f, can be calibrated with the known rate constant 35,36 and the slopes of the lines in Fig. 6(b) through the equation E4'.

1

𝑓[CH 3 ] = 1 𝑓[CH 3 ] 0 + 2k 𝑓 𝑡 (E4')
The slopes of the lines (2k/f) in Fig. 6(b) are fitted and take an average value of 0.007 counts -1 s -1 , and then the scaling factor is determined to be f = 1.6×10 -8 counts cm 3 molecule -1 with the self-reaction rate constant k at 5.52×10 -11 cm 3 molecule -1 s -1 . 35,36 .

So the absolute concentration of CH 3 at each reaction time with different initial concentrations are determined. For example, the concentration of CH 3 in Fig. 6(a) is calculated to be 1.6×10 12 molecule cm -3 from the ratio of its intensity in the mass spectrum and the scaling factor. Then based on the equation E1, the LODs of CH 3 are calculated from the photoionization mass spectra with different conditions and shown in Fig. 6(c). The LOD of CH 3 takes an average value of 1.3 ± 0.2 ppbv, which is 2.3 times lower than our previous result of 3 ppbv, 30 benefiting of the ion focusing electric field and in good agreement with the improvement presented in the Section III. A.

D. Initial applications of VUVPI-TOFMS: The chlorine-initiated oxidation of toluene

As representative examples, the chlorine initiated oxidation of toluene has been investigated in the atmospheric simulation chamber at atmospheric pressure and in the fast flow tube at 2 Torr.

In the simulation chamber, Cl-atoms were generated through the UV photolysis of ~5 ppmv (parts per million by volume) trichloroacetyl chloride (CCl 3 CClO). 37 Clean air was added as bath gas and the initial concentration of toluene was 1.5 ppmv.

The photoionization mass spectra with the UV lamps off and on were measured and are presented in Fig. 7 time resolution of 1 min. We can see that after turning on the UV lamps (at ~4 min), the ion intensity of toluene (in black, left y-axis) decreased and the ion intensity of the products increased within several minutes. For benzaldehyde (in cyan, m/z = 106, right yaxis), the signal increased to the maximum at ~11 minutes and then decreased. This can be explained by the fact that (i) the production of benzaldehyde is exclusively due to the Cl + toluene reaction and thus, as toluene decreases, its production decreases as well; and

(ii) the secondary reaction of Cl + benzaldehyde is faster (k = 8.9×10 -11 cm 3 molecule -1 s - 1 ) 41 than the Cl + toluene reaction (k = 5.8×10 -11 cm 3 molecule -1 s -1 ) 42 , and thus with increasing benzaldehyde and decreasing toluene concentrations the Cl-atoms react increasingly with benzaldehyde rather than with toluene. As the photoionization cross section of benzaldehyde is almost the same as that of toluene at hν = 10.6 eV, 43 the benzaldehyde concentration at ~11 minutes was calculated to be about 50 ppbv based on the calibrated concentration of toluene in the photoionization mass spectra.

For the products glyoxal (m/z = 58), methyl glyoxal (m/z = 72) and benzene (m/z = 78), their intensities increased with the reaction time and approached a steady-state after ~15 minutes, meaning that the production rate of these species is roughly as fast as their consumption. This is a very brief analysis and the inherent mechanism in the atmospheric simulation chamber needs to be simulated to get a more detailed answer, 44,45 which is out of the scope of this manuscript.

The chlorine-initiated oxidation of toluene was also studied in the fast flow tube at 2

Torr. Chlorine atoms were generated by the microwave discharge of 1% Cl 2 diluted in helium and fed into the fast flow tube via the injector. Toluene was added by bubbling helium through a reservoir with liquid toluene, and oxygen gas was added through a calibrated mass flow meter, both were introduced into the main tube through side arms. obtained from the H-abstraction reaction of toluene with chlorine atom, was also detected and can be clearly identified at m/z = 91 in the photoionization mass spectra. In addition, two other peaks at m/z = 126 and 128 can be observed in the mass spectra too and are assigned as C 7 H 7 35 Cl and C 7 H 7 37 Cl, the products of the secondary reaction of the benzyl radical with chlorine atoms, or with remaining Cl 2 . 46,47 The ratio of the C 7 H 7 35 Cl and C 7 H 7 37 Cl peaks is about 3:1, agreeing well with the natural abundances of 35 and 37 Cl isotopes.

After addition of excess oxygen relative to benzyl radicals and chlorine atoms into the flow tube, the peak at m/z = 106 appeared in the photoionization mass spectrum of cm 3 molecule -1 s -1 ), 47 but higher than the literature data of Bartels et al (4.98×10 -13 cm 3 molecule -1 s -1 ) 42 . In our case, neglecting the loss of benzyl radicals through selfreaction, wall and reaction with remaining Cl 2 will lead to somewhat overestimation of the rate constant and can explain the difference.

IV. CONCLUSIONS

Fig. 1 .FIG. 1 .

 11 FIG. 1. Schematic diagram of the VUV-lamp photoionization orthogonal acceleration reflectron time-of-flight mass spectrometer with a capillary effusive beam and a continuous molecular beam connected to an atmospheric simulation chamber and a microwave discharge (MW) flow tube to investigate radical reactions at different pressures and reaction times.

  photoionization region. The distance between the tip of the injector and the skimmer can be adjusted and thus the reaction time can be changed.

FIG. 2 .

 2 FIG. 2. (a) Geometry of the cage-shaped photoionization source and simulated ion trajectories (red solid lines) along the long light path of VUV lamp. The equi-potential curves (blue dotted lines) show the ion focusing effect of electric field in the cage. (b) The ion trajectories and the equi-potential curves of a parallel repelling ionization source.

  FIG. 3. (a)The detailed configuration of the VUVPI-TOFMS, mainly consisting of the cageshaped photoionization source, an Einzel lenses and an orthogonal acceleration reflectron TOF mass analyzer with ion trajectories in red. (b) The 3D layout of the setup with the VUV lamp, the capillary effusive beam inlet and the fast flow tube. (c) The dual pulses and ion signal process of the VUVPI-TOFMS.After being accelerated with 2400 eV kinetic energy within the acceleration region of the reflectron TOF mass analyzer, the ions flied into the field free drift region, entered the reflectron region with a double space-focusing configuration and then backed again into the field free drift region. Finally, the ions were detected by a pair of micro-channel plates (MCPs, 50 mm diameter) with a chevron configuration. The MCPs signals were amplified by a preamplifier and recorded by a multiscaler counter (TA2000B-1 & P7888-2, FAST ComTec, Germany). As shown in Fig.3(b), the VUVPI-TOFMS has a compact design with a length/width of 0.6/0.4 m and its typical parameters are listed in Table2.

FIG. 4 .

 4 FIG. 4. VUV lamp photoionization mass spectra of benzene, toluene and p-xylene, measured with the cage-shaped photoionization source (in black) and the parallel photoionization source (in red) of Fig. 2.

  (a), the ion intensities of benzene (m/z = 78), toluene (m/z = 92) and p-xylene (m/z = 106) integrated from their individual peak areas in the mass spectra decreased linearly with their diluted concentrations. As an example, a typical mass spectrum of the species with diluted concentrations of 22 ppbv, accumulated for 60 s, is presented in Fig. 5(b) and displays a near zero baseline. The ion peaks of benzene, toluene and p-xylene with their individual 13 C isotopic species at m/z = 79, 93 and 107 are observed in the mass spectrum.

FIG. 5 .

 5 FIG. 5. (a) The ion intensities of benzene (m/z = 78), toluene (m/z = 92), and p-xylene (m/z = 106) in photoionization mass spectra change linearly with their concentrations. (b) Photoionization mass spectrum of the mixture at a concentration of 22 ppbv acquired via the capillary effusive beam sampling.

  5(b), σ takes the value of 0.18 between the peaks of m/z = 80 and 90, and B is calculated to be 0.03. The ion peaks of benzene, toluene and p-xylene have the counts of 51, 35 and 24. Then the LODs of the VUVPI-TOFMS with the capillary sampling are calculated at 0.3, 0.4 and 0.5 ppbv for benzene, toluene and p-xylene, respectively, with an accumulation time of 60 s. As listed in Table

CH 3 +

 3 FIG. 6. (a) Photoionization mass spectrum of CH 3 . (b) The inverse intensity of CH 3 in photoionization mass spectra changes linearly with reaction time in the self-reaction with different initial concentrations. (c) The LODs of CH 3 determined from photoionization mass spectra (open symbols) and its average marked with a purple dashed line.

FIG. 7 .

 7 Fig. 7(b). These peaks are assigned to the oxidation products glyoxal (C 2 H 2 O 2 , m/z = 58), methyl glyoxal (C 3 H 4 O 2 , m/z = 72), benzene (C 6 H 6 , m/z = 78) and benzaldehyde (C 7 H 6 O, m/z = 106).38,39 Benzyl alcohol (C 7 H 8 O, m/z = 108) is at the most a very minor product, in contrast to the findings of Nozière et al.40 However, due to the capillary sampling, the reactive intermediates involved in the oxidation process were not transferred into the VUVPI-TOFMS and have thus not been detected.

Fig. 7 (

 7 Fig. 7(c) presents the change of the ion intensities of toluene and of the oxidation products as a function of the reaction time in the atmospheric simulation chamber, with a

FIG. 8 .

 8 FIG. 8. Photoionization mass spectra of the chlorine-initiated oxidation of toluene (a) without and (b) with adding oxygen gas in fast flow tube. (c) Plot of the logarithmic ion intensity of the benzyl radical versus reaction time, with the initial oxygen concentration of [O 2 ] = 8.7 × 10 13 cm -3 .Photoionization mass spectra were obtained without and with adding oxygen in the fast flow tube and are presented in Fig.8(a) and (b). Similar to the above chamber experiments, the most intense ion peak in the photoionization mass spectra is at m/z = 92 and is attributed to toluene, together with its13 C-toluene isotope at m/z = 93. But in this setup with the molecular beam sampling, the reaction intermediate benzyl radical (C 7 H 7 ),

Fig. 8 (

 8 Fig. 8(b) and is assigned as benzaldehyde (C 7 H 6 O) from the oxidation reaction of the benzyl radical, and at the same time the ion intensities of the secondary reaction products C 7 H 7 35 Cl and C 7 H 7 37 Cl decreased. Fig. 8(c) presents the ion intensity of the benzyl radical versus reaction time, obtained through changing the distance between the tip of the injector and the detection region. A good linearity can be found for the decrease of the logarithmic ion intensity of the benzyl radical as a function of the reaction time, in agreement with a pseudo-first-order reaction kinetics. The obtained pseudo-first order rate constant is 142 s -1 , which would result in a rate constant of k C7H7+O2 = 1.6×10 -12 cm 3 molecule -1 s -1 , in accordance with the result of Nelson and McDonald (1.5×10 -12

  VUV photoionization orthogonal acceleration reflectron time-of-flight mass spectrometer (VUVPI-TOFMS) based on a commercial 10.6 eV krypton lamp has been developed to study gas-phase radical reactions. An atmospheric simulation chamber operated at atmospheric pressure and a fast flow tube operated at low pressure (1-10 Torr) are employed as chemical reactors to initiate radical reactions of atmospheric interest at different pressures and allowing studying the time-resolved evolution of the gas mixture over different reaction times, in the minutes to hours range in the atmospheric simulation chamber and in the milliseconds range in the fast flow tube. The atmospheric simulation chamber and the flow tube are connected to the VUVPI-TOFMS for online product 20 analysis via a capillary effusive beam and a continuous molecular beam, respectively, which can be operated independently and switched conveniently in-between within minutes.Several strategies like installing the krypton lamp just above the photoionization region within vacuum and utilizing dual pulses to reduce background signals have been employed in the VUVPI-TOFMS to enhance its detection sensitivity. In particular, a cage-shaped photoionization source has been designed leading to an ion focusing electric field to allow all the ions formed along the long light path of the krypton lamp to be extracted and collected. Thus the novel configuration of the photoionization source could compensate the low photon flux of the krypton lamp and provide a high detection efficiency. Some preliminary experiments have been performed including the selfreaction of methyl radicals and the chlorine-initiated oxidation of toluene in the fast flow tube or the atmospheric simulation chamber. The present experimental results have demonstrated that both stable species and reactive intermediates can be detected with a high sensitivity: the achieved LODs are down to the sub-ppbv concentration range, e.g. 0.3 ppbv for benzene and 1.3 ppbv for the methyl radical, without any pretreatments, providing a great potential to reveal the complex chemical mechanism of atmospherically relevant radical reactions in detail.

Table 1 .

 1 Limits of detection (LODs) of typical time-of-flight mass spectrometers with different photoionization light sources and pressures (P ion ) inside the ionization source.

	Light sources	λ (nm)	P ion. (Pa)	Photon flux (photon s -1 )	LODs a (ppbv) time /Accumulated	Refs.
	ALS	122	~10 -2	5×10 13	20/60s b	7
	NSRL	89	~10 -2	~10 13	<100/5s c	8
	Laser	118	~10 -3	2×10 11	8/5s	16
	Laser	118	~10 -3	---	3/	17
	EBEL	126	~10 -3	1×10 14	27/1s	20
	Kr lamp	117	30	~10 11	3/10s	18
	Kr lamp	117	700	~10 11	0.015/60s	19
	RF Kr lamp	124	1000	5×10 14	0.003/10s	21
	Kr lamp	117	10 -2	~10 11	0.3(1.3)/60s d	This work

a Most of LODs for benzene, b for acetone, c for CO 2 , d 0.3 ppbv for benzene and 1.3 ppbv for the methyl radical. ALS: Advanced Light Source, US; NSRL: National Synchrotron Radiation Laboratory, China; EBEL: Electron beam pumped excimer lamp.

Table 2 .

 2 Typical parameters of the VUV-lamp photoionization orthogonal acceleration time-offlight mass spectrometer.

	Components	Key parameters
	Ionization source	Cage shape with 18 mm diameter; +15 V
	Einzel lens	0 V, -30 V, 0V
	Slit	2 × 6 mm 2 , 0 V
	Extraction region	12 mm length; ±600 V pulse voltage,
		35 kHz, 1 μs width
	Acceleration region	30 mm length; -2400 V
	Electric field-free drift region 320 mm length; -2400V
	Dual-stage ion reflector	-195 V for mirror 1,
		+795 V for mirror 2
	MCPs detector	900 V per piece of MCP

Table 3 .

 3 Calculated ion collection efficiencies of the focusing electric field (E cage ) and the parallel electric field (E para ) in Fig.2.

	Ionization region (mm)	E cage	E para	E cage /E para
	∅1×18	0.84	0.32	2.7
	∅2×18	0.79	0.29	2.7
	∅3×18	0.73	0.26	2.8

ACKNOWLEDGEMENTS

This work was financially supported by the National Natural Science Foundation of China (NSFC, Nos. 21773249, 91961123, 91644109), the National Key R&D Program of China (NO. 2016YFC0200300) and the International Partnership Program of Chinese Academy of Sciences (No. 116134KYSB20170048). C.F. would like to thank the support of the Chinese Academy of Sciences President's International Fellowship Initiative (PIFI, No. 2018VMA0055) and the Xu Guangqi program (No. 43270ZD).