SAPPHIRE: Preconditioned Stochastic Variance Reduction for Faster Large-Scale Statistical Learning — Scientific Code Review

Paper: Jingruo Sun, Zachary Frangella, Madeleine Udell · published paper · arXiv 2501.15941
Code: udellgroup/sapphire · commit 1802bb15b208

This review ran executable checks that compare the released code against what the paper states. Issues come first, grouped by importance: major. Where the released materials are too incomplete to check a claim, the finding says so in its text; such gaps are reproducibility limitations, not demonstrated errors. Findings whose fix is a direct revision of the paper's TeX or the released code include a Prompt to fix panel with a ready-to-use LLM prompt.

Major Findings

R1 [major] — Paper/Code Algorithm Inconsistency

The paper refreshes a scheduled SSN or NySSN preconditioner before sampling the gradient batch, forming v_k = \widehat{\nabla} L(w_k) - \widehat{\nabla} L(\tilde{w}) + \nabla L(\tilde{w}), and applying that iteration's update.

The released step first forms the gradient and updates the parameters, then calls update_precond. The same-batch SVRG estimator and the implemented last-iterate snapshot agree with the paper to error 0.0; the difference is specifically that a scheduled Hessian is formed at the post-step iterate.

Any reported run whose refresh schedule selects a different Hessian point must be rerun after correcting the order.

Suggested fix: Move scheduled preconditioner refreshes before gradient formation and parameter updating. Define the schedule's indexing convention, add an order-sensitive test, and rerun affected experiments.

Repair the released repository for “SAPPHIRE: Preconditioned Stochastic Variance Reduction for Faster Large-Scale Statistical Learning” (the arXiv identifier was not supplied in the review evidence). In sapphire/src/step/step.py, move each scheduled preconditioner refresh before _basic_step forms the stochastic gradient and applies the parameter update. Preserve the same-batch estimator v_k = grad_hat L(w_k) - grad_hat L(w_tilde) + grad L(w_tilde), the full-gradient snapshot refresh, and the valid last-iterate snapshot option. Inspect sapphire/src/updates/precond_updates.py, sapphire/src/gradient_oracles/variance_reduced_oracles.py, sapphire/src/updates/apply_updates.py, and sapphire/src/opt_states/init_interface.py for the associated schedule and state transitions. Explicitly decide and document whether refresh indices are zero-based or one-based. Add an order-sensitive regression test that records the iterate supplied to the Hessian/preconditioner update and verifies that it is the iterate used by that stochastic update. Then rerun every experiment whose refresh schedule can choose a different Hessian point and regenerate the affected results.

R2 [major] — Paper/Code Preconditioner Inconsistency

The paper constructs NySSN with the stabilization \nu = \sqrt{p}\,\varepsilon(\sigma_{\max}(M)), retains the configured damping \rho, and derives the APG step size from the preconditioner's largest eigenvalue. It gives the APG rule as \alpha = (\lambda_1(P) + \rho)^{-1}.

The released SSN operator is internally consistent, with apply/inverse errors 8.881784197001252e-16 and 1.1102230246251565e-16, but stores \alpha = 0.09523809523809522 instead of 1/\lambda_{\max}(P) = 0.10526315789473684. NySSN has apply/inverse errors 8.881784197001252e-16 and 4.926614671774132e-16, but changes requested \rho = 0.05 to 0.7137582940345772 and stores \alpha = 0.11450864745425206 instead of 0.10642007627429591. Its stabilization is spacing(trace(Core))*trace(Core), not the published rule. The standalone APG recurrence itself agrees with error 0.0.

Thus the released preconditioners apply their stored matrices consistently, but their stabilization, damping, and APG step sizes differ from the paper. The effect on convergence curves or timing cannot be determined without rerunning the experiments.

Suggested fix: Either restore the published NySSN stabilization and configured damping or document the implemented alternatives in the paper. Compute APG's step size from the actual returned preconditioner, add spectral tests, and rerun all SSN and NySSN experiments.

Reconcile the paper and released repository for “SAPPHIRE: Preconditioned Stochastic Variance Reduction for Faster Large-Scale Statistical Learning” (the arXiv identifier was not supplied in the review evidence). Inspect sapphire/src/preconditioners/ssn.py, sapphire/src/preconditioners/nyssn.py, sapphire/configs/precond_configs.py, sapphire/src/numpy_linops.py, and sapphire/src/regularizers/_scaled_prox.py. The published NySSN stabilization is nu=sqrt(p)*eps(sigma_max(M)), while the release uses spacing(trace(Core))*trace(Core). The release also changes requested rho=0.05 to 0.7137582940345772. SSN stores alpha=0.09523809523809522 instead of 1/lambda_max(P)=0.10526315789473684, and NySSN stores alpha=0.11450864745425206 instead of 0.10642007627429591. Decide whether the paper or the released behavior is authoritative. If the paper is authoritative, implement its stabilization, preserve configured rho, and compute alpha from the actual returned preconditioner. If the released behavior is intentional, revise the paper’s Sections 3.2.2–3.3 and Algorithm 3 to specify the replacement stabilization, adaptive damping, and resulting APG step-size rule. Do not alter the APG momentum recurrence, which already agrees numerically. Add spectral identity and configured-damping tests for SSN and NySSN, then rerun every SAPPHIRE-SSN and SAPPHIRE-NySSN experiment and regenerate the affected tables and figures.

R3 [major] — Broken SCAD Objective and Proximal Operator

The paper formulates a scalar regularized empirical-risk objective and reports SCAD least-squares experiments on rna-seq and p53.

In the release, public SCAD evaluation raises TypeError: SCADInterface.evaluate() missing 1 required positional argument: 'step_size', and its proximal call raises TypeError: Regularizer.prox() missing 1 required positional argument: 'step_size'. The low-level evaluator returns shape [4] instead of a coordinate-summed scalar, and scad is excluded from the accepted experiment penalties. By contrast, least-squares, logistic, and MCP value/proximal comparisons agree with errors 0.0, 0.0, 2.7755575615628914e-17, and 0.0.

The released SCAD experiments cannot execute or produce the scalar objective plotted in the paper. Corrected rna-seq and p53 results require a rerun.

Suggested fix: Return a coordinate-summed SCAD value, implement its proximal operator, and expose scad through the experiment configuration. Add boundary and interface tests, then rerun both SCAD datasets.

Repair the released repository for “SAPPHIRE: Preconditioned Stochastic Variance Reduction for Faster Large-Scale Statistical Learning” (the arXiv identifier was not supplied in the review evidence). In sapphire/src/regularizers/regularizer_library/scad.py, sapphire/src/regularizers/regularizer_interface.py, sapphire/src/models/glm_api.py, sapphire/src/models/model_factory.py, sapphire/src/models/linear.py, and experiments/configs.py, make SCAD satisfy the same public scalar-objective and proximal interfaces as the other regularizers. Remove the extra undocumented evaluate argument, sum the coordinate penalties to one scalar, and implement and test the SCAD proximal operator. Match the paper’s penalty exactly: lambda|w| for |w|<=lambda; -(|w|^2-2a lambda|w|+lambda^2)/[2(a-1)] for lambda<|w|<a lambda; and (a+1)lambda^2/2 above a lambda, with lambda>0 and a>2. Explicitly decide and document which branch covers |w|=a lambda because the published display leaves that equality uncovered. Add scad to the experiment configuration using the paper’s disclosed parameters; if numerical parameters are unavailable, publish them rather than inventing them. Add scalar-value, boundary, and proximal tests, then rerun and republish the SCAD results for rna-seq and p53.

R4 [major] — Incomplete and Incorrect Dataset Materialization

The paper benchmarks rna-seq, rcv1, p53, yearmsd, covtype, url, and avazu, with medium scale defined by n < 10^6 and large scale by n \ge 10^6.

The released data path cannot materialize any of the seven protocols. rna-seq and p53 have no configuration; url is configured but not downloaded; and avazu, covtype, rcv1, and yearmsd are downloaded only as .bz2 files while their configurations request decompressed files and decompression is disabled. A module imported by every loader is also absent. For RCV1, the official 677399-row test file is configured as training data and the official 20242-row training file as test data, changing the batch rule from 202 to 2048 and iterations per epoch from 100 to 330.

No reported benchmark can be reproduced through this data path. RCV1 additionally exchanges the populations used for training and test loss, so corrected tables and curves require rerunning all ten dataset–objective combinations.

Suggested fix: Publish a complete data manifest and working loaders for all seven datasets, including decompression and the missing configurations. Restore RCV1's official train/test roles, encode every reported objective assignment, and rerun all ten benchmarks.

Repair the released repository for “SAPPHIRE: Preconditioned Stochastic Variance Reduction for Faster Large-Scale Statistical Learning” (the arXiv identifier was not supplied in the review evidence). Update data_processing/configs/dataset_factory.py, data_processing/configs/libsvm_configs.py, data_processing/loaders/libsvm_loader.py, data_processing/loaders/base_loader.py, experiments/experiment.py, and experiments/configs.py, and add the missing data/data_dirs.py module. Publish a working acquisition manifest, enable verified .bz2 decompression, add identified sources and configurations for rna-seq and p53, and make the downloader acquire url. Restore RCV1’s official roles: 20242 training rows and 677399 test rows, not the reverse. Encode the paper’s Table 4 tuples exactly: rna-seq 640/160/20,530/85.83%; rcv1 20,242/677,399/47,236/0.02%; p53 25,136/6,284/5,408/98.51%; yearmsd 463,715/51,630/90/100%; covtype 464,810/116,201/54/22%; url 1,916,904/479,226/3,231,961/0.01%; and avazu 12,642,186/1,719,304/999,990/0.01%. Encode all paper assignments: l1-logistic on rcv1 and covtype; Lasso on rna-seq and yearmsd; elastic-net logistic on avazu and url; SCAD least squares on rna-seq and p53; and MCP logistic on avazu and url. Add materialization and train/test-role tests, then rerun all ten combinations and regenerate Table 4 and Figures 2–7.

R5 [major] — Inconsistent Gradient-Oracle and Budget Accounting

The paper defines effective passes as the total number of full-gradient evaluations. Medium-scale convex and SCAD runs are capped at 120 seconds and 200 full-gradient evaluations; large-scale convex runs are capped at 300 seconds and 200 evaluations. The paper states no corresponding cap for large-scale MCP runs.

The released training loop checks wall time only after a complete epoch and does not enforce the 200-evaluation cap. SVRG makes two minibatch-gradient calls but records one batch, its full snapshot gradient is not counted, and the counters are not exported with the results. The relative-loss formula (R - R_{\mathrm{best}})/R_{\mathrm{best}} and test loss are retained correctly.

The released results therefore cannot reconstruct the paper's effective-pass curves, and every such panel requires a rerun with complete work accounting.

Suggested fix: Instrument every gradient call, export cumulative work counters, and enforce the stated evaluation and time limits during each epoch. Rerun and republish all effective-pass panels.

Repair the released repository for “SAPPHIRE: Preconditioned Stochastic Variance Reduction for Faster Large-Scale Statistical Learning” (the arXiv identifier was not supplied in the review evidence). Update sapphire/src/_training_loop/training_loop.py, sapphire/src/_training_loop/logger.py, sapphire/src/opt_states/opt_states.py, experiments/experiment.py, experiments/results_containers.py, and plotting/plotting.py so that work is counted by one method-independent component-gradient rule. Count both SVRG minibatch-gradient calls, count every full snapshot gradient, export cumulative work counters through OptimizerResult, and convert them to effective passes as the total number of full-gradient evaluations. Enforce both limits during execution rather than only after an epoch: 120 seconds and 200 full-gradient evaluations for medium-scale convex runs, 120 seconds and 200 evaluations for medium-scale SCAD runs, and 300 seconds and 200 evaluations for large-scale convex runs. Do not invent a large-scale MCP cap. Preserve the implemented relative-loss formula (R-R_best)/R_best and exported test loss. Add counter and stopping-boundary regression tests, then rerun every method and regenerate all effective-pass panels.

R6 [major] — Missing Hyperparameter and Tuning Disclosure

The paper distinguishes the conservative convex-theory choice \eta = O(1/(L_P m)) from the practical choice \eta = O(1/L_P). It also claims robust performance without further data-dependent tuning and says that only a few APG iterations are required.

The release exposes an initial step size of 0.25, an APG horizon of 20, preconditioner batch sizes of 256, damping of 0.001, update frequencies of 20, and NySSN rank 10. It also contains rules for a one-percent batch, a large-data batch of 2048, logistic refreshes every 5 epochs, and SSN damping equal to the inverse batch size. It does not provide the numerical L_P estimates or multipliers, per-dataset SAPPHIRE settings, or a tuning grid or audit; nor does the paper give a numerical meaning for “a few.”

The released materials do not allow the practical step-size order or the no-further-tuning claim to be checked. This is a disclosure gap, not a demonstrated error.

Suggested fix: Publish the numerical L_P estimates and multipliers, complete per-dataset configurations, tuning grid or audit, and the APG horizon used for every reported run.

R7 [major] — Missing End-to-End Reproduction Materials

The paper reports that SAPPHIRE is best on every convex benchmark, reaches machine-precision suboptimality on avazu while Catalyst reaches about 10^{-3} to 10^{-4}, and often converges more than 20 times faster. It also reports superior SCAD results and linear convergence for MCP, with a p53 explanation based on NySSN truncating noisy bad-curvature directions.

None of the six convex or four nonconvex benchmarks can be rerun end to end from the released materials. Raw inputs, paper-matched numerical configurations, and the paper-time execution environment are incomplete. Historical result files and generated PDFs are author outputs rather than independent reruns; the current driver does not expose the reported SCAD and MCP benchmark sets, and no p53 ablation tests the bad-curvature explanation.

The released materials therefore do not allow the performance or causal claims to be checked. This is a reproduction gap, not evidence that the reported comparisons are wrong.

Suggested fix: Publish the raw-input acquisition path, complete paper-matched configurations, and a reproducible execution environment for all ten benchmarks. Regenerate Figures 2–7 and provide the p53 ablation needed to test the bad-curvature explanation.