Search

The paper describes a sparse direct solver for the linear systems that arise from the discretization of an elliptic PDE on a two-dimensional domain. The scheme decomposes the domain into thin subdomains, or “slabs” and uses a two-level approach that is designed with parallelization in mind. The scheme takes advantage of H 2 -matrix structure emerging during factorization and utilizes randomized algorithms to efficiently recover this structure. As opposed to multi-level nested dissection schemes that incorporate the use of H or H 2 matrices for a hierarchy of front sizes, SlabLU is a two-level scheme which only uses H 2 -matrix algebra for fronts of roughly the same size. The simplicity allows the scheme to be easily tuned for performance on modern architectures and GPUs. The solver described is compatible with a range of different local discretizations, and numerical experiments demonstrate its performance for regular discretizations of rectangular and curved geometries. The technique becomes particularly efficient when combined with very high-order accurate multidomain spectral collocation schemes. With this discretization, a Helmholtz problem on a domain of size 1000 λ × 1000 λ (for which N = 100 M ) is solved in 15 min to 6 correct digits on a high-powered desktop with GPU acceleration. [ABSTRACT FROM AUTHOR]

We consider the problem of state estimation from a few linear measurements, where the state to recover is an element of the manifold M of solutions of a parameter-dependent equation. The state is estimated using prior knowledge on M coming from model order reduction. Variational approaches based on linear approximation of M , such as PBDW, yield a recovery error limited by the Kolmogorov width of M . To overcome this issue, piecewise-affine approximations of M have also been considered, that consist in using a library of linear spaces among which one is selected by minimizing some distance to M . In this paper, we propose a state estimation method relying on dictionary-based model reduction, where space is selected from a library generated by a dictionary of snapshots, using a distance to the manifold. The selection is performed among a set of candidate spaces obtained from a set of ℓ 1 -regularized least-squares problems. Then, in the framework of parameter-dependent operator equations (or PDEs) with affine parametrizations, we provide an efficient offline-online decomposition based on randomized linear algebra, that ensures efficient and stable computations while preserving theoretical guarantees. [ABSTRACT FROM AUTHOR]

We present an acceleration method for sequences of large-scale linear systems, such as the ones arising from the numerical solution of time-dependent partial differential equations coupled with algebraic constraints. We discuss different approaches to leverage the subspace containing the history of solutions computed at previous time steps in order to generate a good initial guess for the iterative solver. In particular, we propose a novel combination of reduced-order projection with randomized linear algebra techniques, which drastically reduces the number of iterations needed for convergence. We analyze the accuracy of the initial guess produced by the reduced-order projection when the coefficients of the linear system depend analytically on time. Extending extrapolation results by Demanet and Townsend to a vector-valued setting, we show that the accuracy improves rapidly as the size of the history increases, a theoretical result confirmed by our numerical observations. In particular, we apply the developed method to the simulation of plasma turbulence in the boundary of a fusion device, showing that the time needed for solving the linear systems is significantly reduced. [ABSTRACT FROM AUTHOR]

Structure–function relationships occur throughout the sciences. Motivated by optimization of such systems, we develop a framework for estimating the input modes from the singular value decomposition from the action of the forward operator alone. These can then be used to determine the input (structure) changes, which induce the largest output (function) changes. The accuracy of the estimate is determined by reference to the method of snapshots. The proposed method is demonstrated on several example problems, and finally used to approximate the optimal airway treatment set for a problem in respiratory physiology. [ABSTRACT FROM AUTHOR]

Several quantum algorithms for linear algebra problems, and in particular quantum machine learning problems, have been “dequantized” in the past few years. These dequantization results typically hold when classical algorithms can access the data via length-squared sampling. This assumption, which is standard in the field of randomized linear algebra, means that for a unit-norm vector u ∈ C n , we can sample from the distribution p u : { 1 , … , n } → [ 0 , 1 ] defined as p u (i) = | u (i) | 2 for each i ∈ { 1 , … , n } . Since this distribution corresponds to the distribution obtained by measuring the quantum state | u > in the computational basis, length-squared sampling access gives a reasonable classical analogue to the kind of quantum access considered in many quantum algorithms for linear algebra problems. In this work we investigate how robust these dequantization results are. We introduce the notion of approximate length-squared sampling, where classical algorithms are only able to sample from a distribution close to the ideal distribution in total variation distance. While quantum algorithms are natively robust against small perturbations, current techniques in dequantization are not. Our main technical contribution is showing how many techniques from randomized linear algebra can be adapted to work under this weaker assumption as well. We then use these techniques to show that the recent low-rank dequantization framework by Chia, Gilyén, Li, Lin, Tang and Wang (JACM 2022) and the dequantization framework for sparse matrices by Gharibian and Le Gall (STOC 2022), which are both based on the Quantum Singular Value Transformation, can be generalized to the case of approximate length-squared sampling access to the input. We also apply these results to obtain a robust dequantization of many quantum machine learning algorithms, including quantum algorithms for recommendation systems, supervised clustering and low-rank matrix inversion. [ABSTRACT FROM AUTHOR]

The Nyström method is a popular choice for finding a low-rank approximation to a symmetric positive semidefinite matrix. The method can fail when applied to symmetric indefinite matrices for which the error can be unboundedly large. In this work, we first identify the main challenges in finding a Nyström approximation to symmetric indefinite matrices. We then prove the existence of a variant that overcomes the instability, and establish relative-error nuclear norm bounds of the resulting approximation that hold when the singular values decay rapidly. The analysis naturally leads to a practical algorithm, whose robustness is illustrated with experiments. [ABSTRACT FROM AUTHOR]

We present a sequential version of the multilinear Nyström algorithm which is suitable for low-rank Tucker approximation of tensors given in a streaming format. Accessing the tensor A exclusively through random sketches of the original data, the algorithm effectively leverages structures in A , such as low-rankness, and linear combinations. We present a deterministic analysis of the algorithm and demonstrate its superior speed and efficiency in numerical experiments including an application in video processing. • Enhances multilinear Nyström (MLN) with sequential truncation in streaming model. • Undertakes error analysis of the proposed streaming modification of MLN. • Illustrates performance gains of up to 70 percent without compromising accuracy. [ABSTRACT FROM AUTHOR]

Data assimilation methods are commonly used to predict petrophysical properties of deep saline aquifers for carbon dioxide sequestration studies. However, data assimilation is usually computationally challenging for large-scale geological models with a large number of geophysical observations. A computationally efficient method is developed for large-scale data assimilation problems. In the proposed method, model parameters and observations are reduced by randomized tensor decomposition, and the inversion is performed in lower-dimensional spaces. Tensors are multiway arrays that generalize matrices to multiple dimensions and therefore provide a natural way to represent three-dimensional geological models and geophysical measurements. To exploit the inherent multidimensional structure of tensors, high-order singular value decomposition performs orthogonal decomposition of the tensor in the high-order space instead of using a flattened matrix representation. However, for complex geological models, the dimension of tensors becomes very large, making the decomposition computationally unfeasible. To overcome this limitation, the dimension of the tensors is first reduced by a randomized linear algebra algorithm which guarantees that most of the information of the original tensors is preserved with high probability in the reduced tensor, and then high-order singular value decomposition is applied on the smaller tensors with relatively small computational cost. By using randomized tensor decomposition for model and data re-parameterization, data assimilation is efficiently performed in the low-dimensional model and data space. The proposed method is demonstrated and validated on the three-dimensional data set synthesized based on the Johansen formation, a potential carbon dioxide storage site offshore in Norway. [ABSTRACT FROM AUTHOR]

Linear programming (LP) is an extremely useful tool which has been successfully applied to solve various problems in a wide range of areas, including operations research, engineering, economics, or even more abstract mathematical areas such as combinatorics. It is also used in many machine learning applications, such as '1-regularized SVMs, basis pursuit, nonnegative matrix factorization, etc. Interior Point Methods (IPMs) are one of the most popular methods to solve LPs both in theory and in practice. Their underlying complexity is dominated by the cost of solving a system of linear equations at each iteration. In this paper, we consider both feasible and infeasible IPMs for the special case where the number of variables is much larger than the number of constraints. Using tools from Randomized Linear Algebra, we present a preconditioning technique that, when combined with the iterative solvers such as Conjugate Gradient or Chebyshev Iteration, provably guarantees that IPM algorithms (suitably modified to account for the error incurred by the approximate solver), converge to a feasible, approximately optimal solution, without increasing their iteration complexity. Our empirical evaluations verify our theoretical results on both real-world and synthetic data. [ABSTRACT FROM AUTHOR]

Low rank tensor approximation is a fundamental tool in modern machine learning and data science. In this paper, we study the characterization, perturbation analysis, and an efficient sampling strategy for two primary tensor CUR approximations, namely Chidori and Fiber CUR. We characterize exact tensor CUR decompositions for low multilinear rank tensors. We also present theoretical error bounds of the tensor CUR approximations when (adversarial or Gaussian) noise appears. Moreover, we show that low cost uniform sampling is sufficient for tensor CUR approximations if the tensor has an incoherent structure. Empirical performance evaluations, with both synthetic and real-world datasets, establish the speed advantage of the tensor CUR approximations over other state-of-the-art low multilinear rank tensor approximations. [ABSTRACT FROM AUTHOR]

We present a method for converting tensors into the tensor train format based on actions of the tensor as a vector-valued multilinear function. Existing methods for constructing tensor trains require access to "array entries" of the tensor and are therefore inefficient or computationally prohibitive if the tensor is accessible only through its action, especially for high order tensors. Our method permits efficient tensor train compression of large high order derivative tensors for nonlinear mappings that are implicitly defined through the solution of a system of equations. Array entries of these derivative tensors are not directly accessible, but actions of these tensors can be computed efficiently via a procedure that we discuss. Such tensors are often amenable to tensor train compression in theory, but until now no efficient algorithm existed to convert them into tensor train format. We demonstrate our method by compressing a Hilbert tensor of size 41 x 42 x 43 x 44 x 45, and by forming high order (up to fifth order derivatives/sixth order tensors) Taylor series surrogates of the noise-whitened parameter-to-output map for a stochastic partial differential equation with boundary output. [ABSTRACT FROM AUTHOR]

We introduce an algorithm for the least squares solution of a rectangular linear system Ax = b, in which A may be arbitrarily ill-conditioned. We assume that a complementary matrix Z is known such that A â?' AZ* A is numerically low rank. Loosely speaking, Z* acts like a generalized inverse of A up to a numerically low-rank error. We give several examples of (A,Z) combinations in function approximation, where we can achieve high-order approximations in a number of nonstandard settings: the approximation of functions on domains with irregular shapes, weighted least squares problems with highly skewed weights, and the spectral approximation of functions with localized singularities. The algorithm is most efficient when A and Z* have fast matrix-vector multiplication and when the numerical rank of A â?' AZ* A is small. [ABSTRACT FROM AUTHOR]

We study the problem of actively imaging a range-limited far-field scene using an antenna array. We describe how the range limit imposes structure in the measurements across multiple wavelengths. This structure allows us to introduce a novel trade-off: the number of spatial array measurements (i.e., beams that have to be formed) can be reduced to be significantly lower than the number array elements if the scene is illuminated with a broadband source. To take advantage of this trade-off, we use a small number of "generic" linear combinations of the array outputs, instead of the phase offsets used in conventional beamforming. We provide theoretical justification for the proposed trade-off without making any strong structural assumptions on the target scene (such as sparsity) except that it is range-limited. In proving our theoretical results, we take inspiration from the sketching literature. We also provide simulation results to establish the merit of the proposed signal acquisition strategy. Our proposed method results in a reduction in the number of required spatial measurements in an array imaging system and hence can directly impact their speed and cost of operation. [ABSTRACT FROM AUTHOR]

Abstract Projection-cost preservation is a low-rank approximation guarantee which ensures that the cost of any rank- k projection can be preserved using a smaller sketch of the original data matrix. We present a general structural result outlining four sufficient conditions to achieve projection-cost preservation. These conditions can be satisfied using tools from the Randomized Linear Algebra literature. [ABSTRACT FROM AUTHOR]

The singular value decomposition (SVD) has a crucial role in model order reduction. It is often utilized in the offline stage to compute basis functions that project the high-dimensional nonlinear problem into a low-dimensional model which is then evaluated cheaply. It constitutes a building block for many techniques such as the proper orthogonal decomposition (POD) and dynamic mode decomposition (DMD). The aim of this work is to provide an efficient computation of low-rank POD and/or DMD modes via randomized matrix decompositions. This is possible due to the randomized singular value decomposition (rSVD) which is a fast and accurate alternative of the SVD. Although this is considered an offline stage, this computation may be extremely expensive; therefore, the use of compressed techniques drastically reduce its cost. Numerical examples show the effectiveness of the method for both POD and DMD. [ABSTRACT FROM AUTHOR]

Kernel k-means clustering can correctly identify and extract a far more varied collection of cluster structures than the linear k-means clustering algorithm. However, kernel k- means clustering is computationally expensive when the non-linear feature map is highdimensional and there are many input points. Kernel approximation, e.g., the Nyström method, has been applied in previous works to approximately solve kernel learning problems when both of the above conditions are present. This work analyzes the application of this paradigm to kernel k-means clustering, and shows that applying the linear k-means clustering algorithm to k/∊ (1 + o(1)) features constructed using a so-called rank-restricted Nyström approximation results in cluster assignments that satisfy a 1+∊ approximation ratio in terms of the kernel k-means cost function, relative to the guarantee provided by the same algorithm without the use of the Nyström method. As part of the analysis, this work establishes a novel 1+ relative-error trace norm guarantee for low-rank approximation using the rank-restricted Nyström approximation. Empirical evaluations on the 8:1 million instance MNIST8M dataset demonstrate the scalability and usefulness of kernel k-means clustering with Nyström approximation. This work argues that spectral clustering using Nyström approximation|a popular and computationally effcient, but theoretically unsound approach to non-linear clustering| should be replaced with the effcient and theoretically sound combination of kernel k-means clustering with Nyström approximation. The superior performance of the latter approach is empirically veriöed. [ABSTRACT FROM AUTHOR]

We address the statistical and optimization impacts of the classical sketch and Hessian sketch used to approximately solve the Matrix Ridge Regression (MRR) problem. Prior research has quantified the effects of classical sketch on the strictly simpler least squares regression (LSR) problem. We establish that classical sketch has a similar effect upon the optimization properties of MRR as it does on those of LSR: namely, it recovers nearly optimal solutions. By contrast, Hessian sketch does not have this guarantee; instead, the approximation error is governed by a subtle interplay between the "mass" in the responses and the optimal objective value. For both types of approximation, the regularization in the sketched MRR problem results in significantly different statistical properties from those of the sketched LSR problem. In particular, there is a bias-variance trade-off in sketched MRR that is not present in sketched LSR. We provide upper and lower bounds on the bias and variance of sketched MRR; these bounds show that classical sketch significantly increases the variance, while Hessian sketch significantly increases the bias. Empirically, sketched MRR solutions can have risks that are higher by an order-of-magnitude than those of the optimal MRR solutions. We establish theoretically and empirically that model averaging greatly decreases the gap between the risks of the true and sketched solutions to the MRR problem. Thus, in parallel or distributed settings, sketching combined with model averaging is a powerful technique that quickly obtains near-optimal solutions to the MRR problem while greatly mitigating the increased statistical risk incurred by sketching. [ABSTRACT FROM AUTHOR]

In this era of large-scale data, distributed systems built on top of clusters of commodity hardware provide cheap and reliable storage and scalable processing of massive data. With cheap storage, instead of storing only currently relevant data, it is common to store as much data as possible, hoping that its value can be extracted later. In this way, exabytes (10 ^18 bytes) of data are being created on a daily basis. Extracting value from these data, however, requires scalable implementations of advanced analytical algorithms beyond simple data processing, e.g., statistical regression methods, linear algebra, and optimization algorithms. Most such traditional methods are designed to minimize floating-point operations, which is the dominant cost of in-memory computation on a single machine. In parallel and distributed environments, however, load balancing and communication, including disk and network input/output (I/O) , can easily dominate computation. These factors greatly increase the complexity of algorithm design and challenge traditional ways of thinking about the design of parallel and distributed algorithms. Here, we review recent work on developing and implementing randomized matrix algorithms in large-scale parallel and distributed environments. Randomized algorithms for matrix problems have received a great deal of attention in recent years, thus far typically either in theory or in machine learning applications or with implementations on a single machine. Our main focus is on the underlying theory and practical implementation of random projection and random sampling algorithms for very large very overdetermined (i.e., overconstrained) \ell1- and \ell2 -regression problems. Randomization can be used in one of two related ways: either to construct subsampled problems that can be solved, exactly or approximately, with traditional numerical methods; or to construct preconditioned versions of the original full problem that are easier to solve with traditional iterative algorithms. Theoretical results demonstrate that in near input-sparsity time and with only a few passes through the data one can obtain very strong relative-error approximate solutions, with high probability. Empirical results highlight the importance of various tradeoffs (e.g., between the time to construct an embedding and the conditioning quality of the embedding, between the relative importance of computation versus communication, etc.) and demonstrate that \ell1- and \ell2 -regression problems can be solved to low, medium, or high precision in existing distributed systems on up to terabyte-sized data. [ABSTRACT FROM PUBLISHER]

The development of randomized algorithms for numerical linear algebra, e.g. for computing approximate QR and SVD factorizations, has recently become an intense area of research. This paper studies one of the most frequently discussed algorithms in the literature for dimensionality reduction-specifically for approximating an input matrix with a low-rank element. We introduce a novel and rather intuitive analysis of the algorithm in [], which allows us to derive sharp estimates and give new insights about its performance. This analysis yields theoretical guarantees about the approximation error and at the same time, ultimate limits of performance (lower bounds) showing that our upper bounds are tight. Numerical experiments complement our study and show the tightness of our predictions compared with empirical observations. [ABSTRACT FROM AUTHOR]

In this paper, a nonlinear electrostatic surface integral equation is presented that is suitable for predicting corrosion-related fields. Nonlinear behavior arises due to electrochemical reactions at polarized surfaces. Hierarchical H2 matrices are used to compress the discretized integral equation for the fast solution of large problems. A technique based on randomized linear algebra is discussed for the efficient computation of the Jacobian matrix required at each iteration of a nonlinear solution. [ABSTRACT FROM AUTHOR]

• Time autocorrelation functions and power spectral density are common outcomes of large numerical simulations. • When there are many particles or grid points, only very short lag autocorrelations are possible due to the required storage. • Randomized sketching techniques can be used to reduce the cost of storage. • These new techniques have theoretical guarantees and better practical performance than existing techniques. Sketching is a stochastic dimension reduction method that preserves geometric structures of data and has applications in high-dimensional regression, low rank approximation and graph sparsification. In this work, we show that sketching can be used to compress simulation data and still accurately estimate time autocorrelation and power spectral density. For a given compression ratio, the accuracy is much higher than using previously known methods. In addition to providing theoretical guarantees, we apply sketching to a molecular dynamics simulation of methanol and find that the estimate of spectral density is 90% accurate using only 10% of the data. [ABSTRACT FROM AUTHOR]