diff --git a/Project.toml b/Project.toml index b367288..8564c52 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ParameterHandling" uuid = "2412ca09-6db7-441c-8e3a-88d5709968c5" authors = ["Invenia Technical Computing Corporation"] -version = "0.4.10" +version = "0.5.0" [deps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" diff --git a/src/parameters_matrix.jl b/src/parameters_matrix.jl index 9b7e453..f6cab08 100644 --- a/src/parameters_matrix.jl +++ b/src/parameters_matrix.jl @@ -39,31 +39,36 @@ function flatten(::Type{T}, X::Orthogonal) where {T<:Real} end """ - positive_definite(X::AbstractMatrix{<:Real}) + positive_definite(X::AbstractMatrix{<:Real}, ε = eps(T)) -Produce a parameter whose `value` is constrained to be a positive-definite matrix. The argument `X` needs to -be a positive-definite matrix (see https://en.wikipedia.org/wiki/Definite_matrix). +Produce a parameter whose `value` is constrained to be a strictly positive-definite +matrix. The argument `X` minus `ε` times the identity needs to be a positive-definite matrix +(see https://en.wikipedia.org/wiki/Definite_matrix). The optional second argument `ε` must +be a positive real number. The unconstrained parameter is a `LowerTriangular` matrix, stored as a vector. """ -function positive_definite(X::AbstractMatrix{<:Real}) - isposdef(X) || throw(ArgumentError("X is not positive-definite")) - return PositiveDefinite(tril_to_vec(cholesky(X).L)) +function positive_definite(X::AbstractMatrix{T}, ε=eps(T)) where {T<:Real} + ε > 0 || throw(ArgumentError("ε must be positive.")) + _X = X - ε * I + isposdef(_X) || throw(ArgumentError("X-ε*I is not positive-definite for ε=$ε")) + return PositiveDefinite(tril_to_vec(cholesky(_X).L), ε) end -struct PositiveDefinite{TL<:AbstractVector{<:Real}} <: AbstractParameter +struct PositiveDefinite{TL<:AbstractVector{<:Real},Tε<:Real} <: AbstractParameter L::TL + ε::Tε end -Base.:(==)(X::PositiveDefinite, Y::PositiveDefinite) = X.L == Y.L - A_At(X) = X * X' -value(X::PositiveDefinite) = A_At(vec_to_tril(X.L)) +Base.:(==)(X::PositiveDefinite, Y::PositiveDefinite) = X.L == Y.L && X.ε == Y.ε + +value(X::PositiveDefinite) = A_At(vec_to_tril(X.L)) + X.ε * I function flatten(::Type{T}, X::PositiveDefinite) where {T<:Real} v, unflatten_v = flatten(T, X.L) - unflatten_PositiveDefinite(v_new::Vector{T}) = PositiveDefinite(unflatten_v(v_new)) + unflatten_PositiveDefinite(v_new::Vector{T}) = PositiveDefinite(unflatten_v(v_new), X.ε) return v, unflatten_PositiveDefinite end diff --git a/test/parameters_matrix.jl b/test/parameters_matrix.jl index cd12028..717f77c 100644 --- a/test/parameters_matrix.jl +++ b/test/parameters_matrix.jl @@ -32,11 +32,21 @@ using ParameterHandling: vec_to_tril, tril_to_vec end X_mat = ParameterHandling.A_At(rand(3, 3)) # Create a positive definite object X = positive_definite(X_mat) - @test X == X @test value(X) ≈ X_mat @test isposdef(value(X)) @test vec_to_tril(X.L) ≈ cholesky(X_mat).L - @test_throws ArgumentError positive_definite(rand(3, 3)) + + # Zeroing the unconstrained value preserve positivity + X.L .= 0 + @test isposdef(value(X)) + + # Check that zeroing the unconstrained value does not affect the original array + @test X != positive_definite(X_mat) + + @test positive_definite(X_mat, 1e-5) != positive_definite(X_mat, 1e-6) + @test_throws ArgumentError positive_definite(zeros(3, 3)) + @test_throws ArgumentError positive_definite(X_mat, 0.0) + test_parameter_interface(X) x, re = flatten(X) @@ -46,12 +56,6 @@ using ParameterHandling: vec_to_tril, tril_to_vec return logdet(value(X)) end, ) - ΔL = first( - Zygote.gradient(vec_to_tril(X.L)) do L - return logdet(L * L') - end, - ) - @test vec_to_tril(Δl) == tril(ΔL) ChainRulesTestUtils.test_rrule(vec_to_tril, x) end end