High-Performance Algorithms for Numerical Linear Algebra

Eigenstructure of order-one-quasiseparable matrices. Three-term and two-term recurrence relations

This paper derives tradeoffs between three basic costs of a parallel algorithm: synchronization, data movement, and computational cost. These tradeoffs are lower bounds on the execution time of the algorithm which are independent of the number of processors, but dependent on the problem size. Therefore, they provide lower bounds on the parallel execution time of any algorithm computed by a system composed of any number of homogeneous components, each with associated computational, communication, and synchronization payloads. We employ a theoretical model counts the amount of work and data movement as a maximum of any execution path during the parallel computation. By considering this metric, rather than the total communication volume over the whole machine, we obtain new insights into the characteristics of parallel schedules for algorithms with non-trivial dependency structures. We also present reductions from BSP and LogP algorithms to our execution model, extending our lower bounds to these two models of parallel computation. We first develop our results for general dependency graphs and hypergraphs based on their expansion properties, then we apply the theorem to a number of specific algorithms in numerical linear algebra, namely triangular substitution, Gaussian elimination, and Krylov subspace methods. Our lower bound for LU factorization demonstrates the optimality of Tiskin's LU algorithm answering an open question posed in his paper, as well as of the 2.5D LU algorithm which has analogous costs. We treat the computations in a general manner by noting that the computations share a similar dependency hypergraph structure and analyzing the communication requirements of lattice hypergraph structures.

Hardware trends have motivated the development of mixed precision algorithms in numerical linear algebra, which aim to decrease runtime while maintaining acceptable accuracy. One recent development is the development of an adaptive precision sparse matrix–vector produce routine, which may be used to accelerate the solution of sparse linear systems by iterative methods. This approach is also applicable to the application of inexact preconditioners, such as sparse approximate inverse preconditioners used in Krylov subspace methods. In this work, we develop an adaptive precision sparse approximate inverse preconditioner and demonstrate its use within a five-precision GMRES-based iterative refinement method. We call this algorithm variant BSPAI-GMRES-IR. We then analyze the conditions for the convergence of BSPAI-GMRES-IR, and determine settings under which BSPAI-GMRES-IR will produce similar backward and forward errors as the existing SPAI-GMRES-IR method, the latter of which does not use adaptive precision in preconditioning. Our numerical experiments show that this approach can potentially lead to a reduction in the cost of storing and applying sparse approximate inverse preconditioners, although a significant reduction in cost may comes at the expense of increasing the number of GMRES iterations required for convergence.

We investigate randomized distributed algorithms for matrix computations over loosely coupled distributed systems, such as P2P networks or sensor networks. In this poster, we discuss orthogonalization methods and orthogonal iteration. These algorithms are very well understood in the sequential or in the classical parallel context, and they are important building blocks for many algorithms in numerical linear algebra.

Singular value decomposition (SVD) is a key step in many algorithms in statistics, machine learning and numerical linear algebra. While classical singular value decomposition has been made efficient in terms of computational complexity, classical algorithms are not able to fully utilise modern computing environments. The goal of this work is to survey various implementations and applications of randomized algorithms for SVD. Algorithms are compared in terms of accuracy and time of execution. On example of robust principal component analysis (RPCA), it is shown that using randomized algorithms can yield a significant speedup for image processing and similar applications.

This survey highlights the recent advances in algorithms for numerical linear algebra that have come from the technique of linear sketching, whereby given a matrix, one first compresses it to a much smaller matrix by multiplying it by a (usually) random matrix with certain properties. Much of the expensive computation can then be performed on the smaller matrix, thereby accelerating the solution for the original problem. In this survey we consider least squares as well as robust regression problems, low rank approximation, and graph sparsification. We also discuss a number of variants of these problems. Finally, we discuss the limitations of sketching methods.

Classical probabilistic rounding error analysis is particularly well suited to stochastic rounding (SR), and it yields strong results when dealing with floating-point algorithms that rely heavily on summation. For many numerical linear algebra algorithms, one can prove probabilistic error bounds that grow as O(nu), where n is the problem size and u is the unit roundoff. These probabilistic bounds are asymptotically tighter than the worst-case ones, which grow as O(nu). For certain classes of algorithms, SR has been shown to be unbiased. However, all these results were derived under the assumption that SR is implemented exactly, which typically requires a number of random bits that is too large to be suitable for practical implementations. We investigate the effect of the number of random bits on the probabilistic rounding error analysis of SR. To this end, we introduce a new rounding mode, limited-precision SR. By taking into account the number r of random bits used, this new rounding mode matches hardware implementations accurately, unlike the ideal SR operator generally used in the literature. We show that this new rounding mode is biased and that the bias is a function of r. As r approaches infinity, however, the bias disappears, and limited-precision SR converges to the ideal, unbiased SR operator. We develop a novel model for probabilistic error analysis of algorithms employing SR. Several numerical examples corroborate our theoretical findings.

Enhancements to array processors in the form of broadcast buses have been proposed for improving the speeds of algorithms in linear algebra, image processing and computational geometry. In this paper, we consider certain practical issues in such arrays which include fewer processors connected to the broadcast buses and finite time, namely log N time for broadcasting data to N processors. We propose a modified broadcast bus VLSI architecture which consists of a √N × √N mesh connected structure and a hierarchy of broadcast buses in each row and in each column such that every bus has k PE's. That is, in any row or column, in the first level there are \(\frac{{\sqrt N }}{k}\)buses with groups of k PE's connected to each bus. One PE from each of this group is connected to second level of buses in a similar manner. This is recursively done until there are only k PE's left which are connected by a broadcast bus. With this architecture, maximum, minimum, or sum of N values can be found in O(log N), median row of a binary picture in O(log N), shortest distance between two points in O(log N). This architecture is well suited for parallel processing of applications in Linear Algebra, Image Processing, Computational Geometry and Numerical Computations. These restricted connections to buses reduces the I/O ports significantly, and therefore, is more suitable for VLSI implementation.

The advent of the electronic computer has brought the problem of the solution of simultaneous linear equations to the fore, particularly in regard to the solution of symmetrical equations such as normal equations in least squares solutions. The object of this paper will be to remind readers of some of the usual ways of solving simultaneous equations, and then to discuss in greater detail the solution of normal equations by the Gauss–Doolittle, Cholesky and other methods with modifications and explanations.

The efficient solution of large linear matrix problems plays a central role in both linear and nonlinear structural analysis. Accordingly, a substantial effort has been allocated for the design of computer software in order to handle standard tasks like the solution of linear equations or eigenreduction.KeywordsFinite Element TechniqueJacobi MethodOperation CountInverse IterationHessenberg MatrixThese keywords were added by machine and not by the authors. This process is experimental and the keywords may be updated as the learning algorithm improves.

Modular Algorithms for Scientific Computing in Julia

SummaryProper range and precision analysis play an important role in the development of fixed‐point algorithms for embedded system applications. Numerical linear algebra algorithms used to find singular value decomposition of symmetric matrices are suitable for signal and image‐processing applications. These algorithms have not been attempted much in fixed‐point arithmetic. The reason is wide dynamic range of data and vulnerability of the algorithms to round‐off errors. For any real‐time application, the range of the input matrix may change frequently. This poses difficulty for constant and variable fixed‐point formats to decide on integer wordlengths during float‐to‐fixed conversion process because these formats involve determination of integer wordlengths before the compilation of the program. Thus, these formats may not guarantee to avoid overflow for all ranges of input matrices. To circumvent this problem, a novel dynamic fixed‐point format has been proposed to compute integer wordlengths adaptively during runtime. Lanczos algorithm with partial orthogonalization, which is a tridiagonalization step in computation of singular value decomposition of symmetric matrices, has been taken up as a case study. The fixed‐point Lanczos algorithm is tested for matrices with different dimensions and condition numbers along with image covariance matrix. The accuracy of fixed‐point Lanczos algorithm in three different formats has been compared on the basis of signal‐to‐quantization‐noise‐ratio, number of accurate fractional bits, orthogonality and factorization errors. Results show that dynamic fixed‐point format either outperforms or performs on par with constant and variable formats. Determination of fractional wordlengths requires minimization of hardware cost subject to accuracy constraint. In this context, we propose an analytical framework for deriving mean‐square‐error or quantization noise power among Lanczos vectors, which can serve as an accuracy constraint for wordlength optimization. Error is found to propagate through different arithmetic operations and finally accumulate in the last Lanczos vector. It is observed that variable and dynamic fixed‐point formats produce vectors with lesser round‐off error than constant format. All the three fixed‐point formats of Lanczos algorithm have been synthesized on Virtex 7 field‐programmable gate array using Vivado high‐level synthesis design tool. A comparative study of resource usage and power consumption is carried out. Copyright © 2015 John Wiley & Sons, Ltd.

This paper is to introduce a type of full multigrid method for the nonlinear eigenvalue problem. The main idea is to transform the solution of nonlinear eigenvalue problem into a series of solutions of the corresponding linear boundary value problems on the sequence of finite element spaces and nonlinear eigenvalue problems on the coarsest finite element space. The linearized boundary value problems are solved by some multigrid iterations. Besides the multigrid iteration, all other efficient iteration methods for solving boundary value problems can serve as the linear problem solver. We will prove that the computational work of this new scheme is truly optimal, the same as solving the linear corresponding boundary value problem. In this case, this type of iteration scheme certainly improves the overfull efficiency of solving nonlinear eigenvalue problems. Some numerical experiments are presented to validate the efficiency of the new method.

The first Survey and Review article in this issue is “Decay Properties of Spectral Projectors with Applications to Electronic Structure,” by Michele Benzi, Paolo Boito, and Nader Razouk. The linear scaling methods that motivate this work use a reformulation of the conventional electronic structure calculation. Instead of iterating through linear eigenvalue problems, a sequence of spectral projectors, known as density matrices, is computed, from which physical quantities of interest follow directly. Large-scale computations then become feasible, provided that entries in the density matrices are localized; that is, they decay rapidly as we move away from the diagonal. The authors first introduce the basic principles of electronic structure theory and then survey current computational approaches, most of which have arisen in the physical sciences, and their underlying localization requirements. The overview highlights a distinction, often blurred outside the mathematics literature, between physical insight/intuition and mathematical rigor. The key contribution is then to formulate the localization problem mathematically, with transparent assumptions, and develop a unified analytical approach. The tools of matrix analysis and approximation theory are used to justify, where feasible, the localization “results” that are relied upon by today's algorithms. The theory is framed in an asymptotic regime where the system volume increases but the density of particles remains fixed. This type of thermodynamic limit will be familiar to readers who work on multiscale models in physics and chemistry, but may seem strange to numerical analysts who are used to “convergence” in the sense of fixed problem size and arbitrarily fine mesh. The upshot is a sequence of matrices of increasing dimension that look very different from those arising in the discretization of PDEs. After pushing the theory as far as possible, the authors comment on the practical implications of their bounds and raise a number of open questions. They also point out that these decay results may find useful application in other, unrelated, areas: quantum information theory, complex networks, and eigenvalue solvers for tridiagonal matrices. This article illustrates how applied analysis can justify some of the leaps of faith made in the physical sciences and also resolve controversies (section 8.7). It will be of particular interest to readers who work in matrix computation and approximation theory. The second article, “Conditional Gradient Algorithms for Rank-One Matrix Approximations with a Sparsity Constraint,” by Ronny Luss and Marc Teboulle, shares with the first article the theme of studying algorithms based on large-scale matrix computations for which our current analytical understanding fails to explain the performance seen in practice. Given a symmetric matrix, $A$, most applied mathematicians know that dimension reduction can be achieved by using the dominant eigenvectors (singular vectors) of $A$ to form optimal low-rank approximations. In the statistics literature, modulo a centering operation, the same idea is called principal component analysis (PCA). This type of least-squares approach tries to recover all rows and columns of $A$; equivalently, if we think of $A$ as representing pairwise interactions between nodes in a network, it aims to summarize the full set of interactions. It is, of course, possible to impose a postprocessing threshold to cut off the less important contributions, but it is also reasonable to shoot directly for the important rows/columns or the important network nodes. This alternative leads us to the class of problems discussed here---matrix optimization problems where there is a prescribed upper bound on the number of nonzeros in the required “eigenvectors.” The authors discuss a variety of computational approaches that have been proposed in the literature and introduce a unifying framework, which they call ConGradU, to characterize and analyze them. The basic iterative scheme, summarized in Algorithm 1, has an optimization substep with closed-form solutions in a number of important cases. The authors also include an informative computational example based on the textual content of State of the Union addresses from 1790--2011. For example, Table 6.2 shows the key word stems associated with the three principal factors for both thresholded PCA and the more direct sparsity-seeking alternative. It is clear that the sparse version is summarizing different information. This article will appeal to readers who are keen to keep up with developments in modern optimization and to learn about techniques that fall under the current “Big Data” banner.

SummaryIn this work, we consider the numerical solution of large nonlinear eigenvalue problems that arise in thermoacoustic simulations involved in the stability analysis of large combustion devices. We briefly introduce the physical modeling that leads to a nonlinear eigenvalue problem that is solved using a nonlinear fixed point iteration scheme. Each step of this nonlinear method requires the solution of a complex non‐Hermitian linear eigenvalue problem. We review a set of state‐of‐the‐art eigensolvers and discuss strategies to recycle spectral information from one nonlinear step to the next. More precisely, we consider the Jacobi–Davidson algorithm, the Implicitly Restarted Arnoldi method, the Krylov–Schur solver and its block‐variant, and the subspace iteration method with Chebyshev acceleration. On a small test example, we study the relevance of the different approaches and illustrate on a large industrial test case the performance of the parallel solvers best suited to recycle spectral information for large‐scale thermoacoustic stability analysis. Copyright © 2015 John Wiley & Sons, Ltd.