scispace - formally typeset
Search or ask a question
Journal ArticleDOI

Matching-based preprocessing algorithms to the solution of saddle-point problems in large-scale nonconvex interior-point optimization

01 Apr 2007-Computational Optimization and Applications (Kluwer Academic Publishers-Plenum Publishers)-Vol. 36, Iss: 2, pp 321-341
TL;DR: This work presents combinatorial methods to preprocess these matrices to establish more favorable numerical properties for the subsequent factorization in a sparse direct LDLT factorization method where the pivoting is restricted to static supernode data structures.
Abstract: Interior-point methods are among the most efficient approaches for solving large-scale nonlinear programming problems. At the core of these methods, highly ill-conditioned symmetric saddle-point problems have to be solved. We present combinatorial methods to preprocess these matrices in order to establish more favorable numerical properties for the subsequent factorization. Our approach is based on symmetric weighted matchings and is used in a sparse direct LDL T factorization method where the pivoting is restricted to static supernode data structures. In addition, we will dynamically expand the supernode data structure in cases where additional fill-in helps to select better numerical pivot elements. This technique can be seen as an alternative to the more traditional threshold pivoting techniques. We demonstrate the competitiveness of this approach within an interior-point method on a large set of test problems from the CUTE and COPS sets, as well as large optimal control problems based on partial differential equations. The largest nonlinear optimization problem solved has more than 12 million variables and 6 million constraints.

Summary (2 min read)

1 Introduction

  • In recent years, a large amount of work has been devoted to the problem of solving large symmetric indefinite systems in saddle-point form efficiently.
  • For simplicity in the notation the authors assume without loss of generality that all variables have only a lower bound.
  • For a detailed survey on solution techniques for large linear saddle-point systems the interested reader should consult [2].
  • Unfortunately, until recently, sparse symmetric indefinite factorization methods were relatively inefficient compared to symmetric positive definite solvers [22].
  • The matching algorithms used to reorder the matrix before the factorization are discussed in Sect.

2 Optimization algorithm

  • Over the past ten years, a number of primal–dual interior-point methods for general nonlinear continuous optimization problems, such as (1.1), have been proposed (see [19] for a comprehensive survey).
  • Therefore, at least close to the solution, the algorithms derived from both points of views are very similar, and fast local convergence can be achieved by either.
  • To ensure this, SQP-type line search methods usually require that the projection of W̃k (or a modification thereof) onto the null space of ATk is positive definite, or, equivalently, that the matrix in (2.4) is non-singular and has exactly m (the number of constraints) negative eigenvalues [35].
  • The perturbation parameter for the Hessian, δx ≥ 0, is chosen to ensure that the matrix in (2.6) has the correct inertia, and that therefore the search direction has the required descent properties.

3 Supernodal factorizations methods

  • The authors use a Level-3 BLAS supernodal left-looking factorization as described in [27, 33, 38] and discuss two additional pivoting techniques for symmetric indefinite systems that are described below.
  • The first pivoting technique has also been used in [38]—the second variant is novel and is used to improve the accuracy.
  • The element growth in the diagonal block after k steps is bounded by the factor (2.57)k−1.
  • In the extreme case, the supernode exists of only one column and the SBK pivoting can degenerate to diagonal pivoting.
  • In Sect. 4, the authors will discuss these permutations which are based on combinatorial graph algorithms.

4 Symmetric matchings for a-priori pivoting

  • Matching algorithms work on the associate graph representations of the matrices.
  • The use of combinatorial techniques in the solution of indefinite linear systems goes back to the early eighties.
  • Duff and Pralet [13] showed that the problem of finding symmetric maximum weight matchings in a symmetric matrix is equivalent to finding a non-bipartite matching.
  • The matching algorithm also provides, as a byproduct, a row scaling Dr and a column scaling Dc.

5 Numerical experiments

  • In this section the authors present a number of numerical results to explore the robustness and efficiency of different direct linear solvers within the interior point code IPOPT.
  • The solvers used for the comparisons is the PARDISO solver with the different ordering strategies, as well as the Harwell routines MA27 and MA57 (Version 3.0) [10].
  • Further to the right the authors can see how far behind an option is (e.g., in about 55% of the problems, PARDISO1 was not more than 21 = 2 times worse than the best solver for each individual instance).
  • 2 Optimal control problems with discretized PDEs.

5.2.1 Parallel solution

  • In their final experiments the authors explore the parallel scale-up of the PARDISO code (only using the PARDISO2 method), solving the discretized PDE optimal control problems for instances with several million variables.
  • From each problem class the authors picked the first instance, and used the same options for IPOPT as before.
  • The table shows the number of grip points N for each problem, the fraction of runtime for the linear solver PARDISO compared to IPOPT, and the numbers of processors available to the direct solver.

6 Conclusion

  • The authors presented an analysis of direct solution methods based on restricted pivoting and combinatorial preprocessing for the solution of large-scale saddle-point problems stemming from interior-point optimization methods.
  • The results indicate that the combination of the preprocessing step based on weighted symmetric matchings, and the static factorization with Supernodal-Bunch– Kaufman pivoting provides sufficient accuracy both in terms of inertia information (an important requirement for solving nonconvex nonlinear optimization problems), as well as accuracy of the solution.
  • Let us briefly recall the main ingredients necessary for this progress: Furthermore, the approximate constraint matching approach is also of importance for the computation of orderings and scalings at a manageable cost.
  • The authors thank Andrew Conn for valuable suggestions and Iain Duff for a temporary license from HSL.

Did you find this useful? Give us your feedback

Content maybe subject to copyright    Report

Comput Optim Appl (2007) 36: 321–341
DOI 10.1007/s10589-006-9003-y
Matching-based preprocessing algorithms to the
solution of saddle-point problems in large-scale
nonconvex interior-point optimization
Olaf Schenk ·Andreas Wächter ·
Michael Hagemann
Published online: 22 February 2007
© Springer Science+Business Media, LLC 2007
Abstract Interior-point methods are among the most efficient approaches for solving
large-scale nonlinear programming problems. At the core of these methods, highly
ill-conditioned symmetric saddle-point problems have to be solved. We present com-
binatorial methods to preprocess these matrices in order to establish more favorable
numerical properties for the subsequent factorization. Our approach is based on sym-
metric weighted matchings and is used in a sparse direct LDL
T
factorization method
where the pivoting is restricted to static supernode data structures. In addition, we
will dynamically expand the supernode data structure in cases where additional fill-
in helps to select better numerical pivot elements. This technique can be seen as an
alternative to the more traditional threshold pivoting techniques. We demonstrate the
competitiveness of this approach within an interior-point method on a large set of test
problems from the CUTE and COPS sets, as well as large optimal control problems
based on partial differential equations. The largest nonlinear optimization problem
solved has more than 12 million variables and 6 million constraints.
Keywords Nonconvex nonlinear programming · Interior-point method ·
Saddle-point problem · Numerical linear algebra · Maximum weight matching
O. Schenk (
) · M. Hagemann
Departement of Computer Science, University of Basel, Klingelbergstr. 50, 4056 Basel,
Switzerland
e-mail: olaf.schenk@unibas.ch
M. Hagemann
e-mail: michael.hagemann@unibas.ch
A. Wächter
IBM T.J. Watson Research Center, Yorktown Heights, NY, USA
e-mail: andreasw@watson.ibm.com

322 O. Schenk et al.
1 Introduction
In recent years, a large amount of work has been devoted to the problem of solving
large symmetric indefinite systems in saddle-point form efficiently. One reason for
this surge in interest is due the success of interior-point methods in nonlinear pro-
gramming, which at their core require the solution of a series of linear systems in
saddle-point form. Consider a nonlinear programming problem given by an objective
function f(x): R
n
R and constraint functions c(x) :R
n
R
m
, which are both as-
sumed to be twice continuously differentiable. The objective is to find a local solution
of the optimization problem
min
xR
n
f(x) (1.1a)
s.t. c(x) =0, (1.1b)
x 0. (1.1c)
For simplicity in the notation we assume without loss of generality that all variables
have only a lower bound. Also note that problems with general inequality constraints
can be transformed into the above formulation by introducing slack variables.
In an interior-point optimization framework, the solution of (1.1) is found by a se-
ries of Newton-type iterations that require in each step the solution of a linear system
of equations of the form
Kx =
HA
A
T
C
x =b (1.2)
where the n × n matrix H is symmetric and potentially indefinite, C is a diagonal
regularization matrix which is typically small and often zero, and the n × m matrix
A has full column rank. See Sect. 2 for a detailed description of these matrices. For
a detailed survey on solution techniques for large linear saddle-point systems the
interested reader should consult [2].
Commonly, sparse direct factorization methods for symmetric indefinite systems
areusedtosolve(1.2). These methods compute, for example, sparse Bunch–Kaufman
[6]orDuffReid[14] factorizations
K =PQLDL
T
Q
T
P
T
(1.3)
where Q is a pivoting permutation matrix used for numerical stability, P is a fill-in
minimization permutation matrix, L is unit lower triangular, and D is a block diag-
onal matrix with blocks of dimension 1 and 2. Unfortunately, until recently, sparse
symmetric indefinite factorization methods were relatively inefficient compared to
symmetric positive definite solvers [22]. However, with the invention of fast combi-
natorial algorithms [11, 36], which improve the diagonal dominance of the linear sys-
tems, the situation has dramatically changed. The use of matching-based preprocess-
ing steps for (1.2) leads to symmetric indefinite sparse direct solvers that are almost
as efficient as their positive definite counterparts. The key idea is to transform the
matrix K in (1.2)into
ˆ
K =P
M
D
s
KD
T
s
P
T
M
, (1.4)

Matching-based preprocessing algorithms to the solution of saddle-point problems 323
such that the permuted system has a greater block-diagonal dominance than the origi-
nal matrix, and to use factorization methods that honor the block-diagonal dominance
of (1.4). Here, P
M
is a permutation matrix that permutes large-off diagonal entries
to be close to the diagonal, and the diagonal matrix D
s
is a scaling matrix such that
both columns and rows of
ˆ
K have unit infinity-norm.
The paper is organized as follows. Section 2 provides background on interior-
point methods and discusses the requirements for factorization algorithms to solve
(1.2). We briefly review the supernodal factorization approach for the direct solution
of saddle-point matrices in Sect. 3. The matching algorithms used to reorder the ma-
trix before the factorization are discussed in Sect. 4. Section 5 introduces the test
problems and gives various performance statistics and comparisons of the proposed
methods.
2 Optimization algorithm
Over the past ten years, a number of primal–dual interior-point methods for general
nonlinear continuous optimization problems, such as (1.1), have been proposed (see
[19] for a comprehensive survey). They can be derived in roughly two different ways,
which we outline in the following. We contrast the two approaches in order to em-
phasize the importance of a proper handling of the nonconvex case.
The more recent point of view was developed in the 1980’s, initially for problems
in which f(x)and c(x) are linear. Here, the perturbed primal–dual equations
f(x)+∇c(x)λ z = 0, (2.1a)
c(x) = 0, (2.1b)
XZe μe = 0 (2.1c)
play a central role. The vectors λ and z are the Lagrangian multipliers for the equality
and bound constraints, and X and Z denote the diagonal matrices with the vector
elements of x and z on the diagonal. The vector e denotes the vector of all ones,
i.e., e = (1,...,1)
T
, of appropriate dimension. For μ = 0, Eqs. (2.1) together with
x,z 0 are the first-order optimality conditions for (1.1). Methods in this class
obtain a point satisfying those conditions by applying Newton’s method to (2.1),
where μ takes on strictly positive values, and, in a homotopy approach, is eventually
driven to zero. The steps are obtained from
W
k
A
k
I
A
T
k
00
Z
k
0 X
k

x
k
λ
k
z
k
=−
f(x
k
) +A
k
λ
k
z
k
c(x
k
)
X
k
Z
k
e μ
k
e
. (2.2)
Here A
k
=∇c(x
k
), and W
k
denotes the Hessian of the Lagrangian function for (1.1)
w.r.t. x. Once a step has been computed, an appropriate step size is chosen that among
other things ensures that non-negativity conditions (x,z > 0) are satisfied for each
iterate. Using the fact that the Jacobian of the system of Eq. (2.1) is non-singular at
a non-degenerate local solution (x
,z
) of (1.1), it is possible to derive methods
that converge superlinearly toward (x
,z
); see, e.g., [23].

324 O. Schenk et al.
A number of methods have been proposed that apply the primal–dual strategy as
outlined above directly to nonlinear optimization problems (e.g., [15, 40]). Progress
towards the solution is ensured by monitoring the norm of the optimality conditions
((2.1) with μ =0) and driving it to zero. While this is appropriate for problems where
every solution of (2.1) is a solution for (1.1) (for example, when (1.1) corresponds to
a convex optimization problem), we believe it is not suitable in the general nonconvex
case. In a sense, these methods ignore the optimization aspect of the problem, and re-
duce the algorithm to finding a root of the optimality conditions. However, in practice,
this can easily lead to convergence to stationary points that are not minimizers.
A different point of view of interior-point methods, originating from research di-
rectly for general nonlinear optimization problems, has been discussed and analyzed
in detail by Fiacco and McCormick [16] in the 1960s. Given the original problem in
the form (1.1), a sequence of corresponding barrier problems,
min
xR
n
ϕ
μ
(x) =f(x)μ
n
i=1
ln
x
(i)
(2.3a)
s.t. c(x) =0, (2.3b)
is solved to increasingly tighter tolerances, while again the barrier parameter μ is
driven to zero. Here, it is possible to use techniques (for the solution of the individual
barrier problems) that have been developed for general (nonconvex) optimization,
such as SQP-type methods. In particular, we may use globalization approaches that
have been proposed to handle the case where not all stationary points satisfying the
first order optimality conditions are local minimizers. In the following, we restrict
our discussion to line search algorithms, such as the one implemented in I
POPT [41],
the code used for the numerical results presented in Sect. 5.
The search directions are obtained from solving the linear system
˜
W
k
A
k
A
T
k
0

x
k
λ
k
=−
ϕ
μ
(x
k
) +A
k
λ
k
c(x
k
)
, (2.4)
where
˜
W
k
is an approximation of the Hessian of the Lagrangian for the barrier prob-
lem (2.3). One can show that these steps, together with
z
k
=μX
1
k
e z
k
Σ
k
x
k
, (2.5)
correspond to solutions of (2.2), if
˜
W
k
=W
k
+Σ
k
with Σ
k
=X
1
k
Z
k
. Therefore, at
least close to the solution, the algorithms derived from both points of views are very
similar, and fast local convergence can be achieved by either.
However, the crucial difference between the two classes of algorithms lies in how
they behave when the iterates are not close to a solution. SQP-type method use glob-
alization techniques based on a merit function or a filter method (the approach used
in I
POPT,see[42] for details). Here, it is crutial that the search directions gener-
ated from (2.4) have appropriate descent properties to guarantee that progress can
be made for the globalization framework, so that overall convergence to non-optimal
stationary points is less likely or can be avoided. To ensure this, SQP-type line search
methods usually require that the projection of
˜
W
k
(or a modification thereof) onto

Matching-based preprocessing algorithms to the solution of saddle-point problems 325
the null space of A
T
k
is positive definite, or, equivalently, that the matrix in (2.4)is
non-singular and has exactly m (the number of constraints) negative eigenvalues [35].
In contrast to this, a method that uses the unmodified Hessian matrix in (2.2) might
generate steps that are increasing the objective function even for arbitrarily small
step sizes, if the current iterate is feasible and the projection of
˜
W
k
is indefinite. We
therefore strongly believe that it is desirable to include the fact that one is dealing
with a minimization problem in the design of an optimization method for nonconvex
optimization.
2.1 Requirements for direct solvers
A number of algorithms for the factorization of indefinite symmetric systems provide
information about the inertia of the factorized matrix on the fly, e.g., [6]. In the context
of an SQP-type line search optimization method it is therefore possible to ensure the
required descent properties by performing a “trial” factorization for systems of the
form (2.4), where the matrix
˜
W
k
is modified according to some heuristics, until the
inertia of the matrix is correct.
The method implemented in I
POPT is a primal–dual interior-point algorithm that
generates iterates for the variables, (x
k
k
,z
k
), using steps from the linear system
(2.2), but for the reasons just given, instead of solving this nonsymmetric system
directly, it solves a “regularized” form of (2.4), namely,
W
k
+Σ
k
+δ
x
IA
k
A
T
k
δ
c
I

x
k
λ
k
=−
ϕ
μ
(x
k
) +A
k
λ
k
¸(x
k
)
, (2.6)
and obtains the search direction for the bound multipliers from (2.5). The perturbation
parameter for the Hessian, δ
x
0, is chosen to ensure that the matrix in (2.6) has
the correct inertia, and that therefore the search direction has the required descent
properties. In addition, a small perturbation parameter δ
c
10
8
might be chosen
in case the matrix A
k
appears to be column rank-deficient, to ensure that the matrix
in (2.6) is non-singular. The heuristic for choosing δ
x
and δ
c
implemented in IPOPT
is described in [41]; it usually attempts first to work with δ
x
= δ
c
= 0 to obtain the
pure primal–dual Newton direction for fast convergence, and only if the inertia is not
correct with this choice are other values tried.
Another approach for modifying the Hessian matrix is to use an inertia controlling
factorization method (see, e.g., [17]), which modifies the matrix on the fly during the
factorization and ensures the desired descent properties.
Using an iterative solver for the computation of a good search direction from (2.2)
or (2.6) appears difficult, since it seems not possible to ensure the required descent
properties. A hybrid approach has been used in trust-region algorithms (see, e.g., [7]),
where a matrix of the form (2.4), with
˜
W
k
replaced by Σ
k
, is factorized and used to
compute the feasibility component of the search direction, and an iterative method,
working in the null space of the constraint Jacobian, is used to compute the remaining
component [7]. If the inertia of the linear system is not correct, this iterative procedure
will encounter directions of negative curvature, which now can be exploited explic-
itly, so that a good descent search direction is obtained. Other iterative procedures for
the solution of saddle point problems have been proposed, but they usually require

Citations
More filters
01 Jan 2006
TL;DR: Three new variations of a direct factorization scheme to tackle the is- sue of indeniteness in sparse symmetric linear systems, including a reordering that is based on a symmetric weighted matching of the matrix, which is effective for highly indenite symmetric systems.
Abstract: This paper discusses new pivoting factorization methods for solving sparse symmetric indenite sys- tems. As opposed to many existing pivoting methods, our SupernodenBunchnKaufman (SBK) pivoting method dy- namically selects and pivots and may be supplemented by pivot perturbation techniques. We demonstrate the effectiveness and the numerical accuracy of this algorithm and also show that a high performance implementa- tion is feasible. We will also show that symmetric maximum-weighted matching strategies add an additional level of reliability to SBK. These techniques can be seen as a complement to the alternative idea of using more complete pivoting techniques during the numerical factorization. Numerical experiments validate these conclusions. where is a diagonal matrix with and pivot blocks, is a sparse lower triangu- lar matrix, and is a symmetric indenite diagonal matrix that reects small half-machine precision perturbations, which might be necessary to tackle the problem of tiny pivots. is a reordering that is based on a symmetric weighted matching of the matrix , and tries to move the largest off-diagonal elements directly alongside the diagonal in order to form good initial or diagonal block pivots. is a ll reducing reordering which honors the structure of . We will present three new variations of a direct factorization scheme to tackle the is- sue of indeniteness in sparse symmetric linear systems. These methods restrict the pivoting search, to stay as long as possible within predened data structures for efcient Level-3 BLAS factorization and parallelization. On the other hand, the imposed pivoting restrictions can be reduced in several steps by taking the matching permutation into account. The rst al- gorithm uses SupernodenBunchnKaufman (SBK) pivoting and dynamically selects and pivots. It is supplemented by pivot perturbation techniques. It uses no more storage than a sparse Cholesky factorization of a positive denite matrix with the same sparsity structure due to restricting the pivoting to interchanges within the diagonal block associated to a single supernode. The coefcient matrix is perturbed whenever numerically acceptable and pivots cannot be found within the diagonal block. One or two steps of iterative re- nement may be required to correct the effect of the perturbations. We will demonstrate that this restricting notion of pivoting with iterative renement is effective for highly indenite symmetric systems. Furthermore the accuracy of this method is for a large set of matrices from different applications areas as accurate as a direct factorization method that uses com- plete sparse pivoting techniques. In addition, we will discuss two preprocessing algorithms to identify large entries in the coefcient matrix that, if permuted close to the diagonal, permit

474 citations

Journal ArticleDOI
TL;DR: This article gives an algorithm that computes a (1 − 1 − 0))-approximate maximum weight matching in O(i) time, that is, optimal linear time for any fixed ε, and should be appealing in all applications that can tolerate a negligible relative error.
Abstract: The maximum cardinality and maximum weight matching problems can be solved in O(m√n) time, a bound that has resisted improvement despite decades of research. (Here m and n are the number of edges and vertices.) In this article, we demonstrate that this “m√n barrier” can be bypassed by approximation. For any e > 0, we give an algorithm that computes a (1 − e)-approximate maximum weight matching in O(me−1 log e−1) time, that is, optimal linear time for any fixed e. Our algorithm is dramatically simpler than the best exact maximum weight matching algorithms on general graphs and should be appealing in all applications that can tolerate a negligible relative error.

287 citations


Cites methods from "Matching-based preprocessing algori..."

  • ...) Approximate MWM algorithms are used as a heuristic preprocessing step in several sparse linear system solvers [Duff and Gilbert 2002; Hagemann and Schenk 2006; Olschowka and Neumaier 1996; Schenk et al. 2007]....

    [...]

  • ...Approximate MWM algorithms are used as a heuristic preprocessing step in several sparse linear system solvers [Duff and Gilbert 2002; Hagemann and Schenk 2006; Olschowka and Neumaier 1996; Schenk et al. 2007]....

    [...]

Journal ArticleDOI
TL;DR: Efficient preconditioning algorithms for an eigenvalue problem arising in quantum physics, namely, the computation of a few interior eigenvalues and their associated eigenvectors for large-scale sparse real and symmetric indefinite matrices of the Anderson model of localization, are proposed.
Abstract: We propose efficient preconditioning algorithms for an eigenvalue problem arising in quantum physics, namely, the computation of a few interior eigenvalues and their associated eigenvectors for large-scale sparse real and symmetric indefinite matrices of the Anderson model of localization. We compare the Lanczos algorithm in the 1987 implementation by Cullum and Willoughby with the shift-and-invert techniques in the implicitly restarted Lanczos method and in the Jacobi-Davidson method. Our preconditioning approaches for the shift-and-invert symmetric indefinite linear system are based on maximum weighted matchings and algebraic multilevel incomplete $LDL^T$ factorizations. These techniques can be seen as a complement to the alternative idea of using more complete pivoting techniques for the highly ill-conditioned symmetric indefinite Anderson matrices. We demonstrate the effectiveness and the numerical accuracy of these algorithms. Our numerical examples reveal that recent algebraic multilevel preconditioning solvers can accelerate the computation of a large-scale eigenvalue problem corresponding to the Anderson model of localization by several orders of magnitude.

166 citations


Cites background from "Matching-based preprocessing algori..."

  • ...Furthermore, these techniques have been successfully transferred to the symmetric case [24, 26], allowing modern state-of-the-art direct solvers [59, 60, 62] to be orders of magnitudes faster and more memory efficient than ever, finally leading to symmetric indefinite sparse direct solvers that are almost as efficient as their symmetric positive definite counterparts....

    [...]

  • ...With the invention of fast matchings-based algorithms [51], which improve the diagonal dominance of linear systems, the situation has dramatically changed and the impact on preconditioning methods [7], as well as the benefits for sparse direct solvers [60, 62], has been recognized....

    [...]

Journal ArticleDOI
TL;DR: This work revisits the sparse linear algebra computations of the parallel solver PIPS with the goal of improving the shared-memory performance and decreasing the time to solution.
Abstract: We present a scalable approach and implementation for solving stochastic optimization problems on high-performance computers. In this work we revisit the sparse linear algebra computations of the parallel solver PIPS with the goal of improving the shared-memory performance and decreasing the time to solution. These computations consist of solving sparse linear systems with multiple sparse right-hand sides and are needed in our Schur-complement decomposition approach to compute the contribution of each scenario to the Schur matrix. Our novel approach uses an incomplete augmented factorization implemented within the PARDISO linear solver and an outer BiCGStab iteration to efficiently absorb pivot perturbations occurring during factorization. This approach is capable of both efficiently using the cores inside a computational node and exploiting sparsity of the right-hand sides. We report on the performance of the approach on high-performance computers when solving stochastic unit commitment problems of unpre...

150 citations


Cites background or methods from "Matching-based preprocessing algori..."

  • ...A downside of the pivot perturbation that is not negligible is that solving with Ki might require an increased number of refinement steps to achieve the requested accuracy in the interior-point optimization algorithm [32]....

    [...]

  • ...To find an accurate solution during the Gaussian block elimination process we need to find a permutation matrix P that minimizes fill-in and avoids small pivots during the factorization [32]....

    [...]

  • ...These pivoting strategies are mainly based on graph-matching Bunch– Kaufman pivoting that trades fill-in versus stability to make efficient parallelization possible [32]....

    [...]

  • ...In this section we extend a PDE-constrained quadratic program that has been used in [32] to compare the triangular solves approach with the augmented approach....

    [...]

  • ...In many computational science applications, the numerical factorization phase Ki = LiDiLi of forming the partial Schur-complement contributions Si = B T i K −1 i Bi has generally received the most attention, because it is typically the largest component of execution time; most of the algorithmic improvements [7, 30, 32] in the factorization are related to the exploitation of the sparsity structure in Ki....

    [...]

Journal ArticleDOI
TL;DR: It is concluded that for a fixed number of unknowns and polynomial degree of approximation, a higher degree of continuity k drastically increases the CPU time and RAM needed to solve the problem when using a direct solver.

136 citations

References
More filters
Book
01 Nov 2008
TL;DR: Numerical Optimization presents a comprehensive and up-to-date description of the most effective methods in continuous optimization, responding to the growing interest in optimization in engineering, science, and business by focusing on the methods that are best suited to practical problems.
Abstract: Numerical Optimization presents a comprehensive and up-to-date description of the most effective methods in continuous optimization. It responds to the growing interest in optimization in engineering, science, and business by focusing on the methods that are best suited to practical problems. For this new edition the book has been thoroughly updated throughout. There are new chapters on nonlinear interior methods and derivative-free methods for optimization, both of which are used widely in practice and the focus of much current research. Because of the emphasis on practical methods, as well as the extensive illustrations and exercises, the book is accessible to a wide audience. It can be used as a graduate text in engineering, operations research, mathematics, computer science, and business. It also serves as a handbook for researchers and practitioners in the field. The authors have strived to produce a text that is pleasant to read, informative, and rigorous - one that reveals both the beautiful nature of the discipline and its practical side.

17,420 citations

Journal ArticleDOI
TL;DR: A comprehensive description of the primal-dual interior-point algorithm with a filter line-search method for nonlinear programming is provided, including the feasibility restoration phase for the filter method, second-order corrections, and inertia correction of the KKT matrix.
Abstract: We present a primal-dual interior-point algorithm with a filter line-search method for nonlinear programming. Local and global convergence properties of this method were analyzed in previous work. Here we provide a comprehensive description of the algorithm, including the feasibility restoration phase for the filter method, second-order corrections, and inertia correction of the KKT matrix. Heuristics are also considered that allow faster performance. This method has been implemented in the IPOPT code, which we demonstrate in a detailed numerical study based on 954 problems from the CUTEr test set. An evaluation is made of several line-search options, and a comparison is provided with two state-of-the-art interior-point codes for nonlinear programming.

7,966 citations

Journal ArticleDOI
TL;DR: This work presents a new coarsening heuristic (called heavy-edge heuristic) for which the size of the partition of the coarse graph is within a small factor of theSize of the final partition obtained after multilevel refinement, and presents a much faster variation of the Kernighan--Lin (KL) algorithm for refining during uncoarsening.
Abstract: Recently, a number of researchers have investigated a class of graph partitioning algorithms that reduce the size of the graph by collapsing vertices and edges, partition the smaller graph, and then uncoarsen it to construct a partition for the original graph [Bui and Jones, Proc. of the 6th SIAM Conference on Parallel Processing for Scientific Computing, 1993, 445--452; Hendrickson and Leland, A Multilevel Algorithm for Partitioning Graphs, Tech. report SAND 93-1301, Sandia National Laboratories, Albuquerque, NM, 1993]. From the early work it was clear that multilevel techniques held great promise; however, it was not known if they can be made to consistently produce high quality partitions for graphs arising in a wide range of application domains. We investigate the effectiveness of many different choices for all three phases: coarsening, partition of the coarsest graph, and refinement. In particular, we present a new coarsening heuristic (called heavy-edge heuristic) for which the size of the partition of the coarse graph is within a small factor of the size of the final partition obtained after multilevel refinement. We also present a much faster variation of the Kernighan--Lin (KL) algorithm for refining during uncoarsening. We test our scheme on a large number of graphs arising in various domains including finite element methods, linear programming, VLSI, and transportation. Our experiments show that our scheme produces partitions that are consistently better than those produced by spectral partitioning schemes in substantially smaller time. Also, when our scheme is used to compute fill-reducing orderings for sparse matrices, it produces orderings that have substantially smaller fill than the widely used multiple minimum degree algorithm.

5,629 citations

Journal ArticleDOI
TL;DR: It is shown that performance profiles combine the best features of other tools for performance evaluation to create a single tool for benchmarking and comparing optimization software.
Abstract: We propose performance profiles — distribution functions for a performance metric — as a tool for benchmarking and comparing optimization software. We show that performance profiles combine the best features of other tools for performance evaluation.

3,729 citations


Additional excerpts

  • ...1, using Dolan-Moré performance profiles [8]....

    [...]

Book
01 Jan 1993
TL;DR: An efficient translator is implemented that takes as input a linear AMPL model and associated data, and produces output suitable for standard linear programming optimizers.
Abstract: Practical large-scale mathematical programming involves more than just the application of an algorithm to minimize or maximize an objective function. Before any optimizing routine can be invoked, considerable effort must be expended to formulate the underlying model and to generate the requisite computational data structures. AMPL is a new language designed to make these steps easier and less error-prone. AMPL closely resembles the symbolic algebraic notation that many modelers use to describe mathematical programs, yet it is regular and formal enough to be processed by a computer system; it is particularly notable for the generality of its syntax and for the variety of its indexing operations. We have implemented an efficient translator that takes as input a linear AMPL model and associated data, and produces output suitable for standard linear programming optimizers. Both the language and the translator admit straightforward extensions to more general mathematical programs that incorporate nonlinear expressions or discrete variables.

3,176 citations


"Matching-based preprocessing algori..." refers methods in this paper

  • ...The first test set consists of 721 problems from the CUTE [5] collection, as provided by Benson [1] in the AMPL modeling language [20]....

    [...]

Frequently Asked Questions (20)
Q1. What contributions have the authors mentioned in the paper "Matching-based preprocessing algorithms to the solution of saddle-point problems in large-scale nonconvex interior-point optimization" ?

The authors present combinatorial methods to preprocess these matrices in order to establish more favorable numerical properties for the subsequent factorization. Their approach is based on symmetric weighted matchings and is used in a sparse direct LDL factorization method where the pivoting is restricted to static supernode data structures. The authors demonstrate the competitiveness of this approach within an interior-point method on a large set of test problems from the CUTE and COPS sets, as well as large optimal control problems based on partial differential equations. 

The basic idea of the symmetric a-priori pivoting is to form 2 × 2 diagonal blocks based on the matched off-diagonal entries of the matrix. 

In order to maintain the diagonal block structure, the authors compress the matrix by merging the structure of the rows and columns belonging to a 2 × 2 diagonal block. 

The idea introduced in [11] is to determine a bipartite maximum weight matching of the matrix, and to symmetrize it based on the cycle structure of the matrix. 

The vectors λ and z are the Lagrangian multipliers for the equality and bound constraints, and X and Z denote the diagonal matrices with the vector elements of x and z on the diagonal. 

In the pivoting variant II the authors will identify 2 × 2 and enforce the use of these preselected pivot blocks by adding additional zero elements to the structure of K . 

Duff and Pralet [13] showed that the problem of finding symmetric maximum weight matchings in a symmetric matrix is equivalent to finding a non-bipartite matching. 

By matching in AT it is possible to detect certain structural rank deficiencies, which can, for example, stem from redundant constraints in the optimization problem formulation. 

The solvers used for the comparisons is the PARDISO solver with the different ordering strategies, as well as the Harwell routines MA27 and MA57 (Version 3.0) [10]. 

4.2 Matching variant: constraint matching in ATIf the authors focus on finding favorable 2 × 2 diagonal pivots and assume that only diagonal entries are matched in H , the authors can approximate a global symmetric maximum weight matching by determining a row-perfect maximum weight matching in the constraint block AT only. 

The heuristic for choosing δx and δc implemented in IPOPT is described in [41]; it usually attempts first to work with δx = δc = 0 to obtain the pure primal–dual Newton direction for fast convergence, and only if the inertia is not correct with this choice are other values tried. 

In their case the authors are interested in a perfect maximum weight matching, i.e., a matching with maximum weight under the condition that all rows and columns are matched. 

Using the fact that the Jacobian of the system of Eq. (2.1) is non-singular at a non-degenerate local solution (x∗, λ∗, z∗) of (1.1), it is possible to derive methods that converge superlinearly toward (x∗, λ∗, z∗); see, e.g., [23]. 

Since more than 60% of the CUTE problems are solved in less than 0.1 CPU seconds, the authors do not present a comparison of the runtimes. 

One reason for this surge in interest is due the success of interior-point methods in nonlinear programming, which at their core require the solution of a series of linear systems in saddle-point form. 

In the context of an SQP-type line search optimization method it is therefore possible to ensure the required descent properties by performing a “trial” factorization for systems of the form (2.4), where the matrix W̃k is modified according to some heuristics, until the inertia of the matrix is correct. 

closer examination of the results shows that NLPs, which are degenerate so that the iteration matrix is singular in every iteration, are not always handled well, because the heuristics implemented in IPOPT to detect such structural degeneracy rely on the detection of singularity of the matrix by the linear solver. 

In their final experiments the authors explore the parallel scale-up of the PARDISO code (only using the PARDISO2 method), solving the discretized PDE optimal control problems for instances with several million variables. 

In order to enforce the use of these preselected pivot blocks during the LDLT factorization, the authors merge the 2 × 2 column structure in L in such a way that these initial 2 × 2 pivot structure is maintained in the factor L.This is illustrated in Fig. 

Since these general matching algorithms are much more expensive than bipartite matching algorithms, the symmetric maximum weight matchings are typically only approximated.