# A Jacobi-type Newton method for Nash equilibrium problems with descent
guarantees

L. F. Bueno Institute of Science and Technology, Federal University of São
Paulo, São José dos Campos, SP, Brazil<EMAIL_ADDRESS>G. Haeser
Department of Applied Mathematics, University of São Paulo, São Paulo, SP,
Brazil<EMAIL_ADDRESS>O. Kolossoski Department of Applied Mathematics,
University of São Paulo, São Paulo, SP, Brazil<EMAIL_ADDRESS>
###### Abstract

A common strategy for solving an unconstrained two-player Nash equilibrium
problem with continuous variables is applying Newton’s method to the system
obtained by the corresponding first-order necessary optimality conditions.
However, when taking into account the game dynamics, it is not clear what is
the goal of each player when considering they are taking their current
decision following Newton’s iterates. In this paper we provide an
interpretation for Newton’s iterate as follows: instead of minimizing the
quadratic approximation of the objective functions parameterized by the other
player current decision (the Jacobi-type strategy), we show that the Newton
iterate follows this approach but with the objective function parameterized by
a prediction of the other player action. This interpretation allows us to
present a new Newtonian algorithm where a backtracking procedure is introduced
in order to guarantee that the computed Newtonian directions, for each player,
are descent directions for the corresponding parameterized functions. Thus,
besides favoring global convergence, our algorithm also favors true minimizers
instead of maximizers or saddle points, unlike the standard Newton method,
which does not consider the minimization structure of the problem in the non-
convex case. Thus, our method is more robust compared to other Jacobi-type
strategies or the pure Newtonian approach, which is corroborated by our
numerical experiments. We also present a proof of the well-definiteness of the
algorithm under some standard assumptions, together with a preliminary
analysis of its convergence properties taking into account the game dynamics.

Keywords: Nash equilibrium, Newtonian method, Jacobi-type methods, Non-convex
game, Game dynamics.

##  1 Introduction

The Nash Equilibrium Problem (NEP), introduced in [20], models an economic
game where each player aims at minimizing their own objective function, but
the decision taken by each player influences the payoff of the other players.
This is a fundamental problem in economics and social behavior theory that has
been studied extensively. See, for instance, [2], [19] and [1]. For surveys on
the subject, see [10] and [12]. In this paper we consider the unconstrained
two-player NEP given by

$\begin{array}[]{cl}\text{Minimize }&f_{1}(x_{1},x_{2}),\\\\[-1.0pt]
\leavevmode\resizebox{}{7.11317pt}{$x_{1}\in\mathbb{R}^{n_{1}}$}\\\\[4.0pt]
\end{array}$ (1)

and

$\begin{array}[]{cl}\text{Minimize }&f_{2}(x_{1},x_{2}),\\\\[-1.0pt]
\leavevmode\resizebox{}{7.11317pt}{$x_{2}\in\mathbb{R}^{n_{2}}$}\\\\[4.0pt]
\end{array}$ (2)

where $f_{1},f_{2}:\mathbb{R}^{n_{1}}\times\mathbb{R}^{n_{2}}\to\mathbb{R}$
are (possibly non-convex) twice continuously differentiable functions that
describe each player’s objective to be minimized, parameterized by the other
player’s decision. A (local) solution to the NEP (1-2) is a point
$(x_{1}^{*},x_{2}^{*})\in\mathbb{R}^{n_{1}}\times\mathbb{R}^{n_{2}}$ such that
$x_{1}^{*}$ is a (local) minimizer for (1) with $x_{2}=x_{2}^{*}$, and
$x_{2}^{*}$ is a (local) minimizer for (2) with $x_{1}=x_{1}^{*}$.

In order to solve the NEP (1-2), some studies deal with a reformulation of the
problem, such as a variational inequality problem or as a complementarity
problem; see for instance [18] and [9]. One problem with this approach is that
several solutions of the reformulated problem may not be solutions of the
associated equilibrium problem. Another possibility is to reformulate the NEP
as an optimization problem via the Nikaido-Isoda function (see [21]), which
transforms the NEP into a minimax problem. A penalty update scheme was also
proposed in the works of [13] and [11]. Those indirect approaches are somewhat
computationally demanding, as they require solving a nontrivial optimization
problem at each step of the algorithm. In order to avoid such heavy
computations, the method proposed in this paper aims at solving the system
(1-2) without rewriting it as an indirect problem.

One may also consider partitioning methods, namely Jacobi- and Gauss-Seidel-
type splitting methods (see [10] and [28]) that deal with problem (1-2) by
iterating the idea of fixing the variables of one player and taking a step
towards minimizing the problem of the other player. In particular, a standard
Jacobi strategy would fix preliminary decisions for each player, and then
iteratively update them simultaneously by solving a model for the
corresponding optimization problem parameterized by the fixed decision of the
other player. Alternatively, so-called best response methods fix the decision
of each player as the best possible response, see [6] and [7]. The advantage
of these approaches is that they are able to treat problem (1-2) in its
original form, without losing information due to reformulations. However,
these methods also require solving a general optimization problem at each step
for finding the best response for the other player before solving the
corresponding optimization problem for one player, which may also be
computationally expensive, in addition to having very restrictive convergence
conditions.

The minimization structure plays a key role when optimization problems are
non-convex. Several important applications result in such problems, including
NEPs (see [8]). Moreover, the dynamics of the algorithms, namely, the way that
the iterates are computed, may capture some elements of the true dynamics of
the game, reflecting the way a competitive game actually unfolds in practice.
That is, we suspect that economic systems behave in practice as a dynamical
system, where each player considers the current decision of the other players
and a model for predicting their future behavior, before optimizing their
strategy. We expect that Jacobi-type methods should capture this behavior,
predicting accurately the final outcome of a game, even if no equilibrium is
reached. In fact, differently from what is usual in an algorithm for solving a
single optimization problem, where we expect and hope to build a sequence
converging to a solution of the problem, in a NEP viewed as such dynamical
system, the dynamics may converge to an equilibrium state, but other
interesting phenomena may be observed such as orbiting, cycling or divergence.
The reader should have this observation in mind when considering what type of
theoretical property should be expected by an algorithm aiming to capture such
dynamics.

In this paper we introduce a method based on a Jacobi-type approach within a
Newtonian framework. The main novelty of our approach is that instead of
taking a decision based on the current decision of the opponent, our Newton
model is built around a prediction of the opponent step. This should reflect
the reality of some equilibrium games in practice, where each player tries to
make the best decision based on what is the expected response of the other
player. The algorithm proposed in this paper is posed in such a way that each
player aims at decreasing their objective function, with respect to the
predicted decision of the other player, using quadratic approximations. This
avoids the heavy computations of best response methods, while possibly still
reflecting the dynamics of equilibrium games in practice.

Additionally, we note that our method considers the true minimization
structure of the system (1-2) in a Newtonian framework, instead of considering
only the system of equations given by the necessary optimality conditions
associated with each problem. The idea is that one should favor minimizers in
place of maximizers or saddle points when considering the solutions of this
system of equations, which is usually not considered in other Newtonian
algorithms for solving non-convex instances of (1-2).

In a related work, a Jacobi-type approach using quadratic approximations of
players’ objective functions was considered in [28]. The proposed algorithm
consists in a trust region method which does not directly use the Hessian
information on its iterative step; however, it relies on some strong
assumptions on the true Hessians that we do not rely on in our work.
Additionally, for each player, the method in [28] is based on the current
response of the other players, rather than on a predicted one.

The paper is organized as follows: In Section 2 we introduce our Jacobi-type
algorithm in a Newtonian framework. In Section 3 we establish some notations
and summarize the overall assumptions to be considered throughout the paper,
together with some preliminary results. In Section 4 we establish the well-
definiteness of our algorithm. In Section 5 we state some theoretical
properties of the algorithm, establishing convergence to stationarity under
certain conditions. In section 6 we illustrate the behavior and effectiveness
of the proposed method with some numerical experiments. Finally, some
concluding remarks and future prospective works are given in Section 7.

Notation: Given a twice continuously differentiable function
$f:\mathbb{R}^{n_{1}+n_{2}}\to\mathbb{R},$$(x_{1},x_{2})\mapsto
f(x_{1},x_{2})$, for $i=1,2$ we denote by $\nabla_{x_{i}}f(\overline{x})$ its
partial gradient with respect to variable $x_{i}$ evaluated at $\overline{x}$.
The partial Hessian with respect to variables $x_{i}$ and $x_{j}$ evaluated at
$\overline{x}$ is denoted by $\nabla_{x_{i}x_{j}}^{2}f(\overline{x})$, for
$i,j=1,2$. We use $\|\cdot\|$ to denote the $2$-norm of vectors and matrices,
while $u^{T}v$ denotes the canonical inner product of vectors $u$ and $v$ of
the same dimensions. We denote by $\mathbb{N}$ the set of positive integers.

##  2 A Descent Newton Algorithm for the Two-Player Nash Equilibrium Problem

Recall that a solution $(x_{1},x_{2})$ to the equilibrium problem (1-2)
satisfies the following first-order necessary optimality condition:

$\left(\begin{array}[]{c}\nabla_{x_{1}}f_{1}(x_{1},x_{2})\\\\[3.0pt]
\nabla_{x_{2}}f_{2}(x_{1},x_{2})\end{array}\right)=\left(\begin{array}[]{c}0\\\\[3.0pt]
0\end{array}\right).$ (3)

Most algorithms for solving the NEP (1-2) rely on approaching directly the
nonlinear system of equations (3), disregarding the minimization structure of
problem (1-2). This poses no issues when $f_{1}$ and $f_{2}$ are convex with
respect to their decision variables, however, in the presence of non-
convexities, the algorithm may find local maximizers or saddle points as often
as local minimizers. Our goal in this paper is to propose an algorithm that
takes into account the minimization structure of the problem in order to find
local minimizers more often, thus better reflecting the behavior of two-player
games in practice.

Given a current iterate $(x_{1}^{k},x_{2}^{k})$, the classical Newtonian step
$(d_{1}^{k},d_{2}^{k})$ for solving (3) is given by solving the following
linear system:

$\left(\begin{array}[]{cc}\nabla^{2}_{x_{1}x_{1}}f_{1}(x_{1}^{k},x_{2}^{k})&\nabla^{2}_{x_{1}x_{2}}f_{1}(x_{1}^{k},x_{2}^{k})\\\\[3.0pt]
\nabla^{2}_{x_{1}x_{2}}f_{2}(x_{1}^{k},x_{2}^{k})&\nabla^{2}_{x_{2}x_{2}}f_{2}(x_{1}^{k},x_{2}^{k})\end{array}\right)\left(\begin{array}[]{c}d_{1}^{k}\\\
d_{2}^{k}\end{array}\right)=-\left(\begin{array}[]{c}\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k})\\\\[3.0pt]
\nabla_{x_{2}}f_{2}(x_{1}^{k},x_{2}^{k})\end{array}\right).$ (4)

It is well known that, for optimization problems, applying Newton’s method to
the corresponding first-order necessary optimality conditions can be better
interpreted as minimizing a quadratic approximation of the objective function.
With this interpretation, we can extend Newton’s method to non-convex
problems, and we can also provide elements for globalization schemes. Two
well-stablished strategies for doing this are the use of trust regions and
line search techniques, ensuring progress near the current point. The first
one guarantees the existence of a minimizer for the quadratic approximation of
the objective function by establishing a compact domain for the subproblem. In
line search approaches, one should use some positive definite approximation
for the Hessian of the objective function to obtain a coercive model and
descent directions. In this case, one must carefully choose a convenient step-
size along the Newton direction.

In order to exploit the minimization structure of the NEP (1-2), we shall
consider an idea which resembles a standard line-search globalization
described above. Given an approximate solution $(x_{1}^{k},x_{2}^{k})$ at some
iteration $k$, the most natural step to be made by the first player would be
associated with minimizing, with respect to variable $d_{1}$, a simplified
model for the function $f_{1}(x_{1}^{k}+d_{1},x_{2}^{k})$. Namely, the first
player would select its response by computing

$\begin{array}[]{cl}\text{Minimize}&\phi_{1}(x_{1}^{k}+d_{1},x_{2}^{k}),\\\\[-5.0pt]
\leavevmode\resizebox{}{5.12128pt}{$d_{1}$}\\\\[2.0pt] \end{array}$ (5)

where $\phi_{1}$ is an approximation of $f_{1}$.

Similarly, the second player would seek to minimize at the same time a model
for the function $f_{2}(x_{1}^{k},x_{2}^{k}+d_{2})$ with respect to $d_{2}$,
leading to

$\begin{array}[]{cl}\text{Minimize}&\phi_{2}(x_{1}^{k},x_{2}^{k}+d_{2}).\\\\[-5.0pt]
\leavevmode\resizebox{}{5.12128pt}{$d_{2}$}\\\\[2.0pt] \end{array}$ (6)

This approach is related to Jacobi’s method and its theoretical investigation
has been conducted in [28] under a trust-region globalization. The main
concern that we want to address is to base the minimization model for each
player on a prediction of the other player’s action, rather than on its actual
decision. Thus, the main contribution of this work is considering the case
when the direction $d_{1}^{k}$, the decision of the first player, is obtained
by a Newton step with respect to a predicted behavior of the second player.
Namely, given $x_{2}^{k}+d_{2}^{k}$, a predicted decision for player two, the
problem we aim to solve for player one is:

$\begin{array}[]{cl}\text{Minimize}&\phi_{1}(x_{1}^{k}+d_{1},x_{2}^{k}+d_{2}^{k}).\\\\[-5.0pt]
\leavevmode\resizebox{}{5.12128pt}{$d_{1}$}\\\\[2.0pt] \end{array}$ (7)

Inspired by a Newtonian appeal, let us consider the following convex quadratic
approximation:

$f_{1}(x_{1}^{k}+d_{1},x_{2}^{k}+d_{2}^{k})\approx
f_{1}(x_{1}^{k},x_{2}^{k}+d_{2}^{k})+\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+d_{2}^{k})^{T}d_{1}+\frac{1}{2}d_{1}^{T}H_{1}^{k}d_{1},$
(8)

where $H_{1}^{k}$ is a positive definite approximation for the Hessian
$\nabla^{2}_{x_{1}x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+d_{2}^{k})$. Since
$d_{2}^{k}$ needs to be found together with $d_{1}^{k}$, suppose that the same
approach is done simultaneously by the other player in order to make their
decision $d_{2}^{k}$. Then, given $(x_{1}^{k},x_{2}^{k})$, we would still have
a Nash equilibrium problem in the variables $(d_{1},d_{2})$. However, the
resulting NEP is simpler to solve, since the objective functions are
parameterized convex quadratics in the corresponding decision variables. But,
since the term $\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+d_{2}^{k})^{T}d_{1}$
might still combine the variables in a non-linear way, the approximated
problem might still be hard to solve. Thinking in eliminating this
inconvenience with a typical idea from Newton’s method, we use linear
approximations of the gradients in the following way:

$\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+d_{2}^{k})\approx\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k})+\nabla^{2}_{x_{1}x_{2}}f_{1}(x_{1}^{k},x_{2}^{k})d_{2}^{k}.$
(9)

In this way, our approximated problem (7) is now clearly defined by combining
approximations (8) and (9), which has a simple structure and its solution is
readily given by solving the linear system:

$H_{1}^{k}d_{1}+\nabla^{2}_{x_{1}x_{2}}f_{1}(x_{1}^{k},x_{2}^{k})d_{2}^{k}+\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k})=0.$

Repeating the approach for the second player, with a positive definite
$H_{2}^{k}\approx\nabla^{2}_{x_{2}x_{2}}f_{2}(x_{1}^{k}+d_{1}^{k},x_{2}^{k})$,
our iteration is based on the solution of the following linear system:

$\left(\begin{array}[]{cc}H_{1}^{k}&\nabla^{2}_{x_{1}x_{2}}f_{1}(x_{1}^{k},x_{2}^{k})\\\\[3.0pt]
\nabla^{2}_{x_{1}x_{2}}f_{2}(x_{1}^{k},x_{2}^{k})&H_{2}^{k}\end{array}\right)\left(\begin{array}[]{c}d_{1}\\\
d_{2}\end{array}\right)=-\left(\begin{array}[]{c}\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k})\\\\[3.0pt]
\nabla_{x_{2}}f_{2}(x_{1}^{k},x_{2}^{k})\end{array}\right).$ (10)

Given the equivalence with the standard Newtonian system (4), it is natural to
choose $H_{1}^{k}\approx\nabla^{2}_{x_{1}x_{1}}f_{1}(x_{1}^{k},x_{2}^{k})$ and
$H_{2}^{k}\approx\nabla^{2}_{x_{2}x_{2}}f_{2}(x_{1}^{k},x_{2}^{k})$. Thus, we
have associated the standard Newtonian system (4) with the original NEP, where
every player makes its decision considering the predicted decision of the
other player. It now seems natural to consider descent conditions associated
with the predicted functions $f_{1}(x_{1},x_{2}^{k}+d_{2}^{k})$ and
$f_{2}(x_{1}^{k}+d_{1}^{k},x_{2})$. However, we cannot guarantee that the
resulting directions from (10) are descent directions for these functions, as
is shown in the following example.

###### Example 2.1

Consider the equilibrium problem where
$f_{1}(x_{1},x_{2}):=\frac{1}{2}x_{1}^{2}+(x_{2}^{2}+2x_{2}+1)x_{1}$ and
$f_{2}(x_{1},x_{2}):=\frac{1}{2}(x_{2}+2)^{2}$. In this case, we have

  * •

$\nabla_{x_{1}}f_{1}(x_{1},x_{2})=x_{1}+x_{2}^{2}+2x_{2}+1$,
$\nabla^{2}_{x_{1}x_{1}}f_{1}(x_{1},x_{2})=1$,
$\nabla^{2}_{x_{1}x_{2}}f_{1}(x_{1},x_{2})=2x_{2}+2$

  * •

$\nabla_{x_{2}}f_{2}(x_{1},x_{2})=x_{2}+2$,
$\nabla^{2}_{x_{2}x_{2}}f_{2}(x_{2},x_{2})=1$ and
$\nabla^{2}_{x_{1}x_{2}}f_{2}(x_{2},x_{2})=0$.

Taking $(x_{1}^{k},x_{2}^{k}):=(0,0)$, the Newton direction
$(d_{1}^{k},d_{2}^{k})$ obtained from (10) is the solution of:

$\left(\begin{array}[]{cc}1&2\\\\[3.0pt]
0&1\end{array}\right)\left(\begin{array}[]{c}d_{1}\\\
d_{2}\end{array}\right)=-\left(\begin{array}[]{c}1\\\\[3.0pt]
2\end{array}\right).$

So $d_{1}^{k}=3$ and $d_{2}^{k}=-2$. Note that
$\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k})=\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+d_{2}^{k})=1$
and therefore $d_{1}^{k}$ is not a descent direction for
$f_{1}(x_{1},x_{2}^{k}+d_{2}^{k})$ (nor for $f_{1}(x_{1},x_{2}^{k})$) at
$x_{1}^{k}$.

This means that we do not expect the descent condition
$f_{1}(x_{1}^{k}+td_{1}^{k},x_{2}^{k}+d_{2}^{k})<f_{1}(x_{1}^{k},x_{2}^{k}+d_{2}^{k})$
to hold even for arbitrarily small $t$. Therefore, an algorithm that requires
an Armijo-like line search could be not well-defined. The reason for not
obtaining a descent direction when minimizing the convex quadratic
approximation of the function comes from the fact that we did not use its true
gradient. That is, the approximation of
$\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+d_{2}^{k})$ in (9) might not be
accurate enough when $\|d_{2}^{k}\|$ is large. To improve this approximation,
we should consider also taking small steps along the direction $d_{2}^{k}$. In
this way, the line search should be done simultaneously for $f_{1}$ and
$f_{2}$, considering a simultaneous backtracking in $d_{1}^{k}$ and
$d_{2}^{k}$. Thus, replacing $d_{2}^{k}$ by $td_{2}^{k}$ in the previous
derivation of (10), we arrive at the following system of equations for
computing $(d_{1},d_{2})$, given a tentative step length $t$:

$\left(\begin{array}[]{cc}H_{1}^{k}&t\nabla^{2}_{x_{1}x_{2}}f_{1}(x_{1}^{k},x_{2}^{k})\\\\[3.0pt]
t\nabla^{2}_{x_{1}x_{2}}f_{2}(x_{1}^{k},x_{2}^{k})&H_{2}^{k}\end{array}\right)\left(\begin{array}[]{c}d_{1}\\\
d_{2}\end{array}\right)=-\left(\begin{array}[]{c}\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k})\\\\[3.0pt]
\nabla_{x_{2}}f_{2}(x_{1}^{k},x_{2}^{k})\end{array}\right).$

Note that for $t>0$ sufficiently small, we expect the solution $(d_{1},d_{2})$
to be such that $d_{1}$ is a descent direction for problem (1) while $d_{2}$
is a descent direction for problem (2) considering the predicted decision for
the other player. Hence, our method favors local minimizers for solving
problem (1-2). Notice, however, that the solution $(d_{1},d_{2})$ depends on
the choice of the step length $t$, thus, once a step length is rejected by the
descent condition, a new direction $(d_{1},d_{2})$ must be recomputed,
similarly to a trust region approach. Finally, once the descent condition is
met for some $t$, the new iterate is defined as
$(x_{1}^{k+1},x_{2}^{k+1}):=(x_{1}^{k},x_{2}^{k})+t(d_{1},d_{2})$.

Our algorithm selects the stepsizes based on a standard Armijo-type condition,
namely, checking whether the decrease along a component (with the other
component fixed at the predicted point) is proportional to what is predicted
by the first order approximation of the function. In mathematical terms, for
some $\alpha\in(0,1)$, we perform the check

$\displaystyle f_{1}(x_{1}^{k}+td_{1},x_{2}^{k}+td_{2})\leq
f_{1}(x_{1}^{k},x_{2}^{k}+td_{2})+\alpha
t\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+td_{2})^{T}d_{1},$ $\displaystyle
f_{2}(x_{1}^{k}+td_{1},x_{2}^{k}+td_{2})\leq
f_{2}(x_{1}^{k}+td_{1},x_{2}^{k})+\alpha
t\nabla_{x_{2}}f_{2}(x_{1}^{k}+td_{1},x_{2}^{k})^{T}d_{2},$

decreasing $t$ and recomputing the directions whenever one of these
inequalities do not hold.

Before formally presenting the algorithm, let us discuss some other conditions
needed in the line search procedure. Differently from optimization problems,
since the directions are not necessarily descent directions when $t$ is large,
we must ensure a negative slope over the search direction when computing the
step-size $t$. In doing so, we are lead to also avoiding directions that are
near orthogonal to the gradient of the predicted function. Additionally, in
order to avoid stagnation when far from stationary points, we also demand that
the direction is not too small compared with the corresponding gradient of the
predicted function.

Mathematically, for $\theta\in(0,1)$ and $\gamma>0$, our line search checks
also if

$\displaystyle\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+td_{2})^{T}d_{1}$
$\displaystyle\leq-\theta\|\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+td_{2})\|\cdot\|d_{1}\|,$
$\displaystyle\|d_{1}\|$
$\displaystyle\geq\gamma\|\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+td_{2})\|,$
(11)

with a similar requirement for the second player. Notice that, in
optimization, similar inequalities are standard requirements on the direction
for guaranteeing global convergence of a descent method using an Armijo line
search.

Guided by all these ideas, we formalize below the algorithm for computing a
sequence $\\{(x_{1}^{k},x_{2}^{k})\\}_{k\in\mathbb{N}}$ for solving the NEP
(1-2), where we note that the algorithm includes a safeguarding strategy that
replaces the mixed hessians with zeroes whenever stationarity for the
corresponding player has been reached and $t$ is small enough, which will be
explained later.

1 Jacobi-type descent Newton algorithm

Step 0.

Given $(x_{1}^{0},x_{2}^{0})\in\mathbb{R}^{n_{1}}\times\mathbb{R}^{n_{2}}$,
let $\alpha\in(0,1)$, $\tau\in(0,1]$, $\theta>0,\gamma>0$ and set $k:=0$.

Step 1.

Compute $g_{1}^{k}:=\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k})$ and
$g_{2}^{k}:=\nabla_{x_{2}}f_{2}(x_{1}^{k},x_{2}^{k})$, choose symmetric
positive definite matrices $H_{1}^{k}$ and $H_{2}^{k}$, and set $t:=1$.

Step 2.

Set

$\displaystyle
M_{1}:=\left\\{\begin{array}[]{cc}\nabla^{2}_{x_{1}x_{2}}f_{1}(x_{1}^{k},x_{2}^{k}),&\text{
if }\|g_{1}^{k}\|>0\text{ or }t>\tau,\\\ 0,&\text{
otherwise;}\end{array}\right.$ (12) $\displaystyle
M_{2}:=\left\\{\begin{array}[]{cc}\nabla^{2}_{x_{1}x_{2}}f_{2}(x_{1}^{k},x_{2}^{k}),&\text{
if }\|g_{2}^{k}\|>0\text{ or }t>\tau,\\\ 0,&\text{
otherwise.}\end{array}\right.$

Step 3.

Check if the matrix

$\left[\begin{array}[]{cc}H_{1}^{k}&tM_{1}\\\\[3.0pt]
tM_{2}&H_{2}^{k}\end{array}\right]$

is non-singular. If not, set $t:=t/2$ and repeat Step 3.

Step 4.

Find $(d_{1},d_{2})$ by solving the linear system

$\left(\begin{array}[]{cc}H_{1}^{k}&tM_{1}\\\\[3.0pt]
tM_{2}&H_{2}^{k}\end{array}\right)\left(\begin{array}[]{c}d_{1}\\\
d_{2}\end{array}\right)=-\left(\begin{array}[]{c}g_{1}^{k}\\\\[3.0pt]
g_{2}^{k}\end{array}\right).$ (13)

Step 5.

Check the inequalities

$\displaystyle f_{1}(x_{1}^{k}+td_{1},x_{2}^{k}+td_{2})\leq
f_{1}(x_{1}^{k},x_{2}^{k}+td_{2})+\alpha
t\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+td_{2})^{T}d_{1},$ (14)
$\displaystyle\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+td_{2})^{T}d_{1}\leq-\theta\|\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+td_{2})\|\cdot\|d_{1}\|,$
(15)
$\displaystyle\gamma\|\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+td_{2})\|\cdot\|g_{1}^{k}\|\leq\|d_{1}\|\cdot\|g_{1}^{k}\|;$
(16) $\displaystyle f_{2}(x_{1}^{k}+td_{1},x_{2}^{k}+td_{2})\leq
f_{2}(x_{1}^{k}+td_{1},x_{2}^{k})+\alpha
t\nabla_{x_{2}}f_{2}(x_{1}^{k}+td_{1},x_{2}^{k})^{T}d_{2},$ (17)
$\displaystyle\nabla_{x_{2}}f_{2}(x_{1}^{k}+td_{1},x_{2}^{k})^{T}d_{2}\leq-\theta\|\nabla_{x_{2}}f_{2}(x_{1}^{k}+td_{1},x_{2}^{k})\|\cdot\|d_{2}\|,$
(18)
$\displaystyle\gamma\|\nabla_{x_{2}}f_{2}(x_{1}^{k}+td_{1},x_{2}^{k})\|\cdot\|g_{2}^{k}\|\leq\|d_{2}\|\cdot\|g_{2}^{k}\|.$
(19)

If one of the inequalities (14)-(19) do not hold, set $t:=t/2$ and go to Step
2.

Step 6.

Set $t^{k}:=t,d_{1}^{k}:=d_{1},d_{2}^{k}:=d_{2}$ and update
$x_{1}^{k+1}:=x_{1}^{k}+t^{k}d_{1}^{k}$,
$x_{2}^{k+1}:=x_{2}^{k}+t^{k}d_{2}^{k}$, and $k:=k+1$. Go to Step 1.

Notice that for finding the direction $(d_{1},d_{2})$, we use the linear
system (13), which for $i=1,2$ replaces the mixed Hessian matrix
$t\nabla^{2}_{x_{1}x_{2}}f_{i}(x_{1}^{k},x_{2}^{k})$, used in the previous
discussion when deducing (10), with zeroes whenever $g_{i}^{k}=0$ and
$t\leq\tau$. This is done in order to force $d_{i}$ to be zero if stationarity
is already reached for player $i$ and the step size is sufficiently small. As
a consequence, a standard Newton step is set for the other player. The reason
for doing this is that in these situations it may not be possible to obtain a
function decrease for player $i$. We also replaced the aforementioned
condition (11), which controls the ratio between gradient and direction, by
(16). Both conditions are the same when $g_{1}^{k}\neq 0$, but using (16)
ensures that when stationarity is reached for the first player, inequality
(16) is automatically satisfied. A similar situation occurs for the second
player with respect to condition (19).

On the other hand, taking $d_{i}^{k}=0$ may not be the most adequate strategy,
even if $g_{i}^{k}=0$, since stationarity for player $i$ may be lost after a
step is made for the other player. Thus, this is done only as a safeguarding
procedure after the stepsize is at least as small as a threshold value
$\tau\in(0,1]$. Notice that when $\tau:=1$, the safeguarding procedure is
always activated whenever $g_{i}^{k}=0$, while when $\tau<1$, it is activated
only after rejecting the Newtonian step at least once. This is done in order
to try using the Newtonian direction as much as possible due to its fast local
properties.

##  3 General Assumptions

The assumptions stated in this section are supposed to hold throughout the
paper without specific mentioning. Before stating them, we establish some
notation that shall be used in the sequel. We denote

$\displaystyle g^{k}:=\left[\begin{array}[]{c}g_{1}^{k}\\\
g_{2}^{k}\end{array}\right]=\left[\begin{array}[]{c}\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k})\\\
\nabla_{x_{2}}f_{2}(x_{1}^{k},x_{2}^{k})\end{array}\right],\quad{\cal
H}_{t,k}:=\left[\begin{array}[]{cc}H_{1}^{k}&tM_{1}\\\
tM_{2}&H_{2}^{k}\end{array}\right],$ $\displaystyle
d^{k}:=\left[\begin{array}[]{c}d_{1}^{k}\\\ d_{2}^{k}\end{array}\right],\quad
d:=\left[\begin{array}[]{c}d_{1}\\\ d_{2}\end{array}\right],\quad\text{ and
}\quad x^{k}:=\left[\begin{array}[]{c}x_{1}^{k}\\\
x_{2}^{k}\end{array}\right].$

Additionally, for each player $i=1,2$, we denote the index of the other player
by $\neg i$. The general assumptions that we use in all our analysis are the
following:

Assumption A  There exists $L>0$ such that $\nabla_{x_{i}}f_{i}(\cdot,x_{\neg
i})$ are $L$-Lipschitz continuous for every fixed $x_{\neg i}$ and
$i\in\\{1,2\\}$. As a consequence, we have that

$f_{i}(x_{i}+td_{i},x_{\neg i})\leq
f_{i}(x)+t\nabla_{x_{i}}f_{i}(x)^{T}d_{i}+\frac{t^{2}L}{2}\|d_{i}\|^{2}$

for every $x=(x_{i},x_{\neg i})\in\mathbb{R}^{n_{i}+n_{\neg
i}},d_{i}\in\mathbb{R}^{n_{i}}$, and $t>0$.

Assumption B  There exists $C_{H}>0$ such that
$\|\nabla^{2}_{x_{1}x_{2}}f_{i}(x)\|\leq C_{H}$ for all
$x\in\mathbb{R}^{n_{1}+n_{2}}$ and $i\in\\{1,2\\}$.

Assumption C  There exists $C_{R}>0$ and functions $r_{1}$ and $r_{2}$ such
that for $i\in\\{1,2\\}$ and for all $x=(x_{i},x_{\neg
i})\in\mathbb{R}^{n_{i}+n_{\neg i}}$, all $d_{\neg i}\in\mathbb{R}^{n_{\neg
i}}$ and $t>0$ it holds

$\nabla_{x_{i}}f_{i}(x_{i},x_{\neg i}+td_{\neg
i})=\nabla_{x_{i}}f_{i}(x)+t\nabla_{x_{1}x_{2}}^{2}f_{i}(x)d_{\neg
i}+r_{i}(t,x,d_{\neg i}),$ (22)

and

$\|r_{i}(t,x,d_{\neg i})\|\leq C_{R}t^{2}\|d_{\neg i}\|^{2}.$

Assumption D  There exist $\lambda_{max}\geq\lambda_{min}>0$ such that for all
$i\in\\{1,2\\}$ and $k\in\mathbb{N}$ the eigenvalues of $H_{i}^{k}$ lie in
$[\lambda_{min},\lambda_{max}]$. As a consequence, for all
$x_{i}\in\mathbb{R}^{n_{i}}$, we have that

$\begin{split}\lambda_{min}\leq\|H_{i}^{k}\|\leq\lambda_{max},\quad\frac{1}{\lambda_{max}}\leq\|(H_{i}^{k})^{-1}\|\leq\frac{1}{\lambda_{min}},\quad\text{
and}\\\ \lambda_{min}\|x_{i}\|^{2}\leq
x_{i}^{T}H_{i}^{k}x_{i}\leq\lambda_{max}\|x_{i}\|^{2},\quad
i\in\\{1,2\\}.\end{split}$

We highlight that all Assumptions A to C are simple assumptions on the
smoothness and boundedness of the high order derivatives of the functions
$f_{1}$ and $f_{2}$, while Assumption D can be easily guaranteed when choosing
the matrices $H_{i}^{k}$.

We end this section with a few results that arise from Assumptions A-D. The
first result guarantees uniform boundedness of the matrices ${\cal H}_{t,k}$
and ${\cal H}^{-1}_{t,k}$ for small enough values of $t$.

###### Lemma 3.1

For every $k\in\mathbb{N}$ and $t\in(0,1]$, the matrix ${\cal H}_{t,k}$
satisfies

$\|{\cal
H}_{t,k}\|\leq\mu_{max}:=\sqrt{\lambda_{max}^{2}+4\lambda_{max}C_{H}+C_{H}^{2}}.$
(23)

Moreover, if $t\leq\frac{\lambda_{min}^{2}}{8\lambda_{max}C_{H}}$, then ${\cal
H}_{t,k}$ is non-singular and $\|{\cal H}_{t,k}^{-1}\|\leq\mu_{min}$, where
the constant $\mu_{min}$ is given by

$\mu_{min}:=\frac{\sqrt{2}}{\lambda_{min}}.$ (24)

Proof: By the definition of ${\cal H}_{t,k}$, for a vector
$z=(z_{1},z_{2})\in\mathbb{R}^{n_{1}}\times\mathbb{R}^{n_{2}}$, we have that

${\cal H}_{t,k}z=\left[\begin{array}[]{cc}H_{1}^{k}&tM_{1}\\\
tM_{2}&H_{2}^{k}\end{array}\right]\cdot\left[\begin{array}[]{c}z_{1}\\\
z_{2}\end{array}\right]=\left[\begin{array}[]{c}H_{1}^{k}z_{1}+tM_{1}z_{2}\\\
tM_{2}z_{1}+H_{2}^{k}z_{2}\end{array}\right].$

Therefore, we can write

$\displaystyle\|{\cal H}_{t,k}z\|^{2}$
$\displaystyle=(H_{1}^{k}z_{1}+tM_{1}z_{2})^{T}(H_{1}^{k}z_{1}+tM_{1}z_{2})+(tM_{2}z_{1}+H_{2}^{k}z_{2})^{T}(tM_{2}z_{1}+H_{2}^{k}z_{2})$
$\displaystyle=\|H_{1}^{k}z_{1}\|^{2}+2tz_{1}^{T}H_{1}^{k}M_{1}z_{2}+\|tM_{1}z_{2}\|^{2}+\|H_{2}^{k}z_{2}\|^{2}+2tz_{2}^{T}H_{2}^{k}M_{2}z_{1}$
$\displaystyle+\|tM_{2}z_{1}\|^{2}.$ (25)

Now, by the definition of $M_{i}$ on (12), we see that
$\|M_{i}\|\leq\|\nabla_{x_{i}x_{\neg i}}^{2}f_{i}(x^{k})\|,i=1,2$. Hence, by
the Cauchy-Schwarz inequality, we have that

$\displaystyle\|{\cal
H}_{t,k}z\|^{2}\leq\|H_{1}^{k}z_{1}\|^{2}+2t\|H_{1}^{k}\|\cdot\|\nabla_{x_{1}x_{2}}^{2}f_{1}(x^{k})\|\cdot\|z_{1}\|\cdot\|z_{2}\|$
$\displaystyle+t^{2}\|\nabla_{x_{1}x_{2}}^{2}f_{1}(x^{k})\|^{2}\cdot\|z_{2}\|^{2}+\|H_{2}^{k}z_{2}\|^{2}+2t\|H_{2}^{k}\|\cdot\|\nabla_{x_{2}x_{1}}^{2}f_{2}(x^{k})\|\cdot\|z_{1}\|\cdot\|z_{2}\|$
$\displaystyle+t^{2}\|\nabla_{x_{2}x_{1}}^{2}f_{2}(x^{k})\|^{2}\cdot\|z_{1}\|^{2}.$

Using Assumptions B and D, and the fact that $t\leq 1$, we have that

$\displaystyle\|{\cal
H}_{t,k}z\|^{2}\leq\lambda_{max}^{2}\|z_{1}\|^{2}+2t\lambda_{max}C_{H}\cdot\|z_{1}\|\cdot\|z_{2}\|+t^{2}C_{H}^{2}\cdot\|z_{2}\|^{2}$
$\displaystyle+\lambda_{max}^{2}\|z_{2}\|^{2}+2t\lambda_{max}C_{H}\cdot\|z_{1}\|\cdot\|z_{2}\|+t^{2}C_{H}^{2}\|z_{1}\|^{2}$
$\displaystyle\leq\lambda_{max}^{2}\left(\|z_{1}\|^{2}+\|z_{2}\|^{2}\right)+4\lambda_{max}C_{H}\cdot\|z_{1}\|\cdot\|z_{2}\|+C_{H}^{2}\left(\|z_{1}\|^{2}+\|z_{2}\|^{2}\right)$
$\displaystyle=\left(\lambda_{max}^{2}+C_{H}^{2}\right)\left(\|z_{1}\|^{2}+\|z_{2}\|^{2}\right)+4\lambda_{max}C_{H}\|z_{1}\|\cdot\|z_{2}\|.$

Hence, observing that $\|z_{i}\|\leq\|z\|$ for $i=1,2$, we have that

$\|{\cal
H}_{t,k}z\|^{2}\leq\left(\lambda_{max}^{2}+4\lambda_{max}C_{H}+C_{H}^{2}\right)\|z\|^{2},$

and thus (23) holds. To prove the bound on $\|{\cal H}_{t,k}\|$, recalling
(25), we have that

$\displaystyle\|{\cal H}_{t,k}z\|^{2}$
$\displaystyle=\|H_{1}^{k}z_{1}\|^{2}+2tz_{1}^{T}H_{1}^{k}M_{1}z_{2}+\|tM_{1}z_{2}\|^{2}+\|H_{2}^{k}z_{2}\|^{2}+2tz_{2}^{T}H_{2}^{k}M_{2}z_{1}$
$\displaystyle+\|tM_{2}z_{1}\|^{2}\geq\lambda_{min}^{2}\|z_{1}\|^{2}-2t\lambda_{max}\|M_{1}\|\cdot\|z_{1}\|\cdot\|z_{2}\|+\|tM_{1}z_{2}\|^{2}$
$\displaystyle+\lambda_{min}^{2}\|z_{2}\|^{2}-2t\lambda_{max}\|M_{2}\|\cdot\|z_{1}\|\cdot\|z_{2}\|+\|tM_{2}z_{1}\|^{2}.$

Since $\|tM_{i}z_{\neg i}\|^{2}\geq 0$, $\|M_{i}\|\leq C_{H}$,
$\|z\|^{2}=\|z_{1}\|^{2}+\|z_{2}\|^{2}$, and $\|z_{i}\|\leq\|z\|$, we have

$\displaystyle\|{\cal H}_{t,k}z\|^{2}$
$\displaystyle\geq\left(\lambda_{min}^{2}-4t\lambda_{max}C_{H}\right)\|z\|^{2}.$

Thus, for $t\leq\frac{\lambda_{min}^{2}}{8\lambda_{max}C_{H}}$, we have
$\|{\cal H}_{t,k}z\|^{2}\geq\frac{\lambda_{min}^{2}}{2}\|z\|^{2}$, so ${\cal
H}_{t,k}$ is non-singular and $\|{\cal
H}_{t,k}^{-1}\|\leq\frac{\sqrt{2}}{\lambda_{min}}$. $\;\;\;\Box$

Next we prove a technical result that shall be used later.

###### Lemma 3.2

For every $k\in\mathbb{N}$, if
$t\leq\frac{\lambda_{min}^{2}}{8\lambda_{max}C_{H}}$ and $d$ is as in Step 4
of Algorithm 2, then $\|d\|\leq
C_{k}:=\mu_{min}(\|g_{1}^{k}\|+\|g_{2}^{k}\|)$, where $\mu_{min}$ is given by
(24).

Proof: By the system equation (13) and using Lemma 3.1, we have that

$\|d\|=\|{\cal H}^{-1}_{t,k}g^{k}\|\leq\|{\cal
H}^{-1}_{t,k}\|\cdot(\|g_{1}^{k}\|+\|g_{2}^{k}\|)\leq\mu_{min}(\|g_{1}^{k}\|+\|g_{2}^{k}\|)=C_{k}.$
(26)

$\;\;\;\Box$
The next result shows that for each component $i=1,2$ and each step $k$, the
gradients $g_{i}^{k}$ are comparable with $d_{i}$ provided that $t$ is small
enough.

###### Lemma 3.3

For every $k\in\mathbb{N}$, $i=1,2$, if $t$ satisfies

$t\leq\min\left\\{\frac{\lambda_{min}^{2}}{8\lambda_{max}C_{H}},\frac{\|g_{i}^{k}\|}{2C_{H}C_{k}}\right\\},$
(27)

where $C_{k}$ is given as in Lemma 3.2, then

$\frac{1}{2}\|g_{i}^{k}\|\leq\|H_{i}^{k}d_{i}\|\leq\frac{3}{2}\|g_{i}^{k}\|.$
(28)

Proof: We shall prove the result for $i=1$, being the result analogous for
$i=2$. By the first equation of system (13), the definition of $M_{1}$,
Assumption B and Lemma 3.2, we have that

$\|H_{1}^{k}d_{1}\|=\|-tM_{1}d_{2}-g_{1}^{k}\|\leq tC_{H}C_{k}+\|g_{1}^{k}\|.$

Thus, from (27), we conclude the right-hand side inequality of (28).
Similarly, to prove the left-hand side inequality in (28), considering again
the first equation of system (13), we have
$\|g_{1}^{k}\|=\|-H_{1}^{k}d_{1}-tM_{1}d_{2}\|\leq\|H_{1}^{k}d_{1}\|+tC_{H}C_{k}.$
Thus, the result follows from (27). $\;\;\;\Box$

##  4 Well Definiteness of the Algorithm

This section is devoted to proving well-definiteness of Algorithm 2. As
previously mentioned, by the way Algorithm 2 is built, in cases where one of
the gradients is zero, say $\|g_{2}^{k}\|=0$, and the step-size $t$ has
sufficiently decreased ($t\leq\tau$), it turns out that the direction for the
second player is null while for the first player the direction is computed as
a standard Newton direction for optimization with Armijo line search, which is
known to be well-defined and satisfy conditions (14)-(16) with
$t~{}\leq~{}\frac{2(1-\alpha)\lambda_{min}^{2}}{L\lambda_{max}}$,
$\theta=\frac{\lambda_{min}}{\lambda_{max}}$ and
$\gamma=\frac{1}{\lambda_{max}}$ (see for instance the discussion in [3],
pages 29-36, and Lemma 2.20 in [14]). Thus, we focus on the case where
$\|g_{i}^{k}\|\neq 0,i=1,2$, and we show that the inequalities (14)-(19) are
satisfied for sufficiently small $t$, meaning that the stepsize does not
decrease indefinitely.

We begin with the gradient/direction ratio inequalities (16) and (19).
However, we show a stronger version bounding the directions by the gradient on
the points $(x_{i}^{k},x_{\neg i}^{k}+td_{\neg i})$ from both above and below.

###### Lemma 4.1

For every $k\in\mathbb{N}$ and $i=1,2$, if $t$ satisfies

$t\leq\min\left\\{\frac{\|g_{i}^{k}\|}{8C_{H}C_{k}},\frac{\lambda_{min}^{2}}{8\lambda_{max}C_{H}}\right\\},$
(29)

where $C_{k}$ is given as in Lemma 3.2, then

$\gamma\left\|\nabla_{x_{i}}f_{i}(x_{i}^{k},x_{\neg i}^{k}+td_{\neg
i})\right\|\leq\|d_{i}\|\leq\beta\left\|\nabla_{x_{i}}f_{i}(x_{i}^{k},x_{\neg
i}^{k}+td_{\neg i})\right\|,$ (30)

for every $\gamma$ and $\beta$ such that

$\gamma\leq\frac{2}{3\lambda_{max}},\quad\beta\geq\frac{2}{\lambda_{min}}.$
(31)

Proof: We prove the result only for $i=1$, as the proof for $i=2$ is similar.
By (13), we have that

$-g_{1}^{k}=H_{1}^{k}d_{1}+tM_{1}d_{2}.$ (32)

Also, we have that

$\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+td_{2})=g_{1}^{k}+t\left[\displaystyle\int_{0}^{1}\nabla_{x_{1}x_{2}}^{2}f_{1}(x_{1}^{k},x_{2}^{k}+\xi
td_{2})d\xi\right]d_{2}.$ (33)

Combining (32) and (33) we obtain

$\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+td_{2})=-H_{1}^{k}d_{1}+t\left[-M_{1}+\displaystyle\int_{0}^{1}\nabla_{x_{1}x_{2}}^{2}f_{1}(x_{1}^{k},x_{2}^{k}+\xi
td_{2})d\xi\right]d_{2}.$ (34)

By the triangular inequality, we have that

$\|\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+td_{2})\|\geq\Bigg{|}\|H_{1}^{k}d_{1}\|-t\bigg{\lVert}\left[-M_{1}+\displaystyle\int_{0}^{1}\nabla_{x_{1}x_{2}}^{2}f_{1}(x_{1}^{k},x_{2}^{k}+\xi
td_{2})d\xi\right]d_{2}\bigg{\rVert}\Bigg{|}$ (35)

and

$\|\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+td_{2})\|\leq\|H_{1}^{k}d_{1}\|+t\bigg{\lVert}\left[-M_{1}+\displaystyle\int_{0}^{1}\nabla_{x_{1}x_{2}}^{2}f_{1}(x_{1}^{k},x_{2}^{k}+\xi
td_{2})d\xi\right]d_{2}\bigg{\rVert}.$ (36)

Next, we prove that for $t$ sufficiently small it holds

$\frac{1}{2}\|H_{1}^{k}d_{1}\|\leq\|\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+td_{2})\|\leq\frac{3}{2}\|H_{1}^{k}d_{1}\|.$
(37)

Let us suppose that
$\left[-M_{1}+\displaystyle\int_{0}^{1}\nabla_{x_{1}x_{2}}^{2}f_{1}(x_{1}^{k},x_{2}^{k}+\xi
td_{2})d\xi\right]d_{2}\neq 0$, since otherwise this is trivial from (34).
Hence, if we select $t$ satisfying

$t\leq\frac{\|H_{1}^{k}d_{1}\|}{2\bigg{\lVert}\left[-M_{1}+\displaystyle\int_{0}^{1}\nabla_{x_{1}x_{2}}^{2}f_{1}(x_{1}^{k},x_{2}^{k}+\xi
td_{2})d\xi\right]d_{2}\bigg{\rVert}},$ (38)

we see from (35) that the left-hand side inequality in (37) holds and, by
(36), we obtain the right-hand side inequality. Now, using that
$\lambda_{min}\|d_{1}\|\leq\|H_{1}^{k}d_{1}\|\leq\lambda_{max}\|d_{1}\|$, (37)
implies that

$\frac{\lambda_{min}}{2}\|d_{1}\|\leq\|\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+td_{2})\|\leq\frac{3\lambda_{max}}{2}\|d_{1}\|,$

which implies the desired bounds in (30).

Finally, let us show that $t$ satisfying (29) implies the bound (38). By
Assumption B and Lemma 3.2, we see that

$\displaystyle\bigg{\lVert}\left[-M_{1}+\displaystyle\int_{0}^{1}\nabla_{x_{1}x_{2}}^{2}f_{1}(x_{1}^{k},x_{2}^{k}+\xi
td_{2})d\xi\right]d_{2}\bigg{\rVert}$
$\displaystyle\leq\left[\|M_{1}\|+\bigg{\lVert}\displaystyle\int_{0}^{1}\nabla_{x_{1}x_{2}}^{2}f_{1}(x_{1}^{k},x_{2}^{k}+\xi
td_{2})d\xi\bigg{\rVert}\right]\cdot\|d_{2}\|$
$\displaystyle\leq\left[\|M_{1}\|+\displaystyle\int_{0}^{1}\lVert\nabla_{x_{1}x_{2}}^{2}f_{1}(x_{1}^{k},x_{2}^{k}+\xi
td_{2})\rVert d\xi\right]\cdot\|d_{2}\|\leq 2C_{H}C_{k},$

therefore, we have that

$\frac{\|H_{1}^{k}d_{1}\|}{4C_{H}C_{k}}\leq\frac{\|H_{1}^{k}d_{1}\|}{2\bigg{\lVert}\left[-M_{1}+\displaystyle\int_{0}^{1}\nabla_{x_{1}x_{2}}^{2}f_{1}(x_{1}^{k},x_{2}^{k}+\xi
td_{2})d\xi\right]d_{2}\bigg{\rVert}},$

thus, in order to attain (38), it suffices to select

$t\leq\frac{\|H_{1}^{k}d_{1}\|}{4C_{H}C_{k}}.$ (39)

Now, notice that inequalities (28) guaranteed by Lemma 3.3 ensure (39) holds
with $t$ satisfying the hypothesis (29). $\;\;\;\Box$

The next result ensures the angle inequalities (15) and (18) hold for $t$
sufficiently small.

###### Lemma 4.2

For every $k\in\mathbb{N}$ and $i=1,2$, if $g_{i}^{k}\neq 0$, $t$ satisfies

$t\leq\min\left\\{\displaystyle\sqrt{\frac{\lambda_{min}\|g_{i}^{k}\|}{8\lambda_{max}C_{R}(C_{k})^{2}}},\frac{\|g_{i}^{k}\|}{8C_{H}C_{k}},\frac{\lambda_{min}^{2}}{8\lambda_{max}C_{H}}\right\\},$
(40)

where $C_{k}$ is given as Lemma 3.2, and $\theta$ satisfies

$\theta\leq\frac{\lambda_{min}}{4\lambda_{max}},$ (41)

then,

$\nabla_{x_{i}}f_{i}(x^{k}_{i},x^{k}_{\neg i}+td_{\neg
i})^{T}d_{i}\leq-\theta\|\nabla f_{i}(x^{k}_{i},x^{k}_{\neg i}+td_{\neg
i})\|\cdot\|d_{i}\|.$ (42)

Proof: Once again, we work the proof of the main claim for the first
component. For the second one, the proof follows similarly. Since we are
considering $g_{1}^{k}\neq 0$, we have that
$M_{1}=\nabla_{x_{1}x_{2}}^{2}f_{1}(x^{k})$. Hence, Algorithm 2 yields

$H_{1}^{k}d_{1}=-g_{1}^{k}-t\nabla_{x_{1}x_{2}}^{2}f_{1}(x^{k})d_{2}.$ (43)

By Assumption C, we have that

$\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+td_{2})=g_{1}^{k}+t\nabla^{2}_{x_{1}x_{2}}f_{1}(x^{k})d_{2}+r_{1}(t,x^{k},d_{2}),$
(44)

where the remainder $r_{1}(t,x^{k},d_{2})$ is such that

$\|r_{1}(t,x^{k},d_{2})\|\leq C_{R}t^{2}\|d_{2}\|^{2}.$ (45)

Therefore, adding (44) to (43), we have

$H_{1}^{k}d_{1}=-\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+td_{2})+r_{1}(t,x^{k},d_{2}).$
(46)

Hence, since $\|H_{1}^{k}d_{1}\|\leq\lambda_{max}\|d_{1}\|$, we have that

$\displaystyle\|d_{1}\|\geq\frac{1}{\lambda_{max}}\|-\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+td_{2})+r_{1}(t,x^{k},d_{2})\|,$

and using the triangular inequality at the right-hand side we obtain

$\|d_{1}\|\geq\frac{1}{\lambda_{max}}\Big{|}\|r_{1}(t,x^{k},d_{2})\|-\|\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+td_{2})\|\Big{|}.$
(47)

On the other hand, multiplying (46) by $d_{1}^{T}$ and using the Cauchy-
Schwarz inequality, Assumptions C and Lemma 3.2, we have that

$\displaystyle d_{1}^{T}\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+td_{2})$
$\displaystyle\leq-
d_{1}^{T}H_{1}^{k}d_{1}+\|r_{1}(t,x^{k},d_{2})\|\cdot\|d_{1}\|$
$\displaystyle\leq-\lambda_{min}\|d_{1}\|^{2}+C_{R}C_{k}^{2}t^{2}\|d_{1}\|.$
(48)

Therefore, if we choose $t$ in (4) satisfying

$t\leq\displaystyle\sqrt{\frac{\lambda_{min}\|d_{1}\|}{2C_{R}C_{k}^{2}}},$

we obtain

$d_{1}^{T}\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+td_{2})\leq\frac{-\lambda_{min}\|d_{1}\|^{2}}{2}.$
(49)

Thus, we can combine (47) and (49), to write

$d_{1}^{T}\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+td_{2})\leq-\frac{\lambda_{min}}{2\lambda_{max}}\|d_{1}\|\Big{|}\|\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+td_{2})\|-\|r_{1}(t,x^{k},d_{2})\|\Big{|}.$
(50)

Next, to bound the right-hand side of (50), for $t$ satisfying (40) we have by
(28) that

$\displaystyle
t^{2}\leq\frac{\lambda_{min}\|g_{1}^{k}\|}{8\lambda_{max}C_{R}C_{k}^{2}}\leq\frac{\lambda_{min}\|H_{1}^{k}d_{1}\|}{4\lambda_{max}C_{R}C_{k}^{2}}\leq\frac{\lambda_{min}\|d_{1}\|}{4C_{R}C_{k}^{2}}.$
(51)

Hence, using (51) on (45), we have that

$\|r_{1}(t,x^{k},d_{2})\|\leq C_{R}t^{2}\|d_{2}\|^{2}\leq
C_{R}C_{k}^{2}t^{2}\leq\frac{\lambda_{min}\|d_{1}\|}{4}\leq\frac{\|\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+td_{2})\|}{2},$

where we used Lemma 4.1 in the last inequality. Therefore, we have that

$\displaystyle\Big{|}\|\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+td_{2})\|-\|r_{1}(t,x^{k},d_{2})\|\Big{|}$
$\displaystyle=\|\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+td_{2})\|-\|r_{1}(t,x^{k},d_{2})\|$
$\displaystyle\geq\frac{\|\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+td_{2})\|}{2}.$

Then, (50) gives us the desired inequality

$d_{1}^{T}\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+td_{2})\leq-\frac{\lambda_{min}}{4\lambda_{max}}\|d_{1}\|\cdot\|\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+td_{2})\|.$

$\;\;\;\Box$

Now we are ready to finish establishing the well-definiteness of Algorithm 2,
where we prove the remaining relations (14) and (17).

###### Proposition 4.1

For $i=1,2$, suppose that $g_{i}^{k}\neq 0$ for some $k\in\mathbb{N}$, and
assume that the stepsize $t$ satisfies

$t\leq\min\left\\{\displaystyle\sqrt{\frac{\lambda_{min}\|g_{i}^{k}\|}{8\lambda_{max}C_{R}C_{k}^{2}}},\frac{\|g_{i}^{k}\|}{8C_{H}C_{k}},\frac{\lambda_{min}^{2}}{8\lambda_{max}C_{H}},\frac{(1-\alpha)2\theta}{\beta
L}\right\\},$ (52)

where $\beta$ and $\theta$ are constants satisfying (31) and (41)
respectively. Then all inequalities (14)-(19) are satisfied. As a consequence,
iteration $k$ of Algorithm 2 is well-defined.

Proof: By Lemmas 4.1 and 4.2, it remains to prove (14), since (17) is
analogous. Fix $k\in\mathbb{N}$ and $i=1$. By Assumption A, we have that

$f_{1}(x_{1}^{k}+td_{1},x_{2}^{k}+td_{2})-f_{1}(x_{1}^{k},x_{2}^{k}+td_{2})\leq
t\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+td_{2})^{T}d_{1}+\frac{t^{2}L_{1}}{2}\|d_{1}\|^{2}.$

If $d_{1}=0$, then the result is immediate. If $d_{1}\neq 0$ and $t>0$, we
have that

$\displaystyle
t\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+td_{2})^{T}d_{1}+\frac{t^{2}L}{2}\|d_{1}\|^{2}\leq\alpha
t\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+td_{2})^{T}d_{1}$
$\displaystyle\Leftrightarrow
t\leq\frac{2(1-\alpha)\left[-\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+td_{2})^{T}d_{1}\right]}{L\|d_{1}\|^{2}},$

which can be guaranteed by

$t\leq\frac{2(1-\alpha)\theta\|\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+td_{2})\|}{L\|d_{1}\|}$

due to (42). Also, by Lemma 4.1, we have that
$\frac{\|\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k}+td_{2})\|}{\|d_{1}\|}\geq\frac{1}{\beta}$,
and thus, it suffices to select $t$ satisfying (52) to ensure (14).
$\;\;\;\Box$

Notice that our well-definiteness results assume that $t$ is smaller than a
positive constant which depends on the gradients $g_{i}^{k}$ at iteration $k$.
This is the technical reason why we must switch to the standard Newton method
for optimization whenever one player has already reached stationarity. We
could not prove or find a counter-example for the well-definiteness of the
method that never switches to standard Newton, that is, considering $M_{i}$
always equals to $\nabla^{2}_{x_{1}x_{2}}f_{i}(x^{k}),i=1,2,$ in place of
(12). With the use of this safeguarding procedure we have that any iteration
of the algorithm is well defined.

We end this section with a particular result when the iterates all lie on a
bounded set. The result states we can infer the boundedness of the stepsize
sequence from the boundedness of both gradients from below. This guarantees
that in this situation the line search criteria will not get increasingly hard
to be satisfied, forcing small steps and stalling the computation of the
iterates, whenever the decisions of both players are far from stationarity.

###### Proposition 4.2

Suppose that the sequence $\\{x^{k}\\}$ generated by Algorithm 2 lie on a
bounded set. If there exists a constant $c>0$ and $k_{0}\in\mathbb{N}$ such
that $\|g_{i}^{k}\|>c$ for each component $i=1,2$ and all $k>k_{0}$, then the
stepsize sequence $\\{t^{k}\\}$ is bounded away from zero.

Proof: Let $0<\bar{t}:=\min\\{t^{k}\\}$ for $k\leq k_{0}$. For $k>k_{0}$,
since $\\{x^{k}\\}$ is bounded, there exists an universal constant $C_{d}$
such that $0<C_{k}\leq C_{d}$, where $C_{k}$ is defined in Lemma 3.2.
Moreover, using that $\|g_{i}^{k}\|>c$, Proposition 4.1 ensures that if

$t\leq
t^{\prime}:=\min\left\\{\displaystyle\sqrt{\frac{\lambda_{min}c}{8\lambda_{max}C_{R}C_{d}^{2}}},\frac{c}{8C_{H}C_{d}},\frac{\lambda_{min}^{2}}{8\lambda_{max}C_{H}},\frac{(1-\alpha)2\theta}{\beta
L}\right\\},$

then conditions (14)-(19) hold. Thus, by the stepsize update rule we have
that, for all $k$,

$t^{k}\geq\min\left\\{1,\frac{t^{\prime}}{2},\bar{t}\right\\}>0.$

$\;\;\;\Box$

##  5 Convergence Analysis

In this section we establish our convergence results for Algorithm 2. Before
starting, we point out that equilibrium points are key elements in the study
of a game, but they certainly do not capture everything one may need in order
to understand the modeled problem. Unlike optimization problems, where the
solution represents the goal for the decision maker, NEP solutions may not be
desirable by the players. This happens, for example, in the famous Prisoner’s
Dilemma and in congestion problems (see [24] and [25]). In this way, an
algorithm that captures the dynamics involved in the game, even if it does not
obtain convergence, can sometimes be even more important than the
identification of the solution itself. A more concrete example comes from the
so-called evolutionary game theory (see [27]), which is used to explain the
basis of evolutionary behavior. This field differs from standard game theory
in the sense that it focus on how a change of strategy affects the dynamics of
the game, rather than finding equilibrium points, since a change of strategies
describes competition between populations. In arbitrary dynamical systems, one
can obtain stable equilibrium points, where trajectories in their neighborhood
converge to them. However, it is known that one can also observe trajectories
that may orbit or cycle, or that may simply diverge from an unstable
equilibrium point or even have a chaotic behavior, see [23]. Since our
algorithm in some sense tries to emulate the game dynamics, we believe that
all these situations may occur with the sequences it generates. Therefore, it
does not seem reasonable to show convergence to equilibrium points in general
problems.

On the other hand, Nash equilibrium situations are observed in many real-world
applications. Therefore, it is natural to imagine that the dynamics involved
in these applications favors the stability of the system. Moreover, when
convergence occurs, the limit point must be an attractor and therefore it is
expected to be an equilibrium point. In other words, we believe that the
process of evolution of the generated sequence should privilege convergence to
equilibrium points, even if they are not beneficial for the agents. However,
in situations where conditions change very quickly it may be hard to achieve
an equilibrium, because, in general, players tend to respond to outdated
information. Traditional best-response algorithms suffer from this problem and
therefore need very restrictive assumptions in order to achieve convergence.
Our method is constructed on basing the iteration on a predicted response
rather than on the current one, and therefore generates more stability in the
process. This type of approach may be more related to the real dynamics of the
modeled problem. Considering this, we expect to obtain equilibrium points more
frequently than other Jacobi-type methods in the literature.

We also would like to mention that the proper mathematical study of such
dynamical systems arising from realistic behavior of agents in a NEP together
with iterative schemes meant to simulate their behavior is an interesting
subject that is rarely discussed in the literature, with several nuances that
must be addressed. In this aspect, our results are definitely incomplete, as
we are mainly concerned with stationarity of a limit point of the sequence
generated, as usual in the optimization practice. That is, since we believe
that cycling or orbiting involving non-stationary points are reasonable
outcomes of the dynamics, our assumptions are presented in order to force
convergence of the algorithm, leaving the other possibilities to a later
study. With all this in mind we start stating our main convergence results.

###### Proposition 5.1

Consider the sequences $\\{x^{k}\\}$ and $\\{t^{k}\\}$ generated by Algorithm
2. Then, the following holds:

  * i)

If $\lim_{\mathbb{N}^{\prime}}d^{k}=0$ on an infinite subset
$\mathbb{N}^{\prime}\subset\mathbb{N}$, then the sequence of gradients
$\\{g_{1}^{k}\\}$ and $\\{g_{2}^{k}\\}$ converge to zero on the same subset.

  * ii)

If the whole sequence $\\{x^{k}\\}$ converges to a point $x^{*}$, then either
$\nabla_{x_{1}}f_{1}(x^{*})=0$ or $\nabla_{x_{2}}f_{2}(x^{*})=0$.

  * iii)

If $\\{x^{k}\\}$ converges to a point $x^{*}$ and if in addition the sequence
of step sizes $\\{t^{k}\\}$ is bounded away from zero, then $x^{*}$ is
stationary for the NEP (1-2), that is, $\nabla_{x_{1}}f_{1}(x^{*})=0$ and
$\nabla_{x_{2}}f_{2}(x^{*})=0$.

Proof: By (13) and (23), we have that $\|g^{k}\|=\|{\cal
H}_{t,k}d^{k}\|\leq\|{\cal H}_{t,k}\|\cdot\|d^{k}\|\leq\mu_{max}\|d^{k}\|$.
Therefore, if $\lim_{k\in\mathbb{N}^{\prime}}d^{k}=0$, then
$\lim_{k\in\mathbb{N}^{\prime}}g^{k}=0$, which proves item (i).

To prove item (iii), assume that the sequence converges. Then, we have that

$0=\lim_{k\to\infty}x^{k+1}-x^{k}=\lim_{k\to\infty}t^{k}d^{k}.$

Since $t^{k}\geq\bar{t}$ for some constant $\overline{t}>0$, we have that
$d^{k}$ converges to zero. Thus, by item (i) the whole gradient sequence
converges to zero implying that the limit point $x^{*}$ is stationary.

Finally, to prove (ii), assume that it does not hold. Thus, there exists
$\epsilon>0$ such that $\|g_{i}^{k}\|\geq\epsilon$ for $i=1,2$ and $k$ large
enough, say $k>k_{0}\in\mathbb{N}$. Additionally, since $\\{x^{k}\\}$ is
convergent, the sequence lies on a bounded set. Hence, by Proposition 4.2, we
have that $\\{t^{k}\\}$ is bounded away from zero, which by item (iii) implies
that $g^{k}\rightarrow 0$, a contradiction. $\;\;\;\Box$

Differently from the case of optimization, when the sequence $\\{x^{k}\\}$ is
not necessarily convergent, one can not expect to prove stationarity of its
limit points. However, if there is some indication of asymptotic stationarity
for one of the players, one can infer stationarity/optimality for the other
player, in some sense. In the next proposition we show that when the sequence
of iterates varies only with respect to the first player (that is,
$d_{2}^{k}=0$) and $\\{t^{k}\\}$ is bounded, then either $\\{f_{1}(x^{k})\\}$
is unbounded from below, which would somehow characterize the ill-posedness of
the problem, or all limit points of $\\{x^{k}\\}$ are stationary for the NEP
(1-2). Clearly, the role of each player may be exchanged in the statement of
this result.

###### Proposition 5.2

Consider $\\{x^{k}\\}$ and $\\{t^{k}\\}$ given by Algorithm 2 and suppose that
there exists a positive real number $\overline{t}$ such that
$t^{k}>\overline{t}$ for all $k\in\mathbb{N}$.

  * i)

If there exists $k_{0}$ such that $d_{2}^{k}=0$ for all $k\geq k_{0}$ then
either $\\{f_{1}(x^{k})\\}$ is unbounded from below, or
$\lim_{k\to+\infty}g_{1}^{k}=0$.

  * ii)

If $\sum_{k=1}^{\infty}\|d_{2}^{k}\|=\bar{d}$ for some $\bar{d}\in\mathbb{R}$,
then either $\\{f_{1}(x^{k})\\}$ is unbounded from below, or there exists an
infinite subset $\mathbb{N}^{\prime}\subset\mathbb{N}$ such that
$\lim_{k\in\mathbb{N}^{\prime}}g_{1}^{k}=0$.

As a consequence, considering that $f_{1}$ is bounded, item i) states that
every limit point of $\\{x^{k}\\}$ is stationary for the NEP (1-2) while item
ii) states that there exists a limit point that is stationary.

Proof: By conditions (14) and (15), for $k>k_{0}$ we have that

$\displaystyle f_{1}(x_{1}^{k+1},x_{2}^{k+1})$ $\displaystyle\leq
f_{1}(x_{1}^{k},x_{2}^{k+1})+t^{k}\alpha\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k+1})^{T}d_{1}^{k}$
$\displaystyle\leq
f_{1}(x_{1}^{k},x_{2}^{k+1})-t^{k}\alpha\theta\|\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k+1})\|\cdot\|d_{1}^{k}\|.$

Summing and subtracting $f(x^{k})$, using that $t^{k}\geq\bar{t}$ and the mean
value theorem, we have for some $\xi_{k}\in[0,1]$ that

$\displaystyle f_{1}(x^{k+1})-f_{1}(x^{k})$ $\displaystyle\leq
f_{1}(x_{1}^{k},x_{2}^{k+1})-f_{1}(x_{1}^{k},x_{2}^{k})-t^{k}\alpha\theta\|\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k+1})\|\cdot\|d_{1}^{k}\|$
$\displaystyle\leq
t^{k}\nabla_{x_{2}}f_{1}(x_{1}^{k},x_{2}^{k}+\xi_{k}t^{k}d_{2}^{k})^{T}d_{2}^{k}-t^{k}\alpha\theta\|\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k+1})\|\cdot\|d_{1}^{k}\|$
$\displaystyle\leq\|\nabla_{x_{2}}f_{1}(x_{1}^{k},x_{2}^{k}+\xi_{k}t^{k}d_{2}^{k})\|\|d_{2}^{k}\|-\bar{t}\alpha\theta\|\nabla_{x_{1}}f_{1}(x_{1}^{k},x_{2}^{k+1})\|\cdot\|d_{1}^{k}\|.$
(53)

Using (16) and a telescoping sum on (53), if $d_{2}^{k}=0$ for $k>k_{0}$ then

$f_{1}(x^{k+1})\leq
f_{1}(x^{k_{0}})-\bar{t}\alpha\theta\gamma\sum_{j=k_{0}}^{n}\|g_{1}^{k}\|^{2}.$

So, if $\\{f_{1}(x^{k})\\}$ is bounded, we must have that
$\lim_{k\to\infty}g_{1}^{k}=0$. Since $d_{2}^{k}=0$ implies that
$\|d_{1}^{k}\|\leq\frac{\|g_{1}^{k}\|}{\lambda_{min}}$, the result follows
from Proposition 5.1.

Let us consider now the second case. Given $\epsilon>0$, since
$d_{2}^{k}\rightarrow 0$, Assumptions B and C ensure that there exists
$k_{1}\in\mathbb{N}$ such that if $\|g_{1}^{k}\|>\epsilon$ and $k>k_{1}$ then

$\|\nabla_{x_{2}}f_{1}(x_{1}^{k},x_{2}^{k}+\xi_{k}t^{k}d_{2}^{k})\|\leq
2\|g_{1}^{k}\|,$
$\frac{1}{2}\|g_{1}^{k}\|\leq\|\nabla_{x_{2}}f_{1}(x_{1}^{k},x_{2}^{k}+t^{k}d_{2}^{k})\|\leq
2\|g_{1}^{k}\|,$
$\|d_{2}^{k}\|\leq\frac{\overline{t}\alpha\theta\gamma}{16}\|g_{1}^{k}\|.$

So, in this situation we have by (53) and (16) that

$\displaystyle f_{1}(x^{k+1})-f_{1}(x^{k})$ $\displaystyle\leq
2\|g_{1}^{k}\|\cdot\|d_{2}^{k}\|-\frac{\overline{t}\alpha\theta}{2}\|g_{1}^{k}\|\cdot\|d_{1}^{k}\|$
$\displaystyle\leq-\frac{\overline{t}\alpha\theta\gamma}{8}\|g_{1}^{k}\|^{2}\leq-\frac{\overline{t}\alpha\theta\gamma\epsilon^{2}}{8}.$
(54)

If (54) holds for all $k$ large enough we have that $\\{f_{1}(x^{k})\\}$ is
unbounded. Thus, assuming that it is bounded, there exist an infinite subset
$\mathbb{N}^{\prime}\subset\mathbb{N}$ such that
$\lim_{k\in\mathbb{N}^{\prime}}g_{1}^{k}=0$. Using the first equation of (13),
Assumptions B and D, and the fact that $d_{2}^{k}\to 0$ we conclude that
$\lim_{k\in\mathbb{N}^{\prime}}d_{1}^{k}=0$ and so the result follows from
Proposition 5.1.

$\;\;\;\Box$

The results of Proposition 5.2 are, in some sense, a generalization of what
happens when Algorithm 2 is applied to a problem where optimality for one of
the player is always satisfied, and only optimality of the other player is to
be sought. In this case, $t^{k}$ is always bounded away from zero and the
Newton direction has a norm greater than a fraction of the gradient.
Therefore, our result recovers the standard result from Newton’s method that
either the function is unbounded from below or all limit points are
stationary. Note that, by Lemma 3.3, it would be possible to include in
Algorithm 2 a test to ensure that $\gamma\|g_{i}^{k}\|\leq\|d_{i}^{k}\|$, and
therefore $\lim g_{2}^{k}=0$ would be a consequence of $\lim d_{2}^{k}=0$. We
chose not to do this explicitly since in some sense a similar requirement is
already guaranteed by (16) and (19). Finally, we note that this type of
analysis is fairly non-standard, and even a full analysis of the possible
outcomes of the NEP viewed as a dynamical system is still lacking; this
analysis would be very welcomed, as it would provide new elements for the
analysis of algorithms of the type we propose.

We conclude our convergence results by showing that when the objective
functions are strictly convex quadratic functions, the full step with $t=1$ is
always accepted by the line search and one iteration is enough to solve the
problem. This is similar to what is known for the standard Newton method for
optimization. This type of result corroborates the plausibility of the
assumption on the boundedness away from zero of the stepsize, which was very
important in the previous results.

###### Proposition 5.3

Let $A_{1}\in\mathbb{R}^{n_{1}\times n_{1}}$, $A_{2}\in\mathbb{R}^{n_{2}\times
n_{2}}$, $B_{1}\in\mathbb{R}^{n_{1}\times n_{2}}$,
$B_{2}\in\mathbb{R}^{n_{2}\times n_{1}}$, $c_{1}\in\mathbb{R}^{n_{1}}$, and
$c_{2}\in\mathbb{R}^{n_{2}}$ with $A_{1}$ and $A_{2}$ being positive definite
matrices. Consider the NEP defined by the functions
$f_{1}(x_{1},x_{2})=\frac{1}{2}(x_{1})^{T}A_{1}x_{1}+(B_{1}x_{2}-c_{1})^{T}x_{1}$
and
$f_{2}(x_{1},x_{2})=\frac{1}{2}(x_{2})^{T}A_{2}x_{2}+(B_{2}x_{1}-c_{2})^{T}x_{2}$.
If $\alpha\leq\frac{1}{2}$, $H_{1}^{0}=A_{1}$, $H_{2}^{0}=A_{2}$, and the
matrix $\left(\begin{array}[]{cc}A_{1}&B_{1}\\\ B_{2}&A_{2}\end{array}\right)$
is non-singular, then Algorithm 2 finds a solution in a single iteration
without performing backtracking; ie, the solution of system (13) with $t=1$
satisfies all inequalities (14)-(19) with constants $\gamma$ and $\theta$
satisfying (31) and (41) respectively and the new iterate $x^{1}$ is a
solution of the problem.

Proof: We prove the result for $i=1$. For $i=2$ the result is similar. Since
$f_{1}(x)$ is quadratic with respect to $x_{1}$ and linear with respect to
$x_{2}$, we have that $\nabla_{x_{1}}f_{1}(x)$ is linear. Thus, given an
arbitrary $x$ and $d$,

$\displaystyle\nabla_{x_{1}}f_{1}(x_{1}+d_{1},x_{2})=\nabla_{x_{1}}f_{1}(x)+A_{1}d_{1},\quad\text{
and }$ (55)
$\displaystyle\nabla_{x_{1}}f_{1}(x_{1},x_{2}+d_{2})=\nabla_{x_{1}}f_{1}(x)+B_{1}d_{2}.$
(56)

So, using (56) in (13) we obtain that

$A_{1}d_{1}=-\nabla_{x_{1}}f_{1}(x_{1}^{0},x_{2}^{0}+d_{2}).$ (57)

Then, by (55), $\nabla_{x_{1}}f_{1}(x^{0}+d)=0$. Since
$\nabla_{x_{1}x_{1}}^{2}f_{1}(x)=A_{1}>0$ for all $x$, we have that $x^{0}+d$
is a solution of the NEP.

Next, we show that the inequalities (14)-(19) hold, which means that the point
$x^{1}:=x^{0}+d$ is accepted by the algorithm. By (57) we have that

$\|\nabla_{x_{1}}f_{1}(x_{1}^{0},x_{2}^{0}+d_{2})\|\leq\|A_{1}\|\|d_{1}\|\leq\lambda_{max}\|d_{1}\|.$

So we obtain (16) and (19) for $\gamma\leq\frac{1}{\lambda_{max}}$. As a
consequence of this and (57)

$\displaystyle d_{1}^{T}\nabla_{x_{1}}f_{1}(x_{1}^{0},x_{2}^{0}+d_{2})$
$\displaystyle=-d_{1}^{T}A_{1}d_{1}\leq-\lambda_{min}\|d_{1}\|^{2}$
$\displaystyle\leq-\lambda_{min}\gamma\|d_{1}\|\|\nabla_{x_{1}}f_{1}(x_{1}^{0},x_{2}^{0}+d_{2})\|.$

Thus, (15) and (18) hold. Finally, to prove (14) and (17), we use that
$f_{1}(\cdot,x_{2}^{0}+d_{2})$ is quadratic and (57) to obtain

$\displaystyle f_{1}(x+d)$
$\displaystyle=f_{1}(x_{1}^{0},x_{2}^{0}+d_{2})+\nabla_{x_{1}}f_{1}(x_{1}^{0},x_{2}^{0}+d_{2})^{T}d_{1}+\frac{1}{2}d_{1}^{T}A_{1}d_{1}$
$\displaystyle=f_{1}(x_{1}^{0},x_{2}^{0}+d_{2})+\frac{1}{2}d_{1}^{T}\nabla_{x_{1}}f_{1}(x_{1}^{0},x_{2}^{0}+d_{2}).$

$\;\;\;\Box$

Proposition 5.3 guarantees convergence, in one iteration, asking only that
$A_{1}$ and $A_{2}$ are positive definite together with the non-singularity of
the matrix $M:=\left(\begin{array}[]{cc}A_{1}&B_{1}\\\
B_{2}&A_{2}\end{array}\right)$, without requiring its positive definiteness.
In particular, without any special requirement on the matrices $B_{1}$ and
$B_{2}$. In contrast, a standard Jacobi method for the system of equations (3)
would require the norm of the iteration matrix to be less than $1$ (which can
be assured, for instance, when $M$ is strictly diagonally dominant), while the
Jacobi trust-region method proposed by [28] would require for convergence,
among other things, the positive definiteness of $M$.

Additionally, the fact that Proposition 5.3 ensures convergence for quadratic
problems in one iteration, hints at the behavior of the method close to the
solution when applied to non-quadratic functions; namely, for sufficiently
well-behaved functions, it is expected that the behavior of the algorithm
would be similar to its behavior on quadratic functions. Thus, we would expect
that close to the solution the method would not need to decrease the stepsize,
and local convergence of the pure Newtonian strategy would be recovered,
however the pursue of a result of this type is out of the scope of this paper.

##  6 Numerical Experiments

In this section we present some numerical experiments depicting the behavior
of Algorithm 2. First we test it on five selected illustrative example
problems and then we test the Algorithm on an application on facility location
problems.

###  6.1 Illustrative Examples

For the first experiments, we chose the initial point as $x^{0}:=(-5,1)^{T}$,
an error tolerance of $10^{-4}$ for the norm of the gradients $\|g^{k}\|$ as
stopping criterion, together with a maximum number of iterations of $k=1000$.
In Table 1 we show the convergence point, the error, and the number of
iterations for the five selected equilibrium problems to be defined below. We
compared Algorithm 2 (abbreviated by Alg1 in the tables) with the Jacobi-type
trust-region method of [28] (abbreviated Yuan), an exact Jacobi-type method
that deals with system (3) directly at each step (abbreviated Jacobi), and
Newton’s method for system (3) with unitary step, abbreviated Newton. We also
display each problem’s actual solution for comparison.

Jacobi’s method for system (3) is done as follows: given
$(x_{1}^{k},x_{2}^{k})$ we algebrically compute the next iterate
$(x_{1}^{k+1},x_{2}^{k+1})$ by solving the equation
$\nabla_{x_{1}}f_{1}(x_{1},x_{2}^{k})=0$ for $x_{1}$ and
$\nabla_{x_{2}}f_{2}(x_{1}^{k},x_{2})=0$ for $x_{2}$, and then update
$(x_{1}^{k+1},x_{2}^{k+1}):=(x_{1},x_{2})$. As for Newton’s method, given an
iterate $(x_{1}^{k},x_{2}^{k})$, we compute $(d_{1}^{k},d_{2}^{k})$ by solving
system (4) and then compute $(x_{1}^{k+1},x_{2}^{k+1})$ as
$(x_{1}^{k},x_{2}^{k})+(d_{1}^{k},d_{2}^{k})$.

Yuan’s algorithm, a Jacobi-type method for NEPs that also uses Newtonian
ideas, depends on several parameters, for details see [28]. In all our
implementations we chose the parameters
$\delta_{v}:=0.01,\tau_{v}:=1,t_{v,1}:=1$, $G_{i}=1$, and
$\beta_{i}:=\frac{1}{2},i=1,2$.

For our method, we chose $\alpha:=10^{-6},\theta=0.01,\gamma=10^{-6},$ and
$\tau=0.99$. Recall that according to Lemmas 4.1 and 4.2, $\gamma$ is
proportional to $\frac{1}{\lambda_{max}}$ and $\theta$ to the ratio
$\frac{\lambda_{min}}{\lambda_{max}}$. Selecting parameters $\gamma$ and
$\theta$ in this way, we enforce a large bounding interval for
$\|H_{i}^{k}\|$, which guarantees the validity of Assumption D. In order to
choose the matrices $H_{i}^{k}$, we computed the modified Cholesky
decomposition of $\nabla^{2}_{x_{i}x_{i}}f_{i}(x^{k})$ as described in [26],
which resumes to this matrix whenever it is positive definite, and, when it is
not, the procedure returns a positive definite matrix by adding a small
diagonal perturbation to the matrix.

All tests were run in Matlab R2017a in an AMD Ryzen 5 2400G with 3.6Ghz
graphics and 8Gb RAM processor. Our illustrative problem set is given by the
following simple one-dimensional NEPs:

###### Example 6.1

$f_{1}(x):=x_{1}^{2}+x_{1}x_{2}-5x_{1},\quad
f_{2}(x):=\frac{3x_{2}^{2}}{2}-x_{1}x_{2}-x_{2}.$

This example defines positive definite quadratic functions such that solving
its first-order necessary conditions is sufficient for solving the problem, in
addition, the underlying linear system is strictly diagonally dominant; thus,
we expect the problem to be solved properly even by the exact Jacobi-type
method.

###### Example 6.2

$f_{1}(x):=\frac{x_{1}^{2}}{4}+x_{1}x_{2}-5x_{1},\quad
f_{2}(x):=\frac{x_{2}^{2}}{6}-x_{1}x_{2}-x_{2}.$

This example is also a positive definite quadratic, however the spectral
radius of the Jacobi iteration matrix is greater than 1, which may hinder the
convergence of the exact Jacobi-type method.

###### Example 6.3

$f_{1}(x):=x_{1}^{2}+x_{1}x_{2}-5x_{1},\quad
f_{2}(x):=-\frac{3x_{2}^{2}}{2}-x_{1}x_{2}-x_{2}.$

In this example we considered a game without an equilibrium solution, given
that $f_{2}$ has no local minimizer for all $x_{1}$ fixed. However, this
problem has a non-equilibrium stationary point at $(3.2,-1.4)$.

###### Example 6.4

$f_{1}:=-x_{1}(0.6-x_{2}),\quad f_{2}:=x_{2}(0.7-x_{1}).$

This example comes from a mixed strategies formulation on a discrete
equilibrium game, where the government needs to decide between the
distribution of two different types of vaccines for a viral disease, for more
details see [4]. This kind of problem has the special property that the
Hessians are null, meaning that Jacobi’s method is undefined.

###### Example 6.5

$f_{1}(x):=\frac{x_{1}^{3}x_{2}^{2}}{3}+\frac{x_{1}^{2}}{2},\quad
f_{2}(x):=\frac{x_{1}^{2}x_{2}^{3}}{3}+\frac{x_{2}^{2}}{2}.$

In the fifth example, we chose a cubic problem with an equilibrium point in
$(0,0)$ and a stationary point in $(-1,-1)$, which is not an equilibrium.

The results are given in Table 1. For each method we report the last point
obtained, the norm of $g^{k}$, and the number of iterations used by each
method.

Table 1: Test results

Example 6.1: Quadratic with matrix with spectral radius $<1$ \- Solution:
$(2,1)$
---
Method | Alg1 | Yuan | Jacobi | Newton
Point | $(2,1)$ | $(1.99997,0.99999)$ | $(2,1)$ | $(2,1)$
$\|g^{k}\|$ | 0 | $9.85431\cdot 10^{-5}$ | $0$ | 0
# iter | 1 | 1 | 2 | 1
Example 6.2: Quadratic with matrix with spectral radius $>1$ \- Solution:
$(0.57142,4.71428)$
Method | Alg1 | Yuan | Jacobi | Newton
Point | $(0.57142,4.71428)$ | $(0.57145,4.71419)$ | $(\infty,\infty)$ | $(0.57142,4.71428)$
$\|g^{k}\|$ | $2.22009\cdot 10^{-16}$ | $9.71940\cdot 10^{-5}$ | $\infty$ | $2.22009\cdot 10^{-16}$
# iter | 1 | 158 | 398 | 1
Example 6.3: Quadratic problem without a solution.
Method | Alg1 | Yuan | Jacobi | Newton
Point | $(-\infty,\infty)$ | $10^{3}\cdot(-0.49401,0.99900)$ | $(3.19999,-1.39999)$ | $(3.2,-1.4)$
$\|g^{k}\|$ | $\infty$ | $2.50109\cdot 10^{3}$ | $4.28624\cdot 10^{-5}$ | $1.77615\cdot 10^{-16}$
# iter | 372 | 1000 | 8 | 1
Example 6.4: Problem with null Hessians - Solution: $(0.7,0.6)$
Method | Alg1 | Yuan | Jacobi | Newton
Point | $(0.7000,0.60001)$ | $10^{9}\cdot(1.86823,-0.62812)$ | - | $(0.7,0.6)$
$\|g^{k}\|$ | $7.00011\cdot 10^{-5}$ | $2.59973\cdot 10^{9}$ | - | $2.22009\cdot 10^{-16}$
# iter | 1 | 1000 | - | 1
Example 6.5: Cubic problem - Solution: $(0,0)$
Method | Alg1 | Yuan | Jacobi | Newton
Point | $10^{-8}\cdot(0.27025,0.27025)$ | $10^{-4}\cdot(0.64551,0.31008)$ | $(-10^{-154}\cdot 0.74020,-\infty)$ | $(-1,-1)$
$\|g^{k}\|$ | $5.40980\cdot 10^{-9}$ | $8.50096\cdot 10^{-5}$ | $\infty$ | $5.04872\cdot 10^{-16}$
# iter | 9 | 25 | 6 | 7

In Examples 6.1 and 6.2 we obtained the expected results: Jacobi managed to
converge for the first one but not for the second one, while all the other
three methods found the solution. Algorithm 2 and Newton’s method converged in
one iteration since the problems are quadratic with positive definite
Hessians.

In Example 6.3, which does not have a solution, the exact Jacobi and Newton’s
methods converged to the stationary point $(3.2,-1.4)$, meaning that they did
not aim at solving the actual problem. Our algorithm and Yuan’s algorithm, on
the other hand, did not converge, as should be expected. Yuan’s method
performed the maximum number of iterations, while Algorithm 2 found
$(-\infty,+\infty)$ in $372$ iterations.

In Example 6.4 the exact Jacobi method is undefined, since the objective
functions’ Hessians $\nabla_{x_{i}x_{i}}^{2}f_{i}(x^{k})$ are null, meaning
that the component-wise systems in (3) are undefined. Since the problem is
quadratic, Newton’s method managed to solve it in one step, as expected. Our
method does not use the Newton step here, as the Hessians are not positive
definite, thus implying it had to compute a positive definite approximation.
Even so, the method still managed to find a solution with the required
precision in just one step. We highlight that Yuan’s method could not find the
solution in this example, due to the fact that Yuan’s method rely on strong
hypotheses on the system matrix that demand the positive definiteness of the
jacobian of the function that defines system (3), which does not occur here.
On the other hand, our method managed to reach the solution, since its theory
is based on positive definite approximated Hessians instead.

As for Example 6.5, since the gradients are given by $x_{i}(x_{i}x_{\neg
i}^{2}+1)$ for each component $i$, the stationary points can be achieved
either if $x_{i}=0$ or if $x_{i}x_{\neg i}^{2}+1=0$. In order to implement the
exact Jacobi method for system (3), we need to choose how we are going to
satisfy the system directly, either by selecting $x_{i}=0$ (the equilibrium
point) or selecting $x_{i}$ satisfying $x_{i}x_{\neg i}^{2}+1=0$. If we choose
the former, Jacobi’s method unsurprisingly takes one step. If instead we
choose the latter, Jacobi’s method behaves as follows: one of the components
diverges decreasing the function value indefinitely and the other component
manages to find a point that zeroes the respective gradient. On the other
hand, Algorithm 2, Yuan’s and Newton’s method successfully managed to find a
stationary point in this example. However, Newton’s method converged to
$(-1,-1)$, the closest stationary point, instead of to the actual solution.
Algorithm 2 and Yuan’s method, on the other hand, converged to the actual
solution $(0,0)$.

We highlight that Algorithm 2 never reduced the stepsize in order to achieve
the conditions (14)-(19) in the quadratic Examples 6.1-6.4. This is expected
and consistent with Proposition 5.3 for the strictly convex cases. In Example
6.5 the stepsize was reduced at some iterations, however, at most one
reduction of the stepsize was needed for the direction computed to be accepted
by the descent criteria (14)-(19). Even when the stepsize was reduced, this
occurred only far from the solution, which is also consistent with a Newtonian
approach and indicates that the hypothesis of the sequence $\\{t^{k}\\}$ being
bounded away from zero is plausible.

###  6.2 Facility location problem

In the facility location problem, first studied in [15], one may choose the
best place to open facilities. In general, servers seek to be close to their
customers due to the various commercial and logistical advantages that this
can entail. Here we consider the problem of two players who each want to open
a facility in order to obtain the highest expected return on trading with $N$
customers. A player’s chances of serving a certain customer depends on the
proximity of their facility with the position of the customer in relation to
the competitor’s proximity. Denoting by $x_{i}$ the position of the facility
$i$, $i\in\\{1,2\\}$, and by $z_{j}$ the position of the client $j$,
$j\in\\{1,\cdots,N\\}$, let us consider here that the probability of player
$i$ serving client $j$ is
$1-\frac{\|x_{i}-z_{j}\|^{2}}{\|x_{i}-z_{j}\|^{2}+\|x_{\neg i}-z_{j}\|^{2}}$.
Thus, if $b_{j}^{i}$ is the profit of player $i$ in serving customer $j$, the
function to be minimized by each player to obtain the best expected revenue is

$f_{i}(x_{i},x_{\neg
i})=\displaystyle\sum_{j=1}^{N}\frac{b^{i}_{j}\|x_{i}-z_{j}\|^{2}}{\|x_{i}-z_{j}\|^{2}+\|x_{\neg
i}-z_{j}\|^{2}},\quad i=1,2,$ (58)

resulting in a NEP.

Notice that the functions are not defined when both facilities are located
exactly at the position of one of the customers, a very pathological situation
which was never encountered by any of the algorithms we run. Anyway, these
problems are complex, since the functions are non-convex. In Figure 1 we
illustrate the objective function of the first player in a one-dimensional
problem, with the decision of the second player fixed at $0.915$ and three
clients located at $z_{1}=1,z_{2}=-1$ and $z_{3}=3$, with a uniform profit of
$1$. For this problem an equilibrium solution is $(1.901,0.915)$, which is
found by Algorithm 2, with the same parameters used in Subsection 6.1, from
the starting point $(2,1)$.

Figure 1: Unidimensional objective function for the facility location problem.

For a more extensive analysis, let us now consider a two-dimensional problem.
For our test we selected a problem in the form (58) with the clients’
coordinates $z_{1}=(1,0),z_{2}=(0,1),z_{3}=(-1,0)$, and $z_{4}=(0,-1)$. We
chose the profit vector for each player as $b^{1}=(1,2,1,1)$ and
$b^{2}=(1,2,2,3)$. We run Algorithm 2 for $100$ random initial points with
coordinates generated in the interval $[-2,2]$, comparing with Yuan’s, Jacobi
and Newton’s method. Here, for the Jacobi method we solved problems (5) and
(6) with $\phi_{i}=f_{i}$ using Matlab’s fminunc solver instead of dealing
with system (3). In Figure 2 we report the final outcome of each algorithm
according to the stopping criterion $\|g^{k}\|\leq 10^{-6}$. We report simply
whether the method converged to an equilibrium point, a non-equilibrium
stationary point or if it diverged. We see that Jacobi and Yuan’s method
perform poorly due to the properties of the Jacobian of the correspondent
first-order conditions. Jacobi diverged in $88$ out of $100$ runs while Yuan’s
method diverged in $97$ runs, reaching the maximum number of iterations.
Algorithm 2 converged to a true equilibrium point in all cases, taking on
average $36$ iterations and $1.08$ seconds of CPU time, while Newton’s method
mostly converged to a non-equilibrium stationary point, having attained the
true equilibrium point in only $13$ of the $100$ runs. Newton’s method, on
average, took $19$ iterations and $0.76$ seconds to converge, however, for the
$13$ initial points where both algorithms converged to the equilibrium, they
were able to perform similarly as Algorithm 2 took $0.82$ seconds to perform
$19$ iterations while Newton’s method took $0.75$ seconds to perform $16$
iterations, on average.

In Figure 3 we see the behavior of Algorithm 2 for this problem for the
initial point $x^{0}_{1}=(2,3)$ and $x^{0}_{2}=(-3,2)$. The filled circles in
the plane represent the optimal facility locations found for the agents, while
an empty circle represents the iterates. We see that Algorithm 2 forces the
trajectory to the optimal placement. Running Newton’s method with this same
data yields a sequence that diverges, forcing the second facility to
$(+\infty,-\infty)$, yet still stopping declaring success since the gradients
tend to zero.

Figure 2: Facility allocation problem: quality of convergence for different
methods. Figure 3: Bidimensional behavior of convergence of Algorithm 2.

In conclusion, our numerical experiments attest that our method may avoid
converging to a non-equilibrium stationary point, differently from methods
that deal only with the nonlinear system of equations obtained by the first
order necessary optimality conditions as the standard Newton method. Also,
differently from Jacobi and Yuan’s method, our method is based on a predicted
behavior of the other player, which seems to considerably speed up
convergence. Moreover, Jacobi and Yuan’s method seem to struggle when the
Jacobian of the system of equations is not well-behaved, which does not occur
with our method.

##  7 Final remarks

The formulation of Nash Equilibrium Problems (NEPs) is one of the most
important concepts in modeling the behavior of social and economic agents in
practice. In particular, in this paper, we considered the unconstrained NEP
with two players and continuous variables, where we focus on the possibility
of non-convexities in the objective functions of each player. Usually, for
this kind of problem in the optimization literature, one is interested only in
computing a Nash equilibrium, or, more accurately, a point satisfying the
system of nonlinear equations that are satisfied by all equilibrium points.
This approach may be misleading as this system of equations may include also
other undesired solutions.

The contribution of this paper is to propose an algorithm for this problem
with two main characteristics: Firstly, our algorithm addresses the NEP in its
entirety, considering its minimization structure and avoiding convergence to
an undesired stationary point. Secondly, the iterates of our algorithm are
computed in such a way that they can be interpreted as mimicking the behavior
of the agents in practice, namely, the iterates computed by the method
approximate in some sense the true decision made by the players in different
time periods. That is, a new iterate is computed by player one by minimizing a
convex quadratic model of their objective function parameterized by a
prediction of the action that will be taken by the second player. The
prediction is computed simultaneously, assuming that the second player is
adopting the same strategy. This gives rise to a Jacobi-type Newton method,
globalized by an Armijo-type linesearch, which guarantees that at the new
iterate both players will decrease their predicted objectives; this avoids, as
much as possible, convergence to undesired stationary points. On the other
hand, the method can also be seen as a new interpretation of Newton’s method,
which is known to be intrinsically linked to fast local convergence.

Both characteristics of our algorithm are new and relevant for the literature
of NEPs, especially when non-convexities are considered, favoring convergence
to a true equilibrium rather than simply a stationary point. Moreover, our
approach opens the path to new and interesting investigations of the behavior
of Newtonian iterates built with the goal of mimicking the behavior of true
economic agents. Clearly, this is a starting point for this type of
investigation and many aspects still need to be elucidated. For example, it
would be interesting to interpret the many parameters of the algorithm with
the true agents’ behavior. We believe that the synergy between the
interpretation of the game’s dynamics and the sequence generated by an
algorithm can inspire both better techniques for solving NEPs and a better
understanding of practical situations associated with models of this type.

In this paper, we focused on the analysis considering convergence of the
iterative sequence when the step size sequence $\\{t^{k}\\}$ is bounded away
from zero. An important open problem for a better understanding of the
algorithm would be to prove whether or not it is possible to guarantee the
boundedness away from zero of $\\{t^{k}\\}$. Additionally, it would be
interesting to investigate whether or not the step $t=1$ is always accepted
when close enough to solutions that satisfy the second-order sufficient
optimality condition. This would be associated with both the boundedness away
from zero of $\\{t^{k}\\}$ and the quadratic convergence of the algorithm.
However, more refined analysis should be carried out to also capture other
types of phenomena that may occur to the sequence besides convergence. For
instance, in a general dynamical system, the interest goes far beyond mere
convergence to an equilibrium, since several other interesting phenomena may
occur. We expect that subsequent research on this topic should focus on
understanding and classifying such phenomena, with the analysis carried out to
the algorithms built for mimicking the dynamics.

Taking into account that our algorithm converges more often to a true
equilibrium, we envision several possible extensions of this work. One may
consider $N>2$ players, while also considering that each player’s decision
variable $x_{i}$, $i=1,2,\dots,N$, is constrained to a feasible set
$X_{i}\subset\mathbb{R}^{n_{i}}$, or, more generally, the feasible set may
also be parameterized by the other players’ decision $x_{\neg i}$, namely,
$X_{i}:=X_{i}(x_{\neg i})$, as in what is known as a generalized NEP. Our
intention is to extend the ideas of this paper to consider such parameterized
constrained problems, borrowing ideas from Sequential Quadratic Programming,
where one considers a linearization of the constraints in a Newtonian
framework. In addition, our unconstrained algorithm may be employed in an
Augmented Lagrangian framework for the general constrained problem (see [16]
and [5]). In this framework, the constraints are penalized and at each
iteration of the algorithm an unconstrained NEP must be solved, and finding
true solutions of the subproblems rather than mere stationary points should
favor the algorithm in finding true solutions of the original problem.

Another interesting point to be investigated in the future is the connection
with dynamic learning, reinforcing decisions that have better objective
function values. Since the philosophy of using best response functions
together with predicted decisions for the other player is similar to our
approach, this may generate relevant ramifications when substituting the
standard two step process of independent decision making and probability
distribution updated considered in [22, 17] for ours. Furthermore, in
stochastic versions of the algorithm, the notions of cycling and orbiting may
be more naturally understood as generalized convergence, simplifying the
analysis.

## Acknowledgments

This work was supported by FAPESP (São Paulo Research Foundation) grants
2013/07375-0, 2015/02528-8, 2018/24293-0, and 2021/05007-0, and CNPq.

## References

  * [1] Aubin, J. P. (1993). Optima and equilibria. Berlin, Heidelberg, New York: Springer.
  * [2] Arrow, K. J., Debreu, G. (1954). Existence of an equilibrium for a competitive economy. Econometrica, 22 (pp. 265-290).
  * [3] Bertsekas, D. P. (2016). Nonlinear Programming. Athena Scientific.
  * [4] Bortolossi, H., Garbagio, G., Sartini, B. (2007). Uma Introdução à Teoria Econômica dos Jogos. 26th Colóquio Brasileiro de Matemática.
  * [5] Bueno, L. F., Haeser, G., Rojas, F. N. (2019). Optimality Conditions and Constraint Qualifications for Generalized Nash Equilibrium Problems and their Practical Implications. SIAM Journal on Optimization, 29 (1) (pp. 31-54).
  * [6] Caruso, F., Ceparano, M., Morgan, J. (2018). Uniqueness of Nash equilibrium in continuous two-player weighted potential games. J. Math. Anal. Appl., 459 (pp. 1208-1221). for Nash equilibria. SIAM Journal on Optimization, 30 (01) (pp. 1638-1663).
  * [7] Caruso, F., Ceparano, M., Morgan, J. (2020). An inverse-adjusted best response algorithm for Nash equilibria. SIAM Journal on Optimization, 30 (01) (pp. 1638-1663).
  * [8] Chui, Y., Pang, J. S. (2021). Modern Nonconvex Nondifferentiable Optimization. MOS-SIAM Series on Optimization.
  * [9] Facchinei, F., Pang, J. S. (2003). Finite-Dimensional Variational Inequality and Complementarity Problems. Springer, vols. 1 and 2.
  * [10] Facchinei, F., Kanzow, C. (2007). Generalized Nash equilibrium problems. 4OR 5 (pp. 173-210).
  * [11] Facchinei, F., Kanzow, C. (2009). Penalty methods for the solution of generalized Nash equilibrium problems. SIAM Journal on Optimization 20 (5) (pp. 2228-2253). Recent Advances and Challenges. Pesquisa Operacional, 34 (3) (pp. 521-558).
  * [12] Fischer, A., Herrich, H., Schonefeld, K. (2014). Generalized Nash Equilibrium Problems - Recent Advances and Challenges. Pesquisa Operacional, 34 (3) (pp. 521-558).
  * [13] Fukushima, M., Pang, J. S. (1998). Some feasibility issues in mathematical programs with equilibrium constraints. SIAM Journal on Optimization 8 (pp. 673-681).
  * [14] Izmailov, F. A., Solodov, M. V. (2014). Newton-type Methods for Optimization and Variational Problems. Springer.
  * [15] Hotelling, H. (1929). Stability in competition. The Economic Journal, 39 (pp. 41-57). of generalized Nash equilibrium problems. SIAM Journal on Optimization, 26 (4) (pp. 2034-2058).
  * [16] Kanzow, C., Steck, D. (2016). Augmented lagrangian methods for the solution of generalized Nash equilibrium problems. SIAM Journal on Optimization, 26 (4) (pp. 2034-2058).
  * [17] Leslie, D. D., Perkins, S., Xuc, Z. (2020). Best-response dynamics in zero-sum stochastic games. Journal of Economic Theory, 189.
  * [18] Majig, M-A., Fukushima, M. (2008). Restricted-step Josephy-Newton method for general variational inequalities with polyhedral constraints. Technical report, Kyoto University.
  * [19] Myerson, R. B. (1991). Game theory. Analysis of conflict. Harvard University Press.
  * [20] Nash, J. (1951). Non-cooperative games. Annals of mathematics, 54 (2) (pp. 286-295).
  * [21] Nikaido, H., Isoda, K. (1955). Note on noncooperative convex games. Pacific Journal of Mathematics, 5 (pp. 807-815).
  * [22] Ozdaglar, A., Sayin, M. O., Zhang, K. (2021). Independent Learning in Stochastic Games. Preprint available at https://doi.org/10.48550/arXiv.2111.11743. Palis, J., Melo, W. (1982). Geometric Theory of Dynamical Systems. Springer.
  * [23] Palis, J., Melo, W. (1982). Geometric Theory of Dynamical Systems. Springer. Rapoport, A., Chammah, A. (1965). Prisoner’s Dilemma. University of Michigan Press.
  * [24] Rapoport, A., Chammah, A. (1965). Prisoner’s Dilemma. University of Michigan Press.
  * [25] Roughgarden, T., Tardos, E. Bounding the inefficiency of equilibria in nonatomic congestion games. Games and Economic Behavior 47 (2) (pp. 389 - 403).
  * [26] Schnabel, R. B., Eskow, E. (1990). A New Modified Cholesky Factorization. SIAM J. Sci. Stat. Comput, 11 (6) (pp. 1136-1158).
  * [27] Weibull, J. W. Evolutionary Game Theory (1997). MIT Press. of Optimization, 7 (1) (pp. 125-138).
  * [28] Yuan, Y-X. (2011). A trust region algorithm for Nash equilibrium problems. Pacific Journal of Optimization, 7 (1) (pp. 125-138).