Pith. sign in

REVIEW 3 major objections 6 minor 6 references

Physics-Informed Linear Model (PILM): Analytical Representations and Application to Crustal Strain Rate Estimation

T0 review · 3 major / 6 minor · reviewed 2026-08-06 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗ Review JSON

Pith's one-line read This paper establishes a physics-informed linear model (PILM) in which linear PDE problems with constant coefficients admit closed-form least-squares solutions, and reports that on Japanese GNSS data, mathematical smoothing is…

desk verdict The analytical PILM framework is solid and reproducible, but the headline claim that mathematical regularization beats physical regularization is undermined by an unnormalized improper-prior comparison. read the letter →

arxiv 2507.12218 v2 pith:DUXY6MO3 submitted 2025-07-16 cs.LG physics.geo-ph

Tomohisa Okazaki This is my paper · ORCID
classification cs.LGphysics.geo-ph MSC 65D0765K1062F1586A22
keywords physics-informedmachinelearningbasisfunctionexpansionB-splinesclosed-formleastsquarescrustalstrainrateGNSSBayesianregularizationpartialdifferentialequations
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

The paper sets out to show that a linear, basis-function version of physics-informed machine learning can make the usual PDE residual loss exactly solvable: when the solution is a linear combination of fixed basis functions whose derivative products can be integrated analytically, the entire loss is a quadratic form and the optimal coefficients are given by a closed-form least-squares formula. This matters because it removes the two numerical error sources, collocation-point quadrature and stochastic optimization, that plague neural-network formulations, leaving only the choice of basis as the modeling decision. The method is demonstrated on damped oscillation, the diffusion equation with unknown coefficient and unknown initial/boundary conditions, and then on crustal strain-rate estimation from 458 GNSS stations in central Japan. In the application, mathematical regularization that penalizes second derivatives of the velocity field is compared, in a Bayesian marginal-likelihood framework, with physical regularization enforcing thin-elastic-sheet force balance; the log marginal likelihood prefers mathematical smoothing by roughly 105 to 113 nats for every Poisson ratio tested, and hybrid regularization does not beat it. If correct, PILM offers an exact, reproducible reference class for linear inverse problems and a statistical argument that elasticity-based force-balance priors are not the right regularizer for these GNSS data.

What carries the argument

The key object is the definite-integral matrix of B-spline derivatives, $\mathbf{R}_{ab} = \int \boldsymbol{\Phi}^{(a)} \boldsymbol{\Phi}^{(b)\mathrm{T}} dx$, which carries the argument. Because cubic B-splines are piecewise cubic, these matrices for derivative orders 0, 1, and 2 are known rational numbers; they are tabulated in the supplement and reused for every problem. The PDE operator is encoded linearly: for the diffusion equation, $\mathbf{G} = \mathbf{R}_{11}\otimes\mathbf{R}_{00} + k^2 \mathbf{R}_{00}\otimes\mathbf{R}_{22} - k(\mathbf{R}_{10}\otimes\mathbf{R}_{02} + \mathbf{R}_{01}\otimes\mathbf{R}_{20})$; for the elastic-sheet equations, the same blocks form $\mathbf{G}_{\mathrm{phys}}$ through similar Kronecker combinations. All that changes between forward, inverse, and geodetic problems is the observation matrix $\mathbf{H}$; the optimizer remains the closed-form least-squares expression. This object matters because it is what makes the method analytical and lets regularization terms be written as prior covariance matrices in a Bayesian model.

What would settle it

Re-run the marginal-likelihood comparison on synthetic velocity fields generated by a known thin-elastic-sheet model with noise; if mathematical smoothing still achieves a higher marginal likelihood than physical regularization on data that truly came from the elastic model, the model-selection criterion is misleading. Alternatively, recompute Table 2 with station-specific or full GNSS covariance in place of the common-$\sigma$ likelihood; if the gap drops below the noise level or reverses, the paper's applied conclusion fails.

Watch

Extended reading notes

Core claim

The central discovery is that for linear PDEs with constant coefficients, the whole physics-informed loss reduces to a matrix quadratic form. With cubic B-splines as basis functions, the definite-integral matrices $\mathbf{R}_{ab} = \int \boldsymbol{\Phi}^{(a)} \boldsymbol{\Phi}^{(b)\mathrm{T}} dx$ are piecewise-polynomial integrals that can be written down exactly, so the PDE loss, boundary/initial-condition loss, and data misfit combine into $L(\mathbf{a}) = (\mathbf{d}-\mathbf{H}\mathbf{a})^{\mathrm{T}}(\mathbf{d}-\mathbf{H}\mathbf{a}) + \mathbf{a}^{\mathrm{T}}\mathbf{G}\mathbf{a}$, whose minimizer is $\mathbf{a}^{*} = (\mathbf{H}^{\mathrm{T}}\mathbf{H}+\mathbf{G})^{-1}\mathbf{H}^{\mathrm{T}}\mathbf{d}$. In two dimensions the PDE matrix $\mathbf{G}$ is assembled from Kronecker products of one-dimensional $\mathbf{R}_{ab}$ blocks, so arbitrary second-order linear PDEs with fixed coefficients fit the same machinery. The paper verifies this on forward and inverse examples, estimating a damping coefficient and a diffusion coefficient, reconstructing an initial temperature, and recovering harmonic solutions with error scaling about $M^{-3}$. It then applies the framework to GNSS velocities. The applied result is that the log marginal likelihood of the Bayesian model is higher for smoothness regularization ($-1842.75$) than for thin-elastic-sheet physical regularization ($-1947$ to $-1956$ depending on Poisson ratio), and adding the physical penalty to the mathematical one never improves the marginal likelihood. The paper also states its scope explicitly: PILM applies to linear PDEs with fixed coefficients on rectangular domains and is practical only in low dimensions, while nonlinear, variable-coefficient, complex-geometry, or high-dimensional problems remain the province of physics-informed neural networks.

Load-bearing premise

The headline comparison rests on the statistical noise model in Eq. (24), which treats every GNSS velocity component as independent with one common standard deviation $\sigma$; real GNSS velocities carry heterogeneous, correlated station-specific uncertainties, and if this assumption is wrong the roughly 105-nat marginal-likelihood advantage of mathematical smoothing could shrink or change sign.

Editorial extensions

If this is right

  • Linear forward problems for constant-coefficient PDEs can be solved without collocation points or iterative training; the only approximation is the basis-function representation itself.
  • Unknown PDE coefficients or boundary and initial values can be estimated by profiling: for each candidate coefficient, compute the optimal $\mathbf{a}$ in closed form and scan a one-dimensional loss curve, as done for the damping coefficient $c$ and the diffusion coefficient $k$.
  • The Bayesian marginal-likelihood criterion selects hyperparameters objectively; the paper finds mathematical smoothing strictly preferred over thin-elastic-sheet equilibrium for the Japanese GNSS velocity field.
  • Hybrid regularization that adds the physical penalty to the smoothness penalty does not increase marginal likelihood over pure mathematical regularization, so the physical constraint is statistically redundant or harmful for this dataset.
  • The same precomputed B-spline integral matrices are reusable for any linear PDE in one or two dimensions, giving an exact reference solution within a chosen model class against which neural-network solvers can be checked.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • Because the machinery only needs known integrals of basis functions, the same closed-form treatment should extend to other piecewise-polynomial or Fourier bases and to other constant-coefficient linear PDE families such as Helmholtz, Stokes, or poroelasticity; the paper itself demonstrates only ODE, diffusion, and thin-sheet elasticity.
  • The applied comparison assumes independent velocity components with a single common noise variance $\sigma$; re-estimating with station-specific or correlated GNSS covariance could change the roughly 105-nat marginal-likelihood gap, so that test is needed before reading the superiority of mathematical regularization as a general geophysical conclusion.
  • The author's own caveat suggests that point-force elastic interpolation might outperform the continuously distributed force-balance prior used here; PILM could test this directly by replacing $\mathbf{G}_{\mathrm{phys}}$ with a penalty derived from point forces.
  • The preference for smoothness over force balance may reflect the elastic model's inability to represent localized deformation rather than a failure of physics-based priors in general; testing on synthetic data generated from known elastic fields would separate the two.
Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

3 major / 6 minor

Summary. The paper introduces the Physics-Informed Linear Model (PILM), a method that represents PDE solutions as linear combinations of fixed basis functions (cubic B-splines) and computes the loss integrals analytically, yielding a closed-form least-squares solution. The formulation is demonstrated on a damped-oscillation ODE (forward and inverse problems) and a 1-D diffusion equation (inverse problem with unknown boundary/initial conditions), and is then applied to estimate crustal strain rates from GNSS velocities in central Japan. In that application, the author compares two regularizations—a 'mathematical' smoothness penalty and a 'physical' elastic-equilibrium penalty—within a Bayesian framework, using the log marginal likelihood for hyperparameter selection and model comparison. The paper's central methodological claim is that PILM provides an analytically solvable, reproducible alternative to PINNs for linear problems; its applied claim is that mathematical regularization is statistically superior to physical regularization for this GNSS dataset.

Significance. If the statistical comparison were valid, the paper would make a useful contribution at the intersection of physics-informed machine learning and geodetic inversion: it provides a transparent, exactly solvable baseline for PINNs, and it tackles a practical question about whether elasticity-based regularization is preferable to smoothness constraints in strain-rate estimation. The paper's strengths include analytical derivations with machine-checkable expressions, numerical verification against analytical solutions (Figures 2, 3, S3, S4), explicit scaling analysis in M, and publicly available code. However, the headline applied conclusion rests on a model comparison that is not well-defined as presented, and the noise model used in the comparison is unrealistic. These issues are load-bearing for the application, and therefore the applied claim cannot be accepted as established, although the methodological core appears sound.

major comments (3)
  1. [Section 3.2, Eq. (29)] The log marginal likelihood for the improper prior p(a) ∝ exp[-1/(2ρ_i^2) a^T G_i a] is defined only up to an additive constant that depends on the dimension and geometry of the null space of G_i. Since G_math and G_phys have different null spaces (for ν = 0.5, G_phys annihilates the mode u = x^2 − y^2, v = −2xy, so rank(G_phys) < rank(G_math)), the constant does not cancel when comparing LL values across rows of Table 2. The reported gap of about 105 nats is therefore not a well-defined Bayesian model comparison unless a proper prior is placed on the null-space components or another identifiability device is introduced. Please specify such a prior or replace the marginal-likelihood comparison with a criterion that is invariant to the choice of null-space parametrization (e.g., K-fold cross-validation of predictive performance).
  2. [Section 3.2, Eq. (24)] The likelihood assumes all GNSS velocity components are independent with a single common variance σ^2. GNSS velocity uncertainties are spatially heterogeneous and correlated, and the choice of noise model is known to affect strain-rate estimates and hyperparameter optimization. Because σ is optimized jointly with α_i^2 in the marginal-likelihood comparison, the ranking between mathematical and physical regularization in Table 2 could change if per-station uncertainties or a more realistic covariance were used. Please justify the common-σ assumption or repeat the comparison using reported per-station uncertainties.
  3. [Section 3.4, Eq. (34)] The hybrid regularization inherits the same improper-prior issue: the expression LL(α_i^2; G~_i) is used to compare models with different effective null spaces, and the conclusion that hybrid regularization never exceeds LL_math^∗ is therefore not well-defined as a Bayesian statement. In addition, the grid of α_phys^2 values shown in Figure 8b appears sparse (only four fixed values); a more thorough exploration or a formal marginalization over α_phys would be needed to support the claim of no improvement.
minor comments (6)
  1. [Abstract] The abstract contains typographical artifacts ('solving the se equations', 't he', 'understand ing'); these should be corrected in a final version.
  2. [Eq. (20)] The text 'Heere' should read 'Here'.
  3. [Table 1] The 'Solution' column describes PILM as 'Analytical', but for inverse problems with unknown coefficients (Sections 2.1, 2.2, 3.3), the coefficient is determined by a grid search over the loss L^*(c) or L^*(k), not by a fully closed-form expression. Please clarify that the closed-form solution applies to fixed model coefficients.
  4. [Eq. (29)] The notation log‖G_i‖ is nonstandard for a singular matrix; it should be defined as the sum of the logarithms of the nonzero eigenvalues (often written log |G_i|_+).
  5. [Figure S4] The residual magnitude for the quartic verification (3×10^−5) is much larger than for the quadratic verification (10^−13); a one-sentence explanation of this difference (e.g., due to boundary truncation or conditioning) would be helpful.
  6. [Section 3.3] The physical regularization for ν = −1 is described as 'decoupled' and for ν = 0.5, 0 as 'coupled', but the marginal likelihood values in Table 2 show ν = 0 is worse than ν = 0.5; a brief discussion of the non-monotonic dependence on ν would aid interpretation.

Circularity Check

0 steps flagged · score 0.0 of 10

No significant circularity: the PILM solutions and strain-rate comparisons follow from the stated loss functions and data, with self-citations used only for data, code, and auxiliary integrals.

full rationale

The derivation chain is self-contained. The closed-form solution (Eq. 10) is the analytic minimizer of the stated quadratic loss (Eq. 9); no target answer is used to construct it. The inverse estimates of the damping and diffusion coefficients are found by minimizing the same data-plus-physics loss over a grid and then evaluated against synthetic truth, which is standard empirical validation rather than circular reasoning. In the crustal strain-rate application, the two regularization matrices G_math and G_phys are inserted into the same marginal-likelihood formula (Eq. 29), and hyperparameters are optimized on the same GNSS data; this is model selection by the paper's own Bayesian criterion, not a fitted quantity renamed as a prediction. Self-citations (Okazaki et al. 2021 for GNSS velocities, Okazaki 2025 for code) and Nozue and Fukahata 2022 for R_ab integrals provide data and auxiliary values, but the derivation of the PILM solution does not depend on those citations as a load-bearing premise. The paper explicitly acknowledges limitations such as linear PDEs, rectangular domains, and the alternative localized-force approach, which further supports the independence of the results. The skeptical concern about missing null-space normalization in Eq. (29) is a statistical-modeling issue, not a circularity, and no quotable reduction of a claim to its own input exists in the manuscript.

Assumptions & free parameters 3 free parameters · 5 assumptions · 0 invented entities

The central method rests on the standard property of cubic B-splines and the algebraic form of the quadratic loss. The applied comparison rests on the Gaussian common-variance noise model and on the thin-elastic-sheet regularization. No new physical entities are introduced.

free parameters (3)
  • alpha_math^2 and alpha_phys^2 (regularization hyperparameters) = not reported
    Optimized by grid search over the marginal likelihood (Eqs. 28-29); their values are not tabulated, only the resulting log marginal likelihood and sigma.
  • Poisson's ratio nu = 0.5, 0, -1
    Chosen by hand in Section 3.3; controls the coupling terms G_xy and G_yx in the physical regularization.
  • Basis spacing and number of basis functions M = 20 km spacing, M=529 per component for GNSS; Delta t = 1 or 0.1 in synthetic examples
    Chosen by hand; determines resolution and edge artifacts, and affects the marginal likelihood comparison.
assumptions (5)
  • standard math Cubic B-splines are twice continuously differentiable and their definite integrals R_ab are analytically computable.
    Needed for Eqs. (7)-(8), (16), and (22); stated in Section 2.1 and Figure S1.
  • standard math The L2 loss in Eq. (2b) has global minimum zero if and only if u solves the PDE and BCs, and with a finite basis the closed-form solution is optimal within the chosen model class.
    Basis of Eqs. (9)-(10); the representational limitation is acknowledged in Sections 2.1 and 4.
  • domain assumption Gaussian likelihood with a single common variance sigma and Gaussian priors, with marginal likelihood as the model comparison criterion.
    Eqs. (23)-(29); load-bearing for the math-versus-phys comparison and never justified against heterogeneous GNSS errors.
  • domain assumption The crust behaves as a thin elastic sheet under plane stress, and minimizing the integral of the squared external force is a meaningful prior.
    Eqs. (20)-(21); the paper itself calls this a regularization rather than a strict governing equation (Sections 1 and 3.3).
  • domain assumption Basis functions truncated at domain boundaries, with edge effects treated as non-physical artifacts.
    Section 3.3 acknowledges artificial high strain rates at the margins, especially in the ocean and top-right corner.

how reviews work

0 comments
Cite this review

Pith. "Pith review of Physics-Informed Linear Model (PILM): Analytical Representations and Application to Crustal Strain Rate Estimation." pith.science (2026). https://pith.science/paper/DUXY6MO3

@misc{pith2026250712218,
  author       = {Pith},
  title        = {Pith review of: Physics-Informed Linear Model (PILM): Analytical Representations and Application to Crustal Strain Rate Estimation},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/DUXY6MO3}},
  note         = {Machine review of arXiv:2507.12218}
}
read the original abstract

Many physical systems are described by partial differential equations (PDEs), and solving these equations and estimating their coefficients or boundary conditions (BCs) from observational data play a crucial role in understanding the associated phenomena. Recently, a machine learning approach known as physics-informed neural network, which solves PDEs using neural networks by minimizing the sum of residuals from the PDEs, BCs, and data, has gained significant attention in the scientific community. In this study, we investigate a physics-informed linear model (PILM) that uses linear combinations of basis functions to represent solutions, thereby enabling an analytical representation of optimal solutions. The PILM was formulated and verified for illustrative forward and inverse problems including cases with uncertain BCs. Furthermore, the PILM was applied to estimate crustal strain rates using geodetic data. Specifically, physical regularization that enforces elastic equilibrium on the velocity fields was compared with mathematical regularization that imposes smoothness constraints. From a Bayesian perspective, mathematical regularization exhibited superior performance. The PILM provides an analytically solvable framework applicable to linear forward and inverse problems, underdetermined systems, and physical regularization.

Figures

Figures reproduced from arXiv: 2507.12218 by the authors.

Figure 1
Figure 1. Shape of the basis functions for 𝑇 = 10 and 𝑀 = 13 (Δ𝑡 = 1). Cubic B-spline functions are placed at equal intervals and are truncated at the domain boundaries. First, the damped oscillation was solved for 𝑇 = 10. The ODE coefficients were fixed at 𝑚 = 1 and 𝑘 = 1, and the damping coefficient was varied as 𝑐 = 0 (harmonic oscillation), 1 (underdamping), 2 (critical damping), and 3 (overdamping). Three ICs, (𝑢0, 𝑣0 ) … view at source ↗

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

6 extracted references · 5 canonical work pages

  1. [1]

    Akaike, H. (1980). Likelihood and the Bayes procedure. In Bayesian Statistics, University Press, Valencia, 143–166. https://doi.org/10.1007/978-1-4612-1694-0_24 Aki, K. & Richards, P. G. (2002). Quantitative Seismology, Second Edition, University Science Books. Bishop, C. M. (2006). Pattern Recognition and Machine Learning, Springer. 15 Brunton, S. L., & ...

    work page Pith review arXiv 1980

  2. [52]

    https://doi.org/10.1186/s40645-024-00654-7 Okazaki, T., Hirahara, K., Ito, T., Kano, M. and Ueda, N., (2025), Physics-informed deep learning for forward and inverse modeling of inplane crustal deformation, Journal of Geophysical Research: Machine Learning and Computation, 2, e2024JH000474. https://doi.org/10.1029/2024JH000474 Pagani, C., Bodin, T., Métois...

    work page Pith review arXiv 2025

  3. [53]

    https://doi.org/10.1186/s40623-021-01474-5 Okazaki, T., Ito, T., Hirahara, K., & Ueda, N. (2022). Physics -informed deep learning approach for modeling crustal deformation. Nature Communications, 13(1),

    work page doi:10.1186/s40623-021-01474-5 2022

  4. [58]

    https://doi.org/10.1186/s40623-022-01622-5 16 Nozue, Y ., & Fukahata, Y . (2022). Comparison between the boxcar and cubic B -spline functions in estimating displacement fields by basis function expansion. Zisin 2, 75, 1 –12. https://doi.org/10.4294/zisin.20218 (in Japanese with English abstract) Nozue, Y ., & Fukahata, Y . (2025). Influence of the 2011 To...

    work page doi:10.1186/s40623-022-01622-5 2022

  5. [80]

    https://doi.org/10.1186/s40623-025- 02190-0 Okazaki, T. (2025). Program codes for crustal strain rate estimation using basis function expansion . Zenodo. https://doi.org/10.5281/zenodo.16014772 Okazaki, T., Fukahata, Y ., & Nishimura, T. (2021). Consistent estimation of strain-rate fields from GNSS velocity data using basis function expansion with ABIC. E...

    work page doi:10.1186/s40623-025- 2025

  6. [7092]

    https://doi.org/10.1038/s41467- 022-34922-1 Okazaki, T., Hirahara, K., & Ueda, N. (2024). Fault geometry invariance and dislocation potential in antiplane crustal deformation: physics -informed simultaneous solutions. Progress in Earth and Planetary Science, 11,

    doi:10.1038/s41467- 2024

Pith tools

Reviewed August 6, 2026 · model on record in the stance chip above.