# Variational Approach for Efficient KL Divergence Estimation in Dirichlet
Mixture Models

Samyajoy Pal  Department of Statistics, LMU Munich  Christian Heumann
Department of Statistics, LMU Munich

###### Abstract

This study tackles the efficient estimation of Kullback-Leibler (KL)
Divergence in Dirichlet Mixture Models (DMM), crucial for clustering
compositional data. Despite the significance of DMMs, obtaining an
analytically tractable solution for KL Divergence has proven elusive. Past
approaches relied on computationally demanding Monte Carlo methods, motivating
our introduction of a novel variational approach. Our method offers a closed-
form solution, significantly enhancing computational efficiency for swift
model comparisons and robust estimation evaluations. Validation using real and
simulated data showcases its superior efficiency and accuracy over traditional
Monte Carlo-based methods, opening new avenues for rapid exploration of
diverse DMM models and advancing statistical analyses of compositional data.

##  1 Introduction

The Kullback-Leibler (KL) divergence (Csiszár, 1975) stands as a fundamental
measure in statistics, quantifying the statistical distance between
probability distributions. Defined as the expectation of the logarithmic
difference between two probability distributions, KL Divergence is pivotal in
diverse fields. In the context of Mixture Models (McLachlan and Peel, 2000), a
powerful statistical tool with applications ranging from pattern recognition
(Youbi et al., 2016) to bioinformatics (Liu and Yao, 2022), KL divergence
assumes heightened importance. The Kullback-Leibler (KL) divergence (also
recognized as relative entropy) between two probability density functions
$f(x)$ and $g(x)$ is defined by the integral expression:

$D(f\|g)\ \displaystyle\mathop{=}^{{\rm def}}\ \int
f(x)\log\left(\frac{f(x)}{g(x)}\right){\rm dx}\leavevmode\nobreak\ .$

It operates as a measure of the dissimilarity between the probability
distributions encoded by $f(x)$ and $g(x)$. The KL divergence adheres to
fundamental properties, denoted as the divergence properties:

  * •

Self Similarity: $D(f\|f)=0$

  * •

Self Identification: $D(f\|g)=0$ if and only if $f=g$

  * •

Positivity: $D(f\|g)\geq 0$ for all $f,g$.

These properties underscore the significance of KL divergence in capturing the
nuances of distributional disparities, making it a cornerstone in statistical
analyses.

Particularly in the realm of Dirichlet Mixture Models (DMM), the utility of KL
Divergence becomes pronounced. DMM has demonstrated superior performance over
other popular methods such as Gaussian Mixture Models (GMM) for compositional
data (Pal and Heumann, 2022). Compositional data (Aitchison, 1982) is
characterized by elements expressing proportions or percentages within a
whole. This data form is ubiquitous in fields such as genomics and ecology.

While closed-form solutions for KL divergence exist for the Dirichlet
distribution, extending this analytical tractability to Dirichlet Mixture
Models has posed a significant challenge. Past research (Ma et al., 2014)
predominantly turned to Monte Carlo methods to approximate the KL divergence
in DMMs. However, these methods are computationally intensive and time
consuming, which presents a substantial hurdle.

This study addresses these challenges by proposing a variational approach to
approximate KL divergence in Dirichlet Mixture Models. Unlike previous
methods, our approach provides a closed-form solution, substantially enhancing
computational efficiency. This advancement not only facilitates rapid model
comparisons but also ensures a robust evaluation of estimation quality.

The validation of our proposed method using both real and simulated data
demonstrates its superiority in efficiency and accuracy over traditional Monte
Carlo-based approaches. The results underscore the potential of our
variational approximation to accelerate the estimation process while improving
the quality of the estimates. This innovation marks a significant stride in
expediting the exploration of diverse DMM models, thereby fostering
advancements in statistical analyses of compositional data.

##  2 Methods

The Kullback-Leibler (KL) divergence serves as a crucial metric for assessing
the quality of parameter estimation. It quantifies the statistical distance
between the distribution with the true parameter value and the one with the
estimated parameter value, providing a valuable measure of goodness of fit. A
lower KL divergence indicates a better fit, emphasizing the significance of
obtaining accurate parameter estimates. In the subsequent sections, we provide
a concise overview of the estimation process for Dirichlet Mixture Models
(DMM) and introduce various approximations of KL divergence, essential for
quantifying the divergence between distributions.

###  2.1 Estimation of Dirichlet Mixture Model Parameters

Let $\boldsymbol{X}_{1},\boldsymbol{X}_{2},\ldots,\boldsymbol{X}_{N}$ denote a
random sample of size $N$, where $\boldsymbol{X}_{i}$ is a $p$ dimensional
random vector with probability density function $f(\boldsymbol{x}_{i})$ on
$\mathbb{R}^{p}$. An observed random sample is denoted by
$\boldsymbol{x}=(\boldsymbol{x}_{1}^{T},\ldots,\boldsymbol{x}_{N}^{T})^{T}$,
where $\boldsymbol{x}_{i}$ is the observed value of the random vector
$\boldsymbol{X}_{i}$.

The Dirichlet density is given by

$f(\boldsymbol{x}_{i})=\frac{\Gamma(\sum_{m=1}^{p}\alpha_{m})}{\prod_{m=1}^{p}\Gamma(\alpha_{m})}\prod_{m=1}^{p}x_{im}^{\alpha_{m}-1}\leavevmode\nobreak\
,$ (1)

where $\sum_{m=1}^{p}x_{im}=1,x_{im}\text{'s}>0\leavevmode\nobreak\
,\alpha_{jm}\text{'s}>0$ and $\Gamma(\cdot)$ denotes a gamma function. Thus,
the parameter $\boldsymbol{\alpha}=(\alpha_{1},\ldots,\alpha_{p})$ is a
$p$-dimensional vector.
The density of a mixture model with $k$ mixture components for one observation
$\boldsymbol{x}_{i}$ is given by the mixture density

$p(\boldsymbol{x}_{i})=\sum_{j=1}^{k}\pi_{j}f(\boldsymbol{x}_{i}\mid\boldsymbol{\alpha}_{j})\leavevmode\nobreak\
,$ (2)

where $\boldsymbol{\pi}=(\pi_{1},\ldots,\pi_{k})$ contains the corresponding
mixture proportions with $\sum_{j=1}^{k}\pi_{j}=1$, $0<\pi_{j}<1$. The density
component of mixture $j$ is given by
$f(\boldsymbol{x}_{i}\mid\boldsymbol{\alpha}_{j})$ and
$\boldsymbol{\alpha}_{j}$, $j=1,2,...,k$ is the vector of component specific
parameters for each density. Then
$\boldsymbol{\alpha}=(\boldsymbol{\alpha}_{1},\ldots,\boldsymbol{\alpha}_{k})$
denotes the vector of all parameters of the model.
The parameters of Dirichlet mixture model, can be estimated using an EM
algorithm (Dempster et al., 1977). In our previous work (Pal and Heumann,
2022) we presented a Hard version (Celeux and Govaert, 1992) of the EM
algorithm to estimate the parameters. The log likelihood of the model for a
sample of size $N$ is then given by

$\log
p(\boldsymbol{x_{1}},\ldots,\boldsymbol{x_{N}}|\boldsymbol{\alpha},\pi)=\sum_{i=1}^{N}\log\left[\sum_{j=1}^{k}\pi_{j}f(\boldsymbol{x_{i}}|\boldsymbol{\alpha_{j}})\right]\leavevmode\nobreak\
.$ (3)

Latent variables $Z_{i}$’s are introduced, which are categorical variables
taking on values $1,\ldots,k$ with probabilities $\pi_{1},\ldots,\pi_{k}$ such
that $Pr(\boldsymbol{X_{i}}|Z_{i}=j)=f(\boldsymbol{x_{i}})$, $j=1,\ldots,k$.
The cluster membership probabilities of data point $i$ for cluster $j$ can be
obtained by $\gamma_{ij}$, where,

$\gamma_{ij}(x_{i})=Pr(Z_{i}=j|\boldsymbol{x_{i}},\boldsymbol{\alpha_{j}})=\frac{\pi_{j}f(\boldsymbol{x_{i}}|\boldsymbol{\alpha_{j}})}{\sum_{j=1}^{k}\pi_{j}f_{(}\boldsymbol{x_{i}}|\boldsymbol{\alpha_{j})}}\leavevmode\nobreak\
.$ (4)

Hard EM maximizes the classification log likelihood; it applies a delta
function approximation to the posterior probabilities
$Pr(Z_{i}=j|X=x,\alpha)$, where $Z_{i},i=1,\ldots,N$ are the latent variables
representing class labels. The approximation changes the E step as follows,

$Pr(Z_{i}=j|\boldsymbol{x_{i}},\boldsymbol{\alpha_{j}})\approx\mathbb{I}(j=z_{i}^{*}),$
(5)

where, $z_{i}^{*}=\operatorname*{argmax}\limits_{j}\gamma_{ij}$.
$\gamma_{ij}$’s are nothing but the responsibilities (probabilities) for each
data point that belongs to different clusters. In other words, it introduces a
classification step, where all the data points are classified into different
clusters based on the posterior probability. Let $N_{j}$ be the number of data
points in cluster $j$. Then, the ML estimates of $\pi_{j}$ is obtained as
$\displaystyle\frac{N_{j}}{N}$. And estimates of $\boldsymbol{\alpha_{j}}$’s
are obtained by finding the MLE using $N_{j}$ data points in cluster $j$.
However, as the solution is not analytically tractable, some numerical methods
need to be used. It is generally done by fixed-point iteration. The iterative
equation that needs to be solved is given by

$\Psi(\alpha_{jm}^{new})=\Psi(\sum_{m=1}^{p}\alpha_{jm}^{old})+\frac{1}{N_{j}}\sum_{i=1}^{N_{j}}\log(x_{im})\leavevmode\nobreak\
,$ (6)

where, $\Psi$ is the digamma function. The inverse of diggama function is used
to get the estimates of $\alpha_{jm}$. Now, let’s explore various
approximations to compute the KL distance for DMM.

###  2.2 KL Divergence: Monte Carlo Approach

Monte Carlo sampling is frequently employed by researchers to compute KL
divergence. Let $f_{a}$ and $g_{b}$ be two DMMs such that,

$\displaystyle f_{a}$
$\displaystyle=f(\boldsymbol{x})=\sum_{a}\pi_{a}Dir(\boldsymbol{x};\boldsymbol{\alpha}_{a})$
$\displaystyle g_{b}$
$\displaystyle=g(\boldsymbol{x})={\sum}_{b}\omega_{b}Dir(\boldsymbol{x};\boldsymbol{\alpha}_{b})\leavevmode\nobreak\
.$

The approach involves drawing a sample $\boldsymbol{x_{i}}$ from the
probability density function $f_{a}$ such that,

$E_{f_{a}}[\log
f_{a}(\boldsymbol{x_{i}})/g_{b}(\boldsymbol{x_{i}})]=D(f_{a}\|g_{b})\leavevmode\nobreak\
.$

Utilizing $N$ independent and identically distributed samples
$\\{\boldsymbol{x_{i}}\\}_{i=1}^{N}$, the estimation is expressed as follows:

$D_{{\rm MC}}(f_{a}\|g_{b})=\frac{1}{N}\sum_{i=1}^{N}\log
f_{a}(x_{i})/g_{b}(x_{i})\rightarrow D(f_{a}\|g_{b})\leavevmode\nobreak\ ,$
(7)

as $N\rightarrow\infty$. The variance of the estimation error is
$\frac{1}{N}{\rm Var}_{f_{a}}[\log f_{a}/g_{b}]$. To compute $D_{{\rm
MC}}(f\|g)$, it is necessary to generate the independent and identically
distributed samples $\\{\boldsymbol{x_{i}}\\}_{i=1}^{N}$ from $f_{a}$. Drawing
a sample $\boldsymbol{x_{i}}$ from the DMM $f_{a}$ involves initially drawing
a discrete sample $a_{i}$ according to the probabilities $\pi_{a}$.
Subsequently, a continuous sample $\boldsymbol{x_{i}}$ is drawn from the
resulting Dirichlet component $f_{a_{i}}(\boldsymbol{x})$.

The Monte Carlo method represents a convergent technique. It satisfies the
similarity property; however, the positivity property does not hold. The
identification property is likely to fail only under highly artificial
circumstances and with very low probability.

###  2.3 KL Divergence: Variational Approach

Now, we delve into the variational approximation to compute the KL divergence
for Dirichlet Mixture Models (DMM). Let us denote $\left\langle.\right\rangle$
as the inner product. In our context
$\left\langle.\right\rangle_{\boldsymbol{x}}$ means expectation with respect
to x.

######  Theorem 1.

Let $\boldsymbol{X}$ be an $p\times 1$ random vector. Assume two Dirichlet
distributions $u$ and $v$ specifying the probability distribution of
$\boldsymbol{X}$ as,

$\begin{split}u:\;\boldsymbol{X}&\sim\mathrm{Dir}(\alpha_{11},\ldots,\alpha_{1p})\\\
v:\;\boldsymbol{X}&\sim\mathrm{Dir}(\alpha_{21},\ldots,\alpha_{2p})\;.\end{split}$

Then, the Kullback-Leibler divergence of $u$ from $v$ is given by,

$\mathrm{D}(u\,||\,v)=\log\frac{\Gamma\left(\sum_{i=1}^{p}\alpha_{1i}\right)}{\Gamma\left(\sum_{i=1}^{p}\alpha_{2i}\right)}+\sum_{i=1}^{p}\log\frac{\Gamma(\alpha_{2i})}{\Gamma(\alpha_{1i})}+\sum_{i=1}^{p}\left(\alpha_{1i}-\alpha_{2i}\right)\left[\psi(\alpha_{1i})-\psi\left(\sum_{i=1}^{p}\alpha_{1i}\right)\right]\;.$
(8)

###### Proof.

For Dirichlet distributions KL divergence is,

$\begin{split}\mathrm{D}(u\,||\,v)&=\int_{\mathcal{\boldsymbol{X}}^{p}}\mathrm{Dir}(\boldsymbol{x};\alpha_{11},\ldots,\alpha_{1p})\,\log\frac{\mathrm{Dir}(\boldsymbol{x};\alpha_{11},\ldots,\alpha_{1p})}{\mathrm{Dir}(\boldsymbol{x};\alpha_{21},\ldots,\alpha_{2p})}\,\mathrm{d}\boldsymbol{x}\\\
&=\left\langle\log\frac{\mathrm{Dir}(x;\alpha_{11},\ldots,\alpha_{1p})}{\mathrm{Dir}(x;\alpha_{21},\ldots,\alpha_{2p})}\right\rangle_{u(\boldsymbol{x})}\end{split}$
(9)

We can do some algebraic manipulations using the probability density function
of the Dirichlet distribution to get,

$\begin{split}\mathrm{D}(u\,||\,v)&=\left\langle\log\frac{\frac{\Gamma\left(\sum_{i=1}^{p}\alpha_{1i}\right)}{\prod_{i=1}^{p}\Gamma(\alpha_{1i})}\,\prod_{i=1}^{p}{x_{i}}^{\alpha_{1i}-1}}{\frac{\Gamma\left(\sum_{i=1}^{p}\alpha_{2i}\right)}{\prod_{i=1}^{p}\Gamma(\alpha_{2i})}\,\prod_{i=1}^{p}{x_{i}}^{\alpha_{2i}-1}}\right\rangle_{u(\boldsymbol{x})}\\\
&=\left\langle\log\left(\frac{\Gamma\left(\sum_{i=1}^{p}\alpha_{1i}\right)}{\Gamma\left(\sum_{i=1}^{p}\alpha_{2i}\right)}\cdot\frac{\prod_{i=1}^{p}\Gamma(\alpha_{2i})}{\prod_{i=1}^{p}\Gamma(\alpha_{1i})}\cdot\prod_{i=1}^{p}{x_{i}}^{\alpha_{1i}-\alpha_{2i}}\right)\right\rangle_{u(\boldsymbol{x})}\\\
&=\left\langle\log\frac{\Gamma\left(\sum_{i=1}^{p}\alpha_{1i}\right)}{\Gamma\left(\sum_{i=1}^{p}\alpha_{2i}\right)}+\sum_{i=1}^{p}\log\frac{\Gamma(\alpha_{2i})}{\Gamma(\alpha_{1i})}+\sum_{i=1}^{p}(\alpha_{1i}-\alpha_{2i})\cdot\log(x_{i})\right\rangle_{u(\boldsymbol{x})}\\\
&=\log\frac{\Gamma\left(\sum_{i=1}^{p}\alpha_{1i}\right)}{\Gamma\left(\sum_{i=1}^{p}\alpha_{2i}\right)}+\sum_{i=1}^{p}\log\frac{\Gamma(\alpha_{2i})}{\Gamma(\alpha_{1i})}+\sum_{i=1}^{p}(\alpha_{1i}-\alpha_{2i})\cdot\left\langle\log
x_{i}\right\rangle_{u(\boldsymbol{x})}\;.\end{split}$

Moreover,

$\left\langle\log\boldsymbol{x}_{i}\right\rangle=\psi(\alpha_{i})-\psi\left(\sum_{i=1}^{p}\alpha_{i}\right)\;.$

Thus the Kullback-Leibler divergence of $u$ from $v$ becomes:

$\mathrm{D}(u\,||\,v)=\log\frac{\Gamma\left(\sum_{i=1}^{p}\alpha_{1i}\right)}{\Gamma\left(\sum_{i=1}^{p}\alpha_{2i}\right)}+\sum_{i=1}^{p}\log\frac{\Gamma(\alpha_{2i})}{\Gamma(\alpha_{1i})}+\sum_{i=1}^{p}\left(\alpha_{1i}-\alpha_{2i}\right)\left[\psi(\alpha_{1i})-\psi\left(\sum_{i=1}^{p}\alpha_{1i}\right)\right]\;.$

∎

######  Proposition 1.

Let $f_{a}$ and $g_{b}$ be two DMMs such that,

$\displaystyle f_{a}$
$\displaystyle=f(x)=\sum_{a}\pi_{a}Dir(\boldsymbol{x};\boldsymbol{\alpha}_{a})$
$\displaystyle g_{b}$
$\displaystyle=g(x)={\sum}_{b}\omega_{b}Dir(\boldsymbol{x};\boldsymbol{\alpha}_{b})\leavevmode\nobreak\
.$

Then using a variational approach, an approximated KL divergence can be
expressed as,

$D_{{\rm var}{\rm ia}{\rm ti}{\rm on}{\rm
al}}(f\|g)=\sum_{a}\pi_{a}\log{\sum_{a^{\prime}}\pi_{a^{\prime}}e^{-D(f_{a}\|f_{a^{\prime}})}\over\sum_{b}\omega_{b}e^{-D(f_{a}\|g_{b})}}.$
(10)

###### Proof.

Hershey and Olsen (2007) has provided a similar proof for Gaussian Mixture
Models. That approach can be adopted for DMMs as well. We can define the log-
likelihood,$L_{f}(g)$ as, $L_{f}(g)=\left\langle[\log
g(\boldsymbol{x})]\right\rangle_{f(\boldsymbol{x})}$. Then KL divergence can
be written in terms of log-likelihood in the following way.

$D(f\|g)=L_{f}(f)-L_{f}(g)\leavevmode\nobreak\ .$ (11)

We define variational parameters,$\phi_{b|a}>0$ such that
$\Sigma_{b}\phi_{b|a}=1$. Using Jensen’s inequality we can show,

$\displaystyle L_{f}(g)\ $ $\displaystyle\displaystyle\mathop{=}^{{\rm
def}}\left\langle\log g(x)\right\rangle_{f(\boldsymbol{x})}$
$\displaystyle=\left\langle\log\sum_{b}\omega_{b}g_{b}(x)\right\rangle_{f(\boldsymbol{x})}$
$\displaystyle=\left\langle\log\sum_{b}\phi_{b|a}{\omega_{b}g_{b}(x)\over\phi_{b|a}}\right\rangle_{f(\boldsymbol{x})}$
$\displaystyle\geq\left\langle\sum_{b}\phi_{b|a}\log{\omega_{b}g_{b}(x)\over\phi_{b|a}}\right\rangle_{f(\boldsymbol{x})}$
$\displaystyle\displaystyle\mathop{=}^{{\rm def}}\ {\cal
L}_{f}(g,\phi)\leavevmode\nobreak\ .$ (12)

The above is a lower bound on $L_{f}(g)$. We can get the best bound by
maximizing ${\cal L}_{f}(g,\phi)$ with respect to $\phi$. The maximum value is
obtained with:

$\hat{\phi}_{b|a}={\omega_{b}e^{-D(f_{a}\|g_{b})}\over\sum_{b}^{\prime}\pi_{b}^{\prime}e^{-D(f_{a}\|g_{b^{\prime}})}}\leavevmode\nobreak\
.$ (13)

Likewise, we define ,

${\cal L}_{f}(f,\psi)\ \ \displaystyle\mathop{=}^{{\rm def}}\ \
\left\langle\sum_{a^{\prime}}\psi_{a^{\prime}|a}\log{\pi_{a^{\prime}}f_{a^{\prime}}(x)\over\psi_{a^{\prime}|a}}\right\rangle_{f(\boldsymbol{x})}\leavevmode\nobreak\
.$ (14)

The optimal $\psi_{a^{\prime}|a}$ is given by,

$\hat{\psi}_{a^{\prime}|a}={\pi_{a^{\prime}}e^{-D(f_{a}\|f_{a^{\prime}})}\over\sum_{\hat{a}}\pi_{\hat{a}}e^{-D(f_{a}\|f_{\hat{a}})}}$
(15)

Now, like equation 11, we can define $D_{{\rm var}{\rm ia}{\rm ti}{\rm on}{\rm
al}}(f\|g)={\cal L}_{f}(f,\hat{\psi})-{\cal L}_{f}(g,\hat{\phi})$. After
substituting $\hat{\phi}_{b|a}$ and $\hat{\psi}_{a^{\prime}|a}$, we finally
get,

$D_{{\rm var}{\rm ia}{\rm ti}{\rm on}{\rm
al}}(f\|g)=\sum_{a}\pi_{a}\log{\sum_{a^{\prime}}\pi_{a^{\prime}}e^{-D(f_{a}\|f_{a^{\prime}})}\over\sum_{b}\omega_{b}e^{-D(f_{a}\|g_{b})}}\leavevmode\nobreak\
.$

∎

$D_{{\rm va}{\rm ri}{\rm at}{\rm io}{\rm nal}}(f\|g)$ satisfies the self
similarity and self identification property. However, it may not always
satisfy the positivity property. In that case we may use the absolute value of
the divergence.

##  3 Results

We conducted extensive experiments using both simulated and real data sets,
including four simulated data sets (Data Set 1, Data Set 2, Data Set 3, Data
Set 4) drawn from different Dirichlet distributions and a real gene expression
data set (Modencodefly Data). Data Set 3, Data Set 4 and Modencodefly Data are
highdimensional in nature. Some details of the datasets are given in table 1.
The primary metrics for comparison are the Kullback-Leibler (KL) Divergence
and the time taken for each approach. For real data analysis, KL divergence is
typically not used as the true parameter values of the mixture distribution
are unknown. However, in cases where true class labels are available, it
becomes feasible to approximate the true parameter values. In our analysis,
the dataset is classified into the true clusters, and for each cluster, we fit
a Dirichlet distribution to obtain the approximated true $\boldsymbol{\alpha}$
parameter values. The ratio of the number of data points in each cluster to
the total number of observations serves as the estimates for the true
$\boldsymbol{\pi}$ values. This approach enables us to leverage the known
class labels to derive estimations for the true parameters of the mixture
distribution, facilitating a more accurate assessment of KL divergence in
real-world datasets.

Table 1: Details of the data sets

Data Set Sample Size ($N$) Dimension($p$) No. of Clusters ($k$) Data Set 1
4200 3 3 Data Set 2 4200 3 3 Data Set 3 230 75 6 Data Set 4 260 90 6
Modencodefly 147 75 6

The Modencodefly dataset, derived from a comprehensive study of Drosophila
melanogaster (Graveley et al., 2011), serves as a rich resource for exploring
transcriptome dynamics throughout development. This dataset, obtained through
RNA-Seq, tiling microarrays, and cDNA sequencing across 30 distinct
developmental stages, contains 111,195 newly identified elements, including
genes, transcripts, exons, splicing events, and RNA editing sites. In our
study, we utilized the normalized gene profiles (Rau and Maugis-Rabusseau,
2018) derived from this dataset, transforming raw read counts into
compositional data. After filtering for genes exhibiting low variation among
different biological samples, we retained 75 genes for our analysis. The
Modencodefly data, available from ReCount (Frazee et al., 2011), provides a
high-resolution view of transcriptome dynamics and serves as a valuable
resource for investigating gene expression patterns in Drosophila
melanogaster.

Table 2: Comparison of KL Divergence and Time Taken for Different Data Sets

Data Set Metric Variational Monte Carlo $n=10000$ Monte Carlo $n=100000$ Monte
Carlo $n=1000000$ Data Set 1 KL Div. 0.0017 0.0024 0.0021 0.0017 Time 0.0004
1.7 17.3391 181.3682 Data Set 2 KL Div. 0.0022 0.0024 0.0024 0.0022 Time
0.0005 1.7283 17.4883 180.9482 Data Set 3 KL Div. 6.2723 6.1737 6.1815 6.1824
Time 0.0014 3.6702 37.4568 384.2183 Data Set 4 KL Div. 8.1268 7.7653 7.8025
7.8277 Time 0.0015 3.9773 40.5630 418.5154 Modencodefly KL Div. 34.7652
33.5581 33.5471 33.5579 Time 0.0014 3.6353 36.0671 375.7836

From Table 2, it is evident that our proposed variational approach
consistently provided faster solutions for all data sets. Notably, the Monte-
Carlo method yielded KL divergence values closer to variational method as we
increase the number of samples, which increases the execution time to a great
extent and indicating that our variational solution is more robust and
applicable for practical applications.

##  4 Conclusion

In this study, we addressed the challenging task of efficiently estimating
Kullback-Leibler (KL) divergence within the context of Dirichlet Mixture
Models (DMM). KL Divergence is a fundamental measure for quantifying the
statistical distance between probability distributions, playing a pivotal role
in various fields. Our investigation focused on the utility of DMMs, a
powerful statistical tool for clustering compositional data.

Despite the analytical tractability of KL Divergence for Dirichlet
distributions, extending this to DMMs proved to be a formidable challenge.
Past research primarily relied on Monte Carlo methods for approximating KL
divergence. However, these methods are computationally demanding and time-
consuming, presenting significant resource challenges. In response to these
challenges, we proposed a novel variational approach for efficiently
estimating KL divergence in DMMs. Our method provides a closed-form solution,
markedly enhancing computational efficiency compared to existing techniques.
This advancement facilitates swift model comparisons and robust evaluation of
estimation quality.

We validated our approach using a comprehensive set of both simulated and
real-world datasets, which includes a gene expression data. The results,
meticulously detailed in Table 2, underscore the superior efficiency and
accuracy of our novel variational approach when juxtaposed with traditional
Monte Carlo-based methods. Intriguingly, our analysis reveals a nuanced trend:
as the number of random samples escalates, the Monte Carlo method converges to
a KL divergence value remarkably proximate to our proposed variational
approach. However, this convergence comes at a significant cost - the
execution time expands exponentially with larger sample sizes. Consequently,
our study compellingly advocates for the adoption of the variational approach,
offering not only reliable KL divergence but also unparalleled computational
expediency.

In conclusion, our variational approach offers a transformative solution for
the efficient estimation of KL divergence in Dirichlet Mixture Models. This
innovation opens avenues for more rapid and reliable exploration of diverse
DMM models, fostering advancements in statistical analyses of compositional
data across various domains.

## References

  * Csiszár [1975] Imre Csiszár.  I-divergence geometry of probability distributions and minimization problems.  _The annals of probability_ , pages 146–158, 1975.
  * McLachlan and Peel [2000] Geoffrey McLachlan and David Peel.  _Finite Mixture Models_.  Wiley, September 2000.  ISBN 9780471721185.  doi: 10.1002/0471721182.  URL http://dx.doi.org/10.1002/0471721182.
  * Youbi et al. [2016] Zineb Youbi, Larbi Boubchir, Meriem D Bounneche, Arab Ali-Chérif, and Abdelhani Boukrouche.  Human ear recognition based on multi-scale local binary pattern descriptor and kl divergence.  In _2016 39th International conference on telecommunications and signal processing (TSP)_ , pages 685–688. IEEE, 2016.
  * Liu and Yao [2022] Suli Liu and Wu Yao.  Prediction of lung cancer using gene expression and deep learning with kl divergence gene selection.  _BMC bioinformatics_ , 23(1):175, 2022.
  * Pal and Heumann [2022] Samyajoy Pal and Christian Heumann.  Clustering compositional data using dirichlet mixture model.  _Plos one_ , 17(5):e0268438, 2022.
  * Aitchison [1982] John Aitchison.  The statistical analysis of compositional data.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 44(2):139–160, 1982.
  * Ma et al. [2014] Zhanyu Ma, Pravin Kumar Rana, Jalil Taghia, Markus Flierl, and Arne Leijon.  Bayesian estimation of dirichlet mixture model with variational inference.  _Pattern Recognition_ , 47(9):3143–3157, 2014.
  * Dempster et al. [1977] Arthur P Dempster, Nan M Laird, and Donald B Rubin.  Maximum likelihood from incomplete data via the em algorithm.  _Journal of the royal statistical society: series B (methodological)_ , 39(1):1–22, 1977.
  * Celeux and Govaert [1992] Gilles Celeux and Gérard Govaert.  A classification em algorithm for clustering and two stochastic versions.  _Computational statistics & Data analysis_, 14(3):315–332, 1992.
  * Hershey and Olsen [2007] John R Hershey and Peder A Olsen.  Approximating the kullback leibler divergence between gaussian mixture models.  In _2007 IEEE International Conference on Acoustics, Speech and Signal Processing-ICASSP’07_ , volume 4, pages IV–317. IEEE, 2007.
  * Graveley et al. [2011] Brenton R Graveley, Angela N Brooks, Joseph W Carlson, Michael O Duff, Jane M Landolin, Li Yang, Carlo G Artieri, Marijke J Van Baren, Nathan Boley, Benjamin W Booth, et al.  The developmental transcriptome of drosophila melanogaster.  _Nature_ , 471(7339):473–479, 2011.
  * Rau and Maugis-Rabusseau [2018] Andrea Rau and Cathy Maugis-Rabusseau.  Transformation and model choice for rna-seq co-expression analysis.  _Briefings in bioinformatics_ , 19(3):425–436, 2018.
  * Frazee et al. [2011] Alyssa C Frazee, Ben Langmead, and Jeffrey T Leek.  Recount: a multi-experiment resource of analysis-ready rna-seq gene count datasets.  _BMC bioinformatics_ , 12:1–5, 2011.


# On Zero-Sum Two Person Perfect Information Stochastic Games

K. G. Bakshia,1 111Corresponding author, S. Sinhaa,2
_a,1 Department Of Mathematics, Jadavpur University, Kolkata, 700 032, India._
_a,2 Department Of Mathematics, Jadavpur University, Kolkata , 700 032,
India._
_a,1 Email<EMAIL_ADDRESS>
_a,2 Email<EMAIL_ADDRESS>
###### Abstract

A zero-sum two person Perfect Information Stochastic game (PISG) under
limiting average payoff has a value and both the maximiser and the minimiser
have optimal pure stationary strategies. Firstly we form the matrix of
undiscounted payoffs corresponding to each pair of pure stationary strategies
(for each initial state) of the two players and prove that this matrix has a
pure saddle point. Then by using the results by Derman [1] we prove the
existence of optimal pure stationary strategy pair of the players. A crude but
finite step algorithm is given to compute such an optimal pure stationary
strategy pair of the players.

Keywords: Stochastic games, Markov Decision Processes, Perfect Information,
Stationary Strategies, Linear Programming.

AMS subject classifications: 90C40, 91A15, 90C05.

##  1 Introduction

Stochastic games are generalizations of Markov decision processes (MDPs) to
the case of two or more players. Shapley (1953) [12] introduced ’Stochastic
games’ in his paper, which is known as Markov games these days. If two players
play a matrix game repeatedly over the infinite time horizon and the limiting
average payoff is considered, then the value of this infinitely repeated game
coincides with the value of the one shot game (by Folk Theorem [3]). Shapley
[12] introduced the idea of not playing the same matrix game everyday (i.e.,
in every stage of the game), but playing one among finitely many matrix games,
with a motion among them governed by the present game and the actions chosen
there in such a manner that the game is certain to stop in finite time. Then
the payoffs of the players can be formulated as the ratio of two bilinear
forms. Neumann [10] established the minimax theorem for such games and Loomis
[9] gave an elementary proof of this theorem. The case of non-terminating
limiting average Stochastic games were studied by Gillette [4], Hoffman and
Karp [5] . By undiscounted pay-off we mean limiting average pay-off in this
paper. Liggett and Lippman [8] previously proved the existence of pure
stationary optimal strategy pair of the players in an undiscounted perfect
information stochastic game. We propose an alternative proof (with less
complexity) of the result by Liggett and Lippman [8]. By forming the matrix of
undiscounted payoffs corresponding to each pair of pure stationary strategies
(for each initial state) of the two players we prove that this matrix has a
pure saddle point, which is esentially a pure semi-stationary strategy pair of
the players. Then we prove the existence of optimal pure stationary strategy
pair of the players by using the results by Derman [1]. We consider the
policy-improvement algorithm to compute optimal pure stationary strategy pair
of the players. This is a best response algorithm, in which each player looks
for his own Blackwell optimal strategy. It is obvious that this is a finite
step algorithm and it terminates in finite time by the conjecture $8.1$ of
Raghavan and Syed (2002) [11]. The paper is organized as follows. Section 2
contains definitions and properties of an undiscounted two person zero-sum
Stochastic games considered under limiting average pay-off. Section 3 contains
main result of this paper. In section 4 we propose a policy improvement
algorithm to compute an optimal stationary strategy pair for the players of
such perfect information undiscounted Stochastic games. Section 5 contains
some numerical examples illustrating our theorem and proposed algorithm.

##  2 Preliminaries

###  2.1 Finite two preson zero-sum Stochastic games

A zero-sum two person finite stochastic game is described by a collection of
five objects $\Gamma=<S,\\{A(s):s\in S\\},\\{B(s):s\in S\\},q,r>$, where
$S=\\{0,1,\cdots,z\\}$ is the finite non-empty state space and
$A(s)=\\{0,1,\cdots,m_{s}\\},B(s)=\\{0,1,\cdots,n_{s}\\}$ are respectively the
non-empty sets of admissible actions of the players I and II respectively in
the state $s$. Let us denote $K=\\{(s,i,j):s\in S,i\in A(s),j\in B(s)\\}$ to
be the set of admissible triplets. For each $(s,i,j)\in K$, we denote $q(.\mid
s,i,j)$ to be the transition law of the game. Finally $r$ is the real valued
functions on $K$, which represents the immediate (expected) reward for the
player-I (whereas -$r$ is the reward for the player-II). Let us consider
player I as the maximiser and player II as the minimiser in the zero-sum two
person stochastic game.
The Stochastic game over infinite time is played as follows. At the $0$th
decision epoch, the game strats at $s_{0}\in S$ and the players I and II
simultaneously and independently choose actions $i_{0}\in A(s_{0})$ and
$j_{0}\in B(s_{0})$ respectively. Consequently player I and II get immediate
rewards $r(s_{0},i_{0},j_{0})$ and $-r(s_{0},i_{0},j_{0})$ respectively and
the game moves to the state $s_{1}$ with probability $q(s_{1}\mid
s_{0},i_{0},j_{0})$. After reaching the state $s_{1}$ on the next decision
epoch, the game is repeated over infinite time with the state $s_{0}$ replaced
by $s_{1}$. Shapley extended the idea of defining SGs where
$\sum_{s^{{}^{\prime}}\in S}q(s^{{}^{\prime}}\mid s,i,j)<1$ for all
$(s,i,j)\in K$ and the play terminates with probability
$1-\sum_{s^{{}^{\prime}}\in S}q(s^{{}^{\prime}}\mid s,i,j)<1$. Such games are
called ‘stopping SGs’. The ‘non-stopping SGs’ are those where
$\sum_{s^{{}^{\prime}}\in S}q(s^{{}^{\prime}}\mid s,i,j)=1$ for all
$(s,i,j)\in K$, i.e., the play never terminates.

By a strategy (behavioural) $\pi_{1}$ of the player I, we mean a sequence
$\\{(\pi_{1})_{n}(.\mid hist_{n})\\}_{n=1}^{\infty}$, where $(\pi_{1})_{n}$
specifies which action is to be chosen on the $n$-th decision epoch by
associating with each history $hist_{n}$ of the system up to $n$th decision
epoch (where $hist_{n}$=$(s_{0},a_{0},b_{0},s_{1},a_{1}\\\
,b_{1}\cdots,s_{n-1},a_{n-1},b_{n-1},s_{n})$ for $n\geq 2$, $hist_{1}=(s_{0})$
and $(s_{k},a_{k},j_{k})\in K$ are respectively the state and actions of the
players at the $k$-th decision epoch) a probability distribution
$(\pi_{1})_{n}(.\mid hist_{n})$ on $A(s_{n})$. Behavioural strategy $\pi_{2}$
for player II can be defined analogously. Generally by any unspecified
strategy, we mean behavioural strategy here. We denote $\Pi_{1}$ and $\Pi_{2}$
to be the sets of strategy (behavioural) spaces of the players I and II
respectively. A strategy
$f^{{}^{\prime}}=\\{f^{{}^{\prime}}_{n}\\}_{n=1}^{\infty}$ for the player I is
called semi-Markov if for each $n$, $f^{{}^{\prime}}_{n}$ depends on
$s_{1},s_{n}$ and the decision epoch number $n$. Similarly we can define a
semi-Markov strategy
$g^{{}^{\prime}}=\\{g^{{}^{\prime}}_{n}\\}_{n=1}^{\infty}$ for the player II.
A strategy $\pi_{1}=\\{\pi_{1n}\\}_{n=1}^{\infty}$ is called a stationary
strategy if $\exists$ a map
$f:S\rightarrow\mathbb{P}(A)=\\{\mathbb{P}(A(s)):s\in S\\}$, where
$\mathbb{P}(A(s))$ is the set of probability distribution on $A(s)$ such that
$\pi_{1n}=f$ for all $n$ and $f(s)\in\mathbb{P}(A(s))$. A stationary strategy
for player I is defined as $z$ tuple $f=(f(1),f(2),\cdots,f(z))$, where each
$f(s)$ is the probability distribution on $A(s)$ given by
$f(s)=(f(s,1),f(s,2),\cdots,f(s,m_{s}))$. $f(s,i)$ denotes the probability of
choosing action $i$ in the state $s$ by player-I. By similar manner, one can
define a stationary strategy $g$ for player II as $g=(g(1),g(2),\cdots,g(z))$
where each $g(s)$ is the probability distribution on $B(s)$. Let us denote
$F_{1}^{s}$ and $F_{2}^{s}$ to be the set of stationary strategies for player
I and II respectively. A semi-stationary strategy is a semi-Markov strategy
which is independent of the decision epoch $n$, i.e., for a initial state
$s_{1}$ and present state $s_{2}$, if a semi-Markov strategy
$f^{{}^{\prime}}(s_{1},s_{2},n)$ turns out to be independent of $n$, then we
call it a semi-stationary strategy. Let us denote $\xi_{1}$ and $\xi^{2}$ to
be the set of semi-stationary strategies for player-I and II respectively.
A stationary strategy is called pure if any player selects a particular action
with probability $1$ while visiting a state $s$. We denote $F_{1}^{sp}$ and
$F_{2}^{sp}$ to be the set of pure stationary strategies of the players I and
II respectively. Also $\xi_{1}^{sp}$ and $\xi_{2}^{sp}$ are denoted as the set
of pure semi-stationary strategies for the player-I and II respectively.
Definition 1 A zero-sum two person SG $\Gamma=<S,\\{A(s):s\in
S\\},\\{B(s):s\in S\\},q,r>$ is called a perfect information stochastic game
(PISG) if the following properties hold
(i)$S=S_{1}\cup S_{2},S_{1}\cap S_{2}=\phi$.
(ii)$\mid B(s)\mid=1$, for all $s\in S_{1}$, i.e., on $S_{1}$ player-II is a
dummy.
(iii)$\mid A(s)\mid=1$, for all $s\in S_{2}$, i.e., on $S_{2}$ player-I is a
dummy.

###  2.2 Undiscounted zero-sum two person stochastic games

Let $(X_{1},A_{1},B_{1},X_{2},A_{2},B_{2}\cdots)$ be a co-ordinate sequence in
$S\times(A\times B\times S)^{\infty}$. Given behavioural strategy pair
$(\pi_{1},\pi_{2})\in\Pi_{1}\times\Pi_{2}$, initial state $s\in S$, there
exists a unique probability measure $P_{\pi_{1}\pi_{2}}(.\mid X_{0}=s)$ (hence
an expectation $E_{\pi_{1}\pi_{2}}(.\mid X_{0}=s)$) on the product $\sigma$\-
field of $S\times(A\times B\times S)^{\infty}$ by Kolmogorov’s extension
theorem. For a pair of strategies $(\pi_{1},\pi_{2})\in\Pi_{1}\times\Pi_{2}$
for the players I and II respectively, the limiting average (undiscounted)
pay-off for player I, starting from a state $s\in S$ is defined by:

$\phi(s,\pi_{1},\pi_{2})=\liminf_{n\to\infty}\frac{1}{n}E_{\pi_{1}\pi_{2}}\sum_{m=1}^{n}[r(X_{m},A_{m},B_{m})\mid
X_{0}=s]$ (2.1)

Alternatively, for any pair of stationary strategies $(f_{1},f_{2})\in
F_{1}^{s}\times F_{2}^{s}$ of player I and II, we write the undiscounted pay-
off for player I as:

$\phi(s,f_{1},f_{2})=\liminf_{n\to\infty}\frac{1}{n}\sum_{m=1}^{n}r^{m}(s,f_{1},f_{2})$
(2.2)

for all $s\in S$. Where $r^{m}(s,f_{1},f_{2})$ is the respectively the
expected reward for player I at the $m$ th decision epoch, when player I
chooses $f_{1}$ and player II chooses $f_{2}$ respectively and the initial
state is $s$.
Definition 2 For a pair of strategies $(f_{1},f_{2})\in F_{1}^{s}\times
F_{2}^{s}$, we define the transition probability matrix by:

$Q(f_{1},f_{2})=[q(s^{{}^{\prime}}\mid
s,f_{1}(s),f_{2}(s))]_{s,s^{{}^{\prime}}=1}^{z}$,

where $q(s^{{}^{\prime}}\mid s,f_{1}(s),f_{2}(s))=\sum_{i\in A(s)}\sum_{j\in
B(s)}q(s^{{}^{\prime}}\mid s,i,j)f_{1}(s,i)f_{2}(s,j)$ is the probability is
that the system jumps to the state $s^{{}^{\prime}}$ from given state $s$ when
the players play the stationary strategies $f_{1}$ and $f_{2}$.
Lemma 1(Kemeney and Snell, 1976, [7]) Let $Q$ be any $z\times z$ Markov
matrix, then the sequence
$\lim_{n\to\infty}\frac{1}{n+1}\sum_{m=0}^{n}Q^{m}(f_{1},f_{2})$ converges as
$n\rightarrow\infty$ to a Markov matrix $Q^{\ast}$ (the cesaro limiting
matrix) such that $QQ^{\ast}=Q^{\ast}Q=Q^{\ast}Q^{\ast}=Q^{\ast}$.
For each $(f_{1},f_{2})\in F_{1}\times F_{2}$, we define
$r(f_{1},f_{2})=[r(s,f_{1},f_{2})]_{z\times 1}$ as the expected reward, where
for each $s\in S$,

$r(s,f_{1},f_{2})=\sum_{i\in A(s)}\sum_{j\in
B(s)}r(s,i,j)f_{1}(s,i)f_{2}(s,j)$.

Now we have the following result:
Proposition 1 For each player of pure stationary strategies $(f_{1},f_{2})\in
F_{1}^{sp}\times F_{2}^{sp}$,

$\phi(s,f_{1},f_{2})=[Q^{\ast}(f_{1},f_{2})r(f_{1},f_{2})](s)\forall s\in S.$

Where $Q^{\ast}(f_{1},f_{2})$ is the cesaro limiting matrix of
$Q(f_{1},f_{2})$.
Definition 3 A zero-sum two person undiscounted stochastic game is said to
have a value vector $\phi=[\phi(s)]_{N\times 1}$ if
$\sup_{\pi_{1}\in\Pi_{1}}\inf_{\pi_{2}\in\Pi_{2}}\phi(s,\pi_{1},\pi_{2})=\phi(s)=\inf_{\pi_{2}\in\Pi_{2}}\sup_{\pi_{1}\in\Pi_{1}}\phi(s,\pi_{1},\pi_{2})$
for all $s\in S$. A pair of strategies
$(\pi_{1}^{\ast},\pi_{2}^{\ast})\in\Pi_{1},\times\Pi_{2}$ is said to be an
optimal strategy pair for the players if
$\phi(s,\pi_{1}^{\ast},\pi_{2})\geq\phi(s)\geq\phi(s,\pi_{1},\pi_{2}^{\ast})$
for all $s\in S$ and all $(\pi_{1},\pi_{2})\in\Pi_{1}\times\Pi_{2}$. A finite
(state and action spaces) Markov decision process is defined by a collection
of four objects $\hat{\Gamma}=<S,\hat{A}=\\{A(s):s\in S\\},\hat{q},\hat{r}>$,
where $S=\\{0,1,\cdots,z\\}$ is the finite state space,
$\hat{A}(s)=\\{1,2,\cdots,d\\}$ is the finite set of admissible actions in the
state $s$. $\hat{q}(s^{{}^{\prime}}\mid s,a)$ is the transition probabilty
(i.e., $\hat{q}(s^{{}^{\prime}}\mid s,a)\geq 0$ and $\sum_{s^{{}^{\prime}}\in
S}\hat{q}(s^{{}^{\prime}}\mid s,a)=1$) that the next state is
$s^{{}^{\prime}}$, where $s$ is the initial state and the decision maker
chooses action $a$ in the state $s$. The decision process proceeds over
infinite time just as stochastic game, where instead of two players we
consider a single decision maker. The definition of strategy spaces for the
decision maker is same as in the case of stochastic games. Let us denote
$\Pi$, $F$, $F_{s}$ as the set of behavioural, stationary, pure-stationary
strategies respectively of the decision maker. Let
$(X_{1},A_{1},X_{2},A_{2},\cdots)$ be a coordinate sequence in
$S\times(\hat{A}\times S)^{\infty}$. Given a behavioural strategy $\pi\in\Pi$,
initial state $s\in S$, there exists a unique probability measure
$P_{\pi}(.\mid X_{0}=s)$ (hence an expectation $E_{\pi}(.\mid X_{0}=s)$) on
the product $\sigma$\- field of $S\times(\hat{A}\times S)^{\infty}$ by
Kolmogorov’s extension theorem.

For a behavioural strategy $\pi\in\Pi$, the expected limiting average pay-off
is defined by

$\hat{\phi}(s,\pi)=\liminf_{n\to\infty}\frac{1}{n}\sum_{m=1}^{n}E_{\pi}[\hat{r}(X_{m},A_{m})\mid
X_{0}=s].$ (2.3)

for all $s\in S$.

##  3 Main result

Theorem 2 Any zero-sum two person undiscounted perfect information Stochastic
game has a solution in pure stationary strategies.

###### Proof.

Let $\Gamma=<S=S_{1}\cup S_{2},A=\\{A(s):s\in S_{1}\\},B=\\{B(s):s\in
S_{2}\\},q,r>$ be a zero-sum two person perfect information Stochastic game
under limiting average pay-off, where $S=\\{0,1,,\cdots,z\\}$ is the finite
state space. We assume that in $\mid S_{1}\mid$ number of states, player-II is
a dummy and for states $\\{\mid S_{1}\mid+1,\cdots,\mid S_{1}\mid+\mid
S_{2}\mid\\}$ player-I is a dummy. We assume that in this perfect information
game, each player has $d$ number of pure actions in each state where they are
non-dummy. Thus, player-I has $\mid S_{1}\mid.d$ number of pure actions
available in each state $s\in S$, where he/she is non-dummy and player-II has
$\mid S_{2}\mid.d$ number of pure actions where he/she is non-dummy in the
PISG $\Gamma$. Let us the consider the pay-off matrix

$A_{\mid S_{1}\mid.d\times\mid
S_{2}\mid.d}=\left[{\begin{array}[]{cccc}\phi(s,f_{0},g_{0})&\phi(s,f_{0},g_{1})&\cdots&\phi(s,f_{0},g_{\mid
S_{2}\mid.d})\\\
\phi(s,f_{1},g_{0})&\phi(s,f_{1},g_{1})&\cdots&\phi(s,f_{1},g_{\mid
S_{2}\mid.d})\\\ \vdots&\vdots&\ddots&\vdots\\\ \phi(s,f_{\mid
S_{1}\mid.d},g_{0})&\phi(s,f_{\mid S_{1}\mid.d},g_{1})&\cdots&\phi(s,f_{\mid
S_{1}\mid.d},g_{\mid S_{2}\mid.d})\\\ \end{array}}\right]$

Where $(f_{0},f_{1},\cdots,f_{\mid S_{1}\mid.d})$ and
$(g_{0},g_{1},\cdots,g_{\mid S_{2}\mid.d})$ are the pure stationary strategies
chosen by player-I and II repsectively. In order to prove the existence of a
pure semi-stationary strategy, we have to prove that this matrix has a pure
saddle point for each initial state $s\in S$. Now by Shapley [2], if A is the
matrix of a two-person zero-sum game and if every $2\times 2$ submatrix of $A$
has a saddle point, then A has a saddle point. So, we concentrate only on a
$2\times 2$ matrix and observe if it has a saddle point or not. We consider
the $2\times 2$ submatrix:

$\left[{\begin{array}[]{cccc}\phi(s,f_{i},g_{j})&\phi(s,f_{i},g_{j^{{}^{\prime}}})\\\
\phi(s,f_{i^{{}^{\prime}}},g_{j})&\phi(s,f_{i^{{}^{\prime}}},g_{j^{{}^{\prime}}})\\\
\end{array}}\right]$

Where $i^{{}^{\prime}},i\in\\{0,1,\cdots,\mid S_{1}\mid.d\\},(i\neq
i^{{}^{\prime}})$ and $j,j^{{}^{\prime}}\in\\{0,1,\cdots,\mid
S_{2}\mid.d\\},(j\neq j^{{}^{\prime}})$. Now, by suitably renumbering the
strategies, we can write the above sub-matrix as:

$\left[{\begin{array}[]{cccc}\phi(s,f_{1},g_{1})&\phi(s,f_{1},g_{2})\\\
\phi(s,f_{2},g_{1})&\phi(s,f_{2},g_{2})\\\ \end{array}}\right]$

Using the definition of $\phi(s,f_{1},f_{2})$ in section $2$, we get that

${\begin{array}[]{cc}\phi(s,f_{i},g_{j})=\sum_{s^{{}^{\prime}}\in
S}q^{\ast}(s^{{}^{\prime}}\mid s,f_{i},g_{j})r(s^{{}^{\prime}},f_{i},g_{j})\\\
=\sum_{t=1}^{S_{1}}[q^{\ast}(t\mid
s,f_{i.})r(t,f_{i.})]+\sum_{v=S_{1}+1}^{S_{1}+S_{2}}[q^{\ast}(v\mid
s,g_{.j})r(v,g_{.j})]\end{array}}$

Where

$f_{i}(s,.)=\left\\{\begin{array}[]{lr}f_{i.}(s,.)&s\in S_{1}\\\ 1&s\in
S_{2}\\\ \end{array}\right.$

and

$g_{j}(s,.)=\left\\{\begin{array}[]{lr}1&s\in S_{1}\\\ g_{.j}(s,.)&s\in
S_{2}\\\ \end{array}\right.$

We replace $\phi(s,f_{i},g_{j})$ by the expression above in the matrix $A$. We
consider the following two cases when $A$ can not have a pure saddle point.
Case-1: $\phi(s,f_{1},g_{1})$ is row minimum and column minimum,
$\phi(s,f_{1},g_{2})$ is row maximum and column maximum, $\phi(s,f_{2},g_{1})$
is row-maximum and column maximum and $\phi(s,f_{2},g_{2})$ is row-minimum and
column-minimum. These four conditions can be written as:
$\phi(s,f_{1},g_{1})<\phi(s,f_{1},g_{2})$,
$\phi(s,f_{1},g_{1})<\phi(s,f_{2},g_{1}),$
$\phi(s,f_{2},g_{2})<\phi(s,f_{2},g_{1})$,
$\phi(s,f_{2},g_{2})<\phi(s,f_{1},g_{2})$. Thus we get the following
inequalities:

$\displaystyle\sum_{t=1}^{S_{1}}[q^{\ast}(t\mid
s,f_{1.})r(t,f_{1.})]+\sum_{v=S_{1}+1}^{S_{1}+S_{2}}[q^{\ast}(v\mid
s,g_{.1})r(v,g_{.1})]$ (3.1) $\displaystyle<\sum_{t=1}^{S_{1}}[q^{\ast}(t\mid
s,f_{1.})r(t,f_{1.})]+\sum_{v=S_{1}+1}^{S_{1}+S_{2}}[q^{\ast}(v\mid
s,g_{.2})r(v,g_{.2})]$
$\displaystyle\noindent\sum_{t=1}^{S_{1}}[q^{\ast}(t\mid
s,f_{1.})r(t,f_{1.})]+\sum_{v=S_{1}+1}^{S_{1}+S_{2}}[q^{\ast}(v\mid
s,g_{.1})r(v,g_{.1})]$ (3.2) $\displaystyle<\sum_{t=1}^{S_{1}}[q^{\ast}(t\mid
s,f_{2.})r(t,f_{2.})]+\sum_{v=S_{1}+1}^{S_{1}+S_{2}}[q^{\ast}(v\mid
s,g_{.1})r(v,g_{.1})]$
$\displaystyle\noindent\sum_{t=1}^{S_{1}}[q^{\ast}(t\mid
s,f_{2.})r(t,f_{2.})]+\sum_{v=S_{1}+1}^{S_{1}+S_{2}}[q^{\ast}(v\mid
s,g_{.2})r(v,g_{.2})]$ (3.3) $\displaystyle<\sum_{t=1}^{S_{1}}[q^{\ast}(t\mid
s,f_{2.})r(t,f_{2.})]+\sum_{v=S_{1}+1}^{S_{1}+S_{2}}[q^{\ast}(v\mid
s,g_{.1})r(v,g_{.1})]$
$\displaystyle\noindent\sum_{t=1}^{S_{1}}[q^{\ast}(t\mid
s,f_{2.})r(t,f_{2.})]+\sum_{v=S_{1}+1}^{S_{1}+S_{2}}[q^{\ast}(v\mid
s,g_{.2})r(v,g_{.2})]$ (3.4) $\displaystyle<\sum_{t=1}^{S_{1}}[q^{\ast}(t\mid
s,f_{1.})r(t,f_{1.})]+\sum_{v=S_{1}+1}^{S_{1}+S_{2}}[q^{\ast}(v\mid
s,g_{.2})r(v,g_{.2})]$

Hence, $(3.1)$ yields

$\sum_{v=S_{1}+1}^{S_{1}+S_{2}}q^{\ast}(v\mid
s,g_{.2})r(v,g_{.2})-q^{\ast}(v\mid s,g_{.1})r(v,g_{.1})\textgreater 0$ (3.5)

$(3.3)$ yields

$\sum_{v=S_{1}+1}^{S_{1}+S_{2}}q^{\ast}(v\mid
s,g_{.1})r(v,g_{.1})-q^{\ast}(v\mid s,g_{.2})r(v,g_{.2})\textgreater 0$ (3.6)

From $(3.5)$ and $(3.6)$ we clearly get a contradiction. Now we consider the
next case:

Case-2: $\phi(s,f_{1},g_{1})$ is row maximum and column maximum,
$\phi(s,f_{1},g_{2})$ is row minimum and column minimum, $\phi(s,f_{2},g_{1})$
is row-minimum and column minimum and $\phi(s,f_{2},g_{2})$ is row-maximum and
column-maximum. These four conditions can be written as: These four conditions
can be written as: $\phi(s,f_{1},g_{1})>\phi(s,f_{1},g_{2})$,
$\phi(s,f_{1},g_{1})>\phi(s,f_{2},g_{1})$,
$\phi(s,f_{2},g_{2})>\phi(s,f_{2},g_{1})$,
$\phi(s,f_{2},g_{2})>\phi(s,f_{1},g_{2})$. We can re-write them as follows:

$\displaystyle\sum_{t=1}^{S_{1}}[q^{\ast}(t\mid
s,f_{1.})r(t,f_{1.})]+\sum_{v=S_{1}+1}^{S_{1}+S_{2}}[q^{\ast}(v\mid
s,g_{.1})r(v,g_{.1})]$ (3.7) $\displaystyle>\sum_{t=1}^{S_{1}}[q^{\ast}(t\mid
s,f_{1.})r(t,f_{1.})]+\sum_{v=S_{1}+1}^{S_{1}+S_{2}}[q^{\ast}(v\mid
s,g_{.2})r(v,g_{.2})]$
$\displaystyle\noindent\sum_{t=1}^{S_{1}}[q^{\ast}(t\mid
s,f_{1.})r(t,f_{1.})]+\sum_{v=S_{1}+1}^{S_{1}+S_{2}}[q^{\ast}(v\mid
s,g_{.1})r(v,g_{.1})]$ (3.8) $\displaystyle>\sum_{t=1}^{S_{1}}[q^{\ast}(t\mid
s,f_{2.})r(t,f_{2.})]+\sum_{v=S_{1}+1}^{S_{1}+S_{2}}[q^{\ast}(v\mid
s,g_{.1})r(v,g_{.1})]$
$\displaystyle\noindent\sum_{t=1}^{S_{1}}[q^{\ast}(t\mid
s,f_{2.})r(t,f_{2.})]+\sum_{v=S_{1}+1}^{S_{1}+S_{2}}[q^{\ast}(v\mid
s,g_{.2})r(v,g_{.2})]$ (3.9) $\displaystyle>\sum_{t=1}^{S_{1}}[q^{\ast}(t\mid
s,f_{2.})r(t,f_{2.})]+\sum_{v=S_{1}+1}^{S_{1}+S_{2}}[q^{\ast}(v\mid
s,g_{.1})r(v,g_{.1})]$
$\displaystyle\noindent\sum_{t=1}^{S_{1}}[q^{\ast}(t\mid
s,f_{2.})r(t,f_{2.})]+\sum_{v=S_{1}+1}^{S_{1}+S_{2}}[q^{\ast}(v\mid
s,g_{.2})r(v,g_{.2})]$ (3.10) $\displaystyle>\sum_{t=1}^{S_{1}}[q^{\ast}(t\mid
s,f_{1.})r(t,f_{1.})]+\sum_{v=S_{1}+1}^{S_{1}+S_{2}}[q^{\ast}(v\mid
s,g_{.2})r(v,g_{.2})].$

Hence, $(3.7)$ yields

$\sum_{v=S_{1}+1}^{S_{1}+S_{2}}q^{\ast}(v\mid
s,g_{.1})r(v,g_{.1})-q^{\ast}(v\mid s,g_{.2})r(v,g_{.2})\textgreater 0$ (3.11)

$(3.9)$ yields

$\sum_{v=S_{1}+1}^{S_{1}+S_{2}}q^{\ast}(v\mid
s,g_{.2})r(v,g_{.2})-q^{\ast}(v\mid s,g_{.1})r(v,g_{.1})\textgreater 0$ (3.12)

From $(3.11)$ and $(3.12)$ we clearly get a contradiction. Thus, every
$2\times 2$ submatrix has a pure saddle point and by Shapley [2], we claim
that the matrix $A$ has a pure saddle point, namely $(F^{\ast},G^{\ast})$. Now
$F^{\ast}=(f_{0},f_{1},\cdots,f_{t},\cdots,f_{z})$ and
$G^{\ast}=(g_{0},g_{1},\cdots,g_{t},\cdots,g_{z})$ where $f_{t}$ and $g_{t}$
are the pure stationary strategies for the initial state $t$ chosen by
player-I and II respectively. Now we prove the following lemma to prove the
existence of pure stationary strategy pair which is optimal for the players:

######  Lemma 3.1:

Let us fix an initial state $t\in S$ in the PISG $\Gamma$. Suppose
$(f_{t},g_{t})\in F_{1}^{sp}\times F_{2}^{sp}$ be an optimal pure stationary
strategy pair of the players satisfying:

$\phi(t,f_{t},g_{t})\leq\phi(t,f_{t},g)$ $\forall g\in F_{2}^{sp}$ and for
some initial state $t\in S$.

Let us denote $D_{t}$ to be the $t$-th row of the bi-matrix identifying the
strategy pair $(f_{t},g_{t})$, i.e.,
$D_{t}=((f_{t}(t,0),g_{t}(t,0))\cdots,(f_{t}(t,d),g_{t}(t,d))$, where $d$ is
the total number of pure actions in state $t$ for both the players. Then
$(f^{\ast},g^{\ast})\in F_{1}^{sp}\times F_{2}^{sp}$ is a pure stationary
strategy pair of the players identified by the bi-matrix $D^{\ast}$ having
$D_{t}$ as its $t$-th row. We can write the bi-matrix $D^{\ast}$ as:

$D^{\ast}_{(z+1)\times
d}=\left[{\begin{array}[]{cccc}(f_{0}(0,0),g_{0}(0,0))&(f_{0}(0,1),g_{0}(0,1))&\cdots&(f_{0}(0,d),g_{0}(0,d))\\\
(f_{1}(1,0),g_{1}(1,0)))&(f_{1}(1,1),g_{1}(1,1))&\cdots&(f_{1}(1,d),g_{1}(1,d))\\\
\vdots&\vdots&\ddots&\vdots\\\
(f_{z}(z,0),g_{z}(z,0))&(f_{z}(z,1),g_{z}(z,1))&\cdots&(f_{z}(z,d),g_{z}(z,d))\\\
\end{array}}\right]$

and the pair $(f^{\ast},g^{\ast})$ satisfies:

$\phi(t,f^{\ast},g^{\ast})\leq\phi(t,f^{\ast},g)\forall g\in
F_{2}^{sp},\forall t\in S.$ (3.13)

###### Proof.

For an initial state $t\in S(=\\{0,1,\cdots,z\\})$ and a pair of behavioural
strategy $(\pi_{1},\pi_{2})\in\Pi_{1}\times\Pi_{2}$ of the players, we
consider the $(z+1)d^{2}$ component vector:

$\xi_{n}^{\pi_{1}\pi_{2}}=\\{x_{n000}^{t},x_{n001}^{t},\cdots,x_{ns^{{}^{\prime}}ab}^{t},\cdots,x_{nzd^{2}}^{t}\\}$

where
$x_{ns^{{}^{\prime}}ab}^{t}=\frac{1}{n}\sum_{m=1}^{n}P_{\pi_{1}\pi_{2}}(X_{m}=s^{{}^{\prime}},A_{m}=a,B_{m}=b\mid
X_{0}=t)$. Let
$\xi^{\pi_{1}\pi_{2}}(t)=\lim_{n\to\infty}\xi_{n}^{\pi_{1}\pi_{2}}(t)$,
whenever the limit exists and
$\lim_{n\to\infty}x_{ns^{{}^{\prime}}ab}^{t}=x_{s^{{}^{\prime}}ab}^{t}$.
Denote $\Theta(\xi_{n}^{\pi_{1}\pi_{2}}(t))=\sum_{s^{{}^{\prime}}\in
S}\sum_{a\in A(s^{{}^{\prime}})}\sum_{b\in
B(s^{{}^{\prime}})}x_{ns^{{}^{\prime}}ab}^{t}.r(s^{{}^{\prime}},a,b)$. Then

$\displaystyle\phi(t,\pi_{1},\pi_{2})=\liminf_{n\to\infty}\sum_{s^{{}^{\prime}}\in
S}\sum_{a\in A(s^{{}^{\prime}})}\sum_{b\in
B(s^{{}^{\prime}})}x_{ns^{{}^{\prime}}ab}^{t}.r(s^{{}^{\prime}},a,b)$
$\displaystyle=$
$\displaystyle\liminf_{n\to\infty}[\xi_{n}^{\pi_{1}\pi_{2}}(t)].\bar{r}.$
(3.14) $\displaystyle=$
$\displaystyle\liminf_{n\to\infty}\Theta(\xi_{n}^{\pi_{1}\pi_{2}})(t)$

Where $\bar{r}$ is the reward vector of order $zd^{2}$. Define
$\Theta(\xi^{fg}(t))=\lim_{n\to\infty}\Theta(\xi_{n}^{fg}(t))$, considering
that the limit exists for all pure stationary strategy pair $(f,g)\in
F_{1}^{sp}\times F_{2}^{sp}$.

Let $p(s^{{}^{\prime}}\mid t,f^{\ast},g^{\ast})$ be the transition probability
from the state $t$ to $s^{{}^{\prime}}$ defined for the strategy pair
$(f^{\ast},g^{\ast})$. As this is a stochastic game, we can apply the Markov
property that for any two states $x,y\in S$ and $m,n\in\mathbb{N}$,

$\displaystyle p^{n+m}(x,y)$ $\displaystyle=$ $\displaystyle P(X_{n+m}=y\mid
X_{0}=x)$ (3.15) $\displaystyle=$ $\displaystyle\sum_{z\in S}P(X_{n}=z\mid
X_{0}=x)P(X_{n+m}=y\mid X_{0}=x,X_{n}=z)$ $\displaystyle=$
$\displaystyle\sum_{z\in S}P^{n}(x,z)P(X_{n+m}=y\mid X_{0}=x,X_{n}=z)$
$\displaystyle=$ $\displaystyle p^{n}(x,z)p^{m}(z,y)$

where $p^{m}(z,y)$ is the $m$-th step transition probability from the state
$z$ to $y$. Now using the above property and using the definition of
$\xi^{f_{t}g_{t}}(t)$ we have

$\displaystyle\Theta(\xi^{f_{t}g_{t}}(t))$ $\displaystyle=$
$\displaystyle\Theta(\sum_{s^{{}^{\prime}}\in S}p(s^{{}^{\prime}}\mid
t,f^{\ast},g^{\ast})\xi^{f_{t}g_{t}}(s^{{}^{\prime}}))$ (3.16)
$\displaystyle=$ $\displaystyle\sum_{s^{{}^{\prime}}\in
S}p(s^{{}^{\prime}}\mid
t,f^{\ast},g^{\ast})\Theta(\xi^{f_{t}g_{t}}(s^{{}^{\prime}}))\text{[as
$\Theta$ is a continuous function]}$

Now, as ($f_{t},g_{t}$) is an optimal pure stationary strategy pair for the
players when the initial state is t, we can write $(3.16)$ as

$\Theta(\xi^{f_{t}g_{t}}(t))=\sum_{s^{{}^{\prime}}\in S}p(s^{{}^{\prime}}\mid
t,f^{\ast},g^{\ast})\Theta(\xi^{f_{s^{{}^{\prime}}}g_{s^{{}^{\prime}}}}(s^{{}^{\prime}}))$
(3.17)

Now iterating $(3.17)$ $l$ times we get

$\displaystyle\Theta(\xi^{f_{t}g_{t}}(t))$ $\displaystyle=$
$\displaystyle\sum_{s^{{}^{\prime}}\in S}p^{l}(s^{{}^{\prime}}\mid
t,f^{\ast},g^{\ast})\Theta(\xi^{f_{s^{{}^{\prime}}}g_{s^{{}^{\prime}}}}(s^{{}^{\prime}}))$
(3.18)

If we expand the right hand side of the above expression, the right hand side
becomes:

$p^{l}(0\mid t,f^{\ast},g^{\ast})[\sum_{s\in S}\sum_{a\in A(s)}\sum_{b\in
B(s)}r_{0}(s,a,b)x_{sab}^{0}]+\cdots+p^{l}(z\mid
t,f^{\ast},g^{\ast})[\sum_{s\in S}\sum_{a\in A(s)}\\\ \sum_{b\in
B(s)}r_{z}(s,a,b)x^{z}_{sab}]$ (3.19)

Let $r^{{}^{\prime}}(s,a,b)=r_{0}(s,a,b)+r_{1}(s,a,b)+\cdots+r_{z}(s,a,b)$.
Then we can write $(3.19)$ as:

$\displaystyle\Theta(\xi^{f_{t}g_{t}}(t))$ $\displaystyle=$
$\displaystyle\sum_{s\in S}\sum_{a\in A(s)}\sum_{b\in
B(s)}r^{{}^{\prime}}(s,a,b).x_{sab}^{t}$ (3.20) $\displaystyle=$
$\displaystyle\Theta(\xi^{f^{\ast}g^{\ast}}(t))$

Thus form $(3.20)$ and $(3.13)$ we get that

$\phi(t,f^{\ast},g^{\ast})\leq\phi(t,f^{\ast},g)\forall t\in S$ and $\forall
g\in F_{2}^{sp}.$

∎

By similar manner we can show that
$\phi(t,f^{\ast},g^{\ast})\geq\phi(t,f,g^{\ast})\forall t\in S$ and $\forall
f\in F_{1}^{sp}$. Thus the pair ($f^{\ast},g^{\ast})$ is the optimal pure
stationary strategy pair of the players in the PISG $\Gamma$. ∎

##  4 Algorithm to solve a zero-sum two person perfect information stochastic
game

Let $\Gamma$ be a zero-sum two person perfect information stochastic game. We
consider the following policy-improvement algorithm to compute optimal
stationary strategy of the players. This is a best response algorithm, in
which each player looks for his own Blackwell optimal strategy. The algorithm
is stated below:
Step 1: Choose a random pure strategy for player-II $g_{0}$ and set $k=0$.
Step 2: Find the Blackwell optimal strategy $f_{k}$ for player-I in the MDP
$\Gamma(g_{k})$.
Step 3: if $g_{k}$ is blackwell optimal strategy for player-II in
$\Gamma(f_{k})$, set $(f^{\ast},g^{\ast})=(f_{k},g_{k})$ and stop.
Step 4: else find the blackwell optimal strategy $g_{k+1}$ for player-II in
the MDP $\Gamma(f_{k})$, set $k=k+1$ and go to step $2$. It is obvious that
this is a finite step algorithm and it terminates in finite time by the
conjecture $8.1$ of Raghavan and Syed (2002) [11]. The process of finding a
Blackwell optimal strategy for an undiscounred MDP was proposed by Hordijk et
al.(1985) [6]. It consists of a linear programming problem with several
parameters as given below:

$\max\sum_{s=1}^{z}\sum_{a\in A(s)}r(s,a)w_{sa}$

subject to:

$\sum_{s=1}^{z}\sum_{a\in
A(s)}(\delta(s,s^{{}^{\prime}})-q(s^{{}^{\prime}}\mid s,a))w_{sa}=0$,
$s^{{}^{\prime}}\in S$

$\sum_{a\in A(s)}w_{sa}+\sum_{s=1}^{z}\sum_{a\in
A(s)}(\delta(s,s^{{}^{\prime}})-q(s^{{}^{\prime}}\mid s,a))y_{sa}=\beta_{s}$,
$s^{{}^{\prime}}\in S$

$w_{sa}\geq 0$

where $\beta_{s}>0$ are given numbers for each $s\in S$, such that $\sum_{s\in
S}\beta_{s}=1$. The Blackwell optimal pure stationary strategy is computed as:

$f^{\ast}(s)=\frac{w^{\ast}_{sa}}{\sum_{a\in A(s)w^{\ast}_{sa}}}$

where $w^{\ast}_{sa}$ is the optimal solutionof the above LP. By Hordijk et
al.[6], this pure stationary strategy is average optimal as well. We elaborate
the above algorithm by following examples:

##  5 Numerical examples

Example 1: Consider a PISG $\Gamma$ with three states $S=\\{1,2,3\\}$,
$A(1)=\\{1,2\\}=A(2)$, $A(3)=\\{1\\}$, $B(1)=\\{1\\}=B(2)$ AND
$B(3)=\\{1,2,3\\}$. In this example player-I is a dummy player here for the
state $3$ and player-II is dummy for states $1$ and $2$. Rewards and
transition probabilities for the players are given below

State-1:  5 ($\frac{1}{2}$, $\frac{1}{2}$, 0) 7 (0,1,0) State-2:  1
($\frac{1}{3}$,0, $\frac{2}{3}$) 0.5 (0,0,1) State-3:  3 (0,
$\frac{1}{2}$,$\frac{1}{2}$) 4 (1,0,0) 2 ($\frac{1}{2}$, $\frac{1}{4}$,
$\frac{1}{4}$)
where a cell  r ($q_{1}$, $q_{2}$, $q_{3}$) represents that $r$ is the
immediate reward function and $(q_{1},q_{2},q_{3})$ are the transition
probabilities that the next states are $1$, $2$ and $3$ respectively if this
cell is chosen at present state. The pure strategies for player-I are:
$f_{0}=\\{(1,0),(1,0),1\\}$, $f_{1}=\\{(1,0),(0,1),1\\}$,
$f_{2}=\\{(0,1),(1,0)\\}$,$f_{3}=\\{(0,1),(0,1)\\}$. The pure strategies of
player-II are: $g_{0}=\\{1,1,(1,0,0)\\}$, $g_{1}=\\{1,1,(0,1,0)\\}$,
$g_{2}=\\{1,1,(0,0,1)\\}$. Firstly set $k=0$ and we fix the strategy $g_{0}$
of the player-II in $\Gamma$. Thus we get a reduced MDP $\Gamma(g_{0})$ given
below:

State-1:  5 ($\frac{1}{2}$, $\frac{1}{2}$, 0) 7 (0,1,0) State-2:  1
($\frac{1}{3}$,0, $\frac{2}{3}$) 0.5 (0,0,1) State-3:  3 (0,
$\frac{1}{2}$,$\frac{1}{2}$)
Now we formulate the following linear programming problem inthe variables
$x=(x_{11},x_{12},x_{21},\\\ x_{22},x_{31})$ and
$y=(y_{11},y_{12},y_{21},y_{22},y_{31})$ to obtain player-I’s Blackwell
optimal strategy:

$\max R=5x_{11}+7x_{12}+x_{21}+0.5x_{22}+3x_{31}$

subject to

$\displaystyle 3x_{11}+6x_{12}-2_{x21}=0$ (5.1)
$\displaystyle-3x_{11}-6x_{12}+6x_{21}+6x_{22}-3x_{31}=0$ (5.2)
$\displaystyle-8x_{21}-12x_{22}+6x_{31}=0$ (5.3) $\displaystyle
6x_{11}+6x_{12}+3y_{11}+6y_{12}-2y_{21}=6\beta_{1}$ (5.4) $\displaystyle
2x_{21}+2x_{22}-y_{11}-2y_{12}+2y_{21}+2.y_{22}-y_{31}=2\beta_{2}$ (5.5)
$\displaystyle 12x_{31}-8y_{21}-12y_{22}+6y_{31}=12\beta_{3}$ (5.6)
$\displaystyle x,y\geq 0.$ (5.7)

We fix $\beta_{1}=\beta_{2}=\beta_{3}=\frac{1}{3}$. The solution of the above
linear programming problem bt dual-simplex method is given below:
$\max R=2.778$, $x=(0.222,0,0.333,0,0.444)$, $y=(0,0.111,0,0.111,0)$.
Now by the method to compute optimal pure stationary strategy described in
section $4$, we get that $f_{0}=\\{(1,0),(0,1),1\\}$ is the optimal pure
stationary strategy for player-I in $\Gamma(g_{0})$. Now we fix this strategy
for player-I. Thus we get a resultant MDP as follows:
State-1:  5 ($\frac{1}{2}$, $\frac{1}{2}$, 0) State-2:  1 ($\frac{1}{3}$,0,
$\frac{2}{3}$) State-3:  3 (0, $\frac{1}{2}$,$\frac{1}{2}$) 4 (1,0,0) 2
($\frac{1}{2}$, $\frac{1}{4}$, $\frac{1}{4}$) We formulate the linear
programming problem of the above MDP for the variables
$x=(x_{11},x_{21},x_{31},x_{32},x_{33})$ and
$y=(y_{11},y_{21},y_{31},y_{32},y_{33})$ as follows:

$\min R^{{}^{\prime}}=5x_{11}+x_{21}+3x_{31}+4x_{32}+2x_{33}$

subject to

$\displaystyle 3x_{11}-2x_{21}-6x_{32}-3x_{33}=0$ (5.8)
$\displaystyle-2x_{11}+4x_{21}-2x_{31}-x_{33}=0$ (5.9)
$\displaystyle-8x_{21}+6x_{31}+12x_{32}+9x_{33}=0$ (5.10) $\displaystyle
6x_{11}+3y_{11}-2y_{21}-6y_{32}-3y_{33}=6\beta_{1}$ (5.11) $\displaystyle
4x_{21}-2y_{11}+4y_{21}-2y_{31}-y_{33}=4\beta_{2}$ (5.12) $\displaystyle
12x_{31}+12x_{32}+12x_{33}-8y_{21}+6y_{31}+12y_{32}+9y_{33}=12\beta_{3}$
(5.13) $\displaystyle x,y\geq 0.$ (5.14)

The solution of the above LP by dual-simplex method is given below:
$\min R^{{}^{\prime}}=2.778$, $x=(0.222,0.333,0.444,0,0)$,
$y=(0.333,0.1667,0,0,0)$. So by the same method described in section $4$, we
compute the optimal pure stationary strategy for player-II as:
$g_{0}=\\{1,1,(1,0,0)\\}$. Thus the algorithm stops in this step and we get
the optimal pure (limiting average) stationary strategy
$(F^{\ast},G^{\ast})=(f_{0},g_{0})$.

Example 2: Consider a PISG $\Gamma$ with four states $S=\\{1,2,3,4\\}$,
$A(1)=\\{1,2\\}=A(2)$, $A(3)=\\{1\\}$, $B(1)=\\{1\\}=B(2)$ AND
$B(3)=\\{1,2,3\\}$. In this example player-I is a dummy player here for the
state $3$ and player-II is dummy for states $1$ and $2$. Rewards and
transition probabilities for the players are given below

State-1:  2 ($\frac{1}{2}$, $\frac{1}{2}$, 0,0) 3 (0,1,0,0) State-2:  1
($\frac{1}{3}$,0, $\frac{2}{3}$,0) 0.5 (0,0,1,0) State-3:  $5$ $0$
$(0,0,\frac{1}{2},\frac{1}{2})$ $(0,0,1,0)$

State 4:  $11$ $12$ $(\frac{1}{2},0,\frac{1}{2},0)$ $(1,0,0,0)$

where a cell  $r$ ($q_{1}$, $q_{2}$, $q_{3}$, $q_{4}$) represents that $r$ is
the immediate reward function and $(q_{1},q_{2},q_{3},q_{4})$ are the
transition probabilities that the next states are $1$, $2$, $3$ and $4$
respectively if this cell is chosen at present state. The pure strategies for
player-I are: $f_{0}=\\{(1,0),(1,0),1,1\\}$, $f_{1}=\\{(1,0),(0,1),1,1\\}$,
$f_{2}=\\{(0,1),(1,0),1,1\\}$,$f_{3}=\\{(0,1),(0,1),1,1\\}$. The pure
strategies of player-II are: $g_{0}=\\{1,1,(1,0),(1,0)\\}$,
$g_{1}=\\{1,1,(1,0),(0,1)\\}$, $g_{2}=\\{1,1,(0,1),(1,0)\\}$ and
$g_{3}=\\{1,1,(0,1),(0,1)\\}$. Firstly set $k=0$ and we fix the strategy
$g_{0}$ of the player-II in $\Gamma$. Thus we get a reduced MDP
$\Gamma(g_{0})$ given below:

State-1:  2 ($\frac{1}{2}$, $\frac{1}{2}$, 0,0) 3 (0,1,0,0) State-2:  1
($\frac{1}{3}$,0, $\frac{2}{3}$,0) 0.5 (0,0,1,0) State-3:  $5$
($0,0,\frac{1}{2},\frac{1}{2}$)

State 4:  $11$ ($\frac{1}{2},0,\frac{1}{2},0$)

Now we formulate the following linear programming problem in the variables
$x=(x_{11},x_{12},x_{21},x_{22},x_{31},x_{32},x_{41},x_{42})$ and
$y=(y_{11},y_{12},y_{21},y_{22},y_{31},y_{32},y_{41},y_{42})$ to obtain
player-I’s Blackwell optimal strategy:

$\max R=2x_{11}+3x_{12}+0.5x_{22}+5x_{31}+11x_{41}$

subject to

$\displaystyle 3x_{11}+6x_{12}-2x_{21}-3x_{41}=0$ (5.15)
$\displaystyle-3x_{11}-6x_{12}+6x_{21}+6x_{22}=0$ (5.16)
$\displaystyle-4x_{21}-6x_{22}+3x_{31}-3x_{41}=0$ (5.17)
$\displaystyle-3x_{31}+6x_{41}=0$ (5.18) $\displaystyle
6x_{11}+6x_{12}+3y_{11}+6y_{12}-2_{y}{21}-3y_{41}=6\beta_{1}$ (5.19)
$\displaystyle 2x_{21}+2x_{22}-3y_{11}-6y_{12}+6y_{21}+6y_{22}=2\beta_{2}$
(5.20) $\displaystyle 6x_{31}--4y_{21}-6y_{22}+3y_{31}-3y_{41}=12\beta_{3}$
(5.21) $\displaystyle 2x_{41}-x_{31}+2x_{41}=2\beta_{4}$ (5.22) $\displaystyle
x,y\geq 0.$ (5.23)

We fix $\beta_{1}=\beta_{2}=\beta_{3}=\beta_{4}=\frac{1}{4}$. The solution of
the above linear programming problem bt dual-simplex method is given below:
$\max R=5.6875$, $x=(0,0.1250,0,0.1250,0.5000,0.2500)$,
$y=(0,0.1250,0,0.2500,0,0)$.
Now by the method to compute optimal pure stationary strategy described in
section $4$, we get that $f_{0}=\\{(0,1),(0,1),1,1\\}$ is the optimal pure
stationary strategy for player-I in $\Gamma(g_{0})$. Now we fix this strategy
for player-I. Thus we get a resultant MDP as follows:
State-1:  3 $(0,1,0,0)$ State-2:  0.5 (0,0,1,0) State-3:  $5$ $0$
$(0,0,\frac{1}{2},\frac{1}{2})$ $(0,0,1,0)$

State 4:  $11$ $12$ $(\frac{1}{2},0,\frac{1}{2},0)$ $(1,0,0,0)$

Now we formulate the following linear programming problem in the variables
$x=(x_{11},x_{12},x_{21},x_{22},\\\ x_{31},x_{32},x_{41},x_{42})$ and
$y=(y_{11},y_{12},y_{21},y_{22},y_{31},y_{32},y_{41},y_{42})$ to obtain
player-I’s Blackwell optimal strategy:

$\min R=2x_{11}+3x_{21}+0.5x_{31}+5x_{41}+11x_{42}$

subject to

$\displaystyle 2x_{11}-x_{41}-2x_{42}=0$ (5.24) $\displaystyle-
x_{11}+x_{21}=0$ (5.25) $\displaystyle-2x_{21}+x_{31}-x_{41}=0$ (5.26)
$\displaystyle-x_{31}+2x_{41}+2x_{42}=0$ (5.27) $\displaystyle
2x_{11}+2y_{11}-y_{41}-2y_{42}=2\beta_{1}$ (5.28) $\displaystyle
x_{21}-y_{11}+y_{21}=\beta_{2}$ (5.29) $\displaystyle
2x_{31}+2x_{32}-2y_{21}+y_{31}-y_{41}=2\beta_{3}$ (5.30) $\displaystyle
x_{41}+x_{42}-y_{31}+2y_{41}+y_{42}=2\beta_{4}$ (5.31) $\displaystyle x,y\geq
0.$ (5.32)

The solution of the above LP by dual-simplex method is given below:
$\min R^{{}^{\prime}}=2.778$, $x=(0.1250,0.1250,0.5,0,0.25,0)$,
$y=(0.1250,0.2500,0,0,0,0)$. So by the same method described in section $4$,
we compute the optimal pure stationary strategy for player-II as:
$g_{0}=\\{1,1,(1,0),(0,1))\\}$. Thus the algorithm stops in this step and we
get the optimal pure (limiting average) stationary strategy
$(F^{\ast},G^{\ast})=(f_{0},g_{0})$.

## References

  * [1] Cyrus Derman.  On sequential control processes.  The Annals of Mathematical Statistics, 35(1):341–349, 1964.
  * [2] Melvin Dresher, Lloyd S Shapley, and Albert William Tucker.  Advances in Game Theory.(AM-52), Volume 52, volume 52.  Princeton University Press, 2016.
  * [3] Drew Fudenberg and Eric Maskin.  The folk theorem in repeated games with discounting or with incomplete information.  In A long-run collaboration on long-run games, pages 209–230. World Scientific, 2009.
  * [4] Dean Gillette.  Stochastic games with zero stop probabilities.  Contributions to the Theory of Games, 3:179–187, 1957.
  * [5] Alan J Hoffman and Richard M Karp.  On nonterminating stochastic games.  Management Science, 12(5):359–370, 1966.
  * [6] Arie Hordijk, Rommert Dekker, and Lodewijk Cornelis Maria Kallenberg.  Sensitivity-analysis in discounted markovian decision problems.  Operations-Research-Spektrum, 7(3):143–151, 1985.
  * [7] John G Kemeny and J Laurie Snell.  Finite Markov chains: with a new appendix” Generalization of a fundamental matrix”.  Springer, 1983.
  * [8] Thomas M Liggett and Steven A Lippman.  Stochastic games with perfect information and time average payoff.  Siam Review, 11(4):604–607, 1969.
  * [9] Lynn H Loomis.  On a theorem of von neumann.  Proceedings of the National Academy of Sciences, 32(8):213–215, 1946.
  * [10] J von Neumann.  A model of general economic equilibrium.  In Readings in the Theory of Growth, pages 1–9. Springer, 1971\.
  * [11] TES Raghavan and Zamir Syed.  A policy-improvement type algorithm for solving zero-sum two-person stochastic games of perfect information.  Mathematical Programming, 95(3):513–532, 2003.
  * [12] Lloyd S Shapley.  Stochastic games.  Proceedings of the national academy of sciences, 39(10):1095–1100, 1953.