From 2dd3dbcfed5daa43409aefa1bee60a3870b0d7ec Mon Sep 17 00:00:00 2001 From: TymoteuszTula Date: Mon, 13 Jul 2026 14:00:26 +0200 Subject: [PATCH 1/3] Added MPO compression for Finite and Infinite MPO Hamiltonians --- Project.toml | 2 + src/MPSKit.jl | 5 +- src/operators/mpocompression.jl | 634 +++++++++++++++++++++ test/operators/finite_mpo_compression.jl | 375 ++++++++++++ test/operators/infinite_mpo_compression.jl | 388 +++++++++++++ 5 files changed, 1403 insertions(+), 1 deletion(-) create mode 100644 src/operators/mpocompression.jl create mode 100644 test/operators/finite_mpo_compression.jl create mode 100644 test/operators/infinite_mpo_compression.jl diff --git a/Project.toml b/Project.toml index 2a2d97ebc..61f60f476 100644 --- a/Project.toml +++ b/Project.toml @@ -10,6 +10,7 @@ projects = ["test", "docs", "examples"] Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" BlockTensorKit = "5f87ffc2-9cf1-4a46-8172-465d160bd8cd" Compat = "34da2185-b29b-5c13-b0c7-acf172513d20" +Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" HalfIntegers = "f0d1745a-41c9-11e9-1dd9-e5d34d218721" KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" @@ -38,6 +39,7 @@ Accessors = "0.1" Adapt = "4" BlockTensorKit = "0.3.14" Compat = "3.47, 4.10" +Dictionaries = "0.4.6" DocStringExtensions = "0.9.3" HalfIntegers = "1.6.0" KrylovKit = "0.8.3, 0.9.2, 0.10" diff --git a/src/MPSKit.jl b/src/MPSKit.jl index b7387f3da..3c8da7226 100644 --- a/src/MPSKit.jl +++ b/src/MPSKit.jl @@ -49,6 +49,7 @@ export entropy, entanglement_spectrum export open_boundary_conditions, periodic_boundary_conditions export entanglementplot, transferplot export r_LL, l_LL, r_RR, l_RR, r_RL, r_LR, l_RL, l_LR # TODO: rename +export mpo_compression # unexported using Compat: @compat @@ -79,13 +80,14 @@ using Accessors using HalfIntegers using DocStringExtensions -using LinearAlgebra: diag, Diagonal +using LinearAlgebra: diag, Diagonal, diagm using LinearAlgebra: LinearAlgebra using Random using Base: @kwdef, @propagate_inbounds using LoggingExtras using OhMyThreads using TimerOutputs: TimerOutput, @timeit, timeit, reset_timer!, disable_timer!, enable_timer! +using Dictionaries # TODO: change mpocompression to not use Dictionaries # Includes # -------- @@ -130,6 +132,7 @@ include("operators/projection.jl") include("operators/timedependence.jl") include("operators/multipliedoperator.jl") include("operators/lazysum.jl") +include("operators/mpocompression.jl") include("transfermatrix/transfermatrix.jl") include("transfermatrix/transfer.jl") diff --git a/src/operators/mpocompression.jl b/src/operators/mpocompression.jl new file mode 100644 index 000000000..9691abf3f --- /dev/null +++ b/src/operators/mpocompression.jl @@ -0,0 +1,634 @@ +""" + mpo_compression(H::FiniteMPOHamiltonian; η::Number=10^-8) -> Hᵪ, Rs + mpo_compression(H::InfiniteMPOHamiltonian; η::Number=10^-8) -> Qs, Ps + +Returns the compressed version of the FiniteMPOHamiltonian or InfiniteMPOHamiltonian, +as described in: https://arxiv.org/pdf/1909.06341. + +### Arguments +- `H`: Hamiltonian to be compressed. Either FiniteMPOHamiltonian or InfiniteMPOHamiltonian. + +### Keyword arguments +- `η::Number=10^-8`: precision of the compression. The values smaller than η in the SVD + decomposition during compression will be disregarded. + +### Returns +- mpo_compression(H::FiniteMPOHamiltonian; η::Number=10^-8) + `Hᵪ::FiniteMPOHamiltonian`: compressed Hamiltonian in the left canonical form. + `Rs::Vector{BlockTensorMap}`: matrices that can be applied to change between left and + right canonical forms of Hᵪ. +- mpo_compression(H::InfiniteMPOHamiltonian; η::Number=10^-8) + `Qs::InfiniteMPOHamiltonian`: compressed Hamiltonian in the left canonical form. + `Ps::InfiniteMPOHamiltonian`: compressed Hamiltonian in the right canonical form. +""" + +function mpo_compression end + +function mpo_compression(H::FiniteMPOHamiltonian, η::Number=10^-8) + t = @elapsed begin + # make sure that the MPO have at least three sites + N = length(H) + @assert N ≥ 3 + # change MPO from N x 1 x 1 x M blocks to 3 x 1 x 1 x 3 block structure + TT = TensorMap{scalartype(H), spacetype(H[1]), 2, 2, Vector{scalartype(H)}} + Hnew = Vector{SparseBlockTensorMap{TT, scalartype(H), spacetype(H[1]), 2, 2, 4}}( + undef, N) + for n in 1:N + Hnew[n] = SparseBlockTensorMap(reduce_blocks_mpo(H[n])) + end + + Qs, Ls = right_canonical_mpo_finite(Hnew) + TT = AbstractTensorMap{scalartype(H), spacetype(H), 2, 2} + Cs = Vector{SparseBlockTensorMap{TT, scalartype(H), spacetype(H), 2, 2, 4}}(undef, N) + Cs[1] = Qs[1] + Rs = Vector{AbstractBlockTensorMap}(undef, N-2) + R = id(codomain(Qs[2])[1]) + for n in 2:N-1 + @tensor Qnew[a,i;j,b] := R[a,c] * Qs[n][c,i;j,b] + Q, R = qr_block_respecting(Qnew) + # For now it is safe to assume that R will be 3x3 blocked BlockTensorMap, + # since it went through right_canonical_mpo_finite first + M, R′ = decompose_R(R) + U, S, V′ = block_respecting_svd(M, η) + @tensor C[a,i;j,b] := Q[a,i;j,c] * U[c,b] + @tensor R[a;b] := S[a,c] * V′[c,d] * R′[d,b] + Cs[n] = C + Rs[n-1] = R + end + @tensor Cnew[a,i;j,b] := R[a,c] * Qs[N][c,i;j,b] + Cs[N] = Cnew + # Last step: Make new Cs into FiniteMPOHamiltonian + # TA = AbstractTensorMap{scalartype(H), spacetype(H), 2, 2} + # TB = AbstractTensorMap{scalartype(H), spacetype(H), 2, 1} + # TC = AbstractTensorMap{scalartype(H), spacetype(H), 1, 2} + # TD = AbstractTensorMap{scalartype(H), spacetype(H), 1, 1} + Hᵪtype = JordanMPOTensor{scalartype(H), spacetype(H), Vector{scalartype(H)}} + Hᵪ = Vector{Hᵪtype}(undef, N) + for n in 1:N + W = SparseBlockTensorMap(Cs[n]) + A = W[2:(end-1), 1, 1, 2:(end-1)] + B = removeunit(W[2:(end-1), 1, 1, end], 4) + C = removeunit(W[1,1,1,2:(end-1)], 1) + D = removeunit(removeunit(W[1,1,1,end:end], 4), 1) + W = JordanMPOTensor(space(W), A, B, C, D) + Hᵪ[n] = W + end + Hᵪ = FiniteMPOHamiltonian(Hᵪ) + end + # Print reduction data + tot_dim_original = total_virt_dimension(H) + tot_dim_compressed = total_virt_dimension(Hᵪ) + printstyled("┌───────────────────────────────────────\n", color=:cyan) + printstyled("| MPO Compression performed succesfully\n", color=:cyan) + printstyled("│ ⏱ Time: ", color=:cyan) + printstyled("$(round(t, digits=1)) s\n", color=:yellow) + printstyled("| Initial total virtual dimensions: ", color=:cyan) + printstyled("$(tot_dim_original)\n", color=:yellow) + printstyled("| Final total virtual dimensions: ", color=:cyan) + printstyled("$(tot_dim_compressed)\n", color=:yellow) + printstyled("└───────────────────────────────────────\n") + + return Hᵪ, Rs +end + +function mpo_compression(H::InfiniteMPOHamiltonian, η::Number=10^-8) + tot_dim_original = total_virt_dimension(H) + t = @elapsed begin + N = length(H) # Number of sites + # change MPO from N x 1 x 1 x M blocks to 3 x 1 x 1 x 3 block structure + TT = TensorMap{scalartype(H), spacetype(H[1]), 2, 2, Vector{scalartype(H)}} + Hnew = Vector{SparseBlockTensorMap{TT, scalartype(H), spacetype(H[1]), 2, 2, 4}}( + undef, N) + for n in 1:N + Hnew[n] = SparseBlockTensorMap(reduce_blocks_mpo(H[n])) + end + + HL, _ = left_canonical_mpo_infinite_iter_msites(Hnew) + _, HR = right_canonical_mpo_infinite_iter_msites(HL) + HL, Cs = left_canonical_mpo_infinite_iter_msites(HR) + Us = Vector{BlockTensorMap}(undef, N) + Ss = Vector{BlockTensorMap}(undef, N) + Vs = Vector{BlockTensorMap}(undef, N) + for n = 1:N + Us[n], Ss[n], Vs[n] = block_respecting_svd(Cs[n], η) + end + Qs = Vector{BlockTensorMap}(undef, N) + Ps = Vector{BlockTensorMap}(undef, N) + for n = 1:N + # Left canonical version of compressed iMPO + @tensor Q[a,i;j,b] := Us[mod1(n-1,N)]'[a,c] * HL[n][c,i;j,d] * Us[n][d,b] + # Right canonical version of compressed iMPO + @tensor P[a,i;j,b] := Vs[mod1(n-1,N)][a,c] * HR[n][c,i;j,d] * Vs[n]'[d,b] + Qs[n] = Q + Ps[n] = P + end + + Q̂s = create_impoham_from_mpos_msites(Qs) + P̂s = create_impoham_from_mpos_msites(Ps) + end + + # Print reduction data + tot_dim_compressed = total_virt_dimension(Q̂s) + printstyled("┌───────────────────────────────────────\n", color=:cyan) + printstyled("| MPO Compression performed succesfully\n", color=:cyan) + printstyled("│ ⏱ Time: ", color=:cyan) + printstyled("$(round(t, digits=1)) s\n", color=:yellow) + printstyled("| Initial total virtual dimensions: ", color=:cyan) + printstyled("$(tot_dim_original)\n", color=:yellow) + printstyled("| Final total virtual dimensions: ", color=:cyan) + printstyled("$(tot_dim_compressed)\n", color=:yellow) + printstyled("└───────────────────────────────────────\n") + + return Q̂s, P̂s +end + +# utility +# ------- +function trace_single_block(d::AbstractTensorMap) + return @tensor trd[a;b] := d[a,i;i,b] +end + +function change_value_at_fusiontree!(T::AbstractTensorMap, i::Integer, v::Number) + # change value of the fusiontree at index i of a TensorMap T with new value v + f = fusiontrees(T) + if f isa Vector + T[f[i]...] .= v + elseif f isa Indices + T[f.values[i]...] .= v + else + throw(ArgumentError("fusiontrees of T is neither a Vector or a Indices. It is + a $(typeof(f))")) + end +end + +function total_virt_dimension(H::FiniteMPOHamiltonian) + # calculate the sum of bond dimensions for entire H MPO + N = length(H) + tvd = sum([dims(H[i])[4] for i in 1:N-1]) + return tvd +end + +function total_virt_dimension(H::InfiniteMPOHamiltonian) + # calculate the sum of bond dimensions for entire H MPO + N = length(H) + tvd = sum([dims(H[i])[4] for i in 1:N]) + return tvd +end + +# functions +# --------- +function qr_block_respecting(W::AbstractBlockTensorMap) + # Block respecting qr (qr applied only to the upper-left part of the matrix) + bl_st = size(W) + dim_st = dims(W) + + # sanity checks + @assert length(bl_st) == length(dim_st) == 4 + blN = bl_st[1] + blM = bl_st[4] + + # First, for each nonzero fusiontree from blocks in the first row, make them + # traceless + id1 = id(codomain(W)[2]) + + t = zeros(ComplexF64, domain(W)[2][1] ← domain(W)[2][2:end-1]) + for block in 2:blM-1 + d = W[1,1,1,block] + trd = trace_single_block(d) + @tensor W[1,1,1,block][a,i;j,b] = d[a,i;j,b] - trd[a;b] * id1[i;j] + # Here goes the creation of matrix R' with a t' tensor + t[block-1] = trd / dim_st[2] + end + # Perform QR of the 2:blM-1 block + V = W[1:blN-1, :, :, 2:blM-1] + V = permute(V, ((1,2,3), (4,))) + + V = convert(SparseBlockTensorMap{TensorMap{eltype(V).parameters...}, + eltype(storagetype(V)), spacetype(V), 3, 1, 4,}, V) + Q, R = qr_compact(V) + + Q = permute(Q, ((1,2),(3,4))) + + # Put Q into larger Q̂ + Q̂ = zeros(scalartype(W), codomain(W) ← domain(W)[1] ⊗ (domain(W)[2][1] ⊞ domain(Q)[2] ⊞ + domain(W)[2][end])) + Q̂[1,1,1,1] = Q̂[end,1,1,end] = W[1,1,1,1] + Q̂[1:end-1,1,1,2:end-1] = Q + Q̂[1:end-1,1,1,end] = W[1:end-1,1,1,end] + + # Put R and t into larger R̂ + R̂ = zeros(scalartype(W), domain(W)[2][1] ⊞ codomain(R)[1] ⊞ domain(W)[2][end] ← + domain(W)[2][1] ⊞ domain(R)[1] ⊞ domain(W)[2][end]) + R̂[2:end-1,2:end-1] = R + R̂[1,2:end-1] = t + change_value_at_fusiontree!(R̂[1,1],1,1) + change_value_at_fusiontree!(R̂[end,end],1,1) + + return Q̂, R̂ +end + +function lq_block_respecting(W::AbstractBlockTensorMap) + # Block respecting lq (lq applied only to the bottom-right part of the matrix) + bl_st = size(W) + dim_st = dims(W) + + # sanity checks + @assert length(bl_st) == length(dim_st) == 4 + blN = bl_st[1] + blM = bl_st[4] + # First, for each nonzero fusiontree from blocks in the last column, make them + # traceless + id1 = id(codomain(W)[2]) + + t = zeros(ComplexF64, codomain(W)[1][2:end-1] ← codomain(W)[1][end]) + for block in 2:blN-1 + d = W[block,1,1,end] + trd = trace_single_block(d) + @tensor W[block,1,1,end][a,i;j,b] = d[a,i;j,b] - trd[a;b] * id1[i;j] + # Here goes the creation of matrix R' with a t' tensor + t[block-1] = trd / dim_st[2] + end + # Perform QR of the 2:blM-1 block (with virtual indices transposed) + V = W[2:blN-1, :, :, 2:blM] + V = permute(V, ((1,), (2,3,4))) + + V = convert(SparseBlockTensorMap{TensorMap{eltype(V).parameters...}, + eltype(storagetype(V)), spacetype(V), 1, 3, 4,}, V) + # to perform QR + L, Q = lq_compact(V) + # transpose vitrual indices back (now to reconstruct initial matrix + # one has to apply R on the left) + Q = permute(Q, ((1,2),(3,4))) + + # Put Q into larger Q̂ + Q̂ = zeros(scalartype(W), (codomain(W)[1][1] ⊞ codomain(Q)[1] ⊞ codomain(W)[1][end]) ⊗ + codomain(W)[2] ← domain(W)) + Q̂[1,1,1,1] = Q̂[end,1,1,end] = W[1,1,1,1] + Q̂[2:end-1,1,1,2:end] = Q + Q̂[1,1,1,2:end] = W[1,1,1,2:end] + + # Put L and t into larger L̂ + L̂ = zeros(scalartype(W), domain(W)[2][1] ⊞ codomain(L)[1] ⊞ domain(W)[2][end] ← + domain(W)[2][1] ⊞ domain(L)[1] ⊞ domain(W)[2][end]) + L̂[2:end-1,2:end-1] = L + L̂[2:end-1,end] = t + change_value_at_fusiontree!(L̂[1,1],1,1) + change_value_at_fusiontree!(L̂[end,end],1,1) + + return L̂, Q̂ +end + +function left_canonical_mpo_finite( + H::Union{FiniteMPOHamiltonian, V} + ) where {T<:SparseBlockTensorMap{<:Any, <:Any, <:Any, 2, 2, 4}, V<:Vector{T}} + # change finite Hamiltonian MPO into its Left canonical form. The 'H' is a + # FiniteMPOHamiltonian from MPSKit (l and r boundary conditions are already + # part of the left-most and right-most MPOs, they can not be compressed, and + # they can be excluded from left canonicalisation algorithm). + + # make sure that the MPO have at least three sites + N = length(H) + @assert N ≥ 3 + + Q, R = qr_block_respecting(H[2]) + TT = AbstractTensorMap{scalartype(H), spacetype(H[1]), 2, 2} + Qs = Vector{SparseBlockTensorMap{TT, scalartype(H), spacetype(H[1]), 2, 2, 4}}(undef, N) + Qs[1] = H[1] + Qs[2] = Q + Rs = Vector{AbstractBlockTensorMap}(undef, N-2) + Rs[1] = R + for n in 3:N-1 + @tensor Wnew[a,i;j,b] := R[a,c] * H[n][c,i;j,b] + Q, R = qr_block_respecting(Wnew) + Qs[n] = Q + Rs[n-1] = R + end + @tensor Wlast[a,i;j,b] := R[a,c] * H[N][c,i;j,b] + Qs[N] = Wlast + return Qs, Rs +end + +function right_canonical_mpo_finite( + H::Union{FiniteMPOHamiltonian, V} + ) where {T<:SparseBlockTensorMap{<:Any, <:Any, <:Any, 2, 2, 4}, V<:Vector{T}} + # change finite Hamiltonian MPO into its Right canonical form. The 'H' is a + # FiniteMPOHamiltonian from MPSKit (l and r boundary conditions are already + # part of the left-most and right-most MPOs, they can not be compressed, and + # they can be excluded from right canonicalisation algorithm). + + # make sure that the MPO have at least three sites + N = length(H) + @assert N ≥ 3 + + L, Q = lq_block_respecting(H[end-1]) + TT = AbstractTensorMap{scalartype(H), spacetype(H[1]), 2, 2} + Qs = Vector{SparseBlockTensorMap{TT, scalartype(H), spacetype(H[1]), 2, 2, 4}}(undef, N) + Qs[end] = H[end] + Qs[end-1] = Q + Ls = Vector{AbstractBlockTensorMap}(undef, N-2) + Ls[end] = L + for n in N-2:-1:2 + @tensor Wnew[a,i;j,b] := H[n][a,i;j,c] * L[c,b] + L, Q = lq_block_respecting(Wnew) + Qs[n] = Q + Ls[n-1] = L + end + @tensor Wlast[a,i;j,b] := H[1][a,i;j,c] * L[c,b] + Qs[1] = Wlast + return Qs, Ls +end + +function decompose_R(R::AbstractBlockTensorMap) + M = zeros(ComplexF64, codomain(R) ← domain(R)) + R′ = zeros(ComplexF64, domain(R) ← domain(R)) + change_value_at_fusiontree!(M[1,1], 1, 1) + change_value_at_fusiontree!(M[end,end], 1, 1) + change_value_at_fusiontree!(R′[1,1], 1, 1) + change_value_at_fusiontree!(R′[end,end], 1, 1) + M[2:end-1,2:end-1] = R[2:end-1,2:end-1] + R′[1,2:end-1] = R[1,2:end-1] + # R′[2:end-1,2:end-1] = id(codomain(R[2,2])) + # Make diagonal matrix even with codomain and domain being different spaces + # Right now I'm assuming that each of the fusiontree is a square block + for (s, f) in fusiontrees(R′[2:end-1,2:end-1]) + if s.uncoupled[1] == f.uncoupled[1] + dim = size(R′[2:end-1,2:end-1][s,f])[1] + R′[2:end-1,2:end-1][s,f] .= diagm(ones(dim)) + end + end + return M, R′ +end + +function block_respecting_svd(M::AbstractBlockTensorMap, η::Number; perform_truncation=true) + # Perform block-respecting svd on matrix M, with precision η + # TODO Change it so that it works when sectors are numbered by indices rather than + # enumerated + M2 = M[2:end-1,2:end-1] + U, S, V′ = svd_compact(M2) + + num_sectors = length(blocksectors(S)) + new_dims = Vector{Int32}(undef, num_sectors) + + # Old code - not working if blocksectors(S) returned Indices rather than Vector + # for n in 1:num_sectors + # csec = blocksectors(S)[n] + # new_dims[n] = blockdim(S.domain, csec) + # if perform_truncation + # for diag_el in block(S, csec).diag + # if diag_el ≤ η + # new_dims[n] -= 1 + # end + # end + # end + # end + + # New code - works for both cases + for (n, csec) in enumerate(blocksectors(S)) + new_dims[n] = blockdim(S.domain, csec) + if perform_truncation + for diag_el in block(S, csec).diag + if diag_el ≤ η + new_dims[n] -= 1 + end + end + end + end + + fdims = Vector{Pair{sectortype(S), Int32}}() + # Old code + # for n in 1:num_sectors + # if new_dims[n] > 0 + # push!(fdims, blocksectors(S)[n] => new_dims[n]) + # end + # end + + # New code - works for both cases + for (n, csec) in enumerate(blocksectors(S)) + if new_dims[n] > 0 + push!(fdims, csec => new_dims[n]) + end + end + red_space = spacetype(S)(fdims...) + S_new = DiagonalTensorMap(undef, red_space) + U_new = zeros(ComplexF64, codomain(U) ← red_space) + V′_new = zeros(ComplexF64, red_space ← domain(V′)) + + #num_sectors2 = length(blocksectors(S_new)) + + # Old code + # for n in 1:num_sectors2 + # csec = blocksectors(S_new)[n] + # cl = length(block(S_new, csec).diag) + # block(S_new, csec).diag[1:end] = block(S, csec).diag[1:cl] + # block(U_new, csec)[:,:] = block(U, csec)[:,1:cl] + # block(V′_new, csec)[:,:] = block(V′, csec)[1:cl,:] + # end + + # New code - works for both cases + for csec in blocksectors(S_new) + cl = length(block(S_new, csec).diag) + block(S_new, csec).diag[1:end] = block(S, csec).diag[1:cl] + block(U_new, csec)[:,:] = block(U, csec)[:,1:cl] + block(V′_new, csec)[:,:] = block(V′, csec)[1:cl,:] + end + + Û = zeros(ComplexF64, codomain(M) ← domain(M)[1][1] ⊞ domain(U_new)[1] ⊞ domain(M)[1][end]) + Ŝ = zeros(ComplexF64, codomain(M)[1][1] ⊞ codomain(S_new)[1] ⊞ codomain(M)[1][end] ← + domain(M)[1][1] ⊞ domain(S_new)[1] ⊞ domain(M)[1][end]) + V̂′ = zeros(ComplexF64, codomain(M)[1][1] ⊞ codomain(V′_new)[1] ⊞ codomain(M)[1][end] ← domain(M)) + + change_value_at_fusiontree!(Û[1,1], 1, 1) + change_value_at_fusiontree!(Û[end,end], 1, 1) + change_value_at_fusiontree!(Ŝ[1,1], 1, 1) + change_value_at_fusiontree!(Ŝ[end,end], 1, 1) + change_value_at_fusiontree!(V̂′[1,1], 1, 1) + change_value_at_fusiontree!(V̂′[end,end], 1, 1) + + Û[2:end-1,2:end-1] = U_new + Ŝ[2,2] = S_new + V̂′[2:end-1,2:end-1] = V′_new + + return Û, Ŝ, V̂′ +end + +function reduce_blocks_mpo(W::AbstractBlockTensorMap) + # Function that copies BlockTensorMap W from N x 1 x 1 x M blocked matrix to a + # 3 x 1 x 1 x 3 blocked matrix, with the middle block being joined dimensions of N-2, + # M-2 blocks from the original matrix + # TODO Make exception if either N or M is smaller than 3 + + # If the matrix is already blocked properly, return it + if size(W)[1] ≤ 3 && size(W)[4] ≤ 3 + return W + end + + old_cod = codomain(W)[1][2:end-1] + old_dom = domain(W)[2][2:end-1] + + new_cod = spacetype(W)([sec => blockdim(old_cod, sec) for sec in sectors(old_cod)]...) + new_dom = spacetype(W)([sec => blockdim(old_dom, sec) for sec in sectors(old_dom)]...) + + # Create new matrix with zeros + if size(W)[1] ≤ 3 + Vcod = codomain(W) + else + Vcod = (codomain(W)[1][1] ⊞ new_cod ⊞ codomain(W)[1][end]) ⊗ codomain(W)[2] + end + if size(W)[4] ≤ 3 + Vdom = domain(W) + else + Vdom = domain(W)[1] ⊗ (domain(W)[2][1] ⊞ new_dom ⊞ domain(W)[2][end]) + end + W2 = zeros(scalartype(W), Vcod ← Vdom) + + # Corners can be copied trivially + W2[1,1,1,1] = W[1,1,1,1] + W2[end,1,1,end] = W[end,1,1,end] + W2[1,1,1,end] = W[1,1,1,end] + W2[end,1,1,1] = W[end,1,1,1] + # Treat the cases with smaller codomain or domain separately + if size(W)[1] ≤ 2 + # Middle upperblock + for (s,f) in fusiontrees(W[1,1,1,2:end-1]) + W2[1,1,1,2][s,f] = W[1,1,1,2:end-1][s,f] + end + # Middle bottomblock (if size(W)[1] is 1, this is the same block as in + # previous loop) + for (s,f) in fusiontrees(W[end,1,1,2:end-1]) + W2[end,1,1,2][s,f] = W[end,1,1,2:end-1][s,f] + end + elseif size(W)[4] ≤ 2 + # Middle leftblock + for (s,f) in fusiontrees(W[2:end-1,1,1,1]) + W2[2,1,1,1][s,f] = W[2:end-1,1,1,1][s,f] + end + # Middle rightblock (if size(W)[4] is 1, this is the same block as in previous loop) + for (s,f) in fusiontrees(W[2:end-1,1,1,end]) + W2[2,1,1,end][s,f] = W[2:end-1,1,1,end][s,f] + end + else + # All other blocks can be dealt with fusiontrees + # Middle block (most important) + for (s,f) in fusiontrees(W[2:end-1,1,1,2:end-1]) + W2[2,1,1,2][s,f] = W[2:end-1,1,1,2:end-1][s,f] + end + # Upper block + for (s,f) in fusiontrees(W[1,1,1,2:end-1]) + W2[1,1,1,2][s,f] = W[1,1,1,2:end-1][s,f] + end + # Bottom block + for (s,f) in fusiontrees(W[end,1,1,2:end-1]) + W2[3,1,1,2][s,f] = W[end,1,1,2:end-1][s,f] + end + # Left block + for (s,f) in fusiontrees(W[2:end-1,1,1,1]) + W2[2,1,1,1][s,f] = W[2:end-1,1,1,1][s,f] + end + # Right block + for (s,f) in fusiontrees(W[2:end-1,1,1,end]) + W2[2,1,1,3][s,f] = W[2:end-1,1,1,end][s,f] + end + end + return W2 +end + +function check_how_far_from_id(R::AbstractBlockTensorMap) + # Function checking how far the diagonal elements of R are from either + # +1 or -1. Assume R is 3x3 and [1,1] and [3,3] blocks are always id + tdev = 0 + for (s,f) in fusiontrees(R[2,2]) + if s.uncoupled[1] == f.uncoupled[1] + for d in size(R[2,2][s,f])[1] + tdev += min(abs(R[2,2][s,f][d,d] - 1), abs(R[2,2][s,f][d,d] + 1)) + end + end + end + return tdev +end + +function create_impoham_from_mpos_msites(Qs::Vector{BlockTensorMap}) + # Given the Vector{BlockTensorMap} Qs, return the multisite InfiniteMPOHamiltonian + N = length(Qs) + sctype = scalartype(Qs[1]) + sptype = spacetype(Qs[1]) + # TA = AbstractTensorMap{sctype, sptype, 2, 2} + # TB = AbstractTensorMap{sctype, sptype, 2, 1} + # TC = AbstractTensorMap{sctype, sptype, 1, 2} + # TD = AbstractTensorMap{sctype, sptype, 1, 1} + CType = JordanMPOTensor{sctype, sptype, Vector{sctype}} + Cs = Vector{CType}(undef, N) + for n in 1:N + W = SparseBlockTensorMap(Qs[n]) + A = W[2:(end-1), 1, 1, 2:(end-1)] + B = removeunit(W[2:(end-1), 1, 1, end], 4) + C = removeunit(W[1,1,1,2:(end-1)], 1) + D = removeunit(removeunit(W[1,1,1,end:end], 4), 1) + W = JordanMPOTensor(space(W), A, B, C, D) + Cs[n] = W + end + Hnew = InfiniteMPOHamiltonian(Cs) + return Hnew +end + +function left_canonical_mpo_infinite_iter_msites( + H::Union{InfiniteMPOHamiltonian, V}, + η=10^-10 + ) where {T<:SparseBlockTensorMap{<:Any, <:Any, <:Any, 2, 2, 4}, V<:Vector{T}} + # Return the left_canonical form of the iMPO, H can have multiple sites in the unit cell + N = length(H) # Number of sites + + εs = fill(Inf, N) + Ls = fill(id(ComplexF64, codomain(H[1])[1]), N) + Rs = fill(id(ComplexF64, codomain(H[1])[1]), N) + Q = similar(H[1]) + printstyled("┌─── Started iterative left orthogonalization ────\n", color=:cyan) + while sum(εs)/N > η + # First store the R matrices and change them H into Qs + for n in 1:N + Q, R = qr_block_respecting(H[n]) + H[n] = Q + Rs[n] = R + Ls[n] = R * Ls[n] + εs[n] = check_how_far_from_id(R) + end + # Then assign R[(n+1)//N] * Q[n] to H[n] + for n in 1:N + @tensor H[n][a,i;j,b] = Rs[mod1(n-1,N)][a;c] * H[n][c,i;j,b] + end + printstyled("| total err = $(sum(εs)/N)\n", color=:white) + end + printstyled("└─── Left orthogonalization finished correctly ──\n", color=:cyan) + return H, Ls +end + +function right_canonical_mpo_infinite_iter_msites( + H::Union{InfiniteMPOHamiltonian, V}, + η=10^-10 + ) where {T<:SparseBlockTensorMap{<:Any, <:Any, <:Any, 2, 2, 4}, V<:Vector{T}} + # Return the right_canonical form of the iMPO, H can have multiple sites in the unit + # cell + N = length(H) # Number of sites + + εs = fill(Inf, N) + Rs = fill(id(ComplexF64, codomain(H[1])[1]), N) + Ls = fill(id(ComplexF64, codomain(H[1])[1]), N) + Q = similar(H[N]) + printstyled("┌─── Started iterative right orthogonalization ────\n", color=:cyan) + while sum(εs)/N > η + # First store the R matrices and change them H into Qs + for n in N:-1:1 + L, Q = lq_block_respecting(H[n]) + H[n] = Q + Ls[n] = L + Rs[n] = Rs[n] * L + εs[n] = check_how_far_from_id(L) + end + # Then assign R[(n+1)//N] * Q[n] to H[n] + for n in N:-1:1 + @tensor H[n][a,i;j,b] = H[n][a,i;j,c] * Ls[mod1(n+1, N)][c,b] + end + printstyled("| total err = $(sum(εs)/N)\n", color=:white) + end + printstyled("└─── Right orthogonalization finished correctly ──\n", color=:cyan) + return Rs, H +end \ No newline at end of file diff --git a/test/operators/finite_mpo_compression.jl b/test/operators/finite_mpo_compression.jl new file mode 100644 index 000000000..07c259cf6 --- /dev/null +++ b/test/operators/finite_mpo_compression.jl @@ -0,0 +1,375 @@ +using Test, Random +using TensorKit +using MPSKit +using MPSKitModels + +function max_dim_ham(H) + N = length(H) + max_dim = 1 + for i in 1:N + if max_dim < dims(H[i])[4] + max_dim = dims(H[i])[4] + end + end + return max_dim +end + +# Utility functions to create Hamiltonians for testing + +function create_long_range_ising_symmetries_finite(M, k) + chain = fill(Z2Space(0=>1, 1=>1), M) + ZZ = S_zz(Z2Irrep) + X = S_x(Z2Irrep) + single_site_operators = [i => X for i in 1:M] + two_site_operators = [(i,i+j) => ZZ for i in 1:M for j in 1:k if i+j ≤ M] + H = FiniteMPOHamiltonian(chain, single_site_operators..., + two_site_operators...) + return H +end + +function create_long_range_ising_symmetries_random_finite(M, k, σⱼ, σₕ) + chain = fill(Z2Space(0=>1, 1=>1), M) + ZZ = S_zz(Z2Irrep) + X = S_x(Z2Irrep) + single_site_operators = [i => (σₕ * randn()) * X for i in 1:M] + two_site_operators = [(i,i+j) => (1 + σⱼ * randn()) * ZZ for i in 1:M for j in 1:k if i+j ≤ M] + H = FiniteMPOHamiltonian(chain, single_site_operators..., + two_site_operators...) + return H +end + +function transverse_field_ising_fermion_parity_finite(M) + chain = FiniteChain(M) + return transverse_field_ising(ComplexF64, FermionParity, chain) +end + +function transverse_field_ising_z2_finite(M) + chain = FiniteChain(M) + return transverse_field_ising(ComplexF64, Z2Irrep, chain) +end + +function transverse_field_ising_trivial_finite(M) + chain = FiniteChain(M) + return transverse_field_ising(ComplexF64, Trivial, chain) +end + +function kitaev_model_finite(M) + chain = FiniteChain(M) + return kitaev_model(ComplexF64, chain) +end + +function heisenberg_XXX_trivial_finite(M) + chain = FiniteChain(M) + return heisenberg_XXX(Trivial, chain) +end + +function heisenberg_XXX_U1_finite(M) + chain = FiniteChain(M) + return heisenberg_XXX(U1Irrep, chain) +end + +function heisenberg_XXX_SU2_finite(M) + chain = FiniteChain(M) + return heisenberg_XXX(SU2Irrep, chain) +end + +function bilinear_biquadratic_model_trivial_finite(M) + chain = FiniteChain(M) + return bilinear_biquadratic_model(ComplexF64, Trivial, chain) +end + +function bilinear_biquadratic_model_U1_finite(M) + chain = FiniteChain(M) + return bilinear_biquadratic_model(ComplexF64, U1Irrep, chain) +end + +function bilinear_biquadratic_model_SU2_finite(M) + chain = FiniteChain(M) + return bilinear_biquadratic_model(ComplexF64, SU2Irrep, chain) +end + +function quantum_potts_trivial_finite(M) + chain = FiniteChain(M) + return quantum_potts(ComplexF64, Trivial, chain) +end + +function quantum_potts_Z3_finite(M) + chain = FiniteChain(M) + return quantum_potts(ComplexF64, Z3Irrep, chain) +end + +# Tests + +@testset "finite_mpo_compression_energy" begin + # test the finite compression for a long range Ising model + # check if the energy after the compression is the same as before + D = 20 # max bond dimensions + L = 10 # number of sites + range = 4 # cutoff for the maximum range between interactions + trunc_st = notrunc() + + mps = FiniteMPS(L, Z2Space(0=>1, 1=>1), Z2Space(0=>D, 1=>D)) + H = create_long_range_ising_symmetries_finite(L, range) + find_groundstate!(mps, H, DMRG2(; maxiter=100, trscheme=trunc_st)) + E0 = expectation_value(mps, H) + #println(" = $real(E0)") + + mps2 = FiniteMPS(L, Z2Space(0=>1, 1=>1), Z2Space(0=>D, 1=>D)) + H2, Rs = mpo_compression(H) + find_groundstate!(mps2, H2, DMRG2(; maxiter=100, trscheme=trunc_st)) + E1 = expectation_value(mps2, H2) +end + +@testset "finite_mpo_compression_transverse_ising_model_energy" begin + M = 10 + D = 10 # Max bond dimension + trunc_st = trunctol() + + mps1 = FiniteMPS(M, ℂ^2, ℂ^D) + mps2 = FiniteMPS(M, Z2Space(0=>1, 1=>1), Z2Space(0=>D, 1=>D)) + mps3 = FiniteMPS(M, Vect[FermionParity](0=>1, 1=>1), Vect[FermionParity](0=>D, 1=>D)) + + H1 = transverse_field_ising_trivial_finite(M) + H2 = transverse_field_ising_z2_finite(M) + H3 = transverse_field_ising_fermion_parity_finite(M) + + # Find groundstate of non-compressed Hamiltonian + find_groundstate!(mps1, H1, DMRG2(; maxiter=100, trscheme=trunc_st)) + E0_1 = expectation_value(mps1, H1) + find_groundstate!(mps2, H2, DMRG2(; maxiter=100, trscheme=trunc_st)) + E0_2 = expectation_value(mps2, H2) + find_groundstate!(mps3, H3, DMRG2(; maxiter=100, trscheme=trunc_st)) + E0_3 = expectation_value(mps3, H3) + + # Compress MPO + mps1c = FiniteMPS(M, ℂ^2, ℂ^D) + mps2c = FiniteMPS(M, Z2Space(0=>1, 1=>1), Z2Space(0=>D, 1=>D)) + mps3c = FiniteMPS(M, Vect[FermionParity](0=>1, 1=>1), Vect[FermionParity](0=>D, 1=>D)) + + H1c, _ = mpo_compression(H1, 10^-10) + H2c, _ = mpo_compression(H2, 10^-10) + H3c, _ = mpo_compression(H3, 10^-10) + + # Find groundstate of compressed Hamiltonian + find_groundstate!(mps1c, H1c, DMRG2(; maxiter=100, trscheme=trunc_st)) + E0_1c = expectation_value(mps1c, H1c) + find_groundstate!(mps2c, H2c, DMRG2(; maxiter=100, trscheme=trunc_st)) + E0_2c = expectation_value(mps2c, H2c) + find_groundstate!(mps3c, H3c, DMRG2(; maxiter=100, trscheme=trunc_st)) + E0_3c = expectation_value(mps3c, H3c) + + # Assert that the ground state energies are equal up to a precision, display the + # compression + @assert E0_1 ≈ E0_1c + @assert E0_2 ≈ E0_2c + @assert E0_3 ≈ E0_3c + + println("Max dim uncompressed H1: $(max_dim_ham(H1)) || Max dim compressed H1c: + $(max_dim_ham(H1c))" ) + println("Max dim uncompressed H2: $(max_dim_ham(H2)) || Max dim compressed H2c: + $(max_dim_ham(H2c))" ) + println("Max dim uncompressed H3: $(max_dim_ham(H3)) || Max dim compressed H3c: + $(max_dim_ham(H3c))" ) + +end + +@testset "finite_mpo_compression_kitaev_model_energy" begin + M = 10 + D = 10 # Max bond dimension + trunc_st = trunctol() + + mps1 = FiniteMPS(M, Vect[FermionParity](0=>1, 1=>1), Vect[FermionParity](0=>D, 1=>D)) + + H1 = kitaev_model_finite(M) + + # Find groundstate of non-compressed Hamiltonian + find_groundstate!(mps1, H1, DMRG2(; maxiter=100, trscheme=trunc_st)) + E0_1 = expectation_value(mps1, H1) + + # Compress MPO + mps1c = FiniteMPS(M, Vect[FermionParity](0=>1, 1=>1), Vect[FermionParity](0=>D, 1=>D)) + + H1c, _ = mpo_compression(H1, 10^-10) + + # Find groundstate of compressed Hamiltonian + find_groundstate!(mps1c, H1c, DMRG2(; maxiter=100, trscheme=trunc_st)) + E0_1c = expectation_value(mps1c, H1c) + + # Assert that the ground state energies are equal up to a precision, display the + # compression + @assert E0_1 ≈ E0_1c + + println("Max dim uncompressed H1: $(max_dim_ham(H1)) || Max dim compressed H1c: + $(max_dim_ham(H1c))" ) + +end + +@testset "finite_mpo_compression_heisenberg_XXX_energy" begin + M = 10 + D = 6 # Max bond dimension + trunc_st = trunctol() + + mps1 = FiniteMPS(M, ℂ^3, ℂ^D) + mps2 = FiniteMPS(M, U1Space(0=>1, 1=>1, -1=>1), U1Space(0=>D, 1=>D, -1=>D)) + mps3 = FiniteMPS(M, SU2Space(1 => 1), SU2Space(1=>D)) + + H1 = heisenberg_XXX_trivial_finite(M) + H2 = heisenberg_XXX_U1_finite(M) + H3 = heisenberg_XXX_SU2_finite(M) + + # Find groundstate of non-compressed Hamiltonian + find_groundstate!(mps1, H1, DMRG2(; maxiter=100, trscheme=trunc_st)) + E0_1 = expectation_value(mps1, H1) + find_groundstate!(mps2, H2, DMRG2(; maxiter=100, trscheme=trunc_st)) + E0_2 = expectation_value(mps2, H2) + find_groundstate!(mps3, H3, DMRG2(; maxiter=100, trscheme=trunc_st)) + E0_3 = expectation_value(mps3, H3) + + # Compress MPO + mps1c = FiniteMPS(M, ℂ^3, ℂ^D) + mps2c = FiniteMPS(M, U1Space(0=>1, 1=>1, -1=>1), U1Space(0=>D, 1=>D, -1=>D)) + mps3c = FiniteMPS(M, SU2Space(1 => 1), SU2Space(1=>D)) + + H1c, _ = mpo_compression(H1, 10^-10) + H2c, _ = mpo_compression(H2, 10^-10) + H3c, _ = mpo_compression(H3, 10^-10) + + # Find groundstate of compressed Hamiltonian + find_groundstate!(mps1c, H1c, DMRG2(; maxiter=100, trscheme=trunc_st)) + E0_1c = expectation_value(mps1c, H1c) + find_groundstate!(mps2c, H2c, DMRG2(; maxiter=100, trscheme=trunc_st)) + E0_2c = expectation_value(mps2c, H2c) + find_groundstate!(mps3c, H3c, DMRG2(; maxiter=100, trscheme=trunc_st)) + E0_3c = expectation_value(mps3c, H3c) + + # Assert that the ground state energies are equal up to a precision, display the + # compression + @assert E0_1 ≈ E0_1c + @assert E0_2 ≈ E0_2c + @assert E0_3 ≈ E0_3c + + println("Max dim uncompressed H1: $(max_dim_ham(H1)) || Max dim compressed H1c: + $(max_dim_ham(H1c))" ) + println("Max dim uncompressed H2: $(max_dim_ham(H2)) || Max dim compressed H2c: + $(max_dim_ham(H2c))" ) + println("Max dim uncompressed H3: $(max_dim_ham(H3)) || Max dim compressed H3c: + $(max_dim_ham(H3c))" ) + +end + +@testset "finite_mpo_compression_bilinear_biquadratic_model_energy" begin + M = 10 + D = 6 # Max bond dimension + trunc_st = trunctol() + + mps1 = FiniteMPS(M, ℂ^3, ℂ^D) + mps2 = FiniteMPS(M, U1Space(0=>1, 1=>1, -1=>1), U1Space(0=>D, 1=>D, -1=>D)) + mps3 = FiniteMPS(M, SU2Space(1 => 1), SU2Space(1=>D)) + + H1 = bilinear_biquadratic_model_trivial_finite(M) + H2 = bilinear_biquadratic_model_U1_finite(M) + H3 = bilinear_biquadratic_model_SU2_finite(M) + + # Find groundstate of non-compressed Hamiltonian + find_groundstate!(mps1, H1, DMRG2(; maxiter=100, trscheme=trunc_st)) + E0_1 = expectation_value(mps1, H1) + find_groundstate!(mps2, H2, DMRG2(; maxiter=100, trscheme=trunc_st)) + E0_2 = expectation_value(mps2, H2) + find_groundstate!(mps3, H3, DMRG2(; maxiter=100, trscheme=trunc_st)) + E0_3 = expectation_value(mps3, H3) + + # Compress MPO + mps1c = FiniteMPS(M, ℂ^3, ℂ^D) + mps2c = FiniteMPS(M, U1Space(0=>1, 1=>1, -1=>1), U1Space(0=>D, 1=>D, -1=>D)) + mps3c = FiniteMPS(M, SU2Space(1 => 1), SU2Space(1=>D)) + + H1c, _ = mpo_compression(H1, 10^-10) + H2c, _ = mpo_compression(H2, 10^-10) + H3c, _ = mpo_compression(H3, 10^-10) + + # Find groundstate of compressed Hamiltonian + find_groundstate!(mps1c, H1c, DMRG2(; maxiter=100, trscheme=trunc_st)) + E0_1c = expectation_value(mps1c, H1c) + find_groundstate!(mps2c, H2c, DMRG2(; maxiter=100, trscheme=trunc_st)) + E0_2c = expectation_value(mps2c, H2c) + find_groundstate!(mps3c, H3c, DMRG2(; maxiter=100, trscheme=trunc_st)) + E0_3c = expectation_value(mps3c, H3c) + + # Assert that the ground state energies are equal up to a precision, display the + # compression + @assert E0_1 ≈ E0_1c + @assert E0_2 ≈ E0_2c + @assert E0_3 ≈ E0_3c + + println("Max dim uncompressed H1: $(max_dim_ham(H1)) || Max dim compressed H1c: + $(max_dim_ham(H1c))" ) + println("Max dim uncompressed H2: $(max_dim_ham(H2)) || Max dim compressed H2c: + $(max_dim_ham(H2c))" ) + println("Max dim uncompressed H3: $(max_dim_ham(H3)) || Max dim compressed H3c: + $(max_dim_ham(H3c))" ) + +end + +@testset "finite_mpo_compression_quantum_potts_model_energy" begin + M = 10 + D = 6 # Max bond dimension + trunc_st = trunctol() + + mps1 = FiniteMPS(M, ℂ^3, ℂ^D) + mps2 = FiniteMPS(M, Z3Space(0=>1, 1=>1, 2=>1), Z3Space(0=>1, 1=>1, 2=>1)) + + H1 = quantum_potts_trivial_finite(M) + H2 = quantum_potts_Z3_finite(M) + + # Find groundstate of non-compressed Hamiltonian + find_groundstate!(mps1, H1, DMRG2(; maxiter=100, trscheme=trunc_st)) + E0_1 = expectation_value(mps1, H1) + find_groundstate!(mps2, H2, DMRG2(; maxiter=100, trscheme=trunc_st)) + E0_2 = expectation_value(mps2, H2) + + # Compress MPO + mps1c = FiniteMPS(M, ℂ^3, ℂ^D) + mps2c = FiniteMPS(M, Z3Space(0=>1, 1=>1, 2=>1), Z3Space(0=>1, 1=>1, 2=>1)) + + H1c, _ = mpo_compression(H1, 10^-10) + H2c, _ = mpo_compression(H2, 10^-10) + + # Find groundstate of compressed Hamiltonian + find_groundstate!(mps1c, H1c, DMRG2(; maxiter=100, trscheme=trunc_st)) + E0_1c = expectation_value(mps1c, H1c) + find_groundstate!(mps2c, H2c, DMRG2(; maxiter=100, trscheme=trunc_st)) + E0_2c = expectation_value(mps2c, H2c) + + # Assert that the ground state energies are equal up to a precision, display the + # compression + @assert E0_1 ≈ E0_1c + @assert E0_2 ≈ E0_2c + + println("Max dim uncompressed H1: $(max_dim_ham(H1)) || Max dim compressed H1c: + $(max_dim_ham(H1c))" ) + println("Max dim uncompressed H2: $(max_dim_ham(H2)) || Max dim compressed H2c: + $(max_dim_ham(H2c))" ) + +end + +# Test random finite Hamiltonian +@testset "finite_random_long_range_ham" begin + # test the finite compression for a long range random Ising model + # check if the energy after the compression is the same as before + Random.seed!() + D = 20 # max bond dimensions + L = 10 # number of sites + range = 4 # cutoff for the maximum range between interactions + trunc_st = notrunc() + + mps = FiniteMPS(L, Z2Space(0=>1, 1=>1), Z2Space(0=>D, 1=>D)) + H = create_long_range_ising_symmetries_random_finite(L, range, 0.3, 0.4) + find_groundstate!(mps, H, DMRG2(; maxiter=100, trscheme=trunc_st)) + E0 = expectation_value(mps, H) + #println(" = $real(E0)") + + mps2 = FiniteMPS(L, Z2Space(0=>1, 1=>1), Z2Space(0=>D, 1=>D)) + H2, Rs = mpo_compression(H) + find_groundstate!(mps2, H2, DMRG2(; maxiter=100, trscheme=trunc_st)) + E1 = expectation_value(mps2, H2) +end \ No newline at end of file diff --git a/test/operators/infinite_mpo_compression.jl b/test/operators/infinite_mpo_compression.jl new file mode 100644 index 000000000..caebc96c6 --- /dev/null +++ b/test/operators/infinite_mpo_compression.jl @@ -0,0 +1,388 @@ +using Test, Random +using TensorKit +using MPSKit +using MPSKitModels + +function max_dim_ham(H) + N = length(H) + max_dim = 1 + for i in 1:N + if max_dim < dims(H[i])[4] + max_dim = dims(H[i])[4] + end + end + return max_dim +end + +# Utility functions to create Hamiltonians for testing + +function create_long_range_ising_symmetries_infinite(k) + infinite_chain = PeriodicVector([Z2Space(0=>1, 1=>1)]) + ZZ = S_zz(Z2Irrep) + X = S_x(Z2Irrep) + two_site_operators = [(1,j) => -ZZ for j in 1:k] + H = InfiniteMPOHamiltonian(infinite_chain, two_site_operators...) + return H +end + +function create_long_range_ising_symmetries_infinite_msite_random(M, k, σⱼ, σₕ) + V = Z2Space(0=>1, 1=>1) + infinite_chain = PeriodicVector(fill(V, M)) + ZZ = S_zz(Z2Irrep) + X = S_x(Z2Irrep) + single_site_operators = [i => (σₕ * randn()) * X for i in 1:M] + two_site_operators = [(i,i+j) => (1 + σⱼ * randn()) * ZZ for i in 1:M for j in 1:k] + H = InfiniteMPOHamiltonian(infinite_chain, single_site_operators..., + two_site_operators...) + return H +end + +function create_long_range_ising_symmetries_infinite_twosite_diff_int_strength(k, J, h) + infinite_chain = PeriodicVector([Z2Space(0=>1, 1=>1), Z2Space(0=>1, 1=>1)]) + ZZ = S_zz(Z2Irrep) + X = S_x(Z2Irrep) + two_site_operators_J12 = [(a,a+j) => J[1] * ZZ for a in [1 2] for j in 1:2:2k] + two_site_operators_J11 = [(1,j) => J[2] * ZZ for j in 3:2:2k] + two_site_operators_J22 = [(2,j) => J[3] * ZZ for j in 4:2:2k] + single_site_operators = [1 => h[1] * X, 2 => h[2] * X] + H = InfiniteMPOHamiltonian(infinite_chain, single_site_operators..., + two_site_operators_J11..., two_site_operators_J12..., + two_site_operators_J22...) + return H +end + +function transverse_field_ising_fermion_parity_infinite(M) + chain = InfiniteChain(M) + return transverse_field_ising(ComplexF64, FermionParity, chain) +end + +function transverse_field_ising_z2_infinite(M) + chain = InfiniteChain(M) + return transverse_field_ising(ComplexF64, Z2Irrep, chain) +end + +function transverse_field_ising_trivial_infinite(M) + chain = InfiniteChain(M) + return transverse_field_ising(ComplexF64, Trivial, chain) +end + +function kitaev_model_infinite(M) + chain = InfiniteChain(M) + return kitaev_model(ComplexF64, chain) +end + +function heisenberg_XXX_trivial_infinite(M) + chain = InfiniteChain(M) + return heisenberg_XXX(Trivial, chain) +end + +function heisenberg_XXX_U1_infinite(M) + chain = InfiniteChain(M) + return heisenberg_XXX(U1Irrep, chain) +end + +function heisenberg_XXX_SU2_infinite(M) + chain = InfiniteChain(M) + return heisenberg_XXX(SU2Irrep, chain) +end + +function bilinear_biquadratic_model_trivial_infinite(M) + chain = InfiniteChain(M) + return bilinear_biquadratic_model(ComplexF64, Trivial, chain) +end + +function quantum_potts_trivial_infinite(M) + chain = InfiniteChain(M) + return quantum_potts(ComplexF64, Trivial, chain) +end + +function quantum_potts_Z3_infinite(M) + chain = InfiniteChain(M) + return quantum_potts(ComplexF64, Z3Irrep, chain) +end + +@testset "infinite_mpo_compression_energy" begin + # test the infinite compression for a long range Ising model + # check if the energy after the compression is the same as before + D = 4 # max bond dimensions + range = 2 # cutoff for the maximum range between interactions + trunc_st = notrunc() + + mps = InfiniteMPS(Z2Space(0=>1, 1=>1), Z2Space(0=>D, 1=>D)) + H = create_long_range_ising_symmetries_infinite(range) + mps, = find_groundstate(mps, H, VUMPS(;maxiter=10)) + E0 = expectation_value(mps, H) + #println(" = $real(E0)") + + # mps2 = InfiniteMPS(Z2Space(0=>1, 1=>1), Z2Space(0=>D, 1=>D)) + # H2, Rs = mpo_finite_compression(H) + # find_groundstate!(mps2, H2, DMRG2(; maxiter=100, trscheme=trunc_st)) + # E1 = expectation_value(mps2, H2) + mps2 = InfiniteMPS(Z2Space(0=>1, 1=>1), Z2Space(0=>D, 1=>D)) + Q, P = mpo_compression(H) + H2 = Q + mps2, = find_groundstate(mps2, H2, VUMPS(;maxiter=100)) + E0 = expectation_value(mps2, H2) +end + +@testset "infinite_mpo_compression_energy_two_site" begin + # test the infinite compression for a long range Ising model + # check if the energy after the compression is the same as before + D = 10 # max bond dimensions + range = 4 # cutoff for the maximum range between interactions + trunc_st = notrunc() + J = [1, 1, 1] + h = [0.001, 0.001] + + spacephys = fill(Z2Space(0=>1, 1=>1), 2) + space_virt = fill(Z2Space(0=>D, 1=>D), 2) + + mps = InfiniteMPS(spacephys, space_virt) + H = create_long_range_ising_symmetries_infinite_twosite_diff_int_strength(range, J, h) + mps, = find_groundstate(mps, H, VUMPS(;maxiter=100)) + E0 = expectation_value(mps, H) + #println(" = $real(E0)") + + # mps2 = InfiniteMPS(Z2Space(0=>1, 1=>1), Z2Space(0=>D, 1=>D)) + # H2, Rs = mpo_finite_compression(H) + # find_groundstate!(mps2, H2, DMRG2(; maxiter=100, trscheme=trunc_st)) + # E1 = expectation_value(mps2, H2) + mps2 = InfiniteMPS(spacephys, space_virt) + Qs, Ps = mpo_compression(H, 0) + H2 = Qs + mps2, = find_groundstate(mps2, H2, VUMPS(;maxiter=100)) + E0 = expectation_value(mps2, H2) +end + +@testset "infinite_mpo_compression_transverse_ising_model_energy" begin + M = 4 + D = 10 # Max bond dimension + + mps1 = InfiniteMPS(fill(ℂ^2, M), fill(ℂ^D, M)) + mps2 = InfiniteMPS(fill(Z2Space(0=>1, 1=>1), M), fill(Z2Space(0=>D, 1=>D), M)) + mps3 = InfiniteMPS(fill(Vect[FermionParity](0=>1, 1=>1), M), fill(Vect[FermionParity]( + 0=>D, 1=>D), M)) + + H1 = transverse_field_ising_trivial_infinite(M) + H2 = transverse_field_ising_z2_infinite(M) + H3 = transverse_field_ising_fermion_parity_infinite(M) + + # Find groundstate of non-compressed Hamiltonian + mps1, = find_groundstate(mps1, H1, VUMPS(; maxiter=100)) + E0_1 = expectation_value(mps1, H1) + mps2, = find_groundstate(mps2, H2, VUMPS(; maxiter=100)) + E0_2 = expectation_value(mps2, H2) + mps3, = find_groundstate(mps3, H3, VUMPS(; maxiter=100)) + E0_3 = expectation_value(mps3, H3) + + # Compress MPO + mps1c = InfiniteMPS(fill(ℂ^2, M), fill(ℂ^D, M)) + mps2c = InfiniteMPS(fill(Z2Space(0=>1, 1=>1), M), fill(Z2Space(0=>D, 1=>D), M)) + mps3c = InfiniteMPS(fill(Vect[FermionParity](0=>1, 1=>1), M), fill(Vect[FermionParity]( + 0=>D, 1=>D), M)) + + H1c, _ = mpo_compression(H1, 10^-10) + H2c, _ = mpo_compression(H2, 10^-10) + H3c, _ = mpo_compression(H3, 10^-10) + + # Find groundstate of compressed Hamiltonian + mps1c, = find_groundstate(mps1c, H1c, VUMPS(; maxiter=100)) + E0_1c = expectation_value(mps1c, H1c) + mps2c, = find_groundstate(mps2c, H2c, VUMPS(; maxiter=100)) + E0_2c = expectation_value(mps2c, H2c) + mps3c, = find_groundstate(mps3c, H3c, VUMPS(; maxiter=100)) + E0_3c = expectation_value(mps3c, H3c) + + # Assert that the ground state energies are equal up to a precision, display the + # compression + @assert E0_1 ≈ E0_1c + @assert E0_2 ≈ E0_2c + @assert E0_3 ≈ E0_3c + + println("Max dim uncompressed H1: $(max_dim_ham(H1)) || Max dim compressed H1c: + $(max_dim_ham(H1c))" ) + println("Max dim uncompressed H2: $(max_dim_ham(H2)) || Max dim compressed H2c: + $(max_dim_ham(H2c))" ) + println("Max dim uncompressed H3: $(max_dim_ham(H3)) || Max dim compressed H3c: + $(max_dim_ham(H3c))" ) + +end + +@testset "infinite_mpo_compression_kitaev_model_energy" begin + M = 10 + D = 10 # Max bond dimension + + mps1 = InfiniteMPS(fill(Vect[FermionParity](0=>1, 1=>1), M), + fill(Vect[FermionParity](0=>D, 1=>D), M)) + + H1 = kitaev_model_infinite(M) + + # Find groundstate of non-compressed Hamiltonian + mps1, = find_groundstate(mps1, H1, VUMPS(; maxiter=100)) + E0_1 = expectation_value(mps1, H1) + + # Compress MPO + mps1c = InfiniteMPS(fill(Vect[FermionParity](0=>1, 1=>1), M), + fill(Vect[FermionParity](0=>D, 1=>D), M)) + + H1c, _ = mpo_compression(H1, 10^-10) + + # Find groundstate of compressed Hamiltonian + mps1c, = find_groundstate(mps1c, H1c, VUMPS(; maxiter=100)) + E0_1c = expectation_value(mps1c, H1c) + + # Assert that the ground state energies are equal up to a precision, display the + # compression + @assert E0_1 ≈ E0_1c + + println("Max dim uncompressed H1: $(max_dim_ham(H1)) || Max dim compressed H1c: + $(max_dim_ham(H1c))" ) + +end + +@testset "infinite_mpo_compression_heisenberg_XXX_energy" begin + M = 10 + D = 6 # Max bond dimension + + mps1 = InfiniteMPS(fill(ℂ^3, M),fill(ℂ^D, M)) + mps2 = InfiniteMPS(fill(U1Space(0=>1, 1=>1, -1=>1), M), fill(U1Space(0=>D, 1=>D, -1=>D), M)) + mps3 = InfiniteMPS(fill(SU2Space(1 => 1), M), fill(SU2Space(1=>D), M)) + + H1 = heisenberg_XXX_trivial_infinite(M) + H2 = heisenberg_XXX_U1_infinite(M) + H3 = heisenberg_XXX_SU2_infinite(M) + + # Find groundstate of non-compressed Hamiltonian + mps1, = find_groundstate(mps1, H1, VUMPS(; maxiter=100)) + E0_1 = expectation_value(mps1, H1) + mps2, = find_groundstate(mps2, H2, VUMPS(; maxiter=100)) + E0_2 = expectation_value(mps2, H2) + mps3, = find_groundstate(mps3, H3, VUMPS(; maxiter=100)) + E0_3 = expectation_value(mps3, H3) + + # Compress MPO + mps1c = InfiniteMPS(fill(ℂ^3, M),fill(ℂ^D, M)) + mps2c = InfiniteMPS(fill(U1Space(0=>1, 1=>1, -1=>1), M), fill(U1Space(0=>D, 1=>D, -1=>D), M)) + mps3c = InfiniteMPS(fill(SU2Space(1 => 1), M), fill(SU2Space(1=>D), M)) + + H1c, _ = mpo_compression(H1, 10^-10) + H2c, _ = mpo_compression(H2, 10^-10) + H3c, _ = mpo_compression(H3, 10^-10) + + # Find groundstate of compressed Hamiltonian + mps1c, = find_groundstate(mps1c, H1c, VUMPS(; maxiter=100)) + E0_1c = expectation_value(mps1c, H1c) + mps2c, = find_groundstate(mps2c, H2c, VUMPS(; maxiter=100)) + E0_2c = expectation_value(mps2c, H2c) + mps3c, = find_groundstate(mps3c, H3c, VUMPS(; maxiter=100)) + E0_3c = expectation_value(mps3c, H3c) + + # Assert that the ground state energies are equal up to a precision, display the + # compression + @assert E0_1 ≈ E0_1c + @assert E0_2 ≈ E0_2c + @assert E0_3 ≈ E0_3c + + println("Max dim uncompressed H1: $(max_dim_ham(H1)) || Max dim compressed H1c: + $(max_dim_ham(H1c))" ) + println("Max dim uncompressed H2: $(max_dim_ham(H2)) || Max dim compressed H2c: + $(max_dim_ham(H2c))" ) + println("Max dim uncompressed H3: $(max_dim_ham(H3)) || Max dim compressed H3c: + $(max_dim_ham(H3c))" ) + +end + +@testset "infinite_mpo_compression_bilinear_biquadratic_model_energy" begin + M = 10 + D = 6 # Max bond dimension + + mps1 = InfiniteMPS(fill(ℂ^3, M),fill(ℂ^D, M)) + + + H1 = bilinear_biquadratic_model_trivial_infinite(M) + + # Find groundstate of non-compressed Hamiltonian + mps1, = find_groundstate(mps1, H1, VUMPS(; maxiter=100)) + E0_1 = expectation_value(mps1, H1) + + # Compress MPO + mps1c = InfiniteMPS(fill(ℂ^3, M),fill(ℂ^D, M)) + + H1c, _ = mpo_compression(H1, 10^-10) + + # Find groundstate of compressed Hamiltonian + mps1c, = find_groundstate(mps1c, H1c, VUMPS(; maxiter=100)) + E0_1c = expectation_value(mps1c, H1c) + + # Assert that the ground state energies are equal up to a precision, display the + # compression + @assert E0_1 ≈ E0_1c + + println("Max dim uncompressed H1: $(max_dim_ham(H1)) || Max dim compressed H1c: + $(max_dim_ham(H1c))" ) + +end + +@testset "infinite_mpo_compression_quantum_potts_model_energy" begin + M = 10 + D = 6 # Max bond dimension + + mps1 = InfiniteMPS(fill(ℂ^3, M), fill(ℂ^D, M)) + mps2 = InfiniteMPS(fill(Z3Space(0=>1, 1=>1, 2=>1), M), fill(Z3Space(0=>1, 1=>1, 2=>1), M)) + + H1 = quantum_potts_trivial_infinite(M) + H2 = quantum_potts_Z3_infinite(M) + + # Find groundstate of non-compressed Hamiltonian + mps1, = find_groundstate(mps1, H1, VUMPS(; maxiter=100)) + E0_1 = expectation_value(mps1, H1) + mps2, = find_groundstate(mps2, H2, VUMPS(; maxiter=100)) + E0_2 = expectation_value(mps2, H2) + + # Compress MPO + mps1c = InfiniteMPS(fill(ℂ^3, M), fill(ℂ^D, M)) + mps2c = InfiniteMPS(fill(Z3Space(0=>1, 1=>1, 2=>1), M), fill(Z3Space(0=>1, 1=>1, 2=>1), M)) + + H1c, _ = mpo_compression(H1, 10^-10) + H2c, _ = mpo_compression(H2, 10^-10) + + # Find groundstate of compressed Hamiltonian + mps1c, = find_groundstate(mps1c, H1c, VUMPS(; maxiter=100)) + E0_1c = expectation_value(mps1c, H1c) + mps2c, = find_groundstate(mps2c, H2c, VUMPS(; maxiter=100)) + E0_2c = expectation_value(mps2c, H2c) + + # Assert that the ground state energies are equal up to a precision, display the + # compression + @assert E0_1 ≈ E0_1c + @assert E0_2 ≈ E0_2c + + println("Max dim uncompressed H1: $(max_dim_ham(H1)) || Max dim compressed H1c: + $(max_dim_ham(H1c))" ) + println("Max dim uncompressed H2: $(max_dim_ham(H2)) || Max dim compressed H2c: + $(max_dim_ham(H2c))" ) + +end + +# Test random infinite Hamiltonian +@testset "infinite_random_long_range_ham" begin + # test the finite compression for a long range random Ising model + # check if the energy after the compression is the same as before + Random.seed!() + D = 20 # max bond dimensions + L = 10 # number of sites + range = 4 # cutoff for the maximum range between interactions + + mps = InfiniteMPS(fill(Z2Space(0=>1, 1=>1),L), fill(Z2Space(0=>D, 1=>D),L)) + H = create_long_range_ising_symmetries_infinite_msite_random(L, range, 0.1, 0.001) + mps, = find_groundstate(mps, H, VUMPS(; maxiter=100)) + E0 = expectation_value(mps, H) + #println(" = $real(E0)") + + mps2 = InfiniteMPS(fill(Z2Space(0=>1, 1=>1),L), fill(Z2Space(0=>D, 1=>D),L)) + H2, Rs = mpo_compression(H, 0) + mps2, = find_groundstate(mps2, H2, VUMPS(; maxiter=100)) + E1 = expectation_value(mps2, H2) + + @assert E0 ≈ E1 +end \ No newline at end of file From 9da38cf15ed12720af33eef1d32fcb268871e983 Mon Sep 17 00:00:00 2001 From: herviou Date: Wed, 15 Jul 2026 13:53:44 +0200 Subject: [PATCH 2/3] Fixing format, add some basic verbosity control, correct a few names --- src/operators/mpocompression.jl | 470 ++++++++++++++++---------------- 1 file changed, 236 insertions(+), 234 deletions(-) diff --git a/src/operators/mpocompression.jl b/src/operators/mpocompression.jl index 9691abf3f..3867702d1 100644 --- a/src/operators/mpocompression.jl +++ b/src/operators/mpocompression.jl @@ -1,6 +1,6 @@ """ - mpo_compression(H::FiniteMPOHamiltonian; η::Number=10^-8) -> Hᵪ, Rs - mpo_compression(H::InfiniteMPOHamiltonian; η::Number=10^-8) -> Qs, Ps + mpo_compression(H::FiniteMPOHamiltonian; truncbelow::Number=10^-8) -> Hᵪ, Rs + mpo_compression(H::InfiniteMPOHamiltonian; truncbelow::Number=10^-8) -> Qs, Ps Returns the compressed version of the FiniteMPOHamiltonian or InfiniteMPOHamiltonian, as described in: https://arxiv.org/pdf/1909.06341. @@ -9,22 +9,22 @@ as described in: https://arxiv.org/pdf/1909.06341. - `H`: Hamiltonian to be compressed. Either FiniteMPOHamiltonian or InfiniteMPOHamiltonian. ### Keyword arguments -- `η::Number=10^-8`: precision of the compression. The values smaller than η in the SVD +- `truncbelow::Number=10^-8`: precision of the compression. The values smaller than truncbelow in the SVD decomposition during compression will be disregarded. ### Returns -- mpo_compression(H::FiniteMPOHamiltonian; η::Number=10^-8) +- mpo_compression(H::FiniteMPOHamiltonian; truncbelow::Number=10^-8) `Hᵪ::FiniteMPOHamiltonian`: compressed Hamiltonian in the left canonical form. `Rs::Vector{BlockTensorMap}`: matrices that can be applied to change between left and right canonical forms of Hᵪ. -- mpo_compression(H::InfiniteMPOHamiltonian; η::Number=10^-8) +- mpo_compression(H::InfiniteMPOHamiltonian; truncbelow::Number=10^-8) `Qs::InfiniteMPOHamiltonian`: compressed Hamiltonian in the left canonical form. `Ps::InfiniteMPOHamiltonian`: compressed Hamiltonian in the right canonical form. """ function mpo_compression end -function mpo_compression(H::FiniteMPOHamiltonian, η::Number=10^-8) +function mpo_compression(H::FiniteMPOHamiltonian; truncbelow::Number = 10^-8, verbosity::Int = 0) t = @elapsed begin # make sure that the MPO have at least three sites N = length(H) @@ -32,7 +32,8 @@ function mpo_compression(H::FiniteMPOHamiltonian, η::Number=10^-8) # change MPO from N x 1 x 1 x M blocks to 3 x 1 x 1 x 3 block structure TT = TensorMap{scalartype(H), spacetype(H[1]), 2, 2, Vector{scalartype(H)}} Hnew = Vector{SparseBlockTensorMap{TT, scalartype(H), spacetype(H[1]), 2, 2, 4}}( - undef, N) + undef, N + ) for n in 1:N Hnew[n] = SparseBlockTensorMap(reduce_blocks_mpo(H[n])) end @@ -41,21 +42,21 @@ function mpo_compression(H::FiniteMPOHamiltonian, η::Number=10^-8) TT = AbstractTensorMap{scalartype(H), spacetype(H), 2, 2} Cs = Vector{SparseBlockTensorMap{TT, scalartype(H), spacetype(H), 2, 2, 4}}(undef, N) Cs[1] = Qs[1] - Rs = Vector{AbstractBlockTensorMap}(undef, N-2) + Rs = Vector{AbstractBlockTensorMap}(undef, N - 2) R = id(codomain(Qs[2])[1]) - for n in 2:N-1 - @tensor Qnew[a,i;j,b] := R[a,c] * Qs[n][c,i;j,b] + for n in 2:(N - 1) + @tensor Qnew[a, i; j, b] := R[a, c] * Qs[n][c, i; j, b] Q, R = qr_block_respecting(Qnew) # For now it is safe to assume that R will be 3x3 blocked BlockTensorMap, # since it went through right_canonical_mpo_finite first M, R′ = decompose_R(R) - U, S, V′ = block_respecting_svd(M, η) - @tensor C[a,i;j,b] := Q[a,i;j,c] * U[c,b] - @tensor R[a;b] := S[a,c] * V′[c,d] * R′[d,b] + U, S, V′ = block_respecting_svd(M, truncbelow) + @tensor C[a, i; j, b] := Q[a, i; j, c] * U[c, b] + @tensor R[a;b] := S[a, c] * V′[c, d] * R′[d, b] Cs[n] = C - Rs[n-1] = R + Rs[n - 1] = R end - @tensor Cnew[a,i;j,b] := R[a,c] * Qs[N][c,i;j,b] + @tensor Cnew[a, i; j, b] := R[a, c] * Qs[N][c, i; j, b] Cs[N] = Cnew # Last step: Make new Cs into FiniteMPOHamiltonian # TA = AbstractTensorMap{scalartype(H), spacetype(H), 2, 2} @@ -66,59 +67,61 @@ function mpo_compression(H::FiniteMPOHamiltonian, η::Number=10^-8) Hᵪ = Vector{Hᵪtype}(undef, N) for n in 1:N W = SparseBlockTensorMap(Cs[n]) - A = W[2:(end-1), 1, 1, 2:(end-1)] - B = removeunit(W[2:(end-1), 1, 1, end], 4) - C = removeunit(W[1,1,1,2:(end-1)], 1) - D = removeunit(removeunit(W[1,1,1,end:end], 4), 1) + A = W[2:(end - 1), 1, 1, 2:(end - 1)] + B = removeunit(W[2:(end - 1), 1, 1, end], 4) + C = removeunit(W[1, 1, 1, 2:(end - 1)], 1) + D = removeunit(removeunit(W[1, 1, 1, end:end], 4), 1) W = JordanMPOTensor(space(W), A, B, C, D) Hᵪ[n] = W end Hᵪ = FiniteMPOHamiltonian(Hᵪ) end # Print reduction data - tot_dim_original = total_virt_dimension(H) - tot_dim_compressed = total_virt_dimension(Hᵪ) - printstyled("┌───────────────────────────────────────\n", color=:cyan) - printstyled("| MPO Compression performed succesfully\n", color=:cyan) - printstyled("│ ⏱ Time: ", color=:cyan) - printstyled("$(round(t, digits=1)) s\n", color=:yellow) - printstyled("| Initial total virtual dimensions: ", color=:cyan) - printstyled("$(tot_dim_original)\n", color=:yellow) - printstyled("| Final total virtual dimensions: ", color=:cyan) - printstyled("$(tot_dim_compressed)\n", color=:yellow) - printstyled("└───────────────────────────────────────\n") - + if verbosity > 0 + tot_dim_original = total_virt_dimension(H) + tot_dim_compressed = total_virt_dimension(Hᵪ) + printstyled("┌───────────────────────────────────────\n", color = :cyan) + printstyled("| MPO Compression performed succesfully\n", color = :cyan) + printstyled("│ ⏱ Time: ", color = :cyan) + printstyled("$(round(t, digits = 1)) s\n", color = :yellow) + printstyled("| Initial total virtual dimensions: ", color = :cyan) + printstyled("$(tot_dim_original)\n", color = :yellow) + printstyled("| Final total virtual dimensions: ", color = :cyan) + printstyled("$(tot_dim_compressed)\n", color = :yellow) + printstyled("└───────────────────────────────────────\n") + end return Hᵪ, Rs end -function mpo_compression(H::InfiniteMPOHamiltonian, η::Number=10^-8) +function mpo_compression(H::InfiniteMPOHamiltonian; truncbelow::Number = 10^-8, verbosity::Int = 0, canonical_precision = 1.0e-10) tot_dim_original = total_virt_dimension(H) t = @elapsed begin N = length(H) # Number of sites # change MPO from N x 1 x 1 x M blocks to 3 x 1 x 1 x 3 block structure TT = TensorMap{scalartype(H), spacetype(H[1]), 2, 2, Vector{scalartype(H)}} Hnew = Vector{SparseBlockTensorMap{TT, scalartype(H), spacetype(H[1]), 2, 2, 4}}( - undef, N) + undef, N + ) for n in 1:N Hnew[n] = SparseBlockTensorMap(reduce_blocks_mpo(H[n])) end - HL, _ = left_canonical_mpo_infinite_iter_msites(Hnew) - _, HR = right_canonical_mpo_infinite_iter_msites(HL) - HL, Cs = left_canonical_mpo_infinite_iter_msites(HR) + HL, _ = left_canonical_mpo_infinite_iter_msites(Hnew; precision = canonical_precision, verbosity) + _, HR = right_canonical_mpo_infinite_iter_msites(HL; precision = canonical_precision, verbosity) + HL, Cs = left_canonical_mpo_infinite_iter_msites(HR; precision = canonical_precision, verbosity) Us = Vector{BlockTensorMap}(undef, N) Ss = Vector{BlockTensorMap}(undef, N) Vs = Vector{BlockTensorMap}(undef, N) - for n = 1:N - Us[n], Ss[n], Vs[n] = block_respecting_svd(Cs[n], η) + for n in 1:N + Us[n], Ss[n], Vs[n] = block_respecting_svd(Cs[n]; truncbelow) end Qs = Vector{BlockTensorMap}(undef, N) Ps = Vector{BlockTensorMap}(undef, N) - for n = 1:N + for n in 1:N # Left canonical version of compressed iMPO - @tensor Q[a,i;j,b] := Us[mod1(n-1,N)]'[a,c] * HL[n][c,i;j,d] * Us[n][d,b] + @tensor Q[a, i; j, b] := Us[mod1(n - 1, N)]'[a, c] * HL[n][c, i; j, d] * Us[n][d, b] # Right canonical version of compressed iMPO - @tensor P[a,i;j,b] := Vs[mod1(n-1,N)][a,c] * HR[n][c,i;j,d] * Vs[n]'[d,b] + @tensor P[a, i; j, b] := Vs[mod1(n - 1, N)][a, c] * HR[n][c, i; j, d] * Vs[n]'[d, b] Qs[n] = Q Ps[n] = P end @@ -128,30 +131,31 @@ function mpo_compression(H::InfiniteMPOHamiltonian, η::Number=10^-8) end # Print reduction data - tot_dim_compressed = total_virt_dimension(Q̂s) - printstyled("┌───────────────────────────────────────\n", color=:cyan) - printstyled("| MPO Compression performed succesfully\n", color=:cyan) - printstyled("│ ⏱ Time: ", color=:cyan) - printstyled("$(round(t, digits=1)) s\n", color=:yellow) - printstyled("| Initial total virtual dimensions: ", color=:cyan) - printstyled("$(tot_dim_original)\n", color=:yellow) - printstyled("| Final total virtual dimensions: ", color=:cyan) - printstyled("$(tot_dim_compressed)\n", color=:yellow) - printstyled("└───────────────────────────────────────\n") - + if verbosity > 0 + tot_dim_compressed = total_virt_dimension(Q̂s) + printstyled("┌───────────────────────────────────────\n", color = :cyan) + printstyled("| MPO Compression performed succesfully\n", color = :cyan) + printstyled("│ ⏱ Time: ", color = :cyan) + printstyled("$(round(t, digits = 1)) s\n", color = :yellow) + printstyled("| Initial total virtual dimensions: ", color = :cyan) + printstyled("$(tot_dim_original)\n", color = :yellow) + printstyled("| Final total virtual dimensions: ", color = :cyan) + printstyled("$(tot_dim_compressed)\n", color = :yellow) + printstyled("└───────────────────────────────────────\n") + end return Q̂s, P̂s end # utility # ------- function trace_single_block(d::AbstractTensorMap) - return @tensor trd[a;b] := d[a,i;i,b] + return @tensor trd[a;b] := d[a, i; i, b] end function change_value_at_fusiontree!(T::AbstractTensorMap, i::Integer, v::Number) # change value of the fusiontree at index i of a TensorMap T with new value v f = fusiontrees(T) - if f isa Vector + return if f isa Vector T[f[i]...] .= v elseif f isa Indices T[f.values[i]...] .= v @@ -164,7 +168,7 @@ end function total_virt_dimension(H::FiniteMPOHamiltonian) # calculate the sum of bond dimensions for entire H MPO N = length(H) - tvd = sum([dims(H[i])[4] for i in 1:N-1]) + tvd = sum([dims(H[i])[4] for i in 1:(N - 1)]) return tvd end @@ -191,39 +195,49 @@ function qr_block_respecting(W::AbstractBlockTensorMap) # traceless id1 = id(codomain(W)[2]) - t = zeros(ComplexF64, domain(W)[2][1] ← domain(W)[2][2:end-1]) - for block in 2:blM-1 - d = W[1,1,1,block] + t = zeros(ComplexF64, domain(W)[2][1] ← domain(W)[2][2:(end - 1)]) + for block in 2:(blM - 1) + d = W[1, 1, 1, block] trd = trace_single_block(d) - @tensor W[1,1,1,block][a,i;j,b] = d[a,i;j,b] - trd[a;b] * id1[i;j] + @tensor W[1, 1, 1, block][a, i; j, b] = d[a, i; j, b] - trd[a;b] * id1[i;j] # Here goes the creation of matrix R' with a t' tensor - t[block-1] = trd / dim_st[2] + t[block - 1] = trd / dim_st[2] end # Perform QR of the 2:blM-1 block - V = W[1:blN-1, :, :, 2:blM-1] - V = permute(V, ((1,2,3), (4,))) - - V = convert(SparseBlockTensorMap{TensorMap{eltype(V).parameters...}, - eltype(storagetype(V)), spacetype(V), 3, 1, 4,}, V) + V = W[1:(blN - 1), :, :, 2:(blM - 1)] + V = permute(V, ((1, 2, 3), (4,))) + + V = convert( + SparseBlockTensorMap{ + TensorMap{eltype(V).parameters...}, + eltype(storagetype(V)), spacetype(V), 3, 1, 4, + }, V + ) Q, R = qr_compact(V) - Q = permute(Q, ((1,2),(3,4))) + Q = permute(Q, ((1, 2), (3, 4))) # Put Q into larger Q̂ - Q̂ = zeros(scalartype(W), codomain(W) ← domain(W)[1] ⊗ (domain(W)[2][1] ⊞ domain(Q)[2] ⊞ - domain(W)[2][end])) - Q̂[1,1,1,1] = Q̂[end,1,1,end] = W[1,1,1,1] - Q̂[1:end-1,1,1,2:end-1] = Q - Q̂[1:end-1,1,1,end] = W[1:end-1,1,1,end] + Q̂ = zeros( + scalartype(W), codomain(W) ← domain(W)[1] ⊗ ( + domain(W)[2][1] ⊞ domain(Q)[2] ⊞ + domain(W)[2][end] + ) + ) + Q̂[1, 1, 1, 1] = Q̂[end, 1, 1, end] = W[1, 1, 1, 1] + Q̂[1:(end - 1), 1, 1, 2:(end - 1)] = Q + Q̂[1:(end - 1), 1, 1, end] = W[1:(end - 1), 1, 1, end] # Put R and t into larger R̂ - R̂ = zeros(scalartype(W), domain(W)[2][1] ⊞ codomain(R)[1] ⊞ domain(W)[2][end] ← - domain(W)[2][1] ⊞ domain(R)[1] ⊞ domain(W)[2][end]) - R̂[2:end-1,2:end-1] = R - R̂[1,2:end-1] = t - change_value_at_fusiontree!(R̂[1,1],1,1) - change_value_at_fusiontree!(R̂[end,end],1,1) - + R̂ = zeros( + scalartype(W), domain(W)[2][1] ⊞ codomain(R)[1] ⊞ domain(W)[2][end] ← + domain(W)[2][1] ⊞ domain(R)[1] ⊞ domain(W)[2][end] + ) + R̂[2:(end - 1), 2:(end - 1)] = R + R̂[1, 2:(end - 1)] = t + change_value_at_fusiontree!(R̂[1, 1], 1, 1) + change_value_at_fusiontree!(R̂[end, end], 1, 1) + return Q̂, R̂ end @@ -240,47 +254,55 @@ function lq_block_respecting(W::AbstractBlockTensorMap) # traceless id1 = id(codomain(W)[2]) - t = zeros(ComplexF64, codomain(W)[1][2:end-1] ← codomain(W)[1][end]) - for block in 2:blN-1 - d = W[block,1,1,end] + t = zeros(ComplexF64, codomain(W)[1][2:(end - 1)] ← codomain(W)[1][end]) + for block in 2:(blN - 1) + d = W[block, 1, 1, end] trd = trace_single_block(d) - @tensor W[block,1,1,end][a,i;j,b] = d[a,i;j,b] - trd[a;b] * id1[i;j] + @tensor W[block, 1, 1, end][a, i; j, b] = d[a, i; j, b] - trd[a;b] * id1[i;j] # Here goes the creation of matrix R' with a t' tensor - t[block-1] = trd / dim_st[2] + t[block - 1] = trd / dim_st[2] end # Perform QR of the 2:blM-1 block (with virtual indices transposed) - V = W[2:blN-1, :, :, 2:blM] - V = permute(V, ((1,), (2,3,4))) - - V = convert(SparseBlockTensorMap{TensorMap{eltype(V).parameters...}, - eltype(storagetype(V)), spacetype(V), 1, 3, 4,}, V) + V = W[2:(blN - 1), :, :, 2:blM] + V = permute(V, ((1,), (2, 3, 4))) + + V = convert( + SparseBlockTensorMap{ + TensorMap{eltype(V).parameters...}, + eltype(storagetype(V)), spacetype(V), 1, 3, 4, + }, V + ) # to perform QR L, Q = lq_compact(V) # transpose vitrual indices back (now to reconstruct initial matrix # one has to apply R on the left) - Q = permute(Q, ((1,2),(3,4))) + Q = permute(Q, ((1, 2), (3, 4))) # Put Q into larger Q̂ - Q̂ = zeros(scalartype(W), (codomain(W)[1][1] ⊞ codomain(Q)[1] ⊞ codomain(W)[1][end]) ⊗ - codomain(W)[2] ← domain(W)) - Q̂[1,1,1,1] = Q̂[end,1,1,end] = W[1,1,1,1] - Q̂[2:end-1,1,1,2:end] = Q - Q̂[1,1,1,2:end] = W[1,1,1,2:end] + Q̂ = zeros( + scalartype(W), (codomain(W)[1][1] ⊞ codomain(Q)[1] ⊞ codomain(W)[1][end]) ⊗ + codomain(W)[2] ← domain(W) + ) + Q̂[1, 1, 1, 1] = Q̂[end, 1, 1, end] = W[1, 1, 1, 1] + Q̂[2:(end - 1), 1, 1, 2:end] = Q + Q̂[1, 1, 1, 2:end] = W[1, 1, 1, 2:end] # Put L and t into larger L̂ - L̂ = zeros(scalartype(W), domain(W)[2][1] ⊞ codomain(L)[1] ⊞ domain(W)[2][end] ← - domain(W)[2][1] ⊞ domain(L)[1] ⊞ domain(W)[2][end]) - L̂[2:end-1,2:end-1] = L - L̂[2:end-1,end] = t - change_value_at_fusiontree!(L̂[1,1],1,1) - change_value_at_fusiontree!(L̂[end,end],1,1) + L̂ = zeros( + scalartype(W), domain(W)[2][1] ⊞ codomain(L)[1] ⊞ domain(W)[2][end] ← + domain(W)[2][1] ⊞ domain(L)[1] ⊞ domain(W)[2][end] + ) + L̂[2:(end - 1), 2:(end - 1)] = L + L̂[2:(end - 1), end] = t + change_value_at_fusiontree!(L̂[1, 1], 1, 1) + change_value_at_fusiontree!(L̂[end, end], 1, 1) return L̂, Q̂ end function left_canonical_mpo_finite( H::Union{FiniteMPOHamiltonian, V} - ) where {T<:SparseBlockTensorMap{<:Any, <:Any, <:Any, 2, 2, 4}, V<:Vector{T}} + ) where {T <: SparseBlockTensorMap{<:Any, <:Any, <:Any, 2, 2, 4}, V <: Vector{T}} # change finite Hamiltonian MPO into its Left canonical form. The 'H' is a # FiniteMPOHamiltonian from MPSKit (l and r boundary conditions are already # part of the left-most and right-most MPOs, they can not be compressed, and @@ -295,22 +317,22 @@ function left_canonical_mpo_finite( Qs = Vector{SparseBlockTensorMap{TT, scalartype(H), spacetype(H[1]), 2, 2, 4}}(undef, N) Qs[1] = H[1] Qs[2] = Q - Rs = Vector{AbstractBlockTensorMap}(undef, N-2) + Rs = Vector{AbstractBlockTensorMap}(undef, N - 2) Rs[1] = R - for n in 3:N-1 - @tensor Wnew[a,i;j,b] := R[a,c] * H[n][c,i;j,b] + for n in 3:(N - 1) + @tensor Wnew[a, i; j, b] := R[a, c] * H[n][c, i; j, b] Q, R = qr_block_respecting(Wnew) Qs[n] = Q - Rs[n-1] = R + Rs[n - 1] = R end - @tensor Wlast[a,i;j,b] := R[a,c] * H[N][c,i;j,b] + @tensor Wlast[a, i; j, b] := R[a, c] * H[N][c, i; j, b] Qs[N] = Wlast return Qs, Rs end function right_canonical_mpo_finite( H::Union{FiniteMPOHamiltonian, V} - ) where {T<:SparseBlockTensorMap{<:Any, <:Any, <:Any, 2, 2, 4}, V<:Vector{T}} + ) where {T <: SparseBlockTensorMap{<:Any, <:Any, <:Any, 2, 2, 4}, V <: Vector{T}} # change finite Hamiltonian MPO into its Right canonical form. The 'H' is a # FiniteMPOHamiltonian from MPSKit (l and r boundary conditions are already # part of the left-most and right-most MPOs, they can not be compressed, and @@ -320,20 +342,20 @@ function right_canonical_mpo_finite( N = length(H) @assert N ≥ 3 - L, Q = lq_block_respecting(H[end-1]) + L, Q = lq_block_respecting(H[end - 1]) TT = AbstractTensorMap{scalartype(H), spacetype(H[1]), 2, 2} Qs = Vector{SparseBlockTensorMap{TT, scalartype(H), spacetype(H[1]), 2, 2, 4}}(undef, N) Qs[end] = H[end] - Qs[end-1] = Q - Ls = Vector{AbstractBlockTensorMap}(undef, N-2) + Qs[end - 1] = Q + Ls = Vector{AbstractBlockTensorMap}(undef, N - 2) Ls[end] = L - for n in N-2:-1:2 - @tensor Wnew[a,i;j,b] := H[n][a,i;j,c] * L[c,b] + for n in (N - 2):-1:2 + @tensor Wnew[a, i; j, b] := H[n][a, i; j, c] * L[c, b] L, Q = lq_block_respecting(Wnew) Qs[n] = Q - Ls[n-1] = L + Ls[n - 1] = L end - @tensor Wlast[a,i;j,b] := H[1][a,i;j,c] * L[c,b] + @tensor Wlast[a, i; j, b] := H[1][a, i; j, c] * L[c, b] Qs[1] = Wlast return Qs, Ls end @@ -341,53 +363,39 @@ end function decompose_R(R::AbstractBlockTensorMap) M = zeros(ComplexF64, codomain(R) ← domain(R)) R′ = zeros(ComplexF64, domain(R) ← domain(R)) - change_value_at_fusiontree!(M[1,1], 1, 1) - change_value_at_fusiontree!(M[end,end], 1, 1) - change_value_at_fusiontree!(R′[1,1], 1, 1) - change_value_at_fusiontree!(R′[end,end], 1, 1) - M[2:end-1,2:end-1] = R[2:end-1,2:end-1] - R′[1,2:end-1] = R[1,2:end-1] + change_value_at_fusiontree!(M[1, 1], 1, 1) + change_value_at_fusiontree!(M[end, end], 1, 1) + change_value_at_fusiontree!(R′[1, 1], 1, 1) + change_value_at_fusiontree!(R′[end, end], 1, 1) + M[2:(end - 1), 2:(end - 1)] = R[2:(end - 1), 2:(end - 1)] + R′[1, 2:(end - 1)] = R[1, 2:(end - 1)] # R′[2:end-1,2:end-1] = id(codomain(R[2,2])) # Make diagonal matrix even with codomain and domain being different spaces # Right now I'm assuming that each of the fusiontree is a square block - for (s, f) in fusiontrees(R′[2:end-1,2:end-1]) + for (s, f) in fusiontrees(R′[2:(end - 1), 2:(end - 1)]) if s.uncoupled[1] == f.uncoupled[1] - dim = size(R′[2:end-1,2:end-1][s,f])[1] - R′[2:end-1,2:end-1][s,f] .= diagm(ones(dim)) - end + dim = size(R′[2:(end - 1), 2:(end - 1)][s, f])[1] + R′[2:(end - 1), 2:(end - 1)][s, f] .= diagm(ones(dim)) + end end return M, R′ end -function block_respecting_svd(M::AbstractBlockTensorMap, η::Number; perform_truncation=true) - # Perform block-respecting svd on matrix M, with precision η - # TODO Change it so that it works when sectors are numbered by indices rather than +function block_respecting_svd(M::AbstractBlockTensorMap, truncbelow::Number; perform_truncation = true) + # Perform block-respecting svd on matrix M, with precision truncbelow + # TODO Change it so that it works when sectors are numbered by indices rather than # enumerated - M2 = M[2:end-1,2:end-1] + M2 = M[2:(end - 1), 2:(end - 1)] U, S, V′ = svd_compact(M2) num_sectors = length(blocksectors(S)) new_dims = Vector{Int32}(undef, num_sectors) - # Old code - not working if blocksectors(S) returned Indices rather than Vector - # for n in 1:num_sectors - # csec = blocksectors(S)[n] - # new_dims[n] = blockdim(S.domain, csec) - # if perform_truncation - # for diag_el in block(S, csec).diag - # if diag_el ≤ η - # new_dims[n] -= 1 - # end - # end - # end - # end - - # New code - works for both cases for (n, csec) in enumerate(blocksectors(S)) new_dims[n] = blockdim(S.domain, csec) if perform_truncation for diag_el in block(S, csec).diag - if diag_el ≤ η + if diag_el ≤ truncbelow new_dims[n] -= 1 end end @@ -395,14 +403,6 @@ function block_respecting_svd(M::AbstractBlockTensorMap, η::Number; perform_tru end fdims = Vector{Pair{sectortype(S), Int32}}() - # Old code - # for n in 1:num_sectors - # if new_dims[n] > 0 - # push!(fdims, blocksectors(S)[n] => new_dims[n]) - # end - # end - - # New code - works for both cases for (n, csec) in enumerate(blocksectors(S)) if new_dims[n] > 0 push!(fdims, csec => new_dims[n]) @@ -413,47 +413,38 @@ function block_respecting_svd(M::AbstractBlockTensorMap, η::Number; perform_tru U_new = zeros(ComplexF64, codomain(U) ← red_space) V′_new = zeros(ComplexF64, red_space ← domain(V′)) - #num_sectors2 = length(blocksectors(S_new)) - - # Old code - # for n in 1:num_sectors2 - # csec = blocksectors(S_new)[n] - # cl = length(block(S_new, csec).diag) - # block(S_new, csec).diag[1:end] = block(S, csec).diag[1:cl] - # block(U_new, csec)[:,:] = block(U, csec)[:,1:cl] - # block(V′_new, csec)[:,:] = block(V′, csec)[1:cl,:] - # end - # New code - works for both cases for csec in blocksectors(S_new) cl = length(block(S_new, csec).diag) block(S_new, csec).diag[1:end] = block(S, csec).diag[1:cl] - block(U_new, csec)[:,:] = block(U, csec)[:,1:cl] - block(V′_new, csec)[:,:] = block(V′, csec)[1:cl,:] + block(U_new, csec)[:, :] = block(U, csec)[:, 1:cl] + block(V′_new, csec)[:, :] = block(V′, csec)[1:cl, :] end Û = zeros(ComplexF64, codomain(M) ← domain(M)[1][1] ⊞ domain(U_new)[1] ⊞ domain(M)[1][end]) - Ŝ = zeros(ComplexF64, codomain(M)[1][1] ⊞ codomain(S_new)[1] ⊞ codomain(M)[1][end] ← - domain(M)[1][1] ⊞ domain(S_new)[1] ⊞ domain(M)[1][end]) + Ŝ = zeros( + ComplexF64, codomain(M)[1][1] ⊞ codomain(S_new)[1] ⊞ codomain(M)[1][end] ← + domain(M)[1][1] ⊞ domain(S_new)[1] ⊞ domain(M)[1][end] + ) V̂′ = zeros(ComplexF64, codomain(M)[1][1] ⊞ codomain(V′_new)[1] ⊞ codomain(M)[1][end] ← domain(M)) - change_value_at_fusiontree!(Û[1,1], 1, 1) - change_value_at_fusiontree!(Û[end,end], 1, 1) - change_value_at_fusiontree!(Ŝ[1,1], 1, 1) - change_value_at_fusiontree!(Ŝ[end,end], 1, 1) - change_value_at_fusiontree!(V̂′[1,1], 1, 1) - change_value_at_fusiontree!(V̂′[end,end], 1, 1) + change_value_at_fusiontree!(Û[1, 1], 1, 1) + change_value_at_fusiontree!(Û[end, end], 1, 1) + change_value_at_fusiontree!(Ŝ[1, 1], 1, 1) + change_value_at_fusiontree!(Ŝ[end, end], 1, 1) + change_value_at_fusiontree!(V̂′[1, 1], 1, 1) + change_value_at_fusiontree!(V̂′[end, end], 1, 1) - Û[2:end-1,2:end-1] = U_new - Ŝ[2,2] = S_new - V̂′[2:end-1,2:end-1] = V′_new + Û[2:(end - 1), 2:(end - 1)] = U_new + Ŝ[2, 2] = S_new + V̂′[2:(end - 1), 2:(end - 1)] = V′_new return Û, Ŝ, V̂′ end function reduce_blocks_mpo(W::AbstractBlockTensorMap) - # Function that copies BlockTensorMap W from N x 1 x 1 x M blocked matrix to a - # 3 x 1 x 1 x 3 blocked matrix, with the middle block being joined dimensions of N-2, + # Function that copies BlockTensorMap W from N x 1 x 1 x M blocked matrix to a + # 3 x 1 x 1 x 3 blocked matrix, with the middle block being joined dimensions of N-2, # M-2 blocks from the original matrix # TODO Make exception if either N or M is smaller than 3 @@ -461,9 +452,9 @@ function reduce_blocks_mpo(W::AbstractBlockTensorMap) if size(W)[1] ≤ 3 && size(W)[4] ≤ 3 return W end - - old_cod = codomain(W)[1][2:end-1] - old_dom = domain(W)[2][2:end-1] + + old_cod = codomain(W)[1][2:(end - 1)] + old_dom = domain(W)[2][2:(end - 1)] new_cod = spacetype(W)([sec => blockdim(old_cod, sec) for sec in sectors(old_cod)]...) new_dom = spacetype(W)([sec => blockdim(old_dom, sec) for sec in sectors(old_dom)]...) @@ -482,51 +473,51 @@ function reduce_blocks_mpo(W::AbstractBlockTensorMap) W2 = zeros(scalartype(W), Vcod ← Vdom) # Corners can be copied trivially - W2[1,1,1,1] = W[1,1,1,1] - W2[end,1,1,end] = W[end,1,1,end] - W2[1,1,1,end] = W[1,1,1,end] - W2[end,1,1,1] = W[end,1,1,1] + W2[1, 1, 1, 1] = W[1, 1, 1, 1] + W2[end, 1, 1, end] = W[end, 1, 1, end] + W2[1, 1, 1, end] = W[1, 1, 1, end] + W2[end, 1, 1, 1] = W[end, 1, 1, 1] # Treat the cases with smaller codomain or domain separately if size(W)[1] ≤ 2 # Middle upperblock - for (s,f) in fusiontrees(W[1,1,1,2:end-1]) - W2[1,1,1,2][s,f] = W[1,1,1,2:end-1][s,f] + for (s, f) in fusiontrees(W[1, 1, 1, 2:(end - 1)]) + W2[1, 1, 1, 2][s, f] = W[1, 1, 1, 2:(end - 1)][s, f] end - # Middle bottomblock (if size(W)[1] is 1, this is the same block as in + # Middle bottomblock (if size(W)[1] is 1, this is the same block as in # previous loop) - for (s,f) in fusiontrees(W[end,1,1,2:end-1]) - W2[end,1,1,2][s,f] = W[end,1,1,2:end-1][s,f] + for (s, f) in fusiontrees(W[end, 1, 1, 2:(end - 1)]) + W2[end, 1, 1, 2][s, f] = W[end, 1, 1, 2:(end - 1)][s, f] end elseif size(W)[4] ≤ 2 # Middle leftblock - for (s,f) in fusiontrees(W[2:end-1,1,1,1]) - W2[2,1,1,1][s,f] = W[2:end-1,1,1,1][s,f] + for (s, f) in fusiontrees(W[2:(end - 1), 1, 1, 1]) + W2[2, 1, 1, 1][s, f] = W[2:(end - 1), 1, 1, 1][s, f] end # Middle rightblock (if size(W)[4] is 1, this is the same block as in previous loop) - for (s,f) in fusiontrees(W[2:end-1,1,1,end]) - W2[2,1,1,end][s,f] = W[2:end-1,1,1,end][s,f] + for (s, f) in fusiontrees(W[2:(end - 1), 1, 1, end]) + W2[2, 1, 1, end][s, f] = W[2:(end - 1), 1, 1, end][s, f] end else # All other blocks can be dealt with fusiontrees # Middle block (most important) - for (s,f) in fusiontrees(W[2:end-1,1,1,2:end-1]) - W2[2,1,1,2][s,f] = W[2:end-1,1,1,2:end-1][s,f] + for (s, f) in fusiontrees(W[2:(end - 1), 1, 1, 2:(end - 1)]) + W2[2, 1, 1, 2][s, f] = W[2:(end - 1), 1, 1, 2:(end - 1)][s, f] end # Upper block - for (s,f) in fusiontrees(W[1,1,1,2:end-1]) - W2[1,1,1,2][s,f] = W[1,1,1,2:end-1][s,f] + for (s, f) in fusiontrees(W[1, 1, 1, 2:(end - 1)]) + W2[1, 1, 1, 2][s, f] = W[1, 1, 1, 2:(end - 1)][s, f] end # Bottom block - for (s,f) in fusiontrees(W[end,1,1,2:end-1]) - W2[3,1,1,2][s,f] = W[end,1,1,2:end-1][s,f] + for (s, f) in fusiontrees(W[end, 1, 1, 2:(end - 1)]) + W2[3, 1, 1, 2][s, f] = W[end, 1, 1, 2:(end - 1)][s, f] end # Left block - for (s,f) in fusiontrees(W[2:end-1,1,1,1]) - W2[2,1,1,1][s,f] = W[2:end-1,1,1,1][s,f] + for (s, f) in fusiontrees(W[2:(end - 1), 1, 1, 1]) + W2[2, 1, 1, 1][s, f] = W[2:(end - 1), 1, 1, 1][s, f] end # Right block - for (s,f) in fusiontrees(W[2:end-1,1,1,end]) - W2[2,1,1,3][s,f] = W[2:end-1,1,1,end][s,f] + for (s, f) in fusiontrees(W[2:(end - 1), 1, 1, end]) + W2[2, 1, 1, 3][s, f] = W[2:(end - 1), 1, 1, end][s, f] end end return W2 @@ -536,10 +527,10 @@ function check_how_far_from_id(R::AbstractBlockTensorMap) # Function checking how far the diagonal elements of R are from either # +1 or -1. Assume R is 3x3 and [1,1] and [3,3] blocks are always id tdev = 0 - for (s,f) in fusiontrees(R[2,2]) + for (s, f) in fusiontrees(R[2, 2]) if s.uncoupled[1] == f.uncoupled[1] - for d in size(R[2,2][s,f])[1] - tdev += min(abs(R[2,2][s,f][d,d] - 1), abs(R[2,2][s,f][d,d] + 1)) + for d in size(R[2, 2][s, f])[1] + tdev += min(abs(R[2, 2][s, f][d, d] - 1), abs(R[2, 2][s, f][d, d] + 1)) end end end @@ -547,22 +538,19 @@ function check_how_far_from_id(R::AbstractBlockTensorMap) end function create_impoham_from_mpos_msites(Qs::Vector{BlockTensorMap}) - # Given the Vector{BlockTensorMap} Qs, return the multisite InfiniteMPOHamiltonian + # Given the Vector{BlockTensorMap} Qs, return the multisite InfiniteMPOHamiltonian N = length(Qs) sctype = scalartype(Qs[1]) sptype = spacetype(Qs[1]) - # TA = AbstractTensorMap{sctype, sptype, 2, 2} - # TB = AbstractTensorMap{sctype, sptype, 2, 1} - # TC = AbstractTensorMap{sctype, sptype, 1, 2} - # TD = AbstractTensorMap{sctype, sptype, 1, 1} + CType = JordanMPOTensor{sctype, sptype, Vector{sctype}} Cs = Vector{CType}(undef, N) for n in 1:N W = SparseBlockTensorMap(Qs[n]) - A = W[2:(end-1), 1, 1, 2:(end-1)] - B = removeunit(W[2:(end-1), 1, 1, end], 4) - C = removeunit(W[1,1,1,2:(end-1)], 1) - D = removeunit(removeunit(W[1,1,1,end:end], 4), 1) + A = W[2:(end - 1), 1, 1, 2:(end - 1)] + B = removeunit(W[2:(end - 1), 1, 1, end], 4) + C = removeunit(W[1, 1, 1, 2:(end - 1)], 1) + D = removeunit(removeunit(W[1, 1, 1, end:end], 4), 1) W = JordanMPOTensor(space(W), A, B, C, D) Cs[n] = W end @@ -571,9 +559,10 @@ function create_impoham_from_mpos_msites(Qs::Vector{BlockTensorMap}) end function left_canonical_mpo_infinite_iter_msites( - H::Union{InfiniteMPOHamiltonian, V}, - η=10^-10 - ) where {T<:SparseBlockTensorMap{<:Any, <:Any, <:Any, 2, 2, 4}, V<:Vector{T}} + H::Union{InfiniteMPOHamiltonian, V}; + precision = 10^-10, + verbosity::Int = 0 + ) where {T <: SparseBlockTensorMap{<:Any, <:Any, <:Any, 2, 2, 4}, V <: Vector{T}} # Return the left_canonical form of the iMPO, H can have multiple sites in the unit cell N = length(H) # Number of sites @@ -581,8 +570,10 @@ function left_canonical_mpo_infinite_iter_msites( Ls = fill(id(ComplexF64, codomain(H[1])[1]), N) Rs = fill(id(ComplexF64, codomain(H[1])[1]), N) Q = similar(H[1]) - printstyled("┌─── Started iterative left orthogonalization ────\n", color=:cyan) - while sum(εs)/N > η + if verbosity > 1 + printstyled("┌─── Started iterative left orthogonalization ────\n", color = :cyan) + end + while sum(εs) / N > precision # First store the R matrices and change them H into Qs for n in 1:N Q, R = qr_block_respecting(H[n]) @@ -593,28 +584,35 @@ function left_canonical_mpo_infinite_iter_msites( end # Then assign R[(n+1)//N] * Q[n] to H[n] for n in 1:N - @tensor H[n][a,i;j,b] = Rs[mod1(n-1,N)][a;c] * H[n][c,i;j,b] + @tensor H[n][a, i; j, b] = Rs[mod1(n - 1, N)][a;c] * H[n][c, i; j, b] + end + if verbosity > 1 + printstyled("| total err = $(sum(εs) / N)\n", color = :white) end - printstyled("| total err = $(sum(εs)/N)\n", color=:white) end - printstyled("└─── Left orthogonalization finished correctly ──\n", color=:cyan) + if verbosity > 1 + printstyled("└─── Left orthogonalization finished correctly ──\n", color = :cyan) + end return H, Ls end function right_canonical_mpo_infinite_iter_msites( - H::Union{InfiniteMPOHamiltonian, V}, - η=10^-10 - ) where {T<:SparseBlockTensorMap{<:Any, <:Any, <:Any, 2, 2, 4}, V<:Vector{T}} - # Return the right_canonical form of the iMPO, H can have multiple sites in the unit + H::Union{InfiniteMPOHamiltonian, V}; + precision = 10^-10, + verbosity::Int = 0 + ) where {T <: SparseBlockTensorMap{<:Any, <:Any, <:Any, 2, 2, 4}, V <: Vector{T}} + # Return the right_canonical form of the iMPO, H can have multiple sites in the unit # cell N = length(H) # Number of sites - + εs = fill(Inf, N) Rs = fill(id(ComplexF64, codomain(H[1])[1]), N) Ls = fill(id(ComplexF64, codomain(H[1])[1]), N) Q = similar(H[N]) - printstyled("┌─── Started iterative right orthogonalization ────\n", color=:cyan) - while sum(εs)/N > η + if verbosity > 1 + printstyled("┌─── Started iterative right orthogonalization ────\n", color = :cyan) + end + while sum(εs) / N > precision # First store the R matrices and change them H into Qs for n in N:-1:1 L, Q = lq_block_respecting(H[n]) @@ -625,10 +623,14 @@ function right_canonical_mpo_infinite_iter_msites( end # Then assign R[(n+1)//N] * Q[n] to H[n] for n in N:-1:1 - @tensor H[n][a,i;j,b] = H[n][a,i;j,c] * Ls[mod1(n+1, N)][c,b] + @tensor H[n][a, i; j, b] = H[n][a, i; j, c] * Ls[mod1(n + 1, N)][c, b] end - printstyled("| total err = $(sum(εs)/N)\n", color=:white) + if verbosity > 1 + printstyled("| total err = $(sum(εs) / N)\n", color = :white) + end + end + if verbosity > 1 + printstyled("└─── Right orthogonalization finished correctly ──\n", color = :cyan) end - printstyled("└─── Right orthogonalization finished correctly ──\n", color=:cyan) return Rs, H -end \ No newline at end of file +end From 296354260bea54c036b3c7f1628383e08f548c23 Mon Sep 17 00:00:00 2001 From: herviou Date: Wed, 15 Jul 2026 14:26:39 +0200 Subject: [PATCH 3/3] Actually fixing formatting. --- test/algorithms/dynamical_dmrg.jl | 1 - test/operators/finite_mpo_compression.jl | 152 +++++++-------- test/operators/infinite_mpo_compression.jl | 206 +++++++++++---------- 3 files changed, 190 insertions(+), 169 deletions(-) diff --git a/test/algorithms/dynamical_dmrg.jl b/test/algorithms/dynamical_dmrg.jl index 146c43af3..8e07757f4 100644 --- a/test/algorithms/dynamical_dmrg.jl +++ b/test/algorithms/dynamical_dmrg.jl @@ -60,4 +60,3 @@ end @test data ≈ predicted atol = 1.0e-8 end end - diff --git a/test/operators/finite_mpo_compression.jl b/test/operators/finite_mpo_compression.jl index 07c259cf6..8b4cf3dec 100644 --- a/test/operators/finite_mpo_compression.jl +++ b/test/operators/finite_mpo_compression.jl @@ -17,24 +17,28 @@ end # Utility functions to create Hamiltonians for testing function create_long_range_ising_symmetries_finite(M, k) - chain = fill(Z2Space(0=>1, 1=>1), M) + chain = fill(Z2Space(0 => 1, 1 => 1), M) ZZ = S_zz(Z2Irrep) X = S_x(Z2Irrep) single_site_operators = [i => X for i in 1:M] - two_site_operators = [(i,i+j) => ZZ for i in 1:M for j in 1:k if i+j ≤ M] - H = FiniteMPOHamiltonian(chain, single_site_operators..., - two_site_operators...) + two_site_operators = [(i, i + j) => ZZ for i in 1:M for j in 1:k if i + j ≤ M] + H = FiniteMPOHamiltonian( + chain, single_site_operators..., + two_site_operators... + ) return H end function create_long_range_ising_symmetries_random_finite(M, k, σⱼ, σₕ) - chain = fill(Z2Space(0=>1, 1=>1), M) + chain = fill(Z2Space(0 => 1, 1 => 1), M) ZZ = S_zz(Z2Irrep) X = S_x(Z2Irrep) single_site_operators = [i => (σₕ * randn()) * X for i in 1:M] - two_site_operators = [(i,i+j) => (1 + σⱼ * randn()) * ZZ for i in 1:M for j in 1:k if i+j ≤ M] - H = FiniteMPOHamiltonian(chain, single_site_operators..., - two_site_operators...) + two_site_operators = [(i, i + j) => (1 + σⱼ * randn()) * ZZ for i in 1:M for j in 1:k if i + j ≤ M] + H = FiniteMPOHamiltonian( + chain, single_site_operators..., + two_site_operators... + ) return H end @@ -108,15 +112,15 @@ end range = 4 # cutoff for the maximum range between interactions trunc_st = notrunc() - mps = FiniteMPS(L, Z2Space(0=>1, 1=>1), Z2Space(0=>D, 1=>D)) + mps = FiniteMPS(L, Z2Space(0 => 1, 1 => 1), Z2Space(0 => D, 1 => D)) H = create_long_range_ising_symmetries_finite(L, range) - find_groundstate!(mps, H, DMRG2(; maxiter=100, trscheme=trunc_st)) + find_groundstate!(mps, H, DMRG2(; maxiter = 100, trscheme = trunc_st)) E0 = expectation_value(mps, H) #println(" = $real(E0)") - mps2 = FiniteMPS(L, Z2Space(0=>1, 1=>1), Z2Space(0=>D, 1=>D)) + mps2 = FiniteMPS(L, Z2Space(0 => 1, 1 => 1), Z2Space(0 => D, 1 => D)) H2, Rs = mpo_compression(H) - find_groundstate!(mps2, H2, DMRG2(; maxiter=100, trscheme=trunc_st)) + find_groundstate!(mps2, H2, DMRG2(; maxiter = 100, trscheme = trunc_st)) E1 = expectation_value(mps2, H2) end @@ -126,50 +130,50 @@ end trunc_st = trunctol() mps1 = FiniteMPS(M, ℂ^2, ℂ^D) - mps2 = FiniteMPS(M, Z2Space(0=>1, 1=>1), Z2Space(0=>D, 1=>D)) - mps3 = FiniteMPS(M, Vect[FermionParity](0=>1, 1=>1), Vect[FermionParity](0=>D, 1=>D)) + mps2 = FiniteMPS(M, Z2Space(0 => 1, 1 => 1), Z2Space(0 => D, 1 => D)) + mps3 = FiniteMPS(M, Vect[FermionParity](0 => 1, 1 => 1), Vect[FermionParity](0 => D, 1 => D)) H1 = transverse_field_ising_trivial_finite(M) H2 = transverse_field_ising_z2_finite(M) H3 = transverse_field_ising_fermion_parity_finite(M) # Find groundstate of non-compressed Hamiltonian - find_groundstate!(mps1, H1, DMRG2(; maxiter=100, trscheme=trunc_st)) + find_groundstate!(mps1, H1, DMRG2(; maxiter = 100, trscheme = trunc_st)) E0_1 = expectation_value(mps1, H1) - find_groundstate!(mps2, H2, DMRG2(; maxiter=100, trscheme=trunc_st)) + find_groundstate!(mps2, H2, DMRG2(; maxiter = 100, trscheme = trunc_st)) E0_2 = expectation_value(mps2, H2) - find_groundstate!(mps3, H3, DMRG2(; maxiter=100, trscheme=trunc_st)) + find_groundstate!(mps3, H3, DMRG2(; maxiter = 100, trscheme = trunc_st)) E0_3 = expectation_value(mps3, H3) # Compress MPO mps1c = FiniteMPS(M, ℂ^2, ℂ^D) - mps2c = FiniteMPS(M, Z2Space(0=>1, 1=>1), Z2Space(0=>D, 1=>D)) - mps3c = FiniteMPS(M, Vect[FermionParity](0=>1, 1=>1), Vect[FermionParity](0=>D, 1=>D)) + mps2c = FiniteMPS(M, Z2Space(0 => 1, 1 => 1), Z2Space(0 => D, 1 => D)) + mps3c = FiniteMPS(M, Vect[FermionParity](0 => 1, 1 => 1), Vect[FermionParity](0 => D, 1 => D)) H1c, _ = mpo_compression(H1, 10^-10) H2c, _ = mpo_compression(H2, 10^-10) H3c, _ = mpo_compression(H3, 10^-10) # Find groundstate of compressed Hamiltonian - find_groundstate!(mps1c, H1c, DMRG2(; maxiter=100, trscheme=trunc_st)) + find_groundstate!(mps1c, H1c, DMRG2(; maxiter = 100, trscheme = trunc_st)) E0_1c = expectation_value(mps1c, H1c) - find_groundstate!(mps2c, H2c, DMRG2(; maxiter=100, trscheme=trunc_st)) + find_groundstate!(mps2c, H2c, DMRG2(; maxiter = 100, trscheme = trunc_st)) E0_2c = expectation_value(mps2c, H2c) - find_groundstate!(mps3c, H3c, DMRG2(; maxiter=100, trscheme=trunc_st)) + find_groundstate!(mps3c, H3c, DMRG2(; maxiter = 100, trscheme = trunc_st)) E0_3c = expectation_value(mps3c, H3c) - # Assert that the ground state energies are equal up to a precision, display the + # Assert that the ground state energies are equal up to a precision, display the # compression @assert E0_1 ≈ E0_1c @assert E0_2 ≈ E0_2c @assert E0_3 ≈ E0_3c println("Max dim uncompressed H1: $(max_dim_ham(H1)) || Max dim compressed H1c: - $(max_dim_ham(H1c))" ) + $(max_dim_ham(H1c))") println("Max dim uncompressed H2: $(max_dim_ham(H2)) || Max dim compressed H2c: - $(max_dim_ham(H2c))" ) + $(max_dim_ham(H2c))") println("Max dim uncompressed H3: $(max_dim_ham(H3)) || Max dim compressed H3c: - $(max_dim_ham(H3c))" ) + $(max_dim_ham(H3c))") end @@ -178,29 +182,29 @@ end D = 10 # Max bond dimension trunc_st = trunctol() - mps1 = FiniteMPS(M, Vect[FermionParity](0=>1, 1=>1), Vect[FermionParity](0=>D, 1=>D)) + mps1 = FiniteMPS(M, Vect[FermionParity](0 => 1, 1 => 1), Vect[FermionParity](0 => D, 1 => D)) H1 = kitaev_model_finite(M) # Find groundstate of non-compressed Hamiltonian - find_groundstate!(mps1, H1, DMRG2(; maxiter=100, trscheme=trunc_st)) + find_groundstate!(mps1, H1, DMRG2(; maxiter = 100, trscheme = trunc_st)) E0_1 = expectation_value(mps1, H1) # Compress MPO - mps1c = FiniteMPS(M, Vect[FermionParity](0=>1, 1=>1), Vect[FermionParity](0=>D, 1=>D)) + mps1c = FiniteMPS(M, Vect[FermionParity](0 => 1, 1 => 1), Vect[FermionParity](0 => D, 1 => D)) H1c, _ = mpo_compression(H1, 10^-10) # Find groundstate of compressed Hamiltonian - find_groundstate!(mps1c, H1c, DMRG2(; maxiter=100, trscheme=trunc_st)) + find_groundstate!(mps1c, H1c, DMRG2(; maxiter = 100, trscheme = trunc_st)) E0_1c = expectation_value(mps1c, H1c) - # Assert that the ground state energies are equal up to a precision, display the + # Assert that the ground state energies are equal up to a precision, display the # compression @assert E0_1 ≈ E0_1c println("Max dim uncompressed H1: $(max_dim_ham(H1)) || Max dim compressed H1c: - $(max_dim_ham(H1c))" ) + $(max_dim_ham(H1c))") end @@ -210,50 +214,50 @@ end trunc_st = trunctol() mps1 = FiniteMPS(M, ℂ^3, ℂ^D) - mps2 = FiniteMPS(M, U1Space(0=>1, 1=>1, -1=>1), U1Space(0=>D, 1=>D, -1=>D)) - mps3 = FiniteMPS(M, SU2Space(1 => 1), SU2Space(1=>D)) + mps2 = FiniteMPS(M, U1Space(0 => 1, 1 => 1, -1 => 1), U1Space(0 => D, 1 => D, -1 => D)) + mps3 = FiniteMPS(M, SU2Space(1 => 1), SU2Space(1 => D)) H1 = heisenberg_XXX_trivial_finite(M) H2 = heisenberg_XXX_U1_finite(M) H3 = heisenberg_XXX_SU2_finite(M) # Find groundstate of non-compressed Hamiltonian - find_groundstate!(mps1, H1, DMRG2(; maxiter=100, trscheme=trunc_st)) + find_groundstate!(mps1, H1, DMRG2(; maxiter = 100, trscheme = trunc_st)) E0_1 = expectation_value(mps1, H1) - find_groundstate!(mps2, H2, DMRG2(; maxiter=100, trscheme=trunc_st)) + find_groundstate!(mps2, H2, DMRG2(; maxiter = 100, trscheme = trunc_st)) E0_2 = expectation_value(mps2, H2) - find_groundstate!(mps3, H3, DMRG2(; maxiter=100, trscheme=trunc_st)) + find_groundstate!(mps3, H3, DMRG2(; maxiter = 100, trscheme = trunc_st)) E0_3 = expectation_value(mps3, H3) # Compress MPO mps1c = FiniteMPS(M, ℂ^3, ℂ^D) - mps2c = FiniteMPS(M, U1Space(0=>1, 1=>1, -1=>1), U1Space(0=>D, 1=>D, -1=>D)) - mps3c = FiniteMPS(M, SU2Space(1 => 1), SU2Space(1=>D)) + mps2c = FiniteMPS(M, U1Space(0 => 1, 1 => 1, -1 => 1), U1Space(0 => D, 1 => D, -1 => D)) + mps3c = FiniteMPS(M, SU2Space(1 => 1), SU2Space(1 => D)) H1c, _ = mpo_compression(H1, 10^-10) H2c, _ = mpo_compression(H2, 10^-10) H3c, _ = mpo_compression(H3, 10^-10) # Find groundstate of compressed Hamiltonian - find_groundstate!(mps1c, H1c, DMRG2(; maxiter=100, trscheme=trunc_st)) + find_groundstate!(mps1c, H1c, DMRG2(; maxiter = 100, trscheme = trunc_st)) E0_1c = expectation_value(mps1c, H1c) - find_groundstate!(mps2c, H2c, DMRG2(; maxiter=100, trscheme=trunc_st)) + find_groundstate!(mps2c, H2c, DMRG2(; maxiter = 100, trscheme = trunc_st)) E0_2c = expectation_value(mps2c, H2c) - find_groundstate!(mps3c, H3c, DMRG2(; maxiter=100, trscheme=trunc_st)) + find_groundstate!(mps3c, H3c, DMRG2(; maxiter = 100, trscheme = trunc_st)) E0_3c = expectation_value(mps3c, H3c) - # Assert that the ground state energies are equal up to a precision, display the + # Assert that the ground state energies are equal up to a precision, display the # compression @assert E0_1 ≈ E0_1c @assert E0_2 ≈ E0_2c @assert E0_3 ≈ E0_3c println("Max dim uncompressed H1: $(max_dim_ham(H1)) || Max dim compressed H1c: - $(max_dim_ham(H1c))" ) + $(max_dim_ham(H1c))") println("Max dim uncompressed H2: $(max_dim_ham(H2)) || Max dim compressed H2c: - $(max_dim_ham(H2c))" ) + $(max_dim_ham(H2c))") println("Max dim uncompressed H3: $(max_dim_ham(H3)) || Max dim compressed H3c: - $(max_dim_ham(H3c))" ) + $(max_dim_ham(H3c))") end @@ -263,50 +267,50 @@ end trunc_st = trunctol() mps1 = FiniteMPS(M, ℂ^3, ℂ^D) - mps2 = FiniteMPS(M, U1Space(0=>1, 1=>1, -1=>1), U1Space(0=>D, 1=>D, -1=>D)) - mps3 = FiniteMPS(M, SU2Space(1 => 1), SU2Space(1=>D)) + mps2 = FiniteMPS(M, U1Space(0 => 1, 1 => 1, -1 => 1), U1Space(0 => D, 1 => D, -1 => D)) + mps3 = FiniteMPS(M, SU2Space(1 => 1), SU2Space(1 => D)) H1 = bilinear_biquadratic_model_trivial_finite(M) H2 = bilinear_biquadratic_model_U1_finite(M) H3 = bilinear_biquadratic_model_SU2_finite(M) # Find groundstate of non-compressed Hamiltonian - find_groundstate!(mps1, H1, DMRG2(; maxiter=100, trscheme=trunc_st)) + find_groundstate!(mps1, H1, DMRG2(; maxiter = 100, trscheme = trunc_st)) E0_1 = expectation_value(mps1, H1) - find_groundstate!(mps2, H2, DMRG2(; maxiter=100, trscheme=trunc_st)) + find_groundstate!(mps2, H2, DMRG2(; maxiter = 100, trscheme = trunc_st)) E0_2 = expectation_value(mps2, H2) - find_groundstate!(mps3, H3, DMRG2(; maxiter=100, trscheme=trunc_st)) + find_groundstate!(mps3, H3, DMRG2(; maxiter = 100, trscheme = trunc_st)) E0_3 = expectation_value(mps3, H3) # Compress MPO mps1c = FiniteMPS(M, ℂ^3, ℂ^D) - mps2c = FiniteMPS(M, U1Space(0=>1, 1=>1, -1=>1), U1Space(0=>D, 1=>D, -1=>D)) - mps3c = FiniteMPS(M, SU2Space(1 => 1), SU2Space(1=>D)) + mps2c = FiniteMPS(M, U1Space(0 => 1, 1 => 1, -1 => 1), U1Space(0 => D, 1 => D, -1 => D)) + mps3c = FiniteMPS(M, SU2Space(1 => 1), SU2Space(1 => D)) H1c, _ = mpo_compression(H1, 10^-10) H2c, _ = mpo_compression(H2, 10^-10) H3c, _ = mpo_compression(H3, 10^-10) # Find groundstate of compressed Hamiltonian - find_groundstate!(mps1c, H1c, DMRG2(; maxiter=100, trscheme=trunc_st)) + find_groundstate!(mps1c, H1c, DMRG2(; maxiter = 100, trscheme = trunc_st)) E0_1c = expectation_value(mps1c, H1c) - find_groundstate!(mps2c, H2c, DMRG2(; maxiter=100, trscheme=trunc_st)) + find_groundstate!(mps2c, H2c, DMRG2(; maxiter = 100, trscheme = trunc_st)) E0_2c = expectation_value(mps2c, H2c) - find_groundstate!(mps3c, H3c, DMRG2(; maxiter=100, trscheme=trunc_st)) + find_groundstate!(mps3c, H3c, DMRG2(; maxiter = 100, trscheme = trunc_st)) E0_3c = expectation_value(mps3c, H3c) - # Assert that the ground state energies are equal up to a precision, display the + # Assert that the ground state energies are equal up to a precision, display the # compression @assert E0_1 ≈ E0_1c @assert E0_2 ≈ E0_2c @assert E0_3 ≈ E0_3c println("Max dim uncompressed H1: $(max_dim_ham(H1)) || Max dim compressed H1c: - $(max_dim_ham(H1c))" ) + $(max_dim_ham(H1c))") println("Max dim uncompressed H2: $(max_dim_ham(H2)) || Max dim compressed H2c: - $(max_dim_ham(H2c))" ) + $(max_dim_ham(H2c))") println("Max dim uncompressed H3: $(max_dim_ham(H3)) || Max dim compressed H3c: - $(max_dim_ham(H3c))" ) + $(max_dim_ham(H3c))") end @@ -316,39 +320,39 @@ end trunc_st = trunctol() mps1 = FiniteMPS(M, ℂ^3, ℂ^D) - mps2 = FiniteMPS(M, Z3Space(0=>1, 1=>1, 2=>1), Z3Space(0=>1, 1=>1, 2=>1)) + mps2 = FiniteMPS(M, Z3Space(0 => 1, 1 => 1, 2 => 1), Z3Space(0 => 1, 1 => 1, 2 => 1)) H1 = quantum_potts_trivial_finite(M) H2 = quantum_potts_Z3_finite(M) # Find groundstate of non-compressed Hamiltonian - find_groundstate!(mps1, H1, DMRG2(; maxiter=100, trscheme=trunc_st)) + find_groundstate!(mps1, H1, DMRG2(; maxiter = 100, trscheme = trunc_st)) E0_1 = expectation_value(mps1, H1) - find_groundstate!(mps2, H2, DMRG2(; maxiter=100, trscheme=trunc_st)) + find_groundstate!(mps2, H2, DMRG2(; maxiter = 100, trscheme = trunc_st)) E0_2 = expectation_value(mps2, H2) # Compress MPO mps1c = FiniteMPS(M, ℂ^3, ℂ^D) - mps2c = FiniteMPS(M, Z3Space(0=>1, 1=>1, 2=>1), Z3Space(0=>1, 1=>1, 2=>1)) + mps2c = FiniteMPS(M, Z3Space(0 => 1, 1 => 1, 2 => 1), Z3Space(0 => 1, 1 => 1, 2 => 1)) H1c, _ = mpo_compression(H1, 10^-10) H2c, _ = mpo_compression(H2, 10^-10) # Find groundstate of compressed Hamiltonian - find_groundstate!(mps1c, H1c, DMRG2(; maxiter=100, trscheme=trunc_st)) + find_groundstate!(mps1c, H1c, DMRG2(; maxiter = 100, trscheme = trunc_st)) E0_1c = expectation_value(mps1c, H1c) - find_groundstate!(mps2c, H2c, DMRG2(; maxiter=100, trscheme=trunc_st)) + find_groundstate!(mps2c, H2c, DMRG2(; maxiter = 100, trscheme = trunc_st)) E0_2c = expectation_value(mps2c, H2c) - # Assert that the ground state energies are equal up to a precision, display the + # Assert that the ground state energies are equal up to a precision, display the # compression @assert E0_1 ≈ E0_1c @assert E0_2 ≈ E0_2c println("Max dim uncompressed H1: $(max_dim_ham(H1)) || Max dim compressed H1c: - $(max_dim_ham(H1c))" ) + $(max_dim_ham(H1c))") println("Max dim uncompressed H2: $(max_dim_ham(H2)) || Max dim compressed H2c: - $(max_dim_ham(H2c))" ) + $(max_dim_ham(H2c))") end @@ -362,14 +366,14 @@ end range = 4 # cutoff for the maximum range between interactions trunc_st = notrunc() - mps = FiniteMPS(L, Z2Space(0=>1, 1=>1), Z2Space(0=>D, 1=>D)) + mps = FiniteMPS(L, Z2Space(0 => 1, 1 => 1), Z2Space(0 => D, 1 => D)) H = create_long_range_ising_symmetries_random_finite(L, range, 0.3, 0.4) - find_groundstate!(mps, H, DMRG2(; maxiter=100, trscheme=trunc_st)) + find_groundstate!(mps, H, DMRG2(; maxiter = 100, trscheme = trunc_st)) E0 = expectation_value(mps, H) #println(" = $real(E0)") - mps2 = FiniteMPS(L, Z2Space(0=>1, 1=>1), Z2Space(0=>D, 1=>D)) + mps2 = FiniteMPS(L, Z2Space(0 => 1, 1 => 1), Z2Space(0 => D, 1 => D)) H2, Rs = mpo_compression(H) - find_groundstate!(mps2, H2, DMRG2(; maxiter=100, trscheme=trunc_st)) + find_groundstate!(mps2, H2, DMRG2(; maxiter = 100, trscheme = trunc_st)) E1 = expectation_value(mps2, H2) -end \ No newline at end of file +end diff --git a/test/operators/infinite_mpo_compression.jl b/test/operators/infinite_mpo_compression.jl index caebc96c6..e9b7ce082 100644 --- a/test/operators/infinite_mpo_compression.jl +++ b/test/operators/infinite_mpo_compression.jl @@ -17,37 +17,41 @@ end # Utility functions to create Hamiltonians for testing function create_long_range_ising_symmetries_infinite(k) - infinite_chain = PeriodicVector([Z2Space(0=>1, 1=>1)]) + infinite_chain = PeriodicVector([Z2Space(0 => 1, 1 => 1)]) ZZ = S_zz(Z2Irrep) X = S_x(Z2Irrep) - two_site_operators = [(1,j) => -ZZ for j in 1:k] + two_site_operators = [(1, j) => -ZZ for j in 1:k] H = InfiniteMPOHamiltonian(infinite_chain, two_site_operators...) return H end function create_long_range_ising_symmetries_infinite_msite_random(M, k, σⱼ, σₕ) - V = Z2Space(0=>1, 1=>1) + V = Z2Space(0 => 1, 1 => 1) infinite_chain = PeriodicVector(fill(V, M)) ZZ = S_zz(Z2Irrep) X = S_x(Z2Irrep) single_site_operators = [i => (σₕ * randn()) * X for i in 1:M] - two_site_operators = [(i,i+j) => (1 + σⱼ * randn()) * ZZ for i in 1:M for j in 1:k] - H = InfiniteMPOHamiltonian(infinite_chain, single_site_operators..., - two_site_operators...) + two_site_operators = [(i, i + j) => (1 + σⱼ * randn()) * ZZ for i in 1:M for j in 1:k] + H = InfiniteMPOHamiltonian( + infinite_chain, single_site_operators..., + two_site_operators... + ) return H end function create_long_range_ising_symmetries_infinite_twosite_diff_int_strength(k, J, h) - infinite_chain = PeriodicVector([Z2Space(0=>1, 1=>1), Z2Space(0=>1, 1=>1)]) + infinite_chain = PeriodicVector([Z2Space(0 => 1, 1 => 1), Z2Space(0 => 1, 1 => 1)]) ZZ = S_zz(Z2Irrep) X = S_x(Z2Irrep) - two_site_operators_J12 = [(a,a+j) => J[1] * ZZ for a in [1 2] for j in 1:2:2k] - two_site_operators_J11 = [(1,j) => J[2] * ZZ for j in 3:2:2k] - two_site_operators_J22 = [(2,j) => J[3] * ZZ for j in 4:2:2k] + two_site_operators_J12 = [(a, a + j) => J[1] * ZZ for a in [1 2] for j in 1:2:2k] + two_site_operators_J11 = [(1, j) => J[2] * ZZ for j in 3:2:2k] + two_site_operators_J22 = [(2, j) => J[3] * ZZ for j in 4:2:2k] single_site_operators = [1 => h[1] * X, 2 => h[2] * X] - H = InfiniteMPOHamiltonian(infinite_chain, single_site_operators..., - two_site_operators_J11..., two_site_operators_J12..., - two_site_operators_J22...) + H = InfiniteMPOHamiltonian( + infinite_chain, single_site_operators..., + two_site_operators_J11..., two_site_operators_J12..., + two_site_operators_J22... + ) return H end @@ -108,9 +112,9 @@ end range = 2 # cutoff for the maximum range between interactions trunc_st = notrunc() - mps = InfiniteMPS(Z2Space(0=>1, 1=>1), Z2Space(0=>D, 1=>D)) + mps = InfiniteMPS(Z2Space(0 => 1, 1 => 1), Z2Space(0 => D, 1 => D)) H = create_long_range_ising_symmetries_infinite(range) - mps, = find_groundstate(mps, H, VUMPS(;maxiter=10)) + mps, = find_groundstate(mps, H, VUMPS(; maxiter = 10)) E0 = expectation_value(mps, H) #println(" = $real(E0)") @@ -118,10 +122,10 @@ end # H2, Rs = mpo_finite_compression(H) # find_groundstate!(mps2, H2, DMRG2(; maxiter=100, trscheme=trunc_st)) # E1 = expectation_value(mps2, H2) - mps2 = InfiniteMPS(Z2Space(0=>1, 1=>1), Z2Space(0=>D, 1=>D)) + mps2 = InfiniteMPS(Z2Space(0 => 1, 1 => 1), Z2Space(0 => D, 1 => D)) Q, P = mpo_compression(H) H2 = Q - mps2, = find_groundstate(mps2, H2, VUMPS(;maxiter=100)) + mps2, = find_groundstate(mps2, H2, VUMPS(; maxiter = 100)) E0 = expectation_value(mps2, H2) end @@ -134,12 +138,12 @@ end J = [1, 1, 1] h = [0.001, 0.001] - spacephys = fill(Z2Space(0=>1, 1=>1), 2) - space_virt = fill(Z2Space(0=>D, 1=>D), 2) + spacephys = fill(Z2Space(0 => 1, 1 => 1), 2) + space_virt = fill(Z2Space(0 => D, 1 => D), 2) mps = InfiniteMPS(spacephys, space_virt) H = create_long_range_ising_symmetries_infinite_twosite_diff_int_strength(range, J, h) - mps, = find_groundstate(mps, H, VUMPS(;maxiter=100)) + mps, = find_groundstate(mps, H, VUMPS(; maxiter = 100)) E0 = expectation_value(mps, H) #println(" = $real(E0)") @@ -148,9 +152,9 @@ end # find_groundstate!(mps2, H2, DMRG2(; maxiter=100, trscheme=trunc_st)) # E1 = expectation_value(mps2, H2) mps2 = InfiniteMPS(spacephys, space_virt) - Qs, Ps = mpo_compression(H, 0) + Qs, Ps = mpo_compression(H; truncbelow = 0) H2 = Qs - mps2, = find_groundstate(mps2, H2, VUMPS(;maxiter=100)) + mps2, = find_groundstate(mps2, H2, VUMPS(; maxiter = 100)) E0 = expectation_value(mps2, H2) end @@ -159,52 +163,62 @@ end D = 10 # Max bond dimension mps1 = InfiniteMPS(fill(ℂ^2, M), fill(ℂ^D, M)) - mps2 = InfiniteMPS(fill(Z2Space(0=>1, 1=>1), M), fill(Z2Space(0=>D, 1=>D), M)) - mps3 = InfiniteMPS(fill(Vect[FermionParity](0=>1, 1=>1), M), fill(Vect[FermionParity]( - 0=>D, 1=>D), M)) + mps2 = InfiniteMPS(fill(Z2Space(0 => 1, 1 => 1), M), fill(Z2Space(0 => D, 1 => D), M)) + mps3 = InfiniteMPS( + fill(Vect[FermionParity](0 => 1, 1 => 1), M), fill( + Vect[FermionParity]( + 0 => D, 1 => D + ), M + ) + ) H1 = transverse_field_ising_trivial_infinite(M) H2 = transverse_field_ising_z2_infinite(M) H3 = transverse_field_ising_fermion_parity_infinite(M) # Find groundstate of non-compressed Hamiltonian - mps1, = find_groundstate(mps1, H1, VUMPS(; maxiter=100)) + mps1, = find_groundstate(mps1, H1, VUMPS(; maxiter = 100)) E0_1 = expectation_value(mps1, H1) - mps2, = find_groundstate(mps2, H2, VUMPS(; maxiter=100)) + mps2, = find_groundstate(mps2, H2, VUMPS(; maxiter = 100)) E0_2 = expectation_value(mps2, H2) - mps3, = find_groundstate(mps3, H3, VUMPS(; maxiter=100)) + mps3, = find_groundstate(mps3, H3, VUMPS(; maxiter = 100)) E0_3 = expectation_value(mps3, H3) # Compress MPO mps1c = InfiniteMPS(fill(ℂ^2, M), fill(ℂ^D, M)) - mps2c = InfiniteMPS(fill(Z2Space(0=>1, 1=>1), M), fill(Z2Space(0=>D, 1=>D), M)) - mps3c = InfiniteMPS(fill(Vect[FermionParity](0=>1, 1=>1), M), fill(Vect[FermionParity]( - 0=>D, 1=>D), M)) - - H1c, _ = mpo_compression(H1, 10^-10) - H2c, _ = mpo_compression(H2, 10^-10) - H3c, _ = mpo_compression(H3, 10^-10) + mps2c = InfiniteMPS(fill(Z2Space(0 => 1, 1 => 1), M), fill(Z2Space(0 => D, 1 => D), M)) + mps3c = InfiniteMPS( + fill(Vect[FermionParity](0 => 1, 1 => 1), M), fill( + Vect[FermionParity]( + 0 => D, 1 => D + ), M + ) + ) + + H1c, _ = mpo_compression(H1; truncbelow = 10^-10) + H2c, _ = mpo_compression(H2; truncbelow = 10^-10) + H3c, _ = mpo_compression(H3; truncbelow = 10^-10) # Find groundstate of compressed Hamiltonian - mps1c, = find_groundstate(mps1c, H1c, VUMPS(; maxiter=100)) + mps1c, = find_groundstate(mps1c, H1c, VUMPS(; maxiter = 100)) E0_1c = expectation_value(mps1c, H1c) - mps2c, = find_groundstate(mps2c, H2c, VUMPS(; maxiter=100)) + mps2c, = find_groundstate(mps2c, H2c, VUMPS(; maxiter = 100)) E0_2c = expectation_value(mps2c, H2c) - mps3c, = find_groundstate(mps3c, H3c, VUMPS(; maxiter=100)) + mps3c, = find_groundstate(mps3c, H3c, VUMPS(; maxiter = 100)) E0_3c = expectation_value(mps3c, H3c) - # Assert that the ground state energies are equal up to a precision, display the + # Assert that the ground state energies are equal up to a precision, display the # compression @assert E0_1 ≈ E0_1c @assert E0_2 ≈ E0_2c @assert E0_3 ≈ E0_3c println("Max dim uncompressed H1: $(max_dim_ham(H1)) || Max dim compressed H1c: - $(max_dim_ham(H1c))" ) + $(max_dim_ham(H1c))") println("Max dim uncompressed H2: $(max_dim_ham(H2)) || Max dim compressed H2c: - $(max_dim_ham(H2c))" ) + $(max_dim_ham(H2c))") println("Max dim uncompressed H3: $(max_dim_ham(H3)) || Max dim compressed H3c: - $(max_dim_ham(H3c))" ) + $(max_dim_ham(H3c))") end @@ -212,31 +226,35 @@ end M = 10 D = 10 # Max bond dimension - mps1 = InfiniteMPS(fill(Vect[FermionParity](0=>1, 1=>1), M), - fill(Vect[FermionParity](0=>D, 1=>D), M)) + mps1 = InfiniteMPS( + fill(Vect[FermionParity](0 => 1, 1 => 1), M), + fill(Vect[FermionParity](0 => D, 1 => D), M) + ) H1 = kitaev_model_infinite(M) # Find groundstate of non-compressed Hamiltonian - mps1, = find_groundstate(mps1, H1, VUMPS(; maxiter=100)) + mps1, = find_groundstate(mps1, H1, VUMPS(; maxiter = 100)) E0_1 = expectation_value(mps1, H1) # Compress MPO - mps1c = InfiniteMPS(fill(Vect[FermionParity](0=>1, 1=>1), M), - fill(Vect[FermionParity](0=>D, 1=>D), M)) + mps1c = InfiniteMPS( + fill(Vect[FermionParity](0 => 1, 1 => 1), M), + fill(Vect[FermionParity](0 => D, 1 => D), M) + ) - H1c, _ = mpo_compression(H1, 10^-10) + H1c, _ = mpo_compression(H1; truncbelow = 10^-10) # Find groundstate of compressed Hamiltonian - mps1c, = find_groundstate(mps1c, H1c, VUMPS(; maxiter=100)) + mps1c, = find_groundstate(mps1c, H1c, VUMPS(; maxiter = 100)) E0_1c = expectation_value(mps1c, H1c) - # Assert that the ground state energies are equal up to a precision, display the + # Assert that the ground state energies are equal up to a precision, display the # compression @assert E0_1 ≈ E0_1c println("Max dim uncompressed H1: $(max_dim_ham(H1)) || Max dim compressed H1c: - $(max_dim_ham(H1c))" ) + $(max_dim_ham(H1c))") end @@ -244,51 +262,51 @@ end M = 10 D = 6 # Max bond dimension - mps1 = InfiniteMPS(fill(ℂ^3, M),fill(ℂ^D, M)) - mps2 = InfiniteMPS(fill(U1Space(0=>1, 1=>1, -1=>1), M), fill(U1Space(0=>D, 1=>D, -1=>D), M)) - mps3 = InfiniteMPS(fill(SU2Space(1 => 1), M), fill(SU2Space(1=>D), M)) + mps1 = InfiniteMPS(fill(ℂ^3, M), fill(ℂ^D, M)) + mps2 = InfiniteMPS(fill(U1Space(0 => 1, 1 => 1, -1 => 1), M), fill(U1Space(0 => D, 1 => D, -1 => D), M)) + mps3 = InfiniteMPS(fill(SU2Space(1 => 1), M), fill(SU2Space(1 => D), M)) H1 = heisenberg_XXX_trivial_infinite(M) H2 = heisenberg_XXX_U1_infinite(M) H3 = heisenberg_XXX_SU2_infinite(M) # Find groundstate of non-compressed Hamiltonian - mps1, = find_groundstate(mps1, H1, VUMPS(; maxiter=100)) + mps1, = find_groundstate(mps1, H1, VUMPS(; maxiter = 100)) E0_1 = expectation_value(mps1, H1) - mps2, = find_groundstate(mps2, H2, VUMPS(; maxiter=100)) + mps2, = find_groundstate(mps2, H2, VUMPS(; maxiter = 100)) E0_2 = expectation_value(mps2, H2) - mps3, = find_groundstate(mps3, H3, VUMPS(; maxiter=100)) + mps3, = find_groundstate(mps3, H3, VUMPS(; maxiter = 100)) E0_3 = expectation_value(mps3, H3) # Compress MPO - mps1c = InfiniteMPS(fill(ℂ^3, M),fill(ℂ^D, M)) - mps2c = InfiniteMPS(fill(U1Space(0=>1, 1=>1, -1=>1), M), fill(U1Space(0=>D, 1=>D, -1=>D), M)) - mps3c = InfiniteMPS(fill(SU2Space(1 => 1), M), fill(SU2Space(1=>D), M)) + mps1c = InfiniteMPS(fill(ℂ^3, M), fill(ℂ^D, M)) + mps2c = InfiniteMPS(fill(U1Space(0 => 1, 1 => 1, -1 => 1), M), fill(U1Space(0 => D, 1 => D, -1 => D), M)) + mps3c = InfiniteMPS(fill(SU2Space(1 => 1), M), fill(SU2Space(1 => D), M)) - H1c, _ = mpo_compression(H1, 10^-10) - H2c, _ = mpo_compression(H2, 10^-10) - H3c, _ = mpo_compression(H3, 10^-10) + H1c, _ = mpo_compression(H1; truncbelow = 10^-10) + H2c, _ = mpo_compression(H2; truncbelow = 10^-10) + H3c, _ = mpo_compression(H3; truncbelow = 10^-10) # Find groundstate of compressed Hamiltonian - mps1c, = find_groundstate(mps1c, H1c, VUMPS(; maxiter=100)) + mps1c, = find_groundstate(mps1c, H1c, VUMPS(; maxiter = 100)) E0_1c = expectation_value(mps1c, H1c) - mps2c, = find_groundstate(mps2c, H2c, VUMPS(; maxiter=100)) + mps2c, = find_groundstate(mps2c, H2c, VUMPS(; maxiter = 100)) E0_2c = expectation_value(mps2c, H2c) - mps3c, = find_groundstate(mps3c, H3c, VUMPS(; maxiter=100)) + mps3c, = find_groundstate(mps3c, H3c, VUMPS(; maxiter = 100)) E0_3c = expectation_value(mps3c, H3c) - # Assert that the ground state energies are equal up to a precision, display the + # Assert that the ground state energies are equal up to a precision, display the # compression @assert E0_1 ≈ E0_1c @assert E0_2 ≈ E0_2c @assert E0_3 ≈ E0_3c println("Max dim uncompressed H1: $(max_dim_ham(H1)) || Max dim compressed H1c: - $(max_dim_ham(H1c))" ) + $(max_dim_ham(H1c))") println("Max dim uncompressed H2: $(max_dim_ham(H2)) || Max dim compressed H2c: - $(max_dim_ham(H2c))" ) + $(max_dim_ham(H2c))") println("Max dim uncompressed H3: $(max_dim_ham(H3)) || Max dim compressed H3c: - $(max_dim_ham(H3c))" ) + $(max_dim_ham(H3c))") end @@ -296,30 +314,30 @@ end M = 10 D = 6 # Max bond dimension - mps1 = InfiniteMPS(fill(ℂ^3, M),fill(ℂ^D, M)) + mps1 = InfiniteMPS(fill(ℂ^3, M), fill(ℂ^D, M)) H1 = bilinear_biquadratic_model_trivial_infinite(M) # Find groundstate of non-compressed Hamiltonian - mps1, = find_groundstate(mps1, H1, VUMPS(; maxiter=100)) + mps1, = find_groundstate(mps1, H1, VUMPS(; maxiter = 100)) E0_1 = expectation_value(mps1, H1) # Compress MPO - mps1c = InfiniteMPS(fill(ℂ^3, M),fill(ℂ^D, M)) + mps1c = InfiniteMPS(fill(ℂ^3, M), fill(ℂ^D, M)) - H1c, _ = mpo_compression(H1, 10^-10) + H1c, _ = mpo_compression(H1; truncbelow = 10^-10) # Find groundstate of compressed Hamiltonian - mps1c, = find_groundstate(mps1c, H1c, VUMPS(; maxiter=100)) + mps1c, = find_groundstate(mps1c, H1c, VUMPS(; maxiter = 100)) E0_1c = expectation_value(mps1c, H1c) - # Assert that the ground state energies are equal up to a precision, display the + # Assert that the ground state energies are equal up to a precision, display the # compression @assert E0_1 ≈ E0_1c println("Max dim uncompressed H1: $(max_dim_ham(H1)) || Max dim compressed H1c: - $(max_dim_ham(H1c))" ) + $(max_dim_ham(H1c))") end @@ -328,39 +346,39 @@ end D = 6 # Max bond dimension mps1 = InfiniteMPS(fill(ℂ^3, M), fill(ℂ^D, M)) - mps2 = InfiniteMPS(fill(Z3Space(0=>1, 1=>1, 2=>1), M), fill(Z3Space(0=>1, 1=>1, 2=>1), M)) + mps2 = InfiniteMPS(fill(Z3Space(0 => 1, 1 => 1, 2 => 1), M), fill(Z3Space(0 => 1, 1 => 1, 2 => 1), M)) H1 = quantum_potts_trivial_infinite(M) H2 = quantum_potts_Z3_infinite(M) # Find groundstate of non-compressed Hamiltonian - mps1, = find_groundstate(mps1, H1, VUMPS(; maxiter=100)) + mps1, = find_groundstate(mps1, H1, VUMPS(; maxiter = 100)) E0_1 = expectation_value(mps1, H1) - mps2, = find_groundstate(mps2, H2, VUMPS(; maxiter=100)) + mps2, = find_groundstate(mps2, H2, VUMPS(; maxiter = 100)) E0_2 = expectation_value(mps2, H2) # Compress MPO mps1c = InfiniteMPS(fill(ℂ^3, M), fill(ℂ^D, M)) - mps2c = InfiniteMPS(fill(Z3Space(0=>1, 1=>1, 2=>1), M), fill(Z3Space(0=>1, 1=>1, 2=>1), M)) + mps2c = InfiniteMPS(fill(Z3Space(0 => 1, 1 => 1, 2 => 1), M), fill(Z3Space(0 => 1, 1 => 1, 2 => 1), M)) - H1c, _ = mpo_compression(H1, 10^-10) - H2c, _ = mpo_compression(H2, 10^-10) + H1c, _ = mpo_compression(H1; truncbelow = 10^-10) + H2c, _ = mpo_compression(H2; truncbelow = 10^-10) # Find groundstate of compressed Hamiltonian - mps1c, = find_groundstate(mps1c, H1c, VUMPS(; maxiter=100)) + mps1c, = find_groundstate(mps1c, H1c, VUMPS(; maxiter = 100)) E0_1c = expectation_value(mps1c, H1c) - mps2c, = find_groundstate(mps2c, H2c, VUMPS(; maxiter=100)) + mps2c, = find_groundstate(mps2c, H2c, VUMPS(; maxiter = 100)) E0_2c = expectation_value(mps2c, H2c) - # Assert that the ground state energies are equal up to a precision, display the + # Assert that the ground state energies are equal up to a precision, display the # compression @assert E0_1 ≈ E0_1c @assert E0_2 ≈ E0_2c println("Max dim uncompressed H1: $(max_dim_ham(H1)) || Max dim compressed H1c: - $(max_dim_ham(H1c))" ) + $(max_dim_ham(H1c))") println("Max dim uncompressed H2: $(max_dim_ham(H2)) || Max dim compressed H2c: - $(max_dim_ham(H2c))" ) + $(max_dim_ham(H2c))") end @@ -373,16 +391,16 @@ end L = 10 # number of sites range = 4 # cutoff for the maximum range between interactions - mps = InfiniteMPS(fill(Z2Space(0=>1, 1=>1),L), fill(Z2Space(0=>D, 1=>D),L)) + mps = InfiniteMPS(fill(Z2Space(0 => 1, 1 => 1), L), fill(Z2Space(0 => D, 1 => D), L)) H = create_long_range_ising_symmetries_infinite_msite_random(L, range, 0.1, 0.001) - mps, = find_groundstate(mps, H, VUMPS(; maxiter=100)) + mps, = find_groundstate(mps, H, VUMPS(; maxiter = 100)) E0 = expectation_value(mps, H) #println(" = $real(E0)") - mps2 = InfiniteMPS(fill(Z2Space(0=>1, 1=>1),L), fill(Z2Space(0=>D, 1=>D),L)) - H2, Rs = mpo_compression(H, 0) - mps2, = find_groundstate(mps2, H2, VUMPS(; maxiter=100)) + mps2 = InfiniteMPS(fill(Z2Space(0 => 1, 1 => 1), L), fill(Z2Space(0 => D, 1 => D), L)) + H2, Rs = mpo_compression(H; truncbelow = 0) + mps2, = find_groundstate(mps2, H2, VUMPS(; maxiter = 100)) E1 = expectation_value(mps2, H2) @assert E0 ≈ E1 -end \ No newline at end of file +end