Department of Computer Science, James Parsons Building, Byrom Street,
Liverpool John Moores University, Liverpool, L3 3AF,
UKp.c.bell@ljmu.ac.ukhttps://orcid.org/0000-0003-2620-635XDepartment of
Computer Science, University of Oxford, Wolfson Building, Parks Road, Oxford,
OX1 3QD, UKpavel.semukhin@cs.ox.ac.ukhttps://orcid.org/0000-0002-7547-6391
Paul C. Bell and Pavel Semukhin <ccs2012> <concept>
<concept_id>10003752.10003766</concept_id> <concept_desc>Theory of computation
Formal languages and automata theory</concept_desc>
<concept_significance>500</concept_significance> </concept> <concept>
<concept_id>10003752.10003753.10003754</concept_id> <concept_desc>Theory of
computation Computability</concept_desc>
<concept_significance>500</concept_significance> </concept> <concept>
<concept_id>10003752.10003753.10003757</concept_id> <concept_desc>Theory of
computation Probabilistic computation</concept_desc>
<concept_significance>500</concept_significance> </concept> </ccs2012>
[500]Theory of computation Formal languages and automata theory [500]Theory of
computation Computability [500]Theory of computation Probabilistic computation

# Decision Questions for Probabilistic Automata on Small Alphabets

Paul C. Bell  Pavel Semukhin

###### Abstract

We study the emptiness and $\lambda$-reachability problems for unary and
binary Probabilistic Finite Automata (PFA) and characterise the complexity of
these problems in terms of the degree of ambiguity of the automaton and the
size of its alphabet. Our main result is that emptiness and
$\lambda$-reachability are solvable in EXPTIME for polynomially ambiguous
unary PFA and if, in addition, the transition matrix is over $\\{0,1\\}$, we
show they are in NP. In contrast to the Skolem-hardness of the
$\lambda$-reachability and emptiness problems for exponentially ambiguous
unary PFA, we show that these problems are NP-hard even for finitely ambiguous
unary PFA. For binary polynomially ambiguous PFA with fixed and commuting
transition matrices, we prove NP-hardness of the $\lambda$-reachability
(dimension $9$), nonstrict emptiness (dimension $37$) and strict emptiness
(dimension $40$) problems.

###### keywords:

Probabilistic finite automata; unary alphabet; emptiness problem; bounded
ambiguity

##  1 Introduction

There are many possible extensions of the fundamental notion of a
nondeterministic finite automaton. One such notion is that of a Probabilistic
Finite Automata (PFA) which was first introduced by Rabin [19]. In a PFA
$\mathcal{P}$ over a (finite) input alphabet $\Sigma$ the outgoing transitions
from a state, for each input letter of $\Sigma$, form a probability
distribution, as does the initial state vector. Thus, a word $w\in\Sigma^{*}$
is accepted with a certain probability, which we denote $\mathcal{P}(w)$.

There are a variety of interesting questions that one may ask about a PFA
$\mathcal{P}$ over an alphabet $\Sigma$. In this article we focus on two
decision questions, that of _$\lambda$ -reachability_ and _emptiness_. The
$\lambda$-reachability problem is stated thus: given a probability
$\lambda\in[0,1]$, does there exist some word $w\in\Sigma^{*}$ such that
$\mathcal{P}(w)=\lambda$? In the (strict) emptiness problem, we ask if there
exists _any_ word $w\in\Sigma^{*}$ such that $P(w)\geq\lambda$ (resp.
$P(w)>\lambda$). We also mention the related _cutpoint isolation_ problem — to
determine if for each $\epsilon>0$, there exists a word $w\in\Sigma$ such that
$|\mathcal{P}(w)-\lambda|<\epsilon$.

In general, the emptiness problem is undecidable for PFA [18], even over a
binary alphabet when the automaton has $25$ states [12]. The cutpoint
isolation problem is undecidable [4] even for PFA with $420$ states over a
binary alphabet [5]. The problem is especially interesting given the seminal
result of Rabin that if a cutpoint $\lambda$ is isolated, then the cutpoint
language associated with $\lambda$ is necessarily regular [19].

We may ask which restrictions of PFA may lead to decidability of the previous
problems. In this paper we are interested in PFA of _bounded ambiguity_ ,
where the ambiguity of a word denotes the number of accepting runs of that
word in the PFA. A PFA $\mathcal{P}$ is $f$-ambiguous, for a function
$f:\mathbb{N}\to\mathbb{N}$, if every word of length $n$ has at most $f(n)$
accepting runs. A run is _accepting_ if the probability of that run ending in
a final state is strictly positive. The degree of ambiguity is thus a property
of the NFA underlying a PFA (i.e. the NFA produced by setting all nonzero
transition probabilities to $1$). We may consider the notions of finite,
polynomial or exponential ambiguity of $\mathcal{P}$ based on whether $f$ is
bounded by a constant, is a polynomial or else is exponential, respectively.
Characterisations of the degree of ambiguity of NFA are given by Weber and
Seidel [23].

The authors of [8] show that emptiness for PFA remains undecidable even for
polynomially ambiguous automata (quadratic ambiguity), show PSPACE-hardness
results for finitely ambiguous PFA and that emptiness is in NP for the class
of $k$-ambiguous PFA for every $k>0$. The emptiness problem for PFA was later
shown to be undecidable for linearly ambiguous automata [7].

Another restriction is to constrain input words of the PFA to come from a
given language $\mathcal{L}$. If $\mathcal{L}$ is a _letter-bounded_ language,
then the emptiness and $\lambda$-reachability problems remain undecidable for
polynomially ambiguous PFA, even when all transition matrices commute [2]. In
contrast, the _cutpoint-isolation_ problem is decidable even for exponentially
ambiguous PFA when inputs are constrained to come from a given letter-bounded
context-free language, although it is NP-hard for $3$-state PFA on letter-
bounded inputs [3].

Our main results are as follows. We show that the $\lambda$-reachability and
emptiness problems for probabilistic finite automata are:

  * •

In EXPTIME for the class of polynomially ambiguous unary PFA and are NP-
complete if, in addition, the transition matrix is over $\\{0,1\\}$ [Theorem
3.3 and Corollary 1].

  * •

NP-hard for polynomially ambiguous PFA over a binary alphabet with _fixed_ and
_commuting_ transition matrices of dimension $40$ (strict emptiness problem),
$37$ (nonstrict emptiness problem) and $9$ ($\lambda$-reachability problem)
[Theorem 4.1].

We also show NP-hardness for the class of finitely ambiguous unary PFA with
$\\{0,1\\}$ transition matrix [Theorem 3.15]. Our hardness results rely on the
NP-hardness of solving binary quadratic equations and the universality problem
for unary regular expressions. An interesting question, that is left open, is
to find out the exact computational complexity of the above problems in the
case of polynomially ambiguous unary PFA, i.e. to close the gap between the
EXPTIME upper bound and NP lower bound.

##  2 Probabilistic Finite Automata and Notation

We denote by $\mathbb{Q}^{n\times n}$ the set of all $n\times n$ matrices over
$\mathbb{Q}$. Given two column vectors $u\in\mathbb{Q}^{n}$ and
$v\in\mathbb{Q}^{m}$, we denote by $[u|v]$ the column vector
$(u_{1},\ldots,u_{n},v_{1},\ldots,v_{m})^{T}\in\mathbb{Q}^{n+m}$. For a
sequence of vectors $u_{1},u_{2},\ldots,u_{k}$, we write
$[u_{1}|u_{2}|\cdots|u_{k}]$ for the column vector which stacks the vectors on
top of each other.

Given $A=(a_{ij})\in{\mathbb{Q}}^{m\times m}$ and $B\in{\mathbb{Q}}^{n\times
n},$ we define the direct sum $A\oplus B$ and Kronecker product $A\otimes B$
of $A$ and $B$ by:

$A\oplus B=\left[\begin{array}[]{@{}c|l@{}}A&\mathbf{0}_{m,n}\\\
\hline\cr\mathbf{0}_{n,m}&B\end{array}\right],\quad A\otimes
B=\left[\begin{array}[]{cccc}a_{11}B&a_{12}B&\cdots&a_{1m}B\\\
a_{21}B&a_{22}B&\cdots&a_{2m}B\\\ \vdots&\vdots&&\vdots\\\
a_{m1}B&a_{m2}B&\cdots&a_{mm}B\\\ \end{array}\right],$

where $\mathbf{0}_{i,j}$ denotes the zero matrix of dimension $i\times j$.
Note that neither $\oplus$ nor $\otimes$ are commutative in general. The
following useful properties of $\oplus$ and $\otimes$ are well known.

######  Lemma 2.1.

Let $A,B,C,D\in\mathbb{Q}^{n\times n}$. Then we have:

  * •

Associativity: $(A\otimes B)\otimes C=A\otimes(B\otimes C)$ and $(A\oplus
B)\oplus C=A\oplus(B\oplus C)$, thus $A\otimes B\otimes C$ and $A\oplus
B\oplus C$ are unambiguous.

  * •

Mixed product properties: $(A\otimes B)(C\otimes D)=(AC\otimes BD)$ and
$(A\oplus B)(C\oplus D)=(AC\oplus BD)$.

  * •

If $A$ and $B$ are stochastic matrices, then so are $A\oplus B$ and $A\otimes
B$.

  * •

If $A,B\in\mathbb{Q}^{n\times n}$ are both upper-triangular, then so are
$A\oplus B$ and $A\otimes B$.

See [13] for proofs of the first three properties of Lemma 2.1. The fourth
property follows directly from the definition of the Kronecker sum and product
and is not difficult to prove.

A Probabilistic Finite Automaton (PFA) $\mathcal{P}$ with $n$ states over an
alphabet $\Sigma$ is defined as $\mathcal{P}=(u,\\{M_{a}|a\in\Sigma\\},v)$
where $u\in\mathbb{Q}^{n}$ is the _initial probability distribution_ ;
$v\in\\{0,1\\}^{n}$ is the _final state vector_ and each
$M_{a}\in\mathbb{Q}^{n\times n}$ is a (row) stochastic matrix. We will
primarily be interested in _unary_ and _binary_ PFA, for which $|\Sigma|=1$
and $|\Sigma|=2$ respectively. For a word $w=a_{1}a_{2}\cdots
a_{k}\in\Sigma^{*}$, we define the acceptance probability
${\mathcal{P}}(w):\Sigma^{*}\to\mathbb{Q}$ of ${\mathcal{P}}$ as:

${\mathcal{P}}(w)=u^{T}M_{a_{1}}M_{a_{2}}\cdots M_{a_{k}}v\in[0,1],$

which denotes the acceptance probability of $w$.111Some authors interchange
the order of ${u}$ and ${v}$ and use column stochastic matrices, although the
two definitions are trivially equivalent.

For a given cutpoint $\lambda\in[0,1]$, we define the following languages:
$L_{\geq\lambda}({\mathcal{P}})=\\{w\in\Sigma^{*}\ |\
{\mathcal{P}}(w)\geq\lambda\\}$, a nonstrict cutpoint language, and
$L_{>\lambda}({\mathcal{P}})=\\{w\in\Sigma^{*}\ |\
{\mathcal{P}}(w)>\lambda\\}$, a strict cutpoint language. The (strict)
emptiness problem for a cutpoint language is to determine if
$L_{\geq\lambda}({\mathcal{P}})=\emptyset$ (resp.
$L_{>\lambda}({\mathcal{P}})=\emptyset$). We are also interested in the
_$\lambda$ -reachability_ problem, for which we ask if there exists a word
$w\in\Sigma^{*}$ such that $\mathcal{P}(w)=\lambda$.

###  2.1 PFA Ambiguity

The degree of ambiguity of a finite automaton is a structural parameter,
roughly indicating the number of accepting runs for a given input word. See
[23] for a thorough discussion of ambiguity for nondeterministic automata and
[2, 3, 7, 8] for connections to PFA.

Let $w\in\Sigma^{*}$ be an input word of an NFA
$\mathcal{N}=(Q,\Sigma,\delta,Q_{I},Q_{F})$, with $Q$ the set of states,
$\Sigma$ the input alphabet, $\delta\subset Q\times\Sigma\times Q$ the
transition function, $Q_{I}$ the set of initial states and $Q_{F}$ the set of
final states. For each $(p,w,q)\in Q\times\Sigma^{*}\times Q$, define
$\textrm{da}_{\mathcal{N}}(p,w,q)$ as the number of paths for $w$ in
$\mathcal{N}$ leading from state $p$ to $q$. The _degree of ambiguity_ of $w$
in $\mathcal{N}$, denoted $\textrm{da}_{\mathcal{N}}(w)$, is defined as the
number of all _accepting paths_ for $w$ (starting from an initial and ending
in a final state). The _degree of ambiguity_ of $\mathcal{N}$, denoted
$\textrm{da}(\mathcal{N})$, is the supremum of the set
$\\{\textrm{da}_{\mathcal{N}}(w)\ |\ w\in\Sigma^{*}\\}$. $\mathcal{N}$ is
called infinitely ambiguous if $\textrm{da}(\mathcal{N})=\infty$, finitely
ambiguous if $\textrm{da}(\mathcal{N})<\infty$, and unambiguous if
$\textrm{da}(\mathcal{N})\leq 1$. The _degree of growth_ of the ambiguity of
$\mathcal{N}$, denoted $\textrm{deg}(\mathcal{N})$, is defined as the minimum
degree of a univariate polynomial $h$ with positive integral coefficients such
that for all $w\in\Sigma^{*}$, $\textrm{da}_{\mathcal{N}}(w)\leq h(|w|)$ (if
such a polynomial exists, in which case $\mathcal{N}$ is called polynomially
ambiguous, otherwise the degree of growth is infinite and $\mathcal{N}$ is
called exponentially ambiguous).

The above notions relate to NFA. We may derive an analogous notion of
ambiguity for PFA by considering an embedding of a PFA ${\mathcal{P}}$ to an
NFA $\mathcal{N}$ in such a way that for each letter $a\in\Sigma$, if the
probability of transitioning from a state $i$ to state $j$ is nonzero under
${\mathcal{P}}$, then there is an edge from state $i$ to $j$ under
$\mathcal{N}$ for letter $a$. The initial states of $\mathcal{N}$ are those of
${\mathcal{P}}$ having nonzero initial probability and the final states of
$\mathcal{N}$ and ${\mathcal{P}}$ coincide. We then say that ${\mathcal{P}}$
is finitely/polynomially/exponentially ambiguous if $\mathcal{N}$ is
(respectively).

A state $q\in Q$ in an NFA (resp. PFA) is called _useful_ if there exists an
accepting path which visits $q$ (resp. an accepting path of nonzero
probability which visits $q$). We can characterise whether an NFA
$\mathcal{A}$ (and thus a PFA by the above embedding) has finite, polynomial
or exponential ambiguity using the following properties:

EDA \- There is a useful state $q\in Q$ such that, for some word
$v\in\Sigma^{*}$, $da_{\mathcal{A}}(q,v,q)\geq 2$.

$\text{IDA}_{d}$ \- There are useful states $r_{1},s_{1},\ldots,r_{d},s_{d}\in
Q$ and words $v_{1},u_{2},v_{2},\ldots,u_{d},v_{d}\in\Sigma^{*}$ such that for
all $1\leq i\leq d$, $r_{i}$ and $s_{i}$ are distinct and
$(r_{i},v_{i},r_{i}),(r_{i},v_{i},s_{i}),(s_{i},v_{i},s_{i})\in\delta$ and for
all $2\leq i\leq d$, $(s_{i-1},u_{i},r_{i})\in\delta$.

######  Theorem 2.2 ([14, 20, 23]).

An NFA (or PFA) $\mathcal{A}$ having the EDA property is equivalent to it
being exponentially ambiguous. For any $d\in\mathbb{N}$, an NFA (or PFA)
$\mathcal{A}$ having property $\text{IDA}_{d}$ is equivalent to
$\textrm{deg}(\mathcal{A})\geq d$.

Clearly, if $\mathcal{N}$ agrees with $\text{IDA}_{d}$ for some $d>0$, then it
also agrees with $\text{IDA}_{1},\ldots,\text{IDA}_{d-1}$. An NFA (or PFA) is
thus finitely ambiguous if it does not possess property $\text{IDA}_{1}$.

##  3 Unary PFA

Our main focus is on unary automata. We begin by giving a simple folklore
proof that the $\lambda$-reachability and emptiness problems are as
computationally difficult as the famous Skolem problem, which is only know to
be decidable for instances of depth $4$ [22]. See also [1] for connections to
reachability problems for Markov chains.

######  Theorem 3.1.

The $\lambda$-reachability and emptiness problems for unary exponentially
ambiguous Probabilistic Finite Automata are Skolem-hard.

######  Proof 3.2.

_(Folklore)._ The $\lambda$-reachability problem for unary exponentially
ambiguous PFA can be shown Skolem-hard based on the well known matrix
formulation of Skolem’s problem [11] and Turakainen’s technique showing the
equivalence of (strict) cutpoint language acceptance of generalised automata
and exponentially ambiguous probabilistic automata [21].

The emptiness problem can be shown Skolem-hard by encoding the positivity
problem which is known to be Skolem-hard, see [17] for example.

We now move to prove our main result, specifically that the emptiness and
$\lambda$-reachability problems for polynomially ambiguous unary probabilistic
finite automata are in EXPTIME. Note again that without the restriction of
polynomial ambiguity the problem is Skolem-hard by Theorem 3.1 and thus not
even known to be decidable.

######  Theorem 3.3.

The $\lambda$-reachability and (strict) emptiness problems for unary
polynomially ambiguous Probabilistic Finite Automata are decidable in EXPTIME.

In order to establish Theorem 3.3, we need to prove a series of lemmas.

The next lemma states that we may consider a unary polynomially ambiguous PFA
whose transition matrix is upper-triangular. This will prove useful since in
that case the eigenvalues of the transition matrix are rational nonnegative.
In general, a polynomially ambiguous unary PFA may have a transition matrix
with complex eigenvalues as illustrated here: $A=\begin{pmatrix}0&1&0\\\
0&0&1\\\ 1&0&0\end{pmatrix}$ whose eigenvalues are
$\left\\{1,-\frac{1}{2}\pm{\bf i}{\frac{\sqrt{3}}{2}}\right\\}$.

######  Lemma 3.4.

Let $\mathcal{P}=(u,A,v)$ be a polynomially ambiguous unary Probabilistic
Finite Automaton with acceptance function $\mathcal{P}(a^{k})=u^{T}A^{k}v$.
Then we can compute in EXPTIME a set of $d$ polynomially ambiguous unary PFAs
$\\{\mathcal{P}_{s}=(u_{s},U,v^{\prime})\ |\ 0\leq s\leq d-1\\}$ such that $U$
is rational upper-triangular and
$\mathcal{P}(a^{k})=\mathcal{P}_{s}(a^{r})=u_{s}^{T}U^{r}v^{\prime}$, where
$k=rd+s$ with $0\leq s\leq d-1$.

######  Proof 3.5.

We will identify $\mathcal{P}$ and its underlying graph in which an edge
$(p,q)$ exists iff $A_{p,q}\neq 0$. Two states $p,q$ of a PFA are said to be
_connected_ if there exists a path from $p$ to $q$ and from $q$ to $p$. We
partition the set of states into Strongly Connected Components (SCC) denoted
$S_{1},S_{2},\ldots,S_{\ell}$ so that for any SCC $S_{j}$, either $|S_{j}|=1$,
or else any two states in $S_{j}$ are connected. These SCCs can be computed in
linear time.

A polynomially ambiguous PFA does not have the EDA property (see Sec. 2.1).
This implies that every $S_{j}$, with $|S_{j}|>1$, consists of a single
directed cycle, possibly with transitions to other SCCs. To see this, suppose
there are two different directed cycles inside $S_{j}$ of lengths $m$ and $n$
and a common vertex $p$. Then one can construct two different paths of length
$mn$ from $p$ to $p$ by going $m$ times along the first cycle and $n$ time
along the second cycle, respectively, contradicting the assumption that
${\mathcal{P}}$ does not have the EDA property.

Note that if there exists a path from a state $p\in S_{j_{1}}$ to some $q\in
S_{j_{2}}$, then there does not exist any path from any state in $S_{j_{2}}$
to a state in $S_{j_{1}}$, otherwise $S_{j_{1}}$ and $S_{j_{2}}$ would merge
to a single SCC (since all vertices are then connected). This implies that the
connected components $S_{1},S_{2},\ldots,S_{\ell}$ can be reordered in such a
way that there are no transitions from $S_{j}$ to $S_{i}$ for $i<j$. Hence
there exists a permutation matrix $P$ such that the following matrix is
stochastic block upper-triangular:

$B=PAP^{-1}=\begin{pmatrix}B_{1}&*&\cdots&*\\\ 0&B_{2}&\ddots&*\\\
\vdots&\ddots&\ddots&\vdots\\\ 0&0&\cdots&B_{\ell}\end{pmatrix},$

such that each $B_{j}\in{\mathbb{Q}}^{d_{j}\times d_{j}}$, where $d_{j}$ is
the size of $S_{j}$, and $B_{j}\preceq P_{j}$, where
$P_{j}\in{\mathbb{N}}^{d_{j}\times d_{j}}$ is a permutation matrix, and the
entries $*$ are arbitrary. Here $M\preceq N$ means that $M$ is entrywise less
than $N$, i.e. $M_{i,j}\leq N_{i,j}$.

Let $d=\textrm{lcm}\\{d_{j}\ |\ 1\leq j\leq\ell\\}$ (in fact, we can simply
take $d=\prod_{j=1}^{\ell}d_{j}$). We then see that:

$U:=B^{d}=PA^{d}P^{-1}=\begin{pmatrix}B_{1}^{d}&*&\cdots&*\\\
0&B_{2}^{d}&\ddots&*\\\ \vdots&\ddots&\ddots&\vdots\\\
0&0&\cdots&B_{\ell}^{d}\end{pmatrix}.$

Note that each $B_{j}^{d}\preceq P_{j}^{d}=I_{j}$, where
$I_{j}\in{\mathbb{N}}^{d_{j}\times d_{j}}$ is the identity matrix, and the
entries $*$ are arbitrary. Therefore, each $B_{j}^{d}$ is diagonal, and so $U$
is clearly upper-triangular.

We then define $\mathcal{P}_{s}=(u_{s},U,v^{\prime})$, for $0\leq s\leq d-1$,
with $u_{s}^{T}=u^{T}A^{s}P^{-1}$ and $v^{\prime}=Pv$ noting that $Pv$ is a
binary vector as required of a final state vector. We now see that:

$\mathcal{P}(a^{k})=u^{T}A^{k}v=u^{T}A^{s}A^{rd}v=u^{T}A^{s}P^{-1}(PA^{rd}P^{-1})Pv=u_{s}^{T}U^{r}v^{\prime}=\mathcal{P}_{s}(a^{r})$

for $k=rd+s$ with $0\leq s\leq d-1$ as required. Here we used the identity
$U^{r}=PA^{rd}P^{-1}$.

Finally, note that $d$ can be exponential in the number of states of
$\mathcal{P}$, which in turn is bounded by the input size. Hence computing $U$
and all $u_{s}$, for $0\leq s\leq d-1$, takes exponential time.

The next lemma gives us an efficient method to compute an explicit formula for
the acceptance probability function of a unary PFA with upper-triangular
transition matrix.

######  Lemma 3.6.

Let $\mathcal{P}=(u,A,v)$ be a unary probabilistic finite automaton such that
$A$ is rational upper-triangular, and let
$\lambda_{0}=1>\lambda_{1}>\cdots>\lambda_{m}\geq 0$ be distinct eigenvalues
of $A$. Then there exist a constant $c\in{\mathbb{Q}}$ and univariate
polynomials $p_{1},\ldots,p_{m}$ over ${\mathbb{Q}}$, all of which can be
computed in polynomial time, such that

$\mathcal{P}(a^{k})=c+\sum_{i=1}^{m}p_{i}(k)\lambda_{i}^{k}.$

######  Proof 3.7.

First, we write $A$ in _Jordan normal form_ $A=S^{-1}JS$, where $S$ is a
nonsingular (det$(S)\neq 0$) matrix consisting of the _generalised_
eigenvectors of $A$. Note that computing $J$ and $S$ involves solving systems
of linear equations over ${\mathbb{Q}}$, and hence it can be done in
polynomial time. Computing $S^{-1}$ also requires polynomial time. Matrix $J$
has the form
$J=\bigoplus_{i=0}^{m}\bigoplus_{j=1}^{n_{i}}J_{\ell_{i,j}}(\lambda_{i})$,
where $J_{\ell_{i,j}}(\lambda_{i})$ is a $\ell_{i,j}\times\ell_{i,j}$ _Jordan
block_ and $n_{i}$ is the geometric multiplicity of $\lambda_{i}$ (hence
$\sum_{j=1}^{n_{i}}\ell_{i,j}$ is the algebraic multiplicity of
$\lambda_{i}$). Recall that a Jordan block $J_{\ell}(\lambda)$ of size
$\ell\times\ell$ that corresponds to an eigenvalue $\lambda$ has the form:

$J_{\ell}(\lambda)=\begin{pmatrix}\lambda&1&0&\cdots&0\\\
0&\lambda&1&\cdots&0\\\ 0&0&\lambda&\cdots&0\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
0&0&0&\cdots&\lambda\end{pmatrix}\in\mathbb{Q}^{\ell\times\ell}.$

Noting that $\binom{x}{y}=0$ if $y>x$, we see that

$J_{\ell}(\lambda)^{k}=\begin{pmatrix}\lambda^{k}&\binom{k}{1}\lambda^{k-1}&\binom{k}{2}\lambda^{k-2}&\cdots&\binom{k}{\ell-1}\lambda^{k-(\ell-1)}\\\
0&\lambda^{k}&\binom{k}{1}\lambda^{k-1}&\cdots&\binom{k}{\ell-2}\lambda^{k-(\ell-2)}\\\
0&0&\lambda^{k}&\cdots&\binom{k}{\ell-3}\lambda^{k-(\ell-3)}\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\ 0&0&0&\cdots&\lambda^{k}\end{pmatrix}.$
(1)

Note that the entries of $J_{\ell}(\lambda)^{k}$ have the form
$q_{i,j}(k)\lambda^{k}$, where $q_{i,j}(k)$ are polynomials over
${\mathbb{Q}}$ that can be computed in polynomial time. Namely,
$q_{i,i+p}(k)=\binom{k}{p}\lambda^{-p}$ for $0\leq p\leq\ell-i$, and
$q_{i,j}(k)=0$ for $i>j$. Note that even though $p$ appears in the exponent of
$\lambda^{-p}$ and as $p!$ in $\binom{k}{p}$, these values are still
computable in PTIME from the input data because $p$ is bounded by the
dimension of the matrix, which in turn is bounded by the input size.

Next, we note that
$J^{k}=\bigoplus_{i=0}^{m}\bigoplus_{j=1}^{n_{i}}J_{\ell_{i,j}}(\lambda_{i})^{k}$.
Hence the entries of $J^{k}$ have the form $p_{s,t}(k)\lambda_{i}^{k}$, where
$p_{s,t}(k)$ are polynomials over ${\mathbb{Q}}$. So we can write the function
$\mathcal{P}(a^{k})$ as follows:

$\mathcal{P}(a^{k})=u^{T}A^{k}v=(u^{T}S^{-1})J^{k}(Sv).$

Note that in the above equation, $u^{T}S^{-1}$ and $Sv$ are rational vectors.
It follows that

$\mathcal{P}(a^{k})=\sum_{i=0}^{m}p_{i}(k)\lambda_{i}^{k}$

for some polynomials $p_{i}(k)$ over ${\mathbb{Q}}$. In fact, these
polynomials are rational linear combinations of those $p_{s,t}(k)$ that
multiply $\lambda_{i}^{k}$ in the expression for $J^{k}$, and so they can be
computed in polynomial time.

Finally, recall that $\lambda_{0}=1$ and note that the Jordan blocks that
correspond to the dominant eigenvalues of a stochastic matrix have size
$1\times 1$ (for the proof of this fact see, e.g. [9, Theorem 6.5.3]). It
follows from (1) that the terms $\lambda_{0}^{k}$ in the formula for $J^{k}$
are multiplied by constant polynomials $p_{s,t}(k)=1$. Hence $p_{0}(k)=c$ for
some constant $c\in{\mathbb{Q}}$.

The next technical lemma is crucial in our later analysis of the running time
of the algorithms for the emptiness and $\lambda$-reachability problems
presented in Lemmas 3.10 and 3.12.

######  Lemma 3.8.

Let $D\in{\mathbb{R}}$ be such that $\ln D>2$. Then for all $x>3D\ln D$, we
have $D\ln x<x$.

######  Proof 3.9.

Our goal is to find $x_{0}>0$ such that every $x>x_{0}$ satisfies $D\ln x<x$.
First, let us make a substitution $x=Dt$, where $t>1$. Then we can rewrite
$D\ln x<x$ as follows

$\displaystyle D\ln(Dt)$ $\displaystyle<Dt,$ $\displaystyle\ln t+\ln D$
$\displaystyle<t.$

We want to find $t_{0}>1$ such that every $t>t_{0}$ satisfies $\ln t+\ln D<t$.
Let us make another substitution $t=\ln D+u\ln\ln D$, where $u>0$. Then we can
write the previous inequality as

$\displaystyle\ln(\ln D+u\ln\ln D)+\ln D$ $\displaystyle<\ln D+u\ln\ln D,$
$\displaystyle\ln\left(\ln D\left(1+u\frac{\ln\ln D}{\ln D}\right)\right)$
$\displaystyle<u\ln\ln D,$ $\displaystyle\ln\ln D+\ln\left(1+u\frac{\ln\ln
D}{\ln D}\right)$ $\displaystyle<u\ln\ln D.$ (2)

So we need to find $u_{0}>0$ such that for all $u>u_{0}$, the inequality (2)
holds. In order to do this, we can replace the left-hand side of (2) with a
larger value using $\ln\left(1+u\frac{\ln\ln R}{\ln R}\right)<u\frac{\ln\ln
R}{\ln R}$. Thus we obtain

$\displaystyle\ln\ln D+u\frac{\ln\ln D}{\ln D}$ $\displaystyle<u\ln\ln D,$
$\displaystyle 1+\frac{u}{\ln D}$ $\displaystyle<u,\qquad\ln D+u<u\ln
D,\qquad\frac{\ln D}{\ln D-1}<u.$

Recall that by our assumption $\ln D>2$. In this case, $\frac{\ln D}{\ln
D-1}<2$, and hence we can choose $u_{0}=2$. This gives us the values
$t_{0}=\ln D+u_{0}\ln\ln D=\ln D+2\ln\ln D$ and $x_{0}=Dt_{0}=D(\ln D+2\ln\ln
D)$. Since $\ln\ln D<\ln D$, we can choose $x_{0}$ to be $x_{0}=3D\ln D$.

We now proceed to the proof of our main result. We split the analysis into two
cases depending on whether or not the cutpoint $\lambda$ coincides with the
limit $\lim_{k\to\infty}\mathcal{P}(a^{k})$, which is unique by Lemma 3.6.

######  Lemma 3.10.

Let $\mathcal{P}=(u,A,v)$ be a unary probabilistic finite automaton such that
$A$ is rational upper-triangular, and let $\lambda\in[0,1]\cap{\mathbb{Q}}$ be
a cutpoint. Assuming that $\lambda\neq\lim_{k\to\infty}\mathcal{P}(a^{k})$,
the (strict) emptiness and $\lambda$-reachability problems for $\mathcal{P}$
and $\lambda$ are decidable in EXPTIME.

######  Proof 3.11.

By Lemma 3.6, we can write
$\mathcal{P}(a^{k})=c+\sum_{i=1}^{m}p_{i}(k)\lambda_{i}^{k}$, where
$1>\lambda_{1}>\cdots>\lambda_{m}$ are the eigenvalues of $A$ and $c$ and the
coefficients of $p_{i}$ are rational numbers that can be computed in
polynomial time. By assumption,
$\lim_{k\to\infty}\mathcal{P}(a^{k})=c\neq\lambda$. Let
$\epsilon=\frac{|c-\lambda|}{2}$. We now determine a natural number $k_{0}>0$
such that $\mathcal{P}(a^{k})\in(c-\epsilon,c+\epsilon)$ for all $k>k_{0}$.

Let each $p_{i}(k)$ have the form
$p_{i}(k)=a_{i,s}k^{s}+a_{i,s-1}k^{s-1}+\cdots+a_{i,0}$, where $s\leq n$ is
the size of the largest Jordan block in the Jordan normal form of $A$ (we do
not assume here that $a_{i,s}\neq 0$). Then for all $k>0$ we have

$\left|\sum_{i=1}^{m}p_{i}(k)\lambda_{i}^{k}\right|\leq\lambda_{1}^{k}\sum_{i=1}^{m}\left|p_{i}(k)\right|\leq\lambda_{1}^{k}k^{s}\sum_{i=1}^{m}\sum_{j=0}^{s}\left|a_{i,j}\right|=d\,k^{s}\lambda_{1}^{k},$

where $d=\sum_{i=1}^{m}\sum_{j=0}^{s}\left|a_{i,j}\right|\in{\mathbb{Q}}$ can
be computed in polynomial time by Lemma 3.6.

Let $k_{1}>0$ be a number to be defined later such that for all $k>k_{1}$,

$k^{s}<\left(\frac{1}{\sqrt{\lambda_{1}}}\right)^{k}=\lambda_{1}^{-\frac{k}{2}}.$

Then for all $k>k_{1}$, we have
$d\,k^{s}\lambda_{1}^{k}<d\lambda_{1}^{\frac{k}{2}}$. Thus we need to find
$k_{0}\geq k_{1}$ such that for all $k>k_{0}$, we have
$\lambda_{1}^{\frac{k}{2}}<\epsilon/d$. Note that if $\epsilon/d\geq 1$, then
we can take $k_{0}=k_{1}$. Hence we assume that $\epsilon/d<1$.

The inequality $\lambda_{1}^{\frac{k}{2}}<\epsilon/d$ is equivalent to
$k\ln\lambda_{1}<2\ln(\epsilon/d)$. Since $\ln\lambda_{1}<0$, the previous
inequality is equivalent to

$k>\frac{2\ln(\epsilon/d)}{\ln\lambda_{1}}=\frac{2\ln(d/\epsilon)}{-\ln\lambda_{1}}.$
(3)

To determine $k_{0}$, we need an upper bound on the right-hand side of (3). We
will use the fact that for any rational $r>1$, $\ln
r<\log_{2}r\leq\log_{2}\lceil r\rceil<\mathrm{bins}(\lceil r\rceil)$, where
$\mathrm{bins}(n)$ is the size of the binary representation of $n$. Thus
$\mathrm{bins}(\lceil r\rceil)$ gives a polynomially computable integer upper
bound for $\ln r$.

Next, using the fact that $\ln(1+x)<x$ for $x\neq 0$, we obtain

$\ln\lambda_{1}=\ln(1+(\lambda_{1}-1))<\lambda_{1}-1,$

which gives $-\ln\lambda_{1}>1-\lambda_{1}$. Hence a polynomially computable
upper bound on the right-hand side of (3) is

$\frac{2\ln(d/\epsilon)}{-\ln\lambda_{1}}<\frac{2\,\mathrm{bins}(\lceil
d/\epsilon\rceil)}{1-\lambda_{1}}.$ (4)

Next we compute a value $k_{1}$ such that for all $k>k_{1}$:

$k^{s}<\lambda_{1}^{-\frac{k}{2}}\qquad\text{or, equivalently,}\qquad C\ln
k<k,$ (5)

where $C=\dfrac{2s}{-\ln\lambda_{1}}$. Using the fact that $\ln(1+x)<x$ for
$x\neq 0$, we obtain $C<\dfrac{2s}{1-\lambda_{1}}$. Hence in order to find
$k_{1}$, we can replace $C$ in (5) with $D=\dfrac{2s}{1-\lambda_{1}}$. In
addition, we can assume that $\ln D>2$, since otherwise we can replace $D$
with a larger value that satisfies this condition, e.g. with $D=9$. Now, Lemma
3.8 implies that every $k>3D\ln D$ satisfies $D\ln k<k$. To make this value
polynomially computable, we can choose it to be

$k_{1}=3\lceil D\rceil\mathrm{bins}(\lceil D\rceil),\text{ where
}D=\max\left\\{\dfrac{2s}{1-\lambda_{1}},\,9\right\\}.$

Finally, combining the right-hand side of (4) with the above formula, we can
define

$k_{0}=\max\left\\{\dfrac{2\,\mathrm{bins}(\lceil
d/\epsilon\rceil)}{1-\lambda_{1}},\ 3\lceil D\rceil\mathrm{bins}(\lceil
D\rceil)\right\\}.$

Note that all the values that appear in the above formula, e.g. $\epsilon$,
$d$ and $D$, can be computed in polynomial time from the input data.

At this point we have derived a polynomially computable $k_{0}$ such that
$\mathcal{P}(a^{k})=u^{T}A^{k}v\in(c-\epsilon,c+\epsilon)$ and, in particular,
$\mathcal{P}(a^{k})\neq\lambda$ for all $k>k_{0}$. Now, to decide the
$\lambda$-reachability problem, we need to check for each integer
$k\in[0,k_{0}]$ whether $u^{T}A^{k}v=\lambda$. Note that the number of
integers in $[0,k_{0}]$ is equal to $2^{\mathrm{bins}(k_{0})}$, which is
exponential in the instance size. Also, computing $A^{k}$ for a given
$k\in[0,k_{0}]$ takes exponential time because
$\mathrm{bins}(A^{k})=\mathcal{O}(2^{\mathrm{bins}(k_{0})}\mathrm{bins}(A))$.
So, the overall algorithm is in EXPTIME.

In a similar way, we can decide the (strict) emptiness problem in EXPTIME. For
instance, suppose $\lambda>c$. Then for all $k>k_{0}$, we have
$\mathcal{P}(a^{k})<c+\epsilon<\lambda$. Thus deciding whether there exists
$k$ such that $\mathcal{P}(a^{k})<\lambda$ is trivial. Suppose we want to know
if there exists $k$ such that $\mathcal{P}(a^{k})\geq\lambda$. In this case,
we need to check for each integer $k\in[0,k_{0}]$ whether
$u^{T}A^{k}v\geq\lambda$. By the same argument as before, this can be done in
EXPTIME.

######  Lemma 3.12.

Let $\mathcal{P}=(u,A,v)$ be a unary polynomially ambiguous probabilistic
finite automaton such that $A$ is upper-triangular and let
$\lambda\in[0,1]\cap{\mathbb{Q}}$ be a cutpoint. Assuming that
$\lambda=\lim_{k\to\infty}\mathcal{P}(a^{k})$, the (strict) emptiness and
$\lambda$-reachability problems for $\mathcal{P}$ and $\lambda$ are decidable
in EXPTIME.

######  Proof 3.13.

Recall that by Lemma 3.6, we can write
$\mathcal{P}(a^{k})=c+\sum_{i=1}^{m}p_{i}(k)\lambda_{i}^{k}$, where
$1>\lambda_{1}>\cdots>\lambda_{m}$ are the eigenvalues of $A$ and $c$ and the
coefficients of $p_{i}$ are rational numbers computable in polynomial time. By
our assumption, $\lambda=\lim_{k\to\infty}\mathcal{P}(a^{k})=c$. As before,
let each $p_{i}(k)$ have the form
$p_{i}(k)=a_{i,s}k^{s}+a_{i,s-1}k^{s-1}+\cdots+a_{i,0}$, where $s\leq n$ (we
do not assume here that $a_{i,s}\neq 0$).

In addition, assume that the leading coefficient of $p_{1}(k)$ is $a_{1,t}$,
for some $t\leq s$. Without loss of generality, suppose $a_{1,t}>0$; the case
when $a_{1,t}<0$ is similar. First, we compute $k_{0}$ such that
$p_{1}(k)>\frac{1}{2}a_{1,t}k^{t}$ for all $k>k_{0}$. To do this, we will use
the following inequalities:

$\displaystyle
a_{1,t}k^{t}+a_{1,t-1}k^{t-1}+\cdots+a_{1,0}>\frac{1}{2}a_{1,t}k^{t}\ \
\Longleftrightarrow\ \
\frac{1}{2}a_{1,t}k^{t}+a_{1,t-1}k^{t-1}+\cdots+a_{1,0}>0$
$\displaystyle\text{and }|a_{1,t-1}k^{t-1}+\cdots+a_{1,0}|\leq
k^{t-1}(|a_{1,t-1}|+\cdots+|a_{1,0}|)=k^{t-1}\sum_{j=0}^{t-1}|a_{1,j}|\quad\text{if
}k\geq 1.$

So, the inequality $p_{1}(k)>\frac{1}{2}a_{1,t}k^{t}$ follows from
$\frac{1}{2}a_{1,t}k^{t}>k^{t-1}\sum_{j=0}^{t-1}|a_{1,j}|$, which is
equivalent to $k>\frac{2}{a_{1,t}}\sum_{j=0}^{t-1}|a_{1,j}|$. Therefore, we
conclude that

$p_{1}(k)>\frac{1}{2}a_{1,t}k^{t}\quad\text{for all $k$ such that
}k>k_{0}:=\max\left\\{1,\
\frac{2}{a_{1,t}}\sum_{j=0}^{t-1}|a_{1,j}|\right\\}.$ (6)

Now we want to find $k_{1}\geq k_{0}$ such that for all $k>k_{1}$, we have

$\lambda_{1}^{k}p_{1}(k)+\lambda_{2}^{k}p_{2}(k)+\cdots+\lambda_{m}^{k}p_{m}(k)>0.$
(7)

Note that

$|\lambda_{2}^{k}p_{2}(k)+\cdots+\lambda_{m}^{k}p_{m}(k)|\leq\lambda_{2}^{k}(|p_{2}(k)|+\cdots+|p_{m}(k)|)\leq
dk^{s}\lambda_{2}^{k},$ (8)

where $d=\sum_{i=2}^{m}\sum_{j=0}^{s}|a_{i.j}|$. Using (6) and (8), we see
that (7) holds whenever $k>k_{0}$ and
$dk^{s}\lambda_{2}^{k}<\frac{1}{2}a_{1,t}k^{t}\lambda_{1}^{k}$, which is
equivalent to

$\displaystyle\frac{2dk^{s-t}}{a_{1,t}}<{\left(\frac{\lambda_{1}}{\lambda_{2}}\right)}^{k}\qquad\text{or}\qquad$
$\displaystyle\ln\frac{2d}{a_{1,t}}+(s-t)\ln
k<k\ln\frac{\lambda_{1}}{\lambda_{2}}$
$\displaystyle\frac{1}{\ln\lambda_{1}/\lambda_{2}}\left(\ln\frac{2d}{a_{1,t}}+(s-t)\ln
k\right)<k.$ (9)

We will use the following inequality

$\ln\frac{\lambda_{1}}{\lambda_{2}}=-\ln\frac{\lambda_{2}}{\lambda_{1}}=-\ln\left(1+\frac{\lambda_{2}-\lambda_{1}}{\lambda_{1}}\right)>-\frac{\lambda_{2}-\lambda_{1}}{\lambda_{1}}>\lambda_{1}-\lambda_{2}.$

Then we can replace (9) with a stronger inequality

$\frac{1}{\lambda_{1}-\lambda_{2}}\left(\ln\frac{2d}{a_{1,t}}+(s-t)\ln
k\right)<k.$ (10)

In the following, we will assume $t<s$ since otherwise (10) simplifies to
$\frac{1}{\lambda_{1}-\lambda_{2}}\ln\frac{2d}{a_{1,t}}<k$. Let us make the
substitution $k=t{\left(\frac{2d}{a_{1,t}}\right)}^{-\frac{1}{s-t}}$, where
$t>0$. Then (10) can be written as

$\displaystyle\frac{1}{\lambda_{1}-\lambda_{2}}\left(\ln\frac{2d}{a_{1,t}}+(s-t)\ln
t+(s-t)\frac{-1}{s-t}\ln\frac{2d}{a_{1,t}}\right)$
$\displaystyle<t{\left(\frac{2d}{a_{1,t}}\right)}^{-\frac{1}{s-t}}$
$\displaystyle{\left(\frac{2d}{a_{1,t}}\right)}^{\frac{1}{s-t}}\frac{s-t}{\lambda_{1}-\lambda_{2}}\ln
t$ $\displaystyle<t.$

Let $D=\max\left\\{9,\
{\left(\frac{2d}{a_{1,t}}\right)}^{\frac{1}{s-t}}\frac{s-t}{\lambda_{1}-\lambda_{2}}\right\\}$.
Here $9$ is needed to satisfy the requirement $\ln D>2$ in Lemma 3.8. Then by
Lemma 3.8, the above inequality holds when $t>3D\ln D$. Therefore, (10) and
hence (9) holds when
$k>3{\left(\frac{2d}{a_{1,t}}\right)}^{-\frac{1}{s-t}}D\ln D$. To make this
bound polynomially computable, we can simplify it as follows. Suppose that
$2d\geq a_{1,t}$. Then (9) holds when

$k>k_{1}:=3\lceil E\rceil\mathrm{bins}(\lceil E\rceil),\quad\text{where
}E=\max\left\\{9,\
\frac{2d}{a_{1,t}}\cdot\frac{s-t}{\lambda_{1}-\lambda_{2}}\right\\}$

because in this case
${\left(\frac{2d}{a_{1,t}}\right)}^{\frac{1}{s-t}}\leq\frac{2d}{a_{1,t}}$ and
${\left(\frac{2d}{a_{1,t}}\right)}^{-\frac{1}{s-t}}\leq 1$. On the other hand,
if $2d<a_{1,t}$, then (9) holds when

$k>k_{1}:=3\left\lceil\frac{a_{1,t}}{2d}E\right\rceil\mathrm{bins}(\lceil
E\rceil),\quad\text{where }E=\max\left\\{9,\
\frac{s-t}{\lambda_{1}-\lambda_{2}}\right\\}$

because in this case
${\left(\frac{2d}{a_{1,t}}\right)}^{-\frac{1}{s-t}}<{\left(\frac{2d}{a_{1,t}}\right)}^{-1}$
and ${\left(\frac{2d}{a_{1,t}}\right)}^{\frac{1}{s-t}}<1$.

Finally, we conclude that (7) holds for all $k>k_{2}:=\max\\{k_{0},k_{1}\\}$,
where both $k_{0}$ and $k_{1}$ are computable in PTIME. In other words, we
obtained a polynomially computable value $k_{2}$ such that
$\mathcal{P}(a^{k})>c=\lambda$ for all $k>k_{2}$. Using the same argument as
at the end of the proof of Lemma 3.10, we can show that the (strict) emptiness
and $\lambda$-reachability problems are decidable in EXPTIME.

We are now ready to give a proof of Theorem 3.3.

######  Proof 3.14 (Proof of Theorem 3.3).

Let $\mathcal{P}=(u,A,v)$ be a polynomially ambiguous unary PFA. By Lemma 3.4,
we can compute in EXPTIME a set of $d$ polynomially ambiguous unary PFAs
$\\{\mathcal{P}_{s}=(u_{s},U,v^{\prime})\ |\ 0\leq s\leq d-1\\}$ such that $U$
is rational upper-triangular and

$\mathcal{P}(a^{rd+s})=\mathcal{P}_{s}(a^{r})=u_{s}^{T}U^{r}v^{\prime},$

where $0\leq s\leq d-1$.

Suppose $\lambda$ is a given cutpoint. If we want to decide whether there
exists $k$ such that $\mathcal{P}(a^{k})=\lambda$ (or
$\mathcal{P}(a^{k})\geq\lambda$), we can check for every $s$ from $0$ to $d-1$
whether there exists $r$ such that $\mathcal{P}_{s}(a^{r})=\lambda$ (or
$\mathcal{P}_{s}(a^{r})\geq\lambda$, respectively), which can be done in
EXPTIME using Lemmas 3.10 and 3.12. Namely, we will use Lemma 3.10 if
$\lambda\neq c_{s}$ and Lemma 3.12 if $\lambda=c_{s}$ for the current values
of $s\in[0,d-1]$. Finally, we note that even though the value of $d$ can be
exponential in the input size, the whole procedure can still be done in
EXPTIME.

Skolem’s problem is at least NP-hard [6] implying that the
$\lambda$-reachability and emptiness problems are also NP-hard, at least for
PFA of exponential ambiguity. Our next result shows that NP-hardness can be
established even for unary PFAs of _finite ambiguity_.

######  Theorem 3.15.

The $\lambda$-reachability and emptiness problems for unary finitely ambiguous
Probabilistic Finite Automata $\mathcal{P}=(u,A,v)$ with $\\{0,1\\}$-matrix
$A$ are NP-hard.

######  Proof 3.16.

The NP-hardness of Skolem’s problem was established in [6]. Specifically,
Corollary 1.3 of [6] states that the problem of determining, for a given
matrix $A\in\\{0,1\\}^{n\times n}$ and row vectors $b,c\in\\{0,1\\}^{n}$, if
$b^{T}A^{k}c=0$ for some $k\geq 0$ is NP-hard. Examination of the proof of
this corollary shows that in fact $\mathcal{P}$ is finitely ambiguous as we
shall show.

The proof of Theorem 1.1 of [6] shows a reduction of 3SAT on $m$ clauses with
$n$ letters to a unary rational expression $E$ of the form:

$E=\bigcup_{j=0}^{k}a^{z_{j}}(a^{r_{j}})^{*},$

where $k=\mathcal{O}(n^{3}m)$ and $z_{j},r_{j}=\mathcal{O}(n^{6})$ as is not
difficult to see from the proof in [6]. Notice then that each $z_{j},r_{j}$
represented in unary has a polynomial size in terms of the 3SAT instance and
thus $E$ also has a polynomial representation size.

We may then invoke Kleene’s theorem [15] to state that the language recognised
by $E$ is also recognised by an NFA $\mathcal{P}=(b,\\{A\\},c)$ which thus
allows the derivation of Corollary 1.3 of [6]. Note that $E$ is simply the
union of rational expressions of the form $E_{j}=a^{z_{j}}(a^{r_{j}})^{*}$.
Each $E_{j}$ can be transformed to an NFA $\mathcal{N}_{j}$ with
$z_{j}+r_{j}+1$ states
$S_{j}=n_{0,j},\ldots,n_{z_{j},j},n_{z_{j}+1,j},\ldots,n_{z_{j}+r_{j},j}$ with
initial state $n_{0,j}$, final state $n_{z_{j}+1,j}$ and transition function
$\delta:S_{j}\times\\{a\\}\to S_{j}$ given by $\delta(n_{i,j},a)=n_{i+1,j}$
for $0\leq i\leq z_{j}+r_{j}-1$ and $\delta(n_{z_{j}+r_{j}},a)=n_{z_{j}+1,j}$.

We may then form an NFA $\mathcal{N}$ by
$\mathcal{N}=\bigcup_{j=0}^{k}\mathcal{N}_{j}$ with the usual construction. In
this case, $\mathcal{N}$ has set of initial states $\\{n_{0,j}\ |\ 1\leq j\leq
k\\}$, set of final states $\\{n_{z_{j}+1,j}\ |\ 1\leq j\leq k\\}$ and states
in disjoint subsets $S_{j}$ and $S_{j^{\prime}}$ with $j\neq j^{\prime}$ are
not connected. This implies by the IDA property of [23] that $\mathcal{N}$ is
finitely ambiguous since there does not exist any state with two outgoing
transitions (by which reasoning we also know that each row of $\mathcal{N}$’s
transition matrix has exactly one entry $1$ with all others $0$). In fact one
may see that $\mathcal{N}$ is $k$-ambiguous with $k=\mathcal{O}(n^{3}m)$. The
number of states of $\mathcal{N}$ is
$d=\sum_{j=0}^{k}z_{j}+r_{j}+1=\mathcal{O}(n^{9}m)$ which is polynomial in the
3SAT instance representation size.

We note that actually $\mathcal{N}$ is already close to a PFA. Since each row
is zero except for exactly one entry $1$, matrix $A$ is stochastic. We thus
consider Probabilistic Finite Automaton $\mathcal{P}=(u,\\{A\\},c)$ where
$u=\frac{b}{|b|}$ is the initial (stochastic) vector. $\mathcal{P}$ has
polynomial ambiguity since $\mathcal{N}$ does. Therefore, deciding if there
exists $k\geq 0$ such that $\mathcal{P}(a^{k})=0$ or $\mathcal{P}(a^{k})\leq
0$ is NP-hard to determine, proving NP-hardness of the $\lambda$-reachability
and emptiness problems. Since we did not modify $\mathcal{N}$ to derive
$\mathcal{P}$ other than to scale the initial vector, the degree of ambiguity
is retained.

######  Corollary 1.

The $\lambda$-reachability and emptiness problems for unary polynomially
ambiguous PFA $\mathcal{P}=(u,A,v)$ with $\\{0,1\\}$-matrix $A$ are NP-
complete.

######  Proof 3.17.

NP-hardness follows from Theorem 3.15 since finite ambiguity is a stronger
property than polynomially ambiguity. To prove the NP upper bound, we will
show that the algorithm in the proof of Theorem 3.3 can be done in NP. We
again use Lemmas 3.4, 3.6, 3.10 and 3.12. Note that the value $d$ from Lemma
3.4 can be exponential. However, its binary presentation has polynomial size.
So, instead of cycling though all $s$ from $0$ to $d-1$, we can
nondeterministically guess in polynomial time a value $s\in[0,d-1]$.

Next, we note that the values of $k_{0}$ in Lemma 3.10 and $k_{2}$ in Lemma
3.12 also have binary representations of polynomial size. Again, instead of
checking every $k$ in $[0,k_{0}]$ or $[0,k_{2}]$, we can nondeterministically
guess $k$ in polynomial time.

Finally, in the verification step of our algorithm we need to compute the
matrices $A^{d}$, $A^{s}$ and $(A^{d})^{k}$. This can be done in polynomial
time using exponentiation by squaring. Indeed, the exponentiation by squaring
requires polynomially many steps. Also, any power of a stochastic
$\\{0,1\\}$-matrix is also a stochastic $\\{0,1\\}$-matrix, so the entries of
the power matrices do not grow in size.

##  4 Binary PFA

The following theorem shows that the $\lambda$-reachability and emptiness
problems are NP-hard for binary PFA of _polynomial ambiguity_ with a _fixed_
and _commuting_ transition matrices. The emptiness problem for non-commutative
binary PFA over $25$ states is known to be undecidable, at least over
exponentially ambiguous PFA [12]. Emptiness is also undecidable for
exponentially ambiguous commutative PFA, although with many more states and a
larger alphabet [2].

######  Theorem 4.1.

The $\lambda$-reachability and emptiness problems are NP-hard for binary
probabilistic finite automata of polynomial ambiguity with fixed and commuting
matrices of dimension $9$ for $\lambda$-reachability, $37$ for nonstrict
emptiness, and $40$ for strict emptiness.

######  Proof 4.2.

We use a reduction from the solvability of binary quadratic Diophantine
equations. Namely, given an equation of the form $ax^{2}+by-c=0$, where
$a,b,c\in\mathbb{N}$, it is NP-hard to determine if there exists
$x,y\in\mathbb{N}$ satisfying the equation [16]. We begin with the
$\lambda$-reachability problem before considering the emptiness problem.

$\lambda$-Reachability reduction. Let $A=\begin{pmatrix}1&1\\\
0&1\end{pmatrix}$ and note that $A^{k}=\begin{pmatrix}1&k\\\ 0&1\end{pmatrix}$
and that $(A\otimes A)^{k}_{1,4}=(A^{k}\otimes A^{k})_{1,4}=k^{2}$. We form a
weighted automaton222For our purposes here, by a _weighted automaton_ we
simply mean an automaton whose initial vector, final vector, and transition
matrices are over nonnegative integers. $W_{1}$ on binary alphabet
$\Sigma=\\{h,g\\}$ in the following way to encode $ax^{2}+by$ (we will deal
with $c$ later). Let $W_{1}=(u_{1},\phi,v_{1})$ where
$u_{1},v_{1}\in\mathbb{N}^{7}$ and $\phi:\Sigma^{*}\to\mathbb{N}^{7\times 7}$.
We define $u_{1}=(a,0,0,0,b,0,0)^{T}$, $v_{1}=(0,0,0,1,0,1,0)^{T}$ and
$\phi(\ell)=\frac{1}{4}\phi^{\prime}(\ell)$ for $\ell\in\\{h,g\\}$ with

$\phi^{\prime}(h)=\left(\begin{array}[]{@{}c|l@{}}(A\otimes A)\oplus
I_{2}&t_{1}\\\ \hline\cr{\bf
0}^{6}&4\end{array}\right),\quad\phi^{\prime}(g)=\left(\begin{array}[]{@{}c|l@{}}I_{4}\oplus
A&t_{2}\\\ \hline\cr{\bf 0}^{6}&4\end{array}\right),$

with ${\bf 0}^{k}=(0,0,\ldots,0)\in\mathbb{N}^{k}$, $t_{1}=(0,2,2,3,3,3)^{T}$
and $t_{2}=(3,3,3,3,2,3)^{T}$. We see then that each row of
$\phi^{\prime}(\ell)$ is nonnegative and sums to $4$, thus $\phi(\ell)$ is
stochastic for $\ell\in\\{g,h\\}$. Furthermore, by the mixed product property
of the Kronecker product, we see that $((A\otimes A)\oplus
I_{2})^{x}=(A^{x}\otimes A^{x})\oplus I_{2}$ and $(I_{4}\oplus
A)^{y}=I_{4}\oplus A^{y}$ for $x,y\in{\mathbb{N}}$ and thus by the block upper
triangular structure of $\phi^{\prime}(h),\phi^{\prime}(g)$, we see that

$\phi^{\prime}(h^{x}g^{y})=\left(\begin{array}[]{@{}c|l@{}}(A^{x}\otimes
A^{x})\oplus A^{y}&t_{xy}\\\ \hline\cr{\bf 0}^{6}&4^{x+y}\end{array}\right),$

where $t_{xy}$ is a nonnegative vector maintaining the row sum at $4^{x+y}$.
We now see that

$\displaystyle u_{1}^{T}\phi(h^{x}g^{y})v_{1}=\frac{ax^{2}+by}{4^{x+y}}$ (11)

We define a second weighted automaton $W_{2}=(u_{2},\psi,v_{2})$ with
$u_{2}=(c,0)^{T}$, $v_{2}=(0,1)^{T}$ and
$\psi:\Sigma^{*}\to{\mathbb{N}}^{2\times 2}$ with
$\psi(\ell)=\frac{1}{4}\psi^{\prime}(\ell)$ for $\ell\in\\{h,g\\}$ defined
thus: $\psi^{\prime}(h)=\psi^{\prime}(g)=\begin{pmatrix}1&3\\\
0&4\end{pmatrix}$. We therefore see that

$\displaystyle
u_{2}^{T}\psi(h^{x}g^{y})v_{2}=\frac{c(4^{x+y}-1)}{4^{x+y}}=c(1-\frac{1}{4^{x+y}})$
(12)

We now join $W_{1}$ and $W_{2}$ into a $9$-state PFA
$\mathcal{P}=(u,\gamma,v)$ where $u=\frac{1}{a+b+c}[u_{1}|u_{2}]$,
$v=[v_{1}|v_{2}]$ and $\gamma(\ell)=\phi(\ell)\oplus\psi(\ell)$. Combining
Eqns (11) and (12) we see that

$\displaystyle u^{T}\gamma(h^{x}g^{y})v$ $\displaystyle=$
$\displaystyle\frac{1}{a+b+c}\left(\frac{ax^{2}+by}{4^{x+y}}+c(1-\frac{1}{4^{x+y}})\right)$
(13) $\displaystyle=$
$\displaystyle\frac{1}{a+b+c}\left(c+\frac{ax^{2}+by-c}{4^{x+y}}\right)$

which equals $\frac{c}{a+b+c}$ if and only if $ax^{2}+by-c=0$. Note that
$\gamma(h)$ and $\gamma(g)$ commute by their structure since clearly
$(A\otimes A)\oplus I$ and $I_{4}\oplus A$ commute, giving $(A\otimes A)\oplus
A$ in both cases (as a consequence of the mixed product properties of Lemma
2.1) and the rightmost vector of the matrix simply retains the row sum at $1$
for such a product since the matrices are stochastic. Both $\gamma(h)$ and
$\gamma(g)$ are upper-triangular thus $\mathcal{P}$ is polynomially ambiguous.

Nonstrict Emptiness reduction. We now show the proof of the emptiness problem.
We showed that the $\lambda$-reachability problem is NP-hard by deriving a PFA
${\mathcal{P}}$ over the binary alphabet $\\{h,g\\}$ such that
${\mathcal{P}}(h^{x}g^{y})$ is given by Eqn. 13. We note however that a non
solution to $ax^{2}+by-c=0$ can be positive or negative and thus we may be
above or below the threshold $\frac{c}{a+b+c}$. This encoding thus cannot be
used to show the NP-hardness of the _emptiness_ problem.

Instead, we can use a similar encoding of the quartic polynomial given by
$(ax^{2}+by-c)^{2}=a^{2}x^{4}+2abx^{2}y+b^{2}y^{2}+c^{2}-2acx^{2}-2bcy$ with
$a,b,c\in{\mathbb{N}}$. Note that we arranged the four positive terms first,
followed by the two negative terms. Clearly $(ax^{2}+by-c)^{2}$ is nonnegative
and equals zero if and only if $ax^{2}+by-c=0$. We will derive a PFA
${\mathcal{P}}_{2}$ such that

${\mathcal{P}}_{2}(h^{x}g^{y})=\frac{1}{z}\left((2ac+2bc)+\frac{1}{16^{x+y}}(ax^{2}+by+c)^{2}\right),$

where $z=a^{2}+2ab+b^{2}+c^{2}+2ac+2bc$, with the property that
${\mathcal{P}}_{2}(h^{x}g^{y})\geq\frac{2ac+2bc}{z}$ with equality if and only
if $(ax^{2}+by-c)^{2}=0$ which is NP-hard to determine. To this end, we
compute the following four matrices $\\{H_{+},G_{+},H_{-},G_{-}\\}$, the idea
being that $H_{+}$ and $G_{+}$ will be used to compute the positive four terms
and $H_{-}$ and $G_{-}$ will compute the negative terms:

$\displaystyle H_{+}$ $\displaystyle=$ $\displaystyle\underbrace{(A\otimes
A\otimes A\otimes A)}_{x^{4}}\oplus\underbrace{(A\otimes A\otimes
I_{2})}_{x^{2}y}\oplus\underbrace{(I_{2}\otimes
I_{2})}_{y^{2}}\oplus\underbrace{1}_{1}$ $\displaystyle G_{+}$
$\displaystyle=$ $\displaystyle\underbrace{(I_{2}\otimes I_{2}\otimes
I_{2}\otimes I_{2})}_{x^{4}}\oplus\underbrace{(I_{2}\otimes I_{2}\otimes
A)}_{x^{2}y}\oplus\underbrace{(A\otimes A)}_{y^{2}}\oplus\underbrace{1}_{1}$
$\displaystyle H_{-}$ $\displaystyle=$ $\displaystyle\underbrace{(A\otimes
A)}_{x^{2}}\oplus\underbrace{I_{2}}_{y}$ $\displaystyle G_{-}$
$\displaystyle=$ $\displaystyle\underbrace{(I_{2}\otimes
I_{2})}_{x^{2}}\oplus\underbrace{A}_{y}$

and by the mixed product property of Kronecker products of Lemma 2.1),

$\displaystyle H_{+}^{x}G_{+}^{y}$ $\displaystyle=$
$\displaystyle(A^{x}\otimes A^{x}\otimes A^{x}\otimes
A^{x})\oplus(A^{x}\otimes A^{x}\otimes A^{y})\oplus(A^{y}\otimes A^{y})\oplus
1$ $\displaystyle H_{-}^{x}G_{-}^{y}$ $\displaystyle=$
$\displaystyle(A^{x}\otimes A^{x})\oplus A^{y}$

Note that $H_{+}^{x}G_{+}^{y}$ and $H_{-}^{x}G_{-}^{y}$ each contain the
positive and negative (respectively) term of $(ax^{2}+by-c)^{2}$, excluding
the coefficients, e.g. $(H_{+}^{x}G_{+}^{y})_{1,16}=x^{4}$ and
$(H_{+}^{x}G_{+}^{y})_{17,24}=x^{2}y$ etc. Note also that
$H_{+}G_{+}=G_{+}H_{+}$ and $H_{-}G_{-}=G_{-}H_{-}$ which also follows from
the mixed product properties and thus matrices $\\{H_{+},G_{+}\\}$ and
$\\{H_{-},G_{-}\\}$ commute.

As before, we may now increase the dimension of each matrix
$\\{H_{+},H_{-},G_{+},G_{-}\\}$ by $1$ to ensure a common row sum (of $16$ in
this case) by adding a new column on the right hand side of each matrix, and
then divide each matrix by this common value to give
$\\{H^{\prime}_{+},H^{\prime}_{-},G^{\prime}_{+},G^{\prime}_{-}\\}$ so that
each of these matrices is row stochastic. Matrices
$\\{H^{\prime}_{+},G^{\prime}_{+}\\}$ and
$\\{H^{\prime}_{-},G^{\prime}_{-}\\}$ still commute since this change only has
an effect on the final column of the matrix.

We now show how to handle each term of $(ax^{2}+by-c)^{2}$. We first handle
the positive terms. We define
$u_{1}=(a^{2},0,\ldots,0)^{T}\in{\mathbb{Q}}^{16}$,
$u_{2}=(2ab,0,\ldots,0)^{T}\in{\mathbb{Q}}^{8}$,
$u_{3}=(b^{2},0,0,0)^{T}\in{\mathbb{Q}}^{4}$ and $u_{4}=c^{2}$ and then let
$u_{+}=[u_{1}|u_{2}|u_{3}|u_{4}|0]\in{\mathbb{Q}}^{30}$. We let
$v_{1}=(0,\ldots,0,1)^{T}\in{\mathbb{Q}}^{16}$,
$v_{2}=(0,\ldots,0,1)^{T}\in{\mathbb{Q}}^{8}$,
$v_{3}=(0,0,0,1)^{T}\in{\mathbb{Q}}^{4}$ and $v_{4}=1$, and let
$v_{+}=[v_{1}|v_{2}|v_{3}|v_{4}|0]\in{\mathbb{Q}}^{30}$. We then see that

$\displaystyle u_{+}^{T}(H^{\prime}_{+})^{x}(G^{\prime}_{+})^{y}v_{+}$ (14)
$\displaystyle=$ $\displaystyle\frac{1}{16^{x+y}}\left(u_{1}^{T}(A^{x}\otimes
A^{x}\otimes A^{x}\otimes A^{x})v_{1}+u_{2}^{T}(A^{x}\otimes A^{x}\otimes
A^{y})v_{2}+u_{3}^{T}(A^{y}\otimes A^{y})v_{3}+u_{4}^{T}v_{4}\right)$
$\displaystyle=$
$\displaystyle\frac{1}{16^{x+y}}\left(a^{2}x^{4}+2abx^{2}y+b^{2}y^{2}+c^{2}\right)$

We next handle the negative terms, which is essentially accomplished by
switching final and non-final states in the final state vectors to follow.
Define $u_{5}=(2ac,0,0,0)^{T}\in{\mathbb{Q}}^{4}$ and
$u_{6}=(2bc,0)^{T}\in{\mathbb{Q}}^{2}$ and let
$u_{-}=[u_{5}|u_{6}|0]\in{\mathbb{Q}}^{7}$. We let
$v_{5}=(0,0,0,1)^{T}\in{\mathbb{Q}}^{4}$ and
$v_{6}=(0,1)^{T}\in{\mathbb{Q}}^{2}$. Define
$v_{-}=[v_{5}|v_{6}|0]\in{\mathbb{Q}}^{7}$. We then see that

$\displaystyle u_{-}^{T}(H^{\prime}_{-})^{x}(G^{\prime}_{-})^{y}({\bf
1}-v_{-})$ (15) $\displaystyle=$
$\displaystyle(2ac+2bc)-\frac{1}{16^{x+y}}\left(u_{5}^{T}(A^{x}\otimes
A^{x})v_{5}+u_{6}^{T}A^{y}v_{6}+0\right)$ $\displaystyle=$
$\displaystyle(2ac+2bc)-\frac{1}{16^{x+y}}\left(2acx^{2}+2bcy\right),$

where ${\bf 1}=(1,1,\ldots,1)^{T}\in{\mathbb{Q}}^{7}$. We used here the fact
that $X{\bf 1}={\bf 1}$ for a row stochastic matrix $X$. We finally define
that $H=H^{\prime}_{+}\oplus H^{\prime}_{-}\in\mathbb{Q}^{37\times 37}$ and
$G=G^{\prime}_{+}\oplus G^{\prime}_{-}\in\mathbb{Q}^{37\times 37}$, both of
which are row stochastic and commute, and let
$u_{\star}=\frac{[u_{+}|u_{-}]}{z}\in\mathbb{Q}^{37}$ and
$v_{\star}=[v_{+}|({\bf 1}-v_{-})]\in\mathbb{Q}^{37}$, with
$z=a^{2}+2ab+b^{2}+c^{2}+2ac+2bc$ to normalise vector $u_{\star}$. We see then
that $u_{\star}$ is a stochastic vector as required. We define the PFA
${\mathcal{P}}_{2}=(u_{\star},\\{H,G\\},v_{\star})$ and we can now compute
that

$\displaystyle\mathcal{P}_{2}(h^{x}g^{y})=u_{\star}^{T}H^{x}G^{y}v_{\star}$
$\displaystyle=$ $\displaystyle u_{\star}^{T}(H_{+}^{\prime x}G_{+}^{\prime
y}\oplus H_{-}^{\prime x}G_{-}^{\prime y})v_{\star}$ $\displaystyle=$
$\displaystyle\frac{1}{16^{x+y}}\left(\frac{[u_{+}|u_{-}]}{z}^{T}\left(\begin{array}[]{c|c}\begin{array}[]{c|c}H_{+}^{x}G_{+}^{y}&*\\\
\hline\cr{\bf 0}&16^{x+y}\end{array}&{\bf 0}\\\ \hline\cr{\bf
0}&\begin{array}[]{c|c}H_{-}^{x}G_{-}^{y}&*\\\ \hline\cr{\bf
0}&16^{x+y}\end{array}\end{array}\right)[v_{+}|({\bf 1}-v_{-})]\right)$

$\displaystyle=$
$\displaystyle\frac{1}{z16^{x+y}}\left(u_{+}^{T}H_{+}^{x}G_{+}^{y}v_{+}+u_{-}^{T}H_{-}^{x}G_{-}^{y}({\bf
1}-v_{-})\right)$ (19) $\displaystyle=$
$\displaystyle\frac{1}{z}\left(u_{+}^{T}(H^{\prime}_{+})^{x}(G^{\prime}_{+})^{y}v_{+}+u_{-}^{T}(H^{\prime}_{-})^{x}(G^{\prime}_{-})^{y}({\bf
1}-v_{-})\right)$ $\displaystyle=$
$\displaystyle\frac{1}{z}\left((2ac+2bc)+\frac{1}{16^{x+y}}\left(a^{2}x^{4}+2abx^{2}y+b^{2}y^{2}+c^{2}\right)-\frac{1}{16^{x+y}}\left(2acx^{2}+2bcy\right)\right)$
$\displaystyle=$
$\displaystyle\frac{1}{z}\left((2ac+2bc)+\frac{1}{16^{x+y}}(ax^{2}+by-c)^{2}\right)$

where $*$ denote the column vectors used to ensure row sums of $16^{x+y}$ and
${\bf 0}$ denotes zero matrices of appropriate sizes. We also used Eqns (14)
and (15).

Since $(ax^{2}+by-c)^{2}$ is nonnegative, we see that
$u_{\star}^{T}H^{x}G^{y}v_{\star}\geq\frac{2ac+2bc}{z}$ with equality if and
only if $(ax^{2}+by-c)^{2}=0$, which is NP-hard to determine. Therefore using
cutpoint $\lambda=\frac{2ac+2bc}{z}\in{\mathbb{Q}}\cap[0,1]$ means the
(nonstrict) emptiness problem is NP-hard (i.e. does there exist
$x,y\in{\mathbb{N}}$ such that $u_{\star}^{T}H^{x}G^{y}v_{\star}\leq\lambda$
is NP-hard). As before, matrices $H$ and $G$ are upper-triangular and commute
by their structure, and therefore the result holds.

Strict Emptiness reduction. Finally we show how to handle the strict emptiness
problem. We proceed with a technique inspired by [10]. By (19), if
$P_{2}(h^{x}g^{y})=u_{\star}^{T}H^{x}G^{y}v_{\star}\neq\frac{1}{z}(2ac+2bc)$,
then
$u_{\star}^{T}H^{x}G^{y}v_{\star}\geq\frac{1}{z}\left((2ac+2bc)+\frac{1}{16^{x+y}}\right)$
therefore $P_{2}(h^{x}g^{y})\leq\frac{1}{z}(2ac+2bc)$ if and only if
$P_{2}(h^{x}g^{y})<\frac{1}{z}\left((2ac+2bc)+\frac{1}{16^{x+y}}\right)$.

Let us adapt ${\mathcal{P}}_{2}$ in the following way to create a new PFA
${\mathcal{P}}_{3}$. Note that ${\mathcal{P}}_{2}$ has $6$ initial states (by
$u_{\star}$). We add three new states to ${\mathcal{P}}_{3}$, denoted
$q_{0},q_{F}$ and $q_{*}$. State $q_{0}$ is a new initial state of
${\mathcal{P}}_{3}$ which, for any input letter, has probability
$\frac{1}{2\cdot 6}$ of moving to each of the $6$ initial states of
${\mathcal{P}}_{2}$ and probability $\frac{1}{2}$ to move to new state
$q_{F}$. State $q_{F}$ is a new final state that remains in $q_{F}$ for any
input letter with probability $1-\frac{1}{16z}$ and moves to a new non-
accepting absorbing sink state $q_{*}$ with probability $\frac{1}{16z}$. We
now see that for any $a\in\\{h,g\\}$:

${\mathcal{P}}_{3}(aw)=\frac{1}{2}{\mathcal{P}}_{2}(w)+\frac{1}{2}\left(1-\frac{1}{16^{|w|}z^{|w|}}\right)$

If there exists $w_{1}=h^{x}g^{y}$ with $x,y\geq 0$ such that
${\mathcal{P}}_{2}(w_{1})\leq\frac{1}{z}(2ac+2bc)$ then
${\mathcal{P}}_{2}(w_{1})=\frac{1}{z}(2ac+2bc)$ and thus:

${\mathcal{P}}_{3}(aw_{1})=\frac{1}{2}\left(\frac{1}{z}(2ac+2bc)\right)+\frac{1}{2}\left(1-\frac{1}{16^{|w_{1}|}z^{|w_{1}|}}\right)<\frac{1}{2}\left(\frac{1}{z}(2ac+2bc)+1\right).$

For any $w_{2}=h^{x}g^{y}$ with $x,y\geq 0$ such that
${\mathcal{P}}_{2}(w_{2})>\frac{1}{z}(2ac+2bc)$ then
${\mathcal{P}}_{2}(w_{2})\geq\frac{1}{z}(2ac+2bc)+\frac{1}{16^{x+y}}$ by (19).
Thus:

${\mathcal{P}}_{3}(aw_{2})\geq\frac{1}{2}\left(\frac{1}{z}(2ac+2bc)+\frac{1}{16^{|w_{2}|}}\right)+\frac{1}{2}\left(1-\frac{1}{16^{|w_{2}|}z^{|w_{2}|}}\right)>\frac{1}{2}\left(\frac{1}{z}(2ac+2bc)+1\right).$

Thus determining if there exists $w=h^{x}g^{y}$ such that
${\mathcal{P}}_{3}(w)<\frac{1}{2}\left(\frac{1}{z}(2ac+2bc)+1\right)$, i.e.
the strict emptiness problem for ${\mathcal{P}}_{3}$ on cutpoint
$\frac{1}{2}\left(\frac{1}{z}(2ac+2bc)+1\right)$, is NP-hard. The
modifications to ${\mathcal{P}}_{2}$ retain polynomial ambiguity since $q_{0}$
and $q_{F}$ have no incoming (non self looping) edges and $q_{*}$ has no
outgoing edges, therefore property EDA does not hold. Commutativity of the PFA
is unaffected since $\mathcal{P}_{3}$ is identical to $\mathcal{P}_{2}$ except
for adding three new states, behaving identically for both input letters. Note
that ${\mathcal{P}}_{3}$ has $37+3=40$ states.

## References

  * [1] S. Akshay, T. Antonopoulos, J. Ouaknine, and J. Worrell.  Reachability problems for markov chains.  Information Processing Letters, 115(2):155–158, 2015.
  * [2] P. C. Bell.  Polynomially ambiguous probabilistic automata on restricted languages.  In International Colloquium on Automata, Languages, and Programming, number 105 in ICALP’19, pages 1–14, 2019.
  * [3] P. C. Bell and P. Semukhin.  Decidability of cutpoint isolation for probabilistic finite automata on letter-bounded inputs.  In International Conference on Concurrency Theory, number 22 in CONCUR’20, pages 1–16, 2020.
  * [4] A. Bertoni, G. Mauri, and M. Torelli.  Some recursively unsolvable problems relating to isolated cutpoints in probabilistic automata.  In Automata, Languages and Programming, volume 52, pages 87–94, 1977.
  * [5] V. Blondel and V. Canterini.  Undecidable problems for probabilistic automata of fixed dimension.  Theory of Computing Systems, 36:231–245, 2003.
  * [6] V. D. Blondel and N. Portier.  The presence of a zero in an integer linear recurrent sequence is np-hard to decide.  Linear Algebra and its Applications, pages 91–98, 2002.
  * [7] L. Daviaud, M. Jurdzinski, R. Lazic, F. Mazowiecki, G. A. Pérez, and J. Worrell.  When is containment decidable for probabilistic automata?  In International Colloquium on Automata, Languages, and Programming, number 121 in ICALP’18, pages 1–14, 2018.
  * [8] N. Fijalkow, C. Riveros, and J. Worrell.  Probabilistic automata of bounded ambiguity.  In 28th International Conference on Concurrency Theory (CONCUR), pages 19:1–19:14, 2017.
  * [9] S. Friedland.  Matrices: Algebra, Analysis and Applications.  World Scientific Publishing Company Pte Limited, 2015.
  * [10] H. Gimbert and Y. Oualhadj.  Probabilistic automata on finite words: decidable and undecidable problems.  In International Colloquium on Automata, Languages and Programming (ICALP’10), volume 2, pages 527–538, 2010.
  * [11] V. Halava, T. Harju, M. Hirvensalo, and J. Karhumäki.  Skolem’s problem — on the border between decidability and undecidability.  In TUCS Technical Report Number 683, 2005.
  * [12] M. Hirvensalo.  Improved undecidability results on the emptiness problem of probabilistic and quantum cut-point languages.  SOFSEM 2007: Theory and Practice of Computer Science, Lecture Notes in Computer Science, 4362:309–319, 2007.
  * [13] R. A. Horn and C. R. Johnson.  Topics in matrix analysis.  Cambridge University Press, 1991.
  * [14] O. Ibarra and B. Ravikumar.  On sparseness, ambiguity and other decision problems for acceptors and transducers.  In Proc. STACS 1986, volume 210, pages 171–179, 1986.
  * [15] S. C. Kleene.  Representation of events in nerve nets and finite automata.  Automata Studies, Annals of Mathematical Studies, 34, 1956.
  * [16] K. L. Manders and L. Adleman.  NP-complete decision problems for binary quadratics.  Journal of Computer and System Sciences, 16:168–184, 1978.
  * [17] J. Ouaknine and J. Worrell.  Positivity problems for low-order linear recurrence sequences.  In Proceedings of the Twenty-Fifth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA’14, pages 366–379, SODA, 2014.
  * [18] A. Paz.  Introduction to Probabilistic Automata.  Academic Press, 1971.
  * [19] M. O. Rabin.  Probabilistic automata.  Information and Control, 6:230–245, 1963.
  * [20] C. Reutenauer.  Propertiétés arithmétiques et topologiques de séries rationnelles en variables non commutatives.  Thèse troisième cycle, Université Paris VI, 1977.
  * [21] P. Turakainen.  Generalized automata and stochastic languages.  Proceedings of the American Mathematical Society, 21:303–309, 1969\.
  * [22] N. K. Vereshchagin.  The problem of appearance of a zero in a linear recurrence sequence (in russian).  Mat. Zametki, 38(2), 1985.
  * [23] A. Weber and H. Seidl.  On the degree of ambiguity of finite automata.  Theoretical Computer Science, 88(2):325–349, 1991.


#  Ergodic Rate Analysis of Reconfigurable Intelligent Surface-Aided Massive
MIMO Systems with ZF Detectors

Kangda Zhi, Cunhua Pan, Hong Ren and Kezhi Wang (Corresponding author: Cunhua
Pan.)K. Zhi, C. Pan are with the School of Electronic Engineering and Computer
Science at Queen Mary University of London, London E1 4NS, U.K. (e-mail:
k.zhi, c.pan@qmul.ac.uk). H. Ren is with the National Mobile Communications
Research Laboratory, Southeast University, Nanjing 210096, China.
(hren@seu.edu.cn). K. Wang is with Department of Computer and Information
Sciences, Northumbria University, UK. (e-mail: kezhi.wang@northumbria.ac.uk).

###### Abstract

This letter investigates the reconfigurable intelligent surface (RIS)-aided
massive multiple-input multiple-output (MIMO) systems with a two-timescale
design. First, the zero-forcing (ZF) detector is applied at the base station
(BS) based on instantaneous aggregated channel state information (CSI), which
is the superposition of the direct channel and the cascaded user-RIS-BS
channel. Then, by leveraging the channel statistical property, we derive the
closed-form ergodic achievable rate expression. Using a gradient ascent
method, we design the RIS passive beamforming relying only on the long-term
statistical CSI. We prove that the ergodic rate scales on the order of
$\mathcal{O}\left(\log_{2}\left(MN\right)\right)$, where $M$ and $N$ denote
the number of BS antennas and RIS elements, respectively. We also prove the
striking superiority of the considered RIS-aided system with ZF detectors over
the RIS-free systems and RIS-aided systems with maximum-ratio combining (MRC).

###### Index Terms:

Reconfigurable intelligent surface (RIS), intelligent reflecting surface
(IRS), statistical CSI, massive MIMO, ZF.

##  I Introduction

As an emerging technique, reconfigurable intelligent surface (RIS) has been
widely investigated and recognized as a cost-effective complement for future
systems[1, 2]. The RIS mainly consists of a large number of passive reflecting
elements that have low hardware cost and energy consumption. Besides, RIS can
help conventional systems overcome the blockage issue and assist the
transmission by creating high-quality transmission paths.

A well-acknowledged challenge for the RIS is that it may introduce heavy
channel estimation overhead. Fortunately, a novel and more practical
countermeasure, named as two-timescale beamforming design, has been proposed
and validated by some contributions[3, 4, 5, 6, 7, 8, 9]. On one hand, the
two-timescale scheme aims at designing the passive RIS beamforming based on
purely statistical channel state information (CSI), and this could effectively
reduce the overhead and the energy consumption in the operation of the RIS[5].
On the other hand, the two-timescale scheme designs the BS beamforming based
on the instantaneous aggregated channel, which is the superposition of the
direct channel and cascaded user-RIS-BS channel. As a result, the aggregated
channel has the same dimension as that in conventional systems. Thus, two-
timescale schemes possess the same channel estimation overhead as conventional
systems.

Inspired by the above benefits, the two-timescale design has been recently
exploited in RIS-aided massive MIMO systems[10]. It has demonstrated that by
integrating an RIS into conventional massive MIMO systems, the rate
performance can be significantly improved, especially when the original direct
links are weak due to the blockage. However, only the simple maximum-ratio
combining (MRC) detector was considered in [10], and it was revealed that the
achieved gains are limited by the multi-user interference. Therefore, it is
expected that the zero-forcing (ZF) detector, which can effectively mitigate
this interference, is more suitable for RIS-aided massive MIMO systems.
Different from [10], when using ZF detectors, matrix inversion operator
introduces the additional technical challenges of deriving the ergodic
capacity.

Against the above background, in this letter, we consider an RIS-aided massive
MIMO system with ZF detectors. We derive the closed-form expression for the
ergodic rate, which only relies on the long-term CSI. We then design the RIS
based on a gradient ascent algorithm. By analyzing the rate expression, we
find that it scales on the order of
$\mathcal{O}\left(\log_{2}\left(MN\right)\right)$, which indicates that the
RIS-aided massive MIMO system with ZF detectors has the ability to achieve
ultra-high system capacity.

##  II System Model

Figure 1: Uplink transmission in RIS-aided massive MIMO systems.

As illustrated in Fig. 1, the uplink transmission of a massive MIMO system is
considered. Different from conventional systems, an RIS is introduced and
equipped at the facade of a tall building close to $K$ single-antenna users to
improve their channel conditions. The considered model is especially suitable
for the scenario where some cell-edge users suffer from service degradation.
Denote the number of BS antennas and RIS elements as $M$ and $N$,
respectively, where $M>K$. Then, we can define the channel between the users
and the RIS, the channel between the RIS and the BS, the direct channel
between the users and the BS as $\mathbf{H}_{1}\in\mathbb{C}^{N\times K}$,
$\mathbf{H}_{2}\in\mathbb{C}^{M\times N}$ and
$\mathbf{D}\in\mathbb{C}^{M\times K}$, respectively.

Define the phase shift matrix of the RIS as $\mathbf{\Phi}={\rm
diag}\left\\{e^{j\theta_{1}},\ldots,e^{j\theta_{N}}\right\\}$ where
$\theta_{n}$ is the phase shift of the $n$-th RIS element. Herein, we can
express the cascaded user-RIS-BS channel as
$\mathbf{G}=\mathbf{H}_{2}\mathbf{\Phi}\mathbf{H}_{1}\in\mathbb{C}^{M\times
K}$, and then express the aggregated channel from users to the BS as
$\mathbf{Q}=\mathbf{G}+\mathbf{D}\in\mathbb{C}^{M\times K}$. It is worth
noting that this aggregated channel $\mathbf{Q}$ possesses the same dimension
as conventional massive MIMO systems.

Based on the above definitions, we next present the detailed channel model for
$\mathbf{Q}$. Firstly, considering that the direct links may be easily
blocked[3], we adopt the Rayleigh channel model for $\mathbf{D}$ as follows

$\displaystyle\mathbf{D}=\tilde{\mathbf{D}}\bm{\Omega}_{d}^{1/2},$ (1)

where
$\bm{\Omega}_{d}=\operatorname{diag}\left\\{\gamma_{1},\ldots,\gamma_{K}\right\\}$,
and $\gamma_{k}$ denotes the distance-dependent path-loss factor. Each element
of matrix $\tilde{\mathbf{D}}\in\mathbb{C}^{M\times K}$ is independent and
identically distributed (i.i.d.) complex Guassian random variables, whose mean
is zero and variance is unit.

Next, since we consider that the RIS is deployed close to the users, and
according to the fact that the RIS is often installed above the ground, we
assume that the user-RIS channels have purely line-of-sight (LoS) paths. Then,
we denote

$\displaystyle\mathbf{H}_{1}=\left[\sqrt{\alpha_{1}}\,\overline{\mathbf{h}}_{1},\ldots,\sqrt{\alpha_{K}}\,\overline{\mathbf{h}}_{K}\right],$
(2)

where ${\alpha_{k}}$ denotes the path loss. To specify the LoS channel
$\overline{\mathbf{h}}_{k}$, we utilize the two-dimensional uniform squared
planar array (USPA) model[11]. Then, the array response vector for a
$\sqrt{X}\times\sqrt{X}$ USPA can be expressed as follows

$\displaystyle{{\bf
a}_{X}}\left({\vartheta^{a},\vartheta^{e}}\right)=\left[1,...,{e^{j2\pi\frac{d}{\lambda}\left({x\sin\vartheta^{e}\sin\vartheta^{a}+y\cos\vartheta^{e}}\right)}},\right.$
$\displaystyle\qquad\;\;\left.\ldots,{e^{j2\pi\frac{d}{\lambda}\left({\left({\sqrt{X}-1}\right)\sin\vartheta^{e}\sin\vartheta^{a}+\left({\sqrt{X}-1}\right)\cos\vartheta^{e}}\right)}}\right]^{T},$
(3)

where $0\leq x,y\leq\sqrt{X}-1$ are element indices in the two-dimensional
planar array, $d$ and $\lambda$ denote the element spacing and wavelength,
$\vartheta^{a}$ and $\vartheta^{e}$ are azimuth and elevation angles in the
propagation path, respectively. Therefore, denoting by the azimuth and
elevation AoA of user $k$ as $\varphi_{kr}^{a}$ and $\varphi_{kr}^{e}$, we can
now express $\overline{\mathbf{h}}_{k}={{\bf
a}_{N}}\left({\varphi_{kr}^{a},\varphi_{kr}^{e}}\right)$.

Since the RIS is placed near the users, the distance between the RIS and the
BS could be a bit large. Even though both the RIS and the BS have certain
heights, it is still not guaranteed that the RIS-BS channel is purely LoS. As
a result, the Rician model is suitable for the considered RIS-BS channel.
Besides, by adjusting the value of Rician factors, we can study the impacts of
scatterers in RIS-aided systems. This feature is important, since many works
have proven that rich scattering environment is beneficial in conventional
massive MIMO systems[12], while the corresponding impact in RIS-aided massive
MIMO systems with ZF detector is still unknown. Thus, we define

$\displaystyle\begin{array}[]{l}\mathbf{H}_{2}=\sqrt{\frac{\beta\delta}{\delta+1}}\overline{\mathbf{H}}_{2}+\sqrt{\frac{\beta}{\delta+1}}\tilde{\mathbf{H}}_{2},\end{array}$
(5)

where $\beta$ is the path loss, $\delta$ represents the Rician factor. Note
that the Rician factor, varying from $0$ to $\infty$, characterizes the
strength ratio between LoS and non-LoS (NLoS) paths. The NLoS path
$\tilde{\mathbf{H}}_{2}$ contains i.i.d. complex Gaussian random variables
with zero mean and unit variance. Recalling USPA model (II), the LoS path,
$\overline{\mathbf{H}}_{2}$, is written as

$\displaystyle\overline{\mathbf{H}}_{2}=\mathbf{a}_{M}\left(\phi_{r}^{a},\phi_{r}^{e}\right)\mathbf{a}_{N}^{H}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right)\triangleq\mathbf{a}_{M}\mathbf{a}_{N}^{H},$
(6)

where a notational simplification
$\overline{\mathbf{H}}_{2}\triangleq\mathbf{a}_{M}\mathbf{a}_{N}^{H}$ is
applied in the sequel of this paper. Note that the rank of matrix
$\overline{\mathbf{H}}_{2}$ is one. This means that when $\delta\to\infty$,
the cascaded channel $\mathbf{G}$ may become rank-deficient, and then the
achievable spatial multiplexing gains may degrade.

We can now express the $M\times 1$ received signal vector at the BS as

$\displaystyle\mathbf{y}=\sqrt{p}\mathbf{Q}\mathbf{x}+\mathbf{n}=\sqrt{p}\left(\mathbf{H}_{2}\mathbf{\Phi}\mathbf{H}_{1}+\mathbf{D}\right)\mathbf{x}+\mathbf{n},$
(7)

where
$\mathbf{x}=\left[x_{1},\ldots,x_{K}\right]^{T}\sim\mathcal{CN}\left(\mathbf{0},\mathbf{I}_{K}\right)$
includes the transmit symbols from $K$ users, and
$\mathbf{n}\sim\mathcal{CN}\left(\mathbf{0},\sigma^{2}\mathbf{I}_{M}\right)$
is the noise vector. For simplicity, we assume that all users transmit with
the same power $p$.

To facilitate the analysis, we assume that in each channel coherence time, the
instantaneous aggregated channel $\mathbf{Q}$ is perfectly known at the BS,
which serves as an upper bound for practical systems. Based on the two-
timescale design framework, we need to design the BS beamforming based on
instantaneous aggregated CSI, i.e., $\mathbf{Q}$. Thus, the ZF detector at the
BS is designed as
$\mathbf{A}=\mathbf{Q}\left(\mathbf{Q}^{H}\mathbf{Q}\right)^{-1}$, which
results in $\mathbf{A}^{H}\mathbf{Q}=\mathbf{I}_{K}$. Thus, the detected
symbol vector is given by

$\displaystyle\begin{aligned}
\mathbf{r}&=\mathbf{A}^{H}\mathbf{y}=\sqrt{p}\mathbf{x}+\left(\mathbf{Q}^{H}\mathbf{Q}\right)^{-1}\mathbf{Q}^{H}\mathbf{n}.\end{aligned}$
(8)

For ZF, there is no multi-user interference. As a result, the signal-to-
interference-plus-noise ratio (SINR) reduces to the ratio of transmit power
and noise. Based on (8), the SINR of user $k$ is given by

$\displaystyle\mathrm{SINR}_{k}$
$\displaystyle=\frac{p\left[\mathbb{E}_{\mathbf{x}}\left\\{\mathbf{x}\mathbf{x}^{H}\right\\}\right]_{kk}}{[\mathbb{E}_{\mathbf{n}}\\{(\mathbf{Q}^{H}\mathbf{Q})^{-1}\mathbf{Q}^{H}\mathbf{n}\mathbf{n}^{H}\mathbf{Q}(\mathbf{Q}^{H}\mathbf{Q})^{-1}\\}]_{kk}}$
$\displaystyle=\frac{p}{\sigma^{2}\left[(\mathbf{Q}^{H}\mathbf{Q})^{-1}\right]_{kk}}.$
(9)

Then, the $k$-th user’s ergodic rate is lower bounded by

$\displaystyle R_{k}$
$\displaystyle=\mathbb{E}\left\\{\log_{2}\left(1+\mathrm{SINR}_{k}\right)\right\\}$
(10)
$\displaystyle\overset{(a)}{\geq}\log_{2}\left(1+\frac{p}{\sigma^{2}\mathbb{E}\left\\{\left[(\mathbf{Q}^{H}\mathbf{Q})^{-1}\right]_{kk}\right\\}}\right),$
(11)

where $(a)$ utilizes the Jensen’s inequality based on the fact that function
$f(x)\\!=\\!\log_{2}\left(1+\frac{1}{x}\right)$ is convex with respect to $x$.

##  III Rate Analysis and RIS Design

In this section, we first derive the closed-form expression for the rate
$R_{k}$, and then use the derived expression to propose a statistical CSI-
based RIS design.

To derive $R_{k}$, we need to compute
$\mathbb{E}\left\\{\left[(\mathbf{Q}^{H}\mathbf{Q})^{-1}\right]_{kk}\right\\}$.
To this end, we expand matrix $\mathbf{Q}^{H}$ as

$\displaystyle\begin{array}[]{l}\mathbf{Q}^{H}=\sqrt{\frac{\beta\delta}{\delta+1}}\mathbf{H}_{1}^{H}\mathbf{\Phi}^{H}\overline{\mathbf{H}}_{2}^{H}+\sqrt{\frac{\beta}{\delta+1}}\mathbf{H}_{1}^{H}\mathbf{\Phi}^{H}\tilde{\mathbf{H}}_{2}^{H}+\mathbf{\Omega}_{d}^{1/2}\tilde{\mathbf{D}}^{H}.\end{array}$
(13)

For ease of exposition, we define
$\mathbf{Q}^{H}\triangleq\left[\mathbf{q}_{1},\ldots,\mathbf{q}_{M}\right]$,
$\overline{\mathbf{H}}_{2}^{H}\triangleq\left[\overline{\mathbf{c}}_{1},\ldots,\overline{\mathbf{c}}_{M}\right]$,
$\tilde{\mathbf{H}}_{2}^{H}\triangleq\left[\tilde{\mathbf{c}}_{1},\ldots,\tilde{\mathbf{c}}_{M}\right]$,
and
$\tilde{\mathbf{D}}^{H}\triangleq\left[\tilde{\mathbf{d}}_{1},\ldots,\tilde{\mathbf{d}}_{M}\right]$.
Recalling that $\tilde{\mathbf{H}}_{2}$ and $\tilde{\mathbf{D}}$ are all
comprised of i.i.d. complex Gaussian variables, and $\tilde{\mathbf{H}}_{2}$
and $\tilde{\mathbf{D}}$ are mutual independent, we therefore have

$\displaystyle\tilde{\mathbf{c}}_{m}\sim\mathcal{CN}\left(\mathbf{0},\mathbf{I}_{N}\right),1\leq
m\leq M,$ (14)
$\displaystyle\tilde{\mathbf{d}}_{m}\sim\mathcal{CN}\left(\mathbf{0},\mathbf{I}_{K}\right),1\leq
m\leq M,$ (15)

where $\tilde{\mathbf{c}}_{i}$ and $\tilde{\mathbf{c}}_{j}$ are mutual
independent for $i\neq j$; $\tilde{\mathbf{d}}_{i}$ and
$\tilde{\mathbf{d}}_{j}$ are mutual independent, for $i\neq j$;
$\tilde{\mathbf{c}}_{i}$ and $\tilde{\mathbf{d}}_{j}$ are mutual independent
for all $i$ and $j$. Then, since the linear transformation for a standard
Gaussian random vector is still a Gaussian random vector[13], we obtain

$\displaystyle\begin{array}[]{l}\sqrt{\frac{\beta\delta}{\delta+1}}\mathbf{H}_{1}^{H}\mathbf{\Phi}^{H}\overline{\mathbf{c}}_{m}+\sqrt{\frac{\beta}{\delta+1}}\mathbf{H}_{1}^{H}\mathbf{\Phi}^{H}\tilde{\mathbf{c}}_{m}\\\
\sim\mathcal{CN}\left(\sqrt{\frac{\beta\delta}{\delta+1}}\mathbf{H}_{1}^{H}\mathbf{\Phi}^{H}\overline{\mathbf{c}}_{m},\frac{\beta}{\delta+1}\mathbf{H}_{1}^{H}\mathbf{H}_{1}\right),\forall
m\end{array}$ (18)

and
$\mathbf{\Omega}_{d}^{1/2}\tilde{\mathbf{d}}_{m}\sim\mathcal{CN}\left(\mathbf{0},\mathbf{\Omega}_{d}\right),\forall
m$, where the facts $\mathbf{\Phi}^{H}\mathbf{\Phi}=\mathbf{I}_{N}$ and
$\mathbf{\Omega}_{d}=\mathbf{\Omega}_{d}^{H}$ were used.

Next, taking into account that the sum of independent Gaussian vectors is
still Gaussian distributed[13, Theorem 1.2.14], we can obtain the statistics
of the $m$-th column of aggregated channel $\mathbf{Q}^{H}$ as follows

$\displaystyle\begin{array}[]{l}\mathbf{q}_{m}\sim\mathcal{C}\mathcal{N}\left(\sqrt{\frac{\beta\delta}{\delta+1}}\mathbf{H}_{1}^{H}\mathbf{\Phi}^{H}\overline{\mathbf{c}}_{m},\;\frac{\beta}{\delta+1}\mathbf{H}_{1}^{H}\mathbf{H}_{1}+\bm{\Omega}_{d}\right),\end{array}$
(20)

where $\mathbf{q}_{m}$, $1\leq m\leq M$, are mutual independent. Therefore,
$\mathrm{vec}\left(\mathbf{Q}^{H}\right)$ is a complex Gaussian vector with
the following mean and covariance matrices

$\displaystyle\begin{array}[]{l}\mathbb{E}\left\\{{\mathrm{vec}}\left(\mathbf{Q}^{H}\right)\right\\}={\mathrm{vec}}\left(\sqrt{\frac{\beta\delta}{\delta+1}}\mathbf{H}_{1}^{H}\mathbf{\Phi}^{H}\overline{\mathbf{H}}_{2}^{H}\right),\\\
\mathrm{Cov}\left\\{{\mathrm{vec}}\left(\mathbf{Q}^{H}\right)\right\\}=\mathbf{I}_{M}\otimes\left(\frac{\beta}{\delta+1}\mathbf{H}_{1}^{H}\mathbf{H}_{1}+\mathbf{\Omega}_{d}\right),\end{array}$
(23)

where $\mathrm{vec}$ and $\otimes$ denote the vectorization by column stacking
and Kronecker product, respectively.

Then, using the distribution of $\mathrm{vec}\left(\mathbf{Q}^{H}\right)$ and
following the notations in [12, Page 2], matrix $\mathbf{Q}$ is a complex
Gaussian distributed matrix, written as

$\displaystyle\begin{array}[]{l}\mathbf{Q}\sim\mathcal{C}\mathcal{N}\\!\left(\\!\sqrt{\frac{\beta\delta}{\delta+1}}\overline{\mathbf{H}}_{2}\mathbf{\Phi}\mathbf{H}_{1},\mathbf{I}_{M}\\!\otimes\\!\left(\frac{\beta}{\delta+1}\mathbf{H}_{1}^{H}\mathbf{H}_{1}\\!+\\!\mathbf{\Omega}_{d}\\!\right)\right)\\!.\\!\\!\\!\end{array}$
(25)

Therefore, the product $\mathbf{Q}^{H}\mathbf{Q}$ has a complex non-central
Wishart distribution[13, Definition 10.3.1], which can be repressed as

$\displaystyle\begin{array}[]{l}\mathbf{Q}^{H}\mathbf{Q}\sim\mathcal{W}_{K}\left(M,\frac{\beta}{\delta+1}\mathbf{H}_{1}^{H}\mathbf{H}_{1}+\mathbf{\Omega}_{d}\right.\\\
\;\;\quad,\left.(\frac{\beta}{\delta+1}\mathbf{H}_{1}^{H}\mathbf{H}_{1}+\mathbf{\Omega}_{d})^{-1}\\!\frac{\beta\delta}{\delta+1}\mathbf{H}_{1}^{H}\mathbf{\Phi}^{H}\overline{\mathbf{H}}_{2}^{H}\overline{\mathbf{H}}_{2}\mathbf{\Phi}\mathbf{H}_{1}\\!\right)\\!.\\!\\!\\!\end{array}$
(28)

Even though the non-central Wishart distribution (28) is accurate, its
statistics are very complicated, and then we cannot obtain a tractable
expression for insightful analysis. To facilitate the analysis, as in
contributions [14, 15, 16], we next approximate the non-central Wishart
distribution (28) as a central Wishart distribution with the same first-order
moment.

To begin with, the first-order moment for the considered non-central Wishart
distribution is[15, Eq. (45)]

$\displaystyle\begin{aligned}
&\mathbb{E}\left\\{\mathbf{Q}^{H}\mathbf{Q}\right\\}\\\
&=M\left(\frac{\beta}{\delta+1}\mathbf{H}_{1}^{H}\mathbf{H}_{1}+\mathbf{\Omega}_{d}\right)+\frac{\beta\delta}{\delta+1}\mathbf{H}_{1}^{H}\mathbf{\Phi}^{H}\overline{\mathbf{H}}_{2}^{H}\overline{\mathbf{H}}_{2}\mathbf{\Phi}\mathbf{H}_{1}\\\
&=M\left(\frac{\beta}{\delta+1}\mathbf{H}_{1}^{H}\mathbf{H}_{1}+\mathbf{\Omega}_{d}\right)+M\frac{\beta\delta}{\delta+1}\mathbf{H}_{1}^{H}\mathbf{\Phi}^{H}\mathbf{a}_{N}\mathbf{a}_{N}^{H}\mathbf{\Phi}\mathbf{H}_{1},\end{aligned}$
(29)

where the last equality is obtained by using (6) and
$\mathbf{a}_{M}^{H}\mathbf{a}_{M}=M$.

Therefore, a virtual central Wishart distribution with this moment is given
by[15, Sec. $\rm V$. A]

$\displaystyle\begin{array}[]{l}\mathbf{Q}^{H}\mathbf{Q}\\!\sim\\!\mathcal{W}_{K}\\!\\!\left(\\!M,\frac{\beta}{\delta+1}\mathbf{H}_{1}^{H}\mathbf{H}_{1}\\!+\\!\mathbf{\Omega}_{d}\\!+\\!\frac{\beta\delta}{\delta+1}\mathbf{H}_{1}^{H}\mathbf{\Phi}^{H}\mathbf{a}_{N}\mathbf{a}_{N}^{H}\mathbf{\Phi}\mathbf{H}_{1}\\!\right)\\!.\end{array}$
(31)

Based on the obtained complex central Wishart distribution (31), with the help
of [17, Table $\rm I$], we can obtain the expectation of the matrix inverse as
follows

$\displaystyle\mathbb{E}\\!\\!\left\\{\\!\\!\left(\mathbf{Q}^{H}\mathbf{Q}\right)^{-1}\\!\right\\}\\!\\!=\\!\frac{\\!\\!\left(\\!\frac{\beta}{\delta+1}\mathbf{H}_{1}^{H}\mathbf{H}_{1}\\!+\\!\mathbf{\Omega}_{d}\\!+\\!\frac{\beta\delta}{\delta+1}\mathbf{H}_{1}^{H}\mathbf{\Phi}^{H}\mathbf{a}_{N}\mathbf{a}_{N}^{H}\mathbf{\Phi}\mathbf{H}_{1}\\!\\!\right)^{-1}\\!}{M-K}\\!.$
(32)

Substituting (32) into (11), we obtain the lower bound of the ergodic rate of
user $k$ as follows

$\displaystyle
R_{k}\\!\geq\\!\log_{2}\\!\\!\left(\\!1\\!+\\!\frac{p\left(M-K\right)}{\sigma^{2}(\delta\\!+\\!1)\\!\left[(\mathbf{\Lambda}+\beta\delta\mathbf{H}_{1}^{H}\mathbf{\Phi}^{H}\mathbf{a}_{N}\mathbf{a}_{N}^{H}\mathbf{\Phi}\mathbf{H}_{1})^{-1}\right]_{kk}}\right)\\!,$
(33)

where
$\mathbf{\Lambda}=\beta\mathbf{H}_{1}^{H}\mathbf{H}_{1}+(\delta+1)\mathbf{\Omega}_{d}$.

Note that the derived expression, (33), depends only on the statistical CSI,
since the instantaneous CSI-related variables have been averaged out.
Therefore, based on the two-timescale design framework, we can use (33) to
design the phase shifts of the RIS only relying on statistical CSI. Since the
statistical CSI-based phase shifts design only needs to be done on a large
time-scale, the overhead can be effectively reduced. Besides, it is clear that
(33) is an increasing function of $p$ and the RIS-BS channel strength $\beta$,
but it is a decreasing function of noise power $\sigma^{2}$.

######  Corollary 1.

As $M\to\infty$, the rate can maintain non-zero when the power is scaled down
proportionally to $p=1/M$.

Proof: It can be proved by noticing that all the matrices in the denominator
of (33) do not depend on $M$. $\blacksquare$

######  Corollary 2.

When $M\to\infty$ or $p\to\infty$, RIS-aided massive MIMO systems with ZF
detectors perform much better than that with MRC detectors.

Proof: Based on (33), when $M\to\infty$ or $p\to\infty$, we have
$R_{k}\to\infty$, while the rate in RIS-aided massive MIMO systems with MRC
detectors is still bounded due to the multi-user interference, as proved in
[10, Eq. (7)]. $\blacksquare$

######  Corollary 3.

When $\beta=0$, i.e., without the existence of the RIS, the rate of user $k$
reduces to

$\displaystyle
R_{k}\geq\log_{2}\left(1+{p{\left(M-K\right)\gamma_{k}}}/{\sigma^{2}}\right),$
(34)

which is the same rate as [18, Eq. (20)], and scales on the order of
$\mathcal{O}\left(\log_{2}\left(M\right)\right)$.

######  Corollary 4.

The ergodic rate of user $k$ in (33) is further lower bounded by

$\displaystyle R_{k}$
$\displaystyle\geq\log_{2}\left(1+\frac{p\left(M-K\right)}{\sigma^{2}(\delta+1)\left[\mathbf{\Lambda}^{-1}\right]_{kk}}\right)$
(35)
$\displaystyle\approx\log_{2}\left(1+\frac{p\left(M-K\right)}{\sigma^{2}}\Big{(}\frac{N\alpha_{k}\beta}{\delta+1}+\gamma_{k}\Big{)}\right),\text{
as }N\to\infty.$ (36)

which scales on the order of
$\mathcal{O}\left(\log_{2}\left(MN\right)\right)$.

Proof: Since we consider the existence of direct links, there is
$\mathbf{\Omega}_{d}\succ\mathbf{0}$ and then
$\mathbf{\Lambda}\succ\mathbf{0}$ and $\mathbf{\Lambda}^{-1}\succ\mathbf{0}$.
Besides, we have $\mathbf{\Lambda}^{H}=\mathbf{\Lambda}$. Based on the
Woodbury’s identity, we have

$\displaystyle{\left[\left(\mathbf{\Lambda}+\beta\delta\mathbf{H}_{1}^{H}\mathbf{\Phi}^{H}\mathbf{a}_{N}\mathbf{a}_{N}^{H}\mathbf{\Phi}\mathbf{H}_{1}\right)^{-1}\right]_{kk}}\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad$
$\displaystyle\\!=\\!\\!\left[\mathbf{\Lambda}^{-1}\right]_{kk}-\frac{\beta\delta\left[\mathbf{\Lambda}^{-1}\mathbf{H}_{1}^{H}\mathbf{\Phi}^{H}\mathbf{a}_{N}\mathbf{a}_{N}^{H}\mathbf{\Phi}\mathbf{H}_{1}\mathbf{\Lambda}^{-1}\right]_{kk}}{1+\beta\delta\mathbf{a}_{N}^{H}\mathbf{\Phi}\mathbf{H}_{1}\mathbf{\Lambda}^{-1}\mathbf{H}_{1}^{H}\mathbf{\Phi}^{H}\mathbf{a}_{N}}$
(37)
$\displaystyle\\!=\\!\\!\left[\mathbf{\Lambda}^{-1}\\!\right]_{kk}\\!\\!-\\!\frac{\beta\delta\left|\left[\mathbf{\Lambda}^{-1}\mathbf{H}_{1}^{H}\mathbf{\Phi}^{H}\mathbf{a}_{N}\right]_{k}\right|^{2}}{1\\!+\\!\beta\delta\mathbf{a}_{N}^{H}\mathbf{\Phi}\mathbf{H}_{1}\mathbf{\Lambda}^{-1}\mathbf{H}_{1}^{H}\mathbf{\Phi}^{H}\mathbf{a}_{N}}\\!\leq\\!\\!\left[\mathbf{\Lambda}^{-1}\right]_{kk}\\!.$
(38)

Substituting (38) into (33), we arrive at (35). Note that the $k$-th diagonal
element of $\mathbf{H}_{1}^{H}\mathbf{H}_{1}$ equals ${\alpha_{k}}N$, while
the non-diagonal elements are not proportional to $N$. When $N\to\infty$, we
can approximate $\mathbf{H}_{1}^{H}\mathbf{H}_{1}$ as
$N\mathrm{diag}\left\\{\alpha_{1},\ldots,\alpha_{K}\right\\}$, which results
in the approximation in (36). $\blacksquare$

Corollary 4 reveals a very promising capacity gain. It is well-known that the
ergodic rate of RIS-aided systems scales as
$\mathcal{O}\left(\log_{2}\left(MN^{2}\right)\right)$ in the single-user
scenario[3]. Here, we prove that by using ZF detectors in setup of multiple
users, the rate of each user could still scale as
$\mathcal{O}\left(\log_{2}\left(MN\right)\right)$, which demonstrates that the
considered systems can achieve a promising sum user rate. Besides, comparing
(36) with (34), it is shown that the RIS-aided massive MIMO systems with ZF
detectors always outperform RIS-free massive MIMO systems. Meanwhile, it can
be observed that the lower bound (36) tends to (33) when $\delta\to 0$ and
tends to (34) when $\delta\to\infty$.

Next, we design the RIS phase shifts based on (33), which depends only on the
statistical CSI. The sum-rate maximization problem can be formulated as
follows

$\max\limits_{\mathbf{\Phi}}\;\;R^{s}=\sum\nolimits_{k=1}^{K}R_{k},\qquad\quad$
(39a) $\text{s.t. }\;\;\left|[\mathbf{\Phi}]_{nn}\right|=1,\;\;1\leq n\leq N.$
(39b)

Problem (39) is non-convex due to the non-convex unit modulus constraint.
However, we can still obtain a sub-optimal solution based on the gradient
ascent method. For tractability, we rewrite $\mathbf{\Phi}={\rm
diag}\left\\{\mathbf{v}^{H}\right\\}$, where
$\mathbf{v}=[e^{j\theta_{1}},\ldots,e^{j\theta_{N}}]^{H}$. Then, we provide
the gradient vector with respect to $\mathbf{v}$ in the following lemma.

######  Lemma 1.

The gradient of the objective function in (39) is

$\displaystyle\frac{\partial
R^{s}(\mathbf{v})}{\partial\mathbf{v}^{*}}=\sum_{k=1}^{K}\frac{\frac{\mathbf{B}\mathbf{v}}{\mathbf{v}^{H}\mathbf{A}_{k}\mathbf{v}}-\frac{\mathbf{v}^{H}\mathbf{B}\mathbf{v}\mathbf{A}_{k}\mathbf{v}}{\left(\mathbf{v}^{H}\mathbf{A}_{k}\mathbf{v}\right)^{2}}}{\ln(2)\left(1+\frac{\mathbf{v}^{H}\mathbf{B}\mathbf{v}}{\mathbf{v}^{H}\mathbf{A}_{k}\mathbf{v}}\right)},$
(40)

where

$\displaystyle\mathbf{A}_{k}=\frac{\sigma^{2}(\delta+1)}{p(M-K)}\left(\left[\mathbf{\Lambda}^{-1}\right]_{kk}\mathbf{B}-\beta\delta\mathbf{s}_{k}\mathbf{s}_{k}^{H}\right),$
(41)
$\displaystyle\mathbf{B}=\frac{1}{N}\mathbf{I}_{N}+\beta\delta\operatorname{diag}\left(\mathbf{a}_{N}^{H}\right)\mathbf{H}_{1}\mathbf{\Lambda}^{-1}\mathbf{H}_{1}^{H}\operatorname{diag}\left(\mathbf{a}_{N}\right),$
(42)

with
$\mathbf{s}_{k}^{H}\triangleq\left[\mathbf{\Lambda}^{-1}\mathbf{H}_{1}^{H}\operatorname{diag}\left(\mathbf{a}_{N}\right)\right]_{({k,:})}$
corresponds to the $k$-th row vector.

Proof: Substituting
$\mathbf{\Phi}^{H}\mathbf{a}_{N}=\operatorname{diag}\left(\mathbf{a}_{N}\right)\mathbf{v}$
into (III) and utilize (42), we have

$\displaystyle{\left[\left(\mathbf{\Lambda}+\beta\delta\mathbf{H}_{1}^{H}\mathbf{\Phi}^{H}\mathbf{a}_{N}\mathbf{a}_{N}^{H}\mathbf{\Phi}\mathbf{H}_{1}\right)^{-1}\right]_{kk}}$
$\displaystyle=\\!\left[\\!\mathbf{\Lambda}^{-1}\\!\right]_{kk}\\!-\\!\frac{\beta\delta\\!\left[\\!\mathbf{\Lambda}^{-1}\mathbf{H}_{1}^{H}\operatorname{diag}\\!\left(\mathbf{a}_{N}\right)\\!\mathbf{v}\mathbf{v}^{H}\operatorname{diag}\\!\left(\mathbf{a}_{N}^{H}\right)\mathbf{H}_{1}\mathbf{\Lambda}^{-1}\\!\right]_{kk}}{1+\beta\delta\mathbf{v}^{H}\operatorname{diag}\\!\left(\mathbf{a}_{N}^{H}\right)\mathbf{H}_{1}\mathbf{\Lambda}^{-1}\mathbf{H}_{1}^{H}\operatorname{diag}\left(\mathbf{a}_{N}\right)\mathbf{v}}$
$\displaystyle=\\!\frac{\mathbf{v}^{H}\left\\{\left[\mathbf{\Lambda}^{-1}\right]_{kk}\mathbf{B}-\beta\delta\mathbf{s}_{k}\mathbf{s}_{k}^{H}\right\\}\mathbf{v}}{\mathbf{v}^{H}\mathbf{B}\mathbf{v}}.$
(43)

Then, substituting (III) into (33), the sum rate can be rewritten as
$R^{s}=\sum_{k=1}^{K}\log_{2}\left(1+\frac{\mathbf{v}^{H}\mathbf{B}\mathbf{v}}{\mathbf{v}^{H}\mathbf{A}_{k}\mathbf{v}}\right)$.
Based on the chain rule, the gradient of a real function with respect to
complex vector variable is given by[19]

$\displaystyle\frac{\partial
R^{s}\left(\mathbf{v}\right)}{\partial\mathbf{v}^{*}}=\sum\limits_{k=1}^{K}\frac{1}{\ln(2)\left(1+\frac{\mathbf{v}^{H}\mathbf{B}\mathbf{v}}{\mathbf{v}^{H}\mathbf{A}_{k}\mathbf{v}}\right)}\frac{\partial\left(\frac{\mathbf{v}^{H}\mathbf{B}\mathbf{v}}{\mathbf{v}^{H}\mathbf{A}_{k}\mathbf{v}}\right)}{\partial\mathbf{v}^{*}}.$
(44)

Using
$\frac{\partial\left\\{\mathbf{v}^{H}\mathbf{B}\mathbf{v}\right\\}}{\partial\mathbf{v}^{*}}=\mathbf{B}\mathbf{v}$,
and
$\frac{\partial\left\\{\mathbf{v}^{H}\mathbf{A}_{k}\mathbf{v}\right\\}}{\partial\mathbf{v}^{*}}=\mathbf{A}_{k}\mathbf{v}$,
we have

$\displaystyle\begin{aligned}
&\frac{\partial\left(\frac{\mathbf{v}^{H}\mathbf{B}\mathbf{v}}{\mathbf{v}^{H}\mathbf{A}_{k}\mathbf{v}}\right)}{\partial\mathbf{v}^{*}}=\frac{\left\\{\frac{\partial\left(\mathbf{v}^{H}\mathbf{B}\mathbf{v}\right)}{\partial\mathbf{v}^{*}}\right\\}\mathbf{v}^{H}\mathbf{A}_{k}\mathbf{v}-\mathbf{v}^{H}\mathbf{B}\mathbf{v}\left\\{\frac{\partial\left(\mathbf{v}^{H}\mathbf{A}_{k}\mathbf{v}\right)}{\partial\mathbf{v}^{*}}\right\\}}{\left(\mathbf{v}^{H}\mathbf{A}_{k}\mathbf{v}\right)^{2}}\\\
&=\frac{\mathbf{Bv}\mathbf{v}^{H}\mathbf{A}_{k}\mathbf{v}-\mathbf{v}^{H}\mathbf{Bv}\mathbf{A}_{k}\mathbf{v}}{\left(\mathbf{v}^{H}\mathbf{A}_{k}\mathbf{v}\right)^{2}}=\frac{\mathbf{Bv}}{\mathbf{v}^{H}\mathbf{A}_{k}\mathbf{v}}-\frac{\mathbf{v}^{H}\mathbf{B}\mathbf{v}\mathbf{A}_{k}\mathbf{v}}{\left(\mathbf{v}^{H}\mathbf{A}_{k}\mathbf{v}\right)^{2}}.\end{aligned}$
(45)

Substituting (45) into (44) completes the proof. $\blacksquare$

Assume the variable in the $t$-th iteration is $\mathbf{v}^{t}$. Then, the
next variable $\mathbf{v}^{t+1}$ in the $(t+1)$-th iteration is given by

$\displaystyle\tilde{\mathbf{v}}^{t+1}=\mathbf{v}^{t}+\mu\left.\frac{\partial
R^{s}\left(\mathbf{v}\right)}{\partial\mathbf{v}^{*}}\right|_{\mathbf{v}=\mathbf{v}^{t}},$
(46)
$\displaystyle\mathbf{v}^{t+1}=\exp\left(j\arg\left(\tilde{\mathbf{v}}^{t+1}\right)\right),$
(47)

where $\mu$ is the step size which can be chosen by using backtracking line
search[4, 9]. (47) is a projection operation for meeting the unit modulus
constraint (39b).

##  IV Simulation Results

Unless otherwise stated, we consider $K=4$ users evenly located on the half-
circle centered of an RIS with a radius $d_{UI}=20$ m. The distance between
the RIS and the BS is $d_{IB}=700$ m. Using $d_{UI}$ and $d_{IB}$, the
distance between the users and the BS can be calculated by their geometric
relationship as [10]. Based on the distances, the path loss factors
$\alpha_{k}$, $\beta$ and $\gamma_{k}$ are calculated the same as [10].
Besides, we set $M=N=64$, $p=30$ dBm, $\delta=1$ and $\sigma^{2}=-104$ dBm.
The angles in the LoS channels are generated randomly from $[0,2\pi]$. The
Monte Carlo simulations are obtained based on (10) with $10^{4}$ times
average.

In Fig. 4, it can be observed that ZF-based RIS design outperforms the random
phase shifts-based design, the MRC-based design, and the RIS-free systems. The
superiority of ZF over MRC lies in the fact that RIS-aided systems suffer from
severe multi-user interference[10]. This is because users share the same RIS-
BS channel and then their cascaded channels are highly correlated. Therefore,
by effectively eliminating the interference, ZF can achieve a higher ergodic
rate than MRC. In addition, Fig. 4 validates the accuracy of the approximate
in (36). In Fig. 4, we verify the power scaling law as expected in Corollary
1. It again emphases the advantages of ZF-based RIS systems. Besides, all
numerical results show that the derived lower bound (33) is very tight with
the Monte Carlo simulations.

Figure 2: Rate versus the number of RIS
elements $N$.

Figure 3: Rate versus $M$, where power is scaled
down as $p=10/M$.

Figure 4: Rate versus the Rician factor $\delta$.

Fig. 4 plots the rate versus Rician factor $\delta$. It can be seen that when
$d_{IB}=700$ m, the rate with large $\delta$ performs worse, while a contrary
result is observed when $d_{IB}=300$ m. This is because when $d_{IB}$ is
large, the direct channel becomes weak, and then the cascaded channel
$\mathbf{G}$ becomes a dominant factor. In this case, when $\delta$ is large,
the channel $\mathbf{G}$ becomes rank-deficient, which degrades the rate.
However, when $d_{IB}$ is small, the direct links are strong. Since the direct
links have full-rank, the aggregate channel could always have full-rank. As
shown in (33), the variable $\mathbf{\Phi}$ can play more roles when $\delta$
is large, which results in large performance gains. Furthermore, Fig. 4
validates the tightness of (36) under all Rician factors.

##  V Conclusion

An RIS-aided massive MIMO system with ZF detectors was considered in this
paper. We first derived the closed-form ergodic rate expressions, whose lower
bound demonstrates that the rate can scale on the order of
$\mathcal{O}\left(\log_{2}\left(MN\right)\right)$. Then, using the gradient
ascent algorithm, we optimized the phase shifts of the RIS based on
statistical CSI. Finally, simulation results validated the correctness of our
analytical results.

## References

  * [1] C. Pan, H. Ren, K. Wang, M. Elkashlan, A. Nallanathan, J. Wang, and L. Hanzo, “Intelligent reflecting surface aided MIMO broadcasting for simultaneous wireless information and power transfer,” _IEEE J. Sel. Areas Commun._ , vol. 38, no. 8, pp. 1719–1734, Aug. 2020.
  * [2] C. Pan _et al._ , “Multicell MIMO communications relying on intelligent reflecting surfaces,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 8, pp. 5218–5233, Aug. 2020.
  * [3] Y. Han, W. Tang, S. Jin, C.-K. Wen, and X. Ma, “Large intelligent surface-assisted wireless communication exploiting statistical CSI,” _IEEE Trans. Veh. Technol._ , vol. 68, no. 8, pp. 8238–8242, Aug. 2019.
  * [4] Q. U. A. Nadeem, A. Kammoun, A. Chaaban, M. Debbah, and M. S. Alouini, “Asymptotic max-min SINR analysis of reconfigurable intelligent surface assisted MISO systems,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 12, pp. 7748–7764, Dec. 2020.
  * [5] M. M. Zhao, Q. Wu, M. J. Zhao, and R. Zhang, “Intelligent reflecting surface enhanced wireless networks: Two-timescale beamforming optimization,” _IEEE Trans. Wireless Commun._ , vol. 20, no. 1, pp. 2–17, Jan. 2021.
  * [6] Y. Jia, C. Ye, and Y. Cui, “Analysis and optimization of an intelligent reflecting surface-assisted system with interference,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 12, pp. 8068–8082, Dec. 2020.
  * [7] Y. Gao, J. Xu, W. Xu, D. W. K. Ng, and M. S. Alouini, “Distributed IRS with statistical passive beamforming for MISO communications,” _IEEE Wireless Commun. Lett._ , vol. 10, no. 2, pp. 221–225, Feb. 2021.
  * [8] A. Papazafeiropoulos, C. Pan, P. Kourtessis, S. Chatzinotas, and J. M. Senior, “Intelligent reflecting surface-assisted MU-MISO systems with imperfect hardware: Channel estimation and beamforming design,” _IEEE Trans. Wireless Commun._ , early access, 2021.
  * [9] A. Papazafeiropoulos _et al._ , “Asymptotic analysis of max-min weighted SINR for IRS-assisted MISO systems with hardware impairments,” _IEEE Wireless Commun. Lett._ , early access, 2021.
  * [10] K. Zhi, C. Pan _et al._ , “Statistical CSI-based design for reconfigurable intelligent surface-aided massive MIMO systems with direct links,” _IEEE Wireless Commun. Lett._ , vol. 10, no. 5, pp. 1128–1132, May 2021.
  * [11] S. Zhou, W. Xu _et al._ , “Spectral and energy efficiency of IRS-assisted MISO communication with hardware impairments,” _IEEE Wireless Commun. Lett._ , vol. 9, no. 9, pp. 1366–1369, Sep. 2020.
  * [12] S. Jin, X. Gao, and X. You, “On the ergodic capacity of rank-$1$ ricean-fading MIMO channels,” _IEEE Trans. Inf. Theory_ , vol. 53, no. 2, pp. 502–517, Feb, 2007.
  * [13] R. J. Muirhead, _Aspects of multivariate statistical theory_. John Wiley & Sons, 2009, vol. 197.
  * [14] C. Siriteanu, Y. Miyanaga, S. D. Blostein _et al._ , “MIMO zero-forcing detection analysis for correlated and estimated rician fading,” _IEEE Trans. Veh. Tech._ , vol. 61, no. 7, pp. 3087–3099, Sep. 2012.
  * [15] C. Siriteanu, A. Takemura, S. Kuriki _et al._ , “Schur complement based analysis of MIMO zero-forcing for rician fading,” _IEEE Trans. Wireless Commun._ , vol. 14, no. 4, pp. 1757–1771, Apr. 2015.
  * [16] Q. Zhang, S. Jin, K. Wong, H. Zhu, and M. Matthaiou, “Power scaling of uplink massive MIMO systems with arbitrary-rank channel means,” _IEEE J. Sel. Topics Signal Process._ , vol. 8, no. 5, pp. 966–981, Oct. 2014\.
  * [17] J. A. Tague and C. I. Caldwell, “Expectations of useful complex wishart forms,” _Multidimensional Syst. Signal Process._ , vol. 5, no. 3, pp. 263–279, Jul. 1994.
  * [18] H. Q. Ngo, E. G. Larsson, and T. L. Marzetta, “Energy and spectral efficiency of very large multiuser MIMO systems,” _IEEE Trans. Commun._ , vol. 61, no. 4, pp. 1436–1449, Apr. 2013.
  * [19] R. Hunger, “An introduction to complex differentials and complex differentiability,” 2007.