Search

This paper considers the problem of solving a special quartic–quadratic optimization problem with a single sphere constraint, namely, finding a global and local minimizer of $$\frac{1}{2}\mathbf {z}^{*}A\mathbf {z}+\frac{\beta }{2}\sum _{k=1}^{n}|z_{k}|^{4}$$ such that $$\Vert \mathbf {z}\Vert _{2}=1$$ . This problem spans multiple domains including quantum mechanics and chemistry sciences and we investigate the geometric properties of this optimization problem. Fourth-order optimality conditions are derived for characterizing local and global minima. When the matrix in the quadratic term is diagonal, the problem has no spurious local minima and global solutions can be represented explicitly and calculated in $$O(n\log {n})$$ operations. When A is a rank one matrix, the global minima of the problem are unique under certain phase shift schemes. The strict-saddle property, which can imply polynomial time convergence of second-order-type algorithms, is established when the coefficient $$\beta $$ of the quartic term is either at least $$O(n^{3/2})$$ or not larger than O(1). Finally, the Kurdyka–Łojasiewicz exponent of quartic–quadratic problem is estimated and it is shown that the largest exponent is at least 1/4 for a broad class of stationary points.

Many iterative methods in applied mathematics can be thought of as fixed-point iterations, and such algorithms are usually analyzed analytically, with inequalities. In this paper, we present a geometric approach to analyzing contractive and nonexpansive fixed point iterations with a new tool called the scaled relative graph (SRG). The SRG provides a correspondence between nonlinear operators and subsets of the 2D plane. Under this framework, a geometric argument in the 2D plane becomes a rigorous proof of convergence., Published in Mathematical Programming

Mean-field games arise in various fields including economics, engineering, and machine learning. They study strategic decision making in large populations where the individuals interact via certain mean-field quantities. The ground metrics and running costs of the games are of essential importance but are often unknown or only partially known. In this paper, we propose mean-field game inverse-problem models to reconstruct the ground metrics and interaction kernels in the running costs. The observations are the macro motions, to be specific, the density distribution, and the velocity field of the agents. They can be corrupted by noise to some extent. Our models are PDE constrained optimization problems, which are solvable by first-order primal-dual methods. Besides, we apply Bregman iterations to find the optimal model parameters. We numerically demonstrate that our model is both efficient and robust to noise.

The method of block coordinate gradient descent (BCD) has been a powerful method for large-scale optimization. This paper considers the BCD method that successively updates a series of blocks selected according to a Markov chain. This kind of block selection is neither i.i.d. random nor cyclic. On the other hand, it is a natural choice for some applications in distributed optimization and Markov decision process, where i.i.d. random and cyclic selections are either infeasible or very expensive. By applying mixing-time properties of a Markov chain, we prove convergence of Markov chain BCD for minimizing Lipschitz differentiable functions, which can be nonconvex. When the functions are convex and strongly convex, we establish both sublinear and linear convergence rates, respectively. We also present a method of Markov chain inertial BCD. Finally, we discuss potential applications.

We design fast numerical methods for Hamilton–Jacobi equations in density space (HJD), which arises in optimal transport and mean field games. We proposes an algorithm using a generalized Hopf formula in density space. The formula helps transforming a problem from an optimal control problem in density space, which are constrained minimizations supported on both spatial and time variables, to an optimization problem over only one spatial variable. This transformation allows us to compute HJD efficiently via multi-level approaches and coordinate descent methods. Rigorous derivation of the Hopf formula is provided under restricted assumptions and for a relatively narrow case; meanwhile our practical investigation allows us to conjecture that the actual range of applicability should be wider, and therefore we conjecture the formula can be applied to a wider class of practical examples.

It is known that operator splitting methods based on forward–backward splitting, Douglas–Rachford splitting, and Davis–Yin splitting decompose difficult optimization problems into simpler subproblems under proper convexity and smoothness assumptions. In this paper, we identify an envelope (an objective function), whose gradient descent iteration under a variable metric coincides with Davis–Yin splitting iteration. This result generalizes the Moreau envelope for proximal-point iteration and the envelopes for forward–backward splitting and Douglas–Rachford splitting iterations identified by Patrinos, Stella, and Themelis. Based on the new envelope and the stable–center manifold theorem, we further show that, when forward–backward splitting or Douglas–Rachford splitting iterations start from random points, they avoid all strict saddle points with probability one. This result extends the similar results by Lee et al. from gradient descent to splitting methods.

Primal–dual hybrid gradient (PDHG) and alternating direction method of multipliers (ADMM) are popular first-order optimization methods. They are easy to implement and have diverse applications. As first-order methods, however, they are sensitive to problem conditions and can struggle to reach the desired accuracy. To improve their performance, researchers have proposed techniques such as diagonal preconditioning and inexact subproblems. This paper realizes additional speedup about one order of magnitude. Specifically, we choose general (non-diagonal) preconditioners that are much more effective at reducing the total numbers of PDHG/ADMM iterations than diagonal ones. Although the subproblems may lose their closed-form solutions, we show that it suffices to solve each subproblem approximately with a few proximal-gradient iterations or a few epochs of proximal block-coordinate descent, which are simple and have closed-form steps. Global convergence of this approach is proved when the inner iterations are fixed. Our method opens the choices of preconditioners and maintains both low per-iteration cost and global convergence. Consequently, on several typical applications of primal–dual first-order methods, we obtain 4–95 $$\times $$ speedup over the existing state-of-the-art.

In this paper, we analyze the convergence of the alternating direction method of multipliers (ADMM) for minimizing a nonconvex and possibly nonsmooth objective function, $\phi(x_0,\ldots,x_p,y)$, subject to coupled linear equality constraints. Our ADMM updates each of the primal variables $x_0,\ldots,x_p,y$, followed by updating the dual variable. We separate the variable $y$ from $x_i$'s as it has a special role in our analysis. The developed convergence guarantee covers a variety of nonconvex functions such as piecewise linear functions, $\ell_q$ quasi-norm, Schatten-$q$ quasi-norm ($0, Comment: 33 pages, 1 figure, Accepted by Journal of Scientific Computing

In this paper, we present a method for identifying infeasible, unbounded, and pathological conic programs based on Douglas–Rachford splitting. When an optimization program is infeasible, unbounded, or pathological, the iterates of Douglas–Rachford splitting diverge. Somewhat surprisingly, such divergent iterates still provide useful information, which our method uses for identification. In addition, for strongly infeasible problems the method produces a separating hyperplane and informs the user on how to minimally modify the given problem to achieve strong feasibility. As a first-order method, the proposed algorithm relies on simple subroutines, and therefore is simple to implement and has low per-iteration cost.

Recent years have witnessed the surge of asynchronous parallel (async-parallel) iterative algorithms due to problems involving very large-scale data and a large number of decision variables. Because of asynchrony, the iterates are computed with outdated information, and the age of the outdated information, which we call delay, is the number of times it has been updated since its creation. Almost all recent works prove convergence under the assumption of a finite maximum delay and set their stepsize parameters accordingly. However, the maximum delay is practically unknown. This paper presents convergence analysis of an async-parallel method from a probabilistic viewpoint, and it allows for large unbounded delays. An explicit formula of stepsize that guarantees convergence is given depending on delays' statistics. With $p+1$ identical processors, we empirically measured that delays closely follow the Poisson distribution with parameter $p$, matching our theoretical model, and thus the stepsize can be set accordingly. Simulations on both convex and nonconvex optimization problems demonstrate the validness of our analysis and also show that the existing maximum-delay induced stepsize is too conservative, often slowing down the convergence of the algorithm., Comment: accepted to JORSC

We propose a new algorithm to approximate the Earth Mover’s distance (EMD). Our main idea is motivated by the theory of optimal transport, in which EMD can be reformulated as a familiar $$L_1$$ type minimization. We use a regularization which gives us a unique solution for this $$L_1$$ type problem. The new regularized minimization is very similar to problems which have been solved in the fields of compressed sensing and image processing, where several fast methods are available. In this paper, we adopt a primal-dual algorithm designed there, which uses very simple updates at each iteration and is shown to converge very rapidly. Several numerical examples are provided.

Image restoration is one of the essential tasks in image processing. In order to restore images from blurs and noise while also preserving their edges, one often applies total variation (TV) minimization. Cauchy noise, which frequently appears in engineering applications, is a kind of impulsive and non-Gaussian noise. Removing Cauchy noise can be achieved by solving a nonconvex TV minimization problem, which is difficult due to its nonconvexity and nonsmoothness. In this paper, we adapt recent results in the literature and develop a specific alternating direction method of multiplier to solve this problem. Theoretically, we establish the convergence of our method to a stationary point. Experimental results demonstrate that the proposed method is competitive with other methods in visual and quantitative measures. In particular, our method achieves higher PSNRs for 0.5 dB on average.

In this paper we develop a parallel method for solving possibly non-convex time-dependent Hamilton---Jacobi equations arising from optimal control and differential game problems. The subproblems are independent so they can be implemented in an embarrassingly parallel fashion, which usually has an ideal parallel speedup. The algorithm is proposed to overcome the curse of dimensionality (Bellman in Adaptive control processes: a guided tour. Princeton University Press, Princeton, 1961; Dynamic programming. Princeton University Press, Princeton, 1957) when solving HJ PDE . We extend previous work Chow et al. (Algorithm for overcoming the curse of dimensionality for certain non-convex Hamilton---Jacobi equations, Projections and differential games, UCLA CAM report, pp 16---27, 2016) and Darbon and Osher (Algorithms for overcoming the curse of dimensionality for certain Hamilton---Jacobi equations arising in control theory and elsewhere, UCLA CAM report, pp 15---50, 2015) and apply a generalized Hopf formula to solve HJ PDE involving time-dependent and perhaps non-convex Hamiltonians. We suggest a coordinate descent method for the minimization procedure in the Hopf formula. This method is preferable when even the evaluation of the function value itself requires some computational effort, and also when we handle higher dimensional optimization problem. For an integral with respect to time inside the generalized Hopf formula, we suggest using a numerical quadrature rule. Together with our suggestion to perform numerical differentiation to minimize the number of calculation procedures in each iteration step, we are bound to have numerical errors in our computations. These errors can be effectively controlled by choosing an appropriate mesh-size in time and the method does not use a mesh in space. The use of multiple initial guesses is suggested to overcome possibly multiple local extrema in the case when non-convex Hamiltonians are considered. Our method is expected to have wide application in control theory and differential game problems, and elsewhere.

This paper introduces a parallel and distributed extension to the alternating direction method of multipliers (ADMM) for solving convex problem: minimize $\sum_{i=1}^N f_i(x_i)$ subject to $\sum_{i=1}^N A_i x_i=c, x_i\in \mathcal{X}_i$. The algorithm decomposes the original problem into N smaller subproblems and solves them in parallel at each iteration. This Jacobian-type algorithm is well suited for distributed computing and is particularly attractive for solving certain large-scale problems. This paper introduces a few novel results. Firstly, it shows that extending ADMM straightforwardly from the classic Gauss-Seidel setting to the Jacobian setting, from 2 blocks to N blocks, will preserve convergence if matrices $A_i$ are mutually near-orthogonal and have full column-rank. Secondly, for general matrices $A_i$, this paper proposes to add proximal terms of different kinds to the N subproblems so that the subproblems can be solved in flexible and efficient ways and the algorithm converges globally at a rate of o(1/k). Thirdly, a simple technique is introduced to improve some existing convergence rates from O(1/k) to o(1/k). In practice, some conditions in our convergence theorems are conservative. Therefore, we introduce a strategy for dynamically tuning the parameters in the algorithm, leading to substantial acceleration of the convergence in practice. Numerical results are presented to demonstrate the efficiency of the proposed method in comparison with several existing parallel algorithms. We implemented our algorithm on Amazon EC2, an on-demand public computing cloud, and report its performance on very large-scale basis pursuit problems with distributed data., Comment: The second version has a new title

The article “On the Convergence of Asynchronous Parallel Iteration with Unbounded Delays”, written by Zhimin Peng, Yangyang Xu, Ming Yan, Wotao Yin was originally published electronically on the publisher’s Internet portal (currently SpringerLink) on 9th December 2017 without open access.

The formulation $$\begin{aligned} \min _{x,y} ~f(x)+g(y),\quad \text{ subject } \text{ to } Ax+By=b, \end{aligned}$$minx,yf(x)+g(y),subjecttoAx+By=b,where f and g are extended-value convex functions, arises in many application areas such as signal processing, imaging and image processing, statistics, and machine learning either naturally or after variable splitting. In many common problems, one of the two objective functions is strictly convex and has Lipschitz continuous gradient. On this kind of problem, a very effective approach is the alternating direction method of multipliers (ADM or ADMM), which solves a sequence of f/g-decoupled subproblems. However, its effectiveness has not been matched by a provably fast rate of convergence; only sublinear rates such as O(1 / k) and $$O(1/k^2)$$O(1/k2) were recently established in the literature, though the O(1 / k) rates do not require strong convexity. This paper shows that global linear convergence can be guaranteed under the assumptions of strong convexity and Lipschitz gradient on one of the two functions, along with certain rank assumptions on A and B. The result applies to various generalizations of ADM that allow the subproblems to be solved faster and less exactly in certain manners. The derived rate of convergence also provides some theoretical guidance for optimizing the ADM parameters. In addition, this paper makes meaningful extensions to the existing global convergence theory of ADM generalizations.

Surface conformal maps between genus-0 surfaces play important roles in applied mathematics and engineering, with applications in medical image analysis and computer graphics. Previous work (Gu and Yau in Commun Inf Syst 2(2):121---146, 2002) introduces a variational approach, where global conformal parameterization of genus-0 surfaces was addressed through minimizing the harmonic energy, with two weaknesses: its gradient descent iteration is slow, and its solutions contain undesired parameterization foldings when the underlying surface has long sharp features. In this paper, we propose an algorithm that significantly accelerates the harmonic energy minimization and a method that iteratively removes foldings by taking advantages of the weighted Laplace---Beltrami eigen-projection. Experimental results show that the proposed approaches compute genus-0 surface harmonic maps much faster than the existing algorithm in Gu and Yau (Commun Inf Syst 2(2):121---146, 2002) and the new results contain no foldings.

Based on the classic augmented Lagrangian multiplier method, we propose, analyze and test an algorithm for solving a class of equality-constrained non-smooth optimization problems (chiefly but not necessarily convex programs) with a particular structure. The algorithm effectively combines an alternating direction technique with a nonmonotone line search to minimize the augmented Lagrangian function at each iteration. We establish convergence for this algorithm, and apply it to solving problems in image reconstruction with total variation regularization. We present numerical results showing that the resulting solver, called TVAL3, is competitive with, and often outperforms, other state-of-the-art solvers in the field.

Minimization with orthogonality constraints (e.g., $$X^\top X = I$$) and/or spherical constraints (e.g., $$\Vert x\Vert _2 = 1$$) has wide applications in polynomial optimization, combinatorial optimization, eigenvalue problems, sparse PCA, p-harmonic flows, 1-bit compressive sensing, matrix rank minimization, etc. These problems are difficult because the constraints are not only non-convex but numerically expensive to preserve during iterations. To deal with these difficulties, we apply the Cayley transform--a Crank-Nicolson-like update scheme--to preserve the constraints and based on it, develop curvilinear search algorithms with lower flops compared to those based on projections and geodesics. The efficiency of the proposed algorithms is demonstrated on a variety of test problems. In particular, for the maxcut problem, it exactly solves a decomposition formulation for the SDP relaxation. For polynomial optimization, nearest correlation matrix estimation and extreme eigenvalue problems, the proposed algorithms run very fast and return solutions no worse than those from their state-of-the-art algorithms. For the quadratic assignment problem, a gap 0.842 % to the best known solution on the largest problem "tai256c" in QAPLIB can be reached in 5 min on a typical laptop.

The matrix completion problem is to recover a low-rank matrix from a subset of its entries. The main solution strategy for this problem has been based on nuclear-norm minimization which requires computing singular value decompositions—a task that is increasingly costly as matrix sizes and ranks increase. To improve the capacity of solving large-scale problems, we propose a low-rank factorization model and construct a nonlinear successive over-relaxation (SOR) algorithm that only requires solving a linear least squares problem per iteration. Extensive numerical experiments show that the algorithm can reliably solve a wide range of problems at a speed at least several times faster than many nuclear-norm minimization algorithms. In addition, convergence of this nonlinear SOR algorithm to a stationary point is analyzed.

This short article analyzes an interesting property of the Bregman iterative procedure, which is equivalent to the augmented Lagrangian method, for minimizing a convex piece-wise linear function J(x) subject to linear constraints Ax=b. The procedure obtains its solution by solving a sequence of unconstrained subproblems of minimizing $J(x)+\frac{1}{2}\|Ax-b^{k}\|_{2}^{2}$ , where b k is iteratively updated. In practice, the subproblem at each iteration is solved at a relatively low accuracy. Let w k denote the error introduced by early stopping a subproblem solver at iteration k. We show that if all w k are sufficiently small so that Bregman iteration enters the optimal face, then while on the optimal face, Bregman iteration enjoys an interesting error-forgetting property: the distance between the current point $\bar{x}^{k}$ and the optimal solution set X ? is bounded by ?w k+1?w k ?, independent of the previous errors w k?1,w k?2,?,w 1. This property partially explains why the Bregman iterative procedure works well for sparse optimization and, in particular, for ? 1-minimization. The error-forgetting property is unique to J(x) that is a piece-wise linear function (also known as a polyhedral function), and the results of this article appear to be new to the literature of the augmented Lagrangian method.

This paper introduces an algorithm for the nonnegative matrix factorization-and-completion problem, which aims to find nonnegative low-rank matrices X and Y so that the product XY approximates a nonnegative data matrix M whose elements are partially known (to a certain accuracy). This problem aggregates two existing problems: (i) nonnegative matrix factorization where all entries of M are given, and (ii) low-rank matrix completion where nonnegativity is not required. By taking the advantages of both nonnegativity and low-rankness, one can generally obtain superior results than those of just using one of the two properties. We propose to solve the non-convex constrained least-squares problem using an algorithm based on the classical alternating direction augmented Lagrangian method. Preliminary convergence properties of the algorithm and numerical simulation results are presented. Compared to a recent algorithm for nonnegative matrix factorization, the proposed algorithm produces factorizations of similar quality using only about half of the matrix entries. On tasks of recovering incomplete grayscale and hyperspectral images, the proposed algorithm yields overall better qualities than those produced by two recent matrix-completion algorithms that do not exploit nonnegativity.

We present an alternating direction dual augmented Lagrangian method for solving semidefinite programming (SDP) problems in standard form. At each iteration, our basic algorithm minimizes the augmented Lagrangian function for the dual SDP problem sequentially, first with respect to the dual variables corresponding to the linear constraints, and then with respect to the dual slack variables, while in each minimization keeping the other variables fixed, and then finally it updates the Lagrange multipliers (i.e., primal variables). Convergence is proved by using a fixed-point argument. For SDPs with inequality constraints and positivity constraints, our algorithm is extended to separately minimize the dual augmented Lagrangian function over four sets of variables. Numerical results for frequency assignment, maximum stable set and binary integer quadratic programming problems demonstrate that our algorithms are robust and very efficient due to their ability or exploit special structures, such as sparsity and constraint orthogonality in these problems.

In this paper, we present a novel image-based technique that transfers illumination from a source face image to a target face image based on the Logarithmic Total Variation (LTV) model. Our method does not require any prior information regarding the lighting conditions or the 3D geometries of the underlying faces. We first use a Radial Basis Functions (RBFs)-based deformation technique to align key facial features of the reference 2D face with those of the target face. Then, we employ the LTV model to factorize each of the two aligned face images to an illumination-dependent component and an illumination-invariant component. Finally, illumination transferring is achieved by replacing the illumination-dependent component of the target face by that of the reference face. We tested our technique on numerous grayscale and color face images from various face datasets including the Yale face Database, as well as the application of illumination-preserved face coloring.

Decentralized optimization has attracted much research interest for resource-limited networked multi-agent systems in recent years. Decentralized consensus optimization, which is one of the decentralized optimization problems of great practical importance, minimizes an objective function that is the sum of the terms from individual agents over a set of variables on which all the agents should reach a consensus. This problem can be reformulated into an equivalent model with two blocks of variables, which can then be solved by the alternating direction method (ADM) with only communications between neighbor nodes. Motivated by a recently emerged class of so-called multi-block ADMs, this article demonstrates that it is more natural to reformulate a decentralized consensus optimization problem to one with multiple blocks of variables and solve it by a multi-block ADM. In particular, we focus on the multi-block ADM with parallel splitting, which has easy decentralized implementation. Convergence rate is analyzed in the setting of average consensus, and the relation between two-block and multi-block ADMs are studied. Numerical experiments demonstrate the effectiveness of the multi-block ADM with parallel splitting in terms of speed and communication cost and show that it has better network scalability.