Search

The efficient utilization of mixed-precision numerical linear algebra algorithms can offer attractive acceleration to scientific computing applications. Especially with the hardware integration of low-precision special-function units designed for machine learning applications, the traditional numerical algorithms community urgently needs to reconsider the floating point formats used in the distinct operations to efficiently leverage the available compute power. In this work, we provide a comprehensive survey of mixed-precision numerical linear algebra routines, including the underlying concepts, theoretical background, and experimental results for both dense and sparse linear algebra problems.

We show that the H2-norm model reduction problem for single-input/single-output (SISO) linear time-invariant (LTI) systems is essentially an eigenvalue problem (EP), from which the globally optimal solution(s) can be retrieved. The first-order optimality conditions of this model reduction problem constitute a system of multivariate polynomial equations that can be converted to an affine (or inhomogeneous) multiparameter eigenvalue problem (AMEP). We solve this AMEP by using the so-called augmented block Macaulay matrix, which is introduced in this paper and has a special (block) multi-shift invariant null space. The set of all stationary points of the optimization problem, i.e., the (2r)-tuples (r is the order of the reduced model) of affine eigenvalues and eigenvectors of the AMEP, follows from a standard EP related to the structure of that null space. At least one of these (2r)-tuples corresponds to the globally optimal solution of the H2-norm model reduction problem. We present a simple numerical example to illustrate our approach.

This survey describes probabilistic algorithms for linear algebraic computations, such as factorizing matrices and solving linear systems. It focuses on techniques that have a proven track record for real-world problems. The paper treats both the theoretical foundations of the subject and practical computational issues.Topics include norm estimation, matrix approximation by sampling, structured and unstructured random embeddings, linear regression problems, low-rank approximation, subspace iteration and Krylov methods, error estimation and adaptivity, interpolatory and CUR factorizations, Nyström approximation of positive semidefinite matrices, single-view (‘streaming’) algorithms, full rank-revealing factorizations, solvers for linear systems, and approximation of kernel matrices that arise in machine learning and in scientific computing.

Randomized numerical linear algebra - RandNLA, for short - concerns the use of randomization as a resource to develop improved algorithms for large-scale linear algebra computations. The origins of contemporary RandNLA lay in theoretical computer science, where it blossomed from a simple idea: randomization provides an avenue for computing approximate solutions to linear algebra problems more efficiently than deterministic algorithms. This idea proved fruitful in the development of scalable algorithms for machine learning and statistical data analysis applications. However, RandNLA's true potential only came into focus upon integration with the fields of numerical analysis and "classical" numerical linear algebra. Through the efforts of many individuals, randomized algorithms have been developed that provide full control over the accuracy of their solutions and that can be every bit as reliable as algorithms that might be found in libraries such as LAPACK. Recent years have even seen the incorporation of certain RandNLA methods into MATLAB, the NAG Library, NVIDIA's cuSOLVER, and SciKit-Learn. For all its success, we believe that RandNLA has yet to realize its full potential. In particular, we believe the scientific community stands to benefit significantly from suitably defined "RandBLAS" and "RandLAPACK" libraries, to serve as standards conceptually analogous to BLAS and LAPACK. This 200-page monograph represents a step toward defining such standards. In it, we cover topics spanning basic sketching, least squares and optimization, low-rank approximation, full matrix decompositions, leverage score sampling, and sketching data with tensor product structures (among others). Much of the provided pseudo-code has been tested via publicly available MATLAB and Python implementations., Comment: v1: this is the first arXiv release of LAPACK Working Note 299. v2: complete rewrite of the subsection on trace estimation, among other changes. See frontmatter page ii (pdf page 5) for revision history

Randomized Numerical Linear Algebra (RandNLA) is a powerful class of methods, widely used in High Performance Computing (HPC). RandNLA provides approximate solutions to linear algebra functions applied to large signals, at reduced computational costs. However, the randomization step for dimensionality reduction may itself become the computational bottleneck on traditional hardware. Leveraging near constant-time linear random projections delivered by LightOn Optical Processing Units we show that randomization can be significantly accelerated, at negligible precision loss, in a wide range of important RandNLA algorithms, such as RandSVD or trace estimators., Comment: Add "This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 860830"

Understanding dynamics in complex systems is challenging because there are many degrees of freedom, and those that are most important for describing events of interest are often not obvious. The leading eigenfunctions of the transition operator are useful for visualization, and they can provide an efficient basis for computing statistics such as the likelihood and average time of events (predictions). Here we develop inexact iterative linear algebra methods for computing these eigenfunctions (spectral estimation) and making predictions from a data set of short trajectories sampled at finite intervals. We demonstrate the methods on a low-dimensional model that facilitates visualization and a high-dimensional model of a biomolecular system. Implications for the prediction problem in reinforcement learning are discussed., Comment: 24 pages, 16 figures

The solution of sparse triangular linear systems (sptrsv) is the bottleneck of many numerical methods. Thus, it is crucial to count with efficient implementations of such kernel, at least for commonly used platforms. In this sense, Field–Programmable Gate Arrays (FPGAs) have evolved greatly in the last years, entering the HPC hardware ecosystem largely due to their superior energy–efficiency relative to more established accelerators. Up until recently, the design for FPGAs implied the use of low–level Hardware Description Languages (HDL) such as VHDL or Verilog. Nowadays, manufacturers are making a large effort to adopt High–Level Synthesis languages like C/C++ or OpenCL, but the gap between their performance and that of HDLs is not yet fully studied. This work focuses on the performance offered by FPGAs to compute the sptrsv using OpenCL. For this purpose, we implement different parallel variants of this kernel and experimentally evaluate several setups, varying among others the work–group size, the number of compute units, the unroll–factor and the vectorization–factor.

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

Matrix computations lie at the heart of many scientific computations. While sophisticated algorithms have been established for various numerical linear algebra problems such as the solution of linear simultaneous equations and eigenvalue problems, they require considerable modification with the advent of exaFLOPS- scale supercomputers, which are expected to have a huge number of computing cores, deep memory hierarchy, and increased probability of hardware errors. In this chapter, we discuss projected hardware characteristics of exaFLOPS machines and summarize the challenges to be faced by numerical linear algebra algorithms in the near future. Based on these preparations, we present a brief survey of recent research efforts in the field of numerical linear algebra targeted at meeting these challenges.

This survey provides an introduction to the use of randomization in the design of fast algorithms for numerical linear algebra. These algorithms typically examine only a subset of the input to solve basic problems approximately, including matrix multiplication, regression and low-rank approximation. The survey describes the key ideas and gives complete proofs of the main results in the field. A central unifying idea is sampling the columns (or rows) of a matrix according to their squared lengths.

My dissertation project addresses the integration of numerical linear algebra approaches to the field of computer graphics and visualization, especially physically-based Fluid Animations. I mainly focus on (matrix) decomposition techniques with regard to spectral theory and related concepts. In many applications of Fluid Animations, in particular particle-based ones, these techniques are not consulted due to its irregular nature. However, the use of spectral-theoretic approaches leads to both theoretical insights and practical enhancements of the algorithm in terms of precision, efficiency, and stability. My research is conducted at the University of Stuttgart and at the Hochschule der Medien, Stuttgart, where I am a member of the joint graduate school Digital Media.

Large-scale problems have attracted much attention in the last decades since they arise from different applications in several fields. Moreover, the matrices that are involved in those problems are often sparse, this is, the majority of their entries are zero. Around 40 years ago, the most common problems related to large-scale and sparse matrices consisted in solving linear systems, finding eigenvalues and/or eigenvectors, solving least square problems or computing singular value decompositions. However, in the last years, large-scale and sparse problems of different natures have appeared, motivating and challenging numerical linear algebra to develop effective and efficient algorithms to solve them. Common difficulties that appear during the development of algorithms for solving modern large-scale problems are related to computational costs, storage issues and CPU time, given the large size of the matrices, which indicate that direct methods can not be used. This suggests that projection methods based on Krylov subspaces are a good option to develop procedures for solving large-scale and sparse modern problems. In this PhD Thesis we develop novel and original algorithms for solving two large-scale modern problems in numerical linear algebra: first, we introduce the R-CORK method for solving rational eigenvalue problems and, second, we present projection methods to compute the solution of T-Sylvester matrix equations, both based on Krylov subspaces. The R-CORK method is an extension of the compact rational Krylov method (CORK) [104] introduced to solve a family of nonlinear eigenvalue problems that can be expressed and linearized in certain particular ways and which include arbitrary polynomial eigenvalue problems, but not arbitrary rational eigenvalue problems. The R-CORK method exploits the structure of the linearized problem by representing the Krylov vectors in a compact form in order to reduce the cost of storage, resulting in a method with two levels of orthogonalization. The first level of orthogonalization works with vectors of the same size as the original problem, and the second level works with vectors of size much smaller than the original problem. Since vectors of the size of the linearization are never stored or orthogonalized, R-CORK is more efficient from the point of view of memory and orthogonalization costs than the classical rational Krylov method applied to the linearization. Moreover, since the R-CORK method is based on a classical rational Krylov method, the implementation of implicit restarting is possible and we present an efficient way to do it, that preserves the compact representation of the Krylov vectors. We also introduce in this dissertation projection methods for solving the TSylvester equation, which has recently attracted considerable attention as a consequence of its close relation to palindromic eigenvalue problems and other applications. The theory concerning T-Sylvester equations is rather well understood, and before the work in this thesis, there were stable and efficient numerical algorithms to solve these matrix equations for small- to medium- sized matrices. However, developing numerical algorithms for solving large-scale T-Sylvester equations was a completely open problem. In this thesis, we introduce several projection methods based on block Krylov subspaces and extended block Krylov subspaces for solving the T-Sylvester equation when the right-hand side is a low-rank matrix. We also offer an intuition on the expected convergence of the algorithm based on block Krylov subspaces and a clear guidance on which algorithm is the most convenient to use in each situation. All the algorithms presented in this thesis have been extensively tested, and the reported numerical results show that they perform satisfactorily in practice. Adicionalmente se recibió ayuda parcial de los proyectos de investigación: “Structured Numerical Linear Algebra: Matrix Polynomials, Special Matrices, and Conditioning” (Ministerio de Economía y Competitividad de España, Número de proyecto: MTM2012-32542) y “Structured Numerical Linear Algebra for Constant, Polynomial and Rational Matrices” (Ministerio de Economía y Competitividad de España, Número de proyecto: MTM2015-65798-P), donde el investigador principal de ambos proyectos fue Froilán Martínez Dopico. Programa Oficial de Doctorado en Ingeniería Matemática Presidente: José Mas Marí.- Secretario: Fernando de Terán Vergara.- Vocal: José Enrique Román Moltó

Over the last ten years, tremendous speedups for problems in randomized numerical linear algebra such as low rank approximation and regression have been made possible via the technique of randomized data dimensionality reduction, also known as sketching. In theory, such algorithms have led to optimal input sparsity time algorithms for a wide array of problems. While several scattered implementations of such methods exist, the goal of this work is to provide a comprehensive comparison of such methods to alternative approaches. We investigate least squares regression, iteratively reweighted least squares, logistic regression, robust regression with Huber and Bisquare loss functions, leverage score computation, Frobenius norm low rank approximation, and entrywise $\ell_1$-low rank approximation. We give various implementation techniques to speed up several of these algorithms, and the resulting implementations demonstrate the tradeoffs of such techniques in practice.

We propose and simulate integrated optical devices for accelerating numerical linear algebra (NLA) calculations. Data is modulated on chirped optical pulses and these propagate through a multimode waveguide where speckle provides the random projections needed for NLA dimensionality reduction.

A thorough survey of energy efficient techniques for linear algebra operations.A complete list of valuable research efforts.A deeply review of all aspects of feasible solutions.Detailed discussion and comprehensive comparison.Thought-provoking directions of future research are provided. Extreme scale supercomputers available before the end of this decade are expected to have 100 million to 1billion computing cores. The power and energy efficiency issue has become one of the primary concerns of extreme scale high performance scientific computing. This paper surveys the research on saving power and energy for numerical linear algebra algorithms in high performance scientific computing on supercomputers around the world. We first stress the significance of numerical linear algebra algorithms in high performance scientific computing nowadays, followed by a background introduction on widely used numerical linear algebra algorithms and software libraries and benchmarks. We summarize commonly deployed power management techniques for reducing power and energy consumption in high performance computing systems by presenting power and energy models and two fundamental types of power management techniques: static and dynamic. Further, we review the research on saving power and energy for high performance numerical linear algebra algorithms from four aspects: profiling, trading off performance, static saving, and dynamic saving, and summarize state-of-the-art techniques for achieving power and energy efficiency in each category individually. Finally, we discuss potential directions of future work and summarize the paper.

Generating sparsity patterns for effective block-Jacobi preconditioners is a challenging and computationally expensive task, in particular for problems with unknown origin. In this paper we design a convolutional neural network (CNN) to detect natural block structures in matrix sparsity patterns. For test matrices where a natural block structure is complemented with a random distribution of nonzeros (noise), we show that a trained network succeeds in identifying strongly connected components with more than 95% prediction accuracy, and the resulting block-Jacobi preconditioner effectively accelerating an iterative GMRES solver. Segmenting a matrix into diagonal tiles of size 128x128, for each tile the sparsity pattern of an effective block-Jacobi preconditioner can be generated in less than a millisecond when using a production-line GPU.

Summary Proper 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 chapter is based on lectures on Randomized Numerical Linear Algebra from the 2016 Park City Mathematics Institute summer school on The Mathematics of Data., Comment: To appear in the edited volume of lectures from the 2016 PCMI summer school

We establish a new perturbation theory for orthogonal polynomials using a Riemann--Hilbert approach and consider applications in numerical linear algebra and random matrix theory. This new approach shows that the orthogonal polynomials with respect to two measures can be effectively compared using the difference of their Stieltjes transforms on a suitably chosen contour. Moreover, when two measures are close and satisfy some regularity conditions, we use the theta functions of a hyperelliptic Riemann surface to derive explicit and accurate expansion formulae for the perturbed orthogonal polynomials. In contrast to other approaches, a key strength of the methodology is that estimates can remain valid as the degree of the polynomial grows. The results are applied to analyze several numerical algorithms from linear algebra, including the Lanczos tridiagonalization procedure, the Cholesky factorization and the conjugate gradient algorithm. As a case study, we investigate these algorithms applied to a general spiked sample covariance matrix model by considering the eigenvector empirical spectral distribution and its limits. For the first time, we give precise estimates on the output of the algorithms, applied to this wide class of random matrices, as the number of iterations diverges. In this setting, beyond the first order expansion, we also derive a new mesoscopic central limit theorem for the associated orthogonal polynomials and other quantities relevant to numerical algorithms.

In recent years, a variety of randomized constructions of sketching matrices have been devised, that have been used in fast algorithms for numerical linear algebra problems, such as least squares regression, low-rank approximation, and the approximation of leverage scores. A key property of sketching matrices is that of subspace embedding. In this paper, we study sketching matrices that are obtained from bipartite graphs that are sparse, i.e., have left degree s that is small. In particular, we explore two popular classes of sparse graphs, namely, expander graphs and magical graphs. For a given subspace ${\mathcal{U}} \subseteq {{\mathbb{R}}^n}$ of dimension k, we show that the magical graph with left degree s = 2 yields a (1 ± ϵ) l 2 -subspace embedding for ${\mathcal{U}}$, if the number of right vertices (the sketch size) $m = {\mathcal{O}}\left( {{k^2}/{\varepsilon ^2}} \right)$. The expander graph with $s = {\mathcal{O}}(\log k/\varepsilon )$ yields a subspace embedding for $m = {\mathcal{O}}\left( {k\log k/{\varepsilon ^2}} \right)$. We also discuss the construction of sparse sketching matrices with reduced randomness using expanders based on error-correcting codes. Empirical results on various synthetic and real datasets show that these sparse graph sketching matrices work very well in practice.

We create classical (non-quantum) dynamic data structures supporting queries for recommender systems and least-squares regression that are comparable to their quantum analogues. De-quantizing such algorithms has received a flurry of attention in recent years; we obtain sharper bounds for these problems. More significantly, we achieve these improvements by arguing that the previous quantum-inspired algorithms for these problems are doing leverage or ridge-leverage score sampling in disguise; these are powerful and standard techniques in randomized numerical linear algebra. With this recognition, we are able to employ the large body of work in numerical linear algebra to obtain algorithms for these problems that are simpler or faster (or both) than existing approaches. Our experiments demonstrate that the proposed data structures also work well on real-world datasets., Comment: Minor edits to exposition

The traditional metric for the efficiency of a numerical algorithm has been the number of arithmetic operations it performs. Technological trends have long been reducing the time to perform an arithmetic operation, so it is no longer the bottleneck in many algorithms; rather, communication, or moving data, is the bottleneck. This motivates us to seek algorithms that move as little data as possible, either between levels of a memory hierarchy or between parallel processors over a network. In this paper we summarize recent progress in three aspects of this problem. First we describe lower bounds on communication. Some of these generalize known lower bounds for dense classical (O(n3)) matrix multiplication to all direct methods of linear algebra, to sequential and parallel algorithms, and to dense and sparse matrices. We also present lower bounds for Strassen-like algorithms, and for iterative methods, in particular Krylov subspace methods applied to sparse matrices. Second, we compare these lower bounds to widely used versions of these algorithms, and note that these widely used algorithms usually communicate asymptotically more than is necessary. Third, we identify or invent new algorithms for most linear algebra problems that do attain these lower bounds, and demonstrate large speed-ups in theory and practice.

The statistical analysis of Randomized Numerical Linear Algebra (RandNLA) algorithms within the past few years has mostly focused on their performance as point estimators. However, this is insufficient for conducting statistical inference, e.g., constructing confidence intervals and hypothesis testing, since the distribution of the estimator is lacking. In this article, we develop an asymptotic analysis to derive the distribution of RandNLA sampling estimators for the least-squares problem. In particular, we derive the asymptotic distribution of a general sampling estimator with arbitrary sampling probabilities. The analysis is conducted in two complementary settings, i.e., when the objective of interest is to approximate the full sample estimator or is to infer the underlying ground truth model parameters. For each setting, we show that the sampling estimator is asymptotically normally distributed under mild regularity conditions. Moreover, the sampling estimator is asymptotically unbiased in both settings. Based on our asymptotic analysis, we use two criteria, the Asymptotic Mean Squared Error (AMSE) and the Expected Asymptotic Mean Squared Error (EAMSE), to identify optimal sampling probabilities. Several of these optimal sampling probability distributions are new to the literature, e.g., the root leverage sampling estimator and the predictor length sampling estimator. Our theoretical results clarify the role of leverage in the sampling process, and our empirical results demonstrate improvements over existing methods., Comment: 33 pages, 13 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.

SUMMARY The mixed sensitivity minimization problem is revisited for a class of single-input-single-output unstable infinite dimensional plants with low order weights. It is shown that H ∞ controllers can be computed from the singularity conditions of a parameterized matrix whose dimension is the same as the order of the sensitivity weight. The result is applied to the design of H ∞ controllers with integral action. Connections with the so-called Hamiltonian approach are also established. Copyright © 2012 John Wiley & Sons, Ltd.

This survey describes probabilistic algorithms for linear algebra computations, such as factorizing matrices and solving linear systems. It focuses on techniques that have a proven track record for real-world problem instances. The paper treats both the theoretical foundations of the subject and the practical computational issues. Topics covered include norm estimation; matrix approximation by sampling; structured and unstructured random embeddings; linear regression problems; low-rank approximation; subspace iteration and Krylov methods; error estimation and adaptivity; interpolatory and CUR factorizations; Nystr��m approximation of positive-semidefinite matrices; single view ("streaming") algorithms; full rank-revealing factorizations; solvers for linear systems; and approximation of kernel matrices that arise in machine learning and in scientific computing.

This work explores how randomization can be exploited to deliver sophisticatedalgorithms with provable bounds for: (i) The approximation of matrix functions, suchas the log-determinant and the Von-Neumann entropy; and (ii) The low-rank approximationof matrices. Our algorithms are inspired by recent advances in RandomizedNumerical Linear Algebra (RandNLA), an interdisciplinary research area that exploitsrandomization as a computational resource to develop improved algorithms forlarge-scale linear algebra problems. The main goal of this work is to encourage thepractical use of RandNLA approaches to solve Big Data bottlenecks at industriallevel. Our extensive evaluation tests are complemented by a thorough theoreticalanalysis that proves the accuracy of the proposed algorithms and highlights theirscalability as the volume of data increases. Finally, the low computational time andmemory consumption, combined with simple implementation schemes that can easilybe extended in parallel and distributed environments, render our algorithms suitablefor use in the development of highly efficient real-world software.