Showing papers in "arXiv: Numerical Analysis in 2006"

Tensor rank and the ill-posedness of the best low-rank approximation problem

Abstract: There has been continued interest in seeking a theorem describing optimal low-rank approximations to tensors of order 3 or higher, that parallels the Eckart-Young theorem for matrices. In this paper, we argue that the naive approach to this problem is doomed to failure because, unlike matrices, tensors of order 3 or higher can fail to have best rank-r approximations. The phenomenon is much more widespread than one might suspect: examples of this failure can be constructed over a wide range of dimensions, orders and ranks, regardless of the choice of norm (or even Bregman divergence). Moreover, we show that in many instances these counterexamples have positive volume: they cannot be regarded as isolated phenomena. In one extreme case, we exhibit a tensor space in which no rank-3 tensor has an optimal rank-2 approximation. The notable exceptions to this misbehavior are rank-1 tensors and order-2 tensors. In a more positive spirit, we propose a natural way of overcoming the ill-posedness of the low-rank approximation problem, by using weak solutions when true solutions do not exist. In our work we emphasize the importance of closely studying concrete low-dimensional examples as a first step towards more general results. To this end, we present a detailed analysis of equivalence classes of 2-by-2-by-2 tensors, and we develop methods for extending results upwards to higher orders and dimensions. Finally, we link our work to existing studies of tensors from an algebraic geometric point of view. The rank of a tensor can in theory be given a semialgebraic description; i.e., can be determined by a system of polynomial inequalities. In particular we make extensive use of the 2-by-2-by-2 hyperdeterminant.

Recovery algorithms for vector valued data with joint sparsity constraints

Abstract: Vector valued data appearing in concrete applications often possess sparse expansions with respect to a preassigned frame for each vector component individually. Additionally, different components may also exhibit common sparsity patterns. Recently, there were introduced sparsity measures that take into account such joint sparsity patterns, promoting coupling of non-vanishing components. These measures are typically constructed as weighted $\ell_1$ norms of componentwise $\ell_q$ norms of frame coefficients. We show how to compute solutions of linear inverse problems with such joint sparsity regularization constraints by fast thresholded Landweber algorithms. Next we discuss the adaptive choice of suitable weights appearing in the definition of sparsity measures. The weights are interpreted as indicators of the sparsity pattern and are iteratively up-dated after each new application of the thresholded Landweber algorithm. The resulting two-step algorithm is interpreted as a double-minimization scheme for a suitable target functional. We show its $\ell_2$-norm convergence. An implementable version of the algorithm is also formulated, and its norm convergence is proven. Numerical experiments in color image restoration are presented.

Fast linear algebra is stable

Abstract: In an earlier paper, we showed that a large class of fast recursive matrix multiplication algorithms is stable in a normwise sense, and that in fact if multiplication of $n$-by-$n$ matrices can be done by any algorithm in $O(n^{\omega + \eta})$ operations for any $\eta > 0$, then it can be done stably in $O(n^{\omega + \eta})$ operations for any $\eta > 0$. Here we extend this result to show that essentially all standard linear algebra operations, including LU decomposition, QR decomposition, linear equation solving, matrix inversion, solving least squares problems, (generalized) eigenvalue problems and the singular value decomposition can also be done stably (in a normwise sense) in $O(n^{\omega + \eta})$ operations.

Fast Computation of Fourier Integral Operators

Abstract: We introduce a general purpose algorithm for rapidly computing certain types of oscillatory integrals which frequently arise in problems connected to wave propagation and general hyperbolic equations. The problem is to evaluate numerically a so-called Fourier integral operator (FIO) of the form $\int e^{2\pi i \Phi(x,\xi)} a(x,\xi) \hat{f}(\xi) \mathrm{d}\xi$ at points given on a Cartesian grid. Here, $\xi$ is a frequency variable, $\hat f(\xi)$ is the Fourier transform of the input $f$, $a(x,\xi)$ is an amplitude and $\Phi(x,\xi)$ is a phase function, which is typically as large as $|\xi|$; hence the integral is highly oscillatory at high frequencies. Because an FIO is a dense matrix, a naive matrix vector product with an input given on a Cartesian grid of size $N$ by $N$ would require $O(N^4)$ operations. This paper develops a new numerical algorithm which requires $O(N^{2.5} \log N)$ operations, and as low as $O(\sqrt{N})$ in storage space. It operates by localizing the integral over polar wedges with small angular aperture in the frequency plane. On each wedge, the algorithm factorizes the kernel $e^{2 \pi i \Phi(x,\xi)} a(x,\xi)$ into two components: 1) a diffeomorphism which is handled by means of a nonuniform FFT and 2) a residual factor which is handled by numerical separation of the spatial and frequency variables. The key to the complexity and accuracy estimates is that the separation rank of the residual kernel is \emph{provably independent of the problem size}. Several numerical examples demonstrate the efficiency and accuracy of the proposed methodology. We also discuss the potential of our ideas for various applications such as reflection seismology.

Analysis of a Force-Based Quasicontinuum Approximation

Abstract: We analyze a force-based quasicontinuum approximation to a one-dimensional system of atoms that interact by a classical atomistic potential. This force-based quasicontinuum approximation is derived as the modification of an energy-based quasicontinuum approximation by the addition of nonconservative forces to correct nonphysical ``ghost'' forces that occur in the atomistic to continuum interface. We prove that the force-based quasicontinuum equations have a unique solution under suitable restrictions on the loads. For Lennard-Jones next-nearest-neighbor interactions, we show that unique solutions exist for loads in a symmetric region extending nearly to the tensile limit. We give an analysis of the convergence of the ghost force iteration method to solve the equilibrium equations for the force-based quasicontinuum approximation. We show that the ghost force iteration is a contraction and give an analysis for its convergence rate.

Adaptive, Fast and Oblivious Convolution in Evolution Equations with Memory

Abstract: To approximate convolutions which occur in evolution equations with memory terms, a variable-stepsize algorithm is presented for which advancing N steps requires only O(N log(N)) operations and O(log(N)) active memory, in place of O(N^2) operations and O(N) memory for a direct implementation. A basic feature of the fast algorithm is the reduction, via contour integral representations, to differential equations which are solved numerically with adaptive step sizes. Rather than the kernel itself, its Laplace transform is used in the algorithm. The algorithm is illustrated on three examples: a blow-up example originating from a Schr\"odinger equation with concentrated nonlinearity, chemical reactions with inhibited diffusion, and viscoelasticity with a fractional order constitutive law.

Sparse reconstruction by convex relaxation: Fourier and Gaussian measurements

Abstract: We want to exactly reconstruct a sparse signal f (a vector in R^n of small support) from few linear measurements of f (inner products with some fixed vectors). A nice and intuitive reconstruction by Linear Programming has been advocated since 80-ies by Dave Donoho and his collaborators. Namely, one can relax the reconstruction problem, which is highly nonconvex, to a convex problem -- and, moreover, to a linear program. However, when is exactly the reconstruction problem equivalent to its convex relaxation is an open question. Recent work of many authors shows that the number of measurements k(r,n) needed to exactly reconstruct any r-sparse signal f of length n (a vector in R^n of support r) from its linear measurements with the convex relaxation method is usually O(r polylog(n)). However, known estimates of the number of measurements k(r,n) involve huge constants, in spite of very good performance of the algorithms in practice. In this paper, we consider random Gaussian measurements and random Fourier measurements (a frequency sample of f). For Gaussian measurements, we prove the first guarantees with reasonable constants: k(r,n) < 12 r (2 + log(n/r)), which is optimal up to constants. For Fourier measurements, we prove the best known bound k(r,n) = O(r log(n) . log^2(r) log(r log n)), which is optimal within the log log n and log^3 r factors. Our arguments are based on the technique of Geometric Functional Analysis and Probability in Banach spaces.

High order relaxation schemes for non linear degenerate diffusion problems

Abstract: Several relaxation approximations to partial differential equations have been recently proposed. Examples include conservation laws, Hamilton-Jacobi equations, convection-diffusion problems, gas dynamics problems. The present paper focuses onto diffusive relaxation schemes for the numerical approximation of nonlinear parabolic equations. These schemes are based on suitable semilinear hyperbolic system with relaxation terms. High order methods are obtained by coupling ENO and WENO schemes for space discretization with IMEX schemes for time integration. Error estimates and convergence analysis are developed for semidiscrete schemes with numerical analysis for fully discrete relaxed schemes. Various numerical results in one and two dimension illustrate the high accuracy and good properties of the proposed numerical schemes. These schemes can be easily implemented for parallel computer and applied to more general system of nonlinear parabolic equations in two- and three-dimensional cases.

Adaptive Weak Approximation of Diffusions with Jumps

Abstract: This work develops Monte Carlo Euler adaptive time stepping methods for the weak approximation problem of jump diffusion driven stochastic differential equations. The main result is the derivation of a new expansion for the omputational error, with computable leading order term in a posteriori form, based on stochastic flows and discrete dual backward problems which extends the results in [STZ]. These expansions lead to efficient and accurate computation of error estimates. Adaptive algorithms for either stochastic time steps or quasi-deterministic time steps are described. Numerical examples show the performance of the proposed error approximation and of the described adaptive time-stepping methods.

A family of C^0 finite elements for Kirchhoff plates I: Error analysis

Abstract: A new finite element formulation for the Kirchhoff plate model is presented. The method is a displacement formulation with the deflection and the rotation vector as unknowns and it is based on ideas stemming from a stabilized method for the Reissner--Mindlin model and a method to treat a free boundary. Optimal a-priori and a-posteriori error estimates are derived.

Smoothed analysis of complex conic condition numbers

Abstract: Smoothed analysis of complexity bounds and condition numbers has been done, so far, on a case by case basis. In this paper we consider a reasonably large class of condition numbers for problems over the complex numbers and we obtain smoothed analysis estimates for elements in this class depending only on geometric invariants of the corresponding sets of ill-posed inputs. These estimates are for a version of smoothed analysis proposed in this paper which, to the best of our knowledge, appears to be new. Several applications to linear and polynomial equation solving show that estimates obtained in this way are easy to derive and quite accurate.

Numerical solutions to boundary value problem for anomalous diffusion equation with Riesz-Feller fractional operator

Abstract: In this paper, we present a numerical solution to an ordinary differential equation of a fractional order in one-dimensional space. The solution to this equation can describe a steady state of the process of anomalous diffusion. The process arises from interactions within complex and non-homogeneous background. We present a numerical method which is based on the finite differences method. We consider a boundary value problem (Dirichlet conditions) for an equation with the Riesz-Feller fractional derivative. In the final part of this paper, same simulation results are shown. We present an example of non-linear temperature profiles in nanotubes which can be approximated by a solution to the fractional differential equation.

Coarse-graining schemes and a posteriori error estimates for stochastic lattice systems

Abstract: The primary objective of this work is to develop coarse-graining schemes for stochastic many-body microscopic models and quantify their effectiveness in terms of a priori and a posteriori error analysis. In this paper we focus on stochastic lattice systems of interacting particles at equilibrium. %such as Ising-type models. The proposed algorithms are derived from an initial coarse-grained approximation that is directly computable by Monte Carlo simulations, and the corresponding numerical error is calculated using the specific relative entropy between the exact and approximate coarse-grained equilibrium measures. Subsequently we carry out a cluster expansion around this first-and often inadequate-approximation and obtain more accurate coarse-graining schemes. The cluster expansions yield also sharp a posteriori error estimates for the coarse-grained approximations that can be used for the construction of adaptive coarse-graining methods. We present a number of numerical examples that demonstrate that the coarse-graining schemes developed here allow for accurate predictions of critical behavior and hysteresis in systems with intermediate and long-range interactions. We also present examples where they substantially improve predictions of earlier coarse-graining schemes for short-range interactions.

Solution of the time-harmonic Maxwell equations using discontinuous Galerkin methods

Abstract: We present numerical results concerning the solution of the time-harmonic Maxwell's equations discretized by discontinuous Galerkin methods. In particular, a numerical study of the convergence, which compares different strategies proposed in the literature for the elliptic Maxwell equations, is performed in the two-dimensional case.

The algebra of the box-spline

Abstract: In this paper we want to revisit results of Dahmen and Micchelli on box-splines which we reinterpret and make more precise. We compare these ideas with the work of Brion, Szenes, Vergne and others on polytopes and partition functions.

Optimal Approximation of Elliptic Problems by Linear and Nonlinear Mappings II

Abstract: We study the optimal approximation of the solution of an operator equation by certain n-term approximations with respect to specific classes of frames. We study worst case errors and the optimal order of convergence and define suitable nonlinear frame widths. The main advantage of frames compared to Riesz basis, which were studied in our earlier papers, is the fact that we can now handle arbitrary bounded Lipschitz domains--also for the upper bounds. Key words: elliptic operator equation, worst case error, frames, nonlinear approximation, best n-term approximation, manifold width, Besov spaces on Lipschitz domains

Optimized Schwarz Methods for Maxwell equations

Abstract: Over the last two decades, classical Schwarz methods have been extended to systems of hyperbolic partial differential equations, and it was observed that the classical Schwarz method can be convergent even without overlap in certain cases. This is in strong contrast to the behavior of classical Schwarz methods applied to elliptic problems, for which overlap is essential for convergence. Over the last decade, optimized Schwarz methods have been developed for elliptic partial differential equations. These methods use more effective transmission conditions between subdomains, and are also convergent without overlap for elliptic problems. We show here why the classical Schwarz method applied to the hyperbolic problem converges without overlap for Maxwell's equations. The reason is that the method is equivalent to a simple optimized Schwarz method for an equivalent elliptic problem. Using this link, we show how to develop more efficient Schwarz methods than the classical ones for the Maxwell's equations. We illustrate our findings with numerical results.

The Spectral Basis and Rational Interpolation

Abstract: The Euclidean Algorithm is the often forgotten key to rational approximation techniques, including Taylor, Lagrange, Hermite, osculating, cubic spline, Chebyshev, Pade and other interpolation schemes A unified view of these various interpolation techniques is eloquently expressed in terms of the concept of the spectral basis of a factor ring of polynomials When these methods are applied to the minimal polynomial of a matrix, they give a family of rational forms of functions of that matrix

Spectral algorithms for reaction-diffusion equations

Abstract: A collection of codes (in MATLAB & Fortran 77), and examples, for solving reaction-diffusion equations in one and two space dimensions is presented. In areas of the mathematical community spectral methods are used to remove the stiffness associated with the diffusive terms in a reaction-diffusion model allowing explicit high order timestepping to be used. This is particularly valuable for two (and higher) space dimension problems. Our aim here is to provide codes, together with examples, to allow practioners to easily utilize, understand and implement these ideas; we incorporate recent theoretical advances such as exponential time differencing methods and provide timings and error comparisons with other more standard approaches. The examples are chosen from the literature to illustrate points and queries that naturally arise.

On Wu and Schaback's Error Bound

Abstract: In this paper we point out that the commonly used error bound in the theory of radial basis functions contains an important error It doesn't apply for derivatives Thus one should be very careful when dealing with differential equations

Iterative solution of linear equations with unbounded operators

Abstract: A convergent iterative process is constructed for solving any solvable linear equation in a Hilbert space.

Discrete momentum maps for lattice EPDiff

Abstract: We focus on the spatial discretization produced by the Variational Particle-Mesh (VPM) method for a prototype fluid equation the known as the EPDiff equation}, which is short for Euler-Poincar\'e equation associated with the diffeomorphism group (of $\mathbb{R}^d$, or of a $d$-dimensional manifold $\Omega$). The EPDiff equation admits measure valued solutions, whose dynamics are determined by the momentum maps for the left and right actions of the diffeomorphisms on embedded subspaces of $\mathbb{R}^d$. The discrete VPM analogs of those dynamics are studied here. Our main results are: (i) a variational formulation for the VPM method, expressed in terms of a constrained variational principle principle for the Lagrangian particles, whose velocities are restricted to a distribution $D_{\VPM}$ which is a finite-dimensional subspace of the Lie algebra of vector fields on $\Omega$; (ii) a corresponding constrained variational principle on the fixed Eulerian grid which gives a discrete version of the Euler-Poincar\'e equation; and (iii) discrete versions of the momentum maps for the left and right actions of diffeomorphisms on the space of solutions.

Analysing singularities of a benchmark problem

Abstract: The purpose of this paper is to analyze the singularities of a well known benchmark problem ``Andrews' squeezing mechanism''. We show that for physically relevant parameter values this system admits singularities. The method is based on Gr\"obner bases computations and ideal decomposition. It is algorithmic and can thus be applied to study constraint singularities which arise in more general situations.

Spectrum of the Jacobi tau approximation for the second derivative operator

Abstract: It is proved that the eigenvalues of the Jacobi Tau method for the second derivative operator with Dirichlet boundary conditions are real, negative and distinct for a range of the Jacobi parameters. Special emphasis is placed on the symmetric case of the Gegenbauer Tau method where the range of parameters included in the theorems can be extended and characteristic polynomials given by successive order approximations interlace. This includes the common Chebyshev and Legendre, Tau and Galerkin methods. The characteristic polynomials for the Gegenbauer Tau method are shown to obey three term recurrences plus a constant term which vanishes for the Legendre Tau and Galerkin cases. These recurrences are equivalent to a tridiagonal plus one row matrix structure. The spectral integration formulation of the Gegenbauer Tau method is shown to lead directly to that fundamental and well-conditioned tridiagonal plus one row matrix structure. A Matlab code is provided.

Stability estimates for h-p spectral element methods for general elliptic problems on curvilinear domains

Abstract: In this paper we show that the h-p spectral element method developed in \cite{pdstrk1,tomar-01,tomar-dutt-kumar-02} applies to elliptic problems in curvilinear polygons with mixed Neumann and Dirichlet boundary conditions provided that the Babuska--Brezzi inf--sup conditions are satisfied.

Two-step hybrid methods adapted to the numerical integration of perturbed oscillators

Abstract: Two-step hybrid methods specially adapted to the numerical integration of perturbed oscillators are obtained. The formulation of the methods is based on a refinement of classical Taylor expansions due to Scheifele [{\em Z. Angew. Math. Phys.}, {\bf 22}, 186--210 (1971)]. The key property is that those algorithms are able to integrate exactly harmonic oscillators with frequency $\omega$ and that, for perturbed oscillators, the local error contains the (small) perturbation parameter as a factor. The methods depend on a parameter $ u=\omega h$, where $h$ is the stepsize. Based on the B2-series theory of Coleman [{\em IMA J. Numer. Anal.}, {\bf 23}, 197--220 (2003)] we derive the order conditions of this new type of methods. The linear stability and phase properties are examined. The theory is illustrated with some fourth- and fifth-order explicit schemes. Numerical results carried out on an assortment of test problems (such as the integration of the orbital motion of earth satellites) show the relevance of the theory.

Exit from a basin of attraction for stochastic weakly damped nonlinear Schr\"odinger equations

Abstract: We consider weakly damped nonlinear Schr\"odinger equations perturbed by a noise of small amplitude. The small noise is either complex and of additive type or real and of multiplicative type. It is white in time and colored in space. Zero is an asymptotically stable equilibrium point of the deterministic equations. We study the exit from a neighborhood of zero, invariant by the flow of the deterministic equation, in $\xLtwo$ or in $\xHone$. Due to noise, large fluctuations off zero occur. Thus, on a sufficiently large time scale, exit from these domains of attraction occur. A formal characterization of the small noise asymptotic of both the first exit times and the exit points is given.

Uncertainty Principles and Vector Quantization

Abstract: Given a frame in C^n which satisfies a form of the uncertainty principle (as introduced by Candes and Tao), it is shown how to quickly convert the frame representation of every vector into a more robust Kashin's representation whose coefficients all have the smallest possible dynamic range O(1/\sqrt{n}). The information tends to spread evenly among these coefficients. As a consequence, Kashin's representations have a great power for reduction of errors in their coefficients, including coefficient losses and distortions.

Grid Generation and Adaptation by Functionals

Abstract: Accuracy of a simulation is strongly depend on the grid quality. Here, quality means orthogonality at the boundaries and quasi-orthogonality within the critical regions, smoothness, bounded aspect ratios, solution adaptive behaviour, etc. We review various functionals for generating high quality structured quadrilateral meshes in two dimensional domains. Analysis of Winslow and Modified Liao functionals are presented. Numerical examples are also presented to support our theoretical analysis. We will demonstrate the use of the Area functional for generating adaptive quadrilateral meshes.

Remarks on the Cayley Representation of Orthogonal Matrices and on Perturbing the Diagonal of a Matrix to Make it Invertible

Abstract: This note contains two remarks. The first remark concerns the extension of the well-known Cayley representation of rotation matrices by skew symmetric matrices to rotation matrices admitting -1 as an eigenvalue and then to all orthogonal matrices. We review a method due to Hermann Weyl and another method involving multiplication by a diagonal matrix whose entries are +1 or -1. The second remark has to do with ways of flipping the signs of the entries of a diagonal matrix, C, with nonzero diagonal entries, obtaining a new matrix, E, so that E + A is invertible, where A is any given matrix (invertible or not).