\section{Introduction}

Some classes of materials exhibit behavior that is intrinsically random. A compelling example is set forth by brittle materials with random tensile strength \cite{Bower:2010}. Structural members made of such materials may fail---and lose their load-bearing capacity altogether---at tensile stresses that are sensitive to critical flaws in the material. Since such flaws have random geometry, the tensile strength of the material is likewise random (cf., e.~g., refs.~\cite{Reynolds:1974, Bennett:1983} concerning the example of carbon fibers). In addition, empirical material data may be noisy due to experimental measurement error and scatter. This observational noise may be reduced by improving the precision of the experimental measurements. By contrast, the intrinsic randomness of the material is not eliminated by removing observational noise.

\begin{figure}[H]
	\begin{subfigure}{0.45\textwidth}\caption{} \includegraphics[width=0.99\linewidth]{Figures/Data4.pdf}
	\end{subfigure}
    $\quad$
	\begin{subfigure}{0.45\textwidth}\caption{} \includegraphics[width=0.99\linewidth]{Figures/Data5.pdf}
	\end{subfigure}
    \caption{Brittle material with random tensile strength, member $e$ of the structure.
    a) Local material likelihood function $L_{D,e}$ providing the likelihood that the local material state $(\epsilon_e,\sigma_e)$ be on the unfailed and failed branches of the local material set, respectively. b) Empirical local likelihood measure $\mu_{D,e}$ sampled from $L_{D,e}$.} \label{uCXq6P}
\end{figure}

The behavior of random materials is characterized by a likelihood measure instead of a constitutive relation. In principle, any local material state $(\epsilon_e,\sigma_e)$ is possible, albeit with varying likelihood. For instance, for a brittle material with random tensile strength undergoing monotonically increasing strains, the stress may take values in anyone of two branches, corresponding to a failed or un-failed material, with respective likelihoods determined by the distribution of tensile strength, Fig.~\ref{uCXq6P}a.

In addition, local stresses and strains in solids and structures are subject to compatibility and equilibrium constraints. These field constraints are known exactly, though the loading may also be random in some applications. The inference problem is then to determine the likelihood of structural outcomes of interest when the states of the structure are required to be admissible, in the sense of satisfying compatibility and equilibrium, and the material is characterized by a given material likelihood measure.

However, the requisite material likelihood measures are often unknown exactly and characterized only partially through empirical material data, Fig.~\ref{uCXq6P}b. A common response to that challenge is to model the empirical material data, e.~g., by assuming the existence of an underlying reduced manifold in phase space masked by noise obeying a prespecified prior distribution. However, in general no such reduced manifold may exist and the form of the prior may be unknown and subject to surmise, which inevitably introduces modeling bias and uncertainty and may prevent convergence with respect to the empirical data altogether \cite{Owhadi:2015}.

For instance, for brittle materials with random tensile strength the material states may occupy one of two branches, corresponding to a failed or un-failed material, and no single-valued stress-strain relation exists, Fig.~\ref{uCXq6P}a. In addition, in approaches based on Bayesian statistics (cf., e.~g., \cite{Dashti:2017} for a review) it is common to presume a certain class of priors, most commonly Gaussian \cite{Knapik:2011, Knapik:2016}, which may not contain the actual likelihood function of the material. For instance, if the tensile strength of a brittle material obeys a Weibull distribution (cf., e.~g., \cite{Bader:1993, Naito:2008} for the example of carbon fibers), methods assuming a Gaussian prior inevitably converge to the wrong answer with increasing sample size.

In this work, we present an alternative method of inference that determines likelihoods of outcomes directly from the empirical material data and requires no material or prior-distribution {\sl ans\"atze}. We begin by 'thermalizing' the data, which effectively replaces the empirical material likelihood by a sum of Gaussians of a certain width \cite{Conti:2021}. The thermalized likelihoods also follow variationally by minimizing a regularized Kullback-Leibler divergence \cite{Kullback:1951, Kullback:1959, Pinski:2015} from the empirical material likelihood. Likewise, an optimal annealing schedule follows variationally by minimizing the regularized Kullback-Leibler divergence with respect to temperature. The sequence of thermalized likelihoods thus obtained is then used to approximate expectations of outcomes. In this manner, the computation of expectations is reduced to explicit sums over local material data sets and quadratures over admissible states, i.~e., states satisfying compatibility and equilibrium. When restricted to the computation of maximum-likelihood outcomes, the present approach reduces to the minimum-distance (min-dist) Data-Driven method of \cite{kirchdoerfer2016data} in the deterministic case with uniform data convergence, and to the maximum-entropy (max-ent) Data-Driven method of \cite{kirchdoerfer2017data} in the case of noisy data with outliers.

It bears emphasis that, in the present Data-Driven approach, expectations are computed directly from the material data and that no material model or prior distribution need to be contrived at any point of the calculations. The complexity of the Data-Driven inference calculations is linear in the number of material data points and in the number of members in the structure or Gauss points in the solid. We formulate efficient population annealing procedures and fast-search algorithms for accelerating calculations further. Specifically, the population annealing algorithm reduces the complexity of the calculations and it renders it linear in the size of the material data-point population. The calculations are further accelerated by recourse to a hierarchical $k$-means algorithm for the computation of the annealing energy with limits on the number backtracking operations.

The convergence properties of the method with respect to the empirical data are assessed with the aid of selected applications and benchmark tests, including a standard verification test based on Gaussian material data, brittle materials exhibiting random tensile strength, and a simple lightweight space structure. Robust convergence of posterior distributions is obtained in these tests, consistent with quantitative error estimates derived from analysis \cite{Conti:2021, Conti:2022}. The ability of the proposed Data-Driven method of inference to deal effectively with general material data sets and complex behavior without need for models, hypotheses or assumptions is quite remarkable. The numerical tests additionally attest to feasibility of the calculations and the effectiveness of acceleration methods such as population annealing and hierarchical $k$-means searches.

\section{The inference problem}
\label{M2Xt7P}

We consider finite-dimensional systems comprising $m$ components whose state is characterized by two work-conjugate fields $\epsilon \equiv \{\epsilon_e \in \mathbb{R}^d,\ e=1,\dots,m\}$ and $\sigma \equiv \{\sigma_e \in \mathbb{R}^d,\ e=1,\dots,m\}$. We refer to the space of pairs $Z_e = \{z_e \equiv (\epsilon_e, \sigma_e) \in \mathbb{R}^d \times \mathbb{R}^d\}$ as the {\sl local phase space} of component $e$, and $Z = Z_1 \times \cdots \times Z_m = \mathbb{R}^{N\times N}$, $N = m d$, as the {\sl global phase space} of the system. We suppose that a suitable norm is defined in $Z$, e.~g.,
\begin{equation}\label{Poyet2}
    \| z \|
    =
    \| (\epsilon,\sigma) \|
    =
    \left(
    \sum_{e=1}^m
        w_e
        \Big(
            {\mathbb{C}}_e \epsilon_e \cdot \epsilon_e
            +
            {\mathbb{C}}_e^{-1} \sigma_e \cdot \sigma_e
        \Big)
    \right)^{1/2} ,
\end{equation}
where $w_e > 0$ are weights and $\mathbb C_e\in \mathbb R^{d\times d}_{{\rm sym},+}$ are positive definite symmetric matrices, $e=1,\dots,m$.

\begin{example}[Trusses]\label{nBpq7d}
We illustrate the essential structure of discrete field theories by means of the simple example of truss structures. Trusses are assemblies of bars that deform in uniaxial tension or compression. The bars are articulated at common joints, or nodes, that act as hinges, i.~e., cannot transmit moments. Trusses are examples of connected networks that obey conservation laws. Other examples in the same class include electrical circuits, pipeline networks, traffic networks, and others.

The material behavior of a bar $e$ is characterized by a particularly simple relation between uniaxial strain $\epsilon_e$ and uniaxial stress $\sigma_e$. Thus, in this case the local phase spaces are $Z_e = \mathbb{R} \times \mathbb{R}$. These local states are subject to the following laws:

i) Compatibility: Suppose that bar $e$ is connected to nodes $a$ and $b$. Then, the strain in the bar is
\begin{equation}
    \epsilon_e = \frac{u_b - u_a}{L_e} \cdot d_e ,
\end{equation}
where $L_e$ is the length of the bar and $d_e$ is the unit vector pointing from $a$ to $b$.

ii) Equilibrium: Let $S_a$ be the star of an unconstrained node $a$, i.~e., the collection of members connected to $a$. Then, we must have
\begin{equation}
    \sum_{e\in S_a} \sigma_e d_e A_e + f_a = 0,
\end{equation}
where $A_e$ is the cross-sectional area of bar $e$ and $f_a$ is the force applied to node $a$. \hfill$\square$
\end{example}

\subsection{The constraint set of admissible states}\label{sec:Cset}

The above example shows that the state of the system is subject to linear constraints of the general form
\begin{subequations}\label{9qHWzU}
\begin{align}
    &
    \sum_{e=1}^m w_e B_e^T \sigma_e = f ,
    \\ &
    \epsilon_e = B_e u + g_e , \quad e = 1,\dots m ,
\end{align}
\end{subequations}
where $u \in \mathbb{R}^{n}$ is the array of degrees of freedom of the system, $w_e$ are positive weights, $B_e \in \mathbb{R}^{d \times n}$ is a discrete gradient operator, $B_e^T$ is a discrete divergence operator, $f \in \mathbb{R}^n$ is a force array resulting from distributed sources and Neumann boundary conditions and the arrays $g_e \in \mathbb{R}^{d}$ follow from Dirichlet boundary conditions. The constraints (\ref{9qHWzU}) are material independent and define an affine subspace $E$ of $Z$, the {\sl constraint set}. The constraint set $E$ encodes all the data of the problem, including geometry, loading and boundary conditions. The constraints (\ref{9qHWzU}) can also be expressed in matrix form as
\begin{subequations}\label{vZftZT}
\begin{align}
    &
    B^T \tau = f
    \\ &
    \epsilon = B u + g ,
\end{align}
\end{subequations}
with $B = (B_1, \dots, B_m) \in \mathbb{R}^{N \times n}$, $\epsilon = (\epsilon_1, \dots, \epsilon_m)  \in \mathbb{R}^N$, $\tau = (w_1 \sigma_1,\dots,$ $w_m \sigma_m)$  $\in$ $\mathbb{R}^N$ and $g = (g_1, \dots, g_m)  \in \mathbb{R}^N$.

We note that the affine space $E$ defined by the constraints (\ref{vZftZT}) is a translate of the linear space $E_0$ defined by the homogeneous constraints
\begin{subequations}
\begin{align}
    &   \label{qmS3Q5}
    B^T \tau = 0
    \\ &  \label{jJ9B6T}
    \epsilon = B u .
\end{align}
\end{subequations}
Evidently, $E_0 = E_\epsilon \times E_\sigma$, where $E_\epsilon$ is the linear space defined by (\ref{jJ9B6T}) and $E_\sigma$ is the linear space defined by (\ref{qmS3Q5}). Therefore, we have
\begin{equation}
    {\rm dim}(E_0)
    =
    {\rm dim}(E_\epsilon)
    +
    {\rm dim}(E_\sigma)
    =
    {\rm dim}({\rm Im}(B)) + {\rm dim}({\rm Ker}(B^T))
    =
    N .
\end{equation}
Since $Z = \mathbb{R}^N \times \mathbb{R}^N$, it follows that the constraint set $E$ is an affine subspace of $Z$ of dimension $N$ and co-dimension $N$.

This observation  characterizes the structure and dimensionality of the phase space $Z$ and of the subspace $E$ of all admissible states in $Z$, or constraint set, i.~e., the set of all the states that satisfy the conservation laws (\ref{vZftZT}). To wit, the phase space has the conventional dual structure $Z = \mathbb{R}^N \times \mathbb{R}^N$ of work-conjugate pairs, whereas the constraint set $E$ is a linear subspace of $Z$ of dimension $N$ and co-dimension $N$. This structure is particular to systems obeying conservation laws and sets them apart from other applications in data science where the data is heterogeneous and unstructured.

\subsection{Material characterization}

We assume that the material behavior is inherently random and characterized by a positive, continuous, likelihood function $L_D(y)$ representing the likelihood of observing a material state $y\in Z$, cf.~Fig.~\ref{uCXq6P}a. Since material behavior is not localized to a bounded region of phase space in general, the material likelihood function $L_D$ contains infinite mass and cannot be normalized to define a probability density. The function $\Phi_D(y) = - \log L_D(y)$ is the corresponding {\sl potential}.

\begin{example}[Local material behavior] \label{dzeXXX}{\rm
Material behavior is often {\sl local} and can be characterized over each member phase space $Z_e = \mathbb{R}^d \times \mathbb{R}^d$ by a local material likelihood function $L_{D,e} \in \mathcal{M}(Z_e)$. Assuming that the behavior of the members is independent, then the {\sl global} likelihood function is
\begin{equation}
    L_D = L_{D,1} \times \cdots \times L_{D,m} ,
\end{equation}
and the {\sl global} potential is
\begin{equation}
    \Phi_D = \Phi_{D,1} + \cdots + \Phi_{D,m} .
\end{equation}
Likelihood functions need not be integrable and do not define a probability density in general. For instance, suppose that the local behavior is characterized by a sliding Gaussian of the form, cf.~Fig.~\ref{uCXq6P}a.,
\begin{equation}
    L_{D,e}(y_e)
    =
    \exp\Big( - \frac{w_e}{2s^2} \| \sigma_e - \mathbb{C}_e \epsilon_e \|^2\Big)
\end{equation}
where $y_e = (\epsilon_e, \sigma_e) \in Z_e$ and the parameter $s>0$ measures the width of the distribution. The global likelihood function is, then,
\begin{equation}\label{nJFGU9}
    L_{D}(y)
    =
    \exp
    \Big(
        -
        \frac{1}{2s^2}
        \sum_{e=1}^m w_e \| \sigma_e - \mathbb{C}_e \epsilon_e \|^2
    \Big) ,
    \quad
    \| \sigma_e \|^2 = \mathbb{C}_e^{-1} \sigma_e \cdot \sigma_e ,
\end{equation}
$y = (\epsilon,\sigma)$. Equivalently,
\begin{equation}\label{0VUrvV}
    L_{D}(y)
    =
    \exp
    \Big(
        -
        \frac{1}{2s^2}
        \Big(
            \| \epsilon \|^2 + \| \sigma \|^2 - 2 \sigma \cdot \epsilon
        \Big)
    \Big) ,
\end{equation}
with
\begin{equation}
    \| \epsilon \|^2
    =
    \sum_{e=1}^m
         w_e \mathbb{C}_e \epsilon_e \cdot \epsilon_e ,
    \quad
    \| \sigma \|^2
    =
    \sum_{e=1}^m
         w_e \mathbb{C}_e^{-1} \sigma_e \cdot \sigma_e ,
    \quad
    \sigma \cdot \epsilon = \sum_{e=1}^m w_e \sigma_e \cdot \epsilon_e ,
\end{equation}
or
\begin{equation}
    L_{D}(y)
    =
    \exp
    \Big(
        -
        \frac{1}{2s^2}
        \Big(
            \| y \|^2 - 2 \sigma \cdot \epsilon
        \Big)
    \Big)
\end{equation}
in the norm (\ref{Poyet2}). Evidently, $L_{D,e}$ decays away from the centerline $D_e = \{ \sigma_e - \mathbb{C}_e \epsilon_e \}$ but is constant along the line. Hence,
\begin{equation}
    \int_{Z_e} L_{D,e}(y_e) \, dy_e = +\infty ,
\end{equation}
and neither $L_{D,e}(y_e)$ nor $L_D(y)$ can be normalized.
\hfill$\square$}
\end{example}

\subsection{Classical inference}

If the constraint set $E$ and the material likelihood function $L_D$ are fully known, , cf.~Fig.~\ref{uCXq6P}a, the expectation of any quantity of interest, represented by a bounded continuous function $f \in C_b(Z)$, is simply given by
\begin{equation}\label{K2mXSg}
    \mathbb{E}[f]
    =
    \frac{
        \int_E
            f(z) L_D(z)
        \, dz
    }
    {
        \int_E
            L_D(z)
        \, dz
    } ,
\end{equation}
which fully characterizes the response of the system in the sense of probability.

\begin{example}[Random trusses]\label{suO4Ci}
Consider a truss such as defined in Example~\ref{nBpq7d}, with material behavior characterized by the likelihood function (\ref{nJFGU9}). In addition, we metrize phase space by means of the Euclidean norm (\ref{Poyet2}). We begin by parameterizing the constraint space $E$. Recall that $m$ is the number of members of the truss, $Z = \mathbb{R}^m \times \mathbb{R}^m$ is the phase space and $n < m$ the number of unconstrained degrees of freedom. Let $l = m-n$ and $A \in \mathbb{R}^{m\times l}$ the matrix whose columns define a basis of ${\rm Ker}(B^T)$. Further, we define the matrix of Gauss point volumes as $\mathbb{W}={\rm diag}(w_1,\dots,w_m)$. Suppose non-degeneracy in the sense that $\sigma$ satisfies the equilibrium condition $B^T\mathbb{W}\sigma = 0$ if and only if there is $v \in \mathbb{R}^l$ such that
\begin{equation}\label{RYaYcU}
    \sigma = \mathbb{W}^{-1}A v .
\end{equation}
We may thus regard $v$ as a discrete Airy potential and $A$ as a discrete Airy operator. In addition, suppose that $\epsilon$ satisfies the compatibility condition $A^T \epsilon = 0$ if and only if there are displacements $u$ such that $\epsilon = B u$. The, the constraint set $E_0$ through the origin, corresponding to $f=0$ and $g=0$, admits the representation
\begin{equation}
    E_0
    =
    \{
    (\epsilon,\sigma) \in Z = \mathbb{R}^m \times \mathbb{R}^m
    \, : \,
    \epsilon = B u,\ u \in \mathbb{R}^n;
    \ \sigma = \mathbb{W}^{-1}A v,\ v \in \mathbb{R}^l \} .
\end{equation}
The general constraint set is then the translation
\begin{equation}
    E = z_0 + E_0,
    \quad
    z_0 = (\epsilon_0,\sigma_0) ,
    \quad
    \epsilon_0 = g,
    \quad
    B^T \mathbb{W}\sigma_0 = f .
\end{equation}
We verify that the dimension of $Z$ is $2N = 2m$ and the dimension of $E$ is $N = m = n + l$. We may regard $(u,v) \in \mathbb{R}^n \times \mathbb{R}^l$ as a set of coordinates parameterizing $E$. In this representation, the posterior likelihood function takes the form
\begin{equation}\label{nx9SKf}
    L(u,v)
    =
    L_D(Bu+\epsilon_0, \mathbb{W}^{-1}Av+\sigma_0) ,
\end{equation}
where we write $L_D(\epsilon,\sigma)$ for the material likelihood function. The corresponding posterior probability functions then follows by normalization of $L(u,v)$. For instance, if $L_D$ is of the form (\ref{0VUrvV}), representation (\ref{nx9SKf}) gives
\begin{equation}
    L(u,v)
    =
    \exp
    \Big(
        -
        \frac{1}{2s^2}
        \Big(
            \| Bu+\epsilon_0 \|^2 + \| \mathbb{W}^{-1}Av+\sigma_0 \|^2 - 2 (\mathbb{W}^{-1}Av+\sigma_0) \cdot (Bu+\epsilon_0)
        \Big)
    \Big) ,
\end{equation}
or, using the orthogonality properties of $A$ and $B$,
\begin{equation}
    L(u,v)
    =
    \exp
    \Big(
        -
        \frac{1}{2s^2}
        \Big(
            \| Bu+\epsilon_0 \|^2
            +
            \| \mathbb{W}^{-1}Av+\sigma_0 \|^2
            -
            2 \mathbb{W}^{-1}Av \cdot \epsilon_0
            -
            2 \sigma_0 \cdot Bu
        \Big)
    \Big) .
\end{equation}
From the properties of Gaussian integrals, the expected values $(\bar{u},\bar{v})$ of the coordinates $(u,v)$ are obtained by minimizing the potential $\Phi=-\log(L)$, with the result,
\begin{subequations}
\begin{align}
    &
    \bar{u} = (B^T\mathbb{W}\mathbb{C}B)^{-1} B^T\mathbb{W}(\sigma_0 - \mathbb{C} \epsilon_0) ,
    \\ &
    \bar{v} = (A^T\mathbb{W}^{-1}\mathbb{C}^{-1}A)^{-1} A^T(\epsilon_0 - \mathbb{C}^{-1} \sigma_0) ,
\end{align}
\end{subequations}
which are computed by inverting the stiffness and compliance matrices of the truss, $B^T \mathbb{W}\mathbb{C} B$ and $A^T\mathbb{W}^{-1}\mathbb{C}^{-1}A$, respectively, with $\mathbb{C}={\rm diag}(\mathbb{C}_1, \cdots, \mathbb{C}_m)$.
The posterior probability density then follows as
\begin{equation}
    {L}(u,v)
    =
    \sqrt{\det(B^T \mathbb{W}\mathbb{C} B /2\pi s^2 )}
    \sqrt{\det(A^T \mathbb{W}^{-1}\mathbb{C}^{-1} A /2\pi s^2 )}
    \times
\end{equation}
\begin{equation*}
    \exp
    \Big(
        -
        \frac{1}{2s^2}
        \Big(
            (B^T\mathbb{W}\mathbb{C}B)(u-\bar{u}) \cdot (u-\bar{u})
            +
            (A^T\mathbb{W}^{-1}\mathbb{C}^{-1}A)(v-\bar{v}) \cdot (v-\bar{v})
        \Big)
    \Big) ,
\end{equation*}
which provides a full account of the probability of outcomes. \hfill$\square$
\end{example}

\subsection{Approximation by empirical data} \label{8yK9TU}

Suppose, contrariwise, that, as in commonly the case in practice, the material likelihood function $L_D$ is known only approximately through a sequence of empirical samples, e.~g., drawn from experimental measurements, consisting of point-data sets $P_h = \{ y_{h,i},\ i=1,\dots,M_h \} \subset Z$, $h=0,1,\dots$, of size $M_h$, cf.~Fig.~\ref{uCXq6P}b. We additionally associate with each discrete data point a certain likelihood $c_{h,i} \in [0,1]$ accounting, e.~g., for experimental reliability and scatter \cite{kirchdoerfer2017data}. Thus, $c_{h,i} = 1$ for a data point $y_{h,i}$ that is sure to be free of experimental error, whereas $c_{h,i} < 1$ if the experimental measurements show scatter or are of questionable provenance. For every material point-set $P_h$, we consider discrete expectations of the general form
\begin{equation}\label{a2VR0G}
    \mathbb{E}_h[f]
    =
    \sum_{i=1}^{M_h}
    \Big( \int p_{h,i}(z) f(z) \, dz \Big),
\end{equation}
where $f \in C_b(Z)$ is a continuous function representing a quantity of interest,
and we require
\begin{equation}\label{fREO9W}
    p_{h,i}(z) \geq 0 ,
    \quad
    \sum_{i=1}^{M_h}
    \Big( \int_E p_{h,i}(z) \, dz \Big)
    =
    1 .
\end{equation}
We envision an experimental campaign resulting in point-data samples $P_h$ of increasing size, $M_h \uparrow \infty$ and increasing fidelity $c_{h,i} \uparrow 1$. The central question is, then, under what conditions on the data, if any, and what choices of probabilities $p_{h,i}(z)$ the approximate expectations $\mathbb{E}_h[f]$ {\sl converge} to the exact ones $\mathbb{E}[f]$ in the limit as $h\to\infty$.

We specifically consider probabilities $p_h \equiv \{p_{h,i}(z),\ i=1,\dots,M_h\}$ that minimize the regularized Kullback-Leibler (KL) divergence \cite{Kullback:1951, Kullback:1959}
\begin{equation}\label{Qg3xQI}
    G_\beta(p_h)
    =
    \sum_{i=1}^{M_h}
    \Big(
        \int_E
            \Big(
                \beta \| y_{h,i} - z \|^2
                +
                \log\dfrac{p_{h,i}(z)}{c_{h,i}}
                -
                N \log\frac{\beta}{\beta_0}
            \Big)
            \, p_{h,i}(z)
        \, dz
    \Big) ,
\end{equation}
for some reference $\beta_0 > 0$.
Evidently, the divergence function (\ref{Qg3xQI}) penalizes distance to the constraint set $E$, i.~e., it assigns an increasing discrepancy to material data points $y_{h,i}$ that are far from $E$, while simultaneously attempting to match the a priori likelihood $c_{h,i}$ of the point. The optimal regularized discrete probability $p_h$ then follows by minimization of $G_\beta$, with the explicit result
\begin{equation}\label{zNR2jN}
    p_{h,i}(z)
    =
    \frac
    {
        c_{h,i}
        {\rm e}^{-\beta\|y_{h,i}-z\|^2}
    }
    {
        \sum_{i=1}^{M_h} \int_E c_{h,i}
        {\rm e}^{-\beta\|y_{h,i}-z\|^2} \, dz
    } ,
    \qquad
    z \in E .
\end{equation}
We see from this expression that material data points $y_{h,i}$ that are far from the constraint set $E$ contribute less than points are are close, as expected. The parameter $\beta$ sets the scale of the distance as $1/\sqrt{\beta}$.

\subsection{Annealing schedule}

We note that, for $\beta \to +\infty$, the posterior probabilities (\ref{zNR2jN}) degenerate in general. This degeneracy simply reflects the fact that the point sets $P_h$ are unlikely to have a consequential intersection with the admissible set $E$, i.~e., the likelihood that a point $y_{h,i}$ be admissible with respect to a given deterministic loading is zero. As $\beta \to 0$, the posterior probabilities (\ref{zNR2jN}) also degenerate. There must therefore be an optimal value of $\beta$ where the probabilities are least degenerate. Conveniently, this optimal value $\beta_h$ of $\beta$ again follows by minimization of the KL divergence (\ref{Qg3xQI}). Thus, evaluating $G_\beta(p_h)$ at the minimizer (\ref{zNR2jN}), we obtain
\begin{equation}
    \min G_\beta
    =
    -
    \log
    \Big(
        \sum_{i=1}^{M_h} \int_E c_{h,i}
        {\rm e}^{-\beta\|y_{h,i}-z\|^2} \, dz
    \Big)
    -
    N \log \frac{\beta}{\beta_0} .
\end{equation}
Minimizing with respect to $\beta$ gives the optimality condition
\begin{equation}\label{cDTXAV}
    \frac{1}{\beta_h}
    =
    \frac{1}{N}
    \frac
    {
        \sum_{i=1}^{M_h}
        \int_E
            c_{h,i}
            \|y_{h,i}-z\|^2
            {\rm e}^{-\beta_h \|y_{h,i}-z\|^2}
        \, dz
    }
    {
        \sum_{i=1}^{M_h}
        \int_E
            c_{h,i}
            {\rm e}^{-\beta_h \|y_{h,i}-z\|^2}
        \, dz
    } ,
\end{equation}
which defines $\beta_h$ implicitly. Finally, inserting (\ref{cDTXAV}) into (\ref{a2VR0G}) we obtain
\begin{equation}\label{Pa55uo}
    \mathbb{E}_h[f]
    =
    \frac
    {
        \sum_{i=1}^{M_h}
        \Big(
            \int_E
                c_{h,i}
                {\rm e}^{-\beta_h \|y_{h,i}-z\|^2}
                f(z)
            \, dz
        \Big)
    }
    {
        \sum_{i=1}^{M_h}
        \Big(
            \int_E
                c_{h,i}
                {\rm e}^{-\beta_h \|y_{h,i}-z\|^2}
            \, dz
        \Big)
    } ,
\end{equation}
which supplies an expression for the approximate expectations that is explicit up to quadratures.

\subsection{Factorization and localization}

By exchanging the integrals over $E$ and the sums over data in (\ref{Pa55uo}), the approximate expectations have the equivalent form
\begin{equation}\label{tQ6bZL}
\begin{split}
    \mathbb{E}_h[f]
    =
    \frac
    {
        \int_E
        \sum_{i=1}^{M_h}
            f(z)
            {\rm e}^{-\beta_h \|y_{h,i}-z\|^2}
            c_{h,i}
        \, dz
    }
    {
        \int_E
        \sum_{i=1}^{M_h}
            {\rm e}^{-\beta_h \|y_{h,i}-z\|^2}
            c_{h,i}
        \, dz
    } .
\end{split}
\end{equation}
The great advantage of this alternative form is that it factorizes with respect to local data, which eliminates the combinatorial complexity of the first form. Thus, suppose that
\begin{equation}
    D_h = \prod_{e=1}^m D_{e,h} ,
\end{equation}
where $D_{e,h} = \{ y_{e,h,i},\ i=1,\dots,M_{e,h} \}$  are local data files. Then,
\begin{equation}\label{0lgoUa}
    \int_E
    \sum_{i=1}^{M_h}
        f(z)
        {\rm e}^{-\beta_h \|y_{h,i}-z\|^2}
        c_{h,i}
    \, dz
    =
    \int_E
        f(z) {L}_h(z)
    \, dz ,
\end{equation}
where
\begin{equation}
    {L}_h(z)
    =
    \sum_{i=1}^{M_h}
        {\rm e}^{-\beta_h \|y_{h,i}-z\|^2}
        c_{h,i}
\end{equation}
is an inferred likelihood density. For local data, the sum over data factorizes, with the result
\begin{equation}
    {L}_h(z)
    =
    \prod_{e=1}^m {L}_{e,h}(z_e) ,
\end{equation}
where
\begin{equation}\label{6r8WgM}
    {L}_{e,h}(z_e)
    =
    \sum_{i=1}^{M_{e,h}}
        {\rm e}^{-\beta_h \|y_{e,h,i}-z_e\|^2}
        c_{e,h,i}
\end{equation}
are local likelihood density functions. We note that the computational complexity of the expectation is $\sum_{e=1}^m M_{e,h}$, instead of $\prod_{e=1}^m M_{e,h}$, as before, which renders the evaluation of approximate expectations computationally feasible.

We further note that, for local data, we can generalize (\ref{6r8WgM}) to
\begin{equation}
    {L}_{e,h}(z_e)
    =
    \sum_{i=1}^{M_{e,h}}
        {\rm e}^{-\beta_{e,h} \|y_{e,h,i}-z_e\|^2}
        c_{e,h,i} ,
\end{equation}
where we localize temperature. We can then estimate $1/\beta_{e,h}$ as the average minimum distance squared between points in $D_{e,h}$, i.~e.,
\begin{equation}\label{eq:beta_avg}
    \frac{1}{\beta_{e,h}}
    \sim
    \frac{1}{M_{e,h}}
    \sum_{i=1}^{M_{e,h}}
    \min_{j \neq i} \| y_{e,h,i} - y_{e,h,j} \|^2 .
\end{equation}
Finally, the integrals over $E$ can be conveniently computed by stochastic quadrature. These and other matters of implementation are discussed next.

\section{Numerical implementation}

\subsection{Problem statement}
In evaluating (\ref{tQ6bZL}), two main operations need to be performed: the integration over the constraint set $E$ and the summation over the global material data set for the evaluation of the material likelihood. The main computational challenges pertaining to these operations are:
\begin{enumerate}
    \item The evaluation of the local material likelihood involves a summation over local data sets for each material point that can be computationally demanding. However, most data points have negligible influence on the total sum, which suggests the use of efficient, carefully designed {\sl a propos} data structures.
    \item The Monte Carlo integration over the constraint set $E$ requires the generation of random points on $E$, which requires the repeated solution of global systems of equations if the sampling is done by projection. It is therefore advantageous to identify sampling strategies that reduce the number of projections.
    \item The thermalized likelihood function implied by (\ref{Pa55uo}) is exceedingly multi-modal, featuring peaks at every data point (cf.~\cite{kanno2019mixed} for a discussion in the context of {\sl min-dist} DD). This intricate likelihood landscape arises even if the underlying likelihood function is smooth and compounds the evaluation of the integrals over $E$ in (\ref{Pa55uo}).
\end{enumerate}
These challenges are addressed next in turn.

\subsection{Population annealing}
The complexity of the thermalized likelihood landscape renders the direct application of standard Markov chain importance sampling strategies, such as the classical Metropolis--Hastings algorithm \cite{metropolis1953equation}, difficult. Therefore, a more tailored sampling strategy is required, as discussed next.

\subsubsection{Population Annealing}
{\sl Population Annealing} (PA) \cite{iba2001population, machta2010population, weigel2021understanding} is a sequential Monte Carlo method designed for complex likelihood functions and scalable parallel computations. A set of states, the population, is maintained in order to explore different regions of phase space that are possibly separated by energy barriers resulting from local minima. Each population member performs a relatively short Markov chain importance sampling. In addition, the population is resampled by a combination of pruning and enriching, in analogy to {\sl Genetic Algorithms}. The process is driven by a {\sl Simulated Annealing} schedule, which slowly increases the inverse local temperature $\beta_e$. For every new temperature, the population is re-equilibrated by resampling and a specified number of Markov chain moves. In this manner, at low $\beta_e$ the population can escape readily from local minima and perform a more exhaustive exploration. At high $\beta_e$, the population is increasingly driven towards the constraint set $E$.

\subsubsection{Initialization and resampling step}
The principle idea of the resampling step is to remove and duplicate the population samples according to their energy while increasing $\beta_e$ in order to accelerate the convergence to low-energy states. Resampling is done in such a way that the population size fluctuates around a target population size $N_P^*$ and equilibrium of the distribution is maintained for every new temperature.

Initially, at $\beta_e=0$, a population $\mathcal{P} = \{z_p\in E \, : \, p=1\dots N_{P}\}$ of $N_{P}$ states is generated either randomly by the closest-point projection of random global data points onto $E$ or using a {\sl min-dist} DD solver \cite{kirchdoerfer2016data}. At initialization, the population size $N_P$ may be chosen equal to the target size $N_P^*$, or set to a different value, in which case the size gets adapted to the target size during resampling. The energy of the population is initially set to infinity. For $\beta_e>0$, the energy is defined as
\begin{equation}\label{eq:energy}
    e_p
    =
    -
    \frac{1}{\beta_e} \,
    \sum_{e=1}^{m}
        \log
        \Big(
            \frac{1}{M_e}
            \sum_{i=1}^{M_e}
            {\rm e}^{-\beta_e\|y_{e,i}-z_{e,p}\|^2}
        \Big).
\end{equation}
where we assume local data. The efficient evaluation of (\ref{eq:energy}) is discussed in Section \ref{sec:energy}.

At each iteration, or {\sl quench}, $\beta_e$ is increased according to
\begin{equation}
    \beta_e = \beta_e+\Delta\beta_e
\end{equation}
for a specified number of quenches $N_Q$. For simplicity, in calculations we choose a constant increment
\begin{equation}
    \Delta\beta_e = \beta_{e,f}/N_Q
\end{equation}
for a pre-computed target $\beta_{e,f}$. The simple choice of $\beta_{e,f}$ given by the inverse average nearest neighbor distance (\ref{eq:beta_avg}) is used in the numerical examples presented in Sections~\ref{sec:three_bars_Gaussian} to \ref{sec:pactruss_study}.

For resampling, the relative weights
\begin{equation}
    \hat\tau_p = N_P^*\frac{{\rm e}^{-\Delta\beta_e e_p}}{\sum_{i=1}^{N_{P}}{\rm e}^{-\Delta\beta_e e_i}}
\end{equation}
are computed for each population member. The required number of copies of population member $p$ is then conveniently estimated as \cite{wang2015comparing,barash2017gpu}
\begin{align}
    l_p
    =
    \begin{cases}
        \lfloor\hat\tau_p\rfloor, &  \text{if}~ u>\hat\tau_p-\lfloor\hat\tau_p\rfloor ,\\
        \lfloor\hat\tau_p\rfloor+1, & \text{otherwise} ,
    \end{cases}
\end{align}
where $u$ is a random nuber drawn uniformly from $[0,1)$ and the floor function $\lfloor x\rfloor$ is the largest integer less or equal $x$. If $l_p=0$, member $p$ is removed from the population. The resampling step incurs a modest computational overhead that is offset by the overall performance gain resulting from reconfiguring the population to the relevant regions of the phase space.

\subsubsection{Markov chain importance sampling step}\label{sec:pca}
Resampling introduces a bias since population members may be duplicated. In order to eliminate this bias, a Markov chain importance sampling step is performed. In this step, every population member executes a random walk, thereby exploring new states and decorrelating the population. The random walk consists of repeated random moves in the constraint set $E$. We recall that the constraint set $E$ is a hyper-plane of dimension $N=md$ (number of material points times dimension of stress/strain) in an $2N-$dimensional global phase space $Z$. We can conveniently exploit this affine structure and represent $E$ in terms of an $N$-dimensional basis.

We compute the requisite basis by means of a principal component analysis (PCA). To this end, we generate $K>N$ random global data points $z_i \in Z$ and project them onto $E$. For consistency with the norm (\ref{Poyet2}), we further weight the strains and stresses by $\sqrt{w_e}\mathbb{C}_e^{1/2}$ and $\sqrt{w_e}\mathbb{C}_e^{-1/2}$, respectively. The weighted states $z_i^w$ (here and subsequently the superscript $w$ indicates weighting) are collected in a matrix $T\in\mathbb{R}^{K\times 2N}$ in order to compute its covariance matrix $C_E\in\mathbb{R}^{2N\times 2N}$. Computing the eigenvectors of $C_E$ associated to the $N$ largest eigenvalues then gives the basis vectors of the constraint set as a matrix $A_E\in\mathbb{R}^{2N\times N}$.

To effect a move, $N$ random numbers are drawn from a standard normal distribution and multiplied by a step size $s_p$, transformed back by $A_E^T$ and added to the weighted state $z^w_p$
\begin{equation}
    z^w_{p,\rm trial} = z^w_p + s_p~ g~A_E^T, \quad g\sim\mathcal{N}(0,I_N) .
\end{equation}
The proposed move to the trial state $z_{p,\rm trial}$ (obtained by reverse weighting of $z^w_{p,\rm trial}$) is then either accepted or rejected according to the Metropolis-Hastings acceptance probability
\begin{equation}
    P_{\beta_e}=\min(1,\exp(-\beta_e(e_{p,\rm trial}-e_p))) ,
\end{equation}
where $e_{p,\rm trial}$ is the energy of $z_{p,\rm trial}$. This acceptance criterion allows moves that increase the energy slightly, in particular at the initial stages of annealing characterized by low $\beta_e$. In addition, the step size $s_p$ influences the acceptance rate and, therefore, the efficiency of the algorithm. In calculations, we chose the simple adaption scheme
\begin{equation}
    s_p = s_p + (r_p-r^*)s_p ,
\end{equation}
with a target acceptance rate $r^*\in[0,1]$. The acceptance rate $r_p$ is measured as the ratio of accepted moves to the total number of trials $N_T$. The adaption scheme has the effect of reducing the step size if the acceptance rate is smaller than $r^*$, which signals an overly-aggressive exploration strategy, and increasing it if the acceptance rate is too small, which contrariwise signals an overly-timid exploration strategy.

\subsubsection{Energy computation}\label{sec:energy}
The complexity of the evaluation of the energy (\ref{eq:energy}) is $\sum_e M_e$. However,  for sufficiently large $\beta_e$ a vast majority of the data points contribute neglibly to(\ref{eq:energy}), which suggests restricting the sum with the aid of supporting data structures. In \cite{eggersmann2021efficient}, a number of data structures for local data sets were investigated in the framework of min-dist DD solvers and the hierarchical k-means tree \cite{fukunaga1975branch} was found to be particular well-suited for executing approximate nearest neighbor searches.

In this work, we likewise use the hierarchical $k$-means implementation of the FLANN library \cite{muja2014scalable}. The approximate nearest neighbor search provided by the library is performed using a priority queue of maximum length $N_{\rm checks}$ in order to reduce backtracking operations. We employ a radius search technique, which returns all data points $y_{e,i}\in D_e$ within a specified radius
\begin{equation}
    r^2_{\rm TOL}=-\frac{1}{\beta_e}\log({\rm TOL}).
\end{equation}
around the query point $z_{p,e}$. This technique has the advantage over a $k$-nearest-neighbor ($k$NN) search that the number of nearest neighbors $k$ need not be specified. Instead, the sum is restricted to points whose contribution exceeds the tolerance ${\rm TOL}$.

Our numerical tests reveal that using a radius search increases performance markedly and that, remarkably, the performance gains are insensitive to the choice of tolerance. Thus, in tests values of ${\rm TOL} = $  $10^{-8},10^{-16},10^{-32}$ delivered nearly identical performance. In view of this tradeoff, in all calculations presented subsequently we use ${\rm TOL}=10^{-16}$. By contrast, the performance of the algorithm is strongly dependent on the degree of backtracking as set by $N_{\rm checks}$, cf.~Section \ref{sec:pactruss_study}.

It bears emphasis that only one material data structure needs to be set up for each specific material in the structure, since all members made of the same material access the same local data set. Likewise, only one inverse local temperature $\beta_{e,f}$ and one quenching scheme is required for each specific material. In particular, if all structural members are made of the same material, only one material data structure and quenching scheme need to be maintained and run.

\subsubsection{Summary of the algorithm}
The described algorithm is summarized in Algorithm \ref{alg:PA}. It returns the final population $\mathcal{P}$ of size $N_P$. The expectation of a quantity of interest $f$ is evaluated as
\begin{equation}
    \mathbb{E}[f]=\frac{1}{N_{P}}\sum_{p=1}^{N_{P}} f(z_p), \quad z_p\in \mathcal{P}.
\end{equation}
One advantage of the algorithm is that does not require extensive fine-tuning of the input variables. Thus, the target population size $N_{P}^*$, the number of quenches $N_Q$ and the number of Markov chain moves $N_T$ are parameters that are chosen as large as permitted by computational resources. In all calculations, the target acceptance rate is set to $r^*=0.25$. The initial step size $s_{0}$ influences only the first few iterations as it is adapted during the quenching schedule. The step size tens to decrease with increasing dimensionality \cite{graves2011automatic}.
Lastly, we note that the population member perform their moves independently of each other, which renders the algorithm {\sl embarrassingly parallel}.

\begin{algorithm}[H]
\caption{Population annealing algorithm}
\label{alg:PA}
\begin{algorithmic}
\REQUIRE Constraint set $E \subset Z$, basis vectors $A_E$.
\REQUIRE Local data set(s) $D_e$
\REQUIRE $N_Q,N_T,N_{P}^*,s_{0},r^*,\beta_{e,f}$
\STATE Initialize: $\mathcal{P}\leftarrow\{z_p\in E\}$ randomly or heuristically
\STATE Initialize: step size $s_p\leftarrow s_{0}$
\STATE Initialize: energy $e_p\leftarrow \infty$
\STATE Initialize: $\beta_e\leftarrow0$
\STATE Set: $\Delta\beta_e\leftarrow\beta_{e,f}/N_Q$

\FOR{$q=0,\dots,N_Q-1$}
\STATE Set: $\beta_e\leftarrow\beta_e+\Delta\beta_e$
\STATE Set: $\mathcal{P}\leftarrow\texttt{Resampling}(\mathcal{P},\Delta\beta_e)$

\FOR{\textbf{each} $z_p\in \mathcal{P}$}
\FOR{$t=0,\dots,N_T-1$}
\STATE Set: $z_{p,\rm trial}\leftarrow \texttt{RandomMove}(z_p,s_p,A_E)$
\STATE Set: $e_{p,\rm trial}\leftarrow \texttt{Energy}(D_e,\beta_e,z_{p,\rm trial})$

\STATE Set: $u\leftarrow \texttt{UniformRand}(0,1)$
\STATE Set: $w \leftarrow \exp(- \beta_e(e_{p,\rm trial} -e_p))$
\IF {$u\leq w$}

\STATE Set: $z_p\leftarrow z_{p,\rm trial}$
\STATE Set: $e_p\leftarrow e_{p,\rm trial}$
\ENDIF
\ENDFOR
\ENDFOR
\STATE Set: $r_p\leftarrow\texttt{AcceptanceRate}$
\STATE Set: $s_p\leftarrow s_p+(r_p-r^*)s_p$
\ENDFOR

\RETURN $\mathcal{P}$
\end{algorithmic}
\pagebreak
\end{algorithm}

\section{Examples of application}

We conclude with selected examples of application that illustrate simply the range, scope and convergence properties of the approach.

\subsection{Three-bar truss with sliding-Gaussian material data set}\label{sec:three_bars_Gaussian}

We begin with a simple verification example concerned with the three-bar truss shown in Figure \ref{fig:truss_gauss}a loaded by a force $P$ and undergoing a corresponding displacement $\Delta$. In calculations, we set $P=100$, $L=1$ and all cross-sectional areas $A=1$ (here and subsequently units are omitted). The outcome quantity of interest is chosen to be the displacement $\Delta$.

In order to have an analytical reference solution, we assume that the material data set obeys a sliding-Gaussian distribution, in which case the exact solution of the inference problem is given in Section~\ref{suO4Ci}. The local material data set is generated according to Section \ref{dzeXXX} with $\mathbb{C}_e=10000$ and standard deviation $s=5 \times 10^{-4}$. A typical data set is shown in Fig.~\ref{fig:truss_gauss}b. We emphasize that such material data sets are the only input to the calculations and, once generated, the underlying distribution whence they are sampled is jettisoned altogether and presumed unknown.

\begin{figure}[ht]
	\begin{subfigure}{0.45\textwidth}\caption{} \includegraphics[width=0.99\linewidth]{Figures/HW8_Problem1_Figure1.pdf}
	\end{subfigure}
    $\quad$
	\begin{subfigure}{0.40\textwidth}\caption{} \includegraphics[width=0.99\linewidth]{Figures/dataset_1000b.pdf}
	\end{subfigure}
    \caption{a) Truss with prescribed applied force $P$ undergoing a displacement $\Delta$. b) Local material data set with $1000$ data points sampled from a sliding-Gaussian likelihood function.} \label{fig:truss_gauss}
\end{figure}

Calculations are carried out using Algorithm~\ref{alg:PA} with number of quenches $N_Q=100$, number of trials $N_T=20$ and initial step size $s_{0}=1.0$. Nearest-neighbor searches with radius ${\rm TOL}=10^{-16}$ are performed exactly, i.~e., with an unlimited number of backtracking checks $N_{\rm checks}$. As a measure of the error we choose the Kolmogorov–Smirnov statistic $e_{\rm KS}$, which measures the maximum difference between the exact and computed cumulative distributions.

\begin{figure}[ht]
	\begin{subfigure}{0.30\textwidth}\caption{} \includegraphics[width=0.99\linewidth]{Figures/data_dist_0b.pdf}
	\end{subfigure}
	\begin{subfigure}{0.30\textwidth}\caption{} \includegraphics[width=0.99\linewidth]{Figures/data_dist_1.pdf}
	\end{subfigure}
    	\begin{subfigure}{0.30\textwidth}\caption{} \includegraphics[width=0.99\linewidth]{Figures/data_dist_2.pdf}
	\end{subfigure}
    \caption{Three-bar truss with sliding-Gaussian material data. Computed histograms {\sl vs.} exact distribution (red) of displacement $\Delta$ for material data-sets of sizes: a) $M_e=10^3$; b) $M_e=10^4$; c) $M_e=10^5$.} \label{fig:gauss_pdf_study}
\end{figure}

Fig.~\ref{fig:gauss_pdf_study} collects three computed histograms of displacement $\Delta$ for three different material data-set sizes $\{10^3,10^4,10^5\}$. The exact Gaussian distribution is shown in red for comparison. A clear trend towards convergence of the computed histogram to the exact distribution is apparent from the figure. It bears emphasis that, as expected, both the exact distribution and the computed histograms of displacement are strictly {\sl unimodal}. Unimodality is indicative of a lack of complexity of the behavior of a system. In particular, the maximum-likelihood state of the system---equivalently, the minimum free-energy solution---is uniquely defined \cite{kirchdoerfer2017data} and the outcomes of the system cluster around the maximum-likelihood state.

\begin{figure}[ht]
	\begin{subfigure}{0.45\textwidth}\caption{}
	\centering
	\includegraphics[height=5cm]{Figures/study_beta_2.pdf}
	\end{subfigure}

	\begin{subfigure}{0.45\textwidth}\caption{}
	\centering
	\includegraphics[height=5cm]{Figures/data_conv2.pdf}
	\end{subfigure}

	\begin{subfigure}{0.45\textwidth}\caption{}
	\centering
	\includegraphics[height=5cm]{Figures/population_conv.pdf}
	\end{subfigure}
    \caption {Three-bar truss with sliding-Gaussian material data. Convergence of temperatures and computed histograms of displacement $\Delta$. a) Evolution of $\beta_e$ during annealing process. b) Kolmogorov–Smirnov (KS) error {\sl vs.} material data-set size, $N_P=10^4$. c) KS error {\sl vs.} population size, $M_e=10^5$. The lines show the mean of $10$ data set realizations, the shaded area indicate the maximal deviation from mean.} \label{fig:data_set_error}
\end{figure}

Fig.~\ref{fig:data_set_error}a shows the evolution of $\beta_e$ during the annealing process. As may be seen from the figure, $\beta_e$ increases monotonically with the number of quenches and tends towards a final value. This final value itself increases with the material data-set size as the average minimum distance between material data points decreases. Fig.~\ref{fig:data_set_error}b shows the dependence of the Kolmogorov–Smirnov error $e_{KS}$ on the number of material data points, which again exhibits a clear trend towards convergence. An analysis of the convergence plot reveals that, asymptotically for large material data sets, $e_{KS} \sim M_e^{-\alpha}$, where $\alpha \sim 0.36$ is the rate of convergence. These results are consistent with quantitative error estimates derived from analysis \cite{Conti:2022}. Fig.~\ref{fig:data_set_error}c shows the dependence of $e_{KS}$ on the population size $N_P$.

\subsection{Three-bar truss with Weibull tensile-strength distribution}\label{sec:three_bars_Weibull}

Next, we present a simple example that is designed to illustrate the ability of the method to deal effectively with complex data and behavior. The example concerns the same structure as in Section~\ref{sec:three_bars_Gaussian}, Fig.~\ref{fig:dataset_fail}a, but the material is now assumed to be brittle with random tensile strength.

\begin{figure}[ht]
	\begin{subfigure}{0.45\textwidth}\caption{} \includegraphics[width=0.99\linewidth]{Figures/HW8_Problem1_Figure1.pdf}
	\end{subfigure}
    $\quad$
	\begin{subfigure}{0.40\textwidth}\caption{} \includegraphics[width=0.99\linewidth]{Figures/dataset_failureb.pdf}
	\end{subfigure}
    \caption{a) Truss with prescribed displacement $\Delta$ and reaction force $P$. b) Local material data set with $1000$ data points sampled from a bimodal likelihood function with failed and unfailed branches selected according to Weibull strength statistics, plus superimposed Gaussian noise.} \label{fig:dataset_fail}
\end{figure}

The control parameter is now the displacement $\Delta$, which is increased monotonically. The outcome of interest is the corresponding the reaction force $P$. By virtue of the monotonicity of $\Delta$, the strains in all bars also increase monotonically. The behavior of the bars in compression is linear elastic. By contrast, in bars in tension the stress increases linearly up to the tensile strength $\sigma_t > 0$ of the material and subsequently drops to zero at larger strains. In addition, the tensile strength of the material is assumed to be random and to obey Weibull statistics, with distribution
\begin{equation}
    W(\sigma_t)
    =
    1-{\rm e}^{-(\sigma_t/\sigma_0)^p} ,
\end{equation}
where $\sigma_0$ and $p$ are material-specific parameters. Under these conditions, the likelihood of observing a local tensile state $(\epsilon_e,\sigma_e)$ is {\sl bimodal} and consists of: a branch $\sigma_e = \mathbb{C}_e \epsilon_e$  with likelihood $1-W(\mathbb{C}_e \epsilon_e)$, corresponding to an unfailed bar; and a branch $\sigma_e = 0$ with likelihood $W(\mathbb{C}_e \epsilon_e)$, corresponding to a failed bar.

A typical material data set sampled from this likelihood function, with superimposed Gaussian noise of standard deviation $s=10^{-4}$, is shown in Fig.~\ref{fig:dataset_fail}b. The data set is intended to mimic empirical data sets such as might be generated experimentally, including experimental scatter. Again, we emphasize that such material data sets are the only input to the calculations and, once generated, the underlying distribution whence they are sampled is jettisoned altogether and presumed unknown.

\begin{figure}[ht]
    \centering
	\begin{subfigure}{0.45\textwidth}\caption{}
    \centering
    \includegraphics[height=5.0cm]{Figures/rf.pdf}
	\end{subfigure}%

	\begin{subfigure}{0.45\textwidth}\caption{}
    \centering
    \includegraphics[height=5.0cm]{Figures/rf_dist.pdf}
	\end{subfigure}
    \caption{Three-bar truss with Weibull tensile strength data. a) Computed mean and standard deviation of the reaction force $P$ for prescribed displacement $\Delta$. b) Computed distribution of $P$ at $\Delta=2\times 10^{-2}$.} \label{fig:rf}
\end{figure}

For the simulation parameters, we choose $N_Q=50, N_T=10$ and $N_P=2\times 10^4$ with data-set size of $M_e=10^3$. Fig.~\ref{fig:rf}a shows the computed mean and standard deviation of the reaction force $P$ over a range of  prescribed displacements $\Delta$. Interestingly, the figure reveals three regimes: i) and initial linear elastic response with low uncertainty, as measured by the standard deviation; ii) an intermediate regime combining failed and unfailed bars characterized by a decreasing reaction and high uncertainty; and iii) full failure, or {\sl collapse}, with low uncertainty.

Fig.~\ref{fig:rf}b in turn shows the computed distribution of $P$ at $\Delta=2\times 10^{-2}$. Remarkably, the distribution is strongly {\sl multimodal}, a telltale hallmark of complexity. The various modes correspond to several different possible combinations of failed and unfailed bars and their attendant likelihoods. For the data used in calculations, the right bar is computed to remain unfailed with probability $1$ or, equivalently, the probability of the right bar failing is computed to be $0$. By contrast, the remaining two bars can be failed or unfailed with varying likelihoods: no failure (right peak), double failure (left peak), single failures (two middle peaks). From a global perspective, the various modes in the posterior likelihood function of $P$ correspond to a {\sl clustering}, or {\sl classification}, of global material data points close to the constraint set. The ability of the method to classify outcomes automatically and effectively bears remark.

\subsection{Lightweight space structure}\label{sec:pactruss_study}
By way of demonstration, we consider the lightweight structure shown in Fig.~\ref{fig:pactruss} designed to be part of a modular space telescope \cite{lee2016architecture}. We specifically aim to assess the feasibility of the proposed approach for higher-dimensional problems (constraint set dimension $N=39$) and its computational performance in terms of resource requirements and execution times.

\begin{figure}[ht]
    \includegraphics[width=0.40\linewidth]{Figures/pactruss.pdf}
    \caption{Lightweight space structure \cite{lee2016architecture}. The red arrow indicates the quantity of interest (joint eccentricity $u_e$) to be inferred from the analysis.} \label{fig:pactruss}
\end{figure}

In calculations, we set the cross-sectional area of the bars to $A=0.1$ and the length of the non-diagonal bars to $1$. The structure is constrained throughout its base and is loaded by longitudinal joint forces $P=5$ as shown in Fig.~\ref{fig:pactruss}. As a quantity of interest, we choose the top central joint eccentricity $u_e=\sqrt{u_x^2+u_y^2}$ whereby $u_x$ and $u_y$ are the in-plane displacements of the joint. We specifically consider a material data set with $M_e=10^4$ data points sampled from a sliding Gaussian as in Section~\ref{sec:three_bars_Gaussian}, cf.~Fig.~ \ref{fig:truss_gauss}. Conveniently, this assumption affords a closed-form analytical solution for the distribution of $u_e$ that can be used by way of reference, cf.~Section~\ref{suO4Ci}.

For algorithm parameters, we choose $N_P=2\times 10^3, N_T=10, N_Q=100$ and $ s_{0}=0.01$. For faster convergence to low-energy states, we initialize the population by means of $100$ random min-dist solutions \cite{kirchdoerfer2016data} for noise-free data. As expected, the posterior distribution derived from the initial population is far from the exact distribution, Fig.~\ref{fig:pactruss_pdf_study}a. By contrast, a close agreement is achieved following the application of the population annealing algorithm, Fig.~\ref{fig:pactruss_pdf_study}b, which is remarkable in view of the relatively small population and material-data sizes used in the calculations.

\begin{figure}[ht]
	\begin{subfigure}{0.5\textwidth}\caption{} \includegraphics[height=4.8cm]{Figures/pactruss_mindist_1000.pdf}
	\end{subfigure}%
	\begin{subfigure}{0.5\textwidth}\caption{} \includegraphics[height=4.8cm]{Figures/pactruss_distr_2000.pdf}
	\end{subfigure}
    \caption {Lightweight space structure \cite{lee2016architecture}. Computed posterior distributions of the joint eccentricity $u_e$ and exact reference distribution (red curves). a) Posterior distribution obtained from an initial population of $1000$ random min-dist solutions. b) Posterior distribution following application of the population annealing algorithm. }
    \label{fig:pactruss_pdf_study}
\end{figure}

Finally, we assess the impact of the approximate nearest neighbor search on accuracy and CPU times for the population size $N_P=2\times 10^3$. The absolute CPU times are, of course, implementation and machine dependent. In this work, the computations were performed on a 12-Core AMD Ryzen 9 3900X machine with Python bindings \cite{jakob2019pybind11} to a C++ implementation using Standard Template Library Parallel Algorithms.

The CPU times for a direct evaluation of Eq. (\ref{eq:energy}) without any data structuring or acceleration is $\sim 662$ seconds. The radius search with full backtracking terminates in $\sim 92$ seconds, or a seven-fold acceleration. Fig.~\ref{fig:pactruss_cpu_study} shows the effect of the degree of backtracking measured in terms of $N_{\rm checks}$. We see from the figure that, whereas the computation time increases rapidly with the degree of backtracking, accuracy is not substantially improved above $N_{\rm checks}= 2^3$, at which limit the CPU time is reduced to $~8$ seconds. This tradeoff further suggests that the performance of the algorithm can benefit greatly from data structuring techniques such as approximate nearest-neighbor searches \cite{eggersmann2021efficient}, especially for large systems and material data sets. Furthermore, we recall that the algorithm is embarrassingly parallel, which opens the way for efficient parallel implementations, e.~g., on GPU machines \cite{barash2017gpu}.

\begin{figure}[ht]
    \centering
	\begin{subfigure}{0.45\textwidth}\caption{}
    \centering
    \includegraphics[height=4.5cm]{Figures/cpu.pdf}
	\end{subfigure}%

	\begin{subfigure}{0.45\textwidth}\caption{}
    \centering
    \includegraphics[height=4.5cm]{Figures/error.pdf}
	\end{subfigure}
    \caption {Lightweight space structure \cite{lee2016architecture}. Influence of the number of checks for backtracking approximate in the nearest-neighbor search. a) CPU times in seconds on a 12-Core AMD Ryzen 9 3900X machine. b) Error $e_{KS}$ of node eccentricity $u_e$ with error bars derived from a sample of $25$ calculations.}
    \label{fig:pactruss_cpu_study}
\end{figure}

\section{Summary and concluding remarks}

The behavior of many materials is characterized by a {\sl likelihood measure} instead of a constitutive relation. Specifically, the material likelihood measure expresses the likelihood of observing a material state, understood as a point in stress-strain space, or {\sl phase space}, in the laboratory. We have assumed that the material likelihood measure is known only partially through an empirical point-data set in phase space. The state of the solid or structure is additionally subject to compatibility and equilibrium constraints that encode, in particular, all information pertaining to geometry, boundary conditions and loading. The problem is then to infer the {\sl posterior likelihood} of a given outcome of interest.

We have presented an {\sl ansatz}-free Data-Driven method of inference that determines likelihoods of outcomes directly from the empirical material data. In particular, the method of inference requires no material or prior modeling. The computation of expectations is reduced to explicit sums over points in local material data sets and quadratures over admissible states, i.~e., states satisfying compatibility and equilibrium. The method reduces to the max-ent Data-Driven approach proposed in
\cite{kirchdoerfer2017data} when the analysis is restricted to the computation of maximum-likelihood outcomes. In its present form, the Data-Driven method of inference delivers a full characterization of expectations and distributions of outcomes from general material data sets.

The complexity of the Data-Driven inference calculations is linear in the number of material data points and in the number of members in the structure or Gauss points in the solid. Efficient population annealing procedures and fast-search algorithms for accelerating calculations have been presented. Specifically, the population annealing algorithm reduces the complexity of the calculations and it renders it linear in the size of the material data-point population. The calculations are further accelerated by recourse to a hierarchical $k$-means algorithm for the computation of the annealing energy with limits on the number backtracking operations.

The convergence properties of the method with respect to the empirical data have been assessed with the aid of a standard verification test based on Gaussian material data. Robust convergence of posterior distributions is obtained, consistent with quantitative error estimates derived from analysis \cite{Conti:2021, Conti:2022}. It bears emphasis that, in these tests and in all subsequent calculations, the empirical material-data sets are the only input to the calculations and, once generated, the underlying distribution whence they are sampled is discarded altogether and plays no subsequent role.

The full power of the method becomes apparent in applications to strongly non-Gaussian material data sets. Thus, we have presented an example concerned with brittle materials exhibiting random tensile strength. Specifically, the tensile strength of the material is assumed to obey Weibull statistics. The computed posterior distributions in this case are strongly multimodal and encode complex behavior characterized by multiple alternatives. It bears emphasis that the material data-sets arising in this example challenge conventional methods in that: i) there is no reduced manifold in stress-strain space, and ii) the prior distribution is strongly non-Gaussian. Indeed, owing to the random character of the tensile strength of the solid, the feasible material region in stress-strain space is not a manifold, but a full set. In addition, in Bayesian approaches to inference it is common to hypothesize Gaussian priors for the observational noise as a matter of convenience. Evidently, such  hypothesis is not born out in the Weibull strength statistics scenario of the current example and may lead to lack or convergence or, worse still, to convergence to the wrong limit (cf.~\cite{Owhadi:2015}). Against this backdrop, the ability of the proposed Data-Driven method of inference to deal effectively with general material data sets and complex behavior without need for models, hypothesis or assumptions bears remark.

Finally, we have presented a simple benchmark case, concerned with the deformations of a lightweight space structure, that demonstrates the expected computational performance of the method. The calculations attest to the effectiveness of acceleration methods such as population annealing and hierarchical $k$-means searches. While the benchmark calculations presented here are of modest size, the linear scaling of the algorithm with system size and population size, as well as the embarrassingly parallel property of the algorithm, bode well for applications at scale.

\section{Acknowledgements}
All authors gratefully acknowledge the financial support of the Deutsche Forschungsgemeinschaft (DFG) and French Agence Nationale de la Recherche (ANR) through the project "Direct Data-Driven Computational Mechanics for Anelastic Material Behaviours" (ANR-19-CE46-0012-01, RE 1057/47-1, project number 431386925) within the French-German Collaboration for Joint Projects in Natural, Life and Engineering (NLE) Sciences. SR gratefully acknowledges the funding of the DFG projects "Model order reduction in space and parameter dimension - towards damage-based modeling of polymorphic uncertainty in the context of robustness and reliability" (SPP 1886, project number 312911604) and "Coupling of intrusive and non-intrusive locally decomposed model order reduction techniques for rapid simulations of road systems" (TP B05/TRR 339, project number 453596084). MO is greatful for support from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) {\sl via} project 211504053 - SFB 1060; project 441211072 - SPP 2256; and project 390685813 -  GZ 2047/1 - HCM.
\bibliographystyle{unsrt}


\section{Introduction}

The application of deep reinforcement learning to solve single player Markov Decision Processes are relatively well-researched in today's machine learning literature, however, the application of novel deep reinforcement learning methods to solve multi-agent competitive MDP's are relatively few. Particularly, the challenge revolves around competing objectives and solving for, often computationally intractable, Nash equilibria \citep{Nash:1950} in a competitive setting, where agents cannot merely maximize only their respective Q functions. Nash equilibria in Markov Games are particularly useful in the area of \textit{dynamic pricing}, a well-known problem in today's online eCommerce marketplaces.  In a \textit{dynamic pricing} game, multiple firms compete for market share of their product on the same online platform. To model this market, we adopt a modified version of Bertrand oligopoly \citep{bertrand:1883}, where the marginal cost of production is 0, and where purchasing behaviour is probabilistic.

\subsection{Market Oligopoly}

Especially in today's online commerce industries, dynamic pricing provides platform firms with an adversarial advantage and increases profit margins in a competitive economy. In an ideal case, a firm's pricing strategy should consider the actions of other firms, providing the justification for Nash equilibrium policies over unilateral reward maximization policies. In previous approaches, \citep{Liu:2019-dyp} modelled the dynamic pricing problem using deep reinforcement learning, however, focused primarily on single agent revenue maximization and not solving for multi-agent Nash equilibria. \citep{Schlosser:2018} apply stochastic dynamic programming in a simulated environment to estimate consumer demand and maximize expected profit, however, do not consider the equilibrium solution of the market. \citep{vandergeer:2018-dpc} model market demand using parametric functions, with unknown parameters, yet they do not consider a Nash equilibrium nor the actions of previous pricing on demand. Thus, in almost all recent research on the dynamic pricing problem under uncertainty, there exists a desideratum to compute and promote Nash equilibria for market pricing strategies. Such strategies can be applied to eCommerce platforms such as Amazon auto pricing, allowing firms to opt for an automated pricing algorithm.

A large body of literature on dynamic pricing has assumed an unknown demand function and tried to learn the intrinsic  relationship between price and demand. \citep{harrison2012bayesian} showed that a myopic Bayesian policy may lead to incomplete learning and poor profit performance. \citep{liu2020dynamic} studied a joint pricing and inventory problem while learning the price-sensitive demand using a Bayesian dynamic program. Many recent studies have revolved around non-parametric methods. \citep{besbes2009dynamic} developed a pricing strategy that balances the trade-off between exploration and exploitation. Furthermore, in \citep{Liu:2019-dyp} the demand function was a black box system where the demand was approximated from experience replay and live online testing in a reinforcement learning framework.

Dynamic pricing in oligopolies presents multiple problems. One aspect is the optimization problem, where we learn a Nash equilibrium policy when presented with the parameters of the market demand function, and second is the learning of the demand parameters itself. For example, \citep{den2014simultaneously} and \citep{keskin2014dynamic} studied parametric approaches, using maximum likelihood to estimate unknown parameters of the demand function. Recent research has been conducted in demand learning, for example in cases where multiple products exist, and firms can learn to price on an assortment of items \citep{ferreira:2022demand} using multinomial logit regression, or when multiple products exist in a demand learning scenario, and a limited inventory exists, Multi-Armed Bandit approaches have shown to be effective at devising profitable policies \citep{ferreira:2018thompson}. However, these approaches do not consider past pricing actions influencing the future market demand parameters. Furthermore, they do not consider competing firms aiming for a Nash equilibrium.

We propose a demand market where past pricing affects the future market demand in a market competition. To search for a Nash equilibrium, we apply a multi-agent Markov Decision Process, where the state transition is driven by the joint pricing actions of all players (or firms). To simplify the problem, we assume that each firm has unlimited inventory, and there is only a single product with no substitutes.

\subsection{Multi-agent Markov Decision Processes}

Multi-agent Markov Decision Processes, or \textit{Markov Games} (MG's), constitute an important area of research, especially when multiple agents are self-interested, and an exact or approximate Nash equilibrium for the system is desired. The computation of Nash equilibria additionaly presents great difficulty when searching over the enormous joint state action space of the problem, and approximations to this search problem do exist \citep{Porter:2008}. Moreover, the existence of multiple Nash equilibria can further complicate the solution, as some solutions may be more desirable than others.

Yet approximate search algorithms require prior knowledge of the joint reward function and are often limited to two players, modelled by a \textit{best response function}. We treat our payoff function as an \textit{oracle} \citep{Vorobeychik:2008}, meaning we have knowledge regarding the parametric form of the payoff function, however, there is no knowledge to the agents in the MG regarding how the reward function generates rewards for the players. The state of the game is visible to the agents, yet the agents have no knowledge regarding the state transition function. Nevertheless, the market demand parameters and Nash equilibria are known for purposes of experimentation. This eliminates the need for pre-period model fitting to pre-estimate market demand parameters \citep{ferreira2016:6analytics} allowing us to compare our solution to the theoretical solution, as opposed to an empirical estimate.

Model-based approaches furthermore exist to provably find the existence of a Nash equilibrium \citep{Zhang:2020}. However, we aim for a Nash equilibrium solver which is model-free, where no knowledge regarding the reward or transition function is known to the agents. In recent literature, \citep{Raponi:2021} present a model-free algorithm to solve competitive MDP's when the environment parameters can be altered by an external party - but this is not always the case in market economics. Furthermore, \citep{Sayin:2021} propose a model-free decentralized deep learning framework, where the agent is blind to the other agents actions, and convergence to a Nash equilibrium is proven. \citep{Guan:2021} propose a convergent solution for zero-sum MG's via entropy-regularization. However, both \citep{Guan:2021} and \citep{Sayin:2021} impose a number of theoretical restrictions on the MG in order for this convergence to occur. Moreover, \citep{Kozuno:2021} present a provably convergent MG solver for imperfect information restricted to two agents. Nevertheless, we are concerned with an approximate solution to a full information MG for N agents.

\textbf{Contributions:} This work outlines a methodology for Deep Q Learning, as introduced in \citep{mnih:2015}, by extending the framework to multi-agent reinforcement learning (MARL) with a Nash equilibrium objective based on the methods in \citep{hu:2003} and \citep{wang:2003-mgames}. Although the framework presented in \citep{hu:2003} is theoretically sound, the solution to the Nash equilibrium function is often intractable. We therefore apply approximation methods to compute a Nash equilibrium. Black box optimization is applied to estimate an $\epsilon$-Nash Equilibrium policy, and this approximation is learned as a series of neural networks. This MARL model is then applied to the domain of dynamic pricing.

\section{Dynamic Pricing Game} \label{sec:dynamic-pricing-game}

An oligopoly across $N$ firms is modelled as a multi-agent Markov Decision Process with a fully observable state action space. The state space is represented by the fully-observable demand-influencing reference price $\bar{p_t}$ being the state of the game $s^t$. In discrete intervals at time $t$, each agent issues a price for the item to be sold to the general market, $x^t_n$. The reference price, of the market at time $t$ is determined by a stochastic function of all the agents actions at time $t-1$. To focus the problem strictly on dynamic pricing, for each firm we assume a sufficiently large inventory, with no marginal costs of production, no possibility for inventory stockouts, no holding costs, and the capacity to always meet any realization of market demand.

\subsection{Markov Game Parameters} \label{sec:markov-game-params}

The joint action space $\{x_1^t \ ... \ x_n^t\} \in \mathbf{X}$ constitutes the current actions of all agents at time $t$, which drive the state transition by setting prices at time $t$ and represent the set price of the item from agent $n$ at time $t$. The state-action-reward space is defined as a tuple $(s^{t}, x_1^{t}, . . . , x_N^{t}, r_1^{t}, . . . , r_N^{t}, s^{t+1})$ where $r_n^{t}$ is the reward for agent $n$ at time $t$. The joint reward can be written as $\bar{r}_t = (r_1^{t}, . . . , r_N^{t})$, and the joint action can be written as $\bar{x}_t = (x_1^{t}, . . . , x_N^{t})$. $N$ denotes the number of agents. The exact transition probabilities are not known to the agents. Each agent must learn the demand function and optimal strategy as the MG progresses. A state $s^{t}$ is determined by the reference price of the market observable to all agents, $s \in \mathbb{R}$. We discretize the action space into even segments representing a firm's price. Where each segment represents an action $x_n \in X_n$, and $X_n$ is the action space for any agent $n \in N$.

\subsection{Reference Pricing} \label{sec:refernce-pricing}

The first pillar of any market dynamic is constituted by the \textit{demand function}, dictated by both the historical and contemporary prices, of the firms in the market. Demand processes can be ideally linear and continuous, however may not be guaranteed to be stationary or continuous \citep{vdenboer:2020-discon}, and can be subject to various exogenous factors such as environmental carryover effects \citep{rao1993pricing}.

We create an idealization of a market based on \citep{taudes2012integrating}, where a dynamic pricing problem, with reference price effects and a linear demand function, is imposed. The expected current demand of a product is a function of the average price set by all firms, denoted as $\tilde{x}$, and a historical reference price, $\bar{p}$. The reference price $\bar{p}$ is given, and cannot be modified during the current time $t$, moreover, it is a function of the immediate past, $t-1$. Although \citep{taudes2012integrating} is not necessarily the only model that incorporates reference pricing \citep{mazumdar2005reference}, \citep{fibich2003explicit}, \citep{popescu2007dynamic}, \citep{heidhues2014regular}, we adopt it due to its simplicity and the existence of provable Nash equilibria in an oligopolistic framework.

\begin{align}
    f(\tilde{x}) = \beta_0 + \beta_1 \tilde{x} + \beta_2(\tilde{x} - \bar{p}), \quad \beta_0  &\geq 0, \ \beta_1 < 0, \ \beta_2 \leq 0, \ f(\tilde{x}) \geq  0, \tilde{x} \geq 0  \label{eq:taudes-demand}

\end{align}

$\epsilon_D$ is defined as the noise from a Poisson process. Such a Poisson process has the arrival rate $\lambda_t(\tilde{x}) = f(\tilde{x})$ from Eq. \eqref{eq:taudes-demand}, and standard deviation $\sigma = \sqrt{\lambda_t(\tilde{x})}$. Furthermore, \citep{taudes2012integrating} make the stipulation of decreasing demand with respect to increasing price, as illustrated in Condition \eqref{eq:taudes-demand}.

This demand-influencing reference price can be affected by a variety of factors, such as inflation, unemployment, psychological perception \citep{raman:2002}. Moreover, in many proposed oligopoly models, as in \citep{janiszewski:1999range} and \citep{briesch:1997comparative}, the reference price is dictated by a historical joint market price of multiple firms. However, modelling a competitive market oligopoly with an autocorrelated reference price in a MG setting is not heavily investigated until now. In our model, we focus on designing a market whose reference price is driven by the historical $t-1$ joint market price, and additional factors which also affect the reference price are represented as noise $\epsilon_D$. Thus the transition of the reference price is determined by the average of the previous joint prices plus some Guassian randomness, $\epsilon_D$.

In our experiment, the reference price of a product is determined by the previous actions of the firms. We express the reference price function as $\bar{p}(t): \ \mathbb{R}^{N} \xrightarrow[]{} \mathbb{R}$, mapping a vector of historical pricing actions $\mathbf{x}^{t-1}$ to $\mathbb{R}$.  In the beginning of the Markovian process, the reference price is randomly generated within the action space of the pricing range. The stochastic nature of the market price transition entail the Markov property of the MG.

\begin{align}

  \bar{p}(t) &\vcentcolon= \frac{1}{N} \sum^{N}_{n=1} x^{t-1}_n + \epsilon_D, \quad \text{where} \quad \bar{p}(t=0) = \frac{1}{N} \sum^{\mathbf{x^0}}_{n=1} x^{0}_n + \epsilon_D
\end{align}

\subsection{Probabilistic Demand Function} \label{sec:profit-func}

The expected profit for player $n$ is $\mathbf{E}[\Pi_n(x_n)]$, which equals revenue as the marginal cost is of production is assumed to be 0, is defined as,

\begin{align}
    \mathbf{E}[\Pi_n(x_n)] &= \Phi_n(x_n, x_{-n}) x_n f(\tilde{x}) \label{eq:profit-func1} \\
    \Phi_n(x_n, x_{-n}) &= \frac{ e ^{\phi(x_n)} }{\sum^{N}_{i=1} e^{\phi(x_i)}  } \label{eq:softmax-elas}
\end{align}

$\Phi_n(x_n, x_{-n})$ in Eq. \eqref{eq:softmax-elas} represents the probability that a customer purchases the item sold by firm $n$, where the player $n$ prices their merchandise with price $x_n$. $f(x_n)$ is the expected demand during the single stage game for a fixed time period. Following the quantal response function from \citep{luce:1959} and \citep{goeree:2020}, we define a \textit{purchase elasticity function}, $\phi_n(x_n, x_{-n})$,

\begin{align}
  \phi(x) = &
    \begin{cases}
      b_n - a_n x
        &\text{for \quad} 0 < x < b_n / a_n \\
      0 &\text{for \quad}  x > b_n / a_n
    \end{cases} \label{eq:puchase-elas-function} \\
    \text{where} \quad &  x \in \mathbb{R}, x > 0, a \geq 0, b_n = 1 \label{eq:purchase_elas_price_constraints}
\end{align}

In $\phi(x)$, $b_n$ represents the maximum weighted contribution to the probability of an item being purchased by a customer given a price. When the price $x$ is 0, this measure is $b_n$. For simplification, we assume the linear marginal decline of this measure of all players $a_n$ as equal, $a_n = a$, that is the market has the same elasticity regarding the marginal decrease of the customer willingness to purchase from any firm as price increases, with a negative slope $a$. The probability of a customer choosing to purchase from firm $n$ at price $x_n$ among prices set by other firms $x_{-n}$ is defined by combining a softmax function and purchase elasticity function, $\phi(x)$. This mechanism will prevent a solution where each firm undercuts each other to set their prices to 0, as the lowest price need not guarantee that a consumer will buy their product.

\section{Nash Equilibrium Conditions} \label{sec:eq-setting-main}

Computation of exact Nash, or approximate $\epsilon$-Nash, equilibria is an NP-hard problem and notoriously difficult to compute \citep{Daskalakis:2006}. This involves a search over the entire joint state action space, computing the value of each state under a candidate policy. Furthermore, it involves knowledge of the joint Q function and the transition probability of the system. In our scenario, we have the condition that all agents are identical, therefore the solution of one agent can apply to the solution of another.

\subsection{Theoretical $\epsilon$-Nash Equilibria}

Under the market outlined in Section \ref{sec:dynamic-pricing-game}, given the market inputs $\bar{p}$ and $\tilde{x}$, in any pricing strategy, there exists an optimal deviation $d^* \in \mathbf{R}$, such that $\mathbf{E}[\Pi_n(x_n - d^*)]$ will yield the maximum profit advantage $\epsilon$.

\begin{align}
    \epsilon &= \underset{d^* \in \mathbf{R} }{\mathrm{max}} \ \Bigg(  \mathbf{E}[\Pi_n(x_n - d^*)] - \mathbf{E}[\Pi_n(x_n)] \Bigg) \label{eq:nash-def-diff}

\end{align}

An increase in the individual profit of an agent $n$ from unilaterally deviating is denoted as $\epsilon$. Given this optimal deviation, $d^*$, a maximum theoretical upper bound on the profit advantage $\epsilon$ that can be obtained. Eq. \eqref{eq:d-star-main} provides the theoretical value of $d^*$ which is a function of both reference price $\bar{p}$ and market price $\tilde{x}$, as well as the fixed market parameters $\beta_0, \beta_1, \beta_2$. A derivation of $d^*$ can be found in Appendix \ref{sec:eq-setting}.

\begin{align}
    d^* &= \frac{ \sqrt{c_1^2 - c_1 + 4(c_2-1)c_2 - 2c_2} }{2c_2}  \label{eq:d-star-main} \\
    \text{where} \quad  c_1 &= \frac{- (\beta_1 + \beta_2)}{f(\tilde{x})} - \frac{1}{\tilde{x}}, \quad c_2 = \frac{- (\beta_1 + \beta_2)}{f(\tilde{x})\tilde{x}} \nonumber
\end{align}

\begin{figure}[!htb]
\minipage{0.5\textwidth}
  \includegraphics[width=\linewidth]{plots/market2.png}
  \caption*{\textbf{Market Scenario 1:} \\ $\beta_0 = 15, \beta_1 = -1.05, \beta_2 = -3.1, a = 0.1$
\endminipage\hfill
\minipage{0.5\textwidth}
  \includegraphics[width=\linewidth]{plots/market3.png}
  \caption*{\textbf{Market Scenario 2:} \\ $\beta_0 = 25, \beta_1 = -0.6, \beta_2 = -6.1, a = 0.1$
\endminipage\hfill
\caption{Surface plot of potential advantage from deviation ($\epsilon$-deviation) from Nash equilibrium with respect to market price $\tilde{x}$, and reference price $\bar{p}$, with their respective market parameters $\beta_0, \beta_1, \beta_2, a$ for arrival of a single sales event. } \label{fig:market-topo-main}
\end{figure}

Using the Eq. \eqref{eq:d-star-main}, as visualized in Fig. \ref{fig:market-topo-main}, multiple Nash equilibria can exist when $\epsilon = 0$, or when $\epsilon$ is sufficiently small (see Section \ref{sec:results}). Such behaviour occurs when one agent unilaterally deviates from $\tilde{x}$, under reference price $\bar{p}$. Thus, the lower plateau from Fig. \ref{fig:market-topo-main} represents the region where a theoretical Nash Equilibrium exists. Bounds on the state space and action space are limited to $|S| = |X_n| = 10$, that is there are 10 discrete intervals for each state $s^t$ and action $x_n^t$, representing price values.

\subsection{Nash Equilibrium in Markov Games ($\delta$)}

The value of a policy, at each state, can be represented by a joint action which maximizes the joint Q function \citep{Watkins:1992} under policy $\pi(s)$. The probability of  agent $n$ taking action $x_n$ is defined as $\pi_n(s, x_n)$.

\begin{align}
    v(s, \pi_n, \pi_{-n}) = \underset{\bar{x} \in \mathbf{X} }{\mathrm{max}} \ \prod_{n=1}^N \pi_n(s, x_n) Q_n(s, x_{n}, x_{-n}) \label{eq:val-joint-policy}
\end{align}

An $\epsilon$-Nash equilibrium is defined as a joint policy such that the reward of a single stage game will not result in a greater payoff to the agent $n$ by more than $\epsilon$, when any agent unilaterally deviates from said joint policy. Provided $\epsilon$, we consider a bound on the corresponding MG, $\delta$, which defines a bound on the gain of the \textit{value} of a policy, $v(s, {\pi_{n}}^*, \pi^*_{-n})$, should agent $n$ unilaterally deviate with an alternative policy. The solution to Eq. \eqref{eq:val-joint-policy} can be computed by searching over the joint action space $\{ x_n, x_{-n} \}$.

\begin{align}
    \underset{\bar{x} \in \mathbf{X} }{\mathrm{argmax}} \prod_{n=1}^N \pi_n(s, x_n) Q_n(s, x_{n}, x_{-n}^*) &\leq \underset{\bar{x} \in \mathbf{X} }{\mathrm{argmax}} \prod_{n=1}^N \pi_n(s, x_n) Q_n(s, x_{n}^*, x_{-n}^*) + \delta \label{eq:nash-cond-q} \\
    v(s, \pi_n, \pi_{-n}^*) &\leq v(s, \pi_{n}^*, \pi_{-n}^*)  + \delta, \quad \forall s \in \mathbf{S} \label{eq:nash-cond-v}
\end{align}

$\delta$ serves as an upper bound on $\epsilon$, therefore the minimization of $\delta$ will also minimize any possible existence of $\epsilon$ for each single stage game in the MG. In fact, the existence of $\epsilon$ implies the existence of upper bound $\delta$, we provide a proof of this existence in Appendix \ref{sec:eps-delta-exis}.

\section{Multi-Agent Nash Q Learning} \label{sec:ma-ql}

In our model, the game provides full information, where information regarding the state and the actions of all agents are visible to any agent. This allows for \textit{experience replay} for Q learning \citep{mnih:2015} to occur, and the Q function can be shared for all agents as they are assumed to be identical. The joint action space $\bar{x} = (x_1, ..., x_N)$ is defined as the combined actions of all agents at time $t$. Normally, in Q learning, the update mechanism searches for the action that maximizes the Q function, however, in \textit{Nash Q learning}, we must update the Q function with the solution to the Nash equilibrium. In our model, we represent the tabular Q function via a function approximator, which is denoted simply as $Q(s, \bar{x})$. As new experiences are obtained when the MG is played, the Nash Q value is updated. We utilize the Q update mechanism defined in \citep{hu:2003} as the update mechanism for the Nash Q estimator for Q learning. Given a Q function and a Nash policy, we search over the joint action that maximizes the scaled Q function \citep{laumonier:2005-multiagent}, as in Eq. \eqref{eq:nash-q-max}. In our representation, the Q function is a vector, returning the respective Q value for each agent, based on the joint probability input for the joint policy.

\begin{align}
    Q'(s, \bar{x}) &\xleftarrow[]{} (1-\alpha) Q(s, \bar{x}) + \alpha [r + \gamma \mathcal{N}(s') ] \label{eq:q-update} \\
    \bar{x}^* &= \underset{\bar{x}}{\mathrm{argmax}} \ Q(s', \bar{x}) \prod_{i=1}^N \pi_n^*(s', x_n) \label{eq:nash-q-max}\\
    \mathcal{N}(s') &= Q(s', \bar{x}^*) \prod_{i=1}^N \pi_n^*(s', x^*_n)  \label{eq:nash-operator}

\end{align}

$\mathcal{N}(s')$ represents a \textit{Nash Operator}, indicating the Q value of a Nash equilibrium, this is approximated by a scaling factor, according to the Nash equilibrium policy, $\pi^*$, multiplied by the Q value of the joint action, $Q(s, \bar{x})$. Extending from the work of \citep{hu:2003}, we introduce a neural network for mapping any state action combination to its Nash policy, where the Nash scaling factor derives from and is used to compute $\mathcal{N}(s')$ in conjunction with the Q function. The Q-function used in Eq. \eqref{eq:nash-operator} and \eqref{eq:q-update} is represented by a deep neural network, simply referred to as the Nash Q Net $Q(s)$. In a full information game with homogenous agents, a shared neural network representing $Q(s, \bar{x})$ can be used for all agents in the Markov Game (if not homogeneous a separate Q network must be retained and updated separately for each agent). The Q Network parameters representing $\mathcal{N}(s')$ is learned in the same method as Deep Q learning \citep{mnih:2015}, the key innovation, is that the scaling factor to compute $\mathcal{N}(s')$ is obtained via a Nash Policy Net, $\Psi(s)$ (defined later in Section \ref{sec:nash-pol-learn}).

\subsection{Estimating Value Advantage via Black Box Optimization}

We apply a deep neural network to represent a mapping of a joint policy to its respective $\delta$ from Eq. \eqref{eq:nash-cond-v}, designated as $\Gamma(\pi) \xrightarrow[]{} \bm{\delta}$, where $\bm{\delta}$ is a vector containing the $\delta_s$ of each state. As the gradient of $\Gamma(\pi)$ cannot be easily evaluated, we apply gradient free black box optimization for $\delta$ minimization. Trust Region Optimization (TRO) has been shown to be effective for solving high-dimensional optimization problems \citep{wang:2021-lamcts} \citep{ericksson:2020-turbo} \citep{Diouane:2021-trego} \citep{Regis:2016}, particularly when the computational resources are available. To compute the existence of an $\epsilon$-Nash equilibrium in the high-dimensional policy space efficiently, we apply model-based Bayesian Optimization via Trust Region optimization (\texttt{TuRBO}) from \citep{ericksson:2020-turbo}. \texttt{TuRBO} combines standard TRO with a multi-armed bandit system via Thompson Sampling \citep{Thompson:1933} to optimize for local multiple trust regions simultaneously, with sufficient convergence time. However, the candidate generation step in \texttt{TuRBO} is not constrained to account for the valid joint probabilities of each agent, in which the sum of probabilities for each agent in its respective action space must sum to 1. We alter \texttt{TuRBO} by simply normalizing the original candidate probabilities over each set of probabilities belonging to an agent $n$. The resulting modified algorithm is denoted \texttt{TuRBO-p}, where any candidate value generated by \texttt{TuRBO-p} will have joint probabilities that sum to 1 for policies corresponding to each agent.

\begin{align}
    V(s, \pi) &= \underset{\bar{x}}{\mathrm{max}} \ Q(s, \bar{x}) \prod_{i=1}^N \pi_n(s, x_n) \\
    \delta &= \underset{\pi_n'}{\mathrm{max}} \ \Bigg( V(s, \pi_n', \pi_{-n}) - V(s, \pi_{n}, \pi_{-n}) \Bigg) \quad \forall s \in \mathbf{S} \label{eq:delta-max}

\end{align}

The maximization problem is formulated in Eq. \eqref{eq:delta-max} representing the maximum gain in value of a policy, where agent $n$ deviates from policy $\pi_{n}$ with $\pi_n^{'}$.

\subsection{Nash Policy Learning} \label{sec:nash-pol-learn}

The $\epsilon$-Nash policy is defined as a joint policy $(\pi^*_n, \pi^*_{-n})$ that minimizes $\delta$ from Eq. \eqref{eq:nash-cond-v}. Drawing inspiration from \citep{Ceppi:2010}, where a Nash equilibrium is found by effectively minimizing for $\epsilon$ via linear programming, we apply a similar technique of $\delta$-minimization. However, in \citep{Ceppi:2010}, the model parameters are known to the agents, and the Markov game was constrained to two players. In our MG, the joint reward function must be learned. Therefore, we perform approximate search over the policy space instead of using exact methods to discover any possible approximate Nash equilibrium policies. In principle, each state has maps to a corresponding probability in accordance with the Nash Policy which minimizes $\delta$ in Eq. \ref{eq:delta-max}. However, a table keeping track of such an enormous state policy space is not feasible. Therefore, we implement a Nash Policy Net $\Psi(s) \xrightarrow[]{} \pi^*$ to learn the state to policy mapping, which is the joint policy $\pi$ producing an approximate Nash equilibrium as approximated via \texttt{TuRBO-p}, designated as $\hat{\pi}^*(s)$.

\begin{align}
    \hat{\pi}^*(s) = \underset{\pi}{\mathrm{argmin}} \ \Gamma(\pi)_s \label{eq:gamma-min}
\end{align}

\begin{algorithm}[h!]
\caption{Nash equilibrium learning}\label{algo:nash-eq-learning}
\begin{algorithmic}[1]

    \State Initialize state joint policy $\pi$.
    \State Initialize random parameters for $\Psi(s_t)$, $Q(s_t)$ and $\Gamma(\pi)$, as $\theta_{Q}$, $\theta_{\Psi}$, and $\theta_{\Gamma}$.
    \State Initialize MDP Environment, $\mathcal{M}$.

    \For {$e \in episodes $}:
        \State Get initial state $s_0$ from $\mathcal{M}$.
        \For {$t \in (0, t_{max})$}: Iterate until end of episode.

            \State Get action probabilities $\pi(s_t, x_n)$ from Nash Policy Net $\Psi(s_t)$.
                \For {$n \in (1, N)$}: Iterate through agents.

                    \State Sample selected action from Nash policy. $x_n \thicksim Multinomial(s_t, \pi(\overline{x}))$
                \EndFor
            \State Obtain joint action $\overline{x} \xleftarrow[]{}  \{ x_1, ... , x_N \}$
            \State Assign  $\delta_s \xleftarrow[]{} \underset{\pi_n'}{\mathrm{max}} \ \Bigg( V(s, \pi_n', \pi_{-n}) - V(s, \pi_{n}, \pi_{-n}) \Bigg)$ via \texttt{TuRBO-p}.
            \State Find the Nash policy $\pi^* \xleftarrow[]{}  \underset{\pi}{\mathrm{argmin}} \ \Gamma(\pi, s_t)$ via \texttt{TuRBO-p}.
            \State Execute joint action $\overline{x}$ in $(\mathcal{M}, \ s_{t})$ to obtain $s_{t+1}, \overline{r}_t$.
            \State Append experience $(s_t, \overline{x}_t, \overline{r}_t, \epsilon_s, \pi, s_{t+1} )$ to experience replay $\mathcal{D}$.

            \State Update State $s_{t} \xleftarrow[]{} s_{t+1}$.

            \If {$|\mathcal{D}| > batchsize$}:
                \State Sample minibatch $d$ from $\mathcal{D}$.
                \State Set $Q'(s_t, \bar{x}) \xleftarrow[]{}
                    \begin{cases}
                      r_t &\text{for \ } s_t \ \text{terminal} \\
                      (1-\alpha) Q(s, \bar{x}) + \alpha [r + \gamma \mathcal{N}(s') ] &\text{for \ }  s_t \ \text{not terminal}
                    \end{cases} $

                    \State $\mathcal{L}_{Q} \xleftarrow[]{} (Q'(S_t, \bar{x}) - Q(S_t, \bar{x}))^2$ Nash Q Update.
                \State $\mathcal{L}_{\Psi} \xleftarrow[]{} (\hat{\pi} - \pi^*)^2$ Nash policy update, where $\hat{\pi}$ is from $\Psi(s)$.
                \State $\mathcal{L}_{\Gamma} \xleftarrow[]{} (\hat{\epsilon} - \epsilon)^2$ Nash $\epsilon$ update, where $\hat{\epsilon}$ is from $\Gamma(\pi)$
                \State Backpropagate loss functions $\mathcal{L}_{Q}$, $\mathcal{L}_{\Psi}$ and $\mathcal{L}_{\Gamma}$ on $\mathcal{Q}$, $\Psi$, and $\Gamma$.
                \State Update model parameters $\theta_{Q}$, $\theta_{\Psi}$, and $\theta_{\Gamma}$ via gradient descent.

            \EndIf

            \State Update joint policy $\overline{\pi} \xleftarrow[]{} \overline{\pi}^*$.
        \EndFor
    \EndFor

\end{algorithmic}
\end{algorithm}

\section{Results} \label{sec:results}

\begin{figure}[!htb]
\minipage{0.5\textwidth}
  \includegraphics[width=\linewidth]{plots/sc1_loss.png}
  \caption*{Loss behaviour of Scenario 1 with $n=3$.
\endminipage\hfill
\minipage{0.5\textwidth}
  \includegraphics[width=\linewidth]{plots/market2_loss_5n.png}
  \caption*{Loss behaviour of Scenario 2 with $n=5$.
\endminipage\hfill
\caption{Decreasing loss behaviour of batch update during Deep Q Learning $Q$ (left y axis) and Nash Net update $\Psi(s)$ (right y axis) demonstrating that $\Psi(s)$ has learned a representation of the Nash policy (x-axis represents iterations multiplied by episode).} \label{fig:loss_plot_main}
\end{figure}

\textbf{Stabilization of realized market rewards to Nash equilibrium:} Fig. \ref{fig:reward_plots_main} represents the convergence of the average market reward of a single agent (randomly selected as agent 0) towards a Nash equilibrium. We superimpose the boundaries of the theoretical Nash equilibrium reward over the plot. The reward is obtained from the boundaries of theoretical Nash equilibria, where $\epsilon < \epsilon_l$ ($\epsilon_l = 0.0001$). The topology of this function is illustrated in Fig. \ref{fig:market-topo-main}, and for each state, or reference price, there exists a boundary where the policy deviation of any agent can occur without significant unilateral reward gain. The expected reward is computed using the market model parameters per Scenario ($\beta_0, \beta_1, \beta_2, a$). From the expected reward, we compare it with the positive difference of theoretical equilibrium reward $\Pi' = f(x_n') \ x_n'$ where $x_n'$. This exists as two boundary points indicating the area where $\epsilon < \epsilon_l$. For each episode the average reward per timestep over the episode lengths is recorded. From Fig. \ref{fig:reward_plots_main}, we observe the average reward per agent of the system per episode (blue dashed line), and the average reward of a single agent per episode (orange dashed line), to converge within the boundary of Nash equilibrium (blue shade). However, due to randomization from black box optimization and/or stochastic behaviour of the Markov Game this expected reward can sometimes fall outside of the NE region.

\textbf{Diminishing loss function:} We demonstrate empirically, as shown in Fig. \ref{fig:loss_plot_main}, that the loss function of $Q(s)$ and $\Psi(s)$ decrease with each RL episode. This merely indicates that a representation of the function mapping state to its corresponding Nash policy based on policy-to-$\epsilon$ model $\Gamma(\pi)$ is being learned, and a stable solution is proposed. Hyperparameters are presented in Appendix \ref{sec:hardware-hyperparam}.

\begin{figure}[!htb]
\minipage{0.5\textwidth}
  \includegraphics[width=\linewidth]{plots/sc1_rewards.png}
  \caption*{Market rewards per episode Scenario 1, $n=3$.}\label{fig:market3_reward}
\endminipage\hfill
\minipage{0.5\textwidth}
  \includegraphics[width=\linewidth]{plots/market2_rewards_5n.png}
  \caption*{Market rewards per episode Scenario 2, $n=5$.}\label{fig:market2_reward}
\endminipage\hfill
\caption{In a Nash equilibrium, both the market average reward (blue), and single agent reward (orange) should fall within the Nash equilibria boundary (blue shade), as the MG progresses. }\label{fig:reward_plots_main}
\end{figure}

\textbf{Comparison with baseline cooperative Q-learning model:} We compare our NE learning model wit a naive cooperative Q learning model to serve as a baseline measure of model performance. In the naive Q learning model, the objective is simply to maximize each agent's respective rewards, and avoid any kind of NE computation. The Q update function in Eq. \eqref{eq:q-update} is simply a maximization of the joint Q function, instead of using $\mathcal{N}(s')$. From Fig.\ref{fig:baseline_comp}, we generally observe lower average reward values from the baseline model compared to the NE learning model. The NE reward ranges remain very similar. From an economic perspective, the Nash Q learning agents take into account that other agents will compete against them, and the solver optimizes accordingly, whereas, in cooperative Q learning no consideration for the other agents' strategies are taken, and agents may to over-compete, driving average rewards lower in an oligopoly.

\begin{figure}[!htb]
\minipage{0.5\textwidth}
  \includegraphics[width=\linewidth]{plots/sc1_baseline_comp_reward.png}
  \caption*{Scenario 1 - NE Learning vs Q learning, $n=3$.}\label{fig:market3_max_gain}
\endminipage\hfill
\minipage{0.5\textwidth}
  \includegraphics[width=\linewidth]{plots/sc2_baseline_comp_reward.png}
  \caption*{Scenario 2 - NE Learning vs Q learning, $n=5$.}\label{fig:market2_reward}
\endminipage\hfill
\caption{Average revenue per episode from a NE learning model (blue) versus a baseline Q learning model (green), with NE reward boundaries overlayed in their respective colours.  }\label{fig:baseline_comp}
\end{figure}

\section{Conclusion} \label{sec:conclusion}

We created a Markov Game that represents an n-firm oligopoly, based on previously established market models where theoretical $\epsilon$ bounds for a Nash equilibrium policy exist. A black box algorithm is applied to estimate an upper bound on the $\epsilon$ value of a joint policy, represented by neural network $\Gamma(\pi)$. Similarly a Nash Policy Net $\Psi(s)$ is learned to represent the $\epsilon$-minimizing policy from black box optimization, constituting an $\epsilon$-Nash policy. Thus $\Psi(s)$ can be used together with traditional Deep Q learning, to solve for a Nash equilibrium. Empirically, we show that the average reward of all agents, and the reward of a respective single agent, converges to an approximate Nash equilibrium. The limitations of this research are the limited action space of the agents, as well as the identical nature of the agents. Larger scale experimentation under this framework is suggested. Furthermore, the proposed market model could be enhanced by implementing more complex non-linear market oligopoly environments.