diff --git a/Project.toml b/Project.toml
index bf587bfee..2dec89104 100644
--- a/Project.toml
+++ b/Project.toml
@@ -28,22 +28,24 @@ Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 [compat]
 Accessors = "0.1"
 ChainRulesCore = "1.0"
+ChainRulesTestUtils = "1.13"
+SafeTestsets = "0.1"
 Compat = "3.46, 4.2"
 DocStringExtensions = "0.9.3"
 FiniteDifferences = "0.12"
 KrylovKit = "0.9.5, 0.10"
 LinearAlgebra = "1"
 LoggingExtras = "1"
-MPSKit = "0.13.7"
+MPSKit = "0.13.9"
 MPSKitModels = "0.4"
-MatrixAlgebraKit = "0.5.0"
+MatrixAlgebraKit = "0.6"
 OhMyThreads = "0.7, 0.8"
 OptimKit = "0.4"
 Printf = "1"
 QuadGK = "2.11.1"
 Random = "1"
 Statistics = "1"
-TensorKit = "0.15"
+TensorKit = "0.16.2"
 TensorOperations = "5"
 TestExtras = "0.3"
 VectorInterface = "0.4, 0.5"
diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl
index c8ee35996..9c31b7476 100644
--- a/src/PEPSKit.jl
+++ b/src/PEPSKit.jl
@@ -8,9 +8,10 @@ using VectorInterface
 import VectorInterface as VI
 
 using MatrixAlgebraKit
-using MatrixAlgebraKit: TruncationStrategy, LAPACK_DivideAndConquer, LAPACK_QRIteration
-using MatrixAlgebraKit: NoTruncation, LAPACK_EighAlgorithm, truncate
-using MatrixAlgebraKit: eigh_pullback!, eigh_trunc_pullback!, findtruncated, diagview
+using MatrixAlgebraKit: LAPACK_DivideAndConquer, LAPACK_QRIteration
+using MatrixAlgebraKit:
+    TruncationStrategy, NoTruncation, truncate, findtruncated, truncation_error, diagview
+using MatrixAlgebraKit: LAPACK_EighAlgorithm, eigh_pullback!, eigh_trunc_pullback!
 
 using TensorKit
 using TensorKit: AdjointTensorMap, SectorDict
diff --git a/src/algorithms/bp/beliefpropagation.jl b/src/algorithms/bp/beliefpropagation.jl
index 9178fc16f..0246d7e20 100644
--- a/src/algorithms/bp/beliefpropagation.jl
+++ b/src/algorithms/bp/beliefpropagation.jl
@@ -126,11 +126,7 @@ function tr_distance(A::BPEnv, B::BPEnv)
     end / length(A.messages)
 end
 
-function trnorm(M::AbstractTensorMap, p::Real = 1)
-    return TensorKit._norm(svdvals(M), p, zero(real(scalartype(M))))
-end
-function trnorm!(M::AbstractTensorMap, p::Real = 1)
-    return TensorKit._norm(svdvals!(M), p, zero(real(scalartype(M))))
-end
+trnorm(M::AbstractTensorMap, p::Real = 1) = norm(svdvals(M), p)
+trnorm!(M::AbstractTensorMap, p::Real = 1) = norm(svdvals!(M), p)
 
 project_hermitian!!(t) = add(t, t', 1 / 2, 1 / 2)
diff --git a/src/algorithms/bp/gaugefix.jl b/src/algorithms/bp/gaugefix.jl
index 2e6ce4982..561dd2183 100644
--- a/src/algorithms/bp/gaugefix.jl
+++ b/src/algorithms/bp/gaugefix.jl
@@ -117,7 +117,7 @@ function SUWeight(env::BPEnv)
     wts = map(Iterators.product(1:2, axes(env, 2), axes(env, 3))) do (dir′, row, col)
         I = CartesianIndex(mod1(dir′ + 1, 2), row, col)
         sqrtM12, _, sqrtM21, _ = _sqrt_bp_messages(I, env)
-        Λ = svd_vals!(sqrtM12 * sqrtM21)
+        Λ = DiagonalTensorMap(svd_vals!(sqrtM12 * sqrtM21))
         return isdual(space(sqrtM12, 1)) ? _fliptwist_s(Λ) : Λ
     end
     return SUWeight(wts)
diff --git a/src/algorithms/contractions/bondenv/gaugefix.jl b/src/algorithms/contractions/bondenv/gaugefix.jl
index c49030c13..d30c1b2ce 100644
--- a/src/algorithms/contractions/bondenv/gaugefix.jl
+++ b/src/algorithms/contractions/bondenv/gaugefix.jl
@@ -61,8 +61,8 @@ function fixgauge_benv(
                             ↓
                             2
     =#
-    QL, L = left_orth(permute(Z, ((2, 1), (3,))))
-    QR, R = left_orth(permute(Z, ((3, 1), (2,))))
+    QL, L = left_orth(permute(Z, ((2, 1), (3,))); positive = true)
+    QR, R = left_orth(permute(Z, ((3, 1), (2,))); positive = true)
     @debug "cond(L) = $(LinearAlgebra.cond(L)); cond(R) = $(LinearAlgebra.cond(R))"
     # TODO: find a better way to fix gauge that avoids `inv`
     Linv, Rinv = inv(L), inv(R)
diff --git a/src/algorithms/contractions/ctmrg/renormalize_edge.jl b/src/algorithms/contractions/ctmrg/renormalize_edge.jl
index 8e005285e..3a4921787 100644
--- a/src/algorithms/contractions/ctmrg/renormalize_edge.jl
+++ b/src/algorithms/contractions/ctmrg/renormalize_edge.jl
@@ -91,7 +91,7 @@ function renormalize_east_edge(
     )
     return renormalize_east_edge(
         env.edges[EAST, row, _next(col, end)],
-        P_left[EAST, row, col, end],
+        P_left[EAST, row, col],
         P_right[EAST, _prev(row, end), col],
         network[row, col],
     )
diff --git a/src/algorithms/ctmrg/ctmrg.jl b/src/algorithms/ctmrg/ctmrg.jl
index 3006dc7b0..5709c0025 100644
--- a/src/algorithms/ctmrg/ctmrg.jl
+++ b/src/algorithms/ctmrg/ctmrg.jl
@@ -164,27 +164,31 @@ end
 @non_differentiable ctmrg_logfinish!(args...)
 @non_differentiable ctmrg_logcancel!(args...)
 
-# TODO: we might want to consider embedding the smaller tensor into the larger space and then compute the difference
 """
-    _singular_value_distance((S₁, S₂))
+    _singular_value_distance(S₁, S₂)
 
 Compute the singular value distance as an error measure, e.g. for CTMRG iterations.
 To that end, the singular values of the current iteration `S₁` are compared with the
 previous one `S₂`. When the virtual spaces change, this comparison is not directly possible
 such that both tensors are projected into the smaller space and then subtracted.
 """
-function _singular_value_distance((S₁, S₂))
-    V₁ = space(S₁, 1)
-    V₂ = space(S₂, 1)
-    if V₁ == V₂
-        return norm(S₁ - S₂)
-    else
-        V = infimum(V₁, V₂)
-        e1 = isometry(V₁, V)
-        e2 = isometry(V₂, V)
-        return norm(e1' * S₁ * e1 - e2' * S₂ * e2)
+function _singular_value_distance(S₁::SV, S₂::SV) where {SV <: TensorKit.SectorVector}
+    # allocate vector for difference - possibly grow
+    V₁ = Vect[sectortype(S₁)](c => length(v) for (c, v) in blocks(S₁))
+    V₂ = Vect[sectortype(S₂)](c => length(v) for (c, v) in blocks(S₂))
+    diff = zerovector!(SV(undef, supremum(V₁, V₂)))
+
+    for (c, b) in blocks(S₁)
+        diff[c][1:length(b)] .= b
     end
+    for (c, b) in blocks(S₂)
+        diff[c][1:length(b)] .-= b
+    end
+
+    return norm(diff)
 end
+_singular_value_distance(S₁::DiagonalTensorMap, S₂::DiagonalTensorMap) =
+    _singular_value_distance(diagview(S₁), diagview(S₂))
 
 """
     calc_convergence(env, CS_old, TS_old)
@@ -196,10 +200,10 @@ This determined either from the previous corner and edge singular values
 """
 function calc_convergence(env, CS_old, TS_old)
     CS_new = map(svd_vals, env.corners)
-    ΔCS = maximum(_singular_value_distance, zip(CS_old, CS_new))
+    ΔCS = maximum(splat(_singular_value_distance), zip(CS_old, CS_new))
 
     TS_new = map(svd_vals, env.edges)
-    ΔTS = maximum(_singular_value_distance, zip(TS_old, TS_new))
+    ΔTS = maximum(splat(_singular_value_distance), zip(TS_old, TS_new))
 
     @debug "maxᵢ|Cⁿ⁺¹ - Cⁿ|ᵢ = $ΔCS   maxᵢ|Tⁿ⁺¹ - Tⁿ|ᵢ = $ΔTS"
 
diff --git a/src/algorithms/ctmrg/gaugefix.jl b/src/algorithms/ctmrg/gaugefix.jl
index 72f10bf63..f3d869d3b 100644
--- a/src/algorithms/ctmrg/gaugefix.jl
+++ b/src/algorithms/ctmrg/gaugefix.jl
@@ -80,8 +80,8 @@ function gauge_fix(envfinal::CTMRGEnv{C, T}, envprev::CTMRGEnv{C, T}, ::Scrambli
         ρfinal = transfermatrix_fixedpoint(Tsfinal, M, ρinit)
 
         # Decompose and multiply
-        Qprev, = left_orth!(ρprev)
-        Qfinal, = left_orth!(ρfinal)
+        Qprev, = left_orth!(ρprev; positive = true)
+        Qfinal, = left_orth!(ρfinal; positive = true)
 
         return Qprev * Qfinal'
     end
@@ -115,8 +115,8 @@ function gauge_fix(envfinal::CTMRGEnv{C, T}, envprev::CTMRGEnv{C, T}, ::Scrambli
     ρfinal = c4v_transfermatrix_fixedpoint(Tfinal, M, ρinit)
 
     # Decompose and multiply
-    Qprev, = left_orth!(ρprev)
-    Qfinal, = left_orth!(ρfinal)
+    Qprev, = left_orth!(ρprev; positive = true)
+    Qfinal, = left_orth!(ρfinal; positive = true)
 
     σ = Qprev * Qfinal'
 
diff --git a/src/algorithms/ctmrg/sequential.jl b/src/algorithms/ctmrg/sequential.jl
index afb24f540..1c9631675 100644
--- a/src/algorithms/ctmrg/sequential.jl
+++ b/src/algorithms/ctmrg/sequential.jl
@@ -86,10 +86,7 @@ function sequential_projectors(col::Int, network, env::CTMRGEnv, alg::ProjectorA
     )
     proj_and_info = similar(coordinates, T_dst)
     proj_and_info′::typeof(proj_and_info) = dtmap!!(proj_and_info, coordinates) do (r, c)
-        trunc = truncation_strategy(alg, env.edges[WEST, _prev(r, size(env, 2)), c])
-        proj, info = sequential_projectors(
-            (WEST, r, c), network, env, @set(alg.trunc = trunc)
-        )
+        proj, info = sequential_projectors((WEST, r, c), network, env, alg)
         return proj, info
     end
     return _split_proj_and_info(proj_and_info′)
@@ -99,9 +96,11 @@ function sequential_projectors(
     )
     _, r, c = coordinate
     r′ = _prev(r, size(env, 2))
+    trunc = truncation_strategy(alg, env.edges[WEST, r′, c])
+    alg′ = @set alg.trunc = trunc
     Q1 = TensorMap(EnlargedCorner(network, env, (SOUTHWEST, r, c)))
     Q2 = TensorMap(EnlargedCorner(network, env, (NORTHWEST, r′, c)))
-    return compute_projector((Q1, Q2), coordinate, alg)
+    return compute_projector((Q1, Q2), coordinate, alg′)
 end
 function sequential_projectors(
         coordinate::NTuple{3, Int}, network, env::CTMRGEnv, alg::FullInfiniteProjector
@@ -110,13 +109,15 @@ function sequential_projectors(
     coordinate_nw = _next_coordinate(coordinate, rowsize, colsize)
     coordinate_ne = _next_coordinate(coordinate_nw, rowsize, colsize)
     coordinate_se = _next_coordinate(coordinate_ne, rowsize, colsize)
+    trunc = truncation_strategy(alg, env.edges[WEST, coordinate_nw[2:3]...])
+    alg′ = @set alg.trunc = trunc
     ec = (
         TensorMap(EnlargedCorner(network, env, coordinate_se)),
         TensorMap(EnlargedCorner(network, env, coordinate)),
         TensorMap(EnlargedCorner(network, env, coordinate_nw)),
         TensorMap(EnlargedCorner(network, env, coordinate_ne)),
     )
-    return compute_projector(ec, coordinate, alg)
+    return compute_projector(ec, coordinate, alg′)
 end
 
 """
diff --git a/src/algorithms/ctmrg/simultaneous.jl b/src/algorithms/ctmrg/simultaneous.jl
index 34ac95a16..731ef57de 100644
--- a/src/algorithms/ctmrg/simultaneous.jl
+++ b/src/algorithms/ctmrg/simultaneous.jl
@@ -102,13 +102,12 @@ end
 function simultaneous_projectors(
         coordinate, enlarged_corners::Array{E, 3}, env, alg::FullInfiniteProjector
     ) where {E}
-    coordinate′ = _next_coordinate(coordinate, size(env)[2:3]...)
-    trunc = truncation_strategy(alg, env.edges[coordinate[1], coordinate′[2:3]...])
-    alg′ = @set alg.trunc = trunc
-    rowsize, colsize = size(enlarged_corners)[2:3]
+    rowsize, colsize = size(env)[2:3]
     coordinate2 = _next_coordinate(coordinate, rowsize, colsize)
     coordinate3 = _next_coordinate(coordinate2, rowsize, colsize)
     coordinate4 = _next_coordinate(coordinate3, rowsize, colsize)
+    trunc = truncation_strategy(alg, env.edges[coordinate[1], coordinate2[2:3]...])
+    alg′ = @set alg.trunc = trunc
     ec = (
         enlarged_corners[coordinate4...],
         enlarged_corners[coordinate...],
diff --git a/src/algorithms/time_evolution/evoltools.jl b/src/algorithms/time_evolution/evoltools.jl
index 146ede367..8c9d34945 100644
--- a/src/algorithms/time_evolution/evoltools.jl
+++ b/src/algorithms/time_evolution/evoltools.jl
@@ -143,8 +143,8 @@ function _qr_bond(A::PT, B::PT; gate_ax::Int = 1) where {PT <: Union{PEPSTensor,
             ((1, 3, 5, 6), (2, 4)), ((1, 3, 4, 5), (2, 6)), (1, 2, 5, 3, 4), Tuple(1:5)
         end
     end
-    X, a = left_orth(permute(A, permA))
-    Y, b = left_orth(permute(B, permB))
+    X, a = left_orth(permute(A, permA); positive = true)
+    Y, b = left_orth(permute(B, permB); positive = true)
     # no longer needed after TensorKit 0.15
     # @assert !isdual(space(a, 1))
     # @assert !isdual(space(b, 1))
@@ -226,11 +226,7 @@ function _apply_gate(
     else
         trunc
     end
-
-    # TODO: replace this with actual truncation error once TensorKit is updated
-    ac, sc, bc = svd_compact!(a2b2; alg = LAPACK_QRIteration())
-    a, s, b, ϵ = _truncate_compact((ac, sc, bc), trunc)
-
+    a, s, b, ϵ = svd_trunc!(a2b2; trunc, alg = LAPACK_QRIteration())
     a, b = absorb_s(a, s, b)
     if need_flip
         a, s, b = flip(a, numind(a)), _fliptwist_s(s), flip(b, 1)
diff --git a/src/algorithms/time_evolution/simpleupdate3site.jl b/src/algorithms/time_evolution/simpleupdate3site.jl
index 895df8dca..bbe432807 100644
--- a/src/algorithms/time_evolution/simpleupdate3site.jl
+++ b/src/algorithms/time_evolution/simpleupdate3site.jl
@@ -125,7 +125,7 @@ function qr_through(
     @assert !isdual(codomain(R0, 1))
     @assert !isdual(domain(M, 1)) && !isdual(codomain(M, 1))
     @tensor A[-1 -2 -3 -4; -5] := R0[-1; 1] * M[1 -2 -3 -4; -5]
-    _, r = left_orth!(A)
+    _, r = left_orth!(A; positive = true)
     normalize && normalize!(r, Inf)
     return r
 end
@@ -134,7 +134,7 @@ function qr_through(
         ::Nothing, M::GenericMPSTensor{S, 4}; normalize::Bool = true
     ) where {S <: ElementarySpace}
     @assert !isdual(domain(M, 1))
-    _, r = left_orth(M)
+    _, r = left_orth(M; positive = true)
     normalize && normalize!(r, Inf)
     return r
 end
@@ -153,7 +153,7 @@ function lq_through(
     @assert !isdual(domain(L1, 1))
     @assert !isdual(codomain(M, 1)) && !isdual(domain(M, 1))
     @tensor A[-1; -2 -3 -4 -5] := M[-1 -2 -3 -4; 1] * L1[1; -5]
-    l, _ = right_orth!(A)
+    l, _ = right_orth!(A; positive = true)
     normalize && normalize!(l, Inf)
     return l
 end
@@ -163,7 +163,7 @@ function lq_through(
     ) where {S <: ElementarySpace}
     @assert !isdual(codomain(M, 1))
     A = permute(M, ((1,), (2, 3, 4, 5)))
-    l, _ = right_orth!(A)
+    l, _ = right_orth!(A; positive = true)
     normalize && normalize!(l, Inf)
     return l
 end
@@ -207,11 +207,7 @@ function _proj_from_RL(
     @assert isdual(domain(r, 1)) == isdual(codomain(r, 1)) == false
     @assert isdual(domain(l, 1)) == isdual(codomain(l, 1)) == false
     rl = r * l
-
-    # TODO: replace this with actual truncation error once TensorKit is updated
-    uc, sc, vhc = svd_compact!(rl)
-    u, s, vh, ϵ = _truncate_compact((uc, sc, vhc), trunc)
-
+    u, s, vh, ϵ = svd_trunc!(rl; trunc)
     sinv = sdiag_pow(s, -1 / 2)
     Pa, Pb = l * vh' * sinv, sinv * u' * r
     return Pa, s, Pb, ϵ
diff --git a/src/algorithms/truncation/bond_truncation.jl b/src/algorithms/truncation/bond_truncation.jl
index 9d70947b1..d8feb8cb3 100644
--- a/src/algorithms/truncation/bond_truncation.jl
+++ b/src/algorithms/truncation/bond_truncation.jl
@@ -112,7 +112,7 @@ function bond_truncate(
         _, 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
+        Δs = _singular_value_distance(s, s0) / s_nrm
         Δcost = abs(cost - cost0) / cost00
         Δfid = abs(fid - fid0)
         cost0, fid0, s0 = cost, fid, s
@@ -160,8 +160,8 @@ function bond_truncate(
         --- a == b ---   ==>   - Qa ← Ra == Rb ← Qb -
             ↓    ↓               ↓               ↓
     =#
-    Qa, Ra = left_orth(a)
-    Rb, Qb = right_orth(b)
+    Qa, Ra = left_orth(a; positive = true)
+    Rb, Qb = right_orth(b; positive = true)
     @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`
diff --git a/src/algorithms/truncation/fullenv_truncation.jl b/src/algorithms/truncation/fullenv_truncation.jl
index 797d9532a..ea339c8d8 100644
--- a/src/algorithms/truncation/fullenv_truncation.jl
+++ b/src/algorithms/truncation/fullenv_truncation.jl
@@ -256,7 +256,7 @@ function fullenv_truncate(
         u, s, vh = svd_trunc!(b1; trunc = alg.trunc)
         # determine convergence
         s_nrm = norm(s0, Inf)
-        Δs = ((space(s) == space(s0)) ? _singular_value_distance((s, s0)) : NaN) / s_nrm
+        Δs = _singular_value_distance(s, s0) / s_nrm
         Δfid = fid - fid0
         s0 = deepcopy(s)
         fid0 = fid
diff --git a/src/environments/suweight.jl b/src/environments/suweight.jl
index ea4491aec..2966d4155 100644
--- a/src/environments/suweight.jl
+++ b/src/environments/suweight.jl
@@ -151,7 +151,7 @@ end
 
 function compare_weights(wts1::SUWeight, wts2::SUWeight)
     @assert size(wts1) == size(wts2)
-    return sum(_singular_value_distance, zip(wts1.data, wts2.data)) / length(wts1)
+    return sum(splat(_singular_value_distance), zip(wts1.data, wts2.data)) / length(wts1)
 end
 
 function Base.show(io::IO, ::MIME"text/plain", wts::SUWeight)
diff --git a/src/utility/eigh.jl b/src/utility/eigh.jl
index 05a10784d..fab3c0f94 100644
--- a/src/utility/eigh.jl
+++ b/src/utility/eigh.jl
@@ -153,7 +153,8 @@ function _eigh_trunc!(
         trunc::TruncationStrategy,
     )
     D, V = eigh_full!(t; alg)
-    D̃, Ṽ, truncerror = _truncate_eigh((D, V), trunc)
+    (D̃, Ṽ), ind = truncate(eigh_trunc!, (D, V), trunc)
+    truncerror = truncation_error(diagview(D), ind)
 
     # construct info NamedTuple
     condnum = cond(D)
@@ -163,18 +164,6 @@ function _eigh_trunc!(
     return D̃, Ṽ, info
 end
 
-# hacky way of computing the truncation error for current version of eigh_trunc!
-# TODO: replace once TensorKit updates to new MatrixAlgebraKit which returns truncation error as well
-function _truncate_eigh((D, V), trunc::TruncationStrategy)
-    if !(trunc isa NoTruncation) && !isempty(blocksectors(D))
-        D̃, Ṽ = truncate(eigh_trunc!, (D, V), trunc)[1]
-        truncerror = sqrt(abs(norm(D)^2 - norm(D̃)^2))
-        return D̃, Ṽ, truncerror
-    else
-        return D, V, zero(real(scalartype(D)))
-    end
-end
-
 """
 $(TYPEDEF)
 
diff --git a/src/utility/svd.jl b/src/utility/svd.jl
index b2e8e6892..6166f0c96 100644
--- a/src/utility/svd.jl
+++ b/src/utility/svd.jl
@@ -142,7 +142,8 @@ function _svd_trunc!(
         trunc::TruncationStrategy,
     )
     U, S, V⁺ = svd_compact!(t; alg)
-    Ũ, S̃, Ṽ⁺, truncerror = _truncate_compact((U, S, V⁺), trunc)
+    (Ũ, S̃, Ṽ⁺), ind = truncate(svd_trunc!, (U, S, V⁺), trunc)
+    truncerror = truncation_error(diagview(S), ind)
 
     # construct info NamedTuple
     condnum = cond(S)
@@ -152,18 +153,6 @@ function _svd_trunc!(
     return Ũ, S̃, Ṽ⁺, info
 end
 
-# hacky way of computing the truncation error for current version of svd_trunc!
-# TODO: replace once TensorKit updates to new MatrixAlgebraKit which returns truncation error as well
-function _truncate_compact((U, S, V⁺), trunc::TruncationStrategy)
-    if !(trunc isa NoTruncation) && !isempty(blocksectors(S))
-        Ũ, S̃, Ṽ⁺ = truncate(svd_trunc!, (U, S, V⁺), trunc)[1]
-        truncerror = sqrt(norm(S)^2 - norm(S̃)^2)
-        return Ũ, S̃, Ṽ⁺, truncerror
-    else
-        return U, S, V⁺, zero(real(scalartype(S)))
-    end
-end
-
 """
 $(TYPEDEF)
 
diff --git a/test/ctmrg/flavors.jl b/test/ctmrg/flavors.jl
index 023d985aa..a39af1985 100644
--- a/test/ctmrg/flavors.jl
+++ b/test/ctmrg/flavors.jl
@@ -33,20 +33,12 @@ eigh_algs = [:qriteration, :lanczos]
     # compare singular values
     CS_sequential = map(svd_vals, env_sequential.corners)
     CS_simultaneous = map(svd_vals, env_simultaneous.corners)
-    ΔCS = maximum(zip(CS_sequential, CS_simultaneous)) do (c_lm, c_as)
-        smallest = infimum(MPSKit._firstspace(c_lm), MPSKit._firstspace(c_as))
-        e_old = isometry(MPSKit._firstspace(c_lm), smallest)
-        e_new = isometry(MPSKit._firstspace(c_as), smallest)
-        return norm(e_new' * c_as * e_new - e_old' * c_lm * e_old)
-    end
+    ΔCS = maximum(splat(PEPSKit._singular_value_distance), zip(CS_sequential, CS_simultaneous))
     @test ΔCS < 1.0e-2
 
     TS_sequential = map(svd_vals, env_sequential.edges)
     TS_simultaneous = map(svd_vals, env_simultaneous.edges)
-    ΔTS = maximum(zip(TS_sequential, TS_simultaneous)) do (t_lm, t_as)
-        MPSKit._firstspace(t_lm) == MPSKit._firstspace(t_as) || return scalartype(t_lm)(Inf)
-        return norm(t_as - t_lm)
-    end
+    ΔTS = maximum(splat(PEPSKit._singular_value_distance), zip(TS_sequential, TS_simultaneous))
     @test ΔTS < 1.0e-2
 
     # compare Heisenberg energies