Realization Calculus

Download original file

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

Read the typeset paper

1 \section{Introduction}\label{main:section:001}\label{sec:introduction}
2 
3 A matrix can look dense while its action is not. Independent diagonal blocks, sparse transfers, low rank couplings, tensor factors, and transform coordinates can remove most of the arithmetic and storage before a kernel is launched. The gain is mathematical first. A device can exploit only the structure that the representation exposes.
4 
5 This article begins with a finite direct sum
6 \[
7 X=\bigoplus_{\alpha\in I}X_\alpha
8 \]
9 and a complete family of coordinate projectors
10 \[
11 P_\alpha P_\beta=\delta_{\alpha\beta}P_\alpha,
12 \qquad
13 \sum_{\alpha\in I}P_\alpha=I.
14 \]
15 An operator has packet blocks
16 \[
17 A_{\beta\alpha}=P_\beta AP_\alpha.
18 \]
19 The active addresses form a directed graph. Block dimensions, block types, approximation tolerances, and output requirements turn this graph into a computational intermediate representation.
20 
21 The underlying linear algebra is classical. Direct and iterative methods are treated in \cite{GolubVanLoan2013,Saad2003,Davis2006}. Low rank approximation and updates are governed by singular value decomposition and the Sherman--Morrison--Woodbury identity \cite{EckartYoung1936,Mirsky1960,Hager1989,BunchNielsen1978}. Floating point analysis follows the standard model in \cite{Higham2002,IEEE7542019}. Randomized discovery of low rank structure is developed in \cite{HalkoMartinssonTropp2011,MartinssonTropp2020}. Communication and attainable performance are shaped by data movement and arithmetic intensity \cite{HongKung1981,WilliamsWatermanPatterson2009,BallardDemmelHoltzSchwartz2011,FrigoLeisersonProkopRamachandran1999}.
22 
23 The first task is semantic. A compressed representation must specify the operator that it denotes. It must also specify what an algorithm returns. A product represented by block and low rank factors is not the same output as an explicit matrix. A timing comparison is valid only after the output contract and timing scope agree.
24 
25 The second task is algebraic. Block diagonal plus low rank operators are closed under multiplication before recompression. Their unrecompressed rank grows additively. Woodbury inversion reduces a large inverse to packet solves and one small dense inverse. The low rank factors can then be recompressed by orthogonalizing them and taking the singular value decomposition of a small core. This produces the best truncated approximation in spectral and Frobenius norm.
26 
27 The third task is combinatorial. Sparse transfer graphs predict possible support under powers, polynomials, and Schur elimination. A directed acyclic graph gives a finite path expansion. A block triangular solve can be scheduled by topological level. Graph regularity matters because equal shapes can be batched, while irregular edges require indexing and accumulation.
28 
29 The fourth task is numerical. Packet perturbations can be certified by their block sensitivities. A continuous model assigns tolerances through a convex calculation. Upper bounds on available tolerances produce a clipped allocation. A finite set of data types gives a multiple choice budget problem that can be solved exactly after the errors are quantized.
30 
31 The central computational result concerns compilation. Suppose a finite rooted tree records candidate refinements. Each node may be kept as one representation or split into its children. If cost and error are additive across the selected leaves, the optimal pruning is found exactly by dynamic programming. A quantized error budget has a second exact dynamic program. The decision is stable when the score gap exceeds the perturbation in the local models. A change inside one subtree requires recomputation only along the ancestor chain. These results prove a restricted packet compiler. They do not solve unrestricted discovery or arbitrary device scheduling.
32 
33 The measured part uses a sanitized corpus collected on one recorded CPU and one NVIDIA GeForce RTX 3070 Ti Laptop GPU. The preserved corpus contains 154 recorded rows from sixteen suites. The timings are inherited observations. No GPU timing is rerun for this edition. Regular equal block operators and block diagonal plus low rank operators show large gains. Recursive HODLR, Strassen multiplication, iterative refinement, a small Walsh--Hadamard transform, and an attempted two stream overlap provide negative controls. These rows are retained because structure does not guarantee speed.
34 
35 The article separates four sources of gain:
36 \[
37 \text{representation},
38 \qquad
39 \text{device},
40 \qquad
41 \text{kernel},
42 \qquad
43 \text{amortization}.
44 \]
45 A single speedup number can mix all four. The benchmark protocol therefore requires case identity, output contract, timing scope, source precision, correctness reference, discovery time, compilation time, and reuse count. The archived rows do not satisfy every field of this ideal contract: source-value hashes and explicit correctness-reference fields are absent, while setup, discovery, and reuse metadata are sparse. Recorded ratios are consequently descriptive comparisons under their stated scopes, rather than fair matches with complete provenance.
46 
47 \begin{figure}[t]
48 \centering
49 \includegraphics[width=\textwidth]{01_packet_compilation_architecture.pdf}
50 \caption{The representation and the output contract are fixed before the device schedule is chosen.}
51 \label{fig:architecture}
52 \end{figure}
53 
54 The mathematical statements are finite. The empirical statements apply only to the recorded environment and generated cases. The article does not claim a new generic matrix multiplication exponent. Fast matrix multiplication and its stability form a separate theory \cite{Strassen1969,DemmelHoltzSchwartz2007}. It also does not claim that every structured matrix has a profitable GPU realization.
55 
56 \section{Packet operator representations}\label{main:section:002}\label{sec:ir}
57 
58 \subsection{Exact semantics}\label{main:subsection:001}
59 
60 \begin{definition}[Packet layout]\label{main:definition:001}\label{def:packet_layout}
61 A packet layout is a finite family of finite dimensional Hilbert spaces
62 \[
63 \Pcal=(X_\alpha)_{\alpha\in I}
64 \]
65 with direct sum
66 \[
67 X=\bigoplus_{\alpha\in I}X_\alpha.
68 \]
69 Write
70 \[
71 n_\alpha=\dim X_\alpha,
72 \qquad
73 n=\sum_{\alpha\in I}n_\alpha.
74 \]
75 The coordinate inclusion and projection are $\iota_\alpha:X_\alpha\to X$ and $\pi_\alpha:X\to X_\alpha$. The ambient coordinate projector is $P_\alpha=\iota_\alpha\pi_\alpha$.
76 \end{definition}
77 
78 \begin{definition}[Packet operator representation]\label{main:definition:002}\label{def:packet_ir}
79 A finite packet operator representation consists of
80 \[
81 \IR=(\Pcal,E,(B_{\beta\alpha})_{(\beta,\alpha)\in E}),
82 \]
83 where
84 \[
85 E\subseteq I\times I
86 \]
87 and
88 \[
89 B_{\beta\alpha}:X_\alpha\longrightarrow X_\beta.
90 \]
91 Local blocks are identified with their ambient extensions by inclusion and projection. In particular, the notation $P_\beta AP_\alpha$ denotes either the ambient block or its restriction $X_\alpha\to X_\beta$ according to the displayed domain, and $\iota_\beta B_{\beta\alpha}\pi_\alpha$ is the corresponding ambient extension. Its evaluation is
92 \[
93 \operatorname{eval}(\IR)
94 =
95 \sum_{(\beta,\alpha)\in E}
96 \iota_\beta B_{\beta\alpha}\pi_\alpha,
97 \]
98 where \(\iota_\beta:X_\beta\to X\) is the coordinate inclusion.
99 \end{definition}
100 
101 \begin{theorem}[Exact evaluation]\label{main:theorem:001}\label{thm:exact_evaluation}
102 Let \(A\in\End(X)\). Set
103 \[
104 B_{\beta\alpha}=\pi_\beta A\iota_\alpha
105 \]
106 and let \(E\) contain exactly the nonzero blocks. Then
107 \[
108 \boxed{
109 A=
110 \sum_{(\beta,\alpha)\in E}
111 P_\beta AP_\alpha.
112 }
113 \]
114 For every matrix of right hand sides \(Z\),
115 \[
116 AZ
117 =
118 \sum_{(\beta,\alpha)\in E}
119 \iota_\beta B_{\beta\alpha}\pi_\alpha Z.
120 \]
121 \end{theorem}
122 
123 \begin{proof}\phantomsection\label{main:proof:001}
124 Insert the identity on both sides of \(A\):
125 \[
126 A
127 =
128 \left(\sum_\beta P_\beta\right)
129 A
130 \left(\sum_\alpha P_\alpha\right)
131 =
132 \sum_{\alpha,\beta}P_\beta AP_\alpha.
133 \]
134 The zero blocks can be removed. Applying the identity to \(Z\) gives the second formula.
135 \phantomsection\label{main:proof:001:end}\end{proof}
136 
137 \begin{definition}[Support graph]\label{main:definition:003}\label{def:support_graph}
138 The support graph \(G_A\) has vertex set \(I\). It has a directed edge
139 \[
140 \alpha\longrightarrow\beta
141 \]
142 when \(P_\beta AP_\alpha\neq0\).
143 \end{definition}
144 
145 \begin{figure}[t]
146 \centering
147 \includegraphics[width=.88\textwidth]{02_packet_ir_graph.pdf}
148 \caption{A packet intermediate representation stores active block addresses together with dimensions, block types, tolerances, and scheduling metadata.}
149 \label{fig:packet_ir}
150 \end{figure}
151 
152 \begin{theorem}[Block composition]\label{main:theorem:002}\label{thm:block_composition}
153 Let \(A,B\in\End(X)\). Then
154 \[
155 \boxed{
156 (AB)_{\gamma\alpha}
157 =
158 \sum_{\beta\in I}
159 A_{\gamma\beta}B_{\beta\alpha}.
160 }
161 \]
162 Consequently,
163 \[
164 \supp(AB)
165 \subseteq
166 \supp(A)\circ\supp(B),
167 \]
168 where relational composition is defined by
169 \[
170 (\gamma,\alpha)\in E\circ F
171 \]
172 when some \(\beta\) satisfies
173 \[
174 (\gamma,\beta)\in E,
175 \qquad
176 (\beta,\alpha)\in F.
177 \]
178 \end{theorem}
179 
180 \begin{proof}\phantomsection\label{main:proof:002}
181 Insert the packet identity between \(A\) and \(B\):
182 \[
183 P_\gamma ABP_\alpha
184 =
185 \sum_\beta
186 P_\gamma AP_\beta BP_\alpha.
187 \]
188 A summand can be nonzero only when both block addresses are active.
189 \phantomsection\label{main:proof:002:end}\end{proof}
190 
191 \begin{corollary}[Adjoint and tensor semantics]\label{main:corollary:001}\label{cor:adjoint_tensor}
192 The adjoint representation has blocks
193 \[
194 (A^*)_{\alpha\beta}=A_{\beta\alpha}^*.
195 \]
196 If \(A\) and \(B\) are packet resolved on \(X\) and \(Y\), then
197 \[
198 (A\otimes B)_{(\beta,\delta),(\alpha,\gamma)}
199 =
200 A_{\beta\alpha}\otimes B_{\delta\gamma}.
201 \]
202 \end{corollary}
203 
204 \begin{proof}\phantomsection\label{main:proof:003}
205 Both identities follow by inserting the corresponding coordinate projectors and using
206 \[
207 (P_\beta AP_\alpha)^*=P_\alpha A^*P_\beta.
208 \]
209 The tensor identity follows from
210 \[
211 (P_\beta\otimes Q_\delta)
212 (A\otimes B)
213 (P_\alpha\otimes Q_\gamma)
214 =
215 A_{\beta\alpha}\otimes B_{\delta\gamma}.
216 \]
217 \phantomsection\label{main:proof:003:end}\end{proof}
218 
219 \subsection{Output contracts}\label{main:subsection:002}
220 
221 \begin{definition}[Output contract]\label{main:definition:004}\label{def:output_contract}
222 An output contract specifies the mathematical object returned by an algorithm. The contract may require an explicit vector, an explicit matrix, block factors, low rank factors, spectral factors, or an implicit operator with a declared application rule.
223 \end{definition}
224 
225 \begin{definition}[Timing scope]\label{main:definition:005}\label{def:timing_scope}
226 A timing scope states which costs are included. Possible costs include structure discovery, compilation, host transfer, device transfer, kernel execution, graph capture, graph replay, factor reconstruction, and explicit materialization.
227 \end{definition}
228 
229 \begin{proposition}[Contract conversion]\label{main:proposition:001}\label{prop:contract_conversion}
230 Suppose one algorithm returns a compressed representation \(R\) and another returns an explicit object \(Y\). A direct timing comparison is valid only if
231 \[
232 \operatorname{decode}(R)=Y
233 \]
234 and the time required by \(\operatorname{decode}\) is included, or if the downstream task consumes \(R\) directly.
235 \end{proposition}
236 
237 \begin{proof}\phantomsection\label{main:proof:004}
238 Without decoding, the two procedures return different mathematical objects. A timing ratio between different outputs does not compare implementations of the same map. Including the decoding cost gives two procedures with the same output. If the next operation consumes \(R\), decoding is not part of the declared computation.
239 \phantomsection\label{main:proof:004:end}\end{proof}
240 
241 \begin{remark}\label{main:remark:001}\label{rem:semantic_boundary}
242 A compressed product and an explicit matrix product may both be exact. They still have different contracts. This distinction becomes decisive when materialization costs \(n^2\) entries.
243 \end{remark}
244 
245 \section{Arithmetic, storage, and data movement}\label{main:section:003}\label{sec:cost}
246 
247 We count one scalar multiplication and one scalar addition as two modeled floating point operations. These counts describe a declared event model for a representation. They do not predict wall time without a memory and launch model. A multiplication with accumulation charges two scalar operations even for the first coefficient. Thin corrections accumulate into an already computed diagonal or sparse result. Index metadata, temporary workspace, launch overhead, discovery, transfer, and materialization are accounted for separately.
248 
249 \begin{theorem}[Dense and packet application counts]\label{main:theorem:003}\label{thm:packet_counts}
250 A dense \(n\) by \(n\) matrix applied to \(q\) right hand sides has modeled cost
251 \[
252 \boxed{
253 \flop_{\rm dense}=2n^2q.
254 }
255 \]
256 A packet representation with active edge set \(E\) has modeled cost
257 \[
258 \boxed{
259 \flop_{\rm packet}
260 =
261 2q
262 \sum_{(\beta,\alpha)\in E}
263 n_\beta n_\alpha.
264 }
265 \]
266 It stores
267 \[
268 \boxed{
269 M_{\rm packet}
270 =
271 \sum_{(\beta,\alpha)\in E}
272 n_\beta n_\alpha
273 }
274 \]
275 scalar entries.
276 \end{theorem}
277 
278 \begin{proof}\phantomsection\label{main:proof:005}
279 The active block \(B_{\beta\alpha}\) is an \(n_\beta\) by \(n_\alpha\) matrix. Applying it to an \(n_\alpha\) by \(q\) input costs \(2n_\beta n_\alpha q\) modeled operations and stores \(n_\beta n_\alpha\) scalars. Summation over the active addresses gives both formulas.
280 \phantomsection\label{main:proof:005:end}\end{proof}
281 
282 \begin{corollary}[Equal diagonal packets]\label{main:corollary:002}\label{cor:equal_blocks}
283 Let $p\ge1$ and $b\ge1$ be integers and $n=pb$. Suppose $n$ is divided into $p$ equal diagonal packets of width $b$. Then
284 \[
285 \flop_{\rm block}=2qnb,
286 \qquad
287 M_{\rm block}=nb.
288 \]
289 The following factor describes the dense-to-block storage reduction and, when $q\ge1$, the arithmetic reduction:
290 \[
291 \boxed{
292 \frac{n}{b}.
293 }
294 \]
295 \end{corollary}
296 
297 \begin{proof}\phantomsection\label{main:proof:006}
298 There are \(n/b\) active diagonal blocks. Each contains \(b^2\) entries.
299 \phantomsection\label{main:proof:006:end}\end{proof}
300 
301 \begin{definition}[Block diagonal plus low rank operator]\label{main:definition:006}\label{def:bdlr}
302 A block diagonal plus low rank operator has the form
303 \[
304 A=D+UV^*,
305 \]
306 where
307 \[
308 D=\bigoplus_{\alpha\in I}D_\alpha,
309 \qquad
310 U,V\in\CC^{n\times r}.
311 \]
312 The abbreviation BDLR will be used below.
313 \end{definition}
314 
315 \begin{theorem}[BDLR application and storage]\label{main:theorem:004}\label{thm:bdlr_counts}
316 For \(q\) right hand sides,
317 \[
318 AZ=DZ+U(V^*Z)
319 \]
320 has modeled cost
321 \[
322 \boxed{
323 2q
324 \left(
325 \sum_\alpha n_\alpha^2+2nr
326 \right).
327 }
328 \]
329 The representation stores
330 \[
331 \boxed{
332 \sum_\alpha n_\alpha^2+2nr
333 }
334 \]
335 scalars. For $n=pb$, with integers $p,b\ge1$, the following factor describes the ideal dense storage reduction and, when $q\ge1$, the arithmetic reduction:
336 \[
337 \boxed{
338 \frac{n}{b+2r}.
339 }
340 \]
341 \end{theorem}
342 
343 \begin{proof}\phantomsection\label{main:proof:007}
344 The block diagonal product costs
345 \[
346 2q\sum_\alpha n_\alpha^2.
347 \]
348 The products \(V^*Z\) and \(U(V^*Z)\) cost \(2nrq\) each. Storage consists of the diagonal blocks and the two low rank factors. For equal packets,
349 \[
350 \sum_\alpha n_\alpha^2=nb.
351 \]
352 \phantomsection\label{main:proof:007:end}\end{proof}
353 
354 \begin{figure}[t]
355 \centering
356 \includegraphics[width=.76\textwidth]{03_ideal_arithmetic_reduction.pdf}
357 \caption{The ideal reduction grows with dimension only when packet width and transfer rank remain controlled.}
358 \label{fig:ideal_reduction}
359 \end{figure}
360 
361 \begin{definition}[One pass streaming model]\label{main:definition:007}\label{def:streaming_model}
362 Fix a scalar size \(s\) bytes. For every active block, the one pass streaming model reads the block once, reads its source packet once, and reads and writes its target packet once. It does not assume cache reuse across edges.
363 \end{definition}
364 
365 \begin{proposition}[Streaming byte count]\label{main:proposition:002}\label{prop:streaming_bytes}
366 Under \cref{def:streaming_model}, applying a packet operator to \(q\) right hand sides moves
367 \[
368 \boxed{
369 B_{\rm stream}
370 =
371 s
372 \sum_{(\beta,\alpha)\in E}
373 \left(
374  n_\beta n_\alpha
375  +n_\alpha q
376  +2n_\beta q
377 \right)
378 }
379 \]
380 bytes.
381 \end{proposition}
382 
383 \begin{proof}\phantomsection\label{main:proof:008}
384 The three terms count the matrix block, the source values, and the read plus write traffic for the target accumulation. The model counts these quantities independently for every active edge.
385 \phantomsection\label{main:proof:008:end}\end{proof}
386 
387 \begin{definition}[Arithmetic intensity]\label{main:definition:008}\label{def:intensity}
388 When the modeled byte traffic is strictly positive, the modeled arithmetic intensity is
389 \[
390 \Ical
391 =
392 \frac{\flop}{\bytes}.
393 \]
394 \end{definition}
395 
396 \begin{theorem}[Roofline lower bound]\label{main:theorem:005}\label{thm:roofline}
397 Suppose a device sustains at most $P>0$ modeled operations per second and at most $W>0$ bytes per second. A computation that requires \(F\) operations and \(B\) bytes satisfies
398 \[
399 \boxed{
400 T\geq
401 \max\left\{
402 \frac{F}{P},
403 \frac{B}{W}
404 \right\}.
405 }
406 \]
407 For positive elapsed time and positive byte traffic, its work rate is at most
408 \[
409 \min\{P,W\Ical\}.
410 \]
411 \end{theorem}
412 
413 \begin{proof}\phantomsection\label{main:proof:009}
414 The operation capacity gives \(T\geq F/P\). The bandwidth capacity gives \(T\geq B/W\). Both restrictions hold at once.
415 \phantomsection\label{main:proof:009:end}\end{proof}
416 
417 The roofline model is a guide rather than an exact execution law \cite{WilliamsWatermanPatterson2009}. Packet structure changes both \(F\) and \(B\). An irregular representation may remove arithmetic while lowering intensity enough to become bandwidth or launch limited.
418 
419 \begin{figure}[t]
420 \centering
421 \includegraphics[width=.76\textwidth]{04_roofline_schematic.pdf}
422 \caption{A representation can move the same mathematical problem between memory limited and compute limited regimes.}
423 \label{fig:roofline}
424 \end{figure}
425 
426 \begin{corollary}[Capacity threshold]\label{main:corollary:003}\label{cor:capacity}
427 Let a device have available memory $M\ge0$ bytes and scalar size $s>0$, and let $n\ge0$. The matrix representation alone can fit only when
428 \[
429 n\leq\sqrt{M/s}.
430 \]
431 The equal packet BDLR representation alone fits when
432 \[
433 n(b+2r)s\leq M.
434 \]
435 \end{corollary}
436 
437 \begin{proof}\phantomsection\label{main:proof:010}
438 Substitute the dense and BDLR storage counts into the capacity inequality.
439 \phantomsection\label{main:proof:010:end}\end{proof}
440 
441 \begin{proposition}[Materialization cost]\label{main:proposition:003}\label{prop:materialization_cost}
442 Materializing \(D+UV^*\) as an explicit matrix requires storage for \(n^2\) scalars and modeled arithmetic
443 \[
444 2n^2r+O(n^2)
445 \]
446 when the low rank correction is formed by dense multiplication.
447 \end{proposition}
448 
449 \begin{proof}\phantomsection\label{main:proof:011}
450 The product \(UV^*\) multiplies an \(n\) by \(r\) matrix by an \(r\) by \(n\) matrix. Adding the diagonal part changes only the lower order term.
451 \phantomsection\label{main:proof:011:end}\end{proof}
452 
453 \section{Structured packet primitives}\label{main:section:004}\label{sec:primitives}
454 
455 The packet representation can store different exact local operator types. This section records the finite formulas used by several measured suites. The results concern representation and arithmetic. They do not assert that one implementation is fastest at every size.
456 
457 \begin{proposition}[Equal block banded count]\label{main:proposition:004}\label{prop:block_banded_count}
458 Let $p,b\ge1$ be integers, $n=pb$, and let $p=n/b$ equal packets have width $b$. Suppose a block is active precisely when
459 \[
460 \abs{i-j}\leq w,
461 \qquad
462 0\leq w\leq p-1.
463 \]
464 The number of active blocks is
465 \[
466 \boxed{
467 N_{p,w}
468 =
469 p(2w+1)-w(w+1).
470 }
471 \]
472 For \(q\) right hand sides, the modeled application cost and scalar storage are
473 \[
474 \boxed{
475 2qb^2N_{p,w}
476 \quad\text{and}\quad
477 b^2N_{p,w}.
478 }
479 \]
480 \end{proposition}
481 
482 \begin{proof}\phantomsection\label{main:proof:012}
483 The diagonal contributes \(p\) blocks. At offset \(k\), where \(1\leq k\leq w\), the upper and lower block diagonals contribute \(2(p-k)\) blocks. Therefore
484 \[
485 p+2\sum_{k=1}^{w}(p-k)
486 =
487 p(2w+1)-w(w+1).
488 \]
489 The cost and storage follow from \cref{thm:packet_counts}.
490 \phantomsection\label{main:proof:012:end}\end{proof}
491 
492 \begin{proposition}[Sparse plus low rank application]\label{main:proposition:005}\label{prop:sparse_low_rank}
493 Let
494 \[
495 A=S+UV^*,
496 \]
497 where \(S\) has \(z\) stored scalar nonzeros and \(U,V\in\CC^{n\times r}\). In the model that charges two operations per stored sparse entry and right hand side, applying \(A\) to \(q\) right hand sides costs
498 \[
499 \boxed{
500 2q(z+2nr)
501 }
502 \]
503 modeled operations and stores
504 \[
505 \boxed{
506 z+2nr
507 }
508 \]
509 scalar values, apart from sparse index metadata.
510 \end{proposition}
511 
512 \begin{proof}\phantomsection\label{main:proof:013}
513 The sparse product costs \(2zq\). The two dense thin products cost \(2nrq\) each. The stored values are the nonzeros of \(S\) and the entries of the two factors.
514 \phantomsection\label{main:proof:013:end}\end{proof}
515 
516 \begin{theorem}[Circulant spectral representation]\label{main:theorem:006}\label{thm:circulant}
517 Let $n\ge1$ be an integer and $c\in\CC^n$, with indices taken modulo $n$, and define
518 \[
519 C(c)_{jk}=c_{j-k}.
520 \]
521 Let \(F\) be the unitary discrete Fourier matrix
522 \[
523 F_{jk}
524 =
525 \frac{1}{\sqrt n}
526 \exp\left(-\frac{2\pi\ii jk}{n}\right).
527 \]
528 Then
529 \[
530 \boxed{
531 C(c)
532 =
533 F^*\diag(\sqrt n\,Fc)F.
534 }
535 \]
536 The operator is determined by \(n\) scalars and can be applied through two Fourier transforms and one pointwise product.
537 \end{theorem}
538 
539 \begin{proof}\phantomsection\label{main:proof:014}
540 For each frequency \(k\), let
541 \[
542 v^{(k)}_j
543 =
544 \frac{1}{\sqrt n}
545 \exp\left(\frac{2\pi\ii jk}{n}\right).
546 \]
547 Then
548 \[
549 (C(c)v^{(k)})_j
550 =
551 \sum_{r=0}^{n-1}c_rv^{(k)}_{j-r}
552 =
553 \left(
554 \sum_{r=0}^{n-1}
555 c_r\exp\left(-\frac{2\pi\ii rk}{n}\right)
556 \right)v^{(k)}_j.
557 \]
558 The eigenvalue is \(\sqrt n(Fc)_k\). The vectors \(v^{(k)}\) are the columns of \(F^*\), which gives the formula.
559 \phantomsection\label{main:proof:014:end}\end{proof}
560 
561 The fast Fourier transform evaluates the two transforms in \(O(n\log n)\) arithmetic for the usual composite lengths \cite{CooleyTukey1965}. This is an algorithmic factorization of the exact identity in \cref{thm:circulant}.
562 
563 \begin{corollary}[Toeplitz circulant embedding]\label{main:corollary:004}\label{cor:toeplitz_embedding}
564 Let $n\ge1$ be an integer and $T\in\CC^{n\times n}$ be Toeplitz, with first column
565 \[
566 (c_0,c_1,\ldots,c_{n-1})^T
567 \]
568 and first row
569 \[
570 (c_0,r_1,\ldots,r_{n-1}).
571 \]
572 Set
573 \[
574 g
575 =
576 (c_0,c_1,\ldots,c_{n-1},0,r_{n-1},\ldots,r_1)^T
577 \in\CC^{2n}
578 \]
579 and append \(n\) zeros to \(x\in\CC^n\). The first \(n\) entries of the circulant product \(C(g)(x,0)^T\) equal \(Tx\).
580 \end{corollary}
581 
582 \begin{proof}\phantomsection\label{main:proof:015}
583 For \(0\leq j<n\), the circular convolution sum separates into indices \(k\leq j\), where the coefficient is \(c_{j-k}\), and indices \(k>j\), where the wrapped coefficient is \(r_{k-j}\). The padded coordinates do not contribute. The resulting sum is exactly the \(j\)th Toeplitz row applied to \(x\).
584 \phantomsection\label{main:proof:015:end}\end{proof}
585 
586 \begin{theorem}[Kronecker modewise application]\label{main:theorem:007}\label{thm:kronecker_apply}
587 Let \(A\in\CC^{m\times m}\), \(B\in\CC^{n\times n}\), and \(X\in\CC^{n\times m}\). With column major vectorization,
588 \[
589 \boxed{
590 (A\otimes B)\operatorname{vec}(X)
591 =
592 \operatorname{vec}(BXA^T).
593 }
594 \]
595 The right side has modeled cost
596 \[
597 \boxed{
598 2n^2m+2nm^2
599 }
600 \]
601 and stores \(n^2+m^2\) factor entries instead of \(n^2m^2\) entries.
602 \end{theorem}
603 
604 \begin{proof}\phantomsection\label{main:proof:016}
605 The \((i,j)\) entry of \(BXA^T\) is
606 \[
607 \sum_{k=1}^{n}\sum_{\ell=1}^{m}
608 B_{ik}X_{k\ell}A_{j\ell}.
609 \]
610 This is the coordinate formula for \((A\otimes B)\operatorname{vec}(X)\). Multiplying first by \(B\) and then by \(A^T\) gives the stated count.
611 \phantomsection\label{main:proof:016:end}\end{proof}
612 
613 \begin{corollary}[Equal tensor factors]\label{main:corollary:005}\label{cor:tensor_modewise}
614 For integers $k,d\ge0$, a product of $k$ square factors of order $d$, acting on a tensor with $d^k$ entries, can be applied by $k$ modewise contractions, reusing intermediate slices, with modeled cost
615 \[
616 \boxed{
617 2k d^{k+1}
618 }
619 \]
620 and factor storage \(kd^2\).
621 \end{corollary}
622 
623 \begin{proof}\phantomsection\label{main:proof:017}
624 For $k=0$ there are no contractions and both modeled work and factor storage are zero. For $k\ge1$, one mode has $d^{k-1}$ fibers of length $d$. Multiplying each fiber by a \(d\) by \(d\) matrix costs \(2d^2\), so one mode costs \(2d^{k+1}\). Sum over the \(k\) modes.
625 \phantomsection\label{main:proof:017:end}\end{proof}
626 
627 \begin{theorem}[Walsh--Hadamard recursion]\label{main:theorem:008}\label{thm:walsh_hadamard}
628 Let \(N=2^m\), and define
629 \[
630 H_1=[1],
631 \qquad
632 H_{2N}
633 =
634 \begin{bmatrix}
635 H_N&H_N\\
636 H_N&-H_N
637 \end{bmatrix}.
638 \]
639 Then
640 \[
641 H_NH_N^T=NI.
642 \]
643 The product \(H_Nx\) can be evaluated in exactly \(N\log_2N\) additions or subtractions through the recursion.
644 \end{theorem}
645 
646 \begin{proof}\phantomsection\label{main:proof:018}
647 The orthogonality identity follows by induction from block multiplication. One recursive stage forms \(x_0+x_1\) and \(x_0-x_1\), which uses \(N\) additions or subtractions, and then applies two transforms of size \(N/2\). Thus
648 \[
649 T(N)=2T(N/2)+N,
650 \qquad
651 T(1)=0,
652 \]
653 whose solution is \(T(N)=N\log_2N\).
654 \phantomsection\label{main:proof:018:end}\end{proof}
655 
656 \section{Exact low rank algebra and recompression}\label{main:section:005}\label{sec:bdlr}
657 
658 \begin{figure}[t]
659 \centering
660 \includegraphics[width=.68\textwidth]{05_bdlr_structure.pdf}
661 \caption{A regular diagonal packet structure is coupled by global low rank factors.}
662 \label{fig:bdlr_structure}
663 \end{figure}
664 
665 \subsection{Products and inverses}\label{main:subsection:003}
666 
667 \begin{theorem}[Exact BDLR product]\label{main:theorem:009}\label{thm:bdlr_product}
668 Let
669 \[
670 A=D_A+U_AV_A^*,
671 \qquad
672 B=D_B+U_BV_B^*
673 \]
674 use the same diagonal packet layout. Then
675 \[
676 AB=D_AD_B+U_CV_C^*,
677 \]
678 where
679 \[
680 U_C
681 =
682 \begin{bmatrix}
683 D_AU_B+U_A(V_A^*U_B)&U_A
684 \end{bmatrix}
685 \]
686 and
687 \[
688 V_C
689 =
690 \begin{bmatrix}
691 V_B&D_B^*V_A
692 \end{bmatrix}.
693 \]
694 The unrecompressed correction rank is at most \(r_A+r_B\).
695 \end{theorem}
696 
697 \begin{proof}\phantomsection\label{main:proof:019}
698 Expand the product:
699 \[
700 AB
701 =
702 D_AD_B
703 +D_AU_BV_B^*
704 +U_AV_A^*D_B
705 +U_A(V_A^*U_B)V_B^*.
706 \]
707 The second and fourth terms share the right factor \(V_B^*\). The third term equals
708 \[
709 U_A(D_B^*V_A)^*.
710 \]
711 Concatenation gives the stated factors.
712 \phantomsection\label{main:proof:019:end}\end{proof}
713 
714 \begin{corollary}[Rank growth in a product chain]\label{main:corollary:006}\label{cor:rank_growth}
715 Let
716 \[
717 A_j=D_j+U_jV_j^*,
718 \qquad
719 j=1,\ldots,m,
720 \]
721 share one diagonal packet layout. Before recompression, the product
722 \[
723 A_m\cdots A_1
724 \]
725 has a BDLR representation whose correction rank is at most
726 \[
727 \boxed{
728 \sum_{j=1}^m r_j.
729 }
730 \]
731 \end{corollary}
732 
733 \begin{proof}\phantomsection\label{main:proof:020}
734 Apply \cref{thm:bdlr_product} inductively. Each new factor adds at most its own correction rank.
735 \phantomsection\label{main:proof:020:end}\end{proof}
736 
737 \begin{theorem}[Woodbury inversion]\label{main:theorem:010}\label{thm:woodbury}
738 Assume \(D\) and
739 \[
740 S=I+V^*D^{-1}U
741 \]
742 are invertible. Then
743 \[
744 \boxed{
745 (D+UV^*)^{-1}
746 =
747 D^{-1}
748 -D^{-1}US^{-1}V^*D^{-1}.
749 }
750 \]
751 \end{theorem}
752 
753 \begin{proof}\phantomsection\label{main:proof:021}
754 Multiply the proposed inverse by \(D+UV^*\). The coefficient of \(V^*D^{-1}\) is
755 \[
756 U-US^{-1}(I+V^*D^{-1}U)=0.
757 \]
758 The remaining term is the identity. The right inverse calculation is the same.
759 \phantomsection\label{main:proof:021:end}\end{proof}
760 
761 The formula and its history are discussed in \cite{Hager1989}. It reduces one large solve to packet solves with \(D\) and one dense solve with the small matrix \(S\).
762 
763 \begin{corollary}[Inverse norm bound]\label{main:corollary:007}\label{cor:inverse_norm}
764 Under the hypotheses of \cref{thm:woodbury},
765 \[
766 \norm{(D+UV^*)^{-1}}
767 \leq
768 \norm{D^{-1}}
769 +
770 \norm{D^{-1}U}
771 \norm{S^{-1}}
772 \norm{V^*D^{-1}}.
773 \]
774 \end{corollary}
775 
776 \begin{proof}\phantomsection\label{main:proof:022}
777 Apply the triangle inequality and submultiplicativity to the Woodbury formula.
778 \phantomsection\label{main:proof:022:end}\end{proof}
779 
780 \subsection{Optimal factor recompression}\label{main:subsection:004}
781 
782 \begin{standardfact}[Eckart--Young--Mirsky theorem]\label{main:standardfact:001}\label{fact:eym}
783 If \(M\) has singular values
784 \[
785 \sigma_1\geq\sigma_2\geq\cdots,
786 \]
787 pad this finite sequence with zeros. For each integer $k\ge0$, its truncated singular value decomposition $M_k$ is a best rank-at-most-$k$ approximation. It satisfies
788 \[
789 \norm{M-M_k}_2=\sigma_{k+1}
790 \]
791 and
792 \[
793 \norm{M-M_k}_{\rm F}
794 =
795 \left(
796 \sum_{j>k}\sigma_j^2
797 \right)^{1/2}.
798 \]
799 \end{standardfact}
800 
801 This is classical \cite{EckartYoung1936,Mirsky1960}. The next theorem shows that an \(n\) dimensional singular value decomposition is unnecessary when a correction is already factored.
802 
803 \begin{theorem}[Small core recompression]\label{main:theorem:011}\label{thm:small_core_recompression}
804 Let
805 \[
806 L\in\CC^{m\times r},
807 \qquad
808 R\in\CC^{n\times r}.
809 \]
810 Take compact QR factorizations
811 \[
812 L=Q_LR_L,
813 \qquad
814 R=Q_RR_R,
815 \]
816 where the columns of \(Q_L\) and \(Q_R\) are orthonormal. Form the small core
817 \[
818 G=R_LR_R^*
819 \]
820 and let
821 \[
822 G=W\Sigma Z^*
823 \]
824 be an ordered thin singular value decomposition, with nonnegative nonincreasing singular values. For an integer $k\ge0$, retain $\min\{k,d\}$ columns, where $d$ is the number of supplied thin singular directions. Tail singular values beyond $d$ are zero. Define
825 \[
826 L_k=Q_LW_k\Sigma_k^{1/2},
827 \qquad
828 R_k=Q_RZ_k\Sigma_k^{1/2}.
829 \]
830 Then
831 \[
832 L_kR_k^*
833 \]
834 is a best rank at most \(k\) approximation to \(LR^*\). Its spectral and Frobenius errors are the singular value tails in \cref{fact:eym}. The core dimensions equal the column ranks of the two input factors. A noncompact reduced QR may retain extra zero singular directions without changing the conclusion.
835 \end{theorem}
836 
837 \begin{proof}\phantomsection\label{main:proof:023}
838 The exact factor product is
839 \[
840 LR^*
841 =
842 Q_LGQ_R^*.
843 \]
844 The two outer maps are isometries on their column spaces. They preserve the nonzero singular values and the spectral and Frobenius norms of the core. For every ambient rank-at-most-$k$ competitor $B$, its compression $Q_L^*BQ_R$ has rank at most $k$, and compression cannot increase either norm. Thus the core best-approximation inequality also bounds the error of every ambient competitor from below. Truncating the singular value decomposition of \(G\) and applying \cref{fact:eym} gives the result. The compact QR dimensions equal the ranks of the factor column spaces.
845 \phantomsection\label{main:proof:023:end}\end{proof}
846 
847 \begin{figure}[t]
848 \centering
849 \includegraphics[width=.76\textwidth]{06_recompression_spectrum.pdf}
850 \caption{The approximation decision is made from the singular values of a small core whose dimension is the unrecompressed rank.}
851 \label{fig:recompression}
852 \end{figure}
853 
854 \begin{remark}\label{main:remark:002}\label{rem:recompression_work}
855 The large matrices enter through QR factorizations of the factors. The singular value problem has dimension at most the unrecompressed rank. Updating and recompressing low rank factors is standard numerical linear algebra \cite{GolubKahan1965,BunchNielsen1978}.
856 \end{remark}
857 
858 \begin{theorem}[Product chain perturbation]\label{main:theorem:012}\label{thm:product_chain_error}
859 Let \(A_1,\ldots,A_m\) be exact factors and let
860 \[
861 \widehat A_j=A_j+E_j.
862 \]
863 Then
864 \[
865 \boxed{
866 \left\|
867 \widehat A_m\cdots\widehat A_1
868 -
869 A_m\cdots A_1
870 \right\|
871 \leq
872 \sum_{j=1}^m
873 \left(
874 \prod_{i=j+1}^m\norm{\widehat A_i}
875 \right)
876 \norm{E_j}
877 \left(
878 \prod_{i=1}^{j-1}\norm{A_i}
879 \right).
880 }
881 \]
882 Empty products are one.
883 \end{theorem}
884 
885 \begin{proof}\phantomsection\label{main:proof:024}
886 Use the telescoping identity
887 \[
888 \widehat A_m\cdots\widehat A_1
889 -A_m\cdots A_1
890 =
891 \sum_{j=1}^m
892 \widehat A_m\cdots\widehat A_{j+1}
893 E_j
894 A_{j-1}\cdots A_1.
895 \]
896 Take norms term by term.
897 \phantomsection\label{main:proof:024:end}\end{proof}
898 
899 \begin{theorem}[Certified repeated recompression]\label{main:theorem:013}\label{cor:recompression_chain}
900 Let
901 \[
902 P_0=I,
903 \qquad
904 P_j=A_jP_{j-1},
905 \]
906 be exact partial products. Suppose an implemented chain satisfies
907 \[
908 \widehat P_0=I,
909 \qquad
910 \widehat P_j
911 =
912 A_j\widehat P_{j-1}+E_j,
913 \]
914 where \(E_j\) is the additive error introduced after the \(j\)th multiplication, for example by factor recompression. Then
915 \[
916 \boxed{
917 \norm{\widehat P_m-P_m}
918 \leq
919 \sum_{j=1}^{m}
920 \left(
921 \prod_{i=j+1}^{m}\norm{A_i}
922 \right)
923 \norm{E_j}.
924 }
925 \]
926 If \(\norm{E_j}\leq\delta_j\), the same formula with \(\delta_j\) gives a certified final error.
927 \end{theorem}
928 
929 \begin{proof}\phantomsection\label{main:proof:025}
930 Set \(D_j=\widehat P_j-P_j\). Then
931 \[
932 D_j=A_jD_{j-1}+E_j,
933 \qquad
934 D_0=0.
935 \]
936 Repeated substitution gives
937 \[
938 D_m
939 =
940 \sum_{j=1}^{m}
941 A_mA_{m-1}\cdots A_{j+1}E_j.
942 \]
943 Taking norms yields the result.
944 \phantomsection\label{main:proof:025:end}\end{proof}
945 
946 \section{Sparse transfer graphs and schedules}\label{main:section:006}\label{sec:graphs}
947 
948 \subsection{Support under algebraic operations}\label{main:subsection:005}
949 
950 \begin{theorem}[Path support of powers]\label{main:theorem:014}\label{thm:path_support}
951 For every positive integer \(m\),
952 \[
953 P_\gamma A^mP_\alpha
954 \]
955 is the sum over directed walks (vertices may repeat)
956 \[
957 \alpha=\alpha_0
958 \longrightarrow
959 \alpha_1
960 \longrightarrow\cdots\longrightarrow
961 \alpha_m=\gamma
962 \]
963 in \(G_A\), with contribution
964 \[
965 A_{\alpha_m\alpha_{m-1}}
966 \cdots
967 A_{\alpha_1\alpha_0}.
968 \]
969 If no such path exists, the block vanishes.
970 \end{theorem}
971 
972 \begin{proof}\phantomsection\label{main:proof:026}
973 Insert the packet identity between every two consecutive factors of \(A\). Every nonzero term determines one path of length \(m\).
974 \phantomsection\label{main:proof:026:end}\end{proof}
975 
976 \begin{corollary}[Polynomial support]\label{main:corollary:008}\label{cor:polynomial_support}
977 If
978 \[
979 p(z)=\sum_{j=0}^dc_jz^j,
980 \]
981 then
982 \[
983 \supp(p(A))
984 \]
985 is contained in the union of the support relations of \(A^j\) for \(0\leq j\leq d\).
986 \end{corollary}
987 
988 \begin{proof}\phantomsection\label{main:proof:027}
989 Apply \cref{thm:path_support} to every power in the polynomial.
990 \phantomsection\label{main:proof:027:end}\end{proof}
991 
992 \begin{corollary}[Acyclic transfer nilpotence]\label{main:corollary:009}\label{cor:acyclic_nilpotence}
993 Let \(K\) contain only strict transfer blocks. If its support graph is acyclic and its longest directed path has \(h\) edges, then
994 \[
995 K^{h+1}=0.
996 \]
997 \end{corollary}
998 
999 \begin{proof}\phantomsection\label{main:proof:028}
1000 There is no directed path of length \(h+1\). Every block of \(K^{h+1}\) therefore vanishes by \cref{thm:path_support}.
1001 \phantomsection\label{main:proof:028:end}\end{proof}
1002 
1003 \begin{figure}[t]
1004 \centering
1005 \includegraphics[width=.9\textwidth]{07_path_support.pdf}
1006 \caption{Composition creates possible support along directed paths. Cancellation can remove an allowed block, so the graph gives an inclusion rather than an equality.}
1007 \label{fig:path_support}
1008 \end{figure}
1009 
1010 \subsection{Block triangular solves}\label{main:subsection:006}
1011 
1012 \begin{definition}[Topological level]\label{main:definition:009}\label{def:topological_level}
1013 Let \(G\) be a directed acyclic graph. The first level contains its source vertices. Remove them. The next level contains the new sources. Continue until every vertex has been assigned.
1014 \end{definition}
1015 
1016 \begin{theorem}[Level schedule for a triangular solve]\label{main:theorem:015}\label{thm:level_schedule}
1017 Let \(A=D+K\) be block lower triangular after a topological ordering. Assume every diagonal block \(D_\alpha\) is invertible. The solve
1018 \[
1019 Ax=b
1020 \]
1021 can be completed in one sequential stage per topological level. Blocks in the same level may be solved independently after all earlier levels are known.
1022 \end{theorem}
1023 
1024 \begin{proof}\phantomsection\label{main:proof:029}
1025 The equation in packet \(\beta\) is
1026 \[
1027 D_\beta x_\beta
1028 =
1029 b_\beta
1030 -
1031 \sum_{\alpha\to\beta}K_{\beta\alpha}x_\alpha.
1032 \]
1033 Every predecessor of \(\beta\) lies in an earlier level. Once those values are available, \(x_\beta\) is obtained from the local solve with \(D_\beta\). Vertices in the same level have no directed edge between them, so their solves are independent.
1034 \phantomsection\label{main:proof:029:end}\end{proof}
1035 
1036 \begin{corollary}[Sequential depth]\label{main:corollary:010}\label{cor:sequential_depth}
1037 For a nonempty finite graph, the number of sequential stages in the level schedule equals one plus the length of the longest directed path. The empty graph has depth zero.
1038 \end{corollary}
1039 
1040 \begin{proof}\phantomsection\label{main:proof:030}
1041 The recursive level of a vertex is one plus the largest level of its predecessors. For a nonempty graph, the largest level therefore counts the vertices on a longest path. For an empty graph, no task is scheduled and the depth is zero.
1042 \phantomsection\label{main:proof:030:end}\end{proof}
1043 
1044 \begin{figure}[t]
1045 \centering
1046 \includegraphics[width=.82\textwidth]{08_topological_schedule.pdf}
1047 \caption{Packets in the same topological level can be processed concurrently in a triangular solve.}
1048 \label{fig:topological_schedule}
1049 \end{figure}
1050 
1051 \subsection{Fill under elimination}\label{main:subsection:007}
1052 
1053 \begin{theorem}[Inverse support lies in the transitive closure]\label{main:theorem:016}\label{thm:inverse_support}
1054 Let \(A\) be an invertible packet operator. Reachability here includes a zero-length path from a packet to itself. If there is no directed path from packet \(\alpha\) to packet \(\beta\) in \(G_A\), then
1055 \[
1056 \boxed{
1057 P_\beta A^{-1}P_\alpha=0.
1058 }
1059 \]
1060 \end{theorem}
1061 
1062 \begin{proof}\phantomsection\label{main:proof:031}
1063 Let \(R\) be the set of packets reachable from \(\alpha\), including \(\alpha\), and set
1064 \[
1065 X_R=\bigoplus_{\gamma\in R}X_\gamma.
1066 \]
1067 No active edge leaves \(R\), so \(AX_R\subseteq X_R\). The restriction of \(A\) to the finite dimensional space \(X_R\) is injective. It is therefore onto \(X_R\). Hence \(A^{-1}X_R\subseteq X_R\). Since \(\beta\notin R\), the stated block of \(A^{-1}\) vanishes.
1068 \phantomsection\label{main:proof:031:end}\end{proof}
1069 
1070 \begin{proposition}[Schur support inclusion]\label{main:proposition:006}\label{prop:schur_fill}
1071 Split the packets into retained indices \(R\) and hidden indices \(H\). Suppose the hidden block \(A_{HH}\) is invertible. The retained Schur complement is
1072 \[
1073 S
1074 =
1075 A_{RR}
1076 -
1077 A_{RH}A_{HH}^{-1}A_{HR}.
1078 \]
1079 A retained block $({\beta},{\alpha})$ can be active only if it was already active in $A_{RR}$, or if there is a directed walk from $\alpha$ through hidden packets to $\beta$. The internal hidden path may have length zero, so entrance and exit may use the same hidden packet.
1080 \end{proposition}
1081 
1082 \begin{proof}\phantomsection\label{main:proof:032}
1083 The correction block is
1084 \[
1085 P_\beta A_{RH}A_{HH}^{-1}A_{HR}P_\alpha.
1086 \]
1087 It can be nonzero only when \(A_{HR}P_\alpha\) enters a hidden packet, the inverse hidden operator connects hidden packets, and \(P_\beta A_{RH}\) exits to \(\beta\). By \cref{thm:inverse_support}, the hidden inverse can connect only packets joined by a hidden directed path. Cancellation may still remove an allowed block.
1088 \phantomsection\label{main:proof:032:end}\end{proof}
1089 
1090 Sparse elimination and fill are treated systematically in \cite{Davis2006}. The packet graph records the same structural issue at block scale.
1091 
1092 \section{Krylov methods and packet preconditioning}\label{main:section:007}\label{sec:krylov}
1093 
1094 \begin{definition}[Krylov space]\label{main:definition:010}\label{def:krylov}
1095 For a linear operator \(A\), a vector \(v\), and a positive integer \(m\), define
1096 \[
1097 \Kcal_m(A,v)
1098 =
1099 \operatorname{span}
1100 \{v,Av,\ldots,A^{m-1}v\}.
1101 \]
1102 \end{definition}
1103 
1104 \begin{proposition}[Invariant packet Krylov decomposition]\label{main:proposition:007}\label{prop:krylov_invariant}
1105 Assume
1106 \[
1107 [A,P_\alpha]=0
1108 \]
1109 for every packet. If
1110 \[
1111 v=\sum_\alpha v_\alpha,
1112 \qquad
1113 v_\alpha=P_\alpha v,
1114 \]
1115 then
1116 \[
1117 A^jv
1118 =
1119 \sum_\alpha A_\alpha^jv_\alpha
1120 \]
1121 for every \(j\geq0\), where
1122 \[
1123 A_\alpha=P_\alpha AP_\alpha.
1124 \]
1125 Hence
1126 \[
1127 \Kcal_m(A,v)
1128 \subseteq
1129 \bigoplus_\alpha
1130 \Kcal_m(A_\alpha,v_\alpha).
1131 \]
1132 \end{proposition}
1133 
1134 \begin{proof}\phantomsection\label{main:proof:033}
1135 Commutation gives
1136 \[
1137 A^jP_\alpha=P_\alpha A^j.
1138 \]
1139 Apply this identity to each packet component of \(v\).
1140 \phantomsection\label{main:proof:033:end}\end{proof}
1141 
1142 \begin{theorem}[Approximate packet Krylov bound]\label{main:theorem:017}\label{thm:approx_krylov}
1143 Write
1144 \[
1145 A=D+K,
1146 \]
1147 where \(D\) is packet diagonal. For every positive integer \(j\),
1148 \[
1149 \boxed{
1150 \norm{A^j-D^j}
1151 \leq
1152 jM^{j-1}\norm{K},
1153 }
1154 \]
1155 where
1156 \[
1157 M=\max\{\norm{A},\norm{D}\}.
1158 \]
1159 Consequently,
1160 \[
1161 \dist(A^jv,\bigoplus_\alpha\Kcal_{j+1}(D_\alpha,v_\alpha))
1162 \leq
1163 jM^{j-1}\norm{K}\norm{v}.
1164 \]
1165 \end{theorem}
1166 
1167 \begin{proof}\phantomsection\label{main:proof:034}
1168 Use the telescoping identity
1169 \[
1170 A^j-D^j
1171 =
1172 \sum_{k=0}^{j-1}
1173 A^{j-1-k}KD^k.
1174 \]
1175 Each term has norm at most \(M^{j-1}\norm{K}\). The vector \(D^jv\) lies in the direct sum of local Krylov spaces.
1176 \phantomsection\label{main:proof:034:end}\end{proof}
1177 
1178 \begin{theorem}[Block Jacobi condition bound]\label{main:theorem:018}\label{thm:block_jacobi}
1179 Let
1180 \[
1181 A=D+K
1182 \]
1183 act on a nonempty finite coordinate space and be Hermitian positive definite with positive definite block diagonal part $D$. Here $D^{-1/2}$ is the actual inverse of the positive square root of $D$. Suppose
1184 \[
1185 E=D^{-1/2}KD^{-1/2}
1186 \]
1187 is Hermitian and
1188 \[
1189 \norm{E}_2\leq\eta<1.
1190 \]
1191 Then
1192 \[
1193 \boxed{
1194 \cond_2(D^{-1/2}AD^{-1/2})
1195 \leq
1196 \frac{1+\eta}{1-\eta}.
1197 }
1198 \]
1199 \end{theorem}
1200 
1201 \begin{proof}\phantomsection\label{main:proof:035}
1202 For every unit vector \(x\),
1203 \[
1204 \abs{x^*Ex}\leq\eta.
1205 \]
1206 The Rayleigh quotient of \(I+E\) lies in \([1-\eta,1+\eta]\). The quotient of the extreme eigenvalue bounds gives the result.
1207 \phantomsection\label{main:proof:035:end}\end{proof}
1208 
1209 \begin{standardfact}[Preconditioned conjugate gradient convergence]\label{main:standardfact:002}\label{fact:cg}
1210 Let $A$ and $M$ be Hermitian positive definite on a nonempty finite coordinate space, and set
1211 \[
1212 \kappa
1213 =
1214 \cond_2(M^{-1/2}AM^{-1/2}).
1215 \]
1216 In exact arithmetic, the energy norm error of preconditioned conjugate gradient satisfies
1217 \[
1218 \norm{x_m-x_*}_A
1219 \leq
1220 2
1221 \left(
1222 \frac{\sqrt\kappa-1}{\sqrt\kappa+1}
1223 \right)^m
1224 \norm{x_0-x_*}_A.
1225 \]
1226 \end{standardfact}
1227 
1228 This standard estimate is discussed in \cite{Saad2003}. Combining it with \cref{thm:block_jacobi} yields an explicit iteration certificate when the normalized cross packet coupling is small.
1229 
1230 \begin{figure}[t]
1231 \centering
1232 \includegraphics[width=.75\textwidth]{09_packet_preconditioner.pdf}
1233 \caption{The block diagonal packet operator serves as a reusable local preconditioner in a weakly coupled positive system.}
1234 \label{fig:preconditioner}
1235 \end{figure}
1236 
1237 \begin{corollary}[Certified packet preconditioning]\label{main:corollary:011}\label{cor:certified_pcg}
1238 Under the hypotheses of \cref{thm:block_jacobi}, the convergence factor in \cref{fact:cg} is bounded by the value obtained from
1239 \[
1240 \kappa
1241 =
1242 \frac{1+\eta}{1-\eta}.
1243 \]
1244 \end{corollary}
1245 
1246 \begin{proof}\phantomsection\label{main:proof:036}
1247 Insert the condition bound from \cref{thm:block_jacobi} into the monotone convergence factor in \cref{fact:cg}.
1248 \phantomsection\label{main:proof:036:end}\end{proof}
1249 
1250 \section{Packetwise mixed precision}\label{main:section:008}\label{sec:precision}
1251 
1252 \subsection{One operator}\label{main:subsection:008}
1253 
1254 \begin{theorem}[Packet perturbation certificate]\label{main:theorem:019}\label{thm:packet_error}
1255 Let
1256 \[
1257 A=
1258 \sum_{e\in E}A_e,
1259 \qquad
1260 \widehat A
1261 =
1262 \sum_{e\in E}(A_e+E_e),
1263 \]
1264 where edge \(e\) starts at packet \(\alpha(e)\), and assume
1265 \[
1266 A_e=A_eP_{\alpha(e)},
1267 \qquad
1268 E_e=E_eP_{\alpha(e)}.
1269 \]
1270 Then
1271 \[
1272 \boxed{
1273 \norm{(\widehat A-A)x}
1274 \leq
1275 \sum_{e\in E}
1276 \norm{E_e}
1277 \norm{P_{\alpha(e)}x}.
1278 }
1279 \]
1280 If
1281 \[
1282 \norm{E_e}
1283 \leq
1284 \eps_e\norm{A_e},
1285 \]
1286 then
1287 \[
1288 \boxed{
1289 \norm{(\widehat A-A)x}
1290 \leq
1291 \sum_{e\in E}
1292 \eps_e\norm{A_e}
1293 \norm{P_{\alpha(e)}x}.
1294 }
1295 \]
1296 \end{theorem}
1297 
1298 \begin{proof}\phantomsection\label{main:proof:037}
1299 Exact packet evaluation gives
1300 \[
1301 (\widehat A-A)x
1302 =
1303 \sum_eE_eP_{\alpha(e)}x.
1304 \]
1305 Apply the triangle inequality and the operator norm bound to each term.
1306 \phantomsection\label{main:proof:037:end}\end{proof}
1307 
1308 \begin{standardfact}[Floating point matrix product model]\label{main:standardfact:003}\label{fact:rounding}
1309 For real $A\in\RR^{m\times p}$ and $B\in\RR^{p\times n}$, consider ordinary dot-product matrix multiplication in uniform precision with unit roundoff $u$, no overflow or underflow, and $pu<1$. Its computed output $\widehat C$ satisfies the componentwise forward-error bound
1310 \[
1311 |\widehat C-AB|\le\gamma_p|A||B|,
1312 \qquad\gamma_p=\frac{pu}{1-pu}.
1313 \]
1314 Absolute values and the inequality are entrywise. Complex arithmetic, fused kernels, mixed accumulation, and tensor kernels require their own stated constants. This fact does not assert one common small input perturbation for every column, and an input-dependent rounded action is not automatically a fixed linear operator.
1315 \end{standardfact}
1316 
1317 Detailed forms and their hypotheses are given in \cite{Higham2002,HighamMary2021}. For a fixed computed matrix, the bound can certify its error matrix. A rounded application to varying inputs needs a separately justified uniform certificate before it can be used as one fixed operator error in \cref{thm:packet_error}.
1318 
1319 \subsection{Several operators}\label{main:subsection:009}
1320 
1321 \begin{corollary}[Mixed precision product chain]\label{main:corollary:012}\label{cor:mixed_chain}
1322 If each fixed linear factor $A_j$ is replaced by a fixed packetwise linear approximation $\widehat A_j$, write $E_j=\widehat A_j-A_j=\sum_e\Delta_{j,e}P_{\alpha(j,e)}$. A uniform certificate is
1323 \[
1324 \|E_j\|\le\delta_j:=\sum_e\|\Delta_{j,e}\|\,\|P_{\alpha(j,e)}\|.
1325 \]
1326 For orthogonal packet projectors, $\delta_j\le\sum_e\|\Delta_{j,e}\|$. The final operator error is bounded by \cref{thm:product_chain_error} with these uniform $\delta_j$. The corresponding vector error is this operator bound multiplied by the input norm.
1327 \end{corollary}
1328 
1329 \begin{proof}\phantomsection\label{main:proof:038}
1330 Take the supremum over unit inputs in \cref{thm:packet_error}, using $\|P_{\alpha}x\|\le\|P_{\alpha}\|\|x\|$. Substitute the resulting uniform operator error bounds into \cref{thm:product_chain_error}. The action estimate follows from submultiplicativity.
1331 \phantomsection\label{main:proof:038:end}\end{proof}
1332 
1333 \subsection{Continuous allocation}\label{main:subsection:010}
1334 
1335 \begin{theorem}[Continuous precision allocation]\label{main:theorem:020}\label{thm:continuous_precision}
1336 Let $E$ be a finite nonempty edge set, let $s_e>0$ be sensitivities and $c_e>0$ be costs, and let $\tau>0$. Consider
1337 \[
1338 \min_{\eps_e>0}
1339 \sum_e c_e\log\frac1{\eps_e}
1340 \]
1341 subject to
1342 \[
1343 \sum_es_e\eps_e=\tau.
1344 \]
1345 The unique minimizer is
1346 \[
1347 \boxed{
1348 \eps_e
1349 =
1350 \frac{\tau c_e}{s_e\sum_fc_f}.
1351 }
1352 \]
1353 \end{theorem}
1354 
1355 \begin{proof}\phantomsection\label{main:proof:039}
1356 The Lagrangian derivative gives
1357 \[
1358 -\frac{c_e}{\eps_e}+\lambda s_e=0.
1359 \]
1360 Thus
1361 \[
1362 \eps_e=\frac{c_e}{\lambda s_e}.
1363 \]
1364 The constraint gives
1365 \[
1366 \lambda=\frac{\sum_fc_f}{\tau}.
1367 \]
1368 The objective is strictly convex on the positive orthant, so the stationary point is the unique minimizer.
1369 \phantomsection\label{main:proof:039:end}\end{proof}
1370 
1371 \begin{theorem}[Clipped precision allocation]\label{main:theorem:021}\label{thm:clipped_precision}
1372 Retain the finite nonempty edge set, strictly positive sensitivities and costs, and add strictly positive caps $m_e>0$ with upper bounds
1373 \[
1374 0<\eps_e\leq m_e.
1375 \]
1376 Assume
1377 \[
1378 0<\tau\leq\sum_es_em_e.
1379 \]
1380 There is a unique \(\lambda>0\), except at the maximal tolerance endpoint, such that
1381 \[
1382 \boxed{
1383 \eps_e
1384 =
1385 \min\left\{
1386 m_e,
1387 \frac{c_e}{\lambda s_e}
1388 \right\}
1389 }
1390 \]
1391 and
1392 \[
1393 \sum_es_e\eps_e=\tau.
1394 \]
1395 This vector minimizes the same objective over the clipped feasible set.
1396 \end{theorem}
1397 
1398 \begin{proof}\phantomsection\label{main:proof:040}
1399 The Karush--Kuhn--Tucker conditions give the unconstrained form on every nonactive coordinate and equality with \(m_e\) on every active upper bound. The weighted sum
1400 \[
1401 \phi(\lambda)
1402 =
1403 \sum_es_e
1404 \min\left\{
1405 m_e,
1406 \frac{c_e}{\lambda s_e}
1407 \right\}
1408 \]
1409 is continuous and nonincreasing. It decreases from the maximal feasible value to zero. For a strictly submaximal budget, once any coordinate is unsaturated, $\phi$ is strictly decreasing, which gives a unique multiplier. At maximal budget the only feasible vector is $(m_e)_e$ and every $0<\lambda\le\min_e c_e/(s_em_e)$ produces it. Thus the minimizing vector is unique at this endpoint while its multiplier need not be. Strict convexity gives uniqueness of the minimizer in every case.
1410 \phantomsection\label{main:proof:040:end}\end{proof}
1411 
1412 \begin{figure}[t]
1413 \centering
1414 \includegraphics[width=.75\textwidth]{10_precision_allocation.pdf}
1415 \caption{More sensitive packet edges receive tighter continuous tolerances.}
1416 \label{fig:precision}
1417 \end{figure}
1418 
1419 \subsection{Discrete data types}\label{main:subsection:011}
1420 
1421 \begin{definition}[Discrete precision option]\label{main:definition:011}\label{def:dtype_option}
1422 For edge \(e\), a discrete precision option \(k\) has a measured or modeled cost \(c_{e,k}\) and a nonnegative integer error weight \(w_{e,k}\). Exactly one option must be selected on every edge.
1423 \end{definition}
1424 
1425 \begin{theorem}[Exact discrete precision program]\label{main:theorem:022}\label{thm:discrete_precision}
1426 Let the total integer error budget be \(B\). Define
1427 \[
1428 F_j(b)
1429 \]
1430 as the minimum cost for assigning precisions to the first \(j\) edges with total weight exactly \(b\). Then
1431 \[
1432 F_0(0)=0,
1433 \qquad
1434 F_0(b)=+\infty
1435 \quad\text{for }b\neq0.
1436 \]
1437 Set \(F_{j-1}(u)=+\infty\) whenever \(u<0\) or the state is unreachable. Then
1438 \[
1439 \boxed{
1440 F_j(b)
1441 =
1442 \min_k
1443 \left\{
1444 F_{j-1}(b-w_{j,k})+c_{j,k}
1445 \right\}.
1446 }
1447 \]
1448 The minimum of \(F_m(b)\) over \(0\leq b\leq B\) is the exact minimum cost for the quantized model.
1449 \end{theorem}
1450 
1451 \begin{proof}\phantomsection\label{main:proof:041}
1452 Every feasible assignment for the first \(j\) edges ends with one option \(k\) on edge \(j\). Removing it leaves an optimal or suboptimal assignment for the first \(j-1\) edges with residual weight \(b-w_{j,k}\). Replacing a suboptimal prefix by an optimal one can only lower cost. The recurrence therefore has optimal substructure. Induction on \(j\) proves the result.
1453 \phantomsection\label{main:proof:041:end}\end{proof}
1454 
1455 \begin{remark}\label{main:remark:003}\label{rem:dtype_boundary}
1456 The discrete program is exact only for the supplied options and quantized error weights. Real kernels also depend on accumulation type, scaling, overflow, underflow, tensor core availability, and shape restrictions. Mixed precision iterative refinement has its own conditioning requirements \cite{CarsonHigham2018}.
1457 \end{remark}
1458 
1459 \section{Discovery and exact packet compilation}\label{main:section:009}\label{sec:compiler}
1460 
1461 \subsection{Exact spectral discovery}\label{main:subsection:012}
1462 
1463 \begin{standardfact}[Simultaneous spectral decomposition]\label{main:standardfact:004}\label{fact:simultaneous}
1464 A finite family of pairwise commuting normal matrices admits a common orthonormal basis of joint eigenvectors. The joint eigenspaces associated with distinct eigenvalue tuples give the coarsest orthogonal packet system on which every matrix in the family acts as a scalar.
1465 \end{standardfact}
1466 
1467 \begin{proposition}[Exact block discovery certificate]\label{main:proposition:008}\label{prop:block_discovery_certificate}
1468 For an orthogonal packet system \((P_\alpha)\), define the block diagonal compression
1469 \[
1470 \Pi_{\Pcal}(A)
1471 =
1472 \sum_\alpha P_\alpha AP_\alpha.
1473 \]
1474 Then
1475 \[
1476 \boxed{
1477 \norm{A-\Pi_{\Pcal}(A)}_{\rm F}^2
1478 =
1479 \sum_{\alpha\neq\beta}
1480 \norm{P_\beta AP_\alpha}_{\rm F}^2.
1481 }
1482 \]
1483 In operator norm,
1484 \[
1485 \boxed{
1486 \norm{A-\Pi_{\Pcal}(A)}_2
1487 \leq
1488 \sum_{\alpha\neq\beta}
1489 \norm{P_\beta AP_\alpha}_2.
1490 }
1491 \]
1492 Thus any proposed packet discovery has an exact finite residual certificate after its projectors are known.
1493 \end{proposition}
1494 
1495 \begin{proof}\phantomsection\label{main:proof:042}
1496 Distinct matrix blocks are orthogonal in the Hilbert--Schmidt inner product. Their squared Frobenius norms therefore add. The operator norm estimate follows from the triangle inequality.
1497 \phantomsection\label{main:proof:042:end}\end{proof}
1498 
1499 This supplies one exact discovery regime and one exact a posteriori certificate. Perturbation of spectral packets requires a gap. The Davis--Kahan theorem gives quantitative control for Hermitian spectral subspaces under such a gap \cite{DavisKahan1970}. Randomized low rank discovery supplies another controlled regime \cite{HalkoMartinssonTropp2011,MartinssonTropp2020}. Neither result gives globally optimal packet discovery for arbitrary operators.
1500 
1501 \subsection{Candidate packet trees}\label{main:subsection:013}
1502 
1503 \begin{definition}[Candidate representation tree]\label{main:definition:012}\label{def:candidate_tree}
1504 A candidate representation tree is a finite rooted tree. Every node \(v\) has
1505 \[
1506 \kappa_v\geq0,
1507 \qquad
1508 \epsilon_v\geq0,
1509 \]
1510 representing the cost and error of keeping \(v\) as one leaf. Every internal node also has a split overhead \(\sigma_v\geq0\). An admissible pruning contains exactly one selected node on every root to terminal path.
1511 \end{definition}
1512 
1513 The cost and error of a pruning \(\Lcal\) are
1514 \[
1515 \cost(\Lcal)
1516 =
1517 \sum_{v\in\Lcal}\kappa_v
1518 +
1519 \sum_{u\text{ split}}\sigma_u
1520 \]
1521 and
1522 \[
1523 \err(\Lcal)=\sum_{v\in\Lcal}\epsilon_v.
1524 \]
1525 The model can include arithmetic, storage, launch, transfer, and precision costs inside \(\kappa_v\). It remains exact only when the declared totals are additive.
1526 
1527 \begin{theorem}[Exact Lagrangian packet compiler]\label{main:theorem:023}\label{thm:tree_compiler}
1528 Fix \(\lambda\geq0\) and define
1529 \[
1530 J_\lambda(\Lcal)
1531 =
1532 \cost(\Lcal)
1533 +
1534 \lambda\err(\Lcal).
1535 \]
1536 For every node \(v\), set
1537 \[
1538 K_\lambda(v)
1539 =
1540 \kappa_v+\lambda\epsilon_v
1541 \]
1542 and
1543 \[
1544 S_\lambda(v)
1545 =
1546 \sigma_v
1547 +
1548 \sum_{w\in\operatorname{child}(v)}F_\lambda(w).
1549 \]
1550 For a terminal node, set
1551 \[
1552 F_\lambda(v)=K_\lambda(v).
1553 \]
1554 For an internal node, the recurrence
1555 \[
1556 \boxed{
1557 F_\lambda(v)
1558 =
1559 \min\{K_\lambda(v),S_\lambda(v)\}
1560 }
1561 \]
1562 computes the exact minimum Lagrangian objective over every admissible pruning of the subtree rooted at \(v\).
1563 \end{theorem}
1564 
1565 \begin{proof}\phantomsection\label{main:proof:043}
1566 A pruning of the subtree rooted at \(v\) has two possible forms. It keeps \(v\), in which case its objective is \(K_\lambda(v)\). Otherwise it splits \(v\) and independently chooses an admissible pruning in every child subtree. By induction, the optimal child values are \(F_\lambda(w)\). The split objective is therefore \(S_\lambda(v)\). Taking the smaller value is exact.
1567 \phantomsection\label{main:proof:043:end}\end{proof}
1568 
1569 This recurrence is a computational version of best basis selection and Bellman optimality \cite{Bellman1957,CoifmanWickerhauser1992}. The node costs here may include operator application, transfer, launch, precision, and output conversion rather than only a signal entropy.
1570 
1571 \begin{figure}[t]
1572 \centering
1573 \includegraphics[width=.82\textwidth]{11_tree_compiler.pdf}
1574 \caption{The thick leaves are selected by the exact Lagrangian recurrence.}
1575 \label{fig:tree_compiler}
1576 \end{figure}
1577 
1578 \begin{corollary}[Shared compilation for an operator family]\label{main:corollary:013}\label{cor:shared_compiler}
1579 Let operators \(A^{(1)},\ldots,A^{(m)}\) use the same candidate tree. Give operator \(j\) a nonnegative reuse weight \(r_j\). Replace the local keep cost and error by
1580 \[
1581 \kappa_v
1582 =
1583 \kappa_v^{\rm setup}
1584 +
1585 \sum_{j=1}^m r_j\kappa_v^{(j)}
1586 \]
1587 and
1588 \[
1589 \epsilon_v
1590 =
1591 \sum_{j=1}^m r_j\epsilon_v^{(j)}.
1592 \]
1593 Then \cref{thm:tree_compiler} gives the exact optimal shared pruning for the additive weighted objective.
1594 \end{corollary}
1595 
1596 \begin{proof}\phantomsection\label{main:proof:044}
1597 The weighted family objective is still the sum of local leaf and split terms. Apply \cref{thm:tree_compiler} to the aggregated quantities.
1598 \phantomsection\label{main:proof:044:end}\end{proof}
1599 
1600 \begin{proposition}[Monotone Lagrangian frontier]\label{main:proposition:009}\label{prop:lagrangian_frontier}
1601 Let \(0\leq\lambda_1<\lambda_2\). Choose optimal prunings \(\Lcal_1\) and \(\Lcal_2\) for the two Lagrangian parameters. Then
1602 \[
1603 \boxed{
1604 \err(\Lcal_2)
1605 \leq
1606 \err(\Lcal_1)
1607 }
1608 \]
1609 and
1610 \[
1611 \boxed{
1612 \cost(\Lcal_2)
1613 \geq
1614 \cost(\Lcal_1).
1615 }
1616 \]
1617 \end{proposition}
1618 
1619 \begin{proof}\phantomsection\label{main:proof:045}
1620 Optimality gives
1621 \[
1622 \cost(\Lcal_1)+\lambda_1\err(\Lcal_1)
1623 \leq
1624 \cost(\Lcal_2)+\lambda_1\err(\Lcal_2)
1625 \]
1626 and
1627 \[
1628 \cost(\Lcal_2)+\lambda_2\err(\Lcal_2)
1629 \leq
1630 \cost(\Lcal_1)+\lambda_2\err(\Lcal_1).
1631 \]
1632 Adding the inequalities yields
1633 \[
1634 (\lambda_2-\lambda_1)
1635 \bigl(
1636 \err(\Lcal_2)-\err(\Lcal_1)
1637 \bigr)
1638 \leq0.
1639 \]
1640 The error ordering follows. Substitution into the first inequality gives the cost ordering.
1641 \phantomsection\label{main:proof:045:end}\end{proof}
1642 
1643 \begin{corollary}[Supported Pareto points]\label{main:corollary:014}\label{cor:supported_pareto}
1644 For \(\lambda>0\), every Lagrangian optimum is Pareto optimal for the pair
1645 \[
1646 (\cost,\err).
1647 \]
1648 The quantized budget compiler may recover Pareto points that are not optimal for any Lagrangian parameter.
1649 \end{corollary}
1650 
1651 \begin{proof}\phantomsection\label{main:proof:046}
1652 If another pruning had no larger cost and no larger error, with one strict inequality, its Lagrangian objective would be smaller. This contradicts optimality. A finite nonconvex frontier can contain points that are not supported by a line of slope \(-\lambda\), so the budget formulation is genuinely additional.
1653 \phantomsection\label{main:proof:046:end}\end{proof}
1654 
1655 \subsection{An exact error budget compiler}\label{main:subsection:014}
1656 
1657 \begin{theorem}[Quantized budget dynamic program]\label{main:theorem:024}\label{thm:budget_compiler}
1658 Assume every keep error is a nonnegative integer \(w_v\). Let \(B\in\NN\). For each node \(v\), define \(F_v(b)\) as the minimum cost of a pruning of the subtree at \(v\) with total error exactly \(b\). The keep option contributes
1659 \[
1660 F_v(w_v)\leq\kappa_v.
1661 \]
1662 The split option is obtained by min plus convolution of the child tables, followed by addition of \(\sigma_v\). The smallest value among
1663 \[
1664 F_{\rm root}(b),
1665 \qquad
1666 0\leq b\leq B,
1667 \]
1668 is the exact minimum cost among prunings whose quantized error is at most \(B\).
1669 \end{theorem}
1670 
1671 \begin{proof}\phantomsection\label{main:proof:047}
1672 The keep case is immediate. In the split case, an admissible pruning is the union of one admissible pruning from each child. Its error is the sum of the child errors and its cost is the sum of their costs plus \(\sigma_v\). Repeated min plus convolution enumerates every possible distribution of the integer budget across the child subtrees. Induction on the tree proves exactness.
1673 \phantomsection\label{main:proof:047:end}\end{proof}
1674 
1675 \begin{corollary}[Pseudo polynomial complexity]\label{main:corollary:015}\label{cor:budget_complexity}
1676 If the maximum number of children is bounded and the budget is \(B\), the quantized compiler has pseudo polynomial complexity in the number of tree nodes and \(B\). A direct pairwise convolution implementation has quadratic dependence on \(B\) at each merge.
1677 \end{corollary}
1678 
1679 \begin{proof}\phantomsection\label{main:proof:048}
1680 Every table has at most \(B+1\) entries. Merging two tables by direct convolution examines at most \((B+1)^2\) pairs. A bounded number of merges occurs at every node.
1681 \phantomsection\label{main:proof:048:end}\end{proof}
1682 
1683 \begin{figure}[t]
1684 \centering
1685 \includegraphics[width=.75\textwidth]{12_budget_compiler_frontier.pdf}
1686 \caption{The quantized dynamic program produces the exact cost and error frontier for its finite model.}
1687 \label{fig:budget_frontier}
1688 \end{figure}
1689 
1690 \subsection{Decision stability and local updates}\label{main:subsection:015}
1691 
1692 \begin{definition}[Decision gap]\label{main:definition:013}\label{def:decision_gap}
1693 At an internal node \(v\), let
1694 \[
1695 g_v
1696 =
1697 \abs{K_\lambda(v)-S_\lambda(v)}.
1698 \]
1699 This is the gap between the keep and split scores after the child subproblems have been solved.
1700 \end{definition}
1701 
1702 \begin{theorem}[Compiler stability]\label{main:theorem:025}\label{thm:compiler_stability}
1703 Suppose the solved keep score and split score at every internal node change in absolute value by less than $g_v/2$. These are scores after child subproblems have been solved, rather than unconstrained changes in raw local costs. Then every strict keep or split decision in the original optimal pruning remains unchanged.
1704 \end{theorem}
1705 
1706 \begin{proof}\phantomsection\label{main:proof:049}
1707 If the original keep score is lower by \(g_v\), then after perturbation it can increase by less than \(g_v/2\), while the split score can decrease by less than \(g_v/2\). The order remains strict. The split case is the same. Apply the argument from the leaves upward.
1708 \phantomsection\label{main:proof:049:end}\end{proof}
1709 
1710 \begin{theorem}[Local recompilation]\label{main:theorem:026}\label{thm:local_recompilation}
1711 Suppose local costs change only inside the subtree rooted at $u$, while ancestor and sibling local models are unchanged. Use correct cached sibling values from the old solve. Then every dynamic programming value outside that subtree and outside the ancestor chain from \(u\) to the root remains unchanged. Recomputing the subtree at \(u\) and that ancestor chain gives the exact new global optimum.
1712 \end{theorem}
1713 
1714 \begin{proof}\phantomsection\label{main:proof:050}
1715 A node value depends only on its own local data and the values of its children. Nodes outside the changed subtree and outside the ancestor chain have the same local data and the same child values. Their values remain unchanged. Every affected dependency lies inside the changed subtree or on the unique path from \(u\) to the root.
1716 \phantomsection\label{main:proof:050:end}\end{proof}
1717 
1718 \begin{remark}[Boundary of the compiler theorem]\label{main:remark:004}\label{rem:compiler_boundary}
1719 The exact compiler requires a finite candidate tree and additive local models. Arbitrary overlapping decompositions, nonlocal cache interactions, global rank coupling, and unrestricted structure discovery can destroy optimal substructure. No universal approximation factor is proved for those settings.
1720 \end{remark}
1721 
1722 \section{Device scheduling and amortization}\label{main:section:010}\label{sec:scheduling}
1723 
1724 \subsection{Shape classes}\label{main:subsection:016}
1725 
1726 \begin{definition}[Shape class]\label{main:definition:014}\label{def:shape_class}
1727 Two active packet blocks lie in the same shape class when their input dimension, output dimension, right hand side count, scalar type, accumulation type, and kernel contract agree.
1728 \end{definition}
1729 
1730 \begin{proposition}[Launch reduction by shape]\label{main:proposition:010}\label{prop:shape_launch}
1731 Suppose \(B\) independent active blocks fall into \(g\) shape classes, and suppose the backend accepts one batched array or pointer list for each class. Grouped execution then requires at most \(g\) launches instead of \(B\), without changing the arithmetic or the output.
1732 \end{proposition}
1733 
1734 \begin{proof}\phantomsection\label{main:proof:051}
1735 Place all blocks in one class along a batch dimension and invoke one batched kernel for that class. Every block receives the same local computation. One launch per class is sufficient.
1736 \phantomsection\label{main:proof:051:end}\end{proof}
1737 
1738 CUDA provides graph capture and replay for repeated dependency graphs \cite{NVIDIA2026}. These mechanisms reduce repeated launch work, but they do not change the operator or its arithmetic count.
1739 
1740 \subsection{Critical path}\label{main:subsection:017}
1741 
1742 \begin{theorem}[Acyclic schedule depth]\label{main:theorem:027}\label{thm:dag_depth}
1743 Let a packet task graph be acyclic. Assume every task takes one normalized stage, every edge is a dependency, and an unlimited number of independent tasks may run concurrently. For a nonempty finite graph, the minimum possible number of sequential stages equals the number of vertices on a longest directed path. The empty task graph has depth zero.
1744 \end{theorem}
1745 
1746 \begin{proof}\phantomsection\label{main:proof:052}
1747 Every path gives a lower bound because its tasks must execute in order. The topological level schedule in \cref{thm:level_schedule} achieves exactly this depth. With no tasks, the empty schedule has depth zero.
1748 \phantomsection\label{main:proof:052:end}\end{proof}
1749 
1750 \subsection{One time costs}\label{main:subsection:018}
1751 
1752 \begin{theorem}[Amortization threshold]\label{main:theorem:028}\label{thm:amortization}
1753 Let \(s\geq0\) be a one time setup cost. Let $t_b,t_c\ge0$ be the baseline and candidate times per use and let $R\in\NN$ be a reuse count. If \(t_b>t_c\), the candidate is faster after \(R\) uses exactly when
1754 \[
1755 R>
1756 \frac{s}{t_b-t_c}.
1757 \]
1758 The smallest integer reuse count is
1759 \[
1760 \boxed{
1761 R_{\min}
1762 =
1763 \left\lfloor
1764 \frac{s}{t_b-t_c}
1765 \right\rfloor+1.
1766 }
1767 \]
1768 If \(t_b\leq t_c\), no positive reuse count amortizes a positive setup cost.
1769 \end{theorem}
1770 
1771 \begin{proof}\phantomsection\label{main:proof:053}
1772 Compare
1773 \[
1774 s+Rt_c
1775 \]
1776 with
1777 \[
1778 Rt_b.
1779 \]
1780 Rearranging gives the stated inequality.
1781 \phantomsection\label{main:proof:053:end}\end{proof}
1782 
1783 \begin{corollary}[Graph replay threshold]\label{main:corollary:016}\label{cor:graph_threshold}
1784 If graph capture costs $c\ge0$, ordinary execution costs $t_o\ge0$, and replay costs $0\le t_r<t_o$, then capture is profitable when
1785 \[
1786 R>
1787 \frac{c}{t_o-t_r}.
1788 \]
1789 \end{corollary}
1790 
1791 \begin{proof}\phantomsection\label{main:proof:054}
1792 Apply \cref{thm:amortization} with setup cost \(c\).
1793 \phantomsection\label{main:proof:054:end}\end{proof}
1794 
1795 \begin{corollary}[Materialization threshold]\label{main:corollary:017}\label{cor:materialization_threshold}
1796 Suppose a compressed representation costs $t_c\ge0$ per downstream use, an explicit representation costs $t_e\ge0$ per use, and one time materialization costs $m\ge0$, with reuse count $R\in\NN$. Materialization is profitable exactly when
1797 \[
1798 R>
1799 \frac{m}{t_c-t_e}
1800 \]
1801 and \(t_c>t_e\).
1802 \end{corollary}
1803 
1804 \begin{proof}\phantomsection\label{main:proof:055}
1805 Compare repeated compressed use with materialization followed by repeated explicit use.
1806 \phantomsection\label{main:proof:055:end}\end{proof}
1807 
1808 \section{Benchmark contracts and measurement protocol}\label{main:section:011}\label{sec:benchmark_protocol}
1809 
1810 Performance data are meaningful only after the mathematical and timing contracts are fixed. The protocol follows standard concerns from numerical linear algebra, floating point analysis, and machine performance modeling \cite{Higham2002,GolubVanLoan2013,WilliamsWatermanPatterson2009}.
1811 
1812 \begin{definition}[Case signature]\label{main:definition:015}\label{def:case_signature}
1813 A case signature records the dimensions, packet layout, rank, right hand side count, random seed, source values, source type, reference computation, and output contract.
1814 \end{definition}
1815 
1816 \begin{definition}[Fair match]\label{main:definition:016}\label{def:fair_match}
1817 Two measured rows are a fair match when they share one case signature, return the same declared output, use the same correctness reference, and have identical timing scopes. If one row includes an extra conversion or setup cost, it is a different scope and must be reported separately or added to both rows before comparison.
1818 \end{definition}
1819 
1820 \begin{proposition}[Fair matching is an equivalence relation]\label{main:proposition:011}\label{prop:fair_equivalence}
1821 On rows with complete metadata, equality of the case signature, output contract, correctness reference, and timing scope is an equivalence relation.
1822 \end{proposition}
1823 
1824 \begin{proof}\phantomsection\label{main:proof:056}
1825 Equality of each metadata field is reflexive, symmetric, and transitive. Their conjunction has the same properties.
1826 \phantomsection\label{main:proof:056:end}\end{proof}
1827 
1828 \begin{definition}[Representation gain]\label{main:definition:017}\label{def:representation_gain}
1829 Representation gain compares dense and structured work or storage for the same operator, independently of the processor.
1830 \end{definition}
1831 
1832 \begin{definition}[Device gain]\label{main:definition:018}\label{def:device_gain}
1833 Device gain compares processors for the same representation, output, and timing scope.
1834 \end{definition}
1835 
1836 \begin{definition}[Kernel gain]\label{main:definition:019}\label{def:kernel_gain}
1837 Kernel gain compares two implementations on the same device for the same representation, output, and timing scope.
1838 \end{definition}
1839 
1840 \begin{definition}[Amortization gain]\label{main:definition:020}\label{def:amortization_gain}
1841 Amortization gain records the effect of one time discovery, compilation, transfer, capture, or materialization at a declared reuse count.
1842 \end{definition}
1843 
1844 A measured total speedup can mix these four effects. A multiplicative decomposition is only exact when matched rows isolate the factors. Otherwise the decomposition is interpretive.
1845 
1846 \begin{figure}[t]
1847 \centering
1848 \includegraphics[width=.86\textwidth]{22_benchmark_contracts.pdf}
1849 \caption{The corpus records timing scope and output contract as separate fields.}
1850 \label{fig:contracts}
1851 \end{figure}
1852 
1853 \begin{remark}[Logical work rate]\label{main:remark:005}\label{rem:logical_gflops}
1854 A structured method may report a dense equivalent logical work rate greater than the physical dense throughput because it does not execute the dense operations. This statistic describes avoided work. It is not a claim that the device performed those operations.
1855 \end{remark}
1856 
1857 \section{Recorded CPU and GPU corpus}\label{main:section:012}\label{sec:measured}
1858 
1859 The empirical corpus contains 154 measured rows from sixteen suites. It has 74 CPU rows and 80 GPU rows. The GPU was an NVIDIA GeForce RTX 3070 Ti Laptop GPU with compute capability 8.6 and 8 GB nominal memory. The recorded software used PyTorch 2.11.0 with a CUDA 12.8 build. These facts identify the environment. They do not predict another machine.
1860 
1861 \begin{table}[t]
1862 \centering
1863 \small
1864 \begin{tabular}{>{\raggedright\arraybackslash}p{.24\textwidth}r>{\raggedright\arraybackslash}p{.24\textwidth}r}
1865 \toprule
1866 Suite & Rows & Suite & Rows \\
1867 \midrule
1868 Dense multiplication & 28 & BDLR application & 15 \\
1869 Block diagonal & 12 & BDLR product & 12 \\
1870 Circulant & 12 & Toeplitz & 12 \\
1871 Kronecker & 12 & Block banded & 8 \\
1872 Sparse plus low rank & 8 & HODLR & 8 \\
1873 Strassen & 6 & Transfer pipeline & 5 \\
1874 Walsh Hadamard & 4 & Hidden blocks & 4 \\
1875 Randomized low rank & 4 & Iterative refinement & 4 \\
1876 \bottomrule
1877 \end{tabular}
1878 \caption{Coverage of the sanitized measured corpus.}
1879 \label{tab:suites}
1880 \end{table}
1881 
1882 The 154 archived rows contain recorded case identifiers, seeds, types, output contracts, and timing scopes. Source-value hashes and explicit correctness-reference fields are absent. Setup time, discovery time, and reuse count are nonnull in only two rows each. Correctness errors are present in 149 rows. Equality of available metadata cannot establish the full fair-match relation in \cref{def:fair_match}. In particular, CPU-resident and GPU-device-resident scopes differ, and graph replay excludes capture while a device-resident baseline has a different raw scope. The two-stream example also compares wall time with a serial device-time scope. The values below are retained as historical observations with these limitations.
1883 
1884 \subsection{Dense precision frontier}\label{main:subsection:019}
1885 
1886 At dimension \(4096\), the device resident dense rows gave the values in \cref{tab:dense_frontier}.
1887 
1888 \begin{table}[t]
1889 \centering
1890 \small
1891 \begin{tabular}{lrrr}
1892 \toprule
1893 Method & Median time in ms & Logical GFLOP per second & Relative error \\
1894 \midrule
1895 FP64 & 541.2803 & 253.91 & \(2.18\times10^{-15}\) \\
1896 FP32 & 18.1309 & 7580.35 & \(5.76\times10^{-7}\) \\
1897 TF32 & 8.6999 & 15797.75 & \(2.94\times10^{-4}\) \\
1898 FP16 & 4.0264 & 34134.72 & \(3.60\times10^{-4}\) \\
1899 BF16 & 4.1482 & 33132.00 & \(2.88\times10^{-3}\) \\
1900 \bottomrule
1901 \end{tabular}
1902 \caption{Recorded dense accuracy and throughput at one dimension.}
1903 \label{tab:dense_frontier}
1904 \end{table}
1905 
1906 \begin{figure}[t]
1907 \centering
1908 \includegraphics[width=.76\textwidth]{13_dense_gpu_throughput.pdf}
1909 \caption{Recorded device resident dense throughput across dimensions and scalar types.}
1910 \label{fig:dense_throughput}
1911 \end{figure}
1912 
1913 \begin{figure}[t]
1914 \centering
1915 \includegraphics[width=.75\textwidth]{14_accuracy_frontier.pdf}
1916 \caption{The speed and error tradeoff depends on both precision and problem size.}
1917 \label{fig:accuracy_frontier}
1918 \end{figure}
1919 
1920 The frontier is not a universal ordering of scalar types. It records one generated family. Packet precision assignment should use local sensitivities and supported kernels rather than one global type.
1921 
1922 \subsection{Regular packet operators}\label{main:subsection:020}
1923 
1924 For an \(8192\) by \(8192\) equal block operator with block width \(64\) and \(32\) right hand sides, the batched packet calculation took
1925 \[
1926 0.046032\text{ ms}.
1927 \]
1928 It was
1929 \[
1930 14.8919
1931 \]
1932 times faster than the recorded device resident dense GPU row with the same recorded case, output, and timing scope and
1933 \[
1934 235.7404
1935 \]
1936 times faster than the recorded CPU dense row (resident CPU scope versus device-resident GPU scope). The relative error was
1937 \[
1938 8.17\times10^{-8}.
1939 \]
1940 
1941 \begin{figure}[t]
1942 \centering
1943 \includegraphics[width=.74\textwidth]{15_block_diagonal_speedup.pdf}
1944 \caption{Regular equal shapes expose both representation reduction and batched execution.}
1945 \label{fig:block_speedup}
1946 \end{figure}
1947 
1948 For BDLR application with the same dimension and right hand side count, block width \(64\), and rank \(8\), CUDA Graph replay took
1949 \[
1950 0.079872\text{ ms}.
1951 \]
1952 It was
1953 \[
1954 10.9355
1955 \]
1956 times faster than the recorded dense GPU row (graph-replay scope versus device-resident scope) and
1957 \[
1958 153.6527
1959 \]
1960 times faster than the recorded CPU dense row (resident CPU scope versus device-resident GPU scope). Capture time was excluded from the replay scope and must be amortized separately.
1961 
1962 \begin{figure}[t]
1963 \centering
1964 \includegraphics[width=.77\textwidth]{16_bdlr_timings.pdf}
1965 \caption{Representation and launch strategy both affect BDLR application time.}
1966 \label{fig:bdlr_timings}
1967 \end{figure}
1968 
1969 \subsection{Compressed products}\label{main:subsection:021}
1970 
1971 A compressed BDLR product returns block and low rank factors. Materialization returns all matrix entries. The historical product rows illustrate the cost of changing the output contract. Compressed and materialized outputs are separate computations. The GPU materialized BDLR rows at dimensions 512 and 1024 are sums of separately measured factor and expansion medians, as the raw notes state. They are not timing medians for complete workflows.
1972 
1973 \begin{figure}[t]
1974 \centering
1975 \includegraphics[width=.82\textwidth]{17_output_contract_timings.pdf}
1976 \caption{Compressed and explicit products answer different computational questions.}
1977 \label{fig:output_contracts}
1978 \end{figure}
1979 
1980 The algebra in \cref{thm:bdlr_product} explains why the compressed product can remain small. The materialization cost in \cref{prop:materialization_cost} explains why explicit output can erase the advantage.
1981 
1982 \subsection{Discovery and transfer}\label{main:subsection:022}
1983 
1984 For hidden permuted blocks at dimension \(512\), structure discovery took \(2.1102\) ms, dense application took \(0.0969\) ms, and packet application took \(0.0151\) ms. The exact threshold in \cref{thm:amortization} is \(26\) uses. At dimension \(1024\), the corresponding values were \(5.2726\), \(0.1896\), and \(0.0241\) ms, giving \(32\) uses.
1985 
1986 \begin{figure}[t]
1987 \centering
1988 \includegraphics[width=.75\textwidth]{18_discovery_amortization.pdf}
1989 \caption{A profitable repeated representation can remain unprofitable for one use.}
1990 \label{fig:discovery}
1991 \end{figure}
1992 
1993 For one \(2048\) dense FP32 product, device resident graph replay took \(1.9993\) ms. The end to end pageable path took \(7.3036\) ms. The pinned path took \(6.3084\) ms. Two serial device resident products took \(3.6557\) ms, while the tested two stream path took \(4.2368\) ms.
1994 
1995 \begin{figure}[t]
1996 \centering
1997 \includegraphics[width=.78\textwidth]{19_transfer_pipeline.pdf}
1998 \caption{Device residence, transfer, graph replay, and a multioutput pipeline have different timing scopes.}
1999 \label{fig:transfer_pipeline}
2000 \end{figure}
2001 
2002 \subsection{Negative controls}\label{main:subsection:023}
2003 
2004 The corpus retains methods that lost in the recorded baseline comparisons. These are descriptive ratios. Scope and complete-metadata qualifications apply.
2005 
2006 \begin{itemize}
2007 \item At dimension \(512\), recursive Strassen multiplication with leaf size \(128\) took \(7.0950\) ms, while the recorded OpenBLAS dense multiplication took \(1.9597\) ms.
2008 \item At dimension \(2048\), the tested recursive HODLR GPU path took \(5.6079\) ms, while dense GPU multiplication took \(0.0860\) ms.
2009 \item The tested FP32 iterative refinement solve took \(22.0848\) ms, while the direct FP64 GPU solve took \(12.6286\) ms with a comparable final error.
2010 \item The two stream dense path was slower than the serial two product path.
2011 \item The GPU Walsh--Hadamard transform was slower than the CPU implementation at length \(65536\), then became faster at length \(262144\).
2012 \end{itemize}
2013 
2014 Fast multiplication can be asymptotically valuable and still lose at measured sizes because of recursion, memory traffic, and kernel granularity. Hierarchical matrix methods have similar crossover issues \cite{AmbikasaranDarve2013}. Transform algorithms and tensor factorizations form additional structure classes \cite{CooleyTukey1965,Oseledets2011}.
2015 
2016 \begin{figure}[t]
2017 \centering
2018 \includegraphics[width=.77\textwidth]{20_negative_controls.pdf}
2019 \caption{Values above one are slowdowns against the matched baseline.}
2020 \label{fig:negative_controls}
2021 \end{figure}
2022 
2023 \begin{figure}[t]
2024 \centering
2025 \includegraphics[width=.78\textwidth]{21_best_gpu_speedups.pdf}
2026 \caption{The best recorded gain differs strongly across suites. Regular packet structure dominates the largest measured rows.}
2027 \label{fig:best_speedups}
2028 \end{figure}
2029 
2030 \section{Implementation and verification}\label{main:section:013}\label{sec:implementation}
2031 
2032 The source corpus contains a small reference implementation of the exact operators and optimization algorithms used in the paper. The mathematical interface is deliberately narrow. This Python implementation is a numerical prototype, not extracted Lean code. Its stored-key support can retain zero-valued blocks, whereas the exact support graph uses nonzero blocks. Its reduced numerical QR is not the exact rank-compact range construction proved in Appendix F. Prototype tests and theorem proofs therefore have different targets.
2033 
2034 \begin{enumerate}
2035 \item A packet layout records dimensions and coordinate slices.
2036 \item A packet operator stores active blocks and evaluates them without materializing inactive blocks.
2037 \item BDLR factors support exact application, product, solve, and small core recompression.
2038 \item Graph utilities compute support products, path closure, Schur fill inclusion, and topological levels.
2039 \item Tree dynamic programs solve the Lagrangian and quantized budget models.
2040 \item A multiple choice dynamic program assigns discrete precision options.
2041 \item Benchmark utilities audit case identity, output contracts, timing scopes, and retained negative controls.
2042 \end{enumerate}
2043 
2044 The accompanying verification record reports a deterministic verifier with $37392$ checks in nine groups. It covers packet evaluation, adjoints, composition, structured primitives, recompression optimality, product chain error, graph support, Schur fill, tree compilation, precision assignment, and the benchmark corpus. It also reports eleven passing unit tests and twenty-three passing original-source tests for BDLR operations, Krylov methods, mixed precision, measured data schema, and figure generation.
2045 
2046 \begin{table}[t]
2047 \centering
2048 \small
2049 \begin{tabular}{>{\raggedright\arraybackslash}p{.34\textwidth}>{\raggedright\arraybackslash}p{.55\textwidth}}
2050 \toprule
2051 Layer & Verification target \\
2052 \midrule
2053 Packet semantics & materialization, application, adjoint, composition \\
2054 Low rank algebra & exact product, Woodbury solve, optimal core recompression \\
2055 Graph algorithms & path support, topological levels, Schur fill inclusion \\
2056 Optimization & Lagrangian tree compiler, budget compiler, precision assignment \\
2057 Numerical bounds & product chain error, block perturbation, clipped allocation \\
2058 Empirical corpus & row count, suite count, device split, status, error ceiling \\
2059 Document & proofs, labels, citations, figures, font embedding, render audit \\
2060 \bottomrule
2061 \end{tabular}
2062 \caption{Verification layers for the supplied implementation.}
2063 \label{tab:verification_layers}
2064 \end{table}
2065 
2066 \section{Hardware implications and exact boundaries}\label{main:section:014}\label{sec:boundaries}
2067 
2068 The finite results suggest a device interface with grouped block kernels, low rank contractions, packet metadata, a transfer graph scheduler, and per edge precision tags. Such an interface may reduce materialization and repeated launch work when the operator family remains closed under the chosen representation.
2069 
2070 The exact results stop at a finite candidate tree and a declared device model. Three broader problems remain open.
2071 
2072 \begin{enumerate}
2073 \item A universal discovery algorithm would need to find a profitable packet structure from an arbitrary operator and data distribution. The simultaneous spectral theorem and randomized low rank approximation cover special cases only.
2074 \item A universal compiler would need a controlled approximation factor when costs are nonlocal and candidate decompositions overlap. The dynamic programs in \cref{thm:tree_compiler,thm:budget_compiler} do not cover that setting.
2075 \item A packet native memory hierarchy would need a communication theorem for a family of operators and a specified architecture. The one pass byte model in \cref{prop:streaming_bytes} is not such a lower bound.
2076 \end{enumerate}
2077 
2078 Communication lower bounds require explicit machine models \cite{HongKung1981,BallardDemmelHoltzSchwartz2011}. CUDA Graph behavior and kernel constraints depend on the installed CUDA platform \cite{NVIDIA2026}. The measured results in this paper do not prove an optimal accelerator.
2079 
2080 The article also makes no claim about the exponent of generic dense matrix multiplication. Packet algorithms answer a different question. They exploit a restricted operator class and may return a compressed output.
2081 
2082 \section{Conclusion}\label{main:section:015}\label{sec:conclusion}
2083 
2084 A packet operator representation records which block maps exist, how they are stored, and what output they denote. Exact block evaluation and composition turn the support graph into algebra. Arithmetic, storage, and one pass data movement follow from the active block dimensions. Block diagonal plus low rank operators remain closed under multiplication before recompression. Their inverse is reduced by Woodbury, and their correction can be recompressed optimally through a small singular value problem.
2085 
2086 Sparse transfer graphs predict path support, triangular schedule depth, and possible fill under elimination. Packet diagonal structure gives local Krylov spaces and a certified block Jacobi condition bound. Packetwise perturbations give a direct mixed precision certificate. Continuous tolerances follow from a convex allocation, while finite data type choices lead to an exact multiple choice budget program.
2087 
2088 A finite candidate refinement tree supports a genuine packet compiler. The Lagrangian optimum, the quantized error budget optimum, and local recompilation are exact under additive hypotheses. This is narrower than unrestricted automatic discovery, but it turns one broad compiler conjecture into a proved algorithmic statement.
2089 
2090 The measured corpus shows the practical distinction between regular and irregular structure. Equal block and BDLR operators gained substantially on the recorded GPU. Recursive HODLR, Strassen, iterative refinement, a small transform, and an attempted overlap lost to recorded baselines. Discovery, transfer, capture, and materialization changed the break even point.
2091 
2092 The computational principle is therefore precise. Remove work through a verified representation. Preserve the declared output. Compile only inside a model whose costs and errors are stated. Assign precision by sensitivity. Amortize one time costs. Keep negative controls. Treat the final timing as a property of the representation, device, kernel, and reuse regime together.
2093 
2094 

Back to source controls

Back to Realization Calculus