Abstract
ITensor is a system for programming tensor network calculations with an interface modeled on tensor diagram notation, which allows users to focus on the connectivity of a tensor network without manually bookkeeping tensor indices. The ITensor interface rules out common programming errors and enables rapid prototyping of tensor network algorithms. After discussing the philosophy behind the ITensor approach, we show examples of each part of the interface including Index objects, the ITensor product operator, tensor factorizations, tensor storage types, algorithms for matrix product state (MPS) and matrix product operator (MPO) tensor networks, quantum number conserving blocksparse tensors, and the NDTensors library. We also review publications that have used ITensor for quantum manybody physics and for other areas where tensor networks are increasingly applied. To conclude we discuss promising features and optimizations to be added in the future.
Contents
 1 Introduction
 2 Interface Examples
 3 Index Objects
 4 The ITensor Product Operator
 5 Tensor Decompositions
 6 Tensor Storage Layer
 7 High Level Features: MPS and MPO Algorithms
 8 Quantum Number BlockSparse ITensors
 9 NDTensors Library
 10 Other Features of ITensor
 11 Applications of ITensor
 12 Future Directions
 A Full Code Examples
 B ITensor Implementation and Interface in the C++ Language
1 Introduction
Tensor networks are a technique for working with tensors which have many indices [Orus:2014a, Cichocki:2016, Schollwoeck:2011, Ostlund:1995, Hackbusch:2009, Oseledets:2011]. The cost of naively storing and computing with a tensor having indices (an order tensor) scales exponentially with . A tensor network is a representation of a large, highorder tensor as the contracted product of many loworder tensors. When all of the tensors in the network are loworder, a tensor network can make it efficient to perform operations such as summing two highorder tensors or computing their inner product. These operations remain efficient whether the highorder tensor represented by the network has hundreds, thousands, or even an infinite number of indices.
Describing tensor networks can be difficult when using traditional notation: one must come up with distinct names for indices and it can be hard to see the connectivity pattern of the network. An elegant alternative is tensor diagram notation [Penrose]. In diagram notation, tensors are shapes and indices are depicted as lines emanating from them. Connecting two index lines means they are contracted or summed over. For example, the following diagram is equivalent to the traditional expression below it:
Diagram notation is enormously helpful for expressing tensor networks, as it emphasizes key aspects of tensor algorithms while suppressing implementation details such as the ordering of tensor indices. It is just as rigorous as traditional index notation.
ITensor, short for intelligent tensor, is a software library inspired by tensor diagram notation. Its goal is enabling users to translate a tensor diagram into code without reintroducing concepts not expressed by tensor diagrams. For example, when summing two ITensors the only requirement is that they have the same set of indices in any order; the ITensor system handles all other details of performing the sum.
Two “philosophical” principles guided the design of ITensor. The first was that in using the library, any implementation details which are not a part of the conceptual algorithms should be kept hidden from the user as much as possible. Not having to think about these details allows one to focus more clearly on the essentials. A key early insight was that this principle could be tied to the ordering of indices in an ITensor. In typical tensor software, the user is constantly thinking about the order of the indices. However, tensor diagrams do not have any ordering of the indices, just labels that keep track of the relevant information. Thus ITensors have the ordering of their indices abstracted away as a hidden implementation detail, using intelligent indices. The second key principle is that the software should allow one to interact with it at a variety of levels. At a high level, for calculations done in a standard way one can call functions encapsulating a sophisticated algorithm (say, the density matrix renormalization group, DMRG) without understanding much of the implementation details. At an intermediate level, one can gain flexibility by working with moderately sophisticated routines, such as adding MPS. And finally, to do something more novel, one can work at the lower level of ITensors. This multilevel access mandated that ITensor be a library, not a single executable program with complicated input files.
These initial principles led to other important design choices over time. One consequence of having intelligent tensor indices as distinct data objects (of type Index
) is that they can store extra information about themselves. A key use case is indices which have internal subspaces labeled by conserved quantum numbers. Storing this information in the Index
objects ensures that when contracting two quantumnumber conserving ITensors, both are guaranteed to use the same ordering of the subspaces for storing their data. Another consequence is that both dense and various kinds of sparse ITensors can all be of the same type and have essentially the same interface, because the implementation can inspect the indices and internal storage type to determine the actual ‘type’ of any ITensor. Thus users can write very generic code that works for any type of ITensor.
ITensor was first implemented in C++ and extensively developed and refined through three major releases over ten years.^{1}^{1}1ITensor Github Repository (C++): https://github.com/ITensor/ITensor Recently, ITensor has been ported to the Julia language.^{2}^{2}2ITensor Github Repository (Julia): https://github.com/ITensor/ITensors.jl In what follows we show examples in Julia, though we emphasize that the highlevel C++ and Julia interfaces are quite similar (see the Appendix for full code examples in each language). Both versions are full implementations of ITensor in each language: the Julia version is not a wrapper around the C++ version.
The goal of this article is to provide a highlevel overview of the ITensor system, its design goals, and its main features. Much more information including detailed documentation of the ITensor interface, code examples, and tutorials can be found on the ITensor website: https://itensor.org.
2 Interface Examples
We first introduce ITensor by giving examples as an informal overview. In later sections, we will discuss many more details of the individual elements making up the ITensor system such as “intelligent” tensor indices, tensor factorizations, and blocksparse ITensors.
2.1 Installing ITensor
Julia features a builtin package manager that makes installing libraries simple. To install the ITensor library, all a user has to do is issue the following commands, starting from their terminal:
The first line starts an interactive Julia session and the second enters package
manager mode. The command add ITensors
downloads and installs all the
dependencies of the ITensors.jl package
then finally the ITensor library itself.
2.2 Basic ITensor Usage
To begin using the ITensors package in a Julia session or script, input the line
Before creating an ITensor, one first creates its indices. The line of code
creates a tensor Index of dimension 3 and assigns this Index object to the reference
i
.
Upon creation, this Index is stamped with an immutable, unique id number which allows copies of
the Index to be compared and matched to one another.
A portion of this id is shown when printing the Index, with typical example output
of the command @show i
being:
After making a few Index objects i,j,k,l
one can define ITensors:
Because matching Index pairs can automatically recognize each other through their id numbers, tensor contraction can be carried out as:
The *
operator finds all matching indices between two ITensors and
sums over or contracts these indices. The i
Index is summed in the first
contraction above and j
in the second, leaving D
with indices l
and k
.
The ITensor product operator “*
” can also be used for outer products and scalar products,
and is discussed in more detail in Section 4.
2.3 Setting ITensor Elements
Setting an element of an ITensor A = ITensor(i,j,k)
is done by
which assigns the value 0.837 to the element of A
for which index i
has the value 2, j
has value 3, and k
has value 1.
(Note that in Julia, the notation x=>y
makes a Pair(x,y)
object.) ITensor indices are 1indexed, similar to Julia arrays.
Because the Index objects are provided along with their values, they can be passed in any order. Thus the following lines of code
have exactly the same effect on the ITensor A
.
To create an ITensor with normallydistributed random elements instead of specific values, one can use the constructor
to make a realvalued random tensor or
to construct a complexvalued random ITensor.
2.4 Matrix Example
To illustrate the usefulness of the ITensor approach involving Index objects
and the *
operator,
consider a pair of order2 tensors (matrices)
In a typical matrix or tensor library, to contract A
with B
and sum
over their shared index j
, one would need to write code similar to
Note that the above line is not ITensor code!
Within ITensor, all one needs to do is to write
and the *
operator handles the transposition of B
automatically.
If B
is redefined with the ordering of its indices reversed, the operation
A * B
continues to give the correct result. This type of behavior makes ITensor
applications robust to changes in the code that may modify the ordering of tensor indices
or the layout of tensors in memory.
2.5 Summing ITensors
ITensors can be added and subtracted as long as they have the same set of indices. Even if the indices are in a different order, addition always works straightforwardly because the ITensor system is able to internally deduce the data permutation required:
ITensors may also be subtracted and multiplied by scalars, including complex scalars, for example:
2.6 Priming Indices
Sometimes it is not desirable to contract all of the indices shared between two tensors. Consider two ITensors
and say we want to contract only over the index j
leaving the i
indices uncontracted.
A convenient way to achieve this while still using the *
operator is to prime one of the i
indices
The ITensor Ap
has the same elements as A
but has indices (i’,j)
. When contracting Ap
with B
, now only the j
indices will match or compare equal, so it will
be the only Index contracted
Diagrammatically we can notate the above contraction as:
2.7 Compiling ITensor
Although the experience of using Julia is similar to using an interpreted language, it is actually a just in time compiled language. As of version 1.5 of the Julia language, the initial compilation overhead when starting a Julia program can be rather significant. To reduce this overhead, we offer a convenient way for users to compile most of the ITensors.jl code ahead of time with the following commands within an interactive Julia session:
The compilation process can take many minutes, but only has to be done once for each
version of ITensors.jl. After the command is run, it will suggest commandline arguments
that can be passed to the julia
language program
that will load a precompiled ITensors.jl system image when running Julia. Running ITensor
code this way typically reduces startup times from tens of seconds to about a second
or less.
3 Index Objects
A core concept of the ITensor system is that tensor indices carry information beyond just their dimension. Mathematically, this corresponds to the notion that an index labels the basis of a vector space, and that two vector spaces may be different from each other despite having the same dimension.
The notion that a tensor index corresponds to a specific vector space is encoded in the unique id number assigned to an Index object when it is constructed:
Printing an Index as in the code above shows a portion of the (64 bit) id number.
Because a new id is assigned each time an Index is constructed,
other separately constructed Index objects
will not compare equal to i
even if they have the same dimension
In other words, comparison operations (==
,!=
) require two Index objects to have the same id for them to compare equal.
To enrich the Index system one may also add tags to indices
The Index s
above has a dimension 3, as well as two
tags "s"
and "Site"
.
For efficiency reasons, tags can have a maximum of eight characters and indices can have a maximum of four tags.
Tags can serve multiple purposes: helping to identify Index objects when printing them; collecting subsets of indices sharing a common tag or tags; and preventing certain Index pairs from contracting with each other. This last use of tags extends the rule for Index comparisons: for Index objects to compare equal they must have the same tags as well as the same id number.
As discussed in the previous section, one other way to prevent Index objects from comparing equal is to change their prime level. Every Index carries an integer prime level which defaults to zero.
A copy of Index i
but with
a prime level of 1 can be created by calling
or for convenience by writing
Two copies of the same Index which have different prime levels do not compare equal
Because both primes and tags can be used to
prevent Index objects from comparing equal
to each other and being contracted by the
*
operator, some experience is needed
to choose the best approach. Primes are
useful when indices are only distinguished
temporarily; it is easy afterward to call
noprime(T)
on an ITensor to reset
the prime levels of all of its indices.
On the other hand, tags should be used when
there is some applicationspecific understanding
of why certain indices are distinguished,
such as in a tensor network with a square
lattice structure, where all indices linking the
tensors together may describe the same vector space,
but there may be "Up"
, "Down"
, "Left"
,
and "Right"
copies of the index. This is
particularly useful in applications involving translational
invariance, where many copies of the same Index can appear
in different contexts and it can become cumbersome to distinguish
them by prime levels alone.
4 The ITensor Product Operator
Just as tensor diagrams unify many concepts, the ITensor product operator *
likewise unifies many operations into a single operation:

the
*
product of ITensors with no indices in common computes an outer product 
the
*
product of ITensors with all the same indices computes an inner product, resulting in a scalar ITensor 
otherwise, for a pair of ITensors having just some indices in common, the
*
operator computes a tensor contraction
A simple example of an outer product is the product of two vectors which do not share a common index:
Using the *
operator to compute an inner product results in a scalar ITensor
with no indices as in the following example (note that the indices do not need to be in the same order for the result to be correct):
The scalar
function can be called to convert
a scalar ITensor into a real or complex number
or alternatively one can call x = C[]
.
Finally, to illustrate the case of a tensor contraction where only some of the indices are summed, we can use the following example which was also shown at the beginning of this article:
In the diagram above, we have omitted the names of the indices to emphasize
the typical user experience: all that a user needs to know to get a correct
result in the above example is that T
and M
share one Index.
Keeping track of the ordering of the uncontracted indices, which become the
indices of R
, is not necessary.
Besides contracting regular tensors, the *
operator can also be used
in conjunction with specially constructed tensors to manipulate tensor indices.
One example of such a special tensor type is a delta tensor,
also known as a copy tensor, which has all diagonal
elements equal to one and other elements equal to zero, and is
often shown as a solid black circle in tensor diagrams.
In the ITensor library, a delta tensor uses special diagonalsparse storage
internally, not only to save memory but also to ensure that the contraction of delta
tensors with other tensors is performed using specially optimized routines.
A delta tensor can be used to replace an Index with another Index of the same dimension:
or to duplicate (or split) an index as follows:
Note that in Julia, one can use the unicode character
to write the code above as B = B *
(k,i,j)
.
Another example of a special tensor type is a combiner ITensor. When contracted with another ITensor, a combiner merges multiple indices into a single Index.
The Index c
shown in the diagram above can be retrieved by
calling combinedind(C)
on the combiner ITensor.
Taking the product with the combiner a second time reverses this operation.
Like delta tensors, combiners also use a special storage type with a negligible memory footprint and optimized contraction algorithms for combining and uncombining indices.
5 Tensor Decompositions
Commonly used tensor network decompositions are built from more fundamental matrix decompositions such as the QR and singular value decompositions (SVD) known from linear algebra. Despite being defined in terms of matrices, these factorizations can be straightforwardly defined for tensors too. The only additional step needed is a mapping from a tensor to a matrix, defined by specifying a certain group of indices as row indices and the rest as column indices, then treating each group as a single larger index when computing the decomposition. ITensor automates the tedious and errorprone process of converting tensors to matrices and back, providing a tensorlevel interface for various decompositions.
Consider an ITensor T
with indices i,j,k
. We can compute a QR decomposition of
T
by just specifying that
i,k
are the row indices as follows:
A new Index is generated by the qr
function which
links the Q
tensor to the R
tensor as shown
above. This makes it straightforward to recover the tensor
T
just by using the *
operator:
(In Julia, the operator is overloaded to compute the relative difference
between the two sides of an equation, and return true if it is below a prescribed
threshold.)
Note that when computing the product Q*R
one does not
need to know any details of the new Index introduced by the
factorization, such as whether it is the first or second index
of R, or its dimension. However, in situations where one wants
to retrieve this Index, a convenient way to do it is as follows:
where the commonind
function returns the first Index found that
is shared by the two ITensors.
The SVD plays a key role in tensor network calculations, and is implemented as
In the example above, j,i
were specified as
the row indices, leaving m,k
as the column indices.
An important feature of certain decompositions such as the SVD is that they allow controlled truncation of the tensors resulting from the factorization. By default, ITensor decompositions do not truncate, though they do always compute the “thin” version of a decomposition when available. A truncated decomposition can be computed by specifying truncation keyword arguments. In the following example
the truncation will be determined by summing the squares of the singular values from smallest to largest until the truncation error reaches while also ensuring that the maximum number of singular values kept is less than or equal to 10.
6 Tensor Storage Layer
A powerful feature of ITensor is that ITensors can have a wide variety of storage formats while offering the same user interface. Users can mix sparse and dense tensors together in calculations and manipulate any kind of tensor using identical highlevel code.
In most cases users do not set the storage type manually; instead special storage types occur automatically when using other features: after computing the singular value decomposition of an ITensor, the singular values are returned as an ITensor with diagonalsparse storage; constructing an ITensor from indices with quantum number subspaces makes the storage automatically blocksparse.
Importantly, because the storage types used by an ITensor are
distinct types, each one can use the most optimal memory layout possible,
and performancecritical algorithms such as tensor
contraction and factorization can be specialized for each
storage type or combination of storage types.
For this purpose, we take full advantage of Julia’s multiple
dispatch mechanism to organize specialized algorithms into
separate code pathways to keep the library code simple.
These optimizations are hidden from the user, who can just
contract ITensors together using the *
operation and
automatically get the best possible performance available.
Some of the most common storage types available in ITensor are:

Dense storage: this is the default storage type when constructing an ITensor from regular Index objects and setting elements. The Dense storage type is parameterized over its element type, so that
Dense{Float64}
(realvalued dense storage) andDense{ComplexF64}
(complexvalued dense storage) are actually different storage types. The container type used for Dense storage can also be changed through a second, optional type parameter. 
Diagonal storage
: diagonalsparse tensors occur naturally in algorithms such as the singular value decomposition and eigenvalue decomposition. In such settings, all of the diagonal elements can be different and so an array of the diagonal elements is stored. A special case of diagonal storage is uniform diagonal storage, where all of the elements of the diagonal are constrained to be the same. For this special storage only the value of the repeated, identical diagonal element is stored and speciallyoptimized contraction algorithms are invoked. If the uniform diagonal value is equal to
1.0
then such a diagonal tensor can be used to replace one Index with another under the contraction or*
operation, or as a “copy” or “delta” () tensor as used in certain tensor network algorithms. 
Combiner storage: this storage type uses essentially no memory and stores no tensor components. Rather, it stands for a tensor which conceptually merges two or more indices into one larger index. A combiner tensor
C
can be created asC = combiner(i,j,k)
wherei,j,k
are the indices one wants to combine together. Contracting the combiner ITensor with an ITensor having these indices results in a new ITensor where the indices are merged into the Indexcind = combinedind(C)
. The new combined Index is created automatically by the combiner. 
Blocksparse storage: blocksparse storage is automatically used when an ITensor is created from Index objects with quantumnumber subspaces. This is an important case for quantum physics calculations, where the sparsity enforces symmetries or conservation laws and allows calculations to be performed more efficiently. The block sparse and quantum number system is discussed in more detail in Section 8.

GPU storage: GPU (graphics processing unit) storage is an experimental feature supported by the ITensorsGPU.jl package. An ITensor with GPU storage stores its elements in GPU memory, and calls specialized routines for operations including tensor contraction and tensor factorizations. Taking advantage of the parallel processing capabilities of GPUs can give speedups ranging from two to a hundred times the speed of CPU calculations. Because different storage types are handled automatically behind the same ITensor interface, GPU ITensors can take advantage of the same set of highlevel algorithms available in the ITensor library written originally for regular tensors stored in host memory.
The flexibility of the ITensor storage system will let us explore other interesting
possibilities in the future. One direction we plan to pursue is further sparsity patterns, including
fully general sparsity. Another direction is support for lowerprecision floatingpoint data,
which can significantly speed up calculations such as when using GPU hardware.
A particularly interesting storage type to be investigated is
lazy storage. For an operation such as A*B*C
, a tensor with lazy
storage would simply hold references to A
, B
, C
and only evaluate the result
when triggered by an event such as accessing a tensor element. This would allow
optimizations such as determining the optimal contraction order to be automatically performed
before contractions are evaluated. Note that some lazy ITensor operations are already available
through native Julia features, such as parsing and broadcasting syntax.
7 High Level Features: MPS and MPO Algorithms
To make ITensor a productive system for rapidly prototyping tensor network algorithms, it provides the most common and welldeveloped tensor network formats and algorithms. The two most well developed formats are the matrix product state (MPS) tensor network [Ostlund:1995, Vidal:2003, PerezGarcia:2007], also known as the tensor train [Oseledets:2011], and the matrix product operator (MPO) tensor network [Verstraete:2004d, McCulloch:2007].
Algorithms included with the core ITensor library include summation of MPS and MPO; truncation of MPS and of MPO; optimization of MPS through the DMRG algorithm; and multiplication of an MPS by an MPO. These algorithms offer a high degree of customizability: the multiplication of an MPS by an MPO can be performed using at least three different algorithms (selected by a keyword argument), with each algorithm offering tradeoffs in terms of scaling, performance, and controllability. The DMRG code offers different modes, including finding the ground state (dominant eigenvector) of an implied sum of multiple MPOs or finding excited states (subdominant eigenvectors).
7.1 AutoMPO
A very useful and popular feature of ITensor is the AutoMPO system, which is a domainspecific language for constructing MPO tensor networks from sums of products of local linear operators. Constructing such sums is particularly important for physics applications, where one studies Hamiltonian operators, a typical example being the Heisenberg Hamiltonian:
(1) 
This particular Hamiltonian can be exactly written as an MPO of bond dimension 5, but the construction is technical and tedious to program by hand.
The AutoMPO system automates the construction of this Hamiltonian MPO:
Comparing the lines of code in the for loop above to the Hamiltonian definition in Eq. (1) one can observe a close similarity. The AutoMPO system is very powerful: following a major enhancement of the backend code by Anna Keselman based on Ref. [Chan:2016], AutoMPO can accept terms with more than two local operators and local operators separated by arbitrary distances, and uses an SVDbased compression algorithm to obtain a nearlyoptimal MPO bond dimension. Improvements to this system in the future could include better compression techniques and handling of infinite, translationinvariant systems along the lines of Refs. [Parker:2019, Vanhecke:2019] as well as betterscaling algorithms generalized from techniques used for nonlocal Hamiltonians arising in quantum chemistry [Stoudenmire:2017s].
7.2 DMRG Algorithm
One of the most heavily used highlevel algorithms included with ITensor is the density matrix renormalization group (DMRG) [White:1992, Schollwoeck:2011]. The DMRG algorithm computes lowenergy states of quantum systems, or in mathematical terms, dominant eigenvectors of very large Hermitian linear operators.
The main inputs to a DMRG calculation is a Hamiltonian and an initial guess for its ground state . The ITensor DMRG implementation works generically for any Hamiltonian which can be represented as an MPO tensor network, so that the same code can be applied not only to onedimensional systems, but also quasitwodimensional systems and systems with longrange interactions. By taking advantage of the AutoMPO system discussed above, users can rapidly set up DMRG calculations of complicated Hamiltonians.
Given a Hamiltonian MPO constructed as in Section 7.1 above, one can prepare an initial product state, schedule of sweeps (DMRG algorithm iterations) and accuracy parameters, then run the DMRG algorithm:
For Hamiltonians defined as the sum of different sets of terms one can run a DMRG calculation as:
where H1,H2,H3
are separate MPOs. Instead of
summing these MPOs explicitly, which can be costly and inaccurate,
the algorithm loops over them internally as if they were summed.
This technique can be very helpful in applications such as quantum
chemistry where Hamiltonians can become large and complex, yet have
a nearly blockdiagonal MPO form if represented as a single MPO.
To compute an excited state of a Hamiltonian (subdominant eigenvector)
with ITensor DMRG having first computed both the ground state MPS psi0
,
and first excited state psi1
, say,
one provides [psi0,psi1]
as an extra argument to DMRG, meaning that
the next state computed should be constrained to be orthogonal these previous ones:
In the implementation of this particular DMRG routine, projectors onto the previous states psi0
and psi1
are effectively added to the Hamiltonian times an “energy penalty”, pushing up the energy
of these states in the eigenvalue spectrum so they are no longer part of
the lowenergy subspace. Other techniques for computing excited states might be provided in the future.
7.3 MPS and MPO Operations
Far from being blackbox software for performing calculations with MPS, ITensor provides many elementary building blocks for creating custom algorithms involving MPS, MPOs, and other tensor networks built from these components such as projected entangled pair states (PEPS).
The most elementary interface to MPS and MPO tensor networks involves
retrieving and updating individual factor tensors making up the network.
An MPS is a factorization of a tensor psi
of the following form
(2) 
where we have omitted an explicit sitelabel on each tensor for compactness. The factor tensor on site can be obtained as
and updated as
An important technical step involving an MPS is bringing it into an orthogonal form, where all of the factor tensors to the left or right of the center tensor at a site are equivalent to partial isometries (i.e. either their rows or their columns are orthogonal). To bring an MPS into orthogonal form efficiently in ITensor, one calls:
where we note the convention adopted in Julia programming that functions whose
name end with !
may modify their first argument.
An interesting feature of ITensor MPS objects is that they store information about
which tensors are known to be orthogonal, so that calling
orthogonalize!(psi,j)
repeatedly for the same value of j
does no extra work, and shifting
the orthogonality center of an already partially orthogonalized MPS
can be done with the minimum amount of computation.
Another fundamental operation is truncating an MPS: computing another MPS of
a smaller bond dimension which is as close to the original MPS as possible.
For MPS such a truncation can be done optimally through various deterministic algorithms.
Truncating an MPS psi
in ITensor can be done by calling:
where for the sake of example we have shown specific values of the two most commonly
used truncation parameters. The maxdim
parameter sets an upper limit on the
bond dimension of the MPS after the truncation, whereas the cutoff
parameter
allows the new bond dimension to be determined adaptively as long as the resulting
truncation error remains below the value provided. Using a cutoff can allow the
bond dimension to fall below the maxdim
when possible while still ensuring an accurate
approximation of the original MPS.
ITensor supports arithmetic involving MPS and MPOs to be performed using the
add
function. Performing exact sums can lead to quickly growing costs,
so that one normally truncates the result by providing a truncationerror cutoff.
For example, to add two MPS psi
and phi
one can call:
and similarly for adding two MPOs.
Algorithms such as timeevolving quantum states or contracting twodimensional “PEPS” tensor
networks can be formulated in terms of products of an MPO with MPS or with another MPO.
To approximately multiply an MPS psi
by an MPO W
, one can call the function
with example parameters controlling the truncation shown.
The product of two MPOs R
and W
can also be computed:
Importantly, these functions provide multiple backend algorithm implementations with various tradeoffs in terms of the cost, accuracy, and control offered. For example, to select the accurate yet expensive “naive” algorithm for multiplying an MPS by an MPO one may call
8 Quantum Number BlockSparse ITensors
An important technique used in stateoftheart physics calculations is enforcing constraints on tensors arising from conserved quantities. These are quantities such as total particle number or spin rotation symmetry which become conserved due to symmetries of the Hamiltonian operator. The value of each conserved quantity is known as a quantum number.
Quantum number conservation can be important since physical systems commonly preserve symmetries such as rotational symmetry or particle number conservation symmetry, making it necessary for simulations to conserve these to be comparable to experimental results. And conserving quantum numbers allows calculations to run much faster and use less memory because of a blocksparse structure that is naturally imposed on the tensors in a tensor network.
The power of the ITensor approach to conserving quantum numbers is that quantumnumber conserving, blocksparse ITensors offer nearly the same interface as regular dense ITensors. Algorithms can be written generically for dense ITensors and automatically work for the blocksparse case too.
The design of the ITensor quantum number (QN) system is that QN information is stored in Index objects in a fixed order. This information is queried when an ITensor is constructed to determine whether the storage should be blocksparse and the layout of the blocks. When such ITensors are summed, contracted, or factorized, optimized routines are used and the QN information is propagated to the indices of the resulting ITensor.
As an illustrative example of ITensor’s QN blocksparse system, say we have defined two indices with information about their QN subspaces:
The Index i
has a total dimension of 5 because it has two subspaces, one carrying a quantum
number QN(0)
and of dimension of 2; the other carrying a quantum number QN(1)
and of dimension 3. Similarly j
has a total dimension of 3, coming from its two subspaces.
Using these indices, we can define an ITensor T
as
where emptyITensor
creates an ITensor with unallocated storage and an asyet unspecified pattern of nonzero blocks. Then, we set an element of T
as
Note that this element corresponds to the QN(1)
subspace of i
and the QN(1)
subspace of j
, for a combined “QN flux” of flux(T)==QN(2)
.
When setting any further elements, only those elements of T
consistent with a
flux QN(2)
will be allowed to be nonzero.
This imposes a blocksparse structure on T
, since most values of the indices combine
to form fluxes other than QN(2)
and thus remain zero. Only allowed blocks consistent
with the total flux are stored in memory. Blocksparse computations can then be much more
efficient than with dense tensors because fewer nonzero elements have to be handled and
the presence of disjoint blocks allows major parts of calculations to be performed in parallel.
For the rest of this section, we discuss in more detail the different types composing the QN blocksparse ITensor system.
8.1 QN Objects
Blocksparse ITensors arise from vector spaces which are a direct sum of smaller subspaces. In physics calculations, these subspaces are associated with different quantum numbers. In ITensor, sets of quantum numbers are stored in QN objects as a collection of namevalue pairs, where the value is always an integer. Different values may be combined according to the usual rules of integer addition and subtraction, or according to addition modulo .
QN objects carrying a single quantum number, such as total zcomponent spin "Sz"
, may
be constructed as:
QNs may be added, subtracted, and compared:
QN objects can also carry multiple quantum numbers as follows:
Because the quantum numbers are named, they can be provided to the QN constructor in any order and are sorted internally. For convenience when there is only one quantum number, its name can be omitted; this is equivalent to choosing the name to be the empty string.
Some quantum numbers of physical systems obey a addition rule. A key example is fermion parity, which is only conserved modulo two in systems such as superconductors. A addition rule for a quantum number can be specified by providing as the third entry of the tuple defining that quantum number:
The 2
following the quantum number values above specifies that the "P"
quantum
number obeys addition.
The reason quantum numbers have name strings and are not just distinguished positionally is that having names allows QNs containing different quantum numbers to be combined. This is important to allow users to define different local physical spaces (such as spin versus particle degrees of freedom) in separate codes, then combine them together later. A key example of a physical model combining two otherwise separate types of physical spaces is the HubbardHolstein model, where electron sites are intermixed with boson sites.
8.2 QN Index
As discussed above, the blocksparse structure of quantum number
conserving tensors arises from the directsum structure of the
vector spaces over which they are defined. To specify additional
information about directsum subspaces, an Index
object
can be constructed from QNinteger pairs, as follows:
where note that (a=>b) == Pair(a,b)
is builtin Julia notation for
constructing a pair of values a
and b
.
In the example above, the Index i
has three subspaces,
of dimensions 1, 3, and 2 respectively. Therefore the total dimension
of i
is six, or dim(i)==6
.
The subspaces are associated with the quantum numbers QN("N",0)
,
QN("N",1)
, and QN("N",2)
respectively.
A crucial aspect of QN Index objects not yet discussed is that they have
an Arrow
direction, which can be Out
or In
, with Out
being the default. Mathematically, the direction of an Index says whether it is
covariant (In
) or contravariant (Out
).
The arrows of QN indices play two important roles in working with QN ITensors:

A pair of QN indices must have opposite arrow directions to be contracted.

When computing the QN flux of an ITensor block, QNs corresponding to an
Out
Index are added and QNs corresponding to anIn
Index are subtracted.
Examples of these arrow and flux rules will be given in the next section on QN ITensors.
8.3 QN ITensor
Constructing an ITensor from QN Indices makes it a QN ITensor, with a blocksparse storage type. In addition to the blocksparse real and complex storage types, there are also diagonal blocksparse storage types which are usually obtained from factorizations such as the SVD of blocksparse ITensors.
In most respects, working with QN ITensors is very similar to working with dense ITensors.
Operations like adding QN ITensors or multiplying them by scalars work
in a straightforward way.
However, one small but important difference from dense ITensors arises when contracting QN ITensors:
matching QN indices must have opposite arrow directions to be contracted. This rule is important for
consistent bookkeeping of QN flux under Hermitian conjugation of ITensors.
But because computing the
Hermitian conjugate dag(T)
of a QN ITensor T
is defined to reverse all of the
arrows of its indices, code which is already written correctly for complex, dense ITensors
(with proper use of dag
to handle complex conjugation) will automatically be correct
in terms of QN conservation too.
Having discussed all of the types involved in the QN ITensor system, let us discuss some examples which integrate all of these elements. An example motivated by physics is the Hilbert space of a single “hardcore” boson: a type of particle which cannot share an orbital or site with another boson. Such bosons can be used to model atoms which have large, shortrange repulsive interactions. The Hilbert space of a single hardcore boson is spanned by two basis states and , representing no particle and one particle. Along with these basis states, one can define the elementary operators , , , which lower, raise, or measure the number of particles:
(3) 
Diagrammatically the equation can be expressed as
where in the diagram note that the tensors now have arrows on their indices,
with contracted indices having opposite arrow directions (In
versus Out
).
Within ITensor, we can represent the Hilbert space of this boson as an Index
This Index
is the representation in code of the index lines
in the diagram above.
We can next construct the operator as the following ITensor
The first line constructs a
as an ITensor with indices s’
and dag(s)
with elements all zero, and the second
line sets the only nonzero element of a
. We can visualize the resulting tensor
as follows
The small, red labels above denote the subspaces of the Index s
by labelling them according to the value of the "N"
quantum number.
We can also see the single nonzero element corresponding to the (1,2) entry
of the tensor and having the value 1.0.
A key point about the example of the ITensor for the operator is that the only element stored in memory is the one shown above. All other entries shown in light gray are assumed zero and not stored in memory. To see why this is the case, let us label each block of the tensor (or any tensor having the same indices as ) by its quantumnumber flux:
The nonzero element of the tensor is in the block with flux QN("N",1)
and physically
means this operator always reduces the particle number by 1. Because the convention in ITensor
is that QNconserving ITensors must have a welldefined flux, only blocks with the same flux
are stored in memory and the rest are assumed to be zero and not stored.
In contrast, the operator has a flux of zero, and therefore will have
two allowed blocks: the blocks labeled and shown in blue in the diagram above.
Unlike the examples above, general QNconserving ITensors will have many blocks which can be nonzero, and having block sizes greater than . Summations, factorizations, and especially contractions of general blocksparse tensors can be much faster than for dense tensors with the same index dimensions, not only because the zero (unallocated) blocks can be skipped over, but also because operations on nonzero blocks can be performed in parallel. The C++ implementation of ITensor currently uses multicore parallelism within the algorithm for blocksparse tensor contraction, with speedups of up to observed in practical physics applications, though the speedups vary depending on the application.
Finally, all other operations available for dense tensors work for QNconserving ITensors too, with exactly the same interface. This includes the use of combiner ITensors, factorizations such as the SVD and QR, and higherlevel algorithms involving matrix product states and operators.
9 NDTensors Library
Early on in the design of the ITensor library, a conscious decision was made to separate the highlevel ITensor interface, involving “intelligent” Index objects and related features, from the lowerlevel parts of the code focusing on efficient contraction routines and sparse tensor storage layouts. With the port of the ITensor library to Julia, we have taken this design one step further by making the lowerlevel part of the library a separate package known as NDTensors.jl (Dimensional Tensors) which can be used and developed separately from ITensors.jl.
Some of the goals of developing NDTensors as a separate package include:

separating lowlevel NDTensors algorithms from highlevel ITensor logic simplifies and modularizes the code and prevents bugs

encouraging more community contributions to the ITensor project, since some community members may find the NDTensors interface and features more familiar and appealing, and may not prefer to work with the ITensor layer when making contributions

other software besides ITensor could use NDTensors as a backend, which would promote community efforts to improve tensor software and share resources
9.1 Basic Interface
The NDTensors library is a fullfeatured, standalone library emphasizing generic, highperformance algorithms and support for a variety of sparse tensor types. Unlike the ITensor library, NDTensors has a more traditional interface where users must keep track of the ordering of tensor indices. For example, one can construct a dense tensor with dimensions as
which by default is filled with all zeros and then set its elements as
Tensor objects are 1indexed, similar to Julia arrays. A tensor with complex entries can be constructed as
Contracting two Tensors is done by specifying temporary labels for tensor indices; matching labels indicate two indices are contracted while unique labels denote uncontracted indices. In the following example:
the label 1
of the first index of A
matches the 1
label of the third index of B
, so those two indices are contracted with each
other. Likewise the third index of A
and second of B
share
the label 2
and are contracted. The use of negative integers to label contracted indices
is not required, but is just a convention to make the code more readable.
9.2 Blocksparse Tensors
NDTensors provides sparse tensor types as well. An important example is block sparsity. One way to construct a block sparse tensor is as follows:
The code above specifies that the tensor A
has two indices of dimension 5 (= 2+3) and
9 (= 4+5) respectively, with the first index having two subspaces of dimensions 2 and 3 and the
second index having two subspaces of dimensions 4 and 5. Thus A
has four blocks overall
(because its two indices each have two subspaces). The array nzblocks
lists which
blocks of A
can be nonzero and will be actually allocated in memory, with each tuple giving a subspace number for each index.
9.3 Generic Index Types
A crucial feature of NDTensor is that tensors are allowed to represent their indices not just as a collection of integers or block dimensions (specifying each the dimension of each index), but as any object providing a certain index interface. This generic design allows seamless interoperation between the NDTensors library and the ITensor library, as well as making it easy to provide features such as tensor slicing.
For dense and diag storage, essentially all that is required of the container inds
representing the indices of a Tensor
is that one can call the function dim
on its n
th element.
Examples of valid inds
objects are collections of integers, collections of ITensor Index
objects (provided an overload of the method dim
is provided), Dims
objects
provided by the Julia Base library for indexing builtin Julia tensors,
and BlockDims
objects defined by NDTensors for indexing blocksparse tensors. By default, the strides are determined by the dimensions of the indices, but can be overloaded if needed such as for tensor slicing applications. Unless an explicit set of indices is provided,
Tensor
objects default to using the Dims
type (a tuple of integers) to represent its indices and BlockSparseTensor
objects default to using BlockDims
.
For blocksparse storage types, an overload of the blockdim
function is required for any block index, which is used to query the size of a specified block in a specified dimension.
9.4 Tensor Contraction Backend
Tensor contraction is often the computational bottleneck of tensornetwork algorithms. Thus implementing it as efficiently as possible is critical for performance.
For contracting two dense tensors, NDTensors currently uses a strategy of permuting and reshaping the
tensors into matrices, so that the contraction maps to a matrix multiplication. The motivation
behind this strategy is that BLAS libraries such as Intel MKL offer such high performance that
the extra overhead of permuting the tensors is worthwhile. It is also important to note that
the tensor permutation has a subleading scaling relative to the matrix multiplication, so that in
the limit of large tensors the computation is dominated by the BLAS dgemm
or zgemm
routines.
Though this strategy is a common one for tensor libraries, its implementation in NDTensors is done
carefully to ensure that every case where permutation can be avoided is taken advantage of. Also
if two equivalent strategies exist to permute the contracted tensors to matrices where one of the
permutations is trivial, the code chooses to permute the smaller of the two tensors.
The case of blocksparse tensor contraction reduces to doing a set of smaller, dense tensor contractions on various pairs of blocks from the tensors being contracted. Thus it is built on top of the dense contraction layer of NDTensors, but also offers an excellent opportunity to exploit parallelism, since contraction of the blocks can be done independently, although one does have to handle cases where multiple block pairs contribute to the same block of the resulting tensor. By exploiting multicore parallelism for the same algorithm within the C++ implementation of ITensor we have observed speedups of for DMRG and related MPS calculations (depending on the sparsity and block sizes involved, which varies strongly based on the symmetries used), and up to for tree tensor network calculations. We plan to similarly parallelize the blocksparse contraction algorithm in NDTensors.
Looking ahead, a key improvement to NDTensors will be to offer support for more advanced tensor contraction algorithms that have been recently developed. These algorithms build on sophisticated research into BLAS software, where it was realized that modern BLAS implementations could apply to the case of tensors of arbitrary order, and not just matrices. The two implementations of this type we are aware of are TBLIS [TBLIS] and TCL/GETT [TCL]. These libraries significantly reduce, if not totally eliminate, the permutation overhead inherent to the permutetomatrix strategy discussed above, offering superior performance to the default contraction algorithm of C++ ITensor and NDTensors [tn_benchmarks].
10 Other Features of ITensor
The ITensor library has many other features which are important for productive programming, developing new algorithms or treating new problem domains, but whose precise details are somewhat beyond this highlevel introduction. In this section we briefly highlight these features.
10.1 Writing and Reading ITensor Objects with the HDF5 Format
One important feature is that nearly every type involved in a tensor network, from Index objects to IndexSet’s to ITensors, MPS, and MPOs can be written into and read from HDF5 files. The HDF5 format is a widely used and standardized format for writing large datasets and heterogenous data. It offers portability across operating systems with different binary formats; metadata and a filesystem structure for organizing and retrieving data; and efficient use of memory including compression of numerical data. ITensor objects written to HDF5 files can be both written to and read from both the C++ and Julia versions of ITensor, allowing users with large C++ codes to use the new Julia version for tasks such as performing analysis of simulation results.
10.2 Defining Custom Local Hilbert Spaces
An important feature for physics applications is the ability to define custom “degrees of freedom” or local Hilbert spaces and associated local operators to allow users to implement their own systems of interest within highlevel tools like AutoMPO. ITensor includes builtin definitions for only a handful of common cases such as and spin degrees of freedom, spinless and spinful fermions, and the Hilbert space of the
model. But physics applications of ITensor often call for other definitions, such as of local Hilbert spaces for bosons, higher spin moments such as
, and more exotic degrees of freedom such as parafermions. Users may also want to extend builtin Hilbert space types by defining additional local operators. The C++ version of ITensor already lets users define custom local Hilbert spaces and operators, but due to limitations of the C++ language the customization process has remained cumbersome and users have often had trouble mastering the necessary tasks of defining C++ types, constructors, and methods.Fortunately, in the Julia version of ITensor we have been able to streamline the process of defining and using custom Hilbert spaces. The key innovation is that certain Index tags can be designated as special by defining associated “site types”. For example, say a user wants any Index carrying the tag "S=3/2"
to be interpreted as a spin (the Index should also have the appropriate dimension of 4). Practically this means we want systems such as AutoMPO to know how to make the appropriate local operators such as "Sz"
, "S+"
, and "S"
which act on the Hilbert space of this Index. To tell ITensor how these operators should be defined, a user can create overloads of the ITensors.op!
method which accept a special type: SiteType"S=3/2"
. Examples of such overloads are shown in Listing 1 and can be defined outside the ITensor library in user code.
The notation SiteType"S=3/2"
is a convenient Julia macro syntax which is used to create a
unique type parameterized by a string. Creating types out of values
allows one to effectively overload functions over
different values, even though technically functions can only be overloaded over different types.
After defining these functions, the following code will return "Sz"
, "S+"
, and "S"
operators as ITensors given an Index s
which has the "S=3/2"
tag
The ITensor library reads the tags of the Index passed as the second argument to op
, then checks if any of these tags have an associated SiteType
overload of ITensors.op!
. If exactly one tag and operator name pair does have an ITensors.op!
method defined for it, such as the ::SiteType"S=3/2"
, ::OpName"Sz"
overload in Listing 1 above, then that overload is called to produce the operator corresponding to the requested name as an ITensor. Users can also overload a function called op
which has a similar purpose but which both constructs and returns the operator ITensor, giving more control over the whole process.
What makes this system powerful is that the same op
method and its overloads are called by the AutoMPO system and various MPS and MPO constructors within ITensor library code. So after defining the SiteType"S=3/2"
overloads of the op!
or op
functions above, the following code “just works” and correctly makes an MPO of the Heisenberg Hamiltonian for an site system of spins:
Various special tags with associated SiteType
operator definitions can even be mixed together in Index arrays like the sites
array above, permitting easy setup of calculations for mixed systems such as spin chains of alternating and sites or models of alternating spin and boson sites.
10.3 DMRG Observer System
The DMRG code within ITensor is the most heavily used highlevel feature of the library due to the continued popularity and staying power of the DMRG algorithm. Although ITensor’s implementation of DMRG prints some useful details about the results of each sweep, such as the estimated energy (dominant eigenvalue) and typical bond dimension of the MPS being optimized, there are many situations where a user would like to customize the code further, such as to measure local observables throughout each sweep.
To make this customization process as easy as possible, the ITensor DMRG code accepts an optional observer
keyword argument which allows users to pass any object which is a subtype of AbstractObserver
. This type should also have an overload of at least one of the methods measure!
and checkdone!
defined for it too. These methods can be defined in any way the user sees fit and have minimal requirements. Both are called by the ITensor DMRG code at each step of the DMRG algorithm.
The measure!
method gets passed a variety of properties describing the current state of the DMRG calculation, such as the number of the current sweep and location of the site(s) of the MPS whose local tensors are currently being optimized, and even the entire MPS itself. A customized measure!
function can use this information to produce a detailed snapshot of how the optimization is proceeding. One such use of the observer system in the past was to make animated movies of a DMRG calculation to be used in lectures.
The checkdone!
method can be defined if the user wants to set some criterion for the DMRG calculation to stop before all of the requested sweeps have been performed. Example criteria could include some measure of convergence, such as the energy variance, or an external signal from the user.
11 Applications of ITensor
ITensor has been cited in approximately 310 research articles from 2009 through the first half of 2020.^{3}^{3}3List of papers citing ITensor: https://itensor.org/papers Below we highlight papers which show the diverse applications of ITensor. We expect to see ever wider applications in the future as physics tensor network algorithms become more powerful for two and threedimensional systems, abinitio Hamiltonians, and longtime dynamics [White:2018q, Rakovszky:2020]
, and as more applications of tensor networks are developed in applied mathematics, computer science, and machine learning
[Oseledets:2010, Novikov:2015, pmlrv108rakhshan20a].11.1 Equilibrium Quantum Systems
The most common application area of tensor networks and the ITensor software to date has been equilibrium quantum systems. A common starting point for understanding such systems is through their ground state, and the DMRG algorithm which launched the field of tensor networks is primarily a method for finding ground states. More recently, tensor network methods have been extended to study finitetemperature systems. Another important area of development in the field has been extending DMRG and MPS methods to handle ab initio systems such as in quantum chemistry, where details of continuum, atomic physics must be treated.
An excellent example of a groundstate study using ITensor is that of Keselman and Berg [Keselman:2015], who used ITensor’s DMRG algorithm to compute properties of a onedimensional model of superconducting electrons. A detailed study of properties of finitesize systems, including of quantities at the edge of open systems, supports the existence of a topological state of matter even in the absence of a gap in the excitation spectrum.
The stateoftheart efficiency of ITensor’s DMRG codes make it a powerful tool for studying twodimensional systems as well. DMRG remains one of the leading methods for studying twodimensional quantum systems even though it scales exponentially in the transverse system size. In Refs. [Kallin:2014, Stoudenmire:2014], Kallin, Gustainis, Johal, Stoudenmire, Melko, et al. used a combination of exact diagonalization, numerical linked cluster methods, and ITensor DMRG to obtain entanglement entropies associated with sharp corners in the subsystem geometry for various quantum systems at their critical points. Based on the numerical results, a conjecture was put forward for a universal scaling of this corner entanglement which was afterward supported by field theoretic methods [Bueno:2015].
An exemplary study using ITensor DMRG for a twodimensional system of stronglycorrelated electrons is the work by Venderley and Kim [Venderley:2019], who studied the holedoped Hubbard model on the triangular lattice, finding a transition from pwave to dwave superconductivity as the strength of onsite interactions increase.
ITensor has also been used for studying continuum electronic systems such as quantum chemistry calculations of hydrogen chains [Stoudenmire:2017s, Motta:2017t, White:2019m]; for finitetemperature studies, primarily in the context of the minimally entangled typical thermal state (METTS) algorithm [White:2009, Stoudenmire:2010, Bruognolo:2017, Chen:2019]; and for calculations involving PEPS twodimensional tensor networks [Jiang:2019, Hyatt:2019].
11.2 Dynamics of Quantum Systems
Dynamical behavior of quantum systems or quantum systems outofequilibrium is currently a very active research area, where the flexibility and customizability offered by ITensor has been an excellent fit. Such customizability is important because there are many algorithms available for timeevolving quantum states [Paeckel:2019], most of which are not totally blackbox and require some care to use well. Frontier research problems also involve a variety of settings, such as closed versus open systems, or evolution via Hamiltonians versus circuits, as well as a wide range of measurements to be made of the state.
One paper typifying the use of ITensor for dynamics research, blending numerical results with theoretical predictions, is that of Alba and Calabrese Ref. [Alba:2017], who showed that for integrable systems, such as the XXZ spin chain, one can accurately predict the entanglement entropy at both short and long times.
Nahum, Ruhman, Vijay, and Haah used ITensor in Ref. [Nahum:2017] to simulate dynamics of quantum states evolved by random unitary circuits, supporting their prediction that the growth of entanglement entropy is governed by the KPZ universality class related to the classical statistical physics of surface growth.
Schreiber et al. used ITensor to simulate the dynamics of cold atom experiments in Ref. [Schreiber:2015]
, obtaining good agreement with experimental observations of the difference between the number of atoms in even versus odd minima of the external potential.
A rather different application of dynamical tensor network methods are as “solver” subroutines for the dynamical mean field theory (DMFT) algorithm, which can treat infinitesize systems in two and three dimensions. A novel DMFT solver based on fork tensor network states was proposed and demonstrated using ITensor by Bauernfeind, Zingl, et al. in Ref. [Bauernfeind:2017], allowing DMFT methods to achieve greater resolution for electron spectral functions and other benefits.
11.3 Other Application Areas
Historically tensor network methods have mainly been developed and applied within condensed matter physics. But the recent decade has seen a major broadening in applications of tensor networks inside and outside of physics. These newer applications range from studying holographic dualities between physical theories [Swingle:2012m, Pastawski_2015] to computing highdimensional integrals [Oseledets:2010, Dolgov:2020].
An area where tensor network methods are becoming increasingly important is quantum computing, where they can be used to perform efficient classical simulations of quantum devices. Using matrix product states to represent the state of the device leads to simulations with only polynomial cost for shallow enough circuits or by incurring a gradual loss of fidelity. One effort in this area is the XACC software by McCaskey et al., which offers a matrix product state simulator backend based on ITensor [McCaskey_2018, nguyen2020enabling]. ITensor was also used in a recent study by Zhou et al. which shows that if a quantum computing device is sufficiently noisy, experiments performed with it can be efficiently simulated classically [zhou2020limits].
A rather different application area of tensor networks is applied mathematics and machine learning. Here tensor decomposition methods have found many different uses, from compressing weight layers of neural networks
[Novikov:2015], to recovering missing or corrupted data using partial information [huang2019lowrank]. Machine learning is an area where ITensor has potential to be used much more in the future, and ITensor has already been used to investigate new models and algorithms for machine learning, including supervised [NIPS2016_6211, reyes2020multiscale] and unsupervised [Stoudenmire_2018] learning using models parameterized by tensor networks, and to investigate generalization of these models by studying synthetic data [bradley2019modeling].12 Future Directions
Although the ITensor library already offers high performance and powerful features for implementing any tensor network algorithm, we still envision making many improvements and optimizations in upcoming years. Here we briefly discuss the main features we plan to add, though we emphasize that some may take quite a different form when implemented.
A highpriority new feature is support for automatic differentiation (AD). This technique has been popularized for applications in machine learning and neural networks, but has recently been demonstrated to work well for tensor networks too. For example, AD can be used for stateoftheart infinite PEPS calculations and for calculating critical properties of classical systems [Liao:2019]. The flexible storagetype system already built in to ITensor is ideal for defining differentiation through ITensor types, including MPS and other tensor networks. Existing Julia AD frameworks such as Zygote [Innes:2018] can be used to define a basic ITensor AD system involving ITensors in just a few lines of code, though more work is needed to avoid certain computational or memory bottlenecks specific to tensor algorithms, and we plan to investigate other AD frameworks too.
Another highpriority feature is automatic support for fermionic Hilbert spaces. Systems of fermions are foundational for physics applications of tensor networks, and are the most common type of system studied in condensed matter physics. Currently, the only automatic support for fermions in ITensor is within AutoMPO, which relies on lookup tables of operator names designated as anticommuting. This system works well, but leads to a somewhat confusing experience for users when AutoMPO handles fermions automatically, yet they must then manually introduce JordanWigner string operators for computing certain measurements after a DMRG calculation. We are therefore experimenting with a system that introduces fermionic properties instead at the level of Index objects, where QN subspaces carrying oddparity fermionic quantum numbers have the property that index permutations yield a minus sign if the permutation is odd. Our ambitious goal is for calculations involving fermions to work with exactly the same code as for bosonic degrees of freedom. Even if some manual steps end up being required, this new fermion system could be very useful.
More sophisticated optimizations of tensor contractions are another future direction for ITensor.
One of these is lazy, or delayed evaluation of tensor network expressions such as A*B*C*D
, which can facilitate optimization of the order of the contractions e.g. determining that the optimal order of contraction would be (A*(B*C))*D
automatically without any user input. Another optimization is to support general sparsity, which could be useful for large MPOs or tensors determined from samples of data. Finally, support for nonAbelian quantum numbers would be a very powerful feature for symmetric spin models or quantum chemistry Hamiltonians.
We also plan to offer first class support for the most modern infinite MPS and MPO algorithms. This will include the latest developments in obtaining dominant and subdominant eigenvalues and MPS eigenvectors of infinite MPOs (such as the VUMPS algorithm [ZaunerStauber:2018]), obtaining canonical forms of infinite MPS and MPOs, and applying infinite MPOs to infinite MPS. This will all be offered with the same level of convenience as the currently available finite MPS and MPO methods, including an infinite version of AutoMPO.
We plan to continue to develop GPU support throughout the library. Currently, only dense tensor operations can be performed on GPU, so an initial goal will be to support block sparse tensor operations on GPU. In addition, we plan to make GPU support a first class feature, with convenient syntax for indicating that blocks of code should be performed on GPU. In addition, we plan to port the multithreading of block sparse tensor contractions that are currently available in the C++ version of ITensor to Julia, and additionally expand multithreading support to other block sparse tensor operations, such reshaping and decomposing block sparse tensors.
Last but not least, we hope to offer more highlevel features for PEPS (twodimensional tensor network) calculations. Algorithms and methods for optimizing PEPS have reached a point of maturity such that there are now a handful of essentially standard approaches, such as variational iPEPS [corboz2016variational] and fixedpoint methods for computing PEPS environment tensors [Fishman:2018]. Many of these algorithms could be provided with ITensor in the future.
Acknowledgements
We would like to acknowledge key contributors to ITensor: Katharine Hyatt for developing a GPUaccelerated backend for the ITensors.jl package; ^{4}^{4}4ITensorsGPU: https://github.com/ITensor/ITensorsGPU.jl Anna Keselman for contributing a major improvement to the AutoMPO system which handles longrange interactions and multisite operators; Thomas E. Baker for expanding and improving the ITensor documentation, in particular the tutorials. Thanks also to Jing Chen, YingJer Kao, John Terilla, and Tyler Bryson for discussions about automatically handling fermion signs.
Significant contributions and bug fixes to the C++ version of ITensor were made by Anna Keselman, Mingru Yang, Jack Kemp, Kyungmin Lee, Tatsuto Yamamoto, Juraj Hasik, Benedikt Bruognolo, Jose Lado, Hoi Hui, LarsHendrik Frahm, Lucas Vieira, Markus Wallerberger, Miles Chen, Yevgeny BarLev, Jessica Alfonsi, Chuang Xi, and Andrey Antipov. We would also like to thank Nils Wentzell, Alex Wietek, and Daniel Bauernfeind for their help designing and testing blocksparse multithreading with OpenMP.
Significant features and bug fixes to the initial release of ITensors.jl (the Julia version of ITensor) were contributed by Katharine Hyatt, Ori Alberton, Christopher White, Jan Schneider, Alvaro RubioGarcia, Yiqing Zhou, Michael Abbott, Nicolau Werneck, Michael Sven Ferguson, Nick Robinson, and Amartya Bose.
Funding information
SRW acknowledges the support of the U.S. Department of Energy under grant DESC0008696. The Flatiron Institute is a division of the Simons Foundation.
Appendix A Full Code Examples
a.1 Contraction Example
To show a fully working example of contracting two ITensors with a complicated index structure, consider the following code:
The contraction computed by this code can be expressed by the following diagram:
Note that the Index
tags such as "a"
,"b"
,"c"
, etc. are not required for
this code to function properly, but in this context are just for making the indices
easier to identify when printed.
The line of code
shows the output of the hasinds
function which checks that the ITensor C
has all of the indices a,b,c,i,j
. The code above will output
a.2 DMRG Example
The following code example shows the use of higherlevel features of the ITensor Library to compute the groundstate wavefunction of the Heisenberg quantum spin chain model using the density matrix renormalization group (DMRG) algorithm:
A typical output of this code is:
Brief explanations of the major steps of the above code are:

Construct an array of
Index
objects corresponding to spins (which are dimension2Index
objects labeled by the tag"S=1/2"
). 
Input the terms of the onedimensional Heisenberg Hamiltonian into an
AutoMPO
object. 
Construct an MPO
H
out of theAutoMPO
. 
Construct a random MPS
psi0
of bond dimension 10. 
Create a
Sweeps
struct which indicates that five sweeps of the DMRG algorithm are to be performed, with various maximum bond dimensions allowed for each sweep and a truncation error cutoff of throughout. 
Run the DMRG algorithm, which returns the groundstate energy and groundstate wavefunction MPS.
Appendix B ITensor Implementation and Interface in the C++ Language
In this appendix, we give code examples for the C++ version of ITensor to show the similarities to and differences from the Julia version.
b.1 C++ Contraction Example
Here we show the same example of contracting two ITensors with a complicated index structure as in the previous Appendix section A.1. Consider the following code:
By comparing to the Julia language example A.1, one can see that the C++ code
above is very similar with the main differences being the use of include
statements to import the
library headers, the use of the C++ keyword auto
on lines
of code that result in the definition of a new variable, and semicolons terminating each line of
procedural code.
The last line uses a macro Print
provided by ITensor, which has a similar behavior to the Julia @show
macro and
which in this case generates the output:
b.2 C++ DMRG Example
Here we show the same example of contracting two ITensors with a complicated index structure as in the previous Appendix section A.2. Consider the following code:
By comparing to the Julia language example A.2, one can again see that the codes
are rather similar overall. Some key differences beyond the ones mentioned for the contraction example
include that the site Index
arrays (“site sets”) in the C++ version include quantum number
information by default, which we turn off in this example, and the dmrg
routine outputs
much more information by default, so we pass the named argument {"Quiet=",true}
. These
two parts of the code highlight a custom namedargument system developed for the C++ version of ITensor
which could be more generally useful in other C++ codes and which we plan to release as a separate
library in the future.