Skip to content
Merged
Show file tree
Hide file tree
Changes from all 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
45 changes: 20 additions & 25 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,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...)
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