Advanced ITensor Usage Guide
Installing and updating ITensors.jl
The ITensors package can be installed with the Julia package manager. Assuming you have already downloaded Julia, which you can get here, from the Julia REPL, type ]
to enter the Pkg REPL mode and run:
$ julia
julia> ]
+
+pkg> add ITensors
Or, equivalently, via the Pkg
API:
julia> import Pkg; Pkg.add("ITensors")
We recommend using ITensors.jl with Intel MKL in order to get the best possible performance. If you have not done so already, you can replace the current BLAS and LAPACK implementation used by Julia with MKL by using the MKL.jl package. Please follow the instructions here.
To use the latest registered (stable) version of ITensors.jl, use update ITensors
in Pkg
mode or import Pkg; Pkg.update("ITensors")
. We will commonly release new patch versions (such as updating from v0.1.12
to v0.1.13
) with bug fixes and improvements. However, make sure to double check before updating between minor versions (such as from v0.1.41
to v0.2.0
) because new minor releases may be breaking.
Remember that if you are compiling system images of ITensors.jl, such as with the ITensors.compile()
command, you will need to rerurn this command to compile the new version of ITensor after an update.
To try the "development branch" of ITensors.jl (for example, if there is a feature or fix we added that hasn't been released yet), you can do add ITensors#main
. You can switch back to the latest released version with add ITensors
. Using the development/main branch is generally not encouraged unless you know what you are doing.
Using ITensors.jl in the REPL
There are many ways you can write code based on ITensors.jl, ranging from using it in the REPL to writing a small script to making a package that depends on it.
For example, you can just start the REPL from your command line like:
$ julia
assuming you have an available version of Julia with the ITensors.jl package installed. Then just type:
julia> using ITensors
and start typing ITensor commands. For example:
julia> i = Index(2, "i")
+(dim=2|id=355|"i")
+
+julia> A = random_itensor(i, i')
+ITensor ord=2 (dim=2|id=355|"i") (dim=2|id=355|"i")'
+NDTensors.Dense{Float64,Array{Float64,1}}
+
+julia> @show A;
+A = ITensor ord=2
+Dim 1: (dim=2|id=355|"i")
+Dim 2: (dim=2|id=355|"i")'
+NDTensors.Dense{Float64,Array{Float64,1}}
+ 2×2
+ 1.2320011464276275 1.8504245734277216
+ 1.0763652402177477 0.030353720156277037
+
+julia> (A*dag(A))[]
+3.9627443142240617
Note that there are some "gotchas" with working in the REPL like this. Technically, all commands in the REPL are in the "global scope". The global scope might not work as you would expect, for example:
julia> for _ in 1:3
+ A *= 2
+ end
+ERROR: UndefVarError: A not defined
+Stacktrace:
+ [1] top-level scope at ./REPL[12]:2
+ [2] eval(::Module, ::Any) at ./boot.jl:331
+ [3] eval_user_input(::Any, ::REPL.REPLBackend) at /home/mfishman/software/julia-1.4.0/share/julia/stdlib/v1.4/REPL/src/REPL.jl:86
+ [4] run_backend(::REPL.REPLBackend) at /home/mfishman/.julia/packages/Revise/AMRie/src/Revise.jl:1023
+ [5] top-level scope at none:0
since the A
inside the for-loop introduces a new local variable. Some alternatives are to wrap that part of the code in a let-block or a function:
julia> function f(A)
+ for _ in 1:3
+ A *= 2
+ end
+ A
+ end
+f (generic function with 1 method)
+
+julia> A = f(A)
+ITensor ord=2 (dim=2|id=355|"i") (dim=2|id=355|"i")'
+NDTensors.Dense{Float64,Array{Float64,1}}
+
+julia> @show A;
+A = ITensor ord=2
+Dim 1: (dim=2|id=355|"i")
+Dim 2: (dim=2|id=355|"i")'
+NDTensors.Dense{Float64,Array{Float64,1}}
+ 2×2
+ 9.85600917142102 14.803396587421773
+ 8.610921921741982 0.2428297612502163
In this particular case, you can alternatively modify the ITensor in-place:
julia> for _ in 1:3
+ A ./= 2
+ end
+
+julia> @show A;
+A = ITensor ord=2
+Dim 1: (dim=2|id=355|"i")
+Dim 2: (dim=2|id=355|"i")'
+NDTensors.Dense{Float64,Array{Float64,1}}
+ 2×2
+ 1.2320011464276275 1.8504245734277216
+ 1.0763652402177477 0.030353720156277037
A common place you might accidentally come across this is when you are creating a Hamiltonian with OpSum
:
julia> using ITensors, ITensorMPS
+
+julia> N = 4;
+
+julia> sites = siteinds("S=1/2",N);
+
+julia> os = OpSum();
+
+julia> for j=1:N-1
+ os += "Sz", j, "Sz", j+1
+ end
+ERROR: UndefVarError: os not defined
+Stacktrace:
+ [1] top-level scope at ./REPL[16]:2
+ [2] eval(::Module, ::Any) at ./boot.jl:331
+ [3] eval_user_input(::Any, ::REPL.REPLBackend) at /home/mfishman/software/julia-1.4.0/share/julia/stdlib/v1.4/REPL/src/REPL.jl:86
+ [4] run_backend(::REPL.REPLBackend) at /home/mfishman/.julia/packages/Revise/AMRie/src/Revise.jl:1023
+ [5] top-level scope at none:0
In this case, you can use os .+= ("Sz", j, "Sz", j+1)
, add!(os, "Sz", j, "Sz", j+1)
, or wrap your code in a let-block or function.
Take a look at Julia's documentation here for rules on scoping. Also note that this behavior is particular to Julia v1.4 and below, and is expected to change in v1.5.
Note that the REPL is very useful for prototyping code quickly, but working directly in the REPL and outside of functions can cause sub-optimal performance. See Julia's performance tips for more information.
We recommend the package OhMyREPL which adds syntax highlighting to the Julia REPL.
Finding documentation interactively
Julia provides many tools for searching for documentation interactively at the REPL. Say that you want to learn more about how to use an ITensor from the command line. You can start by typing ?
followed by ITensor
:
julia> using ITensors
+
+julia> ?ITensor
+search: ITensor ITensors itensor random_itensor
+
+ An ITensor is a tensor whose interface is independent of its
+ memory layout. Therefore it is not necessary to know the ordering
+ of an ITensor's indices, only which indices an ITensor has.
+ Operations like contraction and addition of ITensors automatically
+ handle any memory permutations.
+
+ Examples
+ ≡≡≡≡≡≡≡≡≡≡
+
+ julia> i = Index(2, "i")
+ (dim=2|id=287|"i")
+
+ julia> A = random_itensor(i', i)
+ ITensor ord=2 (dim=2|id=287|"i")' (dim=2|id=287|"i")
+ NDTensors.Dense{Float64,Array{Float64,1}}
+
+ julia> @show A;
+ A = ITensor ord=2
+ Dim 1: (dim=2|id=287|"i")'
+ Dim 2: (dim=2|id=287|"i")
+ NDTensors.Dense{Float64,Array{Float64,1}}
+ 2×2
+ 0.28358594718392427 1.4342219756446355
+ 1.6620103556283987 -0.40952231269251566
+
+ julia> @show inds(A);
+ inds(A) = IndexSet{2} (dim=2|id=287|"i")' (dim=2|id=287|"i")
+[...]
(the specific output may be different for different versions of ITensors.jl as we update the docs). You can use the help prompt (which you get by typing ?
at the REPL) to print out documentation for types and methods.
Another way to get information about types is with the function fieldnames
:
julia> fieldnames(ITensor)
+(:store, :inds)
which shows the fields of a type. Note that in general the specific names of the fields and structures of types may change (we consider those to be internal details), however we often make functions to access the fields of a type that have the same name as the field, so it is a good place to get started. For example, you can access the storage and indices of an ITensor A
with the functions store(A)
and inds(A)
.
Another helpful function is apropos
, which search through all documentation for a string (ignoring the case) and prints a list of all types and methods with documentation that contain the string.
Based on the apropos
function, we can make some helper functions that may be useful. For example:
using ITensors
+
+function finddocs(s)
+ io = IOBuffer()
+ apropos(io, s)
+ v = chomp(String(take!(io)))
+ return split(v, "\n")
+end
+
+function finddocs(s...)
+ intersect(finddocs.(s)...)
+end
+
+found_methods = finddocs("indices", "set difference")
+display(found_methods)
returns:
3-element Array{SubString{String},1}:
+ "ITensors.noncommoninds"
+ "Base.setdiff"
+ "ITensors.uniqueinds"
which are the functions that have docs that contain the strings "indices"
and "set difference"
. We can print the docs for uniqueinds
to find:
help?> uniqueinds
+search: uniqueinds unique_siteinds uniqueind uniqueindex
+
+ uniqueinds(A, B; kwargs...)
+ uniqueinds(::Order{N}, A, B; kwargs...)
+
+
+ Return an IndexSet with indices that are unique to the set of
+ indices of A and not in B (the set difference).
+
+ Optionally, specify the desired number of indices as Order(N),
+ which adds a check and can be a bit more efficient.
We can also filter the results to only specify functions from certain modules, for example:
julia> filter(x -> startswith(x, "ITensors"), finddocs("indices", "set difference"))
+2-element Array{SubString{String},1}:
+ "ITensors.noncommoninds"
+ "ITensors.uniqueinds"
+
+julia> filter(x -> !startswith(x, "ITensors"), finddocs("indices", "set difference"))
+1-element Array{SubString{String},1}:
+ "Base.setdiff"
Ideally we could have apropos
do a "smart" Google-like search of the appropriate docstrings, but this is a pretty good start.
Additionally, the names
function can be useful, which prints the names of all functions and types that are exported by a module. For example:
julia> names(ITensors)
+264-element Array{Symbol,1}:
+ Symbol("@OpName_str")
+ Symbol("@SiteType_str")
+ Symbol("@StateName_str")
+ Symbol("@TagType_str")
+ Symbol("@disable_warn_order")
+ Symbol("@reset_warn_order")
+ Symbol("@set_warn_order")
+ Symbol("@ts_str")
+ :AbstractObserver
+ :OpSum
+ :DMRGObserver
+ :ITensor
+ :ITensors
+ :Index
+[...]
Of course this is a very long list (and the methods are returned as Symbol
s, which are like strings but not as easy to work with). However, we can convert the list to strings and filter the strings to find functions we are interested in, for example:
julia> filter(x -> contains(x, "common") && contains(x, "ind"), String.(names(ITensors)))
+8-element Array{String,1}:
+ "common_siteind"
+ "common_siteinds"
+ "commonind"
+ "commonindex"
+ "commoninds"
+ "hascommoninds"
+ "noncommonind"
+ "noncommoninds"
Julia types do not have member functions, so people coming from object oriented programming languages may find that at first it is more difficult to find methods that are applicable to a certain type. However, Julia has many fantastic tools for introspection that we can use to make this task easier.
Make a small project based on ITensors.jl
Once you start to have longer code, you will want to put your code into one or more files. For example, you may have a short script with one or more functions based on ITensors.jl:
# my_itensor_script.jl
+using ITensors
+
+function norm2(A::ITensor)
+ return (A*dag(A))[]
+end
Then, in the same directory as your script my_itensor_script.jl
, just type:
julia> include("my_itensor_script.jl");
+
+julia> i = Index(2; tags="i");
+
+julia> A = random_itensor(i', i);
+
+julia> norm2(A)
+[...]
As your code gets longer, you can split it into multiple files and include
this files into one main project file, for example if you have two files with functions in them:
# file1.jl
+
+function norm2(A::ITensor)
+ return (A*dag(A))[]
+end
and
# file2.jl
+
+function square(A::ITensor)
+ return A .^ 2
+end
# my_itensor_project.jl
+
+using ITensors
+
+include("file1.jl")
+
+include("file2.jl")
Then, as before, you can use your functions at the Julia REPL by just including the file my_itensor_project.jl
:
julia> include("my_itensor_project.jl");
+
+julia> i = Index(2; tags="i");
+
+julia> A = random_itensor(i', i);
+
+julia> norm2(A)
+[...]
+
+julia> square(A)
+[...]
As your code gets more complicated and has more files, it is helpful to organize it into a package. That will be covered in the next section.
Make a Julia package based on ITensors.jl
In this section, we will describe how to make a Julia package based on ITensors.jl. This is useful to do when your project gets longer, since it helps with:
- Code organization.
- Adding dependencies that will get automatically installed through Julia's package system.
- Versioning.
- Automated testing.
- Code sharing and easier package installation.
- Officially registering your package with Julia.
and many more features that we will mention later.
Start up Julia and install PkgTemplates
$ julia
+
+julia> ]
+
+pkg> add PkgTemplates
then press backspace and type:
julia> using PkgTemplates
+
+julia> t = Template(; user="your_github_username", plugins=[Git(; ssh=true),])
+
+julia> t("MyITensorsPkg")
You should put your Github account name instead of "your_github_username"
, if you want to use Github to host your package. The option plugins=[Git(; ssh=true),]
sets the Github authentication to use ssh, which is generally more convenient. You can switch to https (where you have to type your username and password to push changes) by setting ssh=false
or leaving off plugins=[...]
. By default, the package will be located in the directory ~/.julia/dev
, you can change this with the keyword argument dir=[...]
. However, ~/.julia/dev
is recommended since that is the directory Julia's package manager (and other packages like Revise
) will look for development packages. Please see the PkgTemplate
documentation for more customization options.
Then, we want to tell Julia about our new package. We do this as follows:
julia> ]
+
+pkg> dev ~/.julia/dev/MyITensorsPkg
then you can do:
julia> using MyITensorsPkg
from any directory to use your new package. However, it doesn't have any functions available yet. Additionally, there should be an empty test file already set up here:
~/.julia/dev/MyITensorsPkg/test/runtests.jl
which you can run from any directory like:
julia> ]
+
+pkg> test MyITensorsPkg
It should show something like:
[...]
+Test Summary: |
+MyITensorsPkg.jl | No tests
+ Testing MyITensorsPkg tests passed
since there are no tests yet.
First we want to add ITensors as a dependency of our package. We do this by "activating" our package environment and then adding ITensors:
julia> ]
+
+pkg> activate MyITensorsPkg
+
+(MyITensorsPkg) pkg> add ITensors
This will edit the file ~/.julia/dev/MyITensorsPkg/Project.toml
and add the line
[deps]
+ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
Because your package is under development, back in the main Pkg environment you should type resolve
:
(MyITensorsPkg) pkg> activate
+
+pkg> resolve
Now, if you or someone else uses the package, it will automatically install ITensors.jl for you.
Now your package is set up to develop! Try editing the file ~/.julia/dev/MyITensorsPkg/src/MyITensorsPkg.jl
and add the norm2
function, which calculates the squared norm of an ITensor:
module MyITensorsPkg
+
+using ITensors
+
+export norm2
+
+norm2(A::ITensor) = (A*dag(A))[]
+
+end
The export command makes norm2
available in the namespace without needing to type MyITensorsPkg.norm2
when you do using MyITensorsPkg
. Now in a new Julia session you can do:
julia> using ITensors
+
+julia> i = Index(2)
+(dim=2|id=263)
+
+julia> A = random_itensor(i)
+ITensor ord=1 (dim=2|id=263)
+NDTensors.Dense{Float64,Array{Float64,1}}
+
+julia> norm(A)^2
+6.884457016011188
+
+julia> norm2(A)
+ERROR: UndefVarError: norm2 not defined
+[...]
+
+julia> using MyITensorsPkg
+
+julia> norm2(A)
+6.884457016011188
Unfortunately, if you continue to edit the file MyITensorsPkg.jl
, even if you type using MyITensorsPkg
again, if you are in the same Julia session the changes will not be reflected, and you will have to restart your Julia session. The Revise package will allow you to edit your package files and have the changes reflected in real time in your current Julia session, so you don't have to restart the session.
Now, we can add some tests for our new functionality. Edit the file ~/.julia/dev/MyITensorsPkg/test/runtests.jl
to look like:
using MyITensorsPkg
+using ITensors
+using Test
+
+@testset "MyITensorsPkg.jl" begin
+ i = Index(2)
+ A = random_itensor(i)
+ @test isapprox(norm2(A), norm(A)^2)
+end
Now when you test your package you should see:
pkg> test MyITensorsPkg
+[...]
+Test Summary: | Pass Total
+MyITensorsPkg.jl | 1 1
+ Testing MyITensorsPkg tests passed
Your package should already be set up as a git repository by the PkgTemplates
commands we started with. We recommend using Github or similar versions control systems for your packages, especially if you plan to make them public and officially register them as Julia packages.
You can set up your local package as a Github repository by following the steps here. Many of the steps may be unnecessary since they were already set up by PkgTemplates
. You should be able to go to the website here, create a new Github repository with the name MyITensorsPkg.jl
, and then following the instructions under "push an existing repository from the command line".
You may also want to switch between HTTPS and SSH authentication as described here, if you didn't choose your preferred authentication protocol with PkgTemplates.
There are many more features you can add to your package through various Julia packages and Github, for example:
- Control of precompilation with tools like SnoopCompile.
- Automatic testing of your package at every pull request/commit with Github Actions, Travis, or similar services.
- Automated benchmarking of your package at every pull request with BenchmarkTools, PkgBenchmark and BenchmarkCI.
- Automated building of your documentation with Documenter.
- Compiling your package with PackageCompiler.
- Automatically check what parts of your code your tests check with code coverage.
- Officially register your Julia package so that others can easily install it and follow along with updated versions using the Registrator.
You can take a look at the ITensors Github page for inspiration on setting up some of these services and ideas for organizing your package.
Developing ITensors.jl
This section is for someone who is interested in modifying the source code of ITensors.jl, and then possibly contribute you changes to the official ITensors.jl package.
This should not be necessary for most people. If for whatever reason you think that the functionality of ITensors.jl needs to be modified, oftentimes you can add new functions outside of ITensors.jl or directly overload a function of ITensors.jl (for example with the import keyword).
However, if you would like to only modify parts of the internals of an ITensors.jl function, and/or plan to contribute changes like bug fixes or new features to the official ITensors.jl package, this section is for you.
If you install a package like ITensors with the package manager using the standard Pkg.add
command:
julia> using Pkg
+
+julia> Pkg.add("ITensors")
it will automatically clone the latest registered/tagged version of ITensors
in a randomly generated directory inside ~/.julia/packages
. You can find out what version you are using with Pkg.status
:
julia> Pkg.status("ITensors")
+ Status `~/.julia/environments/v1.7/Project.toml`
+ [9136182c] ITensors v0.2.16
and you can use pkgdir
to find out the directory of the source code of a package that you have loaded:
julia> using ITensors
+
+julia> pkgdir(ITensors)
+"/home/mfishman/.julia/packages/ITensors/cu9Bo"
The source code of a package loaded in this way is read-only, so you won't be able to modify it.
If you want to modify the source code of ITensors.jl
, you should check out the packages NDTensors.jl
and ITensors.jl
in development mode with Pkg.develop
:
julia> Pkg.develop(["NDTensors", "ITensors"])
+Path `/home/mfishman/.julia/dev/ITensors` exists and looks like the correct repo. Using existing path.
+ Resolving package versions...
+ Updating `~/.julia/environments/v1.7/Project.toml`
+ [9136182c] ~ ITensors v0.2.16 ⇒ v0.2.16 `~/.julia/dev/ITensors`
+ [23ae76d9] ~ NDTensors v0.1.35 ⇒ v0.1.35 `~/.julia/dev/ITensors/NDTensors`
+ Updating `~/.julia/environments/v1.7/Manifest.toml`
+ [9136182c] ~ ITensors v0.2.16 ⇒ v0.2.16 `~/.julia/dev/ITensors`
+ [23ae76d9] ~ NDTensors v0.1.35 ⇒ v0.1.35 `~/.julia/dev/ITensors/NDTensors`
+
+julia> Pkg.status(["NDTensors", "ITensors"])
+ Status `~/.julia/environments/v1.7/Project.toml`
+ [9136182c] ITensors v0.2.16 `~/.julia/dev/ITensors`
+ [23ae76d9] NDTensors v0.1.35 `~/.julia/dev/ITensors/NDTensors`
Then, Julia will use the version of ITensors.jl
living in the directory ~/.julia/dev/ITensors
and the version of NDTensors.jl
living in the directory ~/.julia/dev/ITensors/NDTensors
, though you may need to restart Julia for this to take affect.
We recommend checking out the development versions of both NDTensors.jl
and ITensors.jl
since we often develop both packages tandem, so the development branch of ITensors.jl
may rely on changes we make in NDTensors.jl
.
By default, when you modify code in ~/.julia/dev/ITensors
or ~/.julia/dev/ITensors/NDTensors
you will need to restart Julia for the changes to take affect. A way around this issue is the Revise package. We highly recommend using the Revise package when you are developing packages, which automatically detects changes you are making to a package you have checked out for development and edit code and not have to restart your Julia session. In short, if you have Revise.jl
loaded, you can edit the code in ~/.julia/dev/ITensors
or ~/.julia/dev/ITensors/NDTensors
and the changes you make will be reflected on the fly as you use the package (there are some limitations, for example you will need to restart Julia if you change the definitions of types).
Note that the code in ~/.julia/dev/ITensors
is just a git repository cloned from the repository https://github.com/ITensor/ITensors.jl, so you can do anything that you would with any other git repository (use forks of the project, check out branches, push and pull changes, etc.).
The standard procedure for submitting a bug fix or new feature to ITensors.jl would then be to first fork the ITensors.jl repository. Then, check out your fork for development with:
julia> using Pkg
+
+julia> Pkg.develop(url="https://github.com/mtfishman/ITensors.jl")
where you would replace mtfishman
with your own Github username. Make the changes to the code in ~/.julia/dev/ITensors
, push the changes to your fork, and then make a pull request to the ITensors.jl Github repository.
To go back to the official version of the NDTensors.jl
and ITensors.jl
packages, you can use the command Pkg.free(["NDTensors", "ITensors"])
:
julia> Pkg.free(["NDTensors", "ITensors"])
+ Resolving package versions...
+ Updating `~/.julia/environments/v1.7/Project.toml`
+ [9136182c] ~ ITensors v0.2.16 `~/.julia/dev/ITensors` ⇒ v0.2.16
+ [23ae76d9] ~ NDTensors v0.1.35 `~/.julia/dev/ITensors/NDTensors` ⇒ v0.1.35
+ Updating `~/.julia/environments/v1.7/Manifest.toml`
+ [9136182c] ~ ITensors v0.2.16 `~/.julia/dev/ITensors` ⇒ v0.2.16
+ [23ae76d9] ~ NDTensors v0.1.35 `~/.julia/dev/ITensors/NDTensors` ⇒ v0.1.35
+
+julia> Pkg.status(["NDTensors", "ITensors"])
+ Status `~/.julia/environments/v1.7/Project.toml`
+ [9136182c] ITensors v0.2.16
+ [23ae76d9] NDTensors v0.1.35
so it returns to the version of the package you would have just after installing with Pkg.add
.
Some of the Julia package development workflow definitely takes some getting used to, but once you figure out the "flow" and have a picture of what is going on there are only a small set of commands you really need to use.
A small note is that we follow the Blue style guide for formatting the source code in ITensors.jl. To make this more automated, we use the wonderful package JuliaFormatter.jl. To format your developed version of ITensors.jl, all you have to do is change your directory to ~/.julia/dev/ITensors
and run the command format(".")
after loading the JuliaFormatter
package:
julia> using Pkg
+
+julia> Pkg.status("ITensors")
+ Status `~/.julia/environments/v1.7/Project.toml`
+ [9136182c] ITensors v0.2.16 `~/.julia/dev/ITensors`
+
+julia> using ITensors
+
+julia> pkgdir(ITensors)
+"/home/mfishman/.julia/dev/ITensors"
+
+julia> cd(pkgdir(ITensors))
+
+julia> using JuliaFormatter
+
+julia> format(".")
+false
+
+julia> format(".") # Check the formatting succeeded
+true
This will automatically change the style of the code according to the Blue
style guide. The format
command returns false
if the code was not already formatted (and therefore if the command made changes to the source code to follow the style guide), and returns true
otherwise.
If you make changes to ITensors that you think will be useful to others, such as fixing bugs or adding new features, please consider making a pull request. However, please ask us first before doing so – either by raising an issue on Github or asking a question on the ITensor support forum – to make sure it is a change or addition that we will want to include or to check that it is not something we are currently working on. Coordinating with us in that way will help save your time and energy as well as ours!
Here is a great introduction to Julia package development as well as making pull requests to existing Julia packages by the irreplacable Chris Rackauckas.
Compiling ITensors.jl
You might notice that the time to load ITensors.jl (with using ITensors
) and the time to run your first few ITensor commands is slow. This is due to Julia's just-in-time (JIT) compilation. Julia is compiling special versions of each function that is being called based on the inputs that it gets at runtime. This allows it to have fast code, often nearly as fast as fully compiled languages like C++, while still being a dynamic language.
However, the long startup time can still be annoying. In this section, we will discuss some strategies that can be used to minimize this annoyance, for example:
- Precompilation.
- Staying in the same Julia session with Revise.
- Using PackageCompiler to compile ITensors.jl ahead of time.
Precompilation is performed automatically when you first install ITensors.jl or update a version and run the command using ITensors
for the first time. For example, when you first use ITensors after installation or updating, you will see:
julia> using ITensors
+[ Info: Precompiling ITensors [9136182c-28ba-11e9-034c-db9fb085ebd5]
The process is done automatically, and puts some compiled binaries in your ~/.julia
directory. The goal is to decrease the time it takes when you first type using ITensors
in your next Julia session, and also the time it takes for you to first run ITensor functions in a new Julia session. This helps the startup time, but currently doesn't help enough. This is something both ITensors.jl and the Julia language will try to improve over time.
To avoid this time, it is recommended that you work as much as you can in a single Julia session. You should not need to restart your Julia session very often. For example, if you are writing code in a script, just include
the file again which will pull in the new changes to the script (the exception is if you change the definition of a type you made, which would requiring restarting the REPL).
If you are working on a project, we highly recommend using the Revise package which automatically detects changes you are making in your packages and reflects them real-time in your current REPL session. Using these strategies should minimize the number of times you need to restart your REPL session.
If you plan to use ITensors.jl directly from the command line (i.e. not from the REPL), and the startup time is an issue, you can try compiling ITensors.jl using PackageCompiler.
Before using PackageCompiler to compile ITensors, when we first start using ITensors.jl we might see:
julia> @time using ITensors
+ 3.845253 seconds (10.96 M allocations: 618.071 MiB, 3.95% gc time)
+
+julia> @time i = Index(2);
+ 0.000684 seconds (23 allocations: 20.328 KiB)
+
+julia> @time A = random_itensor(i', i);
+ 0.071022 seconds (183.24 k allocations: 9.715 MiB)
+
+julia> @time svd(A, i');
+ 5.802053 seconds (24.56 M allocations: 1.200 GiB, 7.83% gc time)
+
+julia> @time svd(A, i');
+ 0.000177 seconds (450 allocations: 36.609 KiB)
ITensors provides the command ITensors.compile()
to create what is called a "custom system image", a custom version of Julia that includes a compiled version of ITensors (see the PackageCompiler documentation for more details). Just run the command:
julia> ITensors.compile()
+[...]
By default, this will create the file sys_itensors.so
in the directory ~/.julia/sysimages
. Then if we start julia with:
$ julia --sysimage ~/.julia/sysimages/sys_itensors.so
then you should see something like:
julia> @time using ITensors
+ 0.330587 seconds (977.61 k allocations: 45.807 MiB, 1.89% gc time)
+
+julia> @time i = Index(2);
+ 0.000656 seconds (23 allocations: 20.328 KiB)
+
+julia> @time A = random_itensor(i', i);
+ 0.000007 seconds (7 allocations: 576 bytes)
+
+julia> @time svd(A, i');
+ 0.263526 seconds (290.02 k allocations: 14.220 MiB)
+
+julia> @time svd(A, i');
+ 0.000135 seconds (350 allocations: 29.984 KiB)
which is much better.
Note that you will have to recompile ITensors with the command ITensors.compile()
any time that you update the version of ITensors in order to keep the system image updated. We hope to make this process more automated in the future.
Benchmarking and profiling
Julia has great built-in tools for benchmarking and profiling. For benchmarking fast code at the command line, you can use BenchmarkTools:
julia> using ITensors;
+
+julia> using BenchmarkTools;
+
+julia> i = Index(100, "i");
+
+julia> A = random_itensor(i, i');
+
+julia> @btime 2*$A;
+ 4.279 μs (8 allocations: 78.73 KiB)
We recommend packages like ProfileView to get detailed profiles of your code, in order to pinpoint functions or lines of code that are slower than they should be.
ITensor type design and writing performant code
Advanced users might notice something strange about the definition of the ITensor type, that it is often not "type stable". Some of this is by design. The definition for ITensor is:
mutable struct ITensor
+ inds::IndexSet
+ store::TensorStorage
+end
These are both abstract types, which is something that is generally discouraged for peformance.
This has a few disadvantages. Some code that you might expect to be type stable, like getindex
, is not, for example:
julia> i = Index(2, "i");
+
+julia> A = random_itensor(i, i');
+
+julia> @code_warntype A[i=>1, i'=>2]
+Variables
+ #self#::Core.Compiler.Const(getindex, false)
+ T::ITensor
+ ivs::Tuple{Pair{Index{Int64},Int64}}
+ p::Tuple{Union{Nothing, Int64}}
+ vals::Tuple{Any}
+
+Body::Number
+1 ─ %1 = NDTensors.getperm::Core.Compiler.Const(NDTensors.getperm, false)
+│ %2 = ITensors.inds(T)::IndexSet{1,IndexT,DataT} where DataT<:Tuple where IndexT<:Index
+│ %3 = Base.broadcasted(ITensors.ind, ivs)::Base.Broadcast.Broadcasted{Base.Broadcast.Style{Tuple},Nothing,typeof(ind),Tuple{Tuple{Pair{Index{Int64},Int64}}}}
+│ %4 = Base.materialize(%3)::Tuple{Index{Int64}}
+│ (p = (%1)(%2, %4))
+│ %6 = NDTensors.permute::Core.Compiler.Const(NDTensors.permute, false)
+│ %7 = Base.broadcasted(ITensors.val, ivs)::Base.Broadcast.Broadcasted{Base.Broadcast.Style{Tuple},Nothing,typeof(val),Tuple{Tuple{Pair{Index{Int64},Int64}}}}
+│ %8 = Base.materialize(%7)::Tuple{Int64}
+│ (vals = (%6)(%8, p))
+│ %10 = Core.tuple(T)::Tuple{ITensor}
+│ %11 = Core._apply_iterate(Base.iterate, Base.getindex, %10, vals)::Number
+│ %12 = Core.typeassert(%11, ITensors.Number)::Number
+└── return %12
+
+julia> typeof(A[i=>1, i'=>2])
+Float64
Uh oh, that doesn't look good! Julia can't know ahead of time, based on the inputs, what the type of the output is, besides that it will be a Number
(though at runtime, the output has a concrete type, Float64
).
So why is it designed this way? The main reason is to allow more generic and dynamic code than traditional, statically-typed Arrays. This allows us to have code like:
julia> i = Index(2, "i")
+(dim=2|id=811|"i")
+
+julia> A = ITensor(i', i);
+
+julia> @show A;
+A = ITensor ord=2
+Dim 1: (dim=2|id=811|"i")'
+Dim 2: (dim=2|id=811|"i")
+NDTensors.Empty{Float64,NDTensors.Dense{Float64,Array{Float64,1}}}
+ 2×2
+
+
+
+julia> A[i' => 1, i => 2] = 1.2;
+
+julia> @show A;
+A = ITensor ord=2
+Dim 1: (dim=2|id=811|"i")'
+Dim 2: (dim=2|id=811|"i")
+NDTensors.Dense{Float64,Array{Float64,1}}
+ 2×2
+ 0.0 1.2
+ 0.0 0.0
Here, the type of the storage of A is changed in-place. It starts as an Empty
storage, a special trivial storage. When we set an element, we then allocate the appropriate storage. Allocations are performed only when needed, so if another element is set then no allocation is performed. More generally, this allows ITensors to have more generic in-place functionality, so you can write code where you don't know what the storage is until runtime.
This can lead to certain types of code having perfomance problems, for example looping through ITensors with many elements can be slow:
julia> function myscale!(A::ITensor, x::Number)
+ for n in 1:dim(A)
+ A[n] = x * A[n]
+ end
+ end;
+
+julia> d = 10_000;
+
+julia> i = Index(d);
+
+julia> @btime myscale!(A, 2) setup = (A = random_itensor(i));
+ 2.169 ms (117958 allocations: 3.48 MiB)
However, this is fast:
julia> function myscale!(A::Array, x::Number)
+ for n in 1:length(A)
+ A[n] = x * A[n]
+ end
+ end;
+
+julia> @btime myscale!(A, 2) setup = (A = randn(d));
+ 3.451 μs (0 allocations: 0 bytes)
+
+julia> myscale2!(A::ITensor, x::Number) = myscale!(array(A), x)
+myscale2! (generic function with 1 method)
+
+julia> @btime myscale2!(A, 2) setup = (A = random_itensor(i));
+ 3.571 μs (2 allocations: 112 bytes)
How does this work? It relies on a "function barrier" technique. Julia compiles functions "just-in-time", so that calls to an inner function written in terms of a type-stable type are still fast. That inner function is compiled to very fast code. The main overhead is that Julia has to determine which function to call at runtime.
Therefore, users should keep this in mind when they are writing ITensors.jl code, and we warn that explicitly looping over large ITensors by individual elements should be done with caution in performance critical sections of your code. However, be sure to benchmark and profile your code before prematurely optimizing, since you may be surprised about what are the fast and slow parts of your code.
Some strategies for avoiding ITensor loops are:
- Use broadcasting and other built-in ITensor functionality that makes use of function barriers.
- Convert ITensors to type-stable collections like the Tensor type of NDTensors.jl and write functions in terms of the Tensor type (i.e. the function barrier techique that is used throughout ITensors.jl).
- When initializing very large ITensors elementwise, use built-in ITensor constructors, or first construct an equivalent tensor as an Array or Tensor and then convert it to an ITensor.
ITensor in-place operations
In-place operations can help with optimizing code, when the memory of the output tensor of an operation is preallocated.
The main way to access this in ITensor is through broadcasting. For example:
A = random_itensor(i, i')
+B = random_itensor(i', i)
+A .+= 2 .* B
Internally, this is rewritten by Julia as a call to broadcast!
. ITensors.jl overloads this call (or more specifically, a lower level function copyto!
written in terms of a special lazy type that saves all of the objects and operations). Then, this call is rewritten as
map!((x,y) -> x+2*y, A, A, B)
This is mostly an optimization to use when you can preallocate storage that can be used multiple times.
Additionally, ITensors makes the unique choice that:
C .= A .* B
is interpreted as an in-place tensor contraction. What this means is that this calls a function:
mul!(C, A, B)
(likely to be given an alternative name contract!
) which contracts A
and B
into the pre-allocated memory C
.
Because of the design of the ITensor type (see the section above), there is some flexibility we take in allocating memory for users. For example, if the storage type is more narrow than the result, for convenience we might expand it in-place. If you are worried about memory allocations, we recommend using benchmarking and profiling to pinpoint slow parts of your code (often times, you may be surprised by what is actually slow).
NDTensors and ITensors
ITensors.jl is built on top of another, more traditional tensor library called NDTensors. NDTensors implements AbstractArrays with a variety of sparse storage types, with more to come in the future.
NDTensors implements functionality like permutation of dimensions, fast get and set index, broadcasting, and tensor contraction (where labels of the dimensions must be specified).
For example:
using ITensors
+using NDTensors
+
+T = Tensor(2,2,2)
+T[1,2,1] = 1.3 # Conventional element setting
+
+i = Index(2)
+T = Tensor((i,i',i')) # The identifiers are ignored, just interpreted as above
+T[1,2,1] = 1.3
To make performant ITensor code (refer to the the previous section on type stability and function barriers), ITensor storage data and indices are passed by reference into Tensors, where the performance critical operations are performed.
An example of a function barrier using NDTensors is the following:
julia> using NDTensors
+
+julia> d = 10_000;
+
+julia> i = Index(d);
+
+julia> function myscale!(A::Tensor, x::Number)
+ for n in 1:dim(A)
+ A[n] = x * A[n]
+ end
+ end;
+
+julia> @btime myscale!(A, 2) setup = (A = Tensor(d));
+ 3.530 μs (0 allocations: 0 bytes)
+
+julia> myscale2!(A::ITensor, x::Number) = myscale!(tensor(A), x)
+myscale2! (generic function with 1 method)
+
+julia> @btime myscale2!(A, 2) setup = (A = random_itensor(i));
+ 3.549 μs (2 allocations: 112 bytes)
A very efficient function is written for the Tensor type. Then, the ITensor version just wraps the Tensor function by calling it after converting the ITensor to a Tensor (without any copying) with the tensor
function. This is the basis for the design of all performance critical ITensors.jl functions.