Highs_getRowsByRange is slow

[jd-lara]

Recently I noticed that the numerical checks in PowerSimulations.jl since it takes a long time to run with HiGHS w.r.t to Gurobi or Xpress.

This is the loop we use is something like this

function get_constraint_numerical_bounds(model::OperationModel)
    if !is_built(model)
        error("Model not built, can't calculate constraint numerical bounds")
    end
    bounds = ConstraintBounds()
    for (const_key, constraint_array) in get_constraints(get_optimization_container(model))
        # TODO: handle this at compile and not at run time
        if isa(constraint_array, SparseAxisArray)
            for idx in eachindex(constraint_array)
                constraint_array[idx] == 0.0 && continue
                con_obj = JuMP.constraint_object(constraint_array[idx])
                update_coefficient_bounds(bounds, con_obj, (const_key, idx))
                update_rhs_bounds(bounds, con_obj, (const_key, idx))
            end
        else
            for idx in Iterators.product(constraint_array.axes...)
                !isassigned(constraint_array, idx...) && continue
                con_obj = JuMP.constraint_object(constraint_array[idx...])
                update_coefficient_bounds(bounds, con_obj, (const_key, idx))
                update_rhs_bounds(bounds, con_obj, (const_key, idx))
            end
        end
    end
    return bounds
end

[odow]

To clarify, what type of constraints is this iterating over?

[jd-lara]

To clarify, what type of constraints is this iterating over?

All supported ones in direct mode. Nothing esoteric

[odow]

I think we decided that this was because it queries the constraint matrix:

HiGHS.jl/src/MOI_wrapper.jl

Lines 1662 to 1708 in 1d192e2

	function MOI.get(
	model::Optimizer,
	::MOI.ConstraintFunction,
	c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},S},
	) where {S<:_SCALAR_SETS}
	r = row(model, c)
	num_row = Ref{HighsInt}(0)
	num_nz = Ref{HighsInt}(0)
	matrix_start = Ref{HighsInt}(0)
	lower, upper = Ref{Cdouble}(), Ref{Cdouble}()
	_ = Highs_getRowsByRange(
	model,
	r,
	r,
	num_row,
	lower,
	upper,
	num_nz,
	matrix_start,
	C_NULL,
	C_NULL,
	)
	matrix_index = Vector{HighsInt}(undef, num_nz[])
	matrix_value = Vector{Cdouble}(undef, num_nz[])
	ret = Highs_getRowsByRange(
	model,
	r,
	r,
	num_row,
	lower,
	upper,
	num_nz,
	matrix_start,
	matrix_index,
	matrix_value,
	)
	_check_ret(ret)
	return MOI.ScalarAffineFunction(
	MOI.ScalarAffineTerm{Float64}[
	MOI.ScalarAffineTerm(
	val,
	model.variable_info[CleverDicts.LinearIndex(col + 1)].index,
	) for (col, val) in zip(matrix_index, matrix_value) if !iszero(val)
	],
	0.0,
	)
	end

One option would be to use a cache instead of direct_model, so something like model = Model(HiGHS.Optimizer; add_bridges = false).

I don't think there's an easy way to speed this up in HiGHS.jl.

[odow]

Confirm that the issue is O(N) behavior in Highs_getRowsByRange and that we need to loop over every constraint:

import HiGHS
import Gurobi
import MathOptInterface as MOI
import Plots

function main(model, N)
    x = MOI.add_variables(model, N)
    ci = Any[]
    for i in 2:N
        f = 1.0 * x[i] + 1.0 * x[i-1]
        push!(ci, MOI.add_constraint(model, f, MOI.GreaterThan(0.0)))
    end
    return [MOI.get(model, MOI.ConstraintFunction(), ci_i) for ci_i in ci]
end

function time(model, x)
    main(model(), first(x))  # Precompile
    y = Float64[]
    for N in x
        GC.gc()
        push!(y, @elapsed main(model(), N))
    end
    return y
end
    
x = [1_000, 2_000, 4_000, 8_000, 10_000, 12_000, 14_000, 16_000];
y_grb = time(Gurobi.Optimizer, x);
y_highs = time(HiGHS.Optimizer, x);
y_moi = time(MOI.Utilities.Model{Float64}, x);
plt = Plots.plot(; xlabel = "N", ylabel = "Time [sec]")
Plots.plot!(plt, x, y_grb; label = "Gurobi")
Plots.plot!(plt, x, y_highs; label = "HiGHS")
Plots.plot!(plt, x, y_moi; label = "Utilities.Model")

[odow]

@jajhall shall we mark this as expected behavior, or can we make querying rows more efficient? I assume that you store things column-wise.

[jajhall]

Once Highs::run() is called, the matrix is held column-wise, so extracting rows one-by-one will be slow. In this case I'm slightly puzzled since, as the constraints are being added row-by-row, the matrix will be held row-wise. I don't see anything being done to cause HiGHS to transpose the representation. I'll try to reproduce it sometime