{'fraction_non_alphanumeric': 0.030280857770023078, 'fraction_numerical': 0.008611470704114369, 'mean_word_length': 4.74256625727214, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Simulation of Quantum Many-Body Systems on Amazon Cloud

Quantum many-body systems (QMBs) are some of the most challenging physical systems to simulate numerically. Methods involving tensor networks (TNs) have proven to be viable alternatives to algorithms such as quantum Monte Carlo or simulated annealing, but have been applicable only for systems of either small size or simple geometry due to the NP-hardness of TN contraction. In this paper, we present a heuristic improvement of TN contraction that reduces the computing time, the amount of memory, and the communication time. We demonstrate our heuristic with the Ising model on powerful memory optimized Amazon Web Services (AWS) x1.32x large EC2 instances, showing the viability of cloud computing for scientific applications.


Introduction and Motivation
Quantum many-body (QMB) physics is concerned with the study of microscopic systems involving a large number of interacting particles [1]. Studies of QMB have applicability across a wide range of physics, chemistry, and material science problems [2,3,4], such as the study of magnetism, superconductivity, topological order, and spin glasses [5,6].
Several tensor network algorithms have been formulated to implement computer simulation of these systems [5,7,8,9,10,11,12], but so far such simulations have been limited to systems of modest size and particular geometries due to the computational hardness of the tensor contraction [13,14]. Approximate methods to overcome these limitations have been proposed, for instance, by simplifying the tensor environment thus, avoiding a full contraction [15,16,17,18]. However, often the focus of the literature has been on infinite, translation-invariant systems and no particular attention has been paid to adapt the methodology to distributed computing. In this paper, we investigate a QMB spin system and explore the limitations of a cloud computing environment using on a heuristic for parallel TN contractions without approximations. For this study we use the Amazon Web Services (AWS) Elastic Compute Cloud (EC2) instances.
Tensors are multi-indexed data structures that can be viewed as multi-dimensional arrays. The order of a tensor is defined by the number of indices needed to specify a tensor element. For instance, scalars are tensor of zero order, while vectors and matrices are tensors of order one and two, respectively. A tensor network (TN) is a decomposition of a high-order tensor into a set of low-order tensors which share indices under a specified geometry. To extract information from a TN, it is necessary to perform a summation over all shared indices, a procedure termed tensor contraction.
TNs are now ubiquitous in QMB simulations because they provide a systematic way to represent and approximate quantum wave functions. These wave functions can have exponentially many components. For instance, the Hilbert space of a QMB system of N interacting spin-1/2 particles (each being a two-state subsystem) has dimension 2 N on the spin sector. A tensor representation for such a system requires a tensor of order N where each index represents a spin in the lattice. A more compact and often more efficient representation is obtained by decomposing this high-order tensor into a TN, where shared tensor indices are representative of the entanglement between particles generated by their interactions.
Even though the Hilbert space of a QMB system grows exponentially with the number of degrees of freedom, it is also often the case that only a small region of that space contains useful information about the state of the system [19]. This is particularly the case for systems containing local interactions and whose ground-state energy has a gap separating it from excited states; this gap does not scale with system size. Thus, to find the ground state of such gapped systems involving local interactions, only a small region of the Hilbert space needs to be included in the system's ground state wave function. This phenomenon is captured by how the amount of entanglement between one subsystem and another scales with their sizes. If the entanglement entropy scales only with the size of the boundary separating the two subsystems, the so-called "area law" [20], an efficient TN representation of the ground state is possible, with tensors being polynomially bounded. In some cases, particularly near critical points of QMB systems, the entanglement scales with the volume of the smallest subsystem. In these cases, a TN is known to require an exponential number of resources for an accurate representation of the wave function.
To advance the study of QMB systems, even away from critical points, it is necessary to optimize TN contractions. The optimization of TN contractions is a currently extensive area of study [21,22,23,24,4,25]. Of particular importance, it is noted that the order of summations taken is critical to the determination of the computational time necessary for a contraction. Orus demonstrates the importance of this by comparing two different orderings for the contraction of three tensors, each with indices taking D values [26]. Finding an optimal order is NP-hard problem [27]. Equally important, as we will show in this paper, is the partitioning of the system into various independent parallel contractions.
The development of new processors such as the Tensor Processing Units (TPU) and the creation of libraries such as the Cyclops Tensor Framework (CTF) [28] provides minimal help in the advancement of QMB simulations. TPUs are optimized for matrix-matrix operations, but do not accommodate the exponential growth of tensors throughout the contraction. The partitioning of tensor across multiple processors in CTF is too general, failing to take advantage of the geometry of a physical system. What is needed for efficient contraction without approximation is a set of processors with a large cache and physical memory and a partitioning of the TN across processors in a systematic and geometrically advantageous way. We focus our attention on implementing parallel geometric TN partitioning on AWS EC2 instances with large memories.
This paper is organized as follows: In Sec. 2, we discuss progresses in cloud services that have motivated their use for the simulation of QMB systems. In Sec. 3, we describe the QMB system model used for the numerical calculations. In Sec. 4, we examine the computational difficulties involved in the contraction of a TN representing a many-body spin lattice system and present a heuristic for parallel contraction. In Sec. 5, we compare various algorithms for the contraction of TNs of the spin system, and in Sec. 6 we present the performance analysis of our algorithm on the specific EC2 x1e instance type, along with some benchmark results for the spin system. The conclusions are given in Sec. 7.

Cloud QMB Simulation
The largest QMB simulation problems solved in recent years have typically been relegated to supercomputers. For instance, a group from ETH Zurich was able to perform a 45-qubit simulation using a supercomputer at the Lawrence Berkeley National Laboratory [29]. However, cloud computation has made significant advances recently, making it a viable cost-effective alternative for QMB simulations.
The cloud computing infrastructure is designed to perform optimally for Big Data, online transaction processing, and data streaming applications. Such applications exploit data-level, task-level, and thread-level parallelism. The Warehouse scale computers (WSCs), which serve as the backbone of the cloud infrastructure, host tens to hundreds of thousands processors communicating through networks with sufficient bandwidth but also with a relatively high latency. This infrastructure is advantageous for the enterprise cloud applications, but the communication latency significantly affects the performance of applications in science and engineering. Proce-dures such as Remote Procedure Calls (RPCs), serialization, deserialization, and compression of buffers account for 22-27% of the CPU cycles. This is an unavoidable "WSC architecture tax" [30,31] QMB applications typically exhibit fine grained parallelism, deploying many parallel threads communicating frequently with each other, and using barrier synchronization to transit from one stage of computation to the next. Because of the communication latency of current cloud services, scientific applications are commonly performed on supercomputers with fast interconnection networks such as Mirinet, Infiniband, or some custom designed network.
In recent years, Cloud Service Providers (CSPs) have narrowed the performance gap vis-a-vis supercomputers. Clusters with faster interconnects are now offered. Instances with physical memory on the order of hundreds of GiB are being provided. Faster processors and coprocessors are being developed in new instance types. For example, Amazon provides Graphics Processing Units (GPUs), optimal for linear algebra operations, while Google has pioneered the Domain Specific Architectures (DSAs) with the creation and provision of Tensor Processing Units (TPUs), which are optimal for the many small matrix-matrix operations needed in deep learning applications. However, the challenge of managing high order tensor contractions in QMB applications is still very much present. Both the number of operations and the memory footprint involved in the computation grow exponentially with the system size.
Our approach for the implementation of QMB simulation on AWS instances includes selecting the instances providing the largest storage, allowing contractions to be carried out concurrently, and reducing the amount of communication between parallel threads. We avoid distributing a single tensor across multiple threads and choosing instead to distribute groups of tensors according to the geometry of the system considered.

The QMB model system
For the numerical computations, we adopt the spin-1/2 Ising model in the presense of a transverse field as our QMB system. The Ising model is used in statistical mechanics to describe magnetic lattice systems with strong anisotropy [36]. Recently, it has been used as a paradigm for the study of quantum phase transitions [37]. The model consists of discrete variables S z i that represent magnetic dipole moments of atomic spins that can be in one of two states (+1) or (−1). The spins are arranged in a lattice, allowing each spin to interact with its nearest neighbors (spin-spin interaction). For spatial dimensions larger than one, the model has a finite-temperature phase transition and critical behavior.
When a transverse magnetic field is present, the model yields a zero-temperature phase transition (i.e., a quantum phase transition) driven by the competition between the spin-spin interaction, which favors ferromagnetism (if J > 0) or antiferromagnetism (if J < 0), and the external field, which favors paramagnetism. In this paper, we consider the case where the spins are located on a rectangular lattice. Mathematically, the model is defined by the Hamiltonian (total energy) of the system, with the two termŝ describing, respectively, the spin-spin interactions between nearest-neighbor sites and the coupling of the spins to a transverse field (here denoted by Γ). The constant J quantifies the spin-spin interaction. Different components of the on-site spin operator do not commute, namely,Ŝ z iŜ x i =Ŝ x iŜ z i , lending the two terms in the Hamiltonian non commuting. It is this noncommutability that yields quantum critical behavior and entangled many-body states for the spin system. At zero temperature, one finds a critical point when Γ ≈ 3J [38]. At this point, the spins in the system are highly entangled.
The ground state energy E 0 corresponds to the lowest eigenvalue of the operationĤ. All other eigenvalues are associated to excited states.
The vector describing the wave function of the system can be written as where σ k = ±1, with k = 1, . . . , N indicating the N spin degrees of freedom. The connection between these variables and conventional binary ones is straightforward: x k = (σ k + 1)/2. Notice that there are 2 N basis vectors |σ 1 · · · σ N . The amplitudes A({σ k }) are in general complex numbers; however, for the model in consideration they can always be defined as real. The eigenvector |Ψ 0 associate to E 0 yields the ground state wave function of the system, namely,Ĥ|Ψ 0 = E 0 |Ψ 0 . (In the absence of the external field, the ground state is two-fold degenerate due to the spin inversion symmetry ofĤ J .) An important quantity associate to the state vector is its norm, There are also a number of physical quantities of importance that can be obtained from the state vector. The expectation value of the total energy of the system is defined as The local transverse and longitudinal magnetizations of the system at site k are given by and respectively. Finally, the longitudinal spin-spin correlation between spins at sites i and j is equal to By representing the amplitudes A({σ k }) as a tensor network, all the physical quantities above can be computed via suitable tensor network contractions.

Tensor Network Contraction
As previously mentioned, a TN is a decomposition of a high-order tensor into a set of low-order tensors sharing internal indices under a specific geometry. Consider, for instance, a tensor with index set {x i } k=1,...,N , with x k = 0, 1 for all k = 1, . . . , N , whose elements are expressed as [A] x 1 ,x 2 ,...,x N . This tensor has 2 N elements thus require an exponential amount of storaga memtory. A possible decomposition of this tensor is given by the expression where {α k } k=1,...,N , are the internal (repeated) indices of the network, with each index α k = 1, . . . , χ i for a suitable χ k (often referred as bond dimension). This particular chain decomposition is known as a matrix product state (MPS) [35], see Fig. 1a. An MPS is obtained by repeated applications of singular-value decomposition operations on the original tensor A. Notice that each index x i now resides on an individual tensor A i of order two (for k = 1, N ) or three (k = 2, . . . , N − 1). If χ k ∼ poly(N ), Equation (10) provides a compact decomposition of tensor A, requiring only a polynomial amount of storage space.
A matrix product state is not the only possible decomposition of a tensor. Consider, for instance, the case when N = N h N w . One can then decompose the tensor A into a N w × N h rectangular lattice, see Fig. 1b. As mentioned in Sec. 1, for representative classes of QMB systems, it is indeed the case that bond dimensions are polynomially bounded, and the tensor network decomposition of a groundstate wave function provides a compact representation when performed appropriately. The most suitable decomposition minimizes the bond dimensions χ k and is determined by the interactions present and the system geometry. For instance, the decomposition in Fig. 1b is particularly useful for the representation of the quantum amplitude A({σ k }) of the Ising model wave function, see Eq. (4), as each tensor in the lattice can be associated to one physical spin.
To contract a tensor decomposed into a network, it is necessary to perform a summation over all the internal (repeated) indices in the network. Tensor contraction can be done in different ways depending upon a number of factors including the network topology.
By virtue of our choice of QMB system, in this paper we consider only planar, rectangular TNs, as shown in the example of Fig. 1b. In addition, we focus on the computation of scalar quantities such as those defined in Eqs. (5) to (9), which can be cast as the contraction of two planar TNs (one for A * and another for A) into a single planar TN with no external indices. Thus, the computation of phys- ical quantities requires only the sum over all internal indices of a planar TN. More specifically, consider a square lattice of tensors with size L × L. The tensors on the four corners of the lattice are of order two, those along the edge are of order three, and those within the bulk are of order four. In practice, the full contraction of a TN occurs by contracting tensors pairwise and sequentially. The number of tensor elements (we call dimension), an important quantity for determining memory requirements, evolve as follows. Consider the contraction of two tensors, A 1 and A 2 , with dimensions d(A 1 ) and d(A 2 ), respectively, into a tensor B. The dimension d(B) of the resulting tensor satisfies where x is the set of shared index between A 1 and A 2 . Thus, the tensor dimension can increase substantially after a pairwise contraction. For every full TN contraction there exists a bottleneck contraction, after which every tensor pair contraction no longer increases the memory footprint. It is our desire to minimize the size of this bottleneck, which in turn optimizes the memory requirements and the number of floating point operations (FLOPs) necessary to perform the computation. Naively, we could consider simply contracting the square lattice along the rows from one edge to the other. Assuming the dimension of each initial index was χ, this would leave a chain of tensors each hav-ing χ L elements, which could then be contracted as an MPS. However, this is not the only option. If instead we define four quadrants for the lattice and contract from each edge to the midpoint of the lattice, we end with a ring of tensors each having χ L/2 elements, which can then be fully contracted as an MPS. Another option yet is to contract from the edge somewhere between these two previous approaches, ending up with at least one tensor in the final ring with more than χ L/2 elements. Thus, as mentioned earlier, the order of the contraction matters.
In addition to selecting an appropriate contraction order, we partition the lattice for parallel computation according to both the geometry of the lattice and contraction order. Figure 2 demonstrates this for both of the previously discussed contraction orders. Notice that the partitioning of the quadrant scheme allows for the fewest number of messages between processes.
Our choice of contraction ordering and geometryspecific partitioning can be extended to multiple geometries. The heuristic is as follows: given a TN geometry, select the order which minimizes the bottleneck contraction size and partition the lattice around the ring of tensors for parallel computation with minimal communication. In Sec. 6, we compare these parallel contractions to the contraction of square lattice where individual tensors are cyclically parallelized, in contrast to our parallelization of groups of tensors within the lattice.

QMB Ground State Computation
In this section we discuss a few of the fundamental concepts behind the determination of the ground state of a QMB systems using TNs. A particularly important family of TN states are projected entangled pair states (PEPS) [11]. In contrast to the one-dimensional nature of MPS, PEPS correspond to higher dimensional tensor networks, e.g., two-dimensional TNs [7,15,17,18,33,34]. For gapped systems with local interactions, PEPS is particulaly suitable since: (a) bond dimensions of tensors in the TN -the range χ of which the indexes run -can be kept small and yet provide an accurate description of low-energy states; PEPS naturally satisfy the area-law scaling of the entanglement entropy; and (d) PEPS can handle reasonably well two-point correlation functions that decay polynomially with the separation distance.
Following the scheme shown in Fig. 1b, for our simulations we associate to each spin in the bulk of the two-dimensional lattice a tensor of rank 5. Four indices of this tensor account for the bonds to nearest-neighbor spins. The fifth index accounts for the binary nature of the spin variable. For spins located at the edges and corners of the lattice, only three and two bonds are required, respectively. When all bond (internal) indices are contracted, the resulting scalar quantity yields the probability amplitude A({σ k }) of finding the spin system in the particular basis state {σ k }, k = 1, ..., N,, see Eq. (4).
A standard algorithm for the determination of the ground state energy of a PEPS is the imaginary time evolution (ITE) algorithm. This algorithm is defined by the iterative application of an incremental imaginary time evolution operator U δτ = e −iĤ (−iδτ ) = e −Ĥ δτ to an initial quantum state |Ψ initial , over m steps. This process evolves the system's state in imaginary time in incremental steps δτ = τ /m, where τ is the total time of the evolution. For τ max{J −1 , Γ −1 }, the state vector |Ψ = m i=1Û δτ |Ψ initial becomes exponentially close to the ground state |Ψ 0 , provided that the initial state had a nonzero overlap with the ground state, namely, Ψ initial |Ψ 0 = 0. Typically, one chooses as the initial state |Ψ initial a product state that is a random superposition of individual spin states. We initialize our TN to a uniform superposition of product states, (i.e., the uniform bond dimension is initially χ = 1). Each iteration of the time evolution can be divided into subintervals defined by the application of noncommuting operators in the HamiltonianĤ. For our example, the HamiltonianĤ in Eq. (1) is a summation of the noncommuting operatorsĤ J andĤ Γ . If we consider the second-order Trotter-Suzuki approximation [39] ofÛ δτ , then we obtain whereÛ Γ δτ /2 = e −Ĥ Γ (δτ /2) andÛ J δτ = e −Ĥ J δτ . These operators naturally define three subintervals according to the operators involved within each time step δτ .
To accurately monitor the convergence of the algorithm to the optimal state, the expectation value of the total energy of the system is periodically evaluated every two δτ steps, see Eq. (6). For the Ising model, this results in the calculation of the local magnetization m z i , see Eq. (8), and the and twospin correlation c ij , see Eq. (9), for each spin site and site pair, respectively. This too defines two additional sub-steps during the algorithm.
After applying the sets of operators {Û Γ δτ /2 } to every site or {Û J δτ } to every pair, the system is normalized, |Ψ → |Ψ / √ N , maintaining both the stability of the algorithm and the probabilistic interpretation of |Ψ .
Within a single time step, the actual order of events is as follows: a) update the system by applying the one-body operators to every site, and then normalize. b) update the system by applying the two-body operators to every pair of sites, and then normalize. c) after every two iterations of applying operators to the system, calculate the expectation value of the total energy of the system. With all of this in mind, it is now understood that each stage of the imaginary time evolution algorithm allows for the modification or evaluation of the system in three distinct ways: 1. a tensor can be locally updated by the application an operator; 2. the norm of the system can be calculated, or 3. an expectation value can be calculated.
We now examine how each of these three tasks are accomplished. To update the system and push it towards the ground state, the operatorsÛ Γ δτ /2 andÛ J δτ are applied iteratively as previously mentioned. These operators are themselves composed of a series of local operators applied sequentially to either each single site or to pair of sites in the lattice. For the onebody operatorsÛ δτ /2 acting on site i, a tensorization of the operator is constructed and then contracted along the physical index of site i. Full contraction of the lattice is unnecessary at this stage.
For the two-body operatorsÛ acting on a pair of sites ij, the update includes an additional step. First, the tensorization in the tensor product basis of the pair's physical indices is constructed. Then the contraction over the physical indices is performed. The resultant tensor is then decomposed back into the two individual sites in their respective basis by performing a singular value decomposition, as shown in Fig 3. The singular value decomposition is a necessary step in allowing the exchange of information throughout the lattice, which, in turn enables an increase in the entanglement entropy. This increase in the entanglement entropy leads to the growth of bond dimensions. If the system is far away from the critical point, i.e., the gap between the ground state and the first excited state is finite and system-size independent, the imaginary-time evolution encounters only a low entanglement growth (i.e. small bond dimensions are maintained as the system size grows).
However, as the critical point is approached, there is a significant increase in the entanglement. Entanglement scales up significantly with the system size in this case, significantly increasing the bond dimension necessary for accurate simulation. The interplay between bond dimension and lattice size is of fundamental importance to the determination of what types of problems are practically solvable. As mentioned earlier, the quantum critical point for the twodimensional Ising model with transverse field occurs where the quantity Γ/J ≈ 3. Figure 4 shows that for small lattices far from this point, the final state obtained while allowing only minimal bond growth is in good agreement with the solutions obtained from exact diagonalization. We therefore primarily focus on systems where Γ/J is far from the critical point. After applying the nonunitary operatorsÛ Γ δτ /2 and U J δτ to the two-dimensional tensor lattice, the state must be renormalized to keep the algorithm numerically stable. The norm is calculated simply by contracting the tensor network state |Ψ with its conjugate Ψ| along every spin index. The network must be similarly fully contracted for every site and pair in the lattice to calculate E, m x i , and c ij . For each lattice site i the stateŜ x i |Ψ is contracted with Ψ|. For each pair ij , the stateŜ z iŜ z j |Ψ is contracted with Ψ|. Figure 5 shows that the most computationally demanding portion of the algorithm is the calculation of expectation values. For a L×L lattice, the number of contractions are L 2 and 2L(L − 1) for the calculation of m x i and c ij , respectively. Optimization of this portion of the algorithm is the chief motivation behind the previously discussed geometry specific lattice partitioning and quadrant contraction ordering.

Experimental Results
The code implementing the ITE algorithm discussed in the previous section was run on the AWS Elastic Computer Cluster (EC2). We used the X1.32x large EC2 instance. This instance runs four Intel Xeon E7 8880v3 processors, offering 1952 GiB of DRAM and up to 128 virtual CPUs (vCPUs). All contractions between pairs of tensors were performed by folding the tensors into matrices and utilizing the optimized Basic Linear Algebra Subprogram (BLAS) Doubleprecision General Matrix Multiplication (DGEMM) routine [40].
As a first test of the algorithm performance, we compare each contraction scheme previously mentioned in order to verify whether an analysis of the bottleneck tensors is beneficial in determining the optimal algorithm for tensor contraction. We then determine the limits of the algorithm and the types of lattices it can handle. We also compare this to a lattice contracted using the Cyclops Tensor Framework library (CTF) [28]. This library focuses on the parallelization of individual tensors by distributing tensor elements cyclically across processors for generic lattice geometries, as opposed to our algorithm which is optimized for a specific lattice geom- etry. Figure 6 shows a comparison of the execution time of the quadrant contraction scheme, the row contraction scheme, and the CTF contraction for a square lattice size L × L, every element set to value 1, and uniform bond dimension χ . Quadrant contraction is more efficient than the row contraction, as anticipated by the analysis of each algorithm's bottleneck tensors. For the lattice sizes considered, the quadrant contraction is also demonstrated to be favorable over the CTF, proving the importance of latticespecific optimizations in the cloud environment. Focusing therefore on the quadrant contraction, we also determined the average memory and communication costs for various system sizes. Figure 7 shows that the memory footprint and the communication costs grow exponentially with the system size.
To analyze the interplay between the bond dimension and the lattice size, we allowed for larger bond dimensions representing systems in regimes where entanglement is expected to be large (e.g., near critical points). Figure 8 shows the execution time for systems of various sizes, up to a uniform bond dimension χ = 4, again contracting lattices with each element set to value 1. The maximum lattice sizes for χ = 2 and χ = 4 are L = 12 and L = 6, respectively, based on the L3 cache of the AWS instanced used.
We performed the imaginary time evolution simulation of the Ising model with a transverse field strength of Γ = 1 and coupling strength J = 1 for varying lattice sizes L, measuring the time to convergence. Time steps were fixed to δτ = 3/100. The singular value cutoff parameter was fixed to = 0.01, where λ k /λ 1 ≥ , to temper the growth of the bonds in the system. As anticipated, the time for completion of the ITE scales exponentially with the system size L, as shown in Fig. 9.
Moving away from systems of low entanglement, we also performed simulations of the ITE at a transverse field strength of Γ = 3 and coupling strength J = 1 (i.e., near the critical point). The lattice lengths tested were L = 6 and L = 8. For L = 6 and L = 8, each time step δ was set to δτ = 3/75 and δτ = 4/250, respectively. The ground states were reached with a computational runtime of 12.5 hours and 293.3 hours. As anticipated, bond dimensions grew larger in this regime, reaching its maximal  Figs. 10 and 11. The observables we calculate are not in exact agreement, which we hypothesize to be indicative of the necessity of larger bond dimensions χ. They do, however, qualitatively exhibit the behavior indicative of a phase transition near Γ = 3. At the paramagnetic phase (Γ > 3J), the on-site transversed magnetization σ x is maximum, while in the ferromagnetic phase (Γ < 3J) the on-site longitudinal spin-spin local correlator σ z σ z is maximum instead.

Conclusions
At this time supercomputers built around lowlatency interconnection networks are still the best and often the only option for running scientific and engineering codes. Tensor network contraction is a CPU-intensive and memory-intensive application with a large memory footprint exhibiting fine- grained parallelism; this explains why in the past TN-contraction was mostly done on supercomputers.
Distributing tensors over the multiple nodes of a cloud instance increases the communication intensity, therefore slows down the execution on a cloud with relatively high communication latency. Using TPUs will further amplify the gap between computing and communication speed and will shorten the execution time of small size networks; whenever tensors are distributed across multiple nodes only a network with low latency will show significant performance advantage.
This motivated our approach of using instances with a large memory rather than partition and distribute tensors among multiple nodes. Cloud supplies significantly cheaper computing resources than supercomputers and our results show that this approach can be used for a range of problems of interest. Several EC2 instances types offer cost effective alternatives to simulating quantum many-body systems with tensor networks, such as the x1.32x large instance used in our analysis.
Our tensor network contraction implementation minimizes the communication costs but is limited by the memory size of the vCPU. This limitation is determined by the largest single matrix that can be stored in the L3 cache hosting several vCPUs, a trade off we make to minimize communication. For the particular example of the two-dimensional Ising model in the presence of a transverse field, this limitation in memory size appears for lattice sizes of 6×6 (i.e., for a 2 36 dimensional Hilbert space), causing some deviation from the expected results when the quantum entanglement in the simulated system is high (i.e. at Γ J = 3). Yet, the results are encouraging given the affordability and accessibility of EC2 instances.
The tensor contraction procedures discussed in Sections 4 and 6 perform better that the Cyclops library for square lattices with L ≤ 10. Imminent advancements in the cloud infrastructure will benefit applications exhibiting fine-grained parallelism including QMB simulations. Increased memory and larger L1 and L2 caches of individual cores, as well as larger last-level cache of multicore processors will increase the range of problems that can be solved without the benefit of tensor partitioning, as described in this paper. It is unclear if faster networks will make tensor partitioning more appealing as the communication complexity of very large problems is likely to increase faster than the benefits due to lower communication latency.
Profiling the code for imaginary time evolution shows that the runtime of the algorithm is dominated by two in house procedures, contract tVtl par and contract lattice mpi. The first, responsible for the contraction of a tensor pair, dominates  10, these results are also qualitatively in agreement with those obtained with tree tensor networks [41]. The expectation value saturates to a minimal value near Γ = 3, which is indicative of a phase transition.
the execution time and is called 31,347 times for the 6x6 ITE run, accounting for 97.18 % of the execution time. contract lattice mpi, a subroutine responsible for the parallelized contraction of a lattice, is the parent function of contract tVtl par and it is called a total of 3,843 times. Within this subroutine itself, the time per call is dominated by the calls made to the contract tVtl par procedure, taking 1.07 ms of the total 1.1 ms per call. This validates the assumption that the algorithm bottleneck is the contraction of the two largest tensors. Further confirmation of this is found in the low communication cost, evidenced by each call to contract lattice mpi requiring 4 calls to MPI Send and 2 calls to MPI Recv, and yet only accounting for a combined 2.7 % of the time per call for contract lattice mpi. This analysis shows that the search for optimal contraction algorithms for tensor networks is critical for solving increasingly large problems [42]. An optimal tensor network contraction algorithm could reduce significantly the communication complexity and the memory footprint required for the contraction of the largest tensors, the bottleneck of the process discussed in Sections 4 and 6.
In addition, it shows that by properly partitioning a given lattice (so as to minimize inter-process communication), and selecting appropriate contraction orderings, cloud services can be effectively utilized for QMB simulation applications. This is particularly true in the regime where bond dimensions remain manageable (i.e. away from critical points). Further work will investigate hybrid algorithms involving both individual tensor parallelization (as in CTF) and the geometry specific parallelization approach (as in our algorithm); we will also explore different lattice geometries under the heuristic of the minimal bottleneck tensor contraction.