scispace - formally typeset
Search or ask a question

Showing papers in "Optimization Methods & Software in 2012"


Journal ArticleDOI
TL;DR: The basic components and the underlying philosophy of ADMB are described, with an emphasis on functionality found in no other statistical software, and the main advantages are flexibility, speed, precision, stability and built-in methods to quantify uncertainty.
Abstract: Many criteria for statistical parameter estimation, such as maximum likelihood, are formulated as a nonlinear optimization problem. Automatic Differentiation Model Builder (ADMB) is a programming framework based on automatic differentiation, aimed at highly nonlinear models with a large number of parameters. The benefits of using AD are computational efficiency and high numerical accuracy, both crucial in many practical problems. We describe the basic components and the underlying philosophy of ADMB, with an emphasis on functionality found in no other statistical software. One example of such a feature is the generic implementation of Laplace approximation of high-dimensional integrals for use in latent variable models. We also review the literature in which ADMB has been used, and discuss future development of ADMB as an open source project. Overall, the main advantages of ADMB are flexibility, speed, precision, stability and built-in methods to quantify uncertainty.

1,753 citations


Journal ArticleDOI
TL;DR: The robust optimization framework in the modelling language YALMIP is presented, which carries out robust modelling and uncertainty elimination automatically and allows the user to concentrate on the high-level model.
Abstract: This paper presents the robust optimization framework in the modelling language YALMIP, which carries out robust modelling and uncertainty elimination automatically and allows the user to concentrate on the high-level model. While introducing the software package, a brief summary of robust optimization is given, as well as some comments on modelling and tractability of complex convex uncertain optimization problems.

210 citations


Journal ArticleDOI
TL;DR: This paper studies a generalization of the distributed optimization problem, and discusses an algorithm and proves its convergence, and then discusses extensions to more general and complex distributed optimization problems.
Abstract: In a distributed optimization problem, the complete problem information is not available at a single location but is rather distributed among different agents in a multi-agent system. In the problems studied in the literature, each agent has an objective function and the network goal is to minimize the sum of the agents’ objective functions over a constraint set that is globally known. In this paper, we study a generalization of the above distributed optimization problem. In particular, the network objective is to minimize a function of the sum of the individual objective functions over the constraint set. The ‘outer’ function and the constraint set are known to all the agents. We discuss an algorithm and prove its convergence, and then discuss extensions to more general and complex distributed optimization problems. We provide a motivation for our algorithms through the example of distributed regression of distributed data.

117 citations


Journal ArticleDOI
TL;DR: A new interior-point method is proposed, which is based on an extension of the ideas of self-scaled optimization to the general cases, using the primal correction process to find a scaling point and using this point to compute a strictly feasible primal–dual pair by simple projection.
Abstract: In this paper, we propose a new interior-point method, which is based on an extension of the ideas of self-scaled optimization to the general cases. We suggest using the primal correction process to find a scaling point . This point is used to compute a strictly feasible primal–dual pair by simple projection. Then, we define an affine-scaling direction and perform a prediction step. This is the only moment when the dual barrier is used. Thus, we need only to compute its value, which can even be done approximately. In the second part of the paper, we develop a 4 n -self-concordant barrier for an n -dimensional p -cone, which can be used for numerical testing of the proposed technique.

79 citations


Journal ArticleDOI
TL;DR: Regular control problems in the sense of the Legendre condition are defined, and second-order necessary and sufficient optimality conditions in this class are reviewed.
Abstract: Regular control problems in the sense of the Legendre condition are defined, and second-order necessary and sufficient optimality conditions in this class are reviewed. Adding a scalar homotopy parameter, differential pathfollowing is introduced. The previous sufficient conditions ensure the definiteness and regularity of the path. The role of AD for the implementation of this approach is discussed, and two examples excerpted from quantum and space mechanics are detailed.

66 citations


Journal ArticleDOI
TL;DR: NMLS is able to identify the zero components of a stationary point after a finite number of steps under some mild conditions and the global convergence of FPC_AS is established based on the properties of NMLS.
Abstract: We analyse an abridged version of the active-set algorithm FPC_AS proposed in Wen et al . [A fast algorithm for sparse reconstruction based on shrinkage, subspace optimisation and continuation, SIAM J. Sci. Comput. 32 2010, pp. 1832–1857] for solving the l 1-regularized problem, i.e. a weighted sum of the l 1-norm | x |1 and a smooth function f x . The active-set algorithm alternatively iterates between two stages. In the first ‘non-monotone line search NMLS’ stage, an iterative first-order method based on ‘shrinkage’ is used to estimate the support at the solution. In the second ‘subspace optimization’ stage, a smaller smooth problem is solved to recover the magnitudes of the non-zero components of x . We show that NMLS itself is globally convergent and the convergence rate is at least R-linearly. In particular, NMLS is able to identify the zero components of a stationary point after a finite number of steps under some mild conditions. The global convergence of FPC_AS is established based on the properties of NMLS.

58 citations


Journal ArticleDOI
TL;DR: A hierarchy of increasingly better outer polyhedral approximations to the copositive cone is proposed, and it is established that the sequence of approxIMations is exact in the limit.
Abstract: We consider linear optimization problems over the cone of copositive matrices. Such conic optimization problems, called copositive programs, arise from the reformulation of a wide variety of difficult optimization problems. We propose a hierarchy of increasingly better outer polyhedral approximations to the copositive cone. We establish that the sequence of approximations is exact in the limit. By combining our outer polyhedral approximations with the inner polyhedral approximations due to de Klerk and Pasechnik [SIAM J. Optim. 12 (2002), pp. 875–892], we obtain a sequence of increasingly sharper lower and upper bounds on the optimal value of a copositive program. Under primal and dual regularity assumptions, we establish that both sequences converge to the optimal value. For standard quadratic optimization problems, we derive tight bounds on the gap between the upper and lower bounds. We provide closed-form expressions of the bounds for the maximum stable set problem. Our computational results shed light...

57 citations


Journal ArticleDOI
TL;DR: The interval inequality relations and the concept of solutions of the matrix game with interval payoffs are defined and it is shown that the models proposed extend those of the classical matrix game.
Abstract: The aim of this paper is to study how to solve a type of matrix game with interval payoffs. In this paper, the interval inequality relations and the concept of solutions of the matrix game with interval payoffs are defined. Based on the fuzzy ranking index defined, the solution of the matrix game with interval payoffs can be obtained through solving a pair of bi-objective linear programming models derived from the constructed auxiliary interval programming models. It is shown that the models proposed in this paper extend those of the classical matrix game. The validity and applicability of the proposed methodology are illustrated with a numerical example.

50 citations


Journal ArticleDOI
TL;DR: A GPU-accelerated PSO (GPSO) algorithm is proposed by using a thread pool model and implement GPSO on a GPU, a promising method for tackling high-dimensional and difficult optimization problems using a low-cost and many-core GPU system.
Abstract: Particle swarm optimization (PSO) is a population-based stochastic and derivative-free method that has been used to solve various optimization problems due to its simplicity and efficiency. While solving high-dimensional or complicated problems, PSO requires a large number of particles to explore the problem domains and consequently introduces high computational costs. In this paper, we focus on the acceleration of PSO for solving box-constrained, load-balanced optimization problems by parallelization on a graphics processing unit (GPU). We propose a GPU-accelerated PSO (GPSO) algorithm by using a thread pool model and implement GPSO on a GPU. Numerical results show that the GPU architecture fits the PSO framework well by reducing computational timing, achieving high parallel efficiency and finding better optimal solutions by using a large number of particles. For example, while solving the 100-dimensional test problems with 65,536 (16×212) particles, GPSO has achieved up to 280X and 83X speedups on a NVI...

47 citations


Journal ArticleDOI
TL;DR: This work test and compare different methods from both subgradient methods and bundle methods as well as some methods which may be considered as hybrids of these two and/or some others to get some insight on which method to select for certain types of problems.
Abstract: Most nonsmooth optimization (NSO) methods can be divided into two main groups: subgradient methods and bundle methods. In this paper, we test and compare different methods from both groups as well as some methods which may be considered as hybrids of these two and/or some others. All the solvers tested are so-called general black box methods which, at least in theory, can be applied to solve almost all NSO problems. The test set includes a large number of unconstrained nonsmooth convex and nonconvex problems of different size. In particular, it includes piecewise linear and quadratic problems. The aim of this work is not to foreground some methods over the others but to get some insight on which method to select for certain types of problems.

47 citations


Journal ArticleDOI
TL;DR: The aim of this paper is to improve the convergence theory from Steffensen and Ulbrich and translate this local regularization idea to MPVCs and obtain a new solution method for this class of optimization problems for which several convergence results are given.
Abstract: Mathematical programmes with equilibrium or vanishing constraints MPECs or MPVCs are both known to be difficult optimization problems which typically violate all standard constraint qualifications. A number of methods try to exploit the particular structure of MPECs and MPVCs in order to overcome these difficulties. In a recent paper by Steffensen and Ulbrich S. Steffensen and M. Ulbrich, A new regularization scheme for mathematical programs with equilibrium constraints , SIAM J. Optim. 2010., this was done for MPECs by a local regularization idea that may be viewed as a modification of the popular global regularization technique by Scholtes S. Scholtes, Convergence properties of a regularization scheme for mathematical programs with complementarity constraints , SIAM J. Optim. 11 2001, pp. 918–936.. The aim of this paper is twofold. First, we improve the convergence theory from S. Steffensen and M. Ulbrich, A new regularization scheme for mathematical programs with equilibrium constraints , SIAM J. Optim. 2010. in the MPEC setting, and second we translate this local regularization idea to MPVCs and obtain a new solution method for this class of optimization problems for which several convergence results are given.

Journal ArticleDOI
TL;DR: A randomized initialization technique that involves a nonlinear version of the classical multidimensional scaling, and a dimensionality relaxation scheme with optional weighting is developed.
Abstract: Our experiments show that the method easily solves the artificial problems introduced by More and Wu. It also solves the 12 much more difficult protein fragment problems introduced by Hendrickson, and the six larger protein problems introduced by Grooms, Lewis and Trosset.

Journal ArticleDOI
TL;DR: This work combines the numerical concept of variational discretization and semi-smooth Newton methods for the numerical solution of pde-constrained optimization with control constraints and proves fast local convergence of a globalized algorithm.
Abstract: Combining the numerical concept of variational discretization and semi-smooth Newton methods for the numerical solution of pde-constrained optimization with control constraints, we place special emphasis on the implementation and globalization of the numerical algorithm. We prove fast local convergence of a globalized algorithm and illustrate our analytical and algorithmical findings by numerical experiments.

Journal ArticleDOI
TL;DR: The authors' complexity upper bounds match in order (as a function of the accuracy of approximation), and sometimes even improve, those obtained by Nesterov and Polyak for these same problem classes, without requiring exact Hessians or exact or global solution of the subproblem.
Abstract: The adaptive cubic regularization algorithms described in Cartis, Gould and Toint [Adaptive cubic regularisation methods for unconstrained optimization Part II: Worst-case function- and derivative-evaluation complexity, Math. Program. (2010), doi:10.1007/s10107-009-0337-y (online)]; [Part I: Motivation, convergence and numerical results, Math. Program. 127(2) (2011), pp. 245–295] for unconstrained (nonconvex) optimization are shown to have improved worst-case efficiency in terms of the function- and gradient-evaluation count when applied to convex and strongly convex objectives. In particular, our complexity upper bounds match in order (as a function of the accuracy of approximation), and sometimes even improve, those obtained by Nesterov [Introductory Lectures on Convex Optimization, Kluwer Academic Publishers, Dordrecht, 2004; Accelerating the cubic regularization of Newton's method on convex problems, Math. Program. 112(1) (2008), pp. 159–181] and Nesterov and Polyak [Cubic regularization of Newton's m...

Journal ArticleDOI
TL;DR: It is proved that under more natural assumptions than the ones employed until now, penalty parameters are bounded and the original algorithm has been slightly modified for proving the new boundedness result.
Abstract: Augmented Lagrangian methods are effective tools for solving large-scale nonlinear programming problems. At each outer iteration, a minimization subproblem with simple constraints, whose objective function depends on updated Lagrange multipliers and penalty parameters, is approximately solved. When the penalty parameter becomes very large, solving the subproblem becomes difficult; therefore, the effectiveness of this approach is associated with the boundedness of the penalty parameters. In this paper, it is proved that under more natural assumptions than the ones employed until now, penalty parameters are bounded. For proving the new boundedness result, the original algorithm has been slightly modified. Numerical consequences of the modifications are discussed and computational experiments are presented.

Journal ArticleDOI
TL;DR: The mathematical foundations of SOA sensitivity analysis are discussed and it is shown that it provides an efficient approach to obtain Hessian-vector products in large-scale data assimilation problems.
Abstract: In this paper, we discuss the mathematical foundations of SOA sensitivity analysis and show that it provides an efficient approach to obtain Hessian-vector products. We study the benefits of using second-order information in the numerical optimization process for data assimilation applications. The numerical studies are performed in a twin experiment setting with a two-dimensional shallow water model. Different scenarios are considered with different discretization approaches, observation sets, and noise levels. Optimization algorithms that employ second-order derivatives are tested against widely used methods that require only first-order derivatives. Conclusions are drawn regarding the potential benefits and the limitations of using high-order information in large-scale data assimilation problems.

Journal ArticleDOI
TL;DR: By exploiting the properties of a (minimum) flow network, and identifying pleasing properties of some maximum cuts, an elegant single-pass optimal ordinal relabelling algorithm is formulated.
Abstract: Noise in multi-criteria data sets can manifest itself as non-monotonicity. Work on the remediation of such non-monotonicity is rather scarce. Nevertheless, errors are often present in real-life data sets, and several monotone classification algorithms are unable to use such partially non-monotone data sets. Fortunately, as we will show here, it is possible to restore monotonicity in an optimal way, by relabelling part of the data set. By exploiting the properties of a (minimum) flow network, and identifying pleasing properties of some maximum cuts, an elegant single-pass optimal ordinal relabelling algorithm is formulated.

Journal ArticleDOI
TL;DR: A non-convex loss function is proposed to construct a robust support vector regression (SVR) and a Newton-type algorithm is developed to solve it, which can both retain the sparseness of SVR and oppress outliers in the training samples.
Abstract: The classical support vector machines are constructed based on convex loss functions. Recently, support vector machines with non-convex loss functions have attracted much attention for their superiority to the classical ones in generalization accuracy and robustness. In this paper, we propose a non-convex loss function to construct a robust support vector regression SVR. The introduced non-convex loss function includes several truncated loss functions as its special cases. The resultant optimization problem is a difference of convex functions program. We employ the concave–convex procedure and develop a Newton-type algorithm to solve it, which can both retain the sparseness of SVR and oppress outliers in the training samples. The experiments on both synthetic and real-world benchmark data sets confirm the robustness and effectiveness of the proposed method.

Journal ArticleDOI
TL;DR: This paper presents as global optimization algorithm a specific version of the monotonic basin hopping method which employs, as local minimizer, an efficient versions of the Frank–Wolfe method which shows the effectiveness of the presented methodology in terms of global optimization.
Abstract: This paper considers a portfolio selection problem in which portfolios with minimum number of active assets are sought. This problem is motivated by the need of inducing sparsity on the selected portfolio to reduce transaction costs, complexity of portfolio management, and instability of the solution. The resulting problem is a difficult combinatorial problem. We propose an approach based on the definition of an equivalent smooth concave problem. In this way, we move the difficulty of the original problem to that of solving a concave global minimization problem. We present as global optimization algorithm a specific version of the monotonic basin hopping method which employs, as local minimizer, an efficient version of the Frank–Wolfe method. We test our method on various data sets of small, medium, and large dimensions involving real-world capital market from major stock markets. The obtained results show the effectiveness of the presented methodology in terms of global optimization. Furthermore, also the out-of-sample performances of the selected portfolios, as measured by Sharpe ratio, appear satisfactory.

Journal ArticleDOI
TL;DR: This interdisciplinary paper reviews higher-order numerical methods for the solution of nonlinear problems, and proposes a synthesis of two different conceptual frameworks, namely automatic differentiation and the asymptotic numerical method.
Abstract: Modelling often involves nonlinear parametric problems and bifurcation analysis. This interdisciplinary paper reviews higher-order numerical methods for the solution of nonlinear problems, and proposes a synthesis of two different conceptual frameworks, namely automatic differentiation and the asymptotic numerical method. Various mechanical problems and references illustrate the presentation.

Journal ArticleDOI
TL;DR: In this article, the robust optimization of simulation-based functions is studied, where the robust counterpart of an optimization problem without derivatives falls into the category of the bilevel problems under consideration here.
Abstract: One important application of our work appears in the robust optimization of simulation-based functions, which may arise due to implementation variables or uncertain parameters. The robust counterpart of an optimization problem without derivatives falls into the category of the bilevel problems under consideration here. We provide numerical illustrations of the application of our algorithmic framework to such robust optimization examples.

Journal ArticleDOI
TL;DR: An interior-point method (IPM) for Cartesian P *(κ)- linear complementarity problems over symmetric cones (SCLCP) is analysed and the complexity results are presented.
Abstract: An interior-point method IPM for Cartesian P *κ-linear complementarity problems over symmetric cones SCLCP is analysed and the complexity results are presented. The Cartesian P *κ-SCLCPs have been recently introduced as the generalization of the more commonly known and more widely used monotone-SCLCPs. The IPM is based on the barrier functions that are defined by a large class of univariate functions called eligible kernel functions, which have recently been successfully used to design new IPMs for various optimization problems. Eligible barrier kernel functions are used in calculating the Nesterov–Todd search directions and the default step-size which lead to very good complexity results for the method. For some specific eligible kernel functions, we match the best-known iteration bound for the long-step methods while for the short-step methods the best iteration bound is matched for all cases.

Journal ArticleDOI
TL;DR: This work investigates the computation of Hessian matrices via Automatic Differentiation, using a graph model and an algebraic model, and develops edge_pushing, a new truly reverse Hessian computation algorithm that fully exploits the Hessian's symmetry.
Abstract: We investigate the computation of Hessian matrices via Automatic Differentiation, using a graph model and an algebraic model. The graph model reveals the inherent symmetries involved in calculating the Hessian. The algebraic model, based on Griewank and Walther's [Evaluating derivatives, in Principles and Techniques of Algorithmic Differentiation, 2nd ed., Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2008] state transformations synthesizes the calculation of the Hessian as a formula. These dual points of view, graphical and algebraic, lead to a new framework for Hessian computation. This is illustrated by developing edge_pushing, a new truly reverse Hessian computation algorithm that fully exploits the Hessian's symmetry. Computational experiments compare the performance of edge_pushing on 16 functions from the CUTE collection [I. Bongartz et al. Cute: constrained and unconstrained testing environment, ACM Trans. Math. Softw. 21(1) (1995), pp. 123–160] against two algorithms av...

Journal ArticleDOI
TL;DR: It is proved that a numerical trajectory converges to an optimal trajectory of the continuous-time control problem as the time step goes to zero, with both the limiting optimal state and costate trajectories being absolutely continuous.
Abstract: This paper presents a numerical scheme for solving the continuous-time convex linear-quadratic LQ optimal control problem with mixed polyhedral state and control constraints. Unifying a discretization of this optimal control problem as often employed in model predictive control and that obtained through time-stepping methods based on the differential variational inequality reformulation, the scheme solves a sequence of finite-dimensional convex quadratic programs QPs whose optimal solutions are employed to construct a sequence of discrete-time trajectories dependent on the time step. Under certain technical primal–dual assumptions primarily to deal with the algebraic constraints involving the state variable, we prove that such a numerical trajectory converges to an optimal trajectory of the continuous-time control problem as the time step goes to zero, with both the limiting optimal state and costate trajectories being absolutely continuous. This provides a constructive proof of the existence of a solution to the optimal control problem with such regularity properties. Additional properties of the optimal solutions to the LQ problem are also established that are analogous to those of the finite-dimensional convex QP. Our results are applicable to problems with convex but not necessarily strictly convex objective functions and with possibly unbounded mixed state–control constraints.

Journal ArticleDOI
TL;DR: A specialized parallel factorization procedure for dense saddle-point linear systems in a distributed-memory environment, together with a streamlined method for assembling the distributed dense matrix is developed.
Abstract: We present a novel approach for solving dense saddle-point linear systems in a distributed-memory environment. This work is motivated by an application in stochastic optimization problems with recourse, but the proposed approach can be used for a large family of dense saddle-point systems, in particular, for those arising in convex programming. Although stochastic optimization problems have many important applications, they can present serious computational difficulties. In particular, sample average approximation SAA problems with a large number of samples are often too big to solve on a single shared-memory system. Recent work has developed interior-point methods and specialized linear algebra to solve these problems in parallel, using a scenario-based decomposition that distributes the data, and work across computational nodes. Even for sparse SAA problems, the decomposition produces a dense and possibly very large saddle-point linear system that must be solved repeatedly. We developed a specialized parallel factorization procedure for these systems, together with a streamlined method for assembling the distributed dense matrix. Strong scaling tests indicate over 90% efficiency on 1024 cores on a stochastic unit commitment problem with 57 million variables. Stochastic unit commitment problems with up to 189 million variables are solved efficiently on up to 2048 cores.

Journal ArticleDOI
TL;DR: This paper considers an extension of ordinary linear programming that adds weighted logarithmic barrier terms for some variables and shows that the problem has a dual of the same form and gives complexity results for several different interior-point algorithms.
Abstract: In this paper, we consider an extension of ordinary linear programming LP that adds weighted logarithmic barrier terms for some variables. The resulting problem generalizes both LP and the problem of finding the weighted analytic centre of a polytope. We show that the problem has a dual of the same form and give complexity results for several different interior-point algorithms. We obtain an improved complexity result for certain cases by utilizing a combination of the volumetric and logarithmic barriers. As an application, we consider the complexity of solving the Eisenberg–Gale formulation of a Fisher equilibrium problem with linear utility functions.

Journal ArticleDOI
TL;DR: Two strategies are developed to show the trade-off of the optimal objective function value with tightening or loosening general constraints, and a simple method which may be performed immediately after a single optimization and a detailed method performing biobjective optimization on the minimization of the objective versus a constraint of interest.
Abstract: This paper proposes a framework for trade-off analyses of blackbox constrained optimization problems. Two strategies are developed to show the trade-off of the optimal objective function value with tightening or loosening general constraints. These are a simple method which may be performed immediately after a single optimization and a detailed method performing biobjective optimization on the minimization of the objective versus a constraint of interest. The detailed method provides points on the Pareto front, the trade-off curve, of the objective versus a chosen constraint. The simple method provides points near the trade-off curve, which may be all the designer needs. The trade-off information is generally used by engineers rather than the first-order sensitivity estimates provided by the Lagrange multipliers, which only provide the tangent to the Pareto front at the solution found. The proposed methods are tested on an academic test case and on an engineering problem using the mesh-adaptive direct search algorithm.

Journal ArticleDOI
TL;DR: This paper introduces data structures for storing the third derivative (tensor) of a multivariate scalar function when the second and third derivatives of the function are sparse and gives the number of arithmetic operations of the methods of the Halley class compared to Newton's method.
Abstract: In this paper, we introduce data structures for storing the third derivative (tensor) of a multivariate scalar function when the second and third derivatives of the function are sparse. We consider solving unconstrained optimization problems using methods of the Halley class. The gradient, Hessian matrix and the tensor of the objective function are required for the methods and they have a third-order rate of convergence. However, these methods require the solution of two linear systems of equations at each iteration. Solving the linear systems using a factorization will cause fill-ins. The impact of these fill-ins will be discussed and we give the number of arithmetic operations of the methods of the Halley class compared to Newton's method. Finally, we will give some preliminary numerical results showing that including the time to compute the derivatives, the ratios of elapsed time of methods in the Halley class and Newton's method are close to the ratios of arithmetic operations.

Journal ArticleDOI
TL;DR: The interior-point algorithm described in Peng et al. is adapted to symmetric optimization for linear optimization based on kernel functions and a complexity bound is derived.
Abstract: We present a generalization to symmetric optimization of interior-point methods for linear optimization based on kernel functions. Symmetric optimization covers the three most common conic optimization problems: linear, second-order cone and semi-definite optimization problems. Namely, we adapt the interior-point algorithm described in Peng et al . [ Self-regularity: A New Paradigm for Primal–Dual Interior-point Algorithms . Princeton University Press, Princeton, NJ, 2002.] for linear optimization to symmetric optimization. The analysis is performed through Euclidean Jordan algebraic tools and a complexity bound is derived.

Journal ArticleDOI
TL;DR: A pivoting algorithm for solving linear programs with linear complementarity constraints, which generalizes the simplex method for linear programming to deal with complementarity conditions, and an anticycling scheme that can verify Bouligand stationarity is developed.
Abstract: We present a pivoting algorithm for solving linear programs with linear complementarity constraints. Our method generalizes the simplex method for linear programming to deal with complementarity conditions. We develop an anticycling scheme that can verify Bouligand stationarity. We also give an optimization-based technique to find an initial feasible vertex. Starting with a feasible vertex, our algorithm always finds a minimizer or an unbounded descent search direction in a finite number of pivoting steps.