Skip to content

Commit 72737da

Browse files
authored
Merge pull request #30 from ITensor/docs-update
Updated docs with block2 comparison.
2 parents 1a566b1 + d8915ca commit 72737da

4 files changed

Lines changed: 122 additions & 103 deletions

File tree

README.md

Lines changed: 65 additions & 83 deletions
Original file line numberDiff line numberDiff line change
@@ -127,7 +127,7 @@ Additionally, a further constructor is provided which takes in a modifying `func
127127

128128
## Haldane-Shasty Hamiltonian and truncation
129129

130-
The Haldane-Shasty Hamiltonian defined on $N$ spin-$\frac{1}{2}$ particles is
130+
The Haldane-Shasty Hamiltonian defined on $N$ spin-half particles is
131131

132132
$$
133133
H = \frac{J \pi^2}{N^2} \sum_{n = 1}^N \sum_{m = n + 1}^N \frac{\mathbf{S}_m \cdot \mathbf{S}_n}{\sin^2 \left( \pi \frac{n - m}{N} \right)} \ .
@@ -143,100 +143,82 @@ Starting with the MPO from ITensorMPOConstruction obtained with the standard `to
143143

144144
## Benchmarks: Fermi-Hubbard Hamiltonian in Momentum Space
145145

146-
We constructed the momentum space Fermi-Hubbard Hamiltonian using ITensorMPS, ITensorMPOConstruction and [block2](https://github.com/block-hczhai/block2-preview) which has one of the most sophisticated MPO construction algorithms.
147-
148-
For block2, we used the `FastBlockedDisjointSVD` algorithm for MPO construction.
149-
150-
For even $N$, the Hamiltonian can be represented exactly as an MPO of bond dimension $10 N - 4$, and all the algorithms achieve this minimal bond dimension. ITensorMPOConstruction is also not only able to construct this particular MPO much faster than the competition, but the sparsity of the resulting MPO is much higher.
146+
We constructed the momentum space Fermi-Hubbard Hamiltonian using ITensorMPS and ITensorMPOConstruction. For even $N$, the Hamiltonian can be represented exactly as an MPO of bond dimension $10 N - 4$, and both algorithms achieve this minimal bond dimension. ITensorMPOConstruction is also not only able to construct this particular MPO much faster, but the sparsity of the resulting MPO is much higher.
151147

152148
### Bond Dimension
153-
| $N$ | ITensorMPS | block2 | ITensorMPOConstruction |
154-
|-----|------------|--------|------------------------|
155-
| 10 | 96 | 96 | 96 |
156-
| 20 | 196 | 196 | 196 |
157-
| 30 | 296 | 296 | 296 |
158-
| 40 | N/A | 396 | 396 |
159-
| 50 | N/A | 496 | 496 |
160-
| 100 | N/A | 996 | 996 |
161-
| 200 | N/A | 1996 | 1996 |
162-
| 300 | N/A | 2996 | 2996 |
163-
| 400 | N/A | 3996 | 3996 |
164-
| 500 | N/A | ???? | 4996 |
149+
| $N$ | ITensorMPS | ITensorMPOConstruction |
150+
|-----|------------|------------------------|
151+
| 10 | 96 | 96 |
152+
| 20 | 196 | 196 |
153+
| 30 | 296 | 296 |
154+
| 40 | N/A | 396 |
155+
| 50 | N/A | 496 |
156+
| 100 | N/A | 996 |
157+
| 200 | N/A | 1996 |
158+
| 300 | N/A | 2996 |
159+
| 400 | N/A | 3996 |
160+
| 500 | N/A | 4996 |
165161

166162
### Sparsity
167-
| $N$ | ITensorMPS | block2 | ITensorMPOConstruction |
168-
|-----|------------|--------|------------------------|
169-
| 10 | 92.7% | 96.7% | 99.32% |
170-
| 20 | 92.6% | 98.4% | 99.70% |
171-
| 30 | 92.6% | 99.0% | 99.81% |
172-
| 40 | N/A | 99.2% | 99.86% |
173-
| 50 | N/A | 99.4% | 99.89% |
174-
| 100 | N/A | 99.7% | 99.94% |
175-
| 200 | N/A | 99.85% | 99.97% |
176-
| 300 | N/A | 99.90% | 99.982% |
177-
| 400 | N/A | 99.92% | 99.986% |
178-
| 500 | N/A | N/A | 99.999% |
163+
164+
Sparsity of the `ITensorMPS` MPO with the default `splitblocks=true`, and the `ITensorMPOConstruction` MPO with the less aggressive `combine_qn_sectors::Bool=false`.
165+
166+
| $N$ | ITensorMPS | ITensorMPOConstruction |
167+
|-----|------------|------------------------|
168+
| 10 | 92.7% | 99.32% |
169+
| 20 | 92.6% | 99.70% |
170+
| 30 | 92.6% | 99.81% |
171+
| 40 | N/A | 99.86% |
172+
| 50 | N/A | 99.89% |
173+
| 100 | N/A | 99.94% |
174+
| 200 | N/A | 99.97% |
175+
| 300 | N/A | 99.982% |
176+
| 400 | N/A | 99.986% |
177+
| 500 | N/A | 99.999% |
179178

180179
### Runtime
181-
| $N$ | ITensorMPS | block2 | ITensorMPOConstruction |
182-
|-----|------------|--------|------------------------|
183-
| 10 | 0.32s | 0.016s | 0.009s |
184-
| 20 | 30.6s | 0.089s | 0.052s |
185-
| 30 | 792s | 0.30s | 0.14s |
186-
| 40 | N/A | 0.72s | 0.38s |
187-
| 50 | N/A | 1.5s | 0.63s |
188-
| 100 | N/A | 20s | 7.7s |
189-
| 200 | N/A | 489s | 103s |
190-
| 300 | N/A | 3711s | 500s |
191-
| 400 | N/A | 18373s | 1554s |
192-
| 500 | N/A | N/A | 3802s |
180+
| $N$ | ITensorMPS | ITensorMPOConstruction |
181+
|-----|------------|------------------------|
182+
| 10 | 0.32s | 0.009s |
183+
| 20 | 30.6s | 0.052s |
184+
| 30 | 792s | 0.14s |
185+
| 40 | N/A | 0.38s |
186+
| 50 | N/A | 0.63s |
187+
| 100 | N/A | 7.7s |
188+
| 200 | N/A | 103s |
189+
| 300 | N/A | 500s |
190+
| 400 | N/A | 1554s |
191+
| 500 | N/A | 3802s |
192+
193+
## A note on sparsity
194+
195+
The core component of pretty many MPO construction algorithms is to take an operator $\hat{O}$ defined on a bipartite system
193196

197+
$$
198+
\hat{O} = \sum_{a = 1}^{N_A} \sum_{b = 1}^{N_B} \gamma_{ab} \hat{A}_a \otimes \hat{B}_b \ ,
199+
$$
194200

195-
## Benchmarks: Electronic Structure Hamiltonian
201+
and turn it into a two site MPO
196202

197-
We also construct a particle number and spin preserving two-electron Hamiltonian with random coefficients, to mock the performance of our algorithm on constructing the electronic structure Hamiltonian. In this case we compare against the `FastBipartite` algorithm from block2. This bipartite algorithm essentially performs no compression of the resulting MPO, since it ignores the specific value of the coefficients, but it produces MPOs of high sparsity. In this case, since the coefficients are random, there is no underlying pattern to compress and so the bipartite algorithm works well.
203+
$$
204+
\hat{O} = \sum_{\chi = 1}^w \left( \sum_{a = 1}^{N_A} \alpha_{\chi a} \hat{A}_a \right) \otimes \left( \sum_{b = 1}^{N_B} \beta_{\chi b} \hat{B}_b \right) \ .
205+
$$
198206

199-
We see that for $N \geq 40$, where $N$ is the number of spin-orbitals, ITensorMPOConstruction constructs MPOs of bond dimension slightly larger than block2. This increase in bond dimension can be explained by the attempt to perform compression on an incompressible Hamiltonian. For $N \geq 70$, ITensorMPOConstruction is also slower than block2.
207+
The MPO bond dimension is $w$ and the MPO tensors are essentially operator valued vectors who's $\chi$-th entry is $\left( \sum_{a = 1}^{N_A} \alpha_{\chi a} \hat{A}_a \right)$ for the left tensor. This vector of operators can be reshaped into the standard matrix of operators if $a$ is a combined incoming link and onsite index.
200208

201-
By default, the MPO from ITensorMPOConstruction is also denser than the MPOs from ITensorMPS and block2. However, both ITensorMPS and block2 create blocks of size one, whereas ITensorMPOConstruction creates larger blocks. Using `ITensorMPS.splitblocks` we can split the larger blocks in the MPO from ITensorMPOConstruction up into blocks of size one. After this, the MPO from ITensorMPOConstruction is sparser than the competition.
209+
In the case of operators with a global $U(1)$ symmetry the matrix $\gamma$ can be permuted into a block diagonal form. This form has many benefits since each block can be decomposed into a MPO independently and the $\alpha$ and $\beta$ matrices inherit the block diagonal nature. The way that $\gamma$ is brought into block diagonal form in the original `ITensorMPS` algorithm is to use the quantum numbers associated with the operators $\hat{A}_a$ and $\hat{B}_b$. This has two problems; first it requires that the user provide the symmetry information, and second there may be other block diagonal blocks unrelated to the symmetries. The approach taken in this library addresses both these issues.
202210

203-
In some sense, this Hamiltonian is a poor match for ITensorMPOConstruction due to its relatively low sparsity and lack of compression. Nevertheless, ITensorMPOConstruction is competitive with the better suited bipartite algorithm. A more direct comparison would be with block2's `FastBlockedDisjointSVD` algorithm, which for $N = 70$ constucts an MPO of bond dimension 10117 that is 86% sparse in 1431 seconds. While the bond dimension is similar to the other methods, the sparsity and construction time are worse.
211+
In `ITensorMPOConstruction` $\gamma$ is brought into block diagonal form by viewing it as a bipartite graph adjacency matrix ($\gamma_{a b} \neq 0$ implies there is an edge between left-vertex $a$ and right-vertex $b$) and finding the connected components. Each connected component is then a block in the block diagonal representation. This does not require the use of any symmetry information and is guaranteed to produce the maximum possible number of blocks. However, although `ITensorMPOConstruction` will produce a matrix $\alpha$ of minimal bond dimension $w$ and greatest number of diagonal blocks this does not mean that the overall sparsity of the MPO tensor is maximized. This is because `ITensors`' sparse format is more flexible, and most of the time each diagonal block in the $\alpha$ matrix winds up being stored itself in a sparse format. We use the sparse QR decomposition to decompose each block of $\gamma$, and while the resulting matrices (which become the blocks in $\alpha$ and $\beta$) are sparse, their sparsity is not necessarily optimal.
204212

205-
### Bond Dimension
206-
| $N$ | ITensorMPS | block2 | ITensorMPOConstruction |
207-
|-----|------------|--------|------------------------|
208-
| 10 | 227 | 227 | 227 |
209-
| 20 | 852 | 852 | 852 |
210-
| 30 | N/A | 1877 | 1877 |
211-
| 40 | N/A | 3302 | 3355 |
212-
| 50 | N/A | 5127 | 5134 |
213-
| 60 | N/A | 7352 | 7364 |
214-
| 70 | N/A | 9977 | 9985 |
215-
| 80 | N/A | 13002 | 13006 |
216-
| 90 | N/A | 16427 | 16473 |
213+
To illustrate the suboptimal sparsity, we turn to [Block2](https://github.com/block-hczhai/block2-preview) which has a sophisticated set of MPO construction algorithms. Specifically we will use the `FastBipartite` algorithm, based on the bipartite graph algorithm from [RenLi2020](https://doi.org/10.1063/5.0018149). The bipartite algorithm is very efficient and also produces MPO tensors of exceptional sparsity. The drawback is that it is unable to compress the MPO bond dimension. For example, for the momentum space Fermi-Hubbard Hamiltonian the bond dimension it produces is $O(N^2)$. However, for some operators such as the electronic structure Hamiltonian the bipartite algorithm and the rank decomposition algorithm (used in `ITensorMPS` and here) produce similar MPO bond dimensions. In these cases, the bipartite algorithm will likely produce MPOs of greater sparsity.
217214

218-
### Sparsity
219-
| $N$ | ITensorMPS | block2 | ITensorMPOConstruction | After `splitblocks` |
220-
|-----|------------|--------|------------------------|---------------------|
221-
| 10 | 94.5% | 89.2% | 93.4% | 95.7% |
222-
| 20 | 95.4% | 93.1% | 94.1% | 97.1% |
223-
| 30 | N/A | 95.0% | 94.4% | 97.6% |
224-
| 40 | N/A | 96.0% | 94.4% | 97.8% |
225-
| 50 | N/A | 96.7% | 94.6% | 97.9% |
226-
| 60 | N/A | 97.2% | 94.6% | 97.9% |
227-
| 70 | N/A | 97.6% | 94.7% | N/A |
228-
| 80 | N/A | 97.8% | 94.7% | N/A |
229-
| 90 | N/A | 98.1% | 94.7% | N/A |
215+
In the table below we present data from constructing two different electronic structure Hamiltonians, the second of which is from [ZhaiLee2023](https://doi.org/10.1021/acs.jpca.3c06142). Our rank decomposition algorithm only slightly reduces the bond dimension compared to the bipartite MPO from `Block2`, but it results in a much denser MPO. This increase in sparsity has a significant impact on the subsequent DMRG performance, which is larger by 75% for our rank decomposition MPO (timings and sparsities taken from `Block2`'s by transferring over the MPO from `ITensorMPOConstruction`).
230216

231-
### Runtime
232-
| $N$ | ITensorMPS | block2 | ITensorMPOConstruction |
233-
|-----|------------|--------|------------------------|
234-
| 10 | 3.65s | 0.166s | 0.052s |
235-
| 20 | 605s | 2.67s | 1.21s |
236-
| 30 | N/A | 15.0s | 7.00s |
237-
| 40 | N/A | 50.4s | 29.5s |
238-
| 50 | N/A | 137s | 104s |
239-
| 60 | N/A | 332s | 284s |
240-
| 70 | N/A | 619s | 625s |
241-
| 80 | N/A | 1220s | 1545s |
242-
| 90 | N/A | 1944s | 3968s |
217+
| system | algorithm | bond dimension | sparsity |
218+
|-------------------------------------------------------------------------------------|-----------|----------------|----------|
219+
| C <sub>2</sub> | rank | 704 | 97.21% |
220+
| C <sub>2</sub> | bipartite | 722 | 98.68% |
221+
| [Fe <sub>2</sub> S(C H <sub>3</sub>) (S C H <sub>3</sub>)<sub>4</sub>]<sup>3-</sup> | rank | 2698 | 88.70% |
222+
| [Fe <sub>2</sub> S(C H <sub>3</sub>) (S C H <sub>3</sub>)<sub>4</sub>]<sup>3-</sup> | bipartite | 2738 | 95.64% |
223+
224+
To further complicate matters `Block2` has a different, less flexible, sparse storage format from `ITensors`. Specifically, they store MPO tensors in a "matrix of operators" format, where the onsite operator is always dense. Essentially this format stores a sparse representation of $\alpha$, while maintaining a dense form for $\hat{A}_a$. To facilitate comparisons between the two libraries without having to convert the MPOs between them, we provide the function `block2_nnz(mpo::MPO)::Tuple{Int, Int}` which returns the total number of blocks (the size of $\alpha$ summed across each site) and the number of non-zero blocks.

docs/plot-generators/block2-plots.py

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -115,10 +115,10 @@ def electronic_structure(N, alg):
115115
print(f"N = {N}, time = {stop - start}")
116116

117117

118-
for N in [10]:
119-
fermi_hubbard(N, MPOAlgorithmTypes.FastBlockedDisjointSVD)
118+
for N in [10, 26]:
119+
fermi_hubbard(N, MPOAlgorithmTypes.FastBlockedDisjointSVD, J=-0.5)
120120
print()
121121

122-
for N in [10, 10, 20, 30, 40, 50, 60, 70]:
123-
electronic_structure(N, MPOAlgorithmTypes.FastBipartite)
124-
print()
122+
# for N in [10, 10, 20, 30, 40, 50, 60, 70]:
123+
# electronic_structure(N, MPOAlgorithmTypes.FastBipartite)
124+
# print()

examples/fermi-hubbard.jl

Lines changed: 17 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -183,30 +183,32 @@ let
183183
end
184184
end
185185

186-
let
187-
for N in [10]
188-
println(
189-
"Constructing the Fermi-Hubbard momentum space MPO for $N sites using ITensorMPS"
190-
)
191-
@time "Total construction time" mpo = Fermi_Hubbard_momentum_space(
192-
N; useITensorsAlg=true
193-
)
194-
println("The maximum bond dimension is $(maxlinkdim(mpo))")
195-
println("The sparsity is $(ITensorMPOConstruction.sparsity(mpo))")
196-
println()
197-
end
198-
end
186+
# let
187+
# for N in [10]
188+
# println(
189+
# "Constructing the Fermi-Hubbard momentum space MPO for $N sites using ITensorMPS"
190+
# )
191+
# @time "Total construction time" mpo = Fermi_Hubbard_momentum_space(
192+
# N; useITensorsAlg=true
193+
# )
194+
# println("The maximum bond dimension is $(maxlinkdim(mpo))")
195+
# println("The sparsity is $(ITensorMPOConstruction.sparsity(mpo))")
196+
# @show ITensorMPOConstruction.block2_nnz(mpo)
197+
# println()
198+
# end
199+
# end
199200

200201
let
201-
for N in [10]
202+
for N in [10, 26]
202203
println(
203204
"Constructing the Fermi-Hubbard momentum space MPO for $N sites using ITensorMPOConstruction",
204205
)
205206
@time "Total construction time" mpo = transcorrelated_Fermi_Hubbard_momentum_space_OpIDSum(
206-
N
207+
N, 1, 4, -0.5
207208
)
208209
println("The maximum bond dimension is $(maxlinkdim(mpo))")
209210
println("The sparsity is $(ITensorMPOConstruction.sparsity(mpo))")
211+
@show ITensorMPOConstruction.block2_nnz(mpo)
210212
println()
211213
end
212214
end

src/MPOConstruction.jl

Lines changed: 35 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -202,3 +202,38 @@ function sparsity(mpo::MPO)::Float64
202202

203203
return num_zeros / num_entries
204204
end
205+
206+
function block2_nnz(mpo::MPO)::Tuple{Int, Int}
207+
total_blocks = 0
208+
nnz_blocks = 0
209+
sites = noprime(siteinds(first, mpo))
210+
for (i, t) in enumerate(mpo)
211+
t = t.tensor
212+
total_blocks += prod(size(t)) ÷ dim(sites[i])^2
213+
214+
link_locs = [j for j in 1:ndims(t) if inds(t)[j] (dag(sites[i]), sites[i]')]
215+
@assert length(link_locs) (1, 2)
216+
217+
nz_link_coords = Set{Tuple{Int, Int}}()
218+
219+
for b in ITensors.eachnzblock(t)
220+
block = ITensors.blockview(t, b)
221+
blockStart = NDTensors.blockstart(t, b) .- 1
222+
223+
for coords in CartesianIndices(block)
224+
value = block[coords]
225+
value == 0 && continue
226+
227+
coords = coords.I .+ blockStart
228+
229+
link_coords = coords[link_locs]
230+
length(link_coords) == 1 && (link_coords = (only(link_coords), 0))
231+
push!(nz_link_coords, link_coords)
232+
end
233+
end
234+
235+
nnz_blocks += length(nz_link_coords)
236+
end
237+
238+
return total_blocks, nnz_blocks
239+
end

0 commit comments

Comments
 (0)