Merge pull request #72 from JuliaGaussianProcesses/positive-semidefinite

Add epsilon to `PositiveDefinite`

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"

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