# Fixed-label sensor interpolation: exact expectation and tail control

This is the mathematical and numerical justification for Sensors Part 2's expected minimum-norm error. It replaces the large-sensor approximation. No simulated average is used to compute the displayed expectation.

## Scope

Let the training matrix be \(X=ya^T+\sqrt v Z\), with \(Z\) having independent standard Gaussian entries, \(\|y\|^2=n\), \(\|a\|^2=p\), \(v>0\), and \(p\ge n\). Here \(v=\sigma^2\) is the measurement-noise variance used on the Sensors page. Both vectors may be fixed: only their norms matter by rotational invariance. These are the normalized-label comparisons in Part 2. They are different from the unnormalized-label cohort in Part 1.

For \(w=X^\dagger y\), its fresh population risk is

\[
R=(a^Tw-1)^2+v\|w\|^2.
\]

We average this risk over training measurement noise. The formula below does not apply to a fixed realized rack, to unnormalized random labels without another average, or to the underparameterized least-squares case.

## 1. Reduction to one projected Gaussian vector

Rotate sample coordinates so that \(y=\sqrt n e_1\). The last \(n-1\) rows of \(X\) are centered Gaussian noise. Let \(P\) project onto their nullspace, whose dimension is \(k=p-n+1\), almost surely. It is a uniform random subspace independent of the first noise row. Put \(b=Pa\), and express vectors in an orthonormal basis of that subspace. Then

\[
q=\sqrt n b+\sqrt v z,\quad z\sim N(0,I_k),\qquad
w=\frac{\sqrt n q}{\|q\|^2}.
\]

With \(x=q/\sqrt v\), \(\mu=\sqrt{n/v}\,b\), and \(Q=\|x\|^2\),

\[
R=\left(\frac{\mu^Tx}{Q}-1\right)^2+\frac nQ.
\]

Here \(\lambda=\|\mu\|^2=\lambda_0T\), \(\lambda_0=np/v\), and
\(T\sim\operatorname{Beta}(k/2,(n-1)/2)\); when \(n=1\), \(T=1\).

## 2. Tail theorem (including the interpolation boundary)

Conditional on \(P\), \(Q\) is noncentral chi-square with \(k\) degrees of freedom and finite noncentrality \(\lambda\). Near zero its density is a positive constant times \(q^{k/2-1}\): this follows directly by integrating the positive continuous Gaussian density over a sphere of radius \(\sqrt q\), including the radial Jacobian. Uniform upper bounds hold because \(0\le\lambda\le\lambda_0<\infty\).

Cauchy–Schwarz gives

\[
\frac nQ\le R\le2+\frac{2\lambda+n}{Q}
\le2+\frac{2\lambda_0+n}{Q}.
\]

Consequently, for every real \(r>0\),

\[
\mathbb E[R^r]<\infty\quad\Longleftrightarrow\quad k>2r.
\]

Necessity follows from the lower bound and the positive density at the origin; sufficiency follows from the uniform upper bound, integration of \(q^{k/2-r-1}\) near zero, and bounded inverse powers away from zero.

Thus the mean is infinite at \(p=n\) and \(p=n+1\). At \(p=n+2\) and \(p=n+3\), the mean is finite but the variance is infinite. These are polynomial tails, not light tails in the exponential-moment sense. An ordinary sample mean with a normal standard-error bar is not justified at those latter two widths. Even where the variance exists, its mere existence is not a numerical accuracy certificate for a chosen sample size.

## 3. Exact finite expectation

Assume \(k>2\). Write \(H_d(t)=\mathbb E[1/\chi_d^2(t)]\). The Gaussian Laplace transform gives

\[
H_d(t)=\frac12\int_0^1 u^{d/2-2}e^{-t(1-u)/2}\,du.
\]

Conditional on \(\lambda\), integration of the first and second moments against that same Laplace transform yields

\[
\mathbb E[R\mid\lambda]
=\frac{k-3}{2}
 +\left[n-\frac{(k-5)(k-2)}2-\frac{(k-3)\lambda}2\right]H_k(\lambda).
\]

For completeness, the two identities used here are
\(\mathbb E[\mu^Tx/Q]=1-(k-2)H_k(\lambda)\) and
\(\mathbb E[(\mu^Tx)^2/Q^2]=(k-1)/2-[(k-3)\lambda+(k-1)(k-2)]H_k(\lambda)/2\).
They follow by substituting \(u=(1+2s)^{-1}\) into
\(1/Q=\int_0^\infty e^{-sQ}ds\) and \(1/Q^2=\int_0^\infty s e^{-sQ}ds\); the tilted Gaussian has mean \(u\mu\) and covariance \(uI\).

The beta-product identity reduces the remaining subspace average. If \(U\sim\operatorname{Beta}(1,k/2-1)\) independently of \(T\), then \(V=TU\sim\operatorname{Beta}(1,p/2-1)\), and

\[
\mathbb E[T\mid V]=\frac{k-2}{p-2}+\frac{n-1}{p-2}V.
\]

Both facts follow by splitting a Dirichlet vector with parameters \((1,k/2-1,(n-1)/2)\); the \(n=1\) case follows directly with \(T=1\). It follows that

\[
\mathbb E H_k(\lambda)=\frac{p-2}{k-2}H_p(\lambda_0),
\]

\[
\mathbb E[\lambda H_k(\lambda)]
=\lambda_0H_p(\lambda_0)
 +\frac{n-1}{k-2}\big[(\lambda_0+p-2)H_p(\lambda_0)-1\big].
\]

Substitution gives an exact one-dimensional integral expression:

\[
\mathbb E R=\frac{p-2}{2(p-n-1)}
\left[(p-n-2)+\{2(p-1)-(p-n-2)(p-2+\lambda_0)\}H_p(\lambda_0)\right].
\]

Direct evaluation of this expression can cancel large terms. The implementation instead uses the following nonnegative series.

## 4. Nonnegative Poisson series and deterministic error bound

Let \(J\sim\operatorname{Poisson}(\lambda_0/2)\). Expanding the exponential in the integral (equivalently, the noncentral chi-square mixture) gives
\(H_p(\lambda_0)=\mathbb E[1/(p+2J-2)]\).
Using \(\mathbb E[Jf(J)]=(\lambda_0/2)\mathbb E[f(J+1)]\) in the preceding formula gives

\[
\boxed{\quad\mathbb ER=\frac{p-2}{p-n-1}\mathbb E[g(J)]\quad}
\]

where

\[
g(0)=\frac{p-1}{p-2},\qquad
g(j)=\frac{(p-1)(p-4)+2j(n+1)}{(p+2j-4)(p+2j-2)}\quad(j\ge1).
\]

The separate zero term matters at \(p=4\). For all admissible integers \(p\ge n+2\), including \(p=3,n=1\),

\[
0\le\frac{p-2}{p-n-1}g(j)\le B:=\frac{p-1}{p-n-1}.
\]

This follows by writing, for \(j\ge1\),
\(g(j)=(p-1)/(p+2j-2)-(p-n-2)2j/[(p+2j-4)(p+2j-2)]\), and inspecting the nonnegative numerator in the boxed definition. The upper bound also holds at zero.

This bounded auxiliary Poisson average is an analytic representation of the heavy-tailed risk expectation; **we do not sample it**. For Poisson mean \(m\), the implementation sums relative weights from
\(L=\max(0,\lfloor m-12\sqrt{m+1}-10\rfloor)\) to
\(U=\lceil m+12\sqrt{m+1}+60\rceil\), centered at the mode. Each weight is obtained by multiplying the preceding weight by \(m/j\) or \((j+1)/m\). All sums are nonnegative.

For \(t<m\), the Chernoff bound is
\(P(J\le t)\le\exp[-m+t(1+\log(m/t))]\), with the value \(e^{-m}\) at \(t=0\); for \(t>m\) the same expression bounds \(P(J\ge t)\). These follow by Markov's inequality applied to \(e^{sJ}\), using its elementary generating function \(\exp[m(e^s-1)]\) and choosing \(s=\log(t/m)\).

Apply these at \(L-1\) and \(U+1\) to bound omitted mass by \(\delta\). Normalizing the retained weights gives \(\mathbb E[F(J)\mid L\le J\le U]\), where \(0\le F\le B\). Its difference from \(\mathbb EF\) is therefore at most \(B\delta\). This bound is unrelated to, and does not require, finite variance of the original risk distribution.

The implementation returns an `errorBound` that combines this truncation bound with a conservative floating-point budget. With \(N=U-L+1\), it uses \(\gamma_{8N+64}\), where \(\gamma_q=q\epsilon/(1-q\epsilon)\) and \(\epsilon=2^{-52}\), to cover positive recurrence, summation, rational-term, and normalization arithmetic. This uses twice the usual unit roundoff. It adds \(4\epsilon\max(1,m)\) to cover rounding the mean, using the bound \(|d\mathbb EF/dm|=|\mathbb E[F(J+1)-F(J)]|\le B\). The Chernoff tail allowance is doubled. This is an engineering roundoff budget under normal IEEE-754 arithmetic, not an interval-arithmetic certificate for every JavaScript engine's elementary functions. The rigorous series truncation bound and the exact expectation identity do not depend on that engineering allowance.

The site uses \(1\le n\le20\), \(4\le p\le10000\), and \(0.5\le v\le16\). Unsupported underparameterized inputs return unavailable; the old unconditional-label least-squares formula is not retained as a conditional expectation. A numerical-domain guard likewise returns unavailable instead of silently substituting an asymptotic formula.

## Validation

`verify-sensor-expectation.py` checks the production implementation against high-precision deterministic integration of the independent Laplace expression, including the first finite-mean widths, the exceptional low-dimensional terms, small noise, and 10,000 sensors. Tests also check divergent means and infinite variances at their exact boundaries. The displayed expectation is computed from the exact formula; its numerical evaluation is checked against independent high-precision integration, rather than a Monte Carlo average.


Download the [verifier](verify-sensor-expectation.py) and [production model](../assets/sensor-models.js) into the same folder. The check requires Python 3 with mpmath and Node.js. In a terminal in that folder, run:

```sh
python3 -m venv .venv
.venv/bin/python -m pip install mpmath
.venv/bin/python verify-sensor-expectation.py --model sensor-models.js --output sensor-expectation-check.json
```

The verifier uses 65-digit deterministic Laplace integration, independent of the production Poisson-series evaluation, for 45 parameter settings. It checks each discrepancy against the reported numerical error budget and writes the individual results to the requested JSON file. In a repository checkout, run `python3 blog/data/verify-sensor-expectation.py`; the model path is found automatically and the summary is printed without writing a file.