Introduce strategy for stable Cholesky factorizations

Project.toml

@@ -8,6 +8,7 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+PositiveFactorizations = "85a6dd25-e78a-55b7-8502-1745935b8125"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
@@ -18,9 +19,9 @@ StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
-julia = "1"
 Documenter = "0.25.2"
 ForwardDiff = "0.10.12"
 SpecialFunctions = "0.8.0, 0.10.3"
 StatsBase = "0.32.2, 0.33.1"
 StatsFuns = "0.9.5"
+julia = "1"

src/ForneyLab.jl

@@ -11,10 +11,11 @@ using Printf: @sprintf
 using StatsFuns: logmvgamma, betainvcdf, gammainvcdf, poisinvcdf
 using ForwardDiff
 using StatsBase: Weights
+using PositiveFactorizations
 
 import Statistics: mean, var, cov
 import Base: +, -, *, ^, ==, exp, convert, show, prod!
-import LinearAlgebra: dot
+import LinearAlgebra: dot, cholesky
 import StatsBase: sample
 
 # Helpers

src/factor_nodes/gaussian.jl

@@ -108,13 +108,13 @@ end
 function sample(dist::ProbabilityDistribution{Multivariate, F}) where F<:Gaussian
     isProper(dist) || error("Cannot sample from improper distribution")
     (m,V) = unsafeMeanCov(dist)
-    return (cholesky(V)).U' *randn(dims(dist)) + m
+    return (cholesky(default_cholesky_mode, V)).U' *randn(dims(dist)) + m
 end
 
 function sample(dist::ProbabilityDistribution{Multivariate, F}, n_samples::Int64) where F<:Gaussian
     isProper(dist) || error("Cannot sample from improper distribution")
     (m,V) = unsafeMeanCov(dist)
-    U = (cholesky(V)).U
+    U = (cholesky(default_cholesky_mode, V)).U
     d = dims(dist)
 
     return [U' *randn(d) + m for i in 1:n_samples]

src/factor_nodes/sample_list.jl

@@ -244,7 +244,7 @@ function bootstrap(dist_mean::ProbabilityDistribution{Multivariate, SampleList},
     s_m = dist_mean.params[:s] # Samples representing the mean
     N = length(s_m)
     V = dist_var.params[:m] # Fixed variance
-    U = (cholesky(V)).U # Precompute Cholesky
+    U = (cholesky(default_cholesky_mode, V)).U # Precompute Cholesky
 
     return [U' *randn(d) + s_m[i] for i in 1:N] # New samples
 end
@@ -263,7 +263,7 @@ function bootstrap(dist_mean::ProbabilityDistribution{Multivariate, <:Gaussian},
     s_V = dist_var.params[:s] # Samples representing the covariance
     N = length(s_V)
     (m, V) = unsafeMeanCov(dist_mean)
-    s_U = [(cholesky(s_V[i] + V)).U for i in 1:N] # Precompute Cholesky for each covariance sample; this can be expensive
+    s_U = [(cholesky(default_cholesky_mode, s_V[i] + V)).U for i in 1:N] # Precompute Cholesky for each covariance sample; this can be expensive
 
     return [s_U[i]' *randn(d) + m for i in 1:N] # New samples
 end

src/helpers.jl

@@ -1,16 +1,19 @@
 export huge, tiny, cholinv, diageye, eye, format, *, ^, mat, step!, init
 
+# Constants to define smallest/largest supported numbers.
+# Used for clipping quantities to ensure numerical stability.
+const huge = 1e12
+const tiny = 1e-12
+
+# Forces cholesky argument to positive definite; replace with "Nothing" for cholesky to error on non-pd argument
+const default_cholesky_mode = Positive
+
 """Cast input to a `Matrix` if necessary"""
 ensureMatrix(arr::AbstractMatrix{T}) where T<:Number = arr
 ensureMatrix(arr::Vector{T}) where T<:Number = Diagonal(arr)
 ensureMatrix(n::Number) = fill!(Array{typeof(n)}(undef,1,1), n)
 ensureMatrix(n::Nothing) = nothing
 
-# Constants to define smallest/largest supported numbers.
-# Used for clipping quantities to ensure numerical stability.
-const huge = 1e12
-const tiny = 1e-12
-
 """
 Wrapper for `logabsgamma` function that returns first element of its output
 """
@@ -25,29 +28,15 @@ function labsbeta(x::Number, y::Number)
     return logabsbeta(x, y)[1]
 end
 
+cholesky(::Type{Nothing}, M::AbstractMatrix) = cholesky(Hermitian(Matrix(M))) # No strategy for enforcing PD-ness of M
 
 """
-Matrix inversion using Cholesky decomposition,
-attempts with added regularization (1e-8*I) on failure.
+Matrix inversion using Cholesky decomposition
 """
-function cholinv(M::AbstractMatrix)
-    try
-        return inv(cholesky(Hermitian(Matrix(M))))
-    catch
-        try
-            return inv(cholesky(Hermitian(Matrix(M) + 1e-8*I)))
-        catch exception
-            if isa(exception, PosDefException)
-                error("PosDefException: Matrix is not positive-definite, even after regularization. $(typeof(M)):\n$M")
-            else
-                println("cholinv() errored when inverting $(typeof(M)):\n$M")
-                rethrow(exception)
-            end
-        end
-    end
-end
+cholinv(M::AbstractMatrix) = inv(cholesky(default_cholesky_mode, M))
 cholinv(m::Number) = 1.0/m
 cholinv(D::Diagonal) = Diagonal(1 ./ D.diag)
+
 eye(n::Number) = Diagonal(I,n)
 diageye(dims::Int64) = Diagonal(ones(dims))
 

test/test_helpers.jl

@@ -39,13 +39,9 @@ using LinearAlgebra: Diagonal, isposdef, I, Hermitian
         A = [2.0 1.0; 1.0 2.0]
         @test cholinv(A) ≈ inv(A)
         B = [2.0 1.0; 1.0 1.0] #PD
-        C = [1.0 1.0-1e-18; 1.0-1e-18 1.0] #non-PD due to numerical precision
-        D = [5.000000030387355e7 -4.9999999803873554e7;
-            -4.9999999803873554e7 5.0000000303873554e7] # computed inverse of C
         E = [1.0 1.0+1e8; 1.0+1e8 1.0] #non-PD not due to numerical precision
         @test isApproxEqual(cholinv(B), [1.0 -1.0; -1.0 2.0])
-        @test isApproxEqual(cholinv(C), D)
-        @test_throws Exception cholinv(E)
+        @test_throws Exception cholesky(Nothing, E)
         F = Diagonal([2.0, 3.0])
         @test cholinv(F) == inv(F)
         G = [Diagonal([2.0]) reshape([1.0],1,1); reshape([1.0],1,1) Diagonal([2.0])]