Search

This paper introduces a new sweeping preconditioner for the iterative solution of the variable coefficient Helmholtz equation in two and three dimensions. The algorithms follow the general structure of constructing an approximate LDLI factorization by eliminating the unknowns layer by layer starting from an absorbing layer or boundary condition. The central idea of this paper is to approximate the Schur complement matrices of the factorization using moving perfectly matched layers (PMLs) introduced in the interior of the domain. Applying each Schur complement matrix is equivalent to solving a quasi-1D problem with a banded LU factorization in the 2D case and to solving a quasi-2D problem with B multifrontal method in the 3D case. The resulting preconditioner has linear application cost, and the preconditioned iterative solver converges in a number of iterations that is essentially independent of the number of unknowns or the frequency. Numerical results are presented in both two and three dimensions to demonstrate the efficiency of this new preconditioner. [ABSTRACT FROM AUTHOR]

This paper concerns numerical methods for computing complex geometrical optics (CGO) solutions to the conductivity equation ∇ ⋅ σ∇u (⋅, k) = 0 in ℝ² for piecewise smooth conductivities σ, where k is a complex parameter. The key is to solve an R-linear singular integral equation defined in the unit disk. Recently, Astala et al. [Appl. Comput. Harmon. Anat., 29 (2010), pp. 2-17] proposed a complicated method for numerical computation of CGO solutions by solving a periodic version of the ℝ-linear integral equation in a rectangle containing the unit disk. In this paper, based on the fast algorithms in [P. Daripa and D. Mashat, Numer. Algorithms, 18 (1998), pp. 133-157] for singular integral transforms, we propose a simpler numerical method which solves the ℝ-linear integral equation in the unit disk directly. For the resulting ℝ-linear operator equation, a minimal residual iterative method is proposed. Numerical examples illustrate the accuracy and efficiency of the new method. [ABSTRACT FROM AUTHOR]

An earlier paper proved the convergence of a variable stepsize Bregman operator splitting algorithm (BOSVS) for minimizing ø(Bu) + H(u), where H and ø are convex functions, and ø is possibly nonsmooth. The algorithm was shown to be relatively efficient when applied to partially parallel magnetic resonance image reconstruction problems. In this paper, the convergence rate of BOSVS is analyzed. When H(u) = ‖Au - f ‖2, where A is a matrix, it is shown that for an ergodic approximation uk obtained by averaging k BOSVS iterates, the error in the objective value ø(Buk)+H(uk) is O(1/k). When the optimization problem has a unique solution u*, we obtain the estimate ‖uk-u*‖ = O(1/ √ k). The theoretical analysis is compared to observed convergence rates for partially parallel magnetic resonance image reconstruction problems where A is a large dense ill-conditioned matrix. [ABSTRACT FROM AUTHOR]

A framework, PhaseLift, was recently proposed to solve the phase retrieval problem. In this framework, the problem is solved by optimizing a cost function over the set of complex Hermitian positive semidefinite matrices. This approach to phase retrieval motivates a more general consideration of optimizing cost functions on semidefinite Hermitian matrices where the desired minimizers are known to have low rank. This paper considers an approach based on an alternative cost function defined on a union of appropriate manifolds. It is related to the original cost function in a manner that preserves the ability to find a global minimizer and is significantly more efficient computationally. A rank-based optimality condition for stationary points is given and optimization algorithms based on state-of-the-art Riemannian optimization and dynamically reducing rank are proposed. Empirical evaluations are performed using the PhaseLift problem. The new approach is shown to be an effective method of phase retrieval with computational efficiency increased substantially compared to a state-of-the-art algorithm, the Wirtinger flow algorithm. A preliminary version of this paper can be found in [W. Huang, K. A. Gallivan, and X. Zhang, Procedia Comput. Sci., 80 (2016), pp. 1125-1134]. [ABSTRACT FROM AUTHOR]

The purpose of this paper is to present a theoretic and numerical study of utilizing squeezing and phase shift in coherent feedback control of linear quantum optical systems. A quadrature representation with built-in phase shifters is proposed for such systems. Fundamental structural characterizations of linear quantum optical systems are derived in terms of the new quadrature representation. These results reveal considerable insights into the issue of the physical realizability of such quantum systems. The problem of coherent quantum linear quadratic Gaussian (LQG) feedback control studied in H. I. Nurdin, M. R. James, and I. R. Petersen, Automatica, IFAC, 45 (2009), pp. 1837-1846; G. Zhang and M. R. James, IEEE Trans. Automat. Control, 56 (2011), pp. 1535-1550 is reinvestigated in depth. First, the optimization methods in these papers are extended to a multistep optimization algorithm which utilizes ideal squeezers. Second, a two-stage optimization approach is proposed on the basis of controller parametrization. Numerical studies show that closed-loop systems designed via the second approach may offer LQG control performance even better than that when the closed-loop systems are in the vacuum state. When ideal squeezers in a closed-loop system are replaced by (more realistic) degenerate parametric amplifiers, a sufficient condition is derived for the asymptotic stability of the resultant new closed-loop system; the issue of performance convergence is also discussed in the LQG control setting. [ABSTRACT FROM AUTHOR]

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. [ABSTRACT FROM AUTHOR]

Let A ∈ ℂN x N be a nonsingular complex matrix and b and c be complex vectors of length N. The goal of this paper is to investigate approaches for efficient approximations of the bilinear form c*A-1b. Equivalently, we wish to approximate the scalar value c*x, where x solves the linear system Ax = b. Here the matrix A can be very large or its elements can be too costly to compute so that A is not explicitly available and it is used only in the form of the matrix-vector product. Therefore a direct method is not an option. For A Hermitian positive definite, b*A-1b can be efficiently approximated as a by-product of the conjugate-gradient iterations, which is mathematically equivalent to the matching moment approximations computed via the Gauss-Christoffel quadrature. In this paper we propose a new method using the biconjugate gradient iterations which is applicable to the general complex case. The proposed approach will be compared with existing ones using analytic arguments and numerical experiments. [ABSTRACT FROM AUTHOR]

A systematic mathematical framework for the study of numerical algorithms would allow comparisons, facilitate conjugacy arguments, and enable the discovery of improved, accelerated, data-driven algorithms. Over the course of the past century, the Koopman operator has provided a mathematical framework for the study of dynamical systems which facilitates conjugacy arguments and can provide efficient reduced descriptions. More recently, numerical approximations of the operator have enabled the analysis of a large number of deterministic and stochastic dynamical systems in a completely data-driven, essentially equation-free pipeline. Discrete- or continuous-time numerical algorithms (integrators, nonlinear equation solvers, optimization algorithms are themselves dynamical systems. In this paper, we use this insight to leverage the Koopman operator framework in the data-driven study of such algorithms and discuss benefits for analysis and acceleration of numerical computation. For algorithms acting on high-dimensional spaces by quickly contracting them toward low-dimensional manifolds, we demonstrate how basis functions adapted to the data help to construct efficient reduced representations of the operator. Our illustrative examples include the gradient descent and Nesterov optimization algorithms as well as the Newton-Raphson algorithm. [ABSTRACT FROM AUTHOR]

In this paper we use a displacement structure approach to design a class of new preconditioners for the conjugate gradient method applied to the solution of large Toeplitz linear equations. Explicit formulas are suggested for the G.Strang-type and for the T.Chan-type preconditioners belonging to any of eight classes of matrices diagonalized by the corresponding discrete cosine or sine transforms. Under the standard Wiener class assumption the clustering property is established for all of these preconditioners, guaranteeing rapid convergence of the preconditioned con- jugate gradient method. All the computations related to the new preconditioners can be done in real arithmetic, and to fully exploit this advantageous property one has to suggest a fast real-arithmetic algorithm for multiplication of a Toeplitz matrix by a vector. It turns out that the obtained formulas for the Strang-type preconditioners allow a number of representations for Toeplitz matrices leading to a wide variety of real-arithmetic multiplication algorithms based on any of eight discrete cosine or sine transforms. Recently, transformations of Toeplitz matrices to Vandermonde-like or Cauchy-like matrices have been found to be useful in developing accurate direct methods for Toeplitz linear equations. In this paper we suggest further extending the range of the transformation approach by exploring it for iterative methods; this technique allowed us to reduce the complexity of each iteration of the preconditioned conjugate gradient method. The results of this paper were announced in [T. Kailath and V. Olshevsky, Calcolo, 33 (1996), pp. 191-208]. [ABSTRACT FROM AUTHOR]

In (SIAM J. Matrix Anal. Appl., 23 (2002), pp. 706–727, we developed numerical algorithms for computing sparse low-rank approximations of matrices, and we also provided a detailed error analysis of the proposed algorithms together with some numerical experiments. The low-rank approximations are constructed in a certain factored form with the degree of sparsity of the factors controlled by some user-specified parameters. In this paper, we cast the sparse low-rank approximation problem in the framework of penalized optimization problems. We discuss various approximation schemes for the penalized optimization problem which are more amenable to numerical computations. We also include some analysis to show the relations between the original optimization problem and the reduced one. We then develop a globally convergent discrete Newton-like iterative method for solving the approximate penalized optimization problems. We also compare the reconstruction errors of the sparse low-rank approximations computed by our new methods with those obtained using the methods in the earlier paper and several other existing methods for computing sparse low-rank approximations. Numerical examples show that the penalized methods are more robust and produce approximations with factors which have fewer columns and are sparser. [ABSTRACT FROM AUTHOR]

This paper discusses convergence behavior of the nonlinear inexact Uzawa algorithm for solving saddle-point problems presented in a recent paper of Bramble, Pasciak, and Vassilev (SIAM J. Numer. Anal., 34 (1997), pp. 1072­1092). It is shown that this algorithm converges under a condition weaker than that stated in their paper. Key words: saddle-point problems; nonlinear inexact Uzawa algorithm; convergence behavior [ABSTRACT FROM AUTHOR]

Let F be a compact set of a Banach space X. This paper analyzes the "generalized empirical interpolation method," which, given a function f ∊ F, builds an interpolant Jn[f] in an n dimensional subspace Xn ⊂ X with the knowledge of n outputs (σi(f))ni=1, where σi ∊ X' and X' is the dual space of X. The space Xn is built with a greedy algorithm that is adapted to F in the sense that it is generated by elements of F itself. The algorithm also selects the linear functionals (σi)ni=1 from a dictionary Σ ⊂ X'. In this paper, we study the interpolation error maxf∊F ‖f - Jn[f]‖X by comparing it with the best possible performance on an n dimensional space, i.e., the Kolmogorov n-width of F in X, dn(F,X). For polynomial or exponential decay rates of dn(F,X), we prove that the interpolation error has the same behavior modulo the norm of the interpolation operator. Sharper results are obtained in the case where X is a Hilbert space. [ABSTRACT FROM AUTHOR]

This paper considers an implicitly restarted Krylov subspace method that approximates a stable, linear transfer function f(s) of order n by one of order m, where n >> m. It is well known that oblique projections onto a Krylov subspace may generate unstable partial realizations. To remedy this situation, the oblique projectors obtained via classical Krylov subspace methods are supplemented with further projectors which enable the formation of stable partial realizations directly from f(s). A key feature of this process is that it may be incorporated into an implicit restart scheme. A second difficulty arises from the fact that Krylov subspace methods often generate partial realizations that contain nonessential modes. To this end, balanced truncation may be employed to discard the unwanted part of the reduced-order model. This paper proposes oblique projection methods for large-scale model reduction that simultaneously compute stable reduced-order models while discarding all nonessential modes. It is shown that both of these tasks may be effected by a single oblique projection process. Furthermore, the process is shown to naturally fit into an implicit restart framework. The theoretical properties of these methods are thoroughly investigated, and exact low-dimensional expressions for the L∞-norm of the residual errors are derived. Finally, the behavior of the algorithm is illustrated on two large-scale examples. [ABSTRACT FROM AUTHOR]

The generalized null space decomposition (GNSD) is a unitary reduction of a general matrix A of order n to a block upper triangular form that reveals the structure of the Jordan blocks of A corresponding to a zero eigenvalue. The reduction was introduced by Kublanovskaya. It was extended first by Ruhe and then by Golub and Wilkinson, who based the reduction on the singular value decomposition. Unfortunately, if A has large Jordan blocks, the complexity of these algorithms can approach the order of n4. This paper presents an alternative algorithm, based on repeated updates of a QR decomposition of A, that is guaranteed to be of order n³. Numerical experiments confirm the stability of this algorithm, which turns out to produce essentially the same form as that of Golub and Wilkinson. The effect of errors in A on the ability to recover the original structure is investigated empirically. Several applications are discussed, including the computation of the Drazin inverse. [ABSTRACT FROM AUTHOR]

We introduce a new algorithm of proper generalized decomposition (PGD) for parametric symmetric elliptic partial differential equations. For any given dimension, we prove the existence of an optimal subspace of at most that dimension which realizes the best approximation---in the mean parametric norm associated to the elliptic operator---of the error between the exact solution and the Galerkin solution calculated on the subspace. This is analogous to the best approximation property of the proper orthogonal decomposition (POD) subspaces, except that in our case the norm is parameter-dependent. We apply a deflation technique to build a series of approximating solutions on finite-dimensional optimal subspaces, directly in the online step, and we prove that the partial sums converge to the continuous solution in the mean parametric elliptic norm. We show that the standard PGD for the considered parametric problem is strongly related to the deflation algorithm introduced in this paper. This opens the possibility of computing the PGD expansion by directly solving the optimization problems that yield the optimal subspaces. [ABSTRACT FROM AUTHOR]

The Euler genus of a graph is a fundamental and well-studied parameter in graph theory and topology. Computing it has been shown to be NP-hard by Thomassen [J. Algorithms, 10 (1989), pp. 568{576; J. Combin. Theory, Ser. B, 57 (1993), pp. 196{206], and it is known to be fixed-parameter tractable. However, the approximability of the Euler genus is wide open. While the existence of an O(1)-approximation is not ruled out, only an O( p n)-approximation [J. Chen, S. P. Kanchi, and A. Kanevsky, Inform. Process. Lett., 61 (1997), pp. 317{322] is known even in bounded-degree graphs. In this paper we give a polynomial-time algorithm which, given a boundeddegree graph of Euler genus g, computes a drawing in a surface of Euler genus gO(1) · logO(1) n. Combined with the upper bound from [J. Chen, S. P. Kanchi, and A. Kanevsky, Inform. Process. Lett., 61 (1997), pp. 317{322], our result also implies a O(n1/2-α)-approximation for some constant α > 0. Using our algorithm for approximating the Euler genus as a subroutine, we obtain, in a uniform fashion, algorithms with approximation ratios of the form OPTO(1) · logO(1) n for several related problems on bounded-degree graphs. These include the problems of orientable genus, crossing number, and planar edge and vertex deletion. Our algorithm and proof of correctness for the crossing number problem are simpler compared to the long and difficult proof in the recent breakthrough by Chuzhoy [Proceedings of the ACM Symposium on Theory of Computing, 2011, pp. 303{312], while essentially obtaining a qualitatively similar result. For planar edge and vertex deletion problems our results are the first to obtain a bound of the form poly(OPT; log n). We also highlight some further applications of our results in the design of algorithms for graphs with small genus. Many such algorithms require that a drawing of the graph is given as part of the input. Our results imply that in several interesting cases, we can implement such algorithms even when the drawing is unknown. [ABSTRACT FROM AUTHOR]

Threshold dynamics and its extensions have proven useful in computing interfacial motions in applications as diverse as materials science and machine learning. Certain desirable properties of the algorithm, such as unconditional monotonicity in two-phase ows and gradient stability more generally, hinge on positivity properties of the convolution kernel and its Fourier transform. Recent developments in the analysis of this class of algorithms indicate that sometimes, as in the case of certain anisotropic curvature ows arising in materials science, these properties of the convolution kernel cannot be expected. Other applications, such as machine learning, would benefit from as great a level of exibility in choosing the convolution kernel as possible. In this paper, we establish certain desirable properties of threshold dynamics, such as gamma convergence of its associated energy, for a substantially wider class of kernels than has been hitherto possible. We also present variants of the algorithm that extend some of these properties to even wider classes of convolution kernels. [ABSTRACT FROM AUTHOR]

The roundoff-error-free (REF) LU and Cholesky factorizations, combined with the REF substitution algorithms, allow rational systems of linear equations to be solved exactly and efficiently by working entirely in integer arithmetic. These REF computational tools share two key properties: their constituent divisions are exact, and their matrix entries' bit-lengths (i.e., required number of storage bits) are bounded polynomially. However, as this paper explains, updating the REF factorizations via the common delete-insert-reduce approach turns out to be inefficient in terms of operand bit-length growth and increased computational effort. In fact, we prove that this inefficiency can eliminate the computational savings expected of the factorization update process. Hence, the current work develops REF update algorithms that differ significantly from their traditional counterparts. To achieve this, it introduces the frame matrix methodology for representing and applying the REF factorizations, and then it devises new REF operations and shortcuts for transitioning between "nearby" frame matrices. The presented push-and-swap update approach utilizes these transitions to preserve the special structure of the REF factorizations and to avoid additional growth in their matrix entries' bit-lengths; this update process requires the expected O(n²) operations. The featured updates are column addition, deletion, and replacement with respect to the REF LU factorization. We also prove that analogous row updates can be performed via the column updates and discuss special considerations for updating the REF Cholesky factorization. An accompanying set of computational tests attests that for fully dense random instances, the push-and-swap column replacement update for the REF LU factorization is on average an order of magnitude faster than the corresponding delete-insert-reduce update for the exact rational arithmetic LU factorization. Altogether, the REF update algorithms represent a promising preliminary step toward enhancing the basic solution validation subroutines employed by state-of-the-art mixed-precision methodologies for exact linear and mixed-integer programming. [ABSTRACT FROM AUTHOR]

We seek discrete approximations to solutions u : Ω → ℝ of semilinear elliptic PDE of the form Δu + fs (u) = 0, where fs is a one-parameter family of nonlinear functions and Ω is a domain in ℝd. The main achievement of this paper is the approximation of solutions to the PDE on the cube Ω = (0,π)³ ⊆ ℝ³. There are 323 possible isotropy subgroups of functions on the cube, which fall into 99 conjugacy classes. The bifurcations with symmetry in this problem are quite interesting, including many with 3-dimensional critical eigenspaces. Our automated symmetry analysis is necessary with so many isotropy subgroups and bifurcations among them, and it allows our code to follow one branch in each equivalence class that is created at a bifurcation point. Our most complicated result is the complete analysis of a degenerate bifurcation with a 6-dimensional critical eigenspace. This article extends our work in [Internat. J. Bifur. Chaos Appl. Sci. Engrg., 19 (2009), pp. 2531-2556], wherein we combined symmetry analysis with modified implementations of the gradient Newton-Galerkin algorithm (GNGA [J. M. Neuberger and J. W. Swift, Internat. J. Bifur. Chaos Appl. Sci. Engrg., 11 (2001), pp. 801-820]) to automatically generate bifurcation diagrams and solution graphics for small, discrete problems with large symmetry groups. The code described in the current paper is efficiently implemented in parallel, allowing us to investigate a relatively fine-mesh discretization of the cube. We use the methodology and corresponding library presented in our previous paper [Int. J. Parallel Program., 40 (2012), pp. 443-464]. [ABSTRACT FROM AUTHOR]

We present two iterative algorithms for solving real nonsymmetric and complex nonHermitian 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. [ABSTRACT FROM AUTHOR]

This paper presents an efficient multiscale stochastic framework for uncertainty quantification in modeling of flow through porous media with multiple rock types. The governing equations are based on Darcy's law with nonstationary stochastic permeability represented as a sum of local Karhunen-Loève expansions. The approximation uses stochastic collocation on either a tensor product or a sparse grid, coupled with a domain decomposition algorithm known as the multiscale mortar mixed finite element method. The latter method requires solving a coarse scale mortar interface problem via an iterative procedure. The traditional implementation requires the solution of local fine scale linear systems on each iteration. We employ a recently developed modification of this method that precomputes a multiscale flux basis to avoid the need for subdomain solves on each iteration. In the stochastic setting, the basis is further reused over multiple realizations, leading to collocation algorithms that are more efficient than the traditional implementation by orders of magnitude. Error analysis and numerical experiments are presented. [ABSTRACT FROM AUTHOR]

In this paper we show that recently introduced quantics tensor train (QTT) decomposition call be considered as an algebraic wavelet transform with adaptively determined filters. The main algorithm for obtaining QTT decomposition can be reformulated as a method seeking "good subspaces" or "good bases" and considered as a parameterized transformation of an initial tensor into a sparse tensor. This interpretation allows us to introduce a modification of the tensor train-SVD (TT-SVD) algorithm to make it work in cases where the original algorithm does not work; it results in the new wavelet-like transforms called wavelet tensor train (WTT) transform. Properties of WTT transforms are studied numerically, and a theoretical conjecture on the number of vanishing moments is proposed. It is shown that WTT transforms are orthogonal by construction, and the efficiency of WTT is compared with and often outperforms Daubechies wavelet transforms on certain classes of function-related vectors and matrices. [ABSTRACT FROM AUTHOR]

In this paper, we propose and study the use of alternating direction algorithms for several ℓ1-norm minimization problems arising from sparse solution recovery in compressive sensing, including the basis pursuit problem, the basis pursuit denoising problems of both unconstrained and constrained forms, and others. We present and investigate two classes of algorithms derived from either the primal or the dual form of ℓ1-problems. The construction of the algorithms consists of two main steps: (1) to reformulate an ℓ1-problem into one having blockwise separable objective functions by adding new variables and constraints; and (2) to apply an exact or inexact alternating direction method to the augmented Lagrangian function of the resulting problem. The derived alternating direction algorithms can be regarded as first-order primal-dual algorithms because both primal and dual variables are updated at every iteration. Convergence properties of these algorithms are established or restated when they already exist. Extensive numerical experiments are performed, using randomized partial Walsh Hadamard sensing matrices, to demonstrate the versatility and effectiveness of the proposed approach. Moreover, we present numerical results to emphasize two practically important but perhaps overlooked points: (i) that algorithm speed should be evaluated relative to appropriate solution accuracy; and (ii) that when erroneous measurements possibly exist, the ℓ1-fidelity should generally be preferable to the ℓ2-fidelity. [ABSTRACT FROM AUTHOR]

This paper analyzes and improves the linearized Bregman method for solving the basis pursuit and related sparse optimization problems. The analysis shows that the linearized Bregman method has the exact regularization property; namely, it converges to an exact solution of the basis pursuit problem whenever its smooth parameter a is greater than a certain value. The analysis is based on showing that the linearized Bregman algorithm is equivalent to gradient descent applied to a certain dual formulation. This result motivates generalizations of the algorithm enabling the use of gradient-based optimization techniques such as line search, Barzilai--Borwein, limited memory BFGS (L-BFGS), nonlinear conjugate gradient, and Nesterov's methods. In the numerical simulations, the two proposed implementations, one using Barzilai--Borwein steps with nonmonotone line search and the other using L-BFGS, gave more accurate solutions in much shorter times than the basic implementation of the linearized Bregman method with a so-called kicking technique. [ABSTRACT FROM AUTHOR]

This work arises from the purpose of applying new tools in dynamical systems to time problems in biological systems. The main aim of this paper is to develop a numerical method to perform the effective computation of the phase advancement when we stimulate an oscillator which has not yet reached the asymptotic state (a limit cycle). That is, we want to extend the computation of the phase resetting curves (PRCs) (the classical tool to compute the phase advancement) to a neighborhood of the limit cycle, obtaining what we call the phase resetting surfaces (PRSs). To achieve this goal we first perform a careful study of the theoretical grounds (the parameterization method for invariant manifolds and another approach using Lie symmetries), which allows us to describe the isochronous sections of the limit cycle and, from them, to obtain the PRSs. In order to make this theoretical framework applicable, we use the numerical algorithms of the parameterization method and other semianalytical tools to extend invariant manifolds; as a result, we design a numerical scheme to compute both the isochrons and the PRSs of a given oscillator. Finally, to illustrate this algorithm, we apply it to some well-known biological models and we include a discussion on different biological and numerical aspects suggested by these examples. [ABSTRACT FROM AUTHOR]

In Part II of this paper we first refine the analysis of error-free vector transformations presented in Part I. Based on that we present an algorithm for calculating the rounded-to-nearest result of s := Σ pi for a given vector of floating-point numbers pi, as well as algorithms for directed rounding. A special algorithm for computing the sign of s is given, also working for huge dimensions. Assume a floating-point working precision with relative rounding error unit eps. We define and investigate a K-fold faithful rounding of a real number r. Basically the result is stored in a vector Resν of K nonoverlapping floating-point numbers such that Σ Resν approximates r with relative accuracy epsK, and replacing ResK by its floating-point neighbors in Σ Resν forms a lower and upper bound for r. For a given vector of floating-point numbers with exact sum s, we present an algorithm for calculating a K-fold faithful rounding of s using solely the working precision. Furthermore, an algorithm for calculating a faithfully rounded result of the sum of a vector of huge dimension is presented. Our algorithms are fast in terms of measured computing time because they allow good instruction-level parallelism, they neither require special operations such as access to mantissa or exponent, they contain no branch in the inner loop, nor do they require some extra precision. The only operations used are standard floating-point addition, subtraction, and multiplication in one working precision, for example, double precision. Certain constants used in the algorithms are proved to be optimal. [ABSTRACT FROM AUTHOR]

We describe a class of logically rectangular quadrilateral and hexahedral grids for solving PDEs in circular and spherical domains, including grid mappings for the circle, the surface of the sphere, and the three-dimensional ball. The grids are logically rectangular and the computational domain is a single Cartesian grid. Compared to alternative approaches based on a multiblock data structure or unstructured triangulations, this approach simplifies the implementation of numerical methods and the use of adaptive refinement. A more general domain with a smooth boundary can be gridded by composing one of the mappings from this paper with another smooth mapping from the circle or sphere to the desired domain. Although these grids are highly nonorthogonal, we show that the high-resolution wave-propagation algorithm implemented in CLAWPACK can be used effectively to approximate hyperbolic problems on these grids. Since the ratio between the largest and smallest grids is below 2 for most of our grid mappings, explicit finite volume methods such as the wave-propagation algorithm do not suffer from the center or pole singularities that arise with polar or latitude-longitude grids. Numerical test calculations illustrate the potential use of these grids for a variety of applications including Euler equations, shallow water equations, and acoustics in a heterogeneous medium. Pattern formation from a reaction-diffusion equation on the sphere is also considered. All examples are implemented in the CLAWPACK software package and full source code is available on the web, along with MATLAB routines for the various mappings. [ABSTRACT FROM AUTHOR]

As a solution to the frequently reported problem of the high frequencies reconstructed sometimes many iterations after the low frequencies in tomography, we proposed in a recent paper [V. Israel-Jost, Ph. Choquet, and A. Constantinesco, Int. J. Biomed. Imaging, 2006, paper 34043.] a general adaptation that can apply to many algorithms: the incomplete backprojection. We stressed the fact that this type of frequency adapted (FA) algorithm works only when the linear problem of tomography to be solved involves deblurring, i.e., when the point spread functions (PSFs) associated with each voxel are large. Thus, depending on the desired level of detail, a parameter is set to achieve the reconstruction within a small number of iterations. The aim of this paper is to give a mathematical analysis of both the acceleration obtained and the convergence of an FA algorithm derived from the simultaneous algebraic reconstruction technique (SART) algorithm of Andersen and Kak, the so-called FA-SART algorithm. [ABSTRACT FROM AUTHOR]

This paper is concerned with the development of well-balanced high order numerical schemes for systems of balance laws with a linear flux function, whose coefficients may be variable. First, well-balanced first order numerical schemes are obtained based on the use of exact solvers of Riemann problems that include both the flux and the source terms. Godunov's methods so obtained are extended to higher order schemes by using a technique of reconstruction of states. The main contribution of this paper is to introduce a reconstruction technique that preserves the well-balanced property of Godunov's methods. Some numerical experiments are presented to verify in practice the properties of the developed numerical schemes. [ABSTRACT FROM AUTHOR]

In this paper we consider the numerical solution of discrete boundary integral equations on polyhedral surfaces in three dimensions. When the solution contains typical edge singularities, highly stretched meshes are preferred to uniform meshes, since they reduce the number of degrees of freedom needed to obtain a fixed accuracy. The classical panel-clustering method can still be applied in the presence of such highly stretched meshes. However, we will show that the savings in computation time and storage become suboptimal because the nearfield matrix arising in the panel-clustering algorithm is no longer as sparse as it is in the case of uniform meshes. Hence, a natural question arises as to whether a new enhanced panel-clustering algorithm can be designed which performs efficiently even in the presence of highly stretched meshes. The main result of this paper is to formulate such an enhanced version of the panel-clustering algorithm. The key features of the algorithm are (i) the employment of partial analytic integration in the direction of stretching, yielding a new kernel function on a one-dimensional manifold where the influence of high aspect ratios in the stretched elements is removed, and (ii) the introduction of a generalized admissibility condition with respect to the partially integrated kernel, which ensures that certain stretched clusters which are inadmissible in the classical sense now become admissible. In the context of a model problem, we prove that our algorithm yields an accurate (up to the discretization error) matrix-vector multiplication which requires ...(NlogκN) operations, where N is the number of degrees of freedom and k is small and independent of the aspect ratio. The generalized admissibility condition can be viewed as an addition to the classical method which may be useful in general when stretched meshes are present. We also have performed a numerical experiment which shows that the sparsity of the nearfield matrix for the enhanced panel-clustering method is not negatively affected by stretched elements, and the method will perform optimally. [ABSTRACT FROM AUTHOR]

This paper approaches optimal control problems for discrete-time controlled Markov processes by representing the value of the problem in a dual Lagrangian form; the value is expressed as an infimum over a family of Lagrangian martingales of an expectation of a pathwise supremum of the objective adjusted by the Lagrangian martingale term. This representation opens up the possibility of numerical methods based on Monte Carlo simulation, which may be advantageous in high-dimensional problems or in problems with complicated constraints. [ABSTRACT FROM AUTHOR]

In this paper, we show how to perform the Hessenberg reduction of a rank structured matrix under unitary similarity operations in an efficient way, using the Givens-weight representation. This reduction can be used as a first step for eigenvalue computation. We also show how the algorithm can be modified to compute the bidiagonal reduction of a rank structured matrix, this latter method being a preprocessing step for computing the singular values of the matrix. For the main cases of interest, the algorithms we describe in this paper are of complexity O((ar + bs)n2), where n is the matrix size, r is some measure for the average rank index of the rank structure, s is some measure for the bandwidth of the unstructured matrix part around the main diagonal, and a, b ? R are certain weighting parameters. Numerical experiments demonstrate the stability of this approach. [ABSTRACT FROM AUTHOR]

In this paper we study a discrete version of the classical team decision problem. It has been shown previously that the general discrete team decision problem is NP-hard. Here we present an efficient approximation algorithm for this problem. For the maximization version of this problem with nonnegative rewards, this algorithm computes decision rules which are guaranteed to be within a fixed bound of optimal. [ABSTRACT FROM AUTHOR]

In this paper a new stochastic algorithm for the solution of population balance equations is presented. The population balance equations have the form of extended Smoluchowski equations which include linear and source terms. The new algorithm, called the linear process deferment algorithm (LPDA), is used for solving a detailed model describing the formation of soot in premixed laminar flames. A measure theoretic formulation of a stochastic jump process is developed and the corresponding generator presented. The numerical properties of the algorithm are studied in detail and compared to the direct simulation algorithm and various splitting techniques. LPDA is designed for all kinds of population balance problems where nonlinear processes cannot be neglected but are dominated in rate by linear ones. In those cases the LPDA is seen to reduce run times for a population balance simulation by a factor of up to 1000 with a negligible loss of accuracy. [ABSTRACT FROM AUTHOR]

This paper uses periodic matrix multiplication to improve the time complexities for a number of graph problems. The time for finding an asteroidal triple is reduced from O(nm) to O(n2.82), and the time for finding a star cutset, a two-pair, and a dominating pair is reduced from O(nm) to O(n2.79). It is also shown that each of these problems is at least as hard as one of three basic graph problems for which the best known algorithms run in times O(nm) and O(nα). We note that the fast matrix multiplication algorithms do not seem to be practical because of the enormous constants needed to achieve the asymptotic time bounds. These results are important theoretically for breaking the n³ barrier rather than giving efficient algorithms for a user. [ABSTRACT FROM AUTHOR]

We propose a new algebraic multigrid (AMG) method for solving the eddy current approximations to Maxwell's equations. This AMG method has its roots in an algorithm proposed by Reitzinger and Schöberl. The main focus in the Reitzinger and Schöberl method is to maintain null-space properties of the weak ∇ × ∇ × operator on coarse grids. While these null-space properties are critical, they are not enough to guarantee h-independent convergence rates of the overall multigrid scheme. We present a new strategy for choosing intergrid transfers that not only maintains the important null-space properties on coarse grids but also yields significantly improved multigrid convergence rates. This improvement is related to those we explored in a previous paper, but is fundamentally simpler, easier to compute, and performs better with respect to both multigrid operator complexity and convergence rates. The new strategy builds on ideas in smoothed aggregation to improve the approximation property of an existing interpolation operator. By carefully choosing the smoothing operators, we show how it is sometimes possible to achieve h-independent convergence rates with a modest increase in multigrid operator complexity. Though this ideal case is not always possible, the overall algorithm performs significantly better than the original scheme in both iterations and run time. Finally, the Reitzinger and Schöberl method, as well as our previous smoothed method, are shown to be special cases of this new algorithm. [ABSTRACT FROM AUTHOR]

This paper presents a multigrid method for solving variable coefficient Maxwell's equations. The novelty in this method is the use of interpolation operators that do not produce multilevel commutativity complexes that lead to multilevel exactness. Rather, the effects of multilevel exactness are built into the level equations themselves—on the finest level using a discrete T – V formulation and on the coarser grids through the Galerkin coarsening procedure of a T – V formulation. These built-in structures permit the levelwise use of an effective hybrid smoother on the curl-free near-nullspace components and permit the development of interpolation operators for handling the curl-free and divergence-free error components separately. The resulting block-diagonal interpolation operator does not satisfy multilevel commutativity but has good approximation properties for both of these error components. Applying operator-dependent interpolation for each of these error components leads to an effective multigrid scheme for variable coefficient Maxwell's equations, where multilevel commutativity-based methods can degrade. Numerical results are presented to verify the effectiveness of this new scheme. [ABSTRACT FROM AUTHOR]

In this paper, an efficient and stable numerical algorithm for computing multiple saddle points with symmetries is developed by modifying the local minimax method established in [Y. Li and J. Zhou, SIAM J. Sci. Comput. 23 (2001), pp. 840--865; Y. Li and J. Zhou, SIAM J. Sci. Comput., 24 (2002), pp. 840--865]. First an invariant space is defined in a more general sense and a principle of invariant criticality is proved for the generalization. Then the orthogonal projection to the invariant space is used to preserve the invariance and to reduce computational error across iterations. Simple averaging formulas are used for the orthogonal projections. Numerical computations of examples with various symmetries, of which some can and others cannot be characterized by a compact group of linear isomorphisms, are carried out to confirm the theory and to illustrate applications. The mathematical features of various problems demonstrated in these examples fall into two categories: nodal solutions of saddle-point type with large Morse indices and nonradial positive solutions via symmetry breaking in radially symmetric elliptic problems. The new numerical algorithm generates these rather unstable solutions in an efficient and stable way. The existence of many unstable solutions and their behavior found in this paper remain to be investigated. [ABSTRACT FROM AUTHOR]

Local mesh-refining algorithms known from adaptive finite element methods are adopted for locally conservative and monotone finite volume discretizations of boundary value problems for steady-state convection-diffusion-reaction equations. The paper establishes residual-type explicit error estimators and averaging techniques for a posteriori finite volume error control with and without upwind in global H1- and L2-norms. Reliability and efficiency are verified theoretically and confirmed empirically with experimental support for the superiority of the suggested adaptive mesh-refining algorithms over uniform mesh refining. A discussion of adaptive computations in the simulation of contaminant concentration in a nonhomogeneous water reservoir concludes the paper. [ABSTRACT FROM AUTHOR]

For several graph-theoretic parameters such as vertex cover and dominating set, it is known that if their sizes are bounded by k, then the treewidth of the graph is bounded by some function of k. This fact is used as the main tool for the design of several fixed-parameter algorithms on minor-closed graph classes such as planar graphs, single-crossing-minor-free graphs, and graphs of bounded genus. In this paper we examine whether similar bounds can be obtained for larger minor-closed graph classes and for general families of graph parameters, including all those for which such behavior has been reported so far. Given a graph parameter P, we say that a graph family F has the parameter-treewidth property for P if there is an increasing function t such that every graph G ∈ F has treewidth at most t(P(C)). We prove as our main result that, for a large family of graph parameters called contraction-bidimensional, a minor-closed graph family F has the parameter-treewidth property if F has bounded local treewidth. We also show "if and only if" for some graph parameters, and thus, this result is in some sense tight. In addition we show that, for a slightly smaller family of graph parameters called minor-bidimensional, all minor-closed graph families F, excluding some fixed graphs, have the parameter-treewidth property. The contraction-bidimensional parameters include many domination and covering graph parameters such as vertex cover, feedback vertex set, dominating set, edge-dominating set, and q-dominating set (for fixed q). We use our theorems to develop new fixed-parameter algorithms in these contexts. [ABSTRACT FROM AUTHOR]

This paper provides a nice example to illustrate that without supporting analysis, a numerical simulation may lead to incorrect conclusions. We explore a pedagogical example of failure of Runge-Kutta-Felberg (RKF) algorithms for a simple dynamical system that models the coupling of two oscillators. Although the system appears to be well-behaved, the explicit RKF solvers provide erratic numerical solutions. The mode of failure is based in a period-doubling route to chaos due to the existence of stable linear solutions in the problem. [ABSTRACT FROM AUTHOR]

The main purpose of this paper is to give numerical algorithms and the error analysis for delay quadratic problems in the calculus of variations. These methods are new, efficient, and accurate and have a global a priori error of O(h²), where h is the distance between any two successive node points. We also derive the results for the general, numerical, delay problem, but we focus on proving our results in the simpler quadratic case since the extra technical numerical details have been given previously by the second author. In addition, the authors have previously shown how to reformulate general delay constrained problems in optimal control theory/constrained calculus of variations as unconstrained delay problems. Thus our numerical results and methods will hold for these general constrained problems also. Finally, we note that our algorithm, which solves the stationary condition(s) numerically, avoids the more difficult problems of piecing the solution of second order equations together and requires less smoothness in the solution. Thus we replace difficult second order boundary value problems with the easier task of approximating definite integrals involving first order derivatives. [ABSTRACT FROM AUTHOR]

In this paper we show how to compute recursively an approximation of the left and right dominant singular subspaces of a given matrix. In order to perform as few as possible operations on each column of the matrix, we use a variant of the classical Gram­Schmidt algorithm to estimate this subspace. The method is shown to be particularly suited for matrices with many more rows than columns. Bounds for the accuracy of the computed subspace are provided. Moreover, the analysis of error propagation in this algorithm provides new insights in the loss of orthogonality typically observed in the classical Gram­Schmidt method. [ABSTRACT FROM AUTHOR]

This paper presents numerical results for the asynchronous version of the fast adaptive composite-grid algorithm (AFACx). These results confirm the level-independent convergence bounds established theoretically in a companion paper. These numerical results include the case of AFACx applied to first-order system least-squares finite element discretizations of the stationary Stokes equations on curvilinear adaptive mesh refinement grids. [ABSTRACT FROM AUTHOR]

In an earlier paper, we developed a new formulation of the shallow-water equations on the sphere, based on the semi-Lagrangian advection of potential vorticity. An important feature of our two-time-level model is that forward extrapolation of the nonlinear terms is avoided. In this paper, we develop a multigrid algorithm for solving the resulting system of three coupled nonlinear equations. We include a discussion of basic multigrid methods, along with some more advanced techniques required for the development of the multigrid scheme for this system. Results are presented showing that usual optimal multigrid efficiency is obtained. [ABSTRACT FROM AUTHOR]

Dynamical systems can often be decomposed into loosely coupled subsystems. The system of ordinary differential equations (ODEs) modelling such a problem can then be partitioned corresponding to the subsystems, and the loose couplings can be exploited by special integration methods to solve the problem using a parallel computer or just solve the problem more efficiently than by standard methods. This paper presents accuracy analysis of methods for the numerical integration of stiff partitioned systems of ODEs. The discretization formulas are based on the implicit Euler formula and the second order implicit backward differentiation formula (BDF2). Each subsystem of the partitioned problem is discretized independently, and the couplings to the other subsystems are based on solution values from previous time steps. Applied this way, the discretization formulas are called decoupled. The stability properties of the decoupled implicit Euler formula are well understood. This paper presents error bounds and asymptotic error expansions to be used in controlling step size, relaxation between subsystems and the validity of the partitioning. The decoupled BDF2 formula is analyzed within the same framework. Finally, the analysis is used in the design of a decoupled numerical integration algorithm with variable step size to control the local error and adaptive selection of partitionings. Two versions of the algorithm with decoupled implicit Euler and BDF2, respectively, are used in examples where a realistic problem is solved. The examples compare the results from the decoupled implicit Euler and BDF2 formulas, and compare them with results from the corresponding classical formulas. [ABSTRACT FROM AUTHOR]

In this paper we show that the use of spectral Galerkin methods for the approximation of the Boltzmann equation in the velocity space permits one to obtain spectrally accurate numerical solutions at a reduced computational cost. We prove that the spectral algorithm preserves the total mass and approximates with infinite-order accuracy momentum and energy. Consistency of the method is also proved, and a stability result for a smoothed positive scheme is given. We demonstrate that the Fourier coefficients associated with the collision kernel of the equation have a very simple structure and in some cases can be computed explicitly. Numerical examples for homogeneous test problems in two and three dimensions confirm the advantages of the method. [ABSTRACT FROM AUTHOR]

This paper presents a new algorithm for accurate floating-point computation of the singular value decomposition (SVD) of the product A = BτSC, where B ∈ Rp×m, C ∈ Rq×n, S ∈ Rp×q, and p ≤ m, q ≤ n. The new algorithm uses diagonal scalings, the LU factorization with complete pivoting, the QR factorization with column pivoting, and matrix multiplication to replace A by A′ = B′τS′C′, where A and A′ have the same singular values and the matrix A′ is computed explicitly. The singular values of A′ are computed using the Jacobi SVD algorithm. It is shown that the accuracy of the new algorithm is determined by (i) the accuracy of the QR factorizations of Bτ and Cτ; (ii) the accuracy of the LU factorization with complete pivoting of S; and (iii) the accuracy of the computation of the SVD of a matrix A′ with moderate minD=diagκ2(A′D). Theoretical analysis and numerical evidence show that, in the case of rank(B) = rank(C) = p and full rank S, the accuracy of the new algorithm is unaffected by replacing B, S, C with, respectively, D1B, D2SD3, D4C, where Di, i = 1,…,4, are arbitrary diagonal matrices. As an application, the paper proposes new accurate algorithms for computing the (H,K)­SVD and (H-1,K)­SVD of S. [ABSTRACT FROM AUTHOR]

The development of strong linearizations preserving whatever structure a matrix polynomial might possess has been a very active area of research in the last years, since such linearizations are the starting point of numerical algorithms for computing eigenvalues of structured matrix polynomials with the properties imposed by the considered structure. In this context, Hermitian matrix polynomials are one of the most important classes of matrix polynomials arising in applications and their real eigenvalues are of great interest. The sign characteristic is a set of signs attached to these real eigenvalues which is crucial for determining the behavior of systems described by Hermitian matrix polynomials and, therefore, it is desirable to develop linearizations that preserve the sign characteristic of these polynomials, but, at present, only one such linearization is known. In this paper, we present a complete characterization of all the Hermitian strong linearizations that preserve the sign characteristic of a given Hermitian matrix polynomial and identify several families of such linearizations that can be constructed very easily from the coefficients of the polynomial. [ABSTRACT FROM AUTHOR]

Solving feasibility problems is a central task in mathematics and the applied sciences. One particularly successful method is the Douglas--Rachford algorithm. In this paper, we provide many new conditions sufficient for finite convergence. Numerou s examples illustrate our results. [ABSTRACT FROM AUTHOR]