Linearly scaling direct method for accurately inverting sparse banded matrices

Abstract

Abstract. In many problems in Computational Physics and Chemistry, one finds a special kind of sparse matrices, called banded matrices. These matrices, which are defined as having non-zero entries only within a given distance from the main diagonal, need often to be inverted in order to solve the associated linear system of equations. In this work, we introduce a new O(n) algorithm for solving such a system, with the size of the matrix being n × n. We derive analytical recursive expressions that allow us to directly obtain the solution. In addition, we describe the extension to deal with matrices that are banded plus a small number of non-zero entries outside the band, and we use the same ideas to produce a method for obtaining the full inverse matrix. Finally, we show that our new algorithm is competitive, both in accuracy and in numerical efficiency, when compared with a standard method based on Gaussian elimination. We do this using sets of large random banded matrices, as well as the ones that appear in the calculation of Lagrange multipliers in proteins.

Exact and efficient calculation of lagrange multipliers in biological polymers with constrained bond lengths and bond angles: Proteins and nucleic acids as example cases

“…To show this, first, we will prove that the value of the vectors p and q can be obtained in O(N c ) operations. Then, we will show that the same is true for all the non-zero entries of matrix R, and finally we will briefly discuss the results in [18], where we introduced an algorithm to solve the system in (2.6) also in O(N c ) operations.…”

“…When the semi-band width m is not constant along the whole matrix, things are more complicated and the cost is always between O(N c m 2 min ) and O(N c m 2 max ), depending on how the different rows are arranged. In general, we want to minimize the number of zero fillings in the process of Gaussian elimination (see [18] for further details), which is achieved by not having zeros below non-zero entries. This is easier to understand with an example: Consider the following matrices, where Ω and ω represent different non-zero values for every entry (i.e., not all ω, nor all Ω must take the same value, and different symbols have been chosen only to highlight the main diagonal):…”

“…However, as we showed in the previous sections, the fact that many biological molecules, and proteins in particular, are essentially linear, allows to index the constraints in such a way that the matrix R in eq. (2.6) is banded and use different techniques for solving the problem which require only O(N c ) floating point operations [18].…”

“…In this work, we show that the fact that many interesting biological molecules are esentially linear polymers allows to calculate the Lagrange multipliers in order N c operations (for a molecule where N c constraints are imposed) in an exact (up to machine precision), non-iterative way. Moreover, we provide a method to do so which is based in a clever ordering of the constraints indices, and in a recently introduced algorithm for solving linear banded systems [18]. It is worth mentioning that, in the specialized literature, this possibility has not been considered as far as we are aware; with some works commenting that solving this kind of linear problems (or related ones) is costly (but not giving further details) [19][20][21], and some other works explicitly stating that such a computation must take O(N 3 c ) [22] or O(N 2 c ) [23,24] operations.…”

A survey of the parallel performance and accuracy of Poisson solvers for electronic structure calculations

“…whereL is a sparse matrix (taking advantage of the sparsity of a system of equations can greatly reduce the numerical complexity of its solution 57 ), y y y is known (y = −ρ/ε 0 , in this case) and x x x is the quantity we are solving for, in this case the electrostatic potential. Equation (Eq.…”