DanceQ: High-performance library for number conserving bases

Detailed discussion of efficient implementation of number-conserving basis.

No spatial symmetries.

The authors present a careful discussion and implementation of a number-conserving basis, to be used in the context of exact diagonalization. As indicated in the manuscript, this is close to methods used by practitioners in the field. Nevertheless, this is often poorly documented. It is therefore valuable to have a detailed discussion that includes generalizations and an optimized implementation strategy.

The main shortcoming that I see is the absence of spatial symmetries. For example, at the end of the discussion (i.e., on page 14) there is a comment that 120 TByte of main memory would be required for a spin-1/2 system with 46 sites. If this system has translational symmetry, it should be possible to reduce this memory requirement to below 3 TByte, i.e., a much more reasonable amount. Implementing spatial symmetries efficiently is in fact the bigger challenge than number conservation. ${\cal H}\Phi$ [41] may not use spatial symmetries either, but there are packages implementing these, such as SPINPACK, https://www-e.uni-magdeburg.de/jschulen/spin/ (already mentioned by Referee 2) and Xdiag by Alexander Wietek, https://github.com/awietek/xdiag. To be quite clear, I am not saying that inclusion of spatial symmetries is a condition for publication of the present manuscript and library, but I think that the authors should add more comments on this issue than the half sentence that they devote to it on the right column of page 1.

A more technical comment concerns the parallelization strategy outlined in sections II B and IV A. If I understood correctly, they propose to parallelize over target. I believe that parallelization over source would be the better strategy since a fixed target element avoids write conflicts. Given that the Hamiltonian is a Hermitean matrix, the two strategies are formally equivalent, but e.g. CPU cache performance can differ when looping over source or target. I encourage the authors to think about this issue, and at least reconsider their related comments.

One part that I really enjoyed reading is the historical Introduction. I nevertheless permit myself to propose complementing the well-known Ref. [24] by the older but lesser known publication by J. C. Bonner and M. E. Fisher in Proc. Phys. Soc. 80, 508 (1962).

A final detail: I believe that typesetting of the $\vert\vec{\sigma^{(0)}} \rangle$ in Eq. (11) should be corrected.

1- Expand comments on spatial symmetries.
2- Cite further exact diagonalization packages.
3- Revisit parallelization strategy.
4- Possibly add a citation of J. C. Bonner and M. E. Fisher, Proc. Phys. Soc. 80, 508 (1962).
5- Correct typesetting of $\vert \vec{\sigma^{(0)} }\rangle$ in Eq. (11).
6- I was wondering why the authors did not use the SciPost style file. This is not mandatory, but would help to appraise formatting and may speed up production.

Ask for minor revision

1- Efficient and elegant numerical library to generate and access qubit configurations with a fixed U(1) charge (e.g. magnetization for spin system)
2- Open source code
3- Very generic code able to tackle any local Hilbert space

This paper provides a numerical library to generate and access generic qubit configurations for any quantum model with a conserved U(1) quantity. This is very often the case, e.g. total Sz magnetization for a spin system, and hence relevant to most numerical studies. Due to the exponentially large Hilbert space, generating the proper configurations (among the $q^N$ where $q$ is the local dimension and $N$ the number of sites) is not straightforward. The reverse question, regarding the index number of a given configuration, is needed when computing a matrix-vector multiplication on-the-fly, as done in state-of-the-art exact diagonalization (ED) methods. Again it is nontrivial and requires lookup tables that can become quite large and limit the computations.

This numerical library is very generic and provides an implementation for both mappings between physical configurations and indices in the full list. Thus, it is a very useful tool for any ED-like code. As reminded in the paper, a solution was provided by Lin long ago by dividing the system in two parts. The main achievement here is to generalize this idea to an arbitrary decomposition as well as local Hilbert space $q$.
Basically, it is possible to divide the $N$ sites into various subsystems such that the lookup tables can fit in memory (obviously if one uses less subsystems, computations will be more efficient but this is more costly in memory).

I find the paper very well written and the results convincing.
For all these reasons, I recommend the publication in SciPost Physics Codebases.

1- In many applications, $q=2$ (qubit or spin 1/2) and it is possible to generate valid spin configurations using a Gray-like code, which should be mentioned. Moreover, since they are ordered by construction, finding the index of a given configuration can be done with a simple binomial search. Why is it written that it is "prohibitevely expensive" ? It seems to me that it could be done for large systems as well ?

2- Minor changes: there are missing capital letters in some names, e.g. "bose", "fermi", "hubbard" etc.

Publish (surpasses expectations and criteria for this Journal; among top 10%)

I think the presented manuscript could be published in SciPost Physics Codebases after revision. Overall I think that the manuscript provides valuable insight for numerically oriented research groups dealing with spin or Hubbard systems. I have a few recommendations:

1. The authors use "particle number" in order to name the deviation from some
basis state that could be termed vacuum. However, since they use spin models as examples most of the time I suggest to explain this a bit more carefully in the abstract and introduction.

2. In the introduction the authors mention state of the art numerical work in Refs. 30-32. I suggest to add some work of Kristel Michielsen who for a long time
has been holding records in terms of number of spins for dynamical problems, that by the way can be treated in similar ways.

3. A proof of Eq. (7) was presented in [51], not only the result. I leave it to the authors whether they want to present their proof in Sec. A; maybe it is more insightful or more easily accessible.

4. Either above Algorithm 1 or in the discussion it should be mentioned that the scheme developed in [50] applies to mixed spin systems not just Q=2.

5. Discussion on page 14: The use of "half filling" for spin systems is again (see 1) unusual or misleading. At least, the authors should add that this is equivalent to M=0.

6. If the author want to refer to one of the most potent free programs in the field, they could refer to spinpack by Schulenburg that uses the method of Lin. It was e.g. used for finite temperature Lanczos calculations for a kagome system with N=42 spins addressing the full Hilbert space as well as for N up to 72 addressing few magnon spaces in order to describe magnon condensation related to flat bands.
spinpack works in MPI/openMP hybrid mode.

Ask for minor revision

1 - Provides rich context of the lexicographic problem and clearly explains the proposed solution and its implementation.
2 - The solution proposed in the basis enumeration algorithm has immediate and wide application.

1 - An explanation of the titular library is lacking.
2 - The sectioning is sometimes unclear.

"DanceQ: High-performance library for number conserving bases" proposes a scheme for enumerating the computational basis for a generic quantum system, such as spin chains, with possibly larger local Hilbert space. My general impression is that the proposed algorithm represents a valuable contribution to the community and probably deserves to have a software paper associated with it available in a publication, such as SciPost Physics Codebases. However, the submitted article mainly concerns the basis construction algorithm while being presented as being about a library.

The majority of the paper, section III, is dedicated to outlining an algorithm for constructing lookup tables which is well-explained. I appreciate the connection to the most common sectioning of the system into two subsystems as a limit of the algorithm.

I question whether, assuming the paper's topic is the library, the submitted paper meets point 4 of the Scipost Physics Codebases acceptance criteria. A description of the "logic of its workings and its added value" is provided for the algorithm (in great detail in Section III and Appendix B), however, I had to spend several hours trying to understand the library itself from the otherwise rich and detailed documentation available on GitLab. For instance, even an example program for the code is not provided. Section III F exemplifies this problem. There, the authors briefly introduce the Operator class, notably in a tone disconnected from the rest of the paper. This section clearly concerns the library design and not the basis construction algorithm. Further, the State, Basis, Vector, or variants of SparseMatrix and ShellMatrix classes are never mentioned but appear important to the design of the library.

As a consequence of this tension, I see some minor internal inconsistencies. For example: at the beginning of section IV, the authors write that "for optimal performance, the required memory should not exceed [the L3 cache] limit" while in section IV-A they conclude that "the lookup and retrieval of the indices plays only a secondary role during the full matrix-vector multiplication." This observation to me somewhat invalidates the conclusions drawn from Fig. 4.

I ran some example programs provided in the repository and found that the installation process is generally straightforward however the library has important dependencies which need to be configured separately. The instructions are provided however CMake is not configured to attempt to find some libraries automatically. Some further minor observations are included in the following section.

1 - An explanation of the core design philosophy of the package and some example code should probably be added.

2 - I found the discussion about matrix-vector multiplication (MVM) very insightful albeit short, focusing on the memory requirements of the lookup tables. I would be interested in, e.g., more precise timing for a single MVM which the authors characterize as "only slightly dependent on $N$". The "several MPI threads" should also be changed to a number. It would also be interesting to know how MVM with fixed $L$ and $N$ scales with increasing the number of cores. An important example could be comparing the $L=46$ spin system for $N=2$ and $N=3$ regarding memory requirement and time consumption.

3 - In the discussion of Fig. 5, the authors note "a clear trend that larger systems and larger filling fractions require mode subsystems for ideal [lookup table] performance." Yet some outliers are noticed in particular for $Q>1$ for larger $L=32$ in disagreement with this general trend. Can this be explained somehow?

4- Sections III E and IV could be improved by clearer subsectioning. I.e. a single subsection does not make sense to me.

5 - The symbols used in equations should be checked, e.g., the pairs $q-Q$ or $c_0-c_{P_0}$ are sometimes used interchangeably. The caption of Fig. 6 contains several typos.

6 - The authors are aware that Krylow space algorithms are useful beyond (imaginary) time evolution and ground state searches even without storing additional wavefunctions (e.g., Phys. Rev. B 102, 054408), so I am surprised by this omission, e.g., when they write "a Krylov subspace which can used to perform time evolution or to compute ground states and excitations..., or other eigenvectors using spectral transformations."

Ask for minor revision