Skip to content

Commit

Permalink
Corbett/types and docs (#9)
Browse files Browse the repository at this point in the history 
* Specify type of MPO and pass OpCacheVec as strings

* Docs and examples
  • Loading branch information
@corbett5
corbett5 authored Feb 27, 2024
1 parent e0cb5bf commit 637dd24
Show file tree
Hide file tree
Showing 7 changed files with 423 additions and 150 deletions.
135 changes: 60 additions & 75 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@
[![Coverage](https://codecov.io/gh/ITensor/ITensorMPOConstruction.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/ITensor/ITensorMPOConstruction.jl)
[![Code Style: Blue](https://img.shields.io/badge/code%20style-blue-4495d1.svg)](https://github.com/invenia/BlueStyle)

A fast algorithm for constructing a Matrix Product Operator (MPO) from a sum of local operators. This is a replacement for `MPO(os::OpSum, sites::Vector{<:Index})`. In all cases examined so far this algorithm constructs an MPO with a smaller (or equal) bond dimension faster than the competition.
A fast algorithm for constructing a Matrix Product Operator (MPO) from a sum of local operators. This is a replacement for `MPO(os::OpSum, sites::Vector{<:Index})`. In all cases examined so far this algorithm constructs an MPO with a smaller (or equal) bond dimension faster than the competition. All runtimes below are taken from a single sample on a 2021 MacBook Pro with the M1 Max CPU and 32GB of memory.

## Installation

Expand All @@ -17,15 +17,15 @@ julia> using Pkg; Pkg.add(url="https://github.com/ITensor/ITensorMPOConstruction

## Constraints

This algorithm shares same constraints as ITensor's default algorithm.
This algorithm shares the same constraints as ITensors' default algorithm.

* The operator must be expressed as a sum of products of single site operators. For example a CNOT could not appear in the sum since it is a two site operator.
* When dealing with Fermionic systems the parity of each term in the sum must be even. That is the combined number of creation and annihilation operators must be even.
1. The operator must be expressed as a sum of products of single site operators. For example a CNOT could not appear in the sum since it is a two site operator.
2. When dealing with Fermionic systems the parity of each term in the sum must be even. That is the combined number of creation and annihilation operators in each term in the sum must be even.

There are also two additional constraints:

* Each term in the sum of products representation can only have a single operator acting on a site. For example a term such as $\mathbf{X}^{(1)} \mathbf{X}^{(1)}$ is not allowed. However, there is a pre-processing utility that can automatically replace $\mathbf{X}^{(1)} \mathbf{X}^{(1)}$ with $\mathbf{I}^{(1)}$. This is not a hard requirement for the algorithm, just a simplification to improve performance. Many operators of interest can be easily expressed in a form where only a single operator acts on each site in a term.
* When constructing a quantum number (QN) conserving operator the total flux of the operator must be zero. It would be trivial to remove this constraint.
3. Each term in the sum of products representation can only have a single operator acting on a site. For example a term such as $\mathbf{X}^{(1)} \mathbf{X}^{(1)}$ is not allowed. However, there is a pre-processing utility that can automatically replace $\mathbf{X}^{(1)} \mathbf{X}^{(1)}$ with $\mathbf{I}^{(1)}$. This is not a hard requirement for the algorithm but a simplification to improve performance.
4. When constructing a quantum number conserving operator the total flux of the operator must be zero. It would be trivial to remove this constraint.

## `MPO_new`

Expand All @@ -36,8 +36,22 @@ function MPO_new(os::OpSum, sites::Vector{<:Index}; kwargs...)::MPO
```

The optional keyword arguments are
* `tol::Real`: The tolerance used in the sparse QR decomposition. The default value is calculated separately for each QR decomposition, it is almost always a good value.
* `basisOpCacheVec::OpCacheVec`: A list of operators to use as a basis for each site. The operators on at each site are expressed as one of these basis operators.
* `basisOpCacheVec`: A list of operators to use as a basis for each site. The operators on each site are expressed as one of these basis operators.
* `tol::Real`: The tolerance used in the sparse QR decomposition (which is done by SPQR). The default value is the SPQR default which is calculated separately for each QR decomposition. If you want a MPO that is accurate up to floating point errors the default tolerance should work well. If instead you want to compress the MPO the value `tol` will differ from the `cutoff` passed to `ITensor.MPO` since the truncation method is completely different. If you want to replicate the same truncation behavior first construct the MPO with a suitably small (or default) `tol` and then use `ITensors.truncate!`.

## Examples: Fermi-Hubbard Hamiltonian in Real Space

The one dimensional Fermi-Hubbard Hamiltonian with periodic boundary conditions on $N$ sites can be expressed in real space as

$$
\mathcal{H} = -t \sum_{i = 1}^N \sum_{\sigma \in (\uparrow, \downarrow)} \left( c^\dagger_{i, \sigma} c_{i + 1, \sigma} + c^\dagger_{i, \sigma} c_{i - 1, \sigma} \right) + U \sum_{i = 1}^N n_{i, \uparrow} n_{i, \downarrow} \ ,
$$

where the periodic boundary conditions enforce that $c_k = c_{k + N}$. For this Hamiltonian all that needs to be done to switch over to using ITensorMPOConstruction is switch `MPO(os, sites)` to `MPO_New(os, sites)`.

https://github.com/ITensor/ITensorMPOConstruction.jl/blob/main/examples/fermi-hubbard.jl#L4-L24

For $N = 1000$ both ITensors and ITensorMPOConstruction can construct an MPO of bond dimension 10 in under two seconds. For a compelling reason to use ITensorMPOConstruction we need to look at a more complicated Hamiltonian.

## Examples: Fermi-Hubbard Hamiltonian in Momentum Space

Expand All @@ -47,86 +61,57 @@ $$
\mathcal{H} = \sum_{k = 1}^N \epsilon(k) \left( n_{k, \downarrow} + n_{k, \uparrow} \right) + \frac{U}{N} \sum_{p, q, k = 1}^N c^\dagger_{p - k, \uparrow} c^\dagger_{q + k, \downarrow} c_{q, \downarrow} c_{p, \uparrow}
$$

where $\epsilon(k) = -2 t \cos(\frac{2 \pi k}{N})$ and $c_k = c_{k + N}$. Below is a plot of the bond dimension of the MPO produced by ITensors' default algorithm, [Renormalizer](https://github.com/shuaigroup/Renormalizer) which uses the [bipartite-graph algorithm](https://doi.org/10.1063/5.0018149), and `ITensorMPOConstruction`.
where $\epsilon(k) = -2 t \cos(\frac{2 \pi k}{N})$ and $c_k = c_{k + N}$. The code to construct the `OpSum` is shown below.

![](./docs/plot-generators/fh.png)
https://github.com/ITensor/ITensorMPOConstruction.jl/blob/main/examples/fermi-hubbard.jl#L26-L49

Of note is that the bond dimension of the MPO produced by Renormalizer scales as $O(N^2)$, both ITensors and ITensorMPOConstruction however produce an MPO with a bond dimension that scales as $O(N)$. Below is a table of the time it took to construct the MPO for various number of sites. Some warm up was done for the Julia calculations to avoid measuring compilation overhead. Data recorded on 2021 MacBook Pro with the M1 Max CPU and 32GB of memory.
Unlike the previous example, some more involved changes will be required to use ITensorMPOConstruction. This is because the `OpSum` has multiple operators acting on the same site, violating constraint #3. For example when $k = 0$ in the second loop we have terms of the form $c^\dagger_{p, \uparrow} c^\dagger_{q, \downarrow} c_{q, \downarrow} c_{p, \uparrow}$. You could always create a special case for $k = 0$ and rewrite it as $n_{p, \uparrow} n_{q, \downarrow}$. However if using "Electron" sites then you would also need to consider other cases such as when $p = q$, this would introduce a lot of extraneous complication. Luckily ITensorMPOConstruction provides a method to automatically perform these transformations. If you provide a set of operators for each site to `MPO_new` it will attempt to express the operators acting on each site as a single one of these "basis" operators. The code to do this is shown below.

| $N$ | ITensors | Renormalizer | ITensorMPOConstruction |
|-----|----------|--------------|------------------------|
| 10 | 0.35s | 0.26 | 0.03s |
| 20 | 27s | 3.4s | 0.15s |
| 30 | N/A | 17s | 0.41s |
| 40 | N/A | 59s | 0.93s |
| 50 | N/A | 244s | 1.8s |
| 100 | N/A | N/A | 20s |
| 200 | N/A | N/A | 310s |
https://github.com/ITensor/ITensorMPOConstruction.jl/blob/main/examples/fermi-hubbard.jl#L51-L76

### `OpIDSum`

The code for this example can be found in [examples/fermi-hubbard.jl](https://github.com/ITensor/ITensorMPOConstruction.jl/blob/main/examples/fermi-hubbard.jl). Just like with ITensors, the terms in the Hamiltonian are put into an `OpSum` and then the `OpSum` is transformed into an `MPO`. The code to produce the `OpSum` is
For $N = 200$ constructing the `OpSum` takes 42s and constructing the MPO from the `OpSum` with ITensorMPOConstruction takes another 306s. For some systems constructing the `OpSum` can actually be the bottleneck. In these cases you can construct an `OpIDSum` instead.

```julia
os = OpSum{Float64}()
for k in 1:N
epsilon = -2 * t * cospi(2 * k / N)
os .+= epsilon, "Nup", k
os .+= epsilon, "Ndn", k
end

for p in 1:N
for q in 1:N
for k in 1:N
os .+= U / N, "Cdagup", mod1(p - k, N), "Cdagdn", mod1(q + k, N), "Cdn", q, "Cup", p
end
end
end
```
`OpIDSum` plays the same roll as `OpSum` but in a much more efficient manner. To specify an operator in a term of an `OpSum` you specify a string (the operator's name) and a site index, whereas to specify an operator in a term of an `OpIDSum` you specify an `OpID` which contains an operator index and a site. The operator index is the index of the operator in the provided basis for the site.

The astute reader will notice that this `OpSum` has multiple operators acting on the same site. For example, it contains the term `"Cdagup", 1, "Cdagdn", 1, "Cdn", 1, "Cup", 1`. If we try and pass this `OpSum` to `MPO_New` directly it will throw an error. However, `"Cdagup", 1, "Cdagdn", 1, "Cdn", 1, "Cup", 1` is equivalent to `"Nup * Ndn", 1` and by passing a set of basis operators to `MPO_New` we can automatically convert any product of operators acting on a single site into one of these single basis operators. This is accomplished by the following
For $N = 200$ constructing an `OpIDSum` takes only 0.4s. Shown below is code to construct the Hamiltonian using an `OpIDSum`.

```julia
sites = siteinds("Electron", N; conserve_qns=true)

operatorNames = [
"I",
"Cdn",
"Cup",
"Cdagdn",
"Cdagup",
"Ndn",
"Nup",
"Cup * Cdn",
"Cup * Cdagdn",
"Cup * Ndn",
"Cdagup * Cdn",
"Cdagup * Cdagdn",
"Cdagup * Ndn",
"Nup * Cdn",
"Nup * Cdagdn",
"Nup * Ndn",
]

opCacheVec = [
[OpInfo(ITensors.Op(name, n), sites[n]) for name in operatorNames] for
n in eachindex(sites)
]

return MPO_new(os, sites; basisOpCacheVec=opCacheVec)
```
https://github.com/ITensor/ITensorMPOConstruction.jl/blob/main/examples/fermi-hubbard.jl#L79-L130

## Benchmarks: Fermi-Hubbard Hamiltonian in Momentum Space

Below is a plot of the bond dimension of the MPO produced by ITensors' default algorithm, [Renormalizer](https://github.com/shuaigroup/Renormalizer) which uses the [bipartite-graph algorithm](https://doi.org/10.1063/5.0018149), and ITensorMPOConstruction.

![](./docs/plot-generators/fh.png)

Of note is that the bond dimension of the MPO produced by Renormalizer scales as $O(N^2)$, both ITensors and ITensorMPOConstruction however produce an MPO with a bond dimension that scales as $O(N)$.

Below is a table of the time it took to construct the MPO (including the time it took to specify the Hamiltonian) for various number of sites. For ITensorMPOConstruction an `OpIDSum` was used. Some warm up was done for the Julia calculations to avoid measuring compilation overhead.

| $N$ | ITensors | Renormalizer | ITensorMPOConstruction |
|-----|----------|--------------|------------------------|
| 10 | 0.35s | 0.26 | 0.04s |
| 20 | 27s | 3.4s | 0.10s |
| 30 | N/A | 17s | 0.24s |
| 40 | N/A | 59s | 0.59s |
| 50 | N/A | 244s | 1.2s |
| 100 | N/A | N/A | 16s |
| 200 | N/A | N/A | 283s |

## Examples: Electronic Structure Hamiltonian
## Benchmarks: Electronic Structure Hamiltonian

After looking at the previous example one might assume that the relative speed of ITensorMPOConstruction over Renormalizer might be due to the fact that for the Fermi-Hubbard Hamiltonian ITensorMPOConstruction is able to construct a much more compact MPO. In the case of the electronic structure Hamiltonian all algorithms produce MPOs of similar bond dimensions.
After looking at the previous example you might assume that the relative speed of ITensorMPOConstruction over Renormalizer might be due to the fact that for the Fermi-Hubbard Hamiltonian ITensorMPOConstruction is able to construct a more compact MPO. In the case of the electronic structure Hamiltonian all algorithms produce MPOs of similar bond dimensions.

![](./docs/plot-generators/es.png)

However ITensorMPOConstruction is still an order of magnitude faster than Renormalizer. The code for this example can be found in [examples/electronic-structure.jl](https://github.com/ITensor/ITensorMPOConstruction.jl/blob/main/examples/electronic-structure.jl). The run time to generate these MPOs are shown below (again on a 2021 MacBook Pro with the M1 Max CPU and 32GB of memory).
However, ITensorMPOConstruction is still an order of magnitude faster than Renormalizer. The code for this example can be found in [examples/electronic-structure.jl](https://github.com/ITensor/ITensorMPOConstruction.jl/blob/main/examples/electronic-structure.jl). The run time to generate these MPOs (including the time it took to specify the Hamiltonian) are shown below. For ITensorMPOConstruction an `OpIDSum` was used.

| $N$ | ITensors | Renormalizer | ITensorMPOConstruction |
|-----|----------|--------------|------------------------|
| 10 | 3.0s | 2.1s | 0.37s |
| 20 | 498s | 61s | 7.1s |
| 30 | N/A | 458s | 46s |
| 40 | N/A | 2250s | 183 |
| 50 | N/A | N/A | 510s |
| 10 | 3.0s | 2.1s | 0.31s |
| 20 | 498s | 61s | 5.9s |
| 30 | N/A | 458s | 36s |
| 40 | N/A | 2250s | 162s |
| 50 | N/A | N/A | 504s |
| 60 | N/A | N/A | 1323s |
152 changes: 119 additions & 33 deletions examples/electronic-structure.jl
View file
Original file line number Diff line number Diff line change
@@ -1,41 +1,43 @@
using ITensorMPOConstruction
using ITensors

function electronic_structure_example(
function electronic_structure(
N::Int, complexBasisFunctions::Bool; useITensorsAlg::Bool=false
)::MPO
coeff() = rand() + 1im * complexBasisFunctions * rand()
coeff() = !complexBasisFunctions ? rand() : rand() + 1im * rand()

os = complexBasisFunctions ? OpSum{ComplexF64}() : OpSum{Float64}()
for a in 1:N
for b in a:N
epsilon_ab = coeff()
os .+= epsilon_ab, "Cdagup", a, "Cup", b
os .+= epsilon_ab, "Cdagdn", a, "Cdn", b

a == b && continue
os .+= conj(epsilon_ab), "Cdagup", b, "Cup", a
os .+= conj(epsilon_ab), "Cdagdn", b, "Cdn", a
@time "\tConstructing OpSum" let
for a in 1:N
for b in a:N
epsilon_ab = coeff()
os .+= epsilon_ab, "Cdagup", a, "Cup", b
os .+= epsilon_ab, "Cdagdn", a, "Cdn", b

a == b && continue
os .+= conj(epsilon_ab), "Cdagup", b, "Cup", a
os .+= conj(epsilon_ab), "Cdagdn", b, "Cdn", a
end
end
end

for j in 1:N
for s_j in ("dn", "up")
for k in 1:N
s_k = s_j
(s_k, k) <= (s_j, j) && continue
for j in 1:N
for s_j in ("dn", "up")
for k in 1:N
s_k = s_j
(s_k, k) <= (s_j, j) && continue

for l in 1:N
for s_l in ("dn", "up")
(s_l, l) <= (s_j, j) && continue
for l in 1:N
for s_l in ("dn", "up")
(s_l, l) <= (s_j, j) && continue

for m in 1:N
s_m = s_l
(s_m, m) <= (s_k, k) && continue
for m in 1:N
s_m = s_l
(s_m, m) <= (s_k, k) && continue

value = coeff()
os .+= value, "Cdag$s_j", j, "Cdag$s_l", l, "C$s_m", m, "C$s_k", k
os .+= conj(value), "Cdag$s_k", k, "Cdag$s_m", m, "C$s_l", l, "C$s_j", j
value = coeff()
os .+= value, "Cdag$s_j", j, "Cdag$s_l", l, "C$s_m", m, "C$s_k", k
os .+= conj(value), "Cdag$s_k", k, "Cdag$s_m", m, "C$s_l", l, "C$s_j", j
end
end
end
end
Expand All @@ -47,7 +49,7 @@ function electronic_structure_example(

## The only additional step required is to provide an operator basis in which to represent the OpSum.
if useITensorsAlg
return MPO(os, sites)
return @time "\tConstrucing MPO" MPO(os, sites)
else
operatorNames = [
"I",
Expand All @@ -73,18 +75,102 @@ function electronic_structure_example(
n in eachindex(sites)
]

return MPO_new(os, sites; basisOpCacheVec=opCacheVec)
return @time "\tConstrucing MPO" MPO_new(os, sites; basisOpCacheVec=opCacheVec)
end
end

let
N = 25
function electronic_structure_OpIDSum(N::Int, complexBasisFunctions::Bool)::MPO
operatorNames = [
[
"I",
"Cdn",
"Cup",
"Cdagdn",
"Cdagup",
"Ndn",
"Nup",
"Cup * Cdn",
"Cup * Cdagdn",
"Cup * Ndn",
"Cdagup * Cdn",
"Cdagup * Cdagdn",
"Cdagup * Ndn",
"Nup * Cdn",
"Nup * Cdagdn",
"Nup * Ndn",
] for _ in 1:N
]

↓ = false
↑ = true

opC(k::Int, spin::Bool) = OpID(2 + spin, k)
opCdag(k::Int, spin::Bool) = OpID(4 + spin, k)
opN(k::Int, spin::Bool) = OpID(6 + spin, k)

coeff() = !complexBasisFunctions ? rand() : rand() + 1im * rand()

# Ensure compilation
electronic_structure_example(5, false)
os = complexBasisFunctions ? OpIDSum{ComplexF64}() : OpIDSum{Float64}()
@time "\tConstructing OpIDSum" let
for a in 1:N
for b in a:N
epsilon_ab = coeff()
push!(os, epsilon_ab, opCdag(a, ↑), opC(b, ↑))
push!(os, epsilon_ab, opCdag(a, ↓), opC(b, ↓))

a == b && continue
push!(os, conj(epsilon_ab), opCdag(b, ↑), opC(a, ↑))
push!(os, conj(epsilon_ab), opCdag(b, ↓), opC(a, ↓))
end
end

for j in 1:N
for s_j in (↓, ↑)
for k in 1:N
s_k = s_j
(s_k, k) <= (s_j, j) && continue

for l in 1:N
for s_l in (↓, ↑)
(s_l, l) <= (s_j, j) && continue

for m in 1:N
s_m = s_l
(s_m, m) <= (s_k, k) && continue

value = coeff()
push!(os, value, opCdag(j, s_j), opCdag(l, s_l), opC(m, s_m), opC(k, s_k))
push!(
os, conj(value), opCdag(k, s_k), opCdag(m, s_m), opC(l, s_l), opC(j, s_j)
)
end
end
end
end
end
end
end

sites = siteinds("Electron", N; conserve_qns=true)
return @time "\tConstructing MPO" MPO_new(
os, sites, operatorNames; basisOpCacheVec=operatorNames
)
end

let
N = 20
useITensorsAlg = false

println("Constructing the Electronic Structure MPO for $N orbitals")
@time mpo = electronic_structure_example(N, false)
@time "Total" mpo = electronic_structure(N, false; useITensorsAlg=useITensorsAlg)
println("The maximum bond dimension is $(maxlinkdim(mpo))")
end

let
N = 20

println("Constructing the Electronic Structure MPO for $N orbitals using OpIDSum")
@time "Total" mpo = electronic_structure_OpIDSum(N, false)
println("The maximum bond dimension is $(maxlinkdim(mpo))")
end

Expand Down
Loading

0 comments on commit 637dd24

Please sign in to comment.