Search

QR decomposition is one of the most useful factorization kernels in modern numerical linear algebra algorithms. In particular, the decomposition of tall-and-skinny matrices (TSMs) has major applications in areas including scientific computing, machine learning, image processing, wireless networks, and numerical methods. Traditionally, CPUs and GPUs have achieved better throughput on these applications by using large cache hierarchies and compute cores running at a high frequency, leading to high power consumption. With the advent of heterogeneous platforms, however, FPGAs are emerging as a promising viable alternative. In this work, we propose a high-throughput FPGA-based engine that has a very high computational efficiency (ratio of achieved to peak throughput) compared to similar QR solvers running on FPGAs. Although comparable QR solvers achieve an efficiency of 36%, our design exhibits an efficiency of 54%. For TSMs, our experimental results show that our design can outperform highly optimized QR solvers running on CPUs and GPUs. For TSMs with more than 50K rows, our design outperforms the Intel MKL solver running on an Intel quad-core processor by a factor of 1.5×. For TSMs containing 256 columns or less, our design outperforms the NVIDIA CUBLAS solver running on a K40 GPU by a factor of 3.0×. In addition to being fast, our design is energy efficient—competing platforms execute up to 0.6 GFLOPS/Joule, whereas our design executes more than 1.0 GFLOPS/Joule.

Rational Krylov subspaces have become a fundamental ingredient in numerical linear algebra methods associated with reduction strategies. Nonetheless, many structural properties of the reduced matrices in these subspaces are not fully understood. We advance in this analysis by deriving bounds on the entries of rational Krylov reduced matrices and of their functions, that ensure an a-priori decay of their entries as we move away from the main diagonal. As opposed to other decay pattern results in the literature, these properties hold in spite of the lack of any banded structure in the considered matrices. Numerical experiments illustrate the quality of our results.

We present a perturbed subspace iteration algorithm to approximate the lowermost eigenvalue cluster of an elliptic eigenvalue problem. As a prototype, we consider the Laplace eigenvalue problem posed in a polygonal domain. The algorithm is motivated by the analysis of inexact (perturbed) inverse iteration algorithms in numerical linear algebra. We couple the perturbed inverse iteration approach with mesh refinement strategy based on residual estimators. We demonstrate our approach on model problems in two and three dimensions., Comment: 30 pages, 15 figures

In this expository paper we illustrate the role that the field of values (or numerical range) of a matrix plays in connection with certain problems of numerical analysis. These include the approximation of matrix functions and the convergence of preconditioned Krylov subspace methods for solving large systems of equations arising from the discretization of partial differential equations.

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.

This work develops validated numerical methods for linear stability analysis at an equilibrium solution of a system of delay differential equations (DDEs). In addition to providing mathematically rigorous bounds on the locations of eigenvalues, our method leads to validated counts. For example we obtain the computer assisted theorems about Morse indices (number of unstable eigenvalues). The case of a single constant delay is considered. The method downplays the role of the scalar transcendental characteristic equation in favor of a functional analytic approach exploiting the strengths of numerical linear algebra/techniques of scientific computing. The idea is to consider an equivalent implicitly defined discrete time dynamical system which is projected onto a countable basis of Chebyshev series coefficients. The projected problem reduces to questions about certain sparse infinite matrices, which are well approximated by \begin{document}$ N \times N $\end{document} matrices for large enough \begin{document}$ N $\end{document} . We develop the appropriate truncation error bounds for the infinite matrices, provide a general numerical implementation which works for any system with one delay, and discuss computer-assisted theorems in a number of example problems.

This thesis concerns the optimization and application of low-rank methods, with a special focus on tensor trains (TTs). In particular, we develop methods for computing TT approximations of a given tensor in a variety of low-rank formats and we show how to solve the tensor completion problem for TTs using Riemannian methods. This is then applied to train a machine learning (ML) estimator based on discretized functions. We also study randomized methods for obtaining low-rank approximations of matrices and tensors. Finally, we consider how such randomized methods can be used to solve general linear matrix and tensor equations.

These lecture notes accompany a refresher course in applied mathematics with a focus on numerical concepts (Part I), numerical linear algebra (Part II), numerical analysis, Fourier series and Fourier transforms (Part III), and differential equations (Part IV). Several numerical projects for group work are provided in Part V. In these projects, the tasks are threefold: mathematical modeling, algorithmic design, and implementation. Therein, it is important to draw interpretations of the obtained results and provide measures (Parts I-IV) how to build confidence into numerical findings such intuition, error analysis, convergence analysis, and comparison to manufactured solutions. Both authors have been jointly teaching over several years this class and bring in a unique mixture of their respective teaching and research fields.

Numerical analysis of truncated periodic structures by utilizing conventional method of moments (MoM) technique is challenging because the number of unknowns involved is often very large; hence, the computational cost is high. This letter presents a novel approach, based on the dipole moment (DM) method, for analyzing electromagnetic scattering from truncated periodic structures whose geometrical and/or material properties vary locally from those of a strictly doubly infinite periodic structure. The main advantage of using the DM method is that it provides convenient closed-form expressions for the scattered fields, facilitating the matrix computation in the context of the MoM much more efficiently than when the conventional Rao–Wilton–Glisson (RWG) type of basis functions are used for the same computation. Numerical results obtained by using the proposed approach are compared with those derived from a commercial MoM software package, and good agreement is found. The time advantage is clearly evident from this comparison.

The story of linear programming is one with all the elements of a grand historical drama. The original idea of testing if a polyhedron is non-empty by using a variable elimination to project down one dimension at a time until a tautology emerges dates back to a paper by Fourier in 1823. This gets re-invented in the 1930s by Motzkin. The real interest in linear programming happens during World War II when mathematicians ponder best ways of utilising resources at a time when they are constrained. The problem of optimising a linear function over a set of linear inequalities becomes the focus of the effort. Dantzig’s Simplex Method is announced and the Rand Corporation becomes a hot bed of computational mathematics. The range of applications of this modelling approach grows and the powerful machinery of numerical analysis and numerical linear algebra becomes a major driver for the advancement of computing machines. In the 1970s, constructs of theoretical computer science indicate that linear programming may in fact define the frontier of tractable problems that can be solved effectively on large instances. This raised a series of questions and answers: Is the Simplex Method a polynomial-time method and if not can we construct novel polynomial time methods, etc. And that is how the Ellipsoid Method from the Soviet Union and the Interior Point Method from Bell Labs make their way into this story as the heroics of Khachiyan and Karmarkar. We have called this paper a primer on linear programming since it only gives the reader a quick narrative of the grand historical drama. Hopefully it motivates a young reader to delve deeper and add another chapter.

In this thesis, we use statistical inference and competitive games to design algorithms for computational mathematics. In the first part, comprising chapters two through six, we use ideas from Gaussian process statistics to obtain fast solvers for differential and integral equations. We begin by observing the equivalence of conditional (near-)independence of Gaussian processes and the (near-)sparsity of the Cholesky factors of its precision and covariance matrices. This implies the existence of a large class of dense matrices with almost sparse Cholesky factors, thereby greatly increasing the scope of application of sparse Cholesky factorization. Using an elimination ordering and sparsity pattern motivated by the screening effect in spatial statistics, we can compute approximate Cholesky factors of the covariance matrices of Gaussian processes admitting a screening effect in near-linear computational complexity. These include many popular smoothness priors such as the Matérn class of covariance functions. In the special case of Green's matrices of elliptic boundary value problems (with possibly unknown elliptic operators of arbitrarily high order, with possibly rough coefficients), we can use tools from numerical homogenization to prove the exponential accuracy of our method. This result improves the state-of-the-art for solving general elliptic integral equations and provides the first proof of an exponential screening effect. We also derive a fast solver for elliptic partial differential equations, with accuracy-vs-complexity guarantees that improve upon the state-of-the-art. Furthermore, the resulting solver is performant in practice, frequently beating established algebraic multigrid libraries such as AMGCL and Trilinos on a series of challenging problems in two and three dimensions. Finally, for any given covariance matrix, we obtain a closed-form expression for its optimal (in terms of Kullback-Leibler divergence) approximate inverse-Cholesky factorization subject to a sparsity constraint, recovering Vecchia approximation and factorized sparse approximate inverses. Our method is highly robust, embarrassingly parallel, and further improves our asymptotic results on the solution of elliptic integral equations. We also provide a way to apply our techniques to sums of independent Gaussian processes, resolving a major limitation of existing methods based on the screening effect. As a result, we obtain fast algorithms for large-scale Gaussian process regression problems with possibly noisy measurements. In the second part of this thesis, comprising chapters seven through nine, we study continuous optimization through the lens of competitive games. In particular, we consider competitive optimization, where multiple agents attempt to minimize conflicting objectives. In the single-agent case, the updates of gradient descent are minimizers of quadratically regularized linearizations of the loss function. We propose to generalize this idea by using the Nash equilibria of quadratically regularized linearizations of the competitive game as updates (linearize the game). We provide fundamental reasons why the natural notion of linearization for competitive optimization problems is given by the multilinear (as opposed to linear) approximation of the agents' loss functions. The resulting algorithm, which we call competitive gradient descent, thus provides a natural generalization of gradient descent to competitive optimization. By using ideas from information geometry, we extend CGD to competitive mirror descent (CMD) that can be applied to a vast range of constrained competitive optimization problems. CGD and CMD resolve the cycling problem of simultaneous gradient descent and show promising results on problems arising in constrained optimization, robust control theory, and generative adversarial networks. Finally, we point out the GAN-dilemma that refutes the common interpretation of GANs as approximate minimizers of a divergence obtained in the limit of a fully trained discriminator. Instead, we argue that GAN performance relies on the implicit competitive regularization (ICR) due to the simultaneous optimization of generator and discriminator and support this hypothesis with results on low-dimensional model problems and GANs on CIFAR10.

The notion of a tensor captures three great ideas: equivariance, multilinearity, separability. But trying to be three things at once makes the notion difficult to understand. We will explain tensors in an accessible and elementary way through the lens of linear algebra and numerical linear algebra, elucidated with examples from computational and applied mathematics., Comment: 208 pages, 5 figures

Under the direction of the Institute for Informatics and Telematics in Pisa, Calcolo publishes original contributions on numerical analysis and its applications, and on the theory of computation. The main focus of the journal is on numerical linear algebra, approximation theory and its applications, numerical solutions of differential and integral equations, computational complexity, algorithmics, mathematical aspects of computer science, and optimization theory. Calcolo contains expository papers which introduce emerging topics as well as a abstracts of PhD theses, news and reports from conferences, and book reviews

Purpose This paper aims to propose a novel direct method for indefinite algebraic linear systems. It is well adapted for sparse linear systems, such as those of two-dimensional (2-D) finite elements problems, especially for coupled systems. Design/methodology/approach The proposed method is developed on an example of an indefinite symmetric matrix. The algorithm of the method is given next, and a comparison between the numbers of operations required by the method and the Cholesky method is also given. Finally, an application on a magnetostatic problem for classical methods (Gauss and Cholesky) shows the relative efficiency of the proposed method. Findings The proposed method can be used advantageously for 2-D finite elements in stepping methods without using a block decomposition of matrices. Research limitations/implications This method is advantageous for direct linear solving for 2-D problems, but it is not recommended at this time for three-dimensional problems. Originality/value The proposed method is the first direct solver for algebraic linear systems proposed since more than a half century. It is not limited for symmetric positive systems such as many of direct and iterative methods.

Funder: University of Cambridge, In this paper we study a class of fast geometric image inpainting methods based on the idea of filling the inpainting domain in successive shells from its boundary inwards. Image pixels are filled by assigning them a color equal to a weighted average of their already filled neighbors. However, there is flexibility in terms of the order in which pixels are filled, the weights used for averaging, and the neighborhood that is averaged over. Varying these degrees of freedom leads to different algorithms, and indeed the literature contains several methods falling into this general class. All of them are very fast, but at the same time all of them leave undesirable artifacts such as “kinking” (bending) or blurring of extrapolated isophotes. Our objective in this paper is to build a theoretical model in order to understand why these artifacts occur and what, if anything, can be done about them. Our model is based on two distinct limits: a continuum limit in which the pixel width h→0 and an asymptotic limit in which h>0 but h≪1. The former will allow us to explain “kinking” artifacts (and what to do about them) while the latter will allow us to understand blur. Both limits are derived based on a connection between the class of algorithms under consideration and stopped random walks. At the same time, we consider a semi-implicit extension in which pixels in a given shell are solved for simultaneously by solving a linear system. We prove (within the continuum limit) that this extension is able to completely eliminate kinking artifacts, which we also prove must always be present in the direct method. Finally, we show that although our results are derived in the context of inpainting, they are in fact abstract results that apply more generally. As an example, we show how our theory can also be applied to a problem in numerical linear algebra.

A determinantal point process is a stochastic point process that is commonly used to capture negative correlations. It has become increasingly popular in machine learning in recent years. Sampling a determinantal point process however remains a computationally intensive task. This note introduces a heuristic independent particle approximation to determinantal point processes. The approximation is based on the physical intuition of fermions and is implemented using standard numerical linear algebra routines. Sampling from this independent particle approximation can be performed at a negligible cost. Numerical results are provided to demonstrate the performance of the proposed algorithm.

Under the direction of the Institute for Informatics and Telematics in Pisa, Calcolo publishes original contributions on numerical analysis and its applications, and on the theory of computation. The main focus of the journal is on numerical linear algebra, approximation theory and its applications, numerical solutions of differential and integral equations, computational complexity, algorithmics, mathematical aspects of computer science, and optimization theory. Calcolo contains expository papers which introduce emerging topics as well as a abstracts of PhD theses, news and reports from conferences, and book reviews.

In this article we consider the numerical approximation of lower bounds of Structured Singular Values, SSV. The SSV is a wellknown mathematical quantity which is widely used to analyse and syntesize the robust stability and instability analysis of linear feedback systems in control theory. It links a bridge between numerical linear algebra and system theory. The computation of lower bounds of SSV by means of ordinary differential equations based technique is presented. The obtained numerical results for the lower bounds of SSV are compared with the well-known MATLAB function mussv available in MATLAB control toolbox.

International audience; The acoustic study of propellers in a hydrodynamic tunnel is of paramount importance during the design process, but can involve significant difficulties due to the boundary layer noise (BLN). Indeed, advanced denoising methods are needed to recover the acoustic signal in case of poor signal-to-noise ratio. The technique proposed in this paper is based on the decomposition of the wall-pressure cross-spectral matrix (CSM) by taking advantage of both the low-rank property of the acoustic CSM and the sparse property of the BLN CSM. Thus, the algorithm belongs to the class of robust principal component analysis (RPCA), which derives from the widely used principal component analysis. If the BLN is spatially decorrelated, the proposed RPCA algorithm can blindly recover the acoustical signals even for negative signal-to-noise ratio. Unfortunately, in a realistic case, acoustic signals recorded in a hydrodynamic tunnel show that the noise may be partially correlated. A prewhitening strategy is then considered in order to take into account the spatially coherent background noise. Numerical simulations and experimental results show an improvement in terms of BLN reduction in the large hydrodynamic tunnel. The effectiveness of the denoising method is also investigated in the context of acoustic source localization.

In this paper, we study the ensemble Kalman filter (EnKF) method for chemical species simulation in air quality forecast data assimilation. The main contribution of this paper is that we study the sparse observation data and make use of the matrix structure of the EnKF update equations to design an algorithm for the purpose of computing the analysis of chemical species in an air quality forecast system efficiently. The proposed method can also handle the combined observations from multiple chemical species together. We applied the proposed method and tested its performance in real air quality data assimilation. Numerical examples are presented to demonstrate the efficiency of the proposed computation method for EnKF updating and the effectiveness of the proposed method for NO2, NO, CO, SO2, O3, PM2.5, and PM10 prediction in air quality forecast data assimilation.

In this paper, the resonant cavity problem with anisotropic fully conducting media, complex geometrical structure and perfect electric conductor walls is investigated. We solve this problem based on the finite element method (FEM) with tangential and linear normal (CT/LN) element and standard linear element. An effective numerical method is proposed by us such that it is free of nonphysical modes. After the FEM discretization, we need to solve a quadratic algebraic eigenvalue problem with a linear constraint condition. In order to overcome this difficulty in the field of numerical algebra, we change this algebraic eigenvalue problem into a generalized eigenvalue problem by introducing an auxiliary zero eigenvector. Moreover, when the permittivity and conductivity are two constants, both the eigenmodes of infinite algebraic multiplicity and all the nonphysical modes are also removed by linearization method. Several numerical experiments show that computational method in this paper can suppress all the spurious modes.

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.

The numerical rank determination frequently occurs in matrix computation when the conventional exact rank of a hidden matrix is desired to be recovered. This paper presents a Matlab package RankRev that implements two efficient algorithms for computing the numerical rank and numerical subspaces of a matrix along with updating/downdating capabilities for making adjustment to the results when a row or column is inserted/deleted. The package and the underlying algorithms are accurate, reliable, and much more efficient than the singular value decomposition when the matrix is of low rank or low nullity.