diff --git a/Project.toml b/Project.toml
index 2a2d97ebc..61f60f476 100644
--- a/Project.toml
+++ b/Project.toml
@@ -10,6 +10,7 @@ projects = ["test", "docs", "examples"]
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
 BlockTensorKit = "5f87ffc2-9cf1-4a46-8172-465d160bd8cd"
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
+Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 HalfIntegers = "f0d1745a-41c9-11e9-1dd9-e5d34d218721"
 KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
@@ -38,6 +39,7 @@ Accessors = "0.1"
 Adapt = "4"
 BlockTensorKit = "0.3.14"
 Compat = "3.47, 4.10"
+Dictionaries = "0.4.6"
 DocStringExtensions = "0.9.3"
 HalfIntegers = "1.6.0"
 KrylovKit = "0.8.3, 0.9.2, 0.10"
diff --git a/src/MPSKit.jl b/src/MPSKit.jl
index b7387f3da..3c8da7226 100644
--- a/src/MPSKit.jl
+++ b/src/MPSKit.jl
@@ -49,6 +49,7 @@ export entropy, entanglement_spectrum
 export open_boundary_conditions, periodic_boundary_conditions
 export entanglementplot, transferplot
 export r_LL, l_LL, r_RR, l_RR, r_RL, r_LR, l_RL, l_LR # TODO: rename
+export mpo_compression
 
 # unexported
 using Compat: @compat
@@ -79,13 +80,14 @@ using Accessors
 using HalfIntegers
 using DocStringExtensions
 
-using LinearAlgebra: diag, Diagonal
+using LinearAlgebra: diag, Diagonal, diagm
 using LinearAlgebra: LinearAlgebra
 using Random
 using Base: @kwdef, @propagate_inbounds
 using LoggingExtras
 using OhMyThreads
 using TimerOutputs: TimerOutput, @timeit, timeit, reset_timer!, disable_timer!, enable_timer!
+using Dictionaries # TODO: change mpocompression to not use Dictionaries
 
 # Includes
 # --------
@@ -130,6 +132,7 @@ include("operators/projection.jl")
 include("operators/timedependence.jl")
 include("operators/multipliedoperator.jl")
 include("operators/lazysum.jl")
+include("operators/mpocompression.jl")
 
 include("transfermatrix/transfermatrix.jl")
 include("transfermatrix/transfer.jl")
diff --git a/src/operators/mpocompression.jl b/src/operators/mpocompression.jl
new file mode 100644
index 000000000..9691abf3f
--- /dev/null
+++ b/src/operators/mpocompression.jl
@@ -0,0 +1,634 @@
+"""
+    mpo_compression(H::FiniteMPOHamiltonian; η::Number=10^-8) -> Hᵪ, Rs
+    mpo_compression(H::InfiniteMPOHamiltonian; η::Number=10^-8) -> Qs, Ps
+
+Returns the compressed version of the FiniteMPOHamiltonian or InfiniteMPOHamiltonian,
+as described in: https://arxiv.org/pdf/1909.06341.
+
+### Arguments
+- `H`: Hamiltonian to be compressed. Either FiniteMPOHamiltonian or InfiniteMPOHamiltonian.
+
+### Keyword arguments
+- `η::Number=10^-8`: precision of the compression. The values smaller than η in the SVD 
+  decomposition during compression will be disregarded. 
+
+### Returns
+- mpo_compression(H::FiniteMPOHamiltonian; η::Number=10^-8)
+    `Hᵪ::FiniteMPOHamiltonian`: compressed Hamiltonian in the left canonical form.
+    `Rs::Vector{BlockTensorMap}`: matrices that can be applied to change between left and 
+      right canonical forms of Hᵪ.
+- mpo_compression(H::InfiniteMPOHamiltonian; η::Number=10^-8)
+    `Qs::InfiniteMPOHamiltonian`: compressed Hamiltonian in the left canonical form. 
+    `Ps::InfiniteMPOHamiltonian`: compressed Hamiltonian in the right canonical form.
+"""
+
+function mpo_compression end
+
+function mpo_compression(H::FiniteMPOHamiltonian, η::Number=10^-8)
+    t = @elapsed begin
+        # make sure that the MPO have at least three sites
+        N = length(H)
+        @assert N ≥ 3
+        # change MPO from N x 1 x 1 x M blocks to 3 x 1 x 1 x 3 block structure
+        TT = TensorMap{scalartype(H), spacetype(H[1]), 2, 2, Vector{scalartype(H)}}
+        Hnew = Vector{SparseBlockTensorMap{TT, scalartype(H), spacetype(H[1]), 2, 2, 4}}(
+                    undef, N)
+        for n in 1:N
+            Hnew[n] = SparseBlockTensorMap(reduce_blocks_mpo(H[n]))
+        end
+
+        Qs, Ls = right_canonical_mpo_finite(Hnew)
+        TT = AbstractTensorMap{scalartype(H), spacetype(H), 2, 2}
+        Cs = Vector{SparseBlockTensorMap{TT, scalartype(H), spacetype(H), 2, 2, 4}}(undef, N)
+        Cs[1] = Qs[1]
+        Rs = Vector{AbstractBlockTensorMap}(undef, N-2)
+        R = id(codomain(Qs[2])[1])
+        for n in 2:N-1
+            @tensor Qnew[a,i;j,b] := R[a,c] * Qs[n][c,i;j,b]
+            Q, R = qr_block_respecting(Qnew)
+            # For now it is safe to assume that R will be 3x3 blocked BlockTensorMap,
+            # since it went through right_canonical_mpo_finite first
+            M, R′ = decompose_R(R)
+            U, S, V′ = block_respecting_svd(M, η)
+            @tensor C[a,i;j,b] := Q[a,i;j,c] * U[c,b]
+            @tensor R[a;b] := S[a,c] * V′[c,d] * R′[d,b]
+            Cs[n] = C
+            Rs[n-1] = R
+        end
+        @tensor Cnew[a,i;j,b] := R[a,c] * Qs[N][c,i;j,b]
+        Cs[N] = Cnew
+        # Last step: Make new Cs into FiniteMPOHamiltonian
+        # TA = AbstractTensorMap{scalartype(H), spacetype(H), 2, 2}
+        # TB = AbstractTensorMap{scalartype(H), spacetype(H), 2, 1}
+        # TC = AbstractTensorMap{scalartype(H), spacetype(H), 1, 2}
+        # TD = AbstractTensorMap{scalartype(H), spacetype(H), 1, 1}
+        Hᵪtype = JordanMPOTensor{scalartype(H), spacetype(H), Vector{scalartype(H)}}
+        Hᵪ = Vector{Hᵪtype}(undef, N)
+        for n in 1:N
+            W = SparseBlockTensorMap(Cs[n])
+            A = W[2:(end-1), 1, 1, 2:(end-1)]
+            B = removeunit(W[2:(end-1), 1, 1, end], 4)
+            C = removeunit(W[1,1,1,2:(end-1)], 1)
+            D = removeunit(removeunit(W[1,1,1,end:end], 4), 1)
+            W = JordanMPOTensor(space(W), A, B, C, D)
+            Hᵪ[n] = W
+        end
+        Hᵪ = FiniteMPOHamiltonian(Hᵪ)
+    end
+    # Print reduction data
+    tot_dim_original = total_virt_dimension(H)
+    tot_dim_compressed = total_virt_dimension(Hᵪ)
+    printstyled("┌───────────────────────────────────────\n", color=:cyan)
+    printstyled("| MPO Compression performed succesfully\n", color=:cyan)
+    printstyled("│ ⏱  Time:   ", color=:cyan)
+    printstyled("$(round(t, digits=1)) s\n", color=:yellow)
+    printstyled("| Initial total virtual dimensions:    ", color=:cyan)
+    printstyled("$(tot_dim_original)\n", color=:yellow)
+    printstyled("| Final total virtual dimensions:    ", color=:cyan)
+    printstyled("$(tot_dim_compressed)\n", color=:yellow)
+    printstyled("└───────────────────────────────────────\n")
+
+    return Hᵪ, Rs
+end
+
+function mpo_compression(H::InfiniteMPOHamiltonian, η::Number=10^-8)
+    tot_dim_original = total_virt_dimension(H)
+    t = @elapsed begin
+        N = length(H)   # Number of sites
+        # change MPO from N x 1 x 1 x M blocks to 3 x 1 x 1 x 3 block structure
+        TT = TensorMap{scalartype(H), spacetype(H[1]), 2, 2, Vector{scalartype(H)}}
+        Hnew = Vector{SparseBlockTensorMap{TT, scalartype(H), spacetype(H[1]), 2, 2, 4}}(
+                    undef, N)
+        for n in 1:N
+            Hnew[n] = SparseBlockTensorMap(reduce_blocks_mpo(H[n]))
+        end
+
+        HL, _ = left_canonical_mpo_infinite_iter_msites(Hnew)
+        _, HR = right_canonical_mpo_infinite_iter_msites(HL)
+        HL, Cs = left_canonical_mpo_infinite_iter_msites(HR)
+        Us = Vector{BlockTensorMap}(undef, N)
+        Ss = Vector{BlockTensorMap}(undef, N)
+        Vs = Vector{BlockTensorMap}(undef, N)
+        for n = 1:N
+            Us[n], Ss[n], Vs[n] = block_respecting_svd(Cs[n], η)
+        end
+        Qs = Vector{BlockTensorMap}(undef, N)
+        Ps = Vector{BlockTensorMap}(undef, N)
+        for n = 1:N
+            # Left canonical version of compressed iMPO
+            @tensor Q[a,i;j,b] := Us[mod1(n-1,N)]'[a,c] * HL[n][c,i;j,d] * Us[n][d,b]
+            # Right canonical version of compressed iMPO
+            @tensor P[a,i;j,b] := Vs[mod1(n-1,N)][a,c] * HR[n][c,i;j,d] * Vs[n]'[d,b]
+            Qs[n] = Q
+            Ps[n] = P
+        end
+
+        Q̂s = create_impoham_from_mpos_msites(Qs)
+        P̂s = create_impoham_from_mpos_msites(Ps)
+    end
+
+    # Print reduction data
+    tot_dim_compressed = total_virt_dimension(Q̂s)
+    printstyled("┌───────────────────────────────────────\n", color=:cyan)
+    printstyled("| MPO Compression performed succesfully\n", color=:cyan)
+    printstyled("│ ⏱  Time:   ", color=:cyan)
+    printstyled("$(round(t, digits=1)) s\n", color=:yellow)
+    printstyled("| Initial total virtual dimensions:    ", color=:cyan)
+    printstyled("$(tot_dim_original)\n", color=:yellow)
+    printstyled("| Final total virtual dimensions:    ", color=:cyan)
+    printstyled("$(tot_dim_compressed)\n", color=:yellow)
+    printstyled("└───────────────────────────────────────\n")
+
+    return Q̂s, P̂s
+end
+
+# utility
+# -------
+function trace_single_block(d::AbstractTensorMap)
+    return @tensor trd[a;b] := d[a,i;i,b]
+end
+
+function change_value_at_fusiontree!(T::AbstractTensorMap, i::Integer, v::Number)
+    # change value of the fusiontree at index i of a TensorMap T with new value v
+    f = fusiontrees(T)
+    if f isa Vector
+        T[f[i]...] .= v
+    elseif f isa Indices
+        T[f.values[i]...] .= v
+    else
+        throw(ArgumentError("fusiontrees of T is neither a Vector or a Indices. It is
+                            a $(typeof(f))"))
+    end
+end
+
+function total_virt_dimension(H::FiniteMPOHamiltonian)
+    # calculate the sum of bond dimensions for entire H MPO
+    N = length(H)
+    tvd = sum([dims(H[i])[4] for i in 1:N-1])
+    return tvd
+end
+
+function total_virt_dimension(H::InfiniteMPOHamiltonian)
+    # calculate the sum of bond dimensions for entire H MPO
+    N = length(H)
+    tvd = sum([dims(H[i])[4] for i in 1:N])
+    return tvd
+end
+
+# functions
+# ---------
+function qr_block_respecting(W::AbstractBlockTensorMap)
+    # Block respecting qr (qr applied only to the upper-left part of the matrix)
+    bl_st = size(W)
+    dim_st = dims(W)
+
+    # sanity checks
+    @assert length(bl_st) == length(dim_st) == 4
+    blN = bl_st[1]
+    blM = bl_st[4]
+
+    # First, for each nonzero fusiontree from blocks in the first row, make them
+    # traceless
+    id1 = id(codomain(W)[2])
+
+    t = zeros(ComplexF64, domain(W)[2][1] ← domain(W)[2][2:end-1])
+    for block in 2:blM-1
+        d = W[1,1,1,block]
+        trd = trace_single_block(d)
+        @tensor W[1,1,1,block][a,i;j,b] = d[a,i;j,b] - trd[a;b] * id1[i;j]
+        # Here goes the creation of matrix R' with a t' tensor
+        t[block-1] = trd / dim_st[2]
+    end
+    # Perform QR of the 2:blM-1 block
+    V = W[1:blN-1, :, :, 2:blM-1]
+    V = permute(V, ((1,2,3), (4,)))
+    
+    V = convert(SparseBlockTensorMap{TensorMap{eltype(V).parameters...}, 
+                eltype(storagetype(V)), spacetype(V), 3, 1, 4,}, V)
+    Q, R = qr_compact(V)
+
+    Q = permute(Q, ((1,2),(3,4)))
+
+    # Put Q into larger Q̂
+    Q̂ = zeros(scalartype(W), codomain(W) ← domain(W)[1] ⊗ (domain(W)[2][1] ⊞ domain(Q)[2] ⊞ 
+                domain(W)[2][end]))
+    Q̂[1,1,1,1] = Q̂[end,1,1,end] = W[1,1,1,1]
+    Q̂[1:end-1,1,1,2:end-1] = Q
+    Q̂[1:end-1,1,1,end] = W[1:end-1,1,1,end]
+
+    # Put R and t into larger R̂
+    R̂ = zeros(scalartype(W), domain(W)[2][1] ⊞ codomain(R)[1] ⊞ domain(W)[2][end] ←
+              domain(W)[2][1] ⊞ domain(R)[1] ⊞ domain(W)[2][end])
+    R̂[2:end-1,2:end-1] = R
+    R̂[1,2:end-1] = t
+    change_value_at_fusiontree!(R̂[1,1],1,1)
+    change_value_at_fusiontree!(R̂[end,end],1,1)
+    
+    return Q̂, R̂
+end
+
+function lq_block_respecting(W::AbstractBlockTensorMap)
+    # Block respecting lq (lq applied only to the bottom-right part of the matrix)
+    bl_st = size(W)
+    dim_st = dims(W)
+
+    # sanity checks
+    @assert length(bl_st) == length(dim_st) == 4
+    blN = bl_st[1]
+    blM = bl_st[4]
+    # First, for each nonzero fusiontree from blocks in the last column, make them
+    # traceless
+    id1 = id(codomain(W)[2])
+
+    t = zeros(ComplexF64, codomain(W)[1][2:end-1] ← codomain(W)[1][end])
+    for block in 2:blN-1
+        d = W[block,1,1,end]
+        trd = trace_single_block(d)
+        @tensor W[block,1,1,end][a,i;j,b] = d[a,i;j,b] - trd[a;b] * id1[i;j]
+        # Here goes the creation of matrix R' with a t' tensor
+        t[block-1] = trd / dim_st[2]
+    end
+    # Perform QR of the 2:blM-1 block (with virtual indices transposed)
+    V = W[2:blN-1, :, :, 2:blM]
+    V = permute(V, ((1,), (2,3,4)))
+
+    V = convert(SparseBlockTensorMap{TensorMap{eltype(V).parameters...}, 
+                eltype(storagetype(V)), spacetype(V), 1, 3, 4,}, V)
+    # to perform QR
+    L, Q = lq_compact(V)
+    # transpose vitrual indices back (now to reconstruct initial matrix
+    # one has to apply R on the left)
+    Q = permute(Q, ((1,2),(3,4)))
+
+    # Put Q into larger Q̂
+    Q̂ = zeros(scalartype(W), (codomain(W)[1][1] ⊞ codomain(Q)[1] ⊞ codomain(W)[1][end]) ⊗ 
+                codomain(W)[2] ← domain(W))
+    Q̂[1,1,1,1] = Q̂[end,1,1,end] = W[1,1,1,1]
+    Q̂[2:end-1,1,1,2:end] = Q
+    Q̂[1,1,1,2:end] = W[1,1,1,2:end]
+
+    # Put L and t into larger L̂
+    L̂ = zeros(scalartype(W), domain(W)[2][1] ⊞ codomain(L)[1] ⊞ domain(W)[2][end] ←
+              domain(W)[2][1] ⊞ domain(L)[1] ⊞ domain(W)[2][end])
+    L̂[2:end-1,2:end-1] = L
+    L̂[2:end-1,end] = t
+    change_value_at_fusiontree!(L̂[1,1],1,1)
+    change_value_at_fusiontree!(L̂[end,end],1,1)
+
+    return L̂, Q̂
+end
+
+function left_canonical_mpo_finite(
+        H::Union{FiniteMPOHamiltonian, V}
+    ) where {T<:SparseBlockTensorMap{<:Any, <:Any, <:Any, 2, 2, 4}, V<:Vector{T}}
+    # change finite Hamiltonian MPO into its Left canonical form. The 'H' is a
+    # FiniteMPOHamiltonian from MPSKit (l and r boundary conditions are already
+    # part of the left-most and right-most MPOs, they can not be compressed, and
+    # they can be excluded from left canonicalisation algorithm).
+
+    # make sure that the MPO have at least three sites
+    N = length(H)
+    @assert N ≥ 3
+
+    Q, R = qr_block_respecting(H[2])
+    TT = AbstractTensorMap{scalartype(H), spacetype(H[1]), 2, 2}
+    Qs = Vector{SparseBlockTensorMap{TT, scalartype(H), spacetype(H[1]), 2, 2, 4}}(undef, N)
+    Qs[1] = H[1]
+    Qs[2] = Q
+    Rs = Vector{AbstractBlockTensorMap}(undef, N-2)
+    Rs[1] = R
+    for n in 3:N-1
+        @tensor Wnew[a,i;j,b] := R[a,c] * H[n][c,i;j,b]
+        Q, R = qr_block_respecting(Wnew)
+        Qs[n] = Q
+        Rs[n-1] = R
+    end
+    @tensor Wlast[a,i;j,b] := R[a,c] * H[N][c,i;j,b]
+    Qs[N] = Wlast
+    return Qs, Rs
+end
+
+function right_canonical_mpo_finite(
+        H::Union{FiniteMPOHamiltonian, V}
+    ) where {T<:SparseBlockTensorMap{<:Any, <:Any, <:Any, 2, 2, 4}, V<:Vector{T}}
+    # change finite Hamiltonian MPO into its Right canonical form. The 'H' is a
+    # FiniteMPOHamiltonian from MPSKit (l and r boundary conditions are already
+    # part of the left-most and right-most MPOs, they can not be compressed, and
+    # they can be excluded from right canonicalisation algorithm).
+
+    # make sure that the MPO have at least three sites
+    N = length(H)
+    @assert N ≥ 3
+
+    L, Q = lq_block_respecting(H[end-1])
+    TT = AbstractTensorMap{scalartype(H), spacetype(H[1]), 2, 2}
+    Qs = Vector{SparseBlockTensorMap{TT, scalartype(H), spacetype(H[1]), 2, 2, 4}}(undef, N)
+    Qs[end] = H[end]
+    Qs[end-1] = Q
+    Ls = Vector{AbstractBlockTensorMap}(undef, N-2)
+    Ls[end] = L
+    for n in N-2:-1:2
+        @tensor Wnew[a,i;j,b] := H[n][a,i;j,c] * L[c,b]
+        L, Q = lq_block_respecting(Wnew)
+        Qs[n] = Q
+        Ls[n-1] = L
+    end
+    @tensor Wlast[a,i;j,b] := H[1][a,i;j,c] * L[c,b]
+    Qs[1] = Wlast
+    return Qs, Ls
+end
+
+function decompose_R(R::AbstractBlockTensorMap)
+    M = zeros(ComplexF64, codomain(R) ← domain(R))
+    R′ = zeros(ComplexF64, domain(R) ← domain(R))
+    change_value_at_fusiontree!(M[1,1], 1, 1)
+    change_value_at_fusiontree!(M[end,end], 1, 1)
+    change_value_at_fusiontree!(R′[1,1], 1, 1)
+    change_value_at_fusiontree!(R′[end,end], 1, 1)
+    M[2:end-1,2:end-1] = R[2:end-1,2:end-1]
+    R′[1,2:end-1] = R[1,2:end-1]
+    # R′[2:end-1,2:end-1] = id(codomain(R[2,2]))
+    # Make diagonal matrix even with codomain and domain being different spaces
+    # Right now I'm assuming that each of the fusiontree is a square block
+    for (s, f) in fusiontrees(R′[2:end-1,2:end-1])
+        if s.uncoupled[1] == f.uncoupled[1]
+            dim = size(R′[2:end-1,2:end-1][s,f])[1]
+            R′[2:end-1,2:end-1][s,f] .= diagm(ones(dim))
+        end 
+    end
+    return M, R′
+end
+
+function block_respecting_svd(M::AbstractBlockTensorMap, η::Number; perform_truncation=true)
+    # Perform block-respecting svd on matrix M, with precision η
+    # TODO Change it so that it works when sectors are numbered by indices rather than 
+    #      enumerated
+    M2 = M[2:end-1,2:end-1]
+    U, S, V′ = svd_compact(M2)
+
+    num_sectors = length(blocksectors(S))
+    new_dims = Vector{Int32}(undef, num_sectors)
+
+    # Old code - not working if blocksectors(S) returned Indices rather than Vector
+    # for n in 1:num_sectors
+    #     csec = blocksectors(S)[n]
+    #     new_dims[n] = blockdim(S.domain, csec)
+    #     if perform_truncation
+    #         for diag_el in block(S, csec).diag
+    #             if diag_el ≤ η
+    #                 new_dims[n] -= 1
+    #             end
+    #         end
+    #     end
+    # end
+    
+    # New code - works for both cases
+    for (n, csec) in enumerate(blocksectors(S))
+        new_dims[n] = blockdim(S.domain, csec)
+        if perform_truncation
+            for diag_el in block(S, csec).diag
+                if diag_el ≤ η
+                    new_dims[n] -= 1
+                end
+            end
+        end
+    end
+
+    fdims = Vector{Pair{sectortype(S), Int32}}()
+    # Old code 
+    # for n in 1:num_sectors
+    #     if new_dims[n] > 0
+    #         push!(fdims, blocksectors(S)[n] => new_dims[n])
+    #     end
+    # end
+
+    # New code - works for both cases
+    for (n, csec) in enumerate(blocksectors(S))
+        if new_dims[n] > 0
+            push!(fdims, csec => new_dims[n])
+        end
+    end
+    red_space = spacetype(S)(fdims...)
+    S_new = DiagonalTensorMap(undef, red_space)
+    U_new = zeros(ComplexF64, codomain(U) ← red_space)
+    V′_new = zeros(ComplexF64, red_space ← domain(V′))
+
+    #num_sectors2 = length(blocksectors(S_new))
+
+    # Old code
+    # for n in 1:num_sectors2
+    #     csec = blocksectors(S_new)[n]
+    #     cl = length(block(S_new, csec).diag)
+    #     block(S_new, csec).diag[1:end] = block(S, csec).diag[1:cl]
+    #     block(U_new, csec)[:,:] = block(U, csec)[:,1:cl]
+    #     block(V′_new, csec)[:,:] = block(V′, csec)[1:cl,:]
+    # end
+
+    # New code - works for both cases
+    for csec in blocksectors(S_new)
+        cl = length(block(S_new, csec).diag)
+        block(S_new, csec).diag[1:end] = block(S, csec).diag[1:cl]
+        block(U_new, csec)[:,:] = block(U, csec)[:,1:cl]
+        block(V′_new, csec)[:,:] = block(V′, csec)[1:cl,:]
+    end
+
+    Û = zeros(ComplexF64, codomain(M) ← domain(M)[1][1] ⊞ domain(U_new)[1] ⊞ domain(M)[1][end])
+    Ŝ = zeros(ComplexF64, codomain(M)[1][1] ⊞ codomain(S_new)[1] ⊞ codomain(M)[1][end] ← 
+              domain(M)[1][1] ⊞ domain(S_new)[1] ⊞ domain(M)[1][end])
+    V̂′ = zeros(ComplexF64, codomain(M)[1][1] ⊞ codomain(V′_new)[1] ⊞ codomain(M)[1][end] ← domain(M))
+
+    change_value_at_fusiontree!(Û[1,1], 1, 1)
+    change_value_at_fusiontree!(Û[end,end], 1, 1)
+    change_value_at_fusiontree!(Ŝ[1,1], 1, 1)
+    change_value_at_fusiontree!(Ŝ[end,end], 1, 1)
+    change_value_at_fusiontree!(V̂′[1,1], 1, 1)
+    change_value_at_fusiontree!(V̂′[end,end], 1, 1)
+
+    Û[2:end-1,2:end-1] = U_new
+    Ŝ[2,2] = S_new
+    V̂′[2:end-1,2:end-1] = V′_new
+
+    return Û, Ŝ, V̂′
+end
+
+function reduce_blocks_mpo(W::AbstractBlockTensorMap)
+    # Function that copies BlockTensorMap W from N x 1 x 1 x M blocked matrix to a 
+    # 3 x 1 x 1 x 3 blocked matrix, with the middle block being joined dimensions of N-2, 
+    # M-2 blocks from the original matrix
+    # TODO Make exception if either N or M is smaller than 3
+
+    # If the matrix is already blocked properly, return it
+    if size(W)[1] ≤ 3 && size(W)[4] ≤ 3
+        return W
+    end
+    
+    old_cod = codomain(W)[1][2:end-1]
+    old_dom = domain(W)[2][2:end-1]
+
+    new_cod = spacetype(W)([sec => blockdim(old_cod, sec) for sec in sectors(old_cod)]...)
+    new_dom = spacetype(W)([sec => blockdim(old_dom, sec) for sec in sectors(old_dom)]...)
+
+    # Create new matrix with zeros
+    if size(W)[1] ≤ 3
+        Vcod = codomain(W)
+    else
+        Vcod = (codomain(W)[1][1] ⊞ new_cod ⊞ codomain(W)[1][end]) ⊗ codomain(W)[2]
+    end
+    if size(W)[4] ≤ 3
+        Vdom = domain(W)
+    else
+        Vdom = domain(W)[1] ⊗ (domain(W)[2][1] ⊞ new_dom ⊞ domain(W)[2][end])
+    end
+    W2 = zeros(scalartype(W), Vcod ← Vdom)
+
+    # Corners can be copied trivially
+    W2[1,1,1,1] = W[1,1,1,1]
+    W2[end,1,1,end] = W[end,1,1,end]
+    W2[1,1,1,end] = W[1,1,1,end]
+    W2[end,1,1,1] = W[end,1,1,1]
+    # Treat the cases with smaller codomain or domain separately
+    if size(W)[1] ≤ 2
+        # Middle upperblock
+        for (s,f) in fusiontrees(W[1,1,1,2:end-1])
+            W2[1,1,1,2][s,f] = W[1,1,1,2:end-1][s,f]
+        end
+        # Middle bottomblock (if size(W)[1] is 1, this is the same block as in 
+        # previous loop)
+        for (s,f) in fusiontrees(W[end,1,1,2:end-1])
+            W2[end,1,1,2][s,f] = W[end,1,1,2:end-1][s,f]
+        end
+    elseif size(W)[4] ≤ 2
+        # Middle leftblock
+        for (s,f) in fusiontrees(W[2:end-1,1,1,1])
+            W2[2,1,1,1][s,f] = W[2:end-1,1,1,1][s,f]
+        end
+        # Middle rightblock (if size(W)[4] is 1, this is the same block as in previous loop)
+        for (s,f) in fusiontrees(W[2:end-1,1,1,end])
+            W2[2,1,1,end][s,f] = W[2:end-1,1,1,end][s,f]
+        end
+    else
+        # All other blocks can be dealt with fusiontrees
+        # Middle block (most important)
+        for (s,f) in fusiontrees(W[2:end-1,1,1,2:end-1])
+            W2[2,1,1,2][s,f] = W[2:end-1,1,1,2:end-1][s,f]
+        end
+        # Upper block
+        for (s,f) in fusiontrees(W[1,1,1,2:end-1])
+            W2[1,1,1,2][s,f] = W[1,1,1,2:end-1][s,f]
+        end
+        # Bottom block
+        for (s,f) in fusiontrees(W[end,1,1,2:end-1])
+            W2[3,1,1,2][s,f] = W[end,1,1,2:end-1][s,f]
+        end
+        # Left block
+        for (s,f) in fusiontrees(W[2:end-1,1,1,1])
+            W2[2,1,1,1][s,f] = W[2:end-1,1,1,1][s,f]
+        end
+        # Right block
+        for (s,f) in fusiontrees(W[2:end-1,1,1,end])
+            W2[2,1,1,3][s,f] = W[2:end-1,1,1,end][s,f]
+        end
+    end
+    return W2
+end
+
+function check_how_far_from_id(R::AbstractBlockTensorMap)
+    # Function checking how far the diagonal elements of R are from either
+    # +1 or -1. Assume R is 3x3 and [1,1] and [3,3] blocks are always id
+    tdev = 0
+    for (s,f) in fusiontrees(R[2,2])
+        if s.uncoupled[1] == f.uncoupled[1]
+            for d in size(R[2,2][s,f])[1]
+                tdev += min(abs(R[2,2][s,f][d,d] - 1), abs(R[2,2][s,f][d,d] + 1))
+            end
+        end
+    end
+    return tdev
+end
+
+function create_impoham_from_mpos_msites(Qs::Vector{BlockTensorMap})
+    # Given the Vector{BlockTensorMap} Qs, return the multisite InfiniteMPOHamiltonian 
+    N = length(Qs)
+    sctype = scalartype(Qs[1])
+    sptype = spacetype(Qs[1])
+    # TA = AbstractTensorMap{sctype, sptype, 2, 2}
+    # TB = AbstractTensorMap{sctype, sptype, 2, 1}
+    # TC = AbstractTensorMap{sctype, sptype, 1, 2}
+    # TD = AbstractTensorMap{sctype, sptype, 1, 1}
+    CType = JordanMPOTensor{sctype, sptype, Vector{sctype}}
+    Cs = Vector{CType}(undef, N)
+    for n in 1:N
+        W = SparseBlockTensorMap(Qs[n])
+        A = W[2:(end-1), 1, 1, 2:(end-1)]
+        B = removeunit(W[2:(end-1), 1, 1, end], 4)
+        C = removeunit(W[1,1,1,2:(end-1)], 1)
+        D = removeunit(removeunit(W[1,1,1,end:end], 4), 1)
+        W = JordanMPOTensor(space(W), A, B, C, D)
+        Cs[n] = W
+    end
+    Hnew = InfiniteMPOHamiltonian(Cs)
+    return Hnew
+end
+
+function left_canonical_mpo_infinite_iter_msites(
+        H::Union{InfiniteMPOHamiltonian, V},
+        η=10^-10
+    ) where {T<:SparseBlockTensorMap{<:Any, <:Any, <:Any, 2, 2, 4}, V<:Vector{T}}
+    # Return the left_canonical form of the iMPO, H can have multiple sites in the unit cell
+    N = length(H)       # Number of sites
+
+    εs = fill(Inf, N)
+    Ls = fill(id(ComplexF64, codomain(H[1])[1]), N)
+    Rs = fill(id(ComplexF64, codomain(H[1])[1]), N)
+    Q = similar(H[1])
+    printstyled("┌─── Started iterative left orthogonalization ────\n", color=:cyan)
+    while sum(εs)/N > η
+        # First store the R matrices and change them H into Qs
+        for n in 1:N
+            Q, R = qr_block_respecting(H[n])
+            H[n] = Q
+            Rs[n] = R
+            Ls[n] = R * Ls[n]
+            εs[n] = check_how_far_from_id(R)
+        end
+        # Then assign R[(n+1)//N] * Q[n] to H[n]
+        for n in 1:N
+            @tensor H[n][a,i;j,b] = Rs[mod1(n-1,N)][a;c] * H[n][c,i;j,b]
+        end
+        printstyled("| total err = $(sum(εs)/N)\n", color=:white)
+    end
+    printstyled("└─── Left orthogonalization finished correctly ──\n", color=:cyan)
+    return H, Ls
+end
+
+function right_canonical_mpo_infinite_iter_msites(
+        H::Union{InfiniteMPOHamiltonian, V},
+        η=10^-10
+    ) where {T<:SparseBlockTensorMap{<:Any, <:Any, <:Any, 2, 2, 4}, V<:Vector{T}}
+    # Return the right_canonical form of the iMPO, H can have multiple sites in the unit 
+    # cell
+    N = length(H)       # Number of sites
+    
+    εs = fill(Inf, N)
+    Rs = fill(id(ComplexF64, codomain(H[1])[1]), N)
+    Ls = fill(id(ComplexF64, codomain(H[1])[1]), N)
+    Q = similar(H[N])
+    printstyled("┌─── Started iterative right orthogonalization ────\n", color=:cyan)
+    while sum(εs)/N > η
+        # First store the R matrices and change them H into Qs
+        for n in N:-1:1
+            L, Q = lq_block_respecting(H[n])
+            H[n] = Q
+            Ls[n] = L
+            Rs[n] = Rs[n] * L
+            εs[n] = check_how_far_from_id(L)
+        end
+        # Then assign R[(n+1)//N] * Q[n] to H[n]
+        for n in N:-1:1
+            @tensor H[n][a,i;j,b] = H[n][a,i;j,c] * Ls[mod1(n+1, N)][c,b]
+        end
+        printstyled("| total err = $(sum(εs)/N)\n", color=:white)
+    end
+    printstyled("└─── Right orthogonalization finished correctly ──\n", color=:cyan)
+    return Rs, H
+end
\ No newline at end of file
diff --git a/test/operators/finite_mpo_compression.jl b/test/operators/finite_mpo_compression.jl
new file mode 100644
index 000000000..07c259cf6
--- /dev/null
+++ b/test/operators/finite_mpo_compression.jl
@@ -0,0 +1,375 @@
+using Test, Random
+using TensorKit
+using MPSKit
+using MPSKitModels
+
+function max_dim_ham(H)
+    N = length(H)
+    max_dim = 1
+    for i in 1:N
+        if max_dim < dims(H[i])[4]
+            max_dim = dims(H[i])[4]
+        end
+    end
+    return max_dim
+end
+
+# Utility functions to create Hamiltonians for testing
+
+function create_long_range_ising_symmetries_finite(M, k)
+    chain = fill(Z2Space(0=>1, 1=>1), M)
+    ZZ = S_zz(Z2Irrep)
+    X = S_x(Z2Irrep)
+    single_site_operators = [i => X for i in 1:M]
+    two_site_operators = [(i,i+j) => ZZ for i in 1:M for j in 1:k if i+j ≤ M]
+    H = FiniteMPOHamiltonian(chain, single_site_operators...,
+                                    two_site_operators...)
+    return H
+end
+
+function create_long_range_ising_symmetries_random_finite(M, k, σⱼ, σₕ)
+    chain = fill(Z2Space(0=>1, 1=>1), M)
+    ZZ = S_zz(Z2Irrep)
+    X = S_x(Z2Irrep)
+    single_site_operators = [i => (σₕ * randn()) * X for i in 1:M]
+    two_site_operators = [(i,i+j) => (1 + σⱼ * randn()) * ZZ for i in 1:M for j in 1:k if i+j ≤ M]
+    H = FiniteMPOHamiltonian(chain, single_site_operators...,
+                                    two_site_operators...)
+    return H
+end
+
+function transverse_field_ising_fermion_parity_finite(M)
+    chain = FiniteChain(M)
+    return transverse_field_ising(ComplexF64, FermionParity, chain)
+end
+
+function transverse_field_ising_z2_finite(M)
+    chain = FiniteChain(M)
+    return transverse_field_ising(ComplexF64, Z2Irrep, chain)
+end
+
+function transverse_field_ising_trivial_finite(M)
+    chain = FiniteChain(M)
+    return transverse_field_ising(ComplexF64, Trivial, chain)
+end
+
+function kitaev_model_finite(M)
+    chain = FiniteChain(M)
+    return kitaev_model(ComplexF64, chain)
+end
+
+function heisenberg_XXX_trivial_finite(M)
+    chain = FiniteChain(M)
+    return heisenberg_XXX(Trivial, chain)
+end
+
+function heisenberg_XXX_U1_finite(M)
+    chain = FiniteChain(M)
+    return heisenberg_XXX(U1Irrep, chain)
+end
+
+function heisenberg_XXX_SU2_finite(M)
+    chain = FiniteChain(M)
+    return heisenberg_XXX(SU2Irrep, chain)
+end
+
+function bilinear_biquadratic_model_trivial_finite(M)
+    chain = FiniteChain(M)
+    return bilinear_biquadratic_model(ComplexF64, Trivial, chain)
+end
+
+function bilinear_biquadratic_model_U1_finite(M)
+    chain = FiniteChain(M)
+    return bilinear_biquadratic_model(ComplexF64, U1Irrep, chain)
+end
+
+function bilinear_biquadratic_model_SU2_finite(M)
+    chain = FiniteChain(M)
+    return bilinear_biquadratic_model(ComplexF64, SU2Irrep, chain)
+end
+
+function quantum_potts_trivial_finite(M)
+    chain = FiniteChain(M)
+    return quantum_potts(ComplexF64, Trivial, chain)
+end
+
+function quantum_potts_Z3_finite(M)
+    chain = FiniteChain(M)
+    return quantum_potts(ComplexF64, Z3Irrep, chain)
+end
+
+# Tests
+
+@testset "finite_mpo_compression_energy" begin
+    # test the finite compression for a long range Ising model
+    # check if the energy after the compression is the same as before
+    D = 20 # max bond dimensions
+    L = 10 # number of sites
+    range = 4 # cutoff for the maximum range between interactions
+    trunc_st = notrunc()
+
+    mps = FiniteMPS(L, Z2Space(0=>1, 1=>1), Z2Space(0=>D, 1=>D))
+    H = create_long_range_ising_symmetries_finite(L, range)
+    find_groundstate!(mps, H, DMRG2(; maxiter=100, trscheme=trunc_st))
+    E0 = expectation_value(mps, H)
+    #println("<mps|H|mps> = $real(E0)")
+
+    mps2 = FiniteMPS(L, Z2Space(0=>1, 1=>1), Z2Space(0=>D, 1=>D))
+    H2, Rs = mpo_compression(H)
+    find_groundstate!(mps2, H2, DMRG2(; maxiter=100, trscheme=trunc_st))
+    E1 = expectation_value(mps2, H2)
+end
+
+@testset "finite_mpo_compression_transverse_ising_model_energy" begin
+    M = 10
+    D = 10 # Max bond dimension
+    trunc_st = trunctol()
+
+    mps1 = FiniteMPS(M, ℂ^2, ℂ^D)
+    mps2 = FiniteMPS(M, Z2Space(0=>1, 1=>1), Z2Space(0=>D, 1=>D))
+    mps3 = FiniteMPS(M, Vect[FermionParity](0=>1, 1=>1), Vect[FermionParity](0=>D, 1=>D))
+
+    H1 = transverse_field_ising_trivial_finite(M)
+    H2 = transverse_field_ising_z2_finite(M)
+    H3 = transverse_field_ising_fermion_parity_finite(M)
+
+    # Find groundstate of non-compressed Hamiltonian
+    find_groundstate!(mps1, H1, DMRG2(; maxiter=100, trscheme=trunc_st))
+    E0_1 = expectation_value(mps1, H1)
+    find_groundstate!(mps2, H2, DMRG2(; maxiter=100, trscheme=trunc_st))
+    E0_2 = expectation_value(mps2, H2)
+    find_groundstate!(mps3, H3, DMRG2(; maxiter=100, trscheme=trunc_st))
+    E0_3 = expectation_value(mps3, H3)
+
+    # Compress MPO
+    mps1c = FiniteMPS(M, ℂ^2, ℂ^D)
+    mps2c = FiniteMPS(M, Z2Space(0=>1, 1=>1), Z2Space(0=>D, 1=>D))
+    mps3c = FiniteMPS(M, Vect[FermionParity](0=>1, 1=>1), Vect[FermionParity](0=>D, 1=>D))
+
+    H1c, _ = mpo_compression(H1, 10^-10)
+    H2c, _ = mpo_compression(H2, 10^-10)
+    H3c, _ = mpo_compression(H3, 10^-10)
+
+    # Find groundstate of compressed Hamiltonian
+    find_groundstate!(mps1c, H1c, DMRG2(; maxiter=100, trscheme=trunc_st))
+    E0_1c = expectation_value(mps1c, H1c)
+    find_groundstate!(mps2c, H2c, DMRG2(; maxiter=100, trscheme=trunc_st))
+    E0_2c = expectation_value(mps2c, H2c)
+    find_groundstate!(mps3c, H3c, DMRG2(; maxiter=100, trscheme=trunc_st))
+    E0_3c = expectation_value(mps3c, H3c)
+
+    # Assert that the ground state energies are equal up to a precision, display the 
+    # compression
+    @assert E0_1 ≈ E0_1c
+    @assert E0_2 ≈ E0_2c
+    @assert E0_3 ≈ E0_3c
+
+    println("Max dim uncompressed H1: $(max_dim_ham(H1)) || Max dim compressed H1c: 
+            $(max_dim_ham(H1c))" )
+    println("Max dim uncompressed H2: $(max_dim_ham(H2)) || Max dim compressed H2c: 
+            $(max_dim_ham(H2c))" )
+    println("Max dim uncompressed H3: $(max_dim_ham(H3)) || Max dim compressed H3c: 
+            $(max_dim_ham(H3c))" )
+
+end
+
+@testset "finite_mpo_compression_kitaev_model_energy" begin
+    M = 10
+    D = 10 # Max bond dimension
+    trunc_st = trunctol()
+
+    mps1 = FiniteMPS(M, Vect[FermionParity](0=>1, 1=>1), Vect[FermionParity](0=>D, 1=>D))
+
+    H1 = kitaev_model_finite(M)
+
+    # Find groundstate of non-compressed Hamiltonian
+    find_groundstate!(mps1, H1, DMRG2(; maxiter=100, trscheme=trunc_st))
+    E0_1 = expectation_value(mps1, H1)
+
+    # Compress MPO
+    mps1c = FiniteMPS(M, Vect[FermionParity](0=>1, 1=>1), Vect[FermionParity](0=>D, 1=>D))
+
+    H1c, _ = mpo_compression(H1, 10^-10)
+
+    # Find groundstate of compressed Hamiltonian
+    find_groundstate!(mps1c, H1c, DMRG2(; maxiter=100, trscheme=trunc_st))
+    E0_1c = expectation_value(mps1c, H1c)
+
+    # Assert that the ground state energies are equal up to a precision, display the 
+    # compression
+    @assert E0_1 ≈ E0_1c
+
+    println("Max dim uncompressed H1: $(max_dim_ham(H1)) || Max dim compressed H1c: 
+            $(max_dim_ham(H1c))" )
+
+end
+
+@testset "finite_mpo_compression_heisenberg_XXX_energy" begin
+    M = 10
+    D = 6 # Max bond dimension
+    trunc_st = trunctol()
+
+    mps1 = FiniteMPS(M, ℂ^3, ℂ^D)
+    mps2 = FiniteMPS(M, U1Space(0=>1, 1=>1, -1=>1), U1Space(0=>D, 1=>D, -1=>D))
+    mps3 = FiniteMPS(M, SU2Space(1 => 1), SU2Space(1=>D))
+
+    H1 = heisenberg_XXX_trivial_finite(M)
+    H2 = heisenberg_XXX_U1_finite(M)
+    H3 = heisenberg_XXX_SU2_finite(M)
+
+    # Find groundstate of non-compressed Hamiltonian
+    find_groundstate!(mps1, H1, DMRG2(; maxiter=100, trscheme=trunc_st))
+    E0_1 = expectation_value(mps1, H1)
+    find_groundstate!(mps2, H2, DMRG2(; maxiter=100, trscheme=trunc_st))
+    E0_2 = expectation_value(mps2, H2)
+    find_groundstate!(mps3, H3, DMRG2(; maxiter=100, trscheme=trunc_st))
+    E0_3 = expectation_value(mps3, H3)
+
+    # Compress MPO
+    mps1c = FiniteMPS(M, ℂ^3, ℂ^D)
+    mps2c = FiniteMPS(M, U1Space(0=>1, 1=>1, -1=>1), U1Space(0=>D, 1=>D, -1=>D))
+    mps3c = FiniteMPS(M, SU2Space(1 => 1), SU2Space(1=>D))
+
+    H1c, _ = mpo_compression(H1, 10^-10)
+    H2c, _ = mpo_compression(H2, 10^-10)
+    H3c, _ = mpo_compression(H3, 10^-10)
+
+    # Find groundstate of compressed Hamiltonian
+    find_groundstate!(mps1c, H1c, DMRG2(; maxiter=100, trscheme=trunc_st))
+    E0_1c = expectation_value(mps1c, H1c)
+    find_groundstate!(mps2c, H2c, DMRG2(; maxiter=100, trscheme=trunc_st))
+    E0_2c = expectation_value(mps2c, H2c)
+    find_groundstate!(mps3c, H3c, DMRG2(; maxiter=100, trscheme=trunc_st))
+    E0_3c = expectation_value(mps3c, H3c)
+
+    # Assert that the ground state energies are equal up to a precision, display the 
+    # compression
+    @assert E0_1 ≈ E0_1c
+    @assert E0_2 ≈ E0_2c
+    @assert E0_3 ≈ E0_3c
+
+    println("Max dim uncompressed H1: $(max_dim_ham(H1)) || Max dim compressed H1c: 
+            $(max_dim_ham(H1c))" )
+    println("Max dim uncompressed H2: $(max_dim_ham(H2)) || Max dim compressed H2c: 
+            $(max_dim_ham(H2c))" )
+    println("Max dim uncompressed H3: $(max_dim_ham(H3)) || Max dim compressed H3c: 
+            $(max_dim_ham(H3c))" )
+
+end
+
+@testset "finite_mpo_compression_bilinear_biquadratic_model_energy" begin
+    M = 10
+    D = 6 # Max bond dimension
+    trunc_st = trunctol()
+
+    mps1 = FiniteMPS(M, ℂ^3, ℂ^D)
+    mps2 = FiniteMPS(M, U1Space(0=>1, 1=>1, -1=>1), U1Space(0=>D, 1=>D, -1=>D))
+    mps3 = FiniteMPS(M, SU2Space(1 => 1), SU2Space(1=>D))
+
+    H1 = bilinear_biquadratic_model_trivial_finite(M)
+    H2 = bilinear_biquadratic_model_U1_finite(M)
+    H3 = bilinear_biquadratic_model_SU2_finite(M)
+
+    # Find groundstate of non-compressed Hamiltonian
+    find_groundstate!(mps1, H1, DMRG2(; maxiter=100, trscheme=trunc_st))
+    E0_1 = expectation_value(mps1, H1)
+    find_groundstate!(mps2, H2, DMRG2(; maxiter=100, trscheme=trunc_st))
+    E0_2 = expectation_value(mps2, H2)
+    find_groundstate!(mps3, H3, DMRG2(; maxiter=100, trscheme=trunc_st))
+    E0_3 = expectation_value(mps3, H3)
+
+    # Compress MPO
+    mps1c = FiniteMPS(M, ℂ^3, ℂ^D)
+    mps2c = FiniteMPS(M, U1Space(0=>1, 1=>1, -1=>1), U1Space(0=>D, 1=>D, -1=>D))
+    mps3c = FiniteMPS(M, SU2Space(1 => 1), SU2Space(1=>D))
+
+    H1c, _ = mpo_compression(H1, 10^-10)
+    H2c, _ = mpo_compression(H2, 10^-10)
+    H3c, _ = mpo_compression(H3, 10^-10)
+
+    # Find groundstate of compressed Hamiltonian
+    find_groundstate!(mps1c, H1c, DMRG2(; maxiter=100, trscheme=trunc_st))
+    E0_1c = expectation_value(mps1c, H1c)
+    find_groundstate!(mps2c, H2c, DMRG2(; maxiter=100, trscheme=trunc_st))
+    E0_2c = expectation_value(mps2c, H2c)
+    find_groundstate!(mps3c, H3c, DMRG2(; maxiter=100, trscheme=trunc_st))
+    E0_3c = expectation_value(mps3c, H3c)
+
+    # Assert that the ground state energies are equal up to a precision, display the 
+    # compression
+    @assert E0_1 ≈ E0_1c
+    @assert E0_2 ≈ E0_2c
+    @assert E0_3 ≈ E0_3c
+
+    println("Max dim uncompressed H1: $(max_dim_ham(H1)) || Max dim compressed H1c: 
+            $(max_dim_ham(H1c))" )
+    println("Max dim uncompressed H2: $(max_dim_ham(H2)) || Max dim compressed H2c: 
+            $(max_dim_ham(H2c))" )
+    println("Max dim uncompressed H3: $(max_dim_ham(H3)) || Max dim compressed H3c: 
+            $(max_dim_ham(H3c))" )
+
+end
+
+@testset "finite_mpo_compression_quantum_potts_model_energy" begin
+    M = 10
+    D = 6 # Max bond dimension
+    trunc_st = trunctol()
+
+    mps1 = FiniteMPS(M, ℂ^3, ℂ^D)
+    mps2 = FiniteMPS(M, Z3Space(0=>1, 1=>1, 2=>1), Z3Space(0=>1, 1=>1, 2=>1))
+
+    H1 = quantum_potts_trivial_finite(M)
+    H2 = quantum_potts_Z3_finite(M)
+
+    # Find groundstate of non-compressed Hamiltonian
+    find_groundstate!(mps1, H1, DMRG2(; maxiter=100, trscheme=trunc_st))
+    E0_1 = expectation_value(mps1, H1)
+    find_groundstate!(mps2, H2, DMRG2(; maxiter=100, trscheme=trunc_st))
+    E0_2 = expectation_value(mps2, H2)
+
+    # Compress MPO
+    mps1c = FiniteMPS(M, ℂ^3, ℂ^D)
+    mps2c = FiniteMPS(M, Z3Space(0=>1, 1=>1, 2=>1), Z3Space(0=>1, 1=>1, 2=>1))
+
+    H1c, _ = mpo_compression(H1, 10^-10)
+    H2c, _ = mpo_compression(H2, 10^-10)
+
+    # Find groundstate of compressed Hamiltonian
+    find_groundstate!(mps1c, H1c, DMRG2(; maxiter=100, trscheme=trunc_st))
+    E0_1c = expectation_value(mps1c, H1c)
+    find_groundstate!(mps2c, H2c, DMRG2(; maxiter=100, trscheme=trunc_st))
+    E0_2c = expectation_value(mps2c, H2c)
+
+    # Assert that the ground state energies are equal up to a precision, display the 
+    # compression
+    @assert E0_1 ≈ E0_1c
+    @assert E0_2 ≈ E0_2c
+
+    println("Max dim uncompressed H1: $(max_dim_ham(H1)) || Max dim compressed H1c: 
+            $(max_dim_ham(H1c))" )
+    println("Max dim uncompressed H2: $(max_dim_ham(H2)) || Max dim compressed H2c: 
+            $(max_dim_ham(H2c))" )
+
+end
+
+# Test random finite Hamiltonian
+@testset "finite_random_long_range_ham" begin
+    # test the finite compression for a long range random Ising model
+    # check if the energy after the compression is the same as before
+    Random.seed!()
+    D = 20 # max bond dimensions
+    L = 10 # number of sites
+    range = 4 # cutoff for the maximum range between interactions
+    trunc_st = notrunc()
+
+    mps = FiniteMPS(L, Z2Space(0=>1, 1=>1), Z2Space(0=>D, 1=>D))
+    H = create_long_range_ising_symmetries_random_finite(L, range, 0.3, 0.4)
+    find_groundstate!(mps, H, DMRG2(; maxiter=100, trscheme=trunc_st))
+    E0 = expectation_value(mps, H)
+    #println("<mps|H|mps> = $real(E0)")
+
+    mps2 = FiniteMPS(L, Z2Space(0=>1, 1=>1), Z2Space(0=>D, 1=>D))
+    H2, Rs = mpo_compression(H)
+    find_groundstate!(mps2, H2, DMRG2(; maxiter=100, trscheme=trunc_st))
+    E1 = expectation_value(mps2, H2)
+end
\ No newline at end of file
diff --git a/test/operators/infinite_mpo_compression.jl b/test/operators/infinite_mpo_compression.jl
new file mode 100644
index 000000000..caebc96c6
--- /dev/null
+++ b/test/operators/infinite_mpo_compression.jl
@@ -0,0 +1,388 @@
+using Test, Random
+using TensorKit
+using MPSKit
+using MPSKitModels
+
+function max_dim_ham(H)
+    N = length(H)
+    max_dim = 1
+    for i in 1:N
+        if max_dim < dims(H[i])[4]
+            max_dim = dims(H[i])[4]
+        end
+    end
+    return max_dim
+end
+
+# Utility functions to create Hamiltonians for testing
+
+function create_long_range_ising_symmetries_infinite(k)
+    infinite_chain = PeriodicVector([Z2Space(0=>1, 1=>1)])
+    ZZ = S_zz(Z2Irrep)
+    X = S_x(Z2Irrep)
+    two_site_operators = [(1,j) => -ZZ for j in 1:k]
+    H = InfiniteMPOHamiltonian(infinite_chain, two_site_operators...)
+    return H
+end
+
+function create_long_range_ising_symmetries_infinite_msite_random(M, k, σⱼ, σₕ)
+    V = Z2Space(0=>1, 1=>1)
+    infinite_chain = PeriodicVector(fill(V, M))
+    ZZ = S_zz(Z2Irrep)
+    X = S_x(Z2Irrep)
+    single_site_operators = [i => (σₕ * randn()) * X for i in 1:M]
+    two_site_operators = [(i,i+j) => (1 + σⱼ * randn()) * ZZ for i in 1:M for j in 1:k]
+    H = InfiniteMPOHamiltonian(infinite_chain, single_site_operators...,
+                               two_site_operators...)
+    return H
+end
+
+function create_long_range_ising_symmetries_infinite_twosite_diff_int_strength(k, J, h)
+    infinite_chain = PeriodicVector([Z2Space(0=>1, 1=>1), Z2Space(0=>1, 1=>1)])
+    ZZ = S_zz(Z2Irrep)
+    X = S_x(Z2Irrep)
+    two_site_operators_J12 = [(a,a+j) => J[1] * ZZ for a in [1 2] for j in 1:2:2k]
+    two_site_operators_J11 = [(1,j) => J[2] * ZZ for j in 3:2:2k]
+    two_site_operators_J22 = [(2,j) => J[3] * ZZ for j in 4:2:2k]
+    single_site_operators = [1 => h[1] * X, 2 => h[2] * X]
+    H = InfiniteMPOHamiltonian(infinite_chain, single_site_operators...,
+                                two_site_operators_J11..., two_site_operators_J12...,
+                                two_site_operators_J22...)
+    return H
+end
+
+function transverse_field_ising_fermion_parity_infinite(M)
+    chain = InfiniteChain(M)
+    return transverse_field_ising(ComplexF64, FermionParity, chain)
+end
+
+function transverse_field_ising_z2_infinite(M)
+    chain = InfiniteChain(M)
+    return transverse_field_ising(ComplexF64, Z2Irrep, chain)
+end
+
+function transverse_field_ising_trivial_infinite(M)
+    chain = InfiniteChain(M)
+    return transverse_field_ising(ComplexF64, Trivial, chain)
+end
+
+function kitaev_model_infinite(M)
+    chain = InfiniteChain(M)
+    return kitaev_model(ComplexF64, chain)
+end
+
+function heisenberg_XXX_trivial_infinite(M)
+    chain = InfiniteChain(M)
+    return heisenberg_XXX(Trivial, chain)
+end
+
+function heisenberg_XXX_U1_infinite(M)
+    chain = InfiniteChain(M)
+    return heisenberg_XXX(U1Irrep, chain)
+end
+
+function heisenberg_XXX_SU2_infinite(M)
+    chain = InfiniteChain(M)
+    return heisenberg_XXX(SU2Irrep, chain)
+end
+
+function bilinear_biquadratic_model_trivial_infinite(M)
+    chain = InfiniteChain(M)
+    return bilinear_biquadratic_model(ComplexF64, Trivial, chain)
+end
+
+function quantum_potts_trivial_infinite(M)
+    chain = InfiniteChain(M)
+    return quantum_potts(ComplexF64, Trivial, chain)
+end
+
+function quantum_potts_Z3_infinite(M)
+    chain = InfiniteChain(M)
+    return quantum_potts(ComplexF64, Z3Irrep, chain)
+end
+
+@testset "infinite_mpo_compression_energy" begin
+    # test the infinite compression for a long range Ising model
+    # check if the energy after the compression is the same as before
+    D = 4 # max bond dimensions
+    range = 2 # cutoff for the maximum range between interactions
+    trunc_st = notrunc()
+
+    mps = InfiniteMPS(Z2Space(0=>1, 1=>1), Z2Space(0=>D, 1=>D))
+    H = create_long_range_ising_symmetries_infinite(range)
+    mps, = find_groundstate(mps, H, VUMPS(;maxiter=10))
+    E0 = expectation_value(mps, H)
+    #println("<mps|H|mps> = $real(E0)")
+
+    # mps2 = InfiniteMPS(Z2Space(0=>1, 1=>1), Z2Space(0=>D, 1=>D))
+    # H2, Rs = mpo_finite_compression(H)
+    # find_groundstate!(mps2, H2, DMRG2(; maxiter=100, trscheme=trunc_st))
+    # E1 = expectation_value(mps2, H2)
+    mps2 = InfiniteMPS(Z2Space(0=>1, 1=>1), Z2Space(0=>D, 1=>D))
+    Q, P = mpo_compression(H)
+    H2 = Q
+    mps2, = find_groundstate(mps2, H2, VUMPS(;maxiter=100))
+    E0 = expectation_value(mps2, H2)
+end
+
+@testset "infinite_mpo_compression_energy_two_site" begin
+    # test the infinite compression for a long range Ising model
+    # check if the energy after the compression is the same as before
+    D = 10 # max bond dimensions
+    range = 4 # cutoff for the maximum range between interactions
+    trunc_st = notrunc()
+    J = [1, 1, 1]
+    h = [0.001, 0.001]
+
+    spacephys = fill(Z2Space(0=>1, 1=>1), 2)
+    space_virt = fill(Z2Space(0=>D, 1=>D), 2)
+
+    mps = InfiniteMPS(spacephys, space_virt)
+    H = create_long_range_ising_symmetries_infinite_twosite_diff_int_strength(range, J, h)
+    mps, = find_groundstate(mps, H, VUMPS(;maxiter=100))
+    E0 = expectation_value(mps, H)
+    #println("<mps|H|mps> = $real(E0)")
+
+    # mps2 = InfiniteMPS(Z2Space(0=>1, 1=>1), Z2Space(0=>D, 1=>D))
+    # H2, Rs = mpo_finite_compression(H)
+    # find_groundstate!(mps2, H2, DMRG2(; maxiter=100, trscheme=trunc_st))
+    # E1 = expectation_value(mps2, H2)
+    mps2 = InfiniteMPS(spacephys, space_virt)
+    Qs, Ps = mpo_compression(H, 0)
+    H2 = Qs
+    mps2, = find_groundstate(mps2, H2, VUMPS(;maxiter=100))
+    E0 = expectation_value(mps2, H2)
+end
+
+@testset "infinite_mpo_compression_transverse_ising_model_energy" begin
+    M = 4
+    D = 10 # Max bond dimension
+
+    mps1 = InfiniteMPS(fill(ℂ^2, M), fill(ℂ^D, M))
+    mps2 = InfiniteMPS(fill(Z2Space(0=>1, 1=>1), M), fill(Z2Space(0=>D, 1=>D), M))
+    mps3 = InfiniteMPS(fill(Vect[FermionParity](0=>1, 1=>1), M), fill(Vect[FermionParity](
+        0=>D, 1=>D), M))
+
+    H1 = transverse_field_ising_trivial_infinite(M)
+    H2 = transverse_field_ising_z2_infinite(M)
+    H3 = transverse_field_ising_fermion_parity_infinite(M)
+
+    # Find groundstate of non-compressed Hamiltonian
+    mps1, = find_groundstate(mps1, H1, VUMPS(; maxiter=100))
+    E0_1 = expectation_value(mps1, H1)
+    mps2, = find_groundstate(mps2, H2, VUMPS(; maxiter=100))
+    E0_2 = expectation_value(mps2, H2)
+    mps3, = find_groundstate(mps3, H3, VUMPS(; maxiter=100))
+    E0_3 = expectation_value(mps3, H3)
+
+    # Compress MPO
+    mps1c = InfiniteMPS(fill(ℂ^2, M), fill(ℂ^D, M))
+    mps2c = InfiniteMPS(fill(Z2Space(0=>1, 1=>1), M), fill(Z2Space(0=>D, 1=>D), M))
+    mps3c = InfiniteMPS(fill(Vect[FermionParity](0=>1, 1=>1), M), fill(Vect[FermionParity](
+        0=>D, 1=>D), M))
+
+    H1c, _ = mpo_compression(H1, 10^-10)
+    H2c, _ = mpo_compression(H2, 10^-10)
+    H3c, _ = mpo_compression(H3, 10^-10)
+
+    # Find groundstate of compressed Hamiltonian
+    mps1c, = find_groundstate(mps1c, H1c, VUMPS(; maxiter=100))
+    E0_1c = expectation_value(mps1c, H1c)
+    mps2c, = find_groundstate(mps2c, H2c, VUMPS(; maxiter=100))
+    E0_2c = expectation_value(mps2c, H2c)
+    mps3c, = find_groundstate(mps3c, H3c, VUMPS(; maxiter=100))
+    E0_3c = expectation_value(mps3c, H3c)
+
+    # Assert that the ground state energies are equal up to a precision, display the 
+    # compression
+    @assert E0_1 ≈ E0_1c
+    @assert E0_2 ≈ E0_2c
+    @assert E0_3 ≈ E0_3c
+
+    println("Max dim uncompressed H1: $(max_dim_ham(H1)) || Max dim compressed H1c: 
+            $(max_dim_ham(H1c))" )
+    println("Max dim uncompressed H2: $(max_dim_ham(H2)) || Max dim compressed H2c: 
+            $(max_dim_ham(H2c))" )
+    println("Max dim uncompressed H3: $(max_dim_ham(H3)) || Max dim compressed H3c: 
+            $(max_dim_ham(H3c))" )
+
+end
+
+@testset "infinite_mpo_compression_kitaev_model_energy" begin
+    M = 10
+    D = 10 # Max bond dimension
+
+    mps1 = InfiniteMPS(fill(Vect[FermionParity](0=>1, 1=>1), M), 
+                       fill(Vect[FermionParity](0=>D, 1=>D), M))
+
+    H1 = kitaev_model_infinite(M)
+
+    # Find groundstate of non-compressed Hamiltonian
+    mps1, = find_groundstate(mps1, H1, VUMPS(; maxiter=100))
+    E0_1 = expectation_value(mps1, H1)
+
+    # Compress MPO
+    mps1c = InfiniteMPS(fill(Vect[FermionParity](0=>1, 1=>1), M), 
+                       fill(Vect[FermionParity](0=>D, 1=>D), M))
+
+    H1c, _ = mpo_compression(H1, 10^-10)
+
+    # Find groundstate of compressed Hamiltonian
+    mps1c, = find_groundstate(mps1c, H1c, VUMPS(; maxiter=100))
+    E0_1c = expectation_value(mps1c, H1c)
+
+    # Assert that the ground state energies are equal up to a precision, display the 
+    # compression
+    @assert E0_1 ≈ E0_1c
+
+    println("Max dim uncompressed H1: $(max_dim_ham(H1)) || Max dim compressed H1c: 
+            $(max_dim_ham(H1c))" )
+
+end
+
+@testset "infinite_mpo_compression_heisenberg_XXX_energy" begin
+    M = 10
+    D = 6 # Max bond dimension
+
+    mps1 = InfiniteMPS(fill(ℂ^3, M),fill(ℂ^D, M))
+    mps2 = InfiniteMPS(fill(U1Space(0=>1, 1=>1, -1=>1), M), fill(U1Space(0=>D, 1=>D, -1=>D), M))
+    mps3 = InfiniteMPS(fill(SU2Space(1 => 1), M), fill(SU2Space(1=>D), M))
+
+    H1 = heisenberg_XXX_trivial_infinite(M)
+    H2 = heisenberg_XXX_U1_infinite(M)
+    H3 = heisenberg_XXX_SU2_infinite(M)
+
+    # Find groundstate of non-compressed Hamiltonian
+    mps1, = find_groundstate(mps1, H1, VUMPS(; maxiter=100))
+    E0_1 = expectation_value(mps1, H1)
+    mps2, = find_groundstate(mps2, H2, VUMPS(; maxiter=100))
+    E0_2 = expectation_value(mps2, H2)
+    mps3, = find_groundstate(mps3, H3, VUMPS(; maxiter=100))
+    E0_3 = expectation_value(mps3, H3)
+
+    # Compress MPO
+    mps1c = InfiniteMPS(fill(ℂ^3, M),fill(ℂ^D, M))
+    mps2c = InfiniteMPS(fill(U1Space(0=>1, 1=>1, -1=>1), M), fill(U1Space(0=>D, 1=>D, -1=>D), M))
+    mps3c = InfiniteMPS(fill(SU2Space(1 => 1), M), fill(SU2Space(1=>D), M))
+
+    H1c, _ = mpo_compression(H1, 10^-10)
+    H2c, _ = mpo_compression(H2, 10^-10)
+    H3c, _ = mpo_compression(H3, 10^-10)
+
+    # Find groundstate of compressed Hamiltonian
+    mps1c, = find_groundstate(mps1c, H1c, VUMPS(; maxiter=100))
+    E0_1c = expectation_value(mps1c, H1c)
+    mps2c, = find_groundstate(mps2c, H2c, VUMPS(; maxiter=100))
+    E0_2c = expectation_value(mps2c, H2c)
+    mps3c, = find_groundstate(mps3c, H3c, VUMPS(; maxiter=100))
+    E0_3c = expectation_value(mps3c, H3c)
+
+    # Assert that the ground state energies are equal up to a precision, display the 
+    # compression
+    @assert E0_1 ≈ E0_1c
+    @assert E0_2 ≈ E0_2c
+    @assert E0_3 ≈ E0_3c
+
+    println("Max dim uncompressed H1: $(max_dim_ham(H1)) || Max dim compressed H1c: 
+            $(max_dim_ham(H1c))" )
+    println("Max dim uncompressed H2: $(max_dim_ham(H2)) || Max dim compressed H2c: 
+            $(max_dim_ham(H2c))" )
+    println("Max dim uncompressed H3: $(max_dim_ham(H3)) || Max dim compressed H3c: 
+            $(max_dim_ham(H3c))" )
+
+end
+
+@testset "infinite_mpo_compression_bilinear_biquadratic_model_energy" begin
+    M = 10
+    D = 6 # Max bond dimension
+
+    mps1 = InfiniteMPS(fill(ℂ^3, M),fill(ℂ^D, M))
+
+
+    H1 = bilinear_biquadratic_model_trivial_infinite(M)
+
+    # Find groundstate of non-compressed Hamiltonian
+    mps1, = find_groundstate(mps1, H1, VUMPS(; maxiter=100))
+    E0_1 = expectation_value(mps1, H1)
+
+    # Compress MPO
+    mps1c = InfiniteMPS(fill(ℂ^3, M),fill(ℂ^D, M))
+
+    H1c, _ = mpo_compression(H1, 10^-10)
+
+    # Find groundstate of compressed Hamiltonian
+    mps1c, = find_groundstate(mps1c, H1c, VUMPS(; maxiter=100))
+    E0_1c = expectation_value(mps1c, H1c)
+
+    # Assert that the ground state energies are equal up to a precision, display the 
+    # compression
+    @assert E0_1 ≈ E0_1c
+
+    println("Max dim uncompressed H1: $(max_dim_ham(H1)) || Max dim compressed H1c: 
+            $(max_dim_ham(H1c))" )
+
+end
+
+@testset "infinite_mpo_compression_quantum_potts_model_energy" begin
+    M = 10
+    D = 6 # Max bond dimension
+
+    mps1 = InfiniteMPS(fill(ℂ^3, M), fill(ℂ^D, M))
+    mps2 = InfiniteMPS(fill(Z3Space(0=>1, 1=>1, 2=>1), M), fill(Z3Space(0=>1, 1=>1, 2=>1), M))
+
+    H1 = quantum_potts_trivial_infinite(M)
+    H2 = quantum_potts_Z3_infinite(M)
+
+    # Find groundstate of non-compressed Hamiltonian
+    mps1, = find_groundstate(mps1, H1, VUMPS(; maxiter=100))
+    E0_1 = expectation_value(mps1, H1)
+    mps2, = find_groundstate(mps2, H2, VUMPS(; maxiter=100))
+    E0_2 = expectation_value(mps2, H2)
+
+    # Compress MPO
+    mps1c = InfiniteMPS(fill(ℂ^3, M), fill(ℂ^D, M))
+    mps2c = InfiniteMPS(fill(Z3Space(0=>1, 1=>1, 2=>1), M), fill(Z3Space(0=>1, 1=>1, 2=>1), M))
+
+    H1c, _ = mpo_compression(H1, 10^-10)
+    H2c, _ = mpo_compression(H2, 10^-10)
+
+    # Find groundstate of compressed Hamiltonian
+    mps1c, = find_groundstate(mps1c, H1c, VUMPS(; maxiter=100))
+    E0_1c = expectation_value(mps1c, H1c)
+    mps2c, = find_groundstate(mps2c, H2c, VUMPS(; maxiter=100))
+    E0_2c = expectation_value(mps2c, H2c)
+
+    # Assert that the ground state energies are equal up to a precision, display the 
+    # compression
+    @assert E0_1 ≈ E0_1c
+    @assert E0_2 ≈ E0_2c
+
+    println("Max dim uncompressed H1: $(max_dim_ham(H1)) || Max dim compressed H1c: 
+            $(max_dim_ham(H1c))" )
+    println("Max dim uncompressed H2: $(max_dim_ham(H2)) || Max dim compressed H2c: 
+            $(max_dim_ham(H2c))" )
+
+end
+
+# Test random infinite Hamiltonian
+@testset "infinite_random_long_range_ham" begin
+    # test the finite compression for a long range random Ising model
+    # check if the energy after the compression is the same as before
+    Random.seed!()
+    D = 20 # max bond dimensions
+    L = 10 # number of sites
+    range = 4 # cutoff for the maximum range between interactions
+
+    mps = InfiniteMPS(fill(Z2Space(0=>1, 1=>1),L), fill(Z2Space(0=>D, 1=>D),L))
+    H = create_long_range_ising_symmetries_infinite_msite_random(L, range, 0.1, 0.001)
+    mps, = find_groundstate(mps, H, VUMPS(; maxiter=100))
+    E0 = expectation_value(mps, H)
+    #println("<mps|H|mps> = $real(E0)")
+
+    mps2 = InfiniteMPS(fill(Z2Space(0=>1, 1=>1),L), fill(Z2Space(0=>D, 1=>D),L))
+    H2, Rs = mpo_compression(H, 0)
+    mps2, = find_groundstate(mps2, H2, VUMPS(; maxiter=100))
+    E1 = expectation_value(mps2, H2)
+
+    @assert E0 ≈ E1
+end
\ No newline at end of file