Basis assembly using Proper Orthogonal Decomposition

In Greedy basis assembly using DMRG we have seen that we can customize the snapshot solvers as well as the compression methods during RB assembly. What we want to demonstrate in this example is that we can also use different strategies for basis assembly altogether. In particular, we will show how to use the Proper Orthogonal Decomposition (POD) technique in the offline stage.

First we provide the setup already discussed in The reduced basis workflow:

using LinearAlgebra
using SparseArrays
using ReducedBasis

σx = sparse([0.0 1.0; 1.0 0.0])
σy = sparse([0.0 -im; im 0.0])
σz = sparse([1.0 0.0; 0.0 -1.0])

function to_global(op::M, L::Int, i::Int) where {M<:AbstractMatrix}
    d = size(op, 1)
    if i == 1
        kron(op, M(I, d^(L - 1), d^(L - 1)))
    elseif i == L
        kron(M(I, d^(L - 1), d^(L - 1)), op)
    else
        kron(kron(M(I, d^(i - 1), d^(i - 1)), op), M(I, d^(L - i), d^(L - i)))
    end
end

function xxz_chain(L)
    H1 = 0.25 * sum(1:L-1) do i
        to_global(σx, L, i) * to_global(σx, L, i + 1) +
        to_global(σy, L, i) * to_global(σy, L, i + 1)
    end
    H2 = 0.25 * sum(1:L-1) do i
        to_global(σz, L, i) * to_global(σz, L, i + 1)
    end
    H3 = 0.5  * sum(1:L) do i
        to_global(σz, L, i)
    end
    AffineDecomposition([H1, H2, H3], μ -> [1.0, μ[1], -μ[2]])
end;

The conceptual difference between POD and the greedy assembly strategy is that with POD, a truth solve is performed at all parameter points in the selected training grid, followed by a singular value decomposition of the snapshot matrix. In this way, we obtain an orthogonal basis by using the singular vectors as our RB. While this procedure is less complex than the greedy strategy, it comes with the significantly increased cost of having to solve snapshots at all grid points. Nonetheless, it can be useful, e.g., to obtain a reference RB and to compare against a greedy basis.

So let us stay with the example of the XXZ spin chain and initialize the Hamiltonian as before (using the functions defined in the first example) and choose a grid as well as a solver method:

L = 6
H = xxz_chain(L)

Δ = range(-1.0, 2.5; length=20)
hJ = range(0.0, 3.5; length=20)
grid_train = RegularGrid(Δ, hJ)

lobpcg = LOBPCG(; tol_degeneracy=1e-4);

Notice that we now use a coarser 20 × 20 grid since we perform truth solves on all parameter points and want to keep the computational effort low. Moreover, we are restricted to exact diagonalization solvers, since we need to explicitly construct the snapshot matrix in order to be able to perform an SVD on it. To assemble using POD, we create a POD object where we specify the number of retained columns, i.e., singular vectors of the snapshot matrix:

pod = POD(; n_vectors=24);

We then call assemble using our parameters, including pod, which selects the POD assembly method:

rbres = assemble(H, grid_train, pod, lobpcg);

Truth solving on (20, 20) grid...   1%|▏                 |  ETA: 0:00:11
Truth solving on (20, 20) grid...   2%|▌                 |  ETA: 0:00:08
Truth solving on (20, 20) grid...   4%|▊                 |  ETA: 0:00:07
Truth solving on (20, 20) grid...   6%|█▏                |  ETA: 0:00:07
Truth solving on (20, 20) grid...   8%|█▌                |  ETA: 0:00:06
Truth solving on (20, 20) grid...  10%|█▊                |  ETA: 0:00:06
Truth solving on (20, 20) grid...  11%|██                |  ETA: 0:00:06
Truth solving on (20, 20) grid...  13%|██▎               |  ETA: 0:00:06
Truth solving on (20, 20) grid...  14%|██▋               |  ETA: 0:00:06
Truth solving on (20, 20) grid...  16%|██▉               |  ETA: 0:00:06
Truth solving on (20, 20) grid...  18%|███▎              |  ETA: 0:00:06
Truth solving on (20, 20) grid...  20%|███▌              |  ETA: 0:00:06
Truth solving on (20, 20) grid...  21%|███▊              |  ETA: 0:00:05
Truth solving on (20, 20) grid...  23%|████▏             |  ETA: 0:00:05
Truth solving on (20, 20) grid...  24%|████▍             |  ETA: 0:00:05
Truth solving on (20, 20) grid...  26%|████▋             |  ETA: 0:00:05
Truth solving on (20, 20) grid...  28%|█████             |  ETA: 0:00:05
Truth solving on (20, 20) grid...  29%|█████▎            |  ETA: 0:00:05
Truth solving on (20, 20) grid...  30%|█████▌            |  ETA: 0:00:05
Truth solving on (20, 20) grid...  32%|█████▊            |  ETA: 0:00:05
Truth solving on (20, 20) grid...  34%|██████            |  ETA: 0:00:05
Truth solving on (20, 20) grid...  35%|██████▎           |  ETA: 0:00:05
Truth solving on (20, 20) grid...  36%|██████▋           |  ETA: 0:00:04
Truth solving on (20, 20) grid...  38%|██████▉           |  ETA: 0:00:04
Truth solving on (20, 20) grid...  40%|███████▏          |  ETA: 0:00:04
Truth solving on (20, 20) grid...  41%|███████▍          |  ETA: 0:00:04
Truth solving on (20, 20) grid...  43%|███████▊          |  ETA: 0:00:04
Truth solving on (20, 20) grid...  44%|████████          |  ETA: 0:00:04
Truth solving on (20, 20) grid...  46%|████████▎         |  ETA: 0:00:04
Truth solving on (20, 20) grid...  48%|████████▌         |  ETA: 0:00:04
Truth solving on (20, 20) grid...  49%|████████▉         |  ETA: 0:00:04
Truth solving on (20, 20) grid...  50%|█████████▏        |  ETA: 0:00:03
Truth solving on (20, 20) grid...  52%|█████████▍        |  ETA: 0:00:03
Truth solving on (20, 20) grid...  54%|█████████▋        |  ETA: 0:00:03
Truth solving on (20, 20) grid...  55%|██████████        |  ETA: 0:00:03
Truth solving on (20, 20) grid...  57%|██████████▎       |  ETA: 0:00:03
Truth solving on (20, 20) grid...  58%|██████████▌       |  ETA: 0:00:03
Truth solving on (20, 20) grid...  60%|██████████▊       |  ETA: 0:00:03
Truth solving on (20, 20) grid...  62%|███████████▏      |  ETA: 0:00:03
Truth solving on (20, 20) grid...  63%|███████████▍      |  ETA: 0:00:03
Truth solving on (20, 20) grid...  64%|███████████▋      |  ETA: 0:00:03
Truth solving on (20, 20) grid...  66%|███████████▉      |  ETA: 0:00:02
Truth solving on (20, 20) grid...  68%|████████████▎     |  ETA: 0:00:02
Truth solving on (20, 20) grid...  69%|████████████▌     |  ETA: 0:00:02
Truth solving on (20, 20) grid...  71%|████████████▊     |  ETA: 0:00:02
Truth solving on (20, 20) grid...  72%|█████████████     |  ETA: 0:00:02
Truth solving on (20, 20) grid...  74%|█████████████▎    |  ETA: 0:00:02
Truth solving on (20, 20) grid...  75%|█████████████▌    |  ETA: 0:00:02
Truth solving on (20, 20) grid...  77%|█████████████▉    |  ETA: 0:00:02
Truth solving on (20, 20) grid...  78%|██████████████▏   |  ETA: 0:00:02
Truth solving on (20, 20) grid...  80%|██████████████▍   |  ETA: 0:00:01
Truth solving on (20, 20) grid...  81%|██████████████▋   |  ETA: 0:00:01
Truth solving on (20, 20) grid...  83%|██████████████▉   |  ETA: 0:00:01
Truth solving on (20, 20) grid...  84%|███████████████▏  |  ETA: 0:00:01
Truth solving on (20, 20) grid...  86%|███████████████▍  |  ETA: 0:00:01
Truth solving on (20, 20) grid...  87%|███████████████▊  |  ETA: 0:00:01
Truth solving on (20, 20) grid...  89%|████████████████  |  ETA: 0:00:01
Truth solving on (20, 20) grid...  90%|████████████████▎ |  ETA: 0:00:01
Truth solving on (20, 20) grid...  92%|████████████████▌ |  ETA: 0:00:01
Truth solving on (20, 20) grid...  93%|████████████████▊ |  ETA: 0:00:01
Truth solving on (20, 20) grid...  94%|█████████████████ |  ETA: 0:00:00
Truth solving on (20, 20) grid...  96%|█████████████████▎|  ETA: 0:00:00
Truth solving on (20, 20) grid...  98%|█████████████████▌|  ETA: 0:00:00
Truth solving on (20, 20) grid...  99%|█████████████████▉|  ETA: 0:00:00
Truth solving on (20, 20) grid... 100%|██████████████████| Time: 0:00:07

Since we do not compute any Hamiltonian compressions during POD, we need to compute them afterwards using the HamiltonianCache constructor (recall that h is needed in the online stage):

h_cache = HamiltonianCache(H, rbres.basis);

Again, we arrive at the online phase which is performed analogously to The reduced basis workflow.

M = AffineDecomposition([H.terms[3]], [2 / L])
m, _ = compress(M, rbres.basis)
m_reduced = m()
Δ_online = range(first(Δ), last(Δ); length=100)
hJ_online = range(first(hJ), last(hJ); length=100)
grid_online = RegularGrid(Δ_online, hJ_online)
fulldiag = FullDiagonalization(lobpcg)

using Statistics
magnetization = map(grid_online) do μ
    _, φ_rb = solve(h_cache.h, rbres.basis.metric, μ, fulldiag)
    mean(u -> abs(dot(u, m_reduced, u)), eachcol(φ_rb))
end;

With this we can again produce a magnetization diagram:

using Plots
hm = heatmap(grid_online.ranges[1], grid_online.ranges[2], magnetization';
             xlabel=raw"$\Delta$", ylabel=raw"$h/J$", title="magnetization",
             colorbar=true, clims=(0.0, 1.0), leg=false)
plot!(hm, grid_online.ranges[1], x -> 1 + x; lw=2, ls=:dash, legend=false, color=:green)
Example block output

The magnetization phase diagram is correctly reproduced, however this time without the parameter point markers being plotted. This is due to the fact that all points in grid_train have been solved and incorporated into the RB according to the POD procedure. The number of retained singular vectors therefore does not correspond directly to snapshots at certain parameter points but to linear combinations of them.