Matrix-vector Operations Research Articles

Acoustic 3D modeling by the method of integral equations

This paper presents a parallel algorithm for frequency-domain acoustic modeling by the method of integral equations (IE). The algorithm is applied to seismic simulation. The IE method reduces the size of the problem but leads to a dense system matrix. A tolerable memory consumption and numerical complexity were achieved by applying an iterative solver, accompanied by an effective matrix-vector multiplication operation, based on the fast Fourier transform (FFT). We demonstrate that, the IE system matrix is better conditioned than that of the finite-difference (FD) method, and discuss its relation to a specially preconditioned FD matrix. We considered several methods of matrix-vector multiplication for the free-space and layered host models. The developed algorithm and computer code were benchmarked against the FD time-domain solution. It was demonstrated that, the method could accurately calculate the seismic field for the models with sharp material boundaries and a point source and receiver located close to the free surface. We used OpenMP to speed up the matrix-vector multiplication, while MPI was used to speed up the solution of the system equations, and also for parallelizing across multiple sources. The practical examples and efficiency tests are presented as well.

Delayed Slater determinant update algorithms for high efficiency quantum Monte Carlo.

Within ab initio Quantum Monte Carlo simulations, the leading numerical cost for large systems is the computation of the values of the Slater determinants in the trial wavefunction. Each Monte Carlo step requires finding the determinant of a dense matrix. This is most commonly iteratively evaluated using a rank-1 Sherman-Morrison updating scheme to avoid repeated explicit calculation of the inverse. The overall computational cost is, therefore, formally cubic in the number of electrons or matrix size. To improve the numerical efficiency of this procedure, we propose a novel multiple rank delayed update scheme. This strategy enables probability evaluation with an application of accepted moves to the matrices delayed until after a predetermined number of moves, K. The accepted events are then applied to the matrices en bloc with enhanced arithmetic intensity and computational efficiency via matrix-matrix operations instead of matrix-vector operations. This procedure does not change the underlying Monte Carlo sampling or its statistical efficiency. For calculations on large systems and algorithms such as diffusion Monte Carlo, where the acceptance ratio is high, order of magnitude improvements in the update time can be obtained on both multi-core central processing units and graphical processing units.

FAST ALGORITHMS: SCIENCE, ART, CRAFT

State contains the results of a study of the specifics, possibilities and advantages of fast algorithms. This paper focuses on the description of the proposed approach for the development of fast algorithms using vector-matrix operations.

Scalable domain decomposition solvers for stochastic PDEs in high performance computing

Stochastic spectral finite element models of practical engineering systems may involve solutions of linear systems or linearized systems for non-linear problems with billions of unknowns. For stochastic modeling, it is therefore essential to design robust, parallel and scalable algorithms that can efficiently utilize high-performance computing to tackle such large-scale systems. Domain decomposition based iterative solvers can handle such systems. Although these algorithms exhibit excellent scalabilities, significant algorithmic and implementational challenges exist to extend them to solve extreme-scale stochastic systems using emerging computing platforms. Intrusive polynomial chaos expansion based domain decomposition algorithms are extended here to concurrently handle high resolution in both spatial and stochastic domains using an in-house implementation. Sparse iterative solvers with efficient preconditioners are employed to solve the resulting global and subdomain level local systems through multi-level iterative solvers. Parallel sparse matrix–vector operations are used to reduce the floating-point operations and memory requirements. Numerical and parallel scalabilities of these algorithms are presented for the diffusion equation having spatially varying diffusion coefficient modeled by a non-Gaussian stochastic process. Scalability of the solvers with respect to the number of random variables is also investigated.

A numerical investigation of matrix-free implicit time-stepping methods for large CFD simulations

This paper is concerned with development and testing of advanced time-stepping methods for large unsteady CFD problems in the method of lines approach, where the semi-discretization in space is performed first. The performance of several time discretization methods is studied numerically with regards to computational efficiency, order of accuracy, and stability, as well as the ability to effectively treat stiff problems. We consider matrix-free implementations, a popular approach for time-stepping methods applied to large CFD applications due to its adherence to scalable matrix-vector operations and a small memory footprint. We compare explicit methods with matrix-free implementations of implicit, linearly-implicit, as well as Rosenbrock–Krylov methods. We show that Rosenbrock–Krylov methods are competitive with existing techniques excelling for a number of problem types and settings.

Detecting buried wave-penetrable scatterers in a two-layered medium

In a two-layered medium, we prove that a buried inhomogeneous scatterer is uniquely determined from the wave field data measured in the upper half-space with respect to many incident point sources. Moreover, we extend the multilevel sampling method in Liu and Zou (2013) to numerically reconstruct the buried scatterer applying only few incident fields and partial scattered data. The extended recovery scheme only involves matrix–vector operations and does not need to solve any large-scale ill-posed linear systems or any optimization process. It is feasible to deal with the scatterers of different features and easy to implement, highly tolerant to noise and computationally quite cheap. We can regard it as an effective yet simple computational method to provide a reliable initial guess for the implementation in existing more accurate and refined optimization-type reconstruction algorithms.

Scaling Machine Learning via Compressed Linear Algebra

Large-scale machine learning (ML) algorithms are often iterative, using repeated read-only data access and I/Obound matrix-vector multiplications to converge to an optimal model. It is crucial for performance to fit the data into single-node or distributed main memory and enable very fast matrix-vector operations on in-memory data. Generalpurpose, heavy- and lightweight compression techniques struggle to achieve both good compression ratios and fast decompression speed to enable block-wise uncompressed operations. Compressed linear algebra (CLA) avoids these problems by applying lightweight lossless database compression techniques to matrices and then executing linear algebra operations such as matrix-vector multiplication directly on the compressed representations. The key ingredients are effective column compression schemes, cache-conscious operations, and an efficient sampling-based compression algorithm. Experiments on an initial implementation in SystemML show in-memory operations performance close to the uncompressed case and good compression ratios.We thereby obtain significant end-to-end performance improvements up to 26x or reduced memory requirements.

A new approach to integrate GPU-based Monte Carlo simulation into inverse treatment plan optimization for proton therapy

Monte Carlo (MC)-based spot dose calculation is highly desired for inverse treatment planning in proton therapy because of its accuracy. Recent studies on biological optimization have also indicated the use of MC methods to compute relevant quantities of interest, e.g. linear energy transfer. Although GPU-based MC engines have been developed to address inverse optimization problems, their efficiency still needs to be improved. Also, the use of a large number of GPUs in MC calculation is not favorable for clinical applications. The previously proposed adaptive particle sampling (APS) method can improve the efficiency of MC-based inverse optimization by using the computationally expensive MC simulation more effectively. This method is more efficient than the conventional approach that performs spot dose calculation and optimization in two sequential steps. In this paper, we propose a computational library to perform MC-based spot dose calculation on GPU with the APS scheme. The implemented APS method performs a non-uniform sampling of the particles from pencil beam spots during the optimization process, favoring those from the high intensity spots. The library also conducts two computationally intensive matrix-vector operations frequently used when solving an optimization problem. This library design allows a streamlined integration of the MC-based spot dose calculation into an existing proton therapy inverse planning process. We tested the developed library in a typical inverse optimization system with four patient cases. The library achieved the targeted functions by supporting inverse planning in various proton therapy schemes, e.g. single field uniform dose, 3D intensity modulated proton therapy, and distal edge tracking. The efficiency was 41.6 ± 15.3% higher than the use of a GPU-based MC package in a conventional calculation scheme. The total computation time ranged between 2 and 50 min on a single GPU card depending on the problem size.

Parallel algorithm for 3D modeling of monochromatic acoustic field using integral equations

We present a parallel algorithm for solution of the three-dimensional Helmholtz equation in the frequency domain by the method of volume integral equations. The algorithm is applied to seismic forward modeling. The method of integral equations reduces the size of the problem by dividing the geologic model into the anomalous and background parts, but leads to a dense system matrix. A tolerable memory consumption and numerical complexity were achieved by applying an iterative solver, accompanied by an effective matrix-vector multiplication operation, based on the fast Fourier transform. We used OpenMP to speed up the matrix-vector multiplication, while MPI was used to speed up the equation system solver, and also for parallelizing across multiple sources. Practical examples and efficiency tests are presented.

Reducing latency cost in 2D sparse matrix partitioning models

Abstract Sparse matrix partitioning is a common technique used for improving performance of parallel linear iterative solvers. Compared to solvers used for symmetric linear systems, solvers for nonsymmetric systems offer more potential for addressing different multiple communication metrics due to the flexibility of adopting different partitions on the input and output vectors of sparse matrix-vector multiplication operations. In this regard, there exist works based on one-dimensional (1D) and two-dimensional (2D) fine-grain partitioning models that effectively address both bandwidth and latency costs in nonsymmetric solvers. In this work, we propose two new models based on 2D checkerboard and jagged partitioning. These models aim at minimizing total message count while maintaining a balance on communication volume loads of processors; hence, they address both bandwidth and latency costs. We evaluate all partitioning models on two nonsymmetric system solvers implemented using the widely adopted PETSc toolkit and conduct extensive experiments using these solvers on a modern system (a BlueGene/Q machine) successfully scaling them up to 8K processors. Along with the proposed models, we put practical aspects of eight evaluated models (two 1D- and six 2D-based) under thorough analysis. To the best of our knowledge, this is the first work that analyzes practical performance of 2D models on this scale. Among evaluated models, the models that rely on 2D jagged partitioning obtain the most promising results by striking a balance between minimizing bandwidth and latency costs.

Fast Hadamard transforms for compressive sensing of joint systems: measurement of a 3.2 million-dimensional bi-photon probability distribution.

We demonstrate how to efficiently implement extremely high-dimensional compressive imaging of a bi-photon probability distribution. Our method uses fast-Hadamard-transform Kronecker-based compressive sensing to acquire the joint space distribution. We list, in detail, the operations necessary to enable fast-transform-based matrix-vector operations in the joint space to reconstruct a 16.8 million-dimensional image in less than 10 minutes. Within a subspace of that image exists a 3.2 million-dimensional bi-photon probability distribution. In addition, we demonstrate how the marginal distributions can aid in the accuracy of joint space distribution reconstructions.

A FAST AND ROBUST NUMERICAL METHOD FOR OPTION PRICES AND GREEKS IN A JUMP-DIFFUSION MODEL

We propose a fast and robust finite difference method for Merton``s jump diffusion model, which is a partial integro-differential equation. To speed up a computational time, we compute a matrix so that we can calculate the non-local integral term fast by a simple matrix-vector operation. Also, we use non-uniform grids to increase efficiency. We present numerical experiments such as evaluation of the option prices and Greeks to demonstrate a performance of the proposed numerical method. The computational results are in good agreements with the exact solutions of the jump-diffusion model.

Compression approaches for the regularized solutions of linear systems from large-scale inverse problems

We introduce and compare new compression approaches to obtain regularized solutions of large linear systems which are commonly encountered in large scale inverse problems. We first describe how to approximate matrix vector operations with a large matrix through a sparser matrix with fewer nonzero elements, by borrowing from ideas used in wavelet image compression. Next, we describe and compare approaches based on the use of the low rank singular value decomposition (SVD), which can result in further size reductions. We describe how to obtain the approximate low rank SVD of the original matrix using the sparser wavelet compressed matrix. Some analytical results concerning the various methods are presented and the results of the proposed techniques are illustrated using both synthetic data and a very large linear system from a seismic tomography application, where we obtain significant compression gains with our methods, while still resolving the main features of the solutions.

A Simple Observer for Gyro and Accelerometer Biases in Land Navigation Systems

In various applications of land vehicle navigation and automatic guidance systems, Global Navigation Satellite System/Inertial Measurement Unit (GNSS/IMU) positioning performance crucially depends on the attitude determination accuracy affected by gyro and accelerometer bias instabilities. Traditional bias estimation approaches based on the Kalman filter suffer from implementation complexity and require non-intuitive tuning procedures. In this paper we propose, as an alternative, a simple observer that estimates inertial sensor biases exclusively in terms of quantities with obvious geometrical meaning. By this, any multidimensional vector-matrix operations are avoided and observer tuning is substantially simplified. The observer has been successfully tested in a farming vehicle navigation system.

Increasing the Performance of the Jacobi--Davidson Method by Blocking

Block variants of the Jacobi--Davidson method for computing a few eigenpairs of a large sparse matrix are known to improve the robustness of the standard algorithm when it comes to computing multiple or clustered eigenvalues. In practice, however, they are typically avoided because the total number of matrix-vector operations increases. In this paper we present the implementation of a block Jacobi--Davidson solver. By detailed performance engineering and numerical experiments we demonstrate that the increase in operations is typically more than compensated by performance gains through better cache usage on modern CPUs, resulting in a method that is both more efficient and robust than its single vector counterpart. The steps to be taken to achieve a block speedup involve both kernel optimizations for sparse matrix and block vector operations, and algorithmic choices to allow using blocked operations in most parts of the computation. We discuss the aspect of avoiding synchronization in the algorithm and show by numerical experiments with our hybrid parallel implementation that a significant speedup through blocking can be achieved for a variety of matrices on up to 5120 CPU cores as long as at least about 20 eigenpairs are sought.

A Compact Numerical Implementation for Solving Stokes Equations Using Matrix-vector Operations

Abstract In this work, a numerical scheme is implemented to solve Stokes equations based on cell-centered finite difference over staggered grid. In this scheme, all the difference operations have been vectorized thereby eliminating loops. This is particularly important when using programming languages that require interpretations, e.g., MATLAB and Python. Using this scheme, the execution time becomes significantly smaller compared with non-vectorized operations and also become comparable with those languages that require no repeated interpretations like FORTRAN, C, etc. This technique has also been applied to Navier-Stokes equations under laminar flow conditions.

Cloud data integrity checking protocol from lattice

Cloud computing is a distributed computing model, which provides services and resources. As a new type of storage services, cloud storage provides a new method for mass data storage management. Nowadays, more and more users upload their data to cloud storage servers and retrieve data when needed. One of the basic tasks of cloud storage is to guarantee the integrity between retrieved and uploaded data. Therefore, data integrity checking of cloud storage becomes a key issue. This paper proposes a privacy-preserving public verifying cloud data integrity checking protocol. To achieve security against quantum computer attacks, we resort to constructing homomorphic linear authenticator based on a lattice-based homomorphic signature algorithm. To eliminate data content exposure to achieve privacy, we use the blind checking technique. Theoretical analysis indicates that our scheme achieves security against malicious cloud server and privacy against third party auditor TPA, and resists the known-proof forgery attack. Compared with other protocols, our protocol does not resort to rank-based authenticated skip list, saving storage and communication overhead. Furthermore, it makes use of matrix-vector multiplication operation, avoiding great amount of module exponential operation, which reduces the computational labour for cloud server and TPA.

SlashBurn: Graph Compression and Mining beyond Caveman Communities

Given a real world graph, how should we lay-out its edges? How can we compress it? These questions are closely related, and the typical approach so far is to find clique-like communities, like the cavemen graph', and compress them. We show that the block-diagonal mental image of the cavemen graph' is the wrong paradigm, in full agreement with earlier results that real world graphs have no good cuts. Instead, we propose to envision graphs as a collection of hubs connecting spokes, with super-hubs connecting the hubs, and so on, recursively. Based on the idea, we propose the SLASHBURN method to recursively split a graph into hubs and spokes connected only by the hubs. We also propose techniques to select the hubs and give an ordering to the spokes, in addition to the basic SLASHBURN. We give theoretical analysis of the proposed hub selection methods. Our view point has several advantages: (a) it avoids the no good cuts' problem, (b) it gives better compression, and (c) it leads to faster execution times for matrix-vector operations, which are the back-bone of most graph processing tools. Through experiments, we show that SLASHBURN consistently outperforms other methods for all data sets, resulting in better compression and faster running time. Moreover, we show that SLASHBURN with the appropriate spokes ordering can further improve compression while hardly sacrificing the running time.

Distributed power-line outage detection based on wide area measurement system.

In modern power grids, the fast and reliable detection of power-line outages is an important functionality, which prevents cascading failures and facilitates an accurate state estimation to monitor the real-time conditions of the grids. However, most of the existing approaches for outage detection suffer from two drawbacks, namely: (i) high computational complexity; and (ii) relying on a centralized means of implementation. The high computational complexity limits the practical usage of outage detection only for the case of single-line or double-line outages. Meanwhile, the centralized means of implementation raises security and privacy issues. Considering these drawbacks, the present paper proposes a distributed framework, which carries out in-network information processing and only shares estimates on boundaries with the neighboring control areas. This novel framework relies on a convex-relaxed formulation of the line outage detection problem and leverages the alternating direction method of multipliers (ADMM) for its distributed solution. The proposed framework invokes a low computational complexity, requiring only linear and simple matrix-vector operations. We also extend this framework to incorporate the sparse property of the measurement matrix and employ the LSQRalgorithm to enable a warm start, which further accelerates the algorithm. Analysis and simulation tests validate the correctness and effectiveness of the proposed approaches.

On Application of Matrix-Vector Operations for Solving Constrained Optimization Problems

This paper discusses the application of matrix-vector operations technique for solving constrained optimization problems. The method is aimed at circumventing the computational rigours undergone using the simplex and revised simplex method in solving this class of problems with aim of reducing the computer memory space occupied by the methods. In other to achieve this, a straightforward and simple to handle matrix-vector operations algorithm has been developed to solve the same constrained problems. Numerical results show some improvements compared with the classical method.