Merge branch 'main' into jh/truncerr

Author: lkdvos

ext/MatrixAlgebraKitChainRulesCoreExt.jl

@@ -14,7 +14,7 @@ MatrixAlgebraKit.iszerotangent(::AbstractZero) = true
 @non_differentiable MatrixAlgebraKit.select_algorithm(args...)
 @non_differentiable MatrixAlgebraKit.initialize_output(args...)
 @non_differentiable MatrixAlgebraKit.check_input(args...)
-@non_differentiable MatrixAlgebraKit.isisometry(args...)
+@non_differentiable MatrixAlgebraKit.isisometric(args...)
 @non_differentiable MatrixAlgebraKit.isunitary(args...)
 
 function ChainRulesCore.rrule(::typeof(copy_input), f, A)

src/MatrixAlgebraKit.jl

@@ -9,7 +9,7 @@ using LinearAlgebra: Diagonal, diag, diagind, isdiag
 using LinearAlgebra: UpperTriangular, LowerTriangular
 using LinearAlgebra: BlasFloat, BlasReal, BlasComplex, BlasInt
 
-export isisometry, isunitary, ishermitian, isantihermitian
+export isisometric, isunitary, ishermitian, isantihermitian
 
 export project_hermitian, project_antihermitian, project_isometric
 export project_hermitian!, project_antihermitian!, project_isometric!
@@ -65,6 +65,7 @@ export notrunc, truncrank, trunctol, truncerror, truncfilter
             :svd_pullback!, :svd_trunc_pullback!
         )
     )
+    eval(Expr(:public, :is_left_isometric, :is_right_isometric))
 end
 
 include("common/defaults.jl")

src/common/matrixproperties.jl

@@ -1,5 +1,5 @@
 """
-    isisometry(A; side=:left, isapprox_kwargs...) -> Bool
+    isisometric(A; side=:left, isapprox_kwargs...) -> Bool
 
 Test whether a linear map is an isometry, where the type of isometry is controlled by `kind`:
 
@@ -8,13 +8,14 @@ Test whether a linear map is an isometry, where the type of isometry is controll
 
 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).
+New specializations should overload [`MatrixAlgebraKit.is_left_isometric`](@ref) and
+[`MatrixAlgebraKit.is_right_isometric`](@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...)
+function isisometric(A; side::Symbol = :left, isapprox_kwargs...)
+    side === :left && return is_left_isometric(A; isapprox_kwargs...)
+    side === :right && return is_right_isometric(A; isapprox_kwargs...)
 
     throw(ArgumentError(lazy"Invalid isometry side: $side"))
 end
@@ -25,48 +26,42 @@ end
 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).
+See also [`isisometric`](@ref).
 """
 function isunitary(A; isapprox_kwargs...)
-    return is_left_isometry(A; isapprox_kwargs...) &&
-        is_right_isometry(A; isapprox_kwargs...)
+    return is_left_isometric(A; isapprox_kwargs...) &&
+        is_right_isometric(A; isapprox_kwargs...)
 end
 function isunitary(A::AbstractMatrix; isapprox_kwargs...)
     size(A, 1) == size(A, 2) || return false
-    return is_left_isometry(A; isapprox_kwargs...)
+    return is_left_isometric(A; isapprox_kwargs...)
 end
 
 @doc """
-    is_left_isometry(A; isapprox_kwargs...) -> Bool
+    is_left_isometric(A; isapprox_kwargs...) -> Bool
 
-Test whether a linear map is a left isometry, i.e. `A' * A ≈ I`.
+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
+See also [`isisometric`](@ref) and [`MatrixAlgebraKit.is_right_isometric`](@ref).
+""" is_left_isometric
 
-function is_left_isometry(A::AbstractMatrix; atol::Real = 0, rtol::Real = defaulttol(A), norm = LinearAlgebra.norm)
+function is_left_isometric(A::AbstractMatrix; atol::Real = 0, rtol::Real = defaulttol(A), norm = LinearAlgebra.norm)
     P = A' * A
     nP = norm(P) # isapprox would use `rtol * max(norm(P), norm(I))`
     diagview(P) .-= 1
     return norm(P) <= max(atol, rtol * nP) # assume that the norm of I is `sqrt(n)`
 end
 
 @doc """
-    is_right_isometry(A; isapprox_kwargs...) -> Bool
+    is_right_isometric(A; isapprox_kwargs...) -> Bool
 
-Test whether a linear map is a right isometry, i.e. `A * A' ≈ I`.
+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; atol::Real = 0, rtol::Real = defaulttol(A), norm = LinearAlgebra.norm)
-    P = A * A'
-    nP = norm(P) # isapprox would use `rtol * max(norm(P), norm(I))`
-    diagview(P) .-= 1
-    return norm(P) <= max(atol, rtol * nP) # assume that the norm of I is `sqrt(n)`
-end
+See also [`isisometric`](@ref) and [`MatrixAlgebraKit.is_left_isometric`](@ref).
+""" is_right_isometric
+is_right_isometric(A; kwargs...) = is_left_isometric(A'; kwargs...)
 
 """
     ishermitian(A; isapprox_kwargs...)

test/amd/eigh.jl

@@ -48,14 +48,14 @@ end
 
         D1, V1, ϵ1 = @constinferred eigh_trunc(A; alg, trunc=truncrank(r))
         @test length(diagview(D1)) == r
-        @test isisometry(V1)
+        @test isisometric(V1)
         @test A * V1 ≈ V1 * D1
         @test LinearAlgebra.opnorm(A - V1 * D1 * V1') ≈ D₀[r + 1]
 
         trunc = trunctol(; atol=s * D₀[r + 1])
         D2, V2, ϵ2 = @constinferred eigh_trunc(A; alg, trunc)
         @test length(diagview(D2)) == r
-        @test isisometry(V2)
+        @test isisometric(V2)
         @test A * V2 ≈ V2 * D2
 
         # test for same subspace

test/cuda/eigh.jl

@@ -43,14 +43,14 @@ end
 
         D1, V1, ϵ1 = @constinferred eigh_trunc(A; alg, trunc=truncrank(r))
         @test length(diagview(D1)) == r
-        @test isisometry(V1)
+        @test isisometric(V1)
         @test A * V1 ≈ V1 * D1
         @test LinearAlgebra.opnorm(A - V1 * D1 * V1') ≈ D₀[r + 1]
 
         trunc = trunctol(; atol = s * D₀[r + 1])
         D2, V2, ϵ2 = @constinferred eigh_trunc(A; alg, trunc)
         @test length(diagview(D2)) == r
-        @test isisometry(V2)
+        @test isisometric(V2)
         @test A * V2 ≈ V2 * D2
 
         # test for same subspace

test/eigh.jl

@@ -52,15 +52,15 @@ end
 
         D1, V1, ϵ1 = @constinferred eigh_trunc(A; alg, trunc = truncrank(r))
         @test length(diagview(D1)) == r
-        @test isisometry(V1)
+        @test isisometric(V1)
         @test A * V1 ≈ V1 * D1
         @test LinearAlgebra.opnorm(A - V1 * D1 * V1') ≈ D₀[r + 1]
         @test ϵ1 ≈ norm(view(D₀, (r + 1):m))
 
         trunc = trunctol(; atol = s * D₀[r + 1])
         D2, V2, ϵ2 = @constinferred eigh_trunc(A; alg, trunc)
         @test length(diagview(D2)) == r
-        @test isisometry(V2)
+        @test isisometric(V2)
         @test A * V2 ≈ V2 * D2
         @test ϵ2 ≈ norm(view(D₀, (r + 1):m))
 

test/lq.jl

@@ -17,11 +17,11 @@ eltypes = (Float32, Float64, ComplexF32, ComplexF64)
         @test L isa Matrix{T} && size(L) == (m, minmn)
         @test Q isa Matrix{T} && size(Q) == (minmn, n)
         @test L * Q ≈ A
-        @test isisometry(Q; side = :right)
+        @test isisometric(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 isisometry(Nᴴ; side = :right)
+        @test isisometric(Nᴴ; side = :right)
 
         Ac = similar(A)
         L2, Q2 = @constinferred lq_compact!(copy!(Ac, A), (L, Q))
@@ -51,12 +51,12 @@ eltypes = (Float32, Float64, ComplexF32, ComplexF64)
         # unblocked algorithm
         lq_compact!(copy!(Ac, A), (L, Q); blocksize = 1)
         @test L * Q ≈ A
-        @test isisometry(Q; side = :right)
+        @test isisometric(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 isisometry(Nᴴ; side = :right)
+        @test isisometric(Nᴴ; side = :right)
         if m <= n
             lq_compact!(copy!(Q2, A), (noL, Q2); blocksize = 1) # in-place Q
             @test Q ≈ Q2
@@ -66,12 +66,12 @@ eltypes = (Float32, Float64, ComplexF32, ComplexF64)
         end
         lq_compact!(copy!(Ac, A), (L, Q); blocksize = 8)
         @test L * Q ≈ A
-        @test isisometry(Q; side = :right)
+        @test isisometric(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 isisometry(Nᴴ; side = :right)
+        @test isisometric(Nᴴ; side = :right)
         @test Nᴴ * Nᴴ' ≈ I
 
         qr_alg = LAPACK_HouseholderQR(; blocksize = 1)
@@ -91,15 +91,15 @@ eltypes = (Float32, Float64, ComplexF32, ComplexF64)
         # positive
         lq_compact!(copy!(Ac, A), (L, Q); positive = true)
         @test L * Q ≈ A
-        @test isisometry(Q; side = :right)
+        @test isisometric(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 isisometry(Q; side = :right)
+        @test isisometric(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