Efficient Nash Equilibria using JuMP

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) :

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).

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

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

@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.

2 Likes

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

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.

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.

1 Like

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).

6 Likes

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

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.

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

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.

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!

1 Like

Use an affine expression as the objective instead:

    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.

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
2 Likes

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.

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

2 Likes