# Connecting a simple first-order solver to solve standard form linear program to JuMP

**URL:** <https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694>\
**Category:** Optimization (Mathematical)\
**Tags:** jump\
**Created:** [March 7, 2023, 8:56pm UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694 "2023-03-07T20:56:35Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Shuvomoy\_Das\_Gupta](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shuvomoy_das_gupta/32/10069_2.png) [@Shuvomoy\_Das\_Gupta](https://discourse.julialang.org/u/Shuvomoy_Das_Gupta)\
**Post date:** [March 7, 2023, 8:56pm UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/1 "2023-03-07T20:56:35Z")

</div>

Dear All,

To demonstrate the usage of first-order algorithms to solve linear programs in a graduate course that I am TAing (MIT 15.084: Nonlinear Optimization, link:[https://gitlab.com/vanparys-courses/15.084-sp23](https://gitlab.com/vanparys-courses/15.084-sp23)) I have written a simple Julia package called `SimplePDHG.jl` (https://github.com/Shuvomoy/SimplePDHG.jl) that is used to solve standard form linear programs of the form

\textrm{min}\_x c^{\top} x \quad \textrm{subject to } Ax = b, x\geq 0

using the vanilla PDHG algorithm (https://odlgroup.github.io/odl/math/solvers/nonsmooth/pdhg.html).

I was wondering if there is a way to create a generic linear programming model in `JuMP`, call the `SimplePDHG.jl` solver within `JuMP`, and then get the solution? Any documentation or other resources that I can look into regarding this will be very useful.

**Update**. Thanks to @odow’s help, the package is working fine now with model created in `JuMP`! Package available at: [https://github.com/Shuvomoy/SimplePDHG.jl](https://github.com/Shuvomoy/SimplePDHG.jl)

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [March 7, 2023, 8:59pm UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/2 "2023-03-07T20:59:31Z")

</div>

> I was wondering if there is a way to create a generic linear programming model in `JuMP` , call the `SimplePDHG.jl` solver within `JuMP` , and then get the solution?

The short answer is “yes.”

The longer answer is, “but you’ll need to write an MOI interface.”

Let me see if I can quickly mock something up.

The documentation to read is: [Implementing a solver interface · JuMP](https://jump.dev/JuMP.jl/stable/moi/tutorials/implementing/). But we have some tooling that should make this easier.

---

<div class="post-metadata">

**Author:** ![Shuvomoy\_Das\_Gupta](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shuvomoy_das_gupta/32/10069_2.png) [@Shuvomoy\_Das\_Gupta](https://discourse.julialang.org/u/Shuvomoy_Das_Gupta)\
**Post date:** [March 7, 2023, 9:06pm UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/3 "2023-03-07T21:06:06Z")

</div>

Great, thanks so much @odow , please let me know!

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [March 7, 2023, 9:38pm UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/4 "2023-03-07T21:38:26Z")

</div>

Try this:

```julia
module SimplePDHG

import MathOptInterface as MOI
import SparseArrays

"""
    solve_pdhg(
        A::SparseArrays.SparseMatrixCSC{Float64,Int},
        b::Vector{Float64},
        c::Vector{Float64},
    )::Tuple{MOI.TerminationStatusCode,Vector{Float64}}
"""
function solve_pdhg(
    A::SparseArrays.SparseMatrixCSC{Float64,Int},
    b::Vector{Float64},
    c::Vector{Float64},
)::Tuple{MOI.TerminationStatusCode,Vector{Float64}}
    x = fill(NaN, size(A, 2))
    return MOI.OTHER_ERROR, x
end

MOI.Utilities.@product_of_sets(RHS, MOI.Zeros)

const OptimizerCache = MOI.Utilities.GenericModel{
    Float64,
    MOI.Utilities.ObjectiveContainer{Float64},
    MOI.Utilities.VariablesContainer{Float64},
    MOI.Utilities.MatrixOfConstraints{
        Float64,
        MOI.Utilities.MutableSparseMatrixCSC{
            Float64,
            Int,
            MOI.Utilities.OneBasedIndexing,
        },
        Vector{Float64},
        RHS{Float64},
    },
}

function MOI.add_constrained_variables(
    model::OptimizerCache,
    set::MOI.Nonnegatives,
)
    x = MOI.add_variables(model, MOI.dimension(set))
    MOI.add_constraint.(model, x, MOI.GreaterThan(0.0))
    ci = MOI.ConstraintIndex{MOI.VectorOfVariables,MOI.Nonnegatives}(x[1].value)
    return x, ci
end

mutable struct Optimizer <: MOI.AbstractOptimizer
    x_primal::Dict{MOI.VariableIndex,Float64}
    termination_status::MOI.TerminationStatusCode

    function Optimizer()
        return new(Dict{MOI.VariableIndex,Float64}(), MOI.OPTIMIZE_NOT_CALLED)
    end
end

function MOI.is_empty(model::Optimizer)
    return isempty(model.x_primal) &&
        model.termination_status == MOI.OPTIMIZE_NOT_CALLED
end

function MOI.empty!(model::Optimizer)
    empty!(model.x_primal)
    model.termination_status = MOI.OPTIMIZE_NOT_CALLED
    return
end

function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.VectorAffineFunction{Float64}},
    ::Type{MOI.Zeros},
)
    return true
end

MOI.supports_add_constrained_variables(::Optimizer, ::Type{MOI.Reals}) = false

function MOI.supports_add_constrained_variables(
    ::Optimizer,
    ::Type{MOI.Nonnegatives},
)
    return true
end

function MOI.supports(
    ::Optimizer,
    ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}},
)
    return true
end

function MOI.optimize!(dest::Optimizer, src::MOI.ModelLike)
    cache = OptimizerCache()
    index_map = MOI.copy_to(cache, src)
    @assert all(iszero, cache.variables.lower)
    @assert all(==(Inf), cache.variables.upper)
    A = convert(
        SparseArrays.SparseMatrixCSC{Float64,Int},
        cache.constraints.coefficients,
    )
    b = cache.constraints.constants
    c = zeros(size(A, 2))
    offset = cache.objective.scalar_affine.constant
    for term in cache.objective.scalar_affine.terms
        c[term.variable.value] += term.coefficient
    end
    if cache.objective.sense == MOI.MAX_SENSE
        c *= -1
    end
    dest.termination_status, x_primal = solve_pdhg(A, b, c)
    for x in MOI.get(src, MOI.ListOfVariableIndices())
        dest.x_primal[x] = x_primal[index_map[x].value]
    end
    return index_map, false
end

function MOI.get(model::Optimizer, ::MOI.VariablePrimal, x::MOI.VariableIndex)
    return model.x_primal[x]
end

function MOI.get(model::Optimizer, ::MOI.ResultCount)
    return model.termination_status == MOI.OPTIMAL ? 1 : 0
end

function MOI.get(model::Optimizer, ::MOI.RawStatusString)
    return "$(model.termination_status)"
end

MOI.get(model::Optimizer, ::MOI.TerminationStatus) = model.termination_status

function MOI.get(model::Optimizer, ::MOI.PrimalStatus)
    if model.termination_status == MOI.OPTIMAL
        return MOI.FEASIBLE_POINT
    else
        return MOI.NO_SOLUTION
    end
end

MOI.get(model::Optimizer, ::MOI.DualStatus) = MOI.NO_SOLUTION

MOI.get(::Optimizer, ::MOI.SolverName) = "PDHG"

end

using JuMP
model = Model(SimplePDHG.Optimizer)
@variable(model, x >= 1)
@constraint(model, x <= 2)
@objective(model, Max, x + 1)
optimize!(model)
solution_summary(model; verbose = true)

```

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [March 7, 2023, 9:49pm UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/5 "2023-03-07T21:49:29Z")

</div>

cc @blegat, perhaps I should add this as a tutorial, with an implementation of simplex?

---

<div class="post-metadata">

**Author:** ![Shuvomoy\_Das\_Gupta](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shuvomoy_das_gupta/32/10069_2.png) [@Shuvomoy\_Das\_Gupta](https://discourse.julialang.org/u/Shuvomoy_Das_Gupta)\
**Post date:** [March 7, 2023, 10:03pm UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/6 "2023-03-07T22:03:17Z")

</div>

Great thanks so much @odow ! I will test it and let you know. Also, if I get time, I will try to do a step-by-step tutorial showing how to do that for my students and upload the video on Youtube.

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [March 7, 2023, 10:05pm UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/7 "2023-03-07T22:05:32Z")

</div>

Ah. There’s probably a little bug. I think I setup `Ax + b in {0}`, so the `b` might be the negative of what you expect.

---

<div class="post-metadata">

**Author:** ![Shuvomoy\_Das\_Gupta](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shuvomoy_das_gupta/32/10069_2.png) [@Shuvomoy\_Das\_Gupta](https://discourse.julialang.org/u/Shuvomoy_Das_Gupta)\
**Post date:** [March 7, 2023, 10:08pm UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/8 "2023-03-07T22:08:54Z")

</div>

Sounds good, I will try to implement it step by step and keep you posted if I run into issues. But the code that you provided is very helpful for me to get started. Thanks again @odow

---

<div class="post-metadata">

**Author:** ![blegat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/blegat/32/217090_2.png) [@blegat](https://discourse.julialang.org/u/blegat)\
**Post date:** [March 8, 2023, 8:59am UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/9 "2023-03-08T08:59:49Z")

</div>

Thanks, adding a tutorial for this would indeed be helpful. The code is short and concise but it’s difficult to get it right the first time.

---

<div class="post-metadata">

**Author:** ![Shuvomoy\_Das\_Gupta](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shuvomoy_das_gupta/32/10069_2.png) [@Shuvomoy\_Das\_Gupta](https://discourse.julialang.org/u/Shuvomoy_Das_Gupta)\
**Post date:** [March 8, 2023, 7:01pm UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/10 "2023-03-08T19:01:38Z")

</div>

Hi @odow , the code you provided is working great with minor modifications regarding the sign of `b`! The package (less than 350 lines of code) is available at [https://github.com/Shuvomoy/SimplePDHG.jl](https://github.com/Shuvomoy/SimplePDHG.jl), and can be installed by running

```julia
] add https://github.com/Shuvomoy/SimplePDHG.jl.git

```

Thanks again for your help, really appreciate it. The students will like it a lot!

---

<div class="post-metadata">

**Author:** ![A\_M](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/a_m/32/46038_2.png) [@A\_M](https://discourse.julialang.org/u/A_M)\
**Post date:** [August 10, 2023, 8:27pm UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/11 "2023-08-10T20:27:57Z")

</div>

@odow

I have an LP solver in C++, which can already run in Julia via this command.

objective, solution = sh.lp\_solve(solver, aᵢ, aⱼ, aᵥ, b, c)  
solver: solver options (dictionary)  
aᵢ, aⱼ, aᵥ: CSR representation for matrix A

The form that I solve is

\textrm{min}\_x c^{\top} x \quad \textrm{subject to } Ax \leq b

no x\>= 0, How can I modify your code snippet to consider this solver

Thanks

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [August 10, 2023, 8:35pm UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/12 "2023-08-10T20:35:30Z")

</div>

I haven’t tested, but you’d need something like:

```julia
MOI.Utilities.@product_of_sets(RHS, MOI.Nonpositives) # <-- change

const OptimizerCache = MOI.Utilities.GenericModel{
    Float64,
    MOI.Utilities.ObjectiveContainer{Float64},
    MOI.Utilities.VariablesContainer{Float64},
    MOI.Utilities.MatrixOfConstraints{
        Float64,
        MOI.Utilities.MutableSparseMatrixCSC{
            Float64,
            Int,
            MOI.Utilities.OneBasedIndexing,
        },
        Vector{Float64},
        RHS{Float64},
    },
}

mutable struct Optimizer <: MOI.AbstractOptimizer
    x_primal::Dict{MOI.VariableIndex,Float64}
    termination_status::MOI.TerminationStatusCode

    function Optimizer()
        return new(Dict{MOI.VariableIndex,Float64}(), MOI.OPTIMIZE_NOT_CALLED)
    end
end

function MOI.is_empty(model::Optimizer)
    return isempty(model.x_primal) &&
        model.termination_status == MOI.OPTIMIZE_NOT_CALLED
end

function MOI.empty!(model::Optimizer)
    empty!(model.x_primal)
    model.termination_status = MOI.OPTIMIZE_NOT_CALLED
    return
end

function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.VectorAffineFunction{Float64}},
    ::Type{MOI.Nonpostives}, # <-- change
)
    return true
end

function MOI.supports_constraint( # <-- new. No variable bounds
    ::Optimizer,
    ::Type{MOI.VariableIndex},
    ::Type{
        <:Union{
            MOI.LessThan{Float64},
            MOI.GreaterThan{Float64},
            MOI.EqualTo{Float64},
            MOI.Interval{Float64},
        },
    }
)
    return false
end

function MOI.supports(
    ::Optimizer,
    ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}},
)
    return true
end

function MOI.optimize!(dest::Optimizer, src::MOI.ModelLike)
    cache = OptimizerCache()
    index_map = MOI.copy_to(cache, src)
    @assert all(==(-Inf), cache.variables.lower) # <-- change
    @assert all(==(Inf), cache.variables.upper)
    A = convert(
        SparseArrays.SparseMatrixCSC{Float64,Int},
        cache.constraints.coefficients,
    )
    b = cache.constraints.constants
    c = zeros(size(A, 2))
    offset = cache.objective.scalar_affine.constant
    for term in cache.objective.scalar_affine.terms
        c[term.variable.value] += term.coefficient
    end
    if cache.objective.sense == MOI.MAX_SENSE
        c *= -1
    end
    dest.termination_status, x_primal = solve_pdhg(A, b, c)
    for x in MOI.get(src, MOI.ListOfVariableIndices())
        dest.x_primal[x] = x_primal[index_map[x].value]
    end
    return index_map, false
end

```

---

<div class="post-metadata">

**Author:** ![A\_M](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/a_m/32/46038_2.png) [@A\_M](https://discourse.julialang.org/u/A_M)\
**Post date:** [August 10, 2023, 8:39pm UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/13 "2023-08-10T20:39:35Z")

</div>

Thanks, I will test it out.

---

<div class="post-metadata">

**Author:** ![A\_M](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/a_m/32/46038_2.png) [@A\_M](https://discourse.julialang.org/u/A_M)\
**Post date:** [August 12, 2023, 9:50pm UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/15 "2023-08-12T21:50:35Z")

</div>

Thanks @odow; it worked well with a few modifications. How can I extract the Q matrix from this objective?

\min\_{x} \frac{1}{2} x^T Q x + c^T x \quad \text{subject to} \quad A x \leq b

Considering I have two solvers, lp\_solve and qp\_solve, I must differentiate LP vs QP and extract c for LP and (c, Q) for QP

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [August 13, 2023, 9:06pm UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/16 "2023-08-13T21:06:56Z")

</div>

You’ll need to add

```plaintext
function MOI.supports(
    ::Optimizer,
    ::MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}},
)
    return true
end

```

and then the function is in `cache.objective.scalar_quadratic`. But just as a list of terms, so you’ll need to loop over the `cache.objective.scalar_quadratic.quadratic_terms` to extract the coefficients.

---

<div class="post-metadata">

**Author:** ![A\_M](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/a_m/32/46038_2.png) [@A\_M](https://discourse.julialang.org/u/A_M)\
**Post date:** [August 14, 2023, 8:21pm UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/18 "2023-08-14T20:21:21Z")

</div>

I get this error

```julia
ERROR: LoadError: MathOptInterface.UnsupportedConstraint{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.LessThan{Float64}}: `MathOptInterface.ScalarAffineFunction{Float64}`-in-`MathOptInterface.LessThan{Float64}` constraint is not supported by the model.

```

but I added

```julia
function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.ScalarAffineFunction{Float64}},
    ::Type{MOI.LessThan{Float64}},
)
    return true
end

```

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [August 14, 2023, 8:30pm UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/19 "2023-08-14T20:30:52Z")

</div>

> [@odow](#):
>
> `MOI.Utilities.@product_of_sets(RHS, MOI.Nonpositives) # <-- change`

The `RHS` was setup such that it supported only VectorAffineFunction-in-Nonnegatives. Normally bridges should be able to reformulate that for you though. How did you initialize the model?

---

<div class="post-metadata">

**Author:** ![A\_M](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/a_m/32/46038_2.png) [@A\_M](https://discourse.julialang.org/u/A_M)\
**Post date:** [August 14, 2023, 8:46pm UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/20 "2023-08-14T20:46:31Z")

</div>

When I am trying to run, it takes forever from when I call !Optimize to get into it (I put a print statements with time stamp.)

```julia
println(Dates.format(now(), "HH:MM:SS"), " Solving")
optimize!(model)

function MOI.optimize!(dest::Optimizer, src::MOI.ModelLike)
    println(Dates.format(now(), "HH:MM:SS"), " Solver Caching ...")

```

I thought If I would add `add_bridges = false` it would help. Anyhow I changed my constraints so add\_bridges flag is accepted but still takes very long. I used Mosek and that time was very small.

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [August 14, 2023, 9:08pm UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/21 "2023-08-14T21:08:23Z")

</div>

Do you have a reproducible example?

---

<div class="post-metadata">

**Author:** ![A\_M](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/a_m/32/46038_2.png) [@A\_M](https://discourse.julialang.org/u/A_M)\
**Post date:** [August 14, 2023, 10:55pm UTC](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694/22 "2023-08-14T22:55:55Z")

</div>

Let me work on that.

If I also have \underline{x} \leq x \leq \overline{x} , what should I add in addition to `MOI.Utilities.@product_of_sets(RHS, MOI.Nonpositives)`

because the following didn’t work

```julia
function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.VariableIndex{Float64}},
    ::Type{MOI.Interval},
)
    return true
end

```

[Next page](https://discourse.julialang.org/t/connecting-a-simple-first-order-solver-to-solve-standard-form-linear-program-to-jump/95694.md?page=2)