Linear Algebra Algorithms Research Articles

Approximation and sampling of multivariate probability distributions in the tensor train decomposition

General multivariate distributions are notoriously expensive to sample from, particularly the high-dimensional posterior distributions in PDE-constrained inverse problems. This paper develops a sampler for arbitrary continuous multivariate distributions that is based on low-rank surrogates in the tensor train format, a methodology that has been exploited for many years for scalable, high-dimensional density function approximation in quantum physics and chemistry. We build upon recent developments of the cross approximation algorithms in linear algebra to construct a tensor train approximation to the target probability density function using a small number of function evaluations. For sufficiently smooth distributions, the storage required for accurate tensor train approximations is moderate, scaling linearly with dimension. In turn, the structure of the tensor train surrogate allows sampling by an efficient conditional distribution method since marginal distributions are computable with linear complexity in dimension. Expected values of non-smooth quantities of interest, with respect to the surrogate distribution, can be estimated using transformed independent uniformly-random seeds that provide Monte Carlo quadrature or transformed points from a quasi-Monte Carlo lattice to give more efficient quasi-Monte Carlo quadrature. Unbiased estimates may be calculated by correcting the transformed random seeds using a Metropolis–Hastings accept/reject step, while the quasi-Monte Carlo quadrature may be corrected either by a control-variate strategy or by importance weighting. We show that the error in the tensor train approximation propagates linearly into the Metropolis–Hastings rejection rate and the integrated autocorrelation time of the resulting Markov chain; thus, the integrated autocorrelation time may be made arbitrarily close to 1, implying that, asymptotic in sample size, the cost per effectively independent sample is one target density evaluation plus the cheap tensor train surrogate proposal that has linear cost with dimension. These methods are demonstrated in three computed examples: fitting failure time of shock absorbers; a PDE-constrained inverse diffusion problem; and sampling from the Rosenbrock distribution. The delayed rejection adaptive Metropolis (DRAM) algorithm is used as a benchmark. In all computed examples, the importance weight-corrected quasi-Monte Carlo quadrature performs best and is more efficient than DRAM by orders of magnitude across a wide range of approximation accuracies and sample sizes. Indeed, all the methods developed here significantly outperform DRAM in all computed examples.

Hemolysis Detection in Sub-Microliter Volumes of Blood Plasma

Hemolysis is one of the main reasons for blood sample rejection in clinical laboratories. The scope of this study is to develop a simple, sensitive, and cost-effective method for hemolysis detection in presence of interferants like bilirubin, and lipids in sub-microliter volumes of blood plasma. Experimental samples were prepared in plasma extracted from whole blood. These samples were typically mixtures of hemoglobin, bilirubin, and lipids of varying concentrations. A multispectral optical setup was devised to probe these samples without any optical filters or moving parts. The hemoglobin content in the samples were measured within the range of 0mg/dL to 400 mg/dL concentration using <1μL of detection volume. Detection sensitivity of >90% accuracy and ∼10% coefficient of variation across 27 unknown samples. The optical attenuation path length used was merely ∼150 micrometers. The spectral interference due to overlapping absorption spectrum of bilirubin and scattering spectrum of lipids are resolved using linear matrix algebra algorithms. Hemolysis can lead to spurious measurements of key clinical parameters such as Potassium, Lactate Dehydrogenase, Aspartic Acid Transferase, and other diagnostic biomarkers. Commonly used visual inspection of blood plasma coloration is prone to variability. Description of the sample preparation, calibration and verification of the experimental setup and linear matrix algebra algorithm for analyte interference determination is reported here. Owing to the sub-microliter detection volume and high sensitivity, the system has realistic potential to be implemented in point of care medical devices that demands such low volumes of clinical specimen.

The Parallel Tiled WZ Factorization Algorithm for Multicore Architectures

Abstract The aim of this paper is to investigate dense linear algebra algorithms on shared memory multicore architectures. The design and implementation of a parallel tiled WZ factorization algorithm which can fully exploit such architectures are presented. Three parallel implementations of the algorithm are studied. The first one relies only on exploiting multithreaded BLAS (basic linear algebra subprograms) operations. The second implementation, except for BLAS operations, employs the OpenMP standard to use the loop-level parallelism. The third implementation, except for BLAS operations, employs the OpenMP task directive with the depend clause. We report the computational performance and the speedup of the parallel tiled WZ factorization algorithm on shared memory multicore architectures for dense square diagonally dominant matrices. Then we compare our parallel implementations with the respective LU factorization from a vendor implemented LAPACK library. We also analyze the numerical accuracy. Two of our implementations can be achieved with near maximal theoretical speedup implied by Amdahl’s law.

Scalable diagnostics for global atmospheric chemistry using Ristretto library (version 1.0)

Abstract. We introduce a new set of algorithmic tools capable of producing scalable, low-rank decompositions of global spatiotemporal atmospheric chemistry data. By exploiting emerging randomized linear algebra algorithms, a suite of decompositions are proposed that extract the dominant features from big data sets (i.e., global atmospheric chemistry at longitude, latitude, and elevation) with improved interpretability. Importantly, our proposed algorithms scale with the intrinsic rank of the global chemistry space rather than the ever increasing spatiotemporal measurement space, thus allowing for the efficient representation and compression of the data. In addition to scalability, two additional innovations are proposed for improved interpretability: (i) a nonnegative decomposition of the data for improved interpretability by constraining the chemical space to have only positive expression values (unlike PCA analysis); and (ii) sparse matrix decompositions, which threshold small weights to zero, thus highlighting the dominant, localized spatial activity (again unlike PCA analysis). Our methods are demonstrated on a full year of global chemistry dynamics data, showing the significant improvement in computational speed and interpretability. We show that the decomposition methods presented here successfully extract known major features of atmospheric chemistry, such as summertime surface pollution and biomass burning activities.

Application of Computer Algebra Methods to Investigate the Dynamics of the System of Two Connected Bodies Moving along a Circular Orbit

Computer algebra methods are used to investigate properties of a nonlinear algebraic system that determines the equilibrium orientations for a system of two bodies connected by a spherical hinge that moves along a circular orbit under the action of gravitational torque. To determine the equilibrium orientations for the system of two bodies, the system of 12 stationary algebraic equations is decomposed using linear algebra methods and algorithms for Grobner basis construction. Depending on the parameters of the problem, the number of equilibria is found by analyzing the real roots of the algebraic equations from the Grobner basis constructed. Evolution of the conditions for equilibria existence in the dimensionless parameter space of the problem is investigated. The effectiveness of the algorithms for Grobner basis construction is analyzed depending on the number of parameters for the problem under consideration.

Performance and energy analysis of OpenMP runtime systems with dense linear algebra algorithms

In this article, we analyze performance and energy consumption of five OpenMP runtime systems over a non-uniform memory access (NUMA) platform. We also selected three CPU-level optimizations or techniques to evaluate their impact on the runtime systems: processors features Turbo Boost and C-States, and CPU Dynamic Voltage and Frequency Scaling through Linux CPUFreq governors. We present an experimental study to characterize OpenMP runtime systems on the three main kernels in dense linear algebra algorithms (Cholesky, LU, and QR) in terms of performance and energy consumption. Our experimental results suggest that OpenMP runtime systems can be considered as a new energy leverage, and Turbo Boost, as well as C-States, impacted significantly performance and energy. CPUFreq governors had more impact with Turbo Boost disabled, since both optimizations reduced performance due to CPU thermal limits. An LU factorization with concurrent-write extension from libKOMP achieved up to 63% of performance gain and 29% of energy decrease over original PLASMA algorithm using GNU C compiler (GCC) libGOMP runtime.

Untangling the nonlinearity in inverse scattering with data-driven reduced order models

The motivation of this work is an inverse problem for the acoustic wave equation, where an array of sensors probes an unknown medium with pulses and measures the scattered waves. The goal of the inversion is to determine from these measurements the structure of the scattering medium, modeled by a spatially varying acoustic impedance function. Many inversion algorithms assume that the mapping from the unknown impedance to the scattered waves is approximately linear. The linearization, known as the Born approximation, is not accurate in strongly scattering media, where the waves undergo multiple reflections before they reach the sensors in the array. Thus, the reconstructions of the impedance have numerous artifacts. The main result of the paper is a novel, linear-algebraic algorithm that uses a reduced order model (ROM) to map the data to those corresponding to the single scattering (Born) model. The ROM construction is based only on the measurements at the sensors in the array. The ROM is a proxy for the wave propagator operator, that propagates the wave in the unknown medium over the duration of the time sampling interval. The output of the algorithm can be input into any off-the-shelf inversion software that incorporates state of the art linear inversion algorithms to reconstruct the unknown acoustic impedance.

A sparse linear algebra algorithm for fast computation of prediction variances with Gaussian Markov random fields

Gaussian Markov random fields are used in a large number of disciplines in machine vision and spatial statistics. The models take advantage of sparsity in matrices introduced through the Markov assumptions, and all operations in inference and prediction use sparse linear algebra operations that scale well with dimensionality. Yet, for very high-dimensional models, exact computation of predictive variances of linear combinations of variables is generally computationally prohibitive, and approximate methods (generally interpolation or conditional simulation) are typically used instead. A set of conditions isestablished under which the variances of linear combinations of random variables can be computed exactly using the Takahashi recursions. The ensuing computational simplification has wide applicability and may be used to enhance several software packages where model fitting is seated in a maximum-likelihood framework. The resulting algorithm is ideal for use in a variety of spatial statistical applications, including LatticeKrig modelling, statistical downscaling, and fixed rank kriging. It can compute hundreds of thousands exact predictive variances of linear combinations on a standard desktop with ease, even when large spatial GMRF models are used.

An FGLM-like algorithm for computing the radical of a zero-dimensional ideal

We present a linear algebra algorithm, which, given a zero-dimensional ideal via a “good” representation, produces a separating linear form, its radical, again given via the same “good” representation and also a Kronecker parametrization of it. As an irrelevant byproduct, if you like, you can also read a Gröbner basis of the radical.

Estimation of DOA for Noncircular Signals via Vandermonde Constrained Parallel Factor Analysis

We provide a complete study on the direction-of-arrival (DOA) estimation of noncircular (NC) signals for uniform linear array (ULA) via Vandermonde constrained parallel factor (PARAFAC) analysis. By exploiting the noncircular property of the signals, we first construct an extended matrix which contains two times sampling number of the received signal. Then, taking the Vandermonde structure of the array manifold matrix into account, the extended matrix can be turned into a tensor model which admits the Vandermonde constrained PARAFAC decomposition. Based on this tensor model, an efficient linear algebra algorithm is applied to obtain the DOA estimation via utilizing the rotational invariance between two submatrices. Compared with some existing algorithms, the proposed method has a better DOA estimation performance. Meanwhile, the proposed method consistently has a higher estimation accuracy and a much lower computational complexity than the trilinear alternating least square- (TALS-) based PARAFAC algorithm. Finally, numerical examples are conducted to demonstrate the effectiveness of the proposed approach in terms of estimation accuracy and computational complexity.

Optimization of Triangular and Banded Matrix Operations Using 2d-Packed Layouts

Over the past few years, multicore systems have become increasingly powerful and thereby very useful in high-performance computing. However, many applications, such as some linear algebra algorithms, still cannot take full advantage of these systems. This is mainly due to the shortage of optimization techniques dealing with irregular control structures. In particular, the well-known polyhedral model fails to optimize loop nests whose bounds and/or array references are not affine functions. This is more likely to occur when handling sparse matrices in their packed formats. In this article, we propose using 2d-packed layouts and simple affine transformations to enable optimization of triangular and banded matrix operations. The benefit of our proposal is shown through an experimental study over a set of linear algebra benchmarks.

An efficient solution of real-time data processing for multi-GNSS network

Global navigation satellite systems (GNSS) are acting as an indispensable tool for geodetic research and global monitoring of the Earth, and they have been rapidly developed over the past few years with abundant GNSS networks, modern constellations, and significant improvement in mathematic models of data processing. However, due to the increasing number of satellites and stations, the computational efficiency becomes a key issue and it could hamper the further development of GNSS applications. In this contribution, this problem is overcome from the aspects of both dense linear algebra algorithms and GNSS processing strategy. First, in order to fully explore the power of modern microprocessors, the square root information filter solution based on the blocked QR factorization employing as many matrix–matrix operations as possible is introduced. In addition, the algorithm complexity of GNSS data processing is further decreased by centralizing the carrier-phase observations and ambiguity parameters, as well as performing the real-time ambiguity resolution and elimination. Based on the QR factorization of the simulated matrix, we can conclude that compared to unblocked QR factorization, the blocked QR factorization can greatly improve processing efficiency with a magnitude of nearly two orders on a personal computer with four 3.30 GHz cores. Then, with 82 globally distributed stations, the processing efficiency is further validated in multi-GNSS (GPS/BDS/Galileo) satellite clock estimation. The results suggest that it will take about 31.38 s per epoch for the unblocked method. While, without any loss of accuracy, it only takes 0.50 and 0.31 s for our new algorithm per epoch for float and fixed clock solutions, respectively.

Universality laws for randomized dimension reduction, with applications

AbstractDimension reduction is the process of embedding high-dimensional data into a lower dimensional space to facilitate its analysis. In the Euclidean setting, one fundamental technique for dimension reduction is to apply a random linear map to the data. This dimension reduction procedure succeeds when it preserves certain geometric features of the set. The question is how large the embedding dimension must be to ensure that randomized dimension reduction succeeds with high probability. This article studies a natural family of randomized dimension reduction maps and a large class of data sets. It proves that there is a phase transition in the success probability of the dimension reduction map as the embedding dimension increases. For a given data set, the location of the phase transition is the same for all maps in this family. Furthermore, each map has the same stability properties, as quantified through the minimum RSV. These results can be viewed as new universality laws in high-dimensional stochastic geometry. Universality laws for randomized dimension reduction have many applications in applied mathematics, signal processing and statistics. They yield design principles for numerical linear algebra algorithms, for compressed sensing measurement ensembles and for random linear codes. Furthermore, these results have implications for the performance of statistical estimation methods under a large class of random experimental designs.

Enhancing the scalability of selected inversion factorization algorithms in genomic prediction

A parallel distributed-memory approach for the exact calculation of selected entries of the inverse of a matrix arising in a Best Linear Unbiased Estimation (BLUE) problem in genomic prediction is presented. The particular structure of the matrices involved in this stochastic process, consisting of sparse and dense blocks, requires a framework coupling sparse and dense linear algebra algorithms. Our approach exploits direct sparse techniques based on the Takahashi equations, coupled with distributed LU dense factorizations and Schur-complement computations. The algorithm is validated on several matrices on a Cray XC40 supercomputer.

Adjusting data sparsity problem using linear algebra and machine learning algorithm

Data sparsity is one of the most important challenges in data in which each user only rates a small set of items. This problem is critical with increasing dimensions of data. We present an idea based on linear algebra and machine learning to solve this problem. This research applies a framework to cluster users and items in similar groups simultaneously. This method imputes appropriate values for missing data based on similar ratings in each cluster. This has the advantages of more accurate process results in each cluster due to users’ similarity of interests, and the reduction of sparsity negative effect. This approach is represented on 3dimensional data of users, items and times. The experimental results on MovieLense datasets, show that the method can help to overcome data sparsity, and increase the accuracy of prediction.

Randomized algorithms in numerical linear algebra

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 Matlab Radial Basis Function Toolbox

Radial Basis Function (RBF) methods are important tools for scattered data interpolation and for the solution of Partial Differential Equations in complexly shaped domains. The most straight forward approach used to evaluate the methods involves solving a linear system which is typically poorly conditioned. The Matlab Radial Basis Function toolbox features a regularization method for the ill-conditioned system, extended precision floating point arithmetic, and symmetry exploitation for the purpose of reducing flop counts of the associated numerical linear algebra algorithms.

Vandermonde Constrained Tensor Decomposition Based Blind Carrier Frequency Synchronization for OFDM Transmissions

To overcome the inter-carrier interference (ICI) of orthogonal frequency division multiplexing (OFDM) transmission subject to unknown carrier frequency offset and multipath channel, Vandermonde constrained tensor decomposition based blind approach is investigated in this paper. The main idea is that signal reception model is formulated as a tensor decomposition framework involving a Vandermonde structure factor matrix. Then an efficient linear algebra algorithm is resorted to the computation of the factors for ICI suppression and symbol recovery. Taking the Vandermonde structure into account in tensor model, not only the accuracy of the estimate is increased, but also the computational cost is reduced. Compared with the existing methods reported in the literatures, theoretical analysis and simulation results illustrate that the proposed algorithm provides lower computation complexity and better performance than ESPRIT and PARAFAC algorithms for OFDM system synchronous receiving.

Weighted SGD for ℓp regression with randomized preconditioning

In recent years, stochastic gradient descent (SGD) methods and randomized linear algebra (RLA) algorithms have been applied to many large-scale problems in machine learning and data analysis. SGD m...

Analog Computing in a Modern Context: A Linear Algebra Accelerator Case Study

Approaching the post-Moore's law era, researchers are looking for scalable ways to get useful computation from existing silicon technology. This article presents a programmable analog accelerator for solving systems of linear equations. The authors compensate for commonly perceived downsides of analog computing, such as low precision and accuracy, limited problem sizes, and difficulty applying it to different workloads. On the basis of a prototyped analog accelerator chip, they compare the analog solver's performance and energy consumption against an efficient digital algorithm running on a general-purpose processor. The analog accelerator approach is 10 times faster and provides 33 percent energy savings. Owing to the speed and efficiency of linear algebra algorithms running on digital computers, an analog accelerator that matches digital performance needs a large silicon footprint, which limits scalability. The authors conclude that problem classes outside of systems of linear equations could hold more promise for analog acceleration.