Add GammaLoss

src/GCPDecompositions.jl

@@ -21,6 +21,7 @@ export AbstractLoss,
     NonnegativeLeastSquaresLoss,
     PoissonLoss,
     PoissonLogLoss,
+    GammaLoss,
     UserDefinedLoss
 
 include("type-cpd.jl")

src/type-losses.jl

@@ -90,6 +90,28 @@ value(::PoissonLogLoss, x, m) = exp(m) - x * m
 deriv(::PoissonLogLoss, x, m) = exp(m) - x
 domain(::PoissonLogLoss) = Interval(-Inf, +Inf)
 
+"""
+    GammaLoss(eps::Real = 1e-10)
+
+Loss corresponding to a statistical assumption of Gamma-distributed data `X`
+with scale given by the low-rank model tensor `M`.
+
+- **Distribution:** ``x_i \\sim \\operatorname{Gamma}(k_i, \\theta_i)``
+- **Link function:** ``m_i = k_i \\sigma_i``
+- **Loss function:** ``f(x,m) = \\frac{x}{m + \\epsilon} + \\log(m + \\epsilon)``
+- **Domain:** ``m \\in [0, \\infty)``
+"""
+struct GammaLoss{T<:Real} <: AbstractLoss 
+  eps::T
+  GammaLoss{T}(eps::T) where {T<:Real} =
+    eps >= zero(eps) ? new(eps) :
+    throw(DomainError(eps, "Gamma loss requires nonnegative `eps`"))
+end
+GammaLoss(eps::T = 1e-10) where {T<:Real} = GammaLoss{T}(eps)
+value(loss::GammaLoss, x, m) = log(m + loss.eps) + x / (m + loss.eps)
+deriv(loss::GammaLoss, x, m) = -1 * (x / (m^2 + loss.eps)) + (1 / (m + loss.eps))
+domain(::GammaLoss) = Interval(0.0, +Inf)
+
 # User-defined loss
 
 """

test/items/gcp-opt.jl

@@ -91,6 +91,33 @@ end
     end
 end
 
+@testitem "GammaLoss" begin
+    using Random
+    using Distributions
+
+    @testset "size(X)=$sz, rank(X)=$r" for sz in [(15, 20, 25), (30, 40, 50)], r in 1:2
+        Random.seed!(0)
+        M = CPD(ones(r), randn.(sz, r))
+        X = [rand(Gamma(exp(M[I]))) for I in CartesianIndices(size(M))]
+
+        # Compute reference
+        Random.seed!(0)
+        Mr = GCPDecompositions._gcp(
+            X,
+            r,
+            (x, m) -> log(m + 1e-10) + x / (m + 1e-10),
+            (x, m) -> -1 * (x / (m^2 + 1e-10)) + (1 / (m + 1e-10)),
+            0.0,
+            (;),
+        )
+
+        # Test 
+        Random.seed!(0)
+        Mh = gcp(X, r, GammaLoss())
+        @test maximum(I -> abs(Mh[I] - Mr[I]), CartesianIndices(X)) <= 1e-5
+    end
+end
+
 @testitem "UserDefinedLoss" begin
     using Random, Distributions, IntervalSets