Inverting a huge sparse banded matrix

Question

I have a matrix of $9,200 \times 9,200$ elements. I have approximately $90$ of these matrices to invert.

The reason for this is I am running a nonlinear regression on a problem with significant errors in the independent and dependent variables. Each run, I update the covariance matrix of the dependent variable, and therefore requires inversion to include in the weighted regression problem.

The matrix has a sparse banded structure. The structure of the first $1,000$ values looks like this and the structure of the first $20$ values looks like this. That is, sets of $2 \times 2$ blocks running down the main diagonal, and then also off diagonals seperated by $368$ zeros. The matrix is also symmetric positive definite.

I believe this is called a banded matrix. Given only 0.5% of the matrix is non-zero, I believe there is a faster way of computing this, but am struggling.

Answer 1

To expand on George's answer, if we have a matrix consisting of diagonal matrices $D_{ij}$ like this :

$$M = \left[\begin{array}{ccc} D_{11}&\cdots &D_{1n}\\\vdots&\ddots&\vdots\\D_{n1}&\cdots&D_{nn} \end{array}\right]$$

We can always find a permutation similarity

$$A = PMP^{t}$$

So that $A = \left[\begin{array}{ccc} B_{1}&0 &0\\0&\ddots&0\\0&0&B_{m} \end{array}\right]$

In Gnu Octave or Matlab code this permutation transformation seems to be possible to express as a sparse matrix

P = sparse(((1:n)'+(0:(m-1))*n)(:),((0:(n-1))'*m+(1:m))(:),ones(n*m,1))'

Now inverting $A$ will be a simple case of inverting each of the $B$:s separately.

Answer 2

I actually believe I have answered my own question.

The key to shortening this time down it recognizing the repeating structure of the matrix will lead to most row-column multiplications going to $0$.

Note: The matrix is repeating units of $2\times2$ blocks.

We can consider the symmetric matrix to be a set of $2\times9200$ block matrices arranged as columns. These will only be non-zero when multiplied by themselves, or the $i\times368$ column (there is a diagonal of $2\times2$ block matrices every 368 values).

We can arrange these columns into 184 matrices of 25-stacked block matrices. We can then realise that each of these columns have a bunch of bands of zeros too that can be removed.

We can therefore treat this covariance matrix as a set of 184 $50\times50$ matrices, whose inversion is far quicker.

Currently on python I can get 21 seconds to invert the $9200\times9200$ matrix, and 0.75 seconds using the above method.