Add isisometry function and integrate into tests

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"

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
+
 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")

src/common/matrixproperties.jl

@@ -0,0 +1,7 @@
+"""
+    isisometry(A; kwargs...) -> Bool
+
+Test whether a linear map is an isometry, ie `A' * A ≈ I`. 
+The `kwargs` are passed on to `isapprox` to control the tolerances.
+"""
+isisometry(A::AbstractMatrix; kwargs...) = isapprox(A' * A, LinearAlgebra.I; kwargs...)

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 isisometry(V) && isisometry(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

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')
         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ᴴ')
 
         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')
         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ᴴ')
         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')
         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ᴴ')
         # 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')
         @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')
         @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 isisometry(Q) && isisometry(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 isisometry(Q) && isisometry(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 isisometry(Q) && isisometry(Q')
         lq_full!(copy!(Ac, A), (noL, Q2); blocksize=1)
         @test Q == Q2
         if n == m
@@ -112,22 +112,22 @@ end
         # # other blocking
         lq_full!(copy!(Ac, A), (L, Q); blocksize=18)
         @test L * Q ≈ A
-        @test Q * Q' ≈ I
+        @test isisometry(Q) && isisometry(Q')
         lq_full!(copy!(Ac, A), (noL, Q2); blocksize=18)
         @test Q == Q2
         # pivoted
         @test_throws ArgumentError lq_full!(copy!(Ac, A), (L, Q); pivoted=true)
         # positive
         lq_full!(copy!(Ac, A), (L, Q); positive=true)
         @test L * Q ≈ A
-        @test Q * Q' ≈ I
+        @test isisometry(Q) && isisometry(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 isisometry(Q) && isisometry(Q')
         @test all(>=(zero(real(T))), real(diag(L)))
         lq_full!(copy!(Ac, A), (noL, Q2); positive=true, blocksize=1)
         @test Q == Q2

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ᴴ')
         @test LinearAlgebra.norm(A * adjoint(Nᴴ)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T)
-        @test Nᴴ * Nᴴ' ≈ I
+        @test isisometry(Nᴴ')
         @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')
         @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T)
-        @test Nᴴ2 * Nᴴ2' ≈ I
+        @test isisometry(Nᴴ')
         @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')
         @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T)
-        @test Nᴴ2 * Nᴴ2' ≈ I
+        @test isisometry(Nᴴ')
         @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')
         @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T)
-        @test Nᴴ2 * Nᴴ2' ≈ I
+        @test isisometry(Nᴴ2')
         @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')
             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')
                 @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')
                 @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T)
-                @test Nᴴ2 * Nᴴ2' ≈ I
+                @test isisometry(Nᴴ2')
                 @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')
                 @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T)
-                @test Nᴴ2 * Nᴴ2' ≈ I
+                @test isisometry(Nᴴ2')
                 @test Vᴴ2' * Vᴴ2 + Nᴴ2' * Nᴴ2 ≈ I
             else
                 @test_throws ArgumentError right_orth!(copy!(Ac, A), (C, Vᴴ); kind=kind,

test/polar.jl

@@ -24,15 +24,15 @@ 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)
             W2, P2 = @constinferred left_polar!(copy!(Ac, A), (W, P), alg)
             @test W2 === W
             @test P2 === P
             @test W * P ≈ A
-            @test W' * W ≈ I
+            @test isisometry(W)
             @test isposdef(P)
         end
     end
@@ -52,15 +52,15 @@ 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ᴴ')
             @test isposdef(P)
 
             Ac = similar(A)
             P2, Wᴴ2 = @constinferred right_polar!(copy!(Ac, A), (P, Wᴴ), alg)
             @test P2 === P
             @test Wᴴ2 === Wᴴ
             @test P * Wᴴ ≈ A
-            @test Wᴴ * Wᴴ' ≈ I
+            @test isisometry(Wᴴ')
             @test isposdef(P)
         end
     end