From c706f6fbcc1ca2c94f94d604c41708ddba36d1e9 Mon Sep 17 00:00:00 2001 From: David Hong Date: Thu, 12 Oct 2023 16:13:35 -0400 Subject: [PATCH 1/4] Handle LossFunctions.jl via value, deriv, domain Reduces duplication and is probably a more natural way to implement this package extension. --- ext/LossFunctionsExt.jl | 17 +++++------------ src/gcp-opt.jl | 2 +- 2 files changed, 6 insertions(+), 13 deletions(-) diff --git a/ext/LossFunctionsExt.jl b/ext/LossFunctionsExt.jl index 4fd2ff4..46778f1 100644 --- a/ext/LossFunctionsExt.jl +++ b/ext/LossFunctionsExt.jl @@ -1,20 +1,13 @@ module LossFunctionsExt using GCPDecompositions, LossFunctions -import GCPDecompositions: _factor_matrix_lower_bound +using IntervalSets const SupportedLosses = Union{LossFunctions.DistanceLoss,LossFunctions.MarginLoss} -GCPDecompositions.gcp(X::Array, r, loss::SupportedLosses) = GCPDecompositions._gcp( - X, - r, - (x, m) -> loss(m, x), - (x, m) -> LossFunctions.deriv(loss, m, x), - _factor_matrix_lower_bound(loss), - (;), -) - -_factor_matrix_lower_bound(::LossFunctions.DistanceLoss) = -Inf -_factor_matrix_lower_bound(::LossFunctions.MarginLoss) = -Inf +GCPDecompositions.value(loss::SupportedLosses, x, m) = loss(m, x) +GCPDecompositions.deriv(loss::SupportedLosses, x, m) = LossFunctions.deriv(loss, m, x) +GCPDecompositions.domain(::LossFunctions.DistanceLoss) = Interval(-Inf, Inf) +GCPDecompositions.domain(::LossFunctions.MarginLoss) = Interval(-Inf, Inf) end diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl index 1b8a06a..95a03ab 100644 --- a/src/gcp-opt.jl +++ b/src/gcp-opt.jl @@ -15,7 +15,7 @@ to see what losses are supported. See also: `CPD`, `AbstractLoss`. """ -gcp(X::Array, r, loss::AbstractLoss = LeastSquaresLoss()) = _gcp( +gcp(X::Array, r, loss = LeastSquaresLoss()) = _gcp( X, r, (x, m) -> value(loss, x, m), From 1c19c6adf28548eb157bfa6d174813f6fe0460f5 Mon Sep 17 00:00:00 2001 From: David Hong Date: Thu, 12 Oct 2023 16:35:04 -0400 Subject: [PATCH 2/4] New main signature for _gcp Provide richer info to _gcp so we can dispatch at that level. Preparing the way for having various algorithms. --- src/gcp-opt.jl | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl index 95a03ab..d3834fc 100644 --- a/src/gcp-opt.jl +++ b/src/gcp-opt.jl @@ -15,14 +15,7 @@ to see what losses are supported. See also: `CPD`, `AbstractLoss`. """ -gcp(X::Array, r, loss = LeastSquaresLoss()) = _gcp( - X, - r, - (x, m) -> value(loss, x, m), - (x, m) -> deriv(loss, x, m), - _factor_matrix_lower_bound(loss), - (;), -) +gcp(X::Array, r, loss = LeastSquaresLoss()) = _gcp(X, r, loss, (;)) # Choose lower bound on factor matrix entries based on the domain of the loss function _factor_matrix_lower_bound(loss) @@ -52,6 +45,14 @@ function _factor_matrix_lower_bound(loss) return min end +_gcp(X::Array{TX,N}, r, loss, lbfgsopts) where {TX,N} = _gcp( + X, + r, + (x, m) -> value(loss, x, m), + (x, m) -> deriv(loss, x, m), + _factor_matrix_lower_bound(loss), + lbfgsopts, +) function _gcp(X::Array{TX,N}, r, func, grad, lower, lbfgsopts) where {TX,N} # T = promote_type(nonmissingtype(TX), Float64) T = Float64 # LBFGSB.jl seems to only support Float64 From 9647b7523c3f29400edd122832c0d6ab1f7b5164 Mon Sep 17 00:00:00 2001 From: David Hong Date: Thu, 12 Oct 2023 16:40:03 -0400 Subject: [PATCH 3/4] Add todo --- src/gcp-opt.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl index d3834fc..72f5639 100644 --- a/src/gcp-opt.jl +++ b/src/gcp-opt.jl @@ -45,6 +45,8 @@ function _factor_matrix_lower_bound(loss) return min end +# TODO: remove the older `func, grad, lower` signature +# will require reworking how we do testing _gcp(X::Array{TX,N}, r, loss, lbfgsopts) where {TX,N} = _gcp( X, r, From 4b5431a4248ef9b50dbe20b9fafb0f16350246cc Mon Sep 17 00:00:00 2001 From: David Hong Date: Thu, 12 Oct 2023 19:19:08 -0400 Subject: [PATCH 4/4] Make newer `_gcp` signature the main one Implement fall-back using `UserDefinedLoss`. --- src/gcp-opt.jl | 25 +++++++++++++------------ 1 file changed, 13 insertions(+), 12 deletions(-) diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl index 72f5639..5114ddf 100644 --- a/src/gcp-opt.jl +++ b/src/gcp-opt.jl @@ -45,20 +45,21 @@ function _factor_matrix_lower_bound(loss) return min end -# TODO: remove the older `func, grad, lower` signature +# TODO: remove this `func, grad, lower` signature # will require reworking how we do testing -_gcp(X::Array{TX,N}, r, loss, lbfgsopts) where {TX,N} = _gcp( +_gcp(X::Array{TX,N}, r, func, grad, lower, lbfgsopts) where {TX,N} = _gcp( X, r, - (x, m) -> value(loss, x, m), - (x, m) -> deriv(loss, x, m), - _factor_matrix_lower_bound(loss), + UserDefinedLoss(func; deriv = grad, domain = Interval(lower, +Inf)), lbfgsopts, ) -function _gcp(X::Array{TX,N}, r, func, grad, lower, lbfgsopts) where {TX,N} +function _gcp(X::Array{TX,N}, r, loss, lbfgsopts) where {TX,N} # T = promote_type(nonmissingtype(TX), Float64) T = Float64 # LBFGSB.jl seems to only support Float64 + # Choose lower bound on factor matrix entries based on the domain of the loss + lower = _factor_matrix_lower_bound(loss) + # Random initialization M0 = CPD(ones(T, r), rand.(T, size(X), r)) M0norm = sqrt(sum(abs2, M0[I] for I in CartesianIndices(size(M0)))) @@ -73,12 +74,12 @@ function _gcp(X::Array{TX,N}, r, func, grad, lower, lbfgsopts) where {TX,N} vec_ranges = ntuple(k -> vec_cutoffs[k]+1:vec_cutoffs[k+1], Val(N)) function f(u) U = map(range -> reshape(view(u, range), :, r), vec_ranges) - return gcp_func(CPD(ones(T, r), U), X, func) + return gcp_func(CPD(ones(T, r), U), X, loss) end function g!(gu, u) U = map(range -> reshape(view(u, range), :, r), vec_ranges) GU = map(range -> reshape(view(gu, range), :, r), vec_ranges) - gcp_grad_U!(GU, CPD(ones(T, r), U), X, grad) + gcp_grad_U!(GU, CPD(ones(T, r), U), X, loss) return gu end @@ -90,18 +91,18 @@ function _gcp(X::Array{TX,N}, r, func, grad, lower, lbfgsopts) where {TX,N} end # Objective function and gradient (w.r.t. `M.U`) -function gcp_func(M::CPD{T,N}, X::Array{TX,N}, func) where {T,TX,N} - return sum(func(X[I], M[I]) for I in CartesianIndices(X) if !ismissing(X[I])) +function gcp_func(M::CPD{T,N}, X::Array{TX,N}, loss) where {T,TX,N} + return sum(value(loss, X[I], M[I]) for I in CartesianIndices(X) if !ismissing(X[I])) end function gcp_grad_U!( GU::NTuple{N,TGU}, M::CPD{T,N}, X::Array{TX,N}, - grad, + loss, ) where {T,TX,N,TGU<:AbstractMatrix{T}} Y = [ - ismissing(X[I]) ? zero(nonmissingtype(eltype(X))) : grad(X[I], M[I]) for + ismissing(X[I]) ? zero(nonmissingtype(eltype(X))) : deriv(loss, X[I], M[I]) for I in CartesianIndices(X) ]