Level-dependent QBD processes with finitely many levels: the rate matrices R^(n) by the backward matrix continued fraction, and the stationary distribution. More...

#include <cstddef>
#include <string>
#include <vector>
#include "line/api/mam/qbd_r.h"
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/linalg.h"
#include "line/util/lu.h"
#include "line/util/matrix.h"
#include "line/util/svd.h"

Include dependency graph for ldqbd.h:

Classes
struct	line::mam::LdqbdPi< T >
	Stationary distribution of a level-dependent QBD, per level and per phase. More...
struct	line::mam::LdqbdResult< T >
	R and the stationary distribution together (ldqbd.m). More...

Namespaces
namespace	line
namespace	line::mam

Functions
template<class T>
std::vector< Matrix< T > >	line::mam::ldqbd_R (const std::vector< Matrix< T > > &q0, const std::vector< Matrix< T > > &q1, const std::vector< Matrix< T > > &q2)
	Rate matrices R^(1), ..., R^(N) of a level-dependent QBD (ldqbd_R.m).
template<class T>
LdqbdPi< T >	line::mam::ldqbd_pi (const std::vector< Matrix< T > > &R, const std::vector< Matrix< T > > &q0, const std::vector< Matrix< T > > &q1, const std::vector< Matrix< T > > &q2)
	Stationary distribution of a level-dependent QBD given its rate matrices (ldqbd_pi.m).
template<class T>
LdqbdResult< T >	line::mam::ldqbd (const std::vector< Matrix< T > > &q0, const std::vector< Matrix< T > > &q1, const std::vector< Matrix< T > > &q2)
	Solve a level-dependent QBD: rate matrices and stationary law (ldqbd.m).

Detailed Description

Level-dependent QBD processes with finitely many levels: the rate matrices R^(n) by the backward matrix continued fraction, and the stationary distribution.

Templated port of matlab/src/api/mam/ldqbd_R.m, ldqbd_pi.m and ldqbd.m, which implement Algorithms 1 and 3 of T. Phung-Duc, H. Masuyama, S. Kasahara, Y. Takahashi, "A Simple Algorithm for the Rate Matrices of Level-Dependent QBD Processes", QTNA 2010.

The generator is block tridiagonal over levels 0..N,

Q = | Q1(0)  Q0(0)   0      ...            |
    | Q2(1)  Q1(1)  Q0(1)   ...            |
    |  0     Q2(2)  Q1(2)   ...            |
    | ...                   Q2(N)   Q1(N)  |

with Q0[n] the up-block from level n, Q1[n] the local block at level n and Q2[n] the down-block from level n. The dimensions may vary with the level, so R^(n) is (order of level n-1) x (order of level n). The recursion runs downwards from the top level,

R^(N) = Q0[N-1] (-Q1[N])^-1
R^(n) = Q0[n-1] (-Q1[n] - R^(n+1) Q2[n+1])^-1,   n = N-1, ..., 1

and the stationary vectors follow from the level-0 balance pi_0 (Q1[0] + R^(1) Q2[1]) = 0 and pi_n = pi_{n-1} R^(n), normalized over all levels.

INDEXING. The reference passes Q0, Q1, Q2 as MATLAB cell arrays with three different origins: Q0{k} is Q0^(k-1), Q1{k} is Q1^(k-1) and Q2{k} is Q2^(k). The port takes three std::vector, all indexed by LEVEL: q0[n] for n = 0..N-1, q1[n] for n = 0..N, q2[n] for n = 1..N, with q2[0] required to be present but unused (a level-0 down-block does not exist), so q2.size() == q1.size(). This removes the origin mismatch that makes the MATLAB call sites hard to read; the tests cross-check against MATLAB with the shift applied explicitly.

ARITHMETIC. The recursion is FINITE – N matrix inverses and N products – with no tolerance and no iteration, so ldqbd_R and ldqbd_pi instantiate at Rational and return the rate matrices and the stationary law as exact fractions. That is unusual for a QBD: the level-independent case needs cyclic reduction and can only ever be approximate (see qbd_r.h), whereas a level-dependent chain with a finite top level is a finite linear algebra problem. The one exception is the singular fallback, see below.

SINGULAR LEVELS. The reference tests abs(det(U)) > 1e-14 and falls back to Q0 * pinv(U) when it fails. An absolute determinant threshold is not meaningful at exact arithmetic (a rational matrix is singular or it is not, and det scales like the n-th power of the entries), so the port splits:

inexact T keeps the reference behaviour, |det(U)| <= 1e-14 routes to the LAPACK pseudo-inverse of util/svd.h, in double, with the result lifted back into T;
exact T tests for an exactly singular U and raises NumericError naming the level, because a pseudo-inverse of a rational matrix is not rational and silently returning a rounded one would make the "exact" label false.

Definition in file ldqbd.h.

Classes

Namespaces

Functions

Detailed Description