SIAM Digital Library
 
 
 

SIAM J. on Scientific Computing

All Journal content published prior to 1997 is part of LOCUS.

Search Issue | RSS Feeds RSS
Previous Issue Next Issue

2011

Volume 33, Issue 5, pp. 2115-3086

† Special Section: 2010 Copper Mountain Conference


3D Composite Finite Elements for Elliptic Boundary Value Problems with Discontinuous Coefficients

Tobias Preusser, Martin Rumpf, Stefan Sauter, and Lars Ole Schwen

SIAM J. Sci. Comput. 33, pp. 2115-2143 (29 pages)

Online Publication Date: September 01, 2011

Full Text: | Download PDF

Show Abstract
For scalar and vector-valued elliptic boundary value problems with discontinuous coefficients across geometrically complicated interfaces, a composite finite element approach is developed. Composite basis functions are constructed, mimicking the expected jump condition for the solution at the interface in an approximate sense. The construction is based on a suitable local interpolation on the space of admissible functions. We study the order of approximation and the convergence properties of the method numerically. As applications, heat diffusion in an aluminum foam matrix filled with polymer and linear elasticity of microstructured materials, in particular, specimens of trabecular bone, are investigated. Furthermore, a numerical homogenization approach is developed for periodic structures and real material specimens which are not strictly periodic but are considered as statistical prototypes. Thereby, effective macroscopic material properties can be computed.

Methods for Pricing American Options under Regime Switching

Y. Huang, P. A. Forsyth, and G. Labahn

SIAM J. Sci. Comput. 33, pp. 2144-2168 (25 pages) | Cited 1 time

Online Publication Date: September 01, 2011

Full Text: | Download PDF

Show Abstract
We analyze a number of techniques for pricing American options under a regime switching stochastic process. The techniques analyzed include both explicit and implicit discretizations with the focus being on methods which are unconditionally stable. In the case of implicit methods we also compare a number of iterative procedures for solving the associated nonlinear algebraic equations. Numerical tests indicate that a fixed point policy iteration, coupled with a direct control formulation, is a reliable general purpose method. Finally, we remark that we formulate the American problem as an abstract optimal control problem; hence our results are applicable to more general problems as well.

An Online Method for Interpolating Linear Parametric Reduced-Order Models

David Amsallem and Charbel Farhat

SIAM J. Sci. Comput. 33, pp. 2169-2198 (30 pages)

Online Publication Date: September 01, 2011

Full Text: | Download PDF

Show Abstract
A two-step online method is proposed for interpolating projection-based linear parametric reduced-order models (ROMs) in order to construct a new ROM for a new set of parameter values. The first step of this method transforms each precomputed ROM into a consistent set of generalized coordinates. The second step interpolates the associated linear operators on their appropriate matrix manifold. Real-time performance is achieved by precomputing inner products between the reduced-order bases underlying the precomputed ROMs. The proposed method is illustrated by applications in mechanical and aeronautical engineering. In particular, its robustness is demonstrated by its ability to handle the case where the sampled parameter set values exhibit a mode veering phenomenon.

Duality Based A Posteriori Error Estimation for Quasi-Periodic Solutions Using Time Averages

M. Braack, E. Burman, and N. Taschenberger

SIAM J. Sci. Comput. 33, pp. 2199-2216 (18 pages)

Online Publication Date: September 06, 2011

Full Text: | Download PDF

Show Abstract
We propose an a posteriori error estimation technique for the computation of average functionals of solutions for nonlinear time dependent problems based on duality techniques. The exact solution is assumed to have a periodic or quasi-periodic behavior favoring a fixed mesh strategy in time. We show how to circumvent the need of solving time dependent dual problems. The estimator consists of an averaged residual weighted by sensitivity factors coming from a stationary dual problem and an additional averaging error term coming from nonlinearities of the operator considered. In order to illustrate this technique the resulting adaptive algorithm is applied to several model problems: a linear scalar parabolic problem with known exact solution, the nonsteady Navier–Stokes equations with known exact solution, and finally to the well-known benchmark problem for Navier–Stokes (flow behind a cylinder) in order to verify the modeling assumptions.

Time Discretization Schemes for Poincaré Waves in Finite-Element Shallow-Water Models

Daniel Y. Le Roux, Michel Dieme, and Abdou Sene

SIAM J. Sci. Comput. 33, pp. 2217-2246 (30 pages)

Online Publication Date: September 13, 2011

Full Text: | Download PDF

Show Abstract
The finite-element spatial discretization of the linear shallow-water equations is examined in the context of several temporal discretization schemes. Three finite-element pairs are considered, namely, the $P^{}_{0}-P^{}_{1}$, $P^{NC}_{1}-P^{}_{1}$, and $RT^{}_{0}-P^{}_{0}$ schemes, and the backward and forward Euler, Crank–Nicolson, and second and third order Adams–Bashforth time stepping schemes are employed. A Fourier analysis is performed at the discrete level for the Poincaré waves, and it determines the stability limit of the schemes and the error in wave amplitude and phase that can be expected. Numerical solutions of test problems to simulate Poincaré waves illustrate the analytical results.

Nonadiabatic Transitions Through Tilted Avoided Crossings

Volker Betz and Benjamin D. Goddard

SIAM J. Sci. Comput. 33, pp. 2247-2276 (30 pages)

Online Publication Date: September 13, 2011

Full Text: | Download PDF

Show Abstract
We investigate the transition of a quantum wave-packet through a one-dimensional avoided crossing of molecular energy levels when the energy levels at the crossing point are tilted. Using superadiabatic representations, and an approximation of the dynamics near the crossing region, we obtain an explicit formula for the transition wave function. Our results agree extremely well with high precision ab-initio calculations.

A Higher Order Scheme for a Tangentially Stabilized Plane Curve Shortening Flow with a Driving Force

Martin Balažovjech and Karol Mikula

SIAM J. Sci. Comput. 33, pp. 2277-2294 (18 pages)

Online Publication Date: September 15, 2011

Full Text: | Download PDF

Show Abstract
We introduce a new higher order scheme for computing a tangentially stabilized curve shortening flow with a driving force represented by an intrinsic partial differential equation for an evolving curve position vector. Our new scheme is a combination of the explicit forward Euler and the fully implicit backward Euler schemes. At any discrete time step, the solution is found efficiently using a few semi-implicit iterations. Basic properties of the new scheme are proved in the paper, and its precision is tested by comparing the results with known analytical solutions. For any choice of the time step, the new higher order scheme gives exact radius of evolving uniformly discretized circles in case of flow by curvature and in case of rotation by a constant tangential velocity. Such properties do not hold for other schemes solving flow by mean curvature like the classical explicit, semi-implicit, or fully implicit schemes. In general, the scheme is second order accurate, which is shown by comparing a numerically evolving encompassed area with known analytical expression. The behavior of the scheme is discussed on representative examples, and its advantages with respect to the balance between CPU time and precision are shown.

Tensor-Train Decomposition

I. V. Oseledets

SIAM J. Sci. Comput. 33, pp. 2295-2317 (23 pages) | Cited 3 times

Online Publication Date: September 22, 2011

Full Text: | Download PDF

Show Abstract
A simple nonrecursive form of the tensor decomposition in $d$ dimensions is presented. It does not inherently suffer from the curse of dimensionality, it has asymptotically the same number of parameters as the canonical decomposition, but it is stable and its computation is based on low-rank approximation of auxiliary unfolding matrices. The new form gives a clear and convenient way to implement all basic operations efficiently. A fast rounding procedure is presented, as well as basic linear algebra operations. Examples showing the benefits of the decomposition are given, and the efficiency is demonstrated by the computation of the smallest eigenvalue of a 19-dimensional operator.

A General Framework for Deriving Integral Preserving Numerical Methods for PDEs

Morten Dahlby and Brynjulf Owren

SIAM J. Sci. Comput. 33, pp. 2318-2340 (23 pages)

Online Publication Date: September 22, 2011

Full Text: | Download PDF

Show Abstract
A general procedure for constructing conservative numerical integrators for time-dependent partial differential equations is presented. In particular, linearly implicit methods preserving a time discretized version of the invariant are developed for systems of partial differential equations with polynomial nonlinearities. The framework is rather general and allows for an arbitrary number of dependent and independent variables with derivatives of any order. It is proved formally that second order convergence is obtained. The procedure is applied to a test case, and numerical experiments are provided.

Fast Evaluation of Multiquadric RBF Sums by a Cartesian Treecode

Robert Krasny and Lei Wang

SIAM J. Sci. Comput. 33, pp. 2341-2355 (15 pages) | Cited 1 time

Online Publication Date: September 22, 2011

Full Text: | Download PDF

Show Abstract
A treecode is presented for evaluating sums defined in terms of the multiquadric radial basis function (RBF), $\phi({\bf x}) = (|{\bf x}|^2+c^2)^{1/2}$, where ${\bf x} \in \mathbb{R}^3$ and $c \ge 0$. Given a set of $N$ nodes, evaluating an RBF sum directly requires CPU time that scales like $O(N^2)$. For a given level of accuracy, the treecode reduces the CPU time to $O(N\log N)$ using a far-field expansion of $\phi({\bf x})$. We consider two options for the far-field expansion: (1) a Laurent series previously used in applications of the Fast Multipole Method to multiquadric RBFs, and (2) a certain Taylor series previously used in treecode particle simulations, but not yet in the context of multiquadric RBFs. It is known that the Laurent series converges when the RBF parameter $c$ lies in an interval $0 \le c \le \bar{c}$, where $\bar{c}$ is proportional to the minimum node spacing, but here we show that the Taylor series converges uniformly for $c \ge 0$. We implement the treecode in Cartesian coordinates and use a recurrence relation to compute the Taylor coefficients. Numerical results exhibit the treecode error, CPU time, and memory usage in two test cases, random nodes in a cube and on the surface of a sphere. The treecode approach presented here is applicable to generalized multiquadrics in any dimension.

A Well-Balanced Symplecticity-Preserving Gas-Kinetic Scheme for Hydrodynamic Equations under Gravitational Field

Jun Luo, Kun Xu, and Na Liu

SIAM J. Sci. Comput. 33, pp. 2356-2381 (26 pages)

Online Publication Date: September 27, 2011

Full Text: | Download PDF

Show Abstract
A well-balanced scheme for an isolated gravitational hydrodynamic system is defined as a scheme which exactly preserves an isothermal hydrostatic solution. In this paper, a well-balanced gas-kinetic symplecticity-preserving BGK (SP-BGK) scheme is developed. In the construction of the scheme, the gravitational potential is modeled as a piecewise constant function inside each cell with a potential jump at the cell interface. In the process of designing such a scheme, the energy conservation, Liouville's theorem, and the symplecticity-preserving property of a Hamiltonian flow play important roles in the description of particles penetration and reflection from a potential barrier. More importantly, the use of the symplecticity-preserving property is crucial in the evaluation of the moments of a postinteraction gas distribution function with a potential jump in terms of the moments of preinteraction distribution function. The SP-BGK method is the first well-balanced shock-capturing gas-kinetic scheme for the Navier–Stokes equation. A few theorems are proved for this scheme, which include the necessity to use an exact Maxwellian for keeping the isothermal hydrostatic state, the total mass and energy (the sum of kinetic, thermal, and gravitational ones) conservation, and the well-balanced property of the SP-BGK scheme to keep an isothermal hydrostatic state during the process of particle transport and collision. Many numerical examples are presented to validate the SP-BGK scheme.

Algorithms for Area Preserving Flows

Catherine Kublik, Selim Esedoḡlu, and Jeffrey A. Fessler

SIAM J. Sci. Comput. 33, pp. 2382-2401 (20 pages)

Online Publication Date: September 27, 2011

Full Text: | Download PDF

Show Abstract
We propose efficient and accurate algorithms for computing certain area preserving geometric motions of curves in the plane, such as area preserving motion by curvature. These schemes are based on a new class of diffusion generated motion algorithms using signed distance functions. In particular, they alternate two very simple and fast operations, namely convolution with the Gaussian kernel and construction of the distance function, to generate the desired geometric flow in an unconditionally stable manner. We present applications of these area preserving flows to large scale simulations of coarsening.

A Third Order Accurate Fast Marching Method for the Eikonal Equation in Two Dimensions

Shahnawaz Ahmed, Stanley Bak, Joyce McLaughlin, and Daniel Renzi

SIAM J. Sci. Comput. 33, pp. 2402-2420 (19 pages)

Online Publication Date: September 29, 2011

Full Text: | Download PDF

Show Abstract
In this paper, we develop a third order accurate fast marching method for the solution of the eikonal equation in two dimensions. There have been two obstacles to extending the fast marching method to higher orders of accuracy. The first obstacle is that using one-sided difference schemes is unstable for orders of accuracy higher than two. The second obstacle is that the points in the difference stencil are not available when the gradient is closely aligned with the grid. We overcome these obstacles by using a two-dimensional (2D) finite difference approximation to improve stability, and by locally rotating the grid 45 degrees (i.e., using derivatives along the diagonals) to ensure all the points needed in the difference stencil are available. We show that in smooth regions the full difference stencil is used for a suitably small enough grid size and that the difference scheme satisfies the von Neumann stability condition for the linearized eikonal equation. Our method reverts to first order accuracy near caustics without developing oscillations by using a simple switching scheme. The efficiency and high order of the method are demonstrated on a number of 2D test problems.

Stable and Spectrally Accurate Schemes for the Navier–Stokes Equations

Jun Jia and Jie Liu

SIAM J. Sci. Comput. 33, pp. 2421-2439 (19 pages)

Online Publication Date: September 29, 2011

Full Text: | Download PDF

Show Abstract
In this paper, we present an accurate, efficient and stable numerical method for the incompressible Navier–Stokes equations (NSEs). The method is based on (1) an equivalent pressure Poisson equation formulation of the NSE with proper pressure boundary conditions, which facilitates the design of high-order and stable numerical methods, and (2) the Krylov deferred correction (KDC) accelerated method of lines transpose ($\mbox{MoL}^{\it T}$), which is very stable, efficient, and of arbitrary order in time. Numerical tests with known exact solutions in three dimensions show that the new method is spectrally accurate in time, and a numerical order of convergence 9 was observed. Two-dimensional computational results of flow past a cylinder and flow in a bifurcated tube are also reported.

Schur Complement Preconditioners for Surface Integral-Equation Formulations of Dielectric Problems Solved with the Multilevel Fast Multipole Algorithm

Tahi̇r Malas and Levent Gürel

SIAM J. Sci. Comput. 33, pp. 2440-2467 (28 pages)

Online Publication Date: October 04, 2011

Full Text: | Download PDF

Show Abstract
Surface integral-equation methods accelerated with the multilevel fast multipole algorithm (MLFMA) provide a suitable mechanism for electromagnetic analysis of real-life dielectric problems. Unlike the perfect-electric-conductor case, discretizations of surface formulations of dielectric problems yield $2 \times 2$ partitioned linear systems. Among various surface formulations, the combined tangential formulation (CTF) is the closest to the category of first-kind integral equations, and hence it yields the most accurate results, particularly when the dielectric constant is high and/or the dielectric problem involves sharp edges and corners. However, matrix equations of CTF are highly ill-conditioned, and their iterative solutions require powerful preconditioners for convergence. Second-kind surface integral-equation formulations yield better conditioned systems, but their conditionings significantly degrade when real-life problems include high dielectric constants. In this paper, for the first time in the context of surface integral-equation methods of dielectric objects, we propose Schur complement preconditioners to increase their robustness and efficiency. First, we approximate the dense system matrix by a sparse near-field matrix, which is formed naturally by MLFMA. The Schur complement preconditioning requires approximate solutions of systems involving the (1,1) partition and the Schur complement. We approximate the inverse of the (1,1) partition with a sparse approximate inverse (SAI) based on the Frobenius norm minimization. For the Schur complement, we first approximate it via incomplete sparse matrix-matrix multiplications, and then we generate its approximate inverse with the same SAI technique. Numerical experiments on sphere, lens, and photonic crystal problems demonstrate the effectiveness of the proposed preconditioners. In particular, the results for the photonic crystal problem, which has both surface singularity and a high dielectric constant, shows that accurate CTF solutions for such problems can be obtained even faster than with second-kind integral equation formulations, with the acceleration provided by the proposed Schur complement preconditioners.

A Fast Iterative Method for Solving the Eikonal Equation on Triangulated Surfaces

Zhisong Fu, Won-Ki Jeong, Yongsheng Pan, Robert M. Kirby, and Ross T. Whitaker

SIAM J. Sci. Comput. 33, pp. 2468-2488 (21 pages)

Online Publication Date: October 06, 2011

Full Text: | Download PDF

Show Abstract
This paper presents an efficient, fine-grained parallel algorithm for solving the Eikonal equation on triangular meshes. The Eikonal equation, and the broader class of Hamilton–Jacobi equations to which it belongs, have a wide range of applications from geometric optics and seismology to biological modeling and analysis of geometry and images. The ability to solve such equations accurately and efficiently provides new capabilities for exploring and visualizing parameter spaces and for solving inverse problems that rely on such equations in the forward model. Efficient solvers on state-of-the-art, parallel architectures require new algorithms that are not, in many cases, optimal, but are better suited to synchronous updates of the solution. In previous work [W. K. Jeong and R. T. Whitaker, SIAM J. Sci. Comput., 30 (2008), pp. 2512–2534], the authors proposed the fast iterative method (FIM) to efficiently solve the Eikonal equation on regular grids. In this paper we extend the fast iterative method to solve Eikonal equations efficiently on triangulated domains on the CPU and on parallel architectures, including graphics processors. We propose a new local update scheme that provides solutions of first-order accuracy for both architectures. We also propose a novel triangle-based update scheme and its corresponding data structure for efficient irregular data mapping to parallel single-instruction multiple-data (SIMD) processors. We provide detailed descriptions of the implementations on a single CPU, a multicore CPU with shared memory, and SIMD architectures with comparative results against state-of-the-art Eikonal solvers.

Interpolatory Projection Methods for Parameterized Model Reduction

Ulrike Baur, Christopher Beattie, Peter Benner, and Serkan Gugercin

SIAM J. Sci. Comput. 33, pp. 2489-2518 (30 pages)

Online Publication Date: October 13, 2011

Full Text: | Download PDF

Show Abstract
We provide a unifying projection-based framework for structure-preserving interpolatory model reduction of parameterized linear dynamical systems, i.e., systems having a structured dependence on parameters that we wish to retain in the reduced-order model. The parameter dependence may be linear or nonlinear and is retained in the reduced-order model. Moreover, we are able to give conditions under which the gradient and Hessian of the system response with respect to the system parameters is matched in the reduced-order model. We provide a systematic approach built on established interpolatory $\mathcal{H}_2$ optimal model reduction methods that will produce parameterized reduced-order models having high fidelity throughout a parameter range of interest. For single input/single output systems with parameters in the input/output maps, we provide reduced-order models that are optimal with respect to an $\mathcal{H}_2\otimes\mathcal{L}_2$ joint error measure. The capabilities of these approaches are illustrated by several numerical examples from technical applications.

A Sinc Function Analogue of Chebfun

Mark Richardson and Lloyd N. Trefethen

SIAM J. Sci. Comput. 33, pp. 2519-2535 (17 pages)

Online Publication Date: October 13, 2011

Full Text: | Download PDF

Show Abstract
Chebfun is an established software system for computing with functions of a real variable, but its capabilities for handling functions with singularities are limited. Here an analogous system is described based on sinc function expansions instead of Chebyshev series. This experiment sheds light on the strengths and weaknesses of sinc function techniques. It also serves as a review of some of the main features of sinc methods, including construction, evaluation, zerofinding, optimization, integration, and differentiation.

Geometric Properties of the Icosahedral-Hexagonal Grid on the Two-Sphere

Ning Wang and Jin-Luen Lee

SIAM J. Sci. Comput. 33, pp. 2536-2559 (24 pages)

Online Publication Date: October 18, 2011

Full Text: | Download PDF

Show Abstract
An icosahedral-hexagonal grid on the two-sphere is created by dividing the faces of an icosahedron and projecting the vertices onto the sphere. This grid and its Voronoi tessellation have several desirable features for numerical simulations of physical processes on the sphere. While several methods to construct the icosahedral grid mesh have been proposed over the past decades, and empirical data have been collected to understand and help improve the grid, rarely have analytical analyses been done to investigate the basic geometric properties of the grid. In this paper, we present an analytical analysis of several geometric properties of the icosahedral grids based on two basic constructions: recursive and nonrecursive construction. We point out that these geometric properties can be improved with modified construction procedures. We demonstrate how these improvements impact the numerical integration of PDEs over the sphere.

A Comparison of the Dispersion and Dissipation Errors of Gauss and Gauss–Lobatto Discontinuous Galerkin Spectral Element Methods

Gregor Gassner and David A. Kopriva

SIAM J. Sci. Comput. 33, pp. 2560-2579 (20 pages)

Online Publication Date: October 18, 2011

Full Text: | Download PDF

Show Abstract
We examine the dispersion and dissipation properties of the Gauss and Gauss–Lobatto discontinuous Galerkin spectral element methods (DGSEMs), for linear wave propagation problems. We show that the inherent underintegration in the Gauss–Lobatto variant can be interpreted as a modal filtering of the highest polynomial mode. This in turn has a drastic impact on the dispersion and dissipation relations of the Gauss–Lobatto DGSEM compared to the Gauss variant. We show that the Gauss DGSEM is typically more accurate than the Gauss–Lobatto variant, needing fewer points per wavelength for a given accuracy while on the other hand being more restricted in the explicit time step choice. We show that the spectra of the DGSEM operators depend on the boundary conditions applied and that the ratio of the time step restrictions of the two schemes depends on the choice of boundary conditions.

An Algorithm for the Principal Component Analysis of Large Data Sets

Nathan Halko, Per-Gunnar Martinsson, Yoel Shkolnisky, and Mark Tygert

SIAM J. Sci. Comput. 33, pp. 2580-2594 (15 pages) | Cited 1 time

Online Publication Date: October 18, 2011

Full Text: | Download PDF

Show Abstract
Recently popularized randomized methods for principal component analysis (PCA) efficiently and reliably produce nearly optimal accuracy—even on parallel processors—unlike the classical (deterministic) alternatives. We adapt one of these randomized methods for use with data sets that are too large to be stored in random-access memory (RAM). (The traditional terminology is that our procedure works efficiently out-of-core.) We illustrate the performance of the algorithm via several numerical examples. For example, we report on the PCA of a data set stored on disk that is so large that less than a hundredth of it can fit in our computer's RAM.

Nonsymmetric Preconditioner Updates in Newton–Krylov Methods for Nonlinear Systems

Stefania Bellavia, Daniele Bertaccini, and Benedetta Morini

SIAM J. Sci. Comput. 33, pp. 2595-2619 (25 pages)

Online Publication Date: October 20, 2011

Full Text: | Download PDF

Show Abstract
Newton–Krylov methods, a combination of Newton-like methods and Krylov subspace methods for solving the Newton equations, often need adequate preconditioning in order to be successful. Approximations of the Jacobian matrices are required to form preconditioners, and this step is very often the dominant cost of Newton–Krylov methods. Therefore, working with preconditioners may destroy the “Jacobian-free” (or matrix-free) setting where the single Jacobian-vector product can be provided without forming and storing the element of the true Jacobian. In this paper, we propose and analyze a preconditioning technique for sequences of nonsymmetric Jacobian matrices based on the update of an earlier preconditioner. The proposed strategy can be implemented in a matrix-free manner. Numerical experiments on popular test problems confirm the effectiveness of the approach in comparison with the standard ILU-preconditioned Newton–Krylov approaches.

Interpolating Irregularly Spaced Observations for Filtering Turbulent Complex Systems

John Harlim

SIAM J. Sci. Comput. 33, pp. 2620-2640 (21 pages)

Online Publication Date: October 20, 2011

Full Text: | Download PDF

Show Abstract
We present a numerically fast reduced filtering strategy, the Fourier domain Kalman filter with appropriate interpolations to account for irregularly spaced observations of complex turbulent signals. The design of such a reduced filter involves: (i) interpolating irregularly spaced observations to the model regularly spaced grid points, (ii) understanding under which situation the small scale oscillatory artifact from such interpolation won't degrade the filtered solutions, (iii) understanding when the interpolated covariance structure can be approximated by its diagonal terms when observations are corrupted by independent Gaussian noise, and (iv) applying a scalar Kalman filter formula on each Fourier component independently with an approximate diagonal interpolated covariance matrix. From the practical point of view, there is an emerging need to understand the effect of (i) toward the filtered solutions, for example, in utilizing the data produced from various interpolation techniques that merge multiple satellite measurements of atmospheric and ocean dynamical quantities. To understand point (iii) above and to see how many of nondiagonal terms are effectively ignored in (iv), we compute a ratio $\Lambda$ between the largest nondiagonal components and the smallest diagonal components of the interpolated covariance matrix. We find that for piecewise linear interpolation, this ratio, $\Lambda$, is always smaller than that of the alternative interpolation schemes such as trigonometric and cubic spline for any irregularly spaced observation networks. When observations are not so sparse with small noise, we find that the small scale oscillatory artifact in (ii) above is negligible when piecewise linear interpolation is used whereas for the other schemes such as the nearest neighbor, trigonometric, and cubic spline interpolation, the oscillatory artifact degrades the filtered solutions significantly. Finally, we also find that the reduced filtering strategy with piecewise linear interpolation produces more accurate filtered solutions than conventional approaches when observations are extremely irregularly spaced (such that the ratio $\Lambda$ is not so small) and very sparse.

Optimal Control of Stochastic Flow over a Backward-Facing Step Using Reduced-Order Modeling

Max Gunzburger and Ju Ming

SIAM J. Sci. Comput. 33, pp. 2641-2663 (23 pages)

Online Publication Date: October 20, 2011

Full Text: | Download PDF

Show Abstract
We explore a stochastic optimal control problem for incompressible flow in a channel by using reduced-order modeling (ROM) based on proper orthogonal decomposition (POD). The control acts on a portion of the boundary and has a stochastic perturbation. The objective of control is to minimize the enstrophy, i.e., the $L^2$ norm of the vorticity. Owing to the high computational requirements, the classical Monte Carlo method is not an appropriate option. The POD basis is constructed in the usual manner from a set of “snapshots” obtained by sampling the solution of the Navier–Stokes system at several time instants and over several independent realizations of the random data. Then, the stochastic ROM is constructed via a Galerkin method with respect to the POD basis. The optimal control problem is solved using the stochastic ROM and compared to results obtained using a full finite element discretization that employs many more degrees of freedom. These comparisons show the effectiveness of the use of the stochastic ROM for obtaining accurate optimal solutions of the stochastic control problem at low cost.

Superconvergence of Discontinuous Galerkin Solutions for Delay Differential Equations of Pantograph Type

Qiumei Huang, Hehu Xie, and Hermann Brunner

SIAM J. Sci. Comput. 33, pp. 2664-2684 (21 pages)

Online Publication Date: October 25, 2011

Full Text: | Download PDF

Show Abstract
This paper is concerned with the superconvergence of the discontinuous Galerkin solutions for delay differential equations with proportional delays vanishing at $t = 0$. Two types of superconvergence are analyzed here. The first is based on interpolation postprocessing to improve the global convergence order by finding the superconvergence points of discontinuous Galerkin solutions. The second type follows from the integral iteration which just requires a local integration procedure applied to the discontinuous Galerkin solution, thus increasing the order of convergence. The theoretical results are illustrated by a broad range of numerical examples.
back to top
RSS Feeds
FREE

Special Section: 2010 Copper Mountain Conference

Ray Tuminaro, Guest Editor , Michele Benzi, Xiao-Chuan Cai, Iain Duff, Howard Elman, Roland Freund, Kirk Jordan, Tim Kelley, David Keyes, Misha Kilmer, Sven Leyffer, Tom Manteuffel, Steve McCormick, David Silvester, Homer Walker, et al.

SIAM J. Sci. Comput. 33, pp. 2685-2685 (1 page)

Full Text: | Download PDF

Show Abstract
The biennial Copper Mountain Conference on Iterative Methods was held April 4–9, 2010. This meeting included more than 140 presentations covering many scientific computing areas, such as uncertainty quantification, optimization, Markov chains, saddle-point systems, inverse problems, direct factorizations, Krylov methods, algebraic multigrid, software frameworks, and advanced computer architectures. A partial list of the represented application domains included acoustic scattering, material science, circuit simulation, neutron transport, image processing, magnetohydrodynamics, and porous media flow.
Howard Elman (University of Maryland) and Ray Tuminaro (Sandia National Laboratories) served as co-chairs. Management services were provided by Annette Anthony of Front Range Scientific Computations, Inc. The Conference is gratefully for financial support from the Department of Energy, IBM, Lawrence Livermore National Laboratory, Los Alamos National Laboratory, Sandia National Laboratories, and the National Science Foundation. The Copper Mountain Conference is organized in cooperation with the Society for Industrial and Applied Mathematics (SIAM). Submissions to this special section were open to the public. This was advertised in advance on web sites of the SIAM Journal on Scientific Computing (SISC) and the Copper Mountain Conference.
The collection of papers in this section all center on iterative methods. However, many diverse topics are covered spanning contributions with a theoretical orientation to those that are more applied or driven by a specific application area. Additionally, several papers are highly influenced by current computer architecture trends. Some of the submissions also appeared in the student paper competition, a unique event at the Copper Mountain conferences that has proven to be a successful forum for presenting the latest results in the field. We have had increasingly strong participation from talented recent Ph.D. holders and current Ph.D. students in this competition.
I wish to thank the guest editorial board, the SIAM/SISC staff, and the referees for their hard work to help meet deadlines. Special acknowledgment is due to Mitch Chernoff (SIAM Publications Manager) and Heather Blythe (SIAM Senior Publications Coordinator) for their efforts on this special section.
We hope that this special section will continue to foster new developments in this important area of research. The diversity of topics clearly demonstrates the significant impact of iterative methods across many scientific computing areas and the need for continued research.

Estimates of the Norm of the Error in Solving Linear Systems with FOM and GMRES

Gérard Meurant

SIAM J. Sci. Comput. 33, pp. 2686-2705 (20 pages)

Online Publication Date: October 27, 2011

Full Text: | Download PDF

Show Abstract
We provide formulas for the norm of the error when solving nonsymmetric linear systems with the full orthogonalization method (FOM) and the generalized minimum residual method (GMRES) as well as relations between the error norm and the residual norm. From these formulas we are able to compute estimates of the norm of the error during the iterations. Since stopping criteria based on the norm of the residual may sometimes be misleading, such estimates could lead to a more robust way to stop the iterations. Numerical experiments show that the proposed norm estimates work nicely on difficult linear systems.

Algebraic Multigrid for Linear Systems Obtained by Explicit Element Reduction

Thomas A. Brunner and Tzanio V. Kolev

SIAM J. Sci. Comput. 33, pp. 2706-2731 (26 pages)

Online Publication Date: October 27, 2011

Full Text: | Download PDF

Show Abstract
We consider sparse linear systems, where a set of “interior” degrees of freedom have been eliminated in order to reduce the problem size. This elimination is assumed to be local, so the “interior” principal submatrix is block-diagonal, and the resulting Schur complement is still sparse. For it to be beneficial, the elimination process should lead to reduced memory requirements, but equally importantly, it should also result in an algebraic problem that can be solved efficiently. In this paper we propose a general element reduction approach and show how the elimination process can exploit a particular “subzonal” discretization to maintain the sparsity of the Schur complement. We also investigate algebraic multigrid (AMG) solution algorithms applied to the reduced problem, and we discuss the influence of the local elimination on solver-related properties of the matrix, such as near-nullspace preservation and the availability of stable subspace decompositions. We focus on BoomerAMG, a parallel variant of classical Ruge–Stüben AMG, applied to scalar diffusion problems [V. Henson and U. Yang, Appl. Numer. Math., 41 (2002), pp. 155–177], and the auxiliary-space Maxwell solver (AMS) for electromagnetic diffusion applications [T. Kolev and P. Vassilevski, J. Comput. Math., 27 (2009), pp. 604–623]. In the electromagnetic case, we establish algebraically a reduced version of the Hiptmair–Xu decomposition from [R. Hiptmair and J. Xu, SIAM J. Numer. Anal., 45 (2007), pp. 2483–2509] and consider a modification of the reduction process that targets the singular problems arising in simulations with pure void (zero conductivity) regions. For scalar diffusion problems, our particular stencil analysis shows that the reduction has a positive effect on meshes with stretched elements. We present a number of two-dimensional, three-dimensional, and axisymmetric numerical experiments, which demonstrate that the combination of an appropriately chosen local elimination with the use of the BoomerAMG and AMS solvers can lead to significant improvements in the overall solution time.

Peano—A Traversal and Storage Scheme for Octree-Like Adaptive Cartesian Multiscale Grids

Tobias Weinzierl and Miriam Mehl

SIAM J. Sci. Comput. 33, pp. 2732-2760 (29 pages)

Online Publication Date: October 27, 2011

Full Text: | Download PDF

Show Abstract
Almost all approaches to solving partial differential equations (PDEs) are based upon a spatial discretization of the computational domain—a grid. This paper presents an algorithm to generate, store, and traverse a hierarchy of $d$-dimensional Cartesian grids represented by a $(k=3)$-spacetree, a generalization of the well-known octree concept, and it also shows the correctness of the approach. These grids may change their adaptive structure throughout the traversal. The algorithm uses $2d+4$ stacks as data structures for both cells and vertices, and the storage requirements for the pure grid reduce to one bit per vertex for both the complete grid connectivity structure and the multilevel grid relations. Since the traversal algorithm uses only stacks, the algorithm's cache hit rate is continually higher than 99.9 percent, and the runtime per vertex remains almost constant; i.e., it does not depend on the overall number of vertices or the adaptivity pattern. We use the algorithmic approach as the fundamental concept for a mesh management for $d$-dimensional PDEs and for a matrix-free PDE solver represented by a compact discrete $3^d$-point operator. In the latter case, one can implement a Jacobi smoother, a Krylov solver, or a geometric multigrid scheme within the presented traversal scheme which inherits the low memory requirements and the good memory access characteristics directly.

Analysis of Augmented Lagrangian-Based Preconditioners for the Steady Incompressible Navier–Stokes Equations

Michele Benzi and Zhen Wang

SIAM J. Sci. Comput. 33, pp. 2761-2784 (24 pages)

Online Publication Date: October 27, 2011

Full Text: | Download PDF

Show Abstract
We analyze a class of modified augmented Lagrangian-based preconditioners for both stable and stabilized finite element discretizations of the steady incompressible Navier–Stokes equations. We study the eigenvalues of the preconditioned matrices obtained from Picard linearization, and we devise a simple and effective method for the choice of the augmentation parameter $\gamma$ based on Fourier analysis. Numerical experiments on a wide range of model problems demonstrate the robustness of these preconditioners, yielding fast convergence independent of mesh size and only mildly dependent on viscosity on both uniform and stretched grids. Good results are also obtained on linear systems arising from Newton linearization. We also show that performing inexact preconditioner solves with an algebraic multigrid algorithm results in excellent scalability. Comparisons of the modified augmented Lagrangian preconditioners with other state-of-the-art techniques show the competitiveness of our approach.

Iterative Methods for Neutron Transport Eigenvalue Problems

Fynn Scheben and Ivan G. Graham

SIAM J. Sci. Comput. 33, pp. 2785-2804 (20 pages)

Online Publication Date: October 27, 2011

Full Text: | Download PDF

Show Abstract
We discuss iterative methods for computing criticality in nuclear reactors. In general this requires the solution of a generalized eigenvalue problem for an unsymmetric integro-differential operator in six independent variables, modeling transport, scattering, and fission, where the dependent variable is the neutron angular flux. In engineering practice this problem is often solved iteratively, using some variant of the inverse power method. Because of the high dimension, matrix representations for the operators are often not available and the inner solves needed for the eigenvalue iteration are implemented by matrix-free inner iterations. This leads to technically complicated inexact iterative methods, for which there appears to be no published rigorous convergence theory. For the monoenergetic homogeneous model case with isotropic scattering and vacuum boundary conditions, we show that, before discretization, the general nonsymmetric eigenproblem for the angular flux is equivalent to a certain related eigenproblem for the scalar flux, involving a symmetric positive definite weakly singular integral operator (in space only). This correspondence to a symmetric problem (in a space of reduced dimension) permits us to give a convergence theory for inexact inverse iteration and related methods. In particular this theory provides rather precise criteria on how accurate the inner solves need to be in order for the whole iterative method to converge. We also give examples of discretizations which have a corresponding symmetric finite-dimensional reduced form. The theory is illustrated with numerical examples for several test problems of physical relevance, using GMRES as the inner solver.

A Second Order Discretization of Maxwell's Equations in the Quasi-Static Regime on OcTree Grids

Lior Horesh and Eldad Haber

SIAM J. Sci. Comput. 33, pp. 2805-2822 (18 pages)

Online Publication Date: October 27, 2011

Full Text: | Download PDF

Show Abstract
In this study we consider adaptive mesh refinement for the solution of Maxwell's equations in the quasi-static or diffusion regime. We propose a new finite volume OcTree discretization for the problem and show how to construct second order stencils on Yee grids, extending the known first order discretization stencils. We then develop an effective preconditioner to the problem. We show that our preconditioner performs well for discontinuous conductivities as well as for a wide range of frequencies.

Domain-Decomposition-Type Methods for Computing the Diagonal of a Matrix Inverse

Jok M. Tang and Yousef Saad

SIAM J. Sci. Comput. 33, pp. 2823-2847 (25 pages)

Online Publication Date: October 27, 2011

Full Text: | Download PDF

Show Abstract
This paper presents two methods based on domain decomposition concepts for determining the diagonal of the inverse of specific matrices. The first uses a divide-and-conquer principle and the Sherman–Morrison–Woodbury formula and assumes that the matrix can be decomposed into a $2 \times 2$ block-diagonal matrix and a low-rank matrix. The second method is a standard domain decomposition approach in which local solves are combined with a global correction. Both methods can be successfully combined with iterative solvers and sparse approximation techniques. The efficiency of the methods usually depends on the specific implementation, which should be fine-tuned for different test problems. Preliminary results for some two-dimensional (2D) problems are reported to illustrate the proposed methods.

Efficient Parallel Nonnegative Least Squares on Multicore Architectures

Yuancheng Luo and Ramani Duraiswami

SIAM J. Sci. Comput. 33, pp. 2848-2863 (16 pages)

Online Publication Date: October 27, 2011

Full Text: | Download PDF

Show Abstract
We parallelize a version of the active-set iterative algorithm derived from the original works of Lawson and Hanson [Solving Least Squares Problems, Prentice-Hall, 1974] on multicore architectures. This algorithm requires the solution of an unconstrained least squares problem in every step of the iteration for a matrix composed of the passive columns of the original system matrix. To achieve improved performance, we use parallelizable procedures to efficiently update and downdate the $QR$ factorization of the matrix at each iteration, to account for inserted and removed columns. We use a reordering strategy of the columns in the decomposition to reduce computation and memory access costs. We consider graphics processing units (GPUs) as a new mode for efficient parallel computations and compare our implementations to that of multicore CPUs. Both synthetic and nonsynthetic data are used in the experiments.

Multigrid Smoothers for Ultraparallel Computing

Allison H. Baker, Robert D. Falgout, Tzanio V. Kolev, and Ulrike Meier Yang

SIAM J. Sci. Comput. 33, pp. 2864-2887 (24 pages) | Cited 1 time

Online Publication Date: October 27, 2011

Full Text: | Download PDF

Show Abstract
This paper investigates the properties of smoothers in the context of algebraic multigrid (AMG) running on parallel computers with potentially millions of processors. The development of multigrid smoothers in this case is challenging, because some of the best relaxation schemes, such as the Gauss–Seidel (GS) algorithm, are inherently sequential. Based on the sharp two-grid multigrid theory from [R. D. Falgout and P. S. Vassilevski, SIAM J. Numer. Anal., 42 (2004), pp. 1669–1693] and [R. D. Falgout, P. S. Vassilevski, and L. T. Zikatanov, Numer. Linear Algebra Appl., 12 (2005), pp. 471–494] we characterize the smoothing properties of a number of practical candidates for parallel smoothers, including several $C$-$F$, polynomial, and hybrid schemes. We show, in particular, that the popular hybrid GS algorithm has multigrid smoothing properties which are independent of the number of processors in many practical applications, provided that the problem size per processor is large enough. This is encouraging news for the scalability of AMG on ultraparallel computers. We also introduce the more robust $\ell_1$ smoothers, which are always convergent and have already proven essential for the parallel solution of some electromagnetic problems [T. Kolev and P. Vassilevski, J. Comput. Math., 27 (2009), pp. 604–623].

Algebraic Multigrid for High-Order Hierarchical $H(curl)$ Finite Elements

James H. Lai and Luke N. Olson

SIAM J. Sci. Comput. 33, pp. 2888-2902 (15 pages)

Online Publication Date: October 27, 2011

Full Text: | Download PDF

Show Abstract
Classic multigrid methods are often not directly applicable to nonelliptic problems such as curl-type partial differential equations (PDEs). Curl-curl PDEs require specialized smoothers that are compatible with the gradient-like (near) null space. Moreover, recent developments have focused on replicating the grad-curl-div de Rham complex in a multilevel hierarchy through smoothed aggregation based algebraic multigrid. These approaches have been successful for Nédélec finite elements (i.e., $H(curl)$ edge elements), but do not extend naturally to high-order representations. In this paper we consider hierarchical high-order Whitney elements for the curl-curl eddy current problem and devise a scalable multilevel approach. Our method generates a hierarchy similar to $p$-type multigrid, which results in a coarse level that is amenable to further coarsening through the established process of a multilevel complex. The natural hierarchy of the elements induces an effective interpolation operator and motivates the construction of a compatible gradient smoothing process. We detail the multilevel solver for a hierarchical $H(curl)$ basis and present numerical results in support of the method.

Preconditioning Iterative Methods for the Optimal Control of the Stokes Equations

Tyrone Rees and Andrew J. Wathen

SIAM J. Sci. Comput. 33, pp. 2903-2926 (24 pages)

Online Publication Date: October 27, 2011

Full Text: | Download PDF

Show Abstract
Solving problems regarding the optimal control of partial differential equations (PDEs)—also known as PDE-constrained optimization—is a frontier area of numerical analysis. Of particular interest is the problem of flow control, where one would like to effect some desired flow by exerting, for example, an external force. The bottleneck in many current algorithms is the solution of the optimality system—a system of equations in saddle point form that is usually very large and ill conditioned. In this paper we describe two preconditioners—a block diagonal preconditioner for the minimal residual method and a block lower-triangular preconditioner for a nonstandard conjugate gradient method—which can be effective when applied to such problems where the PDEs are the Stokes equations. We consider only distributed control here, although we believe other problems could be treated in the same way. We give numerical results, and we compare these with those obtained by solving the equivalent forward problem using similar techniques.

On-the-Fly Adaptive Smoothed Aggregation Multigrid for Markov Chains

Eran Treister and Irad Yavneh

SIAM J. Sci. Comput. 33, pp. 2927-2949 (23 pages)

Online Publication Date: October 27, 2011

Full Text: | Download PDF

Show Abstract
A new adaptive algebraic multigrid scheme is developed for the solution of Markov chains, where the hierarchy of operators is adapted on-the-fly in a setup process that is interlaced with the solution process. The setup process feeds the solution process with improved operators, while the solution process provides the adaptive setup process with better approximations on which to base further-improved operators. The approach is demonstrated using Petrov–Galerkin smoothed aggregation where only the prolongation operator is smoothed, while the restriction remains of low order. Results show that the on-the-fly adaptive scheme can improve the performance of multigrid solvers that require extensive setup computations, in both serial and parallel environments.

LSMR: An Iterative Algorithm for Sparse Least-Squares Problems

David Chin-Lung Fong and Michael Saunders

SIAM J. Sci. Comput. 33, pp. 2950-2971 (22 pages)

Online Publication Date: October 27, 2011

Full Text: | Download PDF

Show Abstract
An iterative method LSMR is presented for solving linear systems $Ax=b$ and least-squares problems $\min \|Ax-b\|_2$, with $A$ being sparse or a fast linear operator. LSMR is based on the Golub–Kahan bidiagonalization process. It is analytically equivalent to the MINRES method applied to the normal equation $A^T\! Ax = A^T\! b$, so that the quantities $\|A^T\! r_k\|$ are monotonically decreasing (where $r_k = b - Ax_k$ is the residual for the current iterate $x_k$). We observe in practice that $\|r_k\|$ also decreases monotonically, so that compared to LSQR (for which only $\|r_k\|$ is monotonic) it is safer to terminate LSMR early. We also report some experiments with reorthogonalization.

Iterative Parameter-Choice and Multigrid Methods for Anisotropic Diffusion Denoising

D. Chen, S. MacLachlan, and M. Kilmer

SIAM J. Sci. Comput. 33, pp. 2972-2994 (23 pages)

Online Publication Date: October 27, 2011

Full Text: | Download PDF

Show Abstract
Anisotropic diffusion methods are well known for giving good qualitative results for image denoising. This paper gives a review of the anisotropic diffusion methodology and its application to image denoising. We propose a fixed-point iteration using a multigrid solver to solve a regularized anisotropic diffusion equation, which is not only well-posed, but also has a nontrivial steady-state solution. A new regularization parameter-choice method (Brent-NCP), combining Brent's method and the normalized cumulative periodogram information of the misfit, is also introduced. We test our algorithm on several common test images with different noise levels. The experimental results demonstrate the effectiveness of the anisotropic diffusion with a multigrid approach and the broad applicability of the Brent-NCP parameter-choice algorithm.

A Factorization of the Spectral Galerkin System for Parameterized Matrix Equations: Derivation and Applications

Paul G. Constantine, David F. Gleich, and Gianluca Iaccarino

SIAM J. Sci. Comput. 33, pp. 2995-3009 (15 pages)

Online Publication Date: October 27, 2011

Full Text: | Download PDF

Show Abstract
Recent work has explored solver strategies for the linear system of equations arising from a spectral Galerkin approximation of the solution of PDEs with parameterized (or stochastic) inputs. We consider the related problem of a matrix equation whose matrix and right-hand side depend on a set of parameters (e.g., a PDE with stochastic inputs semidiscretized in space) and examine the linear system arising from a similar Galerkin approximation of the solution. We derive a useful factorization of this system of equations, which yields bounds on the eigenvalues, clues to preconditioning, and a flexible implementation method for a wide array of problems. We complement this analysis with (i) a numerical study of preconditioners on a standard elliptic PDE test problem and (ii) a fluids application using existing CFD codes; the MATLAB codes used in the numerical studies are available online.

Accelerating the Explicitly Restarted Arnoldi Method with GPUs Using an Autotuned Matrix Vector Product

Jérôme Dubois, Christophe Calvin, and Serge Petiton

SIAM J. Sci. Comput. 33, pp. 3010-3019 (10 pages)

Online Publication Date: October 27, 2011

Full Text: | Download PDF

Show Abstract
This paper presents a parallelized hybrid single-vector Arnoldi algorithm for computing approximations to eigenpairs of a nonsymmetric matrix. We are interested in the use of accelerators and multicore units to speed up the Arnoldi process. The main goal is to propose a parallel version of the Arnoldi solver, which can efficiently use multiple multicore processors or multiple graphics processing units (GPUs) in a mixed coarse and fine grain fashion. In the proposed algorithms, this is achieved by an autotuning of the matrix vector product before starting the Arnoldi eigensolver as well as the reorganization of the data and global communications so that communication time is reduced. The execution time, performance, and scalability are assessed with well-known dense and sparse test matrices on multiple Nehalems, GT200 NVidia Tesla, and next generation Fermi Tesla. With one processor, we see a performance speedup of 2 to 3x when using all the physical cores, and a total speedup of 2 to 8x when adding a GPU to this multicore unit, and hence a speedup of 4 to 24x compared to the sequential solver.

The BiCOR and CORS Iterative Algorithms for Solving Nonsymmetric Linear Systems

B. Carpentieri, Y.-F. Jing, and T.-Z. Huang

SIAM J. Sci. Comput. 33, pp. 3020-3036 (17 pages)

Online Publication Date: October 27, 2011

Full Text: | Download PDF

Show Abstract
We present two iterative algorithms for solving real nonsymmetric and complex non-Hermitian linear systems of equations and that were developed from variants of the nonsymmetric Lanczos method. In this paper, we give the theoretical background of the two iterative methods and discuss their main computational aspects. Using a large number of numerical experiments, we analyze their convergence properties, and we also compare them with other popular nonsymmetric iterative solvers in use today.

Restarting the Nonsymmetric Lanczos Algorithm for Eigenvalues and Linear Equations Including Multiple Right-Hand Sides

Ronald B. Morgan and Dywayne A. Nicely

SIAM J. Sci. Comput. 33, pp. 3037-3056 (20 pages)

Online Publication Date: October 27, 2011

Full Text: | Download PDF

Show Abstract
A restarted nonsymmetric Lanczos algorithm is given for computing eigenvalues and both right and left eigenvectors. The restarting limits the storage so that finding eigenvectors is practical. Restarting also makes it possible to deal with roundoff error in new ways. We give a scheme for avoiding near-breakdown and discuss maintaining biorthogonality. A system of linear equations can be solved simultaneously with the eigenvalue computations. Deflation from the presence of the eigenvectors allows the linear equations to generally have good convergence in spite of the restarting. The right and left eigenvectors generated while solving the linear equations can be used to help solve systems with multiple right-hand sides.

A Regularized Gauss–Newton Trust Region Approach to Imaging in Diffuse Optical Tomography

Eric de Sturler and Misha E. Kilmer

SIAM J. Sci. Comput. 33, pp. 3057-3086 (30 pages)

Online Publication Date: October 27, 2011

Full Text: | Download PDF

Show Abstract
We present a new algorithm for the solution of nonlinear least squares problems arising from parameterized imaging problems with diffuse optical tomographic data [D. Boas et al., IEEE Signal Process. Mag., 18 (2001), pp. 57–75]. The parameterization arises from the use of parametric level sets for regularization [M. E. Kilmer et al., Proc. SPIE, 5559 (2004), pp. 381–391], [A. Aghasi, M. E. Kilmer, and E. L. Miller, SIAM J. Imaging Sci., 4 (2011), pp. 618–650]. Such problems lead to Jacobians that have relatively few columns, a relatively modest number of rows, and are ill-conditioned. Moreover, such problems have function and Jacobian evaluations that are computationally expensive. Our optimization algorithm is appropriate for any inverse or imaging problem with those characteristics. In fact, we expect our algorithm to be effective for more general problems with ill-conditioned Jacobians. The algorithm aims to minimize the total number of function and Jacobian evaluations by analyzing which spectral components of the Gauss–Newton direction should be discarded or damped. The analysis considers for each component the reduction of the objective function and the contribution to the search direction, restricting the computed search direction to be within a trust region. The result is a truncated SVD-like approach to choosing the search direction. However, we do not necessarily discard components in order of decreasing singular value, and some components may be scaled to maintain fidelity to the trust region model. Our algorithm uses the Basic Trust Region Algorithm from [A. R. Conn, N. I. M. Gould, and Ph. L. Toint, Trust-Region Methods, SIAM, Philadelphia, 2000]. We prove that our algorithm is globally convergent to a critical point. Our numerical results show that the new algorithm generally outperforms competing methods applied to the DOT imaging problem with parametric level sets, and regularly does so by a significant factor.
Close

close