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