Explore PytOptInferface for writing matrix over native C API

[pz-max]

Background

Linopy either constructs a big matrix itself and stores it in a standardized file format such as the “LP” or the “MPS” file format or writes the matrix through direct Python solver interfaces. Linopy uses xarray and, therefore, has good memory and speed performance in creating the matrix/lp files. The syntax is also easier to read compared to Jump.

PyOptInterface was just released now (April 2024) by @metab0t and its benchmark shows that the speed in constructing the matrix is comparable to Jump.

Explore

Can we use PyOptInterface within Linopy to speed up the matrix construction?

[metab0t]

Thanks for your attention! @pz-max
The design of PyOptInterface is similar to the direct mode of JuMP.jl. It does not store the model internally and only records the mapping between variable/constraint and column/row index.
Do you plan to use PyOptInterface in place of direct Python solver interfaces in linopy?

[metab0t]

I made a minimal example here to build the model described in link.

On my computer, the output is:

Time to create linopy model: 0.29837608337402344
Time to convert to gurobipy model: 10.099567651748657
Time to convert to poi model: 10.64040493965149
Time to create poi model: 3.972501039505005

I am afraid that using PyOptInterface will not accelerate linopy because PyOptInterface does not have shortcut to create variables or constraints in bulk. I haven't profiled the code to find the time-consuming step, do you have any clues?

[pz-max]

@metab0t We would be keen to explore possibilities. Would you be free for a 30-minute exchange with @FabianHofmann (linopy maintainer) and me? If yes, could you ping me per email: max.parzen(at)openenergytransition.org with few time suggestions that are EU friendly? :)

[metab0t]

@pz-max I have sent you an email for arrangement of time. Have you received it?

[aurelije]

So this has potential of abstracting things in io module of linopy (where magic of translating linopy model into solver specific one happens)?

[metab0t]

After profiling, I think that the bottleneck is in the MatrixAccessor class to build matrix instead of gurobipy.

import time
from numpy import arange
from linopy import Model

import pyoptinterface as poi
from pyoptinterface import gurobi

def access_matrix(m):
    M = m.matrices
    vlabels = M.vlabels
    names = "x" + vlabels.astype(str).astype(object)
    vtypes = M.vtypes
    lb = M.lb
    ub = M.ub

    if m.is_quadratic:
        Q = M.Q

    c = M.c

    A = M.A
    clabels = M.clabels
    sense = M.sense
    b = M.b
    names = "c" + clabels.astype(str).astype(object)

def create_model(N):
    m = Model()
    x = m.add_variables(coords=[arange(N), arange(N)])
    y = m.add_variables(coords=[arange(N), arange(N)])
    m.add_constraints(x - y >= arange(N))
    m.add_constraints(x + y >= 0)
    m.add_objective((x * x).sum() + y.sum())
    return m

def create_poi_model(N):
    m = gurobi.Model()
    x = m.add_variables(range(N), range(N))
    y = m.add_variables(range(N), range(N))
    for i in range(N):
        for j in range(N):
            m.add_linear_constraint(x[i, j] - y[i, j], poi.Geq, i)
            m.add_linear_constraint(x[i, j] + y[i, j], poi.Geq, 0)
    expr = poi.ExprBuilder()
    poi.quicksum_(expr, x, lambda x: 2 * x)
    poi.quicksum_(expr, y)
    m.set_objective(expr)
    return m


N = 1000

t0 = time.time()
model = create_model(N)
t1 = time.time()
print('Time to create linopy model:', t1 - t0)

t0 = time.time()
access_matrix(model)
t1 = time.time()
print('Time to access matrices:', t1 - t0)

t0 = time.time()
raw_model = model.to_gurobipy()
t1 = time.time()
print('Time to convert to gurobipy model:', t1 - t0)

t0 = time.time()
raw_model = create_poi_model(N)
t1 = time.time()
print('Time to create poi model:', t1 - t0)

The output:

Time to create linopy model: 0.3041038513183594
Time to access matrices: 8.096210241317749
Time to convert to gurobipy model: 14.520848751068115
Time to create poi model: 6.902266025543213

More than half of the time is spent on constructing the matrix. The time spent on gurobipy (14.52 - 8.10 = 6.42s) is quite close to the time to build model in POI, so it is unlikely to be reduced significantly.

It also reflects the performance of POI where the sequential construction of POI is only slightly slower than the bulk construction with addMVar/addMConstr of gurobipy.

[pz-max]

Summary of exchange:

• POI has a very interesting design philosophy using C-APIs and schemas 🥇
• Constructing the Linopy model is fast, translating it to Gurobi via gurobipy is multiple times slower than construction
• With some work in POI, it could be possible to use the POI C API for the "translation" to remove the Python solver interface, increasing speed. This is useful not only for Linopy but also for many other tools.
• @metab0t will test stuff over the next 2 weeks

[metab0t]

@pz-max @FabianHofmann
Currently, the translation of linopy model to solver model includes two steps:

1. linopy model -> matrix form
2. matrix form -> gurobi model

In fact, the bottleneck is in the step 1. Gurobipy has done a good job in step 2 and swapping it with POI will not have significant influence on the performance even if POI implements adding variables/constraints in bulk.

What do you think about it?

[odow]

See also some related discussion in #207