Report Approximate Truncation Error from Apply Function (#232)

Author: JoeyT1994

Project.toml

@@ -1,7 +1,7 @@
 name = "ITensorNetworks"
 uuid = "2919e153-833c-4bdc-8836-1ea460a35fc7"
 authors = ["Matthew Fishman <[email protected]>, Joseph Tindall <[email protected]> and contributors"]
-version = "0.13.4"
+version = "0.13.5"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"

src/apply.jl

@@ -34,7 +34,7 @@ function full_update_bp(
   nfullupdatesweeps=10,
   print_fidelity_loss=false,
   envisposdef=false,
-  (singular_values!)=nothing,
+  callback=Returns(nothing),
   symmetrize=false,
   apply_kwargs...,
 )
@@ -65,21 +65,23 @@ function full_update_bp(
     apply_kwargs...,
   )
   if symmetrize
-    Rᵥ₁, Rᵥ₂ = factorize_svd(
+    singular_values! = Ref(ITensor())
+    Rᵥ₁, Rᵥ₂, spec = factorize_svd(
       Rᵥ₁ * Rᵥ₂,
       inds(Rᵥ₁);
       ortho="none",
       tags=edge_tag(v⃗[1] => v⃗[2]),
       singular_values!,
       apply_kwargs...,
     )
+    callback(; singular_values=singular_values![], truncation_error=spec.truncerr)
   end
   ψᵥ₁ = Qᵥ₁ * Rᵥ₁
   ψᵥ₂ = Qᵥ₂ * Rᵥ₂
   return ψᵥ₁, ψᵥ₂
 end
 
-function simple_update_bp_full(o, ψ, v⃗; envs, (singular_values!)=nothing, apply_kwargs...)
+function simple_update_bp_full(o, ψ, v⃗; envs, callback=Returns(nothing), apply_kwargs...)
   cutoff = 10 * eps(real(scalartype(ψ)))
   envs_v1 = filter(env -> hascommoninds(env, ψ[v⃗[1]]), envs)
   envs_v2 = filter(env -> hascommoninds(env, ψ[v⃗[2]]), envs)
@@ -116,9 +118,11 @@ function simple_update_bp_full(o, ψ, v⃗; envs, (singular_values!)=nothing, ap
   v1_inds = [v1_inds; siteinds(ψ, v⃗[1])]
   v2_inds = [v2_inds; siteinds(ψ, v⃗[2])]
   e = v⃗[1] => v⃗[2]
-  ψᵥ₁, ψᵥ₂ = factorize_svd(
+  singular_values! = Ref(ITensor())
+  ψᵥ₁, ψᵥ₂, spec = factorize_svd(
     oψ, v1_inds; ortho="none", tags=edge_tag(e), singular_values!, apply_kwargs...
   )
+  callback(; singular_values=singular_values![], truncation_error=spec.truncerr)
   for inv_sqrt_env_v1 in inv_sqrt_envs_v1
     ψᵥ₁ *= dag(inv_sqrt_env_v1)
   end
@@ -129,7 +133,7 @@ function simple_update_bp_full(o, ψ, v⃗; envs, (singular_values!)=nothing, ap
 end
 
 # Reduced version
-function simple_update_bp(o, ψ, v⃗; envs, (singular_values!)=nothing, apply_kwargs...)
+function simple_update_bp(o, ψ, v⃗; envs, callback=Returns(nothing), apply_kwargs...)
   cutoff = 10 * eps(real(scalartype(ψ)))
   envs_v1 = filter(env -> hascommoninds(env, ψ[v⃗[1]]), envs)
   envs_v2 = filter(env -> hascommoninds(env, ψ[v⃗[2]]), envs)
@@ -164,14 +168,16 @@ function simple_update_bp(o, ψ, v⃗; envs, (singular_values!)=nothing, apply_k
   rᵥ₂ = commoninds(Qᵥ₂, Rᵥ₂)
   oR = apply(o, Rᵥ₁ * Rᵥ₂)
   e = v⃗[1] => v⃗[2]
-  Rᵥ₁, Rᵥ₂ = factorize_svd(
+  singular_values! = Ref(ITensor())
+  Rᵥ₁, Rᵥ₂, spec = factorize_svd(
     oR,
     unioninds(rᵥ₁, sᵥ₁);
     ortho="none",
     tags=edge_tag(e),
     singular_values!,
     apply_kwargs...,
   )
+  callback(; singular_values=singular_values![], truncation_error=spec.truncerr)
   Qᵥ₁ = contract([Qᵥ₁; dag.(inv_sqrt_envs_v1)])
   Qᵥ₂ = contract([Qᵥ₂; dag.(inv_sqrt_envs_v2)])
   ψᵥ₁ = Qᵥ₁ * Rᵥ₁
@@ -188,7 +194,7 @@ function ITensors.apply(
   nfullupdatesweeps=10,
   print_fidelity_loss=false,
   envisposdef=false,
-  (singular_values!)=nothing,
+  callback=Returns(nothing),
   variational_optimization_only=false,
   symmetrize=false,
   reduced=true,
@@ -224,15 +230,15 @@ function ITensors.apply(
         nfullupdatesweeps,
         print_fidelity_loss,
         envisposdef,
-        singular_values!,
+        callback,
         symmetrize,
         apply_kwargs...,
       )
     else
       if reduced
-        ψᵥ₁, ψᵥ₂ = simple_update_bp(o, ψ, v⃗; envs, singular_values!, apply_kwargs...)
+        ψᵥ₁, ψᵥ₂ = simple_update_bp(o, ψ, v⃗; envs, callback, apply_kwargs...)
       else
-        ψᵥ₁, ψᵥ₂ = simple_update_bp_full(o, ψ, v⃗; envs, singular_values!, apply_kwargs...)
+        ψᵥ₁, ψᵥ₂ = simple_update_bp_full(o, ψ, v⃗; envs, callback, apply_kwargs...)
       end
     end
     if normalize

test/test_apply.jl

@@ -11,7 +11,7 @@ using ITensorNetworks:
   random_tensornetwork,
   siteinds,
   update
-using ITensors: ITensors, inner, op
+using ITensors: ITensors, ITensor, inner, op
 using NamedGraphs.NamedGraphGenerators: named_grid
 using NamedGraphs.PartitionedGraphs: PartitionVertex
 using SplitApplyCombine: group
@@ -39,9 +39,15 @@ using Test: @test, @testset
   envsGBP = environment(bp_cache, [(v1, "bra"), (v1, "ket"), (v2, "bra"), (v2, "ket")])
   inner_alg = "exact"
   ngates = 5
+  truncerr = 0.0
+  singular_values = ITensor()
+  function callback(; singular_values, truncation_error)
+    truncerr = truncation_error
+    singular_values = singular_values
+  end
   for i in 1:ngates
     o = op("RandomUnitary", s[v1]..., s[v2]...)
-    ψOexact = apply(o, ψ; cutoff=1e-16)
+    ψOexact = apply(o, ψ; cutoff=nothing)
     ψOSBP = apply(
       o,
       ψ;
@@ -50,6 +56,7 @@ using Test: @test, @testset
       normalize=true,
       print_fidelity_loss=true,
       envisposdef=true,
+      callback,
     )
     ψOv = apply(o, ψv; maxdim=χ, normalize=true)
     ψOVidal_symm = ITensorNetwork(ψOv)
@@ -73,6 +80,7 @@ using Test: @test, @testset
     fGBP =
       inner(ψOGBP, ψOexact; alg=inner_alg) /
       sqrt(inner(ψOexact, ψOexact; alg=inner_alg) * inner(ψOGBP, ψOGBP; alg=inner_alg))
+    @test !iszero(truncerr)
     @test real(fGBP * conj(fGBP)) >= real(fSBP * conj(fSBP))
     @test isapprox(real(fSBP * conj(fSBP)), real(fVidal * conj(fVidal)); atol=1e-3)
   end