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?