REVIEW 2 major objections 5 minor 1 cited by
Hybrid Least Squares/Gradient Descent Methods for DeepONets
T0 review · 2 major / 5 minor · reviewed 2026-08-05 · deepseek-v4-flash
Pith's one-line read A hybrid least-squares/gradient-descent method cuts DeepONet training compute by a factor of ten.
desk verdict Solid math, honest limitations, but the headline acceleration claim needs wall-clock timings before I'd trust the 10x speedup. read the letter →
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
The reading
What carries the argument
The factorization identity A_k = K_{P,Q_k}(T_k ⊗ B) (Theorem 3), where B is the branch pre-output matrix, T_k is the trunk output matrix under the loss operator L_k, and K_{P,Q_k} is a commutation matrix rearranging row order. This identity reduces the large least-squares problem eq. (3.11) to the generalized Sylvester equation (3.12), whose solution is given in closed form (Proposition 6) via spectral decompositions of the two Gram matrices, making the 'LS step' of hybrid training a small dense linear algebra solve.
What would settle it
Profile the per-work-unit wall-clock time of LS+Adam against Adam-only on the advection benchmark: if a single LS step (one spectral decomposition of B^T B and Σ ε_k T_k^T T_k plus the entrywise solve) costs more than about 5% of the five Adam epochs it replaces, the claimed speedup at equal error will not materialize in wall-clock time even though the work-unit count to a fixed error drops. Also, trying to apply the closed-form solve to a physics loss whose operator depends on the input function (the variable-coefficient Poisson residual) should fail condition (3.5), so the method's reported
Extended reading notes
Core claim
The central claim: for vanilla DeepONets, the branch network's last-layer weights can be optimized by an exact least-squares solve, not backpropagation, and this is feasible because the full LS matrix factors as A_k = K_{P,Q_k}(T_k ⊗ B) when the data is a Cartesian product and each loss operator acts only on trunk outputs. The normal equations then reduce to the matrix equation B^T B C^T (Σ ε_k T_k^T T_k) + λ C^T = B^T (Σ ε_k F_k^T T_k), a generalized Sylvester equation solved in closed form by spectral decomposition. The paper further claims LS+Adam—Adam on hidden layers plus this LS solve—matches or beats Adam-only's validation L2 error at 10,000 work units versus 100,000, across advection
Load-bearing premise
Everything rests on the loss data being a Cartesian product of input functions and evaluation points with each linear loss operator acting only on trunk outputs; if either fails, the closed-form least-squares update cannot be formed, and the claimed speedup also assumes the LS step's wall-clock overhead stays at only a few percent per work unit.
Editorial extensions
If this is right
- On the tested PDE benchmarks, LS+Adam reaches the same or lower validation L2 error at 10,000 work units that Adam-only reaches only at 100,000 work units.
- The hybrid method extends to physics-informed losses for linear PDEs, where the residual and initial/boundary operators satisfy the factorization conditions for several tested problems.
- The regularization term λ‖C‖² makes the LS system full-rank and can be annealed over training, for both supervised and unsupervised settings.
- The factorization splits the LS problem into branch and trunk subproblems, so the cost of an LS step is governed by the small matrices B and T_k rather than the full P·Q_k × I·J system.
- The scheme is limited to losses of the sum-of-squared-l2 form with linear operators acting on trunk outputs; the nonlinear reaction term and the variable-coefficient Poisson physics loss are explicitly outside its scope.
Reading between the lines
- This factorization is effectively an exact coordinate-descent step on the DeepONet loss; one could analyze hybrid training as alternating minimization and potentially prove convergence-rate statements under strong convexity of the last-layer subproblem — an analysis the paper does not attempt.
- The same Kronecker-product trick should transfer to MIONet-style multi-input operators, where the tensor-product structure of the branch outputs would likely yield an analogous commutation-matrix factorization.
- The advertised speedup is stated in work units; converting it into wall-clock speedup depends on the LS solve's overhead, so a lightweight implementation that reuses spectral decompositions across LS steps could make the tenfold claim hold in real time as well.
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
Summary. The paper proposes a hybrid least-squares/gradient-descent (LSGD) training scheme for DeepONets. Because the DeepONet output is linear in the last-layer matrix C of the branch network, the authors fix the hidden parameters and update C by solving least-squares normal equations. Their main theoretical contribution (Theorem 3) is a factorization of the large LS matrix as A_k = K_{P,Q_k}(T_k ⊗ B) under a Cartesian-product data assumption and a linear-operator locality assumption (3.4)-(3.5). This reduces the normal equations to a Sylvester-type matrix equation whose solution is given in closed form by Proposition 6 via spectral decompositions. Algorithm 2 (LS+Adam) alternates blocks of Adam steps on hidden parameters with full-batch LS solves of the last layer. Experiments on advection, diffusion-reaction, and Poisson problems with supervised and physics-informed losses report that LS+Adam reaches at 10,000 work units the validation L2 error that Adam-only reaches at 100,000 work units.
Significance. The algebraic core is clean and appears correct; the factorization and closed-form solve are nontrivial and could be useful beyond this paper. The scope conditions (3.4)-(3.5) are stated explicitly, and the limitations for nonlinear operators and input-dependent linear operators are acknowledged in Sections 4.2 and 4.3.1. The paper honestly separates the exact LS step from the heuristic alternation with Adam. The main weakness is that the headline 'acceleration' claim rests on an unmeasured conversion from work units to wall-clock time. If the per-WU overhead is as stated, the 10x WU improvement is a substantial practical advance; without timing data, the paper demonstrates only faster convergence per WU.
major comments (2)
- [§4, Work Unit definition and Algorithm 2] The sentence 'the training time for a single WU in the hybrid stage is 2–5% longer than that of Adam's' is the only quantitative link between work units and wall-clock time, but no timing measurement is reported. One hybrid WU performs a full-batch LS solve (forward pass over all P=1000 branch inputs, formation of B^T B and T_k^T T_k, and two eigendecompositions), while an Adam-only WU uses batches of 50 functions; the 2–5% overhead is therefore not self-evident. Because the abstract promises accelerated DeepONet training, the experiments should report wall-clock runtimes (or at least measured overhead with variance) for the main benchmarks. Otherwise the 10x WU gap demonstrates faster convergence per WU, not necessarily per second.
- [§4, Figures 2, 5, 7] The convergence curves and the comparisons 'LS+Adam at 10,000 WU vs Adam-only at 100,000 WU' appear to come from single runs, with no number of random seeds, confidence intervals, or variance. Given the stochasticity of Adam and initialization, a few independent runs (e.g., 5) with mean and standard deviation of validation L2 error at the comparison checkpoints would make the empirical claim robust. This is secondary to the wall-clock issue but is needed for a quantitative evaluation.
minor comments (5)
- [§3, Proposition 6] The proposition states λ≥0, but the entrywise division in (3.14) requires dA,r dB,s + λ > 0; with λ=0 and singular A or B the formula is undefined. The text later sets λ>0, so the proposition should state that condition.
- [§3, after Proposition 6] The sentence 'we can find the last layer parameter C ... that minimizes eq. (3.6)' should likely refer to the LS problem (3.11), not the factorization (3.6).
- [References] Reference [17] is incomplete: the Adam reference should include J. Ba as a coauthor.
- [Abstract] The phrase 'two smaller subproblems — one for the branch network and one for the trunk network' is a bit stronger than what is shown: the method solves a structured matrix equation using spectral decompositions of both Gram matrices, rather than literally decoupling the two networks.
- [§4, Algorithm 2 and its description] The paper says Adam momentums are maintained across the LS step; since θL is overwritten by the exact least-squares solution, the stale Adam moments for θL may be inconsistent with the new parameter value. A sentence explaining why this does not destabilize training would help.
Circularity Check
No circularity: the LS factorization and closed-form solve are derived from stated assumptions with standard linear algebra; the acceleration evidence is empirical, not definitional.
full rationale
The paper's derivation chain is self-contained. The LS problem (3.3) is formed directly from the loss (3.1) and the DeepONet representation (2.1). The factorization theorem (Theorem 3, Eq. 3.6) is proved from the explicit assumptions (3.4) and (3.5) using the definitions of the commutation matrix and Kronecker product; it is not assumed or fitted. The normal equation (3.12) follows from Lemmas 4 and 5, and Proposition 6 gives a closed-form solution by spectral decomposition. No fitted constant is renamed as a prediction. The reported WU-based speedups are empirical comparisons; the assertion that a hybrid WU costs only 2–5% more wall-clock time is not backed by timing data, but that is an evidentiary gap about cost, not a circular step. The load-bearing citations are to external prior work (Cyr et al., Chen & Chen, Lu et al., Magnus & Neudecker) and standard linear algebra; the only self-citations involving author Chang-Ock Lee ([18,19]) appear in the introduction as examples of acceleration methods and are not used to justify the paper's core claims. Therefore no step reduces to its own inputs, and there is no significant circularity.
Assumptions & free parameters
free parameters (5)
- λ (last-layer regularization weight) =
10^-6, 10^-9, or 10^-14 per experiment; decays from 10^-9 to 10^-14 for the Poisson-source case
- ϵ2 (physics-loss weight) =
10^-1, 10^-4, or 10^-14 depending on problem
- Number of Adam epochs per LS step (J0) =
5
- Initial Adam stage length (I0) =
500 epochs
- λ decay schedule =
10^-9 to 10^-14 between 100 and 1000 WU for Poisson-source case
assumptions (5)
- domain assumption Data for each loss term has the Cartesian product structure χ_k = β × τ_k (eq. 3.4).
- domain assumption Each linear loss operator L_k acts only on trunk outputs: L_k[\tilde b_j t_i](u,y) = \tilde b_j(u) L_k[t_i](y) (eq. 3.5).
- domain assumption DeepONet with a bias-free, activation-free linear last layer in the branch possesses the universal approximation property (Theorem 2 and Section 2.1).
- standard math Commutation and Kronecker product identities from Magnus and Neudecker [28] hold.
- standard math Symmetric positive semi-definite matrices have spectral decompositions.
Cite this review
Pith. "Pith review of Hybrid Least Squares/Gradient Descent Methods for DeepONets." pith.science (2026). https://pith.science/paper/GXF3TFFM
@misc{pith2026250815394,
author = {Pith},
title = {Pith review of: Hybrid Least Squares/Gradient Descent Methods for DeepONets},
year = {2026},
howpublished = {\url{https://pith.science/paper/GXF3TFFM}},
note = {Machine review of arXiv:2508.15394}
}
abstract
We propose an efficient hybrid least squares/gradient descent method to accelerate DeepONet training. Since the output of DeepONet can be viewed as linear with respect to the last layer parameters of the branch network, these parameters can be optimized using a least squares (LS) solve, and the remaining hidden layer parameters are updated by means of gradient descent form. However, building the LS system for all possible combinations of branch and trunk inputs yields a prohibitively large linear problem that is infeasible to solve directly. To address this issue, our method decomposes the large LS system into two smaller, more manageable subproblems $\unicode{x2014}$ one for the branch network and one for the trunk network $\unicode{x2014}$ and solves them separately. This method is generalized to a broader type of $L^2$ loss with a regularization term for the last layer parameters, including the case of unsupervised learning with physics-informed loss.
Figures
Figures from the paper (6 more)
Forward citations
Cited by 1 Pith paper
-
Multilateralism in the Global Governance of Artificial Intelligence
AI multilateralism looks multi-stakeholder but operates in the shadow of state hierarchy, framed by principles of epochal change, determinism, and dialectical understanding.
Reference graph
Works this paper leans on
-
[1]
M. T. Augustine , A survey on universal approximation theorems , arXiv preprint arXiv:2407.12895, (2024)
arXiv 2024
-
[2]
A. R. Barron, Universal approximation bounds for superpositions of a sigmoidal function , IEEE Transactions on Information Theory, 39 (1993), pp. 930–945
work page 1993
-
[3]
R. H. Bartels and G. W. Stewart, Algorithm 432 [C2]: solution of the matrix equation AX+XB=C [F4] , Communications of the ACM, 15 (1972), pp. 820–826
work page 1972
-
[4]
A. G. Baydin, B. A. Pearlmutter, A. A. Radul, and J. M. Siskind, Automatic differentiation in machine learning: a survey , Journal of Machine Learning Research, 18 (2018), pp. 1–43
work page 2018
-
[5]
J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. V anderPlas, S. W anderman-Milne, and Q. Zhang, JAX: composable transformations of Python+NumPy programs , 2018
work page 2018
-
[6]
T. Chen and H. Chen, Universal approximation to nonlinear operators by neural net- works with arbitrary activation functions and its application to dynamical systems , IEEE Transactions on Neural Networks, 6 (1995), pp. 911–917
work page 1995
-
[7]
K.-w. E. Chu, The solution of the matrix equations AXB-CXD=E and (YA-DZ, YC- BZ)=(E, F) , Linear Algebra and its Applications, 93 (1987), pp. 93–105
work page 1987
-
[8]
G. Cybenko, Approximation by superpositions of a sigmoidal function , Mathematics of Control, Signals and Systems, 2 (1989), pp. 303–314. 23
work page 1989
Show all 37 references
-
[9]
E. C. Cyr, M. A. Gulian, R. G. Patel, M. Perego, and N. A. Trask, Robust training and initialization of deep neural networks: An adaptive basis viewpoint , in Mathematical and Scientific Machine Learning, PMLR, 2020, pp. 512–536
2020
-
[10]
Dettmers and L
T. Dettmers and L. Zettlemoyer, Sparse networks from scratch: Faster training with- out losing performance , arXiv preprint arXiv:1907.04840, (2019)
1907 arXiv
-
[11]
G. H. Golub, P. C. Hansen, and D. P. O’Leary, Tikhonov regularization and total least squares, SIAM Journal on Matrix Analysis and Applications, 21 (1999), pp. 185–194
1999
-
[12]
K. He, X. Zhang, S. Ren, and J. Sun, Delving deep into rectifiers: Surpassing human- level performance on imagenet classification , in Proceedings of the IEEE International Con- ference on Computer Vision, 2015, pp. 1026–1034
2015
-
[13]
Heinecke, J
A. Heinecke, J. Ho, and W.-L. Hwang, Refinement and universal approximation via sparsely connected ReLU convolution nets , IEEE Signal Processing Letters, 27 (2020), pp. 1175–1179
2020
-
[14]
Hornik, M
K. Hornik, M. Stinchcombe, and H. White, Multilayer feedforward networks are uni- versal approximators, Neural Networks, 2 (1989), pp. 359–366
1989
-
[15]
D.-K. Jang, K. Kim, and H. H. Kim, Partitioned neural network approximation for partial differential equations enhanced with Lagrange multipliers and localized loss functions, Computer Methods in Applied Mechanics and Engineering, 429 (2024), p. 117168
2024
-
[16]
P. Jin, S. Meng, and L. Lu, MIONET: Learning multiple-input operators via tensor product, SIAM Journal on Scientific Computing, 44 (2022), pp. A3490–A3514
2022
-
[17]
D. P. Kingma , Adam: A method for stochastic optimization , arXiv preprint arXiv:1412.6980, (2014)
2014 arXiv
-
[18]
Y. Lee, J. Park, and C.-O. Lee, Two-level group convolution , Neural Networks, 154 (2022), pp. 323–332
2022
-
[19]
, IEEE Transactions on Neural Networks and Learning Systems, 35 (2024), pp
, Parareal neural networks emulating a parallel-in-time algorithm. , IEEE Transactions on Neural Networks and Learning Systems, 35 (2024), pp. 6353–6364
2024
-
[20]
Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar, Fourier neural operator for parametric partial differential equations, arXiv preprint arXiv:2010.08895, (2020)
2010 arXiv
-
[21]
, Neural operator: Graph kernel network for partial differential equations , arXiv preprint arXiv:2003.03485, (2020)
2003 arXiv
-
[22]
Lin and S
H. Lin and S. Jegelka, Resnet with one-neuron hidden layers is a universal approximator, Advances in Neural Information Processing Systems, 31 (2018)
2018
-
[23]
D. C. Liu and J. Nocedal, On the limited memory BFGS method for large scale opti- mization, Mathematical Programming, 45 (1989), pp. 503–528
1989
-
[24]
L. Lu, P. Jin, G. Pang, Z. Zhang, and G. E. Karniadakis, Learning nonlinear op- erators via DeepONet based on the universal approximation theorem of operators , Nature Machine Intelligence, 3 (2021), pp. 218–229. 24
2021
-
[25]
L. Lu, X. Meng, S. Cai, Z. Mao, S. Goswami, Z. Zhang, and G. E. Karniadakis, A comprehensive and fair comparison of two neural operators (with practical extensions) based on fair data , Computer Methods in Applied Mechanics and Engineering, 393 (2022), p. 114778
2022
-
[26]
Z. Lu, H. Pu, F. W ang, Z. Hu, and L. W ang, The expressive power of neural networks: A view from the width , Advances in Neural Information Processing Systems, 30 (2017)
2017
-
[27]
Z. Lu, Y. Zhou, Y. Zhang, X. Hu, Q. Zhao, and X. Hu, A fast general thermal simu- lation model based on multi-branch physics-informed deep operator neural network , Physics of Fluids, 36 (2024), p. 037142
2024
-
[28]
J. R. Magnus and H. Neudecker, The Commutation Matrix: Some Properties and Applications, The Annals of Statistics, 7 (1979), pp. 381 – 394
1979
-
[29]
Neudecker, The Kronecker matrix product and some of its applications in econometrics, Statistica Neerlandica, 22 (1968), pp
H. Neudecker, The Kronecker matrix product and some of its applications in econometrics, Statistica Neerlandica, 22 (1968), pp. 69–82
1968
-
[30]
Raissi, P
M. Raissi, P. Perdikaris, and G. E. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, Journal of Computational Physics, 378 (2019), pp. 686–707
2019
-
[31]
Ramachandran, B
P. Ramachandran, B. Zoph, and Q. V. Le, Searching for activation functions , arXiv preprint arXiv:1710.05941, (2017)
2017 arXiv
-
[32]
Salimans and D
T. Salimans and D. P. Kingma, Weight normalization: A simple reparameterization to accelerate training of deep neural networks , Advances in Neural Information Processing Systems, 29 (2016)
2016
-
[33]
Son, ELM-DeepONets: Backpropagation-free training of deep operator networks via extreme learning machines , IEEE Access, 13 (2025), pp
H. Son, ELM-DeepONets: Backpropagation-free training of deep operator networks via extreme learning machines , IEEE Access, 13 (2025), pp. 86927–86934
2025
-
[34]
A. N. Tikhonov, On the solution of ill-posed problems and the method of regularization , in Doklady Akademii Nauk, vol. 151, Russian Academy of Sciences, 1963, pp. 501–504
1963
-
[35]
W ang, H
S. W ang, H. W ang, and P. Perdikaris, Learning the solution operator of parametric partial differential equations with physics-informed DeepONets , Science Advances, 7 (2021), p. eabi8605
2021
-
[36]
Zheng, A
S. Zheng, A. Vishnu, and C. Ding, Accelerating deep learning with shrinkage and recall, in 2016 IEEE 22nd International Conference on Parallel and Distributed Systems (ICPADS), IEEE, 2016, pp. 963–970
2016
-
[37]
Zhou, Universality of deep convolutional neural networks , Applied and Computa- tional Harmonic Analysis, 48 (2020), pp
D.-X. Zhou, Universality of deep convolutional neural networks , Applied and Computa- tional Harmonic Analysis, 48 (2020), pp. 787–794. 25
2020
Reviewed August 5, 2026 · model on record in the stance chip above.
Discussion (0). Sign in to comment.