Skip to content

FunNystrom++ V2 - #138

Open
mmelnich wants to merge 31 commits into
mainfrom
funnystrompp-v2-baseline
Open

FunNystrom++ V2#138
mmelnich wants to merge 31 commits into
mainfrom
funnystrompp-v2-baseline

Conversation

@mmelnich

Copy link
Copy Markdown
Contributor

WIP

mmelnich added 2 commits May 12, 2026 10:47
…entize

Reference-aligned Persson-Kressner port with bit-exact MATLAB↔C++
cross-validation at LAPACK noise floor.

- Phase 5: vendored RandLAPACK::LanczosFA + BlockLanczosFA from
  funnystrompp; util::upsize raw-pointer overload, util::symmetrize.
- Phase 6: SASO sketch (RandBLAS::SparseSkOp) + SkOp-aware first matvec
  via right_spmm through ExplicitSymLinOp overload.
- Phase 7a: Cholesky-fast / eig(YᵀY) HMT §5.1 path in NystromEVD_v2.
  5-rep median speedup vs SVD-pinv fall-back: 2.71×-3.80× across
  k∈{500,1000,1500,2000} at n=2000 (3.80× at k=n=2000, exceeds the
  ~2× claim in funnystrompp PR #132).
- Phase 7a-perf: force_fallback knob + tight t_specrec_ms inner timer.
- Gap port (#1,#2,#3,#5): optional f_zero correction, Phase 2 skip at
  k==n, gesdd(G)→syevd(G) in fall-back, SkOp overload on
  FunNystromPP_v2::call.
- Componentize: NystromEVD_v2 + workspace moved to drivers/rl_nystrom_evd_v2.hh.

12 cross-validation fixtures pass at tol=1e-11 (10/12 at strict
tol=1e-12; 2 misses are alg_poly.t2 ~1.4e-12 long-standing cancellation
noise). v2 unit tests 4/4. Full RandLAPACK suite 309/309.
@mmelnich mmelnich mentioned this pull request May 15, 2026
mmelnich and others added 11 commits May 25, 2026 12:32
…pace, 11-slot timing

Bundle #2 of the v2-polish cycle:

  * Rename rl_nystrom_evd_v2.hh → rl_nystrom_evd.hh; NystromEVD_v2 free
    function → NystromEVD; NystromEVD_v2_workspace → NystromEVD_workspace.
  * Convert workspace from std::vector<T> members to raw T* + int64_t
    pairs grown via util::upsize and freed in a destructor, matching
    PR #132 / rl_revd2.hh patterns. Output params now use
    (T*&, int64_t&) references-to-pointer.
  * Add 11-slot timing vector (alloc / syrf / matvec / gram / potrf
    / trsm / svd / post_svd / error_est=0 / rest / total) with per-
    phase chrono instrumentation; gated on ws.times_enabled. Layout
    mirrors PR #132's NystromEVD for cross-driver tooling parity.
    Slot 8 (error_est) is unused in this fixed-k variant; kept at
    index 8 to preserve the V1 layout.
  * FunNystromPP_v2: U, lambda, Y_2, fAOmega, Y0 promoted from stack
    std::vectors to persistent raw-pointer class members; nystrom_ws
    promoted to a class member so allocations are reused across
    repeated call() invocations. Destructor frees all five buffers.

FunNystromPP_v2's public call signature is unchanged.

Verification:
  * make -j8 + make install -j8 + benchmark rebuild: clean
  * v2 unit tests: 4/4 PASS (FullRankCapture err_t1 = 2.426e-16)
  * 12-fixture cross-validation at tol=1e-12: 10/12 pass, same 2
    alg_poly.t2 cancellation-regime cells as pre-refactor
    (rel_Δ = 4.241e-13, 1.399e-12). est at noise floor across all
    12 fixtures (rel_Δ ∈ [0.0, 3.835e-15]).
**Silence all RandLAPACK compiler warnings under `-Wall -Wextra`**

Resolves all 22 warnings emitted when building the workshop docker image
(`ghcr.io/rileyjmurray/randlapack-python-workshop:latest`).

**Categories fixed** (-Wunused-but-set-variable, -Wunused-parameter,
-Wsign-compare, -Wvolatile, -Wdangling-else):
- Removed dead code: `prevnormNR` in `rl_determiter.hh`; unused
`[rows_B, cols_B]` structured bindings in `rl_dense_linop.hh` /
`rl_composite_linop.hh`.
- Removed genuinely vestigial parameters: `lda` from sparse
`ABRIK::call` (LDA is meaningless for sparse — `test_ABRIK_general` uses
`if constexpr` to dispatch dense vs. sparse), `k` from
`test_csr_col_split`, `norm_A` from `test_CQRRPT_orthogonalization`. All
callers updated.
- For parameters that are unused only in one specialization but required
by template-overload symmetry (`density`/`layout` in `make_operator`),
used the existing `UNUSED(x)` macro from `RandBLAS/exceptions.hh` with
explanatory comments.
- Misc: braced `ASSERT_FALSE` in `test_rpchol.cc` for `-Wdangling-else`;
rewrote `x += i` on `volatile int` as `x = x + i` to avoid C++20
deprecation in `test_memory_tracker.cc`.

Added issue #137 to deal with clang warnings/errors (will involve
updating RandBLAS submodule).
)

Two RandLAPACK sources triggered clang `-Wunused-but-set-variable` and
`-Wmisleading-indentation` under `-Wall -Wextra`. GCC 13.3 heuristics do
not flag either, which is why they passed through the prior
warnings-cleanup PR #136.

`RandLAPACK/testing/rl_gen.hh:174`: removed dead `int idx = offset;` and
`++idx;` inside `gen_exp_singvals`; the variable was written but never
read.

`RandLAPACK/drivers/rl_abrik.hh:281, 286, 415, 419`: de-indented the
bodies of four `#ifdef RandBLAS_HAS_OpenMP` blocks so they align with
the surrounding function-body indent. The previous extra indent visually
chained off the preceding single-statement `if(this->timing)`, which is
what triggers
`-Wmisleading-indentation`. No control-flow change.

Verified inside ghcr.io/rileyjmurray/randlapack-python-workshop with
clang 20.1.8: 0 warnings under `-Wall -Wextra` (down from 46), 708/708
ctest pass. GCC 13.3 build unchanged.
Two related bugs were present in both the CPU (rl_bqrrp.hh) and GPU
(rl_bqrrp_gpu.hh) implementations of BQRRP. They are independent and the
fix is in two parts:

1. Termination condition (CPU and GPU) The loop termination used
`curr_sz >= n`. For wide A (m < n), curr_sz tops out at min(m, n) = m <
n, so this condition never fires on full-rank wide input. The main loop
exhausts maxiter iterations and falls through `return 0` without setting
`this->rank` or running the termination handler. Result: `rank` is left
at its default (0). Fix: use `curr_sz >= std::min(m, n)`.

2. reconcile_and_cleanup wide-tail copy (GPU only) The GPU's
pointer-swap dance (Algorithm 7 of the BQRRP paper, used instead of the
CPU's in-place column permutation) needs an explicit reconciliation step
at termination to bring data home from whichever buffer holds the latest
writes. The existing reconciliation: * Walks block-aligned column ranges
[0, processed) where processed = (iter + 1) * b_sz_const (or smaller if
block_rank != b_sz_const). * Conditionally copies the trailing region
past `processed` only when `iter != maxiter - 1` (intended to handle
early termination). For wide A at full completion: processed = m < n, so
a valid R12 tail exists in [m, n) but the existing condition `iter !=
maxiter - 1` skips the copy (full completion has iter == maxiter - 1).
Additionally, the J pivot vector has no trailing copy at all on the
odd-iter path, leaving stale iter-(k-1) values in J[m:n] for wide A.
Since iter k's col_swap pulls values from those positions into earlier
slots, the final J ends up with duplicates. Fix: replace the `iter !=
maxiter - 1` gate with `processed < n`, which handles both early
termination and wide A at full completion. Add a symmetric J wide-tail
copy on the odd-iter path. The same fix applies in both the block_zero
early-exit and the normal termination paths.

The CPU implementation does not need the second fix because it uses
in-place column permutation -- the factorization data ends up directly
in caller's A across iterations, so loop fallthrough only loses the
`rank` field. This explains why the wide CPU performance results in the
BQRRP paper (Fig 3.12) were unaffected by the bug -- timing didn't
depend on the rank output.

Tests added:
  * test_bqrrp.cc::BQRRP_wide_aspect            (CPU, m=2000, n=5000)
  * test_bqrrp_gpu.cu::BQRRP_GPU_wide_aspect    (GPU, m=1000, n=2000)

The upstream test suite previously had no wide-A coverage (all tests
were tall m=5000, n=2800 or square 1000 x 1000), which is why this bug
went undetected.

Verified via an equivalent fix in a standalone CUDA port of BQRRP:
residual went from ~5e-1 (with both bugs) to ~1.5e-15 (machine epsilon)
on a m=256, n=512 Gaussian input. The reverse-engineered fix in that
port matches the patches here.

Co-authored-by: Maksim Melnichenko <mmelnichenko@icsi.berkeley.edu>
Route every PR to @mmelnich as code owner. Pairs with forthcoming branch
protection on main that requires Code Owner review.
Converted 150+ printf statements to std::cout across library, tests, and
benchmarks (not using streams where applicable yet).

Issue #75.
…ion (#140)

Four coordinated changes preparing RandLAPACK for the bindings repo:

install.sh dependency discovery via BLASPP_INSTALL_DIR /
LAPACKPP_INSTALL_DIR / RANDOM123_INSTALL_DIR env vars; skip clones and
reuse existing installs when set. Also fixes a CUDA flag conflict when
external blaspp depends on CUDAToolkit.

New RandLAPACK/rl_exceptions.hh defines RandLAPACK::Error plus a single
`randlapack_require(cond) << "msg " << x;` macro backed by a RAII
StreamThrower whose destructor throws at end of full expression.
Registered in the umbrella header and CMakeLists.txt. The macro uses
`for(; !(cond); )` rather than `if(cond);else` to short circuit without
triggering -Wdangling-else inside unbraced outer ifs.

211 call sites across 19 RandLAPACK source files migrate from
randblas_require to randlapack_require with hand readable messages that
include the actual parameter values, e.g. `randlapack_require(b_sz >= 1)
<< "BQRRP::call: b_sz=" << b_sz << " must be positive"`. Polarity flips
from "error if BAD" to "require GOOD"; (long long) casts drop since
operator<< handles int64_t directly.

Driver entry point validation in every driver's call() and constructor,
closing the segfault paths that bindings would expose: BQRRP (8 checks),
RSVD (5), CQRRT (7), CQRRPT (8), REVD2 (6 across 2 overloads), ABRIK (9
across 3 overloads), krill_full_rpchol (6 new + 1 retained mu_size).
Coordinates RandLAPACK with RandBLAS PRs #158, #169, and #172. Each
change stands alone; the commit history reflects that grouping.

### Submodule bump

`6abdfcf` → `8417f4b` (current `origin/main`), covering #158 (test
infrastructure refactor), #160, #164, #169 (CMake improvements), #170,
#171, and #172.

### Required follow-on changes from #169

**`Random123::Random123` rename.** RandBLAS #169 retired the bare
`Random123` INTERFACE target in favor of the namespaced
`Random123::Random123 IMPORTED GLOBAL`. The `target_link_libraries`
entry in `RandLAPACK/CMakeLists.txt` is updated to match.

**`RandBLAS_DIR` layout.** RandBLAS #169 moved the installed config from
`lib/cmake/RandBLASConfig.cmake` to
`lib/cmake/RandBLAS/RandBLASConfig.cmake`. The fallback in
`RandLAPACKConfig.cmake.in` now points at
`${CMAKE_CURRENT_LIST_DIR}/../RandBLAS` instead of
`${CMAKE_CURRENT_LIST_DIR}/..`.

### Test migration (forced by #158)

RandBLAS #158 promoted the test-only comparison and sparse-data
utilities from `RandBLAS/test/` headers into the library proper under
`RandBLAS/RandBLAS/testing/`. They install with the library; consumers
reach them via the `RandLAPACK::RandLAPACK` transitive include path.
Eleven files updated (10 RandLAPACK tests + 1 extras test):

| Before | After |
|---|---|
| `"../../RandBLAS/test/comparison.hh"` |
`<RandBLAS/testing/comparison.hh>` |
| `"../../RandBLAS/test/test_datastructures/test_spmats/common.hh"` |
`<RandBLAS/testing/sparse_data.hh>` |
| `test::comparison::` | `RandBLAS::testing::` |
| `test::test_datastructures::test_spmats::` | `RandBLAS::testing::` |

`extras/test/CMakeLists.txt` drops the now-obsolete probe that added the
RandBLAS source tree to the include path to resolve a test-only header
that wasn't installed. Post-#158 those headers are part of the library
install and reach consumers through the package's normal include path;
re-pointing the probe at the new source location would re-introduce the
same anti-pattern, so it's deleted.

### Consumer-facing CMake target improvements (ported from #169)

- Build-time-only flags wrapped in `$<BUILD_INTERFACE:...>` so `-Wall
-Wextra` (and `-fsanitize=address` when `SANITIZE_ADDRESS=ON`) apply
in-tree but disappear from the installed export.
- `target_compile_features(RandLAPACK INTERFACE cxx_std_20)` advertises
the language requirement on the exported target.
- `RandLAPACKConfig.cmake` calls `find_dependency(OpenMP COMPONENTS
CXX)` directly when RandLAPACK was built with OpenMP, rather than
relying on RandBLAS' Config to incidentally pull it into scope.
- `EXPORT_LINK_INTERFACE_LIBRARIES` (deprecated since CMake 2.8.12)
removed from `install(EXPORT)`.

### `RandLAPACK_BUILD_TESTS` option

New `option(RandLAPACK_BUILD_TESTS "Build RandLAPACK's test suite" ON)`
gates `add_subdirectory(test)` at the top level. The name is prefixed so
it can coexist with RandBLAS's own `BUILD_TESTS` when RandBLAS is
embedded as a subdirectory; the two options control their respective
test subdirs independently.

`install.sh` is updated to pass `-DBUILD_TESTS=OFF
-DRandLAPACK_BUILD_TESTS=ON` to the RandLAPACK cmake invocation: the
end-user installation path doesn't run RandBLAS's tests anyway (they're
covered by RandBLAS's own CI) and skipping them shaves measurable time
off the install.
@CLAassistant

CLAassistant commented May 30, 2026

Copy link
Copy Markdown

CLA assistant check
All committers have signed the CLA.

@rileyjmurray rileyjmurray left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I recommend you have Claude pull all changes from the benchmark/ files on this branch and move them to a new PR. Then you can reduce the scope of this PR to actually be about FunNystrom++.

mmelnich added 14 commits June 12, 2026 20:39
…e_bin helpers to testing utils, add LanczosFA/NystromEVD timing + reorth flag
…Alg 1 / Alg 2 / Lanczos-FA listings); code simplifications
…o-copy q=1 sketch, ping-pong subspace iters, lower-only block tridiag), pseudocode-tagged comments
…-sketch mode and the FunNystromPP SkOp overload; document fscalar/fAfun contract and symmetrize rationale
…-sketch mode and the FunNystromPP SkOp overload; RandBLAS-only RNG + raw buffers in test/benchmark; document fscalar/fAfun contract and symmetrize rationale
… at q=1, no densification), ws.Q rename kills the sketch alias, SkOp matvec required at compile time; drop wasted zero-inits
… ws.Q rename, SkOp matvec required; BlockLanczosFA: assemble T_k in place during the recurrence (drop A_blk/B_blk staging); drop wasted zero-inits
…ernalize normalized-Gaussian Omega2; FunNystromPP use_qfa switch + QFA/FA consistency test
…&C gesdd bug per BALLISTIC thread), loosen CQRRT orthogonality atol eps^0.7->eps^0.65
…owed relative-change stop, d_used); run_lanczos gains an optional per-step stop callback + steps_run; adaptive test
…inal syevd on adaptive stop and the (d*s)^2 copy in the fixed path (preserve_T flag); drop prose em-dashes in QFA/driver comments
…ccuracy, ~3 fewer steps; floor adaptive_min+delay = 4)
…tive-QFA depth probe + budget-closing allocation, rank cap n/2, s <= n-k, batched Phase-2 QFA); BlockLanczosFA: throw on block wider than n; BlockLanczosQFA: adaptive_delay default 5 -> 2
…stop plus fixed-depth mode; audit fixes (scalar-QFA include, non-adaptive pivot guard, auto-tier doc, k->n potrf message)
mmelnich added 4 commits July 28, 2026 13:41
…e-2 skip; name the block-QFA adaptive defaults; time reorthogonalization explicitly

Four fixes, three of them found while auditing for solver reuse.

1. NystromEVD was always called with vec_nnz (default 8) regardless of k. A
   SASO with more nonzeros per column than the column has entries is rejected
   by RandBLAS ("(vec_nnz <= dim_major) was required"), so ANY call with k < 8
   threw instead of degrading. This is not an exotic corner: the benchmark's
   smallest budgets give k = B/2 = 5 at B = 10, and the knob-free auto tier
   picks its own k and lands below 8 at small budgets (47 of 221 rungs in a
   full-campaign rehearsal died this way). Clamp to min(vec_nnz, k); vec_nnz = k
   is a dense column, the correct degenerate limit of a SASO.

2. t_fafun_ms was assigned only inside the k < n branch, so on a k == n call
   (Phase 2 skipped) it kept the previous call's value and consumers computing
   assembly = t_phase2_ms - t_fafun_ms got a negative number. Latent while every
   call constructed a fresh driver; a wrong-answer bug under reuse.

3. BlockLanczosQFA::adaptive_delay / adaptive_min are now exposed as named
   default_adaptive_* constants so a caller reusing an oracle can restore them
   explicitly. Relying on construction to supply the defaults leaks the previous
   call's convergence window forward, changing d_used, the estimate and the
   matvec count.

4. LanczosFA / BlockLanczosFA record _t_reorth_us and widen `times` to six slots
   (matvec, run_lanczos, apply_f, rest, total, reorth). Full reorthogonalization
   is O(d^2 n s) against the recurrence's O(d n s) and is the entire difference
   between this oracle and the basis-free LanczosQFA, so it is now a reported
   number rather than something inferred by differencing two tiers. The QFA
   classes emit 0 in that slot by design (no reorth).

Tests: SmallRankBelowVecNnzDoesNotThrow, FafunTimerResetAtKEqualsN and
ReuseAcrossCallsIsBitIdentical (a grow/shrink/grow sequence of (k, s) against
one reused driver, compared bit-for-bit against fresh instances) are new; the
second fails against the pre-fix driver. 18/18 TestFunNystromPPv2 pass.
# Conflicts:
#	RandLAPACK/CMakeLists.txt
#	RandLAPACK/drivers/rl_abrik.hh
#	install.sh
#	test/comps/test_determiter.cc
#	test/comps/test_qb.cc
#	test/comps/test_rpchol.cc
#	test/drivers/test_cqrrt.cc
#	test/drivers/test_krill.cc
…oes not exist on MSVC)

BinaryIoRoundTrip hardcoded /tmp and POSIX mkstemps/close, the only
compile failure in the whole branch under MSVC. save_dense_bin creates
the file itself, so no create-then-close dance is needed. Verified:
Linux suite 331/331; Windows (MSVC 19.44 + oneMKL ILP64) 767/767 via
install.ps1.
…eans budget/n-bounded)

The fixed default silently floored the oracle bias whenever the certified
depth exceeded 200: in the 2026-07-29 campaign (kappa >= 1e6, certified
depth 700-2200) the auto tier stagnated 10-38x above the 0.1 percent
target across the entire budget ladder. The probe is now bounded only by
n and half the budget; a positive auto_depth_cap remains an opt-in fixed
cap. Regression test pins the old failure (probe must exceed 200 on a
geometric kappa=1e6 spectrum at eps=1e-6; budget closure unchanged).
Suite 332/332.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants