From b2cd36c6e03682175afef1da917224dfbb557f66 Mon Sep 17 00:00:00 2001 From: TymoteuszTula Date: Mon, 13 Jul 2026 14:00:26 +0200 Subject: [PATCH 1/8] 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 6ea720fb6..671f81b5f 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" @@ -39,6 +40,7 @@ Accessors = "0.1" Adapt = "4" BlockTensorKit = "0.3.19" 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 ea209f449..be82b29d8 100644 --- a/src/MPSKit.jl +++ b/src/MPSKit.jl @@ -52,6 +52,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 @@ -84,13 +85,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, NoTimerOutput, @timeit, timeit, print_timer +using Dictionaries # TODO: change mpocompression to not use Dictionaries # Includes # -------- @@ -137,6 +139,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 9b5e571b11009cc0183210f3a8900ed66f3d6af6 Mon Sep 17 00:00:00 2001 From: herviou Date: Wed, 15 Jul 2026 13:53:44 +0200 Subject: [PATCH 2/8] 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 eaf9ecb4741a1972eac1e707f1210cddd20e4ffa Mon Sep 17 00:00:00 2001 From: herviou Date: Wed, 15 Jul 2026 14:26:39 +0200 Subject: [PATCH 3/8] Actually fixing formatting. --- test/operators/finite_mpo_compression.jl | 152 +++++++-------- test/operators/infinite_mpo_compression.jl | 206 +++++++++++---------- 2 files changed, 190 insertions(+), 168 deletions(-) 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 From 716f50ad56782cea572445ce4078ae2e7077d62e Mon Sep 17 00:00:00 2001 From: lkdvos Date: Mon, 12 Oct 2026 11:43:48 -0400 Subject: [PATCH 4/8] Infinite MPOHamiltonian compression via `SvdCut` Replace the standalone `mpo_compression` implementation by an integration with the existing `left_canonicalize!`/`right_canonicalize!` and `SvdCut` machinery: - factor the finite canonicalization steps into `_left_orth_jordan`/`_absorb_left_gauge` and their mirrors, which are shared with the infinite case - `left_canonicalize!`/`right_canonicalize!(::InfiniteMPOHamiltonian)` iterate rank-revealing block QR/LQ decompositions in place until the gauge transformation is unitary - `changebonds(::InfiniteMPOHamiltonian, ::SvdCut)` implements Algorithm 5 of Parker et al., truncating the singular values of the gauge between the right- and left-canonical forms This also fixes the canonical form of (finite) MPOHamiltonians when `C`/`B` contain identity components: the projection onto the identity was missing a `1 / dim(P)` normalization, and rank-deficient blocks could leave directions in `Q` that overlap with the identity. Both are resolved by decomposing the full block with the identity as its first column/row. Co-Authored-By: Claude Opus 5.5 --- Project.toml | 2 - docs/src/changelog.md | 3 + src/MPSKit.jl | 5 +- src/algorithms/changebonds/svdcut.jl | 52 ++ src/operators/mpocompression.jl | 636 --------------------- src/operators/ortho.jl | 344 +++++++---- test/algorithms/changebonds.jl | 31 + test/hamiltonian/finite.jl | 36 ++ test/hamiltonian/infinite.jl | 33 ++ test/operators/finite_mpo_compression.jl | 379 ------------ test/operators/infinite_mpo_compression.jl | 406 ------------- 11 files changed, 386 insertions(+), 1541 deletions(-) delete mode 100644 src/operators/mpocompression.jl delete mode 100644 test/operators/finite_mpo_compression.jl delete mode 100644 test/operators/infinite_mpo_compression.jl diff --git a/Project.toml b/Project.toml index 671f81b5f..6ea720fb6 100644 --- a/Project.toml +++ b/Project.toml @@ -10,7 +10,6 @@ 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" @@ -40,7 +39,6 @@ Accessors = "0.1" Adapt = "4" BlockTensorKit = "0.3.19" 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/docs/src/changelog.md b/docs/src/changelog.md index 6005693f5..e6523fd19 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -24,6 +24,9 @@ When releasing a new version, move the "Unreleased" changes to a new version sec - `ClusterExpansion` time-evolution algorithm for `make_time_mpo`, constructing evolution MPOs for finite and infinite nearest-neighbor Hamiltonians through a nonperturbative cluster expansion that matches exact exponentials on clusters of up to `N` sites. ([#252](https://github.com/QuantumKitHub/MPSKit.jl/pull/252)) +- `changebonds(::InfiniteMPOHamiltonian, ::SvdCut)` compresses infinite Hamiltonians with the + almost-Schmidt decomposition of [Parker et al.](https://arxiv.org/abs/1909.06341), and + `left_canonicalize!`/`right_canonicalize!` bring an `InfiniteMPOHamiltonian` into canonical form. ### Changed diff --git a/src/MPSKit.jl b/src/MPSKit.jl index be82b29d8..ea209f449 100644 --- a/src/MPSKit.jl +++ b/src/MPSKit.jl @@ -52,7 +52,6 @@ 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 @@ -85,14 +84,13 @@ using Accessors using HalfIntegers using DocStringExtensions -using LinearAlgebra: diag, Diagonal, diagm +using LinearAlgebra: diag, Diagonal using LinearAlgebra: LinearAlgebra using Random using Base: @kwdef, @propagate_inbounds using LoggingExtras using OhMyThreads using TimerOutputs: TimerOutput, NoTimerOutput, @timeit, timeit, print_timer -using Dictionaries # TODO: change mpocompression to not use Dictionaries # Includes # -------- @@ -139,7 +137,6 @@ 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/algorithms/changebonds/svdcut.jl b/src/algorithms/changebonds/svdcut.jl index 79ad973ba..ee8ef9fdb 100644 --- a/src/algorithms/changebonds/svdcut.jl +++ b/src/algorithms/changebonds/svdcut.jl @@ -178,3 +178,55 @@ function changebonds!(H::FiniteMPOHamiltonian, alg::SvdCut) return H end + +changebonds(H::InfiniteMPOHamiltonian, alg::SvdCut) = changebonds!(copy(H), alg) +function changebonds!(H::InfiniteMPOHamiltonian, alg::SvdCut) + channels_before = [dim(right_virtualspace(H[i])) - 2 for i in 1:length(H)] + + # the interaction channels of the gauge between the right- and left-canonical forms + right_canonicalize!(H) + _, Cs = uniform_leftorth!(H) + + Us = map(C -> first(svd_trunc(C; trunc = alg.trunc, alg = alg.alg_svd)), Cs) + for i in 1:length(H) + H[i] = _project_channels(Us[i - 1], H[i], Us[i]) + end + + emptied = findall( + i -> channels_before[i] > 0 && dim(right_virtualspace(H[i])) - 2 == 0, + eachindex(channels_before) + ) + isempty(emptied) || @warn "`SvdCut` discarded every interaction crossing bond(s) $emptied; those terms are gone from the compressed operator. Loosen `trunc` or check the scale of the input Hamiltonian." + + return H +end + +# Compute `diag(1, Uₗ', 1) * W * diag(1, Uᵣ, 1)`. +function _project_channels(Uₗ, W::JordanMPOTensor, Uᵣ) + P = physicalspace(W) + Vₗ, Vᵣ = _channelspace(Uₗ), _channelspace(Uᵣ) + WA, WB, WC = W.A, W.B, W.C + Uₗᴴ = adjoint(Uₗ) + + if dim(Vₗ) == 0 || dim(Vᵣ) == 0 + A = similar(WA, Vₗ ⊗ P ← P ⊗ Vᵣ) + else + @plansor A[l p; p' r] := Uₗᴴ[l; l'] * WA[l' p; p' r'] * Uᵣ[r'; r] + end + if dim(Vₗ) == 0 + B = similar(WB, Vₗ ⊗ P ← P) + else + @plansor B[l p; p'] := Uₗᴴ[l; l'] * WB[l' p; p'] + end + if dim(Vᵣ) == 0 + C = similar(WC, P ← P ⊗ Vᵣ) + else + @plansor C[p; p' r] := WC[p; p' r'] * Uᵣ[r'; r] + end + + unit = rightunitspace(P) + V = ⊞(unit, Vₗ, unit) ⊗ P ← P ⊗ ⊞(unit, Vᵣ, unit) + return JordanMPOTensor(V, A, B, C, W.D) +end + +_channelspace(U) = dim(U) == 0 ? SumSpace{spacetype(U)}() : domain(U)[1] diff --git a/src/operators/mpocompression.jl b/src/operators/mpocompression.jl deleted file mode 100644 index 3867702d1..000000000 --- a/src/operators/mpocompression.jl +++ /dev/null @@ -1,636 +0,0 @@ -""" - 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. - -### Arguments -- `H`: Hamiltonian to be compressed. Either FiniteMPOHamiltonian or InfiniteMPOHamiltonian. - -### Keyword arguments -- `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; 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; 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; truncbelow::Number = 10^-8, verbosity::Int = 0) - 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, 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 - 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 - 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; 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 - ) - for n in 1:N - Hnew[n] = SparseBlockTensorMap(reduce_blocks_mpo(H[n])) - end - - 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 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 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] - # 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 - 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] -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) - return 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, 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)] - U, S, V′ = svd_compact(M2) - - num_sectors = length(blocksectors(S)) - new_dims = Vector{Int32}(undef, num_sectors) - - 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 ≤ truncbelow - new_dims[n] -= 1 - end - end - end - end - - fdims = Vector{Pair{sectortype(S), Int32}}() - 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′)) - - # 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]) - - 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}; - 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 - - ε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]) - 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]) - 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 - if verbosity > 1 - printstyled("| total err = $(sum(εs) / N)\n", color = :white) - end - end - 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}; - 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]) - 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]) - 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 - 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 - return Rs, H -end diff --git a/src/operators/ortho.jl b/src/operators/ortho.jl index 9ea1f64a5..c164e6024 100644 --- a/src/operators/ortho.jl +++ b/src/operators/ortho.jl @@ -2,75 +2,191 @@ _dense(t::SparseBlockTensorMap) = BlockTensorMap(t) _dense(t) = t +# Finite MPOHamiltonians +# ---------------------- + function left_canonicalize!( H::FiniteMPOHamiltonian, i::Int; alg::MatrixAlgebraKit.AbstractAlgorithm = Defaults.alg_orth() ) 1 ≤ i < length(H) || throw(ArgumentError("Bounds error in canonicalize")) + H[i], t, R = _left_orth_jordan(H[i]; alg) + H[i + 1] = _absorb_left_gauge(H[i], t, R, H[i + 1]) + return H +end + +function right_canonicalize!( + H::FiniteMPOHamiltonian, i::Int; + alg::MatrixAlgebraKit.AbstractAlgorithm = Defaults.alg_orth() + ) + 1 < i ≤ length(H) || throw(ArgumentError("Bounds error in canonicalize")) + R, t, H[i] = _right_orth_jordan(H[i]; alg) + H[i - 1] = _absorb_right_gauge(H[i - 1], R, t, H[i]) + return H +end + +# Infinite MPOHamiltonians +# ------------------------ + +""" + left_canonicalize!(H::InfiniteMPOHamiltonian; kwargs...) -> H + +Bring `H` into left-canonical form, where the upper-left block `[1 C; 0 A]` of every site has +orthonormal columns, by iterating block-respecting QR decompositions until the gauge transformation +is unitary. Redundant interaction channels are removed along the way. + +# Keyword arguments + +- `tol = Defaults.tolgauge`: convergence tolerance on the deviation of the gauge transformation from + a unitary one, which is also the relative tolerance for discarding redundant channels +- `maxiter = Defaults.maxiter`: maximal number of sweeps over the unit cell +- `verbosity = VERBOSE_WARN`: setting for how much information is displayed + +# References + +* [Parker et al. Phys. Rev. B 102 (2020)](@cite parker2020) +""" +function left_canonicalize!(H::InfiniteMPOHamiltonian; kwargs...) + uniform_leftorth!(H; kwargs...) + return H +end + +""" + right_canonicalize!(H::InfiniteMPOHamiltonian; kwargs...) -> H + +Bring `H` into right-canonical form, where the lower-right block `[A B; 0 1]` of every site has +orthonormal rows, by iterating block-respecting LQ decompositions until the gauge transformation +is unitary. Redundant interaction channels are removed along the way. + +See [`left_canonicalize!`](@ref) for the keyword arguments. +""" +function right_canonicalize!(H::InfiniteMPOHamiltonian; kwargs...) + uniform_rightorth!(H; kwargs...) + return H +end + +# Bring `H` into left-canonical form in place, and return the middle blocks `L` of the accumulated +# gauge transformations `diag(1, L[i - 1], 1) * H₀[i] = H[i] * diag(1, L[i], 1)` up to the gauging +# away of identity components. +function uniform_leftorth!( + H::InfiniteMPOHamiltonian; + tol::Real = Defaults.tolgauge, maxiter::Int = Defaults.maxiter, verbosity::Int = VERBOSE_WARN + ) + alg_trunc = MatrixAlgebraKit.TruncatedAlgorithm(Defaults.alg_svd(), trunctol(; rtol = tol)) + alg = MatrixAlgebraKit.LeftOrthViaSVD(alg_trunc) + Ls = PeriodicVector(Vector{Any}(nothing, length(H))) + return with_verbosity(; verbosity) do + log = IterLog("MPO LC") + @log_initialization loginit!(log, Inf) + for iter in 1:maxiter + local t, R + for i in 1:length(H) + H[i], t, R = _left_orth_jordan(H[i]; alg) + H[i + 1] = _absorb_left_gauge(H[i], t, R, H[i + 1]) + Ls[i] = isnothing(Ls[i]) || dim(R) == 0 ? R : R * Ls[i] + end + ϵ = _gauge_distance(t, R) + if ϵ <= tol + @log_convergence logfinish!(log, iter, ϵ) + break + elseif iter == maxiter + @log_nonconvergence logcancel!(log, iter, ϵ) + else + @log_iteration logiter!(log, iter, ϵ) + end + end + return H, map(identity, Ls) + end +end + +# Bring `H` into right-canonical form in place. +function uniform_rightorth!( + H::InfiniteMPOHamiltonian; + tol::Real = Defaults.tolgauge, maxiter::Int = Defaults.maxiter, verbosity::Int = VERBOSE_WARN + ) + alg_trunc = MatrixAlgebraKit.TruncatedAlgorithm(Defaults.alg_svd(), trunctol(; rtol = tol)) + alg = MatrixAlgebraKit.RightOrthViaSVD(alg_trunc) + return with_verbosity(; verbosity) do + log = IterLog("MPO RC") + @log_initialization loginit!(log, Inf) + for iter in 1:maxiter + local t, R + for i in length(H):-1:1 + R, t, H[i] = _right_orth_jordan(H[i]; alg) + H[i - 1] = _absorb_right_gauge(H[i - 1], R, t, H[i]) + end + ϵ = _gauge_distance(t, R) + if ϵ <= tol + @log_convergence logfinish!(log, iter, ϵ) + break + elseif iter == maxiter + @log_nonconvergence logcancel!(log, iter, ϵ) + else + @log_iteration logiter!(log, iter, ϵ) + end + end + return H + end +end - W = H[i] +# distance of the gauge transformation `[1 t 0; 0 R 0; 0 0 1]` (or its mirror) to a unitary one +function _gauge_distance(t, R) + codomain(R) == domain(R) || return Inf + return norm(t) + norm(R' * R - id(storagetype(R), domain(R))) +end + +# Implementation +# -------------- + +# Decompose `W = Q * G` with `G = [1 t 0; 0 R 0; 0 0 1]` and `Q` left-canonical. +function _left_orth_jordan(W::JordanMPOTensor; alg = Defaults.alg_orth()) P = physicalspace(W) _rightunit = rightunitspace(P) S = spacetype(W) d = sqrt(dim(P)) + WA, WB, WD = W.A, W.B, W.D - # block accessors recompute a fresh `SparseBlockTensorMap` on every access, so bind - # them once and reuse the locals throughout - WA, WB, WC, WD = W.A, W.B, W.C, W.D - - # orthogonalize second column against first - Wi = W[1, 1, 1, 1] - WI = removeunit(Wi, 1) - @plansor t[l; r] := conj(WI[p; p' l]) * WC[p; p' r] - # TODO: the following is currently broken due to a TensorKit bug - # @plansor C′[p; p' r] := WC[p; p' r] - WI[p; p' l] * t[l; r] - @plansor C′[p; p' r] := -WI[p; p' l] * t[l; r] - add!(C′, _dense(WC)) - - # QR of second column - if size(W, 1) == 1 - tmp = transpose(C′, ((2, 1), (3,))) - Q, R = left_orth!(tmp; alg) - - if dim(R) == 0 # fully truncated - V = _rightunit ⊞ _rightunit - Q1 = typeof(WA)(undef, SumSpace{S}() ⊗ P ← P ⊗ SumSpace{S}()) - Q2 = typeof(WC)(undef, P ← P ⊗ SumSpace{S}()) - else - V = ⊞(_rightunit, space(R, 1), _rightunit) - scale!(Q, d) - scale!(R, inv(d)) - Q1 = typeof(WA)(undef, SumSpace{S}() ⊗ P ← P ⊗ space(R, 1)) - Q2 = transpose(Q, ((2,), (1, 3))) - end - H[i] = JordanMPOTensor(codomain(W) ← physicalspace(W) ⊗ V, Q1, WB, Q2, WD) + # QR of the upper-left block `[1 C; 0 A]` with the identity as first column, which keeps all other + # columns of `Q` orthogonal to it, even if the block is rank-deficient + rows = 1:max(size(W, 1) - 1, 1) + V = BlockTensorMap(SparseBlockTensorMap(W)[rows, 1, 1, 1:(end - 1)]) + Q, R = left_orth!(transpose(V, ((3, 1, 2), (4,))); alg = Defaults.alg_orth()) + Π = _split_unit(domain(Q)[1], _rightunit, storagetype(W)) + Q = (Q * Π')[:, :, :, 2:2] + R = Π * R + t = scale!(R[1:1, 2:end], inv(d)) + R = R[2:2, 2:end] + + # other (possibly truncating) algorithms only act on the interaction channels + if alg != Defaults.alg_orth() + U, R = left_orth!(R; alg) + Q = Q * U + end + scale!(R, inv(d)) + + if dim(R) == 0 # fully truncated + V = _rightunit ⊞ _rightunit + Q1 = typeof(WA)(undef, codomain(WA) ← P ⊗ SumSpace{S}()) + Q2 = typeof(W.C)(undef, P ← P ⊗ SumSpace{S}()) else - tmp = transpose(cat(insertleftunit(C′, 1), WA; dims = 1), ((3, 1, 2), (4,))) - Q, R = left_orth!(tmp; alg) - - if dim(R) == 0 # fully truncated - V = _rightunit ⊞ _rightunit - Q1 = typeof(WA)(undef, codomain(WA) ← P ⊗ SumSpace{S}()) - Q2 = typeof(WC)(undef, P ← P ⊗ SumSpace{S}()) - else - scale!(Q, d) - scale!(R, inv(d)) - Q′ = transpose(Q, ((2, 3), (1, 4))) - Q1 = Q′[2:end, 1, 1, 1] - Q2 = removeunit(SparseBlockTensorMap(Q′[1:1, 1, 1, 1]), 1) - V = ⊞(_rightunit, right_virtualspace(Q′), _rightunit) - end - H[i] = JordanMPOTensor(codomain(W) ← P ⊗ V, Q1, WB, Q2, WD) + Q = scale!(transpose(Q, ((2, 3), (1, 4))), d) + Q1 = size(W, 1) == 1 ? typeof(WA)(undef, SumSpace{S}() ⊗ P ← P ⊗ right_virtualspace(Q)) : + SparseBlockTensorMap(Q[2:end, 1, 1, 1]) + Q2 = removeunit(SparseBlockTensorMap(Q[1:1, 1, 1, 1]), 1) + V = ⊞(_rightunit, right_virtualspace(Q), _rightunit) end + return JordanMPOTensor(codomain(W) ← P ⊗ V, Q1, WB, Q2, WD), t, R +end - # absorb into next site - W′ = H[i + 1] +# Compute `G * W′` for the gauge `G = [1 t 0; 0 R 0; 0 0 1]` obtained from `_left_orth_jordan`, +# where `Q` is the tensor to the left of `W′`. +function _absorb_left_gauge(Q::JordanMPOTensor, t, R, W′::JordanMPOTensor) W′A, W′B, W′C, W′D = W′.A, W′.B, W′.C, W′.D if size(W′, 4) > 1 && dim(R) != 0 @plansor A′[l p; p' r] := R[l; r'] * W′A[r' p; p' r] else - A′ = similar(W′A, right_virtualspace(H[i].A) ⊗ physicalspace(W′) ← domain(W′A)) + A′ = similar(W′A, right_virtualspace(Q.A) ⊗ physicalspace(W′) ← domain(W′A)) end if size(W′, 4) > 1 @@ -83,88 +199,67 @@ function left_canonicalize!( if dim(R) != 0 @plansor B′[l p; p'] := R[l; r] * W′B[r p; p'] else - B′ = similar(W′B, right_virtualspace(H[i].A) ⊗ physicalspace(W′) ← domain(W′B)) + B′ = similar(W′B, right_virtualspace(Q.A) ⊗ physicalspace(W′) ← domain(W′B)) end @plansor D′[l p; p'] := t[l; r] * W′B[r p; p'] D′ = add!(removeunit(D′, 1), _dense(W′D)) - H[i + 1] = JordanMPOTensor( - right_virtualspace(H[i]) ⊗ physicalspace(W′) ← domain(W′), - A′, B′, C′, D′ + return JordanMPOTensor( + right_virtualspace(Q) ⊗ physicalspace(W′) ← domain(W′), A′, B′, C′, D′ ) - return H end -function right_canonicalize!( - H::FiniteMPOHamiltonian, i::Int; - alg::MatrixAlgebraKit.AbstractAlgorithm = Defaults.alg_orth() - ) - 1 < i ≤ length(H) || throw(ArgumentError("Bounds error in canonicalize")) - - W = H[i] +# Decompose `W = G * Q` with `G = [1 0 0; 0 R t; 0 0 1]` and `Q` right-canonical. +function _right_orth_jordan(W::JordanMPOTensor; alg = Defaults.alg_orth()) P = physicalspace(W) _rightunit = rightunitspace(P) S = spacetype(W) d = sqrt(dim(P)) + WA, WC, WD = W.A, W.C, W.D - # block accessors recompute a fresh `SparseBlockTensorMap` on every access, so bind - # them once and reuse the locals throughout - WA, WB, WC, WD = W.A, W.B, W.C, W.D - - # orthogonalize second row against last - Wi = W[end, 1, 1, end] - WI = removeunit(Wi, 4) - @plansor t[l; r] := conj(WI[r p; p']) * WB[l p; p'] - # TODO: the following is currently broken due to a TensorKit bug - # @plansor B′[l p; p'] := WB[l p; p'] - WI[r p; p'] * t[l; r] - @plansor B′[l p; p'] := -WI[r p; p'] * t[l; r] - add!(B′, _dense(WB)) - - # LQ of second row - if size(W, 4) == 1 - tmp = transpose(B′, ((1,), (3, 2))) - R, Q = right_orth!(tmp; alg) - - if dim(R) == 0 - V = _rightunit ⊞ _rightunit - Q1 = typeof(WA)(undef, SumSpace{S}() ⊗ P ← P ⊗ SumSpace{S}()) - Q2 = typeof(WB)(undef, SumSpace{S}() ⊗ P ← P) - else - V = ⊞(_rightunit, space(Q, 1), _rightunit) - scale!(Q, d) - scale!(R, inv(d)) - Q1 = typeof(WA)(undef, space(Q, 1) ⊗ P ← P ⊗ SumSpace{S}()) - Q2 = transpose(Q, ((1, 3), (2,))) - end - H[i] = JordanMPOTensor(V ⊗ P ← domain(W), Q1, Q2, WC, WD) + # LQ of the lower-right block `[A B; 0 1]` with the identity as first row, which keeps all other + # rows of `Q` orthogonal to it, even if the block is rank-deficient + cols = min(2, size(W, 4)):size(W, 4) + V = SparseBlockTensorMap(W)[2:end, 1, 1, cols] + V = BlockTensorMap(cat(V[end:end, 1, 1, :], V[1:(end - 1), 1, 1, :]; dims = 1)) + R, Q = right_orth!(transpose(V, ((1,), (3, 4, 2))); alg = Defaults.alg_orth()) + Π = _split_unit(codomain(Q)[1], _rightunit, storagetype(W)) + Q = (Π * Q)[2:2, :, :, :] + R = R * Π' + t = scale!(R[2:end, 1:1], inv(d)) + R = R[2:end, 2:2] + + # other (possibly truncating) algorithms only act on the interaction channels + if alg != Defaults.alg_orth() + R, U = right_orth!(R; alg) + Q = U * Q + end + scale!(R, inv(d)) + + if dim(R) == 0 # fully truncated + V = _rightunit ⊞ _rightunit + Q1 = typeof(WA)(undef, SumSpace{S}() ⊗ P ← domain(WA)) + Q2 = typeof(W.B)(undef, SumSpace{S}() ⊗ P ← P) else - B′′ = insertleftunit(B′, 4) - tmp = transpose(cat(B′′, WA; dims = 4), ((1,), (3, 4, 2))) - R, Q = right_orth!(tmp; alg) - if dim(R) == 0 - V = _rightunit ⊞ _rightunit - Q1 = typeof(WA)(undef, SumSpace{S}() ⊗ P ← domain(WA)) - Q2 = typeof(WB)(undef, SumSpace{S}() ⊗ P ← P) - else - scale!(Q, d) - scale!(R, inv(d)) - Q′ = transpose(Q, ((1, 4), (2, 3))) - Q1 = SparseBlockTensorMap(Q′[1, 1, 1, 2:end]) - Q2 = removeunit(SparseBlockTensorMap(Q′[1, 1, 1, 1:1]), 4) - V = ⊞(_rightunit, left_virtualspace(Q′), _rightunit) - end - H[i] = JordanMPOTensor(V ⊗ P ← domain(W), Q1, Q2, WC, WD) + Q = scale!(transpose(Q, ((1, 4), (2, 3))), d) + Q1 = size(W, 4) == 1 ? typeof(WA)(undef, left_virtualspace(Q) ⊗ P ← P ⊗ SumSpace{S}()) : + SparseBlockTensorMap(Q[1, 1, 1, 1:(end - 1)]) + Q2 = removeunit(SparseBlockTensorMap(Q[1, 1, 1, end:end]), 4) + V = ⊞(_rightunit, left_virtualspace(Q), _rightunit) end + return R, t, JordanMPOTensor(V ⊗ P ← domain(W), Q1, Q2, WC, WD) +end - # absorb into previous site - W′ = H[i - 1] +# Compute `W′ * G` for the gauge `G = [1 0 0; 0 R t; 0 0 1]` obtained from `_right_orth_jordan`, +# where `Q` is the tensor to the right of `W′`. +function _absorb_right_gauge(W′::JordanMPOTensor, R, t, Q::JordanMPOTensor) W′A, W′B, W′C, W′D = W′.A, W′.B, W′.C, W′.D if size(W′, 1) > 1 && dim(R) != 0 @plansor A′[l p; p' r] := W′A[l p; p' r'] * R[r'; r] else - A′ = similar(W′A, codomain(W′A) ← physicalspace(W′A) ⊗ left_virtualspace(H[i].A)) + A′ = similar(W′A, codomain(W′A) ← physicalspace(W′A) ⊗ left_virtualspace(Q.A)) end if size(W′, 1) > 1 @@ -177,11 +272,32 @@ function right_canonicalize!( if dim(R) != 0 @plansor C′[p; p' r] := W′C[p; p' r'] * R[r'; r] else - C′ = similar(W′C, codomain(W′C) ← physicalspace(W′) ⊗ left_virtualspace(H[i].A)) + C′ = similar(W′C, codomain(W′C) ← physicalspace(W′) ⊗ left_virtualspace(Q.A)) end @plansor D′[p; p' r] := W′C[p; p' r'] * t[r'; r] D′ = add!(removeunit(D′, 3), _dense(W′D)) - H[i - 1] = JordanMPOTensor(codomain(W′) ← physicalspace(W′) ⊗ V, A′, B′, C′, D′) - return H + return JordanMPOTensor( + codomain(W′) ← physicalspace(W′) ⊗ left_virtualspace(Q), A′, B′, C′, D′ + ) +end + +# Isometry that splits the first unit-sector basis vector of `V = ⊞(X)` off into its own block. +function _split_unit(V::SumSpace, unit, ::Type{A}) where {A} + X = only(V) + c₀ = only(sectors(unit)) + rest = X ⊖ unit + T = scalartype(A) + TT = tensormaptype(spacetype(X), 1, 1, A) + Π₁ = TT(undef, unit ← X) + copyto!(block(Π₁, c₀), [one(T) zeros(T, 1, blockdim(X, c₀) - 1)]) + Π₂ = TT(undef, rest ← X) + for (c, b) in blocks(Π₂) + Id = Matrix{T}(LinearAlgebra.I, size(b, 1), size(b, 1)) + copyto!(b, c == c₀ ? [zeros(T, size(b, 1)) Id] : Id) + end + Π = BlockTensorMap{TT}(undef, (unit ⊞ rest) ← V) + Π[1, 1] = Π₁ + Π[2, 1] = Π₂ + return Π end diff --git a/test/algorithms/changebonds.jl b/test/algorithms/changebonds.jl index ab451e990..93b894b84 100644 --- a/test/algorithms/changebonds.jl +++ b/test/algorithms/changebonds.jl @@ -249,6 +249,37 @@ end @test expectation_value(ψ, H′) ≈ expectation_value(ψ, H) rtol = 1.0e-10 end +@testset "SvdCut on InfiniteMPOHamiltonian $(spacetype(pspace))" for (pspace, Dspace) in spacelist + Random.seed!(2468) + nn = rand(ComplexF64, pspace * pspace, pspace * pspace) + nn += nn' + L = 3 + lattice = PeriodicVector(fill(pspace, L)) + H0 = InfiniteMPOHamiltonian(lattice, (i, i + 1) => nn for i in 1:L) + + InfiniteMPOHamiltonian(lattice, (i, i + 2) => 0.5^i * nn for i in 1:L) + # adding the terms twice duplicates every interaction channel + H = H0 + 2 * H0 + ψ = InfiniteMPS(rand, ComplexF64, fill(pspace, L), fill(Dspace, L)) + + H′ = changebonds(H, SvdCut(; trunc = trunctol(; atol = 1.0e-12))) + @test all(i -> dim(right_virtualspace(H′[i])) ≤ dim(right_virtualspace(H0[i])), 1:L) + @test expectation_value(ψ, H′) ≈ expectation_value(ψ, H) rtol = 1.0e-10 + + # exponentially decaying interactions, written out as explicit terms + Z = rand(ComplexF64, pspace, pspace) + Z += Z' + H_exp = InfiniteMPOHamiltonian(PeriodicVector([pspace]), (1, 1 + r) => 0.5^r * (Z ⊗ Z) for r in 1:20) + H_exp′ = changebonds(H_exp, SvdCut(; trunc = trunctol(; rtol = 1.0e-5))) + @test dim(right_virtualspace(H_exp′[1])) < dim(right_virtualspace(H_exp[1])) + ψ1 = InfiniteMPS(rand, ComplexF64, [pspace], [Dspace]) + @test expectation_value(ψ1, H_exp′) ≈ expectation_value(ψ1, H_exp) rtol = 1.0e-4 + + # purely on-site terms have no channels to compress + H_onsite = InfiniteMPOHamiltonian(PeriodicVector([pspace]), (1,) => Z) + @test expectation_value(ψ1, changebonds(H_onsite, SvdCut(; trunc = truncrank(0)))) ≈ + expectation_value(ψ1, H_onsite) +end + # Regression: RandExpand used to be able to introduce a NaN entanglement entropy. @testset "RandExpand does not introduce NaN entropy" begin ψ = InfiniteMPS([ℂ^2], [ℂ^5]) diff --git a/test/hamiltonian/finite.jl b/test/hamiltonian/finite.jl index 14f7c9794..30abf990c 100644 --- a/test/hamiltonian/finite.jl +++ b/test/hamiltonian/finite.jl @@ -153,6 +153,42 @@ end end end +# Regression: the identity component of `C`/`B` was projected out without normalizing by `dim(P)`, +# and rank-deficient blocks left `Q` with directions that overlap with the identity. +@testset "Finite MPOHamiltonian canonical forms" begin + for V in (ℂ^2, U1Space(-1 => 1, 0 => 1, 1 => 1), Rep[SU₂](1 // 2 => 1)) + O₂ = randn(ComplexF64, V^2 ← V^2) + O₂ += O₂' + L = 4 + H = FiniteMPOHamiltonian(fill(V, L), (i, i + 1) => O₂ for i in 1:(L - 1)) + H_tm = convert(TensorMap, H) + + HL = copy(H) + for i in 1:(L - 1) + MPSKit.left_canonicalize!(HL, i) + end + @test convert(TensorMap, HL) ≈ H_tm + for W in HL[2:(L - 1)] + Vb = MPSKit.SparseBlockTensorMap(W)[1:(end - 1), 1, 1, 1:(end - 1)] + @plansor G[r; r′] := conj(Vb[l p; p′ r]) * Vb[l p; p′ r′] + G = TensorMap(G) / dim(physicalspace(W)) + @test G ≈ id(domain(G)) + end + + HR = copy(H) + for i in L:-1:2 + MPSKit.right_canonicalize!(HR, i) + end + @test convert(TensorMap, HR) ≈ H_tm + for W in HR[2:(L - 1)] + Vb = MPSKit.SparseBlockTensorMap(W)[2:end, 1, 1, 2:end] + @plansor G[l; l′] := Vb[l p; p′ r] * conj(Vb[l′ p; p′ r]) + G = TensorMap(G) / dim(physicalspace(W)) + @test G ≈ id(codomain(G)) + end + end +end + @testset "Finite MPOHamiltonian repeated indices" begin X = randn(ComplexF64, ℂ^2, ℂ^2) X += X' diff --git a/test/hamiltonian/infinite.jl b/test/hamiltonian/infinite.jl index f44ed42b4..1aa601901 100644 --- a/test/hamiltonian/infinite.jl +++ b/test/hamiltonian/infinite.jl @@ -68,3 +68,36 @@ end @test expectation_value(ψ1, Hreal_long_range + Hcomplex_local) ≈ expectation_value(ψ1, Hreal_long_range) + expectation_value(ψ1, Hcomplex_local) atol = 1.0e-10 end + +@testset "InfiniteMPOHamiltonian canonical forms $(sectortype(pspace))" for (pspace, Dspace) in zip(pspaces, vspaces) + O₂ = rand(ComplexF64, pspace^2, pspace^2) + O₂ += O₂' + O₃ = rand(ComplexF64, pspace^3, pspace^3) + O₃ += O₃' + lattice = PeriodicVector(fill(pspace, 2)) + H = InfiniteMPOHamiltonian(lattice, (i, i + 1) => O₂ for i in 1:2) + + InfiniteMPOHamiltonian(lattice, (i, i + 1, i + 2) => O₃ for i in 1:2) + + ψ = InfiniteMPS(rand, ComplexF64, fill(pspace, 2), fill(Dspace, 2)) + E = expectation_value(ψ, H) + + # the upper-left block `[1 C; 0 A]` has orthonormal columns + HL = MPSKit.left_canonicalize!(copy(H)) + @test expectation_value(ψ, HL) ≈ E + for W in parent(HL) + V = MPSKit.SparseBlockTensorMap(W)[1:(end - 1), 1, 1, 1:(end - 1)] + @plansor G[r; r′] := conj(V[l p; p′ r]) * V[l p; p′ r′] + G = TensorMap(G) / dim(physicalspace(W)) + @test G ≈ id(domain(G)) + end + + # the lower-right block `[A B; 0 1]` has orthonormal rows + HR = MPSKit.right_canonicalize!(copy(H)) + @test expectation_value(ψ, HR) ≈ E + for W in parent(HR) + V = MPSKit.SparseBlockTensorMap(W)[2:end, 1, 1, 2:end] + @plansor G[l; l′] := V[l p; p′ r] * conj(V[l′ p; p′ r]) + G = TensorMap(G) / dim(physicalspace(W)) + @test G ≈ id(codomain(G)) + end +end diff --git a/test/operators/finite_mpo_compression.jl b/test/operators/finite_mpo_compression.jl deleted file mode 100644 index 8b4cf3dec..000000000 --- a/test/operators/finite_mpo_compression.jl +++ /dev/null @@ -1,379 +0,0 @@ -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 diff --git a/test/operators/infinite_mpo_compression.jl b/test/operators/infinite_mpo_compression.jl deleted file mode 100644 index e9b7ce082..000000000 --- a/test/operators/infinite_mpo_compression.jl +++ /dev/null @@ -1,406 +0,0 @@ -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; truncbelow = 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; 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)) - 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; truncbelow = 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; 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)) - 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; truncbelow = 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; truncbelow = 10^-10) - H2c, _ = mpo_compression(H2; truncbelow = 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; truncbelow = 0) - mps2, = find_groundstate(mps2, H2, VUMPS(; maxiter = 100)) - E1 = expectation_value(mps2, H2) - - @assert E0 ≈ E1 -end From e3103fb81fe068b8b402b2b7a92985d4474d6496 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Mon, 12 Oct 2026 11:51:59 -0400 Subject: [PATCH 5/8] Merge infinite canonicalization into `left/right_canonicalize!` Drop the separate `uniform_leftorth!`/`uniform_rightorth!` helpers: the infinite `left_canonicalize!`/`right_canonicalize!` now run the iteration themselves and return the accumulated gauge alongside `H`. Fold the finite canonical-form tests into the existing canonicalization testset. Co-Authored-By: Claude Opus 5.5 --- src/algorithms/changebonds/svdcut.jl | 2 +- src/operators/ortho.jl | 56 +++++++++++++--------------- test/hamiltonian/finite.jl | 13 ++----- test/hamiltonian/infinite.jl | 4 +- 4 files changed, 32 insertions(+), 43 deletions(-) diff --git a/src/algorithms/changebonds/svdcut.jl b/src/algorithms/changebonds/svdcut.jl index ee8ef9fdb..a49211655 100644 --- a/src/algorithms/changebonds/svdcut.jl +++ b/src/algorithms/changebonds/svdcut.jl @@ -185,7 +185,7 @@ function changebonds!(H::InfiniteMPOHamiltonian, alg::SvdCut) # the interaction channels of the gauge between the right- and left-canonical forms right_canonicalize!(H) - _, Cs = uniform_leftorth!(H) + _, Cs = left_canonicalize!(H) Us = map(C -> first(svd_trunc(C; trunc = alg.trunc, alg = alg.alg_svd)), Cs) for i in 1:length(H) diff --git a/src/operators/ortho.jl b/src/operators/ortho.jl index c164e6024..28626fb15 100644 --- a/src/operators/ortho.jl +++ b/src/operators/ortho.jl @@ -29,11 +29,15 @@ end # ------------------------ """ - left_canonicalize!(H::InfiniteMPOHamiltonian; kwargs...) -> H + left_canonicalize!(H::InfiniteMPOHamiltonian; kwargs...) -> H, L -Bring `H` into left-canonical form, where the upper-left block `[1 C; 0 A]` of every site has -orthonormal columns, by iterating block-respecting QR decompositions until the gauge transformation -is unitary. Redundant interaction channels are removed along the way. +Bring `H` into left-canonical form in place, where the upper-left block `[1 C; 0 A]` of every site +has orthonormal columns, by iterating block-respecting QR decompositions until the gauge +transformation is unitary. Redundant interaction channels are removed along the way. + +Also returns the interaction-channel blocks `L` of the accumulated gauge transformations, such that +`diag(1, L[i - 1], 1) * H₀[i] = H[i] * diag(1, L[i], 1)` up to the identity components that are +gauged away. # Keyword arguments @@ -46,29 +50,7 @@ is unitary. Redundant interaction channels are removed along the way. * [Parker et al. Phys. Rev. B 102 (2020)](@cite parker2020) """ -function left_canonicalize!(H::InfiniteMPOHamiltonian; kwargs...) - uniform_leftorth!(H; kwargs...) - return H -end - -""" - right_canonicalize!(H::InfiniteMPOHamiltonian; kwargs...) -> H - -Bring `H` into right-canonical form, where the lower-right block `[A B; 0 1]` of every site has -orthonormal rows, by iterating block-respecting LQ decompositions until the gauge transformation -is unitary. Redundant interaction channels are removed along the way. - -See [`left_canonicalize!`](@ref) for the keyword arguments. -""" -function right_canonicalize!(H::InfiniteMPOHamiltonian; kwargs...) - uniform_rightorth!(H; kwargs...) - return H -end - -# Bring `H` into left-canonical form in place, and return the middle blocks `L` of the accumulated -# gauge transformations `diag(1, L[i - 1], 1) * H₀[i] = H[i] * diag(1, L[i], 1)` up to the gauging -# away of identity components. -function uniform_leftorth!( +function left_canonicalize!( H::InfiniteMPOHamiltonian; tol::Real = Defaults.tolgauge, maxiter::Int = Defaults.maxiter, verbosity::Int = VERBOSE_WARN ) @@ -99,13 +81,26 @@ function uniform_leftorth!( end end -# Bring `H` into right-canonical form in place. -function uniform_rightorth!( +""" + right_canonicalize!(H::InfiniteMPOHamiltonian; kwargs...) -> H, L + +Bring `H` into right-canonical form in place, where the lower-right block `[A B; 0 1]` of every site +has orthonormal rows, by iterating block-respecting LQ decompositions until the gauge transformation +is unitary. Redundant interaction channels are removed along the way. + +Also returns the interaction-channel blocks `L` of the accumulated gauge transformations, such that +`H₀[i] * diag(1, L[i], 1) = diag(1, L[i - 1], 1) * H[i]` up to the identity components that are +gauged away. + +See [`left_canonicalize!`](@ref) for the keyword arguments. +""" +function right_canonicalize!( H::InfiniteMPOHamiltonian; tol::Real = Defaults.tolgauge, maxiter::Int = Defaults.maxiter, verbosity::Int = VERBOSE_WARN ) alg_trunc = MatrixAlgebraKit.TruncatedAlgorithm(Defaults.alg_svd(), trunctol(; rtol = tol)) alg = MatrixAlgebraKit.RightOrthViaSVD(alg_trunc) + Ls = PeriodicVector(Vector{Any}(nothing, length(H))) return with_verbosity(; verbosity) do log = IterLog("MPO RC") @log_initialization loginit!(log, Inf) @@ -114,6 +109,7 @@ function uniform_rightorth!( for i in length(H):-1:1 R, t, H[i] = _right_orth_jordan(H[i]; alg) H[i - 1] = _absorb_right_gauge(H[i - 1], R, t, H[i]) + Ls[i - 1] = isnothing(Ls[i - 1]) || dim(R) == 0 ? R : Ls[i - 1] * R end ϵ = _gauge_distance(t, R) if ϵ <= tol @@ -125,7 +121,7 @@ function uniform_rightorth!( @log_iteration logiter!(log, iter, ϵ) end end - return H + return H, map(identity, Ls) end end diff --git a/test/hamiltonian/finite.jl b/test/hamiltonian/finite.jl index 30abf990c..7f4645001 100644 --- a/test/hamiltonian/finite.jl +++ b/test/hamiltonian/finite.jl @@ -113,8 +113,7 @@ using TensorKit: ℙ end end -@testset "Finite MPOHamiltonian canonicalization with decoupled sites" begin - # regression test for https://github.com/QuantumKitHub/MPSKit.jl/issues/522 +@testset "Finite MPOHamiltonian canonicalization" begin T = ComplexF64 for V in (ℂ^2, U1Space(-1 => 1, 0 => 1, 1 => 1), Rep[SU₂](1 // 2 => 1)) O₁ = randn(T, V, V) @@ -123,6 +122,7 @@ end O₂ += O₂' E = id(storagetype(O₁), domain(O₁)) + # regression test for https://github.com/QuantumKitHub/MPSKit.jl/issues/522 # bulk block of the middle site is empty on the left but not on the right H_right = FiniteMPOHamiltonian(fill(V, 3), (1,) => O₁, (2, 3) => O₂) H_right_tm = O₁ ⊗ E ⊗ E + E ⊗ O₂ @@ -150,15 +150,8 @@ end H5_trunc = changebonds(H5, SvdCut(; trunc = truncrank(16))) @test convert(TensorMap, H5_trunc) ≈ H5_tm - end -end -# Regression: the identity component of `C`/`B` was projected out without normalizing by `dim(P)`, -# and rank-deficient blocks left `Q` with directions that overlap with the identity. -@testset "Finite MPOHamiltonian canonical forms" begin - for V in (ℂ^2, U1Space(-1 => 1, 0 => 1, 1 => 1), Rep[SU₂](1 // 2 => 1)) - O₂ = randn(ComplexF64, V^2 ← V^2) - O₂ += O₂' + # canonical forms, also when `C`/`B` have identity components and the blocks are rank-deficient L = 4 H = FiniteMPOHamiltonian(fill(V, L), (i, i + 1) => O₂ for i in 1:(L - 1)) H_tm = convert(TensorMap, H) diff --git a/test/hamiltonian/infinite.jl b/test/hamiltonian/infinite.jl index 1aa601901..4dc2f4c51 100644 --- a/test/hamiltonian/infinite.jl +++ b/test/hamiltonian/infinite.jl @@ -82,7 +82,7 @@ end E = expectation_value(ψ, H) # the upper-left block `[1 C; 0 A]` has orthonormal columns - HL = MPSKit.left_canonicalize!(copy(H)) + HL, = MPSKit.left_canonicalize!(copy(H)) @test expectation_value(ψ, HL) ≈ E for W in parent(HL) V = MPSKit.SparseBlockTensorMap(W)[1:(end - 1), 1, 1, 1:(end - 1)] @@ -92,7 +92,7 @@ end end # the lower-right block `[A B; 0 1]` has orthonormal rows - HR = MPSKit.right_canonicalize!(copy(H)) + HR, = MPSKit.right_canonicalize!(copy(H)) @test expectation_value(ψ, HR) ≈ E for W in parent(HR) V = MPSKit.SparseBlockTensorMap(W)[2:end, 1, 1, 2:end] From 61895a88f3713906ea45f18d23462e6e5220324e Mon Sep 17 00:00:00 2001 From: lkdvos Date: Mon, 12 Oct 2026 12:40:42 -0400 Subject: [PATCH 6/8] Add PR link to changelog Co-Authored-By: Claude Opus 5.5 --- docs/src/changelog.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/changelog.md b/docs/src/changelog.md index e6523fd19..f39bd9275 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -26,7 +26,7 @@ When releasing a new version, move the "Unreleased" changes to a new version sec that matches exact exponentials on clusters of up to `N` sites. ([#252](https://github.com/QuantumKitHub/MPSKit.jl/pull/252)) - `changebonds(::InfiniteMPOHamiltonian, ::SvdCut)` compresses infinite Hamiltonians with the almost-Schmidt decomposition of [Parker et al.](https://arxiv.org/abs/1909.06341), and - `left_canonicalize!`/`right_canonicalize!` bring an `InfiniteMPOHamiltonian` into canonical form. + `left_canonicalize!`/`right_canonicalize!` bring an `InfiniteMPOHamiltonian` into canonical form. ([#539](https://github.com/QuantumKitHub/MPSKit.jl/pull/539)) ### Changed From 77e07c70c3dcd7c04674b80013a655842ae1b6b1 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Mon, 12 Oct 2026 14:02:05 -0400 Subject: [PATCH 7/8] Document finite and infinite `left/right_canonicalize!` together Co-Authored-By: Claude Opus 5.5 --- src/operators/ortho.jl | 86 +++++++++++++++++++++++++----------------- 1 file changed, 51 insertions(+), 35 deletions(-) diff --git a/src/operators/ortho.jl b/src/operators/ortho.jl index 28626fb15..898c0f37e 100644 --- a/src/operators/ortho.jl +++ b/src/operators/ortho.jl @@ -5,6 +5,38 @@ _dense(t) = t # Finite MPOHamiltonians # ---------------------- +""" + left_canonicalize!(H::FiniteMPOHamiltonian, i::Int; alg = Defaults.alg_orth()) -> H + left_canonicalize!(H::InfiniteMPOHamiltonian; tol, maxiter, verbosity) -> H, L + +Bring `H` into left-canonical form in place, where the upper-left block `[1 C; 0 A]` of a site has +orthonormal columns with respect to the trace inner product normalized by the physical dimension. +The identity components of `C` are gauged away, and the decompositions respect the block structure. + +For a `FiniteMPOHamiltonian`, this performs a single step: site `i` is brought into left-canonical +form with a QR decomposition, and the gauge transformation is absorbed into site `i + 1`. Any other +(possibly truncating) algorithm `alg` is applied to the interaction channels only, as used by +[`SvdCut`](@ref). + +For an `InfiniteMPOHamiltonian`, the decompositions are iterated over the unit cell until the gauge +transformation is unitary, removing redundant interaction channels along the way. This additionally +returns the interaction-channel blocks `L` of the accumulated gauge transformations, such that +`diag(1, L[i - 1], 1) * H₀[i] = H[i] * diag(1, L[i], 1)` up to the identity components that are +gauged away. + +# Keyword arguments + +- `alg = Defaults.alg_orth()`: (finite) algorithm for the decomposition of the interaction channels +- `tol = Defaults.tolgauge`: (infinite) convergence tolerance on the deviation of the gauge + transformation from a unitary one, which is also the relative tolerance for discarding redundant + channels +- `maxiter = Defaults.maxiter`: (infinite) maximal number of sweeps over the unit cell +- `verbosity = VERBOSE_WARN`: (infinite) setting for how much information is displayed + +# References + +* [Parker et al. Phys. Rev. B 102 (2020)](@cite parker2020) +""" function left_canonicalize!( H::FiniteMPOHamiltonian, i::Int; alg::MatrixAlgebraKit.AbstractAlgorithm = Defaults.alg_orth() @@ -15,6 +47,25 @@ function left_canonicalize!( return H end +""" + right_canonicalize!(H::FiniteMPOHamiltonian, i::Int; alg = Defaults.alg_orth()) -> H + right_canonicalize!(H::InfiniteMPOHamiltonian; tol, maxiter, verbosity) -> H, L + +Bring `H` into right-canonical form in place, where the lower-right block `[A B; 0 1]` of a site has +orthonormal rows with respect to the trace inner product normalized by the physical dimension. +The identity components of `B` are gauged away, and the decompositions respect the block structure. + +For a `FiniteMPOHamiltonian`, this performs a single step: site `i` is brought into right-canonical +form with an LQ decomposition, and the gauge transformation is absorbed into site `i - 1`. + +For an `InfiniteMPOHamiltonian`, the decompositions are iterated over the unit cell until the gauge +transformation is unitary, removing redundant interaction channels along the way. This additionally +returns the interaction-channel blocks `L` of the accumulated gauge transformations, such that +`H₀[i] * diag(1, L[i], 1) = diag(1, L[i - 1], 1) * H[i]` up to the identity components that are +gauged away. + +See [`left_canonicalize!`](@ref) for the keyword arguments. +""" function right_canonicalize!( H::FiniteMPOHamiltonian, i::Int; alg::MatrixAlgebraKit.AbstractAlgorithm = Defaults.alg_orth() @@ -28,28 +79,6 @@ end # Infinite MPOHamiltonians # ------------------------ -""" - left_canonicalize!(H::InfiniteMPOHamiltonian; kwargs...) -> H, L - -Bring `H` into left-canonical form in place, where the upper-left block `[1 C; 0 A]` of every site -has orthonormal columns, by iterating block-respecting QR decompositions until the gauge -transformation is unitary. Redundant interaction channels are removed along the way. - -Also returns the interaction-channel blocks `L` of the accumulated gauge transformations, such that -`diag(1, L[i - 1], 1) * H₀[i] = H[i] * diag(1, L[i], 1)` up to the identity components that are -gauged away. - -# Keyword arguments - -- `tol = Defaults.tolgauge`: convergence tolerance on the deviation of the gauge transformation from - a unitary one, which is also the relative tolerance for discarding redundant channels -- `maxiter = Defaults.maxiter`: maximal number of sweeps over the unit cell -- `verbosity = VERBOSE_WARN`: setting for how much information is displayed - -# References - -* [Parker et al. Phys. Rev. B 102 (2020)](@cite parker2020) -""" function left_canonicalize!( H::InfiniteMPOHamiltonian; tol::Real = Defaults.tolgauge, maxiter::Int = Defaults.maxiter, verbosity::Int = VERBOSE_WARN @@ -81,19 +110,6 @@ function left_canonicalize!( end end -""" - right_canonicalize!(H::InfiniteMPOHamiltonian; kwargs...) -> H, L - -Bring `H` into right-canonical form in place, where the lower-right block `[A B; 0 1]` of every site -has orthonormal rows, by iterating block-respecting LQ decompositions until the gauge transformation -is unitary. Redundant interaction channels are removed along the way. - -Also returns the interaction-channel blocks `L` of the accumulated gauge transformations, such that -`H₀[i] * diag(1, L[i], 1) = diag(1, L[i - 1], 1) * H[i]` up to the identity components that are -gauged away. - -See [`left_canonicalize!`](@ref) for the keyword arguments. -""" function right_canonicalize!( H::InfiniteMPOHamiltonian; tol::Real = Defaults.tolgauge, maxiter::Int = Defaults.maxiter, verbosity::Int = VERBOSE_WARN From 9d3a5ed28b1f5ad3da862b038f0443e9124c4298 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Mon, 12 Oct 2026 14:40:08 -0400 Subject: [PATCH 8/8] Whole-chain `left/right_canonicalize!` for finite MPOHamiltonians The finite methods now sweep over the whole chain and return the gauge transformations alongside the operator, matching the infinite versions, which simplifies the finite `SvdCut`. Also add changelog entries for the API change and the canonical-form fix. Co-Authored-By: Claude Opus 5.5 --- docs/src/changelog.md | 6 +++ src/algorithms/changebonds/svdcut.jl | 15 ++---- src/operators/ortho.jl | 70 ++++++++++++++-------------- test/hamiltonian/finite.jl | 27 ++++------- 4 files changed, 56 insertions(+), 62 deletions(-) diff --git a/docs/src/changelog.md b/docs/src/changelog.md index f39bd9275..64bb570b0 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -30,6 +30,10 @@ When releasing a new version, move the "Unreleased" changes to a new version sec ### Changed +- `left_canonicalize!`/`right_canonicalize!(::FiniteMPOHamiltonian)` now sweep over the whole chain + instead of a single site, and return the gauge transformations alongside the operator, matching + the infinite versions. ([#539](https://github.com/QuantumKitHub/MPSKit.jl/pull/539)) + ### Deprecated ### Removed @@ -38,6 +42,8 @@ When releasing a new version, move the "Unreleased" changes to a new version sec - `convert(TensorMap, H)` for single-site `FiniteMPOHamiltonian`s no longer contracts the site tensor with itself. ([#252](https://github.com/QuantumKitHub/MPSKit.jl/pull/252)) - `show` methods for finite MPS/MPOs no longer check unit space condition through `oneunit`, which is not defined for multifusion spaces. ([#538](https://github.com/QuantumKitHub/MPSKit.jl/pull/538)) +- The canonical forms of `MPOHamiltonian`s, and therefore `SvdCut` compression, are now correct when + interactions contain identity components or redundant channels. ([#539](https://github.com/QuantumKitHub/MPSKit.jl/pull/539)) ### Performance diff --git a/src/algorithms/changebonds/svdcut.jl b/src/algorithms/changebonds/svdcut.jl index a49211655..611ff87a6 100644 --- a/src/algorithms/changebonds/svdcut.jl +++ b/src/algorithms/changebonds/svdcut.jl @@ -153,22 +153,15 @@ end changebonds(mpo::FiniteMPOHamiltonian, alg::SvdCut) = changebonds!(copy(mpo), alg) function changebonds!(H::FiniteMPOHamiltonian, alg::SvdCut) # orthogonality center to the left - for i in length(H):-1:2 - H = right_canonicalize!(H, i) - end + right_canonicalize!(H) # a Jordan bond has a start and a finish level on top of its interaction channels channels_before = [dim(left_virtualspace(H, i)) - 2 for i in 2:length(H)] - # swipe right alg_trunc = MatrixAlgebraKit.TruncatedAlgorithm(alg.alg_svd, alg.trunc) - for i in 1:(length(H) - 1) - H = left_canonicalize!(H, i; alg = MatrixAlgebraKit.LeftOrthViaSVD(alg_trunc)) - end - # swipe left -- TODO: do we really need this double sweep? - for i in length(H):-1:2 - H = right_canonicalize!(H, i; alg = MatrixAlgebraKit.RightOrthViaSVD(alg_trunc)) - end + left_canonicalize!(H; alg = MatrixAlgebraKit.LeftOrthViaSVD(alg_trunc)) + # TODO: do we really need this double sweep? + right_canonicalize!(H; alg = MatrixAlgebraKit.RightOrthViaSVD(alg_trunc)) emptied = findall( i -> channels_before[i] > 0 && dim(left_virtualspace(H, i + 1)) - 2 == 0, diff --git a/src/operators/ortho.jl b/src/operators/ortho.jl index 898c0f37e..5a4779a51 100644 --- a/src/operators/ortho.jl +++ b/src/operators/ortho.jl @@ -6,23 +6,24 @@ _dense(t) = t # ---------------------- """ - left_canonicalize!(H::FiniteMPOHamiltonian, i::Int; alg = Defaults.alg_orth()) -> H + left_canonicalize!(H::FiniteMPOHamiltonian; alg = Defaults.alg_orth()) -> H, L left_canonicalize!(H::InfiniteMPOHamiltonian; tol, maxiter, verbosity) -> H, L -Bring `H` into left-canonical form in place, where the upper-left block `[1 C; 0 A]` of a site has -orthonormal columns with respect to the trace inner product normalized by the physical dimension. -The identity components of `C` are gauged away, and the decompositions respect the block structure. +Bring `H` into left-canonical form in place, where the upper-left block `[1 C; 0 A]` of every site +has orthonormal columns with respect to the trace inner product normalized by the physical +dimension. The identity components of `C` are gauged away, and the decompositions respect the block +structure. -For a `FiniteMPOHamiltonian`, this performs a single step: site `i` is brought into left-canonical -form with a QR decomposition, and the gauge transformation is absorbed into site `i + 1`. Any other -(possibly truncating) algorithm `alg` is applied to the interaction channels only, as used by +Also returns the interaction-channel blocks `L` of the gauge transformations on every bond, such +that `diag(1, L[i - 1], 1) * H₀[i] = H[i] * diag(1, L[i], 1)` up to the identity components that +are gauged away, with `L[0]` the identity for a `FiniteMPOHamiltonian`. + +For a `FiniteMPOHamiltonian`, this is a single sweep of QR decompositions from left to right. Any +other (possibly truncating) algorithm `alg` is applied to the interaction channels only, as used by [`SvdCut`](@ref). -For an `InfiniteMPOHamiltonian`, the decompositions are iterated over the unit cell until the gauge -transformation is unitary, removing redundant interaction channels along the way. This additionally -returns the interaction-channel blocks `L` of the accumulated gauge transformations, such that -`diag(1, L[i - 1], 1) * H₀[i] = H[i] * diag(1, L[i], 1)` up to the identity components that are -gauged away. +For an `InfiniteMPOHamiltonian`, the sweeps over the unit cell are repeated until the gauge +transformation is unitary, removing redundant interaction channels along the way. # Keyword arguments @@ -38,42 +39,43 @@ gauged away. * [Parker et al. Phys. Rev. B 102 (2020)](@cite parker2020) """ function left_canonicalize!( - H::FiniteMPOHamiltonian, i::Int; - alg::MatrixAlgebraKit.AbstractAlgorithm = Defaults.alg_orth() + H::FiniteMPOHamiltonian; alg::MatrixAlgebraKit.AbstractAlgorithm = Defaults.alg_orth() ) - 1 ≤ i < length(H) || throw(ArgumentError("Bounds error in canonicalize")) - H[i], t, R = _left_orth_jordan(H[i]; alg) - H[i + 1] = _absorb_left_gauge(H[i], t, R, H[i + 1]) - return H + Ls = map(1:(length(H) - 1)) do i + H[i], t, R = _left_orth_jordan(H[i]; alg) + H[i + 1] = _absorb_left_gauge(H[i], t, R, H[i + 1]) + return R + end + return H, Ls end """ - right_canonicalize!(H::FiniteMPOHamiltonian, i::Int; alg = Defaults.alg_orth()) -> H + right_canonicalize!(H::FiniteMPOHamiltonian; alg = Defaults.alg_orth()) -> H, L right_canonicalize!(H::InfiniteMPOHamiltonian; tol, maxiter, verbosity) -> H, L -Bring `H` into right-canonical form in place, where the lower-right block `[A B; 0 1]` of a site has -orthonormal rows with respect to the trace inner product normalized by the physical dimension. +Bring `H` into right-canonical form in place, where the lower-right block `[A B; 0 1]` of every site +has orthonormal rows with respect to the trace inner product normalized by the physical dimension. The identity components of `B` are gauged away, and the decompositions respect the block structure. -For a `FiniteMPOHamiltonian`, this performs a single step: site `i` is brought into right-canonical -form with an LQ decomposition, and the gauge transformation is absorbed into site `i - 1`. +Also returns the interaction-channel blocks `L` of the gauge transformations on every bond, such +that `H₀[i] * diag(1, L[i], 1) = diag(1, L[i - 1], 1) * H[i]` up to the identity components that +are gauged away, with `L[length(H)]` the identity for a `FiniteMPOHamiltonian`. -For an `InfiniteMPOHamiltonian`, the decompositions are iterated over the unit cell until the gauge -transformation is unitary, removing redundant interaction channels along the way. This additionally -returns the interaction-channel blocks `L` of the accumulated gauge transformations, such that -`H₀[i] * diag(1, L[i], 1) = diag(1, L[i - 1], 1) * H[i]` up to the identity components that are -gauged away. +For a `FiniteMPOHamiltonian`, this is a single sweep of LQ decompositions from right to left. For an +`InfiniteMPOHamiltonian`, the sweeps over the unit cell are repeated until the gauge transformation +is unitary, removing redundant interaction channels along the way. See [`left_canonicalize!`](@ref) for the keyword arguments. """ function right_canonicalize!( - H::FiniteMPOHamiltonian, i::Int; - alg::MatrixAlgebraKit.AbstractAlgorithm = Defaults.alg_orth() + H::FiniteMPOHamiltonian; alg::MatrixAlgebraKit.AbstractAlgorithm = Defaults.alg_orth() ) - 1 < i ≤ length(H) || throw(ArgumentError("Bounds error in canonicalize")) - R, t, H[i] = _right_orth_jordan(H[i]; alg) - H[i - 1] = _absorb_right_gauge(H[i - 1], R, t, H[i]) - return H + Ls = map(length(H):-1:2) do i + R, t, H[i] = _right_orth_jordan(H[i]; alg) + H[i - 1] = _absorb_right_gauge(H[i - 1], R, t, H[i]) + return R + end + return H, reverse!(Ls) end # Infinite MPOHamiltonians diff --git a/test/hamiltonian/finite.jl b/test/hamiltonian/finite.jl index 7f4645001..3fe4cd3ef 100644 --- a/test/hamiltonian/finite.jl +++ b/test/hamiltonian/finite.jl @@ -126,26 +126,21 @@ end # bulk block of the middle site is empty on the left but not on the right H_right = FiniteMPOHamiltonian(fill(V, 3), (1,) => O₁, (2, 3) => O₂) H_right_tm = O₁ ⊗ E ⊗ E + E ⊗ O₂ - H = MPSKit.right_canonicalize!(copy(H_right), 2) + H, = MPSKit.right_canonicalize!(copy(H_right)) @test convert(TensorMap, H) ≈ H_right_tm # bulk block of the middle site is empty on the right but not on the left H_left = FiniteMPOHamiltonian(fill(V, 3), (1, 2) => O₂, (3,) => O₁) H_left_tm = O₂ ⊗ E + E ⊗ E ⊗ O₁ - H = MPSKit.left_canonicalize!(copy(H_left), 2) + H, = MPSKit.left_canonicalize!(copy(H_left)) @test convert(TensorMap, H) ≈ H_left_tm L = 5 H5 = FiniteMPOHamiltonian(fill(V, L), (1, 2) => O₂, (3,) => O₁, (4, 5) => O₂) H5_tm = O₂ ⊗ E ⊗ E ⊗ E + E ⊗ E ⊗ O₁ ⊗ E ⊗ E + E ⊗ E ⊗ E ⊗ O₂ - H = copy(H5) - for i in L:-1:2 - MPSKit.right_canonicalize!(H, i) - end + H, = MPSKit.right_canonicalize!(copy(H5)) @test convert(TensorMap, H) ≈ H5_tm - for i in 1:(L - 1) - MPSKit.left_canonicalize!(H, i) - end + H, = MPSKit.left_canonicalize!(H) @test convert(TensorMap, H) ≈ H5_tm H5_trunc = changebonds(H5, SvdCut(; trunc = truncrank(16))) @@ -156,10 +151,9 @@ end H = FiniteMPOHamiltonian(fill(V, L), (i, i + 1) => O₂ for i in 1:(L - 1)) H_tm = convert(TensorMap, H) - HL = copy(H) - for i in 1:(L - 1) - MPSKit.left_canonicalize!(HL, i) - end + HL, Ls = MPSKit.left_canonicalize!(copy(H)) + @test length(Ls) == L - 1 + @test all(i -> codomain(Ls[i])[1] == right_virtualspace(HL[i].A), 1:(L - 1)) @test convert(TensorMap, HL) ≈ H_tm for W in HL[2:(L - 1)] Vb = MPSKit.SparseBlockTensorMap(W)[1:(end - 1), 1, 1, 1:(end - 1)] @@ -168,10 +162,9 @@ end @test G ≈ id(domain(G)) end - HR = copy(H) - for i in L:-1:2 - MPSKit.right_canonicalize!(HR, i) - end + HR, Rs = MPSKit.right_canonicalize!(copy(H)) + @test length(Rs) == L - 1 + @test all(i -> domain(Rs[i])[1] == left_virtualspace(HR[i + 1].A), 1:(L - 1)) @test convert(TensorMap, HR) ≈ H_tm for W in HR[2:(L - 1)] Vb = MPSKit.SparseBlockTensorMap(W)[2:end, 1, 1, 2:end]