COS Method

Download original file

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

Read the typeset paper

Read this chapter with formatted equations

1 \chapter{Exact streamed representations}\label{ch:streaming}
2 \chapterlead{Large factors need not be resident all at once. The identities in this
3 chapter separate total data size from the working set without changing the
4 operator. A storage partition is therefore a choice about execution, not a
5 statistical approximation. Weighted features provide a concrete construction,
6 but every subsequent identity applies to a supplied real factor.}
7 \section{Weighted-feature specialization}\label{sec:final-7}
8 
9 A prominent application constructs the general factor from weighted,
10 covariate-residualized features.  The sequence kernel association test (SKAT) provides one application
11 \cite{Wu2011}, but the same formulas apply to any feature matrix with
12 nonnegative column weights.
13 
14 Let
15 \[
16 G\in\mathbb R^{n\times m},
17 \qquad
18 W=\operatorname{diag}(w_1,\ldots,w_m),
19 \quad w_j\ge 0,
20 \]
21 and let \(Q_X\in\mathbb R^{n\times p}\) have orthonormal columns spanning
22 the null-model covariates.  Define
23 \[
24 M_X=I-Q_XQ_X^\top,
25 \qquad
26 A=M_XGW^{1/2},
27 \qquad
28 K=AA^\top.
29 \]
30 
31 The central implementation difficulty is not the algebraic size of the
32 model but the mismatch between the total factor size \(n\times m\) and the
33 working memory of an accelerator.  The key point of this book is that
34 several decompositions act on different mathematical coordinates and therefore
35 compose:
36 \[
37 \boxed{
38 \text{storage blocks}
39 \circ
40 \text{structural phases}
41 \circ
42 \text{Krylov/range compression}
43 \circ
44 \text{spectral tail packets}.
45 }
46 \]
47 
48 Storage blocking controls residence.  Structural phases encode linkage-disequilibrium (LD) windows,
49 annotations, frequency bands, or other coordinate groups.  Krylov methods
50 compress the operator to the rank actually requested.  Spectral packets
51 compress the final generalized-\(\chi^2\) transform.
52 
53 
54 \begin{figure}[htbp]
55 \centering
56 \includegraphics[width=\linewidth]{figures_production/figure_01_residualization_geometry.pdf}
57 \caption{Orthogonal residualization separates a feature into a component in the covariate span and a component retained by $M_X$. The diagram is an illustrative coordinate model. The identity behind it is the exact projector decomposition.}
58 \label{fig:residualization-geometry}
59 \end{figure}
60 
61 
62 \section{Variant-block Gram decomposition}\label{sec:final-8}
63 
64 Partition the variant indices into disjoint blocks
65 \[
66 \{1,\ldots,m\}
67 =
68 J_1\sqcup\cdots\sqcup J_B.
69 \]
70 Let
71 \[
72 G_b=G_{[:,J_b]},
73 \qquad
74 W_b=W_{J_b,J_b},
75 \qquad
76 A_b=M_XG_bW_b^{1/2}.
77 \]
78 Then
79 \[
80 A=[A_1\ \cdots\ A_B].
81 \]
82 
83 \begin{theorem}[Additive block Gram identity]\label{claim:final-8-1}
84 \[
85 \boxed{
86 K
87 =
88 AA^\top
89 =
90 \sum_{b=1}^B A_bA_b^\top
91 =
92 M_X
93 \left(
94 \sum_{b=1}^B G_bW_bG_b^\top
95 \right)
96 M_X.
97 }
98 \]
99 \end{theorem}
100 
101 \begin{proof}
102 Horizontal concatenation gives
103 \[
104 [A_1\ \cdots\ A_B]
105 \begin{bmatrix}
106 A_1^\top\\
107 \vdots\\
108 A_B^\top
109 \end{bmatrix}
110 =
111 \sum_bA_bA_b^\top.
112 \]
113 The final form follows from
114 \(A_b=M_XG_bW_b^{1/2}\) and
115 \(M_X^\top=M_X=M_X^2\).
116 \end{proof}
117 
118 \begin{remark}\label{claim:final-8-2}
119 The theorem is exact for every partition.  Block boundaries are execution
120 parameters, not statistical parameters.
121 \end{remark}
122 
123 
124 \begin{figure}[htbp]
125 \centering
126 \includegraphics[width=\linewidth]{figures_production/figure_02_exact_block_kernel_assembly.pdf}
127 \caption{Each resident column block contributes a complete sample-space Gram summand $A_bA_b^{\mathsf T}$. Adding the summands recovers $AA^{\mathsf T}$ exactly. Storage partitioning changes residence and evaluation order, not the represented operator.}
128 \label{fig:block-kernel-assembly}
129 \end{figure}
130 
131 
132 \section{Exact streamed kernel actions}\label{sec:final-9}
133 
134 \begin{theorem}[One-pass sample-space streaming action]\label{claim:final-9-1}
135 For every block vector \(V\in\mathbb R^{n\times s}\),
136 \[
137 \boxed{
138 KV
139 =
140 M_X
141 \sum_{b=1}^B
142 G_bW_bG_b^\top
143 M_XV.
144 }
145 \]
146 \end{theorem}
147 
148 \begin{proof}
149 Apply the additive block Gram identity to \(KV\), then associate each block
150 product as
151 \[
152 G_bW_b(G_b^\top M_XV).
153 \]
154 \end{proof}
155 
156 \begin{algorithm}[Streamed \(KV\)]\label{claim:final-9-2}
157 \begin{enumerate}[label=\arabic*.]
158 \item Compute \(U=M_XV\).
159 \item Set \(Y=0\).
160 \item For each resident block \(b\), compute
161 \[
162 T_b=G_b^\top U,
163 \qquad
164 Y\leftarrow Y+G_b(W_bT_b).
165 \]
166 \item Return \(M_XY\).
167 \end{enumerate}
168 \end{algorithm}
169 
170 \begin{corollary}[Working memory independent of total variant count]\label{claim:final-9-3}
171 If at most \(c\) variants are resident and the block-vector width is \(s\),
172 then the algebraic working set is
173 \[
174 \boxed{
175 O(nc+ns+cs+np),
176 }
177 \]
178 independent of total variant count \(m\).
179 \end{corollary}
180 
181 
182 \section{Streamed factor actions and Golub--Kahan}\label{sec:final-10}
183 
184 Let \(V\in\mathbb R^{m\times s}\) be partitioned conformally:
185 \[
186 V=
187 \begin{bmatrix}
188 V_1\\
189 \vdots\\
190 V_B
191 \end{bmatrix}.
192 \]
193 
194 \begin{theorem}[Streamed \(AV\)]\label{claim:final-10-1}
195 \[
196 \boxed{
197 AV
198 =
199 M_X
200 \sum_{b=1}^B
201 G_bW_b^{1/2}V_b.
202 }
203 \]
204 \end{theorem}
205 
206 \begin{theorem}[Streamed \(A^\top U\)]\label{claim:final-10-2}
207 For \(U\in\mathbb R^{n\times s}\),
208 \[
209 \boxed{
210 (A^\top U)_b
211 =
212 W_b^{1/2}G_b^\top M_XU.
213 }
214 \]
215 \end{theorem}
216 
217 \begin{proof}
218 Both identities follow by substituting
219 \(A=M_XGW^{1/2}\) and partitioning the variant coordinate.
220 \end{proof}
221 
222 \begin{corollary}[Streaming Golub--Kahan equivalence]\label{claim:final-10-3}
223 With the same initial state, algorithmic choices, and random draws, any
224 exact-arithmetic algorithm whose access to \(A\) occurs only through
225 \(AV\) and \(A^\top U\) produces the same iterates when those products are
226 implemented by the preceding streamed formulas.
227 \end{corollary}
228 
229 This applies directly to block Golub--Kahan bidiagonalization, matrix-free
230 singular-value iteration, range finding, and restarted variants.  The
231 factor need not be resident.
232 
233 
234 \section{Two-pass dual action}\label{sec:final-11}
235 
236 The variant-space normal operator is
237 \[
238 H=A^\top A.
239 \]
240 
241 \begin{theorem}[Two-pass streamed dual action]\label{claim:final-11-1}
242 For every \(V\in\mathbb R^{m\times s}\),
243 \[
244 \boxed{
245 HV=A^\top(AV).
246 }
247 \]
248 The product is evaluated exactly by a streamed pass forming
249 \(Y=AV\in\mathbb R^{n\times s}\), followed by a second streamed pass
250 emitting the blocks
251 \[
252 (HV)_b=W_b^{1/2}G_b^\top M_XY.
253 \]
254 \end{theorem}
255 
256 This form is particularly natural for variant-space polyphase Krylov
257 iterations because the right-coordinate blocks can be produced and consumed
258 sequentially.
259 
260 
261 \section{Two-dimensional sample/variant tiling}\label{sec:final-12}
262 
263 A variant block may itself exceed device memory.  Let
264 \[
265 H_b=G_bW_b^{1/2}
266 \]
267 and partition samples into row tiles \(I_1,\ldots,I_R\).
268 
269 \begin{theorem}[Two-pass residualized tile construction]\label{claim:final-12-1}
270 Define
271 \[
272 R_b=Q_X^\top H_b.
273 \]
274 Then
275 \[
276 \boxed{
277 R_b
278 =
279 \sum_r
280 Q_{X,I_r}^\top H_{b,I_r},
281 }
282 \]
283 and after \(R_b\) is known,
284 \[
285 \boxed{
286 A_{b,I_r}
287 =
288 H_{b,I_r}-Q_{X,I_r}R_b.
289 }
290 \]
291 Consequently,
292 \[
293 \boxed{
294 K_{I_r,I_s}
295 =
296 \sum_b
297 A_{b,I_r}A_{b,I_s}^\top.
298 }
299 \]
300 \end{theorem}
301 
302 \begin{proof}
303 The first identity is the row-partitioned matrix product
304 \(Q_X^\top H_b\).  The second is the definition
305 \(M_XH_b=H_b-Q_XQ_X^\top H_b\).  The last identity is the corresponding
306 sample-tile expansion of \(AA^\top\).
307 \end{proof}
308 
309 With sample tile height \(t\) and resident variant width \(c\), the dominant
310 working storage, including a resident covariate tile, is
311 \[
312 O(tc+pc+tp+t^2).
313 \]
314 This count excludes the backing input store, completed output blocks,
315 allocator overhead, and library workspaces.
316 
317 
318 \section{Partition refinement and permutation invariance}\label{sec:final-13}
319 
320 \begin{theorem}[Partition-refinement invariance]\label{claim:final-13-1}
321 Let \(\mathcal J\) be any partition of the variants and let
322 \(\mathcal J'\) refine \(\mathcal J\).  Then
323 \[
324 \boxed{
325 \sum_{J\in\mathcal J}A_JA_J^\top
326 =
327 \sum_{J'\in\mathcal J'}A_{J'}A_{J'}^\top
328 =
329 AA^\top.
330 }
331 \]
332 The analogous streamed operator action is also invariant.
333 \end{theorem}
334 
335 \begin{proof}
336 Each coarse block is a disjoint union of refined blocks.  Apply the additive
337 Gram theorem within each coarse block and sum.
338 \end{proof}
339 
340 \begin{corollary}[Adaptive memory splitting]\label{claim:final-13-2}
341 A runtime may recursively subdivide any block that exceeds available device
342 memory without changing the exact operator.
343 \end{corollary}
344 
345 \begin{proposition}[Variant permutation invariance]\label{claim:final-13-3}
346 If the columns of \(G\) and the diagonal entries of \(W\) are permuted by the
347 same permutation, \(K\) and every streamed sample-space action are unchanged.
348 \end{proposition}
349 
350 
351 \begin{figure}[htbp]
352 \centering
353 \includegraphics[width=\linewidth]{figures_production/figure_08_partition_invariance.pdf}
354 \caption{Gram and operator-action discrepancies for five retained partition controls. Horizontal connectors pair the two measurements for each partition. The small floating discrepancies support the implementation check, while the partition identity itself is exact algebra.}
355 \label{fig:partition-invariance}
356 \end{figure}
357 
358 
359 \section{Memory planning law}\label{sec:final-27}
360 
361 Let \(d\) be the number of bytes per floating scalar, \(c\) the resident
362 variant count, \(s\) the Krylov width, and \(p\) the covariate rank.
363 Ignoring library workspace and allocator overhead, a streamed sample-space
364 action has working memory on the order of
365 \[
366 \boxed{
367 d\,[nc+2ns+cs+np+O(ps+c)].
368 }
369 \]
370 
371 The full factor costs
372 \[
373 dnm.
374 \]
375 
376 For fixed \(n,p,s,c\),
377 \[
378 \boxed{
379 \frac{M_{\rm stream}}{M_{\rm full}}
380 =
381 O\!\left(
382 \frac{c}{m}
383 +
384 \frac{s}{m}
385 +
386 \frac{p}{m}
387 \right),
388 }
389 \]
390 so device memory remains bounded as total variant count grows.
391 
392 \begin{table}[htbp]
393 \centering
394 \caption{Algebraic float64 memory model for \(n=2000\), \(p=4\),
395 Krylov width \(s=16\), resident variant block \(c=1024\), and
396 sample tile height \(t=256\).}
397 \begin{tabular}{rrrr}
398 \toprule
399 \(m\) & full factor (MiB) & streamed \(KV\) (MiB) & tiled Gram (MiB)\\
400 \midrule
401 10,000 & 152.6 & 16.308 & 8.555\\
402 50,000 & 762.9 & 16.308 & 8.555\\
403 100,000 & 1,525.9 & 16.308 & 8.555\\
404 500,000 & 7,629.4 & 16.308 & 8.555\\
405 1,000,000 & 15,258.8 & 16.308 & 8.555\\
406 \bottomrule
407 \end{tabular}
408 \end{table}
409 
410 
411 \begin{figure}[htbp]
412 \centering
413 \includegraphics[width=\linewidth]{figures_production/figure_03_memory_law.pdf}
414 \caption{Binary64 payload as total factor width increases. The complete factor grows with $m$, while the streamed and tiled payloads remain constant at fixed resident dimensions. These are algebraic sizing data, not measurements of total device allocation.}
415 \label{fig:memory-law}
416 \end{figure}
417 
418 
419 
420 \section{The running example in three storage blocks}\label{sec:worked-stream}
421 Partition the factor of Section~\ref{sec:worked-law} into consecutive pairs of columns,
422 \[
423 A=[A_1\ A_2\ A_3],\qquad A_i\in\mathbb R^{4\times2}.
424 \]
425 Although $A_2=A_3$, both column groups must be counted. Their Gram contributions are
426 \[
427 A_1A_1^\top=\frac12
428 \begin{pmatrix}1&0&1&0\\0&1&0&1\\1&0&1&0\\0&1&0&1\end{pmatrix},\qquad
429 A_2A_2^\top=A_3A_3^\top=\frac12
430 \begin{pmatrix}1&0&-1&0\\0&1&0&-1\\-1&0&1&0\\0&-1&0&1\end{pmatrix}.
431 \]
432 The sum is $K$. Merging the last two blocks preserves it. Discarding $A_3$ instead changes the operator and its probability law.
433 
434 For the declared vector $v=(1,2,-1,0)^\top$, the action can be evaluated without forming any block Gram matrix:
435 \[
436 \sum_{j=1}^3 A_j(A_j^\top v)=Kv=(2,3,-2,-1)^\top.
437 \]
438 The symbolic check confirms both storage partitions. Every column remains represented. Dropping an operator coupling instead changes the operator and requires the error analysis developed next.

Back to source controls

Back to COS Method