Abstract / alternative linear operator storage

[matbesancon]

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]

[blegat]

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

[chriscoey]

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

[matbesancon]

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

[blegat]

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 ?

[matbesancon]

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

[mtanneau]

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.

[matbesancon]

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]

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]

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

[ericphanson]

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 jump-dev/Convex.jl#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).

[odow]

There are some related points in #525