{"id":"63dd80c4-3d33-416a-91a3-8cee6b97dced","arxiv_id":"2411.18026","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":7.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":4,"one_line_summary":"A proxy-surface-based fast direct boundary element solver for 2D elastic wave scattering by cavities achieves O(N) complexity in the low-frequency range, with fast multiple right-hand-side solves.","lead":"This paper presents a fast direct solver for two-dimensional elastic wave scattering by a cavity, reaching near-linear time scaling and efficient reuse for many incident waves. It extends the proxy surface method to vector elastodynamic fields and avoids explicit matrix inverses in the Martinsson-Rokhlin scheme.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Multi-level accuracy degrades with N at fixed ID tolerance, so the O(N) claim lacks a demonstrated accuracy-controlled rank bound; a fixed epsilon does not yield a convergent solver as N grows.","rationale":"The paper has genuine strengths: it extends the proxy-surface direct solver to a vector elastodynamic Galerkin BEM, avoids explicit inverses, and reports reproducible-looking scaling data. The reader's conditional verdict is appropriate. My concern sharpens the reader's weakest assumption: Table 1 shows not just missing rank bounds, but actual error growth with N at fixed ε. Single-level results in Table 2 show the single-level compression is stable, isolating the multi-level recursive skeleton selection as the cause. This makes the numerical O(N) scaling in Figure 5 potentially an artifact of insufficient rank at large N, since all timings used ε=10^-8. The proposed test settles this by checking whether accuracy can be controlled without breaking linear scaling. If it cannot, the method is still a useful engineering tool, but the advertised complexity claim would need qualification.","tokens_in":17612,"tokens_out":11507,"duration_ms":104617,"concrete_test":"Run the multi-level FDS for N=12,800, 25,600, and 51,200 at ε=10^-10 and ε=10^-12, recording the relative error versus the conventional BEM and the rank used at each level. If the relative error does not decrease or stay bounded as ε tightens, or if the elapsed time at the tightest ε needed to keep error below a fixed threshold departs from O(N), then the O(N) complexity claim is not accuracy-controlled and the verdict should remain conditional.","verdict_should_be":"UNCHANGED","load_bearing_attack":"Table 1 of Section 8.1 is the most load-bearing evidence problem. At fixed interpolative-decomposition tolerance ε=10^-8, the relative error grows from 8.479×10^-9 (N=400) to 8.404×10^-7 (N=6400), roughly linearly in N; at ε=10^-10 it grows from 2.790×10^-11 to 1.037×10^-8. The authors attribute this to skeletons clustering near cell boundaries in the recursive multi-level compression and state that a more accurate multi-level version is future work. The O(N) timing in Figure 5 was obtained at ε=10^-8, so the timings may not represent an accuracy-controlled solve for large N. The central claim therefore depends on the unproven premise that ranks k=O(1) suffice with error controlled uniformly through all levels. Without a rank bound or an accuracy-controlled scaling test, the numerical O(N) behavior does not establish that the method solves the Galerkin system to a prescribed tolerance at large N.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The manuscript presents a fast direct solver for two-dimensional elastic wave scattering by a cavity in an unbounded elastic solid. The boundary integral equation is the Burton-Miller combination of double-layer and hypersingular operators, discretized by a Galerkin method with piecewise-linear basis functions. The solver is a Martinsson-Rokhlin-type hierarchical direct solver in which shared low-rank factors U_i and V_j for off-diagonal blocks are obtained by interpolative decomposition of proxy-surface interaction matrices, with a variant that avoids explicit inversion of the coefficient matrix and uses fewer matrix-matrix multiplications. Numerical experiments report O(N) elapsed time in the low-frequency range, strong scaling efficiency around 70%, fast multiple-right-hand-side solves, and absence of fictitious eigenfrequencies on the real axis. The paper also includes a comparison with a conventional LU-based BEM and a table of relative errors.","tokens_in":17812,"tokens_out":8417,"duration_ms":75470,"significance":"If the claims are fully substantiated, the solver would be a useful contribution to elastodynamic BEM: it combines the Galerkin Burton-Miller formulation with a linearly scaling direct solver, which is not available in the cited literature for 2D elastic waves. The extension of the proxy surface method to vector wave fields through displacement-direction-wise shared bases is natural and clearly explained. The paper is commendably explicit about the algorithm and about its own limitations, including the multi-level accuracy degradation. However, the current numerical evidence does not by itself establish accuracy-controlled O(N) scalability, because the multi-level error grows with N at fixed truncation tolerance and the timing data were obtained at that same tolerance.","major_comments":[{"comment":"The multi-level fast direct solver shows a systematic accuracy loss as N grows: at ε=10^{-8}, the relative error increases from 8.479×10^{-9} at N=400 to 8.404×10^{-7} at N=6400, and at ε=10^{-10} from 2.790×10^{-11} to 1.037×10^{-8}. The O(N) timing in Figure 5 was generated with ε=10^{-8}, so the large-N timings correspond to solves that are progressively less accurate, not to solves meeting a fixed tolerance. Since the central claim is that the solver is O(N) for solving the Galerkin system, the authors need either (i) a proof or quantitative argument that the interpolative-decomposition ranks remain O(1) with uniformly controlled error through all levels, or (ii) an accuracy-controlled scaling test, e.g., decreasing ε with N or improving the multi-level skeleton selection so that the relative error is held below a prescribed threshold while wall-clock time is measured. The single-level data in Table 2 show that such a test is practicable.","section":"§8.1, Table 1; §8.2, Figure 5"},{"comment":"The Burton-Miller coupling constant α is introduced in Eq. (7) and used throughout the numerical experiments, but no value or rule for choosing it is given anywhere in §8. The behavior of the combined-field formulation, including the relative errors in Table 1 and the absence of fictitious-frequency spikes in Figure 9, can depend strongly on α. Please report the exact value (and, if it varies with ω or the discretization, the formula) used in every experiment; otherwise the numerical results are not reproducible.","section":"§8, numerical setup; Eq. (7)"},{"comment":"The claimed O(N) complexity rests on the assumption that off-diagonal blocks admit a shared low-rank factorization A_{ij}=U_i R_{ij} V_j with rank k=O(1) at every level of the binary tree, but the manuscript provides no bound on k and no error analysis for the recursive compression. The proxy-surface construction in §7 justifies computing U_i and V_j from local interaction matrices, but it does not prove that the resulting rank is bounded as N grows or that the truncation error accumulates benignly through the recursion; the parent-level low-rank claim is deferred to [3]. In view of the observed multi-level accuracy degradation, I ask for either a concrete rank/error estimate for the elastodynamic kernels, or a systematic numerical study reporting the ID ranks and per-level truncation errors as N and ω vary. Without one of these, 'O(N)' is an empirical statement for the tested configurations rather than a property of the method.","section":"§6.2 and §7"}],"minor_comments":[{"comment":"In the display of \\tilde A^p, the lower-left entry is printed as \\tilde A^p_{12}; it should be \\tilde A^p_{21}.","section":"Eq. (39)"},{"comment":"The conclusion refers to the 'Barton–Miller-type boundary integral equation'; this should read 'Burton–Miller'.","section":"§9"},{"comment":"The statement that the 112-core curve 'seems to have a complexity of better than O(N)' as a result of parallelization is confusing. Parallelization changes wall-clock time by a constant factor (ideally) and should not affect the asymptotic complexity order. Please rephrase to separate the measured scaling of the algorithm from the effect of parallel speedup.","section":"§8.2, Figure 5 caption"},{"comment":"The assertion that M_{V_i}^{left} 'has no singularity' and therefore 'can be low-rank approximated' is not a valid justification: a full-rank non-singular matrix need not be low-rank. The actual reason is the smoothness/oscillation of the elastodynamic kernel and the separation of the proxy surface from the cell; the text should say this directly.","section":"§7.1"}],"recommendation":"major_revision","confidential_remarks":"To the editor: The manuscript is within the scope of a numerical analysis/computational engineering journal. The contribution is incremental but solid: it adapts an existing direct-solver framework to 2D elastodynamics with a Galerkin Burton-Miller formulation. The main issue is not correctness of the algebra but the strength of the empirical evidence for accuracy-controlled O(N) scaling; I believe this can be fixed within the scope of a revision, so I recommend major revision rather than rejection. The authors are transparent about limitations, which I found helpful. The paper cites several Japanese-language articles for closely related work; if a broad international audience is intended, a more detailed account of those methods in the introduction would be useful."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Good paper to know about if you work on fast BEM. The genuinely new piece is the extension of the proxy-surface method to 2D elastodynamics: shared direction-wise coefficients for the vector wave field, an inverse-free variant of the Martinsson-Rokhlin sweep, and a Galerkin discretization of the Burton-Miller equation. The derivations look consistent, and the authors are transparent that their sweep algebra follows Gillman et al.'s proof. The multiple-RHS performance is excellent: 180 incident waves at N=12800, with per-vector time from the second vector about 28,900 times less than the first. Strong scaling at 70% is also solid. The single-level version holds accuracy up to N=51200, which is good evidence that the underlying low-rank structure is there.\n\nThe stress-test note lands. Table 1 shows the multi-level solver losing accuracy as N grows at fixed interpolative decomposition tolerance: at eps=1e-8, relative error goes from 8.5e-9 at N=400 to 8.4e-7 at N=6400. The authors attribute this to skeleton clustering near cell boundaries and list a fix as future work. That is honest, but it means the O(N) timing in Figure 5, run at eps=1e-8, is for a solve whose accuracy is not controlled as N grows. The O(N) claim therefore needs qualification: linear scaling for a fixed rank budget, not a proof that the Galerkin system is solved to a prescribed tolerance. A rank bound or an accuracy-controlled scaling test would close the gap. This is the main issue.\n\nTwo smaller points. The Burton-Miller parameter alpha is never specified anywhere; that is a reproducibility gap because alpha affects conditioning and accuracy. And there is no comparison with the H-matrix direct solvers they cite [14,15], which would help calibrate the contribution. No code is shipped, which is common but still limiting.\n\nBottom line: the paper deserves a serious referee. The method is real, the numerics are honestly reported, and the central idea is not wrong. The revision should nail down alpha, add an accuracy-controlled scaling study, and ideally compare against an existing fast elastodynamic solver. I would cite the paper if I were working on fast direct BEM, though the accuracy issue keeps it from being a definitive reference yet.","headline":"Genuine new extension of proxy-surface direct solvers to 2D elastodynamics with honest numerics, but multi-level accuracy degrades with N, so the O(N) claim is not yet accuracy-controlled.","tokens_in":18351,"tokens_out":3721,"would_cite":true,"duration_ms":32396,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["65N38","65F55","74J20"],"pacs":[],"model":"deepseek-v4-flash","headline":"A proxy-surface fast direct solver for 2D elastic wave scattering by a cavity runs in $O(N)$ time in the low-frequency range and avoids fictitious eigenfrequencies on the real axis.","keywords":["elastodynamics","boundary element method","fast direct solver","proxy surface method","Galerkin method","Burton-Miller formulation","low-rank approximation","two-dimensional wave scattering"],"falsifier":"Run the solver at a fixed frequency and geometry while doubling $N$ from, say, 6,400 to 102,400 at a fixed truncation error, and record the rank $k$ chosen by interpolative decomposition at each level; if $k$ grows systematically with $N$, or if the multi-level relative error in Table 1 keeps rising rather than flattening, the constant-rank premise and hence the $O(N)$ claim is falsified.","tokens_in":17394,"feed_emoji":"🌊","tokens_out":8197,"duration_ms":71108,"temperature":0.7,"pith_summary":"Wave scattering by a cavity in an elastic solid leads to dense boundary-element matrices whose direct solution costs $O(N^3)$ in the number of unknowns. This paper claims that a proxy-surface fast direct solver can solve the Galerkin-discretized Burton-Miller boundary integral equation in $O(N)$ time at low frequencies, with no fictitious eigenfrequencies on the real frequency axis. The payoff is practical: after the first right-hand side, additional incident waves become dramatically cheaper, and large two-dimensional scattering simulations become feasible on parallel hardware. If the claimed scaling holds, this gives an asymptotically optimal direct boundary-element method for this class of problems.","feed_headline":"Cavity wave scattering solved in O(N) time","feed_subtitle":"A proxy-surface BEM solver keeps accuracy while handling many incident waves without fictitious eigenfrequencies.","key_machinery":"The machine is the proxy surface method extended to elastodynamics. For each cell of a binary tree decomposition of the boundary, one builds a small local 'proxy' boundary around the cell, discretizes the interaction between the cell and that proxy, and applies interpolative decomposition to the resulting block rows and columns, separately for each displacement direction, to obtain shared coefficients $U_i$ and $V_j$; off-diagonal blocks then factor as $A_{ij}=U_iR_{ij}V_j$ with rank $k$. The fast direct solver alternates an upward step, which compresses the system level by level without explicitly computing $A^{-1}$, and a downward step, which reconstructs the solution from the compressed system using $x_i=A_i^{-1}f_i-A_i^{-1}U_i\\tilde f_i+A_i^{-1}U_i\\tilde A_i y_i$. The Galerkin-discretized Burton-Miller equation supplies the hypersingular operator whose regularization lets piecewise-linear boundary elements work.","core_discovery":"The paper's central claim is that the off-diagonal blocks of the elastodynamic boundary-element matrix admit shared low-rank factorizations $A_{ij}=U_iR_{ij}V_j$, with $U_i$ and $V_j$ obtained by interpolative decomposition on proxy-surface interaction matrices computed separately for the two displacement directions. Feeding these factorizations into a multi-level direct solver that avoids explicitly forming the inverse yields a method whose total work is $O(N)$ in the low-frequency range. Because the solver is built on a Galerkin-discretized Burton-Miller equation, it is free of spurious resonances on the real frequency axis, so frequency sweeps do not hit nonphysical spikes. Numerical results show single-core times consistent with $O(N)$, strong scaling efficiency near 70% on 112 cores, and per-vector times for the 179 additional right-hand sides about 28,900 times smaller than the first solve.","pith_inferences":["Editorial inference: the $O(N)$ claim rests on the rank of the proxy-surface factorizations staying bounded as $N$ grows; the paper's multi-level accuracy degrades with increasing $N$ (its Table 1), which suggests the current interpolative-decomposition skeleton selection does not yet fully enforce this bound.","Editorial extension: a distance-aware cell-size rule or a skeleton-selection criterion that avoids clustering at cell boundaries could plausibly restore the multi-level accuracy that the single-level version already shows.","Editorial extension: comparing the solver against a fast-multipole iterative BEM on the same cavity problems, including setup cost and many right-hand sides, would identify the problem sizes and condition numbers where this direct solver is the faster choice; the paper compares only with conventional dense BEM."],"forward_implications":["With $O(N)$ scaling, two-dimensional elastic wave scattering problems with tens of thousands of boundary nodes become solvable directly, where conventional dense BEM would require $O(N^3)$ work.","Because the solver is direct, it avoids the convergence problems of iterative Krylov solvers for ill-conditioned systems; the paper demonstrates this advantage for problems with many right-hand sides.","Because the formulation is Burton-Miller, scans over real angular frequency show no fictitious eigenfrequency spikes, unlike the non-Burton-Miller BEM compared in the paper.","The upward step is independent of the right-hand side, so once it is done, additional incident waves cost very little; the paper reports about 28,900 times faster per-vector time for the 179 subsequent right-hand sides."],"supporting_citations":[{"why":"Supplies the original proxy-surface fast direct solver that this paper modifies and whose $O(N)$ complexity target it inherits.","marker":"[3]"},{"why":"Extends the proxy surface method to a Helmholtz transmission problem; the elastodynamic extension here follows this template.","marker":"[11]"},{"why":"Contains the equivalent formulas for a direct solver that avoids explicit inversion; the paper's variant uses these formulas in the numerical method.","marker":"[23]"},{"why":"Introduces the Burton-Miller boundary integral equation that removes fictitious eigenfrequencies; the paper's formulation is based on it.","marker":"[25]"},{"why":"Provides the Galerkin and regularization treatment for hypersingular kernels that allows piecewise-linear boundary elements in this method.","marker":"[22]"},{"why":"Defines interpolative decomposition, the linear-algebra tool used to extract shared low-rank coefficients from the proxy-surface interaction matrices.","marker":"[28]"},{"why":"Represents prior H-matrix elastodynamic solvers whose boundary integral formulation retains fictitious eigenfrequencies, motivating the Burton-Miller choice here.","marker":"[14]"},{"why":"Represents another prior H-matrix elastodynamic solver with fictitious eigenfrequencies, further motivating the need for a fast direct solver with the Burton-Miller formulation.","marker":"[15]"}],"fun_headline_variants":["O(N) elastic wave scattering via proxy-surface solver","Proxy surfaces make elastic cavity scattering O(N)","Fast direct solver cracks 2D elastic wave scattering","Elastic scattering solver hits O(N) with proxy surfaces","Linearly scalable solver for elastic wave cavity scattering"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The load-bearing premise is that the interaction between well-separated pieces of the boundary can always be captured by a fixed small set of proxy modes, so the compressed block rank stays bounded as the number of unknowns grows.","fun_headline_variants_meta":{"raw":{"variants":["O(N) elastic wave scattering via proxy-surface solver","Proxy surfaces make elastic cavity scattering O(N)","Fast direct solver cracks 2D elastic wave scattering","Elastic scattering solver hits O(N) with proxy surfaces","Linearly scalable solver for elastic wave cavity scattering"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.00019,"raw_usage":{"total_tokens":1311,"prompt_tokens":889,"completion_tokens":422,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":505,"completion_tokens_details":{"reasoning_tokens":348}},"tokens_in":505,"tokens_out":422,"duration_ms":4024,"temperature":1.0,"reasoning_tokens":348,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-12T11:35:58.148156+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run the solver at a fixed frequency and geometry while doubling $N$ from, say, 6,400 to 102,400 at a fixed truncation error, and record the rank $k$ chosen by interpolative decomposition at each level; if $k$ grows systematically with $N$, or if the multi-level relative error in Table 1 keeps rising rather than flattening, the constant-rank premise and hence the $O(N)$ claim is falsified.","supporting_citations":[{"cited_title":"A fast direct solver for boundary in- tegral equations in two dimensions,","cited_arxiv_id":null,"evidence_quote":"Supplies the original proxy-surface fast direct solver that this paper modifies and whose $O(N)$ complexity target it inherits."},{"cited_title":"A fast direct solver using the Burton- Miller formulation for one-periodic transmission problems of Helmholt’ equation in 2D,","cited_arxiv_id":null,"evidence_quote":"Extends the proxy surface method to a Helmholtz transmission problem; the elastodynamic extension here follows this template."},{"cited_title":"A direct solver with O(N) complexity for integral equations on one-dimensional domains,","cited_arxiv_id":null,"evidence_quote":"Contains the equivalent formulas for a direct solver that avoids explicit inversion; the paper's variant uses these formulas in the numerical method."},{"cited_title":"The application of integral equation methods to the numerical solution of some exterior boundary-value problems,","cited_arxiv_id":null,"evidence_quote":"Introduces the Burton-Miller boundary integral equation that removes fictitious eigenfrequencies; the paper's formulation is based on it."},{"cited_title":"Application of fast multipole galerkin boundary integral equation method to elastostatic crack problems in 3d,","cited_arxiv_id":null,"evidence_quote":"Provides the Galerkin and regularization treatment for hypersingular kernels that allows piecewise-linear boundary elements in this method."},{"cited_title":"On the compression of low rank matrices,","cited_arxiv_id":null,"evidence_quote":"Defines interpolative decomposition, the linear-algebra tool used to extract shared low-rank coefficients from the proxy-surface interaction matrices."},{"cited_title":"Theory and implementation of H-matrix based iterative and direct solvers for helmholtz and elastody- namic oscillatory kernels,","cited_arxiv_id":null,"evidence_quote":"Represents prior H-matrix elastodynamic solvers whose boundary integral formulation retains fictitious eigenfrequencies, motivating the Burton-Miller choice here."},{"cited_title":"Improvement of hierarchical matrices for 3d elastodynamic problems with a complex wavenumber,","cited_arxiv_id":null,"evidence_quote":"Represents another prior H-matrix elastodynamic solver with fictitious eigenfrequencies, further motivating the need for a fast direct solver with the Burton-Miller formulation."}],"review_version":1}