From c83da376eb5074589cc4097d5445ecd8c6422e95 Mon Sep 17 00:00:00 2001
From: Paul Brehmer <paul.brehmer@univie.ac.at>
Date: Thu, 11 Sep 2025 18:31:56 +0200
Subject: [PATCH 01/24] Add `maximize_fidelity!` (not yet working)

---
 src/PEPSKit.jl                                |  2 +
 .../time_evolution/fidelity_maximization.jl   | 90 +++++++++++++++++++
 2 files changed, 92 insertions(+)
 create mode 100644 src/algorithms/time_evolution/fidelity_maximization.jl

diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl
index 27c09f395..1e8fd2459 100644
--- a/src/PEPSKit.jl
+++ b/src/PEPSKit.jl
@@ -71,6 +71,7 @@ include("algorithms/truncation/bond_truncation.jl")
 include("algorithms/time_evolution/evoltools.jl")
 include("algorithms/time_evolution/simpleupdate.jl")
 include("algorithms/time_evolution/simpleupdate3site.jl")
+include("algorithms/time_evolution/fidelity_maximization.jl")
 
 include("algorithms/toolbox.jl")
 include("algorithms/correlators.jl")
@@ -99,6 +100,7 @@ export fixedpoint
 export absorb_weight
 export ALSTruncation, FullEnvTruncation
 export su_iter, su3site_iter, simpleupdate, SimpleUpdate
+export maximize_fidelity!
 
 export InfiniteSquareNetwork
 export InfinitePartitionFunction
diff --git a/src/algorithms/time_evolution/fidelity_maximization.jl b/src/algorithms/time_evolution/fidelity_maximization.jl
new file mode 100644
index 000000000..3ccb798a4
--- /dev/null
+++ b/src/algorithms/time_evolution/fidelity_maximization.jl
@@ -0,0 +1,90 @@
+"""
+    maximize_fidelity!(
+        pepsdst::InfinitePEPS, pepssrc::InfinitePEPS, env₀::CTMRGEnv;
+        maxiter = 5, tol = 1.0e-3, kwargs...
+    )
+
+Iteratively maximize the fidelity of `pepssrc` and `pepsdst` where the contents of `pepssrc`
+are embedded into `pepsdst`. To contract the respective networks, the virtual spaces of
+`env₀` are used and kept fixed. The CTMRG contraction algorithm is specified using the
+`kwargs...` which are passed onto the `leading_boundary` calls.
+"""
+function maximize_fidelity!(
+        pepsdst::InfinitePEPS, pepssrc::InfinitePEPS, env₀::CTMRGEnv;
+        maxiter = 5, tol = 1.0e-3, kwargs...
+    )
+    @assert size(pepsdst) == size(pepssrc) "incompatible unit cell sizes"
+
+    for (pdst, psrc) in zip(unitcell(pepsdst), unitcell(pepssrc))
+        embed!(pdst, psrc)
+    end
+
+    peps = pepssrc
+    peps′ = pepsdst
+    env, = leading_boundary(env₀, peps; kwargs...)
+    peps /= sqrt(norm(peps, env)) # normalize PEPSs to ensure that fidelity is bounded by 1
+    fid = 0.0
+    for i in 1:maxiter
+        # normalize updated PEPS and contract ⟨ψ₁|ψ₂⟩
+        env′, = leading_boundary(env, peps′; kwargs...)
+        peps′ /= sqrt(norm(peps′, env′))
+        nw = InfiniteSquareNetwork(peps′, peps)
+        envnw, = leading_boundary(env′, nw; kwargs...)
+
+        # remove peps′ from fidelity network and compute fidelity
+        ∂fid = _∂contract_site(peps′, envnw)
+        fid′ = real(_fidelity(peps, ∂fid))
+        @info @sprintf("Fidmax. iter %d:   fid = %.4e   Δfid = %.4e", i, fid′, fid′ - fid)
+        1 - fid′ ≤ tol && break
+
+        # update PEPSs
+        peps = peps′
+        peps′ = ∂fid
+        fid = fid′
+    end
+
+    return peps′
+end
+
+function _fidelity(peps::InfinitePEPS, ∂fid::InfinitePEPS)
+    return prod(eachcoordinate(peps)) do (r, c)
+        @tensor conj(peps[r, c][d; D_N D_E D_S D_W]) * ∂fid[r, c][d; D_N D_E D_S D_W]
+    end
+end
+
+function _∂contract_site(peps::InfinitePEPS, env::CTMRGEnv)
+    return InfinitePEPS(
+        map(eachcoordinate(peps)) do ind
+            _∂contract_site(ind, peps, env) * _contract_corners(ind, env) /
+                _contract_vertical_edges(ind, env) / _contract_horizontal_edges(ind, env)
+        end
+    )
+end
+function _∂contract_site(ind::Tuple{Int, Int}, peps::InfinitePEPS, env::CTMRGEnv)
+    r, c = ind
+    return _∂contract_site(
+        env.corners[NORTHWEST, _prev(r, end), _prev(c, end)],
+        env.corners[NORTHEAST, _prev(r, end), _next(c, end)],
+        env.corners[SOUTHEAST, _next(r, end), _next(c, end)],
+        env.corners[SOUTHWEST, _next(r, end), _prev(c, end)],
+        env.edges[NORTH, _prev(r, end), c], env.edges[EAST, r, _next(c, end)],
+        env.edges[SOUTH, _next(r, end), c], env.edges[WEST, r, _prev(c, end)],
+        peps[r, c],
+    )
+end
+function _∂contract_site(
+        C_northwest, C_northeast, C_southeast, C_southwest,
+        E_north::CTMRG_PEPS_EdgeTensor, E_east::CTMRG_PEPS_EdgeTensor,
+        E_south::CTMRG_PEPS_EdgeTensor, E_west::CTMRG_PEPS_EdgeTensor, ψ,
+    )
+    return @autoopt @tensor peps′[d; D_N_below D_E_below D_S_below D_W_below] :=
+        E_west[χ_WSW D_W_above D_W_below; χ_WNW] *
+        C_northwest[χ_WNW; χ_NNW] *
+        E_north[χ_NNW D_N_above D_N_below; χ_NNE] *
+        C_northeast[χ_NNE; χ_ENE] *
+        E_east[χ_ENE D_E_above D_E_below; χ_ESE] *
+        C_southeast[χ_ESE; χ_SSE] *
+        E_south[χ_SSE D_S_above D_S_below; χ_SSW] *
+        C_southwest[χ_SSW; χ_WSW] *
+        ψ[d; D_N_above D_E_above D_S_above D_W_above]
+end

From e2280726e4cbb8e97113a694f5bd0123160ce5d1 Mon Sep 17 00:00:00 2001
From: Paul Brehmer <paul.brehmer@univie.ac.at>
Date: Fri, 12 Sep 2025 11:13:10 +0200
Subject: [PATCH 02/24] Fix environment initialization for differently sized
 PEPS and add docs

---
 .../time_evolution/fidelity_maximization.jl   | 41 +++++++++++--------
 1 file changed, 24 insertions(+), 17 deletions(-)

diff --git a/src/algorithms/time_evolution/fidelity_maximization.jl b/src/algorithms/time_evolution/fidelity_maximization.jl
index 3ccb798a4..3479a868e 100644
--- a/src/algorithms/time_evolution/fidelity_maximization.jl
+++ b/src/algorithms/time_evolution/fidelity_maximization.jl
@@ -1,58 +1,64 @@
 """
     maximize_fidelity!(
-        pepsdst::InfinitePEPS, pepssrc::InfinitePEPS, env₀::CTMRGEnv;
-        maxiter = 5, tol = 1.0e-3, kwargs...
+        pepsdst::InfinitePEPS, pepssrc::InfinitePEPS, envspace;
+        maxiter = 5, tol = 1.0e-3, boundary_alg=(; verbosity=1)
     )
 
 Iteratively maximize the fidelity of `pepssrc` and `pepsdst` where the contents of `pepssrc`
-are embedded into `pepsdst`. To contract the respective networks, the virtual spaces of
-`env₀` are used and kept fixed. The CTMRG contraction algorithm is specified using the
-`kwargs...` which are passed onto the `leading_boundary` calls.
+are embedded into `pepsdst`. To contract the respective networks, the specified `envspace`
+is used on the environment bonds and kept fixed. The CTMRG contraction algorithm is specified
+via the `boundary_alg` `NamedTuple`.
 """
 function maximize_fidelity!(
-        pepsdst::InfinitePEPS, pepssrc::InfinitePEPS, env₀::CTMRGEnv;
-        maxiter = 5, tol = 1.0e-3, kwargs...
+        pepsdst::InfinitePEPS, pepssrc::InfinitePEPS, envspace;
+        maxiter = 5, tol = 1.0e-3, boundary_alg = (; verbosity = 1)
     )
     @assert size(pepsdst) == size(pepssrc) "incompatible unit cell sizes"
 
+    # absorb src PEPS tensors into dst tensors in-place
     for (pdst, psrc) in zip(unitcell(pepsdst), unitcell(pepssrc))
         embed!(pdst, psrc)
     end
 
     peps = pepssrc
     peps′ = pepsdst
-    env, = leading_boundary(env₀, peps; kwargs...)
+    env, = leading_boundary(CTMRGEnv(peps, envspace), peps; boundary_alg...)
     peps /= sqrt(norm(peps, env)) # normalize PEPSs to ensure that fidelity is bounded by 1
     fid = 0.0
     for i in 1:maxiter
         # normalize updated PEPS and contract ⟨ψ₁|ψ₂⟩
-        env′, = leading_boundary(env, peps′; kwargs...)
+        env′, = leading_boundary(CTMRGEnv(peps′, envspace), peps′; boundary_alg...)
         peps′ /= sqrt(norm(peps′, env′))
         nw = InfiniteSquareNetwork(peps′, peps)
-        envnw, = leading_boundary(env′, nw; kwargs...)
+        envnw, = leading_boundary(CTMRGEnv(nw, envspace), nw; boundary_alg...)
 
         # remove peps′ from fidelity network and compute fidelity
-        ∂fid = _∂contract_site(peps′, envnw)
-        fid′ = real(_fidelity(peps, ∂fid))
+        ∂nval = ∂network_value(peps′, envnw)
+        fid′ = abs2(network_value(peps, ∂nval))
         @info @sprintf("Fidmax. iter %d:   fid = %.4e   Δfid = %.4e", i, fid′, fid′ - fid)
-        1 - fid′ ≤ tol && break
+        abs(1 - fid′) ≤ tol && break
 
         # update PEPSs
         peps = peps′
-        peps′ = ∂fid
+        peps′ = ∂nval
         fid = fid′
     end
 
     return peps′
 end
 
-function _fidelity(peps::InfinitePEPS, ∂fid::InfinitePEPS)
+function network_value(peps::InfinitePEPS, ∂nval::InfinitePEPS)
     return prod(eachcoordinate(peps)) do (r, c)
-        @tensor conj(peps[r, c][d; D_N D_E D_S D_W]) * ∂fid[r, c][d; D_N D_E D_S D_W]
+        @tensor conj(peps[r, c][d; D_N D_E D_S D_W]) * ∂nval[r, c][d; D_N D_E D_S D_W]
     end
 end
 
-function _∂contract_site(peps::InfinitePEPS, env::CTMRGEnv)
+"""
+    ∂network_value(peps::InfinitePEPS, env::CTMRGEnv)
+
+Compute the `InfinitePEPS` resulting from removing the bra PEPS tensors in `network_value`.
+"""
+function ∂network_value(peps::InfinitePEPS, env::CTMRGEnv)
     return InfinitePEPS(
         map(eachcoordinate(peps)) do ind
             _∂contract_site(ind, peps, env) * _contract_corners(ind, env) /
@@ -60,6 +66,7 @@ function _∂contract_site(peps::InfinitePEPS, env::CTMRGEnv)
         end
     )
 end
+
 function _∂contract_site(ind::Tuple{Int, Int}, peps::InfinitePEPS, env::CTMRGEnv)
     r, c = ind
     return _∂contract_site(

From c04ce21758715e4896a2e6bdc01762357dfaaae2 Mon Sep 17 00:00:00 2001
From: Paul Brehmer <paul.brehmer@univie.ac.at>
Date: Fri, 12 Sep 2025 15:52:22 +0200
Subject: [PATCH 03/24] Add `single_site_fidelity_intiailize`

---
 src/PEPSKit.jl                                |  2 +-
 .../time_evolution/fidelity_maximization.jl   | 33 +++++++++++++++++++
 2 files changed, 34 insertions(+), 1 deletion(-)

diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl
index 1e8fd2459..d642330b3 100644
--- a/src/PEPSKit.jl
+++ b/src/PEPSKit.jl
@@ -100,7 +100,7 @@ export fixedpoint
 export absorb_weight
 export ALSTruncation, FullEnvTruncation
 export su_iter, su3site_iter, simpleupdate, SimpleUpdate
-export maximize_fidelity!
+export maximize_fidelity!, single_site_fidelity_initialize
 
 export InfiniteSquareNetwork
 export InfinitePartitionFunction
diff --git a/src/algorithms/time_evolution/fidelity_maximization.jl b/src/algorithms/time_evolution/fidelity_maximization.jl
index 3479a868e..8d7c4a2b4 100644
--- a/src/algorithms/time_evolution/fidelity_maximization.jl
+++ b/src/algorithms/time_evolution/fidelity_maximization.jl
@@ -1,3 +1,36 @@
+"""
+    single_site_fidelity_initialize(
+        peps₀::InfinitePEPS, envspace, [bondspace = _maxspace(peps₀)];
+        noise_amp = 1.0e-2, kwargs...
+    )
+
+Generate a single-site unit cell PEPS from a (possibly) multi-site `peps₀` by maximizing
+the fidelity w.r.t. `peps₀`. Here, `envspace` determines the virtual environment spaces for
+CTMRG contractions. By default, the maximal bond space of `peps₀` is used for all virtual
+legs of the single-site PEPS.
+
+The single-site PEPS is intialized with Gaussian nosie and multiplied by `noise_amp`.
+The `kwargs...` are passed onto the [`maximize_fidelity!`](@ref) call, refer to the docs
+for further details.
+"""
+function single_site_fidelity_initialize(
+        peps₀::InfinitePEPS, envspace, bondspace = _maxspace(peps₀);
+        noise_amp = 1.0e-2, kwargs...
+    )
+    @assert allequal(map(p -> space(p, 1), unitcell(peps₀))) "PEPS must have uniform physical spaces"
+
+    physspace = space(unitcell(peps₀)[1], 1)
+    peps_single = noise_amp * InfinitePEPS(randn, scalartype(peps₀), physspace, bondspace) # single-site unit cell with random noise
+
+    peps_uc = InfinitePEPS(fill(only(unitcell(peps_single)), size(peps₀))) # fill peps₀ unit cell with peps_singles
+    maximize_fidelity!(peps_uc, peps₀, envspace; kwargs...) # modifies peps_single in-place
+
+    return peps_single
+end
+
+# maximal virtual space over unit cell
+_spacemax(peps::InfinitePEPS) = argmax(dim, map(p -> argmax(dim, domain(p)), unitcell(peps)))
+
 """
     maximize_fidelity!(
         pepsdst::InfinitePEPS, pepssrc::InfinitePEPS, envspace;

From cc2d4acf90543a9504d6a6cf9aee406b271919d8 Mon Sep 17 00:00:00 2001
From: Paul Brehmer <paul.brehmer@univie.ac.at>
Date: Fri, 12 Sep 2025 16:38:44 +0200
Subject: [PATCH 04/24] Replace network_value with simple norm

---
 .../time_evolution/fidelity_maximization.jl   | 38 ++++++++++---------
 1 file changed, 21 insertions(+), 17 deletions(-)

diff --git a/src/algorithms/time_evolution/fidelity_maximization.jl b/src/algorithms/time_evolution/fidelity_maximization.jl
index 8d7c4a2b4..81d69e615 100644
--- a/src/algorithms/time_evolution/fidelity_maximization.jl
+++ b/src/algorithms/time_evolution/fidelity_maximization.jl
@@ -56,48 +56,52 @@ function maximize_fidelity!(
     peps = pepssrc
     peps′ = pepsdst
     env, = leading_boundary(CTMRGEnv(peps, envspace), peps; boundary_alg...)
-    peps /= sqrt(norm(peps, env)) # normalize PEPSs to ensure that fidelity is bounded by 1
+    peps /= sqrt(_local_norm(peps, peps, env)) # normalize PEPSs to ensure that fidelity is bounded by 1
     fid = 0.0
     for i in 1:maxiter
         # normalize updated PEPS and contract ⟨ψ₁|ψ₂⟩
         env′, = leading_boundary(CTMRGEnv(peps′, envspace), peps′; boundary_alg...)
-        peps′ /= sqrt(norm(peps′, env′))
+        peps′ /= sqrt(_local_norm(peps′, peps′, env′))
         nw = InfiniteSquareNetwork(peps′, peps)
         envnw, = leading_boundary(CTMRGEnv(nw, envspace), nw; boundary_alg...)
 
         # remove peps′ from fidelity network and compute fidelity
-        ∂nval = ∂network_value(peps′, envnw)
-        fid′ = abs2(network_value(peps, ∂nval))
+        ∂norm = _∂local_norm(peps′, envnw)
+        fid′ = abs2(_local_norm(peps, ∂norm))
         @info @sprintf("Fidmax. iter %d:   fid = %.4e   Δfid = %.4e", i, fid′, fid′ - fid)
         abs(1 - fid′) ≤ tol && break
 
         # update PEPSs
         peps = peps′
-        peps′ = ∂nval
+        peps′ = ∂norm
         fid = fid′
     end
 
     return peps′
 end
 
-function network_value(peps::InfinitePEPS, ∂nval::InfinitePEPS)
-    return prod(eachcoordinate(peps)) do (r, c)
-        @tensor conj(peps[r, c][d; D_N D_E D_S D_W]) * ∂nval[r, c][d; D_N D_E D_S D_W]
+"""
+$(SIGNATURES)
+
+Sum over `contract_local_norm` values of all unit cell entries.
+"""
+function _local_norm(ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv)
+    return sum(ind -> contract_local_norm((ind,), ket, bra, env), eachcoordinate(ket))
+end
+function _local_norm(peps::InfinitePEPS, ∂norm::InfinitePEPS)
+    return sum(eachcoordinate(peps)) do (r, c)
+        @tensor conj(peps[r, c][d; D_N D_E D_S D_W]) * ∂norm[r, c][d; D_N D_E D_S D_W]
     end
 end
 
+
 """
-    ∂network_value(peps::InfinitePEPS, env::CTMRGEnv)
+    _∂local_norm(peps::InfinitePEPS, env::CTMRGEnv)
 
-Compute the `InfinitePEPS` resulting from removing the bra PEPS tensors in `network_value`.
+Compute the `InfinitePEPS` resulting from removing the bra PEPS tensors in `_local_norm`.
 """
-function ∂network_value(peps::InfinitePEPS, env::CTMRGEnv)
-    return InfinitePEPS(
-        map(eachcoordinate(peps)) do ind
-            _∂contract_site(ind, peps, env) * _contract_corners(ind, env) /
-                _contract_vertical_edges(ind, env) / _contract_horizontal_edges(ind, env)
-        end
-    )
+function _∂local_norm(peps::InfinitePEPS, env::CTMRGEnv)
+    return InfinitePEPS(map(ind -> _∂contract_site(ind, peps, env), eachcoordinate(peps)))
 end
 
 function _∂contract_site(ind::Tuple{Int, Int}, peps::InfinitePEPS, env::CTMRGEnv)

From ce581c34baa0b834aa901f3e6a54b1ea380fe654 Mon Sep 17 00:00:00 2001
From: Paul Brehmer <paulbrehmer8@googlemail.com>
Date: Wed, 17 Sep 2025 14:56:02 +0200
Subject: [PATCH 05/24] Change filelocation and rename to MPSKit.approximate!

---
 src/PEPSKit.jl                                     |  2 +-
 .../fidelity_initialize.jl}                        | 14 +++++++-------
 2 files changed, 8 insertions(+), 8 deletions(-)
 rename src/algorithms/{time_evolution/fidelity_maximization.jl => optimization/fidelity_initialize.jl} (92%)

diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl
index d642330b3..cacb9057e 100644
--- a/src/PEPSKit.jl
+++ b/src/PEPSKit.jl
@@ -71,7 +71,6 @@ include("algorithms/truncation/bond_truncation.jl")
 include("algorithms/time_evolution/evoltools.jl")
 include("algorithms/time_evolution/simpleupdate.jl")
 include("algorithms/time_evolution/simpleupdate3site.jl")
-include("algorithms/time_evolution/fidelity_maximization.jl")
 
 include("algorithms/toolbox.jl")
 include("algorithms/correlators.jl")
@@ -80,6 +79,7 @@ include("utility/symmetrization.jl")
 
 include("algorithms/optimization/fixed_point_differentiation.jl")
 include("algorithms/optimization/peps_optimization.jl")
+include("algorithms/optimization/fidelity_initialize.jl")
 
 include("algorithms/select_algorithm.jl")
 
diff --git a/src/algorithms/time_evolution/fidelity_maximization.jl b/src/algorithms/optimization/fidelity_initialize.jl
similarity index 92%
rename from src/algorithms/time_evolution/fidelity_maximization.jl
rename to src/algorithms/optimization/fidelity_initialize.jl
index 81d69e615..237714364 100644
--- a/src/algorithms/time_evolution/fidelity_maximization.jl
+++ b/src/algorithms/optimization/fidelity_initialize.jl
@@ -32,17 +32,17 @@ end
 _spacemax(peps::InfinitePEPS) = argmax(dim, map(p -> argmax(dim, domain(p)), unitcell(peps)))
 
 """
-    maximize_fidelity!(
+    approximate!(
         pepsdst::InfinitePEPS, pepssrc::InfinitePEPS, envspace;
         maxiter = 5, tol = 1.0e-3, boundary_alg=(; verbosity=1)
     )
 
-Iteratively maximize the fidelity of `pepssrc` and `pepsdst` where the contents of `pepssrc`
-are embedded into `pepsdst`. To contract the respective networks, the specified `envspace`
-is used on the environment bonds and kept fixed. The CTMRG contraction algorithm is specified
-via the `boundary_alg` `NamedTuple`.
+Approximate `pepssrc` with `pepsdst` by iteratively maximizing their fidelity where the
+contents of `pepssrc` are embedded into `pepsdst`. To contract the respective networks, the
+specified `envspace` is used on the environment bonds and kept fixed. The CTMRG contraction
+algorithm is specified via the `boundary_alg` `NamedTuple`.
 """
-function maximize_fidelity!(
+function MPSKit.approximate!(
         pepsdst::InfinitePEPS, pepssrc::InfinitePEPS, envspace;
         maxiter = 5, tol = 1.0e-3, boundary_alg = (; verbosity = 1)
     )
@@ -50,7 +50,7 @@ function maximize_fidelity!(
 
     # absorb src PEPS tensors into dst tensors in-place
     for (pdst, psrc) in zip(unitcell(pepsdst), unitcell(pepssrc))
-        embed!(pdst, psrc)
+        absorb(pdst, psrc)
     end
 
     peps = pepssrc

From 83049b14442e1a46e5e763830e7735abf3e0d66d Mon Sep 17 00:00:00 2001
From: Paul Brehmer <paulbrehmer8@googlemail.com>
Date: Wed, 17 Sep 2025 15:18:38 +0200
Subject: [PATCH 06/24] Set up small test

---
 src/PEPSKit.jl                                     |  2 +-
 src/algorithms/optimization/fidelity_initialize.jl |  2 +-
 test/runtests.jl                                   |  3 +++
 test/toolbox/fidelity_initialize.jl                | 10 ++++++++++
 4 files changed, 15 insertions(+), 2 deletions(-)
 create mode 100644 test/toolbox/fidelity_initialize.jl

diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl
index cacb9057e..3c528f253 100644
--- a/src/PEPSKit.jl
+++ b/src/PEPSKit.jl
@@ -96,11 +96,11 @@ export correlator, correlation_length
 export leading_boundary
 export PEPSOptimize, GeomSum, ManualIter, LinSolver, EigSolver
 export fixedpoint
+export single_site_fidelity_initialize
 
 export absorb_weight
 export ALSTruncation, FullEnvTruncation
 export su_iter, su3site_iter, simpleupdate, SimpleUpdate
-export maximize_fidelity!, single_site_fidelity_initialize
 
 export InfiniteSquareNetwork
 export InfinitePartitionFunction
diff --git a/src/algorithms/optimization/fidelity_initialize.jl b/src/algorithms/optimization/fidelity_initialize.jl
index 237714364..e96fa08cc 100644
--- a/src/algorithms/optimization/fidelity_initialize.jl
+++ b/src/algorithms/optimization/fidelity_initialize.jl
@@ -23,7 +23,7 @@ function single_site_fidelity_initialize(
     peps_single = noise_amp * InfinitePEPS(randn, scalartype(peps₀), physspace, bondspace) # single-site unit cell with random noise
 
     peps_uc = InfinitePEPS(fill(only(unitcell(peps_single)), size(peps₀))) # fill peps₀ unit cell with peps_singles
-    maximize_fidelity!(peps_uc, peps₀, envspace; kwargs...) # modifies peps_single in-place
+    approximate!(peps_uc, peps₀, envspace; kwargs...) # modifies peps_single in-place
 
     return peps_single
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index 2dc9ec4c0..9234c2b15 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -69,6 +69,9 @@ end
         @time @safetestset "Density matrix from double-layer PEPO" begin
             include("toolbox/densitymatrices.jl")
         end
+        @time @safetestset "Fidelity initialization" begin
+            include("toolbox/fidelity_initialize.jl")
+        end
     end
     if GROUP == "ALL" || GROUP == "UTILITY"
         @time @safetestset "LocalOperator" begin
diff --git a/test/toolbox/fidelity_initialize.jl b/test/toolbox/fidelity_initialize.jl
new file mode 100644
index 000000000..b858da42a
--- /dev/null
+++ b/test/toolbox/fidelity_initialize.jl
@@ -0,0 +1,10 @@
+using TensorKit
+using PEPSKit
+using Test
+
+@testset "Fidelity initialization with approximate!" begin
+    pepssrc = InfinitePEPS(2, 2)
+    peps_single = single_site_fidelity_initialize(pepssrc, ℂ^6, ℂ^3)
+    @test peps_single isa InfinitePEPS
+end
+

From e8e8ce6188475958976349c42c68850319ea3786 Mon Sep 17 00:00:00 2001
From: Paul Brehmer <paulbrehmer8@googlemail.com>
Date: Thu, 18 Sep 2025 10:31:57 +0200
Subject: [PATCH 07/24] Fix formatting

---
 test/toolbox/fidelity_initialize.jl | 1 -
 1 file changed, 1 deletion(-)

diff --git a/test/toolbox/fidelity_initialize.jl b/test/toolbox/fidelity_initialize.jl
index b858da42a..bae3955ce 100644
--- a/test/toolbox/fidelity_initialize.jl
+++ b/test/toolbox/fidelity_initialize.jl
@@ -7,4 +7,3 @@ using Test
     peps_single = single_site_fidelity_initialize(pepssrc, ℂ^6, ℂ^3)
     @test peps_single isa InfinitePEPS
 end
-

From 5f621061780466de576b1a863436c1b1b91ac624 Mon Sep 17 00:00:00 2001
From: Paul Brehmer <paulbrehmer8@googlemail.com>
Date: Thu, 18 Sep 2025 16:55:20 +0200
Subject: [PATCH 08/24] Rewrite `approximate!` and add non-uniform unit cell
 test

---
 .../optimization/fidelity_initialize.jl       | 49 +++++++++++--------
 test/toolbox/fidelity_initialize.jl           | 14 +++++-
 2 files changed, 41 insertions(+), 22 deletions(-)

diff --git a/src/algorithms/optimization/fidelity_initialize.jl b/src/algorithms/optimization/fidelity_initialize.jl
index e96fa08cc..854bf9e4b 100644
--- a/src/algorithms/optimization/fidelity_initialize.jl
+++ b/src/algorithms/optimization/fidelity_initialize.jl
@@ -14,8 +14,8 @@ The `kwargs...` are passed onto the [`maximize_fidelity!`](@ref) call, refer to
 for further details.
 """
 function single_site_fidelity_initialize(
-        peps₀::InfinitePEPS, envspace, bondspace = _maxspace(peps₀);
-        noise_amp = 1.0e-2, kwargs...
+        peps₀::InfinitePEPS, envspace, bondspace = _spacemax(peps₀);
+        noise_amp = 1.0e-1, kwargs...
     )
     @assert allequal(map(p -> space(p, 1), unitcell(peps₀))) "PEPS must have uniform physical spaces"
 
@@ -23,9 +23,9 @@ function single_site_fidelity_initialize(
     peps_single = noise_amp * InfinitePEPS(randn, scalartype(peps₀), physspace, bondspace) # single-site unit cell with random noise
 
     peps_uc = InfinitePEPS(fill(only(unitcell(peps_single)), size(peps₀))) # fill peps₀ unit cell with peps_singles
-    approximate!(peps_uc, peps₀, envspace; kwargs...) # modifies peps_single in-place
+    peps_single, env_single = approximate!(peps_uc, peps₀, envspace; kwargs...) # modifies peps_single in-place
 
-    return peps_single
+    return peps_single, env_single
 end
 
 # maximal virtual space over unit cell
@@ -44,40 +44,49 @@ algorithm is specified via the `boundary_alg` `NamedTuple`.
 """
 function MPSKit.approximate!(
         pepsdst::InfinitePEPS, pepssrc::InfinitePEPS, envspace;
-        maxiter = 5, tol = 1.0e-3, boundary_alg = (; verbosity = 1)
+        maxiter = 10, tol = 1.0e-3, boundary_alg = (; verbosity = 1)
     )
     @assert size(pepsdst) == size(pepssrc) "incompatible unit cell sizes"
 
     # absorb src PEPS tensors into dst tensors in-place
     for (pdst, psrc) in zip(unitcell(pepsdst), unitcell(pepssrc))
-        absorb(pdst, psrc)
+        absorb!(pdst, psrc)
     end
 
-    peps = pepssrc
-    peps′ = pepsdst
+    # normalize reference PEPS
+    peps₀ = pepssrc # smaller bond spaces
+    env₀, = leading_boundary(CTMRGEnv(peps₀, envspace), peps₀; boundary_alg...)
+    peps₀ /= sqrt(abs(_local_norm(peps₀, peps₀, env₀))) # normalize to ensure that fidelity is bounded by 1
+
+    # normalize maximizer PEPS
+    peps = pepsdst
     env, = leading_boundary(CTMRGEnv(peps, envspace), peps; boundary_alg...)
-    peps /= sqrt(_local_norm(peps, peps, env)) # normalize PEPSs to ensure that fidelity is bounded by 1
+    peps /= sqrt(abs(_local_norm(peps, peps, env)))
+
+    nw₀ = InfiniteSquareNetwork(peps₀, peps)
+    envnw, = leading_boundary(CTMRGEnv(nw₀, envspace), nw₀; boundary_alg...)
+    peps′ = _∂local_norm(peps₀, envnw)
     fid = 0.0
     for i in 1:maxiter
-        # normalize updated PEPS and contract ⟨ψ₁|ψ₂⟩
-        env′, = leading_boundary(CTMRGEnv(peps′, envspace), peps′; boundary_alg...)
-        peps′ /= sqrt(_local_norm(peps′, peps′, env′))
-        nw = InfiniteSquareNetwork(peps′, peps)
-        envnw, = leading_boundary(CTMRGEnv(nw, envspace), nw; boundary_alg...)
-
-        # remove peps′ from fidelity network and compute fidelity
-        ∂norm = _∂local_norm(peps′, envnw)
-        fid′ = abs2(_local_norm(peps, ∂norm))
+        # compute fidelity from ∂norm
+        fid′ = abs2(_local_norm(peps, peps′))
         @info @sprintf("Fidmax. iter %d:   fid = %.4e   Δfid = %.4e", i, fid′, fid′ - fid)
         abs(1 - fid′) ≤ tol && break
 
-        # update PEPSs
+        # contract boundary of fidelity network
+        envnw, = leading_boundary(env, InfiniteSquareNetwork(peps, peps′); boundary_alg...)
+        ∂norm = _∂local_norm(peps, envnw)
+
+        # renormalize current PEPS
         peps = peps′
+        env, = leading_boundary(env, peps; boundary_alg...)
+        peps /= sqrt(abs(_local_norm(peps, peps, env)))
+
         peps′ = ∂norm
         fid = fid′
     end
 
-    return peps′
+    return peps, env
 end
 
 """
diff --git a/test/toolbox/fidelity_initialize.jl b/test/toolbox/fidelity_initialize.jl
index bae3955ce..046c93de8 100644
--- a/test/toolbox/fidelity_initialize.jl
+++ b/test/toolbox/fidelity_initialize.jl
@@ -1,9 +1,19 @@
 using TensorKit
 using PEPSKit
+using Random
 using Test
 
 @testset "Fidelity initialization with approximate!" begin
-    pepssrc = InfinitePEPS(2, 2)
-    peps_single = single_site_fidelity_initialize(pepssrc, ℂ^6, ℂ^3)
+    Random.seed!(18092025)
+    unitcell = (3, 3)
+    Pspaces = fill(2, unitcell...)
+    Nspaces = rand(2:4, unitcell...)
+    Espaces = rand(2:4, unitcell...)
+    pepssrc = InfinitePEPS(randn, ComplexF64, Pspaces, Nspaces, Espaces)
+
+    peps_single, = single_site_fidelity_initialize(
+        pepssrc, ℂ^6; maxiter = 10, noise_amp = 0.1,
+        boundary_alg = (; tol = 1.0e-6, maxiter = 20, verbosity = 1)
+    )
     @test peps_single isa InfinitePEPS
 end

From 03f3b2167d13871715c0c79b8d6df5ac53024044 Mon Sep 17 00:00:00 2001
From: Paul Brehmer <paulbrehmer8@googlemail.com>
Date: Thu, 18 Sep 2025 17:29:41 +0200
Subject: [PATCH 09/24] Fix _spacemax with reduce and add test

---
 .../optimization/fidelity_initialize.jl         |  6 ++++--
 test/toolbox/fidelity_initialize.jl             | 17 +++++++++++++++++
 2 files changed, 21 insertions(+), 2 deletions(-)

diff --git a/src/algorithms/optimization/fidelity_initialize.jl b/src/algorithms/optimization/fidelity_initialize.jl
index 854bf9e4b..963a23824 100644
--- a/src/algorithms/optimization/fidelity_initialize.jl
+++ b/src/algorithms/optimization/fidelity_initialize.jl
@@ -29,7 +29,9 @@ function single_site_fidelity_initialize(
 end
 
 # maximal virtual space over unit cell
-_spacemax(peps::InfinitePEPS) = argmax(dim, map(p -> argmax(dim, domain(p)), unitcell(peps)))
+function _spacemax(peps::InfinitePEPS)
+    return reduce(supremum, map(p -> supremum(domain(p)[1], domain(p)[2]), unitcell(peps)))
+end
 
 """
     approximate!(
@@ -105,7 +107,7 @@ end
 
 
 """
-    _∂local_norm(peps::InfinitePEPS, env::CTMRGEnv)
+$(SIGNATURES)
 
 Compute the `InfinitePEPS` resulting from removing the bra PEPS tensors in `_local_norm`.
 """
diff --git a/test/toolbox/fidelity_initialize.jl b/test/toolbox/fidelity_initialize.jl
index 046c93de8..54197299f 100644
--- a/test/toolbox/fidelity_initialize.jl
+++ b/test/toolbox/fidelity_initialize.jl
@@ -1,5 +1,6 @@
 using TensorKit
 using PEPSKit
+using PEPSKit: _spacemax
 using Random
 using Test
 
@@ -17,3 +18,19 @@ using Test
     )
     @test peps_single isa InfinitePEPS
 end
+
+_spacemax(peps::InfinitePEPS) = reduce(supremum, map(p -> (dim, domain(p)), PEPSKit.unitcell(peps)))
+
+@testset "_spacemax for non-uniform unit cell and symmetry sectors" begin
+    Pspaces = fill(Z2Space(0 => 1, 1 => 1), (2, 2))
+    Nspaces = [
+        Z2Space(0 => 2, 1 => 3) Z2Space(0 => 2, 1 => 5)
+        Z2Space(0 => 4, 1 => 3) Z2Space(0 => 2, 1 => 1)
+    ]
+    Espaces = [
+        Z2Space(0 => 6, 1 => 3) Z2Space(0 => 2, 1 => 1)
+        Z2Space(0 => 2, 1 => 3) Z2Space(0 => 2, 1 => 4)
+    ]
+    peps = InfinitePEPS(Pspaces, Nspaces, Espaces)
+    @test _spacemax(peps) == Z2Space(0 => 6, 1 => 5)
+end

From 10e6cc548043d4dc82f8b97677ccdd3c92362ba4 Mon Sep 17 00:00:00 2001
From: Paul Brehmer <paulbrehmer8@googlemail.com>
Date: Thu, 18 Sep 2025 17:31:26 +0200
Subject: [PATCH 10/24] Remove +spacemax from test

---
 test/toolbox/fidelity_initialize.jl | 2 --
 1 file changed, 2 deletions(-)

diff --git a/test/toolbox/fidelity_initialize.jl b/test/toolbox/fidelity_initialize.jl
index 54197299f..0c0724cae 100644
--- a/test/toolbox/fidelity_initialize.jl
+++ b/test/toolbox/fidelity_initialize.jl
@@ -19,8 +19,6 @@ using Test
     @test peps_single isa InfinitePEPS
 end
 
-_spacemax(peps::InfinitePEPS) = reduce(supremum, map(p -> (dim, domain(p)), PEPSKit.unitcell(peps)))
-
 @testset "_spacemax for non-uniform unit cell and symmetry sectors" begin
     Pspaces = fill(Z2Space(0 => 1, 1 => 1), (2, 2))
     Nspaces = [

From 2866c36be04a83dcd3767905557479e02571e5e3 Mon Sep 17 00:00:00 2001
From: Paul Brehmer <paulbrehmer8@googlemail.com>
Date: Fri, 19 Sep 2025 10:10:42 +0200
Subject: [PATCH 11/24] Fix docs reference

---
 src/algorithms/optimization/fidelity_initialize.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/algorithms/optimization/fidelity_initialize.jl b/src/algorithms/optimization/fidelity_initialize.jl
index 963a23824..05d542aac 100644
--- a/src/algorithms/optimization/fidelity_initialize.jl
+++ b/src/algorithms/optimization/fidelity_initialize.jl
@@ -10,7 +10,7 @@ CTMRG contractions. By default, the maximal bond space of `peps₀` is used for
 legs of the single-site PEPS.
 
 The single-site PEPS is intialized with Gaussian nosie and multiplied by `noise_amp`.
-The `kwargs...` are passed onto the [`maximize_fidelity!`](@ref) call, refer to the docs
+The `kwargs...` are passed onto the [`approximate!`](@ref) call, refer to the docs
 for further details.
 """
 function single_site_fidelity_initialize(

From 24417cf0e59d851037c12e2efa0963b1787233c8 Mon Sep 17 00:00:00 2001
From: leburgel <lander.burgelman@gmail.com>
Date: Fri, 19 Sep 2025 11:34:11 +0200
Subject: [PATCH 12/24] Bump TensorKit compat

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index cf50d3418..36559cd78 100644
--- a/Project.toml
+++ b/Project.toml
@@ -41,7 +41,7 @@ Printf = "1"
 QuadGK = "2.11.1"
 Random = "1"
 Statistics = "1"
-TensorKit = "0.14.9"
+TensorKit = "0.14.11"
 TensorOperations = "5"
 TestExtras = "0.3"
 VectorInterface = "0.4, 0.5"

From 384b900c1ac6c331145d3e8fffe8af8ab7e1f23d Mon Sep 17 00:00:00 2001
From: Paul Brehmer <paul.brehmer@univie.ac.at>
Date: Mon, 22 Sep 2025 10:53:19 +0200
Subject: [PATCH 13/24] Add `approximate` method

---
 .../optimization/fidelity_initialize.jl       | 19 +++++++++++++++----
 1 file changed, 15 insertions(+), 4 deletions(-)

diff --git a/src/algorithms/optimization/fidelity_initialize.jl b/src/algorithms/optimization/fidelity_initialize.jl
index 05d542aac..ce59dbc1c 100644
--- a/src/algorithms/optimization/fidelity_initialize.jl
+++ b/src/algorithms/optimization/fidelity_initialize.jl
@@ -33,17 +33,24 @@ function _spacemax(peps::InfinitePEPS)
     return reduce(supremum, map(p -> supremum(domain(p)[1], domain(p)[2]), unitcell(peps)))
 end
 
-"""
-    approximate!(
+@doc """
+    approximate(pepsdst::InfinitePEPS, pepssrc::InfinitePEPS, envspace; kwargs...)
+    approximate!(pepsdst::InfinitePEPS, pepssrc::InfinitePEPS, envspace; kwargs...)
         pepsdst::InfinitePEPS, pepssrc::InfinitePEPS, envspace;
         maxiter = 5, tol = 1.0e-3, boundary_alg=(; verbosity=1)
     )
 
 Approximate `pepssrc` with `pepsdst` by iteratively maximizing their fidelity where the
 contents of `pepssrc` are embedded into `pepsdst`. To contract the respective networks, the
-specified `envspace` is used on the environment bonds and kept fixed. The CTMRG contraction
-algorithm is specified via the `boundary_alg` `NamedTuple`.
+specified `envspace` is used on the environment bonds and kept fixed.
+
+## Keyword arguments
+- `maxiter=5` : Maximal number of maximization iterations
+- `tol=1.0e-3` : Absolute convergence tolerance for the infidelity
+- `boundary_alg=(; verbosity=2)` : CTMRG contraction algorithm, either specified as a `NamedTuple` or `CTMRGAlgorithm`
 """
+approximate, approximate!
+
 function MPSKit.approximate!(
         pepsdst::InfinitePEPS, pepssrc::InfinitePEPS, envspace;
         maxiter = 10, tol = 1.0e-3, boundary_alg = (; verbosity = 1)
@@ -57,6 +64,7 @@ function MPSKit.approximate!(
 
     # normalize reference PEPS
     peps₀ = pepssrc # smaller bond spaces
+    boundary_alg = _alg_or_nt(CTMRGAlgorithm, boundary_alg)
     env₀, = leading_boundary(CTMRGEnv(peps₀, envspace), peps₀; boundary_alg...)
     peps₀ /= sqrt(abs(_local_norm(peps₀, peps₀, env₀))) # normalize to ensure that fidelity is bounded by 1
 
@@ -90,6 +98,9 @@ function MPSKit.approximate!(
 
     return peps, env
 end
+function MPSKit.approximate(pepsdst::InfinitePEPS, pepssrc::InfinitePEPS, envspace; kwargs...)
+    return approximate!(deepcopy(pepsdst), pepssrc, envspace; kwargs...)
+end
 
 """
 $(SIGNATURES)

From a9fcd7326c15e2f3881dc27e213f1191124d7831 Mon Sep 17 00:00:00 2001
From: Paul Brehmer <paul.brehmer@univie.ac.at>
Date: Mon, 22 Sep 2025 12:00:51 +0200
Subject: [PATCH 14/24] Add `LoggingExtras.withlevel` logging

---
 .../optimization/fidelity_initialize.jl       | 91 ++++++++++++-------
 1 file changed, 58 insertions(+), 33 deletions(-)

diff --git a/src/algorithms/optimization/fidelity_initialize.jl b/src/algorithms/optimization/fidelity_initialize.jl
index ce59dbc1c..380b3eefd 100644
--- a/src/algorithms/optimization/fidelity_initialize.jl
+++ b/src/algorithms/optimization/fidelity_initialize.jl
@@ -53,7 +53,7 @@ approximate, approximate!
 
 function MPSKit.approximate!(
         pepsdst::InfinitePEPS, pepssrc::InfinitePEPS, envspace;
-        maxiter = 10, tol = 1.0e-3, boundary_alg = (; verbosity = 1)
+        maxiter = 10, tol = 1.0e-3, verbosity = 3, boundary_alg = (; verbosity = 1)
     )
     @assert size(pepsdst) == size(pepssrc) "incompatible unit cell sizes"
 
@@ -62,46 +62,71 @@ function MPSKit.approximate!(
         absorb!(pdst, psrc)
     end
 
-    # normalize reference PEPS
-    peps₀ = pepssrc # smaller bond spaces
-    boundary_alg = _alg_or_nt(CTMRGAlgorithm, boundary_alg)
-    env₀, = leading_boundary(CTMRGEnv(peps₀, envspace), peps₀; boundary_alg...)
-    peps₀ /= sqrt(abs(_local_norm(peps₀, peps₀, env₀))) # normalize to ensure that fidelity is bounded by 1
-
-    # normalize maximizer PEPS
-    peps = pepsdst
-    env, = leading_boundary(CTMRGEnv(peps, envspace), peps; boundary_alg...)
-    peps /= sqrt(abs(_local_norm(peps, peps, env)))
-
-    nw₀ = InfiniteSquareNetwork(peps₀, peps)
-    envnw, = leading_boundary(CTMRGEnv(nw₀, envspace), nw₀; boundary_alg...)
-    peps′ = _∂local_norm(peps₀, envnw)
-    fid = 0.0
-    for i in 1:maxiter
-        # compute fidelity from ∂norm
-        fid′ = abs2(_local_norm(peps, peps′))
-        @info @sprintf("Fidmax. iter %d:   fid = %.4e   Δfid = %.4e", i, fid′, fid′ - fid)
-        abs(1 - fid′) ≤ tol && break
-
-        # contract boundary of fidelity network
-        envnw, = leading_boundary(env, InfiniteSquareNetwork(peps, peps′); boundary_alg...)
-        ∂norm = _∂local_norm(peps, envnw)
-
-        # renormalize current PEPS
-        peps = peps′
-        env, = leading_boundary(env, peps; boundary_alg...)
+    log = MPSKit.IterLog("Approx.")
+    return LoggingExtras.withlevel(; verbosity) do
+        # normalize reference PEPS
+        peps₀ = pepssrc # smaller bond spaces
+        boundary_alg = _alg_or_nt(CTMRGAlgorithm, boundary_alg)
+        env₀, = leading_boundary(CTMRGEnv(peps₀, envspace), peps₀, boundary_alg)
+        peps₀ /= sqrt(abs(_local_norm(peps₀, peps₀, env₀))) # normalize to ensure that fidelity is bounded by 1
+
+        # normalize maximizer PEPS
+        peps = pepsdst
+        env, = leading_boundary(CTMRGEnv(peps, envspace), peps, boundary_alg)
         peps /= sqrt(abs(_local_norm(peps, peps, env)))
 
-        peps′ = ∂norm
-        fid = fid′
+        nw₀ = InfiniteSquareNetwork(peps₀, peps)
+        envnw, = leading_boundary(CTMRGEnv(nw₀, envspace), nw₀, boundary_alg)
+        peps′ = _∂local_norm(peps₀, envnw)
+        approximate_loginit!(log, one(real(scalartype(peps))), zero(real(scalartype(peps))))
+        for iter in 1:maxiter
+            # compute fidelity from ∂norm
+            fid = abs2(_local_norm(peps, peps′))
+            infid = 1 - fid
+            if abs(infid) ≤ tol
+                approximate_logfinish!(log, iter, infid, fid)
+                break
+            end
+            if iter == maxiter
+                approximate_logcancel!(log, iter, infid, fid)
+            else
+                approximate_logiter!(log, iter, infid, fid)
+            end
+
+            # contract boundary of fidelity network
+            envnw, = leading_boundary(env, InfiniteSquareNetwork(peps, peps′), boundary_alg)
+            ∂norm = _∂local_norm(peps, envnw)
+
+            # renormalize current PEPS
+            peps = peps′
+            env, = leading_boundary(env, peps, boundary_alg)
+            peps /= sqrt(abs(_local_norm(peps, peps, env)))
+
+            peps′ = ∂norm
+        end
+
+        return peps, env
     end
-
-    return peps, env
 end
 function MPSKit.approximate(pepsdst::InfinitePEPS, pepssrc::InfinitePEPS, envspace; kwargs...)
     return approximate!(deepcopy(pepsdst), pepssrc, envspace; kwargs...)
 end
 
+# custom fidelity maximization logging
+function approximate_loginit!(log, infid, fid)
+    return @infov 2 loginit!(log, infid, fid)
+end
+function approximate_logiter!(log, iter, infid, fid)
+    return @infov 3 logiter!(log, iter, infid, fid)
+end
+function approximate_logfinish!(log, iter, infid, fid)
+    return @infov 2 logfinish!(log, iter, infid, fid)
+end
+function approximate_logcancel!(log, iter, infid, fid)
+    return @warnv 1 logcancel!(log, iter, infid, fid)
+end
+
+
 """
 $(SIGNATURES)
 

From d1faf8d67228414b9f58f9ad7d219e69a2d29e9a Mon Sep 17 00:00:00 2001
From: Paul Brehmer <paul.brehmer@univie.ac.at>
Date: Mon, 22 Sep 2025 12:18:12 +0200
Subject: [PATCH 15/24] Update tests

---
 .../optimization/fidelity_initialize.jl           | 13 +++++--------
 test/toolbox/fidelity_initialize.jl               | 15 +++++++++++----
 2 files changed, 16 insertions(+), 12 deletions(-)

diff --git a/src/algorithms/optimization/fidelity_initialize.jl b/src/algorithms/optimization/fidelity_initialize.jl
index 380b3eefd..2c516b4e9 100644
--- a/src/algorithms/optimization/fidelity_initialize.jl
+++ b/src/algorithms/optimization/fidelity_initialize.jl
@@ -1,7 +1,6 @@
 """
     single_site_fidelity_initialize(
-        peps₀::InfinitePEPS, envspace, [bondspace = _maxspace(peps₀)];
-        noise_amp = 1.0e-2, kwargs...
+        peps₀::InfinitePEPS, envspace, [bondspace = _maxspace(peps₀)]; kwargs...
     )
 
 Generate a single-site unit cell PEPS from a (possibly) multi-site `peps₀` by maximizing
@@ -9,9 +8,10 @@ the fidelity w.r.t. `peps₀`. Here, `envspace` determines the virtual environme
 CTMRG contractions. By default, the maximal bond space of `peps₀` is used for all virtual
 legs of the single-site PEPS.
 
-The single-site PEPS is intialized with Gaussian nosie and multiplied by `noise_amp`.
-The `kwargs...` are passed onto the [`approximate!`](@ref) call, refer to the docs
-for further details.
+## Keyword arguments
+- `noise_amp=1.0-1` : Gaussian noise amplitude of initial single-site PEPS
+
+Additionally, all keyword arguments of [`approximate!`](@ref) can be passed.
 """
 function single_site_fidelity_initialize(
         peps₀::InfinitePEPS, envspace, bondspace = _spacemax(peps₀);
@@ -36,9 +36,6 @@ end
 @doc """
     approximate(pepsdst::InfinitePEPS, pepssrc::InfinitePEPS, envspace; kwargs...)
     approximate!(pepsdst::InfinitePEPS, pepssrc::InfinitePEPS, envspace; kwargs...)
-        pepsdst::InfinitePEPS, pepssrc::InfinitePEPS, envspace;
-        maxiter = 5, tol = 1.0e-3, boundary_alg=(; verbosity=1)
-    )
 
 Approximate `pepssrc` with `pepsdst` by iteratively maximizing their fidelity where the
 contents of `pepssrc` are embedded into `pepsdst`. To contract the respective networks, the
diff --git a/test/toolbox/fidelity_initialize.jl b/test/toolbox/fidelity_initialize.jl
index 0c0724cae..10b2c7120 100644
--- a/test/toolbox/fidelity_initialize.jl
+++ b/test/toolbox/fidelity_initialize.jl
@@ -6,17 +6,24 @@ using Test
 
 @testset "Fidelity initialization with approximate!" begin
     Random.seed!(18092025)
+
+    pepssrc_1x1 = InfinitePEPS(randn, ComplexF64, 2, 2)
+    peps_single1, = single_site_fidelity_initialize(
+        pepssrc_1x1, ℂ^6, ℂ^3; noise_amp = 0.3, tol = 1.0e-2,
+    )
+    @test peps_single1 isa InfinitePEPS
+
     unitcell = (3, 3)
     Pspaces = fill(2, unitcell...)
     Nspaces = rand(2:4, unitcell...)
     Espaces = rand(2:4, unitcell...)
-    pepssrc = InfinitePEPS(randn, ComplexF64, Pspaces, Nspaces, Espaces)
+    pepssrc_3x3 = InfinitePEPS(randn, ComplexF64, Pspaces, Nspaces, Espaces)
 
-    peps_single, = single_site_fidelity_initialize(
-        pepssrc, ℂ^6; maxiter = 10, noise_amp = 0.1,
+    peps_single2, = single_site_fidelity_initialize(
+        pepssrc_3x3, ℂ^6; maxiter = 10, noise_amp = 0.1,
         boundary_alg = (; tol = 1.0e-6, maxiter = 20, verbosity = 1)
     )
-    @test peps_single isa InfinitePEPS
+    @test peps_single2 isa InfinitePEPS
 end
 
 @testset "_spacemax for non-uniform unit cell and symmetry sectors" begin

From 23d55a4dfa400db17375b077674e4e7eaea918db Mon Sep 17 00:00:00 2001
From: Paul Brehmer <paul.brehmer@univie.ac.at>
Date: Mon, 22 Sep 2025 16:56:49 +0200
Subject: [PATCH 16/24] Actually return single-site PEPS

---
 src/algorithms/optimization/fidelity_initialize.jl | 4 ++--
 test/toolbox/fidelity_initialize.jl                | 6 ++++--
 2 files changed, 6 insertions(+), 4 deletions(-)

diff --git a/src/algorithms/optimization/fidelity_initialize.jl b/src/algorithms/optimization/fidelity_initialize.jl
index 2c516b4e9..d0059f295 100644
--- a/src/algorithms/optimization/fidelity_initialize.jl
+++ b/src/algorithms/optimization/fidelity_initialize.jl
@@ -23,9 +23,9 @@ function single_site_fidelity_initialize(
     peps_single = noise_amp * InfinitePEPS(randn, scalartype(peps₀), physspace, bondspace) # single-site unit cell with random noise
 
     peps_uc = InfinitePEPS(fill(only(unitcell(peps_single)), size(peps₀))) # fill peps₀ unit cell with peps_singles
-    peps_single, env_single = approximate!(peps_uc, peps₀, envspace; kwargs...) # modifies peps_single in-place
+    peps_single, = approximate!(peps_uc, peps₀, envspace; kwargs...) # modifies peps_single in-place
 
-    return peps_single, env_single
+    return InfinitePEPS([peps_single[1];;])
 end
 
 # maximal virtual space over unit cell
diff --git a/test/toolbox/fidelity_initialize.jl b/test/toolbox/fidelity_initialize.jl
index 10b2c7120..def304846 100644
--- a/test/toolbox/fidelity_initialize.jl
+++ b/test/toolbox/fidelity_initialize.jl
@@ -8,10 +8,11 @@ using Test
     Random.seed!(18092025)
 
     pepssrc_1x1 = InfinitePEPS(randn, ComplexF64, 2, 2)
-    peps_single1, = single_site_fidelity_initialize(
+    peps_single1 = single_site_fidelity_initialize(
         pepssrc_1x1, ℂ^6, ℂ^3; noise_amp = 0.3, tol = 1.0e-2,
     )
     @test peps_single1 isa InfinitePEPS
+    @test size(peps_single1) == (1, 1)
 
     unitcell = (3, 3)
     Pspaces = fill(2, unitcell...)
@@ -19,11 +20,12 @@ using Test
     Espaces = rand(2:4, unitcell...)
     pepssrc_3x3 = InfinitePEPS(randn, ComplexF64, Pspaces, Nspaces, Espaces)
 
-    peps_single2, = single_site_fidelity_initialize(
+    peps_single2 = single_site_fidelity_initialize(
         pepssrc_3x3, ℂ^6; maxiter = 10, noise_amp = 0.1,
         boundary_alg = (; tol = 1.0e-6, maxiter = 20, verbosity = 1)
     )
     @test peps_single2 isa InfinitePEPS
+    @test size(peps_single2) == (1, 1)
 end
 
 @testset "_spacemax for non-uniform unit cell and symmetry sectors" begin

From 2d7f26c04f1fee2bd0e7ca75b714015015863985 Mon Sep 17 00:00:00 2001
From: Paul Brehmer <paul.brehmer@univie.ac.at>
Date: Mon, 22 Sep 2025 17:28:18 +0200
Subject: [PATCH 17/24] Add missing break

---
 src/algorithms/optimization/fidelity_initialize.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/algorithms/optimization/fidelity_initialize.jl b/src/algorithms/optimization/fidelity_initialize.jl
index d0059f295..2486d7954 100644
--- a/src/algorithms/optimization/fidelity_initialize.jl
+++ b/src/algorithms/optimization/fidelity_initialize.jl
@@ -72,10 +72,10 @@ function MPSKit.approximate!(
         env, = leading_boundary(CTMRGEnv(peps, envspace), peps, boundary_alg)
         peps /= sqrt(abs(_local_norm(peps, peps, env)))
 
+        approximate_loginit!(log, one(real(scalartype(peps))), zero(real(scalartype(peps))))
         nw₀ = InfiniteSquareNetwork(peps₀, peps)
         envnw, = leading_boundary(CTMRGEnv(nw₀, envspace), nw₀, boundary_alg)
         peps′ = _∂local_norm(peps₀, envnw)
-        approximate_loginit!(log, one(real(scalartype(peps))), zero(real(scalartype(peps))))
         for iter in 1:maxiter
             # compute fidelity from ∂norm
             fid = abs2(_local_norm(peps, peps′))
@@ -86,6 +86,7 @@ function MPSKit.approximate!(
             end
             if iter == maxiter
                 approximate_logcancel!(log, iter, infid, fid)
+                break
             else
                 approximate_logiter!(log, iter, infid, fid)
             end

From 38b009e9d9466134e9b742e8b5cbef1726cf5101 Mon Sep 17 00:00:00 2001
From: Paul Brehmer <paul.brehmer@univie.ac.at>
Date: Tue, 23 Sep 2025 11:43:05 +0200
Subject: [PATCH 18/24] Remove export of `single_site_fidelity_initialize`

---
 src/PEPSKit.jl | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl
index 4a9ba92d1..ff31064c1 100644
--- a/src/PEPSKit.jl
+++ b/src/PEPSKit.jl
@@ -97,7 +97,6 @@ export correlator, correlation_length
 export leading_boundary
 export PEPSOptimize, GeomSum, ManualIter, LinSolver, EigSolver
 export fixedpoint
-export single_site_fidelity_initialize
 
 export absorb_weight
 export ALSTruncation, FullEnvTruncation

From cffbdd5f9326757fed8424f50488228c6b912be0 Mon Sep 17 00:00:00 2001
From: Paul Brehmer <paul.brehmer@univie.ac.at>
Date: Tue, 23 Sep 2025 12:00:24 +0200
Subject: [PATCH 19/24] Remove `absorb!` from `approximate!`

---
 .../optimization/fidelity_initialize.jl        | 18 +++++++++---------
 test/toolbox/fidelity_initialize.jl            |  2 +-
 2 files changed, 10 insertions(+), 10 deletions(-)

diff --git a/src/algorithms/optimization/fidelity_initialize.jl b/src/algorithms/optimization/fidelity_initialize.jl
index 2486d7954..27b7dcef0 100644
--- a/src/algorithms/optimization/fidelity_initialize.jl
+++ b/src/algorithms/optimization/fidelity_initialize.jl
@@ -9,6 +9,7 @@ CTMRG contractions. By default, the maximal bond space of `peps₀` is used for
 legs of the single-site PEPS.
 
 ## Keyword arguments
+
 - `noise_amp=1.0-1` : Gaussian noise amplitude of initial single-site PEPS
 
 Additionally, all keyword arguments of [`approximate!`](@ref) can be passed.
@@ -22,8 +23,10 @@ function single_site_fidelity_initialize(
     physspace = space(unitcell(peps₀)[1], 1)
     peps_single = noise_amp * InfinitePEPS(randn, scalartype(peps₀), physspace, bondspace) # single-site unit cell with random noise
 
+    # absorb peps₀ tensors into single-site tensors in-place
     peps_uc = InfinitePEPS(fill(only(unitcell(peps_single)), size(peps₀))) # fill peps₀ unit cell with peps_singles
-    peps_single, = approximate!(peps_uc, peps₀, envspace; kwargs...) # modifies peps_single in-place
+    absorb!(peps_uc[1], peps₀[1]) # absorb (1, 1) tensor of peps₀ (applies to all peps_uc entries since absorb! is mutating)
+    peps_single, = approximate!(peps_uc, peps₀, envspace; kwargs...)
 
     return InfinitePEPS([peps_single[1];;])
 end
@@ -37,11 +40,12 @@ end
     approximate(pepsdst::InfinitePEPS, pepssrc::InfinitePEPS, envspace; kwargs...)
     approximate!(pepsdst::InfinitePEPS, pepssrc::InfinitePEPS, envspace; kwargs...)
 
-Approximate `pepssrc` with `pepsdst` by iteratively maximizing their fidelity where the
-contents of `pepssrc` are embedded into `pepsdst`. To contract the respective networks, the
-specified `envspace` is used on the environment bonds and kept fixed.
+Approximate `pepssrc` with `pepsdst` by iteratively maximizing their fidelity, using
+`pepsdst` as an initial guess. To contract the respective networks, the specified `envspace`
+is used on the environment bonds and kept fixed.
 
 ## Keyword arguments
+
 - `maxiter=5` : Maximal number of maximization iterations
 - `tol=1.0e-3` : Absolute convergence tolerance for the infidelity
 - `boundary_alg=(; verbosity=2)` : CTMRG contraction algorithm, either specified as a `NamedTuple` or `CTMRGAlgorithm`
@@ -53,11 +57,7 @@ function MPSKit.approximate!(
         maxiter = 10, tol = 1.0e-3, verbosity = 3, boundary_alg = (; verbosity = 1)
     )
     @assert size(pepsdst) == size(pepssrc) "incompatible unit cell sizes"
-
-    # absorb src PEPS tensors into dst tensors in-place
-    for (pdst, psrc) in zip(unitcell(pepsdst), unitcell(pepssrc))
-        absorb!(pdst, psrc)
-    end
+    @assert all(map((pdst, psrc) -> space(pdst, 1) == space(psrc, 1), unitcell(pepsdst), unitcell(pepssrc))) "incompatible physical spaces"
 
     log = MPSKit.IterLog("Approx.")
     return LoggingExtras.withlevel(; verbosity) do
diff --git a/test/toolbox/fidelity_initialize.jl b/test/toolbox/fidelity_initialize.jl
index def304846..0bc759eac 100644
--- a/test/toolbox/fidelity_initialize.jl
+++ b/test/toolbox/fidelity_initialize.jl
@@ -1,6 +1,6 @@
 using TensorKit
 using PEPSKit
-using PEPSKit: _spacemax
+using PEPSKit: single_site_fidelity_initialize, _spacemax
 using Random
 using Test
 

From 7b565560ffa4cf8e0f37e45c6a639adce6022adf Mon Sep 17 00:00:00 2001
From: Paul Brehmer <paul.brehmer@univie.ac.at>
Date: Tue, 23 Sep 2025 12:05:17 +0200
Subject: [PATCH 20/24] Add comments

---
 src/algorithms/optimization/fidelity_initialize.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/algorithms/optimization/fidelity_initialize.jl b/src/algorithms/optimization/fidelity_initialize.jl
index 27b7dcef0..fddbff846 100644
--- a/src/algorithms/optimization/fidelity_initialize.jl
+++ b/src/algorithms/optimization/fidelity_initialize.jl
@@ -73,7 +73,7 @@ function MPSKit.approximate!(
         peps /= sqrt(abs(_local_norm(peps, peps, env)))
 
         approximate_loginit!(log, one(real(scalartype(peps))), zero(real(scalartype(peps))))
-        nw₀ = InfiniteSquareNetwork(peps₀, peps)
+        nw₀ = InfiniteSquareNetwork(peps₀, peps) # peps₀ has different virtual spaces than peps
         envnw, = leading_boundary(CTMRGEnv(nw₀, envspace), nw₀, boundary_alg)
         peps′ = _∂local_norm(peps₀, envnw)
         for iter in 1:maxiter
@@ -92,6 +92,7 @@ function MPSKit.approximate!(
             end
 
             # contract boundary of fidelity network
+            # initialize CTMRG on environment of peps′ (must have matching virtual spaces!)
             envnw, = leading_boundary(env, InfiniteSquareNetwork(peps, peps′), boundary_alg)
             ∂norm = _∂local_norm(peps, envnw)
 

From 9b12a60c38c9eb752835be9c69bfa001c8f0aa2e Mon Sep 17 00:00:00 2001
From: Paul Brehmer <paul.brehmer@univie.ac.at>
Date: Wed, 24 Sep 2025 17:11:35 +0200
Subject: [PATCH 21/24] Add `ApproximateAlgorithm` abstract type and
 `FidelityMaxCrude` algorithm struct

---
 src/Defaults.jl                               | 15 ++++
 src/PEPSKit.jl                                |  1 +
 .../optimization/fidelity_initialize.jl       | 88 ++++++++++++++++---
 3 files changed, 92 insertions(+), 12 deletions(-)

diff --git a/src/Defaults.jl b/src/Defaults.jl
index b1f901d85..e33f93da9 100644
--- a/src/Defaults.jl
+++ b/src/Defaults.jl
@@ -73,6 +73,14 @@ Module containing default algorithm parameter values and arguments.
 * `ls_maxfg=$(Defaults.ls_maxfg)` : Maximum number of function evaluations for the line search in each step of the optimization.
 * `lbfgs_memory=$(Defaults.lbfgs_memory)` : Size of limited memory representation of BFGS Hessian matrix.
 
+## Approximate
+
+* `approximate_alg=:$(Defaults.approximate_alg)` : Default `ApproximateAlgorithm` for PEPS approximation.
+* `approximate_tol=$(Defaults.approximate_tol)` : Infidelity tolerance of the approximation iteration.
+* `approximate_maxiter=$(Defaults.approximate_maxiter)` : Maximal number of approximation steps.
+* `approximate_miniter=$(Defaults.approximate_miniter)` : Minimal number of approximation steps.
+* `approximate_verbosity=$(Defaults.approximate_verbosity)` : Approximator output information verbosity.
+
 ## OhMyThreads scheduler
 
 - `scheduler=Ref{Scheduler}(...)` : Multithreading scheduler which can be accessed via `set_scheduler!`.
@@ -125,6 +133,13 @@ const ls_maxiter = 10
 const ls_maxfg = 20
 const lbfgs_memory = 20
 
+# Approximate
+const approximate_alg = :fidelitymaxcrude
+const approximate_tol = 1.0e-2
+const approximate_maxiter = 10
+const approximate_miniter = 3
+const approximate_verbosity = 3
+
 # OhMyThreads scheduler defaults
 const scheduler = Ref{Scheduler}()
 
diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl
index ff31064c1..29ce65754 100644
--- a/src/PEPSKit.jl
+++ b/src/PEPSKit.jl
@@ -97,6 +97,7 @@ export correlator, correlation_length
 export leading_boundary
 export PEPSOptimize, GeomSum, ManualIter, LinSolver, EigSolver
 export fixedpoint
+export FidelityMaxCrude
 
 export absorb_weight
 export ALSTruncation, FullEnvTruncation
diff --git a/src/algorithms/optimization/fidelity_initialize.jl b/src/algorithms/optimization/fidelity_initialize.jl
index fddbff846..fa0801008 100644
--- a/src/algorithms/optimization/fidelity_initialize.jl
+++ b/src/algorithms/optimization/fidelity_initialize.jl
@@ -1,3 +1,64 @@
+"""
+$(TYPEDEF)
+
+Abstract super type for tensor approximation algorithms.
+"""
+abstract type ApproximateAlgorithm end
+
+const APPROXIMATE_SYMBOLS = IdDict{Symbol, Type{<:ApproximateAlgorithm}}()
+
+"""
+$(TYPEDEF)
+
+PEPS approximation algorithm maximizing the fidelity by successively applying the fidelity
+derivative with respect to the approximator PEPS.
+
+## Constructors
+
+    FidelityMaxCrude(; kwargs...)
+
+Construct the approximation algorithm from the following keyword arguments:
+
+* `tol::Float64=$(Defaults.approximate_tol)` : Infidelity tolerance of the approximation iteration.
+* `maxiter::Int=$(Defaults.approximate_maxiter)` : Maximal number of approximation steps.
+* `miniter::Int=$(Defaults.approximate_miniter)` : Minimal number of approximation steps.
+* `verbosity::Int=$(Defaults.approximate_verbosity)` : Approximator output information verbosity.
+* `boundary_alg::Union{<:CTMRGAlgorithm,NamedTuple}` : CTMRG algorithm used for contracting the norm and fidelity networks.
+"""
+struct FidelityMaxCrude <: ApproximateAlgorithm
+    tol::Float64
+    maxiter::Int
+    miniter::Int
+    verbosity::Int
+    boundary_alg::CTMRGAlgorithm
+end
+function FidelityMaxCrude(; kwargs...)
+    return ApproximateAlgorithm(; alg = :fidelitymaxcrude, kwargs...)
+end
+
+APPROXIMATE_SYMBOLS[:fidelitymaxcrude] = FidelityMaxCrude
+
+"""
+    ApproximateAlgorithm(; kwargs...)
+
+Keyword argument parser returning the appropriate `ApproximateAlgorithm` algorithm struct.
+"""
+function ApproximateAlgorithm(;
+        alg = Defaults.approximate_alg,
+        tol = Defaults.approximate_tol,
+        maxiter = Defaults.approximate_maxiter, miniter = Defaults.approximate_miniter,
+        verbosity = Defaults.approximate_verbosity,
+        boundary_alg = (; verbosity = maximum(1, verbosity - 2)),  # shouldn't be smaller than one by default
+    )
+    # replace symbol with projector alg type
+    haskey(APPROXIMATE_SYMBOLS, alg) || throw(ArgumentError("unknown approximate algorithm: $alg"))
+    alg_type = APPROXIMATE_SYMBOLS[alg]
+
+    boundary_algorithm = _alg_or_nt(CTMRGAlgorithm, boundary_alg)
+
+    return alg_type(tol, maxiter, miniter, verbosity, boundary_algorithm)
+end
+
 """
     single_site_fidelity_initialize(
         peps₀::InfinitePEPS, envspace, [bondspace = _maxspace(peps₀)]; kwargs...
@@ -52,39 +113,42 @@ is used on the environment bonds and kept fixed.
 """
 approximate, approximate!
 
+
+function MPSKit.approximate!(pepsdst::InfinitePEPS, pepssrc::InfinitePEPS, envspace; kwargs...)
+    alg = ApproximateAlgorithm(; kwargs...)
+    return approximate!(pepsdst, pepssrc, envspace, alg)
+end
 function MPSKit.approximate!(
-        pepsdst::InfinitePEPS, pepssrc::InfinitePEPS, envspace;
-        maxiter = 10, tol = 1.0e-3, verbosity = 3, boundary_alg = (; verbosity = 1)
+        pepsdst::InfinitePEPS, pepssrc::InfinitePEPS, envspace, alg::ApproximateAlgorithm
     )
     @assert size(pepsdst) == size(pepssrc) "incompatible unit cell sizes"
     @assert all(map((pdst, psrc) -> space(pdst, 1) == space(psrc, 1), unitcell(pepsdst), unitcell(pepssrc))) "incompatible physical spaces"
 
     log = MPSKit.IterLog("Approx.")
-    return LoggingExtras.withlevel(; verbosity) do
+    return LoggingExtras.withlevel(; alg.verbosity) do
         # normalize reference PEPS
         peps₀ = pepssrc # smaller bond spaces
-        boundary_alg = _alg_or_nt(CTMRGAlgorithm, boundary_alg)
-        env₀, = leading_boundary(CTMRGEnv(peps₀, envspace), peps₀, boundary_alg)
+        env₀, = leading_boundary(CTMRGEnv(peps₀, envspace), peps₀, alg.boundary_alg)
         peps₀ /= sqrt(abs(_local_norm(peps₀, peps₀, env₀))) # normalize to ensure that fidelity is bounded by 1
 
         # normalize maximizer PEPS
         peps = pepsdst
-        env, = leading_boundary(CTMRGEnv(peps, envspace), peps, boundary_alg)
+        env, = leading_boundary(CTMRGEnv(peps, envspace), peps, alg.boundary_alg)
         peps /= sqrt(abs(_local_norm(peps, peps, env)))
 
         approximate_loginit!(log, one(real(scalartype(peps))), zero(real(scalartype(peps))))
         nw₀ = InfiniteSquareNetwork(peps₀, peps) # peps₀ has different virtual spaces than peps
-        envnw, = leading_boundary(CTMRGEnv(nw₀, envspace), nw₀, boundary_alg)
+        envnw, = leading_boundary(CTMRGEnv(nw₀, envspace), nw₀, alg.boundary_alg)
         peps′ = _∂local_norm(peps₀, envnw)
         for iter in 1:maxiter
             # compute fidelity from ∂norm
             fid = abs2(_local_norm(peps, peps′))
             infid = 1 - fid
-            if abs(infid) ≤ tol
+            if abs(infid) ≤ alg.tol && iter ≥ alg.miniter
                 approximate_logfinish!(log, iter, infid, fid)
                 break
             end
-            if iter == maxiter
+            if iter == alg.maxiter
                 approximate_logcancel!(log, iter, infid, fid)
                 break
             else
@@ -93,7 +157,7 @@ function MPSKit.approximate!(
 
             # contract boundary of fidelity network
             # initialize CTMRG on environment of peps′ (must have matching virtual spaces!)
-            envnw, = leading_boundary(env, InfiniteSquareNetwork(peps, peps′), boundary_alg)
+            envnw, = leading_boundary(env, InfiniteSquareNetwork(peps, peps′), alg.boundary_alg)
             ∂norm = _∂local_norm(peps, envnw)
 
             # renormalize current PEPS
@@ -107,8 +171,8 @@ function MPSKit.approximate!(
         return peps, env
     end
 end
-function MPSKit.approximate(pepsdst::InfinitePEPS, pepssrc::InfinitePEPS, envspace; kwargs...)
-    return approximate!(deepcopy(pepsdst), pepssrc, envspace; kwargs...)
+function MPSKit.approximate(pepsdst::InfinitePEPS, args...; kwargs...)
+    return approximate!(deepcopy(pepsdst), args...; kwargs...)
 end
 
 # custom fidelity maximization logging

From fc93f50a778d7361733ff427d068ec0cf0ba3939 Mon Sep 17 00:00:00 2001
From: Paul Brehmer <paul.brehmer@univie.ac.at>
Date: Wed, 24 Sep 2025 17:52:29 +0200
Subject: [PATCH 22/24] Match `approximate!` function signature with MPSKit

---
 .../optimization/fidelity_initialize.jl       | 86 +++++++++++--------
 1 file changed, 48 insertions(+), 38 deletions(-)

diff --git a/src/algorithms/optimization/fidelity_initialize.jl b/src/algorithms/optimization/fidelity_initialize.jl
index fa0801008..17bed867f 100644
--- a/src/algorithms/optimization/fidelity_initialize.jl
+++ b/src/algorithms/optimization/fidelity_initialize.jl
@@ -30,6 +30,7 @@ struct FidelityMaxCrude <: ApproximateAlgorithm
     maxiter::Int
     miniter::Int
     verbosity::Int
+    envspace::ElementarySpace
     boundary_alg::CTMRGAlgorithm
 end
 function FidelityMaxCrude(; kwargs...)
@@ -61,33 +62,32 @@ end
 
 """
     single_site_fidelity_initialize(
-        peps₀::InfinitePEPS, envspace, [bondspace = _maxspace(peps₀)]; kwargs...
+        peps::InfinitePEPS, [bondspace = _maxspace(peps₀)]; kwargs...
     )
 
-Generate a single-site unit cell PEPS from a (possibly) multi-site `peps₀` by maximizing
-the fidelity w.r.t. `peps₀`. Here, `envspace` determines the virtual environment spaces for
-CTMRG contractions. By default, the maximal bond space of `peps₀` is used for all virtual
-legs of the single-site PEPS.
+Generate a single-site unit cell PEPS from a (possibly) multi-site `peps` by approximating
+the respective entries using [`approximate!`](@ref). By default, the maximal bond space of
+`peps₀` is used for all virtual legs of the single-site PEPS.
 
 ## Keyword arguments
 
-- `noise_amp=1.0-1` : Gaussian noise amplitude of initial single-site PEPS
+* `noise_amp=1.0-1` : Gaussian noise amplitude of initial single-site PEPS
 
-Additionally, all keyword arguments of [`approximate!`](@ref) can be passed.
+All additional keyword arguments will be passed to the [`approximate!`](@ref) call.
 """
 function single_site_fidelity_initialize(
-        peps₀::InfinitePEPS, envspace, bondspace = _spacemax(peps₀);
+        peps::InfinitePEPS, bondspace = _spacemax(peps);
         noise_amp = 1.0e-1, kwargs...
     )
-    @assert allequal(map(p -> space(p, 1), unitcell(peps₀))) "PEPS must have uniform physical spaces"
+    @assert allequal(map(p -> space(p, 1), unitcell(peps))) "PEPS must have uniform physical spaces"
 
-    physspace = space(unitcell(peps₀)[1], 1)
-    peps_single = noise_amp * InfinitePEPS(randn, scalartype(peps₀), physspace, bondspace) # single-site unit cell with random noise
+    physspace = space(unitcell(peps)[1], 1)
+    peps_single = noise_amp * InfinitePEPS(randn, scalartype(peps), physspace, bondspace) # single-site unit cell with random noise
 
     # absorb peps₀ tensors into single-site tensors in-place
-    peps_uc = InfinitePEPS(fill(only(unitcell(peps_single)), size(peps₀))) # fill peps₀ unit cell with peps_singles
-    absorb!(peps_uc[1], peps₀[1]) # absorb (1, 1) tensor of peps₀ (applies to all peps_uc entries since absorb! is mutating)
-    peps_single, = approximate!(peps_uc, peps₀, envspace; kwargs...)
+    peps_uc = InfinitePEPS(fill(only(unitcell(peps_single)), size(peps))) # fill unit cell with peps_single tensors
+    absorb!(peps_uc[1], peps[1]) # absorb (1, 1) tensor of peps₀ (applies to all peps_uc entries since absorb! is mutating)
+    peps_single, = approximate!(peps_uc, peps; kwargs...)
 
     return InfinitePEPS([peps_single[1];;])
 end
@@ -98,48 +98,57 @@ function _spacemax(peps::InfinitePEPS)
 end
 
 @doc """
-    approximate(pepsdst::InfinitePEPS, pepssrc::InfinitePEPS, envspace; kwargs...)
-    approximate!(pepsdst::InfinitePEPS, pepssrc::InfinitePEPS, envspace; kwargs...)
+    approximate(ψ₀::InfinitePEPS, ψ::InfinitePEPS; kwargs...)
+    approximate!(ψ₀::InfinitePEPS, ψ::InfinitePEPS; kwargs...)
+    # expert versions
+    approximate(ψ₀::InfinitePEPS, ψ::InfinitePEPS, alg::ApproximateAlgorithm)
+    approximate!(ψ₀::InfinitePEPS, ψ::InfinitePEPS, alg::ApproximateAlgorithm)
 
-Approximate `pepssrc` with `pepsdst` by iteratively maximizing their fidelity, using
-`pepsdst` as an initial guess. To contract the respective networks, the specified `envspace`
-is used on the environment bonds and kept fixed.
+Approximate `ψ` from the initial guess `ψ₀`. The approximation algorithm is specified via
+the keyword arguments or directly by passing an [`ApproximateAlgorithm`](@ref) struct.
 
 ## Keyword arguments
 
-- `maxiter=5` : Maximal number of maximization iterations
-- `tol=1.0e-3` : Absolute convergence tolerance for the infidelity
-- `boundary_alg=(; verbosity=2)` : CTMRG contraction algorithm, either specified as a `NamedTuple` or `CTMRGAlgorithm`
+* `alg::Symbol=:$(Defaults.approximate_alg)` : Approximation algorithm, which can be one of the following:
+    - `:fidelitymaxcrude` : Maximize the fidelity from the fidelity gradient in a power-method fashion.
+* `tol::Float64=$(Defaults.approximate_tol)` : Infidelity tolerance of the approximation iteration.
+* `maxiter::Int=$(Defaults.approximate_maxiter)` : Maximal number of approximation steps.
+* `miniter::Int=$(Defaults.approximate_miniter)` : Minimal number of approximation steps.
+* `verbosity::Int=$(Defaults.approximate_verbosity)` : Approximator output information verbosity.
+* `boundary_alg::Union{<:CTMRGAlgorithm,NamedTuple}` : CTMRG algorithm used for contracting the norm and fidelity networks.
+
+## Return values
+
+The final approximator and its environment are returned.
 """
 approximate, approximate!
 
 
-function MPSKit.approximate!(pepsdst::InfinitePEPS, pepssrc::InfinitePEPS, envspace; kwargs...)
+function MPSKit.approximate!(ψ₀::InfinitePEPS, ψ::InfinitePEPS; kwargs...)
     alg = ApproximateAlgorithm(; kwargs...)
-    return approximate!(pepsdst, pepssrc, envspace, alg)
+    return approximate!(ψ₀, ψ, alg)
 end
-function MPSKit.approximate!(
-        pepsdst::InfinitePEPS, pepssrc::InfinitePEPS, envspace, alg::ApproximateAlgorithm
-    )
-    @assert size(pepsdst) == size(pepssrc) "incompatible unit cell sizes"
-    @assert all(map((pdst, psrc) -> space(pdst, 1) == space(psrc, 1), unitcell(pepsdst), unitcell(pepssrc))) "incompatible physical spaces"
+function MPSKit.approximate!(ψ₀::InfinitePEPS, ψ::InfinitePEPS, alg::FidelityMaxCrude)
+    @assert size(ψ₀) == size(ψ) "incompatible unit cell sizes"
+    @assert all(map((p₀, p) -> space(p₀, 1) == space(p, 1), unitcell(ψ₀), unitcell(ψ))) "incompatible physical spaces"
 
-    log = MPSKit.IterLog("Approx.")
+    log = MPSKit.IterLog("Approximate")
     return LoggingExtras.withlevel(; alg.verbosity) do
         # normalize reference PEPS
-        peps₀ = pepssrc # smaller bond spaces
-        env₀, = leading_boundary(CTMRGEnv(peps₀, envspace), peps₀, alg.boundary_alg)
-        peps₀ /= sqrt(abs(_local_norm(peps₀, peps₀, env₀))) # normalize to ensure that fidelity is bounded by 1
+        peps_init = ψ # smaller bond spaces
+        envspace = domain(ψ₀[1])
+        env₀, = leading_boundary(CTMRGEnv(peps_init, envspace), peps_init, alg.boundary_alg)
+        peps_init /= sqrt(abs(_local_norm(peps_init, peps_init, env₀))) # normalize to ensure that fidelity is bounded by 1
 
         # normalize maximizer PEPS
-        peps = pepsdst
+        peps = ψ₀
         env, = leading_boundary(CTMRGEnv(peps, envspace), peps, alg.boundary_alg)
         peps /= sqrt(abs(_local_norm(peps, peps, env)))
 
         approximate_loginit!(log, one(real(scalartype(peps))), zero(real(scalartype(peps))))
-        nw₀ = InfiniteSquareNetwork(peps₀, peps) # peps₀ has different virtual spaces than peps
+        nw₀ = InfiniteSquareNetwork(peps_init, peps) # peps₀ has different virtual spaces than peps
         envnw, = leading_boundary(CTMRGEnv(nw₀, envspace), nw₀, alg.boundary_alg)
-        peps′ = _∂local_norm(peps₀, envnw)
+        peps′ = _∂local_norm(peps_init, envnw)
         for iter in 1:maxiter
             # compute fidelity from ∂norm
             fid = abs2(_local_norm(peps, peps′))
@@ -171,8 +180,8 @@ function MPSKit.approximate!(
         return peps, env
     end
 end
-function MPSKit.approximate(pepsdst::InfinitePEPS, args...; kwargs...)
-    return approximate!(deepcopy(pepsdst), args...; kwargs...)
+function MPSKit.approximate(ψ₀::InfinitePEPS, args...; kwargs...)
+    return approximate!(deepcopy(ψ₀), args...; kwargs...)
 end
 
 # custom fidelity maximization logging
@@ -214,6 +223,7 @@ function _∂local_norm(peps::InfinitePEPS, env::CTMRGEnv)
     return InfinitePEPS(map(ind -> _∂contract_site(ind, peps, env), eachcoordinate(peps)))
 end
 
+# contract CTMRG environment leaving open the bra-PEPS virtual and physical bonds
 function _∂contract_site(ind::Tuple{Int, Int}, peps::InfinitePEPS, env::CTMRGEnv)
     r, c = ind
     return _∂contract_site(

From 4bef8fb9d54e4e41115ede76dbb7240290b4111b Mon Sep 17 00:00:00 2001
From: Paul Brehmer <paul.brehmer@univie.ac.at>
Date: Wed, 24 Sep 2025 18:48:22 +0200
Subject: [PATCH 23/24] Remove `approximate!` and fix errors

---
 src/PEPSKit.jl                                |  2 +-
 .../optimization/fidelity_initialize.jl       | 26 +++++++------------
 test/toolbox/fidelity_initialize.jl           | 11 ++++----
 3 files changed, 17 insertions(+), 22 deletions(-)

diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl
index 29ce65754..7da8c0e96 100644
--- a/src/PEPSKit.jl
+++ b/src/PEPSKit.jl
@@ -97,7 +97,7 @@ export correlator, correlation_length
 export leading_boundary
 export PEPSOptimize, GeomSum, ManualIter, LinSolver, EigSolver
 export fixedpoint
-export FidelityMaxCrude
+export FidelityMaxCrude, approximate
 
 export absorb_weight
 export ALSTruncation, FullEnvTruncation
diff --git a/src/algorithms/optimization/fidelity_initialize.jl b/src/algorithms/optimization/fidelity_initialize.jl
index 17bed867f..40301bce2 100644
--- a/src/algorithms/optimization/fidelity_initialize.jl
+++ b/src/algorithms/optimization/fidelity_initialize.jl
@@ -30,7 +30,6 @@ struct FidelityMaxCrude <: ApproximateAlgorithm
     maxiter::Int
     miniter::Int
     verbosity::Int
-    envspace::ElementarySpace
     boundary_alg::CTMRGAlgorithm
 end
 function FidelityMaxCrude(; kwargs...)
@@ -49,7 +48,7 @@ function ApproximateAlgorithm(;
         tol = Defaults.approximate_tol,
         maxiter = Defaults.approximate_maxiter, miniter = Defaults.approximate_miniter,
         verbosity = Defaults.approximate_verbosity,
-        boundary_alg = (; verbosity = maximum(1, verbosity - 2)),  # shouldn't be smaller than one by default
+        boundary_alg = (; verbosity = max(1, verbosity - 2)),  # shouldn't be smaller than one by default
     )
     # replace symbol with projector alg type
     haskey(APPROXIMATE_SYMBOLS, alg) || throw(ArgumentError("unknown approximate algorithm: $alg"))
@@ -87,7 +86,7 @@ function single_site_fidelity_initialize(
     # absorb peps₀ tensors into single-site tensors in-place
     peps_uc = InfinitePEPS(fill(only(unitcell(peps_single)), size(peps))) # fill unit cell with peps_single tensors
     absorb!(peps_uc[1], peps[1]) # absorb (1, 1) tensor of peps₀ (applies to all peps_uc entries since absorb! is mutating)
-    peps_single, = approximate!(peps_uc, peps; kwargs...)
+    peps_single, = approximate(peps_uc, peps; kwargs...)
 
     return InfinitePEPS([peps_single[1];;])
 end
@@ -99,10 +98,8 @@ end
 
 @doc """
     approximate(ψ₀::InfinitePEPS, ψ::InfinitePEPS; kwargs...)
-    approximate!(ψ₀::InfinitePEPS, ψ::InfinitePEPS; kwargs...)
-    # expert versions
+    # expert version
     approximate(ψ₀::InfinitePEPS, ψ::InfinitePEPS, alg::ApproximateAlgorithm)
-    approximate!(ψ₀::InfinitePEPS, ψ::InfinitePEPS, alg::ApproximateAlgorithm)
 
 Approximate `ψ` from the initial guess `ψ₀`. The approximation algorithm is specified via
 the keyword arguments or directly by passing an [`ApproximateAlgorithm`](@ref) struct.
@@ -121,14 +118,14 @@ the keyword arguments or directly by passing an [`ApproximateAlgorithm`](@ref) s
 
 The final approximator and its environment are returned.
 """
-approximate, approximate!
+approximate
 
 
-function MPSKit.approximate!(ψ₀::InfinitePEPS, ψ::InfinitePEPS; kwargs...)
+function MPSKit.approximate(ψ₀::InfinitePEPS, ψ::InfinitePEPS; kwargs...)
     alg = ApproximateAlgorithm(; kwargs...)
-    return approximate!(ψ₀, ψ, alg)
+    return approximate(ψ₀, ψ, alg)
 end
-function MPSKit.approximate!(ψ₀::InfinitePEPS, ψ::InfinitePEPS, alg::FidelityMaxCrude)
+function MPSKit.approximate(ψ₀::InfinitePEPS, ψ::InfinitePEPS, alg::FidelityMaxCrude)
     @assert size(ψ₀) == size(ψ) "incompatible unit cell sizes"
     @assert all(map((p₀, p) -> space(p₀, 1) == space(p, 1), unitcell(ψ₀), unitcell(ψ))) "incompatible physical spaces"
 
@@ -136,7 +133,7 @@ function MPSKit.approximate!(ψ₀::InfinitePEPS, ψ::InfinitePEPS, alg::Fidelit
     return LoggingExtras.withlevel(; alg.verbosity) do
         # normalize reference PEPS
         peps_init = ψ # smaller bond spaces
-        envspace = domain(ψ₀[1])
+        envspace = domain(ψ₀[1])[1]
         env₀, = leading_boundary(CTMRGEnv(peps_init, envspace), peps_init, alg.boundary_alg)
         peps_init /= sqrt(abs(_local_norm(peps_init, peps_init, env₀))) # normalize to ensure that fidelity is bounded by 1
 
@@ -149,7 +146,7 @@ function MPSKit.approximate!(ψ₀::InfinitePEPS, ψ::InfinitePEPS, alg::Fidelit
         nw₀ = InfiniteSquareNetwork(peps_init, peps) # peps₀ has different virtual spaces than peps
         envnw, = leading_boundary(CTMRGEnv(nw₀, envspace), nw₀, alg.boundary_alg)
         peps′ = _∂local_norm(peps_init, envnw)
-        for iter in 1:maxiter
+        for iter in 1:(alg.maxiter)
             # compute fidelity from ∂norm
             fid = abs2(_local_norm(peps, peps′))
             infid = 1 - fid
@@ -171,7 +168,7 @@ function MPSKit.approximate!(ψ₀::InfinitePEPS, ψ::InfinitePEPS, alg::Fidelit
 
             # renormalize current PEPS
             peps = peps′
-            env, = leading_boundary(env, peps, boundary_alg)
+            env, = leading_boundary(env, peps, alg.boundary_alg)
             peps /= sqrt(abs(_local_norm(peps, peps, env)))
 
             peps′ = ∂norm
@@ -180,9 +177,6 @@ function MPSKit.approximate!(ψ₀::InfinitePEPS, ψ::InfinitePEPS, alg::Fidelit
         return peps, env
     end
 end
-function MPSKit.approximate(ψ₀::InfinitePEPS, args...; kwargs...)
-    return approximate!(deepcopy(ψ₀), args...; kwargs...)
-end
 
 # custom fidelity maximization logging
 function approximate_loginit!(log, infid, fid)
diff --git a/test/toolbox/fidelity_initialize.jl b/test/toolbox/fidelity_initialize.jl
index 0bc759eac..6596881a0 100644
--- a/test/toolbox/fidelity_initialize.jl
+++ b/test/toolbox/fidelity_initialize.jl
@@ -4,12 +4,13 @@ using PEPSKit: single_site_fidelity_initialize, _spacemax
 using Random
 using Test
 
-@testset "Fidelity initialization with approximate!" begin
+@testset "Fidelity initialization with single_site_fidelity_initialize" begin
     Random.seed!(18092025)
 
-    pepssrc_1x1 = InfinitePEPS(randn, ComplexF64, 2, 2)
+    peps_1x1 = InfinitePEPS(randn, ComplexF64, 2, 2)
     peps_single1 = single_site_fidelity_initialize(
-        pepssrc_1x1, ℂ^6, ℂ^3; noise_amp = 0.3, tol = 1.0e-2,
+        peps_1x1, ℂ^3; noise_amp = 0.3, tol = 1.0e-2, miniter = 5,
+        boundary_alg = (; tol = 1.0e-6, trscheme = truncdim(10), verbosity = 1)
     )
     @test peps_single1 isa InfinitePEPS
     @test size(peps_single1) == (1, 1)
@@ -18,10 +19,10 @@ using Test
     Pspaces = fill(2, unitcell...)
     Nspaces = rand(2:4, unitcell...)
     Espaces = rand(2:4, unitcell...)
-    pepssrc_3x3 = InfinitePEPS(randn, ComplexF64, Pspaces, Nspaces, Espaces)
+    peps_3x3 = InfinitePEPS(randn, ComplexF64, Pspaces, Nspaces, Espaces)
 
     peps_single2 = single_site_fidelity_initialize(
-        pepssrc_3x3, ℂ^6; maxiter = 10, noise_amp = 0.1,
+        peps_3x3; maxiter = 10, noise_amp = 0.1,
         boundary_alg = (; tol = 1.0e-6, maxiter = 20, verbosity = 1)
     )
     @test peps_single2 isa InfinitePEPS

From 61ba01f8f436b51448bc4f6869610542fdc64d8c Mon Sep 17 00:00:00 2001
From: Paul Brehmer <paul.brehmer@univie.ac.at>
Date: Wed, 24 Sep 2025 18:53:29 +0200
Subject: [PATCH 24/24] Add small `approximate` test

---
 test/toolbox/fidelity_initialize.jl | 19 +++++++++++++++++++
 1 file changed, 19 insertions(+)

diff --git a/test/toolbox/fidelity_initialize.jl b/test/toolbox/fidelity_initialize.jl
index 6596881a0..05030fb16 100644
--- a/test/toolbox/fidelity_initialize.jl
+++ b/test/toolbox/fidelity_initialize.jl
@@ -29,6 +29,25 @@ using Test
     @test size(peps_single2) == (1, 1)
 end
 
+@testset "approximate methods with FidelityMaxCrude" begin
+    peps1 = InfinitePEPS(2, 2)
+    peps2 = InfinitePEPS(2, 2)
+    Random.seed!(1234)
+    peps_approx1, = approximate(
+        peps1, peps2; maxiter = 3,
+        boundary_alg = (; trscheme = truncdim(6), verbosity = 1)
+    )
+    Random.seed!(1234)
+    peps_approx2, = approximate(
+        peps1, peps2, FidelityMaxCrude(;
+            maxiter = 3,
+            boundary_alg = (; trscheme = truncdim(6), verbosity = 1)
+        )
+    )
+
+    @test peps_approx1 ≈ peps_approx2
+end
+
 @testset "_spacemax for non-uniform unit cell and symmetry sectors" begin
     Pspaces = fill(Z2Space(0 => 1, 1 => 1), (2, 2))
     Nspaces = [