Search

We review the basic outline of the highly successful diffusion Monte Carlo technique commonly used in contexts ranging from electronic structure calculations to rare event simulation and data assimilation, and propose a new class of randomized iterative algorithms based on similar principles to address a variety of common tasks in numerical linear algebra. From the point of view of numerical linear algebra, the main novelty of the Fast Randomized Iteration schemes described in this article is that they work in either linear or constant cost per iteration (and in total, under appropriate conditions) and are rather versatile: we will show how they apply to solution of linear systems, eigenvalue problems, and matrix exponentiation, in dimensions far beyond the present limits of numerical linear algebra. While traditional iterative methods in numerical linear algebra were created in part to deal with instances where a matrix (of size $\mathcal{O}(n^2)$) is too big to store, the algorithms that we propose are effective even in instances where the solution vector itself (of size $\mathcal{O}(n)$) may be too big to store or manipulate. In fact, our work is motivated by recent DMC based quantum Monte Carlo schemes that have been applied to matrices as large as $10^{108} \times 10^{108}$. We provide basic convergence results, discuss the dependence of these results on the dimension of the system, and demonstrate dramatic cost savings on a range of test problems., 44 pages, 7 figures

We discuss and compare two greedy algorithms that compute discrete versions of Fekete-like points for multivariate compact sets by basic tools of numerical linear algebra. The first gives the so-called approximate Fekete points by QR factorization with column pivoting of Vandermonde-like matrices. The second computes discrete Leja points by LU factorization with row pivoting. Moreover, we study the asymptotic distribution of such points when they are extracted from weakly admissible meshes.

In 1981 Hong and Kung proved a lower bound on the amount of communication (amount of data moved between a small, fast memory and large, slow memory) needed to perform dense, n-by-n matrix multiplication using the conventional O(n3) algorithm, where the input matrices were too large to fit in the small, fast memory. In 2004 Irony, Toledo, and Tiskin gave a new proof of this result and extended it to the parallel case (where communication means the amount of data moved between processors). In both cases the lower bound may be expressed as Ω(#arithmetic_operations/M), where M is the size of the fast memory (or local memory in the parallel case). Here we generalize these results to a much wider variety of algorithms, including LU factorization, Cholesky factorization, LDLT factorization, QR factorization, the Gram–Schmidt algorithm, and algorithms for eigenvalues and singular values, i.e., essentially all direct methods of linear algebra. The proof works for dense or sparse matrices and for sequential or para...

Probability one homotopy algorithms are a class of methods for solving nonlinear systems of equations that are globally convergent with probability one. These methods are theoretically powerful, and if constructed and implemented properly, are robust, numerically stable, accurate, and practical. The concomitant numerical linear algebra problems deal with rectangular matrices, and good algorithms require a delicate balance (not always achieved) of accuracy, robustness, and efficiency in both space and time. The author's experience with globally convergent homotopy algorithms is surveyed here, and some of the linear algebra difficulties for dense and sparse problems are discussed.

The existence of parallel and pipeline computers has inspired a new approach to algorithmic analysis. Classical numerical methods are generally unable to exploit multiple processors and powerful vector-oriented hardware. Efficient parallel algorithms can be created by reformulating familiar algorithms or by discovering new ones, and the results are often surprising. A comprehensive survey of parallel techniques for problems in linear algebra is given. Specific topics include: relevant computer models and their consequences for programs, evaluation of arithmetic expressions, solution of general and special linear systems of equations, and computation of eigenvalues.

International audience; Traditional rounding error analysis in numerical linear algebra leads to backward error bounds involving the constant γ n = nu/(1 − nu), for a problem size n and unit roundoff u. In the light of large-scale and possibly low-precision computations, such bounds can struggle to provide any useful information. We develop a new probabilistic rounding error analysis that can be applied to a wide range of algorithms. By using a concentration inequality and making probabilistic assumptions about the rounding errors, we show that in several core linear algebra computations γ n can be replaced by a relaxed constant γ n proportional to √ n log n u with a probability bounded below by a quantity independent of n. The new constant γ n grows much more slowly with n than γn. Our results have three key features: they are backward error bounds; they are exact, not first order; and they are valid for any n, unlike results obtained by applying the central limit theorem, which apply only as n → ∞. We provide numerical experiments that show that for both random and real-life matrices the bounds can be much smaller than the standard deterministic bounds and can have the correct asymptotic growth with n. We also identify two special situations in which the assumptions underlying the analysis are not valid and the bounds do not hold. Our analysis provides, for the first time, a rigorous foundation for the rule of thumb that "one can take the square root of an error constant because of statistical effects in rounding error propagation".

We present first-order perturbation analysis of a simple eigenvalue and the corresponding right and left eigenvectors of a general square matrix, not assumed to be Hermitian or normal. The eigenvalue result is well known to a broad scientific community. The treatment of eigenvectors is more complicated, with a perturbation theory that is not so well known outside a community of specialists. We give two different proofs of the main eigenvector perturbation theorem. The first, a block-diagonalization technique inspired by the numerical linear algebra research community and based on the implicit function theorem, has apparently not appeared in the literature in this form. The second, based on complex function theory and on eigenprojectors, as is standard in analytic perturbation theory, is a simplified version of well-known results in the literature. The second derivation uses a convenient normalization of the right and left eigenvectors defined in terms of the associated eigenprojector, but although this dates back to the 1950s, it is rarely discussed in the literature. We then show how the eigenvector perturbation theory is easily extended to handle other normalizations that are often used in practice. We also explain how to verify the perturbation results computationally. We conclude with some remarks about difficulties introduced by multiple eigenvalues and give references to work on perturbation of invariant subspaces corresponding to multiple or clustered eigenvalues. Throughout the paper we give extensive bibliographic commentary and references for further reading.

We show that the spline kernel of order $p$ is a so-called semiseparable function with semiseparability rank $p$. A consequence of this is that kernel matrices generated by the spline kernel are rank structured matrices that can be stored and factorized efficiently. We use this insight to derive new recursive algorithms with linear complexity in the number of knots for various kernel matrix computations. We also discuss applications of these algorithms, including smoothing spline regression, Gaussian process regression, and some related hyperparameter estimation problems.

Quantum embedding theories are playing an increasingly important role in bridging different levels of approximation to the many-body Schrodinger equation in physics, chemistry, and materials scienc...

Every student of numerical linear algebra is familiar with block matrices and vectors. The same ideas can be applied to the continuous analogues of operators, functions, and functionals. It is shown here how the explicit consideration of block structures at the continuous level can be a useful tool. In particular, block operator diagrams lead to templates for spectral discretization of differential and integral equation boundary-value problems in one space dimension by the rectangular differentiation, identity, and integration matrices introduced recently by Driscoll and Hale. The templates are so simple that we are able to present them as executable Matlab codes just a few lines long, developing ideas through a sequence of 12 increasingly advanced examples. The notion of the rectangular shape of a linear operator is made mathematically precise by the theory of Fredholm operators and their indices, and the block operator formulations apply to nonlinear problems too. We propose the convention of representing nonlinear blocks as shaded. At each step of a Newton iteration for a nonlinear problem, the structure is linearized and the blocks become unshaded, representing Fréchet derivative operators, square or rectangular.

The generalized singular value decomposition (GSVD) is one of the essential tools in numerical linear algebra. This paper proposes a regularization method, combining Tikhonov regularization in general form with the truncated GSVD. Then the randomized algorithms are adopted to implement the truncation process. This randomized GSVD for the regularization of the large-scale ill-posed problems can achieve good accuracy with less computational time and memory requirement than the classical regularization methods. Finally, we present the error analyses for the randomized algorithms. Some illustrative numerical examples are provided.

Finding a fixed point to a nonexpansive operator, i.e., $x^*=Tx^*$, abstracts many problems in numerical linear algebra, optimization, and other areas of scientific computing. To solve fixed-point problems, we propose ARock, an algorithmic framework in which multiple agents (machines, processors, or cores) update $x$ in an asynchronous parallel fashion. Asynchrony is crucial to parallel computing since it reduces synchronization wait, relaxes communication bottleneck, and thus speeds up computing significantly. At each step of ARock, an agent updates a randomly selected coordinate $x_i$ based on possibly out-of-date information on $x$. The agents share $x$ through either global memory or communication. If writing $x_i$ is atomic, the agents can read and write $x$ without memory locks. Theoretically, we show that if the nonexpansive operator $T$ has a fixed point, then with probability one, ARock generates a sequence that converges to a fixed points of $T$. Our conditions on $T$ and step sizes are weaker than comparable work. Linear convergence is also obtained. We propose special cases of ARock for linear systems, convex optimization, machine learning, as well as distributed and decentralized consensus problems. Numerical experiments of solving sparse logistic regression problems are presented., updated the linear convergence proofs

We design a distributed-memory randomized structured multifrontal solver for large sparse matrices. Two layers of hierarchical tree parallelism are used. A sequence of innovative parallel methods are developed for randomized structured frontal matrix operations, structured update matrix computation, skinny extend-add operation, selected entry extraction from structured matrices, etc. Several strategies are proposed to reuse computations and reduce communications. Unlike an earlier parallel structured multifrontal method that still involves large dense intermediate matrices, our parallel solver performs the major operations in terms of skinny matrices and fully structured forms. It thus significantly enhances the efficiency and scalability. Systematic communication cost analysis shows that the numbers of words are reduced by factors of about $O(\sqrt{n}/r)$ in two dimensions and about $O(n^{2/3}/r)$ in three dimensions, where $n$ is the matrix size and $r$ is an off-diagonal numerical rank bound of the int...

Several innovative random-sampling and random-mixing techniques for solving problems in linear algebra have been proposed in the last decade, but they have not yet made a significant impact on numerical linear algebra. We show that by using a high-quality implementation of one of these techniques, we obtain a solver that performs extremely well in the traditional yardsticks of numerical linear algebra: it is significantly faster than high-performance implementations of existing state-of-the-art algorithms, and it is numerically backward stable. More specifically, we describe a least-squares solver for dense highly overdetermined systems that achieves residuals similar to those of direct QR factorization-based solvers (lapack), outperforms lapack by large factors, and scales significantly better than any QR-based solver.

Appearing frequently in applications, generalized eigenvalue problems represent one of the core problems in numerical linear algebra. The QZ algorithm of Moler and Stewart is the most widely used algorithm for addressing such problems. Despite its importance, little attention has been paid to the parallelization of the QZ algorithm. The purpose of this work is to fill this gap. We propose a parallelization of the QZ algorithm that incorporates all modern ingredients of dense eigensolvers, such as multishift and aggressive early deflation techniques. To deal with (possibly many) infinite eigenvalues, a new parallel deflation strategy is developed. Numerical experiments for several random and application examples demonstrate the effectiveness of our algorithm on two different distributed memory HPC systems.

We show that if the nearly-linear time solvers for Laplacian matrices and their generalizations can be extended to solve just slightly larger families of linear systems, then they can be used to quickly solve all systems of linear equations over the reals. This result can be viewed either positively or negatively: either we will develop nearly-linear time algorithms for solving all systems of linear equations over the reals, or progress on the families we can solve in nearly-linear time will soon halt.

In order to compute polarization energy of biomolecules, we describe a boundary element approach to solving the linearized Poisson-Boltzmann equation. Our approach combines several important features, including the derivative boundary formulation of the problem and a smooth approximation of the molecular surface based on the algebraic spline molecular surface. State of the art software for numerical linear algebra and the kernel independent fast multipole method is used for both simplicity and efficiency of our implementation. We perform a variety of computational experiments, testing our method on a number of actual proteins involved in molecular docking and demonstrating the effectiveness of our solver for computing molecular polarization energy.

We consider low-rank reconstruction of a matrix using its columns and we present asymptotically optimal algorithms for both spectral norm and Frobenius norm reconstruction. The main tools we introduce to obtain our r esults are: (i) the use of fast approximate SVD-like decompositions for column reconstruction, and (ii) two deter ministic algorithms for selecting rows from matrices with orthonormal columns, building upon the sparse represen tation theorem for decompositions of the identity that appeared in \cite{BSS09}., SIAM Journal on Computing (SICOMP), invited to special issue of FOCS 2011

In the mathematical modeling of real-life applications, systems of equations with complex coefficients often arise. While many techniques of numerical linear algebra, e.g., Krylov-subspace methods, extend directly to the case of complex-valued matrices, some of the most effective preconditioning techniques and linear solvers are limited to the real-valued case. Here, we consider the extension of the popular algebraic multigrid method to such complex-valued systems. The choices for this generalization are motivated by classical multigrid considerations, evaluated with the tools of local Fourier analysis, and verified on a selection of problems related to real-life applications.

The upwind discontinuous Galerkin method is an attractive method for solving time-dependent hyperbolic conservation laws. It is possible to use high-order explicit time-stepping methods and high-order spatial approximations without incurring heavy numerical linear algebra overheads. However, the Courant-Friedrichs-Lewy (CFL) condition for these methods depends on the polynomial order used, and there is a somewhat excessive cost for using very high order spatial approximation. We discuss the impact of a covolume mesh based filter on the CFL number for these methods and present an algorithm which has a CFL number independent of the spatial order of approximation. We present computational results for the advection equation and the wave equation on one-dimensional meshes using up to tenth order in space and time.

The growth factor of a matrix quantifies the amount of potential error growth possible when a linear system is solved using Gaussian elimination with row pivoting. While it is an easy matter [N. J. Higham and D. J. Higham, SIAM J. Matrix Anal. Appl., 10 (1989), pp. 155-164] to construct examples of $n\times n$ matrices having any growth factor up to the maximum of $2^{n-1}$, the weight of experience and analysis [N. J. Higham, Accuracy and Stability of Numerical Algorithms, SIAM, Philadelphia, 1996], [L. N. Trefethen and R. S. Schreiber, SIAM J. Matrix Anal. Appl., 11 (1990), pp. 335-360], [L. N. Trefethen and I. D. Bau, Numerical Linear Algebra, SIAM, Philadelphia, 1997] suggest that matrices with exponentially large growth factors are exceedingly rare. Here we show how to conduct numerical experiments on random matrices using a multicanonical Monte Carlo method to explore the tails of growth factor probability distributions. Our results suggest, for example, that the occurrence of an $8\times 8$ matrix with a growth factor of 40 is on the order of a once-in-the-age-of-the-universe event.

We consider the problem of performing accurate computations with rectangular $(m\times n)$ totally nonnegative matrices. The matrices under consideration have the property of having a unique representation as products of nonnegative bidiagonal matrices. Given that representation, one can compute the inverse, LDU decomposition, eigenvalues, and SVD of a totally nonnegative matrix to high relative accuracy in $\mathcal{O}(\max(m^3,n^3))$ time—much more accurately than conventional algorithms that ignore that structure. The contribution of this paper is to show that the high relative accuracy is preserved by operations that preserve the total nonnegativity—taking a product, re-signed inverse (when $m=n$), converse, Schur complement, or submatrix of a totally nonnegative matrix, any of which costs at most $\mathcal{O}(\max(m^3,n^3))$. In other words, the class of totally nonnegative matrices for which we can do numerical linear algebra very accurately in $\mathcal{O}(\max(m^3,n^3))$ time (namely, those for which we have a product representation via nonnegative bidiagonals) is closed under the operations listed above.

We describe how classical Floquet theory may be utilized, in a continuation framework, to construct an efficient Fourier spectral algorithm for approximating periodic orbits. At each continuation step, only a single square matrix, whose size equals the dimension of the phase-space, needs to be factorized; the rest of the required numerical linear algebra just consists of back-substitutions with this matrix. The eigenvalues of this key matrix are the Floquet exponents, whose crossing of the imaginary axis indicates bifurcation and change-in-stability. Hence we also describe how the new periodic orbits created at a period-doubling bifurcation point may be efficiently computed using our approach.

The stable similarity reduction of a nonsymmetric square matrix to tridiagonal form has been a long-standing problem in numerical linear algebra. The biorthogonal Lanczos process is in principle a candidate method for this task, but in practice it is confined to sparse matrices and is restarted periodically because roundoff errors affect its three-term recurrence scheme and degrade the biorthogonality after a few steps. This adds to its vulnerability to serious breakdowns or near-breakdowns, the handling of which involves recovery strategies such as the look-ahead technique, which needs a careful implementation to produce a block-tridiagonal form with unpredictable block sizes. Other candidate methods, geared generally towards full matrices, rely on elementary similarity transformations that are prone to numerical instabilities. Such concomitant difficulties have hampered finding a satisfactory solution to the problem for either sparse or full matrices. This study focuses primarily on full matrices. After outlining earlier tridiagonalization algorithms from within a general framework, we present a new elimination technique combining orthogonal similarity transformations that are stable. We also discuss heuristics to circumvent breakdowns. Applications of this study include eigenvalue calculation and the approximation of matrix functions.

An object-oriented MATLAB system is described for performing numerical linear algebra on continuous functions and operators rather than the usual discrete vectors and matrices. About eighty MATLAB functions from plot and sum to svd and cond have been overloaded so that one can work with our "chebfun" objects using almost exactly the usual MATLAB syntax. All functions live on [-1,1] and are represented by values at sufficiently many Chebyshev points for the polynomial interpolant to be accurate to close to machine precision. Each of our overloaded operations raises questions about the proper generalization of familiar notions to the continuous context and about appropriate methods of interpolation, differentiation, integration, zerofinding, or transforms. Applications in approximation theory and numerical analysis are explored, and possible extensions for more substantial problems of scientific computing are mentioned.

We propose a new interpretation of the generalized overdetermined eigenvalue problem $\left({\bf A} - \lambda {\bf B}\right){\bf v} \approx 0$ for two $m\times n$ ($m > n$) matrices ${\bf A}$ and ${\bf B}$, its stability analysis, and an efficient algorithm for solving it. Usually, the matrix pencil $\{{\bf A} - \lambda {\bf B}\}$ does not have any rank deficient member. Therefore we aim to compute $\lambda$ for which ${\bf A} - \lambda {\bf B}$ is as close as possible to rank deficient; i.e., we search for $\lambda$ that locally minimize the smallest singular value over the matrix pencil $\{{\bf A} - \lambda {\bf B}\}$. Practically, the proposed algorithm requires ${\cal O}(mn^2)$ operations for computing all the eigenpairs. We also describe a method to compute practical starting eigenpairs. The effectiveness of the new approach is demonstrated with numerical experiments. A MATLAB-based implementation of the proposed algorithm can be found at http://www.mat.univie.ac.at/ neum/software/oeig/.

This is the first of two papers on the row by row decoupling problem and the triangular decoupling problem. In this paper we study the row by row decoupling problem with stability in control theory. We first prove a nice reduction property for the row by row decoupling problem with stability and then develop a numerically reliable method for solving it. The basis of our main results is some condensed forms under orthogonal transformations, which can be implemented in numerically stable ways. Hence our results lead to numerically reliable methods for solving the studied problem using existing numerical linear algebra software such as MATLAB. In the sequel [SIAM J. Matrix Anal. Appl., 23 (2002), pp. 1171--1182], we will consider a related problem---the triangular decoupling problem---and parameterize all solutions for it.

Matrix factorization in numerical linear algebra (NLA) typically serves the purpose of restating some given problem in such a way that it can be solved more readily; for example, one major application is in the solution of a linear system of equations. In contrast, within applied statistics/psychometrics (AS/P), a much more common use for matrix factorization is in presenting, possibly spatially, the structure that may be inherent in a given data matrix obtained on a collection of objects observed over a set of variables. The actual components of a factorization are now of prime importance and not just as a mechanism for solving another problem. We review some connections between NLA and AS/P and their respective concerns with matrix factorization and the subsequent rank reduction of a matrix. We note in particular that several results available for many decades in AS/P were more recently (re)discovered in the NLA literature. Two other distinctions between NLA and AS/P are also discussed briefly: how a generalized singular value decomposition might be defined, and the differing uses for the (newer) methods of optimization based on cyclic or iterative projections.

We introduce a methodology for obtaining inventories of error results for families of numerical dense linear algebra algorithms. The approach for deriving the analyses is goal-oriented, systematic, and layered. The presentation places the analysis side-by-side with the algorithm so that it is obvious where roundoff error is introduced. The approach supports the analysis of more complex algorithms, such as the blocked LU factorization. For this operation we derive a tighter error bound than has been previously reported in the literature.

In this paper, the ill-posed problem of computation of the Fourier transform $\hat{f}$ for a noisy bandlimited signal $f$ is discussed. Here, $\hat{f}$ is computed by a regularized Fourier series which is obtained by a regularization of the Fourier transformed sampling formula of Shannon. If the regularization parameter tends to zero, the uniform and $L^2$-convergence of this method is proved. Numerical examples show the performance of this approach.

There is a need for flexible iterative solvers that can solve large-scale ($>10^6$ unknowns) nonsymmetric sparse linear systems to a small tolerance. Among flexible solvers, flexible GMRES (FGMRES) is attractive because it minimizes the residual norm over a particular subspace. In practice, FGMRES is often restarted periodically to keep memory and work requirements reasonable; however, like restarted GMRES, restarted FGMRES can suffer from stagnation. This has led us to develop a flexible variant of the Krylov linear solver GCROT (generalized conjugate residual with inner orthogonalization and outer truncation). Unlike the original GCROT algorithm, the proposed GCROT variant uses a simplified truncation strategy similar to loose GMRES (LGMRES). This modification is motivated by numerical experiments that suggest the specific subspace retained in the outer iteration of GCROT is less important than its size. The flexible GCROT variant appears to be well suited for advection-dominated problems. In particular, when applied to an adjoint problem from computational aerodynamics, the proposed GCROT variant is robust and efficient compared with several popular truncated Krylov subspace methods. Finally, a flexible version of LGMRES is easily constructed by recognizing algorithmic similarities to GCROT.

Given a set of pairwise object distances and a dimension $k$, FastMap and RobustMap algorithms compute a set of $k$-dimensional coordinates for the objects. These metric space embedding methods implicitly assume a higher-dimensional coordinate representation and are a sequence of translations and orthogonal projections based on a sequence of object pair selections (called pivot pairs). We develop a matrix computation viewpoint of these algorithms that operates on the coordinate representation explicitly using Householder reflections. The resulting coordinate mapping algorithm (CMA) is a fast approximate alternative to truncated principal component analysis (PCA), and it brings the FastMap and RobustMap algorithms into the mainstream of numerical computation where standard BLAS building blocks are used. Motivated by the geometric nature of the embedding methods, we further show that truncated PCA can be computed with CMA by specific pivot pair selections. Describing FastMap, RobustMap, and PCA as CMA computations with different pivot pair choices unifies the methods along a pivot pair selection spectrum. We also sketch connections to the semidiscrete decomposition and the QLP decomposition.

In this paper, we describe block matrix algorithms for the iterative solution of a large-scale linear-quadratic optimal control problem involving a parabolic partial differential equation over a finite control horizon. We consider an “all at once” discretization of the problem and formulate three iterative algorithms. The first algorithm is based on preconditioning a symmetric positive definite reduced linear system involving only the unknown control variables; however inner-outer iterations are required. The second algorithm modifies the first algorithm to avoid inner-outer iterations by introducing an auxiliary variable. It yields a symmetric indefinite system with a positive definite block preconditioner. The third algorithm is the central focus of this paper. It modifies the preconditioner in the second algorithm by a parallel-in-time preconditioner based on the parareal algorithm. Theoretical results show that the preconditioned algorithms have optimal convergence properties and parallel scalability. Numerical experiments confirm the theoretical results.

We present a new iterative scheme for PageRank computation. The algorithm is applied to the linear system formulation of the problem, using inner-outer stationary iterations. It is simple, can be easily implemented and parallelized, and requires minimal storage overhead. Our convergence analysis shows that the algorithm is effective for a crude inner tolerance and is not sensitive to the choice of the parameters involved. The same idea can be used as a preconditioning technique for nonstationary schemes. Numerical examples featuring matrices of dimensions exceeding 100,000,000 in sequential and parallel environments demonstrate the merits of our technique. Our code is available online for viewing and testing, along with several large scale examples.

An algebraic multigrid (AMG) method is presented for the calculation of the stationary probability vector of an irreducible Markov chain. The method is based on standard AMG for nonsingular linear systems, but in a multiplicative, adaptive setting. A modified AMG interpolation formula is proposed that produces a nonnegative interpolation operator with unit row sums. We show how the adoption of a previously described lumping technique maintains the irreducible singular M-matrix character of the coarse-level operators on all levels. Together, these properties are sufficient to guarantee the well-posedness of the algorithm. Numerical results show how it leads to nearly optimal multigrid efficiency for a representative set of test problems.

We describe how a time-discretized wave equation in a homogeneous medium can be solved by boundary integral methods. The time discretization can be a multistep, Runge-Kutta, or a more general multistep-multistage method. The resulting convolutional system of boundary integral equations belongs to the family of convolution quadratures of Lubich. The aim of this work is twofold. It describes an efficient, robust, and easily parallelizable method for solving the semidiscretized system. The resulting algorithm has the main advantages of time-stepping methods and of Fourier synthesis: at each time-step, a system of linear equations with the same system matrix needs to be solved, yet computations can easily be done in parallel; the computational cost is almost linear in the number of time-steps; and only the Laplace transform of the time-domain fundamental solution is needed. The new aspect of the algorithm is that all this is possible without ever explicitly constructing the weights of the convolution quadrature. This approach also readily allows the use of modern data-sparse techniques to efficiently perform computation in space. We investigate theoretically and numerically to which extent hierarchical matrix ($\mathcal H$-matrix) techniques can be used to speed up the space computation. The second aim of this article is to perform a series of large-scale 3D experiments with a range of multistep and multistage time discretization methods: the backward difference formula of order 2 (BDF2), the Trapezoid rule, and the 3-stage Radau IIA methods are investigated in detail. One of the conclusions of the experiments is that the Radau IIA method often performs overwhelmingly better than the linear multistep methods, especially for problems with many reflections, yet, in connection with hyperbolic problems, BDFs have so far been predominant in the literature on convolution quadrature.

We develop and compare multilevel algorithms for solving constrained nonlinear variational problems via interior point methods. Several equivalent formulations of the linear systems arising at each iteration of the interior point method are compared from the point of view of conditioning and iterative solution. Furthermore, we show how a multilevel continuation strategy can be used to obtain good initial guesses ("hot starts") for each nonlinear iteration. Some minimal surface and volume-constrained image registration problems are used to illustrate the various approaches. © 2009 Society for Industrial and Applied Mathematics.

We present a line-search algorithm for large-scale continuous optimization. The algorithm is matrix-free in that it does not require the factorization of derivative matrices. Instead, it uses iterative linear system solvers. Inexact step computations are supported in order to save computational expense during each iteration. The algorithm is an interior-point approach derived from an inexact Newton method for equality constrained optimization proposed by Curtis, Nocedal, and Wachter [SIAM J. Optim., 20 (2009), pp. 1224-1249], with additional functionality for handling inequality constraints. The algorithm is shown to be globally convergent under loose assumptions. Numerical results are presented for nonlinear optimization test set collections and a pair of PDE-constrained model problems.

The innermost computational kernel of many large-scale scientific applications is often a large set of linear equations of the form $Ax=b$ which typically consumes a significant portion of the overall computational time required by the simulation. The traditional approach for solving this problem is to use direct methods. This approach is often preferred in industry because direct solvers are robust and effective for moderate size problems. However, direct methods can consume a huge amount of memory, and CPU time, in large-scale cases. In these cases, iterative techniques are the only viable alternative. Unfortunately, iterative methods lack the robustness of direct methods. The situation is especially difficult when the matrix is nonsymmetric. A lot of research has been devoted to trying to develop a robust iterative algorithm for nonsymmetric systems. The present paper describes a new robust and efficient algorithm aimed at solving iteratively nonsymmetric linear systems. It is based on looking for an approximation to the “optimal” polynomial $P_m(z)$ which satisfies $||P_m(z)||_{\infty}=\min_{Q\in\Pi_m}||Q(z)||_{\infty}$, $z\in D$, where $\Pi_m$ is the set of all polynomials of degree $m$ which satisfies $Q_m(0)=1$ and $D$ is a domain in the complex plane which includes all the eigenvalues of $A$. The resulting algorithm is an efficient one, especially in the case where we have a set of linear systems which share the same matrix $A$. We present several applications, including the exterior Helmholtz problem, which leads to a large indefinite, nonsymmetric, and complex system.

For the finite volume discretization of a second-order elliptic model problem, we derive a posteriori error estimates which take into account an inexact solution of the associated linear algebraic system. We show that the algebraic error can be bounded by constructing an equilibrated Raviart-Thomas-Nedelec discrete vector field whose divergence is given by a proper weighting of the residual vector. Next, claiming that the discretization error and the algebraic one should be in balance, we construct stopping criteria for iterative algebraic solvers. An attention is paid, in particular, to the conjugate gradient method which minimizes the energy norm of the algebraic error. Using this convenient balance, we also prove the efficiency of our a posteriori estimates; i.e., we show that they also represent a lower bound, up to a generic constant, for the overall energy error. A local version of this result is also stated. This makes our approach suitable for adaptive mesh refinement which also takes into account the algebraic error. Numerical experiments illustrate the proposed estimates and construction of efficient stopping criteria for algebraic iterative solvers.

We develop in this paper some efficient algorithms which are essential to implementations of spectral methods on the sparse grid by Smolyak's construction based on a nested quadrature. More precisely, we develop a fast algorithm for the discrete transform between the values at the sparse grid and the coefficients of expansion in a hierarchical basis; and by using the aforementioned fast transform, we construct two very efficient sparse spectral-Galerkin methods for a model elliptic equation. In particular, the Chebyshev-Legendre-Galerkin method leads to a sparse matrix with a much lower number of nonzero elements than that of low-order sparse grid methods based on finite elements or wavelets, and can be efficiently solved by a suitable sparse solver. Ample numerical results are presented to demonstrate the efficiency and accuracy of our algorithms.

Three techniques for improving the robustness and performance of incomplete factorization preconditioners for sparse systems with symmetric positive definite or mildly indefinite coefficient matrices are introduced. The primary contribution is two new block algorithms for incomplete factorization that result in an improvement in the performance of both the preconditioner generation and the iterative solution phases. One of the algorithms applies to matrices that have a natural block structure in their original form, and the other applies to matrices without natural dense blocks. Additionally, two relatively simple but highly effective heuristic strategies are introduced. These include selecting the solver based on the definiteness properties of the preconditioner and automatic selection and tuning of incomplete factorization parameters. All three techniques have adaptive components; i.e., the preconditioner-solver combination chooses parameters or algorithmic components based on the properties of the coefficient matrix and its incomplete factors. Two of the three techniques are applicable to incomplete LU factorization for unsymmetric systems as well.