Updated docs for new memory usage behaviour.

CITATION.cff

@@ -39,5 +39,5 @@ abstract: >-
 keywords:
   - analog quantum computing
   - emulation
-version: 3.0.0
+version: 2.10.0rc1
 date-released: '2026-07-31'

ci/emu_mps/pyproject.toml

@@ -24,7 +24,7 @@ classifiers=[
 ]
 dynamic = ["version"]
 dependencies = [
-  "emu-base==3.0.0"]
+  "emu-base==2.10.0rc1"]
 
 [project.urls]
 Documentation = "https://pasqal-io.github.io/emulators/"

ci/emu_sv/pyproject.toml

@@ -24,7 +24,7 @@ classifiers=[
 ]
 dynamic = ["version"]
 dependencies = [
-  "emu-base==3.0.0"]
+  "emu-base==2.10.0rc1"]
 
 [project.urls]
 Documentation = "https://pasqal-io.github.io/emulators/"

docs/emu_mps/advanced/resource_estimation.md

@@ -44,10 +44,10 @@ $$
 |\mathrm{bath_{CPU}}| < 2s\chi^2\sum_{i=0}^{N/2-1}(3+i) = 4\chi^2N(N+10).
 $$
 
-Just as for the state this is a strict over-estimation that is asymptotically correct, because it assumes all the bonds in the state are of size $\chi$. Note that the baths take up more memory than the state by a factor $O(N)$, where the prefactor depends on the interaction type (for Rydberg it is $1/4$), making this term the most dominant for large qubit numbers. For this reason those baths which are not used in the Lanczos algorithm are always stored on CPU, even when the emulator is configured to use GPU. By overlapping the data tranfers with other computations, we can avoid runtime impact at the cost of occasionally having a third bath tensor present on the GPU during the Lanczos algorithm. In this case we get
+Just as for the state this is a strict over-estimation that is asymptotically correct, because it assumes all the bonds in the state are of size $\chi$. Note that the baths take up more memory than the state by a factor $O(N)$, where the prefactor depends on the interaction type (for Rydberg it is $1/8$), making this term the most dominant for large qubit numbers. For this reason those baths which are not used in the Lanczos algorithm are always stored on CPU, even when the emulator is configured to use GPU. By overlapping the data transfers with other computations, we can avoid runtime impact at the cost of occasionally having a third bath tensor present on the GPU during the Lanczos algorithm. In this case we get
 
 $$
-|\mathrm{bath_GPU}| \leq 2s\chi^2(3h_{max}-1) = 16\chi^2(\frac{3N}{2}+5).
+|\mathrm{bath_GPU}| \leq s\chi^2(3h_{max}-1) = 16\chi^2(\frac{3N}{2}+5).
 $$
 
 ### Contribution from the Krylov space
@@ -62,10 +62,10 @@ where $k$ is the value of `max_krylov_dim`. Recall that the default value of $k=
 
 ### Contribution from temporary tensors
 
-Finally, to compute the above Krylov vectors, the effective two-site Hamiltonian has to be applied to the previous Krylov vector to obtain the next one. The resulting tensor network contraction cannot be done in-place, so it has to store two intermediate results that get very large. The intermediate results take the most memory at the center qubit, where the bond dimension of the Hamiltonian becomes $h$. At this point,
+Finally, to compute the above Krylov vectors, the effective two-site Hamiltonian has to be applied to the previous Krylov vector to obtain the next one. The resulting tensor network contraction cannot be done in-place, so it has to store two intermediate results that get very large. Additionally, the contraction requires storage of a permuted copy of one the two bath tensors. The intermediate results take the most memory at the center qubit, where the bond dimension of the Hamiltonian becomes $h$. At this point,
 
 $$
-|\mathrm{intermediate}| = 2*shp^2\chi^2 = 128h\chi^2
+|\mathrm{intermediate}| = s\chi^2 (2hp^2+h) = h\chi^2(144)
 $$
 
 It should be noted that the value of $h$ cited above assumes that all qubits in the system interact via a two-body term, which is technically true for the Rydberg interaction. However, small terms in the interaction matrix can be truncated away via a configuration option, and this will save memory.
@@ -75,34 +75,37 @@ It should be noted that the value of $h$ cited above assumes that all qubits in
 Putting all of this together, for the total memory consumption $m$ of the program, we can write the following bound for the CPU:
 
 $$
- m_{CPU}(N,\chi,k) = |\psi| + |\mathrm{bath_{CPU}}| + |\mathrm{krylov}| + |\mathrm{intermediate}| < 32N\chi^2 + 4\chi^2N(N+10) + 64*k*\chi^2 + 64(N+4)\chi^2 = 4\chi^2[N(N+34) + 16k + 64].
+ m_{CPU}(N,\chi,k) = |\psi| + |\mathrm{bath_{CPU}}| + |\mathrm{krylov}| + |\mathrm{intermediate}| < 32N\chi^2 + 4\chi^2N(N+10) + 64*k*\chi^2 + 72(N+4)\chi^2 = 4\chi^2[N(N+36) + 16k + 72].
 $$
 
 The quadratic scaling in $N$ is due to the bath contribution. On the GPU, we get a linear scaling instead:
 
 $$
- m_{GPU}(N,\chi,k) = |\psi| + |\mathrm{bath_{GPU}}| + |\mathrm{krylov}| + |\mathrm{intermediate}| < 32N\chi^2 + 16\chi^2(\frac{3N}{2}+5) + 64*k*\chi^2 + 64(N+4)\chi^2 = 8\chi^2[15N + 8k+42]
+ m_{GPU}(N,\chi,k) = |\psi| + |\mathrm{bath_{GPU}}| + |\mathrm{krylov}| + |\mathrm{intermediate}| < 32N\chi^2 + 16\chi^2(\frac{3N}{2}+5) + 64*k*\chi^2 + 72(N+4)\chi^2 = 8\chi^2[16N + 8k+46]
 $$
 
-Note that this estimate is __pessimistic__, since not all $k$ Krylov vectors are likely to be needed, and not all tensors in $\psi$ and the baths have the maximum bond dimension $d$. On the other hand, the estimate for $|intermediate|$ is likely to be accurate, since the bond dimension of $\chi$ is probably attained at the center qubit if the sequence is long enough.
+In this case, there is of course a quadratic scaling in the consumption of CPU RAM. For large system sizes, this will outstrip the memory requirements on the GPU, and it becomes important to ensure enough RAM is available for the program to run.
+
+Note that these estimates are __pessimistic__, since not all $k$ Krylov vectors are likely to be needed, and not all tensors in $\psi$ and the baths have the maximum bond dimension $d$. On the other hand, the estimate for $|intermediate|$ is likely to be accurate, since the bond dimension of $\chi$ is probably attained at the center qubit if the sequence is long enough.
 
 Both TDVP and DMRG rely on the same bath construction and effective Hamiltonian machinery, so their memory requirements are expected to scale similarly with $N$ and $\chi$.
 
-To test the accuracy of the above memory estimations, we benchmarked the __TDVP__ algorithm by fixing the bond dimension to a particular desired value.
-For different combinations of the number of atoms in a register $N$ and the fixed bond dimension $\chi$, we collect the maximum resident size, or RSS, which is expected to capture the maximum memory needed to run the emulation. We plot the RSS in the following picture (left), as a function of the number of qubits and for different bond dimensions. Notice that, once the RSS is normalized by $\chi^2$, as suggested by our estimate above, all the points fall into the same functional dependency on the number of atoms. Moreover, as we plot the normalized function $m(N,\chi,k)/\chi^2$, for a reasonable estimate of the size of the Krylov subspace ($k=30$), it is clear that our upper bound on memory occupation can be reasonably trusted on a wide range of system sizes and bond dimensions.
+To test the accuracy of the above memory estimations, we benchmarked the __TDVP__ algorithm by fixing the bond dimension to a particular desired value. We ran the benchmark on GPU, and consider the memory burden on the GPU.
+
+For different combinations of the number of atoms in a register $N$ and the fixed bond dimension $\chi$, we collect the maximum resident set size, or RSS,  on the GPU, which is expected to capture the maximum memory needed to run the emulation. We plot the RSS in the following picture (left), as a function of the number of qubits and for different bond dimensions. Notice that, once the RSS is normalized by $\chi^2$, as suggested by our estimate above, all the points fall into the same functional dependency on the number of atoms. Moreover, as we plot the normalized function $m(N,\chi,k)/\chi^2$, for a reasonable estimate of the size of the Krylov subspace ($k=10$), it is clear that our upper bound on memory occupation can be reasonably trusted on a wide range of system sizes and bond dimensions.
 
 <img src="../../benchmarks/benchmark_plots/RSS_vs_N.png"  width="49.7%">
 <img src="../../benchmarks/benchmark_plots/emumps_maxRSS_map.png"  width="49.7%">
 
 Finally, having established an estimate for the memory consumption, it makes sense to explore what are the available regimes of qubits/bond dimension can be reached for a given hardware capability.
 Since all heavy simulations will be run on an NVIDIA A100 (on Pasqal's DGX cluster), we have $40$ GB of available memory.
-Therefore, above, we show (right image) the contour lines of the RSS estimate $m(N,\chi,k=30) < 40$ GB for particular useful values of the total memory, allowing to quickly estimate the memory footprint of an emu-mps emulation.
+Therefore, above, we show (right image) the contour lines of the RSS estimate $m(N,\chi,k=10) < 40$ GB for particular useful values of the total memory, allowing to quickly estimate the memory footprint of an emu-mps emulation.
 
 Although these results are shown for __TDVP__, a similar analysis could be performed for __DMRG__, and we expect the resulting memory scaling to be comparable.
 
 ### An example
 
-For example, the results from the [case study](convergence.md) were obtained using $N=49$ and $d=1600$ on 2 GPUs. Taking the above formula, and halving the contributions from $\psi$ and $|\mathrm{bath}|$ since they are split evenly on the GPUs, we reproduce the memory consumption of the program for $k=13$. Notice that the actual number of Krylov vectors required to reach convergence is likely larger than 13, but here we underestimate it since the contributions of $\psi$ and $|\mathrm{bath}|$ are over-estimated.
+For example, the results from the [case study](convergence.md) were obtained using $N=49$ and $d=1600$. Taking the above formula for GPU, we find a memory consumption of $16.998$ GB even for $k=0$, while the actual cited memory consumption on GPU is $14.659$ GB. The actual number of Krylov vectors required to reach convergence is likely around 10, but this contribution is not visible since $\psi$ and $|\mathrm{bath}|$ are over-estimated.
 
 ## Estimating the runtime of a simulation
 

docs/emu_mps/benchmarks/index.md

@@ -19,9 +19,9 @@ These quantities are represented, in the following plots:
 <img src="./benchmark_plots/emumps_maxRSS_map.png"  width="49.7%">
 <img src="./benchmark_plots/emumps_runtime_map.png"  width="49.7%">
 
-The RSS plot (left) shows the peak memory cost of the emulation.
+The RSS plot (left) shows the peak device-dependent memory cost of the emulation.
 It is expected to stay constant at fixed bond dimension and thus represent the total memory occupation of the emulation of a sequence.
-As evident, the emulator is mostly limited by the available memory (40 GB on NVIDIA A100), as it restricts the maximum number of qubits/bond dimension pair allowed.
+As evident, the emulator is mostly limited by the available memory (40 GB on NVIDIA A100), as it restricts the maximum number of qubits/bond dimension pair allowed. Note that in addition to this device-dependent contribution, there is always a contribution in RAM coming from the [bath tensors](../advanced/resource_estimation.md#contribution-from-the-baths). When the simulation is run on CPU, the two contributions add. When running on GPU, the graph can be used to determine what system size can fit on the GPU, and the formula for the bath contribution can be used to determine how much RAM to allocate for the work (assuming access to computing resources where RAM is abundant).
 To get the total estimated runtime instead, one should simply multiply the time estimate in the timing plot (right) by the number of steps in the emulated sequence.
 Finally, given the technical nature of these estimates, they rely on some previous knowledge about [matrix product states](../advanced/mps/index.md) and the [TDVP](../advanced/algorithms.md) algorithm. We encourage anyone who might be interested into the derivation, to have a look at the [resource estimation](../advanced/resource_estimation.md) page in advanced topic section of this documentation.
 

docs/emu_mps/benchmarks/performance.md

@@ -23,15 +23,15 @@ In the right column (both CPU and GPU figure), we explore the available register
 
 ### Performance for a 7x7 grid
 
-Let us now analyze the performance of the simulation for a much larger system. We show results for the adiabatic sequence with 49 qubits arranged in a 7x7 grid. The parameters of the simulation were `max_bond_dim = 1600`, `extra_krylov_tolerance=1e-5` and `precision=1e-6`, and we ran the simulation on 2 GPUs. The maximum bond dimension of the state, its size in memory, the total memory consumption of the program on GPU 1, and the time taken per emulation time step (there are 390 time steps of `dt=10 ns` each) are shown in the graph below.
+Let us now analyze the performance of the simulation for a much larger system. We show results for the adiabatic sequence with 49 qubits arranged in a 7x7 grid. The parameters of the simulation were `max_bond_dim = 1600`, `extra_krylov_tolerance=1e-5` and `precision=1e-6`, and we ran the simulation on GPU. The maximum bond dimension of the state, its size in memory, the total memory consumption of the program on the GPU, and the time taken per emulation time step (there are 450 time steps of `dt=10 ns` each) are shown in the graph below.
 
 <div style="text-align:center;">
 <img src="../benchmark_plots/afm_performance_big.png" width="50%" style="background-color:white;">
 </div>
 
-First off, note that the peak memory consumption on GPU 1 reaches almost 30 GB at the end of the simulation and the memory profile on GPU 2 will be very similar. Note that this memory consumption can be estimated ([see here](../advanced/resource_estimation.md#an-example)), and that the simulation would not have fit on a single GPU. Next, the memory consumption stops increasing as quickly when the maximum bond dimension plateaus, but it does not stop increasing entirely. This is because when the maximum bond dimension reaches the cutoff value of `1600`, most of the tensors in the MPS will not have reached maximum size yet. However, the rate of memory consumption growth will decrease as more of the tensors reach this maximum size.
+First off, note that the peak memory consumption on the GPU reaches almost 15 GB at the end of the simulation. Note that this memory consumption can be estimated ([see here](../advanced/resource_estimation.md#an-example)), and that a significantly larger amount of RAM is used to store bath tensors. Next, the memory consumption stops increasing as quickly when the maximum bond dimension plateaus, but it does not stop increasing entirely. This is because when the maximum bond dimension reaches the cutoff value of `1600`, most of the tensors in the MPS will not have reached maximum size yet. However, the rate of memory consumption growth will decrease as more of the tensors reach this maximum size.
 
-Finally, it can be seen that the time taken per time step scales roughly linearly with the memory consumption of the quantum state. This proposes an obvious method for speeding up the simulation. We've mentioned above, that the qubit density results were insensitive to changes in `max_bond_dim`, and that a value of `max_bond_dim` smaller than `1600` would likely still generate good results. As a consequence of the above graph, we ran the simulation with `max_bond_dim=1000` and all other parameters the same, and indeed, the qubit density was identical up to the 2 digits of precision used above. As can be read off from the above graph, this corresponds to a final time per step of `100s`, yielding a significant reduction in simulation time while keeping the results at the desired precision.
+Finally, it can be seen that the time taken per time step scales roughly linearly with the memory consumption of the quantum state. This proposes an obvious method for speeding up the simulation. We've mentioned above that the qubit density results were insensitive to changes in `max_bond_dim`, and that a value of `max_bond_dim` smaller than `1600` would likely still generate good results. As a consequence of the above graph, we ran the simulation with `max_bond_dim=1000` and all other parameters the same, and indeed, the qubit density was identical up to the 2 digits of precision used above. Furthermore, running the simulation at this reduced bond dimension yielded a final time per step of `55s`, a significant reduction from the final time per step of `195s` from the graph above.
 
 ## Quench
 

docs/emu_mps/index.md

@@ -29,20 +29,32 @@ The following features are currently supported:
     - computational statistics: each time step during the simulation will generate the following information:
         - $\chi$ : is the maximum bond dimension of the MPS
         - $|\Psi|$: MPS (the state) memory footprint
-        - RSS: max memory allocation
+        - RSS: max memory allocation. If run on GPU this is cited both for RAM and GPU.
         - $\triangle t$: time that the step took to run (given in seconds)
 - Specification of:
     - Initial state ( as `MPS` or `MPS._from_state_amplitudes`)
     - Various precision parameters
-    - Whether to run on CPU or GPU(s)
+    - Whether to run on CPU or GPU
     - The interaction coefficients $U_{ij}$ from [here](./advanced/hamiltonian.md#qpu-hamiltonian)
     - A cutoff below which $U_{ij}$ are set to 0 (this makes the computation more memory efficient)
 
 ## Planned features
 
-- More efficient use of GPU by storing tensors on CPU where possible.
 - Differentiability.
 
+## Environment variables
+
+The workflow in emu-mps is not typical of machine learning workloads, and as a consequence certain default behaviours of torch are not optimal. Specifically
+
+- In emu-mps tensor sizes tend to be very unpredictable, in contrast to typical machine learning models. The default torch behaviour of caching GPU memory allocations in deterministic block sizes must be overridden to avoid memory fragmentation that reduces the effectively available memory on the GPU.
+- In emu-mps, we use page-locked RAM for asynchronously transferring tensors from RAM to GPU. Since the tensor sizes are unpredictable, the default torch caching behaviour of allocating memory in powers of 2 works poorly. It is not needed for performance, and it causes the application to use on average 1.5 times more RAM than needed. The host allocator cache should be disabled. This is only possible in torch 2.13 and newer.
+
+Both these things can be configured via an environment variable that should be set before running your python script. Concretely, define the following for optimal memory usage:
+
+`PYTORCH_CUDA_ALLOC_CONF=expandable_segments:True,pinned_max_cached_size_mb:0`
+
+On torch versions lower than 2.13 only use everything before the comma.
+
 ## More Info
 
 Please see the API specification for a list of available config options ([see here](api.md)).