ALGORITHMS FOR REACHABILITY PROBLEMS ON STOCHASTIC MARKOV REWARD MODELS

IRFAN MUHAMMAD

A thesis submitted to the University of Birmingham for the degree of

DOCTOR OF PHILOSOPHY

School of Computer Science

College of Engineering and Physical Sciences

University of Birmingham

September 2020

Probabilistic model-checking is a field which seeks to automate the formal analysis of probabilistic models such as Markov chains. In this thesis, we study and develop the stochastic Markov reward model (sMRM) which extends the Markov chain with rewards as random variables. The model recently being introduced, does not have much in the way of techniques and algorithms for their analysis. The purpose of this study is to derive such algorithms that are both scalable and accurate.

Additionally, we derive the necessary theory for probabilistic model-checking of sMRMs against existing temporal logics such as PRCTL. We present the equations for computing first-passage reward densities, expected value problems, and other reachability problems. Our focus however is on finding strictly numerical solutions for first-passage reward densities. We solve for these by firstly adapting known direct linear algebra algorithms such as Gaussian elimination, and iterative methods such as the power method, Jacobi and Gauss-Seidel. We provide solutions for both discrete-reward sMRMs, where all rewards discrete (lattice) random variables. And also for continuous-reward sMRMs, where all rewards are strictly continuous random variables, but not necessarily having continuous probability density functions (pdfs). Our solutions involve the use of fast Fourier transform (FFT) for faster computation, and we adapted existing quadrature rules for convolution to gain more accurate solutions, rules such as the trapezoid rule, Simpson's rule or Romberg's method.

In the discrete-reward setting, existing solutions are either derived by hands, or a combination of graph-reduction algorithms and symbolically solving them via computer algebra systems. The symbolic approach is not scalable, and for this we present strictly numerical but relatively more scalable algorithms. We found each - direct and iterative - capable of solving problems with larger state spaces. The best performer was the power method, owed partially to its simplicity, leading to easier vectorization of its implementation. Whilst, the Gauss-Seidel method was shown to converge with fewer iterations, it was slower due to costs of deconvolution. The Gaussian Elimination algorithm performed poorly relative to these.

In the continuous-reward setting, existing solutions are adaptable from literature on semi-Markov processes. However, it appears that other algorithms should still be researched for the cases where rewards have discontinuous pdfs. The algorithm we have developed has the ability to resolve such a case, albeit the solution does not appear as scalable as the discrete-reward setting.

Much thanks goes to my supervisor Professor David Parker, who opened for me the doors to academia. Again thanks goes to him for his guidance, and his reviewing of my work. I also thank my RSMG group members, Professor Peter Hancox and Dr. Shan He, both of whom motivated me and trusted that I could accomplish this. And also thanks are to be given to Associate Professor Nick Hawes, and Professor Jeremy Wyatt, my previous co-supervisors who still benefited me significantly in the little of what they said.

I also would like to gratefully acknowledge my family and friends who helped me pull through during times difficult.

Finally, I am gratefully to DARPA for funding me this opportunity through the PRINCESS project (contract FA8750-16-C-0045), via the DARPA BRASS programme.

CHAPTER: INTRODUCTION
Probabilistic models have been utilized effectively in many areas of research and have numerous industrial applications. By abstracting real world systems, and making some sound assumptions of how they behave, we can use such models to represent a simplification of a natural phenomenon, or an artificial or engineered one. Doing so grants us a less expensive representation, one that we can study for insights as to how the true system actually behaves. The understanding we obtain is useful for many reasons. If the system being modelled is artificial, we can use the knowledge to determine whether the true system meets safety specifications. If it is a real world phenomena, we can use it to make predictions and forecasts, or perhaps to develop or prove scientific theories.

A probabilistic model is one which models systems predictable only through relative ratios of outcomes, i.e. systems best characterized as random. These models can then be used to determine the likelihood of certain critical events, such as failure rates of an industrial facility <cit.>, the likelihood of cancer re-emission or death <cit.> or the survivability of coronary patients <cit.>. Other applications of these models could be for weather-forecasting, speech-recognition, in the field of computational biology and robotics. A sufficiently flexible model that is well developed has the capacity to model complex real world problems. A flexible model is one that imposes fewer requirements on real-world systems, to be representable.

The field of probabilistic verification comprises various formal techniques to investigate behaviours of probabilistic models. These models range from simple structures such as Markov chains and probabilistic Petri nets to complicated models such as those written with probabilistic programming languages. A contribution to this field would be for example designing methods for investigating behaviours of models not previously analysable. Or alternatively, presenting new models. One new class of model in literature is the stochastic Markov reward model (sMRM) <cit.>, a model which captures the accumulation of rewards over a random (Markov) process. These rewards are random themselves, and are accumulated as independent increments. Its development is recent, and we do not find much work directly on them. In this thesis, we develop algorithms for investigating a set of questions concerning them. We consider both the discrete-reward and continuous-reward variants of these models, and we focus mostly on a class of problems known as reachability problems.

The techniques investigated in probabilistic verification are various, however one popular family of methods is probabilistic model checking. The method allows us to investigate both qualitative and quantitative  properties of particular probabilistic models, via a combination of specification languages of which are usually temporal logics - to write what we intend to investigate of a model, and a set of model checking algorithms that consists of logic translators, model transformers, and numerical algorithms such as Gaussian elimination, or the Gauss-Seidel methods and their variants. We focus on designing scalable practical numerical algorithms for sMRMs, implementing and experimenting with them.

§ MOTIVATION
Let us consider for the time being application in robotics. The ability to automate robots such that they succeed at their tasks without risk of failure opens the door for many an industrial application. For example, they may be used to investigate the sea floors for undetonated mines left from previous wars or used in nuclear facilities where humans are at risk from radiation poisoning. Or mundane tasks within factories or farming.

We motivate our research into a class of probabilistic models - stochastic Markov reward models - with the following problem:  A robot leaves a charging station to perform a set of tasks in its surrounding environment. The tasks consist of performing inspections on physical structures. The order in which it performs its tasks is determined by some planner; it visits a set of objects and inspects them before returning to recharge. Each movement or inspection consumes energy and time of the robot, both being random variables.

Concerning such a robot, we can ask some questions about its behaviour: 1) What is the probability that the robot runs out of energy during any of its missions? 2) Is it possible for the robot to complete the mission under some time bound $t$? 3) What is the likelihood that it fails to complete more than $x$ number of tasks due to running out of energy?

Problems like these can be modelled with discrete-time Markov chains (DTMCs), a model which allows us to approximate all possible behaviour of the robot with a finite set of configurations known as states, and which captures its behaviour evolution through time by using transitions between these states. For example, the mission can be modelled as a DTMC where the state of a robot is two variables: location, and whether the robot is inspecting or not. The state transitions represent movements of the robot from one location to another, and they also represent the starting, continuation and ending of an inspection task. To represent the energy and time costs, we have the choice of integrating these costs into the state space of the DTMC, that is the state of the robot will also include the amount of energy the robot currently has and how much time it has spent away from the recharging station. For complicated models however, this approach would blow up the memory requirements exponentially, thus leading us to considering better techniques.

Another solution is to use a related model known as the Markov reward model (MRM). In these models, rewards (or synonymously costs) can be represented in a manner whereby they are not integrated into the state space. Thus, the memory requirements are generally lesser than using DTMCs. An MRM is essentially a DTMC with a reward structure connected to it. Both the DTMC and the reward structure can be represented as individual matrices of equal size. We can represent time/energy costs more efficiently with these models. However, the MRM is more restrictive than a regular DTMC allowing only non-negative rewards/costs over transitions. For example with MRMs, we cannot transition to states which would reduce time spent on the mission or increase the robot's current energy amount (by recharging for example). This is not a problem we focus on at this time. Note that modelling both time/energy costs leads to a bi-dimensional (bi-variate) reward structure, thus two matrices (of equal size) are required for it instead of one.

A drawback with MRMs is that it only deals with deterministic rewards, i.e. it only allows us to model the energy between movements as fixed amounts. A more recent model has been introduced into literature at the time of this work <cit.>, which allows these rewards to be random variables. These models are referred to as stochastic Markov reward models (sMRM). It is this class of model that is the focus of this thesis  and what it builds upon.

However, a challenge with sMRMs is that multivariate random variable rewards (or random vector rewards) are not generally composed of mutually independent univariate random variables, unlike MRMs with deterministic rewards (i.e. degenerate random variables), where they are always independent of each other. Unfortunately, this blow up in complexity whilst important to us, we have no general solution for. The focus therefore is solely on univariate reward structures, or random vectors with mutually independent components.

One might propose another set of models instead of sMRMs, known as (transition-based) hidden Markov models (HMMs). This model extends DTMCs, by allowing there to be a further set of states called hidden states. Hence the original set of states of the underlying DTMC are known as observable states.  The complete state space of the model is the Cartesian product $\{\textit{hidden states} \times \textit{observable states}\}$. A required property of the model is that the likelihood of being in a hidden state at any particular time step is solely dependent on the transition made at the previous time step (between two observable states) by the HMM, and not on any other (previous) step.[In a state-based HMM, the likelihood of being in a hidden state at a particular time instance solely depends on the state the HMM is in at that time instance.] One might try to model question 1 above using an HMM by first using the observable states to represent the location of the robot and whether the robot is inspecting or not, just like with the DTMC. Then, introducing two hidden states for the model: the robot still has energy, and the robot has run out of energy. When the robot transits between observable states, e.g. moves location, it has a likelihood of running out of energy, or still having energy, i.e. it has a likelihood of ending up in one of the two hidden states. If the process is to be modelled as a HMM, the likelihood of being in a hidden state (e.g. having no energy) at any time must be modelled as a function of the most recent transition the robot took. However, this is theoretically unsound as the likelihood of having no energy is dependent on all previous transitions  and not just the most recent. As for a formal explanation, we delay this until after we have defined sMRMs properly (see Section <ref>).

If the robot does not yet have a plan of action (or mission), and we have a set of actions we can program the robot to do at any given state, if we were to ask a fourth question - for any state the robot is in, what is the best action the robot can take, in terms of preventing the robot from running out of energy and increasing the likelihood of completing a set of inspections? Then, a more general model or framework is required to answer this. One way to handle this is by extending previous models (DTMCs, MRMs, or sMRMs) with an action set which is used to annotate each state with actions we want the robot to be able to make. State transitions are now dependent on the action chosen by the robot at a given state. These modified models can then be studied to find the best action to make at any state given a particular problem. Extending DTMCs this way gives us the class of models known as Markov decision processes (MDPs). In the case of MRMs or sMRMs, we can call their extensions Markov reward decision processes (MRDPs) or stochastic MRDPs (sMRDPs). However, we have explained that the HMM is not a valid choice of model for our problem, and therefore its extension, called the partially observable Markov decision process (POMDP) cannot be considered.  However, an alternative manner of solving this is by using a black-box (or off-the-shelf) optimization algorithm that produces a plan of action. The black-box algorithm chooses a set of actions for the robot, one for each state. This set of actions then translates to a set of behaviours for the robot which can be modelled as previously via a DTMC, MRM or sMRM. The models are then investigated to see if for example the robot will run out of energy. If it is proven that the plan is sufficient (e.g. the robot will not run out of energy with high probability), the robot can be configured to run the plan. If not, the black-box algorithm will attempt to find a better plan through a new set of actions. Therefore an iteration between plan generation (via the black box algorithm), and plan verification (via the sMRM or other) is needed to determine a good set of actions that would be sufficient. In this thesis, we will not focus on planning problems. Rather the focus is strictly on sMRMs. However, as just explained, sMRMs can still be used as part of an optimization procedure to generate usable plans for a robot.

The robot problem serves to introduce the class of models being investigated in this thesis. It is not however our main focus, but it does clarify how these models can be used.

§ OBJECTIVES AND CONTRIBUTIONS

The stochastic Markov reward models (sMRMs) were recently introduced by <cit.>.  The authors would analyse their behaviour via simulation techniques. We seek to derive numerical algorithms for reachability problems defined over sMRMs, that are scalable and accurate. We avoid simulation techniques which offer generally statistical guarantees on the accuracy of the result, in favour of numerical approaches which can give formal guarantees. Additionally, simulation techniques are considered slow when accuracy is of concern.

Our contributions are as follows:

* We tie together model-checking with sMRMs using a temporal logic called PRCTL - probabilistic reward control tree logic. <cit.>. We lay a foundation for sMRM theory and present theoretical solutions for several problems.

* We present algorithms for our main problem, the computation of the the passage-time reward mass functions (described later), over the discrete-reward sMRMs that are fast and scalable. The algorithms involve adapting existing solutions for solving systems of linear equations: The Gaussian elimination algorithm, the power, Jacobi and Gauss Seidel methods.

* Algorithms for the continuous-reward sMRMs were also derived, that are fast but only slightly scalable. We discovered new quadrature rules for convolution that are amenable for use within the sMRMs that can provide more accurate answers. The solution we provide is perhaps a more general solution relative to some algorithms in literature which are adaptable for use, allowing us to solve problems that other algorithms cannot do without further development.

§ THESIS LAYOUT

The remainder of this thesis is organized as follows. The following chapter discusses relevant literature to the sMRM models. It discusses probabilistic model checking in the context of such models, as well as existing numerical algorithms utilizable for solving some problems of these models. Chapter 3 presents the necessary background for the thesis. It discusses Markov models, temporal logics and the algorithms for resolving properties over Markov models. Chapter 4 presents the theoretical foundations for sMRMs, introducing the system of convolution equations and provides proofs for the basis of our work. Chapter 5 begins with deriving an exact solution for discrete-reward sMRM problems via the Gaussian elimination algorithm. Chapter 6 introduces iterative methods for solving discrete-reward sMRMs via the power, Jacobi and Gauss-Seidel method. These methods are more scalable relative to the direct Gaussian elimination method. Chapter 7 concerns an attempt to resolve continuous-reward sMRMs. Finally, Chapter 8 concludes this thesis and presents direction for future work.

§ COMPUTER DETAILS

For the experiments found in this thesis, two different computers were used. Additionally, the algorithms we present are implemented in python. We made heavy use of the following python libraries: 1) numpy <cit.> for their numerical algorithms and matrix operations. 2) fftw, a python wrapper for the FFTW library <cit.> for FFT operations. 3) PaCal <cit.> - a probabilistic arithmetic calculator. This was used heavily as a benchmarking tool for the algorithms we develop.

The specs. of the two computers are as follows.

* Computer 1: A laptop with 7.7GB of RAM, and an Intel Core i7-5500U CPU (@ 2.40GHz x 4). Here, we are using the OPENBLAS <cit.> package as a back-end for numpy, but multi-threading was turned off.

* Computer 2: A desktop PC  with 31.9GB of RAM, and an AMD Ryzen 5 3600 6-Core CPU (@ 3.593GHz x 6).

Computer 1 is our default computer, and unless mentioned otherwise, can be assumed to be the computer used for a particular experiment in this thesis.

CHAPTER: RELATED WORK

§ INTRODUCTION
In this chapter, we discuss models related to the sMRM such as the Markov reward model (MRM) with deterministic rewards, and semi-Markov processes (SMP). We also delve into probabilistic model checking with temporal logics, an area concerned with solving general problems over probabilistic models such as DTMCs, MRMs, and others. Doing so will give us a reference as to what questions may be asked of sMRM behaviour. Then, we also discuss our focus and compare the algorithms derived here with possible algorithms that exist in literature.

§ PROBABILISTIC MODEL CHECKING
Probabilistic model checking is a wide field covering formal methods used to verify, prove or investigate the behaviour of probabilistic models. Models can be as specific as software code, or an abstraction of real-word systems. These models can be developed by hand, or generated automatically by software. When it is the latter, work to prevent memory overflow include using symbolic representations <cit.>, and SAT methods <cit.>. Once these models have been generated, there are a host of algorithms to analyse them, with results from fields including logic theory, automata theory, numerical methods and graph theory. The set of techniques we focus on are discussed in the book <cit.>, involving state spaces and temporal logics (or more generally specification languages).

Probabilistic model checking with temporal logics consists of a combination of three separate components, the first two being: A probabilistic model, and a specification language (e.g. a temporal logic). The model is used to represent the phenomena or system at hand, whilst the language allows us to write questions (synonymous to statements or properties or problems) we want to investigate concerning the model. We have seen briefly what the models and questions could be in the introduction. For those questions that are decidable and solutions exist for them which are unique, the algorithms that resolve them are an integral part of model checking, and form the third component.

A more formal understanding is that these logics present a way to express particular behaviours of the model by allowing us to write properties it may potentially exhibit. Then, the goal is to determine whether the model has these properties, or formally that the model satisfies these properties. The action of determining so, is termed model-checking.

We summarize below previously investigated problems on probabilistic specific models related to sMRMs, and also detail some existing temporal logics defined over them. The relevant logics are detailed in the following chapter and specific algorithms related to our work for proving satisfaction of their properties on models are presented.

§.§ DTMCs
The discrete time Markov chain (DTMC) forms the foundation for sMRMs and has been studied extensively. Problems for sMRMs that are independent of notions of rewards reduces to problems for DTMCs. Such problems can be partitioned into two: the long term behaviour of a Markov chain, and the short-term (or bounded) behaviour. The former includes investigating (i) expected first passage/arrival times - the expected number of steps to reach a state $j$ beginning from another state $i$, (ii) equilibrium/steady state distribution - the  probability distribution over being in each state of a Markov chain as time tends to infinity, (iii) and the mean recurrence time - the expected amount of steps it takes for a Markov chain to return to state $i$ beginning from the same state. See the textbook <cit.> for further details. For bounded behaviour, we have for example (iv) transient state probabilities - the probability distribution over being in a particular state at some given time $n$.

Probabilistic model-checking is a field which also studies DTMCs. The authors Hansson and Jonsson <cit.> introduced a temporal logic PCTL - probabilistic computation tree logic, which allows the expression of particular properties for DTMCs. All properties written in PCTL are entirely computable, and the algorithms developed to solve PCTL problems cover the entire class. The logic allows expression for a wide range of properties, which in turns allows us to solve problems including determining (i) first passage (or reachability) probabilities - the probability of first entering a state $j$ from each state of the process, (ii) step-bounded reachability probabilities - the probability of first entering a state $j$ from each state under $n$ steps, transient state probabilities, and (iii) repeated reachability and persistence probabilities - the probability that the DTMC repeatedly enters a set of states, and the probability of a DTMC transitioning only within a particular set of states respectively. See the textbook <cit.> for an exposition to the subject.

A property written in such specification languages is to be resolved over its respective DTMCs. If a DTMC satisfies a property, then we mean by this that the DTMC is guaranteed to exhibit such a behaviour. The algorithms involved for resolving properties over DTMCs can be categorized into three groups: 1) Translating the logical property into an equivalent form amenable to simpler computation. This would involve logic theory. 2) Transforming the DTMC with respect to the new property, preparing it for computation. This yields a simpler property as well. This is a combination of results and algorithms from graph theory and automata theory. 3) Solving the transformed DTMC with respect to the remaining property. As for algorithms for solving PCTL statements, they include those that range from being numerical or symbolic, global or local (on-the-fly), and deterministic or statistical. See for example <cit.>. If numerical approaches are used, then part of the solution may involve solving a system of linear equations, or repeated matrix-vector multiplications depending on the statement.

Whilst PCTL is our language of focus, there are other languages which allow resolving of other properties over DTMCs. For example there is PCTL* - probabilistic computation tree logic star by <cit.>. It is a logic that includes as subsets PCTL, and LTL - linear temporal logic introduced by <cit.>, which allows resolving of $\omega$-regular properties over DTMCs.

Additionally from model checking literature, is the problem of parametric model checking, the problem where models are not completely described and have parameters instead, and the goal is to determine if such models satisfy particular logical properties for a range of different values of these parameters. For DTMCs, see <cit.>. Additionally the problem of repairing Markov models has also been studied for controllable DTMCs <cit.>, where if such a DTMC does not satisfy a logical property, another closely related DTMC is sought to ensure the satisfaction of that property. PCTL has also been studied for DTMCs with continuous state spaces, where transition matrices are replaced with kernels instead <cit.>. Solutions may involve analogues of existing numerical solutions of finite-space DTMCs, see <cit.> as a guide. Available general surveys on probabilistic model checking include <cit.>.

§.§ Discrete-time MRMs and sMRMs
Markov reward models (MRMs) are essentially (discrete-time) Markov chains extended with a reward structure, allowing us to weight the occurrence of events of the process with a reward (or dually, cost) function. Then, we can ask not only for probabilities of events as with DTMCs but also the expected reward accumulated for these events. These rewards are either attached to states <cit.>, or attached to transitions of the process <cit.>, and the process accumulates these rewards when being in a state or transitioning between them respectively. If the rewards are deterministic, then we have regular MRMs. If they are random variables, then we have stochastic MRMs (sMRMs).

The study of Markov reward models similarly includes determining long-term and short-term behaviour. As for the former, it includes (i) the expected cumulated reward for reaching a set of states $B$ beginning in any particular state, (ii) and a conditional variant where we have the expected reward to reach $B$ under the condition that $B$ is eventually reached, called conditional expected cumulated reward. (iii) Also, we have long running averages of states - the expected cumulated reward earned when beginning in a state $s$, and (iv) quantile probabilities - the minimum reward bound such that the probability of reaching a set of states $B$ from each state $s$ whilst earning a reward less than the bound, is greater than a pre-specified probability $p$ <cit.>, and this is generalized for multivariate rewards by <cit.>.

One existing specification language for MRM is PRCTL - probabilistic reward computation tree logic, introduced by <cit.>. It is an extension of PCTL, but allows the computation of various expected value problems over MRMs that includes those above (except for quantile probabilities) but provides algorithms to other properties also. See the paper for details.

As for existing work directly on sMRMs, the variance of the cumulated reward has already been studied, as has the covariance of cumulated reward between two sMRMs with algorithms for both found in <cit.>. As for those who introduced the sMRM and gave it its name <cit.>, they presented Monte-Carlo algorithms for model-checking a class of (dependent) multivariate-reward  sMRMs with the logic PRCTL. They presented an example problem where they solved for the probability of reaching a set of states $B$ from a particular state $i$, with the mean cumulated reward being less than or equal to $r$. We however chose to focus on deriving numerical algorithms instead and focused on multivariate-rewards that are mutually independent as a first. Additionally, it would appear that expected value problems for sMRMs can draw ideas from regular MRMs, as we will show in Chapter 4, two expected value problems including the one above can use solutions similar to that for regular MRMs.

§.§ Continuous-time Markov models
Other related models include continuous-time Markov chains (CTMCs) and semi-Markov processes (SMPs).

As for CTMCs, this is an extension of DTMCs, where transition times are now random and distributed by exponential distributions. CTMCs have been studied to solve for the steady state distribution, and for the Kolmogorov forward and backward equations, with the latter being considered a major goal for CTMCs <cit.>. Roughly speaking, the forward equation describes the probability of being in a state $j$ at some time $t+h$, given that it was in state $i$ at time zero. The backwards equation is the reverse of this, it gives the probability of being in state $i$ at time zero, given that the process is in state $j$ at time $t+h$.  Within CTMCs are other popular models, such as birth-processes and birth-death processes, both of which forward and backward equations are studied for.

In probabilistic model checking, the temporal logic CSL - continuous stochastic logic has been introduced for CTMCs by <cit.>, which is a continuous-time variant of PCTL, hence similar properties can be model-checked. It is extended by <cit.> who also presented approximate model checking algorithms for the logic. The logic allows writing properties for (i) determining the probability of reaching a set of states (from a particular initial state) within a specific time interval, (ii)  finding the probability of remaining within a set of states within a specific time interval, and others.

If CTMCs are extended with (deterministic) rewards, then the new model is called a continuous MRM (CMRM). CMRMs are a subset of univariate sMRMs, where each reward is distributed as an exponential distribution, and is related to time. Here, CSRL - continuous stochastic reward logic  has been introduced for model-checking CMRMs by <cit.>. It includes as sub-logics, both CSL and CRL - continuous reward logic. CSRL allows writing properties similar to CSL, they can be used additionally for (i) determining the probability of reaching a set of states (from a particular initial state) within a specified time interval and within a specified reward interval, (ii) finding the probability of remaining within a set of states within a specific time interval and within a specified reward interval.

If we remove the model's restriction on transition time distributions being exponential distributions, we obtain the semi-Markov process (SMP) model. Hence, CTMCs are subsets of SMPs. Further, SMPs are syntactically univariate sMRMs (with non-negative rewards [i.e. time]) and are one of the closest models to it. Therefore much of their theory can be borrowed. Problems for SMPs include computing (i) first-passage time densities - the probability distribution function over the reward accumulated when starting from a state $s$ and reaching a set of states $B$ <cit.>, (ii) their moments <cit.>, (iii) the cumulative distributions for said densities and (iv) the hazard functions derived from these densities <cit.>. Additionally, the (v) mean recurrence time (as defined earlier for DTMCs) and
(vi) asymptotic state probabilities - the probability an SMP is in state $j$ given that it began in state $i$, when time tends to infinity  <cit.>.

Summarily, these various models presented help lay the foundation for sMRM theory. Not just their theory can be borrowed but also the practical algorithms developed for them for problem solving.

More generalized models can be studied that are beyond the scope of this thesis. For example by adding actions to DTMCs, we arrive at an important model, the Markov decision process (MDP). This can be annotated with rewards, yielding Markov reward decision processes (MRDPs or sMRDPs). SMPs have a parametrized variant where covariates have been introduced to them <cit.>. There is also work on stochastic hybrid systems <cit.>  and probabilistic programming languages <cit.>, both being quite general models. There are other probabilistic models which exist that have been extended with reward structures. For example, Markov automata <cit.> and stochastic Petri nets <cit.>.

§ ALGORITHMS FOR SMRMS
Our main focus is in resolving first passage reward densities or reward reachability densities - the distribution function over cumulated reward for first reaching a set of states $B$, having begun from any state of the sMRM. The reason for this focus is that reachability problems are cornerstones of probabilistic model checking, and solving them enables us to solve a large class of problems. An understanding of this can be grasped after reading the next two chapters.

We intend to determine these densities for two types of sMRMs: Continuous-reward and discrete-reward sMRMs. Their definitions will be given formally in Chapter 4, however the distinction is that continuous-reward sMRMs have their rewards all characterized by continuous random variables, whilst discrete-reward sMRMs have rewards characterized by discrete (lattice) random variables. These two classes are chosen as we have found them to have forms amenable to faster computations.

§.§ Continuous-Reward sMRM and SMPs

Syntactically, a semi-Markov process (SMP) is a stochastic Markov reward model, where the rewards are univariate random variables and represent time. Therefore many problems for semi-Markov processes are analogous to problems for sMRMs. Likewise are we able to adapt their solutions for use.

When solving for first passage time densities in an SMP or equivalently, first passage reward densities in an sMRM, we are generally confronted with a system of equations to solve <cit.>. Then, we have found that traditional algorithms apply: Direct numerical methods such as Gaussian elimination can be used, iterative numerical methods such as the power, Jacobi or Gauss-Seidel methods, or symbolic approaches for small problems.

In the original form, the system to solve is a system of convolution equations (described later in Chapter 4). This system is transformable into a particular system of linear equations, where each term is a function. A solution to the system is usually computed at samples of these functions. Thus, a numerical approach generally requires solving a system of equations for each sample. This space-complexity blow up means that symbolic approaches for small problems are useful, as solving the system avoids sampling the functions.

One algorithm for solving first passage time densities consists of three sub-algorithms: 1) Transforming the system into a set of linear equations. 2) Solving the system. 3) Inverse transforming the solution to obtain the densities.

Transforming the system can be done via the continuous Fourier transform, discrete Fourier, or Laplace. This is done exactly either algebraically by hand, or via a computer algebra system (CAS). If neither are possible, it can be approximated numerically. Note that the Laplace transform does not exist for every random variable, whilst the continuous Fourier does. The discrete Fourier is used as an approximation to the problem and is always done numerically.

Solving the system can be done as previously mentioned, either numerically, or algebraically (i.e. symbolically). However, combining numerical and symbolic approaches is also possible.

Inverse transforming the solution is less straightforward than the initial transformation of the system. This is because inversion problems can be ill-conditioned. The inverse Laplace transform is considered one of them <cit.>. However, for SMPs the inversion procedure appears unaffected since the transform is applied to a special class of functions - non-negative random variables, that are absolutely continuous with respect to the Lebesgue measure, this stated in <cit.>. However, in general sMRMs may contain random variables that are not non-negative. Nevertheless, there appears to be quite a few different algorithms for inverse transforms of the continuous Fourier, for example <cit.> and Laplace, where we have a multitude <cit.>.

As for existing solutions to SMPs, we find <cit.> using the Laplace transform, the power method, and two inversions transforms: Laplace-Euler <cit.> or Laplace-Laguerre <cit.>, as their three sub-algorithms. The paper details solving a system with more than a million states (under 10 minutes), hence showing the scalability of their approach. Another solution is by <cit.> who also uses the Laplace transform, a perhaps direct numerical approach for solving the system, and the Laplace-Euler inversion. Their problems included one with 9 states, and is mentioned to be resolvable under a second. An earlier paper by the same authors <cit.>, used the Laplace transform, a graph reduction algorithm which they have described for solving the system and experimented with two inversion algorithms: the Laplace-Euler and a saddlepoint approximation algorithm. For each of the inversion procedures, the paper describes a case where they would perform poorly. It would appear that the Laplace-Euler does not work for empirical distributions naively. They did propose a fix for this, however it is specifically for empirical Laplace transforms (ELTs) <cit.>[This paper may not be visible to the public, and we could not find access to it.]. Additionally from paper <cit.>, the Laplace transforms of the empirical distributions are numerically derived, and for a problem with around three states (and therefore up to nine transitions/pmfs) would take several minutes, which is slow. Another approach is <cit.> which uses the continuous Fourier transform, and experimented with two inverse transforms: The algorithm by <cit.>, and a saddlepoint approximation via <cit.>. The paper showed that the former inversion to be better relative to the saddlepoint approximation. The models they presented were small and algebraic solutions were presented for the first passage time densities.

Another algorithm from SMP literature available for finding first passage reward densities involves using the moments of these densities to infer the density itself. This is done using a vector of (non-negative integer moments) via the method of moments algorithm presented by <cit.>. It appears that whilst these algorithms may be useful as an approximation if using a few moments, may lead to representation explosion if high precision is required, this stated in  <cit.>. As a vector is required for each state of the system to store the moments, these vectors will grow as long as the required accuracy of the density has not been achieved.

Yet the most recent work <cit.>, that of which introduced the sMRM model into literature, resorted to sampling. This choice is perhaps due to sMRMs having generally $n$-dimensional reward random vectors, and therefore the complexity of such algorithms would be too high to resolve $n$-dimensional pdfs numerically. However, if the reward random vector consists of mutually independent random variables, then this reduces to solving 1-dimensional sMRMs $n$ times, which is significantly more tractable. This was the case with regular MRMs, where the reward random vector is solely composed of degenerate random variables (or constants), and hence are always independent of each other. Therefore independence between non-degenerate random variables is still a step up, and not sidewards. Sampling is generally regarded as slow, when high precision is required.

As for general n-dimensional reward random vectors. This is presumably future work. Whatever future technique is to be considered, one has to keep in mind the dangers of representation explosion or time-complexity growth when moving away from lattice representations.

In our thesis, we focus on univariate sMRMs. When rewards are continuous, we approximate the density via the discrete Fourier transform (DFT), which benefits from some quadrature rules we have developed. The system is solved via the power method only.  The transform is numerical and fast, and the inverse appears to work over discontinuous distributions, which was a problem with the popular Laplace-Euler technique for the Laplace inversion transform. In fact, <cit.> suggested to consider the DFT for the case of discontinuous distributions, or perhaps more generally, for problems where the Laplace inversion is generally ill-posed for. The power method was a good choice, as <cit.> showed how scalable it could be in solving problems with millions of states, although it is not sure if our algorithm scales as well with the DFT. Additionally the DFT requires just the pdf of these random variables, whether analytical or empirical. This is unlike the Laplace transform of a pdf which does not always exist. And the DFT is strictly numerical and can be computed quickly without requiring algebraic derivations or precise numerical integrations.

§.§ Discrete-reward sMRMs and SMPs
As for work on solving first passage time mass functions for discrete-time SMPs numerically, the only work we are familiar with is that by <cit.>, which used the continuous Fourier transform of a function when known, and the discrete Fourier when not known. The inversion was done by the inverse discrete Fourier transform. They presented a problem with 3 states, which was solvable algebraically by hand. In this case, the discrete Fourier transform is used as an approximation to the solution. We will reproduce their problem later in Section <ref>. The authors stated that they were deterred from the Laplace transform (perhaps due to the difficulty of inversion with discrete random variables).

For the case where rewards are discrete, we develop solutions for solving first passage reward mass functions, using the discrete Fourier transform and inverse transform. We are able to develop algorithms obtaining machine precision, via iterative methods such as the power method, Jacobi and Gauss-Seidel. Secondly, we also present a direct (exact) solution via the Gaussian elimination algorithm adapted for the system of convolution equations, which may prove useful for slowly converging problems that are not too big (due to its time complexity), if they occur. Thus, we move away from algebraic solutions found by hand (or computers) towards numerical algorithms.

CHAPTER: PRELIMINARIES

§ INTRODUCTION
We introduce here the theoretical foundations for probabilistic model checking for sMRMs with univariate rewards. To do so, we first define DTMCs, and use it to introduce the language PCTL - probabilistic tree control logic. Then MRMs are introduced as well as PRCTL - probabilistic reward control tree logic.

Furthermore, theory is presented on the topics of summations of random variables, and characteristic functions due to their relevance to sMRMs. Within the topic of summation of random variables, is the topic of convolution (as is seen later). We find it necessary to also introduce the inverse operation of (discrete) convolution; deconvolution, which will be used for subsequent chapters.

This chapter partially summarizes the book <cit.> on the topics of probabilistic model checking. Thus further results and more in-depth explanations can be found there. We will also adopt their notation quite considerably. Additionally, some main derivations presented in this thesis will be adaptations of a strategy found in <cit.>. The technique they presented is helpful in deriving solutions. Their notation will also be in this thesis.

§ DTMCS AND PCTL
Discrete time Markov chains are a class of probabilistic models that is used to represent a system consisting of a (finite) number of states, and where the system can only be in one state at any given time. As time passes, the system can transition between states randomly. Randomness is represented quantitatively as probabilities. Additionally, the model being Markovian assumes that if the system is in a particular state, the probability of transiting into another state must not be dependent on the previous states that the model was in.

Formally, A DTMC $\mathcal{M}$ is a tuple $(S, \textbf{P}, i_{init})$, where:

* $S$ is the state space, with the size being finite, i.e. $|S| < \infty$,

* $\textbf{P}: S \times S \rightarrow [0,1]$ is a map for the transition probabilities between any two states of the Markov chain. We set $\textbf{P}$ to be constant with respect to time, hence it is stationary.

* $i_{init}: S \rightarrow [0,1]$ is the initial distribution function of the chain, with $\sum_{s \in S}i_{init}(s) = 1$. This is to specify in which state the Markov chain (the system) is likely to have started in.

§.§.§ Sample space, events, and a probability measure

Let us define a sequence of states as a path, e.g. a path could be $s_0,s_1,s_2,\cdots s_n$ with each $s_i \in S$. A path can then be used to denote a possible outcome of the DTMC starting from state $s_0$, transiting consecutively to states $s_1,s_2,\cdots,s_{n-1}$, and ending in $s_n$. Then the set of all unique infinite length paths a DTMC exhibits defines a sample space. Any measurable subset of the sample space is generally called an event.

The manner in which probabilities are assigned to events must satisfy Kolmogorov's axioms, to be a valid classical probability theory. However, we can define the probability of a finite path $s_0,s_1,s_2,\cdots s_n$ without much complication,  written as $Pr(s_0,s_1,s_2,\cdots s_n)$, which is equal to the probability that the system began in $s_0$, and performs transitions until it reaches $s_n$. Since the process is Markovian, we have

= Pr(s_0)Pr(s_1|s_0)Pr(s_2|s_1,s_0)⋯·Pr(s_n|s_n-1⋯,s_0)

= Pr(s_0)Pr(s_1|s_0)Pr(s_2|s_1)Pr(s_3|s_2)⋯·Pr(s_n|s_n-1)

= i_init(s_0)P(s_0,s_1)P(s_1,s_2)⋯P(s_n-1,s_n)

Let $\pi$ be an arbitrary finite path, and $\Omega_{\pi}$ be the set of all infinite paths that begin with $\pi$. Then it must be the case that
$$Pr(\Omega_{\pi}) = Pr(\pi)$$

Therefore, we are now able to talk about the probability of certain events.

In literature, $\Omega_{\pi}$ is called the cylinder set of $\pi$.  Let $\Pi$ be the set of all finite paths possible through a DTMC. Then the set $\{\Omega_{\pi}| \pi \in \Pi \}$ characterises the set of all basis events, and the smallest $\sigma$-algebra derived from it represents the event space.

Then let $\mathscr{S}$ be the sample space, and $\mathscr{E}$ be the event space (or $\sigma$-algebra). We can define a complete probability space for DTMCs as $(\mathscr{S},\mathscr{E},Pr)$.

By definition, the empty path fragment $\pi = \{\} = \mathscr{S}$ (the sample space), occurs with probability 1, i.e. $Pr(\pi) = 1$.

§.§.§ PCTL notation for events and probability
We can define a grammar to specify events using probabilistic computation tree logic (PCTL) <cit.>. The logic introduces the following operators:$\ \Diamond, \Square, \Circle, \text{U}$, which translates loosely as eventually, always, next, and until.

Define $B,C$, both to be arbitrary subsets of $S$. Then for example:

* $\Diamond B$ - means the set of paths of a DTMC which eventually reaches $B$, i.e. all outcomes where a transition into a state of $B$ occurred,

* $C\ \text{U}\ B$ -  denotes the event, the collection of all outcomes, where the DTMC was in $C$ (until) before directly transiting into $B$. This excludes all outcomes not beginning in $C$, or not eventually transiting into $B$.

* $C\ \text{U}^{\leq n}\ B$ - denotes an event similar to above, however the set of outcomes is restricted only to those that transit into $B$ under $n$ steps, where $n \in \mathbb{N}$.

* $\Square B$ - the set of outcomes that never transition out of the states of $B$, i.e. always remaining within these states.

* $\Circle B$ - the set of outcomes that enter $\next B$ after the next (immediate) transition.

Let $s_i \in S$ be some state. Additionally from PCTL, the notation
$$s_i \vDash \Diamond B$$

means that the event $\Diamond B$ holds for $s_i$. In terms of events, this is conditioning the sample space to the set of paths beginning strictly from $s_i$, and restricting the set of paths to those satisfying $\Diamond B$. Thus in <cit.>, this is creating a new DTMC with $s_i$ as a deterministic initial state, which leads us to a new sample space $\mathscr{S}_s$ and then determining the event $\Diamond B$. An alternative representation could be
$$\Diamond B\ |\ (\text{initial state is $s_i$})$$

Finally, PCTL is specifically used to ask for the probability of these events, therefore we can ask for $Pr(s_i \vDash \Diamond B)$. Then

$$Pr(s_i \vDash \Diamond B) = Pr(\Diamond B | \ \text{initial state is $s_i$})$$

§.§.§ PCTL grammar
Valid statements of PCTL are determined by the following grammar. For the probability of the event $\psi$ beginning in a state $s$, i.e. $Pr(s \vDash \psi)$, $\psi$ can be defined from any of:

\begin{align*}
\psi \Coloneqq \Circle\Phi\ |\ \Phi_1\ \text{U}\ \Phi_2 \ |\  \Phi_1\ \text{U}^{\leq n}\ \Phi_2
\end{align*}

\begin{align*}
\Phi \Coloneqq \text{true}\ |\ C\ |\ \Phi_1 \cap \Phi_2\ |\ \neg\Phi\ |\ \mathbb{P}_{[a,b]}(\psi)
\end{align*}

where $C \subseteq S$ and
$$\mathbb{P}_{[a,b]}(\psi) = Pr(\psi) \in [a,b]$$
denoting whether or not the probability of the event $\psi$ is within some interval $[a,b] \subseteq [0,1]$ (excluding the empty-set $\{\}$), and where $n \in \mathbb{N}$. In the literature, $\Phi$ is termed a state formula, whilst $\psi$ is termed a path formula. A state $s$ satisfies a state formula $\Phi$ when $s \in \Phi$, and a path $\pi$ satisfies a path formula $\psi$ when $\pi \vDash \psi$, or  $\pi \in \psi$ if we interpret $\psi$ as an event in a probability space.

The grammar is written quite concisely, hiding the inclusion of the other operators previously mentioned. For example, the remaining logical operators: OR $\vee$, and implication $\rightarrow$, can be derived using $\neg,\wedge$ alone. Additionally, the temporal-logic operators $\Square$, $\Diamond$, are both derivable using the until operator $\text{U}$. For example, (perhaps intuitively) we have $\Diamond B \equiv \text{true U}\ B$, and $\square B \equiv \neg \Diamond \neg B$. The operator $\text{U}^{\leq n}$ is called constrained-until or step-bounded-until. Other constrained operators can be defined, e.g. $\Diamond^{\leq n}$, or $\Square^{\leq n}$.

Reachability problems The event $\Diamond \Phi$ characterizes the class of reachability  problems. The phrase $\Diamond \Phi$ is synonymous to the set of paths where the set of states satisfying the state formula $\Phi$ is eventually reached. This is important since by logic theory, all unconstrained statements above reduces to mostly resolving strictly $\Diamond B_1$ or $\Circle B_2$ statements, i.e. the complexity of all non-constrained PCTL statements is not significantly much more than being able to resolve reachability statements and next statements.

Notice that the state formula $\Phi$ always reduces to a set of states $B \subseteq S$. In total therefore, there are only three main algorithms required for resolving PCTL statements. An algorithm for next statements, e.g. $\Circle B$, until statements $B_2\ \text{U}\ B_2$, and constrained-until statements e.g. $B_2\ \text{U}^{\leq n} \ B_2$, where $B,B_1,B_2$ are arbitrary subsets of $S$. Repeating from above, it is provable that until, and constrained-until statements for a DTMC $\mathcal{M}$ reduces to reachability and constrained reachability of a new DTMC $\mathcal{M}'$ <cit.>. Aside from the algorithm needed for this transformation, we mostly only need to solve statements of the form $\Circle B, \Diamond B_1$ and $\Diamond^{\leq r} B_2$.

Thesis focus In this thesis, we focus solely on solving unconstrained reachability problems. Hence, we will not delve much into $\Circle \Phi$ properties, nor constrained properties such as $\Phi_1\ \text{U}^{\leq n}\ \Phi_2$. For model checking algorithms for these, please refer to <cit.>. Reachability problems for sMRMs will be defined in the next chapter.  We detail below the means of computing reachability probabilities - $Pr(s \vDash \Diamond B)$.

§.§.§ Solution to reachability statements: $Pr(s \vDash \Diamond B)$

Given the event $\psi = \Diamond B$, then the set of all paths in $\psi$ can be uniquely defined as the set consisting of all cylinder sets of all unique finite paths $\pi= s_0,s_1,\cdots,s_n$, where $s_0 = s$, with $s_1,\cdots,s_{n-1} \notin B$ and $s_n \in B$. Then let this set of finite paths be $\Pi_\psi$, and the set $\Omega_{\psi} = \{\Omega_{\pi} | \pi \in \Pi_\psi \}$. In this case, the elements of $\Omega_{\pi}$ do not overlap , i.e. there are no two cylinder sets $\Omega_{\pi_1}, \Omega_{\pi_2}$ such that $\Omega_{\pi_1} \in \Omega_{\pi_2}$. Then,

$$Pr(\psi)  = \sum_{\Omega_{\pi} \in \Omega_{\psi} } Pr(\Omega_{\pi}) = \sum_{\pi \in \Pi_\psi} Pr(\pi)$$

Notation Before we proceed further, for the derivations to come, we will typically reserve letters $s,t$ to denote states, and symbols $\hat{\pi}, \hat{\psi}$ to denote paths. Additionally, when we write $\Pi.t$, we mean the set of paths of $\Pi$ that begin from state $t$, i.e. $\{\hat{\pi} \in \Pi | \ \hat{\pi}[0] = t \}$. We can concatenate states and paths via dot notation, for example $P(s.\hat{\pi})$ denotes the probability of a path event, one that begins with state $s$ and is followed by path $\hat{\pi}$.

Given a DTMC $\mathcal{M}$, let $B$ be a set of target states we are interested in. Given an initial state $s$, the probability of the event that $B$ is eventually reached is defined to be the collection of finite paths $\Pi_{\Diamond B}$, where each path $s_0,s_1,\cdots ,s_{n-1},s_n$ has $s_n \in B$, and $s_0,\cdots s_{n-1} \notin B$. Let $\Pi_{\Diamond B}.s$ denote the set of paths satisfying $\Diamond B$ beginning from $s$. Then,

Pr(s B)
= ∑_π∈Π_B.s Pr(π)

= ∑_t ∈S ∑_π∈Π_B.t Pr(s.π)

= ∑_t ∈S ∑_π∈Π_B.t P(s,t)Pr(π)

= ∑_t ∈S P(s,t) ∑_π∈Π_B.t  Pr(π)

= ∑_t ∈S P(s,t) Pr(t B)

We also have $Pr(s \vDash \Diamond B) = 1$ for all $s \in B$. Thus, the above is equivalent to

Pr(s B)
= ∑_q ∈S_? P(s,q) Pr(q B) + ∑_u ∈BP(s,u)

with $S_?$ being the set of states that reach $B$ with non-zero probability but are not within it.

One way to solve this for DTMCs with finite state spaces, is to represent the problem as a system of linear equations. Let $x_s \triangleq Pr(s \vDash \Diamond B)$ for all $s \in S$, then
$$x_s = \sum_{t \in S_?}\textbf{P}(s,t)x_t + \sum_{u \in B}\textbf{P}(s,u)$$

Then let x be the vector $(x)_{s \in S_?}$, and $\textbf{b} = (\sum_{u \in B}\textbf{P}(s,u))_{s \in S_?}$, and we define a matrix A where $A_{i,j} = \textbf{P}(i,j)$. Now we can write the system as
$$\textbf{x} = \textbf{Ax} + \textbf{b}$$
with the solution being
$$\textbf{x} = (\textbf{I}-\textbf{A})^{-1}\textbf{b}$$

Methods available for solving the system could be Gaussian-elimination, the power method, Jacobi or Gauss-Seidel methods.

§ MRMS AND PRCTL
Markov reward models (MRMs) extend DTMCs, by allowing transitions to be detailed with rewards or costs.[Rewards and costs are synonymous in this thesis. However, in the context of a problem, the use of one may be more suitable than the other.] Thus, this process not only transitions from state to state as time progresses (like a DTMC), but also accumulates reward along these transitions. MRMs can then be used to study the expected rewards of reaching a particular state beginning in another, for example.

An MRM is a tuple ($\mathcal{M}$, rew), where:

* $\mathcal{M}$ is a DTMC, i.e. $\mathcal{M} = (S, \textbf{P}, i_{init})$.

* rew: $S \times S \rightarrow \mathbb{R_{+}}$ is a reward function, that assigns a reward to each transition $s \rightarrow t$ in $\mathcal{M}$ for all $s,t \in S$. Note that the rewards are strictly non-negative for each transition.

Let $\pi = s_0,s_1,\cdots,s_n$ be a finite path, then the cumulated reward over $\pi$, denoted $Rew(\pi)$ is defined as
$$ Rew(\pi) = rew(s_0,s_1)+rew(s_1,s_2)+rew(s_2,s_3)+\cdots +rew(s_{n-1},s_n) $$

The probability of this path occurring, i.e. $Pr(\pi)$, is defined as it was for DTMCs. The same probability space as DTMCs can be used for MRMs.

§.§.§ Contingent expected rewards
One important measure for MRMs is to compute the expected (or mean) cumulated reward earned for a particular event, beginning in some state $s$. Let us denote this as $ExpRew(\psi)$, for the event $\psi$. Recall that the definition of the expected value of a discrete random variable defined over the real line is typically:

$$\mathbb{E}[X] = \sum_x x\cdot f(x)$$

In the context of MRMs, the expected (cumulated) reward with respect to the sample space $\mathscr{S}$ has to be defined in a manner that avoids $\infty$ being introduced unnecessarily. One consideration is that like DTMCs we can make use of cylinder sets. Let $\Omega_{\pi}$ be an event, where $\pi$ is a finite path. Then define the cumulated reward earned by the event $\Omega_{\pi}$ as $Rew(\pi)$. Thus all basis events $\Omega_{\pi}$ are well defined and less than $\infty$.

An event of interest can be decomposed into a set of non-overlapping cylinder sets (or basis events). Let this event be $\psi$, and $\Omega_{\psi}$ be the set of non-overlapping cylinder sets that describes all paths of the event. Define also $\Pi_\psi = \{\pi |\ \Omega_{\pi} \in \Omega_{\psi} \}$ . Then we define the expected reward contingent on this event as:

$$\mathbb{E}[\psi] = \sum_{\Omega_{\pi} \in \Omega_{\psi}} Rew(\pi)\cdot Pr(\Omega_{\pi})
=  \sum_{\pi \in \Pi_\psi} Rew(\pi)\cdot Pr(\pi) $$

However, $\psi$ being decomposable into cylinder sets is not enough to guarantee that $\mathbb{E}[\psi] \neq \infty$. This is since an event can characterize a single infinite path $\pi \in \mathscr{S}$, of which generally accumulates infinite reward. Thus, a possible consideration (in certain applications) is to assign zero-reward to such events.

Expected rewards of reachability problems The class of expected rewards contingent on the event $s \vDash \Diamond \Phi$ with $\Phi$ being a PCTL state formula and satisfying the condition that $Pr(s \vDash \Diamond \Phi) = 1$ for all $s \in S$, is known to always be finite. I.e. $\mathbb{E}[s \vDash \Diamond \Phi] < \infty$ for any $\Phi$. This class of reward problems is added to PCTL in an extension known as PRCTL - probabilistic reward computation tree logic <cit.>.

§.§.§ PRCTL grammar
The grammar for PRCTL extends PCTL with two terms:

* $\mathbb{E}(\Phi)$ - the expected reward for reaching a set of states determined by the state formula $\Phi$. Its definition is
$$\mathbb{E}_{[c,d]}(\Phi) = ExpRew(s \vDash \Diamond \Phi)$$

* $\Phi_1\ \text{U}_{\leq r}\ \Phi_2$ - reachability probabilities constrained on bounded rewards. Elaborating, this is the event where $\Phi_2$ is reached without accumulating more than $r \in \mathbb{N}$ in reward, with the DTMC being in states $\Phi_1$ up till that point. From this formula, we can derive reward-bounded reachability problems of the form $\Diamond_{\leq r} \Phi$ -  the event where $\Phi$ is reached without accumulating more than $r$ reward. These formulas are used to compute the probabilities of such events, i.e. $Pr(s \vDash \Diamond_{\leq r} \Phi)$.

The full PRCTL grammar is as follows:

\begin{align*}
\psi \Coloneqq \Circle\Phi\ |\ \Phi_1\ \text{U}\ \Phi_2 \ |\  \Phi_1\ \text{U}^{\leq n}\ \Phi_2  |\ \Phi_1\ \text{U}_{\leq r}\ \Phi_2 \
\end{align*}

\begin{align*}
\Phi \Coloneqq \text{true}\ |\ B\ |\ \Phi_1 \cap \Phi_2\ |\ \neg\Phi\ |\ \mathbb{P}_{[a,b]}(\psi)|\ \mathbb{E}_{[c,d]}(\Phi)
\end{align*}

where $[c,d] \subseteq [0,\infty)$, and $ [a,b] \subseteq [0,1]$ excluding the empty-set, and both $n,r \in \mathbb{N}$ and

$$ \mathbb{E}_{[c,d]}(\Phi) = ExpRew(s \vDash \Diamond B) \in [c,d]$$

The grammar above can be made to include the term $\Phi_1\ \text{U}_{[a,b]}\ \Phi_2$ - the event where $\Phi_2$ is reached in accumulating reward only between $a$ and $b$ (both in $\mathbb{N}$), but the DTMC also remaining within $\Phi_1$ until transitioning into $\Phi_2$. Then the computation $Pr(\Phi_1\ \text{U}_{[a,b]}\ \Phi_2)$ is equivalent to
$$Pr(\Phi_1\ \text{U}_{\leq b}\ \Phi_2) - Pr(\Phi_1\ \text{U}_{\leq a}\ \Phi_2)$$

Hence knowing how to solve for the reward-bounded reachability is all that is required.

Reachability problems Note that $ExpRew(s \vDash \Diamond \Phi)$  reduces to $ExpRew(s \vDash \Diamond C)$, for some $C \subseteq S$. Also, determining the probability of $Pr(s \vDash \Phi_1\ \text{U}_{\leq r}\ \Phi_2)$ is not much harder than determining $Pr(s \vDash \Diamond_{\leq r}\ B)$ for some $B \subseteq S$ <cit.>. We proceed to present algorithms for these two reachability reward problems.

§.§.§ Solution to expected reward of reachability statements: $\mathbb{E}(B) = ExpRew(s \vDash \Diamond B)$
Given an MRM $\mathcal{R}$, let $B \subseteq S$. Let us denote $\Pi_{\Diamond B}.s$ as the set of paths of the event $s \vDash \Diamond B$. Then

ExpRew(s B)
= ∑_π∈Π_B.s Rew(π) Pr(π)

= ∑_t ∈S ∑_π∈Π_B.t Rew(s.π) Pr(s.π)

= ∑_t ∈S ∑_π∈Π_B.t (rew(s,t) + Rew(π))·P(s,t)Pr(π)

= ∑_t ∈S P(s,t) ∑_π∈Π_B.t rew(s,t) Pr(π) + Rew(π) Pr(π)

= ∑_t ∈S P(s,t) (∑_π∈Π_B.t rew(s,t) Pr(π) + ∑_π∈Π_B.tRew(π) Pr(π) )

= ∑_t ∈S P(s,t) (rew(s,t) ∑_π∈Π_B.t  Pr(π) + ExpRew(t B) )

= ∑_t ∈S P(s,t)(rew(s,t)Pr(t B) + ExpRew(t B)  )

where we have also defined $ExpRew(t \vDash \Diamond B) = 0$, for all $t \in B$. I.e. the expected cumulated reward for states already in $B$ is zero. Note that the result above is strange in that typically we do not find the term $Pr(t \vDash \Diamond B)$ within the derivation. For example <cit.> both do not present such a term. However, in both, their derivations assumed that $Pr(s \vDash \Diamond B)=1$.

If $Pr(s \vDash \Diamond B)=1 $ for all $s \in S$, then the above is equivalent to:

ExpRew(s B)
= ∑_t ∈S P(s,t)(rew(s,t) + ExpRew(t B)  )

= ∑_t ∈S_? P(s,t)(rew(s,t) + ExpRew(t B)  ) +  ∑_u ∈B P(s,u)rew(s,u)

with $S_?$ being the set of states that can reach $B$ with non-zero probability, but not in $B$, i.e. $S \setminus B$ since $Pr(s \vDash \Diamond B)=1 $ for all $s \in S$.

Let $r_s \triangleq ExpRew(s \vDash \Diamond B)$, and $x_s \triangleq Pr(s \vDash \Diamond B)$ for all $s \in S$. Then the general form (<ref>) can be written as

$$r_s = \sum_{t \in S}\textbf{P}(s,t)(rew(s,t)\cdot x_t + r_t)$$

Then define

\begin{align*}
&\textbf{r} \triangleq      (r_s)_{s \in S_?} &\\
&\textbf{b} \triangleq (\sum_{t \in S_?}P(s,t)rew(s,t)x_s + \sum_{u \in B}\textbf{P}(s,u)rew(s,u)_{s \in S_?} &\\
&\textbf{A} \triangleq (A_{i,j})_{i,j \in S_?^2} \triangleq  (\textbf{P}(i,j)))_{i,j \in S_?^2}
\end{align*}

Now we can write a system of equations in the form $\textbf{r} = \textbf{Ar} + \textbf{b}$. And the solution is
$$\textbf{r} = (\textbf{I}-\textbf{A})^{-1}\textbf{b}$$

§.§.§ Solution to reward-bounded reachability probabilities: $Pr(s \vDash \Diamond_{\leq r} B)$

The approach in <cit.> solves this problem as follows:

Pr(s _≤r B)
= ∑_π∈Π_B.sPr(Rew ≤r| π)Pr(π)

= ∑_t ∈S∑_π∈Π_B.tPr(Rew ≤r| s.π)Pr(s.π)

= ∑_t ∈S∑_π∈Π_B.tPr(Rew ≤r| s.π)P(s,t)Pr(π)

= ∑_t ∈SP(s,t)∑_π∈Π_B.tPr(Rew ≤r - rew(s,t)| π)Pr(π)

= ∑_t ∈SP(s,t)Pr(t _≤(r - rew(s,t)) B)

Also, we have $Pr(s \vDash \Diamond_{\leq r} B)  = 1$ if $s \in B$. This is the same for $Pr(s \vDash \Diamond_{\leq q} B)$ for any $q \in \mathbb{N}$, when $s \in B$. This is since no reward is accumulated for these states, thus they satisfy the inequality $(Rew \leq r)$ immediately. If $s \notin B$, then $Pr(Rew \leq r - rew(s,t)| \pi) = 0$ if $rew(s,t) > r$, for any $t \in S$, since it is impossible to reach $B$ through $t$ satisfying $(Rew \leq r)$.  We also assume that $Pr(s \vDash \Diamond B) = 1$ for all $s \in S$.

Let $x_{s,r} \triangleq Pr(s \vDash \Diamond_{\leq r} B)$ for all $s \in S$. Then the results above can be written as

$$x_{s,r} = \sum_{t \in S}\textbf{P}(s,t)x_{t,(r-rew(s,t))}

Notice that $x_{s,p}$ can be computed independently from all $x_{t,q}$ where $q > p$, and $t \in S$. Thus, $x_{s,r}$ is solved by computing successively, the terms
$$(x_{s,0})_{s \in S_?}, (x_{s,1})_{s \in S_?},\cdots,(x_{s,r-1})_{s \in S_?},(x_{s,r})_{s \in S_?}$$
where each $(x_{s,p})_{s \in S_?}$ uses the previous terms as is seen in the equation above.

The above solution requires recursive computations, but yields a system of equations to solve when zero-rewards exist, i.e. when there exists $rew(s,t) = 0$, for some $s,t \in S^2$. For a state $t$, let $S_{t,0}$ be the set of states such that $rew(t,u) = 0$ for all $u \in S$.  Then, we can write the solution above as:

$$x_{s,r} = \sum_{t \in S\setminus S_{t,0}}\textbf{P}(s,t)x_{t,(r-rew(s,t))}
+ \sum_{q \in S_{t,0}}\textbf{P}(s,t)x_{t,r}

\begin{align*}
&\textbf{x}_p \triangleq      (x_{s,p})_{s,p \in S \times \{0,1,\cdots,r\}} &\\
&\textbf{b} \triangleq (\sum_{t \in S\setminus S_{t,0}}\textbf{P}(s,t)x_{t,(r-rew(s,t))})_{s \in S}     &\\
&\textbf{A} \triangleq (A_{i,j})_{i,j \in S^2} \triangleq  \Bigg(\begin{cases}
\textbf{P}(i,j)) & rew(i,j) = 0\\
0 & otherwise
\end{cases}
\Bigg)_{i,j \in S^2}
\end{align*}

Now we can write a system of equations in the form $\textbf{x}_p = \textbf{A}\textbf{x}_{p} + \textbf{b}$ with the solution being just $\textbf{x}_p = (\textbf{I}-\textbf{A})^{-1}\textbf{b}$. In this way we can compute
consecutively, to solve the problem.

§ SMPS
A semi-Markov process is a DTMC, except that transitions do not occur deterministically with respect to time, but rather by random. This randomness is characterized by transition time distributions.
Semi-Markov processes define a rich class of probabilistic models, including the DTMC and CTMC.

A semi-Markov process is a tuple ($\mathcal{M},\textbf{G}$) where:

* $\mathcal{M}$ is a DTMC, i.e. $\mathcal{M} = (S, \textbf{P}, i_{init})$.

* $\textbf{G}$: $(S \times S) \rightarrow f$ is a map between every state transition $t \in S \times S$ of the DTMC, and a probability distribution $f$ that determines the transition time distribution.

The difference between an SMP and a univariate sMRM is that $\textbf{G}$ is generalized to not be related to time. It is just a map over state transitions to reward distributions. Due to the closeness of these models, we choose not to present reachability problems and their theoretical solutions for SMPs here to prevent overlapping results. Instead, the details are presented only for sMRMs in the next chapter.

§ RANDOM VARIABLES AND CHARACTERISTIC FUNCTIONS

§.§ Representation of random variables

Given two random variables X, Y, they can be represented in a variety of ways, e.g. using their densities (pdf) $f_X(x), f_Y(y)$, or their cumulative distributions $F_X(x), F_Y(y)$, or their characteristic functions $\phi_X(\tau), \phi_Y(\tau)$. For any random variable $X$, its probability density function $f_X(x)$ can be transformed to a characteristic function $\phi_X(\tau)$, and then inverse transformed back into $f_X(x)$.

A characteristic function of a random variable $X$, written as $\phi_X(\tau)$ can be derived via the formula: $\phi_X(\tau) = E[e^{i \tau X}]$. Therefore if $X$ is continuous, then
$E[e^{i \tau X}]= \int_{- \infty}^{\infty}e^{\tau ix}f_X(x)dx$
and if $X$ is a discrete (lattice) random variable, for example if it has $\mathbb{N}$ as a support, then
$E[e^{i \tau X}]= \sum_{x=- \infty}^{\infty}e^{\tau ix}f_X(x)$

The characteristic functions can equivalently be derived by applying the correct Fourier transforms to either the pdf or pmf of $X$.

The continuous Fourier transform (FT) applied to a function $f$ is written as $\mathcal{F}\{f\}$ and defined as
$\mathcal{F}\{f\}(\tau) = \int_{- \infty}^{\infty}e^{\tau ix}f(x)dx$. The inverse continuous FT applied to a function $\phi$ is written as $\mathcal{F}^{-1}\{\phi\}$ with definition
$\mathcal{F}^{-1}\{\phi\}(x) = \int_{-\infty}^{\infty}\phi(\tau)e^{2\pi ix\tau}d\tau$.

The discrete-time FT applied to a function $g$ is written $\mathcal{F}_d\{g\}$ and defined as
$\mathcal{F}_d\{g\}(\tau) = \sum_{x=- \infty}^{\infty}e^{\tau ix}g(x)$. The inverse discrete-time FT applied to a function $\psi$ is written as $\mathcal{F}_d^{-1}\{\psi\}$ with definition
$\mathcal{F}_d^{-1}\{\psi\}(x) = \sum_{-\infty}^{\infty}\psi(\tau)e^{2\pi ix\tau}$.

Let $X$ be a random variable. If $X$ is continuous, then $\phi_X(\tau) = E[e^{i \tau X}] = \mathcal{F}\{f_X\}(\tau)$. If $X$ is a discrete lattice defined over $\mathbb{N}$, then $E[e^{i \tau X}] = \mathcal{F}_d\{f_X\}(\tau)$. Each of these characteristic functions can then be transformed back to the pdf or pmf with the respective inverse Fourier transform. Note that in the future we may drop the subscript $d$ from $F_d$, and therefore the type of Fourier transform is to be inferred from the context.

Note however that there are random variables that have characteristic functions but no analytical expressions are known for their probability density function (e.g. the Stable distribution). Nevertheless, every random variable  has a characteristic function. This is due to the proposition below.

For any random variable $X$, we have that $|\phi_X(\tau)| \leq 1$ for all $\tau$. Hence, the integral (or summation) above in $\mathcal{F}\{f_X\}(\tau)$ absolutely converges and always exists. See <cit.> for proof and details (of the continuous case).

A random variable $Z$ defined as a result of a sum of two other (independent) random variables $X,Y$, e.g $Z=X+Y$, can be represented as a density function, derived from a convolution operation on the other two respective probability density functions: $f_Z(z) = \int_{-\infty}^{\infty}f_X(y)f_Y(y-z)dy$. We will denote this operation as $f_Z(z) = f_X\ast f_Y$ where we use $\ast$ as the symbol for the linear convolution operator.

If we use characteristic functions instead, then the summation of these random variables can be performed via multiplication instead: $\phi_Z(\tau) = \phi_X(\tau)\phi_Y(\tau)$ .

§.§.§ Convolution and deconvolution

Convolution and deconvolution are two operators, each the inverse of the other. In this work we will denote convolution as $\ast$ and deconvolution as $\deconv$. The properties of convolution and deconvolution are similar to that of multiplication and division.  Let $f,g,h$ denote probability density functions or discrete lattice functions (i.e. arrays or vectors). Then, it is known that convolution is:

* Commutative: $f \ast g = g \ast f$

* Associative: $f \ast (g \ast h) = (f \ast g) \ast h$

* Distributive: $f \ast (g + h) = (f\ast g) + (f \ast h)$

Additionally, we have that $f \ast \delta_{x,0} = f$ and that $a(f \ast g) = (af)\ast g$, where $a$ is a constant, and $\delta_{x,0}$ is the Dirac delta if $f,g,h$ are continuous, or the Kronecker delta if discrete (lattices).

With respect to equations involving deconvolutions, then deconvolution has the properties:

* Right distributive: $((f + g) \deconv h) = (f \deconv h) + (g \deconv h)$, but not left distributive: $h \deconv (f+g) \neq  (f \deconv h) + (g \deconv h)$. This is like division.

* Yields the identity: $f \deconv f = \delta_{k,0}$, where $\delta_{k,0}$ is the Dirac or Kronecker delta, depending whether $f$ is continuous or discrete respectively.

Therefore for example we have that $((f + g) \pentagon h) - ((f + k) \pentagon h) \equiv (g-k) \pentagon h$.

Additionally, the following properties hold between convolution and deconvolution:

* $(f \ast g) \pentagon h = (f \pentagon h) \ast g$.

Proof:  Let $\mathcal{F}\{f\} = F$ be the Fourier transform of the pdf $f$. Then $\mathcal{F}\{(f \ast g) \pentagon h\} = \frac{F\cdot G}{H} = \frac{F}{G}H = \mathcal{F}\{(f \pentagon h) \ast g\}$. Applying the inverse Fourier transform both sides of this yields the result.
* $|f\ast g| \leq |f|\ast|g|$ and $|f \pentagon g| \geq |f|\pentagon|g|$.

Proof:  For convolution we have $|f \ast g|(t) = |\int_{x}dx\cdot f(t-x)g(x)| \leq \int_{x}dx\cdot|f(t-x)||g(x)| = (|f| \ast |g|)(t)$. For deconvolution, then since we have $|f| = |g \ast h| \leq |g| \ast |h|$, deconvolving both sides by $|h|$ yields, $|f| \deconv |h| \leq |g|$. Also we have $|g| = |f \deconv h|$. Therefore $|g| = |f \deconv h| \geq |f| \deconv |h|$.

§ SUMMARY
In this chapter we defined DTMCs and MRMs, and presented an introduction to model checking with temporal logics, more specifically PCTL and PRCTL. For particular problems, we explained their solutions, and by doing so we introduced much of the notation we will be using in this thesis.

More importantly, we highlighted that reachability problems are one of the main problems for model-checking DTMCs as the event $\Phi_1 \ \text{U}\ \Phi_2$ for a DTMC $\mathcal{M}_1$ can be reduced to some an event $\Diamond B$ in a transformed DTMC $\mathcal{M}_2$. This holds true for sMRMs too.

In the next chapter, we lay the theoretical foundations of sMRMs and define several problems of interest over them.

CHAPTER: STOCHASTIC MARKOV REWARD MODELS (SMRMS)

§ INTRODUCTION

In this chapter, we introduce the theory on stochastic MRMs (sMRMs), an extension of the traditional MRMs which allows rewards to be random variables (or random vectors). The theory on sMRM was introduced into probabilistic model checking recently by <cit.>. In the literature, sMRMs may be known previously as Markov processes with random rewards <cit.>, statistical flowgraphs <cit.>, simply rewards defined over Markov chains <cit.>. Whilst we present here derivations and theory of our own, literature on Markov chains with rewards and SMPs exists that share a similar theory, for example see the previously cited articles  and <cit.>.

The Markov chain in Figure <ref> captures the movement of a robot within a building. The robot begins in position $l_0$ of the building, and by moving probabilistically between places, it ends up eventually in $l_2$ or $l_4$. These states are its final destination, and we want to determine if it is capable of reaching such positions without running out of energy.

The energy cost of each transition is a random variable. Assuming we know these random variables in advance, we are able to incorporate them into our Markov chain to form stochastic Markov reward models.

[shorten >=2pt,node distance=3cm,auto]
[state,initial] (l0) at (0,0) $l_0$;
[state, right of=l0] (l1)  $l_1$;
[state,accepting, above of=l1]  (l2) $l_2$;
[state, right of=l1] (l3)  $l_3$;
[state,accepting, right of=l3] (l4)  $l_4$;

[->] (l0) edge node $0.2$ (l1);
[->] (l0) [bend left] edge node $0.8$ (l2);
[->] (l1) edge node $0.9$ (l3);
[->] (l1) [bend right] edge node $0.1$ (l2);
[->] (l2) edge [loop above] node $1$ (l2);
[->] (l3) edge[bend left] node $0.4$ (l1);
[->] (l3) edge node $0.6$ (l4);
[->] (l4) edge [loop above] node $1$ (l4);

A Markov chain with two absorbing states $l_2,l_4$.

A stochastic Markov reward model (sMRM) is a tuple (M, rew) where:

* $\mathcal{M}$ is a DTMC, a tuple $(S, \textbf{P}, i_{init})$.

* rew: $(S \times S) \rightarrow f$ is a map between every state transition $t \in S \times S$ of the DTMC, and a probability distribution $f$.

If the rewards are discrete random variables, then $f$ is a probability mass function (pmf), defined over lattices i.e. $f: \mathbb{L} \rightarrow [0,1]$, where $ \mathbb{L}$ is a countable lattice subset of  $\mathbb{R}$. $f$ is a pdf if the rewards are continuous instead.

$f$ may be a joint distribution in the case where the rewards are random vectors (or multivariate). Whilst we focus however on the univariate case, mutually independent piecewise multivariate rewards can be solved by considering each variable separately.

The sample space of the sMRM is represented as the product $\Pi \times \mathbb{R}$, where $\Pi$ is the set of all paths that the underlying DTMC can generate and $\mathbb{R}$ is the set of values for the reward. The basis events for $\Pi$ are the cylinder sets, whereas the basis events for rewards is $\mathcal{E}= \{(-\infty,x] | x \in \mathbb{R} \}$. The event space for rewards is the Borel $\sigma$-algebra, whilst the event space for paths is the smallest $\sigma$-algebra over the set of all cylinder sets. Let us denote them as $\Sigma_1$ and $\Sigma_2$ respectively. Then, we can define the product $\sigma$-algebra as $\Sigma_1 \times \Sigma_2$, or the product event space.

The probability space can be defined as  ($\mathscr{S} = \Pi \times \mathbb{R}, \mathscr{E} = \Sigma_1 \times \Sigma_2, Pr$).

For the majority of this work, we enforce that for all states in $s \in S$, the probability of reaching $B$ is one, i.e.
$$Pr(s \vDash \Diamond B) = 1$$
A manner to remove this restriction will be detailed later (see Sec. <ref>). Also, for the remainder of the thesis, we only solve for the case where each reward random variable $rew(s,t)$ is strictly non-negative, for all pairs $s,t \in S^2$.

§ REACHABILITY PROBLEMS

Given an sMRM $R_M$, there are several reachability problems we can try to solve:

* $Pr(s \vDash \Diamond_{=r} B)$ - The probability of accumulating $r$ reward, and eventually reaching $B \subseteq S$ starting from a state $s \in S$. However, we will use the notation $Pr(Rew = r\ \cap\ s \vDash \Diamond B)$ instead, and it will frequently be seen shortened as
$$Pr(r\ \cap\ s \vDash \Diamond B)$$

Let $\Pi$ be the set of paths starting from $s$ and ending in $B$, (i.e. those that satisfy $s \vDash \Diamond B$) then we have that

= Pr(r ∩ Π) = ∑_π̂ ∈ΠPr(r  ∩π̂) = ∑_π̂ ∈ΠPr(π̂)Pr(r | π̂)

where $Pr(r | \hat{\pi})$ is the probability density (or mass) function of the accumulated reward given a particular path $\hat{\pi} = s_0,\cdots,s_n$, i.e. it is the pdf of the random variable $Rew(r;\hat{\pi}) =  rew(s_0,s_1)+rew(s_1,s_2)+rew(s_2,s_3)+\cdots +rew(s_{n-1},s_n)$.

If $\sum_{\hat{\pi} \in \Pi}Pr(\hat{\pi}) = 1$, then $\mrmpdf$ is a convex combination of probability density functions; a mixture distribution. This is what we call the first-passage reward density (or mass function), the focus of our work. In the model checking literature, this may be better termed as reachability reward density. This may be a pmf or pdf depending on whether it is a discrete-reward sMRM or continuous-reward. The semantics is slightly different to PRCTL for MRMs as now we can interpret $Pr(Rew = r\ \cap\ s \vDash \Diamond B)$  as a function where $r$ is the variable of the function, i.e. we can denote it as $f_{s \vDash \Diamond B}(r)$ instead, or $f_s(r)$ if the property is obvious.

* $\mrmpdfcostbounded$ - The probability of reaching a set $B$ from $s$ with the cumulated reward being less than or equal to $r$. This is equivalent to reward-bounded reachability probability from PRCTL. Alternative notation would be $Pr(Rew \leq r\ \cap\ s \vDash \Diamond B)$. If interpreted as a function, this is the cumulative distribution function (cdf) of $Pr(Rew = r\ \cap\ s \vDash \Diamond B)$, and is denoted as $F_{s \vDash \Diamond B}(r)$ or $F_{s}(r)$.

* $Pr(s \vDash \Diamond_{\mathbb{E}[\leq r]} B)$ - the probability of reaching a set $B$ from $s$, with the mean cumulated reward being less than or equal to $r$. An example of this appeared in <cit.> for sMRMs.

* $\mrmexpected$ - The expected amount of reward accrued when starting from a state $s$ before reaching the target set $B$. This value is a scalar in $\mathbb{R}$. This is the expected reward for reachability from PRCTL.

* $Pr(s \vDash \Diamond_{\leq ?} B) > p$ - asks for the minimal reward bound $r$ such that the probability of reaching $B$ from $s$ is greater than $p$. This is known as a quantile query <cit.>.

* $Pr(s \vDash \Diamond^{\leq t}_{=r} B)$ - The constrained (or step-bounded) variant of $Pr(s \vDash \Diamond_{=r} B)$. Alternative notation for this is $Pr(r\ \cap\ s \vDash \Diamond^{\leq t} B)$, which is a pdf/pmf. A cdf variant of this can be constructed.

* $Pr(s \vDash \Circle_{=r} B)$ - The probability of reaching $B$ in the next step, beginning from $s$, whilst accumulating $r$ reward. The alternative notation for this is $Pr(r\ \cap\ s \vDash \Circle B)$, again a pdf/pmf, and we can have a cdf variant of this as well.

Hence, after learning how to solve for these properties, we can define PRCTL for sMRMs. The grammar is identical as for regular MRMs (see Section <ref>). The PCTL subset of PRCTL can be resolved with DTMC algorithms. Resolution of the expected value problems  (3, 4) for sMRMs is almost identical for those of MRMs. Only, (1, 2, 5, 6, 7) are differently computed, with (1) being able to borrow solution ideas from SMP theory.

The focus of this thesis will be on problem (1) above. As for (2), then this is just the numerical integration of (1).  As for (5), the quantile query is the smallest $r$ of the cdf $\mrmpdfcostbounded$, such that the probability is greater than $p$. Hence having solved (1), we have the ability to derive this. The solution for (6, 7) can be derived indirectly from the power method algorithm for solving (1). This is shown in Theorem <ref>.

We proceed to derive the solutions for the reachability problems of (3,4,2,5,1) in that respective order. We have left (1) for last as this is the main focus and will be elaborated.

§.§ $\mrmexpected$ properties
This is computed as

ExpRew(s B)
= 𝔼[∑_π∈Π_B.s Rew(π) Pr(π)]

= ∑_π∈Π_B.s 𝔼[Rew(π)] Pr(π)

= ∑_t ∈S ∑_π∈Π_B.t 𝔼[Rew(s.π)] Pr(s.π)

= ∑_t ∈S ∑_π∈Π_B.t 𝔼[(rew(s,t) + Rew(π))]·P(s,t)Pr(π)

= ∑_t ∈S P(s,t) ∑_π∈Π_B.t 𝔼[rew(s,t)] Pr(π) + 𝔼[Rew(π)] Pr(π)

= ∑_t ∈S P(s,t) (∑_π∈Π_B.t 𝔼[rew(s,t)] Pr(π) + ∑_π∈Π_B.t𝔼[Rew(π)] Pr(π) )

= ∑_t ∈S P(s,t) (𝔼[rew(s,t)] ∑_π∈Π_B.t  Pr(π) + ExpRew(t B) )

= ∑_t ∈S P(s,t)(𝔼[rew(s,t)]Pr(t B) + ExpRew(t B)  )

or equivalently since $Pr(s \vDash \Diamond B) = 1$ for all $s \in S$,

ExpRew(s B)
= ∑_t ∈S P(s,t)(𝔼[rew(s,t)] + ExpRew(t B)  )

which is just a system of linear equations as before.

Therefore, all that is needed is to be able to compute the expected value of all rewards, e.g. $\mathbb{E}[rew(s,t)]$, then the computation is identical to that of regular MRMs.

§.§ $Pr(s \vDash \Diamond_{\mathbb{E}[\leq r]} B)$ properties
It can be shown that

Pr(s _𝔼[≤r] B)
= ∑_t ∈SP(s,t)Pr(t _𝔼[≤(r - rew(s,t))] B)

using a similar derivation for MRMs described in Section (<ref>). Therefore the algorithm for solving this property can borrow the algorithm from Section (<ref>). Like the previous property, we are required first to compute the expected value of all rewards $\mathbb{E}[rew(s,t)]$.

§.§ $Pr(s \vDash \Diamond_{\leq r} B)$ properties
Assume for now we can already compute $Pr(r \cap s \vDash \Diamond B)$. Let us denote this pdf/pmf as $f_s(x)$ for short, and let $\Pi$ be the set of finite paths starting from $s$ and ending in $B$, (i.e. those that satisfy $s \vDash \Diamond B$). Then we have

Pr(s \vDash \Diamond_{\leq r} B) = \int_{x=0}^{r}f_{s}(x)dx

Pr(s _≤r B)    = ∑_π̂Pr(π̂ ∩(Rew ≤r))= ∑_π̂Pr(π̂)Pr(Rew ≤r | π̂)

= ∑_π̂Pr(π̂)∫_x=0^rf_Rew(π̂)(x)dx= ∫_x=0^r∑_π̂Pr(π̂)f_Rew(π̂)(x)dx

= ∫_x=0^rf_s(x)dx

Hence, if we have already computed $f_s(x)$, we can compute reward-bounded reachability probabilities by integration.

Note that if we compute the cdf

$$F_s(r) = \int_{x=0}^{r}f_{s}(x)dx$$

then we can obtain

$$Pr(s \vDash \Diamond_{[a,b]} B) = F_s(b) - F_s(a)$$

the probability that $B$ is reached from $s$ with reward accumulated only within the interval $[a,b]$.

§.§.§ Multivariate mutually-independent rewards
Consider an sMRM problem where the rewards are random vectors of dimension $n$. Then let Rew be the random vector, denoting the multivariate accumulated reward. Firstly, we have  $Pr(s \vDash \Diamond_{\leq r} B) = Pr(Rew \leq r \cap s \vDash \Diamond B)$. Then define

$$Pr(\textbf{Rew} \leq \vec{r}, s_0 \vDash B)$$
to be the probability that we can reach the set of states B from $s_0$ with reward accumulated under or equal to $\vec{r}$, where $\textbf{Rew} \leq \vec{r} = \cap_{i =1}^{n} Rew_i \leq r_i$. If the components of the random vectors are mutually independent of each other, then

$$Pr(\textbf{Rew} \leq \vec{r}, s_0 \vDash B)  = Pr(\cap_{i=1}^{n}Rew_i \leq r_i, s_0 \vDash B) = \prod_{i=1}^{n}Pr(Rew_i \leq r_i, s_0 \vDash B)$$

where for each $i$, $Pr(Rew_i \leq r_i, s_0 \vDash B)$ can be solved independently.

The solution above also gives us a means to compute quantile queries for multivariate independent rewards, which was studied for regular MRMs in <cit.>.

§.§ $Pr(s \vDash \Diamond_{\leq ?} B) > p$ quantile queries
The quantile query can be solved via interpolation of $F_s(x)$. Given $F_s(x)$, we are to find the smallest $r$ such that $F_s(r) > p$. A naive algorithm based on trial-and-error is to sample the distribution and compare it to $p$ and then move towards the direction of $p$ knowing that $F$ is monotonic. The algorithms for computing $F_s(x)$ in our thesis only computes $F_s(x)$ for values $x \in [0,k]$. If $ F_s(k)< p$, we have to recompute the problem with larger $k$, otherwise an algorithm can be used to find the point. However, whilst typically $\lim_{r \rightarrow \infty}F_s(r) = 1$, if $Pr(s \vDash \Diamond B) < 1$, then this no longer holds. See Section <ref>.

However, in the setting where this is true, i.e. when $Pr(s \vDash \Diamond B) = 1$, a guide to arrive at a sufficiently large interval $[0,k]$ at which $F_s(k) = 1$ is to use the expected value and variances of the first passage-reward density, for example
$$k = ExpRew(s \vDash \Diamond B) + \gamma VarRew(s \vDash \Diamond B)$$
where $\gamma$ is a parameter used to increase the range involved. The computation for $VarRew(s \vDash \Diamond B)$ can be found in <cit.>.

§.§ $\mrmpdf$ properties

System of convolution equations We present a computation of $\mrmpdf$ for each state of a finite sMRM with a common graph-based technique. Let us first state the solution, and leave the derivation till later.

Firstly for any states in $B$, then their cumulated rewards are assigned to zero. We also make them entirely self-absorbent, and we fix the reward transition distributions (for these self-loops) as the zero distribution, expressed by the Dirac delta (or Kronecker in the discrete case) $\delta_{x,0}(x)$. Consequently then, for any path absorbed by $B$, further reward is not cumulated in the sMRM process after entering $B$.

We can shorten $\mrmpdf$ to $f_{s,B}(r)$. However, if we write $f_s(r)$ it is assumed that the set of states to reach is $B$. Then  define $f_{rew(s,t)}(r)$ as the reward pdf for the transition $s \rightarrow t$. Let $Pre^*(B)$ denotes the set of states that reach $B$ with non-zero probability. Then define $S_? = Pre^*(B)\char`\\B$, to be the set of states that reach $B$ but are not in it. Firstly, since for each state $u \in B$ its reward  is zero, we have

$$f_u(r) = \delta_{x,0}(r)$$

which is the zero distribution (returns zero with probability one). Then for the remaining states, we have

f_s(r)    = ∑_t ∈S(P(s,t)f_rew(s,t) ∗f_t)(r)

= ∑_t ∈S_?(P(s,t)f_rew(s,t) ∗f_t)(r) + ∑_u ∈BP(s,u)f_rew(s,u)(r) ∗f_u(r)

= ∑_t ∈S_?(P(s,t)f_rew(s,t) ∗f_t)(r) + ∑_u ∈BP(s,u)f_rew(s,u)(r) ∗δ_x,0(r)

= ∑_t ∈S_?(P(s,t)f_rew(s,t) ∗f_t)(r) + ∑_u ∈BP(s,u)f_rew(s,u)(r)

for all $s \in S_?$. This then yields what we call a system of convolution equations.

Then since $S_? \cup B = S$ (due to all states reaching $B$ with probability 1), we have a solution for all states. Intuitively, the expression means that the reward accumulated from a state $s$ before arriving at $B$, is the convex combination of the rewards accumulated by all states that $s$ can immediately transition into (the combination weighted by the probability of entering these states).

Instead of solving a system of convolution equations, we can convert it to a set of systems of linear equations via the Fourier transform. Let the characteristic function transform operator (or Fourier transform) be represented as $\mathcal{F}$, and the transform applied to a function $f(x)$ be denoted as $\mathcal{F}\{f(x)\}$. Then, the Fourier transform of $\mrmpdf$ (which is $f_s$ from (<ref>)) is

ℱ{} = ℱ{}
= ∑_π̂ ∈ΠPr(π̂)ℱ{(r | π̂)}

which always exists, since Fourier transforms of pdfs and pmfs always exist.

From now on, we will denote $\mathcal{F}\{\mrmpdf\}$ as $\phi_{s}(\tau)$. And $\mathcal{F}\{\mrmpdfsub(r | \hat{\pi})\}$ as $\phi_{Rew(\hat{\pi})}(\tau)$. We also assume that the temporal logic property is always $\Diamond B$. Now we rewrite the equation (<ref>) above as

$$ \mrmfouriers =  \sum_{\hat{\pi} \in \Pi} Pr(\hat{\pi})\mrmfourierrews  $$

The results above can be derived knowing the properties for Fourier transforms over functions; its linearity with respect to constants and addition. Alternatively, using the law of total expectation: using traditional probability notation, then for two dependent random variables $X,Y$ we have $\mathcal{F}\{X\} = E[e^{itX}] = \mathbb{E}[\mathbb{E}[e^{itX}|Y] = \sum_{\hat{y} \in Y}P(\hat{y})*\mathbb{E}[e^{itX}| \hat{y}] = \sum_{\hat{y} \in Y}P(\hat{y})*\mathcal{F}\{X|y\}$. If this is not clear, let $X \sim Pr(Rew=r)$ with the domain being the real line, and $Y \sim Pr(\hat{\pi})$ with the domain being $\Pi$, then the above holds.

If using characteristic functions instead, i.e. using $\mrmfouriers =     \mathcal{F}\{\mrmpdf\}$, we have

ϕ_s(τ) = ∑_t ∈S_?P(s,t)ϕ_rew(s,t)(τ)ϕ_t(τ) + ∑_u ∈BP(s,u)ϕ_rew(s,u)(τ)

providing us with a system of linear equations, for each $\tau$.  And for each $u \in B$,

$$\phi_u(r) = \mathcal{F}\{\delta_{r,0}(r)\}=1$$

There is a one-to-one correspondence between the two formulas (<ref>) and (<ref>). This is due to the Fourier transform being bijective.

We now solve the problem presented in Sec. <ref>. The characteristic functions will be used instead of pdfs.

For each transition $s\rightarrow t$  of the system we assign the symbolic reward $rew(s,t)(r)$, of which the Fourier transform is $\mathcal{F}\{rew(s,t)(r)\} = \phi_{rew(s,t)}(\tau)$. The goal is to compute $\mrmfouriers$, for the temporal property $\Diamond B$, where $B =\{l_2, l_4 \}$. In the following, we may drop $(\tau)$ from our notation for a simpler representation.

Firstly, since $B = \{l_2,l_4\}$, then we make them self-looping and assign the loop reward transition the characteristic function of the Dirac delta. Then, using equation (<ref>), we arrive at the following linear system of equations for this sMRM:

[shorten >=2pt,node distance=3cm,auto]
[state,initial] (l0) at (0,0) $l_0$;
[state, right of=l0] (l1)  $l_1$;
[state,accepting, above of=l1]  (l2) $l_2$;
[state, right of=l1] (l3)  $l_3$;
[state,accepting, right of=l3] (l4)  $l_4$;

[->] (l0) edge node $0.2\phi_{rew(0,1)}$ (l1);
[->] (l0) edge[bend left]node $0.8\phi_{rew(0,2)}$ (l2);
[->] (l1) edge node $0.9\phi_{rew(1,3)}$ (l3);
[->] (l1) edge [bend right] node $0.1\phi_{rew(1,2)}$ (l2);
[->] (l2) edge [loop above] node $\phi_{\delta}$ (l2);
[->] (l3) edge[bend left] node $0.4\phi_{rew(3,1)}$ (l1);
[->] (l3) edge node $0.6\phi_{rew(3,4)}$ (l4);
[->] (l4) edge [loop above] node $\phi_{\delta}$ (l4);

A stochastic Markov reward model (sMRM) with two absorbing states $l_2,l_4$. Each transition has been annotated with the probabilities multiplied by the Fourier transform of the reward/cost distribution (which we have left as variables rather than define them explicitly).

ϕ_0    = 0.2ϕ_rew(0,1)ϕ_1 + 0.8ϕ_rew(0,2)ϕ_2

ϕ_1    = 0.1ϕ_rew(1,2)ϕ_2 + 0.9ϕ_rew(1,3)ϕ_3

ϕ_2    = ϕ_δ

ϕ_3    = 0.4ϕ_rew(3,1)ϕ_1 + 0.6ϕ_rew(3,4)ϕ_4

ϕ_4    = ϕ_δ

One can then derive manually that the solution for $\phi_0,\phi_1,\phi_3$ (i.e. characteristic functions for state $s_0,s_1,s_3$) is

ϕ_0    = 0.2ϕ_rew(0,1)ϕ_1 +  0.8ϕ_rew(0,2)ϕ_δ

ϕ_3    = 0.4ϕ_rew(3,1)ϕ_1 + 0.6ϕ_rew(3,4)ϕ_δ

ϕ_1    = 0.1ϕ_rew(1,2) + 0.9ϕ_rew(1,3)0.6ϕ_rew(3,4)/1 - (0.9ϕ_rew(1,3)0.4ϕ_rew(3,1))

The manual solution here can be automated via a symbolic solver such as Sympy, but as we shall see later, we will need to move away from solving the solution symbolically due to symbolic solvers not being very scalable. The reason is simply the intractability that comes with symbolic solvers.

We proceed to detail how equation (<ref>) was derived. It derives indirectly the equation (<ref>).

We are given an sMRM  with state space $S$, with a set of goal states $B \subset S$. Let it be the case that every state in $S$ can reach $B$. Then we assign to each state $u \in B$ zero rewards, such that the first-passage reward density is the (Dirac or Kronecker) delta $\delta(r)$. Let $S_? = S\char`\\B$ be the set of states not in $B$ (but can reach $B$).  For each $s \in S$, let $\mrmfouriers$ be the Fourier transform of $\mrmpdf$, the first-passage reward density, i.e. $\mrmfouriers = \mathcal{F}\{\mrmpdf\}$. Thus for each state $u \in B$, we immediately have $\phi_u{(\tau)} =\mathcal{F}\{\delta\}(\tau) = \phi_{\delta}{(\tau)} = 1 $.

Then, for each state $s \in S_?$, we have the equivalence

$$\phi_s(\tau) = \sum_{t \in S_?}\textbf{P}(s,t)\phi_{rew(s,t)}(\tau)\phi_t(\tau) + \sum_{u \in B}\textbf{P}(s,u)\phi_{rew(s,u)}(\tau)$$

Proof: Our derivation methodology is analogous to that found in <cit.>, except there the author derived solutions for the expected reachability reward (see Question <ref>), and the variance of the reachability reward. This type of derivation allows us to obtain the set equations via paths.

Let us define $\Pi.s$ to be the set of (finite) paths starting in $s$ ending at a state in $B$. This is a shortened variant of the earlier notation which would require us to write $\Pi_{\Diamond B}.s$ instead.

Then for each $s \in S_?$:

= {  definition of $\phi_s(\tau)$ }

∑_π̂ ∈Π.s Pr(π̂) ϕ_Rew(π̂)(τ)

= { letting $\hat{\pi} = s.\hat{\psi}$. Then since $s \notin B$, $\hat{\psi} \neq \emptyset$.  }

∑_t ∈S  ∑_ψ̂ ∈Π.t Pr(s.ψ̂) ϕ_rew(s.ψ̂)(τ)

= { definition of path probability: $Pr(s.\hat{\psi}) = \textbf{P}(s,t)Pr(\hat{\psi})$  }

∑_t ∈S  ∑_ψ̂ ∈Π.t P(s,t)Pr(ψ̂) ϕ_rew(s.ψ̂)(τ)

= { sum of 2 rvs using their cfs:  ϕ_rew(s. ψ̂)(τ) =ϕ_rew(s,t)(τ)ϕ_rew(ψ̂)(τ) }

∑_t ∈S  ∑_ψ̂ ∈Π.t P(s,t)Pr(ψ̂)ϕ_rew(s,t)(τ)ϕ_rew(ψ̂)(τ)

= { distribute P(s,t) and $\phi_{rew(s,t)}(\tau)$ outside since they are independent of $\hat{\psi}$  }

∑_t ∈SP(s,t)ϕ_rew(s,t)(τ)∑_ψ̂ ∈Π.t Pr(ψ̂)ϕ_rew(ψ̂)(τ)

= { simplify since by definition$\sum_{\hat{\psi} \in \Pi.t} Pr(\hat{\psi})\phi_{rew(\hat{\psi})}(\tau) = \phi_t(\tau)$ }

∑_t ∈SP(s,t)ϕ_rew(s,t)(τ)ϕ_t(τ)

However, since we have that $\phi_u(\tau) \triangleq \phi_\delta(\tau)$ for all states $u \in B$, we  now have that

ϕ_s(τ)    = ∑_t ∈S_?P(s,t)ϕ_rew(s,t)(τ)ϕ_t(τ) +
∑_u ∈BP(s,u)ϕ_rew(s,u)(τ)ϕ_u(τ)

= ∑_t ∈S_?P(s,t)ϕ_rew(s,t)(τ)ϕ_t(τ) +
∑_u ∈BP(s,u)ϕ_rew(s,u)(τ)

which completes the proof.

Remark:  If we apply the inverse Fourier transform to the formula above, we arrive at
$$f_s(r) = \sum_{t \in S_?}((\textbf{P}(s,t)f_{rew(s,t)}) \ast f_{t})(r) + \sum_{u \in B}\textbf{P}(s,u)f_{rew(s,u)}(r)$$

which provides us with a system of convolution equations.

§.§.§ Matrix notation

We now present the solution to $\mrmpdf$ in matrix form. Define

\begin{align*}
\textbf{f} &\triangleq (f_s)_{s \in S_?} &\\
A &\triangleq (\textbf{P}(s,t))_{s,t \in S_?}&\\
\textbf{G} &\triangleq (G_{s,t})_{s,t \in S_?} \triangleq (f_{(rew(s,t))})_{s,t \in S_?}&\\
\textbf{h} &\triangleq  (h_s)_{s \in S_?} \triangleq (\sum_{u \in B}\textbf{P}(s,u)f_{rew(s,u)})_{s \in S_?}
\end{align*}

where each pdf of $f_s$, $f_{rew(s,i)}$ with $s \in S_?,\ i \in S$ can be represented as an analytical function, denoted symbolically (or algebraically). Alternatively, we can use vectors to represent them, such that

\begin{align*}
f_s &= (f_s(r))_{r \in \mathbb{R}} &\\
f_{rew(s,i)} &= (f_{rew(s,i)}(r))_{r \in \mathbb{R}}
\end{align*}

This implies that $\textbf{G}$ is a `three-dimensional matrix', i.e. a hypermatrix, and both f,h are two-dimensional vectors, or hypervectors. We will use the latter representation for the majority of our work however we do avoid it in the examples of this chapter. Using this representation, in the continuous case, G has dimensions $(\mathbb{R} \times S_? \times S_?)$, whilst f,h have dimensions $(\mathbb{R} \times S_?)$.

Notation From now on, we will write regular vectors and matrices in the ordinary manner e.g. $f,g$ for vectors and $A,B$ as matrices. Hypervectors and hypermatrices are similar except bold instead, e.g. h,z for hypervectors and C,D for hypermatrices. This also applies to vectors and matrices of functions. There will be some exceptions to this, for example P is the probability matrix, or map over Markov chain transitions. Additionally, when we index hypermatrices or hypervectors, we generally make their letters light again, e.g. $\textbf{G}$, but $G_{s,t}$, or $G_{s,t}(r)$. We do not change cases for subscripted terms, e.g. not $a_{i,j}$, but $A_{i,j}$ for a matrix $A$.

Now we can rewrite (<ref>) in matrix form as

f = (A∘G)$\ast$ f + h

where $\ast$ denotes hypermatrix convolution, analogous to matrix multiplication, and $A\circ\textbf{G}$ is the element-wise hypermatrix-matrix multiplication between the matrix $A$ and hypermatrix $\textbf{G}$. Let us now explicate the definition of the operators on hypermatrices.

Let  $\textbf{A},\textbf{B}, \textbf{C}$ be three hypermatrices, with dimensions $(\mathbb{R} \times N \times M), (\mathbb{R} \times M \times L)$ and $(\mathbb{R} \times M \times L)$ respectively. Then,

* The hypermatrix convolution product $\textbf{D} = \textbf{A} \textcircled{$\ast$} \textbf{B}$, is defined with dimensions $(\mathbb{R} \times N \times L)$, of which values $D_{i,j}(r)$, for all $i \in N, j \in L, r \in \mathbb{R}$, are determined by

$$D_{i,j}(r) =  \sum_{q=0}^{N-1}(A_{i,q} \ast B_{q,j})(r) $$

* The hypermatrix multiplication product $\textbf{U} = \textbf{AB}$, is defined with dimensions $(\mathbb{R} \times N \times L)$ of which values $U_{i,j}(r)$, for all $i \in N, j \in L, r \in \mathbb{R}$, are determined by

$$U_{i,j}(r) = \sum_{q=0}^{N-1}A_{i,q}(r)B_{q,j}(r)$$

or $U(r) = A(r)B(r)$ where the right-hand side is a matrix multiplication.

* The resulting sum/subtraction $\textbf{S} = \textbf{B} \pm \textbf{C}$, is defined with dimensions equal to $\textbf{B}$ of which values $S_{i,j}(r)$, for all $i \in M, j \in L, r \in \mathbb{R}$, are determined by

$$S_{i,j}(r) =  A_{i,j}(r) \pm B_{i,j}(r) $$

or $S(r) = A(r) \pm B(r)$.

We now define operators between hypermatrices and matrices.

Let $D$ be a matrix of dimensions $(N \times M)$ and $\textbf{A}$ the hypermatrix with dimensions $(\mathbb{R} \times N \times M)$. Then,

* The hypermatrix-matrix Hadamard product $\textbf{H} = D\circ \textbf{A} = \textbf{A} \circ D$ is a hypermatrix with dimensions equal to A, with values $H_{i,j}(r)$ determined by

$$ H_{i,j}(r) =  D_{i,j}A_{i,j}(r) = A_{i,j}(r)D_{i,j}$$

for $i,j \in (N \times M)$ and $r \in \mathbb{R}$, or alternatively  $ H(r) = D \circ A(r) =  A(r) \circ D $
where  $ D \circ A(r)$ is the Hadamard product (or element-wise product) between the two matrices.

* Let $E$ be a matrix with dimensions $(O \times N)$. The hypermatrix-matrix multiplication product $\textbf{Z} = E\textbf{A}$ is a hypermatrix with dimensions equal to ($\mathbb{R} \times O \times M$), with values $Z_{i,j}(r)$ determined by

$$Z_{i,j}(r) = \sum_{q=0}^{N-1}E_{i,q}A_{q,j}(r)$$

or $Z(r) = EA(r)$.

* The hypermatrix-matrix sum/subtraction $\textbf{F} = \textbf{A} \pm D = \mp D + \textbf{A}$, is defined with size equal to $\textbf{A}$, with values $F_{i,j}(r)$ determined by

$$F_{i,j}(r) =  A_{i,j}(r) \pm  D_{i,j} =  \mp D_{i,j} + A_{i,j}(r)$$

or $F(r) =  A(r) \pm  D =  \mp D + A(r)$.

Moving forward, if we use characteristic functions to represent our problem, then firstly define

\begin{align}
\label{eqn:definitions_terms_DFT}
\textbf{x} &\triangleq (x_s)_{s \in S_?} \triangleq (\phi_s)_{s \in S_?} &\\
\textbf{C} &\triangleq  (C_{s,t})_{s,t \in S_?} \triangleq (\phi_{(rew(s,t))})_{s,t \in S_?} \notag &\\
\textbf{d} &\triangleq  (d_s)_{s \in S_?} \triangleq (\sum_{u \in B}\textbf{P}(s,u)\phi_{rew(s,u)})_{s \in S_?} \notag
\end{align}

Then the set of equations (<ref>) has the hypermatrix form,

x(τ) =(A∘C(τ))x(τ) + d(τ)

The solution of which is just

x = (I - (A∘C))^-1d

with $I$ being the identity matrix of size $|S_?| \times |S_?|$. The inverse of a square hypermatrix $\textbf{B}^{-1}$ (having an equal number of rows and columns) is defined such that $\textbf{B}^{-1}\textbf{B} = I\circ \textbf{1}$ with 1 having size equal to B. If we use the analytical (or symbolic) form of these characteristic functions, the general solution to this system will give us rational characteristic functions as we have seen in an earlier example and was shown earlier by for example <cit.>.

The system in (<ref>) can be rewritten in the form

(I - (A∘C))x = d

Then, the matrix $Z(\tau) = (I - A\circ \textbf{C}(\tau))$ is weak and irreducibly diagonally dominant, i.e. there exists a set of states $q \in Q \subseteq S_?$, such that ${Z(\tau)}_{q,q}(\tau)$ exhibits strict diagonal dominance, i.e.

∑_t ∈S_?∖q|Z(τ)_q,t(τ)| < |Z(τ)_q,q(τ)|

and for any remaining states $s \in S_?\setminus Q$, we have weak diagonal dominance, i.e.

∑_t ∈S_?s|Z(τ)_s,t(τ)| ≤|Z(τ)_s,s(τ)|

Proof:  Making use of the fact that for all $\tau$, $|C_{s,t}(\tau)| \leq 1$ (see Definition <ref>), and that $A_{s,t} = \textbf{P}(s,t) \leq 1$, for all $s,t \in S_?^2$, we can easily deduce that $Z(\tau)$ is weakly diagonally dominant. To show that it is also irreducibly diagonally dominant, then for some row $i$ of $Z$, we have that $\sum_{t \in S_?} A_{i,t} < 1$. This is since there must be some state that reaches $B$ with non-zero probability. Knowing this, it is trivial to show that for any row such as $i$, $Z$ is strictly diagonally dominant.

The system of equations (<ref>) has a unique solution, and therefore (<ref>) is resolvable.

Proof: From Proposition <ref>, we know that $(I - A\circ \textbf{C}(\tau))$ is weak and irreducibly diagonally dominant. Then, this is a sufficient condition to prove that a unique solution exists, a proof of which can be found in <cit.>. Hence, we have a unique solution for all $\tau$.

Corollary, this theorem proves that reward random variables positive and negative can be used together within a single sMRM and a unique solution will still exist. As another corollary, this implies that a unique solution exists in the time domain (i.e. without having transformed the equations), i.e. since the Fourier transform is bijective and always exists for random variables. [Matrix (symbolic) approach]

Continuing with the study problem in Example <ref>, instead of solving the system manually, we can now solve it automatically via symbolic matrix solvers. One such solver is Sympy <cit.>.  Note that in this example, we use characteristic functions once more and also forgo the vector representation and use symbolic terms, as previously done in Example  <ref>.

Firstly, let us denote terms in the form $p_{i,j}\phi_{rew(i,j)}$ as $\sigma_{i,j}$ for each pair $i,j \in S_?^2$. Then we rewrite the set of linear equations in matrix form  (<ref>) for states in $S_? = {l_0,l_1,l_3}$  giving:

\phi_0 \\
\phi_1 \\
\phi_3
\end{pmatrix} =  \begin{pmatrix}
0 & \sigma_{0,1} & 0 \\
0 & 0 & \sigma_{1,3}\\
0 & \sigma_{3,1} & 0
\end{pmatrix}
\begin{pmatrix}
\phi_0 \\
\phi_1 \\
\phi_3
\end{pmatrix}    +
\begin{pmatrix}
\sigma_{0,2} \\
\sigma_{1,2} \\
\sigma_{3,4}
\end{pmatrix}

Then by (<ref>),

$$    \begin{pmatrix}
\phi_0 \\
\phi_1 \\
\phi_3
\end{pmatrix} =
\begin{pmatrix}
1 & -\sigma_{0,1} & 0 \\
0 & 1 & -\sigma_{1,3}\\
0 & -\sigma_{3,1} & 1
\end{pmatrix}^{-1}
\begin{pmatrix}
\sigma_{0,2} \\
\sigma_{1,2} \\
\sigma_{3,4}
\end{pmatrix}

Solving this system (symbolically) via     (<ref>), gives:

\phi_0 \\
\phi_1 \\
\phi_3
\end{pmatrix}
\alpha
\begin{pmatrix}
\sigma_{0,1}(- \sigma_{3,4}\sigma_{1,3} - \sigma_{1,2}) +
\sigma_{0,2}(\sigma_{3,1}\sigma_{1,3} - 1) \\
\sigma_{3,4}\sigma_{1,3} + \sigma_{1,2} \\
\sigma_{3,1}\sigma_{1,2} +         \sigma_{3,4} \\
\end{pmatrix}

where $\alpha \triangleq \frac{1}{\sigma_{3,1}\sigma_{1,3} - 1}$.

These equations are equal to the equations found manually earlier. This was proven using Sympy <cit.> to show that these equations have the same solution set with the corresponding equations manually derived earlier.

§.§.§ Partial passage-reward densities

When presenting the solution to find the first-passage reward density  $f_s(r) = Pr(s \vDash \Diamond_{=r} B)$ (Question <ref>),
we enforced earlier that all states $s \in S$ could reach $B$. And thus for the set of paths $\Pi$ that reach $B$ beginning from $s_0$, we had
$$\sum_{\hat{\pi} \in \Pi.s_0}Pr(\hat{\pi}) = 1$$

Let us relax this assumption in two ways. Firstly, let there always be absorbing states within our sMRM which every state will eventually reach. Then we define this set of states as $Abs$, a subset of $S$, and let $\Pi_{\alpha}.s$ be the set of paths that end in a state $\alpha \in Abs$, beginning in $s$. We now get

$$\sum_{\alpha \in Abs}\sum_{\hat{\pi} \in \Pi_\alpha.s}Pr(\hat{\pi}) = 1 = \sum_{\alpha \in Abs} Pr(s \vDash \Diamond \alpha)$$

If we set $B = Abs$, we have $Pr(s \vDash \Diamond B) = \sum_{\alpha \in Abs} Pr(s \vDash \Diamond \alpha) = 1$ and if $B \subset Abs$, then $Pr(s \vDash \Diamond B) \leq 1$. Therefore with respect to accumulated rewards, we have in the case where $B= Abs$,

$$\int dr \cdot  f_s(r) = \int dr \cdot \mrmpdf   = 1$$

And if $ B \subset Abs$,

∫_-∞^∞dr ·f_s(r)    = ∫_-∞^∞dr ·

= ∫_-∞^∞dr ·∑_α∈B ∑_π̂ ∈Π_α.sPr(π̂)Pr(r | π̂)

= ∑_s ∈B ∑_π̂_̂ŝ ∈Π_α.sPr(π̂) ∫_-∞^∞dr ·Pr(r | π̂)

=  ∑_s ∈B ∑_π̂_̂ŝ ∈Π_α.sPr(π̂) ∫_-∞^∞dr ·f_rew(π̂)(r)

= ∑_s ∈B ∑_π̂_̂ŝ ∈Π_α.sPr(π̂) ·1

= Pr(s B)

of which implies that $f_s(r)$ generally represents a partial probability density function (seeing that it no longer always integrates to 1).

To obtain the first-passage reward density for such a system above, it is sufficient to remove all states that do not reach $B$ first and then solve it (e.g. via matrix inversion in Example <ref>). More generally if $B \subset Abs$, which includes the case where $B \cap Abs  =\emptyset$, then to determine the first-passage reward density  $Pr(s \vDash \Diamond_{=r} B)$, we first modify the sMRM such that states in $B$ are absorbing artificially. After this we remove any states that do not reach $B$ from the system and then solve.

§ EXPERIMENTING WITH SYMBOLIC INVERSION

We have found that symbolic matrix inversion (using the Sympy package for python <cit.>), is generally intractable for complete (square) matrices with more than six dimensions (or states in our case).

For the problem
$$(I - (A\circ \textbf{C}))\textbf{x} = \textbf{d}$$

where $\textbf{C},\textbf{x},\textbf{d}$ are matrices (or vectors) of functions. Let  $(I - (A\circ \textbf{C}))$ be completely symbolic (hence we are not using the vector representation of characteristic functions). For example for a $2 \times 2$ matrix, we have
a & b \\
c & d \\
\end{pmatrix}$$

where $a,b,c,d$ are strictly variables. Then our goal is to find $(I - (A\circ \textbf{C}))^{-1}$ symbolically.  From some basic experiments, we found the following results:

xlabel=Number of states,
ylabel=Time taken for matrix inversion (secs.),
xmin=0, xmax=8,
ymin=0, ymax=80,
grid style=dashed,
(1,    0.0004210472107)(2,    0.1848955154)(3,0.2642114162)(4,1.207051277)(5,51.1838026)

The time taken to invert the (dense) matrix $\textbf{I} - (\textbf{A} \circ \textbf{C})$ blows up drastically.

where we have given up trying to solve problems larger than six states. Whilst the above shows a bottleneck, it is not a fair analysis alone as usually there is only one variable $x$ of a univariate-reward sMRM.

x   2+x   3x +5   2x + 10    11x-9    8x-20

2+x   x   2+x   x    27x - 10    15x - 20

x   2+x   5x- 3   30x + 1    x+50 -10    11x + 12x^2

x+5 - 10x   20x+ x^2 +x^3   x+8 -10   2+x    42    36x

x+5 - 10x + 10x^3    6x   x+8 -10 - 25*3 + x^2   10+x    -15x    16x-4

x + x^2 + 3   6*x + 12   x- 10   2x   x^2   x-15

When using a single variable $x$, we find that solving times does get easier. For example, inverting the $6\times 6$ matrix above (<ref>) took around 8.5 seconds.  However, the resulting matrix from the inversion was unwieldy, being over hundreds of thousands of characters long. A simplification procedure exists in Sympy which took roughly 215 seconds to complete. The result however was still quite large being around 11000 characters long. One might argue that the matrix above is not a true representation of an sMRM problem as the functions above are arbitrary and not characteristic functions. Whilst true, certain characteristic functions have more complicated forms, such as those of the Weibull distribution, Gaussian or log-normal distribution. Thus with such functions, it may be better to use a separate variable to represent them. This however leads to the blow-up previously shown.

Therefore our current bottleneck is in the inversion process (sympy uses Gaussian elimination as the default algorithm <cit.> which is known to have a complexity of $O(n^3)$ for a matrix of dimensions $n\times n$)). Whilst better symbolic algorithms could be researched, to mitigate this problem, we proceed to replacing the symbolic step with a numerical one; the solutions presented in the next two chapters for discrete-reward sMRMs, and one subsequent chapter for continuous-reward sMRMs.

§ SUMMARY
In this chapter, we presented the theoretical foundations for sMRM model checking. We explained that many problems for sMRMs reduce to solving for the first-passage reward density $Pr(Rew = r\ \cap\ s \vDash \Diamond B)$. Expected value problems were found to reduce to algorithms similar to that for MRMs.

In the subsequent chapters, we move forward to deriving numericals solutions for the first-passage reward density. We first present a direct approach for solving the system of convolution equations, and in the following chapter we present iterative algorithms.

Note that whilst we have presented techniques in the following chapters for solving the system of convolution equations previously derived, there are other ways to derive the first-passage reward densities, for example see <cit.>.

CHAPTER: DIRECT METHODS FOR SMRMS

§ INTRODUCTION
In this chapter, direct algorithms for solving $\mrmpdf$ numerically are presented. A direct algorithm is one that terminates under a finite number of steps, and provides a full solution only upon termination. The first algorithm we present is the naive Gauss elimination algorithm, but adapted for systems of convolution equations. The second is an approximation algorithm. To summarize this chapter, our major finding is that Gaussian elimination can be adapted for the system of convolution equations by carefully replacing multiplication with convolution and division by deconvolution. Addition and subtraction amounts to piecewise addition or subtraction of vectors.

We shall restrict ourselves to using probability mass functions, i.e. discrete reward random variables (defined only over $\mathbb{N}$). Continuous rewards are treated in Chapter 7. Secondly, we present a computation of $\mrmpdf$ only for $r \in \{0,1,2,\cdots,k\}$, for some $k \in \mathbb{N}$. For any sMRM, we assume that all states $s \in S$ can reach $B$ with probability one. We also continue to focus on univariate reward random variable. We first present the necessary definitions for this chapter (of which will be used in subsequent chapters). We then define the algorithms aforementioned.

§ MOTIVATION
As previously discussed, the reason for finding new algorithms is due to the symbolic approach presented in the previous chapter being intractable. The papers <cit.> provide for us the basis of the following work. They suggest replacing characteristic functions of probability density functions with the discrete Fourier transform (DFT) instead via discretization of the characteristic function. This is useful as there exists a fast numerical implementation of the DFT known as the fast Fourier transform (FFT) with which we can perform linear convolutions. However in this case we restrict ourselves to probability mass functions (e.g. discrete R.V.s) by choice, and therefore any distribution to be used that is continuous should be discretized first, or a discrete analogue found for it. For example see <cit.> for a survey of analogues.

Before we present the Gaussian elimination algorithm, we will first define the convolution and deconvolution operations that will be used.

§ DISCRETE CONVOLUTION AND DECONVOLUTION

Given a non-negative discrete lattice random variable $X$ defined over $\mathbb{N}$, let $p_X[x]$ be its probability mass function (pmf). From now on, we represent pmfs as vectors (arrays) (with indices in $\mathbb{N}$), and are not expressed analytically (algebraically).  Additionally, for the computation of $\mrmpdf$, we will restrict ourselves to $r \in \{0,1,2,\cdots,k\} \in \mathbb{N}$ only. Therefore for example, $f_s[r]$ will be computed only for such a range.

The discrete Fourier transform operator $\mathcal{D}$ applied to the pmf of a random variable $X$ is written as $\mathcal{D}\{p_X[x]\}$, with its result being denoted as $\varphi_X[\tau]$. Then when applying $\mathcal{D}$ to $\mrmpdf$, we write this transform as

$$\mathcal{D}\{\mrmpdf\} = \varphi_{s \vDash \Diamond B}[\tau]$$

or $\varphi_{s}[\tau]$ for short.

§.§.§ Sums of random variables & linear convolution via DFTs
Given two independent random variables $X,Y$ defined over $\mathbb{N}$ (for simplicity), with pmfs $p_X, p_Y$ respectively, then the pmf of the operation $X+Y$ is determined via a discrete linear convolution, i.e.

p_Z(k) = (X+Y)(k) = ∑_x=0^kp_X(x)p_Y(k-x)

If we intend to compute the first $k$ values of $p_Z(x)$, i.e. for all $x \in \{0,1,\cdots,k-1\}$, then using DFTs to perform the convolution, the pmf of $Z$ can be computed as

p_Z(k)= 𝒟^-1{𝒟{p_X[x]}𝒟{p_Y[x]}}(k)

where $\mathcal{D}^{-1}$ is the inverse DFT operator, and $\widetilde{p_i}$ is defined as

p_i[x] ≜
p_i[x]     0 ≤x ≤k

0      k < x ≤2k-1

for each $i \in \{X,Y\}$. Note that $\widetilde{p_i}$ is just $p_i$ truncated and zero padded (with $k-1$ zeroes). The zero-padding is necessary with DFTs to prevent a known problem called time-aliasing. Ultimately, the above states that the convolution requires two vectors of length $2k-1$, to determine the first $k$ values of $p_Z$ exactly.

Let us introduce the function
$$\mathtt{conv}_k(v_1,v_2) \triangleq \mathcal{D}^{-1}\{\mathcal{D}\{\widetilde{v}_1[x]\}\mathcal{D}\{\widetilde{v}_2[x]\}\}$$
as notation to specify the convolution of vector $v_1$ with $v_2$ via FFTs for up to $k$ points only. If $v_1,v_2$ are derived from analytical pmfs, only the first $k$ values are sampled, and the resulting vector padded with $k-1$ zeroes. The length of both $\widetilde{v_1},\widetilde{v_2}$ is $2k-1$.

Efficient algorithms exist to compute DFTs (and their inverses) via the fast Fourier transform (FFT) with complexity $O(nlog_2n)$ where $n$ is the size of the vector being transformed. For example see <cit.>. Thus the computation of (<ref>) has complexity $3\cdot O(nlog_2n) + O(n)$ where $n=2k-1$. This is since $\mathcal{D}$ is used twice, $\mathcal{D}^{-1}$ once and element-wise vector multiplication is performed once. Asymptotically, this is equivalent to $O(klog_2k)$.

§.§.§ Linear deconvolution via polynomial division

Given three non-negative discrete lattice random variables $X,Y,Z$ defined over $\mathbb{N}$, where $Z = X+Y$, and $X$ is independent of $Y$. Then the deconvolution operation is sought when we know for example the pmfs of $Z$ and $Y$ only but we would like to find the pmf of $X$. Algebraically, then the pmf $p_X$ is computed as

\begin{align}
p_X[k] =\frac{1}{p_Y[0]}(p_Z[k] - \sum_{x=0}^{k-1}p_X[x]p_Y[k-x])
\end{align}

and for clarity, $p_X[0] = \frac{p_Z[0]}{p_Y[0]}$.

This can be validated by simply rearranging the definition of convolution. Or if starting from the result above, we obtain

\begin{align*}
p_X[k] &=\frac{1}{p_Y[0]}(p_Z[k] - \sum_{x=0}^{k-1}p_X[x]p_Y[k-x]) &\\
p_X[k]p_Y[0] &=p_Z[k] - \sum_{x=0}^{k-1}p_X[x]p_Y[k-x] &\\
p_Z[k] &= p_X[k]p_Y[0] + \sum_{x=0}^{k-1}p_X[x]p_Y[k-x]  &\\
p_Z[k] &= \sum_{x=0}^{k}p_X[x]p_Y[k-x]  &\\
p_Z[k] &= (X+Y)[k]
\end{align*}

which is the original definition of $p_Z[k]$.

Note that deconvolution is computed recursively as seen above (unlike convolution). As such, computing $p_X[k]$ requires an ordered computation, i.e. computing (in order) $p_X[0],p_X[1],\cdots,p_X[k-1]$ before $p_X[k]$. In another sense, deconvolution is a recursive filter, but convolution is non-recursive.

An important thing to note is that if $p_Y[0] = 0$, the computation above will lead to a division by zero, i.e. $1/p_Y[0]$. Generally, it is only possible to resolve this if $p_X,p_Y$ are not of fixed length, e.g. infinitely long. Starting from the third equation above, then

\begin{align*}
p_Z[k] &= p_X[k]p_Y[0] + \sum_{x=0}^{k-1}p_X[x]p_Y[k-x]  &\\
&= 0 + \sum_{x=0}^{k-1}p_X[x]p_Y[k-x]  &\\
&= p_X[k-1]p_Y[1] + \sum_{x=0}^{k-2}p_X[x]p_Y[k-x]
\end{align*}

Then, we can rearrange the result to obtain $p_X[r]$:

\begin{align*}
p_Z[k] &= p_X[k-1]p_Y[1] + \sum_{x=0}^{k-2}p_X[x]p_Y[k-x]  &\\
p_X[k-1] &= \frac{1}{p_Y[1]}(p_Z[k] - \sum_{x=0}^{k-2}p_X[x]p_Y[k-x]) &\\
p_X[k] &= \frac{1}{p_Y[1]}(p_Z[k+1] - \sum_{x=0}^{k-1}p_X[x]p_Y[k+1-x]) &\\
p_X[k] &= \frac{1}{p_Y[1]}(p_Z[k+1] - \sum_{x=0}^{k-1}p_X[x]p_Y[k+1-x])
\end{align*}

where $p_X[0] = \frac{p_Z[1]}{p_Y[1]}$. From a computational point of view, this result is equal to computing $p_X  =p_Z[1:] \deconv p_Y[1:]$, where for example $p_Z[1:]$ is the vector $p_Z$ with each element shifted left by one place. However this latter computation requires the values of $p_Z[1:k+1]$ to compute the first $k$ values of $p_X$. If we only have access to the first k values of $p_Z$ this computation will yield therefore only $k-1$ values which is problematic for our system of convolution equations since the vectors are fixed to some finite length $k$. In the context of finitely long vectors, let us define a full deconvolution $f \deconv g$ to be where $g[0] \neq 0$. We will call  $g$ a full deconvolutor.

Another area where convolution arises is with polynomial multiplication. This is since the coefficients of the polynomial resulting from multiplication are derived via convolution of the coefficients of the multiplicands. For example, we can encode a pmf $p_X[r]$ (with $r = 0,1,\cdots,k-1$) as a polynomial of degree $k-1$[Note that the degree of the polynomial is one less than the length of $p_X$.] by
$$P_{k-1}(x) = p_X[0] + p_X[1]x + p_X[2]x^2 + \cdots + p_X[k]x^{k-1}$$

and we will denote $P_{k-1}(x;p_X)$ to mean that the coefficients of the polynomial is encoded via the vector $p_X$. Then for example, we have the multiplication $P_{k-1}(x;p_X)\times P_{k-1}(x;p_Y) = P_{2k-2}(x;p_Z)$ where $p_Z = p_X * p_Y$.

Likewise, the inverse operation - polynomial division - can be used to reverse the process, e.g. $ P_{2k-2}(x;p_Z)/P_{k-1}(x;p_Y) = P_{k-1}(x;p_X) + R(x)$, where $R$ is the remainder polynomial and $P_{k-1}(x;p_X)$ is termed the quotient. Hence, we can use polynomial division for deconvolution since the quotient obtains $p_X = p_Z \deconv p_Y$. Note that the remainder $R$ is not required and can be discarded. The division above assumes $p_Z$ is known up to $2k-1$ values, and $p_Y$ up to $k$ values. However, if we are given only $k$ values of $p_Z$ and $p_Y$ we can still obtain $P_{k-1}(x;p_X)$ (and therefore the first $k$ values of $p_X$) by zero-padding $p_Z$ with k-1 zeroes. This is provided that we have a full deconvolutor.

Let $v_1,v_2$ be vectors of length $k$. Then define the polynomial

$$ P_{k-1}(x;v_3) + R = P_{2k-2}(x;\widetilde{v_1})/P_{k-1}(x;v_2)$$

where $\widetilde{v_1}$ is equal to $v_1$ zero padded with $k-1$ zeroes. Then, we introduce the function

$$\mathtt{deconv}_k(v_1,v_2) \triangleq \mathtt{coeffs}(P_{k-1}(x;v_3))$$

where $\mathtt{coeffs}(P_{k-1}(x;v_3)) = v_3$, i.e. $\mathtt{coeffs}$ returns the coefficients of the polynomial. The subscript $k$ in front of $\mathtt{deconv}$ denotes that the first $k$ points of deconvolution is obtained.

Let $p_X,p_Y$ be pmfs. Then to obtain the quotient $p_Z[x] = p_X \deconv p_Y$ for $x = 0,1,\cdots,k-1$, define $v_1[x] = p_X[x], v_2[x] = p_Y[x]$ for $x = 0,1,\cdots,k-1$. Then $p_Z[x] = \mathtt{deconv}_k(v_1,v_2)$.

There appears to be algorithms that compute polynomial divisions using FFTs, for example see the lecture notes <cit.>. They are stated to have the same worst-case complexity as the convolution case with the FFT - $O(nlog_2n)$ where $n=2k-1$ for the computation of $\mathtt{deconv}_k$. Note however that the deconvolution operation is generally numerically unstable.

§ NAIVE GAUSSIAN ELIMINATION

We can rearrange (<ref>) into the form

((I ∘Δ_x,0) -A∘G) $\ast$ f = h

where $\bm{\Delta}_{x,0}$ is a hypermatrix where each element is the Kronecker delta $\delta_{x,0}$ with mass 1 over zero (and thus a hypermatrix overall). Thus $(I \circ \bm{\Delta}_{x,0})$ returns a diagonal hypermatrix instead, which represents the equivalent of the identity matrix in traditional linear algebra. Using vectors for pmfs (of discrete (lattice) random variables), then each of $\textbf{f}$, $(A \circ \textbf{G})$, and $\textbf{h}$ is a three dimensional hypermatrix or hypervector. Again, we will restrict ourselves to $r \in \{0,1,2,\cdots,k\} \in \mathbb{N}$ only, e.g. $(G_{s,t}[r])_{s,t \in S_?^2}$ is defined for $r=0,1,2,\cdots$. Since, we are interested in computing $\mrmpdf$ only for $r = 0,1,2,\cdots,k-1$, i.e. just $k$ values, it suffices to fix the sizes of each term as

* $(1 \times {S_?| \times |S_?|})$ for $A$ and $I$.

* $({k \times |S_?| \times |S_?|})$ for $\textbf{G}$ and $\bm{\Delta}_{x,0}$.

* $({k \times |S_?| \times 1})$ for $\textbf{f}$ and $\textbf{h}$.

Let us denote $((I \circ \bm{\Delta}_{x,0}) -A\circ \textbf{G})$ as $\mathbf{\mathcal{A}}$, leaving us with
$$\mathbf{\mathcal{A}}\ \textcircled{$\ast$}\ \textbf{f} = \textbf{h}$$

We refer to the second and third dimension as the rows and columns of the hypermatrix respectively, i.e. $\mathcal{A}$ has size $({k\times \text{no. rows} \times \text{no. columns}})$. Hence $\mathcal{A}_{i,j}(r)$ refers to indexing the hypermatrix by the $i^{th}$ row, the $j^{th}$ column, and the $r^{th}$ value. The indices start from zero, i.e. $i=0,1,\cdots,|S_?|-1$. The $i^{th}$ equation  of the system is denoted as

$$(\mathbf{\mathcal{A}}\ \textcircled{$\ast$}\ \textbf{f})_i =  \sum_{j} \mathcal{A}_{i,j} \ast {f}_{j} =  h_i

We now present a Gaussian elimination (GE) type algorithm to solve for $\textbf{f}$ in the system above. It is helpful to realise at this stage that convolution and deconvolution are used for the same purposes as multiplication and division in the traditional algorithm (for linear equations), whilst addition and subtraction are just element-wise vector addition and subtraction.

Description of GE:
Gaussian elimination solves the system in two stages: The first stage involves eliminating from each equation $i$ where $i \geq 1$, all terms $f_j$ where $j \leq i$. The procedure is done algorithmically: Starting from equation $0$, this equation is used to eliminate from equation $i$ for all $i\geq 1$, the term $f_0$. This can be done by subtracting from each equation $i$ a convolution (multiple) of equation 0. That is, we replace equation $i$ with

(∑_j 𝒜_i,j ∗f_j) - σ_i * (∑_j 𝒜_0,j ∗f_j) =  h_i - (σ_i * h_0)

where $\sigma_i \triangleq \mathcal{A}_{i,0} \deconv \mathcal{A}_{0,0}$ is the convolution (multiple) term. and it is assumed that $\mathcal{A}_{0,0}$ (called a pivot in the literature on Gaussian elimination) is a full deconvolutor (described in Section <ref>). This assumption is proved true in Theorem <ref>. The update for equation $i$ (<ref>) can also be written as

$$ \sum_{j} \mathcal{A}^{'}_{i,j} \ast {f}_{j} =  h^{'}_i$$

where $\mathcal{A}^{'}_{i,j}  = (\mathcal{A}_{i,j} - \sigma_{i} * \mathcal{A}_{0,j})$, and $h^{'}_{i} = h_i - (\sigma_i * h_0)$. Note then that

$$\mathcal{A}^{'}_{i,0} = \mathcal{A}_{i,0} - (\mathcal{A}_{i,0} \deconv \mathcal{A}_{0,0}) * \mathcal{A}_{0,0} = \mathcal{A}_{i,0} - \mathcal{A}_{i,0} = \vec{0}$$

Thus, $f_0$ is eliminated from equation $i$.

Given a system of convolution equations, substituting an equation by a sum of itself and a multiple of another does not change the solution set (see Proposition <ref>). Hence, the system after being repeatedly updated (via (<ref>)) will still have the same solution.

The procedure above ensures $f_0$ is eliminated from equations $i$, for all  $i\geq 1$. Then, re-assign $\mathcal{A} \coloneqq \mathcal{A}'$, and $\textbf{h} \coloneqq \textbf{h}'$, and the elimination procedure repeats as before, but now starting with equation 1, and eliminating $f_1$ from equations $i,$ for all $i \geq 2$. This repetition continues, e.g. starting with equation $j=2,3,\cdots,|S_?|-2$, and eliminating $f_j$, from all equations $i$ where $i>j$. The end product of all of this is the system:

\begin{align}
\label{eqn:row_echelon_form}
\begin{pmatrix}
\mathcal{A}'_{0,0} &  \mathcal{A}'_{0,1} & \mathcal{A}'_{0,2} & \cdots &\cdots & \mathcal{A}'_{0,|S_?|-1}\\
\vec{0} &   \mathcal{A}'_{1,1} & \mathcal{A}'_{1,2} & \cdots &\cdots& \mathcal{A}'_{1,|S_?|-1} \\
\vec{0} &  \vec{0} & \mathcal{A}'_{2,2} & \cdots &\cdots& \mathcal{A}'_{2,|S_?|-1} \\
\vdots & \vdots & \ddots & \ddots & \ddots & \vdots  \\
\vec{0} & \vec{0} & \vec{0} & \cdots & \mathcal{A}'_{|S_?|-2,|S_?|-2} & \mathcal{A}'_{|S_?|-1,|S_?|-2}\\
\vec{0} & \vec{0} & \vec{0} & \cdots & \vec{0} & \mathcal{A}'_{|S_?|-1,|S_?|-1}
\end{pmatrix}\textcircled{$*$}
\begin{pmatrix}
f_0 \\
f_1 \\
f_2 \\
\vdots \\
f_{|S_?|-2} \\
\end{pmatrix}
\begin{pmatrix}
h'_0 \\
h'_1 \\
h'_2 \\
\vdots \\
h'_{|S_?|-2} \\
\end{pmatrix}
\end{align}

or written in matrix form,

$$\mathcal{A}'\textcircled{$*$}\textbf{f} = \textbf{h}'$$

The resulting hypermatrix $\mathcal{A}'$ is upper-triangular, that is the entries $\mathcal{A}'_{i,j} = \vec{0}$ for all $i >j$.

Next, the second stage of Gauss elimination (known as back-substitution) begins. Note from (<ref>) that the solution hypervector f can be determined via the reverse order $f_{|S_?|-1}$, $f_{|S_?|-1}$, $\cdots$, $f_0$. Starting from the last equation of the system (<ref>) i.e. equation $|S_?|-1$,  we have

$$\mathcal{A}^{'}_{|S_?|-1,|S_?|-1}* f_{|S_?|-1} = h^{'}_{|S_?|-1}$$

and deconvolving both sides by $\mathcal{A}^{'}_{|S_?|-1,|S_?|-1}$  gives

$$f_{|S_?|-1} = h^{'}_i \deconv \mathcal{A}^{'}_{|S_?|-1,|S_?|-1}$$