linear algebra - How to find selected elements of inverse of a banded matrix without inverting it?

Is it possible to find some selected elements of the inverse of a large sparse matrix without inverting it?

For example consider this Hermitian matrix (as a general case).

n1 = 20; n2 = 50;
n3 = n1 n2;

a1 = RandomComplex[{0. + 0. I, 1. + 1. I}, {n1, n1}];
a2 = ConjugateTranspose[a1];

a3 = RandomComplex[{0. + 0. I, 1. + 1. I}, {n1, n1}];
a3 = (a3 + ConjugateTranspose[a3])/2.;


a0 = SparseArray[{Band[{1, 1}, {n3, n3}] -> {a3}, 
      Band[{1, n1 + 1}, {n3 - n1, n3}] -> {a1}, 
      Band[{n1 + 1, 1}, {n3, n3 - n1}] -> {a2}}, {n3, n3}]

b0 = Inverse[a0]; // AbsoluteTiming

{0.240657, Null}

and I am interested only in a small portion of b0 say

b0[[(n2 - 1) n1 + 1 ;; n2 n1, 1 ;; n1]] // MatrixForm

In case of banded matrix, 'LinearSolve` works much faster.

f = LinearSolve[a0, SparseArray[{i_, i_} -> 1, {n3, n3}]]; // AbsoluteTiming
f[[(n2 - 1) n1 + 1 ;; n2 n1, 1 ;; n1]] // MatrixForm

{0.068483, Null}

The result is same (I didn't put it here). As you can see the execution is significantly faster. But f is using same memory as b0 where the memory required to store a single block (a1,a2) is very little.

Map[ByteCount, {a1, a2, a0, b0, f, b0[[(n2 - 1) n1 + 1 ;; n2 n1, 1 ;; n1]]}]
%/Min[%] // N

{6776, 6776, 6165680, 16000152, 16000152, 6776}

{1., 1., 909.929, 2361.3, 2361.3, 1.}

Is there any way to reduce both time and memory in this case?

Answer

One possibility is to use a rectangular matrix with one block being the identity matrix. Let us consider the example in block matrix form.

$\begin{pmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \\ \end{pmatrix} \cdot \begin{pmatrix} b_{11} & b_{12} & b_{13} \\ b_{21} & b_{22} & b_{23} \\ b_{31} & b_{32} & b_{33} \\ \end{pmatrix} = \begin{pmatrix} I & 0 & 0 \\ 0 & I & 0 \\ 0 & 0 & I \\ \end{pmatrix} $

Each element here is a $n \times n$ block matrix as per our definition. Now let say we want to find $b_{12}$. For that only three scalar equations will be sufficient, which in matrix form is

$\begin{pmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \\ \end{pmatrix} \cdot \begin{pmatrix} b_{12} \\ b_{22} \\ b_{32} \\ \end{pmatrix} = \begin{pmatrix} 0 \\ I \\ 0 \\ \end{pmatrix} $

So the second argument of LinearSolve will be a $3n \times n$ rectangular matrix.

Lets use this with the main problem.

n1 = 20; n2 = 50;
n3 = n1 n2;

a1 = RandomComplex[{0. + 0. I, 1. + 1. I}, {n1, n1}];
a2 = ConjugateTranspose[a1];


a3 = RandomComplex[{0. + 0. I, 1. + 1. I}, {n1, n1}];
a3 = (a3 + ConjugateTranspose[a3])/2.;

a0 = SparseArray[{Band[{1, 1}, {n3, n3}] -> {a3}, 
      Band[{1, n1 + 1}, {n3 - n1, n3}] -> {a1}, 
      Band[{n1 + 1, 1}, {n3, n3 - n1}] -> {a2}}, {n3, n3}]

b0 = Inverse[a0]; // AbsoluteTiming

{0.223080, Null}

f = LinearSolve[a0, SparseArray[{i_, i_} -> 1, {n3, n3}]]; // AbsoluteTiming
f[[(n2 - 1) n1 + 1 ;; n2 n1, 1 ;; n1]] // MatrixForm

{0.068368, Null}

Now using the rectangular argument

icol = SparseArray[Band[{1, 1}] -> IdentityMatrix[n1], {n3, n1}]
f1 = LinearSolve[a0, icol]; // AbsoluteTiming
f1[[(n2 - 1) n1 + 1 ;; n2 n1, 1 ;; n1]] // MatrixForm

{0.009034,Null}

icol is here $(I,0,0 ...)^T$ because the element we want to calculate is in first column ($b_{31}$).

Comparing the ByteCount

Map[ByteCount, {a1, f, f1, f1[[(n2 - 1) n1 + 1 ;; n2 n1, 1 ;; n1]]}]

%/Min[%] // N

{6776, 16000152, 320152, 6776}

{1., 2361.3, 47.2479, 1.}

So this is a bit improvement. The method should also work on dense matrices. Although there must be a even better way to do this.

Time Comparison

I tried with different system size and check the time (in seconds) for full inversion and partial inversion (evaluating one column of the inverse).

$n1=100$

$\begin{array}{rrr} n2 & full & col. \\ \hline 20 & 0.55 & 0.15 \\ 40 & 2.26 & 0.25 \\ 80 & 8.72 & 0.55 \end{array}$

As you can see $time_{col.} > time_{full}/n2$. If you need to calculate multiple blocks, using LinearSolve for full matrix would be better.