A survey of model reduction methods for large-scale systems

An overview of model reduction methods and a comparison of the resulting algorithms are presented. These approaches are divided into two broad categories, namely SVD based and moment matching based methods. It turns out that the approximation error in the former case behaves better globally in frequency while in the latter case the local behavior is better.

December 2000

Restricted 2-factors in Bipartite Graphs

The k-restricted 2-factor problem is that of finding a spanning subgraph consisting of disjoint cycles with no cycle of length less than or equal to k. It is a generalization of the well known Hamilton cycle problem and is equivalent to this problem when n/2<=k<=n-1. This paper considers necessary and sufficient conditions, algorithms, and polyhedral conditions for 2-factors in bipartite graphs and restricted 2-factors in bipartite graphs. We introduce a generalization of the necessary and sufficient condition for 4-restricted 2-factors in bipartite graphs to 2k-restricted 2-factors in bipartite graphs of a particular form.

October 2000

Optimization for Parameter Estimation with Application to Transmission Electron Microscopy

We consider a parameter estimation problem for an important model in structural molecular biology. We propose two new mathematical formulations for the problem as constrained nonlinear least-squares problems, develop a numerical algorithm for solving this problem using interior-point methodology, and prove the convergence results for nonlinear least-squares problems with general constraints. Through numerical experimentation we show that our approach to the parameter determination problem is more effective than previous methods.

July 2000

Surrogate Modelling and Space Mapping for Engineering Optimization: A Summary of the Danish Technical University November 2000 Workshop

This is intended to be an outline of the topics presented in the November 16-18, 2000 workshop held at the Danish Technical University in Lyngby outside Copenhagen. The focus of this workshop was on the construction of inexpensive models and their use with optimization to support engineering decision making. The purpose of this summary is to organize the various topics presented there into a unified whole. An additional purpose is to give an intuitive introduction to the promising concept of space mapping as a way to leverage a cheap coarse model into an effective optimization surrogate for a more detailed model.

Key words: response surface methodology, kriging, DACE, surrogate-based optimization, space mapping, pattern search methods.

December 2000

Maximum Stable Set Formulations and Heuristics Based on Continuous Optimization

The stability number for a given graph G is the size of a maximum stable set in G. The Lovasz theta number provides an upper bound on the stability number and can be computed as the optimal value of the Lovasz semidefinite program. In this paper, we show that restricting the matrix variable in the Lovasz semidefinite program to be rank-one or rank-two yields a pair of continuous, nonlinear optimization problems each having the global optimal value equal to the stability number. We propose heuristics for obtaining large stable sets in G based on these new formulations and present computational results indicating the effectiveness of the heuristics.

December 2000

Rank-Two Relaxation Heuristics for Max-Cut and Other Binary Quadratic Programs

Semidefinite relaxation for certain discrete optimization problems involves replacing a vector-valued variable by a matrix-valued one, producing a convex program while increasing the number of variables by an order of magnitude. As useful as it is in theory, this approach encounters difficulty in practice as problem size increases. In this paper, we propose a rank-two relaxation approach and construct continuous optimization heuristics applicable to some binary quadratic programs, including primarily the Max-Cut problem but also others such as the Max-Bisection problem. A computer code based on our rank-two relaxation heuristics is compared with two state-of-the-art semidefinite programming codes. We will report some rather intriguing computational results on a large set of test problems and discuss their ramifications.

November 2000

A Fast, Fully Implicit Backward Euler Solver for Dendritic Neurons

We develop and test a C++ implementation of a discretization of the Hodgkin-Huxley equations for dendritic neurons which employs backward Euler in time and finite differences in space. We make use of the sparse analytical Jacobian matrix to perform the nonlinear solve required at each time step via Newton's method.

September 2000

Time-Domain Decomposition Iterative Methods for the Solution of Distributed Linear Quadratic Optimal Control Problems

We study a class of time-domain decomposition based methods for the solution of distributed linear quadratic optimal control problems. Our methods are based on a multiple shooting reformulation of the distributed linear quadratic optimal control problem as a discrete-time optimal control (DTOC) problem in Hilbert space. The optimality conditions for this DTOC problem lead to a linear system with block structure. This motivates the application of block Gauss-Seidel methods for its solution. We show that certain instantaneous control techniques can be viewed as the application of one step of the forward block Gauss-Seidel method applied to the DTOC optimality system. To obtain better convergence properties, we imbed the block Gauss-Seidel methods as preconditioners in a Krylov-subspace method.

Key words: linear quadratic optimal control problems, instantaneous control, suboptimal control, multiple shooting, discrete-time optimal control problem, Gauss-Seidel method, Krylov subspace methods.

November 2000

Understanding code generated by TAMC

The Tangent linear and Adjoint Model Compiler (TAMC) is an automatic differentiation package designed and implemented by Ralf Giering. The code generated by TAMC can be understood in terms of a careful mathematical model for a subroutine implementing an operator.

June 2000

Identification of Regional Variation in the Constitutive Response of Axisymmetric Membranes

We demonstrate that the equilibrium equations for an axisymmetric, nonlinear, anisotropic membrane under hydrostatic pressure allow explicit representation of the longitudinal and azimuthal stresses in terms of the associated longitudinal and azimuthal strains. We apply this result in a numerical simulation of the canine diaphragm. More precisely, we compute the deformation of the membrane under a quasi static increase in the intensity of the applied hydrostatic load. The associated strains are easily estimated via finite differences. As the membrane is inflated the set of achieved strains grows and as a result of our explicit representation formula, we recover larger and larger patches of the associated stress surfaces.

Key words: material identification, constitutive response, membrane, diaphragm.

September 2000

Correcting an Inconsistent System of Linear Inequalities by Nonlinear Programming

We consider the problem of correcting an inconsistent system of linear inequalities, Ax <= b, subject to nonnegativity constraints, x >= 0. We formulate this problem as a nonlinear program and derive the corresponding Karush-Kuhn-Tucker conditions. These conditions reveal several interesting properties that solutions must satisfy and allow us to derive several equivalent problems that involve fewer decision variables and are more amenable to solution. We propose using a gradient projection method to minimize an objective function Ø(x) subject only to x >= 0. We also propose a hybrid approach that exploits an interesting relation between the correction problem and the method of total least squares.

July 2000

Simultaneous Structure Factor and Contrast Transfer Function Parameter Determination in Transmission Electron Microscopy

We present a new method that allows a fully automated simultaneous determination of the structure factor and the parameters of the Contrast Transfer Function (CTF) and noise function. No previous knowledge of the structure factor or of the CTF parameters is assumed. Our approach is based on the new precise mathematical formulations of the problem as constrained nonlinear least squares that treats the structure factor as a set of undetermined variables, as well as the CTF parameters, and on the interior-point algorithm that satisfies the inequality constraints on the bounds.

Key words: electron microscopy, structure factor, contrast transfer function, CTF, nonlinear least-squares problem.

August 2000

Designing and Analyzing Computational Experiments for Global Optimization

We consider a variety of issues that arise when designing and analyzing computational experiments for global optimization. We describe a probability model for objective functions and a method for generating pseudorandom objective functions. We argue in favor of evaluating the performance of global optimization algorithms by measuring the depth of the objective function achieved with a fixed number of function evaluations. We emphasize the importance of replication in computational experiments and describe some useful statistical techniques for assimilating results. We illustrate our methods by performing a small study that compares two multistart strategies for global optimization.

July 2000

Mathematical Modeling of Protein Structure Using Distance Geometry

This paper reviews methods for structure determination with interatomic distances and explores possible improvement of the methods and ways of combining them with potential energy minimization.

Key words: distance geometry, potential energy minimization, nuclear magnetic resonance (NMR) spectroscopy, singular value decomposition (SVD), graph embedding, global optimization.

July 2000

On Eliminating Square Paths in a Square Lattice

Removing the minimum number of vertices or points from a square lattice such that no square path exists is known as the square path problem. Finding this number as the size of the lattice increases is not so trivial. Results provided by Erdös-Pósa and Bienstock-Dean provides an upper bound for eliminating all cycles from a planar graph but sheds little light on the case of the square lattice. This paper provides several values for the minimum number of vertices needed to be removed such that no square path exists.

April 2000

SMV--An Object-Oriented Sparse Matrix-Vector Computation Library: Programming Manual

This is a technical document for a matrix-vector class library we have developed recently to support dense and sparse matrix-vector computation in optimization applications. The library is written in C++ and runs on single processor platforms. A parallel version is planed for future development.

June 2000

Mixed Variable Optimization of the Number and Composition of Heat Intercepts in a Thermal Insulation System

In the literature, thermal insulation systems with a fixed number of heat intercepts have been optimized with respect to intercept locations and temperatures. The number of intercepts and the types of insulators that surround them were chosen by parametric studies. This was because the optimization methods used could not treat such categorical variables. Discrete optimization variables are categorical if the objective function or the constraints can not be evaluated unless the variables take one of a prescribed enumerable set of values. The key issue is that categorical variables can not be treated as ordinary discrete variables are treated by relaxing them to continuous variables with a side constraint that they be discrete at the solution.

A new mixed variable programming (MVP) algorithm makes it possible to optimize directly with respect to mixtures of discrete, continuous, and categorical decision variables. The result of applying MVP is shown here to give a 65% reduction in the objective function over the previously published result for a thermal insulation model from the engineering literature. This reduction is largely because MVP optimizes simultaneously with respect to the number of heat intercepts and the choices from a list of insulator types as well as intercept locations and temperatures. The main purpose of this paper is to show that the mixed variable optimization algorithm can be applied effectively to a broad class of optimization problems in engineering that could not be easily solved with earlier methods.

Key words: optimization, thermal insulation, heat intercepts, categorical variables, mixed variable programming (MVP), pattern search algorithm.

June 2000 (revised May 2001)

On the Use of Direct Search Methods for Stochastic Optimization

We examine the conventional wisdom that commends the use of directe search methods in the presence of random noise. To do so, we introduce new formulations of stochastic optimization and direct search. These formulations suggest a natural strategy for constructing globally convergent direct search algorithms for stochastic optimization by controlling the error rates of the ordering decisions on which direct search depends. This strategy is successfully applied to the class of generalized pattern search methods. However, a great deal of sampling is required to guarantee convergence with probability one.

June 2000

The Inverse Problem of Neuron Identification

Depending on the state of neuron membrane, the inverse problem of neuron identification is divided into two categories: the passive neuron identification and the active neuron identification. In the first category, we provided a more efficient way to recover neuron parameters than the traditional approach. By exploring the impedance function meticulously, our method reveals a clean and analytical relation between the electrical properties of neurons and their response to sub-threshold current stimulation.

Mathematical equations like the Hodgkin-Huxley equations and the Fitzhugh-Nagumo equations that model active neurons have been established for many years.  However, the inverse problem in this category has barely started.  Our research in this direction attempts to establish a proper formulation of the inverse problem and to investigate possible mathematical techniques that are needed to solve it. For the relatively simple Fitzhugh-Nagumo equations, we successfully reconstructed the nonlinear membrane conductance function and the coefficients of the recovering variable. The method is then extended to a more realistic neuron model, the Morris-Lecar model. We provide a computational strategy for systematically recovering the nonlinearity of both calcium and the potassium channels.

May 2000

The Effect of the Separation of Variables on the Molecular Replacement Method

Traditional approaches for solving the molecular replacement problem separate a six-dimensional optimization problem into two three-dimensional ones in order to reduce the computational cost. There are, however, serious drawbacks in such a separation of the rotational and translational degrees of freedom. In this paper, we present computational experiments indicating that even under ideal conditions the separation can fail to preserve the correspondence between the global minima of a target function and the correct rotations when low resolution data are used. This phenomenon is a reason why only high resolution data are used in traditional approaches for solving the molecular replacement problem. In this paper, we provide a theoretical explanation for this phenomenon.

In order to solve difficult molecular replacement problems, we believe that low resolution terms should be utilized because they generate smooth, shape-defining components in a target function, making it more amenable to global optimization. This study indicates that in order to utilize low resolution data in the molecular replacement method, we need to consider all degrees of freedom simultaneously. The full-dimensional optimization formulation, once a prohibitive procedure due to its high computational cost, should now be feasible given the current state of computational resources and algorithms.

Key words: global minimization, continuation method.

June 2000

Branch Decompositions and their Applications

Many real-life problems can be modeled as optimization or decision problems on graphs. Also, many of those real-life problems are NP-hard. One traditional method to solve these problems is by branch and bound while another method is by graph decompositions. In the 1980's, Robertson and Seymour conceived of two new ways to decompose the graph in order to solve these problems. These ingenious ideas were only byproducts of their work proving Wagner's Conjecture. A branch decomposition is one of these ideas. A paper by Arnborg, Lagergren and Seese showed that many NP-complete problems can be solved in polynomial time using divide and conquer techniques on input graphs with bounded branchwidth, but a paper by Seymour and Thomas proved that computing an optimal branch decomposition is also NP-complete. Although computing optimal branch decompositions is NP-complete, there is a plethora of theory about branchwidth and branch decompositions. For example, a paper by Seymour and Thomas offered a polynomial time algorithm to compute the branchwidth and optimal branch decomposition for planar graphs.

This doctoral research is concentrated on constructing branch decompositions for graphs and using branch decompositions to solve NP-complete problems modeled on graphs. In particular, a heuristic to compute near-optimal branch decompositions is presented and the heuristic is compared to previous heuristics in the subject. Furthermore, a practical implementation of an algorithm given in a paper by Seymour and Thomas for computing optimal branch decompositions of planar graphs is implemented with the addition of heuristics to give the algorithm a "divide and conquer" design. In addition, this work includes a theoretical result relating the branchwidth of planar graphs to their duals, characterizations of branchwidth for Halin and chordal graphs. Also, this work presents an algorithm for minor containment using a branch decomposition and a parallel implementation of the heuristic for general graphs using pthreads.

April 2000

Joint Inversion Using the Convolutional Model

A detailed imaging of an acoustic medium via a seismic experiment requires an accurate representation of the source. A joint source and reflectivity inversion may provide the means to obtain the desired detail. Joint inversion using the acoustic wave equation is computationally expensive. The convolutional model of the seismogram in the offset-time domain provides a computationally cheaper means of approximating the wave equation.

We rigorously derive the convolutional model of the seismogram as an approximation to the linearized acoustic wave equation where the medium to be imaged is layered and has a constant density. Using the Hilbert Class Library (a mathematical C++ library that among other things provides access to optimization algorithms), we implement the convolutional model and its inverse. Using linear and non-linear optimization methods applied to a least squares data residual objective function we test the robustness of inversion results from the convolutional model in the offset-time domain.

April 2000

Imaging Complex Structures With Semi-recursive Kirchhoff Migration

Attempting to image the subsurface in areas of complex geology and rapid lateral velocity variation is a challenging problem. In particular, using prestack Kirchhoff migration in conjunction with first arrival travel times produces poor subsurface images with increasing depth. This problem is not a limitation of Kirchhoff migration, but it is the failure of the finite difference method to compute travel times which correspond to the most energetic arrivals. Dimitri Bevc in 1997 proposed a technique that combines the use of the wave equation datuming with dividing the velocity model into subsets. In each subset, calculation of the travel times with the finite differencing eikonal equation is valid. This thesis applies Bevc's technique using a software package (ProMAX) that is widely used among the academic and the industrial communities. We not only get a superior image at depth but we also enjoy the simplicity and the computational efficiency of using the finite difference method. We use the Marmousi synthetic dataset as input, which satisfies the definitions of structural complexity and rapid lateral velocity variation. To demonstrate the effectiveness of our approach, tests were performed on the Marmousi dataset before and after the application of the semi-recursive Kirchhoff migration.

May 2000

Geometrical Optics for Quasi-P Waves: Theories and Numerical Methods

The quasi-P wave in anisotropic solids is of practical importance in obtaining maximal imaging resolution in seismic exploration. The geometrical optics term in the asymptotic expansion for the wave characterizes the high frequency part of the quasi-P wave by using two functions:& a phase (traveltime) function satisfying an eikonal equation and an amplitude function satisfying a transport equation.

I develop theories and numerical methods for constructing the geometrical optics term of quasi-P waves in general anisotropic solids. The traveltime corresponding to the downgoing wave satisfies a paraxial eikonal equation, an evolution equation in depth. This paraxial eikonal equation takes into account the convexity of the quasi-P slowness surface and thus has a built-in reliable indicator of the ray velocity direction. Therefore, high-order finite-difference eikonal solvers are easily constructed by utilizing Weighted Essentially NonOscillating (WENO) schemes.

Because the amplitude function is related to second-order derivatives of the traveltime, a third-order accurate eikonal solver for traveltimes is necessary to get a first-order accurate amplitude. However, the eikonal equation with a point source has an upwind singularity at the source which renders all finite-difference eikonal solvers to be first-order accurate near the source. A new approach combining an adaptive-gridding strategy with WENO schemes can treat this singularity efficiently and can yield highly accurate traveltimes and amplitudes for both isotropic and anisotropic solids.

A variety of numerical experiments verify that the new paraxial eikonal solver and adaptive-gridding-WENO approach are accurate and efficient for capturing the anisotropy. Therefore, the two new methods provide tools for constructing the geometrical optics term of the quasi-P wave in general anisotropic solids.

April 2000

Lyapunov, Lanczos, and Inertia

We present a new proof of the inertia result associated with Lyapunov equations. Furthermore we present a connection between the Lyapunov equation and the Lanczos process which is closely related to the Schwarz form of a matrix. We provide a method for reducing a general matrix to Schwarz form in a finite number of steps (O(n³)). Hence, we provide a finite method for computing inertia without computing eigenvalues. This scheme is unstable numerically and hence is primarily of theoretical interest.

Key words: Lanczos method, Lyapunov equation, inertia, stability, control.

May 2000

A Globally Convergent Primal-Dual Interior-Point Filter Method for Nonconvex Nonlinear Programming

In this paper, the filter technique of Fletcher and Leyffer (1997) is used to globalize the primal-dual interior-point algorithm for nonlinear programming, avoiding the use of merit functions and the updating of penalty parameters.

The new algorithm decomposes the primal-dual step obtained from the perturbed first-order necessary conditions into a normal and a tangential step, whose sizes are controlled by a trust-region type parameter. Each entry in the filter is a pair of coordinates: one resulting from feasibility and centrality, and associated with the normal step; the other resulting from optimality (complementarity and duality), and related with the tangential step.

Global convergence to first-order critical points is proved for the new primal-dual interior-point filter algorithm.

Key words: interior-point methods, primal-dual, filter, global convergence.

April 2000

Semismooth Newton Methods for Operator Equations in Function Spaces

We develop a semismoothness concept for nonsmooth superposition operators in function spaces. The considered class of operators includes NCP-function-based reformulations of infinite-dimensional nonlinear complementarity problems, and thus covers a very comprehensive class of applications. Our results generalize semismoothness and alpha-order semismoothness from finite-dimensional spaces to a Banach space setting. Hereby, a new generalized differential is used that can be seen as an extension of Qi's finite-dimensional C-subdifferential to our infinite-dimensional framework. We apply these semismoothness results to develop a Newton-like method for nonsmooth operator equations and prove its local q-superlinear convergence to regular solutions. If the underlying operator is alpha-order semismoothness, convergence of q-order 1+alpha is proved. We also establish the semismoothness of composite operators and develop corresponding chain rules. The developed theory is accompanied by illustrating examples and by applications to nonlinear complementarity problems.

Key words: Newton-like methods, semismoothness, superposition operators, generalized differentials, nonlinear complementarity problems, superlinear convergence.

April 2000

A Sensitivity and Adjoint Calculus for Discontinuous Solutions of Hyperbolic Conservation Laws with Source Terms

We present a sensitivity and adjoint calculus for the control of entropy solutions of scalar conservation laws with controlled initial data and source term.  The sensitivity analysis is based on shift-variations which are the sum of a standard variation and suitable corrections by weighted indicator functions approximating the movement of the shock locations. Based on a first order approximation by shift-variations in L¹ we introduce the concept of shift-differentiability which is applicable to operators having functions with moving discontinuities as images and implies differentiability for a large class of tracking-type functionals. In the main part of the paper we show that entropy solutions are generically shift-differentiable at almost all times t > 0 with respect to the control. Hereby we admit shift-variations for the initial data which allows to use the shift-differentiability result repeatedly over time slabs. This is useful for the design of optimization methods with time domain decomposition. Our analysis, especially of the shock sensitivity, combines structural results by using generalized characteristics and an adjoint argument. Our adjoint based shock sensitivity analysis does not require to restrict the richness of the shock structure a priori and admits shock generation points. The analysis is based on stability results for the adjoint transport equation with discontinuous coefficients satisfying a one-sided Lipschitz condition. As a further main result we derive and justify an adjoint representation for the derivative of a large class of tracking-type functionals.

Key words: optimal control, scalar conservation law, shock sensitivity, adjoint state, Fréchet differentiability.

March 2000

A Pattern Search Filter Method for Nonlinear Programming without Derivatives

This paper formulates and analyzes a pattern search method for general constrained optimization based on filter methods for step acceptance. Roughly, a filter method accepts a step that either improves the objective function value or the value of some function that measures the constraint violation. The new algorithm does not compute or approximate any derivatives, penalty constants or Lagrange multipliers. A key feature of the new algorithm is that it preserves the useful division into global SEARCH and a local POLL steps. It is shown here that the algorithm identifies limit points at which optimality conditions depend on local smoothness of the functions. Stronger optimality conditions are guaranteed for smoother functions. In the absence of general constraints, the proposed algorithm and its convergence analysis generalizes the previous work on unconstrained, bound constrained and linearly constrained generalized pattern search. The algorithm is illustrated on some test examples and on an industrial wing planform engineering design application.

Key words: pattern search algorithm, filter algorithm, surrogate-based optimization, derivative-free convergence analysis, constrained optimization, nonlinear programming.

March 2000 (revised November 2003)

What is Simulated Annealing?

Beginning in 1983, simulated annealing was marketed as a global optimization methodology that mimics the physical annealing process by which molten substances cool to crystalline lattices of minimal energy. This marketing strategy had a polarizing effect, attracting those who delighted in metaphor and alienating others who found metaphor insufficient at best and facile at worst. In fact, the emotional outbursts that accompany many discussions of simulated annealing are an unfortunate distraction. Whatever its pros and cons, simulated annealing can be grounded in rigorous mathematics. Here we provide an elementary, self-contained introduction to simulated annealing in terms of Markov chains.

February 2000

Analysis of Generalized Pattern Searches

This paper contains a new convergence analysis for the Lewis and Torczon generalized pattern seach (GPS) class of methods for unconstrained and linearly constrained optimization. This analysis is motivated by a desire to understand the successful behavior of the algorithm under hypotheses that are satisfied by many practical problems. Specifically, even if the objective function is discontinuous or extended valued, the methods find a limit point with some minimizing properties. Simple examples show that the strength of the optimality conditions at a limit point does not depend only on the algorithm, but also on the directions it uses, and on the smoothness of the objective at the limit point in question. This contribution of this paper is to provide a simple convergence analysis that supplies detail about the relation of optimality conditions to objective smoothness properties and to the defining directions for the algorithm, and it gives previous results as corollaries.

Key words: Pattern search algorithm, linearly constrained optimization, surrogate-based optimization, nonsmooth optimization, derivative-free convergence analysis.

February 2000 (revised June 2002)

A Study of the Stationary Configurations of the SStress Criterion for Metric Multidimensional Scaling

It is widely believed that both the stress and the sstress criteria for metric multidimensional scaling are plagued by the existence of nonglobal minimizers. At present, there is little theory that enlightens this belief. Trosset and Mathar (1997) established that nonglobal minimizers of the stress criterion can exist, while Glunt, Hayden, and Liu (1991) demonstrated that the distance matrices of all configurations for which the gradient of the sstress criterion vanishes lie on a certain sphere. This report extends existing theory in several directions. Emphasis is placed on the more tractable case of the sstress criterion. Because the stress and sstress criteria must be minimized by numerical optimization, one result that is of immediate practical value is a simple device for improving the quality of the initial configurations from which optimization commences.

January 2000

Parallel Continuous Optimization

Parallel continuous optimization methods are motivated here by applications in science and engineering. The key issues are addressed at different computational levels including local and global optimization as well as strategies for large, sparse versus small but expensive problems. Topics covered include global optimization, direct search with and without surrogates, optimization of linked subsystems, and variable and constraint distribution. Finally, there is a discussion of future research directions.

January 2000 (revised May 2000)

Go to 1999 reports

Go to 2001 reports