{"id":"caea4c24-1b45-4fab-a470-600139e2f0d6","arxiv_id":"2608.07017","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":7.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":3,"one_line_summary":"An extended Walk-on-Spheres method builds a reusable operator toolbox to solve general linear and nonlinear elliptic equations without a mesh.","lead":"The authors propose a new Monte Carlo framework that precomputes the inverse Laplacian and harmonic extension on a domain, then reuses these operators to solve variable-coefficient and nonlinear elliptic equations without a mesh. A generalist might read it because it offers a mesh-free, parallel path to hard PDEs, a long-standing bottleneck for probabilistic solvers.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"No convergence or stability proof for the algebraic reduction (3.5): entrywise LOT consistency (Prop. 2.13) does not control the ill-conditioned inversion/quadrature, so the central claim remains unverified.","rationale":"The reader's weakest assumption—that the reduction to (3.5) lacks convergence and stability analysis, especially for the Voronoi quadrature and the inversion of Λ—is exactly the load-bearing gap I identified. The paper proves consistency for the LOT building blocks (Prop. 2.13), but the PDE-solving step requires uniform operator convergence, inf-sup stability, and control of regularization, none of which is provided. The reported condition numbers (σ_max≈1e-2, σ_min≈1e-7) confirm that Λ is nearly singular, and TSVD is used without error analysis. This is not an internal contradiction; it is an unverified central claim. A deterministic reassembly experiment would settle whether the core discretization is convergent, separating quadrature/trial-space errors from Monte Carlo and conditioning effects. Because the paper is a numerical-methods contribution with plausible experiments and reproducible matrix data promised in the appendix, conditional acceptance remains the right call, contingent on the missing error analysis or a definitive numerical convergence study.","tokens_in":36829,"tokens_out":8632,"duration_ms":88918,"concrete_test":"Deterministic reassembly test: For Example 3 (disk), compute reference matrices U(1_{V_j})(x_i), ∇U(1_{V_j})(x_i), H(1_{Γ_j})(x_i), ∇H(1_{Γ_j})(x_i) using a high-accuracy deterministic solver (e.g., fine FEM) for m=n=100, 200, 400, 800, p=4√m. Assemble Λ, Σ, A^λ exactly and solve (3.5) with correct and with TSVD (rcond=1e-3). Compute L2 and H1 errors vs u_ex. If the exact-assembly error does not decrease with m,n, the point-evaluation quadrature / trial-space discretization is inconsistent; if it does decrease but the Monte Carlo result in Table 2 does not track it, the noise/conditioning is the limiting factor. This isolates the missing analysis from sampling error.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The central claim—that the LOT reduces arbitrary divergence-type Dirichlet BVPs to the finite-dimensional system (3.5) whose solution yields an approximate PDE solution—is not supported by the paper's convergence results. Proposition 2.13 proves only pointwise consistency of individual matrix entries in a triple limit (m→∞, then M→∞, then N→∞) for fixed data f,g. It gives no norm convergence of the random matrices A_M,N, B_M,N, H_M,N, K_M,N as operators, no quadrature-error estimate for replacing integrals over W_i by λ(W_i) times center values (Remark 3.3), and no inf-sup or conditioning analysis for Λ. The assembled matrix Λ is a quadratic function of these random matrices, and its inversion is the step that actually produces qbar; entrywise consistency of the building blocks does not imply that the solution of (3.5) is close to the true u. Indeed, Appendix A.1 reports σ_max(Λ)≈10^-2, σ_min(Λ)≈10^-7 for Examples 3–4, so Λ is severely ill-conditioned and the paper resorts to TSVD with ad hoc thresholds rcond∈{10^-3, 10^-2}. No error bound connects the TSVD truncation or the Monte Carlo noise to the final L2/H1 error. Consequently, the statement that the method can solve 'arbitrary Dirichlet boundary value problems' of divergence type is currently a heuristic claim supported by encouraging benchmarks, not a proven numerical method.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper proposes a Walk-on-Spheres-based Laplacian Operator Toolbox (LOT): random matrices A_{M,N}, B_{M,N}, H_{M,N}, K_{M,N} that approximate the inverse Dirichlet Laplacian U, its gradient, the harmonic extension H, and its gradient. Section 2 develops probabilistic representations, Monte Carlo estimators, and consistency results (Propositions 2.7 and 2.13). Section 3 uses the LOT to reduce linear divergence-type Dirichlet problems to the finite-dimensional system (3.5), and Euler-Lagrange/energy minimization problems to the optimization of J_{n,m,p,M,N} in (3.14). Numerical experiments on Laplace/Poisson problems, non-symmetric/anisotropic linear equations, p-Laplace, minimal surface, and semilinear equations are reported. The central claim is that the precomputed LOT transforms arbitrary divergence-type Dirichlet problems into an algebraic system whose solution yields an approximate PDE solution, and for linear operators gives estimators of the Green and elliptic-measure operators.","tokens_in":37219,"tokens_out":6437,"duration_ms":64236,"significance":"The probabilistic representations and LOT consistency proofs are nontrivial and appear correct; the paper supplies full proofs in Appendix A.2, and the numerical benchmark suite is broader than in most WoS papers. If the algebraic reduction were backed by convergence and stability analysis, the method would be a significant contribution: it offers a mesh-free, parallel building block for variable-coefficient and nonlinear elliptic problems without custom stochastic representations. The downloadable matrices and pseudocode strengthen reproducibility. However, the current manuscript does not establish the central approximation claim; the missing step is exactly the passage from consistent LOT matrices to controlled errors in the solution of (3.5) or of the minimization problem (3.14).","major_comments":[{"comment":"The system (3.5) is the load-bearing step of the paper, yet no theorem states that its solution \\bar q produces \\bar u close to the weak solution u. Proposition 2.13 proves only entrywise limits for A_{M,N}, B^{(j)}_{M,N}, H_{M,N}, K^{(j)}_{M,N} in the triple limit m→∞, M→∞, N→∞ for fixed data and fixed output points; it does not provide norm convergence of these matrices, and it says nothing about the spectrum of Λ, which is a quadratic function of the random matrices. Entrywise consistency of the building blocks does not imply that the pseudoinverse solution of (3.5) is close to any q for which Uq+Hg approximates u. A convergence theorem for \\bar q and for \\bar u, under explicit conditions on Λ, is required for the advertised claim.","section":"§3.1, Eq. (3.5)"},{"comment":"The derivation of (3.2) replaces integrals over W_i by λ(W_i) times integrand values at the centers x_i, including f(x_i), b(x_i), c(x_i). The PDE assumptions in the same section only require bounded measurable coefficients and f∈L^2(D); pointwise evaluation is not defined for such data. No quadrature error estimate is given, and the regularity needed to justify the quadrature is not stated as a hypothesis of any theorem. Since (3.2) determines Λ and Σ in (3.5), this unquantified approximation is load-bearing; moreover, the limit n→∞ of the output partition W_n is not analyzed anywhere in the consistency results.","section":"§3.1, Remark 3.3 and Eq. (3.2)"},{"comment":"The reported spectra σ_max(Λ)≈10^-2, σ_min(Λ)≈10^-7 (Example 3) and σ_max≈10^-3, σ_min≈10^-7 (Example 4) show that Λ is severely ill-conditioned. The paper uses TSVD with relative thresholds rcond∈{10^-3,10^-2} chosen per example, but provides no error bound for the effect of truncation or of Monte Carlo noise on \\bar q, nor on the final estimator (3.6) of G_L f and H_L g. Without such a bound, the interpretation of \\bar u as an approximate PDE solution, and of the matrices in (3.6) as estimators of the Green and elliptic-measure operators, remains heuristic.","section":"§A.1, Table 7 and Eq. (3.6)"},{"comment":"For the nonlinear problems the paper solves the minimization of J_{n,m,p,M,N} or its critical-point equations with Adam or p-HiSD; Table 7 lists learning rates, initial guesses, and saddle indices. The reported quantity ∥∇J(\\bar f_opt)∥_{2,V} is the gradient norm of the discretized algebraic objective, not the PDE residual, and no theorem relates stationary points of J_{n,m,p,M,N} to solutions of the Euler-Lagrange equation (3.9) as n,m,M,N grow. This is an additional unquantified layer between the LOT and the claimed approximate solution.","section":"§3.2, Eq. (3.14)"}],"minor_comments":[{"comment":"In the last displayed limit, [K^{(j)}_{M,N} \\bar f]_i should be [K^{(j)}_{M,N} \\tilde g]_i, since K acts on boundary values; as written it conflicts with the definition (2.26).","section":"§2.3, Proposition 2.13"},{"comment":"K^{(j)}_{M,N} is initialized as an n×m zero matrix in both pseudocode listings, but from (2.26) and the subsequent loops it should be n×p; both occurrences should be corrected.","section":"§2.4, Algorithms 4 and 5"},{"comment":"The heading 'T aming the singular terms' contains a spacing typo and should read 'Taming the singular terms'.","section":"§2.1, heading"},{"comment":"Remark 2.15 refers to 'A_{m,n}' while the toolbox definition (2.27) and Algorithms use A_{M,N}; please harmonize the notation.","section":"§2.3, Remark 2.15"}],"recommendation":"major_revision","confidential_remarks":"The stress-test concern lands: the gap between Proposition 2.13 and the algebraic reduction (3.5) is real and load-bearing. I do not see circular reasoning or missing citations; the paper is novel and potentially valuable, but the abstract's claim of solving 'arbitrary Dirichlet boundary value problems' is not supported by the current analysis. With convergence and stability results added, or with claims softened to a numerical framework validated by benchmarks, the paper could be suitable for publication."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"The LOT idea is real: they precompute Monte Carlo estimators of the inverse Dirichlet Laplacian U, its gradient, and the harmonic extension H and its gradient once per domain, then reuse them to set up an algebraic system for an unknown source term. As far as I can tell from the cited prior art, nobody does that. The closest works are recursive fixed-point WoS for variable coefficients and neural-operator training with WoS supervision; this is different. The paper also gives clean probabilistic representations and proves pointwise consistency of the estimators (Prop. 2.7, 2.13). The numerical section is broad: linear, semilinear, p-Laplace, minimal surface, and they recover the known Breuer-McKenna-Plum multiplicity results, which is a good external check.\n\nThe soft spot is exactly where the stress-test puts it: the step from entrywise consistency of LOT matrices to the solution of the finite-dimensional system (3.5) has no convergence or stability analysis. The Voronoi-center quadrature that replaces integrals is unquantified, the matrix Lambda is severely ill-conditioned in the reported examples (sigma_min/sigma_max around 1e-7/1e-2), and the TSVD threshold is chosen ad hoc. Entrywise convergence under the triple limit m, M, N does not imply that solving the random system approximates the PDE solution; that would need norm-convergence or an inf-sup condition plus quadrature error bounds. So the abstract's claim about reducing arbitrary divergence-type problems to a solvable algebraic system is currently a heuristic supported by encouraging benchmarks.\n\nThat said, the paper does not try to hide the gap: Remark 3.3 and Appendix A.1 spell out the quadrature assumption and the TSVD choices. The authors also share the estimated matrices via a link, which is a step toward reproducibility, though there is no code for the PDE assembly.\n\nI would send this to a serious referee. The framework is novel, the probabilistic part is rigorous, and the experiments are wide enough to be worth refereeing. The referee should press for an error analysis for (3.5) or at least a systematic empirical convergence study (error vs. n, m, N and conditioning of Lambda), plus fuller reproducibility notes. If those come back well, this could become a useful mesh-free tool.","headline":"The LOT precomputation is a genuinely new idea for WoS-type solvers, and the experiments are encouraging, but the missing error analysis for the algebraic solve step leaves the 'general solver' claim heuristic.","tokens_in":37683,"tokens_out":2189,"would_cite":true,"duration_ms":24639,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["65N75","35J25","65C05","65N80"],"pacs":[],"model":"deepseek-v4-flash","headline":"The paper proposes a reusable Monte Carlo toolbox that turns divergence-type elliptic Dirichlet problems into small algebraic systems.","keywords":["Walk-on-Spheres","Monte Carlo","elliptic PDE","Dirichlet problem","harmonic extension","Green function","Voronoi","probabilistic numerical method"],"falsifier":"Take a domain with a known solution, such as the Poisson problem on the unit disk, precompute the LOT matrices with large Monte Carlo sample counts, and solve (3.5) for increasing numbers of interior and boundary points. If the approximate solution does not converge to the exact one as the Voronoi cell diameters tend to zero, or if ever stronger truncation of the singular values is needed as the discretization refines, then the point-evaluation quadrature or the invertibility step behind (3.5) is not sound.","tokens_in":36638,"feed_emoji":"🎲","tokens_out":5897,"duration_ms":62197,"temperature":0.7,"pith_summary":"The paper claims that the classical Walk-on-Spheres method can be extended from the Laplace equation to general divergence-type elliptic problems, linear and nonlinear, without deriving a new stochastic representation for each operator. Its strategy is to precompute, once per domain, random matrices that approximate the inverse Dirichlet Laplacian, the harmonic extension operator, and their gradients, forming a Laplacian Operator Toolbox. Any Dirichlet boundary value problem is then approximated by writing the solution as the toolbox applied to an unknown source term plus the harmonic extension of the boundary data, and solving a finite-dimensional algebraic system for that source term. If this reduction is correct, the same precomputed toolbox solves many different PDEs on the same domain, and for linear operators it also supplies estimators of the Green integral operator and the elliptic measure operator.","feed_headline":"One Monte Carlo toolbox solves elliptic Dirichlet problems","feed_subtitle":"Precomputed Laplace and harmonic-extension matrices handle linear and nonlinear cases alike.","key_machinery":"The central object is the Laplacian Operator Toolbox: random matrices that approximate $U$, $\\nabla U$, $H$, and $\\nabla H$ after a Voronoi discretization of the domain and boundary, defined in (2.27). The matrices are Monte Carlo estimates of occupation times and exit locations of Walk-on-Spheres chains; entries of $A_{M,N}$ and $B^{(j)}_{M,N}$ accumulate weighted visits to spatial cells, while entries of $H_{M,N}$ and $K^{(j)}_{M,N}$ record where chains first approach the boundary. The argument then uses the representation $u\\approx Uq+Hg$ to reduce the linear PDE to the algebraic system (3.5), built from the toolbox matrices and the coefficients, and reduces Euler-Lagrange problems to minimizing or finding saddle points of the discrete functional (3.12).","core_discovery":"The central claim is that a precomputed Laplacian Operator Toolbox (LOT), consisting of the random matrices $A_{M,N}$, $B^{(j)}_{M,N}$, $H_{M,N}$, and $K^{(j)}_{M,N}$ that approximate $U=(-\\frac{1}{2}\\Delta_0)^{-1}$, its gradient, the harmonic extension $H$, and its gradient, as collected in (2.27), is consistent (Proposition 2.13) and can be used to solve arbitrary Dirichlet boundary value problems of divergence type. For linear operators, the solution is represented as $u\\approx Uq+Hg$ with an unknown source $q$; testing the PDE against functions $U\\varphi$ and using the toolbox matrices turns the weak form into the linear algebraic system $\\Lambda\\bar q = -A^\\lambda_{M,N}\\tilde f - \\Sigma\\tilde g$ (3.5), whose solution yields the approximate solution at the discretization points. The same formula also gives Monte Carlo estimators of the Green integral operator $G_L$ and the elliptic measure operator $H_L$. For nonlinear problems arising as Euler-Lagrange equations, the toolbox converts the energy functional into a finite-dimensional functional $J_{n,m,p,M,N}$ on $\\mathbb{R}^m$, which is minimized or whose saddle points are searched. The numerical experiments cover non-symmetric and anisotropic linear equations, $p$-Laplace, minimal surface, and semilinear equations with multiple solutions.","pith_inferences":["The paper leaves open a convergence proof for the algebraic reduction itself: Proposition 2.13 proves consistency of the LOT matrices, but not that solutions of (3.5) converge to the PDE solution as the Voronoi cells refine, a gap the ill-conditioning of $\\Lambda$ in the examples makes concrete.","If the reduction step can be rigorously stabilized, the same toolbox could plausibly be applied in higher dimensions, where mesh-based solvers degrade, but whether the conditioning of $\\Lambda$ worsens with dimension is untested.","The Green-operator estimator in (3.6) hints at a possible one-step propagator for parabolic or time-dependent problems, though the paper does not pursue that direction."],"forward_implications":["For any single domain, one precomputed toolbox can be reused to solve many different linear and nonlinear Dirichlet problems without rebuilding a mesh or deriving new exit laws.","For a general linear second-order elliptic operator, the same matrices provide estimators of the Green integral operator and the elliptic measure operator through formula (3.6).","Nonlinear Euler-Lagrange problems, including $p$-Laplace and minimal surface equations, are reduced to finite-dimensional optimization problems whose solutions are then mapped back to approximate PDE solutions.","The method inherits the mesh-free and parallel-friendly character of Walk-on-Spheres, so the toolbox construction scales favorably with dimension compared to mesh-based solvers."],"supporting_citations":[{"why":"Introduces the Walk-on-Spheres algorithm whose exit-from-sphere sampling the entire toolbox is built on.","marker":"[40]"},{"why":"Supplies the elliptic regularity and Green-function facts used to justify the operators $U$ and $H$ and the density of the representing space.","marker":"[24]"},{"why":"Provides the probabilistic representation of the Poisson Dirichlet problem used for $U$ and its gradient.","marker":"[23]"},{"why":"Recent WoS-based probabilistic representations and approximate distance-to-boundary guarantees used in constructing the toolbox.","marker":"[10]"},{"why":"Recursive WoS for variable-coefficient operators, the approach this paper generalizes by precomputing a reusable toolbox.","marker":"[47]"},{"why":"Truncated SVD regularization used to solve the ill-conditioned system matrix $\\Lambda$ in Examples 3 and 4.","marker":"[30]"},{"why":"Semilinear problem with multiple solutions used as a benchmark for the saddle-point search.","marker":"[13]"},{"why":"Preconditioned high-index saddle dynamics used to locate multiple solutions of the semilinear example.","marker":"[32]"},{"why":"Grid-free WoS geometry processing whose nearest-neighbor and Voronoi techniques support the efficient implementation.","marker":"[45]"}],"fun_headline_variants":["Walk-on-spheres toolbox solves diverse elliptic PDEs","One Monte Carlo toolbox for linear and nonlinear elliptic PDEs","Precompute Laplace, then solve any divergence-type elliptic problem","Reusable walk-on-spheres toolbox for elliptic PDEs","Elliptic solver from precomputed Laplace toolbox"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The reduction of the PDE to the algebraic system assumes that evaluating functions at representative points inside Voronoi cells gives accurate integrals, and that the resulting matrix $\\Lambda$ is invertible or can be safely regularized; the paper proves consistency only for the Monte Carlo matrices, not for this reduction, and its examples require truncated SVD because $\\Lambda$ is ill-conditioned.","fun_headline_variants_meta":{"raw":{"variants":["Walk-on-spheres toolbox solves diverse elliptic PDEs","One Monte Carlo toolbox for linear and nonlinear elliptic PDEs","Precompute Laplace, then solve any divergence-type elliptic problem","Reusable walk-on-spheres toolbox for elliptic PDEs","Elliptic solver from precomputed Laplace toolbox"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000657,"raw_usage":{"total_tokens":3101,"prompt_tokens":1132,"completion_tokens":1969,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":748,"completion_tokens_details":{"reasoning_tokens":1891}},"tokens_in":748,"tokens_out":1969,"duration_ms":14596,"temperature":1.0,"reasoning_tokens":1891,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-10T16:33:15.373089+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Take a domain with a known solution, such as the Poisson problem on the unit disk, precompute the LOT matrices with large Monte Carlo sample counts, and solve (3.5) for increasing numbers of interior and boundary points. If the approximate solution does not converge to the exact one as the Voronoi cell diameters tend to zero, or if ever stronger truncation of the singular values is needed as the discretization refines, then the point-evaluation quadrature or the invertibility step behind (3.5) is not sound.","supporting_citations":[{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Introduces the Walk-on-Spheres algorithm whose exit-from-sphere sampling the entire toolbox is built on."},{"cited_title":"Gilbarg and N","cited_arxiv_id":null,"evidence_quote":"Supplies the elliptic regularity and Green-function facts used to justify the operators $U$ and $H$ and the density of the representing space."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Provides the probabilistic representation of the Poisson Dirichlet problem used for $U$ and its gradient."},{"cited_title":"From Monte Carlo to neural networks approximations of boundary value problems","cited_arxiv_id":"2209.01432","evidence_quote":"Recent WoS-based probabilistic representations and approximate distance-to-boundary guarantees used in constructing the toolbox."},{"cited_title":"Sawhney, D","cited_arxiv_id":null,"evidence_quote":"Recursive WoS for variable-coefficient operators, the approach this paper generalizes by precomputing a reusable toolbox."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Truncated SVD regularization used to solve the ill-conditioned system matrix $\\Lambda$ in Examples 3 and 4."},{"cited_title":"Breuer, P","cited_arxiv_id":null,"evidence_quote":"Semilinear problem with multiple solutions used as a benchmark for the saddle-point search."},{"cited_title":"Huang, H","cited_arxiv_id":null,"evidence_quote":"Preconditioned high-index saddle dynamics used to locate multiple solutions of the semilinear example."},{"cited_title":"Sawhney and K","cited_arxiv_id":null,"evidence_quote":"Grid-free WoS geometry processing whose nearest-neighbor and Voronoi techniques support the efficient implementation."}],"review_version":1}