Inverse computation with a starting guess

6 visualizzazioni (ultimi 30 giorni)
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.

Risposte (1)

John D'Errico
John D'Errico il 28 Mar 2023
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.
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?
  1 Commento
Shlok Vaibhav Singh
Shlok Vaibhav Singh il 17 Apr 2023
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

Accedi per commentare.

Categorie

Scopri di più su General Applications in Help Center e File Exchange

Prodotti


Release

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by