FEAT: Add tools for repeated games

[QBatista]

Adds tools for repeated games including the outer hyperplane approximation described by Judd, Yeltekin, Conklin 2002.

The implementation is currently much slower than the one in Julia.

# Setup

pd_payoff = np.array([[9.0, 1.0], [10.0, 3.0]])

A = Player(pd_payoff)
B = Player(pd_payoff)

nfg = NormalFormGame((A, B))
rpd = RepeatedGame(nfg, 0.75)

• In Julia: 1.532440 seconds (1.17 M allocations: 87.199 MiB, 0.92% gc time)
• In Python: 1 loop, best of 3: 1min 8s per loop

[coveralls]

Coverage increased (+0.09%) to 93.75% when pulling fe1c996 on QBatista:repeated-games into f48e7d7 on QuantEcon:master.

[oyamad]

I haven't followed the details. Only small comments at the moment.

General comments:

• Try %prun to detect bottlenecks.
• Try method='interior-point' (ENH: added "interior-point" method for scipy.optimize.linprog scipy/scipy#7123) as a linprog option, which is available with the latest dev version of scipy.
• Consider providing a linprog_method option to outerapproximation, to be passed to linprog.
• Modification of the docstring in ce_util.py should belong to a separate PR.
• Did you compare gridmake with cartesian?

[oyamad]

I don't think all the routines should be imported. Just import important routines (perhaps RepeatedGame and outerapproximation?).

[oyamad]

The author (Chase Coleman) of the original code also should be here.

[oyamad]

Remove "form".

[oyamad]

Wouldn't C.copy() work?

[oyamad]

Wouldn't C[:] = Cnew work?

[oyamad]

Use and instead of &.

[oyamad]

Are these two blocks really necessary?

[QBatista]

They seem to be -- here is the output I get without them:

array([[ 10.00000001,   3.97266052],
       [ 10.00000001,   2.99999998],
       [  2.99999998,  10.00000001],
       [  3.97266052,  10.00000001],
       [  9.00000001,   8.99999999],
       [  8.99999999,   9.00000001],
       [  2.99999998,   3.00000001],
       [  2.99999998,   3.        ],
       [  2.99999999,   2.99999999],
       [  2.99999998,   3.        ],
       [  3.00000001,   2.99999998],
       [  3.        ,   2.99999998],
       [  2.99999999,   2.99999999],
       [  3.        ,   2.99999998],
       [  9.00000001,   9.        ],
       [  9.00000001,   9.        ],
       [  9.        ,   9.00000001],
       [  9.        ,   9.00000001]])

[oyamad]

I wonder why we have these duplications. We should look into the algorithm.

[cc7768]

In numpy version 1.13.1 they have updated the np.unique function to accept an axis argument -- Once this is released, we could do something like:

_, inds = np.unique(np.round(vertices, tol_int), axis=0, return_index=True)
vertices = vertices[inds, :]

or just

vertices = np.unique(np.round(vertices, tol_int), axis=0)

depending on whether we want the returned values to be rounded or not.

[oyamad]

What's wrong with using np.linspace?

[oyamad]

I would put this in repeated_game.py, as it's very specific to the code there.

[oyamad]

No need to compute all the pure Nash equilibria.

Try:

try:
    next(pure_nash_brute_gen(sg))
except StopIteration:
    raise ValueError('No pure action Nash equilibrium exists in stage game')

[coveralls]

Coverage increased (+0.003%) to 94.378% when pulling a60656a on QBatista:repeated-games into 21c1a94 on QuantEcon:master.

[coveralls]

Coverage increased (+0.02%) to 94.396% when pulling 4afe764 on QBatista:repeated-games into 21c1a94 on QuantEcon:master.

[QBatista]

Here is the comparison between cartesian and gridmake -- there does not seem to be much difference.

• With cartesian: 1min 22s ± 3.32 s per loop (mean ± std. dev. of 7 runs, 1 loop each)
• With gridmake : 1min 18s ± 1.7 s per loop (mean ± std. dev. of 7 runs, 1 loop each)

[oyamad]

It is interesting that cartesian is slower than gridmake (if the difference is not insignificant).
There has been a discussion that gridmake would be removed in favor of cartesian (sometime in the future); see #35.

[albop]

@oyamad and @QBatista : I refreshed my memories by looking at the implementation of cartesian and gridmake.
cartesian uses a function repeat_1d which was at the time faster than numpy's tile and repeat used in gridmake. Maybe this is not true anymore, which would explain why both methods have similar timings.
Actually, cartesian still does many memory allocating operations, so it is probably not near the optimal performance (would be easy to write a 2d benchmark). It doesn't do recursive calls like gridmake though, so I would expect it to be faster for higher dimensional products.

[albop]

Here a small gist to test the claims above, https://gist.github.com/albop/a4e6af9311fe9a9a392462ed757018bb
I compared cartesian, gridmake and a custom implementation in numba (2d and 3d only). I also added one version which assumes intervals are regularly spaced (so doesn't read memory).
Results below seem consistent: with what I said before. cartesian is faster than gridmake (but result is not the same as it is C-order instead of F-order), and the numba implementation is 4 times faster. Apparently computing coordinates doesn't improve anything.

# the product grids of K elements along each of the d dimensions is performed N times
test_1(N=100, d=2, K=1000) # cartesian : 0.499 seconds
test_2(N=100, d=2, K=1000) # gridmake : 0.609 seconds
test_3(N=100, d=2, K=1000) # numba : 0.114 seconds
test_4(N=100, d=2, K=1000) # numba regular : 0.141 seconds
test_1(N=100, d=3, K=100) # cartesian : 0.828 seconds
test_2(N=100, d=3, K=100) # gridmake : 1.012 seconds
test_3(N=100, d=3, K=100) # numba : 0.166 seconds
test_4(N=100, d=3, K=100) # numba regular : 0.164 seconds

[cc7768]

I haven't had a time to properly review this, but it looks like @oyamad has done a pretty thorough review. The code looks well organized and nicely written in the 10 minutes I spent reading through it.

One "Python" vs "Julia" comment that I have is in Julia you write functions that take a type as an argument, but in Python you typically attach these functions to the class itself as methods -- This means it might make sense to have some of these functions (in particular, any of the functions that take the RepeatedGame class as the first argument) become methods of the RepeatedGame class. This might mean that it makes sense to have nH be an argument to RepeatedGame. One argument against this is that it would make it so that the interface in Julia and Python be slightly different, so I'm not 100% convinced but wanted to bring it up. A possible way around this is just to have outerapproximation be the only additional method RepeatedGame and that would take the nH as an argument. Anyways, just thinking out loud. Am keen to hear what others think.

I'm not surprised that this code is a bit slower (this is precisely the type of example where one would expect Julia to perform better). It would be nice to investigate whether this could be sped up a little, but I don't think that is a first order priority for now.

[oyamad]

@cc7768 I agree that outerapproximation should better be a method of RepeatedGame (but with nH as an argument of the method).

For speedups we might implement a LP solver (by a simple simplex method) in Numba (as a medium term project).

[oyamad]

Should be something like:

def outerapproximation(..., linprog_method='simplex'):
    ...
    lpout = linprog(c, A_ub=A, b_ub=b, bounds=(lb, ub), method=linprog_method)

(scipy version 1.0.0 has not been released.)

[coveralls]

Changes Unknown when pulling 48021ff on QBatista:repeated-games into ** on QuantEcon:master**.

[QBatista]

@oyamad Here is a visualization of how the algorithm works on a few experiments: https://nbviewer.jupyter.org/github/QBatista/Notebooks/blob/master/JYC_Algo_Visualization.ipynb

An important observation is that the quality of the approximation is not weakly increasing with the number of gradients. For this game, choosing 32 gradients seems to give a better approximation than choosing 127 gradients which I suspect is because of the symmetry of payoffs. Additionally, it appears that the approximation is very sensitive to the geometry of the initial guess.
Absent a very good initial guess, the number of gradients will need to be much higher than the number of vertices to get a good approximation. I suspect this is because the algorithm runs into problems when the initial guess is not very good because of a combination of its monotonicity and outer approximation properties: once convex combinations of the approximation vertices are close to the true vertices, then even if the vertices are far from the boundary of the solution set, these vertices will not be able to move much. As such, a high number of gradients allows to increase the flexibility of the shape of the approximation, and therefore mitigate this issue.
Therefore, I think it would be better to intelligently pop similar vertices as we iterate once they reach a certain tolerance level by periodically checking for them in order to be able to set a large default number of gradients.
There are different options for doing this. A naive approach would be to substitute the two vertices by a convex combination of these two. The problem is that the outer approximation property would not necessarily be preserved. A better approach would be to substitute the approximation by a new polytope that ignores the edge connecting two similar vertices. This would imply losing the monotonicity property of the algorithm, which is not a real issue because we can only pop a finite number of vertices. Another point to consider is how to pop chains of vertices. One way of doing it would be to pop this group in one go, but this might be too aggressive if the true set has multiple close vertices. Instead, we could increase the frequency of checks depending on how many vertices were popped on the last iteration on the basis that the probability conditional on the number of vertices popped is likely an increasing function of the latter. Alternatively, we could increase the frequency based on the length and number of chains of similar vertices.
Based on the behavior of the algorithm with 128 gradients, it seems we could probably eliminate between 1/3 to 1/2 of the computations by avoiding unnecessary calls to the linear programming solver by implementing this type of feature while also improving the quality of the approximation.

[cc7768]

These visualizations are very cool. Nice work @QBatista. I think @thomassargent30 would be quite interested in seeing these visualizations.

Is there a reason you think that the value set corresponding to 32 subgradients looks much better than the one with 127? I kind of see what you mean since there are a few extra points between (3, 3) and the other two vertices, but they don't seem far off that line.

I agree the solution will have some dependence on the geometry of the initial guess (I think the difference between the set with 127 points and with 128 points actually illustrates this nicely -- The 127 point version doesn't necessarily have (0, 1), (1, 0), (0, -1), and (-1, 0) to work with which seem to be important components of this value set so it has to place more points along what should be the vertical line). This algorithm finds the smallest convex set (for a given set of subgradients!) that contains the fixed point of the B operator. A natural way to investigate this further would be to work on writing up the inner approximation which is also described by JYC. The fixed point of the B operator should lie between the inner and outer approximations -- My suspicion is that there are some very cool graphs you could draw that show the dependence on the initial geometry and how the inner and outer approximations differ for different games/geometries.algorithmically

I'm hesitant to be too aggressive with popping vertices. It does seem that not all vertices end up mattering very much, but I suspect it is hard to determine algorithmically which are the ones that we should keep. For example, in this game (-1, 0) and (0, -1) seem to be important vertices. I would be interested in seeing some of the tests you described above though as a proof-of-concept.

[oyamad]

@QBatista Animations look very nice.

Maybe we need "a better understanding of the manner in which extreme points of the equilibrium payoff set are generated" (Abreu and Sannikov, 2014). We should study Abreu and Sannikov (as we talked).

[mmcky]

@oyamad Is your last comment suggesting further development of this PR or a new project?

[oyamad]

@mmcky I suggested a new project. (I am afraid the JYC algorithm is too inefficient for pure Python/NumPy.)

[mmcky]

Hi @oyamad and @QBatista. Shall I close this PR? There is a lot of interesting work in this branch so I don't want to remove it without that work being documented. From reading this thread - we won't merge this at this time.