Subproblem Solver Research Articles

Fast convergence of trust-regions for non-isolated minima via analysis of CG on indefinite matrices

AbstractTrust-region methods (TR) can converge quadratically to minima where the Hessian is positive definite. However, if the minima are not isolated, then the Hessian there cannot be positive definite. The weaker Polyak–Łojasiewicz (PŁ) condition is compatible with non-isolated minima, and it is enough for many algorithms to preserve good local behavior. Yet, TR with an exact subproblem solver lacks even basic features such as a capture theorem under PŁ. In practice, a popular inexact subproblem solver is the truncated conjugate gradient method (tCG). Empirically, TR-tCG exhibits superlinear convergence under PŁ. We confirm this theoretically. The main mathematical obstacle is that, under PŁ, at points arbitrarily close to minima, the Hessian has vanishingly small, possibly negative eigenvalues. Thus, tCG is applied to ill-conditioned, indefinite systems. Yet, the core theory underlying tCG is that of CG, which assumes a positive definite operator. Accordingly, we develop new tools to analyze the dynamics of CG in the presence of small eigenvalues of any sign, for the regime of interest to TR-tCG.

Fast convergence to non-isolated minima: four equivalent conditions for $${\textrm{C}^{2}}$$ functions

AbstractOptimization algorithms can see their local convergence rates deteriorate when the Hessian at the optimum is singular. These singularities are inescapable when the optima are non-isolated. Yet, under the right circumstances, several algorithms preserve their favorable rates even when optima form a continuum (e.g., due to over-parameterization). This has been explained under various structural assumptions, including the Polyak–Łojasiewicz condition, Quadratic Growth and the Error Bound. We show that, for cost functions which are twice continuously differentiable ($$\textrm{C}^2$$ C 2 ), those three (local) properties are equivalent. Moreover, we show they are equivalent to the Morse–Bott property, that is, local minima form differentiable submanifolds, and the Hessian of the cost function is positive definite along its normal directions. We leverage this insight to improve local convergence guarantees for safe-guarded Newton-type methods under any (hence all) of the above assumptions. First, for adaptive cubic regularization, we secure quadratic convergence even with inexact subproblem solvers. Second, for trust-region methods, we argue capture can fail with an exact subproblem solver, then proceed to show linear convergence with an inexact one (Cauchy steps).

A proximal-gradient method for problems with overlapping group-sparse regularization: support identification complexity

We consider the proximal-gradient method for minimizing the sum of a smooth function and a convex non-smooth overlapping group- ℓ 1 regularizer, which is known to promote sparse solutions with respect to its groups. A feature that distinguishes our work from most in the literature is that the proximal operator for the overlapping group- ℓ 1 function does not admit a closed-form solution, which introduces challenges when proving convergence of the iterates and especially when proving that the iterates correctly identify the group-sparse structure of the optimal solution. To address these challenges, we present an implementable termination condition for the proximal-gradient algorithm, and a specialized primal-dual subproblem solver that is designed to ensure that a group-sparse structure identification property holds. In particular, we give an upper bound on the maximum number of iterations before the correct support (i.e. group-sparse structure) of an optimal solution is identified.

Constructing Three-Phase Envelopes Using a Trust-Region-Based Algorithm

Summary We develop a new trust-region (TR)-based algorithm to construct complete three-phase envelopes for reservoir fluid mixtures. The new algorithm is developed based on a basic algorithm for two-phase envelope constructions (Xu and Li 2023). A state-of-the-art TR method with a realistic exact subproblem solver is implemented in the algorithm, and an integrated strategy is adopted to construct complete phase envelopes with two-phase and three-phase branches. We test the performance of the TR-based algorithm by constructing multiphase envelopes for a total of 15 fluid mixtures that include hydrocarbons-CO2, hydrocarbons-asphaltene, and hydrocarbons-water mixtures. Comparison against the conventional Newton-based algorithm indicates that the TR-based algorithm leads to a much higher computational efficiency and an enhanced robustness.

A trust-region-based phase envelope construction algorithm

Pressure-temperature (PT) and pressure-composition (Px) phase diagrams play an important role in various chemical and petroleum engineering applications. Lots of algorithms based on the Newton-Raphson method have been proposed to numerically construct PT/Px phase envelopes. No attempt is made to improve the computational efficiency of such Newton-Raphson solution method. This work develops a new trust-region-based algorithm to replace the Newton-Raphson method to provide more robust and efficient phase envelope construction calculations. We first convert the phase envelope construction problem into a minimization problem by defining a least-squares objective function. We then apply the trust-region optimization method with an exact subproblem solver to minimize the objective function. We follow the general strategy proposed by Michelsen (1980) to develop the revised phase envelope construction algorithm. We demonstrate the good performance of the trust-region-based algorithm by comparing it with the conventional Newton's method. The computational results indicate that the new trust-region-based algorithm is more robust and cost-effective than the Newton-based algorithm.

Inexact penalty decomposition methods for optimization problems with geometric constraints

This paper provides a theoretical and numerical investigation of a penalty decomposition scheme for the solution of optimization problems with geometric constraints. In particular, we consider some situations where parts of the constraints are nonconvex and complicated, like cardinality constraints, disjunctive programs, or matrix problems involving rank constraints. By a variable duplication and decomposition strategy, the method presented here explicitly handles these difficult constraints, thus generating iterates which are feasible with respect to them, while the remaining (standard and supposingly simple) constraints are tackled by sequential penalization. Inexact optimization steps are proven sufficient for the resulting algorithm to work, so that it is employable even with difficult objective functions. The current work is therefore a significant generalization of existing papers on penalty decomposition methods. On the other hand, it is related to some recent publications which use an augmented Lagrangian idea to solve optimization problems with geometric constraints. Compared to these methods, the decomposition idea is shown to be numerically superior since it allows much more freedom in the choice of the subproblem solver, and since the number of certain (possibly expensive) projection steps is significantly less. Extensive numerical results on several highly complicated classes of optimization problems in vector and matrix spaces indicate that the current method is indeed very efficient to solve these problems.

MINRES: From Negative Curvature Detection to Monotonicity Properties

The conjugate gradient method (CG) has long been the workhorse for inner-iterations of second-order algorithms for large-scale nonconvex optimization. Prominent examples include line-search based algorithms, e.g., Newton-CG, and those based on a trust-region framework, e.g., CG-Steihaug. This is mainly thanks to CG's several favorable properties, including certain monotonicity properties and its inherent ability to detect negative curvature directions, which can arise in nonconvex optimization. This is despite the fact that the iterative method-of-choice when it comes to real symmetric but potentially indefinite matrices is arguably the celebrated minimal residual (MINRES) method. However, limited understanding of similar properties implied by MINRES in such settings has restricted its applicability within nonconvex optimization algorithms. We establish several such nontrivial properties of MINRES, including certain useful monotonicity as well as an inherent ability to detect negative curvature directions. These properties allow MINRES to be considered as a potentially superior alternative to CG for all Newton-type nonconvex optimization algorithms that employ CG as their subproblem solver.

Newton-MR: Inexact Newton Method with minimum residual sub-problem solver

We consider a variant of inexact Newton Method [20] , [40] , called Newton-MR, in which the least-squares sub-problems are solved approximately using Minimum Residual method [79] . By construction, Newton-MR can be readily applied for unconstrained optimization of a class of non-convex problems known as invex, which subsumes convexity as a sub-class. For invex optimization, instead of the classical Lipschitz continuity assumptions on gradient and Hessian, Newton-MR's global convergence can be guaranteed under a weaker notion of joint regularity of Hessian and gradient. We also obtain Newton-MR's problem-independent local convergence to the set of minima. We show that fast local/global convergence can be guaranteed under a novel inexactness condition, which, to our knowledge, is much weaker than the prior related works. Numerical results demonstrate the performance of Newton-MR as compared with several other Newton-type alternatives on a few machine learning problems.

Variable-Metric Localization of Occupied and Virtual Orbitals.

The key idea of the variable-metric approach to orbital localization is to allow nonorthogonality between orbitals while, at the same time, preventing them from becoming linearly dependent. The variable-metric localization has been shown to improve the locality of occupied nonorthogonal orbitals relative to their orthogonal counterparts. In this work, numerous localization algorithms are designed and tested to exploit the conceptual simplicity of the variable-metric approach with the goal of creating a straightforward and reliable localization procedure for virtual orbitals. The implemented algorithms include the steepest descent, conjugate gradient (CG), limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS), and hybrid procedures as well as trust-region (TR) methods based on the CG and Cauchy-point subproblem solvers. Comparative analysis shows that the CG-based TR algorithm is the best overall method to obtain nonorthogonal localized molecular orbitals (NLMOs), occupied or virtual. The L-BFGS and CG algorithms can also be used to obtain NLMOs reliably but often at higher computational cost. Extensive tests demonstrate that the implemented methods allow us to obtain well-localized Boys-Foster (i.e., maximally localized Wannier functions) and Pipek-Mezey, orthogonal and nonorthogonal, and occupied and virtual orbitals for a variety of gas-phase molecules and periodic materials. The tests also show that virtual NLMOs, which have not been described before, are, on average, 13% (Boys-Foster) and 18% (Pipek-Mezey) more localized than their orthogonal counterparts.

Performance Analysis of Trust Region Subproblem Solvers for Limited-Memory Distributed BFGS Optimization Method

The limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) optimization method performs very efficiently for large-scale problems. A trust region search method generally performs more efficiently and robustly than a line search method, especially when the gradient of the objective function cannot be accurately evaluated. The computational cost of an L-BFGS trust region subproblem (TRS) solver depend mainly on the number of unknown variables (n) and the number of variable shift vectors and gradient change vectors (m) used for Hessian updating, withm << nfor large-scale problems. In this paper, we analyze the performances of different methods to solve the L-BFGS TRS. The first method is the direct method using the Newton-Raphson (DNR) method and Cholesky factorization of a densen×nmatrix, the second one is the direct method based on an inverse quadratic (DIQ) interpolation, and the third one is a new method that combines the matrix inversion lemma (MIL) with an approach to update associated matrices and vectors. The MIL approach is applied to reduce the dimension of the original problem withnvariables to a new problem withmvariables. Instead of directly using expensive matrix-matrix and matrix-vector multiplications to solve the L-BFGS TRS, a more efficient approach is employed to update matrices and vectors iteratively. The L-BFGS TRS solver using the MIL method performs more efficiently than using the DNR method or DIQ method. Testing on a representative suite of problems indicates that the new method can converge to optimal solutions comparable to those obtained using the DNR or DIQ method. Its computational cost represents only a modest overhead over the well-known L-BFGS line-search method but delivers improved stability in the presence of inaccurate gradients. When compared to the solver using the DNR or DIQ method, the new TRS solver can reduce computational cost by a factor proportional ton2/mfor large-scale problems.

An accelerated first-order method with complexity analysis for solving cubic regularization subproblems

We propose a first-order method to solve the cubic regularization subproblem (CRS) based on a novel reformulation. The reformulation is a constrained convex optimization problem whose feasible region admits an easily computable projection. Our reformulation requires computing the minimum eigenvalue of the Hessian. To avoid the expensive computation of the exact minimum eigenvalue, we develop a surrogate problem to the reformulation where the exact minimum eigenvalue is replaced with an approximate one. We then apply first-order methods such as the Nesterov’s accelerated projected gradient method (APG) and projected Barzilai-Borwein method to solve the surrogate problem. As our main theoretical contribution, we show that when an $$\epsilon$$ -approximate minimum eigenvalue is computed by the Lanczos method and the surrogate problem is approximately solved by APG, our approach returns an $$\epsilon$$ -approximate solution to CRS in $${\tilde{O}}(\epsilon ^{-1/2})$$ matrix-vector multiplications (where $${\tilde{O}}(\cdot )$$ hides the logarithmic factors). Numerical experiments show that our methods are comparable to and outperform the Krylov subspace method in the easy and hard cases, respectively. We further implement our methods as subproblem solvers of adaptive cubic regularization methods, and numerical results show that our algorithms are comparable to the state-of-the-art algorithms.

Simultaneous iterative solutions for the trust-region and minimum eigenvalue subproblem

ABSTRACT Given the inability to foresee all possible scenarios, it is justified to desire an efficient trust-region subproblem solver capable of delivering any desired level of accuracy on demand; that is, the accuracy obtainable for a given trust-region subproblem should not be partially dependent on the problem itself. Current state-of-the-art iterative eigensolvers all fall into the class of restarted Lanczos methods; whereas, current iterative trust-region solvers at best reduce to unrestarted Lanczos methods; which in this context are well known to be numerically unstable with impractical memory requirements. In this paper, we present the first iterative trust region subproblem solver that at its core contains a robust and practical eigensolver. Our solver leverages the recently announced work of Stathopoulos and Orginos which has not been noticed by the optimization community and can be utilized because, unlike other restarted Lanczos methods, its restarts do not necessarily modify the current Lanczos sequence generated by Conjugate Gradient methods (CG). This innovated strategy can be utilized in the context of TR solvers as well. Moreover, our TR subproblem solver adds negligible computational overhead compared to existing iterative TR approaches.

TASTE: Temporal and Static Tensor Factorization for Phenotyping Electronic Health Records.

Phenotyping electronic health records (EHR) focuses on defining meaningful patient groups (e.g., heart failure group and diabetes group) and identifying the temporal evolution of patients in those groups. Tensor factorization has been an effective tool for phenotyping. Most of the existing works assume either a static patient representation with aggregate data or only model temporal data. However, real EHR data contain both temporal (e.g., longitudinal clinical visits) and static information (e.g., patient demographics), which are difficult to model simultaneously. In this paper, we propose Temporal And Static TEnsor factorization (TASTE) that jointly models both static and temporal information to extract phenotypes. TASTE combines the PARAFAC2 model with non-negative matrix factorization to model a temporal and a static tensor. To fit the proposed model, we transform the original problem into simpler ones which are optimally solved in an alternating fashion. For each of the sub-problems, our proposed mathematical re-formulations lead to efficient sub-problem solvers. Comprehensive experiments on large EHR data from a heart failure (HF) study confirmed that TASTE is up to 14× faster than several baselines and the resulting phenotypes were confirmed to be clinically meaningful by a cardiologist. Using 60 phenotypes extracted by TASTE, a simple logistic regression can achieve the same level of area under the curve (AUC) for HF prediction compared to a deep learning model using recurrent neural networks (RNN) with 345 features.

Outer approximation with conic certificates for mixed-integer convex problems

A mixed-integer convex (MI-convex) optimization problem is one that becomes convex when all integrality constraints are relaxed. We present a branch-and-bound LP outer approximation algorithm for an MI-convex problem transformed to MI-conic form. The polyhedral relaxations are refined with $\mathcal{K}^*$ cuts derived from conic certificates for continuous primal-dual conic subproblems. Under the assumption that all subproblems are well-posed, the algorithm detects infeasibility or unboundedness or returns an optimal solution in finite time. Using properties of the conic certificates, we show that the $\mathcal{K}^*$ cuts imply certain practically-relevant guarantees about the quality of the polyhedral relaxations, and demonstrate how to maintain helpful guarantees when the LP solver uses a positive feasibility tolerance. We discuss how to disaggregate $\mathcal{K}^*$ cuts in order to tighten the polyhedral relaxations and thereby improve the speed of convergence, and propose fast heuristic methods of obtaining useful $\mathcal{K}^*$ cuts. Our new open source MI-conic solver Pajarito (http://github.com/JuliaOpt/Pajarito.jl) uses an external mixed-integer linear (MILP) solver to manage the search tree and an external continuous conic solver for subproblems. Benchmarking on a library of mixed-integer second-order cone (MISOCP) problems, we find that Pajarito greatly outperforms Bonmin (the leading open source alternative) and is competitive with CPLEX's specialized MISOCP algorithm. We demonstrate the robustness of Pajarito by solving diverse MI-conic problems involving mixtures of positive semidefinite, second-order, and exponential cones, and provide evidence for the practical value of our analyses and enhancements of $\mathcal{K}^*$ cuts.

Stochastic trust region inexact Newton method for large-scale machine learning

Nowadays stochastic approximation methods are one of the major research direction to deal with the large-scale machine learning problems. From stochastic first order methods, now the focus is shifting to stochastic second order methods due to their faster convergence and availability of computing resources. In this paper, we have proposed a novel Stochastic Trust RegiOn Inexact Newton method, called as STRON, to solve large-scale learning problems which uses conjugate gradient (CG) to inexactly solve trust region subproblem. The method uses progressive subsampling in the calculation of gradient and Hessian values to take the advantage of both, stochastic and full-batch regimes. We have extended STRON using existing variance reduction techniques to deal with the noisy gradients and using preconditioned conjugate gradient (PCG) as subproblem solver, and empirically proved that they do not work as expected, for the large-scale learning problems. Finally, our empirical results prove efficacy of the proposed method against existing methods with bench marked datasets.

A Direct Algorithm for Optimization Problems With the Huber Penalty.

We present a direct (noniterative) algorithm for 1-D quadratic data fitting with neighboring intensity differences penalized by the Huber function. Applications of such an algorithm include 1-D processing of medical signals, such as smoothing of tissue time concentration curves in kinetic data analysis or sinogram preprocessing, and using it as a subproblem solver for 2-D or 3-D image restoration and reconstruction. dynamic programming was used to develop the direct algorithm. The problem was reformulated as a sequence of univariate optimization problems, for , where is the number of data points. The solution to the univariate problem at index is parameterized by the solution at , except at . Solving the univariate optimization problem at yields the solution to each problem in the sequence using back-tracking. Computational issues and memory cost are discussed in detail. Two numerical studies, tissue concentration curve smoothing and sinogram preprocessing for image reconstruction, are used to validate the direct algorithm and illustrate its practical applications. In the example of 1-D curve smoothing, the efficiency of the direct algorithm is compared with four iterative methods: the iterative coordinate descent, Nesterov's accelerated gradient descent algorithm, FISTA, and an off-the-shelf second order method. The first two methods were applied to the primal problem, the others to the dual problem. The comparisons show that the direct algorithm outperforms all other methods by a significant factor, which rapidly grows with the curvature of the Huber function. The second example, sinogram preprocessing, showed that robustness and speed of the direct algorithm are maintained over a wide range of signal variations, and that noise and streaking artifacts could be reduced with almost no increase in computation time. We also outline how the proposed 1-D solver can be used for imaging applications.

Numerical studies of a class of linear solvers for fine-scale petroleum reservoir simulation

Numerical simulation based on fine-scale reservoir models helps petroleum engineers in understanding fluid flow in porous media and achieving higher recovery ratio. Fine-scale models give rise to large-scale linear systems, and thus require effective solvers for solving these linear systems to finish simulation in reasonable turn-around time. In this paper, we study convergence, robustness, and efficiency of a class of multi-stage preconditioners accelerated by Krylov subspace methods for solving Jacobian systems from a fully implicit discretization. We compare components of these preconditioners, including decoupling and sub-problem solvers, for fine-scale reservoir simulation. Several benchmark and real-world problems, including a ten-million-cell reservoir problem, were simulated on a desktop computer. Numerical tests show that the combination of the alternating block factorization method and multi-stage subspace correction preconditioner gives a robust and memory-efficient solver for fine-scale reservoir simulation.

A branch and price approach for deployment of multi-tier software services in clouds

This paper considers a service deployment problem that combines service placement and replication level decisions in a cloud computing context. The services are composed of multiple components that are to be placed on nodes in the private cloud of the service provider or, if the private cloud has limited capacity, partly in a public cloud. In the service delivery, the provider has to take into account the quality of service guarantees offered to his end-users. To solve the problem, we develop a branch and price algorithm, where the subproblems both are formulated as a linear mixed integer program and a shortest path problem with resource constraints (SPPRC) on a network with a special structure. The SPPRC can be solved by an exact label-setting algorithm, but to speed up the solution process, we develop a heuristic label-setting algorithm based on a reduced network and simplified dominance rule. Our results show that using the heuristic subproblem solver is efficient. Furthermore, the branch and price algorithm performs better than a previously developed pre-generation algorithm for the same problem. In addition, we analyze and discuss the differences in solutions that utilize resources in a public cloud to different degrees. By conducting this analysis we are able to identify some essential characteristics of good solutions.

Column generation based approaches for a tour scheduling problem with a multi-skill heterogeneous workforce

In this paper, we address a multi-activity tour scheduling problem with time varying demand. The objective is to compute a team schedule for a fixed roster of employees in order to minimize the over-coverage and the under-coverage of different parallel activity demands along a planning horizon of one week. Numerous complicating constraints are present in our problem: all employees are different and can perform several different activities during the same day-shift, lunch breaks and pauses are flexible, demand is given for 15 minutes periods. Employees have feasibility and legality rules to be satisfied, but the objective function does not account for any quality measure associated with each individual’s schedule. More precisely, the problem mixes simultaneously days-off scheduling, shift scheduling, shift assignment, activity assignment, pause and lunch break assignment.To solve this problem, we developed four methods: a compact Mixed Integer Linear Programming model, a branch-and-price like approach with a nested dynamic program to solve heuristically the subproblems, a diving heuristic and a greedy heuristic based on our subproblem solver. The computational results, based on both real cases and instances derived from real cases, demonstrate that our methods are able to provide good quality solutions in a short computing time. Our algorithms are now embedded in a commercial software, which is already in use in a mini-mart company.

Shifting bottleneck scheduling for total weighted tardiness minimization—A computational evaluation of subproblem and re-optimization heuristics

Machine-based decomposition of total weighted tardiness job shops is known to be considerably more complicated than in the makespan case, mainly due to the structure of the underlying graph model and thus the arising one-machine subproblems. In fact, the effectiveness of a shifting bottleneck approach crucially depends on the employed subproblem solver. Although a sophisticated exact algorithm exists, problem instances involving more than 30 jobs are still challenging. In this paper, new heuristic approaches to subproblems of this kind are devised which rely on advanced problem-specific concepts like local optimality and dominance principles. The proposed subproblem solvers are combined with an iterated local search method for re-optimizing already scheduled machines. Computational experiments show that the final enhanced shifting bottleneck algorithms are not only applicable to job shops involving up to 100 jobs and 20 machines, but also able to improve existing results for benchmark instances.