Improve bond truncation algorithms (#303)

* Change bond truncation `tol` to normalized bond SVD diff

* Remove bond s normalization in FET

* Fix formatting

* Reduce repeated normalization of `s` when calculating `Δs`

* Update docstring

Author: Yue-Zhengyuan

src/algorithms/truncation/bond_truncation.jl

@@ -15,22 +15,23 @@ The truncation algorithm can be constructed from the following keyword arguments
 
 * `trunc::TruncationStrategy`: SVD truncation strategy when initilizing the truncated tensors connected by the bond.
 * `maxiter::Int=50` : Maximal number of ALS iterations.
-* `tol::Float64=1e-15` : ALS converges when fidelity change between two iterations is smaller than `tol`.
+* `tol::Float64=1e-9` : ALS converges when the relative change in bond SVD spectrum between two iterations is smaller than `tol`.
 * `check_interval::Int=0` : Set number of iterations to print information. Output is suppressed when `check_interval <= 0`. 
 """
 @kwdef struct ALSTruncation
     trunc::TruncationStrategy
     maxiter::Int = 50
-    tol::Float64 = 1.0e-15
+    tol::Float64 = 1.0e-9
     check_interval::Int = 0
 end
 
 function _als_message(
-        iter::Int, cost::Float64, fid::Float64, Δcost::Float64, Δfid::Float64, time_elapsed::Float64,
+        iter::Int, cost::Float64, fid::Float64, Δcost::Float64,
+        Δfid::Float64, Δs::Float64, time_elapsed::Float64,
     )
     return @sprintf(
         "%5d, fid = %.8e, Δfid = %.8e, time = %.4f s\n", iter, fid, Δfid, time_elapsed
-    ) * @sprintf("      cost = %.3e,   Δcost/cost0 = %.3e", cost, Δcost)
+    ) * @sprintf("      cost = %.3e, Δcost/cost0 = %.3e, |Δs| = %.4e.", cost, Δcost, Δs)
 end
 
 """
@@ -65,14 +66,14 @@ function bond_truncate(
     @assert !isdual(space(a, 2))
     @assert !isdual(space(b, 2))
     @assert codomain(benv) == domain(benv)
+    need_flip = isdual(space(b, 1))
     time00 = time()
     verbose = (alg.check_interval > 0)
     a2b2 = _combine_ab(a, b)
     # initialize truncated a, b
     perm_ab = ((1, 3), (4, 2))
-    a, s, b = svd_trunc(permute(a2b2, perm_ab); trunc = alg.trunc)
-    s /= norm(s, Inf)
-    a, b = absorb_s(a, s, b)
+    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
             ↓
@@ -84,8 +85,8 @@ function bond_truncate(
     # cost function will be normalized by initial value
     cost00 = cost_function_als(benv, ab, a2b2)
     fid = fidelity(benv, ab, a2b2)
-    cost0, fid0, Δfid = cost00, fid, 0.0
-    verbose && @info "ALS init" * _als_message(0, cost0, fid, NaN, NaN, 0.0)
+    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)
         time0 = time()
         #= 
@@ -103,20 +104,27 @@ function bond_truncate(
         Rb = _tensor_Rb(benv, a)
         Sb = _tensor_Sb(benv, a, a2b2)
         b, info_b = _solve_ab(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)
+        # 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
+        s_nrm = norm(s0, Inf)
+        Δs = ((space(s) == space(s0)) ? _singular_value_distance((s, s0)) : NaN) / s_nrm
         Δcost = abs(cost - cost0) / cost00
         Δfid = abs(fid - fid0)
-        cost0, fid0 = cost, fid
+        cost0, fid0, s0 = cost, fid, s
         time1 = time()
-        converge = (Δfid < alg.tol)
+        converge = (Δs < alg.tol)
         cancel = (iter == alg.maxiter)
         showinfo =
             cancel || (verbose && (converge || iter == 1 || iter % alg.check_interval == 0))
         if showinfo
             message = _als_message(
-                iter, cost, fid, Δcost, Δfid, time1 - ((cancel || converge) ? time00 : time0),
+                iter, cost, fid, Δcost, Δfid, Δs,
+                time1 - ((cancel || converge) ? time00 : time0),
             )
             if converge
                 @info "ALS conv" * message
@@ -129,9 +137,11 @@ function bond_truncate(
         converge && break
     end
     a, s, b = svd_trunc!(permute(_combine_ab(a, b), perm_ab); trunc = alg.trunc)
-    # normalize singular value spectrum
-    s /= norm(s, Inf)
-    return a, s, b, (; fid, Δfid)
+    a, b = absorb_s(a, s, b)
+    if need_flip
+        a, s, b = flip_svd(a, s, b)
+    end
+    return a, s, b, (; fid, Δfid, Δs)
 end
 
 function bond_truncate(
@@ -144,18 +154,15 @@ function bond_truncate(
     @assert !isdual(space(a, 2))
     @assert !isdual(space(b, 2))
     @assert codomain(benv) == domain(benv)
+    need_flip = isdual(space(b, 1))
     #= initialize bond matrix using QR as `Ra Lb`
 
-        --- a == b ---   ==>   - Qa - Ra == Rb - Qb -
+        --- a == b ---   ==>   - Qa ← Ra == Rb ← Qb -
             ↓    ↓               ↓               ↓
     =#
     Qa, Ra = left_orth(a)
     Rb, Qb = right_orth(b)
-    # if Qa → Ra, a twist is needed to express a as
-    # contraction of Rb, Qb instead of Qa * Ra
-    isdual(space(Ra, 1)) && twist!(Ra, 1)
-    # similarly if Rb → Qb
-    isdual(space(Qb, 1)) && twist!(Rb, 2)
+    @assert !isdual(space(Ra, 1)) && !isdual(space(Qb, 1))
     @tensor b0[-1; -2] := Ra[-1 1] * Rb[1 -2]
     #= initialize bond environment around `Ra Lb`
 
@@ -174,8 +181,12 @@ function bond_truncate(
     )
     # optimize bond matrix
     u, s, vh, info = fullenv_truncate(b0, benv2, alg)
+    u, vh = absorb_s(u, s, vh)
     # truncate a, b tensors with u, s, vh
     @tensor a[-1 -2; -3] := Qa[-1 -2 3] * u[3 -3]
     @tensor b[-1; -2 -3] := vh[-1 1] * Qb[1 -2 -3]
+    if need_flip
+        a, s, b = flip_svd(a, s, b)
+    end
     return a, s, b, info
 end

src/algorithms/truncation/fullenv_truncation.jl

@@ -15,7 +15,7 @@ The truncation algorithm can be constructed from the following keyword arguments
 
 * `trunc::TruncationStrategy` : SVD truncation strategy when optimizing the new bond matrix.
 * `maxiter::Int=50` : Maximal number of FET iterations.
-* `tol::Float64=1e-15` : FET converges when fidelity change between two FET iterations is smaller than `tol`.
+* `tol::Float64=1e-9` : FET converges when the relative change in bond SVD spectrum between two FET iterations is smaller than `tol`.
 * `trunc_init::Bool=true` : Controls whether the initialization of the new bond matrix is obtained from truncated SVD of the old bond matrix. 
 * `check_interval::Int=0` : Set number of iterations to print information. Output is suppressed when `check_interval <= 0`. 
 
@@ -26,7 +26,7 @@ The truncation algorithm can be constructed from the following keyword arguments
 @kwdef struct FullEnvTruncation
     trunc::TruncationStrategy
     maxiter::Int = 50
-    tol::Float64 = 1.0e-15
+    tol::Float64 = 1.0e-9
     trunc_init::Bool = true
     check_interval::Int = 0
 end
@@ -230,8 +230,6 @@ function fullenv_truncate(
     # initialize u, s, vh with truncated or untruncated SVD
     u, s, vh = svd_trunc(b0; trunc = (alg.trunc_init ? alg.trunc : notrunc()))
     b1 = similar(b0)
-    # normalize `s` (bond matrices can always be normalized)
-    s /= norm(s, Inf)
     s0 = deepcopy(s)
     Δfid, Δs, fid, fid0 = NaN, NaN, 0.0, 0.0
     for iter in 1:(alg.maxiter)
@@ -245,25 +243,25 @@ function fullenv_truncate(
         r, info_r = linsolve(Base.Fix1(*, B), p, r, 0, 1)
         @tensor b1[-1; -2] = u[-1; 1] * r[1 -2]
         u, s, vh = svd_trunc(b1; trunc = alg.trunc)
-        s /= norm(s, Inf)
         # update `- l ←  =  - u ← s ←`
         @tensor l[-1 -2] := u[-1; 1] * s[1; -2]
         @tensor p[-1 -2] := conj(vh[-2; 2]) * benv[-1 2; 3 4] * b0[3; 4]
         @tensor B[-1 -2; -3 -4] := conj(vh[-2; 2]) * benv[-1 2; -3 4] * vh[-4; 4]
         _linearmap_twist!(p)
         _linearmap_twist!(B)
         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)
         u, s, vh = svd_trunc!(b1; trunc = alg.trunc)
-        s /= norm(s, Inf)
         # determine convergence
-        Δs = (space(s) == space(s0)) ? _singular_value_distance((s, s0)) : NaN
+        s_nrm = norm(s0, Inf)
+        Δs = ((space(s) == space(s0)) ? _singular_value_distance((s, s0)) : NaN) / s_nrm
         Δfid = fid - fid0
         s0 = deepcopy(s)
         fid0 = fid
         time1 = time()
-        converge = (Δfid < alg.tol)
+        converge = (Δs < alg.tol)
         cancel = (iter == alg.maxiter)
         showinfo =
             cancel || (verbose && (converge || iter == 1 || iter % alg.check_interval == 0))

test/bondenv/bond_truncate.jl

@@ -7,7 +7,7 @@ using LinearAlgebra
 using KrylovKit
 
 Random.seed!(0)
-maxiter = 500
+maxiter = 600
 check_interval = 20
 trunc = truncerror(; atol = 1.0e-10) & truncrank(8)
 Vext = Vect[FermionParity](0 => 100, 1 => 100)
@@ -35,8 +35,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])
-        a1, b1 = PEPSKit.absorb_s(a1, ss[label], b1)
+        # display(ss[label])
         @test info.fid ≈ PEPSKit.fidelity(benv, PEPSKit._combine_ab(a1, b1), a2b2)
         @test info.fid > fid0
     end