Hilbert–Galerkin Neural Operators

• HGNOs are neural operator ansatzes that represent inputs via finite Hilbert basis coefficients and reconstruct outputs in infinite-dimensional spaces.
• They leverage Galerkin projection to approximate not only target functions but also their Fréchet derivatives, supporting fully nonlinear PDE and HJB equation analysis.
• The architecture combines finite-dimensional latent computation with infinite-dimensional fidelity, enabling error-controlled approximations for complex PDEs and control problems.

Hilbert--Galerkin Neural Operators (HGNOs) are neural-operator ansatzes for separable Hilbert spaces in which an input is represented through finitely many basis coordinates, processed by a finite-dimensional neural network, and re-embedded into a Hilbert-space codomain. In the formulation developed for infinite-dimensional PDEs and optimal control, HGNOs are used to solve a single infinite-dimensional PDE directly on a separable Hilbert space $H$, rather than first discretizing the PDE into a finite-dimensional projected equation. The architecture is proved to approximate not only the target function but also its Fréchet derivative, second Fréchet derivative, and unbounded-operator-composed first-derivative terms required by fully nonlinear second-order PDEs and HJB equations on $H$ (Cohen et al., 19 Mar 2026). The mathematical significance of HGNOs is tied to a broader discretization problem: when an operator $F:X\to X$ on an infinite-dimensional Hilbert space can be approximated by finite-dimensional operators $F_V:V\to V$ in a convergent and representation-stable Galerkin sense, and when such a program is obstructed by infinite-dimensional topology (Furuya et al., 2024).

1. Architectural definition and Galerkin form

Let $H_1,H_2$ be separable Hilbert spaces with orthonormal bases $\{e_i\}_{i=1}^\infty\subset H_1$ and $\{g_i\}_{i=1}^\infty\subset H_2$. For $d,p\in\mathbb N$, the coordinate and embedding maps are

$$\mathcal E_d^{H_1}(x)=\big(\langle x,e_i\rangle\big)_{i=1}^d\in \mathbb R^d, \qquad \widehat{\mathcal E}_p^{H_2}(y_1,\dots,y_p)=\sum_{j=1}^p y_j g_j \in H_2.$$

The associated orthogonal projection is

$$P_d^H x=\sum_{i=1}^d \langle x,e_i\rangle e_i =\widehat{\mathcal E}_d^H(\mathcal E_d^H x).$$

An HGNO is then defined by

$H$0

where $H$1 is a trainable finite-dimensional neural network with parameters $H$2. This is explicitly a Galerkin-type representation: the input is truncated to the first $H$3 basis coordinates, and the output is expanded in the first $H$4 basis elements of the codomain. The finite-dimensional core is typically written as

$H$5

with activation $H$6, affine layers $H$7, and parameter set $H$8 (Cohen et al., 19 Mar 2026).

Several specializations are central. A scalar-valued value-function or PDE-solution HGNO has the form

$H$9

A control or operator HGNO is written

$F:X\to X$0

When $F:X\to X$1, the output may be postprocessed by a map $F:X\to X$2 to enforce $F:X\to X$3; for $F:X\to X$4, one example is $F:X\to X$5.

The architectural interpretation is therefore fixed throughout: analysis into finitely many Hilbert coordinates, nonlinear computation in Euclidean latent space, and synthesis back into the Hilbert codomain. A plausible implication is that the term “Galerkin” in HGNO refers primarily to basis truncation and reconstruction, rather than to a weak-form layer by layer finite-element assembly.

2. Infinite-dimensional PDEs, controls, and derivative representation

The principal PDE class treated with HGNOs is the fully nonlinear second-order equation on a separable Hilbert space $F:X\to X$6,

$F:X\to X$7

where $F:X\to X$8 is a closed, densely defined, possibly unbounded linear operator, $F:X\to X$9 is the Fréchet gradient, and $F_V:V\to V$0 is the bounded self-adjoint Hessian. The controlled setting is the stochastic evolution equation

$F_V:V\to V$1

with value function $F_V:V\to V$2 associated with the infinite-dimensional HJB equation

$F_V:V\to V$3

When $F_V:V\to V$4, this reduces to a deterministic first-order HJB equation; when the control is fixed or the admissible-control set is a singleton, one gets the Kolmogorov PDE (Cohen et al., 19 Mar 2026).

The derivative calculus of scalar HGNOs is explicit. For

$F_V:V\to V$5

one has

$F_V:V\to V$6

$F_V:V\to V$7

If $F_V:V\to V$8, then $F_V:V\to V$9 and

$H_1,H_2$0

These formulas are the basis of the “derivative-informed” character of HGNOs: once the Euclidean core is differentiable, the Hilbert-space derivatives needed by the PDE can be recovered by lifting finite-dimensional derivatives through the chosen basis.

A central obstruction is that cylindrical Hessian truncations generally do not converge in operator norm because

$H_1,H_2$1

Accordingly, the relevant topology for Hessians is weakened. On $H_1,H_2$2, the compact-open topology is generated by seminorms

$H_1,H_2$3

and on $H_1,H_2$4 by

$H_1,H_2$5

For $H_1,H_2$6, the weakened compact-open topology is generated by

$H_1,H_2$7

Within these topologies, HGNOs satisfy several universal approximation theorems. On compact sets they are dense in $H_1,H_2$8 with simultaneous approximation of $H_1,H_2$9, $\{e_i\}_{i=1}^\infty\subset H_1$0, and $\{e_i\}_{i=1}^\infty\subset H_1$1. In weighted Sobolev-type norms they are dense in $\{e_i\}_{i=1}^\infty\subset H_1$2. Under additional consistency assumptions they also approximate $\{e_i\}_{i=1}^\infty\subset H_1$3. These results feed directly into residual approximation theorems for the full PDE operator

$\{e_i\}_{i=1}^\infty\subset H_1$4

yielding HGNOs for which either

$\{e_i\}_{i=1}^\infty\subset H_1$5

or

$\{e_i\}_{i=1}^\infty\subset H_1$6

holds. The paper further proves a bounded-inverse estimate for certain dissipative Kolmogorov-type PDEs, so that a small residual implies a small $\{e_i\}_{i=1}^\infty\subset H_1$7 solution error (Cohen et al., 19 Mar 2026).

3. Discretization theory, no-go phenomena, and monotone resolution

The Hilbert--Galerkin aspect of HGNOs is not exhausted by basis truncation in the ansatz; it also raises a structural discretization question. In a general Hilbert-space setting one considers a directed family $\{e_i\}_{i=1}^\infty\subset H_1$8 of finite-dimensional subspaces with dense union, and a discretization family

$\{e_i\}_{i=1}^\infty\subset H_1$9

subject to the approximation requirement

$\{g_i\}_{i=1}^\infty\subset H_2$0

A stronger form requires

$\{g_i\}_{i=1}^\infty\subset H_2$1

The canonical linear discretization is

$\{g_i\}_{i=1}^\infty\subset H_2$2

with $\{g_i\}_{i=1}^\infty\subset H_2$3 the orthogonal projector onto $\{g_i\}_{i=1}^\infty\subset H_2$4. This is precisely the Hilbert-space Galerkin pipeline of restriction to a trial space, application of the continuum operator, and projection back to the trial space (Furuya et al., 2024).

For bijective neural operators framed as $\{g_i\}_{i=1}^\infty\subset H_2$5-diffeomorphisms, there is a fundamental obstruction: there exists no discretization scheme into finite-dimensional diffeomorphisms that is simultaneously convergent and continuous in the natural sense that $\{g_i\}_{i=1}^\infty\subset H_2$6 implies $\{g_i\}_{i=1}^\infty\subset H_2$7 for every fixed finite-dimensional $\{g_i\}_{i=1}^\infty\subset H_2$8. The source of the obstruction is topological. In finite dimensions, diffeomorphism groups split into orientation-preserving and orientation-reversing components; in infinite-dimensional Hilbert spaces, relevant invertible operator groups are path-connected. The proof uses a path joining identity to a reflection in infinite dimension together with finite-dimensional degree arguments. For HGNOs, the consequence is exact: invertibility or diffeomorphicity of a continuum operator does not by itself imply a well-posed discretization theory.

The positive resolution is strong monotonicity. A map $\{g_i\}_{i=1}^\infty\subset H_2$9 is strongly monotone if there exists $d,p\in\mathbb N$0 such that

$d,p\in\mathbb N$1

For a strongly monotone $d,p\in\mathbb N$2-diffeomorphism, each finite-dimensional restriction

$d,p\in\mathbb N$3

remains strongly monotone and is a $d,p\in\mathbb N$4-diffeomorphism. In finite dimensions the derivative matrices are then strictly positive definite, so determinants stay strictly positive and the orientation-switching obstruction disappears.

The neural-operator layer studied in this framework is

$d,p\in\mathbb N$5

where $d,p\in\mathbb N$6 are compact linear operators and $d,p\in\mathbb N$7. If

$d,p\in\mathbb N$8

then $d,p\in\mathbb N$9 is strongly monotone, hence a diffeomorphism. On this class the projection discretization $$\mathcal E_d^{H_1}(x)=\big(\langle x,e_i\rangle\big)_{i=1}^d\in \mathbb R^d, \qquad \widehat{\mathcal E}_p^{H_2}(y_1,\dots,y_p)=\sum_{j=1}^p y_j g_j \in H_2.$$0 satisfies the approximation property and is continuous. The same paper further shows that a bilipschitz neural-operator layer can be written on bounded balls as

$$\mathcal E_d^{H_1}(x)=\big(\langle x,e_i\rangle\big)_{i=1}^d\in \mathbb R^d, \qquad \widehat{\mathcal E}_p^{H_2}(y_1,\dots,y_p)=\sum_{j=1}^p y_j g_j \in H_2.$$1

where $$\mathcal E_d^{H_1}(x)=\big(\langle x,e_i\rangle\big)_{i=1}^d\in \mathbb R^d, \qquad \widehat{\mathcal E}_p^{H_2}(y_1,\dots,y_p)=\sum_{j=1}^p y_j g_j \in H_2.$$2 is either $$\mathcal E_d^{H_1}(x)=\big(\langle x,e_i\rangle\big)_{i=1}^d\in \mathbb R^d, \qquad \widehat{\mathcal E}_p^{H_2}(y_1,\dots,y_p)=\sum_{j=1}^p y_j g_j \in H_2.$$3 or the reflection

$$\mathcal E_d^{H_1}(x)=\big(\langle x,e_i\rangle\big)_{i=1}^d\in \mathbb R^d, \qquad \widehat{\mathcal E}_p^{H_2}(y_1,\dots,y_p)=\sum_{j=1}^p y_j g_j \in H_2.$$4

and each factor has the form

$$\mathcal E_d^{H_1}(x)=\big(\langle x,e_i\rangle\big)_{i=1}^d\in \mathbb R^d, \qquad \widehat{\mathcal E}_p^{H_2}(y_1,\dots,y_p)=\sum_{j=1}^p y_j g_j \in H_2.$$5

If $$\mathcal E_d^{H_1}(x)=\big(\langle x,e_i\rangle\big)_{i=1}^d\in \mathbb R^d, \qquad \widehat{\mathcal E}_p^{H_2}(y_1,\dots,y_p)=\sum_{j=1}^p y_j g_j \in H_2.$$6, then $$\mathcal E_d^{H_1}(x)=\big(\langle x,e_i\rangle\big)_{i=1}^d\in \mathbb R^d, \qquad \widehat{\mathcal E}_p^{H_2}(y_1,\dots,y_p)=\sum_{j=1}^p y_j g_j \in H_2.$$7.

This decomposition yields an architecture-level finite-rank realization. In a separable Hilbert space with orthonormal basis $$\mathcal E_d^{H_1}(x)=\big(\langle x,e_i\rangle\big)_{i=1}^d\in \mathbb R^d, \qquad \widehat{\mathcal E}_p^{H_2}(y_1,\dots,y_p)=\sum_{j=1}^p y_j g_j \in H_2.$$8,

$$\mathcal E_d^{H_1}(x)=\big(\langle x,e_i\rangle\big)_{i=1}^d\in \mathbb R^d, \qquad \widehat{\mathcal E}_p^{H_2}(y_1,\dots,y_p)=\sum_{j=1}^p y_j g_j \in H_2.$$9

and a residual network on $$P_d^H x=\sum_{i=1}^d \langle x,e_i\rangle e_i =\widehat{\mathcal E}_d^H(\mathcal E_d^H x).$$0 is

$$P_d^H x=\sum_{i=1}^d \langle x,e_i\rangle e_i =\widehat{\mathcal E}_d^H(\mathcal E_d^H x).$$1

Each block is a finite-rank residual perturbation of identity. The inverse of a block is obtained by contraction iteration: $$P_d^H x=\sum_{i=1}^d \langle x,e_i\rangle e_i =\widehat{\mathcal E}_d^H(\mathcal E_d^H x).$$2 This provides a rigorous discretization-safe class of invertible Hilbert-space neural operators and a local inversion mechanism compatible with Galerkin truncation.

4. Whole-space residual minimization and numerical realization

The numerical methodology built around HGNOs is explicitly formulated on the whole Hilbert space rather than on a projected PDE. The target residual norm is

$$P_d^H x=\sum_{i=1}^d \langle x,e_i\rangle e_i =\widehat{\mathcal E}_d^H(\mathcal E_d^H x).$$3

and the basic Monte Carlo loss for samples $$P_d^H x=\sum_{i=1}^d \langle x,e_i\rangle e_i =\widehat{\mathcal E}_d^H(\mathcal E_d^H x).$$4 is

$$P_d^H x=\sum_{i=1}^d \langle x,e_i\rangle e_i =\widehat{\mathcal E}_d^H(\mathcal E_d^H x).$$5

The defining point is that the ansatz is finite-dimensional, but the residual is evaluated at points $$P_d^H x=\sum_{i=1}^d \langle x,e_i\rangle e_i =\widehat{\mathcal E}_d^H(\mathcal E_d^H x).$$6, not merely at $$P_d^H x=\sum_{i=1}^d \langle x,e_i\rangle e_i =\widehat{\mathcal E}_d^H(\mathcal E_d^H x).$$7. In practice, $$P_d^H x=\sum_{i=1}^d \langle x,e_i\rangle e_i =\widehat{\mathcal E}_d^H(\mathcal E_d^H x).$$8 is sampled up to a larger truncation $$P_d^H x=\sum_{i=1}^d \langle x,e_i\rangle e_i =\widehat{\mathcal E}_d^H(\mathcal E_d^H x).$$9, so that terms such as $H$00 are evaluated using the full sampled state rather than the critic’s input truncation (Cohen et al., 19 Mar 2026).

Two training gradients are emphasized. For Deep Hilbert--Galerkin Methods (DHGM),

$H$01

For the QHPDE variant,

$H$02

For HJB equations the paper introduces Hilbert Actor--Critic, with critic $H$03 and actor $H$04. The controlled residual is

$H$05

and the actor objective minimizes the current-value Hamiltonian

$H$06

The numerical experiments reported in the paper concern heat and Burgers control problems on $H$07. For deterministic and stochastic heat equation control, the critic HGNO uses $H$08, one hidden layer, and 600 neurons; the actor uses $H$09, one hidden layer, and 600 neurons; training uses $H$10 iterations, $H$11 Monte Carlo samples per iteration, and state samples represented to $H$12 basis modes, with PyTorch + CUDA. Metrics are

$H$13

computed over $H$14 sampled points. For Kolmogorov problems, both DHGM and QHPDE perform well; for HJB actor-critic, QHPDE dramatically outperforms DHGM. In one trace-class covariance-noise example, QHPDE gives critic RMSE $H$15 and actor RMSE $H$16, whereas DHGM actor-critic errors are reported around order $H$17. The derivative errors under QHPDE satisfy

$H$18

while operator-norm Hessian errors remain much larger, consistent with the COCO-based theory.

For the controlled stochastic Burgers equation, the setup uses critic input dimension $H$19, state samples up to $H$20, an LSTM-like architecture from DGM rather than a simple feedforward network, about $H$21 more parameters than in the heat examples, QHPDE only, $H$22, $H$23, and $H$24. A Monte Carlo finite-difference solver on a 251-point grid with $H$25 simulations per evaluation point is used as baseline. The reported ensemble RMSE on 100-point test sets is around $H$26 in the stochastic case and $H$27 in the deterministic case (Cohen et al., 19 Mar 2026).

Problem	Setup	Reported outcome
Heat HJB/Kolmogorov	$H$28, $H$29, one hidden layer, 600 neurons, $H$30	QHPDE generally superior; actor-critic gains especially strong
Heat trace-class HJB	Same base setup	Critic RMSE $H$31, actor RMSE $H$32
Burgers control	$H$33, $H$34, LSTM-like DGM architecture	Ensemble RMSE around $H$35 stochastic, $H$36 deterministic

5. Relation to adjacent Hilbert-space neural operator theories

Two nearby developments are highly relevant to HGNOs even though neither is presented as an HGNO theory in name. The first studies a Residually-Guided Neural Operator (RNO) for the rules-to-equilibrium map of infinite-dimensional LQ mean field games on Hilbert spaces (Firoozi et al., 22 Oct 2025). The input “rule space” is a Hilbert space

$H$37

and the output space is

$H$38

With orthonormal bases $H$39 of $H$40 and $H$41 of $H$42, the projection and embedding maps are

$H$43

The RNO realization is

$H$44

The paper proves a regular universal approximation theorem in which the approximating operator preserves the same $H$45-Hölder class as the target, proves local Lipschitz continuity of the rules-to-equilibrium map

$H$46

and derives PAC-style generalization bounds under a Karhunen--Loève decomposition with exponentially decaying coefficients. The paper explicitly notes that it is “more ‘Hilbert-projection neural operator’ than ‘Galerkin neural operator’” in the narrow sense, since it does not use an explicit weak formulation or PDE discretization-consistency theorem.

The second adjacent framework is FrameNet, developed in a statistical learning theory for nonlinear operator regression between separable Hilbert spaces (Reinhardt et al., 2024). The architecture factors as

$H$47

where $H$48 encodes the input into frame coefficients, $H$49 rescales to a bounded cube, $H$50 is a finite-dimensional latent neural network, and $H$51 decodes coefficients back into the output Hilbert space. The paper proves general empirical-risk-minimization results for compact operator classes $H$52, including the mean-squared-error bound

$H$53

Under holomorphicity assumptions on the target operator and entropy bounds of the form $H$54, the sparse ReLU FrameNet class achieves

$H$55

A parametric elliptic PDE example on the torus yields algebraic statistical rates for a nonlinear coefficient-to-solution map between Sobolev Hilbert spaces. Structurally, FrameNet is very close to an HGNO: fixed Hilbert encoder, finite-dimensional latent map, and decoder back to the ambient space.

Taken together, these theories show that basis projection, finite-rank truncation, regularity control, and statistical learnability are not specific to one HGNO construction. They also clarify a distinction internal to the literature: some Hilbert-space neural operators are Galerkin-like because they rely on basis truncation, while others add a stronger claim of directly targeting an infinite-dimensional PDE residual or discretization-consistent finite-dimensional realizations.

6. Restrictions, common misconceptions, and open problems

A persistent misconception is that continuum invertibility automatically yields principled Galerkin discretization. The no-go theorem for arbitrary Hilbert-space diffeomorphisms shows that this is false: there is no discretization rule into finite-dimensional diffeomorphisms that is both convergent and continuous in the natural representation sense (Furuya et al., 2024). A second misconception is that using a finite-dimensional ansatz implies that the target problem has been reduced to a projected PDE. In the Deep Hilbert--Galerkin framework, this is explicitly rejected: the ansatz $H$56 is finite-dimensional, but the residual minimized is that of the full PDE on $H$57, evaluated at states represented to a larger truncation $H$58 (Cohen et al., 19 Mar 2026).

The limitations are equally clear. The PDE approximation theory is proved mainly for classical solutions $H$59 with $H$60. Approximation of unbounded-operator terms requires basis-dependent consistency assumptions. The sequential continuity assumption on the nonlinear operator $H$61 in the Hessian argument is nonstandard, and the relevant compact-open and strong-operator topologies are globally non-sequential and non-metrizable. On the discretization side, strong monotonicity is restrictive, bilipschitz factorization is local on bounded balls, the generalized neural-operator layer assumes compact linear maps $H$62, and inverse constructions are local and iterative rather than global closed forms. No general optimization-convergence theorem is proved, and practical performance remains dependent on basis choice, the sampling measure $H$63, and truncation levels.

Several open directions are identified in the source material. One is to determine how broad the class of physically relevant operators is that become strongly monotone after suitable reparameterization. Another is the choice of Galerkin spaces $H$64: the theory is largely basis-agnostic, but it does not prescribe optimal trial or test spaces. Sharper quantitative rates for bilipschitz decomposition and discretization remain open. The discretization of nonlinear activations $H$65 in norm topology is described as subtle unless one arranges invariant finite-dimensional spaces or changes layerwise spaces. A further question is whether one can build a full operator-learning theory for HGNOs in which training is consistent across multiple Galerkin levels under monotonicity constraints. These problems delimit the current state of the subject: HGNOs now have a precise Hilbert-space ansatz, strong approximation results for PDE residuals on $H$66, and a rigorous discretization-safe subclass, but not yet a fully general theory of infinite-dimensional operator learning without structural restrictions.