Basic Linear Algebra Operations Research Articles

Evaluation of the efficiency of parallel algorithms using intel parallel studio

The paper analyzes the means of evaluating the effectiveness of parallel programs. The research focuses on im- proving the performance of basic linear algebra operations performed on parallel computing systems with shared memory. The Intel Parallel Studio environment is used, which provides a powerful set of tools and resources for analyzing and optimizing parallel programs. The results of the study make it possible to solve performance-related problems and improve the operation of parallel programs.

The SHARK integral generation and digestion system.

In this paper, the SHARK integral generation and digestion engine is described. In essence, SHARK is based on a reformulation of the popular McMurchie/Davidson approach to molecular integrals. This reformulation leads to an efficient algorithm that is driven by BLAS level 3 operations. The algorithm is particularly efficient for high angular momentum basis functions (up to L=7 is available by default, but the algorithm is programmed for arbitrary angular momenta). SHARK features a significant number of specific programming constructs that are designed to greatly simplify the workflow in quantum chemical program development and avoid undesirable code duplication to the largest possible extent. SHARK can handle segmented, generally and partially generally contracted basis sets. It can be used to generate a host of one- and two-electron integrals over various kernels including, two-, three-, and four-index repulsion integrals, integrals over Gauge Including Atomic Orbitals (GIAOs), relativistic integrals and integrals featuring a finite nucleus model. SHARK provides routines to evaluate Fock like matrices, generate integral transformations and related tasks. SHARK is the essential engine inside the ORCA package that drives essentially all tasks that are related to integrals over basis functions in version ORCA 5.0 and higher. Since the core of SHARK is based on low-level basic linear algebra (BLAS) operations, it is expected to not only perform well on present day but also on future hardware provided that the hardware manufacturer provides a properly optimized BLAS library for matrix and vector operations. Representative timings and comparisons to the Libint library used by ORCA are reported for Intel i9 and Apple M1 max processors.

Construction and application of provable positive and exact cubature formulas

Abstract Many applications require multi-dimensional numerical integration, often in the form of a cubature formula (CF). These CFs are desired to be positive and exact for certain finite-dimensional function spaces (and weight functions). Although there are several efficient procedures to construct positive and exact CFs for many standard cases, it remains a challenge to do so in a more general setting. Here, we show how the method of least squares (LSs) can be used to derive provable positive and exact formulas in a general multi-dimensional setting. Thereby, the procedure only makes use of basic linear algebra operations, such as solving an LSs problem. In particular, it is proved that the resulting LSs CFs are ensured to be positive and exact if a sufficiently large number of equidistributed data points is used. We also discuss the application of provable positive and exact LSs CFs to construct nested stable high-order rules and positive interpolatory formulas. Finally, our findings shed new light on some existing methods for multi-variate numerical integration and under which restrictions these are ensured to be successful.

異種混合プロセッサにおける並列FEMコードの性能評価

The main computations in finite element analysis are basic linear algebra operations on matrices and vectors, and even load distribution is suitable for these parallel computations. Therefore, parallel FEM codes implemented with OpenMP and simple scheduling often achieve good efficiency in conventional multi-core processors. However, in recent years, desktop CPUs such as Intel's Alder Lake and Apple's M1 have adopted heterogeneous architectures with high-efficiency and high-performance cores. Therefore, it is unclear whether the conventional parallel FEM code can extract the processor's performance on those parallel computers. This research evaluates the performance of a heterogeneous processor using a parallel FEM code that uses the conjugate gradient method as a linear solver and discusses the results.

A majorization-minimization scheme for L 2 support vector regression

In a support vector regression (SVR) model, using the squared ϵ-insensitive loss function makes the optimization problem strictly convex and yields a more concise solution. However, the formulation of -SVR leads to a quadratic programming which is expensive to solve. This paper reformulates the optimization problem of -SVR by absorbing the constraints in the objective function, which can be solved efficiently by a majorization-minimization approach, in which an upper bound for the objective function is derived in each iteration which is easier to be minimized. The proposed approach is easy to implement, without requiring any additional computing package other than basic linear algebra operations. Numerical studies on real-world datasets show that, compared to the alternatives, the proposed approach can achieve similar prediction accuracy with substantially higher time efficiency in training.

Constructing the Singular Roesser State-Space Model Description of 3D Spatio-Temporal Dynamics From the Polynomial System Matrix

This paper considers systems theoretic properties of linear systems defined in terms of spatial and temporal indeterminates. These include physical applications where one of the indeterminates is of finite duration. In some cases, a singular Roesser state-space model representation of the dynamics has found use in characterizing systems theoretic properties. The representation of the dynamics of many linear systems is obtained in terms of transform variables and a polynomial system matrix representation. This paper develops a direct method for constructing the singular Roesser state-space realization from the system matrix description for 3D systems such that relevant properties are retained. Since this method developed relies on basic linear algebra operations, it may be highly effective from the computational standpoint. In particular, spatially interconnected systems of the form of the ladder circuits are considered as the example. This application confirms the usefulness and effectiveness of the proposed method independently of the system spatial order.

Laplacian-optimized diffusion for semi-supervised learning

Semi-supervised learning (SSL) is fundamentally a geometric task: in order to classify high-dimensional point sets when only a small fraction of data points are labeled, the geometry of the unlabeled data points is exploited to gain better classifying accuracy. A number of state-of-the-art SSL techniques rely on label propagation through graph-based diffusion, with edge weights that are evaluated either analytically from the data or through compute-intensive training based on nonlinear and nonconvex optimization. In this paper, we bring discrete differential geometry to bear on this problem by introducing a graph-based SSL approach where label diffusion uses a Laplacian operator learned from the geometry of the input data. From a data-dependent graph of the input, we formulate a biconvex loss function in terms of graph edge weights and inferred labels. Its minimization is achieved through alternating rounds of optimization of the Laplacian and diffusion-based inference of labels. The resulting optimized Laplacian diffusion directionally adapts to the intrinsic geometric structure of the data which often concentrates in clusters or around low-dimensional manifolds within the high-dimensional representation space. We show on a range of classical datasets that our variational classification is more accurate than current graph-based SSL techniques. The algorithmic simplicity and efficiency of our discrete differential geometric approach (limited to basic linear algebra operations) also make it attractive, despite the seemingly complex task of optimizing all the edge weights of a graph.

A Nonparametric Approach to Bond Portfolio Immunization

We consider the problem of short term immunization of a bond-like obligation with respect to changes in interest rates using a portfolio of bonds. In the case that the zero-coupon yield curve belongs to a fixed low-dimensional manifold, the problem is widely known as parametric immunization. Parametric immunization seeks to make the sensitivities of the hedged portfolio price with respect to all model parameters equal to zero. However, within a popular approach of nonparametric (smoothing spline) term structure estimation, parametric hedging is not applicable right away. We present a nonparametric approach to hedging a bond-like obligation allowing for a general form of the term structure estimator with possible smoothing. We show that our approach yields the standard duration based immunization in the limit when the amount of smoothing goes to infinity. We also recover the industry best practice approach of hedging based on key rate durations as another particular case. The hedging portfolio is straightforward to calculate using only basic linear algebra operations.

Iteratively reweighted least square for asymmetric L 2-Loss support vector regression

In support vector regression (SVR) model, using the squared -insensitive loss function makes the objective function of the optimization problem strictly convex and yields a more concise solution. However, the formulation leads to a quadratic programing which is expensive to solve. This paper reformulates the optimization problem by absorbing the constraints in the objective function, and the new formulation shares similarity with weighted least square regression problem. Based on this formulation, we propose an iteratively reweighted least square approach to train the L 2-loss SVR, for both linear and nonlinear models. The proposed approach is easy to implement, without requiring any additional computing package other than basic linear algebra operations. Numerical studies on real-world datasets show that, compared to the alternatives, the proposed approach can achieve similar prediction accuracy with substantially higher time efficiency.

Exact Smooth Term-Structure Estimation

We present a nonparametric method to estimate the discount curve from market quotes based on the Moore--Penrose pseudoinverse. The discount curve reproduces the market quotes perfectly, has maximal smoothness, and is given in closed-form. The method is easy to implement and requires only basic linear algebra operations. We provide a full theoretical framework as well as several practical applications.

A nearest-neighbour discretisation of the regularized stokeslet boundary integral equation

The method of regularized stokeslets is extensively used in biological fluid dynamics due to its conceptual simplicity and meshlessness. This simplicity carries a degree of cost in computational expense and accuracy because the number of degrees of freedom used to discretise the unknown surface traction is generally significantly higher than that required by boundary element methods. We describe a meshless method based on nearest-neighbour interpolation that significantly reduces the number of degrees of freedom required to discretise the unknown traction, increasing the range of problems that can be practically solved, without excessively complicating the task of the modeller. The nearest-neighbour technique is tested against the classical problem of rigid body motion of a sphere immersed in very viscous fluid, then applied to the more complex biophysical problem of calculating the rotational diffusion timescales of a macromolecular structure modelled by three closely-spaced non-slender rods. A heuristic for finding the required density of force and quadrature points by numerical refinement is suggested. Matlab/GNU Octave code for the key steps of the algorithm is provided, which predominantly use basic linear algebra operations, with a full implementation being provided on github. Compared with the standard Nyström discretisation, more accurate and substantially more efficient results can be obtained by de-refining the force discretisation relative to the quadrature discretisation: a cost reduction of over 10 times with improved accuracy is observed. This improvement comes at minimal additional technical complexity. Future avenues to develop the algorithm are then discussed.

Open Cezeri Library: A novel java based matrix and computer vision framework

ABSTRACTIn this paper we introduce the Open Cezeri Library (OCL) framework as a domain specific language (DSL) for researchers, scientists, and engineering students to enable them to develop basic linear algebra operations via simple matrix calculations, image processing, computer vision, and machine learning applications in JAVA programming language. OCL provides a strong intuition of coding for the developer while implementing by means of a fluent interface. The significant aspect of the OCL is to combine the methods of well‐known platforms; MATLAB and JAVA, accordingly. Moreover, OCL supports a fluent interface so that users can extend a single line of codes by putting a dot between the methods because all the methods implemented actually return the host class. It was observed that the learning curve of the OCL is lower than the MATLAB and the native JAVA languages, and makes coding more readable, understandable, traceable, and enjoyable. In addition to this, the experiments revealed that the running performance of the OCL is quite comparable and can be used in a variety of diverse applications. © 2016 Wiley Periodicals, Inc. Comput Appl Eng Educ 24:736–743, 2016; View this article online at wileyonlinelibrary.com/journal/cae; DOI 10.1002/cae.21745

GPU performance analysis of a nodal discontinuous Galerkin method for acoustic and elastic models

Finite element schemes based on discontinuous Galerkin methods possess features amenable to massively parallel computing accelerated with general purpose graphics processing units (GPUs). However, the computational performance of such schemes strongly depends on their implementation. In the past, several implementation strategies have been proposed. They are based exclusively on specialized compute kernels tuned for each operation, or they can leverage BLAS libraries that provide optimized routines for basic linear algebra operations. In this paper, we present and analyze up-to-date performance results for different implementations, tested in a unified framework on a single NVIDIA GTX980 GPU. We show that specialized kernels written with a one-node-per-thread strategy are competitive for polynomial bases up to the fifth and seventh degrees for acoustic and elastic models, respectively. For higher degrees, a strategy that makes use of the NVIDIA cuBLAS library provides better results, able to reach a net arithmetic throughput 35.7% of the theoretical peak value.

MPI-CUDA sparse matrix–vector multiplication for the conjugate gradient method with an approximate inverse preconditioner

The preconditioned conjugate gradient (PCG) is one of the most prominent iterative methods for the solution of sparse linear systems with symmetric and positive definite matrix that arise, for example, in the modeling of incompressible flows. The method relies on a set of basic linear algebra operations which determine the overall performance. Therefore, to achieve improvements in the performance, implementations of these basic operations must be adapted to the changes in the architecture of parallel computing systems. In the last years, one of the strategies to increase the computing power of supercomputers has been the usage of Graphics Processing Units (GPUs) as math co-processors in addition to CPUs. This paper presents a MPI-CUDA implementation of the PCG solver for such hybrid computing systems composed of multiple CPUs and GPUs. Special attention has been paid to the sparse matrix–vector multiplication (SpMV), because most of the execution time of the solver is spent on this operation. The approximate inverse preconditioner, which is used to improve the convergence of the CG solver, is also based on the SpMV operation. An overlapping of data transfer and computations is proposed in order to hide the MPI and the CPU-GPU communications needed to perform parallel SpMVs. This strategy has shown a considerable improvement and, as a result, the hybrid implementation of the PCG solver has demonstrated a significant speedup compared to the CPU-only implementation.

Use of General Purpose Graphics Processing Units with MODFLOW

To evaluate the use of general-purpose graphics processing units (GPGPUs) to improve the performance of MODFLOW, an unstructured preconditioned conjugate gradient (UPCG) solver has been developed. The UPCG solver uses a compressed sparse row storage scheme and includes Jacobi, zero fill-in incomplete, and modified-incomplete lower-upper (LU) factorization, and generalized least-squares polynomial preconditioners. The UPCG solver also includes options for sequential and parallel solution on the central processing unit (CPU) using OpenMP. For simulations utilizing the GPGPU, all basic linear algebra operations are performed on the GPGPU; memory copies between the central processing unit CPU and GPCPU occur prior to the first iteration of the UPCG solver and after satisfying head and flow criteria or exceeding a maximum number of iterations. The efficiency of the UPCG solver for GPGPU and CPU solutions is benchmarked using simulations of a synthetic, heterogeneous unconfined aquifer with tens of thousands to millions of active grid cells. Testing indicates GPGPU speedups on the order of 2 to 8, relative to the standard MODFLOW preconditioned conjugate gradient (PCG) solver, can be achieved when (1) memory copies between the CPU and GPGPU are optimized, (2) the percentage of time performing memory copies between the CPU and GPGPU is small relative to the calculation time, (3) high-performance GPGPU cards are utilized, and (4) CPU-GPGPU combinations are used to execute sequential operations that are difficult to parallelize. Furthermore, UPCG solver testing indicates GPGPU speedups exceed parallel CPU speedups achieved using OpenMP on multicore CPUs for preconditioners that can be easily parallelized.

GPU‐accelerated computation of electron transfer

Electron transfer is a fundamental process that can be studied with the help of computer simulation. The underlying quantum mechanical description renders the problem a computationally intensive application. In this study, we probe the graphics processing unit (GPU) for suitability to this type of problem. Time-critical components are identified via profiling of an existing implementation and several different variants are tested involving the GPU at increasing levels of abstraction. A publicly available library supporting basic linear algebra operations on the GPU turns out to accelerate the computation approximately 50-fold with minor dependence on actual problem size. The performance gain does not compromise numerical accuracy and is of significant value for practical purposes.

Solution of Linear Systems and Matrix Inversion in the TT-Format

Tensors arise naturally in high-dimensional problems in chemistry, financial mathematics, and many other areas. The numerical treatment of such problems is difficult due to the curse of dimensionality: the number of unknowns and the computational complexity grow exponentially with the dimension of the problem. To break the curse of dimensionality, low-parametric representations, or formats, have to be used. In this paper we make use of the TT-format (tensor-train format) which is one of the most effective stable representations of high-dimensional tensors. Basic linear algebra operations in the TT-format are now well developed. Our goal is to provide a “black-box” type of solver for linear systems where both the matrix and the right-hand side are in the TT-format. An efficient DMRG (density matrix renormalization group) method is proposed, and several tricks are employed to make it work. The numerical experiments confirm the effectiveness of our approach.

Tensor-Train Decomposition

A simple nonrecursive form of the tensor decomposition in d dimensions is presented. It does not inherently suffer from the curse of dimensionality, it has asymptotically the same number of parameters as the canonical decomposition, but it is stable and its computation is based on low-rank approximation of auxiliary unfolding matrices. The new form gives a clear and convenient way to implement all basic operations efficiently. A fast rounding procedure is presented, as well as basic linear algebra operations. Examples showing the benefits of the decomposition are given, and the efficiency is demonstrated by the computation of the smallest eigenvalue of a 19-dimensional operator.

A note on the O(n)-storage implementation of the GKO algorithm and its adaptation to Trummer-like matrices

We propose a new O(n)-space implementation of the GKO-Cauchy algorithm for the solution of linear systems with Cauchy-like matrix. Despite its slightly higher computational cost, this new algorithm makes a more efficient use of the processor cache memory. Thus, for matrices of size larger than about 500-1000, it outperforms the existing algorithms. We present an applicative case of Cauchy-like matrices with non-reconstructible main diagonal. In this special instance, the O(n) space algorithms can be adapted nicely to provide an efficient implementation of basic linear algebra operations in terms of the low displacement-rank generators.

Linear algebra for tensor problems

By a tensor problem in general, we mean one where all the data on input and output are given (exactly or approximately) in tensor formats, the number of data representation parameters being much smaller than the total amount of data. For such problems, it is natural to seek for algorithms working with data only in tensor formats maintaining the same small number of representation parameters—by the price of all results of computation to be contaminated by approximation (recompression) to occur in each operation. Since approximation time is crucial and depends on tensor formats in use, in this paper we discuss which are best suitable to make recompression inexpensive and reliable. We present fast recompression procedures with sublinear complexity with respect to the size of data and propose methods for basic linear algebra operations with all matrix operands in the Tucker format, mostly through calls to highly optimized level-3 BLAS/LAPACK routines. We show that for three-dimensional tensors the canonical format can be avoided without any loss of efficiency. Numerical illustrations are given for approximate matrix inversion via proposed recompression techniques.