Towards a multigrid method for the minimumcost flow problem
Abstract.
We present a first step towards a multigrid method for solving the mincost flow problem. Specifically, we present a strategy that takes advantage of existing blackbox fast iterative linear solvers, i.e. algebraic multigrid methods. We show with standard benchmarks that, while less competitive than combinatorial techniques on small problems that run on a single core, our approach scales well with problem size, complexity, and number of processors, allowing for tackling largescale problems on modern parallel architectures. Our approach is based on combining interiorpoint with multigrid methods for solving the nonlinear KKT equations via Newton’s method. However, the Jacobian matrix arising in the Newton iteration is indefinite and its condition number cannot be expected to be bounded. In fact, the eigenvalues of the Jacobian can both vanish and blow up near the solution, leading to a significant slowdown of the convergence speed of iterative solvers  or to the loss of convergence at all. In order to allow for the application of multigrid methods, which have been originally designed for elliptic problems, we furthermore show that the occurring Jacobian can be interpreted as the stiffness matrix of a mixed formulation of the weighted graph Laplacian of the network, whose metric depends on the slack variables and the multipliers of the inequality constraints. Together with our regularization, this allows for the application of a blackbox algebraic multigrid method on the Schurcomplement of the system.
Key words and phrases:
Multigrid methods, interiorpoint methods, parallel algorithms.2010 Mathematics Subject Classification:
65M55, 90C51, 68W10.1. Introduction
This work is concerned with the solution of the minimumcost flow (MCF) problem
(1) 
where and are known as the cost and flow vectors, is a incidence matrix of a network with nodes and arcs, and the constraints have the form:
(2) 
where and are vectors of, respectively, lower and upper bounds on the sought flow . There is a long history concerning the development of combinatorial techniques for the MCF problem. Several classes of methods have appeared over the years, starting from the seminal work by Ford and Fulkerson [14]. Among popular approaches, we recall cyclecanceling algorithms [25, 19], successive shortest paths [12, 35], capacityscaling [12, 32], costscaling [18, 20], and the network simplex [1, 23]. Unfortunately, none of these methods has been shown to be of optimal complexity [1, 26], although they usually behave much better in practice than their theoretical (worst case) estimates may predict.
Although already with the ellipsoid method a solution method with a worst case polynimial complexity existed, the “projective algorithm” introduced by Karmarkar [24] not only provided a better estimate ( instead of (), but also worked well in practice. The original algorithm by Karmarkar has been reformulated later as interior point method (IPM). More elaborated variants of interior point methods, i.e. primaldual methods, nowadays are also applied routinely to MCF problems [8, 9, 10, 11, 16, 31, 33, 34]. In the framework of IPMs, the main computational burden consists in the solution of a linear system of equations, involving the symmetric but possibly indefinite Jacobian arising within Newton’s method applied to the KKT equations from the IPM.
As it turns out, there is no clear “best method”, as the choice of the “best” method will depend on the considered application and the targeted problem size. The practical performance of solution methods for the MCF problem is usually evaluated using benchmarks, as described in [17, 26].
While combinatorial techniques are very efficient on problems of moderate size that can run on a single processor, they are in general not trivial to parallelize. Thus, in terms of large scale problems, at a first glance IPMs seem to be more attractive, since parallel linear solvers for large scale problems are available. However in the case of the IPM, the arising Jacobian matrices can become extremely illconditioned in the neighborhood of the solution. For this reason, the most commonly studied approach is to apply a direct solver, such as the Cholesky factorization [30]. Even with this choice, however, problem (1) needs to be regularized before a factorization can be computed. A more detailed discussion can be found in, e.g., in [15, 21], where two related diagonal perturbations of the system are proposed.
While direct solvers in general might be considered to be more robust, their nonoptimal complexity leads to similar shortcomings to those that plague the combinatorial approaches, in particular the suboptimal scaling with respect to the number of arcs . This fact motivated several authors to consider Krylov subspace methods [8, 16, 31, 33, 34] to solve the MCF problem. However, the efficiency of this approach relies on the spectral properties of the Jacobian matrix. Therefore, most of the works in this class focus on preconditioners for this matrix. While using the diagonal of the Jacobian proved to be very effective in the first iterations, the situation is more complex close to the solution, where a spanningtree preconditioner needs to be applied. In [33, 34] the use of the conjugated gradient method is proposed, employing a heuristic to select between the two options. However, it remains an open question if the complexity of this approach is optimal with respect to the problem size.
Outside the scope of the MCF problem, linear solvers of optimal complexity, such as multigrid (MG) methods [4, 5, 6], are traditionally applied in the context of linear elliptic operators. In extension, multigrid methods can also be applied to constrained minimization problems, i.e. the minimization of convex quadratic functionals subject to box constraints, as it occurs in contact problems [27], or to the IPM, see [11], where the Schur’s complement method is used in presence of additional equality constraints [29]. However, the application of multigrid methods to the Jacobian matrices arising in the IPM is not straightforward. This is caused by the particular structure of the MCF, i.e. the combination of linear and inequality constraints with a linear objective function, leads to several difficulties. In this paper, we discuss these difficulties and present easily implementable and robust solutions. We deal in particular with the following the following aspects:

As discussed in Section 3, the Jacobian of the KKT equations for the MCF problem is indefinite rather than positive definite. Fortunately, it turns out that this matrix is formally equivalent to the discretized operators of the mixed formulation of the Laplace operator in finite elements. As the latter is elliptic, this creates the bridge to the application of multigrid methods.

The Schur complement of the system, encoding the actual graph Laplacian operator of the graph, is applied to the Lagrangian multipliers of the equality constraints, rather than to the flow vector ; thus, box constraints on the solution become general inequality constraints for this variable, hindering the application of standard techniques for constrained elliptic problems, such as, e.g., projected GaussSeidel [27]. Here, we develop and present alternatives.

The Schur complement is plagued by the same illconditioning as the full Jacobian, therefore even MG solvers developed for the graph Laplacian [10, 28] fail to converge in the proximity of the true solution. Hence appropriate regularization techniques need to be studied. Our strategy is presented in Section 4.
Our goal is to fill the gap between the study of linear solvers for the IPM applied to the MCF and the optimalcomplexity methods for elliptic problems. More specifically, we decide to feed the Schur’s complement of the Jacobian of the KKT equations directly into a blackbox algebraic MG (AMG) preconditioner, such as, e.g., the BoomerAMG algorithm from [13]. This in principle could extend the class of problems that can be tackled efficiently using solution methods developed for elliptic problems and improve the IPM complexity for the MCF problem with respect to the number of arcs. However, as explained above, a direct application of (A)MG will lead to several difficulties. In order to overcome those, here we develop appropriate regularizations of the problem that yield a lower and upper bounds on the eigenvalues of the graph Laplacian. Our strategy is to consider the MCF problem as the task of finding the static equilibrium of a dynamic contact problem, by adding damping to the system and combining techniques from the IPM and activeset methods (ASM). While we focus on the MCF problem, we argue that these techniques might help dealing with the illconditioning and nonconvexity of Jacobians in more general applications of the IPM.
We here follow the path of regularizing the problem, combining a diagonal mass matrix arising from a pseudokinetic energy term added to (1) and the activeset method (ASM) to handle inequality constraints, as detailed in Section 4.
1.1. Outline
The paper is organized as follows:

In Section 2, we present the mathematical formulation of the IPM applied to the MCF problem.

In Section 3, we outline the basics of linear solvers applied to the MCF problem, highlighting why outofthebox MG methods cannot be used directly.

In Section 4, we discuss our modifications to the standard IPM to solve the KKT equations, which produce a positive definite Schur’s complement.

In Section 5, we show some numerical results on standard benchmarks.
2. Solving the MCF problem via IPM
The MCF problem (1) yields the following first order optimality conditions
(3)  
(4)  
(5)  
(6) 
where and are the Lagrange multipliers associated to the constraints and , respectively, and where we have denoted the gradient of as
(7) 
For solving the nonlinear KKT conditions (3), usually Newton’s Method is used. The main computational effort in each Newton iteration is the solution of the linear system
(8) 
where and are a diagonal matrices having, respectively, and on their main diagonals. The righthandside is evaluated based on the current Newton iterate, thereby measuring the actual residual of the solution.
In order to find a reasonable stepsize for the Newton correction, line search has to be performed to find and such that
(9)  
(10) 
In general, this straight forward application of Newton’s method will lead to line search increments and , which will be too small for making this strategy practical. To solve this issue, often reformulations based on interior point methods (IPM) are employed, see, e.g. [30]. In IPM, problem (1) is modified with a barrier potential
(11) 
where is a centering parameter to be specified. By defining , this formulation yields a modified third KKT optimality condition via the substitution
(12) 
Clearly, as , these equations tend to the KKT conditions for the original problem (1). Among the several known strategies to combine accuracy (small ) with a small number of Newton iterations (large ), we decided to employ Mehrotra’s predictorcorrector algorithm outlined in [30]. While this choice is obviously well suited for direct solvers, as it requires to solve the system twice with the same matrix, we also obtain a benefit with AMG methods, since the setup phase of the AMG has to be executed only once per Newton iteration.In our experiments, this proved to be a small benefit for graphs with a regular structure, e.g. graphs which would arise from discretizations of the Laplacian on structured meshes, but yielded a significant speedup on complicated graphs, such as the NETGEN data set discussed in the Section on numerical experiments.
3. Mixed formulation and its limitations
In order to be able apply multigrid methods, the Newton iteration (8) must be transformed into a positive definite system. First, let us observe that the last row can be solved for, reducing the system to
(13) 
where . This system is indefinite but exhibits a special structure that, in the context of FE, is known as the mixed formulation of the graph Laplacian induced by the metric . We here show more precisely this equivalence, highlighting that the distortion of close to the solution renders the IPM challenging for linear solvers. We start with the mixed formulation of the Laplace equation
where and are Hilbert spaces. In FE, subspaces and are chosen to approximate the solution. In particular, we select piecewise linear functions with, respectively, vertices and edges of the graph as degreesoffreedom, and we look for a solution . This results in a system of linear equations
The matrix represents a gradient operator, using the fact that the gradient of a linear basis function is constant over the edges, we can write it as , thereby separating a purely combinatoric matrix from the weights matrix , containing the contributions depending on the metric (i.e. the embedding of the graph into an ambient space). This yields
With the substitution , we obtain
whose matrix is identical to our system if . Therefore, the Schur’s complement of (13) is the actual Laplace operator of the graph embedded with the metric
(14) 
Since this operator is elliptic and the resulting matrix has a bounded bandwidth, it is wellknown that the MG methods have optimal complexity. Unfortunately, from the definition of it becomes evident that there are three issues that hinder the straightforward application of linear solvers:

If any of the components of tends to zero, then at least one eigenvalue of blows up.

If any of the components of tends to zero, then at least one eigenvalue of tends to zero.

If the gradient of the active constraints (including the equalities) is not fullrank, then the Schur’s complement becomes singular. This would be equivalent to violating the infsup condition in FE, which results in a nonunique solution for the dual variable .
Note that while issue 3 is problemspecific and may not happen in general, 1 and 2 need to occur near the solution. In the following, we will propose a solution for each one of the above issues, allowing for a direct application of a blackbox AMG solver to the Schur’s complement .
4. Solving the three limitations
As we will see, the main intuition behind the proposed regularization is to treat the MCF as the static equilibrium of a dynamic contact problem. While MG methods are well studied in this context [27], the presence of the linear constraints introduces a coupling between bound constraints on the variable , making techniques such as projected GaussSeidel [27] inapplicable. To see this, recall the fact that projected GS needs to project back to the feasible region one component of the solution at the time. However, box constraints on induce general inequality constraints on
(15) 
Clearly, given a that violates this constraint, it is not trivial nor efficient to correct it while satisfying simultaneously the linear equalities. Therefore, we will introduce appropriate boundary conditions on the variable applied to the system before the Schur’s complement is computed (Section 4.2). These, in turn, will make the system singular and therefore will require additional boundary conditions on the dual variable in order to remove the kernel from the system (Section 4.3). We argue that, while the dual boundary conditions are very specific to the MCF problem, the first two ideas can be applied as well to general nonlinear problems exhibiting a saddlepoint structure with box constraints.
4.1. Solving 1: Pseudodynamic minimization
The first issue stems from the fact that the Hessian of (11) with respect to vanishes. Other authors have recognized this problem and have proposed a similar regularization [15, 21]. However, we here give a physical interpretation of the perturbed problem as the problem of finding the static equilibrium of a dynamic problem, opening up the possibility to tune the method for a faster convergence through an appropriate amount of damping added to the system. We start by transforming the functional in (11) to the Lagrangian
(16) 
where is an arbitrary mass matrix. In practice, we found that the choice works reasonably well. By using a damping parameter and the explicit linear approximation in time , we define an unconstrained step as
(17) 
which will serve as an initial guess of the solution after the Newton step. By doing so, we can use a stepandproject approach [22] to project back the unconstrained solution to the manifold defined by our constraints. The projection problem consists of finding the stationary point of the functional
(18) 
By computing the stationary equations, we obtain a modified first KKT condition
(19) 
which results in the regularized version of (13):
(20) 
It is clear to see that the Schur’s complement is
(21) 
whose eigenvalues are bounded below by as approaches zero. It can also be noted that for the special choice , it follows that , so the respective righthand sides of (13) and (20) would coincide. This would be equivalent to the common choice of only modifying the matrix [15]. However, it is clear that there is a tradeoff at play: bigger and smaller imply a better spectrum for , but also more time steps to reach convergence.
In order to improve convergence, we performed the adaptive adjustment shown in Algorithm 1 with the choice , aimed at providing a lower bound on the eigenvalues of close to the solution, while taking larger steps during the first interations. We did not study further the problem of minimizing the total number of linear solver iterations required by all Newton steps. In fact, we found that for all benchmarks except the GOTO data set, we do not need to use Algorithm 1. In these cases, the time step and the damping parameter can be held fixed, confirming the robustness of the proposed approach.
4.2. Solving 2: Activeset methods (ASM)
As discussed earlier, standard techniques for contact problems capable to handle box constraints are not straightforward to be applied within the MG iterations, due to the presence of linear constraints. However, a similar approach can be taken in between successive Newton steps. At the Newton step , we check the condition
(22) 
If satisfied, it implies that some of the variable are near to contact the boundary of the feasible region and we deem the corresponding constraints to be active, i.e. , , where is the set of active constraints. This translates to applying Dirichlet boundary conditions to the variable in the systems (13) and (20). This is achieved via the substitutions
(23)  
(24)  
(25) 
Moreover, we test for the dual condition
(26) 
This is the deactivation criterion, expressing the fact that if contact forces are not strong enough, then the corresponding constraints are in fact not active. In practice, we found this not to be necessary for the MCF problem and therefore we did not have to deactivate any constraint in the benchmark presented in the next section. With such a treatment of the inequality constraints, our algorithm can be viewed as an activeset method, where rather than testing for a large set of constraints, we use the IPM to find an appropriate guess of the active constraints at the exact solution.
4.3. Solving 3: Connected components of a graph
While solutions 1 and 2 are general and can potentially be applied to nonlinear problems, the solution to the third issue exploits the particular nature of the MCF problem. Let us recall the fact that is a discrete gradient operator, it follows that the Lagrange multipliers are a nodebased function whose edgebased gradient is equal to . Therefore, the function is defined up to a constant. When inequality constraints become active, however, the incidence matrix is modified in a way equivalent to removing the edges corresponding to from the original graph. Recall the following
Definition 1.
A weakly connected component is a maximal subgraph of a directed graph such that for every pair of vertices (u, v) in the subgraph, there is an undirected path from u to v and a directed path from v to u.
It follows that the multipliers are defined up to a constant per weakly connected component of the graph. While graphbased algorithms for computing connected components in linear time do exists, we opt instead for taking advantage of parallel matrixvector multiplication routines. To achieve this, we assemble the adjacency matrix and perform the following fixedpoint iteration:
(27) 
where is initialized as . The fixed point of this iteration is a vector whose entries sharing the same label belong to the same connected component. Then, we pick one label per connected component and set
(28)  
(29) 
While not competitive with graphbased algorithms on a single processor, this approach showed to scale very well with the number of cores and is only needed in the close proximity of the true solution, when constraints become active. Hence, further improvements to the run times presented in the next section are possible with a graphbased algorithm, but not to the scaling, which is the primary concern of this work.
5. Numerical experiments
We tested the proposed method on the following data sets [26]:

ROAD_PATHS: These instances are road networks from different states in the USA, where the edge cost is set to the travel time along the edge. There are approximately randomly selected sources and sinks each with a demand that depends on the maximum flow that can be sent between these sources and sinks.

ROAD_FLOWS: Same as above but the capacities are not set to 1, but depending on the category of the road to either 40, 60, 80 or 100.

NETGEN_8: They are generated using the NETGEN random generator. The subscript 8 indicates that the networks are sparse, i.e., and the capacities and costs are chosen uniformly at random between 1 and 1000 and 1 and 10000, respectively. There are approximately sinks and sources and the total demand is .

GOTO_8: They are generated using the GOTO generator. They are grid instances on tori, known to be rather hard. The subscript 8 indicates that the networks are sparse, i.e., and the capacities and costs are chosen uniformly at random between 1 and 1000 and 1 and 10000, respectively. There is one source and one sink node and the supply is adjusted to the capacity.

Signal processing: These instances are regular quadrilateral grids taken from satellite data that can be measured only up to mod . Costs and flows are specified in order to reconstruct the true values of the function over the grid [7].
For our implementation, we rely on BICGStab from PETSc [2] as the linear solver and the BoomerAMG algebraic MG preconditioner, which is part of the Hypre package [13]. The default Hypre options and parameters have been used in all data sets, with the exception of NETGEN and GOTO, for which we discuss below how AMG has been tuned.
5.1. Summary of results
We are interested in three types of scaling:

Scaling with respect to problem size. The results are summarized in Figure 1, using the exponent of the best fit of run times to the complexity model , where is the total number of arcs. As it can be seen, the scaling of the proposed algorithm is better than the best known MCF algorithms in all cases except for the NETGEN data set. However, as discussed below, there is a moderately large constant complexity factor () with our method compared to combinatorial approaches that is not captured by this estimate, making them more efficient for a small problem with a small number of cores. Ad discussed below, such a constant is heavily dependent on the choice of the MG method and might be improved by using a graph Laplacian AMG.

Scaling with respect to problem complexity. In this case, an exact measure is difficult to be defined, but it is useful to compare the results on the different benchmarks, as done in Figure 2, where a breakdown of run times for the most expensive routines is given for our algorithm. From this comparison, it is clear that NETGEN and GOTO are particularly difficult to be solved. However, such a poor performance can be mostly explained by BoomerAMG, which had to be tuned to avoid producing dense interpolation operators exceeding the memory availability, at the expense of performance. On the other hand, on more regular grids, combinatorial algorithms proved to be much more sensitive to problem complexity, taking 650% more time when solving ROAD_FLOW compared to ROAD_PATHS, while our method showed an increase of only 50%.

Weak scaling tests with respect to number of cores and fixed data set, where the number of arcs per core is kept constant as the problem size is increased. Given that the expected complexity of IPM is , we define the weak scaling efficiency of the th instance as
(30) The results concerning weak scaling are given separately for each data set in the remainder of this section.
5.2. Graph Laplacian AMG
As already discussed, from Figure 1 it is clear that the performance on the NETGEN data set is very poor. In order to confirm that this is a limitation of BoomerAMG, we have tested our approach with the LAMG method [28], which is an algebraic MG solver specifically made for the graph Laplacian, for which a Matlab implementation is available. We have obtained a nearlyoptimal scaling factor of , better than any known combinatorial algorithm, comparable to the BoomerAMG performance on a regular grid network, such as the signal processing data set. This gives a good evidence that an optimalcomplexity IPM solver for the MCF on massively parallel architectures is within reach with our strategy, once an efficient parallel implementation of LAMG becomes available. At the same time, we have also tested the application of LAMG to the IPM without our proposed modification and found that both BoomerAMG and LAMG fail to converge near the solution. This agrees with our discussion in Section 3 and confirms the need to regularize the problem in order to apply efficient solvers.
5.3. Weak scaling tests
We present here the results from weak scaling tests. The hardware used is the Piz Dora supercomputer at the CSCS center. It is a Cray XC40 system with 1256 compute nodes, each equipped with 64 GB of RAM and two 12core Intel Haswell CPUs (Intel Xeon E52690 v3).
5.3.1. Road_paths
As presented in [26], the fastest combinatorial algorithm takes 21.99 seconds on the largest instance (TX_07) of this data set. The results for our method are shown in Figures 3 and 4. As it can be seen, our strategy requires about 200 cores to match the singlecore performance of the best combinatorial algorithm. However, despite the fact that 25872 unknowns per core can be considered too small as a load distribution, the efficiency is above optimal. As can be seen in Figure 1, this is also due to the fact that the method scales with rather than . Therefore, the need for 200 cores is not a symptom of poor scaling but it is a constant complexity factor, so our algorithm would be able to efficiently tackle a larger instance of this data set. From Figure 4, it becomes clear that our modifications to the problem improve the convergence behavior of the linear solver. It is worth remarking that had they been disabled, then rather than decrease, the number of linear iterations would blow up in the last Newton iterations, until the point where BoomerAMG would completely fail to converge.
5.3.2. Road_flow
As presented in [26], the fastest combinatorial algorithm takes 145.7 seconds on the largest instance (TX_07). The results for our method are shown in Figures 5 and 6. By comparing this table with the ROAD_PATHS scenario, we see that as the complexity of the problem increases, run times increase by only with our method, while they increase by more than with combinatorial algorithms. It can also be seen that efficiency is lower in this case, due to a larger increase in the number of Newton iterations required to achieve convergence. We can also notice from Figure 6, that the number of linear iterations has a spike around the Newton iteration 38. We argue that this is an indication that an adaptive selection of the and parameters could give a performance benefit.
5.3.3. Netgen
As presented in [26], the fastest combinatorial algorithm takes 419.4 seconds on the largest instance (22a). The results for the proposed method are shown in Figures 7 and 8. For this data set, we have observed that Newton and MG both converge in a reasonable amount of iterations, however BoomerAMG is struggling in the setup phase, resulting is a very slow assembly of the coarse problems and in a very high memory occupation. We have limited these issues by using 5 levels of aggressive coarsening, at the price of an increased number of linear iterations, as can be seen from Figures 2 and 8. Unfortunately, this is an intrinsic limitation of algebraic MG methods and the implementation of a custom geometric MG would be required to considerably improve efficiency. This will be the subject of future studies.
5.3.4. Goto
As presented in [26], the fastest combinatorial algorithm takes 1564 seconds on the second largest instance (19a), while no data is given for the largest one (22a). It is worth remarking that this is a hard problem and run times for most of the combinatorial algorithms are not given beyond the 16a instance. The results for the proposed method are shown in Figures 9 and 10. For this example, we have observed that the number of Newton iterations grow considerably as the instance size increases. To obtain a reasonable performance, we used the adaptive time step strategy presented in Algorithm 1 with . Despite this, it is clear from the first plot in Figure 10 that the number of linear iterations is far from optimal. This is also due to the fact that, as done in the NETGEN case, we had to limit the memory and computational effort of the BoomerAMG setup phase, by restricting the number of nonzero in the interpolation operator to 5. Therefore, it is possible that also in this case a different MG implementation would improve the convergence behavior. As can be seen from the last plot in Figure 10, this is the only example where the active set becomes already active in the first iterations. In general, this behavior is not ideal, but we found it to be necessary here in order to achieve MG convergence.
5.3.5. Signal processing
The results for the proposed method are shown in Figure 11. If we compare Figures 11 and 7, we immediately see that BoomerAMG is 8x faster on this example. This is not surprising, since the regularity of the mesh makes this example much more similar to FE applications, for which BoomerAMG was specifically developed. This also suggests that there would be a considerable speedup available in the NETGEN case from implementing a custom MG solver for graph Laplacians. However, this example shows that our method is wellsuited to treat problem sizes that are difficult to solve in a reasonable time by a combinatorial algorithm on a single core.
6. Conclusions
We have presented an algorithm to solve the MCF problem on massively parallel architectures. Our approach is based on interpreting the Newton iteration arising from solving the KKT equations as a mixed formulation of the graph Laplacian of the network, with an appropriate metric term that depends on the inequality constraints slack and multipliers. We have proposed an intertwining of IPM and ASM in order to prevent such a metric from vanishing or blowing up, so that it is possible to apply optimalcomplexity linear solvers for elliptic operators, such as AMG. We have shown with standard benchmarks that, while less competitive than combinatorial techniques on small problems that run on a single core, our approach scales well with problem size, complexity, and number of processors, allowing for tackling large problems in a reasonable time.
There are future improvements that can be used to enhance our strategy:

A geometric MG or a parallel graph Laplacian AMG can be used to replace the BoomerAMG solver. This would allow to take advantage of the geometric or graph informations in order to select a more appropriate grid hierarchy for the metric matrix , greatly improving efficiency and memory consumption. Indeed, we have shown that nearly optimalcomplexity is already achievable when the network geometry is regular, such as for the signal processing data, or when the LAMG solver made for graph Laplacians is used. This gives a good indication that an optimalcomplexity parallel solver for the MCF problem is within reach, once a parallel implementation of LAMG is available.

An optimal choice of the parameters , , , and , can be studied for achieving optimal damping ratio of 1, in order to reach the static equilibrium as quickly as possible, and to ensure that the total sum of linear solver iterations over the Newton steps is minimized.
References
 [1] R. K. Ahuja, A. V. Goldberg, J. B. Orlin, and R. E. Tarjan, Finding minimumcost flows by double scaling, Math. Program., Vol. 53, pp. 243–266, 1992.
 [2] S. Balay, S. Abhyankar, M F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, V. Eijkhout, W. D. Gropp, D. Kaushik, M. G. Knepley, L. C. McInnes, K. Rupp, B. F. Smith, S. Zampini, and H. Zhang, PETSc Users Manual, Argonne National Laboratory, ANL95/11  Revision 3.6, 2015.
 [3] D. P. Bertsekas and D. A. Castanon, Parallel primal dual methods for the minimum cost flow problem, Computational Optimization and Applications, Vol. 2, No. 4, pp. 317–336, 1993.
 [4] D. Braess, Finite Elements, Cambridge University Press, New York, 2007.
 [5] W. L. Briggs, V. E. Henson, and S. F. McCormick A Multigrid Tutorial, SIAM Books, Philadelphia, 2000. Second edition.
 [6] W. Hackbusch, MultiGrid Methods and Applications, Springer, 1985
 [7] M. Costantini, A Novel Phase Unwrapping Method Based on Network Programming, IEEE Transactions on Geoscience and Remote Sensing, Vol. 36, No. 3, 1998.
 [8] J. S. Chai and K. C. Toh, Preconditioning and iterative solution of symmetric indefinite linear systems arising from interior point methods for linear programming, Computational Optimization and Applications, Vol. 36, 221–247, 2007.
 [9] S. I. Daitch and D. A. Spielman, Faster Approximate Lossy Generalized Flow via Interior Point Algorithms, in STOC (C. Dwork, ed.), pp. 451–460, ACM, 2008.
 [10] P. Dell’Acqua, A. Frangioni, S. SerraCapizzano, Accelerated multigrid for graph Laplacian operators, Applied Mathematics and Computation, Vol. 270, 193–215, Nov 2015,.
 [11] A. Draganescu and C. Petra, Multigrid preconditioning of linear systems for interior point methods applied to a class of boxconstrained optimal control problems, SIAM Journal on Numerical Analysis, Vol. 50, 328–353, 2012.
 [12] J. Edmonds and R. M. Karp, Theoretical Improvements in Algorithmic Efficiency for Network Flow Problems, J. ACM, Vol. 19, pp. 248–264, Apr. 1972.
 [13] R. D. Falgout, U. M. Yang, hypre: a Library of High Performance Preconditioners, Preconditioners, Lecture Notes in Computer Science, 632–641, 2002.
 [14] L. R. Ford and D. R. Fulkerson, Constructing maximal dynamic flows from static flows., Oper. Res., 6:419–433, 1958.
 [15] M. Friedlander and D. Orban, A primaldual regularized interiorpoint method for convex quadratic programs, Mathematical Programming Computation, Vol. 1, 71–107, 2012.
 [16] A. Frangioni and C. Gentile, Experiments with a hybrid interior point/combinatorial approach for network ow problems, Optimization Methods and Software, Vol. 4, 573–585, 2007.
 [17] H. N. Gabow and R. E. Tarjan, Faster Scaling Algorithms for Network Problems, SIAM J. Comput., Vol. 18, No. 5, pp. 1013–1036, 1989.
 [18] A. Goldberg and R. Tarjan, Solving minimumcost flow problems by successive approximation, in Proc. of the 19th ACM symposium on Theory of computing, STOC’87, (New York, NY, USA), pp. 7–18, ACM press, 1987.
 [19] A. V. Goldberg and R. E. Tarjan., Finding minimumcost circulations by canceling negative cycles, in Proc. of the 20th ACM Symposium on Theory of Computing, STOC?88, (New York, NYm USA), pp. 388–397, ACM Press, 1988.
 [20] A. V. Goldberg and M. Kharitonov, On implementing scaling pushrelabel algorithms for the minimumcost flow problem, Vol. 12. DIMACS Series in Discrete Mathematics and Theoretical Computer Science, 1993.
 [21] J. Gondzio, Matrixfree interior point method. Computational Optimization and Applications 51 (2), 457–480, 2012.
 [22] R. Goldenthal et al, Efficient Simulation of Inextensible Cloth., SIGGRAPH (ACM Transactions on Graphics), 2007.
 [23] D. S. Johnson and C. C. McGeoch, Network Flows and Matching: First DIMACS Implementation Challenge, Boston, MA, USA: American Mathematical Society, 1993.
 [24] N. Karmarkar. A new polynomialtime algorithm for linear programming. Combinatorica, 4(4):373–396, 1984
 [25] M. Klein., A primal method for minimal cost flows with applications to the assignment and transportation problems, Manag. Sci., 14:205–220, 1967.
 [26] P. Kovacs, Minimumcost flow algorithms: an experimental evaluation., Optimization Methods and Software, 30:94127, 2015.
 [27] R. Krause, A Nonsmooth Multiscale Method for Solving Frictional TwoBody Contact Problems in 2D and 3D with Multigrid Efficiency, SIAM J. Sci. Comput., 31(2), 1399–1423, 2008
 [28] O. E. Livne and A. Brandt, Lean Algebraic Multigrid (LAMG): Fast Graph Laplacian Linear Solver, SIAM J. Sci. Comput., 34(4), B499–B522.
 [29] B. Maar and V. Schulz, Interior point multigrid methods for topology optimization, Structural and Multidisciplinary Optimization, Vol, 19, Iss. 3, pp. 214–224, 2000.
 [30] J. Nocedal, S. J. Wright, Numerical Optimization, SpringerVerlag, New York, 1999.
 [31] A. Oliveira and D. Sorensen, A new class of preconditioners for largescale linear systems from interior point methods for linear programming. Linear Algebra and its Applications, Vol. 394, 1–24, 2005.
 [32] J. B. Orlin, A Faster Strongly Polynominal Minimum Cost Flow Algorithm, in STOC (J. Simon,ed.), pp. 377–387, ACM, 1988.
 [33] L. F. Portugal, M. G. C. Resende, G. Veiga, J. Patrício, and J. J. Júdice. PDNET–Software for network flow optimization using an interior point method. 2004.
 [34] L. F. Portugal, M. G. C. Resende, G. Veiga, J. Patrício, and J. J. Júdice. Fortran subroutines for network flow optimization using an interior point algorithm. Pesquisa Operacional, 28:243?261, 2008.
 [35] N. Tomizawa, On some techniques useful for solution of transportation network problems, Networks, 1:173–194, 1971.