Add isisometry function and integrate into tests (#38)

* Add `isisometry`

* Use `isisometry` in tests

* export function

* Bump v0.2.5

* Update interface based on discussions

* Update tests accordingly

* Change `kind` to `side`

* Update src/common/matrixproperties.jl

Co-authored-by: Jutho <Jutho@users.noreply.github.com>

---------

Co-authored-by: Jutho <Jutho@users.noreply.github.com>

Author: lkdvos

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

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

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

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,22 +112,22 @@ 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
         @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 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