Titles By Year
Abstracts By Year
Technical Reports (2010)
TR10-35 PDF File
Mass Lumping for Constant Density Acoustics
William W. Symes & Igor S. 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 Q1 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, in which the stairstep diffraction effect is eliminated by the coefficient averaging implicit in mass lumping.
December 2010
The Effects of Coupling Adaptive Time-Stepping and Adjoint-State Methods for Optimal Control Problems
Marco U. Enriquez
This thesis presents the implications of using adaptive time-stepping schemes with the adjoint-state method, a widely used algorithm for computing derivatives in optimal-control problems. Though we gain control over the accuracy of the time- stepping scheme, the forward and adjoint time grids become mismatched. Despite this fact, I claim using adaptive time-stepping for optimal control problems is advantageous for two reasons. First, taking variable time-steps potentially reduces the computational cost and improves accuracy of the forward and adjoint equations’ numerical solution. Second, by appropriately adjusting the tolerances of the time-stepping scheme, convergence of the optimal control problem can be theoretically guaranteed via inexact Newton theory. I present proofs and computational results to support this claim.
The computational results include an extension of prior work on adaptive checkpointing schemes, enabling checkpointing when solving the reference and adjoint equations adaptively. The numerical results in this thesis feature an optimal control problem with a reservoir simulation constraint.
December 2010
Coherent Interferometric Imaging, Time Gating and Beamforming
Liliana Borcea, Josselin Garnier, George Papanicolaou, & Chrysoula Tsogka
Coherent interferometric imaging is based on the backpropagation of local space-time cross correlations of array data and was introduced in order to improve images when the medium between the array and the object to be imaged is inhomogeneous and unknown [Borcea et al., Inverse Problems, 21 (2005), 1419]. Although this method has been shown to be effective and is well founded theoretically, the coherent interferometric imaging function is computationally expensive and therefore difficult to use. In this paper we show that this function is equivalent to a windowed beamformer energy function, that is, a quadratic function that involves only time gating and time delaying signals in emission and in reception. In this form the coherent interferometric imaging can be implemented efficiently both in hardware and software, that is, at a computational cost that is comparable to the usual beamforming and migration imaging methods. We also revisit the trade-off between enhanced image stability and loss of resolution in coherent interferometry from the point of view of its equivalence to a windowed beamformer energy imaging function.
December 2010
A State Space Error Estimate for POD-DEIM Nonlinear Model Reduction
Saifon Chaturantabut & Danny C. Sorensen
This paper derives state space error bounds for the solutions of reduced systems constructed using Proper Orthogonal Decomposition (POD) together with the Discrete Empirical Interpolation Method (DEIM) recently developed in [S. Chaturantabut and D.C. Sorensen, SISC 2010, pp. 2737-2764 ] for nonlinear dynamical systems. The resulting error bounds are shown to be proportional to the sums of the singular values corresponding to neglected POD basis vectors both in Galerkin projection of the reduced system and in the DEIM approximation of the nonlinear term. The analysis is particularly relevant to ODE systems arising from spatial discretizations of parabolic PDEs. The derivation clearly identifies where the parabolicity is crucial. It also shows exactly how the DEIM approximation error involving the nonlinear term comes into play.
S. Chaturantabut and D.C. Sorensen, Nonlinear model reduction via discrete empirical interpolation, SIAM Journal on Scientific Computing, 32 (2010), pp. 2737-2764.December 2010
A Rapid and Robust Numerical Algorithm for Sensitivity Encoding with Sparsity Constraints: Self-Feeding Sparse SENSE
Feng Huang, Yunmei Chen, Wotao Yin, Wei Lin, Xiaojing Ye, Weihong Guo, & Arne Reykowski
The method of enforcing sparsity during magnetic resonance imaging reconstruction has been successfully applied to partially parallel imaging (PPI) techniques to reduce noise and artifact levels and hence to achieve even higher acceleration factors. However, there are two major problems in the existing sparsity-constrained PPI techniques: speed and robustness. By introducing an auxiliary variable and decomposing the original minimization problem into two subproblems that are much easier to solve, a fast and robust numerical algorithm for sparsity-constrained PPI technique is developed in this work. The specific implementation for a conventional Cartesian trajectory data set is named self-feeding Sparse Sensitivity Encoding (SENSE). The computational cost for the proposed method is two conventional SENSE reconstructions plus one spatially adaptive image denoising procedure.
With reconstruction time approximately doubled, images with a much lower root mean square error (RMSE) can be achieved at high acceleration factors. Using a standard eight-channel head coil, a net acceleration factor of 5 along one dimension can be achieved with low RMSE. Furthermore, the algorithm is insensitive to the choice of parameters. This work improves the clinical applicability of SENSE at high acceleration factors.
November 2010
Trust, But Verify: Fast and Accurate Signal Recovery from 1-bit Compressive Measurements
Jason N. Laska, Zaiwen Wen, Wotao Yin, & Richard G. Baraniuk
The recently emerged compressive sensing (CS) framework aims to acquire signals at reduced sample rates compared to the classical Shannon-Nyquist rate. To date, the CS theory has assumed primarily real-valued measurements; it has recently been demonstrated that accurate and stable signal acquisition is still possible even when each measurement is quantized to just a single bit. This property enables the design of simplified CS acquisition hardware based around a simple sign comparator rather than a more complex analog-to-digital converter; moreover, it ensures robustness to gross non-linearities applied to the measurements. In this paper we introduce a new algorithm — restricted-step shrinkage (RSS) — to recover sparse signals from 1-bit CS measurements. In contrast to previous algorithms for 1-bit CS, RSS has provable convergence guarantees, is about an order of magnitude faster, and achieves higher average recovery signal-to-noise ratio. RSS is similar in spirit to trust-region methods for non-convex optimization on the unit sphere, which are relatively unexplored in signal processing and hence of independent interest.
November 2010
Hyperspectral Data Reconstruction Combining Spatial and Spectral Sparsity
SeyyedMajid Valiollahzadeh & Wotao Yin
This report introduces a novel sparse decomposition model for hyperspectral image reconstruction. The model integrates two well-known sparse structures of hyperspectral images: a small set of signature spectral vectors span all spectral vectors (one at each pixel), and like a standard image, a hyperspectral image is spatially redundant. In our model, a threedimensional hyperspectral cube X is first decomposed into a small number of endmembers by X = H beta + n, where H is the endmember dictionary, beta contains the coefficients, and n are errors and noise. Then beta, which is a 3D cube with the same spatial dimensions as X, is further decomposed into overlapping cubelets {beta_i}, which are sparsely represented by a common dictionary D, i.e., beta_i = D alpha_i where alpha_i is a set of sparse coefficients. This model not only exploits spectral sparsity to the original hyperspectral cube X to the smaller cube beta but also applies latest image sparse representation techniques to beta.
Given a corrupted hyperspectral cube with noise and missing voxels, our method reconstructs the cube by learning H, D, and alpha_i’s from the data itself. These parameters are statistically modeled such that alpha_i’s are sparse and Bayesian inferences have closed-form formulas and are thus easy to compute. Numerical simulations were performed on AVIRIS images. We show that merely 5% randomly selected voxels are enough for the proposed method to returned state-of-theart reconstruction results.
November 2010
Oil Spill Sensor using Multispectral Infrared Imaging via L1 Minimization
Yingying Li, Wei-Chuan Shih, Zhu Han, and Wotao Yin
Early detection of oil spill events is the key to environmental protection and disaster management. Current technology lacks the sensitivity and specificity in detecting the early onset of a small-scale oil spill event. Based on an infrared oil-water contrast model recently developed, we propose a novel nonscanning computational infrared sensor that has the potential to achieve unprecedented detection sensitivity. Such a system can be very low-cost and robust for automated outdoor operations, leading to massive offshore deployment. Taking advantage of the characteristic oil thickness multispectral signatures, we have streamlined an algorithm that incorporates 3D image reconstruction and classification in a single inversion step capitalizing on the benefits of L1 minimization.
November 2010
High Resolution OFDM Channel Estimation with Low Speed ADC using Compressive Sensing
Jia (Jasmine) Meng, Yingying Li, Nam Nguyen, Wotao Yin, & Zhu Han
Orthogonal frequency division multiplexing (OFDM) is a technique that will prevail in the next generation wireless communication. Channel estimation is one of the key challenges in an OFDM system. In this paper, we formulate OFDM channel estimation as a compressive sensing problem, which takes advantage of the sparsity of the channel impulse response and reduces the number of probing measurements, which in turn reduces the ADC speed needed for channel estimation. Specifically, we propose sending out pilots with random phases in order to "spread out" the sparse taps in the impulse response over the uniformly downsampled measurements at the low speed receiver ADC, so that the impulse response can still be recovered by sparse optimization. This contribution leads to high resolution channel estimation with low speed ADCs, distinguishing this paper from the existing attempts of OFDM channel estimation. We also propose a novel estimator that performs better than the commonly used L1 minimization. Specifically, it significantly reduces estimation error by combing L1 minimization with iterative support detection and limited-support least-squares. While letting the receiver ADC run at a speed as low as 1/16 of the speed of the transmitter DAC, we simulated various numbers of multipaths and different measurement SNRs. The proposed system has channel estimation resolution as high as the system equipped with the high speed ADCs, and the proposed algorithm provides additional 6 dB gain for signal to noise ratio.
November 2010
A Feasible Method for Optimization with Orthogonality Constraints
Zaiwen Wen & Wotao Yin
Minimization with orthogonality constraints (e.g., X'X = I) and/or spherical constraints (e.g., ||x||_2 = 1) has wide applications in polynomial optimization, combinatorial optimization, eigenvalue problems, sparse PCA, p-harmonic flows, 1-bit compressive sensing, matrix rank minimization, etc. These problems are difficult because the constraints are not only non-convex but numerically expensive to preserve during iterations. To deal with these difficulties, we propose to use a Crank-Nicholson-like update scheme to preserve the constraints and based on it, develop curvilinear search algorithms with lower per-iteration cost compared to those based on projections and geodesics. The efficiency of the proposed algorithms is demonstrated on a variety of test problems. In particular, for the maxcut problem, it exactly solves a decomposition formulation for the SDP relaxation. For polynomial optimization, nearest correlation matrix estimation and extreme eigenvalue problems, the proposed algorithms run very fast and return solutions no worse than those from their state-of-the-art algorithms. For the quadratic assignment problem, a gap 0.842% to the best known solution on the largest problem "256c" in QAPLIB can be reached in 5 minutes on a typical laptop.
November 2010
Reconstructing an Even Damping from a Single Spectrum
Steven J. Cox & Mark Embree
We consider the wave equation on a finite interval with fixed ends and nonuniform viscous damping. We prove that the spectrum of the associated damped wave operator uniquely determines an even damping. We then develop a refined asymptotic formula for the high frequencies. When the damping is even about the domain's midpoint, terms in this expansion are Fourier coefficients for functions of the damping. Furthermore, the expansion is often quite accurate even at low frequencies, thus suggesting a simple numerical procedure for reconstructing even damping coefficients from measured eigenvalues. Computational examples illustrate the efficacy of this procedure, even in the presence of noise.
August 2010
Convergence of a Class of Stationary Iterative Methods for Saddle Point Problems
Yin Zhang
A unified convergence result is derived for an entire class of stationary iterative methods for solving equality constrained quadratic programs or saddle point problems. This class is constructed from essentially all possible splittings of the n × n submatrix residing in the (1,1)-block of the (n+m)×(n+m) augmented matrix that would generate non-expansive iterations in R^n. The classic multiplier method and augmented Lagrangian alternating direction method are two special members of this class.
August 2010
TR10-23 PDF File
Mathematical Analysis of a Multiphysics Problem
Aycil Cesmelioglu & Beatrice Riviere
This paper analyzes the surface/subsurface flow coupled with transport. The flow is modeled by the coupling of Navier-Stokes and Darcy equations. The transport of a species 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.
July 2010
On the Convergence of an Active Set Method for L1 Minimization
Zaiwen We, Wotao Yin, Hongchao Zhang, & Donald Goldfarb
We analyze an abridged version of the active-set algorithm FPC_AS for solving the L1-regularized least squares problem. The active set algorithm alternatively iterates between two stages. In the first "nonmonotone line search (NMLS)" stage, an iterative first-order method based on "shrinkage" is used to estimate the support at the solution. In the second "subspace optimization"stage, a smaller smooth problem is solved to recover the magnitudes of the nonzero components of the solution x. We show that NMLS itself is globally convergent and the convergence rate is at least R-linearly. In particular, NMLS is able to identify of the zero components of a stationary point after a finite number of steps under some mild conditions. The global convergence of FPC_AS is established based on the properties of NMLS.
July 2010
A Two-grid Method for Coupled Free Flow with Porous Media Flow
Prince Chidyagwai & Beatrice Riviere
This paper presents a two-grid method for solving systems of partial differential equations modelling free flow coupled with porous media flow. This work considers both the coupled Stokes and Darcy as well as the coupled Navier-Stokes and Darcy problems. The numerical schemes proposed are based on combinations of the continuous finite element method and the discontinuous Galerkin method. Numerical errors and convergence rates for solutions obtained from the two-grid method are presented. CPU times for the two-grid algorithm are shown to be significantly less than those obtained by solving the fully coupled problem.
June 2010
Optimization Governed by Stochastic Partial Differential Equations
Drew Kouri
This thesis provides a rigorous framework for the solution of stochastic elliptic partial differential equation (SPDE) constrained optimization problems. In modeling physical processes with differential equations, much of the input data is uncertain (e.g. measurement errors in the diffusivity coefficients). When uncertainty is present, the governing equations become a family of equations indexed by a stochastic variable. Since solutions of these SPDEs enter the objective function, the objective function usually involves statistical moments. These optimization problems governed by SPDEs are posed as a particular class of optimization problems in Banach spaces. This thesis discusses Monte Carlo, stochastic Galerkin, and stochastic collocation methods for the numerical solution of SPDEs and identifies the stochastic collocation method as particularly useful for the optimization of SPDEs. This thesis extends the stochastic collocation method to the optimization context and explores the decoupling nature of this method for gradient and Hessian computations.
June 2010
On the Constants in Inverse Inequalities in L2
Sevtap Ozisik, Beatrice Riviere, & Tim Warburton
In this paper we determine the constants in multivariate Markov inequalities in the L2-norm on an interval, a triangle and a tetrahedron. Using orthonormal polynomials, we derive explicit expression for the constants on that given 1-simplex, 2-simplex and 3-simplex. Accurate values for the constants are crucial for the correct derivation of a priori and a posteriori error estimations in adaptive computation.
May 2010
Resistor Networks and Optimal Grids for the Numerical Solution of Electrical Impedance Tomography with Partial Boundary Measurements
Alexander Vasilyevich Mamonov
The problem of Electrical Impedance Tomography (EIT) with partial boundary measurements is to determine the electric conductivity inside a body from the simultaneous measurements of direct currents and voltages on a subset of its boundary. Even in the case of full boundary measurements the non-linear inverse problem is known to be exponentially ill-conditioned. Thus, any numerical method of solving the EIT problem must employ some form of regularization. We propose to regularize the problem by using sparse representations of the unknown conductivity on adaptive finite volume grids known as the optimal grids, that are computed as part of the problem. Then the discretized partial data EIT problem can be reduced to solving the discrete inverse problems for resistor networks. Two distinct approaches implementing this strategy are presented.
The first approach uses the results for the EIT problem with full boundary measurements, which rely on the use of resistor networks with circular graph topology. The optimal grids for such networks are essentially one dimensional objects, which can be computed explicitly using discrete Fourier transform and rational interpolation with continued fractions. We solve the partial data problem by reducing it to the full data case using the theory of extremal quasiconformal (Teichmuller) mappings.
The second approach is based on resistor networks with the pyramidal graph topology. Such network topology is better suited for the partial data problem, since it allows for explicit treatment of the inaccessible part of the boundary. We present a method of computing the optimal grids for the networks with general topology (including pyramidal), which is based on the sensitivity analysis of both the continuum and the discrete EIT problems. This is the first study of the optimal grids for the case, where reduction to one dimension is not possible.
We present extensive numerical results for the two approaches. We demonstrate both the optimal grids and the reconstructions of smooth and discontinuous conductivities in a variety of domains. The numerical results show several advantages of our approaches compared to the traditional optimization-based methods. First, the inversion based on resistor networks is orders of magnitude faster than any iterative algorithm. Second, our approaches are able to correctly reconstruct both smooth and discontinuous conductivities including those of very high contrast, which usually present a challenge to the iterative or linearization-based inversion methods. Finally, our method does not require any form of artificial regularization via penalty terms. However, our method allows for such regularization to incorporate prior information in the solution.
May 2010
Detection and Imaging in Strongly Backscattering Randomly Layered Media
R. Alonso, L. Borcea, G. Papanicolaou, & C. Tsogka
Echoes from small reflectors buried in heavy clutter are weak and difficult to distinguish from the medium backscatter. Detection and imaging with sensor arrays in such media requires filtering out the unwanted backscatter and emphasizing the echoes from the reflectors that we wish to locate. We consider a detection and filtering approach based on the singular value decomposition of the local cosine transform of the array response matrix. The algorithm applies to general clutter, but its analysis depends on the model of the medium. This paper is concerned with the analysis of the algorithm in finely layered random media. We obtain a detailed characterization of the singular values of the transformed array response matrix and justify the systematic approach of the algorithm for detecting and refining the time windows that contain the echoes that are useful in imaging.
May 2010
Optimization of Shell Structure Acoustics
Sean S. Hardesty
This thesis analyzes a mathematical model for shell structure acoustics, and develops and implements the adjoint equations for this model. The adjoint equations allow the computation of derivatives with respect to large parameter sets in shape optimization problems where the thickness and mid-surface of the shell are computed so as to generate a radiated sound field subject to broad-band design requirements.
The structure and acoustics are modeled, respectively, via the Naghdi shell equations, and thin boundary integral equations, with full coupling at the shell mid-surface. In this way, the three-dimensional structural-acoustic equations can be posed as a problem on the two-dimensional mid-surface of the shell. A wide variety of shapes can thus be explored without re-meshing, and the acoustic field can be computed anywhere in the exterior domain with little additional effort. The problem is discretized using triangular MITC shell elements and piecewise-linear Galerkin boundary elements, coupled with a simple one-to-one scheme.
Prior optimization work on coupled shell-acoustics problems has been focused on applications with design requirements over a small range of frequencies. These problems are amenable to a number of simplifying assumptions. In particular, it is often assumed that the structure is dense enough that the air pressure loading can be neglected, or that the structural motions can be expanded in a basis of low-frequency eigenmodes. Optimization of this kind can be done with reasonable success using a small number of shape parameters because simple modal analysis permits a reasonable knowledge of the parts of the design that will require modification. None of these assumptions are made in this thesis. By using adjoint equations, derivatives of the radiated field can be efficiently computed with respect to large numbers of shape parameters, allowing a much richer space of shapes, and thus, a broader range of design problems to be considered. The adjoint equation approach developed in this thesis is applied to the computation of optimal mid-surfaces and shell thicknesses, using a large shape parameter set, proportional in size to the number of degrees of freedom in the underlying finite element discretization.
April 2010
Complex Flow and Transport Phenomena in Porous Media
Aycil Cesmelioglu
This thesis analyzes partial differential equations related to the coupled surface and subsurface flows and develops efficient high order discontinuous Galerkin (DG) methods to solve them numerically. Specifically, the coupling of the Navier-Stokes and the Darcy's equations, which is encountered in the environmental problem of groundwater contamination through lakes and rivers, is considered. Predicting accurately the transport of contaminants by this coupled flow is of great importance for the remediation strategies. The first part of this thesis analyzes a weak formulation of the time-dependent Navier-Stokes equation coupled with the Darcy's equation through the Beavers-Joseph-Saffman condition. The analysis changes depending on whether the inertial forces are included in the interface conditions or not. The inclusion of the inertial forces (Model I) remedies the difficulty in the analysis caused by the nonlinear convection term; however, it does not reflect the physical interactions on the interface correctly. Hence, I also analyze the weak problem by omitting these forces (Model II) which complicates the analysis and necessitates an extra small data condition. For Model I, a fully discrete scheme based on the DG method and the Crank-Nicolson method is introduced. The convergence of the scheme is proven with optimal error estimates. The second part couples the surface flow and a convection-diffusion type transport with miscible displacement in the subsurface. Initially, I consider the coupled stationary Stokes and Darcy's equations for the flow and establish the existence of a weak solution. Next, imposing additional assumptions on the data, I extend the result to the nonlinear case where the surface flow is given by the Navier-Stokes equation. The analysis also applies to the particular case where the flow is loosely coupled to the transport, that is, the velocity field obtained from the flow is an input for the transport equation. The flow is discretized by combinations of the continuous finite element method and the DG method whereas the discretization of the transport is done by a combined DG and backward Euler methods. The scheme yields optimal error estimates and its robustness for fractured porous media is shown by a numerical example.
April 2010
Model Reduction of Large Spiking Neurons
Anthony Richard Kellems
This thesis introduces and applies model reduction techniques to problems associated with simulation of realistic single neurons. Neurons have complicated dendritic structures and spatially-distributed ionic kinetics that give rise to highly nonlinear dynamics. However, existing model reduction methods compromise the geometry, and thus sacrifice the original input-output relationship. I demonstrate that linear and nonlinear model reduction techniques yield systems that capture the salient dynamics of morphologically accurate neuronal models and preserve the input-output maps while using significantly fewer variables than the full systems. Two main dynamic regimes characterize the voltage response of a neuron, and I demonstrate that different model reduction techniques are well-suited to each regime.
Small perturbations from the neuron's rest state fall into the subthreshold regime, which can be accurately described by a linear system. By applying Balanced Truncation (BT), a model reduction technique for general linear systems, I recover subthreshold voltage dynamics, and I provide an efficient Iterative Rational Krylov Algorithm (IRKA), which makes large problems of interest tractable. However, these approximations are not valid once the input to the neuron is sufficient to drive the voltage into the spiking regime, which is characterized by highly nonlinear behavior. To reproduce spiking dynamics, I use a proper orthogonal decomposition (POD) to reduce the number of state variables and a discrete empirical interpolation method (DEIM) to reduce the complexity of the nonlinear terms.
The techniques described above are successful, but they inherently assume that the whole neuron is either passive (linear) or active (nonlinear). However, in realistic cells the voltage response at distal locations is nearly linear, while at proximal locations it is very nonlinear. With this observation, I fuse the aforementioned models together to create a reduced coupled model in which each reduction technique is used where it is most advantageous, thereby making it possible to more accurately simulate a larger class of cortical neurons.
This thesis contains the following articles, which are used with kind permission from Springer Science+Business Media: Kellems, A.R., Chaturantabut, S., Sorensen, D.C., and Cox, S.J. (2010) Morphologically accurate reduced order modeling of spiking neurons. J. Comput. Neurosci. The original publication is available at www.springerlink.com: DOI: 10.1007/s10827-010-0229-4
Kellems, A.R., Roos, D., Xiao, N., and Cox, S.J. (2009) Low-dimensional, morphologically accurate models of subthreshold membrane potential. J. Comput. Neurosci., 27:161-176 The original publication is available at www.springerlink.com: DOI: 10.1007/s10827-008-0134-2
April 2010
Source Localization in Cluttered Acoustic Waveguides
Leila Issa
Mode coupling due to scattering by weak random inhomogeneities leads to the loss of coherence in the wave field measured a long distances of propagation. This in turn leads to the deterioration of coherent source localization methods such as matched-field. In this dissertation, we study with analysis and numerical simulations how such deterioration occurs and we introduce an original incoherent source localization approach for random waveguides. This approach is based on a special form of transport theory for the incoherent fluctuations of the wave field. The statistical stability of the method is analyzed and its performance is illustrated with numerical simulations. In addition, this method is used to estimate the correlation function of the random fluctuations of the wave speed.
April 2010
Source Synthesis for Waveform Inversion
William W. Symes
Seismologists have long used extended sources, both physical and synthetic, to enhance the effectiveness of data acquisition and processing in various ways. Recently, extended synthetic sources have been suggested as a means to substantially reduce the computational complexity of waveform inversion. Most proposed synthetic sources for waveform inversion are random in some way. In contrast, I suggest a deterministic choice of extended source, one that maximizes the difference between data predicted by the current inversion iterate and data synthesized from the target point source records. Optimal source synthesis in this sense has proven effecive in several non-seismic imaging technologies. I will describe the concept, and present some prelimiary evidence concerning its perfomance in inversion of reflection seismograms.
April 2010
Local Error Analysis of Discontinuous Galerkin Methods for Advection-Dominated Elliptic Linear-Quadratic Optimal Control Problems
Dmitriy Leykekhman & Matthias Heinkenschloss
This paper analyzes the local properties of the symmetric interior penalty upwind discontinuous Galerkin method (SIPG) for the numerical solution of optimal control problems governed by linear reaction-advection-diffusion equations with distributed controls. The theoretical and numerical results presented in this paper show that for advection-dominated problems the convergence properties of the SIPG discretization can be superior to the convergence properties of stabilized finite element discretizations such as the streamline upwind Petrov Galerkin (SUPG) method. For example, we show that for a small diffusion parameter the SIPG method is optimal in the interior of the domain. This is in sharp contrast to SUPG discretizations, for which it is known that the existence of boundary layers can pollute the numerical solution of optimal control problems everywhere even into domains where the solution is smooth and, as a consequence, in general reduces the convergence rates to only first order. In order to prove the nice convergence properties of the SIPG discretization for optimal control problems, we first improve local error estimates of the SIPG discretization for single advection-dominating equations by showing that the size of the numerical boundary layer is controlled not by the mesh size but rather by the size of the diffusion parameter. As a result, for small diffusion, the boundary layers are too "weak" to pollute the SIPG solution into domains of smoothness in optimal control problems. This favorable property of the SIPG method is due to the weak treatment of boundary conditions which is natural for discontinuous Galerkin methods, while for SUPG methods strong imposition of boundary conditions is more conventional. The importance of the weak treatment of boundary conditions for the solution of advection dominated optimal control problems with distributed controls is also supported by our numerical results.
Key Words: Optimal control, advection-diffusion equations, discontinuous Galerkin methods, discretization, local error estimates, distributed control.
March 2010
On the Coupling of Finite Volume and Discontinuous Galerkin Method for Elliptic Problems
Prince Chidyagwai, Ilya Mishev, & Beatrice Riviere
The coupling of cell-centered finite volume method with primal discontinuous Galerkin method is introduced in this paper for elliptic problems. Convergence of the method with respect to the mesh size is proved. Numerical examples confirm the theoretical rates of convergence. Advantages of the coupled scheme are shown for problems with discontinuous coefficients or anisotropic diffusion matrix.
March 2010
Continuous and Discontinuous Finite Element Methods for Coupled Surface-Subsurface Flow and Transport Problems
A. Cesmelioglu, P. Chidyagwai, & B. Riviere
A weak formulation of the coupled problem of flow and transport is discretized and analyzed numerically. The flow problem is characterized by the Navier-Stokes (or Stokes) equations coupled by Darcy equations. The velocity field is obtained by couplings of finite element and discontinuous Galerkin methods. The concentration equation is solved by an improved discontinuous Galerkin method. Convergence of the schemes is obtained. Numerical examples show the robustness of the method for heterogeneous and fractured porous media
March 2010
Efficient Numerical Methods for Least-Norm Regularization
D.C. Sorensen & M. Rojas
The problem
arises in the regularization of discrete forms of ill-posed problems when an estimate of the noise level in the data is available. After deriving necessary and sufficient optimality conditions for this problem, we propose two different classes of algorithms: a factorization-based algorithm for small to medium problems, and matrix-free iterations for the large-scale case. Numerical results illustrating the performance of the methods demonstrate that both classes of algorithms are efficient, robust, and accurate. An interesting feature of our formulation is that there is no situation corresponding to the so-called hard case that occurs in the standard trust-region subproblem. Neither small singular values nor vanishing coefficients present any difficulty to solving the relevant secular equations.
March 2010
Solving a Low-Rank Factorization Model for Matrix Completion by a Non-linear Successive Over-Relaxation Algorithm
Zaiwen Wen, Wotao Yin, & Yin Zhang
The matrix completion problem is to recover a low-rank matrix from a subset of its entries. The main solution strategy for this problem has been based on nuclear-norm minimization which requires computing singular value decompositions -- a task that is increasingly costly as matrix sizes and ranks increase. To improve the capacity of solving large-scale problems, we propose a low-rank factorization model and construct a nonlinear successive over-relaxation (SOR) algorithm that only requires solving a linear least squares problem per iteration. Convergence of this nonlinear SOR algorithm is analyzed. Numerical results show that the algorithm can reliably solve a wide range of problems at a speed at least several times faster than nuclear-norm minimization algorithms.
March 2010
IWAVE Implementation of Adjoint State Method
Dong Sun & William W Symes
Adjoint state method is a well-known method to efficiently compute the gradient of a cost or objective function for a simulation-driven optimization problem. Essentially, it computes the adjoint action of Born operator (the linearized forward map) on any given vector. This report presents a derivation of adjoint state algorithm for an acoustic system discretized by staggered grid finite difference schemes, and discusses its implementation based on the modeling package IWAVE.
Our goal is to construct a C++ wrapper of IWAVE, which fits into a general framework for inversion. This report is the second of several describing an implementation of such a wrapper.
March 2010
IWAVE Implementation of Born Simulation
Dong Sun & William W Symes
The single-scattering (or Born) approximation is the most fundamental assumption shared by all seismic imaging methods, and plays a crutial role in the non-linear waveform inversion, an iterative process of linearized inversions. The Born simulator (linearized forward map) shares a computational core with the corresponding simulator (forward map), which has been well implemented in the modeling package IWAVE. This report focuses on implementing the Born simulator based on IWAVE, and reviews the main adaptations we made in IWAVE to accommodate such an implementation in C++.
Our goal is to construct a C++ wrapper of IWAVE, which fits into a general framework for inversion. This report is the first of several describing an implementation of such a wrapper.
March 2010
Discontinuous Galerkin Time Domain Methods for Acoustics and Comparison with Finite Difference Time Domain
Methods
Xin Wang
This thesis describes an implementation of the discontinuous Galerkin finite element time domain (DGTD) method on unstructured meshes to solve acoustic wave equations in heterogeneous media. In oil industry people use finite difference time domain (FDTD) methods to simulate seismic surveys, the first step to explore oil and gas in the earth's subsurface, conducted either in land or sea. The results in this thesis indicate that the first order time shift effect resulting from misalignment between numerical meshes and material interfaces in the DGTD method occurs in the same way as interface errors in the finite difference simulation of wave propagation. This thesis describes two approaches: interface-fitting mesh and local mesh refinement, without modifying the DGTD scheme, to decrease this troublesome effect with verifications of 2D examples. The comparison of the DGTD method on the piecewise linear interface- fitting mesh and the staggered FDTD method both applied to a square-circle model and a 2D dome model in this thesis confirms the fact that the DGTD method can achieve a suboptimal second order convergence rate while the error in the staggered FDTD method is dominated by the first order interface error when the curved material interfaces are presented. I conclude that the DGTD method is more efficient than the staggered FDTD method for the two solutions to have roughly the same accuracy when the accuracy requirement becomes more and more strict and the model becomes more and more complex.
March 2010
An Alternating Direction Algorithm for Nonnegative Matrix Factorization
Yin Zhang
We extend the classic alternating direction method for convex optimization to solving the non-convex, non- negative matrix factorization problem and conduct several carefully designed numerical experiments to compare the proposed algorithms with the most widely used two algorithms for solving this problem. In addition, the proposed algorithm is also briefly compared with two other more recent algorithms. Numerical evidence shows that the alternating direction algorithm tends to deliver higher-quality solutions with faster computing times on the tested problems. A convergence result is given showing that the algorithm converges to a Karush-Kuhn- Tucker point whenever it converges.
January 2010
EdgeCS: Edge Guided Compressive Sensing Reconstruction
Weihong Guo & Wotao Yin
Compressive sensing (CS) reconstructs images from a small number of projections. We propose EdgeCS - edge guided CS reconstruction - to recover images of higher qualities from fewer measurements than the current state-of-the-art methods.
Accurate edge information can significantly improve image recovery quality and speed, but such information is encoded in the CS measurements of an image. To take advantage of edge information in CS recovery, EdgeCS alternatively performs CS reconstruction and edge detection in a way that each benefits from the latest solution of the other. EdgeCS is fast and returns high-quality images. It exactly recovers the 256 × 256 Shepp-Logan phantom from merely 7 radial lines (or 3.03% k-space), which is impossible for most existing algorithms. It accurately reconstructs a 512×512 magnetic resonance image from 21% noisy samples. Moreover, it is also able to reconstruct complex-valued images. Each took about 30 seconds on an ordinary laptop. The algorithm can be easily ported to GPUs for a speedup of more than 10 folds.
January 2010
Practical Compressive Sensing with Toeplitz and Circulant Matrices
Wotao Yin, Simon Morgan, Junfeng Yang, & Yin Zhang
Compressive sensing encodes a signal into a relatively small number of incoherent linear measurements. In theory, the optimal incoherence is achieved by completely random measurement matrices. However, such matrices are difficult and/or costly to implement in hardware realizations.
After summarizing recent results of how random Toeplitz and circulant matrices can be easily (or even naturally) realized in various applications, we introduce fast algorithms for reconstructing signals from incomplete Toeplitz and circulant measurements. We present computational results showing that Toeplitz and circulant matrices are not only as effective as random matrices for signal encoding, but also permit much faster signal decoding.
January 2010