Added binomial model

[saravkin]

I think a basic binomial classifier is good to have. it is not a nonlinear least squares, but provides function/gradient/hessian.

[codecov]

Codecov Report

❌ Patch coverage is 0% with 29 lines in your changes missing coverage. Please review.
✅ Project coverage is 72.42%. Comparing base (85de184) to head (6490976).
⚠️ Report is 15 commits behind head on main.

Files with missing lines	Patch %	Lines
src/binomial_model.jl	0.00%	29 Missing ⚠️

Additional details and impacted files

@@             Coverage Diff             @@
##             main      #79       +/-   ##
===========================================
- Coverage   99.18%   72.42%   -26.76%     
===========================================
  Files          11       20        +9     
  Lines         369      486      +117     
===========================================
- Hits          366      352       -14     
- Misses          3      134      +131     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[dpo]

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).

[dpo]

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.

[dpo]

@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.

[saravkin]

Just updated the binomial model and added a test that shows it works as expected.

[dpo]

Many thanks! I made a number of suggestions to make your code more generic (it works more generally than with Float64; e.g., it should work in Float32 or Float128).

Also, there's no need to commit the Manifest.toml.

[dpo]

Suggested change

	Minimize
	f(x) = sum( log(1 + exp(a_i' * x)) - b_i * (a_i' * x) )
	where b_i ∈ {0, 1} and `A` is (m features) × (n samples).
	Minimize
	f(x) = sum( log(1 + exp(a_i' * x)) - b_i * (a_i' * x) )
	where b_i ∈ {0, 1} and `A` is (m features) × (n samples).

[dpo]

Not sure I understand this last sentence. There's no z in the description above.

[dpo]

Suggested change

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

[dpo]

Suggested change

	eltype(A) <: Real || throw(ArgumentError("A must have a real element type"))
	eltype(b) <: Real || throw(ArgumentError("b must have a real element type"))

[dpo]

Suggested change

	Ax = zeros(Float64, n) # z = A' * x
	p = zeros(Float64, n) # sigmoid(z)
	w = zeros(Float64, n) # p*(1-p)
	tmp_n = zeros(Float64, n) # sample-space buffer
	tmp_v = zeros(Float64, n) # sample-space buffer
	Ax = similar(b) # z = A' * x
	p = similar(b) # sigmoid(z)
	w = similar(b) # p*(1-p)
	tmp_n = similar(b) # sample-space buffer
	tmp_v = similar(b) # sample-space buffer

[dpo]

Suggested change

	s = 0.0
	@inbounds @simd for i in 1:n
	zi = Ax[i]
	s += softplus(zi) - Float64(b[i]) * zi
	s = zero(T)
	@inbounds @simd for i in 1:n
	zi = Ax[i]
	s += softplus(zi) - b[i] * zi

[dpo]

Suggested change

	tmp_n[i] = p[i] - Float64(b[i])
	tmp_n[i] = p[i] - b[i]

[dpo]

Suggested change

	w[i] = pi * (1.0 - pi)
	w[i] = pi * (1 - pi)

[dpo]

Suggested change

	if obj_weight != 1.0
	if obj_weight != 1

[dpo]

Suggested change

	x0 = zeros(Float64, m)
	x0 = fill!(similar(b, m), zero(T))