mfq_solve.h File Reference

Core of the Markovian fluid queue: the fundamental matrices Psi, K, U and the matrix-exponential stationary solution of a general fluid model. More...

#include <cmath>
#include <cstddef>
#include <vector>
#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 dependency graph for mfq_solve.h:

Go to the source code of this file.

Classes
struct	line::mam::FluidFundamental< T >
	Psi, K and U of a fluid queue with drifts normalized to +-1. More...
struct	line::mam::GeneralFluidSolution< T >
	Stationary matrix-exponential solution of a general Markovian fluid model. More...

Namespaces
namespace	line
namespace	line::mam

Enumerations
enum class	line::mam::RiccatiMethod { line::mam::ADDA , line::mam::SDA }
	Which doubling iteration to run for Psi. More...

Functions
template<class T>
FluidFundamental< T >	line::mam::mfq_fundamental (const Matrix< T > &Fpp, const Matrix< T > &Fpm, const Matrix< T > &Fmp, const Matrix< T > &Fmm, const T &precision, unsigned maxNumIt, RiccatiMethod method)
	Psi, K and U of a fluid queue whose drifts have been normalized to +-1.
template<class T>
GeneralFluidSolution< T >	line::mam::mfq_general_solve (const Matrix< T > &Q, const Matrix< T > &R, const Matrix< T > &Q0, const T &prec)
	Stationary law of a general Markovian fluid model, pi(x) = ini exp(K x) clo above level zero plus the point mass mass0 at zero.
template<class T>
GeneralFluidSolution< T >	line::mam::mfq_general_solve (const Matrix< T > &Q, const Matrix< T > &R)
	mfq_general_solve with the regular boundary and the BUTools default prec = 1e-14.
template<class T>
Matrix< T >	line::mam::mfq_transform_to_ones (const std::vector< T > &v)
	The similarity transformation B with B v = e, for a non-negative column vector v.

Detailed Description

Core of the Markovian fluid queue: the fundamental matrices Psi, K, U and the matrix-exponential stationary solution of a general fluid model.

These are the two BUTools routines that every mfq_* entry point of matlab/src/api/mam/ ultimately calls (FluidFundamentalMatrices and GeneralFluidSolve). They carry no LINE-level name of their own, so they live here as shared infrastructure rather than as a ported API function, and mfq_sojourn.h and mfq_fluflu_sojourn.h are the thin entry points over them.

THE MODEL. A background chain with generator Q modulates a fluid level whose drift in state i is R(i,i), of any sign, zero included. Writing the states in the order (zero drift, positive drift, negative drift), censoring out the zero-drift states and rescaling time by |1/R| turns the model into a fluid queue with drifts +-1, whose first passage matrix Psi solves the Riccati equation

Fpm + Fpp Psi + Psi Fmm + Psi Fmp Psi = 0.

Psi(i,j) is the probability that, starting in an up state i at a level, the process returns to that level in the down state j. From it, K = Fpp + Psi Fmp generates the stationary density and U = Fmm + Fmp Psi is the generator of the chain seen at level zero. The stationary law is then a point mass at level zero plus the density pi(x) = ini exp(K x) clo.

THE RICCATI SOLVER. ADDA (Wang, Wang and Li 2011), the BUTools default: an alternating-directional doubling iteration whose per-step cost is four inverses and whose error SQUARES each step, so the 1e-14 default is reached in a few tens of iterations rather than the hundreds a linearly convergent fixed point would need. SDA is the same iteration with a common shift and is offered because it is the form quoted in most of the literature; the two differ only in the shift and in the scaling of E and F. The reference's third option, cyclic reduction, is NOT ported: it is an alternative with the same convergence order and no accuracy advantage on this problem, and it would duplicate the whole solver for nothing. Callers asking for it get an UnsupportedError naming it, rather than a silent substitution.

WHERE THE PORT IS STRICTER THAN THE REFERENCE. GeneralFluidSolve censors the zero-drift states through pinv(-Qv00), a PSEUDO-inverse. That matrix is singular exactly when the zero-drift states form a closed set, i.e. when the fluid can be trapped at a constant level forever and the model has no stationary fluid law to report. The reference then returns whatever the pseudo-inverse yields, without comment. The port uses the ordinary inverse and raises NumericError naming the condition, because a pseudo-inverse of a singular censoring operator is not the answer to a different question, it is an answer to no question. On every non-degenerate model the two agree exactly, the pseudo-inverse of a non-singular matrix being its inverse.

The two overdetermined normalizations that the reference solves with MATLAB's backslash (an (n+1) x n system that is consistent by construction) are solved here through the normal equations, which return the same vector for a consistent full-column-rank system and need no least-squares factorization.

ARITHMETIC. Gated on num_traits<T>::has_transcendental: ADDA terminates on a tolerance and its scaling step takes a square root. Everything else is finite exact linear algebra.

Definition in file mfq_solve.h.