Theory of the skeleton risk =========================== This chapter develops the asymptotic risk of inductive skeleton fitting from its assumptions to the general scalar formula, then derives the closed forms and comparisons used in :doc:`skeleton_sampling`. The results are organized by mathematical dependence; finite verification and reference tables are collected in :doc:`skeleton_risk_verification`. :cite:t:`tanaka2020asymptotic` derived a Gram recursion and evaluated degrees two and three. The results here give the asymptotic risk and optimal allocation for every dimension and degree, as well as the following explicit evaluations. The conditions in each row are part of the result. .. list-table:: Results and their scope :header-rows: 1 * - Range - Result * - Every :math:`M,D` and active :math:`m` - :prf:ref:`thm-skeleton-scalar-sum`: finite scalar sum, with no Gram inverse. * - :math:`D\geq M` - :prf:ref:`prop-skeleton-top`: :math:`a_{M-1}=\binom{D-1}{M-1}`. * - :math:`m=D-1,D-2,D-3`, where the skeleton exists - :prf:ref:`prop-skeleton-degree-faces`, :prf:ref:`prop-skeleton-second-faces` and :prf:ref:`thm-skeleton-third-faces`: formulas uniform in :math:`M,D`. * - :math:`M=2` or :math:`M=3`, every degree - :prf:ref:`thm-skeleton-curves`, :prf:ref:`cor-skeleton-triangle-edge-risk` and :prf:ref:`thm-skeleton-triangle-vertex-risk`: all coefficients in closed form. * - :math:`m=M-2` and :math:`D\geq M-1` - :prf:ref:`thm-skeleton-codimension-one`: one rational sum. * - Fixed :math:`D`, growing :math:`M` - :prf:ref:`thm-skeleton-rational`: dimension polynomials; :prf:ref:`cor-skeleton-large-dimension`: the limit of every coefficient. * - Risk comparison - Curves win at every degree; :prf:ref:`prop-skeleton-wins` classifies all vertex counts for :math:`D\leq8`; every fixed degree wins for sufficiently large :math:`M`. The proofs are self-contained apart from the cited small-ball inequality. Statements certified by a finite computation are labeled computer-assisted. The final section lists the ranges where further simplification or a complete risk comparison is still open. Setting ------- Let the simplex have :math:`M\geq2` vertices :math:`[M]=\{1,\dots,M\}` and degree :math:`D\geq1`. **Controls and faces.** Controls are indexed by :math:`\alpha \in \mathbb N^M` with :math:`|\alpha| = D`. The Bernstein polynomials are :math:`B_\alpha(t) = \frac{D!}{\alpha!} t^\alpha`, and :math:`\operatorname{supp}\alpha = \{i : \alpha_i > 0\}`. For a nonempty face :math:`I \subseteq [M]`, let .. math:: T_I = \{\alpha : \operatorname{supp}\alpha = I\}, \qquad B_I = \{\alpha : \operatorname{supp}\alpha \subsetneq I\} be the controls interior to :math:`I` and those on its boundary. :math:`T_I` has :math:`\binom{D-1}{|I|-1}` elements, so it is empty when :math:`|I| > D`. On the face :math:`I`, only the polynomials with support in :math:`I` are nonzero. **Gram matrices.** Let :math:`U_I` be the uniform probability measure on the face :math:`I`, let :math:`G^I_{\alpha\beta} = \int B_\alpha B_\beta \, dU_I`, and let :math:`G = G^{[M]}`. The Dirichlet moments .. math:: :label: skeleton-moments \int t^\gamma \, dU_I = \frac{k!\,\gamma!}{(|\gamma| + k)!}, \qquad k = |I| - 1, \quad \operatorname{supp}\gamma \subseteq I, give the entries exactly. Write :math:`(M)_k = M (M+1) \cdots (M+k-1)` for the rising factorial. The block :math:`G^I_{T_I T_I}` is positive definite, because the Bernstein polynomials are linearly independent on the face. **Model and design.** Each response is :math:`x = b(t) + \varepsilon \in \mathbb R^L` with :math:`b(t) = \sum_\alpha B_\alpha(t) b_\alpha`. For every face :math:`I` with :math:`1 \leq |I| \leq \min(D, M)`, :math:`n_I` rows are drawn independently from :math:`U_I`. The noise vectors are independent of the rows and of each other, with mean zero and covariance :math:`\sigma^2 I_L`. Gaussian noise is not needed. **Estimator.** Visit the faces in order of increasing size. On each face :math:`I`, with :math:`T = T_I`, estimate the interior controls by least squares on the rows of :math:`I`, holding the boundary controls at their earlier estimates: .. math:: \hat b_T = \operatorname*{arg\,min}_c \sum_r \Bigl\lVert x_r - \sum_{\beta \in B_I} B_\beta(t_r) \hat b_\beta - \sum_{\alpha \in T} B_\alpha(t_r) c_\alpha \Bigr\rVert^2. This is :func:`~torch_bsf.fitting.fit_skeleton_regression` with unit weights. **Risk and allocation.** The risk is :math:`R_N = E \int \lVert \hat b(t) - b(t) \rVert^2 \, dU_{[M]}(t)` with :math:`N = \sum_I n_I`. Let :math:`N^{(m)}` be the number of rows on faces with :math:`m + 1` vertices. Assume that each such face receives :math:`N^{(m)} / \binom{M}{m+1}` rows, up to one row, and that :math:`N^{(m)} / N \to \rho_m > 0` for every :math:`m < \min(D, M)`. :func:`~torch_bsf.sampling.sample_skeleton` with ``method="random"`` produces exactly this design. **Propagation matrices.** For a face :math:`I` with :math:`T_I \neq \emptyset`, define :math:`E_I` with rows indexed by all controls and columns by :math:`T_I`: - :math:`E_I[T_I]` is the identity; - the row of every control whose support does not contain :math:`I` is zero; - for the faces :math:`J \supsetneq I` with :math:`T_J \neq \emptyset`, in order of increasing size, .. math:: :label: skeleton-propagation E_I[T_J] = -\bigl(G^J_{T_J T_J}\bigr)^{-1} G^J_{T_J B_J} E_I[B_J]. The risk coefficients are .. math:: a_m = \binom{M}{m+1} \sum_{|I| = m+1} \operatorname{tr}\Bigl(\bigl(G^I_{T_I T_I}\bigr)^{-1} E_I^\top G E_I\Bigr), \qquad m = 0, \dots, M - 1. The sum is empty, and :math:`a_m = 0`, when :math:`m \geq D`. :func:`~torch_bsf.sampling.skeleton_risk_coefficients` evaluates this definition in floating point. The asymptotic risk ------------------- .. prf:theorem:: Asymptotic risk :label: thm-skeleton-risk Under the assumptions above, :math:`R_N` is finite for all sufficiently large :math:`N`, and .. math:: \lim_{N \to \infty} N R_N = \sigma^2 L \sum_{m < \min(D, M)} \frac{a_m}{\rho_m}. The proof needs to control the inverse of an empirical Gram matrix. Two lemmas do this. .. prf:lemma:: Small balls of sums :label: lem-skeleton-small-ball Let :math:`X_1, \dots, X_n \geq 0` be independent copies of :math:`X`, and suppose :math:`P(X \leq x) \leq (A x)^a` for all :math:`x > 0` and some :math:`A, a > 0`. Then for every :math:`\delta > 0`, .. math:: P\Bigl(\sum_{i=1}^n X_i \leq n \delta\Bigr) \leq \bigl(e\,\Gamma(a+1) (A \delta)^a\bigr)^n. .. prf:proof:: For :math:`\lambda > 0`, :math:`E e^{-\lambda X} = \int_0^\infty \lambda e^{-\lambda x} P(X \leq x)\,dx \leq \Gamma(a+1) (A / \lambda)^a`. Markov's inequality applied to :math:`e^{-\lambda \sum X_i}` gives :math:`P(\sum X_i \leq n\delta) \leq e^{\lambda n \delta} (E e^{-\lambda X})^n`. Take :math:`\lambda = 1/\delta`. .. prf:lemma:: Inverse moments of empirical Gram matrices :label: lem-skeleton-inverse-moments Let :math:`F` be a face and :math:`\chi = (B_\alpha)_{\alpha \in T}` for a set :math:`T` of :math:`k` controls with supports in :math:`F`, such that :math:`H = \int \chi \chi^\top dU_F` is positive definite. Let :math:`\hat H_n = \frac1n \sum_{r=1}^n \chi(t_r) \chi(t_r)^\top`, where the :math:`t_r` are drawn independently from :math:`U_F`. Then for every :math:`q > 0` there is an :math:`n_0` such that :math:`\sup_{n \geq n_0} E \lVert \hat H_n^{-1} \rVert^q < \infty`. In particular, :math:`\hat H_n` is almost surely invertible for :math:`n \geq n_0`. .. prf:proof:: Let :math:`\psi = H^{-1/2} \chi` and :math:`Q(u) = \frac1n \sum_r (u^\top \psi(t_r))^2`. Let :math:`\lambda` be the minimum of :math:`Q` on the unit sphere. Then :math:`\lVert \hat H_n^{-1} \rVert \leq \lVert H^{-1} \rVert / \lambda`, so it suffices to bound :math:`E \lambda^{-q}`. *Bound for one direction.* For a unit vector :math:`u`, :math:`f = u^\top \psi` is a polynomial of degree at most :math:`D` with :math:`\int f^2 \, dU_F = 1`. The uniform distribution on a convex set is log-concave, so the inequality of :cite:t:`carbery2001distributional` gives :math:`P(|f(t)| \leq s) \leq C D s^{1/D}`, with an absolute constant :math:`C`. Hence :math:`X = f(t)^2` satisfies :prf:ref:`lem-skeleton-small-ball` with :math:`a = 1/(2D)` and :math:`A = (CD)^{2D}`, and .. math:: P\bigl(Q(u) \leq 4 \delta\bigr) \leq (c_1 \delta^a)^n, \qquad c_1 = e\,\Gamma(a+1) (4A)^a. *All directions at once.* Bernstein polynomials take values in :math:`[0, 1]`, so :math:`|\psi| \leq K = \lVert H^{-1/2} \rVert \sqrt k` and :math:`|\sqrt{Q(u)} - \sqrt{Q(v)}| \leq \sqrt{Q(u - v)} \leq K |u - v|`. Take a :math:`(\sqrt\delta / K)`-net of the unit sphere with at most :math:`(3K/\sqrt\delta)^k` points; :math:`\sqrt\delta \leq 1 \leq K` for :math:`\delta \leq 1`, because :math:`\int |\psi|^2 \, dU_F = k`. If :math:`\lambda \leq \delta`, some unit vector :math:`u` has :math:`Q(u) \leq \delta`, and the net point :math:`v` nearest to it has :math:`\sqrt{Q(v)} \leq 2\sqrt\delta`. A union bound over the net gives .. math:: P(\lambda \leq \delta) \leq (3K)^k \delta^{-k/2} (c_1 \delta^a)^n. *Moment bound.* Choose :math:`\delta_1 \leq 1` with :math:`c_1 \delta_1^a \leq 1` and an integer :math:`n_0 > (q + k/2)/a`. For :math:`n \geq n_0` and :math:`\delta \leq \delta_1`, :math:`c_1 \delta^a \leq (\delta / \delta_1)^a \leq 1`, so the right-hand side is at most :math:`(3K)^k \delta_1^{-a n_0} \delta^{a n_0 - k/2}`. Therefore .. math:: E \lambda^{-q} = \int_0^\infty q \delta^{-q-1} P(\lambda \leq \delta) \, d\delta \leq \delta_1^{-q} + q (3K)^k \delta_1^{-a n_0} \int_0^{\delta_1} \delta^{a n_0 - k/2 - q - 1} \, d\delta, which is finite and does not depend on :math:`n`. .. prf:proof:: of :prf:ref:`thm-skeleton-risk` *Exact errors.* On a face :math:`I` with :math:`T = T_I` and :math:`B = B_I`, each row is :math:`x_r = \sum_{\beta \in B} B_\beta(t_r) b_\beta + \sum_{\alpha \in T} B_\alpha(t_r) b_\alpha + \varepsilon_r`. Let :math:`Z_T` and :math:`Z_B` be the design matrices of the rows of :math:`I`, and let :math:`\hat G^I = Z^\top Z / n_I` denote the corresponding blocks of the empirical Gram matrix. When :math:`\hat G^I_{TT}` is invertible, the normal equations give, exactly, the error matrix (one row per control, one column per output) .. math:: e_T = \hat b_T - b_T = \eta_I - \bigl(\hat G^I_{TT}\bigr)^{-1} \hat G^I_{TB} e_B, \qquad \eta_I = (Z_T^\top Z_T)^{-1} Z_T^\top \varepsilon_I. Starting from the vertices, this gives :math:`e = \sum_I \hat E_I \eta_I`, where :math:`\hat E_I` is defined by :eq:`skeleton-propagation` with the empirical blocks :math:`\hat G^J` in place of :math:`G^J`. Condition on all rows. Then the :math:`\eta_I` are independent across faces with mean zero, and each column of :math:`\eta_I` has covariance :math:`\sigma^2 (Z_T^\top Z_T)^{-1} = \sigma^2 (\hat G^I_{TT})^{-1} / n_I`. Since :math:`\int \lVert \hat b - b \rVert^2 dU = \operatorname{tr}(e^\top G e)`, .. math:: N R_N = \sigma^2 L\, E X_N, \qquad X_N = \sum_I \frac{N}{n_I} \operatorname{tr}\Bigl(\bigl(\hat G^I_{TT}\bigr)^{-1} \hat E_I^\top G \hat E_I\Bigr). *Almost sure limit.* :math:`R_N` depends only on the distribution of the design and the noise, so we may take the rows of face :math:`I` to be the first :math:`n_I` terms of one independent sequence per face. Every :math:`n_I` grows without bound. By the strong law of large numbers, :math:`\hat G^J \to G^J` almost surely for every face :math:`J`. The blocks :math:`G^J_{T_J T_J}` are invertible, so :math:`\hat E_I \to E_I`. Since :math:`N / n_I \to \binom{M}{m+1} / \rho_m`, it follows that :math:`X_N \to \sum_m a_m / \rho_m` almost surely. *Uniform integrability.* The entries of :math:`\hat G^J` lie in :math:`[0, 1]`, and :math:`N / n_I` is bounded. Unrolling :eq:`skeleton-propagation` therefore bounds :math:`X_N` by a constant times :math:`\prod_J (1 + \lVert (\hat G^J_{T_J T_J})^{-1} \rVert)^r`, for a fixed :math:`r` and the product over the :math:`F` faces with interior controls. By Hölder's inequality and :prf:ref:`lem-skeleton-inverse-moments` with :math:`q = 2rF`, :math:`\sup_{N \geq N_0} E X_N^2 < \infty`. Hence :math:`R_N` is finite for :math:`N \geq N_0`, the family :math:`X_N` is uniformly integrable, and :math:`E X_N` converges to the almost sure limit. The same argument with a single face, the whole simplex, and all controls gives the all-at-once risk. .. prf:corollary:: All-at-once fitting :label: cor-skeleton-aao-risk For ordinary least squares with all controls on :math:`N` rows drawn independently from :math:`U_{[M]}`, :math:`N R_N \to \sigma^2 L \operatorname{tr}(G^{-1} G) = \sigma^2 L \binom{D+M-1}{D}`, the number of controls. .. prf:theorem:: Optimal ratio :label: thm-skeleton-ratio Over fractions :math:`\rho_m \geq 0` with :math:`\sum_m \rho_m = 1`, :math:`\sum_m a_m / \rho_m` is minimized by :math:`\rho_m = \sqrt{a_m} / \sum_k \sqrt{a_k}`, and the minimum is :math:`(\sum_m \sqrt{a_m})^2`. A zero coefficient receives no rows; the corresponding term is defined as :math:`0/0=0`, while :math:`a_m/0=+\infty` for :math:`a_m>0`. .. prf:proof:: An allocation with :math:`a_m>0` and :math:`\rho_m=0` has infinite risk. For the remaining allocations, apply Cauchy-Schwarz on the indices with positive coefficients: :math:`(\sum_m \sqrt{a_m})^2 = (\sum_m \sqrt{a_m / \rho_m} \sqrt{\rho_m})^2 \leq \sum_m (a_m / \rho_m) \sum_{m:a_m>0} \rho_m \leq \sum_m a_m/\rho_m`. Equality requires zero share on inactive indices and :math:`\sqrt{a_m/\rho_m}` proportional to :math:`\sqrt{\rho_m}` on the active ones, giving the stated ratio. How errors propagate -------------------- Fix a face :math:`I` with :math:`|I| = m + 1 \leq D` and a vector :math:`c \in \mathbb R^{T_I}`. Its *error polynomial* is :math:`p_c(t) = \sum_\alpha B_\alpha(t) (E_I c)_\alpha`. For a face :math:`J`, the *bubble functions* of :math:`J` are the products :math:`t^{\mathbf 1_J} g` with :math:`g` a polynomial of degree at most :math:`D - |J|`. These are exactly the span of :math:`\{B_\alpha : \alpha \in T_J\}` on :math:`J`. .. prf:lemma:: Characterization :label: lem-skeleton-propagation-characterization :math:`p_c` is the unique polynomial of degree at most :math:`D` such that (i) its controls in :math:`T_I` are :math:`c`; (ii) its controls whose support does not contain :math:`I` vanish; (iii) for every face :math:`J \supsetneq I` with :math:`|J| \leq D`, its restriction to :math:`J` is orthogonal in :math:`L^2(U_J)` to the bubble functions of :math:`J`. .. prf:proof:: On :math:`J`, :math:`p_c = \sum_{\operatorname{supp}\beta \subseteq J} B_\beta (E_I c)_\beta`, and the supports contained in :math:`J` are exactly :math:`T_J \cup B_J`. Hence :eq:`skeleton-propagation` reads :math:`\int B_\alpha p_c \, dU_J = 0` for every :math:`\alpha \in T_J`. That is condition (iii). Because :math:`G^J_{T_J T_J}` is invertible, (iii) determines the controls in :math:`T_J` from those in :math:`B_J`. Induction on :math:`|J|` gives uniqueness. Faces with more than :math:`D` vertices have no interior controls and impose no condition. The propagation from :math:`I` to every larger face is determined by the propagation to the faces with one more vertex. .. prf:theorem:: Reduction to one more vertex :label: thm-skeleton-propagation-reduction Fix :math:`j \notin I`. Let :math:`e(\beta, n)` be the control of :math:`p_c` on the face :math:`I \cup \{j\}` with entries :math:`\beta` on :math:`I` and :math:`n` at :math:`j`. By symmetry it does not depend on :math:`j`. Then for every :math:`\alpha` whose support contains :math:`I`, .. math:: :label: skeleton-reduction (E_I c)_\alpha = (n + 1)\, e(\alpha|_I, n) \prod_{l \notin I} \frac{1}{\alpha_l + 1}, \qquad n = \sum_{l \notin I} \alpha_l. Equivalently, let :math:`\phi(x, y) = \sum_{|\beta| + n = D} e(\beta, n) \frac{D!}{\beta!\,n!} x^\beta y^n` be the error polynomial on :math:`I \cup \{j\}` in homogeneous form, and let :math:`\Phi = \partial_y (y \phi)`. Then, with :math:`s_l` independent and uniform on :math:`[0, 1]`, .. math:: :label: skeleton-reduction-integral p_c(t) = E_s\, \Phi\Bigl(t_I, \sum_{l \notin I} s_l t_l\Bigr). .. prf:proof:: *The two forms agree.* :math:`\Phi(x, y) = \sum (n+1) e(\beta, n) \frac{D!}{\beta!\,n!} x^\beta y^n` and :math:`E_s (\sum_l s_l t_l)^n = \sum_{|r| = n} \frac{n!}{r!} t^r \prod_l \frac{1}{r_l + 1}`, because :math:`E s^r = 1/(r+1)`. Multiplying these gives the Bernstein coefficients :eq:`skeleton-reduction`. *Conditions (i) and (ii).* Let :math:`p` be the polynomial defined by :eq:`skeleton-reduction`. If :math:`n = 0`, then :math:`(E_I c)_\alpha = e(\alpha, 0) = c_\alpha`. If :math:`\alpha_i = 0` for some :math:`i \in I`, then :math:`e(\alpha|_I, n)` is a control of :math:`I \cup \{j\}` whose support does not contain :math:`I`, so it vanishes. It remains to check condition (iii) of :prf:ref:`lem-skeleton-propagation-characterization`. *Reducing (iii) to an integral over a cone.* Let :math:`J = I \cup K` with :math:`k = |K| \geq 1` and :math:`|J| \leq D`, and let :math:`g` be a polynomial of degree at most :math:`D - |J|`. Write :math:`x = t_I`, :math:`|x| = \sum_i x_i`, and parametrize :math:`U_J` by :math:`x` and by :math:`t_K` on the slice :math:`\sum t_K = 1 - |x|`. By :eq:`skeleton-reduction-integral` with :math:`w_l = s_l t_l`, .. math:: \prod_{l \in K} t_l \; p = \int_{0 \leq w_l \leq t_l} \Phi\Bigl(x, \sum_l w_l\Bigr) dw. Substitute this into :math:`\int p\, t^{\mathbf 1_J} g \, dU_J`, exchange the integrals, and write :math:`t_K = w + v` with :math:`v \geq 0` and :math:`\sum v = r = 1 - |x| - \sum w`. The integral of :math:`g(x, w + v)` over the simplex of such :math:`v` is :math:`r^{k-1}` times a polynomial in :math:`(x, w, r)` of degree at most :math:`D - |J|`. Integrating over :math:`w` on the slices :math:`\sum w = y` contributes a factor :math:`y^{k-1}`. The integral is therefore a combination of .. math:: \int_{x, y \geq 0,\ |x| + y \leq 1} x^{\mathbf 1_I} \Phi(x, y)\, y^{k-1} r^{k-1} x^a y^b r^c \, dx\, dy, \qquad r = 1 - |x| - y, \quad |a| + b + c \leq D - |J|. *Factoring out the radius.* Substitute :math:`(x, y) = \tau\,\omega`, where :math:`\tau \in (0, 1]` and :math:`\omega` lies in the simplex :math:`S` with :math:`m + 2` vertices, so that :math:`dx\,dy = \tau^{m+1} d\tau\, d\omega`. Since :math:`\Phi` is homogeneous of degree :math:`D`, each integral factors into a :math:`\tau`-integral times .. math:: \int_S \Phi(\omega)\, \omega_I^{\mathbf 1_I}\, \omega_I^a\, \omega_y^{b + k - 1} \, d\omega, and :math:`\omega_I^a \omega_y^{b+k-1}` has degree at most :math:`D - m - 2`. So (iii) holds for every :math:`J` once .. math:: :label: skeleton-reduction-orthogonality \int_S \Phi\, \omega_I^{\mathbf 1_I} h \, d\omega = 0 \quad \text{for every polynomial } h \text{ of degree at most } D - m - 2. *The face with one more vertex.* Run the same computation for :math:`J = I \cup \{j\}`, where :math:`k = 1`. There :math:`p = \phi`, because :math:`\int_0^1 \partial_y(y\phi)(x, sy)\, ds = \phi(x, y)`. The inner integral is :math:`g(x, 1 - |x|)`, an arbitrary polynomial in :math:`x` of degree at most :math:`D - m - 2`. Take :math:`g = x^a`. Then the whole integral is a positive :math:`\tau`-integral times the integral in :eq:`skeleton-reduction-orthogonality` with :math:`h = \omega_I^a`. The orthogonality of :math:`\phi` on :math:`I \cup \{j\}` therefore gives :eq:`skeleton-reduction-orthogonality` for these :math:`h`. These monomials span all polynomials of that degree on :math:`S`, because :math:`\omega_y = 1 - |\omega_I|` there. Scalar evaluation of every coefficient --------------------------------------- The propagation theorem leads to one formula for all dimensions, degrees and skeletons. It separates the angular modes on a face from the moments of the remaining coordinates. Closed forms below are specializations of this formula; a finite sum is available even when no shorter evaluation is known. .. prf:lemma:: The outside-coordinate kernel :label: lem-skeleton-kernel Define the formal power series .. math:: F(x,y)&=\sum_{h,l\geq0}\frac{\binom{h+l}{h}}{(h+1)(l+1)}x^hy^l\\ &=\frac{(1-x-y)\log(1-x-y)-(1-x)\log(1-x) -(1-y)\log(1-y)}{xy},\\ K_q(h,l)&=[x^hy^l]F(x,y)^q,\qquad q,h,l\in\mathbb N. The quotient is interpreted by its formal series, including on the coordinate axes. In particular, :math:`K_0(h,l)=\mathbf1_{h=l=0}` and :math:`K_1(h,l)=\binom{h+l}{h}/((h+1)(l+1))`. For fixed :math:`h,l`, the coefficient is a polynomial in :math:`q` of degree at most :math:`h+l`, given by the finite expansion .. math:: :label: skeleton-kernel-polynomial K_q(h,l)=\sum_{j=0}^{h+l}\binom qj[x^hy^l](F-1)^j. Here :math:`\binom qj=0` for :math:`j>q`. A useful boundary value is :math:`K_2(h,0)=K_2(0,h)=2H_{h+1}/(h+2)`, with :math:`H_j=\sum_{i=1}^j1/i` and :math:`H_0=0`. .. prf:proof:: Integrating :math:`(1-x-y)^{-1}` formally once in each variable gives :math:`xyF` and the displayed logarithmic expression. Expand :math:`F^q=(1+(F-1))^q`; since :math:`F-1` has no constant term, only :math:`j\leq h+l` can contribute to the coefficient. On an axis, :math:`F(x,0)=\sum_{h\geq0}x^h/(h+1)`. Squaring and using :math:`1/((j+1)(h-j+1))=((j+1)^{-1}+(h-j+1)^{-1})/(h+2)` proves the boundary value. We use shifted Jacobi polynomials :math:`P_j^{(\alpha,\beta)}(2v-1)`, orthogonal for :math:`(1-v)^\alpha v^\beta`, with endpoint normalization :math:`P_j^{(\alpha,\beta)}(1)=\binom{j+\alpha}{j}`. .. prf:theorem:: Scalar risk sum :label: thm-skeleton-scalar-sum Let :math:`1\leq s\leq\min(D,M)`, :math:`n=D-s`, and :math:`q=M-s`. Use the kernel of :prf:ref:`lem-skeleton-kernel`. Put .. math:: d_s(k)&=\binom{k+s-2}{s-2}\quad(s\geq2),& d_1(k)&=\mathbf1_{k=0},\\ A_{k,h}&=(-1)^h\binom{n-k}{h}\binom{n+k+3s-1}{h} \quad(0\leq h\leq n-k). Then .. math:: :label: skeleton-scalar-sum a_{s-1}=\frac{\binom Ms^2}{(s-1)!\,(M)_{2D}} \sum_{k=0}^n d_s(k)\sum_{h,l=0}^{n-k} A_{k,h}A_{k,l}(2D+s-h-l-1)!\,h!l!K_q(h,l). All factorial arguments are nonnegative. Each coefficient extraction uses only terms of degree at most :math:`n` in each variable, so every operation is finite, also when :math:`q=0` or :math:`s=1`. .. prf:proof:: *Angular modes.* On the face with :math:`s\geq2` vertices, choose an orthogonal basis :math:`Q_{k,\nu}` of polynomials of exact orthogonal degree :math:`k` for the weight :math:`\prod_i z_i^2` under :math:`U_I`. Here exact orthogonal degree means orthogonal to all polynomials of degree less than :math:`k`; its dimension is :math:`d_s(k)=\binom{k+s-2}{s-2}`. The face bubbles :math:`f_{k,\nu}(z)=\prod_i z_i Q_{k,\nu}(z)` for :math:`0\leq k\leq n` are an orthogonal basis of all face bubbles of degree at most :math:`D`. For :math:`s=1`, use only :math:`Q_{0,1}=1` and :math:`f_{0,1}=1` on the one-point face. *Extension by one vertex.* Write :math:`x=rz`, :math:`r=\sum_i x_i` and let :math:`y` be the new coordinate. The homogeneous polynomial :math:`\Phi` for the boundary mode :math:`f_{k,\nu}` is .. math:: :label: skeleton-mode-extension \Phi_{k,\nu}(x,y)=\prod_i x_i\,r^k Q_{k,\nu}(z) \sum_{h=0}^{n-k} A_{k,h}r^{n-k-h}y^h. Indeed the sum, restricted to :math:`r+y=1`, is :math:`P_{n-k}^{(0,\,3s+2k-1)}(2r-1)`. To verify :eq:`skeleton-reduction-orthogonality`, divide :math:`\Phi` by :math:`\prod_i x_i`. The weighted volume element is :math:`r^{3s-1}\prod_i z_i^2\,dr\,dz`. Expand a test polynomial of degree at most :math:`n-1` into homogeneous angular modes. Angular orthogonality leaves only its degree :math:`k` components. Their remaining radial degrees are at most :math:`n-k-1`, with weight :math:`r^{3s+2k-1}`. Jacobi orthogonality makes these integrals zero. At :math:`y=0` the radial sum is :math:`r^{n-k}`, giving the required boundary mode. This proves the extension by the characterization and reduction theorems. At :math:`n=0` there are no test polynomials. *The remaining coordinates.* By :eq:`skeleton-reduction-integral`, replace each :math:`y^h` by .. math:: T_h(t_{I^c})=E_u\Bigl(\sum_{j\notin I}u_jt_j\Bigr)^h =\sum_{|b|=h}\frac{h!}{b!}\prod_{j\notin I}\frac{t_j^{b_j}}{b_j+1}, where the :math:`u_j` are independent uniform variables on :math:`[0,1]`. If :math:`q\geq1`, write :math:`t_{I^c}=(1-r)v`, with :math:`v` uniform on its :math:`q`-vertex simplex. The Dirichlet moments give .. math:: E_v[T_h(v)T_l(v)]=\frac{h!l!}{(q)_{h+l}}K_q(h,l). The aggregation :math:`t_I=rz` has independent :math:`z\sim U_I` and :math:`r\sim\operatorname{Beta}(s,q)`. Thus the ratio of the squared full-simplex norm of the extended mode to its face norm is .. math:: \frac{1}{(s-1)!\,(M)_{2D}} \sum_{h,l=0}^{n-k} A_{k,h}A_{k,l} (2D+s-h-l-1)!\,h!l!K_q(h,l). Specifically, the radial beta integral has exponents :math:`2D+s-h-l-1` and :math:`q+h+l-1`; its factor :math:`(q+h+l-1)!` cancels the corresponding denominator in the preceding moment. The angular norm cancels as well. Trace is invariant under a change of bubble basis, so summing these norm ratios, with :math:`d_s(k)` modes, gives the trace for one face. Multiply by :math:`\binom Ms^2` for the equal faces and the allocation convention. When :math:`q=0` there is no extension. Only :math:`h=l=0` survives in the stated formula, each norm ratio equals one, and the sum is :math:`\sum_k d_s(k)=\binom{D-1}{s-1}`, including :math:`s=1`. Closed forms near the degree ---------------------------- .. prf:proposition:: Top face :label: prop-skeleton-top For :math:`D \geq M`, :math:`a_{M-1} = \binom{D-1}{M-1}`. .. prf:proof:: For :math:`I = [M]`, :math:`E_I` embeds :math:`T_I` and :math:`G^I = G`, so the trace is :math:`\operatorname{tr}(G_{TT}^{-1} G_{TT}) = |T_I|`. .. prf:proposition:: Faces with D vertices :label: prop-skeleton-degree-faces For :math:`D \leq M`, .. math:: a_{D-1} = \binom{M}{D}^2 \frac{(3D-1)!}{(D-1)!\,(M)_{2D}}. In particular, :math:`a_0 = 2M/(M+1)` for :math:`D = 1`. .. prf:proof:: For :math:`|I| = D`, :math:`T_I = \{\mathbf 1_I\}`. No larger face has interior controls, so :math:`E_I` is a unit vector and the trace is :math:`G_{\mathbf 1_I \mathbf 1_I} / G^I_{\mathbf 1_I \mathbf 1_I}`. By :eq:`skeleton-moments` this equals :math:`\frac{(M-1)!}{(2D+M-1)!} \cdot \frac{(3D-1)!}{(D-1)!}`. Multiply by :math:`\binom{M}{D}` for the faces and by :math:`\binom{M}{D}` from the definition. .. prf:proposition:: Faces with D - 1 vertices :label: prop-skeleton-second-faces For :math:`2 \leq D \leq M + 1`, .. math:: a_{D-2} = \binom{M}{D-1}^2 \frac{(D-1)\,(3D-3)!}{4\,(D-2)!\,(M)_{2D}} \bigl(3M^2 - (6D+3) M + 3D^2 + 15D - 14\bigr). .. prf:proof:: Put :math:`s=D-1`, :math:`n=1` and :math:`q=M-s` in :prf:ref:`thm-skeleton-scalar-sum`. The angular degrees :math:`k=0,1` have multiplicities :math:`1,s-1`; the radial vectors are :math:`(1,-3s)` and :math:`(1)`. With :math:`K_q(0,0)=1`, :math:`K_q(1,0)=q/2` and :math:`K_q(1,1)=q(q+1)/4`, substitution gives .. math:: P_{D,D-2}(M)=\frac{s(3s)!}{4(s-1)!} \bigl(3q^2-9q+12s+4\bigr). Replacing :math:`s=D-1` and :math:`q=M-D+1` gives the stated quadratic. At :math:`q=0` it reduces to :math:`s`, agreeing with :prf:ref:`prop-skeleton-top`, including :math:`s=1`. .. prf:theorem:: Faces with D - 2 vertices :label: thm-skeleton-third-faces Let :math:`D\geq3`, :math:`s=D-2`, and :math:`q=M-s\geq0`. Then .. math:: a_{D-3}=\binom Ms^2\frac{(3s-1)!}{(s-1)!\,(M)_{2D}} \frac{s^2(3s+1)}{64}\,Q_s(q), where the following quartic applies to every such :math:`D,M`: .. math:: Q_s(q)={}&9(3s+1)q^4-78(3s+1)q^3 +(432s^2+1209s+547)q^2\\ &-2(648s^2+1221s+623)q +288(s+1)^2(3s+2). In particular, this includes the full range :math:`D\geq4` requested for :math:`a_{D-3}`, and also :math:`a_0` at :math:`D=3`. .. prf:proof:: Apply :prf:ref:`thm-skeleton-scalar-sum` with :math:`n=2`. The angular degrees :math:`k=0,1,2` have multiplicities :math:`1,s-1,s(s-1)/2`. At :math:`s=1` only the first is present. The radial coefficient vectors :math:`(A_{k,h})_h` are .. math:: (1,-2(3s+1),3s(3s+1)/2),\qquad (1,-(3s+2)),\qquad(1). The needed kernel coefficients, obtained by multiplying the constant, linear and quadratic terms of :math:`F^q`, are .. math:: K_q(0,0)&=1,&K_q(1,0)&=q/2,\\ K_q(2,0)&=q(3q+5)/24,&K_q(1,1)&=q(q+1)/4,\\ K_q(2,1)&=q(q+2)(3q+5)/48,\\ K_q(2,2)&=q(9q^3+66q^2+163q+146)/576. Use symmetry for the reversed arguments. Substitution into the scalar sum gives the displayed quartic after factoring out :math:`(3s-1)!s^2(3s+1)/(64(s-1)!)`. This is a symbolic evaluation with arbitrary :math:`s,q`, rather than an interpolation in the degree. At :math:`q=0` the expression is :math:`s(s+1)/2=\binom{D-1}{M-1}`, as required by :prf:ref:`prop-skeleton-top`. Faces of codimension one ------------------------ .. prf:theorem:: Faces of codimension one :label: thm-skeleton-codimension-one For :math:`s\geq1`, :math:`M=s+1`, and :math:`D\geq s`, put :math:`n=D-s` and use :math:`d_s(k)` from :prf:ref:`thm-skeleton-scalar-sum`. Then .. math:: a_{s-1}(s+1,D)=s(s+1)^2\sum_{k=0}^n \frac{d_s(k)}{(n-k+1)(n+k+3s)}. Thus the two radial sums and the kernel extraction reduce to a single sum of :math:`n+1` rational terms; there are no polynomial integrals, Gram matrices or propagation coefficients to compute. .. prf:proof:: In the mode extension of :prf:ref:`thm-skeleton-scalar-sum`, there is just one outside coordinate :math:`y`. Averaging :math:`y^h` multiplies its coefficient by :math:`1/(h+1)`. Write :math:`L=n-k` and :math:`\beta=3s+2k-1`. The homogeneous radial polynomial therefore becomes, on :math:`r+y=1`, .. math:: \sum_{h=0}^L\frac{(-1)^h}{h+1}\binom Lh\binom{L+\beta}{h} r^{L-h}(1-r)^h =\frac{P_L^{(1,\beta)}(2r-1)}{L+1}. The beta density of :math:`r` is :math:`sr^{s-1}`. After the angular norm cancels, the full-to-face squared-norm ratio of this mode is .. math:: \frac{s}{(L+1)^2}\int_0^1 r^\beta P_L^{(1,\beta)}(2r-1)^2\,dr =\frac{s}{(L+1)(L+\beta+1)}. For completeness, the last integral follows from the reproducing kernel for polynomials of degree at most :math:`L` with weight :math:`r^\beta`. The shifted polynomials :math:`P_j^{(0,\beta)}` have squared norms :math:`1/(2j+\beta+1)` and endpoint value one. Hence its kernel at endpoint one is :math:`C_L(r)=\sum_{j=0}^L(2j+\beta+1)P_j^{(0,\beta)}(2r-1)`, with :math:`\int r^\beta C_L(r)g(r)\,dr=g(1)`. Taking :math:`g=(1-r)h`, :math:`\deg h