From 96f867961daf1bb8034d6b82fe9066fa6ccadb88 Mon Sep 17 00:00:00 2001 From: alexmul1114 <96952477+alexmul1114@users.noreply.github.com> Date: Wed, 11 Oct 2023 19:40:39 -0400 Subject: [PATCH 01/12] Add GammaLoss --- src/GCPDecompositions.jl | 1 + src/type-losses.jl | 22 ++++++++++++++++++++++ test/items/gcp-opt.jl | 27 +++++++++++++++++++++++++++ 3 files changed, 50 insertions(+) diff --git a/src/GCPDecompositions.jl b/src/GCPDecompositions.jl index afbe72f..c1583db 100644 --- a/src/GCPDecompositions.jl +++ b/src/GCPDecompositions.jl @@ -21,6 +21,7 @@ export AbstractLoss, NonnegativeLeastSquaresLoss, PoissonLoss, PoissonLogLoss, + GammaLoss, UserDefinedLoss include("type-cpd.jl") diff --git a/src/type-losses.jl b/src/type-losses.jl index a49e608..cd69362 100644 --- a/src/type-losses.jl +++ b/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 """ diff --git a/test/items/gcp-opt.jl b/test/items/gcp-opt.jl index bb353b3..139b552 100644 --- a/test/items/gcp-opt.jl +++ b/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 From 5ca6aa9297986ccc854ce6fc94363e3a6dcdd674 Mon Sep 17 00:00:00 2001 From: Alex Mulrooney <96952477+alexmul1114@users.noreply.github.com> Date: Thu, 12 Oct 2023 09:16:00 -0400 Subject: [PATCH 02/12] Update src/type-losses.jl Co-authored-by: David Hong --- src/type-losses.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/type-losses.jl b/src/type-losses.jl index cd69362..2bde3c3 100644 --- a/src/type-losses.jl +++ b/src/type-losses.jl @@ -96,8 +96,8 @@ domain(::PoissonLogLoss) = Interval(-Inf, +Inf) 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`` +- **Distribution:** ``x_i \\sim \\operatorname{Gamma}(k, \\sigma_i)`` +- **Link function:** ``m_i = k \\sigma_i`` - **Loss function:** ``f(x,m) = \\frac{x}{m + \\epsilon} + \\log(m + \\epsilon)`` - **Domain:** ``m \\in [0, \\infty)`` """ From 67a2a9ac3bc94b469466cce9eb087e4f85610c83 Mon Sep 17 00:00:00 2001 From: Alex Mulrooney <96952477+alexmul1114@users.noreply.github.com> Date: Thu, 12 Oct 2023 09:16:17 -0400 Subject: [PATCH 03/12] Update src/type-losses.jl Co-authored-by: David Hong --- src/type-losses.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/type-losses.jl b/src/type-losses.jl index 2bde3c3..ee9687c 100644 --- a/src/type-losses.jl +++ b/src/type-losses.jl @@ -108,7 +108,7 @@ struct GammaLoss{T<:Real} <: AbstractLoss 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) +value(loss::GammaLoss, x, m) = x / (m + loss.eps) + log(m + loss.eps) deriv(loss::GammaLoss, x, m) = -1 * (x / (m^2 + loss.eps)) + (1 / (m + loss.eps)) domain(::GammaLoss) = Interval(0.0, +Inf) From c324486cef66f2cd1af08799d0fc014593d9fe99 Mon Sep 17 00:00:00 2001 From: Alex Mulrooney <96952477+alexmul1114@users.noreply.github.com> Date: Thu, 12 Oct 2023 09:16:23 -0400 Subject: [PATCH 04/12] Update src/type-losses.jl Co-authored-by: David Hong --- src/type-losses.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/type-losses.jl b/src/type-losses.jl index ee9687c..575cc83 100644 --- a/src/type-losses.jl +++ b/src/type-losses.jl @@ -109,7 +109,7 @@ struct GammaLoss{T<:Real} <: AbstractLoss end GammaLoss(eps::T = 1e-10) where {T<:Real} = GammaLoss{T}(eps) value(loss::GammaLoss, x, m) = x / (m + loss.eps) + log(m + loss.eps) -deriv(loss::GammaLoss, x, m) = -1 * (x / (m^2 + loss.eps)) + (1 / (m + loss.eps)) +deriv(loss::GammaLoss, x, m) = -x / (m + loss.eps)^2 + inv(m + loss.eps) domain(::GammaLoss) = Interval(0.0, +Inf) # User-defined loss From 1b032e17221282526ee150a3f13f417877cab38d Mon Sep 17 00:00:00 2001 From: David Hong Date: Thu, 12 Oct 2023 15:39:54 -0400 Subject: [PATCH 05/12] Change test problem to better match stat assump --- test/items/gcp-opt.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/items/gcp-opt.jl b/test/items/gcp-opt.jl index 139b552..742ff44 100644 --- a/test/items/gcp-opt.jl +++ b/test/items/gcp-opt.jl @@ -97,8 +97,8 @@ end @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))] + M = CPD(ones(r), rand.(sz, r)) + X = [rand(Gamma(1.5, M[I])) for I in CartesianIndices(size(M))] # Compute reference Random.seed!(0) @@ -106,7 +106,7 @@ end X, r, (x, m) -> log(m + 1e-10) + x / (m + 1e-10), - (x, m) -> -1 * (x / (m^2 + 1e-10)) + (1 / (m + 1e-10)), + (x, m) -> -1 * (x / (m + 1e-10)^2) + (1 / (m + 1e-10)), 0.0, (;), ) From 2f59b37d77776fc2316ce4b86d117d213ceb57ab Mon Sep 17 00:00:00 2001 From: David Hong Date: Thu, 12 Oct 2023 15:44:33 -0400 Subject: [PATCH 06/12] Tweak test case to match a bit better --- test/items/gcp-opt.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/items/gcp-opt.jl b/test/items/gcp-opt.jl index 742ff44..3aa0e71 100644 --- a/test/items/gcp-opt.jl +++ b/test/items/gcp-opt.jl @@ -98,7 +98,8 @@ end @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), rand.(sz, r)) - X = [rand(Gamma(1.5, M[I])) for I in CartesianIndices(size(M))] + k = 1.5 + X = [rand(Gamma(k, M[I]/k)) for I in CartesianIndices(size(M))] # Compute reference Random.seed!(0) From f62fcbf105f00b043e19d7c0597352f3f4413079 Mon Sep 17 00:00:00 2001 From: alexmul1114 <96952477+alexmul1114@users.noreply.github.com> Date: Mon, 16 Oct 2023 09:47:25 -0400 Subject: [PATCH 07/12] Add Rayleigh Loss --- src/GCPDecompositions.jl | 1 + src/type-losses.jl | 22 ++++++++++++++++++++++ test/items/gcp-opt.jl | 36 ++++++++++++++++++++++++++++++++---- 3 files changed, 55 insertions(+), 4 deletions(-) diff --git a/src/GCPDecompositions.jl b/src/GCPDecompositions.jl index c1583db..f74cc6c 100644 --- a/src/GCPDecompositions.jl +++ b/src/GCPDecompositions.jl @@ -22,6 +22,7 @@ export AbstractLoss, PoissonLoss, PoissonLogLoss, GammaLoss, + RayleighLoss, UserDefinedLoss include("type-cpd.jl") diff --git a/src/type-losses.jl b/src/type-losses.jl index 575cc83..9f8f9de 100644 --- a/src/type-losses.jl +++ b/src/type-losses.jl @@ -112,6 +112,28 @@ value(loss::GammaLoss, x, m) = x / (m + loss.eps) + log(m + loss.eps) deriv(loss::GammaLoss, x, m) = -x / (m + loss.eps)^2 + inv(m + loss.eps) domain(::GammaLoss) = Interval(0.0, +Inf) +""" + RayleighLoss(eps::Real = 1e-10) + +Loss corresponding to the statistical assumption of Rayleigh data `X` +with sacle given by the low-rank model tensor `M` + + - **Distribution:** ``x_i \\sim \\operatorname{Rayleigh}(\\theta_i)`` + - **Link function:** ``m_i = \\sqrt{\\frac{\\pi}{2}\\theta_i}`` + - **Loss function:** ``f(x, m) = 2\\log(m + \\epsilon) + \\frac{\\pi}{4}(\\frac{x}{m + \\epsilon})^2`` + - **Domain:** ``m \\in [0, \\infty)`` +""" +struct RayleighLoss{T<:Real} <: AbstractLoss + eps::T + RayleighLoss{T}(eps::T) where {T<:Real} = + eps >= zero(eps) ? new(eps) : + throw(DomainError(eps, "Rayleigh loss requires nonnegative `eps`")) +end +RayleighLoss(eps::T = 1e-10) where {T<:Real} = RayleighLoss{T}(eps) +value(loss::RayleighLoss, x, m) = 2*log(m + loss.eps) + (pi / 4) * ((x/(m + loss.eps))^2) +deriv(loss::RayleighLoss, x, m) = 2/(m + loss.eps) - (pi / 2) * (x^2 / (m + loss.eps)^3) +domain(::RayleighLoss) = Interval(0.0, +Inf) + # User-defined loss """ diff --git a/test/items/gcp-opt.jl b/test/items/gcp-opt.jl index 139b552..d77cb63 100644 --- a/test/items/gcp-opt.jl +++ b/test/items/gcp-opt.jl @@ -97,16 +97,17 @@ end @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))] + M = CPD(ones(r), rand.(sz, r)) + k = 1.5 + X = [rand(Gamma(k, M[I]/k)) 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)), + (x, m) -> 2*log(m + 1e-10) + (pi / 4) * ((x/(m + 1e-10))^2), + (x, m) -> 2/(m + 1e-10) - (pi / 2) * (x^2 / (m + 1e-10)^3), 0.0, (;), ) @@ -118,6 +119,33 @@ end end end +@testitem "RayleighLoss" 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), rand.(sz, r)) + X = [rand(Rayleigh(M[I]/(sqrt(pi/2)))) for I in CartesianIndices(size(M))] + + # Compute reference + Random.seed!(0) + Mr = GCPDecompositions._gcp( + X, + r, + (x, m) -> 2*log(m + 1e-10) + (pi / 4) * ((x/(m + 1e-10))^2), + (x, m) -> 2/(m + 1e-10) - (pi / 2) * (x^2 / (m + 1e-10)^3), + 0.0, + (;), + ) + + # Test + Random.seed!(0) + Mh = gcp(X, r, RayleighLoss()) + @test maximum(I -> abs(Mh[I] - Mr[I]), CartesianIndices(X)) <= 1e-5 + end +end + @testitem "UserDefinedLoss" begin using Random, Distributions, IntervalSets From 3a45a89563e7ed5d5eca3c85004acca2e05e63b6 Mon Sep 17 00:00:00 2001 From: alexmul1114 <96952477+alexmul1114@users.noreply.github.com> Date: Mon, 16 Oct 2023 20:20:03 -0400 Subject: [PATCH 08/12] Attempt NegativeBinomial --- src/GCPDecompositions.jl | 1 + src/type-losses.jl | 59 ++++++++++++++++++++++++++++++++++++---- test/items/gcp-opt.jl | 31 +++++++++++++++++++++ 3 files changed, 85 insertions(+), 6 deletions(-) diff --git a/src/GCPDecompositions.jl b/src/GCPDecompositions.jl index 9664bc3..44b6280 100644 --- a/src/GCPDecompositions.jl +++ b/src/GCPDecompositions.jl @@ -25,6 +25,7 @@ export AbstractLoss, RayleighLoss, BernoulliOddsLoss, BernoulliLogitLoss, + NegativeBinomialOddsLoss, UserDefinedLoss export GCPConstraints diff --git a/src/type-losses.jl b/src/type-losses.jl index a1c217d..5a02193 100644 --- a/src/type-losses.jl +++ b/src/type-losses.jl @@ -138,10 +138,10 @@ domain(::RayleighLoss) = Interval(0.0, +Inf) BernoulliOddsLoss(eps::Real = 1e-10) Loss corresponding to the statistical assumption of Bernouli data `X` -with success rate given by the low-rank model tensor `M` +with odds-sucess rate given by the low-rank model tensor `M` - **Distribution:** ``x_i \\sim \\operatorname{Bernouli}(\\rho_i)`` - - **Link function:** ``m_i = \\rho / (1 - \\rho)`` + - **Link function:** ``m_i = \\frac{\\rho_i}{1 - \\rho_i}`` - **Loss function:** ``f(x, m) = \\log(m + 1) - x\\log(m + \\epsilon)`` - **Domain:** ``m \\in [0, \\infty)`` """ @@ -149,7 +149,7 @@ struct BernoulliOddsLoss{T<:Real} <: AbstractLoss eps::T BernoulliOddsLoss{T}(eps::T) where {T<:Real} = eps >= zero(eps) ? new(eps) : - throw(DomainError(eps, "BernoulliOddsLoss loss requires nonnegative `eps`")) + throw(DomainError(eps, "BernoulliOddsLoss requires nonnegative `eps`")) end BernoulliOddsLoss(eps::T = 1e-10) where {T<:Real} = BernoulliOddsLoss{T}(eps) value(loss::BernoulliOddsLoss, x, m) = log(m + 1) - x * log(m + loss.eps) @@ -161,10 +161,10 @@ domain(::BernoulliOddsLoss) = Interval(0.0, +Inf) BernoulliLogitLoss(eps::Real = 1e-10) Loss corresponding to the statistical assumption of Bernouli data `X` -with log-success rate given by the low-rank model tensor `M` +with log odds-success rate given by the low-rank model tensor `M` - **Distribution:** ``x_i \\sim \\operatorname{Bernouli}(\\rho_i)`` - - **Link function:** ``m_i = \\log(\\rho_i / (1 - \\rho_i))`` + - **Link function:** ``m_i = \\log(\\frac{\\rho_i}{1 - \\rho_i})`` - **Loss function:** ``f(x, m) = \\log(1 + e^m) - xm`` - **Domain:** ``m \\in \\mathbb{R}`` """ @@ -172,7 +172,7 @@ struct BernoulliLogitLoss{T<:Real} <: AbstractLoss eps::T BernoulliLogitLoss{T}(eps::T) where {T<:Real} = eps >= zero(eps) ? new(eps) : - throw(DomainError(eps, "BernoulliLogitsLoss loss requires nonnegative `eps`")) + throw(DomainError(eps, "BernoulliLogitsLoss requires nonnegative `eps`")) end BernoulliLogitLoss(eps::T = 1e-10) where {T<:Real} = BernoulliLogitLoss{T}(eps) value(::BernoulliLogitLoss, x, m) = log(1 + exp(m)) - x * m @@ -180,6 +180,53 @@ deriv(::BernoulliLogitLoss, x, m) = exp(m) / (1 + exp(m)) - x domain(::BernoulliLogitLoss) = Interval(-Inf, +Inf) +""" + NegativeBinomialOddsLoss(r::Integer, eps::Real = 1e-10) + +Loss corresponding to the statistical assumption of Negative Binomial +data `X` with log odds failure rate given by the low-rank model tensor `M` + + - **Distribution:** ``x_i \\sim \\operatorname{NegativeBinomial}(r, \\rho_i) `` + - **Link function:** ``m = \\frac{\\rho}{1 - \\rho}`` + - **Loss function:** ``f(x, m) = (r + x) \\log(1 + m) - x\\log(m + \\epsilon) `` + - **Domain:** ``m \\in [0, \\infty)`` +""" +struct NegativeBinomialOddsLoss{S<:Integer, T<:Real} <: AbstractLoss + r::S + eps::T + function NegativeBinomialOddsLoss{S, T}(r::S, eps::T) where {S<: Integer, T<:Real} + eps >= zero(eps) ? new(eps) : + throw(DomainError(eps, "NegativeBinomialOddsLoss requires nonnegative `eps`")) + r >= zero(r) ? new(r) : + throw(DomainError(r, "NegativeBinomialOddsLoss requires nonnegative `r`")) + end +end +NegativeBinomialOddsLoss(r::S, eps::T = 1e-10) where {S<:Integer, T<:Real} = NegativeBinomialOddsLoss{S, T}(r, eps) +value(loss::NegativeBinomialOddsLoss, x, m) = (loss.r + x) * log(1 + m) - x * log(m + loss.eps) +deriv(loss::NegativeBinomialOddsLoss, x, m) = (loss.r + x) / (1 + m) - x / (m + loss.eps) +domain(::NegativeBinomialOddsLoss) = Interval(0.0, +Inf) + + +""" + HuberLoss(Δ::Real) + + Huber Loss for given Δ + + - **Loss function:** ``f(x, m) = (x - m)^2 if \\abs(x - m)\\leq\\Delta, 2\\Delta\\abs(x - m) - \\Delta^2 otherwise`` + - **Domain:** ``m \\in \\mathbb{R}`` +""" +struct HuberLoss{T<:Real} <: AbstractLoss + Δ::T + HuberLoss{T}(Δ::T) where {T<:Real} = + Δ >= zero(Δ) ? new(Δ) : + throw(DomainError(Δ, "HuberLoss requires nonnegative `Δ`")) +end +Huber(Δ::T) where {T<:Real} = BernoulliLogitLoss{T}(Δ) +value(::HuberLoss, x, m) = abs(x - m) <= Δ ? (x - m)^2 : 2 * Δ * abs(x - m) - Δ^2 +deriv(::HuberLoss, x, m) = abs(x - m) <= Δ ? -2 * (x - m) : +domain(::HuberLoss) = Interval(-Inf, +Inf) + + # User-defined loss """ diff --git a/test/items/gcp-opt.jl b/test/items/gcp-opt.jl index a20175b..4e885cb 100644 --- a/test/items/gcp-opt.jl +++ b/test/items/gcp-opt.jl @@ -232,6 +232,37 @@ end end end +@testitem "NegativeBinomialOddsLoss" 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), rand.(sz, r)) + num_failures = 5 + X = [rand(NegativeBinomial(num_failures, M[I]/(M[I] + 1))) for I in CartesianIndices(size(M))] + + # Compute reference + Random.seed!(0) + Mr = GCPDecompositions._gcp( + X, + r, + (x, m) -> (num_failures + x) * log(1 + m) - x * log(m + 1e-10), + (x, m) -> (num_failures + x) / (1 + m) - x / (m + 1e-10), + 0.0, + (;), + ) + + # Test + Random.seed!(0) + Mh = gcp(X, r, NegativeBinomialOddsLoss(num_failures)) + @test maximum(I -> abs(Mh[I] - Mr[I]), CartesianIndices(X)) <= 1e-5 + end +end + + + + @testitem "UserDefinedLoss" begin using Random, Distributions, IntervalSets From 7f3d83d7744a70f0382cbcae4e7dc203e9368854 Mon Sep 17 00:00:00 2001 From: Alex Mulrooney <96952477+alexmul1114@users.noreply.github.com> Date: Tue, 17 Oct 2023 13:04:39 -0400 Subject: [PATCH 09/12] Modify NegativeBinomialOddsLoss Constructor Co-authored-by: David Hong --- src/type-losses.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/type-losses.jl b/src/type-losses.jl index 5a02193..6c7e0b4 100644 --- a/src/type-losses.jl +++ b/src/type-losses.jl @@ -195,10 +195,11 @@ struct NegativeBinomialOddsLoss{S<:Integer, T<:Real} <: AbstractLoss r::S eps::T function NegativeBinomialOddsLoss{S, T}(r::S, eps::T) where {S<: Integer, T<:Real} - eps >= zero(eps) ? new(eps) : + eps >= zero(eps) || throw(DomainError(eps, "NegativeBinomialOddsLoss requires nonnegative `eps`")) - r >= zero(r) ? new(r) : + r >= zero(r) || throw(DomainError(r, "NegativeBinomialOddsLoss requires nonnegative `r`")) + new(r, eps) end end NegativeBinomialOddsLoss(r::S, eps::T = 1e-10) where {S<:Integer, T<:Real} = NegativeBinomialOddsLoss{S, T}(r, eps) From e777fb4b759c3b45b26d1fa99b85e4afc3ba919f Mon Sep 17 00:00:00 2001 From: alexmul1114 <96952477+alexmul1114@users.noreply.github.com> Date: Tue, 17 Oct 2023 14:17:16 -0400 Subject: [PATCH 10/12] Add HuberLoss and BetaDivergenceLoss --- src/GCPDecompositions.jl | 2 ++ src/type-losses.jl | 47 ++++++++++++++++++++++--- test/items/gcp-opt.jl | 76 ++++++++++++++++++++++++++++++++++++++++ 3 files changed, 121 insertions(+), 4 deletions(-) diff --git a/src/GCPDecompositions.jl b/src/GCPDecompositions.jl index 44b6280..b7c42c6 100644 --- a/src/GCPDecompositions.jl +++ b/src/GCPDecompositions.jl @@ -26,6 +26,8 @@ export AbstractLoss, BernoulliOddsLoss, BernoulliLogitLoss, NegativeBinomialOddsLoss, + HuberLoss, + BetaDivergenceLoss, UserDefinedLoss export GCPConstraints diff --git a/src/type-losses.jl b/src/type-losses.jl index 6c7e0b4..b107f20 100644 --- a/src/type-losses.jl +++ b/src/type-losses.jl @@ -222,14 +222,53 @@ struct HuberLoss{T<:Real} <: AbstractLoss Δ >= zero(Δ) ? new(Δ) : throw(DomainError(Δ, "HuberLoss requires nonnegative `Δ`")) end -Huber(Δ::T) where {T<:Real} = BernoulliLogitLoss{T}(Δ) -value(::HuberLoss, x, m) = abs(x - m) <= Δ ? (x - m)^2 : 2 * Δ * abs(x - m) - Δ^2 -deriv(::HuberLoss, x, m) = abs(x - m) <= Δ ? -2 * (x - m) : +HuberLoss(Δ::T) where {T<:Real} = HuberLoss{T}(Δ) +value(loss::HuberLoss, x, m) = abs(x - m) <= loss.Δ ? (x - m)^2 : 2 *loss. Δ * abs(x - m) - loss.Δ^2 +deriv(loss::HuberLoss, x, m) = abs(x - m) <= loss.Δ ? -2 * (x - m) : -2 * sign(x - m) * loss.Δ * x domain(::HuberLoss) = Interval(-Inf, +Inf) -# User-defined loss +""" + BetaDivergenceLoss(β::Real, eps::Real) + + BetaDivergence Loss for given β + + - **Loss function:** ``f(x, m; β) = \\frac{1}{\\beta}m^{\\beta} - \\frac{1}{\\beta - 1}xm^{\\beta - 1} + if \\beta \\in \\mathbb{R} \\{0, 1\\}, + m - x\\log(m) if \\beta = 1, + \\frac{x}{m} + \\log(m) if \\beta = 0`` + - **Domain:** ``m \\in [0, \\infty)`` +""" +struct BetaDivergenceLoss{S<:Real, T<:Real} <: AbstractLoss + β::T + eps::T + BetaDivergenceLoss{S, T}(β::S, eps::T) where {S<:Real, T<:Real} = + eps >= zero(eps) ? new(β, eps) : + throw(DomainError(eps, "BetaDivergenceLoss requires nonnegative `eps`")) +end +BetaDivergenceLoss(β::S, eps::T = 1e-10) where {S<:Real, T<:Real} = BetaDivergenceLoss{S, T}(β, eps) +function value(loss::BetaDivergenceLoss, x, m) + if loss.β == 0 + return x / (m + loss.eps) + log(m + loss.eps) + elseif loss.β == 1 + return m - x * log(m + loss.eps) + else + return 1 / loss.β * m^loss.β - 1 / (loss.β - 1) * x * m^(loss.β - 1) + end +end +function deriv(loss::BetaDivergenceLoss, x, m) + if loss.β == 0 + return -x / (m + loss.eps)^2 + 1 / (m + loss.eps) + elseif loss.β == 1 + return 1 - x / (m + loss.eps) + else + return m^(loss.β - 1) - x * m^(loss.β - 2) + end +end +domain(::BetaDivergenceLoss) = Interval(0.0, +Inf) + +# User-defined loss """ UserDefinedLoss diff --git a/test/items/gcp-opt.jl b/test/items/gcp-opt.jl index 4e885cb..3cc350b 100644 --- a/test/items/gcp-opt.jl +++ b/test/items/gcp-opt.jl @@ -261,6 +261,82 @@ end end +@testitem "HuberLoss" 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), rand.(sz, r)) + X = [M[I] for I in CartesianIndices(size(M))] + + # Compute reference + Δ = 1 + Random.seed!(0) + Mr = GCPDecompositions._gcp( + X, + r, + (x, m) -> abs(x - m) <= Δ ? (x - m)^2 : 2 * Δ * abs(x - m) - Δ^2, + (x, m) -> abs(x - m) <= Δ ? -2 * (x - m) : -2 * sign(x - m) * Δ * x, + -Inf, + (;), + ) + + # Test + Random.seed!(0) + Mh = gcp(X, r, HuberLoss(Δ)) + @test maximum(I -> abs(Mh[I] - Mr[I]), CartesianIndices(X)) <= 1e-5 + end +end + + +@testitem "BetaDivergenceLoss" 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), rand.(sz, r)) + X = [rand(Poisson(M[I])) for I in CartesianIndices(size(M))] + + function beta_value(β, x, m) + if β == 0 + return x / (m + 1e-10) + log(m + 1e-10) + elseif β == 1 + return m - x * log(m + 1e-10) + else + return 1 / β * m^β - 1 / (β - 1) * x * m^(β - 1) + end + end + function beta_deriv(β, x, m) + if β == 0 + return -x / (m + 1e-10)^2 + 1 / (m + 1e-10) + elseif β == 1 + return 1 - x / (m + 1e-10) + else + return m^(β - 1) - x * m^(β - 2) + end + end + + # Compute reference + β = 1 + Random.seed!(0) + Mr = GCPDecompositions._gcp( + X, + r, + (x, m) -> beta_value(β, x, m), + (x, m) -> beta_deriv(β, x, m), + 0.0, + (;), + ) + + + # Test + Random.seed!(0) + Mh = gcp(X, r, BetaDivergenceLoss(β)) + @test maximum(I -> abs(Mh[I] - Mr[I]), CartesianIndices(X)) <= 1e-5 + end +end @testitem "UserDefinedLoss" begin From 6a97ffeb00b904f784d2e69ccf4b6bed3b5f98bb Mon Sep 17 00:00:00 2001 From: alexmul1114 <96952477+alexmul1114@users.noreply.github.com> Date: Tue, 17 Oct 2023 14:37:36 -0400 Subject: [PATCH 11/12] Add tests of multiple betas for BetaDivergenceLoss --- test/items/gcp-opt.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/test/items/gcp-opt.jl b/test/items/gcp-opt.jl index 3cc350b..dcce16d 100644 --- a/test/items/gcp-opt.jl +++ b/test/items/gcp-opt.jl @@ -294,7 +294,7 @@ end using Random using Distributions - @testset "size(X)=$sz, rank(X)=$r" for sz in [(15, 20, 25), (30, 40, 50)], r in 1:2 + @testset "size(X)=$sz, rank(X)=$r, β" for sz in [(15, 20, 25), (30, 40, 50)], r in 1:2, β in [0, 0.5, 1] Random.seed!(0) M = CPD(ones(r), rand.(sz, r)) X = [rand(Poisson(M[I])) for I in CartesianIndices(size(M))] @@ -319,7 +319,6 @@ end end # Compute reference - β = 1 Random.seed!(0) Mr = GCPDecompositions._gcp( X, @@ -330,7 +329,6 @@ end (;), ) - # Test Random.seed!(0) Mh = gcp(X, r, BetaDivergenceLoss(β)) From 25d728efb67580d4306a2bfa1cca262182029567 Mon Sep 17 00:00:00 2001 From: alexmul1114 <96952477+alexmul1114@users.noreply.github.com> Date: Tue, 17 Oct 2023 17:03:04 -0400 Subject: [PATCH 12/12] Add comment to BetaDivergenceLoss test --- test/items/gcp-opt.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/items/gcp-opt.jl b/test/items/gcp-opt.jl index dcce16d..7dba4e4 100644 --- a/test/items/gcp-opt.jl +++ b/test/items/gcp-opt.jl @@ -297,6 +297,7 @@ end @testset "size(X)=$sz, rank(X)=$r, β" for sz in [(15, 20, 25), (30, 40, 50)], r in 1:2, β in [0, 0.5, 1] Random.seed!(0) M = CPD(ones(r), rand.(sz, r)) + # May want to consider other distributions depending on value of β X = [rand(Poisson(M[I])) for I in CartesianIndices(size(M))] function beta_value(β, x, m)