A code example on primal-dual algorithm for LP

Offtopic
1

Here is a simplest julia realization of the barrier method on solving LP.
The exemplary problem is

Min x[1] + 2 * x[2] + 3 * x[3]
s.t. sum(x) == 1
     x .>= 0

My code is

function gapfun(x, z) return transpose(x) * z end
function bfc(g, N) return g / N^2 end # barrier function coefficient
# data of the problem
N = 3
A = [1 1 1.]
c = [1, 2, 3.]
# initial trial primal-dual interior point (both-side feasible)
x = [0.1, 0.4, 0.5]
l = zeros(1) # it is 1-dim only
z = [1, 2, 3.]
# termination gap tolerance
TGT = 1e-6
for ite in 0:999999
    # initial gap
    g = gapfun(x, z)
    @info "ite = $ite ▶ lb = $(first(l)) < $(c' * x) = ub, gap = $g"
    g < TGT && (@info "▶ ▶ ▶ primal-dual algorithm converges"; break)
    BFC = bfc(g, N)
    # generate dir, i.e. Δ-*
    v = x .* z .- BFC
    ø, ⅁ = x ./ z, v ./ z
    dl = ((A .* transpose(ø)) * transpose(A)) \ (A * ⅁)
    dz = -transpose(A) * dl
    dx = -⅁ - ø .* dz
    # generate stepsize
    bx = dx .< 0
    bz = dz .< 0
    (any(bx) || any(bz)) || error("Abnormal: moving towards interior")
    ax = x[bx] ./ -dx[bx]
    az = z[bz] ./ -dz[bz]
    α = 0.999 * minimum([ax; az])
    # update
    @. x += α * dx
    @. l += α * dl
    @. z += α * dz
end
2 Likes
2

Is there are question here?

If it is intended as a tutorial-type post, I might ask that we spend our effort improving the existing JuMP documentation instead of writing posts like this for discourse. For example, there is:

2 Likes
3

It’s not all about software, sometimes algorithm themselves are worth looking into.
I’ll take a look about the existing framework.

4

This is a forum that is explicitly about the Julia programming language. Reflecting about algorithms and theory is of course useful, but perhaps a personal blog would be a more appropriate venue?

1 Like
5

I put it to the offtopic category.

Are there any well-developed blog platform? please inform me. I’m not familiar with the apps (I’m not a westerner).

6

In my opinion, forem.julialang.org is the exact place for this type of content: The Julia Forem: What it is, why we made one, and how to use it! - Julia Community 🟣

7

It’s fine for this time, but people who follow the Offtopic category do not expect a sudden surge in posts about mathematical programming algorithms or solver shortcomings.
As a forum moderator, I would kindly invite you to think about how you use this specific Discourse platform, and whether it is the best way for you to contribute to the Julia community as a whole (which is of course very happy to have you!).

For instance, documentation PRs may be welcome in the JuMP ecosystem, as pointed out by @odow, but they would need to be focused on what users need. To a certain extend, this forum has the same purpose, helping and orienting people. Unfortunately for me and you, what users need does not necessarily coincide with what we find interesting at a given moment. On the other hand, a blog or social network is a perfect format to tell the world what we find interesting!

Regarding the technical blogging stack, it really depends on your target audience, and I’m not an expert either. You can set up a simple static website with GitHub Pages + Jekyll, or you can publish on websites like https://medium.com/ which has lot of math-adjacent content.

I hope this does not deter you from contributing. It’s great that you have lots of energy and advanced technical knowledge to share. I’m just trying to steer it into something that benefits all of us :slight_smile:

3 Likes