{"id":"ba2c668c-aa20-4735-92df-392f12e51a6b","arxiv_id":"1908.04423","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":5,"one_line_summary":"A Jacobi-Davidson based direct linear method avoids storing Hamiltonian and overlap matrices in variational Monte Carlo optimization, cutting wall time and memory while reaching the same or lower energies on systems up to about 60,000 parameters.","lead":"This paper presents a memory-lean version of the linear method for optimizing non-linear wavefunctions in variational Monte Carlo, using the Jacobi-Davidson eigensolver to avoid building and storing large Hamiltonian and overlap matrices. The resulting accelerated linear method (aLM) reaches energies comparable to exact references on systems with up to tens of thousands of parameters, at a fraction of the wall time of the standard linear method.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"The sampled Hamiltonian is nonsymmetric, but the Jacobi–Davidson correction equation is claimed to be solved by conjugate gradient; the paper never states that H is symmetrized, so the core solver is under-specified as written.","rationale":"The reader's weakest-assumption concern about the correlation length is legitimate for the end-to-end O(N^4) scaling claim, but I judge the unsymmetrized-Hamiltonian and conjugate-gradient issue to be more load-bearing because it affects the correctness of the core solver rather than an asymptotic prefactor. The paper explicitly acknowledges in Section II.B that the sampled H is nonsymmetric, yet Section II.C justifies the Jacobi-Davidson correction solve with an argument that holds for Hermitian A. No symmetrization is stated, so the algorithm as written is incomplete. The numerical agreement with SHCI and GFMC is encouraging evidence that a working variant exists, and the outer-product algebra in Eqs. (24)-(26) is sound, but the results do not identify which solver variant was actually used. This does not overturn the conditional verdict: it reinforces it. The verification step would settle the concern by testing both the unsymmetrized and symmetrized formulations against the published timings and energies.","tokens_in":17972,"tokens_out":18653,"duration_ms":194002,"concrete_test":"Reproduce the Neon or C2 optimization exactly as written: build the unsymmetrized sampled \\bar H from Eq. (25), run Jacobi-Davidson with conjugate gradient for Eq. (30), and record the residual and convergence history. Then repeat with explicit symmetrization \\bar H_sym = (\\bar H + \\bar H^T)/2. If the published energies are reproduced only in the second run, the paper must state that symmetrization is used; if the first run also converges, the authors should explain why conjugate gradient is valid for the nonsymmetric F.","verdict_should_be":"UNCHANGED","load_bearing_attack":"Section II.B states that the Monte Carlo sampled Hamiltonian matrix will not be symmetric when a finite average is taken. Section II.C then sets A to a general square matrix and solves the correction equation F·δ = -r with F = P†(A - ρB)P, asserting that F is positive semi-definite and that only a few conjugate-gradient iterations are needed. For the nonsymmetric sampled H, F is not Hermitian, so conjugate gradient is not mathematically justified and positive semi-definiteness is not guaranteed. The paper does not describe any symmetrization step, such as replacing A by (A + A^T)/2 or averaging G_r and G_c, between Eqs. (24)-(25) and the Jacobi-Davidson iteration. If the implementation silently symmetrizes H, that is a missing algorithmic detail that changes the method and should be reported; if it does not, the correction equation in the paper is not a valid conjugate-gradient system. Since every reported optimization passes through this solve, this is the most load-bearing gap in the central algorithmic claim.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"The manuscript presents an accelerated linear method (aLM) for optimizing non-linearly parametrized variational Monte Carlo wavefunctions. The Hamiltonian and overlap matrices of the linear method are rewritten as sums of outer products of sampled primitives (Eqs. 24-25), so that their action on a vector can be computed in O(ns nvar) time without building and storing the matrices. The generalized eigenvalue problem is then solved with a Jacobi-Davidson iterative solver, and AMSGrad is used far from the minimum before switching to the linear method. The method is applied to real-space Jastrow wavefunctions for Be and Ne, orbital-space Jastrow-GHF wavefunctions for C2, H50, a 98-site Hubbard model, and Fe(II) porphyrin, with up to 57,722 parameters. Energies agree with exact diagonalization, SHCI, and GFMC references. The authors report an empirical solver scaling of O(ns nvar^1.16/nproc) and argue for an overall O(N^4) scaling.","tokens_in":18169,"tokens_out":7848,"duration_ms":84869,"significance":"The central idea of avoiding explicit matrix storage by combining the outer-product representation with an iterative eigensolver is conceptually clean and, if correctly implemented, would remove the main memory bottleneck that limits the linear method to roughly 10^4 parameters. The algebraic identities in Eqs. (24)-(26) are straightforward, and the benchmark results against SHCI and GFMC are credible. The paper also includes a sensible check for stochastic bias using long SGD runs after convergence. These are real strengths. The main uncertainties are the correctness of the iterative solve for the nonsymmetric sampled Hamiltonian and the scope of the scaling claims, which are not fully established by the presented data.","major_comments":[{"comment":"The sampled Hamiltonian is explicitly stated to be nonsymmetric in Section II.B after Eq. (11), yet the Jacobi-Davidson correction equation F = P†(A - rho B)P is claimed to be positive semi-definite and to require only a few conjugate-gradient iterations. For a nonsymmetric A, F is not Hermitian, so conjugate gradient is not mathematically justified and positive semi-definiteness is not defined in the usual sense. The paper does not state any symmetrization of Hbar, such as replacing A by (A + A^T)/2 or averaging Gr and Gc, before entering the Davidson iterations. Because every reported aLM optimization passes through this solve, this missing specification is load-bearing for the central algorithmic claim.","section":"Section II.C, Eqs. (30)-(32)"},{"comment":"The empirical scaling measurement covers only the Jacobi-Davidson solver time, not the end-to-end optimization, while the abstract claims that the method 'improves the scaling of the algorithm' with respect to the number of parameters. Please either provide an end-to-end scaling benchmark or explicitly restrict the claim to the solver. In addition, the O(N^4) scaling in Section III depends on the assumption that the Monte Carlo correlation length grows linearly with system size. That assumption is stated without supporting measurements, and the observed correlation length of about 1.5 for FeP in Section IV.G shows that the correlation length is not universally large; if it grows more slowly than linearly, the number of stored primitives ns and the direct-method cost grow accordingly.","section":"Section IV.E, Fig. 4 and Section III"},{"comment":"The paper's motivation is the memory bottleneck of the linear method, but no memory measurements are reported. Timings alone do not directly demonstrate the memory advantage of aLM over LM. A comparison of peak memory per process for aLM versus LM, at least for the larger systems, would make the central claim concrete and is needed to support the abstract's statement that the method 'greatly lowers the cost'.","section":"Section IV, Tables I-VI"}],"minor_comments":[{"comment":"There are several typos and formatting issues: 'wav efunctions' in the title, 'symmerty' in Section II.A, and 'Lowdin' should be 'Löwdin' in Section IV.D and the references.","section":"Throughout"},{"comment":"The equation A' = V†A V = V†VA is ambiguous because V†VA could be read as (V†V)A; please write V†(VA) or define VA = A V explicitly before the step.","section":"Section II.C, Davidson algorithm step 2"},{"comment":"The comparison energies for Be and Ne are exact diagonalization in aug-cc-pVQZ and aug-cc-pV5Z Gaussian bases, while the VMC calculations use Slater-type bases. The text should state explicitly that these are different basis sets, which explains why the VMC energy can lie below the quoted CI reference.","section":"Tables I and II"},{"comment":"The line-search procedure with slave wavefunctions and Eq. (37) is described, but the statement that 'the median value step size' is used is not completely specified; please clarify how the median is selected among the five candidate step sizes.","section":"Section III, line search"},{"comment":"The y-axis label 'Avg. time/(ns/nproc)' and the caption 'per Monte Carlo sample per process' should be reconciled, and the axis label should state explicitly that the horizontal axis is the number of parameters nvar.","section":"Fig. 4"}],"recommendation":"major_revision","confidential_remarks":"The nonsymmetry issue in Section II.C is likely fixable by either symmetrizing Hbar (a common practice in linear-method implementations) or switching to a non-Hermitian iterative solver such as GMRES; the numerical results would probably be unaffected if symmetrization is already performed silently. The scaling claims need an explicit scoping and additional data but do not appear to require a fundamentally different method. I see no grounds for rejection, provided the authors address the load-bearing specification gap and the scaling evidence."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Bottom line: this is a real contribution. The outer-product rewrite of H and S plus Jacobi-Davidson iterative solution is a sensible way around the O(nvar^2) memory wall of the linear method, and the empirical scaling (solver cost ~ nvar^1.16, with 57k parameters demonstrated) is impressive. The hybrid AMSGrad-to-LM schedule is not brand new—they acknowledge Otis and Neuscamman—but the direct LM solver itself is new and useful. Energies agree with SHCI/GFMC external references, and the SGD bias check on the Hubbard model is a nice touch.\n\nThe algebra in Eqs. (24)–(25) is correct as far as I can tell. The authors also deserve credit for testing on a wide range of systems, including Fe(II) porphyrin, and for being honest about the correlation-length assumption behind the O(N^4) scaling claim.\n\nSoft spots, in order of severity:\n\n1. The Jacobi-Davidson/CG description is under-specified for a nonsymmetric sampled H. Section II.B says H will not be symmetric under finite averaging; Section II.C defines A as a general square matrix, then claims F = P†(A−ρB)P is positive semi-definite and solvable by conjugate gradient. For non-Hermitian A, that is not justified. The paper never states whether they symmetrize the sampled H before solving. If they do (e.g., averaging G_r and G_c), that is a missing algorithmic step that changes the method; if they do not, the correction equation is not a valid CG system. Every reported optimization passes through this solve, so this is a load-bearing reproducibility gap, not a cosmetic one.\n\n2. The scaling plot covers only the Davidson solver, not the whole optimization. The abstract and conclusion imply end-to-end advantage; the timing tables help, but there are no error bars and no comparison with BLM, the closest competing method. A head-to-head benchmark on the same systems would strengthen the positioning.\n\n3. The O(N^4) overall scaling rests on the correlation length growing linearly with N. They state this as an assumption; plausible, but not demonstrated. If it grows more slowly, the memory advantage erodes.\n\n4. No public code or data. For a methods paper, that is a real limitation for independent verification.\n\nNone of this is fatal. The central method is sound and the results are credible. The CG/symmetrization issue is a fixable clarification, not a fundamental flaw—standard practice in the LM literature is to symmetrize, so I suspect the implementation does it, but the paper must say so.\n\nWho is this for? Anyone doing VMC optimization with large parameter sets: quantum chemists with Jastrow-Slater ansatzes and condensed-matter folks with variational wavefunctions. A serious referee will get value from engaging with it. I would send it to peer review, but ask for a revision that addresses the symmetrization, adds error bars, and compares with BLM if possible.","headline":"A genuine algorithmic advance for large-parameter VMC optimization, with a reproducibility gap around the nonsymmetric sampled Hamiltonian that needs fixing before the solver details can be trusted.","tokens_in":18738,"tokens_out":2450,"would_cite":true,"duration_ms":24446,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":[],"pacs":[],"model":"deepseek-v4-flash","headline":"This paper claims that the linear method for optimizing variational Monte Carlo wavefunctions can shed its memory bottleneck by solving its generalized eigenvalue problem iteratively with Jacobi-Davidson, without ever constructing the…","keywords":["variational Monte Carlo","linear method","Jacobi-Davidson","wavefunction optimization","Jastrow factor","AMSGrad","generalized eigenvalue problem","Hubbard model"],"falsifier":"Measure the autocorrelation time of $g_i[n]$, $h_i[n]$, and $E_L[n]$ as the system size grows for a fixed ansatz family; if $n_{\\mathrm{corr}}$ does not grow linearly with $N$, the $O(N^4)$ scaling argument in Section III breaks. Alternatively, run the optimizer on a system with a short, saturating correlation length and check whether per-iteration memory or time departs from the reported $O(n_s n_{\\mathrm{var}}^{1.16}/n_{\\mathrm{proc}})$ scaling.","tokens_in":17682,"feed_emoji":"⚛️","tokens_out":9294,"duration_ms":82777,"temperature":0.7,"pith_summary":"Variational Monte Carlo optimizes wavefunctions that depend non-linearly on parameters, and the linear method is an effective optimizer for this task but has been limited by the need to store the Hamiltonian and overlap matrices at each step. This paper claims that applying the Jacobi-Davidson iterative eigenvector solver to the linear method's generalized eigenvalue problem removes that memory bottleneck, because both matrices can be applied to vectors using stored Monte Carlo samples without ever being constructed. The result is a direct linear method whose per-processor memory and time scale almost linearly with the number of parameters, enabling optimizations of wavefunctions with up to 57,722 parameters. A hybrid schedule that starts with the cheap first-order AMSGrad optimizer and switches to the direct linear method near convergence reaches energies matching high-accuracy benchmarks for systems including the Hubbard model and iron(II) porphyrin.","feed_headline":"Jacobi-Davidson lifts variational Monte Carlo's memory bottleneck","feed_subtitle":"Solving the linear method's eigenproblem iteratively cuts cost and reaches 57,722-parameter wavefunctions.","key_machinery":"The load-bearing identity is the outer-product decomposition of the linear-method matrices: with $\\bar g_i[n]$ and $\\bar h_i[n]$ vanishing for index zero, equations (24) and (25) express $\\bar{\\mathbf{S}}$ and $\\bar{\\mathbf{H}}$ as sums of outer products of the sampled vectors. This turns each matrix-vector product into a few dot products and vector additions. The other half of the machinery is Jacobi-Davidson, an iterative eigensolver that builds a small subspace, solves a projected eigenproblem, and corrects its current best vector by approximately solving a preconditioned equation with the Jacobi orthogonal component correction; the paper uses subspace size 25 with restarts of 5 vectors and loose tolerances for fast progress.","core_discovery":"The central claim is that the linear method's step can be obtained by solving $\\bar{\\mathbf{H}}\\,\\delta p = E\\,\\bar{\\mathbf{S}}\\,\\delta p$ with the Jacobi-Davidson algorithm, using only matrix-vector products. The paper shows that $\\bar{\\mathbf{H}}$ and $\\bar{\\mathbf{S}}$ can be written as sums of outer products of the sampled primitives $g_i[n]=\\langle n|\\Psi_i\\rangle/\\langle n|\\Psi\\rangle$, $h_i[n]=\\langle n|\\hat H|\\Psi_i\\rangle/\\langle n|\\Psi\\rangle$, and the local energy $E_L[n]$, so their action on a vector costs $O(n_s n_{\\mathrm{var}})$ rather than the $O(n_s n_{\\mathrm{var}}^2)$ needed to build the matrices. Storing these primitives instead of the matrices gives a memory per processor of $O(n_s n_{\\mathrm{var}}/n_{\\mathrm{proc}})$, and sampling only every correlation length reduces both memory and time by a factor $n_{\\mathrm{corr}}$. Empirically the solver scales as $O(n_s n_{\\mathrm{var}}^{1.16}/n_{\\mathrm{proc}})$, and the paper reports optimizations of a 98-site Hubbard model with 57,722 parameters, with converged energies matching reference results.","pith_inferences":["If the near-linear solver scaling persists beyond 57,722 parameters, the practical ceiling for this optimization strategy shifts from linear algebra to the cost of computing the primitives $g_i$ and $h_i$, so hand-coded analytic gradients could extend the method further.","The same matrix-free recipe, outer-product-stored samples plus an iterative eigensolver, should transfer to other high-dimensional stochastic optimizations, including neural-network quantum states, where the overlap matrix would otherwise be prohibitively large.","The bias check protocol, running many unbiased SGD steps around the converged minimum and checking the energy histogram, is a reusable diagnostic for any stochastic second-order optimizer.","A direct testable extension is to record the measured correlation length $n_{\\mathrm{corr}}$ for each system; the claimed $O(N^4)$ scaling assumes it grows linearly with system size, and publishing these values would let others verify the scaling claim without reimplementing the code."],"forward_implications":["The linear method can now be applied to wavefunctions with tens of thousands of parameters; the paper demonstrates 57,722 parameters on the 98-site Hubbard model.","Per-iteration wall time for the linear method drops by roughly an order of magnitude compared with explicitly building and diagonalizing the matrices.","Because the matrices are never stored, memory per processor shrinks as parallelism grows, making the optimizer progressively more attractive on many-core machines.","Starting with AMSGrad far from the minimum and switching to the direct linear method near convergence reaches lower energies in comparable or less wall time than either approach alone.","For real-space Jastrow wavefunctions, the direct linear method is reported as the optimizer of choice, since first-order methods converge inconsistently on those systems."],"supporting_citations":[{"why":"Supplies the low-memory prescription for applying Hamiltonian and overlap matrices to a vector in Hilbert-space VMC, which the direct method adopts.","marker":"[6]"},{"why":"Establishes the O(N^4) orbital-space VMC framework whose sampled primitives the present optimizer reuses.","marker":"[8]"},{"why":"Introduces the linear method as a generalized eigenvalue problem in the wavefunction-plus-derivative basis.","marker":"[9]"},{"why":"Gives the uniform-scaling parameter update used to stabilize the linear method's steps.","marker":"[11]"},{"why":"The prior blocked linear method for large parameter sets, the main LM-based baseline the work improves upon.","marker":"[17]"},{"why":"Supplies AMSGrad, the first-order optimizer used far from the minimum.","marker":"[18]"},{"why":"Davidson's iterative method for eigenpairs without matrix construction, the solver family being generalized.","marker":"[48]"},{"why":"Introduces Jacobi-Davidson with the Jacobi orthogonal component correction for generalized eigenproblems.","marker":"[52]"},{"why":"Provides the benchmark Hubbard-model energy used to judge the quality of the 57,722-parameter optimization.","marker":"[64]"},{"why":"Provides the near-exact SHCI reference energies for the Fe(II) porphyrin splitting.","marker":"[66]"}],"fun_headline_variants":["Jacobi-Davidson cuts VMC memory, scales to 57,722 params","Jacobi-Davidson solves VMC eigenproblem without matrices","VMC optimizer skips matrix build via Jacobi-Davidson","Jacobi-Davidson makes VMC optimization memory-efficient"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The claimed $\\mathcal O(N^4)$ overall scaling rests on the assumption that the sampled primitives are serially correlated with a correlation length that grows linearly with system size, so that sampling every $\\mathcal O(N)$ Monte Carlo moves cuts storage and computation by a factor $n_{\\mathrm{corr}}$; if the correlation length grows more slowly, memory and CPU costs grow faster than claimed.","fun_headline_variants_meta":{"raw":{"variants":["Jacobi-Davidson cuts VMC memory, scales to 57,722 params","Jacobi-Davidson solves VMC eigenproblem without matrices","VMC optimizer skips matrix build via Jacobi-Davidson","Jacobi-Davidson makes VMC optimization memory-efficient"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.001343,"raw_usage":{"total_tokens":5500,"prompt_tokens":1034,"completion_tokens":4466,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":650,"completion_tokens_details":{"reasoning_tokens":4390}},"tokens_in":650,"tokens_out":4466,"duration_ms":37705,"temperature":1.0,"reasoning_tokens":4390,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-14T13:43:26.943109+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Measure the autocorrelation time of $g_i[n]$, $h_i[n]$, and $E_L[n]$ as the system size grows for a fixed ansatz family; if $n_{\\mathrm{corr}}$ does not grow linearly with $N$, the $O(N^4)$ scaling argument in Section III breaks. Alternatively, run the optimizer on a system with a short, saturating correlation length and check whether per-iteration memory or time departs from the reported $O(n_s n_{\\mathrm{var}}^{1.16}/n_{\\mathrm{proc}})$ scaling.","supporting_citations":[{"cited_title":"In the Jacobi- Davidson scheme, F = (I − uB · u†)(A − θ · B)(I − u · uB†)","cited_arxiv_id":null,"evidence_quote":"Establishes the O(N^4) orbital-space VMC framework whose sampled primitives the present optimizer reuses."},{"cited_title":"This can be done with a Gram- Schmidt procedure: (a) for v ∈ Cols(V) δ := δ − v · (v† · B · δ) (b) xt = δ/ ‖δ‖, where ‖δ‖ = √ δ† · B · δ","cited_arxiv_id":null,"evidence_quote":"Introduces the linear method as a generalized eigenvalue problem in the wavefunction-plus-derivative basis."},{"cited_title":"To ensure the size of the subspace doesn’t grow too large, restarts should be taken advantage of","cited_arxiv_id":null,"evidence_quote":"Gives the uniform-scaling parameter update used to stabilize the linear method's steps."},{"cited_title":"Sorella, Wave function optimization in the variational monte carlo method, Physical Review B 71, 241103 (2005)","cited_arxiv_id":null,"evidence_quote":"The prior blocked linear method for large parameter sets, the main LM-based baseline the work improves upon."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Supplies AMSGrad, the first-order optimizer used far from the minimum."},{"cited_title":"Neuscamman, Communication: A jastrow factor cou- pled cluster theory for weak and strong electron correla- tion (2013)","cited_arxiv_id":null,"evidence_quote":"Davidson's iterative method for eigenpairs without matrix construction, the solver family being generalized."},{"cited_title":"Tahara and M","cited_arxiv_id":null,"evidence_quote":"Introduces Jacobi-Davidson with the Jacobi orthogonal component correction for generalized eigenproblems."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Provides the benchmark Hubbard-model energy used to judge the quality of the 57,722-parameter optimization."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Provides the near-exact SHCI reference energies for the Fe(II) porphyrin splitting."}],"review_version":1}