From d77ab38b3a70d3819942d9f4523593c925736a22 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Mon, 26 May 2025 12:17:24 +0800 Subject: [PATCH 01/11] Add full update bond environment and gauge fixing --- src/PEPSKit.jl | 2 + .../contractions/bondenv/benv_ctm.jl | 77 ++++++++++++ .../contractions/bondenv/gaugefix.jl | 114 ++++++++++++++++++ test/bondenv/benv_fu.jl | 36 ++++++ test/bondenv/benv_gaugefix.jl | 43 +++++++ test/runtests.jl | 6 + 6 files changed, 278 insertions(+) create mode 100644 src/algorithms/contractions/bondenv/benv_ctm.jl create mode 100644 src/algorithms/contractions/bondenv/gaugefix.jl create mode 100644 test/bondenv/benv_fu.jl create mode 100644 test/bondenv/benv_gaugefix.jl diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index 4259e4b55..02253a719 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -45,7 +45,9 @@ include("algorithms/contractions/ctmrg_contractions.jl") include("algorithms/contractions/localoperator.jl") include("algorithms/contractions/vumps_contractions.jl") include("algorithms/contractions/bondenv/benv_tools.jl") +include("algorithms/contractions/bondenv/gaugefix.jl") include("algorithms/contractions/bondenv/als_solve.jl") +include("algorithms/contractions/bondenv/benv_ctm.jl") include("algorithms/ctmrg/sparse_environments.jl") include("algorithms/ctmrg/ctmrg.jl") diff --git a/src/algorithms/contractions/bondenv/benv_ctm.jl b/src/algorithms/contractions/bondenv/benv_ctm.jl new file mode 100644 index 000000000..6fc856b14 --- /dev/null +++ b/src/algorithms/contractions/bondenv/benv_ctm.jl @@ -0,0 +1,77 @@ +""" +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 +``` +where `XX = X' X` and `YY = Y' Y` (stacked together). + +Axis order: `[DX1 DY1; DX0 DY0]`, as in +``` + ┌---------------------┐ + | ┌----┐ | + └-| |-- DX0 DY0 --┘ + |benv| + ┌-| |-- DX1 DY1 --┐ + | └----┘ | + └---------------------┘ +``` +""" +function bondenv_fu(row::Int, col::Int, X::PEPSOrth, Y::PEPSOrth, env::CTMRGEnv) + Nr, Nc = size(env.corners)[[2, 3]] + cm1 = _prev(col, Nc) + cp1 = _next(col, Nc) + cp2 = _next(cp1, Nc) + rm1 = _prev(row, Nr) + rp1 = _next(row, Nr) + c1 = env.corners[1, rm1, cm1] + 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] + #= index labels + + left half right half + C1 -χ4 - T1 ------- χ6 ------- T1 - χ8 - C2 r-1 + | ‖ ‖ | + χ2 DNX DNY χ10 + | ‖ ‖ | + T4 =DWX= XX = DX = = DY = YY =DEY= T2 r + | ‖ ‖ | + χ1 DSX DSY χ9 + | ‖ ‖ | + C4 -χ3 - T3 ------- χ5 ------- T3 - χ7 - C3 r+1 + c-1 c c+1 c+2 + =# + # left half + @autoopt @tensor lhalf[DX1 DX0 χ5 χ6] := ( + 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] + ) + # right half + @autoopt @tensor rhalf[DY1 DY0 χ5 χ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] + ) + # combine + @autoopt @tensor benv[DX1 DY1; DX0 DY0] := (lhalf[DX1 DX0 χ5 χ6] * rhalf[DY1 DY0 χ5 χ6]) + normalize!(benv, Inf) + return benv +end diff --git a/src/algorithms/contractions/bondenv/gaugefix.jl b/src/algorithms/contractions/bondenv/gaugefix.jl new file mode 100644 index 000000000..00458b7bd --- /dev/null +++ b/src/algorithms/contractions/bondenv/gaugefix.jl @@ -0,0 +1,114 @@ +""" +Replace bond environment `benv` by its positive approximant `Z† Z` +(returns the "half environment" `Z`) +``` + ┌-----------------┐ ┌---------------┐ + | ┌----┐ | | | + └-| |-- 3 4 --┘ └-- Z -- 3 4 --┘ + |benv| = ↓ + ┌-| |-- 1 2 --┐ ┌-- Z†-- 1 2 --┐ + | └----┘ | | | + └-----------------┘ └---------------┘ +``` +""" +function positive_approx(benv::BondEnv) + # eigen-decomposition: benv = U D U' + D, U = eigh((benv + benv') / 2) + # determine if `env` is (mostly) positive or negative + sgn = sign(mean(D.data)) + # When optimizing the truncation of a bond, + # its environment can always be multiplied by a number. + # If `benv` is negative (e.g. obtained approximately from CTMRG), + # we can multiply it by (-1). + (sgn == -1) && (D *= -1) + # set (small) negative eigenvalues to 0 + D.data .= max.(D.data, 0) + Z = sdiag_pow(D, 0.5) * U' + return Z +end + +""" +Use QR decomposition to fix gauge of the half bond environment `Z`. +The reduced bond tensors `a`, `b` and `Z` are arranged as +``` + ┌---------------┐ + | | + └---Z---a---b---┘ + | ↓ ↓ + ↓ +``` +Reference: Physical Review B 92, 035142 (2015) +""" +function fixgauge_benv( + Z::AbstractTensorMap{T,S,1,2}, + a::AbstractTensorMap{T,S,1,2}, + b::AbstractTensorMap{T,S,2,1}, +) where {T<:Number,S<:ElementarySpace} + @assert !isdual(space(Z, 1)) + @assert !isdual(space(a, 2)) + @assert !isdual(space(b, 2)) + #= QR/LQ decomposition of Z + + 3 - Z - 2 = 2 - L - 1 3 - QL - 1 + ↓ ↓ + 1 2 + + = 1 - QR - 3 1 - R - 2 + ↓ + 2 + =# + QL, L = leftorth(Z, ((2, 1), (3,))) + QR, R = leftorth(Z, ((3, 1), (2,))) + Linv, Rinv = inv(L), inv(R) + #= fix gauge of Z, a, b + ┌---------------------------------------┐ + | | + └---Z---Rinv)---(R--a)--(b--L)---(Linv--┘ + | ↓ ↓ + ↓ + + -1 - R - 1 - a - -3 -1 - b - 1 - L - -3 + ↓ ↓ + -2 -2 + + ┌-----------------------------------------┐ + | | + └---Z-- 1 --Rinv-- -2 -3 --Linv-- 2 -┘ + ↓ + -1 + =# + @plansor a[-1; -2 -3] := R[-1; 1] * a[1; -2 -3] + @plansor b[-1 -2; -3] := b[-1 -2; 1] * L[-3; 1] + @plansor Z[-1; -2 -3] := Z[-1; 1 2] * Rinv[1; -2] * Linv[2; -3] + (isdual(space(R, 1)) == isdual(space(R, 2))) && twist!(a, 1) + (isdual(space(L, 1)) == isdual(space(L, 2))) && twist!(b, 3) + return Z, a, b, (Linv, Rinv) +end + +""" +When the (half) bond environment `Z` consists of two `PEPSOrth` tensors `X`, `Y` as +``` + ┌---------------┐ ┌-------------------┐ + | | = | | , + └---Z-- --┘ └--Z0---X-- --Y--┘ + ↓ ↓ +``` +apply the gauge transformation `Linv`, `Rinv` for `Z` to `X`, `Y`: +``` + -1 -1 + | | + -4 - X - 1 - Rinv - -2 -4 - Linv - 1 - Y - -2 + | | + -3 -3 +``` +""" +function _fixgauge_benvXY( + X::PEPSOrth{T,S}, + Y::PEPSOrth{T,S}, + Linv::AbstractTensorMap{T,S,1,1}, + Rinv::AbstractTensorMap{T,S,1,1}, +) where {T<:Number,S<:ElementarySpace} + @plansor X[-1 -2 -3 -4] := X[-1 1 -3 -4] * Rinv[1; -2] + @plansor Y[-1 -2 -3 -4] := Y[-1 -2 -3 1] * Linv[1; -4] + return X, Y +end diff --git a/test/bondenv/benv_fu.jl b/test/bondenv/benv_fu.jl new file mode 100644 index 000000000..c95171c30 --- /dev/null +++ b/test/bondenv/benv_fu.jl @@ -0,0 +1,36 @@ +using Test +using TensorKit +using PEPSKit +using LinearAlgebra +using KrylovKit +using Random + +Nr, Nc = 2, 3 +# create random PEPS +Random.seed!(0) +Pspace = Vect[FermionParity](0 => 1, 1 => 1) +Nspace = Vect[FermionParity](0 => 2, 1 => 2) +peps = InfinitePEPS(randn, ComplexF64, Pspace, Nspace; unitcell=(Nr, Nc)) +normalize!.(peps.A, Inf) +# calculate CTMRG environment +Envspace = Vect[FermionParity](0 => 3, 1 => 3) +ctm_alg = SequentialCTMRG(; tol=1e-10, verbosity=2, trscheme=FixedSpaceTruncation()) +env, = leading_boundary(CTMRGEnv(rand, ComplexF64, peps, Envspace), peps, ctm_alg) +for row in 1:Nr, col in 1:Nc + cp1 = PEPSKit._next(1, Nc) + A, B = peps.A[row, col], peps.A[row, cp1] + X, a, b, Y = PEPSKit._qr_bond(A, B) + # verify that gauge fixing can reduce condition number + benv = PEPSKit.bondenv_fu(row, col, X, Y, env) + @assert [isdual(space(benv, ax)) for ax in 1:numind(benv)] == [0, 0, 1, 1] + @tensor ab[DX DY; da db] := a[DX; da D] * b[D db; DY] + nrm1 = PEPSKit.inner_prod(benv, ab, ab) + # gauge fixing + Z = PEPSKit.positive_approx(benv) + cond = PEPSKit._condition_number(Z' * Z) + Z2, a2, b2, (Linv, Rinv) = PEPSKit.fixgauge_benv(Z, a, b) + benv2 = Z2' * Z2 + cond2 = PEPSKit._condition_number(benv2) + @test 1 <= cond2 < cond + @info "benv cond number: (gauge-fixed) $(cond2) ≤ $(cond) (initial)" +end diff --git a/test/bondenv/benv_gaugefix.jl b/test/bondenv/benv_gaugefix.jl new file mode 100644 index 000000000..bdd5cecf7 --- /dev/null +++ b/test/bondenv/benv_gaugefix.jl @@ -0,0 +1,43 @@ +using Test +using TensorKit +using PEPSKit +using LinearAlgebra +using KrylovKit +using Random + +Vphy = Vect[FermionParity ⊠ U1Irrep]((0, 0) => 1, (1, 1) => 1, (1, -1) => 2) +Vin = Vect[FermionParity ⊠ U1Irrep]((0, 0) => 1, (1, 1) => 3, (1, -1) => 2) +V = Vect[FermionParity ⊠ U1Irrep]((0, 0) => 1, (1, 1) => 2, (1, -1) => 3) +Vs = (V, V') +for V1 in Vs, V2 in Vs, V3 in Vs + #= + ┌---┬---------------┬---┐ + | | | | ┌--------------┐ + ├---X--- -2 -3 ---Y---┤ = | | + | | | | └--Z-- -2 -3 -┘ + └---┴-------Z0------┴---┘ ↓ + ↓ -1 + -1 + =# + X = rand(ComplexF64, Vin ⊗ V1' ⊗ Vin' ⊗ Vin) + Y = rand(ComplexF64, Vin ⊗ Vin ⊗ Vin' ⊗ V3) + Z0 = randn(ComplexF64, Vphy ← Vin ⊗ Vin' ⊗ Vin ⊗ Vin ⊗ Vin ⊗ Vin') + @tensor Z[p; Xe Yw] := Z0[p; Xn Xs Xw Yn Ye Ys] * X[Xn Xe Xs Xw] * Y[Yn Ye Ys Yw] + #= + ┌---------------------------┐ + | | + └---Z-- 1 --a-- 2 --b-- 3 --┘ + ↓ ↓ ↓ + -1 -2 -3 + =# + a = randn(ComplexF64, V1 ← Vphy' ⊗ V2) + b = randn(ComplexF64, V2 ⊗ Vphy ← V3) + @tensor half[:] := Z[-1; 1 3] * a[1; -2 2] * b[2 -3; 3] + Z2, a2, b2, (Linv, Rinv) = PEPSKit.fixgauge_benv(Z, a, b) + @tensor half2[:] := Z2[-1; 1 3] * a2[1; -2 2] * b2[2 -3; 3] + @test half ≈ half2 + # test gauge transformation of X, Y + X2, Y2 = PEPSKit._fixgauge_benvXY(X, Y, Linv, Rinv) + @tensor Z2_[p; Xe Yw] := Z0[p; Xn Xs Xw Yn Ye Ys] * X2[Xn Xe Xs Xw] * Y2[Yn Ye Ys Yw] + @test Z2 ≈ Z2_ +end diff --git a/test/runtests.jl b/test/runtests.jl index 2ea46b4fc..25c6b9774 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -50,6 +50,12 @@ end @time @safetestset "Iterative optimization after truncation" begin include("bondenv/bond_truncate.jl") end + @time @safetestset "Gauge fixing" begin + include("bondenv/benv_gaugefix.jl") + end + @time @safetestset "Full bond environment from CTMRG" begin + include("bondenv/benv_fu.jl") + end end if GROUP == "ALL" || GROUP == "TIMEEVOL" @time @safetestset "Cluster truncation with projectors" begin From b0bed4c65eba49ffb47d85aa72d85b5b16b9b3f7 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Mon, 26 May 2025 15:16:18 +0800 Subject: [PATCH 02/11] Change _condition_number definition --- src/utility/svd.jl | 16 ++++++++++++---- test/runtests.jl | 2 +- 2 files changed, 13 insertions(+), 5 deletions(-) diff --git a/src/utility/svd.jl b/src/utility/svd.jl index d4f454fca..ef6695c13 100644 --- a/src/utility/svd.jl +++ b/src/utility/svd.jl @@ -147,12 +147,20 @@ end ## Forward algorithms # Compute condition number smax / smin for diagonal singular value TensorMap -function _condition_number(S::AbstractTensorMap) +function _condition_number_sdiag(S::AbstractTensorMap) smax = maximum(first ∘ last, blocks(S)) smin = maximum(last ∘ last, blocks(S)) return smax / smin end +# Compute condition number smax / smin for any TensorMap +function _condition_number(a::AbstractTensorMap) + s = tsvd(a)[2] + smax = maximum(s.data) + smin = minimum(s.data) + return smax / smin +end + # Copy code from TensorKit but additionally return full U, S and V to make compatible with :fixed mode function _tsvd!( t::TensorMap{<:RealOrComplexFloat}, @@ -176,7 +184,7 @@ function _tsvd!( end # construct info NamedTuple - condnum = _condition_number(S) + condnum = _condition_number_sdiag(S) info = (; truncation_error=truncerr, condition_number=condnum, U_full=U, S_full=S, V_full=V⁺ ) @@ -217,7 +225,7 @@ end function _tsvd!(_, alg::FixedSVD, ::TruncationScheme, ::Real) info = (; truncation_error=0, - condition_number=_condition_number(alg.S), + condition_number=_condition_number_sdiag(alg.S), U_full=alg.U_full, S_full=alg.S_full, V_full=alg.V_full, @@ -272,7 +280,7 @@ function _tsvd!(f, alg::IterSVD, trunc::TruncationScheme, p::Real) trunc isa NoTruncation ? abs(zero(scalartype(f))) : norm(U * S * V - f, p) # construct info NamedTuple - condition_number = _condition_number(S) + condition_number = _condition_number_sdiag(S) info = (; truncation_error, condition_number, U_full=nothing, S_full=nothing, V_full=nothing ) diff --git a/test/runtests.jl b/test/runtests.jl index 25c6b9774..29f8fc0e3 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -53,7 +53,7 @@ end @time @safetestset "Gauge fixing" begin include("bondenv/benv_gaugefix.jl") end - @time @safetestset "Full bond environment from CTMRG" begin + @time @safetestset "Full update bond environment" begin include("bondenv/benv_fu.jl") end end From 4692ad98c87d3bd9c0bc38053eff23d9a36be663 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Wed, 28 May 2025 11:25:20 +0800 Subject: [PATCH 03/11] Construct bondenv_fu in a single @tensor call --- .../contractions/bondenv/benv_ctm.jl | 35 ++++++++----------- 1 file changed, 14 insertions(+), 21 deletions(-) diff --git a/src/algorithms/contractions/bondenv/benv_ctm.jl b/src/algorithms/contractions/bondenv/benv_ctm.jl index 6fc856b14..549440d74 100644 --- a/src/algorithms/contractions/bondenv/benv_ctm.jl +++ b/src/algorithms/contractions/bondenv/benv_ctm.jl @@ -1,12 +1,12 @@ """ 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 + C1---T1---------T1---C2 r-1 + | ‖ ‖ | + T4===XX== ==YY===T2 r + | ‖ ‖ | + C4---T3---------T3---C3 r+1 + c-1 c c+1 c+2 ``` where `XX = X' X` and `YY = Y' Y` (stacked together). @@ -14,9 +14,9 @@ Axis order: `[DX1 DY1; DX0 DY0]`, as in ``` ┌---------------------┐ | ┌----┐ | - └-| |-- DX0 DY0 --┘ - |benv| - ┌-| |-- DX1 DY1 --┐ + └-| |---DX0 DY0---┘ + |benv| + ┌-| |---DX1 DY1---┐ | └----┘ | └---------------------┘ ``` @@ -38,30 +38,25 @@ function bondenv_fu(row::Int, col::Int, X::PEPSOrth, Y::PEPSOrth, env::CTMRGEnv) t4 = env.edges[4, row, cm1] #= index labels - left half right half - C1 -χ4 - T1 ------- χ6 ------- T1 - χ8 - C2 r-1 + C1--χ4--T1X---------χ6---------T1Y--χ8---C2 r-1 | ‖ ‖ | χ2 DNX DNY χ10 | ‖ ‖ | - T4 =DWX= XX = DX = = DY = YY =DEY= T2 r + T4==DWX==XX===DX== ==DY===YY==DEY==T2 r | ‖ ‖ | χ1 DSX DSY χ9 | ‖ ‖ | - C4 -χ3 - T3 ------- χ5 ------- T3 - χ7 - C3 r+1 + C4--χ3--T3X---------χ5---------T3Y--χ7---C3 r+1 c-1 c c+1 c+2 =# - # left half - @autoopt @tensor lhalf[DX1 DX0 χ5 χ6] := ( + @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] - ) - # right half - @autoopt @tensor rhalf[DY1 DY0 χ5 χ6] := ( + t1X[χ4 DNX0 DNX1 χ6] * c3[χ9 χ7] * t2[χ10 DEY0 DEY1 χ9] * c2[χ8 χ10] * @@ -70,8 +65,6 @@ function bondenv_fu(row::Int, col::Int, X::PEPSOrth, Y::PEPSOrth, env::CTMRGEnv) conj(Y[DNY1 DEY1 DSY1 DY1]) * t1Y[χ6 DNY0 DNY1 χ8] ) - # combine - @autoopt @tensor benv[DX1 DY1; DX0 DY0] := (lhalf[DX1 DX0 χ5 χ6] * rhalf[DY1 DY0 χ5 χ6]) normalize!(benv, Inf) return benv end From c593eca98076853e27800903fb6c53e9ff95d95d Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Wed, 28 May 2025 11:25:41 +0800 Subject: [PATCH 04/11] Improve `positive_approx` --- src/algorithms/contractions/bondenv/gaugefix.jl | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/algorithms/contractions/bondenv/gaugefix.jl b/src/algorithms/contractions/bondenv/gaugefix.jl index 00458b7bd..39c8cbc80 100644 --- a/src/algorithms/contractions/bondenv/gaugefix.jl +++ b/src/algorithms/contractions/bondenv/gaugefix.jl @@ -15,15 +15,17 @@ function positive_approx(benv::BondEnv) # eigen-decomposition: benv = U D U' D, U = eigh((benv + benv') / 2) # determine if `env` is (mostly) positive or negative - sgn = sign(mean(D.data)) + sgn = sign(sum(D.data)) # When optimizing the truncation of a bond, # its environment can always be multiplied by a number. # If `benv` is negative (e.g. obtained approximately from CTMRG), # we can multiply it by (-1). - (sgn == -1) && (D *= -1) - # set (small) negative eigenvalues to 0 - D.data .= max.(D.data, 0) - Z = sdiag_pow(D, 0.5) * U' + data = D.data + @inbounds for i in eachindex(data) + d = sgn == -1 ? -data[i] : data[i] + data[i] = ifelse(d > 0, sqrt(d), zero(d)) + end + Z = D * U' return Z end From 13b7944315ed0c92959d30129b65b59e8c0ef33e Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Wed, 28 May 2025 12:06:08 +0800 Subject: [PATCH 05/11] Expand benv_fu test --- test/bondenv/benv_fu.jl | 38 +++++++++++++++++++++++++------------- 1 file changed, 25 insertions(+), 13 deletions(-) diff --git a/test/bondenv/benv_fu.jl b/test/bondenv/benv_fu.jl index c95171c30..7915b5df5 100644 --- a/test/bondenv/benv_fu.jl +++ b/test/bondenv/benv_fu.jl @@ -5,32 +5,44 @@ using LinearAlgebra using KrylovKit using Random -Nr, Nc = 2, 3 -# create random PEPS Random.seed!(0) -Pspace = Vect[FermionParity](0 => 1, 1 => 1) -Nspace = Vect[FermionParity](0 => 2, 1 => 2) -peps = InfinitePEPS(randn, ComplexF64, Pspace, Nspace; unitcell=(Nr, Nc)) -normalize!.(peps.A, Inf) +Nr, Nc = 2, 2 +# create Hubbard iPEPS using simple update +function get_hubbard_state(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) + wpeps = InfiniteWeightPEPS(rand, ComplexF64, Vphy, Vphy; unitcell=(Nr, Nc)) + alg = SimpleUpdate(1e-2, 1e-7, 10000, truncerr(1e-10) & truncdim(4)) + wpeps, = simpleupdate(wpeps, H, alg; bipartite=false) + peps = InfinitePEPS(wpeps) + normalize!.(peps.A, Inf) + return peps +end + +peps = get_hubbard_state() # calculate CTMRG environment -Envspace = Vect[FermionParity](0 => 3, 1 => 3) -ctm_alg = SequentialCTMRG(; tol=1e-10, verbosity=2, trscheme=FixedSpaceTruncation()) +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=1e-10, verbosity=2, trscheme=truncerr(1e-10) & truncdim(8)) env, = leading_boundary(CTMRGEnv(rand, ComplexF64, peps, Envspace), peps, ctm_alg) for row in 1:Nr, col in 1:Nc - cp1 = PEPSKit._next(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) - # verify that gauge fixing can reduce condition number benv = PEPSKit.bondenv_fu(row, col, X, Y, env) @assert [isdual(space(benv, ax)) for ax in 1:numind(benv)] == [0, 0, 1, 1] - @tensor ab[DX DY; da db] := a[DX; da D] * b[D db; DY] - nrm1 = PEPSKit.inner_prod(benv, ab, ab) - # gauge fixing Z = PEPSKit.positive_approx(benv) + # verify that gauge fixing can greatly reduce + # condition number for physical state bond envs cond = PEPSKit._condition_number(Z' * Z) Z2, a2, b2, (Linv, Rinv) = PEPSKit.fixgauge_benv(Z, a, b) benv2 = Z2' * Z2 cond2 = PEPSKit._condition_number(benv2) @test 1 <= cond2 < cond @info "benv cond number: (gauge-fixed) $(cond2) ≤ $(cond) (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 From 2f635d775924bd85fc6c5777e1853935e713ccb2 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Wed, 28 May 2025 22:43:54 +0800 Subject: [PATCH 06/11] Use `LinearAlgebra.cond` in latest TensorKit --- Project.toml | 2 +- src/algorithms/contractions/bondenv/gaugefix.jl | 6 ++++-- src/utility/svd.jl | 16 ++++------------ test/bondenv/benv_fu.jl | 4 ++-- 4 files changed, 11 insertions(+), 17 deletions(-) diff --git a/Project.toml b/Project.toml index a0e8fe041..cb1675223 100644 --- a/Project.toml +++ b/Project.toml @@ -41,7 +41,7 @@ Printf = "1" QuadGK = "2.11.1" Random = "1" Statistics = "1" -TensorKit = "0.14" +TensorKit = "0.14.6" TensorOperations = "5" VectorInterface = "0.4, 0.5" Zygote = "0.6, 0.7" diff --git a/src/algorithms/contractions/bondenv/gaugefix.jl b/src/algorithms/contractions/bondenv/gaugefix.jl index 39c8cbc80..abf9eece5 100644 --- a/src/algorithms/contractions/bondenv/gaugefix.jl +++ b/src/algorithms/contractions/bondenv/gaugefix.jl @@ -22,8 +22,8 @@ function positive_approx(benv::BondEnv) # we can multiply it by (-1). data = D.data @inbounds for i in eachindex(data) - d = sgn == -1 ? -data[i] : data[i] - data[i] = ifelse(d > 0, sqrt(d), zero(d)) + d = (sgn == -1) ? -data[i] : data[i] + data[i] = (d > 0) ? sqrt(d) : zero(d) end Z = D * U' return Z @@ -61,6 +61,8 @@ function fixgauge_benv( =# QL, L = leftorth(Z, ((2, 1), (3,))) QR, R = leftorth(Z, ((3, 1), (2,))) + @debug "Condition number of L = $(LinearAlgebra.cond(L))" + @debug "Condition number of L = $(LinearAlgebra.cond(R))" Linv, Rinv = inv(L), inv(R) #= fix gauge of Z, a, b ┌---------------------------------------┐ diff --git a/src/utility/svd.jl b/src/utility/svd.jl index ef6695c13..d4f454fca 100644 --- a/src/utility/svd.jl +++ b/src/utility/svd.jl @@ -147,20 +147,12 @@ end ## Forward algorithms # Compute condition number smax / smin for diagonal singular value TensorMap -function _condition_number_sdiag(S::AbstractTensorMap) +function _condition_number(S::AbstractTensorMap) smax = maximum(first ∘ last, blocks(S)) smin = maximum(last ∘ last, blocks(S)) return smax / smin end -# Compute condition number smax / smin for any TensorMap -function _condition_number(a::AbstractTensorMap) - s = tsvd(a)[2] - smax = maximum(s.data) - smin = minimum(s.data) - return smax / smin -end - # Copy code from TensorKit but additionally return full U, S and V to make compatible with :fixed mode function _tsvd!( t::TensorMap{<:RealOrComplexFloat}, @@ -184,7 +176,7 @@ function _tsvd!( end # construct info NamedTuple - condnum = _condition_number_sdiag(S) + condnum = _condition_number(S) info = (; truncation_error=truncerr, condition_number=condnum, U_full=U, S_full=S, V_full=V⁺ ) @@ -225,7 +217,7 @@ end function _tsvd!(_, alg::FixedSVD, ::TruncationScheme, ::Real) info = (; truncation_error=0, - condition_number=_condition_number_sdiag(alg.S), + condition_number=_condition_number(alg.S), U_full=alg.U_full, S_full=alg.S_full, V_full=alg.V_full, @@ -280,7 +272,7 @@ function _tsvd!(f, alg::IterSVD, trunc::TruncationScheme, p::Real) trunc isa NoTruncation ? abs(zero(scalartype(f))) : norm(U * S * V - f, p) # construct info NamedTuple - condition_number = _condition_number_sdiag(S) + condition_number = _condition_number(S) info = (; truncation_error, condition_number, U_full=nothing, S_full=nothing, V_full=nothing ) diff --git a/test/bondenv/benv_fu.jl b/test/bondenv/benv_fu.jl index 7915b5df5..ffa4af8d9 100644 --- a/test/bondenv/benv_fu.jl +++ b/test/bondenv/benv_fu.jl @@ -35,10 +35,10 @@ for row in 1:Nr, col in 1:Nc Z = PEPSKit.positive_approx(benv) # verify that gauge fixing can greatly reduce # condition number for physical state bond envs - cond = PEPSKit._condition_number(Z' * Z) + cond = LinearAlgebra.cond(Z' * Z) Z2, a2, b2, (Linv, Rinv) = PEPSKit.fixgauge_benv(Z, a, b) benv2 = Z2' * Z2 - cond2 = PEPSKit._condition_number(benv2) + cond2 = LinearAlgebra.cond(benv2) @test 1 <= cond2 < cond @info "benv cond number: (gauge-fixed) $(cond2) ≤ $(cond) (initial)" # verify gauge fixing is done correctly From db4d3028657b2a828a34e6640b0513047e4b3655 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Wed, 28 May 2025 22:46:29 +0800 Subject: [PATCH 07/11] Fix `smin` in `_condition_number` --- src/utility/svd.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utility/svd.jl b/src/utility/svd.jl index d4f454fca..d174b461f 100644 --- a/src/utility/svd.jl +++ b/src/utility/svd.jl @@ -149,7 +149,7 @@ end # Compute condition number smax / smin for diagonal singular value TensorMap function _condition_number(S::AbstractTensorMap) smax = maximum(first ∘ last, blocks(S)) - smin = maximum(last ∘ last, blocks(S)) + smin = minimum(last ∘ last, blocks(S)) return smax / smin end From b9acac741087ba6bfc05530cf6d3b253181d6e25 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Thu, 29 May 2025 20:32:53 +0800 Subject: [PATCH 08/11] Update debug messages --- src/algorithms/contractions/bondenv/gaugefix.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/algorithms/contractions/bondenv/gaugefix.jl b/src/algorithms/contractions/bondenv/gaugefix.jl index abf9eece5..70a5b0dc8 100644 --- a/src/algorithms/contractions/bondenv/gaugefix.jl +++ b/src/algorithms/contractions/bondenv/gaugefix.jl @@ -61,8 +61,7 @@ function fixgauge_benv( =# QL, L = leftorth(Z, ((2, 1), (3,))) QR, R = leftorth(Z, ((3, 1), (2,))) - @debug "Condition number of L = $(LinearAlgebra.cond(L))" - @debug "Condition number of L = $(LinearAlgebra.cond(R))" + @debug "cond(L) = $(LinearAlgebra.cond(L)); cond(R) = $(LinearAlgebra.cond(R))" Linv, Rinv = inv(L), inv(R) #= fix gauge of Z, a, b ┌---------------------------------------┐ From 2ac1bf2007793e0b00eb5675644424c2eb3d522f Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Thu, 29 May 2025 20:41:44 +0800 Subject: [PATCH 09/11] Adjust test output --- test/bondenv/benv_fu.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/bondenv/benv_fu.jl b/test/bondenv/benv_fu.jl index ffa4af8d9..92f310c15 100644 --- a/test/bondenv/benv_fu.jl +++ b/test/bondenv/benv_fu.jl @@ -5,15 +5,15 @@ using LinearAlgebra using KrylovKit using Random -Random.seed!(0) +Random.seed!(100) Nr, Nc = 2, 2 # create Hubbard iPEPS using simple update function get_hubbard_state(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) wpeps = InfiniteWeightPEPS(rand, ComplexF64, Vphy, Vphy; unitcell=(Nr, Nc)) - alg = SimpleUpdate(1e-2, 1e-7, 10000, truncerr(1e-10) & truncdim(4)) - wpeps, = simpleupdate(wpeps, H, alg; bipartite=false) + alg = SimpleUpdate(1e-2, 1e-8, 10000, truncerr(1e-10) & truncdim(4)) + wpeps, = simpleupdate(wpeps, H, alg; bipartite=false, check_interval=2000) peps = InfinitePEPS(wpeps) normalize!.(peps.A, Inf) return peps From a38820862d8c6d0e897b5e62ee25fd829fdede80 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Sun, 1 Jun 2025 09:23:52 +0800 Subject: [PATCH 10/11] Avoid `cond` name conflict in test --- test/bondenv/benv_fu.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/bondenv/benv_fu.jl b/test/bondenv/benv_fu.jl index 92f310c15..52431a855 100644 --- a/test/bondenv/benv_fu.jl +++ b/test/bondenv/benv_fu.jl @@ -35,12 +35,12 @@ for row in 1:Nr, col in 1:Nc Z = PEPSKit.positive_approx(benv) # verify that gauge fixing can greatly reduce # condition number for physical state bond envs - cond = LinearAlgebra.cond(Z' * Z) + cond1 = cond(Z' * Z) Z2, a2, b2, (Linv, Rinv) = PEPSKit.fixgauge_benv(Z, a, b) benv2 = Z2' * Z2 - cond2 = LinearAlgebra.cond(benv2) - @test 1 <= cond2 < cond - @info "benv cond number: (gauge-fixed) $(cond2) ≤ $(cond) (initial)" + 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] From 18d91691127c6f566b1f2657b48555204fa3020f Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Sun, 1 Jun 2025 09:24:11 +0800 Subject: [PATCH 11/11] Add TODO to improve gauge fixing method --- src/algorithms/contractions/bondenv/gaugefix.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/algorithms/contractions/bondenv/gaugefix.jl b/src/algorithms/contractions/bondenv/gaugefix.jl index 70a5b0dc8..3a2a17d65 100644 --- a/src/algorithms/contractions/bondenv/gaugefix.jl +++ b/src/algorithms/contractions/bondenv/gaugefix.jl @@ -39,7 +39,9 @@ The reduced bond tensors `a`, `b` and `Z` are arranged as | ↓ ↓ ↓ ``` -Reference: Physical Review B 92, 035142 (2015) +Reference: +- Physical Review B 90, 064425 (2014) +- Physical Review B 92, 035142 (2015) """ function fixgauge_benv( Z::AbstractTensorMap{T,S,1,2}, @@ -62,6 +64,7 @@ function fixgauge_benv( QL, L = leftorth(Z, ((2, 1), (3,))) QR, R = leftorth(Z, ((3, 1), (2,))) @debug "cond(L) = $(LinearAlgebra.cond(L)); cond(R) = $(LinearAlgebra.cond(R))" + # TODO: find a better way to fix gauge that avoids `inv` Linv, Rinv = inv(L), inv(R) #= fix gauge of Z, a, b ┌---------------------------------------┐