Add VectorNonlinearFunction

[odow]

Resurrecting VectorNonlinearFunction from the depths of #2059.

There are two main motivating problems for this:

1. VectorNonlinearFunction-in-Complements constrains for nonlinear complementarity problems
   • Add support for VectorNonlinearFunction chkwon/PATHSolver.jl#82
2. ObjectiveFunction{VectorNonlinearFunction} for nonlinear multi-objective problems
   • Add support for VectorNonlinearFunction objectives MultiObjectiveAlgorithms.jl#68

Documentation

https://jump.dev/MathOptInterface.jl/previews/PR2201/reference/standard_form/#MathOptInterface.VectorNonlinearFunction

I also removed this from #2059 because there we a lot of issues with the bridges, but I have a better handle on what to do now, and we can easily run solver-tests to catch the culprits.

https://github.com/jump-dev/MathOptInterface.jl/actions/runs/5218217869

[odow]

https://github.com/jump-dev/MathOptInterface.jl/actions/runs/5238529424

[odow]

Marking this as blocked by #2203

[blegat]

One thing to consider:

struct VectorFunction{F} <: MOI.AbstractVectorFunction
    elements::Vector{F}
end
const VectorOfVariables = VectorFunction{VariableIndex}
const VectorNonlinearFunction = VectorFunction{Any}

It could simplify because we then don't have to implement operate!, etc... for both. It might be a bit more complicated since it's more abstract but I would expect it to go smoothly since we have the promote_operation and operate API to work abstractly on F.
Another advantage is that we have another scalar function defined in https://github.com/jump-dev/PolyJuMP.jl/blob/master/src/functions.jl.
And that would allow us to do

const VectorPolynomialFunction = Vector{ScalarPolynomialFunction}

and get a working vector polynomial function in one line !

One might indeed wonder why we need to define our own vector type when Julia already has Base.Vector.
It makes sense for VectorAffineFunction since it is implementing it as a matrix and not a vector of vectors but for the other one, it's good to have a default vector function.

One could even define an AbstractFunctionConversionBridge from VectorFunction{ScalarAffineFunction} to VectorAffineFunction even if we probably still want to create a VectorAffineFunction from JuMP directly, even if the JuMP functions are really VectorFunction{ScalarAffineFunction} !

[odow]

This is an interesting suggestion. We'd need backward compat for (::VectorOfVariables).variables, but a generic AbstractVectorFunction that is a vector of scalar terms would be useful.

[odow]

This is not fun. I think I'm going to revert the GenericVectorFunction changes to MOI.VectorOfVariables. There are too many little issues that crop up. Perhaps we can revisit once this PR is merged, and once we have a solution for the operate mess.

[blegat]

There are too many little issues that crop up.

The tests seem to pass except for some printing issues though

[odow]

The tests seem to pass except for some printing issues though

The tests don't cover the code enough. If you try running the tests with Clp/SCS etc lots of things break.

[odow]

This is still blocked by the operate issue, but we're getting closer: #2215

[odow]

https://github.com/jump-dev/MathOptInterface.jl/actions/runs/5341127312

[odow]

If you try running the tests with SCS etc lots of things break.

Okay. One massive source of confusion for me had been failures with SCS. Which were because the bridge in SCS was wrong: jump-dev/SCS.jl#277.

It used to do this:

julia> MOI.Bridges.print_active_bridges(model, MOI.VectorNonlinearFunction, MOI.PositiveSemidefiniteConeTriangle)
 * Unsupported constraint: MOI.GenericVectorFunction{Any}-in-MOI.PositiveSemidefiniteConeTriangle
 |  bridged by:
 |   SCS.ScaledPSDConeBridge{Float64, G} where G
 |  may introduce:
 |   * Supported constraint: MOI.VectorAffineFunction{Float64}-in-SCS.ScaledPSDCone

So no wonder lots of weird unrelated things broke.

But I only dug deep enough once I figured out it wasn't a method error with the various operate methods.

[odow]

So I was looking into VectorNonlinearFunction-in-PositiveSemidefiniteConeTriangle etc. And there are a number of thorny issues:

• Utilities.eachscalar assumes that each row of the vector function is scalar_type(F). But the row of a VectorNonlinearFunction might not be ScalarNonlinearFunction.
• Easy fix is to convert every row, but then what's the difference between Utilities.scalarize and Utilities.eachscalar?
• Operations don't simplify. So f = sin(x) and then -(f, f) doesn't simplify to 0. So we get an expression swelling problem that prevents Utilities.isapprox_zero from working. The main problem with this is in the Bridges.Constraint.SquareBridge, which checks if off-diagonals are equivalent by checking if their difference is zero.

[blegat]

Easy fix is to convert every row, but then what's the difference between Utilities.scalarize and Utilities.eachscalar?

Yes, I think every row should be an SNF. The difference of scalarize is that it allocates the vector of functions while eachscalar does not. If you have a lot or rows but you're only going to get one, eachscalar will be faster.
In fact, scalarize is really like conversion to a GenericVectorFunction

The main problem with this is in the Bridges.Constraint.SquareBridge, which checks if off-diagonals are equivalent by checking if their difference is zero.

Maybe it should use isequal_canonical or something like that

[odow]

Should scalarize be a copy? Currently it isn't always:

MathOptInterface.jl/src/Utilities/functions.jl

Lines 2110 to 2120 in 5afc08b

	"""
	scalarize(func::MOI.VectorOfVariables, ignore_constants::Bool = false)
	Returns a vector of scalar functions making up the vector function in the form
	of a `Vector{MOI.SingleVariable}`.
	See also [`eachscalar`](@ref).
	"""
	function scalarize(f::MOI.VectorOfVariables, ignore_constants::Bool = false)
	return f.variables
	end

Ditto with vectorize

MathOptInterface.jl/src/Utilities/functions.jl

Lines 2044 to 2057 in 5afc08b

	"""
	vectorize(x::AbstractVector{<:Number})
	Returns `x`.
	"""
	vectorize(x::AbstractVector{<:Number}) = x
	"""
	vectorize(x::AbstractVector{MOI.VariableIndex})
	Returns the vector of scalar affine functions in the form of a
	`MOI.VectorAffineFunction{T}`.
	"""
	vectorize(x::AbstractVector{MOI.VariableIndex}) = MOI.VectorOfVariables(x)

[blegat]

I don't think it has to be a copy. MA.operate(scalarize, f) should be a copy 😆

[odow]

https://github.com/jump-dev/MathOptInterface.jl/actions/runs/5383997349

[odow]

PATHSolver.jl and MultiObjectiveAlgorithms.jl are both passing tests with this branch, so I don't think we're far away. Just need https://github.com/jump-dev/MathOptInterface.jl/actions/runs/5383997349 to finish.

[mlubin]

I'll take a look over the weekend