On the solution of large-scale SDP problems by the modified barrier method using iterative solvers

The limiting factors of second-order methods for large-scale semidefinite optimization are the storage and factorization of the Newton matrix. For a particular algorithm based on the modified barrier method, we propose to use iterative solvers instead of the routinely used direct factorization techniques. The preconditioned conjugate gradient method proves to be a viable alternative for problems with a large number of variables and modest size of the constrained matrix. We further propose to avoid explicit calculation of the Newton matrix either by an implicit scheme in the matrix–vector product or using a finite-difference formula. This leads to huge savings in memory requirements and, for certain problems, to further speed-up of the algorithm. Dedicated to the memory of Jos Sturm.     

Exact Max-SAT solvers for over-constrained problems

We present a new generic problem solving approach for over-constrained problems based on Max-SAT. We first define a Boolean clausal form formalism, called soft CNF formulas, that deals with blocks of clauses instead of individual clauses, and that allows one to declare each block either as hard (i.e., must be satisfied by any solution) or soft (i.e., can be violated by some solution). We then present two Max-SAT solvers that find a truth assignment that satisfies all the hard blocks of clauses and the maximum number of soft blocks of clauses. Our solvers are branch and bound algorithms equipped with original lazy data structures, powerful inference techniques, good quality lower bounds, and original variable selection heuristics. Finally, we report an experimental investigation on a representative sample of instances (random 2-SAT, Max-CSP, graph coloring, pigeon hole and quasigroup completion) which provides experimental evidence that our approach is very competitive compared with the state-of-the-art approaches developed in the CSP and SAT communities. Research partially supported by projects TIN2004-07933-C03-03 and TIC2003-00950 funded by the Ministerio de Educación y Ciencia. The second author is supported by a grant Ramón y Cajal.     

A formal approach to subgrammar extraction for NLP

The problem of subgrammar extraction is precisely defined in the following way: Given a grammar $G$ and a set of words $W$, find a smallest subgrammar of $G$ that accepts the same set of sentences from $W^1$ as $G$, and for each of them produces the same parse trees. In practical Natural Language Processing applications, the set of words $W$ is obtained from the text unit. There are practical motivations for doing this operation “just-in-time”, i.e. just before processing the text; hence it is required that this operation be done in an automatic and efficient way. After defining the problem in the general framework, we discuss the problem for context-free grammars (CFG), and give an efficient algorithm for it. We prove that finding the smallest subgrammar for HPSGs is an NP-hard problem, and give an efficient algorithm that solves an easier, approximate version of the problem. We also discuss how the algorithm can be efficiently implemented.     

Parallel iterative solvers for boundary value methods

A parallel variant of the block Gauss-Seidel iteration for the solution of block-banded linear systems is presented. The coefficient matrix is partitioned among the processors as in the domain decomposition methods and then it is split so that the resulting iterative method has the same spectral properties of the block Gauss-Seidel iteration.The parallel algorithm is applied to the solution of block-banded linear systems arising from the numerical discretization of initial value problems by means of Boundary Value Methods (BVMs). BVMs define a new approach for the solution of ordinary differential equations and seem to be attractive for their interesting stability properties and a possible parallel implementation. In this paper, we refer to BVMs based on the extended trapezoidal rules.     

Structured Backward Errors for Structured KKT Systems

In this paper we study structured backward errors for some structured KKT systems. Normwise structured backward errors for structured KKT systems are defined, and computable formulae of the structured backward errors are obtained. Simple numerical examples show that the structured backward errors may be much larger than the unstructured ones in some cases.     

Smooth over-parameterized solvers for non-smooth structured optimization

Mathematical Programming - Non-smooth optimization is a core ingredient of many imaging or machine learning pipelines. Non-smoothness encodes structural constraints on the solutions, such as...     

Multi-stage solvers optimized for damping and propagation

Explicit multi-stage solvers are routinely used to solve the semi-discretized equations that arise in Computational Fluid Dynamics (CFD) problems. Often they are used in combination with multi-grid methods. In that case, the role of the multi-stage solver is to efficiently reduce the high frequency modes on the current grid and is called a smoother. In the past, when optimizing the coefficients of the scheme, only the damping characteristics of the smoother were taken into account and the interaction with the remainder of the multi-grid cycle was neglected. Recently it had been found that coefficients that result in less damping, but allow for a higher Courant-Friedrichs-Lewy (CFL) number are often superior to schemes that try to optimize damping alone. While this is certainly true for multi-stage schemes used as a stand-alone solver, we investigate in this paper if using higher CFL numbers also yields better results in a multi-grid setting. We compare the results with a previous study we conducted and where a more accurate model of the multi-grid cycle was used to optimize the various parameters of the solver.We show that the use of the more accurate model results in better coefficients and that in a multi-grid setting propagation is of little importance.We also look into the gains to be made when we allow the parameters to be different for the pre- and post-smoother and show that even better coefficients can be found in this way.     

Parallel linear system solvers for Runge-Kutta methods

If the nonlinear systems arising in implicit Runge-Kutta methods like the Radau IIA methods are iterated by (modified) Newton, then we have to solve linear systems whose matrix of coefficients is of the form I-A hJ with A the Runge-Kutta matrix and J an approximation to the Jacobian of the righthand side function of the system of differential equations. For larger systems of differential equations, the solution of these linear systems by a direct linear solver is very costly, mainly because of the LU-decompositions. We try to reduce these costs by solving the linear systems by a second (inner) iteration process. This inner iteration process is such that each inner iteration again requires the solution of a linear system. However, the matrix of coefficients in these new linear systems is of the form I - B hJ where B is similar to a diagonal matrix with positive diagonal entries. Hence, after performing a similarity transformation, the linear systems are decoupled into s subsystems, so that the costs of the LU-decomposition are reduced to the costs of s LU-decompositions of dimension d. Since these LU-decompositions can be computed in parallel, the effective LU-costs on a parallel computer system are reduced by a factor s ³ . It will be shown that matrices B can be constructed such that the inner iterations converge whenever A and J have their eigenvalues in the positive and nonpositive halfplane, respectively. The theoretical results will be illustrated by a few numerical examples. A parallel implementation on the four-processor Cray-C98/4256 shows a speed-up ranging from at least 2.4 until at least 3.1 with respect to RADAU5 applied in one-processor mode.     

Direct and iterative solvers for finite-element problems

Parallelization promises to allow increasing performance demand from users of scientific computation codes. In this article we discuss implementation of parallelization techniques in industrial computation codes for applications in Fluid Mechanics (N3S) and Structural Mechanics (Code_ASTER). We outline the gains afforded by parallelization on a Cray C98 machine in production mode, using a single version of the program running in sequential or in parallel mode, with minimum modifications to the initial code. Parallelization is applied at the discretion of the user, for whom no specialist knowledge of this technique is assumed. This revised version was published online in August 2006 with corrections to the Cover Date.     

Hybrid solvers for composition and splitting methods

Composition and splitting are useful techniques for constructing special purpose integration methods for numerically solving many types of differential equations. In this article we will review these methods and summarise the essential ingredients of an implementation that has recently been added to a framework for solving differential equations in Mathematica.     

High-order compact solvers for the three-dimensional Poisson equation

New compact approximation schemes for the Laplace operator of fourth- and sixth-order are proposed. The schemes are based on a Padé approximation of the Taylor expansion for the discretized Laplace operator. The new schemes are compared with other finite difference approximations in several benchmark problems. It is found that the new schemes exhibit a very good performance and are highly accurate. Especially on large grids they outperform noncompact schemes.     

Time-varying Riemann solvers for conservation laws on networks

We consider a conservation law on a network and generic Riemann solvers at nodes depending on parameters, which can be seen as control functions. Assuming that the parameters have bounded variation as functions of time, we prove existence of solutions to Cauchy problems on the whole network.     

A new row ordering strategy for frontal solvers

The frontal method is a variant of Gaussian elimination that has been widely used since the mid 1970s. In the innermost loop of the computation the method exploits dense linear algebra kernels, which are straightforward to vectorize and parallelize. This makes the method attractive for modern computer architectures. However, unless the matrix can be ordered so that the front is never very large, frontal methods can require many more floating‐point operations for factorization than other approaches. We are interested in matrices that have a highly asymmetric structure. We use the idea of a row graph of an unsymmetric matrix combined with a variant of Sloan's profile reduction algorithm to reorder the rows. We also look at applying the spectral method to the row graph. Numerical experiments performed on a range of practical problems illustrate that our proposed MSRO and hybrid MSRO row ordering algorithms yield substantial reductions in the front sizes and, when used with a frontal solver, significantly enhance its performance both in terms of the factorization time and storage requirements. Copyright © 1999 John Wiley & Sons, Ltd.     

Inverse problems for regularization matrices

Discrete ill-posed problems are difficult to solve, because their solution is very sensitive to errors in the data and to round-off errors introduced during the solution process. Tikhonov regularization replaces the given discrete ill-posed problem by a nearby penalized least-squares problem whose solution is less sensitive to perturbations. The penalization term is defined by a regularization matrix, whose choice may affect the quality of the computed solution significantly. We describe several inverse matrix problems whose solution yields regularization matrices adapted to the desired solution. Numerical examples illustrate the performance of the regularization matrices determined.     

Multiscale deflation solvers for flow in porous media

This work describes a novel physics-based deflation preconditioner approach for solving porous media flow problems characterized by highly heterogeneous media. The approach relies on high-conductivity block solutions after rearranging the linear system coefficients into high-conductive and low condutive blocks from a given physically driven threshold value. This rearranging relies on the Hoshen-Kopelman (H-K) algorithm that is commonly used to determine percolation clusters. The resulting preconditioner may alternatively be combined with a deflation preconditioning stage. The proposed approach is coined as a physics-based 2-stage deflation preconditioner (P2SDP). Numerical experiments on different permeability distributions reveal that P2SDP is a powerful means to solve pressure systems when compared to more conventional algebraic approaches such as the incomplete Cholesky factorization. (© 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim)     

Parallel multilevel solvers for the cardiac electro-mechanical coupling

We develop a parallel solver for the cardiac electro-mechanical coupling. The electric model consists of two non-linear parabolic partial differential equations (PDEs), the so-called Bidomain model, which describes the spread of the electric impulse in the heart muscle. The two PDEs are coupled with a non-linear elastic model, where the myocardium is considered as a nearly-incompressible transversely isotropic hyperelastic material. The discretization of the whole electro-mechanical model is performed by Q1 finite elements in space and a semi-implicit finite difference scheme in time. This approximation strategy yields at each time step the solution of a large scale ill-conditioned linear system deriving from the discretization of the Bidomain model and a non-linear system deriving from the discretization of the finite elasticity model. The parallel solver developed consists of solving the linear system with the Conjugate Gradient method, preconditioned by a Multilevel Schwarz preconditioner, and the non-linear system with a Newton–Krylov-Algebraic Multigrid solver. Three-dimensional parallel numerical tests on a Linux cluster show that the parallel solver proposed is scalable and robust with respect to the domain deformations induced by the cardiac contraction.     

Parallel iterative linear solvers for multistep Runge-Kutta methods

This paper deals with solving stiff systems of differential equations by implicit Multistep Runge-Kutta (MRK) methods. For this type of methods, nonlinear systems of dimension sd arise, where s is the number of Runge-Kutta stages and d the dimension of the problem. Applying a Newton process leads to linear systems of the same dimension, which can be very expensive to solve in practice. With a parallel iterative linear system solver, especially designed for MRK methods, we approximate these linear systems by s systems of dimension d, which can be solved in parallel on a computer with s processors. In terms of Jacobian evaluations and LU-decompositions, the k-steps-stage MRK applied with this technique is on s processors equally expensive as the widely used k-step Backward Differentiation Formula on 1 processor, whereas the stability properties are better than that of BDF. A simple implementation of both methods shows that, for the same number of Newton iterations, the accuracy delivered by the new method is higher than that of BDF.     

A sequential cutting plane algorithm for solving convex NLP problems

In this paper we look at a new algorithm for solving convex nonlinear programming optimization problems. The algorithm is a cutting plane-based method, where the sizes of the subproblems remain fixed, thus avoiding the issue with constantly growing subproblems we have for the classical Kelley’s cutting plane algorithm. Initial numerical experiments indicate that the algorithm is considerably faster than Kelley’s cutting plane algorithm and also competitive with existing nonlinear programming algorithms.     

Improve-and-Branch Algorithm for the Global Optimization of Nonconvex NLP Problems

A new algorithm to solve nonconvex NLP problems is presented. It is based on the solution of two problems. The reformulated problem RP is a suitable reformulation of the original problem and involves convex terms and concave univariate terms. The main problem MP is a nonconvex NLP that outer-approximates the feasible region and underestimate the objective function. MP involves convex terms and terms which are the products of concave univariate functions and new variables. Fixing the variables in the concave terms, a convex NLP that overestimates the feasible region and underestimates the objective function is obtained from the MP. Like most of the deterministic global optimization algorithms, bounds on all the variables in the nonconvex terms must be provided. MP forces the objective value to improve and minimizes the difference of upper and lower bound of all the variables either to zero or to a positive value. In the first case, a feasible solution of the original problem is reached and the objective function is improved. In general terms, the second case corresponds to an infeasible solution of the original problem due to the existence of gaps in some variables. A branching procedure is performed in order to either prove that there is no better solution or reduce the domain, eliminating the local solution of MP that was found. The MP solution indicates a key point to do the branching. A bound reduction technique is implemented to accelerate the convergence speed. Computational results demonstrate that the algorithm compares very favorably to other approaches when applied to test problems and process design problems. It is typically faster and it produces very accurate results.