diff --git a/Project.toml b/Project.toml
index 23fc5e15f..d2125adcc 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "MPSKit"
 uuid = "bb1c41ca-d63c-52ed-829e-0820dda26502"
-authors = "Lukas Devos, Maarten Van Damme and contributors"
 version = "0.13.8"
+authors = "Lukas Devos, Maarten Van Damme and contributors"
 
 [deps]
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
@@ -26,7 +26,7 @@ VectorInterface = "409d34a3-91d5-4945-b6ec-7529ddf182d8"
 [compat]
 Accessors = "0.1"
 Aqua = "0.8.9"
-BlockTensorKit = "0.3"
+BlockTensorKit = "0.3.4"
 Combinatorics = "1"
 Compat = "3.47, 4.10"
 DocStringExtensions = "0.9.3"
@@ -34,7 +34,7 @@ HalfIntegers = "1.6.0"
 KrylovKit = "0.8.3, 0.9.2, 0.10"
 LinearAlgebra = "1.6"
 LoggingExtras = "~1.0"
-MatrixAlgebraKit = "0.5.0"
+MatrixAlgebraKit = "0.6"
 OhMyThreads = "0.7, 0.8"
 OptimKit = "0.3.1, 0.4"
 Pkg = "1"
@@ -42,7 +42,7 @@ Plots = "1.40"
 Printf = "1"
 Random = "1"
 RecipesBase = "1.1"
-TensorKit = "0.15.1"
+TensorKit = "0.16"
 TensorKitManifolds = "0.7"
 TensorOperations = "5"
 Test = "1"
diff --git a/src/algorithms/changebonds/randexpand.jl b/src/algorithms/changebonds/randexpand.jl
index 09a6c12da..80003ead2 100644
--- a/src/algorithms/changebonds/randexpand.jl
+++ b/src/algorithms/changebonds/randexpand.jl
@@ -82,7 +82,8 @@ Sample basis states within a given `V::VectorSpace` by creating weights for each
 are distributed uniformly, and then truncating according to the given `strategy`.
 """
 function sample_space(V, strategy)
-    S = TensorKit.SectorDict(c => Random.rand(dim(V, c)) for c in sectors(V))
+    S = TensorKit.SectorVector{Float64}(undef, V)
+    Random.rand!(parent(S))
     ind = MatrixAlgebraKit.findtruncated(S, strategy)
     return TensorKit.Factorizations.truncate_space(V, ind)
 end
diff --git a/src/algorithms/changebonds/svdcut.jl b/src/algorithms/changebonds/svdcut.jl
index 1324cecb0..d5d4ea48a 100644
--- a/src/algorithms/changebonds/svdcut.jl
+++ b/src/algorithms/changebonds/svdcut.jl
@@ -120,12 +120,13 @@ function changebonds!(H::FiniteMPOHamiltonian, alg::SvdCut)
     end
 
     # swipe right
+    alg_trunc = MatrixAlgebraKit.TruncatedAlgorithm(alg.alg_svd, alg.trscheme)
     for i in 1:(length(H) - 1)
-        H = left_canonicalize!(H, i; alg = alg.alg_svd, alg.trscheme)
+        H = left_canonicalize!(H, i; alg = alg_trunc)
     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)
+        H = right_canonicalize!(H, i; alg = alg_trunc)
     end
 
     return H
diff --git a/src/algorithms/groundstate/idmrg.jl b/src/algorithms/groundstate/idmrg.jl
index 6cb71afd3..fab2c96da 100644
--- a/src/algorithms/groundstate/idmrg.jl
+++ b/src/algorithms/groundstate/idmrg.jl
@@ -159,9 +159,9 @@ function _localupdate_sweep_idmrg!(ψ, H, envs, alg_eigsolve)
         _, ψ.AC[pos] = fixedpoint(h, ψ.AC[pos], :SR, alg_eigsolve)
         if pos == length(ψ)
             # AC needed in next sweep
-            ψ.AL[pos], ψ.C[pos] = left_orth(ψ.AC[pos])
+            ψ.AL[pos], ψ.C[pos] = left_orth(ψ.AC[pos]; positive = true)
         else
-            ψ.AL[pos], ψ.C[pos] = left_orth!(ψ.AC[pos])
+            ψ.AL[pos], ψ.C[pos] = left_orth!(ψ.AC[pos]; positive = true)
         end
         transfer_leftenv!(envs, ψ, H, ψ, pos + 1)
     end
@@ -171,7 +171,7 @@ function _localupdate_sweep_idmrg!(ψ, H, envs, alg_eigsolve)
         h = AC_hamiltonian(pos, ψ, H, ψ, envs)
         E, ψ.AC[pos] = fixedpoint(h, ψ.AC[pos], :SR, alg_eigsolve)
 
-        ψ.C[pos - 1], temp = right_orth!(_transpose_tail(ψ.AC[pos]; copy = (pos == 1)))
+        ψ.C[pos - 1], temp = right_orth!(_transpose_tail(ψ.AC[pos]; copy = (pos == 1)); positive = true)
         ψ.AR[pos] = _transpose_front(temp)
 
         transfer_rightenv!(envs, ψ, H, ψ, pos - 1)
diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl
index bf10f06da..0f9f7bacc 100644
--- a/src/algorithms/toolbox.jl
+++ b/src/algorithms/toolbox.jl
@@ -9,7 +9,7 @@ entropy(state::InfiniteMPS) = map(Base.Fix1(entropy, state), 1:length(state))
 function entropy(state::Union{FiniteMPS, WindowMPS, InfiniteMPS}, loc::Int)
     S = zero(real(scalartype(state)))
     tol = eps(typeof(S))
-    for (c, b) in entanglement_spectrum(state, loc)
+    for (c, b) in pairs(entanglement_spectrum(state, loc))
         s = zero(S)
         for x in b
             x < tol && break
diff --git a/src/operators/ortho.jl b/src/operators/ortho.jl
index d4f0f607f..34ce1b159 100644
--- a/src/operators/ortho.jl
+++ b/src/operators/ortho.jl
@@ -1,6 +1,6 @@
 function left_canonicalize!(
         H::FiniteMPOHamiltonian, i::Int;
-        alg = Defaults.alg_qr(), trscheme::TruncationStrategy = notrunc()
+        alg::MatrixAlgebraKit.AbstractAlgorithm = Defaults.alg_qr()
     )
     1 ≤ i < length(H) || throw(ArgumentError("Bounds error in canonicalize"))
 
@@ -19,7 +19,7 @@ function left_canonicalize!(
     # QR of second column
     if size(W, 1) == 1
         tmp = transpose(C′, ((2, 1), (3,)))
-        Q, R = _left_orth!(tmp; alg, trunc = trscheme)
+        Q, R = left_orth!(tmp; alg)
 
         if dim(R) == 0 # fully truncated
             V = oneunit(S) ⊞ oneunit(S)
@@ -35,7 +35,7 @@ function left_canonicalize!(
         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,)))
-        Q, R = _left_orth!(tmp; alg, trunc = trscheme)
+        Q, R = left_orth!(tmp; alg)
 
         if dim(R) == 0 # fully truncated
             V = oneunit(S) ⊞ oneunit(S)
@@ -86,7 +86,7 @@ end
 
 function right_canonicalize!(
         H::FiniteMPOHamiltonian, i::Int;
-        alg = Defaults.alg_lq(), trscheme::TruncationStrategy = notrunc()
+        alg::MatrixAlgebraKit.AbstractAlgorithm = Defaults.alg_lq()
     )
     1 < i ≤ length(H) || throw(ArgumentError("Bounds error in canonicalize"))
 
@@ -105,7 +105,7 @@ function right_canonicalize!(
     # LQ of second row
     if size(W, 4) == 1
         tmp = transpose(B′, ((1,), (3, 2)))
-        R, Q = _right_orth!(tmp; alg, trunc = trscheme)
+        R, Q = right_orth!(tmp; alg)
 
         if dim(R) == 0
             V = oneunit(S) ⊞ oneunit(S)
@@ -121,7 +121,7 @@ function right_canonicalize!(
         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)))
-        R, Q = _right_orth!(tmp; alg, trunc = trscheme)
+        R, Q = right_orth!(tmp; alg)
         if dim(R) == 0
             V = oneunit(S) ⊞ oneunit(S)
             Q1 = typeof(W.A)(undef, SumSpace{S}() ⊗ physicalspace(W) ← physicalspace(W) ⊗ SumSpace{S}())
diff --git a/src/states/abstractmps.jl b/src/states/abstractmps.jl
index 5c1e82c1d..3fa45832d 100644
--- a/src/states/abstractmps.jl
+++ b/src/states/abstractmps.jl
@@ -119,11 +119,11 @@ function makefullrank!(A::PeriodicVector{<:GenericMPSTensor}; alg = Defaults.alg
         i = findfirst(!isfullrank, A)
         isnothing(i) && break
         if !isfullrank(A[i]; side = :left)
-            L, Q = _right_orth!(_transpose_tail(A[i]); alg)
+            L, Q = right_orth!(_transpose_tail(A[i]); alg)
             A[i] = _transpose_front(Q)
             A[i - 1] = A[i - 1] * L
         else
-            A[i], R = _left_orth!(A[i]; alg)
+            A[i], R = left_orth!(A[i]; alg)
             A[i + 1] = _transpose_front(R * _transpose_tail(A[i + 1]))
         end
     end
diff --git a/src/states/ortho.jl b/src/states/ortho.jl
index a29942339..320f90209 100644
--- a/src/states/ortho.jl
+++ b/src/states/ortho.jl
@@ -164,8 +164,8 @@ regauge!
 function regauge!(
         AC::GenericMPSTensor, C::MPSBondTensor; alg = Defaults.alg_qr()
     )
-    Q_AC, _ = _left_orth!(AC; alg)
-    Q_C, _ = _left_orth!(C; alg)
+    Q_AC, _ = left_orth!(AC; alg)
+    Q_C, _ = left_orth!(C; alg)
     return mul!(AC, Q_AC, Q_C')
 end
 function regauge!(AC::Vector{<:GenericMPSTensor}, C::Vector{<:MPSBondTensor}; kwargs...)
@@ -178,8 +178,8 @@ function regauge!(
         CL::MPSBondTensor, AC::GenericMPSTensor; alg = Defaults.alg_lq()
     )
     AC_tail = _transpose_tail(AC)
-    _, Q_AC = _right_orth!(AC_tail; alg)
-    _, Q_C = _right_orth!(CL; alg)
+    _, Q_AC = right_orth!(AC_tail; alg)
+    _, Q_C = right_orth!(CL; alg)
     AR_tail = mul!(AC_tail, Q_C', Q_AC)
     return repartition!(AC, AR_tail)
 end
@@ -240,7 +240,7 @@ function gauge_eigsolve_step!(it::IterativeSolver{LeftCanonical}, state)
     if iter ≥ it.eig_miniter
         alg_eigsolve = updatetol(it.alg_eigsolve, 1, ϵ^2)
         _, vec = fixedpoint(flip(TransferMatrix(A, AL)), C[end], :LM, alg_eigsolve)
-        _, C[end] = _left_orth!(vec; alg = it.alg_orth)
+        _, C[end] = left_orth!(vec; alg = it.alg_orth)
     end
     return C[end]
 end
@@ -251,7 +251,7 @@ function gauge_orth_step!(it::IterativeSolver{LeftCanonical}, state)
         # repartition!(A_tail[i], AL[i])
         mul!(CA_tail[i], C[i - 1], A_tail[i])
         repartition!(AL[i], CA_tail[i])
-        AL[i], C[i] = _left_orth!(AL[i]; alg = it.alg_orth)
+        AL[i], C[i] = left_orth!(AL[i]; alg = it.alg_orth)
     end
     normalize!(C[end])
     return C[end]
@@ -298,7 +298,7 @@ function gauge_eigsolve_step!(it::IterativeSolver{RightCanonical}, state)
     if iter ≥ it.eig_miniter
         alg_eigsolve = updatetol(it.alg_eigsolve, 1, ϵ^2)
         _, vec = fixedpoint(TransferMatrix(A, AR), C[end], :LM, alg_eigsolve)
-        C[end], _ = _right_orth!(vec; alg = it.alg_orth)
+        C[end], _ = right_orth!(vec; alg = it.alg_orth)
     end
     return C[end]
 end
@@ -308,7 +308,7 @@ function gauge_orth_step!(it::IterativeSolver{RightCanonical}, state)
     for i in length(AR):-1:1
         AC = mul!(AR[i], A[i], C[i])   # use AR as temporary storage for A * C
         tmp = repartition!(AC_tail[i], AC)
-        C[i - 1], tmp = _right_orth!(tmp; alg = it.alg_orth)
+        C[i - 1], tmp = right_orth!(tmp; alg = it.alg_orth)
         repartition!(AR[i], tmp)       # TODO: avoid doing this every iteration
     end
     normalize!(C[end])
diff --git a/src/utility/plotting.jl b/src/utility/plotting.jl
index 9265311fa..3e2ce109b 100644
--- a/src/utility/plotting.jl
+++ b/src/utility/plotting.jl
@@ -36,7 +36,7 @@ function entanglementplot end
     spectra = entanglement_spectrum(mps, site)
     sectors = []
     spectrum = []
-    for (c, b) in spectra
+    for (c, b) in pairs(spectra)
         if expand_symmetry # Duplicate entries according to the quantum dimension.
             b′ = repeat(b, dim(c))
             sort!(b′; rev = true)
diff --git a/src/utility/utility.jl b/src/utility/utility.jl
index 418e0cd68..49cb96648 100644
--- a/src/utility/utility.jl
+++ b/src/utility/utility.jl
@@ -149,31 +149,3 @@ function check_unambiguous_braiding(V::VectorSpace)
     return check_unambiguous_braiding(Bool, V) ||
         throw(ArgumentError("cannot unambiguously braid $V"))
 end
-
-# temporary workaround for the fact that left_orth and right_orth are poorly designed:
-function _left_orth!(t; alg::MatrixAlgebraKit.AbstractAlgorithm, trunc::MatrixAlgebraKit.TruncationStrategy = notrunc())
-    if alg isa LAPACK_HouseholderQR
-        return left_orth!(t; kind = :qr, alg_qr = alg, trunc)
-    elseif alg isa LAPACK_HouseholderLQ
-        return left_orth!(t; kind = :qr, alg_qr = LAPACK_HouseholderQR(; alg.kwargs...), trunc)
-    elseif alg isa PolarViaSVD
-        return left_orth!(t; kind = :polar, alg_polar = alg, trunc)
-    elseif alg isa LAPACK_SVDAlgorithm
-        return left_orth!(t; kind = :svd, alg_svd = alg, trunc)
-    else
-        error(lazy"unkown algorithm $alg")
-    end
-end
-function _right_orth!(t; alg::MatrixAlgebraKit.AbstractAlgorithm, trunc::TruncationStrategy = notrunc())
-    if alg isa LAPACK_HouseholderLQ
-        return right_orth!(t; kind = :lq, alg_lq = alg, trunc)
-    elseif alg isa LAPACK_HouseholderQR
-        return right_orth!(t; kind = :lq, alg_lq = LAPACK_HouseholderLQ(; alg.kwargs...), trunc)
-    elseif alg isa PolarViaSVD
-        return right_orth!(t; kind = :polar, alg_polar = alg, trunc)
-    elseif alg isa LAPACK_SVDAlgorithm
-        return right_orth!(t; kind = :svd, alg_svd = alg, trunc)
-    else
-        error(lazy"unkown algorithm $alg")
-    end
-end
diff --git a/test/algorithms.jl b/test/algorithms.jl
index c1d0ad9e9..207857951 100644
--- a/test/algorithms.jl
+++ b/test/algorithms.jl
@@ -1047,7 +1047,10 @@ module TestAlgorithms
         H = XY_model(U1Irrep; L)
 
         H_dense = convert(TensorMap, H)
-        vals_dense = TensorKit.SectorDict(c => sort(v; by = real) for (c, v) in eigvals(H_dense))
+        vals_dense = eigvals(H_dense)
+        for v in values(vals_dense)
+            sort!(v; by = real)
+        end
 
         tol = 1.0e-18 # tolerance required to separate degenerate eigenvalues
         alg = MPSKit.Defaults.alg_eigsolve(; dynamic_tols = false, tol)
@@ -1059,7 +1062,7 @@ module TestAlgorithms
         E₀ = expectation_value(gs, H)
         @test E₀ ≈ first(vals_dense[one(U1Irrep)])
 
-        for (sector, vals) in vals_dense
+        for (sector, vals) in pairs(vals_dense)
             # ED tests
             num = length(vals)
             E₀s, ψ₀s, info = exact_diagonalization(H; num, sector, alg)
@@ -1082,8 +1085,11 @@ module TestAlgorithms
         # so effectively shifting the charges by -1
         H_shift = MPSKit.add_physical_charge(H, U1Irrep.([1, 0, 0, 0]))
         H_shift_dense = convert(TensorMap, H_shift)
-        vals_shift_dense = TensorKit.SectorDict(c => sort(v; by = real) for (c, v) in eigvals(H_shift_dense))
-        for (sector, vals) in vals_dense
+        vals_shift_dense = eigvals(H_shift_dense)
+        for v in values(vals_shift_dense)
+            sort!(v; by = real)
+        end
+        for (sector, vals) in pairs(vals_dense)
             sector′ = only(sector ⊗ U1Irrep(-1))
             @test vals ≈ vals_shift_dense[sector′]
         end