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

40 visualizaciones (últimos 30 días)
I have a large sparse matrix A and have gotten its inverse matrix inv(A) . Then I need to change an element value to get a new matrix, A1. I am trying to get the inverse of A1. Is there any way to do it , rather than recalculate the inv(A1)? Could I get some benefits from A or inv(A)?
e.g. (not a sparse matrix)
A = magic(5);
invA = inv(A);
A1 = A;
A1(1,5) = 1;
invA1 = inv(A1);% is there any way to do this?
I know this may be a math problem. It seems impossible for me. I just would like to know other's viewpoints of it. Thank you.

Respuesta aceptada

John D'Errico
John D'Errico el 1 de Ag. de 2020
Editada: John D'Errico el 1 de Ag. de 2020
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.
  2 comentarios
Bruno Luong
Bruno Luong el 1 de Ag. de 2020
Editada: Bruno Luong el 1 de Ag. de 2020
Explicit compute INV is raeely the good idea. We see that the 1-rank update is reduces to something like this:
iAupdate = iA + a * u * v';
where a is scalar and u and v' are (n x 1) vector
Down of the line, at the end of the day this inverse has to be multiplied with some vector x, so
iAupdate*x = iA*x + a * u * v' * x.
Well the most economic to perform that is in this order:
iAupdate*x = iA*x + (a * (v' * x)) * u.
So actually the update of matrix-vector product need O(n) operations where as update explicitly the inverse-matrix needs O(n^2).
There we go, it explains why it is also not terribly useful for CPU as John has pointed out.
This kind of formula (applied to matrix-vector product) gives rise for BFGS formula, Kamal filtering, Householder transformation, etc... which better to be mastered and applied for each specific problem.
wei zhang
wei zhang el 2 de Ag. de 2020
Thank you for your clarification, especially the time consuming comparision. I would like to give up the idea to update. Direct inverse recalculation seems to be the best solution.

Iniciar sesión para comentar.

Más respuestas (2)

Vladimir Sovkov
Vladimir Sovkov el 1 de Ag. de 2020
I am not sure how much profitable it is numerically, but the Sherman Morrison theorem can be a way, see https://en.wikipedia.org/wiki/Sherman%E2%80%93Morrison_formula
In your example, u vector there should be chosen with only the 1st nonzero element and v vector with only the 5th nonzero element so that u1*v5=(1-A15).
  2 comentarios
Vladimir Sovkov
Vladimir Sovkov el 1 de Ag. de 2020
In your example, it turns to (try to rederive it by yourself to ensure that I have not mistaken somewhere)
wei zhang
wei zhang el 2 de Ag. de 2020
Thank you for your introduction for sherman Morrison theorem and kindly inference about the formula.

Iniciar sesión para comentar.


Bruno Luong
Bruno Luong el 1 de Ag. de 2020
Editada: Bruno Luong el 1 de Ag. de 2020
Apply this Sherman Morisson's formula in the page provided by Vladimir, if one set a single change: A(i,j) to new value newAij, the update should goes like this
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(norm(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
Test script
% testupdate.m
% Random matrix test
A = rand(5);
iA = inv(A);
% Random change
i = randi(size(A,1));
j = randi(size(A,2));
newAij = 10;
% Update matrix and its inverse
[Amod, iAupdate] = updateinv(A, i, j, newAij, iA);
% Check
iAmod = inv(Amod)
iAupdate
relerr = norm(iAupdate-iAmod,'fro')/norm(iAmod,'fro')
Result:
>> testupdate
iAmod =
0.3427 -0.2397 0.0926 -0.1965 -0.0649
-3.1551 3.3784 0.1696 1.3277 -0.5608
-4.7509 1.9002 0.1144 2.6354 0.9978
5.2877 -3.8865 -0.3091 -2.5948 1.2505
-1.4484 1.5561 0.1196 1.7662 -1.5276
iAupdate =
0.3427 -0.2397 0.0926 -0.1965 -0.0649
-3.1551 3.3784 0.1696 1.3277 -0.5608
-4.7509 1.9002 0.1144 2.6354 0.9978
5.2877 -3.8865 -0.3091 -2.5948 1.2505
-1.4484 1.5561 0.1196 1.7662 -1.5276
relerr =
4.0984e-16
From my experience this formula is not terribly useful for 2 reasons
  • I rarely store inverse of matrix, epsecially medium/large size.
  • One can lost precision when apply this formula multiple times.
  1 comentario
wei zhang
wei zhang el 2 de Ag. de 2020
Thank you for your answer very much! It is pity I couldn't set your and John D'Errico's answer both at the answer. Your answer is direct and great. His answer shows some time profiling to clarify this problem more clearly, including your comments there. Forgive me for the accepting.

Iniciar sesión para comentar.

Categorías

Más información sobre Logical en Help Center y File Exchange.

Etiquetas

Productos


Versión

R2019a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by