Optimal skeleton sampling

When responses can be measured at chosen parameters, a simplex can be fitted face by face: vertices first, then edges, then larger faces, each time estimating only the controls interior to the face from observations on that face. This is the inductive skeleton fitting of Tanaka et al. [TSKH20]. Its accuracy depends on how a budget of \(N\) observations is split among the skeletons. This chapter gives the asymptotic risk for every dimension and degree, the split that minimizes it, and the functions that compute, sample and fit with it.

Sample, fit and compare

import torch
from torch_bsf import BezierModel, BezierSpace, Samples, fit_skeleton_regression
from torch_bsf.sampling import optimal_skeleton_ratio, sample_skeleton, skeleton_risk_coefficients

space = BezierSpace.simplex(2, degree=3)
print([round(a, 4) for a in skeleton_risk_coefficients(space)])
print([round(r, 3) for r in optimal_skeleton_ratio(space)])

truth = BezierModel(space, torch.linspace(-1, 1, space.n_controls, dtype=torch.float64)[:, None])
points = sample_skeleton(space, 300, seed=0)
noise = 0.01 * torch.randn(300, 1, generator=torch.Generator().manual_seed(1), dtype=torch.float64)
result = fit_skeleton_regression(space, Samples(points, truth(points).detach() + noise)).require_success()
assert sorted(result.faces)[:3] == [(0,), (0, 1), (0, 1, 2)]
assert (result.model.coefficients - truth.coefficients).abs().max() < 0.1

[0.1464, 3.5357, 1.0]
[0.117, 0.576, 0.306]

sample_skeleton() places 11.7% of the rows on vertices, 57.6% on edges and 30.6% inside the triangle by default, splits each skeleton’s share equally among its faces, and draws uniformly on each face. Pass ratio to choose other shares. fit_skeleton_regression() assigns each row to the face spanned by its nonzero coordinates, so zeros must be exact; it reports one FitResult per face. Rows on faces with more than \(D\) vertices carry no information for this estimator and are ignored, and rows of zero weight do not count towards the row a face needs.

Use fit() instead when observations are not on faces or when all of them should inform every control. Fitting stratified rows with fit is a different estimator, and the ratio below is not optimal for it. fit_skeleton() is the value-only variant for unknown parameters, discussed in Fitting from values alone.

The risk and its coefficients

Let the simplex have \(M\) vertices and degree \(D\), let responses \(x = b(t) + \varepsilon\) have \(L\) outputs with independent noise of variance \(\sigma^2\), and let \(N^{(m)}\) rows lie on the \(m\)-skeleton (the union of faces with \(m + 1\) vertices), split equally among its \(\binom{M}{m+1}\) faces and uniform on each. Measure the error by

\[R = E \int_{\Delta^{M-1}} \lVert \hat b(t) - b(t) \rVert^2 \, dU(t),\]

with \(U\) uniform on the simplex. Then

\[R = \sigma^2 L \sum_{m=0}^{M-1} \frac{a_m}{N^{(m)}} + o\!\left(\frac{1}{N}\right).\]

skeleton_risk_coefficients() returns \(a_0, \dots, a_{M-1}\). They follow from three facts, proved in Theory of the skeleton risk.

  1. On a face \(I\), let \(T\) be the controls whose multi-index has support exactly \(I\) and \(B\) those with smaller support. Only these controls are nonzero on the face. With \(n_I\) rows, least squares for \(T\) with \(B\) fixed gives, to first order, the error

    \[e_T = \eta_I - G^I_{TT}{}^{-1} G^I_{TB}\, e_B,\]

    where \(G^I\) is the Gram matrix of the Bernstein basis under the uniform distribution on the face. The innovations \(\eta_I\) are independent across faces, with covariance \(\sigma^2 G^I_{TT}{}^{-1} / n_I\) to first order. Fluctuations of the empirical Gram matrices contribute only \(o(1/N)\).

  2. Starting from the vertices, the errors of all controls are therefore a linear map \(e = \sum_I E_I \eta_I\) of the independent face noises.

  3. The risk is \(E\,\mathrm{tr}(e^\top G e)\) with \(G\) the Gram matrix on the whole simplex. With \(n_I = N^{(m)} / \binom{M}{m+1}\),

    \[a_m = \binom{M}{m+1} \sum_{|I| = m+1} \mathrm{tr}\!\left(G^I_{TT}{}^{-1} E_I^\top G E_I\right).\]

The Gram entries are exact moments of the uniform distribution on a simplex, \(E[t^\alpha] = k!\,\alpha! / (|\alpha| + k)!\) on a \(k\)-simplex.

The available formulas are summarized below. Their proofs and exact ranges are in Theory of the skeleton risk; independent checks and reference polynomials are in Exact verification of the skeleton risk.

  • \(a_m = 0\) for \(m \geq D\): faces with more than \(D\) vertices have no interior controls, and the optimal ratio gives them no rows.

  • For \(D \geq M\), \(a_{M-1} = \binom{D-1}{M-1}\), the number of controls interior to the whole simplex. Their estimate is measured with the same Gram matrix it was fitted with.

  • For \(D \leq M\), \(a_{D-1} = \binom{M}{D}^2 (3D-1)! / ((D-1)!\,(M)_{2D})\), where \((M)_k = M (M+1) \cdots (M+k-1)\). For \(D = 1\) this is \(a_0 = 2M / (M + 1)\). \(a_{D-2}\) and \(a_{D-3}\) also have closed forms in both \(D\) and \(M\); see Proposition 3 and Theorem 11.

  • On curves (\(M = 2\)), \(a_0 = 4 / (D (D + 2))\) and \(a_1 = D - 1\).

  • On triangles (\(M=3\)), all coefficients are closed for every degree: Theorem 15 gives \(a_0(3,D)\) in harmonic numbers and \(\binom{2D}{D}\) for \(D\geq1\); \(a_1=18(H_{D-1}+H_{2D+2}-H_{D+3})/(2D+3)\) for \(D\geq2\); and \(a_2=\binom{D-1}{2}\) for \(D\geq3\). Here \(H_j=\sum_{i=1}^j1/i\); inactive skeletons have coefficient zero. Hence \((\sqrt{a_0}+\sqrt{a_1}+\sqrt{a_2})^2\) is an explicit formula for the optimal normalized triangle risk at every degree.

  • Every codimension-one coefficient has a single finite rational sum: see Theorem 12.

  • For a fixed degree, every \(a_m\) is a rational function of \(M\) and converges as \(M\) grows. Corollary 5 gives every limit; for example, \(a_0 \to 2 D^2 (D+1)^2 / 4^D\).

These closed forms rest on a structure theorem: the errors that a face passes to larger faces are determined by the errors it passes to the faces with one more vertex. More generally, Theorem 10 gives every coefficient as a finite scalar sum, without Gram-matrix inversion or propagation recursion. Simplification to a bounded number of functions remains open in the ranges listed in the theory chapter.

The optimal ratio

Minimizing \(\sum_m a_m / N^{(m)}\) subject to \(\sum_m N^{(m)} = N\) gives, by the Cauchy-Schwarz inequality,

\[\frac{N^{(m)}}{N} = \frac{\sqrt{a_m}}{\sum_k \sqrt{a_k}}, \qquad R^\ast = \frac{\sigma^2 L}{N} \Bigl(\sum_m \sqrt{a_m}\Bigr)^2.\]

optimal_skeleton_ratio() returns these fractions. For comparison, ordinary least squares on \(N\) uniform rows on the whole simplex (all-at-once fitting) has risk \(\sigma^2 L \binom{D+M-1}{D} / N\), the number of controls over \(N\). The table lists \(N R^\ast / (\sigma^2 L)\) and the all-at-once value; the smaller one is in bold.

Table 10 Optimal inductive skeleton risk versus all-at-once risk

Degree

\(M = 2\)

\(M = 3\)

\(M = 4\)

\(M = 5\)

\(M = 6\)

1

1.333 / 2

1.5 / 3

1.6 / 4

1.667 / 5

1.714 / 6

2

2.914 / 3

5.496 / 6

8.667 / 10

11.99 / 15

15.17 / 21

3

3.727 / 4

10.65 / 10

23.88 / 20

44.75 / 35

73.14 / 56

4

4.581 / 5

15.07 / 15

43.04 / 35

104.9 / 70

220.2 / 126

5

5.467 / 6

20.20 / 21

64.12 / 56

184.3 / 126

471.8 / 252

6

6.374 / 7

26.17 / 28

90.36 / 84

283.9 / 210

824.0 / 462

The degree-two and degree-three rows agree with Tables 1 and 2 of Tanaka et al. [TSKH20], which give these coefficients numerically for \(D \in \{2, 3\}\) and \(2 \leq M \leq 8\). The other rows are new. They show that the paper’s observation that the inductive skeleton fitting wins for \(D < 3\) does not extend to a rule in the degree: it wins for every degree shown on curves (\(M = 2\)), and on triangles at degrees five and six. Theory of the skeleton risk proves that it wins:

  • on curves at every degree;

  • for \(D \leq 2\) in every dimension;

  • for every degree once the dimension is large enough, for example from \(M = 11\) at degree three, from \(M = 20\) at degree four and from \(M = 29\) at degree five, \(M = 39\) at degree six, \(M = 49\) at degree seven and \(M = 60\) at degree eight.

Notes on the original derivation

The numerical results of the paper are correct; the statements around them need three qualifications.

  • The risk is an expectation, so the limit concerns \(N R_N\); writing convergence in probability of \(E[(Z^\top Z)^{-1}]\) mixes two notions.

  • If rows are assigned to faces at random, some face receives no rows with positive probability for every finite \(N\), and the estimator and \(R_N\) are undefined. The analysis needs a deterministic split, which sample_skeleton() uses (up to one row per face).

  • The matrices \(\Lambda^{(m)[k]}\) of the paper’s recursion are indexed inconsistently. The derivation above avoids them by propagating the control errors face by face.

Limits of the result

The coefficients describe the limit as every active face receives more rows; they are not finite-sample error bounds. With few rows per face, a design can be rank deficient or have large inverse-Gram moments. The formulas assume the uniform distribution on each face, an equal split among faces of the same size, and risk averaged uniformly over the simplex; a different weighting of the error changes them. Bernstein Gram matrices become ill-conditioned as the degree grows, and skeleton_risk_coefficients() raises ValueError once the float64 result can no longer be guaranteed to six significant digits (above degree 14 for \(M = 2\), 12 for \(M = 3\), 11 for \(M = 4\), 9 for \(M = 5\) and 8 for \(M = 6\)). In that case, pass ratio to sample_skeleton() explicitly. Only a single simplex is supported, not products or assembled spaces.