Olivier Fercoq Ltci
email: olivier.fercoq@telecom-paristech.fr

A generic coordinate descent solver for nonsmooth convex optimization

Keywords: Coordinate descent, convex optimization, generic solver, efficient implementation

We present a generic coordinate descent solver for the minimization of a nonsmooth convex objective with structure. The method can deal in particular with problems with linear constraints. The implementation makes use of efficient residual updates and automatically determines which dual variables should be duplicated. A list of basic functional atoms is pre-compiled for efficiency and a modelling language in Python allows the user to combine them at run time. So, the algorithm can be used to solve a large variety of problems including Lasso, sparse multinomial logistic regression, linear and quadratic programs.

Introduction

Coordinate descent methods decompose a large optimization problem into a sequence of onedimensional optimization problems. The algorithm was first described for the minimization of quadratic functions by Gauss and Seidel in [START_REF] Seidel | Ueber ein Verfahren, die Gleichungen, auf welche die Methode der kleinsten Quadrate führt, sowie lineäre Gleichungen überhaupt, durch successive Annäherung aufzulösen[END_REF]. Coordinate descent methods have become unavoidable in machine learning because they are very efficient for key problems, namely Lasso [START_REF] Friedman | Pathwise coordinate optimization[END_REF], logistic regression [START_REF] Yu | Dual coordinate descent methods for logistic regression and maximum entropy models[END_REF] and support vector machines [START_REF] John | Fast Training of Support Vector Machines using Sequential Minimal Optimization[END_REF][START_REF] Shalev | Stochastic Dual Coordinate Ascent Methods for Regularized Loss Minimization[END_REF]. Moreover, the decomposition into small subproblems means that only a small part of the data is processed at each iteration and this makes coordinate descent easily scalable to high dimensions.

One of the main ingredients of an efficient coordinate descent solver is its ability to compute efficiently partial derivatives of the objective function [START_REF] Nesterov | Efficiency of coordinate descent methods on huge-scale optimization problems[END_REF]. In the case of least squares for instance, this involves the definition of a vector of residuals that will be updated during the run of the algorithm. As this operation needs to be performed at each iteration, and millions of iterations are usually needed, the residual update and directional derivative computation must be coded in a compiled programming language.

Many coordinate descent solvers have been written in order to solve a large variety of problems. However, most of the existing solvers can only solve problems of the type min

x∈R N J j=1 f (A j x -b j ) + I i=1 g(x (i) )
where x (i) ∈ R Ni is the ith block of x, I i=1 N i = N , A j ∈ R Mj×N is a matrix and b j ∈ R Mj is a vector, and where f is a convex differentiable function and g is a convex lower- semicontinuous function whose proximal operator is easy to compute (a.k.a. a proximal-friendly convex function). Each piece of code usually covers only one type of function [START_REF] Fan | Liblinear: A library for large linear classification[END_REF][START_REF] Pedregosa | Scikit-learn: Machine learning in python[END_REF]. Moreover, even when the user has a choice of objective function, the same function is used for every block [START_REF] Blondel | Lightning: large-scale linear classification, regression and ranking in Python[END_REF].

In this work, we propose a generic coordinate descent method for the resolution of the convex optimization problem

min x∈R N 1 2 x ⊤ Qx + J j=1 c f j f j (A f j x -b f j ) + I i=1 c g i g i (D g i x (i) -b g i ) + L l=1 c h l h l (A h l x -b h l ) . (1)
We shall call f j , g i and h l atom functions. Each of them may be different. We will assume that f j 's are differentiable and convex, g i 's and h l 's are proximal-friendly convex functions. As before A f j ∈ R M f j ×N and A h l ∈ R M h l ×N are matrices, D g i is a multiple of the identity matrix of size

N i , b f j ∈ R M f j , b g i ∈ R
Ni and b h l ∈ R M h l are vectors, c f j , c g i and c h l are positive real numbers, Q is a N × N positive semi-definite matrix.

The algorithm we implemented is described in [START_REF] Fercoq | A coordinate-descent primal-dual algorithm with large step size and possibly nonseparable functions[END_REF] and can be downloaded on https://bitbucket.org/ofercoq/cd_solver. The present paper focuses on important implementation details about residual updates and dual variable duplication. The novelty of our code is that it allows a generic treatment of these algorithmic steps and includes a modelling interface in Python for the definition of the optimization problem. Note that unlike most coordinate descent implementations, it can deal with nonseparable nonsmooth objectives and linear constraints.

Literature review on coordinate descent methods

A thorough review on coordinate descent is beyond the scope of this paper. We shall refer the interested reader to the review papers [START_REF] Stephen | Coordinate descent algorithms[END_REF] and [START_REF] Hao-Jun | A primer on coordinate descent algorithms[END_REF]. Instead, for selected papers dealing with smooth functions, separable non-smooth functions or non-separable non-smooth function, we list their main features. We also quickly review what has been done for non-convex functions. We sort papers in order of publication except when there are an explicit dependency between a paper and a follow-up.

Smooth functions

Smooth objectives are a natural starting point for algorithmic innovations. The optimization problem at stake writes

min x∈R N f (x)
where f is a convex differentiable function with Lipschitz-continuous partial derivatives.

In Table 1, we compare papers that introduced important improvements to coordinate descent methods. We shall in particular stress the seminal paper by Tseng and Yun [START_REF] Tseng | A coordinate gradient descent method for nonsmooth separable minimization[END_REF]. It features coordinate gradient steps instead of exact minimization. Indeed a coordinate gradient steps gives similar iteration complexity both in theory and in practice for a much cheaper iteration cost. Moreover, this opened the door for many innovations: blocks of coordinates and the use of proximal operators were developed in the same paper. Another crucial step was made in [START_REF] Nesterov | Efficiency of coordinate descent methods on huge-scale optimization problems[END_REF]: Nesterov showed how randomization can help finding finite-time complexity bounds and proposed an accelerated version of coordinate descent. He also proposed to use a non-uniform sampling of coordinates depending on the coordinate-wise Lipschitz constants.

Separable non-smooth functions

A large literature has been devoted to composite optimization problems with separable nonsmooth functions:

min x∈R N f (x) + n i=1 g i (x (i) )
where f is a convex differentiable function with Lipschitz-continuous partial derivatives and for all i, g i is a convex function whose proximal operator is easy to compute. Indeed, regularized expected risk minimization problems often fit into this framework and this made the success of coordinate descent methods for machine learning applications. Some papers study min x∈R p f (x) + n i=1 g i ((Ax) (i) ), where f is strongly convex and apply coordinate descent to a dual problem written as

min y∈R N f * (-A ⊤ y) + n i=1 g * i (y (i) ) .
One of the challenges of these works is to show that even though we are solving the dual problem, one can still recover rates for a sequence minimizing the primal objective. We present our selection of papers devoted to this type of problems in Table 2.

Non-separable non-smooth functions

Non-separable non-smooth objective functions are much more challenging to coordinate descent methods. One wishes to solve where f is a convex differentiable function with Lipschitz-continuous partial derivatives, g and h are convex functions whose proximal operator are easy to compute and A is a linear operator. Indeed, the linear operator introduces a coupling between the coordinates and a naive approach leads to a method that does not converge to a minimizer [START_REF] Auslender | Optimisation: Méthodes Numériques[END_REF]. When h = ι {b} , the convex indicator function of the set {b}, we have equality constraints. We present our selection of papers devoted to this type of problems in Table 3.

min x∈R N f (x) + g(x) + h(Ax)

Non-convex functions

The goal here is to find a local minimum to the problem loc-min

x∈R N f (x) + n i=1 g i (x (i) )
without any assumption on the convexity of f nor g i . The function f should be continuously differentiable and the proximal operator of each function g i should be easily computable. Note that in the non-convex setting, the proximal operator may be set-valued.

Description of the Algorithm

General scheme

The algorithm we implemented is a coordinate descent primal-dual method developed in [START_REF] Fercoq | A coordinate-descent primal-dual algorithm with large step size and possibly nonseparable functions[END_REF]. Let introduce the notation

F (x) = 1 2 x ⊤ Qx + J j=1 c f j f j (A f j x -b f j ), G(x) = I i=1 c g i g i (D g i x (i) -b g i ), H(z) = L l=1 c h l h l (z (l) -b h l ), J (i) = {j : A h j,i = 0}, I(j) = {i : A h j,i = 0}, m j = |I(j)| and ρ(A) the spectral radius of matrix A. We shall also denote J f (i) = {j : A f j,i = 0},

Paper

Rate Const PD-CD Notable feature Platt '99 [START_REF] John | Fast Training of Support Vector Machines using Sequential Minimal Optimization[END_REF] × ✓ P for SVM Tseng & Yun '09 [START_REF] Tseng | Block-Coordinate Gradient Descent Method for Linearly Constrained Non Separable Optimization[END_REF] ✓ ✓ P adapts Gauss-Southwell rule Tao et al '12 [START_REF] Tao | Stochastic Coordinate Descent Methods for Regularized Smooth and Nonsmooth Losses[END_REF] ✓ × P uses averages of subgradients Necoara et al '12 [START_REF] Necoara | Efficiency of randomized coordinate descent methods on optimization problems with linearly coupled constraints[END_REF] ✓ ✓ P 2-coordinate descent Nesterov '12 [START_REF] Nesterov | Subgradient Methods for Huge-Scale Optimization Problems[END_REF] ✓ × P uses subgradients Necoara & Clipici '13 [START_REF] Necoara | Efficient parallel coordinate descent algorithm for convex optimization problems with separable constraints: application to distributed MPC[END_REF] ✓ ✓ P coupled constraints Combettes & Pesquet '14 [START_REF] Patrick | Stochastic Quasi-Fejér Block-Coordinate Fixed Point Iterations with Random Sweeping[END_REF] × ✓ ✓ 1st PD-CD, short step sizes Bianchi et al '14 [START_REF] Bianchi | A stochastic coordinate descent primaldual algorithm and applications[END_REF] × ✓ ✓ distributed optimization Hong et al '14 [START_REF] Hong | A block successive upper bound minimization method of multipliers for linearly constrained convex optimization[END_REF] × ✓ × updates all dual variables Fercoq & Richtárik '17 [START_REF] Fercoq | Smooth minimization of nonsmooth functions with parallel coordinate descent methods[END_REF] ✓ × P uses smoothing Alacaoglu et al '17 [START_REF] Alacaoglu | Smooth primaldual coordinate descent algorithms for nonsmooth convex optimization[END_REF] ✓ ✓ ✓ 1st PD-CD w rate for constraints Xu & Zhang '18 [START_REF] Xu | Accelerated primal-dual proximal block coordinate updating methods for constrained convex optimization[END_REF] ✓ ✓ × better rate than [21] Chambolle et al '18 [START_REF] Chambolle | Stochastic primal-dual hybrid gradient algorithm with arbitrary sampling and imaging applications[END_REF] ✓ ✓ × updates all primal variables Fercoq & Bianchi '19 [START_REF] Fercoq | A coordinate-descent primal-dual algorithm with large step size and possibly nonseparable functions[END_REF] ✓ ✓ ✓ 1st PD-CD w long step sizes Gao et al '19 [START_REF] Gao | Randomized primal-dual proximal block coordinate updates[END_REF] ✓ ✓ × 1st primal-dual w rate for constraints Latafat et al '19 [START_REF] Latafat | A new randomized blockcoordinate primal-dual proximal algorithm for distributed optimization[END_REF] ✓ ✓ ✓ linear conv w growth condition Table 3. Selected papers for the minimization of non-smooth non-separable functions. Rate: does the paper prove rates? Const: can the method solve problems with linear equality constraints? PD-CD: the method is purely primal (P), the method updates some primal variables but all the Lagrange multipliers or some dual variables but all the primal variables (×), the method updates some primal and some dual variables at each iteration (✓).

Paper

J Q (i) = {j : Q j,i = 0}, A f ∈ R j M f j ×N
the matrix which stacks the matrices (A f j ) 1≤j≤J and A h the matrix which stacks the matrices (A h l ) 1≤l≤L . The algorithm writes then as Algorithm 1.

Algorithm 1 Coordinate-descent primal-dual algorithm with duplicated variables (PD-CD)

Input: Differentiable function F : R N → R, matrix A h ∈ R M h ×N
, functions G and H whose proximal operators are available.

Initialization: Choose x 0 ∈ R N , y 0 ∈ R nnz(A h ) . Denote J (i) = {j : A h j,i = 0}, I(j) = {i : A h j,i = 0}, m j = |I(j)| and ρ(A) the spectral radius of matrix A. Choose step sizes τ ∈ R I + and σ ∈ R L + such that ∀i ∈ {1, . . . I}, τ i < 1
β i + ρ j∈J (i) m j σ j (A h ) ⊤ j,i A h j,i . (2)
For all i ∈ {1, . . . , I}, set w

(i) 0 = j∈J (i) (A h ) ⊤ j,i y (j)
0 (i). For all j ∈ {1, . . . , J}, set z

(j) 0 = 1 mj i∈I(j) y (j) 0 (i).
Iteration k: Define:

y k+1 = prox σ,H ⋆ z k + D(σ)A h x k x k+1 = prox τ,G x k -D(τ ) ∇F (x k ) + 2(A h ) ⊤ y k+1 -w k .
For i = i k+1 ∼ U ({1, . . . , I}) and for each j ∈ J (i k+1 ), update:

x (i) k+1 = x (i) k+1 y (j) k+1 (i) = y (j) k+1 w (i) k+1 = w (i) k + j∈J(i) (A h ) ⊤ j,i (y (j) k+1 (i) -y (j) k (i)) z (j) k+1 = z (j) k + 1 m j (y (j) k+1 (i) -y (j) k (i)) .
Otherwise, set x

(i ′ ) k+1 = x (i ′ ) k , w (i ′ ) k+1 = w (i ′ ) k , z (j ′ ) k+1 = z (j ′ ) k
and y

(j ′ ) k+1 (i ′ ) = y (j ′ ) k (i ′ ).
We will denote U 1 , . . . , U I the columns of the identity matrix corresponding to the blocks of x = (x (1) , . . . , x (I) ), so that U i x (i) ∈ R N and V 1 , . . . , V J the columns of the identity matrix corresponding to the blocks of

A f x -b f = (A f 1 x -b f 1 , . . . , A f J x -b f J ).

Computation of partial derivatives

For simplicity of implementation, we are assuming that G is separable and the blocks of variable will follow the block structure of G. This implies in particular that at each iteration, only ∇ i F (x k ) needs to be computed. This partial derivative needs to be calculated efficiently because it needs to be performed at each iteration of the algorithm. We now describe the efficient residual update method, which is classically used in coordinate descent implementations [START_REF] Nesterov | Efficiency of coordinate descent methods on huge-scale optimization problems[END_REF].

Denote r f,x k = A f x k -b f and r Q,x k = Qx k .
By the chain rule, we have

∇ i F (x k ) = J j=1 c f j (A f ) ⊤ j,i ∇f j (A f j x k -b f j ) + Qx k = j∈J f (i) c f j (A f ) ⊤ j,i ∇f j ((r f,x k ) j ) + j∈J Q (i) U i r Q,x k If r f,x k and r Q,x k are pre-computed, only O(|J f (i)| + |J Q (i)|
) operations are needed. For an efficient implementation, we will update the residuals r f,x k as follows, using the fact that only the coordinate block i k+1 is updated:

r f,x k+1 = A f x k+1 -b f = A f x k + U ik+1 (x (ik+1) k+1 -x (ik+1) k ) -b f = r f,x k + A f U ik+1 (x (ik+1) k+1 -x (ik+1) k ) = r f,x k + j∈J f (ik+1) V j A f j,ik+1 (x (ik+1) k+1 -x (ik+1) k ) Hence, updating r f,x k+1 also requires only O(|J f (i k+1 )|) iterations. Similarly, updating the residuals r Q,x k , r h,x k = A h x k -b h , w k and z k can be done in O(|J Q (i k+1 )|) and O(|J (i k+1 )|) operations.
Although this technique is well known, it is not trivial how to write it in a generic fashion, since residual updates are needed at each iteration and should be written in a compiled language. We coded the residual update using abstract atom functions in order to achieve this goal.

Computation of proximal operators using atom functions

Another major step in the method is the computation of the i th coordinate of prox τ,G (x) for a given x ∈ R N .

As D g is assumed to be diagonal, G is separable. Hence, by the change of variable z = D g i x-b g i , (prox τ,G (x)) i = arg min

x∈R N i c g i g i (D g i x -b g i ) + 1 2τ i x -x (i) 2 = (D g i ) -1 b g i + arg min z∈R N i c g i g i (z) + 1 2τ i (D g i ) -1 (b g i + z) -x (i) 2 = (D g i ) -1 b g i + arg min z∈R N i g i (z) + 1 2c g i (D g i ) 2 τ i z -(D g i x (i) -b g i ) 2 = (D g i ) -1 b g i + prox c g i (D g i ) 2 τigi (D g i x (i) -b g i )
where we used the abuse of notation that D g i is either the scaled identity matrix or any of its diagonal elements. This derivation shows that to compute (prox τ,G (x)) i we only need linear algebra and the proximal operator of the atom function g i .

We can similarly compute prox H. To compute prox σ,H ⋆ , we use Moreau's formula:

prox σ,H ⋆ (z) = z -D(σ) prox σ -1 ,H (D(σ) -1 z)

Duplication of dual variables

Algorithm 1 maintains duplicated dual variables y k ∈ R nnz(A h ) as well as averaged dual variables

z k ∈ R M h where M h = L l=1 M h l and A h l,i is of size M h l × N i .
The sets J (i) for all i are given by the sparse column format representation of A h . Yet, for all i, we need to construct the set of indices of y k+1 that need to be updated. This is the table dual vars to update in the code. Moreover, as H is not separable in general, in order to compute ȳj k+1 , for j ∈ J (i k+1 ), we need to determine the set of dual indices j ′ that belong to the same block as j with respect to the block decomposition of H. This is the purpose of the tables inv blocks h and blocks h.

The procedure allows us to only compute the entries of ȳk+1 that are required for the update of y k . The code is organized in nine files. The main file is cd solver.pyx. It contains the Python callable and the data structure for the problem definition. The other files are atoms.pyx/pxd, algorithm.pyx/pxd, helpers.pyx/pxd and screening.pyx/pxd. They contain the definition of the atom functions, the algorithms and the functions for computing the objective value. In Figure 1, we show for each subfunction, in which function it is used. The user needs to call the Python class Problem and the Python function coordinate descent. Atom functions can be added by the user without modifying the main algorithm.

Code structure

All tables are defined using Numpy's array constructor in the coordinate descent function. The main loop of coordinate descent and the atom functions are pre-compiled for efficiency.

Atom functions

The code allows us to define atom functions independently of the coordinate descent algorithm. As an example, we provide in Figure 2 As inputs, it gets x (an array of numbers which is the point where the operation takes place), buff (the buffer for vectorial outputs), nb coord (is the size of x), mode, prox param and prox param2 (numbers which are needed when computing the proximal operator). The input mode can be:

• GRAD in order to compute the gradient.

• PROX to compute the proximal operator.

• PROX CONJ uses Moreau's formula to compute the proximal operator of the conjugate function.

• LIPSCHITZ to return the Lipschitz constant of the gradient.

• VAL CONJ to return the value of the conjugate function. As this mode is used only by compute smoothed gap for printing purposes, its implementation is optional and can be approximated using the helper function val conj not implemented. Indeed, for a small

ǫ > 0, h * (y) = sup z z, y -h(z) ≈ sup z z, y -h(z) -ǫ 2 z 2 = p, y -h(p) -ǫ 2 p 2
, where p = prox h/ǫ (y/ǫ).

• VAL to return the value of the function. Some functions naturally require multi-dimensional inputs, like • 2 or the log-sum-exp function. For consistency, we define all the atoms with multi-dimensional inputs: for an atom function f 0 : R → R, we extend it to an atom function f : R Ni → R by f (x) = Ni l=1 f 0 (x l ). For efficiency purposes, we are bypassing the square atom function when computing a gradient and implemented it directly in the algorithm.

Modelling language

In order to use the code in all its generality, we defined a modelling language that can be used to define the optimization problem we want to solve [START_REF] Alacaoglu | Smooth primaldual coordinate descent algorithms for nonsmooth convex optimization[END_REF].

The user defines a problem using the class Problem. Its arguments can be:

• N the number of variables, blocks the blocks of coordinates coded in the same fashion as the indptr index of sparse matrices (default [0, 1, . . . , N]), x init the initial primal point (default 0) and y init the initial duplicated dual variable (default 0) • Lists of strings f, g and h that code for the atom functions used. The function string to func is responsible for linking the atom function that corresponds to the string. Our convention is that the string code is exactly the name of the function in atoms.pyx.

The size of the input of each atom function is defined in blocks f, blocks and blocks h. The function strings f, g or h may be absent, which means that the function does not appear in the problem to solve. • Arrays and matrices cf, Af, bf, cg, Dg, bg, ch, Ah, bh, Q. The class initiator transforms matrices into the sparse column format and checks whether Dg is diagonal.

For instance, in order to solve the Lasso problem, min x 1 2 Ax -b 2 2 + λ x 1 , one can type pb lasso = cd solver.Problem(N=A.shape [START_REF] Alacaoglu | Smooth primaldual coordinate descent algorithms for nonsmooth convex optimization[END_REF], f=["square"]*A.shape[0], Af=A, bf=b, cf=[0.5]*A.shape[0], g=["abs"]*A.shape [START_REF] Alacaoglu | Smooth primaldual coordinate descent algorithms for nonsmooth convex optimization[END_REF], cg=[lambda]*A.shape [START_REF] Alacaoglu | Smooth primaldual coordinate descent algorithms for nonsmooth convex optimization[END_REF]) cd solver.coordinate descent(pb lasso)

Extensions

Non-uniform probabilities

We added the following feature for an improved efficiency. Under the argument sampling='kink half', the algorithms periodically detects the set of blocks I kink such that i ∈ I kink if x (i) is at a kink of g i . Then, block i is selected with probability law

P(i k+1 = i) =      1 n if |I kink | = n 1 2n if |I kink | < n and i ∈ I kink 1 2n + 1 2(n-|Ikink|) if |I kink | < n and i ∈ I kink
The rationale for this probability law is that blocks at kinks are likely to incur no move when we try to update them. We thus put more computational power for non-kinks. On the other hand, we still keep an update probability weight of at least 1 2n for each block, so even in unfavourable cases, we should not observe too much degradation in the performance as compared to the uniform law.

Acceleration

We also coded accelerated coordinate descent [START_REF] Fercoq | Accelerated, parallel and proximal coordinate descent[END_REF], as well as its restarted [START_REF] Fercoq | Restarting the accelerated coordinate descent method with a rough strong convexity estimate[END_REF] and primaldual [START_REF] Alacaoglu | Smooth primaldual coordinate descent algorithms for nonsmooth convex optimization[END_REF] variants. The algorithm is given in Algorithm 2. As before, ȳk+1 and xk+1 should not be computed: only the relevant coordinates should be computed.

The accelerated algorithms improve the worst case guarantee as explained in Table 5: However, accelerated algorithms do not take profit of regularity properties of the objective like strong convexity. Hence, they are not guaranteed to be faster, even though restart may help.

Algorithm 2 Smooth, accelerate, randomize the Coordinate Descent (APPROX/SMART-CD)

Input: Differentiable function F : R N → R, matrix A h ∈ R M h ×N
, functions G and H whose proximal operators are available.

Initialization: Choose x 0 ∈ R N , ẏ0 ∈ R M h . Choose γ 1 > 0 and denote B i 0 = β i + ρ((A h :,i ) ⊤ A h :,i ) γ1 . Set s = 0, θ 0 = 1 n , c 0 = 1, x0 = 0 ∈ R N and x0 = x 0 . Iteration k: Define: y k+1 = prox γ -1 k+1 ,H ⋆ ẏs + D(γ k+1 ) -1 (c k A h xk + A h xk ) x k+1 = prox θ 0 θ k B -1 k ,G xk - θ 0 θ k D(B k ) -1 ∇F (c k xk + xk ) + (A h ) ⊤ y k+1 .
For i = i k+1 ∼ U ({1, . . . , I}), update:

x(i) k+1 = x (i) k+1 x(i) k+1 = x(i) k - 1 -θ k /θ 0 c k (x (i) k+1 - x(i) k ) Otherwise, set x (i ′ ) k+1 = x (i ′ )
k . Compute θ k+1 ∈ (0, 1) as the unique positive root of

θ 3 + θ 2 + θ 2 k θ -θ 2 k = 0 if h = 0 θ 2 + θ 2 k θ -θ 2 k = 0 if h = 0 Update γ k+2 = γk+1 1+θk+1 , c k+1 = (1 -θ k+1 )c k and B i k+1 = β i + ρ((A h :,i ) ⊤ A h :,i ) γk+2 for all i. If Restart(k) is true: Set x k+1 = xk+1 + c k xk+1 . Set ẏs+1 = y k+1 and s ← s + 1. Reset xk+1 ← 0, xk+1 ← x k+1 , c k+1 ← 1, θ k+1 ← θ 0 , β k+1 ← β 1 . h = 0 h = 0 PD-CD Alg. 1 O(1/k) O(1/ √ k) APPROX / SMART-CD Alg. 2 O(1/k 2 ) O(1/k)

Variable screening

The code includes the Gap Safe screening method presented in [START_REF] Ndiaye | Safe Optimization Algorithms for Variable Selection and Hyperparameter Tuning[END_REF]. Note that the method has been studied only for the case where h = 0. Given a non-differentiability point x(i) of the function g i where the subdifferential ∂g i (x (i) ) has a non-empty interior, a test is derived to check whether x(i) is the i th variable of an optimal solution. If this is the case, one can set x (i) = x(i) and stop updating this variable. This may lead to a huge speed up in some cases. As the test relies on the computation of the duality gap, which has a nonnegligible cost, it is only performed from time to time.

In order to state Gap Safe screening in a general setting, we need the concept of polar of a support function. Let C be a convex set. The support function of C is σ C defined by

σ C (x) = sup y∈C y, x . The polar to σ C is σ • C (x) = sup x:σC (x)≤1
x, x .

In particular, if x ∈ C, then σ • C (x) ≤ 1. Denote f (z) = J j=1 c f j f j (z (j) -b f j ), so that (∇f (A f x)) j = c f j ∇f j (A f j x -b f j ), G i (x) = c g i g i (D g i x (i) -b g i ),
x ⋆ a solution to the optimization problem (1) and suppose we have a set R such that (∇f

(A f x ⋆ -b f ), Qx ⋆ ) ∈ R. Gap Safe screening states that max (ζ,ω)∈R σ • ∂Gi(x (i) ) ((A f ) ⊤ i ζ + ω (i) ) < 1 ⇒ x (i) ⋆ = x(i) .
Denote z = ∇f (A f x) and w = Qx1 . We choose R as a sphere centered at

( ζ, ω) = (z, w) max(1, max 1≤i≤I σ • dom G * i ((A f i ) ⊤ z + w (i) ) )
and with radius

r = 2Gap(x, ζ, ω) L f,Q
where

L f,Q = max( max 1≤j≤J L(∇f j ), max 1≤i≤I ρ(Q i,i )) and Gap(x, ζ, ω) = 1 2 x ⊤ Qx + f (A f x) + G(x) + G * (-(A f ) ⊤ ζ -ω) + 1 2 ω⊤ Q † ω + f * ( ζ) .
Note that as ω is a rescaled version of w = Qx, we do not need to know Q † in order to compute ω⊤ Q † ω. It is proved in [START_REF] Ndiaye | Safe Optimization Algorithms for Variable Selection and Hyperparameter Tuning[END_REF] that this set R contains (∇f (A f x ⋆ -b f ), Qx ⋆ ) for any optimal solutions x ⋆ to the primal problem. In the case where G i is a norm and x (i) = 0, these expressions simplify since the σ

• ∂Gi(0) = σ • dom G *
i is nothing else than the dual norm associated to G i .

For the estimation of max

(ζ,ω)∈R σ • ∂Gi(x (i) ) ((A f ) ⊤ i ζ + ω (i)
), we use the fact that the polar of a support function is sublinear and positively homogeneous. Indeed, we have

σ • ∂Gi(x (i) ) ((A f ) ⊤ i ζ + ω (i) ) = sup x : σ ∂G i (x (i) ) (x)≤1 (A f ) ⊤ i ζ + ω (i) , x = sup x : σ c g i D g i ∂g i (D g i x(i) -b g i )≤1 (A f ) ⊤ i ζ + ω (i) , x = 1 c g i D g i σ • ∂gi(D g i x(i) -b g i ) ((A f ) ⊤ i ζ + ω (i) ) = 1 c g i D g i σ • ∂gi(xg i ) ((A f ) ⊤ i ζ + ω (i) ) ≤ 1 c g i D g i σ • ∂gi(xg i ) ((A f ) ⊤ i ζ + ω(i) ) + r sup (u,v): (u,v) =1 1 c g i D g i σ • ∂gi(xg i ) ((A f ) ⊤ i u + v (i) ) .
Here xgi = D g i x(i) -b g i is a point where ∂g i (x gi ) has a nonempty interior. Some care should be taken when ∂G i (x (i) ) is unbounded, so that we first check whether (A f

) ⊤ i ζ + ω (i) ∈ dom σ • ∂gi(xg i )
for all (ζ, ω) ∈ R.

Here also, the novelty lies in the genericity of the implementation.

Numerical validation

Performance

In order to evaluate the performance of the implementation, we compare our implementation with a pure Python coordinate descent solver and code written for specific problems: Scikit learn's Lasso solver and Liblinear's SVM solver. We run the code on an Intel Xeon CPU at 3.07GHz. We can see on Table 6 that our code is hundreds of times faster than the pure Python code. This is due to the compiled nature of our code, that does not suffer from the huge number of iterations required by coordinate descent. On the other hand, our code is about 4 times slower than state-of-the-art coordinate descent implementations designed for a specific problem. We can see it in both examples we chose. This overhead is the price of genericity.

Lasso

We believe that, except for critical applications like Lasso or SVM, a 4 times speed-up does not justify writing a new code from scratch, since a separate piece of code for each problem makes it difficult to maintain and to improve with future algorithmic advances.

Genericity

We tested our algorithm on the following problems:

• Lasso problem min x∈R n 1 2 Ax -b 2 2 + λ x 1 • Binomial logistic regression min x∈R n m i=1 log(1 + exp(b i (Ax) i )) + λ 2 x 2 2
where b i ∈ {-1, 1} for all i. • Sparse binomial logistic regression min

x∈R n m i=1 log(1 + exp(b i (Ax) i )) + λ x 1 • Dual SVM without intercept min x∈R n 1 2α A ⊤ D(b)x 2 2 -e ⊤ x + ι [0,1] n (x)
where ι [0,1] n is the convex indicator function of the set [0, 1] n and encodes the constraint

x ∈ [0, 1] n . • Dual SVM with intercept min x∈R n 1 2α A ⊤ D(b)x 2 2 -e ⊤ x + ι [0,1] (x) + ι {0} (b ⊤ x) • Linearly constrained quadratic program min x∈R n 1 2 (A f ) ⊤ x -b f 2 2 + ι {0} (A h x -b h ) • Linear program min x∈R n c ⊤ x + ι R n + (x) + ι R m -(Ax -b) • TV+ℓ 1 -regularized regression min x∈R n 1 n 2 n 3 1 2 Ax -b 2 2 + α 1 Dx 2,1 + α 2 x 1
where D is the discrete gradient operator and y 2,1 = i,j,k

3 l=1 y 2 i,j,k,l . • Sparse multinomial logistic regression min x∈R n×q m i=1 log q j=1 exp n l=1 A i,l x l,j + n i=1 q j=1 x i,j b i,j + n l=1 q j=1 x 2 l,j
where b i,j ∈ {0, 1} for all i, j. This list demonstrates that the method is able to deal with differentiable functions, separable or nonseparable nondifferentiable functions, as well as use several types of atom function in a single problem.

Benchmarking

In Table 7, we compare the performance of Algorithm 1 with and without screening, Algorithm 2 with and without screening as well as 2 alternative solvers for 3 problems exhibiting various situations:

• Lasso: the nonsmooth function in the Lasso problem is separable;

• the TV-regularized regression problem has a nonsmooth, nonseparable regularizer whose matrix A h is sparse; • the dual SVM with intercept has a single linear nonseparable constraint.

For Algorithm 2, we set the restart with a variable sequence as in [START_REF] Fercoq | Restarting the accelerated coordinate descent method with a rough strong convexity estimate[END_REF]. We did not tune the algorithmic parameters for each instance. We evaluate the precision of a primal-dual pair (x, y) as follows. We define the smoothed gap [START_REF] Tran-Dinh | A smooth primal-dual optimization framework for nonsmooth composite convex minimization[END_REF] as

G β,γ (x, y, ζ, ω) = 1 2 x ⊤ Qx + f (A f x) + G(x) + max y ′ A h x, y ′ -H * (y ′ ) - β 2 y -y ′ 2 +H * (y) + 1 2 ω ⊤ Q † ω + f * (ζ) + max x ′ -(A h ) ⊤ y -(A f ) ⊤ ζ -ω, x ′ -G(x ′ ) - γ 2 x -x ′ 2
and we choose the positive parameters β and γ as

β = dist(A h x, dom H) γ = dist(-(A h ) ⊤ y -(A f ) ⊤ ζ -ω, dom G * ) .
It is shown in [START_REF] Tran-Dinh | A smooth primal-dual optimization framework for nonsmooth composite convex minimization[END_REF] that when G β,γ (x, y, ζ, ω), β and γ are small then the objective value and feasibility gaps are small. For the Lasso problems, we compare our implementations of the coordinate descent method with scikit-learn's coordinate descent [START_REF] Pedregosa | Scikit-learn: Machine learning in python[END_REF] and CVXPY's augmented Lagrangian method [START_REF] Diamond | CVXPY: A Python-embedded modeling language for convex optimization[END_REF] called OSQP. As in Table 6, we have a factor 4 between Alg. 1 without screening and scikitlearn. Acceleration and screening allows us to reduce this gap without sacrificing generality. OSQP is efficient on small problems but is not competitive on larger instances.

For the TV-regularized regression problems, we compare ourself with FISTA where the proximal operator of the total variation is computed inexactly and with LBFGS where the total variation is smoothed with a decreasing smoothing parameter. Those two methods have been implemented for [START_REF] Dohmatob | Benchmarking solvers for TV-l1 least-squares and logistic regression in brain imaging[END_REF]. They manage to solve the problem to an acceptable accuracy in a few hours. As the problem is rather large, we did not run OSQP on it. For our methods, as h is nonzero, we cannot use variable screening with the current theory. Alg. 1 quickly reduces the objective value but fails to get a high precision in a reasonable amount of time. On the other hand, Alg. 2 is the quickest among the four solvers tested here.

The third problem we tried is dual support vector machine with intercept. A very famous solver is libsvm [START_REF] Chang | LIBSVM: A library for support vector machines[END_REF], which implements SMO [START_REF] John | Fast Training of Support Vector Machines using Sequential Minimal Optimization[END_REF], a 2-coordinate descent method that ensures the feasibility of the constraints at each iteration. The conclusions are similar to what we have seen above. The specialized solver remains a natural choice. OSQP can only solve small instances. Alg. 1 has trouble finding high accuracy solutions. Alg. 2 is competitive with respect to the specialized solver. Table 7. Time to reach precision ǫ. scr = with gap safe screening. For SVM, we scale the data so that in each column of the data matrix A, A i,j ∈ [-1, 1].

Conclusion

This paper introduces a generic coordinate descent solver. The technical challenge behind the implementation is the fundamental need for a compiled language in the low-level operations that are partial derivative and proximal operator computations. We solved it using pre-defined atom functions that are combined at run time using a python interface. We show how genericity allows us to decouple algorithm development from a particular application problem. As an example, our software can solve at least 12 types of optimization problems on large instances using primal-dual coordinate descent, momentum acceleration, restart and variable screening.

As future works, apart from keeping the algorithm up to date with the state of the art, we plan to bind our solver with CVXPY in order to simplify further the user experience.

Figure 1 .

 1 Figure 1. Code structure when using no screening and no inertial acceleration

Figure 2 .

 2 Figure 2. Code for the square function atom

Table 1 .

 1 Selected papers for the minimization of smooth functions. Rate: we check whether the paper proves convergence

	Paper	Rate Rand Grad Blck Par Acc	Notable feature
	Seidel '74 [49] Warga '63 [62] Luo & Tseng '92 [34] Leventhal & Lewis '08 [28] Tseng & Yun '09 [60] Nesterov '12 [39] Beck & Tetruashvili '13 [4] Lin & Xiao '13 [33] Lee & Sidford '13 [27] Liu et al '13 [31] Glasmachers & Dogan '13 [22] × × × asymp N N N ✓ Y asymp N ✓ Y ✓ N ✓ Y ✓ Y ✓ Y Y Richtárik & Takáč '16 [46] ✓ Y Allen-Zhu et al '16 [2] ✓ Y Sun et al '17 [55] ✓ Y	× × × × ✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓	× × × × ✓ × ✓ ✓ × × × × ✓ ✓	× × × × × × not eff 1st acc & 1st non-uniform quadratic × strictly convex × rate for weakly convex × quadratic f × line search, proximal operator × × not eff finite time analysis cyclic CD improvements on [39, 47] × not eff ✓ 1st efficient accelerated × ✓ 1st asynchronous × × heuristic sampling × ✓ 1st arbitrary sampling × ✓ non-uniform sampling × ✓ × better asynchrony than [31]
					Blck: the paper considers 1D
	coordinates (×) or blocks of coordinates (✓). Par: Algorithm designed for parallel computing (✓). Acc: no momentum
	(×), accelerated but not efficient in practice (not eff), accelerated algorithm (✓)

(×), an asymptotic rate (asymp) or a finite time iteration complexity (✓). Rand: deterministic (N) or randomized (Y) selection of coordinates. Grad: exact minimization (×) or use of partial derivatives (✓).

Table 2 .

 2 An overview of selected papers proposing and analyzing the iteration complexity of coordinate descent methods for separable non-smooth objectives. Prx: uses a proximal operator to deal with the non-smooth part of the objective. Par: updates several blocks of coordinates in parallel. Acc: uses momentum to obtain an improved rate of convergence. Dual: solves a dual problem but still proves rates in the primal (only relevant for weakly convex duals).

	Paper	Prx Par Acc Dual	Notable feature
	Tseng & Yun '09 [60] S-Shwartz & Tewari '09 [50] Bradley et al '11 [7] Richtárik & Takáč '14 [47] S-Shwartz & Zhang '13 [51] Richtárik & Takáč '15 [48] Takáč et al '13 [56] S-Shwartz & Zhang '14 [53] Yun '14 [68] Fercoq & Richtárik '15 [18] Lin, Lu & Xiao '14 [30] Richtárik & Takáč '16 [45] Fercoq et al '14 [17] Lu & Xiao '15 [33] Li & Lin '18 [29] Fercoq & Qu '18 [16]	✓ ℓ 1 ℓ 1 ✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓	× × ✓ × × ✓ ✓ × × ✓ × ✓ ✓ × × ×	× × × × × × × ✓ × ✓ ✓ × ✓ × ✓ ✓	× × × × ✓ × ✓ ✓ × × × × × × ✓ ×	1st prox, line search, deterministic 1st ℓ 1 -regularized w finite time bound ℓ 1 -regularized parallel 1st proximal with finite time bound 1st dual 1st general parallel 1st dual & parallel acceleration in the primal analysis of cyclic CD 1st proximal and accelerated prox & accelerated on strong conv. 1st distributed distributed computation improved complexity over [47, 39] acceleration in the dual restart for obj with error bound

Table 4 .

 4 Selected papers for the local minimization of nonconvex problems. Conv: study only limit points (×) or convergence of the sequence proved (✓). Smth: can deal with non-convex smooth functions. Nsmth: can deal with non-convex and non-smooth functions.

  the code for the square function atom.

	cdef DOUBLE square(DOUBLE[:] x, DOUBLE[:] buff, int nb_coord, MODE mode,
	DOUBLE prox_param, DOUBLE prox_param2) nogil:
	# Function x -> x**2
	cdef int i
	cdef DOUBLE val = 0.
	if mode == GRAD:
	for i in range(nb_coord):
	buff[i] = 2. * x[i]
	return buff[0]
	elif mode == PROX:
	for i in range(nb_coord):
	buff[i] = x[i] / (1. + 2. * prox_param)
	return buff[0]
	elif mode == PROX_CONJ:
	return prox_conj(square, x, buff, nb_coord, prox_param, prox_param2)
	elif mode == LIPSCHITZ:
	buff[0] = 2.
	return buff[0]
	elif mode == VAL_CONJ:
	return val_conj_not_implemented(square, x, buff, nb_coord)
	else: # mode == VAL
	for i in range(nb_coord):
	val += x[i] * x[i]
	return val

Table 5 .

 5 Convergence speed of the algorithms implemented

Table 6 .

 6 Comparison of our code with a pure Python code and reference implementations for performing 100N coordinate descent iterations for the Lasso problem on the Leukemia dataset with regularization parameter λ = 0.1 (A f ) ⊤ b f ∞, and for 10N coordinate descent iterations for the dual SVM problem on the RCV1 dataset with penalty parameter C = 10.

  >2,000s 267.3s; scr: 114.4s 169.8s; scr: 130.7s TV-regularized regression fMRI[START_REF] Sabrina | The neural basis of loss aversion in decision-making under risk[END_REF][START_REF] Dohmatob | Benchmarking solvers for TV-l1 least-squares and logistic regression in brain imaging[END_REF] (α 1 = α 2 = 5.10 -3 )

	Problem	Dataset	Alternative solver 1	Alternative solver 2	Algorithm 1	Algorithm 2
	Lasso (λ = 0.1 X ⊤ y ∞ )	triazines (J=186, I=60)	Scikit-learn: 0.005s	OSQP: 0.033s 0.107s; scr: 0.101s 0.066s; scr: 0.079s
	ǫ = 10 -6	scm1d (J=9,803, I=280)	Scikit-learn: 24.40s	OSQP: 33.21s	91.97s; scr: 8.73s	24.13; scr: 3.21s
		news20.binary (J=19,996, I=1,355,191)	Scikit-learn: 64.6s OSQP: inexact FISTA: 24,341s LBFGS+homotopy: 6,893s	>25,000s	2,734s
	ǫ = 1	(J=768, I=65,280, L=195,840)
	Dual SVM with intercept ionosphere (J=14, I=351)	libsvm: 0.04s	OSQP: 0.12s	3.23s	0.42s
	ǫ = 10 -3	leukemia (J=7,129, I=72)	libsvm: 0.1s	OSQP: 3.2s	25.5s	0.8s
		madelon (J=500, I=2,600)	libsvm: 50s	OSQP: 37s	3842s	170s
		gisette (J=5,000, I=6,000)	libsvm: 70s	OSQP: 936s	170s	901s
		rcv1 (J=47,236, I=20,242)	libsvm: 195s	OSQP: Memory error	>5000s	63s

We reuse the notation w and z here for the purpose of this section.


Ronan Gaugne

Quentin Petit

Jean-Baptiste Barreau

Valérie Gouranton

Interactive and Immersive Tools for Point Clouds in Archaeology

Figure 1: The immersive environment designed in an archaeological context enables navigation in large and dense point clouds (left) and interactions with the data such as cross-sections manipulation through a tracked tablet (center and right)

Introduction

Throughout the last decade, Cultural Heritage, and in particular archaeology, has benefited from the extensive use of 3D technologies, essentially focusing on reconstitution from modeling or scanning [NPMd12]. Many monuments, buildings and cities were given a new digital existence. Associated with virtual reality and augmented reality technologies, these reconstitution allow us to visualize, simulate and even interact with the past. Virtual reality combined with 3D digitization allows to immerse archaeologists in 1:1 copies of monuments and sites. Most of the time, the graphical environment for these reconstitution is constituted of meshes generated from digitization raw data that can be photos, in the case of A point cloud is the basic raw data obtained when digitizing archaeological sites with laser scans. These data represent a rich and faithful record provided that they have adequate tools to exploit them. They can serve as support for the modelling by snapping edges on existing points [CDD * 17] and their representation allows to distinguish themselves from a hypothetical modelled by geometries [BBP * 14]. They can also be used in BIMs (Building Information Model), without being systematically converted in meshes [ALT15]. More recently, they were use for artistic projects in Cultural Heritage context, as this kind of representation constitutes an aesthetic representation with finesse and fragility of the real world [CMM * 18] or the Palimpsest project http://www.interactivearchitecture.org/ lab-projects/palimpsest .

In this paper, we propose a framework to handle point cloud data in immersive VR, for archaeological purposes.

Related works

Digital point cloud-like data is already widely used in other fields such as medicine, industry, and Geosciences. Each of these areas has its own specificity in their exploitation and use, whether for pathology detection in medicine [AIG * 14], quality control in industry [WMN * 13], or erosion monitoring in Geosciences [START_REF] Lague | Accurate 3d comparison of complex topography with terrestrial laser scanner: Application to the rangitikei canyon (n-z)[END_REF]. It is therefore essential to provide tools adapted to the specific needs related to the field of archaeology. Archaeological material in a context of excavation, restoration or conservation, is generally fragile and any intervention on its structure may damage or even destroy it. This statement fully justifies its systematic digitization, especially endangered material, whether due to climatic, geological or political conditions. Another particular aspect is the uniqueness of the study in archaeological context. It is difficult to provide systematic tools as in medicine and industry where situations of study are essentially repetitive and reproducible. Furthermore, in archaeological context, the case studies are largely based on interpretations and the expression of hypotheses and therefore leave an important part to the archaeologist. It is therefore essential to propose environments in which the archaeologist can project himself/herself and use his/her own perception of a context. This is why the virtual reality that allows to immerse the user in 3D, scale 1, and offer customized interactions, is an interesting answer.

One important challenge is the display of point clouds that can be very large with hundred millions or even billions of points, and cannot benefit from the traditional rendering techniques for 3D meshes. The work presented in [START_REF] Kreylos | Immersive visualization and analysis of lidar data[END_REF] describes an application to visualize LiDAR point cloud data in a CAVE, in a geological context, for quality control, and extraction of survey measurements. The tools developed by the private company Octarina (http://www.octarina.com/ creating-content-experiences-vr-ar/) also propose an interactive environment to manipulate point clouds in a VR headset, in industry context. These two works represent interesting illustrations of the value of immersive point cloud interactions in these two domains. Unfortunately, the proposed interactions are not relevant for archaeological purposes. The work presented in [KKB * ] proposes very interesting interactive tools for collaborative work sessions on a multi-user stereoscopic power wall and a tactile table, applied to very large point clouds obtained from laser scans of a large Cultural Heritage site in Italy. The interactive tools, dedicated to archaeology, are oriented to multiple users sessions for the analysis of the digital data. However, they do not propose a scale 1:1 immersion within the data.

Motivation of the work

Our work contributes with an interactive framework in immersive VR, based on point cloud data, designed specifically to answer archaeological questions. The objective of the study is to explore and illustrate the interest for archaeologists, and more generally Cultural Heritage stakeholders, of manipulating point cloud data in an immersive environment. For this purpose, archaeologists were associated all along the project in order to design, test, and validate different tools to help them in their scientific research.

Our work was performed in the context of a multidisciplinary collaboration between archaeologists and computer scientists started several years ago, with the aim to propose and design new tools and new practices for archaeologists, based on 3D technologies. The work was mainly driven by questions raised by archaeologists during the study of a particular building, an ancient Real Tennis court.

Archaeological context

The building concerned by the study was built at the beginning of the 17th century. It was originally a Real Tennis court, but was transformed into a religious chapel at the end of 17th century, and became property of the army after the French revolution at the end of 18th century [START_REF] Ferrette R | RENNES (35). 10-12 rue Saint-Louis : Les occupations de la parcelle 0234 de l'Antiquité à l'époque moderne ; Le jeu de Paume du Pélican[END_REF]. It is still in elevation and the municipality that owns it since 20th century asked for an archaeological diagnostic in 2014, including a full digitization of its inside and outside (see Figure 2), before launching rehabilitation works. The diagnostic was performed in few weeks, in parallel of the digitization with LiDAR, due to time constraints imposed by the municipality.

In order to pursue their study once the access to the building was no more possible, archaeologists asked computer scientists to propose an application based on the 3D data obtained during the digitization campaign. They needed to understand and document the evolution of the structure of the building. This requires a good representation of the actual volumes, an access to the internal organisation and evolution of the walls and floors, and a highlight of particular details such as the clues of removed structures or recovering of ancient paintings. In order to fulfill these expectations, they requested to be immersed in the digital record of the building, with specific functionalities:

• 1:1 high definition immersion. The first requirement from archaeologists was of course to be able to navigate in the digital copy of the building, at scale 1:1, with a resolution enough faithful and precise to enable further studies. • segmentation of the information delivered by the 3D data.

These segmentation are represented in the form of crosssections, cutting plans, or orthophotos, and depend totally on the different archaeological questions. They are often performed by 3D specialists, often at once on a PC, with no possibility for the archaeologist to superimpose and / or modify them autonomously, due to lack of training on these software, sometimes not very intuitive or expensive. • linking between the 3D environment and 2D documentation, in particular photos taken during the diagnostic. These photos represent an important amount of data and it can be interesting to associate the view of some details highlighted by photos with the gloal 3D data.

Other standard tools such as distance measurement, snapshots of the current point of view, or 3D models integration were also requested by archaeologists.

Contribution of the work

The contribution of our work is a VR application that aims to assist the archaeologists in their process of study, according to their needs. It consists in:

• an immersive VR application dedicated to archaeologists with a display of the dense point cloud at scale 1:1. We chose to work with a point cloud, because it represents the raw data obtained with the laser scan, exempt of any geometrical transformation. However, this kind of 3D data implies some constraints to be used in a VR system, for an archaeological purpose. In particular, the point cloud must be dense enough to allow the archaeologist to observe relevant information. The application was deployed in a CAVE-like structure, to allow the archaeologist to be physically present in the simulation and to have a better perception of the scale (see Figure 3). • a tracked tactile tablet to manage an interactive system of segmentation able to generate cross-sections, cutting plans and orthophotos. We chose this kind of interface as archaeologists often use tactile tablets to record notes in their day to day associated point of view materialized with a sphere. The photo will only be visible when the user's head is co-located with the sphere (see Figure 4, Bottom).

Additional features such as scaling, measuring, dynamic display configuration, and navigation complete the toolbox of the application. Technical details are presented in the next section.

Design and implementation

Point clouds production

The data used in this work were generated using a laser scan technology. The scan of the Real Tennis building was performed in June 2014, with a Faro Focus3D X330. The resulting point cloud is constituted of 781 millions points with both inside and outside the capture of the building (Figure 2). It was sampled to produce two point clouds, one high density point cloud (HDPC) for the close view with 240 millions points and one low density (LDPC) for the distant view with 12 millions points.

Point cloud pre-processing

We used the binary Stanford Polygon format (PLY) for the point cloud. It has the advantage of being a compact and standard format supported by most of the 3D software. PLY files are structured as a list of (x,y,z) points coordinates that can also be enriched by the normal coordinates, (r,g,b) color values, and alpha property. The drawback of this format is that it is not designed to define a structure or hierarchy between the points to help its loading, processing and rendering. In order to optimize the loading and rendering of the point clouds, we used an octree-based octree-based multiresolution approach, using a method similar to [START_REF] Kreylos | Immersive visualization and analysis of lidar data[END_REF]. This kind of structure is common and well-known for point clouds. It allows to perform frustum culling on octree nodes to hide non-visible points using a top down traversal of the octree. The value of the threshold is defined in order to obtain a compromise between the number of draw calls, which depends on the size of the octree, and the frustum culling performance, which depends on the size of the leaves of the octree. The point cloud of the Real Tennis building is constituted of 5744 octree nodes including 4861 octree leaves.

In addition, a subsample is performed on each octree leaf by creating 3D voxel grid over the point cloud leaves (An example of the principle of the algorithm can be found at http://pointclouds.org/documentation/ tutorials/voxel_grid.php). The size of the voxel is user defined and induces grid size and point cloud density. In each voxel only the nearest point to the computed voxel centroid is kept. Each leaf of the octree is stored in a distinct binary file where points coordinates, colors and normals are contiguous in order to match the data sent to the GPU and reduce unnecessary processing overhead. All binaries files of the same resolutions are also packed into one big archive file, basically one file for HDPC and one for LDPC.

Point cloud processing

Architecture

The architecture of the application is designed in order to allow a distributed rendering of the VR part, on a cluster of PCs. Each PC of the cluster runs several instances of the application. The application itself is built using 3D game engine Unity and applicative nodes are synchronized using MiddleVR, a VR plugin for Unity. In addition, the application allows the connection of client applications, which can run on PCs or mobile devices, for specific views such as cross sections. Network clients' connection and synchronization are dynamically managed with Unity Network component.

Point cloud loading

Two different modes can be used to load point cloud. The first mode reads archive files to store all points in RAM at the application startup in order to reduce I/O overhead and lead to better overall performance once point cloud is fully loaded. This mode is particularly useful when working in a distributed architecture, with a clustering mode where the overall performance depends of the worst node performance of the cluster nodes. The drawback of this mode is the limitation according to available RAM. The second mode reads single binary file to only load point clouds that are visible by the camera, in order to reduce RAM usage. This mode enables the handling of very large point clouds. The drawback of this mode is that it induces processing overhead during point cloud navigation.

Runtime display

One issue for the point clouds management is their graphical rendering as they cannot benefit from standard rendering optimization used for 3D surface-based data. We implemented adaptations of frustum culling and lod in order to gain performance on rendering, so we can guarantee an acceptable frame rate for immersive and interactive visualization. Transition between HDPC and LDPC resolutions is done by comparing the distance from the user to the bounding box center to a given threshold that can be dynamically modified (see Figure 5). These two techniques, level of detail and frustum culling, work well together on such data structure: at close range frustum culling discard a large number of points while at long range level of detail heavily reduce the number of points rendered.

Interactions with the point clouds

Various interactions are proposed to the user immersed in the point clouds. First interaction is naturally navigation, in a flying mode, using a tracked wand equipped with a joystick, with four degrees of freedom (translation in the direction pointed by the wand, rotation around the vertical axis).

Interactions with the point clouds are possible in two different ways. The first is through a radial menu used in the application, and controlled with the wand. The second is through a set of icons displayed on a tactile tablet connected to the application. The tablet interface is mainly used for cut plan manipulation.

Point clouds display properties

Points are rendered as a colored square composed of two triangles oriented to face the user like a billboard. Several configurations on the point cloud are proposed to the user. The size of the points can be changed for both LDPC and HDPC. Large point's size leads to a uniform surface representation while a small point size highlights the structure of the current object with transparency rendering. pinpoint the position where the photo was originally taken. The placement of these views is manual and done through five 3D cubic handles on photos, on each corner and in the center (see Figure 4). User can interact with the handles through a tracked controller ray cast; he can grab the photo through the center model and manipulate its position and rotation according to the initial distance between controller and photo. Scale is also modifiable by dragging one of the corner models closer or further of the photo center. During the manipulation, the photo is transparent in order to ease its positioning. Once the position, rotation and scale of the photo match the point cloud they are saved in an external configuration file with the current user head position that represents the new sphere position. When the user puts his/her head inside the sphere, the photo becomes visible, and it disappears when the user quits the position of the point of view. This functionality can also be used to integrate 2D outlines in the 3D environment (see Figure 7).

Scaling

User can modify the scale of the point clouds, either to increase the size of small objects, or to decrease the size of larger environments to get an overview. Scaling is calculated relative to the center of the bounding box of the point cloud using a homothetic transformation directly applied in the rendering for an efficient processing.

Cross section and cutaway plan with tablet

A custom 3D marker placed on a tablet is used to manipulate a cut plane in our VR Cluster. This plane is used as a reference to hide points or make them semi-transparent in order to realize in real time, either a cross section, i.e. a parameterized thickness slice of the point cloud (see Figure 1, center), or a cutaway, i.e. a portion of the point cloud constituted of the points behind the cut plane (see Figure 6, Top). A snap option can be enabled to force the plane to be horizontal or vertical when the plane angle is within a threshold. The tablet shows an orthographic view, orthogonal to the cut plan, of the cross section that can be saved at any time and exported as image. The view on the tablet is dynamically adjusted following its position in the virtual environment (see Figure 1, center). Configuration can be done through the tablet like changing the zoom level, adding a translation or rotation offset, or forcing the color of the points (see Figure 1, right).

Photographic viewpoints

Photographic viewpoints can be added to the point cloud as a complementary visualization tool. A sphere, textured with the associated photo, serves as a landmark to locate this viewpoint and to

Mesh rendering

We can integrate other 3D data inside the VE like CAO 3D model or generated 3D model, that will be combined with the point cloud, to offer different level of information. Figure 8 shows the integra- tion of a graphic modeling of the Real Tennis court at the 18th century that was produced by the archaeologist in charge of the study. In this case, the 3D model was integrated and aligned within the development environment, but it is possible and easy to propose a runtime integration of 3D model, with an interactive positioning of the 3D model similar to the one used for positioning the photos.

Results and discussion

Results

Technical results

The resulting applications were deployed on a large CAVE-like facility constituted of 4 screens, one floor of 9.6m x 2.9m, 6688px x 2000px two lateral screens of 2.9m x 3.1m, 2000px x 2160px, and a main vertical screen of 9.6m x 3.1m, 6688px x 2160px, with images displayed by 14 WQXGA laser projectors, driven by 7 PC, with the following characteristics: 2 Intel Xeon E5-2623v4 2.6 2133 4C, 32GB DDR4-2400, 2 NVIDIA Quadro P6000 24GB.

Each one of the seven PCs of the cluster runs two instances of the application, one for each projector. The tablet used for the interaction is a Microsoft Surface Pro 4, equipped with a tracking geometry built with a 3D printer. Due to the clustering architecture, the overall performance is equal to the performance of the slowest node. The performance on one node depends on the number of points displayed and the size of the points. For the largest point clouds, the Real Tennis building and the Chapel ruin, the performance varies from 10 FPS at full load to 60 FPS, which provides a fluent rendering of the point clouds. The application was also successfully deployed and run on a HTC Vive headset.

Archaeological results

Two archaeologists from outside the project evaluated the application in the large immersive facility. The first one pointed the interest of the 1:1 scale to better perceive the volumes and dimensions of the digitized contexts. He argued that the dimensions are very difficult to evaluate on standard PC screens. He noted that the photographic views allowed to more easily make the connection with the real site. The second one was very enthusiastic about the cutting plan that produces 2D cut views, as this kind of documentation is intensively used when producing excavation reports, especially for buildings. Such 2D images are usually produced by IT technicians from specifications given by archaeologists. In our application the real-time visualization of the 2D cut on the tablet allows the archaeologist to completely master his/her 2D cut and to easily explore alternative 2D cuts during the immersive session. She also pointed out the interest of being immersed at 1:1 scale, in the case of the superposition of the hypothetical modeling and the point cloud, to validate or detect incoherence in the hypothetical model. She also appreciated the high-density area and the point size management that provide a rendering close to a textured mesh. Finally, she noted the interest of such rendering to more easily share her understanding of the building, either in the context of an exchange with another expert, or in the context of scientific mediation, which is a very important activity of archaeology.

Instantiation to three other archaeological contexts

We applied our work to three other archaeological contexts, chosen for their different states of conservation, complexity and archaeological interest. (see Figure 9, Top). The first context was a 13th century ruined chapel, measuring 22m x 13m, with some preserved elevations including the complete pinion with a rose window and arcades. Beyond the classic interest of digital preservation of this unique architectural building [Gou73], the granite weathering due to the temperate oceanic environment [MPS08] accelerates the global erosion of the monument and motivated its digitization. The second context was an 18th-century wreck in a foreshore context, accessible at low tide [RPH * 18]. This context presents interesting characteristics encouraging a digitization project, with a limited size of the site (9m x 3m), a difficult access, as the wreck is under the sea most of the time, an endangered heritage, as the foreshore context is highly destructive. The last context is a cremation urn of Iron Age discovered during an archaeological excavation of a necropolis, among a dozen urns [AL14]. The urns were all CT 5.3.2. Point clouds in VR Some properties of point clouds have a significant impact on their usability in virtual reality. It is important that the point cloud is sufficiently dense and that the points are distributed evenly. The respect of these properties requires a good quality of the stage of digitization. As archaeology teams are increasingly trained in the use of digitization techniques, there is a clear improvement in the quality of point clouds produced, which favors the development of tools for handling this type of data.

Point clouds handling in VR requires to adapt the existing optimization techniques used for 3D meshes. The two techniques we integrated, lod and frustum culling, allowed us to reach comfortable performance for the user with large size point clouds. Another optimization techniques such as occlusion culling could also be adapted to enhance the performances. But this particular technique requires a more complex handling, as the block of points can be transparent depending of the distribution of the points. Additional techniques such as [BWPP04] can be applied in this case.

The use of a CAVE by archaeologists is not always practical or possible, and it is difficult to envisage installing such facility in archaeology laboratories and museums. That is why we have also deployed our point cloud tools on an HMD. This kind of equipment is indeed more affordable and transportable and it is quite possible to include it in the panoply of an archaeological excavation site. However, in this case, the interaction modalities must be adapted, in particular for the manipulation of the cut plane by tablet which will have to be managed completely in the application environment of manipulation of the point cloud since the user does not have real-world visibility when wearing a helmet.

Conclusions

This paper presented a set of tools to manipulate point clouds in a highly immersive and interactive virtual reality environment, in an archaeological context. The tools were developed and evaluated on four different contexts, two buildings and one wreck digitized by LiDAR technology, and one cremation urn digitized by CT scan. The resulting application was successfully deployed in a large CAVE-like facility and was positively evaluated by two archaeologists. The main functionalities of the application are a segmentation tool that allows to produce cross sections and cut-aways, a scaling tool, and an inlay tool that allows to position photographic views inside the point cloud. Further developments are underway to completely adapt the application to a more portable VR headset context. The increasing use of digitization in Cultural Heritage context with the production of more and more digital data justify the development of adapted processing tools, especially around the point clouds which constitute a basic data in this field. Virtual reality, by immersing the user at 1:1 scale in these data is an interesting solution, especially by allowing to propose natural interactions to manipulate them. A video of this work is available at https://vimeo.com/349679831 scanned, and the one presented in this paper revealed two interesting metal objects, a fibula and a knife blade.

The first two contexts represent two endangered sites that are interesting to document for comparative studies and general knowledge, but not enough interesting to justify an expensive restoration. In this case, a laser scan allows to perform a 3D recording of the sites that can be further studied using our framework. For both sites, 2D cuts were performed for archaeological documentation.

For the object-size context, the cremation urn, archaeologists found very interesting the change of scale to get inside the object and obtain an immersive navigation that allows a better understanding of the structure. In particular, they were able to follow the insect galleries between bones fragments and indicates that this could help to better understand the shifting of the fragments and re-associate some of them.

Discussion

Point clouds for archaeology

The digitization of archaeological material makes it possible to meet several objectives. First, a point cloud is a long-term, accurate record of the current state of the site, building, or object. This kind of recording is a major stake of the preservation of our Cultural Heritage. The numerical data obtained also makes it possible to document, study and argue. These data are often collected for a certain purpose, but as pointed out in [CBMM13], it is quite possible to reuse these data later with new objectives. Finally, a cloud of points constitutes an aesthetic data format that can be directly used in the scientific mediation framework, provided that tools are available to handle it.

The tools proposed as part of the work presented in this article have been designed to respond to actual needs of archaeologists. A work session around the point cloud of the chapel made it possible to easily and quickly generate cutting plans allowing documentation of the monument subject to significant erosion. In the same way, the point cloud of the wreck is a working support that can be used while the real wreck, buried under 80cm of sand and covered by the sea most of the time, is very difficult to access. Cutting plans were also generated for this wreck and archaeologists noted that the generation of such cuttings was faster and more accurate than the usual hand-made drawings. The views of the urn and the change of scale have allowed archaeologists to have a unique and innovative view of its internal structure. The Real Tennis building was the richest source of data for the design and evaluation of point cloud manipulation tools. Several work sessions took place in the immersive environment and resulted in these first results. The fact that both the external and internal parts of the building were scanned, and that a 3D reconstitution of the original building was modeled by the archaeologist provided a rich amount of data source that is not always available. However, 3D digitization during archaeological excavations is more and more systematic and multi-source. It is thus important to be able to combine different kinds of data, either produced from different technologies such as LiDAR and CT scan, or of different nature, as in the combination of point cloud data, photographic views and surface meshes.

  Different approaches have been proposed to display large point clouds by structuring the data [GEM * 13], [SW11] or by applying rendering optimization techniques [PGK02]. These rendering tools are designed for visualization and manipulation on a single PC. They generally require software skills to use and can create ambiguities regarding the architectural dimensions. Recent works explored the rendering of large point clouds with immersive first-person view [PMM * 17] or in immersive display structures [BSN14], [TBP15]. However, these works focus on the rendering and do not propose interactive tools or limited ones.

Figure 2 :

 2 Figure 2: Point cloud of the Real Tennis building.

Figure 3 :

 3 Figure 3: Immersive navigation in the point cloud in a large CAVE-like structure.

Figure 4 :

 4 Figure 4: Top: Positioning of a photo within the VR application. Bottom: A floating sphere indicates a point of view of a photo in the application.

Figure 6 :

 6 Figure 6: Cutaway manipulation.

Figure 7 :

 7 Figure 7: Integration of an outline inside the virtual environment.

Figure 5 :

 5 Figure 5: Top: Low density point cloud. Middle: High density point cloud close to the user, Low density at longer distance. Bottom: Full High Density point cloud.

Figure 8 :

 8 Figure 8: Superposition of the graphic modeling of the Real Tennis building.

Figure 9 :

 9 Figure 9: Top: Cremation urn. Bottom, left: Wreck. Bottom, right: 13th century Chapel.

Acknowledgements

We would especially like to thank Yann Bernard and Florian Cousseau for the 3D scans done on the sites, Elen Esnault and Théophane Nicolas for their archaeological feedback.