2
$\begingroup$

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?

$\endgroup$
3
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.