{"id":"d99cdb94-e7a9-464f-acce-35c8e620411a","arxiv_id":"2509.15744","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":7.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":1,"one_line_summary":"A superposition of forward and adjoint wavefields with a scaling factor lets the authors approximate the optimization sensitivity using O(number of grid points) memory instead of O(steps times grid points).","lead":"This paper introduces an approximate adjoint-gradient method that stores only a few wave snapshots instead of the entire wave history, cutting memory use for dynamic optimization dramatically. The approach is demonstrated on full waveform inversion and acoustic topology optimization with up to hundreds of millions of unknowns on a single GPU.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Algorithm 1 omits the adjoint-source time history; for volume-supported sources such as TATO the O(dofs) memory claim fails.","rationale":"The reader's weakest assumption is exactly the load-bearing gap I identify: Algorithm 1 needs the adjoint-source time history in the backward pass, and the paper's memory ledger ignores it. The bilinear expansion in Equation 38 is algebraically sound, and the empirical k-sweep in Section 4.1 gives reasonable evidence that the approximation can be calibrated for the presented examples. However, the core memory claim is not merely an implementation detail: for the TATO objective, the adjoint source is supported on a volume Omega_s, so its time history cannot be stored in O(dofs) memory unless Omega_s has fixed, N-independent size and the asymptotic limit is carefully redefined. The paper neither states a storage scheme nor a regeneration procedure, so the central assertion is currently unsupported for volume-supported adjoint sources. Because the same concern was already the basis of the reader's CONDITIONAL verdict, my read does not move the verdict; it confirms that the condition (clarify or restrict the memory accounting) is necessary before the method can be accepted as a general replacement for stored-forward-wavefield adjoints.","tokens_in":26700,"tokens_out":7479,"duration_ms":70819,"concrete_test":"Instrument the released CUDA/CuPy code (Zenodo record [51]) for the TATO case in Section 5.2 to log peak GPU memory and to trace how f_dagger^n is supplied at Algorithm 1 line 13. If f_dagger over Omega_s is stored across all N steps, or u is recomputed, record the measured memory and compare it with four dofs grids; a term proportional to N*|Omega_s| falsifies the O(dofs) ledger. If the code stores no f_dagger history, identify the mechanism that regenerates u on Omega_s during the backward pass.","verdict_should_be":"UNCHANGED","load_bearing_attack":"Algorithm 1's memory ledger in Section 4 counts only the three solution grids and the kernel, but the backward loop (lines 12-16) requires the adjoint source f_dagger^n for every time step n. During the forward loop, f_dagger^{n+1} is computed at line 5 from the current u^{n+1}; the only scaling step (line 10) sits outside the loop, and no array for the f_dagger history is allocated or described. For FWI (Equation 12), f_dagger is supported on receiver points, so the history is a cheap N_r by N time series. For TATO (Equation 22), f_dagger = -2u/A on Omega_s, a volume-supported quantity that depends on the forward solution at every time. Since u is deleted as the forward loop proceeds, implementing line 13 for TATO requires either storing u (or f_dagger) over Omega_s for all N steps, memory proportional to |Omega_s|*N rather than dofs, or recomputing u, which reintroduces the storage/checkpointing cost the method claims to avoid. The paper does not state which is done. If Omega_s is a fixed fraction of the domain, the asymptotic O(dofs) claim and the derived 2.5e9-parameter limit do not hold for TATO; at minimum the claim must be restricted to point-supported adjoint sources unless a regeneration strategy is supplied.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The manuscript proposes an approximate adjoint sensitivity computation for transient, self-adjoint, undamped wave problems. The central idea is a linear superposition u_s = u + k u_dagger; because the Frechet kernel is bilinear, the identity K(u_s,u_s) = K(u,u) + 2k K(u,u_dagger) + k^2 K(u_dagger,u_dagger) yields the approximation K(u,u_dagger) roughly equal to (K(u_s,u_s) - K(u,u))/(2k), provided the k^2 term is negligible (Eqs. (38)-(39)). Since K(u,u) can be accumulated during the forward simulation and K(u_s,u_s) during a backward simulation of the superposed problem, the forward wavefield need not be stored. The paper presents this as Algorithm 1, studies the sensitivity error as a function of the hyperparameter k, and demonstrates the method on three-dimensional full waveform inversion (up to about 1.2e8 degrees of freedom) and two-dimensional transient acoustic topology optimization (up to about 1.5e7 degrees of freedom), with sensitivity computations reported up to about 1e9 degrees of freedom. The authors also provide a CUDA finite-difference solver and state that the implementation is available online.","tokens_in":27157,"tokens_out":12099,"duration_ms":97575,"significance":"If the central memory claim holds, the method would remove the dominant memory bottleneck in adjoint-based transient optimization on GPUs and would be of practical value in full waveform inversion and topology optimization. The bilinear expansion is exact and the approximation error is characterized honestly through the hockey-stick plots in Fig. 5, including a comparison against standard adjoint gradients over a wide range of k. The released code and the documented solver timings are strengths. However, the claimed asymptotic O(dofs) memory is not established for volume-supported adjoint sources, and the printed Algorithm 1 appears to omit the 1/(2k) normalization required by Eq. (39). These issues are central to the paper's contribution and must be resolved before the memory and scalability claims can be accepted.","major_comments":[{"comment":"Algorithm 1's backward loop (lines 12-16) requires the adjoint source f_dagger at every time step, but f_dagger is computed in the forward loop (line 5) and no storage or regeneration of its time history is specified. The memory discussion in Section 4 counts only 'at most four solution vectors' (u^{n-1}/u_s^{n+1}, u^n/u_s^n, u^{n+1}/u_s^{n+1}, K~gamma), which omits this history. For FWI the adjoint source (Eq. (12)) is supported on receiver points, so storing the residual traces adds only N_r times N floats and the O(dofs) claim can stand if this is stated explicitly. For TATO the adjoint source (Eq. (22)) is supported on the volume Omega_s and depends on the forward solution: implementing line 13 requires either storing u (or f_dagger) on Omega_s for all N time steps, i.e., memory of size |Omega_s| times N, which scales as O(dofs times N) when Omega_s is a fixed fraction of the domain. No such storage is described, and recomputing u would reintroduce the checkpointing cost the method claims to avoid. The paper therefore does not support the 'memory usage is independent of the number of time steps' claim (Section 4.1) or the 2.5e9-parameter limit (Section 4.2) for TATO; the claims must be qualified or a valid O(dofs) strategy for volume-supported adjoint sources must be supplied.","section":"Section 4, Algorithm 1 (lines 5-16), Eq. (22)"},{"comment":"Algorithm 1 accumulates K~gamma as -sum incrementKernel(u) + sum incrementKernel(u_s) (lines 7 and 14) and returns this sum directly. However, Eq. (39) defines K~gamma = (1/(2k))(K(u_s,u_s) - K(u,u)). Because line 10 scales f_dagger by k, u_s = u + k u_dagger and K(u_s,u_s) - K(u,u) equals 2k K(u,u_dagger) + k^2 K(u_dagger,u_dagger); without the factor 1/(2k) the returned quantity is not the stated gradient approximation. Unless the implementation applies an unstated normalization, the error curves in Fig. 5 could not match the reference gradients. Please correct Algorithm 1 (for example, multiply K~gamma by 1/(2k) before the return) or revise Eq. (39) to match the implementation.","section":"Section 4, Algorithm 1 vs. Eq. (39)"},{"comment":"The pseudocode has an inconsistent time index and scaling for the adjoint force: line 5 computes f_dagger^{n+1}, line 13 uses f_dagger^n + f^n, and line 10 multiplies only the final f_dagger^{n+1} by k after the forward loop has ended. If the f_dagger history is kept, the scaling has to be applied to every time level (or inside the backward loop); if it is not kept, line 13 cannot be evaluated. Please specify the intended indexing and move the scaling to a well-defined location.","section":"Section 4, Algorithm 1 (lines 5, 10, 13)"}],"minor_comments":[{"comment":"The list 'at most four solution vectors are needed (u^{n-1}/u_s^{n+1}, u^n/u_s^n, u^{n+1}/u_s^{n+1}, K~gamma)' is garbled: it appears to list the forward and superposed grids at two different times in the same slots; please clarify which three solution grids are reused.","section":"Section 4, paragraph after Algorithm 1"},{"comment":"In the sentence preceding Eqs. (30) and (31), 'The are straightforward' should read 'They are straightforward'.","section":"Section 3.1, before Eqs. (30)-(31)"},{"comment":"The phrase 'In contrast to Figure 18' appears to be a self-reference; it should likely refer to Figure 16.","section":"Appendix B, Figure 18 caption"},{"comment":"The statement that K(u_s,u_s) is 'independent of u' is too strong; u_s is driven by a source that depends on u through f_dagger, so the intended meaning is that K(u_s,u_s) can be accumulated without storing the forward history of u.","section":"Section 4, near Eq. (39)"}],"recommendation":"major_revision","confidential_remarks":"The paper is a promising practical contribution, and the code release is a definite strength. The two blocking issues are the unaccounted adjoint-source time history for volume-supported sources (TATO) and the missing 1/(2k) normalization in the printed Algorithm 1. Both are fixable by revision, but if the O(dofs) memory claim cannot be restored for TATO, the title and abstract should be narrowed to point-supported adjoint sources or the additional |Omega_s| times N storage should be reported transparently."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Two things up front. The superposition identity behind the method is exact, and as a way to avoid storing the forward wavefield it is genuinely new compared with checkpointing, compression, or random boundaries. And the paper's memory claim is not fully supported: the backward loop needs the adjoint-source history, and the paper doesn't say where it lives.\n\nThe core trick is simple and easy to state: for a bilinear kernel K, K(u,u†) is approximated by (K(u+k u†, u+k u†) − K(u,u))/(2k), with the k²K(u†,u†) term dropped. That gives an iterative gradient computation with roughly three solution grids plus the gradient field. The authors implement this in CUDA, provide code on Zenodo, and demonstrate it on full waveform inversion and transient acoustic topology optimization. The 3D FWI at about 1.2e8 degrees of freedom and TATO at about 1.5e7 degrees of freedom are real scale-ups, and the error-versus-k hockey-stick curves are informative and honest. They also state the restrictions up front: self-adjoint problems, no damping, time-reversible integration.\n\nThe main soft spot is the memory ledger. Algorithm 1 computes f† during the forward pass and uses it during the backward pass, but no storage or regeneration of that history is described. For FWI with point receivers, the history is just N_r × N floats, cheap. For TATO, f† = −2u/A on the observation volume Ω_s; that is a volume-supported time series. Unless Ω_s is tiny or checkpointing is used, the O(dofs) memory claim and the derived 2.5e9-parameter limit do not follow. The line that scales the adjoint force outside the loop is also sloppy pseudocode; presumably it means scaling the whole history. This is fixable, but it needs to be stated clearly.\n\nSecond soft spot: dropping k²K(u†,u†) has no error bound. The empirical plateau across many orders of magnitude in k is encouraging, and the calibration procedure (decrease k until single- and double-precision results diverge) is reasonable. Still, a formal estimate would make the method more than a heuristic. Also, the standalone 3D FWI case did not converge in ten iterations; the authors attribute that to FWI well-posedness, which is plausible, but it does limit what that particular demo shows.\n\nWho is this for? People doing transient adjoint-based optimization or FWI on GPUs. I would send it to review. The core idea deserves a serious referee, and the memory-history issue should be caught and fixed in revision. I would want to see a revised Section 4 before relying on the billion-parameter claim.","headline":"A genuinely new superposition trick for adjoint gradients with real GPU scalability, but the O(dofs) memory claim misses the adjoint-source history for volume-supported sources.","tokens_in":27483,"tokens_out":3807,"would_cite":true,"duration_ms":34622,"reading_group":"yes","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 a superposition identity lets an adjoint sensitivity be approximated from two self-terms of a combined forward-plus-adjoint field, so dynamic wave optimization no longer needs to store the forward time history.","keywords":["adjoint optimization","GPU acceleration","finite difference method","dynamic optimization","full waveform inversion","topology optimization","acoustics","memory-efficient sensitivity"],"falsifier":"Run Algorithm 1 with a volume-supported adjoint source that depends on the forward field over the whole domain (a TATO problem with $\\Omega_s = \\Omega$). If the method's peak memory stays at a few solution grids regardless of the number of time steps $N$, the O(dofs) claim holds; if the implementation must store or regenerate the forward field over $\\Omega_s$ for every time step and peak memory grows with $N$, the claim is falsified for that class.","tokens_in":26527,"feed_emoji":"🌊","tokens_out":10199,"duration_ms":81672,"temperature":0.7,"pith_summary":"This paper claims that the gradient of a dynamic optimization objective in wave problems can be computed without ever storing the forward wavefield, by running one forward solve and one backward solve of a superposed wavefield and combining their kernel values. The memory needed drops from the number of grid points times the number of time steps to a few grids of size equal to the number of degrees of freedom, which is what makes billion-parameter problems fit on a single GPU. The price is that the gradient is approximate: the argument relies on the quadratic term in the adjoint field being negligible after choosing a scale factor $k$, and it only works for self-adjoint, damping-free problems with time-reversible integrators. The authors demonstrate the method on full waveform inversion and transient acoustic topology optimization, computing sensitivities for problems up to about $10^9$ unknowns. If the approximation holds, dynamic adjoint optimization on GPUs stops being memory-bound for a large class of wave problems.","feed_headline":"Superposition trick cuts wave-optimization memory to a few snapshots","feed_subtitle":"Sensitivities for wave inversion and acoustic design now need only a few solution grids, not the whole time history.","key_machinery":"The load-bearing object is the bilinear Fr\\'echet kernel $K_\\gamma(u,u^\\dagger)$ — the time-integrated product of forward and adjoint wave fields (velocity and gradient terms) that gives the sensitivity of the cost function to material perturbations. The paper's mechanism is the superposition identity (Equation 39): because the kernel is bilinear, the mixed term $K_\\gamma(u,u^\\dagger)$ is recovered as $\\frac{1}{2k}(K_\\gamma(u_s,u_s)-K_\\gamma(u,u))$ where $u_s = u + k u^\\dagger$, with the $k^2 K_\\gamma(u^\\dagger,u^\\dagger)$ term neglected. This converts the standard adjoint computation, which needs both fields at the same physical time and therefore a stored forward history, into two sequential simulation sweeps — a forward sweep accumulating $K_\\gamma(u,u)$ and an adjoint-source term, followed by a backward sweep of the superposed problem accumulating $K_\\gamma(u_s,u_s)$ — with only a small constant number of solution grids alive at any moment.","core_discovery":"The central discovery is that for a bilinear sensitivity kernel $K_\\gamma$, the mixed term needed for the gradient can be extracted from two self-terms by superposition: with $u_s = u + k u^\\dagger$, bilinearity gives $K_\\gamma(u_s,u_s) = K_\\gamma(u,u) + 2k K_\\gamma(u,u^\\dagger) + k^2 K_\\gamma(u^\\dagger,u^\\dagger)$, and when the adjoint field is much smaller than the forward field the last term is negligible. Hence $K_\\gamma(u,u^\\dagger) \\approx \\frac{1}{2k}\\big(K_\\gamma(u_s,u_s)-K_\\gamma(u,u)\\big)$. The term $K_\\gamma(u,u)$ can be accumulated during the forward solve, and $K_\\gamma(u_s,u_s)$ during a backwards solve of the superposed problem, so the full forward history never has to be stored; only a constant number of solution grids is kept. The cost is an approximate sensitivity with a hyperparameter $k$ whose admissible range spans many orders of magnitude, and the applicability is restricted to self-adjoint problems with time-reversible time integration. The authors demonstrate sensitivity computation up to roughly $10^9$ degrees of freedom on a single 40 GB GPU, and run full waveform inversion and transient acoustic topology optimization at problem sizes that standard adjoint methods cannot reach.","pith_inferences":["For point-sensor adjoint sources, as in the full waveform inversion examples, the required time history of the adjoint source is only the sensor traces, which is cheap to store; the O(dofs) memory claim is therefore solid there. For volume-supported adjoint sources that depend on the forward field over a large region, the time history of that source may reintroduce a memory term proportional to su","The same superposition identity should extend to other bilinear kernels over reversible dynamics, such as elastic wave equations without attenuation or Schr\\\"odinger-type equations, since the argument only requires self-adjointness and time reversibility.","Because the optimizer used in the paper is Adam, which tolerates gradient noise, the approximation's dependence on $k$ may be even less consequential in practice than the error curves suggest; a direct test would be to compare optimization trajectories with the exact adjoint gradient on a moderately sized problem.","The paper's full-structure FWI result failed to converge in ten iterations, attributed to ill-posedness rather than the approximation; a testable extension would be to run the same large problem with the exact adjoint method on a machine with enough memory to isolate the approximation's contribution to the failure."],"forward_implications":["Memory per sensitivity computation drops from $O(\\text{dofs} \\times \\text{time steps})$ to $O(\\text{dofs})$, so the practical limit for dynamic adjoint optimization shifts from memory capacity to wall-clock time.","On a 40 GB GPU the theoretical single-precision limit is about $2.5 \\times 10^9$ parameters, roughly three orders of magnitude beyond what the standard adjoint method reaches with thousands of time steps.","The approximation adds no extra computational work beyond one forward and one backward solve, unlike checkpointing or compression schemes that trade extra compute or accuracy for memory.","A usable $k$ can be calibrated during the first gradient computation by decreasing $k$ until single- and double-precision results diverge, without needing a reference gradient.","The method is limited to self-adjoint, damping-free dynamics and time-reversible time integrators, so dissipative wave problems are outside its scope."],"supporting_citations":[{"why":"Provides the adjoint-state method for computing gradients in geophysical inverse problems, which is the foundation of the sensitivity analysis the paper approximates.","marker":"[24]"},{"why":"Serves as the tutorial reference for the adjoint method and supports the derivation of the adjoint wave equation and Fr\\'echet kernel.","marker":"[27]"},{"why":"Gives the density-scaling full waveform inversion parameterization and the specific Fr\\'echet kernel used in the FWI demonstration.","marker":"[56]"},{"why":"Provides the convolutional neural network finite difference solver and baseline timings used for comparison in the forward and gradient computations.","marker":"[23]"},{"why":"Represents the checkpointing strategy whose memory-versus-compute trade-off the new method aims to avoid.","marker":"[28]"},{"why":"Represents the wavefield compression approach that the new method claims to supersede as a memory-reduction technique.","marker":"[32]"},{"why":"Describes an existing memory-efficient gradient computation by re-simulating the forward wavefield backward in time, an alternative the paper positions its approach against.","marker":"[37, 38]"}],"fun_headline_variants":["Superposition slashes memory for billion-parameter wave inversion","Adjoint memory cut: Superposition trick fits billion-param on one GPU","No history storage: New adjoint method scales to billion DOF on A100","Superposition-based adjoint: memory now just a few time steps"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The backward pass must be able to obtain the adjoint source at every time step, and the paper does not state where that time history lives when the source depends on the forward wave over a large region; if it has to be stored, the memory claim fails.","fun_headline_variants_meta":{"raw":{"variants":["Superposition slashes memory for billion-parameter wave inversion","Adjoint memory cut: Superposition trick fits billion-param on one GPU","No history storage: New adjoint method scales to billion DOF on A100","Superposition-based adjoint: memory now just a few time steps"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000667,"raw_usage":{"total_tokens":3099,"prompt_tokens":1058,"completion_tokens":2041,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":674,"completion_tokens_details":{"reasoning_tokens":1963}},"tokens_in":674,"tokens_out":2041,"duration_ms":11998,"temperature":1.0,"reasoning_tokens":1963,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-15T15:50:40.510185+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run Algorithm 1 with a volume-supported adjoint source that depends on the forward field over the whole domain (a TATO problem with $\\Omega_s = \\Omega$). If the method's peak memory stays at a few solution grids regardless of the number of time steps $N$, the O(dofs) claim holds; if the implementation must store or regenerate the forward field over $\\Omega_s$ for every time step and peak memory grows with $N$, the claim is falsified for that class.","supporting_citations":[{"cited_title":"Achieving logarithmic growth of temporal and spatial complexity in reverse automatic differentiation.Optimization Methods and Software, 1(1):35–54, January 1992","cited_arxiv_id":null,"evidence_quote":"Represents the checkpointing strategy whose memory-versus-compute trade-off the new method aims to avoid."},{"cited_title":"Wavefield compression for 26 adjoint methods in full-waveform inversion.GEOPHYSICS, 81(6):R385–R397, 2016","cited_arxiv_id":null,"evidence_quote":"Represents the wavefield compression approach that the new method claims to supersede as a memory-reduction technique."}],"review_version":2}