Chunlei Yu
email: chunlei.yu@hds.utc.fr

Véronique Cherfaoui
email: veronique.cherfaoui@hds.utc.fr

Philippe Bonnifait
email: philippe.bonnifait@hds.utc.fr

An Evidential Sensor Model for Velodyne Scan Grids

For the development of driving assistance systems and autonomous vehicles, a reference perception equipment including navigable space determination and obstacles detection is a key issue. The Velodyne sensor which provides high definition and omnidirectional information can be used for this purpose. Nevertheless, when scanning around the vehicle, uncertainty necessarily arises due to unperceived areas and noisy measurements. This paper proposes an inverse evidential model for the Velodyne in order to exploit its measurements in a 2D occupancy grid mapping framework. The evidential sensor model interprets the data acquired from the Velodyne and successively maps it to a Carthesian evidential grid using a fusion process based on the least commitment principle to guarantee information integrity. Experimental results prove that this approach can handle efficiently the uncertainties of the sensor and thus a highly reliable local reference map near the vehicle can be built for every timestamped perception system that needs evaluation or calibration.

I. INTRODUCTION

To cope with errors and uncertainties when building occupancy grid maps, Bayesian methods are the foundations of usual frameworks. Many extensions have been published in the literature, like the Bayesian Occupancy Filter (BOF) [START_REF] Coue | Bayesian occupancy filteing for multitarget tracking : an automotive application[END_REF] which estimates both the occupancy and the speed of the cells. [START_REF] Lefaudeux | Extended occupation grids for non-rigid moving objects tracking[END_REF] proposed an extended occupancy grid approach which can be used to track non-rigid moving objects. [START_REF] Weiss | Robust driving path detection in urban and highway scenarios using a laser scanner and online occupancy grids[END_REF] applied the bayesian occupancy grid map to detect road boundaries. [START_REF] Pathak | 3d forward sensor modeling and application to occupancy grid based sensor fusion[END_REF] proposed novel forward model to interpret laser observations into occupancy grids. In this paper, we propose an evidential framework to build an occupancy grid map in the proximity of the host vehicle. The evidential theory [START_REF] Shafer | A Mathematical Theory of Evidence[END_REF] is used as the mathematical basis to build an inverse model. Indeed, by using this formalism, it is efficient to model state cell allocation and avoids introducing arbitrary a priori information in unoberved areas. The accumulation of scan echoes contributes to the detection of the ground and of the obstacles. Evidential grid mapping is not new. [START_REF] Pagac | An evidential approach to map-building for autonomous vehicles[END_REF] built an evidential grid map using the data acquired from a sonar, [START_REF] Yang | Evidential mapping for mobile robots with range sensors[END_REF] also adopted the evidential approach and built an evidential occupancy grid map in a Cartesian coordinates frame. Using the evidential framework to manage fusion is popular and [START_REF] Nashashibi | Vehicle recognition and tracking using a generic multi-sensor and multi-algorithm fusion approach[END_REF] applied evidential fusion rules to manage sensor uncertainties.

The V elodyne lidar [START_REF] Inc | HDL-64E S2 and S2.1 High Definition LiDAR Sensor User's Manual and Programming Guide[END_REF] provides rich and accurate information about the surrounding environment, an adapted sensor model to tackle its uncertainty and to fully profit its rich information is put forward. We propose a tailored sensor model which interprets the V elodyne data into local 2D occupancy grid maps. The least commitment principle is adopted here in order to avoid introducing any prior information. The model makes fully use of the evidential theory to handle the rich information provided by the V elodyne while keeping the processing load reasonable. A fusion process based on Dempster Shafer data fusion enhances the reference map by fusing data acquired at different locations.

The paper is organized as follows: section II presents the evidential framework, the theory basis will be illustrated. Section III details the main contribution of the work. An evidential sensor model is developed and utilized to merge high definition lidar measurements into scan grid maps. Section IV illustrates the fusion scheme based on the evidential framework. Section V shows the implementation details and experimental results.

II. EVIDENTIAL FRAMEWORK

In this framework, a frame of discernment Ω is defined to model the state of each cell. The frame of discernment is a finite discrete set, which contains all mutually exclusively propositions of interest. These basic propositions are called singleton propositions. In our case, the frame of discernment is defined as: Ω = {O, F }, the two singletons are the proposition O and the proposition F, indicating respectively that the specific cell is occupied and free. One has to increase this set, by considering the power set which is defined as 2 Ω = {∅, F, O, Ω}. With the frame of discernment, there are 2 |Ω| possible subsets, where |Ω| is the cardinality of the set. Each subset is a possible proposition and in this manner it is possible to exhaustively propose other more general propositions based on interactions between the singleton propositions.

The meaning of each singleton proposition is detailed below:

• 1) O indicates Occupied cell. • 2) F indicates F ree cell.
• 3) Ω indicates ignorance about the state of the cell (U nknown cell). • 4) ∅ indicates that no proposition fits the cell. In evidential theory, different functions may be used to represent the information in the power set. One specific function which is mostly utilized in our approach is the basic probability assignment (BP A), which is a direct support for a specific proposition. Let function m denote a BP A that maps each proposition in 2 Ω to a numerical measure of direct support for that particular proposition. The function m returns values in the range of [0, 1] and satisfies the conditions:

A⊆Ω, A =∅ m (A) = 1, m (∅) = 0
One powerful application of evidential theory is the fusion of different sources of information. Let m 1 and m 2 be two given mass functions describing the occupancy belief of the same cell. The result of the combination using Dempster's conjunctive rule is computed in the following way:

(m1 ⊕ m2) (A) = K ∀B,C∈2 Ω , B∩C=A, A =∅ m1 (B) m2 (C) K -1 = 1 - ∀B,C∈2 Ω , B∩C=∅ m1 (B) m2 (C) (m1 ⊕ m2) (∅) = 0
In the mapping process, one inherent aspect is to merge sensor data from different locations at different time instants into the built map. This fusion process changes the space states from unknown (initial state) to another applicable state and enables to update and confirm the state of the cells. Typically, the Dempster's conjunctive rule is chosen as the fusion operator.

III. FROM VELODYNE DATA TO EVIDENTIAL SCAN GRID MAP

In this section, an inverse sensor model for the V elodyne is developed. The sensor model has to be chosen with special attention since a well adapted model can greatly increase the perception performance. In our approach, we build a 2D evidential occupancy grid map with data from the Velodyne by making a projection on the ground plane. To develop the sensor model, we suppose that the ground is locally flat for every Velodyne scan. In order to be as close as possible to the sensor's rotating acquisition process, the scan grid map is created in a polar frame. As shown in Figure 1, the whole space around the car, which is placed at the origin, is divided into angular sectors, while each sector in the space is divided into different cells. For the BPA assignment process, we consider the sectors independent from each other. Indeed, if the sampling of the grid is high enough and since the laser beam width is very small, this assumption is well verified. In Figure 1, the colors represent different information. The dark cells mean that there is no information in the corresponding space due to an occlusion because of obstacles for instance. Green and red cells represent respectively F ree and Occupied space.

A. Polar sensor model basic concepts

This state assignment respects the least commitment principle. Velodyne points determine the state of the corresponding cells. Therefore, the space where there is no information is treated as U nknown which avoids adding wrong a priori information.

Although Velodyne is a high definition sensor, there might exist objects that are not detected. In Figure 2, six beams are drawn in the angular sector corresponding to angle λ, among which three hit the ground (G1, G2, G3), and the three other ones hit obstacles (H1, H2, H3). Due to the installation of the sensor on the roof, there is little information provided about the space between the host vehicle and the first impact on the ground, and between the consecutive impacts. As shown in the figure, there exists a non-null possibility that some nondetected objects exist in these areas, like O1, O2, O3 and O4. Therefore, when projecting Velodyne data in 2D, one needs to consider with attention what can happen between two consecutive beams. If the sensor is used for autonomous driving or driving assistance, this is crucial for safety.

B. Thresholding the ground

Let define an elevation threshold denoted H which specifies the elevation of points considered as obstacles. The value of the threshold has to be chosen carefully in order to filter noise. When the elevation of an echo is above H, we consider the cell Occupied.

Figure 3 illustrates the state allocation process, and a specific sector is shown. To differentiate the ground information from the above-ground information, we set up a scene where a human and a car are near our host vehicle. Corresponding to the angular sector of angle λ, the lateral view is shown, where eight beams from Velodyne are drawn for illustrating purpose. The four lowest green beams hit the ground, where G1, G2, G3 and G4 are respectively the intersections of the signal and the ground. The four red beams reach the human and the car in the distance, and their intersections are H1, H2, H3 and H4. The grid in the bottom serves as an illustrating plot of the polar world model shown in Figure 3, in which the horizontal axis shows the range variation and the vertical axis represents the angle variation.

G1 G2 H1 H2 G3 H h H3 G4 H4 Range Angle O U U U F F F U F F U U U U U U 1 2
Based on the least commitment principle, the state allocation process obeys the following rules: As shown in Figure 3, the cells which contain the H2, H3 and H4 are marked O(Occupied), as these points are above the threshold; the cells which contain respectively G1, G2, G3 and G4 are marked as F (F ree), as these points are detected on the ground; as depicted above, although H1 is above the ground and hits on the person, the corresponding cell is allocated F ree because it is under the threshold H. All the other cells are marked U (U nknown). One benefit of defining the threshold H is the extension of the F ree region by making a backward extrapolation to the host vehicle. Figure 4 illustrates the principle: considering the beam which hits the ground, one can deduce that there is no obstacle in the interval L which has an elevation superior to threshold H. In this case, the zone F ree corresponding to G1 is extended towards the host vehicle. Two more cells covered by the extension distance L are set to F ree. In Figure 3, the states of cell 8 cell 11 and cell 13 at angle λ are also set to F ree because we extrapolate at point G3, G4 and H1.

C. Backward free space propagation

G1 H h extrapolation α L Range Angle U U U U F U U U U λ F F U U U U U U U U U

D. Eliminating conflicting impacts

One conflicting situation might happen using the above sensor model, as shown in Figure 5, where the host vehicle detects one bus ahead, four beams are drawn for illustration. Based on the cell state allocation principle, the corresponding cells for G1 and G2 should be allocated F ree, which conflicts with reality. Obviously, cell 10 should be allocated Occupied and cell 11 should be allocated U nknown based on the sensor readings. In our approach, we make an additional modeling: a detected obstacle on the ground is modeled as a vertical surf ace that is linked to the ground. So in cells where both obstacle points and ground points are detected, ground points are ignored. This model brings the convenience of eliminating the conflicting cells, but at the expense of making impossible to differentiate hung obstacles from obstacles linked to the ground. The ground points detected behind the first obstacle surface are ignored and the corresponding cells rest U nknown until future exploration when the host vehicle moves.

E. Grid mass assignment

We need now to assign a BP A to the grid cells to quantify the belief. We propose a grid mass assignment model based on information accumulation. In Figure 3, cell 6, cell 9 and cell 12 are all set to F ree. However, we should have unequal amount of beliefs about their F ree state because in cell 6, there exist two points on the ground to support the state, whereas in cell 9 and cell 12, there exists only one. The same stands for the Occupied cells, cell 15 and cell 19 should have unequal amount of beliefs about the occupied state. M ore points supporting one state should contribute to more belief s on the state. This accumulation concept reinforces the belief assigned to each proposition. The BP A values are based on sensor uncertainties. Let α F A and α M D correspond to the the probability of false alarm and missed-detection. A false alarm is when the sensor issues an impact whereas there is nothing. It depends essentially on the sensor noise and on multipath propagation. A missed detection is mainly related to the reflexivity of the target and to the ratio between the cell size and the beam width. Figure 6 shows how this ratio results in missed-detection. δ s represents the divergence of Lidar (beam width), and δ g represents the angular resolution of the polar grid. In this circumstance, one beam of Lidar can not cover the whole sector. This beam can miss potential obstacles within its blind regions of the cell. The missed-detection effect thus has to be considered.

The proposed model calculates the BP As with probabilistic approach. Based on the definition of false alarm, α F A = P (C = F | ξ 1 ), Where ξ 1 represents one obstacle impact in the cell, C stands for the state of the cell. Supposing that errors are independent, the total false alarm probability in one cell given n O obstacle points are detected in this cell should be

P (C = F | ξ 1 , ξ 2, ... , ξ N ) = α n O F A .
Thus the probability of Occupied can be represented as:

P (C = O | ξ 1 , ξ 2 , ..., ξ N ) = 1 -α n O F A .
Based on the same methodology, for F ree cells, the missed-detection probability α M D = P (C = O | ∆), where ∆ represents no above ground impact is returned to the sensor. If we assume n F ground points are detected in this cell, the total missed-detection probability should be α n F M D . Thus the probability of F ree should be represented as 1 -α n F M D . Based on the principle, the BP A assignment: For a F ree cell:

m(O) = 0, m(F ) = 1 -α n F M D , m(Ω) = 1 -m(F ), m(∅) = 0
For an Occupied cell:

m(O) = 1 -α n O F A , m(F ) = 0, m(Ω) = 1 -m(O), m(∅) = 0
For an U nknown cell, the initial state is kept:

m(O) = 0, m(F ) = 0, m(Ω) = 1, m(∅) = 0
To keep the processing load reasonable, we suggest to extrapolate the free level m(F ) uniformly with no decrease to the cells that have no echoes.

F. From polar to cartesian

The approach merges the Velodyne scan data into occupancy grid map. This map is built in polar coordinate system, but for the fusion purpose, we need to transform it into Cartesian coordinate system. All the information collected has to be transformed into Cartesian coordinates, with the least loss. Several methods exist to do the transform. In our approach, we have adopted the bilinear interpolation algorithm introduced by [START_REF] Moras | Evidential perception grids for robotics navigation in urban environment[END_REF].

IV. EGO-MAP GRID FUSION SCHEME

The scan grid map is not complete, because there exist uncertainties in the map due to unperceived space. The U nknown space between the host vehicle and the obstacles should be eliminated. With Dempster's conjunctive rule, the fusion of several successive scan grids allows to gradually eliminate the uncertainties in the map. In the proposed approach, a map is built by fusing the scan grids centered at different locations. To make this fusion, the ego-motion of the host vehicle has to be compensated and then every new scan grid of the Velodyne is merged into a grid denoted EgoM apGrid. Figure 7 illustrates the whole approach. The fusion process is sequential. At time t, the new scan grid ScanGrid(t) updates EgoM apGrid(t -1) to provide a new EgoM apGrid(t). To accommodate to the dynamic environment, we adopt the approach proposed by [START_REF] Moras | Moving objects detection by conflict analysis in evidential grids[END_REF] and use a decay factor for EgoM apGrid. The information in EgoM apGrid can become aged and not consistent with reality. This effect can be especially important when moving objects are in the scene. Adopting decay factor means that all the information in the past will gradually become less evident. The equations below show how mass functions are discounted with a decay factor denoted β.

β mM (A) = β * mM (A), A ⊂ Ω β mM (Ω) = 1 -β + β * mM (Ω)
The fusion process adopts the Demster-Shafer conjunctive rule, as shown in Equation 1. For denotation purpose, let m M,t and m S,t represent respectively the mass functions of EgoM apGrid and ScanGrid at time t. The approach was tested with a vehicle of the experimental platform PACPUS shown in Figure 8. We have implemented the approach in C++. In the approach, the V elodyne HDL -64E data was acquired at 10Hz frequency. The segmentation process detects the ground plane and above ground points. We have adopted the ground labeling technique proposed by [START_REF] Petrovskaya | Model based vehicle tracking for autonomous driving in urban environments[END_REF], which generally provides satisfactory results. The ego-motion between two scans is estimated using CAN data.

mM,t = β mM,t-1 ⊕ mS,t (1)

V. EXPERIMENTAL IMPLEMENTATION AND RESULTS

A. Experimental implementation

For the purpose of demonstration, the scan grids of (72 * 72) meters are built with uniform cells of size(0.1 * 0.1) meters. In the polar grid map, the angular resolution is 0.5 degrees and the radius resolution is 0.1 meters. For the tuning parameters, we have adopted α M D = 0.66, α F A = 0.15. α M D is based on the ratio of the beam divergence of V elodyne (estimated to 0.17 degrees by [START_REF] Lonjaret | Rapport d'essais de qualification de l'imageur laser[END_REF]) and the resolution of the grid (0.5 degrees). We have tuned empirically α F A to 0.15 in order to consider the sensor noise and the multipath phenomenon.

B. Results

Figure 9: 3D display of the scene

For every instant, a Velodyne scan grid is built. For illustration of the system, we report in this section a typical scene where the host vehicle is in an urban road, as shown in Figure 9. The scene is depicted using Point Cloud Library (P CL) [START_REF] Rusu | 3d is here: Point cloud library (pcl)[END_REF]. Three occupancy maps are shown in Figure 10. Figure 10a shows result with no backward extrapolation while Figure 10b and Figure 10c show scan grids with two different H values for comparison, respectively H = 0.2 and H = 0.4. Conform to the sensor model, with a larger H value, longer distance is extrapolated, so larger F ree zones are observed in the grid map in which H = 0.4. Augmentation of H can reduce the uncertainties in the map, as less unperceived space is contained in the map thanks to backward extrapolation, but it also adds uncertainty to the system, because obstacles points with high elevation (but still beneath H) can be recognized as ground points. In reality, the vehicle can cross speed-bumpers or small slopes, with H set appropriately, these situations pose no problem and the approach can be more robust.

The effect of information accumulation model is reflected by the green level, which corresponds to the F ree mass degree. From the 3D display in Figure 9, the ground points received by V elodyne are typically denser in the space near the host system than farther away, which contributes to higher F ree belief in the nearby space due to more ground points. Thus in the resultant grid map in 10c, the F ree (green) space near the center is brighter than at the edge. Figure 11 shows the fusion result. Corresponding to the scene depicted in Figure 9, H is set to 0.2. One can remark that the green level in the central part of the fusion result map is lower compared to the other surrounding parts of the F ree space. The reason is that the sensor receives no information from this space in the present scan. The resulting scan grid in Figure 10b shows that this space is U nknown, which makes no contribution to state confirmation in the EgoM apGrid, whereas the F ree space surrounding this part accumulates in the F ree mass from the ScanGrid, resulting to a higher F ree mass degree in fusion.

The decay factor was set to 0.98 to slowly discount aged information. This effect is noticeable in the fusion result map shown in Figure 11a: the right bottom part of the map shows darker green which means less evidence to be F ree. This can also be explained by the ScanGrid of Figure 10b: the state of this space is U nknown in the ScanGrid. With no evidence supporting the space state, the system tends to gradually forget its past state. In this case, m(F ) decreases until the system totally forgets the state, and it becomes U nknown again.

Figure 11b shows the result of the decision rule consisting in selecting the cell state which contains the maximum mass. This is the resultant reference map.

VI. CONCLUSION

In this paper, we have proposed a new inverse sensor model to transform raw V elodyne data into local evidential grid maps. Based on the least commitment principle, the proposed model avoids introducing incorrect prior information in the grid state allocation of the space. The proposed way to fill the unobserved areas is based on the principle of information accumulation thanks to the fusion of successive scans taken at different points of view. This is performed by a spatiotemporal fusion process based on evidential theory. Based on the real experiments we have done to test the approach, we have observed that the approach provides good results in comparison to real scenes from a qualitative point of view.

Figure 1 :

 1 Figure 1: Space representation in Polar Coordinates, showing howmeasurements from V elodyne can be interpreted in the evidential framework. Green refers to free space, red refers to occupied space and dark refers to unknown space.

Figure 2 :

 2 Figure 2: Some lidar beams in a particular angular sector S. Top, bird view of the host vehicle, x represents the motion direction. Bottom, lateral view of the scene.

Figure 3 :

 3 Figure 3: Sensor model for the Scan grid construction process. Top, lateral visualization of the threshold scene and of the backward free extrapolation (short vertical lines, explained in section III-C). Bottom, state assignment.

Figure 4 :

 4 Figure 4: Backward extrapolation to the host vehicle: the two cells included in the interval L are extended as F ree space.

Figure 5 :

 5 Figure 5: Conflicting information elimination

δs

  Beam divergence of Lidar projected on ground δg Angle resolution of polar grid

Figure 6 :

 6 Figure 6: Missed-detection illustration

Figure 7 :

 7 Figure 7: Workflow of the scan grid construction and fusion

Figure 8 :

 8 Figure 8: Experimental platform with the Velodyne sensor

4 Figure 10 :

 410 Figure 10: Scan Grids, Green represents F ree space, Red represents Occupied space, Dark represents U nknown space.


Marta Conde Pena

Jean-Charles Faugère
email: jean-charles.faugere@inria.fr

Ludovic Perret
email: ludovic.perret@lip6.fr

Algebraic Cryptanalysis of a Quantum Money Scheme The Noise-Free Case

We investigate the Hidden Subspace Problem (HSPq) over Fq:

Input : p1, . . . , pm, q1, . . . , qm ∈ Fq[x1, . . . , xn] of degree d ≥ 3 (and n ≤ m ≤ 2n). Find : a subspace A ⊂ Fq n of dimension n/2 (n is even) such that pi(A) = 0 ∀i ∈ {1, . . . , m} and qj(A ⊥ ) = 0 ∀j ∈ {1, . . . , m}, where A ⊥ denotes the orthogonal complement of A with respect to the usual scalar product in Fq. This problem underlies the security of the first public-key quantum money scheme that is proved to be cryptographically secure under a non quantum but classic hardness assumption. This scheme was proposed by S. Aaronson and P. Christiano [1] at STOC'12. In particular, it depends upon the hardness of HSP2. More generally, Aaronson and Christiano left as an open problem to study the security of the scheme for a general field Fq. We present a randomized polynomial-time algorithm that solves the HSPq for q > 2 with success probability ≈ 1 -1/q. So, the quantum money scheme extended to Fq is not secure. Finally, based on experimental results and a structural property of the polynomials that we prove, we conjecture that there is also a randomized polynomial-time algorithm solving the HSP2 with high probability. To support our theoretical results, we also present several experimental results confirming that our algorithms are very efficient in practice. We emphasize that [1] proposes a non-noisy and a noisy version of the public-key quantum money scheme. The noisy version of the quantum money scheme remains secure.

Proposition 1 (Section 4). Let p = (p 1 , . . . , p m ), q = (q 1 , . . . , q m ) ∈

then we can construct an algebraic system of equations Sys HSP 2 over F 2 in N = n 2 /4 variables with at most 2m n/2 1 + n/2 2 + . . . + n/2 d such that a systematic basis of A

Introduction

The no-cloning theorem in quantum mechanics states the impossibility of creating identical copies of an unknown arbitrary quantum money state. In [START_REF] Wiesner | Conjugate coding[END_REF], Wiesner suggested to take advantage of this physical law in order to construct a scheme for (quantum) money that could not be counterfeited. The initial work of Wiesner has been then followed by several papers that try to improve the initial idea of [START_REF] Wiesner | Conjugate coding[END_REF], i.e. [START_REF] Bennett | Quantum cryptography, or unforgeable subway tokens[END_REF][START_REF] Gavinsky | Quantum money with classical verification[END_REF][START_REF] Mosca | Quantum coins. Error-Correcting Codes[END_REF]. This line of research culminated with the proposal of Aaronson and Christiano [START_REF] Aaronson | Quantum money from hidden subspaces[END_REF] at STOC '12 who proposed a public-key quantum money scheme.

A public-key quantum money scheme is a scheme in which anyone with a quantum device can verify if a banknote is valid rather than only the bank that issued (in contrast to [START_REF] Wiesner | Conjugate coding[END_REF]). A public-key quantum money scheme based on knot theory was introduced in [START_REF] Farhi | Quantum money from knots[END_REF]. However, its security is not well understood. The scheme proposed by Aaronson and Christiano in [START_REF] Aaronson | Quantum money from hidden subspaces[END_REF] is the first that is public-key and proved to be cryptographically secure under a classical (as in non-quantum) hardness assumption. The scheme is based on hiding two orthogonal subspaces by expressing each of them as the common zeros of a set of appropriate random multivariate non-linear polynomials. In particular, its security relies on the assumption that the following problem is hard:

Hidden Subspaces Problem (HSP q ) Input : polynomials p 1 , . . . , p m , q 1 , . . . , q m ∈ F q [x 1 , . . . , x n ] of degree d ≥ 3, n ≤ m ≤ 2n. Find : a subspace A ⊂ F q n of dimension n/2 (n is even) such that p i (A) = 0 ∀i ∈ {1, . . . , m} and q j (A ⊥ ) = 0 ∀j ∈ {1, . . . , m},

where A ⊥ denotes the orthogonal complement of A with respect to the standard scalar product in F q .

We emphasize that in [START_REF] Aaronson | Quantum money from hidden subspaces[END_REF], the authors propose a non-noisy and a noisy version of the public-key quantum money scheme. In this paper we only consider the noise-free version of the quantum money scheme.

In particular, the non-noisy version of the quantum money scheme relies on the HSP 2 , and Aaronson and Christiano conjecture that it cannot be solved in polynomial-time. They also state as an open problem the study of the scheme extended to a general field F q , which brings up the question of the hardness of HSP q .

We analyze the hardness of the HSP q . The main idea is to model the problem as a set of algebraic equations. Expressing elements as the common zeroes of a set of random multivariate non-linear polynomials is the core of algebraic attacks, e.g. [START_REF] Faugère | Algebraic cryptanalysis of hidden field equation (HFE) cryptosystems using gröbner bases[END_REF][START_REF] Courtois | Algebraic attacks on stream ciphers with linear feedback[END_REF][START_REF] Faugère | Polynomial equivalence problems: Algorithmic and theoretical aspects[END_REF]. However, in this case we can exploit that there are two sets of public polynomials whose sets of zeros are two subspaces orthogonal to each other.

Aside from this quantum money scheme, the HSP q has also interest as a general computer algebra problem closely related to the isomorphism of polynomials [START_REF] Patarin | Hidden fields equations (HFE) and isomorphisms of polynomials (IP): two new families of asymmetric algorithms[END_REF]. Given p = T • p ′ • S, where p = (p 1 , . . . , p m ), p ′ = (p ′ 1 , . . . , p ′ m ) ∈ F q [x 1 , . . . , x n ] m and T, S are affine invertible transformations, the Isomorphism of Polynomials (IP) problem consists on recovering T and S. The HSP q can be seen as a slight modification of the isomorphism of polynomials problem where p = 0, T is the identity transformation and S is linear but not invertible.

Main Results

Our results mostly rely on Gröbner bases and linear algebra techniques. This is because we are capable of identifying the solution of the HSP q as the unique solution of an overdetermined system of multivariate equations in N = n 2 /4 unknowns (Section 3.1, Proposition 5). The properties of this system are different for q = 2 and q > 2, so we study separately both cases.

Our first main result (Section 3) solves an open problem presented in [START_REF] Aaronson | Quantum money from hidden subspaces[END_REF], which is the study of the HSP q for q > 2. From the algebraic equations describing HSP q , we observe that we can extract a set of linear equations (Lemma 3). Due to the shape of the linear equations, we can prove prove that sufficiently many linearly independent ones can be extracted. This gives: Theorem 1 (Section 3.2). Let N = n 2 /4. There is a randomized polynomial-time algorithm solving HSP q , for q > 2, with complexity O(N ω ) = O(n 2ω ), where 2 ≤ ω ≤ 3 is the linear algebra constant, and success probability γ q (n/2)γ q (m) γ q (m -n/2) , γ q (k) being the probability that a random k × k matrix with entries in F q is invertible. For n big enough, the success probability is ≈ 1 -1/q.

In Section 3.3, we report experimental results demonstrating that HSP q , with q > 2, can be solved very efficiently. For n ≤ 20, the algorithm requires less than 0.1 s. for various q. We have implemented the algorithm using the Magma software [START_REF] Bosma | The Magma algebra system I: The user language[END_REF]. The code is provided with the submission so that the results can be reproduced or conducted for bigger values of n.

Our second result is concerned with the HSP 2 (Section 4). In this case our system does not contain, except with a small probability, linear equations and so the approach needs to be different. Still, in the case of the HSP 2 , we have an algebraic system of equations which is very overdetermined.

vanishes Sys HSP 2 with probability γ 2 (n/2), where γ 2 (n/2) denotes the probability of a random matrix with entries in F 2 of being invertible. For n big enough, this is ≈ 1/2. So, we can still hope that the computation of a Gröbner basis of the system will be efficient. In Section 4.1, we run experiments to confirm this intuition. It appears that Sys HSP 2 is much easier to solve than a semi-regular system of the same size. The Magma code of this part is also provided with the submission. Typically, we can solve in practice Sys HSP 2 for n ≤ 18 and d = 3 in less than 3 hours (for smaller n, we can solve in few minutes). For n = 18, we have to solve a system of degree-3 equations with 81 variables. In practice, we observed that the maximum degree reached during the computation of a Gröbner basis of Sys HSP 2 is bounded from above by a small constant. Based on this observation we conjecture then that: Conjecture 1. The degree of regularity is bounded above by d + 1. If this conjecture is true, the following result is obtained:

Theorem 2. Let N = n 2 /4.
There is a randomized polynomial-time algorithm solving HSP 2 with a complexity of O(N ω(d+1) ) = O(n 2ω(d+1) ), where 2 ≤ ω ≤ 3 is the linear algebra constant, and success probability γ 2 (n/2), where γ 2 (k) denotes the probability of a random k × k matrix with entries in F 2 of being invertible. For n big enough, the success probability of the algorithm is ≈ 1/2.

To support our assumption we analyze in Section 5 the structure of Sys HSP 2 . We prove a structural property (due to the orthogonality of the hidden subspaces) that allows to obtain equations of degree lower than d from the public polynomials Sys HSP 2 by performing simple manipulations on the initial system. In particular: Proposition 2. We can easily generate O(m 2 ) equations of degree d -1 which are linear combinations of the equations from Sys HSP 2 .

This means that a Gröbner computation on Sys HSP 2 will generate at the very first step many equations of lower degree. This is known as a fall of degree and it is typically a behaviour which is not occurring in a random (i.e. semi-regular) system of equations. So, it is a first step towards proving our conjecture.

Organization of the Paper

In Section 2 we introduce some notation that will be used throughout the paper, we recall the basics of Gröbner bases and we describe the non-noisy version of the quantum money scheme of [START_REF] Aaronson | Quantum money from hidden subspaces[END_REF]. The first part of Section 3 is concerned with the general modeling of the HSP q as a system of multivariate non-linear equations, and the second part is dedicated to obtain the algorithm of the first Theorem. Sections 4 and 5 are dedicated to the HSP 2 . In Section 4 we explain precisely why the behaviour of HSP q is different for q = 2 and q > 2, we report experimental results and we derive our conjecture which, if true, results in the second Theorem. Section 5 is the most technical one, in which we explain that equations of degree < d can be obtained due to the orthogonality of the hidden subspaces.

Preliminaries

We first recall some basics of Gröbner bases as the main tool to approach non-linear systems. Then we describe precisely our target problem and its relation with the quantum money scheme proposed in [START_REF] Aaronson | Quantum money from hidden subspaces[END_REF]. Before that we fix some general notation: we denote by F q the finite field with q elements, we set x = (x 1 , . . . , x n ) and F q [x] = F q [x 1 , . . . , x n ] to be the polynomial ring over F q in the unknowns x 1 , . . . , x n . We denote by M(F q [x]) the set of monomials in F q [x], M(F 2 [x]) refers to the set of squarefree monomials in F 2 [x] and M s (F q [x]) refers to the set of monomials of degree s in F q [x]. As usual, M k,ℓ (F q ) denotes the set of k × ℓ matrices with entries in F q , M k (F q ) denotes the square matrices of order k with entries in F q , and GL k (F q ) the set of invertible matrices in M k (F q ).

Basics of Computer Algebra

As systems of multivariate non-linear equations are the key component of this work we recall some aspects of Gröbner bases computations [START_REF] Buchberger | Comments on the translation of my phd thesis[END_REF][START_REF] Buchberger | Bruno buchberger's phd thesis 1965: An algorithm for finding the basis elements of the residue class ring of a zero dimensional polynomial ideal[END_REF][START_REF] Buchberger | Ein Algorithmus zum Auffinden der Basiselemente des Restklassenringes nach einem nulldimensionalen Polynomideal[END_REF]. Given a polynomial ideal over F q , say I

= f 1 , . . . , f s = { s i=1 f i h i | h 1 , . . . , h s ∈ F q [x
]}, Gröbner bases provide a way to obtain the variety V q (I) = {x ∈ F q | f i (x) = 0, for all 1 ≤ i ≤ s} by transforming the initial generators of the ideal into new generators with better properties, in the sense that computing the variety becomes simpler (this "better" set of generators is precisely the Gröbner basis).

The classic method to compute Gröbner bases is Buchberger's algorithm [START_REF] Buchberger | Comments on the translation of my phd thesis[END_REF][START_REF] Buchberger | Bruno buchberger's phd thesis 1965: An algorithm for finding the basis elements of the residue class ring of a zero dimensional polynomial ideal[END_REF][START_REF] Buchberger | Ein Algorithmus zum Auffinden der Basiselemente des Restklassenringes nach einem nulldimensionalen Polynomideal[END_REF], but more efficient methods, such as F 4 [START_REF] Faugère | A new efficient algorithm for computing Gröbner bases (F4)[END_REF] and F 5 [START_REF] Faugère | A new efficient algorithm for computing Gröbner bases without reduction to zero (F5)[END_REF], have been proposed. F 5 is considered to be one of the most efficient algorithms up to date for computing Gröbner bases. It uses linear algebra techniques and suppresses useless computations carried out in Buchberger's algorithm.

For increasing values of d, the F 5 algorithm successively reduces to row echelon form matrices of the form

A d =     m 1 f 1 . . . . . . . . . . . . m 2 f 2 . . . . . . . . . . . . m 3 f 3 . . . . . . . . . . . . . . . . . . . . . . . . . . .     (1)
where the columns are indexed by the monomials ordered decreasingly with respect to <, and m i are monomials such that deg(m j f i ) ≤ d. At some point, for some d, the reduced row echelon form of the matrix A d contains a Gröbner basis. This maximum degree d reached during a F 5 computation is called the degree of regularity and it is an important parameter when assessing the running time of a Gröbner basis computation. Systems verifying certain hypotheses are called semi-regular [START_REF] Bardet | On the complexity of Gröbner basis computation of semiregular overdetermined algebraic equations[END_REF][START_REF] Bardet | Asymptotic behaviour of the degree of regularity of semi-regular polynomial systems[END_REF], and they are interesting due to two reasons. Firstly, because if a system is chosen at random it turns out to be semi-regular with high probability, and secondly because the degree of regularity is known for this kind of systems. In fact, if {f 1 , . . . , f s } ⊂ F q [x] is a semi-regular system, where each f i has degree d i , its degree of regularity is given by the index of the first non-positive coefficient of the power series

k≥0 c k z k = s i=1 (1 -z d i ) (1 -z) n . (2)
The time complexity of computing a Gröbner basis [START_REF] Bardet | On the Complexity of the F5 Gr bner basis Algorithm[END_REF] is roughly given by the time spent carrying out the row echelon reduction of

A d , which is O((#A d ) ω ), where 2 ≤ ω ≤ 3 is the linear algebra constant.
Since the size of the matrix A d can be roughly approximated by O n d , this gives an overall complexity of O n ωd .

Definition of the Problem

From now on we will always assume that n is even. Recall that we are focusing on analyzing the hardness of the following problem:

Hidden Subspaces Problem (HSP q ) Input : p 1 , . . . , p m , q 1 , . . . , q m ∈ F q [x] of degree d ≥ 3 (n ≤ m ≤ 2m).
Find : a subspace A ⊂ F q n of dimension n/2 such that p i (A) = 0 ∀i ∈ {1, . . . , m} and q j (A ⊥ ) = 0 ∀j ∈ {1, . . . , m},

where A ⊥ denotes the orthogonal complement of A with respect to the standard scalar product in F q .

As mentioned in the introduction, HSP q arises in relation to the security of the non-noisy version of the quantum money scheme proposed in [START_REF] Aaronson | Quantum money from hidden subspaces[END_REF]. The private key of this scheme is a subspace A ⊂ F q n , and the polynomials p 1 , . . . , p m ∈ F q [x] (vanishing on A) and q 1 , . . . , q m ∈ F q [x] (vanishing on A ⊥ ) are the public key. To output money, the bank queries an oracle to obtain a basis of A and using this description of the subspace it generates a quantum state $ which is the banknote. The verifying process is based on the fact that it is easy to check whether a given element is a zero of a polynomial or not.

The recovery of A compromises the security of the scheme, so it becomes crucial that the HSP q cannot be easily solved. It is conjectured in [START_REF] Aaronson | Quantum money from hidden subspaces[END_REF] that, for big enough d, there is no polynomial-time algorithm that solves HSP 2 with success probability Ω 2 -n/2 .

Before proceeding any further we need to detail how the keys are generated. This is specified in [START_REF] Aaronson | Quantum money from hidden subspaces[END_REF]. The generation of an uniformly random subspace is clear by just choosing a full rank matrix in M n/2,n (F q ). The generation of an uniformly random polynomial vanishing on a given subspace can be done in O n d -time Lemma 1 ( [START_REF] Aaronson | Quantum money from hidden subspaces[END_REF]). Denote by I d,A the set of polynomials of degree d that vanish on A, by e i ∈ F q n the vector that has a 1 in its i-th position and 0 elsewhere, and by E the subspace generated by the vectors e 1 , . . . , e n/2 . We have:

1. A polynomial is in I d,E if and only if each of its monomials is divisible by an element in the set

{x n/2+1 , . . . , x n }. 2. If L is an invertible linear transformation on I d,A , the function p (x) → p (xL) maps I d,A to I d,AL -1 .
Applying this lemma, one can generate polynomials vanishing on the appropriate subspace in the following way: We have preformed all our experiments using Proposition 3.

3 The HSP q , for q > 2

We analyze the hardness of the HSP q for q > 2. This is an open problem in [START_REF] Aaronson | Quantum money from hidden subspaces[END_REF] that arises when studying the security of the quantum money scheme extended to a general field F q . We conclude that the quantum money scheme extended to F q is not secure. First we show that, with a certain probability, the HSP q can be modeled by a suitable set of non-linear equations. Then we prove that, with very high probability, enough linear equations that are linearly independent can be extracted from it. This results in a randomized polynomial-time algorithm for the HSP q (q > 2).

General Modeling of HSP q

In this part we show that the HSP q can be rather naturally modeled as a set of algebraic equations. The first straightforward modeling presented is however not optimal as it includes many equivalent solutions. We show how we can use the structure of our problem to remove the unnecessary solutions.

We abuse notation and denote by A either a subspace of F n q of dimension n/2 or a matrix in M n/2,n (F q ) whose rows are the elements of a basis of the subspace A. Proposition 4. Let p = (p 1 , . . . , p m ), q = (q 1 , . . . , q m ) ∈ F q [x] m × F q [x] m be a degree-d instance of HSP q . Let y 1 , . . . , y n/2 be variables and G (resp. G ⊥ ) be a formal matrix of size n/2×n (resp. n/2×n). We consider the system:

SysNaive HSPq = Coeff(p i , t), Coeff(q j , t) | ∀i, j ∈ {1, . . . , m}, ∀t ∈ M F q [y 1 , . . . , y n/2 ] ( 3
)
where Coeff(p i , t) denotes the coefficient of t ∈ M(F q [g 1 , . . . , g N ]) in p i (y 1 , . . . , y n/2 )•G and Coeff(q j , t)

the coefficient of t ∈ M(F q [g ⊥ 1 , . . . , g ⊥ N ]) in q j (y 1 , . . . , y n/2 ) • G ⊥ . SysNaive HSPq is a system of O(n d )
algebraic equations over F q in n 2 variables (the entries of G and G ⊥ ). Let A ⊂ F n q be a vector subspace of dimension n/2. If A is a solution of HSP q on (p, q) ∈ F q [x] m × F q [x] m then the components of A and A ⊥ vanish all the equations SysNaive HSPq .

Proof. This is an immediate consequence of the fact that every element of the subspace A (resp. A ⊥ ) can be expressed as (y 1 , . . . , y n/2 )A (resp. (y 1 , . . . , y n/2 )A ⊥ ). As a consequence, all the coefficients of the polynomials p i (y 1 , . . . , y n/2 ) • A (resp. p i (y 1 , . . . , y n/2 ) • A ⊥ must be equal to zero.

⊓ ⊔

It is easy to see that SysNaive HSPq has many solutions which are equivalent. If a vector subspace A ⊂ F n q is a solution of HSP q , then any basis of A will be a solution SysNaive HSPq . It is then natural to define a canonical form of the solutions of HSP q .

Lemma 2. Let p = (p 1 , . . . , p m ), q = (q 1 , . . . , q m ) ∈

F q [x] m × F q [x] m be a degree-d instance of HSP q . Let A ⊂ F n q be a vector subspace of dimension n/2. If A is a solution of HSP q on (p, q) ∈ F q [x] m ×F q [x] m , then for any S ∈ GL n/2 (F q ), S • A is a solution of HSP q on (p, q) ∈ F q [x] m × F q [x] m .
Proof. For all i, 1 ≤ i ≤ m, p i (y 1 , . . . , y n/2 )S • A = 0 holds as a consequence of S being invertible. Also, since (SA) ⊥ = A ⊥ (as A ⊥ (SA) T = A ⊥ A T S T = 0 considering that A ⊥ A T = 0), it also holds that q j (y 1 , . . . , y n/2 )(SA) ⊥ = 0 for all j, 1 ≤ j ≤ m.

⊓ ⊔

A direct consequence of Lemma 2 is that we can assume, with high probability, that a vector space solution A of HSP q is given in systematic form. This is, we can suppose that A = (I|G), where G ∈ GL n/2 (F q ) and I is the n/2 × n/2 identity matrix. If A has such a form then A ⊥ = (G T |I).

Fact 1 We recall that the probability that a random matrix in M n (F q ) is invertible is given by:

γ q (n) = n i=1 1 - 1 q i .
It is well known that:

lim n→∞ γ q (n) = 1 - 1 q + O 1 q 2 .
For big values of q, γ q (n/2) is close to 1, which justifies the restriction on the shape of the subspace A.

We can now improve the modeling thanks to a canonical form of the solutions. We remove all the solutions of SysNaive HSPq which correspond to equivalent bases. To do so, we generate a similar system of equations but with a smaller number of variables.

Proposition 5. Let p = (p 1 , . . . , p m ), q = (q 1 , . . . , q m ) ∈ F q [x] m × F q [x] m be a degree-d instance of HSP q . Let y 1 , . . . , y n/2 be variables and G be a formal matrix of size n/2 × n/2 and N = n 2 /4. We consider the system:

Sys HSPq = Coeff(p i , t), Coeff(q j , t) | ∀i, j ∈ {1, . . . , m}, ∀t ∈ M F q [y 1 , . . . , y n/2 ] ( 4
)
where Coeff(p i , t) denotes the coefficient of t ∈ M(F q [g 1 , . . . , g N ]) in p i (y 1 , . . . , y n/2 ) • (I|G) , and Coeff(q j , t) the coefficient of t ∈ M(F q [g 1 , . . . , g N ]) in q j (y 1 , . . . , y n/2 ) • (G T |I) . Sys HSPq is a system of O(n d ) algebraic equations over F q in N variables (the entries of G).

Let A ⊂ F n q be a vector subspace of dimension n/2. If A is a solution of HSP q on (p, q) ∈ F q [x] m × F q [x] m , then A admits with probability γ q (n/2) a basis in systematic form whose components vanish all the equations Sys HSPq .

Proof. This follows easily from Lemma 2 and Proposition 4.

⊓ ⊔

So, Lemma 2 permitted to divide by 4 the number of variables that we have to consider.

3.2 Randomized Polynomial-Time Algorithm for HSP q , with q > 2 According to Proposition 5, solving HSP q is equivalent with high probability (w.h.p) to solve the nonlinear system Sys HSPq . In this part we show that the non-linear system can be solved in polynomial-time. This is due to the fact that we can extract from Sys HSPq sufficiently many linear equations that are linearly independent.

Lemma 3. Let p = (p 1 , . . . , p m ), q = (q 1 , . . . , q m ) ∈ F q [x] m × F q [x] m be a degree-d instance of HSP q . Let p

(1) i (resp. q

(1)

i ) be the homogeneous component of degree 1 of p i (resp. q i ), that is:

p i (1) = n j=1
λ p i,j x j , where λ i,1 , . . . , λ i,n ∈ F q , q i (1) = n j=1 λ q i,j x j , with λ i,1 , . . . , λ i,n ∈ F q .

For i ∈ {1, . . . , m} and k ∈ {1, . . . , n/2}, the linear equations:

         n/2 j=1 λ p i,j+n/2 g j+n/2(k-1) + λ p i,k n/2 j=1
λ q i,j g k+n/2(j-1) + λ q i,k+n/2 are in Sys HSPq .

Proof. Let G be a formal matrix of size n/2 × n/2. If we expand the products (y 1 , . . . , y n/2 )(I|G) and (y 1 , . . . , y n/2 )(G T |I), we can see that Sys HSPq is obtained from the coefficients of:

           p i y 1 , . . . , y n/2 , n/2-1 t=0 g tn/2+1 y t+1 , . . . , n/2-1 t=0 g tn/2+n/2 y t+1 , ∀i, 1 ≤ i ≤ m, q j n/2 t=1 g t y t , . . . , n/2 t=1 g n/2(n/2-1)+t y t , y 1 , . . . , y n/2 , ∀j, 1 ≤ j ≤ m. (5)
It is clear that the equations Coeff(p i , y 1 ), . . . , Coeff(p i , y n/2 ) on the one hand, and the equations Coeff(q i , y 1 ), . . . , Coeff(q i , y n/2 ) on the other hand, are linear for all i ∈ {1, . . . , m}. Taking into account the expressions of p

(1) i and q

(1) i as well as (5), we have that for k = 1, 2, . . . , n/2:

Coeff(p i , y k ) = λ p i,k + n/2 j=1 λ p i,j+n/2 g j+(k-1)n/2 ,
Coeff(q i , y k ) = λ q i,k+n/2 + n/2 j=1 λ q i,j g k+(j-1)n/2 , as required.

⊓ ⊔

Let N = n 2 /4. Since m ≥ n, the system of linear equations in Lemma 3 is already overdetermined with at most 2mn/2 = mn ≥ 4N linear equations versus N unknowns. We show now that among these (at most) mn linear equations there are, with high probability, at least N linearly independent ones, enough to solve it. Lemma 4. Let p = (p 1 , . . . , p m ), q = (q 1 , . . . , q m ) ∈ F q [x] m × F q [x] m be a degree-d instance of HSP q . With probability γq(m) γq(m-n/2) , we can extract from Sys HSPq at least N = n 2 /4 linear equations that are linearly independent.

Proof. The mn × N matrix of coefficients associated to the linear system specified in Lemma 3, whose columns are the unknowns g 1 , . . . , g n/2 , . . . , g N -n/2 , . . . , g N , is the following:

                  . . . .
     . ( 6
)
this is, if the matrix ( 6) is of maximum rank. Since the coefficients of the matrix are uniformly random, the probability that a m × n/2 matrix has maximum rank is, according to [START_REF] Brent | Determinants and rank of random matrices over Zm[END_REF], precisely

(1 -1 q ) . . . (1 -1 q m ) (1 -1 q ) . . . (1 -1 q m-n/2 ) = γ q (m) γ q (m -n/2)
.

⊓ ⊔

Considering that the shape of A is of the restricted form we assumed with probability γ q (n/2) and that the system above can be solved successfully with probability γq(m) γq(m-n/2) , the following theorem sums up the results of this section:

Theorem 3. Let q > 2.
There is a randomized polynomial-time algorithm solving HSP q in:

O(n 2ω ),
where 2 ≤ ω ≤ 3 is the linear algebra constant, and with success probability

γ q (n/2)γ q (m) γ q (m -n/2)
.

The success probability of our algorithm can be asymptotically approximated by 1 -1/q.

Proof. The algorithm to solve HSP q is the following:

Input: p 1 , . . . , p m , q 1 , . . . , q m ∈ F q [x] of degree d ≥ 3.
Construct the linear system of Lemma 3. Solve it. Return this solution.

Taking into account that γ q (n) = n i=1

1 -1 q i and that lim

n→∞ γ q (n) = 1 -1 q + O 1 q 2 , lim n→∞ γ q (n/2)γ q (m) γ q (m -n/2) = 1 - 1 q + O 1 q 2
follows, and so the asymptotic success probability of our algorithm increases as we increase q. ⊓ ⊔

Experimental Results

We report here our experimental results for HSP q , with q > 2, obtained with the algorithm of Theorem 3. We have implemented the algorithm using the Magma software [START_REF] Bosma | The Magma algebra system I: The user language[END_REF]. In the tables below, NextPrime(k) is the Magma function that outputs the least prime number greater than k. Also, Time gen is the timed needed to generate the instances, and Time is the time spent solving the linear system. Finally, N = n 2 /4 is the number of unknowns in the linear system. As expected from Theorem 3 the algorithm is very efficient. Note that even for small q, all experiments performed succeeded as the probability of obtaining sufficiently many linearly independent linear equations, γ q (m)/λ q (m -n/2), tends to 1 very quickly even for small values of q. Note that the running time of our algorithm is clearly dominated by the time spent in generating the instance. This is done in polynomial time 5 , so we can infer from the experiments that the algorithm runs in polynomial time, which is coherent with the theoretical results obtained.

d = 3 n q N Timegen

An Efficient Algorithm for Solving HSP 2

We consider in this part the special case of HSP 2 . As in Proposition 5, we can model HSP 2 by a set of algebraic equations. However, the system for q = 2 will have a different structure than Sys HSPq . In particular, it is no longer possible to extract linear equations (this is due to the field equations). As a consequence we have to adopt a different strategy for solving HSP 2 . First we particularize the modeling for HSP 2 : Proposition 6. Let p = (p 1 , . . . , p m ), q = (q 1 , . . . , q m ) ∈ F 2 [x] m × F 2 [x] m be a degree-d instance of HSP 2 . Let y 1 , . . . , y n/2 be variables and G be a formal matrix of size n/2 × n/2 and N = n 2 /4. We consider the system:

Sys HSP 2 = Coeff(p i , t), Coeff(q j , t) | ∀i, j ∈ {1, . . . , m}, ∀t ∈ M F 2 [y 1 , . . . , y n/2 ] ( 7
)
where

Coeff(p i , t) denotes the coefficient of t ∈ M(F 2 [g 1 , . . . , g N ]) in p i (y 1 , . . . , y n/2 ) • (I|G) , and Coeff(q j , t) the coefficient of t ∈ M(F 2 [g 1 , . . . , g N ]) in q j (y 1 , . . . , y n/2 ) • (G T |I) . Sys HSP 2 is a system of at most 2m n/2 1 + n/2 2 + . . . + n/2 d algebraic equations over F 2 in N variables (the entries of G). Let A ⊂ F n 2 be a vector subspace of dimension n/2. If A is a solution of of HSP 2 on (p, q) ∈ F 2 [x] m × F 2 [x] m ,
then A admits with probability γ 2 (n/2) a basis in systematic form whose components vanish all the equations Sys HSP 2 .

Proof. Direct application of Proposition 5 and the fact that

#(M(F 2 [y 1 , . . . , y n/2 ]) = n/2 1 + . . . + n/2 d .

⊓ ⊔

Note that the equations of Sys HSP 2 are of degree d with high probability for big enough parameters. Indeed, let G be a formal matrix of size n/2×n/2 and recall that Sys HSP 2 is obtained from the coefficients of:

           p i y 1 , . . . , y n/2 , n/2-1 t=0 g tn/2+1 y t+1 , . . . , n/2-1 t=0 g tn/2+n/2 y t+1 , ∀i, 1 ≤ i ≤ m, q j n/2 t=1 g t y t , . . . , n/2 t=1 g n/2(n/2-1)+t y t , y 1 , . . . , y n/2 , ∀j, 1 ≤ j ≤ m.
Since we reduce modulo the field equations, the coefficient of a linear term is obtained from the linear terms of p i and q j but also from the coefficients of higher degree terms reduced modulo the field equations. So, expect with a high probability that Sys HSP 2 has no linear equation. However, although Sys HSP 2 is non-linear it is greatly overdetermined. Thus we can expect that computing a Gröbner basis of Sys HSP 2 can still be done efficiently.

Experimental results and Interpretation

The goal of this part is to show that Sys HSP 2 is indeed much easier to solve than a semi-regular system of the same size. Recall that if a system is semi-regular, its degree of regularity is given by the first non-positive coefficient of the power series specified in [START_REF] Bardet | On the complexity of Gröbner basis computation of semiregular overdetermined algebraic equations[END_REF].

We report experiments run on a 2.93 GHz Intel PC with 128 Gb. of RAM with the Magma software [6] (V2.19-1) for the most disadvantageous choice of parameters (this is, m = n). We recall that Magma implements the F 4 algorithm ( [START_REF] Faugère | A new efficient algorithm for computing Gröbner bases (F4)[END_REF]) for computing Gröbner basis. The notation used in the table is the following: n the number of variables of the public polynomials, N = n 2 /4 is the number of unknowns of the system in proposition 6, U eqs is the upper bound on the number of equations as specified in proposition 6, d reg is the degree of regularity observed in practice, and d sg reg is the theoretical degree of regularity treating the system as if it was semi-regular. The first thing we observe is that the number of equations of our system coincides with the upper bound for the maximum number of possible equations. The experiments show that solving these systems is easier than if they were random: the degree of regularity observed in practice is indeed lower than the expected one, which suggests that there is an underlying structure that can be exploited. Furthermore, the degree of regularity stays bounded, so we conjecture the following: Conjecture 1. The degree of regularity is bounded above by d + 1.

If our conjecture is true, there is a randomized polynomial-time algorithm for HSP 2 as follows:

Theorem 4. Assuming Conjecture 1, there is a randomized polynomial-time algorithm (the computation of a Gröbner basis) solving degree-d instances of HSP 2 with a complexity of

O(n 2ω(d+1) ),
where 2 ≤ ω ≤ 3 is the linear algebra constant, and success probability γ 2 (n/2).

Structural Low Degree Equations

The goal of this part is to provide theoretical arguments supporting Conjecture 1. That is, we show that the system of algebraic equations of Proposition 6 has a very particular structure. We prove that suitable linear combinations of the equations will lead to equations of a lower degree. This is actually the first computation performed by a Gröbner basis algorithm on the system of Proposition 6. As a consequence, solving the system of degree d-equations from Proposition 6 reduces to solve a system with equations of degree d and degree d -1. This is typically a behaviour which is not occurring in a random (i.e. semi-regular) system of equations and so it is a first step towards proving Conjecture 1.

Let (p, q) ∈ F 2 [x] × F 2 [x] be polynomials of degree d such that p vanishes on a vector subspace A ⊂ F n 2 of dimension n/2 and q vanishes on the orthogonal space A ⊥ . Let N d = n/2 d . We order lexicographically the monomials of degree d in the rings F 2 y 1 , . . . , y n/2 , F 2 x n/2+1 , . . . , x n , and F 2 x 1 , . . . , x n/2 .

We then denote by t 1 < . . . < t N d , m 1 < . . . < m N d and m ⊥ 1 < . . . < m ⊥ N d the respective monomials in ascending order. This way we can write

p = α 1 m 1 + . . . + α N d m N d + p, α 1 , . . . , α N d ∈ F 2 , p ∈ F 2 [x] \M d (F 2 [x n 2 +1 , . . . , x n ]), q = β 1 m ⊥ 1 + . . . + β N d m ⊥ N d + q, β 1 , . . . , β N d ∈ F 2 , q ∈ F 2 [x] \M d (F 2 [x 1 , . . . , x n/2 ]). ( 8
)
Recall that the notation Coeff(p, t) refers to the coefficient of t ∈ M(F 2 [g 1 , . . . , g N ]) occurring in p (y 1 , . . . , y n/2 )(I|G) , and Coeff(q, t) refers to the coefficient of t ∈ M(F 2 [g 1 , . . . , g N ]) occurring in p (y 1 , . . . , y n/2 )(G T |I) . We will denote by Coeff(p, t) (d) the homogeneous component of degree d of Coeff(p, t).

Let t ∈ M(F 2 [g 1 , . . . , g N ]), and we can deduce from (5) that Coeff(p, t) (d) = 0 = Coeff(q, t) (d) . Then, the homogeneous component of degree d of Coeff(p, t) (resp. Coeff(q, t)) is equal to the sum of the contributions with terms of degree d that each monomial m i (resp. m ⊥ j ) present in (8) makes to, respectively, Coeff(p, t) (d) and Coeff(q, t) (d) , this is, d) .

Coeff(p, t) (d) = α 1 Coeff(m 1 , t) (d) + . . . + α N d Coeff(m N d , t) (d) , Coeff(q, t) (d) = β 1 Coeff(m ⊥ 1 , t) (d) + . . . + β N d Coeff(m ⊥ N d , t) (
The fact that G and G T have the same entries (in different positions) and are involved in the evaluations of p i ((y 1 , . . . , y n/2 )(I|G)) and q j ((y 1 , . . . , y n/2 )(G T |I)) produces certain relations between expressions of the form Coeff(m i , t 1 ) (d) and expressions of the form Coeff(m ⊥ j , t 2 ) (d) for appropriate

t 1 , t 2 ∈ M d (F 2 [y 1 , . . . , y n/2 ]
). These relations are detailed in the following result:

Proposition 7. Let (p, q) ∈ F 2 [x] × F 2 [
x] be polynomials of degree d such that p vanishes on a vector subspace A ⊂ F n 2 of dimension n/2 and q vanishes on the orthogonal space A ⊥ . Let N d = n/2 d . For all i, j ∈ {1, . . . , N d }, it holds that:

Coeff m i , t j (d) = Coeff m ⊥ j , t i (d) ,
where x i k +n/2 , y j 1 y j 2 . . . y j d (d) .

Observing the system (5), this is the coefficient of y j 1 y j 2 . . . y j d in the product

g i k +(j π(k) -1)n/2 . (9)
On the other hand,

Coeff(m ⊥ j , t i ) (d) = Coeff( d k=1 x j k , y j 1 y j 2 . . . y j d ) (d)
.

Again, observing the system (5), this is the coefficient of y i 1 y i 2 . . . y i d in the product

Now ( 10) and ( 9) clearly coincide since π permutation over {1,...,d}

d k=1 g (j k -1)n/2+i π(k) = π permutation over {1,...,d} d k=1 g i π(k) + j π -1 (π(k)) -1 n/2 = π -1 permutation over {1,...,d} d k=1 g i k + j π -1 (k) -1 n/2 .

⊓ ⊔

We can use this proposition to identify the suitable linear combinations that are of degree d -1:

Theorem 5. Let the notations be as in Proposition 7. There exist i, j ∈ {1, . . . , N d } such that the equation Coeff(p, t j ) + Coeff(q, t i ) +

{k =i|α k =0} Coeff(q, t k ) + {ℓ =j|β ℓ =0} Coeff(p, t ℓ ) is of degree d -1.
Proof. Denote by i, j the smallest indexes such that α i , β j = 0. Using Proposition 7,

Coeff(m i , t j ) (d) = Coeff(m ⊥ j , t i ) (d) .
Now, for every k = i such that α k = 0 and for all ℓ = j such that β ℓ = 0, using Proposition 7 we get the following equalities: which means that we cancelled out the terms of degree d and so the required equation is of degree d -1.

Coeff(m k , t j ) (d) = Coeff(m ⊥ j , t k ) (d) , Coeff(m ⊥ ℓ , t i ) (d) = Coeff(m i , t ℓ ) (d

⊓ ⊔

This result can be used to generate low-degree equations:

Corollary 1 Let p = (p 1 , . . . , p m ), q = (q 1 , . . . , q m ) ∈ F 2 [x] m ×F 2 [x] m be a degree-d instance of HSP 2 . We can easily generate O(m 2 ) equations of degree d -1. These equations are linear combinations of the degree-d equations of Sys HSP 2 .

Proof. We apply simply Theorem 5 to each pair of polynomials (p i , q j ) ∈ F 2 [x] × F 2 [x]. From the proof of Theorem 5, it is clear that these equations are linear combinations of the equations from Sys HSP 2 . ⊓ ⊔

To conclude this part, we include below experimental results about the number of equations of degree d -1 generated as Corollary 1 which are linearly independent. In the table, we denote by #eqs pr the number of linearly independent equations obtained in practice and by #eqs th the maximum number of linearly independent equations that can be obtained, which is m 2 . We observe that the behaviour is unstable for small values of the parameters. This is partially due to the fact that if a polynomial p i (resp. q j ) does not have terms of degree d in F 2 [x n/2+1 , . . . , x n ] (resp. F 2 [x 1 , . . . , x n/2 ]), then we do not get equations of degree d -1 applying Theorem 5 to the pair (p i , q k ) for all k ∈ {1, . . . , m} (resp. from the pair (p k , q j ) for all k ∈ {1, . . . , m}). This happens with probability 1/2 ( n/2 d ) , which is not too small for low parameters. So, for small parameters it is possible that we obtain a number of equations of degree d -1 lower than m 2 . However, if this is the case there still are equations of degree d -1 (or lower) produced by the terms of degree d -1 (or lower). We see that the behavior becomes stable for big enough values of the parameters m, n obtaining as many equations of degree d -1 as possible, this is, m 2 .

Conclusions

In this paper we presented a very efficient attack for HSP q with q > 2. Since the asymptotic probability of success of this algorithm is 1 -1/q, the quantum money scheme extended to F q is badly broken for big q. We also provided some experimental and theoretical arguments that support the conjecture that HSP 2 can be solved in polynomial time.

Proposition 3 .

 3 (Vanishing polynomial) The generation of an uniformly random polynomial of degree d vanishing on a given subspace A consists of the following two steps:1. Generate a polynomial p (x) of degree d vanishing on E: by lemma 1(1), this is done including each monomial of degree d or lower independently and with probability 1/2 if it is divisible by an element in the set {x n/2+1 , . . . , x n }. 2. Transform the polynomial p (x) into one vanishing on A: considering the matrix L of change of basis (i.e., E = AL), the polynomial p (xL) vanishes on A by lemma 1(2).

g

  (j k -1)n/2+ℓ y ℓ which, again after expanding it, equals π permutation over {1,...,d} d k=1 g (j k -1)n/2+i π(k) .

d = 3 d = 4 #

 4 eqs pr # eqs th #eqs pr # eqs th m = n = 10

  . . . . . . . . . . . . . . . . . . . . . . . . . .We see that due to its particular shape it has rank N if there exists an n/2 × n/2 invertible submatrix in the following m × n/2 matrix:

								
	λ p i,n/2+1 . . . λ p i,n 0 . . . 0 λ p i,n/2+1 . . . λ p 0 . . . 0 . . . i,n . . .	0 0	. . . 0 . . . 0	   
		. . . . . . . . . 0 . . . 0 . . . . . . . . . λ q j,1 . . . 0	. . . . . . . . . . . . . . . . . . . . . 0 . . . 0 . . . λ p i,n/2+1 . . . λ p i,n . . . . . . . . . . . . . . . . . . . . . λ q j,2	        
								
		0	. . . λ q j,1	0	. . . λ q j,2 . . .	0	. . . λ q j,n/2	  
	        	. . . . . . . . . i,n/2+1 . . . λ p λ p i,n 0 . . . 0 λ p i,n/2+1 . . . λ p . . . . . . . . . . . . . . . . . . . . . 0 . . . 0 . . . 0 . . . 0 i,n . . . 0 . . . 0 . . . . . . . . . 0 . . . 0 0 . . . 0 . . . λ p i,n i,n/2+1 . . . λ p . . . . . . . . . . . . . . . . . . . . .	        	,
			  	λ p 1,n/2+1 λ p 1,n/2+2 . . . λ p 1,n λ p 2,n/2+1 λ p 2,n/2+2 . . . λ p 2,n
			 	. . . m,n/2+1 λ p λ p m,n/2+2 . . . λ p . . . . . . . . . m,n

. . . 0 . . . λ q j,n/2 . . . 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . We restrict our attention to the following mn/2 × N submatrix containing the equations Coeff(p i , y j ) only, for i ∈ {1, . . . , m} and j ∈ {1, . . . , n/2}: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  ) , andCoeff(m ⊥ ℓ , t k ) (d) = Coeff(m k , t ℓ ) (d) .Now adding up both the left-hand side and the right-hand side of all the four equalities, we obtainCoeff(m i , t j ) (d) + {k =i|α k =0} Coeff(m k , t j ) (d) +The left-hand side of this equality is the homogeneous component of degree d of Coeff(p, t j ) + Coeff(q, t i ) +{k =i|α k =0} Coeff(q, t k ) + {ℓ =j|β ℓ =0}Coeff(p, t ℓ )

						Coeff(m ⊥	ℓ , t i )	(d) +
				{ℓ =j|β ℓ =0}
	+	Coeff(m ⊥	ℓ , t k )	(d) + Coeff(m ⊥	j , t i )	(d) +	Coeff(m ⊥	j , t k )	(d) +
	{ℓ =j,k =i|α k =0,β ℓ =0}					{k =i|α k =0}
	+
	{ℓ =j|β l =0}

Coeff(m i , t ℓ ) (d) + {ℓ =j,k =i|α k =0,β ℓ =0} Coeff(m k , t ℓ ) (d) = 0.

Note that the generation of the instance is rather slow in practice, probably due to a non-optimal implementation of the Evaluate function in Magma for symbolic polynomials.


Piermarco Cannarsa

Marco Mazzola
email: marco.mazzola@imj-prg.fr

Carlo Sinestrari

GLOBAL PROPAGATION OF SINGULARITIES FOR TIME DEPENDENT HAMILTON-JACOBI EQUATIONS

Keywords: 2010 Mathematics Subject Classification. Primary: 58F15, 58F17; Secondary: 53C35 Dimension theory, Poincaré recurrences, multifractal analysis

We investigate the properties of the set of singularities of semiconcave solutions of Hamilton-Jacobi equations of the form ut(t, x) + H(∇u(t, x)) = 0, a.e. (t, x) ∈ (0, +∞)

It is well known that the singularities of such solutions propagate locally along generalized characteristics. Special generalized characteristics, satisfying an energy condition, can be constructed, under some assumptions on the structure of the Hamiltonian H. In this paper, we provide estimates of the dissipative behavior of the energy along such curves. As an application, we prove that the singularities of any viscosity solution of (1) cannot vanish in a finite time.

1. Introduction. It is commonly accepted that, in optimal control, a crucial role is played by the Hamilton-Jacobi equation

u t (t, x) + H x, ∇u(t, x) = 0 (t, x) ∈ (0, T ) × R n u(0, x) = u 0 (x) x ∈ R n (2)
where

• H : R n × R n → R is a C 2 smooth function such that (a) lim |p|→∞ inf x∈R n H(x, p) |p| = +∞ (b) D 2 p H(x, p) > 0, ∀(x, p) ∈ R n × R n • u 0 : R n → R is a Lipschitz function.
Indeed, being able to characterize the value function as the unique solution of (2) is the starting point towards a rigorous approach to dynamic programing.

The notion of viscosity solutions, introduced in the seminal papers [START_REF] Crandall | Viscosity solutions of Hamilton-Jacobi equations[END_REF] and [START_REF] Crandall | Some properties of viscosity solutions of Hamilton-Jacobi equations[END_REF], provides the right class of generalized solutions to study existence, uniqueness, and stability issues for problem [START_REF] Albano | Propagation of singularities for solutions of nonlinear first order partial differential equations[END_REF]. The reader will find an overview of the main features of this theory in [START_REF] Bardi | Optimal Control and Viscosity Solutions of Hamilton-Jacobi Equations[END_REF], for first order equations, and [START_REF] Fleming | Controlled Markov Processes and Viscosity Solutions[END_REF], for second order problems.

On the other hand, it is also well known that Hamilton-Jacobi equations have no global smooth solutions, in general, because solutions may develop singularitiesi.e., discontinuities of the gradient-in finite time due to crossing of characteristics.

Indeed, the maximal regularity one may expect for solutions of ( 2) is that, for any t > 0, u(t, •) is locally semiconcave on R n , that is, u(t, •) can be represented as the sum of a concave and a smooth function on each compact subset of R n . In fact, the notion of semiconcave solution was used in the past even to provide existence and uniqueness results for (2) before the theory of viscosity solution was developed, see [START_REF] Douglis | The continuous dependence of generalized solutions of non-linear partial differential equations upon initial data[END_REF], [START_REF] Kruzhkov | Generalized solutions of the Hamilton-Jacobi equations of the eikonal type I[END_REF], and [START_REF] Krylov | Nonlinear Elliptic and Parabolic Equations of the Second Order[END_REF]. Nowadays, semiconcavity is still an important property in the study of Hamilton-Jacobi when related to optimal control problems in euclidean spaces ( [START_REF] Cannarsa | Semiconcave Functions, Hamilton-Jacobi Equations and Optimal Control[END_REF]) and even on Riemannian manifolds ( [START_REF] Villani | Optimal Transport, Old and New[END_REF]). However, it is rather regarded as a regularizing effect of the nonlinear semigroup associated with (2)-in some sense, a sign of irreversibility in Hamilton-Jacobi dynamics.

Another evidence of irreversibility for the equation

u t (t, x) + H x, ∇u(t, x) = 0 (t, x) ∈ (0, T ) × R n (3)
is the persistence of singularities, that is, the fact that once a singularity is created, it will propagate forward in time up to +∞. Unlike the gain of semiconcavity, such a phenomenon is not well understood so far.

What is sufficiently clear to this date is the local propagating structure of the singular set of a viscosity solution u of (3): if (t 0 , x 0 ) ∈ [0, +∞) × R n is a singular point of u, then there exists a Lipschitz arc γ : [t 0 , t 0 + τ ) → R n such that (t, γ(t)) is singular for all t ∈ [t 0 , t 0 + τ ), see [START_REF] Albano | Propagation of singularities for solutions of nonlinear first order partial differential equations[END_REF], [START_REF] Yu | A simple proof of the propagation of singularities for solutions of Hamilton-Jacobi equations[END_REF], and [START_REF] Cannarsa | Singular dynamics for semiconcave functions[END_REF]. Therefore, the question we are now interested in is to provide conditions to ensure that τ = +∞. We note that, in general, this problem has a negative answer if H is allowed Lipschitz dependence in (t, x) even for n = 1, see [START_REF] Cannarsa | Semiconcave Functions, Hamilton-Jacobi Equations and Optimal Control[END_REF]Example 5.6.7].

A first, simple case where the answer to the above problem is positive is when n = 1 and H is sufficiently smooth. Indeed, the x-derivative of u turns out to be a solution of a conservation law for which the results in [START_REF] Dafermos | Generalized characteristics and the structure of solutions of hyperbolic conservation laws[END_REF] ensure the persistence of singularities. However, the one-dimensional case is very special because topological obstructions prevent singularities from disappearing after their onset.

Another result that guarantees the global propagation of singularities in any space dimension was obtained in [START_REF] Albano | Propagation of singularities for concave solutions of Hamilton-Jacobi equations[END_REF] for concave solutions of the Hamilton-Jacobi equation u t + H(∇u) = 0.

Requiring concavity-unlike semiconcavity-for u is, however, a restrictive assumption because it imposes a global constraint on solutions. A last result, which is strongly related to the above problem, concerns the distance function, d Ω , from the boundary of a bounded open subset Ω of a Riemannian manifold. Such a function is indeed the solution of a well known stationary Hamilton-Jacobi equation, that is, the eikonal equation. In [START_REF] Albano | Singular gradient flow of the distance fundtion and homotopy equivalence[END_REF], it is shown that the singular set of d Ω is invariant under the generalized gradient flow, a property which is crucial to prove that Ω has the same homotopy type as the singular set of d Ω .

In this paper, we address the above problem for solution of the Hamilton-Jacobi equation

u t (t, x) + H(∇u(t, x)) = 0, (t, x) ∈ (0, +∞) × Ω, ( 4
)
where Ω is a bounded domain in R n and H : R n → R is the quadratic form

H(p) = 1 2 Ap • p,
with A a positive definite n × n real matrix. Our main result, Theorem 4.2, ensures that the singularities of any viscosity solution of ( 4) persist for all times, or at least until the singular arc touches the boundary of Ω. More precisely, we show that, if (t 0 , x 0 ) is a singular point of u, then there exist T ∈ (0, +∞] and a Lipschitz continuous arc γ : [t 0 , t 0 + T ) → R n , starting from x 0 , with (s, γ(s)) singular for all s ∈ [t 0 , t 0 + T ) and such that lim s→t0+T γ(s) ∈ ∂Ω whenever T < +∞.

The proof of the above result relies on two main ideas that are converted in two technical results, respectively. In the first one, Lemma 3.2, we obtain, as in [START_REF] Albano | Singular gradient flow of the distance fundtion and homotopy equivalence[END_REF], a sharp semiconcavity estimate for a suitable transform of the solution u. In the second one, Theorem 4.1, we establish an inequality showing that the full Hamiltonian associated with (4), that is,

F (τ, p) = τ + H(p),
decreases along a selection of the superdifferential of u, evaluated at any point of a suitable arc. To be more specific, let (t 0 , x 0 ) ∈ Q and let t < t 0 . Then we prove that there exist T ′ > 0 and a Lipschitz continuous arc γ : [t 0 , t 0 + T ′ ) → Ω, starting from x 0 , such that min (τ,p)∈D + u(s,γ(s))

F (τ, p) ≤ t 0 - t s -t 2 min (τ0,p0)∈D + u(t0,x0) F (τ 0 , p 0 ) (5)
for every s ∈ [t 0 , t 0 + T ′ ). Such a dissipative behavior is essential to deduce persistence of singularities. Indeed, the above inequality yields that, if (t 0 , x 0 ) is a singular point of u, hence the minimum on the right-hand side is strictly negative, then the minimum on the left side must be negative too, thus forcing the point (s, γ(s)) to be singular as well for every s ∈ [t 0 , t 0 + T ′ ). Moreover, the quantitative estimate (5) allows to reproduce the above reasoning starting from t 0 + T ′ as long as the arc γ stays away from the boundary of Ω.

Although the structure of the Hamiltonian in ( 4) is quite special, we believe that our approach can be used to treat more complex classes of equations. In a future work, we will show how to adapt the above ideas to treat time dependent Hamilton-Jacobi equations on Riemannian manifolds.

Preliminaries.

Let Ω ⊂ R n be an open set. We consider the Hamilton-Jacobi equation

u t (t, x) + H(∇u(t, x)) = 0 a.e. (t, x) ∈ (0, +∞) × Ω =: Q u(t, x) = ϕ(t, x) for (t, x) ∈ ∂Q, (6)
where H : R n → R is the quadratic form

H(p) = 1 2 Ap • p,
with A positive definite, and ϕ : Q → R is Lipschitz continuous. Here we define

u t = ∂u ∂t and ∇u = ∂u ∂x 1 , . . . , ∂u ∂x n ,
whereas Du = (u t , ∇u) indicates the gradient of u whenever it exists. Let L denote the Legendre transform of H, i.e. L(q) = 1 2 A -1 q • q. We assume that the data ϕ satisfy the following compatibility condition

ϕ(t, x) -ϕ(s, y) ≤ (t -s)L x -y t -s , (7)
for all (t, x), (s, y) ∈ ∂Q such that t > s ≥ 0. Then, see [START_REF] Lions | Generalized Solutions of Hamilton-Jacobi Equations[END_REF], problem (6) has a unique viscosity solution u ∈ Lip( Q) which is given by the Hopf formula

u(t, x) = min (s, y) ∈ ∂Q s < t (t -s)L x -y t -s + ϕ(s, y) . (8)
Equivalently, see [START_REF] Cannarsa | Semiconcave Functions, Hamilton-Jacobi Equations and Optimal Control[END_REF], u given by ( 8) is the unique Lipschitz function on Q which satisfies ( 6) almost everywhere and is locally semiconcave on Q, that is, for every

convex compact K ⊂ Q there is a constant C ≥ 0 such that u(X + H) + u(X -H) -2u(X) ≤ C|h| 2 for every X, H ∈ R n+1 such that X + H, X -H ∈ K.
Since the function u given by ( 8) is in general non differentiable, we denote by Σ(u) the set of points of non-differentiability for u. For X ∈ Q we introduce the superdifferential of u at X

D + u(X) = P ∈ R n+1 : lim sup Y →X u(Y ) -u(X) -P, Y -X |Y -X| ≤ 0
and the set of the reachable gradients of u at X

D * u(X) = P ∈ R n+1 : Q \ Σ(u) ∋ X i → X, Du(X i ) → P . The directional derivative of u at X in the direction V ∈ R n+1 is defined by ∂u(X, V ) = lim h→0+ u(X + hV ) -u(X) h
and the exposed face of D + u(X) in the direction V by

D + u(X, V ) = P ∈ D + u(X) : P, V ≤ Q, V ∀ Q ∈ D + u(X) .
Since u is a locally semiconcave map, the set D + u(X) is the convex hull of D * u(X). Moreover, see [START_REF] Cannarsa | Semiconcave Functions, Hamilton-Jacobi Equations and Optimal Control[END_REF], if K is a convex compact subset of Q and C ≥ 0 is a semiconcavity constant of u on K, the superdifferential of u satisfies the following monotonicity property:

P -Q, X -Y ≤ C|X -Y | 2 (9)
for every

P ∈ D + u(X), Q ∈ D + u(Y ) and X, Y ∈ K.
The directional derivative of the locally semiconcave map u can be connected to the superdifferential and the reachable gradients in the following way:

∂u(X, V ) = min P ∈D + u(X) P, V = min P ∈D * u(X) P, V (10)
for any X ∈ Q and V ∈ R n+1 . Finally, we recall a result involving the exposed face of u. Its proof can be found in [START_REF] Cannarsa | Semiconcave Functions, Hamilton-Jacobi Equations and Optimal Control[END_REF].

Proposition 2.1. Let X ∈ Q and P ∈ R n+1 . Suppose that there are sequences

{X i } ⊂ Q \ {X} and P i ∈ D + u(X i ) satisfying X i → X, P i → P and lim i→∞ X i -X |X i -X| = V .
Then P ∈ D + u(X, V ).

3.

A sharp monotonicity estimate. We first prove a technical lemma. Here, for every X = (t, x) ∈ R n+1 , 0 < t < t and R > 0, we denote by B(X, R) the ball of centre X and radius R, while by Bt(X, R) we refer to the set

Bt(X, R) = (( t, +∞) × R n ) ∩ B(X, R). Lemma 3.1. For every X ′ = (t ′ , x ′ ) ∈ Q there exist R > 0 and 0 < t < t ′ such that Bt(X ′ , R) ⊂ Q and u(t, x) = min y∈Ω (t -t)L x -y t -t + u( t, y) (11)
for each (t, x) ∈ Bt(X ′ , R).

Proof. Let t > 0 and set Qt = ( t, +∞) × Ω. First we observe that

u(t, x) = min (s, y) ∈ ∂Qt s < t (t -s)L x -y t -s + u(s, y) (12)
for every (t, x) ∈ Qt. Indeed, let (t, x), (s, y) ∈ ∂Qt such that t > s ≥ t. Since t > t and (t, x) ∈ ∂Qt, we have that x ∈ ∂Ω. If y ∈ ∂Ω, by [START_REF] Crandall | Viscosity solutions of Hamilton-Jacobi equations[END_REF] we obtain

u(t, x) -u(s, y) = ϕ(t, x) -ϕ(s, y) ≤ (t -s)L x -y t -s .
Otherwise, let y ∈ Ω, so that s = t, and let (r, z) ∈ ∂Q such that r < s and

u(s, y) = (s -r)L y -z s -r + ϕ(r, z) . (13)
Then, by ( 13), [START_REF] Crandall | Viscosity solutions of Hamilton-Jacobi equations[END_REF] and the structure of L, we have

u(t, x) -u(s, y) = ϕ(t, x) -(s -r)L y -z s -r -ϕ(r, z) ≤ (t -r)L x -z t -r -(s -r)L y -z s -r = (t -r)L x -y t -r + 1 t -r A -1 (x -y) • (y -z) + (t -r)L y -z t -r -(s -r)L y -z s -r = (t -s)L x -y t -s - 1 2(t -r) s -r t -s A -1 (x -y) • (x -y) + 1 t -r A -1 (x -y) • (y -z) - 1 2(t -r) t -s s -r A -1 (y -z) • (y -z) = (t -s)L x -y t -s - 1 t -r L s -r t -s (x -y) - t -s s -r (y -z) ≤ (t -s)L x -y t -s .
In both cases, we obtain the compatibility condition

u(t, x) -u(s, y) ≤ (t -s)L x -y t -s .
Then, see [START_REF] Lions | Generalized Solutions of Hamilton-Jacobi Equations[END_REF], the right hand side of ( 12) is the unique viscosity solution of

v t (t, x) + H(∇v(t, x)) = 0 a.e. (t, x) ∈ Qt v(t, x) = u(t, x) for (t, x) ∈ ∂Qt .
Hence it coincides on Qt with the map u. Now fix X ′ = (t ′ , x ′ ) ∈ Q and let R ′ > 0 be such that B(X ′ , R ′ ) ⊂ Q. For any T > 0 the map u is Lipschitz continuous on [0, T ] × Ω, see for example [START_REF] Cannarsa | Semiconcave Functions, Hamilton-Jacobi Equations and Optimal Control[END_REF]. Let l > 0 be a Lipschitz constant for u on the set [0,

t ′ + R ′ ] × Ω. Let λ = min{A -1 z • z : z ∈ ∂B(0, 1)} > 0 ( 14
)
and let K > R ′ be such that

λ(K -R ′ ) 2 l(1 + K + R ′ ) ≥ 1 . (15)
If max{0, t ′ -1 2 } < t < t ′ , then for every (t, x) ∈ Bt(X ′ , R ′ ) and (s, y) ∈ ∂Qt such that t < s < t and |x ′ -y| ≥ K, by ( 14) we obtain

(t -s)L x -y t -s ≥ λ|x -y| 2 2(t -s) > λ(K -R ′ ) 2 .
On the other hand, by [START_REF] Strömberg | Propagation of singularities along broken characteristics[END_REF] we have

|u(t, x) -u(s, y)| < l(1 + K + R ′ ) ≤ λ(K -R ′ ) 2 .
Therefore, if max{0, t ′ -1 2 } < t < t ′ and (t, x) ∈ Bt(X ′ , R ′ ), the minimum in ( 12) is realized at some (s, y) ∈ B(X ′ , K) ∩ Q, such that (s, y) ∈ ∂Qt and s < t. Let

M = sup { |u(X)| : X ∈ B(X ′ , K) ∩ Q} . If 0 < R < R ′ and max{0, t ′ -1 2 } < t < t ′ satisfy (d ∂Ω (x ′ ) -R) 2 t ′ -t + R > 4M λ ,
then for every (t, x) ∈ Bt(X ′ , R), t < s < t and y ∈ ∂Ω we have

(t -s) L x -y t -s + u(s, y) = 1 2(t -s) A -1 (x -y) • (x -y) + u(s, y) ≥ 1 2(t -s) λ(d ∂Ω (x ′ ) -R) 2 + u(s, y) > M .
We just proved that if (t, x) ∈ Bt(X ′ , R), the minimum in ( 12) is realized at some point of the form ( t, y), y ∈ Ω. Therefore we obtain [START_REF] Fleming | Controlled Markov Processes and Viscosity Solutions[END_REF].

In [START_REF] Albano | Singular gradient flow of the distance fundtion and homotopy equivalence[END_REF], the invariance under the generalized gradient flow of the singular set of a solution u : R m → R of the eikonal equation is proved. The argument of the proof relies on the monotonicity estimate

u(X)P -u(Y )Q, X -Y ≤ |X -Y | 2 for every P ∈ D + u(X), Q ∈ D + u(Y ) and X, Y ∈ R m .
This property is a direct consequence of ( 9) and the global semiconcavity of the square of u with constant C = 2. Consider now the viscosity solution u of (6). Lemma 3.1 can be exploited in order to obtain some semiconcavity estimates for u. For example, fixed X ′ = (t ′ , x ′ ) ∈ Q, let R > 0 and 0 < t < t ′ be associated to X ′ as in Lemma 3.1. Using [START_REF] Fleming | Controlled Markov Processes and Viscosity Solutions[END_REF], it is possible to verify that for every x, h ∈ R n such that x-h, x+h ∈ B(x ′ , R),

u(t ′ , x + h) -u(t ′ , x -h) -2u(t ′ , x) ≤ Λ t ′ - t |h| 2 ,
where Λ = max{A -1 z • z : z ∈ ∂B(0, 1)}. This semiconcavity property and ( 9) imply pq, xy ≤ Λ|x -y| 2 t ′ -t for every p ∈ ∇ + u(t ′ , x), q ∈ ∇ + u(t ′ , y) and x, y ∈ B(x ′ , R). Yet, in order to study the propagation in time of the singularities of u, we need an estimate on the monotonicity of the superdifferential jointly in time and space. For this reason, as in [START_REF] Albano | Singular gradient flow of the distance fundtion and homotopy equivalence[END_REF], a suitable transform of the solution u is introduced: for any t > 0, the map vt : ( t, +∞) × Ω → R is defined by

vt(t, x) = (t -t) u(t, x) . (16)
In a similar way as before, given X ′ = (t ′ , x ′ ) ∈ Q, for suitable 0 < t < t ′ and R > 0 some semiconcavity properties of vt on Bt(X ′ , R) can be obtained, jointly in time and space. In the following lemma, we derive the resulting sharp monotonicity estimate for the superdifferential of vt.

Lemma 3.2. Let X ′ = (t ′ , x ′ ) ∈ Q and let R > 0 and 0 < t < t ′ be associated to X ′ as in Lemma 3.1. Then

P 1 -P 2 , X 1 -X 2 ≤ 2L(x 1 -x 2 )
for every X 1 , X 2 ∈ Bt(X ′ , R) and every

P i ∈ D + vt(X i ), i = 1, 2.
Proof. Fix X ′ = (t ′ , x ′ ) ∈ Q and let R and t be such that [START_REF] Fleming | Controlled Markov Processes and Viscosity Solutions[END_REF] holds true for all (t, x) ∈ Bt(X ′ , R). Let (t, x) ∈ Bt(X ′ , R) \ Σ(u). By [START_REF] Fleming | Controlled Markov Processes and Viscosity Solutions[END_REF] there exists y ∈ Ω such that

u(t, x) = (t -t)L x -y t -t + u( t, y) . Then vt(t, x) = L (x -y) + (t -t)u( t, y) (17)
and it is easy to prove that Dvt(t, x) coincides with the gradient of the right hand side of ( 17) at (t, x), that is

Dvt(t, x) = (u( t, y), ∇L(x -y)) . (18)
In general, when (t, x) ∈ Bt(X ′ , R), any element of the superdifferential D + vt(t, x) is the convex combination of elements of the form (18), since vt is locally semiconcave. Hence, given X i = (t i , x i ) ∈ Bt(X ′ , R) and

P i ∈ D + vt(X i ), i = 1, 2, there exist λ k i ≥ 0 and y k i ∈ Ω for i = 1, 2 and k ∈ {0, . . . , n + 1}, such that n+1 k=0 λ k i = 1, vt(X i ) = L(x i -y k i ) + (t i -t) u( t, y k i ) (19)
and

P i = n+1 k=0 λ k i (u( t, y k i ), ∇L(x i -y k i )) . (20)
By ( 11), ( 16) and ( 19), for every k 1 , k 2 ∈ {0, . . . , n + 1} we obtain

L(x 1 -y k1 1 ) + (t 1 -t) u( t, y k1 1 ) ≤ L(x 1 -y k2 2 ) + (t 1 -t) u( t, y k2 2 ) (21) and L(x 2 -y k2 2 ) + (t 2 -t) u( t, y k2 2 ) ≤ L(x 2 -y k1 1 ) + (t 2 -t) u( t, y k1 1 ) . ( 22
)
Thus, (20), ( 21) and ( 22) imply

P 1 -P 2 , X 1 -X 2 = n+1 k1,k2=0 λ k1 1 λ k2 2 (u( t, y k1 1 ) -u( t, y k2 2 ))(t 1 -t 2 ) + ∇L(x 1 -y k1 1 ) -∇L(x 2 -y k2 2 ), x 1 -x 2 = n+1 k1,k2=0 λ k1 1 λ k2 2 (t 1 -t)(u( t, y k1 1 ) -u( t, y k2 2 )) + (t 2 -t)(u( t, y k2 2 ) -u( t, y k1 1 )) + ∇L(x 1 -y k1 1 ) -∇L(x 2 -y k2 2 ), x 1 -x 2 ≤ n+1 k1,k2=0 λ k1 1 λ k2 2 L(x 1 -y k2 2 ) -L(x 1 -y k1 1 ) + L(x 2 -y k1 1 ) -L(x 2 -y k2 2 ) + ∇L(x 1 -y k1 1 ) -∇L(x 2 -y k2 2 ), x 1 -x 2 . ( 23
)
Observe that the special structure of L yield

L(x 1 -y k2 2 ) -L(x 2 -y k2 2 ) = 1 2 A -1 (x 1 -x 2 ) • (x 1 -x 2 ) + A -1 (x 1 -x 2 ) • (x 2 -y k2 2 ) = L(x 1 -x 2 ) + ∇L(x 2 -y k2 2 ), x 1 -x 2 . ( 24
)
Analogously,

L(x 2 -y k1 1 ) -L(x 1 -y k1 1 ) = L(x 1 -x 2 ) -∇L(x 1 -y k1 1 ), x 1 -x 2 . ( 25
)
Then, combining (23), ( 24) and (25), we obtain

P 1 -P 2 , X 1 -X 2 ≤ 2L(x 1 -x 2 )
concluding the proof.

Propagation of singularities.

In what follows, we denote by F the full Hamiltonian associated with [START_REF] Cannarsa | Singular dynamics for semiconcave functions[END_REF], that is, for every (τ, p) ∈ R × R n we set

F (τ, p) = τ + H(p).
Let u be given by the Hopf formula [START_REF] Crandall | Some properties of viscosity solutions of Hamilton-Jacobi equations[END_REF]. Then, see [START_REF] Cannarsa | Semiconcave Functions, Hamilton-Jacobi Equations and Optimal Control[END_REF], u satisfies (3) at any point (t, x) ∈ Q \ Σ(u). Consequently, for every X ∈ Q and any (τ, p) ∈ D * u(X), we have

F (τ, p) = 0. Since D + u(X) is the convex hull of D * u(X), the special structure of F implies that X ∈ Q is a singular point of u if and only if min (τ,p)∈D + u(X) F (τ, p) (26)
is strictly negative. Sufficient conditions are provided in the literature, see for example [START_REF] Cannarsa | Singular dynamics for semiconcave functions[END_REF] and [START_REF] Strömberg | Propagation of singularities along broken characteristics[END_REF], for the existence of generalized characteristics, whose dynamics are determined by selections of the superdifferentials of u that are "energy minimizing" in the sense of (26). To be more precise, in [START_REF] Cannarsa | Singular dynamics for semiconcave functions[END_REF] it was proved that if for any X 0 = (t 0 , x 0 ) ∈ Q there is a unique generalized characteristic starting from X 0 , then any generalized characteristic ξ : [t 0 , t 0 + T 0 ) → Q admits right derivative d ds + ξ(s) for all s ∈ [t 0 , t 0 + T 0 ), this is right-continuous and is given by d ds + ξ(s) = DF (τ (s), p(s)) , where (τ (s), p(s)) ∈ D + u(ξ(s)) is such that

F (τ (s), p(s)) = min (τ,p)∈D + u(ξ(s)) F (τ, p) . ( 27
)
When H is a quadratic form, the uniqueness of the generalized characteristics, given the initial data, is a consequence of Gronwall's Lemma. Then, we can state the following

Proposition 4.1. Let X ′ ∈ Q and R > 0 be such that B(X ′ , R) ⊂ Q. Then, there exists T R > 0 such that for every (t 0 , x 0 ) ∈ B(X ′ , R) there is a Lipschitz continuous arc ξ : [t 0 , t 0 + T R ) → Q satisfying the following properties: (i) ξ(t 0 ) = (t 0 , x 0 ); (ii) the right derivative d ds + ξ(s) does exist for all s ∈ [t 0 , t 0 + T R ); (iii) d ds + ξ(•) is right-continuous and d ds + ξ(s) = DF (τ (s), p(s)) = 1 ∇H(p(s)) , (28)
where (τ (s), p(s)) ∈ D + u(ξ(s)) satisfies ( 27). Moreover,

d ds + u(ξ(s)) = τ (s) + A p(s) • p(s) ∀ s ∈ [t 0 , t 0 + T R ) . ( 29
)
Proof. The first statement of the Proposition is a consequence of [6, Theorem 3.2, Corollary 3.4] and of Gronwall's Lemma. It remains to verify (29). First observe that by (28) for every s ∈ [t 0 , t 0 + T R ) we have

d ds + u(ξ(s)) = ∂u ξ(s), d ds + ξ(s) = ∂u ξ(s), 1 ∇H(p(s)) . (30)
Let (τ (s), p(s)) ∈ D + u(ξ(s)) satisfy (27) for every s ∈ [t 0 , t 0 + T R ). By (iii) of Proposition 4.1 we obtain

lim h↓0 ξ(s + h) = ξ(s) lim h↓0 (τ (s + h), p(s + h)) = (τ (s), p(s))
and

lim h↓0 ξ(s + h) -ξ(s) |ξ(s + h) -ξ(s)| = (1, ∇H(p(s))) |(1, ∇H(p(s)))| .
Then Proposition 2.1 yields

(τ (s), p(s)) ∈ D + u ξ(s), (1, ∇H(p(s))) |(1, ∇H(p(s)))| . Equivalently, (τ (s), p(s)
) ∈ arg min (τ,p)∈D + u(ξ(s)) τ p , 1 ∇H(p(s)) . (31)
Therefore, (30), ( 10), (31) and the special structure of H yield

d ds + u(ξ(s)) = τ (s) + p(s) • ∇H(p(s)) = τ (s) + Ap(s) • p(s) .
In [START_REF] Cannarsa | Singular dynamics for semiconcave functions[END_REF], it was shown that, given a singular point X 0 = (t 0 , x 0 ) of u, the singularity propagates locally in time following the generalized characteristic

ξ(s) = t 0 x 0 + s t0 1 ∇H(p(r)) dr . (32)
We shall provide an upper bound for the dissipation of the minimal energy (27) along this curve. Eventually, this estimate has the consequence that the singularities cannot actually extinguish in a finite time.

Theorem 4.1. Let X ′ = (t ′ , x ′ ) ∈ Q and let R > 0 and 0 < t < t ′ be associated to X ′ as in Lemma 3.1. Then, there exists T ′ > 0 such that for every (t 0 , x 0 ) ∈ Bt(X ′ , R 2 ) the Lipschitz continuous arc defined by (32) satisfies min (τ,p)∈D + u(ξ(s))

F (τ, p) ≤ t 0 - t s -t 2 min (τ0,p0)∈D + u(t0,x0) F (τ 0 , p 0 ) ( 33
)
for every s ∈ [t 0 , t 0 + T ′ ).

Theorem 4.1 implies the global propagation of the singularities.

Theorem 4.2. Let (t 0 , x 0 ) be a singular point of u. Then there exist T ∈ (0, +∞] and a Lipschitz continuous arc γ : [t 0 , t 0 + T ) → R n starting from x 0 , satisfying

(s, γ(s)) ∈ Σ(u) ∀ s ∈ [t 0 , t 0 + T )
and such that lim s→t0+T γ(s) ∈ ∂Ω whenever T < +∞.

Proof. Let ξ(•) =: (•, γ(•)) be the generalized characteristic starting from (t 0 , x 0 ). Set T = sup{r ≥ 0 :

ξ(t 0 + r) ∈ Σ(u)}. ( 34
)
By Theorem 4.1 we have that T > 0. Then either T = +∞ or 0 < T < +∞.

In the latter case we must have x ′ := lim s→t0+T γ(s) ∈ ∂Ω. Indeed, suppose by contradiction that x ′ ∈ Ω. Let R > 0 and 0 < t < t 0 + T be associated with the point (t ′ , x ′ ) := (t 0 + T, x ′ ) ∈ Q as in Lemma 3.1. Let T ′ > 0 be as provided by Theorem 4.1 and let max{ t, t ′ -T ′ } < s < t ′ be such that ξ(s) ∈ B(ξ(t ′ ), R 2 ). Hence, by Theorem 4.1 we have min

(τ,p)∈D + u(ξ(t ′ )) F (τ, p) ≤ s - t t ′ -t 2 min (τ,p)∈D + u(ξ(s))
F (τ, p) < 0.

Then, ξ(t 0 +T ) ∈ Σ(u)∩Q and Theorem 4.1 contradicts the maximality in (34).

Before proving Theorem 4.1, we provide an example showing that the estimate (33) is somehow sharp. Example. For ε > 0 consider the problem

u t (t, x) + 1 2 u 2 x (t, x) = 0 a.e. (t, x) ∈ (0, +∞) × R u(0, x) = (|x|-1) 2 2ε , (35)
The Hopf formula provides the unique viscosity solution of (35):

u ε (t, x) = 1 2 (|x| -1) 2 t + ε .
The singular set of this map is Σ(u ε ) = (0, +∞) × {0} and the curve ξ : (0, +∞) → R 2 defined by ξ(s) = (s, 0) is a generalized characteristic. If we compute the energy minimizing selection of the superdifferential of u ε along this curve, we obtain arg min{F (τ, p) : (τ, p) ∈ D + u ε (s, 0)} = -1 2(s + ε) 2 , 0 =: {(τ (s), p(s))} .

Then, given 0 < t 0 ≤ s, we have

F (τ (s), p(s)) = t 0 + ε s + ε 2 F (τ (t 0 ), p(t 0 )) .
Letting ε ↓ 0 and considering the associated functions u ε , the inequality (33) turns out to be sharp.

Proof of Theorem 4.1. Fix X ′ = (t ′ , x ′ ) ∈ Q. Let R > 0 and 0 < t < t ′ be associated to X ′ as in Lemma 3.1, and T R be associated to X ′ and R as in Proposition 4.1. If (t 0 , x 0 ) ∈ Bt X ′ , R 2 , let ξ : [t 0 , t 0 + T R ) → Q be the generalized characteristic starting from (t 0 , x 0 ) and, for any s ∈ [t 0 , t 0 + T R ], the vector (τ (s), p(s)) ∈ D + u(ξ(s)) satisfy

F (τ (s), p(s)) ≤ F (τ, p) ∀ (τ, p) ∈ D + u(ξ(s)).
Observe that, by (28), for every (t

0 , x 0 ) ∈ Bt X ′ , R 2
the curve ξ is Lipschitz continuous with constant l ξ ≤ 1 + Λl u , where Λ = max{A -1 z • z : z ∈ ∂B(0, 1)} and l u is a Lipschitz constant of u on the set 0,

t ′ + R 2 + T R × Ω. Setting T ′ = R 2l ξ , for every (t 0 , x 0 ) ∈ Bt X ′ , R 2 we have that ξ(s) ∈ Bt(X ′ , R) ∀ s ∈ [t 0 , t 0 + T ′ ) .
In order to obtain the claim of the Theorem, we need to verify that

F (τ (s), p(s)) ≤ t 0 - t s -t 2 F (τ (t 0 ), p(t 0 )) ∀ s ∈ [t 0 , t 0 + T ′ ) . Fix (t 0 , x 0 ) ∈ Bt X ′ , R 2
. Recalling (28), for all s ∈ [t 0 , t 0 + T ′ ] we can define the arc γ by

ξ(s) = t 0 x 0 + s t0 1 ∇H(p(r)) dr =: s γ(s)
.

Let vt be defined as in [START_REF] Villani | Optimal Transport, Old and New[END_REF] and fix s ∈ (t 0 , t 0 +T ′ ). Since ξ([t 0 , t 0 +T ′ ]) ⊂ Bt(X ′ , R), by Lemma 3.2, for every r ∈ [t 0 , s] and 0 < h < t 0 + T ′s, we have

P h -P, ξ(r + h) -ξ(r) ≤ 2L(γ(r + h) -γ(r)) (36)
for any

P h ∈ D + vt(ξ(r + h)), P ∈ D + vt(ξ(r)). Let (t, x) ∈ Bt(X ′ , R) \ Σ(u). By Lemma 3.2, there is y ∈ Ω such that u(t, x) = (t -t)L x -y t -t + u( t, y) .
Then one can prove that

Du(t, x) = -L x-y t-t 1 t-t ∇L(x -y) and Dvt(t, x) = u( t, y) ∇L(x -y) , so that Dvt(t, x) = u(t, x) 0 + (t -t)Du(t, x) .
This and the properties of the superdifferential of semiconcave functions yield that for every (t, x) ∈ Bt(X ′ , R)

D + vt(t, x) = u(t, x) 0 + (t -t)D + u(t, x) . (37)
By ( 36) and (37), we obtain in particular

u(ξ(r + h)) 0 +(r + h -t) τ (r + h) p(r + h) - u (
Setting

G h (r) = (r + h -t) 2 ρ h (r) + 1 2 A -1 γ h (r) • γ h (r) ,
the left hand side of (41) can be rewritten as (r + h -t) -1 d dr + G h (r). Now consider the right hand side of the inequality. Define

Here l u > 0 is a suitable Lipschitz constant for u. Summarizing, (41), (42), ( 43) and (44) yield d dr

+ G h (r) ≤ C(r + h -t) ω h (r) (45)
for some C > 0 independent on r and h. Integrating both sides of (45) on the interval [t 0 , s], we obtain

G h (s) ≤ G h (t 0 ) + C s t0
(r + h -t) ω h (r) .

Taking the limit as h → 0 and using dominated convergence, we obtain F (τ (s), p(s)) ≤ t 0 -t s -t 2 F (τ (t 0 ), p(t 0 )) .

This concludes the proof of the Theorem.

  + h)γ(r) ≤ 2L (γ(r + h)γ(r)) , that is h [u(ξ(r + h))u(ξ(r))] + h(r + h -t) [τ (r + h)τ (r)] + h 2 τ (r) +(r + h -t) p(r + h)p(r), γ(r + h)γ(r) + h p(r), γ(r + h)γ(r) ≤2L (γ(r + h)γ(r)) .Dividing by h 2 and settingτ h (r) = τ (r + h)τ (r) h , p h (r) = p(r + h)p(r) h , γ h (r) = γ(r + h)γ(r) h we obtain u(ξ(r + h))u(ξ(r)) h +(r + h -t)τ h (r) + τ (r) +(r + h -t) p h (r), γ h (r) + p(r), γ h (r) ≤ 2L (γ h (r)) .(38)Let ρ(r) = r t0 τ (z)dz, and ρ h (r) = ρ(r+h)-ρ(r) h. By Proposition 4.1 we haveu(ξ(r + h))u(ξ(rγ h (r) • γ h (r) .(40)By (38), (39) and (40) we obtain2ρ h (r)+ A -1 γ h (r) • γ h (r) + (r + h -t) d dr + ρ h (r) + 1 2 A -1 γ h (r) • γ h (r) ≤ A -1 γ h (r) • γ h (r) -1 h r+h r A p(z) • p(z)dz+A -1 γ h (r) • γ h (r)p(r), γ h (r) +ρ h (r) -d dr + ρ(r) .

  ω h (r) = sup r≤z≤r+h [|p(z)p(r)| + |τ (z)τ (r)|] .The sequence ω h (•) converges pointwise to 0 as h → 0+, since τ and p are rightcontinuous. We haveA -1 γ h (r) • γ h (r) -1 h r+h r A p(z) • p(z)dz ≤ 0 ,(42)A -1 γ h (r) • γ h (r)p(r), γ h (r) ≤ |γ h (r)| • 1 h r+h r |p(z)p(r)|dz ≤ Al u ω h (r)(43)and ρ h (r) -d dr + ρ(r) ≤ 1 h r+h r |τ (z)τ (r)|dz ≤ ω h (r) .

Acknowledgements. This work was co-funded by the European Union under the 7th Framework Programme "FP7-PEOPLE-2010-ITN", grant agreement number 264735-SADCO.

supported by NSF grant xx-xxxx.