Technical Reports (2009)
Analysis of Weak Solutions For the Fully Coupled Stokes-Darcy-Transport Problem
Aycil Cesmelioglu & Beatrice Riviere
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
TR09-33 PDF
An Overview of Time-Stepping Classes for Optimization (TSOpt)
Marco Enriquez & William W. Symes
Abstract TBA
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
Compressed Sensing via Iterative Support Detection
Yilun Wang & Wotao Yin
We present a new compressive sensing 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 cases of L1-based construction due to insufficient measurements, in which the returned signals are not equal or even close to the true signals. One would normally give up on such returned signals, but ISD will learn from them and come up with new minimization problems that return signals closer to the true ones. Specifically, from an incorrect signal ISD detects an index set I that includes components most likely to be true nonzeros, obtains a new signal x by solving min{sum_{i not in I} |x_i| : Ax = b}, and repeats such support detection and minimization using latest x and I from one another until convergence. ISD differs from the orthogonal matching pursuit (OMP) method and its variants in two aspects.
First, although all these methods generate index sets iteratively, I in ISD is not necessarily nested or increasing over iterations. Second, the OMP-family methods iteratively project certain quantities (typically residuals) over the complements of the index sets but the minimization problem above in ISD is no projection. To analyze the performance of the proposed method, we generalize the Null Space Property to Truncated Null Space Property and base our analysis on the latter.
We introduce an efficient implementation of ISD, called threshold-ISD, for recovering signals with fast decaying distributions of nonzeros from compressive measurements. Numerical experiments show that threshold{ISD has significant overall advantages over the classical L1 minimization approach, as well as two other state-of-the-art algorithms such as the iterative reweighted L1 minimization algorithm (IRL1) and the iterative reweighted least-squares algorithm (IRLS).
September 2009
On Balanced Truncation for Inhomogeneously Initialized Systems
M. Heinkenschloss, T. Reis, & A. C. Antoulas
We introduce a balancing-related method for the reduction of linear time-invariant systems whose initial value is not necessarily vanishing. An error bound for the presented reduction method is provided.
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
TR09-25 PDF File
Application of POD and DEIM to Dimension Reduction of Nonlinear Miscible Viscous Fingering in Porous Media
Saifon Chaturantabut & Danny C. Sorensen
Proper Orthogonal Decomposition (POD) in conjunction with a Discrete Empirical Interpolation Method (DEIM) is applied to construct a reduced-order model of a finite difference discretized system used to simulate nonlinear miscible viscous fingering in a 2-D porous medium. POD is first used to extract a low-dimensional basis that optimally captures the dominant characteristics of the sampled trajectory of the system. This POD basis is truncated according to decay of the singular values and the resulting low dimensional reduced basis is then used in a Galerkin projection to construct a reduced-order system. However, this POD based reduced system still has complexity of the full system present in evaluation of the projected nonlinear term and hence provides no computational savings. DEIM is applied to reduce the complexity of the projected nonlinear term to be proportional to the number of reduced 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
TR09-09 PDF File
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
Subgrid Modeling Via Mass Lumping in Constant Density Acoustics
William W. Syme & Igor Terentyev
Finite difference simulations in discontinuous media are only first-order accurate regardless of the formal order of the method. Stairstep diffractions are a widely known manifestation of this first-order error. To recover the optimal convergence rate, we apply a finite element approach with mass lumping to a constant density acoustic wave equation. This approach results in an explicit second-order accurate difference scheme with specifically averaged sound velocity. Two numerical experiments confirm theoretical convergence rate. In particular, application of finite element discretizations with mass lumping leads to elimination of stairstep diffraction observed in simulations based on standard finite difference discretizations.
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
TR09-04 PDF File
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