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{Gaussian quadratic forms}\label{ch:gaussian}
2 \chapterlead{A quadratic statistic may be presented as a matrix, a factor, or a list of
3 spectral coefficients. The representation determines how the calculation is
4 organized. Under the Gaussian model, the spectrum determines the law.
5 This chapter fixes that distinction before any numerical approximation is
6 introduced.}
7 \label{sec:gaussian-core}
8 
9 The primary mathematical object is a Gaussian quadratic form.
10 
11 \section{Central positive forms}
12 
13 Let
14 \[
15 X\sim N(0,\Sigma),
16 \qquad
17 M\succeq0,
18 \]
19 and write \(S=\Sigma^{1/2}\).  With \(X=SZ\), \(Z\sim N(0,I)\),
20 \[
21 Q=X^\top MX
22 =
23 Z^\top(SMS)Z.
24 \]
25 If the nonzero eigenvalues of the positive semidefinite operator
26 \[
27 H=SMS
28 \]
29 are \(\lambda_1,\ldots,\lambda_r\), then
30 \[
31 \boxed{
32 Q\overset d=
33 \sum_{j=1}^r\lambda_jZ_j^2.
34 }
35 \]
36 
37 Here and throughout the central reductions, the variables $Z_j$ are independent
38 standard normal variables. The numerical tail problem therefore depends only
39 on a positive spectrum, regardless of how that spectrum was generated.
40 If there are no nonzero coefficients, then $Q=0$ almost surely and
41 $P(Q\ge q)=1$ for $q\le0$, while $P(Q\ge q)=0$ for $q>0$.
42 This deterministic case is handled directly, without a saddlepoint or an
43 inversion calculation.
44 
45 \begin{proposition}[Representation equivalence]\label{claim:final-2-1}
46 The same central quadratic-form law may be supplied by any of the following:
47 \begin{enumerate}
48 \item a nonnegative spectrum \(\{\lambda_j\}\).
49 \item a positive semidefinite matrix \(H\).
50 \item a factor \(A\) with \(H=AA^\top\) or \(H=A^\top A\).
51 \item a replayable streamed factor
52 \[
53 A=[A_1\ \cdots\ A_B].
54 \]
55 \end{enumerate}
56 The nonzero spectra of \(AA^\top\) and \(A^\top A\) are identical, and
57 \[
58 AA^\top=\sum_bA_bA_b^\top.
59 \]
60 \end{proposition}
61 
62 \begin{proof}
63 If $AA^\top u=\lambda u$ with $\lambda>0$, then $A^\top u\ne0$ and
64 $A^\top A(A^\top u)=\lambda A^\top u$. Conversely, an eigenvector
65 $v$ of $A^\top A$ with positive eigenvalue maps to the eigenvector $Av$
66 of $AA^\top$. These maps restrict to inverse isomorphisms of the
67 corresponding eigenspaces after scaling, so multiplicities agree.
68 Horizontal concatenation gives $AA^\top=\sum_b A_bA_b^\top$ by direct
69 matrix multiplication. Orthogonal diagonalization then identifies the
70 associated Gaussian quadratic-form laws.
71 \end{proof}
72 
73 This proposition is the domain-independent source of the memory-bounded
74 streaming method developed below.
75 
76 \section{Central indefinite forms}
77 
78 The positive case is not the natural boundary of the central theory.  Let
79 \[
80 X\sim N(0,\Sigma),
81 \qquad
82 M=M^\top,
83 \]
84 where $M$ may be indefinite, and again write $S=\Sigma^{1/2}$.  The symmetric
85 operator
86 \[
87 H=SMS
88 \]
89 has a real signed spectrum.
90 
91 \begin{theorem}[Exact central signed reduction]\label{claim:final-2-2}
92 If $\lambda_1,\ldots,\lambda_r$ are the nonzero eigenvalues of $H$, then
93 \[
94 \boxed{
95 X^\top M X
96 \overset d=
97 \sum_{j=1}^r\lambda_jZ_j^2.
98 }
99 \]
100 \end{theorem}
101 
102 \begin{proof}
103 Write $X=SZ$ with $Z\sim N(0,I)$ and diagonalize
104 \[
105 SMS=V\Lambda V^\top.
106 \]
107 Rotational invariance of the standard normal gives
108 \[
109 Z^\top V\Lambda V^\top Z
110 \overset d=
111 \sum_j\lambda_jZ_j^2.
112 \]
113 No positivity assumption on $M$ is needed for this central reduction.
114 \end{proof}
115 
116 The sign pattern determines the support.  If every $\lambda_j\ge0$, the lower
117 endpoint is zero.  If every $\lambda_j\le0$, the upper endpoint is zero.  If
118 both signs occur, the distribution is unbounded in both directions.
119 
120 For the right-tail path, the cumulant generating function is
121 \[
122 \kappa(s)
123 =
124 -\frac12\sum_j\log(1-2s\lambda_j),
125 \]
126 with real domain
127 \[
128 \boxed{
129 \max_{\lambda_j<0}\frac{1}{2\lambda_j}
130 < s <
131 \min_{\lambda_j>0}\frac{1}{2\lambda_j}.
132 }
133 \]
134 Here a missing negative or positive spectral set gives a domain endpoint
135 \(-\infty\) or \(+\infty\), respectively.
136 When $q>E[Q]=\sum_j\lambda_j$, the saddlepoint is the unique positive solution
137 of
138 \[
139 \kappa'(\theta)=q
140 \]
141 inside this domain whenever a positive spectral component is present.  The
142 negative-only case with at least one nonzero coefficient has a positive
143 saddlepoint only for $E[Q]<q<0$. Indeed, $\kappa'(\theta)$ increases from
144 $E[Q]$ to zero as $\theta\to\infty$. For $q\ge0$ its survival
145 probability is zero, and no finite positive saddlepoint solves the query.
146 
147 The packet representation is sign aware: positive and negative eigenvalues are
148 clustered separately, and each same-sign packet stores a center $\mu$ and
149 centered moments.  The local expansion remains
150 \[
151 -\frac12\log(1-2s(\mu+\delta))
152 =
153 -\frac12\log(1-2s\mu)
154 +
155 \frac12\sum_{n\ge1}
156 \frac{\alpha(s)^n\delta^n}{n},
157 \qquad
158 \alpha(s)=\frac{2s}{1-2s\mu},
159 \]
160 with the same geometric remainder condition
161 \[
162 |\alpha(s)|\max|\delta|<1.
163 \]
164 
165 A useful distinction appears in the Fourier tail.  Under the saddlepoint tilt,
166 \[
167 \left|\phi_\theta(t)\right|
168 =
169 \prod_j
170 \left[
171 1+
172 \left(
173 \frac{2\lambda_j t}{1-2\theta\lambda_j}
174 \right)^2
175 \right]^{-1/4}.
176 \]
177 If $r$ coefficients are nonzero, then asymptotically
178 \[
179 |\phi_\theta(t)|=O(|t|^{-r/2}).
180 \]
181 Thus certification difficulty depends strongly on spectral rank, not only on
182 how small the requested probability is.
183 
184 The signed law has an exact spectral reduction, but its numerical
185 certification requires bounds appropriate to its support and Fourier decay.
186 The finite-endpoint, mixed-sign, and low-rank controls are developed in
187 Chapter~\ref{ch:evidence}.
188 
189 \section{Noncentral positive forms}
190 
191 Let
192 \[
193 X\sim N(\mu,\Sigma),
194 \qquad
195 M\succeq0.
196 \]
197 Again set \(S=\Sigma^{1/2}\) and diagonalize
198 \[
199 SMS=V\Lambda V^\top.
200 \]
201 Writing
202 \[
203 c=V^\top SM\mu,
204 \]
205 each positive-eigenvalue direction satisfies
206 \[
207 \lambda_iY_i^2+2c_iY_i
208 =
209 \lambda_i
210 \left(
211 Y_i+\frac{c_i}{\lambda_i}
212 \right)^2
213 -
214 \frac{c_i^2}{\lambda_i}.
215 \]
216 
217 \begin{theorem}[Exact noncentral reduction]\label{claim:final-2-3}
218 For the positive eigenvalues \(\lambda_i>0\),
219 \[
220 \boxed{
221 Q
222 \overset d=
223 q_0
224 +
225 \sum_i
226 \lambda_i
227 \chi^2_1(\delta_i),
228 }
229 \]
230 where
231 \[
232 \delta_i=
233 \left(\frac{c_i}{\lambda_i}\right)^2
234 \]
235 and
236 \[
237 q_0
238 =
239 \mu^\top M\mu
240 -
241 \sum_i\frac{c_i^2}{\lambda_i}
242 \ge0.
243 \]
244 The noncentral chi-square components are independent. To state the
245 zero-shift condition precisely, put $T=M^{1/2}S$, $b=M^{1/2}\mu$,
246 and let $P_T$ be the orthogonal projector onto $\operatorname{range}(T)$.
247 Then
248 \[
249 q_0=\|(I-P_T)b\|^2.
250 \]
251 Thus $q_0=0$ exactly when $M^{1/2}\mu\in\operatorname{range}(M^{1/2}S)$.
252 In particular, $\Sigma\succ0$ is sufficient. For singular covariance the
253 shift is the squared deterministic component perpendicular to the random
254 range in these weighted coordinates.
255 \end{theorem}
256 
257 \begin{proof}
258 Expand
259 \[
260 (\mu+SZ)^\top M(\mu+SZ)
261 =
262 Z^\top SMSZ
263 +
264 2\mu^\top MSZ
265 +
266 \mu^\top M\mu
267 \]
268 and rotate by \(V\).  Completing the square gives the displayed
269 noncentral components.  For a null vector of \(SMS\),
270 \[
271 M^{1/2}Sv=0,
272 \]
273 hence
274 \[
275 v^\top SM\mu
276 =
277 (M^{1/2}Sv)^\top M^{1/2}\mu
278 =
279 0,
280 \]
281 so there is no unresolved stochastic linear term in a zero-eigenvalue
282 direction. For the positive eigenvectors $v_i$ of $T^\top T=SMS$,
283 the vectors $Tv_i/\sqrt{\lambda_i}$ form an orthonormal basis of
284 $\operatorname{range}(T)$. Their inner products with $b$ are
285 $c_i/\sqrt{\lambda_i}$. Hence the subtracted sum equals $\|P_Tb\|^2$,
286 leaving $q_0=\|(I-P_T)b\|^2\ge0$. Orthogonal rotation preserves
287 independence of the normal coordinates, so the squared components are
288 independent as stated.
289 \end{proof}
290 
291 The numerical interface implements this reduction directly.  In an isotropic
292 control,
293 \[
294 \Sigma=I,\qquad M=aI,\qquad a>0,
295 \]
296 the law collapses to the exact scaled noncentral chi-square distribution
297 \[
298 Q/a\sim\chi^2_d(\|\mu\|^2).
299 \]
300 
301 The numerical controls for this reduction appear in
302 Chapter~\ref{ch:evidence}. They test a floating analytic-ledger engine and
303 are not assigned the directed-enclosure guarantee of the central positive
304 backend.
305 
306 \section{Decomposition axes}
307 
308 The solver separates four logically different representations:
309 \[
310 \boxed{
311 \text{storage blocks}
312 \circ
313 \text{structural phases}
314 \circ
315 \text{Krylov/range compression}
316 \circ
317 \text{spectral tail packets}.
318 }
319 \]
320 Storage blocks answer where data reside.  Structural phases answer how the
321 operator propagates.  Krylov/range compression answers which subspace has
322 actually been resolved.  Spectral packets answer how the final transform is
323 evaluated.  These decompositions should not be conflated.
324 
325 
326 
327 
328 \section{A factor, a spectrum, and an exact law}\label{sec:worked-law}
329 We will follow one small example through the remaining parts of the book. Its entries are exact binary fractions, so the same matrix can be used for symbolic algebra and for a binary64 input without introducing an input-rounding ambiguity. Set
330 \[
331 A=\frac12
332 \begin{pmatrix}
333  1& 1& 1& 1& 1& 1\\
334  1&-1& 1&-1& 1&-1\\
335  1& 1&-1&-1&-1&-1\\
336  1&-1&-1& 1&-1& 1
337 \end{pmatrix}.
338 \]
339 Multiplication gives
340 \[
341 K=AA^\top=\frac12
342 \begin{pmatrix}
343 3&0&-1&0\\0&3&0&-1\\-1&0&3&0\\0&-1&0&3
344 \end{pmatrix}.
345 \]
346 The orthogonal matrix
347 \[
348 V=\frac12
349 \begin{pmatrix}
350 1&1&1&1\\1&-1&1&-1\\1&1&-1&-1\\1&-1&-1&1
351 \end{pmatrix}
352 \]
353 satisfies $V^\top K V=\operatorname{diag}(1,1,2,2)$. For $Z\sim N(0,I_4)$,
354 \[
355 Q=Z^\top KZ\ \overset d=\ Y+2W,
356 \qquad Y,W\text{ independent with law }\chi^2_2.
357 \]
358 The spectrum is therefore $(1,1,2,2)$ and $\mathbb E Q=6$. The six-dimensional feature Gram $A^\top A$ has these same four nonzero eigenvalues and two additional zeros.
359 
360 Here the probability law also has a closed form. The densities of $Y$ and $2W$ are $\tfrac12e^{-x/2}$ and $\tfrac14e^{-x/4}$ on the nonnegative half-line. Their convolution is
361 \[
362 f_Q(x)=\int_0^x\frac12e^{-y/2}\frac14e^{-(x-y)/4}\,dy
363        =\frac12\bigl(e^{-x/4}-e^{-x/2}\bigr).
364 \]
365 Integrating from $q$ to infinity gives
366 \begin{equation}
367 P(Q\ge q)=2e^{-q/4}-e^{-q/2},\qquad q\ge0.
368 \label{eq:worked-exact-tail}
369 \end{equation}
370 At $q=40$ the probability is approximately $9.079779837\times10^{-5}$. This formula will be used as an independent control, not as a shortcut inside the COS calculation. The next part begins with the factor $A$ and asks which operations may change its storage or basis without changing this law.

Back to source controls

Back to COS Method