CQRRT diagnostics#129
Open
mmelnich wants to merge 47 commits into
Open
Conversation
Contributor
Author
|
Closing this & adding a reference to the paper. |
Contributor
Author
|
Temporarily re-opened to add functionality. |
When solving for R_sk^{-1} explicitly to precondition a linop, solve
R_sk * X = I (Side::Left) rather than X * R_sk = I (Side::Right).
The two are mathematically equivalent but the latter exhibits poor
backward error in the subsequent product A * R_sk^{-1} when R_sk is
ill-conditioned. On photogrammetry2 (κ(R_sk) ≈ 2.2e8) this changes
orth_err in CQRRT_linop from 1.6e-4 to 1.5e-9 -- matching the
GEQP3/BQRRP stabilized variants and obviating their need as the
default path.
Same fix applied to the initial M = R_1^{-1} formation in
sCholQR3_linops.
Slim CQRRT_linop_applications to run only the patched CQRRT_linop;
swap GETRI for BQRRP in CQRRT_diagnostic to show the stabilized
counterpart alongside trsm/trtri.
…b bugs
- Replace NMR/Kronecker benchmark with CQRRT_linop_irlsq.cc.
Loads a tall sparse .mtx, wraps as SparseLinOp, generates synthetic
x_true + b = J*x_true + noise. For each Q-less QR variant: draws a
fresh sparse sketch S2 (independent of CQRRT's S1), forms x_0 =
R^{-1} R^{-T} (S2 A)^T (S2 b) (paper Algorithm 1, line 3, Q-less
form), then runs 2-step IR with inner CG.
- Strip Tikhonov from IterRefineLSQ. Drop KroneckerOperator,
RegularizedLinOp, and their tests/includes/CMake entries. Drop the
GEQP3-stabilized variant from the benchmark and rename
CQRRT_linop_stb_bqrrp -> CQRRT_linop_bqrrp.
- Fix IR-LSQ populate_times double-counting bug. Inner CG's
TRSM/fwd/adj contributions were tracked into both t_inner_total
(wallclock of the inner_cg call) and the shared
t_{trsm,fwd,adj}_total counters, so 'other = outer_total - inner -
trsm - fwd - adj' went negative. Split into outer-only and
inner-only counters and report inner_ex = t_inner_total -
t_inner_{trsm,fwd,adj}, total_trsm = t_outer_trsm + t_inner_trsm
(similarly fwd, adj). 'other' is now a non-negative residue
(axpy/copy/nrm2 bookkeeping).
- Fix three analytical_kb formulas in rl_memory_tracker.hh.
cqrrt_linops_bqrrp_analytical_kb captured only the BQRRP-precond
moment (d*n + 4n^2) and missed the later Gram-loop moment (d*n + n
+ n^2 + m*b_eff). For tall inputs the latter dominates. Now
returns the max of the two; signature extended to (m, n, d_factor,
block_size). scholqr3_linops_analytical_kb forgot G3_factor at the
iter-3 peak (5n^2 -> 6n^2). scholqr3_linops_basic_analytical_kb
forgot all three G_i_factor members (3n^2 -> 6n^2).
Per collaborator's correction (Oleg, 2026-05-25), the FEM2 operator is the
Petrov-Galerkin form J = B^{-1} * A * U_h with B = chol(M), A = K, U_h = V
(mass-matrix Cholesky factor, stiffness, prolongation respectively).
CQRRT_linop_applications now takes three .mtx files in FEM mode and builds
a doubly-nested CompositeOperator:
J = CompositeOperator(L_inv_op,
CompositeOperator(K_op, V_op))
with L_inv_op = CholSolverLinOp<T>(M_file, half_solve=true). Because L
factors M (not K), the composite does NOT algebraically collapse; the LS
normal equations land on J' J = V' K M^{-1} K V (Petrov-Galerkin coarse-
grid mass-weighted stiffness-squared).
CLI changes:
FEM mode (NEW): <K_file> <M_file> <V_file> <d_factor> [...] (3 files)
Sparse mode: sparse <A_file> <d_factor> [...] (unchanged)
The previous 2-stage composite (L^{-1} * V with L = chol(K)) is gone --
that formulation was based on a misread of the generator's outputs. The
sparse mode is unchanged.
No driver, linop, test, or CMake changes.
- Single binary now selects between SVD post-processing, IR-LSQ refinement,
or both via a positional <mode> arg.
- 5-method dispatch (CQRRT_linop, CholQR, sCholQR3, sCholQR3_basic,
CQRRT_linop_bqrrp) restored to the FEM/sparse composite paths.
- Add power-iteration estimate of ||A||_2 and replace ls_residual_norm
with the Higham normwise backward error ||Ax-b||/(||A||*||x||+||b||),
drivable to machine epsilon for a backward-stable LS solver.
- Add memlite blocked-compute orth_err (O(n^2 + m*b)), runs for every
selected method in every mode; new orth_error column in irlsq CSV.
- FEM + irlsq: b = L^{-1} * Gaussian random vector (no x_true).
- CQRRT_linop_irlsq.cc removed (CMake target dropped).
…s spec)
Introduces comps/rl_cholqr.hh with three free-function templates:
- blocked_preconditioned_gram(A, R_pre, G, ...) : Layer 0; computes
G = R_pre^T A^T A R_pre via blocked linop calls. nullptr R_pre handles
the M=I case via a small per-block identity scratch.
- cholqr_primitive(A, R, shift_factor, ...) : Algorithm 1 (with optional
shift for sCholQR3 iter 1).
- pcholqr_primitive(A, P, R, method, ...) : Algorithm 2; the unifying
building block. Dispatches the P^{-1} step on PCholQRPrecondMethod
(TRSM_IDENTITY / TRTRI / GEQP3 / BQRRP). Adds the new TRTRI method.
Gram step in pcholqr_primitive dispatches on method:
- TRSM_IDENTITY/TRTRI: per-block writes A^T A R_pre into G, then a
single TRSM(P^T, G) applies the left factor. O(n^3/2) vs O(n^3),
stable since P is preserved (the optimization CQRRT_linops used to
have inline; now shared with sCholQR3 iters 2-3).
- GEQP3/BQRRP: per-block GEMM with explicit R_pre^T, preserving the
QRCP stability advantage.
Driver refactor:
- CholQR_linops, sCholQR3_linops (both variants), CQRRT_linops are now
thin wrappers over the primitives. ~33% LOC reduction across the
family (1946 -> 1310 lines).
- rl_cqrrt.hh now holds both dense CQRRT and CQRRT_linops in one file;
rl_cqrrt_linops.hh deleted (CMake had no separate target; 3 benchmark
.cc files updated to drop the include).
- Backwards-compat alias `using CQRRTLinopPrecond = PCholQRPrecondMethod;`
keeps existing benchmark code unchanged.
Memory tracker formulas updated to reflect the now-freed buffers:
scholqr3_linops_analytical_kb: 6n^2 + (m+n)*b -> 2n^2 + (m+n)*b
scholqr3_linops_basic_analytical_kb: m*n + 6n^2 -> m*n + 2n^2
Tests:
- 3 new tests in test_orth_linop.cc cover TRTRI / GEQP3 / BQRRP precond
paths (existing tests only exercised TRSM_IDENTITY).
- All 26 CholQR/sCholQR3/CQRRT tests pass.
PowerOp<InnerOp> (RandLAPACK/linops/rl_power_linop.hh) — generic wrapper representing A^j for a square base linear operator A. Chains j calls to the base op with two ping-pong scratch buffers; A^j is never materialized. j == 1 takes a no-scratch fast path; j == 2 allocates one scratch; j >= 3 allocates two. Op::Trans dispatches base(Op::Trans, ...) j times. Side::Left only (the only consumer pattern we have so far). sparse_axpby_shared_pattern (extras/misc/ext_sparse_axpy.hh) — computes C := alpha*A + beta*B for two CSRMatrix inputs whose sparsity patterns are bit-identical (rowptr + colidxs equal). O(nnz) value-only path. The target consumer is X = K - omega*M in the reduced-spectral application, where K and M from a single FEM mesh always share sparsity exactly. A general-purpose sparse_axpby for different patterns belongs upstream (RandBLAS issue; MKL has mkl_sparse_d_add, cuSPARSE has cusparseDcsrgeam2). Tests: 6 new PowerOp tests (j=1, j=3, multi-RHS, Op::Trans, alpha/beta, PowerOp wrapping CompositeOperator — the rspec usage pattern). 3 new sparse_axpby tests (tridiagonal, shifted-inverse pattern, float type). All 9 pass; full test suite still passes.
Sibling to CholSolverLinOp. Wraps Eigen::SparseLU with COLAMD ordering,
factor-once / solve-many pattern. Handles both SPD and indefinite sparse
matrices — used by the reduced-spectral application when omega is an
interior shift and X = K - omega*M becomes indefinite (Cholesky fails).
Scope:
- In-memory Eigen::SparseMatrix constructor (rspec mode computes X at
runtime via sparse_axpby; the file-based constructor mirroring
CholSolverLinOp's pattern can be added when a consumer needs it).
- operator(): Side::Left, ColMajor, Op::NoTrans on B, Op::NoTrans /
Op::Trans on A. RowMajor and Side::Right deferred.
Tests cover SPD tridiagonal (1D Laplacian), indefinite tridiagonal,
non-symmetric matrix (exercises trans dispatch, A^{-T} != A^{-1}), and
multi-RHS with alpha/beta accumulation. All 4 pass.
TransposedOp<InnerOp> (RandLAPACK/linops/rl_transposed_linop.hh) — implicit
transpose view of any LinearOperator. One-line dispatch wrapper: forwards
to base() with the trans flag flipped. Generic over the LinearOperator
concept so it composes with DenseLinOp, SparseLinOp, CompositeOperator,
PowerOp, and even itself (double-transpose tests pass).
Use case from the rspec application: build C = L^T * X^{-1} * L as
CompositeOperator(TransposedOp(L_op), CompositeOperator(X_inv_op, L_op))
without materializing a separate L^T sparse matrix. But TransposedOp is
intentionally generic — any future 3+ operand chain in our codebase that
needs a transpose-on-one-operand will reuse it.
Also: added the concept-required 12-arg operator() overload (no Side,
delegates to Side::Left) to both PowerOp and TransposedOp. Without this
overload, nesting these wrappers (e.g., PowerOp around PowerOp, or
TransposedOp around PowerOp) fails the LinearOperator concept check.
Other linops (DenseLinOp, SparseLinOp, CompositeOperator, CholSolverLinOp)
already had this overload; the new wrappers now match the convention.
Tests: 5 new TransposedOp tests (dense, double-transpose==identity,
CompositeOperator inner, transposed-of-transposed, around-PowerOp).
All 11 PowerOp + TransposedOp tests pass.
…mark
Implements Algorithm 4 from the collaborator's pseudocode: reduced-basis
Rayleigh-Ritz approximation of eigenvalues of the symmetric operator
C = L^T * (K - omega*M)^{-1} * L, where L L^T = M
on the subspace range(V_app), with V_app = C^j * V_FEM.
New mode "rspec" alongside the existing svd/irlsq/both. Two new positional
CLI args: <omega> (double, default 0.0) and <power_j> (int, default 1,
constrained to {1, 2, 3} per collaborator spec). FEM input only — sparse
mode rejected for rspec.
Operator chain assembly (rl_cqrrt_applications.cc, run_rspec_benchmark):
X = sparse_axpby_shared_pattern(K, -omega*M) // O(nnz)
X_eigen = convert(X) // triplets
X_inv_op = SparseLUSolverLinOp(X_eigen).factorize() // try/catch
L_op = SparseLinOp(L_inv_op.make_L_csc()) // L from chol(M)
C_op = CompositeOperator(TransposedOp(L_op),
CompositeOperator(X_inv_op, L_op)) // L^T X^{-1} L
Cj_op = PowerOp(C_op, power_j)
V_app_op = CompositeOperator(Cj_op, V_op) // implicit, no mat
Per-(algorithm, run): 5-method dispatch (CQRRT_linop, CholQR, sCholQR3,
sCholQR3_basic, CQRRT_linop_bqrrp) reuses the same QR drivers as the
svd/irlsq paths. After QR returns R, Rayleigh-Ritz forms the small n*n
matrix T = R^{-T} * V_app^T * C * V_app * R^{-1} block-by-block (no mxn
intermediate materialization), symmetrizes against rounding drift, then
syevd gives eigenvalues and eigenvectors. Ritz residuals computed for
top-k pairs as ||K v - lambda M v|| / (||K v|| + |lambda| ||M v||) where
v = V_FEM * R^{-1} * u.
Output: <ts>_rspec_results.csv with columns
algorithm, run, m, n, omega, power_j, qr_status, qr_time_us, peak_rss_kb,
analytical_kb, factor_time_us, rspec_total_us, eig_0..eig_{k-1},
resid_0..resid_{k-1}.
Singular-X handling: when omega is too close to an eigenvalue of (K, M),
SparseLU.factorize raises RandLAPACK::Error. Caught at the top of
run_rspec_benchmark; a single failure row with qr_status=-99 is written
and the run returns 0 (not an error — collaborator flagged this as an
expected case).
Fix to PowerOp + TransposedOp: const-correctness on the dense B input
(T* const B -> const T* B) so they compose under CompositeOperator,
which always passes B as const T* through its body.
ext_cholsolver_linop and ext_sparselu_linop minor cleanup.
Build: clean. 36 existing CholQR/sCholQR3/CQRRT/PowerOp/TransposedOp tests
pass. End-to-end smoke test on the small FEM2 problem (75824 x 8304,
omega=0, j=1, method_mask=1) progresses through matrix load + L factor
(310 ms) + X=K-omega*M assembly + SparseLU factor (660 ms) + composite
operator construction. PCholQR warmup is in progress — execution-time
proof that the data flow is correct; the actual benchmark needs the
ISAAC walltime budget (rspec at this size is dominated by the
SparseLU(75824).solve(8304 RHS) at each C-application).
A unit test that compares Ritz eigenvalues against lapack::sygv on a
small synthetic problem is left as a TODO; the FEM smoke test
exercises every new code path.
- CQRRT_linop_applications: remove GSVD post-processing (gesdd on R for
generalized singular values/vectors), upcast-orth diagnostic, and the
"svd"/"both" modes. Available modes are now {irlsq, rspec}. Drops
~410 lines: result-struct SVD fields, A_materialized + AtA_precomputed
setup, do_svd/do_irlsq branching, GSVD CSV writers, skip_svd/upcast_orth
CLI args, and the K_file/V_file params from the inner runner (unused
after the GSVD writer deletion).
- extras/test: delete test_ext_sparse_axpy.cc and drop it from
CMakeLists.txt; the helper it covered is not part of the API surface
we still care about.
Drop the local orth_error wrapper and condition_number reimplementation in CQRRT_diagnostic.cc; their callers now use RandLAPACK::testing::orthogonality_error and RandLAPACK::util::cond_num_check directly. gram_condition_number is kept as a thin convenience that builds G = A^T A then calls cond_num_check.
- Drop the make_eye<T>(n) and fill_lower_from_upper<T> helpers in cqrrt_bench_common.hh. The 6 make_eye callers now allocate locally and call RandLAPACK::util::eye directly; fill_lower_from_upper had no remaining callers. - gram_condition_number: replace gemm A_pre^T A_pre with syrk (upper) + RandBLAS::symmetrize, since the result is symmetric by construction. - Replace the two hand-rolled upper-triangle extraction loops in the GEQP3 and BQRRP paths of CQRRT_diagnostic with lapack::lacpy(Upper).
ext_sparselu_linop.hh
Replace the per-call Eigen::Matrix X(m, n) allocation with a member
T* x_buf_ + int64_t x_buf_size_ buffer, grown on demand. The solve
writes through an Eigen::Map over x_buf_ and the alpha/beta scale-
and-accumulate into C reads from the same buffer. Destructor frees;
copy is now deleted (Eigen members were already non-trivially copyable).
test_ext_sparselu_linop.cc
- Use setFromTriplets instead of reserve/insert for sparse construction.
- apply_A: replace the Eigen::Map dance with a direct InnerIterator
sweep over the sparse matrix.
- Replace the hand-rolled err/ref accumulator in every test with a
shared rel_err() helper that uses blas::nrm2.
Remove the caller-owned G (n×n) and A_temp (m×b_eff) workspace arguments from cholqr_primitive; the primitive now allocates and frees them itself. The three callers (CholQR_linops, sCholQR3_linops, sCholQR3_linops_basic) drop their explicit alloc/free for these two buffers; the timing accumulator alloc_dur for them goes to zero (other per-driver allocations remain). pcholqr_primitive's caller-owned scratch is unchanged.
IR-LSQ (rl_iter_refine_lsq.hh)
Replace the 10 std::vector<T> workspaces in IterRefineLSQ::call (r, g,
c, z, dx, cg_r, cg_p, cg_Mp, tmp_n, tmp_m) with raw T* + a single
free_workspace() cleanup lambda invoked on both return paths.
Applications benchmark (CQRRT_linop_applications.cc)
- Hoist the per-(alg, run) QR-output buffer R, plus the four IR-LSQ
sketch-and-solve buffers (SA, Sb, x_ls, Ax), out of both loops and
into raw T* allocations at function entry; zero-fill R per iter to
match the prior std::vector<T> R(n*n, 0) behavior. Apply the same
pattern to the warmup R_warm and to the rspec R buffer.
- Replace std::vector<T> with raw T* in compute_AtA_blocked (E_block,
A_block, AtA_block), estimate_op_2norm (v, Av), and
compute_orth_error_memlite (X).
Output-state vectors (bench_result lists, top_eigvals/top_residuals,
selected_algs, run_states, ir.times / inner_iters_per_step, and one-shot
RHS construction at function boundaries) keep their std::vector type per
the project rule that small fixed-size scratch / containers are fine.
- compute_Q_from_R helper: Eye buffer → raw T*. - run_algorithms hot loops: Q_uniform (function-level reusable Q buffer), R_rss (CQRRT/CholQR RSS scopes), R_cqrrt/R_cholqr/R_scholqr3/R_dense (per-algorithm per-run QR output), I_mat → raw T* allocated outside the inner loops and zero-filled per iter where needed. Q_uniform freed at function exit. Container vectors (run_states, results, sizes) and breakdown<long> output state retain std::vector per the project convention for small/fixed-size containers.
Two helpers added to the shared header:
- load_csr_verbose: load_csr + the "Loading <label> from <path>... done
(m x n, nnz=N)" progress messages used everywhere we read an .mtx.
- make_run_timestamp: returns the YYYYMMDD_HHMMSS string used to name
CSV output files.
Applications benchmark switches to both helpers. The three sparse-load
sites (sparse-mode A, FEM K, FEM V) and the three timestamp sites (irlsq
CSV, rspec CSV, SparseLU-failure sentinel CSV) collapse to one-liners.
linop_basic and diagnostic kept their inline timestamp idioms: basic uses
a trailing-underscore date_prefix and diagnostic reuses the time_t value
for a ctime() call — neither fits cleanly through the new helper without
restructuring more than this commit warrants.
Two changes to make the CQRRT applications IR-LSQ benchmark numbers
comparable to the April FEM2 plots and to surface failure causes.
CQRRT_linop_applications.cc
Replace compute_orth_error_memlite (X = R^{-T} A^T A R^{-1} ; subtract I)
with compute_orth_error_explicit: materialize Q = A R^{-1} block-by-block
via the linop, then call testing::orthogonality_error<T>(Q, m, n) for
||Q^T Q - I||_F / sqrt(n). The memlite path was mathematically
equivalent but did two TRSMs against an ill-conditioned R, amplifying
forward error by kappa(R)^2 -- giving meaningless ~1e8 numbers on FEM2
(kappa(A) ~ 1e7) even when the algorithm produced a backward-stable LS
solution (Higham residual ~3.8e-13). Matches the OLD path at
composite_applications.cc (e770d15). Extra peak memory: m*n (FEM2 large
~80 GB on top of 46 GB; fits within the ISAAC node).
Drops now-orphan compute_AtA_blocked.
rl_cholqr.hh
Every failure-return path in pcholqr_primitive now prints a tagged
stderr line identifying which sub-step bailed: TRSM_IDENTITY / TRTRI /
GEQP3 / BQRRP diag-zero, lapack::trtri info, BQRRP null state, or
lapack::potrf info on the preconditioned Gram. Same for the potrf in
cholqr_primitive (with the shift_factor echoed). Next ISAAC run's
slurm-err will say e.g. "[pcholqr_primitive] FAIL: lapack::potrf on
preconditioned Gram returned info=4127 (non-PD pivot at column 4127;
preconditioner P stability margin insufficient for kappa(A))",
letting us distinguish iter-2 Cholesky breakdown (likely) from
literal-zero diagonals (unlikely).
…itives Algorithmic core cholqr_primitive / pcholqr_primitive grow two new params (max_retries, shift_growth) with defaults of 0 / 10. On potrf failure the primitive snapshots the Gram, restores it, bumps shift by ×shift_growth, and retries up to max_retries times. shift_factor=0 + max_retries=0 keeps CholQR's vanilla behavior; shift_factor=eps + max_retries=10 implements Oleg's "start at 10^-16·||A||^2, grow on demand" prescription. New: CholQR2_linops Two-pass CholQR via the shared primitives. Iter 1 = cholqr_primitive with shift_factor=eps and adaptive retries; iter 2 = pcholqr_primitive (TRSM_IDENTITY, P=R_1) also with adaptive retries. Returns status 1 (iter-1 exhaust) or 2 (iter-2 exhaust). Companion cholqr2_linops_analytical_kb helper. sCholQR3_linops (blocked) Initial shift_factor_iter1 dropped from 11·eps·n to plain eps. Adaptive retries enabled on all three iterations. The old shift was ~10^5× larger than σ_min²(A) for FEM2, which collapsed iter-2's G_2 to a rank-deficient matrix (eigenvalues 1 - shift/(σ²+shift) → 0) and triggered potrf bail at column ~σ_min's home in the spectrum. Smaller shift keeps G_2 near identity; adaptive growth covers the rare bad case. sCholQR3_linops_basic (storage-efficient variant) Refactored from inline syrk+potrf on a materialized Q_buf to use cholqr_primitive + pcholqr_primitive (block_size=0). Now algorithmically identical to sCholQR3_linops with no blocking; the materialized-Q syrk optimization is gone (storage advantage was an artifact of that). Gains the new diagnostic prints + adaptive retries for free. Analytical_kb formula updated to reflect the new 6·n² + 2·m·n peak. Benchmark dispatch Drops CQRRT_linop_bqrrp (bit 6 / 64). Adds CholQR2 at bit 4 / 16. Default method_mask 79 → 31 (0b11111). CQRRTLinopPrecond::BQRRP option is retained inside pcholqr_primitive for non-benchmark callers.
…ckup CQRRT_linop: was missing the Cholesky Gram G(n*n) and its backup snapshot G_backup(n*n) used by the adaptive-shift retry. Bumped from n*n to 3*n*n in the formula. CholQR_linops: now includes cholqr_primitive's G_backup(n*n) and blocked_preconditioned_gram's I_block(n*b_eff). Bumped from n*n + m*b_eff to 2*n*n + (m+n)*b_eff. New prediction is an upper bound; the OS may not fully commit G_backup's pages during the Gram-loop peak, so observed RSS can be a bit below this number. sCholQR3_linops: was missing the cholqr_primitive's transient G + G_backup + A_temp during iter 1, which is the dominant peak (driver scratches + primitive scratches live concurrently during iter 1). Bumped from 2*n*n + (m+n)*b_eff to 5*n*n + (2m+n)*b_eff. sCholQR3_linops_basic: formula 6*n*n + 2*m*n unchanged; comment updated to note that plugging b_eff = n into sCholQR3_linops's new formula gives exactly this expression, so the old formula was already correct. CQRRT_linop_bqrrp: not touched; the BQRRP option is no longer used in the benchmark dispatch but the helper is retained for non-benchmark callers.
X is SPD for omega < lambda_min(K, M) (omega = 0 / near-zero in this app), so use CholSolverLinOp (Eigen confined to the factorization, solves via RandBLAS TRSM) instead of the Eigen SparseLU operator. Adds an in-memory CSR constructor to CholSolverLinOp; deletes SparseLUSolverLinOp + its test; adds a non-Eigen tridiagonal generator (gen_tridiag_csr) to rl_gen and an Eigen-free in-memory CholSolver test (97/97 extras tests pass).
…ner) Merge pcholqr_primitive into cholqr_primitive with an optional P (nullptr = unpreconditioned); update CholQR/CholQR2/sCholQR3/CQRRT. Use lacpy/transpose_square for the QRCP-inverse triangle/transpose ops, add missing comments, simplify build_R_from_A, and fold the in-memory CholSolver tests into the unified solver test.
… wrappers Factor the CholQR-family iteration into one cholqr_iterate(num_iters, shift_iter1, shift_iter_rest) engine; CholQR(1)/CholQR2(2)/sCholQR3(3)/sCholQR3_basic become thin wrappers that preserve their timing layouts. Factor test-mode Q into materialize_Q_from_R. Use laset for the identity block and lapmr for the jpiv row permutation; single-line BLAS/LAPACK calls; fix reinvented routines (lacpy in rl_cqrrt, RandBLAS::symmetrize in the benchmark) and stale comments across the PR.
New matrix-free linops to build the mu-regularized augmented operator A_hat = [A; mu*I] by composition: - ScaledIdentityOp<T>: matrix-free mu*I_n. - VStackOp<Top,Bot>: vertical concat of any two linops sharing a column dimension. Dense Gram/adjoint paths plus a blocked sketch overload (forms W = A_hat*I_block per output column block and sketches that slice with the full S, so S is never partitioned and no n_rows x d intermediate is formed). Feeding A_hat to any Q-less QR driver yields R = chol(A^T A + mu^2 I) with no driver source changes: the 4 CholQR-family methods Gram it, and CQRRT sketches + Grams it via the blocked overload. CQRRT itself is unchanged. benchmark: new irlsq_reg mode -- regularized augmented preconditioner with independent preconditioner / solve precisions, FEM-faithful kappa column-scaling of V, and a consistent RHS b = A x_true so the backward error ~ u(solve_prec) (kappa-robust) and the forward error ||x-x_true||/||x_true|| ~ u*kappa exposes the precision x kappa interaction. test: test_vstack_linop.cc -- 5 tests incl. R^T R = A^T A + mu^2 I via both the CholQR and CQRRT-through-A_hat paths.
The regularization was mu = mu_factor * u(precond) -- with no ||A|| or size
scaling it is ~mn times too small past kappa(A) ~ 1/sqrt(u), so A^T A + mu^2 I
stays numerically singular at high kappa and even the double/double cell fails
(observed on the 2026-06-15 FEM2 run). Replace with the textbook shifted-
CholeskyQR shift applied through the augmented operator:
mu = mu_factor * ||A||_2 * sqrt((mn + n^2) * u(precond))
so mu^2 is the Fukaya shift (~ ||A||^2 * u * mn), keeping the Gram PD for
kappa(A) up to ~1/u. mu_factor ~ sqrt(11) ~ 3.3 is the principled value; lower it
to reproduce the under-regularized "mu = 10u" regime on purpose. ||A||_2 (already
estimated for the Higham metric) is hoisted above the mu computation and reused.
Also: usage strings now list the irlsq_reg mode.
The A2norm/sqrt(mn)-scaled "Fukaya" shift (e2eb9fb) was the wrong objective: it sizes mu for Q-factor ORTHOGONALITY, but IR-LSQ needs R to be a good right PRECONDITIONER, which requires mu <~ sigma_min. The Fukaya mu (~8e-4*||A|| at FEM2 scale) is orders above sigma_min, so the preconditioned CG matrix M = R^-T A^T A R^-1 had kappa(M) ~ (mu/sigma_min)^2 ~ 1e7+ and CG stalled at the iteration cap (fwd_err ~ 0.9, all methods bit-identical). Restore the collaborator's spec exactly: mu = mu_factor * u(precond), mu_factor=10 => mu = 10u, no ||A|| or size scaling. The augmented operator A_hat = [A; mu*I], its Q-less CholeskyQR R, and R-as-right-preconditioner are unchanged and match the spec. ||A||_2 is still estimated for the Higham metric only.
…olesky Pass the same unbounded eps*trace(G) shift-retry to cholqr_primitive that CholQR/CholQR2 use (unshifted first attempt, max_retries=-1, growth 10), so an ill-conditioned or single-precision Gram is rescued instead of failing at potrf. Double precision is unchanged (shift never activates).
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.
Adds
CQRRT_orth_gap.cc— a standalone benchmark that compares the orthogonality quality of CQRRT_linop and CQRRT_expl on the same input matrix. Five input modes are supported:generate— synthetic sparse matrix with controlled condition number and densityfile— any Matrix Market (.mtx) file from diskcomposite— composite operator from aK.mtx+V.mtxpair (buildsCholSolver(K) ∘ Sparse(V) = L⁻¹V), for generalized LS / FEM-type problemsdiag— step-by-step diagnostic on a Matrix Market file, reimplements CQRRT manually to pinpoint where numerical divergence between the two paths arisesdiag_gen— same step-by-step diagnostic on a synthetic matrixUsage examples