I have completed my review. Here is my assessment.
Adversarial Review: Ridge GCV → cuSOLVER SPD solve (Phase 1)¶
VERDICT: fix-then-ship¶
The actual Phase‑1 deliverable — spd_solve, the cuSOLVER handle, the gating helper, the 3 call‑site swaps, the CMake link/compile‑def edits, and the equivalence test — is correct and faithful to the approved design. I could not find a correctness, ABI, or fold‑isolation bug in any of it.
There is one real problem: an unrelated, out‑of‑scope rewrite of the hot gemm() dispatch path is bundled into the same diff. It isn’t in the approved design, it touches a far broader blast radius than the isolated solve, and it’s not covered by the equivalence test. Split it out (or justify+test it separately) before merging; the Ridge work itself is shippable.
Findings¶
1. [SCOPE] gemm() was rewritten to use static thread_local ReusableDeviceBuffer scratch — not part of this design. cpp/src/core/cuda_dispatch.cpp:1559-1572 (and the dA/dB/dC deref changes at :1597-1603). The approved design (§3) scopes Phase 1 to spd_solve + the cuSOLVER handle + CMake + the 3 dual call‑site swaps + the test. This gemm change is a separate performance optimization (replacing per‑call DevicePtr with grow‑only per‑thread scratch to avoid cudaMalloc/cudaFree device‑sync stalls). It is plausibly correct — dA/dB are fully overwritten by copy_h2d each call, dC is only read when beta != 0 (same contract as the original uninitialized DevicePtr), grow‑only ensure() is sound, and thread_local gives per‑thread isolation — but:
It modifies the primitive used by every CUDA GEMM in libn4m, including the
K = Xc·Xcᵀbuild for the dual design and all PLS moment builds — much wider than the isolatedspd_solve.It is not exercised by the equivalence test (which only calls
spd_solvedirectly), and the design’s optional moment‑sweep parity sanity (§5) wasn’t added.§4.2 of the design explicitly cautioned against introducing reused/persistent device buffers in Phase 1.
Fix: extract this into its own PR with its own benchmark + a
run_moment_sweepparity check (N4M_CUDA_RIDGE_DISABLE=1vs not, scores agree ≤1e‑8). If it must ship together, call it out explicitly and add the parity gate.
2. [STYLE] static thread_local lifetime footnote (only if finding #1 stays). cuda_dispatch.cpp:1564-1566. The scratch destructors run cudaFree at thread/process exit, where the CUDA runtime may already be torn down (cudaErrorCudartUnloading). The return is ignored, so it’s benign, and it mirrors the existing CublasState dtor pattern — but it’s a new lifetime that the per‑call DevicePtr didn’t have, and the buffers stay resident at high‑water mark for the process lifetime. Acceptable as a documented tradeoff; noting for completeness.
3. [OK] spd_solve correctness (cuda_dispatch.cpp:1156-1260). Verified against §3.1 and the §4 traps:
A uploaded as‑is, no transpose — correct (A symmetric ⇒ row‑major image == column‑major image;
uploimmaterial, documented at:1207).B packed column‑major
Bcm[target*n+i]=B[i*q+target](:1200-1204) and X unpacked inverseX[i*q+target]=Xcm[target*n+i](:1242-1245) — both correct.devInfochecked after bothpotrf(:1219-1228) andpotrs(:1234-1238).Return codes exactly per spec:
potrf info>0 → 1(not‑PD),info<0 → 2,potrs info!=0 → 2, andcheck_cusolverthrows → caught by thebad_alloc/exception/...envelope →2(:1248-1259).Workspace sized via
cusolverDnDpotrf_bufferSize(:1210-1213), no assumed size.Overflow/validity guards are thorough: null/zero,
n,q ≤ INT_MAX,mul_overflows(n,n),mul_overflows(n,q), and byte‑size overflow (:1170-1189).No leaks (RAII
DevicePtrfordA/dB/dInfo/dWork,std::vectorstaging). Default‑stream ordering is correct —copy_d2hofinfois a synchronouscudaMemcpyon the default stream, so it waits forpotrf/potrs(no explicit sync needed).
4. [OK] Scope of the swaps. Exactly the 3 dual sites route to solve_dual_spd: predict_ridge_from_dual_design (sweep.cpp:529), ridge_dual_cross_heldout_sse (:577), fit_ridge_from_dual_design (:623). The k×k PᵀW inverse (sweep.cpp:721), the primal moment solve (:913), and the PLS‑prefix inverse (:1403) all remain on host solve_square_qr. ridge.cpp is untouched (confirmed via diff‑stat). No Phase‑2 eigh code was added.
5. [OK] Fallback + gating (sweep.cpp:435-469). solve_dual_spd falls through to host QR on s==1 and s==2 (the if (s==0) return; leaves both other codes falling out of the #if block to solve_square_qr). ridge_cuda_dual_enabled: n<256 → false, N4M_CUDA_RIDGE_DISABLE honored, cuda_runtime_available() checked, #else → false. The λ=0 generic case is handled correctly: potrf reports not‑PD → 1 → host QR, which itself rejects the singular R with N4M_ERR_NUMERICAL_FAILURE — identical to the pre‑change behavior. The [[maybe_unused]] on ridge_cuda_dual_enabled correctly prevents an unused‑function -Werror failure in the non‑CUDA build.
6. [OK] Fold‑clean. spd_solve uses only per‑call DevicePtrs and the passed‑in A/B; no persistent or cross‑fold device buffer. K/Y_work are fold‑local. No fold‑to‑fold state in the solve path. (Caveat: finding #1’s gemm scratch is reused, but it’s fully overwritten each call so it carries no numerical state — not a leakage vector, but it’s the kind of reuse §4.2 flagged.)
7. [OK] ABI. spd_solve is declared only in the internal cuda_dispatch.hpp, has no N4M_API, and the project compiles with -fvisibility=hidden (CMakeLists.txt:77-79, n4m_targets.cmake:111-112) plus the n4m_linux.map version script — so it stays out of the dynamic symbol table, exactly like the existing gemm/pls1_* internals. No cpp/abi/, cpp/include/, bindings/, or catalog/ changes (confirmed via diff‑stat). cuSOLVER handle lifecycle is correct: created after cuBLAS, available gated on both succeeding, cuBLAS rolled back (cublasDestroy_v2 + handle=nullptr) if cusolverDnCreate fails, and destroyed before cuBLAS in the dtor (cuda_dispatch.cpp:47-76).
8. [OK] The test (test_internal_linalg.cpp:304-371). It’s a real test: builds a deterministic SPD A = (GᵀG)/n + I (eigenvalues ≥ 1, well‑conditioned), solves via spd_solve and a host QR reference (qr_square_multi, byte‑for‑byte the same algorithm as the production solve_square_qr), and asserts max abs diff ≤ 1e-8. It self‑skips with a printed message when no GPU is present. The not‑PD path is covered with a diag(1,…,1,0) matrix (a clean zero pivot → potrf info=n>0), asserting spd_solve(...) == 1. The #if defined(N4M_USE_CUDA) guards (helper, test, and invocation) plus the cpp/tests/CMakeLists.txt:94-95 generator‑expression compile‑def correctly make the TU see the macro that n4m_c_static keeps PRIVATE. Note: it tests n=64 < 256 but calls spd_solve directly, bypassing the n<256 gate — correct, so the device path is actually exercised.
9. [STYLE] N4M_CUDA_RIDGE_DISABLE semantics (sweep.cpp:439-442). Any non‑empty value other than a leading '0' disables — so N4M_CUDA_RIDGE_DISABLE=false would also disable CUDA, diverging from the truthy_env idiom used elsewhere (cuda_dispatch.cpp:386). This matches the design’s §3.4 spec verbatim and is defensible for a kill‑switch, so it’s not a bug — just a minor inconsistency worth a one‑line comment if you care about env‑var uniformity.
Bottom line: the Ridge spd_solve change is clean, correct, ABI‑safe, fold‑isolated, and well‑tested. Pull the gemm() thread_local‑scratch optimization out of this diff (finding #1) and the Phase‑1 work ships as designed. The orchestrator’s GPU run should confirm the ≤1e-8 equivalence and the not‑PD ==1 assertions, which my static read indicates will pass.