# Efficient Nash Equilibria using JuMP

**URL:** <https://discourse.julialang.org/t/efficient-nash-equilibria-using-jump/57616>\
**Category:** Optimization (Mathematical)\
**Tags:** question, jump, performance\
**Created:** [March 20, 2021, 6:06pm UTC](https://discourse.julialang.org/t/efficient-nash-equilibria-using-jump/57616 "2021-03-20T18:06:40Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ian\_Slagle](https://avatars.discourse-cdn.com/v4/letter/i/b5a626/32.png) [@Ian\_Slagle](https://discourse.julialang.org/u/Ian_Slagle)\
**Post date:** [March 20, 2021, 6:06pm UTC](https://discourse.julialang.org/t/efficient-nash-equilibria-using-jump/57616/1 "2021-03-20T18:06:40Z")

</div>

Previously I have been using LRSLib (which calls the C library lrsnash) for computing Nash equilibria of two-player simultaneous move zero-sum games. However, I wanted to implement a pure Julia implementation for the sake of performance (and for ease of use, as LRSLib requires additional setup than Pkg). This is my implementation thus far, which works in the very limited test cases I have run it on (adapted from [https://people.csail.mit.edu/costis/6896sp10/lec2.pdf](https://people.csail.mit.edu/costis/6896sp10/lec2.pdf)) :

```julia
using JuMP
import GLPK

function nash(R::Array{Float64, 2})
    # https://people.csail.mit.edu/costis/6896sp10/lec2.pdf

    # Nash Equilibrium for row player
    model = Model(GLPK.Optimizer)
    @variable(model, z)
    @objective(model, Max, z)
    @variable(model, x[1:size(R)[1]])
    @constraint(model, x .>= 0)
    @constraint(model, x' * R .>= z * ones(size(R)[2])')
    @constraint(model, x' * ones(size(R)[1]) == 1)
    optimize!(model)
    z_return = value(z)
    x_return = JuMP.value.(x)
    
    # Nash Equilibrium for column player
    model = Model(GLPK.Optimizer)
    @variable(model, z_1)
    @objective(model, Min, z_1)
    @variable(model, y[1:size(R)[2]])
    @constraint(model, y .>= 0)
    @constraint(model, -y' * R' .>= z_1 * -ones(size(R)[1])')
    @constraint(model, y' * ones(size(R)[2]) == 1)
    optimize!(model)
    y_return = JuMP.value.(y)
    
    return z_return, x_return, y_return
end

```

Note that while the nash equilibrium value itself is the same for both models, I am computing both to get the mixed strategy for the column player (maybe not necessary?). For large matrices, this is outperforming LRSLib beautifully. The results below are from BenchmarkTools (mean for points, error bars are standard deviations. It’s more obvious in the next plot, but it should not be assumed that the distributions are normal, I just wanted to show some measure of variance).

 ![image](https://global.discourse-cdn.com/julialang/original/3X/e/b/eb92b667131c298d072d4ee366ae14f5dda67034.png)

However, for smaller matrices, the results are less than inspiring:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/1/6/1605e08ca72e312cd48cd360c74f1d3552664c98.png)

Here’s my benchmarks (note that LRSLib also requires integers):

```julia
@benchmark nash(R) setup=(R=rand($i, $i) .- 0.5)
@benchmark LRSLib.nashsolve(trunc.(Int64, R .* 1000), -trunc.(Int64, R .* 1000)) setup=(R=rand($i, $i) .- 0.5

```

Most of the matrices I input are small. How can I make this more efficient, hopefully to at least comparable to the LRSLib times? Is it domain-specific optimizations (in more fleshed-out implementations I have versions that solve 1xm, nx1, and 2x2 matrices without linear programming; maybe there are similar tricks for other sizes?)? Are there faster ways of writing linear programming (I am very new to the area)? Are there ways to write the Julia code specifically faster (ie, type-safety concerns, allocations, etc)? Any help on this would be much appreciated.

---

<div class="post-metadata">

**Author:** ![Ian\_Slagle](https://avatars.discourse-cdn.com/v4/letter/i/b5a626/32.png) [@Ian\_Slagle](https://discourse.julialang.org/u/Ian_Slagle)\
**Post date:** [March 21, 2021, 12:58am UTC](https://discourse.julialang.org/t/efficient-nash-equilibria-using-jump/57616/2 "2021-03-21T00:58:38Z")

</div>

Upon further review, it’s absolutely possible to do this in a single model.

```julia
function nash(R::Array{Float64, 2})
    # Set up model and payoff
    model = Model(GLPK.Optimizer)
    @variable(model, z)
    @objective(model, Max, z)
    
    # Solve for row player
    @variable(model, x[1:size(R)[1]])
    @constraint(model, x .>= 0)
    @constraint(model, x' * R .>= z)
    @constraint(model, x' * ones(size(R)[1]) == 1)
    
    # Solve for column player from row player
    @variable(model, y[1:size(R)[2]])
    @constraint(model, y .>= 0)
    @constraint(model, y' * ones(size(R)[2]) == 1)
    @constraint(model, y' * R' .== z)
    
    optimize!(model)
    
    return value(z), JuMP.value.(x), JuMP.value.(y)
end

```

Probably not the fastest way to do so, but it seems to work significantly better (I’ve also put this first plot in log scale to compare better). I think it might (?) be possible to solve it directly instead.

 ![image](https://global.discourse-cdn.com/julialang/original/3X/8/b/8b891896f96505c9fcad9f372ea77dd041aad13a.png)

 ![image](https://global.discourse-cdn.com/julialang/original/3X/e/6/e65d737bcf7d89ed2d6621d4b07b1618e5ccdc61.png)

I’m wondering if it’s possible to set some of the model up ahead of time as a constant? The size of R is variable, but if it comes to it I could set up all of the reasonable possibilities.

---

<div class="post-metadata">

**Author:** ![miles.lubin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/miles.lubin/32/279_2.png) [@miles.lubin](https://discourse.julialang.org/u/miles.lubin)\
**Post date:** [March 21, 2021, 1:24am UTC](https://discourse.julialang.org/t/efficient-nash-equilibria-using-jump/57616/3 "2021-03-21T01:24:44Z")

</div>

A few quick thoughts:

1. Try `direct_model(GLPK.Optimizer)` instead of `Model(GLPK.Optimizer)`. This disables fancy transformations that JuMP can perform and acts more as if you’re calling the GLPK API directly.

2. It should be possible to extract the solution for the column player by looking at the dual solution to the problem for the row-player problem. This would cut the work in half.

3. Write `@variable(model, x[1:size(R)[1]], lower_bound = 0)` and delete `@constraint(model, x .>= 0)`.

4. Write `x' * ones(size(R)[1])` as `sum(x)`.

---

<div class="post-metadata">

**Author:** ![Ian\_Slagle](https://avatars.discourse-cdn.com/v4/letter/i/b5a626/32.png) [@Ian\_Slagle](https://discourse.julialang.org/u/Ian_Slagle)\
**Post date:** [March 21, 2021, 6:45pm UTC](https://discourse.julialang.org/t/efficient-nash-equilibria-using-jump/57616/4 "2021-03-21T18:45:08Z")

</div>

Thanks so much for your help! I have implemented 3 and 4 and it has indeed gotten faster.

```julia
function nash(R::Array{Float64, 2})
    # Set up model and payoff
    model = Model(GLPK.Optimizer)
    @variable(model, z)
    @objective(model, Max, z)
    
    # Solve for row player
    @variable(model, x[1:size(R)[1]], lower_bound = 0.0)
    @constraint(model, x' * R .>= z)
    @constraint(model, sum(x) == 1.0)
    
    # Solve for column player from row player
    @variable(model, y[1:size(R)[2]], lower_bound = 0.0)
    @constraint(model, sum(y) == 1.0)
    @constraint(model, y' * R' .== z)
    
    optimize!(model)
    
    return value(z), JuMP.value.(x), JuMP.value.(y)
end

```

However, I am running into a few issues here.

> [@miles.lubin](#):
>
> Try `direct_model(GLPK.Optimizer)` instead of `Model(GLPK.Optimizer)` . This disables fancy transformations that JuMP can perform and acts more as if you’re calling the GLPK API directly.

My attempt to do this (singe `direct_model(GLPK.Optimizer())` for the third line returned this error:

```julia
The solver does not support an objective function of type MathOptInterface.SingleVariable.

```

when the objective (maximize z, which is admittedly a single variable) was defined.

> [@miles.lubin](#):
>
> It should be possible to extract the solution for the column player by looking at the dual solution to the problem for the row-player problem. This would cut the work in half.

This sounds great! Unfortunately, while I have come across the term ‘dual’ in my reading on linear programming, I’m still not familiar enough to construct this myself. Is there an example of using JuMP in this way that I could look at?

Again, thanks so much for your help!

---

<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 21, 2021, 7:09pm UTC](https://discourse.julialang.org/t/efficient-nash-equilibria-using-jump/57616/5 "2021-03-21T19:09:19Z")

</div>

Use an affine expression as the objective instead:

```nohighlight
    model = direct_model(GLPK.Optimizer())
    @variable(model, z)
    @objective(model, Max, 1.0 * z)

```

For the dual solution, take a look at the values returned by `shadow_price` and `reduced_cost`.

```nohighlight
function nash(R::Matrix{Float64})
    model = direct_model(GLPK.Optimizer())
    @variable(model, z)
    @objective(model, Max, 1.0 * z)
    @variable(model, x[1:size(R, 1] >= 0)
    @constraint(model, c1, x' * R .>= z)
    @constraint(model, c2, sum(x) == 1.0)
    optimize!(model)
    @show shadow_price.(c1)
    @show shadow_price(c2)
    @show reduced_cost.(x)
    return value(z), value.(x)
end

```

---

<div class="post-metadata">

**Author:** ![Ian\_Slagle](https://avatars.discourse-cdn.com/v4/letter/i/b5a626/32.png) [@Ian\_Slagle](https://discourse.julialang.org/u/Ian_Slagle)\
**Post date:** [March 21, 2021, 8:31pm UTC](https://discourse.julialang.org/t/efficient-nash-equilibria-using-jump/57616/6 "2021-03-21T20:31:19Z")

</div>

This is working beautifully. I went ahead and added some checks to save time on solutions with pure strategies, but even without this its still faster than LRSLib at the nxn matrices tested.

```julia
function strat_vec(l::Int64, i::Int64)
    vec_to_return = vec(zeros(l))
    @inbounds vec_to_return[i] = 1.0
    return vec_to_return
end

minmax(R::Array{Float64, 2}, m::Int64) = @inbounds mapreduce(x -> maximum(R[:, x]), min, 1:m)
maxmin(R::Array{Float64, 2}, n::Int64) = @inbounds mapreduce(x -> minimum(R[x, :]), max, 1:n)
findminmax(R::Array{Float64, 2}, n::Int64) = @inbounds strat_vec(n, argmax(map(x -> minimum(R[x, :]), 1:n)))
findmaxmin(R::Array{Float64, 2}, m::Int64) = @inbounds strat_vec(m, argmin(map(x -> maximum(R[:, x]), 1:m)))

function nash(R::Array{Float64, 2})
    n, m = size(R)
    
    # Check if we have to do linear programming
    n == 1 && return minimum(R), vec([1.0]), strat_vec(m, argmin(vec(R)))
    m == 1 && return maximum(R), strat_vec(n, argmin(vec(R))), vec([1.0])
    minmax(R, m) == maxmin(R, n) && return minmax(R, m), findminmax(R, n), findmaxmin(R, m)
    
    # Set up model and payoff
    model = direct_model(GLPK.Optimizer())
    @variable(model, z)
    @objective(model, Max, 1.0 * z)
    
    # Solve for row player
    @variable(model, x[1:n], lower_bound = 0.0)
    @constraint(model, c1, x' * R .>= z)
    @constraint(model, sum(x) == 1.0)
    
    optimize!(model)
    
    return value(z), value.(x), shadow_price.(c1)
end

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/3/8/38a1149ee2f511d5ab29fe1b0364f252d7a066ed.png)