diff --git a/Project.toml b/Project.toml
index 479074d2..ec99e515 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "MatrixAlgebraKit"
 uuid = "6c742aac-3347-4629-af66-fc926824e5e4"
 authors = ["Jutho <jutho.haegeman@ugent.be> and contributors"]
-version = "0.2.4"
+version = "0.2.5"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
diff --git a/src/MatrixAlgebraKit.jl b/src/MatrixAlgebraKit.jl
index 54207e60..567c1510 100644
--- a/src/MatrixAlgebraKit.jl
+++ b/src/MatrixAlgebraKit.jl
@@ -9,6 +9,8 @@ using LinearAlgebra: Diagonal, diag, diagind
 using LinearAlgebra: UpperTriangular, LowerTriangular
 using LinearAlgebra: BlasFloat, BlasReal, BlasComplex, BlasInt, triu!, tril!, rdiv!, ldiv!
 
+export isisometry, isunitary
+
 export qr_compact, qr_full, qr_null, lq_compact, lq_full, lq_null
 export qr_compact!, qr_full!, qr_null!, lq_compact!, lq_full!, lq_null!
 export svd_compact, svd_full, svd_vals, svd_trunc
@@ -40,6 +42,7 @@ include("common/pullbacks.jl")
 include("common/safemethods.jl")
 include("common/view.jl")
 include("common/regularinv.jl")
+include("common/matrixproperties.jl")
 
 include("yalapack.jl")
 include("algorithms.jl")
diff --git a/src/common/matrixproperties.jl b/src/common/matrixproperties.jl
new file mode 100644
index 00000000..e6a6494e
--- /dev/null
+++ b/src/common/matrixproperties.jl
@@ -0,0 +1,59 @@
+"""
+    isisometry(A; side=:left, isapprox_kwargs...) -> Bool
+
+Test whether a linear map is an isometry, where the type of isometry is controlled by `kind`:
+
+- `side = :left` : `A' * A ≈ I`. 
+- `side = :right` : `A * A' ≈ I`.
+
+The `isapprox_kwargs` are passed on to `isapprox` to control the tolerances.
+
+New specializations should overload [`is_left_isometry`](@ref) and [`is_right_isometry`](@ref).
+
+See also [`isunitary`](@ref).
+"""
+function isisometry(A; side::Symbol=:left, isapprox_kwargs...)
+    side === :left && return is_left_isometry(A; isapprox_kwargs...)
+    side === :right && return is_right_isometry(A; isapprox_kwargs...)
+
+    throw(ArgumentError(lazy"Invalid isometry side: $side"))
+end
+
+"""
+    isunitary(A; isapprox_kwargs...)
+
+Test whether a linear map is unitary, i.e. `A * A' ≈ I ≈ A' * A`.
+The `isapprox_kwargs` are passed on to `isapprox` to control the tolerances.
+
+See also [`isisometry`](@ref).
+"""
+function isunitary(A; isapprox_kwargs...)
+    return is_left_isometry(A; isapprox_kwargs...) &&
+           is_right_isometry(A; isapprox_kwargs...)
+end
+
+@doc """
+    is_left_isometry(A; isapprox_kwargs...) -> Bool
+
+Test whether a linear map is a left isometry, i.e. `A' * A ≈ I`.
+The `isapprox_kwargs` can be used to control the tolerances of the equality.
+
+See also [`isisometry`](@ref) and [`is_right_isometry`](@ref).
+""" is_left_isometry
+
+function is_left_isometry(A::AbstractMatrix; isapprox_kwargs...)
+    return isapprox(A' * A, LinearAlgebra.I; isapprox_kwargs...)
+end
+
+@doc """
+    is_right_isometry(A; isapprox_kwargs...) -> Bool
+
+Test whether a linear map is a right isometry, i.e. `A * A' ≈ I`.
+The `isapprox_kwargs` can be used to control the tolerances of the equality.
+
+See also [`isisometry`](@ref) and [`is_left_isometry`](@ref).
+""" is_right_isometry
+
+function is_right_isometry(A::AbstractMatrix; isapprox_kwargs...)
+    return isapprox(A * A', LinearAlgebra.I; isapprox_kwargs...)
+end
diff --git a/test/eigh.jl b/test/eigh.jl
index 6e785158..26e6ce8a 100644
--- a/test/eigh.jl
+++ b/test/eigh.jl
@@ -17,8 +17,7 @@ using MatrixAlgebraKit: TruncatedAlgorithm, diagview
 
         D, V = @constinferred eigh_full(A; alg)
         @test A * V ≈ V * D
-        @test V' * V ≈ I
-        @test V * V' ≈ I
+        @test isunitary(V)
         @test all(isreal, D)
 
         D2, V2 = eigh_full!(copy(A), (D, V), alg)
@@ -47,14 +46,14 @@ end
 
         D1, V1 = @constinferred eigh_trunc(A; alg, trunc=truncrank(r))
         @test length(diagview(D1)) == r
-        @test V1' * V1 ≈ I
+        @test isisometry(V1)
         @test A * V1 ≈ V1 * D1
         @test LinearAlgebra.opnorm(A - V1 * D1 * V1') ≈ D₀[r + 1]
 
         trunc = trunctol(s * D₀[r + 1])
         D2, V2 = @constinferred eigh_trunc(A; alg, trunc)
         @test length(diagview(D2)) == r
-        @test V2' * V2 ≈ I
+        @test isisometry(V2)
         @test A * V2 ≈ V2 * D2
 
         # test for same subspace
diff --git a/test/lq.jl b/test/lq.jl
index 10375e50..f12ea99a 100644
--- a/test/lq.jl
+++ b/test/lq.jl
@@ -14,11 +14,11 @@ using LinearAlgebra: diag, I
         @test L isa Matrix{T} && size(L) == (m, minmn)
         @test Q isa Matrix{T} && size(Q) == (minmn, n)
         @test L * Q ≈ A
-        @test Q * Q' ≈ I
+        @test isisometry(Q; side=:right)
         Nᴴ = @constinferred lq_null(A)
         @test Nᴴ isa Matrix{T} && size(Nᴴ) == (n - minmn, n)
         @test maximum(abs, A * Nᴴ') < eps(real(T))^(2 / 3)
-        @test Nᴴ * Nᴴ' ≈ I
+        @test isisometry(Nᴴ; side=:right)
 
         Ac = similar(A)
         L2, Q2 = @constinferred lq_compact!(copy!(Ac, A), (L, Q))
@@ -35,12 +35,12 @@ using LinearAlgebra: diag, I
         # unblocked algorithm
         lq_compact!(copy!(Ac, A), (L, Q); blocksize=1)
         @test L * Q ≈ A
-        @test Q * Q' ≈ I
+        @test isisometry(Q; side=:right)
         lq_compact!(copy!(Ac, A), (noL, Q2); blocksize=1)
         @test Q == Q2
         lq_null!(copy!(Ac, A), Nᴴ; blocksize=1)
         @test maximum(abs, A * Nᴴ') < eps(real(T))^(2 / 3)
-        @test Nᴴ * Nᴴ' ≈ I
+        @test isisometry(Nᴴ; side=:right)
         if m <= n
             lq_compact!(copy!(Q2, A), (noL, Q2); blocksize=1) # in-place Q
             @test Q ≈ Q2
@@ -50,25 +50,25 @@ using LinearAlgebra: diag, I
         end
         lq_compact!(copy!(Ac, A), (L, Q); blocksize=8)
         @test L * Q ≈ A
-        @test Q * Q' ≈ I
+        @test isisometry(Q; side=:right)
         lq_compact!(copy!(Ac, A), (noL, Q2); blocksize=8)
         @test Q == Q2
         lq_null!(copy!(Ac, A), Nᴴ; blocksize=8)
         @test maximum(abs, A * Nᴴ') < eps(real(T))^(2 / 3)
-        @test Nᴴ * Nᴴ' ≈ I
+        @test isisometry(Nᴴ; side=:right)
         # pivoted
         @test_throws ArgumentError lq_compact!(copy!(Ac, A), (L, Q); pivoted=true)
         # positive
         lq_compact!(copy!(Ac, A), (L, Q); positive=true)
         @test L * Q ≈ A
-        @test Q * Q' ≈ I
+        @test isisometry(Q; side=:right)
         @test all(>=(zero(real(T))), real(diag(L)))
         lq_compact!(copy!(Ac, A), (noL, Q2); positive=true)
         @test Q == Q2
         # positive and blocksize 1
         lq_compact!(copy!(Ac, A), (L, Q); positive=true, blocksize=1)
         @test L * Q ≈ A
-        @test Q * Q' ≈ I
+        @test isisometry(Q; side=:right)
         @test all(>=(zero(real(T))), real(diag(L)))
         lq_compact!(copy!(Ac, A), (noL, Q2); positive=true, blocksize=1)
         @test Q == Q2
@@ -85,14 +85,14 @@ end
         @test L isa Matrix{T} && size(L) == (m, n)
         @test Q isa Matrix{T} && size(Q) == (n, n)
         @test L * Q ≈ A
-        @test Q * Q' ≈ I
+        @test isunitary(Q)
 
         Ac = similar(A)
         L2, Q2 = @constinferred lq_full!(copy!(Ac, A), (L, Q))
         @test L2 === L
         @test Q2 === Q
         @test L * Q ≈ A
-        @test Q * Q' ≈ I
+        @test isunitary(Q)
 
         noL = similar(A, 0, n)
         Q2 = similar(Q)
@@ -102,7 +102,7 @@ end
         # unblocked algorithm
         lq_full!(copy!(Ac, A), (L, Q); blocksize=1)
         @test L * Q ≈ A
-        @test Q * Q' ≈ I
+        @test isunitary(Q)
         lq_full!(copy!(Ac, A), (noL, Q2); blocksize=1)
         @test Q == Q2
         if n == m
@@ -112,7 +112,7 @@ end
         # # other blocking
         lq_full!(copy!(Ac, A), (L, Q); blocksize=18)
         @test L * Q ≈ A
-        @test Q * Q' ≈ I
+        @test isunitary(Q)
         lq_full!(copy!(Ac, A), (noL, Q2); blocksize=18)
         @test Q == Q2
         # pivoted
@@ -120,14 +120,14 @@ end
         # positive
         lq_full!(copy!(Ac, A), (L, Q); positive=true)
         @test L * Q ≈ A
-        @test Q * Q' ≈ I
+        @test isunitary(Q)
         @test all(>=(zero(real(T))), real(diag(L)))
         lq_full!(copy!(Ac, A), (noL, Q2); positive=true)
         @test Q == Q2
         # positive and blocksize 1
         lq_full!(copy!(Ac, A), (L, Q); positive=true, blocksize=1)
         @test L * Q ≈ A
-        @test Q * Q' ≈ I
+        @test isunitary(Q)
         @test all(>=(zero(real(T))), real(diag(L)))
         lq_full!(copy!(Ac, A), (noL, Q2); positive=true, blocksize=1)
         @test Q == Q2
diff --git a/test/orthnull.jl b/test/orthnull.jl
index b0004739..cc8f2a19 100644
--- a/test/orthnull.jl
+++ b/test/orthnull.jl
@@ -63,9 +63,9 @@ end
         @test C isa Matrix{T} && size(C) == (minmn, n)
         @test N isa Matrix{T} && size(N) == (m, m - minmn)
         @test V * C ≈ A
-        @test V' * V ≈ I
+        @test isisometry(V)
         @test LinearAlgebra.norm(A' * N) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T)
-        @test N' * N ≈ I
+        @test isisometry(N)
         @test V * V' + N * N' ≈ I
 
         M = LinearMap(A)
@@ -80,9 +80,9 @@ end
             @test C isa Matrix{T} && size(C) == (minmn, n)
             @test N isa Matrix{T} && size(N) == (m, nullity)
             @test V * C ≈ A
-            @test V' * V ≈ I
+            @test isisometry(V)
             @test LinearAlgebra.norm(A' * N) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T)
-            @test N' * N ≈ I
+            @test isisometry(N)
         end
 
         for alg_qr in ((; positive=true), (; positive=false), LAPACK_HouseholderQR())
@@ -92,9 +92,9 @@ end
             @test C isa Matrix{T} && size(C) == (minmn, n)
             @test N isa Matrix{T} && size(N) == (m, m - minmn)
             @test V * C ≈ A
-            @test V' * V ≈ I
+            @test isisometry(V)
             @test LinearAlgebra.norm(A' * N) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T)
-            @test N' * N ≈ I
+            @test isisometry(N)
             @test V * V' + N * N' ≈ I
         end
 
@@ -105,9 +105,9 @@ end
         @test C2 === C
         @test N2 === N
         @test V2 * C2 ≈ A
-        @test V2' * V2 ≈ I
+        @test isisometry(V2)
         @test LinearAlgebra.norm(A' * N2) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T)
-        @test N2' * N2 ≈ I
+        @test isisometry(N2)
         @test V2 * V2' + N2 * N2' ≈ I
 
         atol = eps(real(T))
@@ -117,9 +117,9 @@ end
         @test C2 !== C
         @test N2 !== C
         @test V2 * C2 ≈ A
-        @test V2' * V2 ≈ I
+        @test isisometry(V2)
         @test LinearAlgebra.norm(A' * N2) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T)
-        @test N2' * N2 ≈ I
+        @test isisometry(N2)
         @test V2 * V2' + N2 * N2' ≈ I
 
         rtol = eps(real(T))
@@ -131,9 +131,9 @@ end
             @test C2 !== C
             @test N2 !== C
             @test V2 * C2 ≈ A
-            @test V2' * V2 ≈ I
+            @test isisometry(V2)
             @test LinearAlgebra.norm(A' * N2) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T)
-            @test N2' * N2 ≈ I
+            @test isisometry(N2)
             @test V2 * V2' + N2 * N2' ≈ I
         end
 
@@ -143,12 +143,12 @@ end
             @test V2 === V
             @test C2 === C
             @test V2 * C2 ≈ A
-            @test V2' * V2 ≈ I
+            @test isisometry(V2)
             if kind != :polar
                 N2 = @constinferred left_null!(copy!(Ac, A), N; kind=kind)
                 @test N2 === N
                 @test LinearAlgebra.norm(A' * N2) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T)
-                @test N2' * N2 ≈ I
+                @test isisometry(N2)
                 @test V2 * V2' + N2 * N2' ≈ I
             end
 
@@ -175,9 +175,9 @@ end
                 @test C2 !== C
                 @test N2 !== C
                 @test V2 * C2 ≈ A
-                @test V2' * V2 ≈ I
+                @test isisometry(V2)
                 @test LinearAlgebra.norm(A' * N2) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T)
-                @test N2' * N2 ≈ I
+                @test isisometry(N2)
                 @test V2 * V2' + N2 * N2' ≈ I
             else
                 @test_throws ArgumentError left_orth!(copy!(Ac, A), (V, C); kind=kind,
@@ -206,9 +206,9 @@ end
         @test Vᴴ isa Matrix{T} && size(Vᴴ) == (minmn, n)
         @test Nᴴ isa Matrix{T} && size(Nᴴ) == (n - minmn, n)
         @test C * Vᴴ ≈ A
-        @test Vᴴ * Vᴴ' ≈ I
+        @test isisometry(Vᴴ; side=:right)
         @test LinearAlgebra.norm(A * adjoint(Nᴴ)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T)
-        @test Nᴴ * Nᴴ' ≈ I
+        @test isisometry(Nᴴ; side=:right)
         @test Vᴴ' * Vᴴ + Nᴴ' * Nᴴ ≈ I
 
         M = LinearMap(A)
@@ -222,9 +222,9 @@ end
         @test Vᴴ2 === Vᴴ
         @test Nᴴ2 === Nᴴ
         @test C2 * Vᴴ2 ≈ A
-        @test Vᴴ2 * Vᴴ2' ≈ I
+        @test isisometry(Vᴴ2; side=:right)
         @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T)
-        @test Nᴴ2 * Nᴴ2' ≈ I
+        @test isisometry(Nᴴ; side=:right)
         @test Vᴴ2' * Vᴴ2 + Nᴴ2' * Nᴴ2 ≈ I
 
         atol = eps(real(T))
@@ -234,9 +234,9 @@ end
         @test Vᴴ2 !== Vᴴ
         @test Nᴴ2 !== Nᴴ
         @test C2 * Vᴴ2 ≈ A
-        @test Vᴴ2 * Vᴴ2' ≈ I
+        @test isisometry(Vᴴ2; side=:right)
         @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T)
-        @test Nᴴ2 * Nᴴ2' ≈ I
+        @test isisometry(Nᴴ; side=:right)
         @test Vᴴ2' * Vᴴ2 + Nᴴ2' * Nᴴ2 ≈ I
 
         rtol = eps(real(T))
@@ -246,9 +246,9 @@ end
         @test Vᴴ2 !== Vᴴ
         @test Nᴴ2 !== Nᴴ
         @test C2 * Vᴴ2 ≈ A
-        @test Vᴴ2 * Vᴴ2' ≈ I
+        @test isisometry(Vᴴ2; side=:right)
         @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T)
-        @test Nᴴ2 * Nᴴ2' ≈ I
+        @test isisometry(Nᴴ2; side=:right)
         @test Vᴴ2' * Vᴴ2 + Nᴴ2' * Nᴴ2 ≈ I
 
         for kind in (:lq, :polar, :svd)
@@ -257,12 +257,12 @@ end
             @test C2 === C
             @test Vᴴ2 === Vᴴ
             @test C2 * Vᴴ2 ≈ A
-            @test Vᴴ2 * Vᴴ2' ≈ I
+            @test isisometry(Vᴴ2; side=:right)
             if kind != :polar
                 Nᴴ2 = @constinferred right_null!(copy!(Ac, A), Nᴴ; kind=kind)
                 @test Nᴴ2 === Nᴴ
                 @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T)
-                @test Nᴴ2 * Nᴴ2' ≈ I
+                @test isisometry(Nᴴ2; side=:right)
                 @test Vᴴ2' * Vᴴ2 + Nᴴ2' * Nᴴ2 ≈ I
             end
 
@@ -275,9 +275,9 @@ end
                 @test Vᴴ2 !== Vᴴ
                 @test Nᴴ2 !== Nᴴ
                 @test C2 * Vᴴ2 ≈ A
-                @test Vᴴ2 * Vᴴ2' ≈ I
+                @test isisometry(Vᴴ2; side=:right)
                 @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T)
-                @test Nᴴ2 * Nᴴ2' ≈ I
+                @test isisometry(Nᴴ2; side=:right)
                 @test Vᴴ2' * Vᴴ2 + Nᴴ2' * Nᴴ2 ≈ I
 
                 C2, Vᴴ2 = @constinferred right_orth!(copy!(Ac, A), (C, Vᴴ); kind=kind,
@@ -288,9 +288,9 @@ end
                 @test Vᴴ2 !== Vᴴ
                 @test Nᴴ2 !== Nᴴ
                 @test C2 * Vᴴ2 ≈ A
-                @test Vᴴ2 * Vᴴ2' ≈ I
+                @test isisometry(Vᴴ2; side=:right)
                 @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T)
-                @test Nᴴ2 * Nᴴ2' ≈ I
+                @test isisometry(Nᴴ2; side=:right)
                 @test Vᴴ2' * Vᴴ2 + Nᴴ2' * Nᴴ2 ≈ I
             else
                 @test_throws ArgumentError right_orth!(copy!(Ac, A), (C, Vᴴ); kind=kind,
diff --git a/test/polar.jl b/test/polar.jl
index 513654ea..b7512d28 100644
--- a/test/polar.jl
+++ b/test/polar.jl
@@ -24,7 +24,7 @@ using MatrixAlgebraKit: PolarViaSVD
             @test W isa Matrix{T} && size(W) == (m, n)
             @test P isa Matrix{T} && size(P) == (n, n)
             @test W * P ≈ A
-            @test W' * W ≈ I
+            @test isisometry(W)
             @test isposdef(P)
 
             Ac = similar(A)
@@ -32,7 +32,7 @@ using MatrixAlgebraKit: PolarViaSVD
             @test W2 === W
             @test P2 === P
             @test W * P ≈ A
-            @test W' * W ≈ I
+            @test isisometry(W)
             @test isposdef(P)
         end
     end
@@ -52,7 +52,7 @@ end
             @test Wᴴ isa Matrix{T} && size(Wᴴ) == (m, n)
             @test P isa Matrix{T} && size(P) == (m, m)
             @test P * Wᴴ ≈ A
-            @test Wᴴ * Wᴴ' ≈ I
+            @test isisometry(Wᴴ; side=:right)
             @test isposdef(P)
 
             Ac = similar(A)
@@ -60,7 +60,7 @@ end
             @test P2 === P
             @test Wᴴ2 === Wᴴ
             @test P * Wᴴ ≈ A
-            @test Wᴴ * Wᴴ' ≈ I
+            @test isisometry(Wᴴ; side=:right)
             @test isposdef(P)
         end
     end
diff --git a/test/qr.jl b/test/qr.jl
index 81c64f0e..83f3ec6d 100644
--- a/test/qr.jl
+++ b/test/qr.jl
@@ -17,9 +17,9 @@ using LinearAlgebra: diag, I
         @test Q * R ≈ A
         N = @constinferred qr_null(A)
         @test N isa Matrix{T} && size(N) == (m, m - minmn)
-        @test Q' * Q ≈ I
+        @test isisometry(Q)
         @test maximum(abs, A' * N) < eps(real(T))^(2 / 3)
-        @test N' * N ≈ I
+        @test isisometry(N)
 
         Ac = similar(A)
         Q2, R2 = @constinferred qr_compact!(copy!(Ac, A), (Q, R))
@@ -36,13 +36,13 @@ using LinearAlgebra: diag, I
         # unblocked algorithm
         qr_compact!(copy!(Ac, A), (Q, R); blocksize=1)
         @test Q * R ≈ A
-        @test Q' * Q ≈ I
+        @test isisometry(Q)
         qr_compact!(copy!(Ac, A), (Q2, noR); blocksize=1)
         @test Q == Q2
         qr_compact!(copy!(Ac, A), (Q2, noR); blocksize=1)
         qr_null!(copy!(Ac, A), N; blocksize=1)
         @test maximum(abs, A' * N) < eps(real(T))^(2 / 3)
-        @test N' * N ≈ I
+        @test isisometry(N)
         if n <= m
             qr_compact!(copy!(Q2, A), (Q2, noR); blocksize=1) # in-place Q
             @test Q ≈ Q2
@@ -53,12 +53,12 @@ using LinearAlgebra: diag, I
         # other blocking
         qr_compact!(copy!(Ac, A), (Q, R); blocksize=8)
         @test Q * R ≈ A
-        @test Q' * Q ≈ I
+        @test isisometry(Q)
         qr_compact!(copy!(Ac, A), (Q2, noR); blocksize=8)
         @test Q == Q2
         qr_null!(copy!(Ac, A), N; blocksize=8)
         @test maximum(abs, A' * N) < eps(real(T))^(2 / 3)
-        @test N' * N ≈ I
+        @test isisometry(N)
 
         # pivoted
         qr_compact!(copy!(Ac, A), (Q, R); pivoted=true)
@@ -69,21 +69,21 @@ using LinearAlgebra: diag, I
         # positive
         qr_compact!(copy!(Ac, A), (Q, R); positive=true)
         @test Q * R ≈ A
-        @test Q' * Q ≈ I
+        @test isisometry(Q)
         @test all(>=(zero(real(T))), real(diag(R)))
         qr_compact!(copy!(Ac, A), (Q2, noR); positive=true)
         @test Q == Q2
         # positive and blocksize 1
         qr_compact!(copy!(Ac, A), (Q, R); positive=true, blocksize=1)
         @test Q * R ≈ A
-        @test Q' * Q ≈ I
+        @test isisometry(Q)
         @test all(>=(zero(real(T))), real(diag(R)))
         qr_compact!(copy!(Ac, A), (Q2, noR); positive=true, blocksize=1)
         @test Q == Q2
         # positive and pivoted
         qr_compact!(copy!(Ac, A), (Q, R); positive=true, pivoted=true)
         @test Q * R ≈ A
-        @test Q' * Q ≈ I
+        @test isisometry(Q)
         if n <= m
             # the following test tries to find the diagonal element (in order to test positivity)
             # before the column permutation. This only works if all columns have a diagonal
@@ -117,14 +117,14 @@ end
         @test Q2 === Q
         @test R2 === R
         @test Q * R ≈ A
-        @test Q' * Q ≈ I
+        @test isunitary(Q)
         qr_full!(copy!(Ac, A), (Q2, noR))
         @test Q == Q2
 
         # unblocked algorithm
         qr_full!(copy!(Ac, A), (Q, R); blocksize=1)
         @test Q * R ≈ A
-        @test Q' * Q ≈ I
+        @test isunitary(Q)
         qr_full!(copy!(Ac, A), (Q2, noR); blocksize=1)
         @test Q == Q2
         if n == m
@@ -134,33 +134,33 @@ end
         # other blocking
         qr_full!(copy!(Ac, A), (Q, R); blocksize=8)
         @test Q * R ≈ A
-        @test Q' * Q ≈ I
+        @test isunitary(Q)
         qr_full!(copy!(Ac, A), (Q2, noR); blocksize=8)
         @test Q == Q2
         # pivoted
         qr_full!(copy!(Ac, A), (Q, R); pivoted=true)
         @test Q * R ≈ A
-        @test Q' * Q ≈ I
+        @test isunitary(Q)
         qr_full!(copy!(Ac, A), (Q2, noR); pivoted=true)
         @test Q == Q2
         # positive
         qr_full!(copy!(Ac, A), (Q, R); positive=true)
         @test Q * R ≈ A
-        @test Q' * Q ≈ I
+        @test isunitary(Q)
         @test all(>=(zero(real(T))), real(diag(R)))
         qr_full!(copy!(Ac, A), (Q2, noR); positive=true)
         @test Q == Q2
         # positive and blocksize 1
         qr_full!(copy!(Ac, A), (Q, R); positive=true, blocksize=1)
         @test Q * R ≈ A
-        @test Q' * Q ≈ I
+        @test isunitary(Q)
         @test all(>=(zero(real(T))), real(diag(R)))
         qr_full!(copy!(Ac, A), (Q2, noR); positive=true, blocksize=1)
         @test Q == Q2
         # positive and pivoted
         qr_full!(copy!(Ac, A), (Q, R); positive=true, pivoted=true)
         @test Q * R ≈ A
-        @test Q' * Q ≈ I
+        @test isunitary(Q)
         if n <= m
             # the following test tries to find the diagonal element (in order to test positivity)
             # before the column permutation. This only works if all columns have a diagonal
diff --git a/test/schur.jl b/test/schur.jl
index 71404e89..897cf155 100644
--- a/test/schur.jl
+++ b/test/schur.jl
@@ -14,7 +14,7 @@ using LinearAlgebra: I
         TA, Z, vals = @constinferred schur_full(A; alg)
         @test eltype(TA) == eltype(Z) == T
         @test eltype(vals) == Tc
-        @test Z' * Z ≈ I
+        @test isisometry(Z)
         @test A * Z ≈ Z * TA
 
         Ac = similar(A)
diff --git a/test/svd.jl b/test/svd.jl
index 375dc8ea..99dd940f 100644
--- a/test/svd.jl
+++ b/test/svd.jl
@@ -3,7 +3,7 @@ using Test
 using TestExtras
 using StableRNGs
 using LinearAlgebra: LinearAlgebra, Diagonal, I, isposdef
-using MatrixAlgebraKit: TruncatedAlgorithm, TruncationKeepAbove, diagview
+using MatrixAlgebraKit: TruncatedAlgorithm, TruncationKeepAbove, diagview, isisometry
 
 @testset "svd_compact! for T = $T" for T in (Float32, Float64, ComplexF32, ComplexF64)
     rng = StableRNG(123)
@@ -32,8 +32,8 @@ using MatrixAlgebraKit: TruncatedAlgorithm, TruncationKeepAbove, diagview
             @test S isa Diagonal{real(T)} && size(S) == (minmn, minmn)
             @test Vᴴ isa Matrix{T} && size(Vᴴ) == (minmn, n)
             @test U * S * Vᴴ ≈ A
-            @test U' * U ≈ I
-            @test Vᴴ * Vᴴ' ≈ I
+            @test isisometry(U)
+            @test isisometry(Vᴴ; side=:right)
             @test isposdef(S)
 
             Ac = similar(A)
@@ -44,8 +44,8 @@ using MatrixAlgebraKit: TruncatedAlgorithm, TruncationKeepAbove, diagview
             @test S2 === S
             @test V2ᴴ === Vᴴ
             @test U * S * Vᴴ ≈ A
-            @test U' * U ≈ I
-            @test Vᴴ * Vᴴ' ≈ I
+            @test isisometry(U)
+            @test isisometry(Vᴴ; side=:right)
             @test isposdef(S)
 
             Sd = @constinferred svd_vals(A, alg′)
@@ -66,8 +66,8 @@ end
             @test S isa Matrix{real(T)} && size(S) == (m, n)
             @test Vᴴ isa Matrix{T} && size(Vᴴ) == (n, n)
             @test U * S * Vᴴ ≈ A
-            @test U' * U ≈ I
-            @test Vᴴ * Vᴴ' ≈ I
+            @test isunitary(U)
+            @test isunitary(Vᴴ)
             @test all(isposdef, diagview(S))
 
             Ac = similar(A)
@@ -76,9 +76,8 @@ end
             @test S2 === S
             @test V2ᴴ === Vᴴ
             @test U * S * Vᴴ ≈ A
-            @test U' * U ≈ I
-            @test Vᴴ * Vᴴ' ≈ I
-            @test Vᴴ' * Vᴴ ≈ I
+            @test isunitary(U)
+            @test isunitary(Vᴴ)
             @test all(isposdef, diagview(S))
 
             Sc = similar(A, real(T), min(m, n))