How to quickly update the inverse of a lower triangular matrix for a 4 elements changes?

Question

I have the inverse of a lower triangular sparse matrix which is 5000 by 5000. Now only 4 elements at locations (i,k), (i,m), (j,k) and (j,m) have updated values, while the values of all other elements are kept unchanged.

I am looking for a quick approach to obtain the new inverse of this updated lower triangular matrix.

Thanks in advance.

Benson

Answer 1

Your question is vague about what you want to do with $L$. Do you need the inverse of $L$ in each iteration for some specific reason, or just want to solve a bunch of $Lx=b$ systems of equations in each iteration? How many systems per iteration?

In general, the time to multiply $Lx$ for a sparse $L$ matrix should be roughly proportional to the number of nonzero entries in $L$. Similarly, the time to forward solve $x=L\backslash b$ should be roughly proportional to the number of nonzeros in L and roughly equal to the time it takes to perform an $Lx$ multiplication.

Depending on the sparsity pattern, the number of nonzeros in inv(L) may be significantly greater than the number of nonzeros in L. Thus solving $x=L\backslash b$ should typically be faster than $x=L^{-1}b$, even if you've precomputed $L^{-1}$ and amortize this cost over many right-hand sides.

I wrote a simple MATLAB script to generate a random sparse lower triangle 5000 by 5000 matrix and then time the three approaches. A typical result was:

Inverse time=7.451e-0 seconds
L mult time=8.139e-05 seconds
Linv mult time=8.905e-03 seconds
solve time=1.075e-04 seconds

That is, $b=Lx$ and $x=L\backslash b$ both took about 1e-4 seconds ($Lx$ was about 20% faster) while $x=L^{-1}b$ took almost 10 times longer (9e-3 seconds, not counting the time to compute the inverse itself.)

Since $L\backslash b$ is very fast, if you're solving a few thousand $Lx=b$ systems of equations per iteration, then you would better off in MATLAB on this system updating the elements in $L$ and doing forward solves after each update.

Notice that you could use 5000 forward solves to find the inverse of $L$, and this would take around 0.5 seconds, faster then the time it took to compute the inverse of L in MATLAB with the inv command. I don't know why inv was so slow on this example- it's surprising.

From what you've described it seems likely that the software you're using is not implementing these sparse matrix operations in the most efficient manner.

Answer 2

@Brian. I tried several approaches but the CPU time does not reduce effectively. My problem is given as follows:

InvA = A;

for I = 1 : 7000

For a given new deltaA, A1 = A + deltaA.

while {converge?}

   solve A1*x_k+1 = b/x_k.

end

end

In the inner loop, I can use x_k+1 = A1(b/x_k), or x_k+1=InvA1*(b/x_k). Since there are only 4 non-zero elements in deltaA, I divided A1 into 4 sub-matrices and try to make better use of InvA to obtain InvA1. I found to get InvA1 is faster for my case. But it still takes about 18 minutes. I am thinking if there is a faster way to obtain InvA1?

Thanks a lot again. You have a good day! Benson