Dev

[balaji1810]

Closed-form LCD approximation for even-dimensional Gaussians

1. The objective

With samples $s_1 \dots s_L$, weights $w_i$, $W = \sum_i w_i$, and writing $c_i = | s_i |^2$ and $T_{ij} = | s_i - s_j |^2$, the distance splits into

$$ D = D_1 - 2 D_2 + D_3 $$

$$ D_1 = \int_0^{b} \frac{x^{N+1}}{(1+x^2)^{N/2}} dx $$

$$ D_2 = \int_0^{b} x \left( \frac{2x^2}{1+2x^2} \right)^{N/2} \sum_i w_i \exp\left( \frac{-c_i}{2(1+2x^2)} \right) dx $$

$$ D_3 = \int_0^{b} x \sum_{i,j} w_i w_j \exp\left( \frac{-T_{ij}}{4x^2} \right) dx $$

$D_1$ is sample-independent, $D_2$ is attraction, $D_3$ is repulsion. All in units of $\pi^{N/2}$, the same normalisation the quadrature path uses.

The even - $N$ restriction is what makes $D_2$ integrable in closed form: with $N = 2k$ the fractional power $(\cdot)^{N/2}$ becomes the integer power $(\cdot)^k$, and the integral collapses to a finite binomial sum.

---

2. The closed forms

All in lib/gm_to_dirac_even/lcd_even_closed_form.h. Throughout, $U = 1 + 2b^2$, $e_U = \exp(-c/(2U))$, $e_0 = \exp(-c/2)$, and $\Delta E = \text{Ei}(-c/(2U)) - \text{Ei}(-c/2)$.

lcd_a_n(N, b) - the $D_1$ integral

$$ A_N(b) = \int_0^{b} \frac{x^{N+1}}{(1+x^2)^{N/2}} dx $$

$$ A_1(b) = \frac{1}{2}\left( b\sqrt{1+b^2} - \text{asinh}(b) \right), \qquad A_2(b) = \frac{1}{2}\left( b^2 - \ln(1+b^2) \right) $$

$$ A_n(b) = \frac{ n A_{n-2}(b) - b^n (1+b^2)^{-(n-2)/2} }{ n-2 }, \qquad n \ge 3 $$

Valid for all $N$, odd and even.

lcd_c_repulsion(b, c) - the per-pair $D_3$ term

$$ C(b,c) = \int_0^{b} x \exp\left( \frac{-c}{4x^2} \right) dx = \frac{b^2}{2} \exp\left( \frac{-c}{4b^2} \right) + \frac{c}{8} \text{Ei}\left( \frac{-c}{4b^2} \right) $$

with $C(b,0) = b^2/2$.

lcd_delta_bkk(k, b, c) and lcd_delta_bkk1(k, b, c) - the per-sample $D_2$ terms

$$ B_{k,q}(c) = \int_0^{b} \frac{x^{2k+1}}{(1+2x^2)^{q}} \exp\left( \frac{-c}{2(1+2x^2)} \right) dx $$

lcd_delta_bkk is $q = k$ (objective), lcd_delta_bkk1 is $q = k+1$ (gradient). Substituting $u = 1+2x^2$ and expanding binomially gives one formula for both:

$$ B_{k,q}(c) = 2^{-(k+2)} \sum_{j=0}^{k} (-1)^j \binom{k}{j} I_{j + q - k} $$

$$ I_d = \int_1^{U} u^{-d} \exp\left( \frac{-c}{2u} \right) du $$

$$ I_0 = U e_U - e_0 + \frac{c}{2} \Delta E, \qquad I_1 = -\Delta E, \qquad I_d = \frac{2}{c}\left( \frac{e_U}{U^{d-2}} - e_0 \right) + \frac{2(d-2)}{c} I_{d-1} $$

The recursion terminates at $d = 2$ because the coefficient carries a factor $(d-2)$. In the code this is lcd_delta_b0, which returns $I_d / 4$; lcd_delta_bkk_closed walks the binomial sum with the shift parameter selecting $q = k$ or $q = k+1$, so both functions share one implementation.

lcd_delta_bkk_zero(k, b) - the $c = 0$ constant

$$ Z(k,b) = B_{k,k}(0) = \int_0^{b} \frac{x^{2k+1}}{(1+2x^2)^{k}} dx = 2^{-(k+2)} \left[ (U-1) - k \ln U + \sum_{j=2}^{k} (-1)^j \binom{k}{j} \frac{U^{1-j}-1}{1-j} \right] $$

This is $O(b^2)$ by construction, it is the constant that gets pulled out of the objective.

lcd_delta_bkk_reduced(k, b, c) - the cancellation-free attraction term

$$ B_{k,k}(c) - Z(k,b) = \int_0^{b} \frac{x^{2k+1}}{(1+2x^2)^{k}} \text{expm1}\left( \frac{-c}{2(1+2x^2)} \right) dx $$

Computed with expm1 in the leading $I_0$ term so the $O(b^2)$ piece $2^{-k}U/4$ is never formed in floating point; lcd_delta_bkk_zero_minus_leading removes the remaining $O(\ln U)$ part.

lcd_reported_offset(N, b, W) - the constant that restores the true distance

$$ K = A_N(b) - 2^{k+1} W Z(k,b) + \frac{W^2 b^2}{2} $$

collecting, in order, the $D_1$ term, the per-sample attraction constant and the per-pair repulsion constant. Each piece is $O(b^2)$, but the leading parts combine to $\frac{1}{2}(1-W)^2 b^2$, which vanishes for normalised weights.

For $N=2$ the branch is algebraically rearranged so that cancellation is exact rather than asymptotic:

$$ K = \frac{b^2 (1-W)^2}{2} + \frac{W}{2}\ln\left( 1 + \frac{b^2}{1+b^2} \right) - \frac{1-W}{2}\ln(1+b^2) $$

---

3. How they are assembled

Objective (calculateD2 + calculateD3, both in the .tpp). The optimizer minimises a reduced objective in which every $O(b^2)$ constant has been removed analytically:

$$ f = -2 \cdot 2^{k} \sum_i w_i \left[ B_{k,k}(c_i) - Z(k,b) \right] + \sum_{i,j} w_i w_j \left[ C(b,T_{ij}) - \frac{b^2}{2} \right] $$

The first sum is lcd_delta_bkk_reduced; the second is written inline as $\frac{b^2}{2}\text{expm1}(z) + \frac{c}{8}\text{Ei}(z)$, so coincident and self-pairs ($T_{ij}=0$) contribute exactly zero instead of $b^2/2$ each.

Gradient. Since $c_q = | s_q |^2$ and $\partial B_{k,k} / \partial c = -\frac{1}{2} B_{k,k+1}$:

$$ \frac{\partial f}{\partial s_{q,d}} = 2^{k+1} w_q B_{k,k+1}(c_q) s_{q,d} + \frac{1}{2} w_q \sum_{j} w_j \text{Ei}\left( \frac{-T_{qj}}{4b^2} \right) \left( s_{q,d} - s_{j,d} \right) $$

The $D_3$ part is unusually clean because the three pieces of $\partial C / \partial c$ collapse to $\frac{1}{8}\text{Ei}(z)$, so one Ei evaluation per pair serves both objective and gradient.

Reported distance. GMToDiracEvenOptimizationParams computes $K$ once in its constructor, and

$$ D = f + K $$

modified_van_mises_distance_sq always reports this true $D$. The minimizer carries $K$ too by default (includeD1InObjective = true) - since $K$ is $O(1)$ rather than $O(b^2)$, adding it costs no $x$-dependent precision and keeps ApproximateOptions::ftolRel (default 1e-10) meaningful instead of firing against a large constant. This follows the pattern of #4, already on main.