Feasibility Pump (FP) is a well-studied mixed integer linear programming primal heuristic that sequentially guides fractional solutions into integer feasibility by minimizing their distance to given target integer points, while preserving the linear constraints of the original problem. Due to its effectiveness and simplicity, FP has been embedded in most commercial solvers, and several variants of the basic method have been proposed, aiming at improving its convergence and the quality of the solutions obtained. This work introduces Kernel Pump (KP), a novel FP speed-up technique that relies on the intuition that not all binary variables are active (i.e., set to one) in a feasible solution, and can, thus, be sub-divided in a kernel search fashion. Precisely, KP identifies a subset of promising binary variables (the kernel) and distributes the rest into buckets ranked by their likelihood of composing a feasible/optimal solution. FP sub-problems are then solved iteratively, focusing on the current kernel and progressively adding variables from the buckets whenever necessary. Extensive computational experiments on benchmark instances from MIPLIB 2017, as well as two additional problems from the literature of vehicle routing and scheduling, reveal a considerable boost in success rate, also decreasing running times for the majority of instances tested. Notably, aside from outperforming FP, the new method was able to find feasible solutions in cases where CPLEX failed within an hour of search. To our knowledge, this is also the first study to introduce a decomposition approach within the FP framework.
Probing and dual fixing are two powerful presolve techniques in mixed integer programming (MIP) solvers. Probing tentatively sets some binary variables to 0 or 1, applies linear constraint based domain propagation techniques to derive better variable bounds, and extracts useful information such as stronger variable implications and better global variable bounds. Dual fixing attempts to fix variables to lower or upper bounds while ensuring that at least one optimal solution is retained, as long as the problem was feasible. In this paper, we investigate how to combine the two approaches in MIP solvers to achieve a better performance. In particular, we first embed dual fixing into the probing framework, deriving more useful variables' implications for enhancing the capability of probing. Then, we develop an improved dual fixing technique where more variable fixings can be applied, and use the probing framework to detect the reductions. Computational results on the MIPLIB 2017 benchmark instances demonstrate the potential of the two proposed techniques in combining probing and dual fixing on the open-source MIP solver HiGHS.
In this paper, we study a class of structured sparse optimization problems characterized by a convex, possibly non-smooth loss function and a capped- ℓ _1 penalty. This model can provide an exact continuous relaxation for the problems with a cardinality penalty. First, we propose an Accelerated Nested Proximal Gradient (ANPG) algorithm, which employs a nested proximal structure and selective extrapolation to enhance computational efficiency. Under some mild and easily verifiable conditions, we prove the subsequence convergence of the ANPG algorithm to the lifted stationary points of the considered problem, which correspond to the strong local minimizers of the corresponding cardinality penalty problems. Moreover, we establish an O(k^-1) convergence rate on the objective function values, while a refined extrapolation strategy ensures sequence convergence of the iterates, albeit with a potential reduction in the theoretical convergence rate on the objective values. Second, for the cases with a smooth loss function, we further propose an Accelerated Proximal Gradient (APG) algorithm that guarantees the sequence convergence on the iterates and achieves a faster convergence rate of o(k^-2) in terms of objective values. Finally, the effectiveness of the ANPG and APG algorithms, as well as the high quality of the solutions they produce, are verified through numerical experiments on the least absolute deviation regression problem and the sparse logistic regression problem, respectively.
We consider several hierarchical optimization problems: (generalized) semi-infinite and existence-constrained semi-infinite programs, minmax, and bilevel programs. Multiple adaptive discretization-based algorithms have been published for these problem classes in recent decades. However, rigorous numerical performance comparisons between these algorithms are lacking. Indeed, if numerical comparisons are provided at all, they typically compare a small selection of algorithms on small test sets, across different platforms, and with various subsolvers required during the solution. Additionally, some algorithms have parameters that impede a fair comparison. Our contribution is threefold: i) We present an open-source software called libDIPS (Discretization-based semi-Infinite and bilevel Programming Solvers), which implements multiple adaptive discretization-based solvers. The primary benefit of libDIPS is that it enables users to flexibly switch between the implemented solvers within a single problem class and switch between the available subsolvers. ii) We compile an extensive library of test problems for the (generalized) semi-infinite, minmax, and bilevel problem classes, which, in total, contains over 600 problem instances. Our set includes eight merged test sets and additional problem instances from over 80 literature sources. iii) We compare the solvers numerically using our library of test problems and identify tradeoffs in the parameter tuning.
Linear optimization solvers commonly employ presolve techniques to simplify and improve the solution characteristics of models before a simplex or interior point algorithm is used for solution. The Fourier-Motzkin elimination, a well-known mathematical algorithm for removing columns from linear inequality systems, has not been utilized as a presolve method in existing solvers due to its exponential row growth. We propose a novel presolve method based on the Fourier-Motzkin elimination, along with postsolve algorithms to recover the primal and dual optimal solutions. The algorithm controls row growth by judiciously selecting which columns to eliminate after computing an upper bound on model size prior to each variable elimination step. Our computational results demonstrate that the proposed presolve with Fourier-Motzkin elimination effectively reduces the number of rows, columns, and nonzero elements in general linear optimization problems. As a result, the algorithm reduces the total number of iterations of simplex methods in CPLEX. Additionally, when combined with the CPLEX presolve, our FME implementation leads to models smaller than those obtained by CPLEX presolve alone, leading to speed-ups for CPLEX’s primal simplex, dual simplex, and barrier methods.
We present a GPU-accelerated backend for QOCO, a C-based solver for quadratic objective second-order cone programs (SOCPs) based on a primal-dual interior point method. Our backend uses NVIDIA’s cuDSS library to perform a direct sparse LDL^⊤ factorization of the KKT system at each iteration. We also develop custom CUDA kernels for cone operations and show that parallelizing these operations is essential for achieving peak performance. Additionally, we refactor QOCO to introduce a modular backend abstraction that decouples solver logic from the underlying linear algebra implementations, allowing the existing CPU and new GPU backend to share a unified codebase. This GPU backend is accessible through a direct Python interface and through CVXPY, allowing for easy use. Numerical experiments on a range of large-scale quadratic programs and SOCPs with tens to hundreds of millions of nonzero elements in the KKT matrix, demonstrate speedups of up to 50–70 times over the CPU implementation.
In this report, we present the BOBILib, a collection of more than 2600 instances of mixed integer bilevel linear optimization problems (MIBLPs). The goal of this library is to provide a large and well-curated set of test instances freely available for the research community so that new and existing algorithms in bilevel optimization can be tested and compared in a standardized way. The library is sub-divided into instances of different types and also contains different benchmark instance sets. Moreover, we present a new data format for MIBLPs that is less error-prone compared to an older format that will now be deprecated. We provide numerical results for all instances of the library using available bilevel solvers. Based on these numerical results, we select benchmark instance sets, which provide a meaningful basis for experimental comparisons of solution methods in a moderate time. The instances, together with solution files, can be downloaded at https://bobilib.org .
Adjustable robust optimization (ARO) is a powerful tool to model problems that have uncertain data and that feature a two-stage decision-making process. Computationally, they are often addressed using the column-and-constraint generation (CCG) algorithm introduced by Zeng and Zhao [45]. While it was empirically shown that the algorithm scales well if all second-stage decisions are continuous, the presence of integer variables in the second stage rapidly leads to challenging large-scale mixed-integer problems within CCG. These problems can no longer be solved to global optimality within reasonable time limits in general. In this work, we explicitly focus on ARO problems with mixed-integer second-stage decisions and discuss the main difficulties of successfully applying CCG to this problem class. We then introduce, for a large set of problems with specific structural properties, a stronger formulation, which can be used in place of the master problem in the classic CCG algorithm. We show how this model can be effectively solved by column generation (CG). Additionally, we introduce a new CG-based heuristic that is able to generate new feasible points to speed up the overall method. We apply this nested scheme, combining CCG and CG, to three problems from logistics and scheduling to show the applicability of our approach.
We propose a primal heuristic for quadratic mixed-integer problems. Our method extends the Boscia framework – originally a mixed-integer convex solver leveraging a Frank-Wolfe-based branch-and-bound approach – to address nonconvex quadratic objective functions and constraints. We reformulate nonlinear constraints, introduce preprocessing steps, and a suite of heuristics including rounding strategies, gradient-guided selection, and large neighborhood search techniques that exploit integer-feasible vertices generated during the Frank-Wolfe iterations. Computational results demonstrate the effectiveness of our method in solving challenging MIQCQPs, achieving improvements on QPLIB instances within minutes and winning first place in the Land-Doig MIP Computational Competition 2025.
The market split problem was proposed by Cornuéjols and Dawande in 1998 as benchmark problem for algorithms solving linear systems with binary variables. The recent (2025) Quantum Optimization Benchmark Library (QOBLIB) contains a set of feasible instances of the market split problem. In QOBLIB an instance of the market split problem is considered as solved as soon as at least one feasible solution has been found. The market split problem seems to be difficult to solve with the conventional branch-and-cut approach of integer linear programming software which reportedly can handle QOBLIB instances up to m=7 . In contrast, a new GPU implementation of the Schroeppel–Shamir algorithm solves instances up to m=11 . In this note we report about experiments with an algorithm that reduces the market split problem to a lattice problem. With the author’s most recent implementation – named solvediophant – instances of the QOBLIB market split benchmark problems can be solved up to m=14 on a standard computer.
Doubly nonnegative (DNN) relaxation usually provides a tight lower bound for a mixed-binary quadratic program (MBQP). However, solving DNN problems is challenging because: (1) the problem size is Ω((n+l)^2) for an MBQP with n variables and l inequality constraints, and (2) the rank of optimal solutions cannot be estimated a priori due to the absence of theoretical bounds. In this work, we propose RiNNAL+, a Riemannian augmented Lagrangian method (ALM) for solving DNN problems. We prove that the DNN relaxation of an MBQP, with matrix dimension (n+l+1), is equivalent to the SDP-RLT relaxation (based on the reformulation-linearization technique) with a smaller matrix dimension (n+1). In addition, we develop a hybrid method that alternates between two phases to solve the ALM subproblems. In phase one, we apply low-rank matrix factorization and random perturbation to transform the feasible region into a lower-dimensional manifold so that we can use the Riemannian gradient descent method. In phase two, we apply a single projected gradient step to update the rank of the underlying variable and escape from spurious local minima arising in the first phase if necessary. To reduce the computation cost of the projected gradient step, we develop pre-processing and warm-start techniques for acceleration. Unlike traditional rank-adaptive methods that require extensive parameter tuning, our hybrid method requires minimal tuning. Extensive experiments confirm the efficiency and robustness of RiNNAL+ in solving various classes of large-scale DNN problems.
In this work, we address the exact D-optimal experimental design problem when the number of design vectors is large. First, we propose a customized column generation algorithm to solve the continuous relaxation of the problem. In the approach, each restricted master problem is constructed carefully so that the number of variables stays small and therefore the subproblem can be solved efficiently by a Primal-Dual Interior-Point-based Semidefinite Programming solver. The support of this solution provides a subset of design points to be used in a local search algorithm for the solution of the integer problem. We prove that a local search algorithm restricted to points of this subset provides an exact design that is provably close to the exact D-optimal design. Our numerical experiments show that, for large-scale instances in which the number of regression points exceeds by far the number of experiments, our approach achieves superior performance compared to existing branch-and-bound-based algorithms in both computational efficiency and solution quality.
This paper proposes an efficient numerical optimization method based on second-order cone programming (SOCP) to solve dynamic optimal transport (DOT) problems with quadratic costs on staggered grid discretizations. By properly reformulating the discretized DOT problem into an equivalent linear SOCP, we develop a highly efficient implementation of an inexact decomposition-based proximal augmented Lagrangian method to solve it. The proposed approach is provided as an open-source software package to facilitate reproducibility and further research. Numerical experiments on a diverse range of DOT problems demonstrate that our software significantly outperforms several state-of-the-art solvers in terms of both accuracy and computational efficiency. Furthermore, it exhibits robust performance when handling measures that are not strictly positive or are in irregular domains with obstacles.
Quadratically constrained quadratic programming (QCQP) has long been recognized as a computationally challenging problem, particularly in large-scale or high-dimensional settings where solving it directly becomes intractable. The complexity further escalates when a sparsity constraint is involved, giving rise to the problem of sparse QCQP (SQCQP), which makes conventional solution methods even less effective. Existing approaches for solving SQCQP typically rely on mixed-integer programming formulations, relaxation techniques, or greedy heuristics but often suffer from computational inefficiency and limited accuracy. In this work, we introduce a novel paradigm by designing an efficient algorithm that directly addresses SQCQP. To be more specific, we introduce P-stationarity to establish first- and second-order optimality conditions of the original problem, leading to a system of nonlinear equations whose generalized Jacobian is proven to be nonsingular under mild assumptions. Most importantly, these equations facilitate the development of a semismooth Newton-type method that exhibits significantly low computational complexity due to the sparsity constraint, achieves a locally quadratic convergence rate, and performs well on large-scale instances in numerical experiments.
Second-order cone programs (SOCPs) with quadratic objective functions are common in optimal control and other fields. Most SOCP solvers which use interior-point methods are designed for linear objectives and convert quadratic objectives into linear ones via slack variables and extra constraints, despite the computational advantages of handling quadratic objectives directly. In applications like model-predictive control and online trajectory optimization, these SOCPs have known sparsity structures and require rapid solutions. When solving these problems, most solvers use sparse linear algebra routines, which introduce computational overhead and hinder performance. In contrast, custom linear algebra routines can exploit the known sparsity structure of problem data and be significantly faster. This work makes two key contributions: (1) the development of QOCO, an open-source C-based solver for quadratic objective SOCPs, and (2) the introduction of QOCOGEN, an open-source custom solver generator for quadratic objective SOCPs, which generates a solver written in C that leverages custom linear algebra. Both implement a primal-dual interior-point method with Mehrotra’s predictor-corrector. On the benchmark problems we run, QOCO is more robust than many commonly used solvers and is faster on small- to medium-sized problems. Additionally, solvers generated by QOCOGEN are significantly faster than QOCO and are free of dynamic memory allocation making them amenable for use on embedded systems.
We present PDLP, a practical first-order method for linear programming (LP) designed to solve large-scale LP problems. PDLP is based on the primal–dual hybrid gradient (PDHG) method applied to the minimax formulation of LP. PDLP incorporates several enhancements to PDHG, including diagonal preconditioning, presolving, adaptive step sizes, adaptive restarting, and feasibility polishing. Our algorithm is implemented in C++, available in Google’s open-source OR-Tools library, and supports multithreading. To evaluate our method, we introduce a new collection of eleven large-scale LP problems with sizes ranging from 125 million to 6.3 billion nonzeros. PDLP solves eight of these instances to optimality gaps of 1 10^-8 ) within 6 days on a single machine. We also compare PDLP with Gurobi barrier, primal simplex, and dual simplex implementations. These traditional methods are designed to solve linear programs to much tighter optimality gaps but struggle to solve these instances. Gurobi barrier solves only three instances, exceeding our 1TB RAM limit on the other eight. While primal and dual simplex are more memory-efficient than the barrier method, they are slower and solve only three instances within 6 days. Compared with the conference version of this work (in: Advances in Neural Information Processing Systems 34 (NeurIPS 2021)), the key new contributions are: (i) feasibility polishing, a technique that quickly finds solutions that are approximately optimal but almost exactly feasible (without which only two of the eleven problems can be solved); (ii) a multithreaded C++ implementation available in Google OR-Tools; and (iii) a new collection of large-scale LP problems. Note that the conference version should be referred to for comparisons with SCS and ablation studies, which we do not repeat in this paper.
We present a graph-theoretic modeling approach for hierarchical optimization that leverages the OptiGraph abstraction implemented in the Julia package Plasmo.jl. We show that the abstraction is flexible and can effectively capture complex hierarchical connectivity that arises from decision-making over multiple spatial and temporal scales (e.g., integration of planning, scheduling, and operations in manufacturing and infrastructures). We also show that the graph abstraction facilitates the conceptualization and implementation of decomposition and approximation schemes. Specifically, we propose a graph-based Benders decomposition (gBD) framework that enables the exploitation of hierarchical (nested) structures and that uses graph aggregation/partitioning procedures to discover such structures. In addition, we provide a Julia implementation of gBD, which we call PlasmoBenders.jl. We illustrate the capabilities using examples arising in the context of energy and power systems.
High-order tensor methods that employ local Taylor models of degree p within adaptive regularization frameworks (ARp) have recently received significant attention, due to their improved/optimal global and local rates of convergence, for both convex and nonconvex optimization problems. The numerical performance of tensor methods for general unconstrained optimization problems remains insufficiently explored/understood, which we address in this paper, by showcasing the numerical performance of standard second- and third-order variants ( p=2,3 ) and proposing novel techniques for key algorithmic aspects when p≥ 3 , to improve the numerical efficiency of tensor variants. In particular, to improve the adaptive choice of the regularization parameter, we extend the interpolation-based updating strategy introduced in [Gould, Porcelli and Toint, Comput Optim Appl (2012) 53:1–22] for p=2 , to the case when p ≥ 3 . We identify fundamental differences between the different local minima of the regularised subproblems for p=2 and p ≥ 3 and their effect on algorithm performance. Then, when p≥ 3 , we introduce a novel pre-rejection technique that rejects poor/unsuccessful subproblem minimizers (that we refer to as ’transient’) prior to any function evaluation, thereby reducing cost and selecting useful (’persistent’) ones. Numerical studies showcase the efficiency improvements generated by our proposed modifications of the AR3 algorithm. We also assess numerically, the effect of different subproblem termination conditions and choice of initial regularization parameter on the overall algorithm performance. Finally, we benchmark our best-performing AR3 variants, as well as those in [Birgin et al, Optim Lett (2020) 14:815–838], against second-order ones (AR2). Encouraging results on standard test problems are obtained, confirming that AR3 variants can be made to outperform second-order variants in terms of objective evaluations, derivative evaluations, and number of subproblem solves. We provide an efficient, extensive and modular software package in MATLAB that includes many AR2 and AR3 variants, including Hessian- and tensor-free ones, allowing ease of use and experimentation for interested users.
Support matrix machine (SMM) is a successful supervised classification model for matrix-type samples. Unlike support vector machines, it employs low-rank regularization on the regression matrix to effectively capture the intrinsic structure embedded in each input matrix. When solving a large-scale SMM, a major challenge arises from the potential increase in sample size, leading to substantial computational and storage burdens. To address these issues, we design a semismooth Newton-CG (SNCG) based augmented Lagrangian method (ALM) for solving the SMM. The ALM exhibits an asymptotic R-superlinear convergence if a strict complementarity condition is satisfied. The SNCG method is employed to solve the ALM subproblems, achieving at least a superlinear convergence rate under the nonemptiness of an index set. Furthermore, the sparsity of samples and the low-rank nature of solutions enable us to reduce the computational cost and storage demands for the Newton linear systems. Additionally, we develop an adaptive sieving strategy that generates a solution path for the SMM by exploiting sample sparsity. The finite convergence of this strategy is also demonstrated. Numerical experiments on both large-scale real and synthetic datasets validate the effectiveness of the proposed methods.
Sample average approximation (SAA) is a technique for obtaining approximate solutions to stochastic programs that uses the average from a random sample to approximate the expected value that is being optimized. Since the outcome from solving an SAA is random, statistical estimates on the optimal value of the true problem can be obtained by solving multiple SAA replications with independent samples. We study techniques to accelerate the solution of this set of SAA replications, when solving them sequentially via Benders decomposition. We investigate how to exploit similarities in the problem structure, as the replications just differ in the realizations of the random samples. Our extensive computational experiments provide empirical evidence that our techniques for using information from solving previous replications can significantly reduce the solution time of later replications.
