Stability of the retarded Green's function calculation
Hello Kwantoptians, I have a question regarding the speed/stability of computing the retarded Green's function of a transport setup using Kwant. In the Kwant sourcecode [1] it is noted that using `kwant.greens_function` is "often slower and less stable than the scattering matrix calculation". I was wondering if you could provide me with some references for this assertion. I have perused the Kwant paper [2] and (part I of) the thesis of Michael Wimmer [3] and cannot find any mention of a speed/stability comparison of the algorithm implemented by `kwant.greens_function` vs `kwant.smatrix`. My understanding is that in both cases a linear system (LHS) is constructed and solved for different right handsides (RHS) using out of the box linear solvers. The solution for each RHS corresponds to one column of the retarded Green's function / extended scattering matrix respectively. The difference between `kwant.greens_function` and `kwant.smatrix` is then the following. In the former case the leads are taken into account by added the retarded selfenergy to the LHS and the RHSs are indicator vectors for the sites of the lead/scattering region interface, which corresponds to a "unit impulse" boundary condition. In the latter case the linear system is "extended" so as to include extra unknowns that correspond to the scattering amplitudes, and the RHSs are indicator vectors with the 1's in this "extended" part, which corresponds to a "single incoming mode" boundary condition. It seems to me that the salient difference is in the boundary conditions, and I do not have a good intuition as to why one set of boundary conditions would make the linear system easier/harder to solve. Happy Kwanting! Joe [1]: https://gitlab.kwantproject.org/kwant/kwant//blob/master/kwant/solvers/com... [2]: https://iopscience.iop.org/article/10.1088/13672630/16/6/063065/pdf [3]: https://epub.uniregensburg.de/12142/
Hello again, I Noticed a couple of minor mistakes in my previous email: 1. I used the term "out of the box" linear solvers, when I meant "black box" linear solvers. By this I meant that we are not tailoring the algorithm used for solving the linear system based on the properties of the LHS (e.g. as we do in RGF by the choice of scattering region slices), but rather trusting that the applied mathematicians have done their job properly, and any incidental structure in the LHS is found "automatically" by the solver. 2. I claimed that the RHSs for the scattering problem were "indicator vectors with 1s in [the] extended part". This is not true, however the interpretation of the RHS as a "single incoming mode" is correct AFAIK. Happy Kwanting, Joe From: Kwantdiscuss <kwantdiscussbounces@kwantproject.org> On Behalf Of Joseph Weston (Aquent LLC  Canada) Sent: Wednesday, March 18, 2020 12:08 PM To: kwantdiscuss@kwantproject.org Subject: [EXTERNAL] [Kwant] Stability of the retarded Green's function calculation Hello Kwantoptians, I have a question regarding the speed/stability of computing the retarded Green's function of a transport setup using Kwant. In the Kwant sourcecode [1] it is noted that using `kwant.greens_function` is "often slower and less stable than the scattering matrix calculation". I was wondering if you could provide me with some references for this assertion. I have perused the Kwant paper [2] and (part I of) the thesis of Michael Wimmer [3] and cannot find any mention of a speed/stability comparison of the algorithm implemented by `kwant.greens_function` vs `kwant.smatrix`. My understanding is that in both cases a linear system (LHS) is constructed and solved for different right handsides (RHS) using out of the box linear solvers. The solution for each RHS corresponds to one column of the retarded Green's function / extended scattering matrix respectively. The difference between `kwant.greens_function` and `kwant.smatrix` is then the following. In the former case the leads are taken into account by added the retarded selfenergy to the LHS and the RHSs are indicator vectors for the sites of the lead/scattering region interface, which corresponds to a "unit impulse" boundary condition. In the latter case the linear system is "extended" so as to include extra unknowns that correspond to the scattering amplitudes, and the RHSs are indicator vectors with the 1's in this "extended" part, which corresponds to a "single incoming mode" boundary condition. It seems to me that the salient difference is in the boundary conditions, and I do not have a good intuition as to why one set of boundary conditions would make the linear system easier/harder to solve. Happy Kwanting! Joe [1]: https://gitlab.kwantproject.org/kwant/kwant//blob/master/kwant/solvers/common.py#L428<https://nam06.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgitlab.kwantproject.org%2Fkwant%2Fkwant%2F%2Fblob%2Fmaster%2Fkwant%2Fsolvers%2Fcommon.py%23L428&data=02%7C01%7Cvjosewe%40microsoft.com%7Cccd4df85ea5a44fb98e508d7cb6fbfa3%7C72f988bf86f141af91ab2d7cd011db47%7C1%7C0%7C637201553194781503&sdata=jNV1lZwk4xs%2BfX25KOuP0Z992wS3as9Uribwou1P1rk%3D&reserved=0> [2]: https://iopscience.iop.org/article/10.1088/13672630/16/6/063065/pdf<https://nam06.safelinks.protection.outlook.com/?url=https%3A%2F%2Fiopscience.iop.org%2Farticle%2F10.1088%2F13672630%2F16%2F6%2F063065%2Fpdf&data=02%7C01%7Cvjosewe%40microsoft.com%7Cccd4df85ea5a44fb98e508d7cb6fbfa3%7C72f988bf86f141af91ab2d7cd011db47%7C1%7C0%7C637201553194791491&sdata=5ltAdppLxjv%2BAf6fHEb4xTlSrZ%2FcGNY1jPFwnWO932s%3D&reserved=0> [3]: https://epub.uniregensburg.de/12142/<https://nam06.safelinks.protection.outlook.com/?url=https%3A%2F%2Fepub.uniregensburg.de%2F12142%2F&data=02%7C01%7Cvjosewe%40microsoft.com%7Cccd4df85ea5a44fb98e508d7cb6fbfa3%7C72f988bf86f141af91ab2d7cd011db47%7C1%7C0%7C637201553194791491&sdata=B1HTfGTAeQqmX%2BuLFsgxrGdtMYMij%2BbDkgupLe6efOg%3D&reserved=0>
Hi Joe, Computing the Green's function or the selfenergy of the lead fails when a terminated lead has a bound state at the energy we're looking at (so there's an exact pole). For example a lead made out of a topological superconductor always has this problem at zero energy. At the same time the mode solvers face no singularities in this case. In linear algebra terms computing a selfenergy corresponds to partially fixing the order of Gaussian elimination, while the most general case may require pivoting incompatible with this order. That's what concerns stability. The speed difference is, I think, insignificant: the number of right hand sides is smaller with the modes solver since we use one per incoming mode, not one per degree of freedom. However in practice obtaining the LU decomposition causes the largest slowdown. Cheers, Anton On Wed, 18 Mar 2020 at 22:12, Joseph Weston (Aquent LLC  Canada) <vjosewe@microsoft.com> wrote:
Hello again,
I Noticed a couple of minor mistakes in my previous email:
I used the term “out of the box” linear solvers, when I meant “black box” linear solvers. By this I meant that we are not tailoring the algorithm used for solving the linear system based on the properties of the LHS (e.g. as we do in RGF by the choice of scattering region slices), but rather trusting that the applied mathematicians have done their job properly, and any incidental structure in the LHS is found “automatically” by the solver. I claimed that the RHSs for the scattering problem were “indicator vectors with 1s in [the] extended part”. This is not true, however the interpretation of the RHS as a “single incoming mode” is correct AFAIK.
Happy Kwanting,
Joe
From: Kwantdiscuss <kwantdiscussbounces@kwantproject.org> On Behalf Of Joseph Weston (Aquent LLC  Canada) Sent: Wednesday, March 18, 2020 12:08 PM To: kwantdiscuss@kwantproject.org Subject: [EXTERNAL] [Kwant] Stability of the retarded Green's function calculation
Hello Kwantoptians,
I have a question regarding the speed/stability of computing the retarded Green’s function of a transport setup using Kwant.
In the Kwant sourcecode [1] it is noted that using `kwant.greens_function` is “often slower and less stable than the scattering matrix calculation”. I was wondering if you could provide me with some references for this assertion.
I have perused the Kwant paper [2] and (part I of) the thesis of Michael Wimmer [3] and cannot find any mention of a speed/stability comparison of the algorithm implemented by `kwant.greens_function` vs `kwant.smatrix`.
My understanding is that in both cases a linear system (LHS) is constructed and solved for different right handsides (RHS) using out of the box linear solvers. The solution for each RHS corresponds to one column of the retarded Green’s function / extended scattering matrix respectively. The difference between `kwant.greens_function` and `kwant.smatrix` is then the following. In the former case the leads are taken into account by added the retarded selfenergy to the LHS and the RHSs are indicator vectors for the sites of the lead/scattering region interface, which corresponds to a “unit impulse” boundary condition. In the latter case the linear system is “extended” so as to include extra unknowns that correspond to the scattering amplitudes, and the RHSs are indicator vectors with the 1’s in this “extended” part, which corresponds to a “single incoming mode” boundary condition.
It seems to me that the salient difference is in the boundary conditions, and I do not have a good intuition as to why one set of boundary conditions would make the linear system easier/harder to solve.
Happy Kwanting!
Joe
[1]: https://gitlab.kwantproject.org/kwant/kwant//blob/master/kwant/solvers/com...
[2]: https://iopscience.iop.org/article/10.1088/13672630/16/6/063065/pdf
Hi Anton, Thanks for the detailed response. I think I understand most of your answer, but I'm a little confused about one part. If I remember correctly, Kwant computes selfenergies for leads using the mode decomposition via V†ΦΛΦ^1, where V is the hopping, Φ Is the matrix of outgoing modes (Φ^1 its rightinverse) and Λ is the matrix of translation eigenvalues. The only places I can see where a pole may appear infinite translation eigenvalues λ or badly conditioned Φ. On the other hand, when constructing the linear system for the scattering problem we extend the lefthandside with terms like V†ΦΛ. When you say "computing a selfenergy corresponds to partially fixing the order of Gaussian elimination" do you mean that computing the selfenergy forces us to invert Φ, and *this* is what leads to an instability? As I write this I realise that I am just restating what you said, but I would like a bit of clarity on this point. Thanks, Joe Original Message From: Anton Akhmerov <anton.akhmerov+kd@gmail.com> Sent: Sunday, March 22, 2020 8:41 AM To: Joseph Weston (Aquent LLC  Canada) <vjosewe@microsoft.com> Cc: kwantdiscuss@kwantproject.org Subject: [EXTERNAL] Re: [Kwant] Stability of the retarded Green's function calculation Hi Joe, Computing the Green's function or the selfenergy of the lead fails when a terminated lead has a bound state at the energy we're looking at (so there's an exact pole). For example a lead made out of a topological superconductor always has this problem at zero energy. At the same time the mode solvers face no singularities in this case. In linear algebra terms computing a selfenergy corresponds to partially fixing the order of Gaussian elimination, while the most general case may require pivoting incompatible with this order. That's what concerns stability. The speed difference is, I think, insignificant: the number of right hand sides is smaller with the modes solver since we use one per incoming mode, not one per degree of freedom. However in practice obtaining the LU decomposition causes the largest slowdown. Cheers, Anton On Wed, 18 Mar 2020 at 22:12, Joseph Weston (Aquent LLC  Canada) <vjosewe@microsoft.com> wrote:
Hello again,
I Noticed a couple of minor mistakes in my previous email:
I used the term “out of the box” linear solvers, when I meant “black box” linear solvers. By this I meant that we are not tailoring the algorithm used for solving the linear system based on the properties of the LHS (e.g. as we do in RGF by the choice of scattering region slices), but rather trusting that the applied mathematicians have done their job properly, and any incidental structure in the LHS is found “automatically” by the solver. I claimed that the RHSs for the scattering problem were “indicator vectors with 1s in [the] extended part”. This is not true, however the interpretation of the RHS as a “single incoming mode” is correct AFAIK.
Happy Kwanting,
Joe
From: Kwantdiscuss <kwantdiscussbounces@kwantproject.org> On Behalf Of Joseph Weston (Aquent LLC  Canada) Sent: Wednesday, March 18, 2020 12:08 PM To: kwantdiscuss@kwantproject.org Subject: [EXTERNAL] [Kwant] Stability of the retarded Green's function calculation
Hello Kwantoptians,
I have a question regarding the speed/stability of computing the retarded Green’s function of a transport setup using Kwant.
In the Kwant sourcecode [1] it is noted that using `kwant.greens_function` is “often slower and less stable than the scattering matrix calculation”. I was wondering if you could provide me with some references for this assertion.
I have perused the Kwant paper [2] and (part I of) the thesis of Michael Wimmer [3] and cannot find any mention of a speed/stability comparison of the algorithm implemented by `kwant.greens_function` vs `kwant.smatrix`.
My understanding is that in both cases a linear system (LHS) is constructed and solved for different right handsides (RHS) using out of the box linear solvers. The solution for each RHS corresponds to one column of the retarded Green’s function / extended scattering matrix respectively. The difference between `kwant.greens_function` and `kwant.smatrix` is then the following. In the former case the leads are taken into account by added the retarded selfenergy to the LHS and the RHSs are indicator vectors for the sites of the lead/scattering region interface, which corresponds to a “unit impulse” boundary condition. In the latter case the linear system is “extended” so as to include extra unknowns that correspond to the scattering amplitudes, and the RHSs are indicator vectors with the 1’s in this “extended” part, which corresponds to a “single incoming mode” boundary condition.
It seems to me that the salient difference is in the boundary conditions, and I do not have a good intuition as to why one set of boundary conditions would make the linear system easier/harder to solve.
Happy Kwanting!
Joe
[1]: https://nam06.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgitl ab.kwantproject.org%2Fkwant%2Fkwant%2F%2Fblob%2Fmaster%2Fkwant%2Fsol vers%2Fcommon.py%23L428&data=02%7C01%7Cvjosewe%40microsoft.com%7C 636fd73dc8614ac7a0c708d7ce776cf0%7C72f988bf86f141af91ab2d7cd011db47%7C 1%7C0%7C637204884633045315&sdata=6BCOPJsvwdXfsD2Zi1huK8BjB0avZFaV9 rFijI1yUAQ%3D&reserved=0
[2]: https://nam06.safelinks.protection.outlook.com/?url=https%3A%2F%2Fiops cience.iop.org%2Farticle%2F10.1088%2F13672630%2F16%2F6%2F063065%2Fpdf &data=02%7C01%7Cvjosewe%40microsoft.com%7C636fd73dc8614ac7a0c708d 7ce776cf0%7C72f988bf86f141af91ab2d7cd011db47%7C1%7C0%7C637204884633045 315&sdata=rYQD5petvR%2FMr18hXSCs9jAhCDBhfewYxAw5XdF3qU8%3D&res erved=0
[3]: https://nam06.safelinks.protection.outlook.com/?url=https%3A%2F%2Fepub .uniregensburg.de%2F12142%2F&data=02%7C01%7Cvjosewe%40microsoft. com%7C636fd73dc8614ac7a0c708d7ce776cf0%7C72f988bf86f141af91ab2d7cd011db47%7C1%7C0%7C637204884633045315&sdata=lCRUeoPlqwsCreOR01xTCstggBLTPxOPzGOEa22GnOI%3D&reserved=0
Hey Joe,
When you say "computing a selfenergy corresponds to partially fixing the order of Gaussian elimination" do you mean that computing the selfenergy forces us to invert Φ, and *this* is what leads to an instability? As I write this I realise that I am just restating what you said, but I would like a bit of clarity on this point.
Indeed: when there's a lead bound state at the energy of interest, Φ becomes not invertible. More generally, there are no guarantees on the condition number of Φ. Cheers, Anton
Thanks,
Joe
Original Message From: Anton Akhmerov <anton.akhmerov+kd@gmail.com> Sent: Sunday, March 22, 2020 8:41 AM To: Joseph Weston (Aquent LLC  Canada) <vjosewe@microsoft.com> Cc: kwantdiscuss@kwantproject.org Subject: [EXTERNAL] Re: [Kwant] Stability of the retarded Green's function calculation
Hi Joe,
Computing the Green's function or the selfenergy of the lead fails when a terminated lead has a bound state at the energy we're looking at (so there's an exact pole). For example a lead made out of a topological superconductor always has this problem at zero energy. At the same time the mode solvers face no singularities in this case. In linear algebra terms computing a selfenergy corresponds to partially fixing the order of Gaussian elimination, while the most general case may require pivoting incompatible with this order.
That's what concerns stability. The speed difference is, I think, insignificant: the number of right hand sides is smaller with the modes solver since we use one per incoming mode, not one per degree of freedom. However in practice obtaining the LU decomposition causes the largest slowdown.
Cheers, Anton
On Wed, 18 Mar 2020 at 22:12, Joseph Weston (Aquent LLC  Canada) <vjosewe@microsoft.com> wrote:
Hello again,
I Noticed a couple of minor mistakes in my previous email:
I used the term “out of the box” linear solvers, when I meant “black box” linear solvers. By this I meant that we are not tailoring the algorithm used for solving the linear system based on the properties of the LHS (e.g. as we do in RGF by the choice of scattering region slices), but rather trusting that the applied mathematicians have done their job properly, and any incidental structure in the LHS is found “automatically” by the solver. I claimed that the RHSs for the scattering problem were “indicator vectors with 1s in [the] extended part”. This is not true, however the interpretation of the RHS as a “single incoming mode” is correct AFAIK.
Happy Kwanting,
Joe
From: Kwantdiscuss <kwantdiscussbounces@kwantproject.org> On Behalf Of Joseph Weston (Aquent LLC  Canada) Sent: Wednesday, March 18, 2020 12:08 PM To: kwantdiscuss@kwantproject.org Subject: [EXTERNAL] [Kwant] Stability of the retarded Green's function calculation
Hello Kwantoptians,
I have a question regarding the speed/stability of computing the retarded Green’s function of a transport setup using Kwant.
In the Kwant sourcecode [1] it is noted that using `kwant.greens_function` is “often slower and less stable than the scattering matrix calculation”. I was wondering if you could provide me with some references for this assertion.
I have perused the Kwant paper [2] and (part I of) the thesis of Michael Wimmer [3] and cannot find any mention of a speed/stability comparison of the algorithm implemented by `kwant.greens_function` vs `kwant.smatrix`.
My understanding is that in both cases a linear system (LHS) is constructed and solved for different right handsides (RHS) using out of the box linear solvers. The solution for each RHS corresponds to one column of the retarded Green’s function / extended scattering matrix respectively. The difference between `kwant.greens_function` and `kwant.smatrix` is then the following. In the former case the leads are taken into account by added the retarded selfenergy to the LHS and the RHSs are indicator vectors for the sites of the lead/scattering region interface, which corresponds to a “unit impulse” boundary condition. In the latter case the linear system is “extended” so as to include extra unknowns that correspond to the scattering amplitudes, and the RHSs are indicator vectors with the 1’s in this “extended” part, which corresponds to a “single incoming mode” boundary condition.
It seems to me that the salient difference is in the boundary conditions, and I do not have a good intuition as to why one set of boundary conditions would make the linear system easier/harder to solve.
Happy Kwanting!
Joe
[1]: https://nam06.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgitl ab.kwantproject.org%2Fkwant%2Fkwant%2F%2Fblob%2Fmaster%2Fkwant%2Fsol vers%2Fcommon.py%23L428&data=02%7C01%7Cvjosewe%40microsoft.com%7C 636fd73dc8614ac7a0c708d7ce776cf0%7C72f988bf86f141af91ab2d7cd011db47%7C 1%7C0%7C637204884633045315&sdata=6BCOPJsvwdXfsD2Zi1huK8BjB0avZFaV9 rFijI1yUAQ%3D&reserved=0
[2]: https://nam06.safelinks.protection.outlook.com/?url=https%3A%2F%2Fiops cience.iop.org%2Farticle%2F10.1088%2F13672630%2F16%2F6%2F063065%2Fpdf &data=02%7C01%7Cvjosewe%40microsoft.com%7C636fd73dc8614ac7a0c708d 7ce776cf0%7C72f988bf86f141af91ab2d7cd011db47%7C1%7C0%7C637204884633045 315&sdata=rYQD5petvR%2FMr18hXSCs9jAhCDBhfewYxAw5XdF3qU8%3D&res erved=0
[3]: https://nam06.safelinks.protection.outlook.com/?url=https%3A%2F%2Fepub .uniregensburg.de%2F12142%2F&data=02%7C01%7Cvjosewe%40microsoft. com%7C636fd73dc8614ac7a0c708d7ce776cf0%7C72f988bf86f141af91ab2d7cd011db47%7C1%7C0%7C637204884633045315&sdata=lCRUeoPlqwsCreOR01xTCstggBLTPxOPzGOEa22GnOI%3D&reserved=0
participants (2)

Anton Akhmerov

Joseph Weston (Aquent LLC  Canada)