Dawen Wu
email: dawen.wu@centralesupelec.fr

Abdel Lisser
email: abdel.lisser@l2s.centralesupelec.fr

Improved Saddle Point Prediction in Stochastic Two-Player Zero-Sum Games with a Deep Learning Approach

Keywords: Stochastic two-player zero-sum games, Saddle points, Neurodynamic optimization, Deep learning, Ordinary differential equations

In this paper, we propose a novel deep learning approach for predicting saddle points in stochastic twoplayer zero-sum games. Our method combines neurodynamic optimization and deep neural networks. First, we model the stochastic two-player zero-sum game as an ordinary differential equation (ODE) system using neurodynamic optimization. Second, we develop a neural network to approximate the solution to the ODE system, which includes the saddle point prediction for the game problem. Third, we introduce a specialized algorithm for training the neural network to enhance the accuracy of the saddle point prediction. Our experiments demonstrate that our model outperforms existing approaches, yielding faster convergence and more accurate saddle point predictions.

Introduction

In non-cooperative game theory, two-player zero-sum games are among the most fundamental models, featuring two players, each with finite actions, and their payoffs add up to zero. A saddle point of such a game represents a situation in which no player can increase his/her payoff by unilaterally changing his/her strategy. von Neumann (1928) proved that there always exists a saddle point for any two-player zero-sum 5 game. Later, [START_REF] Nash | Equilibrium points in n-person games[END_REF] extended this result and proved that there always exists a Nash equilibrium for any n-player general-sum game with finite actions. [START_REF] Charnes | Constrained games and linear programming[END_REF] studied the two-player zero-sum games with linear constraints, which can be formulated as linear programming problems. Recently, [START_REF] Singh | A second-order cone programming formulation for two player zero-sum games with chance constraints[END_REF] studied a stochastic version of the two-player zero-sum games, namely stochastic two-player zero-sum games.

They demonstrated that a saddle point exists if the random vectors defining stochastic linear constraints 10 follow elliptically symmetric distributions.

Neurodynamic optimization refers to a method that uses a first-order ODE system to model a nonlinear optimization problem. [START_REF] Hopfield | neural" computation of decisions in optimization problems[END_REF] proposed the Hopfield network to solve the well-known traveling salesman problem. The Hopfield network was extended to solve nonlinear convex programming by introducing a penalty parameter [START_REF] Kennedy | Neural networks for nonlinear programming[END_REF]. However, this method struggles to achieve an optimal solution because the true minimizer can only be achieved when the penalty parameter goes to infinity. Since then, various types of neurodynamic optimization methods have been proposed without using the penalty parameter, to solve various optimization problems, e.g., quadratic programming [START_REF] Xia | A recurrent neural network for solving linear projection equations[END_REF], nonlinear projection equations [START_REF] Xia | A new neural network for solving nonlinear projection equations[END_REF], second-order cone programming [START_REF] Nazemi | A new collaborate neuro-dynamic framework for solving convex second order cone programming problems with an application in multi-fingered robotic hands[END_REF], non-smooth optimization problems [START_REF] Qin | A two-layer recurrent neural network for nonsmooth convex optimization problems[END_REF], pseudoconvex optimization problems [START_REF] Xu | A neurodynamic approach to nonsmooth constrained 495 pseudoconvex optimization problem[END_REF]. [START_REF] Wu | A dynamical neural network approach for solving stochastic two-player zerosum games[END_REF] proposed a neurodynamic optimization approach to model the stochastic two-player zero-sum game problem using Karush-Kuhn-Tucker conditions.

Deep learning is a type of machine learning that uses deep neural networks, comprising many layers of interconnected nodes, to learn complex patterns and relationships in data. Thanks to the exponential growth of data and computing resources in recent years, deep learning has been applied successfully to a wide range of fields, including image processing [START_REF] Raza | Deep learning for computer vision: A comprehensive review[END_REF], natural language processing [START_REF] Tan | Deep learning for natural language processing: A review[END_REF], bioinformatics [START_REF] Hu | Deep learning in bioinformatics: A comprehensive survey[END_REF], game theory (Wu & Lisser, 2022d,c), and operations research (Wu & Lisser, 2022a;[START_REF] Nair | Solving mixed integer programs using neural networks[END_REF][START_REF] Bengio | Machine learning for combinatorial optimization: a methodological tour d'horizon[END_REF].

Research on utilizing neural networks to solve differential equations dates back to the 1990s, with the initial concept involving training a neural network to satisfy both a given differential equation and its boundary conditions [START_REF] Dissanayake | Neural-network-based approximations for solving partial differential equations[END_REF]. Advances in this field include the development of construction methods that inherently satisfy initial and boundary conditions, with applications in ordinary and partial differential equations [START_REF] Lagaris | Artificial neural networks for solving ordinary and partial differential equations[END_REF][START_REF] Lagaris | Neural-network methods for boundary value problems with irregular boundaries[END_REF][START_REF] Mcfall | Artificial neural network method for solution of boundary value problems with exact satisfaction of arbitrary boundary conditions[END_REF]. The advent of deep learning has revitalized the use of neural networks for solving differential equations, enabling researchers to address the challenge of solving high-dimensional nonlinear PDEs [START_REF] Han | Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations[END_REF][START_REF] Yu | The deep ritz method: a deep learning-based numerical algorithm for solving variational problems[END_REF][START_REF] Han | Solving high-dimensional partial differential equations using deep learning[END_REF] and leading to significant developments such as the introduction of physics-informed neural networks (PINNs) (Raissi et al., 2019a). PINNs have found successful applications across various fields, including computational mechanics (Anitescu et al., 2019;[START_REF] Samaniego | An energy approach to the solution of partial differential equations in computational mechanics via machine learning: Concepts, implementation and applications[END_REF], and have been modified to accommodate a range of problem scenarios [START_REF] Wang | Understanding and mitigating gradient flow pathologies in physics-informed neural networks[END_REF][START_REF] Lu | Physicsinformed neural networks with hard constraints for inverse design[END_REF][START_REF] Zhang | Learning in modal space: Solving timedependent stochastic pdes using physics-informed neural networks[END_REF]. The rapid progress in this research direction has been facilitated by automatic differentiation tools that simplify the calculation of derivatives, gradients, and Jacobian matrices (Baydin et al., 2018;Paszke et al., 2019a), and it has been demonstrated that the neural network approximator converges to the PDE solution as the number of hidden units increases [START_REF] Sirignano | Dgm: A deep learning algorithm for solving partial differential equations[END_REF]. Currently, software packages are available that utilize deep learning methods to solve differential equations (Lu et al., 2021a;[START_REF] Chen | Neurodiffeq: A python package for solving differential equations with neural networks[END_REF]. This paper addresses the saddle point problem in stochastic two-player zero-sum games [START_REF] Singh | A second-order cone programming formulation for two player zero-sum games with chance constraints[END_REF]. Previous research reformulated this problem as a problem of solving the state solution of an ODE system [START_REF] Wu | A dynamical neural network approach for solving stochastic two-player zerosum games[END_REF]. However, the existing solution method demands solving all intermediate states in the ODE system to arrive at the final prediction, which can lead to computationally intensive and timeconsuming processes. To overcome these challenges, this paper aims to enhance the saddle point prediction problem by incorporating advanced deep learning techniques for solving differential equations [START_REF] Raissi | Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations[END_REF][START_REF] Sirignano | Dgm: A deep learning algorithm for solving partial differential equations[END_REF]. By developing an efficient approach to the stochastic two-player zero-sum game problem, we address the current limitations and open new possibilities for solving various game-theoretic problems. The proposed approach can encourage researchers to explore the use of deep learning techniques in tackling other complex game-theoretic scenarios, contributing to advancements in the field of game theory and optimization.

It is important to note that the previous convolutional neural network (CNN) approach (Wu & Lisser, 2022d,c) employed to address two-player zero-sum games is not directly applicable to the stochastic twoplayer zero-sum game considered in this paper. The aforementioned papers focused on classical two-player zero-sum games without random variables, which could be represented by a matrix denoting payoffs and subsequently used as input for the CNN. In contrast, the stochastic two-player zero-sum game examined in this paper incorporates random variables within the constraints, preventing a direct input of the game into the CNN.

Key contributions

Our key contributions can be summarized as follows:

• We propose a novel method that combines neurodynamic optimization and deep learning to solve the stochastic two-player zero-sum game. To the best of our knowledge, this is the first time deep learning has been applied to this problem.

• Our work transforms the stochastic two-player zero-sum game into a neural network training problem, allowing us to solve the problem without the need for a standard numerical integration solver.

• Our method outperforms the state-of-the-art method presented in [START_REF] Wu | A dynamical neural network approach for solving stochastic two-player zerosum games[END_REF]. In the specific game examples given in the experimental section, our method achieves an objective value of 5.48 and maximum violation of constraints of 0.0827, respectively, outperforming the method of [START_REF] Wu | A dynamical neural network approach for solving stochastic two-player zerosum games[END_REF].

Paper outline

The remaining sections are organized as follows. The background knowledge necessary to understand this paper is provided in Section 2, including an introduction to the stochastic two-player zero-sum game, neurodynamic optimization, and numerical integration methods. Section 3 describes our proposed method and how it solves the game problem. Section 4 presents experimental results and a discussion in which we compare our method with the numerical integration methods. Section 5 summarizes this paper and outlines future directions. Index set of constraints for player 2.

Notation

α 1 = α 1 k k∈J1
Confidence levels for player 1's constraints.

α 2 = α 2 l l∈J2
Confidence levels for player 2's constraints.

G(α)

Stochastic two-player zero-sum game with confidence level α; α = (α 1 , α 2 ).

(x * ∈ R m , y * ∈ R n ) Saddle point of G(α). (x ∈ R m , ŷ ∈ R n ) Predicted saddle point. Ellip d (µ, Σ, φ) d-dimensional

Preliminaries

The stochastic two-player zero-sum game is introduced in Section 2.1. The neurodynamic optimization method, which models the stochastic two-player zero-sum game as an ODE system, is presented in Section 2.2.

Numerical integration methods, which are existing approaches for solving the game problem, are described in Section 2.3.

Stochastic two-player zero-sum game

A two-player zero-sum game problem with linear constraints introduced by [START_REF] Charnes | Constrained games and linear programming[END_REF] is characterized as follows.

• The game has two players in total, and players 1 and 2 have pure strategy sets of M = {1, . . . , m} and N = {1, . . . , n}, respectively. m and n denote the number of pure strategies of players 1 and 2, respectively.

• The payoff of the game is represented by a matrix A = (a ij ) ∈ R m×n . When player 1 chooses action i ∈ M and player 2 chooses action j ∈ N , players 1 and 2 receive payoffs of a ij and -a ij , respectively.

• x ∈ X and y ∈ Y are mixed strategies of players 1 and 2, respectively.

X = x ∈ R m | Bx ≤ b, 1 T m x = 1, x ≥ 0 and Y = y ∈ R n | Dy ≤ d, 1 T n y = 1, y ≥ 0 denote the feasible sets, where B ∈ R p×m , b ∈ R p , D ∈ R q×n , d ∈ R q , 1 T m = [1, 1, . . . , 1] ∈ R m , and 1 T n = [1, 1, . . . , 1] ∈ R n . Let J 1 = {1, 2 .
. . , p} and J 2 = {1, 2 . . . , q} be the index sets of constraints of player 1 and 2, respectively. Given a mixed strategy y of player 2, the objective of player 1 is to find a mixed strategy x that solves the following linear programming problem,

                           max x x T Ay s.t. Bx ≤ b, 1 T m x = 1, x ≥ 0. (1)
Similarly, for a given strategy x of player 1, the aim of player 2 is to find a strategy y that solves the following linear programming problem.

                           min y x T Ay s.t.
Dy ≥ d,

1 T n y = 1, y ≥ 0.
(2)

A strategy profile (x, y) is called a saddle point of this game if x and y are optimal solutions of ( 1) and (2), respectively. [START_REF] Charnes | Constrained games and linear programming[END_REF] proved that there is always a saddle point for such a game.

Definition 1 (Elliptically symmetric distribution). A d-dimensional random vector ξ follows an elliptically symmetric distribution Ellip d (µ, Σ, φ), i.e., ξ ∼ Ellip d (µ, Σ, φ) if its characteristic function is given by Ee iz T ξ = e iz T µ φ z T Σz , where z ∈ R d is the argument of the characteristic function, µ is the location parameter, Σ is the scale matrix, and φ is the characteristic generator function.

In this paper, we aim at solving the stochastic two-player zero-sum game introduced by [START_REF] Singh | A second-order cone programming formulation for two player zero-sum games with chance constraints[END_REF], in which the matrices B and D are treated as random variables. A stochastic two player zero sum game is characterized as follows.

• Let B ω and D ω denote random matrices that defines the constraints of player 1 and 2, respectively, and ω denotes some uncertainty parameter.

• We denote the k-th row vector of B ω as B ω k and the l-th row vector of

D ω as D ω l . Both B ω k and D ω l follow elliptical distributions, i.e., B ω k ∼ Ellip m µ 1 k , Σ 1 k , φ 1 k and D ω l ∼ Ellip n µ 2 l , Σ 2 l , φ 2 l .
The k-th element of b is denoted as b k , and the l-th element of d is denoted as d l .

• α 1 k ∈ [0, 1] denotes the confidence level for the k-th constraint of player 1, and α 2 l ∈ [0, 1] denotes the confidence level for the l-th constraint of player 2. Let α 1 = α 1 k k∈J1 , and α 2 = α 2 l l∈J2 , and

α = α 1 , α 2 .
Therefore, the stochastic two-player zero-sum game can be written as:

                           max x x T Ay s.t. P {B ω k x ≤ b k } ≥ α 1 k , ∀k ∈ J 1 1 T m x = 1 x ≥ 0, (3)
and

                           min y x T Ay s.t. P {D ω l y ≥ d l } ≥ α 2 l , ∀l ∈ J 2 1 T n y = 1 y ≥ 0. (4)
A mixed strategy (x * , y * ) is said to be a saddle point of the stochastic two-player zero-sum game if it simultaneously solves for both players' optimization problems, as described in equations ( 3) and ( 4). In other words, (x * , y * ) represents the optimal strategies for each player, such that neither player can unilaterally improve their outcome by deviating from these strategies.

120

By the second-order cone constraint reformulation [START_REF] Henrion | Structural properties of linear probabilistic constraints[END_REF][START_REF] Van De Panne | Minimum-cost cattle feed under probabilistic protein constraints[END_REF][START_REF] Kataoka | A stochastic programming model[END_REF], the feasible sets of (3) and ( 4) are reformulated as

S 1 α 1 = x ∈ R m | 1 T m x = 1, x ≥ 0, x T µ 1 k + Ψ -1 ξ 1 k (α 1 k )∥(Σ 1 k ) 1 2 x∥ ≤ b k , ∀k ∈ J 1 , (5)
and

S 2 α 2 = y ∈ R n | 1 T n y = 1, y ≥ 0, -y T µ 2 l + Ψ -1 ξ 2 l (α 2 l )∥(Σ 2 l ) 1 2 y∥ ≤ -d l , ∀l ∈ J 2 , (6)
respectively, where

ξ 1 k = B ω k x-x T µ 1 k (Σ 1 k ) 1 2 x , k ∈ J 1 , and ξ 2 l = -D ω l y+y T µ 2 l (Σ 2 l ) 1 2 y , l ∈ J 2 follow univariate standard elliptical distributions, i.e., ξ 1 k ∼ Ellip 0, 1, φ 1 k , and ξ 2 l ∼ Ellip 0, 1, φ 2 l . Ψ -1 ξ 1 k α 1 k and Ψ -1 ξ 2 l α 2 l
are the quantile functions of 1-dimensional distribution functions induced by characteristic functions φ 1 k and φ 2 l , respectively. We denote a stochastic two-player zero-sum game with a confidence level α as G(α).

Assumption 1.

125

• S 1 (α 1 ) is strictly feasible, i.e., there exists an x ∈ R m which is a feasible point of S 1 (α 1 ) and the inequality constraints of S 1 (α 1 ) are strictly satisfied by x.

• S 2 (α 2 ) is strictly feasible, i.e., there exists an y ∈ R n which is a feasible point of S 2 (α 2 ) and the inequality constraints of S 2 (α 2 ) are strictly satisfied by x.

Theorem 1 [START_REF] Singh | A second-order cone programming formulation for two player zero-sum games with chance constraints[END_REF], Theorem 3.5). Consider a stochastic two player zero sum game G(α), 130 where the matrices B w and D w defining the constraints of both the players, respectively, are random. Let

the row vectors B w k ∼ Ellip m µ 1 k , Σ 1 k , φ 1 k , k ∈ J 1 , and D w l ∼ Ellip n µ 2 l , Σ 2 l , φ 2 l , l ∈ J 2 . If all k ∈ J 1 and l ∈ J 2 , Σ 1
k and Σ 2 l are positive definite matrices. Then, there exists a saddle point equilibrium for the game G(α) for all α ∈ (0.5, 1] p × (0.5, 1] q .

Theorem 2 [START_REF] Singh | A second-order cone programming formulation for two player zero-sum games with chance constraints[END_REF], Theorem 3.7). Let Assumption 1 hold. (x * , y * ) is the saddle point of the stochastic two-player zero-sum game G(α) if and only if there exist

v 1 * , δ 1 k * k∈J1 , λ 1 * and v 2 * , δ 2 l * l∈J2 , λ 2 * such that y * , v 1 * , δ 1 k * k∈J1 , λ 1 * and x * , v 2 * , δ 2 l * l∈J2 , λ 2 *
are optimal solutions of the following primaldual pair of nonlinear optimization problems.

min y,v 1 ,(δ 1 k ) k∈J 1 ,λ 1 v 1 + k∈J1 λ 1 k b k s.t. (i)Ay -k∈J1 λ 1 k µ 1 k -k∈J1 Σ 1 k 1 2 δ 1 k ≤ v 1 1 m (ii) -y T µ 2 l + Ψ -1 ξ 2 l α 2 l Σ 2 l 1 2 y ≤ -d l , ∀l ∈ J 2 (iii) δ 1 k ≤ λ 1 k Ψ -1 ξ 1 k α 1 k , ∀k ∈ J 1 (iv)1 T y = 1 (v)y ≥ 0 (vi)λ 1 k ≥ 0, ∀k ∈ J 1                                                  (P) max x,v 2 ,(δ 2 l ) l∈J 2 ,λ 2 v 2 + l∈J2 λ 2 l d l s.t. (i)A T x -l∈J2 λ 2 l µ 2 l -l∈J2 Σ 2 l 1 2 δ 2 l ≥ v 2 1 n (ii)x T µ 1 k + Ψ -1 ξ 1 k α 1 k Σ 1 k 1 2 x ≤ b k , ∀k ∈ J 1 (iii) δ 2 l ≤ λ 2 l Ψ -1 ξ 2 l α 2 l , ∀l ∈ J 2 (iv)1 T x = 1 (v)x ≥ 0 (vi)λ 2 l ≥ 0, ∀l ∈ J 2                                                  (D) 135
We refer readers to [START_REF] Singh | A second-order cone programming formulation for two player zero-sum games with chance constraints[END_REF] for the proof and other details related to Theorems 1 and 2.

Additionally, we refer readers to Wu & Lisser (2022b) (Proposition 1) for the reformulation of a stochastic two-player zero-sum game G(α) to (P) and (D).

Neurodynamic optimization

The neurodynamic optimization method introduced by Wu & Lisser ( 2022b) is used to model the stochastic two-player zero-sum game by a ODE system. Because (P) and (D) are a primal-dual pair, it is only necessary to model (P). The optimal solution for (D) can be found in the dual solution of (P).

2.2.1. Reformulation of (P).

Let us denote the decision variables of (P) as

s = (y ∈ R n , v 1 ∈ R, δ ∈ R n * p , λ ∈ R p ), where δ = [δ 1 1 T , . . . , δ 1 p T ]
T , and λ = [λ 1 1 , . . . , λ 1 p ] T . Define ns = n + 1 + n * p + p as the number of decision variables, and s ∈ R ns . We reformulate the equality constraint 1 T n y = 1 in (P) as two inequalities, i.e., 1 T n y -1 ≤ 0 and 1 -1 T n y ≤ 0. The second-order cone constraints in (P) are not differentiable at the origin, so a smoothness technique must be used to transform them into smooth functions. This technique is given by the following expression

∥s∥ smooth = ∥s∥ 2 + ϵ 2 , ( 7
)
where ϵ is a small positive constant, typically chosen to be on the order of 10 -6 or 10 -8 , and is set to ϵ = 10 -6 in this paper. Then, the optimization problem (P) can be rewritten as

min s f (s) s.t. g(s) ≤ 0, (8)
where the objective function f : R ns → R is denoted as

f (s) = v 1 + k∈J1 λ 1 k b k , (9)
and the constraints are denoted as

g(s) =                  g 1 (s) g 2 (s) g 3 (s) g 41 (s) g 42 (s) g 5 (s) g 6 (s)                  =                  Ay -v1 m -k∈J1 Σ 1 k 1 2 δ k -k∈J1 λ k µ 1 k (-y T µ 2 l + Ψ -1 ξ 2 l α 2 l Σ 2 l 1 2 y smooth + d l ) l∈J2 (∥δ k ∥ smooth -Ψ -1 ξ 1 k α 1 k λ k ) k∈J1 1 T n y -1 -1 T n y + 1 -y -λ                  , ( 10
)
where g 1 : R n+1+n * p+p → R n , g 2 : R n+1+n * p+p → R q , g 3 : R n+1+n * p+p → R p , g 41 : R n+1+n * p+p → R, g 42 : R n+1+n * p+p → R, g 5 : R n+1+n * p+p → R n , and g 6 : R n+1+n * p+p → R q . Let define nu = n+q +p+1+1+n+q

Karush-Kuhn-Tucker conditions

The Karush-Kuhn-Tucker conditions (KKT conditions) for the optimization problem (8) are

∇f (s) + ∂g ∂s (s) T u = 0, g(s) ≤ 0, u T ≥ 0, u T g(s) = 0, (11)
where u ∈ R nu denotes the dual variables. ∇f (s) ∈ R ns denotes the gradient of the objective function at point s. ∂g ∂s (s) ∈ R nu×ns denotes the Jacobian matrix of the constraints at point s. The details of u ∈ R nu , ∇f (s) ∈ R ns , and ∂g ∂s (s) ∈ R nu×ns are given as follows

u =                  u 1 ∈ R n u 2 ∈ R q u 3 ∈ R p u 41 ∈ R u 42 ∈ R u 5 ∈ R n u 6 ∈ R q                  , ∇f (s) =         ∇f y (s) ∈ R n ∇f v (s) ∈ R ∇f δ (s) ∈ R p * n ∇f λ (s) ∈ R p         , ( 12
) ∂g ∂s (s) = ∂g ∂y (s) ∈ R nu×n ∂g ∂v (s) ∈ R nu×1 ∂g ∂δ (s) ∈ R nu×n * p ∂g ∂λ (s) ∈ R nu×p , =                  ∂g1 ∂y (s) ∈ R n×n ∂g1 ∂v (s) ∈ R n×1 ∂g1 ∂δ (s) ∈ R n×n * p ∂g1 ∂λ (s) ∈ R n×p ∂g2 ∂y (s) ∈ R q×n ∂g2 ∂v (s) ∈ R q×1 ∂g2 ∂δ (s) ∈ R q×n * p ∂g2 ∂λ (s) ∈ R q×p ∂g3 ∂y (s) ∈ R p×n ∂g3 ∂v (s) ∈ R p×1 ∂g3 ∂δ (s) ∈ R p×n * p ∂g3 ∂λ (s) ∈ R p×p ∂g41 ∂y (s) ∈ R 1×n ∂g41 ∂v (s) ∈ R 1×1 ∂g41 ∂δ (s) ∈ R 1×n * p ∂g41 ∂λ (s) ∈ R 1×p ∂g42 ∂y (s) ∈ R 1×n ∂g42 ∂v (s) ∈ R 1×1 ∂g42 ∂δ (s) ∈ R 1×n * p ∂g42 ∂λ (s) ∈ R 1×p ∂g5 ∂y (s) ∈ R n×n ∂g5 ∂v (s) ∈ R n×1 ∂g5 ∂δ (s) ∈ R n×n * p ∂g5 ∂λ (s) ∈ R n×p ∂g6 ∂y (s) ∈ R q×n ∂g6 ∂v (s) ∈ R q×1 ∂g6 ∂δ (s) ∈ R q×n * p ∂g6 ∂λ (s) ∈ R q×p                  . ( 13
T (t), u T (t)) T = (y T (t), v T (t), δ T (t), λ T (t), u T (t)) T ,
and r : R → R nr , where nr = ns + nu. The following ODE system models the KKT conditions ( 11)

dr dt = Φ(r) =            dy dt dv dt dδ dt dλ dt du dt            =                -∇f y (s) + ∂g ∂y (s) T (u + g (s)) + -∇f v (s) + ∂g ∂v (s) T (u + g (s)) + -∇f δ (s) + ∂g ∂δ (s) T (u + g (s)) + -∇f λ (s) + ∂g ∂λ (s) T (u + g (s)) + (u + g (s)) + -u                , (14)
where (u + g (s))

+ = max{0, u + g (s)}.
Definition 2 (State solution). Consider an ODE system dr dt = Φ(r), where Φ(r) : R nr → R nr . Given an initial point r 0 ∈ R nr , a vector value function r(t) : R → R nr is called a state solution, if it satisfies the ODE system dr dt = Φ(r) and the initial condition r(0) = r 0 .

Definition 3 (Equilibrium point). Consider an ODE system dr dt = Φ(r). A point r * is called an equilibrium point if it satisfies Φ(r * ) = 0.

Theorem 3 [START_REF] Wu | A dynamical neural network approach for solving stochastic two-player zerosum games[END_REF]). Let r * = (s * , u * ) T . The point r * is an equilibrium point of the ODE system (14) if and only if (s * , u * ) satisfy the KKT conditions (11).

Lemma 4 [START_REF] Wu | A dynamical neural network approach for solving stochastic two-player zerosum games[END_REF]). The equilibrium point r * of the proposed ODE system (14) is unique.

Theorem 5 [START_REF] Wu | A dynamical neural network approach for solving stochastic two-player zerosum games[END_REF]). The equilibrium point r * of the proposed ODE system (14) is globally asymptotically stable, i.e., given any initial point r 0 , the state solution has lim t→∞ r(t) = r * .

Remark 1. The primary objective of the ODE system presented in (14) is to address the KKT conditions given in (11). As stated in Theorem 3, the equilibrium point in the ODE system coincides with to the KKT conditions outlined in (11). As a result, by determining the equilibrium point of the ODE system, we can effectively solve the KKT conditions and subsequently identify the saddle point of the stochastic two-player zero-sum game.

Lemma 4 establishes the uniqueness of the equilibrium point in the ODE system, which is consistent with the notion that the stochastic two-player zero-sum game has a unique saddle point. Furthermore, Theorem 5 ensures that the state solution of the ODE system converges to the equilibrium point of the system as t approaches infinity, implying that it converges to the saddle point of the stochastic two-player zero-sum game.

This means that we can solve the saddle point of the corresponding game by simply solving the state solution of the ODE system.

Numerical integration method

In this subsection, we briefly describe how existing methods solve for the state solution of an ODE system, which will be compared with our proposed method in the experimental section. Since the considered ODE system ( 14) is nonlinear, it cannot be solved analytically.

The existing methods for solving the state solution of an ODE system use numerical integration techniques to approximate the solution by discretizing the domain [START_REF] Butcher | Numerical methods for ordinary differential equations[END_REF]. These methods, such as the Runge-Kutta method, choose multiple time points in the domain, known as collocation points, and find a solution that satisfies the ODE system at these points. However, these methods can be inefficient when only the final state is of interest because the final state is only obtained after calculating all the intermediate collocation points.

Numerical integration methods are further divided into two categories: explicit and implicit methods.

Explicit methods, such as RK45, RK23, and DOP853, determine the state at a later time based on the current state [START_REF] Dormand | A family of embedded runge-kutta formulae[END_REF][START_REF] Bogacki | A 3 (2) pair of runge-kutta formulas[END_REF][START_REF] Hairer | Solving ordinary differential equations[END_REF]. Implicit methods, such as Radau and BDF, find the solution by solving equations involving the current and later states [START_REF] Wanner | Solving ordinary differential equations II volume 375[END_REF][START_REF] Shampine | The matlab ode suite[END_REF]. Additionally, LSODA can switch automatically between stiff and nonstiff methods [START_REF] Petzold | Automatic selection of methods for solving stiff and nonstiff systems of ordinary differential equations[END_REF]. Scipy provides software implementations of these methods to make them easier to use [START_REF] Virtanen | & SciPy 1.0 Contributors[END_REF].

Method

In Section 3.1, we describe the problem setup and use an initial value problem to approximate the stochastic two-player zero-sum game. Section 3.2 presents our proposed approach, which uses a neural network to solve the initial value problem and subsequently predict the saddle point of the game. The loss function for training the neural network is discussed in Section 3.3, and an evaluation metric is given in Section 3.4 to assess the model's effectiveness at predicting saddle points. Finally, Section 3.5 provides a complete pipeline for solving the game problem using our proposed neural network.

Problem setup

This subsection provides a brief summary of how to reformulate a stochastic two-player zero-sum game as an initial value problem based on Sections 2.1 and 2.2. As presented in Figure 1, the reformulation procedure is as follows:

1. Reformulate the original game problem G(α), i.e., (3)-( 4), as a primal-dual pair of nonlinear optimization problems (P) and (D). 7) to derive the KKT conditions (11) of the primal problem (P).

Apply the smoothness technique (

3. Represent the KKT conditions using the ODE system (14), where the equilibrium point is unique and globally asymptotically stable.

4. The desired saddle point (x * , y * ) of the game is included in the equilibrium point r * of the ODE system.

In order to obtain the equilibrium point, we must construct an initial value problem (IVP), which consists of (a) the ODE system ( 14

Optimal solution of KKT conditions

Equilibrium of ODE system

Figure 1: The upper part illustrates the workflow of modeling a stochastic two-player zero-sum game (STPZSG) using an ODE system. The STPZSG is represented by equations ( 3) and ( 4), and the ODE system is represented by equation ( 14). The nonlinear optimization problems (NOPs) in equations (P) and (D) and the KKT conditions in equation ( 11) are also shown.

The lower part demonstrates that the saddle point of the STPZSG is incorporated in the optimal solution of the NOPs and KKT conditions, as well as the equilibrium point of the ODE system.

Fully connected network

Fully connected network

Predicted state solution

When

Predicted saddle point When given any initial point r 0 , the corresponding state solution r(t) will converge to the same equilibrium r * as time t goes to infinity. Therefore, we set the initial point of the ODE system to all zeros, i.e.,

r 0 = 0. ( 15
)
The final state r(T ) approximates the equilibrium point r * , i.e.,

r(T ) ≈ r * . ( 16
)
According to Theorem 5, the larger the provided time range [0, T ], the closer the final state is to r * . Our 205 goal is to solve for the final state r(T ) of the state solution r(t), t ∈ [0, T ], given an adequate size of the time range.

Neural network model

In this subsection, we present a model based on a fully connected neural network that is designed to predict the state solution r(t) and its final state r(T ). The proposed model is given as follows

r(t; w) = 1 -e -t N(t; w), t ∈ [0, T ], ( 17
)
where N(t; w) is a fully connected neural network with trainable parameters w.

(1 -e -t ) is an auxiliary function to ensure the satisfaction of the initial point, i.e., r(t = 0; w) = 0.

210

As shown in Figure 2(Left), the model r(t; w) itself is used as a predicted state solution to the IVP, i.e.,

r(t; w) ≈ r(t), t ∈ [0, T ]. ( 18
)
By choosing the time t = T , we have r(t = T ; w) ≈ r(T ). ( 19)

Combining ( 16) with ( 19), we have

r(t = T ; w) ≈ r * . ( 20
)
As shown in Figure 2(Right), r(t = T ; w) contains the predicted saddle point of the considered game problem.

The effectiveness of the proposed neural network in solving stochastic two-player zero-sum game problems is supported by two theorems: 1) Theorem 5, which ensures the convergence of the state solutions of the ODE system to the saddle point of the game.

2) The universal approximation theorem for neural networks [START_REF] Hornik | Approximation capabilities of multilayer feedforward networks[END_REF], which states that there always exists a neural network as a state solution for solving the ODE system.

Loss function

We define the loss function of the proposed model as

L(t, w) = ∂ r(t; w) ∂t -Φ(r(t; w)) , (21)
where Φ(•) refers to the ODE system of ( 14), which corresponds to the game problem being studied. ∂ r(t;w) ∂t is the derivative of the model output r(t; w) with respect to the input time t, which can be computed analytically by automatic differentiation tools, e.g., PyTorch or JAX (Paszke et al., 2019b;[START_REF] Bradbury | JAX: composable transformations of Python+NumPy programs[END_REF].

We denote the objective function of the proposed model as However, calculating E(w) is computationally difficult because of its integral part. In practice, we train the model by minimizing the following batch loss

E(w) = T 0 L(t, w)dt. ( 22
L(T, w) = 1 |T| t∈T L(t, w), (23)
where T contains randomly sampled time points from the interval [0, T ], and |T| denotes the size of the set.

Following the approach used in previous studies [START_REF] Sirignano | Dgm: A deep learning algorithm for solving partial differential equations[END_REF], we sample the time points in T with a uniform distribution over the interval [0, T ]. These time points are sampled in an independent and identically distributed manner. The choice of probability distribution for sampling time points can potentially affect the training performance, which warrants further investigation in future studies.

Evaluation metric

To evaluate the performance of our model in predicting the final state r(t = T ; w), we propose the following evaluation metric:

ϵ(r(T ; w)) =      f (ŝ) if ∥g(ŝ)∥ ∞ ≤ τ, +∞ otherwise , (24)
where ϵ(r(T ; w)) is called the epsilon value of r(T ; w). ŝ is retrieved from the predicted final state, i.e., r(T ; w) = [ŝ, û]. f (•) and g(•) are defined in ( 9) and ( 10), respectively. τ is a positive constant that represents a threshold parameter. When ∥g(ŝ)∥ ∞ ≤ τ , the epsilon value is set to the objective value f (ŝ); otherwise, it is set to +∞.

Pipeline

Algorithm 1: Enhancing neurodynamic optimization with deep learning for solving stochastic twoplayer zero-sum games Input : A stochastic two player zero sum game G(α), A time range [0, T ] Output: Predicted saddle point (x, ŷ) 1 Function Main: 2 Derive the ODE system Φ(•) corresponding to the game G(α).

3

Initialize a neural network model r(t; w).

4

Initialize ϵ best = ϵ(r(T ; w)). Algorithm 1 integrates all the methods introduced in this section. This algorithm is an optimization procedure that minimizes the objective function E(w) while considering the proposed evaluation metric (24).

First, the algorithm reformulates the game problem as an IVP for a given time range [0, T ]. Next, it initializes a neural network to solve this problem and trains the network to improve its accuracy. In training, each iteration samples a batch of t uniformly from the time range [0, T ] as a training dataset. The neural network then performs gradient descent on the batch loss L(T, w) to update its parameters w. Finally, x and ŷ, which are obtained from the predicted final state, need to be transformed into probability distributions using the softmax function.

ϵ best represents the lowest epsilon value that the model has achieved so far. ϵ temp represents the epsilon value that the model obtained at the most recent iteration. After each update of the model's parameters, the algorithm calculates ϵ temp and compares it with ϵ best . If ϵ temp is lower than ϵ best , it indicates that the model has found a better prediction after this iteration of parameter updates. In this case, the algorithm sets ϵ best = ϵ temp and saves the current model. This idea is similar to the concept of early-stopping in deep learning, but here, we consider the epsilon value rather than the loss value.

Algorithm 1 requires specifying a time range as a hyperparameter. According to Theorem 5, the larger the chosen time range, the closer the corresponding final state will be to the optimal solution. In practice, however, a large time range can lead to training difficulties because the neural network has to approximate the state solution over a larger time range. Therefore, the choice of the time range involves a trade-off. A larger time range offers a higher maximum achievable accuracy but requires more training resources, whereas a smaller time range allows for faster convergence with fewer training iterations but results in a weaker achievable accuracy. One should determine the most appropriate time range that strikes a balance between accuracy and available training resources for their specific application.

In addition, the design of the neural network structure should be guided by two factors: 1) the chosen time range and 2) the size of the stochastic two-player zero-sum game to be solved. The size of the stochastic two-player zero-sum game is determined by the sizes of the action sets for both players and the number of constraints in their respective sets. The complexity of the neural network structure, which includes the number of neurons and the number of hidden layers, should be tailored to accommodate these two factors.

Numerical results

In this section, we conduct experiments to evaluate the performance of our proposed neural network model for solving the stochastic two-person zero-sum game. We first present the problem instance and describe the training process of our proposed model in Sections 4.1 and 4.2. Next, we compare our proposed method with the numerical integration method in Section 4.3, and discuss the advantages and limitations of our method in Section 4.4.

We use the Google Colab Pro+ platform to perform our experiments. The neural network model is implemented using Pytorch 1.12.1 with CUDA 11.2 (Paszke et al., 2019b), and the ODE system is modelled with JAX [START_REF] Bradbury | JAX: composable transformations of Python+NumPy programs[END_REF]. Our proposed method is compared with six numerical integration solvers, namely RK45, RK23, DOP853, Radau, BDF, and LSODA [START_REF] Dormand | A family of embedded runge-kutta formulae[END_REF][START_REF] Bogacki | A 3 (2) pair of runge-kutta formulas[END_REF][START_REF] Hairer | Solving ordinary differential equations[END_REF][START_REF] Wanner | Solving ordinary differential equations II volume 375[END_REF][START_REF] Shampine | The matlab ode suite[END_REF][START_REF] Petzold | Automatic selection of methods for solving stiff and nonstiff systems of ordinary differential equations[END_REF], which can be accessed using Scipy [START_REF] Virtanen | & SciPy 1.0 Contributors[END_REF].

Problem instance 275

We consider a specific instance of the stochastic two-player zero-sum game as shown in (3) and (4). Both players 1 and 2 have four pure strategies and one probabilistic constraint. B ω and D ω are 4-dimensional random vectors following normal distributions, i.e., B ω ∼ N µ 1 , Σ 1 , D ω ∼ N µ 2 , Σ 2 , where µ 1 and µ 2 are the means and Σ 1 and Σ 2 are the variances of the normal distributions. The problem data (A, Σ 1 , Σ 2 , µ 1 , µ 2 , b, d, α 1 , α 2 ) are given as follows:

A =        
8.67 0.17 6.04 3.81 7.09 0.37 9.83 6.89 4.75 7.32 4.06 6.52 2.20 2.22 4.76 3.10

        , Σ 1 =        
1.10 0.00 0.00 0.00 0.00 1.40 0.00 0.00 0.00 0.00 1.38 0.00 0.00 0.00 0.00 1.43

        , (25)
Σ 2 =        
1.70 0.00 0.00 0.00 0.00 1.35 0.00 0.00 0.00 0.00 1.69 0.00 0.00 0.00 0.00 1.28

        , µ 1 =         4.53 0.06 2.53 3.66         , µ 2 =         1.74 3.83 7.37 8.77         , (26

Model training

We construct the proposed model to solve the game instance, r(t; w) = (1 -e -t ) N(t; w), where N(t; w) is a fully connected neural network with a hidden layer of 200 neurons and the activation function is Tanh.

The training hyperparameters are as follows: the maximum number of iterations is 100, 000, the optimizer is 280

Adam [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] with a learning rate of 0.001 and batch size of 128. We apply the ODE system (14) to model the game instance given in Section 4.1. The neural network model is trained using Algorithm 1, with the time range set to [0, 10] and the threshold τ in equation ( 24) Table 2(left) lists the predicted saddle points given by our proposed method at the 1st, 1000th, 10,000th, 30,000th, 50,000th, and 100,000th training iterations. Table 2(right) lists the predicted saddle points given by the RK45 method, a typical numerical integration solver, at the 1st, 1000th, 10,000th, 30,000th, 50,000th, and 100,000th collocation points.

Our method

Numerical integration methods We evaluate the performance of our model using two metrics: the objective value (OBJ, lower is better)

Iteration OBJ ↓ Collocation point RK45 RK23 DOP853 Radau BDF LSODA OBJ ↓ OBJ ↓ OBJ ↓ OBJ ↓ OBJ ↓ OBJ ↓ 1 -1.
and the maximum violation of constraints (MVC, lower is better). OBJ is defined as f (s), where f (•) is given by equation ( 9) and here s represents the prediction from our method or a numerical integration method.

MVC is defined as |g(s)| ∞ , where g(•) is given by equation (10).

Figure 4, Table 3, and Table 4 compare our method with numerical integration methods using the OBJ and MVC metrics. Figure 4 shows the results in a continuous manner, while Tables 3 and4 show the results in a discrete manner at selected iterations or collocation points. Figure 4 only compares our method with the RK45 method, while Tables 3 and4 compare six numerical integration methods.

The comparison results yield the following insights:

• Our method outperforms all numerical integration methods in terms of evaluation metrics OBJ and MVC, as shown in Tables 3 and4. Ultimately, our method achieves an OBJ of 5.48 and an MVC of 0.0827. Among the numerical integration methods, BDF attains the best results, with an OBJ of 5.54 and an MVC of 0.1235.

• Our method converges faster than the numerical integration methods. Figure 4 shows that our method obtains an almost exact solution at the 30,000th iteration, with an OBJ of 5.48 and an MVC of 0.0827.

In comparison, the results of the numerical integration methods at the 30,000th collocation point show OBJ and MVC values of 0.317 and 0.8481, respectively, indicating a significant deviation from accuracy.

• Our proposed method can identify a solution that nearly satisfies the constraints in the early stages of solving. Table 4 shows that at the 1000th iteration, our method finds a solution with an MVC of less than 0.1. In contrast, the numerical integration method fails to find solutions with MVC less than 0.1 even after 100,000 collocation points.

• Figure 3 shows a significant decrease in the MSE loss value during the first 20,000 iterations, with a subsequent levelling off. Figure 4 shows that OBJ and MVC still improve significantly after 20,000 iterations, although the decrease in loss value has become less pronounced. This suggests that the MSE loss of the model is not entirely decisive for the OBJ and MVC indicators, although they are correlated.

• The performance of the six numerical integration methods does not differ significantly from each other.

Discussion

In Sections 4.2 and 4.3, we present experimental results that demonstrate the superiority of our proposed deep learning approach. The observed advantages are mainly due to three key aspects:

(i) During each iteration, the neural network can directly predict the saddle point of the game without having to calculate all the intermediate states in the ODE system. This is in contrast to traditional numerical integration methods, which require all intermediate states to be calculated sequentially in order to make a prediction.

(ii) In the algorithm 1, we incorporate several important designs that exploit the structure of the problem.

These include the use of softmax functions to ensure the feasibility of the mixed strategy and the implementation of a mechanism that retains the best results from each iteration.

(iii) Our approach transforms the stochastic game problem into a neural network training problem, eliminating the dependence on numerical integration solvers. This allows us to solve the problem using only deep learning infrastructure, taking advantage of the latest advances such as GPU parallel computing.

While our method demonstrates significant performance advantages, we must acknowledge its limitations:

• Our approach requires careful hyperparameter tuning of the neural network to achieve optimal performance. In contrast, numerical integration methods are more straightforward and do not require numerous algorithmic settings.

• Our approach relies on training neural networks, which can be unstable. In some cases this instability may lead to rapid convergence to a satisfactory solution, while in other cases it may take many iterations to find an appropriate solution. Possible factors contributing to this inconsistency include the initial state of the neural network and the learning rate setting. The underlying cause is due to the black-box nature of neural networks and deserves further investigation.

Conclusion

In this paper, we presented a novel approach to find the saddle point of stochastic two-player zero-sum games by combining neurodynamic optimization and neural networks. This results in faster convergence and superior solutions compared to traditional methods. Our work establishes a link between stochastic games and deep learning, which can benefit the stochastic game theory given the rapid growth of the deep learning community.

However, there are limitations to our work, such as the stability and robustness of the method. Future

  s(•), y(•), v(•), δ(•), λ(•) and u(•) be time dependent functions, i.e., s : R → R ns , y : R → R n , v : R → R, δ : R → R n * p , and λ : R → R p . Denote r(t) = (s

  ), (b) an initial point r 0 , and (c) a time range [0, T ]. The state solution r(t) of the IVP satisfies the ODE system on the domain t ∈ [0, T ] and the initial condition r(0) = r 0 . By Theorem 5,

Figure 2 :

 2 Figure2: The left part of the figure shows the neural network model, r(t; w) = 1 -e -t N(t; w), t ∈ [0, T ], acting as a predicted state solution to the IVP. The right part of the figure illustrates that when the final time is selected (i.e., t = T ), r(t; w) contains the predicted saddle point.

)

  The objective function E(w) is an integral of the loss function over the pre-given time range [0, T ]. The loss value L(t, w) denotes the error of the model at the time t, and the objective function E(w) denotes the overall error of the proposed model over the time range [0, T ].

5

  while iter ≤ Max iteration do 6 T ∼ U (0, T ): Uniformly sample a batch of t from the interval [0, T ]. 7 Forward propagation: Compute the batch loss L(T, w). 8 Backward propagation: Update w by ∇ w L(T, w). 9 Compute the epsilon value: ϵ temp = ϵ(r(T ; w)). 10 if ϵ temp < ϵ best then 11 ϵ best = ϵ temp 12 Save the model with parameters w 13 end 14 Extract x and ŷ from r(T ; w).

  ) b = 7.03, and d = 4.53. The confidence levels of both players are set to α 1 = α 2 = 80%.

Figure 3 :

 3 Figure 3: The mean square error (M.S.E) loss plotted against the number of iterations.

Figure 4 :

 4 Figure 4: The objective value plotted against the number of iterations (left figure). The maximum violation of constraints plotted against the number of iterations (right figure).

  research should focus on: (a) using state-of-the-art machine learning and deep learning techniques to improve our method, (b) exploring advanced neurodynamic optimization approaches for problem modelling, and (c) extending our method to other stochastic game problems.BibliographyAnitescu, C., Atroshchenko, E., Alajlan, N., & Rabczuk, T. (2019). Artificial neural network methods for the solution of second order boundary value problems. Computers,Materials and Continua, 59 , 345-359. Baydin, A. G., Pearlmutter, B. A., Radul, A. A., & Siskind, J. M. (2018). Automatic differentiation in machine learning: a survey. Journal of machine learning research, 18 .

  Table 1 presents the notations list of this paper.

	Notation	Description
	M = {1, . . . , m}	Action set of player 1; m denotes the number of pure strategies.
	N = {1, . . . , n}	Action set of player 2; n denotes the number of pure strategies.
	x ∈ R m	Mixed strategy of player 1.
	y ∈ R n	Mixed strategy of player 2.
	J 1	Index set of constraints for player 1.
	J 2

Table 1 :

 1 Summary of mathematical notations

		elliptically symmetric distribution.
	s ∈ R ns	Decision variables; ns denotes the number of decision variables.
	u ∈ R nu	Dual variables; nu denotes the number of dual variables.
	Φ(•)	ODE system.
	[0, T ]	Time range.
	r : [0, T ] → R nr	State solution of the ODE system.
	r : [0, T ] → R nr	Approximate state solution.
	∥•∥	Euclidean norm.
	∥•∥ ∞	Infinity norm.

Table 3 :

 3 Comparison of the objective (OBJ) value produced by our method and various numerical integration methods, including RK45, RK23, DOP853, Radau, BDF, and LSODA. The table shows the OBJ values for different iterations and collocation points.

Table 4 :

 4 Comparison of the maximum violation of constraints (MVC) produced by our method and various numerical integration methods, including RK45, RK23, DOP853, Radau, BDF, and LSODA. The table shows the MVC values for different iterations and collocation points.

as the number of constraints, and g : R ns → R nu . set to 0.1. We use the MSE loss as a performance metric to evaluate the model's ability to solve the derived ODE system. Figure 3 shows the evolution of the MSE loss throughout the training process. The MSE loss decreases from an initial value of 4,557 to 53 after 100,000 iterations, indicating the transition of the neural network from a naive approximation to an accurate state solution of the ODE system. It is important to note that the reported loss values do not distinguish between training and testing losses, as the model trains on a completely new randomly generated time set T in each iteration. This training strategy is consistent with a range of deep learning methods for solving differential equations [START_REF] Raissi | Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations[END_REF][START_REF] Sirignano | Dgm: A deep learning algorithm for solving partial differential equations[END_REF].

Comparison with numerical integration methods

In this subsection, we compare the performance of our proposed method with the numerical integration solvers on the stochastic two-player zero-sum game instance. We first show the predicted saddle points given by the two types of methods and then evaluate the performance of these saddle points on two evaluation metrics.

The setup of the numerical solver involves dividing the time range [0, 10] into 100, 000 collocation points, which corresponds to the number of training iterations of our method.