diff --git a/Project.toml b/Project.toml
index f5a00ec63..27b180f8b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -25,7 +25,7 @@ VectorInterface = "409d34a3-91d5-4945-b6ec-7529ddf182d8"
 [compat]
 Accessors = "0.1"
 Aqua = "0.8.9"
-BlockTensorKit = "0.1.6"
+BlockTensorKit = "0.2"
 Combinatorics = "1"
 Compat = "3.47, 4.10"
 DocStringExtensions = "0.9.3"
diff --git a/docs/src/assets/mpskit.bib b/docs/src/assets/mpskit.bib
index 11b22006a..41320c7ed 100644
--- a/docs/src/assets/mpskit.bib
+++ b/docs/src/assets/mpskit.bib
@@ -196,6 +196,22 @@ @article{paeckel2019
   keywords = {Density matrix renormalization group (DMRG),Matrix-product states (MPS),Strongly-correlated systems,Time-dependent variational principle (TDVP),Time-evolution methods,Time-evolving block decimation (TEBD)}
 }
 
+@article{parker2020,
+  title = {Local Matrix Product Operators: {{Canonical}} Form, Compression, and Control Theory},
+  shorttitle = {Local Matrix Product Operators},
+  author = {Parker, Daniel E. and Cao, Xiangyu and Zaletel, Michael P.},
+  year = {2020},
+  month = jul,
+  journal = {Physical Review B},
+  volume = {102},
+  number = {3},
+  pages = {035147},
+  publisher = {American Physical Society},
+  doi = {10.1103/PhysRevB.102.035147},
+  url = {https://link.aps.org/doi/10.1103/PhysRevB.102.035147},
+  abstract = {We present a new method for compressing matrix product operators (MPOs) which represent sums of local terms, such as Hamiltonians. Just as with area law states, such local operators may be fully specified with a small amount of information per site. Standard matrix product state (MPS) tools are ill-suited to this case, due to extensive Schmidt values that coexist with intensive ones, and Jordan blocks in the transfer matrix. We ameliorate these issues by introducing an ``almost Schmidt decomposition'' that respects locality. Our method is ``{$\varepsilon$}-close'' to the accuracy of MPS-based methods for finite MPOs, and extends seamlessly to the thermodynamic limit, where MPS techniques are inapplicable. In the framework of control theory, our method generalizes Kung's algorithm for model order reduction. Several examples are provided, including an all-MPO version of the operator recursion method (Lanczos algorithm) directly in the thermodynamic limit. All results are accompanied by practical algorithms, well-suited for the large MPOs that arise in DMRG for long-range or quasi-2D models.}
+}
+
 @article{pollmann2009,
   title = {Theory of {{Finite-Entanglement Scaling}} at {{One-Dimensional Quantum Critical Points}}},
   author = {Pollmann, Frank and Mukerjee, Subroto and Turner, Ari M. and Moore, Joel E.},
diff --git a/src/MPSKit.jl b/src/MPSKit.jl
index 4bab6bf82..51ee82e32 100644
--- a/src/MPSKit.jl
+++ b/src/MPSKit.jl
@@ -115,6 +115,7 @@ include("operators/abstractmpo.jl")
 include("operators/mpo.jl")
 include("operators/jordanmpotensor.jl")
 include("operators/mpohamiltonian.jl") # the mpohamiltonian objects
+include("operators/ortho.jl")
 include("operators/multilinempo.jl")
 include("operators/projection.jl")
 include("operators/timedependence.jl")
diff --git a/src/algorithms/changebonds/svdcut.jl b/src/algorithms/changebonds/svdcut.jl
index cd60bf1b8..7ad26e706 100644
--- a/src/algorithms/changebonds/svdcut.jl
+++ b/src/algorithms/changebonds/svdcut.jl
@@ -1,11 +1,18 @@
 """
 $(TYPEDEF)
 
-An algorithm that uses truncated SVD to change the bond dimension of a ψ.
+An algorithm that uses truncated SVD to change the bond dimension of a state or operator.
+This is achieved by a sweeping algorithm that locally performs (optimal) truncations in a gauged basis.
+
+See also [`changebonds(!)`](@ref changebonds)
 
 ## Fields
 
 $(TYPEDFIELDS)
+
+## References
+
+* [Parker et al. Phys. Rev. B 102 (2020)](@cite parker2020)
 """
 @kwdef struct SvdCut{S} <: Algorithm
     "algorithm used for the singular value decomposition"
@@ -105,3 +112,22 @@ end
 function changebonds(ψ, H, alg::SvdCut, envs = environments(ψ, H))
     return changebonds(ψ, alg), envs
 end
+
+changebonds(mpo::FiniteMPOHamiltonian, alg::SvdCut) = changebonds!(copy(mpo), alg)
+function changebonds!(H::FiniteMPOHamiltonian, alg::SvdCut)
+    # orthogonality center to the left
+    for i in length(H):-1:2
+        H = right_canonicalize!(H, i)
+    end
+
+    # swipe right
+    for i in 1:(length(H) - 1)
+        H = left_canonicalize!(H, i; alg = alg.alg_svd, alg.trscheme)
+    end
+    # swipe left -- TODO: do we really need this double sweep?
+    for i in length(H):-1:2
+        H = right_canonicalize!(H, i; alg = alg.alg_svd, alg.trscheme)
+    end
+
+    return H
+end
diff --git a/src/operators/jordanmpotensor.jl b/src/operators/jordanmpotensor.jl
index e46356f35..7fcc7c6e8 100644
--- a/src/operators/jordanmpotensor.jl
+++ b/src/operators/jordanmpotensor.jl
@@ -67,26 +67,37 @@ end
 
 function JordanMPOTensor(
         V::TensorMapSumSpace{S, 2, 2},
-        A::SparseBlockTensorMap{TA, E, S, 2, 2},
-        B::SparseBlockTensorMap{TB, E, S, 2, 1},
-        C::SparseBlockTensorMap{TC, E, S, 1, 2},
-        D::SparseBlockTensorMap{TD, E, S, 1, 1}
+        A::SparseBlockTensorMap{TA, E, S, 2, 2}, B::SparseBlockTensorMap{TB, E, S, 2, 1},
+        C::SparseBlockTensorMap{TC, E, S, 1, 2}, D::SparseBlockTensorMap{TD, E, S, 1, 1}
     ) where {E, S, TA, TB, TC, TD}
     allVs = eachspace(V)
     VA = space(allVs[2:(end - 1), 1, 1, 2:(end - 1)])
-    VA == space(A) || throw(SpaceMismatch("A-block has incompatible spaces"))
+    VA == space(A) || throw(SpaceMismatch(lazy"A-block has incompatible spaces:\n$VA\n$(space(A))"))
 
     VB = removeunit(space(allVs[2:(end - 1), 1, 1, end]), 4)
-    VB == space(B) || throw(SpaceMismatch("B-block has incompatible spaces"))
+    VB == space(B) || throw(SpaceMismatch(lazy"B-block has incompatible spaces:\n$VB\n$(space(B))"))
 
     VC = removeunit(space(allVs[1, 1, 1, 2:(end - 1)]), 1)
-    VC == space(C) || throw(SpaceMismatch("C-block has incompatible spaces"))
+    VC == space(C) || throw(SpaceMismatch(lazy"C-block has incompatible spaces:\n$VC\n$(space(C))"))
 
     VD = removeunit(removeunit(space(allVs[1, 1, 1, end:end]), 4), 1)
-    VD == space(D) || throw(SpaceMismatch("D-block has incompatible spaces"))
+    VD == space(D) || throw(SpaceMismatch(lazy"D-block has incompatible spaces:\n$VD\n$(space(D))"))
 
     return JordanMPOTensor{E, S, TA, TB, TC, TD}(V, A, B, C, D)
 end
+function JordanMPOTensor(
+        V::TensorMapSumSpace{S, 2, 2},
+        A::AbstractTensorMap{E, S, 2, 2}, B::AbstractTensorMap{E, S, 2, 1},
+        C::AbstractTensorMap{E, S, 1, 2}, D::AbstractTensorMap{E, S, 1, 1}
+    ) where {E, S}
+    return JordanMPOTensor(
+        V,
+        A isa SparseBlockTensorMap ? A : SparseBlockTensorMap(A),
+        B isa SparseBlockTensorMap ? B : SparseBlockTensorMap(B),
+        C isa SparseBlockTensorMap ? C : SparseBlockTensorMap(C),
+        D isa SparseBlockTensorMap ? D : SparseBlockTensorMap(D)
+    )
+end
 
 function JordanMPOTensor(W::SparseBlockTensorMap{TT, E, S, 2, 2}) where {TT, E, S}
     @assert W[1, 1, 1, 1] isa BraidingTensor && W[end, 1, 1, end] isa BraidingTensor
diff --git a/src/operators/ortho.jl b/src/operators/ortho.jl
new file mode 100644
index 000000000..1b7ecb924
--- /dev/null
+++ b/src/operators/ortho.jl
@@ -0,0 +1,190 @@
+function left_canonicalize!(
+        H::FiniteMPOHamiltonian, i::Int;
+        alg = QRpos(), trscheme::TruncationScheme = notrunc()
+    )
+    1 ≤ i < length(H) || throw(ArgumentError("Bounds error in canonicalize"))
+
+    W = H[i]
+    S = spacetype(W)
+    d = sqrt(dim(physicalspace(W)))
+
+    # orthogonalize second column against first
+    WI = removeunit(W[1, 1, 1, 1], 1)
+    @tensor t[l; r] := conj(WI[p; p' l]) * W.C[p; p' r]
+    # TODO: the following is currently broken due to a TensorKit bug
+    # @plansor C′[p; p' r] := W.C[p; p' r] - WI[p; p' l] * t[l; r]
+    @plansor C′[p; p' r] := -WI[p; p' l] * t[l; r]
+    add!(C′, W.C)
+
+    # QR of second column
+    if size(W, 1) == 1
+        tmp = transpose(C′, ((2, 1), (3,)))
+
+        if trscheme == notrunc()
+            Q, R = leftorth!(tmp; alg)
+        else
+            @assert alg == SVD() || alg == SDD()
+            Q, Σ, Vᴴ = tsvd!(tmp; alg, trunc = trscheme)
+            R = Σ * Vᴴ
+        end
+
+        if dim(R) == 0 # fully truncated
+            V = BlockTensorKit.oplus(oneunit(S), oneunit(S))
+            Q1 = typeof(W.A)(undef, SumSpace{S}() ⊗ physicalspace(W) ← physicalspace(W) ⊗ SumSpace{S}())
+            Q2 = typeof(W.C)(undef, physicalspace(W) ← physicalspace(W) ⊗ SumSpace{S}())
+        else
+            V = BlockTensorKit.oplus(oneunit(S), space(R, 1), oneunit(S))
+            scale!(Q, d)
+            scale!(R, inv(d))
+            Q1 = typeof(W.A)(undef, SumSpace{S}() ⊗ physicalspace(W) ← physicalspace(W) ⊗ space(R, 1))
+            Q2 = transpose(Q, ((2,), (1, 3)))
+        end
+        H[i] = JordanMPOTensor(codomain(W) ← physicalspace(W) ⊗ V, Q1, W.B, Q2, W.D)
+    else
+        tmp = transpose(cat(insertleftunit(C′, 1), W.A; dims = 1), ((3, 1, 2), (4,)))
+        if trscheme == notrunc()
+            Q, R = leftorth!(tmp; alg)
+        else
+            @assert alg == SVD() || alg == SDD()
+            Q, Σ, Vᴴ = tsvd!(tmp; alg, trunc = trscheme)
+            R = Σ * Vᴴ
+        end
+        if dim(R) == 0 # fully truncated
+            V = BlockTensorKit.oplus(oneunit(S), oneunit(S))
+            Q1 = typeof(W.A)(undef, SumSpace{S}() ⊗ physicalspace(W) ← physicalspace(W) ⊗ SumSpace{S}())
+            Q2 = typeof(W.C)(undef, physicalspace(W) ← physicalspace(W) ⊗ SumSpace{S}())
+        else
+            scale!(Q, d)
+            scale!(R, inv(d))
+            Q′ = transpose(Q, ((2, 3), (1, 4)))
+            Q1 = Q′[2:end, 1, 1, 1]
+            Q2 = removeunit(SparseBlockTensorMap(Q′[1:1, 1, 1, 1]), 1)
+            V = BlockTensorKit.oplus(oneunit(S), right_virtualspace(Q′), oneunit(S))
+        end
+        H[i] = JordanMPOTensor(codomain(W) ← physicalspace(W) ⊗ V, Q1, W.B, Q2, W.D)
+    end
+
+    # absorb into next site
+    W′ = H[i + 1]
+
+    if size(W′, 4) > 1 && dim(R) != 0
+        @plansor A′[l p; p' r] := R[l; r'] * W′.A[r' p; p' r]
+    else
+        A′ = similar(W′.A, right_virtualspace(H[i].A) ⊗ physicalspace(W′) ← domain(W′.A))
+    end
+
+    if size(W′, 4) > 1
+        @plansor C′[l p; p' r] := t[l; r'] * W′.A[r' p; p' r]
+        C′ = add!(removeunit(C′, 1), W′.C)
+    else
+        C′ = W′.C # empty
+    end
+
+    if dim(R) != 0
+        @plansor B′[l p; p'] := R[l; r] * W′.B[r p; p']
+    else
+        B′ = similar(W′.B, right_virtualspace(H[i].A) ⊗ physicalspace(W′) ← domain(W′.B))
+    end
+
+    @plansor D′[l p; p'] := t[l; r] * W′.B[r p; p']
+    D′ = add!(removeunit(D′, 1), W′.D)
+
+    H[i + 1] = JordanMPOTensor(
+        right_virtualspace(H[i]) ⊗ physicalspace(W′) ← domain(W′),
+        A′, B′, C′, D′
+    )
+    return H
+end
+
+function right_canonicalize!(
+        H::FiniteMPOHamiltonian, i::Int;
+        alg = LQpos(), trscheme::TruncationScheme = notrunc()
+    )
+    1 < i ≤ length(H) || throw(ArgumentError("Bounds error in canonicalize"))
+
+    W = H[i]
+    S = spacetype(W)
+    d = sqrt(dim(physicalspace(W)))
+
+    # orthogonalize second row against last
+    WI = removeunit(W[end, 1, 1, end], 4)
+    @plansor t[l; r] := conj(WI[r p; p']) * W.B[l p; p']
+    # TODO: the following is currently broken due to a TensorKit bug
+    # @plansor B′[l p; p'] := W.B[l p; p'] - WI[r p; p'] * t[l; r]
+    @plansor B′[l p; p'] := -WI[r p; p'] * t[l; r]
+    add!(B′, W.B)
+
+    # LQ of second row
+    if size(W, 4) == 1
+        tmp = transpose(B′, ((1,), (3, 2)))
+        if trscheme == notrunc()
+            R, Q = rightorth!(tmp; alg)
+        else
+            @assert alg == SVD() || alg == SDD()
+            U, Σ, Q = tsvd!(tmp; alg, trunc = trscheme)
+            R = U * Σ
+        end
+
+        if dim(R) == 0
+            V = BlockTensorKit.oplus(oneunit(S), oneunit(S))
+            Q1 = typeof(W.A)(undef, SumSpace{S}() ⊗ physicalspace(W) ← physicalspace(W) ⊗ SumSpace{S}())
+            Q2 = typeof(W.B)(undef, SumSpace{S}() ⊗ physicalspace(W) ← physicalspace(W))
+        else
+            V = BlockTensorKit.oplus(oneunit(S), space(Q, 1), oneunit(S))
+            scale!(Q, d)
+            scale!(R, inv(d))
+            Q1 = typeof(W.A)(undef, space(Q, 1) ⊗ physicalspace(W) ← physicalspace(W) ⊗ SumSpace{S}())
+            Q2 = transpose(Q, ((1, 3), (2,)))
+        end
+        H[i] = JordanMPOTensor(V ⊗ physicalspace(W) ← domain(W), Q1, Q2, W.C, W.D)
+    else
+        tmp = transpose(cat(insertleftunit(B′, 4), W.A; dims = 4), ((1,), (3, 4, 2)))
+        if trscheme == notrunc()
+            R, Q = rightorth!(tmp; alg)
+        else
+            @assert alg == SVD() || alg == SDD()
+            U, Σ, Q = tsvd!(tmp; alg, trunc = trscheme)
+            R = U * Σ
+        end
+        if dim(R) == 0
+            V = BlockTensorKit.oplus(oneunit(S), oneunit(S))
+            Q1 = typeof(W.A)(undef, SumSpace{S}() ⊗ physicalspace(W) ← physicalspace(W) ⊗ SumSpace{S}())
+            Q2 = typeof(W.B)(undef, SumSpace{S}() ⊗ physicalspace(W) ← physicalspace(W))
+        else
+            scale!(Q, d)
+            scale!(R, inv(d))
+            Q′ = transpose(Q, ((1, 4), (2, 3)))
+            Q1 = SparseBlockTensorMap(Q′[1, 1, 1, 2:end])
+            Q2 = removeunit(SparseBlockTensorMap(Q′[1, 1, 1, 1:1]), 4)
+            V = BlockTensorKit.oplus(oneunit(S), left_virtualspace(Q′), oneunit(S))
+        end
+        H[i] = JordanMPOTensor(V ⊗ physicalspace(W) ← domain(W), Q1, Q2, W.C, W.D)
+    end
+
+    # absorb into previous site
+    W′ = H[i - 1]
+
+    if size(W′, 1) > 1 && dim(R) != 0
+        @plansor A′[l p; p' r] := W′.A[l p; p' r'] * R[r'; r]
+    else
+        A′ = similar(W′.A, codomain(W′.A) ← physicalspace(W′.A) ⊗ left_virtualspace(H[i].A))
+    end
+
+    if size(W′, 1) > 1
+        @plansor B′[l p; p' r] := W′.A[l p; p' r'] * t[r'; r]
+        B′ = add!(removeunit(B′, 4), W′.B)
+    else
+        B′ = W′.B
+    end
+
+    if dim(R) != 0
+        @plansor C′[p; p' r] := W′.C[p; p' r'] * R[r'; r]
+    else
+        C′ = similar(W′.C, codomain(W′.C) ← physicalspace(W′) ⊗ left_virtualspace(H[i].A))
+    end
+
+    @plansor D′[p; p' r] := W′.C[p; p' r'] * t[r'; r]
+    D′ = add!(removeunit(D′, 3), W′.D)
+    H[i - 1] = JordanMPOTensor(codomain(W′) ← physicalspace(W′) ⊗ V, A′, B′, C′, D′)
+    return H
+end
diff --git a/test/operators.jl b/test/operators.jl
index 12f3a9706..3ea1ae003 100644
--- a/test/operators.jl
+++ b/test/operators.jl
@@ -124,11 +124,13 @@ module TestOperators
     @testset "Finite MPOHamiltonian" begin
         L = 3
         T = ComplexF64
-        for V in (ℂ^2, U1Space(0 => 1, 1 => 1))
+        for V in (ℂ^2, U1Space(-1 => 1, 0 => 1, 1 => 1))
             lattice = fill(V, L)
-            O₁ = rand(T, V, V)
+            O₁ = randn(T, V, V)
+            O₁ += O₁'
             E = id(storagetype(O₁), domain(O₁))
-            O₂ = rand(T, V^2 ← V^2)
+            O₂ = randn(T, V^2 ← V^2)
+            O₂ += O₂'
 
             H1 = FiniteMPOHamiltonian(lattice, i => O₁ for i in 1:L)
             H2 = FiniteMPOHamiltonian(lattice, (i, i + 1) => O₂ for i in 1:(L - 1))
@@ -164,6 +166,9 @@ module TestOperators
                 FiniteMPOHamiltonian(lattice, 3 => O₁)
             @test 0.8 * H1 + 0.2 * H1 ≈ H1 atol = 1.0e-6
             @test convert(TensorMap, H1 + H2) ≈ convert(TensorMap, H1) + convert(TensorMap, H2) atol = 1.0e-6
+            H1_trunc = changebonds(H1, SvdCut(; trscheme = truncdim(0)))
+            @test H1_trunc ≈ H1
+            @test all(left_virtualspace(H1_trunc) .== left_virtualspace(H1))
 
             # test dot and application
             state = rand(T, prod(lattice))
@@ -181,12 +186,21 @@ module TestOperators
             vertical_operators = Dict(
                 (I, I + I_vertical) => O₂ for I in eachindex(IndexCartesian(), square) if I[1] < size(square, 1)
             )
-            operators = merge(local_operators, vertical_operators)
+            I_horizontal = CartesianIndex(0, 1)
+            horizontal_operators = Dict(
+                (I, I + I_horizontal) => O₂ for I in eachindex(IndexCartesian(), square) if I[2] < size(square, 1)
+            )
+            operators = merge(local_operators, vertical_operators, horizontal_operators)
             H4 = FiniteMPOHamiltonian(grid, operators)
 
             @test H4 ≈
                 FiniteMPOHamiltonian(grid, local_operators) +
-                FiniteMPOHamiltonian(grid, vertical_operators)
+                FiniteMPOHamiltonian(grid, vertical_operators) +
+                FiniteMPOHamiltonian(grid, horizontal_operators) atol = 1.0e-4
+
+            H5 = changebonds(H4 / 3 + 2H4 / 3, SvdCut(; trscheme = truncbelow(1.0e-12)))
+            psi = FiniteMPS(physicalspace(H5), V ⊕ oneunit(V))
+            @test expectation_value(psi, H4) ≈ expectation_value(psi, H5)
         end
     end