Expensive jac_structure!

[amontoison]

Because of the linear API, the cost of the function jac_structure! is now prohibitive for large scale problems.
If we have as much linear constraints as nonlinear constraints, we need 2 * (ncon / 2)^2 comparisons (hidden in the calls to count).

NLPModels.jl/src/nlp/api.jl

Lines 192 to 221 in ac1213a

	function jac_structure!(
	nlp::AbstractNLPModel,
	rows::AbstractVector{T},
	cols::AbstractVector{T},
	) where {T}
	@lencheck nlp.meta.nnzj rows cols
	lin_ind = 1:(nlp.meta.lin_nnzj)
	if nlp.meta.nlin > 0
	if nlp.meta.nnln == 0
	jac_lin_structure!(nlp, rows, cols)
	else
	jac_lin_structure!(nlp, view(rows, lin_ind), view(cols, lin_ind))
	for i in lin_ind
	rows[i] += count(x < nlp.meta.lin[rows[i]] for x in nlp.meta.nln)
	end
	end
	end
	if nlp.meta.nnln > 0
	if nlp.meta.nlin == 0
	jac_nln_structure!(nlp, rows, cols)
	else
	nln_ind = (nlp.meta.lin_nnzj + 1):(nlp.meta.lin_nnzj + nlp.meta.nln_nnzj)
	jac_nln_structure!(nlp, view(rows, nln_ind), view(cols, nln_ind))
	for i in nln_ind
	rows[i] += count(x < nlp.meta.nln[rows[i]] for x in nlp.meta.lin)
	end
	end
	end
	return rows, cols
	end

[tmigot]

I think we can improve the default implementation but it will be marginal. Essentially, it's up to the implementation of the NLPModel API to offer an improved implementation (if possible). So, we should probably think of ways to really document the potential performance issues?

[amontoison]

So it means that we should always implement both API if possible like I did in NLPModelsJuMP.jl ?
cons!, cons_lin!, cons_nln!, jac_structure!, jac_structure_lin!, jac_structure_nln!, etc...

From my discussion with François, we have worst performance if we use the linear API than without.

[tmigot]

Yes, if there a possible optimization we should reimplement those as well.

By default, it's slightly worse but it helps better implementation as solvers will take advantage of splitting the constraints.

[amontoison]

To the best of my knowledge, KNITRO.jl is the only one that takes advantage of splitting the constraints.

[amontoison]

What could be relevant is to store a boolean in the NLPModelMeta to specify that the constraints are already in the good order (linear and then nonlinear). It will avoid the calls to count.
We have control on that in many models (like NLPModelsJuMP.jl).

[tmigot]

What could be the relevant is store a boolean in the NLPModelMeta to specify that the constraints are already in the good order (linear and then nonlinear). It will avoid the calls to count.

It's a good idea.

[dpo]

It may be simpler to adopt a fixed convention in NLPModels that linear constraints come before nonlinear constraints, or the other way around. For example,

• AMPL always orders nonlinear constraints first: see Table 5 (p.11) of Hooking your solver to AMPL (apparently, there's no option to change that)
• CUTEst has a flag to order nonlinear constraints first: set l_order = 2 in the call to csetup().

I haven't investigated JuMP. Currently, ADNLPModels always puts linear constraints first, but I suppose it wouldn't be hard to put them last.

It feels like such a convention wouldn't be very hard to implement and propagate. What do you think?

@amontoison @tmigot

[amontoison]

I like this idea! 👍👍

In NLPModelsJuMP.jl, the ordering is linear --> quadratic --> nonlinear but we have full freedom to change that. We extract the constraints one by one from the JuMP model.

We can also do what we want for ADNLPModels.jl.

I am also wondering if we can not keep one API if we do that (cons!, jac_coord!, ...) where an add an extra keyword argument to specify if we want all constraints or only a subset ("lp", "nlp", "full"). The default value is "full" and we can add a boolean in NLPModelMetato specify if the two other variants are also supported.

@blegat @cvanaret @frapac As external users of NLPModels.jl, what do you think of the suggestion of Dominique?

[blegat]

Having a convention on the ordering is a good idea. From the point of view of NLPModelsJuMP.jl, it will be easy to adapt to whatever the convention is. JuMP constraint indices only need to be different if they have the same constraint type precisely because solvers usually internally do not mix them anyway so they keep offsets and then the index of constraint within constraints of the same type.

[frapac]

I concur with Benoit, I htink having a convention on the ordering is a good idea.