Search

Convex nonsmooth optimization problems, whose solutions live in very high dimensional spaces, have become ubiquitous. To solve them, the class of first-order algorithms known as proximal splitting algorithms is particularly adequate: they consist of simple operations, handling the terms in the objective function separately. In this overview, we demystify a selection of recent proximal splitting algorithms: we present them within a unified framework, which consists in applying splitting methods for monotone inclusions in primal-dual product spaces, with well-chosen metrics. Along the way, we easily derive new variants of the algorithms and revisit existing convergence results, extending the parameter ranges in several cases. In particular, we emphasize that when the smooth term in the objective function is quadratic, e.g., for least-squares problems, convergence is guaranteed with larger values of the relaxation parameter than previously known. Such larger values are usually beneficial for the convergence speed in practice., To appear in SIAM Review

We present the concept of a Generalized Feedback Nash Equilibrium (GFNE) in dynamic games, extending the Feedback Nash Equilibrium concept to games in which players are subject to state and input constraints. We formalize necessary and sufficient conditions for (local) GFNE solutions at the trajectory level, which enable the development of efficient numerical methods for their computation. Specifically, we propose a Newton-style method for finding game trajectories which satisfy necessary conditions for an equilibrium, which can then be checked against sufficiency conditions. We show that the evaluation of the necessary conditions in general requires computing a series of nested, implicitly-defined derivatives, which quickly becomes intractable. To this end, we introduce an approximation to the necessary conditions which is amenable to efficient evaluation, and in turn, computation of solutions. We term the solutions to the approximate necessary conditions Generalized Feedback Quasi-Nash Equilibria (GFQNE), and we introduce numerical methods for their computation. In particular, we develop a Sequential Linear-Quadratic Game approach, in which a LQ local approximation of the game is solved at each iteration. The development of this method relies on the ability to compute a GFNE to inequality- and equality-constrained LQ games, and therefore specific methods for the solution of these special cases are developed in detail. We demonstrate the effectiveness of the proposed solution approach on a dynamic game arising in an autonomous driving application.

We study the problem of computing weighted sum-of-squares (WSOS) certificates for positive polynomials over a compact semialgebraic set. Building on the theory of interior-point methods for convex optimization, we introduce the concept of dual cone certificates, which allows us to interpret vectors from the dual of the sum-of-squares cone as rigorous nonnegativity certificates of a WSOS polynomial. Whereas conventional WSOS certificates are alternative representations of the polynomials they certify, dual certificates are distinct from the certified polynomials; moreover, each dual certificate certifies a full-dimensional convex cone of WSOS polynomials. As a result, rational WSOS certificates can be constructed from numerically computed dual certificates at little additional cost, without any rounding or projection steps applied to the numerical certificates. As an additional algorithmic application, we present an almost entirely numerical hybrid algorithm for computing the optimal WSOS lower bound of a given polynomial along with a rational dual certificate, with a polynomial-time computational cost per iteration and linear rate of convergence., Comment: Preprint accepted in SIAM Journal on Optimization. Major revision compared to v1

We study a general convex optimization problem, which covers various classic problems in different areas and particularly includes many optimal transport related problems arising in recent years. To solve this problem, we revisit the classic Bregman proximal point algorithm (BPPA) and introduce a new inexact stopping condition for solving the subproblems, which can circumvent the underlying feasibility difficulty often appearing in existing inexact conditions when the problem has a complex feasible set. Our inexact condition also covers several existing inexact conditions as special cases and hence makes our inexact BPPA (iBPPA) more flexible to fit different scenarios in practice. Moreover, inspired by Nesterov's acceleration technique, we develop an inertial variant of our iBPPA, denoted by V-iBPPA, and establish the iteration complexity of $O(1/k^{\lambda})$, where $\lambda\geq1$ is a quadrangle scaling exponent of the kernel function. In particular, when the proximal parameter is a constant and the kernel function is strongly convex with Lipschitz continuous gradient (hence $\lambda=2$), our V-iBPPA achieves a faster rate of $O(1/k^2)$ just as existing accelerated inexact proximal point algorithms. Some preliminary numerical experiments for solving the standard OT problem are conducted to show the convergence behaviors of our iBPPA and V-iBPPA under different inexactness settings. The experiments also empirically verify the potential of our V-iBPPA on improving the convergence speed.

We investigate the feasibility of sample average approximation (SAA) for general stochastic optimization problems, including two-stage stochastic programming without the relatively complete recourse assumption. Instead of analyzing problems with specific structures, we utilize results from the Vapnik-Chervonenkis (VC) dimension and Probably Approximately Correct learning to provide a general framework that offers explicit feasibility bounds for SAA solutions under minimal structural or distributional assumption. We show that, as long as the hypothesis class formed by the feasbible region has a finite VC dimension, the infeasibility of SAA solutions decreases exponentially with computable rates and explicitly identifiable accompanying constants. We demonstrate how our bounds apply more generally and competitively compared to existing results.

We provide a proximal average with repect to a $1$-coercive Legendre function. In the sense of Bregman distance, the Bregman envelope of the proximal average is a convex combination of Bregman envelopes of individual functions. The Bregman proximal mapping of the average is a convex combination of convexified proximal mappings of individual functions. Techniques from variational analysis provide the keys for the Bregman proximal average.

We study solution sensitivity for nonlinear programs (NLPs) whose structures are induced by graphs. These NLPs arise in many applications such as dynamic optimization, stochastic optimization, optimization with partial differential equations, and network optimization. We show that for a given pair of nodes, the sensitivity of the primal-dual solution at one node against a data perturbation at the other node decays exponentially with respect to the distance between these two nodes on the graph. In other words, the solution sensitivity decays as one moves away from the perturbation point. This result, which we call exponential decay of sensitivity, holds under the strong second-order sufficiency condition and the linear independence constraint qualification. We also present conditions under which the decay rate remains uniformly bounded; this allows us to characterize the sensitivity behavior of NLPs defined over subgraphs of infinite graphs. The theoretical developments are illustrated with numerical examples.

The focus of this paper is on stochastic variational inequalities (VI) under Markovian noise. A prominent application of our algorithmic developments is the stochastic policy evaluation problem in reinforcement learning. Prior investigations in the literature focused on temporal difference (TD) learning by employing nonsmooth finite time analysis motivated by stochastic subgradient descent leading to certain limitations. These encompass the requirement of analyzing a modified TD algorithm that involves projection to an a-priori defined Euclidean ball, achieving a non-optimal convergence rate and no clear way of deriving the beneficial effects of parallel implementation. Our approach remedies these shortcomings in the broader context of stochastic VIs and in particular when it comes to stochastic policy evaluation. We developed a variety of simple TD learning type algorithms motivated by its original version that maintain its simplicity, while offering distinct advantages from a non-asymptotic analysis point of view. We first provide an improved analysis of the standard TD algorithm that can benefit from parallel implementation. Then we present versions of a conditional TD algorithm (CTD), that involves periodic updates of the stochastic iterates, which reduce the bias and therefore exhibit improved iteration complexity. This brings us to the fast TD (FTD) algorithm which combines elements of CTD and the stochastic operator extrapolation method of the companion paper. For a novel index resetting policy FTD exhibits the best known convergence rate. We also devised a robust version of the algorithm that is particularly suitable for discounting factors close to 1., Comment: arXiv admin note: text overlap with arXiv:2011.02987

This paper is concerned with the column $\ell_{2,0}$-regularized factorization model of low-rank matrix recovery problems and its computation. The column $\ell_{2,0}$-norm of factor matrices is introduced to promote column sparsity of factors and low-rank solutions. For this nonconvex discontinuous optimization problem, we develop an alternating majorization-minimization (AMM) method with extrapolation, and a hybrid AMM in which a majorized alternating proximal method is proposed to seek an initial factor pair with less nonzero columns and the AMM with extrapolation is then employed to minimize of a smooth nonconvex loss. We provide the global convergence analysis for the proposed AMM methods and apply them to the matrix completion problem with non-uniform sampling schemes. Numerical experiments are conducted with synthetic and real data examples, and comparison results with the nuclear-norm regularized factorization model and the max-norm regularized convex model show that the column $\ell_{2,0}$-regularized factorization model has an advantage in offering solutions of lower error and rank within less time.

Chance constraints provide a principled framework to mitigate the risk of high-impact extreme events by modifying the controllable properties of a system. The low probability and rare occurrence of such events, however, impose severe sampling and computational requirements on classical solution methods that render them impractical. This work proposes a novel sampling-free method for solving rare chance constrained optimization problems affected by uncertainties that follow general Gaussian mixture distributions. By integrating modern developments in large deviation theory with tools from convex analysis and bilevel optimization, we propose tractable formulations that can be solved by off-the-shelf solvers. Our formulations enjoy several advantages compared to classical methods: their size and complexity is independent of event rarity, they do not require linearity or convexity assumptions on system constraints, and under easily verifiable conditions, serve as safe conservative approximations or asymptotically exact reformulations of the true problem. Computational experiments on linear, nonlinear and PDE-constrained problems from applications in portfolio management, structural engineering and fluid dynamics illustrate the broad applicability of our method and its advantages over classical sampling-based approaches in terms of both accuracy and efficiency.

Random projection techniques based on Johnson-Lindenstrauss lemma are used for randomly aggregating the constraints or variables of optimization problems while approximately preserving their optimal values, that leads to smaller-scale optimization problems. D'Ambrosio et al. have applied random projection to a quadratic optimization problem so as to decrease the number of decision variables. Although the problem size becomes smaller, the projected problem will also almost surely be non-convex if the original problem is non-convex, and hence will be hard to solve. In this paper, by focusing on the fact that the level of the non-convexity of a non-convex quadratic optimization problem can be alleviated by random projection, we find an approximate global optimal value of the problem by attributing it to a convex problem with smaller size. To the best of our knowledge, our paper is the first to use random projection for convexification of non-convex optimization problems. We evaluate the approximation error between optimum values of a non-convex optimization problem and its convexified randomly projected problem., 25 pages, 4 figures

We consider optimization problems on Riemannian manifolds with equality and inequality constraints, which we call Riemannian nonlinear optimization (RNLO) problems. Although they have numerous applications, the existing studies on them are limited especially in terms of algorithms. In this paper, we propose Riemannian sequential quadratic optimization (RSQO) that uses a line-search technique with an ell_1 penalty function as an extension of the standard SQO algorithm for constrained nonlinear optimization problems in Euclidean spaces to Riemannian manifolds. We prove its global convergence to a Karush-Kuhn-Tucker point of the RNLO problem by means of parallel transport and the exponential mapping. Furthermore, we establish its local quadratic convergence by analyzing the relationship between sequences generated by RSQO and the Riemannian Newton method. Ours is the first algorithm that has both global and local convergence properties for constrained nonlinear optimization on Riemannian manifolds. Empirical results show that RSQO finds solutions more stably and with higher accuracy compared with the existing Riemannian penalty and augmented Lagrangian methods., Comment: 36 pages, 2 figure

We consider the bilinear optimal control of an advection-reaction-diffusion system, where the control arises as the velocity field in the advection term. Such a problem is generally challenging from both theoretical analysis and algorithmic design perspectives mainly because the state variable depends nonlinearly on the control variable and an additional divergence-free constraint on the control is coupled together with the state equation. Mathematically, the proof of the existence of optimal solutions is delicate, and up to now, only some results are known for a few special cases where additional restrictions are imposed on the space dimension and the regularity of the control. We prove the existence of optimal controls and derive the first-order optimality conditions in general settings without any extra assumption. Computationally, the well-known conjugate gradient (CG) method can be applied conceptually. However, due to the additional divergence-free constraint on the control variable and the nonlinear relation between the state and control variables, it is challenging to compute the gradient and the optimal stepsize at each CG iteration, and thus nontrivial to implement the CG method. To address these issues, we advocate a fast inner preconditioned CG method to ensure the divergence-free constraint and an efficient inexactness strategy to determine an appropriate stepsize. An easily implementable nested CG method is thus proposed for solving such a complicated problem. For the numerical discretization, we combine finite difference methods for the time discretization and finite element methods for the space discretization. Efficiency of the proposed nested CG method is promisingly validated by the results of some preliminary numerical experiments.

We consider the problem of minimizing a high-dimensional objective function, which may include a regularization term, using (possibly noisy) evaluations of the function. Such optimization is also called derivative-free, zeroth-order, or black-box optimization. We propose a new $\textbf{Z}$eroth-$\textbf{O}$rder $\textbf{R}$egularized $\textbf{O}$ptimization method, dubbed ZORO. When the underlying gradient is approximately sparse at an iterate, ZORO needs very few objective function evaluations to obtain a new iterate that decreases the objective function. We achieve this with an adaptive, randomized gradient estimator, followed by an inexact proximal-gradient scheme. Under a novel approximately sparse gradient assumption and various different convex settings, we show the (theoretical and empirical) convergence rate of ZORO is only logarithmically dependent on the problem dimension. Numerical experiments show that ZORO outperforms the existing methods with similar assumptions, on both synthetic and real datasets.

We consider the problem of minimizing an objective function that is the sum of a convex function and a group sparsity-inducing regularizer. Problems that integrate such regularizers arise in modern machine learning applications, often for the purpose of obtaining models that are easier to interpret and that have higher predictive accuracy. We present a new method for solving such problems that utilize subspace acceleration, domain decomposition, and support identification. Our analysis shows, under common assumptions, that the iterate sequence generated by our framework is globally convergent, converges to an $\epsilon$-approximate solution in at most $O(\epsilon^{-(1+p)})$ (respectively, $O(\epsilon^{-(2+p)})$) iterations for all $\epsilon$ bounded above and large enough (respectively, all $\epsilon$ bounded above) where $p > 0$ is an algorithm parameter, and exhibits superlinear local convergence. Preliminary numerical results for the task of binary classification based on regularized logistic regression show that our approach is efficient and robust, with the ability to outperform a state-of-the-art method.

Auxiliary variables are often used to model a convex piecewise linear function in the framework of linear optimization. This work shows that such variables yield a block diagonal plus low rank structure in the reduced KKT system of the dual problem. We show how the structure can be detected efficiently, and derive the linear algebra formulas for an interior point method which exploit such structure. The structure is detected in 36% of the cases in Netlib. Numerical results on the inverse planning problem in radiation therapy show an order of magnitude speed-up compared to the state-of-the-art interior point solver CPLEX, and considerable improvements in dose distribution compared to current algorithms.

Linear complementarity problems are a powerful tool for modeling many practically relevant situations such as market equilibria. They also connect many sub-areas of mathematics like game theory, optimization, and matrix theory. Despite their close relation to optimization, the protection of LCPs against uncertainties -- especially in the sense of robust optimization -- is still in its infancy. During the last years, robust LCPs have only been studied using the notions of strict and $\Gamma$-robustness. Unfortunately, both concepts lead to the problem that the existence of robust solutions cannot be guaranteed. In this paper, we consider affinely adjustable robust LCPs. In the latter, a part of the LCP solution is allowed to adjust via a function that is affine in the uncertainty. We show that this notion of robustness allows to establish strong characterizations of solutions for the cases of uncertain matrix and vector, separately, from which existence results can be derived. Our main results are valid for the case of an uncertain LCP vector. Here, we additionally provide sufficient conditions on the LCP matrix for the uniqueness of a solution. Moreover, based on characterizations of the affinely adjustable robust solutions, we derive a mixed-integer programming formulation that allows to solve the corresponding robust counterpart. If, in addition, the certain LCP matrix is positive semidefinite, we prove polynomial-time solvability and uniqueness of robust solutions. If the LCP matrix is uncertain, characterizations of solutions are developed for every nominal matrix, i.e., these characterizations are, in particular, independent of the definiteness of the nominal matrix. Robust solutions are also shown to be unique for positive definite LCP matrix but both uniqueness and mixed-integer programming formulations still remain open problems if the nominal LCP matrix is not psd., Comment: 20 pages

This paper describes an extension of the BFGS and L-BFGS methods for the minimization of a nonlinear function subject to errors. This work is motivated by applications that contain computational noise, employ low-precision arithmetic, or are subject to statistical noise. The classical BFGS and L-BFGS methods can fail in such circumstances because the updating procedure can be corrupted and the line search can behave erratically. The proposed method addresses these difficulties and ensures that the BFGS update is stable by employing a lengthening procedure that spaces out the points at which gradient differences are collected. A new line search, designed to tolerate errors, guarantees that the Armijo-Wolfe conditions are satisfied under most reasonable conditions, and works in conjunction with the lengthening procedure. The proposed methods are shown to enjoy convergence guarantees for strongly convex functions. Detailed implementations of the methods are presented, together with encouraging numerical results., 27 pages, 13 figures, 2 tables

This work is a follow-up and a complement to arXiv:1912.08899 [math.OC] for solving polynomial optimization problems (POPs). The chordal-TSSOS hierarchy that we propose is a new sparse moment-SOS framework based on term-sparsity and chordal extension. By exploiting term-sparsity of the input polynomials we obtain a two-level hierarchy of semidefinite programming relaxations. The novelty and distinguishing feature of such relaxations is to obtain quasi block-diagonal matrices obtained in an iterative procedure that performs chordal extension of certain adjacency graphs. The graphs are related to the terms arising in the original data and not to the links between variables. Various numerical examples demonstrate the efficiency and the scalability of this new hierarchy for both unconstrained and constrained POPs. The two hierarchies are complementary. While the former TSSOS arXiv:1912.08899 [math.OC] has a theoretical convergence guarantee, the chordal-TSSOS has superior performance but lacks this theoretical guarantee., 29 pages, 14 tables, 5 figures. arXiv admin note: text overlap with arXiv:1912.08899

In this paper, we develop convergence analysis of a modified line search method for objective functions whose value is computed with noise and whose gradient estimates are inexact and possibly random. The noise is assumed to be bounded in absolute value without any additional assumptions. We extend the framework based on stochastic methods from [Cartis and Scheinberg, 2018] which was developed to provide analysis of a standard line search method with exact function values and random gradients to the case of noisy functions. We introduce two alternative conditions on the gradient which when satisfied with some sufficiently large probability at each iteration, guarantees convergence properties of the line search method. We derive expected complexity bounds to reach a near optimal neighborhood for convex, strongly convex and nonconvex functions. The exact dependence of the convergence neighborhood on the noise is specified., Comment: 30 pages. arXiv admin note: text overlap with arXiv:1905.01332

In this paper we give an algorithm to round the floating point output of a semidefinite programming solver to a solution over the rationals or a quadratic extension of the rationals. We apply this to get sharp bounds for packing problems, and we use these sharp bounds to prove that certain optimal packing configurations are unique up to rotations. In particular, we show that the configuration coming from the $\mathsf{E}_8$ root lattice is the unique optimal code with minimal angular distance $\pi/3$ on the hemisphere in $\mathbb R^8$, and we prove that the three-point bound for the $(3, 8, \vartheta)$-spherical code, where $\vartheta$ is such that $\cos \vartheta = (2\sqrt{2}-1)/7$, is sharp by rounding to $\mathbb Q[\sqrt{2}]$. We also use our machinery to compute sharp upper bounds on the number of spheres that can be packed into a larger sphere., Comment: 24 pages

This paper presents an accelerated composite gradient (ACG) variant, referred to as the AC-ACG method, for solving nonconvex smooth composite minimization problems. As opposed to well-known ACG variants that are either based on a known Lipschitz gradient constant or a sequence of maximum observed curvatures, the current one is based on the average of all past observed curvatures. More specifically, AC-ACG uses a positive multiple of the average of all observed curvatures until the previous iteration as a way to estimate the "function curvature" at the current point and then two resolvent evaluations to compute the next iterate. In contrast to other variable Lipschitz estimation variants, e.g., the ones based on the maximum curvature, AC-ACG always accepts the aforementioned iterate regardless how poor the Lipschitz estimation turns out to be. Finally, computational results are presented to illustrate the efficiency of AC-ACG on both randomly generated and real-world problem instances., 28 pages

We introduce BELLA, a locally superlinearly convergent Bregman forward-backward splitting method for minimizing the sum of two nonconvex functions, one of which satisfies a relative smoothness condition and the other one is possibly nonsmooth. A key tool of our methodology is the Bregman forward-backward envelope (BFBE), an exact and continuous penalty function with favorable first- and second-order properties, which enjoys a nonlinear error bound when the objective function satisfies a Lojasiewicz-type property. The proposed algorithm is of linesearch type over the BFBE along user-defined update directions and converges subsequentially to stationary points and globally under the Kurdyka-Lojasiewicz condition. Moreover, when the update directions are superlinear in the sense of Facchinei and Pang [Finite-Dimensional Variational Inequalities and Complementarity Problems, Volume I, Springer, New York, 2003], owing to the given nonlinear error bound unit stepsize is eventually always accepted and the algorithm attains superlinear convergence rates even when the limit point is a nonisolated minimum.

We study a class of countably infinite linear programs (CILPs) whose feasible sets are bounded subsets of appropriately defined spaces of measures. The optimal value, optimal points, and minimal points of these CILPs can be approximated by solving finite-dimensional linear programs. We show how to construct finite-dimensional programs that lead to approximations with easy-to-evaluate error bounds, and we prove that the errors converge to zero as the size of the finite-dimensional programs approaches that of the original problem. We discuss the use of our methods in the computation of the stationary distributions, occupation measures, and exit distributions of Markov chains. Read More:https://epubs.siam.org/doi/10.1137/19M1268847 &nbsp

The linear complementarity problem (LCP) is a general set membership problem that includes quadratic cone programming as a special case. In this work we consider a homogeneous embedding of the LCP, which encodes both the optimality conditions of the original LCP as well as certificates of infeasibility. The resulting problem can be expressed as a monotone inclusion problem involving the sum of two maximal monotone operators, to which we apply Douglas-Rachford splitting. The resulting algorithm has almost identical per-iteration cost as applying the splitting to the LCP directly. Specifically, each iteration requires solving a linear system and projecting onto a convex cone. The resulting algorithm has no explicit hyper-parameters and is able to return a primal-dual solution should one exist or a certificate of infeasibility otherwise. If a sequence of related problems are solved it can easily be warm-started and make use of factorization caching of the linear system. We demonstrate on a range of public and synthetic datasets that for feasible problems our approach tends to be somewhat faster than applying operator splitting directly to the LCP, and in cases of infeasibility our approach can be significantly faster than alternative approaches based on diverging iterates. The algorithm we describe has been implemented in C and is available open-source as SCS-3.0.

Multistage stochastic programming deals with operational and planning problems that involve a sequence of decisions over time while responding to realizations that are uncertain. Algorithms designed to address multistage stochastic linear programming (MSLP) problems often rely upon scenario trees to represent the underlying stochastic process. When this process exhibits stagewise independence, sampling-based techniques, particularly the stochastic dual dynamic programming (SDDP) algorithm, have received wide acceptance. However, these sampling-based methods still operate with a deterministic representation of the problem that uses the so-called sample average approximation. In this work, we present a sequential sampling approach for MSLP problems that allows the decision process to assimilate newly sampled data recursively. We refer to this method as the stochastic dynamic linear programming (SDLP) algorithm. Since we use sequential sampling, the algorithm does not necessitate a priori representation of uncertainty, either through a scenario tree or sample average approximation, both of which require a knowledge/estimation of the underlying distribution. In this regard, SDLP is a sequential sampling approach to address MSLP problems. This method constitutes a generalization of the Stochastic Decomposition (SD) for two-stage SLP models. We employ quadratic regularization for optimization problems in the non-terminal stages. Furthermore, we introduce the notion of basic feasible policies which provide a piecewise-affine solution discovery scheme, that is embedded within the optimization algorithm to identify incumbent solutions used for regularization. Finally, we show that the SDLP algorithm provides a sequence of decisions and corresponding value function estimates along a sequence of state trajectories that asymptotically converge to their optimal counterparts, with probability one.

We propose a unified framework to address a family of classical mixed-integer optimization problems with logically constrained decision variables, including network design, facility location, unit commitment, sparse portfolio selection, binary quadratic optimization, sparse principal analysis and sparse learning problems. These problems exhibit logical relationships between continuous and discrete variables, which are usually reformulated linearly using a big-M formulation. In this work, we challenge this longstanding modeling practice and express the logical constraints in a non-linear way. By imposing a regularization condition, we reformulate these problems as convex binary optimization problems, which are solvable using an outer-approximation procedure. In numerical experiments, we establish that a general-purpose numerical strategy, which combines cutting-plane, first-order and local search methods, solves these problems faster and at a larger scale than state-of-the-art mixed-integer linear or second-order cone methods. Our approach successfully solves network design problems with 100s of nodes and provides solutions up to 40\% better than the state-of-the-art; sparse portfolio selection problems with up to 3,200 securities compared with 400 securities for previous attempts; and sparse regression problems with up to 100,000 covariates., Revised version (including title change). The old title was "A unified approach to mixed-integer optimization: Nonlinear formulations and scalable algorithms"

Clustering is one of the fundamental tasks in data analytics and machine learning. In many situations, different clusterings of the same data set become relevant. For example, different algorithms for the same clustering task may return dramatically different solutions. We are interested in applications in which one clustering has to be transformed into another; e.g., when a gradual transition from an old solution to a new one is required. In this paper, we devise methods for constructing such a transition based on linear programming and network theory. We use a so-called clustering-difference graph to model the desired transformation and provide methods for decomposing the graph into a sequence of elementary moves that accomplishes the transformation. These moves are equivalent to the edge directions, or circuits, of the underlying partition polytopes. Therefore, in addition to a conceptually new metric for measuring the distance between clusterings, we provide new bounds on the circuit diameter of these partition polytopes.

We consider a distributionally robust Partially Observable Markov Decision Process (DR-POMDP), where the distribution of the transition-observation probabilities is unknown at the beginning of each decision period, but their realizations can be inferred using side information at the end of each period after an action being taken. We build an ambiguity set of the joint distribution using bounded moments via conic constraints and seek an optimal policy to maximize the worst-case (minimum) reward for any distribution in the set. We show that the value function of DR-POMDP is piecewise linear convex with respect to the belief state and propose a heuristic search value iteration method for obtaining lower and upper bounds of the value function. We conduct numerical studies and demonstrate the computational performance of our approach via testing instances of a dynamic epidemic control problem. Our results show that DR-POMDP can produce more robust policies under misspecified distributions of transition-observation probabilities as compared to POMDP, but has less costly solutions than robust POMDP. The DR-POMDP policies are also insensitive to varying parameter in the ambiguity set and to noise added to the true transition-observation probability values obtained at the end of each decision period.

Low rank matrix recovery problems appear widely in statistics, combinatorics, and imaging. One celebrated method for solving these problems is to formulate and solve a semidefinite program (SDP). It is often known that the exact solution to the SDP with perfect data recovers the solution to the original low rank matrix recovery problem. It is more challenging to show that an approximate solution to the SDP formulated with noisy problem data acceptably solves the original problem; arguments are usually ad hoc for each problem setting, and can be complex. In this note, we identify a set of conditions that we call simplicity that limit the error due to noisy problem data or incomplete convergence. In this sense, simple SDPs are robust: simple SDPs can be (approximately) solved efficiently at scale; and the resulting approximate solutions, even with noisy data, can be trusted. Moreover, we show that simplicity holds generically, and also for many structured low rank matrix recovery problems, including the stochastic block model, $\mathbb{Z}_2$ synchronization, and matrix completion. Formally, we call an SDP simple if it has a surjective constraint map, admits a unique primal and dual solution pair, and satisfies strong duality and strict complementarity. However, simplicity is not a panacea: we show the Burer-Monteiro formulation of the SDP may have spurious second-order critical points, even for a simple SDP with a rank 1 solution., Comment: 24 pages, 1 figure, and 1 table

Recently, in a series of papers [32,38,39,41], the ratio of $\ell_1$ and $\ell_2$ norms was proposed as a sparsity inducing function for noiseless compressed sensing. In this paper, we further study properties of such model in the noiseless setting, and propose an algorithm for minimizing $\ell_1$/$\ell_2$ subject to noise in the measurements. Specifically, we show that the extended objective function (the sum of the objective and the indicator function of the constraint set) of the model in [32] satisfies the Kurdyka-Lojasiewicz (KL) property with exponent 1/2; this allows us to establish linear convergence of the algorithm proposed in [39, Eq. 11] under mild assumptions. We next extend the $\ell_1$/$\ell_2$ model to handle compressed sensing problems with noise. We establish the solution existence for some of these models under the spherical section property [37,44], and extend the algorithm in [39, Eq. 11] by incorporating moving-balls-approximation techniques [4] for solving these problems. We prove the subsequential convergence of our algorithm under mild conditions, and establish global convergence of the whole sequence generated by our algorithm by imposing additional KL and differentiability assumptions on a specially constructed potential function. Finally, we perform numerical experiments on robust compressed sensing and basis pursuit denoising with residual error measured by $ \ell_2 $ norm or Lorentzian norm via solving the corresponding $\ell_1$/$\ell_2$ models by our algorithm. Our numerical simulations show that our algorithm is able to recover the original sparse vectors with reasonable accuracy.

Recently, there has been a surge of interest in designing variants of the classical Newton-CG in which the Hessian of a (strongly) convex function is replaced by suitable approximations. This is mainly motivated by large-scale finite-sum minimization problems that arise in many machine learning applications. Going beyond convexity, inexact Hessian information has also been recently considered in the context of algorithms such as trust-region or (adaptive) cubic regularization for general non-convex problems. Here, we do that for Newton-MR, which extends the application range of the classical Newton-CG beyond convexity to invex problems. Unlike the convergence analysis of Newton-CG, which relies on spectrum preserving Hessian approximations in the sense of L\"{o}wner partial order, our work here draws from matrix perturbation theory to estimate the distance between the subspaces underlying the exact and approximate Hessian matrices. Numerical experiments demonstrate a great degree of resilience to such Hessian approximations, amounting to a highly efficient algorithm in large-scale problems., Comment: 32 pages, 10 figures

Computing the stationary distributions of a continuous-time Markov chain (CTMC) involves solving a set of linear equations. In most cases of interest, the number of equations is infinite or too large, and the equations cannot be solved analytically or numerically. Several approximation schemes overcome this issue by truncating the state space to a manageable size. In this review, we first give a comprehensive theoretical account of the stationary distributions and their relation to the long-term behaviour of CTMCs that is readily accessible to non-experts and free of irreducibility assumptions made in standard texts. We then review truncation-based approximation schemes for CTMCs with infinite state spaces paying particular attention to the schemes' convergence and the errors they introduce, and we illustrate their performance with an example of a stochastic reaction network of relevance in biology and chemistry. We conclude by elaborating on computational trade-offs associated with error control and several open questions. &nbsp

Smooth convex minimization over the unit trace-norm ball is an important optimization problem in machine learning, signal processing, statistics and other fields, that underlies many tasks in which one wishes to recover a low-rank matrix given certain measurements. While first-order methods for convex optimization enjoy optimal convergence rates, they require in worst-case to compute a full-rank SVD on each iteration, in order to compute the projection onto the trace-norm ball. These full-rank SVD computations however prohibit the application of such methods to large problems. A simple and natural heuristic to reduce the computational cost is to approximate the projection using only a low-rank SVD. This raises the question if, and under what conditions, this simple heuristic can indeed result in provable convergence to the optimal solution. In this paper we show that any optimal solution is a center of a Euclid. ball inside-which the projected-gradient mapping admits rank that is at most the multiplicity of the largest singular value of the gradient vector. Moreover, the radius of the ball scales with the spectral gap of this gradient vector. We show how this readily implies the local convergence (i.e., from a "warm-start" initialization) of standard first-order methods, using only low-rank SVD computations. We also quantify the effect of "over-parameterization", i.e., using SVD computations with higher rank, on the radius of this ball, showing it can increase dramatically with moderately larger rank. We extend our results also to the setting of optimization with trace-norm regularization and optimization over bounded-trace positive semidefinite matrices. Our theoretical investigation is supported by concrete empirical evidence that demonstrates the \textit{correct} convergence of first-order methods with low-rank projections on real-world datasets., Accepted to SIAM Journal on Optimization (SIOPT)

This paper presents smoothing schemes for obtaining approximate stationary points of unconstrained or linearly-constrained composite nonconvex-concave min-max (and hence nonsmooth) problems by applying well-known algorithms to composite smooth approximations of the original problems. More specifically, in the unconstrained (resp. constrained) case, approximate stationary points of the original problem are obtained by applying, to its composite smooth approximation, an accelerated inexact proximal point (resp. quadratic penalty) method presented in a previous paper by the authors. Iteration complexity bounds for both smoothing schemes are also established. Finally, numerical results are given to demonstrate the efficiency of the unconstrained smoothing scheme.

In this work, we present a novel approach for solving stochastic shape optimization problems. Our method is the extension of the classical stochastic gradient method to infinite-dimensional shape manifolds. We prove convergence of the method on Riemannian manifolds and then make the connection to shape spaces. The method is demonstrated on a model shape optimization problem from interface identification. Uncertainty arises in the form of a random partial differential equation, where underlying probability distributions of the random coefficients and inputs are assumed to be known. We verify some conditions for convergence for the model problem and demonstrate the method numerically.

Conditional gradient methods have attracted much attention in both machine learning and optimization communities recently. These simple methods can guarantee the generation of sparse solutions. In addition, without the computation of full gradients, they can handle huge-scale problems sometimes even with an exponentially increasing number of decision variables. This paper aims to significantly expand the application areas of these methods by presenting new conditional gradient methods for solving convex optimization problems with general affine and nonlinear constraints. More specifically, we first present a new constraint extrapolated condition gradient (CoexCG) method that can achieve an ${\cal O}(1/\epsilon^2)$ iteration complexity for both smooth and structured nonsmooth function constrained convex optimization. We further develop novel variants of CoexCG, namely constraint extrapolated and dual regularized conditional gradient (CoexDurCG) methods, that can achieve similar iteration complexity to CoexCG but allow adaptive selection for algorithmic parameters. We illustrate the effectiveness of these methods for solving an important class of radiation therapy treatment planning problems arising from healthcare industry. To the best of our knowledge, all the algorithmic schemes and their complexity results are new in the area of projection-free methods.

In this paper, we provide the first provable linear-time (in the number of non-zero entries of the input) algorithm for approximately solving the generalized trust region subproblem (GTRS) of minimizing a quadratic function over a quadratic constraint under some regular conditions. Our algorithm is motivated by and extends a recent linear-time algorithm for the trust region subproblem by Hazan and Koren [Math. Program., 2016, 158(1-2): 363-381]. However, due to the non-convexity and non-compactness of the feasible region, such an extension is nontrivial. Our main contribution is to demonstrate that under some regular condition, the optimal solution is in a compact and convex set and lower and upper bounds of the optimal value can be computed in linear time. Using these properties, we develop a linear-time algorithm for the GTRS.

We propose a methodology for studying the performance of common splitting methods through semidefinite programming. We prove tightness of the methodology and demonstrate its value by presenting two applications of it. First, we use the methodology as a tool for computer-assisted proofs to prove tight analytical contraction factors for Douglas--Rachford splitting that are likely too complicated for a human to find bare-handed. Second, we use the methodology as an algorithmic tool to computationally select the optimal splitting method parameters by solving a series of semidefinite programs., Comment: Published in the SIAM Journal on Optimization

We analyze the coordinate descent method with a new coordinate selection strategy, called volume sampling. This strategy prescribes selecting subsets of variables of certain size proportionally to the determinants of principal submatrices of the matrix, that bounds the curvature of the objective function. In the particular case, when the size of the subsets equals one, volume sampling coincides with the well-known strategy of sampling coordinates proportionally to their Lipschitz constants. For the coordinate descent with volume sampling, we establish the convergence rates both for convex and strongly convex problems. Our theoretical results show that, by increasing the size of the subsets, it is possible to accelerate the method up to the factor which depends on the spectral gap between the corresponding largest eigenvalues of the curvature matrix. Several numerical experiments confirm our theoretical conclusions.

EXTRA is a popular method for dencentralized distributed optimization and has broad applications. This paper revisits EXTRA. First, we give a sharp complexity analysis for EXTRA with the improved $O\left(\left(\frac{L}{\mu}+\frac{1}{1-\sigma_2(W)}\right)\log\frac{1}{\epsilon(1-\sigma_2(W))}\right)$ communication and computation complexities for $\mu$-strongly convex and $L$-smooth problems, where $\sigma_2(W)$ is the second largest singular value of the weight matrix $W$. When the strong convexity is absent, we prove the $O\left(\left(\frac{L}{\epsilon}+\frac{1}{1-\sigma_2(W)}\right)\log\frac{1}{1-\sigma_2(W)}\right)$ complexities. Then, we use the Catalyst framework to accelerate EXTRA and obtain the $O\left(\sqrt{\frac{L}{\mu(1-\sigma_2(W))}}\log\frac{ L}{\mu(1-\sigma_2(W))}\log\frac{1}{\epsilon}\right)$ communication and computation complexities for strongly convex and smooth problems and the $O\left(\sqrt{\frac{L}{\epsilon(1-\sigma_2(W))}}\log\frac{1}{\epsilon(1-\sigma_2(W))}\right)$ complexities for non-strongly convex ones. Our communication complexities of the accelerated EXTRA are only worse by the factors of $\left(\log\frac{L}{\mu(1-\sigma_2(W))}\right)$ and $\left(\log\frac{1}{\epsilon(1-\sigma_2(W))}\right)$ from the lower complexity bounds for strongly convex and non-strongly convex problems, respectively.

The paper studies generalized differentiability properties of the marginal function of parametric optimal control problems of semilinear elliptic partial differential equations. We establish upper estimates for the regular and the limiting subgradients of the marginal function. With some additional assumptions, we show that the solution map of the perturbed optimal control problems has local upper H\"{o}lderian selections. This leads to a lower estimate for the regular subdifferential of the marginal function.

In this paper we consider the unconstrained minimization problem of a smooth function in ${\mathbb{R}}^n$ in a setting where only function evaluations are possible. We design a novel randomized derivative-free algorithm --- the stochastic three points (STP) method --- and analyze its iteration complexity. At each iteration, STP generates a random search direction according to a certain fixed probability law. Our assumptions on this law are very mild: roughly speaking, all laws which do not concentrate all measure on any halfspace passing through the origin will work. For instance, we allow for the uniform distribution on the sphere and also distributions that concentrate all measure on a positive spanning set. Given a current iterate $x$, STP compares the objective function at three points: $x$, $x+\alpha s$ and $x-\alpha s$, where $\alpha>0$ is a stepsize parameter and $s$ is the random search direction. The best of these three points is the next iterate. We analyze the method STP under several stepsize selection schemes (fixed, decreasing, estimated through finite differences, etc). We study non-convex, convex and strongly convex cases.

In this paper, we propose new accelerated methods for smooth convex optimization, called contracting proximal methods. At every step of these methods, we need to minimize a contracted version of the objective function augmented by a regularization term in the form of Bregman divergence. We provide global convergence analysis for a general scheme admitting inexactness in solving the auxiliary subproblem. In the case of using for this purpose high-order tensor methods, we demonstrate an acceleration effect for both convex and uniformly convex composite objective functions. Thus, our construction explains acceleration for methods of any order starting from one. The augmentation of the number of calls of oracle due to computing the contracted proximal steps is limited by the logarithmic factor in the worst-case complexity bound.

Consider the minimization of a nonconvex differentiable function over a polyhedron. A popular primal-dual first-order method for this problem is to perform a gradient projection iteration for the augmented Lagrangian function and then update the dual multiplier vector using the constraint residual. However, numerical examples show that this approach can exhibit "oscillation" and may not converge. In this paper, we propose a proximal alternating direction method of multipliers for the multi-block version of this problem. A distinctive feature of this method is the introduction of a "smoothed" (i.e., exponentially weighted) sequence of primal iterates, and the inclusion, at each iteration, to the augmented Lagrangian function a quadratic proximal term centered at the current smoothed primal iterate. The resulting proximal augmented Lagrangian function is inexactly minimized (via a gradient projection step) at each iteration while the dual multiplier vector is updated using the residual of the linear constraints. When the primal and dual stepsizes are chosen sufficiently small, we show that suitable "smoothing" can stabilize the "oscillation", and the iterates of the new proximal ADMM algorithm converge to a stationary point under some mild regularity conditions. Furthermore, when the objective function is quadratic, we establish the linear convergence of the algorithm. Our proof is based on a new potential function and a novel use of error bounds., Comment: 31 pages,6 figures. Accepted for publication by SIAM Journal on Optimization

This paper is devoted to the study of generalized subdifferentials of spectral functions over Euclidean Jordan algebras. Spectral functions appear often in optimization problems playing the role of "regularizer", "barrier", "penalty function" and many others. We provide formulae for the regular, approximate and horizon subdifferentials of spectral functions. In addition, under local lower semicontinuity, we also furnish a formula for the Clarke subdifferential, thus extending an earlier result by Baes. As application, we compute the generalized subdifferentials of the function that maps an element to its k-th largest eigenvalue. Furthermore, in connection with recent approaches for nonsmooth optimization, we present a study of the Kurdyka-Lojasiewicz (KL) property for spectral functions and prove a transfer principle for the KL-exponent. In our proofs, we make extensive use of recent tools such as the commutation principle of Ram\'irez, Seeger and Sossa and majorization principles developed by Gowda., Comment: 26 pages. Some minor fixes and trimming. Accepted for publication at the SIAM Journal on Optimization

In this paper we study $p$-order methods for unconstrained minimization of convex functions that are $p$-times differentiable ($p\geq 2$) with $\nu$-H\"{o}lder continuous $p$th derivatives. We propose tensor schemes with and without acceleration. For the schemes without acceleration, we establish iteration complexity bounds of $\mathcal{O}\left(\epsilon^{-1/(p+\nu-1)}\right)$ for reducing the functional residual below a given $\epsilon\in (0,1)$. Assuming that $\nu$ is known, we obtain an improved complexity bound of $\mathcal{O}\left(\epsilon^{-1/(p+\nu)}\right)$ for the corresponding accelerated scheme. For the case in which $\nu$ is unknown, we present a universal accelerated tensor scheme with iteration complexity of $\mathcal{O}\left(\epsilon^{-p/[(p+1)(p+\nu-1)]}\right)$. A lower complexity bound of $\mathcal{O}\left(\epsilon^{-2/[3(p+\nu)-2]}\right)$ is also obtained for this problem class., Comment: arXiv admin note: text overlap with arXiv:1907.07053

We consider distributed smooth nonconvex unconstrained optimization over networks, modeled as a connected graph. We examine the behavior of distributed gradient-based algorithms near strict saddle points. Specifically, we establish that (i) the renowned Distributed Gradient Descent (DGD) algorithm likely converges to a neighborhood of a Second-order Stationary (SoS) solution; and (ii) the more recent class of distributed algorithms based on gradient tracking--implementable also over digraphs--likely converges to exact SoS solutions, thus avoiding (strict) saddle-points. Furthermore, new convergence rate results to first-order critical points is established for the latter class of algorithms., Comment: Final version, to appear on SIAM J. on Optimization

A multiobjective optimization problem is simplicial if the Pareto set and front are homeomorphic to a simplex and, under the homeomorphisms, each face of the simplex corresponds to the Pareto set and front of a subproblem. In this paper, we show that strongly convex problems are simplicial under a mild assumption on the ranks of the differentials of the objective mappings. We further prove that one can make any strongly convex problem satisfy the assumption by a generic linear perturbation, provided that the dimension of the source is sufficiently larger than that of the target. We demonstrate that the location problems, a biological modeling, and the ridge regression can be reduced to multiobjective strongly convex problems via appropriate transformations preserving the Pareto ordering and the topology., Comment: 21 pages. Remarks 4.4 and 4.5 are added. A new application is given in section 5.3. Introduction is also revised accordingly

In this paper, we propose two novel non-stationary first-order primal-dual algorithms to solve nonsmooth composite convex optimization problems. Unlike existing primal-dual schemes where the parameters are often fixed, our methods use pre-defined and dynamic sequences for parameters. We prove that our first algorithm can achieve $\mathcal{O}(1/k)$ convergence rate on the primal-dual gap, and primal and dual objective residuals, where $k$ is the iteration counter. Our rate is on the non-ergodic (i.e., the last iterate) sequence of the primal problem and on the ergodic (i.e., the averaging) sequence of the dual problem, which we call semi-ergodic rate. By modifying the step-size update rule, this rate can be boosted even faster on the primal objective residual. When the problem is strongly convex, we develop a second primal-dual algorithm that exhibits $\mathcal{O}(1/k^2)$ convergence rate on the same three types of guarantees. Again by modifying the step-size update rule, this rate becomes faster on the primal objective residual. Our primal-dual algorithms are the first ones to achieve such fast convergence rate guarantees under mild assumptions compared to existing works, to the best of our knowledge. As byproducts, we apply our algorithms to solve constrained convex optimization problems and prove the same convergence rates on both the objective residuals and the feasibility violation. We still obtain at least $\mathcal{O}(1/k^2)$ rates even when the problem is "semi-strongly" convex. We verify our theoretical results via two well-known numerical examples., Comment: 30 pages, 2 figures