Showing papers in "arXiv: Numerical Analysis in 2009"

Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions

Abstract: Low-rank matrix approximations, such as the truncated singular value decomposition and the rank-revealing QR decomposition, play a central role in data analysis and scientific computing. This work surveys and extends recent research which demonstrates that randomization offers a powerful tool for performing low-rank matrix approximation. These techniques exploit modern computational architectures more fully than classical methods and open the possibility of dealing with truly massive data sets. This paper presents a modular framework for constructing randomized algorithms that compute partial matrix decompositions. These methods use random sampling to identify a subspace that captures most of the action of a matrix. The input matrix is then compressed---either explicitly or implicitly---to this subspace, and the reduced matrix is manipulated deterministically to obtain the desired low-rank factorization. In many cases, this approach beats its classical competitors in terms of accuracy, speed, and robustness. These claims are supported by extensive numerical experiments and a detailed error analysis.

Finite element exterior calculus: from Hodge theory to numerical stability

Abstract: This article reports on the confluence of two streams of research, one emanating from the fields of numerical analysis and scientific computation, the other from topology and geometry. In it we consider the numerical discretization of partial differential equations that are related to differential complexes so that de Rham cohomology and Hodge theory are key tools for the continuous problem. After a brief introduction to finite element methods, the discretization methods we consider, we develop an abstract Hilbert space framework for analyzing stability and convergence. In this framework, the differential complex is represented by a complex of Hilbert spaces and stability is obtained by transferring Hodge theoretic structures from the continuous level to the discrete. We show stable discretization is achieved if the finite element spaces satisfy two hypotheses: they form a subcomplex and there exists a bounded cochain projection from the full complex to the subcomplex. Next, we consider the most canonical example of the abstract theory, in which the Hilbert complex is the de Rham complex of a domain in Euclidean space. We use the Koszul complex to construct two families of finite element differential forms, show that these can be arranged in subcomplexes of the de Rham complex in numerous ways, and for each construct a bounded cochain projection. The abstract theory therefore applies to give the stability and convergence of finite element approximations of the Hodge Laplacian. Other applications are considered as well, especially to the equations of elasticity. Background material is included to make the presentation self-contained for a variety of readers.

ADMiRA: Atomic Decomposition for Minimum Rank Approximation

Abstract: We address the inverse problem that arises in compressed sensing of a low-rank matrix. Our approach is to pose the inverse problem as an approximation problem with a specified target rank of the solution. A simple search over the target rank then provides the minimum rank solution satisfying a prescribed data approximation bound. We propose an atomic decomposition that provides an analogy between parsimonious representations of a sparse vector and a low-rank matrix. Efficient greedy algorithms to solve the inverse problem for the vector case are extended to the matrix case through this atomic decomposition. In particular, we propose an efficient and guaranteed algorithm named ADMiRA that extends CoSaMP, its analogue for the vector case. The performance guarantee is given in terms of the rank-restricted isometry property and bounds both the number of iterations and the error in the approximate solution for the general case where the solution is approximately low-rank and the measurements are noisy. With a sparse measurement operator such as the one arising in the matrix completion problem, the computation in ADMiRA is linear in the number of measurements. The numerical experiments for the matrix completion problem show that, although the measurement operator in this case does not satisfy the rank-restricted isometry property, ADMiRA is a competitive algorithm for matrix completion.

A Krylov subspace algorithm for evaluating the phi-functions appearing in exponential integrators

Abstract: We develop an algorithm for computing the solution of a large system of linear ordinary dierential equations (ODEs) with polynomial in- homogeneity. This is equivalent to computing the action of a certain matrix function on the vector representing the initial condition. The matrix function is a linear combination of the matrix exponential and other functions related to the exponential (the so-called '-functions). Such computations are the ma- jor computational burden in the implementation of exponential integrators, which can solve general ODEs. Our approach is to compute the action of the matrix function by constructing a Krylov subspace using Arnoldi or Lanczos iteration and projecting the function on this subspace. This is combined with time-stepping to prevent the Krylov subspace from growing too large. The algorithm is fully adaptive: it varies both the size of the time steps and the dimension of the Krylov subspace to reach the required accuracy with mini- mal eort. We implement this algorithm in the matlab function phipm and we give instructions on how to obtain and use this function. Various numerical experiments show that the phipm function is often significantly more ecient than the state-of-the-art.

Stable and Accurate Interpolation Operators for High-Order Multi-Block Finite-Difference Methods

Abstract: Block-to-block interface interpolation operators are constructed for several common high-order finite difference discretizations. In contrast to conventional interpolation operators, these new interpolation operators maintain the strict stability, accuracy and conservation of the base scheme even when nonconforming grids or dissimilar operators are used in adjoining blocks. The stability properties of the new operators are verified using eigenvalue analysis, and the accuracy properties are verified using numerical simulations of the Euler equations in two spatial dimensions.

Greedy Solution of Ill-Posed Problems: Error Bounds and Exact Inversion

Abstract: The orthogonal matching pursuit (OMP) is an algorithm to solve sparse approximation problems. Sufficient conditions for exact recovery are known with and without noise. In this paper we investigate the applicability of the OMP for the solution of ill-posed inverse problems in general and in particular for two deconvolution examples from mass spectrometry and digital holography respectively. In sparse approximation problems one often has to deal with the problem of redundancy of a dictionary, i.e. the atoms are not linearly independent. However, one expects them to be approximatively orthogonal and this is quantified by the so-called incoherence. This idea cannot be transfered to ill-posed inverse problems since here the atoms are typically far from orthogonal: The ill-posedness of the operator causes that the correlation of two distinct atoms probably gets huge, i.e. that two atoms can look much alike. Therefore one needs conditions which take the structure of the problem into account and work without the concept of coherence. In this paper we develop results for exact recovery of the support of noisy signals. In the two examples in mass spectrometry and digital holography we show that our results lead to practically relevant estimates such that one may check a priori if the experimental setup guarantees exact deconvolution with OMP. Especially in the example from digital holography our analysis may be regarded as a first step to calculate the resolution power of droplet holography.

Semi-Lagrangian schemes for linear and fully non-linear diffusion equations

Abstract: For linear and fully non-linear diffusion equations of Bellman-Isaacs type, we introduce a class of approximation schemes based on differencing and interpolation. As opposed to classical numerical methods, these schemes work for general diffusions with coefficient matrices that may be non-diagonal dominant and arbitrarily degenerate. In general such schemes have to have a wide stencil. Besides providing a unifying framework for several known first order accurate schemes, our class of schemes includes new first and higher order versions. The methods are easy to implement and more efficient than some other known schemes. We prove consistency and stability of the methods, and for the monotone first order methods, we prove convergence in the general case and robust error estimates in the convex case. The methods are extensively tested.

Multi-scale methods for wave propagation in heterogeneous media

Abstract: Multi-scale wave propagation problems are computationally costly to solve by traditional techniques because the smallest scales must be represented over a domain determined by the largest scales of the problem. We have developed and analyzed new numerical methods for multi-scale wave propagation in the framework of heterogeneous multi-scale method. The numerical methods couples simulations on macro- and micro-scales for problems with rapidly oscillating coefficients. We show that the complexity of the new method is significantly lower than that of traditional techniques with a computational cost that is essentially independent of the micro-scale. A convergence proof is given and numerical results are presented for periodic problems in one, two and three dimensions. The method is also successfully applied to non-periodic problems and for long time integration where dispersive effects occur.

Non-intrusive and structure preserving multiscale integration of stiff ODEs, SDEs and Hamiltonian systems with hidden slow dynamics via flow averaging

Abstract: We introduce a new class of integrators for stiff ODEs as well as SDEs. These integrators are (i) {\it Multiscale}: they are based on flow averaging and so do not fully resolve the fast variables and have a computational cost determined by slow variables (ii) {\it Versatile}: the method is based on averaging the flows of the given dynamical system (which may have hidden slow and fast processes) instead of averaging the instantaneous drift of assumed separated slow and fast processes. This bypasses the need for identifying explicitly (or numerically) the slow or fast variables (iii) {\it Nonintrusive}: A pre-existing numerical scheme resolving the microscopic time scale can be used as a black box and easily turned into one of the integrators in this paper by turning the large coefficients on over a microscopic timescale and off during a mesoscopic timescale (iv) {\it Convergent over two scales}: strongly over slow processes and in the sense of measures over fast ones. We introduce the related notion of two-scale flow convergence and analyze the convergence of these integrators under the induced topology (v) {\it Structure preserving}: for stiff Hamiltonian systems (possibly on manifolds), they can be made to be symplectic, time-reversible, and symmetry preserving (symmetries are group actions that leave the system invariant) in all variables. They are explicit and applicable to arbitrary stiff potentials (that need not be quadratic). Their application to the Fermi-Pasta-Ulam problems shows accuracy and stability over four orders of magnitude of time scales. For stiff Langevin equations, they are symmetry preserving, time-reversible and Boltzmann-Gibbs reversible, quasi-symplectic on all variables and conformally symplectic with isotropic friction.

Variational Theory and Domain Decomposition for Nonlocal Problems

Abstract: In this article we present the first results on domain decomposition methods for nonlocal operators. We present a nonlocal variational formulation for these operators and establish the well-posedness of associated boundary value problems, proving a nonlocal Poincare inequality. To determine the conditioning of the discretized operator, we prove a spectral equivalence which leads to a mesh size independent upper bound for the condition number of the stiffness matrix. We then introduce a nonlocal two-domain variational formulation utilizing nonlocal transmission conditions, and prove equivalence with the single-domain formulation. A nonlocal Schur complement is introduced. We establish condition number bounds for the nonlocal stiffness and Schur complement matrices. Supporting numerical experiments demonstrating the conditioning of the nonlocal one- and two-domain problems are presented.

An Analysis of broken $P_1$-Nonconforming Finite Element Method For Interface Problems

Abstract: We study some numerical methods for solving second order elliptic problem with interface. We introduce an immersed interface finite element method based on the `broken' $P_1$-nonconforming piecewise linear polynomials on interface triangular elements having edge averages as degrees of freedom. This linear polynomials are broken to match the homogeneous jump condition along the interface which is allowed to cut through the element. We prove optimal orders of convergence in $H^1$ and $L^2$-norm. Next we propose a mixed finite volume method in the context introduced in \cite{Kwak2003} using the Raviart-Thomas mixed finite element and this `broken' $P_1$-nonconforming element. The advantage of this mixed finite volume method is that once we solve the symmetric positive definite pressure equation(without Lagrangian multiplier), the velocity can be computed locally by a simple formula. This procedure avoids solving the saddle point problem. Furthermore, we show optimal error estimates of velocity and pressure in our mixed finite volume method. Numerical results show optimal orders of error in $L^2$-norm and broken $H^1$-norm for the pressure, and in $H(\Div)$-norm for the velocity.

Topics in compressed sensing

Abstract: Compressed sensing has a wide range of applications that include error correction, imaging, radar and many more. Given a sparse signal in a high dimensional space, one wishes to reconstruct that signal accurately and efficiently from a number of linear measurements much less than its actual dimension. Although in theory it is clear that this is possible, the difficulty lies in the construction of algorithms that perform the recovery efficiently, as well as determining which kind of linear measurements allow for the reconstruction. There have been two distinct major approaches to sparse recovery that each present different benefits and shortcomings. The first, L1-minimization methods such as Basis Pursuit, use a linear optimization problem to recover the signal. This method provides strong guarantees and stability, but relies on Linear Programming, whose methods do not yet have strong polynomially bounded runtimes. The second approach uses greedy methods that compute the support of the signal iteratively. These methods are usually much faster than Basis Pursuit, but until recently had not been able to provide the same guarantees. This gap between the two approaches was bridged when we developed and analyzed the greedy algorithm Regularized Orthogonal Matching Pursuit (ROMP). ROMP provides similar guarantees to Basis Pursuit as well as the speed of a greedy algorithm. Our more recent algorithm Compressive Sampling Matching Pursuit (CoSaMP) improves upon these guarantees, and is optimal in every important aspect.

A modified phase field approximation for mean curvature flow with conservation of the volume

Abstract: This paper is concerned with the motion of a time dependent hypersurface that evolves by mean curvature flow with a a volume constraint. Phase field approximation of this motion leads to the well known nonlocal Allen--Cahn equation. Here we propose a modified version of this equation, and we show that it has better volume preserving properties than the classical one, even in the presence of an additional forcing term.

A Simple Proof that Random Matrices are Democratic

Abstract: The recently introduced theory of compressive sensing (CS) enables the reconstruction of sparse or compressible signals from a small set of nonadaptive, linear measurements. If properly chosen, the number of measurements can be significantly smaller than the ambient dimension of the signal and yet preserve the significant signal information. Interestingly, it can be shown that random measurement schemes provide a near-optimal encoding in terms of the required number of measurements. In this report, we explore another relatively unexplored, though often alluded to, advantage of using random matrices to acquire CS measurements. Specifically, we show that random matrices are democractic, meaning that each measurement carries roughly the same amount of signal information. We demonstrate that by slightly increasing the number of measurements, the system is robust to the loss of a small number of arbitrary measurements. In addition, we draw connections to oversampling and demonstrate stability from the loss of significantly more measurements.

Pathwise Accuracy and Ergodicity of Metropolized Integrators for SDEs

Abstract: Metropolized integrators for ergodic stochastic differential equations (SDE) are proposed which (i) are ergodic with respect to the (known) equilibrium distribution of the SDE and (ii) approximate pathwise the solutions of the SDE on finite time intervals. Both these properties are demonstrated in the paper and precise strong error estimates are obtained. It is also shown that the Metropolized integrator retains these properties even in situations where the drift in the SDE is nonglobally Lipschitz, and vanilla explicit integrators for SDEs typically become unstable and fail to be ergodic.

A Priori and A Posteriori Analysis of the Quasi-Nonlocal Quasicontinuum Method in 1D

Abstract: For a next-nearest neighbour pair interaction model in a periodic domain, a priori and a posteriori analyses of the quasinonlocal quasicontinuum method (QNL-QC) are presented. The results are valid for large deformations and essentially guarantee a one-to-one correspondence between atomistic solutions and QNL-QC solutions. The analysis is based on truncation error and residual estimates in negative norms and novel a priori and a posteriori stability estimates.

Existence and approximation of a (regularized) Oldroyd-B model

Abstract: Two finite element approximations of the Oldroyd-B model for dilute polymeric fluids are considered, in bounded 2- and 3-dimensional domains, under no flow boundary conditions. The pressure and the symmetric conformation tensor are aproximated by either (a) piecewise constants or (b) continuous piecewise linears, the velocity by (a) continuous piecewise quadratics or a reduced version with linear tangential component on each edge, and (b) by continuous piecewise quadratics or the mini-element. Both schemes (a) and (b) satisfy a free energy bound, which involves the logarithm of the conformation tensor, without any constraint on the time step for the backward Euler type time discretization. This extends the results of [Boyaval et al. M2AN 43 (2009) 523--561], where a piecewise constant approximation of the conformation tensor was necessary to treat the advection term in the stress equation, and a restriction on the time step, based on the initial data, was required to ensure that the approximation to the conformation tensor remained positive definite. Furthermore, for (b) in the presence of an additional dissipative term in the stress equation and a cut-off on the conformation tensor on certain terms like in [Barrett and Suli, M3AS 18 (2008) 935--971] for the FENE dumbbell model, we show (subsequence) convergence towards global-in-time weak solutions (when d=2, cut-offs can be replaced with a time step restriction dependent on the spatial discretization parameter). Hence, we prove existence of global-in-time weak solutions to these regularized models.

An Optimal Order Error Analysis of the One-Dimensional Quasicontinuum Approximation

Abstract: We derive a model problem for quasicontinuum approximations that allows a simple, yet insightful, analysis of the optimal-order convergence rate in the continuum limit for both the energy-based quasicontinuum approximation and the quasi-nonlocal quasicontinuum approximation. The optimal-order error estimates for the quasi-nonlocal quasicontinuum approximation are given for all strains up to the continuum limit strain for fracture. The analysis is based on an explicit treatment of the coupling error at the atomistic to continuum interface, combined with an analysis of the error due to atomistic and continuum schemes using the stability of the quasicontinuum approximation.

Reconstruction of the equilibrium of the plasma in a Tokamak and identification of the current density profile in real time

Abstract: The reconstruction of the equilibrium of a plasma in a Tokamak is a free boundary problem described by the Grad-Shafranov equation in axisymmetric configuration. The right-hand side of this equation is a nonlinear source, which represents the toroidal component of the plasma current density. This paper deals with the identification of this nonlinearity source from experimental measurements in real time. The proposed method is based on a fixed point algorithm, a finite element resolution, a reduced basis method and a least-square optimization formulation. This is implemented in a software called Equinox with which several numerical experiments are conducted to explore the identification problem. It is shown that the identification of the profile of the averaged current density and of the safety factor as a function of the poloidal flux is very robust.

An asymptotic preserving scheme for strongly anisotropic elliptic problems

Abstract: In this article we introduce an asymptotic preserving scheme designed to compute the solution of a two dimensional elliptic equation presenting large anisotropies. We focus on an anisotropy aligned with one direction, the dominant part of the elliptic operator being supplemented with Neumann boundary conditions. A new scheme is introduced which allows an accurate resolution of this elliptic equation for an arbitrary anisotropy ratio.

A Convergent Overlapping Domain Decomposition Method for Total Variation Minimization

Abstract: This paper is concerned with the analysis of convergent sequential and parallel overlapping domain decomposition methods for the minimization of functionals formed by a discrepancy term with respect to data and a total variation constraint. To our knowledge, this is the first successful attempt of addressing such strategy for the nonlinear, nonadditive, and nonsmooth problem of total variation minimization. We provide several numerical experiments, showing the successful application of the algorithm for the restoration of 1D signals and 2D images in interpolation/inpainting problems respectively, and in a compressed sensing problem, for recovering piecewise constant medical-type images from partial Fourier ensembles.

Surface Comparison with Mass Transportation

Abstract: We use mass-transportation as a tool to compare surfaces (2-manifolds) In particular, we determine the "similarity" of two given surfaces by solving a mass-transportation problem between their conformal densities This mass transportation problem differs from the standard case in that we require the solution to be invariant under global M\"obius transformations Our approach provides a constructive way of defining a metric in the abstract space of simply-connected smooth surfaces with boundary (ie surfaces of disk-type); this metric can also be used to define meaningful intrinsic distances between pairs of "patches" in the two surfaces, which allows automatic alignment of the surfaces We provide numerical experiments on "real-life" surfaces to demonstrate possible applications in natural sciences

Error bounds for computing the expectation by Markov chain Monte Carlo

Abstract: We study the error of reversible Markov chain Monte Carlo methods for approximating the expectation of a function. Explicit error bounds with respect to different norms of the function are proven. By the estimation the well known asymptotical limit of the error is attained, i.e. there is no gap between the estimate and the asymptotical behavior. We discuss the dependence of the error on a burn-in of the Markov chain. Furthermore we suggest and justify a specific burn-in for optimizing the algorithm.

Certified numerical homotopy tracking

Abstract: Given a homotopy connecting two polynomial systems we provide a rigorous algorithm for tracking a regular homotopy path connecting an approximate zero of the start system to an approximate zero of the target system. Our method uses recent results on the complexity of homotopy continuation rooted in the alpha theory of Smale. Experimental results obtained with the implementation in the numerical algebraic geometry package of Macaulay2 demonstrate the practicality of the algorithm. In particular, we confirm the theoretical results for random linear homotopies and illustrate the plausibility of a conjecture by Shub and Smale on a good initial pair.

Accelerating gradient projection methods for $\ell_1$-constrained signal recovery by steplength selection rules

Abstract: We propose a new gradient projection algorithm that compares favorably with the fastest algorithms available to date for $\ell_1$-constrained sparse recovery from noisy data, both in the compressed sensing and inverse problem frameworks. The method exploits a line-search along the feasible direction and an adaptive steplength selection based on recent strategies for the alternation of the well-known Barzilai-Borwein rules. The convergence of the proposed approach is discussed and a computational study on both well-conditioned and ill-conditioned problems is carried out for performance evaluations in comparison with five other algorithms proposed in the literature.

Modified energy for split-step methods applied to the linear Schr\"odinger equation

Abstract: We consider the linear Schr\"odinger equation and its discretization by split-step methods where the part corresponding to the Laplace operator is approximated by the midpoint rule. We show that the numerical solution coincides with the exact solution of a modified partial differential equation at each time step. This shows the existence of a modified energy preserved by the numerical scheme. This energy is close to the exact energy if the numerical solution is smooth. As a consequence, we give uniform regularity estimates for the numerical solution over arbitrary long time

Pole-based approximation of Fermi-Dirac function

Abstract: Two approaches for the efficient rational approximation of the Fermi-Dirac function are discussed: one uses the contour integral representation and conformal mapping and the other is based on a version of the multipole representation of the Fermi-Dirac function that uses only simple poles. Both representations have logarithmic computational complexity. They are of great interest for electronic structure calculations.

Iterative Methods for the Force-based Quasicontinuum Approximation

Abstract: Force-based atomistic-continuum hybrid methods are the only known pointwise consistent methods for coupling a general atomistic model to a finite element continuum model. For this reason, and due to their algorithmic simplicity, force-based coupling methods have become a popular class of atomistic-continuum hybrid models as well as other types of multiphysics models. However, the recently discovered unusual stability properties of the linearized force-based quasicontinuum (QCF) approximation, especially its indefiniteness, present a challenge to the development of efficient and reliable iterative methods. We present analytic and computational results for the generalized minimal residual (GMRES) solution of the linearized QCF equilibrium equations. We show that the GMRES method accurately reproduces the stability of the force-based approximation and conclude that an appropriately preconditioned GMRES method results in a reliable and efficient solution method.

Strong and weak divergence in finite time of Euler's method for stochastic differential equations with non-globally Lipschitz continuous coefficients

Abstract: The stochastic Euler scheme is known to converge to the exact solution of a stochastic differential equation with globally Lipschitz continuous drift and diffusion coefficient. Recent results extend this convergence to coefficients which grow at most linearly. For superlinearly growing coefficients finite-time convergence in the strong mean square sense remained an open question according to [Higham, Mao & Stuart (2002); Strong convergence of Euler-type methods for nonlinear stochastic differential equations, SIAM J. Numer. Anal. 40, no. 3, 1041-1063]. In this article we answer this question to the negative and prove for a large class of stochastic differential equations with non-globally Lipschitz continuous coefficients that Euler's approximation converges neither in the strong mean square sense nor in the numerically weak sense to the exact solution at a finite time point. Even worse, the difference of the exact solution and of the numerical approximation at a finite time point diverges to infinity in the strong mean square sense and in the numerically weak sense.

Hopf algebras of formal diffeomorphisms and numerical integration on manifolds

Abstract: B-series originated from the work of John Butcher in the 1960s as a tool to analyze numerical integration of differential equations, in particular Runge-Kutta methods. Connections to renormalization theory in perturbative quantum field theory have been established in recent years. The algebraic structure of classical Runge-Kutta methods is described by the Connes-Kreimer Hopf algebra. Lie-Butcher theory is a generalization of B-series aimed at studying Lie-group integrators for differential equations evolving on manifolds. Lie-group integrators are based on general Lie group actions on a manifold, and classical Runge-Kutta integrators appear in this setting as the special case of R^n acting upon itself by translations. Lie--Butcher theory combines classical B-series on R^n with Lie-series on manifolds. The underlying Hopf algebra combines the Connes-Kreimer Hopf algebra with the shuffle Hopf algebra of free Lie algebras. We give an introduction to Hopf algebraic structures and their relationship to structures appearing in numerical analysis, aimed at a general mathematical audience. In particular we explore the close connection between Lie series, time-dependent Lie series and Lie--Butcher series for diffeomorphisms on manifolds. The role of the Euler and Dynkin idempotents in numerical analysis is discussed. A non-commutative version of a Faa di Bruno bialgebra is introduced, and the relation to non-commutative Bell polynomials is explored.