Luca Castelli Aleardi

Alexandre Nolin
email: alexandre.nolin@polytechnique.edu

Maks Ovsjanikov

Efficient and practical tree preconditioning for solving Laplacian systems

We consider the problem of designing efficient iterative methods for solving linear systems. In its full generality, this is one of the oldest problems in numerical analysis with a tremendous number of practical applications. In this paper, we focus on a particular type of linear systems, associated with Laplacian matrices of undirected graphs, and study a class of iterative methods for which it is possible to speed up the convergence through the combinatorial preconditioning. In particular, we consider a class of preconditioners, known as tree preconditioners, introduced by Vaidya, that have been shown to lead to asymptotic speed-up in certain cases. Rather than trying to improve the structure of the trees used in preconditioning, we propose a very simple modification to the basic tree preconditioner, which can significantly improve the performance of the iterative linear solvers in practice. We show that our modification leads to better conditioning for some special graph structures, and provide extensive experimental evidence for the drastic decrease in the complexity of the preconditioned conjugate gradient method for several classes of graphs, including 3D meshes and complex networks.

Introduction

Solving general linear systems of equations is one of the oldest and best studied areas of numerical analysis, with an abundance of both exact and approximate solutions of varying efficiency (see e.g., [START_REF] Golub | Matrix Computations[END_REF]). In this paper, we focus on iterative methods for solving a particular type of linear systems, associated with Laplacian matrices of undirected graphs. These linear systems arise in a variety of applications, which are related to solving the Poisson equation on discretized domains, including physical (e.g. fluid) simulation, complex system analysis, geometry processing and computer graphics [START_REF] Koutis | A fast solver for a class of linear systems[END_REF][START_REF] Krishnan | Efficient preconditioning of laplacian matrices for computer graphics[END_REF], among many others. One class of techniques, which is especially useful in solving large systems of equations with Laplacians matrices of certain (sparse) graphs is the conjugate gradient method. This method can be classified as an iterative approach, since it only requires the ability to compute matrix-vector products and provides progressively better estimates of the final solution. It is also known to terminate in finite number of steps depending on the quality of the initial guess and the condition number of the matrix in question [START_REF] Vishnoi | Lx = b[END_REF]. The convergence speed of the conjugate gradient method can further be improved significantly using preconditioning, which aims to approximate the given matrix A by another matrix B (a preconditioner), whose inverse can be readily computed. The quality of the improvement provided by the preconditioner is directly related to the difference between B -1 A and identity.

In recent years a particular class of preconditioners has been proposed for solving the Poisson equation on undirected graphs, by using the so-called combinatorial (or geometric) preconditioning [START_REF] Beauwens | Lower eigenvalue bounds for pencils of matrices[END_REF][START_REF] Vaidya | Solving linear equations with symmetric diagonally dominant matrices by constructing good preconditioners[END_REF]. The main idea, first proposed by Vaidya, is to approximate a given graph by its subgraph, on which the system of equations can be solved easily. The canonical example of this type of preconditioner is a spanning tree of the graph. Since the Poisson equation can be easily solved in linear time if the graph is a tree, the main idea in Vaidya's approach is to use the Laplacian of the spanning tree as a preconditioner to improve the convergence of iterative methods, such as the conjugate gradient. This basic framework has been extended significantly to both obtain near-optimal trees that can approximate arbitrary graphs, and to use a recursive approach in which a graph can be approximated by a progressively more accurate subgraphs, which can lead to very significant asymptotic speed-up in solving linear systems on general graphs [START_REF] Spielman | Nearly-linear time algorithms for graph partitioning, graph sparsification, and solving linear systems[END_REF]. While the theoretical framework for combinatorial preconditioners has been developed and in some ways settled, with a few notable exceptions, the practical implementations of these ideas are still largely lacking. This can be attributed, in part, to the highly complex nature of the algorithms for obtaining the optimal preconditioners, with potentially very large constants in the asymptotic analysis on the one hand [START_REF] Spielman | Nearly-linear time algorithms for graph partitioning, graph sparsification, and solving linear systems[END_REF][START_REF] Spielman | Nearly linear time algorithms for preconditioning and solving symmetric, diagonally dominant linear systems[END_REF], and the relatively little improvement provided by the basic tree preconditioner on the other hand [START_REF] Chen | Vaidya's preconditioners: implementation and experimental study[END_REF]. As a result, despite the theoretical appeal and the near-optimality in the asymptotic sense of the resulting algorithms [START_REF] Cohen | Solving SDD linear systems in nearly mlog 1/2 n time[END_REF], the practitioners have not yet fully benefited from the potential practical improvements provided by the combinatorial preconditioners.

Contribution In this paper, we concentrate on the basic setting of Vaidya's preconditioners where the Laplacian matrix of a single spanning tree is used as a preconditioner for the conjugate gradient method. Indeed, by extending the experiments of Chen et al. [START_REF] Chen | Vaidya's preconditioners: implementation and experimental study[END_REF] to a variety of graphs and large networks, we show empirically that in most cases the improvement given by a single preconditioner is either minor or even non-existent compared to the baseline conjugate gradient approach. In this context, our main contribution is to propose a very simple modification to Vaidya's tree preconditioner, which provides significant practical improvements, with minimal implementation effort. Our modification can be seen as a combination of a basic Jacobi (diagonal) preconditioner with a combinatorial (tree) one. Despite its extreme simplicity, we show that on a set of important special cases, our approach can lead to a significant decrease in the condition number of the resulting system, compared to the baseline combinatorial preconditioner. Perhaps more importantly, however, we also show via extensive experimentation, that our modification can also lead to very significant practical speedup in the convergence of the conjugate gradient method compared to both the Jacobi and tree preconditioners for a large number of classes of graphs and different target functions. Our approach is not meant to provide a preconditioner structurally very different to existing ones, or to improve on their asymptotic complexity. Rather, by showing that a simple modification of the tree preconditioner can potentially lead to significant practical improvements, we hope to demonstrate the usefulness of such preconditioners and to help eventually bridge theory and practice in this field.

Related works A tremendous amount of progress has been done in solving linear systems associated to symmetric diagonally dominant matrices in the recent past. Classical iterations, as described in [START_REF] Golub | Matrix Computations[END_REF], were very sensitive to systems of poor condition number, and until quite recently efficient preconditioning was mostly a mat-ter of heuristics. A major step was done by Spielman and Teng [START_REF] Spielman | Nearly-linear time algorithms for graph partitioning, graph sparsification, and solving linear systems[END_REF], presenting the first nearly linear algorithm. Their results have since then been split into three papers [START_REF] Spielman | Nearly linear time algorithms for preconditioning and solving symmetric, diagonally dominant linear systems[END_REF][START_REF] Spielman | A local clustering algorithm for massive graphs and its application to nearly linear time graph partitioning[END_REF][START_REF] Spielman | Spectral sparsification of graphs[END_REF]. The feat was made possible by the introduction and refinement of ideas such as spectral sparsification, ultra-sparsifiers and algorithms for constructing low-stretch spanning trees, which they cleverly assembled together to build a recursively preconditioned iterative solver. The general idea of preconditioning recursively is the basis of today's best solvers [START_REF] Cohen | Solving SDD linear systems in nearly mlog 1/2 n time[END_REF], and its individual parts have been separately improved over the years. For instance, low-stretch spanning trees, first introduced with no link to preconditioning [START_REF] Alon | A graph-theoretic game and its application to the k-server problem[END_REF], have seen an improvement over the years. Their use as preconditioners was suggested in [START_REF] Vaidya | Solving linear equations with symmetric diagonally dominant matrices by constructing good preconditioners[END_REF], in the continuity of the ideas of support theory and combinatorial preconditioning (see [START_REF] Beauwens | Lower eigenvalue bounds for pencils of matrices[END_REF][START_REF] Boman | Support theory for preconditioning[END_REF][START_REF] Bern | Support-graph preconditioners[END_REF] for early work and formalizations of this theory). Interested readers can read the progression of the stretch in [START_REF] Abraham | Using petal-decompositions to build a low stretch spanning tree[END_REF], where an algorithm for computing trees of total stretch O(m log n log log n) in time O(m log n log log n) is given, which is the best we know of today. If no better bound has been found since then, recent works introduced a generalization of stretch [START_REF] Cohen | Stretching stretch[END_REF], creating new possibilities of optimization. Spectral sparsification has seen similar improvements, however its progression is less linear than that of spanning trees. The existence of better sparsifiers has been shown in [START_REF] Batson | Twice-ramanujan sparsifiers[END_REF] and fast construction algorithms have been proposed by recent works [START_REF] Kolla | Subgraph sparsification and nearly optimal ultrasparsifiers[END_REF][START_REF] Koutis | Faster spectral sparsification and numerical algorithms for SDD matrices[END_REF] (we refer to [START_REF] Batson | Spectral sparsification of graphs: theory and algorithms[END_REF] for a detailed presentation of the recent advances on this topic). Other works have departed from the initial recursive framework in an attempt to produce simpler solvers, with [START_REF] Nesterov | Efficiency of coordinate descent methods on huge-scale optimization problems[END_REF] and recently improved in [START_REF] Kelner | A simple, combinatorial algorithm for solving sdd systems in nearly-linear time[END_REF][START_REF] Lee | Efficient accelerated coordinate descent methods and faster algorithms for solving linear systems[END_REF]. While their modified gradient descent method that actualizes small parts of the vector at any time is very different from the other solvers, they use similar tools of graph decomposition and approximation (such as low-stretch spanning trees).

Preliminaries and background

Graph Laplacian

Throughout the paper, we consider simple, undirected, non-weighted graphs G = (V, E) with #V = n and #E = m, and d(i) = #{j, (i, j) ∈ E} the degrees of the vertices. The unweighted (and un-normalized) Laplacian matrix L G is given via its relation to diagonal degree matrix D G and the adjacency matrix A G :

D G = d(i) if i = j 0 o/w , A G = 1 if (i, j) ∈ E 0 o/w , L G = D G -A G
It can be readily verified that the Laplacian matrix L G is symmetric, diagonally dominant, and as such only has non-negative eigenvalues. Indeed, it can be readily seen that the number of connected components of G equals the dimension of the null space of L G . Throughout our paper we assume to be working with a connected graph G. In this case the eigenvalues of L G are given as 0 = λ 1 < λ 2 ≤ . . . ≤ λ n . To simplify the exposition, in this paper we assume uniform weights on all edges of the graph. Nevertheless, most of the material (excluding the analysis of our modification of the tree preconditioner) can be adapted to the case of positively-weighted edges.

Solving Linear Systems

The canonical problem that we consider is to solve linear system of equations of the form Ax = b, where, in our case A = L G for some known vector b. Depending on the domain, a problem of this form may also be known as solving the discrete Poisson equation. In general, although the number n of vertices in the graph can be very large, the matrix L G is typically sparse, which can make direct solvers inefficient or even not applicable, since a full n 2 set of variables can easily exceed the available memory. Instead, iterative solvers have been used to solve this problem, and most notably the Conjugate Gradient (CG) method, which is especially useful in cases with limited memory, since it requires only matrix-vector product computations. This method is applicable to symmetric positive (semi)-definite systems, and computes successive approximations of x by taking a step in a direction conjugate to previously taken steps, where conjugacy between two vectors x 1 , x 2 is defined as x T

1 Ax 2 = 0 (please see Chap. 11 in [START_REF] Golub | Matrix Computations[END_REF] for a full discussion of this method). Due to the classical and well-studied nature of the CG method, we do not describe it in detail, but instead refer the reader to the multiple excellent expositions on this topic (for completeness, the pseudo-code for the method is provided in Figure 1). It is well-known that in the absence of rounding errors, the conjugate gradient method will converge in at most n iterations in the worst case. A more relevant bound, however, can be given by using the condition number κ(A) of the matrix A, given by the ratio of its largest and smallest non-zero eigenvalues. In particular, after t iterations, the error of the algorithm is bounded by:

||x (t) -x|| A ≤ 2 1 - 2 λ n /λ 2 + 1 t ||x|| A (1)
Note that while the Conjugate Gradient method is best suited for positive definite matrices, it can also be easily adapted to positive semi-definite systems, such as the ones including the graph Laplacian. One simply has to make sure that the right hand side of the equation lies in the span of the matrix. For us, this means that the vector b has to sum to zero. Let us also stress that λ 1 in the definition of the condition number is the first non-zero eigenvalue. This will become particularly important when we define and analyze the properties of the preconditioned conjugate gradient.

Preconditioning

Since the condition number gives a simple bound on the efficiency of iterative solvers, and of the conjugate gradient method in particular, it is natural to try to introduce linear systems equivalent to the original one, but with a lower condition number, and therefore better convergence properties. This process, called preconditioning, requires a non-singular matrix M , such that M -1 ≈ A -1 . Then, instead of solving Ax = b directly, we solve :

C -1 AC -1 x = C -1 b (2)
where C 2 = M , and x is found by solving Cx = x. Ideally, the preconditioner M should be a positive (semi)-definite matrix, such that the condition number of M -1 A is significantly smaller than that of A itself. The design of optimal preconditioners Initialization

r (0) = b -Ax (0) p (0) = r (0) Iteration                α (k) = ||r (k) || 2 ||p (k) || 2 A x (k+1) = x (k) + α (k) p (k) r (k+1) = r (k) -α (k) Ap (k) β (k) = ||r (k+1) || 2 ||r (k) || 2 p (k+1) = r (k+1) + β (k) p (k) Initialization    r (0) = b -Ax (0) z (0) = M -1 r (0) p (0) = z (0) Iteration                    α (k) = r (k) ,z (k) ||p (k) || 2 A x (k+1) = x (k) + α (k) p (k) r (k+1) = r (k) -α (k) Ap (k) z (k+1) = M -1 r (k+1) β (k) = z (k+1) ,r (k+1) z (k) ,r (k) p (k+1) = z (k+1) + β (k) p (k)
Fig. 1. Conjugate Gradient and Preconditioned Conjugate Gradient : this pseudocode shows the general idea of orthogonalization in both algorithms as well as how the preconditioning takes place in PCG.

typically involves a trade-off: on the one hand, M -1 A should be as close to identity as possible. On the other hand, it should be possible to solve a linear system of the form M x = b very quickly, since it has to be done at every CG iteration. An example of potentially useful preconditioning is the Jacobi Preconditioner for diagonally dominant systems. This consists in taking the matrix D = (δ ij A (i,j) ) (i,j) , i.e the diagonal of the original matrix, as preconditioner. This is both very easy to compute, and solving Dx = b takes an optimal O(n) operations. When both A and M are symmetric positive definite, then solving Eq. 2 can be done without explicitly computing the matrix C, by modifying the steps taken during the iterations of the Conjugate Gradient method. This results in the Preconditioned Conjugate Gradient for which we provide the pseudo-code in Figure 1.

Note that for positive semi-definite systems, one has to use the pseudo-inverse in Eq. 2 above, and make sure that the kernel of M is contained in the kernel of A, for otherwise the system M x = b, may not have a solution. In most cases, when the preconditioning is applied to positive semi-definite systems the kernels of M and A coincide, although the framework can also be applied in a more general case.

Spanning Trees as Preconditioners for Graphs

While both the basic and the preconditioned Conjugate Gradient method can be applied for any positive (semi)-definite linear system, the design of preconditioners can highly benefit from the knowledge of the structure of the matrix in question. As mentioned above, in this paper, we concentrate on the linear systems arising from the Laplacian matrices of undirected graphs. In this case, a particularly promising idea, first proposed by Vaidya and then extended significantly in the recent years, is to use the Laplacian matrix of a subgraph as a preconditioner to the original system. Note that, if the subgraph is connected over the same set of nodes as the original graph, then the kernels of the Laplacian matrices both have dimension 1, and they contain the constant vector 1 n and all the vectors parallel to it, making the use of preconditioning directly applicable. A particularly appealing candidate for a subgraph to be used as a preconditioner is a spanning tree T of the graph G. This is because if L T is the Laplacian matrix of the tree T , then the problem of type L T x = b can be solved very efficiently, in time O(n), with two tree traversals. This makes spanning trees good candidates for preconditioning, because their use keeps the cost per PCG iteration in O(m). It can be shown [START_REF] Vishnoi | Lx = b[END_REF] that for a spanning tree T of G, κ(L † T L G ) ≤ stretch T (G), where the stretch is defined as the sum of the distances in the tree between any two vertices connected by an edge in G. Together with Eq. 1 this can be used to establish the convergence of the preconditioned Conjugate Gradient for Laplacian matrices.

We note briefly that better bounds can be proved by also looking at the distribution of the eigenvalues. A proof using a lower bound for all eigenvalues and an upper bound on the number of eigenvalues above a certain threshold yields that PCG computes an -approximation in O (stretch T (G))

1/3 log(1/ ) iterations (see Lemma 17.2 in [START_REF] Vishnoi | Lx = b[END_REF]). In the past several years, this basic framework for solving linear systems with Laplacian matrices has been extended significantly, with two major research directions: finding trees that can optimize the stretch with respect to arbitrary large graphs [START_REF] Abraham | Using petal-decompositions to build a low stretch spanning tree[END_REF], and changing this basic framework to use a more sophisticated hierarchical graph approximation scheme in which preconditioners themselves can be solved via iterative (and possibly recursive) schemes [START_REF] Koutis | A fast solver for a class of linear systems[END_REF]. Unfortunately, both of these directions lead to highly complex algorithms (practical performances have been evaluated only very recently [START_REF] Hoske | Is nearly-linear the same in theory and practice? a case study with a combinatorial laplacian solver[END_REF]). Rather than trying to improve either of these two directions, our goal is to show that a simple modification to the tree preconditioner can significantly improve the performance of the iterative solver both in theory (for some restricted case) and in practice (over a large number of empirical tests). Our goal, therefore, is to provide a practical and efficient extension of the basic tree preconditioning framework.

Contribution: enhancing tree-based preconditioners

As mentioned in the introduction, our main goal is to show that a simple modification of the combinatorial (tree) preconditioner can have positive impact on the practical performance of the preconditioned conjugate gradient. Indeed, as has been noted by Chen et al. [START_REF] Chen | Vaidya's preconditioners: implementation and experimental study[END_REF] and we confirm in Section 4, the basic version of Vaidya's approach rarely results in significant practical benefits for PCG. The critical remark behind our work is that it is possible to add positive terms to the diagonal of the preconditioning matrix L T without changing its combinatorial structure that enables the fast resolution of associated linear systems. Thus, we introduce the matrix H T = L T + D G -D T = D G -A T . Note that the matrix H T has the same diagonal as the Laplacian L G , but the same sparsity structure as the Laplacian of the subgraph T . Therefore, solving a linear system of equations of the type H T x = b can still be done in exactly the same time as solving L T x = b. Nevertheless, as we show below theoretically (on some restricted cases) and empirically on a large number of different graphs and linear systems, this simple modification can significantly boost the performance of the PCG method.

Before proceeding to the analysis of our modification to Vaidya's preconditioner, we first note that unless T = G, the matrix H T will be full-rank, unlike the L G which has a kernel consisting of vectors parallel to the constant vector 1 n . While in practice, this does not change the method shown in Figure 1, we note that the analysis needs to be adapted slightly. Namely, since we are operating in the space orthogonal to the constant vector 1 n , we need to make sure that the condition number of the preconditioned system is calculated correctly. For this, the following Lemma, which is readily verified, is useful: Lemma 1. The eigenvalues of the generalized eigenvalue system L G x = λH T x are the same as those of the system L G x = λP H T x, where P = (I n -1 n 1 n 1 T n ) is the projection onto the space of vectors orthogonal to the constant vector.

Therefore, computing the condition number κ(L G , H T ) of the preconditioned system can be done by considering the ratio of the largest to smallest non-zero eigenvalues of the matrix H -1 T L G . Equivalently, one can consider the smallest and largest value c such that x T (L G -cH T )x ≥ 0 for all x, such that x T H T x c = 0.

To motivate the use of our preconditioner as well as to provide some intuition on its behavior we proceed in two stages. First, we show some bounds on the condition number for special graphs, and second, we demonstrate empirically that for a very wide range of large scale graphs and linear systems our approach can significantly outperform other baseline preconditioners (Section 4).

Some bounds for special graphs

In this section we provide bounds on the condition number of the preconditioned system for Laplacians of special graphs and show that we can obtain significant theoretical improvement over Vaidya's preconditioners in some important special cases.

The complete graph We first evaluate the performance of our modified tree preconditioner for the class of complete graphs. Let us consider G = K n , the complete graph on n vertices and let T be a star spanning tree, consisting of one root vertex of degree n -1 which is adjacent to all remaining n -1 vertices. Lemma 2. Given the complete graph G and the tree T described above, then for any

n > 2 we have κ(L G , H T ) = n n-1 < κ(L G , L T ) = n. Note, in particular that κ(L G , H T ) → 1 whereas κ(L G , L T ) grows with n.
The ring graph Another important example is the case of the cycle (ring) graph with n vertices. Here, the tree T differs from G by a single edge. In this case we have the following result: Lemma 3. If G is a cycle and T is a spanning tree of G, then κ(L G , H T ) < 2, while κ(L G , L T ) = n for any n.

Note that again, the system preconditioned with H T remains well-conditioned for all n, unlike the system preconditioned by the tree itself, which has an unbounded condition number. Indeed, a strictly more general result holds: Lemma 4. Let G be any graph and T be a tree on G, such that the edge-complement

T c of T in G is a star. Then: κ(L G , H T ) ≤ 2.
Note that this lemma generalizes the previous one since the complement of the tree in the ring graph is a single edge. The wheel graph Our final example is the wheel graph, consisting of a cycle with n -1 vertices that are all connected to a central vertex s, which does not belong to the cycle. In this case, let T be the star graph centered around s. Lemma 5. Given the graph G and the spanning tree T described above then, for any

n odd, κ(L G , H T ) < κ(L G , L T ) = 5.
This example is instructive since the wheel graph can be considered to be a simple case of a triangle mesh, a class of graphs for which we show empirically a significant improvement over Vaidya's preconditioners in the next section. Fig. 3. We compute condition numbers for (unweighted) grids endowed with different spanning trees: blue and red edges correspond to trees with high and low stretch factors respectively. Our precondition matrix H allows to drastically decrease the condition number in both cases, when compared to standard tree preconditioning (dashed lines).

A counterexample We also note that there exist graphs for which the condition number of our system is worse than that of the unmodified approach. The simplest example of such a graph is a path-graph with 6 nodes, with additional edges between nodes (1, 3) and [START_REF] Batson | Spectral sparsification of graphs: theory and algorithms[END_REF][START_REF] Bern | Support-graph preconditioners[END_REF], and where the tree is the path. In this case, it can be easily seen that κ(L G , L T ) = 3 < κ(L G , H T ). Nevertheless our experiments suggest that such cases are rare, and seem to occur when the graph G is very close to the tree T . We leave the full characterization of such cases as interesting future work.

Experimental results

We provide experimental evaluations of the performances of our preconditioner against CG (conjugate gradient with no preconditioning), the diagonal preconditioner JPCG (Jacobi preconditioned conjugate gradient) and TPCG (tree-based Vaidya's preconditioned conjugate gradient). As our preconditioner is a combination of the tree-based and diagonal approaches, we denote it by JTPCG.

Data sets. We run our experiments on a wide collection of graphs including triangle meshes (obtained from the AIM@SHAPE Shape repository), 2D regular grids, and complex networks (from the Stanford Large Network Dataset Collection). We also consider random planar triangulations, generated by the uniform random sampler by Poulalhon and Schaeffer [START_REF] Poulalhon | Optimal coding and sampling of triangulations[END_REF], as well as graphs randomly generated according to the small-world and preferential attachment models.

Evaluating the condition number

Regular grids. Our first experiments concern the evaluation of the condition numbers for regular grids, for which we know how to construct spanning trees of high and low stretch factors. It is not difficult to see that the total stretch of the blue tree T h in Fig. 3 is Θ(n √ n) (observe that a vertical edge (u i , v i ) ∈ G \ T belonging to the i-th column contributes with a stretch of Θ(i), where i ranges from 1 to √ n). The red edges in Fig. 3 define a spanning tree T l having a low stretch factor, which can be evaluated to be O(n log n) using an inductive argument (we refer to [START_REF] Koutis | A fast solver for a class of linear systems[END_REF] for more details). These bounds reflect the numerical evaluation of the condition numbers for both trees T h and T l (plotted as dashed lines in Fig. 3). Experimental evaluations show that our Jacobi-tree preconditioner allows to drastically decrease the condition numbers for both trees T h and T l . More interestingly, using H T instead of L T we obtain new bounds which are extremely close, despite the very different performances of the corresponding spanning trees. Not surprisingly, this behavior does not only concern regular grids, but it is common to a wide class of graph laplacians, as suggested by experimental evidence provided in next sections. Fig. 4. We compute condition numbers for several 3D surface meshes, comparing tree preconditioning and Jacobi-tree preconditioning. Meshes are endowed with both minimum (blue) and maximum (red) spanning trees (weights correspond to Euclidean edge length).

Mesh graphs and complex networks. We compute numerical evaluations of condition numbers for laplacians corresponding to several 3D meshes, of different sizes and topology (refer to Fig. 4). We test and compare our preconditioner against the CG method (without preconditioning) and tree preconditioning, using as test trees both minimum and maximum spanning trees. We consider min spanning trees because their performances are in general pretty worse than those of maximum spanning trees: weights are computed according to the euclidean edge length of the 3D embedding (in the case of unweighted graphs, in order to compute min and max spanning tree, we reweight edges according to a vertex degree driven strategy). Let us firstly remark that our experiments confirm the intuition of Vaidya's seminal work: maximum spanning trees perform in general better as preconditioners than other trees. Once again, our preconditioner is able to get condition numbers which are significantly lower than the ones obtained with the simple preconditioner L T . As already mentioned, when comparing the behaviors of maximum and minimum spanning trees the condition numbers obtained via H T collapse getting very close. As observed in next section when comparing linear solvers, this phenomenon occurs for all tested meshes and graphs and results in a significant improvement of the performance of iterative solvers.

Counting iterations: comparison of iterative linear solvers

In this section we provide experimental evidence of the improvement achieved by our JTPCG preconditioner. We test it against other linear solvers (CG, JPCG, and TPCG) on a large set of surface meshes and random graphs. In order to obtain a fair comparison of performances, we measure the convergence rates of linear solvers for Lx = b counting the total number of iterations required to achieve a given error: as metrics we use the standard relative residual error. Fluid simulation We use iterative solvers as core linear solvers for simulating fluid diffusion on regular 2D grids of different sizes. We count the number of iterations Fig. 5. Fluid simulation: we compare JTPCG against CG and TPCG (JPCG is not shown, as its performance is not significantly different from that of CG). We plot the proportion of the amount of iterations required by different methods for solving 100 (resp. 200) linear systems with a fixed precision of 1e -5. For instance, JTPCG (colored curves) takes between 51 and 127 iterations per linear system on a grid of size 4096, while CG (black curve) requires between 75 and 184 iterations.

required by different solvers at each time step (we use fixed precision 1e -5). As shown by the plots in Fig. 5, JTPCG is able to drastically decrease the number of iterations, using both the high and low stretch factor spanning trees (red and blue curves). Observe that tree-based preconditioner perform pretty well (even without diagonal modification) when combined with low stretch factors (red curves in Fig. 5). Table 1. Solving linear systems: we compare the JTPCG against the classical CG method, the JPCG and TPCG preconditioners. We count total number of iterations required to achieve fixed precision 1e -7.

Solving mesh laplacians. The results reported in Table 1 concern the resolution of linear systems of the form Ax = b, where the right term is a random vector b orthogonal to the constant vector. We use the same starting vector as initial guess for all linear solvers (tests are repeated several times, in order to take into account the dependency of the convergence speed on the initial guess). As confirmed by the numerical results reported for the few meshes in Table 1, our preconditioner always performs better (some times just slightly better) than other solvers (this behavior has been confirmed for all tested meshes).

Iterative eigensolvers and spectral embeddings. Spectral methods proved their relevance in various and distinct application domains, ranging from graph drawing and data visualization to complex networks analysis and geometry processing. (for more comprehensive readings on these topics we refer to [START_REF] Levy | Spectral geometry processing[END_REF][START_REF] Koren | Drawing graphs by eigenvectors: Theory and practice[END_REF][START_REF] Luxburg | A tutorial on spectral clustering[END_REF]). We also have integrated our preconditioner as the core of an iterative eigensolver (we have implemented a hybrid version of the inverse power iteration method). We evaluate the experimental performances by computing the smallest non-trivial eigenvalues of the Laplacian matrix (this is a base fundamental step in problems such as spectral drawing and

  , LT ) κ(LG, HT ) κ(LG, LT ) κ(LG, HT ) (max tree) (max tree) (min tree) (min tree)

Fig. 2 .

 2 Fig. 2. Condition numbers for unweighted graphs: we consider a few example of complex networks and random triangulations. Left pictures show the metabolic system of the C. elegans worm and a random planar triangulation (picture provided by N. Curien).

Fig. 7 .

 7 Fig. 7. Spectral clustering and complex networks: the picture above shows a into five sets of a social network (facebook, 4k nodes) that we obtained applying the K-means algorithm to the spectral embedding of the graph.

spectral clustering). Tables in Fig. 6 and 7 report the total number of iterations performed by the linear solvers required to compute the first three eigenvalues for 3D surface meshes and complex networks. Fig. 6. The pictures above show the triceratops mesh together with its 3D spectral embedding. This table reports the total number of iterations performed by the iterative linear solvers during the inverse power iteration (our tests are run with fixed precision 1e -5).

Concluding remarks

The main goal of this work is to propose a very simple and efficient modification of combinatorial preconditioners, for which we show theoretical (in some restricted cases) and empirical (on a large collection of various graphs) improvement in performance over the original tree preconditioning. The best performances concern mesh-like structures, for which JTPCG achieves very interesting convergence speed in all test applications. In the case of complex networks, it is more difficult to reveal any weakness or strength of a particular linear solver: the behavior significantly depends on the structural properties of tested graphs (for example, while TPCG performs poorly on meshes, it achieves very interesting convergence speed on the Power Grid network). We hope that this work will motivate further theoretical investigations about the practical relevance of our method revealed by our experimental evaluations.

Appendix

Proof of Lemma 2: To simplify the computation, we consider the basis for vectors of size n orthogonal constant vector given by two sets: {B1} consisting of all vectors with value 0 at the center of the star tree, and the vector b2, whose value at the center is

and at all other vertices is -

. Note that if G is a complete graph then LGx = nx for any x orthogonal to the constant vector. Moreover, for the star tree LT x = x if x ∈ {B1} and LT b2 = nb2. Therefore, κ(LG, LT ) = n. It is also easy to see that HT x = (n -1)x if x ∈ {B1} and HT b2 = nb2 + cHxconst. Thus, κ(LG, HT ) = n n-1 . Proof of Lemmas 3 and 4: To see that κ(LG, LT ) = n if G is a cycle and T a tree obtained by removing one edge (between vertices u, w) from G, consider the vector x given by x(u) = c, and x(vi+1) = x(vi) -d, where the ordering is given by the tree, such that v1 = u, and vn = w, and the constants c, d are chosen such that x = 0 and x = 1. It is easy to see that LGx = nLT x. Moreover, since LGx = LT x for any x s.t. x(u) = x(w) = 0, it follows that κ(LG, LT ) = n.

We now prove Lemma 4, which also shows that κ(LG, HT ) ≤ 2 in the setting described by Lemma 3. For this, first note that λmax(LG, HT ) ≤ 2 for any G and T . This is because 2HT -LG = LG + AT c , where AT c is the adjacency matrix for the edge-complement of the tree T , can easily seen to be positive definite. It remains to show that λmin(LG, HT ) ≥ 1 if T c is a star graph. For this, it suffices to show that for any vector x, s.t. x T HT xconst, where xconst is the constant vector, we have x T LG -Hx ≥ 0. Note that y = HT xconst is a vector, such that y(u) = dT c (u) the degree of vertex u in the edge-complement of the tree T in G. Therefore, x T Hxconst = 0 if and only if u x(u)dT c (u) = 0. Now remark that

LG -H = -AT c and thus x T (LG -H)x = -(u,v)∈T c x(u)x(v). Since by assumption T c is a star, and x T HT xconst = 0 there exists a vertex w, s.t. x(w) =v, s.t. (v,w)∈T c x(v), and therefore, x T (LG -H)x =u s.t. (u,w)∈T c x(u)x(w) = v, s.t. (v,w)∈T c x(v) 2 ≥ 0. It follows that LT -HT is positive definite, λmin(LG, HT ) ≥ 1 and therefore κ(LG, HT ) ≤ 2.

Proof of Lemma 5: First note that if x is constant on the outer (cycle) part of the graph, then LGx = LT x = nx. Moreover, if x is 0 on the center of the star tree and alternating on the cycle: x(ui+1) = -x(ui) then LGx = 5x, whereas LT x = x. Therefore κ(LG, LT ) ≥ 5 (In fact, by using the same argument as in the proof of 2, one can see that κ(LG, LT ) = 5). At the same time, it is easy to see that λmax(LG, HT ) = 5/3 and λmin(LG, HT ) > 1/3 for all n. Therefore κ(LG, HT ) < κ(LG, LT ).


Abdelhamid Makiou
email: makiou@telecom-paristech.fryoucefbegriche@ieee.orgahmed@telecom-paristech.fr

Youcef Begriche

Ahmed Serhrouchni

Hybrid Approach to Detect SQLi Attacks and Evasion Techniques

Keywords:

Injections flaws which include SQL injection are the most prevalent security threats affecting Web applications[1].

To mitigate these attacks, Web Application Firewalls (WAFs) apply security rules in order to both inspect HTTP data streams and detect malicious HTTP transactions. Nevertheless, attackers can bypass WAF's rules by using sophisticated SQL injection techniques. In this paper, we introduce a novel approach to dissect the HTTP traffic and inspect complex SQL injection attacks. Our model is a hybrid Injection Prevention System (HIPS) which uses both a machine learning classifier and a pattern matching inspection engine based on reduced sets of security rules.

I. INTRODUCTION

Structured Query Language (SQL) injection is one of the most devastating vulnerabilities that impacts DataBase Management Systems (DBMS), as it can lead to the exposure of all the sensitive information stored in an application's database [START_REF] Ansari | SQL Injection in Oracle: An exploration of vulnerabilities International Journal on Computer Science and Engineering[END_REF]. In order to confront SQL injection attacks, various methodologies and techniques have been used. On one hand, Web application developers adopted safe coding and applied input validation functions. They developed filters to protect the application's entries from SQL code injections. These filters block inputs that contain SQL keywords or special characters commonly used in malicious SQL code injection. On the other hand, Web Application Firewalls protect application's database by inspecting HTTP traffic and applying a set of security rules. The language used for expressing security rules can explicitly describe a signature of an SQL injection attack, or implicitly describe the way of detecting these attacks. It can also express an anomaly score value which increases every time a malicious pattern appears in an HTTP request. If the anomaly value reaches a predefined threshold, the request will be rejected. In spite of the robustness of the above methods, attackers can bypass them by substituting malicious pattern characters or varying its format. However, a basic solution is to write a specific rule for each type of evasion technique, but it requires a high mastery of both HTTP protocol and regular expressions programming. Furthermore, pattern matching algorithms, used by security rules in order to inspect complex patterns, decrease the overall performances of the detection engine. Since pattern matching algorithms require a lot of resources, some WAFs [START_REF] Ristic | ModSecurity Handbook: The Complete Guide to the Popular Open Source Web Application Firewall[END_REF] [9] are configured to only inspect POST request. In this paper, we propose a hybrid approach to detect SQL injection attacks and their evasion techniques. Our proposal enhances both the inspection process of HTTP streams and security rules management.

II. RELATED WORK

In [START_REF] Tajpour | Comparison of SQL Injection Detection and Prevention Techniques, 2nd International Conforence on Education Technology and Computer (ICETC)[END_REF][4], Kruegel and Vigna propose an anomaly-based intrusion detection system for web applications. It characterizes HTTP requests using a number of statistical characteristics derived from parameter's length, character distribution, structure, presence and order. This method focuses only on the incoming query parameters whereas it ignores the corresponding HTTP response. These results are either causing unnecessary false positives or missing certain attacks. AMNESIA [START_REF] Halfond | AMNESIA: Analysis and Monitoring for NEutralizing SQL-Injection Attacks[END_REF] is an SQL injection detection and prevention system which combines static analysis and run-time monitoring. It uses a model-based approach to detect illegal queries. Nonetheless, it requires web application's source code reviewing. In SQLrand [3][6] instead of normal SQL keywords developers create queries using randomized instructions. In this approach a proxy filter intercepts queries to the database and de-randomizes the keywords. By using the randomized instruction set, attacker's injected code could not have been constructed. As it uses a secret key to modify instructions, security of the approach is dependent on attacker ability to seize the key. It requires the integration of a proxy for the database in the system as the same as developer training. ModSecurity WAF, proposed by Ivan Ristic [START_REF] Ristic | ModSecurity Handbook: The Complete Guide to the Popular Open Source Web Application Firewall[END_REF], is an open source solution based on signature attack detection. ModSecurity is widely used and has medium performances. Though, this system is strongly related to some types of web servers and it only analyses POST queries in order to avoid performance deterioration. In addition, the rules formalism is very complex which requires a high expertise in HTTP protocol and in regular expressions coding. IronBee [START_REF]The IronBee Project[END_REF], a new project similar to ModSecurity, aims to improve detection performance and facilitate the expression of security rules by introducing the LUA scripting language. An other recent open source project NAXSI [START_REF]Naxsi project (Nginx Anti Xss Sql Injection[END_REF] uses a heuristic approach for the detection of XSS and SQL injection attacks. Its performances are acceptable but requires a learning process to define white-lists. SQLi and XSS rules are static and use a simple cumulative scoring system based on the appearance of some special strings.

III. SQLI ATTACKS AND EVASION TECHNIQUES

Halfond, Viegas, and Orso researches [START_REF] Halfond | AMNESIA: Analysis and Monitoring for NEutralizing SQL-Injection Attacks[END_REF] proposed a taxonomy of SQL injection attacks. Depending on the goal, attackers can append a syntactically correct SQL code to the original query, or forge their own malicious SQL commands and introduce them to the DBMS via a vulnerable web application inputs. They will use one of these classes: Tautology, Incorrect Queries, UNION Queries, Piggy-backed Queries, Stored Procedure, Inference (Blind Injection, Timing Attacks) . In the following, we will present two techniques attackers can use in order to bypass the WAF rules.

IV. PROPOSED SOLUTION

A. HTTP Protocol Dissection

In most security solutions, traffic dissection process is the first operation before applying any security control. In HAKA project [START_REF] Denis | Arkoon Network Security), un langage orient rseaux et scurit[END_REF], HTTP stream is divided in two tables, HTTPrequest and HTTP-response, each table is associated with one hook and each hook contains all security rules declared for either requests or responses. In this section we introduce an other manner to dissect HTTP streams.

a) Typical HTTP request: HTTP protocol is expressed in a human-readable ASCII text. Headers use text to describe a request from a client (browser) or a response from the server. An HTTP request begins usually with a GET or POST method, followed by the URL and the protocol version. The following headers provide various information about the client, connection, content, etc. These headers are separated by \r\n to distinguish each header. which are separated by \r\n characters. However, before making the dissection, it has to get information about security rules. Indeed, users are obliged to declare security rules for the body and for each header. With the knowledge of headers involved in the inspection process, the dissector will only extract and parse these headers.

V. THE HYBRID INJECTION PREVENTION SYSTEM (HIPS)

• 1 Request Dissector URL: The information that the URL must be inspected is known by the dissector, since the user has declared a security rule on the hook : HTTP URL. The dissector extracts URL string and passes it to the classifier. A. The Machine Learning Model 1) SQL injections data collection framework: Malicious traffic is collected from an attack platform that includes: -Specific SQLi attack tools using evasion techniques -A Web server -A vulnerable web application with known SQL injection attacks.

B. features vector selection

Feature selection plays an important role in the identification of the potential malicious data used to escape SQLi filters. taken by the random variables X 1 , . . . , X m that are assumed conditionally independent relative to the category c (SQLi, Legitimate). Each random variable gives an information about a pattern type in a dissected header. In this model, all random variables are binary: X i = 1 if pattern of type i noted pa i is present, otherwise X i = 0. Consequently, each random variable X i = 1 follows a Bernoulli distribution with parameter p i = p(pa i ).

D. Classification: Naive Bayesian Classifier model

According to the Bayes theorem [START_REF] Robert | Le choix Baysien. Principes et pratiques[END_REF] and the total probabilities theorem, for a vector -→ x = (x 1 , . . . , x n ), the probability to belong to the class c is defined as follows:

p(C = c/ - → X = - → x )= p(C = c)•p( - → X = - → x /C=c) p( - → X = - → x ) (1)
Using the theorem of the total probabilities, we deduce:

p(C = c/ - → X = - → x )= p(C = c)•p( - → X = - → x /C=c) c {SQLi,Leg} p(C=c)•p( - → X = - → x /C=c) (2)

E. Cost Evaluation

Works of [START_REF] Androutsopoulos | An Evaluation of Naive Bayesian Anti-Spam Filtering[END_REF] on cost-sensitive evaluation measures consist of evaluating the false positives effect on the total cost, in term of time wasting by users to delete spams. But in our case, the cost is the impact of false negatives on the trust of users granted to our classifier.

1) The False Negatives Effect: A false negative is mistakenly classifying an SQLi attack as a legitimate content, and a false positive is a legitimate content mistakenly classified as an SQLi attack. In our model, the cost of a false negative is much higher than the cost of a false positive. Indeed, wasting time in analyzing legitimate requests is more acceptable than passing a malicious code to the Web application. The two error types are defined as such follows:

• Classifying an SQLi attack as legitimate content:

(SQLi → Leg) • Classifying a legitimate content as an SQLi attacks:

(Leg → SQLi) In our classifier model, the first error is more serious than the second error. To illustrate this idea, we introduce the parameter λ, as its objective is to give more importance to the first error by assuming that Leg → SQLi is λ times more costly than SQLi → Leg.

2) Classification criteria: According to the above two error types, the selection criteria is as follows:

The content of a header x is classified as legitimate if and only if:

p(C = Leg/ X = x) > λ.p(C = SQLi/ X = x) (3)
given that p(C = Leg/ X = x) + p(C = SQLi/ X = x) = 1, the selection criteria becomes as follows:

p(C = Leg/ X = x) > α (4)
Where :

α = λ 1 + λ , λ = α 1 -α (5)
3) Method and parameters evaluation: In this section, we define the parameters that allow us to evaluate our filter. To this end, two evaluation parameters are used: accuracy (Acc) and error (Err = 1 -Acc) [START_REF] Androutsopoulos | An Evaluation of Naive Bayesian Anti-Spam Filtering[END_REF]. They are defined as follows:

Acc = n sqli→sqli + n leg→leg N sqli + N leg (6)
Err = n sqli→leg + n sqli→sqli N sqli + N leg
where:

• N sqli = n sqli→leg + n sqli→sqli • N leg = n leg→leg + n leg→sqli
• n y→z denotes the number of patterns of class y that are mistakenly classified in class z. The parameters defined above do not take into consideration the notion of weight for the two error types introduced in the previous paragraph. This leads us to introduce the weighted accuracy (W acc) and weighted error (W err = 1-W acc). We assumed that SQLi → Leg is λ times more devastating for our system than Leg → SQLi. To make accuracy and error rate sensitive to this cost, we should treat each SQL injection as if it was λ inputs; when an SQL injection is misclassified, this counts as λ errors; and when it is classified correctly, this counts as λ successes

W acc = λn sqli→sqli + n leg→leg λN sqli + N leg (7)
W err = λn sqli→leg + n leg→sqli λN sqli + N leg
To have a precise idea of the filter's performance, we compare it to a non-filtered system in which all requests are considered as legitimate. We introduce the definition of the base-line weighted error and the base-line weighted accuracy (respectively noted W acc b and W err b ) which are defined as follows:

W acc b = λN leg λN sqli + N leg (8)
W err b = N sqli λN sqli + N leg

The T CR (Total Cost Ratio) value measures the performance of a machine learning classifier to the same environment without a classifier. In the case where the T CR value is negligible, the best approach is to not use a classifier and to send all requests to the rules based detection engine. An effective filter which could be used in real environments should have a T CR value higher than 1.

The T CR formula is defined as follows:

T CR = W err b W err = N sqli λn leg→sqli + n sqli→leg (9)

VI. RESULTS EVALUATION

We have collected a training data set from the framework previously presented in this paper. The training set is a mixture of SQli attacks and legitimate requests shown in the below table with different proportions (80-20)%,(75-25)%,(50-50)%, then we varied, for each scenario, the λ parameter and finally calculated values of W acc, W acc b and TCR.

By increasing the α (threshold) value from 99% to 50%, we have an increase in number of false positives, which means that the rules engine will analyze legitimate requests, but at the same time, the evaluation has shown an increase in TCR value. However, in practice, false positives (legitimate requests classified as SQLi attacks) will be forwarded to the rules engine that will apply all security rules. This will decrease the overall performances of our system. To find a compromise between lower number of false positives and false negatives, the TCR values should be medium values. Acceptable values of TCR are related to α threshold values higher than 50% and lower than 90%.

VII. COMPARATIVE ANALYSIS

Security Rules based WAFs do not have a module that predicts SQL injection attacks. In this case, our approach will improve the inspection performance because the classifier module will not forward legitimate traffic to the detection engine. Only suspected requests are deeply inspected by applying pattern matching algorithms. In [START_REF] Komiya | Classification of malicious web code by machine learning[END_REF], Authors worked on Classification of Malicious Web Code by Machine Learning. They implemented and evaluated two classifiers both SQLIAs and XSS which can use TF-IDF method for weight calculation, and three machine learning approach among SVM, Naive-Bayes, k-Nearest Neighbor Algorithm. They obtain good precision values for their classifier (99.16 %) by using SVM with Gaussian Kernel. Our classifiers obtains 97.6% by using Bayesian algorithm. However, they did not provide a solution to handle false negatives and false positives. On the other hand, the using of SVM algorithms with Gaussian Kernel may require significant CPU resources and decreases the classifier performances in multi-gigabits rates networks.

VIII. CONCLUSION

In this paper, we focused on the problem of detecting complex SQL injections. We proposed a novel approach to dissect HTTP requests in order to cover most evasion techniques and improve security rules management process. We also provided an Injection Prevention System architecture which includes a machine learning classifier. Based on the TCR results, we have shown the effectiveness of the classifier by tuning its values in order to reduce false negatives. We were also able to show that false positives didn't impact the overall performances of our system. A key element of future work is to apply the same approach in order to develop an anti XSS and SQL attacks solution.

Fig. 1 :

 1 Fig. 1: Raw HTTP Request b) HTTP Request Dissection: Our dissection module is able to recognize request's components (headers and the body)which are separated by \r\n characters. However, before making the dissection, it has to get information about security rules. Indeed, users are obliged to declare security rules for the body and for each header. With the knowledge of headers involved in the inspection process, the dissector will only extract and parse these headers.

Fig. 2 :

 2 Fig. 2: Detection Engine Architecture

Fig. 3 :

 3 Fig. 3: Attack Traces Collection Framework

Fig. 4 :

 4 Fig. 4: Attack Traces Collection Framework

  table shows the result of our features vector selection.

	Variables	Tokens (key words)
	x1	SELECT
	x2	UNION
	x3	UPDATE
	..	Other SQL keywords
	xn	* (Asterisk)
	..	Other Special Symbols
	xm	UNHEX
	..	Other Evasion Keywords
	x45	% (percent)

TABLE I :

 I

A simple table C. Request representation Each request's header is characterized by a vector -→ x defined by -→ x = (x 1 , . . . , x m ) where x 1 , . . . , x m are the values

TABLE II :Table 2

 II2