ABRIK Wrap-Up - #130
Open
mmelnich wants to merge 48 commits into
Open
Conversation
mmelnich
force-pushed
the
winter-2025-abrik-clean
branch
from
May 4, 2026 18:45
124f1a7 to
5888ab1
Compare
ABRIK_speed_comparisons.cc and ABRIK_accuracy_analysis.cc both gain a num_runs
positional argument and an outer loop over runs that uses RNGState<RNG>(run) so
each run draws from a distinct seed. CSVs gain a leading 'run' column tagging
each data row with its run index. GESDD (deterministic) runs once and is
reported under run=0 in both benchmarks.
Speed_comparisons CLI: <prec> <outdir> <input> <target_rank> <run_gesdd>
<budget> <num_runs> <num_b_sz> <b_sz...> [sub_ratio] [use_cqrrt]
Accuracy_analysis CLI: <prec> <outdir> <input> <m> <n> <b_sz> <num_matmuls> <num_runs>
Previously the accuracy benchmark used RandLAPACK::gen::mat_gen with custom_input, which only handles whitespace-delimited text files. Feeding it the .bin matrices produced by gen_mat_alg971_paper crashed at the dimension query because the binary header was interpreted as ASCII rows. Switch to BenchIO::load_matrix (the same auto-detecting loader used by ABRIK_speed_comparisons and ABRIK_runtime_breakdown). It dispatches on file extension: .bin via read_bin_matrix, .txt via read_txt_matrix, .mtx via fast_matrix_market. Sparse input is now rejected explicitly with a clear error since accuracy_analysis runs a full GESDD. A is now owned by the LoadedMatrix and freed automatically; removed the matching new T[m*n] and delete[] A. CMakeLists.txt bumps the target's LINK_LIBS from Benchmark_libs to Benchmark_libs_external because ext_matrix_io.hh pulls in fast_matrix_market and Eigen headers, matching the other two ABRIK benchmark targets.
The summary `end_cols = iter * k / 2` truncates when iter*k is odd, which only
happens for odd block size k. For k=1, calling BK::call with max_krylov_iters=1
finishes the prelude + one odd half-step (Y_od[:,0:1] populated, R[0:1,0:1]
populated) but then returns end_cols=0 due to integer floor.
ABRIK::call_with_checkpoints then sees end_cols==0, emits the placeholder row
(total_matvecs=0, err=1), and breaks the entire checkpoint loop. Every
subsequent checkpoint (mv=2, 4, ..., 4096) is skipped, so the speed-comparison
CSVs contain only one sentinel row per run at b_sz=1.
The runtime-breakdown benchmark drives BK::call with fixed max_krylov_iters
>= 2 in every cell, so this path produced valid data — only call_with_checkpoints
at its first b_sz=1 checkpoint exposed the bug.
Fix: ceiling division `(iter * k + 1) / 2`. Bit-identical for even k (all
historical b_sz ∈ {4, 8, 16, 32, 64, 128} cases), and yields the correct
column count for odd k:
iter | k=1 floor | k=1 ceil | k=2 (no change) | k=4 (no change)
---- | --------- | -------- | --------------- | ---------------
1 | 0 (bug) | 1 | 1 | 2
2 | 1 | 1 | 2 | 4
3 | 1 (bug) | 2 | 3 | 6
4 | 2 | 2 | 4 | 8
No new tests added; the b_sz=1 path is now exercised by the existing
ABRIK_speed_comparisons benchmark at the bench level once Bergamo/SPR
campaigns rerun with this fix.
- BK: fix end_cols undercount for even block size with an odd final iteration (((iter+1)/2)*k instead of (iter*k+1)/2), which silently dropped k/2 singular triplets on rank-deficient / odd-max-iter paths. Bit-identical for even final iterations and for k=1 (the previously covered cases). - rs_linop: free Omega_1 on the final stabilization-failure return path (leak). - svd_residual: guard against k<1 and a zero smallest singular value. - QB/RSVD: document that the dense call() overwrites A (in-place deflation). - Drop redundant value-initialization on fully-overwritten U/V/Sigma outputs. - Benchmarks: remove the std::set_terminate backtrace debug handler; guard coo.reserve(0) (empty submatrix) and the budget_to_restarts divide-by-zero. - rl_downdatable_linop: use new[]/delete[] instead of calloc/free. - LinOp QB/RF/RS: note the dense and LinOp paths must be kept in sync. - Remove em-dashes and prose double-dashes from PR comments; fix ABRIK docstring typos.
mmelnich
force-pushed
the
winter-2025-abrik-clean
branch
from
July 7, 2026 20:36
e3e2d29 to
558f9e0
Compare
…&C gesdd bug per BALLISTIC thread), loosen CQRRT orthogonality atol eps^0.7->eps^0.65 (cherry picked from commit 3b8ba9c)
…t growth The adaptive criterion was assessed over every computed triplet rather than over the leading ones the caller asked for. On a decaying spectrum that cannot terminate early: each restart appends trailing triplets whose relative error is of order one, so the assessment is dominated by exactly the terms the restart just introduced, and it only passes once the Krylov subspace saturates. Measured on a spectrum decaying over six decades, the leading-10 residual fell from 3.5e-1 to 6.6e-15 while the all-triplets figure stayed near 1.7. In practice the driver always ran to end_cols = n and then reported failure, so the adaptive exit had never once fired before saturation. The number of assessed triplets is now derived at entry from the initial budget, k = ceil(max_krylov_iters / 2) * b, which is exactly what that budget produces. This needs no new parameter and makes an over-large request impossible to express: you cannot ask for more triplets than your own starting budget yields. k is fixed before the first restart and does not track the triplets subsequently computed, which is what lets deepening the subspace improve a fixed set. adaptive_increment is replaced by adaptive_growth (default 2.0), applied as p <- max(ceil(growth * p), p + 1). Doubling bounds the overshoot past convergence at 2x in iterations, about 4x in work. A larger ratio buys almost nothing: the per-check costs telescope to r^2/(r^2-1) of one check, 1.33 at r=2 against 1.01 at r=10, while the overshoot penalty grows as r^2. When adaptive is set and no budget was given, the initial budget defaults to 2, the smallest value satisfying Algorithm 1's p > 1 precondition. Adds ABRIKTermination and assessed_rank so callers can read why the loop stopped instead of inferring it from the residual, which cannot separate an exhausted retry budget from a saturated subspace. Adds ABRIK_adaptive_hard, the first benchmark to exercise adaptive mode, and a regression test on a rotated decaying spectrum asserting the iteration count is strictly below saturation. Every pre-existing adaptive test uses a flat Gaussian spectrum, where all triplets converge together and the two behaviors are indistinguishable, which is how this survived. Verified that the new test fails against the old criterion: it runs to iters=60 of 60 and terminates by rank deficiency rather than convergence. Suite: 314/314.
Records the derivation above the main loop: the Frobenius-content criterion (norm_R exceeding sqrt(1 - tol^2)||M||_F is equivalent to the relative residual bound, via the Pythagorean split of the two-sided projection) and the rank-deficiency criterion. Replaces the bare commented-out MATLAB expression on the convergence test with the argument behind it.
…lar value Brings the implementation in line with the residual the paper defines: each triplet's two-sided residual is divided by its own estimated singular value, rather than dividing the whole stacked residual by Sigma[k-1]. Previously the estimate was Sigma-weighted, so it reported the accuracy of the worst absolute residual scaled by the smallest retained singular value. An absolute threshold of that kind certifies only ~eps_mach*sigma_1/sigma_i relative accuracy for triplet i, and accepts triplets with sigma_i below the threshold vacuously. Dividing per triplet makes the quantity a Frobenius stack of relative residuals, so r <= eps bounds the relative backward error of every one of the k triplets. This changes reported accuracy everywhere svd_residual is used, including the err column for ABRIK, Spectra and RSVD in the speed benchmark, and ABRIK's own adaptive stopping test. Since each column is now divided by sigma_i >= sigma_k instead of the whole residual by sigma_k, the new value is bounded above by the old one, by a factor of up to the condition number of the retained block. All figures derived from the previous metric need regenerating.
…metrics in one pass Explicit assessed rank. The rank the error is assessed over is still derived from the initial budget by default, but can now be set explicitly. The derived value is a multiple of the block size, so a fixed count such as ten cannot be expressed at b = 4; an evaluation protocol that holds the assessed rank fixed while sweeping the block size needs the override, since block size is a performance knob while the assessed rank is a problem specification. Setting it too high for the initial budget throws, naming the minimum budget required. Honest under-delivery. A small residual over FEWER triplets than were requested is not convergence. The driver previously reported 'converged' in that case, because the assessment clamps to the triplets that exist. On an identity input, whose Krylov space is span(Omega) and never grows, asking for 100 triplets at b = 10 returned 10 and called it success. The two cases are now distinguished: the subspace may simply not have grown yet, which is benign, or it may be unable to grow at all, which is reported as ABRIKTermination::under_delivered. Three metrics in one pass. svd_residual_all returns the two-sided normalized residual (ours), the one-sided normalized variant, and the two-sided absolute variant, sharing the two operator applications rather than tripling the matvec cost. Needed to show the metrics side by side. On a decaying spectrum at eight Krylov iterations the one-sided variant reports 5.2e-15 while the two-sided reports 4.8e-4 on the same factorization, which is the failure mode a one-sided residual has by construction. Suite: 314/314.
…(C4); adaptive_hard: absolute-tolerance mode (B2)
Brings the branch up to date with everything merged this week: the quick-win batch (col_swap via lapmt with raw int64_t pivots, laset triangle helpers, typed constants, GPU single precision), install.sh 2.0 plus the installer CI, and the packager-only external-RandBLAS gate. Conflict resolutions: rl_qb.hh keeps this branch's removal of QB's internal matrix copy and takes main's typed-constant casts (two casts applied to branch-added reorthogonalization lines main never saw); rl_abrik.hh keeps this branch's restructured driver wholesale, since main's changes there were the mechanical cast sweep over the old structure; install.sh takes main's 2.0 rewrite unchanged; test_qb.cc keeps the gesvd reference SVD this branch already carried.
The two existing helpers aggregate into a single Frobenius norm over k triplets. That answers "how accurate is this set, taken together", which is what the adaptive loop needs, but it cannot answer "how many of these are real". A single norm mixes converged triplets with junk and reports one number. That second question is the one that matters for a block Krylov method that may commit basis columns carrying no operator content: such a column comes back as a triplet with sigma near zero. A fabricated direction cannot pass a two-sided NORMALIZED residual, so counting the triplets that do is an honest measure of delivered content, and an upper bound on what the algorithm may claim. svd_residual_per_triplet writes one residual per triplet; svd_triplets_certified returns how many clear a tolerance. Two operator applications, the same cost as either aggregate helper. A triplet with sigma <= 0 is reported as infinite rather than skipped or bailed on. The aggregate helpers return early in that case, which is right for a loop that only needs "not converged", but wrong here: a returned triplet with no singular value is precisely the failure this is meant to catch, so it must show up as one bad entry rather than poison the whole array.
Three separate concerns were being answered by one expression, abs(R_ii[(n + 1) * (k - 1)]) < sqrt(eps): when to terminate, whether the Krylov subspace can still grow, and which columns of the current block to keep. Only the third is what it was written for, and it was the only one working. In non-adaptive mode the iteration budget defaults to INT_MAX, so that exit never fires and only two remain. The other, the Frobenius criterion, was computing lantr over the UPPER triangle of an R that util::transposition stores LOWER, so norm_R was only ||diag(R)|| and the criterion essentially never fired. Rank deficiency was carrying termination single-handedly. That is why removing it broke nine tests with nothing to do with rank deficiency, and why a better criterion is less safe than a bad one until termination is fixed: a threshold relative to ||A|| is zero for a zero matrix. Termination first: - lantr now reads Uplo::Lower. norm_R climbs monotonically to ||A||_F again. - The norm_converged exit moved ABOVE ++iter. end_cols = ((iter + 1) / 2) * k reads iter as a count of COMPLETED iterations, and max_iters_reached breaks before the increment; norm_converged broke after it, so it reported a block that had never been built and gesdd consumed uninitialized basis columns. This was unreachable while the criterion above was dead, and surfaced immediately once it was restored: six passing tests went to residual 2.4 with unchanged iteration counts. - An explicit saturation guard, BKTermination::saturated, so the right basis can never exceed the ambient dimension. This is also a memory-safety bound: R_ii sits at row k*(iter_ev+1) with leading dimension n, so past saturation the band writes land in the NEXT allocated column. Silent corruption, no segfault, nothing for a sanitizer to see. It guards the odd side only, since even iterations append to the left basis, which lives in R^m. - Preconditions. k > min(m, n) was an out-of-bounds READ, not merely unsupported: the band buffers are n*k and (n+k)*k while the probes index (n + 1) * (k - 1). - The constructor initializes num_krylov_iters, norm_R_end and termination_reason, which were left uninitialized. Then the criterion, as block_numerical_rank, a free function so it can be unit tested. It reads the whole trailing block rather than one entry, measures against tau*||A||_F rather than an absolute constant, and returns a width so the healthy prefix of a partially deficient block survives. Following Balabanov, arXiv:2210.09953 Alg. 7, with the anchor carried across to the blocked setting: Alg. 7 measures against ||R||_2 of the factorization it reveals, and anchoring to a block's own scale does not survive blocking, because a wholly dead block has no healthy reference. The terminal block is truncated accordingly, and the two sides need DIFFERENT adjustments: on an odd terminal iteration the truncated block is a right-basis block, so end_cols shrinks while end_rows keeps the full left width and end_rows == end_cols stops holding; on an even one only end_rows shrinks. tau is user-facing, mirroring CQRRPT's eps, and forwarded through both ABRIK call paths. Also: R_11_trans is allocated once and reused every iteration, so a failed CQRRT left the previous block's healthy diagonal in place and the criterion read stale values and detected nothing. It is now cleared and the status honoured. ABRIK gained the zero-width guard its checkpoint sibling already had. Measured: a rank-5 input now returns 5 triplets rather than 10 of which 5 were noise, and the result no longer changes when the input is scaled by 1e8 (which previously gave 50 triplets and a different termination reason). Ill-conditioned input certifies 155 of 200 and a multiplicity wider than the block 200 of 200; both previously failed under every strategy tested. The test assertions change with the code because the old one could not see any of this: it was the two-sided UNNORMALIZED residual over a LEADING subset, and that variant accepts sigma <= eps vacuously while never inspecting the tail where junk lands. It is replaced by "every delivered triplet certifies" plus a normalized backstop.
BK was exercised only through the ABRIK driver, which hides end_rows/end_cols, the band buffers and BKTermination. Everything that made the rank-deficiency work hard to debug lives at that level. Four groups: Criterion unit tests. block_numerical_rank on synthetic diagonals, with no Krylov iteration, no BLAS and no RNG. Includes the two cases that decided the design: a uniformly dead block, where any rule anchored to the block's own scale sees a perfectly conditioned block and flags nothing, and scale invariance, which the previous absolute threshold failed and which alone would have caught the original bug. An interior dip does NOT truncate, and that is deliberate: the factorization is unpivoted, so a small entry in the middle genuinely does not imply the columns after it are junk, and testing the trailing BLOCK rather than the diagonal is what makes the unpivoted case safe. Liveness. BK must terminate with max_krylov_iters left at its INT_MAX default, on zero, identity, denormal-scaled, rank-one and full-rank input. That default had no coverage at all: every pre-existing test overrides it. Each asserts the band stays inside its buffers. An explicit CTest TIMEOUT is set because gtest has no per-test timeout, and it is the only mechanism that turns a non-termination regression into a failure rather than a hung CI job. Structural invariants. band == X_ev' * A * Y_od, on both the odd and even paths. Two GEMMs on a small matrix, and it catches a transpose slip, a permutation not folded into the band, and truncated columns left unaccounted for, at once. It also settles a documentation ambiguity: the buffer AS STORED is what equals X'AY, which is the orientation ABRIK hands to gesdd, despite the header calling it "stored transposed". Determinism and resume equivalence. Same seed gives bitwise identical output, and call(p) equals call(p1) + resume(p) bitwise. Recorded now because resume reconstructs its state as pure k-arithmetic and will be silently wrong once variable-width blocks land. Also the repo's first EXPECT_THROW: randlapack_require is not NDEBUG-gated, so preconditions are testable in Release.
Brings in the six commits main gained since the 2026-08-05 refresh (5da0674): 1cbe9e0 Native Windows (MSVC) support (#154) 4cc72a8 TEMPORARY: quarantine the macOS Accelerate gesdd test (#157) e4e8565 Windows installation overhaul (#156) 26eaed6 CI: fix the broken oneMKL cache path and drop a dead ctest exclusion (#161) ab05872 Installer overhaul: pins, provenance, run-not-just-link BLAS check (#162) d4ee721 Include the GPU layer from RandLAPACK.hh, guarded on __CUDACC__ (#169) Two conflicts, both in test/, both from main's MSVC work colliding with branch test infrastructure. Neither is library code. test/CMakeLists.txt: kept BOTH sides. Main adds /bigobj and randlapack_stage_runtime_dlls(RandLAPACK_tests); the branch adds PROPERTIES TIMEOUT 300 to gtest_discover_tests. The ordering is forced rather than stylistic: the DLL staging must precede gtest_discover_tests, because that command runs the test binary at build time to enumerate cases and on Windows that run needs the DLLs already staged. test/drivers/test_cqrrt.cc: both sides made the same code change, atol from eps^0.7 to eps^0.65, and differed only in the justification comment. Merged the comment to credit both backends. Two unrelated BLAS implementations independently put norm_0/sqrt(n) over the tighter bound on this case, vcpkg oneMKL sequential on Windows at 4.4e-11 and Apple Accelerate on macOS at ~1.7e-11. That it was found twice independently is what says the bound was wrong, not the backends. Three things worth recording. 1. 26eaed6 removed the dead ctest exclusion --exclude-regex "^TestABRIK\.ABRIK_catch_instability" that this branch was still carrying in six places across three workflow files, so the branch no longer needs to. One correction to that PR's description, which states there is no history of the test under that name: edab935^ contains four such tests (_prelim, _good, _bad, _worse) and the regex was unanchored on the right, so it matched all four until edab935 deleted them on 2026-02-02. The exclusion was live once, not always inert. 2. Main's 4cc72a8 quarantines TestQB.Polynomial_Decay_general1 on macOS, while this branch's 5a6b6d3 fixed that same test on Apple Silicon by switching its reference SVD from gesdd to gesvd. Both survive the merge, so the quarantine's own "if this passes, delete the suppression" warning will now fire on every macOS run. Left in place deliberately: 4cc72a8 is a temporary main-side commit pinned to two open upstream PRs, and reverting it here would fight main. Flagged for the PR thread instead. 3. install.sh shrinks by 484 lines because ab05872 turned it into a wrapper delegating to install/install.sh. Nothing was lost, and the branch's own cluster-install optimisation from 2ca557e survives there as -Dbuild_tests=OFF for both blaspp and lapackpp, using those projects' actual lowercase option name rather than the branch's BUILD_TESTING. The RandBLAS submodule advances from 8417f4b (1.1.0-32) to 952251c (1.1.0-42).
test/main.cc defines its own main() calling InitGoogleTest + RUN_ALL_TESTS, but no CMakeLists references it: it is absent from RandLAPACK_test_srcs, and RandLAPACK_tests takes its main from GTest::Main (test/CMakeLists.txt). So the file has never been compiled. It is not merely redundant: adding it to the source list would duplicate main() at link time against GTest::Main, and meanwhile it misleads anyone reading test/ about where main comes from. Verified dead before removal: git grep -n "main\.cc" -- '*CMakeLists.txt' '*.cmake' (no hits)
Two documentation sites had drifted behind the code, both in ways that read as the opposite of the truth. rl_bk.hh described R as an "Upper band matrix (stored transposed)". Measured, the buffer AS STORED equals X_ev' * A * Y_od to 7e-16 (TestBK.BK_band_equals_XtAY_*), and that is also the orientation ABRIK hands to lapack::gesdd with no transpose. The phrase described the per-block transposition done by util::transposition(..., copy_upper_triangle=1), which leaves each diagonal block lower triangular; it never meant the band as a whole. Both R and S now also state their dimensions and leading dimensions (n and n + k), and why S carries an extra k rows. rl_abrik.hh documented U, V and Sigma as exactly ((num_iters / 2) * k) wide. That was wrong twice over. The integer division is (iter + 1) / 2, not iter / 2 (rl_bk.hh, full_cols), and terminal-block truncation makes even the corrected expression an upper bound rather than a width. What the caller actually gets is singular_triplets_found, which rl_abrik.hh sets to end_cols and which is what U, V and Sigma are allocated against. Also drops two stale line-number citations from the test_bk.cc header comment and rephrases the ambiguity it settled in the past tense, since the header it was describing is fixed by this commit. Cited by symbol rather than line number so the comment does not rot again. No functional change.
test/drivers/test_rsvd.cc allocated A with new double[m * n]() and never freed it; the cleanup block covered A_saved, A_copy and A_check but omitted A. 160000 bytes per run. The sibling TestRSVD.LinOpSparse already cleans up correctly. The delete cannot be hoisted: A is aliased into A_linop and is destructively modified in place by the axpy preservation check, whose result is asserted just before the cleanup block. So the cleanup block is the first safe point. Worth recording why CI never caught this, because the answer is not "the asan job is broken". LeakSanitizer is active in the Debug/asan configuration and the leak is real, but it is masked in exactly the configuration CI runs: ./bin/RandLAPACK_tests --gtest_filter=TestRSVD.LinOpDense -> exit 0, no report LSAN_OPTIONS=use_registers=0 (same filter) -> exit 1, leak reported ./bin/RandLAPACK_tests --gtest_filter='TestRSVD.*' -> exit 1, leak reported At process exit a CPU register still holds the pointer, so LSAN roots it and classifies the block reachable rather than leaked. use_registers=0 unmasks it; use_stacks=0 does not, so it is a register artifact and not a stack one. And gtest_discover_tests gives every gtest case its own process, which is precisely the single-test configuration that masks it. Running more than one test in the process clobbers the register and the leak reports. Consequence: the asan job is green, would have stayed green without this fix, and cannot be relied on for this class of leak. Deliberately NOT addressed by setting LSAN_OPTIONS=use_registers=0 in CI, which is a fragile knob that would surface reachable-at-exit allocations across the whole suite. If the coverage is wanted, the safe form is a separate CI step running the binary once unfiltered as a leak sweep, and that should not be added before confirming the unfiltered run is otherwise leak-clean. Note macOS asan ships no LeakSanitizer, so the two asan jobs were never equivalent for leaks in the first place. Regression check must use a multi-test filter; a single-test filter cannot distinguish fixed from masked.
residual_error_comp assigns into the fixture-owned all_data.U_cpy and V_cpy with new T[...] without freeing any previous allocation, while both fixture destructors free only whatever pointer is currently stored. A second call on one TestData therefore leaks the first pair. No current test does this: every caller invokes it once per TestData instance, so this is a latent-leak shape rather than a live leak. But the signature invites the second call, and the Debug/asan job would report it the moment one appears. Safe as written because both ABRIKTestData and ABRIKTestDataSparse initialise U_cpy and V_cpy to nullptr, so the delete[] is a no-op on the first call.
The early `return i - 1` inside the search loop was wrong in two independent
ways.
It undercounted. s is non-increasing, so the first index at or below the
threshold IS the rank: s[0..i-1] are the i retained values. Returning i - 1
reported one fewer, and for an all-zero matrix (where s[0] is 0 and the very
first index trips the threshold) it returned -1.
And it leaked. The delete[] pair sat after the loop, so the early return skipped
it, leaking m*n + n elements on every rank-deficient input, which is precisely
the class of input this function exists to inspect.
Restructured to a single exit, shaped after cond_num_check directly above it:
that function is the copy-paste twin which got the frees right, which is
presumably how the two diverged.
On safety of the behaviour change. Both callers are in rl_gen.hh (the spiked and
adverserial cases of mat_gen), and both are guarded by `if
(info.check_true_rank)`. A repo-wide grep finds check_true_rank set true
nowhere: the only hits are its declaration (default false) and those two reads.
So the function is unreachable in the current tree and the corrected return
value is inert today. That is the argument for fixing it now, before something
starts depending on the wrong answer. Note also that mat_gen writes info.rank
from it rather than reading it, so any future consumer would have inherited the
off-by-one.
Three tests added in test/misc/test_util.cc, deterministic so the expected
answer is exact rather than threshold-dependent. Demonstrated against the
unfixed function before landing:
test_rank_check_full_rank passed before and after (fall-through path was
already correct)
test_rank_check_exact_rank returned 6 for a rank-7 input, now 7
test_rank_check_zero_matrix returned -1, now 0
Cherry-picked from 764a082 on spring-2026-wip (2026-07-28), then reparameterised to match the convention in this file. The original analysis stands and is restated here because that commit is not on this branch's history. The matrix built to break CholeskyQR was in fact the most benign input possible: gen_bad_cholqr_singvals returned all ones, so the condition number was 1 for every value of cond requested. Three faults compounded. 'int offset = k' made the decay loop 'for (i = offset; i < k; ++i)' empty. The rate log(1e8/cond)/(1-(n-offset)) was written for a block of n-offset entries while the loop wrote into a length-k vector, and it divides by zero when n == k+1. And the signature could not express the intent: the dispatcher passes info.rank as k, and rank defaults to n, so both arguments are equal by default and neither could supply the count of leading ones. There were no callers outside gen_bad_cholqr_mat, no test and no benchmark, which is why it survived. Blame is 9bae33d, which split the sigma helper out of the builder and left the offset behind. What changed relative to 764a082: that commit hardcoded offset = k / 2 and its own message conceded the split was a free choice, since the parameter that would have carried it was unusable. This version takes frac_spectrum_one instead, which is what the rest of the file does. gen_poly_singvals takes the fraction as a parameter and computes floor(k * frac_spectrum_one); gen_exp_singvals and gen_step_singvals hardcode floor(k * 0.1); and mat_gen_info already carries frac_spectrum_one defaulting to 0.1 and already threads it through for the polynomial case. So the plumbing existed and only bad_cholqr was not using it. Note this changes the default spectrum shape from half-and-half to 0.1/0.9, deliberately: it matches every sibling generator and is now a caller-supplied knob rather than a constant buried in the body. The unusable second dimension argument is dropped. Guards throw when the leading block would be empty, when the trailing block has fewer than two entries (the interpolation needs both endpoints), and when cond < 1e8, below which the trailing block would rise from 1e-8 toward 1/cond rather than decay, leaving a non-monotone spectrum whose condition number is 1e8 rather than the request. That threshold is also where the modeled failure begins, since unshifted CholeskyQR loses orthogonality past eps^(-1/2), about 1.5e8 in double. Endpoints are exact by construction: i = 0 gives 1e-8, and i = n_decay-1 gives 1e-8 * (cond*1e-8)^-1 = 1/cond, so s[0]/s[k-1] is exactly cond. Tests follow in the next commit, deliberately separate so that reverting this one leaves a failing test.
Deliberately separate from the fix in the previous commit, so that reverting that
commit alone leaves a failing test rather than silently restoring an all-ones
spectrum.
Four tests in a new TestGenSpectra fixture. The existing fixture in this file,
TestGeneratorsMutateState, is about RNG state advancement and is a different
concern, so this does not extend it.
bad_cholqr_singvals_realises_cond s[0] is exactly 1 and s[0]/s[k-1] is
the requested cond. Load-bearing: it
is the assertion that fails outright
without the fix.
bad_cholqr_singvals_block_sizes exactly floor(k*frac) leading ones,
and the cliff between the blocks is
1 -> 1e-8.
bad_cholqr_singvals_rejects_degenerate_shapes
cond < 1e8, an empty leading block,
and a trailing block below two
entries all throw.
bad_cholqr_singvals_is_monotone non-increasing across the whole
vector at cond 1e8, 1e10, 1e12.
Demonstrated against the unfixed generator before landing: three of the four
fail. The fourth, is_monotone, PASSES on the broken version, because a constant
all-ones vector is non-increasing. That is worth recording as the reason the
original fault survived so long, and the reason monotonicity alone is not an
adequate check for this generator.
Not included, and flagged rather than written blind: the end-to-end orthogonality
contrast that is the actual behaviour this generator exists to produce (764a082
measured CholQR at 9.86e-01 orthogonality error at cond 1e8, with CholQR2 and
sCholQR3 degrading gracefully). That belongs in test_orth.cc and needs its
tolerance chosen from a measured local run, not guessed.
…stride bug
Two changes that belong together because the second was found by testing the
first.
1. The even branch discarded CQRRT's return value while the odd branch checks it
and sets rank_deficient with final_block_width = 0. Now both match. Unlike the
odd branch there is no stale-buffer hazard here: S_ii is a fresh region of S
each iteration, zeroed by the initial calloc or by the fill after realloc, so a
failed factorisation leaves zeros rather than a previous block's diagonal, and
block_numerical_rank would return 0 and exit anyway.
This is DEFENSIVE, not a fix for observed behaviour, and the commit says so
because the measurement says so. A sweep of exact ranks 11 to 32 at block size
10 under cqrrt never produced a zero-width even-side block:
r not a multiple of 10 -> even terminal, final_block_width = r mod 10, in [1,9]
r a multiple of 10 -> stops one iteration EARLIER via norm_converged,
because the preceding odd iteration had already
captured all of A's spectral content
So block_numerical_rank always fires first with a nonzero healthy prefix, and a
fully dead even-side block does not arise for exact-rank input. The check still
earns its place: relying on the probe means relying on the buffer happening to
be zeroed, and once a narrowed block continues rather than terminating, a fully
dead even-side block becomes reachable.
The geometry concern that gated this change is resolved, and resolved in favour
of the simple mirror. final_block_width = 0 on an even terminal makes the band
square, where the ordinary even shape is end_rows = end_cols + k, and the worry
was that a square window would strand the terminal coupling block. It does not:
the block a square window excludes is the REJECTED one, which is the point. At
the first even iteration S_i is column 0 and S_ii is row k, so the window
S[0:k, 0:k] is exactly the coupling block X_0' A Y_0 written by the
reorthogonalisation gemm.
2. check_band_identity computed S's leading dimension as n + (end_rows -
end_cols). That difference is k only when the terminal block is full width; it
is the last accepted width once the rank criterion truncates. So the helper read
the buffer at the wrong stride at any truncated even terminal and reported a
band-identity failure that was not there. Measured on the new test: 1.158
before, 1.732e-15 after, against a 1e-10 tolerance. The library's band write was
correct throughout; only the check was wrong. Now takes k and uses n + k, which
is what its own comment already claimed.
This mattered more than a test-only bug usually would. Prune-and-narrow
continuation makes truncated terminals the normal case rather than the terminal
one, so this helper would have produced a spurious "the band is broken" signal
throughout that work, against the very invariant it exists to protect.
New test BK_even_terminal_band_identity_cqrrt covers the truncated even terminal
under cqrrt, which nothing did: it pins the exit route, the width, the non-square
geometry (25 by 20), and the band identity there.
Full suite 814/814.
ABRIK_catch_instability_bad was deleted in edab935 (2026-02-02) along with its _prelim, _good and _worse siblings. It was the original instability signal: a block size that is a large fraction of the ambient dimension, driven to the full dimension. Worth having again now that BK has direct tests, an explicit saturation guard and a termination-reason enum, so a failure here is diagnosable rather than mysterious. Scaled down from 4000x4000 with b_sz 1000 to 800x800 with b_sz 200. The shape that mattered is preserved (b_sz = n/4, target_rank = n, Gaussian input); what is dropped is a size that allocates two 128 MB buffers before ABRIK runs and drives the Krylov space to 4000 columns, against the TIMEOUT 300 that every discovered test now carries, and under Debug + asan on a two-core CI runner. Measured rather than assumed, which is why the size is what it is: Release 1.27 s Debug/asan 0.92 s Both far inside the 300 s budget. The asan figure being the lower of the two is not an anomaly: the heavy arithmetic goes to Release-built MKL, and asan instruments only RandLAPACK's own code. One API drift fixed on the way in: the historical body set ABRIK.num_threads_max and num_threads_min, both of which were removed from ABRIK after this test was deleted. No current test sets them.
Phase 3 step 1 (with steps 2 and 3 folded in; see below). Pure refactor: behaviour is unchanged and the proof is that all 56 printed diagnostic lines across TestBK and TestABRIK are BYTE-IDENTICAL to the pre-change baseline, with 815/815 green. WHY. Every extent, offset and band pointer in the loop was derived from the iteration count under the assumption that each block is exactly k wide. That assumption is what blocks prune-and-narrow continuation, and it is also what made the original defect hard to see. WHAT REPLACES IT. Two unambiguous accumulators, x_cols and y_cols, holding the ACCEPTED column count of each basis. They advance only on acceptance, after the rank probe. The old curr_X_cols/curr_Y_cols were pre-advanced before the block they reserved was written, so they meant "accepted plus pending" and flipped meaning twice per cycle; deriving a reorthogonalisation extent from them would have reintroduced exactly the off-by-one-block ambiguity behind this bug. Plus w_last, the width of the most recently accepted block, and separate allocation high-water marks. The six band pointers are now DERIVED at the top of each branch rather than advanced at the end of the previous one, collapsing three copies of the same arithmetic (fresh init, resume reconstruction, end-of-branch advance) into one. Verified against the old expressions case by case, for example R_ii lands at &R[n*k + k] at iteration 3 and S_ii at &S[(n+k)*k + 2k] at iteration 4 under both schemes. resume() now RESTORES four scalars instead of reconstructing them from `iter` by fixed-width arithmetic, which was silently wrong the moment any earlier block was narrower than k. That deletes the hazard rather than fencing it off, and it is why refusing a narrowed resume is unnecessary: after continuation lands, a narrowed run can legitimately end at max_iters_reached, which is precisely the state ABRIK resumes from. lantr becomes trapezoidal, (x_cols, y_cols). x_cols >= y_cols always, and after a narrowing they differ, so a square window would silently drop the R_i rows of the last accepted X block. Identical to the old square call whenever nothing has narrowed. The saturation guard is restated as the memory bound it always was, x_cols > n, checked once per even iteration since x_cols changes nowhere else. end_rows and end_cols are read straight off the counters, and the parity-dependent truncation adjustment disappears because the counters are already per-side. final_block_width survives as a diagnostic, joined by narrowed_blocks so a test can assert the mechanism FIRED rather than that the answer came out right. WHY STEPS 1-3 LANDED TOGETHER. The plan sequenced "accept the healthy prefix" as a separate step. It cannot be: the old truncated formulas (end_cols = full_cols - (k - width), end_rows = full_cols + width) equal the new counters only if the prefix IS accepted. Keeping the old epilogue for one commit would have meant expressing "what the counters would be if the block were full width", which is the reconstruction being removed. Accepting the prefix and collapsing the epilogue are the same change, and the byte-identical gate covers both, so nothing is lost by combining them. ONE REAL BUG FOUND IN MY OWN FIRST ATTEMPT, recorded because it is the trap in this refactor. Regrouping buffer growth by index space (odd owns Y_od and R, even owns X_ev and S) is correct, but I first placed the growth AFTER the pointer derivation and the GEMM that writes into the new block. On the non-prealloc path that writes past the end of X_ev at iteration 2, which showed up as "realloc(): invalid next size" in BK_terminates_on_identity. The old code got away with growing X_ev in the ODD branch because that was a pre-advance for the next iteration, keeping the destination one block ahead of the write. Growth must now precede both the derivation and the write. Both branches order it that way.
Phase 3 steps 4 through 7. This is the change the phase exists for.
THE DEFECT. X_ev receives a block at the prologue AND on every even iteration
while Y_od receives one only on odd iterations, so the left basis runs
permanently one block ahead, reaches the numerical rank first, and the run
stopped with the right basis up to k columns short. A Krylov space missing even
one direction of the row space is not invariant, so the leading triplets stopped
converging at all, which is why the loss was total rather than proportional.
THE FIX. When the rank criterion returns 0 < r < w, keep the healthy prefix, set
w_last = r, and CONTINUE with a narrower block. Only a fully dead block (r == 0)
is terminal. No random replacement of pruned columns: measured unnecessary, and
it costs orthogonality.
Landed in four verifiable steps, each with its own measurement:
step 4 zero the rejected part of the diagonal band block (rejected COLUMNS on
the odd side, rejected ROWS on the even side, because R_ii is stored
lower and S_ii upper). Their true value is zero by the band's
block-bidiagonal structure and the error introduced is bounded by
tau*||A||, exactly the order the criterion just called negligible.
Verified a no-op: diagnostics byte-identical to baseline. It stops being
a no-op once a block continues, because those positions then become
interior to the reported band with nothing else overwriting them.
step 5 odd-side continuation. Two changes only, both predicted:
ABRIK_adaptive_rank_deficient takes one more half-step (iters 1 -> 2,
same 5 triplets), and BK_terminates_on_rank_one now exits via
norm_converged rather than rank_deficient, which is the more accurate
reason once the narrowed block continues and captures everything.
step 6 even-side continuation. This is what flips the sweep.
step 7 test tightening.
RESULT, rank sweep at block size 10, certified triplets:
rank before after
20 20 of 20 20 of 20
21-29 0 to 3 r of r
30 30 of 30 30 of 30
31-38 3 to 7 r of r
39 0 0 (see below)
40 40 of 40 40 of 40
So 20 of the 21 ranks now certify r of r, where before only the three multiples of
the block size did. The six named regimes: T2 (exact rank 25) goes from 20
claimed / 0 certified to 25 of 25, which is the case that motivated the whole
phase; T4 (15 numerically dead directions) improves from 180 claimed / 114
certified to 185 of 185; T3, T5, T6 and T1 are unchanged. The T2 trace is exactly
the prediction: iteration 4 narrows the left block to 5, iteration 5 completes the
right basis, iteration 6 probes zero width and stops, terminal even with
end_rows = end_cols = 25.
RANK 39, THE ONE HOLDOUT, AND WHY IT IS NOT THIS PHASE'S PROBLEM. It is threshold
sensitivity in tau, not stranding. With tau defaulting to n*eps (4.44e-14 at
n = 200) the tenth column of the iteration-6 block carries a reorthogonalisation
residual just above tau*||A||, so the block is accepted at full width, the left
basis reaches 40 columns for a rank-39 matrix, and the run exits through
norm_converged on an odd terminal with end_rows = 40 > end_cols = 39. That one
junk left column prevents every triplet from certifying. Measured across tau:
tau = n*eps (default) iters 7 reason norm_converged er 40 ec 39
tau = 1e-12 iters 8 reason rank_deficient er 39 ec 39
tau = 1e-10, 1e-8 identical to 1e-12
so any tau at or above 1e-12 makes rank 39 behave exactly like 37, 38 and 40.
Rank 39 measured 39 claimed / 0 certified BEFORE this phase as well, so
continuation neither caused nor fixed it. The default tau is deliberately left
alone: the ill-conditioned regime (kappa 1e10, 155 of 200) depends on tau being
small enough not to discard genuine trailing directions, which is the same
trade-off this whole thread turned on, and it is what the user-facing tau knob
added in Phase 2 exists for. Recorded as
TestBK.BK_rank_39_is_a_tau_sensitivity_not_a_shortfall, and the sweep asserts
rank 39 at its measured value so that fixing it fails loudly here.
TESTS. The sweep is renamed ABRIK_rank_sweep_certifies_full_rank and now asserts
r of r rather than characterising a shortfall, keeping the never-over-deliver
invariant at every rank. T2 tightened from EXPECT_LE to EXPECT_EQ.
BK_even_terminal_band_identity_cqrrt updated to the fixed geometry, and its band
identity still holds at 1.732e-15, so the band write is correct at the new
non-square shapes. narrowed_blocks is asserted where it matters, so a test cannot
pass by luck without continuation having fired.
VERIFICATION. Release 816/816, Debug+asan 816/816, and the whole-binary asan leak
sweep reports nothing.
Two additions, both closing gaps where nothing was checking a property the code depends on. RESUME ACROSS A NARROWING. This closes the one Phase 3 acceptance criterion that was still open. resume() was changed from reconstructing its state to restoring four saved scalars, precisely because the reconstruction (curr_X_cols = (1 + iter_ev) * k, curr_Y_cols = iter_od * k) is silently wrong once any block is not exactly k wide. But no test exercised that: the existing resume test uses a full-rank Gaussian that never narrows, so the hazard the restore was designed to remove was itself untested. Before continuation the case was unreachable, because narrowing terminated the loop. Continuation makes it reachable: a narrowed run now ends at max_iters_reached, which is exactly the state ABRIK resumes from (rl_abrik.hh, the adaptive path). Exact rank 25 at block size 10 narrows the left block to 5 at iteration 4, so a checkpoint at 5 sits after the narrowing. Measured, bitwise identical between the two legs: RESUME single(8): rows=50 cols=40 | 4 then resume(8): rows=50 cols=40 RESUME-NARROWED single(8): rows=25 cols=25 | 5 then resume(8): rows=25 cols=25 The shared helper asserts narrowed_blocks >= 1 at the checkpoint, so the narrowed variant cannot quietly degrade into a second copy of the full-width one if the configuration ever stops narrowing. MATVEC ACCOUNTING. BK applies the operator in exactly three places: the prologue, the odd branch (Trans), and the even branch (NoTrans). Nothing checked the count, so a stray apply or a reorthogonalisation pass routed through the operator would have cost matvecs unnoticed. A counting wrapper around DenseLinOp makes the assertions exact integers with no tolerance. Worth recording that the first version of these assertions was WRONG and the test caught my error rather than a code bug. I asserted total == num_krylov_iters; measurement gave 41 at iters=40 and 7 at iters=6. The prologue applies A once (X = A*Omega) BEFORE iter is incremented to 1, so it is real work that num_krylov_iters does not count. The correct relations, now asserted: total == iters + 1 n_trans == (iters + 1) / 2 one per odd iteration n_notrans== iters / 2 + 1 one per even iteration, plus the prologue Two cases: one that exhausts its budget, and one that terminates early through the rank criterion, so the count is checked against work actually done rather than against a budget the run happens to exhaust. fro_nrm is forwarded to the inner operator and deliberately not counted; it is not an operator application. Full suite 819/819.
check_band_identity checked only that the band EQUALS X'AY, which is blind to structure: a dense band that happened to equal X'AY would pass, and loss of the block-bidiagonal structure through reorthogonalisation drift is exactly the failure mode the identity cannot see. The structural predicate could not be derived by reading the code, so it was measured first with a block-norm map and only then asserted. At k = 10: R (odd terminal, ld n) S (even terminal, ld n+k) diagonal 2.2e-01 diagonal 3.8e-01 .. 2.2e-01 subdiagonal 1.3e-01 subdiagonal 1.7e-01 .. 6.7e-02 ABOVE diag EXACTLY 0 ABOVE diag ~5e-17, computed BELOW subdiag ~3.5e-17, computed BELOW subdiag EXACTLY 0 Both bands are block-bidiagonal in their significant entries, but they are STRUCTURAL MIRRORS in which side is untouched. The odd branch fills a ROW strip of R (the rows of the current X block, across all previously accepted Y columns) while the even branch fills a COLUMN strip of S. So each band has one side merely at roundoff and one side never written, and they are opposite sides. That asymmetry is why the first version of this assertion was wrong. It applied a single rule to both bands, demanding exact zero below the subdiagonal, and failed on R at entry (20,0) with 3.6e-17. R's below-band region is computed, not untouched. The check is now parity-aware and asserts exactness only where the measurement showed the buffer is genuinely never written. The bounds are stated in ENTRY indices (j >= i + k above, i >= j + 2k below) rather than block indices, so they remain valid once the rank criterion narrows a block: widths only ever shrink, so the true support only tightens and these bounds stay conservative. That matters because prune-and-narrow makes narrowed terminals ordinary, and this assertion now rides along on all three band tests including the narrowed cqrrt one. This also settles a question left open in the Phase 3 design notes, which read the S layout as implying a structurally zero diagonal block per column and could not reconcile that with the identity holding at 7e-16. The map shows the diagonal blocks carry the largest mass; that reading was an off-by-one in block indexing. Full suite 819/819.
…udget
The ABRIK paper withholds every per-matvec claim about Spectra because the
figures plot the budget it was HANDED rather than the work it did. This adds the
measurement that unblocks those claims.
The hard part already existed: BudgetedPartialSVDSolver::num_operations() has
been public since the wrapper was written. Nothing called it. run_svds now
reports it, converted from A'A applications to matvecs with A by the factor of
two that ext_budgeted_svd_solver.hh documents in its header, applied explicitly
rather than left for the reader.
The CSV gains actual_matvecs beside total_matvecs rather than replacing it, so
the existing figures keep rendering and the new claim becomes checkable by
comparison. ABRIK, RSVD and GESDD write the same value in both columns, honestly:
their budgets are expressed directly in matvecs, so requested and actual coincide
by construction. Only Spectra can differ, because its restart schedule is derived
from the budget rather than metered against it, which is precisely why the claim
was withheld.
THE PART THAT WOULD HAVE GONE WRONG SILENTLY. The MATLAB reader
(plotting/parse_abrik_csv.m) is POSITIONAL:
textscan(fid2, '%d %s %d %d %f %d', ...)
Reading a 7-field row with that 6-specifier format does not fail. It slides
actual_matvecs into err and truncates err into an integer elapsed_us, producing
plausible-looking, wrong figures with no error anywhere. Inserting the column
without touching the reader would have corrupted every plot that consumes these
files.
So the reader now detects the schema from the header (contains 'actual_matvecs')
and selects the matching format, alongside the two schemas it already handled.
Files written before this column get actual_matvecs = total_matvecs, so
downstream code can read the field unconditionally.
aggregate_abrik_runs.m carries it through, reduced with the TIME reducer rather
than the error reducer since it is a cost measurement.
Those two MATLAB files are edited on disk under ~/matlab/ABRIK_benchmark/, which
is not a git repository, so they are not part of this commit and are flagged
separately.
Verified by compilation only, per the standing rule that benchmarks are compiled
locally and validated on the cluster: the benchmark is a standalone CMake project
needing an installed RandLAPACK, so this was checked with a -fsyntax-only build of
the translation unit against the worktree headers, which exits 0. It has NOT been
run, so the claim that actual_matvecs lands at or below total_matvecs at every
checkpoint is unverified. A count ABOVE the requested budget would mean the budget
control never worked, and that is a finding to report rather than smooth over.
The adaptive-termination figure overshoots its tolerance, reaching 8.5e-15
against a 1e-5 or 1e-10 request, because the driver grows its retry budget by
doubling. Showing that the overshoot is an artifact of the growth factor rather
than of the criterion needs a finer factor to be requestable, and it was not:
ABRIK has always had adaptive_growth, and this benchmark never set it, so every
cell silently ran at the 2.0 default and there was no way to ask for anything
else from the command line.
Added as an optional argument after sub_ratio, defaulting to 2.0, so every
existing cell is unaffected.
TWO THINGS THAT WOULD HAVE MADE THE RESULT WRONG RATHER THAN ABSENT:
The retry cap hardcoded the doubling assumption:
adaptive_max_retries = ceil(log2(iters_max / iters_start))
At growth 1.25 that allows about 4 retries where the same 16-to-1024 span needs
about 15. The run would have been throttled and reported max_retries, which reads
as the finer factor failing to converge rather than as being cut off. Now
log_growth(iters_max / iters_start).
And the CSV header asserted "the driver doubles its own budget" as a fact about
the data. It now states the actual value, so the file cannot misdescribe itself.
Cluster cells: two new entries in gen_abrik_jobs.sh, adaptive_cc3_tolabs5_g125
and _tolabs10_g125, mirroring the existing pair. The trailing "1.0 1.25" is
sub_ratio then adaptive_growth; sub_ratio is positional and comes first, so it
has to be stated explicitly at its default to reach the growth argument. The
generated scripts inherit campus-bigmem, sapphirerapids and the 24h walltime from
the generator, which is deliberate: ai-tenn reaps at roughly 2h despite
advertising 3 days, and MKL is unstable on the AMD nodes.
NOT SUBMITTED, and not to be submitted without explicit instruction. A peer
session's Phase 6 campaign currently holds 6 of the 8 campus-bigmem slots for
about 12h, and these two cells plus an install need 3, so they would not fit
regardless. Sequencing agreed with that session.
Verified by compilation only (-fsyntax-only against the worktree headers, exit 0)
and by generating the sbatch scripts and reading back the emitted argument list.
The benchmark has not been run, so nothing here is claimed about what growth 1.25
actually produces.
…ually guarantees
The comment on block_numerical_rank cited Balabanov Thm 5.6 (arXiv:2210.09953) as
"the contract this buys": cond(retained) <= 10 n^1.5 r / tau. Read against the
paper, that citation claims more than it can. Nothing about the criterion changes;
only the claim made for it, and the test that now backs the weaker, true one.
THE THEOREM, verbatim from the paper:
Theorem 5.6. Let X have normalized columns. Consider Algorithm 7 using the
strong rank-revealing QR algorithm and tau >= 4 n^{3/2} r u. Let Theta be an
eps-embedding for X_(1:r) with eps <= 1/2. Under Assumptions 5.5 (possibly
excluding 26d and 26e), cond(X_(1:r)) <= 10 n^{3/2} r tau^{-1}.
with X an m-by-n matrix, r the revealed rank, u the unit roundoff.
Three problems with inheriting it, none fatal to the criterion, all fatal to the
claim:
1. THE HYPOTHESIS IS VIOLATED BY OUR OWN DEFAULT. tau >= 4 n^1.5 r u is a LOWER
bound on tau. At u = eps/2 and a k = 10 block that is 1.4e-13; reading n as the
ambient column count gives 1.3e-11 at n = 200. Our default tau_eff = n*eps is
4.4e-14, below both by 3x to 280x. The theorem does not apply at the tau we
ship.
2. IT IS VACUOUS WHERE IT DOES APPLY. At tau = 1.4e-13 the bound reads
cond <= 2.3e+16; at our default, 7.1e+16. Double cannot represent a condition
number above 1/eps = 4.5e+15, so the bound permits more ill-conditioning than
the arithmetic can express. Measured against a factor with a known geometric
diagonal, the realised conditioning is 1/tau exactly while the bound is
4573/tau: it never binds, at any tau, by a constant three and a half orders of
magnitude. An assertion on it would pass on every input including a broken one,
which is worse than no assertion because it reads as coverage.
3. IT IS A DIFFERENT ALGORITHM. Thm 5.6 is stated for Alg. 7 using STRONG
rank-revealing QR on a matrix with NORMALIZED COLUMNS. This criterion pivots
nothing and runs on a band block from geqrf or CQRRT. The shape of the rule is
borrowed from Alg. 7 step 3; the guarantee is not.
WHAT IS ASSERTED INSTEAD. The relative property, which is both true and the reason
tau is user-facing: raising tau retains fewer columns and bounds their
conditioning more tightly. Measured at k = 12:
tau=1e-20 retained=11 cond=1.0e+20
tau=1e-14 retained= 8 cond=1.0e+14
tau=1e-10 retained= 6 cond=1.0e+10
tau=1e-06 retained= 4 cond=1.0e+06
tau=1e-02 retained= 2 cond=1.0e+02
so cond(retained) tracks 1/tau. The test asserts monotonicity in both retained
count and retained conditioning, plus that the knob actually bites across the
range, since monotonicity alone is satisfied trivially by a constant.
The same conclusion arrives independently from measurement: rank 39 at block size
10 is the one rank in 20..40 that fails to certify under the default tau and is
correct from 1e-12 upward. Theory wants at least 1.4e-13, measurement at least
1e-12. Both say the default is small. It is still deliberately left alone, because
the ill-conditioned regime needs a small tau not to discard genuine trailing
directions, and that trade-off is exactly what a user-facing tau is for.
This closes the last open item from the Phase 3 plan, which had deferred the
contract test pending a reading of the paper. The answer to "should we assert the
bound" turned out to be no, with reasons.
Full suite 820/820.
ABRIK CSVs recorded what a run was ON -- precision, input matrix, size, target
rank -- but never WHICH BINARY produced them. A grep for GIT_COMMIT across
benchmark/bench_ABRIK/ and the benchmark CMakeLists returned nothing, and existing
result files confirm it.
The gap is invisible until two eras are compared. Every row is individually
honest, the file looks complete, and nothing says the numbers came from different
builds.
It stops being hypothetical for the very next campaign. The growth-1.25 cells
exist to be compared against the growth-2.0 cells, and that comparison spans two
builds BY CONSTRUCTION: 1.25 cannot run until a build carries the adaptive_growth
argument that the 2.0 cells predate. So the one comparison the campaign is for is
exactly the cross-build case, and there was nothing to name the builds involved.
Three-piece chain, matching the convention the CQRRT campaigns already use so a
single MATLAB checker covers both:
1. the install job records the SHA to
<root>/RandNLA-project/build/RANDLAPACK_COMMIT -- next to the binaries it
describes rather than in a script (all three ABRIK install variants);
2. each generated cell cats that file, exports RANDLAPACK_GIT_COMMIT, and echoes
PROVENANCE commit=... into the job log;
3. the four ABRIK writers echo it into their comment header as
"# RANDLAPACK_GIT_COMMIT=<sha>", via a new shared
benchmark/bench_ABRIK/abrik_bench_provenance.hh.
The binary reads an ENVIRONMENT VARIABLE rather than a compile-time define
deliberately: no CMake plumbing, no rebuild to re-stamp, and a broken chain
surfaces honestly as "unknown" instead of a stale SHA baked in at configure time.
Two mistakes of mine caught while wiring it, both worth recording because both
would have produced a confidently wrong SHA rather than a visible failure:
* I first placed the rev-parse after `git checkout` but BEFORE `git pull`, which
records the previous commit. The pull is what determines the built SHA.
* The generator's heredoc is unquoted, so an unescaped $(cat ...) evaluated at
GENERATION time on the laptop, where the file does not exist, baking the
literal string "unknown" into every cell. It has to be escaped so it evaluates
at run time on the cluster, the same way the existing hardware-provenance line
escapes $(hostname).
Verified end to end by simulating the chain locally: install writes the SHA, the
cell export reads it back, the binary emits
"# RANDLAPACK_GIT_COMMIT=b76edea3ed9447b8352ccdeec00128196e59a519", it degrades to
"unknown" with the variable unset, and the line matches the checker's expected
form of 7-to-40 hex in a leading comment. All four writers compile
(-fsyntax-only, exit 0).
NOTE FOR WHOEVER ADDS THE MATLAB CHECK. The CQRRT helper asserts a single build
per ERA. Copying that to ABRIK unchanged would instantly refuse the growth study,
whose whole point is comparing two builds. For this campaign the assertion belongs
at CELL level -- single build within a cell -- with the cross-build comparison
made explicit between cells rather than forbidden.
The slurm tree is not a git repository, so the install and generator edits are on
disk only and are not part of this commit.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
WIP