GitHub

Basis matrices

This tutorial was generated using Literate.jl. Download the source as a .jl file.

This tutorial explains how to query the basis of a linear program.

Setup

This tutorial uses the following packages:

using JuMP
import HiGHS

Standard form example

Consider the following example, which is from the Wikipedia article on Basic feasible solutions:

\[\begin{aligned} \max \; & 0 \\ \text{s.t.}\; & 1x_1 + 5x_2 + 3x_3 + 4x_4 + 6x_5 = 14 \\ & 0x_1 + 1x_2 + 3x_3 + 5x_4 + 6x_5 = 7 \\ & x_i \ge 0, \forall i = 1,\ldots,5. \end{aligned}\]

The A matrix is:

A = [1 5 3 4 6; 0 1 3 5 6]
2×5 Matrix{Int64}:
 1  5  3  4  6
 0  1  3  5  6

and the right-hand side b vector is:

b = [14, 7]
2-element Vector{Int64}:
 14
  7

We can create and optimize the problem in the standard form:

n = size(A, 2)
model = Model(HiGHS.Optimizer)
set_silent(model)
@variable(model, x[1:n] >= 0)
@constraint(model, A * x == b)
optimize!(model)
assert_is_solved_and_feasible(model)

This has a solution:

value.(x)
5-element Vector{Float64}:
 0.0
 2.0
 0.0
 1.0
 0.0

Query the basis status of a variable using MOI.VariableBasisStatus:

get_attribute(x[1], MOI.VariableBasisStatus())
NONBASIC_AT_LOWER::BasisStatusCode = 2

the result is a MOI.BasisStatusCode. Query all of the basis statuses with the broadcast get_attribute.(:

get_attribute.(x, MOI.VariableBasisStatus())
5-element Vector{MathOptInterface.BasisStatusCode}:
 NONBASIC_AT_LOWER::BasisStatusCode = 2
 BASIC::BasisStatusCode = 0
 NONBASIC_AT_LOWER::BasisStatusCode = 2
 BASIC::BasisStatusCode = 0
 NONBASIC_AT_LOWER::BasisStatusCode = 2

For this problem, the values are either MOI.NONBASIC_AT_LOWER or MOI.BASIC. All of the MOI.NONBASIC_AT_LOWER variables have a value at their lower bound. The MOI.BASIC variables correspond to the columns of the optimal basis.

Get the columns using:

indices = get_attribute.(x, MOI.VariableBasisStatus()) .== MOI.BASIC
5-element BitVector:
 0
 1
 0
 1
 0

Filter the basis matrix from A:

B = A[:, indices]
2×2 Matrix{Int64}:
 5  4
 1  5

Since the basis matrix is non-singular, solving the system Bx = b must yield the optimal primal solution of the basic variables:

B \ b
2-element Vector{Float64}:
 2.0
 0.9999999999999998
value.(x[indices])
2-element Vector{Float64}:
 2.0
 1.0

A more complicated example

Often, you may want to work with the basis of a model that is not in a nice standard form. For example:

model = Model(HiGHS.Optimizer)
set_silent(model)
@variable(model, x >= 0)
@variable(model, 0 <= y <= 3)
@variable(model, z <= 1)
@objective(model, Min, 12x + 20y - z)
@constraint(model, c1, 6x + 8y >= 100)
@constraint(model, c2, 7x + 12y >= 120)
@constraint(model, c3, x + y <= 20)
optimize!(model)
assert_is_solved_and_feasible(model)

A common way to query the basis status of every variable is:

v_basis = Dict(
    xi => get_attribute(xi, MOI.VariableBasisStatus()) for
    xi in all_variables(model)
)
Dict{VariableRef, MathOptInterface.BasisStatusCode} with 3 entries:
  y => BASIC
  x => BASIC
  z => NONBASIC_AT_UPPER

Despite the model having three constraints, there are only two basic variables. Since the basis matrix must be square, where is the other basic variable?

The answer is that solvers will reformulate inequality constraints:

\[A x \le b\]

into the system:

\[A x + Is = b\]

Thus, for every inequality constraint there is a slack variable s.

Query the basis status of the slack variables associated with a constraint using MOI.ConstraintBasisStatus:

c_basis = Dict(
    ci => get_attribute(ci, MOI.ConstraintBasisStatus()) for ci in
    all_constraints(model; include_variable_in_set_constraints = false)
)
Dict{ConstraintRef{Model, C, ScalarShape} where C, MathOptInterface.BasisStatusCode} with 3 entries:
  c1 : 6 x + 8 y ≥ 100  => NONBASIC
  c3 : x + y ≤ 20       => BASIC
  c2 : 7 x + 12 y ≥ 120 => NONBASIC

Thus, the basis is formed by x, y, and the slack associated with c3.

A simple way to get the A matrix of an unstructured linear program is with lp_matrix_data:

matrix = lp_matrix_data(model)
matrix.A
3×3 SparseArrays.SparseMatrixCSC{Float64, Int64} with 6 stored entries:
 6.0   8.0   ⋅ 
 7.0  12.0   ⋅ 
 1.0   1.0   ⋅

You can check the permutation of the rows and columns using

matrix.variables
3-element Vector{VariableRef}:
 x
 y
 z

and

matrix.affine_constraints
3-element Vector{ConstraintRef}:
 c1 : 6 x + 8 y ≥ 100
 c2 : 7 x + 12 y ≥ 120
 c3 : x + y ≤ 20

We can construct the slack column associated with c3 as:

s_column = zeros(size(matrix.A, 1))
s_column[3] = 1.0
1.0

The full basis matrix is therefore:

B = hcat(matrix.A[:, [1, 2]], s_column)
3×3 SparseArrays.SparseMatrixCSC{Float64, Int64} with 7 stored entries:
 6.0   8.0   ⋅ 
 7.0  12.0   ⋅ 
 1.0   1.0  1.0

lp_matrix_data returns separate vectors for the lower and upper row bounds. Convert to a single right-hand side vector by taking the finite elements:

b = ifelse.(isfinite.(matrix.b_lower), matrix.b_lower, matrix.b_upper)
3-element Vector{Float64}:
 100.0
 120.0
  20.0

Solving the Basis system as before yields:

B \ b
3-element Vector{Float64}:
 14.999999999999995
  1.250000000000004
  3.75

which is the value of x, y, and the slack associated with c3.

Identifying degenerate variables

Another common task is identifying degenerate variables. A degenerate variable is a basic variable that has an optimal value at its lower or upper bound.

Here is a function that computes whether a variable is degenerate:

function is_degenerate(x)
    if get_attribute(x, MOI.VariableBasisStatus()) == MOI.BASIC
        return (has_lower_bound(x) && ≈(value(x), lower_bound(x))) ||
               (has_upper_bound(x) && ≈(value(x), upper_bound(x)))
    end
    return false
end
is_degenerate (generic function with 1 method)

A simple example of a linear program with a degenerate solution is:

A, b, c = [1 1; 0 1], [1, 1], [1, 1]
model = Model(HiGHS.Optimizer);
set_silent(model)
@variable(model, x[1:2] >= 0)
@objective(model, Min, c' * x)
@constraint(model, A * x == b)
optimize!(model)
degenerate_variables = filter(is_degenerate, all_variables(model))
1-element Vector{VariableRef}:
 x[1]

The solution is degenerate because:

value(x[1])
-0.0

and

get_attribute(x[1], MOI.VariableBasisStatus())
BASIC::BasisStatusCode = 0