Xped is c++ library for the definition and manipulation of symmetry preserving tensors for tensor network related algorithms.
Theory
Tensor
Multiple definitions exist for the term tensor (see e.g. wikipedia). For the purpose of this library, a tensor \(T\) is a multidimensional linear map between the tensor product of \(d\) vector spaces \(\mathbb{C}^{n_1}\otimes\cdots \otimes\mathbb{C}^{n_d}\) which form the domain \(\mathcal{D}\) and a tensor product of \(c\) vector spaces \(\mathcal{\mathbb{C}^{m_1}\otimes\cdots \otimes\mathbb{C}^{m_c}}\) which form the codomain \(\mathcal{C}\):
\begin{align*} T:\mathbb{C}^{n_1}\otimes\cdots \otimes\mathbb{C}^{n_d} &\to \mathbb{C}^{m_1}\otimes\cdots \otimes\mathbb{C}^{m_c} \\ v_1\otimes\cdots\otimes v_d &\mapsto w_1\otimes\cdots\otimes w_c \end{align*}
The rank of the tensor \(T\) is \(r=d+c\). The dimension of the domain is \(dim_{\mathcal{D}}=n_1\cdots n_d\) and of the codomain \(dim_{\mathcal{C}}=m_1\cdots m_c\). Since the map is multi-linear, its action is entirely defined by the action on the basis vectors. This defines the components of the tensor with respect to given bases:
\[ j_1\otimes\cdots\otimes j_c = \sum_{i_1,\dots,i_d}T_{i_1,\dots,i_d}^{j_1,\dots,j_c} i_1\otimes\cdots\otimes i_d \]
By a canonical isomorphism \(\pi_\mathcal{D}\), one can define the composite vector space \(\mathbb{C}^{dim_\mathcal{D}}\) for the domain of the map:
\[ \mathbb{C}^{dim_\mathcal{D}}\cong\mathbb{C}^{n_1}\otimes\cdots \otimes\mathbb{C}^{n_d} \]
Analogously, one can define a composite space \(\mathbb{C}^{dim_\mathcal{C}}\) for the codomain. The tensor T is then an ordinary linear map:
\[ T: \mathbb{C}^{dim_\mathcal{D}} \to \mathbb{C}^{dim_\mathcal{C}} \]
This map can be represented by a matrix which will be useful when storing the tensor in memory.
Symmetric tensor
A symmetric tensor is a tensor which respects a given symmetry. To formalize this, lets assume a symmetry described by a group \(G\). \(G\) is represented on all individual vector spaces of the domain by \(\mathcal{R}^i_G\) and codomain by \(\mathcal{Q}^i_G\) of a tensor. Then a tensor \(T\) is symmetric if:
\[ T(v_1\otimes\cdots\otimes v_d) = \mathcal{Q}^1_G\otimes\cdots\otimes \mathcal{Q}^c_G\left( T(\mathcal{R}^1_G(v_1)\otimes\cdots\otimes \mathcal{R}^d_G(v_d))\right) \]
That means the action of the tensor stays the same when the symmetry acts on the individual spaces. If the tensor is transformed with the isomorphisms \(\pi_\mathcal{D}\) and \(\pi_\mathcal{C}\) it becomes an ordinary linear map. Then Schur's lemma (wikipedia) implies that this map is a constant if the composite representation of the symmetry on domain and codomain is irreducible. Generally, it will be reducible but may be decomposed into a direct sum of irreducible representations (irreps). The underlying matrix becomes then block diagonal.
Construction of isomorphism
The block diagonal form of the underlying matrix is only manifest if the isomorphism \(\pi_\mathcal{D}\) maps to a basis which is the decomposed into the direct sum of irreps. Hence, the crucial task is to find an appropriate isomorphism. A straight-forward method for obtaining \(\pi_\mathcal{D}\) is the use of fusion trees which perform the mapping by a sequential pairwise fusion of two vector spaces into their product state. For a given binary fusion of two spaces, the Clebsch-Gordan expansion is used to find the decomposition of the product state into irreps. Each intermediate space is then combined with next unpaired space until all spaces are combined to the final domain/codomain of the tensor. It is important to store all the intermediate quantum numbers since for general symmetries, several different paths exist which lead to the same output. An example fusion tree of four quantum numbers \(q_1,\dots,q_4\) coupled to \(q_f\) looks as follows:
The triangular shapes correspond to the binary Clebsch-Gordan decomposition. They carry another (Greek) label which is important for symmetries with outermultiplicity in their Clebsch-Gordan expansion (e.g. SU(3)). For Abelian symmetries, the complete fusion tree is entirely determined by the incoming quantum numbers and the intermediate values do not have to be stored.
The fusion trees construct the isomorphisms \(\pi_\mathcal{D}\) and \(\pi_\mathcal{C}\).
- Note
- The isomorphisms \(\pi_\mathcal{D}\) and \(\pi_\mathcal{C}\) are not constructed explicitly but only implicitly by storing all fusion trees.
Tensor operations
For tensor network related algorithms, several tensor operations are ubiquitous. One needs to contract tensors together to obtain final results or needs to decompose a tensor with a singular value decomposition to truncate certain bonds in the network. A fundamental operation is the permutation of the legs of an individual tensor. It is needed for other high-level operations.
Permutations
A permutation is a reshuffling of the individual legs of the tensor. This corresponds to a reordering of the individual vector spaces which form \(\mathcal{D}\) and \(\mathcal{C}\). For nonsymmetric tensors, this can be achieved by choosing appropriate new isomorphisms \(\pi_\mathcal{D}^\prime\) and \(\pi_\mathcal{C}^\prime\) such that the permuted input is mapped to the correct state. This leads to a new view to the original data which can be left untouched. However, the data layout becomes strided and noncontiguous. For some operations, it is therefore more efficient to copy the underlying data so that the layout is again contiguous.
In the case of symmetric tensors, the situation is more complicated. A permutation requires a reordering of the incoming legs of the fusion trees. If the incoming legs are shuffled, the resultant fusion tree is not canonical anymore and therefore not useful for further operations which expects a canonical fusion tree. In order to get the reordered fusion tree in canonical form, one needs the recoupling symbols \(F\) of the symmetry group \(G\).
The sum goes over all possible intermediate quantum numbers \(q_{23}\) and in the case of outermultiplicity also over \(\gamma\) and \(\delta\). The recoupling symbols define a unitary transformation between the composite spaces of different fusion orders. They are entirely determined by the Clebsch Gordan coefficients.
To obtain the reordered fusion tree in canonical form, one needs a second basis operation which exchanges two incoming quantum numbers on the same node. This operations involves the swap symbol for the group.
With these two operations, one can implement arbitrary permutations on a single fusion tree. A tensor consists of two fusion trees, one for the domain and one for the codomain. For permutations which shuffle across the domain and codomain, a third basic operation is needed. This involves the bending of a line pointing inwards to pointing outwards and is quantified by the turn symbol of the group.
The permutation of the legs of a symmetric tensor is then performed by first permuting the fusion tree pair (consisting of domain tree and codomain tree) and second by assigning the old data to the new tensor with careful attention to the recoupling-, swap- and turn symbols.
Contractions
A tensor contraction boils down to the composition of the associated maps. I.e., if \(T_1: \mathcal{D}\to\mathcal{I}\) and \(T_2: \mathcal{I}\to\mathcal{C}\) than the contraction \(T\) of \(T_1\) and \(T_2\) over all the elementary spaces in \(\mathcal{I}\) is the composition:
\[ T = T_2 \circ T_1 \]
The composition of tensors corresponds to the multiplication of the representing matrices which arise after the application of the isomorphisms \(\pi_\mathcal{D}\) and \(\pi_\mathcal{C}\).
In the generic case, the indices of a tensor which should be contracted with indices of another tensor are located in its codomain. Therefore one has to permute the indices of both tensors so that an ordinary composition performs the contraction. After the composition, another permute might bring the indices in the desired order for the final result.
Decompositions
A tensor decomposition can be seen as the inverse operation to a contraction. Because of that it can be handled similarly. An often encountered decomposition is the singular value decomposition (svd). It can be applied to any linear map. Hence, the tensor is again interpreted as an ordinary linear map with the help of the isomorphisms \(\pi_\mathcal{D}\) and \(\pi_\mathcal{C}\). Then, the svd reads:
\[ T = U \Sigma V^\dagger \]
\(U\) has the original domain of the tensor and \(V\) has the original codomain, i.e. the cut is taken through the domain and codomain.
For a general partition of indices, one can apply a permute operation to get an appropriate domain and codomain.
Implementation
Quick tour
To get a feeling of how the library works, lets look at some simple examples for the core part of the library. All entities of the library are defined in namespace Xped.
Symmetries
Each supported symmetry has its own class which supplies all methods related to the symmetry. Currently implemented are the following groups:
- \(Z_n\) \(\rightarrow\)
template<std::size_t N, Xped::Sym::Kind> struct Xped::Sym::Zn; - \(\text{U}(1)\) \(\rightarrow\)
template<Xped::Sym::Kind> struct Xped::Sym::U1; - \(\text{SU}(2)\) \(\rightarrow\)
template<Xped::Sym::Kind> struct Xped::Sym::SU2; - trivial \(\rightarrow\)
struct Xped::Sym::U0;
Partial support exists also for general \(\text{SU}(N)\) symmetry.
All symmetry classes have a CRTP base class template<typename Derived> struct Xped::Sym::SymBase; to avoid duplicated code. The individual symmetry classes are stateless and methods are static. Some basic objects which are implemented in the symmetry classes are:
using qType = ...the type for the quantum numbers of the symmetrystatic std::vector<qType> basis_combine(qType ql, qType qr)returns the decomposition of the direct product of two irreps into irreps.static constexpr qType qvacuum();returns the quantum number of the trivial representation.
To see the full list of required functions check Symmetry/SU2.hpp
All symmetry classes take Xped::Sym::Kind as a template parameter. It describes to which particles the symmetry relates to. Possible values are:
The difference between the first two is only when formatting related quantum numbers because spin quantum numbers can be half-integers but particle quantum numbers are only integers. The latter however, makes a crucial difference. If a symmetry is set up with Xped::Kind::Fermion, the fermionic commutation rules and respective signs are taken care off automatically for tensor operations. This is controlled by the methods
The first methods adds a minus sign if two odd fermion-parity states are swapped. The second includes a sign if a odd parity state is twisted.
Fusion trees
The struct FusionTree is a simple static storage of all relevant data for a single fusion path of Rank quantum numbers of a given Symmetry:
It has the following members:
- Warning
- This class does not check whether all the quantum numbers do fit together but it expects that it is constructed correctly.
To display a FusionTree, there is the handy function Xped::FusionTree::draw(). It generates an ascii drawing of the tree. Here is an example for a rank=4 tree for an SU(2) symmetry:
Fusion tree manipulations
Fusion tree manipulations build the basis for tensor operations. The fundamental operation is the permutation of two adjacent incoming quantum numbers. Unless it involves the first two incoming quantum numbers which enters to the same binary fusion, this involves recoupling operations. Here, especially the \(6j\)-symbol of the symmetry group enters. At the end, a fusion tree can be written as a weighted sum over fusion trees with two adjacent quantum numbers exchanged. The respective method is Xped::FusionTree::swap:
With the exchange of adjacent quantum numbers, one can implement arbitrary permutations of the incoming quantum numbers by decomposing the permutation into a chain of adjacent swaps. This is performed with Xped::util::Permutation::decompose(). The high-level function which performs the general permutation is Xped::FusionTree::permute().
Hilbert bases
To describe a Hilbert basis, the class Qbasis is used:
This class stores a basis which is properly sorted with respect to a given Symmetry, i.e. it contains several multiplets of the symmetry. Furthermore, the basis can be a composite basis of several individual bases. This is controlled via the template parameter depth. To handle the composition, Qbasis uses FusionTree to remember the exact fusion process. Generally, bases are constructed with depth=1 and composite bases are obtained using the method Xped::Qbasis::combine().
The most important member of Qbasis is:
The actual type is more complicated because of AllocationPolicy but this does not matter here. Each entry of data_ corresponds to a quantum number of type qType (first entry in the tuple). For each quantum number, the global number of the first plain state corresponds to the second element in the tuple. The third element of the tuple belongs to the plain basis states of this quantum number. The plain basis state are described by the simple class Xped::Basis.
A few lines of code illustrate the use of the class:
One of the key methods of Qbasis are Xped::Qbasis::leftOffset() and Xped::Qbasis::rightOffset(). If several bases are combined into a composite basis, each basis state in the composite basis belongs to a specific fusion. A rigorous order of all the fused states with same quantum number is fundamental to construct the isomorphism \(\pi_\mathcal{D}\). The exact signature of leftOffset() is:
It takes a FusionTree as a parameter and returns the number of basis states which come before this tree. These methods are fundamental for getting tensor views from the large block diagonal matrix which is stored in memory. See the implementation of Xped::Tensor::subBlock():
Tensor class
The tensor class is the central object of the library:
It takes 6 template parameters but the last two are optional. The first four describe the underlying Scalar type, the rank and corank of the tensor ( \(d\) and \(c\)) and the symmetry under which the tensor is invariant. The fifth parameter decides whether the tensor can be used for automatic differentiation (AD). It is a bool parameter and both choices (true and false) corresponds to different explicit specializations of the general class template. For the AD part see below. The last parameter controls the allocation procedure.
The tensor class has only a single member variable:
The used storage is currently determined at configuration time by a cmake parameter. In the future, this can become an additional template parameter.
The most important methods are:
Xped::Tensor::permute();Xped::Tensor::contract();Xped::Tensor::tSVD();
Storage
For an efficient storage of the tensor data in the case of a symmetry constrain, the tensor is seen as an ordinary linear map between two vector spaces. This procedure was described in the Theory section above. In the case of a symmetry, the representing matrix of the linear map is block-diagonal. Therefore, the most compact storage is obtained when storing only the nonzero block-matrices of the large block-diagonal matrix. Each block is associated to a quantum number of the symmetry and contains all elements of the symmetry which can be fused to this quantum number.
Two specific storage implementations are available:
VecOfMatsee Core/storage/StorageType_vecofmat.hppContiguoussee Core/storage/StorageType_contiguous.hpp
The first option stores a std::vector of matrices to describe the block-diagonal matrix. For this option the entire tensor is not contiguous in memory.
The second option stores a large contiguous buffer (std::vector<Scalar>). When a specific block of the block-diagonal matrix is requested, a view into this large buffer is returned which corresponds to the requested block.
The Storage class must have two specific constructors.
Both take the individual bases of the domain and codomain as an std::array and the world to which the storage is associated in the case of distributed parallelism. The first constructor does not allocate actual memory but only sets the meta data of the storage. The second constructor allocates the needed storage and copies the contiguous data data into the tensor. If size does not match the total elements in the tensor, an assertion is raised.
Additionally, the storage implementations need to include the following public methods.
const MatrixType& block(std::size_t i) const;returns a const reference to the ith block of the block-diagonal matrixMatrixType& block(std::size_t i);returns a nonconst reference to the ith block of the block-diagonal matrixconst MatrixType& block(qType q) const;returns a const reference to the block belonging to the fused quantum numberq.MatrixType& block(qType q);returns a nonconst reference to the block belonging to the fused quantum numberq.std::size_t plainSize() const;returns the total number of scalars in the tensor summed over all blocks.const mpi::XpedWorld& world() constreturns the mpi world, the data is living in. This is only relevant for distributed tensors.
For the complete interface check out Core/storage/StorageType_contiguous.hpp.
Automatic differentiation
The Tensor class supports build in automatic differentiation (AD), i.e. for algorithms build on top of the Tensor class, the derivative can be computed automatically. To declare a tensor with AD support, one needs to set the respective template argument to true. The source code for the AD support for Tensor can be found in include/Xped/AD/ADTensor.hpp (and not in Tensor.hpp itself). It is an explicit specialization and of the Tensor class template. In order to obtain the AD functionality, core features of the third-party library stan/math are used.
Let's look at an example usage of the AD functionality:
Backends
The Xped library is designed to be a high-level library for symmetric tensors. All operations for plain nonsymmetric tensors are reimplemented. Instead, one can choose between different backends that deliver the functionality for plain tensors (or matrices). It is also possible to add a backend by implementing the required interface.
Currently, three backends are available. The Eigen backend is tested most thoroughly.
Eigen
Eigen is a famous c++ library for linear algebra. It is designed with an elegant API and highly optimized matrix operations. Additionally, it has a tensor module which is the basis for the TensorFlow framework by google. This part supports arbitrary ranked tensors and a wealth of operations including contractions and reductions.
The Eigen library also supports the dispatch of fundamental matrix operations to external BLAS and LAPACK implementations. This feature can be used within Xped at configuration time using the XPED_USE_BLAS and similar options. If the intel math kernel library should be integrated, one can switch on the parameter XPED_USE_MKL.
Cyclops tensor framework (ctf)
ctf is a tensor framework for distributed parallelism. This allows the computations to spread over several different compute cores and allows to use the hardware more efficient. It has an overhead for small tensors but allows to scale the compute power with very large tensors.
- Note
- The ctf backend is most useful for tensors without symmetry or with tensors with small internal symmetries.
array
- Todo:
- Add documentation for array backend.
- Warning
- The array backend needs some adaptions to be consistent with the code base.
Algorithms
On top of the core library, Xped also provides high-level tensor network algorithms.
Matrix product states (MPS)
MPSs represent a one-dimensional tensor network tailored for strongly interacting one-dimensional physical lattice systems. For detailed information, check the following references:
- https://link.springer.com/article/10.1007/BF02099178
- https://arxiv.org/abs/1008.3477
- https://arxiv.org/abs/2011.12127
The MPS code is located in the directory Xped/MPS and currently provides a class MPS and basic operations for them. An implementation of matrix product operators (MPO)s and algorithms like the density matrix renormalization-group (DMRG) or the time-dependent variational principle (TDVP) are future development goals.
A small example demonstrates the current capability:
Projected entangled pair states (PEPS)
PEPS is a generalization of MPS for two-dimensional systems. For detailed information, check the following references:
The PEPS code of Xped is located in Xped/PEPS and has the following capabilities:
- CTMRG algorithm for arbitrary unit cells with custom pattern (AB/BA)
- Imaginary time evolution using the simple update method with nearest and next-nearest Hamiltonian terms
- Direct energy-minimization using automatic differentiation
For PEPS simulations, the configuration is parsed as a toml file. The toml file has sections for CTM, imaginary time evolution and nonlinear optimization. The used model is specified in the section [model] and needs to be present in all simulations: Here is an example for the Hubbard model:
The unit cell and PEPS bond dimension are specified in the section [ipeps]:
CTM
The CTM procedure can be configured using toml configuration section [ctm]
Simple update
Imaginary time evolution can be configured through the toml section [Imag].
Energy minimization
Simulations for direct energy minimization using nonlinear optimization procedures are configured through the toml section [optim].
