Saddle-point approximation of the cache normalizing constant. More...

#include <cmath>
#include <cstddef>
#include <limits>
#include <vector>
#include "line/api/cache/cache_erec.h"
#include "line/api/cache/cache_xi_iter.h"
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/matrix.h"

Include dependency graph for cache_spm.h:

Classes
struct	line::cache::CacheSpmResult< T >
	Return value of cache_spm, mirroring [Z,lZ,xi]. More...

Namespaces
namespace	line
namespace	line::cache

Functions
template<class T>
CacheSpmResult< T >	line::cache::cache_spm (const Matrix< T > &gamma_in, const std::vector< int > &m)
	Saddle-point approximation of the cache normalizing constant.

Detailed Description

Saddle-point approximation of the cache normalizing constant.

Templated port of matlab/src/api/cache/cache_spm.m, cross-checked against jar/src/main/java/jline/api/cache/Cache_spm.java.

The constant E(gamma,m) of cache_erec is the coefficient extraction

E = [prod_l z_l^{m_l}] prod_k (1 + sum_l gamma(k,l) z_l) * prod_l m_l!,

evaluated here by a multidimensional saddle point at the multipliers xi from cache_xi_iter. With S(k) = sum_l gamma(k,l) xi(l) and

phi = sum_k log(1 + S(k)) - sum_l m(l) log(xi(l)), C(j,l) = delta(j,l) sum_k gamma(k,j)/(1+S(k))

xi(j) sum_k gamma(k,j) gamma(k,l)/(1+S(k))^2,

the Gaussian integral around the saddle gives

log E ~ phi + sum_l log(m_l!) - (h/2) log(2 pi) - (1/2) sum_l log xi(l)

(1/2) log det C.

TWO BOUNDARIES, both fixed in all four codebases on 2026-08-29; the notes below say what the behaviour used to be, because saved results predating the fix carry it.

n == sum(m): every item is cached, so the capacity equations force every multiplier to infinity and there is no interior saddle. Z comes from the exact cache_erec and xi is reported as +infinity, its limit, WITHOUT running the iteration. Before the fix MATLAB computed Z exactly but then called cache_xi_iter for the third return value and hung (cache_spm([.5 .25;.4 .2;.3 .15],[2 1]) did not return in two minutes); the JAR had no fallback at all and ran the same non-terminating loop for every value; this port raised NumericError from cache_xi_iter's sweep cap, losing the exact Z with it.

m_l == 0: list l has xi(l)=0, which is a boundary of the Laplace integral rather than a direction of it, so list l is dropped before the saddle solve and h shrinks with it. Dropping is exact: setting z_l=0 in the generating function removes list l from E(m), and prod_l m_l! is unchanged because 0!=1. Before the fix no codebase dropped it – only all-zero ROWS of gamma were filtered, never zero-capacity COLUMNS – so the bisection floored xi(l) at ~2^-50 and the -(1/2) sum_l log xi(l) prefactor gained ~+17 per empty list, silently: n=12 items at gamma=(0.8,0.6,0.4) and m=(0,3,3) returned lZ=24.814 against an exact 9.127, a factor of 6.5e8. cache_prob_spm reaches this on any list with m_l==1, since it evaluates E at oner(m,l). All m_l zero is the empty cache: Z=1, lZ=0, xi=0.

ARITHMETIC: log, exp and sqrt throughout, plus the tolerance-stopped cache_xi_iter, so transcendental arithmetic is required.

MATLAB's lZ=real(lZ) discards the imaginary part that appears when det C or some xi(l) is negative; that is implemented here as taking the modulus inside the logarithm, which is the same real branch. Z is returned on the same branch (MATLAB would return a complex Z there).

The item count is gamma's ROW count in every codebase. MATLAB used to read n = length(gamma), which is max(n,h), and used it both for the n == mt degenerate test and as the loop bound over gamma's rows, reading past the end whenever h > n; that was corrected to size(gamma,1) alongside the two fixes above. cache_erec still carries the same length() misuse.

Definition in file cache_spm.h.

Classes

Namespaces

Functions

Detailed Description