Add epsilon to PositiveDefinite

src/ParameterHandling.jl

@@ -8,7 +8,14 @@ using LinearAlgebra
 using SparseArrays
 
 export flatten,
-    value_flatten, positive, bounded, fixed, deferred, orthogonal, positive_definite
+    value_flatten,
+    positive,
+    bounded,
+    fixed,
+    deferred,
+    orthogonal,
+    positive_definite,
+    positive_semidefinite
 
 include("flatten.jl")
 include("parameters_base.jl")

src/parameters_matrix.jl

@@ -39,31 +39,72 @@ function flatten(::Type{T}, X::Orthogonal) where {T<:Real}
 end
 
 """
-    positive_definite(X::AbstractMatrix{<:Real})
+    positive_semidefinite(X::AbstractMatrix{<:Real})
 
-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 positive-semidefinite matrix. The
+argument `X` needs to be a positive-definite matrix
+(see https://en.wikipedia.org/wiki/Definite_matrix).
 
 The unconstrained parameter is a `LowerTriangular` matrix, stored as a vector.
+
+!!! warning
+    Even though the matrix needs to be positive-definite upon construction, the
+    unconstrained parameter can become zero, which represents a matrix which is merely
+    positive-semidefinite. To get a matrix that is always strictly positive-definite, use
+    `positive_definite`.
 """
-function positive_definite(X::AbstractMatrix{<:Real})
+function positive_semidefinite(X::AbstractMatrix{<:Real})
     isposdef(X) || throw(ArgumentError("X is not positive-definite"))
-    return PositiveDefinite(tril_to_vec(cholesky(X).L))
+    return PositiveSemiDefinite(tril_to_vec(cholesky(X).L))
+end
+
+"""
+    positive_definite(X::AbstractMatrix{<:Real}, ε = eps(T))
+
+Produce a parameter whose `value` is constrained to be a strictly positive-semidefinite
+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{T}, ε=eps(T)) where {T<:Real}
+    ε > 0 || throw(ArgumentError("ε is not positive. Use `positive_semidefinite` instead."))
+    _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 PositiveSemiDefinite{TL<:AbstractVector{<:Real}} <: AbstractParameter
     L::TL
 end
 
-Base.:(==)(X::PositiveDefinite, Y::PositiveDefinite) = X.L == Y.L
+Base.:(==)(X::PositiveSemiDefinite, Y::PositiveSemiDefinite) = X.L == Y.L
 
 A_At(X) = X * X'
 
-value(X::PositiveDefinite) = A_At(vec_to_tril(X.L))
+value(X::PositiveSemiDefinite) = A_At(vec_to_tril(X.L))
+
+function flatten(::Type{T}, X::PositiveSemiDefinite) where {T<:Real}
+    v, unflatten_v = flatten(T, X.L)
+    function unflatten_PositiveSemiDefinite(v_new::Vector{T})
+        return PositiveSemiDefinite(unflatten_v(v_new))
+    end
+    return v, unflatten_PositiveSemiDefinite
+end
+
+struct PositiveDefinite{TL<:AbstractVector{<:Real},Tε<:Real} <: AbstractParameter
+    L::TL
+    ε::Tε
+end
+
+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
 

test/parameters_matrix.jl

@@ -24,19 +24,19 @@ using ParameterHandling: vec_to_tril, tril_to_vec
         end
     end
 
-    @testset "positive_definite" begin
+    @testset "positive_semidefinite" begin
         @testset "vec_tril_conversion" begin
             X = tril!(rand(3, 3))
             @test vec_to_tril(tril_to_vec(X)) == X
             @test_throws ErrorException tril_to_vec(rand(4, 5))
         end
         X_mat = ParameterHandling.A_At(rand(3, 3)) # Create a positive definite object
-        X = positive_definite(X_mat)
+        X = positive_semidefinite(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))
+        @test_throws ArgumentError positive_semidefinite(rand(3, 3))
         test_parameter_interface(X)
 
         x, re = flatten(X)
@@ -54,4 +54,26 @@ using ParameterHandling: vec_to_tril, tril_to_vec
         @test vec_to_tril(Δl) == tril(ΔL)
         ChainRulesTestUtils.test_rrule(vec_to_tril, x)
     end
+
+    @testset "positive_definite" begin
+        X_mat = ParameterHandling.A_At(rand(3, 3)) # Create a positive definite object
+        X = positive_definite(X_mat)
+        @test isposdef(value(X))
+        X.L .= 0 # zero the unconstrained value
+        @test isposdef(value(X))
+        @test X != positive_definite(X_mat)
+        @test positive_definite(X_mat, 1e-3) != positive_definite(X_mat, 1e-2)
+        @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)
+        Δl = first(
+            Zygote.gradient(x) do x
+                X = re(x)
+                return logdet(value(X))
+            end,
+        )
+        ChainRulesTestUtils.test_rrule(vec_to_tril, x)
+    end
 end