diff --git a/src/algorithms/contractions/bondenv/als_solve.jl b/src/algorithms/contractions/bondenv/als_solve.jl index c244c7370..af083f618 100644 --- a/src/algorithms/contractions/bondenv/als_solve.jl +++ b/src/algorithms/contractions/bondenv/als_solve.jl @@ -19,11 +19,9 @@ Construct the tensor └-----------------------------------┘ ``` """ -function _tensor_Ra( - benv::BondEnv{T, S}, b::AbstractTensorMap{T, S, 2, 1} - ) where {T <: Number, S <: ElementarySpace} +function _tensor_Ra(benv::BondEnv, b::MPSTensor) return @autoopt @tensor Ra[DX1 Db1; DX0 Db0] := ( - benv[DX1 DY1; DX0 DY0] * b[Db0 DY0; db] * conj(b[Db1 DY1; db]) + benv[DX1 DY1; DX0 DY0] * b[Db0 db; DY0] * conj(b[Db1 db; DY1]) ) end @@ -44,10 +42,10 @@ Construct the tensor ``` """ function _tensor_Sa( - benv::BondEnv{T, S}, b::AbstractTensorMap{T, S, 2, 1}, a2b2::AbstractTensorMap{T, S, 2, 2} + benv::BondEnv, b::MPSTensor, a2b2::AbstractTensorMap{T, S, 2, 2} ) where {T <: Number, S <: ElementarySpace} - return @autoopt @tensor Sa[DX1 Db1; da] := ( - benv[DX1 DY1; DX0 DY0] * conj(b[Db1 DY1; db]) * a2b2[DX0 DY0; da db] + return @autoopt @tensor Sa[DX1 da; Db1] := ( + benv[DX1 DY1; DX0 DY0] * conj(b[Db1 db; DY1]) * a2b2[DX0 DY0; da db] ) end @@ -67,11 +65,9 @@ Construct the tensor └-----------------------------------┘ ``` """ -function _tensor_Rb( - benv::BondEnv{T, S}, a::AbstractTensorMap{T, S, 2, 1} - ) where {T <: Number, S <: ElementarySpace} +function _tensor_Rb(benv::BondEnv, a::MPSTensor) return @autoopt @tensor Rb[Da1 DY1; Da0 DY0] := ( - benv[DX1 DY1; DX0 DY0] * a[DX0 Da0; da] * conj(a[DX1 Da1; da]) + benv[DX1 DY1; DX0 DY0] * a[DX0 da; Da0] * conj(a[DX1 da; Da1]) ) end @@ -92,10 +88,10 @@ Construct the tensor ``` """ function _tensor_Sb( - benv::BondEnv{T, S}, a::AbstractTensorMap{T, S, 2, 1}, a2b2::AbstractTensorMap{T, S, 2, 2} + benv::BondEnv, a::MPSTensor, a2b2::AbstractTensorMap{T, S, 2, 2} ) where {T <: Number, S <: ElementarySpace} - return @autoopt @tensor Sb[Da1 DY1; db] := ( - benv[DX1 DY1; DX0 DY0] * conj(a[DX1 Da1; da]) * a2b2[DX0 DY0; da db] + return @autoopt @tensor Sb[Da1 db; DY1] := ( + benv[DX1 DY1; DX0 DY0] * conj(a[DX1 da; Da1]) * a2b2[DX0 DY0; da db] ) end @@ -116,30 +112,10 @@ Calculate the inner product ``` """ function inner_prod( - benv::BondEnv{T, S}, a1b1::AbstractTensorMap{T, S, 2, 2}, a2b2::AbstractTensorMap{T, S, 2, 2} + benv::BondEnv, a1b1::AbstractTensorMap{T, S, 2, 2}, a2b2::AbstractTensorMap{T, S, 2, 2} ) where {T <: Number, S <: ElementarySpace} return @autoopt @tensor benv[DX1 DY1; DX0 DY0] * - conj(a1b1[DX1 DY1; da db]) * - a2b2[DX0 DY0; da db] -end - -""" -$(SIGNATURES) - -Calculate the fidelity between two evolution steps -``` - |⟨a1,b1|a2,b2⟩|^2 - -------------------------- - ⟨a1,b1|a1,b1⟩⟨a2,b2|a2,b2⟩ -``` -""" -function fidelity( - benv::BondEnv{T, S}, a1b1::AbstractTensorMap{T, S, 2, 2}, a2b2::AbstractTensorMap{T, S, 2, 2} - ) where {T <: Number, S <: ElementarySpace} - b12 = inner_prod(benv, a1b1, a2b2) - b11 = inner_prod(benv, a1b1, a1b1) - b22 = inner_prod(benv, a2b2, a2b2) - return abs2(b12) / abs(b11 * b22) + conj(a1b1[DX1 DY1; da db]) * a2b2[DX0 DY0; da db] end """ @@ -153,14 +129,12 @@ Contract the axis between `a` and `b` tensors ``` """ function _combine_ab( - a::AbstractTensorMap{T, S, 2, 1}, b::AbstractTensorMap{T, S, 1, 2} + a::MPSTensor, b::AbstractTensorMap{T, S, 1, 2} ) where {T <: Number, S <: ElementarySpace} return @tensor ab[DX DY; da db] := a[DX da; D] * b[D; db DY] end -function _combine_ab( - a::AbstractTensorMap{T, S, 2, 1}, b::AbstractTensorMap{T, S, 2, 1} - ) where {T <: Number, S <: ElementarySpace} - return @tensor ab[DX DY; da db] := a[DX D; da] * b[D DY; db] +function _combine_ab(a::MPSTensor, b::MPSTensor) + return @tensor ab[DX DY; da db] := a[DX da; D] * b[D db; DY] end """ @@ -168,41 +142,40 @@ $(SIGNATURES) Calculate the cost function ``` - f(a,b) = ‖ |a1,b1⟩ - |a2,b2⟩ ‖^2 - = ⟨a1,b1|a1,b1⟩ - 2 Re⟨a1,b1|a2,b2⟩ + ⟨a2,b2|a2,b2⟩ + f(a,b) = ‖ |ψ1⟩ - |ψ2⟩ ‖^2 + = ⟨ψ1|benv|ψ1⟩ - 2 Re⟨ψ1|benv|ψ2⟩ + ⟨ψ2|benv|ψ2⟩ +``` +and the fidelity +``` + |⟨ψ1|benv|ψ2⟩|² + ------------------------ + ⟨ψ1|benv|ψ1⟩⟨ψ2|benv|ψ2⟩ ``` """ -function cost_function_als( - benv::BondEnv{T, S}, a1b1::AbstractTensorMap{T, S, 2, 2}, a2b2::AbstractTensorMap{T, S, 2, 2} - ) where {T <: Number, S <: ElementarySpace} - t1 = inner_prod(benv, a1b1, a1b1) - t2 = inner_prod(benv, a2b2, a2b2) - t3 = inner_prod(benv, a1b1, a2b2) - return real(t1) + real(t2) - 2 * real(t3) +function cost_function_als(benv, ψ1, ψ2) + b12 = inner_prod(benv, ψ1, ψ2) + b11 = inner_prod(benv, ψ1, ψ1) + b22 = inner_prod(benv, ψ2, ψ2) + cost = real(b11) + real(b22) - 2 * real(b12) + fid = abs2(b12) / abs(b11 * b22) + return cost, fid end """ $(SIGNATURES) -Solve the equations `Rx x = Sx` (x = a, b) with initial guess `x0` -``` - ┌---------------------------┐ - | ┌----┐ | - └---| |--- 1 -- x -- 2 --┘ - | | ↓ - | Rx | -3 - | | - ┌---| |--- -1 -2 --┐ - | └----┘ | - └---------------------------┘ -``` -""" -function _solve_ab( - Rx::AbstractTensorMap{T, S, 2, 2}, - Sx::AbstractTensorMap{T, S, 2, 1}, - x0::AbstractTensorMap{T, S, 2, 1}, - ) where {T <: Number, S <: ElementarySpace} - f(x) = (@tensor Sx2[-1 -2; -3] := Rx[-1 -2; 1 2] * x[1 2; -3]) - x1, info = linsolve(f, Sx, x0, 0, 1) +Solve the equations `Rx x = Sx` with initial guess `x0`. +""" +function _solve_als( + Rx::AbstractTensorMap{T, S, N, N}, + Sx::GenericMPSTensor{S, N}, + x0::GenericMPSTensor{S, N}; kwargs... + ) where {T, S, N} + @assert N >= 2 + pR = (codomainind(Rx), domainind(Rx)) + pX = ((1, (3:(N + 1))...), (2,)) + pRX = ((1, N + 1, (2:(N - 1))...), (N,)) + f(x) = tensorcontract(Rx, pR, false, x, pX, false, pRX) + x1, info = linsolve(f, Sx, x0, 0, 1; kwargs...) return x1, info end diff --git a/src/algorithms/contractions/bondenv/benv_ctm.jl b/src/algorithms/contractions/bondenv/benv_ctm.jl index 9b97c1576..45262cc67 100644 --- a/src/algorithms/contractions/bondenv/benv_ctm.jl +++ b/src/algorithms/contractions/bondenv/benv_ctm.jl @@ -1,14 +1,16 @@ """ Construct the environment (norm) tensor ``` - C1---T1---------T1---C2 r-1 - | ‖ ‖ | - T4===XX== ==YY===T2 r - | ‖ ‖ | - C4---T3---------T3---C3 r+1 - c-1 c c+1 c+2 + -1 C1---E1---------E1---C2 + | ‖ ‖ | + 0 E4===XX== ==YY===E2 + | ‖ ‖ | + 1 C4---E3---------E3---C3 + -1 0 1 2 ``` -where `XX = X' X` and `YY = Y' Y` (stacked together). +with `X` at position `[row, col]`. +`X, Y` are unitary tensors produced when finding the reduced site tensors +with `_qr_bond`; and `XX = X' X` and `YY = Y' Y` (stacked together). Axis order: `[DX1 DY1; DX0 DY0]`, as in ``` @@ -21,7 +23,9 @@ Axis order: `[DX1 DY1; DX0 DY0]`, as in └---------------------┘ ``` """ -function bondenv_ctm(row::Int, col::Int, X::PEPSOrth, Y::PEPSOrth, env::CTMRGEnv) +function bondenv_ctm( + row::Int, col::Int, X::T, Y::T, env::CTMRGEnv + ) where {T} Nr, Nc = size(env.corners)[[2, 3]] cm1 = _prev(col, Nc) cp1 = _next(col, Nc) @@ -32,29 +36,29 @@ function bondenv_ctm(row::Int, col::Int, X::PEPSOrth, Y::PEPSOrth, env::CTMRGEnv c2 = env.corners[2, rm1, cp2] c3 = env.corners[3, rp1, cp2] c4 = env.corners[4, rp1, cm1] - t1X, t1Y = env.edges[1, rm1, col], env.edges[1, rm1, cp1] - t2 = env.edges[2, row, cp2] - t3X, t3Y = env.edges[3, rp1, col], env.edges[3, rp1, cp1] - t4 = env.edges[4, row, cm1] + e1X, e1Y = env.edges[1, rm1, col], env.edges[1, rm1, cp1] + e2 = env.edges[2, row, cp2] + e3X, e3Y = env.edges[3, rp1, col], env.edges[3, rp1, cp1] + e4 = env.edges[4, row, cm1] #= index labels - C1--χ4--T1X---------χ6---------T1Y--χ8---C2 r-1 + C1--χ4--E1X---------χ6---------E1Y--χ8---C2 r-1 | ‖ ‖ | χ2 DNX DNY χ10 | ‖ ‖ | - T4==DWX==XX===DX== ==DY===YY==DEY==T2 r + E4==DWX==XX===DX== ==DY===YY==DEY==E2 r | ‖ ‖ | χ1 DSX DSY χ9 | ‖ ‖ | - C4--χ3--T3X---------χ5---------T3Y--χ7---C3 r+1 + C4--χ3--E3X---------χ5---------E3Y--χ7---C3 r+1 c-1 c c+1 c+2 =# + X, Y = _prepare_site_tensor(X), _prepare_site_tensor(Y) @autoopt @tensor benv[DX1 DY1; DX0 DY0] := - c4[χ3 χ1] * t4[χ1 DWX0 DWX1 χ2] * c1[χ2 χ4] * t3X[χ5 DSX0 DSX1 χ3] * - X[DNX0 DX0 DSX0 DWX0] * conj(X[DNX1 DX1 DSX1 DWX1]) * t1X[χ4 DNX0 DNX1 χ6] * - c3[χ9 χ7] * t2[χ10 DEY0 DEY1 χ9] * c2[χ8 χ10] * t3Y[χ7 DSY0 DSY1 χ5] * - Y[DNY0 DEY0 DSY0 DY0] * conj(Y[DNY1 DEY1 DSY1 DY1]) * t1Y[χ6 DNY0 DNY1 χ8] - + c4[χ3 χ1] * e4[χ1 DWX0 DWX1 χ2] * c1[χ2 χ4] * e3X[χ5 DSX0 DSX1 χ3] * + twistdual(X, 1)[pX DNX0 DX0 DSX0 DWX0] * conj(X[pX DNX1 DX1 DSX1 DWX1]) * e1X[χ4 DNX0 DNX1 χ6] * + c3[χ9 χ7] * e2[χ10 DEY0 DEY1 χ9] * c2[χ8 χ10] * e3Y[χ7 DSY0 DSY1 χ5] * + twistdual(Y, 1)[pY DNY0 DEY0 DSY0 DY0] * conj(Y[pY DNY1 DEY1 DSY1 DY1]) * e1Y[χ6 DNY0 DNY1 χ8] normalize!(benv, Inf) return benv end diff --git a/src/algorithms/contractions/bondenv/benv_tools.jl b/src/algorithms/contractions/bondenv/benv_tools.jl index 91ad5958e..38400b919 100644 --- a/src/algorithms/contractions/bondenv/benv_tools.jl +++ b/src/algorithms/contractions/bondenv/benv_tools.jl @@ -1,3 +1,13 @@ const BondEnv{T, S} = AbstractTensorMap{T, S, 2, 2} where {T <: Number, S <: ElementarySpace} +const BondEnv3site{T, S} = AbstractTensorMap{T, S, 4, 4} where {T <: Number, S <: ElementarySpace} +const Hair{T, S} = AbstractTensor{T, S, 2} where {T <: Number, S <: ElementarySpace} +# Orthogonal tensors obtained PEPSTensor/PEPOTensor +# with one physical leg factored out by `_qr_bond` const PEPSOrth{T, S} = AbstractTensor{T, S, 4} where {T <: Number, S <: ElementarySpace} const PEPOOrth{T, S} = AbstractTensor{T, S, 5} where {T <: Number, S <: ElementarySpace} + +"Convert tensor `t` connected by the bond to be truncated to a `PEPSTensor`." +_prepare_site_tensor(t::PEPSTensor) = t +_prepare_site_tensor(t::PEPOTensor) = first(fuse_physicalspaces(t)) +_prepare_site_tensor(t::PEPSOrth) = permute(insertleftunit(t, 1), ((1,), (2, 3, 4, 5))) +_prepare_site_tensor(t::PEPOOrth) = permute(t, ((1,), (2, 3, 4, 5))) diff --git a/src/algorithms/contractions/bondenv/gaugefix.jl b/src/algorithms/contractions/bondenv/gaugefix.jl index bfd0f7f01..0572d4852 100644 --- a/src/algorithms/contractions/bondenv/gaugefix.jl +++ b/src/algorithms/contractions/bondenv/gaugefix.jl @@ -1,5 +1,5 @@ """ -Replace bond environment `benv` by its positive approximant `Z† Z` +Find positive approximant `Z† Z` of a norm tensor `benv` (returns the "half environment" `Z`) ``` ┌-----------------┐ ┌---------------┐ @@ -11,9 +11,9 @@ Replace bond environment `benv` by its positive approximant `Z† Z` └-----------------┘ └---------------┘ ``` """ -function positive_approx(benv::BondEnv) +function positive_approx(benv::AbstractTensorMap{T, S, N, N}) where {T, S, N} # eigen-decomposition: benv = U D U' - D, U = eigh_full((benv + benv') / 2) + D, U = eigh_full!(project_hermitian(benv)) # determine if `env` is (mostly) positive or negative sgn = sign(sum(D.data)) # When optimizing the truncation of a bond, diff --git a/src/algorithms/time_evolution/apply_mpo.jl b/src/algorithms/time_evolution/apply_mpo.jl index b60818530..b7ae6f42f 100644 --- a/src/algorithms/time_evolution/apply_mpo.jl +++ b/src/algorithms/time_evolution/apply_mpo.jl @@ -112,7 +112,7 @@ Then the fidelity is just ``` =# """ -Perform QR decomposition through a PEPS tensor +Perform QR decomposition through a `GenericMPSTensor` ``` ╱ ╱ -←-R0-←-M-←- => ---Q-←-R1-←- @@ -120,19 +120,22 @@ Perform QR decomposition through a PEPS tensor ``` """ function qr_through( - R0::MPSBondTensor, M::GenericMPSTensor{S, 4}; normalize::Bool = true - ) where {S <: ElementarySpace} + R0::MPSBondTensor, M::GenericMPSTensor{S, N}; normalize::Bool = true + ) where {S, N} @assert !isdual(codomain(R0, 1)) @assert !isdual(domain(M, 1)) && !isdual(codomain(M, 1)) - @tensor A[-1 -2 -3 -4; -5] := R0[-1; 1] * M[1 -2 -3 -4; -5] + pR = (codomainind(R0), domainind(R0)) + pM = ((1,), Tuple(2:(N + 1))) + pRM = (codomainind(M), domainind(M)) + A = tensorcontract(R0, pR, false, M, pM, false, pRM) _, r = left_orth!(A; positive = true) normalize && normalize!(r, Inf) return r end # for `M` at the left end of the MPS function qr_through( - ::Nothing, M::GenericMPSTensor{S, 4}; normalize::Bool = true - ) where {S <: ElementarySpace} + ::Nothing, M::GenericMPSTensor{S, N}; normalize::Bool = true + ) where {S, N} @assert !isdual(domain(M, 1)) _, r = left_orth(M; positive = true) normalize && normalize!(r, Inf) @@ -140,7 +143,7 @@ function qr_through( end """ -Perform LQ decomposition through a tensor +Perform LQ decomposition through a `GenericMPSTensor` ``` ╱ ╱ -←-L0-←-Q-←- <= -←-M-←-L1-←- @@ -148,21 +151,24 @@ Perform LQ decomposition through a tensor ``` """ function lq_through( - M::GenericMPSTensor{S, 4}, L1::MPSBondTensor; normalize::Bool = true - ) where {S <: ElementarySpace} + M::GenericMPSTensor{S, N}, L1::MPSBondTensor; normalize::Bool = true + ) where {S, N} @assert !isdual(domain(L1, 1)) @assert !isdual(codomain(M, 1)) && !isdual(domain(M, 1)) - @tensor A[-1; -2 -3 -4 -5] := M[-1 -2 -3 -4; 1] * L1[1; -5] + pM = (codomainind(M), domainind(M)) + pL = (codomainind(L1), domainind(L1)) + pML = ((1,), Tuple(2:(N + 1))) + A = tensorcontract(M, pM, false, L1, pL, false, pML) l, _ = right_orth!(A; positive = true) normalize && normalize!(l, Inf) return l end # for `M` at the right end of the MPS function lq_through( - M::GenericMPSTensor{S, 4}, ::Nothing; normalize::Bool = true - ) where {S <: ElementarySpace} + M::GenericMPSTensor{S, N}, ::Nothing; normalize::Bool = true + ) where {S, N} @assert !isdual(codomain(M, 1)) - A = permute(M, ((1,), (2, 3, 4, 5))) + A = permute(M, ((1,), Tuple(2:(N + 1))); copy = true) l, _ = right_orth!(A; positive = true) normalize && normalize!(l, Inf) return l @@ -171,7 +177,7 @@ end """ Given a cluster `Ms`, find all `R`, `L` matrices on each internal bond """ -function _get_allRLs(Ms::Vector{T}) where {T <: GenericMPSTensor{<:ElementarySpace, 4}} +function _get_allRLs(Ms::Vector{T}) where {T <: GenericMPSTensor} # M1 -- (R1,L1) -- M2 -- (R2,L2) -- M3 N = length(Ms) # get the first R and the last L @@ -214,11 +220,14 @@ function _proj_from_RL( end """ -Given a cluster `Ms` and the pre-calculated `R`, `L` bond matrices, -find all projectors `Pa`, `Pb` and Schmidt weights `wts` on internal bonds. +Given a cluster `Ms`, find all projectors `Pa`, `Pb` +and Schmidt weights `wts` on internal bonds. """ -function _get_allprojs(Ms, Rs, Ls, truncs::Vector{E}) where {E <: TruncationStrategy} +function _get_allprojs( + Ms::Vector{T}, truncs::Vector{E} + ) where {T <: GenericMPSTensor, E <: TruncationStrategy} N = length(Ms) + Rs, Ls = _get_allRLs(Ms) @assert length(truncs) == N - 1 projs_errs = map(1:(N - 1)) do i trunc = if isa(truncs[i], FixedSpaceTruncation) @@ -258,14 +267,16 @@ Find projectors to truncate internal bonds of the cluster `Ms`. """ function _cluster_truncate!( Ms::Vector{T}, truncs::Vector{E} - ) where {T <: GenericMPSTensor{<:ElementarySpace, 4}, E <: TruncationStrategy} - Rs, Ls = _get_allRLs(Ms) - Pas, Pbs, wts, ϵs = _get_allprojs(Ms, Rs, Ls, truncs) + ) where {T <: GenericMPSTensor, E <: TruncationStrategy} + Pas, Pbs, wts, ϵs = _get_allprojs(Ms, truncs) # apply projectors # M1 -- (Pa1,wt1,Pb1) -- M2 -- (Pa2,wt2,Pb2) -- M3 for (i, (Pa, Pb)) in enumerate(zip(Pas, Pbs)) - @tensor (Ms[i])[-1 -2 -3 -4; -5] := (Ms[i])[-1 -2 -3 -4; 1] * Pa[1; -5] - @tensor (Ms[i + 1])[-1 -2 -3 -4; -5] := Pb[-1; 1] * (Ms[i + 1])[1 -2 -3 -4; -5] + Ms[i] = Ms[i] * twistdual(Pa, 1) + pP = ((1,), (2,)) + pM = ((1,), Tuple(2:numind(Ms[i + 1]))) + pPM = (codomainind(Ms[i + 1]), domainind(Ms[i + 1])) + Ms[i + 1] = tensorcontract(Pb, pP, false, Ms[i + 1], pM, false, pPM) end return wts, ϵs, Pas, Pbs end diff --git a/src/algorithms/truncation/bond_truncation.jl b/src/algorithms/truncation/bond_truncation.jl index d8feb8cb3..57fd68cc7 100644 --- a/src/algorithms/truncation/bond_truncation.jl +++ b/src/algorithms/truncation/bond_truncation.jl @@ -74,17 +74,11 @@ function bond_truncate( perm_ab = ((1, 3), (4, 2)) a, s0, b = svd_trunc(permute(a2b2, perm_ab); trunc = alg.trunc) a, b = absorb_s(a, s0, b) - #= temporarily reorder axes of a and b to - 1 -a/b- 2 - ↓ - 3 - =# - perm = ((1, 3), (2,)) - a, b = permute(a, perm), permute(b, perm) + # put b in MPS axis order + b = permute(b, ((1, 2), (3,))) ab = _combine_ab(a, b) # cost function will be normalized by initial value - cost00 = cost_function_als(benv, ab, a2b2) - fid = fidelity(benv, ab, a2b2) + cost00, fid = cost_function_als(benv, ab, a2b2) cost0, fid0, Δcost, Δfid, Δs = cost00, fid, NaN, NaN, NaN verbose && @info "ALS init" * _als_message(0, cost0, fid, Δcost, Δfid, Δs, 0.0) for iter in 1:(alg.maxiter) @@ -99,15 +93,14 @@ function bond_truncate( =# Ra = _tensor_Ra(benv, b) Sa = _tensor_Sa(benv, b, a2b2) - a, info_a = _solve_ab(Ra, Sa, a) + a, info_a = _solve_als(Ra, Sa, a) # Fixing `a`, solve for `b` from `Rb b = Sb` Rb = _tensor_Rb(benv, a) Sb = _tensor_Sb(benv, a, a2b2) - b, info_b = _solve_ab(Rb, Sb, b) + b, info_b = _solve_als(Rb, Sb, b) @debug "Bond truncation info" info_a info_b ab = _combine_ab(a, b) - cost = cost_function_als(benv, ab, a2b2) - fid = fidelity(benv, ab, a2b2) + cost, fid = cost_function_als(benv, ab, a2b2) # TODO: replace with truncated svdvals (without calculating u, vh) _, s, _ = svd_trunc!(permute(ab, perm_ab); trunc = alg.trunc) # fidelity, cost and normalized bond-s change diff --git a/src/algorithms/truncation/fullenv_truncation.jl b/src/algorithms/truncation/fullenv_truncation.jl index ea339c8d8..fd9c685bd 100644 --- a/src/algorithms/truncation/fullenv_truncation.jl +++ b/src/algorithms/truncation/fullenv_truncation.jl @@ -49,24 +49,7 @@ between two states specified by the bond matrices `b1`, `b2` function inner_prod( benv::BondEnv{T, S}, b1::AbstractTensorMap{T, S, 1, 1}, b2::AbstractTensorMap{T, S, 1, 1} ) where {T <: Number, S <: ElementarySpace} - val = @tensor conj(b1[1; 2]) * benv[1 2; 3 4] * b2[3; 4] - return val -end - -""" -$(SIGNATURES) - -Given the bond environment `benv`, calculate the fidelity -between two states specified by the bond matrices `b1`, `b2` -``` - F(b1, b2) = (⟨b1|b2⟩ ⟨b2|b1⟩) / (⟨b1|b1⟩ ⟨b2|b2⟩) -``` -""" -function fidelity( - benv::BondEnv{T, S}, b1::AbstractTensorMap{T, S, 1, 1}, b2::AbstractTensorMap{T, S, 1, 1} - ) where {T <: Number, S <: ElementarySpace} - return abs2(inner_prod(benv, b1, b2)) / - real(inner_prod(benv, b1, b1) * inner_prod(benv, b2, b2)) + return @tensor conj(b1[1; 2]) * benv[1 2; 3 4] * b2[3; 4] end """ @@ -252,7 +235,7 @@ function fullenv_truncate( l, info_l = linsolve(Base.Fix1(*, B), p, l, 0, 1) @debug "Bond truncation info" info_l info_r @tensor b1[-1; -2] = l[-1 1] * vh[1; -2] - fid = fidelity(benv, b0, b1) + _, fid = cost_function_als(benv, b0, b1) u, s, vh = svd_trunc!(b1; trunc = alg.trunc) # determine convergence s_nrm = norm(s0, Inf) diff --git a/test/bondenv/benv_ctm.jl b/test/bondenv/benv_ctm.jl index c6406ba43..7676c5139 100644 --- a/test/bondenv/benv_ctm.jl +++ b/test/bondenv/benv_ctm.jl @@ -7,8 +7,12 @@ using Random Random.seed!(100) Nr, Nc = 2, 2 +Envspace = Vect[FermionParity ⊠ U1Irrep]( + (0, 0) => 4, (1, 1 // 2) => 1, (1, -1 // 2) => 1, (0, 1) => 1, (0, -1) => 1 +) +ctm_alg = SequentialCTMRG(; tol = 1.0e-10, verbosity = 2, trunc = truncerror(; atol = 1.0e-10) & truncrank(8)) # create Hubbard iPEPS using simple update -function get_hubbard_state(t::Float64 = 1.0, U::Float64 = 8.0) +function get_hubbard_peps(t::Float64 = 1.0, U::Float64 = 8.0) H = hubbard_model(ComplexF64, Trivial, U1Irrep, InfiniteSquare(Nr, Nc); t, U, mu = U / 2) Vphy = Vect[FermionParity ⊠ U1Irrep]((0, 0) => 2, (1, 1 // 2) => 1, (1, -1 // 2) => 1) peps = InfinitePEPS(rand, ComplexF64, Vphy, Vphy; unitcell = (Nr, Nc)) @@ -20,30 +24,45 @@ function get_hubbard_state(t::Float64 = 1.0, U::Float64 = 8.0) return peps end -peps = get_hubbard_state() -# calculate CTMRG environment -Envspace = Vect[FermionParity ⊠ U1Irrep]( - (0, 0) => 4, (1, 1 // 2) => 1, (1, -1 // 2) => 1, (0, 1) => 1, (0, -1) => 1 -) -ctm_alg = SequentialCTMRG(; tol = 1.0e-10, verbosity = 2, trunc = truncerror(; atol = 1.0e-10) & truncrank(8)) -env, = leading_boundary(CTMRGEnv(rand, ComplexF64, peps, Envspace), peps, ctm_alg) -for row in 1:Nr, col in 1:Nc - cp1 = PEPSKit._next(col, Nc) - A, B = peps.A[row, col], peps.A[row, cp1] - X, a, b, Y = PEPSKit._qr_bond(A, B) - benv = PEPSKit.bondenv_ctm(row, col, X, Y, env) - @assert [isdual(space(benv, ax)) for ax in 1:numind(benv)] == [0, 0, 1, 1] - Z = PEPSKit.positive_approx(benv) - # verify that gauge fixing can greatly reduce - # condition number for physical state bond envs - cond1 = cond(Z' * Z) - Z2, a2, b2, (Linv, Rinv) = PEPSKit.fixgauge_benv(Z, a, b) - benv2 = Z2' * Z2 - cond2 = cond(benv2) - @test 1 <= cond2 < cond1 - @info "benv cond number: (gauge-fixed) $(cond2) ≤ $(cond1) (initial)" - # verify gauge fixing is done correctly - @tensor half[:] := Z[-1; 1 3] * a[1; -2 2] * b[2 -3; 3] - @tensor half2[:] := Z2[-1; 1 3] * a2[1; -2 2] * b2[2 -3; 3] - @test half ≈ half2 +function get_hubbard_pepo(t::Float64 = 1.0, U::Float64 = 8.0) + H = hubbard_model(ComplexF64, Trivial, U1Irrep, InfiniteSquare(Nr, Nc); t, U, mu = U / 2) + pepo = PEPSKit.infinite_temperature_density_matrix(H) + wts = SUWeight(pepo) + alg = SimpleUpdate(; + trunc = truncerror(; atol = 1.0e-10) & truncrank(4), bipartite = false + ) + pepo, = time_evolve(pepo, H, 2.0e-3, 500, alg, wts; check_interval = 100) + normalize!.(pepo.A, Inf) + return pepo +end + +function test_benv_ctm(state::Union{InfinitePEPS, InfinitePEPO}) + network = isa(state, InfinitePEPS) ? state : InfinitePEPS(state) + env, = leading_boundary(CTMRGEnv(rand, ComplexF64, network, Envspace), network, ctm_alg) + for row in 1:Nr, col in 1:Nc + cp1 = PEPSKit._next(col, Nc) + A, B = state.A[row, col], state.A[row, cp1] + X, a, b, Y = PEPSKit._qr_bond(A, B) + benv = PEPSKit.bondenv_ctm(row, col, X, Y, env) + Z = PEPSKit.positive_approx(benv) + # verify that gauge fixing can greatly reduce + # condition number for physical state bond envs + cond1 = cond(Z' * Z) + Z2, a2, b2, (Linv, Rinv) = PEPSKit.fixgauge_benv(Z, a, b) + benv2 = Z2' * Z2 + cond2 = cond(benv2) + @test 1 <= cond2 < cond1 + @info "benv cond number: (gauge-fixed) $(cond2) ≤ $(cond1) (initial)" + # verify gauge fixing is done correctly + @tensor half[:] := Z[-1; 1 3] * a[1; -2 2] * b[2 -3; 3] + @tensor half2[:] := Z2[-1; 1 3] * a2[1; -2 2] * b2[2 -3; 3] + @test half ≈ half2 + end + return +end + +peps = get_hubbard_peps() +pepo = get_hubbard_pepo() +for state in (peps, pepo) + test_benv_ctm(state) end diff --git a/test/bondenv/bond_truncate.jl b/test/bondenv/bond_truncate.jl index 58b0a4158..25f56c200 100644 --- a/test/bondenv/bond_truncate.jl +++ b/test/bondenv/bond_truncate.jl @@ -5,6 +5,7 @@ using TensorKit using PEPSKit using LinearAlgebra using KrylovKit +using PEPSKit: cost_function_als Random.seed!(0) maxiter = 600 @@ -19,6 +20,7 @@ for Vbondl in (Vint, Vint'), Vbondr in (Vint, Vint') # random positive-definite environment Z = randn(Float64, Vext ← Vbond) benv = Z' * Z + normalize!(benv, Inf) # untruncated bond tensor a2b2 = randn(Float64, Vbondl ⊗ Vbondr ← Vphy' ⊗ Vphy') a2, s, b2 = svd_compact(permute(a2b2, perm_ab)) @@ -26,7 +28,7 @@ for Vbondl in (Vint, Vint'), Vbondr in (Vint, Vint') # bond tensor (truncated SVD initialization) a0, s, b0 = svd_trunc(permute(a2b2, perm_ab); trunc = trunc) a0, b0 = PEPSKit.absorb_s(a0, s, b0) - fid0 = PEPSKit.fidelity(benv, PEPSKit._combine_ab(a0, b0), a2b2) + fid0 = cost_function_als(benv, PEPSKit._combine_ab(a0, b0), a2b2)[2] @info "Fidelity of simple SVD truncation = $fid0.\n" ss = Dict{String, DiagonalTensorMap}() for (label, alg) in ( @@ -36,7 +38,7 @@ for Vbondl in (Vint, Vint'), Vbondr in (Vint, Vint') a1, ss[label], b1, info = PEPSKit.bond_truncate(a2, b2, benv, alg) @info "$label improved fidelity = $(info.fid)." # display(ss[label]) - @test info.fid ≈ PEPSKit.fidelity(benv, PEPSKit._combine_ab(a1, b1), a2b2) + @test info.fid ≈ cost_function_als(benv, PEPSKit._combine_ab(a1, b1), a2b2)[2] @test info.fid > fid0 end @test isapprox(ss["ALS"], ss["FET"], atol = 1.0e-3)