Sarmad Munir
email: sarmad.munir@itk.ntnu.no

Morten Hovd
email: morten.hovd@itk.ntnu.no

Guillaume Sandou
email: guillaume.sandou@supelec.fr

Sorin Olaru
email: sorin.olaru@centralesupelec.fr

Controlled contractive sets for low-complexity constrained control

Explicit constrained control is relatively simple when a controlled contractive set is available. However, the complexity of the explicit controller will depend on the complexity of the controlled contractive set. The ability to design a low complexity controlled contractive set is therefore desirable. Most methods for finding controlled contractive sets either assume the use of a constant linear state feedback, or is based on reachable set computations. In the first case, the assumption of a constant linear state feedback is restrictive (as MPC controllers for linear systems are typically piecewise affine), and in the second case the complexity of the controlled contractive set may be very high.

Initial results on the construction of low complexity controlled contractive sets without assuming linear state feedback were reported at the IFAC World Congress in 2014. The present paper addresses shortcomings in the previous results, including the ability to handle identical modes in series (corresponding to a non-diagonalizable A matrix) and oscillatory modes. The paper suggests a method to find a flexible complexity controlled contractive set.

I. INTRODUCTION Model Predictive Control (MPC) has been used in industry

for decades. There is a wide application of MPC in industries nowadays [START_REF] Ogawa | Practice and Challenges in Chemical Process Control Applications in Japan[END_REF]. From slow speed chemical processes to fast systems in automotive sector, MPC is highly useful [START_REF] Hrovat | The Development of Model Predictive Control in Automotive Industry: A Survey[END_REF]. MPC uses the system model information to obtain control law. It allows operating of a system for long periods without expert intervention [START_REF] Holkar | An Overview of Model Predictive Control[END_REF]. Standard MPC is based on online solutions of optimization problems. Due to its computational complexity, its application is limited to the systems with slow dynamics which are not safety critical [START_REF] Hovd | Low Complexity Constraint Control Using Contractive Sets[END_REF]. Explicit MPC provides an answer to these limitations of standard MPC, by formulating the MPC problem as multi-parametric problem. Instead of solving problem online, it can be solved offline and optimal control law can be given as piecewise affine (PWA) functions of the present state [START_REF] Scibilia | Explicit Model Predictive Control: Solutions Via Computational Geometry[END_REF]. Therefore the online MPC computation is transformed to simple evaluation of a PWA function [START_REF] Bemporad | Pistikopoulos The explicit linear quadratic regulator for constrained systems[END_REF]. This allows implementation of simple hardware with fast sampling rate. However, as the problem size increases, number of the regions of the explicit solution and the memory required for storing the explicit solution increases rapidly. This limits the use of explicit MPC to a system with modest number of states and short prediction horizons. Thus, the complexity reduction in explicit Model Predictive Control becomes next big challenge.

There has been significant research going on to find approximate solutions in order to simplify explicit solutions. A simpler solution by directly approximating the control law is formulated in [START_REF] Jones | Approximate explicit mpc using bilevel optimization[END_REF]. Vertex Control approaches can be utilized for simpler controller designs. In the vertex control approach, an admissible input corresponding to vertices is used for controller design [START_REF] Nguyen | Implicit improved vertex control for uncertain, timevarying linear discrete-time systems with state and control constraints[END_REF]. A common way to do this would be to calculate the feasible region for an MPC formulation with guaranteed stability. However, the calculation of feasible regions for MPC may also be computationally complex [START_REF] Scibilia | On feasible sets for MPC and their approximations[END_REF].

A simpler approach is using explicit constrained control with the help of controlled contractive set. In this case, complexity of the Explicit MPC depends on the contractive set complexity. Therefore, finding a low complexity contractive set is essential for complexity reduction in Explicit MPC. The procedure in [START_REF] Dorea | a,b)-invariant polyhedral sets for linear discrete-time systems[END_REF] converges to the maximal controlled contractive set with a specified contraction factor. However, the complexity of the resulting set may be very high. In [START_REF] Hovd | Low Complexity Constraint Control Using Contractive Sets[END_REF], a non-iterative procedure for obtaining lowcomplexity contractive set is proposed and it is shown how the contractive set can be used for designing a controller for a system with state and input constraints. It also proposed an approach to handle oscillatory modes given sufficient contractive dynamics. However, the approach is not generally applicable to systems with oscillatory modes, and is not applicable to identical modes in series. Furthermore, the contractive set obtained in [START_REF] Hovd | Low Complexity Constraint Control Using Contractive Sets[END_REF] is of fixed complexity, without any means of trading off complexity versus the size of the contractive set. The main purpose of this paper is to describe a flexible approach to obtain a low complexity contractive sets enabling the design of low complexity explicit constrained control, which also handles identical and oscillatory modes in series and allows a trade-off between the size of the controlled contractive set and the complexity of the set.

Section 2 describes how to find a contractive set using an optimization based approach. In section 3, different solutions are merged together in order to find a larger contractive set. In section 4, two methods are introduced to simplify complex contractive sets, first method introduces circumscribed ellipsoid technique to simplify the contractive set while second method directly removes hyperplanes in order to obtain a low complexity contractive set. In section 5, polytopic contractive sets are considered for controller design, as polytopic contractive sets are natural starting points for designing explicit constrained controllers. Section 6 provides examples of these approaches and conclusion is drawn in section 7.

II. OPTIMIZATION BASED CALCULATION FOR CONTROLLED CONTRACTIVE SETS

Consider the discrete-time system

x k+1 = Ax k + Bu k (1)
with state and input constraints given by X = {x|H x x ≤ h x } and U = {u|H u u ≤ h u } respectively.

Def inition 1. A closed set P ∈ X with the origin in its interior is called controlled γ-contractive iff, for a given γ ∈ [0, 1) and for all x ∈ P there exists an u ∈ U such that Ax k + Bu k ∈ γP.

The contractive set P for system 1 can be given as

P = {x ∈ R n |F x k ≤ f }. It is contractive if ∀x k ∈ P ∃ u k such that F (Ax k + Bu k ) ≤ γf (2)
In the following an optimization based approach to find a large contractive set will be described. The approach is based on maximizing the volume of an ellipsoid inscribed in the contractive set. Finding the inscribed ellipsoid of maximum volume will be addressed in the next subsection. For the volume of the inscribed ellipsoid to be bounded, the set P itself needs to be bounded, and additional considerations are necessary to ensure that the set is also controlled contractive with the required contraction factor. These issues will be addressed in subsequent subsections.

A. Maximum Volume Ellipsoid inscribed in a set

Consider the set P described by m linear inequalities as P = {x ∈ R n |F i x k ≤ f i , ∀i = 1...m} and the ellipsoid ε given as ε = {Cy + d| y ≤ 1}, where C = C T > 0.

According to [START_REF] Vandenberghe | Determinant maximization with linear matrix inequality constraints[END_REF], the volume of ellipsoid ε is proportional to det(C). The maximum volume ellipsoid ε inscribed in P can be obtained by solving the optimization problem:

max C,d,F,f log(det(C)) (3)
subject to

C = C T > 0 (4)
(f i -F i d)I CF T i F i C f i -F i d 0∀i = 1, ...m (5)
Alternatively, one may maximize the radius z of a ball B having the center at c inscribed inside the polytope P.

B = {x ∈ R n c -x 2 ≤ z} (6)
Where z and c are known as Chebyshev radius and Chebyshev center respectively.

To make sure that the ball is entirely inside P, we ensure that the distance from center c of circle to any hyperplane is greater than or equal to Chebyshev radius z for each constraint {F i x k ≤ f i , ∀i = 1, ...m}. Therefore it imposes following constraints:

F i c -f i ≥ F i 2 z (7)
As c satisfies inequalities F i c ≤ f i , therefore we have

F i c + F i z ≤ f i (8)
Thus, by maximizing Chebyshev radius z, the largest possible ball within a contractive set can be found.

Using ball B instead of ellipsoid ε is simpler but contractive set obtained by inscribing the maximum volume ellipsoid ε, can be much larger than the one obtained by inscribing largest ball B. From here on, we shall consider a maximum volume ellipsoid inscribed in P.

B. Boundedness Constraints

We find maximum volume ellipsoid ε inscribed in P from [START_REF] Holkar | An Overview of Model Predictive Control[END_REF]. But this does not guarantee that contractive set obtained by such formulation will be bounded as well. If the volume of the inscribed ellipsoid is finite (the maxdet problem is well posed), the set P is also bounded. However, we want to impose conditions to make sure that this is the case. To ensure boundedness for contractive set we introduce bounding constraints on set P so that when P is bounded by a set H, the inequalities defining H are redundant when inequalities defining P are added. H should contain the state constraints and added lax constraints, if the state constraints set X is not bounded. Let H be defined as:

H = {x ∈ R n |Hx k ≤ h} (9)
Then from [START_REF] Bitsoris | On the Positive Invaraince of polyhedral sets for discrete -time systems[END_REF] we know that P is contained in H (and hence is bounded), if and only if there exists a matrix M with non-negative elements such that

M F = H (10)
and

M f ≤ h (11)
We can find such M which bounds set P by solving above mentioned problem.

C. Contractive Constraints

The set P is contractive if it fulfills Definition 1. We want to find P = {x ∈ R n |F x k ≤ f } with maximum guaranteed contraction factor such that eq. ( 2) is satisfied. Maximum contraction factor for a given F and f can be found by solving following bi-level optimization problem:

max x γ * (12a) subject to Hx k ≤ h (12b) F x k ≤ f (12c) γ * = min u,γ γ (12d)
subject to

F (Ax k + Bu k ) ≤ γf (12e) H u u ≤ h u (12f)
Here the solution of the lower level problem imposes constraints on upper level problem. By replacing the lower-level problem by the corresponding KKT conditions [START_REF] Hovd | Multi-level Programming for Designing Penalty Functions for MPC Controllers[END_REF], one obtains:

max x,u,γ γ * (13a) subject to Hx k ≤ h (13b) F x k ≤ f (13c) F (Ax k + Bu k ) -γf ≤ 0 (13d) H u u ≤ h u (13e) λ a (F (Ax k + Bu k ) -γf ) ≤ 0 (13f) λ b (H u u -h u ) ≤ 0 (13g) ∇ u,γ L(u, γ, λ) ≤ 0 (13h)
where the Lagrangian function L is given by

L = γ * + λ T a (F Ax + F Bu -γh x ) + λ T b (H u u -h u ) (14)
Here λ a and λ b are Lagrangian multipliers for the inequality constraints. Equations (13f) and (13g) induces nonlinearity in the system. We use binary variables sε{0, 1} to remove this non-linearity as explained in [START_REF] Fortuny-Amat | McCarl A representation and economic interpretation of a two-level programming problem[END_REF] to obtain following single level optimization problem:

max x,u,γ,λ,s γ * (15a) subject to Hx k ≤ h (15b) F x k ≤ f (15c) λ a ≥ 0 (15d) λ b ≥ 0 (15e) λ a ≤ M λ a s (15f) λ b ≤ M λ b s (15g) F (Ax k + Bu k ) -γf ≤ 0 (15h) λ(F (Ax k + Bu k ) -γf ) ≥ -M u a (1 -s) (15i) H u u ≤ h u (15j) H u u -h u ≥ -M u b (1 -s) (15k) ∇ u,γ L(u, γ) = 0 (15l)
Here M λ a ,M λ b ,M u a and M u b are diagonal matrices with appropriate dimensions and sufficiently large elements on the main diagonal. If a candidate contractive set is given, its maximum contraction factor can be found using the formulation described above.

D. Problem Formulation

Previous subsections imposed the constraints on P. We want to maximize volume of an ellipsoid ε so that it is inscribed in contractive set P such that boundedness and contractive constraints are satisfied. Thus we have a following optimization problem to obtain P:

max C,d,F,f log(det(C)) (16a) subject to C = C T > 0 (16b) (f i -F i d)I CF T i F i C f i -F i d 0∀i = 1, ...m (16c)
M F = H (16d) M f ≤ h (16e) max x,u,γ,λ,s γ * (16f) subject to Hx k ≤ h (16g) F x k ≤ f (16h) λ a ≥ 0 (16i) λ b ≥ 0 (16j) λ a ≤ M λ a s (16k)
λ b ≤ M λ b s (16l)
F (Ax k + Bu k ) -γf ≤ 0 (16m) λ(F (Ax k + Bu k ) -γf ) ≥ -M u a (1 -s) (16n)
H u u ≤ h u (16o) H u u -h u ≥ -M u b (1 -s) (16p) ∇ u,γ L(u, γ) = 0 (16q)
The problem formulated above can be used to find the maximum volume ellipsoid ε and the corresponding contractive set P = {x ∈ R n |F x k ≤ f } inside which the ellipsoid ε resides. The advantage of this method is that the complexity of contractive set is flexible. There is a tradeoff between complexity of the contractive set reduced and volume of the set. The complexity of the contractive set can be selected according to the specific application. Note that problem formulated above can be simplified by removing parameter f by normalizing set P such that P = {x ∈ R n |F x k ≤ 1}. If the problem formulation is symmetric, this can be exploited to further simplify the computations. Example VI-A shows the effectiveness of this algorithm by comparing the results with the method described in [START_REF] Dorea | a,b)-invariant polyhedral sets for linear discrete-time systems[END_REF]. The contractive sets obtained of different complexities for the system with identical modes are shown in example VI-B.

E. Solving the Formulated Problem

For the optimization problem in eq [START_REF] Geyer | Optimal complexity reduction of polyhedral piecewise affine systems[END_REF], the possible simplifications mentioned in the previous subsection notwithstanding, is a large optimization problem with highly nonconvex constraints. Solving this problem to a (provable) global optimality is therefore very difficult and computationally costly. Luckily, for this problem a solution need not be globally optimal in order to be useful. In this work, Particle Swarm Optimization (PSO) is used to find good solutions to the problem in eq [START_REF] Geyer | Optimal complexity reduction of polyhedral piecewise affine systems[END_REF].

PSO technique described in [START_REF] Van Den Bergh | E ngelbrecht A New Locally Convergent Particle Swarm Optimiser[END_REF] is used to solve problem formulated in previous subsection. It is a population based stochastic method. Particles, candidate solution, are moved around search space with certain position and velocity. Each particle is moved towards best position in search space, so we move the particle swarm towards best solutions. Particle Swarm Optimization is an unconstrained optimization. In order to apply PSO, eq [START_REF] Geyer | Optimal complexity reduction of polyhedral piecewise affine systems[END_REF] is converted into an unconstrained problem by appending suitable penalty functions accounting for constraint violations to the optimization criterion.

While exploring the search space, the PSO typically finds multiple feasible solutions. Such solutions can be merged to obtain an enlarged contractive set, as explained in the next section.

III. MERGING OF CONTRACTIVE SETS

Def inition 2: The convex hull of a family of sets S is the smallest convex set containing the entire family of sets S.

The volume of the contractive set can be increased by finding the convex hull of different solutions obtained by the Particle Swarm Optimization. Unfortunately, it will not only increase the volume of the contractive set but will also increase the complexity of set.

T heorem 1: The convex hull of the contractive sets S 1 and S 2 with contraction factors γ 1 and γ 2 respectively, is γ contractive, where γ = max(γ 1 , γ 2 ).

P roof : Let x 1 (k) and x 2 (k) be points on the contractive sets P 1 and P 2 respectively whose convex hull is denoted by P 0 . Assume that contraction factors for P 1 and P 2 are γ 1 and γ 2 respectively. Then for

x 1 (k)∃u 1 (k) ∈ U |x 1 (k+1) ∈ γ 1 P 1 and x 2 (k)∃u 2 ∈ U |x 2 (k + 1) ∈ γ 2 P 2 .
Let x 0 (k) be point obtained outside P 1 and P 2 but inside the convex hull. Then, x 0 (k) can be expressed as:

x 0 (k) = α 1 x 1 (k) + α 2 x 2 (k) ( 17
)
where

α 1 , α 2 ≥ 0, α 1 + α 2 ≤ 1 Convex hull is contractive if and only if there exist {u 0 (k) ∈ U |x 0 (k + 1) ∈ γP 0 }. Applying the input, u 0 (k) = α 1 u 1 (k) + α 2 u 2 (k) (18)
the system dynamics (1) results in

x 0 (k+1) = A(α 1 x 1 (k)+α 2 x 2 (k))+B(α 1 u 1 (k)+α 2 u 2 (k)) (19a) x 0 (k + 1) = α 1 (Ax 1 (k) + Bu 1 (k)) + α 2 (Ax 2 (k) + Bu 2 (k)) (19b) x 0 (k + 1) = α 1 x 1 (k + 1) + α 2 x 2 (k + 1) (19c)
As x 1 (k + 1) ∈ γ 1 P 1 and x 2 (k + 1) ∈ γ 2 P 2 , therefore

x 0 (k + 1) ∈ γ(α 1 P 1 ⊕ α 2 P 2 ) ∈ γP 0 (20)

IV. SIMPLIFICATION OF CONTRACTIVE SETS

Assume that the contractive set obtained after taking the convex hull of different solutions provided by the PSO is P c with m hyperplanes and nv number of vertices. We want to minimize the number of hyperplanes so that the complexity of the formulated contractive set is reduced. As such a convex hull for higher dimensional systems can be complex, we'd like to reduce its complexity. Naturally, this complexity reduction should be obtained without reducing the volume of the contractive set substantially, and without increasing the contraction factor much. We introduce two methods to obtain the simplified contractive sets. The first method obtains the new set of vertices by using a circumscribed ellipsoid while the second method removes the hyperplanes so that the set P c remains contractive with required contraction factor within state constraints.

A. Circumscribed Ellipsoid

In this method we find a minimum volume ellipsoid ε c = { Cx + d ≤ 1} with center at d circumscribing the contractive set P c . The ellipsoid ε c can be found by minimizing the volume of ellipsoid containing all vertices of P c , this results in a maxdet problem as follows:

min C,d -log(det(C)) (21a)
subject to C = C T > 0 (21b)

CV i + d ≤ 1 (21c)
Where V i is a vertex of P c , and nv is the number of vertices of P c . Once the circumscribed ellipsoid is found, we look at which points from V i lie on the ellipsoidal boundary. Let us assume vertices V E lie on the boundary of the ellipsoid such that {V E ∈ V i , ∀i = 1, ..nv}. Then we rotate the ellipsoid such that point p f in V i , which was farthest from ellipsoidal boundary lies on new rotated ellipsoid. This adds the vertex which will give the largest increment in the volume. We find a normal plane to vector passing through p f and center of ellipsoid. There can be many directions in which such a plane can lie. We take projection ε p of the ellipsoid ε c and points V i on the plane and then check which new point is farthest from ε p . We add that point to current point p f and continue with projections. This procedure is repeated for n -1 times to obtain n -1 points. Normal to these n -1 points will give us a specific direction for rotation of ellipsoid (orientation of plane). We add two points on the new plane far away from original points and find a new minimum volume ellipsoid inscribing all the points. It makes sure that the point p f lies on boundary of the new ellipsoid. In this way, the ellipsoid will be rotated to find new set of vertices which originally belongs to V i , common points of the ellipsoid and the contractive set are added to V E . We repeat the rotation procedure until we find a large enough contractive set. V E forms the vertices of new simplified contractive set. There is a trade-off between volume of contractive set obtained and complexity reduced.

B. Removing a Hyperplane

The method described above finds a reduced complexity contractive set. As the method focuses on vertex operation for complexity reduction, there may arise a case when it increases the complexity in terms of hyperplanes, in that case the method described above should be discarded. Here, we propose a method for simplification of a polytope by directly operating on the hyperplanes. We remove the hyperplane if contraction factor of the simplified set doesn't exceed the required contraction constraint. We select the hyperplane H i to be removed by checking which hyperplane has to be pushed outwards least, in order to become redundant. The hyperplane H i can be removed if new vertices V obtained by removing H i are contractive and also fulfills the state constraints. The procedure is repeated for all the hyperplanes. As this method directly removes the hyperplanes, it guarantees complexity reduction and increase in volume of the contractive set. In worst case scenario, the procedure will not have any impact on the set unlike method described in previous sub-section. Results obtained with this procedure along with the formulations of previous method are illustrated in example VI-A.

V. DESIGNING A CONTROLLER WITH GIVEN CONTRACTIVE SET

Assume the low complexity controlled contractive set with contraction factor γ is obtained by the methods mentioned above. Then, the control formulation can be given as:

min u k x T k+1 Qx k+1 + u T k Ru k (22a) subject to x k+1 = Ax k + Bu k (22b) H u u ≤ h u (22c) F x k+1 ≤ γαf (22d)
where

α = max{F i x k,i /f i }, ∀i = 1, ..., m (22e)
The contractiveness ensures that the α will reduce by factor γ at each time step, therefore the state trajectories will converge to origin, which ensures stability of the system. The explicit solution to (22) can be obtained by solving it parametrically, with x k and α as the parameters. The complexity of the explicit solution (in terms of the number of critical regions obtained), is given by the number of different combinations of constraints that may be active at the optimum, when the parameters are allowed to vary throughout the given parameter region. MPC formulations with a long prediction horizon will have a high number of constraints, and thus also typically a high number of possible combinations of different constraints. The prediction horizon for the formulation in ( 22) is 1, and thus the complexity of the explicit solution can be expected to be low.

VI. EXAMPLES

A. Spring Mass Damper system

Consider the spring mass damper system example mentioned in [START_REF] Hovd | Low Complexity Constraint Control Using Contractive Sets[END_REF] with state representation given as: 1). The contractive set obtained by the method described in section III has 20 hyperplanes with volume of 40.9193 and contraction factor 0.9774 (Yellow set in figure 1). By simplifying the larger contractive set using method explained in section IV-A, the set obtained has a volume of 38.2663 with 0.9771 contraction factor and 6 complexity (as shown in green in figure 1). By implementing method described in section IV-B, simplified contractive set obtained has 6 hyperplanes with volume 42.4433 and contraction factor 0.9783 (blue set).

ẋ = 0 1 -7 -7 x + 0 1 u (23)

B. Identical Modes system

Consider a system with identical modes in series, i.e, nondiagonalizable system matrix A given as:

x k+1 = 0.98 1 0 0.98 x + 0 0 0.6 0.9 u (24)
The contractive sets obtained for this system are shown in figure 2. Contractive set obtained by the method explained in [START_REF] Dorea | a,b)-invariant polyhedral sets for linear discrete-time systems[END_REF] is shown in blue, while the green, red and yellow sets are obtained using the optimization based approach in this paper, when specifying different set complexities. It can be clearly seen that by decreasing complexity of a set, we also lose volume. Therefore, there is a trade-off between the volume of the contractive set and its complexity. (25) Input and state constraints are given as -1 ≤ u i ≤ 1, ∀i = 1, 2 and -1 ≤ x j ≤ 1, ∀j = 1, ...3 respectively. The contraction factor is chosen to be 0.9745. The explicit solution for a MPC with prediction horizon of 10 has 862 regions, and the total volume of the feasible region is 6.77. The contractive set for the system as explained in sections II is obtained with complexity of 8 and volume 2.0828. The controller is designed by method explained in section V. By using the set obtained from section II formulations, number of regions comes out to be 21. Clearly, the reduction in volume is significant, but so is the reduction in complexity for the explicit controller.

VII. CONCLUSION

A novel method for formulating controlled contractive sets has been described in this paper. An optimization based approach is used to find the contractive sets, subject to boundedness and contractiveness constraints. Multiple contractive sets can subsequently be merged to obtain a larger contractive set. Two techniques have been discussed to further simplify the set. The first technique utilizes circumscribed ellipsoid to find a reduced complexity contractive set while second technique proposes an operation on the hyperplanes to reduce the complexity. Finally, numerical examples shows the efficiency of the proposed method.

Fig. 1 :

 1 Fig. 1: Complexity Reduction of Contractive Set

Fig. 2 :

 2 Fig. 2: Flexible Complexity Contractive Sets

*The work leading to these results has received funding from the People Programme (Marie Curie Actions) of the European Union's Seventh Framework Programme (FP7/2007-2013) under REA grant agreement no 607957 (TEMPO).


Ying Fu
email: fu@ceremade.dauphine.fr

Herschel Rabitz
email: hrabitz@princeton.edu

Gabriel Turinici
email: gabriel.turinici@dauphine.fr

Hamiltonian identification in presence of large control field perturbations

Keywords:

Hamiltonian identification in presence of large control field perturbations

I. INTRODUCTION

The interaction between a quantum system and an applied field can be used to gather information about the system by measuring suitable observables for various incarnations of the field (see [1-3, 8, 14, 17, 18] for some related works). Such an inversion endeavor, seeking information about the system, has to take into account the possibility that the measurements are contaminated by noise, possibly from multiple sources, or that the model does not adequately describe all relevant characteristics of the system. Many considerations arise, including the following list:

S-1

The number of levels is unknown. Except for spins, the number of levels can generally take on any value, including the prospect of there being a continuum. The nature of the field bandwidth and intensity plays into which aspects of the system dynamics are accessible due to the field. In addition, unless the molecule is aligned (which is also possible) with the control field there will be vibration-rotation transitions. The circumstance does not change the key algorithmic principles set out in this work, but it can make the actual number of levels and dipole elements more involved.

S-2

Intense fields can be important in some problems where excited electronic states play a role, which can bring in non-resonant processes (i.e., virtual states, requiring adequate models). Moreover intense fields may lead to models including the nonlinear field coupling coefficients in the Hamiltonian.

S-3

If the system is in a thermal initial state, then the initial state will be a Boltzmann distribution instead of a single cold molecule in the ground or a pure state. The inversion will then need to deal with the distribution of molecules.

The performance of inversion needs to acknowledge the issues above. In addition, there will always be noise in the observations, but assuming that it is random and symmetric (e.g. Gaussian) then signal averaging may reduce its influence. However, noise may accompany the control field with the following circumstances: F-1 Different molecules (even if they are oriented) will see different fields as the field generally has a spatial (e.g., a Gaussian) lateral spread. F-2 A so-called spatio-temporal "chirp" may occur, where there is some unknown complex spatio-temporal pattern (possibly fluctuating) in the field. When the pattern is not uniform it has to be factored into the model with, possibly nonlinear, additional terms, see [START_REF] Gu | Spatial chirp in ultrafast optics[END_REF].

F-3

Although the field may be known before entering the sample, optical distortion during propagation can lead to additional field uncertainty in the domain where the data is taken.

F-4

The variations in field are stochastic from shot-to-shot with respect to the frequency dependence of the amplitudes and phases. The origin of this behavior is mainly due to the laser source and possibly jitter in the optics. The implications for the control data due to this occurrence is likely more significant at higher intensities.

Collecting all of these items pose a daunting, but not an impossible challenge to treat. In this work we limit the focus on the randomness coming from the last item F-4 above. These are related to so-called "fixed systematic errors", see [12, section VI.A. equations (38) and (40)] or "systematic control error", see [START_REF] Khodjasteh | Dynamically error-corrected gates for universal quantum computation[END_REF] in the quantum computing literature. In [START_REF] Alexandre | Experimental protection of quantum gates against decoherence and control errors[END_REF] "low frequency noise" is used (cf. also [10, section IV. C.]): it is the portion of the (control) amplitude noise that has a correlation time up to 10 3 times longer than the timescale of the dynamics therefore it can be considered as constant in time. The understanding of laser noise and various associated models is an active area of research [START_REF] Dufour | Optical noise reduction in a femtosecond ti: sapphire laser pumped by a passively stabilized argon ion laser[END_REF][START_REF] Ivanov | Experimental study of noise properties of a ti: sapphire femtosecond laser[END_REF][START_REF] Zola | The influence of laser field noise on controlled quantum dynamics[END_REF].

We present the theoretical background and analysis in Section II. The numerical algorithm and simulation results are given in Section III. Final remarks are made in Section IV.

II. THEORETICAL FRAMEWORK

A. The model

We introduce the following notation:

• L M 1 ,M 2 ,••• ,Mm is the Lie algebra spanned by the matrices M 1 , M 2 , • • • , M m ;
• For any matrix or vector X we denote by X * its adjoint (the transpose, complex conjugate);

• H N is the set of all Hermitian matrices H N = {X ∈ C N ×N |X * = X}; • Ψ(t, H, u(•), µ, Ψ 0 )
is the solution of the equation (1) below; to simplify the notation, when there is no ambiguity, we denote it as Ψ(t);

• SU (N ) is the special unitary group of degree N , which is the group of N × N unitary matrices with determinant 1;

• su(N ) is the Lie Algebra of skew-Hermitian matrices (the Lie algebra of SU (N ));

• A set of commuting observables O = {O 1 , ..., O K } (named SCO hereafter) satisfies [O k , O ] = 0, ∀k, ∈ {1, ..., K};
• L(X) is the distribution of the random variable X.

Let us consider the following controlled quantum system with time-dependent wavefunction Ψ(t) satisfying the Schrödinger equation:

       i Ψ(t, H, u(•), µ, Ψ 0 ) = (H -u(t)µ)Ψ(t, H, u(•), µ, Ψ 0 ) Ψ(0, H, u(•), µ, Ψ 0 ) = Ψ 0 , ( 1
)
where H is the internal ("free") Hamiltonian and µ the coupling operator between the control u(t) ∈ L 1 loc (R + ; R) and the system. We work in a finite dimensional framework, therefore H, µ ∈ H N for some N ∈ N * . The free Hamiltonian H is taken as known and the goal is to recover the matrix entries of µ from laboratory measurements of some observables depending on Ψ(t). The control u(t) can be changed in the laboratory to gather sufficient information on the system for the inversion to extract µ.

In the laboratory, the control field is produced by superposition of several frequencies, each with a specific amplitude and phase:

u(t) = σ(t) α =β A αβ sin(ω αβ t + θ αβ ), (2)
where σ(t) is a Gaussian envelope in time and ω αβ = E β -E α is the transition frequency between the eigenvalues E α an E β of H. The amplitudes A αβ and the phases θ αβ are the control parameters. In practice, the control fields cannot be produced perfectly. When we repeat nominally the same experiment, there can be a random shift in the amplitudes, which means for each experiment, there is a multiplicative noise factor on the amplitudes.

Accordingly the perturbation is modeled by a random multiplicative factor Y acting on the control i.e., u(t) is replaced by Y • u(t). This model has the assumption that each A αβ has the same shift; in addition we do not treat here possible noise appearing in the phases θ αβ .

In principle, the same methodology may apply to noise appearing in the phases; however at this time adapted ensemble controllability results are not available and this requires working with a high-dimensional probability distribution (one dimension for each independent noise in a phase). As a final remark, noise in the amplitudes and phases does not imply that the origin is in the pulse shaper but from the driving laser and represented in the fashion above.

The perturbation Y takes values in V = {y , ≤ L} ⊂ R. We denote ξ = P(Y = y ),

where the column vector ξ with entries ξ 1 , • • • , ξ L is a probability distribution on V. The numerical values of the possible perturbations y are known but their occurrence probabilities ξ are unknown and thus Ψ is a random variable, as are all measurements depending on Ψ. Repeating the control experiment several times will yield the distribution of the measurements.

The measurements are of the form

OΨ(T, H, u, µ, Ψ 0 ), Ψ(T, H, u, µ, Ψ 0 ) with O ∈ H N
being a member of a list of possible measurable operators. Often, only one observable operator is readily available, but the experiment can be repeated many times, including with distinct chosen fields. Generally for two observables

O 1 , O 2 ∈ H N no information is available on the joint distribution of the values O 1 Ψ(T, H, u, µ, Ψ 0 ), Ψ(T, H, u, µ, Ψ 0 ) and O 2 Ψ(T, H, u, µ, Ψ 0 ), Ψ(T, H, u, µ, Ψ 0 ) .

B. Theoretical result

The following theorem proves that under certain assumptions on the system and the SCO O, if we obtain the same distributions for all observables in O and for all controls, then the dipole moment µ and the probabilities (ξ ) L =1 can be identified up to some multiplicative phases.

Remark II.1. As the noise (with unknown distribution) Y multiplies the unknown dipole

µ, when the couple (Y, µ) is a solution, any couple (Y /λ, λµ) is a solution too; thus it is only possible to obtain Y and µ up to a multiplicative factor (here λ). If, for instance, some additional information on Y or µ is known the constant λ can be set accordingly. As such, we will suppose from now on that for at least one transition an absorption intensity measurement can be performed which provides the value of |µ k, | 2 for some given k, (see [4, Chapter XIII, Section C. [START_REF] Brif | Control of quantum phenomena: past, present and future[END_REF]

.b Fermi's Golden Rule, p. 1299]). This is recalled in assumption

Hyp-D below.

Theorem II.2. Let H, µ 1 , µ 2 ∈ H N , H diagonal, µ 1 = 0, µ 2 = 0, Y 1 , Y 2 two random
variables with values in the same set V, at least one of which is non-null, Ψ 1 0 , Ψ 2 0 ∈ S N some initial states and denote for a = 1, 2 and

u ∈ L 1 loc (R + , R): Ψ a (t, u) = Ψ(t, H, u(•), µ a , Ψ a 0 ).
Let O be a (non-trivial) SCO. We suppose that N ≥ 3 and:

Hyp-A L iH,iµ 1 = L iH,iµ 2 = su(N ); Hyp-B tr(H) = tr(µ 1 ) = tr(µ 2 ) = 0;
Hyp-C the eigenvalues of H are all of multiplicity one.

Hyp-D |(µ 1 ) k, | 2 = |(µ 2 ) k, | 2 = 0 for some fixed k, .
The final observation time is denoted T (assumed large enough) and we suppose the following equality of distributions:

L( OΨ 1 (T, uY 1 ), Ψ 1 (T, uY 1 ) ) = L( OΨ 2 (T, uY 2 ), Ψ 2 (T, uY 2 ) ) ∀u ∈ L 1 ([0, T ]; R), ∀O ∈ O, (3)
then for some (α i ) N i=1 ∈ R N :        (µ 1 ) jk = ±e i(α j -α k ) (µ 2 ) jk , ∀j, k ≤ N, P(Y 1 = y ) = P(Y 2 = ±y ) ∀ ≤ L. ( 4
)
Remark II.3. The ± signs in equation (4) are due to the hypothesis Hyp-D which leaves the sign undetermined; it may also appear when the sign of the noise is ambiguous i.e., when Y 1 and -Y 1 have the same distribution.

Remark II.4. The SCO O may contain just one observable.

Proof. The proof requires the tools introduced in [START_REF] Fu | Quantum Hamiltonian and dipole moment identification in presence of large control perturbations[END_REF] where the additive noise u(•) + Y was considered. In order to keep it simple we only give the main ideas and the modifications with respect to proof presented there. Using Lemma 5.1 in [START_REF] Fu | Quantum Hamiltonian and dipole moment identification in presence of large control perturbations[END_REF], for all ≤ L and O ∈ O,

there exists κ O ( ) ≤ L such that OΨ 1 (T, u • y ), Ψ 1 (T, u • y ) = OΨ 2 (T, u • y κ O ( ) ), Ψ 2 (T, u • y κ O ( ) ) . ( 5
)
Reasoning as in the proof of theorem 4.1 of [START_REF] Fu | Quantum Hamiltonian and dipole moment identification in presence of large control perturbations[END_REF], this implies that there exists W ∈ SU (N )

diagonal such that y κ O ( ) µ 2 = W y µ 1 W -1 . Since at least one of Y 1 or Y 2 is non-null, we can
suppose, without loss of generality, that y = 0; we deduce the existence of some λ ∈ R \ {0}

such that

       (µ 1 ) jk = λe i(α j -α k ) (µ 2 ) jk , ∀j, k ≤ N, P(Y 1 = y ) = P(Y 2 = λy ) ∀ ≤ L. ( 6
)
However, since µ 1 and µ 2 are fixed for all O and we obtain that λ is independent of , which, using assumption Hyp-D, gives the conclusion.

Remark II.5. The hypothesis Hyp-A is required for identification while Hyp-B is rather a convention. Hypothesis Hyp-C can be relaxed (as in [START_REF] Fu | Quantum Hamiltonian and dipole moment identification in presence of large control perturbations[END_REF]) if the O is a Complete Set of Commuting Observables (CSCO).

III. NUMERICAL RESULTS

A. The algorithm

Several numerical simulations were performed in order to illustrate the theoretical result in Section II. In all cases we simulate the "real" system with Hamiltonian H, dipole moment

µ real , noise Y real with distribution L(Y real ) = L =1 ξ real δ y and observables O = {O 1 , ..., O K } that correspond to the Hamiltonian H, specified by projections {|e 1 e 1 |, ... , |e N e N |}.
Here |e 1 is the j-th eigenstate of H. Note that here we take O = {O 1 , ..., O K } to be the entire set of projectors, but one would be enough for the theoretical result to hold, see the Remark II.4.

The measurements provide the distributions of the observables for each control u:

L k=1 ξ real k δ | Ψ(T,H,u•y k ,µ real ,Ψ 0 1 ),e j | 2 ). ( 7
)
The Hamiltonian H is assumed known to high accuracy. In contrast, we suppose that a priori we only know the order of magnitude of the dipole moment. This information is useful as an initial guess for the inversion procedure.

To find the dipole moment and the noise distribution of the control amplitudes we minimize the difference between (i) the observed distribution calculated, at some suitably large time T , with the current dipole candidate µ and the current noise distribution candidate L k=1 ξ k δ y k and (ii) the real observable distribution calculated with the real dipole moment µ real and the real distribution L k=1 ξ real k δ y k . The difference is summed over several controls u 1 ,...,u Nu and defined as:

J (µ, (ξ k ) L k=1 ; (u i ) Nu i=1 ) = log    1 N u Nu i=1 N j=1 W 1 L k=1 ξ k δ | Ψ(T,H,u i •y k ,µ,Ψ 0 1 ),e j | 2 , L k=1 ξ real k δ | Ψ(T,H,u i •y k ,µ real ,Ψ 0 1 ),e j | 2    . ( 8
)
Here W 1 represents the 1-Wasserstein (also known as Kantorovich-Rubinstein) distance between two distributions (see page 34-35 in [START_REF] Villani | Topics in optimal transportation[END_REF]); for two probability distributions Z 1 , Z 2 having cumulative distribution functions F Z 1 (respectively F Z 2 ) the distance is:

W 1 (Z 1 , Z 2 ) = 1 0 |F -1 Z 1 (x) -F -1 Z 2 (x)|dx. (9)
Other distances could also be used, e.g. W 2 . We start the optimization with an initial guess µ 0 (see ( 14)); the distribution ξ 0 is initialized to be uniform. The iteration n ≥ 1 consists in the following steps:

Algo 1 Randomly choose N u controls u n i , i = 1, ..., N u ; Algo 2 minimize ξ → J (µ n-1 , ξ; (u i ) Nu i=1
) and set ξ n to be a minimizer (in practice a close approximation);

Algo 3 minimize µ → J (µ, ξ n ; (u i ) Nu i=1
) and set µ n to be a minimizer (in practice a close approximation);

For the step Algo 2, denote a ij k = | Ψ(T, H, u n i • y k , µ n , Ψ 0 1 ), e j | 2 and b ij k = | Ψ(T, H, u n i • y k , µ real , Ψ 0 1
), e j | 2 . A part of the algorithm is to minimize the error:

log   1 N u Nu i=1 N j=1 W 1 ( L k=1 ξ k δ a ij k , L k=1 ξ real k δ b ij k )   . ( 10
)
The differential of the distance W 1 is not trivial to compute (see [START_REF] Villani | Topics in optimal transportation[END_REF] for a rigorous mathematical treatment); instead, in this step, we use the L 2 distance between the smoothed densities: each Dirac mass is replaced by a Gaussian distribution with small variance (here

ν 2 = 10 -4 ): δ a ij k is replaced by the normal distribution N (a ij k , ν 2 ) and δ b ij k is replaced by N (b ij k , ν 2 )
). The use of parameter ν does not assume any particular distribution for ξ real and does not bias towards one, it is only a rapid way to obtain a computable gradient.

We compute the

L × L matrix M ij whose entries M ij k, = 1 ν• √ 2π e - (a ij k -a ij ) 2 2•ν 2
account for the density of the distribution N (a ij , ν 2 ) at the point a ij k (and a similar matrix

M ij with M ij k, = 1 ν• √ 2π e - (b ij k -b ij ) 2 2•ν 2
for the real distribution). When ν is small the minimizer of the term in ( 10) is close to the minimizer of Nu i=1 N j=1 M ij ξ -M ij ξ real 2 which is given by the formula:

ξ n,raw =   Nu i=1 N j=1 (M ij ) T M ij   -1   Nu i=1 N j=1 (M ij ) T M ij   ξ real . ( 11
)
The term ξ n,raw is corrected to be a probability distribution and set to ξ n :

ξ n = 1 L =1 |ξ n,raw | (|ξ n,raw 1 |, ..., |ξ n,raw L |) . ( 12
)
For the step Algo 3 there is no explicit solution; a classical unconstrained nonlinear optimization algorithm is employed (we used the Gnu Octave procedure "fminunc").

Once the algorithm finished, the remaining overall multiplicative constant λ (see Remark II.1)) is set consistent with the data (see below the examples).

B. Numerical tests: N observables

We consider the 4-level system (N = 4) in [START_REF] Donovan | Exploring the Hamiltonian inversion landscape[END_REF] having:

H =          
0.0833 -0.0038 -0.0087 0.0041 -0.0038 0.0647 0.0083 0.0038 -0.0087 0.0083 0.0036 -0.0076 0.0041 0.0038 -0.0076 0.0357

          , µ real =           0 5 -1 0 5 0 6 -1.5 -1 6 0 7 0 -1.5 7 0           . ( 13
)
We set N u = 36; the controls are defined by formula (2) with amplitudes (A αβ ) chosen at random, uniformly in [0, 0.1 • H l ∞ µ real l ∞ = 0.0012] and the phases θ αβ in [0, 2π]. For the values ( 13) the initial guesses are:

µ 0 g =           0 3.
          , µ 0 e =           0 10 1 1 10 0 10 1 1 10 0 10 1 1 10 0           , µ 0 b =           0
          . ( 14
)
The average relative errors of these initial guesses are 42%, 70% and 50%. Guesses µ 0 g and µ 0

b are obtained by multiplying element-wise µ real by uniform random variables in [0.4, 1.6] (i.e., up to 60% average relative error); the term µ 0 e is obtained by taking roughly the order of magnitude of the entries of µ real : the values smaller or about 1 are taken to be 1 while the others are set to 10. The observables arising from the Hamiltonian H are projections • Y real = Y g being a Gaussian distribution centered at 1 with variance equal to 0.0025:

ξ real,g = f g L k=1 f g k = 1, • • • , L with f g k = 1 √ 0.0025 √ 2π e -(y k -1) 2 2•0.0025 , (16)
• Y real = Y e being a shifted exponential distribution form:

ξ real,e = f e L k=1 f e k = 1, • • • , L with f e k = 5 • e -5•(y k -ym) , (17)
• Y real = Y b being the bi-modal distribution which is the sum of two Gaussian distributions. We choose the first one centered at 0.8 with variance equal to 0.0025 and the second one centered at 1.2 with variance equal to 0.0049:

ξ real,b = f b L k=1 f b k = 1, • • • , L with f b k = 1 √ 0.0025 √ 2π e -(y k -0.8) 2 2•0.0025 + 1 √ 0.0049 √ 2π e -(y k -1.2) 2 2•0.0049 . ( 18
)
The dipole moments converge in 10 iterations. The numerical values are rescaled in order to use that |(µ real ) 12 | 2 = 25; same is done for the noise distributions. We obtain µ 10 g , µ 10 e and µ 10 b :

µ 10 g -µ real ∞ = 5 • 10 -5 , µ 10 g =           0 5 -1 0 5 0 5.99995 -1.5
-1 5.99995 0 6.99999 0 -1.5 6.99999 0

          , ( 19
)
µ 10 e -µ real ∞ = 10 -4 , µ 10 e =           0 5 -1 0 5 0 6 -1.5
-1 6 0 6.9999 0 -1.5 6.9999 0

          , ( 20
)
µ 10 b -µ real ∞ = 6 • 10 -5 , µ 10 b =           0 5 -1 0 5 0 5.99999 -1.5
-1 5.99999 0 6.99994 0 -1.5 6.99994 0

          . ( 21
)
The error norm is • ∞ is the largest, in absolute value, of the error components. See figures 1,2, 3 for the results. The tables I and II present the match of the probability distributions of the observables.

FIG. 1. Identification of the Gaussian distribution [START_REF] Villani | Topics in optimal transportation[END_REF]. The real distribution is in blue, the numerical result in red. Good agreement with the unknown noise distribution ξ real is obtained; the error on the dipole moment is 5 • 10 -5 , see equation [START_REF] Zola | The influence of laser field noise on controlled quantum dynamics[END_REF].

C. A single measured observable

We use the same system as in the Section III B except that here we consider the extreme case when only one observable is available as function of the control field. The observable is the projection |e 3 e 3 | to the third eigenstate.

The initial guess is:

µ 0 b =           0
          , ( 22
) Iteration n = 1 Real distribution O\u u 1 1 u 1 2 u 1 3 O 1 O 2 O 3 O 4 Initial guess O\u u 1 1 u 1 2 u 1 3 O 1 O 2 O 3 O 4 TABLE I.
The match of the observation distributions for the bi-modal distribution [START_REF] Zhu | Potential surfaces from the inversion of time dependent probability density data[END_REF]. We plot the histograms at the start of the inversion algorithm i.e., iteration step n = 1. In the left sub-table are the histograms of the observations

L k=1 ξ real k δ | Ψ(T,H,u 1 1 •y k ,µ real ,Ψ 0 1
),e j | 2 with the real distribution; in the right sub-table the histograms correspond to the initial guess which is the uniform distribution. Each column corresponds to a control field, here only the first 3 control fields u 1 1 , u 1 2 and u 1 3 are shown. Each line corresponds to a specific observable in the SCO set O. The initial guess is seen to be a poor approximation, as the histograms in left and right sub-tables differ substantially. The coefficient µ 2,1 = µ 1,2 = 5 is fixed, as previously stated, but now it is treated as a constraint by the algorithm (which will thus only optimize the other coefficients). The algorithm minimizes the difference:

J (µ, (ξ k ) L k=1 ; (u i ) Nu i=1 ) = log 1 N u Nu i=1 W 1 L k=1 ξ k δ | Ψ(T,H,u i •y k ,µ,Ψ 0 1 ),e 3 | 2 , L k=1 ξ real k δ | Ψ(T,H,u i •y k ,µ real ,Ψ 0 1 ),e 3 | 2    . ( 23
)
The distribution Y real tested is the bi-modal model. The algorithm converges and after 5,10 and 15 iterations, the dipole moments we obtain respectively have an L

IV. CONCLUSION

The inversion of the dipole moment has been considered in a model where the noise coming from the laser source is non-negligible and of unknown distribution. The model considered here has noise acting multiplicatively on the control intensity and is the same for all frequency components.

First, we proved theoretically that if one can measure repeatedly (at least) one observable for many control fields, the set of probability distributions of this observable is enough to recover both the dipole and the noise distribution. Then, a numerical algorithm based on the Wasserstein distance between the probability distributions was proposed and seen to perform well for several different, non-perturbative, noise distributions and initial guesses for the dipole. As subject for future study is the question of how to treat noise of possibly distinct character reflected in each of the amplitudes and phases or other different noise models as discussed in the Introduction. In addition, there are several further issues to consider in future work for creating a realistic algorithm for quantum system data inversion.

{|e 1 e 1 |. ( 15 )

 115 , |e 2 e 2 |, |e 3 e 3 |, |e 4 e 4 |} to eigenstates: |e 1 = 0.0845 -0.1313 0.9651 0.2101 T , |e 2 = -0.1305 -0.0856 -0.2103 0.9651 T , |e 3 = 0.2118 0.9647 0.0838 0.1325 T , |e 4 = 0.9649 -0.2118 -0.1314 0.0830 T In the basis {e 1 , e 2 , e 3 , e 4 } the Hamiltonian H is diagonal with eigenvalues E 1 = 0, E 2 = 0.0365, E 3 = 0.0651, E 4 = 0.0857. We set the final time T = 3200 which is about 10 periods of the smallest transition frequency 2π/(E 4 -E 3 ) = 314 in H. With respect to Remark II.1 we suppose that the value |(µ real ) 12 | 2 = 25 is known. The support for the distribution Y is known and denoted [y m , y M ]; in the numerical tests we take y m = 0.5 and y M = 1.5. We discretize the set of possible values of the perturbation with L = 51 equidistant points y = y m + ( -1) • y M -ym L-1 , = 1, ..., L. The values y are supposed known but not the probabilities ξ real that define the distribution L(Y real ) = L =1 ξ real δ y of the perturbation Y real . Several distributions are tested; they are constructed by discretizing, truncating and re-normalizing several classical distributions:

Final

FIG. 7 .

 7 FIG. 7. The L 2 error µ kµ real L 2 for k = 1, ..., 15 iterations.

TABLE II .

 II Converged result (iteration n = 10) from initial guess TableI. The identification works well as the left and right sub-tables match. As explained, the control fields are chosen randomly at each iteration and in particular the controls chosen at iteration n = 1 and n = 10 are not the same (otherwise the histograms corresponding to observations with the real noise distribution would be the same as in Table I left sub-table). Here the results for u10 1 , u 10 2 and u 10 3 are displayed.FIG.2. Identification of the exponential distribution[START_REF] Zhang | Quantum hamiltonian identification from measurement time traces[END_REF]. The real distribution is in blue, the numerical result in red. Good agreement with the unknown noise distribution ξ real is obtained; the error on the dipole moment is 10 -4 , see equation (20).FIG.3. Identification of the bi-modal distribution[START_REF] Zhu | Potential surfaces from the inversion of time dependent probability density data[END_REF]. The real distribution is in blue, the numerical result in red. Good agreement with the unknown noise distribution ξ real is obtained; the error on the dipole moment is 6 • 10 -5 , see equation (21).

  2 error of 0.17246, 0.03736 and 4.5652 • 10 -4 respectively. In order to test the robustness of the algorithm with respect to other error norms, we use the L 2 error norm • L 2 (the square root of the sum of squares of the components). The convergence is slower than in Section III B when we treat µ 2,1 = µ 1,2 = 5 as a constraint.FIG.4. Identification of the bi-modal distribution (18) after 5 iterations. The real distribution is in blue, the numerical result in red. The noise distribution starts to have the same qualitative features as ξ real ; the L 2 error on the dipole moment is 0.17246. FIG.5. Identification of the bi-modal distribution (18) after 10 iterations. The real distribution is in blue, the numerical result in red. The noise distribution starts to converge; the L 2 error on the dipole moment is 0.03736. FIG.6. Identification of the bi-modal distribution (18) after 15 iterations. The real distribution is in blue, the numerical result in red. Good agreement with the unknown noise distribution ξ real is obtained; the L 2 error on the dipole moment is 4.5652 • 10 -4 .

	0.08
	0.07
	0.06
	0.05
	0.04
	0.03
	0.02
	0.01
	0
	0.4	0.6	0.8	1	1.2	1.4	1.6
	0.08
	0.07
	0.06
	0.05 0.08
	0.04 0.07
	0.03 0.06
	0.02 0.05
	0.01 0.04
	0 0.03
	0.4	0.6	0.8	1	1.2	1.4	1.6
	0.02
	0.01
	0
	0.4	0.6	0.8	1	1.2	1.4	1.6
	Remark III.1. A legitimate question is related to the scaling of the identification with
	respect to the number N of levels. Note first that the computation in equation (23) only
	depends on the number of observables (and not on N , as equation (10) seemed to indicate).

On the other hand, the computation of the W 1 distance is independent on N and, for one dimensional laws, straightforward.

Of course in order to compute the observables, numerical simulations are performed, and these do depend on N . But such simulations allow for trivial parallelization which could bring down the wall-clock time per iteration to that of a single numerical resolution of the N -level system. Such a time is a lower bound because one needs to check whether a candidate solution is indeed a good solution.

ACKNOWLEDGMENTS

Y.F. and G.T. were partially supported by the "Agence Nationale de la Recherche" (ANR) Projet Blanc EMAQS number ANR-2011-BS01-017-01. H.R. acknowledges support from the U.S. National Science Foundation CHE-1464569.


Antoine Boutet
email: antoine.boutet@insa-lyon.fr

Albin Petit
email: albin.petit@insa-lyon.fr

Sonia Ben Mokhtar
email: sonia.benmokhtar@insa-lyon.fr

Léa Laporte
email: lea.laporte@insa-lyon.fr

Poster Abstract: Leveraging Query Sensitivity for Practical Private Web Search

Keywords:

Several private Web search solutions have been proposed to preserve the user privacy while querying search engines. However, most of these solutions are costly in term of processing, network overhead and latency as they mostly rely on cryptographic techniques and/or the generation of fake requests. Furthermore, all these solutions protect all queries similarly, ignoring whether the original request contains sensitive content (e.g., religious, political or sexual orientation) or not. Based on an analysis of a real dataset of Web search requests, we show that queries related to sensitive matters are in practice a minority. As a consequence, protecting all queries similarly results in poor performance as a large number of queries get overprotected.

In this paper, we propose a request sensitivity assessment module that we use for improving the practicability of existing private web search solutions. We assess the sensitivity of a request in two phases: a semantic sensitivity analysis (based on the topic of the query) and a request linkability analysis (based on the similarity between the current query and the query history of the requester). Finally, the sensitivity assessment is used to adapt the level of protection of a given query according to its identified degree of sensitivity: the more sensitive a query is, the more protected it will be.

Experiments with a real dataset show that our approach can improve the performance of state-of-the-arts private Web search solutions by reducing the number of queries overprotected, while ensuring a similar level of privacy to the users, making them more likely to be used in practice.

INTRODUCTION

Search engines have become an essential service for finding content on the Internet. However, by regularly querying these services, users disclose a large amount of their personal data. To limit the disclosure of users' personal in-Permission to make digital or hard copies of part or all of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. Copyrights for third-party components of this work must be honored. For all other uses, contact the owner/author(s). formation, many private Web search solutions have been proposed. There are mainly two types of approaches, the solutions obfuscating the queries or obfuscating the identity of the requester. For instance, Tor [START_REF] Dingledine | Tor: The second-generation onion router[END_REF] uses anonymous communication to hide the identity of the requester behind a chain of proxies, while GooPIR [START_REF] Domingo-Ferrer | h(k)-private information retrieval from privacy-uncooperative queryable databases[END_REF] generates k fake queries to hide among (k + 1) queries the real query of user. A recent study [START_REF] Petit | PEAS: Private, Efficient and Accurate Web Search[END_REF] show that each mechanism separately are not enough to efficiently protect users, and proposed to combine both mechanisms. However, this solution is costly in term of processing and network overheads. Indeed it relies on cryptographic operations and the generation of number of fake queries, which consequently induces loss of performance.

Middleware Posters and

Furthermore, existing private Web search solutions protect all queries similarly regardless the sensitivity of the query. Based on a real dataset of query log and a crowdsourcing campaign to collect the human judgments regarding the semantically sensitive queries, we showed that overall sensitive queries are in practice a minority. The content oblivious nature of existing private Web search solutions can lead to an over-protection of non sensitive queries and an under-protection or sensitive queries. In this work, we explore an adaptive query protection driven by the sensitivity of the query. Specifically, we have identified two types of sensitive queries. The first one are the semantically sensitive queries. While there is no universally adopted definition of what constitutes sensitive content, the privacy policy of most service providers tend to define sensitive personal information with similar categories (e.g., religious beliefs, political orientation, sexual life, or medical facts). Here, we consider a score to quantify the risk that a query is semantically sensitive compared to categories pre-selected as sensitive by its associated requester. The second type of sensitive queries are those which can be linked to a user profile (i.e., sensitive to a de-anonymisation attack), named linkability sensitive queries. To quantify the risk that a query is linkability sensitive, we consider a score reflecting the similarity between the query and past queries of the associated user.

In this poster, we present a sensitivity assessment module which dynamically protects users' queries according to their level of sensitivity. This sensitivity assessment takes into account two sensitive matters, the semantic content and the linkability of the query. To protect a given query, the proposed module combines obfuscation (i.e., fake queries) and anonymous communication. Each of these mechanisms is controlled by the actual level of sensitivity of the query, the more sensitive it is, the more protected it will get. More precisely, the number and the generation of fake query depends on the semantic risk analysis, while the number of proxies part of the anonymous channel depends on the linkability risk analysis.

SENSITIVITY ASSESSMENT

Our sensitivity assessment module analyses the risk that a query is semantically and linkability sensitive.

Semantic Risk Analysis

The semantic risk analysis provides a score in [0, 1] which determines how the query is related to the topics identified as sensitive by the user. To achieve that we use two libraries: (i) Wordnet, a lexical database, and (ii) eXtended WordNet Domains, a mapping of WordNet synsets to domain labels. First the query q is split into keywords. For each keyword, we collect all possible synsets (i.e., all possible meanings of one word in Wordnet) and return the most probable ones using a graph-based disambiguation method. Then, we map each synset to multiple categories using the eXtended Word-Net Domains library. Finally, we use all these mappings to compute and return the probability that the query is related to the categories defined as sensitive by the requester.

Linkability Risk Analysis

The goal of the linkability risk analysis is to provide a score in [0, 1] to determine if the query is sensitive to a deanonymisation attack. In such an attack, an adversary tries to link an anonymous query to a specific user by measuring the distance between the query and background knowledge collected for each user (part of their query history in our case). This attack is successful if the query is close enough to past queries of the user. To quantify the linkability risk, we first represent the query q in a vector space model that has for values either 0 (the term is not in the query) or 1 (the term is in the query). We then compute and return the cosine similarity between the vector associated to the query and the vector representing past queries.

DYNAMIC QUERY PROTECTION

Our sensitivity assessment module dynamically adapts the query protection. More precisely, this assessment controls both the obfuscation and the anonymous communication.

Query obfuscation is achieved through the generation of k fake queries. To dynamically define the value of k for each query, we use a linear projection between the score returned by the semantic risk analysis in [0, 1] and the number of fake queries in [0, kmax]. To dynamically define the number of proxies in the anonymous channel, named p, we use a linear projection between the score returned by the linkability risk analysis in [0, 1] and p in [1, pmax]. Therefore, depending on the level of sensitivity of the query, the protection can be achieved through only one proxy (i.e., for non sensitive queries), or through multiple fake queries and a chain of multiple proxies (i.e., for the most sensitive queries).

PRELIMINARY RESULTS

We evaluate our sensitive assessment module using a dataset of real queries from AOL search. Our experiments show that the semantic assessment lead to the identification of semantically sensitive queries with a recall of 87.9% while the linkability assessment classifies linkable queries with a recall of 94.9%. Furthermore, we show that adapting our module to existing private Web search solution (i.e., PEAS) reduces by 40.3% the number of queries overprotected and consequently significantly decreases the latency and the bandwidth consumption. Lastly, we show that this performance improvement is achieved while maintaining a similar level of privacy (i.e., sensitive information captured by the search engine) and utility (i.e., quality of the results presented to users).

BACKGROUND ON SENSITIVITY

Privacy policies of online services tend to more and more take into consideration the sensitive nature of the information. For instance, the privacy policy of Quora (a questionand-answer service) offers to users the possibility to choose whether to answer questions anonymously or with their names attached. Specifically, this privacy feature has been recently analysed in [START_REF] Peddinti | Cloak and swagger: Understanding data sensitivity through the lens of user anonymity[END_REF] to develop a classifier to predict the sensitivity of individuals posts. Another recent work [START_REF] Biega | R-susceptibility: An ir-centric approach to assessing privacy risks for users in online communities[END_REF] provides a metric to identify the most susceptible users who are exposed to a personal sensitive state (e.g., pregnant or afflicted by depression) from their posts created in online communities. To identify sensitive topics, LDA models are build from large data collections specialized in the considered sensitive state. The similarity between the information published by users and the most salient words associated to sensitive topics is then performed to identify the more exposed users.

CONCLUSIONS

We propose a sensitivity assessment module to dynamically adapt the protection of each query according to its level of sensitivity. The protection is achieved by combining obfuscation through fake queries and anonymous communications. The number of fake queries is controlled by the sensitivity level of the semantic content of the query, and the number of relay in the proxy chain is controlled by the sensitivity level of linkability of the query to its requester. Preliminary results show that the proposed solution provides an important performance improvement while maintaining a similar privacy and utility to state-of-the-arts solutions.

  Demos '16 December 12-16 2016, Trento, Italy c 2016 Copyright held by the owner/author(s). ACM ISBN 978-1-4503-4666-5/16/12. DOI: http://dx.doi.org/10.1145/3007592.3007595