Skip to content
Merged
Show file tree
Hide file tree
Changes from 5 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion ext/MatrixAlgebraKitChainRulesCoreExt.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
3 changes: 2 additions & 1 deletion src/MatrixAlgebraKit.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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!
Expand Down Expand Up @@ -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")
Expand Down
39 changes: 20 additions & 19 deletions src/common/matrixproperties.jl
Original file line number Diff line number Diff line change
@@ -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`:

Expand All @@ -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
Expand All @@ -25,43 +26,43 @@ 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
See also [`isisometric`](@ref) and [`MatrixAlgebraKit.is_left_isometric`](@ref).
""" is_right_isometric

function is_right_isometry(A::AbstractMatrix; atol::Real = 0, rtol::Real = defaulttol(A), norm = LinearAlgebra.norm)
function is_right_isometric(A::AbstractMatrix; atol::Real = 0, rtol::Real = defaulttol(A), norm = LinearAlgebra.norm)
Comment thread
lkdvos marked this conversation as resolved.
Outdated
P = A * A'
nP = norm(P) # isapprox would use `rtol * max(norm(P), norm(I))`
diagview(P) .-= 1
Expand Down
4 changes: 2 additions & 2 deletions test/amd/eigh.jl
Original file line number Diff line number Diff line change
Expand Up @@ -48,14 +48,14 @@ end

D1, V1 = @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 = @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
Expand Down
4 changes: 2 additions & 2 deletions test/cuda/eigh.jl
Original file line number Diff line number Diff line change
Expand Up @@ -43,14 +43,14 @@ end

D1, V1 = @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 = @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
Expand Down
4 changes: 2 additions & 2 deletions test/eigh.jl
Original file line number Diff line number Diff line change
Expand Up @@ -52,14 +52,14 @@ end

D1, V1 = @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 = @constinferred eigh_trunc(A; alg, trunc)
@test length(diagview(D2)) == r
@test isisometry(V2)
@test isisometric(V2)
@test A * V2 ≈ V2 * D2

s = 1 - sqrt(eps(real(T)))
Expand Down
16 changes: 8 additions & 8 deletions test/lq.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand Down Expand Up @@ -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
Expand All @@ -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)
Expand All @@ -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
Expand Down
Loading
Loading