Solve large number of linear inequalities with dynamic arrays

[LHR]

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

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:

And that all p sum to 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):

# 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”, 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?

[mtanneau]

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 and TimerOutputs. 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 and 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.

[LHR]

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?

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!

[mtanneau]

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.