text:
\begin{document}

\title{Nash Equilibrium Problems of Polynomials}

\author[Jiawang Nie]{Jiawang~Nie}
\address{Jiawang Nie, Xindong Tang, Department of Mathematics,
University of California San Diego,
9500 Gilman Drive, La Jolla, CA, USA, 92093.}
\email{njw@math.ucsd.edu, xit039@ucsd.edu}

\author[Xindong Tang]{Xindong~Tang}

\subjclass[2010]{90C33, 91A10, 90C22, 65K05}

\date{}

\keywords{Nash Equilibrium, Polynomial Optimization, Moment-SOS relaxation,
Lagrange multiplier expression, Tight relaxation}

\begin{abstract}
This paper studies Nash equilibrium problems that are given by polynomial functions.
We formulate efficient polynomial optimization problems
for computing Nash equilibria. The Lasserre type Moment-SOS relaxations
are used to solve them. Under generic assumptions,
the method can find a Nash equilibrium if there is one.
Moreover, it can find all Nash equilibria
if there are finitely many ones of them.
The method can also detect nonexistence if there is no Nash equilibrium.
\end{abstract}

\maketitle

\section{Introduction}

The Nash Equilibrium Problem (NEP) is a kind of games for finding strategies
for a group of players such that each player's objective is optimized,
for given other players' strategies.
Suppose there are $N$ players and the $i$th player's strategy is the variable
$x_{i} \in \re^{n_{i}}$ (the $n_i$-dimensional real Euclidean space).
We denote that
\[
x_i\, := \, (x_{i,1}, \ldots, x_{i,n_i}), \quad
\allx \, := \, (x_1, \ldots, x_N).
\]
The total dimension of all players' strategies is
\[
n : \, = \, n_1 + \cdots + n_N.
\]
When the $i$th player's strategy $x_i$ is being optimized,
we use $x_{-i}$ to denote the subvector of all players' strategies except $x_i$, i.e.,
\[
x_{-i} \, := \, (x_1, \ldots, x_{i-1}, x_{i+1}, \ldots, x_N),
\]
and write $\allx = (x_i, x_{-i})$ accordingly.
When the writing $x_{-i}$ appears,
 the $i$th player's strategy is being considered for optimization,
while the vector of all other players' strategy is fixed to be $\xmi$.
In a NEP, the $i$th player's best strategy $x_i$
is a minimizer for the optimization problem
\be
\label{eq:subprob}

\mbf{F}_i(x_{-i}):
\left\{ \begin{array}{cl}
\min\limits_{\xpi\in \Rni}  &  \fpi(\xpi,\xmi) \\
\st & g_{i,j}(\xpi) = 0 \, (j \in \mc{E}_i), \\
    & g_{i,j}(\xpi) \geq 0 \, (j \in \mc{I}_i),
\end{array} \right.
\ee
for given $\xmi$ of other players' strategies.
In the above, $f_i$ is the $i$th player's objective function,
and $g_{i,j}$ are constraining functions in $\xpi$.
The $\mc{E}_i$ and $\mc{I}_i$ are disjoint labeling sets
of finite cardinalities (possibly empty).
The feasible set of the optimization $\mbf{F}_i(x_{-i})$ in \reff{eq:subprob} is
\be \label{eq:feaset}
X_i \, := \,
\{x_i \in \Rni: \, g_{i,j}(\xpi) = 0\, (j \in \mc{E}_i), \,
g_{i,j}(\xpi) \geq 0 \, (j \in \mc{I}_i)  \} .
\ee
It is called the feasible strategy set for the $i$th player.
For NEPs, each set $X_i$ does not depend on $x_{-i}$.
This is different from generalized Nash Equilibrium problems (GNEPs),
where each player's feasible set depends on other players' strategies.
The entire strategy vector $\allx$ is a feasible point if
\[
 \allx   = (x_1, \ldots, x_N) \, \in \,
X_1 \times \cdots \times X_n,
\]
that is, each $x_i \in X_i$.
The NEP can be formulated as
\be
\label{eq:NEP}

\mbox{find} \quad \allx^*  \in\nfR^{n}
\quad \mbox{such that each}
\,\,\xpi^* \,\, \mbox{is a minimizer of} \, \mbf{F}_i(x_{-i}^*),
\ee
where $\allx^* = (x_1^*, \ldots, x_N^*)$.
A solution of \reff{eq:NEP} is called a {\it Nash Equilibrium} (NE).
When all the defining functions $f_i$ and $g_{i,j}$ are polynomials in $\allx$,
the NEP is then called a Nash Equilibrium Problem of Polynomials (NEPP).
The following is an example.

\begin{example}
\label{ep:isolatedep}  \rm
Consider the $2$-player NEP with the individual optimization
\[
\baray{l}
\mbox{1st player:}  \left\{
  \begin{array}{rl}
    \min\limits_{x_1 \in \re^2} &
     x_{1,1}(x_{1,1}+x_{2,1}+4x_{2,2})+2x_{1,2}^2, \\
 \st & 1-(x_{1,1})^2-(x_{1,2})^2 \ge 0,
  \end{array}\right.
\\
\mbox{2nd player:}   \left\{
  \begin{array}{rl}
    \min\limits_{x_2 \in \re^2} & x_{2,1}(x_{1,1}+2x_{1,2}+x_{2,1})+x_{2,2}(2x_{1,1}+x_{1,2}+x_{2,2}), \\
    \st & 1 -(x_{2,1})^2 -(x_{2,2})^2 \ge 0.
  \end{array}\right.
\earay
\]
This NEP has only $3$ NEs (see Section~\ref{ssc:convex}), which are
\[
\begin{array}{rll}
\mbox{1st NE:} & x_1^* = (0, 0), & x_2^* = (0, 0); \\
\mbox{2nd NE:} & x_1^* = (1, 0), & x_2^* = \frac{1}{\sqrt{5}} (-1, -2 ) ;\\
\mbox{3rd NE:} & x_1^* = (-1, 0), & x_2^* = \frac{1}{\sqrt{5}}(1, 2 ) .
\end{array}
\]
It is interesting to note that each player's objective is strictly convex
with respect to its strategy, because their Hessian's
with respect to their own strategies are positive definite.

\end{example}

When each feasible set $X_i$ is a finite set,
the NEP is called a {\it finite game}. For finite games,
Nash Equilibria typically do not exist.
People are also interested in mixed strategies,
which are probability distributions on the strategy set.

Mixed strategy solutions always exist for finite games \cite{Nash1951}.
For the case that $f_1+\dots+f_N = 0$,
the NEP \reff{eq:NEP} is called a {\it zero-sum game}.
The two-player zero-sum game is equivalent to a saddle point problem \cite{Nie2018saddle}.
NEPs have broad applications in Economics modelling

\cite{Aumann1995,maskin1999,Rubinstein1982,Shaked1982,Wildasin1988}.
Applications outside Economics can be found in
\cite{Bell2000,Bjrnerstedt1994,Garfinkel1994,MacKenzie2001,Schofield2002}.
We refer to \cite{Aubin13,BasOls98,Datta03,Myer13,OsbRub94,Rosen65}
for more general work on NEPs.

\iffalse

Moreover, the Generative Adversarial Network (GAN) \cite{Goodfellow2014}
is a useful class of machine learning frameworks,
which are based on solving game problems.
Training GANs is connected to NEPs
\cite{Goodfellow2014,Mao2017,Radford2015}.

\fi

It is generally hard to solve NEPs \cite{DasGoldPap09,Savani2006}.
For solving finite games or finding their mixed strategy solutions, we refer to the work
\cite{Ahmadi2017,Datta10,Kontogiannis2006,Lemke1964}.
For two-player zero sum games, the NEPs are equivalent to saddle point problems
and there are optimization methods for solving them
\cite{CLOu14,Nie2018saddle}.
More work for solving NEPs can be found in

\cite{GurPang09,Herings2010,krawczyk2000,LarLas12,LiBas87,LipMar04,Ratliff2013,UryRub94}.
Moreover, NEPs are special cases of
generalized Nash Equilibrium problems. We refer to
\cite{DreFacKanSag11,Facchinei2007,FacFisPic09,FacKan10OR,
FisMarSch14,nabetani2011parametrized,Pang2005}
for recent work on GNEPs. In the current state of the art,
it is mostly an open question to solve general NEPs efficiently,
especially those whose individual optimization problems are nonconvex.

\medskip \noindent
{\bf Contributions}
This paper focuses on Nash equilibrium problems that are given by polynomial functions,
i.e., the objectives $f_i$ and constraining functions $g_{i,j}$
are polynomial functions in their variables.
We formulate efficient polynomial optimization
for computing one or more Nash equilibria.
The Lasserre type Moment-SOS hierarchy of semidefinite relaxations
are used to solve the polynomial optimization.
Under generic assumptions, we prove the method can find
a Nash equilibrium if there exists one.
Moreover, we can also find all Nash equilibria
if there are finitely many ones of them.
The method can also detect nonexistence if there is no Nash equilibrium.
For NEPs that are given by generic polynomials,
we further show that there are finitely many complex KKT points.
Therefore, there are at most finitely many NEs for generic NEPPs.
This implies that our method can solve general NEPPs successfully.

The paper is organized as follows.
Some preliminaries about polynomial optimization
are given in Section~\ref{sc:pre}.
We give efficient polynomial optimization formulations
in Section~\ref{sc:finding}.
We show how to solve polynomial optimization problems for their optimizers
by the Moment-SOS hierarchy in Section~\ref{sc:solvepop}.
Numerical experiments and applications are given in Section~\ref{sc:num}.
Conclusions and discussions are proposed in Section~\ref{sc:conc}.
In the appendix, we prove that NEPs have finitely many complex KKT points,
when the functions are given by generic polynomials.

\section{Preliminaries}
\label{sc:pre}

\noindent {\bf Notation}
The symbol $\mathbb N$ (resp., $\mathbb R$, $\mathbb C$) stands for the set of
nonnegative integers (resp., real numbers, complex numbers).

For a positive integer $k$, denote the set $[k] := \{1, \ldots, k\}$.
For a real number $t$, $\lceil t \rceil$ (resp., $\lfloor t \rfloor$)
denotes the smallest integer not smaller than $t$
(resp., the biggest integer not bigger than $t$).
For the $i$th player's strategy variable $x_i \in \re^{n_i}$,
the $x_{i,j}$ denotes the $j$th entry of $x_i$, $j = 1, \ldots, n_i$.
Let $\nfR[\allx]$ (resp., $\cpx[\allx]$) denotes the ring of polynomials
with real (resp., complex) coefficients in $\allx$.
Let $\nfR[\allx]_d$ (resp., $\cpx[\allx]_d$) denotes
its subset of polynomials whose degrees are not greater than $d$.

For the $i$th player's strategy vector $x_i$,
the notation $\re[x_i]$, $\cpx[x_i]$, $\re[x_i]_d$, $\cpx[x_i]_d$
are defined in the same way.
For $i$th player's objective $f_i(x_i,x_{-i})$,
the notation $\nabla_{x_i}f_i$, $\nabla^2_{x_i}f_i$
respectively denote its gradient and Hessian with respect to $x_i$.

In the following, we use the letter $z$ to represent either $\allx$ or $x_i$
for convenience of discussion. Suppose $z := (z_1,\ldots,z_l)$
and $\af := (\af_1, \ldots, \af_l) \in \N^{l}$, denote
\[
z^\alpha := z_1^{\alpha_1} \cdots z_l^{\alpha_l}, \quad
|\alpha|:=\alpha_1+\ldots+\alpha_l.
\]
For an integer $d >0$, denote the monomial power set
\[
{\mathbb{N}}_d^l \, := \,
\{\alpha\in {\mathbb{N}}^l: \, \ |\alpha| \le d \}.
\]
We use $[z]_d$ to denote the vector of all monomials in $z$
and whose degree is at most $d$, ordered in the graded alphabetical ordering.
For example, if $z =(z_1, z_2)$, then
\[
[z]_3 = (1,  z_1, z_2, z_1^2, z_1z_2,
z_2^2, z_1^3, z_1^2z_2, z_1z_2^2, z_2^3).
\]
Throughout the paper, the word “generic” is used for a property
if it holds for all but a set of Lebesgue measure zero
in the space of input data.

\subsection{Ideals and positive polynomials}

Let $\mathbb{F}=\re\ \mbox{or}\ \cpx$. For a polynomial
$p \in\mathbb{F}[z]$ and subsets $I, J \subseteq \mathbb{F}[z]$,
define the product and Minkowski sum
\[
p \cdot I  :=\{ p  q: \, q \in I \}, \quad
I+J  := \{a+b: \, a \in I, b \in J  \}.
\]
The subset $I$ is an ideal if $p \cdot I \subseteq I$ for all $p\in\mathbb{F}[z]$
and $I+I \subseteq I$.
For a tuple of polynomials $q = (q_1, \ldots, q_m)$, the set
\[
\idl[q]:= q_1\cdot\mathbb{F}[z] + \ldots + q_m \cdot \mathbb{F}[z]
\]
is the ideal generated by $q$, which is the smallest ideal
containing each $q_i$.

We review basic concepts in polynomial optimization.
A polynomial $\sigma \in \re[z]$ is said to be a sum of squares (SOS)
if $\sigma=s_1^2+s_2^2+\ldots+s_k^2$ for some polynomials $s_1,\dots,s_k\in\nfR[z]$.
The set of all SOS polynomials in $z$ is denoted as $\Sigma[z]$.
For a degree $d$, we denote the truncation
\[
\Sigma[z]_d \, := \, \Sigma[z] \cap \nfR[z]_d.
\]
For a tuple $g=(g_1,\dots,g_t)$ of polynomials in $z$,
its quadratic module is the set
\[
\qmod[g] \, := \,  \Sigma[z] +  g_1 \cdot \Sigma[z] + \ldots + g_t \cdot  \Sigma[z].
\]
Similarly, we denote the truncation of $\qmod(g)$
\[
\qmod[g]_{2d} \, := \, \Sigma[z]_{2d} + g_1\cdot \Sigma[z]_{2d-\deg(g_1)}
+\cdots+g_t\cdot\Sigma[z]_{2d-\deg(g_t)}.
\]
The tuple $g$ determines the basic closed semi-algebraic set
\begin{equation}
  \label{polyrep}
\mathcal{S}(g) \, := \,  \{z \in \nfR^l: g(z) \ge  0  \}.
\end{equation}
For a tuple $h=(h_1,\dots,h_s)$ of polynomials in $\re[z]$, its real zero set is
\[
\mc{Z}(h) := \{u \in\re^l: h_1(u)=\cdots=h_s(u)=0\}.
\]
The set $\idl[h]+\qmod[g]$ is said to be {\it archimedean}
if there exists $\rho \in \idl[h]+\qmod[g]$ such that the set $\mathcal{S}(\rho)$ is compact.
If $\idl[h]+\qmod[g]$ is archimedean, then
$\mathcal{Z}(h)\cap\mathcal{S}(g)$ must be compact.
Conversely, if $\mathcal{Z}(h)\cap\mathcal{S}(g)$ is compact, say,
$\mathcal{Z}(h)\cap\mathcal{S}(g)$ is contained in the ball $R -\|z\|^2 \ge 0$,
then $\idl(h)+\qmod(g,R -\|z\|^2)$ is archimedean
and $\mathcal{Z}(h)\cap\mathcal{S}(g) = \mathcal{Z}(h)\cap\mathcal{S}(g, R -\|z\|^2)$.
Clearly, if $f \in \idl[h]+\qmod[g]$, then
$f \ge 0$ on $\mathcal{Z}(h) \cap \mathcal{S}(g)$.
The reverse is not necessarily true.
However, when $\idl[h]+\qmod[g]$ is archimedean,
if $f > 0$ on $\mathcal{Z}(h)\cap\mathcal{S}(g)$, then $f \in \idl[h]+\qmod[g]$.
This conclusion is referenced as Putinar's Positivestellensatz \cite{putinar1993positive}.
Interestingly, if $f \ge 0$ on $\mathcal{Z}(h)\cap\mathcal{S}(g)$,
we also have $f\in \idl[h]+\qmod[g]$,
under some standard optimality conditions \cite{nie2014optimality}.

\subsection{Localizing and moment matrices}

Let $\re^{ \N_{2d}^{l} }$ denote the space of all real vectors
that are labeled by $\af \in \N_{2d}^{l}$. Each
$y \in \re^{ \N_{2d}^{l} }$ is labeled as
\[
y \, = \, (y_\af)_{ \af \in \N_{2d}^{l} }.
\]
Such $y$ is called a
{\it truncated multi-sequence} (tms) of degree $2d$.
For a polynomial $f  = \sum_{ \af \in \N^l_{2d} } f_\af z^\af \in  \re[z]_{2d}$,
define the operation
\be \label{<f,y>}
\langle f, y \rangle \,= \, \sum_{ \af \in \N^l_{2d} } f_\af y_\af.
\ee
The operation $\langle f, y \rangle$ is a bilinear function in $(f, y)$.

For a polynomial $q \in \re[z]$  with $\deg(q) \le 2k$ and the integer
\[
t \,= \, k - \lceil \deg(q)/2 \rceil ,
\]
the outer product $q \cdot [z]_t ([z]_t)^T$
is a symmetric matrix polynomial in $z$, with length $\binom{n+t}{t}$.
We write the expansion as
\[
q \cdot [z]_t ([z]_t)^T \, = \, \sum_{ \af \in \N_{2d}^l }
z^\af  Q_\af,
\]
for some symmetric matrices $Q_\af$.
Then we define the matrix function
\be \label{df:Lf[y]}
L_{q}^{(d)}[y] \, := \, \sum_{ \af \in \N_{2d}^l } y_\af  Q_\af.
\ee
It is called the $k$th {\it localizing matrix} of $q$ and generated by $y$.
For given $q$, $L_{q}^{(d)}[y]$ is linear in $y$.
Clearly, if $q(u) \geq 0$ and $y = [u]_{2d}$, then
\[ L_{q}^{(d)}[y] = q(u) [u]_t[u]_t^T \succeq 0. \]
For instance, if $l=d=2$
and $q= 1 - z_1-z_1z_2$, then
\[
L_q^{(2)}[y]=\left [\begin{matrix}
y_{00}-y_{10}-y_{11} &  y_{10}-y_{20}-y_{21} &  y_{01}-y_{11}-y_{12} \\
y_{10}-y_{20}-y_{21} &  y_{20}-y_{30}-y_{31} &  y_{11}-y_{21}-y_{22} \\
y_{01}-y_{11}-y_{12} &  y_{11}-y_{21}-y_{22} &  y_{02}-y_{12}-y_{13} \\
\end{matrix}\right ].
\]
When $q$ is the constant one polynomial,
the localizing matrix $L_{1}^{(d)}[y]$
reduces to a moment matrix, which we denote as
\[
M_d[y] \, := \, L_{1}^{(d)}[y].
\]
For instance, for $n=2$ and $y \in \re^{ \N^2_4 }$, we have
$M_0[y] = [y_{00}]$,
\[
M_1[y]= \left[
\begin{array}{cccccc}
y_{00} & y_{10} & y_{01}  \\
y_{10} & y_{20} & y_{11}  \\
y_{01} & y_{11} & y_{02} \\
\end{array} \right],
\,
M_2[y]= \left[
\begin{array}{cccccc}
y_{00} & y_{10} & y_{01} & y_{20} & y_{11} & y_{02} \\
y_{10} & y_{20} & y_{11} & y_{30} & y_{21} & y_{12} \\
y_{01} & y_{11} & y_{02} & y_{21} & y_{12} & y_{03} \\
y_{20} & y_{30} & y_{21} & y_{40} & y_{31} & y_{22} \\
y_{11} & y_{21} & y_{12} & y_{31} & y_{22} & y_{13} \\
y_{02} & y_{12} & y_{03} & y_{22} & y_{13} & y_{04} \\
\end{array} \right].
\]

\section{Polynomial optimization formulations}
\label{sc:finding}

In this section, we show how to formulate efficient polynomial optimization
for solving the NEP \reff{eq:NEP}.

\subsection{Optimality conditions for individual optimization}

Consider the $i$th player's individual optimization problem $\mbf{F}_i(x_{-i})$ in (\ref{eq:subprob}),
for given $\xmi$. For convenience, we write the constraining functions as
\[
g_i(x_i) \, := \, (g_{i,1}(x_i), \ldots, g_{i,m_i}(x_i)).
\]
Suppose $\allx^* = (x_1^*, \ldots, x_N^*)$ is a NE.
Under linear independence constraint qualification condition (LICQC) at $x_i^*$,
i.e., the set of gradients for active constraining functions are linearly independent,
there exist Lagrange multipliers $\lambda_{i,j}^*$ such that
\be
\label{eq:KKTwithLM}
\left\{
\begin{array}{c}
\sum_{j=1}^{m_i}\lambda_{ij}^*\nabla_{x_i} g_{i,j}(\xpi^*) = \nabla_{x_i} f_i(\allx^*), \\
 0\leq  \lambda_{i,j}^*\perp g_{i,j}(\xpi^*) \ge 0  \, (j \in \mc{I}_i ).
\end{array}
\right.
\ee
The above is the Karush-Kuhn-Tucker (KKT) condition for the optimization $\mbf{F}_i(x_{-i}^*)$.
The $x^*$ satisfying \reff{eq:KKTwithLM} is called a KKT point.
So $\allx^*$ and $\lmd_{i,j}^*$ satisfy the following polynomial system
($i=1,\ldots, N$)
\be
\underbrace{\bbm
\nabla_{x_i} g_{i,1} & \nabla_{x_i} g_{i,2} &  \cdots &  \nabla_{x_i} g_{i,m_i} \\
g_{i,1}(x) & 0  & \cdots & 0 \\
0  & g_{i,2}(x)  & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0  &  0  & \cdots & g_{i,m_i}(x)
\ebm}_{G_i(x_i) }
\underbrace{
\bbm  \lmd_{i,1} \\ \lmd_{i,2} \\ \vdots \\ \lmd_{i,m_i} \ebm
}_{\lmd_i}
=
\underbrace{\bbm  \nabla_{x_i} f_i  \\ 0 \\ \vdots \\ 0 \ebm}_{ \hat{f}_i(\allx)} .
\ee

Therefore, if there exists a matrix polynomial $H_i(x_i)$ such that
\be \label{HiGi=Isi}
H_i(x_i) G_i(x_i)  \, = \, I_{m_i},
\ee
then we can express $\lmd_i$ as
\[
\lmd_i = H_i(x_i) G_i(x_i) \lmd_i = H_i(x_i) \hat{f}_i(\allx).
\]
Interestingly, the matrix polynomial $H_i(x_i)$ satisfying \reff{HiGi=Isi}
exists under the nonsingularity condition on $g_i$.
The polynomial tuple $g_i$ is said to be \textit{nonsingular} if
$G_i(x_i)$ has full column rank for all $x_i \in \cpx^{n_i}$ \cite{Nie2018}.

It is a generic condition. We remark that
if $g_i$ is nonsingular, then the LICQC holds at every minimizer
of \reff{eq:subprob},
so there must exist $\lambda_{i,j}$ satisfying \reff{eq:KKTwithLM}
and we can express $\lmd_{i,j}$ as
\be \label{polyexpr:lmdij}
\lmd_{i,j}  =  \lmd_{i,j}(\allx) :=
\big( H_i(x_i) \hat{f}_i(\allx) \big)_j
\ee
for all NEs.

Throughout the paper, we assume that every constraining polynomial tuple
$g_i$ is nonsingular. This is a generic assumption.
So $\lmd_{i,j}(\allx)$ can be expressed as polynomials
as in \reff{polyexpr:lmdij}. Then, each Nash equilibrium $\allx^*$
satisfies the following polynomial systems ($i=1,\ldots, N$)
\be
\label{eq:criticalpt}
(S_i): \left\{
\begin{array}{l}
\nabla_{x_i} f_i(\allx)
-\sum_{j=1}^{m_i}\lambda_{ij}(\allx) \nabla_{x_i} g_{i,j}(x_i) =0,  \\
  g_{i,j}(x_i) = 0 \, (j \in  \mc{E}_i), \,
\lmd_{i,j}(\allx)  g_{i,j}(x_i) =  0 \, (j \in  \mc{I}_i), \\
  g_{i,j}(x_i) \geq  0 \, (j \in  \mc{I}_i), \,
           \lambda_{i,j} (\allx) \geq 0 \, (j \in  \mc{I}_i) .
\end{array}
\right.
\ee
The above are necessary conditions for NEs.
When every optimization in (\ref{eq:subprob}) is convex,
the \reff{eq:criticalpt} are sufficient conditions for NEs.

\subsection{Optimization based on KKT conditions}
\label{sc:alg}

The main task for NEP is to find a tuple $\allx^* = (x_1^*, \ldots, x_N^*)$
such that each $x_i^*$ is a minimizer for the
optimization problem $\mbf{F}_i(x^*_{-i})$.
We assume each constraining tuple $g_i$ is nonsingular.
Then $\allx^*$ must satisfy the polynomial system \reff{eq:criticalpt}.
Choose a generic positive definite matrix
\[
\Theta \, \in \, \re^{ (n+1) \times  (n+1)   } .
\]
Then we consider the following optimization problem
\be
\label{eq:KKTfeasopt}
\left\{
\begin{array}{rll}
\min\limits_{\allx} & [\allx]_1^T \cdot \Theta \cdot [\allx]_1 \\
\st  &  \nabla_{x_i} f_i(\allx)
-\sum_{j=1}^{m_i}\lambda_{ij}(\allx) \nabla_{x_i} g_{i,j}(x_i) =0 \, (i \in [N]),  \\
  &  g_{i,j}(x_i) = 0 \, (j \in  \mc{E}_i, i \in [N]), \\
  & \lmd_{i,j}(\allx)  g_{i,j}(x_i) =  0 \, (j \in  \mc{I}_i, i \in [N]), \\
  &  g_{i,j}(x_i) \geq  0 \, (j \in  \mc{I}_i, i \in [N]), \\
  & \lambda_{i,j} (\allx) \geq 0 \, (j \in  \mc{I}_i, i \in [N]) .
\end{array}
\right.
\ee
In the above, the vector $[\allx]_1 = \bbm 1 & \allx^T \ebm^T \in \re^{n+1}$.
Under the nonsingularity assumptions on $g_i$, every Nash equilibrium
$\allx^*$ is a feasible point of \reff{eq:KKTfeasopt},
while the converse is typically not true. However,
for every feasible point $\allx$ of (\ref{eq:KKTfeasopt}),
the $x_i$ is a critical point for the optimization $\mbf{F}_i(x_{-i})$.
It is important to observe that if \reff{eq:KKTfeasopt} is infeasible,
then there are no NEs. If \reff{eq:KKTfeasopt} is feasible,
then it must have a minimizer, because its objective

is a positive definite quadratic function.
Moreover, for generic $\Theta$,
the minimizer of \reff{eq:KKTfeasopt} is unique.

Assume that $u := (u_1, \ldots, u_N)$ is an optimizer of (\ref{eq:KKTfeasopt}).
If each $u_i$ is a minimizer for the optimization problem $\mbf{F}_i(u_{-i})$,
then $u$ is a NE. To this end, for each player,
consider the optimization problem:
\be  \label{eq:lowerlvopt}
\left\{
\begin{array}{rl}
\omega_i := \min & f_i(x_i, u_{-i})-f_i(u_i, u_{-i}) \\
            \st  &  g_{i,j}(\xpi) = 0 \, (j \in \mc{E}_i),  \\
                 &  g_{i,j}(\xpi) \ge 0 \, (j \in \mc{I}_i) .
\end{array}
\right.
\ee
If all the optimal values $\omega_i \ge 0$,
then $u$ is a Nash Equilibrium. However, if one of them is negative,
say, $\omega_i < 0$, then $u$ is not a NE. Let $U_i$ be a set of some optimizers of
\reff{eq:lowerlvopt}, then $u$ violates the following inequalities
\be \label{cutineq:NE}
f_i(x_i, x_{-i}) \le f_i(v, x_{-i}) \quad  (v \in U_i) .
\ee
However, every Nash equilibrium must satisfy \reff{cutineq:NE}.

When $u$ is not a NE, we aim at finding a new candidate
by posing the inequalities in \reff{cutineq:NE}.
Therefore, we consider the following optimization problem:
\be  \label{eq:upperlvopt}
\left\{
  \begin{array}{rl}
    \min\limits_{\allx} &  [\allx]_1^T \cdot \Theta \cdot [\allx]_1   \\
    \st  &  \nabla_{x_i} f_i(\allx)
-\sum_{j=1}^{m_i}\lambda_{ij}(\allx) \nabla_{x_i} g_{i,j}(x_i) =0 \, (i \in [N]),  \\
  &  g_{i,j}(x_i) = 0 \, (j \in  \mc{E}_i, i \in [N]), \\
  & \lmd_{i,j}(\allx)  g_{i,j}(x_i) =  0 \, (j \in  \mc{I}_i, i \in [N]), \\
  &  g_{i,j}(x_i) \geq  0 \, (j \in  \mc{I}_i, i \in [N]), \\
  & \lambda_{i,j} (\allx) \geq 0 \, (j \in  \mc{I}_i, i \in [N]),  \\
 & f_i(v, x_{-i}) - f_i(x_i, x_{-i}) \ge 0\ (v \in \mc{K}_i, i\in[N]) .
  \end{array}
  \right.
\ee
In the above, each $\mc{K}_i$ is a set of some optimizers of
\reff{eq:lowerlvopt}.
We can solve \reff{eq:upperlvopt} again for an optimizer.
If it is verified to be a NE, then we are done.
If it is not, we can add more inequalities like \reff{cutineq:NE}.
Repeating this procedure, we get the following algorithm.

\begin{alg}
\label{ag:KKTSDP} \rm
For the NEP given as in \reff{eq:subprob} and \reff{eq:NEP},
do the following

\begin{itemize}

\item [Step~0]
Initialize $\mc{K}_i := \emptyset$ for all $i$ and $l := 0$.
Choose a generic positive definite matrix $\Theta$ of length $n+1$.

\item [Step~1]
Solve the polynomial optimization problem \reff{eq:upperlvopt}.
If it is infeasible, there is no NE and stop;
otherwise, solve it for an optimizer $u$.

\item [Step~2]
For each $i=1,\ldots,N$, solve the optimization \reff{eq:lowerlvopt}.

If all $\omega_i \ge 0$, $u$ is a NE and stop.
If one of $\omega_i$ is negative, go to the next step.

\item [Step~3]
For each $i$ with $\omega_i < 0$,
obtain a set $U_i$ of some (may not all) optimizers of \reff{eq:lowerlvopt};
then update the set $\mc{K}_i := \mc{K}_i \cup U_i$.
Let $l := l+1$, then go to Step~1.

\end{itemize}
\end{alg}

In the Step~0, we can set $\Theta = R^T R$
for a randomly generated matrix $R$ of length $n+1$.
The objective in \reff{eq:upperlvopt} is a positive definite quadratic function,
so it has a minimizer if \reff{eq:upperlvopt} is feasible.
The case is slightly different for (\ref{eq:lowerlvopt}).
If the feasible set $X_i$ is compact or $f_i(x_i, u_{-i})$
is coercive for the given $u_{-i}$,
then (\ref{eq:lowerlvopt}) has a minimizer.
If $X_i$ is unbounded and $f_i(x_i, u_{-i})$ is not coercive,
it may be difficult for computing the optimal value $\omega_i$.
In applications, we are mostly interested in cases that
(\ref{eq:lowerlvopt}) has a minimizer, for the existence of a NE.
In Section~\ref{sc:solvepop}, we will discuss how to solve
the optimization problems in Algorithm \ref{ag:KKTSDP},
by the Lasserre type Moment-SOS hierarchy of semidefinite relaxations.

The following is the convergence theorem.

\begin{theorem}  \label{tm:finiteconv}
Assume each constraining polynomial tuple $g_i$ is nonsingular and let
$\lmd_{i,j}(\allx)$ be Lagrange multiplier polynomials as in \reff{polyexpr:lmdij}.
Let $\mc{G}$ be the feasible set of \reff{eq:KKTfeasopt}
and $\mc{G}^*$ be the set of all NEs. If the complement
$\mc{G} \backslash \mc{G}^*$ is a finite set,
i.e., the cardinality $\ell := |\mc{G} \backslash \mc{G}^*| < \infty$,
then Algorithm \ref{ag:KKTSDP} must terminate within at most $\ell$ loops.
\end{theorem}
\begin{proof}
Under the nonsingularity assumption of polynomial tuples $g_i$,
the Lagrange multipliers $\lmd_{i,j}$ can be expressed as polynomials
$\lmd_{i,j}(\allx)$ as in \reff{polyexpr:lmdij}.
For each $u$ that is a feasible point of \reff{eq:KKTfeasopt},
every NE must satisfy the constraint
\[
f_i(u_i, x_{-i}) - f_i(x_i, x_{-i}) \ge 0.
\]
Therefore, every NE must also be a feasible point of \reff{eq:upperlvopt}.
Since the matrix $\Theta$ is positive definite,
the optimization \reff{eq:upperlvopt} must have a minimizer,
unless it is infeasible.
When Algorithm~\ref{ag:KKTSDP} goes to a newer loop,
say, from the $l$th to the $(l+1)$th,
the optimizer $u$ for \reff{eq:upperlvopt} in the $l$th loop
is no longer feasible for \reff{eq:upperlvopt} in the $(l+1)$th loop.
This means that the feasible set of \reff{eq:upperlvopt}
must loose at least one point after each loop,
unless a NE is reached.
Also note that the feasible set of \reff{eq:upperlvopt}
is contained in $\mc{G}$. If $\mc{G} \backslash \mc{G}^*$
is a finite set, Algorithm \ref{ag:KKTSDP} must terminate after some loops.
The number of loops is at most $\ell$.
\end{proof}

When the polynomials $f_i,g_{i,j}$ are generic,
the NEP~\reff{eq:NEP} has finitely many KKT points,
i.e., the feasible set of \reff{eq:KKTfeasopt} is finite.
This is shown in Theorem~\ref{tm:genfinite} in the Appendix.
The assumptions in Theorem~\ref{tm:finiteconv} are general.

\subsection{Convex NEPs}
\label{ssc:convex}

An important class of NEPs is that each individual optimization
$\mbf{F}_i(x_{-i})$ is a convex optimization problem, i.e.,
each objective $f_i(x_i,\xmi)$ is a convex function in $x_i$ for given $x_{-i}$,
the equality constraining function $g_{i,j}(x_i)$ ($j \in \mc{E}_i$)
is linear in $x_i$, and the inequality constraining function
$g_{i,j}(x_i)$ ($j \in \mc{I}_i$) is concave in $x_i$.
This requires that each player's strategy set $X_i$ is convex and
\begin{equation}
\label{eq:convex}
f_i(\theta a+(1-\theta)b, x_{-i}) \le \theta f_i(a, x_{-i})
+(1-\theta)f_i(b, x_{-i}),
\end{equation}
for all $a, b \in X_i$, $\theta \in [0,1]$,
and for given $x_{-i}$.

For convex optimization, the optimizers are equivalent to the KKT points,
under constraint qualification conditions (e.g., the Slater's condition or the LICQC).

In particular, when the constraining polynomial is nonsingular,
a point is a minimizer if and only if it is a KKT point, which means that
every feasible point of \reff{eq:KKTfeasopt} is a NE.
Therefore, we get the following corollary.

\begin{corollary} \label{cor:cvxnep}
Assume each $g_i$ is a nonsingular tuple of polynomials.
Suppose each $g_{i,j}(x_i)$ ($j \in \mc{E}_i$) is linear,
each $g_{i,j}(x_i)$ ($j \in \mc{I}_i$) is concave,
and each $f_i(x_i,\xmi)$ is convex in $x_i$ for given $x_{-i}$.
Then Algorithm~\ref{ag:KKTSDP} must terminate at the first loop with $l = 0$,
returning a NE or reporting that there is no NE.
\end{corollary}

For convex optimization problems, there are infinitely many optimizers
unless the optimizer is unique.
Moreover, if the objective is strictly convex
(i.e., the inequality \reff{eq:convex} holds strictly for all $a\ne b$, $0<\theta<1$),
the optimizer is always unique, if it exists.
However, these conclusions are not true for convex NEPs.
Even if each player's objective $f_i(x_i, x_{-i})$
is strictly convex in $x_i$ for all given $x_{-i}$,
the NEP might have finitely many NEs.
This is the case for Example~\ref{ep:isolatedep}.

\begin{example}

\rm
Consider the $2$-player NEP in Example~\ref{ep:isolatedep}.

Each individual optimization is strictly convex, because the Hessian's
$\nabla^2_{x_1}f_1$ and $\nabla^2_{x_2}f_2$ are positive definite.
The constraints are the convex ball conditions. The KKT system is
\be
\label{eq:1.1}
\left\{
\begin{array}{l}
 2x_{1,1}+x_{2,1}+4x_{2,2}=-2\lambda_1 x_{1,1},4x_{1,2}=-2\lambda_1 x_{1,2},\\
 x_{1,1}+2x_{1,2}+2x_{2,1}=-2\lambda_2 x_{2,1},2x_{1,1}+x_{1,2}+2x_{2,2}=-2\lambda_2 x_{2,2},\\
 \lambda_1(1-(x_{1,1})^2-(x_{1,2})^2)=0,\lambda_2(1-(x_{2,1})^2-(x_{2,2})^2)=0,\\
 1-(x_{1,1})^2-(x_{1,2})^2\ge0,1-(x_{2,1})^2-(x_{2,2})^2\ge0,\\
 \lambda_1\ge0, \lambda_2\ge 0.
\end{array}\right.
\ee
By solving the above directly, one can show that
this NEP has only $3$ NEs, together with Lagrange multipliers as follows.
\[
\begin{array}{lll}
x_1^* = (0, 0), & x_2^* = (0, 0), & \lambda_1^*=\lambda_2^*=0; \\
x_1^* = (1, 0), & x_2^* = \frac{1}{\sqrt{5}} (-1, -2 ),
& \lambda_1^*=\frac{9\sqrt{5}}{10}-1, \,\, \lambda_2^*=\frac{\sqrt{5}}{2}-1;\\
x_1^* = (-1, 0), & x_2^* = \frac{1}{\sqrt{5}}(1, 2 ), &
   \lambda_1^*=\frac{9\sqrt{5}}{10}-1,\,\, \lambda_2^*=\frac{\sqrt{5}}{2}-1.
\end{array}
\]
\end{example}

\subsection{More Nash Equilibria}
\label{sc: all}

Algorithm~\ref{ag:KKTSDP} aims at finding a single NE.
In some applications, people may be interested in more NEs.
For the case that there is a unique NE,
people are also interested in a certificate for the uniqueness.
Here, we give a procedure for computing more NEs
or verifying the uniqueness.

Assume that $\allx^*$ is a Nash Equilibrium produced by Algorithm~\ref{ag:KKTSDP},
i.e., $\allx^*$ is also a minimizer of \reff{eq:upperlvopt}.
Note that all KKT points $\overline{\allx}$ satisfying
\[
[\overline{\allx}]_1^T \Theta [\overline{\allx}]_1 \, < \,
[\allx^*]_1^T \Theta [\allx^*]_1
\]
are excluded from the feasible set of \reff{eq:upperlvopt} by the constraints
\[
f_i(u_i,\xmi)-f_i(x_i,\xmi) \,\ge \, 0  \quad
(\forall \, u \in \mc{K}_i,\, \forall \, i\in[N]).
\]

If $\allx^*$ is an isolated NE (e.g., this is the case if there are finitely many NEs),
there exists a scalar $\dt >0$ such that
\[
[\allx]_1^T \Theta [\allx]_1 \, \ge \,
[\allx^*]_1^T \Theta [\allx^*]_1 + \dt
\]
for all other NEs $\allx$. For such $\dt$, we can try to find a different NE
by solving the following optimization problem
\be
\label{eq:upperlvoptall}
\left\{ \begin{array}{rl}
\min\limits_{\allx} & [\allx]_1^T \Theta [\allx]_1 \\
\st  &  \nabla_{x_i} f_i(\allx)
-\sum_{j=1}^{m_i}\lambda_{ij}(\allx) \nabla_{x_i} g_{i,j}(x_i) =0 \, (i \in [N]),  \\
  &  g_{i,j}(x_i) = 0 \, (j \in  \mc{E}_i, i \in [N]), \\
  & \lmd_{i,j}(\allx)  g_{i,j}(x_i) =  0 \, (j \in  \mc{I}_i, i \in [N]), \\
  &  g_{i,j}(x_i) \geq  0 \, (j \in  \mc{I}_i, i \in [N]), \\
  & \lambda_{i,j} (\allx) \geq 0 \, (j \in  \mc{I}_i, i \in [N]),  \\
  & f_i(v, x_{-i}) - f_i(x_i, x_{-i}) \ge 0\ (v \in \mc{K}_i, i\in[N]),  \\
  & [\allx]_1^T \Theta [\allx]_1 \, \ge \, [\allx^*]_1^T \Theta [\allx^*]_1 + \dt .
\end{array} \right.
\ee
When an optimizer of (\ref{eq:upperlvoptall}) is computed,
we can check if it is a NE or not by solving (\ref{eq:lowerlvopt}).
If it is, we get a new NE that is different from $\allx^*$.
If it is not, we can union new points to the set $\mc{K}_i$.
Repeating the above process,  we are able to get more Nash equilibria.

A concern in computation is how to choose the constant
$\dt >0$ for \reff{eq:upperlvoptall}.
We want a value $\dt >0$ such that there is no other Nash Equilibrum $u$ such that
$[u]_1^T \Theta [u]_1 \, \le \, [\allx^*]_1^T \Theta [\allx^*]_1 + \dt $.
To this end, we consider the following maximization problem

\be \label{eq:maxxi}
\left\{
\begin{array}{rl}
\max\limits_{\allx} & [\allx]_1^T \Theta [\allx]_1   \\
\st  &  \nabla_{x_i} f_i(\allx)
-\sum_{j=1}^{m_i}\lambda_{ij}(\allx) \nabla_{x_i} g_{i,j}(x_i) =0 \, (i \in [N]),  \\
  &  g_{i,j}(x_i) = 0 \, (j \in  \mc{E}_i, i \in [N]), \\
  & \lmd_{i,j}(\allx)  g_{i,j}(x_i) =  0 \, (j \in  \mc{I}_i, i \in [N]), \\
  &  g_{i,j}(x_i) \geq  0 \, (j \in  \mc{I}_i, i \in [N]), \\
  & \lambda_{i,j} (\allx) \geq 0 \, (j \in  \mc{I}_i, i \in [N]),  \\
  & f_i(v, x_{-i}) - f_i(x_i, x_{-i}) \ge 0\ (v \in \mc{K}_i, i\in[N]),  \\
  & [\allx]_1^T \Theta [\allx]_1 \, \le \, [\allx^*]_1^T \Theta [\allx^*]_1 + \dt .\\
\end{array} \right.
\ee
Interestingly, if $\allx^*$ is also a maximizer of \reff{eq:maxxi},
then the feasible set of \reff{eq:upperlvoptall}
contains all NEs except $\allx^*$,
under some general assumptions.

\begin{prop} \label{prop:choice:dt}
Assume $\Theta$ is a generic positive definite matrix
and $\allx^*$ is a minimizer of \reff{eq:upperlvopt}.
\bit

\item [(i)]
If $\allx^*$ is also a maximizer of \reff{eq:maxxi},
then there is no other Nash Equilibrium $u$ satisfying
$[u]_1^T \Theta [u]_1 \, \le \, [\allx^*]_1^T \Theta [\allx^*]_1 + \dt $.

\item [(ii)]
If $\allx^*$ is an isolated KKT point, then there exists $\dt > 0$
such that $\allx^*$ is also a maximizer of \reff{eq:maxxi}.

\eit
\end{prop}
\begin{proof}
Note that every NE is a feasible point of \reff{eq:upperlvopt}.

(i) If $\allx^*$ is also a maximizer of \reff{eq:maxxi},
then the objective $[\allx]_1^T \Theta [\allx]_1$
achieves a constant value in the following set of \reff{eq:maxxi}.
If $u$ is a Nash equilibrium with
$[u]_1^T \Theta [u]_1 \, \le \, [\allx^*]_1^T \Theta [\allx^*]_1 + \dt$,
then
\[
\, [u]_1^T \Theta [u]_1 = [\allx^*]_1^T \Theta [\allx^*]_1 .
\]
This means that $u$ is also a minimizer of \reff{eq:upperlvopt}.
When $\Theta$ is a generic positive definite matrix,
the optimization \reff{eq:upperlvopt} has a unique optimizer,
so $u=\allx^*$.

(ii) Since $\Theta$ is positive definite, there exists $\eps > 0$ such that
\[
 [\allx]_1^T \Theta [\allx]_1  \, \ge \,  \eps( 1 + \| x \|)^2
\]
for all $\allx$. Let $C = \sqrt{ \big( [\allx^*]_1^T \Theta [\allx^*]_1 \big) / \eps }$,
then the following set
\be \nn
T \, := \,
\left\{ y = [x]_2
\left| \begin{array}{c}
  \nabla_{x_i} f_i(\allx)
-\sum_{j=1}^{m_i}\lambda_{ij}(\allx) \nabla_{x_i} g_{i,j}(x_i) =0 \, (i \in [N]),  \\
    g_{i,j}(x_i) = 0 \, (j \in  \mc{E}_i, i \in [N]), \\
   \lmd_{i,j}(\allx)  g_{i,j}(x_i) =  0 \, (j \in  \mc{I}_i, i \in [N]), \\
     g_{i,j}(x_i) \geq  0 \, (j \in  \mc{I}_i, i \in [N]), \\
   \lambda_{i,j} (\allx) \geq 0 \, (j \in  \mc{I}_i, i \in [N]),  \\
   f_i(v, x_{-i}) - f_i(x_i, x_{-i}) \ge 0\ (v \in \mc{K}_i, i\in[N]),  \\
   \, \| \allx \| \leq C
\end{array} \right. \right\}.
\ee
is compact. Note that $[\allx^*]_2 \in T$.
Let $\theta$ be the vector such that
\[
 [\allx]_1^T \Theta [\allx]_1 \,=\, \theta^T y
\]
for all $y = [x]_2$. Since $\allx^*$ is an isolated KKT point,
the $y^* := [\allx^*]_2$ is also an isolated point of $T$.
Then its subset
\[
T_1 \, := \, T \backslash \{ y^* \}
\]
is also a compact set.
Since $\allx^*$ is a minimizer of \reff{eq:upperlvopt},
the hyperplane $H :=\{ \theta^T y = \theta^T y^* \}$
is a supporting hyperplane for the set $T$.
Since $\Theta$ is generic, the optimization \reff{eq:upperlvopt}
has a unique minimizer, which implies that $y^*$ is the unique minimizer
of the linear function $\theta^Ty$ on $T$.
So, $H$ does not intersect $T_1$ and their distance is positive.
There exists a scalar $\tau > 0$ such that
\[
 [\allx]_1^T \Theta [\allx]_1 \,=\, \theta^T y \, \ge \, \theta^T y^* + \tau
  \,=\, [\allx^*]_1^T \Theta [\allx^*]_1 + \tau
\]
for all $y = [x]_2 \in T_1$.
Then, for the choice $\dt := \tau/2$, the point $\allx^*$
is the only feasible point for \reff{eq:maxxi}.
Hence, $\allx^*$ is also a maximizer of \reff{eq:maxxi}.
\end{proof}

Proposition~\ref{prop:choice:dt} shows the existence of $\dt>0$
such that \reff{eq:upperlvopt} and \reff{eq:maxxi}
have the same optimal value.
However, it does not give a concrete lower bound for $\dt$.
In computational practice, we can first give a priori value for $\dt$.
If it does not work, we can decrease $\delta$
to a smaller value (e.g., let $\delta:= \delta/5$).
By repeating this, the optimization (\ref{eq:maxxi})
will eventually have $\allx^*$ as a maximizer.
The following is the algorithm for an NE that is different from $\allx^*$.

\begin{alg}
\label{ag:allKKTSDP} \rm
Give an initial value for $\delta$ (say, $0.1$).
\begin{description}

\item [Step~1] Solve the maximization problem \reff{eq:maxxi}.
If its optimal value $\eta$ equals
$\upsilon := [\allx^*]_1^T \Theta [\allx^*]_1$,
then go to Step~2.
If $\eta$ is bigger than $\upsilon$,
then let $\delta=\min(\delta/5,\eta - \upsilon)$
and repeat this step.

\item [Step~2]
Solve the optimization problem \reff{eq:upperlvoptall}.
If it is infeasible, then there are no additional NEs;
if it is feasible, solve it for a minimizer $u$.

\item [Step~3]
For each $i=1,\ldots,N$, solve the optimization \reff{eq:lowerlvopt}
for the optimal value $\omega_i$.
If all $\omega_i \ge 0$, stop and $u$ is a NE.
If one of $\omega_i$ is negative, go to Step~4.

\item [Step~4]
For each $i\in[N]$, update the set $\mc{K}_i := \mc{K}_i \cup U_i$,
and then go back to Step~2.

\end{description}
\end{alg}

When $\allx^*$ is not an isolated KKT point,
there may not exist a satisfactory $\dt >0$ for the Step~1.
For such a case, more investigation is required to verify uniqueness of the NE
or to find other NEs.
However, for generic NEPs, there are finitely many KKT points
(see Theorem~\ref{tm:genfinite} in the appendix).
The following is the convergence result for Algorithm~\ref{ag:allKKTSDP}.

\begin{theorem}
\label{tm:allnash}
Under the same assumptions in Theorem~\ref{tm:finiteconv},
if $\Theta$ is a generic positive definite matrix
and $\allx^*$ is an isolated KKT point,
then Algorithm~\ref{ag:allKKTSDP} must terminate after finitely many steps,
returning a NE that is different from $\allx^*$ or
reporting the nonexistence of other NEs.
\end{theorem}
\begin{proof}
Under the given assumptions, Proposition~\ref{prop:choice:dt}(ii)
shows the existence of $\dt > 0$ satisfactory for the Step~1
of Algorithm~\ref{ag:allKKTSDP}.
Again, by Proposition~\ref{prop:choice:dt}(i),
the feasible set of \reff{eq:upperlvoptall} contains all
NEs except $\allx^*$. The finite termination of Algorithm~\ref{ag:allKKTSDP}
can be proved in the same way as for Theorem~\ref{tm:finiteconv}.
\end{proof}

Once a new NE is obtained, we can repeatedly apply Algorithm~\ref{ag:allKKTSDP},
to compute more NEs, if they exist. In particular, if there are finitely many NEs,
we can eventually get all of them. Indeed, for generic NEPPs,
the number of NEs is finite (see Theorem~\ref{tm:genfinite} in the appendix).
We can assume the set of equilibria is
\[
\{\allx^{(1)}, \, \ldots, \, \allx^{(s)} \} .
\]
Without loss of generality, we can assume they are ordered as
\[
[\allx^{(1)}]_1^T \Theta [\allx^{(1)}]_1 < \cdots
< [\allx^{(s)}]_1^T \Theta [\allx^{(s)}]_1 ,
\]
since $\Theta$ is generic. If the first $r$ NEs, say,
$\allx^{(1)},\ldots,\allx^{(r)}$, are obtained, there exists $\delta>0$ such that
\[
[\allx^{(j)}]_1^T \Theta [\allx^{(j)}]_1 >
    [\allx^{(r)}]_1^T \Theta [\allx^{(r)}]_1+\delta
\]
for all $j=r+1,\ldots,s$. Therefore, if we apply Algorithm~\ref{ag:allKKTSDP}
with $\allx^*=\allx^{(r)}$, the next Nash equilibrium $\allx^{(r+1)}$
can be obtained, if it exists.
Therefore, we have the following conclusion.

\begin{corollary} \label{cor:allNE}
Under the assumptions of Theorem~\ref{tm:allnash},
if there are finitely many Nash equilibria,
then all of them can be found by applying Algorithm~\ref{ag:allKKTSDP} repeatedly.
\end{corollary}

\section{The Moment-SOS hierarchy for solving optimization}
\label{sc:solvepop}

In this section, we discuss how to solve the polynomial optimization problems
in Algorithms~\ref{ag:KKTSDP} and \ref{ag:allKKTSDP}
by using the Lasserre type Moment-SOS hierarchy of semidefinite relaxations.
We assume the constraining polynomial tuples $g_i$ are all nonsingular.
Therefore, the Lagrange multipliers $\lmd_{i,j}$
can be expressed as polynomial functions $\lmd_{i,j}(\allx)$ as in \reff{polyexpr:lmdij}
for all Nash equilibria. Note that every Nash equilibrium $\allx^*$
must satisfy the polynomial system \reff{eq:criticalpt}.

\subsection{The optimization for all players}

We discuss how to solve the polynomial optimization problems
(\ref{eq:upperlvopt}), \reff{eq:upperlvoptall} and \reff{eq:maxxi},
by using the Moment-SOS hierarchy of semidefinite programming relaxations
\cite{Las01,LasBk15,LasICM,Lau09,LauICM}.

First, we discuss how to solve \reff{eq:upperlvoptall}.
Suppose the set $\mc{K}_i$ is given, for each player.
For notational convenience, denote the polynomial tuples
\begin{multline} \label{poly:Phi:ith}
\Phi_i := \Big \{\nabla_{x_i} f_i(\allx)
  -\sum_{j=1}^{m_i}\lambda_{ij}(\allx)\nabla_{x_i} g_{ij} \Big \}
\cup \Big \{ g_{i,j}: j \in \mc{E}_i \Big  \}  \\
 \cup \Big \{\lambda_{i,j}(\allx)\cdot g_{i,j}: j \in \mc{I}_i \Big  \} ,
\end{multline}
\begin{multline}  \label{poly:Psi:ith}
\Psi_i : =  \Big \{ g_{i,j}: j \in \mc{I}_i \Big  \} \cup
\Big \{\lambda_{i,j}(\allx): j \in \mc{I}_i \Big  \} \\
\cup \Big \{f_i(v,\xmi)-f_i(\xpi,\xmi): \, v \in \mc{K}_i  \Big \}  .
\end{multline}
In the above, for a vector $p = (p_1,\ldots, p_s)$ of polynomials,
the set $\{ p \}$ stands for $\{p_1, \ldots, p_s\}$, for notational convenience.
Denote the unions
\be \label{union:Phi:Psi}
\Phi \, :=\bigcup_{i=1}^N \Phi_i,  \quad
\Psi \,:=\,\bigcup_{i=1}^N \Psi_i.
\ee
They are both finite sets of polynomials.
Then, the optimization (\ref{eq:upperlvopt})
can be equivalently written as
\be  \label{eq:rewtupper}
\left\{
\baray{rl}
\vartheta_{\min} :=   \min\limits_{\allx \in \re^n} &
     \theta(\allx) := [x]_1^T \Theta [x]_1  \\
\st &  p(\allx) = 0 \,\, (\forall \, p \in \Phi),   \\
    &  q(\allx) \ge 0 \,\, (\forall \, q \in \Psi).
\earay
\right.
\ee
Denote the degree
\[
d_0 \, := \, \max\{\lceil\deg(p)/2\rceil: \, p \in \Phi \cup \Psi \}.
\]
For a degree $k \ge d_0$, consider the
the $k$th order moment relaxation for \reff{eq:rewtupper}
\be
\label{eq:d-mom}
\left\{
\baray{rl}
\vartheta_k  \, :=  \,\min\limits_{ y } & \lip \theta,y \rip\\
 \st & y_0=1,\, L_{p}^{(k)}[y] = 0 \, (p \in \Phi), \\
  & M_d[y] \succeq 0,\, L_{q}^{(k)}[y] \succeq 0 \, (q \in \Psi),  \\
  &  y \in \mathbb{R}^{\mathbb{N}^{n}_{2k}}.
\earay
\right.
\ee
Its dual optimization problem is the $k$th order SOS relaxation
\be
\label{eq:d-sos}
\left\{
\baray{ll}
\max & \gamma\\
\st & \theta -\gamma \in  \idl[\Phi]_{2k}+\qmod[\Psi]_{2k}. \\
\earay
\right.
\ee
For relaxation orders $k=d_0, d_0+1, \ldots$,
we get the Moment-SOS hierarchy of semidefinite relaxations
\reff{eq:d-mom}-\reff{eq:d-sos}.
This produces the following algorithm for solving
the polynomial optimization problem \reff{eq:rewtupper}.

\begin{alg} \label{ag:upsdp} \rm
Let $\theta,\Phi,\Psi$ be as in \reff{eq:rewtupper}.
Initialize $k := d_0$.

\begin{description}

\item [Step~1]
Solve the moment relaxation (\ref{eq:d-mom}). If it is infeasible,
\reff{eq:rewtupper} has no feasible points and stop;
otherwise, solve it for a minimizer $y^*$.

\item [Step~2]
Let $u = (y^*_{e_1}, \ldots, y^*_{e_n})$.
If $u$ is feasible for \reff{eq:rewtupper}
and $\vartheta_k = \theta(u)$,
then $u$ is a minimizer of \reff{eq:rewtupper}.
Otherwise, let $k := k+1$ and go to Step~1.

\end{description}
\end{alg}

In the Step~2, $e_i$ denotes the labeling vector such that
the $i$th entry is $1$ while all other entries are $0$.
For example, when $n=4$, then $y_{e_2} = y_{0100}$.

The conclusions of Algorithm~\ref{ag:upsdp} are justified as follows.
The optimization \reff{eq:d-mom} is a relaxation of \reff{eq:rewtupper}.
This is because if $x$ is a feasible point of \reff{eq:rewtupper},
then $y = [x]_{2k}$ must be feasible for \reff{eq:d-mom}.
Hence, if \reff{eq:d-mom} is infeasible,
then \reff{eq:rewtupper} must be infeasible,
which also implies the nonexistence of a NE.
Moreover, the optimal value $\vartheta_k$ of \reff{eq:d-mom}
is a lower bound for the minimum value of \reff{eq:rewtupper},
i.e., $\vartheta_k \leq \theta(x)$ for all $x$
that is feasible for \reff{eq:rewtupper}.
In the Step~2, if $u$ is feasible for \reff{eq:rewtupper}
and $\vartheta_k = \theta(u)$,
then $u$ must be a minimizer of \reff{eq:rewtupper}.
The convergence of Algorithm~\ref{ag:upsdp} is shown as follows.

\begin{theorem} \label{thmcvg:upconv}
Assume the matrix $\Theta$ is a generic positive definite matrix
and $\idl[\Phi]+\qmod[\Psi]$ is archimedean.
\bit

\item [(i)]
If the optimization \reff{eq:rewtupper} is infeasible,
then the moment relaxation \reff{eq:d-mom}
must be infeasible when the order $k$ is big enough.

\item [(ii)]
Suppose the optimization \reff{eq:rewtupper} is feasible.
Let $u^{(k)}$ be the point $u$ produced in the Step~2 of
Algorithm~\ref{ag:upsdp} in the $k$th loop.
Then $u^{(k)}$ converges to the unique minimizer of \reff{eq:rewtupper}.
In particular, if the real zero set of $\Phi$ is finite,
then $u^{(k)}$ is the unique minimizer of \reff{eq:rewtupper},
when $k$ is sufficiently large.

\eit
\end{theorem}
\begin{proof}
(i) If \reff{eq:rewtupper} is infeasible,  the constant polynomial $-1$
can be viewed as a positive polynomial on the feasible set of \reff{eq:rewtupper}.
Since $\idl[\Phi]+\qmod[\Psi]$ is archimedean,
we have $-1 \in  \idl[\Phi]_{2k}+\qmod[\Psi]_{2k}$,
for $k$ big enough, by the Putinar's Positivstellensatz \cite{putinar1993positive}.
For such big $k$, the SOS relaxation \reff{eq:d-sos} is unbounded from above,
hence the moment relaxation \reff{eq:d-mom} must be infeasible.

(ii) When the optimization \reff{eq:rewtupper} is feasible,
it must have a unique minimizer, say, $\allx^*$,
because its objective is a generic positive definite quadratic polynomial.
The convergence of $u^{(k)}$ to $\allx^*$ is shown in
\cite{Swg05} or \cite[Theorem~3.3]{nie2013certifying}.
For the special case that $\Phi(\allx)=0$ has finitely many real solutions,
the point $u^{(k)}$ must be equal to $\allx^*$,
when $k$ is large enough.
This is shown in \cite{lasserre2008semidef}
(also see \cite{nie2013polynomial}).
\end{proof}

The archimedeanness of $\idl[\Phi]+\qmod[\Psi]$
is essentially requiring that the feasible set of \reff{eq:rewtupper}
is compact. If the real zero set of $\Phi$ is compact,
then $\idl[\Phi]+\qmod[\Psi]$ must be archimedean.
In particular, if the NEPP has finitely many real KKT points,

then $\idl[\Phi]+\qmod[\Psi]$ is archimedean.
Interestingly, when the objective and constraining polynomials are generic,
there are finitely many KKT points.
See Theorem~\ref{tm:genfinite} in the appendix.
In fact, as shown in the proof of Theorem~\ref{tm:genfinite},
the zero set of $\Phi$ is finite for generic NEPPs,
and hence Algorithm~\ref{ag:upsdp} has finite convergence.

The other polynomial optimization problem \reff{eq:upperlvoptall}
can be solved in the same way by the Moment-SOS hierarchy of
semidefinte relaxations like \reff{eq:d-mom}-\reff{eq:d-sos}.
The convergence property is the same.
For cleanness of the paper, we omit the details.

For the maximization \reff{eq:maxxi},
we denote the set of polynomials
\be \label{newdef:Psi}
\Phi \, :=\bigcup_{i=1}^N \Phi_i , \quad
\Psi \, :=\bigcup_{i=1}^N \Psi_i \cup
\Big\{
[\allx^*]_1^T \Theta [\allx^*]_1 + \dt - [\allx]_1^T \Theta [\allx]_1
\Big\} .
\ee
Then \reff{eq:maxxi} can be equivalently written as
\be  \label{eq:rewtuppermax}
\left\{
\baray{rl}
 \vartheta_{\min} := \min\limits_{\allx \in \re^n} &
     \theta(\allx) := -[x]_1^T \Theta [x]_1  \\
\st &  p(\allx) = 0 \,\, (\forall \, p \in \Phi),   \\
    &  q(\allx) \ge 0 \,\, (\forall \, q \in \Psi).
\earay
\right.
\ee
Similarly, the hierarchy of moment relaxations~\reff{eq:d-mom}
can be used to solve \reff{eq:rewtuppermax}.
The following is the algorithm.

\begin{alg} \label{ag:maxsdp} \rm
Let $\Phi,\Psi,\theta$ be as in \reff{newdef:Psi}-\reff{eq:rewtuppermax}.
Initialize $k := d_0$.
\begin{description}

\item [Step~1]
Solve the moment relaxation (\ref{eq:d-mom})
for the minimum value $\vartheta_k$ and a minimizer $y^*$ and let $t:=d_0$.

\item [Step~2]
If $\vartheta_k\ge-[\allx^*]_1^T \Theta [\allx^*]_1$, stop;
otherwise, go to Step~3.

\item [Step~3]
If $y^*$ satisfies the rank condition
\be \label{eq:flatrank}
 \Rank{M_t[y^*]} \,=\, \Rank{M_{t-d_0}[y^*]} ,
\ee
then extract a set $U_i$ of
$r :=\Rank{M_t(y^*)}$ minimizers for \reff{eq:rewtuppermax}) and stop.

\item [Step~4]
If \reff{eq:flatrank} fails to hold and $t < k$,
let $t := t+1$ and then go to Step~3;
otherwise, let $k := k+1$ and go to Step~1.

\end{description}
\end{alg}

The optimization \reff{eq:rewtuppermax} is always feasible because
$\allx^*$ is a feasible point.
Therefore, the moment relaxation \reff{eq:d-mom} is also feasible.
Since the minimum value $\vartheta_k$ is a lower bound of $\vartheta_{\min}$,
if $\vartheta_k\ge -[\allx^*]_1^T \Theta [\allx^*]_1$, then
\[
\vartheta_k=\vartheta_{\min}=-[\allx^*]_1^T \Theta [\allx^*]_1
\]
and hence $\allx^*$ is a maximizer of \reff{eq:maxxi}.

In Step~3, the rank condition~\reff{eq:flatrank} is called \textit{flat truncation} \cite{nie2013certifying}.
It is a sufficient (and almost necessary) condition to check convergence
of moment relaxations.

When \reff{eq:flatrank} holds, the method in \cite{HenLas05}
can be used to extract $r$ minimizers for \reff{eq:rewtuppermax}).
The method is implemented in the software {\tt GloptPoly}~\cite{GloPol3}.
Moreover, Algorithms~\ref{ag:upsdp} and \ref{ag:maxsdp}
can be implemented in {\tt GloptPoly}.

The convergence of Algorithm~\ref{ag:maxsdp} is as follows.
Note that $\mbox{Ideal}({\Phi})+\qmod({\Psi})$ is archimedean,
since it contains the polynomial
$[\allx^*]_1^T \Theta [\allx^*]_1 + \dt - [\allx]_1^T \Theta [\allx]_1$.
Therefore, we always have $\vartheta_k \to \vartheta_{\min}$ as $k \to \infty$ \cite{Las01}.
Under some classical optimality conditions,
we have $\vartheta_k = \vartheta_{\min}$
when $k$ is large enough \cite{nie2014optimality}.

Moreover, if the real zero set of $\Phi$ is finite,
then Algorithm~\ref{ag:maxsdp} has finite convergence \cite{nie2013polynomial}.

\subsection{Checking Nash equilibria}
\label{sc:loweropt}

Suppose $u$ is a minimizer of \reff{eq:upperlvopt}.
To check if $u=(u_i,u_{-i})$ is a NE or not, we need to solve the individual
optimization \reff{eq:lowerlvopt} for each player.
For notational convenience, we denote the polynomial tuples
\begin{multline} \label{polylow:Hi(u)}
H_i(u) := \big \{g_{i,j}: j \in \mc{E}_i \big \} \cup
\big \{ \lambda_{i,j}(\xpi,u_{-i})\cdot g_{i,j}: j \in \mc{I}_i \big \}  \\
\cup \big \{ \nabla_{x_i}f_i(\xpi, u_{-i}) -
 \sum\limits_{j=1}^{m_i}\lambda_{ij}(\xpi, u_{-i} ) \nabla_{x_i} g_{ij} \big \} ,
\end{multline}
\be  \label{polylow:Gi(u)}
G_i(u) := \big \{g_{i,j}: j \in \mc{I}_i \big \} \cup
\big \{ \lambda_{i,j}(\xpi,u_{-i}) : j \in \mc{I}_i \big \}   .
\ee
Like the earlier case, the set $\{ p \}$ stands for $\{p_1, \ldots, p_s\}$,
when $p = (p_1,\ldots, p_s)$ is a vector of polynomial.
The sets $H_i(u), G_i(u)$ are finite collections
of polynomials in $x_i$ and parameterized by $u$.
If the optimization (\ref{eq:lowerlvopt})
has a minimizer, then it is equivalent to the following optimization
\be  \label{eq:rewtlower}
\left\{
\baray{rrl}
\omega_i:= &\min\limits_{x_i \in \re^{n_i} }  & f_i(\xpi, u_{-i}) \\
&\st & p(\xpi)=0 \, (p \in H_i(u) ), \\
    && q(\xpi) \ge 0  \, (q \in G_i(u) ).
\earay
\right.
\ee
The above is a polynomial optimization problem in $x_i$.
Denote the degree for its constraining polynomials
\be \label{hfdg:di:ith}
d_{i}:=\max \big \{ \lceil \deg(p)/2 \rceil :
     p \in H_i(u) \cup G_i(u)  \big \}.
\ee
For a degree $k \ge d_{i}$, the $k$th order moment relaxation for
\reff{eq:rewtupper} is
\be
\label{eq:d-momlow}
\left\{
\baray{rrl}
 \omega_i^{(k)}:=&\min\limits_{ y } & \lip f_i(\xpi, u_{-i}),y \rip\\
&\st & y_0=1,\, L_{p}^{(k)}[y] = 0 \, (p \in H_i(u)),  \\
 & & M_k[y] \succeq 0,\, L_{q}^{(k)}[y] \succeq 0 \, (q \in G_i(u)), \\
 & &  y\in\mathbb{R}^{\mathbb{N}^{n_i}_{2k}}.
\earay
\right.
\ee
Its dual optimization problem is the $k$th order SOS relaxation
\be
\label{eq:d-soslow}
\left\{
\baray{ll}
\max & \gamma\\
\st & f_i(\xpi, u_{-i}) - \gamma \in
        \idl[H_i(u)]_{2k}+\qmod[G_i(u)]_{2k}. \\
\earay
\right.
\ee
By solving the above relaxations for $k =d_i, d_i+1, \ldots$,
we get the Moment-SOS hierarchy of relaxations
\reff{eq:d-momlow}-\reff{eq:d-soslow}.
This gives the following algorithm.

\begin{alg} \label{ag:lowsdp} \rm
For the $i$th player's individual optimization \reff{eq:rewtlower},
assume $u$ is a minimizer of \reff{eq:upperlvopt}.
\begin{description}

\item [Step~0]
Construct the sets $H_i(u)$, $G_i(u)$ of polynomials as in
\reff{polylow:Hi(u)}, \reff{polylow:Gi(u)}.
Initialize $k := d_i$.

\item [Step~1]
Solve the moment relaxation (\ref{eq:d-momlow})
for the minumum value $\omega_i^{(k)}$
and a minimizer $y^*$. If $\omega_i^{(k)}\ge0$,
then $\omega_i=0$ and stop;
otherwise, go to the next step.

\item [Step~2]
Let $t:=d_i$ as in \reff{hfdg:di:ith}.
If $y^*$ satisfies the rank condition
\be \label{eq:flatranklow}
 \Rank{M_t[y^*]} \,=\, \Rank{M_{t-d_i}[y^*]} ,
\ee
then extract a set $U_i$ of
 $r :=\Rank{M_t(y^*)}$ minimizers for \reff{eq:rewtlower}) and stop.

\item [Step~3]
If \reff{eq:flatranklow} fails to hold and $t < k$,
let $t := t+1$ and then go to Step~2;
otherwise, let $k := k+1$ and go to Step~1.

\end{description}
\end{alg}

We would like to remark that the optimization \reff{eq:rewtlower}
is always feasible, because $u$ is a minimizer of \reff{eq:upperlvopt}.
The moment relaxation (\ref{eq:d-momlow}) is also feasible.
Because $\omega_i^{(k)}$ is a lower bound for $\omega_i$,
and $\omega_i\le f_i(u_i,u_{-i})=0$,
if $\omega_i^{(k)}\ge0$, then $\omega_i$ must be $0$.

The Algorithm~\ref{ag:lowsdp} can be implemented in {\tt GloptPoly}.

The following theorem is the convergence for Algorithm~\ref{ag:lowsdp}.
Its proof follows from \cite[Theorem~4.4]{Nie2018saddle}.

\begin{theorem}  \label{thmcvg:lowopt}
Assume the $i$th player's constraining polynomial tuple $g_i$
is nonsingluar and its optimization \reff{eq:lowerlvopt}
has a minimizer for the given $u_{-i}$.
Assume either one of the following conditions hold:
\bit

\item [(i)]
The set $\idl[H_i(u)]+\qmod[G_i(u)]$ is archimedean,

\item [(ii)]
The real zero set of polynomials in $H_i(u)$ is finite.

\eit
If each minimizer of \reff{eq:lowerlvopt} is an isolated critical point,
then all minimizers of \reff{eq:d-momlow}
must satisfy the flat truncation \reff{eq:flatranklow},
for all $k$ big enough.
Therefore, Algorithm~\ref{ag:lowsdp} must terminate within
finitely many loops.
\end{theorem}

We would like to remark the following inclusion
\[
\idl[g_{i,j}: j \in \mc{E}_i] \subseteq \idl[H_i(u)], \quad
\qmod[g_{i,j}: j \in \mc{I}_i] \subseteq \qmod[G_i(u)].
\]
If $\idl[g_{i,j}: j \in \mc{E}_i] + \qmod[g_{i,j}: j \in \mc{I}_i]$
is archimedean, then $\idl[H_i(u)]+\qmod[G_i(u)]$ is also archimedean.
Therefore, if the archimedeanness holds for the $i$th player's optimization
\reff{eq:subprob}, then the condition (i) in Theorem~\ref{thmcvg:lowopt}
is satisfied.

\section{Numerical Experiments}
\label{sc:num}

This section reports numerical experiments for solving NEPs
by using Algorithm~\ref{ag:KKTSDP} and \ref{ag:allKKTSDP}.
We apply the software {\tt GloptiPoly~3} \cite{GloPol3}
and {\tt SeDuMi} \cite{sturm1999using} to solve
the Moment-SOS relaxations for the polynomial optimization
\reff{eq:upperlvopt}, \reff{eq:lowerlvopt},
\reff{eq:upperlvoptall} and \reff{eq:maxxi}.
The computation is implemented in an Alienware Aurora R8 desktop,
with an Intel\textsuperscript \textregistered \
Core(TM) i7-9700 CPU at 3.00GHz$\times $8 and 16GB of RAM,
in a Windows 10 operating system. To implement Algorithm~\ref{ag:KKTSDP},
we need Lagrange multiplier representations as in \cite{Nie2018}.
The following cases are frequently used.

\bit

\item For the constraint $\{x_i\in\re^{n_i}:\sum_{j=1}^{n_i}x_{i,j}\le 1,x_i\ge0\}$,
the constraining polynomials are
$g_{i,0}=1-\sum_{j=1}^{n_i}x_{i,j},g_{i,1}=x_{i,1},\cdots,g_{i,n_i}=x_{i,n_i}.$
The Lagrange multipliers $\lmd_{i,j}$ can be represented as
\[
\lmd_{i,0} \,=\, x_i^T\nabla_{x_i} f_i, \quad
\lambda_{i,j} \, = \, \frac{\partial f_i}{\partial x_{i,j}}-x_i^T\nabla_{x_i} f_i ,
\quad j = 1,\ldots, n_i.
\]

\item For the sphere constraint $1-x_i^Tx_i = 0$
or the ball constraint $1-x_i^Tx_i \ge 0 $,
the constraining polynomial is $g_{i,1}=1-x_i^Tx_i$
and the Lagrange multiplier can be expressed as
$\lambda_{i,1} \, = \, -\frac{1}{2}x_i^T\nabla_{x_i} f_i.$

\eit

In Step~2 of Algorithm~\ref{ag:KKTSDP} and Step~3 of Algorithm~\ref{ag:allKKTSDP},
if the optimal value $\omega_i \ge 0$ for all players,
then the point $u$ is a NE. In numerical computation,
we cannot have $\omega_i \ge 0$ exactly, due to round-off errors.
Therefore, we use the parameter
\[
\omega^* \, := \, \min_{i=1,\ldots,N} \omega_i
\]
to measure the accuracy of the computed NE.
Typically, if $\omega^*$ is small, say, $\omega^* \ge -10^{-6}$,
we regard the computed solution as an accurate NE.

For the convex NEP in Example~\ref{ep:isolatedep},
Algorithm~\ref{ag:allKKTSDP} found all the $3$ NEs correctly
with $\omega^*=-1.9512 \cdot 10^{-9}$.
The following is another example of convex NEPs.

\begin{example}
\label{ep:convex}
\rm
Consider the convex NEP
\begin{equation*}
\mbox{1st player:} \left\{
  \begin{array}{cl}
     \min\limits_{x_1 \in \re^2} &  x_{1,1}(x_{1,1}+
     x_{2,1}+4x_{2,2})+4x_{1,2}^2,  \\
\st & 1-(x_{1,1})^2-(x_{1,2})^2 \ge 0,
  \end{array}\right.
\end{equation*}
\begin{equation*}
\mbox{2nd player:}   \left\{
\begin{array}{cl}
\min\limits_{x_2 \in \re^2}& 2x_{2,1}^2+2x_{2,2}^2+(x_{1,1}-2x_{1,2})x_{2,1}\\
     & \qquad \qquad \qquad +(4x_{1,1}+x_{1,2})x_{2,2}, \\
    s.t. & 1-x_{2,1}-x_{2,2} \ge 0, \, x_{2,1} \ge 0, \, x_{2,2} \ge 0 .
  \end{array}\right.
\end{equation*}
For the second player,
the Lagrange multipliers can be represented as
\[
\lambda_{2,1}=-\nabla_{x_2} f_2^Tx_2,\
\lambda_{2,2}=\frac{\partial f_2}{\partial x_{2,1}}+\lambda_{2,1},\
\lambda_{2,3}=\frac{\partial f_2}{\partial x_{2,2}}+\lambda_{2,1}.
\]
For this NEP, Algorithm~\ref{ag:allKKTSDP}, found two NEs:
\[
x_1^* =(0.0000,-0.0000),\, x_2^* = (0.0000, 0.0000), \,
\omega^* = -7.4772\cdot 10^{-9};
\]
\[
x_1^* =(-1.0000,-0.0000), \,  x_2^* = (0.1250,0.8750), \,
\omega^* = -3.3640 \cdot 10^{-8}.
\]
The computation took about $1$ second.
\end{example}

\begin{example}
\label{ep:nonconvex2x3}
\rm
Consider the $2$-player NEP
\begin{equation*}
\mbox{1st player:} \left\{
  \begin{array}{cl}
    \min\limits_{x_1 \in \re^3} & \sum_{j=1}^3 x_{1,j}(x_{1,j}- j\cdot x_{2,j})\\
    s.t. & 1-x_{1,1}x_{1,2} \ge 0, \,1-x_{1,2}x_{1,3}\ge 0,\, x_{1,1} \ge 0,
  \end{array}\right.
\end{equation*}
\begin{equation*}
\mbox{2nd player:}  \left\{
  \begin{array}{cl}
    \min\limits_{x_2 \in \re^3} &
    \prod_{j=1}^3 x_{2,j} +
    \sum\limits_{\mathclap{\substack{1\le i<j\le 3\\1\le k\le 3}}} x_{1,i}x_{1,j}x_{2,k}+
    \sum\limits_{\mathclap{\substack{1\le i\le 3\\1\le j<k\le 3}}} x_{1,i}x_{2,j}x_{2,k}
    \\
    s.t. & 1-(x_{2,1})^2 - (x_{2,2})^2 = 0.
  \end{array}\right.
\end{equation*}
The first player's optimization is non-convex, with an unbounded feasible set.
The Lagrange multipliers for the first player's optimization are
\be \nn
\begin{array}{lll}
  \lambda_{1,1}=(1-x_{1,1}x_{1,2})\frac{\partial f_1}{\partial x_{1,1}}, &
  \lambda_{1,2}=-x_{1,1}\frac{\partial f_1}{\partial x_{1,2}}, &
  \lambda_{1,3}=x_{1,1}\frac{\partial f_1}{\partial x_{1,1}}
      -x_{1,2}\frac{\partial f_1}{\partial x_{1,2}}.
\end{array}
\ee
Applying Algorithm \ref{ag:allKKTSDP}, we get four NEs:
\[
\begin{array}{ll}
x_1^*=(0.3198,  0.6396, -0.6396 ), &x_2^*=( 0.6396,    0.6396,    -0.4264); \\
x_1^*=(0.0000, 0.3895, 0.5842 ), & x_2^*=(-0.8346, 0.3895, 0.3895); \\
x_1^*=(0.2934, -0.5578, 0.8803 ), & x_2^*=( 0.5869, -0.5578,  0.5869);\\
x_1^*=(0.0000, -0.5774, -0.8660 ), & x_2^*=( -0.5774, -0.5774, -0.5774). \\
\end{array}
\]
Their accuracy parameters are respectively
\[
-7.1879\cdot 10^{-8}, \, -3.5040\cdot 10^{-7}, \, -4.3732\cdot 10^{-7},
-6.4360\cdot 10^{-7}.
\]
It took about $30$ seconds.
If the second player's objective becomes
\[
-\prod_{j=1}^3 x_{2,j} +
\sum\limits_{\mathclap{\substack{1\le i\le 3\\1\le j<k\le 3}}} x_{1,i}x_{2,j}x_{2,k}
-\sum\limits_{\mathclap{\substack{1\le i<j\le 3\\1\le k\le 3}}} x_{1,i}x_{1,j}x_{2,k},
\]
then there is no NE, which is detected by Algorithm~\ref{ag:KKTSDP}.
It took around $16$ seconds.
\end{example}

\begin{example} \rm
\label{ep:plain3x2ball_arc_box}
Consider the $3$-player NEP
\begin{equation*}
\mbox{1st player:}  \left\{
 \begin{array}{cl}
 \min\limits_{x_1 \in \re^2} &  (2x_{1,1}-x_{1,2}+3)x_{1,1}x_{2,1}\\
 							&\qquad \qquad \qquad+[(2x_{1,2})^2+(x_{3,2})^2]x_{1,2}\\

\st  & 1 - x_1^Tx_1 \ge 0,
    \end{array}\right.
\end{equation*}
\begin{equation*}
\mbox{2nd player:}  \left\{
 \begin{array}{cl}
 \min\limits_{x_2 \in \re^2 } &  [(x_{2,1})^2-x_{1,2}]x_{2,1} \\
 &\qquad \quad+[(x_{2,2})^2+2x_{3,2}+x_{1,2}x_{3,1}]x_{2,2}\\

 \st &x_2^Tx_2-1 =0, \, x_{2,1} \ge 0, \, x_{2,2} \ge 0,
 \end{array}\right.
\end{equation*}
\begin{equation*}
\mbox{3rd player:}  \left\{
 \begin{array}{cl}
 \min\limits_{x_3 \in \re^2 }& (x_{1,1}x_{1,2}-1)x_{3,1}- [3(x_{3,2})^2+1]x_{3,2}\\
 &\qquad \qquad \qquad\quad+2[x_{3,1}+x_{3,2}]x_{3,1}x_{3,2} \\

\st & 1-(x_{3,1})^2 \ge 0, \, 1 - (x_{3,2})^2 \ge 0.
 \end{array}\right.
\end{equation*}

The Lagrange multipliers can be represented as
\[
\begin{array}{lll}
  	\lambda_{2,1}=\frac{1}{2}(x_2^T\nabla_{x_2} f_2),
  & \lambda_{2,2}=\frac{\partial f_2}{x_{2,1}}-2x_{2,1}\lambda_{2,1},
  &\lambda_{2,3}=\frac{\partial f_2}{x_{2,2}}-2x_{2,2}\lambda_{2,1},\\
  	\lambda_{3,1}=-\frac{x_{3,1}}{2}\frac{\partial f_3}{\partial x_{3,1}},
  &\lambda_{3,2}=-\frac{x_{3,2}}{2}\frac{\partial f_3}{\partial x_{3,2}}.
\end{array}
\]
Applying Algorithm~\ref{ag:allKKTSDP},  we get the unique NE
\[
x_1^* = (-0.3558, -0.9346), \,
x_2^* = (1.0000,   0.0000), \,
x_3^* = (-0.3331,  1.0000).
\]
The accuracy parameter is $-9.2310\cdot 10^{-9}.$
It took around $9$ seconds.
If the third player's objective becomes $-f_1(\allx)-f_2(\allx)$,
then the NEP becomes a zero-sum game and there is no NE.
Algorithm~\ref{ag:KKTSDP} detected the nonexistence.
It took around $3$ seconds.
\end{example}

\begin{example} \rm
\label{ep:ring}
Consider the $2$-player NEP
\begin{equation*}
\mbox{1st player:} \left\{
\begin{array}{cll}
  \min\limits_{x_1 \in \re^2} &  2x_{1,1}x_{1,2}+
  3x_{1,1}(x_{2,1})^2+3(x_{1,2})^2x_{2,2}\\
 \st &(x_{1,1})^2+(x_{1,2})^2-1\ge 0,  \\
  	 &2-(x_{1,1})^2-(x_{1,2})^2\ge 0
  \end{array}\right.
\end{equation*}
\begin{equation*}
\mbox{2nd player:} \left\{
\begin{array}{cll}
  \min\limits_{x_2 \in \re^2 }& (x_{2,1})^3+(x_{2,2})^3+x_{1,1}(x_{2,1})^2 \\
  & \qquad \qquad +x_{1,2}(x_{2,2})^2 +x_{1,1}x_{1,2}(x_{2,1}+x_{2,2}) \\
  \st &(x_{2,1})^2+(x_{2,2})^2-1 \ge 0, \\
      &2-(x_{2,1})^2+(x_{2,2})^2 \ge 0.
  \end{array}\right.
\end{equation*}
The Lagrange multipliers can be represented as ($i=1,2$):
\[
\lambda_{i,1}=\frac{1}{2}\nabla_{x_i} f_i^Tx_i(2-x_i^Tx_i),\quad
 \lambda_{i,2}=\frac{1}{4}\nabla_{x_i} f_i^Tx_i(1-x_i^Tx_i).
 \]
By Algorithm \ref{ag:allKKTSDP}, we get the unique NE
\[
x_1^* =(-1.3339,0.4698), \quad x_2^* = (-1.4118,0.0820),
\]
with the accuracy parameter $-3.5186\cdot 10^{-8}.$
It took around $5$ seconds.
\end{example}

\begin{example}
\label{ep:bs}
\rm
Consider the NEP
\begin{equation*}
\baray{l}
\mbox{1st player:}  \left\{
\begin{array}{cl}
   \min\limits_{x_1 \in \re^{n_1}} &
   \sum_{1 \le i \le j \le n_1} x_{1,i}x_{1,j}(x_{2,i}+x_{2,j})  \\
\st   & 1-(x_{1,1}^2+\cdots+x_{1,n_1}^2) = 0,
\end{array}\right.
\\
\mbox{2nd player:}   \left\{
\begin{array}{cl}
    \min\limits_{x_2 \in \re^{n_2}} &
    \sum_{1 \le i\le j \le n_2} x_{2,i}x_{2,j}(x_{1,i}+x_{1,j})   \\
\st & 1-(x_{2,1}^2+\cdots+x_{2,n_2}^2) = 0,
\end{array}\right.
\earay
\end{equation*}
where the dimension $n_1 = n_2$.
The computational results for cases
$n_1=n_2=3,4,5,6$ are shown in Table~\ref{tab:exm:bs}.
\begin{table}[htbp]
\centering
\caption{Computational results for Example~\ref{ep:bs}.}
\label{tab:exm:bs}
\begin{tabular}{|l|c|c|r|}  \hline
$n_1$  & $x_1^*$ 	 	& $x_2^*$ 		& time \quad   \\ \hline
$3$  & $(-0.5774,\, -0.5774, \, -0.5774)$ & $(-0.5774,\, -0.5774, \, -0.5774)$ & 1.31  \\ \hline
$4$  & $\begin{array}{l}(0.8381, \,0.5024,\\ \qquad -0.0328 , \,-0.2098)\end{array}$
     & $\begin{array}{l} (-0.1791, \, -0.0683,\\ \qquad 0.4066, \,0.8933) \end{array}$
           & 62.85   \\ \hline
$5$  & $\begin{array}{l}(0.8466,\, 0.4407,\, 0.1744,\\
                 \qquad  \qquad \qquad -0.0101,\, -0.2418) \end{array}$
     & $\begin{array}{l}(-0.1944,\,-0.0512, \, 0.1238,\\
                     \qquad\qquad 0.3370, \, 0.9114) \end{array}$
              & 682.67 \\ \hline
$6$  & $\begin{array}{l}(0.8026,\, 0.4724,\, 0.1799, \\
           \qquad 0.1799, \, -0.0637,\, -0.2527)\end{array}$
     & $\begin{array}{l}(-0.1979,\, -0.0772,\, 0.1091,\\
            \qquad  0.1091, \, 0.4040, \, 0.8762) \end{array}$
                & 18079.99  \\ \hline
\end{tabular}
\label{test}
\end{table}
The time is displayed in seconds.
The accuracy parameter $\omega^*$ is respectively
\[
-1.0689\cdot 10^{-7}, \quad  -1.4459\cdot 10^{-9}, \quad
-2.7551\cdot 10^{-9}, \quad  -7.0354\cdot 10^{-9} .
\]
Because of the relatively large amount of computational time,
we only compute one NE for each case in the above.
\end{example}

We would like to remark that our method can also be applied to solve
NEPs for which the individual optimization has no constraints,
or equivalently, the feasible set $X_i$ for \reff{eq:subprob}
is the entire space $\re^{n_i}$. For unconstrained NEPs,
the KKT system \reff{eq:KKTwithLM} becomes
\[
\nabla_{x_i} f_i(x^*) \, = \, 0, \quad i=1,\ldots, N.
\]
The Algorithms~\ref{ag:KKTSDP} and \reff{ag:allKKTSDP}
can be implemented in the same way.

\begin{example}
\label{ep:unconstrained}
\rm
Consider the unconstrained NEP
\begin{equation*}
\baray{l}
\mbox{1st player:}   \left\{
\begin{array}{cl}
   \min   &  \sum\limits_{i=1}^{n_1} (x_{1,i})^4 +
    \sum\limits_{0\le i\le j \le k  \le n_1}
     \frac{x_{1,i}x_{1,j} ( x_{1,k} + x_{2,i}+ x_{3,j})}{(n_1)^2} \\
 \st & x_1 \in \re^{n_1} ,
\end{array} \right.
\\
\mbox{2nd player:}  \left\{
\begin{array}{cl}
   \min  &  \sum\limits_{i=1}^{n_2} (x_{2,i})^4 +
   \sum\limits_{0\le i\le j \le k  \le n_2}
     \frac{x_{2,i}x_{2,j} ( x_{2,k} + x_{3,i}+ x_{1,j} )}{(n_2)^2} \\
\st   &   x_2 \in \re^{n_2},
\end{array}\right.
\\
\mbox{3rd player:}  \left\{
\begin{array}{cl}
   \min  &  \sum\limits_{i=1}^{n_3} (x_{3,i})^4 +
   \sum\limits_{0\le i\le j \le k  \le n_3}
     \frac{x_{3,i}x_{3,j} ( x_{3,k} + x_{1,i}+ x_{2,j} ) }{(n_3)^2} \\
\st   &   x_3 \in \re^{n_3},
\end{array}\right.
\earay
\end{equation*}
where $x_{1,0}=x_{2,0}=x_{3,0}=1$, and $n_1 = n_2 = n_3$.
We implement Algorithm~\ref{ag:allKKTSDP} for the cases $n_1=n_2=n_3=2,3,4,5,6$.
The computational results are shown in the following table.
For all cases, we computed a NE successfully
and obtained that $x_1^*=x_2^*=x_3^*$ (up to round-off errors).
There is a unique NE for each case.
The computational results are reported in Table~\ref{tab:exm:uncon}.
\begin{table}[htbp]
\centering
\caption{The computational results for Example~\ref{ep:unconstrained}.}
\label{tab:exm:uncon}
\begin{tabular}{|l|c|c|r|}  \hline
$n_1$  & $x_1^* = x_2^* = x_3^*$	& $\omega^*$ & time \quad  \\ \hline
$2$ & $(-0.8410,\quad -0.7125)$  & $-8.8291\cdot 10^{-9}$ & 0.34  \\\hline
$3$ & $(-0.6743,-0.6157,-0.5236)$  & $-6.6507\cdot 10^{-9}$ & 1.58  \\\hline
$4$ & $\begin{array}{ll}    (-0.5950,-0.5606\\
   \qquad\qquad-0.5097,-0.4363)\end{array}$
	& $-1.0577\cdot 10^{-9}$ & 16.86  \\\hline
$5$ & $\begin{array}{l}    (-0.5476,-0.5247,-0.4919,\\\qquad\qquad-0.4472,-0.3860)\end{array}$
	& $-4.4438\cdot 10^{-9}$ & 177.63  \\\hline
$6$ & $\begin{array}{l}    (-0.5157,-0.4992,-0.4762,\\-0.4457,-0.4060,-0.3534)\end{array}$
	& $-3.7536\cdot 10^{-9}$ & 1379.27  \\\hline

\end{tabular}
\end{table}
The time is displayed in seconds.
\end{example}

The following are some examples of NEPs from applications.

\begin{example}\rm
Consider the environmental pollution control problem for three countries
for the case \textit{autarky} \cite{breton2006game}.

Let $ x_{i,1} (i= 1,2,3) $ denote the (gross) emissions
from the $i$th country.

The revenue of the $i$th country depends on $ x_{i,1}$,
e.g., a typically one is $x_{i,1}(b_i- \frac{1}{2} x_{i,1})$.
The variable $x_{i,2} $ represents the investment by the $i$th country
to local environmental projects.

The net emission in country $i$ is $x_{i,1}-\gamma_{i}x_{i,2}$,
which is always nonnegative and must be kept below or equal a certain prescribed level $E_i>0$
under an environmental constraint.
The damage cost of the $i$th country is assumed to be

$d_i( x_{i,1}-\gamma x_{i,2})+\sum_{j\ne i}c_{i,j}x_{i,2}x_{j,1}.$
For given parameters $b_i,c_{i,j},d_i,\gamma_i,E_i$,
the $i$th ($i=1,2,3$) country's optimization problem is
\[
\left\{ \begin{array}{cl}
 \min\limits_{ x_i \in \re^2} &-x_{i,1}(b_i-\frac{1}{2} x_{i,1})+\frac{(x_{i,2})^2}{2}
            +d_i( x_{i,1}-\gamma_i x_{i,2})+\sum\limits_{j\ne i}c_{i,j}x_{i,2}x_{j,1} \\

    \st &x_{i,2}\ge0,\  x_{i,1}\le b_i,\\
    & 0\le x_{i,1}-\gamma_i x_{i,2}\le E_i.
  \end{array} \right.
\]
We consider the general cases that $b_i\ne E_i$.
The Lagrange multipliers can be expressed as
\be \nn
\begin{array}{llr}
	\lambda_{i,4}=&\frac{1}{(b_i-E_i)E_i}\left(\frac{\partial f_i}{\partial x_{i,2}}x_{i,2}(x_{i,1}
         -\gamma_ix_{i,2})-\frac{\partial f_i}{\partial x_{i,1}}(b_i-x_{i,1})(x_{i,1}-\gamma_ix_{i,2})\right),\\
	\lambda_{i,3}=&\frac{1}{b_i}\left((b_i-x_{i,1})(\frac{\partial f_i}{\partial x_{i,1}}+\lambda_{i,4})-x_{i,2}(\frac{\partial f_i}{\partial x_{i,2}}-\gamma_i\lambda_{i,4}) \right),\\
	\lambda_{i,2}=&\lambda_{i,3}-\lambda_{i,4}-\frac{\partial f_i}{\partial x_{i,1}},\\
	\lambda_{i,1}=&\frac{\partial f_i}{\partial x_{i,2}}+\gamma_i\lambda_{i,3}-\gamma_i\lambda_{i,4}.
\end{array}
\ee
We solve the NEP for the following typical parameters:
\[
\begin{array}{llllll}
b_1=1.5, & b_2=2, & b_3=1.8, & c_{1,2}=0.2, & c_{1,3}=0.3, & c_{2,1}=0.4,\\
c_{2,3}=0.2, & c_{3,1}=0.5, & c_{3,2}=0.1, & d_1=0.8, & d_2=1.2, & d_3=1.0,\\
E_1=3,  & E_2=4, & E_3=2,  & \gamma_1=0.7, & \gamma_2=0.5, & \gamma_3=0.9.
\end{array}
\]
By Algorithm~\ref{ag:allKKTSDP}, we get the unique NE
\[
x_1^*=(0.7000,0.1600), \quad x_2^* = (0.8000,0.1600), \quad x_3^* = (0.8000,0.4700),
\]
with the accuracy parameter $-1.1059\cdot 10^{-9}$.
It took about $10$ seconds.
\end{example}

\begin{example}\rm
Consider the NEP of the electricity market problem
\cite{Contreras2004,FacKan10}.
There are three generating companies,
and the $i$th company possesses $s_i$ generating units.
For the $i$th company,
the power generation of his $j$th generating unit is denoted by $x_{i,j}$.
Assume $0\le x_{i,j}\le E_{i,j}$,
where the nonzero parameter $E_{i,j}$ represents its maximum capacity,
and the cost of this generating unit is $\frac{1}{2}c_{i,j}(x_{i,j})^2+d_{i,j}x_{i,j}$,
where $c_{i,j},d_{i,j}$ are parameters.
The electricty price is given by
$$\phi(\allx) \, := \, b-a(\sum_{i=1}^3\sum_{j=1}^{s_i}x_{i,j}).$$
The aim of the each company is to maximize its profits,
that is,
to solve the following optimization problem:
\[
\baray{l}
\mbox{$i$th player:} \left\{
 \begin{array}{cl}
 \min\limits_{x_i \in \re^{s_i}} & \frac{1}{2}\sum_{j=1}^{s_i}(c_{i,j}(x_{i,j})^2+d_{i,j}x_{i,j})-\phi(\allx)(\sum_{j=1}^{s_i}x_{i,j}).\\
\st  & 0\le x_{i,j}\le E_{i,j}\,(j\in[s_i]).
  \end{array}\right.
\earay
\]
The Lagrange multipliers according to the constraints
$g_{i,2j-1}:=E_{i,j}-x_{i,j}\ge0,g_{i,2j}:=x_{i,j}\ge0$
can be represented as
$$\lambda_{i,2j-1}=-\frac{\partial f_i}{\partial x_{i,j}}\cdot x_{i,j}/E_{i,j},\lambda_{i,2j}=\frac{\partial f_i}{\partial x_{i,j}}+\lambda_{i,2j-1}.\,(j\in [s_i])$$
We run Algorithm~\ref{ag:allKKTSDP} for the following setting:
\be \nn
\begin{array}{llllll}
 s_i=i,&&a=1,&& b=10,\\
 c_{1,1}=0.4, & c_{2,1}=0.35, & c_{2,2}=0.35, & c_{3,1}=0.46, & c_{3,2}=0.5, & c_{3,3}=0.5, \\
 d_{1,1}=2,   & d_{2,1}=1.75, & d_{2,2}=1,    & d_{3,1}=2.25, & d_{3,2}=3,   & d_{3,3}=3,  \\
 E_{1,1}=2,   & E_{2,1}=2.5,  & E_{2,2}=0.67, & E_{3,1}=1.2,  & E_{3,2}=1.8, & E_{3,3}=1.6 .
\end{array}
\ee
We found the unique NE
\[
x_1^* = 1.7184, \quad x_2^* = (1.8413,0.6700), \quad
x_3^* = (1.2000,0.0823,0.0823).
\]
The accuracy parameter is $-5.1183\cdot 10^{-7}$.
It took about $8$ seconds.
\end{example}

\section{Conclusions}
\label{sc:conc}

This paper studies Nash equilibrium problems that are given by polynomial functions.
Algorithms \ref{ag:KKTSDP} and \ref{ag:allKKTSDP} are proposed for
computing one or all NEs.
The Lasserre type Moment-SOS hierarchy of semidefinite relaxations
are used to solve the appearing polynomial optimization problems.
Under generic assumptions, we can compute a Nash equilibrium if it exists,
and detect the nonexistence if there is none.
Moreover, we can get all Nash equilibria
if there are finitely many ones of them.

\appendix
\appendixpage
\addappheadtotoc
\section{On the finiteness of KKT points for generic NEPPs}

We discuss the finiteness of KKT points for generic NEPPs.
This implies that Algorithm~\ref{ag:KKTSDP} and \ref{ag:allKKTSDP}
has finite convergence.
After enumeration of all possibilities of active inequality constraints,
we can generally consider the case that
\reff{eq:subprob} has only equality constraints.
Consequently, the length $m_i$ of the $i$th player's constraining polynomials
can be assumed less than or equal to $n_i$,
the dimension of its strategy $x_i$.
To prove the finiteness, we can ignore the sign conditions
$\lmd_{i,j} \ge 0$ for Lagrange multipliers.
The KKT system for all players is
\be
\label{KKTsym:allplayers}
\left\{
\begin{array}{rcll}
\sum_{j=1}^{m_i}\lambda_{ij}  \nabla_{x_i} g_{i,j}(x_i) &=&
        \nabla_{x_i} f_i(\allx) & (i \in [N]),  \\
 g_{i,j}(x_i) & = & 0 & (i \in [N], j \in  [m_i]) .
\end{array}
\right.
\ee
When the objectives $f_i$ are generic polynomials in $\allx$
and each $g_{i,j}$ is a generic polynomial in $x_i$,
we show that \reff{KKTsym:allplayers}
has finitely many complex solutions.

\begin{theorem} \label{tm:genfinite}
Let $d_{i,j}>0$, $a_{i,j} >0$ be degrees, for $i \in [N], j \in  [m_i]$.
If each $g_{i,j}$ is a generic polynomial in $x_i$ of degree $d_{i,j}$
and each $f_i$ is a generic polynomial in $\allx$
and its degree in $x_j$ is $a_{i,j}$, then the KKT system
\reff{KKTsym:allplayers} has finitely many complex solutions
and hence the NEP has finitely many KKT points.
\end{theorem}

\begin{proof}
For each player $i=1, \ldots, N$, denote
\[
b_i := a_{i,i}-1+d_{i,1} + \cdots + d_{i,m_i} - m_i.
\]
\[
\tilde{x}_i :=(x_{i,0},x_{i,1},\ldots,x_{i,n_i}), \quad
\tilde{x} := (\tilde{x}_1, \ldots,  \tilde{x}_N  ) .
\]
The homogenization of $g_{i,j}$
is $\tilde{g}_{i,j}$, a form in $\tilde{x}_i$.
Let $\mathbb{P}^{n_i}$ be the $n_i$ dimensional projective space, over the complex field.

Consider the projective varieties
\[
\mathcal{U}_i := \Big \{
(\tilde{x}_1,\dots,\tilde{x}_N)\in
\mathbb{P}^{n_1}\times\dots\times\mathbb{P}^{n_N}:
\tilde{g}_{i}(\tilde{x}_i)=0
\Big \}, i =1,\ldots, N,
\]
\[
\mathcal{U}:= \mc{U}_1 \cap \cdots \cap \mc{U}_N.
\]
When $g_{i,j}$ are generic polynomials in $x_i$,
the codimension of $\mc{U}_i$ is $m_i$ (see \cite{Harris1992}),
so $\mc{U}$ has the codimension $m_1 + \cdots + m_N$.

The $i$th player's objective $f_i$ is a polynomial in
$\allx=(x_1, \ldots, x_N)$,
we denote the multi-homogenization of $f_i(\xpi,\xmi)$ as
\[
\tilde{f}_i(\tilde{x}_i,\tilde{x}_{-i}):=
f_i(x_1/x_{1,0},\dots,x_N/x_{N,0})\cdot(\prod_{j=1}^N ( x_{i,0} )^{a_{i,j}}).
\]
It is a multi-homogeneous polynomial in $\tilde{\allx}$.
For each $i$, consider the determinantal variety
\[
W_i \,:=\, \left\{ \allx \in  \cpx^n\left|
\begin{array}{cccc}
\rank [ \nabla_{x_i} f_i(\allx) & \nabla_{x_i} g_{i,1}(x_i)
           & \cdots & \nabla_{x_i} g_{i,m_i}(x_i)] \le m_i
\end{array}
\right.\right\}.
\]
(The $\nabla_{x_i}$ denote the gradient with respect to $x_i$.)

Its multi-homogenization is
\[\widetilde{W}_i \,:=\, \left\{ \tilde{\allx}  \left|
\begin{array}{cccc}
  \rank [ \nabla_{x_i} \tilde{f}_i(\tilde{\allx}) & \nabla_{x_i} \tilde{g}_{i,1}(\tilde{x}_i)
              & \cdots   & \nabla_{x_i} \tilde{g}_{i,m_i} (\tilde{x}_i)] \le m_i
\end{array}
\right.\right\}.
\]
The matrix in the above can be explicitly written as
\[
J_i(\tilde{x}_i, \tilde{x}_{-i}) \, := \,
\left[\begin{array}{cccc}
\pt_{x_{i,1}} \tilde{f}_i(\tilde{\allx}) & \pt_{x_{i,1}} \tilde{g}_{i,1}(\tilde{x}_i)
              & \cdots   & \pt_{x_{i,1}} \tilde{g}_{i,m_i}(\tilde{x}_i)  \\
\pt_{x_{i,2}} \tilde{f}_i(\tilde{\allx}) & \pt_{x_{i,2}} \tilde{g}_{i,1}(\tilde{x}_i)
              & \cdots   & \pt_{x_{i,2}} \tilde{g}_{i,m_i} (\tilde{x}_i) \\
 \vdots &  \vdots & \ddots &  \vdots \\
\pt_{x_{i,n_i}} \tilde{f}_i(\tilde{\allx}) & \pt_{x_{i,n_i}} \tilde{g}_{i,1}(\tilde{x}_i)
              & \cdots   & \pt_{x_{i,n_i}} \tilde{g}_{i,m_i}(\tilde{x}_i)  \\
\end{array}
\right] .
\]
The $(m_i+1)$-by-$(m_i+1)$ minors of the matrix $J_i$
are homogeneous in $\tilde{x}_i$ of degree $b_i$.
They are homogeneous in $\tilde{x}_j$ of degree $a_{i,j}$, for $j \ne i$.
By \cite[Proposition~2.1]{nie2009},
when $g_{i,j}$ are generic polynomials in $x_i$,
the right $m_i$ columns of $J_i$ are linearly independent for all
$\tilde{x_i} \in \mc{U}_i$. That is, for every $\tilde{x} \in \mc{U}_i$,
there must exist a nonzero $m_i$-by-$m_i$ minor from the right
$m_i$ columns of $J_i$. In the following,
we consider fixed generic polynomials $g_{i,j}$.

First, we show that $\mc{U} \cap \widetilde{W}_1$
have the codimension $n_1 + m_2 + \cdots + m_N$.
Let $\mc{V}$ be the projective variety consisting of all
equivalent classes of the vectors
\be \label{proj:kron:all:xi}
\mathfrak{m}_1(\tilde{\allx}) :=
[\tilde{x}_1]_{b_1}^{hom} \otimes
[\tilde{x}_2]_{a_{1,2}}^{hom}  \otimes \cdots
\otimes [\tilde{x}_N]_{a_{1,N}}^{hom},
\ee
for equivalent classes of $\tilde{x} \in \mc{U}$.
In the above, $\otimes$ denotes the Kronecker product,
$[u]_d^{hom}$ denotes the vector of
all monomials in $u$ of degree equal to $d$.
In other words, $[u]_d^{hom}$ is the subvector of of $[u]_d$
for monomials of the highest degree $d$.
Note that $\mc{U}$ is birational to $\mc{V}$
(consider the natural embedding $\varphi:\mc{U}\hookrightarrow\mc{V}$ such that $\phi(\tilde{\allx})=\mathfrak{m}_1(\tilde{\allx})$).
So $\mc{U}$ and $\mc{V}$ have the same codimension \cite{Shafarevich2013}.
For each subset $I \subseteq [n_1]$ of cardinality $m_1$,
we use $\det_{I} J_1$ to denote the $m_1$-by-$m_1$ minor of $J_1$
for the submatrix whose row indices are in $I$
and whose columns are the right hand side $m_1$ columns. Then
\[
\widetilde{W}_1 \quad = \bigcup_{I \subseteq [n_1], |I| = m_i}  \mc{X}_{I}
\quad \mbox{where}
\]
\[
\mc{X}_{I} := \{ \tilde{\allx} :
\rank\, J_1(\allx) \le m_1, \, \mbox{det}_{I} J_1(\allx) \ne 0 \}.
\]
For each $I$, we have $\tilde{\allx} \in \mc{X}_{I}$
if and only if the $(m_1+1)$-by-$(m_1+1)$ minors of $J_1$,
corresponding to the row indices $I \cup \{\ell\}$ with $\ell \in [n_1] \backslash I$,
are equal to zeros. There are totally $n_1 - m_1$ such minors.
The vanishing of these $(m_1+1)$-by-$(m_1+1)$ minors of $J_1$
gives $n_1 - m_1$ linear equations in the vector
$\mathfrak{m}_1(\tilde{\allx})$ as in \reff{proj:kron:all:xi}.
The coefficients of these linear equations are linearly
parameterized by coefficients of $f_1$. Therefore,
when $f_1$ has generic coefficients, the set
\[
\mc{Y}_I :=
\{\mathfrak{m}_1(\tilde{\allx}) : \tilde{\allx} \in \mc{X}_{I}  \cap \mc{U}  \}
\]
is the intersection of $\mc{V}$ with $n_1 - m_1$
generic linear equations. Since $\mc{X}_{I}  \cap \mc{U}$
is birational to $\mc{Y}_I$, they have the same codimension,
so the codimension of $\mc{X}_{I}  \cap \mc{U}$
is $n_1 + m_2 + \cdots + m_N$.
This conclusion is true for all the above subsets $I$. Since
\[
\mc{U} \cap \widetilde{W}_1 \quad  =
\bigcup_{I \subseteq [n_1], |I| = m_1}  \mc{X}_{I} \cap \mc{U}  ,
\]
the codimension of $\mc{U} \cap \widetilde{W}_1$
is equal to $n_1 + m_2 + \cdots + m_N$.

Second, we can repeat the above argument to show that
\[
( \mc{U} \cap \widetilde{W}_1 ) \cap \widetilde{W}_2
\]
has codimension $n_1 + n_2 + m_3 + \cdots + m_N$.
Let $\mc{V}'$ be the projective variety
consisting of all equivalent classes of the vectors
\be
\mathfrak{m}_2(\tilde{\allx}) \, := \,
[\tilde{x}_1]_{a_{2,1}}^{hom} \otimes
[\tilde{x}_2]_{b_2}^{hom} \otimes
[\tilde{x}_3]_{a_{2,3}}^{hom} \otimes \cdots
\otimes [\tilde{x}_N]_{a_{2,N}}^{hom}
\ee
for equivalent classes of $\tilde{x} \in \mc{U} \cap \widetilde{W}_1$.
Note that $\mc{U} \cap \widetilde{W}_1$ is birational to $\mc{V}'$.
They have the same codimension. Similarly, we have
\[
\widetilde{W}_2 \quad = \bigcup_{I \subseteq [n_2], |I| = m_2}  \mc{X}'_{I}
\quad \mbox{where}
\]
\[
\mc{X}'_{I} \,:= \, \{ \tilde{\allx} :
\rank\, J_2(\allx) \le m_2, \, \mbox{det}_{I} J_2(\allx) \ne 0 \}.
\]
When $f_2$ has generic coefficients, the set
\[
\mc{Y}'_I :=
\{\mathfrak{m}_2(\tilde{\allx}) : \tilde{\allx} \in \mc{X}'_{I}
\cap \mc{U} \cap \widetilde{W}_1  \}
\]
is the intersection of $\mc{V}'$ with $n_2 - m_2$
generic hyperplanes. Since $\mc{X}'_{I}  \cap \mc{U} \cap \widetilde{W}_1$
is birational to $\mc{Y}'_I$, they have the same dimension,
so the codimension of $\mc{X}'_{I}  \cap \mc{U} \cap \widetilde{W}_1$
is $n_1 + n_2 + m_3 + \cdots + m_N$.
This conclusion is true for all $\mc{Y}'_I$. Since
\[
\mc{U} \cap \widetilde{W}_1 \cap \widetilde{W}_2  \quad =
\bigcup_{I \subseteq [n_2], |I| = m_2}  \mc{X}'_{I} \cap \mc{U} \cap \widetilde{W}_1,
\]
we know $\mc{U} \cap \widetilde{W}_1 \cap \widetilde{W}_2$
has the codimension $n_1 + n_2 + m_3 + \cdots + m_N$.

Similarly, by repeating the above, we can eventually show that
\[
\mc{U} \cap \widetilde{W}_1 \cap \widetilde{W}_2 \cap \cdots  \cap \widetilde{W}_N
\]
has codimension $n_1 + n_2  + \cdots + n_N$.
This implies the KKT system \reff{KKTsym:allplayers}
has codimension $n_1 + n_2 + \cdots + n_N$,
i.e., the dimension of the solution set of \reff{KKTsym:allplayers} is zero.
So, there are finitely many complex KKT points.
\end{proof}