MATLAB: Updating inverse of a lower triangular matrix

lower triangular matrixupdating inverse

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?

Best Answer

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.

Related Solutions

MATLAB: Finding the pseudo inverse of a matrix

If the matrix A does not have full rank, there is no inverse. In consequence you cannot find any B, which satisfies A*B=eye . This is the definition of the rank, of invertible and there cannot be an "alternative".

This means, that the question is not meaningful. It is like asking for the inverse of 0. There is none.

MATLAB: Is there any method to calculate the inverse of matrix which changed a few values

Bruno is correct in how to solve the problem, and also SOME of the issues with it. If you use this tool multiple times, the floating point errors will begin to accumulate.

There is however, still a serious problem. That is time. I recall looking at this formula long ago, when I first learned about it. (35 - 40 years ago.) At the time, I concluded it might not be such a great tool, both due to the numerical issues, as well as the time required to perform the computation.

While the idea of a rank 1 update for the inverse of a matrix seems interesting, it has serious problems in practice. Gosh, it must be more efficient that just recomputing the inverse of the matrix, right? Always test these things.

I ran a test of this code. In my tests, I computed a random matrix of size N, then the inverse of that matrix. Next, I chose a random element to modify to some new random value.

Finally, I used timeit to measure the time required for a matrix inverse of that matrix. Then I used timeit to measure the time required for this update to the matrix inverse. The results are found in the table below.

     N       invtime      updatetime    update ratio
    ____    __________    __________    ____________
      50    4.6589e-05    9.3269e-05       2.0019   
     100    0.00012982    0.00048521       3.7375   
     200     0.0004016     0.0018684       4.6524   
     400     0.0014756     0.0073921       5.0094   
     600     0.0041071      0.017325       4.2183   
     800     0.0070771      0.032447       4.5848   
    1000      0.012438      0.055685       4.4771   
    1300      0.025226       0.11205       4.4419   
    1600      0.045956       0.31944       6.9509   
    1900      0.075096       0.81713       10.881   
    2200       0.11136        1.6726       15.019   
    2500       0.16798        3.0376       18.083   
    3000       0.29141        6.2889       21.581

So column 1 in that table is the order of the matrix used in the test.

invtime is the time used to compute one matrix inverse of a matrix of that size.

updatetime is the time required for the function updateinv to do its work. As you should see, updateinv is considerably SLOWER than a simple, direct inverse of the matrix.

The final column is the relative cost. As you should see, updateinv actually starts to get worse as the matrix size itself grows.

In the loglog plot, you should see that the matrix inverse itself takes consistently LESS time to perform. As well, you can see the differential gets worse for larger matrices.

If you don't trust the data I provide, run a test yourself. Even a 2Kx2K matrix inverse is not that costly to compute.

A = rand(2000);
timeit(@() inv(A))
ans =
     0.084881

As you can see, this mirrors the numbers in column 2 in my table. I stopped at 3K square matrices, because the time required for the update was starting to get expensive, and I wanted to post this response.

I've attached the script I used to do these time tests, in case you wish to verify my results.

So you seriously need to consider if this is a good idea. While I always strongly advise considering if you even want to compute the matrix inverse at all as there are better things to do almost always, updating that inverse using the code posted by Bruno was never a savings in time. If you will perform multiple updates, then you are further behind, since now you will also incur penalties due to accumulation of the errors.

Finlly, should Bruno wish to soend some time in optimizing his code, I would add that the profile tool tells me most of the time required seems to be in just in creating the tolerance.

A = rand(2000);
timeit(@() norm(A))
ans =
       1.0575
       
timeit(@() normest(A))
ans =
     0.011862
     
norm(A)
ans =
       1000.1
normest(A)
ans =
       1000.1

So, if I just modify one line from the code written by Bruno, replacing norm with a call to normest,

function [Amod, iAupdate] = updateinv(A, i, j, newAij, iA)
% INPUTS:
%   A (n x n) matrix
%   i, j, newAij are scalars: A(i,j) will change to newAij
%   iA inverse of A
% OUTPUTS:
%   Amod: n x n matrix subjected to this modification
%   iAupdate: inverse of Amod, updated using Shermanâ€“Morrison formula
Amod = A;
Amod(i,j) = newAij;
d = newAij-A(i,j);
c = 1+d*iA(j,i);
tol = max(size(A)) * eps(normest(A));
if abs(c) < tol
    warning('update matrix might be inaccurate');
    if c > 0
        c = tol;
    else
        c = -tol;
    end
end
iAupdate = iA - d/c * (iA(:,i)*iA(j,:));
end

the times required have now inverted in their behavior.

     N       invtime      updatetime    update ratio
    ____    __________    __________    ____________
      50    4.0205e-05    9.2221e-06      0.22938   
     100     0.0001385     3.005e-05      0.21697   
     200    0.00039433    7.6963e-05      0.19517   
     400     0.0018709    0.00019188      0.10256   
     600     0.0040269     0.0012491      0.31018   
     800     0.0071675     0.0025266       0.3525   
    1000      0.012501     0.0052236      0.41786   
    1300      0.023901      0.011431      0.47824   
    1600      0.044222       0.01764       0.3989   
    1900      0.073962      0.024691      0.33384   
    2200       0.10989      0.033139      0.30156   
    2500       0.17115      0.043476      0.25402   
    3000       0.29513      0.063203      0.21415

The cost of using normest is it is only an approximate estimate of the norm of A. But here we see we could have done several updates to the inverse in the time required for one inverse. The question of accumulation of errors is another issue, worth some study in itself.

Best Answer

Related Solutions

MATLAB: Finding the pseudo inverse of a matrix

MATLAB: Is there any method to calculate the inverse of matrix which changed a few values

Related Question