MathOptInterface.jl icon indicating copy to clipboard operation
MathOptInterface.jl copied to clipboard

Abstract / alternative linear operator storage

Open matbesancon opened this issue 5 years ago • 11 comments

The current way of representing A x is with a VectorAffineFunction, which has the following structure:

MOI.VectorAffineFunction(
    [MOI.VectorAffineTerm(
        row_idx,
        ScalarAffineTerm(
            coeff,
            col_idx,
        )
    )],
    [constants]
)

This has a lot of nested structure and imposes a fixed concrete layout:

  • one cannot take advantage of specific matrix structures
  • data gets moved around twice, from the matrix to a VectorAffineTerm, then from VAT to the solver.

Could we have an abstract linear operator function, VectorLinearOperator, with possibly a bridge to VectorAffineTerm to avoid excessive burden on solvers. It would let solvers query the information as they want, rows(vector_linear.A) or cols(vector_linear.A) or through scalar vector_linear.A[i, j]

matbesancon avatar Aug 30 '19 10:08 matbesancon

Do you know a solver which could take advantage of knowing the matrix structure ?

blegat avatar Aug 30 '19 14:08 blegat

When Hypatia is released, I will chime in here. Not quite ready to yet.

chriscoey avatar Aug 30 '19 14:08 chriscoey

From the Gurobi doc:

A brief note for users of the Gurobi MATLAB® and R interfaces: our interfaces to these languages are built around the assumption that you will use the rich matrix-oriented capabilities of the underlying languages to build your optimization models. Thus, our examples for these languages don’t attempt to show you how to build models. We have instead chosen to provide a few simple examples that demonstrate how to pass matrices into our interface.

With this example in R: https://www.gurobi.com/documentation/8.1/examples/mip_r.html

Where the matrix is passed to the model, instead of constructing row by row

matbesancon avatar Aug 30 '19 15:08 matbesancon

Is Gurobi actually using this matrix or is it just transforming it and for the algorithm it won't change anything whether you give this matrix or a VectorAffineFunction ? Is it just to gain model generation time or also solve time ?

blegat avatar Aug 30 '19 15:08 blegat

Not sure there is a way to know, I would say model generation time only.

matbesancon avatar Aug 30 '19 16:08 matbesancon

Do you know a solver which could take advantage of knowing the matrix structure ?

Tulip can. Actually, knowing the matrix structure is the main reason for having coded Tulip in the first place :)

Is it just to gain model generation time or also solve time ?

I'm only 99.99% sure it would affect model generation time only, not solve time.

mtanneau avatar Sep 07 '19 10:09 mtanneau

Is it just to gain model generation time or also solve time ?

Another point where this can affect performance is callbacks, if constraints defined in callbacks require moving around more data, it affects solving time.

matbesancon avatar Nov 12 '19 21:11 matbesancon

I wanted to get back to this issue and start thinking of possible solutions. Pinging @chriscoey now that Hypatia is released. Many solvers will take a pointer to an array of coefficients and a pointer to variable indices, in this form: https://github.com/scipopt/SCIP.jl/blob/master/src/MOI_wrapper/linear_constraints.jl#L12

or something close. This means that there will systematically be an unpacking of values into two arrays, from the one-array-of-struct current SAF model. If users provide the constraint in a form dot(coeff, vars) in set, then this means moving data around twice to go through the MOI layer.

SAF may be a better design for some operations, but at least having the capacity for solvers to specify that they take affine constraints in the pair-of-array form and not array-of-struct could make MOI closer to them.

Bridges from one form to the other should be doable

matbesancon avatar Jan 28 '21 18:01 matbesancon

An alternative to bridges would be to always use an interface when interacting with abstract linear operators, and both SAF and the pair-of-array type implementing this interface

matbesancon avatar Jan 28 '21 18:01 matbesancon

Relatedly, Convex.jl uses a lot of extended formulations which are similar to bridges in many ways, and I experimented with a bunch of "vector affine function" representations in https://github.com/jump-dev/Convex.jl/pull/393. Looking at that PR, the one I liked best was

struct SparseAffineOperation2{T}
    matrix::SparseMatrixCSC{T}
    vector::Vector{T}
end

struct SparseVAFTape2{T}
    operations::Vector{SparseAffineOperation2{T}}
    variables::Vector{MOI.VariableIndex}
end

where SparseAffineOperation2 is equivalent to a MOI.VectorAffineFunction, and SparseVAFTape2 lazily represents the product of a bunch of them (which can be materialized into a single MOI.VectorAffineFunction left-to-right at the end, allowing use of optimized sparse matmuls, instead of a bunch of unoptimized scalar operations).

ericphanson avatar Jan 28 '21 18:01 ericphanson

There are some related points in https://github.com/jump-dev/MathOptInterface.jl/issues/525

odow avatar Aug 12 '21 22:08 odow

I propose closing this in favor of #2180. I don't see us changing this anytime soon, and the complexity of adding a new function type to MOI makes me not want to try this.

odow avatar Sep 12 '23 22:09 odow

We actually have these implemented in DiffOpt: https://github.com/jump-dev/DiffOpt.jl/blob/master/src/utils.jl It's used to communicate the backward differentiation results

blegat avatar Sep 13 '23 08:09 blegat

Sure, but making them a first-class part of MOI and JuMP is a lot more work. And I don't really want more functions that solvers can opt-into with a bridge between them.

odow avatar Sep 13 '23 22:09 odow

This is an item in https://github.com/jump-dev/MathOptInterface.jl/issues/2180, so I think we can close this issue. I don't see us making progress anytime soon. It would be far too disruptive (even though it would be beneficial).

I'll give this a few days and then close. If anyone disagrees, please comment and we can discuss further. (Or re-open it if you're reading this from the far future.)

odow avatar Oct 23 '23 00:10 odow

Sorry, @matbesancon, yet another issue closed without resolution. At least this one doesn't have >100 comments...

odow avatar Oct 23 '23 19:10 odow

another addition to the List:tm:

matbesancon avatar Oct 23 '23 20:10 matbesancon