mirror of
https://github.com/TREX-CoE/Sherman-Morrison.git
synced 2024-12-25 13:53:56 +01:00
Initial commit.
This commit is contained in:
François Coppens
commit
3704856b10
71
SM-Coppens-1.m
Normal file
71
SM-Coppens-1.m
Normal file
@ -0,0 +1,71 @@
|
||||
## Based on Algorithm 3 from P. Maponi,
|
||||
## p. 283, doi:10.1016/j.laa.2006.07.007
|
||||
|
||||
clc ## Clear the screen
|
||||
## Define the matrix to be inverted. This is example 8 from the paper
|
||||
## In the future this matrix needs to be read from the function call arguments
|
||||
#A=[1,1,-1; ...
|
||||
# 1,1,0; ...
|
||||
# -1,0,-1];
|
||||
#A0=diag(diag(A)); ## The diagonal part of A
|
||||
|
||||
### The modified example that gives all singular updates at some point
|
||||
#A=[1,1,1; ...
|
||||
# 1,1,0; ...
|
||||
# -1,0,-1];
|
||||
#A0=diag(diag(A)); ## The diagonal part of A
|
||||
|
||||
## A uniform distributed random square (5x5) integer matrix with entries in [-1,1]
|
||||
do
|
||||
A=randi([-1,1],5,5);
|
||||
#A=rand(5);
|
||||
A0=diag(diag(A)); ## The diagonal part of A
|
||||
until (det(A)!=0 && det(A0)!=0) ## We need both matrices to be simultaniously non-singular
|
||||
A
|
||||
printf("Determinant of A is: %d\n",det(A))
|
||||
printf("Determinant of A0 is: %d\n\n",det(A0))
|
||||
|
||||
Ar=A-A0; ## The remainder of A
|
||||
nCols=columns(Ar); ## The number of coluns of A (M in accompanying PDF)
|
||||
Id=eye(nCols); ## Identity matrix, used for the Cartesian basis vectors
|
||||
Ainv=zeros(nCols,nCols,nCols+1); ## 3d matrix to store intermediate inverse of A
|
||||
U=zeros(nCols,nCols,nCols); ## 3d matrix to store the nCols rank-1 updates
|
||||
a=zeros(1,nCols); ## Vector containing the break-down values
|
||||
used=zeros(1,nCols); ## Vector to keep track of which updates have been used
|
||||
|
||||
## Loop to calculate A0_inv and the rank-1 updates and store in U
|
||||
Ainv(:,:,1)=eye(nCols);
|
||||
for i=1:nCols
|
||||
Ainv(i,i,1)=1/A0(i,i);
|
||||
U(:,:,i)=Ar(:,i)*transpose(Id(:,i));
|
||||
endfor
|
||||
|
||||
## Here starts the calculation of the inverse of A
|
||||
for i=1:nCols ## Outer loop iterated over the intermediates A_l^-1, M times in total
|
||||
|
||||
## Here the break-down values are calculated for each intermediate inverse
|
||||
for j=1:nCols
|
||||
a(j)=1 + transpose(Id(:,j)) * Ainv(:,:,i) * Ar(:,j); ## First time all elmts 1 because A0_inv is diagonal
|
||||
printf("a(%d) = %f\n",j,a(j))
|
||||
## Select next update ONLY if break-down value != 0 AND not yet used
|
||||
if (a(j)!=0 && used(j)!=1);
|
||||
k=j;
|
||||
break;
|
||||
endif
|
||||
endfor
|
||||
|
||||
## Do the actual S-M update
|
||||
Ainv(:,:,i+1) = Ainv(:,:,i) - Ainv(:,:,i) * U(:,:,k) * Ainv(:,:,i) / a(k);
|
||||
used(k)=1; ## Mark this update as used
|
||||
endfor
|
||||
|
||||
## Test if the inverse found is really an inverse (does not work if values are floats)
|
||||
if (A*Ainv(:,:,nCols+1)==eye(nCols))
|
||||
printf("\n");
|
||||
printf("Inverse found. A^{-1}:\n");
|
||||
Ainv(:,:,nCols+1)
|
||||
else
|
||||
printf("\n");
|
||||
printf("Inverse NOT found! A*A^{-1}:\n");
|
||||
A*Ainv(:,:,nCols+1)
|
||||
endif
|
||||
112
SMMaponiA3.m
Normal file
112
SMMaponiA3.m
Normal file
@ -0,0 +1,112 @@
|
||||
## Algorithm 3 from P. Maponi,
|
||||
## p. 283, doi:10.1016/j.laa.2006.07.007
|
||||
clc ## Clear the screen
|
||||
|
||||
## Define the matrix to be inverted. This is example 8 from the paper
|
||||
## In the future this matrix needs to be read from the function call arguments
|
||||
A=[1,1,-1; ...
|
||||
1,1,0; ...
|
||||
-1,0,-1];
|
||||
A0=diag(diag(A)); ## The diagonal part of A
|
||||
|
||||
### The modified example that gives all singular updates at some point
|
||||
#A=[1,1,1; ...
|
||||
# 1,1,0; ...
|
||||
# -1,0,-1];
|
||||
#A0=diag(diag(A)); ## The diagonal part of A
|
||||
|
||||
### A square uniform distributed random integer matrix with entries in [-1,1]
|
||||
#do
|
||||
# A=randi([-1,1],3,3);
|
||||
# A0=diag(diag(A)); ## The diagonal part of A
|
||||
#until (det(A)!=0 && det(A0)!=0) ## We need both matrices to be simultaniously non-singular
|
||||
|
||||
### A square uniform distributed random float matrix with entries in (0,1)
|
||||
#do
|
||||
# A=rand(5);
|
||||
# A0=diag(diag(A)); ## The diagonal part of A
|
||||
#until (det(A)!=0 && det(A0)!=0) ## We need both matrices to be simultaniously non-singular
|
||||
|
||||
Ar=A-A0; ## The remainder of A
|
||||
nCols=columns(Ar); ## The number of coluns of A (M in accompanying PDF)
|
||||
Id=eye(nCols);
|
||||
A0inv=eye(nCols);
|
||||
Ainv=zeros(nCols,nCols);
|
||||
ylk=zeros(nCols,nCols,nCols);
|
||||
p=zeros(nCols,1);
|
||||
breakdown=zeros(nCols,1);
|
||||
|
||||
A,A0
|
||||
printf("Determinant of A is: %d\n",det(A))
|
||||
printf("Determinant of A0 is: %d\n",det(A0))
|
||||
|
||||
## Calculate the inverse of A0 and populate p-vector
|
||||
for i=1:nCols
|
||||
A0inv(i,i) = 1 / A0(i,i);
|
||||
p(i)=i;
|
||||
endfor
|
||||
|
||||
## Calculate all the y0k in M^2 multiplications instead of M^3
|
||||
for k=1:nCols
|
||||
for i=1:nCols
|
||||
#printf("(i,k,1) = (%d,%d,1)\n",i,k);
|
||||
ylk(i,k,1) = A0inv(i,i) * Ar(i,k);
|
||||
endfor
|
||||
endfor
|
||||
|
||||
## Calculate all the ylk from the y0k calculated previously
|
||||
for l=2:nCols
|
||||
## Calculate break-down conditions and put in a vector
|
||||
for j=l-1:nCols
|
||||
breakdown(j) = abs(1+ylk(p(j),p(j),l-1));
|
||||
#printf("|1 + ylk(%d,%d,%d)|\n", p(j), p(j), l-1);
|
||||
endfor
|
||||
[val, lbar] = max(breakdown); ## Find the index of the max value
|
||||
breakdown=zeros(nCols,1); ## Reset the entries to zero for next l-round
|
||||
## Swap p(l) and p(lbar)
|
||||
tmp=p(l-1);
|
||||
p(l-1)=p(lbar);
|
||||
p(lbar)=tmp;
|
||||
for k=l:nCols
|
||||
for i=1:nCols
|
||||
ylk(i,p(k),l) = ylk(i,p(k),l-1) - (ylk(p(l-1),p(k),l-1)) / (1+ylk(p(l-1),p(l-1),l-1)) * (ylk(i,p(l-1),l-1));
|
||||
#printf("ylk(%d,%d,%d) = ylk(%d,%d,%d) - (ylk(%d,%d,%d) / (1+ylk(%d,%d,%d) * (ylk(%d,%d,%d);\n", i,p(k),l,i,p(k),l-1,p(l-1),p(k),l-1,p(l-1),p(l-1),l-1,i,p(l-1),l-1);
|
||||
endfor
|
||||
endfor
|
||||
endfor
|
||||
|
||||
## Construct A-inverse from A0-inverse and the ylk
|
||||
Ainv=A0inv;
|
||||
for l=1:nCols
|
||||
Ainv=(Id - ylk(:,p(l),l) * transpose(Id(:,p(l))) / (1 + ylk(p(l),p(l),l))) * Ainv;
|
||||
#printf("Ainv=(Id - ylk(:,%d,%d) * transpose(Id(:,%d)) / (1 + ylk(%d,%d,%d))) * Ainv\n",p(l),l,p(l),p(l),p(l),l);
|
||||
endfor
|
||||
|
||||
## Test if the inverse found is really an inverse (does not work if values are floats)
|
||||
IdTest=A*Ainv;
|
||||
if (IdTest==eye(nCols))
|
||||
printf("\n");
|
||||
printf("Inverse of A^{-1} FOUND!\n");
|
||||
Ainv
|
||||
else
|
||||
printf("\n");
|
||||
printf("Inverse of A^{-1} NOT found yet.\nRunning another test...\n");
|
||||
for i=1:nCols
|
||||
for j=1:nCols
|
||||
if (abs(IdTest(i,j))<cutOff)
|
||||
IdTest(i,j)=0;
|
||||
elseif (abs(IdTest(i,j))-1<cutOff)
|
||||
IdTest(i,j)=1;
|
||||
endif
|
||||
endfor
|
||||
endfor
|
||||
if (IdTest==eye(nCols))
|
||||
printf("\n");
|
||||
printf("Inverse of A^{-1} FOUND!\n");
|
||||
Ainv
|
||||
else
|
||||
printf("\n");
|
||||
printf("Still not found. Giving up!\n");
|
||||
IdTest
|
||||
endif
|
||||
endif
|
||||
135
SMMaponiA4.m
Normal file
135
SMMaponiA4.m
Normal file
@ -0,0 +1,135 @@
|
||||
## Algorithm 4 from P. Maponi,
|
||||
## p. 283, doi:10.1016/j.laa.2006.07.007
|
||||
clc ## Clear the screen
|
||||
|
||||
% ## Define the matrix to be inverted. This is example 8 from the paper
|
||||
% ## In the future this matrix needs to be read from the function call arguments
|
||||
% A=[1,1,-1; ...
|
||||
% 1,1,0; ...
|
||||
% -1,0,-1];
|
||||
|
||||
## The modified example that gives all singular updates at some point
|
||||
A=[1,1,1; ...
|
||||
1,1,0; ...
|
||||
-1,0,-1];
|
||||
|
||||
%## A square uniform distributed random integer matrix with entries in [-1,1]
|
||||
%do
|
||||
% A=randi([-5,5],4,4);
|
||||
%until (det(A)!=0) ## We neec matrix non-singular
|
||||
|
||||
% ## A square uniform distributed random float matrix with entries in (0,1)
|
||||
% do
|
||||
% A=rand(5);
|
||||
% until (det(A)!=0) ## We need matrix non-singular
|
||||
|
||||
|
||||
nCols=columns(A); ## The number of coluns of A (M in accompanying PDF)
|
||||
Id=eye(nCols);
|
||||
d=0.1
|
||||
A0=d*eye(nCols);
|
||||
Ar=A-A0; ## The remainder of A
|
||||
A0inv=eye(nCols);
|
||||
Ainv=zeros(nCols,nCols);
|
||||
ylk=zeros(nCols,nCols,nCols);
|
||||
breakdown=zeros(nCols,1);
|
||||
cutOff=1e-10;
|
||||
|
||||
## Calculate the inverse of A0 and populate p-vector
|
||||
for i=1:nCols
|
||||
A0inv(i,i) = 1 / A0(i,i);
|
||||
p(i)=i;
|
||||
endfor
|
||||
A,A0,Ar,A0inv
|
||||
|
||||
printf("Determinant of A is: %d\n",det(A))
|
||||
printf("Determinant of A0 is: %d\n",det(A0))
|
||||
printf("Determinant of A0inv is: %d\n",det(A0inv))
|
||||
|
||||
|
||||
|
||||
## Calculate all the y0k in M^2 multiplications instead of M^3
|
||||
for k=1:nCols
|
||||
for i=1:nCols
|
||||
ylk(i,k,1) = A0inv(i,i) * Ar(i,k);
|
||||
% printf("ylk(%d,%d,1) = A0inv(%d,%d) * Ar(%d,%d)\n",i,k,i,i,i,k);
|
||||
endfor
|
||||
endfor
|
||||
|
||||
|
||||
|
||||
## Calculate all the ylk from the y0k calculated previously
|
||||
for l=2:3#nCols
|
||||
el=Id(:,l-1);
|
||||
## Calculate break-down conditions and put in a vector
|
||||
for j=l-1:nCols
|
||||
#breakdown(j) = abs( el(j) + ylk(j,l-1,l-1) )
|
||||
breakdown(j) = abs( el(j) + ylk(j,j,l-1) )
|
||||
% printf("|el(%d) + ylk(%d,%d,%d)|\n", j, j, l-1, l-1);
|
||||
endfor
|
||||
[val, s] = max(breakdown); ## Find the index of the max value
|
||||
printf("l = %d\ns = %d\n",l-1,s);
|
||||
breakdown=zeros(nCols,1); ## Reset the entries to zero for next l-round
|
||||
if (s!=l-1) ## Apply partial pivoting
|
||||
## Swap yl-1,k(r) and yl-1,k(s) for all k=l,l+1,...,M
|
||||
r=l-1
|
||||
for k=l-1:nCols
|
||||
tmp=ylk(r,k,l-1)
|
||||
ylk(r,k,l-1)=ylk(s,k,l-1);
|
||||
printf("ylk(%d,%d,%d)=ylk(%d,%d,%d)\n",r,k,l-1,s,k,l-1);
|
||||
ylk(s,k,l-1)=tmp;
|
||||
endfor
|
||||
## Modify yl-1,r and yl-1,s
|
||||
er=Id(:,r);
|
||||
es=Id(:,s);
|
||||
ylk(:,r,l-1) = ylk(:,r,l-1) + es - er;
|
||||
ylk(:,s,l-1) = ylk(:,s,l-1) + er - es;
|
||||
endif
|
||||
## Compute finally the yl,k
|
||||
for k=l:nCols
|
||||
for i=1:nCols
|
||||
ylk(i,k,l) = ylk(i,k,l-1) - ylk(l-1,k,l-1) / (1 + ylk(l-1,l-1,l-1)) * ylk(i,l-1,l-1);
|
||||
printf("ylk(%d,%d,%d) = ylk(%d,%d,%d) - (ylk(%d,%d,%d) / (1+ylk(%d,%d,%d) * (ylk(%d,%d,%d)\n",i,k,l,i,k,l-1,l-1,k,l-1,l-1,l-1,l-1,i,l-1,l-1);
|
||||
endfor
|
||||
endfor
|
||||
endfor
|
||||
|
||||
|
||||
|
||||
## Construct A-inverse from A0-inverse and the ylk
|
||||
Ainv=A0inv;
|
||||
for l=1:nCols
|
||||
Ainv=(Id - ylk(:,l,l) * transpose(Id(:,l)) / (1 + ylk(l,l,l))) * Ainv;
|
||||
% printf("Ainv=(Id - ylk(:,%d,%d) * transpose(Id(:,%d)) / (1 + ylk(%d,%d,%d))) * Ainv\n",l,l,l,l,l,l);
|
||||
endfor
|
||||
|
||||
|
||||
|
||||
## Test if the inverse found is really an inverse (does not work if values are floats)
|
||||
IdTest=A*Ainv
|
||||
if (IdTest==eye(nCols))
|
||||
printf("\n");
|
||||
printf("Inverse of A^{-1} FOUND!\n");
|
||||
Ainv
|
||||
else
|
||||
printf("\n");
|
||||
printf("Inverse of A^{-1} NOT found yet.\nRunning another test...\n");
|
||||
for i=1:nCols
|
||||
for j=1:nCols
|
||||
if (abs(IdTest(i,j))<cutOff)
|
||||
IdTest(i,j)=0;
|
||||
elseif (abs(IdTest(i,j))-1<cutOff)
|
||||
IdTest(i,j)=1;
|
||||
endif
|
||||
endfor
|
||||
endfor
|
||||
if (IdTest==eye(nCols))
|
||||
printf("\n");
|
||||
printf("Inverse of A^{-1} FOUND!\n");
|
||||
Ainv
|
||||
else
|
||||
printf("\n");
|
||||
printf("Still not found. Giving up!\n");
|
||||
Ainv,IdTest
|
||||
endif
|
||||
endif
|
||||
Loading…
Reference in New Issue
Block a user