Streaming Polyphase Factor Decomposition / SKAT

Download original file

This is the original TeX manuscript source. Use the paper’s PDF for typeset equations.

Read the typeset paper

1 \documentclass[11pt]{article}
2 
3 \usepackage[letterpaper,margin=0.82in]{geometry}
4 \usepackage[T1]{fontenc}
5 \usepackage[utf8]{inputenc}
6 \usepackage{lmodern}
7 \usepackage{microtype}
8 \usepackage{amsmath,amssymb,amsthm,mathtools}
9 \usepackage{booktabs,longtable,array}
10 \usepackage{enumitem}
11 \usepackage{xcolor}
12 \usepackage{hyperref}
13 \usepackage{fancyhdr}
14 \usepackage{listings}
15 \usepackage{graphicx}
16 \usepackage{float}
17 \usepackage{caption}
18 \captionsetup{font=small,labelfont=bf}
19 \graphicspath{{figures_final/}}
20 
21 \hypersetup{
22   colorlinks=true,
23   linkcolor=blue!45!black,
24   citecolor=blue!45!black,
25   urlcolor=blue!45!black,
26   pdftitle={Streaming Polyphase Factor Decomposition for SKAT-like Tests},
27   pdfauthor={Komorebi P.}
28 }
29 
30 \setlength{\parindent}{0pt}
31 \setlength{\parskip}{6pt}
32 \setlength{\headheight}{14pt}
33 \setlength{\emergencystretch}{2em}
34 
35 \newtheorem{theorem}{Theorem}[section]
36 \newtheorem{proposition}[theorem]{Proposition}
37 \newtheorem{corollary}[theorem]{Corollary}
38 \newtheorem{lemma}[theorem]{Lemma}
39 \theoremstyle{definition}
40 \newtheorem{definition}[theorem]{Definition}
41 \newtheorem{algorithm}[theorem]{Algorithm}
42 \theoremstyle{remark}
43 \newtheorem{remark}[theorem]{Remark}
44 
45 \pagestyle{fancy}
46 \fancyhf{}
47 \fancyhead[L]{Streaming Polyphase Factor Decomposition}
48 \fancyhead[R]{Komorebi P.}
49 \fancyfoot[C]{\thepage}
50 
51 \title{\textbf{Streaming Polyphase Factor Decomposition for SKAT-like Quadratic-Form Tests}\\
52 \large Memory-Bounded Factor Actions, Phase-Local Krylov Structure, and Spectral Tail Compression}
53 \author{Komorebi P.}
54 \date{Illustrated manuscript v0.8 -- 2026-09-23}
55 
56 \begin{document}
57 \maketitle
58 
59 \begin{abstract}
60 Sequence kernel association testing (SKAT) is a variance-component framework for set-based rare-variant association testing \cite{Wu2011}. Large rare-variant and kernel tests often begin with a weighted factor
61 \[
62 A=M_XGW^{1/2},
63 \qquad
64 K=AA^\top,
65 \]
66 where \(G\in\mathbb R^{n\times m}\) is a genotype or feature matrix,
67 \(W\) is diagonal and nonnegative, and \(M_X\) projects away null-model
68 covariates.  When \(m\) is large, the direct factor \(A\) may exceed device
69 memory even though the final statistic requires only a spectrum or repeated
70 operator actions.  This manuscript formalizes a composable decomposition
71 calculus with four independent layers: storage blocks, structural
72 polyphase packets, Krylov/range compression, and final spectral packets for
73 deep-tail evaluation.  The core exact result is that both \(K\) and every
74 Krylov action can be evaluated by streaming variant blocks, so total variant
75 count controls arithmetic and I/O but need not control device working-set
76 size.  Storage refinement is proved invariant, two-dimensional sample/variant
77 tiling is exact, and the streaming layer composes with an exact polyphase
78 light-cone theorem, minimal phase-local rank allocation, low-rank interface
79 messages, and exact QR-core updates.  When phase couplings are weak rather
80 than zero, a nonnegative transfer matrix gives certified phasewise leakage
81 bounds for operator powers and complete Krylov polynomials.  A dual
82 query/spectral budget propagates pruning uncertainty to generalized-\(\chi^2\)
83 tail brackets, while residual trace, square-trace, and spectral-cap invariants
84 supply substantially sharper incomplete-range certificates and an exact
85 moment-triggered residual-completion rule.
86 Reference implementations reproduce
87 explicitly formed operators to roughly machine precision.  For
88 \(n=2000\), Krylov width \(16\), covariate rank \(4\), and resident variant
89 block \(1024\), the algebraic streamed kernel-action working set is about
90 \(16.3\) MiB and remains constant from \(m=10^4\) through \(m=10^6\), while
91 the full float64 factor grows from about \(153\) MiB to \(15.3\) GiB.
92 \end{abstract}
93 
94 \tableofcontents
95 \clearpage
96 
97 \section{Problem and decomposition axes}
98 
99 Let
100 \[
101 G\in\mathbb R^{n\times m},
102 \qquad
103 W=\operatorname{diag}(w_1,\ldots,w_m),
104 \quad w_j\ge 0,
105 \]
106 and let \(Q_X\in\mathbb R^{n\times p}\) have orthonormal columns spanning
107 the null-model covariates.  Define
108 \[
109 M_X=I-Q_XQ_X^\top,
110 \qquad
111 A=M_XGW^{1/2},
112 \qquad
113 K=AA^\top.
114 \]
115 
116 The central implementation difficulty is not the algebraic size of the
117 model but the mismatch between the total factor size \(n\times m\) and the
118 working memory of an accelerator.  The key point of this manuscript is that
119 several decompositions act on different mathematical coordinates and therefore
120 compose:
121 \[
122 \boxed{
123 \text{storage blocks}
124 \circ
125 \text{structural phases}
126 \circ
127 \text{Krylov/range compression}
128 \circ
129 \text{spectral tail packets}.
130 }
131 \]
132 
133 Storage blocking controls residence.  Structural phases encode LD windows,
134 annotations, frequency bands, or other coordinate groups.  Krylov methods
135 compress the operator to the rank actually requested.  Spectral packets
136 compress the final generalized-\(\chi^2\) transform.
137 
138 
139 \begin{figure}[H]
140 \centering
141 \includegraphics[width=.72\linewidth]{figure_01_residualization_geometry.png}
142 \caption{Residualization is the geometric starting point. Every sample-space vector splits into a component in the covariate span and an orthogonal residual component. The factor $A=M_XGW^{1/2}$ therefore acts only on variation remaining after the null-model covariates have been removed.}
143 \label{fig:residualization-geometry}
144 \end{figure}
145 
146 \section{Variant-block Gram decomposition}
147 
148 Partition the variant indices into disjoint blocks
149 \[
150 \{1,\ldots,m\}
151 =
152 J_1\sqcup\cdots\sqcup J_B.
153 \]
154 Let
155 \[
156 G_b=G_{[:,J_b]},
157 \qquad
158 W_b=W_{J_b,J_b},
159 \qquad
160 A_b=M_XG_bW_b^{1/2}.
161 \]
162 Then
163 \[
164 A=[A_1\ \cdots\ A_B].
165 \]
166 
167 \begin{theorem}[Additive block Gram identity]
168 \[
169 \boxed{
170 K
171 =
172 AA^\top
173 =
174 \sum_{b=1}^B A_bA_b^\top
175 =
176 M_X
177 \left(
178 \sum_{b=1}^B G_bW_bG_b^\top
179 \right)
180 M_X.
181 }
182 \]
183 \end{theorem}
184 
185 \begin{proof}
186 Horizontal concatenation gives
187 \[
188 [A_1\ \cdots\ A_B]
189 \begin{bmatrix}
190 A_1^\top\\
191 \vdots\\
192 A_B^\top
193 \end{bmatrix}
194 =
195 \sum_bA_bA_b^\top.
196 \]
197 The final form follows from
198 \(A_b=M_XG_bW_b^{1/2}\) and
199 \(M_X^\top=M_X=M_X^2\).
200 \end{proof}
201 
202 \begin{remark}
203 The theorem is exact for every partition.  Block boundaries are execution
204 parameters, not statistical parameters.
205 \end{remark}
206 
207 
208 \begin{figure}[H]
209 \centering
210 \includegraphics[width=.92\linewidth]{figure_02_exact_block_kernel_assembly.png}
211 \caption{Exact block-kernel assembly. The sample interaction kernel carried by each resident variant block is $K_b=A_bA_b^\top$, and the full kernel is their exact sum. Streaming changes the order of evaluation, not the operator.}
212 \label{fig:block-kernel-assembly}
213 \end{figure}
214 
215 \section{Exact streamed kernel actions}
216 
217 \begin{theorem}[One-pass sample-space streaming action]
218 For every block vector \(V\in\mathbb R^{n\times s}\),
219 \[
220 \boxed{
221 KV
222 =
223 M_X
224 \sum_{b=1}^B
225 G_bW_bG_b^\top
226 M_XV.
227 }
228 \]
229 \end{theorem}
230 
231 \begin{proof}
232 Apply the additive block Gram identity to \(KV\), then associate each block
233 product as
234 \[
235 G_bW_b(G_b^\top M_XV).
236 \]
237 \end{proof}
238 
239 \begin{algorithm}[Streamed \(KV\)]
240 \begin{enumerate}[label=\arabic*.]
241 \item Compute \(U=M_XV\).
242 \item Set \(Y=0\).
243 \item For each resident block \(b\), compute
244 \[
245 T_b=G_b^\top U,
246 \qquad
247 Y\leftarrow Y+G_b(W_bT_b).
248 \]
249 \item Return \(M_XY\).
250 \end{enumerate}
251 \end{algorithm}
252 
253 \begin{corollary}[Working memory independent of total variant count]
254 If at most \(c\) variants are resident and the block-vector width is \(s\),
255 then the algebraic working set is
256 \[
257 \boxed{
258 O(nc+ns+cs+np),
259 }
260 \]
261 independent of total variant count \(m\).
262 \end{corollary}
263 
264 \section{Streamed factor actions and Golub--Kahan}
265 
266 Let \(V\in\mathbb R^{m\times s}\) be partitioned conformally:
267 \[
268 V=
269 \begin{bmatrix}
270 V_1\\
271 \vdots\\
272 V_B
273 \end{bmatrix}.
274 \]
275 
276 \begin{theorem}[Streamed \(AV\)]
277 \[
278 \boxed{
279 AV
280 =
281 M_X
282 \sum_{b=1}^B
283 G_bW_b^{1/2}V_b.
284 }
285 \]
286 \end{theorem}
287 
288 \begin{theorem}[Streamed \(A^\top U\)]
289 For \(U\in\mathbb R^{n\times s}\),
290 \[
291 \boxed{
292 (A^\top U)_b
293 =
294 W_b^{1/2}G_b^\top M_XU.
295 }
296 \]
297 \end{theorem}
298 
299 \begin{proof}
300 Both identities follow by substituting
301 \(A=M_XGW^{1/2}\) and partitioning the variant coordinate.
302 \end{proof}
303 
304 \begin{corollary}[Streaming Golub--Kahan equivalence]
305 Any exact-arithmetic algorithm whose access to \(A\) occurs only through
306 \(AV\) and \(A^\top U\) produces the same iterates when those products are
307 implemented by the preceding streamed formulas.
308 \end{corollary}
309 
310 This applies directly to block Golub--Kahan bidiagonalization, matrix-free
311 singular-value iteration, range finding, and restarted variants.  The
312 factor need not be resident.
313 
314 \section{Two-pass dual action}
315 
316 The variant-space normal operator is
317 \[
318 H=A^\top A.
319 \]
320 
321 \begin{theorem}[Two-pass streamed dual action]
322 For every \(V\in\mathbb R^{m\times s}\),
323 \[
324 \boxed{
325 HV=A^\top(AV).
326 }
327 \]
328 The product is evaluated exactly by a streamed pass forming
329 \(Y=AV\in\mathbb R^{n\times s}\), followed by a second streamed pass
330 emitting the blocks
331 \[
332 (HV)_b=W_b^{1/2}G_b^\top M_XY.
333 \]
334 \end{theorem}
335 
336 This form is particularly natural for variant-space polyphase Krylov
337 iterations because the right-coordinate blocks can be produced and consumed
338 sequentially.
339 
340 \section{Two-dimensional sample/variant tiling}
341 
342 A variant block may itself exceed device memory.  Let
343 \[
344 H_b=G_bW_b^{1/2}
345 \]
346 and partition samples into row tiles \(I_1,\ldots,I_R\).
347 
348 \begin{theorem}[Two-pass residualized tile construction]
349 Define
350 \[
351 R_b=Q_X^\top H_b.
352 \]
353 Then
354 \[
355 \boxed{
356 R_b
357 =
358 \sum_r
359 Q_{X,I_r}^\top H_{b,I_r},
360 }
361 \]
362 and after \(R_b\) is known,
363 \[
364 \boxed{
365 A_{b,I_r}
366 =
367 H_{b,I_r}-Q_{X,I_r}R_b.
368 }
369 \]
370 Consequently,
371 \[
372 \boxed{
373 K_{I_r,I_s}
374 =
375 \sum_b
376 A_{b,I_r}A_{b,I_s}^\top.
377 }
378 \]
379 \end{theorem}
380 
381 \begin{proof}
382 The first identity is the row-partitioned matrix product
383 \(Q_X^\top H_b\).  The second is the definition
384 \(M_XH_b=H_b-Q_XQ_X^\top H_b\).  The last identity is the corresponding
385 sample-tile expansion of \(AA^\top\).
386 \end{proof}
387 
388 With sample tile height \(t\) and resident variant width \(c\), the dominant
389 working storage is
390 \[
391 O(tc+pc+t^2).
392 \]
393 
394 \section{Partition refinement and permutation invariance}
395 
396 \begin{theorem}[Partition-refinement invariance]
397 Let \(\mathcal J\) be any partition of the variants and let
398 \(\mathcal J'\) refine \(\mathcal J\).  Then
399 \[
400 \boxed{
401 \sum_{J\in\mathcal J}A_JA_J^\top
402 =
403 \sum_{J'\in\mathcal J'}A_{J'}A_{J'}^\top
404 =
405 AA^\top.
406 }
407 \]
408 The analogous streamed operator action is also invariant.
409 \end{theorem}
410 
411 \begin{proof}
412 Each coarse block is a disjoint union of refined blocks.  Apply the additive
413 Gram theorem within each coarse block and sum.
414 \end{proof}
415 
416 \begin{corollary}[Adaptive memory splitting]
417 A runtime may recursively subdivide any block that exceeds available device
418 memory without changing the exact operator.
419 \end{corollary}
420 
421 \begin{proposition}[Variant permutation invariance]
422 If the columns of \(G\) and the diagonal entries of \(W\) are permuted by the
423 same permutation, \(K\) and every streamed sample-space action are unchanged.
424 \end{proposition}
425 
426 
427 \begin{figure}[H]
428 \centering
429 \includegraphics[width=.84\linewidth]{figure_08_partition_invariance.png}
430 \caption{Partition-refinement invariance in floating-point tests. Coarse, fine, irregular, singleton-heavy, and permuted blockings reproduce the explicitly formed Gram matrix and its action to machine precision. The block schedule is therefore an execution choice rather than a change in the statistical model.}
431 \label{fig:partition-invariance}
432 \end{figure}
433 
434 \section{Polyphase structure}
435 
436 Let \(P_h\) project onto variant phase \(h\), and again set
437 \[
438 H=A^\top A.
439 \]
440 Define the exact phase graph by
441 \[
442 g\to h
443 \quad\Longleftrightarrow\quad
444 P_hHP_g\ne0.
445 \]
446 
447 \begin{theorem}[Polyphase light cone]
448 If the graph distance from phase \(g\) to phase \(h\) exceeds \(k\), then
449 \[
450 \boxed{
451 P_hH^kP_g=0.
452 }
453 \]
454 \end{theorem}
455 
456 \begin{proof}
457 Insert the projector resolution
458 \(\sum_jP_j=I\) between successive factors of \(H\).  Every nonzero term in
459 \(P_hH^kP_g\) corresponds to a phase path of length \(k\).  If no such path
460 exists, every term vanishes.
461 \end{proof}
462 
463 
464 \begin{figure}[H]
465 \centering
466 \includegraphics[width=.82\linewidth]{figure_04_polyphase_lightcone.png}
467 \caption{The exact polyphase light cone. The support of the Krylov iterate expands by at most one graph edge per operator application. Phases outside graph distance $k$ remain identically inactive at depth $k$ and need not be loaded.}
468 \label{fig:polyphase-lightcone}
469 \end{figure}
470 
471 \begin{theorem}[Streaming-polyphase composition]
472 At Krylov depth \(k\), an exact matrix-free iteration may omit every phase
473 outside the exact light cone.  Within each active phase, arbitrary storage
474 refinement may be used.
475 \end{theorem}
476 
477 \begin{proof}
478 The light-cone theorem removes inactive phase contributions exactly.
479 Partition-refinement invariance shows that the remaining active phase
480 operators may be streamed through any exact subpartition.
481 \end{proof}
482 
483 This is the key composition law:
484 \[
485 \boxed{
486 \text{mathematical locality}
487 \longrightarrow
488 \text{active phases}
489 \longrightarrow
490 \text{active storage blocks}.
491 }
492 \]
493 
494 \section{Certified approximate light cones}
495 \label{sec:approx-lightcone}
496 
497 The exact light-cone theorem requires literal zero phase couplings.  Real LD,
498 annotation, and feature graphs are often only approximately sparse.  This
499 section replaces exact support by a computable phasewise envelope.
500 
501 Let the phase projectors be $P_1,\ldots,P_q$, and let $X$ be a vector or block
502 vector.  Define its phase-norm profile
503 \[
504 \rho(X)_g=\lVert P_gX\rVert_F,
505 \]
506 and define the nonnegative phase-transfer matrix
507 \[
508 \mathcal B_{hg}=\lVert P_hHP_g\rVert_2.
509 \]
510 The matrix $\mathcal B$ records the strength of every one-step phase transfer
511 without assuming that weak couplings vanish.
512 
513 \begin{theorem}[Weighted phase-transfer envelope]
514 \label{thm:weighted-transfer}
515 For every block vector $X$,
516 \[
517 \boxed{
518 \rho(HX)\leq \mathcal B\rho(X)
519 }
520 \]
521 componentwise.  Consequently, for every integer $k\geq0$,
522 \[
523 \boxed{
524 \rho(H^kX)\leq \mathcal B^k\rho(X).
525 }
526 \]
527 \end{theorem}
528 
529 \begin{proof}
530 Using $\sum_gP_g=I$ and submultiplicativity,
531 \[
532 \begin{aligned}
533 \lVert P_hHX\rVert_F
534 &=\left\lVert\sum_gP_hHP_gX\right\rVert_F\\
535 &\leq\sum_g\lVert P_hHP_g\rVert_2\lVert P_gX\rVert_F\\
536 &=(\mathcal B\rho(X))_h.
537 \end{aligned}
538 \]
539 Iteration gives the power bound.
540 \end{proof}
541 
542 Now choose a symmetric set of retained phase blocks and define
543 \[
544 H_0=\sum_{(h,g)\in\mathcal E_0}P_hHP_g.
545 \]
546 Let
547 \[
548 (\mathcal B_0)_{hg}=
549 \begin{cases}
550 \mathcal B_{hg},&(h,g)\in\mathcal E_0,\\
551 0,&\text{otherwise},
552 \end{cases}
553 \qquad
554 \mathcal R=\mathcal B-\mathcal B_0.
555 \]
556 Symmetric pair retention preserves self-adjointness of $H_0$, although the
557 following error theorem does not require positive semidefiniteness of $H_0$.
558 
559 \begin{theorem}[Thresholded-power certificate]
560 \label{thm:thresholded-power}
561 For every $k\geq1$,
562 \[
563 \boxed{
564 \rho\bigl((H^k-H_0^k)X\bigr)
565 \leq
566 (\mathcal B^k-\mathcal B_0^k)\rho(X).
567 }
568 \]
569 In particular,
570 \[
571 \boxed{
572 \lVert(H^k-H_0^k)X\rVert_F
573 \leq
574 \left\lVert(\mathcal B^k-\mathcal B_0^k)\rho(X)\right\rVert_2.
575 }
576 \]
577 \end{theorem}
578 
579 \begin{proof}
580 Write $E=H-H_0$.  The telescoping identity gives
581 \[
582 H^k-H_0^k
583 =
584 \sum_{j=0}^{k-1}H^{k-1-j}EH_0^j.
585 \]
586 The block-norm matrix of $E$ is $\mathcal R$.  Applying
587 Theorem~\ref{thm:weighted-transfer} successively to each term gives
588 \[
589 \rho\bigl((H^k-H_0^k)X\bigr)
590 \leq
591 \sum_{j=0}^{k-1}
592 \mathcal B^{k-1-j}\mathcal R\mathcal B_0^j\rho(X).
593 \]
594 The matrix telescoping identity
595 \[
596 \mathcal B^k-\mathcal B_0^k
597 =
598 \sum_{j=0}^{k-1}
599 \mathcal B^{k-1-j}(\mathcal B-\mathcal B_0)\mathcal B_0^j
600 \]
601 yields the phasewise formula.  Orthogonality of the phase projectors gives
602 \[
603 \lVert Y\rVert_F^2=\sum_h\lVert P_hY\rVert_F^2,
604 \]
605 which proves the global bound.
606 \end{proof}
607 
608 \begin{corollary}[Certified approximate light cone]
609 \label{cor:approx-lightcone}
610 Assume $X$ is supported on a set of start phases $S$.  If phase $h$ is not
611 reachable from $S$ by a retained path of length at most $k$, then
612 \[
613 P_hH_0^kX=0
614 \]
615 and
616 \[
617 \boxed{
618 \lVert P_hH^kX\rVert_F
619 \leq
620 \bigl[(\mathcal B^k-\mathcal B_0^k)\rho(X)\bigr]_h.
621 }
622 \]
623 Thus an exact zero outside the retained light cone is replaced by a certified
624 nonzero leakage budget.
625 \end{corollary}
626 
627 \begin{corollary}[Scalar threshold law]
628 For any induced matrix norm, let
629 \[
630 M=\max\{\lVert\mathcal B\rVert,\lVert\mathcal B_0\rVert\}.
631 \]
632 Then
633 \[
634 \boxed{
635 \lVert\mathcal B^k-\mathcal B_0^k\rVert
636 \leq
637 kM^{k-1}\lVert\mathcal R\rVert.
638 }
639 \]
640 Hence weak discarded coupling produces a first-order error bound proportional
641 to the discarded transfer norm and to the Krylov depth.
642 \end{corollary}
643 
644 \begin{theorem}[Krylov-polynomial certificate]
645 \label{thm:polynomial-certificate}
646 Let
647 \[
648 p(t)=\sum_{j=0}^d a_jt^j.
649 \]
650 Then
651 \[
652 \boxed{
653 \rho\bigl((p(H)-p(H_0))X\bigr)
654 \leq
655 \sum_{j=1}^d
656 |a_j|(\mathcal B^j-\mathcal B_0^j)\rho(X).
657 }
658 \]
659 \end{theorem}
660 
661 \begin{proof}
662 Apply Theorem~\ref{thm:thresholded-power} to every nonconstant power and use
663 the triangle inequality.  The constant terms cancel.
664 \end{proof}
665 
666 The polynomial theorem covers arbitrary vectors in a finite Krylov space,
667 including polynomial filters and restarted combinations, rather than only the
668 monomial $H^kX$.
669 
670 \begin{corollary}[Certified phase pruning]
671 If a thresholded phase graph is chosen so that
672 \[
673 \left\lVert
674 (\mathcal B^k-\mathcal B_0^k)\rho(X)
675 \right\rVert_2
676 \leq\tau,
677 \]
678 then replacing $H^kX$ by $H_0^kX$ incurs Frobenius error at most $\tau$.
679 The same statement holds for a Krylov polynomial with the certificate in
680 Theorem~\ref{thm:polynomial-certificate}.
681 \end{corollary}
682 
683 This turns graph thresholding into an error-budget decision.  A phase or
684 interface is not discarded because it merely appears small; it is discarded
685 only when the weighted sum of all paths using discarded couplings remains
686 inside the requested tolerance.
687 
688 \begin{figure}[H]
689 \centering
690 \includegraphics[width=.84\linewidth]{figure_09_approximate_lightcone_certificate.png}
691 \caption{Weak long-range phase couplings replace the exact zero outside the retained light cone by a nonzero but certified envelope.  In this PSD test, the phasewise transfer-matrix certificate bounds the observed leakage at every depth from one through eight.}
692 \label{fig:approx-lightcone-certificate}
693 \end{figure}
694 
695 \begin{figure}[H]
696 \centering
697 \includegraphics[width=.84\linewidth]{figure_10_threshold_tradeoff.png}
698 \caption{Certified phase-graph sparsification at Krylov depth five.  Retaining fewer off-diagonal phase blocks reduces communication and operator work while the transfer-matrix bound remains above the measured truncation error.  The bound is intentionally conservative because it ignores cancellation and singular-vector alignment.}
699 \label{fig:threshold-tradeoff}
700 \end{figure}
701 
702 \subsection{Numerical audit}
703 
704 The validation operator has twelve phases of width four.  A block-bidiagonal
705 factor generates a positive-definite local operator $H_0$, and weak symmetric
706 long-range blocks are then added to obtain $H$.  Both normalized operators
707 remain positive definite.  A three-column start block is supported in phase
708 zero.
709 
710 For depths one through eight:
711 \begin{itemize}
712 \item every observed phase norm lay below $\mathcal B^k\rho(X)$;
713 \item every observed thresholding error lay below
714 $(\mathcal B^k-\mathcal B_0^k)\rho(X)$;
715 \item the largest total observed-error-to-certificate ratio was approximately
716 $0.669$.
717 \end{itemize}
718 
719 A separate hard-threshold sweep retained between $22$ and $132$ directed
720 off-diagonal phase blocks.  Every measured depth-five error remained below its
721 certificate; the largest observed-error-to-certificate ratio was approximately
722 $0.376$.  Finally, a degree-four signed polynomial in $H$ satisfied the
723 phasewise polynomial certificate in all twelve phases, with maximum ratio
724 approximately $0.068$.
725 
726 These computations validate the implementation and illustrate sharpness.  The
727 universal result is supplied by the preceding proofs, not by the finite test.
728 
729 
730 \section{Certified budget-optimal phase scheduling}
731 \label{sec:budget-scheduler}
732 
733 The preceding section answers a validation question: given a retained phase
734 operator $H_0$, how large can the Krylov error be?  This section turns that
735 certificate into a selection rule.  The goal is to remove as much phase
736 communication as possible while remaining inside a prescribed error budget.
737 
738 Let $\mathcal E$ denote the set of undirected off-diagonal phase pairs.  Because
739 $H$ is self-adjoint, $\mathcal B$ is symmetric.  For a phase pair
740 $e=\{h,g\}$ define the nonnegative matrix $\Delta_e$ by
741 \[
742 (\Delta_e)_{hg}=\mathcal B_{hg},\qquad
743 (\Delta_e)_{gh}=\mathcal B_{gh},
744 \]
745 with every other entry zero.  If $\mathcal O\subseteq\mathcal E$ is the set
746 of omitted pairs, define
747 \[
748 \mathcal B_{\mathcal O}
749 =\mathcal B-\sum_{e\in\mathcal O}\Delta_e.
750 \]
751 The diagonal is always retained.
752 
753 For depth $k\ge1$ and starting phase profile $\rho_0=\rho(X)$, define the
754 \emph{edge risk vector}
755 \[
756 \boxed{
757 q_e^{(k)}
758 =
759 \sum_{j=0}^{k-1}
760 \mathcal B^{k-1-j}\Delta_e\mathcal B^j\rho_0
761 }
762 \]
763 and its scalar risk
764 \[
765 \boxed{
766 \beta_e^{(k)}=\lVert q_e^{(k)}\rVert_2.
767 }
768 \]
769 Unlike the raw block norm $\mathcal B_{hg}$, this quantity depends on the
770 start profile and on Krylov depth.
771 
772 \begin{theorem}[Edge-additive leakage majorant]
773 \label{thm:edge-additive-majorant}
774 For every omitted set $\mathcal O\subseteq\mathcal E$,
775 \[
776 \boxed{
777 0\le
778 (\mathcal B^k-\mathcal B_{\mathcal O}^k)\rho_0
779 \le
780 \sum_{e\in\mathcal O}q_e^{(k)}
781 }
782 \]
783 componentwise.  Consequently,
784 \[
785 \boxed{
786 \lVert(H^k-H_{\mathcal O}^k)X\rVert_F
787 \le
788 \left\lVert
789 (\mathcal B^k-\mathcal B_{\mathcal O}^k)\rho_0
790 \right\rVert_2
791 \le
792 \sum_{e\in\mathcal O}\beta_e^{(k)}.
793 }
794 \]
795 \end{theorem}
796 
797 \begin{proof}
798 The first inequality in the second display is
799 Theorem~\ref{thm:thresholded-power}.  For the transfer matrices, the exact
800 noncommutative telescoping identity is
801 \[
802 \mathcal B^k-\mathcal B_{\mathcal O}^k
803 =
804 \sum_{j=0}^{k-1}
805 \mathcal B^{k-1-j}
806 (\mathcal B-\mathcal B_{\mathcal O})
807 \mathcal B_{\mathcal O}^j.
808 \]
809 All matrices and $\rho_0$ are nonnegative, and
810 $0\le\mathcal B_{\mathcal O}\le\mathcal B$ entrywise.  Hence
811 \[
812 \mathcal B_{\mathcal O}^j\rho_0
813 \le
814 \mathcal B^j\rho_0.
815 \]
816 Substituting
817 \[
818 \mathcal B-\mathcal B_{\mathcal O}
819 =
820 \sum_{e\in\mathcal O}\Delta_e
821 \]
822 and distributing the sum gives the componentwise majorant.  The final scalar
823 bound follows from the Euclidean triangle inequality.
824 \end{proof}
825 
826 \begin{corollary}[Polynomial edge-risk certificate]
827 Let
828 \[
829 p(t)=\sum_{d=0}^{D}a_dt^d.
830 \]
831 Define
832 \[
833 q_e^{(p)}
834 =
835 \sum_{d=1}^{D}|a_d|q_e^{(d)},
836 \qquad
837 \beta_e^{(p)}=\lVert q_e^{(p)}\rVert_2.
838 \]
839 Then
840 \[
841 \boxed{
842 \lVert(p(H)-p(H_{\mathcal O}))X\rVert_F
843 \le
844 \sum_{e\in\mathcal O}\beta_e^{(p)}.
845 }
846 \]
847 Thus the same scheduling machinery applies to finite Krylov combinations and
848 polynomial filters.
849 \end{corollary}
850 
851 \begin{corollary}[Query-invisible phase pair]
852 If
853 \[
854 \Delta_e\mathcal B^j\rho_0=0
855 \qquad
856 \text{for every }0\le j<k,
857 \]
858 then
859 \[
860 \beta_e^{(k)}=0.
861 \]
862 Therefore a phase pair can be large in operator norm and still have zero
863 risk for a particular finite-depth query.  A magnitude-only threshold cannot
864 represent this distinction.
865 \end{corollary}
866 
867 \subsection{Equal-cost optimality}
868 
869 Suppose each retained phase pair has the same communication or scheduling
870 cost.  Under the additive certificate, a set of omitted pairs is feasible if
871 \[
872 \sum_{e\in\mathcal O}\beta_e^{(k)}\le\tau.
873 \]
874 Sort the risks as
875 \[
876 \beta_{(1)}\le\beta_{(2)}\le\cdots\le\beta_{(N)}.
877 \]
878 
879 \begin{theorem}[Cardinality-optimal certified pruning]
880 \label{thm:unit-cost-optimal}
881 Let $r$ be the largest integer satisfying
882 \[
883 \sum_{i=1}^{r}\beta_{(i)}\le\tau.
884 \]
885 Then omitting the $r$ least-risk phase pairs is feasible, and no feasible
886 omitted set can contain more than $r$ pairs under the scalar additive
887 certificate.
888 \end{theorem}
889 
890 \begin{proof}
891 Among all subsets of cardinality $r+1$, the smallest possible sum of risks is
892 \[
893 \sum_{i=1}^{r+1}\beta_{(i)}.
894 \]
895 By maximality of $r$, that sum exceeds $\tau$.  Hence no set of $r+1$ omitted
896 pairs is feasible.  The prefix of length $r$ is feasible by construction.
897 \end{proof}
898 
899 This is a genuine optimization theorem, not a heuristic threshold rule.  The
900 optimization target is the additive certificate itself; the tighter exact
901 transfer envelope can leave additional slack.
902 
903 \subsection{Exact-envelope refinement}
904 
905 After the cardinality-optimal additive schedule is found, recomputing
906 \[
907 \left\lVert
908 (\mathcal B^k-\mathcal B_{\mathcal O}^k)\rho_0
909 \right\rVert_2
910 \]
911 may reveal unused budget.  Additional pairs can then be removed one at a time,
912 provided the exact transfer envelope remains below $\tau$.
913 
914 \begin{proposition}[Certified refinement]
915 Any sequence of additional omissions accepted only when
916 \[
917 \left\lVert
918 (\mathcal B^k-\mathcal B_{\mathcal O}^k)\rho_0
919 \right\rVert_2\le\tau
920 \]
921 remains a valid $\tau$-certificate for the original operator error.
922 \end{proposition}
923 
924 The refinement rule used in the reference implementation chooses, among all
925 currently feasible next omissions, the pair that produces the smallest
926 increase in the exact transfer envelope.  It is certified but is not claimed
927 globally optimal.
928 
929 \subsection{Unequal communication costs}
930 
931 Let $c_e\ge0$ denote the communication or storage cost saved by omitting pair
932 $e$.  Under the scalar additive certificate, the exact discrete selection
933 problem is
934 \[
935 \boxed{
936 \max_{x_e\in\{0,1\}}
937 \sum_e c_ex_e
938 \quad\text{subject to}\quad
939 \sum_e\beta_e^{(k)}x_e\le\tau.
940 }
941 \]
942 This is the $0$--$1$ knapsack problem.  Hence unequal-cost scheduling is
943 combinatorial in general.  The implementation supplies a safe
944 cost-to-risk greedy heuristic; its output remains certified because it never
945 exceeds the additive risk budget.  A fractional-knapsack relaxation supplies
946 an upper bound on the maximum possible saved cost.
947 
948 \begin{figure}[H]
949 \centering
950 \includegraphics[width=.84\linewidth]{figure_11_budget_scheduler_comparison.png}
951 \caption{Depth-five pruning under a fixed certified error budget.  A global block-norm threshold is nearly unable to prune this example without crossing the requested tolerance, whereas query-aware edge risks remove many phase pairs.  The exact-envelope refinement safely uses slack left by the additive majorant.}
952 \label{fig:budget-scheduler-comparison}
953 \end{figure}
954 
955 \begin{figure}[H]
956 \centering
957 \includegraphics[width=.82\linewidth]{figure_13_edge_risk_vs_block_norm.png}
958 \caption{Raw phase-block magnitude and finite-depth query risk are different quantities.  The edge score $\beta_{hg}^{(5)}$ incorporates the start profile and every depth-five transfer path through the phase graph; the block norm $\mathcal B_{hg}$ does not.}
959 \label{fig:edge-risk-vs-norm}
960 \end{figure}
961 
962 \subsection{Numerical audit}
963 
964 The twelve-phase operator from Section~\ref{sec:approx-lightcone} was reused at
965 Krylov depth five.  There are $66$ undirected off-diagonal phase pairs.
966 For error budgets $10^{-3}$, $2\times10^{-3}$, $5\times10^{-3}$, and
967 $10^{-2}$, the additive-optimal scheduler retained respectively
968 $61$, $56$, $45$, and $36$ pairs.  Exact-envelope refinement reduced those
969 counts further to
970 $58$, $50$, $36$, and $30$ pairs while remaining inside budget.
971 
972 In the same test, the best single global threshold on the raw block norm
973 retained $65$ pairs at every one of those four budgets.  The difference is
974 structural: the global threshold is blind to start phase and Krylov depth.
975 
976 A separate six-phase problem contains $15$ candidate phase pairs, so all
977 $2^{15}=32768$ subsets can be enumerated.  At four distinct budgets, the
978 sorted-risk scheduler matched the brute-force cardinality optimum under the
979 additive certificate in every case.  The exact-envelope greedy refinement also
980 matched the brute-force exact-envelope optimum in cardinality for all four
981 budgets in this finite audit.  The latter agreement is empirical and is not
982 promoted to a theorem.
983 
984 For an unequal-cost eight-phase test, costs were proportional to the bytes in
985 symmetric dense phase blocks.  The cost-to-risk heuristic saved $5408$ bytes
986 out of a $9472$-byte communication proxy while using additive risk $0.011632$
987 under budget $0.012$.  The fractional relaxation upper bound was $5647.8$
988 bytes, placing the greedy schedule within about $4.25\%$ of that relaxation
989 bound on this instance.
990 
991 \section{Asynchronous double-buffer execution model}
992 \label{sec:double-buffer}
993 
994 The mathematical streaming identities do not require transfers and GEMMs to be
995 serial.  Consider $B$ storage blocks.  Let $t_b$ be the transfer time for block
996 $b$ and $c_b$ its device compute time.
997 
998 \begin{proposition}[Ideal two-buffer pipeline]
999 In a model with perfect asynchronous overlap, transfer of block $b+1$ may run
1000 while block $b$ is being processed.  The serial one-buffer time is
1001 \[
1002 T_{\rm serial}=\sum_{b=1}^{B}(t_b+c_b),
1003 \]
1004 whereas the two-buffer pipeline time is
1005 \[
1006 \boxed{
1007 T_{\rm pipe}
1008 =
1009 t_1
1010 +
1011 \sum_{b=1}^{B-1}\max\{c_b,t_{b+1}\}
1012 +
1013 c_B.
1014 }
1015 \]
1016 For uniform times $t_b=t$ and $c_b=c$,
1017 \[
1018 T_{\rm pipe}=t+c+(B-1)\max\{t,c\}.
1019 \]
1020 \end{proposition}
1021 
1022 \begin{proof}
1023 After the first block arrives, every interior pipeline stage finishes only when
1024 both the current computation and the next transfer have completed.  Their
1025 joint stage duration is therefore the maximum of those two times.  The initial
1026 transfer and final compute cannot be hidden and form the two boundary terms.
1027 \end{proof}
1028 
1029 If one resident block uses $M_{\rm block}$ bytes and all other fixed working
1030 state uses $M_{\rm fixed}$ bytes, the corresponding memory models are
1031 \[
1032 M_{\rm single}=M_{\rm fixed}+M_{\rm block},
1033 \qquad
1034 \boxed{
1035 M_{\rm double}=M_{\rm fixed}+2M_{\rm block}.
1036 }
1037 \]
1038 The proposition is a scheduling model, not a hardware timing claim.  PCIe or
1039 NVLink contention, kernel launch overhead, synchronization, page pinning, and
1040 library workspaces can reduce realized overlap.
1041 
1042 \begin{figure}[H]
1043 \centering
1044 \includegraphics[width=.82\linewidth]{figure_12_double_buffer_model.png}
1045 \caption{Idealized effect of asynchronous double buffering over $64$ streamed blocks.  The largest modeled speedup occurs when transfer and computation times are comparable, because neither resource is left idle for most of the pipeline.  The additional cost is one extra resident block buffer.}
1046 \label{fig:double-buffer-model}
1047 \end{figure}
1048 
1049 \section{Phase-local rank}
1050 
1051 Let \(R\) be a newly generated Krylov residual and
1052 \[
1053 R_h=P_hR.
1054 \]
1055 
1056 \begin{theorem}[Minimal exact phase width]
1057 Let
1058 \[
1059 r_h=\operatorname{rank}(R_h).
1060 \]
1061 Every phase-supported basis spanning \(\operatorname{range}(R_h)\) has at
1062 least \(r_h\) vectors, and a rank-revealing QR or compact SVD attains exactly
1063 \(r_h\).
1064 \end{theorem}
1065 
1066 This justifies nonuniform Krylov block widths by phase.
1067 
1068 
1069 \begin{figure}[H]
1070 \centering
1071 \includegraphics[width=.82\linewidth]{figure_05_phase_local_rank.png}
1072 \caption{Phase-local rank as a dimension law. The singular values of a phase residual reveal how many genuinely new directions that phase contributes. A rank-revealing factorization keeps exactly those directions instead of allocating a uniform global width.}
1073 \label{fig:phase-local-rank}
1074 \end{figure}
1075 
1076 \section{Low-rank interface messages}
1077 
1078 If an inter-phase coupling has exact factorization
1079 \[
1080 P_hHP_g
1081 =
1082 U_{hg}V_{hg}^\top
1083 \]
1084 of rank \(r_{hg}\), then
1085 \[
1086 \boxed{
1087 P_hHP_gx
1088 =
1089 U_{hg}(V_{hg}^\top x).
1090 }
1091 \]
1092 
1093 This representation reduces storage and communication.  It does not
1094 automatically reduce wall-clock time unless the factors are applied through
1095 batched or fused kernels; a Python reference implementation can lose to
1096 dense BLAS despite lower arithmetic and storage.
1097 
1098 \section{Exact QR-core updates}
1099 
1100 Suppose
1101 \[
1102 A=QR,
1103 \qquad
1104 Q^\top Q=I,
1105 \]
1106 and append a new block \(U\).  Set
1107 \[
1108 C=Q^\top U,
1109 \qquad
1110 F=U-QC,
1111 \qquad
1112 F=Q_\perp R_\perp.
1113 \]
1114 
1115 \begin{theorem}[Exact QR-core update]
1116 \[
1117 [A,U]
1118 =
1119 [Q,Q_\perp]
1120 \begin{pmatrix}
1121 R&C\\
1122 0&R_\perp
1123 \end{pmatrix}.
1124 \]
1125 Therefore the nonzero eigenvalues of
1126 \([A,U][A,U]^\top\) are exactly the eigenvalues of the small core Gram.
1127 \end{theorem}
1128 
1129 This is useful for overlapping genomic windows and small incremental variant
1130 updates.
1131 
1132 
1133 \begin{figure}[H]
1134 \centering
1135 \includegraphics[width=.84\linewidth]{figure_06_qr_core_equivalence.png}
1136 \caption{QR-core spectral equivalence. The singular values of the full updated factor and the compact update core lie on top of one another. This permits exact reuse across overlapping windows or appended variant blocks.}
1137 \label{fig:qr-core-equivalence}
1138 \end{figure}
1139 
1140 \section{Range-core spectral compression}
1141 
1142 Suppose an orthonormal basis \(Q_r\) captures the factor range:
1143 \[
1144 \operatorname{range}(A)\subseteq\operatorname{range}(Q_r).
1145 \]
1146 Then
1147 \[
1148 A=Q_rB,
1149 \qquad
1150 B=Q_r^\top A.
1151 \]
1152 
1153 \begin{theorem}[Exact range core]
1154 The nonzero eigenvalues of
1155 \[
1156 K=AA^\top
1157 \]
1158 are exactly the eigenvalues of
1159 \[
1160 BB^\top.
1161 \]
1162 \end{theorem}
1163 
1164 Thus the computational pipeline can terminate its data-streaming stage after
1165 the range is captured:
1166 \[
1167 \boxed{
1168 \text{streamed data}
1169 \to
1170 \text{small spectral core}
1171 \to
1172 \text{tail evaluation}.
1173 }
1174 \]
1175 
1176 \section{From phase pruning to certified SKAT tails}
1177 
1178 The previous sections certify the effect of phase pruning on Krylov
1179 propagation.  A SKAT-like test ultimately depends on the spectrum itself, so
1180 one more bridge is needed:
1181 \[
1182 \text{phase pruning}
1183 \longrightarrow
1184 \text{operator error}
1185 \longrightarrow
1186 \text{spectral error}
1187 \longrightarrow
1188 \text{tail-probability error}.
1189 \]
1190 
1191 \subsection{Block-norm operator domination}
1192 
1193 Let
1194 \[
1195 E=H-H_0
1196 \]
1197 be the self-adjoint error introduced by phase pruning.  Using the same
1198 orthogonal phase projectors \(P_h\), define
1199 \[
1200 \mathcal C_{hg}
1201 =
1202 \|P_hEP_g\|_2.
1203 \]
1204 
1205 \begin{theorem}[Block-norm operator certificate]
1206 \[
1207 \boxed{
1208 \|E\|_2
1209 \le
1210 \|\mathcal C\|_2.
1211 }
1212 \]
1213 \end{theorem}
1214 
1215 \begin{proof}
1216 For any vector \(x\), let
1217 \[
1218 r_g=\|P_gx\|_2.
1219 \]
1220 Then
1221 \[
1222 \|P_hEx\|_2
1223 \le
1224 \sum_g
1225 \|P_hEP_g\|_2\,\|P_gx\|_2
1226 =
1227 (\mathcal C r)_h.
1228 \]
1229 Orthogonality of the phase decomposition gives
1230 \[
1231 \|Ex\|_2^2
1232 =
1233 \sum_h\|P_hEx\|_2^2
1234 \le
1235 \|\mathcal C r\|_2^2
1236 \le
1237 \|\mathcal C\|_2^2\|r\|_2^2
1238 =
1239 \|\mathcal C\|_2^2\|x\|_2^2.
1240 \]
1241 Taking the supremum over unit vectors proves the claim.
1242 \end{proof}
1243 
1244 The important computational point is that \(\mathcal C\) is only a
1245 phase-by-phase matrix.  The same block norms already used to construct the
1246 polyphase scheduler therefore provide a global operator certificate.
1247 
1248 \subsection{Weyl enclosure of the quadratic-form coefficients}
1249 
1250 Assume
1251 \[
1252 H\succeq0
1253 \]
1254 and let the eigenvalues of \(H\) and \(H_0\), in decreasing order, be
1255 \[
1256 \lambda_i
1257 \qquad\text{and}\qquad
1258 \mu_i.
1259 \]
1260 If
1261 \[
1262 \|H-H_0\|_2\le\varepsilon,
1263 \]
1264 Weyl's inequality gives
1265 \[
1266 |\lambda_i-\mu_i|\le\varepsilon.
1267 \]
1268 
1269 \begin{corollary}[Positive spectral interval]
1270 For every \(i\),
1271 \[
1272 \boxed{
1273 \ell_i
1274 =
1275 \max\{0,\mu_i-\varepsilon\}
1276 \le
1277 \lambda_i
1278 \le
1279 u_i
1280 =
1281 \max\{0,\mu_i+\varepsilon\}.
1282 }
1283 \]
1284 \end{corollary}
1285 
1286 This remains valid even if symmetric phase pruning makes \(H_0\) slightly
1287 indefinite.
1288 
1289 \subsection{Certified generalized-\texorpdfstring{\(\chi^2\)}{chi-square} tail bracket}
1290 
1291 For independent standard normal variables \(Z_i\), define
1292 \[
1293 Q_-=\sum_i\ell_iZ_i^2,
1294 \qquad
1295 Q=\sum_i\lambda_iZ_i^2,
1296 \qquad
1297 Q_+=\sum_i u_iZ_i^2.
1298 \]
1299 
1300 \begin{theorem}[Stochastic tail enclosure]
1301 For every threshold \(q\),
1302 \[
1303 \boxed{
1304 \Pr(Q_-\ge q)
1305 \le
1306 \Pr(Q\ge q)
1307 \le
1308 \Pr(Q_+\ge q).
1309 }
1310 \]
1311 \end{theorem}
1312 
1313 \begin{proof}
1314 Couple all three quadratic forms using the same normal variables \(Z_i\).
1315 Since
1316 \[
1317 0\le\ell_i\le\lambda_i\le u_i
1318 \]
1319 for every \(i\),
1320 \[
1321 Q_-\le Q\le Q_+
1322 \]
1323 almost surely.  Taking upper-tail events yields the result.
1324 \end{proof}
1325 
1326 The two endpoint probabilities can be evaluated by the same K-COS machinery
1327 used for the unpruned spectrum.  Its numerical error certificate is separate
1328 from the spectral enclosure above.
1329 
1330 \subsection{A dual-budget pruning rule}
1331 
1332 The local Krylov certificate from the previous section and the global spectral
1333 certificate answer different questions.  The former asks whether omitted
1334 couplings matter to the current Krylov query; the latter asks whether they can
1335 move the spectrum enough to matter to the final test.
1336 
1337 A production schedule should therefore enforce
1338 \[
1339 \boxed{
1340 \sum_{e\in\mathcal O}\beta_e^{(k)}
1341 \le
1342 \tau_K,
1343 \qquad
1344 \|\mathcal C_{\mathcal O}\|_2
1345 \le
1346 \tau_\lambda.
1347 }
1348 \]
1349 
1350 The implementation uses a conservative greedy rule: phase pairs are considered
1351 in increasing query-risk order, and a pair is omitted only when adding it
1352 preserves both inequalities.  The rule is certified but is not claimed to be
1353 globally optimal under the two simultaneous constraints.
1354 
1355 \begin{figure}[H]
1356 \centering
1357 \includegraphics[width=.84\linewidth]{figure_14_dual_budget_tail_bridge.png}
1358 \caption{Dual-budget phase pruning.  The query-local budget controls Krylov
1359 propagation while the operator budget controls final spectral and tail
1360 uncertainty.  A loose operator budget can admit a strong coupling and widen
1361 the tail interval by several decades even when the Krylov-local budget still
1362 passes.}
1363 \label{fig:dual-budget-tail}
1364 \end{figure}
1365 
1366 In the deterministic ten-phase experiment, fixing the query budget at
1367 \(0.05\) and the operator budget at \(0.005\) removed \(23\) of \(45\)
1368 undirected phase pairs.  At a true probability of \(10^{-8}\), the certified
1369 tail interval was
1370 \[
1371 10^{-8.11143}
1372 \le
1373 p
1374 \le
1375 10^{-7.89036},
1376 \]
1377 a width of about \(0.221\) decades.  With an operator budget of \(0.0075\),
1378 \(34\) pairs were removed and the interval widened to about \(0.320\)
1379 decades.  By contrast, allowing a strong coupling to be removed under a very
1380 loose operator budget produced a multi-decade interval.
1381 
1382 \subsection{Certify after prune}
1383 
1384 There is a second route that does not use the pruned operator as the final
1385 spectral approximation.  Let \(U\) be any orthonormal sample-space basis,
1386 possibly generated cheaply from a pruned or approximate Krylov process, and
1387 define
1388 \[
1389 B=U^\top A,
1390 \qquad
1391 F=(I-UU^\top)A,
1392 \qquad
1393 \delta=\|F\|_2.
1394 \]
1395 
1396 \begin{theorem}[Range-core singular-value certificate]
1397 Let the singular values of \(A\) be \(\sigma_i(A)\), and let
1398 \(\widetilde\sigma_i(B)\) denote the singular values of \(B\) padded by zeros
1399 to the ambient spectral dimension.  Then
1400 \[
1401 \boxed{
1402 \left|
1403 \sigma_i(A)
1404 -
1405 \widetilde\sigma_i(B)
1406 \right|
1407 \le
1408 \delta.
1409 }
1410 \]
1411 Consequently,
1412 \[
1413 \max\{0,\widetilde\sigma_i(B)-\delta\}^2
1414 \le
1415 \lambda_i(AA^\top)
1416 \le
1417 \left(
1418 \widetilde\sigma_i(B)+\delta
1419 \right)^2.
1420 \]
1421 \end{theorem}
1422 
1423 \begin{proof}
1424 Because \(U\) has orthonormal columns,
1425 \[
1426 A=UB+F
1427 \]
1428 and
1429 \[
1430 \|A-UB\|_2=\|F\|_2=\delta.
1431 \]
1432 The nonzero singular values of \(UB\) are exactly those of \(B\).
1433 The standard singular-value perturbation inequality gives the stated bound.
1434 Squaring the nonnegative endpoints gives the eigenvalue interval.
1435 \end{proof}
1436 
1437 This provides a useful architecture:
1438 \[
1439 \boxed{
1440 \text{pruned operator for basis discovery}
1441 \to
1442 \text{one exact streamed factor check}
1443 \to
1444 \text{certified spectral core}.
1445 }
1446 \]
1447 
1448 \subsection{A retained negative result}
1449 
1450 The preceding theorem is exact, but a correct interval need not be useful.
1451 
1452 In a full-rank \(50\times50\) stress test, a pruned Krylov basis of dimensions
1453 \(12\) through \(48\) left range residuals between roughly \(0.87\) and
1454 \(0.99\) times the leading singular value.  The corresponding \(10^{-8}\)
1455 tail intervals were therefore vacuous or thousands of decades wide.  At
1456 dimension \(50\), the range residual fell to approximately
1457 \(1.2\times10^{-15}\), and the tail interval collapsed to the exact value.
1458 
1459 \begin{figure}[H]
1460 \centering
1461 \includegraphics[width=.84\linewidth]{figure_15_range_core_certification.png}
1462 \caption{Certify-after-prune range residual.  The range-core theorem is valid
1463 at every depth, but a deep-tail full-spectrum guarantee can remain useless
1464 until the range residual is genuinely small.  This negative result is retained
1465 as a stopping-rule lesson rather than hidden as an implementation failure.}
1466 \label{fig:range-core-cert}
1467 \end{figure}
1468 
1469 The practical lesson is that range-core certification is a stopping test, not
1470 a promise that a severely incomplete basis automatically yields a useful
1471 full-spectrum tail bound.
1472 
1473 
1474 \section{Trace-aware residual certification and exact completion}
1475 
1476 The spectral bridge in the preceding section used an operator-norm residual as
1477 if every unresolved coefficient could move independently by the same amount.
1478 That is rigorous but unnecessarily pessimistic once a nontrivial range basis
1479 has already been captured.  The residual has considerably more structure.
1480 
1481 Let \(U\in\mathbb R^{n\times r}\) have orthonormal columns and define
1482 \[
1483 B=U^\top A,
1484 \qquad
1485 F=(I-UU^\top)A.
1486 \]
1487 Then
1488 \[
1489 A=UB+F,
1490 \qquad
1491 U^\top F=0.
1492 \]
1493 
1494 \begin{theorem}[Exact core--residual decomposition]
1495 \[
1496 \boxed{
1497 A^\top A
1498 =
1499 B^\top B+F^\top F.
1500 }
1501 \]
1502 In particular the unresolved spectral term
1503 \[
1504 R=F^\top F
1505 \]
1506 is positive semidefinite.
1507 \end{theorem}
1508 
1509 \begin{proof}
1510 Expand \((UB+F)^\top(UB+F)\).  The two cross terms vanish because
1511 \(U^\top F=0\), while \(U^\top U=I\).
1512 \end{proof}
1513 
1514 This changes the certification problem.  Instead of an arbitrary Hermitian
1515 error, the missing term is a positive residual with measurable spectral
1516 moments.  Define
1517 \[
1518 t=\operatorname{tr}R,
1519 \qquad
1520 v=\operatorname{tr}(R^2),
1521 \qquad
1522 \rho=\|R\|_2.
1523 \]
1524 Then
1525 \[
1526 \boxed{
1527  t=\|F\|_F^2,
1528  \qquad
1529  v=\|FF^\top\|_F^2,
1530  \qquad
1531  \rho=\|F\|_2^2.
1532 }
1533 \]
1534 The first quantity is the exact missing mean of the Gaussian quadratic form;
1535 the second records how concentrated the residual spectrum is; the third is its
1536 largest possible coefficient.
1537 
1538 \subsection{A trace--Frobenius--cap MGF envelope}
1539 
1540 Set
1541 \[
1542 C=B^\top B,
1543 \qquad
1544 H=C+R,
1545 \qquad
1546 c=\lambda_{\max}(C).
1547 \]
1548 For \(0<2\theta(c+\rho)<1\), put
1549 \[
1550 \alpha=(1-2\theta c)^{-1}.
1551 \]
1552 
1553 \begin{theorem}[Trace--Frobenius--cap log-MGF certificate]
1554 For the Gaussian quadratic form with coefficient matrix \(H\),
1555 \[
1556 K_H(\theta)
1557 =-\frac12\log\det(I-2\theta H)
1558 \]
1559 satisfies
1560 \[
1561 \boxed{
1562 K_H(\theta)
1563 \le
1564 K_C(\theta)
1565 +
1566 \theta\alpha t
1567 +
1568 \frac{\theta^2\alpha^2v}
1569      {1-2\theta\alpha\rho}.
1570 }
1571 \]
1572 \end{theorem}
1573 
1574 \begin{proof}
1575 Factor the determinant through the known core:
1576 \[
1577 I-2\theta(C+R)
1578 =
1579 (I-2\theta C)^{1/2}
1580 (I-2\theta T)
1581 (I-2\theta C)^{1/2},
1582 \]
1583 where
1584 \[
1585 T=(I-2\theta C)^{-1/2}
1586 R
1587 (I-2\theta C)^{-1/2}.
1588 \]
1589 Since
1590 \[
1591 \|(I-2\theta C)^{-1/2}\|_2^2=\alpha,
1592 \]
1593 we have
1594 \[
1595 \operatorname{tr}T\le\alpha t,
1596 \qquad
1597 \operatorname{tr}(T^2)\le\alpha^2v,
1598 \qquad
1599 \|T\|_2\le\alpha\rho.
1600 \]
1601 For every eigenvalue \(x\) of \(T\),
1602 \[
1603 -\frac12\log(1-2\theta x)
1604 \le
1605 \theta x
1606 +
1607 \frac{\theta^2x^2}
1608      {1-2\theta\alpha\rho}.
1609 \]
1610 The scalar inequality follows from the power series of \(-\log(1-u)\) and a
1611 geometric bound on all terms of degree at least two.  Summing over the
1612 eigenvalues of \(T\) proves the result.
1613 \end{proof}
1614 
1615 \begin{corollary}[Trace-aware Chernoff upper tail]
1616 For \(Q=Z^\top HZ\) with \(Z\sim N(0,I)\),
1617 \[
1618 \boxed{
1619 \Pr(Q\ge q)
1620 \le
1621 \inf_{0<2\theta(c+\rho)<1}
1622 \exp\{\widehat K_H(\theta)-\theta q\},
1623 }
1624 \]
1625 where \(\widehat K_H\) is the preceding MGF upper bound.
1626 \end{corollary}
1627 
1628 A complementary lower endpoint is immediate from
1629 \[
1630 H=C+R\succeq C.
1631 \]
1632 Using the same normal vector,
1633 \[
1634 Z^\top HZ\ge Z^\top CZ
1635 \]
1636 almost surely, and therefore
1637 \[
1638 \boxed{
1639 \Pr(Z^\top CZ\ge q)
1640 \le
1641 \Pr(Z^\top HZ\ge q).
1642 }
1643 \]
1644 The lower endpoint can be evaluated by K-COS on the captured core spectrum,
1645 while the upper endpoint is obtained from the optimized trace-aware Chernoff
1646 bound.
1647 
1648 \subsection{Finite stress test}
1649 
1650 The previous coefficient-wise Weyl interval was intentionally conservative.
1651 On the same full-rank, 50-dimensional stress test used in the preceding
1652 section, its \(10^{-8}\) tail interval remained thousands of decades wide
1653 until essentially the complete range had been captured.  Preserving the exact
1654 residual trace and higher residual moments changes that behavior substantially.
1655 
1656 \begin{table}[H]
1657 \centering
1658 \caption{Tail-bracket width at a true probability of \(10^{-8}\).  The
1659 trace/moment certificate is still conservative, but it is orders of magnitude
1660 narrower than coefficient-wise Weyl inflation before full range capture.}
1661 \begin{tabular}{rrrr}
1662 \toprule
1663 Krylov depth & basis dimension & v0.7 Weyl width & v0.8 moment width\\
1664 \midrule
1665 4 & 24 & 90,097 & 11.86\\
1666 5 & 30 & 16,511 & 9.76\\
1667 6 & 36 & 3,747 & 7.67\\
1668 7 & 42 & 2,275 & 5.70\\
1669 8 & 48 & 1,165 & 3.39\\
1670 \bottomrule
1671 \end{tabular}
1672 \end{table}
1673 
1674 \begin{figure}[H]
1675 \centering
1676 \includegraphics[width=.84\linewidth]{figure_16_trace_moment_vs_weyl.png}
1677 \caption{Incomplete-range tail certification before exact completion.  The
1678 trace/moment bound preserves the amount and concentration of unresolved
1679 spectral mass instead of shifting every coefficient by one worst-case radius.}
1680 \label{fig:trace-moment-vs-weyl}
1681 \end{figure}
1682 
1683 A separate random implementation audit used 120 random factor/range pairs and
1684 four admissible positive tilts per pair, for 480 MGF comparisons.  Both the
1685 trace/cap and trace--Frobenius--cap envelopes dominated the exact log MGF in
1686 every tested case.  The trace--Frobenius--cap form was at least as tight as the
1687 trace/cap form in approximately 83.5\% of those comparisons.  This finite
1688 sweep audits the implementation; the theorem above supplies the universal
1689 argument.
1690 
1691 \subsection{Residual rank from moments}
1692 
1693 The same invariants determine when it is better to stop approximating the
1694 missing spectrum and complete it exactly.  Let
1695 \[
1696 s=\operatorname{rank}F=\operatorname{rank}R.
1697 \]
1698 If the nonzero eigenvalues of \(R\) are \(r_1,\ldots,r_s\), then
1699 \[
1700 t=\sum_i r_i,
1701 \qquad
1702 v=\sum_i r_i^2,
1703 \qquad
1704 r_i\le\rho.
1705 \]
1706 
1707 \begin{theorem}[Residual-rank bracket]
1708 \[
1709 \boxed{
1710 \left\lceil
1711 \max\left\{
1712 \frac{t}{\rho},
1713 \frac{t^2}{v}
1714 \right\}
1715 \right\rceil
1716 \le
1717 \operatorname{rank}F
1718 \le
1719 n-r.
1720 }
1721 \]
1722 The zero-residual case is interpreted separately as rank zero.
1723 \end{theorem}
1724 
1725 \begin{proof}
1726 The inequality \(t\le s\rho\) gives \(s\ge t/\rho\).  Cauchy--Schwarz gives
1727 \[
1728 t^2
1729 =
1730 \left(\sum_{i=1}^s r_i\right)^2
1731 \le
1732 s\sum_{i=1}^s r_i^2
1733 =sv,
1734 \]
1735 so \(s\ge t^2/v\).  Finally, every column of \(F\) lies in
1736 \(\operatorname{range}(U)^\perp\), whose dimension is \(n-r\).
1737 \end{proof}
1738 
1739 Whenever the integer lower bound equals \(n-r\), the missing rank is certified
1740 exactly without computing a residual SVD.
1741 
1742 \begin{figure}[H]
1743 \centering
1744 \includegraphics[width=.84\linewidth]{figure_17_residual_rank_completion.png}
1745 \caption{Residual-rank certification in the full-rank stress test.  From depth
1746 3 onward the moment lower bound meets the orthogonal-complement upper bound.
1747 At depth 8 exactly two residual directions remain.}
1748 \label{fig:residual-rank-completion}
1749 \end{figure}
1750 
1751 \subsection{Exact residual completion}
1752 
1753 Suppose a compact factorization of the residual is
1754 \[
1755 F=WD,
1756 \]
1757 where \(W\) has orthonormal columns.  Since \(F\) lies in the orthogonal
1758 complement of \(U\),
1759 \[
1760 U^\top W=0.
1761 \]
1762 
1763 \begin{theorem}[Exact residual completion]
1764 \[
1765 \boxed{
1766 A
1767 =
1768 [U\;W]
1769 \begin{bmatrix}
1770 B\\D
1771 \end{bmatrix}.
1772 }
1773 \]
1774 Consequently, the nonzero singular values of \(A\) are exactly the singular
1775 values of the stacked completion core
1776 \[
1777 \mathcal C_{\rm complete}
1778 =
1779 \begin{bmatrix}
1780 B\\D
1781 \end{bmatrix}.
1782 \]
1783 \end{theorem}
1784 
1785 This identity admits a column-streamed implementation.  For each residual
1786 column block, project against the residual basis accumulated so far, append
1787 only genuinely new directions through rank-revealing QR or SVD, and extend the
1788 coordinate matrices of earlier blocks with zero rows.  New residual directions
1789 are orthogonal to every earlier residual block by construction, so no earlier
1790 data block must be revisited.
1791 
1792 In the deterministic stress test, exact residual completion reproduced the
1793 full singular spectrum at relative errors between approximately
1794 \(2\times10^{-16}\) and \(6\times10^{-16}\) for all tested incomplete bases.
1795 The resulting \(10^{-8}\) tail differed from the full-spectrum reference by
1796 less than \(10^{-14}\) decades.  At depth 8 the candidate range had dimension
1797 48, the missing rank was certified as exactly two, and appending those two
1798 residual directions closed the spectrum without another approximate Krylov
1799 stage.
1800 
1801 \subsection{Hybrid stopping rule}
1802 
1803 The resulting production logic is not ``approximate forever.''  It is a
1804 certified handoff between two regimes:
1805 \[
1806 \boxed{
1807 \begin{array}{c}
1808 \text{streamed / pruned Krylov}\\[2mm]
1809 \downarrow\\[-1mm]
1810 \text{residual moments }(t,v,\rho)\\[2mm]
1811 \downarrow\\[-1mm]
1812 \begin{cases}
1813 \text{continue Krylov}, & \text{residual range still large},\\
1814 \text{exact streamed residual completion}, & \text{completion rank affordable}.
1815 \end{cases}
1816 \end{array}
1817 }
1818 \]
1819 
1820 The trace-aware tail certificate remains useful as a rigorous intermediate
1821 bound and as a stopping diagnostic.  Once only a small number of residual
1822 directions remain, exact completion is mathematically preferable to forcing a
1823 worst-case tail envelope to become artificially sharp.
1824 
1825 \section{Realistic rare-variant SKAT simulation validation}
1826 
1827 The preceding sections establish operator identities and certified numerical
1828 bounds.  This section tests the same machinery on the project's existing
1829 rare-variant SKAT simulation layer, including explicit LD structure,
1830 covariates, per-variant standardization, rare-variant weighting, and the
1831 phenotype-space score statistic.
1832 
1833 \subsection{Simulation provenance and a corrected historical defect}
1834 
1835 Two established genotype generators were used.
1836 
1837 The first is a latent-Gaussian block-LD generator.  Within a genomic block,
1838 alleles are obtained by thresholding correlated Gaussian latent variables at
1839 variant-specific minor-allele frequencies.
1840 
1841 The second is a copy-LD generator in which haplotype states persist along a
1842 block with high probability and are occasionally redrawn.  During the v0.9
1843 audit, a defect was found in the historical implementation: the genotype
1844 matrix was allocated with \verb|np.empty| and then updated using \verb|+=|.
1845 Because the initial contents of \verb|np.empty| are unspecified, those old
1846 copy-LD benchmark values are not reproducible evidence.  The v0.9 generator
1847 uses an explicitly zero-initialized matrix before haplotype accumulation.
1848 All copy-LD results reported below use the corrected implementation.
1849 
1850 \subsection{Tested LD regimes}
1851 
1852 Four deterministic scenarios were evaluated.
1853 
1854 \begin{table}[H]
1855 \centering
1856 \caption{Realistic rare-variant simulation scenarios.  Adjacent $r^2$ is
1857 computed within the simulated LD blocks after variant standardization.}
1858 \begin{tabular}{lrrrrr}
1859 \toprule
1860 scenario & $n$ & $m$ & rank & median adjacent $r^2$ & effective rank\\
1861 \midrule
1862 latent LD, moderate & 600 & 3000 & 596 & 0.0069 & 441.9\\
1863 latent LD, strong & 600 & 3000 & 596 & 0.1140 & 240.3\\
1864 copy LD, strong & 600 & 3000 & 596 & 0.8298 & 200.1\\
1865 copy LD, large & 1000 & 5000 & 996 & 0.8179 & 353.0\\
1866 \bottomrule
1867 \end{tabular}
1868 \end{table}
1869 
1870 The corresponding spectra therefore span a substantial range of concentration:
1871 stronger LD reduces effective rank even when the nominal quadratic-form rank
1872 remains close to the sample residual dimension.
1873 
1874 \begin{figure}[H]
1875 \centering
1876 \includegraphics[width=.84\linewidth]{figure_18_realistic_ld_effective_rank.png}
1877 \caption{LD structure changes the effective spectral rank.  The corrected
1878 copy-LD models create much stronger local dependence than the latent-Gaussian
1879 moderate-LD case, while all remain within the same streamed SKAT pipeline.}
1880 \label{fig:realistic-ld-rank}
1881 \end{figure}
1882 
1883 \subsection{Streaming directly from raw genotype blocks}
1884 
1885 The raw genotype matrix was streamed in resident variant blocks of widths
1886 31, 127, and 509.  Within each block the implementation performs
1887 
1888 \[
1889 G_b
1890 \to
1891 G_b-\bar G_b
1892 \to
1893 \text{variance scaling}
1894 \to
1895 M_X G_b
1896 \to
1897 W_b^{1/2}
1898 \to
1899 A_bA_b^\top,
1900 \]
1901 
1902 without materializing the full weighted factor.
1903 
1904 Across all four LD scenarios and all three resident block sizes, the largest
1905 relative discrepancy against the fully materialized kernel was
1906 
1907 \[
1908 \boxed{6.61\times10^{-16}}.
1909 \]
1910 
1911 \begin{figure}[H]
1912 \centering
1913 \includegraphics[width=.84\linewidth]{figure_19_realistic_streaming_invariance.png}
1914 \caption{Raw-genotype streaming reproduces the materialized SKAT kernel across
1915 multiple LD regimes and block sizes.  The observed discrepancies are at
1916 floating-point summation scale.}
1917 \label{fig:realistic-streaming}
1918 \end{figure}
1919 
1920 This test is stronger than applying the streamed identity to a precomputed
1921 factor because the centering, variance scaling, covariate projection, and
1922 rare-variant weights are all performed inside the block loop.
1923 
1924 \subsection{Phenotype-space equivalence}
1925 
1926 For each scenario, 100 Gaussian null residual draws were evaluated in three
1927 ways:
1928 
1929 \begin{enumerate}
1930 \item directly through the factor,
1931 \[
1932 Q=\frac12\|A^\top e\|_2^2;
1933 \]
1934 \item through the sample kernel,
1935 \[
1936 Q=\frac12 e^\top Ke;
1937 \]
1938 \item through the eigencoordinates of $K$.
1939 \end{enumerate}
1940 
1941 The worst relative discrepancy across all draws and all scenarios was
1942 
1943 \[
1944 \boxed{8.32\times10^{-16}}.
1945 \]
1946 
1947 Thus the simulated phenotype path, the streamed kernel, and the generalized
1948 $\chi^2$ spectral representation agree to floating-point precision.
1949 
1950 \subsection{Deep tails on realistic LD spectra}
1951 
1952 For every scenario, tail queries were generated at tilt fractions
1953 0.90, 0.97, and 0.99 of the positive singularity.  K-COS was compared with a
1954 higher-resolution full-spectrum tilted reference.
1955 
1956 The large corrected copy-LD case was then pushed further, to fractions
1957 0.995, 0.997, 0.999, and 0.9995.  The deepest tested reference value was
1958 
1959 \[
1960 \boxed{\log_{10}p=-457.973916}.
1961 \]
1962 
1963 At that point the displayed K-COS and high-resolution full-spectrum values
1964 were identical.  Across the complete realistic-LD battery, the largest
1965 absolute discrepancy was
1966 
1967 \[
1968 \boxed{3.21\times10^{-12}\ \text{decades}}.
1969 \]
1970 
1971 \begin{figure}[H]
1972 \centering
1973 \includegraphics[width=.84\linewidth]{figure_20_realistic_deep_tail_accuracy.png}
1974 \caption{K-COS versus the high-resolution full-spectrum tilted reference on
1975 realistic LD spectra.  Accuracy remains at roughly machine-scale in log
1976 probability through the deepest corrected copy-LD query, near $10^{-458}$.}
1977 \label{fig:realistic-deep-tail}
1978 \end{figure}
1979 
1980 For the $n=1000$, $m=5000$ corrected copy-LD spectrum, selected results are
1981 
1982 \begin{center}
1983 \begin{tabular}{rrrr}
1984 \toprule
1985 tilt fraction & $\log_{10}p$ & K-COS $K$ & absolute error (decades)\\
1986 \midrule
1987 0.9950 & -65.950669 & 2048 & $4.69\times10^{-13}$\\
1988 0.9970 & -95.305034 & 2048 & $3.70\times10^{-14}$\\
1989 0.9990 & -240.602047 & 8192 & $4.94\times10^{-14}$\\
1990 0.9995 & -457.973916 & 16384 & below displayed precision\\
1991 \bottomrule
1992 \end{tabular}
1993 \end{center}
1994 
1995 \subsection{Interpretation and boundary}
1996 
1997 This simulation block does not add a new algebraic theorem.  It supplies a
1998 more demanding validation layer for the theorem chain already proved.
1999 Specifically, it shows that the exact streaming identities and the K-COS
2000 deep-tail machinery remain numerically stable under substantially different
2001 rare-variant LD structures and through the complete raw-genotype preprocessing
2002 path.
2003 
2004 The evidence remains synthetic.  These generators do not reproduce every
2005 population-genetic feature of a sequencing cohort, and this block is not an
2006 external-cohort replication.  It also addresses the Gaussian linear-model null;
2007 finite-sample Bernoulli or logistic calibration remains a separate statistical
2008 problem.
2009 
2010 
2011 \section{Interface to spectral packets and K-COS}
2012 
2013 The matrix packets above operate before eigendecomposition.  They should be
2014 distinguished from the spectral packets used after the eigenvalues
2015 \(\lambda_j\) are known.
2016 
2017 The complete stack is
2018 \[
2019 \boxed{
2020 \begin{aligned}
2021 G
2022 &\to \text{storage blocks}\\
2023 &\to \text{matrix phases}\\
2024 &\to \text{streamed }A,A^\top,K\\
2025 &\to \text{polyphase/rank-adaptive Krylov}\\
2026 &\to \text{small spectral core}\\
2027 &\to \{\lambda_j\}\\
2028 &\to \text{centered spectral packets}\\
2029 &\to \text{deep-tail K-COS}.
2030 \end{aligned}
2031 }
2032 \]
2033 
2034 Storage decomposition, structural decomposition, spectral compression, and
2035 tail compression are independent axes and can be active simultaneously.
2036 
2037 
2038 \begin{figure}[H]
2039 \centering
2040 \includegraphics[width=.86\linewidth]{figure_07_spectral_packets.png}
2041 \caption{Spectral packets form a second, distinct decomposition layer. After the large matrix problem has been reduced to eigenvalues, nearby eigenvalues are represented by packet centers and centered moments for K-COS transform evaluation. These spectral packets should not be confused with the storage blocks used before eigendecomposition.}
2042 \label{fig:spectral-packets}
2043 \end{figure}
2044 
2045 \section{Memory planning law}
2046 
2047 Let \(d\) be the number of bytes per floating scalar, \(c\) the resident
2048 variant count, \(s\) the Krylov width, and \(p\) the covariate rank.
2049 Ignoring library workspace and allocator overhead, a streamed sample-space
2050 action has working memory on the order of
2051 \[
2052 \boxed{
2053 d\,[nc+2ns+cs+np+O(ps+c)].
2054 }
2055 \]
2056 
2057 The full factor costs
2058 \[
2059 dnm.
2060 \]
2061 
2062 For fixed \(n,p,s,c\),
2063 \[
2064 \boxed{
2065 \frac{M_{\rm stream}}{M_{\rm full}}
2066 =
2067 O\!\left(
2068 \frac{c}{m}
2069 +
2070 \frac{s}{m}
2071 +
2072 \frac{p}{m}
2073 \right),
2074 }
2075 \]
2076 so device memory remains bounded as total variant count grows.
2077 
2078 \begin{table}[H]
2079 \centering
2080 \caption{Algebraic float64 memory model for \(n=2000\), \(p=4\),
2081 Krylov width \(s=16\), resident variant block \(c=1024\).}
2082 \begin{tabular}{rrrr}
2083 \toprule
2084 \(m\) & full factor (MiB) & streamed \(KV\) (MiB) & tiled Gram (MiB)\\
2085 \midrule
2086 10,000 & 152.6 & 16.308 & 8.555\\
2087 50,000 & 762.9 & 16.308 & 8.555\\
2088 100,000 & 1,525.9 & 16.308 & 8.555\\
2089 500,000 & 7,629.4 & 16.308 & 8.555\\
2090 1,000,000 & 15,258.8 & 16.308 & 8.555\\
2091 \bottomrule
2092 \end{tabular}
2093 \end{table}
2094 
2095 
2096 \begin{figure}[H]
2097 \centering
2098 \includegraphics[width=.84\linewidth]{figure_03_memory_law.png}
2099 \caption{The memory law. The complete factor grows linearly with total variant count $m$, whereas the streamed kernel action and the two-dimensional tiled construction remain bounded when resident block size, Krylov width, and covariate rank are fixed.}
2100 \label{fig:memory-law}
2101 \end{figure}
2102 
2103 \section{Numerical verification}
2104 
2105 All identities were compared against explicitly formed reference matrices.
2106 
2107 \begin{table}[H]
2108 \centering
2109 \caption{Partition-invariance checks.}
2110 \begin{tabular}{lrrr}
2111 \toprule
2112 partition & blocks & Gram relative error & action relative error\\
2113 \midrule
2114 coarse 4 & 4 & \(2.91\times10^{-16}\) & \(5.99\times10^{-16}\)\\
2115 fine 64 & 32 & \(3.39\times10^{-16}\) & \(6.00\times10^{-16}\)\\
2116 irregular & 9 & \(3.34\times10^{-16}\) & \(5.87\times10^{-16}\)\\
2117 singleton-front & 21 & \(2.43\times10^{-16}\) & \(5.82\times10^{-16}\)\\
2118 permuted, block 73 & 29 & \(3.52\times10^{-16}\) & \(5.89\times10^{-16}\)\\
2119 \bottomrule
2120 \end{tabular}
2121 \end{table}
2122 
2123 A hierarchical coarse-to-fine refinement gave relative discrepancy
2124 \[
2125 3.17\times10^{-16}.
2126 \]
2127 
2128 For an explicitly PSD polyphase operator
2129 \[
2130 H=L^\top L
2131 \]
2132 with block-bidiagonal \(L\), the smallest eigenvalue was approximately
2133 \[
2134 2.86,
2135 \]
2136 and the largest norm observed outside the exact \(k\)-step light cone was
2137 \[
2138 \boxed{0}.
2139 \]
2140 
2141 The QR-core update reproduced the full singular spectrum with relative
2142 \(\ell_2\) discrepancy
2143 \[
2144 5.89\times10^{-16}.
2145 \]
2146 
2147 
2148 \section{A visual map of the decomposition calculus}
2149 
2150 Figures~\ref{fig:residualization-geometry}--\ref{fig:partition-invariance} trace the mathematical object through its successive representations. Residualization first restricts the problem to the null-model orthogonal complement. Exact block-Gram additivity then permits memory-bounded streaming. Partition invariance separates physical storage granularity from the operator itself, while the phase graph constrains propagation through an exact light cone. Phase-local rank removes redundant basis directions; QR-core reuse preserves spectra under small updates; and spectral packets compress the final generalized-$\chi^2$ transform. The figures therefore represent a sequence of exact representation changes, each preserving the semantic object needed by the next layer.
2151 
2152 \section{Solver portfolio}
2153 
2154 \begin{longtable}{p{0.25\linewidth}p{0.30\linewidth}p{0.36\linewidth}}
2155 \toprule
2156 Regime & Preferred method & Reason\\
2157 \midrule
2158 factor fits comfortably & direct factor or sample Gram & lowest implementation overhead\\
2159 factor exceeds device memory & streamed \(AV\), \(A^\top U\), or \(KV\) & exact memory-bounded execution\\
2160 single block still too large & 2D sample/variant tiling & exact resident-tile bound\\
2161 low numerical rank & block Golub--Kahan / range core & solve at rank scale\\
2162 sparse exact phase graph & polyphase frontier Krylov & skip unreachable phases\\
2163 nonuniform phase novelty & phase-local rank allocation & dimension-minimal exact width\\
2164 low-rank phase interface & compressed message factors & reduce communication/storage\\
2165 overlapping sets & QR-core update & exact reuse\\
2166 full spectrum known & K-COS spectral packets & deep-tail transform compression\\
2167 \bottomrule
2168 \end{longtable}
2169 
2170 
2171 \section{Relation to prior art and candidate novelty}
2172 
2173 The individual numerical primitives used in this manuscript are not new in
2174 isolation.  Blocked matrix multiplication, out-of-core linear algebra,
2175 Lanczos and Golub--Kahan iteration, QR factorization, and low-rank operator
2176 application are classical tools.  Randomized matrix approximation also has a
2177 well-developed pass-efficient and streaming literature; in particular,
2178 Halko, Martinsson, and Tropp describe algorithms intended for matrices too
2179 large to fit in fast memory and emphasize constant-pass or single-pass
2180 approaches \cite{Halko2011}.  In the SKAT literature, fastSKAT specifically
2181 uses randomized projection and Lanczos-type ideas to avoid a complete dense
2182 eigendecomposition, computing leading eigenvalues and approximating the
2183 remaining spectral contribution \cite{Lumley2018}.  The original SKAT
2184 framework and its weighted quadratic-form null distribution are due to Wu
2185 et al. \cite{Wu2011}.
2186 
2187 Accordingly, the contribution claimed here is \emph{not} the invention of
2188 blocking, streaming, matrix-free multiplication, or Krylov iteration.  The
2189 candidate contribution is the exact composition of four independent
2190 representations:
2191 \[
2192 \boxed{
2193 \mathcal D
2194 =
2195 \mathcal D_{\rm storage}
2196 \circ
2197 \mathcal D_{\rm phase}
2198 \circ
2199 \mathcal D_{\rm rank}
2200 \circ
2201 \mathcal D_{\rm spectrum}.
2202 }
2203 \]
2204 The first controls device residency, the second controls structural
2205 propagation, the third controls spectral dimension, and the fourth controls
2206 deep-tail transform cost.
2207 
2208 \subsection{Storage-phase commutation}
2209 
2210 Fix a structural phase partition
2211 \[
2212 \{1,\ldots,m\}=\mathcal P_1\sqcup\cdots\sqcup\mathcal P_H
2213 \]
2214 and let each phase be subdivided into arbitrary storage blocks
2215 \[
2216 \mathcal P_h
2217 =
2218 J_{h1}\sqcup\cdots\sqcup J_{hB_h}.
2219 \]
2220 Let $S$ denote this physical storage subdivision.  The phase projectors
2221 $P_h$ are defined by the structural phases, not by $S$.
2222 
2223 \begin{theorem}[Storage-phase commutation theorem]
2224 Let
2225 \[
2226 A=M_XGW^{1/2},
2227 \qquad H=A^\top A,
2228 \]
2229 and let $S$ and $S'$ be any two exact storage refinements that respect the
2230 same structural phase partition.  Then:
2231 \begin{enumerate}[label=(\roman*)]
2232 \item the exact operator $H$ is identical under $S$ and $S'$;
2233 \item every phase coupling $P_hHP_g$ is identical under $S$ and $S'$;
2234 \item the exact phase graph and every $k$-step polyphase light cone are
2235 identical under $S$ and $S'$;
2236 \item any exact matrix-free Krylov method applied through streamed products
2237 produces the same Krylov subspaces in exact arithmetic under $S$ and $S'$.
2238 \end{enumerate}
2239 \end{theorem}
2240 
2241 \begin{proof}
2242 Storage refinement changes only the order and grouping of finite sums.
2243 By partition refinement invariance, the streamed products reconstruct the
2244 same $A$, $A^\top$, and therefore the same $H=A^\top A$.  Since each
2245 $P_h$ is fixed by the structural phase partition, every block
2246 $P_hHP_g$ is unchanged.  The phase graph is defined only by these blocks, so
2247 its path structure and light cones are unchanged.  Finally, an exact
2248 matrix-free Krylov method depends only on repeated applications of the same
2249 operator to the same starting block; hence its generated subspaces are
2250 identical.
2251 \end{proof}
2252 
2253 \begin{corollary}[Adaptive resident-block scheduler]
2254 An implementation may change resident block size at runtime---including
2255 recursively splitting a phase block after an allocation failure---without
2256 changing the exact phase graph or exact Krylov trajectory, provided all
2257 required storage blocks are processed exactly.
2258 \end{corollary}
2259 
2260 This separates \emph{mathematical locality} from \emph{physical storage
2261 granularity}.  A phase may remain fixed as an LD window or annotation class
2262 while its physical representation is subdivided according to free GPU
2263 memory, scalar precision, transfer granularity, or stream count.
2264 
2265 \subsection{Frontier-to-I/O theorem}
2266 
2267 \begin{theorem}[Exact frontier-to-I/O reduction]
2268 Suppose the exact phase graph of $H$ is used and the current Krylov block is
2269 supported on a phase set $F_0$.  At depth $k$, let $F_k$ be the graph
2270 $k$-neighborhood of $F_0$.  Then no storage block lying wholly in a phase
2271 outside $F_k$ can affect the $k$-step Krylov vector.  Such a block need not be
2272 loaded during that step.
2273 \end{theorem}
2274 
2275 \begin{proof}
2276 For every phase $h\notin F_k$ and every starting phase $g\in F_0$, the
2277 polyphase light-cone theorem gives $P_hH^kP_g=0$.  Storage-phase commutation
2278 allows each active phase operator to be evaluated through arbitrary exact
2279 storage blocks.  Therefore blocks outside the active structural frontier
2280 contribute identically zero to the requested iterate and may be omitted.
2281 \end{proof}
2282 
2283 The theorem turns an algebraic support statement into an execution law:
2284 \[
2285 \boxed{
2286 \text{Krylov depth}
2287 \longrightarrow
2288 \text{reachable phases}
2289 \longrightarrow
2290 \text{storage blocks that must be transferred}.
2291 }
2292 \]
2293 This is stronger than ordinary chunking, because the scheduler is constrained
2294 by the exact operator graph rather than solely by available memory.
2295 
2296 \subsection{Comparison with fastSKAT}
2297 
2298 fastSKAT attacks the same broad spectral bottleneck but follows a different
2299 approximation strategy: it computes leading eigenvalues using randomized or
2300 Lanczos-type methods and approximates the remainder of the quadratic-form
2301 spectrum \cite{Lumley2018}.  The present architecture instead separates
2302 memory blocking from structural propagation and can continue, when desired,
2303 to exact range capture and an explicit small spectral core before handing the
2304 resulting eigenvalues to a separate deep-tail engine.  Thus the strongest
2305 novelty claim is not ``matrix-free SKAT,'' which is prior art, but the
2306 storage-invariant polyphase scheduler and its composition with local-rank,
2307 core-reuse, and tail-packet layers.
2308 
2309 \subsection{Novelty statement and boundary}
2310 
2311 A targeted literature review conducted for this manuscript identified prior
2312 work on SKAT, fastSKAT, out-of-core randomized matrix decompositions, and
2313 classical Krylov solvers, but did not identify a SKAT solver that combines all
2314 of the following in one exact architecture:
2315 \begin{enumerate}
2316 \item residualization-preserving streamed factor actions;
2317 \item storage refinement invariant under runtime block subdivision;
2318 \item an exact structural phase graph with finite-step light cones;
2319 \item phase-local minimal Krylov width;
2320 \item low-rank inter-phase messages and exact QR-core update reuse;
2321 \item an explicit transition from the small spectral core to a separately
2322 certified deep-tail spectral-packet calculation.
2323 \end{enumerate}
2324 This absence is evidence for a candidate systems/methodological contribution,
2325 not a proof of global novelty.  The appropriate publication claim is that the
2326 \emph{composition and its exact invariance laws appear to be novel}, while
2327 its numerical primitives are classical.
2328 
2329 \section{Claim boundary}
2330 
2331 The manuscript deliberately uses a narrow novelty boundary. Blocking, tiled matrix multiplication, out-of-core linear algebra, randomized low-rank approximation, Golub--Kahan/Lanczos iteration, QR, and low-rank factor application are inherited techniques \cite{Halko2011,Lumley2018}. The contribution should not be stated as the invention of blocking or matrix-free SKAT.
2332 
2333 The candidate Komorebi contribution is the resolved synthesis
2334 \[
2335 \boxed{\begin{aligned}
2336 &\text{exact streamed factor actions}
2337 +\text{ polyphase light cones}
2338 +\text{ minimal phase rank}\\
2339 &+\text{ low-rank messages}
2340 +\text{ exact update cores}
2341 +\text{ deep-tail spectral packets}.
2342 \end{aligned}}
2343 \]
2344 
2345 The conceptual claim is that storage decomposition, structural phase
2346 decomposition, spectral compression, and tail compression are mathematically
2347 distinct layers and can be composed without changing the exact operator.
2348 
2349 \section{Current proof and validation ledger}
2350 
2351 \small
2352 \begin{tabular}{p{0.46\linewidth}p{0.18\linewidth}p{0.26\linewidth}}
2353 \toprule
2354 Item & Status & Evidence\\
2355 \midrule
2356 block Gram additivity & proved & direct matrix identity\\
2357 streamed \(KV\) & proved & direct matrix identity\\
2358 streamed \(AV,A^\top U\) & proved & factor partition\\
2359 2D tiling & proved & row-block expansion\\
2360 partition refinement invariance & proved & disjoint-union argument\\
2361 variant permutation invariance & proved & simultaneous coordinate permutation\\
2362 storage-phase commutation & proved & operator invariance under storage refinement\\
2363 frontier-to-I/O reduction & proved & light-cone support plus storage commutation\\
2364 polyphase light cone & proved & projector-path expansion\\
2365 minimal phase width & proved & dimension argument\\
2366 QR-core update & proved & orthonormal left factor\\
2367 range-core spectrum & proved & singular-value invariance\\
2368 machine-precision implementation & tested & explicit-reference comparisons\\
2369 GPU wall-clock crossover & open & requires production CUDA implementation\\
2370 I/O-optimal scheduling & open & hardware/backing-store dependent\\
2371 approximate phase graph error bound & proved/tested & weighted transfer matrix and thresholded-power certificate\\
2372 budget-optimal equal-cost scheduler & proved/tested & additive edge-risk theorem plus exhaustive audit\\
2373 unequal-cost scheduler & formulated/tested & certified knapsack heuristic with fractional upper bound\\
2374 double-buffer pipeline model & proved/model only & ideal asynchronous-overlap proposition\\
2375 trace-aware residual MGF envelope & proved/tested & exact residual decomposition plus trace--Frobenius--cap bound\\
2376 residual-rank moment bracket & proved/tested & trace/cap and trace-square rank inequalities\\
2377 streamed residual completion & proved/tested & orthogonal residual factor and incremental completion core\\
2378 \bottomrule
2379 \end{tabular}
2380 \normalsize
2381 
2382 \section{Next formal targets}
2383 
2384 The next blocks should address:
2385 
2386 \begin{enumerate}
2387 \item integrate the hybrid Krylov/moment/completion stopping rule into the
2388 production K-COS API;
2389 \item implement rank-revealing residual completion with fused GPU block QR or
2390 TSQR rather than the reference SVD path;
2391 \item derive randomized or block-Krylov convergence bounds under streamed,
2392 phase-pruned actions;
2393 \item benchmark the complete scheduler on native GPU hardware, including
2394 PCIe overlap, allocator workspace, and block-size crossover.
2395 \end{enumerate}
2396 
2397 \appendix
2398 
2399 \section{Reference pseudocode}
2400 
2401 \begin{lstlisting}[language=Python]
2402 def stream_KV(block_source, weights, Qx, V):
2403     U = V - Qx @ (Qx.T @ V)
2404     Y = zeros_like(U)
2405     for G_block, w_block in block_source:
2406         T = G_block.T @ U
2407         Y += G_block @ (w_block[:, None] * T)
2408     return Y - Qx @ (Qx.T @ Y)
2409 \end{lstlisting}
2410 
2411 \section{Running manuscript protocol}
2412 
2413 This file is intended to remain cumulative.  New theorem blocks should be
2414 added with one of three labels:
2415 
2416 \begin{itemize}
2417 \item \textbf{proved}: a complete mathematical argument is present;
2418 \item \textbf{validated}: a numerical experiment supports the statement but
2419 does not prove it;
2420 \item \textbf{open}: the obstruction or missing assumption is recorded.
2421 \end{itemize}
2422 
2423 Negative timings, failed approximations, and counterexamples should be retained
2424 in the running ledger rather than removed from later versions.
2425 
2426 
2427 \begin{thebibliography}{9}
2428 
2429 \bibitem{Wu2011}
2430 M. C. Wu, S. Lee, T. Cai, Y. Li, M. Boehnke, and X. Lin,
2431 ``Rare-Variant Association Testing for Sequencing Data with the Sequence
2432 Kernel Association Test,''
2433 \emph{American Journal of Human Genetics}, 89(1):82--93, 2011.
2434 DOI: 10.1016/j.ajhg.2011.05.029.
2435 
2436 \bibitem{Halko2011}
2437 N. Halko, P. G. Martinsson, and J. A. Tropp,
2438 ``Finding Structure with Randomness: Probabilistic Algorithms for Constructing
2439 Approximate Matrix Decompositions,''
2440 \emph{SIAM Review}, 53(2):217--288, 2011.
2441 DOI: 10.1137/090771806.
2442 
2443 \bibitem{Lumley2018}
2444 T. Lumley, J. Brody, G. Peloso, A. Morrison, and K. Rice,
2445 ``FastSKAT: Sequence Kernel Association Tests for Very Large Sets of Markers,''
2446 \emph{Genetic Epidemiology}, 2018.
2447 PMCID: PMC6129408.
2448 
2449 \end{thebibliography}
2450 
2451 \end{document}

Back to source controls

Back to Streaming Polyphase Factor Decomposition / SKAT