**Historical Development of the BFGS Secant Method and Its Characterization Properties**

Joanna Maria Papakonstantinou

The BFGS secant method is the preferred secant method for finite-dimensional unconstrained optimization. The first part of this research consists of recounting the historical development of secant methods in general and the BFGS secant method in particular. Many people believe that the secant method arose from Newton's method using finite difference approximations to the derivative. We compile historical evidence revealing that a special case of the secant method predated Newton's method by more than 3000 years. We trace the evolution of secant methods from 18th-century B.C. Babylonian clay tablets and the Egyptian Rhind Papyrus. Modifications to Newton's method yielding secant methods are discussed and methods we believe influenced and led to the construction of the BFGS secant method are explored.

In the second part of our research, we examine the construction of several rank-two secant update classes that had not received much recognition in the literature. Our study of the underlying mathematical principles and characterizations inherent in the updates classes led to theorems and their proofs concerning secant updates. One class of symmetric rank-two updates that we investigate is the Dennis class. We demonstrate how it can be derived from the general rank-one update formula in a purely algebraic manner not utilizing Powell's method of iterated projections as Dennis did it. The literature abounds with update classes; we show how some are related and show containment when possible. We derive the general formula that could be used to represent all symmetric rank-two secant updates. From this, particular parameter choices yielding well-known updates and update classes are presented. We include two derivations of the Davidon class and prove that it is a maximal class. We detail known characterization properties of the BFGS secant method and describe new characterizations of several secant update classes known to contain the BFGS update. Included is a formal proof of the conjecture made by Schnabel in his 1977 Ph.D. thesis that the BFGS update is in some asymptotic sense the average of the DFP update and the Greenstadt update.

November 2010

**Alternating Direction Augmented Lagrangian Methods for Semidefinite Programming**

Zaiwen Wen, Donald Goldfarb, & Wotao Yin

We present an alternating direction method based on an augmented Lagrangian framework for solving semidefinite programming (SDP) problems in standard form. At each iteration, the algorithm, also known as a two-splitting scheme, minimizes the dual augmented Lagrangian function sequentially with respect to the Lagrange multipliers corresponding to the linear constraints, then the dual slack variables and finally the primal variables, while in each minimization keeping the other variables fixed. Convergence is proved by using a fixed-point argument. A multiple-splitting algorithm is then proposed to handle SDPs with inequality constraints and positivity constraints directly without transforming them to the equality constraints in standard form. Finally, numerical results for frequency assignment, maximum stable set and binary integer quadratic programming problems are presented to demonstrate the robustness and efficiency of our algorithm.

December 2009

**TR09-41**
PDF File

**Collaborative Spectrum Sensing from Sparse Observations for Cognitive Radio Networks**

Jia Meng, Wotao Yin, Husheng Li, Ekram Houssian, & Zhu Han

Spectrum sensing, which aims at detecting spectrum holes, is the precondition for the implementation of cognitive radio. Collaborative spectrum sensing among the cognitive radio nodes is expected to improve the ability of checking complete spectrum usage. Due to hardware limitations, each cognitive radio node can only sense a relatively narrow band. Consequently, the available channel sensing information is far from being sufficient for recognizing the wide range of unoccupied channels precisely. Aiming at breaking this bottleneck, we propose to apply matrix completion and joint sparsity recovery to reduce sensing and transmitting requirements and improve sensing results. Specifically, equipped with a frequency selective filter, each cognitive radio node senses linear combinations of multiple channel information and reports them the fusion center, where occupied channels are then decoded from the reports by novel matrix completion and joint sparsity recovery algorithms. As a result, the number of reports sent from the CRs to the fusion center are significantly reduced. We propose two decoding approaches, one based on matrix completion and the other based on joint sparsity recovery, both of which allow exact recovery from incomplete reports. The numerical results validate the effectiveness and robustness of our approaches. In particular, in small scale networks, the matrix completion approach achieves exact channel detection with a number of samples no more than 50% of the number of channels in the network; While joint sparsity recovery achieves similar performance in large-scale networks.

December 2009

**Copula Density Estimation by Total Variation Penalized Likelihood with Linear Equality Constraints**

Leming Qu & Wotao Yin

A copula density is the joint probability density function (PDF) of a random vector with uniform marginals. An approach to bivariate copula density estimation is introduced that is based on a maximum penalized likelihood estimation (MPLE) with a total variation (TV) penalty term. The marginal unity and symmetry constraints for copula density are enforced by linear equality constraints. The TV-MPLE subject to linear equality constraints is solved by an augmented Lagrangian and operator-splitting algorithm. It offers an order of magnitude improvement in computational efficiency over another TV-MPLE method without constraints solved by log-barrier method for second order cone program. A data-driven selection of the regularization parameter is through K-fold cross-validation (CV). Simulation and real data application show the effectiveness of the proposed approach. The MATLAB code implementing the methodology is available online.

December 2009

**Collaborative Spectrum Sensing from Sparse Observations Using Matrix Completion**

Jia Meng, Wotao Yin, Husheng Li, Ekram Houssian, & Zhu Han

In cognitive radio, spectrum sensing is a key component to detect spectrum holes (i.e., channels not used by any primary users). Collaborative spectrum sensing among the cognitive radio nodes is expected to improve the ability of checking complete spectrum usage states. Unfortunately, due to power limitation and channel fading, available channel sensing information is far from being sufficient to tell the unoccupied channels directly. Aiming at breaking this bottleneck, we apply recent matrix completion techniques to greatly reduce the sensing information needed. We formulate the collaborative sensing problem as a matrix completion subproblem and a joint-sparsity reconstruction subproblem. Results of numerical simulations that validated the effectiveness and robustness of the proposed approach are presented. In particular, in noiseless cases, when number of primary user is small, exact detection was obtained with no more than 8% of the complete sensing information, whilst as number of primary user increases, to achieve a detection rate of 95.55%, the required information percentage was merely 16.8%.

December 2009

**Analysis of Weak Solutions For the Fully Coupled Stokes-Darcy-Transport Problem**

Aycil Cesmelioglu & Béatrice M. Rivière

This paper analyzes the surface/subsurface flow coupled with transport. The flow is modeled by the coupling of Stokes and Darcy equations. The transport is modeled by a convection-dominated parabolic equation. The two-way coupling between flow and transport is nonlinear and it is done via the velocity field and the viscosity. This problem arises from a variety of natural phenomena such as the contamination of the groundwater through rivers. The main result is existence and stability bounds of a weak solution.

November 2009

**Alternating Direction Algorithms for L1-Problems in Compressive Sensing**

Junfeng Yang & Yin Zhang

In this paper, we propose and study the use of alternating direction algorithms for several L1-norm minimization problems arising from sparse solution recovery in compressive sensing, including the basis pursuit problem, the basis-pursuit denoising problems of both unconstrained and constrained forms, as well as others. We present and investigate two classes of algorithms derived from either the primal or the dual forms of the L1-problems. The construction of the algorithms consists of two main steps: (1) to reformulate an L1-problem into one having partially separable objective functions by adding new variables and constraints; and (2) to apply an exact or inexact alternating direction method to the resulting problem. The derived alternating direction algorithms can be regarded as first-order primal-dual algorithms because both primal and dual variables are updated at each and every iteration. Convergence properties of these algorithms are established or restated when they already exist. Extensive numerical results in comparison with several state-of-the-art algorithms are given to demonstrate that the proposed algorithms are efficient, stable and robust. Moreover, we present numerical results to emphasize two practically important but perhaps overlooked points. One point is that algorithm speed should always be evaluated relative to appropriate solution accuracy; another is that whenever erroneous measurements possibly exist, the l1-norm fidelity should be the fidelity of choice in compressive sensing.

November 2009

**Dimension Reduction for Unsteady Nonlinear Partial Differential Equations via Empirical Interpolation Methods**

Saifon Chaturantabut

This thesis evaluates and compares the efficiencies of techniques for constructing reduced-order models for finite difference (FD) and finite element (FE) discretized systems of unsteady nonlinear partial differential equations (PDEs). With nonlinearity, the complexity for solving the reduced-order system constructed directly from the well-known Proper Orthogonal Decomposition (POD) technique alone still depends on the dimension of the original system. *Empirical Interpolation Method (EIM)*, proposed in [2], and its discrete variation, *Discrete Empirical Interpolation Method (DEIM)*, introduced in this thesis, are therefore combined with the POD technique to remove this inefficiency in the nonlinear terms of FE and FD cases, respectively. Numerical examples demonstrate that both POD-EIM and POD-DEIM approaches not only dramatically reduce the dimension of the original system with high accuracy, but also remove the dependence on the dimension of the original system as reflected in the decrease computational time compared to the POD approach.

October 2009

**A Matlab Implementation of a Flat Norm Motivated Polygonal Edge Matching Method using a Decomposition of Boundary into Four 1-Dimensional Currents**

Simon Morgan, Wotao Yin, & Kevin Vixie

We describe and provide code and examples for a polygonal edge matching method.

September 2009

**Image-Based Face Illumination Transferring Using Logarithmic Total Variation Models**

Qing Li, Wotao Yin, & Zhigang Deng

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

__Key Words:__ BIllumination transferring, Face relighting, Logarithmic total variation model, Radial basis functions, Face decomposition

September 2009

**Time-Stepping Classes for Optimization (TSOpt)**

Marco Enriquez & William W. Symes

This report introduces the "Time Stepping Package for Optimization", or TSOpt, which is an interface for time-stepping simulation written in C++. It packages a simulator together with its derivatives (\sensitivities") and adjoint derivatives with respect to simulation parameters in a single object called a Jet, which can be used in conjunction with an optimization algorithm to solve a simulation-driven optimization problem. Further, TSOpt interfaces with the Rice Vector Library (RVL), allowing Jet objects to define a Operator subclass.

September 2009

**A C++ Class Supporting Adjoint-State Methods**

Marco Enriquez

The adjoint-state method is widely used for computing gradients in simulation- driven optimization problems. The adjoint-state evolution equation requires access to the entire history of the system states. There are instances, however, where the required state for the adjoint-state evolution is not readily accessible; consider large- scale problems, for example, where the entire simulation history is not saved to con- serve memory. This thesis introduces a C++ state-access class, StateHistory, to support a myriad of solutions to this problem. Derived StateHistory classes im- plement a (simulation) time-altering function and data-access functions, which can be used in tandem to access the entire state history. This thesis also presents a derived StateHistory class, GriewankStateHistory, which uses Griewank's opti- mal checkpointing scheme. While only storing a small fraction of simulation states, GriewankStateHistory objects can reconstitute unsaved states for a small computa- tional cost. These ideas were implemented in the context of TSOpt, a time-stepping library for simulation-driven optimization algorithms.

September 2009

**Numerical Solutions of Matrix Equations Arising in Model Reduction of Large-Scale Linear-Time-Invariant Systems**

Hung Dinh Nong

This thesis presents and analyzes new algorithms for matrix equations arising from model reduction of linear-time-invariant (LTI) systems. Such systems arise in a variety of areas, especially in circuit simulation. When an integrated circuit has millions of devices, performing a full-system simulation can be infeasible due to the time required for the nonlinear solver to compute the solution of large linearized systems. Model reduction becomes indispensable since model reduction techniques aim to produce low-dimensional systems that capture the same response characteristics as the originals while enabling substantial improvement in simulation time and resulting in greatly reduced storage requirements.

One of the crucial physical features of LTI systems in circuit simulation is passivity, which is essential to preserve during model reduction. Passivity preserving model reduction using the invariant subspace method in conjunction with positive real balancing for LTI systems provides an error bound and involves the solution of two Riccati equations. The contributions of this thesis are fourfold. First, a generalization of the invariant subspace method for passivity preserving model reduction is presented for LTI systems in descriptor form. Second, two iterative algorithms to compute the solution of the two algebraic Riccati equations are derived. Even though these two iterative algorithms have been in existence, the thesis introduces them here in addition with specific conditions on LTI systems for their convergence.

At each step of either iteration, a pair of Lyapunov equations result. The thesis next investigates a parameter free ADI-like (PFADI) method for the solution of the Lyapunov equation. This PFADI method involves the solution of the Sylvester equation at each step, and its computational complexity depends on how efficiently the solution of the Sylvester equation is obtained. The final part of the thesis is devoted to the study of two techniques to solve the Sylvester equation. The first technique is to obtain the solution of the Sylvester equation via that of an invariant subspace problem. The thesis presents a complete analysis and an efficient algorithm to compute the optimal real shift for the Cayley transformation used in the invariant subspace technique. The analysis is then generalized for multiple optimal real shift selection for the ADI technique for the solution of the Sylvester equation.

September 2009

**Sparse Signal Reconstruction via Iterative Support Detection **

Yilun Wang & Wotao Yin

We present a novel sparse signal reconstruction method "ISD", aiming to achieve fast reconstruction and a reduced requirement on the number of measurements compared to the classical L1 minimization approach. ISD addresses failed reconstructions of L1 minimization due to insufficient measurements. It estimates a support set *I* from a current reconstruction and obtains a new reconstruction by solving the minimization problem min{sum_{*i not in
I*} |*x_i*| : *Ax = b*}, and it iterates these two steps for a small number of times. ISD differs from the orthogonal matching pursuit (OMP) method, as well as its variants, because (i) the index set *I* in ISD is
not necessarily nested or increasing and (ii) the minimization problem above updates all the components of *x* at the same time. We generalize the Null Space Property to Truncated Null Space Property and present our analysis
of ISD based on the latter.

We introduce an efficient implementation of ISD, called threshold-ISD, for recovering signals with fast decaying distributions of nonzeros from compressive sensing measurements. Numerical experiments show that threshold-ISD has significant advantages over the classical L1 minimization approach, as well as two state-of-the-art algorithms: the iterative reweighted L1 minimization algorithm (IRL1) and the iterative reweighted least-squares algorithm (IRLS).

MATLAB code is available for download from: http://www.caam.rice.edu/~optimization/L1/ISD/.

September 2009; Revise June 2010

**Balanced Truncation Model Reduction for Systems with Inhomogeneous Initial Conditions**

M. Heinkenschloss, T. Reis, & A. C. Antoulas

We present a rigorous approach to extend balanced truncation model reduction (BTMR) to systems with inhomogeneous initial conditions, we provide estimates for the error between the input-output maps of the original and of the reduced initial value system, and we illustrate numerically the superiority of our approach over the naive application of BTMR. When BTMR is applied to linear time-invariant systems with inhomogeneous initial conditions, it is crucial that the initial data are well represented by the subspaces generated by BTMR. This requirement is often ignored or it is avoided by making the restrictive assumption that the initial data are zero. To ensure that the initial data are well represented by the BTMR subspaces, we add auxiliary inputs determined by the initial data. To obtain our error estimate, we approximate the contribution of the inhomogeneous initial data by a suitable L^{2} input function.

__Key Words:__ Balanced truncation, initial value problems, error bound, Hankel singular values

August 2009

**User's Guide for LMaFit: Low-rank Matrix Fitting**

Yin Zhang

This User's Guide describes the functionality and basic usage of the Matlab package LMaFit for low-rank matrix optimization. It also briefly explains the formulations and algorithms used.

August 2009

**Adaptive Finite Element Methods for Linear-Quadratic Convection Dominated Elliptic Optimal Control Problems**

Eelco Nederkoorn

The numerical solution of linear-quadratic elliptic optimal control problems requires the solution of a coupled system of elliptic partial differential equations (PDEs), consisting of the so-called state PDE, the adjoint PDE and an algebraic equation. Adaptive finite element methods (AFEMs) attempt to locally refine a base mesh in such a way that the solution error is minimized for a given discretization size. This is particularly important for the solution of convection dominated problems where inner and boundary layers in the solutions to the PDEs need to be sufficiently resolved to ensure that the solution of the discretized optimal control problem is a good approximation of the true solution.

This thesis reviews several AFEMs based on energy norm based error estimates for single convection dominated PDEs and extends them to the solution of the coupled system of convection dominated PDEs arising from the optimality conditions for optimal control problems.

__Key Words:__ Adaptive finite element methods, optimal control problems, convection-diffusion equations, local refinement, error estimation

August 2009

**Accelerating the LSTRS Algorithm**

J. Lampe, M. Rojas, D.C. Sorensen, & H. Voss

In a recent paper [Rojas, Santos, Sorensen: ACM ToMS 34 (2008), Article 11] an efficient method for solvingthe Large-Scale Trust-Region Subproblem was suggested which is based on recasting it in terms of a parameter dependent eigenvalue problem and adjusting the parameter iteratively. The essential work at each iteration is the solution of an eigenvalue problem for the smallest eigenvalue of the Hessian matrix (or two smallest eigenvalues in the potential hard case) and associated eigenvector(s). Replacing the implicitly restarted Lanczos method in the original paper with the Nonlinear Arnoldi method makes it possible to recycle most of the work from previous iterations which can substantially accelerate LSTRS.

__Key Words:__ Constrained quadratic optimization, regularization, trust-region, ARPACK, Non-linear Arnoldi method

July 2009

**Application of POD and DEIM to Dimension Reduction of Nonlinear Miscible Viscous Fingering in Porous Media**

Saifon Chaturantabut & Danny C. Sorensen

A Discrete Empirical Interpolation Method (DEIM) is applied in conjunction with Proper Orthogonal Decomposition (POD) to construct a nonlinear reduced-order model of finite difference discretized system used in the simulation of nonlinear miscible viscous fingering in a 2-D porous medium. POD is first applied to extract a low-dimensional basis that optimally captures the dominant characteristics of the system trajectory. This basis is then used in a Galerkin projection scheme to construct a reduced-order system. DEIM is then applied to greatly improve the efficiency in computing the projected nonlinear terms in the POD reduced system. DEIM achieves a complexity reduction of the nonlinearities which is proportional to the number of reduced variables while POD retains a complexity proportional to the original number of variables. Numerical results demonstrate that the dynamics of the viscous fingering in the full-order system of dimension 15000 can be captured accurately by the POD-DEIM reduced system of dimension 40 with the computational time reduced by factor of O(1000).

July 2009

**Domain Decomposition and Balanced Truncation Model Reduction for Shape Optimization of the Stokes System**

H. Antil, M. Heinkenschloss, & R. H. W. Hoppe

The optimal design of structures and systems described by partial differential equations (PDEs) often gives rise to large-scale optimization problems, in particular if the underlying system of PDEs represents a multi-scale, multi-physics problem. Therefore, reduced order modeling techniques such as balanced truncation model reduction, proper orthogonal decomposition, or reduced basis methods are used to significantly decrease the computational complexity while maintaining the desired accuracy of the approximation. In this paper, we are interested in such shape optimization problems where the design issue is restricted to a relatively small portion of the computational domain. In this case, it appears to be natural to rely on a full order model only in that specific part of the domain and to use a reduced order model elsewhere. A convenient methodology to realize this idea consists in a suitable combination of domain decomposition techniques and balanced truncation model reduction. We will consider such an approach for shape optimization problems associated with the time-dependent Stokes system and derive explicit error bounds for the modeling error. As an application in life sciences, we will be concerned with the optimal design of capillary barriers as part of a network of microchannels and reservoirs on microfluidic biochips that are used in clinical diagnostics, pharmacology, and forensics for high-throughput screening and hybridization in genomics and protein profiling in proteomics.

__Key Words:__ Domain decomposition, balanced truncation model reduction, microfluidic biochips.

July 2009

**Domain Decomposition and Model Reduction for the Numerical Solution of PDE Constrained Optimization Problems with Localized Optimization Variables **

H. Antil, M. Heinkenschloss , R. H. W. Hoppe , & D. C. Sorensen

We introduce a technique for the dimension reduction of a class of PDE constrained optimization problems governed by linear time dependent advection diffusion equations for which the optimization variables are related to spatially localized quantities. Our approach uses domain decomposition applied to the optimality system to isolate the subsystem that explicitly depends on the optimization variables from the remaining linear optimality subsystem. We apply balanced truncation model reduction to the linear optimality subsystem. The resulting coupled reduced optimality system can be interpreted as the optimality system of a reduced optimization problem. We derive estimates for the error between the solution of the original optimization problem and the solution of the reduced problem. The approach is demonstrated numerically on an optimal control problem and on a shape optimization problem.

__Key Words:__ Optimal control, shape optimization, domain decomposition, model reduction.

July 2009

**A Random Force is a Force, of Course, of Coarse: Decomposing Complex Enzyme Kinetics with Surrogate Models**

Christopher P. Calderon

The temporal autocorrelation (AC) function associated with monitoring order parameters characterizing conformational fluctuations of an enzyme is analyzed using a collection of surrogate models. The surrogates considered are phenomenological stochastic differential equation (SDE) models. It is demonstrated how an ensemble of such surrogate models, each surrogate being calibrated from a single trajectory, indirectly contains information about unresolved conformational degrees of freedom. This ensemble can be used to construct complex temporal ACs associated with a ``non-Markovian" process. The ensemble of surrogates approach allows researchers to consider models more flexible than a mixture of exponentials to describe relaxation times and at the same time gain physical information about the system. The relevance of this type of analysis to matching single-molecule experiments to computer simulations and how more complex stochastic processes can emerge from a mixture of simpler processes is also discussed. The ideas are illustrated on a toy SDE model and on molecular dynamics simulations of the enzyme dihydrofolate reductase.

June 2009

**TR09-21**
PDF File

**Imaging in Random Waveguides**

L. Borcea, L. Issa (CAAM, Rice University), & C. Tsogka, University of Crete, Greece.

We present a new method for imaging sources in random waveguides, from the time traces of the signals recorded at remote arrays of passive sensors. By random waveguides we mean that the wave speed has rapid fluctuations, that cause wave scattering and a significant loss of coherence of the signals at the array. Our imaging method is based on a special form of transport equations in random waveguides, and it can determine in a statistically stable manner the location of sources in the waveguide and statistical information about the wave speed fluctuations.

June 2009

**Filtering Random Layering Effects For Imaging and Velocity Estimation**

Fernando Gonzalez del Cueto

Imaging compactly supported reflectors in highly heterogeneous media is a challenging problem due to the significant interaction of waves with the medium which causes considerable delay spread and loss of coherence. The imaging problem consists in finding the support of small reflectors using recorded echoes at an array of sensors. The thesis considers the case of randomly layered media, in which significant multiple scattering by the layered structures and quick loss of coherence is observed. These strong, backscattered echoes from the layers can overwhelm the weaker coherent signals due to the compactly supported reflectors. This signal-to-noise problem must be addressed to image effectively. Using techniques routinely used in exploration seismology, filters (layer annihilators) are designed to remove the primary reflections of the stronger layered features in the medium. However, it observed that these filters also remove the incoherent signal that is due to the fine, random layers. The main achievement of this thesis is the theoretical and numerical analysis of this phenomenon. Additionally, the applicability of the layer annihilators for velocity estimation is presented.

May 2009

**A Simple Algorithm For the Inverse Field of Values Problem**

Russell Carden

The field of values of a matrix is the closed convex subset of the complex plane comprising all Rayleigh quotients, a set of interest in the stability analysis of dynamical systems and convergence theory of matrix iterations, among other applications. Recently, Uhlig proposed the inverse field of values problem: given a point in the field of values, determine a vector for which this point is the corresponding Rayleigh quotient. Uhlig also devised a sophisticated algorithm involving random vectors and the boundaries of ellipses for solving the inverse field of values problem. We propose a simpler deterministic algorithm that must converge (in exact arithmetic), and for most points yields an exact result in only a few iterations. The algorithm builds upon the fact that the inverse field of values problem can be solved exactly in the two dimensional case. We also resolve a conjecture posed by Uhlig concerning the number of linearly independent vectors that generate a point in the field of values, and propose a more challenging inverse field of values problem that is of interest in eigenvalue computations.

May 2009

**Convergence of a Discontinuous Galerkin Method For the Miscible Displacement Under Minimal Regularity**

Béatrice M. Rivière & Noel Walkington

Discontinuous Galerkin time discretizations are combined with the mixed finite element and continuous finite element methods to solve the miscible displacement problem. Stable schemes of arbitrary order in space and time are obtained. Under minimal regularity assumptions on the data, convergence of the scheme is proved by using compactness results for functions that may be discontinuous in time.

May 2009

**User's Guide For YALL1: Your Algorithms for L1 Optimization**

Yin Zhang

This User's Guide describes the functionality and basic usage of the Matlab package YALL1 for L1 minimization. The one-for-six algorithm used in the YALL1 solver is briefly introduced in the appendix.

May 2009

**A Parameter Free ADI-like Method for the Numerical Solution of Large Scale Lyapunov Equations**

Ryan Nong & Danny C. Sorensen

An algorithm is presented for constructing an approximate numerical solution to a large scale Lyapunov equation in low rank factored form. The algorithm is based upon a synthesis of an approximate power method and an alternating direction implicit (ADI) method. The former is parameter free and tends to be efficient in practice, but there is little theoretical understanding of its convergence properties. The ADI method has a well understood convergence theory, but the method relies upon selection of shift parameters and a poor shift selection can lead to very slow convergence in practice. The algorithm presented here uses an approximate power method iteration to obtain a basis update. It then constructs a re-weighting of this basis update to provide a factorization update that satisfies ADI-like convergence properties.

May 2009

**PSQR: A Stable and Efficient Penalized Spline Algorithm**

Christopher P. Calderon, Josue G. Martinez, Raymond J. Carroll, & Danny C. Sorensen

We introduce an algorithm for reliably computing quantities associated with several types of semiparametric mixed models in situations where the condition number on the random effects matrix is large. The algorithm is numerically stable and efficient. It was designed to process penalized spline (P-spline) models without making unnecessary numerical approximations. The algorithm, PSQR (P-splines via QR), is formulated in terms of QR decompositions. PSQR can treat both exactly rank deficient and ill-conditioned matrices. The latter situation often arises in large scale mixed models and/or when a P-spline is estimated using a basis with poor numerical properties, e.g. a truncated power function (TPF) basis. We provide concrete examples where unnecessary numerical approximations introduce both subtle and dramatic errors that would likely go undetected, thus demonstrating the importance of using this reliable numerical algorithm. Simulation results studying a univariate function and a longitudinal data set are used to demonstrate the algorithm. Extensions and the utility of the method in more general semiparametric regression applications are briefly discussed. MATLAB scripts demonstrating implementation are provided in the Supplemental Materials.

May 2009

**Analysis of Two Mathematical Models for the Coupled Navier-Stokes/Darcy Problem**

Prince Chidyagwai & Béatrice M. Rivière

This paper introduces and analyzes two models coupling of incompressible Navier-Stokes equations with the porous media flow equations. A numerical method that uses continuous finite elements in the incompressible flow region and discontinuous finite elements in the porous medium, is proposed. Existence and uniqueness results under small data condition of the numerical solution are proved. Optimal a priori error estimates are derived. Numerical examples comparing the two models are provided.

April 2009

**P-Splines Using Derivative Information**

Christopher P. Calderon, Josue G. Martinez, Raymond J. Carroll, & Danny C. Sorensen

Time series associated with single-molecule experiments and/or simulations contain a wealth of multiscale information about complex biochemical systems. However efficiently extracting and representing useful physical information from these time series measurements can be challenging. We demonstrate how Penalized splines (P-Splines) can be useful in summarizing complex single-molecule time series data using quantities estimated from the observed data. A design matrix that simultaneously uses noisy function and derivative scatterplot information to refine function estimates using P-spline techniques is introduced. The approach is called the PuDI (P-Splines using Derivative Information) method. We show how Generalized Least Squares fits seamlessly into the PuDI method; several applications demonstrating how inclusion of uncertainty information improves the PuDI function estimates are presented. The PuDI design matrix can be used to assist scatterplot smoothing applications where both unbiased function and derivative estimates are available.

April 2009

**Morphologically Accurate Reduced Order Modeling of Spiking Neurons**

Anthony R. Kellems, Saifon Chaturantabut, Danny C. Sorensen, & Steven J. Cox

Accurately simulating neurons with realistic morphological structure and synaptic inputs requires the solution of large systems of nonlinear ordinary differential equations. We apply model reduction techniques to recover the complete nonlinear voltage dynamics of a neuron using a system of much lower dimension. Using a proper orthogonal decomposition, we build a reduced-order system from salient snapshots of the full system output, thus reducing the number of state variables. A discrete empirical interpolation method is then used to reduce the complexity of the nonlinear term to be proportional to the number of reduced variables. Together these two techniques allow for up to two orders of magnitude dimension reduction without sacrificing the spatially-distributed input structure, with an associated order of magnitude speed-up in simulation time. We demonstrate that both nonlinear spiking behavior and subthreshold response of realistic cells are accurately captured by these low-dimensional models.

April 2009

**Getting It Right Without Knowing the Answer: Quality Control in a Large Seismic Modeling Project**

William W. Symes, Igor S. Terentyev, & Tetyana Vdovina

Phase I of the SEAM Project will produce a variable-density acoustic synthetic survey over a 3D geological model simulating a deepwater subsalt exploration target. Due to the intended use of the data, the project places a premium on accuracy. Commercially produced Phase I synthetics will be spot-checked against benchmark simulations to assure quality. Thus the accuracy of the benchmark simulator required careful assessment. The authors designed and implemented the benchmark simulator used in this program, subjected it to verification tests, and assisted in the qualification phase of the Phase I project. The key lessons that we have learned so far from this assessment are that (1) the few verification tools available to us - a few analytic solutions and Richardson extrapolation - seem to be adequate, at least in a rough way, and (2) the standard approach to this type of simulation - finite difference methods on regular grids - requires surprisingly fine grid steps to produce small relative RMS errors for models of the type defined by this project.

April 2009

**Mass Lumping For Constant Density Acoustics**

William W. Symes

Mass lumping provides an avenue for efficient time-stepping of time-dependent problems with conforming finite element discretization. Typical justifications for mass lumping use quadrature error estimates which do not hold for nonsmooth coefficients. In this paper, I show that the mass-lumped semidiscrete system for the constant-density acoustic wave equation with Q1 elements exhibits optimal order convergence eeven when the coefficient (bulk modulus) is merely bounded and measurable, provide that the right-hand side possesses some smoothness in time.

April 2009

**Approximate Inverse Scattering Using Pseudodifferential Scaling**

Rami Nammour

This thesis proposes a computationally efficient method for approximating the inverse of the normal operator arising in the linearized inverse problem for reflection seismology.

The inversion of the normal operator using direct matrix methods is computationally infeasible. Approximate inverses estimate the solution of the inverse problem or precondition iterative methods. Application of the normal operator requires an expensive solution of large scale PDE problems. However, the normal operator approximately commutes with pseudodifferential operators, hence shares their near diagonality in a frame of localized monochromatic pulses. Estimation of a diagonal representation in this frame encoded in the symbol of the normal operator:

- follows from its application to a single input vector;
- suffices to approximate its inverse.

I use an efficient algorithm to apply pseudodifferential operators, given their symbol, to construct a rapidly converging optimization algorithm that estimates the symbol of an inverse for the normal operator, thereby approximately solving the inverse problem.

April 2009

**Approximate Constant Density Acoustic Inverse Scattering Using Dip-Dependent Scaling**

Rami Nammour & William W. Symes

This abstract presents a computationally efficient method to approximate the inverse of the Hessian or normal operator arising in a linearized inverse problem for constant density acoustics model of reflection seismology. Solution of the linearized inverse problem problem involves construction of an image via prestack depth migration, then correction of the image amplitudes via application of the inverse of the normal operator. The normal operator acts by dip-dependent scaling of the amplitudes of its input vector. This property permits us to efficiently approximate the normal operator, and its inverse, from the result of its application to a single input vector, for example the image, and thereby approximately solve the linearized inverse scattering problem.

We validate the method on a 2D section of the Marmousi model to correct the amplitudes of the migrated image.

April 2009

**A Software Framework for Finite Difference Simulation**

Igor Terentyev

This paper describes a software framework for solving time dependent PDEs in simple domains using finite difference (FD) methods. The framework is designed for parallel computations on distributed and shared memory computers, thus allowing for efficient solution of large-scale problems. The framework provides automated data exchange between processors based on stencil information. This automated data exchange allows a user to add FD schemes without knowledge about underlying parallel infrastructure. The framework includes acoustic solver based on staggered second-order in time and various orders in space FD schemes with perfectly matched layer and/or free surface boundary conditions.

April 2009

**Mass Lumping for Constant-Density Acoustics**

William W. Syme & Igor Terentyev

Conforming Galerkin discretization of the constant-density acoustic wave equation provides optimal order convergence even in the presence of very rough coefficients, provided that the time dependence of the data (right-hand side) is minimally smooth. Such discretizations avoid the well-known first-order error ("stairstep diffraction") phenomenon produced by standard finite difference methods. On the other hand, Galerkin methods in themselves are inefficient for high frequency wave simulation, due to the implicit nature of the time step system. Mass lumping renders the time step explicit, and provides an avenue for efficient time-stepping of time-dependent problems with conforming finite element spatial discretization. Typical justifications for mass lumping use quadrature error estimates which do not hold for nonsmooth coefficients. In this paper, we show that the mass-lumped semidiscrete system for the constant-density acoustic wave equation with rectangular multilinear elements exhibits optimal order convergence even when the coefficient (bulk modulus) is merely bounded and measurable, provided that the right-hand side possesses some smoothness in time. We illustrate the theory with numerical examples involving discontinuous, non-grid-aligned bulk moduli, in which the coefficient averaging implicit in mass lumping eliminates the stairstep diffraction effect.

April 2009

**Discrete Empirical Interpolation for Nonlinear Model Reduction**

Saifon Chaturantabut & Danny C. Sorensen

A dimension reduction method called Discrete Empirical Interpolation (DEIM) is proposed and shown to dramatically reduce the computational complexity of the popular Proper Orthogonal Decomposition (POD) method for constructing reduced-order models for unsteady and/or parametrized nonlinear partial differential equations (PDEs). In the presence of a general nonlinearity, the standard POD-Galerkin technique reduces dimension in the sense that far fewer variables are present, but the complexity of evaluating the nonlinear term remains that of the original problem. Here we describe DEIM as a modification of POD that reduces the complexity as well as the dimension of general nonlinear systems of ordinary differential equations (ODEs). It is, in particular, applicable to ODEs arising from finite difference discretization of unsteady time dependent PDE and/or parametrically dependent steady state problems. Our contribution is a greatly simplified description of Empirical Interpolation in a finite dimensional setting. The method possesses an error bound on the quality of approximation. An application of DEIM to a finite difference discretization of the 1-D FitzHugh-Nagumo equations is shown to reduce the dimension from 1024 to order 5 variables with negligible error over a long-time integration that fully captures non-linear limit cycle behavior. We also demonstrate applicability in higher spatial dimensions with similar state space dimension reduction and accuracy results.

March 2009

**The Nonlinear Differential Semblance Algorithm for Plane Waves in Layered Media**

Dong Sun

This paper proposes an alternative approach to the output least-squares (OLS) seismic inversion for layered-media. The latter cannot guarantee a reliable solution for either synthetic or field data, because of the existence of many spurious local minima of the objective function for typical data, which lack low-frequency energy. To recover the low-frequency lacuna of typical data, I formulate waveform inversion into a differential semblance optimization (DSO) problem with artificial low-frequency data as control variables. To my knowledge, this approach is the first version of differential semblance with nonlinear modeling that may properly accounts for nonlinear effects of wave propagation, such as multiple reflections. Numerical experiments with synthetic data indicate the smoothness and convexity of the proposed objective function. These results suggest that gradient-related algorithms may successfully approximate a global minimizer from a crude initial guess for typical band-limited data.

April 2009

**A Nonlinear Differential Semblance Strategy for Waveform Inversion: Experiments in Layered Media**

Dong Sun & William W. Symes

This paper proposes an alternative approach to the output least-squares (OLS) seismic inversion for layered-media. The latter cannot guarantee a reliable solution for either synthetic or field data, because of the existence of many spurious local minima of the objective function for typical data, which lack low-frequency energy. To recover the low-frequency lacuna of typical data, we formulate waveform inversion as a differential semblance optimization (DSO) problem with artificial low-frequency data as control variables. This version of differential semblance with nonlinear modeling properly accounts for nonlinear effects of wave propagation, such as multiple reflections. Numerical experiments with synthetic data indicate the smoothness and convexity of the proposed objective function. These results suggest that gradient-related algorithms may successfully approximate a global minimizer from a crude initial guess for typical band-limited data.

April 2009

**Analysis and Generalizations of the Linearized Bregman Method**

Wotao Yin

This paper reviews the Bregman methods, analyzes the linearized Bregman method, and proposes fast generalization of the latter for solving the basis pursuit and related problems. The analysis shows that the linearized Bregman method has the exact penalty property, namely, it converges to an exact solution of the basis pursuit problem if and only if its regularization parameter a is greater than a certain value. The analysis is based on showing that the linearized Bregman algorithm is equivalent to gradient descent applied to a certain dual formulation. This result motivates generalizations of the algorithm enabling the use of gradient-based optimization techniques such as line search, Barzilai-Borwein steps, L-BFGS, and nonlinear conjugate gradient steps. In addition, the paper discusses the selection and update of alpha.

The analysis and discussions are limited to the L1-norm but can be extended to other L-like functions.

May 2009

**A Fast Algorithm For Sparse Reconstruction Based On Shrinkage Subspace Optimization And Continuation**

Zaiwen Wen, Wotao Yin, Donald Goldfarb, & Yin Zhang

We propose a fast algorithm for solving the l1-regularized minimization problem $\min \mu \|x\|_1 + \|Ax -b\|^2_2$ for recovering sparse solutions to an undetermined system of linear equations $Ax=b$. The algorithm is divided into two stages that are performed repeatedly. In the first stage a first-order iterative method called ``shrinkage'' yields an estimate of the subset of components of $x$ likely to be nonzero in an optimal solution. Restricting the decision variables $x$ to this subset and fixing their signs at their current values reduces the l1-norm $\|x\|_1$ to a linear function of $x$. The resulting subspace problem, which involves the minimization of a smaller and smooth quadratic function, is solved in the second phase. Our code FPC-AS embeds this basic two-stage algorithm in a continuation (homotopy) approach by assigning a decreasing sequence of values to $\mu$. This code exhibits state-of-the-art performance both in terms of its speed and its ability to recover sparse signals. It can even recover signals that are not as sparse as required by current compressive sensing theory.

February 2009