Skip to content

Commit 595fa72

Browse files
pbrehmermfherbst
andauthored
Progress (#34)
Co-authored-by: Michael F. Herbst <info@michael-herbst.com>
1 parent aac8667 commit 595fa72

31 files changed

+1085
-499
lines changed

Project.toml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
name = "ReducedBasis"
22
uuid = "f18afe42-61d0-402c-a0bd-58c747008192"
33
authors = ["Michael F. Herbst <info@michael-herbst.com>", "Paul Brehmer <paul.brehmer@rwth-aachen.de>"]
4-
version = "0.1.0"
4+
version = "0.1.1"
55

66
[deps]
77
DFTK = "acf6eb54-70d9-11e9-0013-234b7a5f5337"

README.md

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -10,7 +10,7 @@ accelerate the solution of a parametrized eigenvalue problems across the
1010
parameter domain.
1111

1212
In the RBM approach, a surrogate model is assembled by projecting the full
13-
problem onto a basis consisting of only a few tens of parameter snapshots, the
13+
problem onto a basis consisting of only a few tens of parameter snapshots.
1414
The package focuses on a greedy strategy that selects snapshots by maximally
1515
reducing the estimated error with each additional snapshot.
1616
Once the RBM surrogate is assembled, observables or post-processing can proceed

docs/make.jl

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -20,6 +20,8 @@ makedocs(;
2020
"Examples" => [
2121
"examples/xxz_ed.md",
2222
"examples/xxz_dmrg.md",
23+
"examples/xxz_pod.md",
24+
"examples/multi_ad.md",
2325
],
2426
"api.md"
2527
],

docs/src/examples/multi_ad.md

Lines changed: 121 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,121 @@
1+
# Affine decompositions with multi-indices and additional parameters
2+
3+
In this example, we want to explore the capabilities of the central [`AffineDecomposition`](@ref) type.
4+
To advance the previous examples where we covered the magnetization — a very simple observable that consists of only one affine term and a parameter-independent coefficient — we now turn to observables where the indices are multi-indices ``r = (r_1, \dots, r_d)`` and the coefficients can depend on additional parameters ``p``, aside from the ``\bm{\mu}`` parameter points that are present in the Hamiltonian:
5+
6+
```math
7+
O(\bm{\mu}, p) = \sum_{q=1}^Q \alpha_q(\bm{\mu}, p)\, O_q
8+
```
9+
10+
To stay within the realm of spin physics, we will consider the so-called *spin structure factor*
11+
12+
```math
13+
\mathcal{S}(k) = \frac{1}{L} \sum_{r,r'=1}^L e^{-i (r - r') k} S^z_r S^z_{r'},\quad
14+
\alpha_{r,r'}(k) = \frac{e^{-i (r - r') k}}{L}, \quad
15+
O_{r,r'} = S^z_r S^z_{r'}
16+
```
17+
18+
with a wavevector parameter ``k``, to discuss the implementation of a more complicated observable.
19+
20+
```@example multi_ad; continued = true
21+
using LinearAlgebra, SparseArrays, Plots, ReducedBasis # hide
22+
23+
σx = sparse([0.0 1.0; 1.0 0.0]) # hide
24+
σy = sparse([0.0 -im; im 0.0]) # hide
25+
σz = sparse([1.0 0.0; 0.0 -1.0]) # hide
26+
27+
function to_global(op::M, L::Int, i::Int) where {M<:AbstractMatrix} # hide
28+
d = size(op, 1) # hide
29+
30+
if i == 1 # hide
31+
return kron(op, M(I, d^(L - 1), d^(L - 1))) # hide
32+
elseif i == L # hide
33+
return kron(M(I, d^(L - 1), d^(L - 1)), op) # hide
34+
else # hide
35+
return kron(kron(M(I, d^(i - 1), d^(i - 1)), op), M(I, d^(L - i), d^(L - i))) # hide
36+
end # hide
37+
end # hide
38+
39+
function xxz_chain(L) # hide
40+
H1 = 0.25 * sum([to_global(σx, L, i) * to_global(σx, L, i + 1) + # hide
41+
to_global(σy, L, i) * to_global(σy, L, i + 1) for i in 1:(L-1)]) # hide
42+
H2 = 0.25 * sum([to_global(σz, L, i) * to_global(σz, L, i + 1) for i in 1:(L-1)]) # hide
43+
H3 = 0.5 * sum([to_global(σz, L, i) for i in 1:L]) # hide
44+
coefficient_map = μ -> [1.0, μ[1], -μ[2]] # hide
45+
AffineDecomposition([H1, H2, H3], coefficient_map) # hide
46+
end # hide
47+
48+
L = 6 # hide
49+
H = xxz_chain(L) # hide
50+
Δ = range(-1.0, 2.5; length=40) # hide
51+
hJ = range(0.0, 3.5; length=40) # hide
52+
grid_train = RegularGrid(Δ, hJ) # hide
53+
greedy = Greedy(; estimator=Residual(), tol=1e-3, init_from_rb=true) # hide
54+
lobpcg = LOBPCG(; n_target=1, tol_degeneracy=1e-4, tol=1e-9) # hide
55+
qrcomp = QRCompress(; tol=1e-9) # hide
56+
57+
info = assemble(H, grid_train, greedy, lobpcg, qrcomp) # hide
58+
basis = info.basis; h = info.h_cache.h; # hide
59+
60+
Δ_online = range(first(Δ), last(Δ); length=100) # hide
61+
hJ_online = range(first(hJ), last(hJ); length=100) # hide
62+
grid_online = RegularGrid(Δ_online, hJ_online) # hide
63+
fulldiag = FullDiagonalization(lobpcg) # hide
64+
```
65+
66+
So let us continue the first example where we have generated an ``L=6`` XXZ surrogate `basis` with a reduced Hamiltonian `h` using an exact diagonalization solver.
67+
Now the task is to implement the double-sum in ``\mathcal{S}``, as well as the ``k``-dependency in the coefficients.
68+
The double-sum can be encoded by putting all ``S^z_r S^z_{r'}`` combinations into a ``L \times L`` matrix:
69+
70+
``` @example multi_ad; continued = true
71+
terms = map(idx -> to_global(σz, L, first(idx.I)) * to_global(σz, L, last(idx.I)),
72+
CartesianIndices((1:L, 1:L)))
73+
```
74+
75+
Correspondingly, the coefficient function now has to map one ``k`` value to a matrix of coefficients of the same size as the `terms` matrix:
76+
77+
``` @example multi_ad; continued = true
78+
coefficient_map = k -> map(idx -> cis(-(first(idx.I) - last(idx.I)) * k) / L,
79+
CartesianIndices((1:L, 1:L)))
80+
```
81+
82+
One feature of the structure factor that also shows up in many other affine decompositions with double-sums is that the term indices commute, i.e. ``O_{r,r'} = O_{r',r}``.
83+
In that case, only the upper triangular matrix has to be computed since ``B^\dagger O_{r,r'} B = B^\dagger O_{r',r} B`` are the same in the compressed affine decomposition.
84+
So let's create the [`AffineDecomposition`](@ref) and compress, exploiting this symmetry using the `symmetric_terms` keyword argument:
85+
86+
``` @example multi_ad; continued = true
87+
SFspin = AffineDecomposition(terms, coefficient_map)
88+
sfspin, _ = compress(SFspin, basis; symmetric_terms=true)
89+
```
90+
91+
In the online evaluation of the structure factor, we then need to define some wavevector values and compute the structure factor at each of them.
92+
With the `grid_online` from before, this reads:
93+
94+
``` @example multi_ad; continued = true
95+
wavevectors = [0.0, π/4, π/2, π]
96+
sf = [zeros(size(grid_online)) for _ in 1:length(wavevectors)]
97+
for (idx, μ) in pairs(grid_online)
98+
_, φ_rb = solve(h, basis.metric, μ, fulldiag)
99+
for (i, k) in enumerate(wavevectors)
100+
sf[i][idx] = sum(eachcol(φ_rb)) do u
101+
real(dot(u, sfspin(k), u))
102+
end / size(φ_rb, 2)
103+
end
104+
end
105+
```
106+
107+
Here we again see the convenience of measuring observables in the online stage;
108+
adding more wavevector values does not significantly increase the computational cost, since it corresponds to a mere reevaluation of the coefficient functions and small vector-matrix products.
109+
Finally, let us see how the structure factor behaves for the different wavevector values:
110+
111+
``` @example multi_ad
112+
kwargs = (; xlabel=raw"$\Delta$", ylabel=raw"$h/J$", colorbar=true, leg=false) # hide
113+
hms = []
114+
for (i, q) in enumerate(wavevectors)
115+
push!(hms, heatmap(grid_online.ranges[1], grid_online.ranges[2], sf[i]';
116+
title="\$k = $(round(q/π; digits=3))\\pi\$", kwargs...))
117+
end
118+
plot(hms...)
119+
```
120+
121+
It can be nicely seen that the spin structure factor indicates the ferromagnetic phase at ``k=0`` and then moves through the magnetization plateaus until it reaches the antiferromagnetic plateau at ``k=\pi``.

docs/src/examples/xxz_dmrg.md

Lines changed: 25 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
# Greedy basis assembly using DMRG
22

3-
As a follow-up example, we now want to showcase how to compute an `RBasis` by means of the [density matrix renormalization group](https://tensornetwork.org/mps/algorithms/dmrg/) (DMRG).
4-
To that end, we utilize the `ITensors.jl` library which, among other things, efficiently implements DMRG.
3+
As a follow-up example, we now want to showcase how to compute a reduced basis by means of the [density matrix renormalization group](https://tensornetwork.org/mps/algorithms/dmrg/) (DMRG).
4+
To that end, we utilize the ITensors.jl package which, among other things, efficiently implements DMRG.
55
We will see that, while we need to adjust the way we set up the model Hamiltonian as well as our solver, most steps stay the same.
66
Again, we treat the one-dimensional ``S=1/2`` XXZ model from the previous example.
77

@@ -19,7 +19,7 @@ using ReducedBasis
1919
```
2020

2121
now featuring `ITensors`.
22-
To build the Hamiltonian terms as MPOs, we make use of the [`OpSum()`](https://itensor.github.io/ITensors.jl/stable/tutorials/DMRG.html) object that automatically produces an MPO from a string of operators.
22+
To build the Hamiltonian terms as MPOs, we make use of the `ITensors.OpSum()` object that automatically produces an MPO from a string of operators.
2323
The affine MPO terms are then stored in an [`AffineDecomposition`](@ref) as [`ApproxMPO`](@ref)s which also include possible truncation keyword arguments:
2424

2525
```@example xxz_dmrg; continued = true
@@ -35,29 +35,27 @@ function xxz_chain(sites::IndexSet; kwargs...)
3535
end
3636
magn_term += "Sz", length(sites) # Add last magnetization term
3737
coefficient_map = μ -> [1.0, μ[1], -μ[2]]
38-
AffineDecomposition(
39-
[ApproxMPO(MPO(xy_term, sites), xy_term; kwargs...),
40-
ApproxMPO(MPO(zz_term, sites), zz_term; kwargs...),
41-
ApproxMPO(MPO(magn_term, sites), magn_term; kwargs...)],
42-
coefficient_map,
43-
)
38+
AffineDecomposition([ApproxMPO(MPO(xy_term, sites), xy_term; kwargs...),
39+
ApproxMPO(MPO(zz_term, sites), zz_term; kwargs...),
40+
ApproxMPO(MPO(magn_term, sites), magn_term; kwargs...)],
41+
coefficient_map)
4442
end
4543
```
4644

4745
So let us instantiate such an MPO Hamiltonian where we also specify a singular value `cutoff`, which is passed to the [`ApproxMPO`](@ref) objects:
4846

4947
```@example xxz_dmrg; continued = true
50-
L = 12
48+
L = 18
5149
sites = siteinds("S=1/2", L)
5250
H = xxz_chain(sites; cutoff=1e-9)
5351
```
5452

55-
We now chose a bigger system size, since the tensor format allows for efficient low rank approximations (hence the `cutoff`) that buy us a substantial performance advantage when going to larger systems.
53+
Notice that we can now choose a bigger system size (which is still very small here), since the tensor format allows for efficient low rank approximations (hence the `cutoff`) that buy us a substantial performance advantage when going to larger systems.
5654

5755
## Using the [`DMRG`](@ref) solver for basis assembly
5856

5957
Having created our Hamiltonian in MPO format, we now need a solver that is able to compute ground states from MPOs.
60-
The corresponding ground state will also be obtained in a tensor format, namely as *matrix product states* (MPS).
58+
The corresponding ground state will also be obtained in a tensor format, namely as a *matrix product state* (MPS).
6159
This is achieved by `ITensors.dmrg` which is wrapped in the [`DMRG`](@ref) solver type:
6260

6361
```@example xxz_dmrg; continued = true
@@ -66,10 +64,11 @@ dm = DMRG(; n_states=1, tol_degeneracy=0.0,
6664
observer=() -> DMRGObserver(; energy_tol=1e-9))
6765
```
6866

69-
While the implemented DMRG solver is capable of solving degenerate ground state, we here opt for non-degenerate settings (i.e. `n_states=1` and `tol_degeneracy=0.0`), since one encounters a ``L+1``-fold degeneracy on the parameter domain, where the degenerate DMRG solver can produce instable results for larger ``L``.
67+
While the implemented DMRG solver is capable of solving degenerate ground states, we here opt for non-degenerate settings, i.e. `n_states=1` and `tol_degeneracy=0.0`.
68+
(We do this due to a ``L+1``-fold degeneracy on the parameter domain, where the degenerate DMRG solver can produce instable results for larger ``L``.)
7069

7170
As discussed in the last example, we need a way to orthogonalize the reduced basis.
72-
Due to the MPS format that the snapshots will have, we cannot use QR decompositions anymore and resort to a different method, [`EigenDecomposition`](@ref), featuring an eigenvalue decomposition of the snapshot overlap matrix:
71+
Due to the MPS format that the snapshots will have, we cannot use QR decompositions anymore and resort to a different method, [`EigenDecomposition`](@ref), featuring an eigenvalue decomposition of the snapshot overlap matrix that can drop insignificant snapshots that fall below a cutoff:
7372

7473
```@example xxz_dmrg; continued = true
7574
edcomp = EigenDecomposition(; cutoff=1e-7)
@@ -79,22 +78,23 @@ edcomp = EigenDecomposition(; cutoff=1e-7)
7978
Δ = range(-1.0, 2.5; length=40) # hide
8079
hJ = range(0.0, 3.5; length=40) # hide
8180
grid_train = RegularGrid(Δ, hJ) # hide
82-
greedy = Greedy(; estimator=Residual(), n_truth_max=22, init_from_rb=true) # hide
81+
greedy = Greedy(; estimator=Residual(), n_truth_max=32, init_from_rb=true) # hide
8382
```
8483

85-
Now with different types for `H`, the solver and the orthogonalizer, we call `assemble` using the `greedy` strategy and training grid from the last example:
84+
Now with different types for the Hamiltonian, the solver and the orthogonalizer, we call `assemble` using the `greedy` strategy and training grid from the last example:
8685

8786
```@example xxz_dmrg; continued = true
88-
basis, h, info = assemble(H, grid_train, greedy, dm, edcomp)
87+
info = assemble(H, grid_train, greedy, dm, edcomp)
88+
basis = info.basis; h = info.h_cache.h;
8989
```
9090

9191
The returned `basis` now has snapshot vectors of `ITensors.MPS` type, which we have to keep in mind when we want to compress observables.
92-
That is to say, the observables have to be constructed as `AffineDecompositions` with `ApproxMPO` terms as for the Hamiltonian.
92+
That is to say, the observables have to be constructed as [`AffineDecomposition`](@ref)s with [`ApproxMPO`](@ref) terms as we did for the Hamiltonian.
9393
Again, we want to compute the magnetization so that we can reuse the third term of `H`:
9494

9595
```@example xxz_dmrg; continued = true
9696
M = AffineDecomposition([H.terms[3]], μ -> [2 / L])
97-
m = compress(M, basis)
97+
m, _ = compress(M, basis)
9898
```
9999

100100
```@example xxz_dmrg; continued = true
@@ -112,8 +112,8 @@ magnetization = map(grid_online) do μ # hide
112112
end # hide
113113
```
114114

115-
And at that point, we can continue as before since we have arrived at the online phase where we only operate in the low-dimensional reduced basis space, agnostic of the previously used solver method.
116-
We have to make sure, however, to choose matching degeneracy settings for the [`FullDiagonalization`](@ref) solver in the online phase, which we do via the matching constructor:
115+
And at that point, we continue as before since we have arrived at the online phase where we only operate in the low-dimensional reduced basis space, agnostic of the snapshot solver method.
116+
We have to make sure, however, to choose matching degeneracy settings for the [`FullDiagonalization`](@ref) solver in the online phase:
117117

118118
```@example xxz_dmrg; continued = true
119119
fulldiag = FullDiagonalization(dm)
@@ -128,18 +128,18 @@ grid_online = RegularGrid(Δ_online, hJ_online) # hide
128128
magnetization = map(grid_online) do μ # hide
129129
_, φ_rb = solve(h, basis.metric, μ, fulldiag) # hide
130130
sum(eachcol(φ_rb)) do u # hide
131-
abs(dot(u, m_reduced, u)) / size(φ_rb, 2) # hide
132-
end # hide
131+
abs(dot(u, m_reduced, u)) # hide
132+
end / size(φ_rb, 2) # hide
133133
end # hide
134134
```
135135

136-
In the same way as before, we perform the online calculations and arrive at the following magnetization plot:
136+
In the same way as before, we perform the online computations and arrive at the following magnetization plot:
137137

138138
```@example xxz_dmrg
139139
hm = heatmap(grid_online.ranges[1], grid_online.ranges[2], magnetization';
140140
xlabel=raw"$\Delta$", ylabel=raw"$h/J$", title="magnetization ",
141141
colorbar=true, clims=(0.0, 1.0), leg=false)
142-
plot!(hm, grid_online.ranges[1], x -> 1 + x; lw=2, ls=:dash, legend=false, color=:fuchsia)
142+
plot!(hm, grid_online.ranges[1], x -> 1 + x; lw=2, ls=:dash, legend=false, color=:green)
143143
params = unique(basis.parameters)
144144
scatter!(hm, [μ[1] for μ in params], [μ[2] for μ in params];
145145
markershape=:xcross, color=:springgreen, ms=3.0, msw=2.0)

0 commit comments

Comments
 (0)