quadrature-free closed forms for the LCD / modified Cramer-von-Mises distance of an even-dimensional standard normal target More...

#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_expint.h>
#include <cassert>
#include <cmath>
#include <cstddef>

Data Structures
struct	LcdDeltaBIntegrandParams
	parameters of the raw b-integrand used by the small-c fallback More...

Functions
double	lcd_ei (double x)
	exponential integral Ei(x) for negative arguments

double	lcd_a_n (size_t N, double b)
	D1, the sample-independent term of the mCvM distance.

double	lcd_c_repulsion (double b, double c)
	C(b, c), the per-pair repulsion term of D3.

double	lcd_delta_b_integrand (double b, void *params)
	b^(2k+1) * (1 + 2b^2)^(-q) * exp(-c / (2 (1 + 2b^2)))

double	lcd_delta_b0 (size_t d, double u, double eULeading, double eU, double e0, double dEi, double c)
	dB_{0,d}(c) = B_{0,d}(bMax, c) - B_{0,d}(0, c)

double	lcd_delta_bkk (size_t k, double bMax, double c)
	dBkk(k, c) = B_{k,k}(bMax, c) - B_{k,k}(0, c), the D2 term

double	lcd_delta_bkk1 (size_t k, double bMax, double c)
	dBkk1(k, c) = B_{k,k+1}(bMax, c) - B_{k,k+1}(0, c), the D2 gradient

double	lcd_delta_bkk_closed (size_t k, double bMax, double c, size_t shift, bool reduced=false)
	closed form of dBkk / dBkk1

double	lcd_delta_bkk_zero (size_t k, double bMax)
	dBkk(k, 0) = integral{0 to bMax} b^(2k+1) / (1 + 2b^2)^k db

double	lcd_delta_bkk_zero_minus_leading (size_t k, double bMax)
	lcd_delta_bkk_zero(k, bMax) - pow(2, -k) * U / 4

double	lcd_delta_bkk_reduced (size_t k, double bMax, double c)
	dBkk(k, c) - dBkk(k, 0), the cancellation-free attraction term

double	lcd_delta_bkk_impl (size_t k, double bMax, double c, size_t shift)
	shared body of lcd_delta_bkk / lcd_delta_bkk1

double	lcd_reported_offset (size_t N, double bMax, double weightSum)
	the constant that turns the reduced objective back into the true distance

Detailed Description

quadrature-free closed forms for the LCD / modified Cramer-von-Mises distance of an even-dimensional standard normal target

The mCvM distance of a Dirac mixture with samples s_i and weights w_i against N(0, I) decomposes into

D = D1 - 2 * D2 + D3

where D1 is sample-independent, D2 is the attraction term (samples towards the mean) and D3 is the repulsion term (samples away from each other).

Note: Preconditions for the D2 helpers (lcd_delta_bkk, lcd_delta_bkk1): the target must be the STANDARD normal N(0, I) and N = 2k must be EVEN. lcd_a_n and lcd_c_repulsion are valid for all N.

Function Documentation

◆ lcd_a_n()

double lcd_a_n	(	size_t	N,
		double	b
	)

inline

D1, the sample-independent term of the mCvM distance.

Evaluates A_N(b) = integral_0^b x^(N+1) / (1 + x^2)^(N/2) dx via the recursion

A_1(b) = 0.5 * (b * sqrt(1 + b^2) - asinh(b))

A_2(b) = 0.5 * (b^2 - log(1 + b^2))

A_n(b) = (n * A_{n-2}(b) - b^n / (1 + b^2)^((n-2)/2)) / (n - 2)

applied iteratively. Valid for all N, odd and even.

Parameters

N	dimension, must be >= 1
b	upper integration bound, must be >= 0

Returns: A_N(b)

◆ lcd_c_repulsion()

double lcd_c_repulsion	(	double	b,
		double	c
	)

inline

C(b, c), the per-pair repulsion term of D3.

C(b, c) = integral_0^b x * exp(-c / (4 x^2)) dx

    = (b^2 / 2) * exp(-c / (4 b^2)) + (c / 8) * Ei(-c / (4 b^2))

Parameters

b	upper integration bound, must be > 0
c	squared distance between the two samples, must be >= 0

Returns: C(b, c); C(b, 0) = b^2 / 2

◆ lcd_delta_b0()

double lcd_delta_b0	(	size_t	d,
		double	u,
		double	eULeading,
		double	eU,
		double	e0,
		double	dEi,
		double	c
	)

inline

dB_{0,d}(c) = B_{0,d}(bMax, c) - B_{0,d}(0, c)

Parameters

d	index of the base case
u	1 + 2 * bMax^2
eULeading	exp(-c / (2u)), or expm1(-c / (2u)) for the reduced form
eU	exp(-c / (2u))
e0	exp(-c / 2)
dEi	Ei(-c / (2u)) - Ei(-c / 2)
c	squared norm of the sample, must be > 0

Returns: dB_{0,d}(c), less u/4 for d == 0 if eULeading is the expm1 form

◆ lcd_delta_b_integrand()

double lcd_delta_b_integrand	(	double	b,
		void *	params
	)

inline

b^(2k+1) * (1 + 2b^2)^(-q) * exp(-c / (2 (1 + 2b^2)))

With reduced set, exp is replaced by expm1, so the integral evaluates dBkk(c) - dBkk(0) directly and never forms the O(bMax^2) intermediate - the same trick gm_to_dirac_short::calculateP2 uses.

Parameters

b	integration variable
params	LcdDeltaBIntegrandParams

Returns: value of the integrand at b

◆ lcd_delta_bkk()

double lcd_delta_bkk	(	size_t	k,
		double	bMax,
		double	c
	)

inline

dBkk(k, c) = B_{k,k}(bMax, c) - B_{k,k}(0, c), the D2 term

D2 = pow(2, k) * sum_i w_i * lcd_delta_bkk(k, bMax, c_i), with c_i = ||s_i||^2.

Parameters

k	half the dimension, N = 2k; must be >= 1
bMax	upper integration bound, must be > 0
c	squared norm of the sample, must be >= 0

Returns: dBkk(k, c)

Note: standard-normal target and even N only

◆ lcd_delta_bkk1()

double lcd_delta_bkk1	(	size_t	k,
		double	bMax,
		double	c
	)

inline

dBkk1(k, c) = B_{k,k+1}(bMax, c) - B_{k,k+1}(0, c), the D2 gradient

Parameters

k	half the dimension, N = 2k; must be >= 1
bMax	upper integration bound, must be > 0
c	squared norm of the sample, must be >= 0

Returns: dBkk1(k, c)

Note: standard-normal target and even N only

◆ lcd_delta_bkk_closed()

double lcd_delta_bkk_closed	(	size_t	k,
		double	bMax,
		double	c,
		size_t	shift,
		bool	reduced = `false`
	)

inline

closed form of dBkk / dBkk1

Parameters

k	half the dimension, N = 2k
bMax	upper integration bound
c	squared norm of the sample, must be > 0
shift	0 for dBkk (uses dB_{0,j}), 1 for dBkk1 (uses dB_{0,j+1})

Returns: pow(2,-k) * sum_{j=0..k} (-1)^j * binom(k,j) * dB_{0,j+shift}(c)

◆ lcd_delta_bkk_impl()

double lcd_delta_bkk_impl	(	size_t	k,
		double	bMax,
		double	c,
		size_t	shift
	)

inline

shared body of lcd_delta_bkk / lcd_delta_bkk1

Parameters

k	half the dimension, N = 2k
bMax	upper integration bound
c	squared norm of the sample
shift	0 for dBkk (uses dB_{0,j}), 1 for dBkk1 (uses dB_{0,j+1})

Returns: dBkk (shift 0) or dBkk1 (shift 1)

◆ lcd_delta_bkk_reduced()

double lcd_delta_bkk_reduced	(	size_t	k,
		double	bMax,
		double	c
	)

inline

dBkk(k, c) - dBkk(k, 0), the cancellation-free attraction term

Equal to integral{0 to bMax} b^(2k+1)/(1+2b^2)^k * expm1(-c/(2(1+2b^2))) db.

The constant removed per sample is lcd_delta_bkk_zero(k, bMax); D2 as a whole loses pow(2, k) * (sum_i w_i) * lcd_delta_bkk_zero(k, bMax).

Parameters

k	half the dimension, N = 2k; must be >= 1
bMax	upper integration bound, must be > 0
c	squared norm of the sample, must be >= 0

Returns: dBkk(k, c) - dBkk(k, 0)

Note: standard-normal target and even N only

◆ lcd_delta_bkk_zero()

double lcd_delta_bkk_zero	(	size_t	k,
		double	bMax
	)

inline

dBkk(k, 0) = integral{0 to bMax} b^(2k+1) / (1 + 2b^2)^k db

The c = 0 value of the attraction term, in elementary closed form. With u = 1 + 2b^2 the substitution gives

Z = 2^-(k+2) * integral{1 to U} (1 - 1/u)^k du

= 2^-(k+2) * [ (U-1) - k*ln(U)

sum_{j=2..k} (-1)^j C(k,j) (U^(1-j) - 1)/(1-j) ]

This is the constant lcd_delta_bkk_reduced() removes. It is O(bMax^2) by construction, so it belongs in the reported offset, never inside the optimizer's objective.

Parameters

k	half the dimension, N = 2k; must be >= 1
bMax	upper integration bound, must be > 0

Returns: dBkk(k, 0)

◆ lcd_delta_bkk_zero_minus_leading()

double lcd_delta_bkk_zero_minus_leading	(	size_t	k,
		double	bMax
	)

inline

lcd_delta_bkk_zero(k, bMax) - pow(2, -k) * U / 4

The same expression with (U - 1) replaced by -1, i.e. with the O(bMax^2) leading term removed. The result is O(ln U), so subtracting it from the expm1 form of the closed dBkk introduces no cancellation.

Parameters

k	half the dimension, N = 2k; must be >= 1
bMax	upper integration bound, must be > 0

Returns: dBkk(k, 0) - pow(2, -k) * (1 + 2 bMax^2) / 4

◆ lcd_ei()

double lcd_ei ( double x )

inline

exponential integral Ei(x) for negative arguments

Parameters

x	argument, must be strictly negative

Returns: Ei(x), or exactly 0.0 in case of underflow, domain error or overflow

◆ lcd_reported_offset()

double lcd_reported_offset	(	size_t	N,
		double	bMax,
		double	weightSum
	)

inline

the constant that turns the reduced objective back into the true distance

The optimizer minimises a reduced objective in which every O(bMax^2) piece has been removed analytically:

f_reduced = -2 * pow(2,k) * sum_i w_i * lcd_delta_bkk_reduced(k, bMax, c_i)

sum_ij w_i w_j * C(bMax, T_ij)

and the true distance is D = f_reduced + K, with

K = A_N(bMax)

pow(2, k+1) * W * lcd_delta_bkk_zero(k, bMax)
W^2 * bMax^2 / 2, W = sum_i w_i

collecting the D1 term, the per-sample attraction constant and the per-pair repulsion constant respectively.

Parameters

N	dimension, must be even and >= 2
bMax	upper integration bound, must be > 0
weightSum	sum of the Dirac weights, normally 1

Returns: K

Data Structures

Functions

Detailed Description

Function Documentation

◆ lcd_a_n()

◆ lcd_c_repulsion()

◆ lcd_delta_b0()

◆ lcd_delta_b_integrand()

◆ lcd_delta_bkk()

◆ lcd_delta_bkk1()

◆ lcd_delta_bkk_closed()

◆ lcd_delta_bkk_impl()

◆ lcd_delta_bkk_reduced()

◆ lcd_delta_bkk_zero()

◆ lcd_delta_bkk_zero_minus_leading()

◆ lcd_ei()

◆ lcd_reported_offset()