Add epsilon to PositiveDefinite

src/ParameterHandling.jl

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

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
 

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