Stability and Bifurcation Analysis of Nonlinear PDEs via Random Projection-based PINNs: A Krylov-Arnoldi Approach

math.NA cs.LG cs.NA math.DS Gianluca Fabiani, Michail E. Kavousanakis, Constantinos Siettos, Ioannis G. Kevrekidis · Mar 23, 2026
Local to this browser
What it does
This paper solves stability and bifurcation analysis for nonlinear PDEs using Physics-Informed Random Projection Neural Networks (PI-RPNNs). The core innovation is a matrix-free shift-invert Krylov-Arnoldi method operating directly in...
Why it matters
The core innovation is a matrix-free shift-invert Krylov-Arnoldi method operating directly in weight space to circumvent the exponential singular value decay of the random collocation matrix $\Psi$. This enables reliable computation of...
Main concern
The paper presents a compelling methodological advance that bridges neural network-based PDE solvers with rigorous numerical bifurcation analysis. The combination of theoretical guarantees—almost sure regularity of the generalized...
Community signal
0
0 up · 0 down
Sign in to vote with arrows
AI Review AI reviewed
Plain-language introduction

This paper solves stability and bifurcation analysis for nonlinear PDEs using Physics-Informed Random Projection Neural Networks (PI-RPNNs). The core innovation is a matrix-free shift-invert Krylov-Arnoldi method operating directly in weight space to circumvent the exponential singular value decay of the random collocation matrix $\Psi$. This enables reliable computation of leading eigenpairs for detecting saddle-node, Hopf, and pitchfork bifurcations without requiring additional PDE solves beyond the initial training.

Critical review
Verdict
Bottom line

The paper presents a compelling methodological advance that bridges neural network-based PDE solvers with rigorous numerical bifurcation analysis. The combination of theoretical guarantees—almost sure regularity of the generalized eigenproblem and exponential singular value decay for analytic activations—with a robust algorithmic solution (truncated SVD with matrix-free Arnoldi) effectively addresses the critical ill-conditioning inherent to random projection methods. Validation on three canonical benchmarks (Liouville-Bratu-Gelfand, FitzHugh-Nagumo, and Allen-Cahn) demonstrates accurate detection of bifurcation points and recovery of eigenfunctions.

“The key advance presented here is the use of a matrix-free shift-invert Krylov–Arnoldi method, which allows the computation of the leading eigenvalues of the physical Jacobian without explicitly inverting the ill-conditioned random collocation matrix $\Psi$”
paper · Section 5 (Conclusion)
What holds up

The mathematical foundation is rigorous and convincing. Proposition 2 establishes that for analytic activation functions, the singular values of the collocation matrix decay exponentially as $\hat{\sigma}_{j} = O(R^{-j/2})$, rigorously justifying the observed numerical rank-deficiency. Proposition 1 guarantees that the generalized eigenvalue problem $J_w v = \lambda \tilde{B}_\Psi v$ is almost surely regular, ensuring compatibility with standard LAPACK/ARPACK solvers. The algorithmic strategy—avoiding explicit formation of $\Psi^\dagger$ via a factored, truncated SVD and applying shift-invert Arnoldi only through matrix-vector products—is well-designed and directly addresses the spurious mode contamination issue.

“the singular values $\hat{\sigma}_{1} \geq \hat{\sigma}_{2} \geq \dots \geq \hat{\sigma}_{M}$ of the collocation matrix $\Psi \in \mathbb{R}^{M \times N}$ ... satisfy the exponential decay estimate $\lim_{N \to \infty} \hat{\sigma}_{j} = O(R^{-j/2})$”
paper · Proposition 2 (Asymptotic singular value decay)
“the generalized eigenvalue problem ... defines a regular matrix pencil $(J_u, B)$ almost surely”
paper · Proposition 1 (Almost sure regularity)
Main concerns

While the benchmarks are standard, the empirical evaluation lacks comprehensive scaling studies or comparisons with established spectral/finite-element bifurcation software (e.g., AUTO, pde2path) regarding computational cost and memory scaling for high-resolution discretizations. The method relies on heuristic hyperparameter selection for random weight sampling (Eq. 41) whose bounds may not generalize to problems with steep gradients or widely separated scales without manual tuning. Furthermore, the framework is currently restricted to steady-state bifurcations; extensions to Floquet multipliers for periodic orbits or time-dependent stability analysis remain future work.

“The current implementation relies on heuristic bounds for the random parameters, which, while effective for the problems considered, may require adjustment for equations with very steep gradients or widely separated scales”
paper · Section 5 (Conclusion)
“extending it to directly compute Floquet multipliers for periodic orbits or to analyze stability of time-dependent PDEs would be a valuable direction”
paper · Section 5 (Conclusion)
Evidence and comparison

The evidence convincingly supports the core claims regarding numerical stability. Figure 2(a) demonstrates the exponential singular value decay of $\Psi$, while Figure 2(c) clearly contrasts the spurious eigenvalue cluster produced by the naive $J_u = J_w \Psi^\dagger$ approach against the clean spectrum obtained via the proposed shift-invert method. Comparisons with finite-difference references in Figures 1, 3, 4, 5, and 6 show good agreement in bifurcation diagrams and dominant eigenvalues. However, the paper does not provide quantitative runtime comparisons against concurrent neural network approaches (e.g., Eig-PIELM for linear PDEs) or traditional numerical continuation packages.

“While all three methods recover the correct dominant eigenvalue, the 'naive' approach also produces a cluster of spurious near-zero eigenvalues originating from the truncated singular modes of $\Psi$”
paper · Section 4.2 / Figure 2
“Concurrently, the Eig-PIELM approach [64] has proposed a similar generalized eigenvalue problem, but its scope is limited to linear PDEs and does not address the numerical conditioning challenges that arise in bifurcation regimes”
paper · Introduction
Reproducibility

The paper provides sufficient algorithmic detail in Section 3.3.1 to reproduce the method, including the matrix-free operator $\mathcal{T}_\sigma$ and the truncation tolerance rule $\tau \approx \sqrt{\epsilon_{\text{mach}}} \approx 10^{-8}$. The authors state that code and data will be made available upon publication. However, exact reproduction may be hindered by the lack of explicit tables listing problem-specific hyperparameters (e.g., random seed, exact collocation point counts $M$, neuron counts $N$, and distribution bounds $\alpha_U$) for each case study. The reliance on MATLAB's `eigs` is clearly stated, though the SVD truncation rank $r$ is determined adaptively based on tolerance rather than fixed a priori.

“we select a threshold $\tau$ (e.g., $\tau = 10^{-8}$), guided by a standard rule of thumb: in double-precision arithmetic (unit roundoff $\epsilon_{\text{mach}} \approx 2.2 \times 10^{-16}$), a tolerance $\tau \approx \sqrt{\epsilon_{\text{mach}}}$ suppresses the amplification of rounding errors”
paper · Remark 3.1
“The data and code supporting the findings of this study will be made available upon publication or upon reasonable request”
paper · Data and code availability
Abstract

We address a numerical framework for the stability and bifurcation analysis of nonlinear partial differential equations (PDEs) in which the solution is sought in the function space spanned by physics-informed random projection neural networks (PI-RPNNs), and discretized via a collocation approach. These are single-hidden-layer networks with randomly sampled and fixed a priori hidden-layer weights; only the linear output layer weights are optimized, reducing training to a single least-squares solve. This linear output structure enables the direct and explicit formulation of the eigenvalue problem governing the linear stability of stationary solutions. This takes a generalized eigenvalue form, which naturally separates the physical domain interior dynamics from the algebraic constraints imposed by boundary conditions, at no additional training cost and without requiring additional PDE solves. However, the random projection collocation matrix is inherently numerically rank-deficient, rendering naive eigenvalue computation unreliable and contaminating the true eigenvalue spectrum with spurious near-zero modes. To overcome this limitation, we introduce a matrix-free shift-invert Krylov-Arnoldi method that operates directly in weight space, avoiding explicit inversion of the numerically rank-deficient collocation matrix and enabling the reliable computation of several leading eigenpairs of the physical Jacobian - the discretized Frechet derivative of the PDE operator with respect to the solution field, whose eigenvalue spectrum determines linear stability. We further prove that the PI-RPNN-based generalized eigenvalue problem is almost surely regular, guaranteeing solvability with standard eigensolvers, and that the singular values of the random projection collocation matrix decay exponentially for analytic activation functions.

Challenge the Review

Pick a starting point or write your own. Challenges run in the background, so you can keep reading while the AI investigates.

No challenges yet. Disagree with the review? Ask the AI to revisit a specific claim.