Merge pull request #117 from tpapp/tp/fix-l2-remainder-singularity

Fix l2 remainder singularity, use tanh for calculations.
tpapp · Dec 13, 2023 · 1771aab · 1771aab · tpapp · Dec 13, 2023
2 parents 4b4d0a6 + 088c11c
commit 1771aab
Show file tree

Hide file tree

Showing 4 changed files with 31 additions and 7 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "TransformVariables"
 uuid = "84d833dd-6860-57f9-a1a7-6da5db126cff"
 authors = ["Tamas K. Papp <tkpapp@gmail.com>"]
-version = "0.8.8"
+version = "0.8.9"
 
 [deps]
 ArgCheck = "dce04be8-c92d-5529-be00-80e4d2c0e197"
@@ -11,7 +11,6 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
-Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SimpleUnPack = "ce78b400-467f-4804-87d8-8f486da07d0a"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
@@ -22,7 +21,9 @@ ChangesOfVariables = "0.1"
 DocStringExtensions = "0.8, 0.9"
 ForwardDiff = "0.10"
 InverseFunctions = "0.1"
+LinearAlgebra = "1.6"
 LogExpFunctions = "0.3"
+Random = "1.6"
 SimpleUnPack = "1"
 StaticArrays = "1"
 julia = "1.6"
diff --git a/src/special_arrays.jl b/src/special_arrays.jl
@@ -4,29 +4,39 @@ export UnitVector, UnitSimplex, CorrCholeskyFactor, corr_cholesky_factor
 #### building blocks
 ####
 
+"""
+$(SIGNATURES)
+
+`log(abs(…))` of the derivative of `tanh`, calculated accurately.
+"""
+function _tanh_logabsderiv(x)
+ d = 2*x
+ log(4) + d - 2 * log1pexp(d)
+end
+
 """
  (y, r, ℓ) = $SIGNATURES
 
 Given ``x ∈ ℝ`` and ``0 ≤ r ≤ 1``, return `(y, r′)` such that
 
-1. ``y² + r′² = r²``,
+1. ``y² + (r′)² = r²``,
 
 2. ``y: |y| ≤ r`` is mapped with a bijection from `x`.
 
 `ℓ` is the log Jacobian (whether it is evaluated depends on `flag`).
 """
 @inline function l2_remainder_transform(flag::LogJacFlag, x, r)
- z = 2*logistic(x) - 1
- (z * √r, r*(1 - abs2(z)),
- flag isa NoLogJac ? flag : log(2) + logistic_logjac(x) + 0.5*log(r))
+ # note that 1-tanh(x)^2 = sech(x)^2
+ (tanh(x) * √r, r*sech(x)^2,
+ flag isa NoLogJac ? flag : _tanh_logabsderiv(x) + 0.5*log(r))
 end
 
 """
  (x, r′) = $SIGNATURES
 
 Inverse of [`l2_remainder_transform`](@ref) in `x` and `y`.
 """
-@inline l2_remainder_inverse(y, r) = logit((y/√r+1)/2), r-abs2(y)
+@inline l2_remainder_inverse(y, r) = atanh(y/√r), r-y^2
 
 ####
 #### UnitVector

diff --git a/test/Project.toml b/test/Project.toml
@@ -8,6 +8,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LogDensityProblems = "6fdf6af0-433a-55f7-b3ed-c6c6e0b8df7c"
 LogDensityProblemsAD = "996a588d-648d-4e1f-a8f0-a84b347e47b1"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
+Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SimpleUnPack = "ce78b400-467f-4804-87d8-8f486da07d0a"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"

diff --git a/test/runtests.jl b/test/runtests.jl
@@ -659,3 +659,15 @@ end
  y.b[3] += 1
  @test x[4] == z + 1
 end
+
+@testset "near-singular Cholesky factor" begin
+ x = [8.348500225024523, -3.80486310849193, -15.115725300837742, 5.840812234057503,
+ 7.548980701857334, 0.6546495312434718, 2.863837638357627, 1.0081703617568052,
+ -38.543769810398466, -14.252165683848483, -22.75952203884357, 1.9543987098768612,
+ 5.415229912144962, -1.4360948924991273, 4.957606068283541, -5.443369115798325,
+ -2.536087079311158, -2.0710241403850635, -0.982209305513312, 6.821758239096414,
+ 5.925173901833287]
+ t = corr_cholesky_factor(7)
+ U = transform(t, x)
+ @test isfinite(logabsdet(U)[1])
+end