# Solve large number of linear inequalities with dynamic arrays

**URL:** <https://discourse.julialang.org/t/solve-large-number-of-linear-inequalities-with-dynamic-arrays/41821>\
**Category:** Optimization (Mathematical)\
**Tags:** question, jump, economics, arrays\
**Created:** [June 21, 2020, 6:46pm UTC](https://discourse.julialang.org/t/solve-large-number-of-linear-inequalities-with-dynamic-arrays/41821 "2020-06-21T18:46:02Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![LHR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lhr/32/15823_2.png) [@LHR](https://discourse.julialang.org/u/LHR)\
**Post date:** [June 21, 2020, 6:46pm UTC](https://discourse.julialang.org/t/solve-large-number-of-linear-inequalities-with-dynamic-arrays/41821/1 "2020-06-21T18:46:02Z")

</div>

Hi,

I am trying to make a script to remove dominated strategies from a pay-off matrix G for a game-theoretical problem that I have (2-player game, finite discrete strategy sets).

The goal is to automate the iterative elimination of dominated strategies, similar to, for example, the illustration here under 5.1: [http://faculty.haas.berkeley.edu/stadelis/Game%20Theory/econ160\_sect2a.pdf](http://faculty.haas.berkeley.edu/stadelis/Game%20Theory/econ160_sect2a.pdf)

At the heart the problem you have a system of linear inequalities and the goal is to check if there is a solution for p:

A\*p \> b

where A is m\*n matrix for which it needs to be checked if the product of A and probabilities p\_1, … ,p\_n is larger than the vector b. The vector b represents the pay-offs of the strategy which is checked for being dominated, the n vectors in A represent the pay-offs of (combinations of) strategies that potentially dominate b.

Since p\_1, … ,p\_n are probabilities we have the additional constraints that each p is:

**![0=pi=1](https://global.discourse-cdn.com/julialang/original/3X/4/a/4af5c6fb3704f75b2c567a87f400c0d10d083513.png "0 space less than equals space straight p subscript straight i space less than equals space 1")**

And that all p sum to 1:

![i=1npi=1](https://global.discourse-cdn.com/julialang/original/3X/3/d/3df8dbd2be408ae0aded3ffd445ca0b619ebf0bb.png "sum from space straight i space equals 1 space to straight n of straight p subscript straight i space equals space 1")

Now, if G is a pay-off matrix with s strategies, this problem has to be solved repeatedly for each strategy 1, … , s. So let’s say G is a 5x5 pay-off matrix, then each pay-off vector defined by strategy 1 to 5 (the vector b in the problem above) has to be checked against each (combination) of the other 4 vectors (the m\*n matrix A).

I implemented a solution for this using loops and JUMP. It works fine if G is a 11x11 or 12x12 matrix but computation time grows exponentially with larger matrizes. I need to be able to solve this for a matrix of at least size 21x21 which, with the current script, would take several weeks.

Here is a simplified working version of the code (_excluding_ a higher order function that takes the strategy sets of both player A and B into account):

```julia
# Packages
using Combinatorics
using JuMP
using GLPK

# Create example pay-off matrix G
G = rand(0:5, 5, 5)

# Loop through each strategy of G, then through (combinations of) possibly dominating strategies, pass to solveSysOfEq().
function returnDominatedStrategies(Mat::Array{Int64,2})
    Mat::Array{Int64,2} = transpose(Mat)
    # We loop through each strategy of the matrix.
    dominatedStrategies = zeros(0)
    for strat = 1:length(Mat[1,:])

        # ... and then again to get a separate count for the length indicator i. We have to do this for each strategy.
        for l in 1:(length(Mat[1,:])-1)

             # Here we calculate all combinations of potentially dominating strategies with length l.
             # [x for x in 1:(length(Mat[:,1])) if x != strat] will give us an integer list with all strategies of Mat except the strategy tested for being dominated (strat);
             for combo in combinations([x for x in 1:(length(Mat[1,:])) if x != strat],l)

               dominatedStrVect::Array{Int64,1} = Mat[:,strat] # Extract the pay-off vector of the strategy 'strat' tested for being dominated.
               dominatingStrVect::Array{Int64,2} = Mat[:,combo] # Extract the pay-off vectors of the (combination of) strategies tested for dominating 'strat' .

               # We call on the function that test if there is any combination of probability weights p_i that fulfills the inequality of the matrix product: dominatingStrVect*p_i > dominatedStrVect.
               solutionstatus = solveSysOfEq(dominatingStrVect, dominatedStrVect)

               # If there is a solution, i.e. a combination that dominates the current strategy, we add this strategy to 'dominatedStrategies', and jump to the loop for the next strategy.
               if solutionstatus == MOI.FEASIBLE_POINT
                   append!(dominatedStrategies, strat)
                   @goto end_l
               end
             end
        end
        @label end_l
    end
    return dominatedStrategies
end

# Check if there is a solution to the system of inequalities using Jump.
function solveSysOfEq(dominatingStrVect::Array{Int64,2}, dominatedStrVect::Array{Int64,1})

    model = Model(GLPK.Optimizer)

    A = dominatingStrVect
    b = dominatedStrVect

    @variable(model, 0 <= p[1:length(A[1,:])] <= 1)
    @constraint(model, sum(p[1:length(A[1,:])]) == 1)
    @constraint(model, A*p - b .>= 0.0000001)
    optimize!(model)

    status::MOI.ResultStatusCode = MOI.get(model, MOI.PrimalStatus())
    return status
end

# Function call
returnDominatedStrategies(G)

```

I know the “[Performance Tips](https://docs.julialang.org/en/v1/manual/performance-tips/index.html#man-performance-tips-1)”, and tried to implement some of the advice. However, there is a lot of information and I am not sure what the most promising way to proceed would be (also I didn’t understand everything tbh). This is the first time I need to optimize a piece of code for performance and the first time using Julia (switched from Python for this problem). What I would like to have is some guidance on how to optimize this. What would be the best strategy? What should I focus on?

**TLDR** : I have a system of linear inequalities with dynamically sized m\*n matrix A and vector b, and some additional constraints. Need to solve iteratively many, many times, at the moment using loops. What is the best strategy to optimize performance?

---

<div class="post-metadata">

**Author:** ![mtanneau](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mtanneau/32/17787_2.png) [@mtanneau](https://discourse.julialang.org/u/mtanneau)\
**Post date:** [June 24, 2020, 2:57pm UTC](https://discourse.julialang.org/t/solve-large-number-of-linear-inequalities-with-dynamic-arrays/41821/2 "2020-06-24T14:57:21Z")

</div>

## Some general tips

- before optimizing performance, you need to know what parts are making it slow (that’s true for any language, not just Julia).  
In Julia, you can use [ProfileView](https://github.com/timholy/ProfileView.jl) and [TimerOutputs](https://github.com/KristofferC/TimerOutputs.jl). If you’re new to Julia I would recommend ProfileView.
- Once you have profiled your code, look at where most time is spent, and try to optimize these parts first.  
They may not be the ones you think.
- There are several occurrences of `length(Mat[1, :])` → this will effectively extract the first row of `Mat`, then compute its length, every time that line is executed.  
Since the size of your matrix does not change, you can just compute `m, n = size(Mat)` at the beginning of the function, and use instead `for strat in 1:m`.
- Similarly, if a variable’s value does not change within a loop, avoid re-computing it. For instance, `dominatedStrVect` only depends on `strat`, so you only to compute it in the outer-most `for` loop.

## About the math

It is not possible to enforce _strict_ inequality constraints. In your code you use a tolerance; another way is to model it as follows:

\begin{array}{rl} \displaystyle \min\_{p, t, y} \ \ \ & t\\ s.t. \ \ \ & \displaystyle \sum\_{j=1}^{n} p\_{j} = 1\\ & y = b - Ap\\ & t \geq y\_{i} \ \ \ \forall i=1, ..., m\\ & p \geq 0 \end{array}

where y \in \mathbb{R}^{m}. Note that the upper bounds p\_{j} \leq 1 are not necessary here, as they follow from p \geq 0, e^{T}p = 1.

A strategy is dominated if and only if the objective value of this problem is negative.  
Indeed, if there exists \bar{p} such that A\bar{p} \> b, then \bar{y} = b - A \bar{p} has all its coordinates negative. Take \bar{t} = \max(\bar{y}\_{i}), and you have a feasible solution \bar{p}, \bar{t}, \bar{y} with objective value \bar{t} \< 0. Reciprocally, assume the optimal value of the above problem is p^{\*}, t^{\*}, y^{\*} with t^{\*} \< 0.  
It follows that A p^{\*} -b = - y^{\*} \> 0, i.e., A p^{\*} \> b.

Finally, if you’re trying to dominate strategy j\_{0}, just add the constraint p\_{j\_{0}} = 0 (by setting the corresponding upper bound to 0) to effectively remove it from the problem.

With this approach, you only need to formulate the LP once, then adjust the upper bounds and the right-hand side as needed.  
To do so, use the [`set_upper_bound`](https://jump.dev/JuMP.jl/stable/variables/#Variable-bounds-1) and [`set_normalized_rhs`](https://jump.dev/JuMP.jl/stable/constraints/#JuMP.set_normalized_rhs) functions in JuMP.  
That is much faster than re-creating the LP from scratch every time.

Finally, is it really necessary to loop through all possible supports (called `combo` in your code)? Isn’t it enough to just check dominance against a combination of all strategies that are not `strat`? You’re only returning a list of dominated strategies, not the dominating ones.  
That would require solving only _one_ LP for each strategy, rather than exponentially many.

---

<div class="post-metadata">

**Author:** ![LHR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lhr/32/15823_2.png) [@LHR](https://discourse.julialang.org/u/LHR)\
**Post date:** [June 30, 2020, 1:52pm UTC](https://discourse.julialang.org/t/solve-large-number-of-linear-inequalities-with-dynamic-arrays/41821/3 "2020-06-30T13:52:39Z")

</div>

Hi mtanneau,

thanks very much for the detailed response!

I am trying to implement your suggestions and hitting a roadblock with building the model. In particular, my code does not seem to be able to enforce the required strict inequality. Below a MWE with a) strict dominance of j\_3 over b, b) equality between j\_3 and b.

Any advice? What am I doing wrong?

```julia

using JuMP
using GLPK

##

# Check if there is a solution to the system of inequalities using Jump.
function solveSysOfEq(dominatingStrVect::Array{Float64,2}, dominatedStrVect::Array{Float64,1})

    model = Model(GLPK.Optimizer)

    m, n = size(dominatingStrVect)
    A = dominatingStrVect
    b = dominatedStrVect

    @variable(model, 0 <= p[1:n])
    @variable(model, y[1:m])
    @variable(model, t)

    @constraint(model, sum(p[1:n]) == 1)
    @constraint(model, y .== b - A*p)
    @constraint(model, [i=1:m], t >= y[i])

    @objective(model, Min, t)
    optimize!(model)

    #print(model)
    #print(value.(p))
    #print(value.(y))
    #print(value.(t))

    status::MOI.ResultStatusCode = MOI.get(model, MOI.PrimalStatus())
    return status
end

##

# a) A and b where j_3 strictly dominates b
function initAandb()
    dominatingStrVect::Array{Float64,2} = [5 2 5; 1 1 5; 3 1 2; 5 3 2; 2 1 2]
    dominatedStrVect::Array{Float64,1} = [4, 4, 1, 1, 1]
    return dominatingStrVect, dominatedStrVect
end

dominatingStrVect, dominatedStrVect = initAandb()
solveSysOfEq(dominatingStrVect, dominatedStrVect)

##

# b) A and b without strict dominance (equality between j_3 and b)
function initAandb()
    dominatingStrVect::Array{Float64,2} = [5 2 4; 1 1 4; 3 1 1; 5 3 1; 2 1 1]
    dominatedStrVect::Array{Float64,1} = [4, 4, 1, 1, 1]
    return dominatingStrVect, dominatedStrVect
end

dominatingStrVect, dominatedStrVect = initAandb()
solveSysOfEq(dominatingStrVect, dominatedStrVect)

##

```

Thanks in advance!

---

<div class="post-metadata">

**Author:** ![mtanneau](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mtanneau/32/17787_2.png) [@mtanneau](https://discourse.julialang.org/u/mtanneau)\
**Post date:** [July 1, 2020, 1:15am UTC](https://discourse.julialang.org/t/solve-large-number-of-linear-inequalities-with-dynamic-arrays/41821/4 "2020-07-01T01:15:59Z")

</div>

I don’t understand what you mean by

> my code does not seem to be able to enforce the required strict inequality

Once you have built the model, check the optimal value, e.g., with `objective_value(model)`, and if it is strictly negative, then the strategy is strictly dominated.  
If the objective is positive, the strategy is not dominated.