Solving sparse linear systems with approximate inverse preconditioners on analog devices

The iterative solution of sparse linear systems [1] is a fundamental task across many applications in science, engineering, and optimization. For decades, progress in speeding up preconditioned iterative solvers was primarily driven by constructing more effective preconditioning algorithms and via steady microprocessor performance gains. The last decade saw the advent of fast power-efficient low-precision hardware, such as Graphics Processing Units (GPUs), which gave rise to new opportunities to accelerate sparse iterative solvers [2, 3, 4, 5, 6, 7, 8]. While these advances have led to remarkable gains in the performance of sparse iterative solvers, conventional CMOS-based digital microprocessors possess inherent scaling and power dissipation limitations. It is therefore imperative that we explore alternative hardware architectures to address current and future challenges of sparse iterative solvers. One such alternative paradigm is based on analog devices configured as crossbar arrays of non-volatile memory. Such arrays can achieve high degrees of parallelism with low energy consumption by mapping matrices onto arrays of memristive elements capable of storing information and executing simple operations such as a multiply-and-add. In particular, these devices can perform Matrix-Vector Multiplication (MVM) in near constant time independent of the number of nonzero entries in the operand matrices [9, 10]. Indeed, considerable computational speed-ups have been achieved with analog crossbar arrays [11, 12, 13, 14, 15, 16, 17, 18, 19], mostly for machine learning applications that can tolerate the noise and lower precision of these devices.

Analog designs with enhanced precision have been proposed [20, 21] to meet the accuracy demands of scientific applications, but the improvement in precision comes at the cost of increased hardware complexity and energy consumption. Memristive hardware with precision enhancement has been exploited in the context of Jacobi iterations to solve linear systems stemming from partial differential equations (PDEs) [22], and as an autonomous linear system solver [23] for small dense systems. The latter approach was used to apply overlapping block-Jacobi preconditioners in the Generalized Minimal RESidual (GMRES) iterative solver [20] with bit-slicing for increasing precision. Inner-outer iterations [24, 25, 26] have been considered for dense systems, where analog arrays were used for the inner solver while iterative refinement was used as an outer solver in the digital space.

The main contribution of this paper is a flexible preconditioning framework for inexpensively solving a large class of sparse linear systems by effectively utilizing simple analog crossbar arrays without precision-enhancing extensions. Our approach combines flexible iterative solvers and approximate inverse preconditioners on a hybrid architecture consisting of a conventional digital processor and an analog crossbar array accelerator. We analytically and experimentally show that this combination can solve certain sparse linear systems at a fraction of the cost of solving them on digital processors. While the relative advantage of analog hardware grows with the size of the linear systems, we demonstrate the potential for an order of magnitude speedup even for systems with just a few hundred equations.

While our method can be used to provide a fast stand-alone sparse solver, it also has exciting applications in the context of the upcoming exascale platforms [27]. Extreme-scale solvers are likely to be highly composable [28] and hierarchical [29] with multiple levels of nested algorithmic components. A solver like the one proposed in this paper, running on analog accelerator-equipped nodes of a massively parallel platform, is an ideal fast and low-power candidate for a local subdomain or coarse-grid solver to tackle the local components of a much larger sparse linear system.

We begin by briefly describing our architecture model and its key properties. Figure 1 illustrates a hybrid digital-analog architecture that combines a conventional microprocessor system with an analog crossbar array for performing MVM. The crossbar consists of rows and columns of conductors with a memristive element at each row-column intersection [30]. The conductance of these elements can be set, reset, or updated in a non-volatile manner. Matrix operations are performed by mapping the $n\times n$ operand matrix $M$ onto the crossbar array, where the conductance $G_{ij}$ at the intersection of row $i$ and column $j$ represents the value $M[i,j]$ (i.e., the $(i,j)$-th entry of the matrix $M$) after a suitable scaling.

The MVM operation $y=Mr$ can be emulated by sending pulses of $V_{in}$ volts along the columns of the crossbar such that the length $t_{j}$ of the pulse along column $j$ is proportional to the $j$th entry $r[j]$ of vector $r$, with suitable normalization. Following Ohm’s Law, this contributes a current equal to $V_{in}G_{ij}$ on the conductor corresponding to row $i$ for a duration $t_{j}$. Following Kirchoff’s Current Law, the currents along each row accumulate and can be integrated over the time period equivalent to the maximum pulse length using capacitors, yielding a charge proportional to $\sum_{j=1}^{n}M[i,j]r[j]$, which in turn is proportional to $y[i]$. This integrated value can be recovered as a digital quantity via an analog-to-digital converter (ADC), and yields an approximation $\widehat{y}[i]$ to $y[i]$.

The above procedure involves multiple sources of nondeterministic noise so that the output $\widehat{y}\in\mathbb{R}^{n}$ is equivalent to a low precision approximation of $y$. Writing $M$ to the crossbar array incurs multiplicative and additive write noises $N^{Wm}\in\mathbb{R}^{n\times n}$ and $N^{Wa}\in\mathbb{R}^{n\times n}$, respectively, and the actual conductance values at the crosspoints of the array reflect $\widehat{M}=M\odot(I+N^{Wm})+N^{Wa}$, where $\odot$ denotes element-wise multiplication. Similarly, digital-to-analog conversion (DAC) of the vector $r$ into voltage pulses suffers from multiplicative and additive input noises $N^{Im}\in\mathbb{R}^{n}$ and $N^{Ia}\in\mathbb{R}^{n}$, respectively. As a result, the matrix $\widehat{M}$ is effectively multiplied by a perturbed version of $r$ given by $\widehat{r}=r\odot(\mathbf{1}+N^{Im})+N^{Ia}$, where $\mathbf{1}$ is a vector of all ones.

A characteristic equation to describe the output $\widehat{y}$ of an analog MVM $Mr$ can be written as

where $N^{Om}\in\mathbb{R}^{n}$ and $N^{Oa}\in\mathbb{R}^{n}$ denote the multiplicative and additive components of the output noise, respectively. These components reflect the inherent inexactness in the multiplication based on circuit laws and current integration, as well as the loss of precision in the ADC conversion of the result vector. Equation (1) can be simplified as

where $\widehat{N}^{Oa}$ is the overall effective additive noise that captures all the lower order noise terms from Equation (1) and the nondetermistic error matrix $E$ is given by

here $\otimes$ denotes outer product. We provide additional technical details in Appendix A.

The exact characteristics and magnitudes of the noises depend on the physical implementation [12, 25] of the analog array and the materials used therein, the details of which are beyond the scope of this paper but they are generally captured by our model. Regardless of the implementation, the noises are always device dependent and stochastic.

Once the matrix $M$ is written to the crossbar array, analog MVMs involving the matrix can be performed in $O(1)$ time, resulting in an $O(n^{2})$ speedup over the equivalent digital computation if $M$ is dense. Therefore, a linear system solver that can overcome the analog MVM error, even if it requires many more MVM operations than its digital counterpart, can substantially reduce the cost of the solution.

The goal of preconditioning is to transform an $n\times n$ system of linear equations $Ax=b$ to enable an iterative procedure to compute a sufficient approximation of $x$ in fewer steps. Algebraically, preconditioning results in solving the system $MAx=Mb$, where $M\in\mathbb{R}^{n\times n}$ is an approximation of $A^{-1}$. Since $A^{-1}$ is dense in general, most popular general-purpose preconditioners are applied implicitly; e.g., the ILU preconditioner generates a lower (upper) triangular matrix $L$ ($U$) such that $A\approx LU$. The preconditioner $M=U^{-1}L^{-1}$ is then applied by forward/backward substitution [1].

Our experiments were conducted in a Matlab environment (version R2020b) with 64-bit arithmetic on a single core of a 2.3 GHz 8-Core Intel i9 machine with 64 GB of system memory. We used a Matlab version of the publicly available simulator [41] with a PyTorch interface for emulating the noise, timing, and energy characteristics of an analog crossbar array. The simulator models all sources of analog noise outlined in Section II as scaled Gaussian processes. Using Matlab notation, the components of the matrix write noise were modeled as $\mathtt{randn(\cdot)\times 5.0\mathtt{e}-3}$, and those of the input and output noises were both modeled as $\mathtt{randn(\cdot)\times 1.0\mathtt{e}-2}$; these are the default settings in the simulator based on currently realizable analog hardware [16]. The number of bits used in the ADC and DAC was set to $7$ and $9$, respectively. Table I shows the default parameters used throughout our experiments to simulate the analog device and to construct the approximate inverse preconditioner.

Our test problems originate from discretizations of model PDEs. The first two matrices are derived from a Finite Element (FE) discretization of Poisson’s equation on $\Omega\subset\mathbb{R}^{2}$ with homogeneous Dirichlet boundary conditions

where $\partial$ denotes the partial derivative with respect to an independent variable and $f(x,y)$ denotes the load (or source) vector. We consider two different domains: a square domain $\Omega:=[0,1]^{2}$ and a circular domain $\Omega=\{[x,y]\in\mathbb{R}^{2}\;\big{|}\;x^{2}+y^{2}=1\}$ of radius one centered at the origin. The Laplace operator $\left(\dfrac{\partial^{2}}{\partial x^{2}}+\dfrac{\partial^{2}}{\partial y^{2}}\right)u(x,y)$ is discretized by linear finite elements while the load vector is set equal to the constant source function $f(x,y)\equiv 1$. Figure 2 plots the sparsity patterns of the discretized Laplacian matrices $A$ and the surfaces of the solutions $u(x,y)$ for the square ($n=625$) and circular ($n=362$) meshes. Our third test matrix stems from a Finite Difference (FD) discretization of the Laplace operator in the unit cube with homogeneous Dirichlet boundary conditions and $n=8^{3}$. The average number of nonzero entries per row in the discretized sparse matrices $A$ and in the corresponding preconditioners $M$ are listed in Table II, which shows that the fill-in factor (i.e., the ratio $\mathrm{nnz}(M)/\mathrm{nnz}(A)$) ranges from $14$ to $24$. Our test matrices are relatively small because we are limited by the speed of the Matlab version of the simulator for the analog crossbar arrays; however, our results are meaningful because the relative advantage of analog hardware, in general, only increases with larger matrices.

In our experiments, we mainly compare the following iterative solvers: ($a$) Richardson iterations without preconditioning; ($b$) Richardson iterations with approximate inverse preconditioning applied in the IEEE 754 standard (denoted by “Richardson+d-ai”); and ($c$) Richardson iterations with approximate inverse preconditioning applied through simulated analog hardware (denoted by “Richardson+h-ai”). We observe from the results in Table II that standard Richardson iterations fail to converge within $m_{\mathrm{it}}=50$ iterations for all three test problems, but converge rapidly with approximate inverse preconditioning. This is not surprising since the spectral radius of the matrix $I-A$, i.e., $\rho(I-A)$, is very close to one in all cases. The convergence behavior of the two preconditioned variants is almost identical for the FE matrices, but hybrid preconditioned Richardson iteration requires more than twice the number iterations compared to the digital variant for the FD matrix. Overall, the hybrid preconditioned Richardson requires $5\times$ to $10\times$ fewer digital FLOPs compared to the purely digital version.

Figure 3 plots the relative residual norm after each Richardson iteration both without preconditioning and with the various preconditioners considered in this paper. As expected, hybrid preconditioning requires more iterations; however, this trade-off is highly beneficial because analog hardware allows for a rapid application of the preconditioner. This is evident in Figure 4 which shows that, on an average, hybrid preconditioning requires substantially fewer digital FLOPs to converge. For reference, we have also included Richardson iterations with an ILU(0) preconditioner [1] applied in the IEEE 754 standard.

Figure 5 plots the number of digital FLOPs versus the relative residual norm for $\mathtt{nnz}_{\mathrm{AI}}$ values of $20\times\mathtt{nnz}(A)$, $40\times\mathtt{nnz}(A)$, and $60\times\mathtt{nnz}(A)$. A denser preconditioner speeds up the convergence rate at the cost of additional FLOPs per iteration in the digital implementation. For the square FE mesh, the hybrid implementation yields a greater improvement in FLOPs compared to its digital counterpart as the preconditioner is made denser because the time to apply $M$ on the analog device does not increase with nnz($M$). However, this trend is not observed for the other matrices and diminishes with increasing preconditioner density because the analog noises limit the accuracy of the MVM operation $Mr$ even if $M$ is highly accurate. As a result, the number of iterations does not always decline as nnz($M$) increases.

Figures 3–5 show that Richardson iterations preconditioned with approximate inverses on the hybrid architecture outperform preconditioned Richardson iterations on the digital architecture by greater margins as the condition number (listed in Table II) of $A$ grows larger. Figure 5 also shows that our method can use denser preconditioners more effectively for matrices with higher condition numbers. These trends bode well for the potential use of our proposed preconditioning framework on a hybrid architecture in real-life problems that are likely to be larger and more ill-conditioned.

Figure 6 plots the impact of varying the fill-in factor for the FE square matrix on the time required by the key operations of Algorithm 1, such as constructing $M$, writing it to the analog crossbar, performing an MVM with it in both the digital and hybrid settings, and performing a digital MVM with the coefficient matrix $A$. The sequential time for the construction of the preconditioner in MATLAB is reported for academic purposes only; in practice, this step would utilize compiler optimizations and the ample availability of parallelism [34, 35] to finish much faster. The time for writing $M$ is roughly equivalent to that of performing seven digital MVMs with it, implying that the benefits of the hybrid architecture start accruing after seven iterations. Note that the size of this matrix $M$ is only $625\times 625$. The current state of analog hardware technology is capable of realizing crossbar arrays of sizes up to $4000\times 4000$ [16]. These can accommodate much larger preconditoner matrices for which the analog MVM times would still be similar to those in Figure 6 and the write times will increase as $O(n)$, but the digital MVM times would be much higher, growing as $O(n^{2})$.

Some hardware parameters, such as the number of DAC and ADC bits, can have a profound impact on the time and energy consumption of the analog device, and therefore on its overall performance. We vary the number of DAC and ADC bits and plot the results in Figure 7 to justify our choice of 7 DAC bits and 9 ADC bits. Fewer bits tend to lead to divergence, while more bits result in increased time and energy consumption with little or no improvement in convergence, as the accuracy of the analog computations is limited by other sources of noise.

With the slowing of Moore’s law [42] for digital microprocessors, there has been considerable interest in exploring alternatives, such as analog computing, for speeding up computationally expensive kernels. While simple analog crossbar arrays have been shown to be effective in many applications involving dense matrices, it has been challenging to address sparse matrix problems, such as solving sparse linear systems, because they do not naturally map to dense crossbar arrays and generally have a low tolerance for the stochastic errors pervasive in analog computing. This paper proposes, analyzes, and experimentally evaluates a preconditioning framework that exploits inexpensive MVM on analog arrays to substantially reduce the time and energy required for solving an important practical class of sparse linear systems. These and similar efforts are critical for meeting the continually growing computational demands of various applications of sparse solvers and for preparing these solvers for the fast but energy-constrained embedded and exascale [27] systems of the future.

We would like to thank Tayfun Gokmen, Malte Rasch, and Shashanka Ubaru for helpful discussions and input that substantially contributed to the content of this paper.