Umit V Çatalyürek

Kamer Kaya

Bora Uçar
email: bora.ucar@ens-lyon.fr

On analysis of partitioning models and metrics in parallel sparse matrix-vector multiplication

Keywords:

Graph/hypergraph partitioning models and methods have been successfully used to minimize the communication requirements among processors in several parallel computing applications. Parallel sparse matrix-vector multiplication (SpMxV) is one of the representative applications that renders these models and methods indispensable in many scientific computing contexts. We investigate the interplay of several partitioning metrics and execution times of SpMxV implementations in three libraries: Trilinos, PETSc, and an in-house one. We design and carry out experiments with up to 512 processors and investigate the results with regression analysis. Our experiments show that the partitioning metrics, although not an exact measure of communication cost, influence the performance greatly in a distributed memory setting. The regression analyses demonstrate which metric is the most influential for the execution time of the three libraries used.

Sur l'analyse des modèles de partitionnement et des métriques lors de la multiplication d'une matrice creuse avec un vecteur dense en parallèle Résumé : Les modèles et méthodes de partitionnement de graphes/hypergraphes ont été utilisés avec succès pour minimiser les besoins de communication entre processeurs dans de nombreuses applications de calcul parallèle. La multiplication d'une matrice creuse avec un vecteur dense (SpMxV) est une des applications représentatives qui rendent ces modèles et méthodes indispensables dans de nombreux contextes de calcul scientifique. Nous nous intéressons aux interactions entre plusieurs métriques de partitionnement et le temps d'exécution des implémentations de SpMxV dans les bibliothèques Trilinos, PETSc, et une bibliothèque à nous. Nous effectuons des expériences avec jusqu'à 512 processeurs, et nous étudions les résultats à l'aide d'analyse par régression. Nos expériences montrent que les métriques de partitionnement, bien qu'elles ne fournissent pas une mesure exacte du coût de communication, influent de façon significative sur la performance dans un système à mémoire distribuée. Les analyses par régression montrent quelle métrique a le plus d'influence sur le temps d'exécution pour les trois bibliothèques utilisées.

Mots-clés : Multiplication d'une matrice creuse avec un vecteur en parallèle, partitionnement d'hypergraphes 1 Introduction

Repeated sparse matrix-vector (SpMxV) and sparse matrix-transpose-vector multiplies that involve the same large, sparse matrix are the kernel operations in various iterative algorithms involving sparse linear systems. Such iterative algorithms include solvers for linear systems, eigenvalues, and linear programs. Efficient parallelization of SpMxV operations is therefore very important in virtually all large scale scientific computing applications. Quite a number of partitioning methods and models have been used to enable efficient parallelization of SpMxV. Some of the earlier methods try to balance the computational load of processors without paying attention to the effect of sparsity in the communication requirements. Some more recent methods use graph and hypergraph models to minimize some aspect of the communication cost while achieving load balance.

Hypergraph model-based partitioning methods have been proved to be very useful in implementing different partitioning schemes. These partitioning schemes address different communication cost metrics for the suitable variants of parallel SpMxV operations. In general, the importance of the communication cost metrics, such as the total volume of communication, the total number of messages and these two quantities on per processor basis, depends on the machine architecture, problem size, and the underlying parallel algorithm. In this study, we investigate the effects of the partitioning methods in order to identify the most relevant metrics in various configurations.

Hypergraph partitioning models try to reduce the total communication volume. Optimizing this metric seems to be more tractable than optimizing the others. However, in general, optimizing only this metric cannot always reduce the parallel execution time. This has been recognized by earlier models that implicitly and/or explicitly address other communication metrics. Here we mention two and defer the discussion of related work to Section 5.

Uçar and Aykanat [START_REF] Uçar | Encapsulating multiple communication-cost metrics in partitioning sparse rectangular matrices for parallel matrix-vector multiplies[END_REF] reduce the total number of messages by allowing controlled increase in the total communication volume with respect to the standard hypergraph models [START_REF] Çatalyürek | Hypergraph-partitioning-based decomposition for parallel sparse-matrix vector multiplication[END_REF], while also trying to achieve balance on other communication cost metrics on per processor basis. They report practical experiments with parallel SpMxV computations showing that the method addressing the multiple communication cost metrics is more helpful than the original method that addresses only the total communication volume. Bisseling and Meesen [START_REF] Bisseling | Communication balancing in parallel sparse matrix-vector multiplication[END_REF] also propose methods addressing different communication cost metrics and report improvements with respect to an original method. Both of these studies demonstrate that different communication cost metrics and their interplay can be important to achieve scalable parallel algorithms. It is therefore important to understand the effects of different metrics (optimized by different partitioning models) on the running time of applications under different configurations.

The contribution of this paper is two-fold. We designed and conducted several experiments in a system with 512 processors to show the effects of partitioning models and metrics on SpMxV performance. As far as we know, this is the first work which compares the existing partitioning models and metrics in modern architectures with modern software following message-passing paradigm. Second, we showed that it is difficult, if not impossible, to define the correct partitioning model and metric without analyzing the characteristics of the input matrices and the SpMxV library being used. We experimented

RR n°8301

with three existing libraries, PETSc [START_REF] Balay | Efficient management of parallelism in object oriented numerical software libraries[END_REF][START_REF] Balay | PETSc users manual[END_REF], Trilinos [START_REF] Heroux | An overview of the trilinos project[END_REF], and an in-house library SpMV [START_REF] Uçar | A library for parallel sparse matrix-vector multiplies[END_REF]. In order to overcome the mentioned difficulty, we carefully analyze the results using regression analysis techniques and relate the execution time of SpMxV implementations to the different partitioning metrics. We portray this analysis in detail so as to suggest improved objective functions for partitioning software. Our experiments show that partitioning with the correct metric and model is important for good performance especially when the number of processors is large. Although, we only had an access to a 512-processor machine, the experiments and their analysis show that to scale larger systems, one needs to be more careful while partitioning the matrix-in our experiments the fact that the communication metrics greatly related to the execution time is observable starting from 64 processors.

The rest of the paper is organized as follows. In Section 2, we describe SpMxV in detail with three parallel algorithms using different approaches. Section 3 defines the partitioning models and metrics investigated in this work. Experimental results are given in Section 4 and some related work is mentioned in Section 5. Section 6 concludes the paper.

Parallel sparse matrix vector multiplication

Consider the sparse matrix-vector multiply of the form y ← Ax, where the nonzeros of the m × n matrix A are partitioned arbitrarily among K processors such that each processor P k owns a mutually disjoint subset of nonzeros, A (k) where A = k A (k) . The vectors y and x are also partitioned among processors, where the processor P k holds x (k) , a dense vector of size n k , and it is responsible for computing y (k) , a dense vector of size m k . Let µ(•) denote the nonzero-toprocessor and vector-entry-to-processor assignments. We note that the vectors x (k) for k = 1, . . . , K are disjoint and hence k n k = n; similarly the vectors y (k) for k = 1, . . . , K are disjoint and hence k m k = m. In this setting, the sparse matrix A (k) owned by processor P k can be permuted and written as

A (k) =          A (k) 11 • • • A (k) 1 • • • A (k) 1K . . . . . . . . . . . . . . . A (k) 1 • • • A (k) • • • A (k) K . . . . . . . . . . . . . . . A (k) K1 • • • A (k) K • • • A (k) KK          . (1)
Here, the blocks in the row-block stripe

A (k) k * = {A (k) k1 , . . . , A (k) kk , . . . , A (k)
kK } have row dimension of m k , and similarly the blocks in the column-block stripe

A (k) * k = {A (k) 1k , . . . , A (k) kk , . . . , A (k)
Kk } have column dimension of n k . The x-vector entries that are needed by processor P k are represented as

x(k) = [x (k) 1 , . . . , x(k) k , . . . , x(k) K ]
, a sparse column vector, where x(k) contains only those entries of x ( ) of processor P corresponding to the nonzero columns in A (k) * . Here, the vector k) , defined according to the given partition on the x-vector (hence the vector x(k) is of size at least n k ). The y-vector entries for which the processor P k computes partial results are represented as a sparse vector ŷ

x(k) k is equivalent to x (
(k) = [ŷ (1) k , . . . , ŷ(k) k , . . . , ŷ(K) k ],
where ŷ( ) k contains only the

RR n°8301

partial results for y ( ) corresponding to the nonzero rows in A (k) * . If a row i in A (k) contains nonzeros only in columns corresponding to x (k) , that is if µ(a ij ) = P k for all nonzero a ij ∈ A (k) , and µ(x i ) = P k , the row i is called a local row from the perspective of the processor P k . Similarly, if a column j in A (k) contains nonzeros only in rows corresponding to ŷ(k) k , that is if µ(a ij ) = P k for all nonzero a ij ∈ A (k) , and µ(y i ) = P k , the column j is called a local column, from the perspective of the processor P k .

The standard parallel SpMxV algorithm [START_REF] Çatalyürek | On two-dimensional sparse matrix partitioning: Models, methods, and a recipe[END_REF][START_REF] Uçar | A library for parallel sparse matrix-vector multiplies[END_REF][START_REF] Vastenhouw | A two-dimensional data distribution method for parallel sparse matrix-vector multiplication[END_REF] based on the described nonzero and vector entry partitioning is called the row-column-parallel algorithm. In this algorithm, each processor P k executes the following steps.

(k) = [x (k) 1 , . . . , x(k) k , . . . , x(k) K ], first define x(k) k = x (k) . Then, for each = k where A (k) * contains nonzeros, receive x(k) from processor P = µ(x (k) ), corresponding to the nonzero columns in A (k) * . 3. Compute ŷ(k) ← A (k) x(k) . 4.
For each = k, send the sparse partial-results vector ŷ( ) k to processor P = µ(ŷ ( ) k ), where ŷ( ) k contains only those partial results for y ( ) corresponding to the nonzero rows in A (k) * . 5. Receive the partial-results vector ŷ(k) from each processor P which has computed a partial result for y (k) , i.e., from each processor P where A ( ) k * has nonzeros. 6. Compute y (k) ← ŷ(k) , adding all the partial-results ŷ(k) received in the previous step to its own partial results for y (k) .

There are two communication phases in this algorithm. The first one is just before the local matrix-vector multiply, and it is due to the communication of the x-vector entries (steps 1 and 2). This communication operations is referred to as expand. If the matrix A is distributed columnwise, and the x-vector entries are partitioned conformably with the column partition of A, then the expand operation is not needed. The resulting algorithm is called the column-parallel algorithm. The second communication phase is just after the local matrix-vector multiply (steps 4 and 5), and it is due to the communication of the partial results on the y-vector entries. This communication operation on the y-vector entries is referred to as fold. If the matrix A is distributed rowwise, and the y-vector entries are partitioned conformably with the rows of A, then the fold operation is not needed. The resulting algorithm is called the row-parallel algorithm.

There are different implementations of the above algorithm. We summarize three implementations with which we have experimented. Two of the implementations are in the well-known general libraries Trilinos [START_REF] Heroux | An overview of the trilinos project[END_REF] and PETSc [START_REF] Balay | Efficient management of parallelism in object oriented numerical software libraries[END_REF][START_REF] Balay | PETSc users manual[END_REF]; the third one, SpMV, is an in-house library [START_REF] Uçar | A library for parallel sparse matrix-vector multiplies[END_REF].

Trilinos provides an implementation which can be described as in Algorithm 1 from the point of view of the processor P k . In this implementation, the expand operations are finished before doing any computation. Then, all the scalar multiply-add operations are performed. Later on, the fold operations are completed. Trilinos uses Irecv/Isend and waitall communication primitives PETSc provides an implementation of the above algorithm only for the rowparallel case. The Algorithm 2 summarizes that implementation from the point of view of P k . First, the expand operation is initiated using Irecv and Isend primitives. Then, instead of waiting the reception of all necessary x-vector entries, it performs some local computations so as to overlap communication and computations. In particular, the processor P k performs scalar multiply-add operations using local a ij 's for which µ(x j ) = P k and there is no a i with µ(x ) = P k . Then, upon verifying the reception of all needed x-vector entries using waitall, P k continues with scalar multiply-add operations with the nonzeros on the rows that has at least one nonzero in a column j for which µ(x ) = P k . The implementation can also be seen in an earlier technical report [START_REF] Saad | P-SPARSLIB: A portable library of distributed memory sparse iterative solvers[END_REF]. Figure 1 describes the algorithm pictorially. After issuing Isends and Irecvs (for x(k) ), processor P k performs the computations associated with the horizontally shaded matrix zone. Then, waitall is executed to have all x (k) before continuing with the rows that are below the horizontal ones. Note that the local matrices are actually permuted into the displayed form (local rows and the interface rows). The advantage of this implementation with respect to the Algorithm 1 is that it allows overlap between the reception of messages for the expand operation and scalar multiply-add operations with the nonzeros in local rows.

Consider again the matrix A (k) of processor P k as shown in Fig. 1. Before executing the waitall operation, there are some more scalar multiply-add operations that P k can perform before the reception of any

x(k) . These operations

x (k) ∧ 1 ... x (k) ∧ k ... x (k) ∧ K y (k) ∧ Figure 1:
The zones of the matrix A (k) of processor P k with respect to the vector

x assuming a row-parallel algorithm. PETSc exploits communication-computation overlap using only the horizontal zone. It is possible to use the hatched zone as well. Furthermore, upon the reception of a message we have association computations to do with the vertical zone.

are related to the nonzeros that are in the hatched zone in the figure. In order to exploit the hatched zone for communication computation overlap, one can store that zone in the compressed column storage (CCS) format. This way, one can delay the invocation of the waitall operation for some more time. In fact, we can get rid of the waitall operation and maximize the communication computation overlap by performing all scalar multiply-operations that involve a received x-vector entry before waiting the reception of any other message. This requires storing the vertically shaded zones of the matrix in Fig. 1 in CCS format: with this, when P k receives x(k) , it can visit the respective column and perform all operations. This way of storing the vertically shaded and hatched zones in CCS maximizes the amount of overlap in the strict sense (optimal amount of overlap) when a processor receives a single message from each sender (as should be the case in a proper SpMxV code). The third library that we investigate in this work SpMV [START_REF] Uçar | A library for parallel sparse matrix-vector multiplies[END_REF] implements this approach for row-parallel and row-column parallel algorithms. The row-column-parallel algorithm's implementation of SpMV is shown in Algorithm 3 for processor P k . A K-way partition of a hypergraph H is denoted as Π = {V 1 , V 2 , . . . , V K } where each part is nonempty subset of V, they are pairwise disjoint, and their union is equal to V. In a partition Π, a net that has at least one pin (vertex) in a part is said to connect that part. The number of parts connected by a net n j , i.e., connectivity, is denoted as λ j . A net n j is said to be uncut (internal) if it connects exactly one part (i.e., λ j = 1), and cut (external), otherwise (i.e., λ j > 1).

Let W k denote the total vertex weight in V k (i.e.,

W k = v∈V k w[v]
) and W avg denote the weight of each part when the total vertex weight is equally distributed (i.e.,

W avg = ( vi∈V w[v i ])/K). If each part V k ∈ Π satisfies the balance criterion W k ≤ W avg (1 + ε), for k = 1, 2, . . . , K (2)
we say that Π is balanced where ε represents the maximum allowed imbalance ratio. The set of external nets of a partition Π is denoted as N E . Let χ(Π) denote the cost, i.e., cutsize, of a partition Π. There are various cutsize definitions [START_REF] Lengauer | Combinatorial Algorithms for Integrated Circuit Layout[END_REF]. The following two are the most relevant ones for this work:

χ(Π) = nj ∈N E c[n j ] (3)
χ(Π) = nj ∈N c[n j ](λ j -1) . (4)
In ( 3) and ( 4), the contributions of each cut net n j to the cutsize are c[n j ] and c[n j ](λ j -1), respectively. The cutsize metric given in (3) will be referred to here as cut-net metric and the one in (4) will be referred as connectivity metric. Given ε and an integer K > 1, the hypergraph partitioning problem can be defined as the task of finding a balanced partition Π with K parts such that χ(Π) is minimized. The hypergraph partitioning problem is NP-hard [START_REF] Lengauer | Combinatorial Algorithms for Integrated Circuit Layout[END_REF].

Hypergraph models and methods for sparse matrix partitioning

There are three basic hypergraph models for sparse matrix partitioning: the column-net hypergraph model, the row-net hypergraph model, and the fine-grain model. Each of these models also directly corresponds to a partitioning method.

RR n°8301

We will briefly describe these models below. It is also possible to develop other partitioning methods, like jagged-like, checker-board, Mondriaan (a variant of orthogonal recursive bisection) using these three basic models [START_REF] Bisseling | Two-dimensional approaches to sparse matrix partitioning[END_REF][START_REF] Çatalyürek | On two-dimensional sparse matrix partitioning: Models, methods, and a recipe[END_REF][START_REF] Uçar | A matrix partitioning interface to PaToH in MATLAB[END_REF][START_REF] Vastenhouw | A two-dimensional data distribution method for parallel sparse matrix-vector multiplication[END_REF] with the standard hypergraph partitioning tools, such as PaToH [START_REF] Çatalyürek | PaToH: A Multilevel Hypergraph Partitioning Tool[END_REF] or Zoltan [START_REF] Boman | Zoltan 3.0: Parallel Partitioning, Load Balancing, and Data-Management Services; User's Guide[END_REF]. In this work, among other partitioning methods we will only use the checker-board method, because it is usually the best method to reduce the maximum number of messages sent and received by a processor.

In the column-net hypergraph model [START_REF] Çatalyürek | Hypergraph-partitioning-based decomposition for parallel sparse-matrix vector multiplication[END_REF] used for 1D rowwise partitioning, an m × n matrix A with τ nonzeros is represented as a unit-cost hypergraph

H R = (V R , N C ) with |V R | = m vertices, |N C | = n
nets, and τ pins. In H R , there exists one vertex v i ∈ V R for each row i of matrix A. Weight w[v i ] of a vertex v i is equal to the number of nonzeros in row i. The name of the model comes from the fact that the columns are represented as nets. That is, there exists one unit-cost net n j ∈ N C for each column j of matrix A. The net n j connects the vertices corresponding to the rows that have a nonzero in column j. That is, v i ∈ n j if and only if a ij = 0.

In the row-net hypergraph model [START_REF] Çatalyürek | Hypergraph-partitioning-based decomposition for parallel sparse-matrix vector multiplication[END_REF] used for 1D columnwise partitioning, an m × n matrix A with τ nonzeros is represented as a unit-cost hypergraph

H C = (V C , N R ) with |V C | = n vertices, |N R | = m
nets, and τ pins. In H C , there exists one vertex v j ∈ V C for each column j of matrix A. Weight w[v j ] of a vertex v j ∈ V R is equal to the number of nonzeros in column j. The name of the model comes from the fact that the rows are represented as nets. That is, there exists one unit-cost net n i ∈ N R for each row i of matrix A. The net n i ⊆ V C connects the vertices corresponding to the columns that have a nonzero in row i. That is, v j ∈ n i if and only if a ij = 0.

In the fine-grain model [START_REF] Çatalyürek | A fine-grain hypergraph model for 2D decomposition of sparse matrices[END_REF], otherwise known as the column-row-net hypergraph model, used for 2D nonzero-based fine-grain partitioning, an m×n matrix A with τ nonzeros is represented as a unit-weight and unit-cost hypergraph

H Z = (V Z , N RC ) with |V Z | = τ vertices, |N RC | = m +
n nets and 2τ pins. In V Z , there exists one unit-weight vertex v ij for each nonzero a ij of matrix A. The name of the model comes from the fact that both rows and columns are represented as nets. That is, in N RC , there exists one unit-cost row-net r i for each row i of matrix A and one unit-cost column-net c j for each column j of matrix A. The row-net r i connects the vertices corresponding to the nonzeros in row i of matrix A, and the column-net c j connects the vertices corresponding to the nonzeros in column j of matrix A. That is, v ij ∈ r i and v ij ∈ c j if and only if a ij = 0.

Three standard partitioning methods are obtained by partitioning the vertices of the hypergraphs generated by one of these models using the connectivity metric (4): Column-Net (CN), Row-Net (RN) and Fine-Grain (FG). Each one of those methods minimizes the total communication volume, while keeping the balanced work.

The checker-board partitioning method (CB) [START_REF] Çatalyürek | A hypergraph-partitioning approach for coarse-grain decomposition[END_REF] is a two-step method, where each step models either the expand phase or the fold phase of the parallel SpMxV. There are two alternative schemes for this partitioning method: the one which models the expands in the first step and the folds in the second step, and vice-versa. We describe the first one below (the other can be described similarly). Given an m×n matrix A and the number K of processors organized as a M ×N mesh, first, A is partitioned row-wise into M parts using the columnnet model, producing Π R = {R 1 , . . . , R M }. In the second step, the matrix A is RR n°8301 partitioned columnwise into N parts by using the multi-constraint partitioning to obtain Π C = {C 1 , . . . , C N }. The rowwise and columnwise partitions Π R and Π C together define a 2D partition of the matrix A, where µ(a ij ) = P k, ⇔ r i ∈ R k and c j ∈ C . The multi-constraint formulation is used to achieve load balance among processors. In this formulation, each vertex v i of H C is assigned M weights w[i, k], for k = 1, . . . , M . Here, w[i, k] is equal to the number of nonzeros of column c i in rows R k . Consider an N -way partitioning of H C with M constraints using the vertex weight definition above. Maintaining the M balance constraints (each individual constraint treated as separate balance constraint) corresponds to maintaining computational load balance on the processors of each row of the processor mesh. Establishing the equivalence between the total communication volume and the sum of the cutsizes of the two partitions is fairly straightforward. The volume of communication for the fold operations corresponds exactly to the χ(Π C ). The volume of communication for the expand operations corresponds exactly to the χ(Π R ).

Partitioning metrics

We investigate the practical effects of partitioning models on several metrics. These metrics include the maximum number of nonzeros assigned to a processor (MaxNnz), the total communication volume (TotVol), the maximum send volume of a processor (MaxSV), the total number of messages (TotMsg), and the maximum number of messages sent by a processor (MaxSM). We will also analyze how these metrics are related with the SpMxV performance.

As can be expected, most of our results are based on experiments and their careful analysis (see the next section). Yet, there are some theoretical results worth mentioning. In 1D partitioning, the maximum number of messages can be as high as K -1, and the total number of messages can be K times larger. In the fine-grain model the maximum number of messages can be as high as 2 × (K -1), and the total number of messages as can be, again, K times larger. When a checker-board partitioning is used for distributing a matrix into K = M N processors, each processor only communicates with up to M -1 and N -1 processors in the expand and fold phases. Since the checker-board model is more restricted, the total communication volume is expected to be larger.

Experimental results

We carried our experiments on a 64-node cluster where each node has a 2.27GHz dual quad-core Intel Xeon (Bloomfield) CPU and 48GB main memory. Each core in a socket has 64KB L1 and 256KB L2 caches, and each socket has an 8MB L3 cache shared by 4 cores. The interconnection network is 20Gbps DDR InfiniBand. For parallelism, mvapich2 version 1.6 is used. We built SpMV, PETSc, and Trilinos with gcc 4.4.4 and used optimization flag -O3. For PETSc experiments, we used the matrix type MPIAIJ and the multiplication routine MatMult. Since each node has 8 cores, we have 512 processors in total. In the experiments, we use K ∈ {1, 8, 16, 32, 64, 128, 256, 512}. For an experiment with K = 1 processors, we fully utilize K/8 nodes of the cluster. To measure the time of one SpMxV operation (in secs), we do 500 multiplications for each execution. The values in the tables and figures are the averages of these 500 runs. We used seven large real-life square matrices from different application domains that are available at the University of Florida (UFL) Sparse Matrix Collection (http://www.cise.ufl.edu/research/sparse/matrices) and three synthetically generated matrices corresponding to 5-point stencil meshes in 2D with sizes 1024 × 1024, 2048 × 2048, and 4096 × 4096. The properties of the matrices are given in Table 1.

To generate the partitions with respect to column-net (CN), row-net (RN), fine grain (FG), and checkerboard (CB) models, we used PaToH [START_REF] Çatalyürek | PaToH: A Multilevel Hypergraph Partitioning Tool[END_REF] with the connectivity metric (4), i.e., aiming to minimize TotVol, and 0.03 as the maximum imbalance. We also used a block partitioning (BL) model to generate rowwise partitions of the matrices. In this model, we traverse the rows from 1 to m, generate a part with approximately τ /K nonzeros, and continue with the next part when this number is exceeded. Experimental results show that the average imbalance ratio of the partitions generated with BL is around 10 -4 which is two order magnitude smaller than what we allowed in other partitioning methods. The model does not take any other partitioning metric into account. For matrices based on 2D meshes, we used another rowwise partitioning model called MP. This model tiles the 2D plane with a diamond-like shape [4, Section 4.8] and associates each shape with a processor. This approach balances the volume and the number of messages the processors send and receive. In the experiments on real-life matrices, we use PETSc with rowwise models CN and BL, and SpMV and Trilinos with all models except MP. For meshes, we added MP to each library's model set.

Let us start with two small experiments on mesh-based matrices to demonstrate that the partitioning metrics are relevant. In the first experiment, we use MP which allows us to obtain the same MaxNnz, MaxSV, and MaxSM values for different number of processors (K) and mesh sizes. We exploit this to observe the effect of TotVol and TotMsg on the execution time as Table 2 shows. In this table, from one row to the other, only TotVol and TotMsg change. Yet, the execution time of each library increases significantly when TotVol increases from 8,200 to 131,296, and TotMsg increases from 40 to 752. Hence, a good partitioning with minimized metrics would reduce the execution time.

In the second experiment, we compared the metrics and execution times of Trilinos on mesh-based matrices when models BL and MP are used. Since MP takes TotVol and MaxNnz into account, and BL only deals with the latter, there is a drastic difference between mean TotVol and MaxSV values of these two models in favor of MP as Figure 2 shows. However, due to the banded structure of mesh-based matrices, BL obtains much less TotMsg and MaxSM values. The mean execution times of Trilinos employing these models are given in Fig. 2e. As the experiment shows, the total latency has a larger performance impact on Trilinos on these matrices in our system for larger processor counts, whereas the total volume of messages seems to be important for smaller number of processors. Clearly, one has to determine which metric is the most important for a given number of processors for a given algorithm. Below, we carry out a careful regression analysis of the metrics for different algorithms.

Regression analysis

To evaluate the performance of the libraries with respect to the partitioning metrics, we use linear regression analysis techniques and solve the nonnegative least squares problem (NNLS). In NNLS, given a variable matrix V and a vector t, we want to find a dependency vector d which minimizes Vdt s.t. d ≥ 0.

In our case, V has five columns which correspond to the partitioning metrics MaxNnz, TotVol, MaxSV, TotMsg, and MaxSM. Each row of V corresponds to an SpMxV execution where the execution time is put to the corresponding entry of t. Hence, we have the same V but a different t for each library. We apply a well-known technique in regression analysis and standardize each entry of V by subtracting its column's mean and dividing it to its column's standard deviation so that the mean and the standard deviation of each column become 0 and 1, respectively. This way, the units are removed and each column becomes equally important throughout the analysis. We then used MATLAB's lsqnonneg to solve NNLS. Each entry of the output d i shows the dependency of the execution time to the partitioning metric corresponding to the ith column of V. Tables 3,4, and 5 show the dependency values found in various settings.

We first apply regression analysis to each library with all matrices and rowwise partitioning models CN (column-net) and BL (block). The analysis shows that when K ≤ 64, SpMxV performance depends rigorously on the maximum number of nonzeros assigned to a processor. In this case, the dependency values for MaxNnz are 8.02, 7.81, and 6.80 for SpMV, PETSc, and Trilinos, respectively. As Table 3 shows, the next important metric is MaxSV with values 1.66, 1.53, and 2.20. The latency-based partitioning metrics (TotMsg and MaxSM) do not effect the performance for K ≤ 64. However, when K gets larger, these metrics are of utmost importance. Furthermore, the importance of MaxNnz decreases drastically for all the libraries. For SpMV and PETSc, MaxNnz becomes the 3rd important variable, whereas for Trilinos, it is the 4th. This shows that SpMV and PETSc handle the increase in the communication metrics better than Trilinos.

When K ≥ 128, the dependency of Trilinos to TotMsg is much less than that of SpMV and PETSc. On the contrary, Trilinos' MaxSM dependency is almost 1.75 times more than SpMV and PETSc. This is expected since Trilinos uses Algorithm 1 which has no communication-computation overlap due to the use of waitall primitive. Such primitives can cause close coupling among the processors. When MaxNnz and the variance on the number of messages per processor are large, the overhead due to the bottleneck processor can result in poor SpMxV performance. Note that the dependency profiles of SpMV and PETSc, which are similar due to the communication-computation overlap in Algorithms 2 and 3, do not point out a similar bottleneck. We extend the regression analysis to all matrices and all partitioning models and show the results in Table 4. The performance of SpMV and Trilinos rigorously depend on MaxNnz if K ≤ 32, and on TotMsg and MaxSM when K ≥ 256. Once again, Trilinos' MaxSM dependency is higher than that of SpMV due to the waitall primitive. To see the effect of matrix structure on regression analysis, we use only mesh-based matrices in the next experiment. As Table 5 shows, we observe that for these matrices, the performance of SpMV and Trilinos mostly depend on MaxNnz even when K ≥ 64. Note that these matrices are banded and the communication metrics have relatively lower values compared to those of real-life matrices. Hence, the most dominant factor is MaxNnz. In the light of the regression analysis experiments, we can argue that partitioning metrics in fact effect the performance of parallel SpMxV libraries. However, the best metric (or function of metrics) that needs to be minimized depends on the number of processors, the size and structure of the matrix, which we are planning to investigate as a future work, and even the library itself. Although some of these variables are known while generating the partitions, predicting the others may need a preprocessing phase. For example, we already know that the libraries in this paper are employing point-to-point communication primitives which makes the connectivity metric suitable. However, if collective communication primitives, e.g., MPI_ALLGATHER, had been used, it would be better to minimize the cut-net metric as the main partitioning objective (however, we should note that such collective operations introduce unnecessary synchronization and messages especially for large K values). On the other hand, the matrix structure can be different for each input and a partitioner needs either a manual direction or a preprocessing to predict the best metric for each matrix.

Effects of vector partitioning

For parallel y ← Ax with nonzero-to-processor and vector-entry-to-processor assignment µ, we assume the following: µ(x j ) (or µ(y i )) is equal to P , then there exists a nonzero a ij assigned to processor P . Although vector-entry-toprocessor assignments affect most of the communication metrics such as MaxSV and TotMsg, under this assumption, TotVol is not one of them. To minimize either of the partitioning metrics MaxSV and TotMsg while keeping TotVol constant, we use a simple vector-partitioning heuristic. The heuristic traverses each vector entry and sets µ(x j ) (or µ(y i )) to P if the assignment causes the minimum increase in the concerned metric and P is assigned at least one nonzero a ij .

Figure 3 shows the mean execution times of the libraries with the vectorpartitioning heuristic minimizing either MaxSV or TotMsg. We observed that nearly for all libraries and partitioning models, minimizing TotMsg is more important than minimizing MaxSV for reducing the execution time, especially when K is large. The difference is more clear in Fig. 3b. Starting from K = 64, the difference becomes obvious in favor of TotMsg which is concordant with the regression analyses. For K ∈ {8, 16}, minimizing MaxSV or TotMsg seem to be equally important.

Effects of models on metrics

All the partitions used in this work are generated with the objective of minimizing TotVol and balancing the load for a good MaxNnz. To partition x and y, we used the greedy heuristic described above. Since minimizing TotMsg during vector partitioning seems more promising, we will use that variant in the rest of the paper. Here, we analyze the effect of the partitioning models on the metrics. We do not include MaxNnz since it is in the (1 + ) margin of W avg .

The communication metrics vary with respect to the partitioning model as Figure 4 shows for K = {32, 128, 512}. As expected, BL (block) partitioning produces high TotVol and MaxSV (maximum send volume). However, due to the irregular nonzero distributions in original matrices and with the help of the vector partitioning heuristic, its mean TotMsg becomes minimum and its MaxSM (maximum number of messages sent) is reduced to a number less than those of CN (column-net), RN (row-net), and FG (fine-grain). The only model which is better than BL is CB (checker-board) which also results in good volume-related metrics. As expected, its mean MaxSV (maximum send volume) and MaxSM values are much better than the others. The FG model has the worst number of messages and the maximum number of messages. Note that the relative performances of the models do not change with K. However, their difference tend to increase and hence, the model used for partitioning becomes more important as the parallel matrix-vector multiplication times of the libraries show in Figure 5. As regression analysis experiments also prove, the models do not effect the execution time when K is low. However, when K increases, the effect becomes significant. In agreement with Fig. 5, for SpMV and Trilinos, the CB model is the best one. This is expected since CB obtains a good model for TotVol and TotMsg and is the best one for MaxSV and MaxSM. For PETSc, CN is better than BL since CN's volume-based performance is much better while CN and BL are comparable for the latency-based metrics (TotMsg and MaxSM). When K is increased from 128 to 256, the only significant reduction on the execution time is obtained by SpMV with the CB model. This shows that minimizing the right metric with the right partitioning model may be crucial if one wants to use the available resources in their limits.

As Figure 6 shows, for mesh-based matrices, the effects of the partitioning models on TotVol and MaxSV are similar to the ones obtained for real-life matrices. The BL model produces large communication volumes and all the other models can obtain comparable TotVol values. But the relative MaxSV, TotMsg, and MaxSM performance of CB gets worse due to the band structure of the mesh-based matrices which is not suitable for a 2D distribution. For the same reason, the rowwise model BL becomes the best one w.r.t. the latencybased metrics (TotMsg and MaxSM). With BL-based partitioning, each processor sends at most 1 message throughout an SpMxV operation. The runtimes in Fig. 7 are in agreement with these observations. Unlike the experiment with the real-life matrices, the mean execution times significantly decrease till K goes to 512. Note that the average number of nonzeros in most of the real-life matrices and mesh-based matrices are comparable to each other. However, the band structure of the mesh-based matrices yields lower metrics as shown in Figure 6. Hence, the communication overhead for these matrices is smaller compared with that of real-life matrices. We observed that this is the main reason of better scalability.

Related work

Partitioning literature is a rich one. Here we provide a small survey of related work that we believe characterizes major partitioning approaches. Earlier approaches to sparse matrix partitioning concern only the computational load of processors, without paying attention to the effect of sparsity in communication. These approaches are discussed in some recent as well as older studies [START_REF] Bisseling | Parallel iterative solution of sparse linear systems on a transputer network[END_REF][START_REF] Hendrickson | An efficient parallel algorithm for matrix-vector multiplication[END_REF][START_REF] Lewis | Distributed memory matrix-vector multiplication and conjugate gradient algorithms[END_REF][START_REF] Ogielski | Sparse matrix computations on parallel processor arrays[END_REF][START_REF] Pınar | Fast optimal load balancing algorithms for 1D partitioning[END_REF][START_REF] Saule | Load-balancing spatially located computations using rectangular partitions[END_REF]. Even though those methods would yield fairly good load balance, they rely on inherent ordering of matrix rows and columns for reducing the communication cost-which can be arbitrary. To alleviate this problem, first graph-based models have been proposed, later hypergraph-based models have been developed to overcome limitations of graph models [START_REF] Hendrickson | Graph partitioning models for parallel computing[END_REF].

The first use of the hypergraph partitioning methods for efficient parallel sparse matrix-vector multiply operations dates from mid 90s [START_REF] Çatalyürek | A hypergraph model for mapping repeated sparse matrix-vector product computations onto multicomputers[END_REF]. A more comprehensive study [START_REF] Çatalyürek | Hypergraph-partitioning-based decomposition for parallel sparse-matrix vector multiplication[END_REF] describes the use of the row-net and column-net hypergraph models in 1D sparse matrix partitioning. The vector partitioning problem that complements the hypergraph partitioning is a recent one [START_REF] Bisseling | Communication balancing in parallel sparse matrix-vector multiplication[END_REF][START_REF] Uçar | Encapsulating multiple communication-cost metrics in partitioning sparse rectangular matrices for parallel matrix-vector multiplies[END_REF][START_REF] Uçar | Revisiting hypergraph models for sparse matrix partitioning[END_REF].

A fair treatment of parallel sparse matrix-vector multiplication, analysis and investigations on certain matrix types along with the use of hypergraph partitioning is given in [4, Chapter 4]. Further analysis of hypergraph partitioning on some model problems is given in [START_REF] Uçar | On the scalability of hypergraph models for sparse matrix partitioning[END_REF]. Some different methods for sparse matrix partitioning using hypergraphs can be found in [START_REF] Uçar | A matrix partitioning interface to PaToH in MATLAB[END_REF], including jagged-like and checkerboard partitioning methods, and in [START_REF] Vastenhouw | A two-dimensional data distribution method for parallel sparse matrix-vector multiplication[END_REF], the orthogonal recursive bisection approach. A recipe to choose a partitioning method for a given matrix is given in [START_REF] Çatalyürek | On two-dimensional sparse matrix partitioning: Models, methods, and a recipe[END_REF].

All of the studies mentioned above uses static partitioning. A recent work [START_REF] Lee | Adaptive runtime tuning of parallel sparse matrix-vector multiplication on distributed memory systems[END_REF] presents a dynamic partitioning but requires the matrix to be duplicated on all processors; which is infeasible in many applications and hence will not scale.

Most of these studies are more theoretical and comparisons presented on these studies are in terms of one or two partitioning metrics. Our work differs from those by mainly being a through experimental analysis paper. To the best of our knowledge, our paper is the first one that actually compares all such partitioning methods and metrics, using state-of-the-art software packages.

We did not take the performance difference between intra-node and internode communication operations. In recent work [START_REF] Çatalyürek | On shared-memory parallelization of a sparse matrix scaling algoritm[END_REF], Çatalyürek et al. discuss the relevance of hypergraph partitioning methods in shared memory parallelization of a SpMxV-like algorithm. In particular, they show that connectivity-1 metric (which corresponds to the total communication volume in a distributed memory SpMxV) corresponds to the total size of the private arrays used in a possible shared memory implementation. In other words, although there are differences between the metrics optimized by partitioning algorithms for different computing environments, the same partitioning make sense for both implementations, and hence for a hybrid implementation. We did not touch this issue yet but there are on going work by, among others, Byun et al. [START_REF] Byun | pOSKI: Parallel Optimized Sparse Kernel Interface Library User[END_REF] and Yzelman and Roose [START_REF] Yzelman | High-level strategies for parallel sharedmemory sparse matrix-vector multiplication[END_REF].

Conclusion

We have carried out a detailed study to understand the importance of partitioning models and their effects in parallel SpMxV operations. As mentioned in the experiments, minimizing the right metric with the right partitioning model may be crucial to increase throughput. For example, for the real-life matrices in our test set, CB model is the only one which can obtain a significant reduction on the SpMxV time when K is increased from 128 to 256 (after that we did not see any speed up). It is obvious that the other models fail to obtain such a reduction since the gain by dividing MaxNnz by two does not compensate the communication overhead induced by multiplying K by two. Hence, assuming the communication overhead is doubled on the average, doubling K increases the relative importance of communication on SpMxV four times.

Matrices from today's scientific and industrial applications can be huge. If one has only a few processors, partitioning may not matter, since the contribution of communication to the execution time will be low and the overall improvement on SpMxV via a good partitioning will be insignificant. However, as the regression analyses of Section 4.1 show, after a number of processors, the communication overhead will start to dominate the SpMxV time. For our experiments, this number is somewhere between 32 and 64, and it depends on the characteristics of the matrix, the library and the architecture used for SpMxV operations. Although it may be more than 64, considering the advancements on CPU hardware, we can easily argue that this number will remain highly practical and partitioning will matter more for systems that are larger than those considered here.

Algorithm 2 :

 2 ParSpMxV-Overlap-PETSc variant Input: A, x, µ Output: y 1 Send local xj (i.e., µ(xj) = P k ) to those processors that have at least one nonzero in column j 2 Compute y k i ← aij xj for the local nonzeros and local xj, i.e., the nonzeros for which µ(aij) = P k and µ(xj) = P k 3 Receive all non-local xj (i.e., µ(xj) = P k ) 4 Compute y k i ← y k i +aij xj for the local nonzeros and non-local xj, i.e., the nonzeros for which µ(aij) = P k and µ(xj) = P k

Algorithm 3 : 5

 35 ParSpMxV-Overlap Input: A, x, µ Output: y 1 Send local xj (i.e., µ(xj) = P k ) to those processors that have at least one nonzero in column j 2 Compute y k i ← aij xj for the local nonzeros and local xj, i.e., the nonzeros for which µ(aij) = P k and µ(xj) = P k 3 while not all non-local xj received do 4 Receive non-local xj from any processor P (i.e., µ(xj) = P ) Compute y k i ← y k i +aij xj for the local nonzeros and process P 's xj, i.e., the nonzeros for which µ(aij) = P k and µ(xj) = P 6 Send local nonzero partial results y k i to the processor µ(yi) = P k , for all nonzero y k i 7 while not all partial results y i received do 8 Receive nonzero partial results set y from any processor P 9 Compute yi ← y k i +y i for each yi ∈ y RR n°8301

Figure 2 :

 2 Figure 2: Partitioning metrics and corresponding mean execution time for an SpMxV operation in Trilinos on mesh-based matrices when models BL and MP are used.

Figure 3 :

 3 Figure 3: Effect of the vector-partitioning heuristic minimizing MaxSV or TotMsg.The horizontal and vertical axes show the #processors and the geometric means of the execution times (secs) across all matrices in log scale, respectively.

Figure 4 :

 4 Figure 4: Effects of partitioning models on metrics for matrices from real-life applications. The charts show the geometric mean of each metric w.r.t. the partitioning model.

Figure 5 :

 5 Figure 5: Mean SpMxV times on real-life matrices in log scale for each library with respect to partitioning model. The vector-partitioning heuristic in this experiment tries to minimize TotMsg.

Figure 6 :

 6 Figure 6: Effects of partitioning models on metrics for mesh-based matrices. The charts show the geometric mean of each metric w.r.t. partitioning model.

Figure 7 :

 7 Figure 7: Mean SpMxV times on mesh-based matrices in log scale for each library with respect to partitioning model. The vector-partitioning heuristic in this experiment tries to minimize TotMsg.

  N ) is defined as a set of vertices V and a set of nets (hyperedges) N among those vertices. A net n j ∈ N is a subset of vertices, and the vertices in n j are called its pins. A graph is a special hypergraph whose nets have size two. We use pins[n j ] and nets[v i ] to represent the pins of a net n j , and the set of nets that contain vertex v i , respectively. Vertices can be associated with weights, denoted with w[•], and nets can be associated with costs, denoted with c[•].

	3 Summary of the existing partitioning models
	and metrics
	3.1 Hypergraph-based models for SpMxV
	3.1.1 Hypergraph partitioning
	A hypergraph H = (V,

Table 1 :

 1 Properties of the experiment matrices.

	Matrix	Description	n	τ
	atmosmodl	Atmospheric model.	1,489,752	10,319,760
	TSOPF_RS Optimal pow. flow	38,120	16,171,169
	Freescale1	Semiconductor sim.	3,428,755	17,052,626
	rajat31	Circuit sim.	4,690,002	20,316,253
	RM07R	Comp. fluid dyn.	381,689	37,464,962
	cage15	DNA electrophoresis	5,154,859	99,199,551
	HV15R	3D engine fan	2,017,169 283,073,458
	mesh-1024	5-point stencil	1,048,576	5,238,784
	mesh-2048	5-point stencil	4,194,304	20,963,328
	mesh-4096	5-point stencil	16,777,216	83,869,696

Table 2 :

 2 Effect of partitioning metrics TotVol and TotMsg on execution times (in sec) for mesh-based matrices using the MP partitioning model. From one row to the other, among the partitioning metrics only TotVol and TotMsg change.

	mesh-	K	TotVol TotMsg SpMV PETSc Trilinos
	1,024	8	8,200	40	2.24	2.03	2.20
	2,048	32	32,816	184	2.26	2.07	2.49
	4,096	128 131,296	752	2.71	2.45	2.75

Table 3 :

 3 Regression analysis of SpMV, PETSc and Trilinos with all matrices and models CN and BL.

			8 ≤ K ≤ 64		128 ≤ K ≤ 512
	Metric	SpMV PETSc Trilinos SpMV PETSc Trilinos
	MaxNnz	8.02	7.81	6.80	0.49	0.44	0.83
	TotVol	0.18	0.38	1.00	0.39	0.36	1.06
	MaxSV	1.66	1.53	2.20	0.00	0.00	0.11
	TotMsg	0.15	0.28	0.00	7.90	8.03	4.51
	MaxSM	0.00	0.00	0.00	1.22	1.18	3.49

Table 4 :

 4 Regression analysis of SpMV and Trilinos with all matrices and partitioning models. PETSc is not shown in this table because it cannot handle all partitioning schemes.

		8 ≤ K ≤ 32	64 ≤ K ≤ 128	256 ≤ K ≤ 512
	Metric	SpMV Trilinos SpMV Trilinos SpMV Trilinos
	MaxNnz	8.43	7.54	2.75	2.52	0.00	0.02
	TotVol	0.23	0.89	0.52	1.94	0.38	0.98
	MaxSV	1.35	1.57	1.57	1.69	0.04	0.50
	TotMsg	0.00	0.00	4.66	2.38	6.24	3.06
	MaxSM	0.00	0.00	0.49	1.47	3.34	5.44

Table 5 :

 5 Regression analysis of SpMV and Trilinos with mesh-based matrices and all partitioning models.

RR n°8301


Gregory Z Gutin

Royal Holloway

Dimitrios M Thilikos

Parameterized Complexity and the Understanding, Design, and Analysis of Heuristicst Organizers

In the parameterized/multivariate framework, NP-hardness is just the beginning: a result about the null-parameterization. What follows is a rich dialogue between theory and practice, with the goals of:

• explaining the effectiveness of established heuristics and

• designing better heuristics in mathematically disciplined ways

The workshop brought together researchers from both universities and industry, who are interested in exploring how a multivariate view of complexity analysis and algorithm design can lead to industrially useful algorithms in new, mathematically systematic ways.

The vast majority of practical problems are NP-hard and thus we cannot hope to obtain algorithms that compute exact or optimal solutions efficiently, at least according to the one-dimensional (n, the input size) worst-case framework that essentially contrasts worst-case running times that are polynomial in n, versus worst-case running times that are exponential in n.

Recent discussions of the future of theoretical computer science have pointed to the importance of heuristic algorithms for NP-hard problems. In 2009, in his plenary talk at the FOCS Theory Day, Richard Karp posed the question "why do heuristics work so well in practice?" as one of the most important challenges in the theory of computing. Applied computing communities have been aware of this issue for decades, and have lamented the relative mission-failure of theoretical computer science. Theoretical computer science, has reciprocally lamented its lack of support, both by core funding agencies, and in Computer Science curriculums. There is a central dysfunction between theory and practice to which Karp points.

In real life, typical inputs to a computational problem are the outputs of other resource-constrained processes, both natural and artificial. While datasets may be large, significant (although often multifaceted) structure is almost always present and worth taking into account in algorithm design and complexity analysis.

Parameterized complexity attends to the underlying structure of real-world input distributions, as well as structure in the computational objectives, such as size of solutions, noise in the data, goodness of the desired approximation, or combinations of these. (See [START_REF] Downey | Parameterized complexity[END_REF][START_REF] Flum | Parameterized Complexity Theory[END_REF][START_REF] Niedermeier | An Invitation to Parameterized Algorithms[END_REF].) Secondary measurements (that may form an aggregate parameter) capture such crucially relevant problem structure and open up the possibility of exploiting this structure in order to design efficient worstcase algorithms -or to identify algorithmic strategies that can usefully be incorportated into practical heuristics. The central notion of parameterized complexity is a multivariate generalization of polynomial time that is termed fixed-parameter tractability (FPT), meaning solvability in time f (k)+n c , where f is some function of the parameter k, n is the input size, and c is a constantin other words, the input size contributes polynomially to the complexity, and any exponential costs are confined to some function of the parameter k.

Parameterized algorithms, both directly and indirectly, have made major contributions to practical computing in such areas as Computational Biology and (increasingly) Artificial Intelligence. The historical record is muddy. Since computing began, applications-determined implementors have of course attended to relevant secondary measurements that allowed their implementations to succeed on real-world data-sets (but without any attention to theoretical frameworks). Meanwhile, onedimensionally framed and trained theorists have crept towards a multivariate outlook, under the rhetoric of exact algorithmics.

Although the historical record is muddied, it is essentially the case that the modern computational revolutions in Biology and Medicine would not have happened without parameterized/multivariate algorithmics. Explicit mathematical formulation of the issues adds power to the algorithm engineering opportunities, not only in these application domains, but wherever NP-hardness of central computational tasks has been identified, and heuristics are currently the main game.

Real-life input distributions have significant structure, but what is it, and how should it be explicitly declared in realistic multivariate complexity bookkeeping? The multivariate approach to complexity analysis and algorithm design often demands significant dialog between theory results, and practical computing experience, including the adaptation of FPT algorithms and mathematically-centered algorithmic strategies into practical-computing-targeted heuristics:

• It is well-known that the input of every FPT parameterized problem can be pre-processed, in polynomial time, such that the reduced instance (losing no information) has size bounded (only) by a function of the parameter (for many NP-hard problems, a polynomial function). Such kernelization subroutines can be highly useful in practical computing, regardless of how the reduced instances are solved. For the first time: a useful mathematical theory of pre-processing.

• Parameterization also offers, for the first time, a useful theory of the complexity of local search subroutines.

• FPT subroutines are also useful in the "genetic recombination" phase of hybrid (e.g., memetic) meta-heuristics.

• Algorithmic strategies for FPT algorithms, such as iterative compression and greedy localization, can be systematically adapted into theoretically organized classes of heuristics (see [Karp11]) that are competitive or superior to previous efforts, or in some cases, exactly the same: the successful heuristic was simply an unrecognized FPT algorithm, according to an unrecognized but relevant parameter (see examples in [DF99]).

• The mathematical determination that a parameterized problem is FPT often involves estimates concerning mathematical relationships that are far from the plausible but unproven truth of the matter, or worst-case assumptions about unparameterized aspects of the input that are far from typical. Systematic algorithm engineering compromises about such issues can lead to novel and effective heuristics derived from FPT algorithms (see [TM09] for a pioneering example).

• The central innovation of the parameterized approach to algorithms and complexity is, of course, simply the explicitly multivariate view of complexity analysis and algorithm design. Beyond that, at the core of the multivariate algorithmics project, there is a richer vision of the workflow between theory and practice. In the classical framework of P versus NP, a problem is NP-hard, and that is almost the end of the story. As Karp wrote in [Karp11]:

"Except in rare cases, the problems are NP-hard, and the performance guarantees provided by polynomial-time approximation algorithms are far too pessimistic to be useful."

In the multivariate workflow there are two complementary creative principles:

Deconstruction of hardness results. Most hardness proofs in both classical and parameterized complexity are unreasonable. Read the proofs and ask why the constructions are unreasonable for real-world inputs. This often points to realistic FPT parameterizations.

Model enrichment. When FPT results are obtained, it is often productive to enrich the definition of the problem, in order to obtain better traction in real-world applications. Parameters can be complex aggregates that wrap onto problem realism.

Because a single classical problem can be parameterized in an unlimited number of ways, this rich mathematical dialog exposes many opportunities for fresh approaches in the design of practical algorithms and heuristics. The goal of the Workshop was to explore, systematize and programmatically articulate these opportunities.

Overview of Talks

Practical Aspects of Fixed-Parameter Algorithms: An Implementations' Viewpoint Faisal Abu Khzam, Lebanese American University The talk will cover some implementation techniques and a discussion of design alternatives, including a comparison between best-worst-case algorithms and algorithms based on simple heuristic priorities. Some parallel implementation strategies will be discussed as well.

The Matrix Representation with Flipping problem and the FlipCut heuristic

Sebastian Boecker, University of Jena, Germany In computational phylogenetics, supertree methods assemble phylogenetic trees with non-identical but overlapping taxon sets, into a larger supertree: These supertrees contains all taxa of all input trees and describes the evolutionary relationship of these taxa. This problem can be formalized in different ways, to cope with contradictory information in the input. I will concentrate on the Matrix Representation with Flipping (MRF) framework: The input trees are encoded in a 0/1/?-matrix, and we search for a minimum set of 0/1-flips such that the resulting matrix admits a directed perfect phylogeny. This is an NP-hard problem, but if all input trees share the same set of taxa, the problem is fixed-parameter tractable. Otherwise, the problem is W[2]-hard and cannot be approximated within a constant factor. I will present a data reduction as well as an ILP formulation for the problem which, unfortunately, both do not work well in practice. Finally, I will present a novel heuristic for the problem, named FlipCut supertrees. Evaluations of the heuristic are very promising.

Cluster Editing in Practice

Sebastian Boecker, University of Jena, Germany Cluster Editing ask to transform an undirected graph into a disjoint union of cliques using a minimum number k of edge modifications. I will talk about practical issues that concern the implementation of algorithms for Cluster Editing, including the growth of the parameter k, running times in practice, bestperforming reduction rules, and others.

Combinatorial optimization using diversities

David Bryant, University of Otago A diversity is a generalisation of metric spaces to beyond pairwise comparisons. Formally, it is a pair (X, δ) where X is a set and δ is defined on finite subsets of X satisfying the two axioms

1. δ(A) = 0 ⇔ |A| ≤ 1; 2. δ(A ∪ C) ≤ δ(A ∪ B) + δ(B ∪ C) whenever B = ∅.
Note that δ restricted to pairs gives a metric. Examples of diversities include the diameter, the length of the minimum tour, and the length of an optimal Steiner tree. Diversities were introduced in a recent paper of Bryant and Tupper, and generate a rich theory which we are only just beginning to explore. We survey some of the main results and then speculate wildly on their use in applied combinatorial optimization and data analysis. In particular we discuss their role in hypergraph algorithms and show how L 1 embedding strategies for graphs and metrics generalize to analogous results for hypergraphs and diversities (at least that is the plan).

Exploring parameterised complexity in computational topology

Benjamin Burton, The University of Queensland Topological decision problems in three dimensions are notoriously difficult: the mere existence of an algorithm is often a major result, and many important algorithms have yet to be studied from a complexity viewpoint. Even "simple" problems, such as unknot recognition (testing whether a loop of string is knotted), or 3-sphere recognition (testing whether a triangulated 3-manifold is topologically trivial), have best-known algorithms that are worst-case exponential time. In practice, however, some of these algorithms run surprisingly well in experimental settings, and we discuss some reasons why parameterised complexity now looks to be the "right" tool to explain this behaviour. We also present some initial forays into the parameterised complexity of topological problems, including both fixed-parameter-tractability and W[P]-hardness results.

What have I been doing recently in parameterized complexity

Rod Downey, Victoria University of Wellington I will mention some of the work in this area I have been looking at. This includes the new book with Mike, plus lots of new projects with Judith Egan, Mike, Fran and Peter Shaw. These concern permissive search, incremental computation, hueristics, and parameterizing by the number of numbers. I will give a brief overview of these.

Inferring Haplotypes from Genotypes on a Pedigree

Patricia Evans, University of New Brunswick In genetics, it is considerably easier to determine genotype information (what values an individual has at each of a number of marker sites) than it is to determine haplotype information (how those values are split between the chromosomes). Since haplotype information can be very useful in pedigree analysis, we would like to be able to infer haplotypes from genotypes for a set of individuals with known family relationships in a pedigree. This work examines how haplotypes including recombination can be inferred as an application of Edge Bipartitization. This is joint work with Peter Doan.

Parameterized k-opt

Fedor V. Fomin, Univ. of Bergen We discuss results on the parameterized complexity about searching in the k-exchange neighborhoods for the following optimization problems -Problems on planar graphs -Feedback Arc Set in Tournaments (FAST) -Chordal completion (fill-in)

Domination Analysis of Heuristics

Gregory Gutin, Royal Holloway, University of London Many heuristics do not provide approximation guarantee in the worst case and so we cannot compare worst case performance of heuristics for the same problem using approximation. Domination analysis is an alternative way of comparing heuristics based on the concept of domination number proposed in 1996 by Glover and Punnen. (Interestingly, some results in the area were obtained already in the 1970s.) Consider the Asymmetric Travelling Salesman problem (defined on weighted complete directed graph). The domination number of an ATSP heuristic H donated dom(H, n) is the maximum integer d such that for every instance of ATSP on n vertices H finds a solution not worse than d solutions. For example, for every n > 1 there are instances on n vertices for which the nearest neighbour heuristic (NN) will find the unique worst possible solution. Thus, the domination number of NN is 1. There are ATSP heuristics with domination number at least (n -2)! and even higher. We'll concentrate mainly on TSP heuristics, but heuristics for other problems will also be mentioned.

Exponential Neighborhoods

Gregory Gutin, Royal Holloway University of London Exponential neighbourhoods contain exponential number of solutions. However, there are some such neighbourhoods for which the best solution can be found in polynomial time. We'll consider such neighbourhoods for TSP and other problems. We should be able to get something much bigger if we replace 'polynomial time' with 'FPT time', but for which problems we can do it and what is the parameter?

Implementing and testing fixed-parameter algorithms

Falk Hüfner, Technische Universität Berlin

The talk is about various aspects of implementing and testing fixed-parameter algorithms, such as ensuring correctness, debugging, speedups, and performance evaluation and visualization.

Parameterized Testability

Kazuo Iwama, School of Informatics, Kyoto University In graph property testing, given a graph represented by an oracle, we want to test whether the graph satisfies some predetermined property P or ǫ-far from satisfying the property. Roughly speaking, a graph is called ǫ-far from a property if we need to modify an ǫ-fraction of the graph to make it satisfy the property. A most important focus in this topic has been on general characterizations for (constant-time) testable properties, which have obtained much progress in the last decade: For the dense graph model, Alon et al. finally obtained a purely combinatorial necessary and sufficient condition for P to be testable that is closely related to the Szemerédi Regularity Lemma. For the boundeddegree graph model, due to Benjamini et al., any property is testable if it is minor-closed. In this talk, we are still interested in a general framework of testable properties, but our approach is a bit different from the qualitative one as above. Rather, we are interested in "parameterized properties" that are quite common in NP optimization problems. A bit curiously, this direction has not been so popular in the context of property testing. The obvious reason is that the problem becomes less interesting in both the dense graph model and the bounded degree model. For instance, consider the problem k-Vertex Cover for a constant k. In the dense graph model, we want to decide whether a graph has a vertex cover of size at most k or we need to remove at least ǫn 2 edges to have the property. It turns out that this property is easily testable just by accepting only graphs with at most ǫ × n 2 edges. To avoid this triviality, we consider the general graph model with an augmentation of random edge sampling capability. It is shown that a variety of such problems, including k-Vertex Cover, k-Feedback Vertex Set, k-Multicut, k-Path Free and k-Dominating Set, are constant-time testable if k is constant. It should be noted that the first four problems are fixed parameter tractable (FPT) and it turns out that algorithmic techniques for their FPT algorithms (boundedbranching tree search, color coding, etc.) are also useful for our testers. k-Dominating Set is W[2]-hard, but we can still test the property in constant time since the definition of ǫ-farness makes the problem trivial for non-sparse graphs that are the source of hardness for the original optimization problem. We also consider k-Odd Cycle Transversal, which is another well-known FPT problem, but we only give a sublinear-time tester when k is a constant. This is a joint work with Yuichi Yoshida.

LP-Based Parameterized Algorithms for Separation Problems

Daniel Lokshtanov, Univ. Bergen We show that the most basic branching algorithm for Vertex Cover turboboosted by a classical linear-programming based reduction rule runs in time 4 k-OP T (LP ) poly(n), where OP T (LP ) is the value of the optimum solution to the LP relaxation of the problem. It turns out that this algorithm can be used to solve Odd Cycle Transversal and Almost 2-Sat, as well as a number of other problems, in time 4 k poly(n). We then give an improved algorithm for Vertex Cover with running time 2.32 k-OP T (LP ) poly(n), this improves considerably over the currently fastest parameterized algorithms for Odd Cycle Transversal and Almost 2-SAT.

Linear kernel via conflict packing -application to FAST and dense RTI problems

Christophe Paul, CNRS-LIRMM We introduce a new technique to design kernelization algorithms for parameterized problems, namely Conflict Packing. We illustrate this technique on two well-studied problems: FAST and RTI. For the former, one is given a tournament T = (V, A) and seeks a set of at most k arcs whose reversal in T leads to an acyclic tournament. While a linear vertex-kernel is already known for this problem, it requires a constant-factor approximation algorithm as a subroutine. Using the Conflict Packing, we show how to avoid this subroutine to find a socalled safe partition. In the RTI problem, one is given a set of vertices V and a dense collection R of rooted binary trees over three vertices of V and seeks a rooted tree over V containing all but at most k triplets from R. Using again the Conflict Packing, we prove that the dense RTI problem admits a linear vertexkernel. This result improves the on the quadratic known kernel. joint work with A. Perez and S. Thomassé

Implementing Courcelle's Theorem

Peter Rossmanith I will make a demonstration of our MSO solver and explain how it works. The solver is based on game trees.

Tractable Parameterizations for the Minimum Linear Arrangement Problem

Hadas Shachnai, Technion, Israel

The Minimum Linear Arrangement (MLA) problem asks to embed a given graph on the integer line so that the sum of the edge lengths of the embedded graph is minimized. Most layout problems are either intractable, or not known to be tractable, parameterized by the treewidth of the input graphs. We investigate MLA with respect to three parameters that provide more structure than treewidth. In particular, we give a factor 1 + ǫ -approximation algorithm for MLA parameterized by (ǫ, k), where k is the vertex cover number of the input graph. By a similar approach, we describe two FPT algorithms that exactly solve MLA parameterized by, respectively, the max-leaf and clique-cover numbers of the input graph. Joint work with Michael Fellows, Danny Hermelin and Frances Rosamond.

Solving Hard Problems Incrementaly

Peter Shaw, Charles Darwin University By utilizing an incremental (or solution repair) formulation, some hard problems, when formulated with appropriate parameters, become FPT. This talk describes how algorithmic solutions of such problems can be applied as effective FPT subroutines to improve many local search problems. To illustrate this we consider the Dominating Set problem. We show that the incremental Dominating set problem is FPT. A number of formulations of this problem were considered, but most of these resulted in hard problems. We show a FPT algorithm for one particular formulation and further prove that this problem has no polynomial kernel. This is a joint work with Rod Downey, Judith Egan, Michael Fellows, and Frances Rosamond.

Some Algorithmic Challenges in Arc Routing

Manuel Sorge, Technische Universität Berlin Arc Routing problems occur naturally in many logistic tasks like delivering mail, garbage disposal, street sweeping and snow plowing. Such a task is formulated, for example, in the Chinese Postman problem where, given an edgeweighted graph, one searches for a minimum-weight tour that visits all edges in the graph. In this simple form, Chinese Postman is solvable in polynomial time, but many variants that arise in applications yield NP-hard problems. Yet, despite the practical relevance, the literature on parameterized complexity of Arc Routing problems is rather sparse. We give a brief survey of the known (parameterized) complexity results and identify some interesting open questions. In particular, we briefly report on our theoretical and empirical studies related to the Rural Postman problem. This is joint work with René van Bevern, Rolf Niedermeier, and Mathias Weller.

Deconstructing Algorithms: From Theory to Engineering

Siamak Tazari, Google

We present the generic idea of deconstructing a sophisticated theoretical algorithm, and engineering its parts to obtain a useful, practical "heuristic" via a case study on Steiner tree in planar graphs. Afterwards, we will briefly talk about the disconnect between theoretical computer science and the state of the industry today and propose possible directions on bringing them closer together.

From parameterized complexity to heuristics

Dimitrios M. Thilikos, CNRS-LIRMM & University of Athens

In this talk we describe some of the the existing links netween parameterized complexity and heuristics. More particular we present how the notion of kernelization can serve as a formal mathematical framework for studying preprocessing and we explain why the multivariate computational complexity viewpoint is necessary for this approach. We also give several research directions for the study of various heuristic approaches using tools and techniqeus from parameterized algorithm design.

FPT for Scalable Mining Algorithms

Takeaki Uno, National Institute of Informatics Data mining algorithms often enumerates patterns/structures included in the given data, for exhaust search of candidates. They usually deal with data that are noisy but have some property (sparse, locally dense, having hubs, etc.) The time complexity of the algorithm are usually exponential in the input size, but practically they are efficient. One of the reason is that, the number of solutions is not exponential, relatively small, and the algorithms are output sensitive (# solutions (resp., # edges) is a kind of "fixed parameter"). We would like to introduce fixed parameters to data mining/enumeration algorithms so that we can solve difficult problems, and can understand the mechanism of the practical efficiency.

Efficient Enumeration of the Directed Binary Perfect Phylogenies from Incomplete Data

Yoshio Okamoto, Japan Advanced Institute of Science And Technology

We study a character-based phylogeny reconstruction problem when an incomplete set of data is given. More specifically, we consider the situation under the directed perfect phylogeny assumption with binary characters in which for some species the states of some characters are missing. Our main object is to give an efficient algorithm to enumerate (or list) all perfect phylogenies that can be obtained when the missing entries are completed. While a simple branch-and-bound algorithm (B&B) shows a theoretically good performance, we propose another approach based on a zero-suppressed binary decision diagram (ZDD). Experimental results on randomly generated data exhibit that the ZDD approach outperforms B&B. We also prove that counting the number of phylogenetic trees consistent with a given data is #P-complete, thus providing an evidence that an efficient random sampling seems hard. Joint work with Masashi Kiyomi and Toshiki Saitoh.

Graph Isomorphism Problem Parameterized by Width Parameters

Yota Otachi, Japan Advanced Institute of Science And Technology

We study the parameterized complexity of the graph isomorphism problem

when parameterized by width parameters related to strong tree decompositions. We first present FPT algorithms for some width parameters, and then provide a unified explanation for various isomorphism results concerned with parameters related to tree decompositions. As a first step towards intractability results for parameterized graph isomorphism we develop an FPT Turing-reduction from strong tree width to the a priori unrelated parameter maximum degree. This is joint work with Pascal Schweitzer.

List of participants