Is python sufficient for serious tensor network or DMRG calculations?

Question

There are enough tensor-network/DMRG packages that support python (ALPS, OpenMPS, pyUni10) that it would seem (to an outsider) that it's possible to do some substantial work with DMRG without messing with C or any other languages. Are these packages sufficiently fast and full-featured that students could do meaningful research work using just python?

Answer 1

There probably are good Python bindings to various DMRG implementations, which allow one to run DMRG from Python. Since the implementations typically rely on lower-level C/C++/Fortran routines, the calculations run quite quickly. E.g. PySCF appears to have bindings to various DMRG programs, see https://sunqm.github.io/pyscf/dmrgscf.html.

If you're talking about developing a new Python-only DMRG code, I doubt you'll succeed: the overhead is just too great. Julia might be an interesting alternative, though.

Answer 2

I would argue that you can use a package like DMRG++, which is written in C++, without knowing much C/C++ or indeed any language. It has its own input format that takes some getting used to (unfortunately it isn't overly well documented), but the outputs are provided in text and hdf5 formats that you can then manipulate with your favorite scripting language - e.g. python. Yet all heavy lifting, so to speak, is done by the provided C++ code. Unless you need to modify or troubleshoot that code, perhaps to add a non-standard model, you can definitively avoid messing with C++. Also, the code is actively maintained so even in such cases it might be possible to avoid messing with C++ yourself.

Answer 3

I've heard people talking about working on tensor networks for electronic structure since the 2000s, but here in the 2020s I'm aware of no code that can help me solve the problems in which I'm interested. I think "serious" tensor network codes for electronic structure are still a "dream", and I'm becoming more and more doubtful that any time in the next few years we'll see (more general) tensor networks become "mainstream" in (almost) the way DMRG did.

I have heard of progress with PEPS and MERA but mainly for spin Hamiltonians. I see that the Wikipedia page for DMRG currently provides links to iTensor, and Uni10 for which Adam Iaizzi (one of our favorite people, here at Matter Modeling SE) is a co-author. I can't quite figure out based on the documentation, how general of a tensor network iTensor supports, but it seems that whatever tensor networks they implement, can only be applied to spin Hamiltonians "models" rather than the "electronic structure Hamiltonian" in which I am personally interested. Uni10 seems to be able to do PEPS and MERA in addition to DMRG and maybe more, but again I cannot figure out to which Hamiltonians it can be applied. Both iTensor and Uni10 are written in C++.

"Is python sufficient for serious tensor network/DMRG calculations?"

In 2011 Martin Plenio told me at a Royal Society meeting that he coded a basic DMRG program himself in MATLAB (which you can imagine, is very doable). This means that NumPy too, can suffice for some purposes, since MATLAB and NumPy both call BLAS to do the heavy lifting (matrix arithmetic).

I wrote a very similar code called FeynDyn (Feynman Dynamics) for calcualting Feynman integrals numerically using the "tensor propagation algorithm", which involves extremely high-rank tensors, organized in such a way that all necessary calculations can be done using matrix arithmetic (so the whole program is written in MATLAB which is about as fast as you can get for matrix arithmetic on either CPUs or GPUs). This is as far as I know, the fastest tensor propagation code for Feynman integrals based on all comparisons I've done to FORTRAN codes written by leading colleagues in the field: Not only does it get the calculations done faster than all FORTRAN codes I've compared it to, when I run both my code and their code on a CPU, but thanks to the GPU programmers working at MathWorks my code benefits from a major speed-up, while none of the FORTRAN or C++ codes do can run on GPUs (they could probably get their codes working on GPUs with OpenACC or some CUDA wrapper but only someone who knows their code well enough can do it, and based on what I understand, those codes were not vectorized like mine, so they'll probably have to write the code from scratch to compete with my GPU performance). I had never heard of NumPy when I started writing FeynDyn, but my understanding is that it could be similar in performance but might not get the same GPU acceleration (at least not nearly as easily!).

Closing thoughts:

• If the fastest tensor propagation code for Feynman integrals is written in MATLAB, and beats the FORTRAN codes written by other world leaders in the field, then NumPy may also do reasonably well on CPUs.
• This is because the tensor propagation algorithm (not too different from what happens in a tensor network code, in terms of core low-level arithmetic used) can be vectorized so that it can pretty much be done using only BLAS functions that were written by people whose job was to make BLAS as fast as possible.
• As far as I know, Python (and Octave!) have not caught up to MATLAB in terms of GPU computing. It can be frustrating to work with MATLAB due to licensing issues and the amount of time it takes to load (it doesn't sound like the most important thing, but it's still a problem for me even in 2020). Nevertheless, even my present MSc student who was born in the mid-90s prefers MATLAB for some of the things we're doing.
• I am one of the most anti-Python people you'll meet, and when I first saw this question I immediately thought of pointing to the fact that all "serious" electronic structure codes use FORTRAN or C++, at least for the numerically expensive parts, but halfway through writing this answer I remembered the similarity between my tensor propagation code and tensor networks, and have come to the conclusion that NumPy is sufficient, but MATLAB is better for the GPU support and Julia might be the best thing in the middle.

Answer 4

For 1D spin Hamiltons it is perfectly doable even if you code from scratch. NumPy should be sufficient. I have written DMRG codes in both Julia and python and I see no difference, this is also due to the fact that I am a newbie in Julia and my codes are far from good. However, if you are willing to use python I will suggest you to do all the algorithm level optimization first, find efficient tensor contraction order for all the contractions, write your own cheap Lanczos/Arnoldo eigensolver that saves you from building big matrices to use a Python eigensolver libraries.

Currently my code for long range transverse field Ising model with 200 lattice sites and bond dimension as high as 100 can do a ground state search in just few minutes and all has been coded in python from the scratch.