Algorithms for Overcoming the Curse of Dimensionality for Certain Hamilton-Jacobi Equations Arising in Control Theory and Elsewhere

It is well known that time dependent Hamilton-Jacobi-Isaacs partial differential equations (HJ PDE), play an important role in analyzing continuous dynamic games and control theory problems. An important tool for such problems when they involve geometric motion is the level set method. This was first used for reachability problems. The cost of these algorithms, and, in fact, all PDE numerical approximations is exponential in the space dimension and time. In this work we propose and test methods for solving a large class of the HJ PDE relevant to optimal control problems without the use of grids or numerical approximations. Rather we use the classical Hopf formulas for solving initial value problems for HJ PDE. We have noticed that if the Hamiltonian is convex and positively homogeneous of degree one that very fast methods exist to solve the resulting optimization problem. This is very much related to fast methods for solving problems in compressive sensing, based on $\ell_1$ optimization. We seem to obtain methods which are polynomial in the dimension. Our algorithm is very fast, requires very low memory and is totally parallelizable. We can evaluate the solution and its gradient in very high dimensions at $10^{-4}$ to $10^{-8}$ seconds per evaluation on a laptop. We carefully explain how to compute numerically the optimal control from the numerical solution of the associated initial valued HJ-PDE for a class of optimal control problems. We show that our algorithms compute all the quantities we need to obtain easily the controller. The term curse of dimensionality, was coined by Richard Bellman in 1957 when considering problems in dynamic optimization.


Introduction to Hopf Formulas, HJ PDEs and Level Set Evolutions
We briefly introduce Hamilton-Jacobi equations with initial data and the Hopf formulas to represent the solution.We give some examples to show the potential of our approach, including examples to perform level set evolutions.
Given a continuous function H : R n → R bounded from below by an affine function, we consider the HJ PDE ∂ϕ ∂t (x, t) + H(∇ x ϕ(x, t)) = 0 in R n × (0, +∞), (1) where ∂ϕ ∂t and ∇ x ϕ respectively denote the partial derivative with respect to t and the gradient vector with respect to x of the function ϕ.We are also given some initial data where J : R n → R is convex.For the sake of simplicity we only consider functions ϕ and J that are finite everywhere.Results presented in this paper can be generalized for H : R n → R ∪ {+∞} and J : R n → Date: September 9, 2015 (Revised on March 10, 2016) R ∪ {+∞} under suitable assumptions.We also extend our results to an interesting class of nonconvex initial data in section 2.2.
We wish to compute the viscosity solution [12], [11] for a given x ∈ R n and t > 0.
Using numerical approximations is essentially impossible for n ≥ 4. The complexity of a finite difference equation is exponential in n because the number of grid points is also exponential in n.This has been found to be impossible, even with the use of sophisticated e.g.ENO, WENO, DG, methods [42], [31], [52].High order accuracy is no cure for this curse of dimensionality.
We propose and test a new approach, borrowing ideas from convex optimization, which arise in the ℓ 1 regularization convex optimization [24], [51] used in compressive sensing [6], [17].It has been shown experimentally that these ℓ 1 based methods converge quickly when we use Bregman and split Bregman iterative methods.These are essentially the same as Augmented Lagrangian methods [25] and Alternating Direction Method of Multipliers methods [23].These and related first order and splitting techniques have enjoyed a renaissance since they were recently used very successfully for these ℓ 1 and related problems [51], [24].One explanation for their rapid convergence is the "error forgetting" property discovered and analyzed in [43] for ℓ 1 regularization.
We will solve the initial value problem (1)-( 2) without discretization, using the Hopf formula [30] ϕ(x, t) = (J * + tH) * (x) (3) where the Fenchel-Legendre transform f * : R n → R ∪ {+∞} of a convex, proper, lower semicontinuous function f : R n → R ∪ {+∞} is defined by [18,28,44] (5) J * (v) = sup where •, • denotes the ℓ 2 (R n ) inner product.We also define for any v ∈ R n , v p = ( n i=1 |v i | p ) 1 p for 1 ≤ p < +∞ and v ∞ = max i=1,...,n |v i |.Note that since J : R n → R is convex we have that J * is 1-coercive [28, Prop.1.3.9,p. 46], i.e., lim x →+∞ J * (x) x 2 = +∞.When the gradient ∇ x ϕ(x, t) exists, then it is precisely the unique minimizer of (4); in addition, the gradient will also provide the optimal control considered in this paper (see Section 2) For instance, this holds for any x ∈ R n and t > 0 when H : R n → R is convex, J : R n → R is strictly convex, differentiable and 1-coercive.The Hopf formula only requires the convexity of J and the continuity of H, but we will often require that H in (1) is also convex.
We note that the case H = • 1 corresponds to H(∇ x ϕ(x, t)) = n i=1 ∂ ∂xi ϕ(x, t) , used, e.g. to compute the Manhattan distance to the zero level set of x → J(x), [15].This optimization is closely related to the ℓ 1 type minimization [6,17], where A is an m × n matrix with real entries, m < n, and λ > 0 arising in compressive sensing.Let Ω ⊂ R n be a closed convex set.We denote by int Ω the interior of Ω.If J(x) < 0 for x ∈ int Ω, J(x) > 0 for x ∈ (R n \ Ω) and J(x) = 0 for x ∈ (Ω \ int Ω), then the solution ϕ of ( 1)-( 2) has the property that for any t > 0 the set {x ∈ R n | ϕ(x, t) = 0} is precisely the set of points in R n \ Ω for which the Manhattan distance to Ω is equal to t > 0. This is an example of the use of the level set method [41].
Similarly, if we take H = • 2 then for any t > 0 the resulting zero level set of x → ϕ(x, t) will be the points in R n \ Ω whose Euclidean distance to Ω is equal to t.This fact will be useful later when we find the projection from a point x ∈ R n to a compact, convex set Ω.
We present here two somewhat simple but illustrative examples to show the potential power of our approach.Time results are presented in Section 4 and show that we can compute solution of some HJ-PDEs in fairly high dimensions at a rate below a millisecond per evaluation on a standard laptop. Let Then, for a given t > 0, the set {x ∈ R n | ϕ(x, t) = 0} will be precisely the set of points at Manhattan distance t outside of the ellipsoid determined by {x ∈ R n | J(x) ≤ 0}.Following [28,Prop. 1.3.1,p. 42] it is easy to see that We note that the function to be minimized decouples into scalar minimizations of the form The unique minimizer is the classical soft thresholding operator [33,21,14] defined for any component i = 1, . . ., n by Therefore, for any x ∈ R n , any t > 0 and any i = 1, . . ., n we have Here B(t) ⊆ {0, . . ., n} consists of indices for which |x i | ≤ t, and thus {0, . . ., n} \ B(t) corresponds to indices for which |x i | > t, and the zero level set moves outwards in this elegant fashion.We note that in the above case we were able to compute the solution analytically and the dimension n played no significant role.Of course this is rather a special problem, but this gives us some idea of what to expect in more complicated cases, discussed in section 3.
We will often need another shrink operator, i.e., when we solve the optimization problem with α > 0 and Its unique minimizer is given by and thus its optimal value corresponds to the Huber function (see [50] for instance) To move the unit sphere outwards with normal velocity 1, we use the following formula and, unsurprisingly, the zero level set of x → ϕ(x, t) is the set x satisfying x 2 = t + 1, for t > 0.
These two examples will be generalized below so that we can, with extreme speed, compute the signed distance, either Euclidean, Manhattan or various generalizations, to the boundary of the union of a finite collection of compact convex sets.
The remainder of this paper is organized as follows: Section 2 contains an introduction to optimal control and its connection to HJ-PDE.Section 3 gives the details of our numerical methods.Section 4 presents numerical results with some details.We draw some concluding remarks and give future plans in Section 5.The appendix links our approach to the concepts of gauge and support functions in convex analysis.

Introduction to Optimal Control
First, we give a short introduction to optimal control and its connection to HJ PDE which is given in (11).We also introduce positively homogeneous of degree one Hamiltonians and describe their relationship to optimal control problems.We explain how to recover the optimal control from the solution of the HJ-PDE.An appendix describes further connections between these Hamiltonians and gauge in convex analysis.Second, we present some extensions of our work.
2.1.Optimal control and HJ-PDE.We are largely following the discussion in [16], see also [19], about optimal control and its link with HJ PDE.We briefly present it formally and we specialize it to the cases considered in this paper.
Suppose we are given a fixed terminal time T ∈ R, an initial time t < T along with an initial x ∈ R n .We consider the Lipschitz solution x : [t, T ] → R n of the following ordinary differential equation ( 9) where f : C → R n is a given bounded Lipschitz function and C some given compact set of R n .The solution of ( 9) is affected by the measurable function β : (−∞, T ] → C which is called a control.We note A = {β : (−∞, T ] → C | β is measurable}.We consider the functional for given initial time t < T , x ∈ R n and control β K(x, t; β) = T t L(β(s))ds + J(x(T )), where x is the solution of (9).We assume that the terminal cost J : R n → R is convex.We also assume that the running cost L : R n → R ∪ {+∞} is proper, lower semicontinuous, convex, 1-coercive and dom L ⊆ C where dom L denotes the domain of L defined by dom L = {x ∈ R n | L(x) < +∞}.The minimization of K among all possible controls defines the value function v : R n × (−∞, T ] → R given for any x ∈ R n and any t < T by (10) v(x, t) = inf β∈A K(x, t; β).
The value function (10) satisfies the dynamic programming principle for any x ∈ R n , any t ≥ T and any τ ∈ (t, T ) The value function v also satisfies the following Hamilton-Jacobi-Bellman equation with terminal value Note that the control We have that ϕ is the viscosity solution of (11) where the Hamiltonian H : R n → R ∪ {+∞} is defined by ( 12) We note that the above HJ-PDE is the same as the one we consider thoughout this paper.In this paper we use the Hopf formula (3) to solve (11).We wish the Hamiltonian H : R n → R to be not only convex but also positively 1-homogeneous, i.e., for any p ∈ R n and any α > 0 We proceed as follows.Let us first introduce the characteristic function I Ω : R n → R ∪ {+∞} of the set Ω which is defined by We recall that C is a compact convex set of R n .We take f (c) = −c for any c ∈ C in (9) and Then, (12) gives the Hamiltonian H : R n → R defined by ( 14) Note that the right-hand side of ( where we recall that the Fenchel-Legendre is defined by (5).The nonnegativity of the Hamiltonian is related to the fact that C contains the origin, i.e., 0 ∈ C, and gauges.This connection is described in the appendix.Note that the controller β(t) for t ∈ (−∞, T ) in ( 9) is recovered for the solution ϕ of (11) since we have whenever ϕ(•, t) is differentiable.For any p ∈ R n such that ∇H(p) exists we also have H(p) = p, ∇H(p) .Thus we obtain that the control is given by β(t) = ∇H(∇ x ϕ(x, T − t)).
We present in Section3 our efficient algorithm that computes not only ϕ(x, t) but also ∇ x ϕ(x, t).We emphasize that we do not need to use some numerical approximations to compute the spatial gradient.In other words our algorithm computes all the quantities we need to get the optimal control without using any approximations.
It is sometimes convenient to use polar coordinates.Let us denote the (n − 1)-sphere by The set C can be described in terms of the Wulff shape [40] by the function W : The Hamiltonian H is then naturally defined via γ : S n−1 → R for any R > 0 and any θ ∈ S n−1 by ( 16) and where γ is defined by The for p = 0 0 otherwise.
In future work, we will also consider Hamiltonians defined as the supremum of linear forms such as those that arise in linear programming.We will devise very fast, low memory, totally parallelizable and apparently low time complexity methods for solving (11) with H given by ( 14) in the next section.
2.2.Some extensions and future work.In this section we show that we can solve the problem for a much more general class of Hamiltonians and initial data which arise in optimal control, including an interesting class of nonconvex initial data.
Let us first consider Hamiltonians that correspond to linear controls.Instead of ( 9), we consider the following ordinary differential equation dx ds where M is a n × n matrices with real entries and N (s) for any s ∈ (−∞, T ] is a n × m matrices with real entries.We can make a change of variables and we have dz ds (s) = e −sM N (s)β(s).
The resulting Hamiltonian now depends on t and ( 12) becomes: If C is a closed convex set, and L = I C we have We intend to program and test this in our next paper.
We now review some well known results about the types of convex initial value problems that yield to max/min-plus algebra for optimal control problems (see e.g.[22,35,1]

for instance). Suppose we have k different initial value problems
where all initial data J i : R n → R are convex, the Hamiltonian H : R n → R is convex and 1-coercive.Then, we may use the Hopf-Lax formula to get, for any x ∈ R n and any t > 0

So we can solve the initial value problem
by simply taking the pointwise minimum over the k solutions φ i (x, t), each of which has convex initial data.See Section 4 for numerical results.As an important example, suppose and each J i is a level set function for a convex compact set Ω i with nonempty interior and where the interiors of each Ω i may overlap with each other.We have J i (x) < 0 inside Ω i , J i (x) > 0 outside Ω i , and J i (x) = 0 at the boundary of Ω i .Then min i=1,...,n J i is also a level set function for the union of the Ω i .Thus we can solve complicated level set motion involving merging fronts and compute a closest point and the associated proximal points to nonconvex sets of this type.See section 3.
For completeness we add the following fact about the minimum of Hamiltonians.Let H i : R n → R, with i = 1, . . ., k, be k continuous Hamiltonians bounded from below by a common affine function.We consider where So we find the solution to by solving k different initial value problems and taking the pointwise maximum.See Section 4 for numerical results.
We end this section by showing that explicit formulas can be obtained for the terminal value x(T ) and the control β(t) for another class of running cost L. Suppose that C is a convex compact set containing the origin and take f (c) = −c for any c ∈ C. Assume also that L : R n → R ∪ {+∞} is strictly convex, differentiable when its subdifferential is nonempty and that dom L has a non-empty interior with dom L ⊆ C. Then the associated Hamiltonian H is defined by H = L * .Then, using the results of [13], we have that the (x, t) → ϕ(x, t) which solves (11) is given by the Hopf-Lax formula ϕ(x, t) = min y∈R n J(y) + tH * ( x−y t ) where the minimizer is unique and denoted by ȳ(x, t).Note that the Hopf-Lax formula corresponds to a convex optimization problem which allows us to compute ȳ(x, t).In addition, we can compute the gradient with respect to x since we have ∇ x ϕ(x, t) = ∇H * x−ȳ(x,t) t ∈ ∂J(y(x, t)) for any given x ∈ R n and t > 0.
For any t ∈ (−∞, T ) and fixed x ∈ R n the control is given by ).Note that both the control and the terminal value can be easily computed.More details about these facts will be given in a forthcoming paper.

Overcoming the Curse of Dimensionality for Convex Initial Data and Convex Homogeneous Degree One Hamiltonians -Optimal Control
We first present our approach for evaluating the solution of the HJ-PDE and its gradient using the Hopf formula [30], Moreau's identity [38] and the split Bregman algorithm [24].We note that the split Bregman algorithm can be replaced by other algorithms which converge rapidly for problems of this type.An example might be the primal-dual hybrid gradient method [53,8].Then, we show that our approach can be adapted to compute a closest point on a closed set, which is the union of disjoint closed convex sets with a non empty interior, to a given point.
We take J : R n → R convex and positively 1-homogeneous.We recall that solving (19), i.e., computing the viscosity solution, for a given x ∈ R n , t > 0 using numerical approximations, is essentially impossible, for n ≥ 4 due to the memory issue, and the complexity is exponential in n.
An evaluation of the solution at x ∈ R n and t > 0 for the examples we consider in this paper is of the order of 10 −8 to 10 −4 seconds on a standard laptop (see Section 4).The apparent time complexity seems to be polynomial in n with remarkably small constants.
We will use the Hopf formula [30]: Note that the infimum is always finite and attained (i.e, it is a minimum) since we have assumed that J is finite everywhere on R n and that H is continuous and bounded from below by an affine function.
The Hopf formula (20) requires only the continuity of H, but we will also require the Hamiltonian H be convex as well.We recall that the previous section shows how to relax this condition.
We will use the split Bregman iterative approach to solve this [24] v k+1 = arg min For simplicity we consider λ = 1 and consider v 0 = x, d 0 = x and b 0 = 0 in this paper.The algorithm still works for any positive λ and any finite values for v 0 , d 0 and b 0 .The sequence (v k ) k∈N and (d k ) k∈N are both converging to the same quantity which is a minimizer of (20).We recall that when the minimizer of ( 20) is unique then it is precisely the ∇ x ϕ(x, t); in other words, both (v k ) k∈N and (v k ) k∈N converge to ∇ x ϕ(x, t) under this uniqueness assumption.If the minimizer is not unique then the sequences (v k ) k∈N and (v k ) k∈N converge to an element of the subdifferential of ∂(y → f (y, t))(x) (see below for the definition of a subdifferential).We need to solve (21) and (22).Note that up to some changes of variables, both optimization problems can be reformulated as finding the unique minimizer of ( 24) where z ∈ R n , α > 0, and f : R n → R ∪ {+∞} is a convex, proper, lower semi-continuous function.Its unique minimizer w satisfies the optimal condition where ∂f (x) denotes the subdifferential (see for instance [27, p. 241], [44,Section 23]) of f at x ∈ R n and is defined by where (I + ∂f ) −1 denotes the "resolvent"operator of f (see [2, Def. 2, chp.3, p. 144], [5, p. 54] for instance).It is also called the proximal map of f following the seminal paper of Moreau [38] (see also [28,Def. 4.1.2,p. 318], [44, p. 339]).This mapping has been extensively studied in the context of optimization (see for instance [10,29,45,48]).
Closed form formulas exist for the proximal of map for some specific cases.For instance, we have seen in the introduction that (I + α ∂ • i ) −1 = shrink i (•, α) for i = 1, 2, where we recall that shrink 1 and shrink 2 are defined by ( 6) and (7), respectively.Another classical example consists of considering a quadratic form 1  2 , with A a symmetric positive definite matrix with real values, which yields , where I n denotes the identity matrix of size n.
Assume f is twice differentiable with a bounded Hessian, then the proximal map can be efficiently computed using Newton's method.Algorithms based on Newton's method require us to solve a linear system that involves an n × n matrix.Note that typical high dimensions for optimal control problems are about n = 10.For computational purposes, these order of values for n are small.
We describe an efficient algorithm to compute the proximal map of • ∞ in Section 4.2 using parametric programming [46,Chap. 11,Section 11.M].An algorithm to compute the proximal map for 1 2 • 2 1 is described in Section 4.4.
The proximal maps for f and f * satisfy the celebrated Moreau identity [38] (see also [44,Thm. 31.5,p. 338]) which reads as follows: for any w ∈ R n and any α > 0 This shows that (I + α ∂f ) −1 (w) can be easily computed from I + 1 α ∂f * −1 w α .In other words, depending on the nature and properties of the mappings f and f * , we choose the one for which the proximal point is "easier"to compute.Section 4.4 describes an algorithm to compute the proximal map of 1 2 • 2 ∞ using only evaluations of −1 using Moreau's identity (25).
We shall see that Moreau's identity ( 25) can be very useful to compute the proximal maps of convex and positively 1-homogeneous functions.
We consider problem ( 22) that corresponds to compute the proximal of a convex positively 1-homogeneous function H : R n → R. (we use H instead f to emphasize that we are considering positively 1-homogeneous functions and we set α = 1 to alleviate notations.)We have that H * is the characteristic function of a closed convex set C ⊆ R n , .i.e., the Wulff shape associated to H [40], In other words, (I + ∂H * ) −1 (z) corresponds to the projection of z on the closed convex set C that we denote by π C (z), that is for any Thus, using the Moreau identity (25), we see that (I + ∂H) −1 can be computed from the projection on its associated Wulff shape and we have for any In other words, computing the proximal map of H can be performed by computing the projection on its associated Wulff shape C.This formula is not new, see e.g.[7,39].Let us consider an example.Consider Hamiltonians of the form where A is a symmetric positive matrix.Here the Wulff shape is the ellipsoid We describe in Section 4.3 an efficient algorithm for computing the projection on an ellipsoid.Thus, this allows us to compute efficiently the proximal map of norms of the form • A using (26).
3.2.Projection on closed convex set with the Level Set Method.We now describe an algorithm based on the level set method [41] to compute the projection π Ω on a compact convex set Ω ⊂ R n with a nonempty interior.This problem appears to be of great interest for its own sake.
Let ψ : R n × [0, +∞) be the viscosity solution of the eikonal equation ( 27) where we recall that ∂ψ ∂t (y, t) and ∇ x ψ(y, t) respectively denote the partial derivatives of φ with respect to the time and space variable at (y, s), and where L : R n → R satisfies for any y ∈ R n (28) where int Ω denotes the interior of Ω.Given s > 0, we consider the set which corresponds to all points that are at a (Euclidean) distance s from Γ(0).Moreover, for a given point y ∈ Γ(s), the closest point to y on Γ(0) is exactly the projection π Ω (y) of y on Ω and we have (29) π In this paper we will assume that L is strictly convex, 1-coercive and differentiable so that ∇ x ψ(s, y) exists for any y ∈ R n and s > 0. We note that if Ω is the finite union of sets of this type then ∇ x ψ(s, y) may have isolated jumps.This presents no serious difficulties.Note that ( 27) takes the form of ( 19) with H = • 2 and J = L.We again use split Bregman to solve the optimization given by the Hopf formula (3).To avoid confusion we respectively replace J, v, d and b, by L, w, e and c in ( 21)-( 23) An important observation here is that e k+1 can be solved explicitly in (31) using the shrink 2 operator defined by (7).Note that the algorithm given by ( 21)-( 23) allows us to evaluate not only ψ(y, s) but also ∇ x ψ(y, s).Indeed for any s > 0 ∇ x ψ(y, s) = arg min v∈R n {L * (v) + sH(v) − y, v }, the minimizer being unique.Thus, the above algorithm (30)- (32) generates sequences (w k ) k∈N and (e k ) k∈N that both converge to ∇ x ψ(z, s).
From ( 27) we have ∂ψ ∂t (z, s) = − ∇ x ψ(z, s) 2 for any s > 0. We can thus compute (34).It remains to choose the initial data L related to the set Ω. We would like it to be smooth so that the proximal point in (30) can be computed efficiently using Newton's method.(If L lacks differentiability the approach can be easily modified using [46,Chap 11,Section 11.M] or using [9].)We consider Ω as Wulff shapes that are expressed, thanks to (15), with the function W : As a simple example, if W ≡ 1, then we might try L = y 2 − 1, the signed distance to Γ(0).This does not suit our purposes, because its Hessian is singular and its dual is the indicator function of the l 2 (R n ) unit ball.Instead we take 2 + 1 , both of these are convex and C 2 functions with non vanishing gradients away from the origin, i.e. near Γ(s).This gives us a hint as how to proceed.
Recall that we need to get initial data which behaves as a level set function should; i.e., as defined by (28).We also want either L or L * to be smooth enough, actually twice differentiable with Lipschitz continuous Hessian, so that Newton's method can be used.
We might take where R ≥ 0, θ ∈ S n−1 and m a positive integer.If we consider the important case of W (θ) = θ , with A a symmetric positive definite matrix, which corresponds to Ω = {x ∈ R n | x, A x ≤ 1}, then m ≥ 2 will have a smooth enough Hessian for L. In fact, m = 2 will lead us to a linear Hessian.There will be situations where using L * is preferable because L cannot be made smooth enough this way.
We can obtain L * using polar coordinates and taking 1 ≥ m > 1 2 , that is for any R ≥ 0 and any θ ∈ S n−1 where we recall that γ is defined by (17).As we shall see, L * is often preferable to L yielding a smooth Hessian for 1 ≥ m > 1 2 with m close to 1 2 .Let us consider some important examples.If the set Ω is defined by For case (a) we need only take If m ≥ 1 it is easy to see that Hessian of L is continuous and bounded for p ≥ 2. So we can use Newton's method for the first choice above.For case (b), we construct the Fenchel-Legendre of the function It is easy to see that where 1 p + 1 q = 1 and if 1 2 < m ≤ 1 the Hessian of L * is continuous.Other interesting examples include the following regions defined by For Ω to be convex we require A to be a positive definite symmetric matrix with real entries and its maximal eigenvalue is bounded by twice the minimal eigenvalue.We can take ( 35) for m ≥ 2 and see that this has a smooth Hessian.

Numerical Results
We will consider the following Hamiltonians A• with A symmetric positive definite matrix, and the following initial data •, A• with A a positive definite diagonal matrix.It will be useful to consider the spectral decomposition of A, i.e., A = P DP † where D is a diagonal matrix, P is an orthogonal matrix and P † denotes the transpose of P .The identity matrix in R n is denoted by I n .
First, we present the algorithms to compute the proximal points for the above Hamiltonians and initial data.We shall describe these algorithms using the following generic formulation for the proximal map Second, we present the time results on a standard laptop.Some time results are also provided for a 16 cores computer which shows that our approach scales very well.The latter is due to our very low memory requirement.Finally some plots that represent the solution of some HJ PDE are presented.
4.1.Some explicit formulas for simple specific cases.We are able to obtain explicit formulas for the proximal map for some specific cases of interest in this paper.For instance, as we have seen, considering f = • 1 gives for any i = 1, . . ., n where sign(β) = 1 if β ≥ 0 and −1 otherwise.The case f = • 2 yields a similar formula with µ = max( w 2 − α, 0).The two above cases are computed in linear time with respect to the dimension n.
Proximal maps for positive definite quadratic forms, i.e., f (w) = 1 2 w, Aw are also easy to compute since for any α > 0 and any z ∈ R n where we recall that A = P DP † with D a diagonal matrix and P an orthogonal matrix.The time complexity is dominated by the evaluation of the matrix-vector product involving P and P † .

The case of •
We use Moreau's identity (26) to compute (I + α • ∞ ) −1 ; that is for any α > 0 and for any z ∈ R n where we recall that π αC denotes the projection operator onto the closed convex set αC.We use a simple variation of parametric approaches that are well-known in graph-based optimization algorithm (see [46,Chap. 11,Section 11.M] for instance).
Let us assume that z / ∈ (αC).The projection corresponds to solve Now we use Lagrange duality (see [27, chap.VII] for instance).The Lagrange dual function g : [0, +∞) → R is defined by that is g(µ) = shrink 1 (z, µ) 1 − µα.We denote by μ the value that realizes the maximum of g.Then, we obtain for any z / ∈ αC that (36) where μ satisfies Computing the projection is thus reduced to computing the optimal value of the Lagrange multiplier μ.Consider the function h : [0, z 1 ] → R defined by h(µ) = shrink 1 (z, µ) 1 .We have that h is continuous, piecewise affine, h(0) = z 1 , h( z 1 ) = 0 and h is decreasing (recall that we assume z / ∈ (αC)).Following [46,Chap. 11,Section 11.M], we call breakpoints the values for which h is not differentiable.The set of breakpoints for h is B = {0} i=1,...,n {|z i |} .We sort all breakpoints in increasing order and we denote this sequence by (l i , . . ., l m ) ∈ B m with l i < l i+1 for i = 1, . . ., (m − 1), where m ≤ n is the number of breakpoints.This operation takes O(n log n).Then, using a bitonic search, we can find a simple interpolation computed in constant time yields μ that satisfies (37).We then use (36) to compute the projection.The overall time complexity is, therefore, O(n log n).

The case •
A and projection on a ellipsoid.We follow the same approach as for • ∞ .We consider Using Moreau's identity (26) we only need to compute the projection π EA (w) of w ∈ R n on the ellipsoid E A .Note that we have π EA (w) = P π ED (P † w) where we recall that A = P DP † with D and P a diagonal and orthogonal matrix, respectively.Thus we only describe the algorithm for the projection on an ellipsoid involving positive definite diagonal matrices.
To simplify notation we take d i = D ii for i = 1, . . ., n.We consider the ellipsoid E D defined by Let w / ∈ E D .We can easily show (see [26,Exercise III.8] for instance) that Π ED (w) satisfies for any i = 1. . . ., n where the Lagrange multiplier μ > 0 is the unique solution of n i=1 µ using Newton's method which generates a sequence (µ k ) k∈N converging to μ.We set the initial value to µ 0 = 0 and we stop Newton's iterations for the first k which statisfies |µ k+1 − µ k | ≤ 10 −8 .Once we have the value for μ we use (38) to obtain the approximate projection.

The cases
We have for any α > 0 and any Thus, assuming there exists β ≥ 0 such that we have for any α > 0 and any z ∈ R n (40) The existence of β and an algorithm to compute it follow.Let us assume that z = 0 (otherwise, β = 0 works and the solution is of course 0).Then, consider the function g : [0, z 1 ] → R defined by It is continuous, and g(0) = α z 1 while g ( z 1 ) = − z 1 .The intermediate value theorem tells us that there exists β such that g( β) = 0, that is, satisfying (39).The function g is decreasing, piecewise affine and the breakpoints of g (i.e, the points where g is not differentiable) are B = {0} ∪ i=1,...,n {|z i |}.We now proceed similarly as for the case • ∞ .We note (l i , . . ., l m ) ∈ B m the breakpoints sorted in increasing order, i.e., such that l i < l i+1 for i = 1, . . ., (m − 1), where m ≤ n is the number of breakpoints.We use a bitonic search to find the two consecutive breakpoints l i and l i+1 , such that g(l i ) ≥ 0 > g(l i+1 ).Since g is affine on [l i , l i+1 ] a simple interpolation yields the value β.We then compute I + α 2 ∂ • 2 1 −1 (z) using (40).
We now consider the case 1 2 • 2 ∞ .We have (for instance [18,Prop. 4 Then Moreau's identity (25) yields for any α > 0 and for any z ∈ R n which can be easily computed using the above algorithm for evaluating 4.5.Time results and illustrations.We now give numerical results for several Hamiltonians and initial data.We present time results on a standard laptop using a single core which show that our approach allows us to evaluate very rapidly HJ PDE solutions.We also present some time results on a 16 cores computer to show that are approach scales very well.We also present some plots that depict the solution of some HJ PDEs.
We recall that we consider the following Hamiltonians •, A• with A symmetric positive definite matrix, and the following initial data • with D a positive definite diagonal matrix, where the matrix D and A are defined follows: D is a diagonal matrix of size n × n defined by D ii = 1 + i−1 n−1 for i = 1, . . ., n.The symmetric positive definite matrix A of size n × n is defined by A ii = 2 for i = 1, . . ., n and A ij = 1 for i, j = 1, . . ., n with i = j.
All computations are performed using IEEE double precision floating-points where denormalized number mode has been disabled.The quanties (x, t) are drawn uniformly in [−10, 10] n × [0, 10].We present the average time to evaluate a solution for 1, 000, 000 runs.
We set λ = 1 in the split Bregman algorithm ( 21)- (23).We stop the iterations when the following stopping criteria is met: We first carry out the numerical experiments on an Intel Laptop Core i5-5300U running at 2.3 GHz.The implementation here is single threaded, i.e., only one core is used.Tables 1, 2, 3 and 4 present time results for several dimensions n = 4, 8, 12, 16 and with initial data •, D respectively.We see that it takes about 10 −8 to 10 −4 seconds per evaluation of the solution.We now consider experiments that are carried out on a computer with 2 Intel Xeon E5-2690 processors running at 2.90GHz.Each processor has 8 cores.Table 5 present the average time to compute the solution with Hamiltonian H = • ∞ and initial data J = • 2 1 for several dimensions and various number of used cores.We see that our approach scales very well.This is due to the fact that our algorithm requires little memory which easily fits in the L1 cache of each processor.Therefore cores are not competing for resources.This suggests that our approach is suitable for low-energy embedded systems.
We now consider solutions of HJ PDEs in dimension n = 8 on a 2-dimensional grid.We evaluate φ(x 1 , x 2 , 0, 0, 0, 0, 0, 0) with x i ∈ ∪ k=0,...,99 {−20 + k 40 99 } for i = 1, 2. Figures 1, 2 and 3 depict the solutions with initial data , and with Hamiltonians H = • 2 , H = • 1 .and H = •, D• for various times, respectively.Figure 4 and Figure 5 illustrate the max/min-plus algebra results described in Section 2.2. Figure 4 depicts the HJ solution for the initial data  5.Time results in seconds for the average time per call for evaluting the solution of the HJ-PDE with the initial data J = 1 2 • 2 1 , and the Hamiltonian H = • ∞ , for various dimensions, and several cores.

Conclusion
We have designed algorithms which enable us to solve certain Hamilton-Jacobi equations very rapidly.Our algorithms not only evaluate the solution but also computes the gradient of the solution.These include equations arising in control theory leading to Hamiltonians which are convex and positively homogeneous of degree 1.We were motivated by ideas coming from compressed sensing; we borrowed algorithms devised to solve ℓ 1 regularized problems which are known to rapidly converge.We apparently extended this fast convergence to include convex positively 1-homogeneous regularized problems.
There are no grids involved.Instead of complexity which is exponential in the dimension of the problems, which is typical of grid based methods, ours appears to be polynomial in the dimension with very small constants.We can evaluate the solution on a laptop at about 10 −4 − 10 −8 seconds per evaluation for fairly high dimensions.Our algorithm requires very low memory and is totally parallelizable which suggest that it is suitable for low energy embedded systems.We have chosen to restrict the presentation of the numerical experiments to norm-based Hamiltonians and we emphasize that our approach naturally extends to more elaborate positively 1-homogeneous Hamiltonians (using the min/max algebra results as we did for instance).
As an important step in this procedure we have also derived an equally fast method to find a closest point lying on Ω, a finite union of compact convex sets Ω i , such that Ω = ∪ k i Ω i has a nonempty interior, to a given point.
We can also solve certain so called fast marching [49] and fast sweeping [47] problems equally rapidly in high dimensions.If we wish to find ψ : R n → R with, say ψ = 0 on the boundary of a set Ω defined above, satisfying and locate the zero level set of u(•, t) = 0 for any given t > 0. Indeed any Of course the same approach could be used for any convex, positively 1-homogeneous Hamiltonian H (instead of • 2 ), e.g., H = • 1 .This will give us results related to computing the Manhattan distance.
We expect to extend our work as follows: (1) We will do experiments involving linear controls, allowing x and t > 0 dependence while the Hamiltonian (p, x, t) → H(p, x, t) is still convex and positively 1-homogeneous in p.The procedure was described in section 2. (2) We will extend our fast computation of the projection in several ways.We will consider in detail the case of polyhedral regions defined by the intersection of sets  problems.We expect to develop alternate approaches to several issues arising in LP, including rapidly finding the existence and location of a feasible point.
(3) We will consider nonconvex but positively 1-homogeneous Hamiltonians.These arise in: differential games as well as in the problem of finding a closest point on the boundary of a given compact convex set Ω, to an arbitrary point in the interior of Ω. (4) As an example of a nonconvex Hamiltonians we consider the following problems arising in differential games [20].We need to solve the following scalar problem for any z ∈ R n and any α > 0 min  It is easy to see that the minimizer is the stretch 1 operator which we define for any i = 1, . . ., n as: (41) (stretch 1 (x, α) We note that the discontinuity in the minimizer will lead to a jump in the derivatives (x, t) → ∂ϕ ∂xi (x, t), which is no surprise, given that this interface associated with the equation   The zero level set disappears when t ≥ max i a i as it should.
For completeness, we also consider the nonconvex optimization problem Its minimizer is given by the stretch 2 operator formally defined by stretch 2 (x, α) = x + α x x 2 if x = 0, αθ with θ 2 = 1 if x = 0.
, Paper accepted in Research in the Mathematical Sciences.Research supported by ONR grants N000141410683, N000141210838 and DOE grant DE-SC00183838.

3. 1 .
Numerical optimization algorithm.We present the steps needed to solve(19) and H corresponds to the support function C, that is for any p ∈ R n H(p) = sup s∈C s, p .Following Moreau [38, Example 3.d], the proximal point of z ∈ R n relative to H * is (I + ∂H * )
t) = 0, and the previous initial data, will move inwards, and characteristics will intersect.The solution ϕ(x, t) will remain locally Lipschitz continuous, even though a point inside the ellipsoid may be equally close to two points on the boundary of the original ellipsoid in the Manhattan metric.So we

Table 1 .
1, 1) † , and H = • 1 for various times.Figure5depicts the HJ solution for various time with J = 1 • 2 2 and H = min Time results in seconds for the average time per call for evaluting the solution of the HJ-PDE with the initial data J = 1 2 • 2 2 , several Hamiltonians and various dimensions n.

Table 2 .
Time results in seconds for the average time per call for evaluting the solution of the HJ-PDE with the initial data J = 1 2 • 2 ∞ , several Hamiltonians and various dimensions n.

Table 3 .
Time results in seconds for the average time per call for evaluting the solution of the HJ-PDE with the initial data J = 1

Table 4 .
Time results in seconds for the average time per call for evaluting the solution of the HJ-PDE with the initial data J = 1 2 •, D −1 • , several Hamiltonians and various dimensions n.