A Low-Rank Augmented Lagrangian Method for Doubly Nonnegative Relaxations of Mixed-Binary Quadratic Programs

In this dissertation we develop numerical algorithms and software tools to facilitate parallel solutions of nonlinear programming (NLP) problems. In particular, we address large-scale, block-structured problems with an intrinsic decomposable configuration. These problems arise in a great number of engineering applications, including parameter estimation, optimal control, network optimization, and stochastic programming. The structure of these problems can be leveraged by optimization solvers to accelerate solutions and overcome memory limitations, and we propose variants to two classes of optimization algorithms: augmented Lagrangian (AL) schemes and Schur-complement interior-point methods. The convergence properties of augmented Lagrangian decomposition schemes like the alternating direction method of multipliers (ADMM) and progressive hedging (PH) are well established for convex optimization but convergence guarantees in non-convex settings are still poorly understood. In practice, however, ADMM and PH often perform satisfactorily in complex non-convex NLPs. In this work, we study connections between the method of multipliers (MM), ADMM, and PH to derive benchmarking metrics that explain why PH and ADMM work in practice. We illustrate the concepts using challenging dynamic optimization problems. Our exposition seeks to establish more formalism in benchmarking ADMM, PH, and AL schemes and to motivate algorithmic improvements. The effectiveness of nonlinear interior-point solvers for solving large-scale problems relies quite heavily on the solution of the underlying linear algebra systems. The schur-complement decomposition is very effective for parallelizing the solution of linear systems with modest coupling. However, for systems with large number of coupling variables the schur-complement method does not scale favorably. We implement an approach that uses a Krylov solver (GMRES) preconditioned with ADMM to solve block-structured linear systems that arise in the interior-point method. We show that this ADMM-GMRES approach overcomes the well-known scalability issues of Schur decomposition. One important drawback of using decomposition approaches like ADMM and PH is their convergence rate. Unlike Schur-complement interior-point algorithms that have super-linear convergence, augmented Lagrangian approaches typically exhibit linear and sublinear rates. We exploit connections between ADMM and the Schur-complement decomposition to derive an accelerated version of ADMM. Specifically, we study the effectiveness of performing a Newton-Raphson algorithm to compute multiplier estimates for augmented Lagrangian methods. We demonstrate using two-stage stochastic programming problems that our multiplier update achieves convergence in fewer iterations for MM on general nonlinear problems. In the case of ADMM, the newton update significantly reduces the number of subproblem solves for convex quadratic programs (QPs). Moreover, we show that using newton multiplier updates makes the method robust to the selection of the penalty parameter. Traditionally, state-of-the-art optimization solvers are implemented in low-level programming languages. In our experience, the development of decomposition algorithms in these frameworks is challenging. They present a steep learning curve and can slow the development and testing of new numerical algorithms. To mitigate these challenges, we developed PyNumero, a new open source framework implemented in Python and C++. The package seeks to facilitate development of optimization algorithms for large-scale optimization within a high-level programming environment while at the same time minimizing the computational burden of using Python. The efficiency of PyNumero is illustrated by implementing algorithms for problems arising in stochastic programming and optimal control. Timing results are presented for both serial and parallel implementations. Our computational studies demonstrate that with the appropriate balance between compiled code and Python, efficient implementations of optimization algorithms are achievable in these high-level languages.

Although frequently encountered in many practical applications, singular nonlinear optimization has been always recognized as a difficult problem. In the last decades, classical numerical techniques have been proposed to deal with the singular problem. However, the issue of numerical instability and high computational complexity has not found a satisfactory solution so far. In this paper, we consider the singular optimization problem with bounded variables constraint rather than the common unconstraint model. A novel neural network model was proposed for solving the problem of singular convex optimization with bounded variables. Under the assumption of rank one defect, the original difficult problem is transformed into nonsingular constrained optimization problem by enforcing a tensor term. By using the augmented Lagrangian method and the projection technique, it is proven that the proposed continuous model is convergent to the solution of the singular optimization problem. Numerical simulation further confirmed the effectiveness of the proposed neural network approach.

The linear quadratic optimal control problem under quadratic constraints is an optimization problem over a generally non-convex set. Yakubovich [8] and Megretski [4] have studied this problem, and they show how it may be translated into a two-stage, convex optimization problem. In this paper we study the linear quadratic control problem under quadratic constraints for generalized first-order systems. As in the state-space case the linear quadratic control problem without quadratic constraints may be solved in terms of a linear matrix inequality. Subsequently, we use the results from [8] to derive a linear matrix inequality characterizing the linear quadratic optimal behaviour under quadratic constraints.

<p style='text-indent:20px;'>Total generalization variation (TGV) is a very powerful and important regularization for various inverse problems and computer vision tasks. In this paper, we propose a semismooth Newton based augmented Lagrangian method for solving this problem. The augmented Lagrangian method (also called as method of multipliers) is widely used for lots of smooth or nonsmooth variational problems. However, its efficiency heavily depends on solving the corresponding coupled and nonlinear system together and simultaneously. With efficient primal-dual semismooth Newton methods for the challenging and highly coupled nonlinear subproblems involving total generalized variation, we develop a highly efficient and competitive augmented Lagrangian method compared with some fast first-order method. With the analysis of the metric subregularities of the corresponding functions, we give both the global convergence and local linear convergence rate for the proposed augmented Lagrangian methods.</p>

The term, "Waveform inversion" (WFI), refers to a collection of techniques that use the information from seismic data to derive high-fidelity earth models for seismic imaging. The attractiveness of WFI lies mainly in its lack of approximations, at least in a theoretical sense, in contrast to other model determination techniques such as semblance or tomography. However, a whole raft of approximations must be made to make the technique viable with today'scomputing technology and restrictions of seismic acquisition. These are collectively referred to as "waveform inversion strategies" and in this paper we mainly discuss regularization and preconditioning strategies. Because the wavefields need to be accurately modeled to represent the kinematics of all the waves during WFI iterations, the effects of anisotropy often help to improve WFI results. In this paper, forward modeling and its adjoint computation are based on acoustic wave equations in vertical transversely isotropic (VTI)media. We introduce a multi-parameter inversion for P-wave velocity and anisotropy parameters. WFI is a highly nonlinear, ill-posed problem. We introduce additional information and turn the unconstrained optimization problem into a constrained optimization problem in order to reduce the ill-posedness. The geophysics of the problem leads to appropriate constraints, such as restriction of model parameters, or information from well logs. In this paper, we use well logs as constraints and solve the problem using the augmented Lagrangian method (ALM), a mathematical method that replaces a constrained optimization problem by a series of unconstrained problems. The ALM with well constraints aims at preserving velocity characteristics from welllogs and providing us with more reliable velocity updates. This paper presents the acoustic anisotropic WFI implementation using ALM with well constraints. It also discusses practical strategies for regularization and preconditioning and their influences on the models that are obtained from WFI. We illustrate these approaches on a 2D synthetic Marmousi example and another application to 3D VSOOBC data from the Green Canyon area of the Gulf of Mexico. From the results, we show that multi-parameter VTI WFI with ALM provides us with more useful and reliable model updates. To further evaluate our WFI results, we also compare offset gathers and RTM images and illustrate their significant improvement susing updated models generated from WFI.

The goal of this software is to solve a doubly nonnegative (DNN) relaxation for the quadratic assignment problem (QAP).

In this paper, we show that the quadratic assignment problem (QAP) can be reformulated to an equivalent rank constrained doubly nonnegative (DNN) problem. Under the framework of the difference of convex functions (DC) approach, a semi-proximal DC algorithm is proposed for solving the relaxation of the rank constrained DNN problem whose subproblems can be solved by the semi-proximal augmented Lagrangian method. We show that the generated sequence converges to a stationary point of the corresponding DC problem, which is feasible to the rank constrained DNN problem under some suitable assumptions. Moreover, numerical experiments demonstrate that for most QAP instances, the proposed approach can find the global optimal solutions efficiently, and for others, the proposed algorithm is able to provide good feasible solutions in a reasonable time.

The calibration of volatility models from observable option prices is a fundamental problem in quantitative finance. The most common approach among industry practitioners is based on the celebrated Dupire's formula [6], which requires the knowledge of vanilla option prices for a continuum of strikes and maturities that can only be obtained via some form of price interpolation. In this paper, we propose a new local volatility calibration technique using the theory of optimal transport. We formulate a time continuous martingale optimal transport problem, which seeks a martingale diffusion process that matches the known densities of an asset price at two different dates, while minimizing a chosen cost function. Inspired by the seminal work of Benamou and Brenier [1], we formulate the problem as a convex optimization problem, derive its dual formulation, and solve it numerically via an augmented Lagrangian method and the alternative direction method of multipliers (ADMM) algorithm. The solution effectively reconstructs the dynamic of the asset price between the two dates by recovering the optimal local volatility function, without requiring any time interpolation of the option prices.

—GeoFEM (Iizuka et al., 1999) is a parallel finite element analysis system intended for multi-physics, multi-scale problems of solid earth field phenomena. Very large linear elastic problems have already been resolved by parallel computation with GeoFEM. The next stage is to examine large-scale nonlinear problems using GeoFEM. The analysis of large-scale contact problems for fault zones is particularly important in the development of models that simulate the occurrence and cycle of earthquakes.¶This paper proposes a parallel FEM using an iterative solver with the augmented Lagrange method for solving large-scale contact problems. Direct solvers are presently applied to contact problem analysis because the matrix for such problems is ill-conditioned. However, direct solvers are not suitable for large-scale matrices. The augmented Lagrange method can improve the matrix conditions. The present study evaluates a parallel FEM using an iterative solver with the augmented Lagrange method. Analysis of a contact problem with the augmented Lagrange method revealed that an optimal penalty parameter exists and that large-scale parallel contact analysis using the iterative solver with localized preconditioning is promising.

In this paper, we denoise a given noisy image by minimizing a smoothness-promoting function over a set of local similarity measures which compare the mean of the given image and some candidate image on a large collection of subboxes. The associated convex optimization problem possesses a huge number of constraints which are induced by extended real-valued functions stemming from the Kullback–Leibler divergence. Alternatively, these nonlinear constraints can be reformulated as affine ones, which makes the model seemingly more tractable. For the numerical treatment of both formulations of the model (i.e., the original one as well as the one with affine constraints), we propose a rather general augmented Lagrangian method which is capable of handling the huge amount of constraints. A self-contained, derivative-free, global convergence theory is provided, allowing an extension to other problem classes. For the solution of the resulting subproblems in the setting of our suggested image denoising models, we make use of a suitable stochastic gradient method. Results of several numerical experiments are presented in order to compare both formulations and the associated augmented Lagrangian methods.

This paper presents a new point-to-face contact algorithm for contacts between two polyhedrons with planar boundaries. A new discrete numerical method called three-dimensional discontinuous deformation analysis (3-D DDA) is used and formulations of normal contact submatrices based on the proposed algorithm are derived. The presented algorithm is a simple and efficient method and it can be easily coded into a computer program. This approach does not need to use an iterative algorithm in each time step to obtain the contact plane, unlike the ‘Common-Plane’ method applied in the existing 3-D DDA. In the present 3-D DDA method, block contact constraints are enforced using the penalty method. This approach is quite simple, but may lead to inaccuracies that may be large for small values of the penalty number. The penalty method also creates block contact overlap, which violates the physical constraints of the problem. These limitations are overcome by using the augmented Lagrangian method that is used for normal contacts in this research. This point-to-face contact model has been programmed and some illustrative examples are provided to demonstrate the new contact rule between two blocks. A comparison between results obtained by using the augmented Lagrangian method and the penalty method is presented as well.

The augmented Lagrangian method (ALM) is a classical optimization tool that solves a given “difficult” (constrained) problem via finding solutions of a sequence of “easier” (often unconstrained) subproblems with respect to the original (primal) variable, wherein constraints satisfaction is controlled via the so-called dual variables. ALM is highly flexible with respect to how primal subproblems can be solved, giving rise to a plethora of different primal–dual methods. The powerful ALM mechanism has recently proved to be very successful in various large-scale and distributed applications. In addition, several significant advances have appeared, primarily on precise complexity results with respect to computational and communication costs in the presence of inexact updates and design and analysis of novel optimal methods for distributed consensus optimization. We provide a tutorial-style introduction to ALM and its variants for solving convex optimization problems in large-scale and distributed settings. We describe control-theoretic tools for the algorithms’ analysis and design, survey recent results, and provide novel insights into the context of two emerging applications: federated learning and distributed energy trading.

We present a local convergence analysis of the method of multipliers for equality-constrained variational problems (in the special case of optimization, also called the augmented Lagrangian method) under the sole assumption that the dual starting point is close to a noncritical Lagrange multiplier (which is weaker than second-order sufficiency). Local $$Q$$ Q -superlinear convergence is established under the appropriate control of the penalty parameter values. For optimization problems, we demonstrate in addition local $$Q$$ Q -linear convergence for sufficiently large fixed penalty parameters. Both exact and inexact versions of the method are considered. Contributions with respect to previous state-of-the-art analyses for equality-constrained problems consist in the extension to the variational setting, in using the weaker noncriticality assumption instead of the usual second-order sufficient optimality condition (SOSC), and in relaxing the smoothness requirements on the problem data. In the context of optimization problems, this gives the first local convergence results for the augmented Lagrangian method under the assumptions that do not include any constraint qualifications and are weaker than the SOSC. We also show that the analysis under the noncriticality assumption cannot be extended to the case with inequality constraints, unless the strict complementarity condition is added (this, however, still gives a new result).

We present an Augmented Lagrangian Method (ALM) for solving image reconstruction problems with a cost function consisting of multiple regularization functions with a data fidelity constraint. The presented technique is used to solve inverse problems related to image reconstruction, including compressed sensing formulations. Our contributions include an improvement for reducing the number of computations required by an existing ALM method, an approach for obtaining the proximal mapping associated with p-norm based regularizers, and lastly a particular ALM for the constrained image reconstruction problem with a hybrid cost function including a weighted sum of the p-norm and the total variation of the image. We present examples from Synthetic Aperture Radar imaging and Computed Tomography.

In this paper, we propose an image segmentation model where an $L^1$ variant of the Euler's elastica energy is used as boundary regularization. An interesting feature of this model lies in its preference for convex segmentation contours. However, due to the high order and non-differentiability of Euler's elastica energy, it is nontrivial to minimize the associated functional. As in recent work on the ordinary $L^2$-Euler's elastica model in imaging, we propose using an augmented Lagrangian method to tackle the minimization problem. Specifically, we design a novel augmented Lagrangian functional that deals with the mean curvature term differently than in previous works. The new treatment reduces the number of Lagrange multipliers employed, and more importantly, it helps represent the curvature more effectively and faithfully. Numerical experiments validate the efficiency of the proposed augmented Lagrangian method and also demonstrate new features of this particular segmentation model, such as shape driven and data driven properties.