ctmc_classical.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_API_MC_CTMC_CLASSICAL_H
#define LINE_API_MC_CTMC_CLASSICAL_H
 
/**
 * @file
 * @ingroup api_mc
 * Classical stationary iterative solvers for pi Q = 0, sum(pi) = 1: the power
 * method on the uniformised DTMC, Jacobi, Gauss-Seidel and SOR (Stewart 1994,
 * ch. 3; Haverkort 1998, section 15.1.3).
 *
 * Templated port of matlab/src/api/mc/ctmc_power.m, ctmc_jacobi.m,
 * ctmc_gaussseidel.m and ctmc_sor.m, which share the loop of
 * matlab/src/api/mc/private/stationary_sweep_iterate.m. The JAR twin is
 * jline.api.mc.StationarySweep and the Python twin
 * line_solver.api.mc.classical. All four codebases run the same loop:
 *
 *   x = p0 / sum(p0) (uniform by default); for k = 0..maxiter:
 *     relres = ||x Q||_1 / max|q_ii|; stop with flag 0 if relres <= tol,
 *     flag 1 if k == maxiter; otherwise sweep, renormalise to sum one,
 *     and stop with flag 4 if the sweep sum is zero or not finite.
 *
 * Sweeps, with A = Q' split as D + L + U:
 *   power        x <- x + (x Q)'/q, q = 1.1 max|q_ii| (aperiodic by the
 *                self-loop of P = I + Q/q, so periodic chains converge)
 *   jacobi       x <- -D^{-1} (L + U) x
 *   gaussseidel  x <- (D + L)^{-1} (-U x)
 *   sor          x <- (D + w L)^{-1} ((1 - w) D - w U) x
 *
 * Flag 2 reports a zero diagonal (an absorbing state) for every method but the
 * power method, which needs no diagonal inverse. Jacobi does not converge when
 * the embedded jump chain is periodic, for instance on any birth-death chain;
 * it then reports flag 1. On a reducible generator the limit depends on p0.
 *
 * SOR with omega <= 0 selects the adaptive rule: start at omega = 1 and every
 * 10 sweeps estimate the Gauss-Seidel spectral radius from the ratio of
 * successive update norms, lambda, then set
 *   mu^2 = (lambda + omega - 1)^2 / (lambda omega^2),
 *   omega = min(2 / (1 + sqrt(1 - mu^2)), 1.95)
 * when that increases omega (Young's formula for consistently ordered
 * matrices, Hageman and Young 1981). If the residual later grows above 100
 * times the reference residual recorded at the last increase, omega reverts to
 * 1 and adaptation stops.
 *
 * Like GMRES these stop on a residual tolerance, so the exact instantiation is
 * refused at compile time: there is nothing exact for the iteration to reach.
 */
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <string>
#include <vector>
 
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/matrix.h"
 
namespace line {
namespace mc {
 
constexpr double CLASSICAL_DEFAULT_TOL = 1e-12;
constexpr long CLASSICAL_DEFAULT_MAXITER = 10000;
 
template <class T>
struct ClassicalResult {
    std::vector<T> p;  ///< stationary distribution estimate, sums to one
    int flag;          ///< 0 converged, 1 iteration cap, 2 zero diagonal, 4 breakdown
    double relres;     ///< ||p Q||_1 / max|q_ii| at return
    long iter;         ///< sweeps performed
    double omega;      ///< relaxation factor at return (1 for Gauss-Seidel, 0 otherwise)
};
 
namespace detail {
 
enum class SweepKind { Power, Jacobi, GaussSeidel, Sor };
 
template <class T>
inline bool classical_isfinite(const T& v) {
    const T a = num_abs(T(v));
    return a < T(a + num_traits<T>::from_int(1));
}
 
template <class T>
T classical_residual_l1(const std::vector<T>& x, const Matrix<T>& Q) {
    const std::size_t n = Q.rows();
    T s = num_traits<T>::from_int(0);
    for (std::size_t j = 0; j < n; ++j) {
        T r = num_traits<T>::from_int(0);
        for (std::size_t i = 0; i < n; ++i) r += x[i] * Q(i, j);
        s += num_abs(r);
    }
    return s;
}
 
template <class T>
ClassicalResult<T> stationary_sweep_iterate(SweepKind kind, const Matrix<T>& Q, double tol, long maxiter,
                                            double omega, const std::vector<T>& p0) {
    static_assert(num_traits<T>::has_transcendental,
                  "the classical iterative solvers stop on a residual tolerance and have no exact "
                  "instantiation; use ctmc_solve for exact arithmetic");
    const std::size_t n = Q.rows();
    if (Q.cols() != n) throw InputError("stationary_sweep_iterate: generator is not square");
    if (n == 0) throw InputError("stationary_sweep_iterate: empty generator");
    if (tol <= 0) tol = CLASSICAL_DEFAULT_TOL;
    if (maxiter <= 0) maxiter = CLASSICAL_DEFAULT_MAXITER;
    const T zero = num_traits<T>::from_int(0);
    const T one = num_traits<T>::from_int(1);
 
    std::vector<T> x(n);
    if (p0.empty()) {
        for (std::size_t i = 0; i < n; ++i) x[i] = one / num_traits<T>::from_int(static_cast<long>(n));
    } else {
        if (p0.size() != n) throw InputError("stationary_sweep_iterate: p0 has the wrong length");
        T s = zero;
        for (const T& v : p0) s += v;
        for (std::size_t i = 0; i < n; ++i) x[i] = p0[i] / s;
    }
 
    const bool adaptive = (kind == SweepKind::Sor) && !(omega > 0);
    if (kind == SweepKind::Sor) {
        if (adaptive) {
            omega = 1.0;
        } else if (!(omega > 0 && omega < 2)) {
            throw InputError("ctmc_sor: relaxation factor omega must lie in (0,2)");
        }
    } else if (kind == SweepKind::GaussSeidel) {
        omega = 1.0;
    } else {
        omega = 0.0;
    }
 
    const double ADAPT_EVERY = 10, OMEGA_CAP = 1.95, DIVERGE_FACTOR = 100;
    T qmax = zero;
    bool zeroDiag = false;
    for (std::size_t i = 0; i < n; ++i) {
        const T a = num_abs(Q(i, i));
        if (a > qmax) qmax = a;
        if (Q(i, i) == zero) zeroDiag = true;
    }
    ClassicalResult<T> res{x, 0, 0.0, 0, omega};
    if (qmax == zero) return res;
    if (kind != SweepKind::Power && zeroDiag) {
        res.flag = 2;
        res.relres = num_traits<T>::to_double(classical_residual_l1(x, Q) / qmax);
        return res;
    }
    const T q = qmax * num_traits<T>::from_rational(11, 10);
 
    bool frozen = false, haveRef = false;
    double dPrev = std::nan(""), dCur = std::nan(""), relresRef = 0.0;
    long iter = 0;
    std::vector<T> r(n), xn(n);
    for (long k = 0; k <= maxiter; ++k) {
        for (std::size_t j = 0; j < n; ++j) {
            T s = zero;
            for (std::size_t i = 0; i < n; ++i) s += x[i] * Q(i, j);
            r[j] = s;
        }
        T rn = zero;
        for (const T& v : r) rn += num_abs(v);
        const double relres = num_traits<T>::to_double(rn / qmax);
        res.relres = relres;
        if (!haveRef) {
            relresRef = relres;
            haveRef = true;
        }
        if (relres <= tol) {
            res.flag = 0;
            res.p = x;
            res.iter = iter;
            res.omega = omega;
            return res;
        }
        if (k == maxiter) break;
        if (adaptive && !frozen && relres > DIVERGE_FACTOR * relresRef) {
            omega = 1.0;
            frozen = true;
        }
        switch (kind) {
            case SweepKind::Power:
                for (std::size_t i = 0; i < n; ++i) xn[i] = x[i] + r[i] / q;
                break;
            case SweepKind::Jacobi:
                for (std::size_t i = 0; i < n; ++i) {
                    T s = zero;
                    for (std::size_t j = 0; j < n; ++j)
                        if (j != i) s += Q(j, i) * x[j];  // (A - D) x with A = Q'
                    xn[i] = -s / Q(i, i);
                }
                break;
            default: {
                // (D + w L) xn = ((1 - w) D - w U) x, by forward substitution.
                const T w = num_traits<T>::from_double(omega);
                for (std::size_t i = 0; i < n; ++i) {
                    T s = zero;
                    for (std::size_t j = 0; j < i; ++j) s += Q(j, i) * xn[j];
                    for (std::size_t j = i + 1; j < n; ++j) s += Q(j, i) * x[j];
                    xn[i] = ((one - w) * Q(i, i) * x[i] - w * s) / Q(i, i);
                }
                break;
            }
        }
        T s = zero;
        bool finite = true;
        for (const T& v : xn) {
            s += v;
            if (!classical_isfinite(v)) finite = false;
        }
        if (!finite || !classical_isfinite(s) || s == zero) {
            res.flag = 4;
            res.p = x;
            res.iter = iter;
            res.omega = omega;
            return res;
        }
        T d = zero;
        for (std::size_t i = 0; i < n; ++i) {
            xn[i] /= s;
            d += num_abs(T(xn[i] - x[i]));
        }
        dPrev = dCur;
        dCur = num_traits<T>::to_double(d);
        x.swap(xn);
        iter = k + 1;
        if (adaptive && !frozen && (iter % static_cast<long>(ADAPT_EVERY)) == 0 && dPrev > 0) {
            const double lambda = dCur / dPrev;
            if (lambda > 0 && lambda < 1) {
                const double mu2 = (lambda + omega - 1) * (lambda + omega - 1) / (lambda * omega * omega);
                if (mu2 >= 0 && mu2 < 1) {
                    const double omegaNew = std::min(2.0 / (1.0 + std::sqrt(1.0 - mu2)), OMEGA_CAP);
                    if (omegaNew > omega) {
                        omega = omegaNew;
                        relresRef = num_traits<T>::to_double(classical_residual_l1(x, Q) / qmax);
                    }
                }
            }
        }
    }
    res.flag = 1;
    res.p = x;
    res.iter = iter;
    res.omega = omega;
    return res;
}
 
}  // namespace detail
 
/**
 * @brief Power method on the uniformised DTMC P = I + Q/q, q = 1.1 max|q_ii|.
 * @param Q generator (row sums zero)
 * @param tol residual tolerance on ||p Q||_1 / max|q_ii| (<= 0 for 1e-12)
 * @param maxiter sweep cap (<= 0 for 10000)
 * @param p0 initial vector, empty for uniform
 */
template <class T>
ClassicalResult<T> ctmc_power(const Matrix<T>& Q, double tol = CLASSICAL_DEFAULT_TOL,
                              long maxiter = CLASSICAL_DEFAULT_MAXITER, const std::vector<T>& p0 = {}) {
    return detail::stationary_sweep_iterate(detail::SweepKind::Power, Q, tol, maxiter, 0.0, p0);
}
 
/** @brief Jacobi iteration on Q' x = 0 with renormalisation (see file comment). */
template <class T>
ClassicalResult<T> ctmc_jacobi(const Matrix<T>& Q, double tol = CLASSICAL_DEFAULT_TOL,
                               long maxiter = CLASSICAL_DEFAULT_MAXITER, const std::vector<T>& p0 = {}) {
    return detail::stationary_sweep_iterate(detail::SweepKind::Jacobi, Q, tol, maxiter, 0.0, p0);
}
 
/** @brief Forward Gauss-Seidel iteration on Q' x = 0 with renormalisation. */
template <class T>
ClassicalResult<T> ctmc_gaussseidel(const Matrix<T>& Q, double tol = CLASSICAL_DEFAULT_TOL,
                                    long maxiter = CLASSICAL_DEFAULT_MAXITER, const std::vector<T>& p0 = {}) {
    return detail::stationary_sweep_iterate(detail::SweepKind::GaussSeidel, Q, tol, maxiter, 1.0, p0);
}
 
/**
 * @brief Successive over-relaxation on Q' x = 0 with renormalisation.
 * @param omega relaxation factor in (0,2); <= 0 selects the adaptive rule
 */
template <class T>
ClassicalResult<T> ctmc_sor(const Matrix<T>& Q, double tol = CLASSICAL_DEFAULT_TOL,
                            long maxiter = CLASSICAL_DEFAULT_MAXITER, double omega = 0.0,
                            const std::vector<T>& p0 = {}) {
    return detail::stationary_sweep_iterate(detail::SweepKind::Sor, Q, tol, maxiter, omega, p0);
}
 
/**
 * Whether a linsolver name selects one of the classical solvers, as accepted
 * by ctmc_solve in all four codebases: "power", "jacobi", "gaussseidel", "sor".
 */
inline bool is_classical_linsolver(const std::string& name) {
    return name == "power" || name == "jacobi" || name == "gaussseidel" || name == "sor";
}
 
/** Run the classical solver named by is_classical_linsolver with defaults. */
template <class T>
ClassicalResult<T> ctmc_classical_by_name(const std::string& name, const Matrix<T>& Q) {
    if (name == "power") return ctmc_power(Q);
    if (name == "jacobi") return ctmc_jacobi(Q);
    if (name == "gaussseidel") return ctmc_gaussseidel(Q);
    if (name == "sor") return ctmc_sor(Q);
    throw InputError("ctmc_classical_by_name: unknown classical solver '" + name + "'");
}
 
}  // namespace mc
}  // namespace line
 
#endif  // LINE_API_MC_CTMC_CLASSICAL_H