Fast Matrix-vector Multiplication Research Articles

PF‐GEMV: Utilization maximizing architecture in fast matrix–vector multiplication for GPT‐2 inference

AbstractOwing to the widespread advancement of transformer‐based artificial neural networks, artificial intelligence (AI) processors are now required to perform matrix–vector multiplication in addition to the conventional matrix–matrix multiplication. However, current AI processor architectures are optimized for general matrix–matrix multiplications (GEMMs), which causes significant throughput degradation when processing general matrix–vector multiplications (GEMVs). In this study, we proposed a port‐folding GEMV (PF‐GEMV) scheme employing multiformat and low‐precision techniques while reusing an outer product‐based processor optimized for conventional GEMM operations. This approach achieves 93.7% utilization in GEMV operations with an 8‐bit format on an 8 8 processor, thus resulting in a 7.5 increase in throughput compared with that of the original scheme. Furthermore, when applied to the matrix operation of the GPT‐2 large model, an increase in speed by 7 is achieved in single‐batch inferences.

Exact novel formulas and fast algorithm of potential for a hammock resistor network

The establishment of a resistor network model has become a sharp edge to solve complex scientific problems. In this paper, we introduce Chebyshev polynomials to express the potential formula of the hammock resistor network and improve the general solution of the hammock resistor network. Moreover, through the change in different parameters, special potential formulas are proposed and displayed in 3D dynamic view. A fast algorithm of the calculating potential is given by using the matrix equation model, discrete cosine transform-II, and the fast matrix-vector multiplication. Finally, we show the advantages of our improved potential formula and fast algorithm by the calculation efficiency of the three methods. The modified potential formula and the presented fast algorithm provide a new tool for the field of science and engineering.

Memristor-Based Spectral Decomposition of Matrices and Its Applications

The recently developed memristor technology allows for extremely fast implementation of a number of important matrix operations and algorithms. Moreover, the existence of fast matrix-vector operations offers the opportunity to design new matrix algorithms that exploit these operations. Here, we focus on the spectral decomposition of matrices, a task that plays an important role in a wide variety of applications from different engineering and scientific fields, including network science, control theory, advanced dynamics, and quantum mechanics. While there are a number of algorithms designed to find eigenvalues and eigenvectors of a matrix, these methods often suffer from poor running time performance. In this work, we present an algorithm for finding eigenvalues and eigenvectors that is designed to be used on memristor crossbar arrays. Although this algorithm can be implemented in a non-memristive system, its fast running time relies on the availability of extremely fast matrix-vector multiplication, as is offered by a memristor crossbar array. In this paper, we (1) show the running time improvements of existing eigendecomposition algorithms when matrix-vector multiplications are performed on a memristor crossbar array, and (2) present <monospace>EigSweep</monospace>, a novel, fully-parallel, fast and flexible eigendecomposition algorithm that gives an improvement in running time over traditional eigendecomposition algorithms when all are accelerated by a memristor crossbar. We discuss algorithmic aspects as well as hardware-related aspects of the implementation of <monospace>EigSweep</monospace>, and perform an extensive experimental analysis on real-world and synthetic matrices.

Analytical potential formulae and fast algorithm for a horn torus resistor network.

In this paper, a (u+1)×v horn torus resistor network with a special boundary is researched. According to Kirchhoff's law and the recursion-transform method, a model of the resistor network is established by the voltage V and a perturbed tridiagonal Toeplitz matrix. We obtain the exact potential formula of a horn torus resistor network. First, the orthogonal matrix transformation is constructed to obtain the eigenvalues and eigenvectors of this perturbed tridiagonal Toeplitz matrix; second, the solution of the node voltage is given by using the famous fifth kind of discrete sine transform (DST-V). We introduce Chebyshev polynomials to represent the exact potential formula. In addition, the equivalent resistance formulae in special cases are given and displayed by a three-dimensional dynamic view. Finally, a fast algorithm of computing potential is proposed by using the mathematical model, famous DST-V, and fast matrix-vector multiplication. The exact potential formula and the proposed fast algorithm realize large-scale fast and efficient operation for a (u+1)×v horn torus resistor network, respectively.

Efficient Identification of Butterfly Sparse Matrix Factorizations

Fast transforms correspond to factorizations of the form , where each factor is sparse and possibly structured. This paper investigates essential uniqueness of such factorizations, i.e., uniqueness up to unavoidable scaling ambiguities. Our main contribution is to prove that any matrix having the so-called butterfly structure admits an essentially unique factorization into butterfly factors (where ), and that the factors can be recovered by a hierarchical factorization method, which consists in recursively factorizing the considered matrix into two factors. This hierarchical identifiability property relies on a simple identifiability condition in the two-layer and fixed-support setting. This approach contrasts with existing ones that fit the product of butterfly factors to a given matrix via gradient descent. The proposed method can be applied in particular to retrieve the factorization of the Hadamard or the discrete Fourier transform matrices of size . Computing such factorizations costs , which is on the order of dense matrix-vector multiplication, while the obtained factorizations enable fast matrix-vector multiplications and have the potential to be applied to compress deep neural networks.

Algorithms for solving a class of real quasi-symmetric Toeplitz linear systems and its applications

&lt;abstract&gt;&lt;p&gt;In this paper, fast numerical methods for solving the real quasi-symmetric Toeplitz linear system are studied in two stages. First, based on an order-reduction algorithm and the factorization of Toeplitz matrix inversion, a sequence of linear systems with a constant symmetric Toeplitz matrix are solved. Second, two new fast algorithms are employed to solve the real quasi-symmetric Toeplitz linear system. Furthermore, we show a fast algorithm for quasi-symmetric Toeplitz matrix-vector multiplication. In addition, the stability analysis of the splitting symmetric Toeplitz inversion is discussed. In mathematical or engineering problems, the proposed algorithms are extraordinarily effective for solving a sequence of linear systems with a constant symmetric Toeplitz matrix. Fast matrix-vector multiplication and a quasi-symmetric Toeplitz linear solver are proven to be suitable for image encryption and decryption.&lt;/p&gt;&lt;/abstract&gt;

An extra-component method for evaluating fast matrix-vector multiplication with special functions

In calculating integral or discrete transforms, use has been made of fast algorithms for multiplying vectors by matrices whose elements are specified as values of special (Chebyshev, Legendre, Jacobi, Laguerre, etc.) functions. To achieve higher efficiency, a convenient general approach for calculating matrix-vector products for some class of problems is proposed. A series of fast simple-structure algorithms developed under this approach can be efficiently implemented with software based on modern microprocessors. Computational experiments show that the method has a precomputation complexity of $$O(N^2 \log N)$$ and an execution complexity of $$O(N \log N)$$ in general. This agrees with the analytical estimates obtained for trigonometric functions. In addition, for trigonometric functions, the number of arithmetic operations at the preliminary stage can be decreased to $$O(N \log N)$$ . The results of computational experiments with the algorithms show that these procedures can decrease the calculation time by several orders of magnitude compared with a conventional direct method of matrix-vector multiplication.

Grouped Transformations and Regularization in High-Dimensional Explainable ANOVA Approximation

In this paper we propose a tool for high-dimensional approximation based on trigonometric polynomials where we allow only low-dimensional interactions of variables. In a general high-dimensional setting, it is already possible to deal with special sampling sets such as sparse grids or rank-1 lattices. This requires black-box access to the function, i.e., the ability to evaluate it at any point. Here, we focus on scattered data points and grouped frequency index sets along the dimensions. From there we propose a fast matrix-vector multiplication, the grouped Fourier transform, for high-dimensional grouped index sets. Those transformations can be used in the application of the previously introduced method of approximating functions with low superposition dimension based on the analysis of variance (ANOVA) decomposition where there is a one-to-one correspondence from the ANOVA terms to our proposed groups. The method is able to dynamically detected important sets of ANOVA terms in the approximation. In this paper, we consider the involved least-squares problem and add different forms of regularization: Classical Tikhonov-regularization, namely, regularized least squares and the technique of group lasso, which promotes sparsity in the groups. As for the latter, there are no explicit solution formulas which is why we applied the fast iterative shrinking-thresholding algorithm to obtain the minimizer. Moreover, we discuss the possibility of incorporating smoothness information into the least-squares problem. Numerical experiments in under-, overdetermined, and noisy settings indicate the applicability of our algorithms. While we consider periodic functions, the idea can be directly generalized to non-periodic functions as well.

Preconditioned Legendre Spectral Galerkin Methods for the Non-separable Elliptic Equation

The Legendre spectral Galerkin method of self-adjoint second order elliptic equations usually results in a linear system with a dense and ill-conditioned coefficient matrix. In this paper, the linear system is solved by a preconditioned conjugate gradient (PCG) method where the preconditioner M is constructed by approximating the variable coefficients with a (T+1)-term Legendre series in each direction to a desired accuracy. A feature of the proposed PCG method is that the iteration step increases slightly with the size of the resulting matrix when reaching a certain approximation accuracy. The efficiency of the method lies in that the system with the preconditioner M is approximately solved by a one-step method based on the incomplete LU factorization technique with no fill-in, denoted by ILU(0). The ILU(0) factorization of \(M\in {\mathbb {R}}^{(N-1)^d\times (N-1)^d}\) can be computed using \({\mathcal {O}}(T^{2d} N^d)\) operations, and the number of nonzeros in the factorization factors is of \({\mathcal {O}}(T^{d} N^d)\), \(d=1,2,3\). A conclusion of the algorithm is to fast solve the resulting system from the Legendre Galerkin spectral method for Poisson equations with Dirichlet boundary conditions, which has a complexity of \({\mathcal {O}}(N^d)\). To further speed up the PCG method, an algorithm is developed for fast matrix-vector multiplications by the resulting matrix of Legendre-Galerkin spectral discretization, without the need to explicitly form it. The complexity of the fast matrix-vector multiplications is of \({\mathcal {O}}(N^d (\log _2 N)^{2})\). In view that T is independent of N in one dimension and is set to be of order \({\mathcal {O}}(\log _2 N)\) in two and three dimensions, the PCG method has a \({\mathcal {O}}(N^d (\log _2 N)^{2d})\) quasi-optimal complexity for a d dimensional domain with \((N-1)^d\) unknows, \(d=1,2,3\). In addition, a fast direct solver for the three-dimensional Poisson equation is developed, which is of \({\mathcal {O}}(N^{3} (\log _2 N)^{2})\) and improves the existing results on the computational complexity. Numerical examples are given to demonstrate the efficiency of proposed preconditioners and the algorithm for fast matrix-vector multiplications.

A Novel Method for Radar Echo Simulation based on Fast-Constructed Database

This paper presents a novel method for fast calculation of radar echo in near-field regions after the equivalent source has been computed by method of moments (MoM). An easy-to-access near-field database (NFDB) is established, which is built on the auxiliary tetrahedral meshes surrounding the near-field regions of interest. The near-fields calculation (NFC) of arbitrary observation points can be expressed explicitly via the NF-DB. An efficient matrix compression scheme named random sampling-based butterfly factorization (RS-BF) is proposed to speed up the construction of NFDB. With this approach, each group of O(N) elements in the database can be calculated through one fast matrix-vector multiplication operation that has a computational complexity below O(Nlog N). The proposed method can avoid time-consuming point-by-point NFC of the traditional methods. Several numerical examples are presented to demonstrate the accuracy and efficiency of this method. In particular, the echo simulation of a missile-target encounter example is presented to illustrate its capability for practical applications.

Learning in high-dimensional feature spaces using ANOVA-based fast matrix-vector multiplication

&lt;p style='text-indent:20px;'&gt;Kernel matrices are crucial in many learning tasks such as support vector machines or kernel ridge regression. The kernel matrix is typically dense and large-scale. Depending on the dimension of the feature space even the computation of all of its entries in reasonable time becomes a challenging task. For such dense matrices the cost of a matrix-vector product scales quadratically with the dimensionality &lt;inline-formula&gt;&lt;tex-math id="M1"&gt;\begin{document}$ N $\end{document}&lt;/tex-math&gt;&lt;/inline-formula&gt;, if no customized methods are applied. We propose the use of an ANOVA kernel, where we construct several kernels based on lower-dimensional feature spaces for which we provide fast algorithms realizing the matrix-vector products. We employ the non-equispaced fast Fourier transform (NFFT), which is of linear complexity for fixed accuracy. Based on a feature grouping approach, we then show how the fast matrix-vector products can be embedded into a learning method choosing kernel ridge regression and the conjugate gradient solver. We illustrate the performance of our approach on several data sets.&lt;/p&gt;

Efficient second-order ADI difference schemes for three-dimensional Riesz space-fractional diffusion equations

In this paper, a three-dimensional time-dependent Riesz space-fractional diffusion equation is considered, and an alternating direction implicit (ADI) difference scheme is proposed, in which a weighted and shifted Grünwald difference scheme (see Tian et al. (2015) [33]) is utilized for the discretizations of space-fractional derivatives, and a fractional-Douglas-Gunn type ADI method is utilized for the discretization of time derivative. The method is proved to be unconditionally stable and convergent with second-order accuracy both in time and space with respect to a weighted discrete energy norm. Efficient implementation of the method is carefully discussed, and then based on fast matrix-vector multiplications, a fast conjugate gradient (FCG) solver for the resulting symmetric positive definite linear algebraic system is developed. Numerical experiments support the theoretical analysis and show strong effectiveness and efficiency of the method for large-scale modeling and simulations. Finally, a linearized ADI scheme based on second-order extrapolation method is developed and tested for the nonlinear Riesz space-fractional diffusion equation.

An efficient contour integral based eigensolver for 3D dispersive photonic crystal

Numerical simulations play a significant role in studying the properties of dispersive metallic photonic crystals. The dispersive photonic crystals are modeled by the Maxwell equations, and the equations are then discretized by the widely-used Yee’s scheme. After applying certain similarity transformations to the discretized system, the original simulation problem becomes a non-Hermitian eigenvalue problem with clustered eigenvalues. An efficient contour integral (CI) based eigensolver is developed to overcome the difficulties of applying existing methods to solve eigenvalues in designated regions. This efficient method combines the contour integral, the fast matrix–vector multiplication, and efficient linear system solver. The numerical results illustrate the efficiency of our algorithm.

Fast and Accurate Gaussian Kernel Ridge Regression Using Matrix Decompositions for Preconditioning

This paper presents a preconditioner-based method for solving a kernel ridge regression problem. In contrast to other methods, which utilize either fast matrix-vector multiplication or a preconditioner, the suggested approach uses randomized matrix decompositions for building a preconditioner with a special structure that can also utilize fast matrix-vector multiplications. This hybrid approach is efficient in reducing the condition number, exact, and computationally efficient, enabling the processing of large datasets with computational complexity linear to the number of data points. Also, a theoretical upper bound for the condition number is provided. For Gaussian kernels, we show that given a desired condition number, the rank of the needed preconditioner can be determined directly from the dataset.

Toeplitz Monte Carlo

Motivated mainly by applications to partial differential equations with random coefficients, we introduce a new class of Monte Carlo estimators, called Toeplitz Monte Carlo (TMC) estimator, for approximating the integral of a multivariate function with respect to the direct product of an identical univariate probability measure. The TMC estimator generates a sequence \(x_1,x_2,\ldots \) of i.i.d. samples for one random variable and then uses \((x_{n+s-1},x_{n+s-2}\ldots ,x_n)\) with \(n=1,2,\ldots \) as quadrature points, where s denotes the dimension. Although consecutive points have some dependency, the concatenation of all quadrature nodes is represented by a Toeplitz matrix, which allows for a fast matrix–vector multiplication. In this paper, we study the variance of the TMC estimator and its dependence on the dimension s. Numerical experiments confirm the considerable efficiency improvement over the standard Monte Carlo estimator for applications to partial differential equations with random coefficients, particularly when the dimension s is large.

Global-Gate Controlled One-Transistor One-Digital-Memristor Structure for Low-Bit Neural Network

Memristor based neuromorphic computing system has recently attracted enormous attention due to its fast and energy-efficient matrix vector multiplication, thus providing a novel approach to implement neural networks for artificial intelligence. However, the widely studied analogue memristors exhibit major flaws in terms of high conductance variation and nonlinear/asymmetric characteristics. In this work, we develop global gate controlled one transistor one digital memristor (1T1DM) architecture as the basic binary electronic synapse. Inspired by the current research highlights about low-bit networks, we further implement low-bit neuromorphic computing onto our 1T1DM systems by simulation. Compared to the classical analogue type of memristor with one transistor one analogue memristor (1T1R) structure, our 1T1DM network is light-weighted, highly robust and can work well on challenging visual tasks. Besides, benefiting from the global gated device structure, the on-state conductance of the digital memristors in the network can be simultaneously modulated by the controlling gate, offering possibility to tune the power consumption and operation speed while will not increase the circuit complexity.

A preconditioned fast collocation method for a linear bond-based peridynamic model

We develop a fast collocation method for a static bond-based peridynamic model. Based on the analysis of the structure of the stiffness matrix, a fast matrix-vector multiplication technique was found, which can be used in the Krylov subspace iteration method. In this paper, we also present an effective preconditioner to accelerate the convergence of the Krylov subspace iteration method. Using the block-Toeplitz–Toeplitz-block (BTTB)-type structure of the stiffness matrix, we give a block-circulant-circulant-block (BCCB)-type preconditioner. The numerical experiments show the utility of the preconditioned fast collocation method.

A Novel ReRAM-Based Architecture of Field Sequential Color Driver for High-Resolution LCoS Displays

Liquid crystal on-silicon (LCoS) display is one of the most representative micro-display technologies, and is widely adopted in virtual reality (VR) and augmented reality (AR) devices thanks to a relatively simple structure using a semiconductor manufacturing process to realize high-resolution displays. However, the structural complexity to handle color frames by field sequential color (FSC) scheme hinders more widespread adoptions of the LCoS displays. In this article, to resolve the problem, we propose a novel FSC driver architecture using resistive random access memory (ReRAM) that diminishes the driver’s structural complexity with matrix-vector multiplications. The proposed architecture leverages fast matrix-vector multiplications with a memristor crossbar array to expedite the FSC operation that extracts the individual red, green, and blue color sub-frames from an entire image. We present the hardware performance of our architecture that is implemented using the crossbar array and peripheral circuits. Compared to the conventional static random access memory (SRAM)-based architecture, we confirm that the proposed design is much superior in terms of chip size, leakage power, and frame rate in various image resolutions. Specifically, the chip size and leakage power are reduced by up to 96% and 99%, respectively, and the frame rate is improved by up to 36%. We also analyze image quality loss caused by ReRAM read and write noise.

Structured FISTA for image restoration

SummaryIn this paper, we propose an efficient numerical scheme for solving some large‐scale ill‐posed linear inverse problems arising from image restoration. In order to accelerate the computation, two different hidden structures are exploited. First, the coefficient matrix is approximated as the sum of a small number of Kronecker products. This procedure not only introduces one more level of parallelism into the computation but also enables the usage of computationally intensive matrix–matrix multiplications in the subsequent optimization procedure. We then derive the corresponding Tikhonov regularized minimization model and extend the fast iterative shrinkage‐thresholding algorithm (FISTA) to solve the resulting optimization problem. Because the matrices appearing in the Kronecker product approximation are all structured matrices (Toeplitz, Hankel, etc.), we can further exploit their fast matrix–vector multiplication algorithms at each iteration. The proposed algorithm is thus called structured FISTA (sFISTA). In particular, we show that the approximation error introduced by sFISTA is well under control and sFISTA can reach the same image restoration accuracy level as FISTA. Finally, both the theoretical complexity analysis and some numerical results are provided to demonstrate the efficiency of sFISTA.

A stable explicitly solvable numerical method for the Riesz fractional advection–dispersion equations

In this paper, we present a finite difference scheme for solving the Riesz fractional advection-dispersion equations (RFADEs). The scheme is obtained by using asymmetric discretization technique and modify the shifted Grünwald approximation to fractional derivative. By calculating the unknowns in differential nodal-point sequences at the odd and even time-levels, the discrete solution of the scheme can be obtained explicitly. The computational cost for the scheme at each time step can be O(Klog K) by using the fast matrix-vector multiplication with the help of Toeplitz structure, where K is the number of unknowns. We prove that the scheme is solvable and unconditionally stable. We derive the error estimates in discrete l2-norm, which is optimal in some cases. Numerical examples are presented to verify our theoretical results.