Overview

This section describes the basic working principles of the cuTensorNet library. For a general introduction to quantum circuits, please refer to Introduction to quantum computing.

Introduction to tensor networks

Tensor networks emerge in many mathematical and scientific domains, ranging from quantum circuit simulations, quantum many-body physics and quantum chemistry, to machine learning and probability theory. As network sizes scale up exponentially, there is an ever increasing need for a high-performance library that can drastically facilitate efficient parallel implementations of tensor network algorithms across multiple domains, which cuTensorNet aims to serve.

Tensor and tensor network

Tensors are a generalization of scalars (0D), vectors (1D), and matrices (2D), to an arbitrary number of dimensions. In the cuTensorNet library we follow cuTENSOR’s nomenclature:

  • A rank (or order) \(N\) tensor has \(N\) modes

  • Each tensor mode has an extent (the size of the mode) and can be assigned a label (index) as part of the tensor network

  • Each tensor mode has a stride, reflecting the distance in physical memory between two logically consecutive elements along that mode, in units of elements

For example, a \(3 \times 4\) matrix \(M\) has two modes (i.e., it’s a rank-2 tensor) of extent 3 and 4, respectively. In the C (row-major) memory layout, it has strides (4, 1). We can assign two mode labels, \(i\) and \(j\), to represent it as \(M_{ij}\).

Note

For NumPy/CuPy users, rank/order translates to the array attribute .ndim, the sequence of extents translates to .shape, and the sequence of strides has the same meaning as .strides but in different units (NumPy/CuPy uses bytes instead of counts).

When two or more tensors are contracted to form another tensor, their shared mode labels are summed over. Diagrammatically, tensors and their contractions can be visualized as follows:

../_images/TN_diagrams.png

Here, each vertex represents one tensor object, and each edge stands for one tensor mode (label). When a mode label is contracted, the tensors are connected by the corresponding edges.

A tensor network is a collection of tensors contracted together to form an output tensor. The contractions between the constituent tensors fully determine the network topology. For example, the tensor \(T\) below is given by contracting the tensors \(A\), \(B\), \(C\), and \(D\):

\[T_{abcd} = A_{aij} B_{bjk} C_{klc} D_{lid},\]

where the modes with the same label are implicitly summed over following the Einstein summation convention. In this example, the mode label \(i\) connects the tensors \(D\) and \(A\), the mode label \(j\) connects the tensors \(A\) and \(B\), the mode label \(k\) connects the tensors \(B\) and \(C\), the mode label \(l\) connects the tensors \(C\) and \(D\). The four uncontracted modes with labels \(a\), \(b\), \(c\), and \(d\) refer to free modes (sometimes also referred to as external modes), indicating the resulting tensor \(T\) is of rank 4. Following the diagrammatic convention above, this contraction can be visualized as follows:

../_images/T_ABCD.png

As we can see from the diagram, this contraction has a square topology where each tensor is connected to its adjacent tensors.

Similarly, a quantum circuit can be viewed as a kind of tensor network. As shown in the figure below, single-qubit and two-qubit gates translate to rank-2 and rank-4 tensors, respectively. Meanwhile, the initial single-qubit states \(|0\rangle\) and single-qubit measurement operations can be viewed as vectors (projectors) of size 2. The contraction of the tensor network on the right yields the wavefunction amplitude of the quantum circuit on the left for a particular basis state (ex: \(|010\rangle\)).

../_images/QC_TN.png

Description of tensor networks

A full description of a given tensor network contraction requires two pieces of information: tensor operands and the topology of the tensor network (tensor connectivity hyper-graph). A tensor network in the cuTensorNet library is represented by the cutensornetNetworkDescriptor_t descriptor that effectively encodes the topology graph and data type of the network. To be precise, this descriptor specifies the number of input tensors in numInputs and the number of modes for each tensor in the array numModesIn, along with the modes of each tensor, mode extents, and strides in the arrays of pointers modesIn, extentsIn, and stridesIn, respectively. The tensor network descriptor also accepts an array of qualifiers (qualifiersIn) corresponding to each input tensor, in which each input tensor can be marked as Conjugate (default non-Conjugate), Constant (default non-Constant), and/or requires Gradient (default false). When an input tensor is marked as Conjugate, its data will be complex conjugated upon tensor contraction. When an input tensor is marked as Constant, i.e., its data will not change across subsequent tensor network contractions, cuTensorNet will invoke a caching mechanism (using the provided CACHE workspace memory) to accelerate the computation of subsequent tensor network contractions. See Intermediate tensor(s) reuse. When an input tensor is marked as requires Gradient, its corresponding gradient tensor will be computed upon backward propagation call. Likewise, the tensor network descriptor holds similar information about the output tensor (e.g., numModesOut, modesOut, extentsOut, stridesOut). Note that there is only one output tensor per tensor network.

It is possible for all these network metadata to reside on the host since construction of a tensor network requires only its topology and the data-access pattern; we do not need to know the actual content of input tensors at the tensor network descriptor creation step.

Internally, cuTensorNet utilizes cuTENSOR to create tensor objects and perform pairwise tensor contractions. cuTensorNet’s APIs are designed such that users can focus on creating the tensor network description without having to manage such “low-level” details themselves. The tensor network contraction can be computed in different precisions, specified by the data type given by a cutensornetComputeType_t constant.

Once a valid tensor network is created, one can

  1. Find a cost-optimal tensor network contraction path, possibly with slicing and additional constraints (currently, the cost function can either be the total Flop count or the estimated time to solution)

  2. Access information concerning the tensor network contraction path

  3. Get the needed workspace size to accommodate intermediate tensors

  4. Create a tensor network contraction plan according to the information collected above

  5. Auto-tune the contraction plan to optimize the run time of the tensor network contraction

  6. Perform the actual tensor network contraction to produce the output tensor, either serially or in parallel

It is users’ responsibility to manage device memory for the workspace (from Step 3) and input/output tensors (for Step 5). See API Reference for the cuTensorNet APIs (section Workspace Management API). Alternatively, the user can provide a stream-ordered memory pool to the library to facilitate workspace memory allocations, see Memory Management API for details.

High-level tensor network specification and processing

To simplify the specification and processing of tensor networks encountered in quantum science and other domains, cuTensorNet provides a set of high-level API functions for users to gradually build a given tensor network state and subsequently compute its properties:

  • cutensornetCreateState() is used to create the initial (vacuum) tensor network state in some direct-product tensor space specified by dimensions of all of its constituent vector spaces, such as qudit spaces or, more commonly, qubit spaces (all vector space dimensions are equal to 2 for qubits).

  • cutensornetStateApplyTensorOperator() (cutensornetStateApplyTensor() deprecated) is used to describe an action of tensor operators (e.g., quantum gates) on a tensor network state so as to gradually build the tensor network of interest (e.g., a given quantum circuit) that defines the final tensor network state (a tensor in the defined direct-product tensor space).

  • Once the tensor network state of interest (e.g., the output state of a quantum circuit) has been constructed, its properties can be computed:

    • cutensornetCreateAccessor() specifies a slice of a tensor network state that needs to be computed,

      to inspect the amplitudes of the tensor network state without producing the full state tensor.

    • cutensornetCreateExpectation() specifies the expectation value of a given tensor network operator

      over a given tensor network state. Currently, supported tensor network operators include those defined as a sum of products of tensors acting on disjoint subsets of tensor network state modes. In particular, the Jordan-Wigner transformation produces a subclass of such operators. See cutensornetNetworkOperator_t for more information.

    • cutensornetCreateMarginal() specifies a marginal distribution tensor to be computed for a given tensor network state.

      In physics literature, a marginal distribution tensor is referred to as the Reduced Density Matrix (RDM). It is the tensor obtained by tracing the direct product of the defined tensor network state with its conjugate state (from the dual direct-product tensor space) over the specified tensor modes.

    • cutensornetCreateSampler() creates a tensor network state sampler that can sample the tensor network state

      over specified (either all or a subset of the) tensor modes. When sampling over all tensor state modes, the tensor network state defines the probability distribution in the corresponding direct-product tensor space. Namely, the squared absolute value of each tensor element defines its probability of appearing as a sample in the sampling procedure. In quantum circuit simulations, this would mean generating the output samples from the final state of the quantum circuit according to their probabilities. Sampling over an incomplete set of tensor network state modes would refer to sampling the corresponding marginal probability distribution.

Approximate tensor network algorithms

While tensor network contraction path optimization can greatly reduce the computational cost of the exact contraction of a tensor network, the cost can still quickly scale beyond classical computer’s limit as the network size and complexity increases. A variety of approximate tensor network methods have been developed to address this challenge, such as algorithms based on the matrix product states (MPS). These methods often aim to exploit the sparsity in the tensor network, using tensor decomposition techniques such as QR and SVD.

The QR and SVD decompositions of a tensor can be visualized by the diagram below and can be viewed as the higher order generalization of matrix QR/SVD.

../_images/tensor_decomp.png

The problem can be fully specified by user-requested partitions of modes in the output tensors. For example, the tensor \(T\) in the figure above is decomposed into output tensors \(U\), \(S\), and \(V\) using tensor SVD:

\[T_{iajb} = U_{ijm} S_{mn} V_{ban}\]

Note

Although the output \(S\) is represented as a matrix in the equation and the diagram, only its diagonal elements are nonzero. In all SVD-related APIs of cuTensorNet, only the non-zero diagonal elements are returned (as a vector).

Once the partition of modes in \(U\) and \(V\) is determined, input tensor \(T\) is first permuted to an intermediate tensor \(\bar{T}\) with the proper mode ordering such that it can be viewed as a matrix for the decomposition:

\[T_{iajb} \rightarrow \text{permute} \rightarrow \bar{T}_{ijba} \rightarrow \text{SVD} \rightarrow \bar{U}_{ijm}S_{mn}\bar{V}_{nba}\]

After the matrix SVD/QR decomposition, permutation may still be needed to transform the output intermediate tensors \(\bar{U}\) and \(\bar{V}\), to the requested output \(U\) and \(V\).

Note

Tensor SVD and QR in cuTensorNet operate in reduced mode, i.e., \(k=min(m, n)\) where k denotes the extent of the shared mode of the output, while m and n represent the effective input matrix row/column sizes. While tensor QR is always exact, cuTensorNet provides different ways for users to truncate singular values, \(\bar{U}\), and \(\bar{V}\). For instance, if the user wishes to perform fixed extent truncation in tensor SVD, \(k\) can be set to be lower than \(min(m, n)\). See SVD options for detailed descriptions.

In practice, tensor decomposition and pairwise tensor contraction are often combined to transform the original network to a new topology that reflects its internal sparsity or entanglement geometry in the field of quantum simulation. For instance, in matrix product states (MPS) simulator, the qubits states are mapped to a 1-dimensional tensor chain. When gate tensors are applied to the MPS, a sequence of contraction and tensor QR/SVD operations are performed to ensure that the output states are still represented as a 1-dimensional tensor chain. The compound operation that is used to maintain the MPS geometry is termed as the gate split process. For example, a two-qubit gate tensor \(G\) is applied to two connecting tensors \(A\) and \(B\) to generate output tensors \(\tilde{A}\) and \(\tilde{B}\) with similar modes as \(A\) and \(B\):

\[A_{ipj} B_{jql} G_{pqrs} \rightarrow \text{gate split} \rightarrow \tilde{A}_{irx}\tilde{B}_{jsl}\]

In the gate split process, the number of singular values may also be truncated to keep the MPS size tractable. By iterating this process throughout the network, the original tensor network can be approximated by the output MPS states. The process can be diagrammatically shown as below:

../_images/gate_split.png

For more detailed explanation of how the transformation is performed, please refer to Gate Split Algorithm below.

Contraction Optimizer

A contraction path is a sequence of pairwise contractions represented in the numpy.einsum_path() format. For a given tensor network, the contraction cost can differ by orders of magnitude, depending on the quality of the contraction path. Therefore, it is crucial to use a path optimizer to find a contraction path that minimizes the total cost of contracting the tensor network. Currently, the contraction cost refers to the floating point operations (FLOP) for the cuTensorNet path optimizer’s purpose, and an optimal contraction path corresponds to the path with the minimal FLOP count.

The cuTensorNet path optimizer is based on a graph-partitioning approach (called phase 1), followed by interleaved slicing and reconfiguration optimization (called phase 2). Practically, experience indicates that finding an optimal contraction path can be sensitive to the choice of configuration parameters, which can be configured via cutensornetContractionOptimizerConfigSetAttribute() and queried via cutensornetContractionOptimizerConfigGetAttribute(). The rest of this section aims to introduce the overall optimization scheme and these configurable parameters.

Once the optimization completes, it returns an opaque object cutensornetContractionOptimizerInfo_t that contains all of the attributes for creating a contraction plan. For example, the optimal path (CUTENSORNET_CONTRACTION_OPTIMIZER_INFO_PATH) and the corresponding FLOP count (CUTENSORNET_CONTRACTION_OPTIMIZER_INFO_FLOP_COUNT) can be queried via cutensornetContractionOptimizerInfoGetAttribute().

Note

The returned FLOP count assumes the input tensors are real-valued; for complex-valued tensors, multiplying it by 4 is a good estimation.

Alternatively, we can bypass the cuTensorNet optimizer entirely by supplying your own path using cutensornetContractionOptimizerInfoSetAttribute(). For the unset cutensornetContractionOptimizerInfo_t attributes (ex: the number of slices and the FLOP) cuTensorNet will either use the default values or compute them on the fly when applicable.

Graph partitioning

In general, searching for the optimal contraction path for a given tensor network is a NP-hard problem. Therefore, it is crucial to approach this problem in a divide-and-conquer spirit. Given an input tensor network, we first perform graph partitioning to split the network into N smaller subgraphs and repeat this process recursively until the size of each subgraph is less than the cutoff size. Once the final partitioning scheme is determined, finding the paths for contracting first within each subgraph and then all subgraphs becomes tractable. The process can be illustrated using the diagram below and we get a coarse contraction path or contraction tree for subsequent optimizations.

../_images/graph_partition.png

The cuTensorNet library provides controls over the graph partition algorithms through following parameters:

Slicing

In order to fit a tensor network contraction into available device memory, as specified by workspaceSizeConstraint, it may be necessary to use slicing (also known as variable projection or bond cutting). By slicing, we split the contraction of the entire tensor network into a number of independent smaller contractions where each contraction considers only one particular position of a certain mode (or combination of modes). The total number of slices that are created equals the product of the extents of the sliced modes. The result for full tensor network contraction can be obtained by summing over the output of each sliced contraction. Taking the above \(T\) tensor as an example, if we slice over the mode i we obtain the following representation:

\[T_{abcd} = A_{aij} B_{bjk} C_{klc} D_{lid} \longrightarrow \sum_{i_s} \left( A_{a {i_s} j} B_{bjk} C_{klc} D_{l {i_s} d} \right),\]
../_images/slicing.png

As one can see from the diagram, the sliced mode \(i_s\) is no longer implicitly summed over as part of tensor contraction, but instead explicitly summed in the last step. As a result, slicing can effectively reduce the memory footprint for contracting large tensor networks, in particular quantum circuits. Besides, since each sliced contraction is independent from the others, the computation can be efficiently parallelized in various distributed settings.

Despite all the benefits above, the downside of slicing is that it often increases the total FLOP count of the entire contraction. The overhead of slicing heavily depends on the contraction path and the modes that are sliced. In general, there is no simple way to determine the best set of modes to slice.

The slice-finding algorithm in the cuTensorNet library can be tuned by following parameters:

The slice-finding algorithm in the cuTensorNet library utilizes a set of heuristics and the process can be coupled with the subtree reconfiguration routine to find the optimal set of modes that both satisfies the memory constraint and incurs the least overhead.

Reconfiguration

At the end of each slice-finding iteration, the quality of the contraction tree has been diminished by the slicing. We can improve the contraction tree at this stage by performing reconfiguration. Reconfiguration considers a number of small subtrees within the full contraction tree and attempts to improve their quality. The repeated process of slicing and reconfiguration can be viewed as illustrated in the diagram below.

../_images/subtree_reconfiguration.png

Although this process is computationally expensive, a sliced contraction tree without reconfiguration may be orders of magnitude more expensive to execute than expected. The cuTensorNet library offers some controls to influence the reconfiguration algorithm:

  • CUTENSORNET_CONTRACTION_OPTIMIZER_CONFIG_RECONFIG_NUM_ITERATIONS: Specifies the number of subtrees to consider during each reconfiguration. The amount of time spent in reconfiguration, which usually dominates the optimizer run time, is linearly proportional to this value. Based on our experiments, values between 500 and 1000 provide very good results. Default is 500. Setting this to 0 will disable reconfiguration.

  • CUTENSORNET_CONTRACTION_OPTIMIZER_CONFIG_RECONFIG_NUM_LEAVES: Specifies the maximum number of leaf nodes in each subtree considered by reconfiguration. Since the time spent is exponential in this quantity for optimal subtree reconfiguration, selecting large values will invoke faster non-optimal algorithms. Nonetheless, the time spent by reconfiguration increases very rapidly as this quantity is increased. Default is 8. Must be at least 2. While using the default value usually produces the best FLOP count, setting it to 6 will speed up the optimizer execution without significant increase in the FLOP count for many problems.

Deferred rank simplification

Since the time taken by the path-finding algorithm increases quickly as the number of tensors increases, it is advantageous to minimize the number of tensors, if possible. Rank simplification removes trivial tensor contractions from the network in order to improve performance. These contractions are those where a tensor is only connected with at most two neighbors, effectively making the contraction a small matrix-vector or matrix-matrix multiplication. One example for rank simplification is given in the figure below:

../_images/rank_simplification.png

This technique is particularly useful for tensor networks with low connectivity; for instance, all single-qubit gates in quantum circuits can be fused with the neighboring multi-qubit gates, thereby reducing the search space for the optimizer. The necessary contractions to perform the simplification are not immediately performed, but, rather, are prepended to the contraction path returned. If, for some reason, such simplification is not desired, it can be disabled:

While simplification helps lower the FLOP count in most cases, it may sometimes (depending on the network topology and other factors) lead to a path with a higher FLOP count. We recommend that users experiment with the impact of turning simplification off (using the option listed above) on the computed path.

Hyper-optimizer

cuTensorNet provides a hyper-optimizer for the path optimization that can automatically generate many instances of contraction paths and return the best of them in terms of the total FLOP count. The number of instances is user-controlled via the use of CUTENSORNET_CONTRACTION_OPTIMIZER_CONFIG_HYPER_NUM_SAMPLES which is set to 0 by default. The idea here is that the hyper-optimizer will create CUTENSORNET_CONTRACTION_OPTIMIZER_CONFIG_HYPER_NUM_SAMPLES instances, each reflecting the use of different parameters within the optimizer algorithm. Each instance will run the full optimizer algorithm including reconfiguration and slicing (if requested). At the end of the hyper-optimizer loop, the best path (in term of FLOP counts) is returned.

The hyper-optimizer runs its instances in parallel. The desired number of threads can be set using CUTENSORNET_CONTRACTION_OPTIMIZER_CONFIG_HYPER_NUM_THREADS and is chosen to be half of the available logical cores by default. The number of threads is limited to the number of the available logical cores to avoid the resource contention that is likely with a larger number of threads.

Internally, cuTensorNet holds a thread pool for multithreading in the hyper-optimizer etc. The thread pool lifetime is bound to the cuTensorNet handle. (In previous releases, OpenMP was used; this is no longer the case since cuTensorNet v2.0.0.)

The configuration parameters that are varied by the hyper-optimizer can be found in the graph partitioning section. Some of these parameters may be fixed to a given value (via cutensornetContractionOptimizerConfigSetAttribute()). When a parameter is fixed, the hyper-optimizer will not randomize it. The randomness can be fixed by setting the seed via the attribute CUTENSORNET_CONTRACTION_OPTIMIZER_CONFIG_SEED.

Intermediate tensor(s) reuse

cuTensorNet currently supports caching/reuse of intermediate tensors for cases where some of the input tensors vary in value, but the network structure does not. In such cases, it can drastically accelerate repeated executions of the same tensor network where some tensors update their values. One example is calculation of the probability amplitudes of individual bit-strings or small batches of them, the procedure often used for the verification/validation of quantum processors. In this case, only the calculation of the first bit-string amplitude incurs the full computational cost whereas the subsequent calculations of bit-string probability amplitudes should run much faster, benefiting from the intermediate tensor caching/reuse. Another example is the sampling of the probability distribution of the final quantum circuit state. The sampling procedure will also benefit from intermediate reuse since the underlying reduced density matrices computed for groups of qubits do not change their structure while only updating the values of a relatively small number of input tensors. Finally, variational quantum algorithms with a relatively small number of variational parameters may also be accelerated with this new feature.

To activate intermediate tensor reuse for constant input tensors, users need to:

Approximation setting

SVD Options

In addition to the standard tensor SVD mentioned above, cuTensorNet allows users to perform truncation, normalization, and re-partition during the decomposition using various SVD algorithms. For truncation, the most common strategy is to keep the largest k number of singular values and trim out the remaining ones. In cuTensorNet, this can be easily specified by modifying the extent of the shared mode in the output cutensornetTensorDescriptor_t object. Alternatively, users may truncate the singular values based on the actual value or the distribution of the eigenvalues during runtime. Such truncation setting, SVD algorithm along with normalization and re-partition options are provided by different attributes of the cutensornetTensorSVDConfig_t object:

Note

The value-based truncation options in cutensornetTensorSVDConfig_t can be used in conjunction with the fixed extent truncation, but since the actual reduced extent found at runtime is unknown, the size of the data allocation for the output tensors should always be based on the assumption of no value-based truncation. After the execution, the potentially reduced extent found at runtime will be stored in the CUTENSORNET_TENSOR_SVD_INFO_REDUCED_EXTENT attribute of cutensornetTensorSVDInfo_t and also reflected in the output cutensornetTensorDescriptor_t. Users can either call cutensornetTensorSVDInfoGetAttribute() or cutensornetGetTensorDetails() to query this information.

Gate-split algorithm

In general there are more than one approach to perform the gate split operation. cuTensorNet provides two algorithms in cutensornetGateSplitAlgo_t for this purpose:

  • CUTENSORNET_GATE_SPLIT_ALGO_DIRECT: The three input tensors will first be contracted to form an intermediate tensor. A subsequent tensor SVD will take place to decompose the intermediate tensor to the desired output. This algorithm is generally more memory demanding but more performant when the tensor sizes are relatively small.

  • CUTENSORNET_GATE_SPLIT_ALGO_REDUCED: The input tensors A and B will first be decomposed into smaller fragments via tensor QR. The two R tensors from the decomposition will then be contracted with input tensor G to form an intermediate tensor. A subsequent tensor SVD will take place to decompose the intermediate tensor, the output of which will then contract with the Q tensors from the first step to form the desired outputs. This algorithm is generally less memory demanding and more performant when the tensors sizes get large.

Note

For CUTENSORNET_GATE_SPLIT_ALGO_REDUCED, cuTensorNet may internally skip QR on tensor A or/and B if no memory size reduction can be achieved in the decomposition.

The two algorithms can be diagrammatically visualized as follows (singular values are partitioned equally onto the two outputs in this example):

../_images/gate_split_algo.png

Supported data types

For tensor network contraction, a valid combination of the data and compute types inherits in a straightforward fashion from that of cuTENSOR. Please refer to cutensornetCreateNetworkDescriptor() and cuTENSOR’s User Guide for details.

For approximate tensor network functionalities including tensor SVD/QR decompositions and gate split operation, following data types are supported:

Note

For gate split operation, the valid compute types for the data types above are consistent with that of tensor network contraction.

References

For a technical introduction to cuTensorNet, please refer to the NVIDIA blog:

For further information about general tensor networks, please refer to the following:

For the application of tensor networks to quantum circuit simulations, please see:

Citing cuQuantum

    1. Bayraktar et al., “cuQuantum SDK: A High-Performance Library for Accelerating Quantum Science,” 2023 IEEE International Conference on Quantum Computing and Engineering (QCE), Bellevue, WA, USA, 2023, pp. 1050-1061, doi: 10.1109/QCE57702.2023.00119.