Showing papers in "Mathematics of Computation in 2009"
TL;DR: This paper analyzes the convergence of linearized Bregman iterations and derives a new algorithm that is proven to be convergent with a rate and can be used as another choice of an efficient tool in compressed sensing.
Abstract: Finding a solution of a linear equation Au = f with various minimization properties arises from many applications. One such application is compressed sensing, where an efficient and robust-to-noise algorithm to find a minimal l 1 norm solution is needed. This means that the algorithm should be tailored for large scale and completely dense matrices A, while Au and A T u can be computed by fast transforms and the solution we seek is sparse. Recently, a simple and fast algorithm based on linearized Bregman iteration was proposed in [28, 32] for this purpose. This paper is to analyze the convergence of linearized Bregman iterations and the minimization properties of their limit. Based on our analysis here, we derive also a new algorithm that is proven to be convergent with a rate. Furthermore, the new algorithm is simple and fast in approximating a minimal l 1 norm solution of Au = f as shown by numerical simulations. Hence, it can be used as another choice of an efficient tool in compressed sensing.
327 citations
TL;DR: In this paper, a projection method for the numerical evaluation of Fredholm determinants is proposed, which is derived from the classical Nystrom method for solving Fredholm equations of the second kind, using Gauss-Legendre or Clenshaw-Curtis as the underlying quadrature rule.
Abstract: Some significant quantities in mathematics and physics are most naturally expressed as the Fredholm determinant of an integral operator, most notably many of the distribution functions in random matrix theory. Though their numerical values are of interest, there is no systematic numerical treatment of Fredholm determinants to be found in the literature. Instead, the few numerical evaluations that are available rely on eigenfunction expansions of the operator, if expressible in terms of special functions, or on alternative, numerically more straightforwardly accessible analytic expressions, e.g., in terms of Painleve transcendents, that have masterfully been derived in some cases. In this paper we close the gap in the literature by studying projection methods and, above all, a simple, easily implementable, general method for the numerical evaluation of Fredholm determinants that is derived from the classical Nystrom method for the solution of Fredholm equations of the second kind. Using Gauss—Legendre or Clenshaw—Curtis as the underlying quadrature rule, we prove that the approximation error essentially behaves like the quadrature error for the sections of the kernel. In particular, we get exponential convergence for analytic kernels, which are typical in random matrix theory. The application of the method to the distribution functions of the Gaussian unitary ensemble (GUE), in the bulk scaling limit and the edge scaling limit, is discussed in detail. After extending the method to systems of integral operators, we evaluate the two-point correlation functions of the more recently studied Airy and Airy 1 processes.
251 citations
TL;DR: With respect to space-time tensor-product wavelet bases, parabolic initial boundary value problems are equivalently formulated as bi-infinite matrix problems and adaptive wavelet methods are shown to yield sequences of approximate solutions which converge at the optimal rate.
Abstract: With respect to space-time tensor-product wavelet bases, parabolic initial boundary value problems are equivalently formulated as bi-infinite matrix problems. Adaptive wavelet methods are shown to yield sequences of approximate solutions which converge at the optimal rate. In case the spatial domain is of product type, the use of spatial tensor product wavelet bases is proved to overcome the so-called curse of dimensionality, i.e., the reduction of the convergence rate with increasing spatial dimension.
212 citations
TL;DR: The purpose of this paper is to prove that the linearized Bregman iteration proposed in [40] for the basis pursuit problem indeed converges.
Abstract: One of the key steps in compressed sensing is to solve the basis pursuit problem min u∈ℝ n{∥u∥: Au = f }. Bregman iteration was very successfully used to solve this problem in [40]. Also, a simple and fast iterative algorithm based on linearized Bregman iteration was proposed in [40], which is described in detail with numerical simulations in [35]. A convergence analysis of the smoothed version of this algorithm was given in [11]. The purpose of this paper is to prove that the linearized Bregman iteration proposed in [40] for the basis pursuit problem indeed converges.
192 citations
TL;DR: The new approximate flux is proven to have normal components continuous across inter-element boundaries, to converge in L2 with order k + 1, and to have a divergence converging in L1 also with orderk+1.
Abstract: We identify discontinuous Galerkin methods for second-order elliptic problems in several space dimensions having superconvergence properties similar to those of the Raviart-Thomas and the Brezzi-Douglas-Marini mixed methods. These methods use polynomials of degree k ≥ 0 for both the potential as well as the flux. We show that the approximate flux converges in L2 with the optimal order of k + 1, and that the approximate potential and its numerical trace superconverge, in L2-like norms, to suitably chosen projections of the potential, with order k + 2. We also apply element-by-element postprocessing of the approximate solution to obtain new approximations of the flux and the potential. The new approximate flux is proven to have normal components continuous across inter-element boundaries, to converge in L2 with order k + 1, and to have a divergence converging in L2 also with order k+1. The new approximate potential is proven to converge with order k+2 in L2. Numerical experiments validating these theoretical results are presented.
181 citations
TL;DR: A decoupling approach based on interface approximation via temporal extrapolation is proposed for devising decoupled marching algorithms for the mixed model that models coupled fluid flow and porous media flow.
Abstract: We study numerical methods for solving a non-stationary mixed Stokes-Darcy problem that models coupled fluid flow and porous media flow. A decoupling approach based on interface approximation via temporal extrapolation is proposed for devising decoupled marching algorithms for the mixed model. Error estimates are derived and numerical experiments are conducted to demonstrate the computational effectiveness of the decoupling approach.
179 citations
TL;DR: Formulas that allow us to decompose a function f of d variables into a sum of 2 d terms fu indexed by subsets u of f1;:::;dg, where each term fu depends only on the variables with indices in u, are presented.
Abstract: We present formulas that allow us to decompose a function f of d variables into a sum of 2 d terms fu indexed by subsets u of f1;:::;dg, where each term fu depends only on the variables with indices in u. The decomposition depends on the choice of d commuting projections fPjg d=1 , where Pj(f) does not depend on the variable xj. We present an explicit formula for fu, which is new even for the anova and anchored decompositions; both are special cases of the general decomposition. We show that the decomposition is minimal in the following sense: if f is expressible as a sum in which there is no term that depends on all of the variables indexed by the subset z then, for every choice of fPjg d=1, the terms fu = 0 for all subsets u containing z. Furthermore, in a reproducing kernel Hilbert space setting, we give su-cient conditions for the terms fu to be mutually orthogonal.
169 citations
TL;DR: In this article, a method for treating general boundary conditions in the finite element method generalizing an approach, due to Nitsche (1971), for approximating Dirichlet boundary conditions was introduced.
Abstract: We introduce a method for treating general boundary conditions in the finite element method generalizing an approach, due to Nitsche (1971), for approximating Dirichlet boundary conditions. We use Poisson's equations as a model problem and prove a priori and a posteriori error estimates. The method is also compared with the traditional Galerkin method. The theoretical results are verified numerically.
156 citations
TL;DR: A homographic best approximation problem is presented and it is shown that the new waveform relaxation algorithms are well posed and converge much faster than the classical one: the number of iterations to reach a certain accuracy can be orders of magnitudes smaller.
Abstract: We present and study a homographic best approximation problem, which arises in the analysis of waveform relaxation algorithms with optimized transmission conditions. Its solution characterizes in each class of transmission conditions the one with the best performance of the associated waveform relaxation algorithm. We present the particular class of first order transmission conditions in detail and show that the new waveform relaxation algorithms are well posed and converge much faster than the classical one: the number of iterations to reach a certain accuracy can be orders of magnitudes smaller. We illustrate our analysis with numerical experiments.
137 citations
TL;DR: A convergence analysis for exponential splitting methods applied to linear evolution equations states that the classical order of the splitting method is retained in a setting of unbounded operators, without requiring any additional order condition.
Abstract: We present a convergence analysis for exponential splitting methods applied to linear evolution equations. Our main result states that the classical order of the splitting method is retained in a setting of unbounded operators, without requiring any additional order condition. This is achieved by basing the analysis on the
abstract framework of (semi)groups. The convergence analysis also includes generalizations to splittings consisting of more then two operators, and to variable time steps. We conclude by illustrating that the abstract results are applicable in the context of the Schrodinger equation with an external magnetic field or with an
unbounded potential.
114 citations
TL;DR: This paper constructs operator-adapted subspaces with a dimension smaller than that of the standard full grid spaces but which have the same approximation order as the standardFull grid spaces, provided that certain additional regularity assumptions on the solution are fulfilled.
Abstract: This paper is concerned with the construction of optimized sparse grid approximation spaces for elliptic pseudodifferential operators of arbitrary order. Based on the framework of tensor-product biorthogonal wavelet bases and stable subspace splittings, we construct operator-adapted subspaces with a dimension smaller than that of the standard full grid spaces but which have the same approximation order as the standard full grid spaces, provided that certain additional regularity assumptions on the solution are fulfilled. Specifically for operators of positive order, their dimension is O(2 J ) independent of the dimension n of the problem, compared to O(2 Jn ) for the full grid space. Also, for operators of negative order the overall cost is significantly in favor of the new approximation spaces. We give cost estimates for the case of continuous linear information. We show these results in a constructive manner by proposing a Galerkin method together with optimal preconditioning. The theory covers elliptic boundary value problems as well as boundary integral equations.
TL;DR: An algorithm is described that determines a set of unramified covers of a given hyperelliptic curve, with the property that any rational point will lift to one of the covers, if the algorithm returns an empty set.
Abstract: We describe an algorithm that determines a set of unramified covers of a given hyperelliptic curve, with the property that any rational point will lift to one of the covers. In particular, if the algorithm returns an empty set, then the hyperelliptic curve has no rational points. This provides a relatively efficiently tested criterion for solvability of hyperelliptic curves. We also discuss applications of this algorithm to curves of genus 1 and to curves with rational points.
TL;DR: In this article, the convergence and optimality of adaptive mixed finite element methods for the Poisson equation were established and a quasi-orthogonality property was proved using the fact that the error is orthogonal to the divergence free subspace, while the part of the error that is not divergence free can be bounded by the data oscillation.
Abstract: The convergence and optimality of adaptive mixed finite element methods for the Poisson equation are established in this paper. The main difficulty for mixed finite element methods is the lack of minimization principle and thus the failure of orthogonality. A quasi-orthogonality property is proved using the fact that the error is orthogonal to the divergence free subspace, while the part of the error that is not divergence free can be bounded by the data oscillation using a discrete stability result. This discrete stability result is also used to get a localized discrete upper bound which is crucial for the proof of the optimality of the adaptive approximation.
TL;DR: In this paper, the authors used high order Runge-Kutta (RK) integrators to construct a family of related methods, which they refer to as integral deferred correction (IDC) methods.
Abstract: Spectral deferred correction (SDC) methods for solving ordinary differential equations (ODEs) were introduced by Dutt, Greengard and Rokhlin (2000). It was shown in that paper that SDC methods can achieve arbitrary high order accuracy and possess nice stability properties. Their SDC methods are constructed with low order integrators, such as forward Euler or backward Euler, and are able to handle stiff and non-stiff terms in the ODEs. In this paper, we use high order Runge-Kutta (RK) integrators to construct a family of related methods, which we refer to as integral deferred correction (IDC) methods. The distribution of quadrature nodes is assumed to be uniform, and the corresponding local error analysis is given. The smoothness of the error vector associated with an IDC method, measured by the discrete Sobolev norm, is a crucial tool in our analysis. The expected order of accuracy is demonstrated through several numerical examples. Superior numerical stability and accuracy regions are observed when high order RK integrators are used to construct IDC methods.
TL;DR: The goal of this article is to design robust and simple first order explicit solvers for one-dimensional nonconservative hyperbolic systems and to split the Roe matrices into two parts which are used to calculate the contributions at the cells to the right and to the left, respectively.
Abstract: The goal of this article is to design robust and simple first order explicit solvers for one-dimensional nonconservative hyperbolic systems. These solvers are intended to be used as the basis for higher order methods for one or multidimensional problems. The starting point for the development of these solvers is the general definition of a Roe linearization introduced by Toumi in 1992 based on the use of a family of paths. Using this concept, Roe methods can be extended to nonconservative systems. These methods have good well-balanced and robustness properties, but they have also some drawbacks: in particular, their implementation requires the explicit knowledge of the eigenstructure of the intermediate matrices. Our goal here is to design numerical methods based on a Roe linearization which overcome this drawback. The idea is to split the Roe matrices into two parts which are used to calculate the contributions at the cells to the right and to the left, respectively. This strategy is used to generate two different one-parameter families of schemes which contain, as particular cases, some generalizations to nonconservative systems of the well-known Lax-Friedrichs, Lax-Wendroff, FORCE, and GFORCE schemes. Some numerical experiments are presented to compare the behaviors of the schemes introduced here with Roe methods.
TL;DR: This paper characterize resonances in terms of (improper) eigen-functions of the Helmholtz operator on an unbounded domain and proves that the first of these steps leads to eigenvalue convergence (to the desired resonance values) which is free from spurious computational eigenvalues provided that the size of computational domain is sufficiently large.
Abstract: In this paper, we consider the problem of computing resonances in open systems. We first characterize resonances in terms of (improper) eigen-functions of the Helmholtz operator on an unbounded domain. The perfectly matched layer (PML) technique has been successfully applied to the computation of scattering problems. We shall see that the application of PML converts the resonance problem to a standard eigenvalue problem (still on an infinite domain). This new eigenvalue problem involves an operator which resembles the original Helmholtz equation transformed by a complex shift in the coordinate system. Our goal will be to approximate the shifted operator first by replacing the infinite domain by a finite (computational) domain with a convenient boundary condition and second by applying finite elements on the computational domain. We shall prove that the first of these steps leads to eigenvalue convergence (to the desired resonance values) which is free from spurious computational eigenvalues provided that the size of computational domain is sufficiently large. The analysis of the second step is classical. Finally, we illustrate the behavior of the method applied to numerical experiments in one and two spatial dimensions.
TL;DR: This paper explicitly determine hypergeometric functions modulo higher powers of p and discusses an application to supercongruences and two non-trivial generalized Harmonic sum identities discovered using the computer summation package Sigma.
Abstract: Let p be an odd prime. In 1984, Greene introduced the notion of hypergeometric functions over finite fields. Special values of these functions have been of interest as they are related to the number of Fp points on algebraic varieties and to Fourier coefficients of modular forms. In this paper, we explicitly determine these functions modulo higher powers of p and discuss an application to supercongruences. This application uses two non-trivial generalized Harmonic sum identities discovered using the computer summation package Sigma. We illustrate the usage of Sigma in the discovery and proof of these two identities.
TL;DR: This paper proposes a discretization for the (nonlinearized) compressible Stokes problem with a linear equation of state $\rho=p$, based on Crouzeix-Raviart elements, and proves estimates for the discrete solution, which yields its existence by a topological degree argument, and the convergence of the scheme to a solution of the continuous problem.
Abstract: In this paper, we propose a discretization for the (nonlinearized) compressible Stokes problem with a linear equation of state $\rho=p$, based on Crouzeix-Raviart elements The approximation of the momentum balance is obtained by usual finite element techniques Since the pressure is piecewise constant, the discrete mass balance takes the form of a finite volume scheme, in which we introduce an upwinding of the density, together with two additional stabilization terms We prove {\em a priori} estimates for the discrete solution, which yields its existence by a topological degree argument, and then the convergence of the scheme to a solution of the continuous problem
TL;DR: This work considers an initial value problem for a class of evolution equations incorporating a memory term with a weakly singular kernel bounded by C(t ― s) α-1, and proves a higher convergence rate at the nodes of the time mesh.
Abstract: We consider an initial value problem for a class of evolution equations incorporating a memory term with a weakly singular kernel bounded by C(t ― s) α-1 , where 0 < α < 1. For the time discretization we apply the discontinuous Galerkin method using piecewise polynomials of degree at most q - 1, for q = 1 or 2. For the space discretization we use continuous piecewise-linear finite elements. The discrete solution satisfies an error bound of order k 9 + h 2 l(k), where k and h are the mesh sizes in time and space, respectively, and l(k) = max(1,log k -1 ). In the case q = 2, we prove a higher convergence rate of order k 3 +h 2 l(k) at the nodes of the time mesh. Typically, the partial derivatives of the exact solution are singular at t = 0, necessitating the use of non-uniform time steps. We compare our theoretical error bounds with the results of numerical computations.
TL;DR: F Fourier expansions for the Apostol-Bernoulli and Apostoli-Euler polynomials are investigated using the Lipschitz summation formula and their integral representations are derived.
Abstract: We investigate Fourier expansions for the Apostol-Bernoulli and Apostol-Euler polynomials using the Lipschitz summation formula and obtain their integral representations We give some explicit formulas at rational arguments for these polynomials in terms of the Hurwitz zeta function We also derive the integral representations for the classical Bernoulli and Euler polynomials and related known results
TL;DR: The Bernstein inequality estimates L p Bessel-potential Sobolev norms of functions in this space in terms of the minimal separation and the L p norm of the function itself, and inverse theorems for SBF approximation in the L P norm are derived.
Abstract: The purpose of this paper is to establish L p error estimates, a Bernstein inequality, and inverse theorems for approximation by a space comprising spherical basis functions located at scattered sites on the unit n-sphere. In particular, the Bernstein inequality estimates L p Bessel-potential Sobolev norms of functions in this space in terms of the minimal separation and the L p norm of the function itself. An important step in its proof involves measuring the L p stability of functions in the approximating space in terms of the l p norm of the coefficients involved. As an application of the Bernstein inequality, we derive inverse theorems for SBF approximation in the L P norm. Finally, we give a new characterization of Besov spaces on the n-sphere in terms of spaces of SBFs.
TL;DR: It is shown that an algorithm relying on floating point evaluation of modular functions and on interpolation, which has received little attention in the literature, has a complexity that is essentially linear in the size of the computed polynomials.
Abstract: We analyse and compare the complexity of several algorithms for computing modular polynomials. We show that an algorithm relying on floating point evaluation of modular functions and on interpolation, which has received little attention in the literature, has a complexity that is essentially (up to logarithmic factors) linear in the size of the computed polynomials. In particular, it obtains the classical modular polynomials $\Phi_\ell$ of prime level $\ell$ in time O (\ell^3 \log^4 \ell \log \log \ell). Besides treating modular polynomials for $\Gamma^0 (\ell)$, which are an important ingredient in many algorithms dealing with isogenies of elliptic curves, the algorithm is easily adapted to more general situations. Composite levels are handled just as easily as prime levels, as well as polynomials between a modular function and its transform of prime level, such as the Schlafli polynomials and their generalisations. Our distributed implementation of the algorithm confirms the theoretical analysis by computing modular equations of record level around $10000$ in less than two weeks on ten processors.
TL;DR: This paper identifies function classes suited to the DE-Sinc formulas in a way compatible with the existing theoretical results for the SE-S inc formulas, as well as identifies alternative function classes for the DE/SE/DE formulas, which are more useful in applications.
Abstract: The DE-Sinc formulas, resulting from a combination of the Sinc approximation formula with the double exponential (DE) transformation, provide a highly efficient method for function approximation. In many cases they are more efficient than the SE-Sinc formulas, which are the Sinc approximation formulas combined with the single exponential (SE) transformations. Function classes suited to the SE-Sinc formulas have already been investigated in the literature through rigorous mathematical analysis, whereas this is not the case with the DE-Sinc formulas. This paper identifies function classes suited to the DE-Sinc formulas in a way compatible with the existing theoretical results for the SE-Sinc formulas. Furthermore, we identify alternative function classes for the DE-Sinc formulas, as well as for the SE-Sinc formulas, which are more useful in applications in the sense that the conditions imposed on the functions are easier to verify.
TL;DR: This work shows by direct computation that the Galois group of the Schubert problem of 3-planes in ℂ 8 meeting 15 fixed 5-planes non-trivially is the full symmetric group S 6006.
Abstract: Numerical homotopy continuation of solutions to polynomial equations is the foundation for numerical algebraic geometry, whose development has been driven by applications of mathematics. We use numerical homotopy continuation to investigate the problem in pure mathematics of determining Galois groups in the Schubert calculus. For example, we show by direct computation that the Galois group of the Schubert problem of 3-planes in ℂ 8 meeting 15 fixed 5-planes non-trivially is the full symmetric group S 6006 .
TL;DR: In this paper, it was shown that computing the number of vertices of the Voronoi cell of a lattice is #P-hard, and an algorithm for low dimensional and highly-symmetric lattices was proposed.
Abstract: In this paper we are concerned with finding the vertices of the Voronoi cell of a Euclidean lattice. Given a basis of a lattice, we prove that computing the number of vertices is a #P-hard problem. On the other hand we describe an algorithm for this problem which is especially suited for low dimensional (say dimensions at most 12) and for highly-symmetric lattices. We use our implementation, which drastically outperforms those of current computer algebra systems, to find the vertices of Voronoi cells and quantizer constants of some prominent lattices.
TL;DR: An efficient Filon-type method for the integration of systems containing Bessel functions with exotic oscillators based on a diffeomorphism transformation is presented and applications to Airy transforms are given.
Abstract: In this paper, we present an efficient Filon-type method for the integration of systems containing Bessel functions with exotic oscillators based on a diffeomorphism transformation and give applications to Airy transforms. Preliminary numerical results show the effectiveness and accuracy of the quadrature for large arguments of integral systems.
TL;DR: It is shown that linear combinations of fundamental solutions of operators of order m > 4, with singularities lying on a prescribed pseudo-boundary, are not in general dense in the corresponding solution space.
Abstract: In the present work, we investigate the applicability of the method of fundamental solutions for the solution of boundary value problems of elliptic partial differential equations and elliptic systems. More specifically, we study whether linear combinations of fundamental solutions can approximate the solutions of the boundary value problems under consideration. In our study, the singularities of the fundamental solutions lie on a prescribed pseudo―boundary ― the boundary of a domain which embraces the domain of the problem under consideration. We extend previous density results of Kupradze and Aleksidze, and of Bogomolny, to more general domains and partial differential operators, and with respect to more appropriate norms. Our domains may possess holes and their boundaries are only required to satisfy a rather weak boundary requirement, namely the segment condition. Our density results are with respect to the norms of the spaces C l (Ω). Analogous density results are obtainable with respect to Holder norms. We have studied approximation by fundamental solutions of the Laplacian, biharmonic and m―harmonic and modified Helmholtz and poly-Helmholtz operators. In the case of elliptic systems, we obtain analogous density results for the Cauchy―Navier operator as well as for an operator which arises in the linear theory of thermo―elasticity. We also study alternative formulations of the method of fundamental solutions in cases when linear combinations of fundamental solutions of the equations under consideration are not dense in the solution space. Finally, we show that linear combinations of fundamental solutions of operators of order m > 4, with singularities lying on a prescribed pseudo-boundary, are not in general dense in the corresponding solution space.
TL;DR: This work formulate the problem of finding optimal monotonicity preserving general linear methods for linear autonomous equations, and proposes an efficient algorithm that reliably finds optimal methods even among classes involving very high order accuracy and that use many steps and/or stages.
Abstract: Monotonicity preserving numerical methods for ordinary differential equations prevent the growth of propagated errors and preserve convex boundedness properties of the solution. We formulate the problem of finding optimal monotonicity preserving general linear methods for linear autonomous equations, and propose an efficient algorithm for its solution. This algorithm reliably finds optimal methods even among classes involving very high order accuracy and that use many steps and/or stages. The optimality of some recently proposed methods is verified, and many more efficient methods are found. We use similar algorithms to find optimal strong stability preserving linear multistep methods of both explicit and implicit type, including methods for hyperbolic PDEs that use downwind-biased operators.