Skip to content

Added binomial model #79

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
saravkin wants to merge 6 commits into
base: main
Choose a base branch
from
Open

Added binomial model #79

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
f08b18b
added binomial model
saravkin Dec 28, 2025
6490976
added binomial model
saravkin Dec 28, 2025
7768d00
updated binomial model to have hprod
saravkin Dec 31, 2025
333934f
fixed binomial model
saravkin Jan 9, 2026
5bce344
Apply suggestion from @dpo
saravkin Feb 12, 2026
e0cba0c
Apply review suggestions: make binomial_model generic over element type
saravkin Feb 12, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions src/RegularizedProblems.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@ include("lrcomp_model.jl")
include("matrand_model.jl")
include("group_lasso_model.jl")
include("nnmf.jl")
include("binomial_model.jl")

function __init__()
@require ProximalOperators = "a725b495-10eb-56fe-b38b-717eba820537" begin
Expand Down
62 changes: 62 additions & 0 deletions src/binomial_model.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,62 @@
export binomial_model

"""
nlp_model, nls_model = binomial_model(A, b)

Return an instance of an `NLPModel` representing the binomial logistic regression
problem, and an `NLSModel` representing the system of equations formed by its gradient.
Copy link
Copy Markdown
Member

@dpo dpo Dec 30, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I understand now, but it’s a bit confusing because the two models don’t represent the same problem. I suggest returning the NLPModel only. For that, you only need implement obj and grad!. If I understand well, the current jacv! and jactv! are actually hprod! (the product between the Hessian and a vector), hence why they are the same.

It would be easy to write a generic wrapper to minimize $\tfrac{1}{2} \Vert \nabla f(x)\Vert^2$ for any given NLPModel.

Copy link
Copy Markdown
Member

@dpo dpo Jan 4, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@saravkin I think what you want to do is the following:

function binomial_model(A, b)
  # preallocate storage, etc...

  function obj(x)
    # return sum(exp(...))
  end

  function grad!(g, x)
    # ... what you put in `resid!()`
    # fill `g` with the gradient at `x`
  end

  function hprod!(hv, x, v; obj_weight = 1)
    # ... what you put in `jacv!()`
    # fill `hv` with the product of the Hessian at `x` with `v`
  end

  return NLPModel(x0, obj, grad = grad!, hprod = hprod!, meta_args = (name = "Binomial"))
end

Your problems isn't a least-squares problem, so you only return the NLPModel.


Minimize f(x) = sum(log(1+exp(a_i^T x)) - y_i a_i^T x)
where y_i in {0, 1}.

The NLS residual is defined as the gradient of f(x):
F(x) = A (p - b)
where p = sigmoid(A' x).
"""
function binomial_model(A, b)
Copy link
Copy Markdown
Member

@dpo dpo Jan 12, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
function binomial_model(A, b)
function binomial_model(A::AbstractMatrix{T}, b::AbstractVector{T}) where T <: Real

m, n = size(A) # m features, n samples

# Pre-allocate buffers
Ax = zeros(n)
p = zeros(n)
w = zeros(n)
tmp_n = zeros(n)

function resid!(r, x)
mul!(Ax, A', x)
@. p = 1.0 / (1.0 + exp(-Ax))
@. tmp_n = p - b
mul!(r, A, tmp_n)
return r
end

function jacv!(Jv, x, v)
mul!(Ax, A', x)
@. w = 1.0 / (1.0 + exp(-Ax))
@. w = w * (1.0 - w)

mul!(tmp_n, A', v)
@. tmp_n *= w
mul!(Jv, A, tmp_n)
return Jv
end

function jactv!(Jtv, x, v)
jacv!(Jtv, x, v)
end

function obj(x)
mul!(Ax, A', x)
# Stable computation of log(1+exp(z)) - b*z
return sum(@. (Ax > 0) * ((1 - b) * Ax + log(1 + exp(-Ax))) + (Ax <= 0) * (log(1 + exp(Ax)) - b * Ax))
Copy link
Copy Markdown
Member

@dpo dpo Dec 29, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I’m a bit confused. If this is a nonlinear least-squares problem, the objective should be $\tfrac{1}{2} \Vert r \Vert^2$, where $r$ results from resid!(r, x), and the gradient should be the vector g that results from jactv!(g, x, r).

end

function grad!(g, x)
resid!(g, x)
end

x0 = zeros(m)

return FirstOrderModel(obj, grad!, x0, name="Binomial"),
FirstOrderNLSModel(resid!, jacv!, jactv!, m, x0)
end
Loading