Search

This work considers the low‐rank approximation of a matrix A(t)$$ A(t) $$ depending on a parameter t$$ t $$ in a compact set D⊂ℝd$$ D\subset {\mathbb{R}}^d $$. Application areas that give rise to such problems include computational statistics and dynamical systems. Randomized algorithms are an increasingly popular approach for performing low‐rank approximation and they usually proceed by multiplying the matrix with random dimension reduction matrices (DRMs). Applying such algorithms directly to A(t)$$ A(t) $$ would involve different, independent DRMs for every t$$ t $$, which is not only expensive but also leads to inherently non‐smooth approximations. In this work, we propose to use constant DRMs, that is, A(t)$$ A(t) $$ is multiplied with the same DRM for every t$$ t $$. The resulting parameter‐dependent extensions of two popular randomized algorithms, the randomized singular value decomposition and the generalized Nyström method, are computationally attractive, especially when A(t)$$ A(t) $$ admits an affine linear decomposition with respect to t$$ t $$. We perform a probabilistic analysis for both algorithms, deriving bounds on the expected value as well as failure probabilities for the approximation error when using Gaussian random DRMs. Both, the theoretical results and numerical experiments, show that the use of constant DRMs does not impair their effectiveness; our methods reliably return quasi‐best low‐rank approximations. [ABSTRACT FROM AUTHOR]

Learning controllers from data for stabilizing dynamical systems typically follows a two-step process of first identifying a model and then constructing a controller based on the identified model. However, learning models means identifying generic descriptions of the dynamics of systems, which can require large amounts of data and extracting information that are unnecessary for the specific task of stabilization. The contribution of this work is to show that if a linear dynamical system has dimension (McMillan degree) n , then there always exist n states from which a stabilizing feedback controller can be constructed, independent of the dimension of the representation of the observed states and the number of inputs. By building on previous work, this finding implies that any linear dynamical system can be stabilized from fewer observed states than the minimal number of states required for learning a model of the dynamics. The theoretical findings are demonstrated with numerical experiments that show the stabilization of the flow behind a cylinder from less data than necessary for learning a model. [ABSTRACT FROM AUTHOR]

This paper develops a class of algorithms for general linear systems and eigenvalue problems. These algorithms apply fast randomized dimension reduction ("sketching") to accelerate standard subspace projection methods, such as GMRES and Rayleigh-Ritz. This modification makes it possible to incorporate nontraditional bases for the approximation subspace that are easier to construct. When the basis is numerically full rank, the new algorithms have accuracy similar to classic methods but run faster and may use less storage. For model problems, numerical experiments show large advantages over the optimized MATLAB routines, including a 70x speedup over gmres and a 10x speedup over eigs. [ABSTRACT FROM AUTHOR]

In a modern observational study based on healthcare databases, the number of observations and of predictors typically range in the order of 105-106 and of 104-105. Despite the large sample size, data rarely provide sufficient information to reliably estimate such a large number of parameters. Sparse regression techniques provide potential solutions, one notable approach being the Bayesian method based on shrinkage priors. In the large n and large p setting, however, the required posterior computation encounters a bottleneck at repeated sampling from a high-dimensional Gaussian distribution, whose precision matrix Φ is expensive to compute and factorize. In this article, we present a novel algorithm to speed up this bottleneck based on the following observation: We can cheaply generate a random vector b such that the solution to the linear system Φβ = b has the desired Gaussian distribution. We can then solve the linear system by the conjugate gradient (CG) algorithm through matrix-vector multiplications by Φ; this involves no explicit factorization or calculation of Φ itself. Rapid convergence of CG in this context is guaranteed by the theory of prior-preconditioning we develop. We apply our algorithm to a clinically relevant large-scale observational study with n = 72,489 patients and p = 22,175 clinical covariates, designed to assess the relative risk of adverse events from two alternative blood anti-coagulants. Our algorithm demonstrates an order of magnitude speed-up in posterior inference, in our case cutting the computation time from two weeks to less than a day. Supplementary materials for this article are available online.

We introduce a mixed precision algorithm for computing sparse matrix-vector products and use it to accelerate the solution of sparse linear systems by iterative methods. Our approach is based on the idea of adapting the precision of each matrix element to their magnitude: we split the elements into buckets and use progressively lower precisions for the buckets of progressively smaller elements. We carry out a rounding error analysis of this algorithm that provides us with an explicit rule to decide which element goes into which bucket and allows us to rigorously control the accuracy of the algorithm. We implement the algorithm on a multicore computer and obtain significant speedups (up to a factor 7\times) with respect to uniform precision algorithms, without loss of accuracy, on a range of sparse matrices from real-life applications. We showcase the effectiveness of our algorithm by plugging it into various Krylov solvers for sparse linear systems and observe that the convergence of the solution is essentially unaffected by the use of adaptive precision. [ABSTRACT FROM AUTHOR]

Block low-rank (BLR) compression can significantly reduce the memory and time costs of parallel sparse direct solvers. In this paper, we investigate the performance of the BLR triangular solve phase, which we observe to be underwhelming when dealing with many right-hand sides (RHS). We explain that this is because the bottleneck of the triangular solve is not in accessing the BLR LU factors, but rather in accessing the RHS, which are uncompressed. Motivated by this finding, we propose several new hybrid variants, which combine the right-looking and left-looking communication patterns to minimize the number of accesses to the RHS. We confirm via a theoretical analysis that these new variants can significantly reduce the total communication volume. We assess the impact of this reduction on the time performance on a range of real-life applications using the MUMPS solver, obtaining up to 20% time reduction. [ABSTRACT FROM AUTHOR]

In this paper, we present a numerical linear algebra analytical study of some schemes for the Bertozzi–Esedoglu–Gillette–Cahn–Hilliard equation. Both 1D and 2D finite difference discretizations in space are proposed with semi-implicit and implicit discretizations on time. We prove that our proposed numerical solutions converge to continuous solutions. [ABSTRACT FROM AUTHOR]

The application of scientific engineering computing based on exascale computing not only offers opportunities but also creates challenges for the development of numerical linear algebra algorithms. Firstly, the characteristics of exascale computing are analyzed, including: parallel programming for large-scale heterogeneous parallel architecture has become the mainstream approach; reducing the extremely high energy costs associated with running large- scale applications is a major concern; multi-precision heterogeneous computing hardware has triggered further research of mixed precision computing. Secondly, the optimization work of mainstream dense and sparse linear algebra solvers for high-performance computing architectures is reviewed, and the characteristics and advantages of each solver are compared. Then, the main technology progress of linear algebra solvers is summarized, mainly including: isolating heterogeneous computing modules and designing a new unified programming framework to achieve performance portability of software algorithms; improving the performance level of numerical computing and data storage using mixed precision methods while ensuring the overall requirements of scientific engineering computing applications; combined with hardware multi-level cache and network communication characteristics, advanced parallel computing algorithms are developed to avoid or reduce inefficient large- scale data communication. Finally, this paper provides an outlook on the future research trends in this direction. [ABSTRACT FROM AUTHOR]

Random butterfly matrices were introduced by Parker in 1995 to remove the need for pivoting when using Gaussian elimination. The growing applications of butterfly matrices have often eclipsed the mathematical understanding of how or why butterfly matrices are able to accomplish these given tasks. To help begin to close this gap using theoretical and numerical approaches, we explore the impact on the growth factor of preconditioning a linear system by butterfly matrices. These results are compared to other common methods found in randomized numerical linear algebra. In these experiments, we show that preconditioning using butterfly matrices has a more significant dampening impact on large growth factors than other common preconditioned and a smaller increase to minimal growth factor systems. Moreover, we are able to determine the full distribution of the growth factors for a subclass of random butterfly matrices. Previous results by Trefethen and Schreiber relating to the distribution of random growth factors were limited to empirical estimates of the first moment for Ginibre matrices. [ABSTRACT FROM AUTHOR]

We introduce a novel approach to exploit mixed precision arithmetic for low-rank approximations. Our approach is based on the observation that singular vectors associated with small singular values can be stored in lower precisions while preserving high accuracy overall. We provide an explicit criterion to determine which level of precision is needed for each singular vector. We apply this approach to block low-rank (BLR) matrices, most of whose off-diagonal blocks have low rank. We propose a new BLR LU factorization algorithm that exploits the mixed precision representation of the blocks. We carry out the rounding error analysis of this algorithm and prove that the use of mixed precision arithmetic does not compromise the numerical stability of the BLR LU factorization. Moreover, our analysis determines which level of precision is needed for each floating-point operation (flop), and therefore guides us toward an implementation that is both robust and efficient. We evaluate the potential of this new algorithm on a range of matrices coming from real-life problems in industrial and academic applications. We show that a large fraction of the entries in the LU factors and flops to perform the BLR LU factorization can be safely switched to lower precisions, leading to significant reductions of the storage and expected time costs, of up to a factor three using fp64, fp32, and bfloat16 arithmetics. [ABSTRACT FROM AUTHOR]

This paper presents an adaptive randomized algorithm for computing the butterfly factorization of an m × n matrix with m ≈ n provided that both the matrix and its transpose can be rapidly applied to arbitrary vectors. The resulting factorization is composed of O(log n) sparse factors, each containing O(n) nonzero entries. The factorization can be attained using O(n3/2 log n) computation and O(n log n) memory resources. The proposed algorithm can be implemented in parallel and can apply to matrices with strong or weak admissibility conditions arising from surface integral equation solvers as well as multi-frontal-based finite-difference, finite-element, or finite-volume solvers. A distributed-memory parallel implementation of the algorithm demonstrates excellent scaling behavior.

Modern GPUs equipped with mixed precision tensor core units present great potential to accelerate dense linear algebra operations such as LU factorization. However, state-of-the-art mixed half/single precision LU factorization algorithms all require the matrix to be stored in single precision, leading to expensive data movement and storage costs. This is explained by the fact that simply switching the storage precision from single to half leads to significant loss of accuracy, forfeiting all accuracy benefits from using tensor core technology. In this article, we propose a new factorization algorithm that is able to store the matrix in half precision without incurring any significant loss of accuracy. Our approach is based on a left-looking scheme employing single precision buffers of controlled size and a mixed precision doubly partitioned algorithm exploiting tensor cores in the panel factorizations. Our numerical results show that compared with the state of the art, the proposed approach is of similar accuracy but with only half the data movement and memory footprint, and hence potentially much faster: it achieves up to 2× and 3.5× speedups on V100 and A100 GPUs, respectively. [ABSTRACT FROM AUTHOR]

The target of this article is to propose a novel numerical scheme that can compete its corresponding solvers of the same type for computing the sign matrix. To do this, a three‐step root solver is designed to gain as much as possible of global convergence when finding the matrix sign. It is proved that the novel numerical scheme is of quartical convergence order. The stability of the scheme is discussed in detail as well. Finally, numerical tests are provided to uphold the theoretical findings. [ABSTRACT FROM AUTHOR]

We propose a new approach to generate a reliable reduced model for a parametric elliptic problem, in the presence of noisy data. The reference model reduction procedure is the directional HiPOD method, which combines Hierarchical Model reduction with a standard Proper Orthogonal Decomposition, according to an offline/online paradigm. In this paper we show that directional HiPOD looses in terms of accuracy when problem data are affected by noise. This is due to the interpolation driving the online phase, since it replicates, by definition, the noise trend. To overcome this limit, we replace interpolation with Machine Learning fitting models which better discriminate relevant physical features in the data from irrelevant unstructured noise. The numerical assessment, although preliminary, confirms the potentialities of the new approach. [ABSTRACT FROM AUTHOR]

In multiword arithmetic, a matrix is represented as the unevaluated sum of two or more lower precision matrices, and a matrix product is formed by multiplying the constituents in low precision. We investigate the use of multiword arithmetic for improving the performance-accuracy tradeoff of matrix multiplication with mixed precision block fused multiply--add (FMA) hardware, focusing especially on the tensor cores available on NVIDIA GPUs. Building on a general block FMA framework, we develop a comprehensive error analysis of multiword matrix multiplication. After confirming the theoretical error bounds experimentally by simulating low precision in software, we use the cuBLAS and CUTLASS libraries to implement a number of matrix multiplication algorithms using double-fp16 (double-binary16) arithmetic. When running the algorithms on NVIDIA V100 and A100 GPUs, we find that double-fp16 is not as accurate as fp32 (binary32) arithmetic despite satisfying the same worst-case error bound. Using probabilistic error analysis, we explain why this issue is likely to be caused by the rounding mode used by the NVIDIA tensor cores, and we propose a parameterized blocked summation algorithm that alleviates the problem and significantly improves the performance-accuracy tradeoff. [ABSTRACT FROM AUTHOR]

In this paper, we present a numerical linear algebra analytical study of some schemes for the Bertozzi–Esedoglu–Gillette–Cahn–Hilliard equation. Both 1D and 2D finite difference discretizations in space are proposed with semi-implicit and implicit discretizations on time. We prove that our proposed numerical solutions converge to continuous solutions.

• The INTERNODES method is applied to problems in contact mechanics. • We propose a highly efficient preconditioner to solve the sequence of linear systems. • An algorithm is suggested to cheaply solve linear systems with the preconditioning matrix. • The quality of the preconditioner is assessed on small to medium size 2D problems. The mortar finite element method is a well-established method for the numerical solution of partial differential equations on domains displaying non-conforming interfaces. The method is known for its application in computational contact mechanics. However, its implementation remains challenging as it relies on geometrical projections and unconventional quadrature rules. The INTERNODES (INTERpolation for NOn-conforming DEcompositionS) method, instead, could overcome the implementation difficulties thanks to flexible interpolation techniques. Moreover, it was shown to be at least as accurate as the mortar method making it a very promising alternative for solving problems in contact mechanics. Unfortunately, in such situations the method requires solving a sequence of ill-conditioned linear systems. In this paper, preconditioning techniques are designed and implemented for the efficient solution of those linear systems. [ABSTRACT FROM AUTHOR]

Point cloud (PC) is a collection of discrete geometric samples of a physical object in 3D space. A PC video consists of temporal frames evenly spaced in time, each containing a static PC at one time instant. PCs in adjacent frames typically do not have point-to-point (P2P) correspondence, and thus exploiting temporal redundancy for PC restoration across frames is difficult. In this paper, we focus on the super-resolution (SR) problem for PC video: increase point density of PCs in video frames while preserving salient geometric features consistently across time. We accomplish this with two ideas. First, we establish partial P2P coupling between PCs of adjacent frames by interpolating interior points in a low-resolution PC patch in frame $t$ and translating them to a corresponding patch in frame $t+1$ , via a motion model computed by iterative closest point (ICP). Second, we promote piecewise smoothness in 3D geometry in each patch using feature graph Laplacian regularizer (FGLR) in an easily computable quadratic form. The two ideas translate to an unconstrained quadratic programming (QP) problem with a system of linear equations as solution—one where we ensure the numerical stability by upper-bounding the condition number of the coefficient matrix. Finally, to improve the accuracy of the ICP motion model, we re-sample points in a super-resolved patch at time $t$ to better match a low-resolution patch at time $t+1$ via bipartite graph matching after each SR iteration. Experimental results show temporally consistent super-resolved PC videos generated by our scheme, outperforming SR competitors that optimized on a per-frame basis, in two established PC metrics. [ABSTRACT FROM AUTHOR]

Randomized singular value decomposition (RSVD) is by now a well-established technique for efficiently computing an approximate singular value decomposition of a matrix. Building on the ideas that underpin RSVD, the recently proposed algorithm "randUTV" computes a full factorization of a given matrix that provides low-rank approximations with near-optimal error. Because the bulk of randUTV is cast in terms of communication-efficient operations such as matrix-matrix multiplication and unpivoted QR factorizations, it is faster than competing rank-revealing factorization methods such as column-pivoted QR in most high-performance computational settings. In this article, optimized randUTV implementations are presented for both shared-memory and distributed-memory computing environments. For shared memory, randUTV is redesigned in terms of an algorithm-by-blocks that, together with a runtime task scheduler, eliminates bottlenecks from data synchronization points to achieve acceleration over the standard blocked algorithm based on a purely fork-join approach. The distributed-memory implementation is based on the ScaLAPACK library. The performance of our new codes compares favorably with competing factorizations available on both shared-memory and distributed-memory architectures. [ABSTRACT FROM AUTHOR]

This paper concerns test matrices for numerical linear algebra using an error-free transformation of floating-point arithmetic. For specified eigenvalues given by a user, we propose methods of generating a matrix whose eigenvalues are exactly known based on, for example, Schur or Jordan normal form and a block diagonal form. It is also possible to produce a real matrix with specified complex eigenvalues. Such test matrices with exactly known eigenvalues are useful for numerical algorithms in checking the accuracy of computed results. In particular, exact errors of eigenvalues can be monitored. To generate test matrices, we first propose an error-free transformation for the product of three matrices YSX. We approximate S by S ′ to compute Y S ′ X without a rounding error. Next, the error-free transformation is applied to the generation of test matrices with exactly known eigenvalues. Note that the exactly known eigenvalues of the constructed matrix may differ from the anticipated given eigenvalues. Finally, numerical examples are introduced in checking the accuracy of numerical computations for symmetric and unsymmetric eigenvalue problems. [ABSTRACT FROM AUTHOR]

In this paper, we argue for the utility of deterministic inference in the classical problem of numerical linear algebra, that of solving a linear system. We show how the Gaussian belief propagation solver, known to work for symmetric matrices, can be modified to handle nonsymmetric matrices. Furthermore, we introduce a new algorithm for matrix inversion that corresponds to the generalized belief propagation derived from the cluster variation method. We relate these algorithms to LU and block LU decompositions and provide certain guarantees based on theorems from the theory of belief propagation. The proposed solvers can be accelerated within a framework of geometric multigrid. We benchmark the resulting multigrid method with a set of classical multigrid solvers based on incomplete LU, Chebyshev, and pattern Gauss-Seidel smoothers. The results indicate that Gaussian belief propagation can be used to construct competitive solvers, smoothers, and preconditioners for selected elliptic test problems. [ABSTRACT FROM AUTHOR]

In this paper, algorithms are presented constructing, for the Hadamard maximum determinant problem of dimension n with elements ± 1. This leads to the construction of efficient algorithms specifying the determinants of n × n matrices with elements ± 1. The notion of lexicographically ordered sequences of integers and matrices is presented and we describe algorithms creating all possible lexicographically ordered vectors of dimension. With the ordered vectors, we create all the alternative determinants and search for their values, so the ones we meet first are the lexically smaller ones. The implementation of the algorithms is introduced, both Exhaustive Search, Brute force Algorithm and Partial Search and we compare these algorithms according to their speed and efficiency. Furthermore, we present some theoretical constructions attaining the evaluation of the spectrum of matrices with elements ± 1. Summarizing in this article, we develop theoretical techniques on which we can rely and build faster algorithms to find the determinants of ± 1 ± 1. matrices of the spectrum. [ABSTRACT FROM AUTHOR]

Sparse matrix-matrix multiplication (or SpGEMM) is a key primitive for many high-performance graph algorithms as well as for some linear solvers, such as algebraic multigrid. The scaling of existing parallel implementations of SpGEMM is heavily bound by communication. Even though 3D (or 2.5D) algorithms have been proposed and theoretically analyzed in the at MPI model on Erd}os{Rffenyi matrices, those algorithms had not been implemented in practice and their complexities had not been analyzed for the general case. In this work, we present the first implementation of the 3D SpGEMM formulation that exploits multiple (intranode and internode) levels of parallelism, achieving significant speedups over the state-of-the-art publicly available codes at all levels of concurrencies. We extensively evaluate our implementation and identify bottlenecks that should be subject to further research.

This article is an extension of work entitled “Computing planetary interior normal modes with a highly parallel polynomial filtering eigensolver.” by Shi et al., originally presented at the SC18 conference. A highly parallel polynomial filtered eigensolver was developed and exploited to calculate the planetary normal modes. The proposed method is ideally suited for computing interior eigenpairs for large-scale eigenvalue problems as it greatly enhances memory and computational efficiency. In this article, the second-order finite element method is used to further improve the accuracy as only the first-order finite element method was deployed in the previous work. The parallel algorithm, its parallel performance up to 20k processors, and the great computational accuracy are illustrated. The reproducibility of the previous work was successfully performed on the Student Cluster Competition at the SC19 conference by several participant teams using a completely different Mars-model dataset on different clusters. Both weak and strong scaling performances of the reproducibility by the participant teams were impressive and encouraging. The analysis and reflection of their results are demonstrated and future direction is discussed. [ABSTRACT FROM AUTHOR]

Repeatedly performing sparse matrix‐vector multiplication (SpMV) followed by transposed sparse matrix‐vector multiplication (SpMTV) with the same matrix is a part of several algorithms, for example, the Lanczos biorthogonalization algorithm and the biconjugate gradient method. Such algorithms can benefit from combining parallel SpMV and SpMTV into a single operation we call joint direct and transposed sparse matrix‐vector multiplication (SpMMTV). In this article, we present a parallel SpMMTV algorithm for shared‐memory CPUs. The algorithm uses a sparse matrix format that divides the stored matrix into sparse matrix blocks and compresses the row and column indices of the matrix. This sparse matrix format can be also used for SpMV, SpMTV, and similar sparse matrix‐vector operations. We expand upon existing research by suggesting new variants of the parallel SpMMTV algorithm and by extending the algorithm to efficiently support symmetric matrices. We compare the performance of the presented parallel SpMMTV algorithm with alternative approaches, which use state‐of‐the‐art sparse matrix formats and libraries, using sparse matrices from real‐world applications. The performance results indicate that the median performance of our proposed parallel SpMMTV algorithm is up to 45% higher than of the alternative approaches. [ABSTRACT FROM AUTHOR]

We investigate a special class of data sparse rank-structured matrices that combine a flat block low-rank (BLR) partitioning with the use of shared (called nested in the hierarchical case) bases. This format is to H 2 matrices what BLR is to H matrices: we therefore call it the BLR² matrix format. We present algorithms for the construction and LU factorization of BLR² matrices,and perform their cost analysis---both asymptotically and for a fixed problem size. With weak admissibility, BLR² matrices reduce to block separable matrices (the flat version of HBS/HSS). Our analysis and numerical experiments reveal some limitations of BLR² matrices with weak admissibility,which we propose to overcome with two approaches: strong admissibility, and the use of multiple shared bases per row and column. [ABSTRACT FROM AUTHOR]

This paper presents an adaptive randomized algorithm for computing the butterfly factorization of an m x n matrix with m ≈ n provided that both the matrix and its transpose can be rapidly applied to arbitrary vectors. The resulting factorization is composed of O(log n) sparse factors, each containing O(n) nonzero entries. The factorization can be attained using O(n3/2 log n) computation and O(n log n) memory resources. The proposed algorithm can be implemented in parallel and can apply to matrices with strong or weak admissibilit y conditions arising from surface integral equation solvers as well as multi-frontal-based finite-difference, finite-element, or finite-volume solvers. A distributed-memory parallel implementation of the algorithm demonstrates excellent scaling behavior. [ABSTRACT FROM AUTHOR]

Stochastic rounding rounds a real number to the next larger or smaller floating-point number with probabilities 1 minus the relative distances to those numbers. It is gaining attention in deep learning because it can increase the success of low precision computations. We compare basic properties of stochastic rounding with those for round to nearest, finding properties in common as well as significant differences. We prove that for stochastic rounding the rounding errors are mean independent random variables with zero mean. We derive a new version of our probabilistic error analysis theorem from [N. J. Higham and T. Mary, SIAM J. Sci. Comput., 41 (2019), pp. A2815-- A2835], weakening the assumption of independence of the random variables to mean independence. These results imply that for a wide range of linear algebra computations the backward error for stochastic rounding is unconditionally bounded by a multiple of √nu to first order, with a certain probability, where n is the problem size and u is the unit roundoff. This is the first scenario where the rule of thumb that one can replace nu by √nu in a rounding error bound has been shown to hold without any additional assumptions on the rounding errors. We also explain how stochastic rounding avoids the phenomenon of stagnation in sums, whereby small addends are obliterated by round to nearest when they are too small relative to the sum. [ABSTRACT FROM AUTHOR]

We give two different and simple constructions for dimensionality reduction in l2 via linear mappings that are sparse: only an (ε)-fraction of entries in each column of our embedding matrices are non-zero to achieve distortion 1+ε with high probability, while still achieving the asymptotically optimal number of rows. These are the first constructions to provide subconstant sparsity for all values of parameters, improving upon previous works of Achlioptas [2003] and Dasgupta et al. [2010]. Such distributions can be used to speed up applications where l2 dimensionality reduction is used. [ABSTRACT FROM AUTHOR]