diff --git a/Project.toml b/Project.toml
index 51e51efdf..f4fe33394 100644
--- a/Project.toml
+++ b/Project.toml
@@ -42,7 +42,7 @@ HalfIntegers = "1.6.0"
 KrylovKit = "0.8.3, 0.9.2, 0.10"
 LinearAlgebra = "1.6"
 LoggingExtras = "~1.0"
-MatrixAlgebraKit = "0.6"
+MatrixAlgebraKit = "0.6.5"
 OhMyThreads = "0.7, 0.8"
 OptimKit = "0.3.1, 0.4"
 ParallelTestRunner = "2"
diff --git a/src/MPSKit.jl b/src/MPSKit.jl
index 30db5aa89..c1b47dc40 100644
--- a/src/MPSKit.jl
+++ b/src/MPSKit.jl
@@ -62,7 +62,7 @@ using TensorKit
 using TensorKit: BraidingTensor
 using TensorKit: TupleTools as TT
 using MatrixAlgebraKit
-using MatrixAlgebraKit: TruncationStrategy, PolarViaSVD, LAPACK_SVDAlgorithm
+using MatrixAlgebraKit: TruncationStrategy
 using BlockTensorKit
 using BlockTensorKit: TensorMapSumSpace
 using TensorOperations
diff --git a/src/algorithms/approximate/vomps.jl b/src/algorithms/approximate/vomps.jl
index 224d77486..8893858c6 100644
--- a/src/algorithms/approximate/vomps.jl
+++ b/src/algorithms/approximate/vomps.jl
@@ -65,7 +65,7 @@ end
 function localupdate_step!(
         ::IterativeSolver{<:VOMPS}, state::VOMPSState{<:Any, <:Tuple}, ::SerialScheduler
     )
-    alg_orth = Defaults.alg_qr()
+    alg_orth = Defaults.alg_orth()
 
     ACs = similar(state.mps.AC)
     dst_ACs = state.mps isa Multiline ? eachcol(ACs) : ACs
@@ -88,7 +88,7 @@ end
 function localupdate_step!(
         ::IterativeSolver{<:VOMPS}, state::VOMPSState{<:Any, <:Tuple}, scheduler
     )
-    alg_orth = Defaults.alg_qr()
+    alg_orth = Defaults.alg_orth()
 
     ACs = similar(state.mps.AC)
     dst_ACs = state.mps isa Multiline ? eachcol(ACs) : ACs
diff --git a/src/algorithms/changebonds/svdcut.jl b/src/algorithms/changebonds/svdcut.jl
index d5d4ea48a..379abcba8 100644
--- a/src/algorithms/changebonds/svdcut.jl
+++ b/src/algorithms/changebonds/svdcut.jl
@@ -122,11 +122,11 @@ function changebonds!(H::FiniteMPOHamiltonian, alg::SvdCut)
     # 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_trunc)
+        H = left_canonicalize!(H, i; alg = MatrixAlgebraKit.LeftOrthViaSVD(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_trunc)
+        H = right_canonicalize!(H, i; alg = MatrixAlgebraKit.RightOrthViaSVD(alg_trunc))
     end
 
     return H
diff --git a/src/algorithms/groundstate/vumps.jl b/src/algorithms/groundstate/vumps.jl
index 47c638db9..2c4054493 100644
--- a/src/algorithms/groundstate/vumps.jl
+++ b/src/algorithms/groundstate/vumps.jl
@@ -106,7 +106,7 @@ function localupdate_step!(
         it::IterativeSolver{<:VUMPS}, state, scheduler = Defaults.scheduler[]
     )
     alg_eigsolve = updatetol(it.alg_eigsolve, state.iter, state.ϵ)
-    alg_orth = Defaults.alg_qr()
+    alg_orth = Defaults.alg_orth()
 
     mps = state.mps
     src_Cs = mps isa Multiline ? eachcol(mps.C) : mps.C
@@ -127,7 +127,7 @@ end
 
 function _localupdate_vumps_step!(
         site, mps, operator, envs, AC₀, C₀;
-        parallel::Bool = false, alg_orth = Defaults.alg_qr(),
+        parallel::Bool = false, alg_orth = Defaults.alg_orth(),
         alg_eigsolve = Defaults.eigsolver, which
     )
     if !parallel
diff --git a/src/algorithms/statmech/vomps.jl b/src/algorithms/statmech/vomps.jl
index ed738ebf5..866d2d75b 100644
--- a/src/algorithms/statmech/vomps.jl
+++ b/src/algorithms/statmech/vomps.jl
@@ -102,7 +102,7 @@ end
 function localupdate_step!(
         ::IterativeSolver{<:VOMPS}, state, scheduler = Defaults.scheduler[]
     )
-    alg_orth = Defaults.alg_qr()
+    alg_orth = Defaults.alg_orth()
     mps = state.mps
     ACs = similar(mps.AC)
     dst_ACs = state.mps isa Multiline ? eachcol(ACs) : ACs
@@ -119,7 +119,7 @@ function localupdate_step!(
 end
 
 function _localupdate_vomps_step!(
-        site, mps, operator, envs; parallel::Bool = false, alg_orth = Defaults.alg_qr()
+        site, mps, operator, envs; parallel::Bool = false, alg_orth = Defaults.alg_orth()
     )
     if !parallel
         AC = AC_projection(site, mps, operator, mps, envs)
diff --git a/src/operators/ortho.jl b/src/operators/ortho.jl
index d21985797..b48d12ef1 100644
--- a/src/operators/ortho.jl
+++ b/src/operators/ortho.jl
@@ -1,6 +1,6 @@
 function left_canonicalize!(
         H::FiniteMPOHamiltonian, i::Int;
-        alg::MatrixAlgebraKit.AbstractAlgorithm = Defaults.alg_qr()
+        alg::MatrixAlgebraKit.AbstractAlgorithm = Defaults.alg_orth()
     )
     1 ≤ i < length(H) || throw(ArgumentError("Bounds error in canonicalize"))
 
@@ -88,7 +88,7 @@ end
 
 function right_canonicalize!(
         H::FiniteMPOHamiltonian, i::Int;
-        alg::MatrixAlgebraKit.AbstractAlgorithm = Defaults.alg_lq()
+        alg::MatrixAlgebraKit.AbstractAlgorithm = Defaults.alg_orth()
     )
     1 < i ≤ length(H) || throw(ArgumentError("Bounds error in canonicalize"))
 
diff --git a/src/states/abstractmps.jl b/src/states/abstractmps.jl
index cbc019dfb..40785968b 100644
--- a/src/states/abstractmps.jl
+++ b/src/states/abstractmps.jl
@@ -110,20 +110,20 @@ function isfullrank(V::TensorKit.TensorMapSpace; side = :both)
 end
 
 """
-    makefullrank!(A::PeriodicVector{<:GenericMPSTensor}; alg=Defaults.alg_qr())
+    makefullrank!(A::PeriodicVector{<:GenericMPSTensor}; alg=Defaults.alg_orth())
 
 Make the set of MPS tensors full rank by performing a series of orthogonalizations.
 """
-function makefullrank!(A::PeriodicVector{<:GenericMPSTensor}; alg_leftorth = Defaults.alg_qr(), alg_rightorth = Defaults.alg_lq())
+function makefullrank!(A::PeriodicVector{<:GenericMPSTensor}; alg_orth = Defaults.alg_orth())
     while true
         i = findfirst(!isfullrank, A)
         isnothing(i) && break
         if !isfullrank(A[i]; side = :left)
-            L, Q = right_orth!(_transpose_tail(A[i]); alg = alg_rightorth)
+            L, Q = right_orth!(_transpose_tail(A[i]); alg = alg_orth)
             A[i] = _transpose_front(Q)
             A[i - 1] = A[i - 1] * L
         else
-            A[i], R = left_orth!(A[i]; alg = alg_leftorth)
+            A[i], R = left_orth!(A[i]; alg = alg_orth)
             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 a8158a325..d1646c7fd 100644
--- a/src/states/ortho.jl
+++ b/src/states/ortho.jl
@@ -21,7 +21,7 @@ $(TYPEDFIELDS)
     verbosity::Int = VERBOSE_WARN
 
     "algorithm used for orthogonalization of the tensors"
-    alg_orth = Defaults.alg_qr()
+    alg_orth = Defaults.alg_orth()
     "algorithm used for the eigensolver"
     alg_eigsolve = _GAUGE_ALG_EIGSOLVE
     "minimal amount of iterations before using the eigensolver steps"
@@ -46,7 +46,7 @@ $(TYPEDFIELDS)
     verbosity::Int = VERBOSE_WARN
 
     "algorithm used for orthogonalization of the tensors"
-    alg_orth = Defaults.alg_lq()
+    alg_orth = Defaults.alg_orth()
     "algorithm used for the eigensolver"
     alg_eigsolve = _GAUGE_ALG_EIGSOLVE
     "minimal amount of iterations before using the eigensolver steps"
@@ -73,31 +73,15 @@ end
 
 function MixedCanonical(;
         tol::Real = Defaults.tolgauge, maxiter::Int = Defaults.maxiter,
-        verbosity::Int = VERBOSE_WARN, alg_orth = Defaults.alg_qr(),
+        verbosity::Int = VERBOSE_WARN, alg_orth = Defaults.alg_orth(),
         alg_eigsolve = _GAUGE_ALG_EIGSOLVE,
         eig_miniter::Int = 10, order::Symbol = :LR
     )
-    if alg_orth isa LAPACK_HouseholderQR
-        alg_leftorth = alg_orth
-        alg_rightorth = LAPACK_HouseholderLQ(; alg_orth.kwargs...)
-    elseif alg_orth isa CUSOLVER_HouseholderQR
-        alg_leftorth = alg_orth
-        alg_rightorth = LQViaTransposedQR(CUSOLVER_HouseholderQR(; alg_orth.kwargs...))
-    elseif alg_orth isa LAPACK_HouseholderLQ
-        alg_leftorth = LAPACK_HouseholderQR(; alg_orth.kwargs...)
-        alg_rightorth = alg_orth
-    elseif alg_orth isa LQViaTransposedQR
-        alg_leftorth = alg_orth
-        alg_rightorth = alg_orth.qr_alg
-    else
-        alg_leftorth = alg_rightorth = alg_orth
-    end
-
     left = LeftCanonical(;
-        tol, maxiter, verbosity, alg_orth = alg_leftorth, alg_eigsolve, eig_miniter
+        tol, maxiter, verbosity, alg_orth, alg_eigsolve, eig_miniter
     )
     right = RightCanonical(;
-        tol, maxiter = maxiter, verbosity, alg_orth = alg_rightorth, alg_eigsolve, eig_miniter
+        tol, maxiter = maxiter, verbosity, alg_orth, alg_eigsolve, eig_miniter
     )
 
     return MixedCanonical(left, right, order)
@@ -168,7 +152,7 @@ performance for accuracy.
 regauge!
 
 function regauge!(
-        AC::GenericMPSTensor, C::MPSBondTensor; alg = Defaults.alg_qr()
+        AC::GenericMPSTensor, C::MPSBondTensor; alg = Defaults.alg_orth()
     )
     Q_AC, _ = left_orth!(AC; alg)
     Q_C, _ = left_orth!(C; alg)
@@ -181,7 +165,7 @@ function regauge!(AC::AbstractVector{<:GenericMPSTensor}, C::AbstractVector{<:MP
     return AC
 end
 function regauge!(
-        CL::MPSBondTensor, AC::GenericMPSTensor; alg = Defaults.alg_lq()
+        CL::MPSBondTensor, AC::GenericMPSTensor; alg = Defaults.alg_orth()
     )
     AC_tail = _transpose_tail(AC)
     _, Q_AC = right_orth!(AC_tail; alg)
diff --git a/src/utility/defaults.jl b/src/utility/defaults.jl
index ba1c5aeec..b49fa1719 100644
--- a/src/utility/defaults.jl
+++ b/src/utility/defaults.jl
@@ -9,7 +9,7 @@ import KrylovKit: GMRES, Arnoldi, Lanczos
 using OhMyThreads
 using ..MPSKit: DynamicTol
 using TensorKit: TensorKit
-using MatrixAlgebraKit: LAPACK_HouseholderQR, LAPACK_HouseholderLQ, LAPACK_DivideAndConquer
+using MatrixAlgebraKit: DefaultAlgorithm, Householder
 
 const VERBOSE_NONE = 0
 const VERBOSE_WARN = 1
@@ -57,9 +57,8 @@ function alg_eigsolve(;
     return dynamic_tols ? DynamicTol(alg, tol_min, tol_max, tol_factor) : alg
 end
 
-alg_svd() = LAPACK_DivideAndConquer()
-alg_qr() = LAPACK_HouseholderQR(; positive = true)
-alg_lq() = LAPACK_HouseholderLQ(; positive = true)
+alg_svd() = DefaultAlgorithm()
+alg_orth() = Householder(; positive = true)
 
 # TODO: make verbosity and maxiter actually do something
 function alg_environments(;
diff --git a/test/cuda/operators.jl b/test/cuda/operators.jl
index b3f03796c..cf42a4e7b 100644
--- a/test/cuda/operators.jl
+++ b/test/cuda/operators.jl
@@ -7,10 +7,6 @@ using MatrixAlgebraKit
 using TensorKit: ℙ, tensormaptype, TensorMapWithStorage
 using Adapt, CUDA, cuTENSOR
 
-# TODO revisit this once https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/issues/176
-# is resolved
-MPSKit.Defaults.alg_svd() = CUSOLVER_QRIteration()
-
 @testset "CuFiniteMPO" for V in (ℂ^2, U1Space(0 => 1, 1 => 1))
     # start from random operators
     L = 4