ABRIK Wrap-Up - #130
Open
mmelnich wants to merge 30 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.
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