Inverse computation with a starting guess

I am runinng that uses a for loop where in each iteration, an inverse of a 1000x1000 matrix is being computed which is slowing the loop. It is not possible to use A\b form since i explicitly need the inverse matrix, so I am using inv(A).

But there is an additional piece of information, that is, in each for loop iteration, i am changing a parameter slightly, that alters A by small amounts, hence the inverse in iteration should be close to the inverse in previous iteration. Is it possible to provide the last inverse as a starting point or guess for the inverse computation in matlab? I dont see such a possiblity using inv(A), but if it can be done, it will save me a lot of time.

Answer 1

[John D'Errico]

It is often the case that someone THINKs they need a matrix inverse, yet they don't. And since you don't know if you can do this computation, that makes it entirely more likely you may not actually really need a matrix inverse.

Anyway, having a previous matrix inverse will not generally help you, although that too is difficult to know. For example IF the change from one matrix to the next was a simple rank 1 update, then yes, there is a tool that can be used, the Sherman-Morrison-Woodbury formula.

https://www.cs.cornell.edu/~bindel/class/cs6210-f09/lec12.pdf

Depending if the change is a rank 1 update, or a higher rank update, there are variations on the formula.

But honestly, if you won't tell us what the perturbations to your matrix entail, and you won't tell us why it is that you think you need the matrix inverse, is there a good way for us to know how to help you?

[Shlok Vaibhav Singh]

Hi John, sorry couldnt follow up due to exams.

I am working with green's function in nanoscale devices, and to get current value throught the device, I have to take trace of Green's function, which is computed as the inverse of (hamiltonian - energy ) of the system. I dont think there is a way to avoid the inverse operator.

A trick that might avoid inverse computation is:

so that we can have as current but this involves two \ operators.

And A is always full rank, G can always be found

A sample code is (Taken from Datta, Atom to Transistor)

clear all

figure

%Constants (all MKS, except energy which is in eV)

hbar=1.06e-34;q=1.6e-19;m=.25*9.1e-31;IE=(q*q)/(2*pi*hbar);

Ef=0.1;kT=.025;

%inputs

a=3e-10;t0=(hbar^2)/(2*m*(a^2)*q);

NS=15;NC=30;ND=15;Np=NS+NC+ND;

%Hamiltonian matrix

%NS=15;NC=20;ND=15;Np=NS+NC+ND;UB=0*ones(Np,1);%no barrier

NS=23;NC=4;ND=23;Np=NS+NC+ND;

UB=[zeros(NS,1);0.4*ones(NC,1);zeros(ND,1);];%tunneling barrier

%NS=15;NC=16;ND=15;Np=NS+NC+ND;

%UB=[zeros(NS,1);.4*ones(4,1);zeros(NC-8,1);.4*ones(4,1);zeros(ND,1)];%RT barrier

T=(2*t0*diag(ones(1,Np)))-(t0*diag(ones(1,Np-1),1))-(t0*diag(ones(1,Np-1),-1));

T=T+diag(UB);

%Bias

V=0;mu1=Ef+(V/2);mu2=Ef-(V/2);

U1=V*[.5*ones(1,NS) linspace(0.5,-0.5,NC) -.5*ones(1,ND)];

U1=0*U1; %U1';%Applied potential profile

%Energy grid for Green’s function method

NE=61;E=linspace(-.2,.8,NE);zplus=i*1e-12;dE=E(2)-E(1);

f1=1./(1+exp((E-mu1)./kT));

f2=1./(1+exp((E-mu2)./kT));

%Transmission

I=0;%Current

for k=1:NE

sig1=zeros(Np);sig2=zeros(Np);sig3=zeros(Np);

ck=1-((E(k)+zplus-U1(1)-UB(1))/(2*t0));ka=acos(ck);

sig1(1,1)=-t0*exp(i*ka);gam1=i*(sig1-sig1');

ck=1-((E(k)+zplus-U1(Np)-UB(Np))/(2*t0));ka=acos(ck);

sig2(Np,Np)=-t0*exp(i*ka);gam2=i*(sig2-sig2');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%% G being computed by inversion here %%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

G=inv(((E(k)+zplus)*eye(Np))-T-diag(U1)-sig1-sig2-sig3);

TM(k)=real(trace(gam1*G*gam2*G'));

I=I+(dE*IE*TM(k)*(f1(k)-f2(k)));

end

V,I

XX=a*1e9*[1:1:Np];

XS=XX([1:NS-4]);XD=XX([NS+NC+5:Np]);

hold on

h = plot( E,TM, 'b')

%h = plot(XX,U1+UB,'b');

%h=plot(XS,mu1*ones(1,NS-4),'b--');

%h=plot(XD,mu2*ones(1,ND-4),'b--');

%axis([0 1.1 -.2 .8])

%axis([0 15 -.2 .8])

set(h,'linewidth',[2.0])

set(gca,'Fontsize',[24])

xlabel(' z ( nm ) --->')

ylabel(' Energy ( eV ) ---> ')

grid on