Updating inverse of a lower triangular matrix

Dear All,

I need to solve a matrix equation Ax=b, where the matrix A is a lower triangular matrix and its dimension is very big (could be 10000 by 10000). Now I need to change a row of A and solve Ax=b again (this change will be many times). In order to speed up the calculation, a good approach is to calculate the inverse of matrix A and use the substitution to solve x. But when a row of matrix A is changed, I need to update the inverse of A. Is there a good idea to update the inverse of matrix A?

Answer 1

[John D'Errico]

NO. Computing the inverse of A, if it is lower triangular is a BAD idea.

Solving the problem x = A\b is a forward substitution, so fast as hell. On the order of a matrix*vector multiply in terms of the computational load. I.e., essentially an O(n^2) operation.

Wasting the time to do do an update of an inverse matrix, so that you can then do a matrix multiply is just silly.

A = tril(rand(10000));
b = rand(10000,1);

timeit(@() A*b)
ans =
     0.067386

timeit(@() A\b)
ans =
     0.090263

So the backslash took only slightly more time than a multiply. The update of a matrix inverse would take way more time than that. And for no good reason at all, only to still need to do a matrix*vector multiply.

For example, one can compute the change to a matrix inverse under a rank 1 update. (which is exactly what you want to do.) The classic solution is due to Sherman & Morrison. (I remember Woodbury in the name when I learned about it a zillion years ago, but Sherman-Morrison is a special case.)

https://en.wikipedia.org/wiki/Sherman–Morrison_formula

Now, look carefully at the computations in the formula shown. While you could do them, you need to recognize that they would involve more work than a simple matrix*vector multiply would entail. And since the forward substitution itself is trivial, you cannot possibly gain here.

Answer 2

[Bei Gou]

Dear John,

Thank you very much for your quick responses. Your idea sounds great. Do you think the calculation of forward substitution could be even faster if we use C code to conduct it instead of Matlab?

[John D'Errico]

Please use the comment capability to ask questions, rather than adding answers.

But surely this depends on your capability in writing C code. :) Note that these tools are already compiled code. Can you do better than someone who seriously does know what they are doing?

There is ONE way you can gain though. Use linsolve, and tell it that your matrix is lower triangular.

A = tril(rand(10000));
b = rand(10000,1);

opts.LT = true;
timeit(@() linsolve(A,b,opts))
ans =
     0.057655

timeit(@() A\b)
ans =
      0.12932

The point is if we tell linsolve in advance not to bother checking if the matrix has some special shape like lower triangular, then we can gain some speed. We can't tell backslash that. So backslash takes a fair amount of time just deciding how to solve the problem. If you know in advance something about the matrix, linsolve exists to let you use that information.

[Bei Gou]

Thanks a lot, John, for your helpful comments.

[Bei Gou]

Hi, John,

I am trying to implement your idea for my problem, but I find there is another issue we have to consider. The linear equation is Not actually a linear one: it is like Ax=b(x), where A is a constant sparse invertible matrix, and b is a column vector which is a function of x. So to solve x, I need to do the following iterations:

1. Initialize x.

2. Solve x_k+1=A/b(x_k).

3. Let x_k+1=x_k.

4. Repeat above until |x_k+1-x_k| is small enough.

So in the above process, we need to use several times of Forward Substitution. If we use the updated inverse of A, we can save time on each iteration.

How do you think?

Thanks a lot again.

Bei