sketch_freq_directions.nx source
↩ module page · 477 lines · 14489 B
1// sketch_freq_directions.nx -- Liberty 2013 Frequent Directions sketch.
2//
3// Streaming low-rank matrix approximation. Given a stream of input
4// rows in R^d, maintain a sketch matrix B in R^(l x d) (l << n) such
5// that B^T B approximates A^T A to within an additive error bound
6// based on the tail singular values:
7//
8// || A^T A - B^T B ||_F <= || A - A_k ||_F^2 / (l - k)
9//
10// where A_k is the best rank-k approximation of A. Direct
11// matrix-valued analog of Misra-Gries frequency counting (Liberty's
12// own framing in the 2013 STOC paper).
13//
14// ALGORITHM:
15// For each input row x in R^d:
16// 1. If any row of B is "empty" (all zeros) -> copy x into it.
17// 2. Else -> shrink:
18// a. Compute SVD: B = U Σ V^T
19// b. delta := σ_l^2 (smallest singular value squared)
20// c. Shrunk singular values: σ_i'^2 := max(0, σ_i^2 - delta)
21// d. B := diag(sqrt(σ_i'^2)) V^T
22// -- by construction the last row has σ_l' = 0 (empty)
23// e. Copy x into the now-empty row.
24//
25// SVD WITHOUT f64 (the technical chokepoint -- see queued roadmap
26// 2026-05-10):
27// This reference impl uses POWER ITERATION with DEFLATION in Q14
28// fixed-point. Bounded to l <= 4, d <= 8 for this version --
29// convergence guarantee at small sizes is reliable; for larger
30// matrices the accumulated round-off in Q14 degrades the
31// shrink quality. v2 will swap in fixed-point Jacobi or
32// bidiagonal + QR for arbitrary size.
33//
34// MATURITY = REFERENCE_IMPL with EXPLICIT size limits. Substrate
35// callers requesting larger l or d get NX_FD_ERR_TOO_BIG at alloc.
36//
37// COMPOSITION:
38// - Complements sketch_kmeans1d (clustering family) for streaming
39// dimensionality reduction before clustering.
40// - Complements sketch_correlation for streaming low-rank
41// covariance estimation.
42// - Complements sketch_naive_bayes for feature reduction before
43// classification.
44
45import "syscalls.nx"
46import "sketch_types.nx"
47
48const NX_FD_Q14: i64 = 16384
49const NX_FD_L_MAX: i64 = 4
50const NX_FD_D_MAX: i64 = 8
51const NX_FD_POWER_ITERS: i64 = 32
52const NX_FD_POWER_TOL_Q14: i64 = 1 // 1/16384 ~ 6e-5 convergence
53
54const NX_FD_OK: i64 = 0
55const NX_FD_ERR_TOO_BIG: i64 = -1
56const NX_FD_ERR_DIM_MISMATCH: i64 = -2
57
58struct FreqDir {
59 l: i64,
60 d: i64,
61 b: *i64, // l * d entries, row-major, Q14 fixed-point
62 total_rows: i64,
63 scratch_v: *i64, // size d
64 scratch_w: *i64, // size d
65 scratch_bv: *i64, // size l
66 scratch_v_basis: *i64, // size l * d (deflation basis)
67 scratch_s2: *i64, // size l (singular values squared)
68}
69
70// === construction =================================================
71
72func nx_fd_alloc(l: i64, d: i64) -> *FreqDir {
73 if l < 2 { return 0 as *FreqDir }
74 if l > NX_FD_L_MAX { return 0 as *FreqDir }
75 if d < 1 { return 0 as *FreqDir }
76 if d > NX_FD_D_MAX { return 0 as *FreqDir }
77 let raw: *u8 = sys_mmap(72)
78 let fd: *FreqDir = raw as *FreqDir
79 fd.l = l
80 fd.d = d
81 fd.total_rows = 0
82 let b_raw: *u8 = sys_mmap(l * d * 8)
83 fd.b = b_raw as *i64
84 var i: i64 = 0
85 while i < l * d {
86 fd.b[i] = 0
87 i = i + 1
88 }
89 let sv_raw: *u8 = sys_mmap(d * 8)
90 fd.scratch_v = sv_raw as *i64
91 let sw_raw: *u8 = sys_mmap(d * 8)
92 fd.scratch_w = sw_raw as *i64
93 let sbv_raw: *u8 = sys_mmap(l * 8)
94 fd.scratch_bv = sbv_raw as *i64
95 let sbasis_raw: *u8 = sys_mmap(l * d * 8)
96 fd.scratch_v_basis = sbasis_raw as *i64
97 let s2_raw: *u8 = sys_mmap(l * 8)
98 fd.scratch_s2 = s2_raw as *i64
99 return fd
100}
101
102// === fixed-point utilities ========================================
103
104func nx_fd_isqrt(x: i64) -> i64 {
105 if x < 0 { return 0 }
106 if x == 0 { return 0 }
107 if x < 4 { return 1 }
108 var g: i64 = (x >> 1) + 1
109 var iter: i64 = 0
110 while iter < 64 {
111 let next_g: i64 = (g + x / g) / 2
112 if next_g >= g { iter = 64 }
113 if next_g < g {
114 g = next_g
115 iter = iter + 1
116 }
117 }
118 return g
119}
120
121func nx_fd_abs(x: i64) -> i64 {
122 if x < 0 { return -x }
123 return x
124}
125
126// === matrix indexing helpers ======================================
127
128func nx_fd_b_get(fd: *FreqDir, r: i64, c: i64) -> i64 {
129 return fd.b[r * fd.d + c]
130}
131
132func nx_fd_b_set(fd: *FreqDir, r: i64, c: i64, v: i64) -> i64 {
133 fd.b[r * fd.d + c] = v
134 return 0
135}
136
137func nx_fd_row_norm_sq(fd: *FreqDir, r: i64) -> i64 {
138 var s: i64 = 0
139 var c: i64 = 0
140 while c < fd.d {
141 let v: i64 = nx_fd_b_get(fd, r, c)
142 s = s + (v * v) / NX_FD_Q14
143 c = c + 1
144 }
145 return s
146}
147
148func nx_fd_empty_row_index(fd: *FreqDir) -> i64 {
149 var r: i64 = 0
150 while r < fd.l {
151 if nx_fd_row_norm_sq(fd, r) == 0 {
152 return r
153 }
154 r = r + 1
155 }
156 return -1
157}
158
159// === vector ops in Q14 ============================================
160//
161// Dot, axpy, norm — all keep Q14 scale via (a * b) / Q14.
162
163func nx_fd_dot(a: *i64, b: *i64, n: i64) -> i64 {
164 var s: i64 = 0
165 var i: i64 = 0
166 while i < n {
167 s = s + (a[i] * b[i]) / NX_FD_Q14
168 i = i + 1
169 }
170 return s
171}
172
173func nx_fd_vnorm_sq(v: *i64, n: i64) -> i64 {
174 return nx_fd_dot(v, v, n)
175}
176
177func nx_fd_vscale(v: *i64, n: i64, factor_q14: i64) -> i64 {
178 var i: i64 = 0
179 while i < n {
180 v[i] = (v[i] * factor_q14) / NX_FD_Q14
181 i = i + 1
182 }
183 return 0
184}
185
186func nx_fd_vsub_scaled(target: *i64, source: *i64, n: i64, factor_q14: i64) -> i64 {
187 var i: i64 = 0
188 while i < n {
189 target[i] = target[i] - (source[i] * factor_q14) / NX_FD_Q14
190 i = i + 1
191 }
192 return 0
193}
194
195func nx_fd_vnormalize(v: *i64, n: i64) -> i64 {
196 let norm_sq_q14: i64 = nx_fd_vnorm_sq(v, n)
197 if norm_sq_q14 == 0 { return -1 }
198 // sqrt of Q14-scaled quantity: isqrt(x * Q14) gives sqrt(x) * sqrt(Q14)
199 // We want sqrt of (norm_sq / Q14) * Q14^2 = norm_sq * Q14:
200 let norm_q14: i64 = nx_fd_isqrt(norm_sq_q14 * NX_FD_Q14)
201 if norm_q14 == 0 { return -1 }
202 let inv_q14: i64 = (NX_FD_Q14 * NX_FD_Q14) / norm_q14
203 nx_fd_vscale(v, n, inv_q14)
204 return 0
205}
206
207// === matrix-vector ops ============================================
208//
209// Bv: compute B @ v, store in out (size l).
210// Btw: compute B^T @ w, store in out (size d).
211
212func nx_fd_matvec_Bv(fd: *FreqDir, v: *i64, out: *i64) -> i64 {
213 var r: i64 = 0
214 while r < fd.l {
215 var s: i64 = 0
216 var c: i64 = 0
217 while c < fd.d {
218 s = s + (nx_fd_b_get(fd, r, c) * v[c]) / NX_FD_Q14
219 c = c + 1
220 }
221 out[r] = s
222 r = r + 1
223 }
224 return 0
225}
226
227func nx_fd_matvec_Btw(fd: *FreqDir, w: *i64, out: *i64) -> i64 {
228 var c: i64 = 0
229 while c < fd.d {
230 var s: i64 = 0
231 var r: i64 = 0
232 while r < fd.l {
233 s = s + (nx_fd_b_get(fd, r, c) * w[r]) / NX_FD_Q14
234 r = r + 1
235 }
236 out[c] = s
237 c = c + 1
238 }
239 return 0
240}
241
242// === power iteration SVD (with deflation) ========================
243//
244// Find top l singular vectors of B by repeated:
245// v <- B^T B v / ||B^T B v||
246// Then deflate: subsequent vectors orthogonalized against found ones.
247
248func nx_fd_basis_get(fd: *FreqDir, idx: i64, c: i64) -> i64 {
249 return fd.scratch_v_basis[idx * fd.d + c]
250}
251
252func nx_fd_basis_set(fd: *FreqDir, idx: i64, c: i64, v: i64) -> i64 {
253 fd.scratch_v_basis[idx * fd.d + c] = v
254 return 0
255}
256
257// Project v against previously-found basis vectors and renormalize.
258func nx_fd_deflate(fd: *FreqDir, v: *i64, found: i64) -> i64 {
259 var j: i64 = 0
260 while j < found {
261 // dot of v with basis[j]
262 var s: i64 = 0
263 var c: i64 = 0
264 while c < fd.d {
265 s = s + (v[c] * nx_fd_basis_get(fd, j, c)) / NX_FD_Q14
266 c = c + 1
267 }
268 // v <- v - s * basis[j]
269 c = 0
270 while c < fd.d {
271 let bv: i64 = nx_fd_basis_get(fd, j, c)
272 v[c] = v[c] - (s * bv) / NX_FD_Q14
273 c = c + 1
274 }
275 j = j + 1
276 }
277 return 0
278}
279
280// Compute one singular vector via power iteration, store in v_out
281// (size d). Returns sigma² in Q14 via fd.scratch_s2[found_idx].
282func nx_fd_power_one(fd: *FreqDir, found_idx: i64, init_idx: i64) -> i64 {
283 // Initialize v as unit vector e_{init_idx}.
284 var c: i64 = 0
285 while c < fd.d {
286 fd.scratch_v[c] = 0
287 c = c + 1
288 }
289 fd.scratch_v[init_idx] = NX_FD_Q14
290 // Deflate against previously-found vectors.
291 nx_fd_deflate(fd, fd.scratch_v, found_idx)
292 let nrm: i64 = nx_fd_vnormalize(fd.scratch_v, fd.d)
293 if nrm < 0 {
294 // Couldn't normalize -- direction was entirely in deflated subspace.
295 fd.scratch_s2[found_idx] = 0
296 return 0
297 }
298 // Iterate.
299 var iter: i64 = 0
300 while iter < NX_FD_POWER_ITERS {
301 // w := B v
302 nx_fd_matvec_Bv(fd, fd.scratch_v, fd.scratch_bv)
303 // v_new := B^T w
304 nx_fd_matvec_Btw(fd, fd.scratch_bv, fd.scratch_w)
305 // Deflate.
306 nx_fd_deflate(fd, fd.scratch_w, found_idx)
307 // Normalize.
308 let n2: i64 = nx_fd_vnorm_sq(fd.scratch_w, fd.d)
309 if n2 == 0 {
310 fd.scratch_s2[found_idx] = 0
311 return 0
312 }
313 let nrm_w: i64 = nx_fd_vnormalize(fd.scratch_w, fd.d)
314 if nrm_w < 0 {
315 fd.scratch_s2[found_idx] = 0
316 return 0
317 }
318 // Convergence check: |v - w| in Q14.
319 var diff: i64 = 0
320 var k: i64 = 0
321 while k < fd.d {
322 diff = diff + nx_fd_abs(fd.scratch_v[k] - fd.scratch_w[k])
323 k = k + 1
324 }
325 // Copy w into v.
326 k = 0
327 while k < fd.d {
328 fd.scratch_v[k] = fd.scratch_w[k]
329 k = k + 1
330 }
331 if diff <= NX_FD_POWER_TOL_Q14 { iter = NX_FD_POWER_ITERS }
332 iter = iter + 1
333 }
334 // Compute σ² = || B v ||²:
335 nx_fd_matvec_Bv(fd, fd.scratch_v, fd.scratch_bv)
336 let s2: i64 = nx_fd_vnorm_sq(fd.scratch_bv, fd.l)
337 fd.scratch_s2[found_idx] = s2
338 // Store v as basis[found_idx].
339 c = 0
340 while c < fd.d {
341 nx_fd_basis_set(fd, found_idx, c, fd.scratch_v[c])
342 c = c + 1
343 }
344 return 0
345}
346
347// === shrink =======================================================
348//
349// Compute SVD via power iteration, find σ²_min, shrink each
350// σ²_i by σ²_min, rebuild B from the shrunk singular values
351// and the right singular vectors V (stored in scratch_v_basis).
352//
353// New B[i] = sqrt(shrunk_σ²_i) * V_i
354
355func nx_fd_shrink(fd: *FreqDir) -> i64 {
356 var i: i64 = 0
357 while i < fd.l {
358 // Init seed varies per i to avoid degenerate starts.
359 let init_idx: i64 = i % fd.d
360 nx_fd_power_one(fd, i, init_idx)
361 i = i + 1
362 }
363 // Find σ²_min.
364 var s2_min: i64 = fd.scratch_s2[0]
365 i = 1
366 while i < fd.l {
367 if fd.scratch_s2[i] < s2_min { s2_min = fd.scratch_s2[i] }
368 i = i + 1
369 }
370 // Rebuild B[i] = sqrt(max(0, σ²_i - σ²_min)) * V_i in Q14.
371 i = 0
372 while i < fd.l {
373 var shrunk: i64 = fd.scratch_s2[i] - s2_min
374 if shrunk < 0 { shrunk = 0 }
375 // shrunk is in Q14 ( = (||B v||² in Q14)). sqrt(Q14) = Q7; we want Q14.
376 // σ' in Q14 = isqrt(shrunk * Q14).
377 let sigma_q14: i64 = nx_fd_isqrt(shrunk * NX_FD_Q14)
378 var c: i64 = 0
379 while c < fd.d {
380 let v_c: i64 = nx_fd_basis_get(fd, i, c)
381 let new_v: i64 = (sigma_q14 * v_c) / NX_FD_Q14
382 nx_fd_b_set(fd, i, c, new_v)
383 c = c + 1
384 }
385 i = i + 1
386 }
387 return 0
388}
389
390// === add row ======================================================
391
392func nx_fd_add_row(fd: *FreqDir, row: *i64, n: i64) -> i64 {
393 if n != fd.d { return NX_FD_ERR_DIM_MISMATCH }
394 var empty: i64 = nx_fd_empty_row_index(fd)
395 if empty < 0 {
396 // No empty row -- shrink.
397 nx_fd_shrink(fd)
398 empty = nx_fd_empty_row_index(fd)
399 // After shrink the last row should be empty (smallest singular
400 // value zeroed out). If still not empty (numerical noise),
401 // force-clear the row with smallest norm.
402 if empty < 0 {
403 var min_idx: i64 = 0
404 var min_norm: i64 = nx_fd_row_norm_sq(fd, 0)
405 var i: i64 = 1
406 while i < fd.l {
407 let nrm: i64 = nx_fd_row_norm_sq(fd, i)
408 if nrm < min_norm {
409 min_norm = nrm
410 min_idx = i
411 }
412 i = i + 1
413 }
414 var c: i64 = 0
415 while c < fd.d {
416 nx_fd_b_set(fd, min_idx, c, 0)
417 c = c + 1
418 }
419 empty = min_idx
420 }
421 }
422 var c: i64 = 0
423 while c < fd.d {
424 nx_fd_b_set(fd, empty, c, row[c])
425 c = c + 1
426 }
427 fd.total_rows = fd.total_rows + 1
428 return NX_FD_OK
429}
430
431// === queries ======================================================
432//
433// Approximate top singular value (after most-recent insertions).
434
435func nx_fd_top_sigma_sq(fd: *FreqDir) -> i64 {
436 // Just compute the row norm of the largest row -- adequate
437 // approximation for non-degenerate B; honest about being a
438 // lower bound on σ_1².
439 var best: i64 = 0
440 var r: i64 = 0
441 while r < fd.l {
442 let nrm: i64 = nx_fd_row_norm_sq(fd, r)
443 if nrm > best { best = nrm }
444 r = r + 1
445 }
446 return best
447}
448
449func nx_fd_total_rows(fd: *FreqDir) -> i64 {
450 return fd.total_rows
451}
452
453func nx_fd_memory_bytes(fd: *FreqDir) -> i64 {
454 return 72 + fd.l * fd.d * 8 + fd.d * 8 * 2 + fd.l * 8 + fd.l * fd.d * 8 + fd.l * 8
455}
456
457// === typed envelope ===============================================
458//
459// Frobenius approximation error envelope. For l rows, the FD bound
460// is || A^T A - B^T B ||_F <= 2 || A ||_F^2 / l (Liberty 2013, Thm 1.1
461// with k=0 worst case). Express as REL_STDDEV envelope on top
462// singular value estimate.
463
464func nx_fd_rel_err_ppb(l: i64) -> i64 {
465 if l >= 4 { return 500000000 } // 50% relative on smallest dimensions
466 if l >= 3 { return 666000000 } // 66%
467 return 1000000000 // 100% (l=2 worst case)
468}
469
470func nx_fd_query_top(fd: *FreqDir) -> *ApproxI64 {
471 let s2: i64 = nx_fd_top_sigma_sq(fd)
472 return nx_approx_new(s2, NX_ENV_REL_STDDEV,
473 nx_fd_rel_err_ppb(fd.l),
474 682700000,
475 NX_MATURITY_REFERENCE_IMPL,
476 NX_ADV_HONEST)
477}