Louis Montaut

Quentin Le Lidec

Vladimir Petrik

Josef Sivic

Justin Carpentier

Collision Detection Accelerated: An Optimization Perspective

Collision detection between two convex shapes is an essential feature of any physics engine or robot motion planner. It has often been tackled as a computational geometry problem, with the Gilbert, Johnson and Keerthi (GJK) algorithm being the most common approach today. In this work we leverage the fact that collision detection is fundamentally a convex optimization problem. In particular, we establish that the GJK algorithm is a specific sub-case of the well-established Frank-Wolfe (FW) algorithm in convex optimization. We introduce a new collision detection algorithm by adapting recent works linking Nesterov acceleration and Frank-Wolfe methods. We benchmark the proposed accelerated collision detection method on two datasets composed of strictly convex and non-strictly convex shapes. Our results show that our approach significantly reduces the number of iterations to solve collision detection problems compared to the state-of-the-art GJK algorithm, leading to up to two times faster computation times.

I. INTRODUCTION

Physics engines designed to simulate rigid bodies are an essential tool used in a wide variety of applications, notably in robotics, video games, and computer graphics [START_REF] Coumans | PyBullet, a Python module for physics simulation for games, robotics and machine learning[END_REF][START_REF]Persistent Contact Manifold[END_REF][START_REF] Todorov | Mujoco: A physics engine for model-based control[END_REF]. Collision detection, a crucial feature of any physics engine or robot motion planer [START_REF] Mirabel | HPP: A new software for constrained motion planning[END_REF][START_REF] Marc A Toussaint | Differentiable physics and stable modes for tool-use and manipulation planning[END_REF][START_REF] Schulman | Motion planning with sequential convex optimization and convex collision checking[END_REF], consists of finding which objects are colliding or not, i.e. are sharing at least one common point or if there exists a separating hyper-plane between both. As simulation often needs to deal with multiple objects and run in real-time (i.e., in video games) or at very high-frequencies (i.e., in robotics), it is crucial for collision detection to be carried out as fast as possible. In order to reduce computational times, collision detection is usually decomposed into two phases thoroughly covered in [START_REF] Ericson | Real-Time Collision Detection[END_REF]. The first phase is the broad phase which consists in identifying which pair of simulated objects are potentially colliding. Objects in the pair selected during the broad phase are close to each other, hence the collision is uncertain and needs to be considered carefully. The second phase is the narrow phase in which each pair identified in the broad phase is tested to check whether a collision is indeed occurring. Collision detection during the narrow phase is the focus of this paper. Problem formulation. We consider two convex shapes A 1 and A 2 in R n (with n = 2 or 3 in common applications). If the shapes are not convex, we use their respective Since A1 and A2 are convex, D is also convex. In (a), the shapes are not in collision hence the origin of the configuration space C, 0C (in red) lies outside the Minkowski difference, 0C ̸ ∈ D. The vector x * = x * 1 -x * 2 separates A1 from A2. It is also equal to the projection of 0C onto the Minkowski difference D, x * = proj D (0C). In (b), the shapes are overlapping, thus 0C ∈ D. In this case, we have x * = proj D (0C) = 0C. convex-hulls or decompose the shapes into a collection of convex sub-shapes [START_REF] Mamou | A Simple and Efficient Approach for 3D Mesh Approximate Convex Decomposition[END_REF]. The separation distance between A 1 and A 2 , denoted by dist(A 1 , A 2 ) ∈ R + , can be formulated as a minimization problem of the form:

d 1,2 = min x1∈A1,x2∈A2 ∥x 1 -x 2 ∥ 2 and dist(A 1 , A 2 ) = d 1,2 , (1)
where x 1 ∈ A 1 and x 2 ∈ A 2 are both vectors in R n , d 1,2 is the optimal value of (1) and ∥•∥ is the Euclidian norm of R n . If A 1 and A 2 are intersecting (i.e. in collision), we have dist(A 1 , A 2 ) = 0. If the two shapes do not intersect, we have dist(A 1 , A 2 ) > 0. These two cases are illustrated in Fig. 1.

Problem [START_REF] Bach | Learning with Submodular Functions: A Convex Optimization Perspective[END_REF] allows us to encapsulate and tackle both the distance computation problem and the computationally cheaper Boolean collision check into one single convex optimization problem. In the distance computation problem, we aim at computing the separation distance between A 1 and A 2 , denoted dist(A 1 , A 2 ), i.e. the distance between their closest points. This distance is useful in some applications such as collision-free path planning [START_REF] Gilbert | Distance Functions and Their Application to Robot Path Planning in the Presence of Obstacles[END_REF][START_REF] Stasse | Real-Time (Self)-Collision Avoidance Task on a Hrp-2 Humanoid Robot[END_REF], especially for pairs of objects entering the narrow phase. If a pair of objects has not been selected by the broad phase, a cheap estimate of dist(A 1 , A 2 ) is usually enough [START_REF] Ericson | Real-Time Collision Detection[END_REF]. In the Boolean collision check we only aim at determining if A 1 and A 2 intersect, and computing dist(A 1 , A 2 ) is unnecessary. However, we will later see that the Boolean collision check is a sub-problem of the distance computation problem: solving (1) can be early-stopped once a separating plane between A 1 and A 2 has been found. In the rest of this paper, we will use the generic term "collision detection" when we mean to encapsulate both distance computation and Boolean collision checking. We will specify when the distinction is needed.

Related work. The most well-known algorithm for collision detection between two convex shapes is the so-called Gilbert-Johnson-Keerthi algorithm (GJK) [START_REF] Gilbert | A fast procedure for computing the distance between complex objects in three-dimensional space[END_REF]. It can handle both the distance computation and the Boolean collision check [START_REF] Van Den Bergen | A Fast and Robust GJK Implementation for Collision Detection of Convex Objects[END_REF]. Most alternatives to GJK found in the literature focus on computing collisions between convex polyhedra such as the Lin-Canny algorithm [START_REF] Lin | A Fast Algorithm for Incremental Distance Calculation[END_REF] or the V-Clip [START_REF] Mirtich | V-Clip: Fast and Robust Polyhedral Collision Detection[END_REF] algorithm. Although GJK is equivalent in performance to these algorithms [START_REF] Cameron | A Comparison of Two Fast Algorithms for Computing the Distance between Convex Polyhedra[END_REF], it is not restricted to convex polyhedra. The strength of GJK is formulating the collision detection problem on the Minkowski difference. The properties of the Minkowski difference are used to cleverly compute support vectors on the Minkowski difference (these notions will be introduced in detail in Sec. II). GJK is thus able to handle collision detection and distance computation for many different shapes such as convex polyhedra and basic primitives (i.e., spheres, ellipsoids, cylinders, capsules etc.) [START_REF] Van Den Bergen | A Fast and Robust GJK Implementation for Collision Detection of Convex Objects[END_REF][START_REF] Van Den Bergen | Collision Detection in Interactive 3D Environments[END_REF][START_REF] Ericson | Real-Time Collision Detection[END_REF]. The generality of GJK as well as its efficiency, good precision and ease of implementation makes it the state-of-the-art algorithm for collision detection between two convex shapes. The computation time of GJK has been reduced thanks to improvements to its sub-operations [START_REF] Montanari | Improving the GJK Algorithm for Faster and More Reliable Distance Queries Between Convex Objects[END_REF][START_REF] Van Den Bergen | A Fast and Robust GJK Implementation for Collision Detection of Convex Objects[END_REF]. However, to the best of our knowledge, no previous work has focused explicitly on lowering the number of iterations of GJK. This work notably shows how to accelerate collision detection by directly lowering the number of iterations needed to solve a collision problem instance compared to the vanilla GJK algorithm.

Over the years, collision detection has often been seen as a computational geometry problem. However we argue that this view has limited the improvement of collision detection and propose to reframe it as a convex optimization problem. As briefly mentioned already in their 1988 paper [START_REF] Gilbert | A fast procedure for computing the distance between complex objects in three-dimensional space[END_REF] and brought-up again by [START_REF] Qin | Smoothing Algorithms for Computing the Projection onto a Minkowski Sum of Convex Sets[END_REF], the ideas developed by Gilbert, Johnson and Keerthi are rooted in convex optimization, notably in the work of Wolfe [START_REF] Wolfe | Finding the Nearest Point in A Polytope[END_REF] and Gilbert [START_REF] Gilbert | An Iterative Procedure for Computing the Minimum of a Quadratic Form on a Convex Set[END_REF]. Similarly to GJK, these works aim at solving a Minimum-Norm Point problem and are sub-cases of a more general algorithm: the Frank-Wolfe method, also known as the conditional gradient descent.

The Frank-Wolfe algorithm (FW) dates back to 1956 and is one of the first convex optimization algorithms. It has been heavily studied over the years by the optimization community. This algorithm iterates over the computation of support points in order to approach the optimal solution. The undesired zig-zagging behavior of FW, already identified by its authors, has been addressed by introducing corrections to the original FW method [START_REF] Garber | Playing Non-linear Games with Linear Oracles[END_REF][START_REF] Gilbert | An Iterative Procedure for Computing the Minimum of a Quadratic Form on a Convex Set[END_REF][START_REF] Guélat | Some Comments on Wolfe's 'Away Step[END_REF][START_REF] Jaggi | Revisiting Frank-Wolfe: Projection-free sparse convex optimization[END_REF][START_REF] Kerdreux | Projection-Free Optimization on Uniformly Convex Sets[END_REF][START_REF] Lacoste-Julien | On the Global Linear Convergence of Frank-Wolfe Optimization Variants[END_REF][START_REF] Wolfe | Finding the Nearest Point in A Polytope[END_REF]. In [START_REF] Jaggi | Revisiting Frank-Wolfe: Projection-free sparse convex optimization[END_REF] and [START_REF] Lacoste-Julien | On the Global Linear Convergence of Frank-Wolfe Optimization Variants[END_REF], widely used corrections of the FW algorithm are analyzed as well as their convergence properties. In this work, we notably show in Sec. II that the GJK algorithm is an instance of the fully-corrective Frank-Wolfe algorithm, covered in [START_REF] Lacoste-Julien | On the Global Linear Convergence of Frank-Wolfe Optimization Variants[END_REF], applied to solving a Minimum-Norm Point problem (MNP). Finally, recent works have also tried to accelerate the FW algorithm by applying the so-called Nesterov acceleration [START_REF] Nesterov | A method for solving the convex programming problem with convergence rate O(1/k 2 )[END_REF], a classic acceleration technique in unconstrained optimization. Nesterov momentum has been successfully added by Li et al. [START_REF] Li | A Momentum-Guided Frank-Wolfe Algorithm[END_REF] in order to accelerate FW. In [START_REF] Qin | Smoothing Algorithms for Computing the Projection onto a Minkowski Sum of Convex Sets[END_REF], Qin and An take a different approach as they are interested in the general problem of projecting a point onto a Minkowski difference in any dimension. In order to accelerate the theoretical convergence of the 1966 Gilbert algorithm, the authors devise the NESMINO algorithm which exploits the classic Nesterov acceleration. However, the minimization problem (1) is modified by introducing a smoothing term. By doing so, the authors rely on successive projections onto the original shapes instead of computing support points. This makes the NESMINO algorithm similar to the projected-gradient descent method. Unlike NESMINO, GJK and our work are sub-cases of the FW algorithm. As pointed out by the authors, although the NESMINO algorithm makes use of the Nesterov acceleration, it does not accelerate over the original 1966 Gilbert algorithm. In Sec. IV, we experimentally show that the NESMINO algorithm is slower when compared to GJK and our method.

Contributions. Our work builds on the seminal works by Frank and Wolfe [START_REF] Frank | An algorithm for quadratic programming[END_REF] and Gilbert et al. [START_REF] Gilbert | A fast procedure for computing the distance between complex objects in three-dimensional space[END_REF] as well as on the work of Li et al. [START_REF] Li | A Momentum-Guided Frank-Wolfe Algorithm[END_REF] to globally accelerate distance computation and collision checking algorithms between convex shapes. We make these three main contributions: → We recast the collision detection problem as a convex optimization problem which can be solved by the FW algorithm. Using the ideas developed by Gilbert, Johnson and Keerthi, we show that GJK is in fact a sub-case of the fully-corrective FW algorithm; → We adapt recent works on Nesterov-accelerated FW in order to accelerate both the distance computation and the Boolean collision check problems; → We empirically analyse the convergence of our proposed approach on two large shape benchmarks. Results

show a faster convergence of our approach leading to a computational time up to two times faster than the state-of-the-art GJK algorithm, on both distance computation and Boolean collision checking. Paper outline. The paper is organized as follows. In Sec. II, we recast the distance computation problem as a Frank-Wolfe instance. We introduce the duality-gap of the FW method, allowing us to bound the distance to the optimal solution of the distance computation problem. We also present the fully-corrective version of FW and show the link between GJK and FW. In Sec. III, we introduce recent work on Nesterov-accelerated FW and show how to adapt it for both distance computation and Boolean collision checking. For distance computation, we adapt the convergence criterion of FW when using Nesterov acceleration in order to retain the bound on the distance to the optimal solution. We also propose to adapt the Nesterov acceleration for non-strictly convex shapes. In Sec. IV, we evaluate our approach against the state-of-the-art GJK algorithm on two benchmarks containing both strictly convex shapes and non-strictly convex shapes.

II. COLLISION DETECTION FROM A FRANK-WOLFE PERSPECTIVE

In this section, we highlight the natural connection between computing the distance between convex shapes and convex optimization, particularly within the frame of the Frank-Wolfe setting. We notably show that the GJK algorithm can be seen as a variant of the Frank-Wolfe algorithm that leverages properties of convex 3D shapes to drastically lower the computational complexity.

Distance computation and Boolean collision checking. As recalled in Sec. I, collision detection is a sub-case of distance computation: dist(A 1 , A 2 ) > 0 means that the two shapes do not overlap while dist(A 1 , A 2 ) = 0 means that the shapes are in collision. In the particular case of dist(A 1 , A 2 ) > 0, it is sufficient to find a strictly positive lower bound on d 1,2 to solve the collision problem. In the context of convex shapes, this is often simpler than computing the distance between the two shapes [START_REF] Stasse | Real-Time (Self)-Collision Avoidance Task on a Hrp-2 Humanoid Robot[END_REF] and can be done by finding a plane separating A 1 from A 2 . In the rest of the paper, we focus on the generic problem of computing the distance between A 1 and A 2 , as it encapsulates the simpler Boolean collision check covered later in this section. Results for the particular Boolean collision checking case are analyzed in the experimental section IV.

Recasting the distance computation problem onto the Minkowski difference. The first important idea of Gilbert et al. [START_REF] Gilbert | A fast procedure for computing the distance between complex objects in three-dimensional space[END_REF] is to recast the distance computation problem onto the Minkowski difference D of the shapes, illustrated in Fig. 1 and defined as follows:

D = A 1 -A 2 = {x = x 1 -x 2 | x 1 ∈ A 1 , x 2 ∈ A 2 } ⊂ C , (2)
where C = R n is the so-called collision space. The shapes A 1 and A 2 lie in the shape space, and the Minkowski difference D lies in the collision space. Although both the shape space and the collision space are in R n , we make a distinction between the two to highlight the change in perspective. We stress that the Minkowski difference D is specific to shapes A 1 and A 2 . If the relative position or relative orientation between A 1 and A 2 changes, their Minkowski difference changes accordingly.

The following properties hold for the Minkowski difference D:

1) Since A 1 and A 2 are convex sets, D is also convex.

2) If A 1 and A 2 are intersecting, the origin of C, denoted as 0 C , lies inside the Minkowski difference D, i.e. 0 C = x 1x 2 for some x 1 ∈ A 1 and x 2 ∈ A 2 . 3) If A 1 and A 2 are not intersecting, the projection of 0 C onto D, x * = proj D (0 C ), corresponds to two vectors x * 1 ∈ A 1 and x * 2 ∈ A 2 , also called witness vectors in the computational geometry literature [START_REF] Ericson | Real-Time Collision Detection[END_REF]. Contrary to x * , these vectors x * 1 and x * 2 are not necessarily unique, as is the case for non-strictly convex shapes such as two parallel boxes. 4) Finally, we always have ∥x * ∥ = dist(A 1 , A 2 ). This final remark allows us to recast the distance computation problem (1) onto the Minkowski difference as follows:

d 1,2 = min x∈D ∥x -0 C ∥ 2 = min x∈D ∥x∥ 2 . ( 3
)
The convex optimization problem (3) is equivalent to [START_REF] Bach | Learning with Submodular Functions: A Convex Optimization Perspective[END_REF] and is known as a Minimum-Norm Point problem in the optimization literature [START_REF] Bach | Learning with Submodular Functions: A Convex Optimization Perspective[END_REF][START_REF] Lacoste-Julien | On the Global Linear Convergence of Frank-Wolfe Optimization Variants[END_REF][START_REF] Wolfe | Finding the Nearest Point in A Polytope[END_REF]. In our case, 0 C ∈ C = R n is the null vector i.e. the origin of the collision space.

We thus aim at finding the point in D with the lowest norm. This vector x * is the optimal solution to (3), given by

d 1,2 = ∥x * ∥ 2 = dist(A 1 , A 2 ) 2 .
Directly computing the Minkowski difference D is neither analytically tractable nor computationally efficient. Most of the first and second order methods for constrained convex optimization problems, such as projected gradient descent or interior point methods [START_REF] Boyd | Convex Optimization[END_REF], are thus sub-optimal choices. However, computing support vectors of the Minkowski difference D, a notion defined hereafter in this section, is relatively simple and largely demonstrated by Gilbert et al. [START_REF] Gilbert | A fast procedure for computing the distance between complex objects in three-dimensional space[END_REF]. Solving convex optimization problems by computing support vectors is the foundation of the Frank-Wolfe algorithm and its variants [START_REF] Jaggi | Revisiting Frank-Wolfe: Projection-free sparse convex optimization[END_REF], as we discuss next.

Distance computation using the Frank-Wolfe algorithm. The Frank-Wolfe algorithm (FW) [START_REF] Frank | An algorithm for quadratic programming[END_REF] is one of the oldest convex optimization methods and solves the following constrained optimization problem:

f (x * ) = min x∈D f (x), (4)
where f : R n → R is a convex and differentiable function and D is a compact convex set. For our distance computation problem (3), we use f (x) = ∥x∥ 2 and the Minkowski difference D as convex constraint set. As a side note, the following discussed algorithms all require an initial starting point x 0 ∈ D. Shapes used in physics engines are usually attached to a frame in order to keep track of their position Algorithm 1 Frank-Wolfe algorithm with linesearch [START_REF] Jaggi | Revisiting Frank-Wolfe: Projection-free sparse convex optimization[END_REF] Let x 0 ∈ D, ϵ > 0 For k=0, 1, ... do

1: d k = ∇f (x k ) ▷ Direction of support 2: s k ∈ arg min s∈D ⟨d k , s⟩(= S D (d k )) ▷ Support (7) 3: If g F W (x k ) ≤ ϵ , return f (x k ) ▷ Duality gap (15) 4: γ k = arg min γ∈[0,1] f (γx k + (1 -γ)s k ) ▷ Linesearch 5: x k+1 = γ k x k + (1 -γ k )s k
▷ Update iterate In the case of the distance computation problem (3), where f (x) = ∥x∥ 2 , line 4-5 correspond to projecting 0 C on the segment [x k , s k ]:

6: x k+1 = proj [x k ,s k ] (0 C ) ▷ Project 0 C on [x k , s k ] Fig. 2: Computing a support vector s k in direction ∇f (x k )
on convex set D. We illustrate with the example of distance computation. On the left, we draw the Minkowski difference D which point of minimum norm (MNP) is

x * i.e. x * is the projection of 0C onto D, x * = proj D (0C)
. The iterate at iteration k of the FW algorithm is x k . In purple we draw the level sets of the function f (x) = ∥x∥ 2 . On the right, we draw in purple the level sets of the linearization of f at iterate x k , h k . The first step of the FW algorithm is to compute support vector s k in the direction of ∇f (x k ) (green arrow), s k ∈ SD(∇f (x k )). In the second step of the FW algorithm, we compute x k+1 as a convex combination of x k and s k i.e. x k+1 is a point on the segment

[x k , s k ].
and orientation in space. We denote c 1 ∈ A 1 and c 2 ∈ A 2 the origins of the frames attached to A 1 and A 2 , respectively. In the rest of this paper, we take

x 0 = c 1 -c 2 .
The FW algorithm, summarized in Alg. 1, is a gradient-descent method. It consists in iteratively applying two steps in order to converge towards the optimal solution x * of ( 4). If we denote by x k the estimate of x * at iteration k, these two steps correspond to:

1) First, we compute a support vector s k in the direction of ∇f (x k ), by solving a linear optimization problem on D. 2) Second, we update our current iterate x k to obtain x k+1 , by taking a convex combination of the current iterate x k and the computed support vector s k .

In the following we detail these steps in the context of distance computation. At iteration k, the current iterate x k is the estimate of the optimal solution x * and f (x k ) is the estimate of the optimal value of (4), f (x * ), at iteration k. We write the linearization of the function f at x k and denote it as h k :

h k (s) = f (x k ) + ⟨∇f (x k ), s -x k ⟩ ( 5
)
where s is a vector of R n , ∇f (x k ) is the gradient of f at x k and ⟨•, •⟩ denotes the dot product between two vectors of R n . → Step 1. The first step of the FW algorithm at iteration k consists of finding a minimizer s k ∈ D of h k on the convex set D (line 2 in Alg. 1). Such a vector s k is called a support vector of D or simply a support and is defined as follows:

s k ∈ arg min s∈D h k (s) = arg min s∈D ⟨∇f (x k ), s⟩ . (6)
Fig. 2 gives a graphical understanding of support s k . The vector s k belongs to D and is in the most opposite direction w.r.t. ∇f (x k ). In order to highlight the importance of the direction in which a support s k is computed, we now introduce the notion of support direction and support function. Given a support direction d ∈ R n , the support function S D returns a set of D and is defined as:

S D (d) = arg min s∈D ⟨d, s⟩ ⊂ D . (7)
The support function S D may return a set with more than one vector. We only need to use one vector of this set. Thinking in terms of direction of support allows us to understand that this direction can be rescaled while preserving the output of the support function:

∀d ∈ R n , ∀α > 0, S D (αd) = S D (d). (8)
A support

s k ∈ D at iteration k is thus computed in the direction d k = ∇f (x k ) and belongs to S D (∇f (x k )), s k ∈ S D (∇f (x k )).
We now explain how to compute the support vector s k in the case of the distance computation problem (3) where we minimize f (x) = ∥x∥ 2 on the Minkowski difference D of A 1 and A 2 . First, we have ∇f (x) = 2x. Therefore, in the case of problem (3), it follows that:

s k ∈ S D (x k ) = arg min s∈D ⟨x k , s⟩ . (9)
As demonstrated by Gilbert et al. [START_REF] Gilbert | A fast procedure for computing the distance between complex objects in three-dimensional space[END_REF], any vector s ∈ S D (d) related to the Minkowski difference can be decomposed as the difference between two support vectors s A1 ∈ S A1 (d) and s A2 ∈ S A2 (-d) over the two individual shapes, leading to the following relation:

s = s A1 -s A2 ∈ S D (d). (10)
Equation [START_REF] Ericson | Real-Time Collision Detection[END_REF] shows that we can construct a support of the Minkowski difference from the supports of the original shapes. This property highlights the powerful change of perspective of working on the Minkowski difference. Indeed, there exists a large number of shapes for which computing supports is simple: spheres, ellipsoids, cylinders, capsules, polytopes etc. [START_REF] Ericson | Real-Time Collision Detection[END_REF][START_REF] Van Den Bergen | A Fast and Robust GJK Implementation for Collision Detection of Convex Objects[END_REF][START_REF] Van Den Bergen | Collision Detection in Interactive 3D Environments[END_REF]. Fig. 3 illustrates the construction of a support of the Minkowski difference D using the supports of the original shapes A 1 and A 2 . → Step 2. Once a support vector s k ∈ S D (x k ) has been computed, we update the iterate x k to obtain x k+1 by taking a convex combination between s k and x k . The original FW algorithm uses a parameter-free update: where γ k = k+1 k+2 ∈ [0, 1] controls the step size. Alternatively, a linesearch can be carried out to find a better iterate x k+1 (line 4 in Alg. 1):

x k+1 = γ k x k + (1 -γ k )s k , (11)
γ k = arg min γ∈[0,1] f (γx k + (1 -γ)s k ) x k+1 = γ k x k + (1 -γ k )s k . ( 12
)
In the distance computation case where f (x) = ∥x∥ 2 , this linesearch ( 12) is equivalent to projecting

0 C onto the segment [x k , s k ], x k = proj [x k ,s k ] (0 C ) (line 4 in Alg. 1).
Since D is convex, both ( 11) and ( 12) updates are guaranteed to remain in D.

Stopping criteria. As Frank-Wolfe deals with convex problems, the duality gap associated to problem (4) can be used as a stopping criterion. Due to its convexity, the function f is always above its linearization. Otherwise said, for any x ∈ R n and any s ∈ R n :

f (s) ≥ f (x) + ⟨∇f (x), s -x⟩ . ( 13
)
Reworking this inequality and applying the min operator enables us to compute the Frank-Wolfe duality gap g FW (x) ∈ R + which gives an upper-bound on the difference f (x)f (x * ):

f (x) -f (x * ) ≤ -min s∈D ⟨∇f (x), s -x⟩ = g FW (x) . (14)
In particular, at iteration k of the FW algorithm, we have:

f (x k ) -f (x * ) ≤ g FW (x k ) = ⟨∇f (x k ), x k -s k ⟩ , ( 15
)
where

s k ∈ S D (∇f (x k ))
is the support vector computed at iteration k in the direction of ∇f (x k ). The duality-gap g FW (x k ) serves as a convergence criterion for the Frank-Wolfe method and is cheap to compute. Applied to the distance computation problem (3), the duality gap at iteration k, g FW (x k ), guarantees that:

∥x k ∥ 2 -∥x * ∥ 2 ≤ g FW (x k ) = 2⟨x k , x k -s k ⟩. ( 16
)
Using the triangular inequality of the Euclidian norm and the convexity of the Minkowski difference D, we can show that:

∥x k -x * ∥ 2 ≤ ∥x k ∥ 2 -∥x * ∥ 2 ≤ g FW (x k ) . ( 17
)
Algorithm 2 Boolean collision checking: separating plane condition Insert after line 2 in Alg. 1:

1: If ⟨d k , s k ⟩ > 0 , return False If after termination d 1,2 = 0, return True
Inequality ( 17) is useful in practice as it allows to finely control the desired tolerance on the distance to the optimal solution x * (line 3 in Alg. 1). Indeed, if ones wants to compute an estimate x of the optimal solution x * at precision ϵ, meaning that ∥xx * ∥ ≤ √ ϵ, it is sufficient to check that g FW (x) ≤ ϵ.

Boolean collision checking. As mentioned earlier, the problem of distance computation encompasses the problem of collision checking. Indeed, in collision checking, we are only interested in finding a separating plane between A 1 and A 2 , if it exists. This is equivalent to finding a separating plane between D and 0 C . For any support direction d, if we have:

⟨d, s⟩ > 0, s ∈ S D (d), (18)
then the plane supported by the vector d separates D from 0 C [START_REF] Van Den Bergen | A Fast and Robust GJK Implementation for Collision Detection of Convex Objects[END_REF]. This also means that, in the case where the two shapes intersect, collision checking has the same computational complexity than distance computation. As shown in Alg. 2, we add this separating plane condition before line 2 in Alg. 1.

The Frank-Wolfe active-set. As many gradient-descent algorithms, the FW method tends to zig-zag towards the optimal solution [START_REF] Lacoste-Julien | On the Global Linear Convergence of Frank-Wolfe Optimization Variants[END_REF], slowing down the convergence to the optimum. This behavior is undesired and amplified if the optimal solution x * lies close to the boundary of the constraint set D. In collision detection, this corresponds to the case where the two shapes are not intersecting. This zig-zagging behavior is due to the way that Frank-Wolfe approaches the set of active constraints [START_REF] Lacoste-Julien | On the Global Linear Convergence of Frank-Wolfe Optimization Variants[END_REF], also called active-set in the optimization literature [START_REF] Boyd | Convex Optimization[END_REF]. In the FW setting, the active-set at iteration k, denoted W k = {s 0 , ..., s r } ⊂ D, is the set of vectors in D used by the algorithm to maintain a convex combination of the iterate x k :

x k = r i=0 λ i s i , r i=0 λ i = 1 with s i ∈ W k ⊂ D and λ i > 0. (19
) In Alg. 3, we rewrite the FW algorithm with linesearch (Alg. 1) in order to highlight the notion of active-set: to projecting 0 C onto the segment [x k , s k ] (line 5 in Alg. 3).

• Finally, the active-set is updated W k+1 = {x k+1 } (line 6 in Alg. 3). In practice, discarding previously computed supports when updating the active-set is inefficient and causes the zig-zagging phenomenon observed in the FW algorithm [START_REF] Lacoste-Julien | On the Global Linear Convergence of Frank-Wolfe Optimization Variants[END_REF]. In the optimization literature, a rich and wide variety of variants of the FW algorithm have been introduced to efficiently cope with the active-set in order to improve the convergence rate of the FW method [START_REF] Garber | Playing Non-linear Games with Linear Oracles[END_REF][START_REF] Guélat | Some Comments on Wolfe's 'Away Step[END_REF][START_REF] Holloway | An Extension of the Frank and Wolfe Method of Feasible Directions[END_REF][START_REF] Kerdreux | Restarting Frank-Wolfe[END_REF][START_REF] Wolfe | Finding the Nearest Point in A Polytope[END_REF]. However, these variants remain too generic and are not suited for the specific problem of collision detection. Instead, we propose next to incorporate the active-set strategy used in GJK within the Frank-Wolfe setting.

Connection between GJK and Frank-Wolfe. In the case of collision detection, Gilbert et al. [START_REF] Gilbert | A fast procedure for computing the distance between complex objects in three-dimensional space[END_REF] developed an efficient strategy to handle the active-set at a minimal cost. To represent the current estimate x k and the optimal solution x * , GJK exploits the concept of simplexes in R 3 . A simplex in R n corresponds to a set containing at most n + 1 vectors of R n and the rank r of a simplex is the number of vectors it contains (0 < r ≤ n + 1). For 3-dimensional spaces, a simplex corresponds either to a point (r = 1), a segment (r = 2), a triangle (r = 3) or a tetrahedron (r = 4). Similarly to the simplex methods for Linear Programming [START_REF] Dantzig | Linear programming and extensions[END_REF], the Carathéodory theorem [START_REF] Caratheodory | Uber den Variabilitatsbereich der Koeffizienten von Potenzreihen, die gegebene Werte nicht annehmen[END_REF] motivates the use of simplexes. Let Y be a set of N ≥ n vectors in R n , Y = {y i ∈ R n } 0≤i≤N . The Carathéodory theorem states that any vector x ∈ conv(Y) can be expressed as the convex combination of at most n + 1 vectors of Y:

x = r j=0 λ j y j , with y j ∈ Y, λ j > 0, r i=0 λ j = 1. (20)
Hence, any vector in D, and particularly the optimal solution x * ∈ D = of the distance computation problem (3), can be identified as a convex combination of the vectors composing a simplex W . Relying on simplexes is attractive as there is no need to run any algorithm to compute the convex-hull of a simplex as they are convex by construction. Frank-Wolf algorithms may operate on more complex active-sets, which might become hard to tackle from a computational point of view [START_REF] Jaggi | Revisiting Frank-Wolfe: Projection-free sparse convex optimization[END_REF][START_REF] Lacoste-Julien | On the Global Linear Convergence of Frank-Wolfe Optimization Variants[END_REF]. In other words, the problem of finding the optimal solution x * can be reformulated as the problem of identifying the optimal simplex W * on which x * can be decomposed into a convex combination. This is exactly the approach followed by GJK that we now detail. At iteration k of GJK, the current iterate x k is a convex combination of the vectors composing the simplex W k of rank r k ≤ n. To update x k and W k , the following procedure is applied:

• After computing support vector s k (line 2 in Alg. 3), we add s k to W k to obtain W k+1 = W k ∪ {s k } (line 4 in Alg. 3). The set W k+1 is now a simplex of rank r k+1 ≤ n + 1.
Algorithm 3 Frank-Wolfe algorithm with line-search (see Alg. 1) rewritten with active-sets and applied to the distance computation problem (3)

Let x 0 ∈ D, W 0 = {x 0 }, ϵ > 0 For k=0, 1, ... do 1: d k = x k ▷ Direction of support 2: s k ∈ S D (d k ) ▷ Support (7) 3: If g F W (x k ) ≤ ϵ , return f (x k ) ▷ Duality gap (15) 4: W k+1 = W k ∪ {s k } ▷ Augment active-set 5: x k+1 = proj conv( W k+1 ) (0 C ) ▷ Project 0 C on conv( W k+1 ) 6: W k+1 = {x k+1 } ▷ Update active-set • We then minimize function f (x) = ∥x∥ 2 on W k+1
to obtain x k+1 , corresponding to projecting 0 C onto W k+1 :

x k+1 = proj conv( W k+1 ) (0 C ) 1 (line 5 in Alg.
3). • We then have two cases, summarized in Alg 4: → If x k+1 = 0 C , the algorithm is stopped. Thus, we have x * = 0 C and d 1,2 = 0 in (3) (line 1 in Alg. 4). → Otherwise, we construct W k+1 from W k+1 . To do so, we retain only the minimal number of vectors in W k+1 needed to express x k+1 as a convex combination (line 2 in Alg. 4). Indeed, as 0 C / ∈ W k+1 , the projection x k+1 of 0 C on W k+1 necessarily lies on a face of W k+1 , and can be expressed as a convex combination of the vectors composing this face. This ensures that W k+1 is necessarily of rank r k+1 < r k+1 ≤ n + 1. Through this discussion, it is clear that GJK is a particular case of Frank-Wolfe. More specifically, it is a sub-case of the fully-corrective Frank-Wolfe algorithm analyzed by Lacoste-Julien and Jaggi [START_REF] Lacoste-Julien | On the Global Linear Convergence of Frank-Wolfe Optimization Variants[END_REF]. The strategy used by GJK to handle the active-set has proved to be very efficient in practice and renders the GJK algorithm state of the art for collision detection. In the next section, we propose to leverage the formulation of collision detection as a Frank-Wolfe sub-case to accelerate its convergence following the well-established Nesterov acceleration paradigm [START_REF] Nesterov | A method for solving the convex programming problem with convergence rate O(1/k 2 )[END_REF].

III. ACCELERATING COLLISION DETECTION

Gradient descent (GD) is the backbone of many convex optimization methods and relies solely on the gradient of the objective function. Second order methods [START_REF] Boyd | Convex Optimization[END_REF] such as Newton methods have faster convergence rates than GD at the price of requiring the computation and the inversion of the Hessian quantities. Momentum methods have thus been introduced in the optimization literature to provide gradient-based methods with improved convergence rates without requiring costly Hessian evaluation. In this section we use recent work linking the Nesterov acceleration of GD to the FW algorithm to globally accelerate collision detection. This global acceleration Algorithm 4 Fully-corrective FW using simplexes, applied to the distance computation problem (3). This algorithm is identical to GJK [START_REF] Gilbert | A fast procedure for computing the distance between complex objects in three-dimensional space[END_REF] In Alg. 3, let W 0 = ∅ and replace line 6 by:

1: If x k+1 = 0 C , return 0
If the algorithm has not terminated, update W k+1 to retain only the smallest number of vectors needed to express x k+1 : 2: W k+1 = {s 1 , ..., s r } where s 1 , ..., s r are the smallest number of vectors in W k+1 such that x k+1 is a convex combination of s 1 , ..., s r .

of collision detection is experimentally evaluated in Sec. IV on several benchmarks.

Nesterov acceleration for unconstrained optimization. We initially consider the following unconstrained minimization problem:

f (x * ) = min x∈R n f (x) , (21)
where f : R n → R is a convex and differentiable function.

The vanilla gradient-descent algorithm consists of following the slope of f given by its gradient ∇f . The following scheme is applied iteratively until a given convergence criterion is met (e.g., ∥∇f (x k )∥ < ϵ, with ϵ the desired precision):

x k+1 = x k + α k ∇f (x k ), (22)
where x k ∈ R n is the current iterate and α k ∈ R is the gradient step. This standard setting leads to a simple implementation with linear convergence rate (O(1/k)).

To go beyond this linear convergence regime, acceleration techniques have been devised in the optimization community to provide quadratic convergence rate (O(1/k 2 )) or more [START_REF] Alexandre D'aspremont | Acceleration methods[END_REF], by relying on relatively cheap gradient evaluations. Among these gradient-descent acceleration techniques, the Nesterov acceleration [START_REF] Nesterov | A method for solving the convex programming problem with convergence rate O(1/k 2 )[END_REF] is one of the better studied and most popular in practice [START_REF] Alexandre D'aspremont | Acceleration methods[END_REF]. It is based on accumulating previously computed gradients in a momentum term d k using this momentum d k to update the current iterate x k as:

y k = x k + δ k d k-1 (23a) d k = δ k d k-1 + α k ∇f (y k ) (23b) x k+1 = x k + d k (23c)
where δ k ∈ R is the momentum parameter and y k ∈ R n is an intermediate quantity. The role of momentum d k is to smooth the trajectory of iterates converging towards the optimum by geometrically averaging previously computed gradients.

The δ k momentum parameter is selected to prevent damping or overshooting of the iterate trajectory when going towards the optimal solution x * . Note that the current iterate x k is extrapolated using the momentum term d k to compute the intermediate vector

y k = x k + δ k d k .
The gradient is then computed at the vector y k . Computing the term y k leads to an anticipatory behavior in similar spirit to extra-gradient methods [START_REF] Alexandre D'aspremont | Acceleration methods[END_REF].

Algorithm 6 Nesterov-accelerated Frank-Wolfe [START_REF] Li | A Momentum-Guided Frank-Wolfe Algorithm[END_REF] In Alg. 1 and Alg. 3, let d -1 = s -1 = x 0 , δ k = k+1 k+3 and replace line 1 by:

1: y k = δ k x k + (1 -δ k )s k-1 2: d k = δ k d k-1 + (1 -δ k )∇f (y k )
The Frank-Wolfe algorithm and Nesterov acceleration. Recent works of Li et al. [START_REF] Li | A Momentum-Guided Frank-Wolfe Algorithm[END_REF][START_REF] Li | Heavy Ball Momentum for Conditional Gradient[END_REF] have proposed to adapt the Nesterov acceleration to the FW setting. We propose to leverage and adapt this FW acceleration scheme to the context of collision detection, by notably extending the FW formulation of collision detection previously developed in Sec. II.

In the original FW algorithm, the support vector at iteration k, s k , is computed in the direction of the gradient ∇f (x k ) (line 1 in Alg. 1). In the Nesterov acceleration of FW proposed by Li et al. [START_REF] Li | A Momentum-Guided Frank-Wolfe Algorithm[END_REF], the direction of support for computing s k is instead defined by:

y k = δ k x k + (1 -δ k )s k-1 (24a) d k = δ k d k-1 + (1 -δ k )∇f (y k ) (24b) s k = S D (d k ), (24c)
where

s k-1 is the support computed at iteration k -1, δ k = k+1 k+3 ∈ [0, 1]
is the momentum parameter and S D is the support function as defined in [START_REF] Coumans | PyBullet, a Python module for physics simulation for games, robotics and machine learning[END_REF]. As in the Nesterov acceleration for GD, y k is an intermediary vector used to evaluate the gradient ∇f (y k ). To ensure y k stays in D, it is a convex combination of x k and s k-1 , both vectors of D. The direction of support is then obtained by taking a convex combination of the previous support direction d k-1 and the gradient ∇f (y k ).

Li et al. [START_REF] Li | A Momentum-Guided Frank-Wolfe Algorithm[END_REF] have experimentally shown that this acceleration strategy leads to a better convergence rate of the FW algorithm when compared to the original FW algorithm. In the following, we explain how to adapt the Nesterov acceleration of FW to collision detection.

Adapting Nesterov fully-corrective Frank-Wolfe to distance computation. Preserving GJK's simplex strategy is crucial for collision detection as it greatly speeds up the vanilla FW algorithm. Therefore, we adapt (24) accordingly as:

y k = δ k x k + (1 -δ k )s k-1 (25a) d k = δ k d k-1 + (1 -δ k )∇f (y k ) (25b) s k = S D (d k ), (25c)
W k+1 = W k ∪ {s k }, (25d) x k+1 = proj conv( W k+1 ) (0 C ). (25e)
These steps are also summarized in Alg. 7. The update of simplex W k+1 from W k+1 is then identical to the one described in Alg. 4. The original duality gap defined in Sec. II (Eq. 15) can no longer be used as a convergence

Algorithm 7 Nesterov-accelerated GJK Let x 0 ∈ D, W 0 = ∅, d -1 = s -1 = x 0 , ϵ > 0 For k=0, 1, ... do 1: δ k = k+1 k+3 ▷ Momentum parameter value 2: y k = δ k x k + (1 -δ k )s k-1 ▷ Intermediary point (25a) 3: d k = δ k d k-1 + (1 -δ k )∇f (y k ) ▷ Support dir. (25b) 4: s k ∈ S D (d k )
▷ Support (7) 5: if g(x k ) ≤ ϵ then ▷ Fixed-point condition [START_REF] Nesterov | A method for solving the convex programming problem with convergence rate O(1/k 2 )[END_REF] 6:

If d k = x k , return f (x k ) ▷ Algorithm terminates 7: s k ∈ S D (∇f (x k )) ▷ Compute s k in dir. ∇f (x k ) Replace line 3 by: d k = x k until termination. 8: W k+1 = W k ∪ {s k } ▷ Augment active-set 9: x k+1 = proj conv( W k+1 ) (0 C ) ▷ Project 0 C on conv( W k+1 ) 10: If x k+1 = 0 C ,
return 0 11: W k+1 = {s 1 , ..., s r } where s 1 , ..., s r are the smallest number of vectors in W k+1 such that x k+1 is a convex combination of s 1 , ..., s r . criterion. Indeed, the following inequality:

∥x k -x * ∥ 2 ≤ g FW (x k ) = 2⟨x k , x k -s k ⟩, s k ∈ S D (x k ),
is no longer valid because the support vector s k is no longer computed in the direction of the gradient ∇f (x k ) = 2x k . Next we will show that the original stopping criterion devised in Sec. II cannot be used and we need to derive a new one.

Stopping criterion. As the number of iteration k increases, δ k → k→∞ 1 in [START_REF] Mamou | A Simple and Efficient Approach for 3D Mesh Approximate Convex Decomposition[END_REF]. Therefore, d k tends to be equal to d k-1 (25b) and thus s k = s k-1 (25c). As a consequence, augmenting W k with s k to construct W k+1 (see (25d)) and then projecting 0 C onto W k+1 (25e) will not result in any progress. Therefore, x k+1 = x k : the algorithm reaches a fixed point and is stuck on constant support direction d.

In order to cope with this issue, we use the following strategy. Suppose

x k ̸ = 0 C . Since x k = proj conv(W k ) (0 C ) we have: ∀s i ∈ W k , ⟨x k , x k -s i ⟩ = 0. ( 26
)
After computing s k ∈ S D (d k ), if we have:

⟨x k , x k -s k ⟩ ̸ = 0, (27)
then s k is not a linear combination of vectors in W k . Therefore, augmenting W k with s k to obtain W k+1 and projecting 0 C onto conv( W k+1 ) to obtain x k+1 will result in the algorithm progressing toward the optimum x * . Suppose on the contrary that:

⟨x k , x k -s k ⟩ = 0, (28)
then s k is a linear combination of vectors in W k . Adding s k to W k will thus not result in any progress towards the optimum. As a consequence, Eq. ( 28) encompasses two cases:

• If the support direction d k is aligned with ∇f (x k ), Eq. ( 28), corresponding to g FW (x k ) = 0, matches the termination criteria of the distance computation problem and therefore we have reached the optimum leading to x k = x * . • Otherwise, if d k is not aligned with ∇f (x k ), the algorithm cannot stop as a null duality gap is not met. The algorithm thus enters a cycle where it iterates until Eq. ( 28) does not hold. To cope with this undesired behavior we simply stop the Nesterov acceleration as soon as Eq. ( 28) is met and switch back to the non-accelerated version Alg. 4. We thus define the function g such that for any s k ∈ D:

g(x k ) = 2⟨x k , x k -s k ⟩, (29)
g is used in Alg. 7 as an optimality criterion ( g ≤ ϵ) either for stopping the Nesterov acceleration in order to continue with the vanilla GJK, or as stopping criteria qualifying an optimal solution, in which case g = g FW and ( 17) holds. The entire algorithm is summarized in Alg. 7.

Nesterov acceleration for non-strictly convex shapes. Let us explain the effect of the Nesterov acceleration on the support direction update (25b) and distinguish between strictly convex and non-strictly convex D: . Therefore, even though δ k is not close to 1, the fixed point condition ( 28) can be verified. The Nesterov acceleration is stopped, possibly prematurely. The latter case is especially problematic when shapes A 1 and A 2 are in close-proximity, which is ultimately the type of collision problems which are commonly encountered in simulation or motion planning with contacts. In (25b), this is due to the norm of d k-1 being greatly predominant over the norm of y k as k increases, ∥d k-1 ∥ ≫ ∥∇f (y k )∥. To prevent this phenomenon on non-strictly convex D, we propose to replace (25b) by a simple heuristic which normalizes the gradient and momentum directions as follows:

• If D is
d k = δ k d k-1 ∥d k-1 ∥ + (1 -δ k ) ∇f (y k ) ∥∇f (y k )∥ , (30)
summarized in Alg. 8. In Sec. IV, we experimentally prove this heuristic to significantly reduce the number of steps for distance computations for non-strictly convex shapes.

Algorithm 8 Normalize direction for non-strictly convex shapes Replace line 3 in Alg. 7 by:

1: d k = δ k d k-1 ∥d k-1 ∥ + (1 -δ k ) ∇f (y k ) ∥∇f (y k )∥

IV. EXPERIMENTS

In this section, we study the performance of Nesterov-accelerated GJK (Alg. 7) against the vanilla GJK (Alg. 4) algorithm. We use the HPP-FCL C++ library [START_REF] Pan | FCL: A General Purpose Library for Collision and Proximity Queries[END_REF][START_REF] Pan | HPP-FCL: an extension of the Flexible Collision Library[END_REF] and its implementation of GJK as a starting point for the implementation of the proposed Nesterov-accelerated GJK (Alg. 7). To distinguish between pairs of strictly convex and non-strictly convex shapes, we build a first benchmark only composed of pairs of ellipsoids (strictly convex shapes) and a second benchmark using pairs of standard meshes (represented by their convex hulls) which are taken from the commonly-used ShapeNet dataset [START_REF] Angel | ShapeNet: An Information-Rich 3D Model Repository[END_REF]. Finally, we empirically show that the simplex strategy used by GJK and our method (discussed in Sec. II) is crucial for efficient collision detection. We show that GJK and our method significantly outperform the original FW algorithm and the recent NESMINO [START_REF] Qin | Smoothing Algorithms for Computing the Projection onto a Minkowski Sum of Convex Sets[END_REF] algorithm. Although it differs from FW algorithms, we include the NESMINO algorithm in this analysis as its projected-gradient descent procedure is accelerated using the classic Nesterov acceleration scheme [START_REF] Nesterov | A method for solving the convex programming problem with convergence rate O(1/k 2 )[END_REF]. The code to run the benchmarks is made freely available at: https://github.com/lmontaut/collision-detection-benchmark. Shape datasets. In the ellipsoids benchmark, the ellipsoids are randomly generated by sampling positive-definite matrices. In total, we generate 1000 random pairs of ellipsoids. Given a pair of ellipsoids, we randomly sample relative poses between the two shapes such that both objects do not intersect. This then allows us to finely control the distance dist(A 1 , A 2 ) between the objects, enabling us to measure the of the separation distance on the performance of the studied algorithms. The values used for dist(A 1 , A 2 ) range from -0.1 m to 1 m. Negative values correspond to scenarii where the shapes are intersecting and dist(A 1 , A 2 ) is the norm of the separating vector. The separating vector is the vector of smallest norm needed to translate one of the two shapes such that the two shapes do not intersect. Therefore, for each pair of ellipsoids, 100 random relative poses are sampled such that the shapes do not intersect. For each relative pose, we translate the shapes along the axis given by their closest-points to study the impact of dist(A 1 , A 2 ). We then set dist(A 1 , A 2 ) to fixed values between -0.1 m to 1 m. On the other hand, the mesh dataset contains roughly 35k shapes but we only use 1.8k randomly selected shapes from 10 categories of the dataset to set-up the mesh benchmark. To generate the mesh benchmark, we sample 10 k random pairs of meshes and represent each shape by its convex hull. The resulting meshes contain between 10 and 5000 vertices. About 50% of meshes contain between 100 and 1000 vertices. The remaining shapes are equally distributed between less than 100 vertices or more than 1000 vertices. As in the ellipsoids benchmark, 100 random relative poses are sampled for each pair such that the shapes do not intersect and then set dist(A 1 , A 2 ) to fixed values between -0.1 m and 1 m. In both benchmarks, the characteristic sizes of the shapes range from a few centimeters up to a meter. Finally, for the distance computation problem, we select a convergence tolerance of ϵ = 10 -8 .

Metrics.

In order to measure the performances of the Nesterov-accelerated GJK and the vanilla GJK algorithms, we measure the number of iterations N k to solve a given collision problem. For the mesh benchmark, we also measure the execution time T µ of both methods. To cope with CPU throttling, we solve generated collision problems 100 times. We then report the average of the 90% lowest computation times. Since both algorithms have identical iterations which differ only in the computation of the support direction, the number of iterations is preferred to study the performance of the algorithms; the execution time serves as a reference but can vary based on the implementation and the hardware used.

A. Strictly convex shapes: ellipsoids

Single collision pair. To qualitatively understand the effect of the Nesterov acceleration, we first study the evolution of the optimality criterion [START_REF] Holloway | An Extension of the Frank and Wolfe Method of Feasible Directions[END_REF] in the case of distance computation on a single pair of ellipsoids (Fig. 4).

In Fig. 4a, the shapes are intersecting and both algorithms converge to x * = 0 C . The optimality criterion ( 29) is not used by any of the two methods to terminate. Indeed, they stop when the projection onto the simplex W k yields 0 C (line 10 in Alg. 7 and line 1 in Alg. 4). However, the criterion is still used in Nesterov-accelerated GJK to switch to GJK. Both methods converge rapidly. As explained in Sec. II, this is expected as the optimal solution x * = 0 C lies inside D.

In Fig. 4b, the shapes are close together with a separating distance of dist(A 1 , A 2 ) = 0.01 m. Along with the previous intersecting case, this figure represents a typical scenario one could encounter in the narrow phase of the collision detection problem. Contrary to the intersecting case, the solution x * now lies on the boundary of D. The convergence of GJK is impaired due to the successive switches in support directions. By smoothing the change of support direction, the Nesterov acceleration requires less iterations and divides by more than two the number of iterations to reach the convergence criterion compared to GJK.

When shapes are distant (Fig. 4c), the problem is easier to solve for GJK. This is due to smaller angles between successive support directions which is a consequence of the large distance between 0 C and D. In such a case, the Nesterov acceleration is now detrimental to the convergence of the algorithm as it induces too little successive changes in the support directions. This is a typical case where smoothing the gradient directions via the Nesterov acceleration does not help and takes more than triple the number of iterations for Nesterov-accelerated GJK to achieve convergence compared For Nesterov-accelerated GJK, the switch from Nesterov acceleration to vanilla GJK is denoted by the switch in color from orange to green. In (a), the optimality criterion defined by ( 29) is not used by any of the two methods to terminate as 0C lies inside D and both algorithms terminate when x k = 0C. It is however used by Nesterov-accelerated GJK to stop the Nesterov acceleration when a fixed-point is met, hence the sudden peak denoting the switch to GJK. Nesterov acceleration has beneficial effects in the intersecting (a) and close-proximity (b) setups which are the cases of interest, for instance, in the context of simulation and planning with contacts. to vanilla GJK. This behavior is however specific to ellipsoids and is not observed in meshes as shown below. Fig. 4 highlights three different behaviors in the convergence of both algorithms: when the shapes are intersecting, in close-proximity and distant. Scenarii where shapes are intersecting or in close proximity are of greater importance in general, as they correspond to contexts which are evaluated during the narrow phase. On the opposite, distant scenarii are rejected by the broad phase, and other approximations of dist(A 1 , A 2 ) could be used if needed [START_REF] Ericson | Real-Time Collision Detection[END_REF].

Statistical validation over the ellipsoids dataset. The convergence graphs in Fig. 4 directly depend on the pair of shapes considered. Therefore, to get a better statistical understanding of the performance of Nesterov-accelerated GJK and vanilla GJK, we focus on the ellipsoid benchmark.

Fig. 5a and Fig. 5b show the mean and standard-deviation of the number of iterations N k of each method for the distance computation and the Boolean collision checking problems, respectively. In the case of distance computation (Fig. 5a), the results confirm the trend shown in the convergence graphs in Fig. 4. When the shapes are shallowly intersecting, Nesterov-accelerated GJK converges with the same or even fewer number of iterations than vanilla GJK. The shallower the penetration, the more Nesterov accelerates over vanilla GJK. The irregularity in standard deviation at -0.01 m is a critical zone for the momentum where the variance increases. This could be fixed thanks to a restart strategy [START_REF] Alexandre D'aspremont | Acceleration methods[END_REF] and will be analysed in a future work. When shapes are in close-proximity, the Nesterov acceleration of GJK significantly reduces the number of iterations compared to vanilla GJK. Finally, when shapes are distant, 1 m ≤ dist(A 1 , A 2 ), Nesterov-accelerated GJK is detrimental to convergence on the distance computation problem. A similar pattern of speed-ups of the Nesterov-accelerated GJK over the vanilla GJK is shown for the collision detection problem in Fig. 5b.

B. Non-strictly convex shapes: meshes

Effect of support direction normalization. For meshes, the importance of normalizing the support direction (see Eq. ( 30)) in the Nesterov-accelerated GJK is highlighted in Fig. 6. For both the distance computation and Boolean collision checking problems, it prevents the Nesterov acceleration from reaching a fixed-point too early and consequently it reduces the overall amount of iterations needed to converge. In the following, we thus focus only on Nesterov-accelerated GJK with support normalization and compare its performance against the vanilla GJK algorithm.

Statistical validation over the ShapeNet dataset. In Fig. 7 and Fig. 8, we report the number of iterations N k and execution time T µ for both Nesterov-accelerated and vanilla GJK. In Fig. 9, we report relative accelerations N k GJK /N k Nesterov and T µ GJK /T µ Nesterov of Nesterov-accelerated GJK compared to GJK on a given collision problem. These relative measures allow to analyze the effects of the two algorithms on exactly the same collision problems, which are not captured when using absolute values. Overall, Nesterov-accelerated GJK significantly reduces the number of iterations and execution time when compared to GJK in cases where shapes are shallowly intersecting or in close-proximity. The number of iterations and the execution times are reduced by a factor ranging from 1.5 up to 2. It is worth recalling at this stage, that when two shapes are relatively far from each other, any broad phase algorithm will automatically discard such pair. Only in a ) on the ShapeNet benchmark for close-by or shallowly intersecting shapes where the distance between shapes is -0.1m ≤ dist(A1, A2) ≤ 0.1m. N denotes the number of vertices for each mesh. For all depicted pairs, Nesterov-accelerated GJK always gives a better (or at least as good) timing for both collision computation and distance evaluation. The graphs show the number of iterations N k (lower is better) for GJK and Nesterov-accelerated GJK with and without support direction normalization when the two shapes are in close-proximity: 0 m < dist(A1, A2) ≤ 0.1 m. Using normalization of support direction is beneficial in Nesterov-accelerated GJK, reducing the overall number of iterations compared to GJK and non-normalized Nesterov-accelerated GJK.

small percentage of cases Nesterov-accelerated GJK is slower than GJK. Momentum restarts could help improve these rare cases and will be considered in future work. In Table . I, we select 4 meshes with an increasing number of vertices to highlight the benefits of the Nesterov acceleration. For each pair, we report the mean and the standard deviation of the execution time for both distance computation and Boolean collision checking. We consider the challenging set-up of close-by or shallowly intersecting shapes in the range of separation distances -0.1 m ≤ dist(A 1 , A 2 ) ≤ 0.1 m. The lower mean and standard deviation show that Nesterov-accelerated GJK is faster than the vanilla GJK and reduces the spread of computation times across the different collision problems in

C. Importance of the simplex strategy in GJK

In the previous sections, we have experimentally shown the improvement of our method, Nesterov-accelerated GJK, over the vanilla GJK algorithm for collision problems which are important in practice, i.e. when collision pairs have not been filtered by the broad phase and are thus overlapping or in close-proximity. We conclude this section by demonstrating the importance of the simplex strategy used in GJK and our method when solving collision problems.

To do so, we evaluate the performance of the Frank-Wolfe algorithm (Alg. 1), the recent NESMINO algorithm [START_REF] Qin | Smoothing Algorithms for Computing the Projection onto a Minkowski Sum of Convex Sets[END_REF], GJK and our method on ellipsoids and cubes and report the results in Table . II. Altough the NESMINO algorithm is similar to projected-gradient descent and strongly differs from Frank-Wolfe like algorithms, it uses the classic Nesterov acceleration which makes it interesting to compare to our method. FW, GJK and our method stop when a tolerance of ϵ = 10 -8 on the FW duality-gap is met. Therefore, to render the NESMINO algorithm comparable to the other considered methods, we run NESMINO until the distance between its solution and the solution found by GJK is less than √ ϵ = 10 -4 . In Table IIa, we consider 1000 collision problems between pairs of ellipsoids for each distance category (overlapping, close-proximity and distant). Shapes are in close-proximity when 0 ≤ dist(A 1 , A 2 ) ≤ 0.1 m. In this first scenario, all algorithms have a comparable number of operations per iteration. Indeed, the projection operation used in NESMINO when the shapes are ellipsoids has the same complexity as the support operation used in the three other algorithms. Although GJK and our method also do a simplex projection at each iteration, this operation has about the same complexity as computing the support point. In the case of strictly-convex shapes such as ellipsoids, GJK and our method significantly outperform the FW and NESMINO algorithms. This is especially the case when the shapes are overlapping or in close-proximity where GJK or our method take 3 to 10 times less iterations compared to FW or NESMINO.

In Table IIb, we repeated the same experiments with collision pairs of cubes. Since cubes are polytopes with a small number of vertices, GJK and our method only take a few iterations to reach a tolerance of ϵ = 10 -8 . However, because cubes are non-strictly convex shapes, the convergence of the FW algorithm is O(1/ϵ), i.e. it takes on the order of 1/ϵ iterations to reach a FW duality-gap of ϵ. The NESMINO algorithm takes less iterations than FW but more than 100 times more iterations than GJK and our method. In the specific case of polytopes, the NESMINO algorithm is also much more costly per iteration than FW, GJK or our method, as it replaces the computation of support points by much more costly projections on the original polytopes.

V. CONCLUSION

In this work, we have first established that the well-known GJK algorithm can be understood as a variant of the Frank-Wolfe method, well studied within the convex optimization community, and more precisely GJK can be identified as a sub-case of fully-corrective Frank-Wolfe. Subsequently, this connection has enabled us to accelerate the GJK algorithm in the sense of Nesterov acceleration, by adapting recent contributions on the application of Nesterov acceleration to the context of Frank-Wolfe. Through extensive benchmarks, we have shown that this acceleration is beneficial for both collision detection and distance computation settings for scenarios where shapes are intersecting or in close proximity, accelerating collision detection by up to a factor of 2. Interestingly, these two scenarios notably encompass the generic contexts of planning and control as well as physical simulation, which are essential areas of modern robotics. Therefore, although the proposed acceleration corresponds to improvements of GJK's execution time on the order of few microseconds, modern robotics applications may solve millions to billions of collision problems e.g. when learning a policy with RL [START_REF] Openai | Solving Rubik's Cube with a Robot Hand[END_REF]. The Nesterov acceleration of GJK is already included in the HPP-FCL library [START_REF] Pan | HPP-FCL: an extension of the Flexible Collision Library[END_REF], notably used by the HPP framework [START_REF] Mirabel | HPP: A new software for constrained motion planning[END_REF] for motion planning and the Pinocchio framework [START_REF] Carpentier | The Pinocchio C++ library -A fast and flexible implementation of rigid body dynamics algorithms and their analytical derivatives[END_REF] dedicated to simulation and modelling. One can expect this work to be largely adopted in the current available GJK implementations, as it only requires minor algorithmic changes. This work should benefit a large audience within robotics (e.g., simulation, planning, control) and beyond by addressing issues shared by other communities, including computer graphics and computational geometry.

  (a) Left: two distant shapes. Right: 0 C lies outside of the Minkowski difference of the shapes. (b) Left: two overlapping shapes. Right: 0 C lies inside of the Minkowski difference of the shapes.

Fig. 1 :

 1 Fig. 1: Distant vs. overlapping pairs of shapes and their respective Minkowski difference. Left column: two convex shapes in 2D. Right column: the Minkowski difference D of A1 and A2. Since A1 and A2 are convex, D is also convex. In (a), the shapes are not in collision hence the origin of the configuration space C, 0C (in red) lies outside the Minkowski difference, 0C ̸ ∈ D. The vector x * = x * 1 -x * 2 separates A1 from A2. It is also equal to the projection of 0C onto the Minkowski difference D, x * = proj D (0C). In (b), the shapes are overlapping, thus 0C ∈ D. In this case, we have x * = proj D (0C) = 0C.

Fig. 3 :

 3 Fig.3: Computing a support vector on the Minkowski difference using support vectors on the individual shapes. The vector sA 1 is a support vector of shape A1 in direction d. The vector sA 2 is a support vector of shape A2 in direction -d. The constructed vector s = sA 1 -sA 2 is a support vector of the Minkowski difference D in the direction D.

  strictly convex, any vector s belonging to the surface of D has a unique corresponding direction d such that s = S D (d). Here, we stress the fact that the support function S D returns only one vector. Consequently, we have d k ̸ = d k-1 and therefore s k ̸ = s k-1 . The fixed point condition (28) is thus not met unless δ k = 1 and Nesterov acceleration continues to be applied in Alg. 7. In practice, the algorithm runs until δ k gets close to 1 or x k gets close to 0 C . The condition (28) is then satisfied as the algorithm starts to cycle. The Nesterov acceleration is thus removed and the algorithm runs until the convergence criteria is satisfied, guaranteed by the Frank-Wolfe algorithm. • Otherwise, if D is non-strictly convex, multiple support directions {d 1 , ..., d m , ...} can yield the same support vector s ∈ S D (d 1 ) = ... = S D (d m ) = ... etc. Consequently, it is possible to have d k-1 ̸ = d k and s k = s k-1

1 ,

 1 A 2 ) = 10 -2 [m] GJK GJK + Nesterov (b) Close-proximity ellipsoids.

Fig. 4 :

 4 Fig. 4: Convergence of GJK and Nesterov-accelerated GJK for ellipsoid collisions. The y-axis measures the optimality criterion (29).

- 10 - 1 -

 101 10 -3 0 10 -3 10 -1 10 0 dist(A 1 , A 2 ) [m]

Fig. 5 :

 5 Fig.5: Comparison of Nesterov-accelerated GJK and vanilla GJK on the ellipsoid benchmark. The graphs show the number of iterations (y-axis) vs. the signed distance between the two shapes (x-axis). The curve shows the mean value over 100,000 random trials. The shaded region corresponds to the standard deviation. The Nesterov-accelerated GJK algorithm requires fewer iterations when the shapes are in close proximity.

Fig. 6 :

 6 Fig. 6: Importance of support direction normalization in Nesterov-accelerated GJK on the ShapeNet benchmark.The graphs show the number of iterations N k (lower is better) for GJK and Nesterov-accelerated GJK with and without support direction normalization when the two shapes are in close-proximity: 0 m < dist(A1, A2) ≤ 0.1 m. Using normalization of support direction is beneficial in Nesterov-accelerated GJK, reducing the overall number of iterations compared to GJK and non-normalized Nesterov-accelerated GJK.

  Number of iterations N k . Lower is better. Execution time T . Lower is better.

Fig. 7 :

 7 Fig. 7: Distance computation on the ShapeNet benchmark. The graphs show the number of iterations (a) and the execution time (b) for Nesterov-accelerated GJK with normalization (red) and vanilla GJK (green) for a range of distances (x-axis) between the shapes.this setting.

  Number of iterations N k . Lower is better.-10 -1 -10 -2 -10 -3 10 -3 10 -2 10 -1 10 0 dist(A 1 , A 2 ) [m] Execution time T µ . Lower is better.

Fig. 8 :

 8 Fig. 8: Boolean collision check on the ShapeNet benchmark. The graphs show the number of iterations (a) and the execution time (b) for the Nesterov-accelerated GJK with normalization (red) and vanilla GJK (green) for a range of distances (x-axis) between the shapes.

Fig. 9 :

 9 Fig. 9: Speed-ups on the ShapeNet benchmark. The plots show ratios of the number of iterations (a) and execution times (b) of vanilla GJK and the Nesterov-accelerated GJK with normalization. Ratios over 1.0 show speed-ups of Nesterov-accelerated GJK over GJK.

TABLE I :

 I Computation times (µs) for distance computation (T µ

D ) and Boolean collision checking (T µ C

TABLE II :

 II Number of iterations for distance computation between ellipsoids (a) and between cubes (b). For each distance category (overlapping, close-proximity and distant), we report the mean and standard deviation of the number of iterations. Shapes are in close-proximity when 0 ≤ dist(A1, A2) ≤ 0.1 m.

• A iteration k, the active-set is only composed of x k , W k = {x k }. • The active-set W k is then augmented by computing a support s k (line 2 in Alg. 3) to obtain W k+1 = {x k , s k } (line 4 in Alg. 3). • We then minimize function f on the convex-hull of W k+1 , conv( W k+1 ), which is simply the segment [x k , s k ]. For the distance computation problem (3), this linesearch operation is equivalent

The efficient projection onto simplexes in R 3 , named the distance sub-algorithm by Gilbert et al.[START_REF] Gilbert | A fast procedure for computing the distance between complex objects in three-dimensional space[END_REF], is thoroughly covered in[START_REF] Ericson | Real-Time Collision Detection[END_REF][START_REF] Van Den Bergen | A Fast and Robust GJK Implementation for Collision Detection of Convex Objects[END_REF] and its robustness is improved in[START_REF] Montanari | Improving the GJK Algorithm for Faster and More Reliable Distance Queries Between Convex Objects[END_REF].

ACKNOWLEDGMENTS

The authors would like to warmly thank Francis Bach for fruitful discussions on Min-Norm Point algorithms. This work was partly supported by the European Regional Development Fund under the project IMPACT (reg. no. CZ.02.1.01/0.0/0.0/15 003/0000468), by the French government under management of Agence Nationale de la Recherche as part of the "Investissements d'avenir" program, reference ANR-19-P3IA-0001 (PRAIRIE 3IA Institute) and the Louis Vuitton ENS Chair on Artificial Intelligence.