pfqn_simplex.h File Reference

Shared machinery for closures of the simplex factor of the McKenna-Mitra integral, used by pfqn_aghq.h. More...

#include <cmath>
#include <cstddef>
#include <vector>
#include "line/api/pfqn/pfqn_asympt_common.h"
#include "line/api/pfqn/pfqn_le.h"
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/matrix.h"

Include dependency graph for pfqn_simplex.h:

Go to the source code of this file.

Classes
struct	line::pfqn::simplex::Radial< T >
	log J(c) and the moments of the tilted law of the radius. More...
struct	line::pfqn::simplex::Mode< T >
	Mode, curvature and log-integrand at the mode of the simplex factor. More...

Namespaces
namespace	line
namespace	line::pfqn
namespace	line::pfqn::simplex

Functions
template<class T>
void	line::pfqn::simplex::sym_eig (Matrix< T > A, std::vector< T > &d, Matrix< T > &V)
	Eigenvalues and eigenvectors of a symmetric matrix by cyclic Jacobi, ascending.
template<class T>
void	line::pfqn::simplex::golub_welsch (const std::vector< T > &off, const T &mu0, std::vector< T > &x, std::vector< T > &w)
	Nodes and weights of the Gauss rule with the given Jacobi off-diagonal.
template<class T>
void	line::pfqn::simplex::gauss_legendre (std::size_t n, std::vector< T > &x, std::vector< T > &w)
	N-point Gauss-Legendre rule on [-1,1].
template<class T>
void	line::pfqn::simplex::gauss_hermite (std::size_t q, std::vector< T > &z, std::vector< T > &w)
	Q-point Gauss-Hermite rule of the probabilists' weight exp(-z^2/2).
template<class T>
T	line::pfqn::simplex::radial_logf (const T &t, const std::vector< T > &c, const std::vector< T > &N, const std::vector< T > &Z, std::size_t M)
	Log-integrand of the radial integral in t = log v, Jacobian included.
template<class T>
Radial< T >	line::pfqn::simplex::radial (const std::vector< T > &c, const std::vector< T > &N, const std::vector< T > &Z, std::size_t M)
	log J(c) = log int_0^inf exp(-v) v^(M-1) prod_r (Z_r + v c_r)^N_r dv, plus the moments of the tilted law of v that the simplex derivatives need: G_r = E[T_r], vbar = E[v] and Lam = cov(T) - diag(E[T^2]/N) = grad^2_c log J, with T_r(v) = N_r v/(Z_r + v c_r).
template<class T>
Mode< T >	line::pfqn::simplex::simplex_mode (const Matrix< T > &L, const std::vector< T > &N, const std::vector< T > &Z)
	Mode and curvature of h(w) = log J(L'x(w)) + sum_i log x_i with J the exact radial integral.
template<class T>
std::vector< T >	line::pfqn::simplex::softmax_gauge (const std::vector< T > &w)
	softmax of [w; 0], the logistic parametrisation of the simplex with gauge w_M = 0.
template<class T, class F>
T	line::pfqn::simplex::aghq_rule (const F &h, const std::vector< T > &w0, const T &h0, const Matrix< T > &A, std::size_t q, std::size_t d)
	Log of the tensor Gauss-Hermite sum, accumulated with a running maximum; the det(A)^(-1/2) of the rule is applied by the caller.

Detailed Description

Shared machinery for closures of the simplex factor of the McKenna-Mitra integral, used by pfqn_aghq.h.

At Z = 0 the integrand is homogeneous of degree sum(N), so y = v*x separates and the radius integrates exactly to gamma(N+M), leaving an integral over the unit simplex in which ALL of the error of the logistic expansion lives. With Z > 0 that factorisation is gone: the radius cannot be marginalised and is integrated numerically here rather than closed, leaving the same M-1 simplex directions to a closure.

ARITHMETIC. Everything here is a Laplace-type approximation or a quadrature of a transcendental integrand, so each entry point is gated on num_traits<T>::has_transcendental exactly as pfqn_le is.

EIGENSOLVER. util/eig.h is LAPACK and double-only, while this family is templated on T, so a cyclic Jacobi eigensolver is carried here instead. It is used for the Golub-Welsch construction of the quadrature nodes and for the principal-axis frame of the adaptive Gauss-Hermite rule. Jacobi is chosen over a tridiagonal QL because the matrices are small, it needs no shift strategy, and it is symmetric-exact by construction.

Definition in file pfqn_simplex.h.