Changes for non-CPU array support

Author: kshyatt

Project.toml

@@ -76,3 +76,7 @@ cuTENSOR = "011b41b2-24ef-40a8-b3eb-fa098493e9e1"
 
 [targets]
 test = ["Aqua", "Adapt", "CUDA", "cuTENSOR", "Pkg", "Test", "TestExtras", "Plots", "Combinatorics", "ParallelTestRunner", "TensorKitTensors"]
+
+[sources]
+TensorKit = {url="https://github.com/QuantumKitHub/TensorKit.jl", rev="ksh/cuda_tweaks"}
+BlockTensorKit = {url="https://github.com/kshyatt/BlockTensorKit.jl", rev="ksh/conversions"}

examples/cu_classic2d/1.hard-hexagon/main.jl

@@ -0,0 +1,104 @@
+using Markdown
+
+md"""
+# [The Hard Hexagon model](@id demo_hardhexagon)
+
+![logo](hexagon.svg)
+
+Tensor networks are a natural way to do statistical mechanics on a lattice.
+As an example of this we will extract the central charge of the hard hexagon model.
+This model is known to have central charge 0.8, and has very peculiar non-local (anyonic) symmetries.
+Because TensorKit supports anyonic symmetries, so does MPSKit.
+To follow the tutorial you need the following packages.
+"""
+
+using MatrixAlgebraKit, MPSKit, MPSKitModels, TensorKit, Plots, Polynomials, Adapt, CUDA, cuTENSOR
+
+#src # for reproducibility:
+#src using Random
+#src Random.seed!(123)
+
+md"""
+The [hard hexagon model](https://en.wikipedia.org/wiki/Hard_hexagon_model) is a 2-dimensional lattice model of a gas, where particles are allowed to be on the vertices of a triangular lattice, but no two particles may be adjacent.
+This can be encoded in a transfer matrix with a local MPO tensor using anyonic symmetries, and the resulting MPO has been implemented in MPSKitModels.
+
+In order to use these anyonic symmetries, we need to generalise the notion of the bond dimension and define how it interacts with the symmetry. Thus, we implement away of converting integers to symmetric spaces of the given dimension, which provides a crude guess for how the final MPS would distribute its Schmidt spectrum.
+"""
+mpo = adapt(CuArray, hard_hexagon())
+P = physicalspace(mpo, 1)
+function virtual_space(D::Integer)
+    _D = round(Int, D / sum(dim, values(FibonacciAnyon)))
+    return Vect[FibonacciAnyon](sector => _D for sector in (:I, :τ))
+end
+
+@assert isapprox(dim(virtual_space(100)), 100; atol = 3)
+
+md"""
+## The leading boundary
+
+One way to study statistical mechanics in infinite systems with tensor networks is by approximating the dominant eigenvector of the transfer matrix by an MPS. 
+This dominant eigenvector contains a lot of hidden information.
+For example, the free energy can be extracted by computing the expectation value of the mpo.
+Additionally, we can compute the entanglement entropy as well as the correlation length of the state:
+"""
+
+D = 10
+V = virtual_space(D)
+ψ₀ = adapt(CuArray, InfiniteMPS(ComplexF64, [P], [V]))
+
+# put this here to avoid interfering with initial InfiniteMPS construction
+MPSKit.Defaults.alg_qr() = CUSOLVER_HouseholderQR(; positive = true)
+MPSKit.Defaults.alg_svd() = CUSOLVER_QRIteration()
+MPSKit.Defaults.alg_lq() = LQViaTransposedQR(CUSOLVER_HouseholderQR(; positive = true))
+
+ψ, envs, = leading_boundary(
+    ψ₀, mpo,
+    VUMPS(; verbosity = 0, alg_eigsolve = MPSKit.Defaults.alg_eigsolve(; ishermitian = false))
+) # use non-hermitian eigensolver
+F = real(expectation_value(ψ, mpo))
+S = real(first(entropy(ψ)))
+ξ = correlation_length(ψ)
+println("F = $F\tS = $S\tξ = $ξ")
+
+md"""
+## The scaling hypothesis
+
+The dominant eigenvector is of course only an approximation. The finite bond dimension enforces a finite correlation length, which effectively introduces a length scale in the system. This can be exploited to formulate a scaling hypothesis [pollmann2009](@cite), which in turn allows to extract the central charge.
+
+First we need to know the entropy and correlation length at a bunch of different bond dimensions. Our approach will be to re-use the previous approximated dominant eigenvector, and then expanding its bond dimension and re-running VUMPS.
+According to the scaling hypothesis we should have ``S \propto \frac{c}{6} log(ξ)``. Therefore we should find ``c`` using
+"""
+
+function scaling_simulations(
+        ψ₀, mpo, Ds; verbosity = 0, tol = 1.0e-6,
+        alg_eigsolve = MPSKit.Defaults.alg_eigsolve(; ishermitian = false)
+    )
+    entropies = similar(Ds, Float64)
+    correlations = similar(Ds, Float64)
+    alg = VUMPS(; verbosity, tol, alg_eigsolve)
+
+    ψ, envs, = leading_boundary(ψ₀, mpo, alg)
+    entropies[1] = real(entropy(ψ)[1])
+    correlations[1] = correlation_length(ψ)
+
+    for (i, d) in enumerate(diff(Ds))
+        ψ, envs = changebonds(ψ, mpo, OptimalExpand(; trscheme = truncrank(d)), envs)
+        ψ, envs, = leading_boundary(ψ, mpo, alg, envs)
+        entropies[i + 1] = real(entropy(ψ)[1])
+        correlations[i + 1] = correlation_length(ψ)
+    end
+    return entropies, correlations
+end
+
+bond_dimensions = 10:5:25
+ψ₀ = InfiniteMPS(CuVector{ComplexF64}, [P], [virtual_space(bond_dimensions[1])])
+Ss, ξs = scaling_simulations(ψ₀, mpo, bond_dimensions)
+
+f = fit(log.(ξs), 6 * Ss, 1)
+c = f.coeffs[2]
+
+#+
+
+p = plot(; xlabel = "logarithmic correlation length", ylabel = "entanglement entropy")
+p = plot(log.(ξs), Ss; seriestype = :scatter, label = nothing)
+plot!(p, ξ -> f(ξ) / 6; label = "fit")

examples/cu_quantum1d/1.ising-cft/main.jl

@@ -0,0 +1,163 @@
+using Markdown
+md"""
+# The Ising CFT spectrum
+
+This tutorial is meant to show the finite size CFT spectrum for the quantum Ising model. We
+do this by first employing an exact diagonalization technique, and then extending the
+analysis to larger system sizes through the use of MPS techniques.
+"""
+
+using MPSKit, MPSKitModels, TensorKit, Plots, KrylovKit
+using LinearAlgebra: eigvals, diagm, Hermitian
+using CUDA, cuTENSOR, Adapt
+
+md"""
+The Hamiltonian is defined on a finite lattice with periodic boundary conditions,
+which can be implemented as follows:
+"""
+
+L = 12
+H = periodic_boundary_conditions(transverse_field_ising(), L)
+
+md"""
+## Exact diagonalisation
+
+In MPSKit, there is support for exact diagonalisation by leveraging the fact that applying
+the Hamiltonian to an untruncated MPS will result in an effective Hamiltonian on the center
+site which implements the action of the entire Hamiltonian. Thus, optimizing the middle
+tensor is equivalent to optimixing a state in the entire Hilbert space, as all other tensors
+are just unitary matrices that mix the basis.
+"""
+
+energies, states = exact_diagonalization(H; num = 18, alg = Lanczos(; krylovdim = 200));
+plot(
+    real.(energies);
+    seriestype = :scatter, legend = false, ylabel = "energy", xlabel = "#eigenvalue"
+)
+cu_states = map(state -> adapt(CuArray, state), states)
+
+md"""
+!!! note "Krylov dimension"
+    Note that we have specified a large Krylov dimension as degenerate eigenvalues are
+    notoriously difficult for iterative methods.
+"""
+
+md"""
+## Extracting momentum
+
+Given a state, it is possible to assign a momentum label
+through the use of the translation operator. This operator can be defined in MPO language
+either diagramatically as
+
+```@raw html
+<img src="translation_mpo.svg" alt="translation operator" class="color-invertible" width="50%"/>
+```
+
+or in the code as:
+"""
+
+function O_shift(L)
+    I = id(ComplexF64, ℂ^2)
+    @tensor O[W S; N E] := I[W; N] * I[S; E]
+    return adapt(CuArray, periodic_boundary_conditions(InfiniteMPO([O]), L))
+end
+
+md"""
+We can then calculate the momentum of the ground state as the expectation value of this
+operator. However, there is a subtlety because of the degeneracies in the energy
+eigenvalues. The eigensolver will find an orthonormal basis within each energy subspace, but
+this basis is not necessarily a basis of eigenstates of the translation operator. In order
+to fix this, we diagonalize the translation operator within each energy subspace.
+The resulting energy levels have one-to-one correspondence to the operators in CFT, where
+the momentum is related to their conformal spin as $P_n = \frac{2\pi}{L}S_n$.
+"""
+
+function fix_degeneracies(basis)
+    L = length(basis[1])
+    M = zeros(ComplexF64, length(basis), length(basis))
+    T = O_shift(L)
+    for j in eachindex(basis), i in eachindex(basis)
+        M[i, j] = dot(basis[i], T, basis[j])
+    end
+
+    vals = eigvals(M)
+    return angle.(vals)
+end
+
+momenta = Float64[]
+append!(momenta, fix_degeneracies(cu_states[1:1]))
+append!(momenta, fix_degeneracies(cu_states[2:2]))
+append!(momenta, fix_degeneracies(cu_states[3:3]))
+append!(momenta, fix_degeneracies(cu_states[4:5]))
+append!(momenta, fix_degeneracies(cu_states[6:9]))
+append!(momenta, fix_degeneracies(cu_states[10:11]))
+append!(momenta, fix_degeneracies(cu_states[12:12]))
+append!(momenta, fix_degeneracies(cu_states[13:16]))
+append!(momenta, fix_degeneracies(cu_states[17:18]))
+
+md"""
+We can compute the scaling dimensions $\Delta_n$ of the operators in the CFT from the
+energy gap of the corresponding excitations as $E_n - E_0 = \frac{2\pi v}{L} \Delta_n$,
+where $v = 2$. If we plot these scaling dimensions against the conformal spin $S_n$ from
+above, we retrieve the familiar spectrum of the Ising CFT.
+"""
+
+v = 2.0
+Δ = real.(energies[1:18] .- energies[1]) ./ (2π * v / L)
+S = momenta ./ (2π / L)
+
+p = plot(
+    S, real.(Δ);
+    seriestype = :scatter, xlabel = "conformal spin (S)", ylabel = "scaling dimension (Δ)",
+    legend = false
+)
+vline!(p, -3:3; color = "gray", linestyle = :dash)
+hline!(p, [0, 1 / 8, 1, 9 / 8, 2, 17 / 8]; color = "gray", linestyle = :dash)
+p
+
+md"""
+## Finite bond dimension
+
+If we limit the maximum bond dimension of the MPS, we get an approximate solution, but we
+can reach higher system sizes.
+"""
+
+L_mps = 20
+H_mps = adapt(CuArray, periodic_boundary_conditions(transverse_field_ising(), L_mps))
+D = 64
+
+ψ, envs, δ = find_groundstate(adapt(CuArray, FiniteMPS(L_mps, ℂ^2, ℂ^D)), H_mps, DMRG());
+
+md"""
+Excitations on top of the ground state can be found through the use of the quasiparticle
+ansatz. This returns quasiparticle states, which can be converted to regular `FiniteMPS`
+objects.
+"""
+
+E_ex, qps = excitations(H_mps, QuasiparticleAnsatz(), ψ, envs; num = 18)
+states_mps = vcat(ψ, map(qp -> convert(FiniteMPS, qp), qps))
+energies_mps = map(x -> expectation_value(x, H_mps), states_mps)
+
+momenta_mps = Float64[]
+append!(momenta_mps, fix_degeneracies(states_mps[1:1]))
+append!(momenta_mps, fix_degeneracies(states_mps[2:2]))
+append!(momenta_mps, fix_degeneracies(states_mps[3:3]))
+append!(momenta_mps, fix_degeneracies(states_mps[4:5]))
+append!(momenta_mps, fix_degeneracies(states_mps[6:9]))
+append!(momenta_mps, fix_degeneracies(states_mps[10:11]))
+append!(momenta_mps, fix_degeneracies(states_mps[12:12]))
+append!(momenta_mps, fix_degeneracies(states_mps[13:16]))
+append!(momenta_mps, fix_degeneracies(states_mps[17:18]))
+
+v = 2.0
+Δ_mps = real.(energies_mps[1:18] .- energies_mps[1]) ./ (2π * v / L_mps)
+S_mps = momenta_mps ./ (2π / L_mps)
+
+p_mps = plot(
+    S_mps, real.(Δ_mps);
+    seriestype = :scatter, xlabel = "conformal spin (S)",
+    ylabel = "scaling dimension (Δ)", legend = false
+)
+vline!(p_mps, -3:3; color = "gray", linestyle = :dash)
+hline!(p_mps, [0, 1 / 8, 1, 9 / 8, 2, 17 / 8]; color = "gray", linestyle = :dash)
+p_mps