REVIEW 3 major objections 5 minor 1 cited by
Polynomial chaos turns operator learning into a linear solve
Reviewed by Pith at T0; open to challenge. T0 means a machine referee read the full paper against a public rubric. the ladder, T0–T4 →
PCE coefficients, solved by least squares or PDE residuals, can serve as an operator learning surrogate that also yields UQ moments.
T0 review reviewed 2026-08-05 challenge →
load-bearing objection PC2 is a clean re-framing of PCE as operator learning, but the physics-informed derivation only works for deterministic operators; the advection–diffusion example uses a random coefficient and the stated equations don't cover it. the 3 major comments →
Polynomial Chaos Expansion for Operator Learning
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
The reading
Core claim
The central claim is that an operator mapping input functions to PDE solutions can be captured by a double spectral expansion where the solution is a sum of products of orthogonal polynomials in the random input parameters and orthogonal polynomials in space-time. This reduces operator learning to recovering a coefficient matrix from the factorization S = Phi C Psi. The paper derives the closed-form minimizer C* = (Phi^T Phi)^{-1} Phi^T S Psi^T (Psi Psi^T)^{-1} for data-driven learning, and for PC2 derives a linear system from the PDE, boundary, and initial residuals, with a Newton iteration for nonlinear problems. It further shows that the first two moments of the predicted solution are obt
What carries the argument
The central object is the factorized bivariate spectral expansion S = Phi C Psi, with Phi a matrix of orthogonal spatio-temporal polynomials and Psi a matrix of orthogonal stochastic polynomials matched to the input distribution. Orthogonality converts the operator-learning objective into least-squares problems for the coefficient matrix C, yielding a unique solution when the number of space-time collocation points is at least the number of physical basis terms and the number of input samples is at least the number of stochastic basis terms. In the physics-informed PC2 variant, the same factorization lets the PDE operator act on the physical basis alone, so the residual is linear in C for li
Load-bearing premise
The physics-informed PC2 derivation assumes that the differential operator can be applied to the spatio-temporal basis alone, with the stochastic polynomials passing through untouched; this holds only when the PDE's random inputs appear as source or initial-condition terms and breaks when a random field multiplies a derivative, as in the advection-diffusion example.
What would settle it
Run PC2 on the paper's advection-diffusion equation with a random velocity field that has nonzero mean, and compare the exact PDE residual applied to a product basis function with the factorized form used in Eq. (15) at interior points. The difference involves the random coefficient multiplying the derivative of the physical basis, so the residual after solving Eq. (20) will not converge to machine precision; measuring that stall gives a concrete test of the derivation.
If this is right
- Operator learning becomes a direct linear algebra procedure: no stochastic gradient descent, no network architecture search, and a unique closed-form solution as soon as sample and collocation counts exceed basis counts.
- Uncertainty quantification comes for free at inference: once C is known, predicted mean and covariance at any query point are matrix products, so ensembles or repeated simulations are unnecessary.
- PC2 removes the need for labeled output data in linear problems; for nonlinear problems, Newton iterations on the polynomial residual provide a data-free surrogate.
- Sparse truncation of the stochastic basis is what makes the 21-dimensional heat example tractable, suggesting the method's practical range is set by basis cardinality rather than by sample availability.
Where Pith is reading between the lines
- Where the input-to-solution map is smooth in the random parameters and the stochastic dimension is small to moderate, this method should beat neural operator learners on cost per accuracy; the curse of dimensionality enters through basis cardinality, so extensions need sparse or adaptive bases.
- The paper's PC2 derivation assumes the differential operator acts only on the physical basis functions, which is violated when a random field multiplies a derivative, as in the advection-diffusion example. A direct extension would expand the product of the random coefficient and the stochastic basis before projecting, turning the linear solve into a larger coupled system.
- Because Sobol sensitivity indices are analytic functions of the coefficient matrix, the same coefficients could yield global sensitivity maps of the learned operator, an output the paper mentions but does not demonstrate.
- A testable extension is to feed PC2 residual norms back into basis selection: adaptive refinement of the physical or stochastic basis based on coefficient magnitudes could automatically handle the steep gradients seen in Burgers and boundary layers.
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
Summary. The manuscript proposes polynomial chaos expansion (PCE) as an operator learning method. The solution operator G:u→s is approximated as Ŝ(x,t,ξ)=Φ(x,t) C Ψ(ξ), with Φ a multivariate polynomial basis in the physical variables, Ψ a PCE basis in the random input variables, and C a matrix of coefficients. In the data-driven setting, C is obtained by linear least squares, giving the closed-form solution in Eq. (14). In the physics-informed setting (PC2), C is obtained by minimizing PDE, boundary, and initial condition residuals; for linear PDEs the manuscript derives the closed-form normal-equation solution in Eq. (20), while nonlinear PDEs are said to be solved by Newton-Raphson iteration as in Eq. (21). Mean and covariance estimates are obtained by post-processing C via Eqs. (22)-(23). Numerical experiments cover an anti-derivative problem, a random-coefficient advection-diffusion equation, Burgers' equation with a random source, and a 2D heat equation with a random source, with comparisons to FEM and Monte Carlo references.
Significance. If the derivations were valid in full generality, the paper would provide a simple, interpretable, and computationally cheap alternative to neural operators for low-dimensional random-parameter problems, with the attractive feature of built-in UQ at no extra cost. The data-driven least-squares formulation is standard and appears correct; the UQ post-processing is elegant; and the anti-derivative and 2D-heat examples are consistent with the stated theory. However, two of the four numerical examples rely on a factorization assumption in Eq. (15) that the manuscript neither states nor justifies, so the central claim that PC2 handles random-coefficient and nonlinear PDEs is not supported by the derivation. The paper is a potentially useful contribution once this is repaired, but in its current form the physics-informed component is not self-consistent.
major comments (3)
- [§2, Eq. (15); §3.2, Eq. (25)] Eq. (15) assumes Nν[Φβ Ψα] = Nν[Φβ] Ψα. This factorization holds for a deterministic linear operator but not for the advection-diffusion problem of Eq. (25), where the velocity v(x) is a GRF realized in ξ. Substituting Eq. (9) into Nν[s]=∂t s+v(x,ξ)∂x s−D∂xx s produces the term v(x,ξ)∂xΦβ(x,t)Ψα(ξ), which cannot be written as a realization-independent function of (x,t) times Ψα(ξ). Consequently, the matrices ΦPDE in Eq. (16) and the normal equations in Eqs. (18)-(20) omit the v-Ψ coupling. The reported PC2 MSE of 4.79e-5 cannot be obtained from the stated minimization unless the implementation uses a different loss or an undocumented treatment of the product term. The authors should derive a coupled system (e.g., expand v in the same PCE basis and account for triple products) or explicitly restrict PC2 to operators with deterministic coefficients.
- [§2, Eqs. (15)-(16); §3.3, Eq. (26)] The same factorization is also invalid for the nonlinear Burgers equation in Eq. (26). For Nν[s]=∂t s+s∂x s−ν∂xx s, replacing s by the expansion in Eq. (9) yields quadratic terms with products of coefficients and products of basis functions ΨαΨγ. The residual expression in Eq. (16) is nevertheless linear in C, so it cannot represent the quadratic term. The text states at Eq. (21) that a nonlinear system is solved, but that system is never written down. Please provide the actual residual, gradient, and Hessian formulas used for the Burgers example, or state the numerical optimization objective and algorithm actually implemented.
- [§2, Eq. (21)] The Newton-Raphson update C(k+1)=C(k)−H_L(C(k))^{-1}∇L(C(k)) is dimensionally inconsistent as written: ∇L is declared to be in R^{Q×P}, so its Hessian is a fourth-order tensor, not a matrix whose inverse can be applied to a Q×P gradient. The text says that 'closed-form expressions' can be derived but provides none. Since the only nonlinear example is Burgers' equation, the reader cannot verify the reported MSE of 4.41e-5. Please provide the vectorized system that is actually solved, or specify the optimizer used.
minor comments (5)
- [§2, Eqs. (22)-(23)] The mean and covariance formulas assume the multivariate polynomials Ψα are orthonormal with E[Ψα^2]=1. The paper only states orthogonality. If a non-normalized Wiener-Askey basis is used, Eq. (23) needs E[Ψα^2] factors in the diagonal matrix.
- [§2, Eq. (17)] The loss in Eq. (17) and the normal equations in Eq. (20) are written for an empirical sample average over N realizations and n virtual points. This should be stated explicitly, since it is not a Galerkin projection in the continuous L2 sense and the distinction matters for the interpretation of the PC2 solution.
- [§3.4, Table 1] The reported MAE of 6.48e-16 for the mean in the 2D heat problem is suspiciously close to machine precision. Please explain whether this is due to the source field having zero mean, or whether it is a statistical artifact of the comparison.
- [General] The paper does not include a code/data availability statement. Given that the PC2 derivation is incomplete for two of the four examples, releasing the implementation would be important for reproducibility.
- [Keywords/Abstract] There is a typo in the keywords: 'Operator leaning' should be 'Operator learning'.
Circularity Check
No significant circularity: the PCE and PC2 derivations are self-contained algebraic constructions, and the self-citations are not load-bearing.
full rationale
The paper's derivation chain is self-contained. The data-driven PCE solution in Eq. (14) follows directly from the least-squares objective in Eq. (13) via the normal equations; no fitted parameter is later relabeled as a prediction. The UQ formulas in Eqs. (22)–(23) are algebraic consequences of the PCE orthogonality and the matrix form Eq. (10), and they are validated against independent Monte Carlo references rather than fitted to them. The physics-informed PC2 derivation in Eqs. (15)–(20) is a direct linear-algebra consequence of minimizing the residual loss in Eq. (17); it does not assume the numerical results it later reports. The paper cites its own prior PC2/m-PCE work (refs. 24–28), but explicitly states that the present formulation differs and is more efficient, so those citations are not load-bearing. The skeptical concern about random-coefficient operators and the factorization in Eq. (15) is a mathematical validity/approximation question, not a circularity: the algebra does not reduce to its inputs by construction. No circular step satisfying the required standard (Eq. X = Eq. Y by construction, or fitted parameter renamed as prediction) is present.
Axiom & Free-Parameter Ledger
free parameters (4)
- Stochastic polynomial order p =
3 for first three examples; 4 for 2D heat
- Spatio-temporal polynomial order q =
10, 14, 23, 16 for the four examples
- Number of KL terms r =
6 for first three examples; 21 for 2D heat (20 source terms plus 1 amplitude)
- Hyperbolic truncation factor =
0.9 for 2D heat
axioms (5)
- domain assumption Input functions are Gaussian random fields representable by a truncated Karhunen-Loeve expansion with independent standard normal variables xi.
- domain assumption The solution operator is smooth enough in xi and (x,t) that low-order polynomial truncation (p,q) is accurate.
- ad hoc to paper In PC2, the differential operator acts only on the spatio-temporal basis functions, i.e., N[Phi_beta Psi_alpha] = N[Phi_beta] Psi_alpha.
- domain assumption Newton-Raphson iteration for nonlinear PC2 converges, and the gradient and Hessian have closed-form expressions.
- domain assumption Monte Carlo simulation with 10,000 samples provides converged reference statistics.
Cite this review
Pith. "Pith review of Polynomial Chaos Expansion for Operator Learning." pith.science (2026). https://pith.science/paper/HWHJJNHS
@misc{pith2026250820886,
author = {Pith},
title = {Pith review of: Polynomial Chaos Expansion for Operator Learning},
year = {2026},
howpublished = {\url{https://pith.science/paper/HWHJJNHS}},
note = {Machine review of arXiv:2508.20886}
}
read the original abstract
Operator learning (OL) has emerged as a powerful tool in scientific machine learning (SciML) for approximating mappings between infinite-dimensional functional spaces. One of its main applications is learning the solution operator of partial differential equations (PDEs). While much of the progress in this area has been driven by deep neural network-based approaches such as Deep Operator Networks (DeepONet) and Fourier Neural Operator (FNO), recent work has begun to explore traditional machine learning methods for OL. In this work, we introduce polynomial chaos expansion (PCE) as an OL method. PCE has been widely used for uncertainty quantification (UQ) and has recently gained attention in the context of SciML. For OL, we establish a mathematical framework that enables PCE to approximate operators in both purely data-driven and physics-informed settings. The proposed framework reduces the task of learning the operator to solving a system of equations for the PCE coefficients. Moreover, the framework provides UQ by simply post-processing the PCE coefficients, without any additional computational cost. We apply the proposed method to a diverse set of PDE problems to demonstrate its capabilities. Numerical results demonstrate the strong performance of the proposed method in both OL and UQ tasks, achieving excellent numerical accuracy and computational efficiency.
Figures
Forward citations
Cited by 1 Pith paper
-
Probabilistic operator learning: generative modeling and uncertainty quantification for foundation models of differential equations
ICON is shown to compute the posterior predictive mean of differential equation solutions, and a generative extension, GenICON, provides samples from this distribution for uncertainty quantification.
Reference graph
Works this paper leans on
-
[1]
Learning nonlinear operators via deeponet based on the universal approximation theorem of operators,
L. Lu, P. Jin, G. Pang, Z. Zhang, and G. E. Karniadakis, “Learning nonlinear operators via deeponet based on the universal approximation theorem of operators,” Nature machine intelligence, vol. 3, no. 3, pp. 218–229, 2021
work page 2021
-
[2]
Fourier neural operator for parametric partial differential equations,
Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anand- kumar, “Fourier neural operator for parametric partial differential equations,”arXiv preprint arXiv:2010.08895, 2020
Pith/arXiv arXiv 2010
-
[3]
T. Tripura and S. Chakraborty, “Wavelet neural operator for solving parametric partial dif- ferential equations in computational mechanics problems,” Computer Methods in Applied Mechanics and Engineering, vol. 404, p. 115783, 2023
work page 2023
-
[4]
Laplace neural operator for solving differential equations,
Q. Cao, S. Goswami, and G. E. Karniadakis, “Laplace neural operator for solving differential equations,” Nature Machine Intelligence, vol. 6, no. 6, pp. 631–640, 2024
work page 2024
-
[5]
Convolutional neural operators,
B. Raonic, R. Molinaro, T. Rohner, S. Mishra, and E. de Bezenac, “Convolutional neural operators,” in ICLR 2023 workshop on physics for machine learning, 2023
work page 2023
-
[6]
Neural operator: Graph kernel network for partial differential equations,
Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandku- mar, “Neural operator: Graph kernel network for partial differential equations,” arXiv preprint arXiv:2003.03485, 2020
Pith/arXiv arXiv 2003
-
[7]
Kernel methods are competitive for operator learning,
P. Batlle, M. Darcy, B. Hosseini, and H. Owhadi, “Kernel methods are competitive for operator learning,” Journal of Computational Physics, vol. 496, p. 112549, 2024
work page 2024
-
[8]
Operator learning with gaussian processes,
C. Mora, A. Yousefpour, S. Hosseinmardi, H. Owhadi, and R. Bostanabad, “Operator learning with gaussian processes,” Computer Methods in Applied Mechanics and Engineering, vol. 434, p. 117581, 2025
work page 2025
-
[9]
Operator learning using random features: A tool for scientific computing,
N. H. Nelsen and A. M. Stuart, “Operator learning using random features: A tool for scientific computing,” SIAM Review, vol. 66, no. 3, pp. 535–571, 2024
work page 2024
-
[10]
M. Lowery, J. Turnage, Z. Morrow, J. D. Jakeman, A. Narayan, S. Zhe, and V . Shankar, “Kernel neural operators (knos) for scalable, memory-efficient, geometrically-flexible operator learning,” arXiv preprint arXiv:2407.00809, 2024
Pith/arXiv arXiv 2024
-
[11]
The Wiener–Askey polynomial chaos for stochastic differential equations,
D. Xiu and G. E. Karniadakis, “The Wiener–Askey polynomial chaos for stochastic differential equations,” SIAM Journal on Scientific Computing, vol. 24, no. 2, pp. 619–644, 2002
work page 2002
-
[12]
C. Soize and R. Ghanem, “Physical systems with random uncertainties: chaos representations with arbitrary probability measure,” SIAM Journal on Scientific Computing, vol. 26, no. 2, pp. 395–410, 2004
work page 2004
-
[13]
Neural chaos: A spectral stochastic neural operator,
B. Bahmani, I. G. Kevrekidis, and M. D. Shields, “Neural chaos: A spectral stochastic neural operator,” arXiv e-prints, pp. arXiv–2502, 2025
work page 2025
-
[14]
R. G. Ghanem and P. D. Spanos, Stochastic finite elements: a spectral approach . Courier Corporation, 2003
work page 2003
-
[15]
Stochastic collocation methods: a survey,
D. Xiu, “Stochastic collocation methods: a survey,” Handbook of Uncertainty Quantification, pp. 699–716, 2016
work page 2016
-
[16]
High-order collocation methods for differential equations with random inputs,
D. Xiu and J. S. Hesthaven, “High-order collocation methods for differential equations with random inputs,” SIAM Journal on Scientific Computing, vol. 27, no. 3, pp. 1118–1139, 2005. 17
work page 2005
-
[17]
Stochastic collocation methods on unstructured grids in high dimen- sions via interpolation,
A. Narayan and D. Xiu, “Stochastic collocation methods on unstructured grids in high dimen- sions via interpolation,” SIAM Journal on Scientific Computing, vol. 34, no. 3, pp. A1729– A1752, 2012
work page 2012
-
[18]
Stochastic finite element: a non intrusive ap- proach by regression,
M. Berveiller, B. Sudret, and M. Lemaire, “Stochastic finite element: a non intrusive ap- proach by regression,” Revue Européenne de Mécanique Numérique/European Journal of Computational Mechanics, vol. 15, no. 1-2-3, pp. 81–92, 2006
work page 2006
-
[19]
Efficient collocational approach for parametric uncertainty analysis,
D. Xiu, “Efficient collocational approach for parametric uncertainty analysis,”Communications in Computational Physics, vol. 2, no. 2, pp. 293–309, 2007
work page 2007
-
[20]
Global sensitivity analysis using polynomial chaos expansions,
B. Sudret, “Global sensitivity analysis using polynomial chaos expansions,” Reliability Engi- neering & System Safety, vol. 93, no. 7, pp. 964–979, 2008
work page 2008
-
[21]
On distribution-based global sensitivity analysis by polynomial chaos expansion,
L. Novák, “On distribution-based global sensitivity analysis by polynomial chaos expansion,” Computers & Structures, vol. 267, p. 106808, 2022
work page 2022
-
[22]
On fractional moment estimation from polynomial chaos expansion,
L. Novák, M. Valdebenito, and M. Faes, “On fractional moment estimation from polynomial chaos expansion,” Reliability Engineering & System Safety, vol. 254, p. 110594, 2025
work page 2025
-
[23]
Data-driven polynomial chaos expansion for machine learning regression,
E. Torre, S. Marelli, P. Embrechts, and B. Sudret, “Data-driven polynomial chaos expansion for machine learning regression,” Journal of Computational Physics, vol. 388, pp. 601–623, 2019
work page 2019
-
[24]
Physics-informed polynomial chaos expansions,
L. Novák, H. Sharma, and M. D. Shields, “Physics-informed polynomial chaos expansions,” Journal of Computational Physics, vol. 506, p. 112926, 2024
work page 2024
-
[25]
H. Sharma, L. Novák, and M. Shields, “Physics-constrained polynomial chaos expansion for scientific machine learning and uncertainty quantification,” Computer Methods in Applied Mechanics and Engineering, vol. 431, p. 117314, 2024
work page 2024
-
[26]
Man- ifold learning-based polynomial chaos expansions for high-dimensional surrogate models,
K. Kontolati, D. Loukrezis, K. R. Dos Santos, D. G. Giovanis, and M. D. Shields, “Man- ifold learning-based polynomial chaos expansions for high-dimensional surrogate models,” International Journal for Uncertainty Quantification, vol. 12, no. 4, 2022
work page 2022
-
[27]
K. Kontolati, D. Loukrezis, D. G. Giovanis, L. Vandanapu, and M. D. Shields, “A survey of unsupervised learning methods for high-dimensional uncertainty quantification in black-box- type problems,” Journal of Computational Physics, vol. 464, p. 111313, 2022
work page 2022
-
[28]
On the influence of over-parameterization in manifold based surrogates and deep neural operators,
K. Kontolati, S. Goswami, M. D. Shields, and G. E. Karniadakis, “On the influence of over-parameterization in manifold based surrogates and deep neural operators,” Journal of Computational Physics, vol. 479, p. 112008, 2023
work page 2023
-
[29]
Adaptive sparse polynomial chaos expansion based on least angle regression,
G. Blatman and B. Sudret, “Adaptive sparse polynomial chaos expansion based on least angle regression,” Journal of Computational Physics, vol. 230, no. 6, pp. 2345–2367, 2011
work page 2011
-
[30]
Adam: A method for stochastic optimization,
D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” arXiv preprint arXiv:1412.6980, 2014
Pith/arXiv arXiv 2014
-
[31]
On sensitivity estimation for nonlinear mathematical models,
I. M. Sobol’, “On sensitivity estimation for nonlinear mathematical models,”Matematicheskoe Modelirovanie, vol. 2, no. 1, pp. 112–118, 1990
work page 1990
-
[32]
Uqpy v4. 1: Uncertainty quantification with python,
D. Tsapetis, M. D. Shields, D. G. Giovanis, A. Olivier, L. Novak, P. Chakroborty, H. Sharma, M. Chauhan, K. Kontolati, L. Vandanapu,et al., “Uqpy v4. 1: Uncertainty quantification with python,” SoftwareX, vol. 24, p. 101561, 2023. 18
work page 2023
-
[33]
Jax: composable transformations of python+ numpy programs,
J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Milne,et al., “Jax: composable transformations of python+ numpy programs,” 2018
work page 2018
-
[34]
The FEniCS project version 1.5,
M. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes, and G. N. Wells, “The FEniCS project version 1.5,”Archive of Numerical Software, vol. 3, no. 100, 2015
2015
-
[35]
Numerical Approximation in CFD Problems Using Physics Informed Machine Learning
S. Rout, V . Dwivedi, and B. Srinivasan, “Numerical approximation in cfd problems using physics informed machine learning,” arXiv preprint arXiv:2111.02987, 2021. 19
work page internal anchor Pith review Pith/arXiv arXiv 2021
This paper was first reviewed by deepseek-v4-flash on August 5, 2026.
discussion (0)
Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.